Skip to content
Snippets Groups Projects
Commit 20bd19a3 authored by Pierre Morfouace's avatar Pierre Morfouace
Browse files

New more realistic resolution for Hira

parent 05f7dd15
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Inputs/DetectorConfiguration/hira.detector
+ 1
− 1
View file @ 20bd19a3
......@@ -9,7 +9,7 @@ GeneralTarget
Target
THICKNESS= 45
ANGLE= 0
RADIUS= 15
RADIUS= 10
MATERIAL= CH2
X= 0
Y= 0
......
NPAnalysis/Hira/Analysis
+ 0
− 0
View file @ 20bd19a3
No preview for this file type
NPAnalysis/Hira/Analysis.cxx
+ 12
− 9
View file @ 20bd19a3
......@@ -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 ;
......
NPAnalysis/Hira/Analysis.h
+ 2
− 1
View file @ 20bd19a3
......@@ -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 );
......
NPAnalysis/Hira/Analysis.o
+ 0
− 0
View file @ 20bd19a3
No preview for this file type
NPAnalysis/Hira/Show.C
+ 79
− 3
View file @ 20bd19a3
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;
}
NPSimulation/Hira/Hira.hh
+ 8
− 8
View file @ 20bd19a3
......@@ -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;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment