diff --git a/Inputs/DetectorConfiguration/hira.detector b/Inputs/DetectorConfiguration/hira.detector
index 09b869fb7116a3347c0caeb8c6726a8e81ccc7b0..0420a979eb3256413e3d08bc8be9064e2948a289 100755
--- a/Inputs/DetectorConfiguration/hira.detector
+++ b/Inputs/DetectorConfiguration/hira.detector
@@ -9,7 +9,7 @@ GeneralTarget
Target
THICKNESS= 45
ANGLE= 0
- RADIUS= 15
+ RADIUS= 10
MATERIAL= CH2
X= 0
Y= 0
diff --git a/NPAnalysis/Hira/Analysis b/NPAnalysis/Hira/Analysis
index 276ccdd7d66539d5c152646f617309d7e4f19c40..a1eda041f6342c1023f4cda8ff998fcfc8212266 100755
Binary files a/NPAnalysis/Hira/Analysis and b/NPAnalysis/Hira/Analysis differ
diff --git a/NPAnalysis/Hira/Analysis.cxx b/NPAnalysis/Hira/Analysis.cxx
index 3416df30edbc3edb676fa8fde72518703dca0e28..20889f2d5abdd4687eb37e6955bc12039de3e65b 100644
--- a/NPAnalysis/Hira/Analysis.cxx
+++ b/NPAnalysis/Hira/Analysis.cxx
@@ -107,20 +107,20 @@ int main(int argc, char** argv)
TVector3 PositionOnHira = TVector3(X,Y,Z);
TVector3 ZUnit = TVector3(0,0,1);
- //TVector3 TelescopeNormal = -Must2->GetTelescopeNormal(countHira);
- //TVector3 TargetNormal = TVector3(-sin(TargetAngle*deg),0,cos(TargetAngle*deg));
+ //TVector3 TelescopeNormal = -Hira->GetTelescopeNormal(countHira);
+ TVector3 TargetNormal = TVector3(0,0,1);//TVector3(-sin(TargetAngle*deg),0,cos(TargetAngle*deg));
double X_target = InitialConditions->GetIncidentPositionX();
double Y_target = InitialConditions->GetIncidentPositionY();
double Z_target = InitialConditions->GetIncidentPositionZ();
- //TVector3 PositionOnTarget = TVector3(X_target,Y_target,Z_target);
- TVector3 PositionOnTarget = TVector3(0,0,0);
+ TVector3 PositionOnTarget = TVector3(X_target,Y_target,Z_target);
+ //TVector3 PositionOnTarget = TVector3(0,0,0);
TVector3 HitDirection = PositionOnHira - PositionOnTarget;
TVector3 HitDirectionUnit = HitDirection.Unit();
//double ThetaNormalHira = HitDirection.Angle(TelescopeNormal);
- //double ThetaNormalTarget = HitDirection.Angle(TargetNormal);
+ double ThetaNormalTarget = HitDirection.Angle(TargetNormal);
//double ThetaBeam = InitialConditions->GetICIncidentAngleTheta(countHira);
//double PhiBeam = InitialConditions->GetICIncidentAnglePhi(countHira);
@@ -145,19 +145,22 @@ int main(int argc, char** argv)
for(int countCsI =0; countCsI<Hira->CsI_E.size(); countCsI++){
E_CsI = Hira->CsI_E[countCsI];
if(E_CsI>EnergyThreshold)ELab += E_CsI;
- //if(E_CsI>EnergyThreshold) ELab = E_ThinSi + E_ThickSi + E_CsI;
}
}
- //else ELab = -1000;
+ else E_CsI = -1000;
+
+ ELab = EL_deuteron_CH2.EvaluateInitialEnergy(ELab,
+ (TargetThickness/2)*micrometer,
+ ThetaNormalTarget);
+
// ********************** Angle in the CM frame *****************************
TransfertReaction -> SetNuclei3(ELab, ThetaLab*deg);
ThetaCM = TransfertReaction->GetThetaCM()/deg;
ExcitationEnergy = TransfertReaction->GetExcitation4();
-
- RootOutput::getInstance()->GetTree()->Fill();
}
}
+ RootOutput::getInstance()->GetTree()->Fill();
}
cout << "A total of " << nentries << " event has been annalysed " << endl ;
diff --git a/NPAnalysis/Hira/Analysis.h b/NPAnalysis/Hira/Analysis.h
index b9ff62b5316a414db4b60458a2b9e4ebdc7f02e0..544b4424344bed22bde358b4e24a3fb2c0f4c84e 100644
--- a/NPAnalysis/Hira/Analysis.h
+++ b/NPAnalysis/Hira/Analysis.h
@@ -124,8 +124,9 @@ using namespace NPL ;
namespace ENERGYLOSS
{
+ EnergyLoss EL_deuteron_CH2 = EnergyLoss("deuteron_CH2.G4table","G4Table",100 );
// Declare your Energy loss here :
- /* EnergyLoss ProtonTarget = EnergyLoss ( "CD2.txt" ,
+ /* EnergyLoss DeuerontTarget = EnergyLoss ( "CD2.txt" ,
100 ,
1,
1 );
diff --git a/NPAnalysis/Hira/Analysis.o b/NPAnalysis/Hira/Analysis.o
index 3092ec6b048c0a7096d829ba14213ee9cd6d35a4..23bcb614e32cd6d9b95ed2fc3a516c031d442eb5 100644
Binary files a/NPAnalysis/Hira/Analysis.o and b/NPAnalysis/Hira/Analysis.o differ
diff --git a/NPAnalysis/Hira/Show.C b/NPAnalysis/Hira/Show.C
index 4b1106049dec09ed6b61d64df55025465fc33202..68361983baba7305f80f0c4bbc15e8d48b70943e 100644
--- a/NPAnalysis/Hira/Show.C
+++ b/NPAnalysis/Hira/Show.C
@@ -1,17 +1,36 @@
void InitChain();
+void LoadCut();
+
TChain *chain;
TGraph *kin;
+int nentries;
+
+TCutG *cut_kine;
+
+
void Show()
{
gStyle->SetOptStat(0);
InitChain();
+ LoadCut();
+
+ nentries = chain->GetEntries();
+
+ //**************** Condition ****************//
+ TString CondCut = "cut_kine";
+ TString CondTot = CondCut;//+"&&"+
+ //*******************************************//
TCanvas* mainC0 = new TCanvas("Hira XY", "Hira XY" , 800,600);
- TCanvas* mainC1 = new TCanvas("PID", "PID" , 800,600);
+ TCanvas* mainC1 = new TCanvas("PID", "PID" , 1200,600);
TCanvas* mainC2 = new TCanvas("Kinematics", "Kinematics" , 800,600);
TCanvas* mainC3 = new TCanvas("Phi-Theta", "Phi-Theta" , 800,600);
TCanvas* mainC4 = new TCanvas("BeamSpot", "BeamSpot" , 800,600);
+ TCanvas* mainC5 = new TCanvas("ThetaCM", "ThetaCM" , 800,600);
+ TCanvas* mainC6 = new TCanvas("Efficiency", "Efficiency" , 800,600);
+ TCanvas* mainC7 = new TCanvas("CrossSection","CrossSection",800,600);
+ mainC1->Divide(2,1);
mainC0->cd();
chain->Draw("Y:X>>h0(300,-300,300,200,-200,200)","","colz");
@@ -20,12 +39,21 @@ void Show()
h0->GetYaxis()->SetTitle("Y (mm)");
h0->SetTitle("");
- mainC1->cd();
+
+ mainC1->cd(1);
chain->Draw("E_ThickSi:E_CsI>>h1(1000,0,200,1000,0,25)","","colz");
TH2F* h1 = (TH2F*)gDirectory->FindObjectAny("h1");
h1->GetXaxis()->SetTitle("E_{CsI} (MeV)");
h1->GetYaxis()->SetTitle("E_{Si} (MeV)");
h1->SetTitle("");
+
+ mainC1->cd(2);
+ chain->Draw("E_ThinSi:E_ThickSi>>h12(1000,0,25,1000,0,5)","","colz");
+ TH2F* h12 = (TH2F*)gDirectory->FindObjectAny("h12");
+ h12->GetXaxis()->SetTitle("E_{Si} (MeV)");
+ h12->GetYaxis()->SetTitle("#Delta E (MeV)");
+ h12->SetTitle("");
+
mainC2->cd();
chain->Draw("ELab:ThetaLab>>h2(1000,0,50,1000,0,200)","","colz");
@@ -55,6 +83,32 @@ void Show()
h4->GetXaxis()->SetTitle("X_{beam} (mm)");
h4->GetYaxis()->SetTitle("Y_{beam} (mm)");
h4->SetTitle("");
+
+ mainC5->cd();
+ chain->Draw("ThetaCM>>h5(25,0,100)",CondTot,"E1");
+ TH1F* h5 = (TH1F*)gDirectory->FindObjectAny("h5");
+ h5->GetXaxis()->SetTitle("#theta_{CM} (deg)");
+ h5->SetTitle("");
+
+
+ //Efficiency calculation
+ mainC6->cd();
+ TH1F *Efficiency = new TH1F("Efficiency", "Efficiency", 25, 0, 100);
+ Efficiency->SetXTitle("#theta_{CM}");
+ Efficiency->SetTitle("Efficiency");
+ Efficiency->Add(h5,4*TMath::Pi()/nentries);
+ Efficiency->Draw();
+
+ // Cross-Section
+ mainC7->cd();
+ mainC7->SetLogy();
+ TH1F *h7 = new TH1F("h7", "h7", 25, 0, 100);
+ h7->SetXTitle("#theta_{CM}");
+ h7->SetYTitle("d#sigma/d#Omega (mb/sr)");
+ h7->SetTitle("");
+ h7->Divide(h5,Efficiency,1,nentries/(4*TMath::Pi()));//262000
+ h7->Draw();
+
return;
}
@@ -69,4 +123,26 @@ void InitChain()
chain->Add("/Users/pierremorfouace/Physics/NPTool/nptool/Outputs/Analysis/46Ar_pd.root");
return;
-}
\ No newline at end of file
+}
+
+
+
+///////////////////////////////////////////////////////////////////////////////
+void LoadCut()
+{
+ TString Cut_Name = "CUT/cut_kine.root";
+ TString Object_Name = "cut_kine";
+ TFile* f_cut_kine = new TFile(Cut_Name,"read");
+ cut_kine = (TCutG*) f_cut_kine->FindObjectAny(Object_Name);
+
+ return;
+}
+
+
+
+
+
+
+
+
+
diff --git a/NPSimulation/Hira/Hira.hh b/NPSimulation/Hira/Hira.hh
index f2a0e2b897a2787256475f982b8bf8cd28a0a30e..f049bbe9f6f00d0cac8a8b4a55803001e6e1196f 100644
--- a/NPSimulation/Hira/Hira.hh
+++ b/NPSimulation/Hira/Hira.hh
@@ -47,10 +47,10 @@ namespace HIRA
{
// Resolution
- const G4double ResoTime = 0.212765957 ;// = 500ps // Unit is ns/2.35
- const G4double ResoCsI = 0.106 ;// = 250 kev of resolution // Unit is MeV/2.35
- const G4double ResoThickSi = 0.051 ;// = 120keV of Resolution // Unit is MeV/2.35
- const G4double ResoThinSi = 0.034 ;// = 80keV of Resolution // Unit is MeV/2.35
+ const G4double ResoTime = 0.212765957; // = 500ps // Unit is ns/2.35
+ const G4double ResoCsI = 0.200/2.35; // = 250 kev of resolution // Unit is MeV/2.35
+ const G4double ResoThickSi = 0.065/2.35; // = 120keV of Resolution // Unit is MeV/2.35
+ const G4double ResoThinSi = 0.050/2.35; // = 80keV of Resolution // Unit is MeV/2.35
const G4double EnergyThreshold = 0.;//100*keV;
@@ -59,7 +59,7 @@ namespace HIRA
const G4double Length = 7.*cm ;
const G4int NumberOfStrip = 32;
- const G4double SiliconFace = 64*mm ;
+ const G4double SiliconFace = 63*mm ;
const G4double AluStripThickness = 0.4*micrometer ;
const G4double ThinSiThickness = 65*micrometer ;
const G4double ThickSiThickness = 1500*micrometer ;
@@ -68,10 +68,10 @@ namespace HIRA
const G4double MylarCsIThickness = 3*micrometer;
const G4double CsIThickness = 4.*cm;// + 2*MylarCsIThickness ;
- const G4double CsIXFront = 32.*mm;
+ const G4double CsIXFront = 33.*mm;
const G4double CsIXBack = 37.*mm;
- const G4double CsIYFront = 32.*mm;
- const G4double CsIYBack = 32.*mm;
+ const G4double CsIYFront = 33.*mm;
+ const G4double CsIYBack = 33.*mm;
const G4double DistInterCsI = 0.01*mm;
const G4double ClusterFaceFront = 7*cm;