diff --git a/Projects/e775s/analysis/Analysis.cxx b/Projects/e775s/analysis/Analysis.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..120b33601acb9415c149b28ce3eb305d7a0583d3
--- /dev/null
+++ b/Projects/e775s/analysis/Analysis.cxx
@@ -0,0 +1,1855 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2014 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Adrien MATTA contact address: a.matta@surrey.ac.uk *
+ * Edited for e793s *
+ * *
+ * Creation Date : *
+ * Last update : May 2021 *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * Class describing the property of an Analysis object *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+#include<iostream>
+#include<string>
+
+#include<fstream>
+
+using namespace std;
+
+#include "Analysis.h"
+#include "NPFunction.h"
+#include "NPAnalysisFactory.h"
+#include "NPDetectorManager.h"
+#include "NPOptionManager.h"
+#include "DefineColours.h"
+
+////////////////////////////////////////////////////////////////////////////////
+Analysis::Analysis(){
+}
+
+////////////////////////////////////////////////////////////////////////////////
+Analysis::~Analysis(){
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::Init() {
+cout << "start_Init" << endl;
+ ///////////////////////////////////////////////////////////////////////////////
+
+ if(NPOptionManager::getInstance()->HasDefinition("Sim")){
+ cout << " == == == == SIMULATION == == == ==" << endl;
+ isSim=true;
+ isPhaseSpace=false;
+ } else if (NPOptionManager::getInstance()->HasDefinition("Exp")) {
+ cout << " == == == == EXPERIMENT == == == ==" << endl;
+ isSim=false;
+ isPhaseSpace=false;
+ } else if (NPOptionManager::getInstance()->HasDefinition("PSp")) {
+ cout << " == == == == PHASE SPACE == == == ==" << endl;
+ isSim=false;
+ //isSim=true;
+ isPhaseSpace=true;
+ } else {
+ cout << " == == == == DEFAULT (EXP) == == == ==" << endl;
+ isSim=false;
+ isPhaseSpace=false;
+ }
+
+ if(NPOptionManager::getInstance()->HasDefinition("Au")){
+ cout << " -- -- -- -- Au Target -- -- -- --" << endl;
+ runningAuTarget=true;
+ } else if(NPOptionManager::getInstance()->HasDefinition("CD2")){
+ cout << " -- -- -- -- CD2 Target -- -- -- --" << endl;
+ runningAuTarget=false;
+ } else {
+ cout << " ERROR: NO TARGET SPECIFIED BY --DEFINE !" << endl;
+ cout << " ERROR: NO TARGET SPECIFIED BY --DEFINE !" << endl;
+ cout << " ERROR: NO TARGET SPECIFIED BY --DEFINE !" << endl;
+ }
+
+// double simCentroid=-10.;
+// double simSigma = 0.104;
+// if(isSim){
+// if(NPOptionManager::getInstance()->HasDefinition("S0000")){ simCentroid = 0.000; }
+// else if(NPOptionManager::getInstance()->HasDefinition("S1675")){ simCentroid = 1.675; }
+// else if(NPOptionManager::getInstance()->HasDefinition("S3572")){ simCentroid = 3.572; }
+// else if(NPOptionManager::getInstance()->HasDefinition("S4071")){ simCentroid = 4.701; }
+// else if(NPOptionManager::getInstance()->HasDefinition("S5004")){ simCentroid = 5.004; }
+// else if(NPOptionManager::getInstance()->HasDefinition("S5228")){ simCentroid = 5.228; }
+// else if(NPOptionManager::getInstance()->HasDefinition("S5629")){ simCentroid = 5.629; }
+// else if(NPOptionManager::getInstance()->HasDefinition("S7622")){ simCentroid = 7.622; }
+// else if(NPOptionManager::getInstance()->HasDefinition("S7754")){ simCentroid = 7.754; }
+// else if(NPOptionManager::getInstance()->HasDefinition("S8554")){ simCentroid = 8.554; }
+// else if(NPOptionManager::getInstance()->HasDefinition("S8962")){ simCentroid = 8.962; }
+// else if(NPOptionManager::getInstance()->HasDefinition("S9770")){ simCentroid = 9.770; }
+// else if(NPOptionManager::getInstance()->HasDefinition("S3575")){
+// simCentroid = 3.575;
+// cout << "TESTING TESTING TESTING" << endl;
+// cout << "TESTING TESTING TESTING" << endl;
+// cout << "TESTING TESTING TESTING" << endl;
+// cout << "TESTING TESTING TESTING" << endl;
+// cout << "TESTING TESTING TESTING" << endl;
+// cout << "TESTING TESTING TESTING" << endl;
+// cout << "TESTING TESTING TESTING" << endl;
+// cout << "TESTING TESTING TESTING" << endl;
+// cout << "TESTING TESTING TESTING" << endl;
+// cout << "TESTING TESTING TESTING" << endl;
+// cout << "TESTING TESTING TESTING" << endl;
+// cout << "TESTING TESTING TESTING" << endl;
+// cout << "TESTING TESTING TESTING" << endl;
+// cout << "TESTING TESTING TESTING" << endl;
+// cout << "TESTING TESTING TESTING" << endl;
+// }
+// else {
+// cout << " ADD STATE ENERGY DEFINITION FOR CLEANUP OF SOLID ANGLE SIMULATION!!! " << endl;
+// cout << " ADD STATE ENERGY DEFINITION FOR CLEANUP OF SOLID ANGLE SIMULATION!!! " << endl;
+// cout << " ADD STATE ENERGY DEFINITION FOR CLEANUP OF SOLID ANGLE SIMULATION!!! " << endl;
+// cout << " ADD STATE ENERGY DEFINITION FOR CLEANUP OF SOLID ANGLE SIMULATION!!! " << endl;
+// cout << " ADD STATE ENERGY DEFINITION FOR CLEANUP OF SOLID ANGLE SIMULATION!!! " << endl;
+// cout << " ADD STATE ENERGY DEFINITION FOR CLEANUP OF SOLID ANGLE SIMULATION!!! " << endl;
+// cout << " ADD STATE ENERGY DEFINITION FOR CLEANUP OF SOLID ANGLE SIMULATION!!! " << endl;
+// }
+// }
+// cout << simCentroid << " +- " << simSigma << endl;
+
+ agata_zShift=(49.0-0.63)*mm;
+ //agata_zShift=51*mm;
+ //BrhoRef=0.65;
+
+ if(isSim && !isPhaseSpace){
+cout << "here_InIsSimLoop" << endl;
+ Initial = new TInitialConditions();
+ ReactionConditions = new TReactionConditions();
+ RootInput::getInstance()->GetChain()->SetBranchAddress("InitialConditions",&Initial);
+
+ // Initilize MGX histograms, because must be inside function
+ string base1 = "ThetaCM_detected_MG";
+ string base2 = "ThetaLab_detected_MG";
+ string base1unb = "ThetaCM_detected_MG_Unb";
+ string base2unb = "ThetaLab_detected_MG_Unb";
+
+ for(int i=0; i<7; i++){
+ int j=i+1;
+ string name1 = base1 + to_string(j);
+ string name2 = base2 + to_string(j);
+ string name1unb = base1unb + to_string(j);
+ string name2unb = base2unb + to_string(j);
+
+ ThetaCM_detected_MGX[i]
+ = new TH1F( name1.c_str() , name1.c_str() , NumThetaAngleBins, 0, 180 );
+ ThetaLab_detected_MGX[i]
+ = new TH1F( name2.c_str() , name2.c_str() , NumThetaAngleBins, 0, 180 );
+ ThetaCM_detected_MGX_Unb[i]
+ = new TH1F( name1unb.c_str(), name1unb.c_str(), NumThetaAngleBins, 0, 180 );
+ ThetaLab_detected_MGX_Unb[i]
+ = new TH1F( name2unb.c_str(), name2unb.c_str(), NumThetaAngleBins, 0, 180 );
+ }
+ }
+
+ // initialize input and output branches
+ InitOutputBranch();
+ InitInputBranch();
+
+//M2 = (TMust2Physics*) m_DetectorManager -> GetDetector("MUST2");
+ MG = (TMugastPhysics*) m_DetectorManager -> GetDetector("Mugast");
+ ML = (TModularLeafPhysics*) m_DetectorManager -> GetDetector("ModularLeaf");
+// CATS = (TCATSPhysics*) m_DetectorManager->GetDetector("CATSDetector");
+
+ // get reaction information
+ reaction.ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
+ OriginalBeamEnergy = reaction.GetBeamEnergy();
+ reaction.Print(); //TESTING PARSER
+
+ // get beam position from .reaction file
+ Beam = (NPL::Beam*) reaction.GetParticle1();
+ XBeam = Beam->GetMeanX();
+ YBeam = Beam->GetMeanY();
+
+ cout << " ---------------- Beam Position ---------------- " << endl;
+ cout << " \tX = " << XBeam << " mm\tY = " << YBeam << " mm" << endl;
+ cout << " ----------------------------------------------- " << endl;
+
+ // target thickness
+ TargetThickness = m_DetectorManager->GetTargetThickness();
+ string TargetMaterial = m_DetectorManager->GetTargetMaterial();
+
+ AuThickness = m_DetectorManager->GetBackThickness();
+ string BackMaterial = m_DetectorManager->GetBackMaterial();
+
+ cout << " ---------------- Target Values ---------------- " << endl;
+ cout << " CD2: \tZ = " << m_DetectorManager->GetTargetZ() << " mm\tT = " << TargetThickness << " mm" << endl;
+ cout << " Au: \tZ = " << m_DetectorManager->GetTargetZ() << " mm\tT = " << AuThickness << " mm" << endl;
+ cout << " ----------------------------------------------- " << endl;
+
+ // energy losses
+ string light = NPL::ChangeNameToG4Standard(reaction.GetParticle3()->GetName());
+ string beam = NPL::ChangeNameToG4Standard(reaction.GetParticle1()->GetName());
+ string heavy = NPL::ChangeNameToG4Standard(reaction.GetParticle4()->GetName());
+
+ if(isPhaseSpace){
+ //then doesn't pull definitions properly
+ // TEMP FIX::
+ light = "proton";
+ beam = "O19";
+ heavy = "O20";
+ }
+
+ LightTarget = NPL::EnergyLoss(light+"_"+TargetMaterial+".G4table","G4Table",100 );
+ LightAl = NPL::EnergyLoss(light+"_Al.G4table" ,"G4Table",100);
+ LightSi = NPL::EnergyLoss(light+"_Si.G4table" ,"G4Table",100);
+ LightAu = NPL::EnergyLoss(light+"_Au.G4table" ,"G4Table",100);
+//HeavyAu = NPL::EnergyLoss(heavy+"_Au.G4table" ,"G4Table",100);
+ BeamTarget = NPL::EnergyLoss(beam+"_"+TargetMaterial+".G4table","G4Table",100);
+//HeavyTarget = NPL::EnergyLoss(heavy+"_"+TargetMaterial+".G4table","G4Table",100);
+
+ FinalBeamEnergy = BeamTarget.Slow(OriginalBeamEnergy, 0.5*TargetThickness, 0);
+ reaction.SetBeamEnergy(FinalBeamEnergy);
+
+ cout << "\033[91m Running for reaction "
+ << reaction.GetParticle1()->GetName() << "("
+ << reaction.GetParticle2()->GetName() << ","
+ << reaction.GetParticle3()->GetName() << ")"
+ << reaction.GetParticle4()->GetName() << endl;
+
+ cout << "\033[36m Beam energy at middle of CD2: " << FinalBeamEnergy << "\033[37m"<< endl;
+
+ // initialize various parameters
+ Rand = TRandom3();
+ ParticleMult = 0;
+ GammaMult = 0;
+ DetectorNumber = 0;
+ ThetaNormalTarget = 0;
+ ThetaM2Surface = 0;
+ ThetaMGSurface = 0;
+ Si_E_M2 = 0;
+ CsI_E_M2 = 0;
+ Energy = 0;
+ ThetaGDSurface = 0;
+ elab_tmp = 0;
+ thetalab_tmp = 0;
+ philab_tmp = 0;
+ ThetaGDSurface = 0;
+ recoilBeta = 0;
+ X.clear();
+ Y.clear();
+ Z.clear();
+ dE = 0;
+ BeamDirection = TVector3(0,0,1);
+ //nbTrack=0;
+ //nbHits=0;
+ count=0;
+ AHeavy=reaction.GetParticle4()->GetA();
+ ALight=reaction.GetParticle3()->GetA();
+ MHeavy=reaction.GetParticle4()->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();
+cout << "end_Init" << endl;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::TreatEvent(){
+//cout << "start_TreatEvent" << endl;
+
+ // Reinitiate calculated variable
+ ReInitValue();
+ CalculateVamos();
+
+ if(isSim && !isPhaseSpace){
+ ThetaCM_emmitted->Fill(ReactionConditions->GetThetaCM());
+ ThetaLab_emmitted->Fill(ReactionConditions->GetTheta(0));
+ }
+
+ 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
+ GATCONF_Counter[i]++ ; // increment the array
+ // cout << std::bitset<16>(GATCONF_MASTER) << " " << i+1 << endl;
+ }
+ }
+
+ if(isSim && !isPhaseSpace){
+ OriginalELab = ReactionConditions->GetKineticEnergy(0);
+ OriginalThetaLab = ReactionConditions->GetTheta(0);
+ BeamEnergy = ReactionConditions->GetBeamEnergy();
+ }
+
+ TVector3 BeamDirection(0.,0.,1.);
+ BeamImpact = TVector3(XBeam,YBeam,m_DetectorManager->GetTargetZ());
+
+ ////////////////////////////////////////////////////////////////////////////
+ //////////////////////////////// LOOP on MUGAST ////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////
+ unsigned int sizeMG = MG->DSSD_E.size();
+ //for(unsigned int countMugast = 0; countMugast<sizeMG; countMugast++){
+ for(unsigned int countMugast = 0; countMugast<MG->DSSD_E.size(); countMugast++){
+
+//cout << "start_MugastEvent" << endl;
+
+
+// MOVING LATER TO ADD Ex LIMITATIONS
+// if(isSim && !isPhaseSpace){
+// int MGX = MG->TelescopeNumber[0];
+// // IGNORE ANNULAR EVENTS?
+// if(MGX==11){ break; }
+// MGX = MGX-1;
+
+// ThetaCM_detected_MG->Fill(ReactionConditions->GetThetaCM());
+// ThetaCM_detected_MGX[MGX]->Fill(ReactionConditions->GetThetaCM());
+
+// ThetaLab_detected_MG->Fill(ReactionConditions->GetTheta(0));
+// ThetaLab_detected_MGX[MGX]->Fill(ReactionConditions->GetTheta(0));
+// }
+
+ // Part 1 : Impact Angle
+ ThetaMGSurface = 0;
+ ThetaNormalTarget = 0;
+ thetalab_tmp = 0;
+ philab_tmp = 0;
+ E4_tmp = 0.0;
+
+ TVector3 HitDirection = MG->GetPositionOfInteraction(countMugast) - BeamImpact;
+ thetalab_tmp = HitDirection.Angle(BeamDirection);
+ philab_tmp = HitDirection.Phi();
+
+ X.push_back( MG -> GetPositionOfInteraction(countMugast).X());
+ Y.push_back( MG -> GetPositionOfInteraction(countMugast).Y());
+ Z.push_back( MG -> GetPositionOfInteraction(countMugast).Z());
+
+//cout << " X = " << MG -> GetPositionOfInteraction(countMugast).X() << "\t Y = " << MG -> GetPositionOfInteraction(countMugast).Y() << endl;
+//cout << MG -> GetPositionOfInteraction(countMugast).X() << "\t" << MG -> GetPositionOfInteraction(countMugast).Y() << endl;
+//cout << "STRIP X & Y \t" << MG -> DSSD_X.back() << "\t" << MG -> DSSD_Y.back() << "\t" << MG->DSSD_X.back() + MG->DSSD_Y.back() << endl;
+
+ ThetaMGSurface = HitDirection.Angle( MG -> GetTelescopeNormal(countMugast) );
+ ThetaNormalTarget = HitDirection.Angle( TVector3(0.0,0.0,1.0) ) ;
+ //ThetaNormalTarget = HitDirection.Angle( TVector3(0.0,0.0,-1.0) ) ;
+
+ MGDet.push_back(MG->TelescopeNumber[0]);
+ MGXStrip.push_back(MG->DSSD_X.back());
+ MGYStrip.push_back(MG->DSSD_Y.back());
+
+ // Part 2 : Impact Energy
+ Energy = elab_tmp = 0;
+
+ /* Below is from IZanon - presumably fixing some issue with the annular? */
+ // removed, because I don't care about annular for unbounds
+/**///if(MGDet.back()==11){
+/**/// if(MGYStrip.back()==1 || MGYStrip.back()==2 || MGYStrip.back()==15 && MGYStrip.back()==16){
+/**/// if(MGXStrip.back()>32 && MGXStrip.back()<49) Energy = MG->GetEnergyDeposit(countMugast);
+/**/// }
+/**/// else if(MGYStrip.back()>2 && MGYStrip.back()<7){
+/**/// if(MGXStrip.back()>48 && MGXStrip.back()<65) Energy = MG->GetEnergyDeposit(countMugast);
+/**/// }
+/**/// else if(MGYStrip.back()>6 && MGYStrip.back()<11){
+/**/// if(MGXStrip.back()>16 && MGXStrip.back()<50) Energy = MG->GetEnergyDeposit(countMugast);
+/**/// }
+/**/// else if(MGYStrip.back()>10 && MGYStrip.back()<15){
+/**/// if(MGXStrip.back()>0 && MGXStrip.back()<17) Energy = MG->GetEnergyDeposit(countMugast);
+/**/// }
+/**/// else Energy = 0*MeV;
+/**///}
+/**///else Energy = MG->GetEnergyDeposit(countMugast);
+ Energy = MG->GetEnergyDeposit(countMugast);
+
+ double rawe = Energy; //for parallel line unbound state calc later, new addiiton
+ RawEnergy.push_back(Energy);
+
+// if(!isSim){
+ elab_tmp = LightAl.EvaluateInitialEnergy(
+ Energy, //particle energy after Al
+ 0.4*micrometer, //thickness of Al
+ //ThetaMGSurface); //angle of impingement
+ 0.); //angle of impingement
+ ELoss_Al.push_back(Energy-elab_tmp);
+ double elab_tmp2 = elab_tmp;
+ elab_tmp = LightTarget.EvaluateInitialEnergy(
+ elab_tmp, //particle energy after leaving target
+ TargetThickness*0.5, //distance passed through target
+ ThetaNormalTarget); //angle of exit from target
+ ELoss_Target.push_back(elab_tmp2-elab_tmp);
+ ELoss.push_back(Energy-elab_tmp);
+// } else { //TESTING DIFFERENT ENERGY LOSSES IN SIMULATION
+// elab_tmp = Energy; //so I can add and remove sections
+// //elab_tmp = LightSi.EvaluateInitialEnergy(
+// // elab_tmp, //particle energy after Si
+// // 0.5*500.*micrometer, //thickness of Si
+// // ThetaMGSurface); //angle of impingement
+// //elab_tmp = LightAl.EvaluateInitialEnergy(
+// // elab_tmp, //particle energy after Al
+// // 0.4*micrometer, //thickness of Al
+// // ThetaMGSurface); //angle of impingement
+// elab_tmp = LightTarget.EvaluateInitialEnergy(
+// elab_tmp, //particle energy after leaving target
+// TargetThickness*0.5, //distance passed through target
+// ThetaNormalTarget); //angle of exit from target
+// }
+
+ ELab.push_back(elab_tmp);
+ hasAuTarget.push_back(runningAuTarget);
+
+///**/if(MG->TelescopeNumber.back()==1) thetalab_tmp = (thetalab_tmp/deg + 0.5)*deg;
+///**/if(MG->TelescopeNumber.back()==2) thetalab_tmp = (thetalab_tmp/deg + 0.5)*deg;
+///**/if(MG->TelescopeNumber.back()==3) thetalab_tmp = (thetalab_tmp/deg + 0.5)*deg;
+///**/if(MG->TelescopeNumber.back()==4) thetalab_tmp = (thetalab_tmp/deg - 0.3)*deg;
+///**/if(MG->TelescopeNumber.back()==5) thetalab_tmp = (thetalab_tmp/deg - 0.5)*deg;
+///**/if(MG->TelescopeNumber.back()==6) thetalab_tmp = (thetalab_tmp/deg - 0.5)*deg;
+///**/if(MG->TelescopeNumber.back()==7) thetalab_tmp = (thetalab_tmp/deg - 0.7)*deg;
+
+// // Part 3 : Excitation Energy Calculation
+// double slope11[16] = {1.02215,0.99258,1.,1.,1.,1.,1.,1.,1.,1., 1.02094, 1.023497, 1.01173, 1.02576, 1.04290, 1.01356};
+// double inter11[16] = {0.17232, 0.72019,0.,0.,0.,0.,0.,0.,0.,0., 0.25035, 0.00398, 0.13719, 0.02731, 0.22406, 0.18893};
+// Ex.push_back(reaction.ReconstructRelativistic(elab_tmp, thetalab_tmp, philab_tmp));
+// if(MGDet.back()==11){
+// ExCalib.push_back(Ex.back()*slope11[MGYStrip.back()-1]+inter11[MGYStrip.back()-1]);
+// }
+// else if(MGDet.back()==3) ExCalib.push_back(Ex.back()-0.08*MeV);
+// else if(MGDet.back()==5) ExCalib.push_back(Ex.back()+0.03*MeV);
+// else if(MGDet.back()==6) ExCalib.push_back(Ex.back()-0.02*MeV);
+// else ExCalib.push_back(Ex.back());
+
+ //Ex.push_back(reaction.ReconstructRelativistic(elab_tmp, thetalab_tmp, philab_tmp));
+
+ double ex_tmp;
+ if(isSim){
+ //No adjustment to Ex
+ ex_tmp = reaction.ReconstructRelativistic(elab_tmp, thetalab_tmp);
+ } else {
+ //Adjustment to Ex, using straight line calibration
+ ex_tmp = (reaction.ReconstructRelativistic(elab_tmp, thetalab_tmp, philab_tmp) / 0.989) - 0.0693;
+ }
+ Ex.push_back(ex_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);
+ ThetaCM.push_back(reaction.EnergyLabToThetaCM(elab_tmp, thetalab_tmp)/deg);
+
+ if(sizeMG==1){
+ MG_T = MG->DSSD_T[0];
+ MG_E = MG->DSSD_E[0];
+ MG_X = MG->DSSD_X[0];
+ MG_Y = MG->DSSD_Y[0];
+ MG_D = MG->TelescopeNumber[0];
+ }
+
+ if(isSim && !isPhaseSpace){
+ int MGX = MG->TelescopeNumber[0];
+ // IGNORE ANNULAR EVENTS?
+ if(MGX==11){ break; }
+ MGX = MGX-1;
+
+ // REJECT POOR EVENTS (PUNCHTHROUGH ETC.)
+ //if(abs(ex_tmp - simCentroid)<3*simSigma){
+ ThetaCM_detected_MG->Fill(ReactionConditions->GetThetaCM());
+ ThetaCM_detected_MGX[MGX]->Fill(ReactionConditions->GetThetaCM());
+
+ ThetaLab_detected_MG->Fill(ReactionConditions->GetTheta(0));
+ ThetaLab_detected_MGX[MGX]->Fill(ReactionConditions->GetTheta(0));
+ //}
+ }
+
+// if(isSim && !isPhaseSpace){
+// if(reaction.ReconstructRelativistic(elab_tmp,thetalab_tmp, philab_tmp)<5.9
+// && elab_tmp>1.5
+// && thetalab_tmp/deg>104.
+// && elab_tmp<((0.07*(thetalab_tmp/deg))-5.39)
+// && elab_tmp>((0.07*(thetalab_tmp/deg))-6.18))
+// {
+// ThetaCM_detected_MG_Unb->Fill(ReactionConditions->GetThetaCM());
+// ThetaLab_detected_MG_Unb->Fill(ReactionConditions->GetTheta(0));
+// }
+// }
+
+ double HeavyThetaLab_tmp = reaction.GetEnergyImpulsionLab_4().Theta();
+ double HeavyPhiLab_tmp = reaction.GetEnergyImpulsionLab_4().Phi();
+ HeavyThetaLab.push_back(HeavyThetaLab_tmp/deg);
+ HeavyPhiLab.push_back(HeavyPhiLab_tmp/deg);
+
+ auto tmpHeavy = new TVector3(1,1,1);
+ tmpHeavy->SetTheta(HeavyThetaLab_tmp*(180./M_PI));
+ tmpHeavy->SetPhi(HeavyPhiLab_tmp*(180./M_PI));
+ tmpHeavy->SetMag(WindowDistance/cos(HeavyThetaLab_tmp));//cos=A/H, using average of MUST2 corners as 'Z'
+//HeavyPosXAtWindow.push_back(tmpHeavy->X());
+//HeavyPosYAtWindow.push_back(tmpHeavy->Y());
+//HeavyAngX.push_back(atan(tmpHeavy->X()/WindowDistance));
+//HeavyAngY.push_back(atan(tmpHeavy->Y()/WindowDistance));
+//HeavyPosXAtMUST2.push_back(MUST2InnerDistance*tan( atan(tmpHeavy->X()/WindowDistance) ));
+//HeavyPosYAtMUST2.push_back(MUST2InnerDistance*tan( atan(tmpHeavy->Y()/WindowDistance) ));
+
+// cout << "========= MG ===========" << endl;
+// cout << tmpHeavy->X() << "\t"
+// << (atan(tmpHeavy->X()/WindowDistance)) << endl;
+// cout << tmpHeavy->Y() << "\t"
+// << (atan(tmpHeavy->Y()/WindowDistance)) << endl;
+
+
+
+/**/short OffToF[12] = {0,-1800,-1400,-1450,-2100,-3580,-1300,-2820,0,0,0,0};
+/**/int NDet;
+/**/NDet = MG->TelescopeNumber.back();
+/**/for(int ii = 0; ii<12 ; ii++) ToF = T_MUGAST_FPMW + OffToF[NDet];
+
+//cout << "end_MugastEvent" << endl;
+ }//end loop Mugast
+
+
+ ////////////////////////////////////////////////////////////////////////////
+ ///////////////////////////////// LOOP on AGATA ////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////
+ //********* Taking this section basically completely from IZanon *********//
+ //********* Taking this section basically completely from IZanon *********//
+ //********* Taking this section basically completely from IZanon *********//
+
+ for(int j = 0; j < nbAdd; j++){ // all multiplicity
+
+ // Measured E
+ double Egamma=AddE[j]/1000.; // From keV to MeV
+ TVector3 GammaHit(AddX[j],AddY[j],AddZ[j]+agata_zShift);
+ //GammaHit.RotateZ(265.401*TMath::DegToRad());
+ GammaHit.RotateZ(263.348*TMath::DegToRad());
+ TVector3 GammaDirection = GammaHit-BeamImpact;
+ GammaDirection = GammaDirection.Unit();
+
+ // Construct LV in lab frame
+ TLorentzVector GammaLVM;
+ GammaLVM.SetPx(Egamma*GammaDirection.X());
+ GammaLVM.SetPy(Egamma*GammaDirection.Y());
+ GammaLVM.SetPz(Egamma*GammaDirection.Z());
+ GammaLVM.SetE(Egamma);
+
+ // Test EDC calculated taking into account the recoil direction (Irene and Alain)
+ recoilBeta = reaction.GetNucleus4()->GetEnergyImpulsion().Beta();
+
+ /* E4_tmp = reaction.GetEnergyImpulsionLab_4().E()-reaction.GetNucleus4()->Mass();
+ TLorentzVector LVector4_tmp=reaction.GetEnergyImpulsionLab_4();
+ std::cout << LVector4_tmp.Beta() << "\t\t" << reaction.GetNucleus4()->GetEnergyImpulsion().Beta() << std::endl;
+ TVector3 imp4=LVector4_tmp.Vect();
+ imp4.SetMag(sqrt((E4_tmp+reaction.GetNucleus4()->Mass())*(E4_tmp+reaction.GetNucleus4()->Mass())-reaction.GetNucleus4()->Mass()*reaction.GetNucleus4()->Mass()));
+ TLorentzVector LVector4(imp4,E4_tmp+reaction.GetNucleus4()->Mass());
+ recoilBeta = LVector4.Beta();
+ recoilDir = LVector4.Vect();*/
+
+ beta_MUGAST.push_back(recoilBeta);
+ recoilDir = reaction.GetNucleus4()->GetEnergyImpulsion().Vect();
+
+ //mRecoilDist->Fill(recoilDir.Phi()*TMath::RadToDeg(), recoilDir.Theta()*TMath::RadToDeg());
+
+ ThetaGR.push_back(GammaDirection.Angle(recoilDir));
+ gammaEDC.push_back(Egamma*((1-recoilBeta*TMath::Cos(ThetaGR.back()))/TMath::Sqrt(1-recoilBeta*recoilBeta)));
+ ThetaG.push_back(GammaDirection.Theta()/deg);
+ PhiG.push_back(GammaDirection.Phi()/deg);
+
+ // EDC Optimization
+ //if(T_MUGAST_FPMW > 13e3 && T_MUGAST_FPMW < 26e3 && T_FPMW_HF > 12.1e3 && T_FPMW_HF < 13.2e3){
+/* if(T_MUGAST_FPMW > 7e3 && T_MUGAST_FPMW < 20e3 && T_FPMW_HF > 13.7e3 && T_FPMW_HF < 14.8e3){
+ TVector3 gammaDir = GammaDirection;
+ double tmpGR = 0;
+ for(UInt_t j = 0 ; j < 360 ; j++){
+ gammaDir.RotateZ(1*TMath::DegToRad());
+ tmpGR = gammaDir.Angle(recoilDir);
+ mEDC_PhiAgata->Fill(j,
+ Egamma*1000*((1-recoilBeta*TMath::Cos(tmpGR))/TMath::Sqrt(1-recoilBeta*recoilBeta)));
+ }
+ // GammaDirection is a unitary vector so offset does not work;
+ gammaDir = GammaHit;
+ TVector3 zOffset ;
+ TVector3 tmpGamma;
+ for(Int_t j = -10 ; j <=10 ; j++){
+ zOffset.SetXYZ(0,0,j);
+ tmpGamma = gammaDir+zOffset;
+ //std::cout << "Z = " << tmpGamma.Z() << " = " << gammaDir.Z() << "+ " << zOffset.Z() << std::endl;
+ tmpGR = tmpGamma.Angle(recoilDir);
+ mEDC_OffsetAgata->Fill(j,
+ Egamma*1000*((1-recoilBeta*TMath::Cos(tmpGR))/TMath::Sqrt(1-recoilBeta*recoilBeta)));
+ }
+
+ }
+*/
+
+ //Test (does not work): Energy loss in Au for 20O
+ E4_tmp = reaction.GetEnergyImpulsionLab_4().E()-reaction.GetNucleus4()->Mass();
+ Energy4 = E4_tmp;
+ //if(T_MUGAST_FPMW > 13e3 && T_MUGAST_FPMW < 26e3 && T_FPMW_HF > 12.1e3 && T_FPMW_HF < 13.2e3){std::cout << E4_tmp;}
+ //cout << E4_tmp << "\t";
+ E4_tmp = HeavyTarget.Slow(E4_tmp, TargetThickness*0.5,Theta4*deg);
+ //if(T_MUGAST_FPMW > 13e3 && T_MUGAST_FPMW < 26e3 && T_FPMW_HF > 12.1e3 && T_FPMW_HF < 13.2e3){std::cout << "\t" << E4_tmp;}
+ //cout << E4_tmp << "\t";
+ //cout << Theta4 << "\t";
+ E4_tmp = HeavyAu.Slow(E4_tmp, AuThickness, Theta4*deg);
+
+ //cout <<Theta4 << "\t\t" << AuThickness << "\t" << AuThickness/TMath::Cos(Theta4*deg) << endl;
+ //if(T_MUGAST_FPMW > 13e3 && T_MUGAST_FPMW < 26e3 && T_FPMW_HF > 12.1e3 && T_FPMW_HF < 13.2e3 && ExCalib.back()>1.5&& ExCalib.back()<5.5)
+ //cout << reaction.GetEnergyImpulsionLab_4().E()-reaction.GetNucleus4()->Mass() << "\t" << E4_tmp << "\n";
+
+ TLorentzVector LVector4_tmp_1=reaction.GetEnergyImpulsionLab_4();
+ TVector3 imp4_1=LVector4_tmp_1.Vect();
+ imp4_1.SetMag(sqrt((E4_tmp+reaction.GetNucleus4()->Mass())*(E4_tmp+reaction.GetNucleus4()->Mass())-reaction.GetNucleus4()->Mass()*reaction.GetNucleus4()->Mass()));
+ TLorentzVector LVector4_1(imp4_1,E4_tmp+reaction.GetNucleus4()->Mass());
+ //double bb;
+ //bb = sqrt(1-(reaction.GetNucleus4()->Mass()*reaction.GetNucleus4()->Mass()/(E4_tmp+reaction.GetNucleus4()->Mass())/(E4_tmp+reaction.GetNucleus4()->Mass())));
+ //std::cout << "Beta1 " << LVector4_1.Beta() << "\tBeta2 " << bb << std::endl;
+ //if(T_MUGAST_FPMW > 13e3 && T_MUGAST_FPMW < 26e3 && T_FPMW_HF > 12.1e3 && T_FPMW_HF < 13.2e3){std::cout << "\t" << E4_tmp << std::endl;}
+ //cout<<"E4_tmp="<<E4_tmp<<" E4 reaction="<<reaction.GetEnergyImpulsionLab_4().E()-reaction.GetNucleus4()->Mass()<<" E4Lorentz="<< LVector4.E()-reaction.GetNucleus4()->Mass();
+ //cout<<LVector4.Beta()<<endl;
+
+ double beta_fix = 0.126;
+ EDC_MUGAST.push_back(Egamma*((1-beta_fix*TMath::Cos(ThetaG.back()*deg))/TMath::Sqrt(1-beta_fix*beta_fix)));
+
+ // Boost back in CM
+ /* TVector3 betaM2=reaction.GetEnergyImpulsionLab_4().BoostVector();
+ TVector3 betaM=LVector4.BoostVector();
+ betaM.SetPhi(HitDirection.Phi()+3.14159);
+ betaM.SetMag(0.125);
+ //cout<<"betaM="<<betaM.Mag()<<" beta original="<<betaM2.Mag()<<endl;
+ beta_MUGAST.push_back(betaM.Mag());
+ betaM=TVector3(0,0,0.125);
+ GammaLVM.Boost(-betaM);
+ ThetaGamma.push_back(betaM.Angle(GammaLVM.BoostVector())/deg);
+ // Get EDC
+ EDC_MUGAST.push_back(GammaLVM.Energy());
+ // Construct LV in lab frame
+ double ThetaGammaV = GammaDirection.Angle(BeamDirection)/deg;
+ TLorentzVector GammaLVV;
+ GammaLVV.SetPx(Egamma*GammaDirection.X());
+ GammaLVV.SetPy(Egamma*GammaDirection.Y());
+ GammaLVV.SetPz(Egamma*GammaDirection.Z());
+ GammaLVV.SetE(Egamma);
+ // Boost back in CM
+ TVector3 betaV(0,0,-mBeta); */
+ beta_VAMOS.push_back(LVector4_1.Beta());
+ //std::cout << beta_VAMOS.back() << "\t\t" << recoilBeta-beta_VAMOS.back() << std::endl;
+ //GammaLVV.Boost(betaV);
+ //EDC_VAMOS.push_back(GammaLVV.Energy());
+
+ //IRENE's TEST WITH ENERGY LOSS IN GOLD
+
+ vamosEDC.push_back(Egamma*((1-beta_VAMOS.back()*TMath::Cos(ThetaGR.back()))/TMath::Sqrt(1-beta_VAMOS.back()*beta_VAMOS.back())));
+
+ //Graph for checking degrader thickness
+ /*
+ if(T_MUGAST_FPMW > 13e3 && T_MUGAST_FPMW < 26e3 && T_FPMW_HF > 12.1e3 && T_FPMW_HF < 13.2e3 && Ex.back()>1.0 && Ex.back()<5.5){
+ double tmpGR = 0;
+ TVector3 gammaDir = GammaDirection;
+
+ E4_tmp = reaction.GetEnergyImpulsionLab_4().E()-reaction.GetNucleus4()->Mass();
+ E4_tmp = HeavyTarget.Slow(E4_tmp, TargetThickness*0.5,Theta4*deg);
+
+ for(Int_t j = 0 ; j < 30 ; j++){
+ double off = 0.;
+ //off = (j/10.)*micrometer;
+ off = j/10.*micrometer;
+ double E4_fix = HeavyAu.Slow(E4_tmp, AuThickness+off, Theta4*deg);
+ TLorentzVector LVector4_tmp=reaction.GetEnergyImpulsionLab_4();
+ TVector3 imp4=LVector4_tmp.Vect();
+ imp4.SetMag(sqrt((E4_fix+reaction.GetNucleus4()->Mass())*(E4_fix+reaction.GetNucleus4()->Mass())-reaction.GetNucleus4()->Mass()*reaction.GetNucleus4()->Mass()));
+ TLorentzVector LVector4(imp4,E4_fix+reaction.GetNucleus4()->Mass());
+ double beta4 = LVector4.Beta();
+
+ tmpGR = gammaDir.Angle(recoilDir);
+ backingThick->Fill((AuThickness+off)/micrometer+0.025,Egamma*1000*((1-beta4*TMath::Cos(tmpGR))/TMath::Sqrt(1-beta4*beta4)));
+ //printf("%13.5f %13.5f\n",(AuThickness+off)/micrometer,Egamma*1000*((1-beta4*TMath::Cos(tmpGR))/TMath::Sqrt(1-beta4*beta4)));
+ }
+ }*/ //fine grafico
+ }
+
+ //////////////////// TRACKING //////////////////////
+
+ for(int j = 0; j < nbTrack; j++){ // all multiplicity
+
+ // Measured E
+ double Etrack=trackE[j]/1000.;
+ TVector3 TrackHit(trackX1[j],trackY1[j],trackZ1[j]+agata_zShift-0.578);
+ TrackHit.RotateZ(272.235*TMath::DegToRad());
+ TVector3 GammaDirection = TrackHit-BeamImpact;
+ GammaDirection = GammaDirection.Unit();
+ if(j==0) radius = TrackHit.Mag();
+
+ // Construct LV in lab frame
+ TLorentzVector TrackLVM;
+ TrackLVM.SetPx(Etrack*GammaDirection.X());
+ TrackLVM.SetPy(Etrack*GammaDirection.Y());
+ TrackLVM.SetPz(Etrack*GammaDirection.Z());
+ TrackLVM.SetE(Etrack);
+
+ // Test EDC calculated taking into account the recoil direction (Irene and Alain)
+ //////recoilBeta = reaction.GetNucleus4()->GetEnergyImpulsion().Beta();
+ //////recoilDir = reaction.GetNucleus4()->GetEnergyImpulsion().Vect();
+ recoilBeta = reaction.GetEnergyImpulsionLab_4().Beta();
+ recoilDir = reaction.GetEnergyImpulsionLab_4().Vect();
+ ThetaGR.push_back(GammaDirection.Angle(recoilDir));
+
+ trackEDC.push_back(Etrack*((1-recoilBeta*TMath::Cos(ThetaGR.back()))/TMath::Sqrt(1-recoilBeta*recoilBeta)));
+
+/* E4_tmp = reaction.GetEnergyImpulsionLab_4().E()-reaction.GetNucleus4()->Mass();
+ //Energy4 = E4_tmp;
+ E4_tmp = HeavyTarget.Slow(E4_tmp, TargetThickness*0.5,Theta4*deg);
+ E4_tmp = HeavyAu.Slow(E4_tmp, AuThickness, Theta4*deg);
+ TLorentzVector LVector4_tmp_1=reaction.GetEnergyImpulsionLab_4();
+ TVector3 imp4_1=LVector4_tmp_1.Vect();
+ imp4_1.SetMag(sqrt((E4_tmp+reaction.GetNucleus4()->Mass())*(E4_tmp+reaction.GetNucleus4()->Mass())-reaction.GetNucleus4()->Mass()*reaction.GetNucleus4()->Mass()));
+ TLorentzVector LVector4_1(imp4_1,E4_tmp+reaction.GetNucleus4()->Mass()); */
+ goldEDC.push_back(Etrack*((1-beta_VAMOS.back()*TMath::Cos(ThetaGR.back()))/TMath::Sqrt(1-beta_VAMOS.back()*beta_VAMOS.back())));
+
+
+ /* double tmpGR = 0;
+ TVector3 gammaDir = GammaDirection;
+ for(UInt_t j = 0 ; j < 360 ; j++){
+ gammaDir.RotateZ(1*TMath::DegToRad());
+ tmpGR = gammaDir.Angle(recoilDir);
+ mEDC_PhiAgata->Fill(j,Etrack*1000*((1-recoilBeta*TMath::Cos(tmpGR))/TMath::Sqrt(1-recoilBeta*recoilBeta)));
+ }
+ // GammaDirection is a unitary vector so offset does not work;
+ gammaDir = TrackHit;
+ TVector3 zOffset ;
+ TVector3 tmpGamma;
+ for(Int_t j = -10 ; j <=10 ; j++){
+ zOffset.SetXYZ(0,0,j);
+ tmpGamma = gammaDir+zOffset;
+ //std::cout << "Z = " << tmpGamma.Z() << " = " << gammaDir.Z() << "+ " << zOffset.Z() << std::endl;
+ tmpGR = tmpGamma.Angle(recoilDir);
+ mEDC_OffsetAgata->Fill(j,
+ Etrack*1000*((1-recoilBeta*TMath::Cos(tmpGR))/TMath::Sqrt(1-recoilBeta*recoilBeta)));
+ }*/
+
+ }
+
+ // radius = TMath::Sqrt(trackX1[0]*trackX1[0]+trackY1[0]*trackY1[0]+(trackZ1[0]+*trackZ1[0]);
+
+
+
+ // Agata add back is not always multiplicity 1 ?? NO, not necessarily!
+ if(nbAdd==1){
+ TLorentzVector GammaLV;
+ // Measured E
+ double Egamma=AddE[0]/1000.; // From keV to MeV
+ // Gamma detection position
+ // TrackZ1 to be corrected there is a shift of +51mm
+ TVector3 GammaHit(AddX[0],AddY[0],AddZ[0]+agata_zShift);
+ // TVector3 GammaHit(trackX1[0],trackY1[0],trackZ1[0]);
+ TVector3 GammaDirection = GammaHit-BeamImpact;
+ GammaDirection = GammaDirection.Unit();
+ // Beta from Two body kinematic
+ //TVector3 beta = reaction.GetEnergyImpulsionLab_4().BoostVector();
+ // Beta from the Beam mid target
+ reaction.GetKinematicLine4();
+ //TVector3 beta(0,0,-reaction.GetNucleus4()->GetBeta());
+ TVector3 beta(0,0,-mBeta);
+
+ double ThetaGamma = GammaDirection.Angle(BeamDirection)/deg;
+
+ // Construct LV in lab frame
+ GammaLV.SetPx(Egamma*GammaDirection.X());
+ GammaLV.SetPy(Egamma*GammaDirection.Y());
+ GammaLV.SetPz(Egamma*GammaDirection.Z());
+ GammaLV.SetE(Egamma);
+ // Boost back in CM
+ GammaLV.Boost(beta);
+ // Get EDC
+ AddBack_EDC.push_back(GammaLV.Energy());
+ }
+
+
+
+
+
+
+
+
+// // Agata by track
+// /*
+// for(int j=0; j<nbTrack; j++){ // all multiplicity
+// TLorentzVector GammaLV;
+// // Measured E
+// double Egamma=trackE[j]/1000.; // From keV to MeV
+// // Gamma detection position
+// // TrackZ1 to be corrected there is a shift of +51mm
+// TVector3 GammaHit(trackX1[j],trackY1[j],trackZ1[j]+agata_zShift);
+// // TVector3 GammaHit(trackX1[0],trackY1[0],trackZ1[0]);
+// // Gamma Direction
+// TVector3 GammaDirection = GammaHit-BeamImpact;
+// GammaDirection = GammaDirection.Unit();
+// // Beta from Two body kinematic
+// //TVector3 beta = reaction.GetEnergyImpulsionLab_4().BoostVector();
+// // Beta from the Beam mid target
+// reaction.GetKinematicLine4();
+// TVector3 beta(0,0,-reaction.GetNucleus4()->GetBeta());
+//
+// double ThetaGamma = GammaDirection.Angle(BeamDirection)/deg;
+// // Construct LV in lab frame
+// GammaLV.SetPx(Egamma*GammaDirection.X());
+// GammaLV.SetPy(Egamma*GammaDirection.Y());
+// GammaLV.SetPz(Egamma*GammaDirection.Z());
+// GammaLV.SetE(Egamma);
+// // Boost back in CM
+// GammaLV.Boost(beta);
+// // Get EDC
+// 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){
+// TLorentzVector GammaLV;
+// // Measured E
+// double Egamma = AddE[i]/1000.; // From keV to MeV
+// // Gamma detection position
+// // TrackZ1 to be corrected there is a shift of +51mm
+// TVector3 GammaHit(AddX[i],AddY[i],AddZ[i]+agata_zShift);
+// // TVector3 GammaHit(trackX1[0],trackY1[0],trackZ1[0]);
+//
+//
+//// cout << "gDir " << GammaHit.X() << " "
+//// << GammaHit.Y() << " "
+//// << GammaHit.Z() << " ";
+// /* TEST adding constant to Phi */
+// //double phitemp = GammaHit.Phi();
+// //double phiconst = -2.0;
+// //GammaHit.SetPhi(phitemp+phiconst);
+//// GammaHit.RotateZ(0.0);
+//
+//// cout << "gDir " << GammaHit.X() << " "
+//// << GammaHit.Y() << " "
+//// << GammaHit.Z() << endl;
+//
+// // Gamma Direction
+// TVector3 GammaDirection = GammaHit-BeamImpact;
+// GammaDirection = GammaDirection.Unit();
+//
+//
+// /* Beta from Two body kinematic */
+// TVector3 beta = reaction.GetEnergyImpulsionLab_4().BoostVector();
+// //TVector3 beta = reaction->GetEnergyImpulsionLab_4().BoostVector();
+//
+//// cout << "bDir " << beta.X() << " "
+//// << beta.Y() << " "
+//// << beta.Z() ;
+//// beta.RotateX(0.002847); beta.RotateY(3.144869); beta.RotateZ(0.095923);
+// beta.RotateY(M_PI);
+//// cout << "bDir " << beta.X() << " "
+//// << beta.Y() << " "
+//// << beta.Z() << endl;
+//
+// /* Original beta */
+//// TVector3 beta(0,0,-0.1257);
+//
+// // For beta rotation minimization
+// AGATA_Theta.push_back(GammaDirection.Theta());
+// AGATA_Phi.push_back(GammaDirection.Phi());
+// AGATA_GammaPx.push_back(Egamma*GammaDirection.X());
+// AGATA_GammaPy.push_back(Egamma*GammaDirection.Y());
+// AGATA_GammaPz.push_back(Egamma*GammaDirection.Z());
+// AGATA_GammaE.push_back(Egamma);
+// AGATA_OrigBetaX.push_back(beta.X());
+// AGATA_OrigBetaY.push_back(beta.Y());
+// AGATA_OrigBetaZ.push_back(beta.Z());
+// AGATA_OrigBetaMag.push_back(beta.Mag());
+//
+// /* Other fills */
+// double ThetaGamma = GammaDirection.Angle(BeamDirection)/deg;
+// // Construct LV in lab frame
+// GammaLV.SetPx(Egamma*GammaDirection.X());
+// GammaLV.SetPy(Egamma*GammaDirection.Y());
+// GammaLV.SetPz(Egamma*GammaDirection.Z());
+// GammaLV.SetE(Egamma);
+// // Boost back in CM
+// GammaLV.Boost(beta);
+// // Get EDC
+// AddBack_EDC.push_back(GammaLV.Energy());
+//
+// if(i==0){
+// //First event in loop
+// AddBack_EDC_Event1.push_back(GammaLV.Energy());
+// } else {
+// //second, third, fourth...
+// AddBack_EDC_Event2.push_back(GammaLV.Energy());
+// }
+//
+// }
+//
+// float offset[21] = {0.,27.4,46.7,41.1,71.2,58.7,32.1,40.4,-46.3,15.1,9.9,-39.1,-39.8,-5.3,18.5,0.7,-28.5,-7,-22,-8.5,-8.5};
+// float slope[21] = {1.,0.96,0.93,0.93,0.89,0.91,0.95,0.93,1.06,0.97,0.99,1.05,1.04,1.01,0.97,1.,1.03,0.99,1.,1.02,1.02};
+//
+// for(unsigned int i=0; i<MW_Nr; i++){
+// MWT[i] = -1500;
+// for(unsigned short j=0; j<21; j++)
+// if(MW_N[i]==j)
+// MWT[i] = (MW_T[i]-offset[j])/slope[j];
+// }
+
+
+//cout << "end_TreatEvent" << endl;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void FillSolidAngles(TH1F* hSA, TH1F* hDet, TH1F* hEmm){
+ if(!filledCline){
+ cout << RED << "FILLING CLINE VECTOR! SHOULD ONLY OCCUR ONCE!" << RESET << endl;
+ for(int t=0; t<NumThetaAngleBins; t++){
+ double angleMin = (t)*(180./NumThetaAngleBins);
+ double angleMax = (t+1)*(180./NumThetaAngleBins);
+ cout << " Angle " << angleMin
+ << " to " << angleMax
+ << endl;
+ clineVal.push_back(2.0*M_PI*(cos(angleMin*degtorad) - cos(angleMax*degtorad)));
+ clineX.push_back(angleMin+((angleMax-angleMin)/2.0));
+ filledCline=true;
+ }
+ }
+
+ for (int i = 0; i < hSA->GetNbinsX(); ++i){
+ double val = hDet->GetBinContent(i) / hEmm->GetBinContent(i);
+ double valerr = val * sqrt(
+ pow(hDet->GetBinError(i) / hDet->GetBinContent(i), 2) +
+ pow(hEmm->GetBinError(i) / hEmm->GetBinContent(i), 2) );
+ if (isnan(val)) { val = 0; valerr = 0; }
+ val *= clineVal.at(i);
+ valerr *= clineVal.at(i);
+
+ hSA->SetBinContent(i, val);
+ hSA->SetBinError(i, valerr);
+ }
+}
+
+////////////////////////////////////////////////////////////////////////////////
+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 ;
+
+// TFile* ppp2 = new TFile("TargetT.root", "Recreate");
+// backingThick->Write();
+// ppp2->Close();
+//
+ if(isSim && !isPhaseSpace){
+
+ //TObjArray HistList(0);
+ TList *HistList = new TList();
+
+// //// MUST 2 DETECTOR #5
+// //auto Efficiency_CM_MM = new TH1F(*ThetaCM_detected_MM);
+// //Efficiency_CM_MM->SetName("Efficiency_CM_MM");
+// //Efficiency_CM_MM->SetTitle("Efficiency_CM_MM");
+// //Efficiency_CM_MM->Sumw2();
+// //auto Efficiency_Lab_MM = new TH1F(*ThetaLab_detected_MM);
+// //Efficiency_Lab_MM->SetName("Efficiency_Lab_MM");
+// //Efficiency_Lab_MM->SetTitle("Efficiency_Lab_MM");
+// //Efficiency_Lab_MM->Sumw2();
+// //
+// //auto SolidAngle_CM_MM = new TH1F(*ThetaCM_detected_MM);
+// //SolidAngle_CM_MM->SetName("SolidAngle_CM_MM");
+// //SolidAngle_CM_MM->SetTitle("SolidAngle_CM_MM");
+// //auto SolidAngle_Lab_MM = new TH1F(*ThetaLab_detected_MM);
+// //SolidAngle_Lab_MM->SetName("SolidAngle_Lab_MM");
+// //SolidAngle_Lab_MM->SetTitle("SolidAngle_Lab_MM");
+//
+// //auto Efficiency_CM_MM_wV = new TH1F(*ThetaCM_detected_MM_wV);
+// //Efficiency_CM_MM_wV->SetName("Efficiency_CM_MM_wV");
+// //Efficiency_CM_MM_wV->SetTitle("Efficiency_CM_MM_wV");
+// //Efficiency_CM_MM_wV->Sumw2();
+// //auto Efficiency_Lab_MM_wV = new TH1F(*ThetaLab_detected_MM_wV);
+// //Efficiency_Lab_MM_wV->SetName("Efficiency_Lab_MM_wV");
+// //Efficiency_Lab_MM_wV->SetTitle("Efficiency_Lab_MM_wV");
+// //Efficiency_Lab_MM_wV->Sumw2();
+// //
+// //auto SolidAngle_CM_MM_wV = new TH1F(*ThetaCM_detected_MM_wV);
+// //SolidAngle_CM_MM_wV->SetName("SolidAngle_CM_MM_wV");
+// //SolidAngle_CM_MM_wV->SetTitle("SolidAngle_CM_MM_wV");
+// //auto SolidAngle_Lab_MM_wV = new TH1F(*ThetaLab_detected_MM_wV);
+// //SolidAngle_Lab_MM_wV->SetName("SolidAngle_Lab_MM_wV");
+// //SolidAngle_Lab_MM_wV->SetTitle("SolidAngle_Lab_MM_wV");
+//
+// //auto Efficiency_CM_MM_Unb = new TH1F(*ThetaCM_detected_MM_Unb);
+// //Efficiency_CM_MM_Unb->SetName("Efficiency_CM_MM_Unb");
+// //Efficiency_CM_MM_Unb->SetTitle("Efficiency_CM_MM_Unb");
+// //Efficiency_CM_MM_Unb->Sumw2();
+// //auto Efficiency_Lab_MM_Unb = new TH1F(*ThetaLab_detected_MM_Unb);
+// //Efficiency_Lab_MM_Unb->SetName("Efficiency_Lab_MM_Unb");
+// //Efficiency_Lab_MM_Unb->SetTitle("Efficiency_Lab_MM_Unb");
+// //Efficiency_Lab_MM_Unb->Sumw2();
+//
+// //auto SolidAngle_CM_MM_Unb = new TH1F(*ThetaCM_detected_MM_Unb);
+// //SolidAngle_CM_MM_Unb->SetName("SolidAngle_CM_MM_Unb");
+// //SolidAngle_CM_MM_Unb->SetTitle("SolidAngle_CM_MM_Unb");
+// //auto SolidAngle_Lab_MM_Unb = new TH1F(*ThetaLab_detected_MM_Unb);
+// //SolidAngle_Lab_MM_Unb->SetName("SolidAngle_Lab_MM_Unb");
+// //SolidAngle_Lab_MM_Unb->SetTitle("SolidAngle_Lab_MM_Unb");
+//
+// //Efficiency_CM_MM->Divide(ThetaCM_emmitted);
+// //Efficiency_Lab_MM->Divide(ThetaLab_emmitted);
+//
+// //Efficiency_CM_MM_wV->Divide(ThetaCM_emmitted);
+// //Efficiency_Lab_MM_wV->Divide(ThetaLab_emmitted);
+// //
+// //Efficiency_CM_MM_Unb->Divide(ThetaCM_emmitted);
+// //Efficiency_Lab_MM_Unb->Divide(ThetaLab_emmitted);
+// //
+// ////double dt_MM = 180./Efficiency_Lab_MM->GetNbinsX();
+// ////cout << "Angular infinitesimal (MM) = " << dt_MM << "deg " << endl;
+// ////auto Cline_MM = new TF1("Cline_MM",Form("1./(2*%f*sin(x*%f/180.)*%f*%f/180.)",M_PI,M_PI,dt_MM,M_PI),0,180);
+//
+// //FillSolidAngles(SolidAngle_CM_MM, ThetaCM_detected_MM, ThetaCM_emmitted);
+// //FillSolidAngles(SolidAngle_Lab_MM, ThetaLab_detected_MM, ThetaLab_emmitted);
+// //
+// //FillSolidAngles(SolidAngle_CM_MM_wV, ThetaCM_detected_MM_wV, ThetaCM_emmitted);
+// //FillSolidAngles(SolidAngle_Lab_MM_wV, ThetaLab_detected_MM_wV, ThetaLab_emmitted);
+// //
+// //FillSolidAngles(SolidAngle_CM_MM_Unb, ThetaCM_detected_MM_Unb, ThetaCM_emmitted);
+// //FillSolidAngles(SolidAngle_Lab_MM_Unb, ThetaLab_detected_MM_Unb, ThetaLab_emmitted);
+// //
+// ////SolidAngle_CM_MM->Divide(ThetaCM_emmitted);
+// ////SolidAngle_CM_MM->Divide(Cline_MM,1);
+// ////SolidAngle_CM_MM->Divide(Cline);
+// ////SolidAngle_Lab_MM->Divide(ThetaLab_emmitted);
+// ////SolidAngle_Lab_MM->Divide(Cline_MM,1);
+// ////SolidAngle_Lab_MM->Divide(Cline);
+//
+// //HistList->Add(ThetaCM_emmitted);
+// //HistList->Add(ThetaLab_emmitted);
+// //HistList->Add(ThetaCM_detected_MM);
+// //HistList->Add(ThetaLab_detected_MM);
+// //HistList->Add(Efficiency_CM_MM);
+// //HistList->Add(Efficiency_Lab_MM);
+// //HistList->Add(SolidAngle_CM_MM);
+// //HistList->Add(SolidAngle_Lab_MM);
+// //HistList->Add(ThetaCM_detected_MM_wV);
+// //HistList->Add(ThetaLab_detected_MM_wV);
+// //HistList->Add(Efficiency_CM_MM_wV);
+// //HistList->Add(Efficiency_Lab_MM_wV);
+// //HistList->Add(SolidAngle_CM_MM_wV);
+// //HistList->Add(SolidAngle_Lab_MM_wV);
+// //HistList->Add(ThetaCM_detected_MM_Unb);
+// //HistList->Add(ThetaLab_detected_MM_Unb);
+// //HistList->Add(Efficiency_CM_MM_Unb);
+// //HistList->Add(Efficiency_Lab_MM_Unb);
+// //HistList->Add(SolidAngle_CM_MM_Unb);
+// //HistList->Add(SolidAngle_Lab_MM_Unb);
+//
+// ////HistList->Add(Cline_MM);
+// ////HistList->Add(Cline);
+//
+//
+ // MUGAST
+ auto Efficiency_CM_MG = new TH1F(*ThetaCM_detected_MG);
+ Efficiency_CM_MG->SetName("Efficiency_CM_MG");
+ Efficiency_CM_MG->SetTitle("Efficiency_CM_MG");
+ Efficiency_CM_MG->Sumw2();
+ auto Efficiency_Lab_MG = new TH1F(*ThetaLab_detected_MG);
+ Efficiency_Lab_MG->SetName("Efficiency_Lab_MG");
+ Efficiency_Lab_MG->SetTitle("Efficiency_Lab_MG");
+ Efficiency_Lab_MG->Sumw2();
+
+ auto SolidAngle_CM_MG = new TH1F(*ThetaCM_detected_MG);
+ SolidAngle_CM_MG->SetName("SolidAngle_CM_MG");
+ SolidAngle_CM_MG->SetTitle("SolidAngle_CM_MG");
+ auto SolidAngle_Lab_MG = new TH1F(*ThetaLab_detected_MG);
+ SolidAngle_Lab_MG->SetName("SolidAngle_Lab_MG");
+ SolidAngle_Lab_MG->SetTitle("SolidAngle_Lab_MG");
+
+ auto Efficiency_CM_MG_wV = new TH1F(*ThetaCM_detected_MG_wV);
+ Efficiency_CM_MG_wV->SetName("Efficiency_CM_MG_wV");
+ Efficiency_CM_MG_wV->SetTitle("Efficiency_CM_MG_wV");
+ Efficiency_CM_MG_wV->Sumw2();
+ auto Efficiency_Lab_MG_wV = new TH1F(*ThetaLab_detected_MG_wV);
+ Efficiency_Lab_MG_wV->SetName("Efficiency_Lab_MG_wV");
+ Efficiency_Lab_MG_wV->SetTitle("Efficiency_Lab_MG_wV");
+ Efficiency_Lab_MG_wV->Sumw2();
+
+ auto SolidAngle_CM_MG_wV = new TH1F(*ThetaCM_detected_MG_wV);
+ SolidAngle_CM_MG_wV->SetName("SolidAngle_CM_MG_wV");
+ SolidAngle_CM_MG_wV->SetTitle("SolidAngle_CM_MG_wV");
+ auto SolidAngle_Lab_MG_wV = new TH1F(*ThetaLab_detected_MG_wV);
+ SolidAngle_Lab_MG_wV->SetName("SolidAngle_Lab_MG_wV");
+ SolidAngle_Lab_MG_wV->SetTitle("SolidAngle_Lab_MG_wV");
+
+ auto Efficiency_CM_MG_Unb = new TH1F(*ThetaCM_detected_MG_Unb);
+ Efficiency_CM_MG_Unb->SetName("Efficiency_CM_MG_Unb");
+ Efficiency_CM_MG_Unb->SetTitle("Efficiency_CM_MG_Unb");
+ Efficiency_CM_MG_Unb->Sumw2();
+ auto Efficiency_Lab_MG_Unb = new TH1F(*ThetaLab_detected_MG_Unb);
+ Efficiency_Lab_MG_Unb->SetName("Efficiency_Lab_MG_Unb");
+ Efficiency_Lab_MG_Unb->SetTitle("Efficiency_Lab_MG_Unb");
+ Efficiency_Lab_MG_Unb->Sumw2();
+
+ auto SolidAngle_CM_MG_Unb = new TH1F(*ThetaCM_detected_MG_Unb);
+ SolidAngle_CM_MG_Unb->SetName("SolidAngle_CM_MG_Unb");
+ SolidAngle_CM_MG_Unb->SetTitle("SolidAngle_CM_MG_Unb");
+ auto SolidAngle_Lab_MG_Unb = new TH1F(*ThetaLab_detected_MG_Unb);
+ SolidAngle_Lab_MG_Unb->SetName("SolidAngle_Lab_MG_Unb");
+ SolidAngle_Lab_MG_Unb->SetTitle("SolidAngle_Lab_MG_Unb");
+
+ Efficiency_CM_MG->Divide(ThetaCM_emmitted);
+ Efficiency_Lab_MG->Divide(ThetaLab_emmitted);
+
+ Efficiency_CM_MG_wV->Divide(ThetaCM_emmitted);
+ Efficiency_Lab_MG_wV->Divide(ThetaLab_emmitted);
+
+ Efficiency_CM_MG_Unb->Divide(ThetaCM_emmitted);
+ Efficiency_Lab_MG_Unb->Divide(ThetaLab_emmitted);
+
+ //double dt_MG = 180./Efficiency_Lab_MG->GetNbinsX();
+ //cout << "Angular infinitesimal (MG) = " << dt_MG << "deg " << endl;
+ //auto Cline_MG = new TF1("Cline_MG",Form("1./(2*%f*sin(x*%f/180.)*%f*%f/180.)",M_PI,M_PI,dt_MG,M_PI),0,180);
+
+ /* Testing method for better errors in SolidAngle histograms */
+ FillSolidAngles(SolidAngle_CM_MG, ThetaCM_detected_MG, ThetaCM_emmitted);
+ FillSolidAngles(SolidAngle_CM_MG_wV, ThetaCM_detected_MG_wV, ThetaCM_emmitted);
+ FillSolidAngles(SolidAngle_CM_MG_Unb, ThetaCM_detected_MG_Unb, ThetaCM_emmitted);
+ //SolidAngle_CM_MG->Divide(Cline_MG,1);
+ //SolidAngle_CM_MG->Divide(Cline);
+ FillSolidAngles(SolidAngle_Lab_MG, ThetaLab_detected_MG, ThetaLab_emmitted);
+ FillSolidAngles(SolidAngle_Lab_MG_wV, ThetaLab_detected_MG_wV, ThetaLab_emmitted);
+ FillSolidAngles(SolidAngle_Lab_MG_Unb, ThetaLab_detected_MG_Unb, ThetaLab_emmitted);
+ //SolidAngle_Lab_MG->Divide(Cline_MG,1);
+ //SolidAngle_Lab_MG->Divide(Cline);
+
+ HistList->Add(ThetaCM_detected_MG);
+ HistList->Add(ThetaLab_detected_MG);
+ HistList->Add(Efficiency_CM_MG);
+ HistList->Add(Efficiency_Lab_MG);
+ HistList->Add(SolidAngle_CM_MG);
+ HistList->Add(SolidAngle_Lab_MG);
+ HistList->Add(ThetaCM_detected_MG_wV);
+ HistList->Add(ThetaLab_detected_MG_wV);
+ HistList->Add(Efficiency_CM_MG_wV);
+ HistList->Add(Efficiency_Lab_MG_wV);
+ HistList->Add(SolidAngle_CM_MG_wV);
+ HistList->Add(SolidAngle_Lab_MG_wV);
+ HistList->Add(ThetaCM_detected_MG_Unb);
+ HistList->Add(ThetaLab_detected_MG_Unb);
+ HistList->Add(Efficiency_CM_MG_Unb);
+ HistList->Add(Efficiency_Lab_MG_Unb);
+ HistList->Add(SolidAngle_CM_MG_Unb);
+ HistList->Add(SolidAngle_Lab_MG_Unb);
+
+ //HistList->Add(Cline_MG);
+
+ // MUGAST INDIVIDUAL
+ TH1F *SolidAngle_CM_MGX[7];
+ string base1 = "SolidAngle_CM_MG";
+ for(int i=0;i<7;i++){
+ int j=i+1;
+ string name = base1 + to_string(j);
+ SolidAngle_CM_MGX[i] = new TH1F(*ThetaCM_detected_MGX[i]);
+ SolidAngle_CM_MGX[i]->SetName(name.c_str());
+ SolidAngle_CM_MGX[i]->SetTitle(name.c_str());
+ FillSolidAngles(SolidAngle_CM_MGX[i], ThetaCM_detected_MGX[i], ThetaCM_emmitted);
+ //SolidAngle_CM_MGX[i]->Divide(Cline_MG,1);
+ //SolidAngle_CM_MGX[i]->Divide(Cline);
+ }
+
+ TH1F *SolidAngle_Lab_MGX[7];
+ string base2 = "SolidAngle_Lab_MG";
+ for(int i=0;i<7;i++){
+ int j=i+1;
+ string name = base2 + to_string(j);
+ SolidAngle_Lab_MGX[i] = new TH1F(*ThetaLab_detected_MGX[i]);
+ SolidAngle_Lab_MGX[i]->SetName(name.c_str());
+ SolidAngle_Lab_MGX[i]->SetTitle(name.c_str());
+ FillSolidAngles(SolidAngle_Lab_MGX[i], ThetaLab_detected_MGX[i], ThetaLab_emmitted);
+ //SolidAngle_Lab_MGX[i]->Divide(Cline_MG,1);
+ //SolidAngle_Lab_MGX[i]->Divide(Cline);
+ }
+
+ for(int i=0; i<7; i++){HistList->Add(ThetaCM_detected_MGX[i]);}
+ for(int i=0; i<7; i++){HistList->Add(ThetaLab_detected_MGX[i]);}
+
+ for(int i=0; i<7; i++){HistList->Add(SolidAngle_CM_MGX[i]);}
+ for(int i=0; i<7; i++){HistList->Add(SolidAngle_Lab_MGX[i]);}
+
+ auto clineValGraph = new TGraph(NumThetaAngleBins);
+ clineValGraph->SetName("clineValGraph");
+ clineValGraph->SetTitle("clineValGraph");
+ for(int b=0; b<NumThetaAngleBins; b++){
+ clineValGraph->SetPoint(b,clineX.at(b),clineVal.at(b));
+ }
+ HistList->Add(clineValGraph);
+
+ auto HistoFile = new TFile("~/Programs/nptoolv3/Projects/e775s/analysis/SolidAngleHistFiles/SAHF_New.root","RECREATE");
+ HistList->Write();
+ HistoFile->Close();
+
+
+ cout << "!!! MAKE SURE YOU RENAME THE SOLID ANGLE FILE! !!!" << endl;
+ }
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::InitOutputBranch(){
+cout << "start_InitOutput" << endl;
+ RootOutput::getInstance()->GetTree()->Branch("hasAuTarget",&hasAuTarget);
+
+
+ //RootOutput::getInstance()->GetTree()->Branch("Ex",&Ex,"Ex/D");
+ RootOutput::getInstance()->GetTree()->Branch("Ex",&Ex);
+ RootOutput::getInstance()->GetTree()->Branch("ExCalib",&ExCalib);
+ //RootOutput::getInstance()->GetTree()->Branch("EDC",&EDC,"EDC/D");
+ RootOutput::getInstance()->GetTree()->Branch("EDC_VAMOS",&EDC_VAMOS);
+ RootOutput::getInstance()->GetTree()->Branch("EDC_MUGAST",&EDC_MUGAST);
+ RootOutput::getInstance()->GetTree()->Branch("beta_VAMOS",&beta_VAMOS);
+ RootOutput::getInstance()->GetTree()->Branch("beta_MUGAST",&beta_MUGAST);
+ RootOutput::getInstance()->GetTree()->Branch("ToF",&ToF);
+
+ RootOutput::getInstance()->GetTree()->Branch("AddBack_EDC",&AddBack_EDC);
+ RootOutput::getInstance()->GetTree()->Branch("AddBack_EDC_Event1",&AddBack_EDC_Event1);
+ RootOutput::getInstance()->GetTree()->Branch("AddBack_EDC_Event2",&AddBack_EDC_Event2);
+ RootOutput::getInstance()->GetTree()->Branch("gammaEDC",&gammaEDC);
+ RootOutput::getInstance()->GetTree()->Branch("trackEDC",&trackEDC);
+ RootOutput::getInstance()->GetTree()->Branch("vamosEDC",&vamosEDC);
+ RootOutput::getInstance()->GetTree()->Branch("goldEDC",&goldEDC);
+ RootOutput::getInstance()->GetTree()->Branch("ThetaG",&ThetaG);
+ RootOutput::getInstance()->GetTree()->Branch("radius",&radius);
+ RootOutput::getInstance()->GetTree()->Branch("PhiG",&PhiG);
+ RootOutput::getInstance()->GetTree()->Branch("AGATA_Theta",&AGATA_Theta);
+ RootOutput::getInstance()->GetTree()->Branch("AGATA_Phi",&AGATA_Phi);
+ RootOutput::getInstance()->GetTree()->Branch("AGATA_GammaPx",&AGATA_GammaPx);
+ RootOutput::getInstance()->GetTree()->Branch("AGATA_GammaPy",&AGATA_GammaPy);
+ RootOutput::getInstance()->GetTree()->Branch("AGATA_GammaPz",&AGATA_GammaPz);
+ RootOutput::getInstance()->GetTree()->Branch("AGATA_GammaE",&AGATA_GammaE);
+ RootOutput::getInstance()->GetTree()->Branch("AGATA_OrigBetaX",&AGATA_OrigBetaX);
+ RootOutput::getInstance()->GetTree()->Branch("AGATA_OrigBetaY",&AGATA_OrigBetaY);
+ RootOutput::getInstance()->GetTree()->Branch("AGATA_OrigBetaZ",&AGATA_OrigBetaZ);
+ RootOutput::getInstance()->GetTree()->Branch("AGATA_OrigBetaMag",&AGATA_OrigBetaMag);
+ 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);
+ RootOutput::getInstance()->GetTree()->Branch("ELoss_Al",&ELoss_Al);
+ RootOutput::getInstance()->GetTree()->Branch("ELoss_Target",&ELoss_Target);
+ RootOutput::getInstance()->GetTree()->Branch("ELoss",&ELoss);
+ RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab);
+ RootOutput::getInstance()->GetTree()->Branch("ThetaGamma",&ThetaGamma);
+ RootOutput::getInstance()->GetTree()->Branch("Theta4",&Theta4);
+ RootOutput::getInstance()->GetTree()->Branch("Phi4",&Phi4);
+ RootOutput::getInstance()->GetTree()->Branch("Energy4",&Energy4);
+ RootOutput::getInstance()->GetTree()->Branch("ThetaGR",&ThetaGR);
+ RootOutput::getInstance()->GetTree()->Branch("PhiLab",&PhiLab);
+ RootOutput::getInstance()->GetTree()->Branch("HeavyThetaLab",&HeavyThetaLab);
+ RootOutput::getInstance()->GetTree()->Branch("HeavyPhiLab",&HeavyPhiLab);
+ RootOutput::getInstance()->GetTree()->Branch("HeavyPosXAtWindow",&HeavyPosXAtWindow);
+ RootOutput::getInstance()->GetTree()->Branch("HeavyPosYAtWindow",&HeavyPosYAtWindow);
+ RootOutput::getInstance()->GetTree()->Branch("HeavyPosXAtMUST2",&HeavyPosXAtMUST2);
+ RootOutput::getInstance()->GetTree()->Branch("HeavyPosYAtMUST2",&HeavyPosYAtMUST2);
+ RootOutput::getInstance()->GetTree()->Branch("HeavyAngX",&HeavyAngX);
+ RootOutput::getInstance()->GetTree()->Branch("HeavyAngY",&HeavyAngY);
+ 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("MGXStrip",&MGXStrip);
+ RootOutput::getInstance()->GetTree()->Branch("MGYStrip",&MGYStrip);
+ RootOutput::getInstance()->GetTree()->Branch("MGDet",&MGDet);
+ 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");
+ RootOutput::getInstance()->GetTree()->Branch("T_FPMW_HF_C",&T_FPMW_HF_C,"T_FPMW_CATS1/F");
+ RootOutput::getInstance()->GetTree()->Branch ("IC",IC, "IC[12]/F");
+ RootOutput::getInstance()->GetTree()->Branch ("ICRawPUMult01",ICRawPUMult01, "ICRawPUMult01[7]/s");
+
+ RootOutput::getInstance()->GetTree()->Branch("T_CATS2_HF",&T_CATS2_HF,"T_CATS2_HF/s");
+ RootOutput::getInstance()->GetTree()->Branch("T_MUGAST_FPMW",&T_MUGAST_FPMW,"T_MUGAST_FPMW/s");
+ RootOutput::getInstance()->GetTree()->Branch("T_FPMW_CATS2",&T_FPMW_CATS2,"T_FPMW_CATS2/s");
+ RootOutput::getInstance()->GetTree()->Branch("T_FPMW_HF",&T_FPMW_HF,"T_FPMW_HF/s");
+
+ // DC0
+ RootOutput::getInstance()->GetTree()->Branch ("DC0_X", &DC0_X, "DC0_X/F");
+ RootOutput::getInstance()->GetTree()->Branch ("DC0_Y", &DC0_Y, "DC0_Y/F");
+ // DC1
+ RootOutput::getInstance()->GetTree()->Branch ("DC1_X", &DC1_X, "DC1_X/F");
+ RootOutput::getInstance()->GetTree()->Branch ("DC1_Y", &DC1_Y, "DC1_Y/F");
+ // DC2
+ RootOutput::getInstance()->GetTree()->Branch ("DC2_X", &DC2_X, "DC2_X/F");
+ RootOutput::getInstance()->GetTree()->Branch ("DC2_Y", &DC2_Y, "DC2_Y/F");
+ // DC3
+ RootOutput::getInstance()->GetTree()->Branch ("DC3_X", &DC3_X, "DC3_X/F");
+ RootOutput::getInstance()->GetTree()->Branch ("DC3_Y", &DC3_Y, "DC3_Y/F");
+
+ RootOutput::getInstance()->GetTree()->Branch ("Xf", &Xf, "Xf/F");
+ RootOutput::getInstance()->GetTree()->Branch ("Tf", &Tf, "Tf/F");
+
+ RootOutput::getInstance()->GetTree()->Branch ("Brho", &Brho, "Brho/F");
+ RootOutput::getInstance()->GetTree()->Branch ("Theta", &Theta, "Theta/F");
+ RootOutput::getInstance()->GetTree()->Branch ("Phi", &Phi, "Phi/F");
+ RootOutput::getInstance()->GetTree()->Branch ("Path", &Path, "Path/F");
+
+ RootOutput::getInstance()->GetTree()->Branch ("EWIRE_1_1", &EWIRE_1_1, "EWIRE_1_1/s");
+ RootOutput::getInstance()->GetTree()->Branch ("EWIRE_1_2", &EWIRE_1_2, "EWIRE_1_2/s");
+ 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 ("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");
+ RootOutput::getInstance()->GetTree()->Branch("mBeta",&mBeta,"mBeta/D");
+ RootOutput::getInstance()->GetTree()->Branch("mGamma",&mGamma,"mGamma/D");
+ RootOutput::getInstance()->GetTree()->Branch("mM_Q",&mM_Q,"mM_Q/D");
+ RootOutput::getInstance()->GetTree()->Branch("mM",&mM,"mM/D");
+ RootOutput::getInstance()->GetTree()->Branch("mE",&mE,"mE/D");
+ RootOutput::getInstance()->GetTree()->Branch("mdE",&mdE,"mdE/D");
+*/
+
+ // DC0->3_X/_Y Xf, Tf, Brho, Path, EWIRE_1_1 all, MW_Nr, MW_T, AGAVA_VAMOSTS,
+ // Agata
+ // Time stamp of the agata trigger
+
+ RootOutput::getInstance()->GetTree()->Branch("TStrack",&TStrack,"TStrack/l");
+
+ // Array of reconstructed tracks
+ RootOutput::getInstance()->GetTree()->Branch("nbTrack",&nbTrack,"nbTrack/I");
+ RootOutput::getInstance()->GetTree()->Branch("trackE",trackE,"trackE[nbTrack]/F");
+ RootOutput::getInstance()->GetTree()->Branch("trackX1",trackX1,"trackX1[nbTrack]/F");
+ RootOutput::getInstance()->GetTree()->Branch("trackY1",trackY1,"trackY1[nbTrack]/F");
+ RootOutput::getInstance()->GetTree()->Branch("trackZ1",trackZ1,"trackZ1[nbTrack]/F");
+ RootOutput::getInstance()->GetTree()->Branch("trackT",trackT,"trackT[nbTrack]/F");
+ 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");
+
+ RootOutput::getInstance()->GetTree()->Branch("EventWeight",&EventWeight,"EventWeight/D");
+
+ //AGATA hits for the tracking optimization
+ RootOutput::getInstance()->GetTree()->Branch("nbHits",&nbHits,"nbHits/I");
+ RootOutput::getInstance()->GetTree()->Branch("hitE",hitE,"hitE[nbHits]/F");
+ RootOutput::getInstance()->GetTree()->Branch("hitX",hitX,"hitX[nbHits]/F");
+ RootOutput::getInstance()->GetTree()->Branch("hitY",hitY,"hitY[nbHits]/F");
+ RootOutput::getInstance()->GetTree()->Branch("hitZ",hitZ,"hitZ[nbHits]/F");
+ RootOutput::getInstance()->GetTree()->Branch("hitId",hitId,"hitId[nbHits]/I");
+ RootOutput::getInstance()->GetTree()->Branch("hitSg",hitSg,"hitSg[nbHits]/I");
+// mEDC_PhiAgata = new TH2I("mEDC_PhiAgata",
+// "Optimization DC for Phi AGATA",
+// 360,-0.5,359.5,
+// 3000,0,3000);
+// RootOutput::getInstance()->GetList()->Add(mEDC_PhiAgata);
+// mEDC_OffsetAgata = new TH2I("mEDC_OffsetAgata",
+// "Optimization DC for Offset AGATA",
+// 21,-10.5,10.5,
+// 3000,0,3000);
+// RootOutput::getInstance()->GetList()->Add(mEDC_OffsetAgata);
+// backingThick = new TH2I("backingThick",
+// "Optimization Backing Thickness",
+// 30,10.0,13.0,
+// 3000,0,3000);
+// RootOutput::getInstance()->GetList()->Add(backingThick);
+// kineShift = new TH3I("kineShift","Kinematic for Annular Shift",
+// 11,-0.5,10.5,
+// 1200,0,12,
+// 360,90,180);
+// RootOutput::getInstance()->GetList()->Add(kineShift);
+// for(int i=0; i<7; i++){
+// beamPos[i] = new TH3I(Form("beamPos_%i",i),
+// "Position of the beam",
+// 21,-10.5,10.5,
+// 21,-10.5,10.5,
+// 280,-1,6);
+// RootOutput::getInstance()->GetList()->Add(beamPos[i]);
+// }
+
+ if(isSim && !isPhaseSpace){
+ RootOutput::getInstance()->GetTree()->Branch("OriginalELab",
+ &OriginalELab,"OriginalELab/D");
+ RootOutput::getInstance()->GetTree()->Branch("OriginalThetaLab",
+ &OriginalThetaLab,"OriginalThetaLab/D");
+ RootOutput::getInstance()->GetTree()->Branch("BeamEnergy",
+ &BeamEnergy,"BeamEnergy/D");
+ }
+cout << "end_InitOutput" << endl;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::InitInputBranch(){
+ SetBranchStatus();
+cout << "start_InitInputBranch (after set status)" << endl;
+ // RootInput:: getInstance()->GetChain()->SetBranchAddress("GATCONF",&vGATCONF);
+ //
+ // Vamos
+ RootInput::getInstance()->GetChain()->SetBranchAddress("LTS",<S);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("VTS",&VTS);
+ 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);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("IC",IC);
+
+ RootInput::getInstance()->GetChain()->SetBranchAddress("T_CATS2_HF",&T_CATS2_HF);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("T_MUGAST_FPMW",&T_MUGAST_FPMW);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("T_FPMW_CATS2",&T_FPMW_CATS2);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("T_FPMW_HF",&T_FPMW_HF);
+
+ // DC0
+ RootInput::getInstance()->GetChain()->SetBranchAddress("DC0_X", &DC0_X);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("DC0_Y", &DC0_Y);
+ // DC1
+ RootInput::getInstance()->GetChain()->SetBranchAddress("DC1_X", &DC1_X);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("DC1_Y", &DC1_Y);
+ // DC2
+ RootInput::getInstance()->GetChain()->SetBranchAddress("DC2_X", &DC2_X);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("DC2_Y", &DC2_Y);
+ // DC3
+ RootInput::getInstance()->GetChain()->SetBranchAddress("DC3_X", &DC3_X);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("DC3_Y", &DC3_Y);
+
+ RootInput::getInstance()->GetChain()->SetBranchAddress("Xf", &Xf);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("Tf", &Tf);
+
+ RootInput::getInstance()->GetChain()->SetBranchAddress("Brho", &Brho);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("Theta", &Theta);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("Phi", &Phi);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("Path", &Path);
+
+ RootInput::getInstance()->GetChain()->SetBranchAddress("EWIRE_1_1", &EWIRE_1_1);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("EWIRE_1_2", &EWIRE_1_2);
+ 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);
+
+ // Agata
+
+ RootInput::getInstance()->GetChain()->SetBranchAddress("TStrack",&TStrack);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("nbTrack",&nbTrack);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("trackE",trackE);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("trackX1",trackX1);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("trackY1",trackY1);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("trackZ1",trackZ1);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("trackT",trackT);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("trackCrystalID",trackCrystalID);
+
+ // We need to add the hits for the tracking optimization
+ RootInput::getInstance()->GetChain()->SetBranchAddress("nbHits",&nbHits);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("hitE",hitE);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("hitX",hitX);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("hitY",hitY);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("hitZ",hitZ);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("hitId",hitId);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("hitSg",hitSg);
+
+ 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(isPhaseSpace){
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("EventWeight",&EventWeight);
+ }
+
+ if(isSim && !isPhaseSpace){
+ //RootInput:: getInstance()->GetChain()->SetBranchStatus("InitialConditions",true );
+ //RootInput:: getInstance()->GetChain()->SetBranchStatus("fIC_*",true );
+ RootInput::getInstance()->GetChain()->SetBranchAddress("InitialConditions",
+ &Initial);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("InteractionCoordinates",
+ &Interaction);
+ //RootInput:: getInstance()->GetChain()->SetBranchStatus("ReactionConditions",true );
+ //RootInput:: getInstance()->GetChain()->SetBranchStatus("fRC_*",true );
+ RootInput::getInstance()->GetChain()->SetBranchAddress("ReactionConditions",
+ &ReactionConditions);
+ }
+cout << "end_InitInputBranch (after set status)" << endl;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::SetBranchStatus(){
+cout << "start_SetBranchStatus" << endl;
+ // Set Branch status
+ RootInput::getInstance()->GetChain()->SetBranchStatus("LTS",true);
+ RootInput::getInstance()->GetChain()->SetBranchStatus("VTS",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 );
+ RootInput::getInstance()->GetChain()->SetBranchStatus("T_FPMW_HF_C",true);
+
+ RootInput::getInstance()->GetChain()->SetBranchStatus("T_CATS2_HF",true );
+ RootInput::getInstance()->GetChain()->SetBranchStatus("T_MUGAST_FPMW",true );
+ RootInput::getInstance()->GetChain()->SetBranchStatus("T_FPMW_CATS2",true );
+
+ // DC0
+ RootInput::getInstance()->GetChain()->SetBranchStatus("DC0_X",true);
+ RootInput::getInstance()->GetChain()->SetBranchStatus("DC0_Y",true);
+ // DC1
+ RootInput::getInstance()->GetChain()->SetBranchStatus("DC1_X",true);
+ RootInput::getInstance()->GetChain()->SetBranchStatus("DC1_Y",true);
+ // DC2
+ RootInput::getInstance()->GetChain()->SetBranchStatus("DC2_X",true);
+ RootInput::getInstance()->GetChain()->SetBranchStatus("DC2_Y",true);
+ // DC3
+ RootInput::getInstance()->GetChain()->SetBranchStatus("DC3_X",true);
+ RootInput::getInstance()->GetChain()->SetBranchStatus("DC3_Y",true);
+
+ RootInput::getInstance()->GetChain()->SetBranchStatus("Xf",true );
+ RootInput::getInstance()->GetChain()->SetBranchStatus("Tf",true );
+
+ RootInput::getInstance()->GetChain()->SetBranchStatus("Brho",true );
+ RootInput::getInstance()->GetChain()->SetBranchStatus("Path",true );
+
+ RootInput::getInstance()->GetChain()->SetBranchStatus("EWIRE_1_1",true );
+ RootInput::getInstance()->GetChain()->SetBranchStatus("EWIRE_1_2",true );
+ 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("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("TStrack",true );
+ RootInput::getInstance()->GetChain()->SetBranchStatus("nbTrack",true );
+ RootInput::getInstance()->GetChain()->SetBranchStatus("trackE",true );
+ RootInput::getInstance()->GetChain()->SetBranchStatus("trackX1",true );
+ RootInput::getInstance()->GetChain()->SetBranchStatus("trackY1",true );
+ RootInput::getInstance()->GetChain()->SetBranchStatus("trackZ1",true );
+ RootInput::getInstance()->GetChain()->SetBranchStatus("trackT",true );
+ RootInput::getInstance()->GetChain()->SetBranchStatus("trackCrystalID",true);
+
+ RootInput::getInstance()->GetChain()->SetBranchStatus("hitE",1);
+ RootInput::getInstance()->GetChain()->SetBranchStatus("hitX",1);
+ RootInput::getInstance()->GetChain()->SetBranchStatus("hitY",1);
+ RootInput::getInstance()->GetChain()->SetBranchStatus("hitZ",1);
+ RootInput::getInstance()->GetChain()->SetBranchStatus("hitId",1);
+ RootInput::getInstance()->GetChain()->SetBranchStatus("hitSg",1);
+ RootInput::getInstance()->GetChain()->SetBranchStatus("nbHits",1);
+
+ 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(isPhaseSpace){
+ RootInput:: getInstance()->GetChain()->SetBranchStatus("EventWeight",true );
+ }
+
+ if(isSim && !isPhaseSpace){
+ RootInput:: getInstance()->GetChain()->SetBranchStatus("InitialConditions",true );
+ RootInput:: getInstance()->GetChain()->SetBranchStatus("fIC_*",true );
+ RootInput:: getInstance()->GetChain()->SetBranchStatus("InteractionCoordinates",true );
+ RootInput:: getInstance()->GetChain()->SetBranchStatus("ReactionConditions",true );
+ RootInput:: getInstance()->GetChain()->SetBranchStatus("fRC_*",true );
+ }
+cout << "end_SetBranchStatus" << endl;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::ReInitValue(){
+//cout << "start_ReInit" << endl;
+ hasAuTarget.clear();
+
+ Ex.clear();
+ ExCalib.clear();
+ Ecm.clear();
+ AddBack_EDC.clear();
+ AddBack_EDC_Event1.clear();
+ AddBack_EDC_Event2.clear();
+ AGATA_Theta.clear();
+ AGATA_Phi.clear();
+ AGATA_GammaPx.clear();
+ AGATA_GammaPy.clear();
+ AGATA_GammaPz.clear();
+ AGATA_GammaE.clear();
+ AGATA_OrigBetaX.clear();
+ AGATA_OrigBetaY.clear();
+ AGATA_OrigBetaZ.clear();
+ AGATA_OrigBetaMag.clear();
+ EAgata = -1000;
+ ELab.clear();
+ RawEnergy.clear();
+ ELoss_Al.clear();
+ ELoss_Target.clear();
+ ELoss.clear();
+ ThetaLab.clear();
+ PhiLab.clear();
+ ThetaCM.clear();
+ HeavyThetaLab.clear();
+ HeavyPhiLab.clear();
+ HeavyPosXAtWindow.clear();
+ HeavyPosYAtWindow.clear();
+ HeavyPosXAtMUST2.clear();
+ HeavyPosYAtMUST2.clear();
+ HeavyAngX.clear();
+ HeavyAngY.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;
+
+ filledCline=false;
+
+ MG_T=-1000;
+ MG_E=-1000;
+ MG_X=-1000;
+ MG_Y=-1000;
+ MG_D=-1000;
+
+ //EDC= -1000;
+ EDC.clear();
+
+ /* Just straight copied from IZanons */
+ Ex.clear();
+ ExCalib.clear();
+ //AddBack_EDC=-1000;
+ EAgata = -1000;
+ ELab.clear();
+ MGDet.clear();
+ MGXStrip.clear();
+ MGYStrip.clear();
+
+ ThetaG.clear();
+ radius = 0;
+ PhiG.clear();
+ ThetaLab.clear();
+ ThetaGamma.clear();
+ Theta4 = 0;
+ Phi4 = 0;
+ Theta3.clear();
+ Phi3.clear();
+ ThetaGR.clear();
+ PhiLab.clear();
+ ThetaCM.clear();
+ gammaEDC.clear();
+ trackEDC.clear();
+ vamosEDC.clear();
+ goldEDC.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;
+
+ EDC_VAMOS.clear();
+ EDC_MUGAST.clear();
+ beta_VAMOS.clear();
+ beta_MUGAST.clear();
+
+//cout << "end_ReInit" << endl;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::CalculateVamos(){
+//cout << "start_CalculateVamos" << endl;
+ mT1 = -1500.;
+ mT = -1500.;
+ mV = -1500.;
+ mD2 = -1500;
+ mBeta = -1500.;
+ mGamma = -1500.;
+ mM_Q = -1500.;
+ mM = -1500.;
+ mQ = -1500.;
+ mMrec = -1500.;
+ mdE = -1500;
+ mE = -1500;
+
+ mT1 = T_FPMW_HF_C;
+ mT = 596.0 - mT1;
+
+ Float_t ICTh = 1;
+ if(MW_Nr == 1){
+ switch (MW_N[0]) {
+ case 0:
+ mD2 = 826.465 ;
+ break;
+ case 1:
+ mD2 = 825.677 ;
+ break;
+ case 2:
+ mT += -2.5;
+ mD2 = 830.012 ;
+ break;
+ case 3:
+ mT += -2.0;
+ mD2 = 825.939 ;
+ break;
+ case 4:
+ mD2 = 825.521 ;
+ break;
+ case 5:
+ mT += 1;
+ mD2 = 827.253 ;
+ break;
+ case 6:
+ mT += 2;
+ mD2 = 826.334 ;
+ break;
+ case 7:
+ mT += 2;
+ mD2 = 827.253 ;
+ break;
+ case 8:
+ mT += 1.5;
+ mD2 = 833.559 ;
+ break;
+ case 9:
+ mT += 0;
+ mD2 = 825.043 ;
+ break;
+ case 10:
+ mT += 0.9;
+ mD2 = 826.728 ;
+ break;
+ case 11:
+ mT += 2.2;
+ mD2 = 825.414 ;
+ break;
+ case 12:
+ mT += 5;
+ mD2 = 817.138 ;
+ break;
+ case 13:
+ mT += 5;
+ mD2 = 823.049 ;
+ break;
+ case 14:
+ mT += 4;
+ mD2 = 827.779 ;
+ break;
+ case 15:
+ mT += 3 ;
+ mD2 = 832.508 ;
+ break;
+ case 16:
+ mT += 2;
+ mD2 = 841.966 ;
+ break;
+ case 17:
+ mD2 = 848.403 ;
+ break;
+ case 18:
+ mD2 = 831.194 ;
+ break;
+ case 19:
+ mD2 = 828.698 ;
+ break;
+ case 20:
+ mD2 = 824.232 ;
+ break;
+ default:
+ mT = -1500.;
+ mD2 = -1500 ;
+ break;
+ }
+ }
+ else{
+ mT = -1500.;
+ mD2 = -1500 ;
+ }
+
+ mV = mD2/mT;
+ mBeta = mV/29.9792;
+ mGamma = 1./sqrt(1.0-mBeta*mBeta);
+ mM_Q = Brho/3.105/mBeta/mGamma;
+
+ mE = 0.86*(2.8*IC[0]+IC[1]*(IC[0]>0.1)+IC[2]*(IC[0]>0.1)*(IC[1]>0.1)+IC[3]*(IC[0]>0.1)*(IC[1]>0.1)*(IC[2]>0.1)+IC[4]*(IC[0]>0.1)*(IC[1]>0.1)*(IC[2]>0.1)*(IC[3]>0.1)+IC[5]*(IC[0]>0.1)*(IC[1]>0.1)*(IC[2]>0.1)*(IC[3]>0.1)*(IC[4]>0.1)+IC[6]*(IC[0]>0.1)*(IC[1]>0.1)*(IC[2]>0.1)*(IC[3]>0.1)*(IC[4]>0.1)*(IC[5]>0.1));
+
+ mdE = IC[0]+IC[1]*(IC[0]>0.1);
+
+ mM = mE/931.494/(mGamma-1.);
+ mQ = mM/mM_Q;
+ mMrec = int(mQ+0.5)*mM_Q;
+
+//cout << "end_CalculateVamos" << endl;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+// 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);
+ }
+ };
+
+ proxy_analysis p_analysis;
+}
+
+
diff --git a/Projects/e775s/analysis/Analysis.h b/Projects/e775s/analysis/Analysis.h
new file mode 100644
index 0000000000000000000000000000000000000000..07ec1224f724d0a7120d6bc473f273a0e24d6535
--- /dev/null
+++ b/Projects/e775s/analysis/Analysis.h
@@ -0,0 +1,426 @@
+#ifndef Analysis_h
+#define Analysis_h
+/*****************************************************************************
+ * Copyright (C) 2009-2014 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Adrien MATTA contact address: a.matta@surrey.ac.uk *
+ * Edited for e793s *
+ * *
+ * Creation Date : *
+ * Last update : May 2021 *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * Class describing the property of an Analysis object *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+#include"NPVAnalysis.h"
+#include"NPEnergyLoss.h"
+#include"NPReaction.h"
+#include"NPBeam.h"
+#include"RootOutput.h"
+#include"RootInput.h"
+#include "TF1.h"
+#include "TInitialConditions.h"
+#include "TReactionConditions.h"
+#include "TInteractionCoordinates.h"
+#include "TMust2Physics.h"
+#include "TMugastPhysics.h"
+#include "TCATSPhysics.h"
+//#include "TPlasticPhysics.h"
+//#include "../../NPLib/Detectors/ModularLeaf/TModularLeafPhysics.h"
+#include "TModularLeafPhysics.h"
+#include <TRandom3.h>
+#include <TVector3.h>
+#include <TMath.h>
+#include <bitset>
+
+int NumThetaAngleBins = 1800;// 180 = 1 deg, 360 = 0.5 deg, 900 = 0.2 deg, 1800 = 0.1 deg
+
+auto ThetaCM_emmitted = new TH1F("ThetaCM_emmitted","ThetaCM_emmitted",NumThetaAngleBins,0,180);
+auto ThetaCM_detected_MG = new TH1F("ThetaCM_detected_MG","ThetaCM_detected_MG",NumThetaAngleBins,0,180);
+auto ThetaCM_detected_MG_wV = new TH1F("ThetaCM_detected_MG_wV","ThetaCM_detected_MG_wV",NumThetaAngleBins,0,180);
+auto ThetaCM_detected_MG_Unb = new TH1F("ThetaCM_detected_MG_UNb","ThetaCM_detected_MG_Unb",NumThetaAngleBins,0,180);
+
+auto ThetaLab_emmitted = new TH1F("ThetaLab_emmitted","ThetaLab_emmitted",NumThetaAngleBins,0,180);
+auto ThetaLab_detected_MG = new TH1F("ThetaLab_detected_MG","ThetaLab_detected_MG",NumThetaAngleBins,0,180);
+auto ThetaLab_detected_MG_wV = new TH1F("ThetaLab_detected_MG_wV","ThetaLab_detected_MG_wV",NumThetaAngleBins,0,180);
+auto ThetaLab_detected_MG_Unb = new TH1F("ThetaLab_detected_MG_Unb","ThetaLab_detected_MG_Unb",NumThetaAngleBins,0,180);
+//auto HistCline_MM = new TH1F("HistCline_MM","HistClint_MM",NumThetaAngleBins,0,180);
+//auto HistCline_MG = new TH1F("HistCline_MG","HistClint_MG",NumThetaAngleBins,0,180);
+
+double degtorad = M_PI/180.;
+vector<double> clineVal, clineX;
+bool filledCline;
+
+TH1F *ThetaCM_detected_MGX[7];
+TH1F *ThetaLab_detected_MGX[7];
+
+TH1F *ThetaCM_detected_MGX_wV[7];
+TH1F *ThetaLab_detected_MGX_wV[7];
+
+TH1F *ThetaCM_detected_MGX_Unb[7];
+TH1F *ThetaLab_detected_MGX_Unb[7];
+
+class Analysis: public NPL::VAnalysis{
+ public:
+ Analysis();
+ ~Analysis();
+
+ public:
+ void Init();
+ void TreatEvent();
+ void End();
+
+ void InitOutputBranch();
+ void InitInputBranch();
+ void SetBranchStatus();
+ void CalculateVamos();
+ void ReInitValue();
+ static NPL::VAnalysis* Construct();
+
+ private:
+ bool runningAuTarget;
+ std::vector<bool> hasAuTarget;
+
+ unsigned int ParticleMult;
+ unsigned int GammaMult;
+ //double EDC;
+ std::vector<double> EDC;
+ vector<double> AddBack_EDC;
+ vector<double> AddBack_EDC_Event1;
+ vector<double> AddBack_EDC_Event2;
+ vector<double> AGATA_Theta;
+ vector<double> AGATA_Phi;
+ vector<double> AGATA_GammaPx;
+ vector<double> AGATA_GammaPy;
+ vector<double> AGATA_GammaPz;
+ vector<double> AGATA_GammaE;
+ vector<double> AGATA_OrigBetaX;
+ vector<double> AGATA_OrigBetaY;
+ vector<double> AGATA_OrigBetaZ;
+ vector<double> AGATA_OrigBetaMag;
+ double EAgata;
+ std::vector<double> ELab;
+ std::vector<double> Ex;
+ std::vector<double> Ecm;
+ std::vector<double> RawEnergy;
+ std::vector<double> ELoss_Al;
+ std::vector<double> ELoss_Target;
+ std::vector<double> ELoss;
+ std::vector<double> ThetaLab;
+ std::vector<double> PhiLab;
+ std::vector<double> ThetaCM;
+ std::vector<double> HeavyThetaLab;
+ std::vector<double> HeavyPhiLab;
+ std::vector<double> HeavyPosXAtWindow;
+ std::vector<double> HeavyPosYAtWindow;
+ std::vector<double> HeavyPosXAtMUST2;
+ std::vector<double> HeavyPosYAtMUST2;
+ std::vector<double> HeavyAngX;
+ std::vector<double> HeavyAngY;
+
+ double WindowDistance = 681.33;
+ double MUST2InnerDistance = 391.4;
+
+ //For use in simulations
+ double OriginalELab;
+ double OriginalThetaLab;
+ double BeamEnergy;
+
+ NPL::Reaction reaction;
+ //NPL::Reaction *reaction;
+ // Energy loss table: the G4Table are generated by the simulation
+ NPL::EnergyLoss LightTarget;
+ NPL::EnergyLoss LightAl;
+ NPL::EnergyLoss LightSi;
+ NPL::EnergyLoss LightAu;
+ NPL::EnergyLoss BeamTarget;
+ NPL::EnergyLoss HeavyTarget;
+
+ double TargetThickness ;
+ double WindowsThickness;
+ // Beam Energy
+ double OriginalBeamEnergy ; // AMEV
+ double FinalBeamEnergy;
+
+ // Beam Position
+ double XBeam;
+ double YBeam;
+
+ // intermediate variable
+ TVector3 BeamDirection;
+ TVector3 BeamImpact;
+ TRandom3 Rand;
+ int Run;
+ int DetectorNumber;
+ double ThetaNormalTarget;
+ double ThetaM2Surface;
+ double ThetaMGSurface;
+ double elab_tmp;
+ double thetalab_tmp;
+ double philab_tmp;
+ double Si_E_M2;
+ double CsI_E_M2;
+ double Energy;
+ #define GATCONF_SIZE 16
+ unsigned int GATCONF_Counter[GATCONF_SIZE];
+ string GATCONF_Label[GATCONF_SIZE] = {
+ "MUVI1",
+ "MUVI2",
+ "MUVI3",
+ "",
+ "CATS2%",
+ "VAMOS%",
+ "MUST2%",
+ "",
+ "",
+ "VAMOS",
+ "ORAGATA",
+ } ;
+
+ double MG_T;
+ double MG_E;
+ int MG_X;
+ int MG_Y;
+ int MG_D;
+
+ double ThetaGDSurface ;
+ std::vector<double> X ;
+ std::vector<double> Y ;
+ std::vector<double> Z ;
+ int AHeavy;
+ int ALight;
+ double MHeavy;
+ double MLight;
+
+ // Vamos Branches
+ unsigned long long int LTS;
+ unsigned short T_FPMW_CATS1;
+ float T_FPMW_CATS2_C;
+ unsigned short T_FPMW_HF;
+ unsigned short T_CATS2_HF;
+ unsigned short T_MUGAST_FPMW;
+ unsigned short T_FPMW_CATS2;
+ float T_FPMW_HF_C;
+ float IC[12];
+ unsigned short ICRawPUMult01[7];
+
+ float DC0_X;
+ float DC0_Y;
+
+ float DC1_X;
+ float DC1_Y;
+
+ float DC2_X;
+ float DC2_Y;
+
+ float DC3_X;
+ float DC3_Y;
+
+ //float Xf;
+ float Tf;
+
+ float Brho;
+ float Theta;
+ float Phi;
+ float Path;
+
+ unsigned short EWIRE_1_1;
+ unsigned short EWIRE_1_2;
+ unsigned short EWIRE_2_1;
+ unsigned short EWIRE_2_2;
+
+ int MW_Nr;
+ float MW_T[1000];
+ float MWT[1000];
+ unsigned short MW_N[1000];
+ unsigned long long int AGAVA_VAMOSTS;
+
+ // Vamos Calculated
+ double BrhoRef;
+ double mT;
+ double mV;
+ double mD;
+ double mBeta;
+ double mGamma;
+ double mM_Q;
+ double mM;
+ double mE;
+ double mdE;
+
+ // Agata branches
+ double agata_zShift;
+ /*
+ unsigned long long int TStrack;
+ int nbHits;
+ int nbTrack;
+ float *trackE = new float(1000);
+ float *trackX1 = new float(1000);
+ float *trackY1 = new float(1000);
+ float *trackZ1 = new float(1000);
+ float *trackT = new float(1000);
+ int *trackCrystalID = new int(1000);
+ */
+ int nbAdd;
+ unsigned long long int TSHit;
+ float AddE[1000] ;
+ float AddX[1000] ;
+ float AddY[1000] ;
+ float AddZ[1000] ;
+
+ double EventWeight; //double or uns long long int?
+
+ //
+ double dE;
+ double dTheta;
+ // Branches and detectors
+ TMust2Physics* M2;
+ TMugastPhysics* MG;
+ //TCATSPhysics* CATS;
+ TModularLeafPhysics* ML;
+ //TPlasticPhysics* VM;
+
+ bool warning=true;
+
+
+ //Simulation Stuff
+ bool isSim;
+ bool isPhaseSpace;
+ bool excludePoor;
+ bool writetoscreen;
+ TInitialConditions* Initial;
+ TInteractionCoordinates* Interaction;
+ TReactionConditions* ReactionConditions;
+ bool hitVAMOS = false;
+
+ // Beam object
+ NPL::Beam* Beam;
+
+ unsigned int GATCONF_MASTER;
+
+ unsigned long long int count ;
+
+
+
+
+ /******** FROM IZANON ********/
+ /******** FROM IZANON ********/
+ /******** FROM IZANON ********/
+ /******** FROM IZANON ********/
+
+ private:
+ //double EDC;
+ std::vector<int> MGDet;
+ std::vector<int> MGXStrip;
+ std::vector<int> MGYStrip;
+ std::vector<double> EDC_VAMOS;
+ std::vector<double> EDC_MUGAST;
+ std::vector<double> beta_VAMOS;
+ std::vector<double> beta_MUGAST;
+ unsigned short ToF;
+
+ std::vector<double> gammaEDC;
+ std::vector<double> trackEDC;
+ std::vector<double> vamosEDC;
+ std::vector<double> goldEDC;
+ double recoilBeta;
+ std::vector<double> ThetaGR;
+ TVector3 recoilDir;
+ TVector3 gammaDir;
+ std::vector<double> ThetaG;
+ double radius;
+ std::vector<double> PhiG;
+ std::vector<double> ExCalib;
+ std::vector<double> ExTmp;
+ std::vector<double> ThetaGamma;
+ double Theta4;
+ double Phi4;
+ double Energy4;
+// std::vector<double> Theta4;
+// std::vector<double> Phi4;
+ std::vector<double> Theta3;
+ std::vector<double> Phi3;
+ std::vector<int> NTelescope;
+ std::vector<int> NXStrip;
+ std::vector<int> NYStrip;
+
+ NPL::Reaction tmp_reaction;
+ // Energy loss table: the G4Table are generated by the simulation
+ NPL::EnergyLoss HeavyAu;
+ NPL::EnergyLoss* BeamWindow;
+ NPL::EnergyLoss* LightWindow;
+
+ //DC Optimization
+ TH2I *mEDC_PhiAgata;
+ TH2I *mEDC_OffsetAgata;
+ TH2I *backingThick;
+//TH3I *beamPos[7];
+//TH3I *kineShift;
+
+ double FrontDeformation ;
+ double FrontRadius ;
+ double AuThickness;
+ // Beam Energy
+
+ // intermediate variable
+ double thetalab_2 ;
+ double E4_tmp;
+
+ double DSSD_T;
+ int TelNbr;
+ bool m_ValidEvent;
+
+ // Vamos Branches
+ unsigned long long int VTS;
+ float Xf;
+
+ // Vamos Calculated
+ double mT1;
+ double mD2;
+ double mQ;
+ double mMrec;
+
+ // Agata branches
+ unsigned long long int TStrack;
+ //int nbHits;
+ /*int nbTrack;
+ float *trackE = new float(1000);
+ float *trackX1 = new float(1000);
+ float *trackY1 = new float(1000);
+ float *trackZ1 = new float(1000);
+ float *trackT = new float(1000);
+ int *trackCrystalID = new int(1000);
+ */
+ int nbTrack;
+ float trackE[1000];
+ float trackX1[1000];
+ float trackY1[1000];
+ float trackZ1[1000];
+ float trackT[1000];
+ Int_t trackCrystalID[1000];
+
+ // Hit info for tracking opt
+ Int_t nbHits;
+ float hitE[1000] ;
+ float hitX[1000] ;
+ float hitY[1000] ;
+ float hitZ[1000] ;
+ Int_t hitId[1000];
+ Int_t hitSg[1000];
+ //
+ // Branches and detectors
+ TCATSPhysics* CATS;
+
+
+};
+#endif
+
diff --git a/Projects/e775s/analysis/InputFiles/All.runtotreat b/Projects/e775s/analysis/InputFiles/All.runtotreat
new file mode 100644
index 0000000000000000000000000000000000000000..eb2f79568b28ed6efd14fd09473dfc6a664d1cc6
--- /dev/null
+++ b/Projects/e775s/analysis/InputFiles/All.runtotreat
@@ -0,0 +1,43 @@
+TTreeName
+ TreeMaster
+RootFileName
+ /home/paxman/Experiment/e775s/RawData/run_0115/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0116/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0117/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0123/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0124/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0125/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0128/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0129/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0130/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0131/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0132/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0133/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0134/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0135/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0137/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0138/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0139/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0140/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0142/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0150/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0151/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0152/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0153/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0163/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0164/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0165/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0166/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0167/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0168/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0169/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0171/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0172/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0173/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0174/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0178/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0179/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0180/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0181/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0184/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0185/Out/Analysis/Tree_0000.root
diff --git a/Projects/e775s/analysis/InputFiles/Au.runtotreat b/Projects/e775s/analysis/InputFiles/Au.runtotreat
new file mode 100644
index 0000000000000000000000000000000000000000..c61ee1f14b06cb5c7a3a58014fee5d62289fedea
--- /dev/null
+++ b/Projects/e775s/analysis/InputFiles/Au.runtotreat
@@ -0,0 +1,35 @@
+TTreeName
+ TreeMaster
+RootFileName
+ /home/paxman/Experiment/e775s/RawData/run_0115/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0116/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0117/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0124/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0125/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0128/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0129/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0130/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0131/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0132/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0134/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0135/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0137/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0138/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0139/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0151/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0152/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0153/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0163/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0164/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0165/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0166/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0167/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0168/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0169/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0171/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0172/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0173/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0180/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0181/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0184/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0185/Out/Analysis/Tree_0000.root
diff --git a/Projects/e775s/analysis/InputFiles/CD.runtotreat b/Projects/e775s/analysis/InputFiles/CD.runtotreat
new file mode 100644
index 0000000000000000000000000000000000000000..d466cdac8310e397b8209fd6aaac63afdadb66d7
--- /dev/null
+++ b/Projects/e775s/analysis/InputFiles/CD.runtotreat
@@ -0,0 +1,9 @@
+TTreeName
+ TreeMaster
+RootFileName
+ /home/paxman/Experiment/e775s/RawData/run_0123/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0140/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0150/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0174/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0178/Out/Analysis/Tree_0000.root
+ /home/paxman/Experiment/e775s/RawData/run_0179/Out/Analysis/Tree_0000.root
diff --git a/Projects/e775s/analysis/InputFiles/Calibration.txt b/Projects/e775s/analysis/InputFiles/Calibration.txt
new file mode 100644
index 0000000000000000000000000000000000000000..8ace5e99ecf13c8cddc23e4217d518632dbae815
--- /dev/null
+++ b/Projects/e775s/analysis/InputFiles/Calibration.txt
@@ -0,0 +1,19 @@
+CalibrationFilePath
+/home/paxman/Programs/nptoolv3/Projects/e775s/analysis/InputFiles/CalibrationFiles/Cal_Str_X_E_MG1.cal
+/home/paxman/Programs/nptoolv3/Projects/e775s/analysis/InputFiles/CalibrationFiles/Cal_Str_X_E_MG2.cal
+/home/paxman/Programs/nptoolv3/Projects/e775s/analysis/InputFiles/CalibrationFiles/Cal_Str_X_E_MG3.cal
+/home/paxman/Programs/nptoolv3/Projects/e775s/analysis/InputFiles/CalibrationFiles/Cal_Str_X_E_MG4.cal
+/home/paxman/Programs/nptoolv3/Projects/e775s/analysis/InputFiles/CalibrationFiles/Cal_Str_X_E_MG5.cal
+/home/paxman/Programs/nptoolv3/Projects/e775s/analysis/InputFiles/CalibrationFiles/Cal_Str_X_E_MG6.cal
+/home/paxman/Programs/nptoolv3/Projects/e775s/analysis/InputFiles/CalibrationFiles/Cal_Str_X_E_MG7.cal
+/home/paxman/Programs/nptoolv3/Projects/e775s/analysis/InputFiles/CalibrationFiles/Cal_Str_X_E_MG11.cal
+/home/paxman/Programs/nptoolv3/Projects/e775s/analysis/InputFiles/CalibrationFiles/Cal_Str_Y_E_MG1.cal
+/home/paxman/Programs/nptoolv3/Projects/e775s/analysis/InputFiles/CalibrationFiles/Cal_Str_Y_E_MG2.cal
+/home/paxman/Programs/nptoolv3/Projects/e775s/analysis/InputFiles/CalibrationFiles/Cal_Str_Y_E_MG3.cal
+/home/paxman/Programs/nptoolv3/Projects/e775s/analysis/InputFiles/CalibrationFiles/Cal_Str_Y_E_MG4.cal
+/home/paxman/Programs/nptoolv3/Projects/e775s/analysis/InputFiles/CalibrationFiles/Cal_Str_Y_E_MG5.cal
+/home/paxman/Programs/nptoolv3/Projects/e775s/analysis/InputFiles/CalibrationFiles/Cal_Str_Y_E_MG6.cal
+/home/paxman/Programs/nptoolv3/Projects/e775s/analysis/InputFiles/CalibrationFiles/Cal_Str_Y_E_MG7.cal
+/home/paxman/Programs/nptoolv3/Projects/e775s/analysis/InputFiles/CalibrationFiles/Cal_Str_Y_E_MG11.cal
+/home/paxman/Programs/nptoolv3/Projects/e775s/analysis/InputFiles/CalibrationFiles/Cal_AllTime.cal
+
diff --git a/Projects/e775s/analysis/InputFiles/e775s_2024-07-30.reaction b/Projects/e775s/analysis/InputFiles/e775s_2024-07-30.reaction
new file mode 100644
index 0000000000000000000000000000000000000000..1b95bb10afc00ecb2475ed6fe0976f12d7825bb5
--- /dev/null
+++ b/Projects/e775s/analysis/InputFiles/e775s_2024-07-30.reaction
@@ -0,0 +1,35 @@
+%%%%%%%%%%%%%%%%%%%%%% E775S %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Beam
+ Particle= 19O
+ ExcitationEnergy= 0 MeV
+ Energy= 151.82 MeV
+ SigmaEnergy= 0. MeV
+ SigmaThetaX= 0.0 deg
+ SigmaPhiY= 0.0 deg
+ SigmaX= 1.5 mm
+ SigmaY= 3 mm
+ MeanThetaX= 0 deg
+ MeanPhiY= 0 deg
+ MeanX= -3.41 mm
+ MeanY= -0.20 mm
+ %EnergyProfilePath=
+ %XThetaXProfilePath=
+ %YPhiYProfilePath=
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+TwoBodyReaction
+ Beam= 19O
+ Target= 2H
+ Light= 1H
+ Heavy= 20O
+ ExcitationEnergyLight= 0.0 MeV
+ ExcitationEnergyHeavy= 0.0 MeV
+%% ExcitationEnergyHeavy= 1.67368 MeV
+%% 3.570 4.069 5.223
+%% Zanon used this, I'm using flatbackwards for now %% LabCrossSectionPath = 24.19Odp0 CSR
+ LabCrossSectionPath = flatbackwards.txt CSR
+
+ ShootLight= 1
+ ShootHeavy= 1
+
diff --git a/Projects/e775s/analysis/InputFiles/mugast_2024-07-30.detector b/Projects/e775s/analysis/InputFiles/mugast_2024-07-30.detector
new file mode 100644
index 0000000000000000000000000000000000000000..f45db46c3ee3d0763fc9d40bf6c838c42097616c
--- /dev/null
+++ b/Projects/e775s/analysis/InputFiles/mugast_2024-07-30.detector
@@ -0,0 +1,97 @@
+%%%%%%%%%%Detector%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Target
+ THICKNESS= 2.80 micrometer
+ ANGLE= 0 deg
+ RADIUS= 10 mm
+ MATERIAL= CD2
+ X= +0.00 mm
+ Y= +0.00 mm
+ Z= -0.77 mm
+ BackThickness= 12.6 micrometer
+ BackMaterial= Au
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Mugast Trapezoid
+ DetectorNumber= 1
+ X128_Y128= -23.05 42.83 -99.66 mm
+ X1_Y128= -40.74 25.21 -99.65 mm
+ X1_Y1= -120.65 58.67 -32.39 mm
+ X128_Y1= -55.66 123.17 -32.55 mm
+ VIS= all
+
+%%%%%%%%%%%%%%%%%%%%%POS TO UPDATE%%%%%%%%%
+Mugast Trapezoid
+ DetectorNumber= 2
+ X128_Y128= -45.68 14.28 -98.74 mm
+ X1_Y128= -45.68 -10.92 -98.74 mm
+ X1_Y1= -125.65 -44.06 -31.63 mm
+ X128_Y1= -125.65 47.42 -31.63 mm
+ VIS= all
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Mugast Trapezoid
+ DetectorNumber= 3
+ X128_Y128= -48.05 -23.23 -95.05 mm
+ X1_Y128= -24.83 -39.35 -99.33 mm
+ X1_Y1= -58.24 -119.46 -32.57 mm
+ X128_Y1= -122.86 -54.55 -32.41 mm
+ VIS= all
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Mugast Trapezoid
+ DetectorNumber= 4
+ X128_Y128= -15.11 -44.73 -99.25 mm
+ X1_Y128= 10.11 -44.71 -99.35 mm
+ X1_Y1= 43.37 -125.07 -32.74 mm
+ X128_Y1= -48.18 -125.07 -32.45 mm
+ VIS= all
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Mugast Trapezoid
+ DetectorNumber= 5
+ X128_Y128= 21.68 -41.48 -100.19 mm
+ X1_Y128= 39.55 -23.81 -100.11 mm
+ X1_Y1= 119.55 -57.08 -33.23 mm
+ X128_Y1= 54.81 -121.7 -33.12 mm
+ VIS= all
+
+%%%%%%%%%%%%%%%%%%%%POS TO UPDATE%%%%%%%%%
+Mugast Trapezoid
+ DetectorNumber= 6
+ X128_Y128= 45.68 -14.28 -98.74 mm
+ X1_Y128= 45.68 10.92 -98.74 mm
+ X1_Y1= 125.65 44.06 -31.63 mm
+ X128_Y1= 125.65 -47.42 -31.63 mm
+ VIS= all
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Mugast Trapezoid
+ DetectorNumber= 7
+ X128_Y128= 40.99 21.69 -98.21 mm
+ X1_Y128= 23.34 39.63 -98.04 mm
+ X1_Y1= 57.53 120.53 -32.58 mm
+ X128_Y1= 122.2 55.73 -32.67 mm
+ VIS= all
+
+%%%%%%%%%%%%%%%%%%%POS TO UPDATE%%%%%%%%%
+Mugast Trapezoid
+ DetectorNumber= 8
+ X128_Y128= -15.11 -44.73 -99.25 mm
+ X1_Y128= 10.11 -44.71 -99.35 mm
+ X1_Y1= 43.37 -125.07 -32.74 mm
+ X128_Y1= -48.18 -125.07 -32.45 mm
+ VIS= all
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Mugast Annular
+ DetectorNumber= 11
+ Center= -0.3 1 -126.9 mm
+ VIS= all
+
+ModularLeaf
+ DefaultValue= -5000
+ Leafs= GATCONF_MASTER GATCONF_SLAVE1 GATCONF_SLAVE2 CONFDEC_AGATA CONFDEC_VAMOS DATATRIG_CATS1 DATATRIG_CATS2 T_MUGAST_HF EVAMOS_AGATA E_AGATA T_FAG_HF T_VAMOS_AGATA PLQ_L PLQ_R T_AGATA_HF TVAMOS_MUGAST_HF T_MUGAST_VAMOS T_MUGAST_AGATA T_VAMOS_HF ADC1_9 ADC1_10 ADC1_11
+ Leafs= GATCONF_MASTER CONFDEC_AGATA CONFDEC_VAMOS T_MUGAST_HF EVAMOS_AGATA E_AGATA T_FAG_HF T_VAMOS_AGATA T_AGATA_HF TVAMOS_MUGAST_HF T_MUGAST_VAMOS T_MUGAST_AGATA T_VAMOS_HF ADC1_9 ADC1_10 ADC1_11
+ X= T_MUGAST_HF EVAMOS_AGATA
+ Y= TVAMOS_MUGAST_HF E_AGATA
diff --git a/Projects/e775s/analysis/InputFiles/runPhaseSpace.mac b/Projects/e775s/analysis/InputFiles/runPhaseSpace.mac
new file mode 100644
index 0000000000000000000000000000000000000000..d9502085bfe4544342ddc937bcdbca74211c96f2
--- /dev/null
+++ b/Projects/e775s/analysis/InputFiles/runPhaseSpace.mac
@@ -0,0 +1 @@
+/run/beamOn 10000000
diff --git a/Projects/e775s/analysis/InputFiles/runSimulation.mac b/Projects/e775s/analysis/InputFiles/runSimulation.mac
new file mode 100644
index 0000000000000000000000000000000000000000..325d88eef6b6df38e8ffac1d8fc402ccc816dc69
--- /dev/null
+++ b/Projects/e775s/analysis/InputFiles/runSimulation.mac
@@ -0,0 +1 @@
+/run/beamOn 1000000
diff --git a/Projects/e775s/analysis/Plots_C2SFunctions.h b/Projects/e775s/analysis/Plots_C2SFunctions.h
new file mode 100644
index 0000000000000000000000000000000000000000..16e4efa1e6d06e80069e01fe5549fb5fcc8141f4
--- /dev/null
+++ b/Projects/e775s/analysis/Plots_C2SFunctions.h
@@ -0,0 +1,1239 @@
+/* Predefine functions */
+void C2S();
+void C2S_Diagnosis();
+vector<vector<double>> GetExpDiffCross(double Energy);
+TH1F* PullThetaLabHist(int i, double minTheta, double gatesize);
+TH1F* PullThetaCMHist(int i, double minTheta, double gatesize);
+TH1F* PullPhaseSpaceHist(int i, double minTheta, double gatesize);
+void Scale(TGraph* g , TGraphErrors* ex);
+TGraph* TWOFNR(double E, double J0, double J, double n, double l, double j);
+double ToMininize(const double* parameter);
+TGraph* FindNormalisation(TGraph* theory, TGraphErrors* experiment);
+TGraph* FindNormalisation(TGraph* theory, TGraph* theory2, TGraphErrors* experiment);
+TList* peakFitList = new TList();
+
+/* Global variables */
+vector<Double_t> anglecentres, anglewidth;
+TGraph* currentThry;
+TGraph* currentThry2;
+TGraphErrors* staticExp;
+int indexE;
+double globalS, globalSerr;
+double globalS2, globalSerr2;
+int numAngleBins;
+double widthAngleBins;
+double firstAngle;
+double CSangleL, CSangleH;
+double CSvertL, CSvertH;
+bool GammaGate = false;
+double globGammaGate;
+bool doMix = false;
+bool doDoublet = false;
+double globalBinning;
+double SnTWOFNR = 7.607;
+
+/* Output volume toggle */
+bool loud = 1;
+
+/* Scale method toggle */
+bool scaleTogether = 1;
+
+/* Strings for image */
+string orbitalname;
+string orbital;
+
+/* Strings for SolidAngle input file */
+string statename;
+string inputdate;
+
+
+////////////////////////////////////////////////////////////////////////////////
+void WriteToFile(TGraphErrors* gr, string filename){
+ int npoints = gr->GetN();
+ ofstream file;
+ file.open (filename.c_str());
+
+ cout << " ============================ " << endl;
+ cout << " Writing to " << filename << " ..." << endl;
+
+ file << " i \t x \t y \t xe \t ye" << endl;
+
+ for(int i = 0; i < npoints; i++){
+ double x = gr->GetPointX(i);
+ double y = gr->GetPointY(i);
+ double xe = gr->GetErrorX(i);
+ double ye = gr->GetErrorY(i);
+
+ file << i << "\t"
+ << x << "\t"
+ << y << "\t"
+ << xe << "\t"
+ << ye << "\t"
+ << endl;
+
+ }
+
+ cout << " Complete! " << endl;
+ cout << " ============================ " << endl;
+
+}
+
+////////////////////////////////////////////////////////////////////////////////
+
+TGraphErrors* GetExpDCS_ForFig(string filename){
+ vector<double> x, dx, y, dy;
+ ifstream f(filename.c_str());
+ string line;
+ while (getline(f, line)){
+ stringstream ss(line);
+ // i x y dx dy
+ double tA, tB, tC, tD, tE;
+ if (f >> tA >> tB >> tC >> tD >> tE){
+ //cout << tA << " " << tB << " " << tC << " " << tD << " " << tE << endl;
+ x.push_back(tB);
+ y.push_back(tC);
+ dx.push_back(0.0);
+ dy.push_back(tE);
+ }
+ }
+ f.close();
+
+ TGraphErrors* g = new TGraphErrors(x.size(), &(x[0]), &(y[0]), &(dx[0]), &(dy[0]));
+ return g;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+TGraph* GetSimDCS_ForFig(string filename){
+ vector<double> x, dx, y, dy;
+ ifstream f(filename.c_str());
+ string line;
+ while (getline(f, line)){
+ stringstream ss(line);
+ // i x y dx dy
+ double tA, tB, tC, tD, tE;
+ if (f >> tA >> tB >> tC >> tD >> tE){
+ //cout << tA << " " << tB << " " << tC << " " << tD << " " << tE << endl;
+ x.push_back(tB);
+ y.push_back(tC);
+ dx.push_back(tD);
+ dy.push_back(tE);
+ }
+ }
+ f.close();
+
+ TGraph* g = new TGraph(x.size(), &(x[0]), &(y[0]));
+ return g;
+}
+
+void CanvasPartition(TCanvas *C,const Int_t Nx = 2,const Int_t Ny = 2,
+ Float_t lMargin = 0.15, Float_t rMargin = 0.05,
+ Float_t bMargin = 0.15, Float_t tMargin = 0.05);
+double XtoPad(double x);
+double YtoPad(double x);
+
+////////////////////////////////////////////////////////////////////////////////
+void CanvasPartition(TCanvas *C,const Int_t Nx,const Int_t Ny, Float_t lMargin, Float_t rMargin, Float_t bMargin, Float_t tMargin){
+ if (!C) return;
+
+ // Setup Pad layout:
+ Float_t vSpacing = 0.0;
+ Float_t vStep = (1.- bMargin - tMargin - (Ny-1) * vSpacing) / Ny;
+
+ Float_t hSpacing = 0.0;
+ Float_t hStep = (1.- lMargin - rMargin - (Nx-1) * hSpacing) / Nx;
+
+ Float_t vposd,vposu,vmard,vmaru,vfactor;
+ Float_t hposl,hposr,hmarl,hmarr,hfactor;
+
+ for (Int_t i=0;i<Nx;i++) {
+
+ if (i==0) {
+ hposl = 0.0;
+ hposr = lMargin + hStep;
+ hfactor = hposr-hposl;
+ hmarl = lMargin / hfactor;
+ hmarr = 0.0;
+ } else if (i == Nx-1) {
+ hposl = hposr + hSpacing;
+ hposr = hposl + hStep + rMargin;
+ hfactor = hposr-hposl;
+ hmarl = 0.0;
+ hmarr = rMargin / (hposr-hposl);
+ } else {
+ hposl = hposr + hSpacing;
+ hposr = hposl + hStep;
+ hfactor = hposr-hposl;
+ hmarl = 0.0;
+ hmarr = 0.0;
+ }
+
+ for (Int_t j=0;j<Ny;j++) {
+
+ if (j==0) {
+ vposd = 0.0;
+ vposu = bMargin + vStep;
+ vfactor = vposu-vposd;
+ vmard = bMargin / vfactor;
+ vmaru = 0.0;
+ } else if (j == Ny-1) {
+ vposd = vposu + vSpacing;
+ vposu = vposd + vStep + tMargin;
+ vfactor = vposu-vposd;
+ vmard = 0.0;
+ vmaru = tMargin / (vposu-vposd);
+ } else {
+ vposd = vposu + vSpacing;
+ vposu = vposd + vStep;
+ vfactor = vposu-vposd;
+ vmard = 0.0;
+ vmaru = 0.0;
+ }
+
+ C->cd(0);
+
+ auto name = TString::Format("pad_%d_%d",i,j);
+ auto pad = (TPad*) C->FindObject(name.Data());
+ if (pad) delete pad;
+ pad = new TPad(name.Data(),"",hposl,vposd,hposr,vposu);
+ pad->SetLeftMargin(hmarl);
+ pad->SetRightMargin(hmarr);
+ pad->SetBottomMargin(vmard);
+ pad->SetTopMargin(vmaru);
+
+ pad->SetFrameBorderMode(0);
+ pad->SetBorderMode(0);
+ pad->SetBorderSize(0);
+
+ pad->Draw();
+ }
+ }
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void WriteToFile(TGraph* gr, string filename){
+ int npoints = gr->GetN();
+ ofstream file;
+ file.open (filename.c_str());
+
+ cout << " ============================ " << endl;
+ cout << " Writing to " << filename << " ..." << endl;
+
+ file << " i \t x \t y \t xe \t ye" << endl;
+
+ for(int i = 0; i < npoints; i++){
+ double x = gr->GetPointX(i);
+ double y = gr->GetPointY(i);
+ double xe = gr->GetErrorX(i);
+ double ye = gr->GetErrorY(i);
+
+ file << i << "\t"
+ << x << "\t"
+ << y << "\t"
+ << xe << "\t"
+ << ye << "\t"
+ << endl;
+
+ }
+
+ cout << " Complete! " << endl;
+ cout << " ============================ " << endl;
+
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void canclone(TCanvas* major, int padNum, string name){
+ string minorName = "./CS2_Figures/" + name + ".root";
+ TFile* minorFile = new TFile(minorName.c_str(),"READ");
+ TCanvas* minor = (TCanvas*) minorFile->FindObjectAny("c_peakFits");
+
+ major->cd(padNum);
+ TPad *pad = (TPad*)gPad;
+ minor->cd();
+ TObject *obj, *clone;
+ pad->Range(minor->GetX1(),minor->GetY1(),minor->GetX2(),minor->GetY2());
+ pad->SetTickx(minor->GetTickx());
+ pad->SetTicky(minor->GetTicky());
+ pad->SetGridx(minor->GetGridx());
+ pad->SetGridy(minor->GetGridy());
+ pad->SetLogx(minor->GetLogx());
+ pad->SetLogy(minor->GetLogy());
+ pad->SetLogz(minor->GetLogz());
+ pad->SetBorderSize(minor->GetBorderSize());
+ pad->SetBorderMode(minor->GetBorderMode());
+ minor->TAttLine::Copy((TAttLine&)*pad);
+ minor->TAttFill::Copy((TAttFill&)*pad);
+ minor->TAttPad::Copy((TAttPad&)*pad);
+
+ TIter next(minor->GetListOfPrimitives());
+ gROOT->SetSelectedPad(pad);
+ while ((obj=next())) {
+ clone = obj->Clone();
+ obj->Draw("SAME");
+ pad->GetListOfPrimitives()->Add(clone,obj->GetDrawOption());
+ }
+ pad->Modified();
+ pad->Update();
+ major->cd(padNum);
+ pad->Draw();
+}
+
+////////////////////////////////////////////////////////////////////////////////
+double GetNodes(double spdf){
+ double nodes;
+ if(spdf==0){ nodes=1; }
+ else if(spdf==1){ nodes=0; }
+ else if(spdf==2){ nodes=0; }
+ else if(spdf==3){ nodes=0; }
+ else{
+ cout << " INPUT NODES::" << endl;
+ cin >> nodes;
+ }
+ return nodes;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void GenerateResidual(TGraphErrors* gdSdO, TGraph* Final){
+
+ TGraphErrors* gResid = new TGraphErrors(*gdSdO);
+ for(int n=0; n < gResid->GetN(); n++){
+ double x = gdSdO->GetPointX(n);
+ double residual = gdSdO->GetPointY(n) - Final->Eval(x);
+ gResid->SetPoint(n,x,residual);
+ gResid->SetPointError(n,0,gdSdO->GetErrorY(n));
+ }
+ TLine* markzero = new TLine(CSangleL,0.,CSangleH,0.);
+ gResid->SetTitle("");
+ gResid->GetXaxis()->SetRangeUser(CSangleL,CSangleH);
+ gResid->GetXaxis()->SetNdivisions(512,kTRUE);
+ gResid->GetYaxis()->SetTitle("Residuals");
+ gResid->GetYaxis()->SetTitleSize(0.15);
+ gResid->GetYaxis()->SetTitleOffset(0.36);
+ gResid->GetYaxis()->SetLabelSize(0.08);
+ gResid->GetYaxis()->SetNdivisions(305);
+ gResid->GetXaxis()->SetTitleSize(0.15);
+ gResid->GetXaxis()->SetTitleOffset(0.8);
+ gResid->GetXaxis()->SetLabelSize(0.1);
+ gResid->GetXaxis()->SetTickLength(0.1);
+ gResid->Draw();
+ markzero->SetLineStyle(2);
+ markzero->Draw("same");
+
+}
+
+////////////////////////////////////////////////////////////////////////////////
+TGraph* AddTGraphs(TGraph* a, TGraph* b){
+ if(a->GetN() != b->GetN()){
+ cout << "CANNOT ADD GRAPHS, UNEQUAL NUMBER OF POINTS" << endl;
+ cout << a->GetN() << " != " << b->GetN() << endl;
+ return 0;
+ }
+
+ TGraph* sum = (TGraph*) a->Clone();
+
+ int maxN = a->GetN();
+ for(int i = 0 ; i < maxN ; i++){
+ double xa = a->GetPointX(i);
+ double xb = b->GetPointX(i);
+ if(xa!=xb){cout << "FAIL! Unequal X!" << endl; return 0;}
+
+ double ya = a->GetPointY(i);
+ double yb = b->GetPointY(i);
+ sum->SetPoint(i,xa,ya+yb);
+ }
+
+ return sum;
+
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void C2S(double Energy, double Spin, double spdf, double angmom, double binning, string option){
+ /* Clean global variables, in case of multiple runs */
+ indexE = 100;
+ anglecentres.clear();
+ anglewidth.clear();
+ globalS=0.;
+ globalSerr=0.;
+ globalS2=0.;
+ globalSerr2=0.;
+ peakFitList->Clear();
+ globalBinning = binning;
+
+ //gROOT->SetBatch(1);
+
+ double spdf2, angmom2 = 0;
+ doMix = false;
+ doDoublet = false;
+
+ if(option=="mixed"){
+ doMix = true;
+ cout << BOLDBLUE
+ << " || MIXING OPTION SELECTED || \n"
+ << " || GIVE DETAILS OF 2ND OPT || \n"
+ << " || L = ... || \n"
+ << endl;
+ cin >> spdf2;
+ cout << " || J = ... || \n"
+ << endl;
+ cin >> angmom2;
+ cout << RESET;
+ } else if (option=="doublet") {
+ doDoublet = true;
+ }
+
+ if(binning < 0.01){
+ cout << BOLDRED
+ << " BINNING IS BELOW 0.01 MeV! Exiting..."
+ << RESET
+ << endl;
+ return;
+ }
+
+ /* Assign local variables */
+ // s1/2 -> spdf = 0, angmom = 0.5
+ // J0 is incident spin, which is 19O g.s. therefore J0 = 5/2
+ double J0 = 2.5;
+ double ElasticNorm = 1.0, ElasticNormErr = 0.0;
+
+ string backupFileName = "SolidAngleHistFiles/";
+ string saFileName = "SolidAngleHistFiles/SAHF_Sim_19Odp_";
+
+ // Define angle variables
+//numAngleBins=26; widthAngleBins=2.0; firstAngle=104.0;
+//numAngleBins=25; widthAngleBins=2.0; firstAngle=106.0;
+ numAngleBins=19; widthAngleBins=2.0; firstAngle=106.0;
+ CSangleL = 100.;
+ CSangleH = 150.;
+ CSvertL = 0.01; //CSvertL = 2e-6;
+ CSvertH = 100.; //CSvertH = 5e-3;
+
+ // Reduce angular range for high energy states
+ if(Energy>7.2){
+ double tmax = (Energy-18.)/(-0.075);
+ numAngleBins = (int)floor((double)(floor(tmax)-firstAngle)/widthAngleBins);
+ cout << "!! TRIMMED NUMBER OF ANGLE BINS" << numAngleBins << endl;
+ cout << "!! THEREFORE ANGE RANGE = "
+ << firstAngle << " to "
+ << firstAngle+(widthAngleBins*numAngleBins)
+ << endl;
+ }
+
+ // Extract spin orbit name
+ orbitalname.clear();
+ orbital.clear();
+
+ switch((int)spdf){
+ case 0:
+ orbitalname.append("s_{");
+ orbital.append("s");
+ break;
+ case 1:
+ orbitalname.append("p_{");
+ orbital.append("p");
+ break;
+ case 2:
+ orbitalname.append("d_{");
+ orbital.append("d");
+ break;
+ case 3:
+ orbitalname.append("f_{");
+ orbital.append("f");
+ break;
+ default:
+ cout << RED
+ <<"ERROR: SPDF SELECTION OUT OF RANGE"
+ << endl;
+ }
+ orbitalname.append(to_string((int)(2*angmom)));
+ orbitalname.append("/2}");
+ orbital.append(to_string((int)(2*angmom)));
+ orbital.append("2");
+
+ // Number of nodes
+ double nodes = GetNodes(spdf);
+ double nodes2;
+ if(doMix){
+ nodes2 = GetNodes(spdf2);
+ }
+
+ /* Retrieve array index of the entered peak energy */
+ /* numpeaks and Energy[] defined globally in KnownPeakFitter.h */
+ bool found = 0;
+ for(int i=0;i<numPeaks;i++){
+ if(abs(Energy-means.at(i))<0.01){
+ cout << "========================================================" << endl;
+ cout << "Identified as state #" << i << ", E = " << means.at(i) << endl;
+ indexE = i;
+ found = 1;
+ stringstream ss;
+ ss << setfill('0') << setw(4) << (int)((means.at(i)+0.0001)*1000.);
+ statename = ss.str();
+ }
+ }
+ if(!found){
+ cout << "========================================================" << endl;
+ cout << "NO STATE AT THAT ENERGY INDENTIFIED!! CHECK KNOWN PEAKS!!" << endl;
+ return;
+ }
+ saFileName.append(statename);
+ saFileName.append(".root");
+
+ /* Solid Angle (from simulation) */
+ TFile* file;
+ ifstream infile(saFileName);
+ ifstream backup(backupFileName);
+ if(infile.good()){
+ cout << GREEN << "Opening file " << saFileName << RESET << endl;
+ file = TFile::Open(saFileName.c_str());
+ } else if (backup.good()){
+ cout << BOLDRED << "FAILED TO OPEN " << saFileName << endl;
+ cout << "Open standard backup file " << backupFileName << endl;
+ cout << RESET;
+ file = TFile::Open(backupFileName.c_str());
+ } else {
+ cout << BOLDRED << "FAILED TO OPEN MAIN OR BACKUP SOLID ANGLE FILE" << endl;
+ cout << RED << "Check SolidAngle file exists..." << RESET << endl;
+ }
+
+ string objName = "SolidAngle_Lab_MG";
+ TH1F* SolidAngle = (TH1F*) file->FindObjectAny(objName.c_str());
+ TCanvas* c_SolidAngle = new TCanvas("c_SolidAngle","c_SolidAngle",1000,1000);
+ SolidAngle->Draw("HIST");
+ SolidAngle->GetXaxis()->SetRangeUser(CSangleL,CSangleH);
+ /* (canvas deleted after Area/SA calculation) */
+
+ /* Area of experimental peaks */
+ TCanvas* c_PeakArea = new TCanvas("c_PeakArea","c_PeakArea",1000,1000);
+ vector<vector<double>> areaArray = GetExpDiffCross(means.at(indexE));
+// delete c_PeakArea;
+
+ // Array: peakenergy, peakarea, areaerror, anglemin, anglemax
+ if(loud){
+ for(int i=0; i<areaArray.size();i++){
+ cout << i << " "
+ << areaArray[i][0] << " "
+ << areaArray[i][1] << " "
+ << areaArray[i][2] << " "
+ << areaArray[i][3] << " "
+ << areaArray[i][4] << endl;
+ }
+ }
+
+ /* AoSA = Area/Solid Angle [#/msr] */
+ /* dSdO = Experimental Diff. Cross Sect. (Area/Solid Angle *Norm) [mb/msr] */
+ vector<Double_t> AoSA, AoSAerr;
+ vector<Double_t> expDCS, expDCSerr;
+ for(int i=0; i<areaArray.size();i++){
+ int binmin = SolidAngle->FindBin(areaArray[i][3]+0.0001);
+ int binmax = SolidAngle->FindBin(areaArray[i][4]-0.0001);
+
+ anglecentres.push_back(((areaArray[i][4]-areaArray[i][3])/2.)+areaArray[i][3]);
+ anglewidth.push_back(areaArray[i][4]-areaArray[i][3]);
+
+ double tempsum=0, tempsumerr=0;
+ for(int x=binmin;x<binmax+1;x++){
+ tempsum += SolidAngle->GetBinContent(x);
+ tempsumerr += SolidAngle->GetBinError(x);
+ if(loud){cout << x << endl;}
+ }
+ if(loud){cout << "TEST CHECK " << tempsum << " +- " << tempsumerr << endl;}
+
+ double SAerr;
+ double SA = SolidAngle->IntegralAndError(binmin,binmax,SAerr);
+ SA = SA*1000.; //sr->msr
+ SAerr = SAerr*1000.; //sr->msr
+
+ /* Area over Solid angle ONLY */
+ AoSA.push_back(areaArray[i][1]/SA);
+ AoSAerr.push_back((areaArray[i][1]/SA)
+ * sqrt( pow(areaArray[i][2]/areaArray[i][1],2)
+ + pow(SAerr/SA,2)));
+
+ /* Differential Cross Section */
+ /* NOTE: DON'T INCLUDE NORM ERROR IN ERR PROPAGATION AS IT IS SYSTEMATIC! */
+ double tempvalue = (areaArray[i][1]/SA)*ElasticNorm;
+ double temperror = tempvalue
+ * sqrt( pow(areaArray[i][2]/areaArray[i][1],2)
+ + pow(SAerr/SA,2)
+ );
+ //if(temperror<tempvalue){ // exclude very large error bars
+ expDCS.push_back(tempvalue);
+ expDCSerr.push_back(temperror);
+ //}
+
+ if(loud){
+ cout << "Angle " << areaArray[i][3] << " to " << areaArray[i][4]
+ << ", centre " << anglecentres[i]
+ << ": Area = " << areaArray[i][1] << " +- " << areaArray[i][2] << " cnts"
+ << " SA = " << SA << " +- " << SAerr << " msr"
+ << " Area/SA = " << AoSA[i] << " +- " << AoSAerr[i] << " counts/msr"
+ << setprecision(5)
+ << " Norm = " << ElasticNorm << " +- " << ElasticNormErr
+ << " (don't include norm err in propagation)"
+ << " dSdO = " << tempvalue << " +- " << temperror
+ << setprecision(3)
+ << endl;
+ }
+
+ }
+ //delete c_SolidAngle;
+
+ gROOT->SetBatch(0);
+
+ /* Graph of Area/Solid Angle*/
+ TCanvas* c_AoSA = new TCanvas("c_AoSA","c_AoSA",1000,1000);
+ c_AoSA->SetLogy();
+ TGraphErrors* gAoSA = new TGraphErrors(
+ anglecentres.size(), //n
+ &(anglecentres[0]), &(AoSA[0]), //x, y
+ 0, &(AoSAerr[0]) ); //errX, errY
+ gAoSA->SetTitle("Area/SolidAngle");
+ gAoSA->GetXaxis()->SetTitle("ThetaLab [deg]");
+ gAoSA->GetYaxis()->SetTitle("Counts/#Omega [counts/msr]");
+ gAoSA->Draw();
+
+ /* Graph of experimental diff. cross sect (dSigma/dOmega) */
+ TCanvas* c_dSdO = new TCanvas("c_dSdO","c_dSdO",1000,1000);
+ c_dSdO->SetLogy();
+ TGraphErrors* gdSdO = new TGraphErrors(
+ anglecentres.size(),
+ &(anglecentres[0]), &(expDCS[0]),
+ 0, &(expDCSerr[0]) );
+ gdSdO->SetTitle("Differential Cross Section");
+ gdSdO->GetXaxis()->SetTitle("#theta_{lab} [deg]");
+ gdSdO->GetXaxis()->SetTitleOffset(1.2);
+ gdSdO->GetXaxis()->SetTitleSize(0.04);
+ gdSdO->GetYaxis()->SetTitle("d#sigma/d#Omega [mb/msr]");
+ gdSdO->GetYaxis()->SetTitleOffset(1.2);
+ gdSdO->GetYaxis()->SetTitleSize(0.04);
+ gdSdO->Draw();
+ c_dSdO->Update();
+
+ /* TWOFNR diff. cross section, in mb/msr */
+ double meanTWOFNR = means.at(indexE);
+ if(meanTWOFNR > SnTWOFNR){meanTWOFNR = SnTWOFNR-0.001;}
+ TCanvas* c_TWOFNR = new TCanvas("c_TWOFNR","c_TWOFNR",1000,1000);
+ c_TWOFNR->SetLogy();
+ TGraph* TheoryDiffCross = TWOFNR(meanTWOFNR, Spin, nodes, spdf, angmom);
+ TheoryDiffCross->GetYaxis()->SetTitle("d#sigma/d#Omega [mb/msr]"); //msr set in func above
+ TheoryDiffCross->GetXaxis()->SetTitle("ThetaLab [deg]");
+ TheoryDiffCross->Draw();
+
+ TGraph* TheoryDiffCross2;
+ if(doMix){
+ TheoryDiffCross2 = TWOFNR(meanTWOFNR, Spin, nodes2, spdf2, angmom2);
+ TheoryDiffCross2->SetLineColor(kBlue);
+ TheoryDiffCross2->Draw("same");
+ }
+
+ /** TEMP **/
+ if(loud){
+ cout << "UNSCALED THEORY DIFF CROSS SECTION EVALUATED AT DATA POINTS:::" << endl;
+ cout << setprecision(6);
+ for(int i=0; i<numAngleBins; i++){
+ cout << anglecentres.at(i) << "\t" << TheoryDiffCross->Eval(anglecentres.at(i)) << endl;
+ }
+ cout << setprecision(3);
+ cout << "......................" << endl;
+ }
+ /** **** **/
+
+ /* Using Chi2 minimizaiton */
+ if(loud){cout << "USING CHI2 MINIMIZAITON..." << endl;}
+ TCanvas* c_Chi2Min = new TCanvas("c_Chi2Min","c_Chi2Min",1000,1000);
+ gStyle->SetPadLeftMargin(0.12);
+ gStyle->SetPadRightMargin(0.03);
+
+ TPad *pad1 = new TPad("pad1","pad1",0,0.25,1,1);
+ TPad *pad2 = new TPad("pad2","pad2",0,0,1,0.25);
+ pad1->SetTopMargin(0.1);
+ pad1->SetBottomMargin(0.00001);
+ pad1->SetBorderMode(0);
+ pad1->SetLogy();
+ pad1->SetTickx();
+ pad1->SetTicky();
+ pad2->SetTopMargin(0.00001);
+ pad2->SetBottomMargin(0.3);
+ pad2->SetBorderMode(0);
+ pad2->SetTickx();
+ pad2->SetTicky();
+ pad1->Draw();
+ pad2->Draw();
+ pad1->cd();
+
+ TGraphErrors* gdSdO2 = new TGraphErrors(*gdSdO);
+ gdSdO2->SetName("gdSdO2");
+ gdSdO2->SetTitle("gdSdO2");
+
+ TGraph* Final;
+
+ if(doMix){
+ Final = FindNormalisation(TheoryDiffCross,TheoryDiffCross2,gdSdO);
+ } else {
+ Final = FindNormalisation(TheoryDiffCross,gdSdO);
+ }
+
+ Final->SetName("Final");
+ Final->SetTitle("Final");
+
+ gdSdO->SetLineColor(kRed);
+ gdSdO->SetMarkerColor(kRed);
+ gdSdO->SetMarkerStyle(21);
+
+
+ /* Construct file name string */
+ /**/ ostringstream tempstream;
+ /**/ if(means.at(indexE)<1.0){tempstream << 0;}
+ /**/ tempstream << (int) (means.at(indexE)*1000);
+ /**/ tempstream << "_" << orbital;
+ /**/ tempstream << "_spin" << Spin;
+ /**/ string tempstr = tempstream.str();
+ /* Construct hist title string */
+ /**/ ostringstream textstream;
+ /**/ textstream << std::fixed << setprecision(3);
+ /**/ textstream << " " << means.at(indexE);
+ /**/ textstream << " MeV, ";
+ /**/ textstream << orbitalname;
+ /**/ textstream << ", spin " << (int)Spin;
+ /**/ textstream << " --> scale (not norm.) = " << globalS
+ /**/ << " +- " << globalSerr;
+ /**/ if(doMix){
+ /**/ textstream << " & " << globalS2
+ /**/ << " +- " << globalSerr2;
+ /**/ }
+ /**/ string textstring = textstream.str();
+
+ gdSdO->SetTitle(textstring.c_str());
+ gdSdO->GetYaxis()->SetTitleOffset(1.3);
+ gdSdO->GetYaxis()->SetTitleSize(0.042);
+ gdSdO->GetXaxis()->SetRangeUser(CSangleL,CSangleH);
+ //gdSdO->GetYaxis()->SetRangeUser(5e-6,5e-3); //Same vertical range for all? WNC suggestion
+ gdSdO->GetYaxis()->SetRangeUser(CSvertL,CSvertH); //Same vertical range for all? WNC suggestion
+ gdSdO->GetXaxis()->SetNdivisions(512,kTRUE);
+ gdSdO->Draw("AP");
+
+ Final->Draw("SAME");
+
+ pad2->cd();
+
+ GenerateResidual(gdSdO,Final);
+
+ string savestring1 = "./C2S_Outputs/"+tempstr+".root";
+ string savestring2 = "./C2S_Outputs/"+tempstr+".pdf";
+ c_Chi2Min->SaveAs(savestring1.c_str());
+ c_Chi2Min->SaveAs(savestring2.c_str());
+
+ cout << YELLOW
+ << " ----- C2S = (2J+1)S = (2*" << (int)Spin << "+1)S = ("
+ << (int)((2.*Spin)+1.) << ")S = "
+ << ((2.*Spin)+1.) * globalS << " ----- "
+ << RESET << endl;
+
+
+ //delete c_AoSA;
+ //delete c_dSdO;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+vector<vector<double>> GetExpDiffCross(double Energy){
+ cout << "========================================================" << endl;
+ vector<vector<double>> AllPeaks_OneGate;
+ vector<vector<double>> OnePeak_AllGates;
+ /****CHANGE ANGLE GATING****/
+ //double widthAngleBins = 2.5;
+ //double firstAngle = 105.;
+ /***************************/
+ double x[numAngleBins], y[numAngleBins];
+ //TList* list = new TList();
+
+ double trackScale = 0.0;
+ /* Determine scaling factor for PhaseSpace */
+ TCanvas* c_ExSubPSpace = new TCanvas("c_ExSubPSpace","c_ExSubPSpace",1000,1000);
+ //if(scaleTogether){
+ // cout << BLUE << endl;
+ // cout << " SUBTRACTING PHASE SPACE TOGETHER! " << endl;
+ // TH1F* baseEx = PullThetaLabHist(0,firstAngle,widthAngleBins);
+ // TH1F* basePS = PullPhaseSpaceHist(0,firstAngle,widthAngleBins);
+ // for(int i=1; i<numAngleBins;i++){
+ // TH1F* addEx = PullThetaLabHist(i,firstAngle,widthAngleBins); baseEx->Add(addEx,1.);
+ // TH1F* addPS = PullPhaseSpaceHist(i,firstAngle,widthAngleBins); basePS->Add(addPS,1.);
+ // }
+ //
+ //// /* Subtract flat background equal to smallest bin in range */
+ //// baseEx->GetXaxis()->SetRange(baseEx->FindBin(-1.),baseEx->FindBin(1.));
+ //// double minValueInRange = baseEx->GetBinContent(baseEx->GetMinimumBin());
+ //// baseEx->GetXaxis()->UnZoom();
+ //// cout << "Subtracting background of " << minValueInRange << endl;
+ //// for(int b=1; b<baseEx->GetNbinsX() ; b++){
+ //// baseEx->SetBinContent(b,baseEx->GetBinContent(b)-minValueInRange);
+ //// }
+ //
+ // // Fix error where basePS = 600, baseEx=150
+ // if(loud){
+ // cout << "baseEx bins: " << baseEx->GetNbinsX() << endl;
+ // cout << "basePS bins: " << basePS->GetNbinsX() << endl;
+ // }
+ // if(basePS->GetNbinsX()!=baseEx->GetNbinsX()){
+ // cout << basePS->GetNbinsX() / baseEx->GetNbinsX() << endl;
+ // int ratio = basePS->GetNbinsX() / baseEx->GetNbinsX();
+ // basePS->Rebin(ratio);
+ // }
+ // if(loud){
+ // cout << "baseEx bins: " << baseEx->GetNbinsX() << endl;
+ // cout << "basePS bins: " << basePS->GetNbinsX() << endl;
+ // }
+
+ // /* Begin scaling within range, track changes */
+ // trackScale = 0.7;
+ // basePS->Scale(trackScale);
+ // int numAngleBinsScale = baseEx->GetNbinsX();
+ // int nbinlow = basePS->FindBin(7.5); int nbinhigh = basePS->FindBin(10.0);
+ // for(int b=nbinlow; b<nbinhigh; b++){
+////cout << "----" << endl;
+ // if(baseEx->GetBinContent(b) > 0.0 && basePS->GetBinContent(b) > baseEx->GetBinContent(b)){
+ // while(basePS->GetBinContent(b) > baseEx->GetBinContent(b)){
+ // basePS->Scale(0.99999);
+ // trackScale *= 0.99999;
+////cout << "Ex: " << baseEx->GetBinContent(b) << " PS: " << basePS->GetBinContent(b) << endl;
+ // }
+ // }
+ // }
+
+ // TH1F* finalEx = (TH1F*) baseEx->Clone();
+ // finalEx->SetName("ExSubPSpace");
+ // finalEx->SetTitle("ExSubPSpace");
+ // finalEx->Add(basePS,-1.);
+ // baseEx->SetLineColor(kBlack); baseEx->Draw("hist");
+ // basePS->SetLineColor(kOrange); basePS->Draw("hist same");
+ // finalEx->SetLineColor(kRed); finalEx->Draw("hist same");
+ // if(loud){cout << "PhaseSpace -> ExpData scaling = " << trackScale << endl;}
+ // cout << RESET << endl;
+ //}
+
+ /* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
+ //for(int i=0; i<1; i++){
+ for(int i=0; i<numAngleBins;i++){
+ double min = firstAngle + (i*widthAngleBins);
+ double max = min + widthAngleBins;
+ cout << GREEN << endl;
+ cout << "===================================" << endl;
+ cout << " RUNNING SUB AND FIT " << endl;
+ cout << "min: " << min << " max: " << max << endl;
+
+ stringstream tmp; tmp << fixed << setprecision(0);
+ tmp << "c_peakFits_" << min << "_" << max;
+ string tmp2 = tmp.str();
+ TCanvas* c_peakFits = new TCanvas("c_peakFits",tmp2.c_str(),1000,1000);
+
+ /* Retrieve theta-gated Ex TH1F from file GateThetaLabHistograms.root */
+ /* To change angle gates, run GateThetaLab_MultiWrite() */
+ TH1F *gate, *pspace;
+ gate = PullThetaLabHist(i,firstAngle,widthAngleBins);
+ pspace = PullPhaseSpaceHist(i,firstAngle,widthAngleBins);
+
+ /* Scale the Phase Space at this angle... */
+ /* ... using DrawExThetaLab_SubPSp value */
+ gate->Add(pspace,-0.134256);
+ for(int b=0; b<gate->GetNbinsX();b++){
+ if(gate->GetBinContent(b)<0){
+ if(loud){cout << "flattening bin " << gate->GetBinCenter(b) << ", from " << gate->GetBinContent(b) << " to 0" << endl;}
+ gate->SetBinContent(b,0);
+ }
+ }
+
+ /* ... for all angles together */
+ //if(scaleTogether){
+ // if(loud){ cout << "gate bins: " << gate->GetNbinsX() << endl;
+ // cout << "pspc bins: " << pspace->GetNbinsX() << endl; }
+ // // Fix error where basePS = 600, baseEx=150
+ // if(pspace->GetNbinsX()!=gate->GetNbinsX()){
+ // if(loud){cout << pspace->GetNbinsX() / gate->GetNbinsX() << endl;}
+ // int ratio = pspace->GetNbinsX() / gate->GetNbinsX();
+ // pspace->Rebin(ratio);
+ // }
+ // if(loud){ cout << "gate bins: " << gate->GetNbinsX() << endl;
+ // cout << "pspc bins: " << pspace->GetNbinsX() << endl; }
+ // cout << "SUMMATION: Subtracting phase space, scaled by " << trackScale << endl;
+ // gate->Add(pspace,-trackScale);
+ //}
+ //
+ ///* ... or seperately for each angular bin */
+ ///* NOTE THAT THIS DOES NOT ACCOUNT FOR FLAT BACKGROUND */
+ //else {
+ // if(pspace->Integral() > 50.){ // Non-garbage histogram
+ // trackScale=0.7;
+ // pspace->Scale(trackScale);
+ // int numAngleBins = gate->GetNbinsX();
+ // for(int b=0; b<numAngleBins; b++){
+ // if(loud){cout << " FROM " << pspace->GetBinContent(b) <<
+ // " > " << gate->GetBinContent(b);
+ // }
+ // while(pspace->GetBinContent(b) > gate->GetBinContent(b)){
+ // pspace->Scale(0.9999);
+ // trackScale*=0.9999;
+ // }
+ // if(loud){cout << " TO " << pspace->GetBinContent(b) <<
+ // " > " << gate->GetBinContent(b) << endl;
+ // }
+ // }
+ // if(loud){cout << " !!! SCALE FOR THIS ANGLE = " << trackScale << endl;}
+ // cout << "INDIVIDUAL: Subtracting phase space, scaled by " << trackScale << endl;
+ // gate->Add(pspace,-1);
+ // }
+ //}
+ // Retrieve array containing all fits, for one angle gate. //
+ // Specific peak of interest selected from the vector by //
+ // global variable indexE //
+ gate->Draw();
+ AllPeaks_OneGate = FitKnownPeaks_RtrnArray(gate);
+
+
+ /* Write PS-subtracted spectrum to list */
+ //list->Add(gate);
+ //list->Add(c_peakFits);
+ gate->GetXaxis()->SetRangeUser(-1.0,+8.0);
+ string filename = "./C2S_Outputs/";
+ filename += tmp2 + ".root";
+ c_peakFits->SaveAs(filename.c_str());
+ //auto tempfile = new TFile(filename.c_str(),"UPDATE"); //Reopen newly made file
+ //auto templist = new TList();
+ //templist->Add();
+
+
+ /* Check correct OneGate vector is selected */
+ if(loud){
+ cout << "area of " << means.at(indexE) << " = "
+ << AllPeaks_OneGate[indexE][1]
+ << " +- " << AllPeaks_OneGate[indexE][2]
+ << endl;
+ }
+
+ /* Add min and max angle to end of relevant OneGate vector */
+ AllPeaks_OneGate[indexE].push_back(min);
+ AllPeaks_OneGate[indexE].push_back(max);
+
+ /* Push relevant OneGate vector to end of AllGates vector */
+ OnePeak_AllGates.push_back(AllPeaks_OneGate[indexE]);
+ //delete c_peakFits;
+ }
+
+ /* Write PS-subtracted spectrum to file */
+ //TFile* outfile = new TFile("GateThetaLab_ExMinusGatePhaseSpace.root","RECREATE");
+ //list->Write("GateThetaLab_ExMinusPhaseSpace",TObject::kSingleKey);
+
+ return OnePeak_AllGates;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+TH1F* PullThetaLabHist(int i, double minTheta, double gatesize){
+ string name = "./C2S_Outputs/GateThetaLabHist_Experiment.root";
+ TFile* file = new TFile(name.c_str(),"READ");
+ string histname = "TLab"
+ + to_string((int)(minTheta+(i*gatesize))) + "to"
+ + to_string((int)(minTheta+((i+1)*gatesize)));
+ cout << "Loading " << name << " -> " << histname << endl;
+
+ TList *list = (TList*)file->Get("GateThetaLabHistograms");
+ TH1F* hist = (TH1F*)list->FindObject(histname.c_str());
+
+ double pulledBin = hist->GetBinWidth(10);
+ int ratio = (int)(globalBinning / pulledBin);
+ hist->Rebin(ratio);
+ file->Close();
+
+ return hist;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+TH1F* PullThetaCMHist(int i, double minTheta, double gatesize){
+ //TFile* file = new TFile("GateThetaCMHistograms_47Kdt_18Oct22_bin0p1.root","READ");
+ //string name = "GateThetaCMHistograms_47Kdt_18Oct22_";
+ string name = "GateThetaCMHistograms_47Kdt_09Jun23";
+ if(!GammaGate){
+ //name = name + "bin0p1.root";
+ name = name + ".root";
+ } else {
+ name = name + to_string((int)(globGammaGate*1000))
+ + "GammaGate.root";
+ }
+ TFile* file = new TFile(name.c_str(),"READ");
+
+ string histname = "cThetaCMGate_"
+ + to_string((int)(minTheta+(i*gatesize))) + "-"
+ + to_string((int)(minTheta+((i+1)*gatesize)));
+ cout << "Loading " << histname << endl;
+ TList *list = (TList*)file->Get("GateThetaCMHistograms");
+ TH1F* hist = (TH1F*)list->FindObject(histname.c_str());
+
+ double pulledBin = hist->GetBinWidth(10);
+ int ratio = (int)(globalBinning / pulledBin);
+ hist->Rebin(ratio);
+ file->Close();
+
+ return hist;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+TH1F* PullPhaseSpaceHist(int i, double minTheta, double gatesize){
+ string name = "./C2S_Outputs/GateThetaLabHist_PhaseSpace.root";
+ TFile* file = new TFile(name.c_str(),"READ");
+ string histname = "TLab"
+ + to_string((int) (minTheta+(i*gatesize))) + "to"
+ + to_string((int) (minTheta+((i+1)*gatesize)));
+ cout << "Loading " << name << " -> " << histname << endl;
+
+ TList *list = (TList*)file->Get("GateThetaLabHistograms");
+ TH1F* hist = (TH1F*)list->FindObject(histname.c_str());
+ file->Close();
+
+ double pulledBin = hist->GetBinWidth(10);
+ int ratio = (int)(globalBinning / pulledBin);
+ hist->Rebin(ratio);
+
+ return hist;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Scale(TGraph* g , TGraphErrors* ex){
+ double scale;
+ double mean = 0 ;
+ double* eX = ex->GetX();
+ double* eY = ex->GetY();
+ double totalW = 0;
+ double W = 0;
+ for(int i = 0 ; i < ex->GetN() ; i++){
+ if(eY[i]>1 && eY[i] <1e4){
+ // Incremental Error weighted average
+ W = 1./ex->GetErrorY(i);
+ scale = eY[i]/g->Eval(eX[i]);
+ totalW +=W;
+ mean = mean + W*(scale - mean)/(totalW);
+ }
+ }
+
+ //scaleTWOFNR = mean;
+ if(loud){
+ cout << "SCALED THEORY BY " << mean << endl;
+ cout << " therefore S = " << 1/mean << " ??" << endl;
+ }
+
+ double* x = g->GetX();
+ double* y = g->GetY();
+ for(unsigned int i = 0 ; i < g->GetN() ; i++)
+ g->SetPoint(i,x[i],y[i]*mean);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+double Chi2(TGraph* theory, TGraphErrors* exper){
+ double Chi2 = 0;
+ double chi = 0;
+
+ //cout << setprecision(8);
+ //for(int i = 1 ; i < exper->GetN() ; i++){
+ for(int i = 0 ; i < exper->GetN() ; i++){
+ if(exper->GetPointY(i)>1.0e-4){//10){ //0){
+// cout << setprecision(6)
+// << "COMPARE::::: "
+// << exper->GetPointY(i) << " +- " << exper->GetErrorY(1)
+// << " TO "
+// << theory->Eval(anglecentres[i])
+// << endl;
+
+ //////////////////////////////////////////////////////////////////////////////////////////
+ // ---------------------------------------------------------------------------------------
+ // | Louis Lyons pg 102:
+ // | Least squares should minimize sum(y_obs - y_th / err)^2
+ // | Err SHOULD in principle be error in THEORY
+ // | However it is common practive to use error in OBS for various reasons (pg 103-104)
+ // ---------------------------------------------------------------------------------------
+ //////////////////////////////////////////////////////////////////////////////////////////
+
+ //Chi2 += pow((exper->GetPointY(i) - theory->Eval(anglecentres[i])), 2 ) / theory->Eval(anglecentres[i]);
+ Chi2 += pow(((exper->GetPointY(i) - theory->Eval(anglecentres[i]))) / exper->GetErrorY(i) , 2 );
+ }
+ }
+ if(loud){cout << setprecision(10) << "Chi2 = " << Chi2 << " #dof = " << exper->GetN()-1 << " Chi2/#dof = " << Chi2/(double)exper->GetN()-1 << setprecision(4)<< endl;}
+ return Chi2;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+double ToMininize(const double* parameter){
+ static int f = 0 ;
+ TGraph* g = new TGraph();
+ if(doMix){
+ double* X = currentThry->GetX(); // gets valies from global g1 = tgraph passed to find norm
+ double* Y1 = currentThry->GetY();
+ double* Y2 = currentThry2->GetY();
+ for(int i = 0 ; i < currentThry->GetN() ; i++){
+ g->SetPoint(g->GetN(),X[i],parameter[0]*Y1[i] + parameter[1]*Y2[i]); // scales this tgraph by parameter
+ }
+ } else {
+ double* X = currentThry->GetX(); // gets valies from global g1 = tgraph passed to find norm
+ double* Y = currentThry->GetY();
+ for(int i = 0 ; i < currentThry->GetN() ; i++){
+ g->SetPoint(g->GetN(),X[i],parameter[0]*Y[i]); // scales this tgraph by parameter
+ }
+ }
+
+ double chi2 = Chi2(g,staticExp); //compares theory tgraph to experimental tgrapherrors by chi2
+ return chi2;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+TGraph* FindNormalisation(TGraph* theory, TGraphErrors* experiment){
+ /* (dSdO)meas = S * (dSdO)calc */
+
+ // Set global variable
+ currentThry = theory;
+ staticExp = experiment;
+
+ // Construct minimiser
+ const char* minName ="Minuit";const char* algoName="Migrad";
+ ROOT::Math::Minimizer* min =
+ ROOT::Math::Factory::CreateMinimizer(minName, algoName);
+ min->SetValidError(true);
+
+ // Number of parameters (should only be 1 for me)
+ int mysize = 1;
+
+ // Create funciton wrapper for minimizer
+ // a IMultiGenFunction type
+ ROOT::Math::Functor f(&ToMininize,mysize);
+ min->SetFunction(f);
+
+ //int SpecFactorUpperLimit = 2;
+ int SpecFactorUpperLimit = 20;
+
+ // Set range of parameter(??)
+ double* parameter = new double[mysize];
+ for(unsigned int i = 0 ; i < mysize ; i++){
+ parameter[i] = 0.5;
+ char name[4];
+ sprintf(name,"S%d",i+1);
+ min->SetLimitedVariable(i,name,parameter[i],0.01,0,SpecFactorUpperLimit);
+ }
+
+ //min->SetPrintLevel(3);
+
+ ///// TO IMPROVE: FIND WAY OF OBTAINING NDF AND PRINT CHI2/NDF /////
+
+ // Minimise
+ min->ProvidesError();
+ min->SetValidError(true);
+ min->Minimize();
+ const double *xs = min->X();
+ const double *err = min->Errors();
+
+ cout << "IS THE MINIMIZATION GOOD? -> " << min->IsValidError() << endl;
+
+ // Write out
+ for(int i = 0 ; i < mysize ; i++){
+ cout << Form("S%d=",i+1) << xs[i] << "(" << err[i] << ")" << endl;
+ }
+
+ /* Store S value in global variable, to access for drawing on plots */
+ globalS = xs[0];
+ globalSerr = err[0];
+
+ // Return the Fitted CS
+ TGraph* g = new TGraph();
+ double* X = theory->GetX();
+ double* Y = theory->GetY();
+ if(loud){
+ cout << setprecision(8);
+ cout << "Start: X[0] = " << theory->GetPointX(4) << " Y[0] = " << theory->GetPointY(4) << endl;
+ cout << "multip by " << xs[0] << endl;
+ }
+
+ for(int i=0; i<theory->GetN(); i++){ g->SetPoint(i,X[i],xs[0]*Y[i]); }
+
+ if(loud){
+ cout << "End: X[0] = " << g->GetPointX(4) << " Y[0] = " << g->GetPointY(4) << endl;
+ cout << setprecision(3);
+ }
+
+ return g;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+TGraph* FindNormalisation(TGraph* theory, TGraph* theory2, TGraphErrors* experiment){
+ /* (dSdO)meas = S * (dSdO)calc */
+
+ // Set global variable
+ currentThry = theory;
+ currentThry2 = theory2;
+ staticExp = experiment;
+
+ // Construct minimiser
+ const char* minName ="Minuit";const char* algoName="Migrad";
+ ROOT::Math::Minimizer* min =
+ ROOT::Math::Factory::CreateMinimizer(minName, algoName);
+ min->SetValidError(true);
+
+ // Number of parameters (should only be 1 for me)
+ int mysize = 1;
+ if(doMix){mysize=2;}
+
+ // Create funciton wrapper for minimizer
+ // a IMultiGenFunction type
+ ROOT::Math::Functor f(&ToMininize,mysize);
+ min->SetFunction(f);
+
+ // Set range of parameter(??)
+ double* parameter = new double[mysize];
+ for(unsigned int i = 0 ; i < mysize ; i++){
+ //parameter[i] = 0.8;
+ parameter[i] = 10.;
+ char name[4];
+ sprintf(name,"S%d",i+1);
+ //min->SetLimitedVariable(i,name,parameter[i],0.01,0,10);
+ min->SetLimitedVariable(i,name,parameter[i],0.01,0,20);
+ }
+
+ // Minimise
+ min->Minimize();
+ const double *xs = min->X();
+ const double *err = min->Errors();
+
+ // Write out
+ for(int i = 0 ; i < mysize ; i++){
+ cout << Form("S%d=",i+1) << xs[i] << "(" << err[i] << ")" << endl;
+ }
+
+ /* Store S value in global variable, to access for drawing on plots */
+ globalS = xs[0];
+ globalSerr = err[0];
+ globalS2 = xs[1];
+ globalSerr2 = err[1];
+
+ // Return the Fitted CS
+ TGraph* g = new TGraph();
+ double* X = theory->GetX();
+ double* Y1 = theory ->GetY();
+ double* Y2 = theory2->GetY();
+ if(loud){
+ cout << setprecision(8);
+ cout << "Start: X[0] = " << theory->GetPointX(4) << " Y[0] = " << theory->GetPointY(4) << endl;
+ cout << "multip by " << xs[0] << endl;
+ cout << "Start: X[1] = " << theory2->GetPointX(4) << " Y[0] = " << theory2->GetPointY(4) << endl;
+ cout << "multip by " << xs[1] << endl;
+ }
+
+ for(int i=0; i<theory->GetN(); i++){
+ g->SetPoint(i, X[i], xs[0]*Y1[i] + xs[1]*Y2[i]);
+ }
+
+ if(loud){
+ cout << "End: X = " << g->GetPointX(4) << " Y = " << g->GetPointY(4) << endl;
+ cout << setprecision(3);
+ }
+
+ return g;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void C2S(double Energy, double Spin, double spdf, double angmom, double binning, double GammaGateEnergy){
+ GammaGate=true;
+ globGammaGate = GammaGateEnergy;
+
+ C2S(Energy, Spin, spdf, angmom, binning, "");
+}
diff --git a/Projects/e775s/analysis/Plots_Diagnostics.C b/Projects/e775s/analysis/Plots_Diagnostics.C
new file mode 100644
index 0000000000000000000000000000000000000000..0a117540fd19b23739cb29dc489ac59ea711822e
--- /dev/null
+++ b/Projects/e775s/analysis/Plots_Diagnostics.C
@@ -0,0 +1,30 @@
+#include "Plots_Functions.h"
+
+void Plots_Diagnostics(){
+
+
+//TCutG* cutGood = LoadCut_Good();
+//TCutG* cutBad = LoadCut_Bad();
+
+
+//TFile *f = TFile::Open("../../../Outputs/Analysis/Testing_NoStripsExcluded.root");
+//TFile *f = TFile::Open("../../../Outputs/Analysis/Testing_StripExclusion.root");
+//TFile *f = TFile::Open("../../../Outputs/Analysis/Testing_StripExclusionUsingMapping.root");
+//TFile *f = TFile::Open("../../../Outputs/Analysis/Testing_StripExclusionUsingMapping_V4.root");
+//TFile *f = TFile::Open("../../../Outputs/Analysis/Testing_StripExclusionUsingMapping_V5.root");
+//TFile *f = TFile::Open("../../../Outputs/Analysis/Testing_Thresholds.root");
+//TFile *f = TFile::Open("../../../Outputs/Analysis/Testing_BSpotMin.root");
+
+ vector<string> files;
+ files.push_back("../../../Outputs/Analysis/Testing_ExCalib.root");
+ tree = Chain("PhysicsTree",files,true);
+
+ vector<string> filesau;
+ filesau.push_back("../../../Outputs/Analysis/Testing_ExCalib_Au.root");
+ tree_au = Chain("PhysicsTree",filesau,true);
+
+ vector<string> filescd;
+ filescd.push_back("../../../Outputs/Analysis/Testing_ExCalib_CD.root");
+ tree_cd = Chain("PhysicsTree",filescd,true);
+
+}
diff --git a/Projects/e775s/analysis/Plots_Experiment.C b/Projects/e775s/analysis/Plots_Experiment.C
new file mode 100644
index 0000000000000000000000000000000000000000..1ed0a488e2fd06d1760d4c920a2cceb245c186b0
--- /dev/null
+++ b/Projects/e775s/analysis/Plots_Experiment.C
@@ -0,0 +1,42 @@
+#include "Plots_Functions.h"
+#include "Plots_C2SFunctions.h"
+#include "Plots_ParryFunctions.h"
+
+
+void C2S(){
+ cout << "======================================" << endl;
+ cout << " Run C2S:" << endl;
+ cout << " -- 0.000 MeV, 0+, d5/2: C2S(0.000,0,2,2.5,0.05,\"\")" << endl;
+ cout << " -- 1.675 MeV, 2+, d5/2: C2S(1.675,2,2,2.5,0.05,\"\")" << endl;
+ cout << " -- 3.572 MeV, 4+, d5/2: C2S(3.572,4,2,2.5,0.05,\"\")" << endl;
+ cout << " -- 4.071 MeV, 2+, mixed?: C2S(4.071,2,0,0.5,0.05,\"mixed\")" << endl;
+ cout << " -- 5.004 MeV, ?+, xx/2: C2S(5.004,2,2,2.5,0.05,\"\")" << endl;
+ cout << " -- 5.228 MeV, ?+, s1/2: C2S(5.228,2,0,0.5,0.05,\"\")" << endl;
+ cout << " -- 5.629 MeV, ?+, xx/2: C2S(5.629,2,2,2.5,0.05,\"\")" << endl;
+ cout << "======================================" << endl;
+}
+
+
+void Plots_Experiment(){
+
+ //Load in files
+//LoadChain_Au();
+//LoadChain_CD2();
+//LoadChain_Both();
+ LoadChain_Both_PaxmanAnalysis();
+// TCutG* cutPrompt = LoadCut_Prompt();
+// TCutG* cutTimeBoth = LoadCut_TimeBoth();
+// TCutG* cutBothBlobs = LoadCut_BothBlobs();
+// TCutG* cutPromptPlus = LoadCut_PromptPlus();
+// TCutG* cutGood = LoadCut_Good();
+// TCutG* cutBad = LoadCut_Bad();
+
+ cout << "======================================" << endl;
+ cout << "Time gate GOOD = abs(T_FPMW_HF-13000)<1000" << endl;
+ cout << "Time gate BAD = abs(T_FPMW_HF-10200)< 600" << endl;
+ cout << "======================================" << endl;
+ cout << "USEFUL COMMANDS:::" << endl;
+ cout << " GateGamma_SeeParticle(0.005,0.10,0.09,0.03,0.14,0.03)" << endl;
+ cout << " --First excited decay of 19O, with background" << endl;
+ C2S();
+}
diff --git a/Projects/e775s/analysis/Plots_Functions.h b/Projects/e775s/analysis/Plots_Functions.h
new file mode 100644
index 0000000000000000000000000000000000000000..f6b86a5382d8472ee49c5914c5e22bfbe3fce0ce
--- /dev/null
+++ b/Projects/e775s/analysis/Plots_Functions.h
@@ -0,0 +1,1083 @@
+#include <math.h>
+#include <string>
+#include <sstream>
+#include "Plots_TWOFNR.h"
+#include "Plots_KnownPeakFitter.h"
+#include "DefineColours.h"
+using namespace std;
+
+//=====================================================================//
+//=====================================================================//
+
+vector<string> files;
+TChain* tree_cd=NULL;
+TChain* tree_au=NULL;
+TChain* tree_both=NULL;
+TChain* tree_raw=NULL;
+TChain* tree=NULL;
+TChain* tree_PSp00=NULL;
+TChain* tree_PSp96=NULL;
+
+NPL::Reaction O19dp("19O(d,p)20O@151"); //TEMPORARY BEAM ENERGY!
+
+auto colourR = new TColor((Int_t)TColor::GetFreeColorIndex(), 220./256., 5./256., 12./256.);
+auto colourB = new TColor((Int_t)TColor::GetFreeColorIndex(), 25./256., 101./256., 176./256.);
+auto colourG = new TColor((Int_t)TColor::GetFreeColorIndex(), 78./256., 178./256., 101./256.);
+auto colourP = new TColor((Int_t)TColor::GetFreeColorIndex(), 136./256., 46./256., 114./256.);
+
+
+//=====================================================================//
+//=====================================================================//
+
+//////////////////////////////////////////////////////////////////////////////////////
+void PlotKinematic(NPL::Reaction r, double Ex, Color_t c, int w, int s)
+{
+ r.SetExcitation4(Ex);
+ TGraph* g = r.GetKinematicLine3();
+ g->SetLineColor(c);
+ g->SetLineWidth(w);
+ g->SetLineStyle(s);
+ g->Draw("c same");
+}
+
+//////////////////////////////////////////////////////////////////////////////////////
+TChain* Chain(string TreeName, vector<string>& file, bool EventList)
+{
+ TChain* chain = new TChain(TreeName.c_str());
+ unsigned int size = file.size();
+ for(unsigned int i=0; i<size; i++){
+ cout << "Adding to chain: " << file[i] << endl;
+ chain->Add(file[i].c_str());
+ }
+ return chain;
+}
+
+//////////////////////////////////////////////////////////////////////////////////////
+void LoadChain_Both_PaxmanAnalysis()
+{
+ files.clear();
+//files.push_back("~/Programs/nptoolv3/Outputs/Analysis/19Odp_2024-07-30_CD.root");
+//files.push_back("~/Programs/nptoolv3/Outputs/Analysis/19Odp_2024-07-30_Au.root");
+ files.push_back("~/Programs/nptoolv3/Outputs/Analysis/19Odp_2024-08-12_CD.root");
+ files.push_back("~/Programs/nptoolv3/Outputs/Analysis/19Odp_2024-08-12_Au.root");
+ tree_both = Chain("PhysicsTree",files,true);
+}
+
+//////////////////////////////////////////////////////////////////////////////////////
+void LoadChain_Raw()
+{
+ files.clear();
+ files.push_back("~/Experiment/e775s/RawData/run_0115/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0116/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0117/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0123/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0124/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0125/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0128/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0129/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0130/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0131/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0132/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0133/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0134/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0135/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0137/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0138/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0139/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0140/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0142/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0150/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0151/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0152/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0153/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0163/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0164/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0165/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0166/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0167/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0168/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0169/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0171/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0172/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0173/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0174/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0178/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0179/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0180/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0181/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0184/Out/Analysis/Tree_0000.root");
+ files.push_back("~/Experiment/e775s/RawData/run_0185/Out/Analysis/Tree_0000.root");
+ tree_raw = Chain("TreeMaster",files,true);
+}
+
+//////////////////////////////////////////////////////////////////////////////////////
+void LoadChain_PhaseSpace()
+{
+ files.clear();
+ files.push_back("../../../Outputs/Analysis/PSp_PhaseSpace_InitialTest_0000.root");
+ tree_PSp00 = Chain("PhysicsTree",files,true);
+
+ files.clear();
+ files.push_back("../../../Outputs/Analysis/PSp_PhaseSpace_InitialTest_0096.root");
+ tree_PSp96 = Chain("PhysicsTree",files,true);
+}
+
+//////////////////////////////////////////////////////////////////////////////////////
+void DrawPlotsToLoad(TChain* tree)
+{
+ // ExEg
+ tree->Draw("Ex:AddBack_EDC>>hExEg(8000,0,8,2000,-5,15)","abs(T_FPMW_HF-13000)<1000 && ThetaLab<158","colz");
+ TH2F* hExEg = (TH2F*) gDirectory->Get("hExEg");
+
+ // ExEt
+ tree->Draw("Ex:trackEDC>>hExEt(8000,0,8,2000,-5,15)","abs(T_FPMW_HF-13000)<1000 && ThetaLab<158","colz");
+ TH2F* hExEt= (TH2F*) gDirectory->Get("hExEt");
+
+ // Ex (cleanup)
+ tree->Draw("Ex>>hExGood(2000,-5,15)","abs(T_FPMW_HF-13000)<1000 && ThetaLab<158","colz");
+ TH1F* hExGood = (TH1F*) gDirectory->Get("hExGood");
+ tree->Draw("Ex>>hExBad(2000,-5,15)","abs(T_FPMW_HF-10200)<600 && ThetaLab<158","colz");
+ TH1F* hExBad = (TH1F*) gDirectory->Get("hExBad");
+ TH1F* hExClean = (TH1F*) hExGood->Clone();
+ hExClean->SetTitle("hExClean");
+ hExClean->SetName("hExClean");
+ hExClean->Add(hExBad,-1);
+ hExClean->Draw();
+
+ // Ex-ThetaLab (cleanup)
+ tree->Draw("Ex:ThetaLab>>hExTLGood(80,100,180,2000,-5,15)","abs(T_FPMW_HF-13000)<1000 && ThetaLab<158","colz");
+ TH2F* hExTLGood = (TH2F*) gDirectory->Get("hExTLGood");
+ tree->Draw("Ex:ThetaLab>>hExTLBad(80,100,180,2000,-5,15)","abs(T_FPMW_HF-10200)<600 && ThetaLab<158","colz");
+ TH2F* hExTLBad = (TH2F*) gDirectory->Get("hExTLBad");
+ TH2F* hExTLClean = (TH2F*) hExTLGood->Clone();
+ hExTLClean->SetTitle("hExTLClean");
+ hExTLClean->SetName("hExTLClean");
+ hExTLClean->Add(hExTLBad,-1);
+ hExTLClean->Draw();
+
+}
+
+//////////////////////////////////////////////////////////////////////////////////////
+void DrawPlotsToLoad()
+{
+ DrawPlotsToLoad(tree_both);
+
+}
+//=====================================================================//
+//=====================================================================//
+
+//////////////////////////////////////////////////////////////////////////////////////
+TH1F* Project(TH2F* h2D, const char * XY, Color_t col)
+{
+ TH1F* h;
+ if(string(XY) == "X"){
+ string name = "projX_" + to_string(col);
+ h2D->ProjectionX(name.c_str(),0,-1,"");
+ h = (TH1F*) gDirectory->Get(name.c_str());
+ h->SetLineColor(col);
+ } else if (string(XY) == "Y"){
+ string name = "projY_" + to_string(col);
+ h2D->ProjectionY(name.c_str(),0,-1,"");
+ h = (TH1F*) gDirectory->Get(name.c_str());
+ h->SetLineColor(col);
+ } else {
+ cout << "!!!! ERROR !!!! XY projection selector is invalid" << endl;
+ return 0;
+ }
+ return (TH1F*) h;
+};
+
+//////////////////////////////////////////////////////////////////////////////////////
+TH1F* Project(TH2F* h2D, const char * XY, double mean, double width, Color_t col)
+{
+ TH1F* h;
+ if(string(XY) == "X"){
+ //X projection is GPSG
+ string name = "projX_" + to_string(col);
+ h2D->ProjectionX(name.c_str(),h2D->GetYaxis()->FindBin(mean-(0.5*width)),
+ h2D->GetYaxis()->FindBin(mean+(0.5*width)),"");
+ h = (TH1F*) gDirectory->Get(name.c_str());
+ h->SetLineColor(col);
+ } else if (string(XY) == "Y"){
+ //Y projection is GGSP
+ string name = "projY_" + to_string(col);
+ h2D->ProjectionY(name.c_str(),h2D->GetXaxis()->FindBin(mean-(0.5*width)),
+ h2D->GetXaxis()->FindBin(mean+(0.5*width)),"");
+ h = (TH1F*) gDirectory->Get(name.c_str());
+ h->SetLineColor(col);
+ } else {
+ cout << "!!!! ERROR !!!! XY projection selector is invalid" << endl;
+ return 0;
+ }
+ return (TH1F*) h;
+};
+
+//////////////////////////////////////////////////////////////////////////////////////
+void Prj_DrawLines(double mean, double width, Color_t col)
+{
+ TLine* line = new TLine();
+ line->SetLineColor(col);
+ line->SetLineStyle(kDashed);
+ line->SetLineWidth(2);
+ line->DrawLine(mean-(0.5*width), 0, mean-(0.5*width), 2000);
+ line->DrawLine(mean+(0.5*width), 0, mean+(0.5*width), 2000);
+}
+
+//////////////////////////////////////////////////////////////////////////////////////
+string AddOrTrack(string input)
+{
+ if(input == "addback"){ return "hExEg"; }
+ else if(input == "track"){ return "hExEt"; }
+ else { cout << "ERROR! not 'addback' or 'track' !!!" << endl; return "hExEg";}
+}
+
+//////////////////////////////////////////////////////////////////////////////////////
+void GateGamma_SeeParticle(string addbackOrTrack, double bing, double binp, double Eg1, double width1)
+{
+ string histname = AddOrTrack(addbackOrTrack);
+ //Check if histogram exists alreado, else draw
+ //in future, add check, then try to load, then draw
+ if(gDirectory->Get(histname.c_str()) == 0){
+ DrawPlotsToLoad();
+ }
+
+ //Load in & rebin
+ TH2F* hCopyMe = (TH2F*) gDirectory->Get(histname.c_str());
+ TH2F* hExEg = (TH2F*) hCopyMe->Clone();
+ hExEg->Rebin2D(bing/hExEg->GetXaxis()->GetBinWidth(3),
+ binp/hExEg->GetYaxis()->GetBinWidth(3));
+
+ //Set up canvas
+ TCanvas *cGGSP = new TCanvas("cGGSP","cGGSP",1000,750);
+ cGGSP->Divide(1,2);
+ cGGSP->cd(1);
+
+ //Axis being projected
+ TH1F* projG = Project(hExEg, "X", kBlack);
+ projG->Draw();
+
+ //Projection gates
+ Prj_DrawLines(Eg1, width1, kRed);
+
+ //Projected axis
+ cGGSP->cd(2);
+ TH1F* projP1 = Project(hExEg, "Y", Eg1, width1, kRed);
+ projP1->Draw();
+}
+
+//////////////////////////////////////////////////////////////////////////////////////
+void GateGamma_SeeParticle(string addbackOrTrack, double bing, double binp, double Eg1, double width1,
+ double Eg2, double width2)
+{
+ string histname = AddOrTrack(addbackOrTrack);
+ //Check if histogram exists alreado, else draw
+ //in future, add check, then try to load, then draw
+ if(gDirectory->Get(histname.c_str()) == 0){
+ DrawPlotsToLoad();
+ }
+
+ //Load in & rebin
+ TH2F* hCopyMe = (TH2F*) gDirectory->Get(histname.c_str());
+ TH2F* hExEg = (TH2F*) hCopyMe->Clone();
+ hExEg->Rebin2D(bing/hExEg->GetXaxis()->GetBinWidth(3),
+ binp/hExEg->GetYaxis()->GetBinWidth(3));
+
+ //Set up canvas
+ TCanvas *cGGSP = new TCanvas("cGGSP","cGGSP",1000,750);
+ cGGSP->Divide(1,2);
+ cGGSP->cd(1);
+
+ //Axis being projected
+ TH1F* projG = Project(hExEg, "X", kBlack);
+ projG->Draw();
+
+ //Projection gates
+ Prj_DrawLines(Eg1, width1, kRed);
+ Prj_DrawLines(Eg2, width2, kBlue);
+
+ //Projected axis
+ cGGSP->cd(2);
+ TH1F* projP1 = Project(hExEg, "Y", Eg1, width1, kRed);
+ TH1F* projP2 = Project(hExEg, "Y", Eg2, width2, kBlue);
+ projP1->Draw();
+ projP2->Draw("same");
+}
+
+//////////////////////////////////////////////////////////////////////////////////////
+void GateGamma_SeeParticle(string addbackOrTrack, double bing, double binp, double Eg1, double width1,
+ double Eg2, double width2,
+ double Eg3, double width3)
+{
+ string histname = AddOrTrack(addbackOrTrack);
+ //Check if histogram exists alreado, else draw
+ //in future, add check, then try to load, then draw
+ if(gDirectory->Get(histname.c_str()) == 0){
+ DrawPlotsToLoad();
+ }
+
+ //Load in & rebin
+ TH2F* hCopyMe = (TH2F*) gDirectory->Get(histname.c_str());
+ TH2F* hExEg = (TH2F*) hCopyMe->Clone();
+ hExEg->Rebin2D(bing/hExEg->GetXaxis()->GetBinWidth(3),
+ binp/hExEg->GetYaxis()->GetBinWidth(3));
+
+ //Set up canvas
+ TCanvas *cGGSP = new TCanvas("cGGSP","cGGSP",1000,750);
+ cGGSP->Divide(1,2);
+ cGGSP->cd(1);
+
+ //Axis being projected
+ TH1F* projG = Project(hExEg, "X", kBlack);
+ projG->Draw();
+
+ //Projection gates
+ Prj_DrawLines(Eg1, width1, kRed);
+ Prj_DrawLines(Eg2, width2, kBlue);
+ Prj_DrawLines(Eg3, width3, kGreen);
+
+ //Projected axis
+ cGGSP->cd(2);
+ TH1F* projP1 = Project(hExEg, "Y", Eg1, width1, kRed);
+ TH1F* projP2 = Project(hExEg, "Y", Eg2, width2, kBlue);
+ TH1F* projP3 = Project(hExEg, "Y", Eg3, width3, kGreen);
+ projP1->Draw();
+ projP2->Draw("same");
+ projP3->Draw("same");
+}
+
+//////////////////////////////////////////////////////////////////////////////////////
+void GateGamma_SeeParticle(string addbackOrTrack, double bing, double binp, double Eg1, double width1,
+ double Eg2, double width2,
+ double Eg3, double width3,
+ double Eg4, double width4)
+{
+ string histname = AddOrTrack(addbackOrTrack);
+ //Check if histogram exists alreado, else draw
+ //in future, add check, then try to load, then draw
+ if(gDirectory->Get(histname.c_str()) == 0){
+ DrawPlotsToLoad();
+ }
+
+ //Load in & rebin
+ TH2F* hCopyMe = (TH2F*) gDirectory->Get(histname.c_str());
+ TH2F* hExEg = (TH2F*) hCopyMe->Clone();
+ hExEg->Rebin2D(bing/hExEg->GetXaxis()->GetBinWidth(3),
+ binp/hExEg->GetYaxis()->GetBinWidth(3));
+
+ //Set up canvas
+ TCanvas *cGGSP = new TCanvas("cGGSP","cGGSP",1000,750);
+ cGGSP->Divide(1,2);
+ cGGSP->cd(1);
+
+ //Axis being projected
+ TH1F* projG = Project(hExEg, "X", kBlack);
+ projG->Draw();
+
+ //Projection gates
+ Prj_DrawLines(Eg1, width1, kRed);
+ Prj_DrawLines(Eg2, width2, kBlue);
+ Prj_DrawLines(Eg3, width3, kGreen);
+ Prj_DrawLines(Eg4, width4, kViolet);
+
+ //Projected axis
+ cGGSP->cd(2);
+ TH1F* projP1 = Project(hExEg, "Y", Eg1, width1, kRed);
+ TH1F* projP2 = Project(hExEg, "Y", Eg2, width2, kBlue);
+ TH1F* projP3 = Project(hExEg, "Y", Eg3, width3, kGreen);
+ TH1F* projP4 = Project(hExEg, "Y", Eg4, width4, kViolet);
+ projP1->Draw();
+ projP2->Draw("same");
+ projP3->Draw("same");
+ projP4->Draw("same");
+}
+
+//////////////////////////////////////////////////////////////////////////////////////
+void GateGamma_SeeParticle(string addbackOrTrack, double bing, double binp, double Eg1, double width1,
+ double Eg2, double width2,
+ double Eg3, double width3,
+ double Eg4, double width4,
+ double Eg5, double width5)
+{
+ string histname = AddOrTrack(addbackOrTrack);
+ //Check if histogram exists alreado, else draw
+ //in future, add check, then try to load, then draw
+ if(gDirectory->Get(histname.c_str()) == 0){
+ DrawPlotsToLoad();
+ }
+
+ //Load in & rebin
+ TH2F* hCopyMe = (TH2F*) gDirectory->Get(histname.c_str());
+ TH2F* hExEg = (TH2F*) hCopyMe->Clone();
+ hExEg->Rebin2D(bing/hExEg->GetXaxis()->GetBinWidth(3),
+ binp/hExEg->GetYaxis()->GetBinWidth(3));
+
+ //Set up canvas
+ TCanvas *cGGSP = new TCanvas("cGGSP","cGGSP",1000,750);
+ cGGSP->Divide(1,2);
+ cGGSP->cd(1);
+
+ //Axis being projected
+ TH1F* projG = Project(hExEg, "X", kBlack);
+ projG->Draw();
+
+ //Projection gates
+ Prj_DrawLines(Eg1, width1, kRed);
+ Prj_DrawLines(Eg2, width2, kBlue);
+ Prj_DrawLines(Eg3, width3, kGreen);
+ Prj_DrawLines(Eg4, width4, kViolet);
+ Prj_DrawLines(Eg5, width5, kOrange);
+
+ //Projected axis
+ cGGSP->cd(2);
+ TH1F* projP1 = Project(hExEg, "Y", Eg1, width1, kRed);
+ TH1F* projP2 = Project(hExEg, "Y", Eg2, width2, kBlue);
+ TH1F* projP3 = Project(hExEg, "Y", Eg3, width3, kGreen);
+ TH1F* projP4 = Project(hExEg, "Y", Eg4, width4, kViolet);
+ TH1F* projP5 = Project(hExEg, "Y", Eg5, width5, kOrange);
+ projP1->Draw();
+ projP2->Draw("same");
+ projP3->Draw("same");
+ projP4->Draw("same");
+ projP5->Draw("same");
+}
+
+//////////////////////////////////////////////////////////////////////////////////////
+void GateParticle_SeeGamma(string addbackOrTrack, double bing, double binp, double Ex1, double width1)
+{
+ string histname = AddOrTrack(addbackOrTrack);
+ //Check if histogram exists alreado, else draw
+ //in future, add check, then try to load, then draw
+ if(gDirectory->Get(histname.c_str()) == 0){
+ DrawPlotsToLoad();
+ }
+
+ //Load in & rebin
+ TH2F* hCopyMe = (TH2F*) gDirectory->Get(histname.c_str());
+ TH2F* hExEg = (TH2F*) hCopyMe->Clone();
+ hExEg->Rebin2D(bing/hExEg->GetXaxis()->GetBinWidth(3),
+ binp/hExEg->GetYaxis()->GetBinWidth(3));
+
+ //Set up canvas
+ TCanvas *cGPSG = new TCanvas("cGPSG","cGPSG",1000,750);
+ cGPSG->Divide(1,2);
+ cGPSG->cd(1);
+
+ //Axis being projected
+ TH1F* projP = Project(hExEg, "Y", kBlack);
+ projP->Draw();
+
+ //Projection gates
+ Prj_DrawLines(Ex1, width1, kRed);
+
+ //Projected axis
+ cGPSG->cd(2);
+ TH1F* projG1 = Project(hExEg, "X", Ex1, width1, kRed);
+ projG1->Draw();
+}
+
+//////////////////////////////////////////////////////////////////////////////////////
+void GateParticle_SeeGamma(string addbackOrTrack, double bing, double binp, double Ex1, double width1,
+ double Ex2, double width2)
+{
+ string histname = AddOrTrack(addbackOrTrack);
+ //Check if histogram exists alreado, else draw
+ //in future, add check, then try to load, then draw
+ if(gDirectory->Get(histname.c_str()) == 0){
+ DrawPlotsToLoad();
+ }
+
+ //Load in & rebin
+ TH2F* hCopyMe = (TH2F*) gDirectory->Get(histname.c_str());
+ TH2F* hExEg = (TH2F*) hCopyMe->Clone();
+ hExEg->Rebin2D(bing/hExEg->GetXaxis()->GetBinWidth(3),
+ binp/hExEg->GetYaxis()->GetBinWidth(3));
+
+ //Set up canvas
+ TCanvas *cGPSG = new TCanvas("cGPSG","cGPSG",1000,750);
+ cGPSG->Divide(1,2);
+ cGPSG->cd(1);
+
+ //Axis being projected
+ TH1F* projP = Project(hExEg, "Y", kBlack);
+ projP->Draw();
+
+ //Projection gates
+ Prj_DrawLines(Ex1, width1, kRed);
+ Prj_DrawLines(Ex2, width2, kBlue);
+
+ //Projected axis
+ cGPSG->cd(2);
+ TH1F* projG1 = Project(hExEg, "X", Ex1, width1, kRed);
+ TH1F* projG2 = Project(hExEg, "X", Ex2, width2, kBlue);
+ projG1->Draw();
+ projG2->Draw("same");
+}
+
+//////////////////////////////////////////////////////////////////////////////////////
+void GateParticle_SeeGamma(string addbackOrTrack, double bing, double binp, double Ex1, double width1,
+ double Ex2, double width2,
+ double Ex3, double width3)
+{
+ string histname = AddOrTrack(addbackOrTrack);
+ //Check if histogram exists alreado, else draw
+ //in future, add check, then try to load, then draw
+ if(gDirectory->Get(histname.c_str()) == 0){
+ DrawPlotsToLoad();
+ }
+
+ //Load in & rebin
+ TH2F* hCopyMe = (TH2F*) gDirectory->Get(histname.c_str());
+ TH2F* hExEg = (TH2F*) hCopyMe->Clone();
+ hExEg->Rebin2D(bing/hExEg->GetXaxis()->GetBinWidth(3),
+ binp/hExEg->GetYaxis()->GetBinWidth(3));
+
+ //Set up canvas
+ TCanvas *cGPSG = new TCanvas("cGPSG","cGPSG",1000,750);
+ cGPSG->Divide(1,2);
+ cGPSG->cd(1);
+
+ //Axis being projected
+ TH1F* projP = Project(hExEg, "Y", kBlack);
+ projP->Draw();
+
+ //Projection gates
+ Prj_DrawLines(Ex1, width1, kRed);
+ Prj_DrawLines(Ex2, width2, kBlue);
+ Prj_DrawLines(Ex3, width3, kGreen);
+
+ //Projected axis
+ cGPSG->cd(2);
+ TH1F* projG1 = Project(hExEg, "X", Ex1, width1, kRed);
+ TH1F* projG2 = Project(hExEg, "X", Ex2, width2, kBlue);
+ TH1F* projG3 = Project(hExEg, "X", Ex3, width3, kGreen);
+ projG1->Draw();
+ projG2->Draw("same");
+ projG3->Draw("same");
+}
+
+//////////////////////////////////////////////////////////////////////////////////////
+void GateParticle_SeeGamma(string addbackOrTrack, double bing, double binp, double Ex1, double width1,
+ double Ex2, double width2,
+ double Ex3, double width3,
+ double Ex4, double width4)
+{
+ string histname = AddOrTrack(addbackOrTrack);
+ //Check if histogram exists alreado, else draw
+ //in future, add check, then try to load, then draw
+ if(gDirectory->Get(histname.c_str()) == 0){
+ DrawPlotsToLoad();
+ }
+
+ //Load in & rebin
+ TH2F* hCopyMe = (TH2F*) gDirectory->Get(histname.c_str());
+ TH2F* hExEg = (TH2F*) hCopyMe->Clone();
+ hExEg->Rebin2D(bing/hExEg->GetXaxis()->GetBinWidth(3),
+ binp/hExEg->GetYaxis()->GetBinWidth(3));
+
+ //Set up canvas
+ TCanvas *cGPSG = new TCanvas("cGPSG","cGPSG",1000,750);
+ cGPSG->Divide(1,2);
+ cGPSG->cd(1);
+
+ //Axis being projected
+ TH1F* projP = Project(hExEg, "Y", kBlack);
+ projP->Draw();
+
+ //Projection gates
+ Prj_DrawLines(Ex1, width1, kRed);
+ Prj_DrawLines(Ex2, width2, kBlue);
+ Prj_DrawLines(Ex3, width3, kGreen);
+ Prj_DrawLines(Ex4, width4, kViolet);
+
+ //Projected axis
+ cGPSG->cd(2);
+ TH1F* projG1 = Project(hExEg, "X", Ex1, width1, kRed);
+ TH1F* projG2 = Project(hExEg, "X", Ex2, width2, kBlue);
+ TH1F* projG3 = Project(hExEg, "X", Ex3, width3, kGreen);
+ TH1F* projG4 = Project(hExEg, "X", Ex4, width4, kViolet);
+ projG1->Draw();
+ projG2->Draw("same");
+ projG3->Draw("same");
+ projG4->Draw("same");
+}
+
+//////////////////////////////////////////////////////////////////////////////////////
+void GateParticle_SeeGamma(string addbackOrTrack, double bing, double binp, double Ex1, double width1,
+ double Ex2, double width2,
+ double Ex3, double width3,
+ double Ex4, double width4,
+ double Ex5, double width5)
+{
+ string histname = AddOrTrack(addbackOrTrack);
+ //Check if histogram exists alreado, else draw
+ //in future, add check, then try to load, then draw
+ if(gDirectory->Get(histname.c_str()) == 0){
+ DrawPlotsToLoad();
+ }
+
+ //Load in & rebin
+ TH2F* hCopyMe = (TH2F*) gDirectory->Get(histname.c_str());
+ TH2F* hExEg = (TH2F*) hCopyMe->Clone();
+ hExEg->Rebin2D(bing/hExEg->GetXaxis()->GetBinWidth(3),
+ binp/hExEg->GetYaxis()->GetBinWidth(3));
+
+ //Set up canvas
+ TCanvas *cGPSG = new TCanvas("cGPSG","cGPSG",1000,750);
+ cGPSG->Divide(1,2);
+ cGPSG->cd(1);
+
+ //Axis being projected
+ TH1F* projP = Project(hExEg, "Y", kBlack);
+ projP->Draw();
+
+ //Projection gates
+ Prj_DrawLines(Ex1, width1, kRed);
+ Prj_DrawLines(Ex2, width2, kBlue);
+ Prj_DrawLines(Ex3, width3, kGreen);
+ Prj_DrawLines(Ex4, width4, kViolet);
+ Prj_DrawLines(Ex5, width5, kOrange);
+
+ //Projected axis
+ cGPSG->cd(2);
+ TH1F* projG1 = Project(hExEg, "X", Ex1, width1, kRed);
+ TH1F* projG2 = Project(hExEg, "X", Ex2, width2, kBlue);
+ TH1F* projG3 = Project(hExEg, "X", Ex3, width3, kGreen);
+ TH1F* projG4 = Project(hExEg, "X", Ex4, width4, kViolet);
+ TH1F* projG5 = Project(hExEg, "X", Ex5, width5, kOrange);
+ projG1->Draw();
+ projG2->Draw("same");
+ projG3->Draw("same");
+ projG4->Draw("same");
+ projG5->Draw("same");
+}
+
+
+
+
+//////////////////////////////////////////////////////////////////////////////////////
+void GateGamma_SeeParticle_HighEnergyStaggered(string addbackOrTrack, double bing, double binp)
+{
+ string histname = AddOrTrack(addbackOrTrack);
+ //Check if histogram exists alreado, else draw
+ //in future, add check, then try to load, then draw
+ if(gDirectory->Get(histname.c_str()) == 0){
+ DrawPlotsToLoad();
+ }
+
+ gStyle->SetPalette(kCool);
+
+ double width = 8.0;
+ //vector<double> energies = {3.5, 4.0, 4.8, 5.0};
+ //vector<double> energies = {3.5,3.6,3.7,3.8,3.9,4.0,4.1,4.2,4.3,4.4,4.5,4.6,4.7,4.8,4.9,5.0};
+ vector<double> energies = {3.4,3.6,3.8,4.0,4.8};
+ auto cols = TColor::GetPalette();
+
+ //Load in & rebin
+ TH2F* hCopyMe = (TH2F*) gDirectory->Get(histname.c_str());
+ TH2F* hExEg = (TH2F*) hCopyMe->Clone();
+ hExEg->Rebin2D(bing/hExEg->GetXaxis()->GetBinWidth(3),
+ binp/hExEg->GetYaxis()->GetBinWidth(3));
+
+ //Set up canvas
+ TCanvas *cGPSG_High = new TCanvas("cGPSG_High","cGPSG_High",1000,750);
+ cGPSG_High->Divide(1,2);
+ cGPSG_High->cd(1);
+
+ //Axis being projected
+ TH1F* projG = Project(hExEg, "X", kBlack);
+ projG->Draw();
+
+ //Projection gatesi
+ for(int i=0; i<energies.size(); i++){
+ Prj_DrawLines(energies.at(i)+width/2., width, cols.At(i));
+ }
+
+ //Projected axis
+ cGPSG_High->cd(2);
+ TH1F* projP1 = Project(hExEg, "Y", energies.at(0)+width/2., width, cols.At(0));
+ projP1->Draw("plc");
+ for(int i=1; i<energies.size(); i++){
+ TH1F* projP = Project(hExEg, "Y", energies.at(i)+width/2., width, cols.At(i));
+ projP->Draw("plc same");
+ }
+
+}
+
+
+
+
+//=====================================================================//
+//=====================================================================//
+
+////////////////////////////////////////////////////////////////////////////////
+void DrawExThetaLab_SubPSp(TChain* tree)
+{
+ LoadChain_PhaseSpace();
+ TCanvas* c = new TCanvas("c","c",500,500);
+
+ // ExTL Good
+ tree->Draw("Ex:ThetaLab>>hGood(26,104,156,150,-2,13)","abs(T_FPMW_HF-13000)<1000 && Mugast.TelescopeNumber<10","colz");
+ TH2F* hGood = (TH2F*) gDirectory->Get("hGood");
+
+ // ExTL Bad
+ tree->Draw("Ex:ThetaLab>>hBad(26,104,156,150,-2,13)","abs(T_FPMW_HF-10200)<600 && Mugast.TelescopeNumber<10","colz");
+ TH2F* hBad = (TH2F*) gDirectory->Get("hBad");
+
+ // ExTL Clean
+ TH2F* hClean = (TH2F*) hGood->Clone();
+ hClean->SetTitle("hClean");
+ hClean->SetName("hClean");
+ hClean->Add(hBad,-1);
+
+ // ExTL PSp
+ tree_PSp00->Draw("Ex:ThetaLab>>hPSp(26,104,156,150,-2,13)","","colz");
+ TH2F* hPSp = (TH2F*) gDirectory->Get("hPSp");
+
+ double scale = 1.0;
+ int minX = hGood->GetXaxis()->FindBin(106);//110);
+ int maxX = hGood->GetXaxis()->FindBin(130);//145);
+ int minY = hGood->GetYaxis()->FindBin(7.0);
+ int maxY = hGood->GetYaxis()->FindBin(9.5);
+ int step = (maxX-minX)/3;
+
+ cout << minX << " " << maxX << " || " << minY << " " << maxY << endl;
+
+ for(int x = minX; x < maxX; x++){
+// cout << " ---- X = " << x << " ---- " << endl;
+ for(int y = minY; y < maxY; y++){
+// cout << " ---- Y = " << y << " ---- " << endl;
+
+ if(hClean->GetYaxis()->GetBinCenter(y) < (-0.0833 * hClean->GetXaxis()->GetBinCenter(x) ) + 17.75){
+// cout << x << " " << y << " || " << hPSp->GetBinContent(x,y) << " vs " << hClean->GetBinContent(x,y) << endl;
+
+ if(hClean->GetBinContent(x,y)>2){
+ while(hPSp->GetBinContent(x,y) * scale > hClean->GetBinContent(x,y)){
+ scale *= 0.999;
+ cout << hPSp->GetXaxis()->GetBinCenter(x) << ", " << hPSp->GetYaxis()->GetBinCenter(y)
+ << " | " << hPSp->GetBinContent(x,y) * scale << " > " << hClean->GetBinContent(x,y)
+ << " | " << scale << endl;
+ }
+ }
+
+ }
+ }
+ }
+
+ cout << " |||| scale = " << scale << " |||| " << endl;
+ delete c;
+
+ TCanvas* cClean = new TCanvas("cClean","cClean",500,500);
+ hClean->Draw("colz");
+
+ TCanvas* cPSp = new TCanvas("cPSp","cPSp",500,500);
+ hPSp->Draw("colz");
+
+ TCanvas* cFinal = new TCanvas("cFinal","cFinal",500,500);
+ TH2F* hFinal = (TH2F*) hClean->Clone();
+ hFinal->SetTitle("hFinal");
+ hFinal->SetName("hFinal");
+ hFinal->Add(hPSp,-scale);
+ hFinal->Draw("colz");
+
+ TCanvas* cProj = new TCanvas("cProj","cProj",1000,500);
+ TH1F* hLowC = (TH1F*) hClean->ProjectionY("LC", minX, minX+step);
+ TH1F* hMidC = (TH1F*) hClean->ProjectionY("MC", minX+step, maxX-step);
+ TH1F* hHghC = (TH1F*) hClean->ProjectionY("HC", maxX-step, maxX );
+
+ hLowC->SetLineColor(kBlack);
+ hMidC->SetLineColor(kBlack);
+ hHghC->SetLineColor(kBlack);
+
+ TH1F* hLowP = (TH1F*) hPSp->ProjectionY("LP", minX, minX+step);
+ TH1F* hMidP = (TH1F*) hPSp->ProjectionY("MP", minX+step, maxX-step);
+ TH1F* hHghP = (TH1F*) hPSp->ProjectionY("HP", maxX-step, maxX );
+
+ hLowP->SetLineColor(kRed); hLowP->Scale(scale);
+ hMidP->SetLineColor(kRed); hMidP->Scale(scale);
+ hHghP->SetLineColor(kRed); hHghP->Scale(scale);
+
+ TH1F* hLowF = (TH1F*) hFinal->ProjectionY("LF", minX, minX+step);
+ TH1F* hMidF = (TH1F*) hFinal->ProjectionY("MF", minX+step, maxX-step);
+ TH1F* hHghF = (TH1F*) hFinal->ProjectionY("HF", maxX-step, maxX );
+
+ hLowF->SetLineColor(kBlue);
+ hMidF->SetLineColor(kBlue);
+ hHghF->SetLineColor(kBlue);
+
+ cProj->Divide(3);
+ cProj->cd(1); hLowC->Draw(); hLowP->Draw("same"); hLowF->Draw("same");
+ cProj->cd(2); hMidC->Draw(); hMidP->Draw("same"); hMidF->Draw("same");
+ cProj->cd(3); hHghC->Draw(); hHghP->Draw("same"); hHghF->Draw("same");
+
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void DrawEx_SubPSp(TChain* tree)
+{
+ LoadChain_PhaseSpace();
+ TCanvas* c = new TCanvas("c","c",1000,1000);
+
+ // ExTL Good
+ tree->Draw("Ex>>hGood(300,-2,13)","abs(T_FPMW_HF-13000)<1000 && Mugast.TelescopeNumber<10","colz");
+ TH1F* hGood = (TH1F*) gDirectory->Get("hGood");
+
+ // ExTL Bad
+ tree->Draw("Ex>>hBad(300,-2,13)","abs(T_FPMW_HF-10200)<600 && Mugast.TelescopeNumber<10","colz");
+ TH1F* hBad = (TH1F*) gDirectory->Get("hBad");
+
+ // ExTL Clean
+ TH1F* hClean = (TH1F*) hGood->Clone();
+ hClean->SetTitle("hClean");
+ hClean->SetName("hClean");
+ hClean->Add(hBad,-1);
+
+ // ExTL PSp
+ tree_PSp00->Draw("Ex>>hPSp(300,-2,13)","","colz");
+ TH1F* hPSp = (TH1F*) gDirectory->Get("hPSp");
+
+ double scale = 0.134256;
+ cout << " || SUBTRACTING PHASE SPACE WITH SCALE = " << scale << endl;
+
+ TH1F* hFinal = (TH1F*) hClean->Clone();
+ hFinal->SetTitle("hFinal");
+ hFinal->SetName("hFinal");
+ hFinal->Add(hPSp,-scale);
+
+ hClean->SetLineColor(kBlack); hClean->Draw();
+ hPSp ->SetLineColor(kRed); hPSp->Scale(scale); hPSp ->Draw("same");
+ hFinal->SetLineColor(kBlue); hFinal->Draw("same");
+
+}
+
+/////////////////////////////////////////////////////////////////////////////////////
+void AddKnownStateLines(TH1F* hist, double ymax)
+{
+ hist->Draw();
+
+ TLine* Sn = new TLine(7.608, 0.0,7.608, ymax);
+ Sn->SetLineColor(kOrange);
+ Sn->SetLineWidth(4);
+ Sn->SetLineStyle(kDashed);
+ Sn->Draw();
+
+ for(int i=0; i<numPeaks; i++){
+ TLine* l = new TLine(means.at(i), 0.0, means.at(i), ymax);
+ l->SetLineColor(kBlack);
+ l->SetLineWidth(2);
+ l->SetLineStyle(kDotted);
+ l->Draw();
+ }
+}
+
+void DrawGammas_AddBackVersusTracked()
+{
+ TCanvas* cGammaCompare = new TCanvas("cGammaCompare","cGammaCompare",1600,800);
+
+ tree_both->Draw("trackEDC>>hTracked(1000,-1,9)",
+ "Mugast.TelescopeNumber<9 && abs(T_FPMW_HF-13000)<1000",
+ "");
+ TH1F* hTracked = (TH1F*) gDirectory->Get("hTracked");
+ hTracked->SetLineColor(kRed);
+ hTracked->SetLineWidth(2);
+
+
+ tree_both->Draw("AddBack_EDC>>hAddBack(1000,-1,9)",
+ "Mugast.TelescopeNumber<9 && abs(T_FPMW_HF-13000)<1000",
+ "");
+ TH1F* hAddBack = (TH1F*) gDirectory->Get("hAddBack");
+ hAddBack->SetLineColor(kBlue);
+ hAddBack->SetLineWidth(2);
+
+ hTracked->Draw(); hAddBack->Draw("same");
+
+}
+
+
+//=====================================================================//
+//=====================================================================//
+
+////////////////////////////////////////////////////////////////////////////////
+void Diagnostics_StripNumber(TTree* tree)
+{
+ TCanvas* c = new TCanvas("c","c",1000,1000);
+
+ c->Divide(4,4);
+
+ c->cd(1);
+ tree->Draw("Mugast.DSSD_X>>hX1(130,0,130)","Mugast.TelescopeNumber==1","");
+ c->cd(2);
+ tree->Draw("Mugast.DSSD_Y>>hY1(130,0,130)","Mugast.TelescopeNumber==1","");
+
+ c->cd(3);
+ tree->Draw("Mugast.DSSD_X>>hX5(130,0,130)","Mugast.TelescopeNumber==5","");
+ c->cd(4);
+ tree->Draw("Mugast.DSSD_Y>>hY5(130,0,130)","Mugast.TelescopeNumber==5","");
+
+ c->cd(5);
+ tree->Draw("Mugast.DSSD_X>>hX2(130,0,130)","Mugast.TelescopeNumber==2","");
+ c->cd(6);
+ tree->Draw("Mugast.DSSD_Y>>hY2(130,0,130)","Mugast.TelescopeNumber==2","");
+
+ c->cd(7);
+ tree->Draw("Mugast.DSSD_X>>hX6(130,0,130)","Mugast.TelescopeNumber==6","");
+ c->cd(8);
+ tree->Draw("Mugast.DSSD_Y>>hY6(130,0,130)","Mugast.TelescopeNumber==6","");
+
+ c->cd(9);
+ tree->Draw("Mugast.DSSD_X>>hX3(130,0,130)","Mugast.TelescopeNumber==3","");
+ c->cd(10);
+ tree->Draw("Mugast.DSSD_Y>>hY3(130,0,130)","Mugast.TelescopeNumber==3","");
+
+ c->cd(11);
+ tree->Draw("Mugast.DSSD_X>>hX7(130,0,130)","Mugast.TelescopeNumber==7","");
+ c->cd(12);
+ tree->Draw("Mugast.DSSD_Y>>hY7(130,0,130)","Mugast.TelescopeNumber==7","");
+
+ c->cd(13);
+ tree->Draw("Mugast.DSSD_X>>hX4(130,0,130)","Mugast.TelescopeNumber==4","");
+ c->cd(14);
+ tree->Draw("Mugast.DSSD_Y>>hY4(130,0,130)","Mugast.TelescopeNumber==4","");
+
+
+
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void GateThetaLab_MultiWrite(int numGates, double binsize)
+{
+ double startTheta= 104.;
+ double finishTheta=156.;
+
+ string drawbase = to_string(numGates)
+ + ",104,156,"
+ + to_string(20.0/binsize)
+ + ",-5,15)";
+ string drawgood = "Ex:ThetaLab>>hGood(" + drawbase;
+ string drawbad = "Ex:ThetaLab>>hBad(" + drawbase;
+
+ tree_both->Draw(drawgood.c_str(),"abs(T_FPMW_HF-13000)<1000","colz");
+ TH2F* hGood = (TH2F*) gDirectory->Get("hGood");
+
+ tree_both->Draw(drawbad.c_str(), "abs(T_FPMW_HF-10200)<600","colz");
+ TH2F* hBad = (TH2F*) gDirectory->Get("hBad");
+
+ TH2F* hClean = (TH2F*) hGood->Clone();
+ hClean->SetTitle("hClean");
+ hClean->SetName("hClean");
+ hClean->Add(hBad,-1);
+ hClean->Draw("colz");
+
+ string ytitle = "Counts / " + to_string(binsize) + " MeV";
+ double gatesize = (finishTheta-startTheta)/numGates;
+ TList* list = new TList();
+
+ for (int i=0; i<numGates; i++){
+ cout << "Writing gate " << i+1 << "/" << numGates << endl;
+ double minTheta = startTheta + (i * gatesize);
+ string title = "TLab"+to_string((int) minTheta)+"to"+to_string((int) (minTheta+gatesize));
+
+ TH1F* Ex_ThetaLabGate = (TH1F*) hClean->ProjectionY(title.c_str(), i, i+1, "");
+ Ex_ThetaLabGate->GetXaxis()->SetTitle("Ex [MeV]");
+ Ex_ThetaLabGate->GetYaxis()->SetTitle(ytitle.c_str());
+ Ex_ThetaLabGate->Sumw2();
+ Ex_ThetaLabGate->SetTitle(title.c_str());
+ list->Add(Ex_ThetaLabGate);
+ }
+
+ TFile* file = new TFile("./C2S_Outputs/GateThetaLabHist_Experiment.root","RECREATE");
+ list->Write("GateThetaLabHistograms",TObject::kSingleKey);
+ file->ls();
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void GatePhaseSpace_MultiWrite(int numGates, double binsize)
+{
+ double startTheta= 104.;
+ double finishTheta=156.;
+
+ string draw = "Ex:ThetaLab>>h("
+ + to_string(numGates)
+ + ",104,156,"
+ + to_string(20.0/binsize)
+ + ",-5,15)";
+
+ tree_PSp00->Draw(draw.c_str(),"","colz");
+ TH2F* h = (TH2F*) gDirectory->Get("h");
+ h->Draw("colz");
+
+ string ytitle = "Counts / " + to_string(binsize) + " MeV";
+ double gatesize = (finishTheta-startTheta)/numGates;
+ TList* list = new TList();
+
+ for (int i=0; i<numGates; i++){
+ cout << "Writing gate " << i+1 << "/" << numGates << endl;
+ double minTheta = startTheta + (i * gatesize);
+ string title = "TLab"+to_string((int) minTheta)+"to"+to_string((int) (minTheta+gatesize));
+
+ TH1F* Ex_ThetaLabGate = (TH1F*) h->ProjectionY(title.c_str(), i, i+1, "");
+ Ex_ThetaLabGate->GetXaxis()->SetTitle("Ex [MeV]");
+ Ex_ThetaLabGate->GetYaxis()->SetTitle(ytitle.c_str());
+ Ex_ThetaLabGate->Sumw2();
+ Ex_ThetaLabGate->SetTitle(title.c_str());
+ list->Add(Ex_ThetaLabGate);
+ }
+
+ TFile* file = new TFile("./C2S_Outputs/GateThetaLabHist_PhaseSpace.root","RECREATE");
+ list->Write("GateThetaLabHistograms",TObject::kSingleKey);
+ file->ls();
+}
+
+//=====================================================================//
+//=====================================================================//
+
+void TWOFNR_ApproachingUnbound()
+{
+ TCanvas* cTWOFNR = new TCanvas("cTWOFNR","cTWOFNR",1000,1000);
+ cTWOFNR->Divide(2,2);
+ cTWOFNR->cd(1)->SetLogy();
+ cTWOFNR->cd(2)->SetLogy();
+ cTWOFNR->cd(3)->SetLogy();
+ cTWOFNR->cd(4)->SetLogy();
+
+ cTWOFNR->cd(1);
+ TGraph* s12_7500 = TWOFNR(7.500,2,1,0,0.5); s12_7500->SetLineColor(colourR->GetNumber()); s12_7500->SetLineWidth(4); s12_7500->SetLineStyle(1); s12_7500->GetXaxis()->SetRangeUser(100,160); s12_7500->GetYaxis()->SetRangeUser(0.004,2); s12_7500->Draw();
+ TGraph* s12_7520 = TWOFNR(7.520,2,1,0,0.5); s12_7520->SetLineColor(colourB->GetNumber()); s12_7520->SetLineWidth(2); s12_7520->SetLineStyle(7); s12_7520->GetXaxis()->SetRangeUser(100,160); s12_7520->GetYaxis()->SetRangeUser(0.004,2); s12_7520->Draw("same");
+ TGraph* s12_7540 = TWOFNR(7.540,2,1,0,0.5); s12_7540->SetLineColor(colourB->GetNumber()); s12_7540->SetLineWidth(2); s12_7540->SetLineStyle(7); s12_7540->GetXaxis()->SetRangeUser(100,160); s12_7540->GetYaxis()->SetRangeUser(0.004,2); s12_7540->Draw("same");
+ TGraph* s12_7560 = TWOFNR(7.560,2,1,0,0.5); s12_7560->SetLineColor(colourB->GetNumber()); s12_7560->SetLineWidth(2); s12_7560->SetLineStyle(7); s12_7560->GetXaxis()->SetRangeUser(100,160); s12_7560->GetYaxis()->SetRangeUser(0.004,2); s12_7560->Draw("same");
+ TGraph* s12_7580 = TWOFNR(7.580,2,1,0,0.5); s12_7580->SetLineColor(colourB->GetNumber()); s12_7580->SetLineWidth(2); s12_7580->SetLineStyle(7); s12_7580->GetXaxis()->SetRangeUser(100,160); s12_7580->GetYaxis()->SetRangeUser(0.004,2); s12_7580->Draw("same");
+ TGraph* s12_7600 = TWOFNR(7.600,2,1,0,0.5); s12_7600->SetLineColor(colourG->GetNumber()); s12_7600->SetLineWidth(4); s12_7600->SetLineStyle(1); s12_7600->GetXaxis()->SetRangeUser(100,160); s12_7600->GetYaxis()->SetRangeUser(0.004,2); s12_7600->Draw("same");
+
+ cTWOFNR->cd(2);
+ TGraph* p12_7500 = TWOFNR(7.500,2,0,1,0.5); p12_7500->SetLineColor(colourR->GetNumber()); p12_7500->SetLineWidth(4); p12_7500->SetLineStyle(1); p12_7500->GetXaxis()->SetRangeUser(100,160); p12_7500->GetYaxis()->SetRangeUser(0.004,2); p12_7500->Draw();
+ TGraph* p12_7520 = TWOFNR(7.520,2,0,1,0.5); p12_7520->SetLineColor(colourB->GetNumber()); p12_7520->SetLineWidth(2); p12_7520->SetLineStyle(7); p12_7520->GetXaxis()->SetRangeUser(100,160); p12_7520->GetYaxis()->SetRangeUser(0.004,2); p12_7520->Draw("same");
+ TGraph* p12_7540 = TWOFNR(7.540,2,0,1,0.5); p12_7540->SetLineColor(colourB->GetNumber()); p12_7540->SetLineWidth(2); p12_7540->SetLineStyle(7); p12_7540->GetXaxis()->SetRangeUser(100,160); p12_7540->GetYaxis()->SetRangeUser(0.004,2); p12_7540->Draw("same");
+ TGraph* p12_7560 = TWOFNR(7.560,2,0,1,0.5); p12_7560->SetLineColor(colourB->GetNumber()); p12_7560->SetLineWidth(2); p12_7560->SetLineStyle(7); p12_7560->GetXaxis()->SetRangeUser(100,160); p12_7560->GetYaxis()->SetRangeUser(0.004,2); p12_7560->Draw("same");
+ TGraph* p12_7580 = TWOFNR(7.580,2,0,1,0.5); p12_7580->SetLineColor(colourB->GetNumber()); p12_7580->SetLineWidth(2); p12_7580->SetLineStyle(7); p12_7580->GetXaxis()->SetRangeUser(100,160); p12_7580->GetYaxis()->SetRangeUser(0.004,2); p12_7580->Draw("same");
+ TGraph* p12_7600 = TWOFNR(7.600,2,0,1,0.5); p12_7600->SetLineColor(colourG->GetNumber()); p12_7600->SetLineWidth(4); p12_7600->SetLineStyle(1); p12_7600->GetXaxis()->SetRangeUser(100,160); p12_7600->GetYaxis()->SetRangeUser(0.004,2); p12_7600->Draw("same");
+
+ cTWOFNR->cd(3);
+ TGraph* d52_7500 = TWOFNR(7.500,2,0,2,2.5); d52_7500->SetLineColor(colourR->GetNumber()); d52_7500->SetLineWidth(4); d52_7500->SetLineStyle(1); d52_7500->GetXaxis()->SetRangeUser(100,160); d52_7500->GetYaxis()->SetRangeUser(0.004,2); d52_7500->Draw();
+ TGraph* d52_7520 = TWOFNR(7.520,2,0,2,2.5); d52_7520->SetLineColor(colourB->GetNumber()); d52_7520->SetLineWidth(2); d52_7520->SetLineStyle(7); d52_7520->GetXaxis()->SetRangeUser(100,160); d52_7520->GetYaxis()->SetRangeUser(0.004,2); d52_7520->Draw("same");
+ TGraph* d52_7540 = TWOFNR(7.540,2,0,2,2.5); d52_7540->SetLineColor(colourB->GetNumber()); d52_7540->SetLineWidth(2); d52_7540->SetLineStyle(7); d52_7540->GetXaxis()->SetRangeUser(100,160); d52_7540->GetYaxis()->SetRangeUser(0.004,2); d52_7540->Draw("same");
+ TGraph* d52_7560 = TWOFNR(7.560,2,0,2,2.5); d52_7560->SetLineColor(colourB->GetNumber()); d52_7560->SetLineWidth(2); d52_7560->SetLineStyle(7); d52_7560->GetXaxis()->SetRangeUser(100,160); d52_7560->GetYaxis()->SetRangeUser(0.004,2); d52_7560->Draw("same");
+ TGraph* d52_7580 = TWOFNR(7.580,2,0,2,2.5); d52_7580->SetLineColor(colourB->GetNumber()); d52_7580->SetLineWidth(2); d52_7580->SetLineStyle(7); d52_7580->GetXaxis()->SetRangeUser(100,160); d52_7580->GetYaxis()->SetRangeUser(0.004,2); d52_7580->Draw("same");
+ TGraph* d52_7600 = TWOFNR(7.600,2,0,2,2.5); d52_7600->SetLineColor(colourG->GetNumber()); d52_7600->SetLineWidth(4); d52_7600->SetLineStyle(1); d52_7600->GetXaxis()->SetRangeUser(100,160); d52_7600->GetYaxis()->SetRangeUser(0.004,2); d52_7600->Draw("same");
+
+ cTWOFNR->cd(4);
+ TGraph* f72_7500 = TWOFNR(7.500,3,0,3,3.5); f72_7500->SetLineColor(colourR->GetNumber()); f72_7500->SetLineWidth(4); f72_7500->SetLineStyle(1); f72_7500->GetXaxis()->SetRangeUser(100,160); f72_7500->GetYaxis()->SetRangeUser(0.004,2); f72_7500->Draw();
+ TGraph* f72_7520 = TWOFNR(7.520,3,0,3,3.5); f72_7520->SetLineColor(colourB->GetNumber()); f72_7520->SetLineWidth(2); f72_7520->SetLineStyle(7); f72_7520->GetXaxis()->SetRangeUser(100,160); f72_7520->GetYaxis()->SetRangeUser(0.004,2); f72_7520->Draw("same");
+ TGraph* f72_7540 = TWOFNR(7.540,3,0,3,3.5); f72_7540->SetLineColor(colourB->GetNumber()); f72_7540->SetLineWidth(2); f72_7540->SetLineStyle(7); f72_7540->GetXaxis()->SetRangeUser(100,160); f72_7540->GetYaxis()->SetRangeUser(0.004,2); f72_7540->Draw("same");
+ TGraph* f72_7560 = TWOFNR(7.560,3,0,3,3.5); f72_7560->SetLineColor(colourB->GetNumber()); f72_7560->SetLineWidth(2); f72_7560->SetLineStyle(7); f72_7560->GetXaxis()->SetRangeUser(100,160); f72_7560->GetYaxis()->SetRangeUser(0.004,2); f72_7560->Draw("same");
+ TGraph* f72_7580 = TWOFNR(7.580,3,0,3,3.5); f72_7580->SetLineColor(colourB->GetNumber()); f72_7580->SetLineWidth(2); f72_7580->SetLineStyle(7); f72_7580->GetXaxis()->SetRangeUser(100,160); f72_7580->GetYaxis()->SetRangeUser(0.004,2); f72_7580->Draw("same");
+ TGraph* f72_7600 = TWOFNR(7.600,3,0,3,3.5); f72_7600->SetLineColor(colourG->GetNumber()); f72_7600->SetLineWidth(4); f72_7600->SetLineStyle(1); f72_7600->GetXaxis()->SetRangeUser(100,160); f72_7600->GetYaxis()->SetRangeUser(0.004,2); f72_7600->Draw("same");
+
+
+
+
+
+
+
+
+
+}
+
diff --git a/Projects/e775s/analysis/Plots_KnownPeakFitter.h b/Projects/e775s/analysis/Plots_KnownPeakFitter.h
new file mode 100644
index 0000000000000000000000000000000000000000..4ba99d84c44bcc9ef6cc341a1d5a9aea8882945e
--- /dev/null
+++ b/Projects/e775s/analysis/Plots_KnownPeakFitter.h
@@ -0,0 +1,138 @@
+#include "TMath.h"
+#include "math.h"
+#include <cmath>
+#include "stdlib.h"
+
+double pi = M_PI;
+
+vector<double> means = { 0.000,
+ 1.675,
+ 3.572,
+ 4.071,
+ 5.004,
+ 5.228,
+ 5.629,
+ 7.622,
+ 7.754,
+ 8.554,
+ 8.96,
+ 9.77,
+ 10.13,
+ };
+
+const int numPeaks = 13;//7;//means.size();
+double meanCorrection_m = 1.0;//0.975;
+double meanCorrection_c = 0.0;//0.090;
+double sgmaCorrection_m = 1.0;//-0.00557;
+double sgmaCorrection_c = 0.0;//+0.123;
+
+vector<vector<double>> FitKnownPeaks_RtrnArray(TH1F* hist){
+ //double minfit = -1.0, maxfit = 7.0;
+ double minfit = -1.0, maxfit = 10.0;
+ double binwidth = hist->GetXaxis()->GetBinWidth(5);
+ double sigma = 0.11;
+
+ // Build individual peak functions
+ string nameBase = "Peak";
+ string function = "gaus";
+ TF1 **allPeaks = new TF1*[numPeaks];
+ for(int i=0; i<numPeaks; i++){
+ string nameHere = nameBase;
+ nameHere += to_string(i);
+ allPeaks[i] = new TF1(nameHere.c_str(), function.c_str(), minfit, maxfit);
+ allPeaks[i]->SetLineColor(4001);
+ allPeaks[i]->SetLineStyle(7);
+ allPeaks[i]->SetNpx(500);
+ }
+
+ // Build background function (SIMPLE FOR NOW!)
+ TF1 *bg = new TF1("bg","[0]",minfit, maxfit);
+ bg->SetLineColor(kGreen);
+ bg->SetLineStyle(9);
+ bg->SetParNames("Background");
+
+ // Build total function
+ int bgRefNum = ((numPeaks-1)*3)+3;
+ string s_full = "gaus(0)";
+ for(int i=1; i<numPeaks; i++){
+ s_full += "+ gaus("+ to_string(3*i) + ")";
+ }
+ TF1 *full = new TF1("fitAllPeaks",s_full.c_str(),minfit, maxfit);
+ full->SetLineColor(kRed);
+ full->SetNpx(500);
+
+ // Set and fix some parameters
+ for(int i=0; i<numPeaks; i++){
+ full->FixParameter((i*3)+1, meanCorrection_m * means.at(i) + meanCorrection_c);
+ //full->FixParameter((i*3)+2, sigma);
+ //full->FixParameter((i*3)+2, sgmaCorrection_m * means.at(i) + sgmaCorrection_c);
+ full->SetParameter((i*3)+2, sigma);
+ full->SetParLimits((i*3)+2, 0.95*sigma, 1.05*sigma);
+ full->SetParameter((i*3)+0, 1e1);
+ full->SetParLimits((i*3)+0, 0.0,1e8);
+ }
+
+ // Fit full function
+ hist->Fit(full, "RWQB", "", minfit, maxfit);
+
+ // Extract parameters
+ const Double_t* finalPar = full->GetParameters();
+ const Double_t* finalErr = full->GetParErrors();
+ for(int i=0; i<numPeaks; i++){
+ allPeaks[i]->SetParameters(finalPar[0+(i*3)],
+ finalPar[1+(i*3)],
+ finalPar[2+(i*3)]);
+ allPeaks[i]->SetLineColor(kOrange);
+ }
+
+ //Draw all peaks
+ hist->Draw();
+ for(int i=0; i<numPeaks; i++){
+ allPeaks[i]->Draw("same");
+ }
+
+ //Write to screen
+ cout << "===========================" << endl;
+ cout << "== PEAK =========== AREA ==" << endl;
+
+ //Loop over every peak
+ vector<vector<double>> allpeaks;
+ for(int i=0; i<numPeaks; i++){
+ //For AREA = HEIGHT * SIGMA * SQRT(2*PI)
+ double A = (finalPar[0+(i*3)] * finalPar[2+(i*3)] * sqrt(2*pi)) / binwidth;
+
+ double deltaA = A * sqrt( pow( finalErr[0+(i*3)]/finalPar[0+(i*3)] ,2) + pow( finalErr[2+(i*3)]/finalPar[2+(i*3)] ,2) );
+ cout << "DELTAAREA = "
+ << finalErr[0+(i*3)] << " / "
+ << finalPar[0+(i*3)] << " and "
+ << finalErr[2+(i*3)] << " / "
+ << finalPar[2+(i*3)] << " = "
+ << deltaA << endl;
+
+ cout << fixed << setprecision(3)
+ << " #" << i << " "
+ << finalPar[(i*3)+1] << "\t" << setprecision(0)
+ << A << "\t+- "
+ << deltaA << setprecision(3);
+ cout << " SQRT: " << sqrt(A) << endl;
+
+ vector<double> onepeak; //energy, area and error for one peak
+ onepeak.push_back(finalPar[(i*3)+1]);
+ onepeak.push_back(A);
+ onepeak.push_back(deltaA);
+ allpeaks.push_back(onepeak);
+ cout << "-------------" << endl;
+ }
+
+ //if(removing){
+ // for(int b=1; b<hist->GetNbinsX()-1; b++){
+ // cout << b << " " << hist->GetBinContent(b) << " " << full->Eval(hist->GetBinCenter(b)) << endl;
+ // fitContByBin.push_back(full->Eval(hist->GetBinCenter(b)));
+ // }
+ //}
+
+ return allpeaks;
+}
+
+int FitKnownPeaks(TH1F* hist){auto dump = FitKnownPeaks_RtrnArray(hist); return 0; }
+
diff --git a/Projects/e775s/analysis/Plots_ParryFunctions.h b/Projects/e775s/analysis/Plots_ParryFunctions.h
new file mode 100644
index 0000000000000000000000000000000000000000..bbcbde7cedbf239f028a7e1bfd463c58e4bd4323
--- /dev/null
+++ b/Projects/e775s/analysis/Plots_ParryFunctions.h
@@ -0,0 +1,285 @@
+vector<TH1F*> fitEx;
+
+vector<int> peaksToFit{7622, 7855, 8313, 8554, 8962, 9200, 9400};
+//vector<int> peaksToFit{7622, 7855, 8313, 8554, 9200, 9400};
+//vector<int> peaksToFit{7622, 7855, 8554, 8962, 9400};
+//vector<int> peaksToFit{7622, 7855, 8554, 9200, 9400};
+
+string basePath = "~/Programs/nptoolv3/Outputs/Analysis/Sim_19Odp_";
+string gateBase = "Mugast.TelescopeNumber<9";
+//string gateBase = "Mugast.TelescopeNumber<9 && ELab<(+0.05458*ThetaLab-3.651) && ELab>(-0.05458*ThetaLab+7.701) && ELab>(+0.05458*ThetaLab-5.501)"; // THIS IS NOT WORKING BECAUSE THE PHASE SPACE IS OVERSUBTRACTED!!
+//string gateBase = "Mugast.TelescopeNumber<9 && ELab<(+0.05458*ThetaLab-4.301) && ELab>(-0.05458*ThetaLab+7.051) && ELab>(+0.05458*ThetaLab-4.851)"; // Ex limit = 10 MeV
+string gateG = gateBase + " && abs(T_FPMW_HF-13000)<1000";
+string gateB = gateBase + " && abs(T_FPMW_HF-10200)< 600";
+
+//int angleMin = 104; int angleMax = 134; double angleBinWidth = 10;
+//int angleMin = 104; int angleMax = 144; double angleBinWidth = 5; double exMin = 6; double exMax = 11; double exBinWidth = 0.05;
+int angleMin = 104; int angleMax = 134; double angleBinWidth = 3; double exMin = 6; double exMax = 11; double exBinWidth = 0.05;
+int angleBinNumb = (angleMax-angleMin)/angleBinWidth; int exBinNumb = (exMax-exMin)/exBinWidth;
+string drawstring = "Ex:ThetaLab>>h(" + to_string(angleBinNumb) + ","
+ + to_string(angleMin) + ","
+ + to_string(angleMax) + ","
+ + to_string(exBinNumb) + ","
+ + to_string(exMin) + ","
+ + to_string(exMax) + ")";
+
+vector<TGraphErrors*> parryPeakAoSA; vector<TGraphErrors*> parryPeakdSdO;
+
+TH2F* Parry_DrawPlots_ExpSubPSp(TChain* treeExp){
+ LoadChain_PhaseSpace();
+ TChain* treePSp = tree_PSp00;
+
+ //Prepare
+ gStyle->SetOptStat(0);
+ TCanvas* cExp = new TCanvas("cExp","cExp",500,500);
+
+ //GOOD EXP
+ treeExp->Draw(drawstring.c_str(),gateG.c_str(),"");
+ TH2F* hG = (TH2F*) gDirectory->Get("h"); hG->SetTitle("hG"); hG->SetName("hG");
+
+ //BAD EXP
+ treeExp->Draw(drawstring.c_str(),gateB.c_str(),"");
+ TH2F* hB = (TH2F*) gDirectory->Get("h"); hB->SetTitle("hB"); hB->SetName("hB");
+
+ //CLEAN EXP
+ TH2F* hC = (TH2F*) hG->Clone(); hC->SetTitle("hC"); hC->SetName("hC");
+ hC->Add(hB,-1);
+ cExp->cd(); hC->Draw("colz");
+
+ //PHASE SPACE
+ TCanvas* cPSp = new TCanvas("cPSp","cPSp",500,500);
+ treePSp->Draw(drawstring.c_str(),gateBase.c_str(),"");
+ TH2F* hP = (TH2F*) gDirectory->Get("h"); hP->SetTitle("hP"); hP->SetName("hP");
+ cPSp->cd(); hP->Draw("colz");
+
+ //SCALE PHASE SPACE
+ double scale = 1.0;
+// //for(int x=0; x<hC->GetNbinsX(); x++){
+// for(int x=hP->GetXaxis()->FindBin(110); x<hP->GetXaxis()->FindBin(122); x++){
+// //for(int y=0; y<hC->GetNbinsY(); y++){
+// for(int y=hP->GetYaxis()->FindBin(7.6); y<hP->GetYaxis()->FindBin(9.7); y++){
+// if(hP->GetYaxis()->GetBinCenter(y) < (-0.08*hP->GetXaxis()->GetBinCenter(x))+18.75){
+//cout << hP->GetXaxis()->GetBinCenter(x) << " " << hP->GetYaxis()->GetBinCenter(y) << endl;
+// //if(hP->GetBinContent(x,y) > 0.002*hP->Integral() && hC->GetBinContent(x,y) > 10){
+// //if(hP->GetBinContent(x,y) > 0.005*hP->Integral() && hC->GetBinContent(x,y) > 30){
+// //if(hP->GetBinContent(x,y) > 0.0010*hP->Integral() && hC->GetBinContent(x,y) > 4){
+// //if(hP->GetBinContent(x,y) > 0.0010*hP->Integral() && hC->GetBinContent(x,y) > 2){
+// while((hP->GetBinContent(x,y))*scale > hC->GetBinContent(x,y)){
+// scale *= 0.999;
+////cout << " scaling... " << endl;
+// }
+// }
+// }
+// }
+// cout << "SCALE = " << scale << " !!! SHOULD BE AROUND 0.14!!!!" << endl;
+ scale = 0.134256; cout << GREEN << " USING SCALING FROM DrawEx_SubPSp(tree_both), scale = " << scale << RESET << endl;
+
+ //CHECK SCALE PHASE SPACE
+ TCanvas* cCheck = new TCanvas("cCheck","cCheck",500,500);
+ TH1F* projC = (TH1F*) hC->ProjectionY(); projC->SetLineColor(kBlue);
+ TH1F* projP = (TH1F*) hP->ProjectionY(); projC->SetLineColor(kRed ); projP->Scale(scale);
+ projC->Draw(); projP->Draw("same");
+
+ //SUBTRACT PHASE SPACE
+ TCanvas* cSub = new TCanvas("cSub","cSub",500,500);
+ TH2F* hS = (TH2F*) hC->Clone(); hS->SetTitle("hS"); hS->SetName("hS");
+ hS->Add(hP,-scale);
+ hS->Draw("colz");
+
+ //RETURN
+ return hS;
+}
+
+vector<TH2F*> Parry_DrawPlots_SimPeaks(vector<int> peaks){
+ vector<TH2F*> hSims;
+
+ for(int p=0; p<peaksToFit.size(); p++){
+ string path = basePath + to_string(peaks.at(p)) + ".root";
+ string name = "peak" + to_string(peaks.at(p));
+ TCanvas* c = new TCanvas(name.c_str(),name.c_str(),500,500);
+
+ TFile *f = TFile::Open(path.c_str());
+ TTree* tPeak=(TTree*)f->Get("PhysicsTree");
+
+ tPeak->Draw(drawstring.c_str(),gateBase.c_str(),"");
+ TH2F* htemp = (TH2F*) gDirectory->Get("h");
+ htemp->SetTitle(name.c_str()); htemp->SetName(name.c_str());
+
+ hSims.push_back(htemp);
+ c->cd();
+ htemp->Draw("colz");
+ }
+
+ return hSims;
+}
+
+//void Parry_FitSingleAngle(TH1F* exp, TObjArray* sim, int angle, TFile* f){
+vector<pair<double,double>> Parry_FitSingleAngle(TH1F* exp, TObjArray* sim, int angle, TFile* f){
+
+////TEST INPUTS
+//TCanvas* test = new TCanvas("test","test",800,800);
+//exp->Draw(); exp->Print();
+//for(int i=0; i<sim->GetEntries(); i++){
+// sim->At(i)->Draw("same");
+// sim->At(i)->Print();
+//}
+
+ //DEFINITIONS
+ double frac[sim->GetEntries()]; double fracerr[sim->GetEntries()];
+ TFractionFitter* fracfit = new TFractionFitter(exp, sim);
+ fracfit->SetRangeX(exp->FindBin(7.4), exp->FindBin(10.5));
+ for(int i=0; i<sim->GetEntries(); i++){ fracfit->Constrain(i,0.0,0.4); }
+//for(int b=exp->FindBin(7.9); b<exp->FindBin(8.2); b++){ fracfit->ExcludeBin(b); }
+
+ //PERFORM FIT
+ Int_t fitstatus = fracfit->Fit();
+ if(fitstatus!=0){ cout << RED << " FAILED TO FIT!" << endl; }
+ else { cout << GREEN << " FRACTION FITTER SUCCESSFUL" << endl; }
+
+ //EXTRACT RESULTS
+ fracfit->ReleaseRangeX();
+ TH1F* result = (TH1F*) fracfit->GetPlot();
+ result->SetLineColor(kRed); result->SetLineWidth(2); result->SetLineStyle(1);
+ for(int i = 0; i<sim->GetEntries(); i++){ fracfit->GetResult(i, frac[i],fracerr[i]); }
+
+ //OUTPUT RESULTS TO SCREEN
+ for(int i = 0; i<sim->GetEntries(); i++){ cout << frac[i] << "+-" << fracerr[i] << " || "; }
+ cout << endl;
+ cout << "Chi2 / NDF = " << fracfit->GetChisquare() << " / "
+ << fracfit->GetNDF() << " = "
+ << fracfit->GetChisquare() / fracfit->GetNDF()
+ << RESET << endl;
+
+ //SEE OUTPUTS
+ string canvname = "canv" + to_string(angle);
+ TCanvas* output = new TCanvas(canvname.c_str(),canvname.c_str(),500,500);
+ exp->Draw("EL");
+ result->Draw("same hist");
+ for(int i=0; i<sim->GetEntries(); i++){
+ double scale = (double)(exp->Integral()/((TH1F*) sim->At(i))->Integral())*(double)frac[i];
+ cout << BLUE << " SCALING PEAK " << i << " BY VALUE " << scale << RESET << endl;
+ ((TH1F*) sim->At(i))->Scale(scale);
+ sim->At(i)->Draw("same hist");
+ }
+
+ f->cd();
+ output->Write();
+
+ vector<pair<double,double>> AllPeaksOneAngle;
+
+ for(int i=0; i<sim->GetEntries(); i++){
+ double tmp = (double)(((double)exp->Integral())/(((TH1F*) sim->At(i))->Integral()*(double)frac[i]));
+ pair<double,double> PointErr ( tmp, (fracerr[i]/frac[i])*tmp );
+ cout << exp->Integral() << " "
+ << ((TH1F*)sim->At(i))->Integral() << " "
+ << frac[i] << " "
+ << tmp << " | "
+ << ((TH1F*) sim->At(i))->Integral() << " "
+ << (double)frac[i] << " "
+ << exp->Integral() << " "
+ << endl;
+ AllPeaksOneAngle.push_back(PointErr);
+ }
+
+//return fracfit->GetChisquare() / fracfit->GetNDF();
+ return AllPeaksOneAngle;
+}
+
+void Parry_RunWholeRange(){
+ //Need to automate filling of vector<int> peaksToFit
+ cout << " Fitting peaks: ";
+ for(int i=0; i<peaksToFit.size(); i++){ cout << peaksToFit[i] << " ";}
+ cout << endl;
+
+ //Draw EXPERIMENTAL plots
+ TH2F* hExp2D = Parry_DrawPlots_ExpSubPSp(tree_both);
+
+ //Draw SIMULATED plots
+ vector<TH2F*> hSim2D = Parry_DrawPlots_SimPeaks(peaksToFit);
+
+ //Output file
+ TFile* f = new TFile("ParryFunctionOutputs.root","RECREATE");
+
+ //Cycle through angles to fit each single angle with peaks, and fill peak DCS plots
+ for(int a=0; a<angleBinNumb; a++){
+ string n = "hAng" + to_string(a);
+ TH1F* hExp1D = (TH1F*) hExp2D->ProjectionY(n.c_str(), a, a+1);
+ for(int bin=0; bin<hExp1D->GetNbinsX(); bin++){
+ if(hExp1D->GetBinContent(bin)<0){
+ hExp1D->SetBinContent(bin,0);
+ cout << "flatten bin " << bin << endl;
+ }
+ }
+ hExp1D->SetLineColor(kBlack); hExp1D->SetLineWidth(2); hExp1D->SetLineStyle(1);
+ hExp1D->GetYaxis()->SetRangeUser(0,500);
+ TObjArray *fitSim = new TObjArray(peaksToFit.size());
+ for(int p=0; p<peaksToFit.size(); p++){
+ string n2 = n + "_Peak" + to_string(p);
+ TH1F* hSim1D = (TH1F*) hSim2D[p]->ProjectionY(n2.c_str(), a, a+1);
+ hSim1D->SetLineColor(kBlue); hSim1D->SetLineWidth(2); hSim1D->SetLineStyle(7);
+ fitSim->Add(hSim1D);
+
+ string graphname = "grPeak" + to_string(p);
+ auto gr = new TGraphErrors();
+ gr->SetTitle(graphname.c_str()); gr->SetName(graphname.c_str());
+ parryPeakAoSA.push_back(gr);
+ }
+
+ //Fit all peaks to single angle
+ vector<pair<double,double>> AllPeaksOneAngle = Parry_FitSingleAngle(hExp1D, fitSim, a, f);
+
+ //Extract all angles for single peaks, fill area (later divide by SA)
+ for(int p=0; p<peaksToFit.size(); p++){
+ parryPeakAoSA[p]->SetPoint(a, angleMin+(a*angleBinWidth), AllPeaksOneAngle[p].first);
+ parryPeakAoSA[p]->SetPointError(a, 0., AllPeaksOneAngle[p].second);
+ }
+ }
+
+ //Divide DCS by SA[sr] to get dSdO
+ string filepathSA = "./SolidAngleHistFiles/SAHF_Sim_19Odp_";
+ for(int p=0; p<peaksToFit.size(); p++){
+ string filenameSA = filepathSA + to_string(peaksToFit[p]) + ".root";
+ TFile* fileSA = TFile::Open(filenameSA.c_str());
+ cout << " Using solid angle file " << filenameSA << endl;
+ TH1F* SolidAngle = (TH1F*) fileSA->FindObjectAny("SolidAngle_Lab_MG");
+
+ for(int a=0; a<angleBinNumb; a++){
+ int binmin = SolidAngle->FindBin(angleMin+(a*angleBinWidth)-(0.5*angleBinWidth));
+ int binmax = SolidAngle->FindBin(angleMin+(a*angleBinWidth)+(0.5*angleBinWidth));
+
+ double SAerr; //[sr]
+ double SA = SolidAngle->IntegralAndError(binmin, binmax, SAerr); //[sr]
+
+ double val = parryPeakAoSA[p]->GetPointY(a);
+ double err = parryPeakAoSA[p]->GetErrorY(a);
+ parryPeakAoSA[p]->SetPointY(a, val/SA);
+ parryPeakAoSA[p]->SetPointError(a, 0., (val/SA) * sqrt( pow( err/val, 2) + pow( SAerr/SA, 2) ));
+ }
+ }
+
+ //Plot AoSA for each peak
+ for(int p=0; p<peaksToFit.size(); p++){
+ string cAoSAn = "cAoSA" + to_string(p);
+ TCanvas* cAoSA = new TCanvas(cAoSAn.c_str(),cAoSAn.c_str(),500,500);
+ cAoSA->SetLogy();
+ parryPeakAoSA[p]->GetYaxis()->SetRangeUser(1e-0,1e7);
+ parryPeakAoSA[p]->SetTitle(to_string(peaksToFit[p]).c_str());
+ parryPeakAoSA[p]->Draw();
+ }
+
+
+ cout << "!!!!!!!!!!!! MADE IT OUT !!!!!!!!!!!!" << endl;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/Projects/e775s/analysis/Plots_Simulation.C b/Projects/e775s/analysis/Plots_Simulation.C
new file mode 100644
index 0000000000000000000000000000000000000000..50c21377fcfd6200ed8e495f65b023a2677d0b45
--- /dev/null
+++ b/Projects/e775s/analysis/Plots_Simulation.C
@@ -0,0 +1,60 @@
+#include "Plots_Functions.h"
+
+void DrawPhaseSpaceCuts_Ex(){
+ TCanvas *c = new TCanvas("c","c",1000,1000);
+
+ tree_PSp00->Draw("Ex>>hNoCut(100,7,11)","","");
+ TH2F* hNoCut = (TH2F*) gDirectory->Get("hNoCut"); hNoCut->SetLineColor(kBlack); hNoCut->SetLineWidth(2);
+
+//tree_PSp00->Draw("Ex>>h1p1(100,7,11)","ELab<(0.05458*ThetaLab)-4.576 && ELab>(-0.05458*ThetaLab)+6.776 && ELab>(0.05458*ThetaLab)-4.576","same");
+ tree_PSp00->Draw("Ex>>h1p2(100,7,11)","ELab<(0.05458*ThetaLab)-4.476 && ELab>(-0.05458*ThetaLab)+6.876 && ELab>(0.05458*ThetaLab)-4.676","same");
+//tree_PSp00->Draw("Ex>>h1p3(100,7,11)","ELab<(0.05458*ThetaLab)-4.376 && ELab>(-0.05458*ThetaLab)+6.976 && ELab>(0.05458*ThetaLab)-4.776","same");
+ tree_PSp00->Draw("Ex>>h1p4(100,7,11)","ELab<(0.05458*ThetaLab)-4.276 && ELab>(-0.05458*ThetaLab)+7.076 && ELab>(0.05458*ThetaLab)-4.876","same");
+//tree_PSp00->Draw("Ex>>h1p5(100,7,11)","ELab<(0.05458*ThetaLab)-4.176 && ELab>(-0.05458*ThetaLab)+7.176 && ELab>(0.05458*ThetaLab)-4.976","same");
+ tree_PSp00->Draw("Ex>>h1p6(100,7,11)","ELab<(0.05458*ThetaLab)-4.076 && ELab>(-0.05458*ThetaLab)+7.276 && ELab>(0.05458*ThetaLab)-5.076","same");
+//tree_PSp00->Draw("Ex>>h1p7(100,7,11)","ELab<(0.05458*ThetaLab)-3.976 && ELab>(-0.05458*ThetaLab)+7.376 && ELab>(0.05458*ThetaLab)-5.176","same");
+ tree_PSp00->Draw("Ex>>h1p8(100,7,11)","ELab<(0.05458*ThetaLab)-3.876 && ELab>(-0.05458*ThetaLab)+7.476 && ELab>(0.05458*ThetaLab)-5.276","same");
+//tree_PSp00->Draw("Ex>>h1p9(100,7,11)","ELab<(0.05458*ThetaLab)-3.776 && ELab>(-0.05458*ThetaLab)+7.576 && ELab>(0.05458*ThetaLab)-5.376","same");
+ tree_PSp00->Draw("Ex>>h2p0(100,7,11)","ELab<(0.05458*ThetaLab)-3.676 && ELab>(-0.05458*ThetaLab)+7.676 && ELab>(0.05458*ThetaLab)-5.476","same");
+//tree_PSp00->Draw("Ex>>h2p1(100,7,11)","ELab<(0.05458*ThetaLab)-3.576 && ELab>(-0.05458*ThetaLab)+7.776 && ELab>(0.05458*ThetaLab)-5.576","same");
+ tree_PSp00->Draw("Ex>>h2p2(100,7,11)","ELab<(0.05458*ThetaLab)-3.476 && ELab>(-0.05458*ThetaLab)+7.876 && ELab>(0.05458*ThetaLab)-5.676","same");
+//tree_PSp00->Draw("Ex>>h2p3(100,7,11)","ELab<(0.05458*ThetaLab)-3.376 && ELab>(-0.05458*ThetaLab)+7.976 && ELab>(0.05458*ThetaLab)-5.776","same");
+ tree_PSp00->Draw("Ex>>h2p4(100,7,11)","ELab<(0.05458*ThetaLab)-3.276 && ELab>(-0.05458*ThetaLab)+8.076 && ELab>(0.05458*ThetaLab)-5.876","same");
+//tree_PSp00->Draw("Ex>>h2p5(100,7,11)","ELab<(0.05458*ThetaLab)-3.176 && ELab>(-0.05458*ThetaLab)+8.176 && ELab>(0.05458*ThetaLab)-5.976","same");
+
+//TH2F* h1p1 = (TH2F*) gDirectory->Get("h1p1"); h1p1->SetLineColor(2); h1p1->SetLineWidth(2); h1p1->SetFillColor(2); h1p1->SetFillStyle(3004);
+ TH2F* h1p2 = (TH2F*) gDirectory->Get("h1p2"); h1p2->SetLineColor(3); h1p2->SetLineWidth(2); h1p2->SetFillColor(3); h1p2->SetFillStyle(3004);
+//TH2F* h1p3 = (TH2F*) gDirectory->Get("h1p3"); h1p3->SetLineColor(4); h1p3->SetLineWidth(2); h1p3->SetFillColor(4); h1p3->SetFillStyle(3004);
+ TH2F* h1p4 = (TH2F*) gDirectory->Get("h1p4"); h1p4->SetLineColor(6); h1p4->SetLineWidth(2); h1p4->SetFillColor(6); h1p4->SetFillStyle(3004);
+//TH2F* h1p5 = (TH2F*) gDirectory->Get("h1p5"); h1p5->SetLineColor(7); h1p5->SetLineWidth(2); h1p5->SetFillColor(7); h1p5->SetFillStyle(3004);
+ TH2F* h1p6 = (TH2F*) gDirectory->Get("h1p6"); h1p6->SetLineColor(2); h1p6->SetLineWidth(2); h1p6->SetFillColor(2); h1p6->SetFillStyle(3004);
+//TH2F* h1p7 = (TH2F*) gDirectory->Get("h1p7"); h1p7->SetLineColor(3); h1p7->SetLineWidth(2); h1p7->SetFillColor(3); h1p7->SetFillStyle(3004);
+ TH2F* h1p8 = (TH2F*) gDirectory->Get("h1p8"); h1p8->SetLineColor(4); h1p8->SetLineWidth(2); h1p8->SetFillColor(4); h1p8->SetFillStyle(3004);
+//TH2F* h1p9 = (TH2F*) gDirectory->Get("h1p9"); h1p9->SetLineColor(6); h1p9->SetLineWidth(2); h1p9->SetFillColor(6); h1p9->SetFillStyle(3004);
+ TH2F* h2p0 = (TH2F*) gDirectory->Get("h2p0"); h2p0->SetLineColor(7); h2p0->SetLineWidth(2); h2p0->SetFillColor(7); h2p0->SetFillStyle(3004);
+//TH2F* h2p1 = (TH2F*) gDirectory->Get("h2p1"); h2p1->SetLineColor(2); h2p1->SetLineWidth(2); h2p1->SetFillColor(2); h2p1->SetFillStyle(3004);
+ TH2F* h2p2 = (TH2F*) gDirectory->Get("h2p2"); h2p2->SetLineColor(3); h2p2->SetLineWidth(2); h2p2->SetFillColor(3); h2p2->SetFillStyle(3004);
+//TH2F* h2p3 = (TH2F*) gDirectory->Get("h2p3"); h2p3->SetLineColor(4); h2p3->SetLineWidth(2); h2p3->SetFillColor(4); h2p3->SetFillStyle(3004);
+ TH2F* h2p4 = (TH2F*) gDirectory->Get("h2p4"); h2p4->SetLineColor(6); h2p4->SetLineWidth(2); h2p4->SetFillColor(6); h2p4->SetFillStyle(3004);
+//TH2F* h2p5 = (TH2F*) gDirectory->Get("h2p5"); h2p5->SetLineColor(7); h2p5->SetLineWidth(2); h2p5->SetFillColor(7); h2p5->SetFillStyle(3004);
+
+}
+
+void DrawPhaseSpaceCuts_ELTL(){
+ TCanvas *c = new TCanvas("c","c",1000,1000);
+
+ c->Divide(2);
+ c->cd(1);
+ tree_PSp00->Draw("ELab:ThetaLab>>h1p2(80,100,140,80,1,3)","ELab<(0.05458*ThetaLab)-4.476 && ELab>(-0.05458*ThetaLab)+6.876 && ELab>(0.05458*ThetaLab)-4.676","colz");
+ c->cd(2);
+ tree_PSp00->Draw("ELab:ThetaLab>>h2p4(80,100,140,80,1,3)","ELab<(0.05458*ThetaLab)-3.276 && ELab>(-0.05458*ThetaLab)+8.076 && ELab>(0.05458*ThetaLab)-5.876","colz");
+
+}
+
+/////////////////////////////////////////////////////////////////////////////////////////////
+
+void Plots_Simulation(){
+
+ LoadChain_PhaseSpace();
+
+}
diff --git a/Projects/e775s/analysis/Plots_TWOFNR.h b/Projects/e775s/analysis/Plots_TWOFNR.h
new file mode 100644
index 0000000000000000000000000000000000000000..968b70d3b15166c4ad5f6185ec5fd72ce9820bf2
--- /dev/null
+++ b/Projects/e775s/analysis/Plots_TWOFNR.h
@@ -0,0 +1,186 @@
+void MoveDirectoryTo(string newDir){
+ cout << "--------------------------------------------------------" << endl;
+ cout << "Moving directory from:" << endl;
+ system("pwd");
+ cout << "to:" << endl;
+ chdir(newDir.c_str());
+ system("pwd");
+ cout << "--------------------------------------------------------" << endl;
+
+
+}
+
+TGraph* TWOFNR(double Ex, double Jf, double particleN, double particleL, double particleJ){
+
+ /* This funciton moves between directories in roder to execute TWOFNR */
+ /* and as such CANNOT be run on another machine! */
+
+
+
+ cout << "========================================================" << endl;
+ cout << " BEGINING TWOFNR AUTOMATED PROCESS " << endl;
+ cout << "========================================================" << endl;
+
+ //double Ex = 0.0;
+ double Ji = 2.5;
+ //double Jf = 0.0;
+ double QValue = 5.383 - Ex; //QValue for 19O(d,p)
+ double BeamEnergy = 8.0;
+ int beamA = 19, beamZ = 8;
+ int model_in_A = 6; //JTandy (for now...)
+ int model_in_B = 1; //Reid soft core (for now...)
+ int model_out_A = 6; //KoningDelaroche
+
+ string njjj;
+
+ string twofnrDir = "/home/paxman/Programs/twofnr";
+
+ char origDirChar[200];
+ getcwd(origDirChar,200);
+ string origDir{origDirChar};
+
+ // Moving to TWOFNR directory
+ MoveDirectoryTo(twofnrDir);
+
+ // Delete existing tran.jjj and 24.jjj files
+ remove("tran.jjj");
+ remove("24.jjj");
+
+ // Fill in.front file
+ ofstream frontinput("in.front");
+ /* name */ frontinput << "jjj" << endl;
+ /* label */ frontinput << "auto" << endl;
+ /* reaction - here (d,p) */ frontinput << 2 << endl;
+ /* calcualte entrance DW */ frontinput << 0 << endl;
+ /* calcualte exit DW */ frontinput << 0 << endl;
+ /* beam energy per nucleon */ frontinput << BeamEnergy << endl;
+ /* beam A, Z */ frontinput << beamA << " " << beamZ << endl;
+ /* default integration range */ frontinput << 1 << endl;
+ /* default number of partial waves */ frontinput << 1 << endl;
+ /* default c.o.m steps */ frontinput << "0 0 0" << endl;
+ /* L and J of transfered nucleon */ frontinput << particleL << " " << particleJ << endl;
+ /* nodes in function (0s,0p,1s...)*/ frontinput << particleN << endl;
+ /* calc using reaction Q value */ frontinput << 2 << endl;
+ /* input reaction Q value */ frontinput << QValue << endl;
+ /* no nonlocality in incident */ frontinput << 1 << endl;
+ /* spin, incident channel */ frontinput << Ji << endl;
+ /* use built-in potentials */ frontinput << 1 << endl;
+ /*** select deut potential ***/ frontinput << model_in_A << endl;
+ if(model_in_A==6){
+ /*** select deut potential ***/ frontinput << model_in_B << endl;
+ }
+ /* no nonlocality in outgoing */ frontinput << 1 << endl;
+ /* spin, outgoing channel */ frontinput << Jf << endl;
+ /* use built-in potentials */ frontinput << 1 << endl;
+ /*** select prot protential ***/ frontinput << model_out_A << endl;
+ /*** switch case is to keep second selection consistent with first ***/
+ switch(model_out_A){
+ case 1:
+ frontinput << 1 << endl;
+ break;
+ case 2:
+ frontinput << 2 << endl;
+ break;
+ case 3:
+ frontinput << 2 << endl;
+ break;
+ case 4:
+ frontinput << 2 << endl;
+ break;
+ case 5:
+ frontinput << 3 << endl;
+ break;
+ case 6:
+ frontinput << 4 << endl;
+ break;
+ }
+ /* default <p|d> vertex value */ frontinput << 1 << endl;
+ /* use that value */ frontinput << 1 << endl;
+ /* use zero range <d|p> vertex */ frontinput << 1 << endl;
+ /* default neutron binding potent. */ frontinput << "1.25 0.65" << endl;
+ /* default spin-orbit strengh */ frontinput << 6.0 << endl;
+ /* default bound st. non-locality */ frontinput << 0 << endl;
+ /* default bound st. s-o radius */ frontinput << 0 << endl;
+
+ // Close in.front file
+ frontinput.close();
+
+ // Execute front20
+ string exec_front = "(exec ";
+ exec_front.append(twofnrDir);
+ exec_front.append("/front20 < in.front > /dev/null)");
+
+ cout << "Executing: " << exec_front << endl;
+ system(exec_front.c_str());
+
+ ifstream checkfront("tran.jjj");
+ if(!checkfront){
+ cout << "!!! ERROR !!! -> front20 failed" << endl;
+ return 0;
+ } else {
+ cout << "front20 execution complete" << endl;
+ checkfront.close();
+ }
+
+
+ // Fill in.twofnr file
+ ofstream twofnrinput("in.twofnr");
+ twofnrinput << "tran.jjj" << endl;
+ twofnrinput.close();
+
+ // Execute twofnr20
+ string exec_twofnr = "(exec ";
+ exec_twofnr.append(twofnrDir);
+ exec_twofnr.append("/twofnr20 < in.twofnr > /dev/null)");
+
+ cout << "Executing: " << exec_twofnr << endl;
+ system(exec_twofnr.c_str());
+
+ ifstream checktwofnr("24.jjj");
+ if(!checktwofnr){
+ cout << "!!! ERROR !!! -> twofnr20 failed" << endl;
+ return 0;
+ } else {
+ cout << "twofnr20 execution complete" << endl;
+ checktwofnr.close();
+ }
+
+ TGraph* DCS_Lab = new TGraph("24.jjj");
+ cout << "Loaded DCS (lab), in units of mb/sr!" << endl;
+
+ // Moving back to working directory
+ MoveDirectoryTo(origDir);
+
+
+ return DCS_Lab;
+
+}
+
+void TWOFNR_CompareShapes(double Ex){
+
+ TCanvas* twofnr = new TCanvas("twofnr","twofnr",1000,1000);
+ twofnr->cd(); twofnr->SetLogy();
+
+ TGraph* s12 = TWOFNR(Ex, 2., 1., 0., 0.5);
+ TGraph* p12 = TWOFNR(Ex, 2., 0., 1., 0.5);
+ TGraph* d52 = TWOFNR(Ex, 0., 0., 2., 2.5);
+ TGraph* f72 = TWOFNR(Ex, 3., 0., 3., 3.5);
+
+ s12->SetLineColor(kRed); s12->SetLineWidth(2); s12->SetLineStyle(1);
+ p12->SetLineColor(kGreen); p12->SetLineWidth(1); p12->SetLineStyle(7);
+ d52->SetLineColor(kBlue); d52->SetLineWidth(2); d52->SetLineStyle(1);
+ f72->SetLineColor(kViolet); f72->SetLineWidth(1); f72->SetLineStyle(7);
+
+ s12->Draw();
+ p12->Draw("same");
+ d52->Draw("same");
+ f72->Draw("same");
+
+ auto legend = new TLegend(0.6,0.1,0.9,0.9);
+ legend->AddEntry(s12,"s1/2","L");
+ legend->AddEntry(p12,"p1/2","L");
+ legend->AddEntry(d52,"d5/2","L");
+ legend->AddEntry(f72,"f7/2","L");
+ legend->Draw();
+
+}
diff --git a/Projects/e775s/analysis/run_autopsp.sh b/Projects/e775s/analysis/run_autopsp.sh
new file mode 100755
index 0000000000000000000000000000000000000000..fd136ee17359fa621cdb9d6f0841411aed08b915
--- /dev/null
+++ b/Projects/e775s/analysis/run_autopsp.sh
@@ -0,0 +1,3 @@
+#!/bin/bash
+./run_psp.sh PhaseSpace_InitialTest 0000
+./run_psp.sh PhaseSpace_InitialTest 0096
diff --git a/Projects/e775s/analysis/run_autosim.sh b/Projects/e775s/analysis/run_autosim.sh
new file mode 100755
index 0000000000000000000000000000000000000000..9bc95b15a9683f54e3fe16e574a38911c8dc2436
--- /dev/null
+++ b/Projects/e775s/analysis/run_autosim.sh
@@ -0,0 +1,30 @@
+#./run_sim.sh 19Odp 7622 # 18O(t,p)
+#./run_sim.sh 19Odp 7754 # 18O(t,p)
+#./run_sim.sh 19Odp 8554 # 18O(t,p)
+#./run_sim.sh 19Odp 8962 # 18O(t,p)
+#./run_sim.sh 19Odp 9770 # 18O(t,p)
+#./run_sim.sh 19Odp 7855 # Bohlen
+ ./run_sim.sh 19Odp 7916 # Bohlen
+ ./run_sim.sh 19Odp 8313 # Bohlen
+ ./run_sim.sh 19Odp 8561 # Bohlen
+ ./run_sim.sh 19Odp 8789 # Bohlen
+ ./run_sim.sh 19Odp 8962 # Bohlen
+
+#./run_sim.sh 19Odp 9200 # spaced out 9 MeV to 10 MeV
+#./run_sim.sh 19Odp 9400 # spaced out 9 MeV to 10 MeV
+#./run_sim.sh 19Odp 9600 # spaced out 9 MeV to 10 MeV
+
+#./run_sim.sh 19Odp 8000
+#./run_sim.sh 19Odp 8200
+#./run_sim.sh 19Odp 8400
+#./run_sim.sh 19Odp 8800
+
+#./run_sim.sh 19Odp 0000
+#./run_sim.sh 19Odp 1675
+#./run_sim.sh 19Odp 3572
+#./run_sim.sh 19Odp 4071
+#./run_sim.sh 19Odp 5004
+#./run_sim.sh 19Odp 5228
+#./run_sim.sh 19Odp 5629
+#./run_sim.sh 19Odp 7610
+#./run_sim.sh 19Odp 7865
diff --git a/Projects/e775s/analysis/run_exp.sh b/Projects/e775s/analysis/run_exp.sh
new file mode 100755
index 0000000000000000000000000000000000000000..03fbddf8d7dfc0cade9080c3005384571b99a5bb
--- /dev/null
+++ b/Projects/e775s/analysis/run_exp.sh
@@ -0,0 +1,52 @@
+#!/bin/bash
+cd ~/Programs/nptoolv3/Projects/e775s/analysis;
+cmake ./;
+make -j6;
+
+#=====================
+#====================================================
+ #rfile='./InputFiles/e775s.reaction'
+ rfile='./InputFiles/e775s_2024-07-30.reaction'
+#----------------------------------------------------
+ #dfile='./InputFiles/mugast.detector'
+ dfile='./InputFiles/mugast_2024-07-30.detector'
+#----------------------------------------------------
+ cfile='./InputFiles/Calibration.txt'
+#====================================================
+#=====================
+
+#=====================
+#====================================================
+# Symbolic link for mugast map
+#----------------------------------------------------
+# Remove old links
+ rm ~/Programs/nptoolv3/NPLib/Detectors/Mugast/MugastMap.h
+ rm ~/Programs/nptoolv3/NPLib/Detectors/Mugast/MugastReverseMap.h
+#----------------------------------------------------
+# Link "DATA"
+#echo "SYMBOLIC LINK: DATA"
+#ln -s ~/Programs/nptoolv3/NPLib/Detectors/Mugast/Zanon_Data_MugastMap.h ~/Programs/nptoolv3/NPLib/Detectors/Mugast/MugastMap.h
+#ln -s ~/Programs/nptoolv3/NPLib/Detectors/Mugast/Zanon_Data_MugastReverseMap.h ~/Programs/nptoolv3/NPLib/Detectors/Mugast/MugastReverseMap.h
+#----------------------------------------------------
+# Link "ORIG"
+ echo "SYMBOLIC LINK: ORIG"
+ ln -s ~/Programs/nptoolv3/NPLib/Detectors/Mugast/Zanon_Orig_MugastMap.h ~/Programs/nptoolv3/NPLib/Detectors/Mugast/MugastMap.h
+ ln -s ~/Programs/nptoolv3/NPLib/Detectors/Mugast/Zanon_Orig_MugastReverseMap.h ~/Programs/nptoolv3/NPLib/Detectors/Mugast/MugastReverseMap.h
+#====================================================
+#=====================
+
+#=====================
+#====================================================
+#Analysis -- CD2 Target
+ inputCD='./InputFiles/CD.runtotreat'
+ outputCD=$1'_CD'
+ #npanalysis --definition Exp -R $inputCD -E $rfile -D $dfile -C $cfile -O $outputCD;
+ npanalysis --definition Exp --definition CD2 -R $inputCD -E $rfile -D $dfile -C $cfile -O $outputCD;
+
+#Analysis -- CD2+Au Target
+ inputAu='./InputFiles/Au.runtotreat'
+ outputAu=$1'_Au'
+ npanalysis --definition Exp --definition Au -R $inputAu -E $rfile -D $dfile -C $cfile -O $outputAu;
+#====================================================
+#=====================
+
diff --git a/Projects/e775s/analysis/run_psp.sh b/Projects/e775s/analysis/run_psp.sh
new file mode 100755
index 0000000000000000000000000000000000000000..764826417408ca72ea6d3c2267a24a2d1ca6ba19
--- /dev/null
+++ b/Projects/e775s/analysis/run_psp.sh
@@ -0,0 +1,91 @@
+#!/bin/bash
+
+cd ~/Programs/nptoolv3/Projects/e775s/analysis;
+rm ./InputFiles/Auto.runtotreat
+rm ./InputFiles/AutoPSpace.reaction
+
+sp=" "
+dfile="./InputFiles/mugast_2024-07-30.detector"
+rfile_exp="./InputFiles/e775s_2024-07-30.reaction"
+rfile_sim="./InputFiles/AutoPSpace.reaction"
+
+#Variables
+E=$(echo "scale=3; $2/1000" | bc -l)
+X="-3.41"
+Y="-0.20"
+
+#====================================================
+# Symbolic link for mugast map
+#----------------------------------------------------
+# Remove old links
+ rm ~/Programs/nptoolv3/NPLib/Detectors/Mugast/MugastMap.h
+ rm ~/Programs/nptoolv3/NPLib/Detectors/Mugast/MugastReverseMap.h
+#----------------------------------------------------
+# Link "SIMU"
+#echo "SYMBOLIC LINK: SIMU"
+#ln -s ~/Programs/nptoolv3/NPLib/Detectors/Mugast/Zanon_Simu_MugastMap.h ~/Programs/nptoolv3/NPLib/Detectors/Mugast/MugastMap.h
+#ln -s ~/Programs/nptoolv3/NPLib/Detectors/Mugast/Zanon_Simu_MugastReverseMap.h ~/Programs/nptoolv3/NPLib/Detectors/Mugast/MugastReverseMap.h
+#----------------------------------------------------
+# Link "ORIG"
+ echo "SYMBOLIC LINK: ORIG"
+ ln -s ~/Programs/nptoolv3/NPLib/Detectors/Mugast/Zanon_Orig_MugastMap.h ~/Programs/nptoolv3/NPLib/Detectors/Mugast/MugastMap.h
+ ln -s ~/Programs/nptoolv3/NPLib/Detectors/Mugast/Zanon_Orig_MugastReverseMap.h ~/Programs/nptoolv3/NPLib/Detectors/Mugast/MugastReverseMap.h
+#====================================================
+
+
+#====================================================
+echo -e "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" >> ./InputFiles/AutoPSpace.reaction
+echo -e "Beam" >> ./InputFiles/AutoPSpace.reaction
+echo -e "$sp""Particle=""$sp""19O" >> ./InputFiles/AutoPSpace.reaction
+echo -e "$sp""ExcitationEnergy=""$sp""0""$sp""MeV" >> ./InputFiles/AutoPSpace.reaction
+echo -e "$sp""Energy=""$sp""151.82""$sp""MeV" >> ./InputFiles/AutoPSpace.reaction
+echo -e "$sp""SigmaEnergy=""$sp""0.0""$sp"" MeV" >> ./InputFiles/AutoPSpace.reaction
+echo -e "$sp""SigmaThetaX=""$sp""0.0""$sp""deg" >> ./InputFiles/AutoPSpace.reaction
+echo -e "$sp""SigmaPhiY=""$sp""0.0""$sp""deg" >> ./InputFiles/AutoPSpace.reaction
+echo -e "$sp""SigmaX=""$sp""1.5""$sp""millimeter" >> ./InputFiles/AutoPSpace.reaction
+echo -e "$sp""SigmaY=""$sp""3.0""$sp""millimeter" >> ./InputFiles/AutoPSpace.reaction
+echo -e "$sp""MeanThetaX=""$sp""0""$sp""deg" >> ./InputFiles/AutoPSpace.reaction
+echo -e "$sp""MeanPhiY=""$sp""0""$sp""deg" >> ./InputFiles/AutoPSpace.reaction
+echo -e "$sp""MeanX=""$sp""$X""$sp""millimeter" >> ./InputFiles/AutoPSpace.reaction
+echo -e "$sp""MeanY=""$sp""$Y""$sp""millimeter" >> ./InputFiles/AutoPSpace.reaction
+echo -e "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" >> ./InputFiles/AutoPSpace.reaction
+echo -e "PhaseSpace" >> ./InputFiles/AutoPSpace.reaction
+echo -e "$sp""Beam=""$sp""19O" >> ./InputFiles/AutoPSpace.reaction
+echo -e "$sp""Target=""$sp""2H" >> ./InputFiles/AutoPSpace.reaction
+echo -e "$sp""Daughters=""$sp""1H""$sp""19O""$sp""n" >> ./InputFiles/AutoPSpace.reaction
+echo -e "$sp""ExcitationEnergies=""$sp""0.0""$sp""$E""$sp""0.0""MeV" >> ./InputFiles/AutoPSpace.reaction
+echo -e "$sp""Fermi=""$sp""1" >> ./InputFiles/AutoPSpace.reaction
+#====================================================
+
+cmake ./;
+make -j6;
+
+nameconstruct=$1"_"$2
+echo "Name: $nameconstruct"
+echo "Reaction file: $rfile_sim"
+echo "Detector file: $dfile"
+
+for x in 1 2 3 4 5 6 7 8 9 10
+do
+ outname=$nameconstruct"-"$x
+# npsimulation -D $dfile -E $rfile_sim -B ./InputFiles/runPhaseSpace.mac --random-seed $RANDOM -O $outname;
+done
+
+outfile="PSp_"$nameconstruct
+
+echo "TTreeName" >> ./InputFiles/Auto.runtotreat
+echo " SimulatedTree" >> ./InputFiles/Auto.runtotreat
+echo "RootFileName" >> ./InputFiles/Auto.runtotreat
+
+directory='/home/paxman/Programs/nptoolv3/Outputs/Simulation/'
+
+for x in 1 2 3 4 5 6 7 8 9 10
+do
+ echo "$sp""$directory""$nameconstruct""-""$x"".root" >> ./InputFiles/Auto.runtotreat
+done
+
+npanalysis --definition PSp -R ./InputFiles/Auto.runtotreat -E $rfile_exp -D $dfile -O $outfile;
+
+#mvname="./SolidAngleHistFiles/SAHF_"$outfile".root"
+#mv ./SolidAngleHistFiles/SAHF_New.root $mvname
+
diff --git a/Projects/e775s/analysis/run_sim.sh b/Projects/e775s/analysis/run_sim.sh
new file mode 100755
index 0000000000000000000000000000000000000000..619e9cc6964318ddf85f30067b2ee92b816a5f8f
--- /dev/null
+++ b/Projects/e775s/analysis/run_sim.sh
@@ -0,0 +1,99 @@
+#!/bin/bash
+
+cd ~/Programs/nptoolv3/Projects/e775s/analysis;
+rm ./InputFiles/Auto.runtotreat
+rm ./InputFiles/Auto.reaction
+
+sp=" "
+dfile="./InputFiles/mugast_2024-07-30.detector"
+rfile="./InputFiles/Auto.reaction"
+
+#Variables
+E=$(echo "scale=3; $2/1000" | bc -l)
+X="-3.41"
+Y="-0.20"
+
+#====================================================
+# Symbolic link for mugast map
+#----------------------------------------------------
+# Remove old links
+ rm ~/Programs/nptoolv3/NPLib/Detectors/Mugast/MugastMap.h
+ rm ~/Programs/nptoolv3/NPLib/Detectors/Mugast/MugastReverseMap.h
+#----------------------------------------------------
+# Link "SIMU"
+#echo "SYMBOLIC LINK: SIMU"
+#ln -s ~/Programs/nptoolv3/NPLib/Detectors/Mugast/Zanon_Simu_MugastMap.h ~/Programs/nptoolv3/NPLib/Detectors/Mugast/MugastMap.h
+#ln -s ~/Programs/nptoolv3/NPLib/Detectors/Mugast/Zanon_Simu_MugastReverseMap.h ~/Programs/nptoolv3/NPLib/Detectors/Mugast/MugastReverseMap.h
+#----------------------------------------------------
+# Link "ORIG"
+ echo "SYMBOLIC LINK: ORIG"
+ ln -s ~/Programs/nptoolv3/NPLib/Detectors/Mugast/Zanon_Orig_MugastMap.h ~/Programs/nptoolv3/NPLib/Detectors/Mugast/MugastMap.h
+ ln -s ~/Programs/nptoolv3/NPLib/Detectors/Mugast/Zanon_Orig_MugastReverseMap.h ~/Programs/nptoolv3/NPLib/Detectors/Mugast/MugastReverseMap.h
+#====================================================
+
+#====================================================
+#Boilerplate
+echo -e "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" >> ./InputFiles/Auto.reaction
+echo -e "Beam" >> ./InputFiles/Auto.reaction
+echo -e "$sp""Particle=""$sp""19O" >> ./InputFiles/Auto.reaction
+echo -e "$sp""ExcitationEnergy=""$sp""0""$sp""MeV" >> ./InputFiles/Auto.reaction
+echo -e "$sp""Energy=""$sp""151.82""$sp""MeV" >> ./InputFiles/Auto.reaction
+echo -e "$sp""SigmaEnergy=""$sp""0.0""$sp"" MeV" >> ./InputFiles/Auto.reaction
+echo -e "$sp""SigmaThetaX=""$sp""0.0""$sp""deg" >> ./InputFiles/Auto.reaction
+echo -e "$sp""SigmaPhiY=""$sp""0.0""$sp""deg" >> ./InputFiles/Auto.reaction
+echo -e "$sp""SigmaX=""$sp""1.5""$sp""millimeter" >> ./InputFiles/Auto.reaction
+echo -e "$sp""SigmaY=""$sp""3.0""$sp""millimeter" >> ./InputFiles/Auto.reaction
+echo -e "$sp""MeanThetaX=""$sp""0""$sp""deg" >> ./InputFiles/Auto.reaction
+echo -e "$sp""MeanPhiY=""$sp""0""$sp""deg" >> ./InputFiles/Auto.reaction
+echo -e "$sp""MeanX=""$sp""$X""$sp""millimeter" >> ./InputFiles/Auto.reaction
+echo -e "$sp""MeanY=""$sp""$Y""$sp""millimeter" >> ./InputFiles/Auto.reaction
+echo -e "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" >> ./InputFiles/Auto.reaction
+echo -e "TwoBodyReaction" >> ./InputFiles/Auto.reaction
+echo -e "$sp""Beam=""$sp""19O" >> ./InputFiles/Auto.reaction
+echo -e "$sp""Target=""$sp""2H" >> ./InputFiles/Auto.reaction
+echo -e "$sp""Light=""$sp""1H" >> ./InputFiles/Auto.reaction
+echo -e "$sp""Heavy=""$sp""20O" >> ./InputFiles/Auto.reaction
+echo -e "$sp""ExcitationEnergyLight=""$sp""0.0""$sp""MeV" >> ./InputFiles/Auto.reaction
+echo -e "$sp""ExcitationEnergyHeavy=""$sp""$E""$sp""MeV" >> ./InputFiles/Auto.reaction
+echo -e "$sp""LabCrossSectionPath=""$sp""flatbackwards.txt""$sp""CSR" >> ./InputFiles/Auto.reaction
+echo -e "$sp""ShootLight=""$sp""1" >> ./InputFiles/Auto.reaction
+echo -e "$sp""ShootHeavy=""$sp""0" >> ./InputFiles/Auto.reaction
+
+#====================================================
+
+cmake ./;
+make -j6;
+
+nameconstruct=$1"_"$2
+echo "Name: $nameconstruct"
+echo "Reaction file: $rfile"
+echo "Detector file: $dfile"
+
+for x in 1 2 3 4 5 6 7 8 9 10
+#for x in 11 12 13 14 15
+do
+ outname=$nameconstruct"-"$x
+ npsimulation -D $dfile -E $rfile -B ./InputFiles/runSimulation.mac --random-seed $RANDOM -O $outname;
+done
+
+outfile="Sim_"$nameconstruct
+
+echo "TTreeName" >> ./InputFiles/Auto.runtotreat
+echo " SimulatedTree" >> ./InputFiles/Auto.runtotreat
+echo "RootFileName" >> ./InputFiles/Auto.runtotreat
+
+directory='/home/paxman/Programs/nptoolv3/Outputs/Simulation/'
+defineState='S'$2
+
+echo "Defining state = "$defineState
+
+for x in 1 2 3 4 5 6 7 8 9 10
+do
+ echo "$sp""$directory""$nameconstruct""-""$x"".root" >> ./InputFiles/Auto.runtotreat
+done
+
+npanalysis --definition Sim --definition $defineState -R ./InputFiles/Auto.runtotreat -E $rfile -D $dfile -O $outfile;
+
+mvname="./SolidAngleHistFiles/SAHF_"$outfile".root"
+mv ./SolidAngleHistFiles/SAHF_New.root $mvname
+