From c57deb54863cfa9575b615f0ba081b180170999b Mon Sep 17 00:00:00 2001
From: adrien-matta <a.matta@surrey.ac.uk>
Date: Tue, 21 Oct 2014 10:19:07 +0100
Subject: [PATCH] * Updating Example1 for nicer presentation and adding
comments to NPA * project
---
NPAnalysis/Example1/Analysis.cxx | 89 ++++++++++++++++++------------
NPAnalysis/Example1/Analysis.h | 24 ++++----
NPAnalysis/Example1/ShowResult.C | 17 ++++--
NPAnalysis/Example1/cuts/EDE.root | Bin 4256 -> 4256 bytes
4 files changed, 76 insertions(+), 54 deletions(-)
diff --git a/NPAnalysis/Example1/Analysis.cxx b/NPAnalysis/Example1/Analysis.cxx
index 1e3b646b0..575a3f970 100755
--- a/NPAnalysis/Example1/Analysis.cxx
+++ b/NPAnalysis/Example1/Analysis.cxx
@@ -1,7 +1,6 @@
#include "Analysis.h"
-int main(int argc, char** argv)
-{
+int main(int argc, char** argv){
// command line parsing
NPOptionManager* myOptionManager = NPOptionManager::getInstance(argc,argv);
@@ -57,43 +56,53 @@ int main(int argc, char** argv)
double Si_Y_M2 = 0;
double ZTarget = 0;
double TargetThickness = 18*micrometer;
- double ThetaCM = 0 ;
// Get number of events to treat
cout << endl << "///////// Starting Analysis ///////// "<< endl;
int nentries = Chain->GetEntries();
- // int nentries = 100000;
cout << " Number of Event to be treated : " << nentries << endl;
clock_t begin = clock();
clock_t end = begin;
cout.precision(5);
- // main loop on entries
+ //////////////////////////////////////////////////////////////////////////////
+ // main loop on entries //
for (int i = 0 ; i < nentries; i++) {
if (i%10000 == 0 && i!=0) {
end = clock();
long double TimeElapsed = (long double) (end-begin) / CLOCKS_PER_SEC;
double percent = (double)i/nentries;
double TimeToWait = (TimeElapsed/percent) - TimeElapsed;
- cout << " "<< flush;
- cout << "\r Progression:" << percent*100 << " % \t | \t Remaining time : ~" << TimeToWait <<"s | Analysis Rate : "<< (double) i/TimeElapsed << flush;
+ cout << " "<< flush;
+ cout << "\r Progression:"
+ << percent*100 << " % \t | \t Remaining time : ~"
+ << TimeToWait <<"s | Analysis Rate : "
+ << (double) i/TimeElapsed << flush;
}
else if (i == nentries-1) cout << "\r Progression:" << " 100% " <<endl;
- // get data
+ // Get the raw Data
Chain -> GetEntry(i);
+ // Clear previous Event
myDetector->ClearEventPhysics();
+ // Build the current event
myDetector->BuildPhysicalEvent();
+ // Reinitiate calculated variable
ReInitValue();
- // Get the Init information on beam positio nand energy
+ // Get the Init information on beam position and energy
// and apply by hand the experimental resolution
- double XTarget = Rand.Gaus(Init->GetIncidentPositionX(),1);
+ // This is because the beam diagnosis are not simulated
+
+ // PPAC position resolution on target is assumed to be 1mm
+ double XTarget = Rand.Gaus(Init->GetIncidentPositionX(),1);
double YTarget = Rand.Gaus(Init->GetIncidentPositionY(),1);
- double BeamEnergy = Rand.Gaus(Init->GetIncidentInitialKineticEnergy(),4.5);
TVector3 BeamDirection = Init->GetBeamDirection();
+
+ // Beam energy is measured using F3 and F2 plastic TOF
+ double BeamEnergy = Rand.Gaus(Init->GetIncidentInitialKineticEnergy(),4.5);
BeamEnergy = Li11CD2.Slow(BeamEnergy,TargetThickness/2.,0);
He10Reaction->SetBeamEnergy(BeamEnergy);
-
+
//////////////////////////// LOOP on MUST2 + SSSD Hit //////////////////
for(unsigned int countSSSD = 0 ; countSSSD < SSSD->Energy.size() ; countSSSD++){
for(unsigned int countMust2 = 0 ; countMust2 < M2->Si_E.size() ; countMust2++){
@@ -145,30 +154,34 @@ int main(int argc, char** argv)
E_SSSD = SSSD->Energy[countSSSD];
// if CsI
- if(CsI_E_M2>0 ){
- // The energy in CsI is calculate form dE/dx Table bevause calibration for 3He is not good
- Energy = He3Si.EvaluateEnergyFromDeltaE(Si_E_M2,300*micrometer,ThetaM2Surface, 0.01*MeV, 450.*MeV,0.001*MeV ,1000);
- E_M2=CsI_E_M2;
- }
-
- else
- Energy = Si_E_M2;
-
- E_M2 += Si_E_M2;
-
- // Evaluate energy using the thickness
- ELab = He3Al.EvaluateInitialEnergy( Energy ,2*0.4*micrometer , ThetaM2Surface);
- ELab = He3Si.EvaluateInitialEnergy( ELab ,20*micrometer , ThetaM2Surface);
- ELab = He3Al.EvaluateInitialEnergy( ELab ,0.4*micrometer , ThetaM2Surface);
- // Target Correction
- ELab = He3CD2.EvaluateInitialEnergy( ELab ,TargetThickness/2., ThetaNormalTarget);
+ if(CsI_E_M2>0 ){
+ // The energy in CsI is calculate form dE/dx Table because
+ // 20um resolution is poor
+ Energy =
+ He3Si.EvaluateEnergyFromDeltaE(Si_E_M2,300*micrometer,
+ ThetaM2Surface, 0.01*MeV,
+ 450.*MeV,0.001*MeV ,1000);
+ E_M2=CsI_E_M2;
+ }
+
+ else
+ Energy = Si_E_M2;
+
+ E_M2 += Si_E_M2;
+
+ // Evaluate energy using the thickness
+ ELab = He3Al.EvaluateInitialEnergy( Energy ,2*0.4*micrometer , ThetaM2Surface);
+ ELab = He3Si.EvaluateInitialEnergy( ELab ,20*micrometer , ThetaM2Surface);
+ ELab = He3Al.EvaluateInitialEnergy( ELab ,0.4*micrometer , ThetaM2Surface);
+ // Target Correction
+ ELab = He3CD2.EvaluateInitialEnergy( ELab ,TargetThickness/2., ThetaNormalTarget);
/************************************************/
- /***********M2->Si_E[countMust2];*************************************/
+ /***********888888888888888888888****************/
// Part 3 : Excitation Energy Calculation
Ex = He10Reaction -> ReconstructRelativistic( ELab , ThetaLab );
ThetaLab=ThetaLab/deg;
- /************************************************/
+ /************************************************/
/************************************************/
// Part 4 : Theta CM Calculation
@@ -187,27 +200,31 @@ int main(int argc, char** argv)
RootOutput::getInstance()->Destroy();
RootInput::getInstance()->Destroy();
NPOptionManager::getInstance()->Destroy();
- /////////////////////////////////////////////////////////////////////
+ /////////////////////////////////////////////////////////////////////////////
return 0 ;
}
-////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
void InitOutputBranch() {
RootOutput::getInstance()->GetTree()->Branch("Ex",&Ex,"Ex/D");
RootOutput::getInstance()->GetTree()->Branch("ELab",&ELab,"ELab/D");
RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab,"ThetaLab/D");
+ RootOutput::getInstance()->GetTree()->Branch("ThetaCM",&ThetaCM,"ThetaCM/D");
}
-////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
void InitInputBranch(){
RootInput:: getInstance()->GetChain()->SetBranchAddress("InitialConditions",&Init );
RootInput:: getInstance()->GetChain()->SetBranchStatus("InitialConditions",true );
RootInput:: getInstance()->GetChain()->SetBranchStatus("fIC_*",true );
}
-////////////////////////////////////////////////////////////////////////////////////////////////////
-////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
void ReInitValue(){
Ex = -1000 ;
ELab = -1000;
ThetaLab = -1000;
+ ThetaCM = -1000;
}
+
+
+////////////////////////////////////////////////////////////////////////////////
diff --git a/NPAnalysis/Example1/Analysis.h b/NPAnalysis/Example1/Analysis.h
index 76c2d8d4e..73dd9916a 100755
--- a/NPAnalysis/Example1/Analysis.h
+++ b/NPAnalysis/Example1/Analysis.h
@@ -40,29 +40,25 @@ void ReInitValue() ;
/////////////////////////////////////////////////////////////////////////////////////////////////
// ----------------------------------- DOUBLE, INT, BOOL AND MORE -------------------------------
namespace VARIABLE{
- double Ex;
+ double Ex;
double ELab;
double ThetaLab;
+ double ThetaCM;
TInitialConditions* Init = new TInitialConditions();
}
-
+
using namespace VARIABLE ;
// ----------------------------------------------------------------------------------------------
/////////////////////////////////////////////////////////////////////////////////////////////////
// -----------------------------------ENERGY LOSS----------------------------------------------
-namespace ENERGYLOSS
- {
-
- // Declare your Energy loss here :
- EnergyLoss He3CD2 = EnergyLoss ("He3_CD2.G4table","G4Table",1 );
-
- EnergyLoss He3Al = EnergyLoss ("He3_Aluminium.G4table","G4Table",1);
-
- EnergyLoss He3Si = EnergyLoss ("He3_Si.G4table","G4Table",1);
-
- EnergyLoss Li11CD2 = EnergyLoss ("Li11[0.0]_CD2.G4table","G4Table",100);
+namespace ENERGYLOSS{
+ // Energy loss table: the G4Table are generated by the simulation
+ EnergyLoss He3CD2 = EnergyLoss("He3_CD2.G4table","G4Table",100 );
+ EnergyLoss He3Al = EnergyLoss("He3_Aluminium.G4table","G4Table",10);
+ EnergyLoss He3Si = EnergyLoss("He3_Si.G4table","G4Table",10);
+ EnergyLoss Li11CD2 = EnergyLoss("Li11[0.0]_CD2.G4table","G4Table",100);
}
-
+
using namespace ENERGYLOSS ;
// ----------------------------------------------------------------------------------------------
/////////////////////////////////////////////////////////////////////////////////////////////////
diff --git a/NPAnalysis/Example1/ShowResult.C b/NPAnalysis/Example1/ShowResult.C
index 6a137a24d..d425115db 100644
--- a/NPAnalysis/Example1/ShowResult.C
+++ b/NPAnalysis/Example1/ShowResult.C
@@ -41,9 +41,10 @@ ETOF->Draw("same");
// Kinematical Line //
c1->cd(3);
-chain->Draw("ELab:ThetaLab>>hKine(1000,0,90,400,0,40)","MUST2.CsI_E<0 && MUST2.TelescopeNumber<5 && EDE && ETOF","colz");
+chain->Draw("ELab:ThetaLab>>hKine(500,0,45,400,0,40)","MUST2.CsI_E<0 && MUST2.TelescopeNumber<5 && EDE && ETOF","colz");
NPL::Reaction r("11Li(d,3He)10He@553");
+r.SetExcitationHeavy(1.4);
TGraph* Kine = r.GetKinematicLine3();
Kine->SetLineWidth(2);
Kine->SetLineColor(kOrange-3);
@@ -51,9 +52,9 @@ Kine->Draw("c");
// Excitation Energy //
c1->cd(4);
-int bin=100;
-double Emin = -10;
-double Emax = 10;
+int bin=50;
+double Emin = -5;
+double Emax = 5;
chain->Draw(Form("Ex>>hEx(%d,%f,%f)",bin,Emin,Emax),"MUST2.CsI_E<0 && MUST2.TelescopeNumber<5 && EDE && ETOF");
TH1F* hEx = (TH1F*) gDirectory->FindObjectAny("hEx");
@@ -62,4 +63,12 @@ hEx->GetXaxis()->SetTitle("E_{10He}");
hEx->SetFillStyle(1001);
hEx->SetLineColor(kAzure+7);
hEx->SetFillColor(kAzure+7);
+
+hEx->Fit("gaus");
+
+TF1* f = hEx->GetFunction("gaus");
+f->SetLineWidth(2);
+f->SetLineColor(kOrange-3);
+f->SetNpx(1000);
+
}
diff --git a/NPAnalysis/Example1/cuts/EDE.root b/NPAnalysis/Example1/cuts/EDE.root
index 1d18123c0647fc8d0dce0661e0b1cabf411b40a0..fea217bf9828f59a729befc535082e5d38399a78 100644
GIT binary patch
delta 500
zcmZ3WxIl4&zQC8)b9y3~1fNV7pYrW-^S<{J9n~Zm7^)ap7#JA5{Z=heVbEt_1hP|r
zm|=3Gfat`@Ldrm4uoMGBEKt-ACe6Sgz{S7_F>&GvMP{I~iEk#?zxe$A&}_l-#8;BC
zGn)L5|DMftqyBSn)P_CveSdeSPI_9^EB|JG8~eFw7Bz?dIiJ1!c9K2cjsF{77w#+j
z!0EG>nYsB-y#mL_lRA2d44+seBNFWxEIB!3BxX!+npz>TtI@((LPAzTV!z=A>jgii
zcQW4HuQ0#yO4Hxd4M+E9I_#|XVe?ba<Ow^?z-7+HcCG73*O{aV$srOE5)u;s8&$Wi
z=|2;fwz_8S>1VGO*8NV+%-m*|evT=1OM_LAjL4+}2M+$>;h4kPmB_&UW?JZ^EXl=x
zH_JAzkA7Nq>whFr*R}-CHeUxunQ5Vk+j&y>>U-2~r1SFdD41PIe{h6l(}911%pFHq
zXT(34<S5>5DfNuAq;kS(!^;etR+*%n)2RN?a-gx(FyZ)t488jo9tbG!eQ>5C{M~~G
z3F5Zb_gq=@#fw30=7OyZ47^^TDAB3~h5P1nOzym((48zHU;qha88KiaSb}4X3lVDq
T>5~-%q&R``C$OO0bP59i%EQVi
delta 500
zcmZ3WxIl4&zQC$YX0sxk1)oe8pYrW-^S<{J9n~Zm7^)ap7#JA5{Z?$d$e_=_2xO-K
zF~j6W0nv$*g_ME9U?~QMSfHpKOqzi~fQx|<V&cRTip)S|6W>g(f41e0+w0OegReb&
zJ_naS{_D&7X8%u3-H3bk$N$|mI$5;$82>l_!_4X4GIojo7e0Ggc2b`G&HjD6o)#Nw
zFfF^+*y#ArUZCYiir91mhRP<MH3sqxb68mTc+Pkq^s?c(>o7->r-zTH=f1=Z*?>Q%
zk2JiCpK#tG>)^lC1CQ>TEQqOBV^*Ia!nQh<p;elhIombKHO=sZp%#x04^Pkk1Cyez
zxu@xG4&8Sy_4C(&{eL%YF^QgUoZe^@<shrc$Ca7zAfb+p<s8#dg9hfZ(^{TedoKRF
z*>~Xj+9K<>|8;=6qBpP{R$IWpcY4)^Xf~sEdpF@Xo7vdd1f*XX7bG!7Cj1X%ba}*d
zM!(?Xg6_jIyq{Ti*_?PfBa<OARLUrQifx6{g9Aq-HY6vQOwZ3K=$IH+khVv=yr5u1
z_q^=5tiUQ&hAM^m@(c{TUZ5z^ss)An=5tK$yr9sXEFoY331t~EU?f<AV~q<DYXa$$
R6$GR>f$=A>pxkr{0{}3?#q9t9
--
GitLab