diff --git a/Projects/e793s/Analysis.cxx b/Projects/e793s/Analysis.cxx
index 302b3276c051d2626c9e7f431ba81adf0b635f97..90db9a8e24d6bc4e141ac306325154ecb6a1b7a4 100755
--- a/Projects/e793s/Analysis.cxx
+++ b/Projects/e793s/Analysis.cxx
@@ -42,11 +42,24 @@ Analysis::~Analysis(){
////////////////////////////////////////////////////////////////////////////////
void Analysis::Init() {
- ///////////////////////////////////////////////////////////////////////////////
+ ///////////////////////////////////////////////////////////////////////////////
+
+// if(NPOptionManager::getInstance()->HasDefinition("simulation")){
+ cout << " == == == == SIMULATION == == == ==" << endl;
+ isSim=true;
+// } else {
+// cout << " == == == == EXPERIMENT == == == ==" << endl;
+// isSim=false;
+// }
agata_zShift=51*mm;
//BrhoRef=0.65;
-
+
+ if(isSim){
+ Initial = new TInitialConditions();
+ ReactionConditions = new TReactionConditions();
+ }
+
// initialize input and output branches
InitOutputBranch();
InitInputBranch();
@@ -73,6 +86,11 @@ void Analysis::Init() {
// target thickness
TargetThickness = m_DetectorManager->GetTargetThickness();
string TargetMaterial = m_DetectorManager->GetTargetMaterial();
+
+ cout << " ---------------- Target Values ---------------- " << endl;
+ cout << " \tZ = " << m_DetectorManager->GetTargetZ() << " mm\tT = " << TargetThickness << " mm" << endl;
+ cout << " ----------------------------------------------- " << endl;
+
// energy losses
string light = NPL::ChangeNameToG4Standard(reaction.GetParticle3()->GetName());
string beam = NPL::ChangeNameToG4Standard(reaction.GetParticle1()->GetName());
@@ -110,26 +128,47 @@ void Analysis::Init() {
//nbTrack=0;
//nbHits=0;
//count=0;
- //AHeavy=reaction.GetNucleus4()->GetA();
AHeavy=reaction.GetParticle4()->GetA();
- //ALight=reaction.GetNucleus3()->GetA();
ALight=reaction.GetParticle3()->GetA();
- //MHeavy=reaction.GetNucleus4()->Mass();
MHeavy=reaction.GetParticle4()->Mass();
- //MLight=reaction.GetNucleus3()->Mass();
MLight=reaction.GetParticle3()->Mass();
+ bool writetoscreen=true;
for(int i=0;i<GATCONF_SIZE;i++){ // loop over the bits
GATCONF_Counter[i] = 0 ;
}
+ ThetaCM_detected->Sumw2();
+ ThetaLab_detected->Sumw2();
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::TreatEvent(){
-
+
// Reinitiate calculated variable
ReInitValue();
+
+ if(isSim){
+ //ThetaCM_emmitted->Fill(Initial->GetThetaCM(0));
+ ThetaLab_emmitted->Fill(Initial->GetThetaLab_WorldFrame(0));
+
+ TVector3 labDir = Initial->GetParticleDirection(0);
+ TVector3 comDir = labDir;
+
+ comDir.SetZ( cos(0.5*Initial->GetThetaLab_WorldFrame(0)*deg) );
+ comDir.SetX( labDir.X() / (2*comDir.Z()) );
+ comDir.SetY( labDir.Y() / (2*comDir.Z()) );
+
+ /*
+ cout << Initial->GetThetaLab_WorldFrame(0) << " "
+ << labDir.Theta()/deg << " "
+ << comDir.Theta()/deg << " "
+ << endl;
+ */
+
+ ThetaCM_emmitted->Fill(comDir.Theta()/deg);
+ }
+
GATCONF_MASTER=ML->GetCalibratedValue("GATCONF_MASTER");
for(int i=0;i<GATCONF_SIZE;i++){ // loop over the bits
if(GATCONF_MASTER & (unsigned int)pow(2,i)){ // test if ith bit is on
@@ -138,6 +177,12 @@ void Analysis::TreatEvent(){
}
}
+// if(isSim){
+// OriginalELab = ReactionConditions->GetKineticEnergy(0);
+// OriginalThetaLab = ReactionConditions->GetTheta(0);
+// BeamEnergy = ReactionConditions->GetBeamEnergy();
+// }
+
TVector3 BeamDirection(0.,0.,1.);
BeamImpact = TVector3(XBeam,YBeam,m_DetectorManager->GetTargetZ());
@@ -146,6 +191,7 @@ void Analysis::TreatEvent(){
// FinalBeamEnergy=BeamCD2.Slow(OriginalBeamEnergy,0.5*TargetThickness,BeamDirection.Angle(TVector(0,0,1)));
//reaction.SetBeamEnergy(FinalBeamEnergy);
+
////////////////////////////////////////////////////////////////////////////
//////////////////////////////// LOOP on MUST2 ////////////////////////////
////////////////////////////////////////////////////////////////////////////
@@ -157,6 +203,31 @@ void Analysis::TreatEvent(){
// MUST2
int TelescopeNumber = M2->TelescopeNumber[countMust2];
+ if(isSim){
+ if(TelescopeNumber==5){
+ //ThetaCM_detected->Fill(Initial->GetThetaCM(0));
+ ThetaLab_detected->Fill(Initial->GetThetaLab_WorldFrame(0));
+
+ TVector3 labDir = Initial->GetParticleDirection(0);
+ TVector3 comDir = labDir;
+
+ comDir.SetZ( cos(0.5*Initial->GetThetaLab_WorldFrame(0)*deg) );
+ comDir.SetX( labDir.X() / (2*comDir.Z()) );
+ comDir.SetY( labDir.Y() / (2*comDir.Z()) );
+
+ /*
+ // Check the angles agree
+ cout << Initial->GetThetaLab_WorldFrame(0) << " "
+ << labDir.Theta()/deg << " "
+ << comDir.Theta()/deg << " "
+ << endl;
+ */
+
+ ThetaCM_detected->Fill(comDir.Theta()/deg);
+
+ }
+ }
+
/************************************************/
// Part 1 : Impact Angle
ThetaM2Surface = 0;
@@ -165,6 +236,36 @@ void Analysis::TreatEvent(){
philab_tmp = 0;
TVector3 HitDirection = M2 -> GetPositionOfInteraction(countMust2) - BeamImpact ;
thetalab_tmp = HitDirection.Angle(BeamDirection);
+
+ /******************************************************************************/
+ /*** THIS SECTION IS FOR USE WHEN CONVERTING ISOTROPIC LAB TO ISOTROPIC CoM ***/
+// if(writetoscreen){
+// cout << "====================================" << endl;
+// cout << " CONVERTING ISO LAB TO ISO CoM! " << endl;
+// cout << "====================================" << endl;
+// writetoscreen=false;
+// }
+ /* Generate a normalized vector/unit vector for use in theta calculation*/
+// TVector3 NormHitDir = HitDirection.Unit();
+
+ /* thetaCM = 0.5 thetaLab = 0.5 acos(zLab) */
+// thetalab_tmp = 0.5 * acos(NormHitDir.Z());
+ /* Now, thetalab_tmp is in CoM frame */
+
+ /* zCM = cos(thetaCM) */
+// HitDirection.SetZ(cos(thetalab_tmp));
+ /* Now, (xLab, yLab, zCM) */
+
+ /* xCM = xLab/(2 zCM) */
+// HitDirection.SetX(HitDirection.X()/(2*HitDirection.Z()));
+ /* Now, (xCM, yLab, zCM) */
+
+ /* yCM = yLab/(2 zCM) */
+// HitDirection.SetY(HitDirection.Y()/(2*HitDirection.Z()));
+ /* Now, (xCM, yCM, zCM) */
+ /******************************************************************************/
+ /******************************************************************************/
+
philab_tmp = HitDirection.Phi();
X.push_back(M2->GetPositionOfInteraction(countMust2).X());
@@ -172,7 +273,7 @@ void Analysis::TreatEvent(){
Z.push_back(M2->GetPositionOfInteraction(countMust2).Z());
ThetaM2Surface = HitDirection.Angle(- M2->GetTelescopeNormal(countMust2) );
- ThetaNormalTarget = HitDirection.Angle( TVector3(XBeam,YBeam,1) ) ;
+ ThetaNormalTarget = HitDirection.Angle( TVector3(0.0,0.0,1) ) ;
/************************************************/
// Part 2 : Impact Energy
@@ -185,19 +286,30 @@ void Analysis::TreatEvent(){
if(CsI_E_M2>0 ){
// The energy in CsI is calculate form dE/dx Table because
Energy = CsI_E_M2;
- Energy = LightAl.EvaluateInitialEnergy(Energy, 0.4*micrometer, ThetaM2Surface);
+ if(!isSim){
+ Energy = LightAl.EvaluateInitialEnergy(Energy, 0.4*micrometer, ThetaM2Surface);
+ }
Energy+=Si_E_M2;
}
-
else
Energy = Si_E_M2;
RawEnergy.push_back(Energy);
- // Evaluate energy using the thickness
- elab_tmp = LightAl.EvaluateInitialEnergy(Energy, 0.4*micrometer, ThetaM2Surface);
- // Target Correction
- elab_tmp = LightTarget.EvaluateInitialEnergy(elab_tmp, 0.5*TargetThickness, ThetaNormalTarget);
+ /*
+ cout << M2->TelescopeNumber[0] << " "
+ << M2->GetTelescopeNormal(countMust2).X() << " "
+ << M2->GetTelescopeNormal(countMust2).Y() << " "
+ << M2->GetTelescopeNormal(countMust2).Z() << endl;
+ */
+
+ if(!isSim){
+ // Evaluate energy using the thickness
+ elab_tmp = LightAl.EvaluateInitialEnergy(Energy, 0.4*micrometer, ThetaM2Surface);
+ // Target Correction
+ elab_tmp = LightTarget.EvaluateInitialEnergy(elab_tmp, 0.5*TargetThickness, ThetaNormalTarget);
+ } else {elab_tmp = Energy;}
+
ELab.push_back(elab_tmp);
/************************************************/
@@ -205,7 +317,6 @@ void Analysis::TreatEvent(){
/************************************************/
// Part 3 : Excitation Energy Calculation
Ex.push_back(reaction.ReconstructRelativistic(elab_tmp,thetalab_tmp));
-
Ecm.push_back(Energy*(AHeavy+ALight)/(4*AHeavy*cos(thetalab_tmp)*cos(thetalab_tmp)));
/************************************************/
@@ -217,9 +328,9 @@ void Analysis::TreatEvent(){
ThetaLab.push_back(thetalab_tmp/deg);
PhiLab.push_back(philab_tmp/deg);
-
}
+/*
////////////////////////////////////////////////////////////////////////////
//////////////////////////////// LOOP on MUGAST ////////////////////////////
////////////////////////////////////////////////////////////////////////////
@@ -239,27 +350,37 @@ void Analysis::TreatEvent(){
Z.push_back( MG -> GetPositionOfInteraction(countMugast).Z());
ThetaMGSurface = HitDirection.Angle( MG -> GetTelescopeNormal(countMugast) );
- ThetaNormalTarget = HitDirection.Angle( TVector3(XBeam,YBeam,-1) ) ;
+ ThetaNormalTarget = HitDirection.Angle( TVector3(0.0,0.0,1.0) ) ;
+
+//cout << MG->TelescopeNumber[0] << "\t"
+// << "\t" << MG->GetTelescopeNormal(countMugast).X()
+// << "\t" << MG->GetTelescopeNormal(countMugast).Y()
+// << "\t" << MG->GetTelescopeNormal(countMugast).Z() << endl;
// Part 2 : Impact Energy
Energy = elab_tmp = 0;
Energy = MG->GetEnergyDeposit(countMugast);
RawEnergy.push_back(Energy);
- elab_tmp = LightAl.EvaluateInitialEnergy(
+ if(!isSim){
+ elab_tmp = LightAl.EvaluateInitialEnergy(
Energy, //particle energy after Al
0.4*micrometer, //thickness of Al
ThetaMGSurface); //angle of impingement
- elab_tmp = LightTarget.EvaluateInitialEnergy(
+ elab_tmp = LightTarget.EvaluateInitialEnergy(
elab_tmp, //particle energy after leaving target
TargetThickness*0.5, //distance passed through target
ThetaNormalTarget); //angle of exit from target
+ } else {elab_tmp = Energy;}
+
ELab.push_back(elab_tmp);
// Part 3 : Excitation Energy Calculation
- Ex.push_back(reaction.ReconstructRelativistic(elab_tmp,thetalab_tmp));
- Ecm.push_back(elab_tmp*(AHeavy+ALight)/(4*AHeavy*cos(thetalab_tmp)*cos(thetalab_tmp)));
-
+ if(!isSim){
+ Ex.push_back(reaction.ReconstructRelativistic(elab_tmp,thetalab_tmp));
+ Ecm.push_back(elab_tmp*(AHeavy+ALight)/(4*AHeavy*cos(thetalab_tmp)*cos(thetalab_tmp)));
+ }
+
// Part 4 : Theta CM Calculation
ThetaLab.push_back(thetalab_tmp/deg);
PhiLab.push_back(philab_tmp/deg);
@@ -275,6 +396,8 @@ void Analysis::TreatEvent(){
}//end loop Mugast
+*/
+
////////////////////////////////////////////////////////////////////////////
///////////////////////////////// LOOP on AGATA ////////////////////////////
////////////////////////////////////////////////////////////////////////////
@@ -309,6 +432,8 @@ void Analysis::TreatEvent(){
EDC.push_back(GammaLV.Energy());
}
*/
+
+/********
// Agata add back is not always multiplicity 1 ?? NO, not necessarily!
for(int i=0; i<nbAdd; i++){
//if(nbAdd==1){
@@ -339,6 +464,7 @@ void Analysis::TreatEvent(){
GammaLV.Boost(beta);
// Get EDC
AddBack_EDC.push_back(GammaLV.Energy());
+ AddBack_EDC2.push_back(GammaLV.Energy());
}
@@ -352,44 +478,103 @@ void Analysis::TreatEvent(){
MWT[i] = (MW_T[i]-offset[j])/slope[j];
}
+********/
}
- ////////////////////////////////////////////////////////////////////////////////
- void Analysis::End(){
- cout << endl << "\e[1;32m GATCONF Statistics " << endl ;
- for(int i=0;i<GATCONF_SIZE;i++){ // loop over the bits
- cout << GATCONF_Label[i] << "\t\t" << GATCONF_Counter[i] << endl ;
- }
- cout << endl ;
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::End(){
+ cout << endl << "\e[1;32m GATCONF Statistics " << endl ;
+ for(int i=0;i<GATCONF_SIZE;i++){ // loop over the bits
+ cout << GATCONF_Label[i] << "\t\t" << GATCONF_Counter[i] << endl ;
}
+ cout << endl ;
+
+ if(isSim){
+
+ TObjArray HistList(0);
+
+ auto Efficiency_CM = new TH1F(*ThetaCM_detected);
+ //auto Efficiency_CM = new TH1F(*ThetaCM_emmitted);
+ Efficiency_CM->SetName("Efficiency_CM");
+ Efficiency_CM->SetTitle("Efficiency_CM");
+ Efficiency_CM->Sumw2();
+ auto Efficiency_Lab = new TH1F(*ThetaLab_detected);
+ //auto Efficiency_Lab = new TH1F(*ThetaLab_emmitted);
+ Efficiency_Lab->SetName("Efficiency_Lab");
+ Efficiency_Lab->SetTitle("Efficiency_Lab");
+ Efficiency_Lab->Sumw2();
+
+ auto SolidAngle_CM = new TH1F(*ThetaCM_detected);
+ //auto SolidAngle_CM = new TH1F(*ThetaCM_emmitted);
+ SolidAngle_CM->SetName("SolidAngle_CM");
+ SolidAngle_CM->SetTitle("SolidAngle_CM");
+ auto SolidAngle_Lab = new TH1F(*ThetaLab_detected);
+ //auto SolidAngle_Lab = new TH1F(*ThetaLab_emmitted);
+ SolidAngle_Lab->SetName("SolidAngle_Lab");
+ SolidAngle_Lab->SetTitle("SolidAngle_Lab");
+
+ Efficiency_CM->Divide(ThetaCM_emmitted);
+ //Efficiency_CM->Divide(ThetaCM_detected);
+ Efficiency_Lab->Divide(ThetaLab_emmitted);
+ //Efficiency_Lab->Divide(ThetaLab_detected);
+
+ double dt = 180./Efficiency_Lab->GetNbinsX();
+ cout << "Angular infinitesimal = " << dt << "deg " << endl;
+ auto Cline = new TF1("Cline",Form("1./(2*%f*sin(x*%f/180.)*%f*%f/180.)",M_PI,M_PI,dt,M_PI),0,180);
+
+ SolidAngle_CM->Divide(ThetaCM_emmitted);
+ //SolidAngle_CM->Divide(ThetaCM_detected);
+ SolidAngle_CM->Divide(Cline,1);
+ SolidAngle_Lab->Divide(ThetaLab_emmitted);
+ //SolidAngle_Lab->Divide(ThetaLab_detected);
+ SolidAngle_Lab->Divide(Cline,1);
+
+ HistList.Add(ThetaCM_emmitted);
+ HistList.Add(ThetaLab_emmitted);
+ HistList.Add(ThetaCM_detected);
+ HistList.Add(ThetaLab_detected);
+ HistList.Add(Efficiency_CM);
+ HistList.Add(Efficiency_Lab);
+ HistList.Add(SolidAngle_CM);
+ HistList.Add(SolidAngle_Lab);
+ HistList.Add(Cline);
+
+ //HistoFile->Write();
+ auto HistoFile = new TFile("SolidAngle_HistoFile.root","RECREATE");
+ HistList.Write();
+ HistoFile->Close();
+ }
+}
////////////////////////////////////////////////////////////////////////////////
- void Analysis::InitOutputBranch(){
- //RootOutput::getInstance()->GetTree()->Branch("Ex",&Ex,"Ex/D");
- RootOutput::getInstance()->GetTree()->Branch("Ex",&Ex);
- //RootOutput::getInstance()->GetTree()->Branch("EDC",&EDC,"EDC/D");
- RootOutput::getInstance()->GetTree()->Branch("EDC",&EDC);
- RootOutput::getInstance()->GetTree()->Branch("AddBack_EDC",&AddBack_EDC);
- RootOutput::getInstance()->GetTree()->Branch("EAgata",&EAgata,"EAgata/D");
- RootOutput::getInstance()->GetTree()->Branch("ELab",&ELab);
- RootOutput::getInstance()->GetTree()->Branch("Ecm",&Ecm);
- RootOutput::getInstance()->GetTree()->Branch("RawEnergy",&RawEnergy); // CPx ADDITION
- RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab);
- RootOutput::getInstance()->GetTree()->Branch("PhiLab",&PhiLab);
- RootOutput::getInstance()->GetTree()->Branch("ThetaCM",&ThetaCM);
- RootOutput::getInstance()->GetTree()->Branch("Run",&Run,"Run/I");
- RootOutput::getInstance()->GetTree()->Branch("X",&X);
- RootOutput::getInstance()->GetTree()->Branch("Y",&Y);
- RootOutput::getInstance()->GetTree()->Branch("Z",&Z);
- RootOutput::getInstance()->GetTree()->Branch("dE",&dE,"dE/D");
- RootOutput::getInstance()->GetTree()->Branch("MG_T",MG_T);
- RootOutput::getInstance()->GetTree()->Branch("MG_E",MG_E);
- RootOutput::getInstance()->GetTree()->Branch("MG_X",MG_X);
- RootOutput::getInstance()->GetTree()->Branch("MG_Y",MG_Y);
- RootOutput::getInstance()->GetTree()->Branch("MG_D",MG_D);
-
- // Vamos
- RootOutput::getInstance()->GetTree()->Branch("LTS",<S,"LTS/l");
+void Analysis::InitOutputBranch(){
+ //RootOutput::getInstance()->GetTree()->Branch("Ex",&Ex,"Ex/D");
+ RootOutput::getInstance()->GetTree()->Branch("Ex",&Ex);
+ //RootOutput::getInstance()->GetTree()->Branch("EDC",&EDC,"EDC/D");
+ RootOutput::getInstance()->GetTree()->Branch("EDC",&EDC);
+ RootOutput::getInstance()->GetTree()->Branch("AddBack_EDC",&AddBack_EDC);
+ RootOutput::getInstance()->GetTree()->Branch("AddBack_EDC2",&AddBack_EDC2);
+ RootOutput::getInstance()->GetTree()->Branch("EAgata",&EAgata,"EAgata/D");
+ RootOutput::getInstance()->GetTree()->Branch("ELab",&ELab);
+ RootOutput::getInstance()->GetTree()->Branch("Ecm",&Ecm);
+ RootOutput::getInstance()->GetTree()->Branch("RawEnergy",&RawEnergy); // CPx ADDITION
+ RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab);
+ RootOutput::getInstance()->GetTree()->Branch("PhiLab",&PhiLab);
+ RootOutput::getInstance()->GetTree()->Branch("ThetaCM",&ThetaCM);
+ RootOutput::getInstance()->GetTree()->Branch("Run",&Run,"Run/I");
+ RootOutput::getInstance()->GetTree()->Branch("X",&X);
+ RootOutput::getInstance()->GetTree()->Branch("Y",&Y);
+ RootOutput::getInstance()->GetTree()->Branch("Z",&Z);
+ RootOutput::getInstance()->GetTree()->Branch("dE",&dE,"dE/D");
+ RootOutput::getInstance()->GetTree()->Branch("MG_T",MG_T);
+ RootOutput::getInstance()->GetTree()->Branch("MG_E",MG_E);
+ RootOutput::getInstance()->GetTree()->Branch("MG_X",MG_X);
+ RootOutput::getInstance()->GetTree()->Branch("MG_Y",MG_Y);
+ RootOutput::getInstance()->GetTree()->Branch("MG_D",MG_D);
+
+ // Vamos
+ RootOutput::getInstance()->GetTree()->Branch("LTS",<S,"LTS/l");
/* RootOutput::getInstance()->GetTree()->Branch("T_FPMW_CATS1",&T_FPMW_CATS1,"T_FPMW_CATS1/s");
RootOutput::getInstance()->GetTree()->Branch("T_FPMW_CATS2_C",&T_FPMW_CATS2_C,"T_FPMW_CATS2_C/F");
RootOutput::getInstance()->GetTree()->Branch("T_FPMW_HF",&T_FPMW_HF,"T_FPMW_CATS1/s");
@@ -428,12 +613,12 @@ void Analysis::TreatEvent(){
RootOutput::getInstance()->GetTree()->Branch ("EWIRE_2_1", &EWIRE_2_1, "EWIRE_2_1/s");
RootOutput::getInstance()->GetTree()->Branch ("EWIRE_2_2", &EWIRE_2_2, "EWIRE_2_2/s");
*/
- RootOutput::getInstance()->GetTree()->Branch ("MW_Nr", &MW_Nr, "MW_Nr/I");
- RootOutput::getInstance()->GetTree()->Branch ("MW_T", MW_T, "MW_T[MW_Nr]/F");
- RootOutput::getInstance()->GetTree()->Branch ("MW_N", MW_N, "MW_N[MW_Nr]/s");
- RootOutput::getInstance()->GetTree()->Branch ("MWT", MWT, "MWT[MW_Nr]/F");
+ RootOutput::getInstance()->GetTree()->Branch ("MW_Nr", &MW_Nr, "MW_Nr/I");
+ RootOutput::getInstance()->GetTree()->Branch ("MW_T", MW_T, "MW_T[MW_Nr]/F");
+ RootOutput::getInstance()->GetTree()->Branch ("MW_N", MW_N, "MW_N[MW_Nr]/s");
+ RootOutput::getInstance()->GetTree()->Branch ("MWT", MWT, "MWT[MW_Nr]/F");
- RootOutput::getInstance()->GetTree()->Branch ("AGAVA_VAMOSTS", &AGAVA_VAMOSTS, "AGAVA_VAMOSTS/l");
+ RootOutput::getInstance()->GetTree()->Branch ("AGAVA_VAMOSTS", &AGAVA_VAMOSTS, "AGAVA_VAMOSTS/l");
/* RootOutput::getInstance()->GetTree()->Branch("mT",&mT,"mT/D");
RootOutput::getInstance()->GetTree()->Branch("mV",&mV,"mV/D");
RootOutput::getInstance()->GetTree()->Branch("mD",&mD,"mD/D");
@@ -461,23 +646,32 @@ void Analysis::TreatEvent(){
RootOutput::getInstance()->GetTree()->Branch("trackCrystalID",trackCrystalID,"trackCrystalID[nbTrack]/I");
*/
- //Agata Addback
- RootOutput::getInstance()->GetTree()->Branch("nbAdd",&nbAdd,"nbAdd/I");
- RootOutput::getInstance()->GetTree()->Branch("TSHit",&TSHit,"TSHit/l");
- RootOutput::getInstance()->GetTree()->Branch("AddE",AddE,"AddE[nbAdd]/F");
- RootOutput::getInstance()->GetTree()->Branch("AddX",AddX,"AddX[nbAdd]/F");
- RootOutput::getInstance()->GetTree()->Branch("AddY",AddY,"AddY[nbAdd]/F");
- RootOutput::getInstance()->GetTree()->Branch("AddZ",AddZ,"AddZ[nbAdd]/F");
+ //Agata Addback
+ RootOutput::getInstance()->GetTree()->Branch("nbAdd",&nbAdd,"nbAdd/I");
+ RootOutput::getInstance()->GetTree()->Branch("TSHit",&TSHit,"TSHit/l");
+ RootOutput::getInstance()->GetTree()->Branch("AddE",AddE,"AddE[nbAdd]/F");
+ RootOutput::getInstance()->GetTree()->Branch("AddX",AddX,"AddX[nbAdd]/F");
+ RootOutput::getInstance()->GetTree()->Branch("AddY",AddY,"AddY[nbAdd]/F");
+ RootOutput::getInstance()->GetTree()->Branch("AddZ",AddZ,"AddZ[nbAdd]/F");
+
+ if(isSim){
+ RootOutput::getInstance()->GetTree()->Branch("OriginalELab",
+ &OriginalELab,"OriginalELab/D");
+ RootOutput::getInstance()->GetTree()->Branch("OriginalThetaLab",
+ &OriginalThetaLab,"OriginalThetaLab/D");
+ RootOutput::getInstance()->GetTree()->Branch("BeamEnergy",
+ &BeamEnergy,"BeamEnergy/D");
}
+}
- ////////////////////////////////////////////////////////////////////////////////
- void Analysis::InitInputBranch(){
- SetBranchStatus();
- // RootInput:: getInstance()->GetChain()->SetBranchAddress("GATCONF",&vGATCONF);
- //
- // Vamos
- RootInput::getInstance()->GetChain()->SetBranchAddress("LTS",<S);
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::InitInputBranch(){
+ SetBranchStatus();
+ // RootInput:: getInstance()->GetChain()->SetBranchAddress("GATCONF",&vGATCONF);
+ //
+ // Vamos
+ RootInput::getInstance()->GetChain()->SetBranchAddress("LTS",<S);
/* RootInput::getInstance()->GetChain()->SetBranchAddress("T_FPMW_CATS1",&T_FPMW_CATS1);
RootInput::getInstance()->GetChain()->SetBranchAddress("T_FPMW_CATS2_C",&T_FPMW_CATS2_C);
RootInput::getInstance()->GetChain()->SetBranchAddress("T_FPMW_HF_C",&T_FPMW_HF_C);
@@ -500,12 +694,12 @@ void Analysis::TreatEvent(){
RootInput::getInstance()->GetChain()->SetBranchAddress("EWIRE_2_1", &EWIRE_2_1);
RootInput::getInstance()->GetChain()->SetBranchAddress("EWIRE_2_2", &EWIRE_2_2);
*/
- RootInput::getInstance()->GetChain()->SetBranchAddress("MWTVM", &MW_Nr);
- RootInput::getInstance()->GetChain()->SetBranchAddress("MWTV", MW_T);
- RootInput::getInstance()->GetChain()->SetBranchAddress("MWTVN", MW_N);
- RootInput::getInstance()->GetChain()->SetBranchAddress("AGAVA_VAMOSTS", &AGAVA_VAMOSTS);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("MWTVM", &MW_Nr);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("MWTV", MW_T);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("MWTVN", MW_N);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("AGAVA_VAMOSTS", &AGAVA_VAMOSTS);
- // Agata
+ // Agata
/*
RootInput::getInstance()->GetChain()->SetBranchAddress("TStrack",&TStrack);
RootInput::getInstance()->GetChain()->SetBranchAddress("nbTrack",&nbTrack);
@@ -516,17 +710,27 @@ void Analysis::TreatEvent(){
RootInput::getInstance()->GetChain()->SetBranchAddress("trackT",trackT);
RootInput::getInstance()->GetChain()->SetBranchAddress("trackCrystalID",trackCrystalID);
*/
- RootInput::getInstance()->GetChain()->SetBranchAddress("TSHit",&TSHit);
- RootInput::getInstance()->GetChain()->SetBranchAddress("nbAdd",&nbAdd);
- RootInput::getInstance()->GetChain()->SetBranchAddress("AddE",AddE);
- RootInput::getInstance()->GetChain()->SetBranchAddress("AddX",AddX);
- RootInput::getInstance()->GetChain()->SetBranchAddress("AddY",AddY);
- RootInput::getInstance()->GetChain()->SetBranchAddress("AddZ",AddZ);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("TSHit",&TSHit);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("nbAdd",&nbAdd);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("AddE",AddE);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("AddX",AddX);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("AddY",AddY);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("AddZ",AddZ);
+
+ if(isSim){
+ //RootInput:: getInstance()->GetChain()->SetBranchStatus("InitialConditions",true );
+ //RootInput:: getInstance()->GetChain()->SetBranchStatus("fIC_*",true );
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("InitialConditions",&Initial);
+ //RootInput:: getInstance()->GetChain()->SetBranchStatus("ReactionConditions",true );
+ //RootInput:: getInstance()->GetChain()->SetBranchStatus("fRC_*",true );
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("ReactionConditions",
+ &ReactionConditions);
}
+}
- void Analysis::SetBranchStatus(){
- // Set Branch status
- RootInput::getInstance()->GetChain()->SetBranchStatus("LTS",true);
+void Analysis::SetBranchStatus(){
+ // Set Branch status
+ RootInput::getInstance()->GetChain()->SetBranchStatus("LTS",true);
/* RootInput::getInstance()->GetChain()->SetBranchStatus("T_FPMW_CATS1",true);
RootInput::getInstance()->GetChain()->SetBranchStatus("T_FPMW_CATS2_C",true);
RootInput::getInstance()->GetChain()->SetBranchStatus("T_FPMW_HF",true );
@@ -547,15 +751,15 @@ void Analysis::TreatEvent(){
RootInput::getInstance()->GetChain()->SetBranchStatus("EWIRE_2_1",true );
RootInput::getInstance()->GetChain()->SetBranchStatus("EWIRE_2_2",true );
*/
- RootInput::getInstance()->GetChain()->SetBranchStatus("MW_Nr",true );
- RootInput::getInstance()->GetChain()->SetBranchStatus("MW_T",true );
+ RootInput::getInstance()->GetChain()->SetBranchStatus("MW_Nr",true );
+ RootInput::getInstance()->GetChain()->SetBranchStatus("MW_T",true );
- RootInput::getInstance()->GetChain()->SetBranchStatus("MWTV",true );
- RootInput::getInstance()->GetChain()->SetBranchStatus("MWTVN",true );
- RootInput::getInstance()->GetChain()->SetBranchStatus("MWTVM",true );
+ RootInput::getInstance()->GetChain()->SetBranchStatus("MWTV",true );
+ RootInput::getInstance()->GetChain()->SetBranchStatus("MWTVN",true );
+ RootInput::getInstance()->GetChain()->SetBranchStatus("MWTVM",true );
- RootInput::getInstance()->GetChain()->SetBranchStatus("AGAVA_VAMOSTS",true );
- // Agata
+ RootInput::getInstance()->GetChain()->SetBranchStatus("AGAVA_VAMOSTS",true );
+ // Agata
/*
RootInput::getInstance()->GetChain()->SetBranchStatus("TStrack",true );
RootInput::getInstance()->GetChain()->SetBranchStatus("nbTrack",true );
@@ -566,73 +770,82 @@ void Analysis::TreatEvent(){
RootInput::getInstance()->GetChain()->SetBranchStatus("trackT",true );
RootInput::getInstance()->GetChain()->SetBranchStatus("trackCrystalID",true);
*/
- RootInput::getInstance()->GetChain()->SetBranchStatus("TSHit",true );
- RootInput::getInstance()->GetChain()->SetBranchStatus("nbAdd",true );
- RootInput::getInstance()->GetChain()->SetBranchStatus("AddE",true );
- RootInput::getInstance()->GetChain()->SetBranchStatus("AddX",true );
- RootInput::getInstance()->GetChain()->SetBranchStatus("AddY",true );
- RootInput::getInstance()->GetChain()->SetBranchStatus("AddZ",true );
+ RootInput::getInstance()->GetChain()->SetBranchStatus("TSHit",true );
+ RootInput::getInstance()->GetChain()->SetBranchStatus("nbAdd",true );
+ RootInput::getInstance()->GetChain()->SetBranchStatus("AddE",true );
+ RootInput::getInstance()->GetChain()->SetBranchStatus("AddX",true );
+ RootInput::getInstance()->GetChain()->SetBranchStatus("AddY",true );
+ RootInput::getInstance()->GetChain()->SetBranchStatus("AddZ",true );
+
+ if(isSim){
+ RootInput:: getInstance()->GetChain()->SetBranchStatus("InitialConditions",true );
+ RootInput:: getInstance()->GetChain()->SetBranchStatus("fIC_*",true );
+ RootInput:: getInstance()->GetChain()->SetBranchStatus("ReactionConditions",true );
+ RootInput:: getInstance()->GetChain()->SetBranchStatus("fRC_*",true );
}
+}
- ////////////////////////////////////////////////////////////////////////////////
- void Analysis::ReInitValue(){
- Ex.clear();
- Ecm.clear();
- AddBack_EDC.clear();
- EAgata = -1000;
- ELab.clear();
- RawEnergy.clear();
- ThetaLab.clear();
- PhiLab.clear();
- ThetaCM.clear();
- //relative_angle = -1000;
- //sum_angle = -1000;
- X.clear();
- Y.clear();
- Z.clear();
- dE= -1000;
- ParticleMult = 0;
- GammaMult = 0;
- DetectorNumber = 0;
- ThetaNormalTarget = 0;
- ThetaM2Surface = 0;
- ThetaMGSurface = 0;
- Si_E_M2 = 0;
- CsI_E_M2 = 0;
- Energy = 0;
- ThetaGDSurface = 0;
- BeamDirection = TVector3(0,0,1);
- elab_tmp = 0;
- thetalab_tmp = 0;
- philab_tmp = 0;
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::ReInitValue(){
+ Ex.clear();
+ Ecm.clear();
+ AddBack_EDC.clear();
+ AddBack_EDC2.clear();
+ AddBack_EDC2.push_back(-1.0); //offset by 1
+ EAgata = -1000;
+ ELab.clear();
+ RawEnergy.clear();
+ ThetaLab.clear();
+ PhiLab.clear();
+ ThetaCM.clear();
+ //relative_angle = -1000;
+ //sum_angle = -1000;
+ X.clear();
+ Y.clear();
+ Z.clear();
+ dE= -1000;
+ ParticleMult = 0;
+ GammaMult = 0;
+ DetectorNumber = 0;
+ ThetaNormalTarget = 0;
+ ThetaM2Surface = 0;
+ ThetaMGSurface = 0;
+ Si_E_M2 = 0;
+ CsI_E_M2 = 0;
+ Energy = 0;
+ ThetaGDSurface = 0;
+ BeamDirection = TVector3(0,0,1);
+ elab_tmp = 0;
+ thetalab_tmp = 0;
+ philab_tmp = 0;
- MG_T=-1000;
- MG_E=-1000;
- MG_X=-1000;
- MG_Y=-1000;
- MG_D=-1000;
+ MG_T=-1000;
+ MG_E=-1000;
+ MG_X=-1000;
+ MG_Y=-1000;
+ MG_D=-1000;
- //EDC= -1000;
- EDC.clear();
- }
+ //EDC= -1000;
+ EDC.clear();
+}
- ////////////////////////////////////////////////////////////////////////////////
- // Construct Method to be pass to the AnalysisFactory //
- ////////////////////////////////////////////////////////////////////////////////
- NPL::VAnalysis* Analysis::Construct(){
- return (NPL::VAnalysis*) new Analysis();
- }
+////////////////////////////////////////////////////////////////////////////////
+// Construct Method to be pass to the AnalysisFactory //
+////////////////////////////////////////////////////////////////////////////////
+NPL::VAnalysis* Analysis::Construct(){
+ return (NPL::VAnalysis*) new Analysis();
+}
- ////////////////////////////////////////////////////////////////////////////////
- // Registering the construct method to the factory //
- ////////////////////////////////////////////////////////////////////////////////
- extern "C"{
- class proxy_analysis{
- public:
- proxy_analysis(){
- NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct);
- }
- };
+////////////////////////////////////////////////////////////////////////////////
+// Registering the construct method to the factory //
+////////////////////////////////////////////////////////////////////////////////
+extern "C"{
+ class proxy_analysis{
+ public:
+ proxy_analysis(){
+ NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct);
+ }
+ };
proxy_analysis p_analysis;
}
diff --git a/Projects/e793s/Analysis.h b/Projects/e793s/Analysis.h
index 3d3ed1bbc99fc02287a55c988593b85b56f06f0a..cbff593d09cc94c34d79ee4fc1b569d5dfe93d9e 100755
--- a/Projects/e793s/Analysis.h
+++ b/Projects/e793s/Analysis.h
@@ -28,6 +28,9 @@
#include"NPBeam.h"
#include"RootOutput.h"
#include"RootInput.h"
+#include "TF1.h"
+#include "TInitialConditions.h"
+#include "TReactionConditions.h"
#include "TMust2Physics.h"
#include "TMugastPhysics.h"
#include "TCATSPhysics.h"
@@ -38,6 +41,11 @@
#include <TMath.h>
#include <bitset>
+auto ThetaCM_emmitted = new TH1F("ThetaCM_emmitted","ThetaCM_emmitted",180,0,180);
+auto ThetaCM_detected = new TH1F("ThetaCM_detected","ThetaCM_detected",180,0,180);
+auto ThetaLab_emmitted = new TH1F("ThetaLab_emmitted","ThetaLab_emmitted",180,0,180);
+auto ThetaLab_detected = new TH1F("ThetaLab_detected","ThetaLab_detected",180,0,180);
+
class Analysis: public NPL::VAnalysis{
public:
Analysis();
@@ -61,6 +69,7 @@ class Analysis: public NPL::VAnalysis{
//double EDC;
std::vector<double> EDC;
vector<double> AddBack_EDC;
+ vector<double> AddBack_EDC2;
double EAgata;
std::vector<double> ELab;
std::vector<double> Ex;
@@ -69,7 +78,12 @@ class Analysis: public NPL::VAnalysis{
std::vector<double> ThetaLab;
std::vector<double> PhiLab;
std::vector<double> ThetaCM;
-
+
+ //For use in simulations
+ double OriginalELab;
+ double OriginalThetaLab;
+ double BeamEnergy;
+
NPL::Reaction reaction;
// Energy loss table: the G4Table are generated by the simulation
NPL::EnergyLoss LightTarget;
@@ -217,6 +231,12 @@ class Analysis: public NPL::VAnalysis{
//TCATSPhysics* CATS;
TModularLeafPhysics* ML;
+ //Simulation Stuff
+ bool isSim;
+ bool writetoscreen;
+ TInitialConditions* Initial;
+ TReactionConditions* ReactionConditions;
+
// Beam object
NPL::Beam* Beam;