diff --git a/Inputs/DetectorConfiguration/hira.detector b/Inputs/DetectorConfiguration/hira.detector
index eef6c620062453f45b48ff57b8bcc2888a64478b..0420a979eb3256413e3d08bc8be9064e2948a289 100755
--- a/Inputs/DetectorConfiguration/hira.detector
+++ b/Inputs/DetectorConfiguration/hira.detector
@@ -7,9 +7,9 @@ GeneralTarget
% Radius in mm
% Temperature in K, Pressure in bar
Target
- THICKNESS= 75
+ THICKNESS= 45
ANGLE= 0
- RADIUS= 5
+ RADIUS= 10
MATERIAL= CH2
X= 0
Y= 0
diff --git a/NPAnalysis/Hira/Analysis b/NPAnalysis/Hira/Analysis
index dd799aa3aea2ede7dcd37fee9934c8883833ed50..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 bbe9ee3598425a55a4d151701c9e4af0ed6d06fd..20889f2d5abdd4687eb37e6955bc12039de3e65b 100644
--- a/NPAnalysis/Hira/Analysis.cxx
+++ b/NPAnalysis/Hira/Analysis.cxx
@@ -33,7 +33,7 @@ int main(int argc, char** argv)
// Instantiate some Reaction
NPL::Reaction* TransfertReaction = new Reaction ;
- TransfertReaction -> ReadConfigurationFile("34Ar_pd.reaction") ;
+ TransfertReaction -> ReadConfigurationFile("46Ar_pd.reaction") ;
//Get Detector pointer :
THiraPhysics* Hira = (THiraPhysics*) myDetector -> GetDetector("HIRAArray") ;
@@ -58,7 +58,7 @@ int main(int argc, char** argv)
RootOutput::getInstance()->GetTree()->Branch( "E_ThinSi" , &E_ThinSi , "E_ThinSi/D" ) ;
RootOutput::getInstance()->GetTree()->Branch( "E_ThickSi" , &E_ThickSi , "E_ThickSi/D" ) ;
RootOutput::getInstance()->GetTree()->Branch( "E_CsI" , &E_CsI , "E_CsI/D" ) ;
- RootOutput::getInstance()->GetTree()->Branch( "Etot" , &Etot , "Etot/D" ) ;
+ RootOutput::getInstance()->GetTree()->Branch( "ELab" , &ELab , "ELab/D" ) ;
RootOutput::getInstance()->GetTree()->Branch("ExcitationEnergy", &ExcitationEnergy,"ExcitationEnergy/D") ;
RootOutput::getInstance()->GetTree()->Branch( "X" , &X , "X/D" ) ;
RootOutput::getInstance()->GetTree()->Branch( "Y" , &Y , "Y/D" ) ;
@@ -72,7 +72,7 @@ int main(int argc, char** argv)
clock_t begin = clock();
clock_t end = begin;
- double EnergyThreshold = 1;
+ double EnergyThreshold = 0.;
// main loop on entries
for (int i = 0; i < nentries; i++) {
@@ -107,19 +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 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);
@@ -139,24 +140,27 @@ int main(int argc, char** argv)
E_ThickSi = Hira->ThickSi_E[countHira];
E_ThinSi = Hira->ThinSi_E[countHira];
- //Etot = E_ThinSi + E_ThickSi;
+ ELab = E_ThinSi + E_ThickSi;
if(Hira->CsI_E.size() == 1){
for(int countCsI =0; countCsI<Hira->CsI_E.size(); countCsI++){
E_CsI = Hira->CsI_E[countCsI];
- //Etot += E_CsI;
- if(E_CsI>EnergyThreshold) Etot = E_ThinSi + E_ThickSi + E_CsI;
+ if(E_CsI>EnergyThreshold)ELab += E_CsI;
}
}
- else Etot = -1000;
+ else E_CsI = -1000;
+
+ ELab = EL_deuteron_CH2.EvaluateInitialEnergy(ELab,
+ (TargetThickness/2)*micrometer,
+ ThetaNormalTarget);
+
// ********************** Angle in the CM frame *****************************
- TransfertReaction -> SetNuclei3(Etot, ThetaLab*deg);
+ 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 ;
@@ -176,7 +180,7 @@ void ReInitOuputValue()
E_ThinSi = -1000;
E_ThickSi = -1000;
E_CsI = -1000;
- Etot = -1000;
+ ELab = -1000;
ThetaLab = -1000;
PhiLab = -1000;
ThetaCM = -1000;
diff --git a/NPAnalysis/Hira/Analysis.h b/NPAnalysis/Hira/Analysis.h
index 118d5dcda277b9d784644d2e02f5bb0555ca7511..544b4424344bed22bde358b4e24a3fb2c0f4c84e 100644
--- a/NPAnalysis/Hira/Analysis.h
+++ b/NPAnalysis/Hira/Analysis.h
@@ -43,7 +43,7 @@ Double_t TargetThickness;
Double_t E_ThinSi = 0;
Double_t E_ThickSi = 0;
Double_t E_CsI = 0;
-Double_t Etot = 0;
+Double_t ELab = 0;
Double_t ThetaLab = 0;
Double_t PhiLab = 0;
Double_t ThetaLab_simu = 0;
@@ -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 8c71ea1604d3b9d97e2c261196785011dd8d2bb4..23bcb614e32cd6d9b95ed2fc3a516c031d442eb5 100644
Binary files a/NPAnalysis/Hira/Analysis.o and b/NPAnalysis/Hira/Analysis.o differ
diff --git a/NPAnalysis/Hira/RunToTreat.txt b/NPAnalysis/Hira/RunToTreat.txt
index d571a10601cbf55fb73c24260b02e28eed845713..706313c9f4cc50aa5fece9fcba99378f4e75187d 100644
--- a/NPAnalysis/Hira/RunToTreat.txt
+++ b/NPAnalysis/Hira/RunToTreat.txt
@@ -1,8 +1,8 @@
TTreeName
SimulatedTree
RootFileName
- /Users/pierremorfouace/Physics/NPTool/nptool/Outputs/Simulation/simu.root
- %/Users/pierremorfouace/Physics/NPTool/nptool/Outputs/Simulation/11Be_d3He.root
+ %/Users/pierremorfouace/Physics/NPTool/nptool/Outputs/Simulation/simu.root
+ /Users/pierremorfouace/Physics/NPTool/nptool/Outputs/Simulation/46Ar_pd.root
diff --git a/NPAnalysis/Hira/Show.C b/NPAnalysis/Hira/Show.C
index 4bb9b5c994ffcc81420337f8c300249f16cbd941..68361983baba7305f80f0c4bbc15e8d48b70943e 100644
--- a/NPAnalysis/Hira/Show.C
+++ b/NPAnalysis/Hira/Show.C
@@ -1,16 +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");
@@ -19,27 +39,36 @@ 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("Etot:ThetaLab>>h2(1000,0,50,1000,0,200)","","colz");
+ chain->Draw("ELab:ThetaLab>>h2(1000,0,50,1000,0,200)","","colz");
h2->SetMinimum(1);
TH2F* h2 = (TH2F*)gDirectory->FindObjectAny("h2");
h2->GetXaxis()->SetTitle("#theta_{lab} (deg)");
h2->GetYaxis()->SetTitle("E (MeV)");
h2->SetTitle("");
- NPL::Reaction *r = new NPL::Reaction("34Ar(p,d)33Ar@2380");
- //NPL::Reaction *r = new NPL::Reaction("11Be(d,3He)10Li@770");
+ //NPL::Reaction *r = new NPL::Reaction("34Ar(p,d)33Ar@2380");
+ NPL::Reaction *r = new NPL::Reaction("46Ar(p,d)45Ar@3220");
kin = r->GetKinematicLine3();
kin->SetLineColor(2);
- kin->Draw("plsame");
+ //kin->Draw("plsame");
mainC3->cd();
chain->Draw("PhiLab:ThetaLab>>h3(100,0,50,720,-180,180)","","colz");
@@ -47,6 +76,39 @@ void Show()
h3->GetXaxis()->SetTitle("#theta_{lab} (deg)");
h3->GetYaxis()->SetTitle("#phi_{lab} (MeV)");
h3->SetTitle("");
+
+ mainC4->cd();
+ chain->Draw("fIC_Incident_Position_Y:fIC_Incident_Position_X>>h4(1000,-10,10,100,-10,10)","","colz");
+ TH2F* h4 = (TH2F*)gDirectory->FindObjectAny("h4");
+ 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;
}
@@ -58,7 +120,29 @@ void InitChain()
if(chain!=NULL){delete chain;}
chain = new TChain("AnalysedTree","");
- chain->Add("/Users/pierremorfouace/Physics/NPTool/nptool/Outputs/Analysis/test.root");
+ chain->Add("/Users/pierremorfouace/Physics/NPTool/nptool/Outputs/Analysis/46Ar_pd.root");
+
+ return;
+}
+
+
+
+///////////////////////////////////////////////////////////////////////////////
+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;
-}
\ No newline at end of file
+}
+
+
+
+
+
+
+
+
+
diff --git a/NPAnalysis/newHira/Analysis.cxx b/NPAnalysis/newHira/Analysis.cxx
index 887d1f0d037f10390ffebb60edf2d8a465e07d25..a4384b8fa057454681917a5ba0bd8e7de8074ca8 100644
--- a/NPAnalysis/newHira/Analysis.cxx
+++ b/NPAnalysis/newHira/Analysis.cxx
@@ -36,11 +36,11 @@ void Analysis::Init(){
InitOutputBranch();
InitInputBranch();
- Hira= (THiraPhysics*) m_DetectorManager -> GetDetector("Sharc");
+ Hira= (THiraPhysics*) m_DetectorManager -> GetDetector("HIRAArray");
LightCD2 = EnergyLoss("proton_CD2.G4table","G4Table",100 );
LightAl = EnergyLoss("proton_Al.G4table","G4Table",100);
LightSi = EnergyLoss("proton_Si.G4table","G4Table",100);
- BeamCD2 = EnergyLoss("Rb88_CD2.G4table","G4Table",100);
+ //BeamCD2 = EnergyLoss("Rb88_CD2.G4table","G4Table",100);
myReaction = new NPL::Reaction();
myReaction->ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
TargetThickness = m_DetectorManager->GetTargetThickness()*micrometer;
diff --git a/NPLib/Physics/NPBeam.cxx b/NPLib/Physics/NPBeam.cxx
index ce69aa19cb7e4f2788bf444dac1778cda971ea69..93c90f27f3a8e0867d01fbc2ca6373b157c57dad 100644
--- a/NPLib/Physics/NPBeam.cxx
+++ b/NPLib/Physics/NPBeam.cxx
@@ -333,14 +333,14 @@ void Beam::GenerateRandomEvent(double& E, double& X, double& Y, double& Z, doubl
if(fSigmaX!=-1){
// Shoot within the target unless target size is null (no limit)
- while(sqrt(X*X+Y*Y>fEffectiveTargetSize) || fEffectiveTargetSize == 0){
+ while(sqrt(X*X+Y*Y)>fEffectiveTargetSize || fEffectiveTargetSize == 0){
NPL::RandomGaussian2D(fMeanX, fMeanThetaX, fSigmaX, fSigmaThetaX, X, ThetaX);
NPL::RandomGaussian2D(fMeanY, fMeanPhiY, fSigmaY, fSigmaPhiY, Y, PhiY);
}
}
else{
- while(sqrt(X*X+Y*Y>fEffectiveTargetSize) || fEffectiveTargetSize == 0){
+ while(sqrt(X*X+Y*Y)>fEffectiveTargetSize || fEffectiveTargetSize == 0){
fXThetaXHist->GetRandom2(X,ThetaX);
fYPhiYHist->GetRandom2(Y,PhiY);
}
diff --git a/NPSimulation/Hira/Hira.hh b/NPSimulation/Hira/Hira.hh
index e7b0d4ebacca217f690f0cb5d748ea035830a1b3..f049bbe9f6f00d0cac8a8b4a55803001e6e1196f 100644
--- a/NPSimulation/Hira/Hira.hh
+++ b/NPSimulation/Hira/Hira.hh
@@ -1,5 +1,5 @@
-#ifndef Sharc_h
-#define Sharc_h 1
+#ifndef Hira_h
+#define Hira_h 1
/*****************************************************************************
* Copyright (C) 2009-2013 this file is part of the NPTool Project *
* *
@@ -14,7 +14,7 @@
* Last update : *
*---------------------------------------------------------------------------*
* Decription: *
- * This class describe the Sharc Silicon detector *
+ * This class describe the Hira Telescops *
* *
*---------------------------------------------------------------------------*
* Comment: *
@@ -47,10 +47,10 @@ namespace HIRA
{
// Resolution
- const G4double ResoTime = 0.212765957 ;// = 500ps // Unit is ns/2.35
- const G4double ResoCsI = 0.08 ;// = 188 kev of resolution // Unit is MeV/2.35
- const G4double ResoThickSi = 0.022 ;// = 52keV of Resolution // Unit is MeV/2.35
- const G4double ResoThinSi = 0.064 ;// = 150keV 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,11 +68,11 @@ 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 DistInterCsI = 0.2*mm;
+ const G4double CsIYFront = 33.*mm;
+ const G4double CsIYBack = 33.*mm;
+ const G4double DistInterCsI = 0.01*mm;
const G4double ClusterFaceFront = 7*cm;
const G4double ClusterFaceBack = 9*cm;