diff --git a/Projects/e775s/analysis/Analysis.cxx b/Projects/e775s/analysis/Analysis.cxx
index 120b33601acb9415c149b28ce3eb305d7a0583d3..8cde7c8079da5076efc5678a1eefad19a92bf03d 100644
--- a/Projects/e775s/analysis/Analysis.cxx
+++ b/Projects/e775s/analysis/Analysis.cxx
@@ -34,6 +34,9 @@ using namespace std;
#include "NPOptionManager.h"
#include "DefineColours.h"
+// Ex Observed:Expected correction factors
+double ExCorrM, ExCorrC;
+
////////////////////////////////////////////////////////////////////////////////
Analysis::Analysis(){
}
@@ -69,14 +72,20 @@ cout << "start_Init" << endl;
if(NPOptionManager::getInstance()->HasDefinition("Au")){
cout << " -- -- -- -- Au Target -- -- -- --" << endl;
runningAuTarget=true;
+ ExCorrM = 0.991; //1.0; //0.990;
+ ExCorrC = 0.074; //0.0; //0.079;
} else if(NPOptionManager::getInstance()->HasDefinition("CD2")){
cout << " -- -- -- -- CD2 Target -- -- -- --" << endl;
runningAuTarget=false;
+ ExCorrM = 0.988; //1.0; //0.977;
+ ExCorrC = 0.124; //0.0; //0.142;
} 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;
}
+ cout << " Ex corr.: y = " << ExCorrM << " x + " << ExCorrC << endl;
+ cout << " -- -- -- -- -- -- -- -- -- -- --" << endl;
// double simCentroid=-10.;
// double simSigma = 0.104;
@@ -136,24 +145,16 @@ cout << "here_InIsSimLoop" << endl;
// 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 );
}
}
@@ -161,10 +162,8 @@ cout << "here_InIsSimLoop" << endl;
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());
@@ -184,8 +183,12 @@ cout << "here_InIsSimLoop" << endl;
TargetThickness = m_DetectorManager->GetTargetThickness();
string TargetMaterial = m_DetectorManager->GetTargetMaterial();
- AuThickness = m_DetectorManager->GetBackThickness();
- string BackMaterial = m_DetectorManager->GetBackMaterial();
+ // FOLLOWING LINES DO NOT WORK!
+ // THE ARE FOR CRYOGENIC TARGETS! NOT SOLID TARGETS!
+ //AuThickness = m_DetectorManager->GetBackThickness();
+ //string BackMaterial = m_DetectorManager->GetBackMaterial();
+ // HENCE, HAVE TO HARD-CODE Au THICKNESS HERE! NOT IDEAL!
+ AuThickness = 0.0126; //12.6 um = 0.0126 mm
cout << " ---------------- Target Values ---------------- " << endl;
cout << " CD2: \tZ = " << m_DetectorManager->GetTargetZ() << " mm\tT = " << TargetThickness << " mm" << endl;
@@ -209,10 +212,13 @@ cout << "here_InIsSimLoop" << endl;
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);
+ HeavyTarget = NPL::EnergyLoss(heavy+"_"+TargetMaterial+".G4table","G4Table",100);
+ HeavyAu = NPL::EnergyLoss(heavy+"_Au.G4table" ,"G4Table",100);
+
+ /////////////////////////////////
+ //Slowing of 19O in first half of CD2
FinalBeamEnergy = BeamTarget.Slow(OriginalBeamEnergy, 0.5*TargetThickness, 0);
reaction.SetBeamEnergy(FinalBeamEnergy);
@@ -224,6 +230,30 @@ cout << "here_InIsSimLoop" << endl;
cout << "\033[36m Beam energy at middle of CD2: " << FinalBeamEnergy << "\033[37m"<< endl;
+
+ //Slowing of 20O in second half of CD2
+ double afterCDBeamEnergy = HeavyTarget.Slow(FinalBeamEnergy, 0.5*TargetThickness, 0);
+ cout << "\033[36m Beam energy after CD2: " << afterCDBeamEnergy << "\033[37m"<< endl;
+ //Slowing of 20O in full Au thickness
+ double afterAuBeamEnergy = HeavyAu.Slow(afterCDBeamEnergy, AuThickness, 0);
+ cout << "\033[36m Beam energy after Au: " << afterAuBeamEnergy << "\033[37m"<< endl;
+
+ /////////////////////////////////
+
+ cout << "\033[36m Beta at middle of CD2: "
+ << reaction.GetEnergyImpulsionLab_4().Beta() << "\033[37m"<< endl;
+
+ NPL::Reaction reactionTest = reaction;
+
+ reactionTest.SetBeamEnergy(afterCDBeamEnergy);
+ cout << "\033[36m Beta after CD2: "
+ << reactionTest.GetEnergyImpulsionLab_4().Beta() << "\033[37m"<< endl;
+ reactionTest.SetBeamEnergy(afterAuBeamEnergy);
+ cout << "\033[36m Beta after Au: "
+ << reactionTest.GetEnergyImpulsionLab_4().Beta() << "\033[37m"<< endl;
+
+ /////////////////////////////////
+
// initialize various parameters
Rand = TRandom3();
ParticleMult = 0;
@@ -432,7 +462,7 @@ void Analysis::TreatEvent(){
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_tmp = (reaction.ReconstructRelativistic(elab_tmp, thetalab_tmp, philab_tmp) / ExCorrM) - ExCorrC;
}
Ex.push_back(ex_tmp);
Ecm.push_back(elab_tmp*(AHeavy+ALight)/(4*AHeavy*cos(thetalab_tmp)*cos(thetalab_tmp)));
@@ -466,18 +496,6 @@ void Analysis::TreatEvent(){
//}
}
-// 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);
@@ -487,20 +505,6 @@ void Analysis::TreatEvent(){
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;
@@ -749,9 +753,12 @@ void Analysis::TreatEvent(){
// 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);
+ //reaction.GetKinematicLine4();
+ TVector3 beta(0,0,-reaction.GetNucleus4()->GetBeta());
+ // THIS IS THE ONE USED BY IZANON? -> // TVector3 beta(0,0,-mBeta);
+
+ RawGammaEnergy.push_back(Egamma);
+ BetaUsedByAddBackEDC.push_back(beta.Z());
double ThetaGamma = GammaDirection.Angle(BeamDirection)/deg;
@@ -767,12 +774,6 @@ void Analysis::TreatEvent(){
}
-
-
-
-
-
-
// // Agata by track
// /*
// for(int j=0; j<nbTrack; j++){ // all multiplicity
@@ -937,118 +938,11 @@ void Analysis::End(){
}
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");
@@ -1066,62 +960,12 @@ void Analysis::End(){
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);
@@ -1129,20 +973,6 @@ void Analysis::End(){
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];
@@ -1154,8 +984,6 @@ void Analysis::End(){
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];
@@ -1167,8 +995,6 @@ void Analysis::End(){
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]);}
@@ -1210,6 +1036,8 @@ cout << "start_InitOutput" << endl;
RootOutput::getInstance()->GetTree()->Branch("beta_MUGAST",&beta_MUGAST);
RootOutput::getInstance()->GetTree()->Branch("ToF",&ToF);
+ RootOutput::getInstance()->GetTree()->Branch("RawGammaEnergy",&RawGammaEnergy);
+ RootOutput::getInstance()->GetTree()->Branch("BetaUsedByAddBackEDC",&BetaUsedByAddBackEDC);
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);
@@ -1588,6 +1416,8 @@ void Analysis::ReInitValue(){
Ex.clear();
ExCalib.clear();
Ecm.clear();
+ RawGammaEnergy.clear();
+ BetaUsedByAddBackEDC.clear();
AddBack_EDC.clear();
AddBack_EDC_Event1.clear();
AddBack_EDC_Event2.clear();
diff --git a/Projects/e775s/analysis/Analysis.h b/Projects/e775s/analysis/Analysis.h
index 07ec1224f724d0a7120d6bc473f273a0e24d6535..5572335293c40c99c362fdf1d08973789b7427a6 100644
--- a/Projects/e775s/analysis/Analysis.h
+++ b/Projects/e775s/analysis/Analysis.h
@@ -47,15 +47,9 @@ int NumThetaAngleBins = 1800;// 180 = 1 deg, 360 = 0.5 deg, 900 = 0.2 deg, 1800
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;
@@ -64,12 +58,6 @@ 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();
@@ -95,6 +83,8 @@ class Analysis: public NPL::VAnalysis{
unsigned int GammaMult;
//double EDC;
std::vector<double> EDC;
+ vector<double> RawGammaEnergy;
+ vector<double> BetaUsedByAddBackEDC;
vector<double> AddBack_EDC;
vector<double> AddBack_EDC_Event1;
vector<double> AddBack_EDC_Event2;
@@ -137,6 +127,7 @@ class Analysis: public NPL::VAnalysis{
double BeamEnergy;
NPL::Reaction reaction;
+ NPL::Reaction reactionTest;
//NPL::Reaction *reaction;
// Energy loss table: the G4Table are generated by the simulation
NPL::EnergyLoss LightTarget;
diff --git a/Projects/e775s/analysis/InputFiles/runPhaseSpace.mac b/Projects/e775s/analysis/InputFiles/runPhaseSpace.mac
index d9502085bfe4544342ddc937bcdbca74211c96f2..506556ed2882f22f4f6013ba592297aba63dec25 100644
--- a/Projects/e775s/analysis/InputFiles/runPhaseSpace.mac
+++ b/Projects/e775s/analysis/InputFiles/runPhaseSpace.mac
@@ -1 +1 @@
-/run/beamOn 10000000
+/run/beamOn 100000
diff --git a/Projects/e775s/analysis/Plots_C2SFunctions.h b/Projects/e775s/analysis/Plots_C2SFunctions.h
index 16e4efa1e6d06e80069e01fe5549fb5fcc8141f4..51164b1d67d10a3bc83ab661ead6b1562d40c3bd 100644
--- a/Projects/e775s/analysis/Plots_C2SFunctions.h
+++ b/Projects/e775s/analysis/Plots_C2SFunctions.h
@@ -2,9 +2,9 @@
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);
+TH1F* PullThetaLabHistSubPSp(int i, double minTheta, double gatesize, int angBinSize);
+//TH1F* PullThetaLabHist(int i, double minTheta, double gatesize, int angBinSize);
+//TH1F* PullPhaseSpaceHist(int i, double minTheta, double gatesize, int angBinSize);
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);
@@ -29,6 +29,7 @@ bool GammaGate = false;
double globGammaGate;
bool doMix = false;
bool doDoublet = false;
+bool comparingTWOFNRs = false;
double globalBinning;
double SnTWOFNR = 7.607;
@@ -346,7 +347,8 @@ TGraph* AddTGraphs(TGraph* a, TGraph* b){
}
////////////////////////////////////////////////////////////////////////////////
-void C2S(double Energy, double Spin, double spdf, double angmom, double binning, string option){
+// Allow detailed choices (i.e. initial nucleus spin J0, and number/width/intial angle bins)
+void C2S(double Energy, double Spin, double spdf, double angmom, double binning, string option,double J0, double nAngB, double wAngB, double fAngB ){
/* Clean global variables, in case of multiple runs */
indexE = 100;
anglecentres.clear();
@@ -380,9 +382,9 @@ void C2S(double Energy, double Spin, double spdf, double angmom, double binning,
doDoublet = true;
}
- if(binning < 0.01){
+ if(binning < 0.005){
cout << BOLDRED
- << " BINNING IS BELOW 0.01 MeV! Exiting..."
+ << " BINNING IS BELOW 0.005 MeV! Exiting..."
<< RESET
<< endl;
return;
@@ -391,31 +393,29 @@ void C2S(double Energy, double Spin, double spdf, double angmom, double binning,
/* 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 J0 = 2.5; // NOW INPUT VARIABLE
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;
+ numAngleBins=nAngB; widthAngleBins=wAngB; firstAngle=fAngB; // NOW INPUT VARIABLE
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;
- }
+ CSangleH = 160.;
+ CSvertL = 0.04; //0.01;
+ CSvertH = 100.;
+
+ //// 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();
@@ -459,6 +459,7 @@ void C2S(double Energy, double Spin, double spdf, double angmom, double binning,
/* numpeaks and Energy[] defined globally in KnownPeakFitter.h */
bool found = 0;
for(int i=0;i<numPeaks;i++){
+ cout << i << " " << Energy << " " << means.at(i) << endl;
if(abs(Energy-means.at(i))<0.01){
cout << "========================================================" << endl;
cout << "Identified as state #" << i << ", E = " << means.at(i) << endl;
@@ -469,6 +470,17 @@ void C2S(double Energy, double Spin, double spdf, double angmom, double binning,
statename = ss.str();
}
}
+// for(int i=0;i<numVoigt;i++){
+// if(abs(Energy-meansVoigt.at(i))<0.01){
+// cout << "========================================================" << endl;
+// cout << "Identified as state #" << i << ", E = " << meansVoigt.at(i) << endl;
+// indexE = numPeaks - 1 + i;
+// found = 1;
+// stringstream ss;
+// ss << setfill('0') << setw(4) << (int)((meansVoigt.at(i)+0.0001)*1000.);
+// statename = ss.str();
+// }
+// }
if(!found){
cout << "========================================================" << endl;
cout << "NO STATE AT THAT ENERGY INDENTIFIED!! CHECK KNOWN PEAKS!!" << endl;
@@ -496,13 +508,13 @@ void C2S(double Energy, double Spin, double spdf, double angmom, double binning,
string objName = "SolidAngle_Lab_MG";
TH1F* SolidAngle = (TH1F*) file->FindObjectAny(objName.c_str());
- TCanvas* c_SolidAngle = new TCanvas("c_SolidAngle","c_SolidAngle",1000,1000);
+ TCanvas* c_SolidAngle = new TCanvas("c_SolidAngle","c_SolidAngle",500,500);
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);
+ TCanvas* c_PeakArea = new TCanvas("c_PeakArea","c_PeakArea",500,500);
vector<vector<double>> areaArray = GetExpDiffCross(means.at(indexE));
// delete c_PeakArea;
@@ -580,7 +592,7 @@ void C2S(double Energy, double Spin, double spdf, double angmom, double binning,
gROOT->SetBatch(0);
/* Graph of Area/Solid Angle*/
- TCanvas* c_AoSA = new TCanvas("c_AoSA","c_AoSA",1000,1000);
+ TCanvas* c_AoSA = new TCanvas("c_AoSA","c_AoSA",500,500);
c_AoSA->SetLogy();
TGraphErrors* gAoSA = new TGraphErrors(
anglecentres.size(), //n
@@ -592,7 +604,7 @@ void C2S(double Energy, double Spin, double spdf, double angmom, double binning,
gAoSA->Draw();
/* Graph of experimental diff. cross sect (dSigma/dOmega) */
- TCanvas* c_dSdO = new TCanvas("c_dSdO","c_dSdO",1000,1000);
+ TCanvas* c_dSdO = new TCanvas("c_dSdO","c_dSdO",500,500);
c_dSdO->SetLogy();
TGraphErrors* gdSdO = new TGraphErrors(
anglecentres.size(),
@@ -611,35 +623,46 @@ void C2S(double Energy, double Spin, double spdf, double angmom, double binning,
/* 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);
+ TCanvas* c_TWOFNR = new TCanvas("c_TWOFNR","c_TWOFNR",500,500);
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");
- }
+
+ TGraph* TheoryDiffCross; TGraph* TheoryDiffCross2;
+ TGraph* TDCS_s; TGraph* TDCS_p; TGraph* TDCS_d; TGraph* TDCS_f; TGraph* TDCS_g; TGraph* TDCS_h;
+ if(comparingTWOFNRs==0){ //normal state
+ TheoryDiffCross = TWOFNR(meanTWOFNR, Spin, nodes, spdf, angmom, J0);
+ TheoryDiffCross->GetYaxis()->SetTitle("d#sigma/d#Omega [mb/msr]"); //msr set in func above
+ TheoryDiffCross->GetXaxis()->SetTitle("ThetaLab [deg]");
+ TheoryDiffCross->Draw();
+
+ if(doMix){
+ TheoryDiffCross2 = TWOFNR(meanTWOFNR, Spin, nodes2, spdf2, angmom2, J0);
+ 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;
+ /** 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;
}
- cout << setprecision(3);
- cout << "......................" << endl;
+ } else { // running for ALL ell-transfers, to compare shape
+ cout << " TWOFNR s-wave..." << endl; TDCS_s = TWOFNR(meanTWOFNR, 2, 1, 0, 0.5, 2.5); TDCS_s->SetTitle("TDCS_s");
+ cout << " TWOFNR p-wave..." << endl; TDCS_p = TWOFNR(meanTWOFNR, 0, 0, 1, 0.5, 0.5); TDCS_p->SetTitle("TDCS_p");
+ cout << " TWOFNR d-wave..." << endl; TDCS_d = TWOFNR(meanTWOFNR, 0, 0, 2, 2.5, 2.5); TDCS_d->SetTitle("TDCS_d");
+ cout << " TWOFNR f-wave..." << endl; TDCS_f = TWOFNR(meanTWOFNR, 1, 0, 3, 3.5, 2.5); TDCS_f->SetTitle("TDCS_f");
+ cout << " TWOFNR g-wave..." << endl; TDCS_g = TWOFNR(meanTWOFNR, 2, 0, 4, 4.5, 2.5); TDCS_g->SetTitle("TDCS_g");
+ cout << " TWOFNR h-wave..." << endl; TDCS_h = TWOFNR(meanTWOFNR, 3, 0, 5, 5.5, 2.5); TDCS_h->SetTitle("TDCS_h");
}
/** **** **/
/* Using Chi2 minimizaiton */
if(loud){cout << "USING CHI2 MINIMIZAITON..." << endl;}
- TCanvas* c_Chi2Min = new TCanvas("c_Chi2Min","c_Chi2Min",1000,1000);
+ TCanvas* c_Chi2Min = new TCanvas("c_Chi2Min","c_Chi2Min",500,500);
gStyle->SetPadLeftMargin(0.12);
gStyle->SetPadRightMargin(0.03);
@@ -660,25 +683,44 @@ void C2S(double Energy, double Spin, double spdf, double angmom, double binning,
pad2->Draw();
pad1->cd();
- TGraphErrors* gdSdO2 = new TGraphErrors(*gdSdO);
- gdSdO2->SetName("gdSdO2");
- gdSdO2->SetTitle("gdSdO2");
-
TGraph* Final;
+ TGraph* Fs; TGraph* Fp; TGraph* Fd; TGraph* Ff; TGraph* Fg; TGraph* Fh;
+ if(comparingTWOFNRs==0){ //normal state
+ TGraphErrors* gdSdO2 = new TGraphErrors(*gdSdO);
+ gdSdO2->SetName("gdSdO2");
+ gdSdO2->SetTitle("gdSdO2");
- 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);
+ 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);
+ } else { // running for ALL ell-transfers, to compare shape
+ Fs = FindNormalisation(TDCS_s,gdSdO);
+ Fp = FindNormalisation(TDCS_p,gdSdO);
+ Fd = FindNormalisation(TDCS_d,gdSdO);
+ Ff = FindNormalisation(TDCS_f,gdSdO);
+ Fg = FindNormalisation(TDCS_g,gdSdO);
+ Fh = FindNormalisation(TDCS_h,gdSdO);
+ Fs->SetLineColor(kViolet);
+ Fp->SetLineColor(kRed);
+ Fd->SetLineColor(kGreen);
+ Ff->SetLineColor(kBlue); Ff->SetLineStyle(7);
+ Fg->SetLineColor(kCyan); Fg->SetLineStyle(7);
+ Fh->SetLineColor(kOrange); Fh->SetLineStyle(7);
+
+ gdSdO->SetLineColor(kBlack);
+ gdSdO->SetMarkerColor(kBlack);
+ gdSdO->SetMarkerStyle(21);
+ }
/* Construct file name string */
/**/ ostringstream tempstream;
@@ -706,16 +748,26 @@ void C2S(double Energy, double Spin, double spdf, double angmom, double binning,
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->GetYaxis()->SetRangeUser(CSvertL,CSvertH);
gdSdO->GetXaxis()->SetNdivisions(512,kTRUE);
gdSdO->Draw("AP");
- Final->Draw("SAME");
+ if(comparingTWOFNRs==0){ //normal state
+ Final->Draw("SAME");
+ } else {
+ Fs->Draw("SAME");
+ Fp->Draw("SAME");
+ Fd->Draw("SAME");
+ Ff->Draw("SAME");
+ Fg->Draw("SAME");
+ Fh->Draw("SAME");
+ }
pad2->cd();
- GenerateResidual(gdSdO,Final);
+ if(comparingTWOFNRs==0){ //normal state
+ GenerateResidual(gdSdO,Final);
+ }
string savestring1 = "./C2S_Outputs/"+tempstr+".root";
string savestring2 = "./C2S_Outputs/"+tempstr+".pdf";
@@ -733,6 +785,21 @@ void C2S(double Energy, double Spin, double spdf, double angmom, double binning,
//delete c_dSdO;
}
+////////////////////////////////////////////////////////////////////////////////
+// For simple states, don't require detailed choices
+void C2S(double Energy, double Spin, double spdf, double angmom, double binning, string option){
+ C2S(Energy, Spin, spdf, angmom, binning, option, 2.5, 19, 2.0, 106.0);
+}
+
+
+////////////////////////////////////////////////////////////////////////////////
+// Comparing one state to many TWOFNRs
+void C2S_CompareTWOFNRs(double Energy, double Spin, double spdf, double angmom, double binning, string option, double J0, double nAngB, double wAngB, double fAngB ){
+ comparingTWOFNRs = true;
+ C2S(Energy, Spin, spdf, angmom, binning, option, J0, nAngB, wAngB, fAngB);
+ comparingTWOFNRs = false;
+}
+
////////////////////////////////////////////////////////////////////////////////
vector<vector<double>> GetExpDiffCross(double Energy){
cout << "========================================================" << endl;
@@ -747,7 +814,7 @@ vector<vector<double>> GetExpDiffCross(double Energy){
double trackScale = 0.0;
/* Determine scaling factor for PhaseSpace */
- TCanvas* c_ExSubPSpace = new TCanvas("c_ExSubPSpace","c_ExSubPSpace",1000,1000);
+ TCanvas* c_ExSubPSpace = new TCanvas("c_ExSubPSpace","c_ExSubPSpace",500,500);
//if(scaleTogether){
// cout << BLUE << endl;
// cout << " SUBTRACTING PHASE SPACE TOGETHER! " << endl;
@@ -822,24 +889,29 @@ vector<vector<double>> GetExpDiffCross(double Energy){
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);
+ TCanvas* c_peakFits = new TCanvas("c_peakFits",tmp2.c_str(),500,500);
/* 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);
- }
- }
-
+ TH1F *gate; // HERE< TAKE VERSION WIH PSP ALREADY REMOVED!
+ gate = PullThetaLabHistSubPSp(i,firstAngle,widthAngleBins,(int) widthAngleBins);
+
+ ///* Retrieve theta-gated Ex TH1F from file GateThetaLabHistograms.root */
+ ///* To change angle gates, run GateThetaLab_MultiWrite() */
+ //TH1F *gate, *pspace;
+ //gate = PullThetaLabHist(i,firstAngle,widthAngleBins,(int) widthAngleBins);
+ //pspace = PullPhaseSpaceHist(i,firstAngle,widthAngleBins,(int) 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;
@@ -890,7 +962,7 @@ vector<vector<double>> GetExpDiffCross(double Energy){
/* Write PS-subtracted spectrum to list */
//list->Add(gate);
//list->Add(c_peakFits);
- gate->GetXaxis()->SetRangeUser(-1.0,+8.0);
+ gate->GetXaxis()->SetRangeUser(-1.0,+10.0);
string filename = "./C2S_Outputs/";
filename += tmp2 + ".root";
c_peakFits->SaveAs(filename.c_str());
@@ -906,26 +978,30 @@ vector<vector<double>> GetExpDiffCross(double Energy){
<< " +- " << AllPeaks_OneGate[indexE][2]
<< endl;
}
-
+cout << "here.... indexE = " << indexE << endl;
/* Add min and max angle to end of relevant OneGate vector */
AllPeaks_OneGate[indexE].push_back(min);
AllPeaks_OneGate[indexE].push_back(max);
+cout << "here...." << endl;
/* Push relevant OneGate vector to end of AllGates vector */
OnePeak_AllGates.push_back(AllPeaks_OneGate[indexE]);
//delete c_peakFits;
}
+cout << "here...." << endl;
/* Write PS-subtracted spectrum to file */
//TFile* outfile = new TFile("GateThetaLab_ExMinusGatePhaseSpace.root","RECREATE");
//list->Write("GateThetaLab_ExMinusPhaseSpace",TObject::kSingleKey);
+cout << "here...." << endl;
return OnePeak_AllGates;
}
////////////////////////////////////////////////////////////////////////////////
-TH1F* PullThetaLabHist(int i, double minTheta, double gatesize){
- string name = "./C2S_Outputs/GateThetaLabHist_Experiment.root";
+TH1F* PullThetaLabHistSubPSp(int i, double minTheta, double gatesize, int angBinSize){
+ string name = "./C2S_Outputs/GateThetaLabHist_Experiment_" + to_string(angBinSize) + "deg_SubPSp.root";
+//string name = "./C2S_Outputs/GateThetaLabHist_Parallel_Experiment_SubPSp_" + to_string(angBinSize) + "deg.root";
TFile* file = new TFile(name.c_str(),"READ");
string histname = "TLab"
+ to_string((int)(minTheta+(i*gatesize))) + "to"
@@ -934,6 +1010,7 @@ TH1F* PullThetaLabHist(int i, double minTheta, double gatesize){
TList *list = (TList*)file->Get("GateThetaLabHistograms");
TH1F* hist = (TH1F*)list->FindObject(histname.c_str());
+ hist->GetXaxis()->SetRangeUser(-2.,10.);
double pulledBin = hist->GetBinWidth(10);
int ratio = (int)(globalBinning / pulledBin);
@@ -943,54 +1020,47 @@ TH1F* PullThetaLabHist(int i, double minTheta, double gatesize){
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;
-}
+//////////////////////////////////////////////////////////////////////////////////
+//TH1F* PullThetaLabHist(int i, double minTheta, double gatesize, int angBinSize){
+// string name = "./C2S_Outputs/GateThetaLabHist_Experiment_" + to_string(angBinSize) + "deg.root";
+// //string name = "./C2S_Outputs/GateThetaLabHist_Experiment_" + to_string(angBinSize) + "deg_SubPSp.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());
+// hist->GetXaxis()->SetRangeUser(-2.,10.);
+//
+// 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, int angBinSize){
+// string name = "./C2S_Outputs/GateThetaLabHist_PhaseSpace_" + to_string(angBinSize) + "deg.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){
diff --git a/Projects/e775s/analysis/Plots_Experiment.C b/Projects/e775s/analysis/Plots_Experiment.C
index 1ed0a488e2fd06d1760d4c920a2cceb245c186b0..b51ab6063cd6c25ac8afa7f925b94d409e70f10c 100644
--- a/Projects/e775s/analysis/Plots_Experiment.C
+++ b/Projects/e775s/analysis/Plots_Experiment.C
@@ -1,21 +1,28 @@
#include "Plots_Functions.h"
#include "Plots_C2SFunctions.h"
#include "Plots_ParryFunctions.h"
-
+#include "Plots_ConvoluteBreitWignerGaus.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 << " C2S( Ex, Spin, L, J, ExBin, \"option\", SpinIn, nAngB, wAngBin, fAngBin) " << endl;
+ cout << " standard values: 2.5, 25, 2, 106 " << endl;
+ cout << "======================================" << endl;
+ cout << " -- 0.000 MeV, 0+, d5/2: C2S(0.000,0,2,2.5,0.02,\" \",2.5,17,2,116) low&high angs bad" << endl;
+ cout << " -- 1.675 MeV, 2+, d5/2: C2S(1.675,2,2,2.5,0.02,\" \",2.5,21,2,106) high angs bad" << endl;
+ cout << " -- 3.572 MeV, 4+, d5/2: C2S(3.572,4,2,2.5,0.02,\" \",2.5,24,2,106) " << endl;
+ cout << " -- 4.071 MeV, 2+, sdmx: C2S(4.071,2,0,0.5,0.02,\"mixed\",2.5,24,2,106)" << endl;
+ cout << " -- 4.456 MeV, 0+, xx/2: C2S(4.456,0,2,2.5,0.02,\" \",2.5, 6,8,106) 6p-2h state" << endl;
+ cout << " -- 5.004 MeV, ?+, xx/2: C2S(5.004,2,2,2.5,0.02,\" \",2.5, 6,8,106) " << endl;
+ cout << " -- 5.228 MeV, ?+, s1/2: C2S(5.228,2,0,0.5,0.02,\"mixed\",2.5,24,2,106)" << endl;
+ cout << " -- 5.629 MeV, ?+, xx/2: C2S(5.629,2,2,2.5,0.02,\"mixed\",2.5,24,2,106)" << endl;
+ cout << " -- 7.252 MeV, 5-? xx/2: C2S(7.252,5,2,1.5,0.02,\" \",2.5,11,2,106) Tilley: h-wave, but we pop via d-wave! 3/2 or 5/2?" << endl;
+ cout << " -- 7.622 MeV, 4+? xx/2: C2S(7.622,4,2,1.5,0.02,\" \",2.5,11,2,106)" << endl;
+ cout << " -- 7.855 MeV, 5-? xx/2: C2S(7.855,2,0,0.5,0.02,\" \",2.5,11,2,106) looks s-wave?? how?? 5- accoring to LaFrance" << endl;
+ cout << " -- 8.313 MeV, ??? xx/2: C2S(8.313,4,2,1.5,0.02,\" \",2.5,11,2,106)" << endl;
cout << "======================================" << endl;
}
-
void Plots_Experiment(){
@@ -31,12 +38,26 @@ void Plots_Experiment(){
// TCutG* cutGood = LoadCut_Good();
// TCutG* cutBad = LoadCut_Bad();
+
+
+ //== Assignments for Plots_KnownPeakFitter() ==
+ means.reserve(meansGauss.size() + meansVoigt.size());
+ means.insert(means.end(), meansGauss.begin(),meansGauss.end());
+ means.insert(means.end(), meansVoigt.begin(),meansVoigt.end());
+
+ numGauss = meansGauss.size();
+ numVoigt = meansVoigt.size();
+ numPeaks = numGauss + numVoigt;
+ //=============================================
+
+
+
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 << "Time gate GOOD = abs(T_FPMW_HF-13300)<1500" << endl;
+ cout << "Time gate BAD = abs(T_FPMW_HF-10500)<1500" << 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();
+ C2S();
}
diff --git a/Projects/e775s/analysis/Plots_Functions.h b/Projects/e775s/analysis/Plots_Functions.h
index f6b86a5382d8472ee49c5914c5e22bfbe3fce0ce..2c8fce13166d7b77f621361542f01903d1f444d0 100644
--- a/Projects/e775s/analysis/Plots_Functions.h
+++ b/Projects/e775s/analysis/Plots_Functions.h
@@ -14,7 +14,7 @@ TChain* tree_cd=NULL;
TChain* tree_au=NULL;
TChain* tree_both=NULL;
TChain* tree_raw=NULL;
-TChain* tree=NULL;
+//TChain* tree=NULL;
TChain* tree_PSp00=NULL;
TChain* tree_PSp96=NULL;
@@ -25,9 +25,7 @@ auto colourB = new TColor((Int_t)TColor::GetFreeColorIndex(), 25./256., 101./25
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)
@@ -56,11 +54,28 @@ TChain* Chain(string TreeName, vector<string>& file, bool EventList)
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");
+//files.push_back("~/Programs/nptoolv3/Outputs/Analysis/19Odp_XYZT_0003_CD.root"); files.push_back("~/Programs/nptoolv3/Outputs/Analysis/19Odp_XYZT_0003_Au.root");
+//files.push_back("~/Programs/nptoolv3/Outputs/Analysis/19Odp_2024-08-28_CD.root"); files.push_back("~/Programs/nptoolv3/Outputs/Analysis/19Odp_2024-08-28_Au.root");
+//files.push_back("~/Programs/nptoolv3/Outputs/Analysis/19Odp_2024-09-05_CD.root"); files.push_back("~/Programs/nptoolv3/Outputs/Analysis/19Odp_2024-09-05_Au.root");
+//files.push_back("~/Programs/nptoolv3/Outputs/Analysis/19Odp_2024-09-06_Au.root"); files.push_back("~/Programs/nptoolv3/Outputs/Analysis/19Odp_2024-09-06_CD.root");
+ files.push_back("~/Programs/nptoolv3/Outputs/Analysis/19Odp_2024-09-10_Au.root"); files.push_back("~/Programs/nptoolv3/Outputs/Analysis/19Odp_2024-09-10_CD.root");
tree_both = Chain("PhysicsTree",files,true);
+
+ files.clear();
+//files.push_back("~/Programs/nptoolv3/Outputs/Analysis/19Odp_XYZT_0003_CD.root");
+//files.push_back("~/Programs/nptoolv3/Outputs/Analysis/19Odp_2024-08-28_CD.root");
+//files.push_back("~/Programs/nptoolv3/Outputs/Analysis/19Odp_2024-09-05_CD.root");
+//files.push_back("~/Programs/nptoolv3/Outputs/Analysis/19Odp_2024-09-06_CD.root");
+ files.push_back("~/Programs/nptoolv3/Outputs/Analysis/19Odp_2024-09-10_CD.root");
+ tree_cd = Chain("PhysicsTree",files,true);
+
+ files.clear();
+//files.push_back("~/Programs/nptoolv3/Outputs/Analysis/19Odp_XYZT_0003_Au.root");
+//files.push_back("~/Programs/nptoolv3/Outputs/Analysis/19Odp_2024-08-28_Au.root");
+//files.push_back("~/Programs/nptoolv3/Outputs/Analysis/19Odp_2024-09-05_Au.root");
+//files.push_back("~/Programs/nptoolv3/Outputs/Analysis/19Odp_2024-09-06_Au.root");
+ files.push_back("~/Programs/nptoolv3/Outputs/Analysis/19Odp_2024-09-10_Au.root");
+ tree_au = Chain("PhysicsTree",files,true);
}
//////////////////////////////////////////////////////////////////////////////////////
@@ -123,48 +138,80 @@ void LoadChain_PhaseSpace()
}
//////////////////////////////////////////////////////////////////////////////////////
-void DrawPlotsToLoad(TChain* tree)
+void DrawPlotsToLoad_ExEg_AddBack(TChain* tree)
{
- // ExEg
- tree->Draw("Ex:AddBack_EDC>>hExEg(8000,0,8,2000,-5,15)","abs(T_FPMW_HF-13000)<1000 && ThetaLab<158","colz");
+ TCanvas *c = new TCanvas("cExEg_AddBack","cExEg_AddBack",500,500);
+ tree->Draw("Ex:AddBack_EDC>>hExEg(8000,0,8,2000,-5,15)","abs(T_FPMW_HF-13300)<1500 && ThetaLab<158","colz");
TH2F* hExEg = (TH2F*) gDirectory->Get("hExEg");
+ if( tree->GetEntries() == tree_both->GetEntries() ) { hExEg->SetName("hExEg_both"); }
+ else if( tree->GetEntries() == tree_cd->GetEntries() ) { hExEg->SetName("hExEg_cd"); }
+ else if( tree->GetEntries() == tree_au->GetEntries() ) { hExEg->SetName("hExEg_au"); }
+}
- // ExEt
- tree->Draw("Ex:trackEDC>>hExEt(8000,0,8,2000,-5,15)","abs(T_FPMW_HF-13000)<1000 && ThetaLab<158","colz");
+//////////////////////////////////////////////////////////////////////////////////////
+void DrawPlotsToLoad_ExEg_Tracked(TChain* tree)
+{
+ TCanvas *c = new TCanvas("cExEg_Tracked","cExEg_Tracked",500,500);
+ tree->Draw("Ex:trackEDC>>hExEt(8000,0,8,2000,-5,15)","abs(T_FPMW_HF-13300)<1500 && ThetaLab<158","colz");
TH2F* hExEt= (TH2F*) gDirectory->Get("hExEt");
+ if( tree->GetEntries() == tree_both->GetEntries() ) { hExEt->SetName("hExEt_both"); }
+ else if( tree->GetEntries() == tree_cd->GetEntries() ) { hExEt->SetName("hExEt_cd"); }
+ else if( tree->GetEntries() == tree_au->GetEntries() ) { hExEt->SetName("hExEt_au"); }
+}
- // Ex (cleanup)
- tree->Draw("Ex>>hExGood(2000,-5,15)","abs(T_FPMW_HF-13000)<1000 && ThetaLab<158","colz");
+//////////////////////////////////////////////////////////////////////////////////////
+void DrawPlotsToLoad_Ex(TChain* tree)
+{
+ TCanvas *c = new TCanvas("cEx","cEx",500,500);
+ tree->Draw("Ex>>hExGood(2000,-5,15)","abs(T_FPMW_HF-13300)<1500 && 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");
+ tree->Draw("Ex>>hExBad(2000,-5,15)","abs(T_FPMW_HF-10500)<1500 && 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");
+ if( tree->GetEntries() == tree_both->GetEntries() ) { hExClean->SetName("hExClean_both"); }
+ else if( tree->GetEntries() == tree_cd->GetEntries() ) { hExClean->SetName("hExClean_cd"); }
+ else if( tree->GetEntries() == tree_au->GetEntries() ) { hExClean->SetName("hExClean_au"); }
+}
+
+//////////////////////////////////////////////////////////////////////////////////////
+void DrawPlotsToLoad_ExThetaLab(TChain* tree)
+{
+ TCanvas *c = new TCanvas("cExThetaLab","cExThetaLab",500,500);
+ tree->Draw("Ex:ThetaLab>>hExTLGood(80,100,180,2000,-5,15)","abs(T_FPMW_HF-13300)<1500 && 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");
+ tree->Draw("Ex:ThetaLab>>hExTLBad(80,100,180,2000,-5,15)","abs(T_FPMW_HF-10500)<1500 && 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();
-
+ if( tree->GetEntries() == tree_both->GetEntries() ) { hExTLClean->SetName("hExTLClean_both"); }
+ else if( tree->GetEntries() == tree_cd->GetEntries() ) { hExTLClean->SetName("hExTLClean_cd"); }
+ else if( tree->GetEntries() == tree_au->GetEntries() ) { hExTLClean->SetName("hExTLClean_au"); }
}
//////////////////////////////////////////////////////////////////////////////////////
-void DrawPlotsToLoad()
+void DrawPlotsToLoad_All(TChain* tree)
{
- DrawPlotsToLoad(tree_both);
-
+ DrawPlotsToLoad_ExEg_AddBack(tree);
+ DrawPlotsToLoad_ExEg_Tracked(tree);
+ DrawPlotsToLoad_Ex(tree);
+ DrawPlotsToLoad_ExThetaLab(tree);
}
-//=====================================================================//
-//=====================================================================//
+
+//////////////////////////////////////////////////////////////////////////////////////
+void DrawPlotsToLoad_ExEg_AddBack(){ DrawPlotsToLoad_ExEg_AddBack(tree_both); }
+void DrawPlotsToLoad_ExEg_Tracked(){ DrawPlotsToLoad_ExEg_Tracked(tree_both); }
+void DrawPlotsToLoad_Ex(){ DrawPlotsToLoad_Ex(tree_both); }
+void DrawPlotsToLoad_ExThetaLab(){ DrawPlotsToLoad_ExThetaLab(tree_both); }
+void DrawPlotsToLoad_All(){ DrawPlotsToLoad_All(tree_both); }
+
+//==================================================================================//
//////////////////////////////////////////////////////////////////////////////////////
TH1F* Project(TH2F* h2D, const char * XY, Color_t col)
@@ -185,7 +232,7 @@ TH1F* Project(TH2F* h2D, const char * XY, Color_t col)
return 0;
}
return (TH1F*) h;
-};
+}
//////////////////////////////////////////////////////////////////////////////////////
TH1F* Project(TH2F* h2D, const char * XY, double mean, double width, Color_t col)
@@ -210,7 +257,7 @@ TH1F* Project(TH2F* h2D, const char * XY, double mean, double width, Color_t col
return 0;
}
return (TH1F*) h;
-};
+}
//////////////////////////////////////////////////////////////////////////////////////
void Prj_DrawLines(double mean, double width, Color_t col)
@@ -232,14 +279,30 @@ string AddOrTrack(string input)
}
//////////////////////////////////////////////////////////////////////////////////////
-void GateGamma_SeeParticle(string addbackOrTrack, double bing, double binp, double Eg1, double width1)
+void CheckExEgExists(string histname, string tree, string addbackOrTrack)
{
- 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();
+ if(addbackOrTrack == "track") {
+ if(tree == "cd") { DrawPlotsToLoad_ExEg_Tracked(tree_cd); }
+ if(tree == "au") { DrawPlotsToLoad_ExEg_Tracked(tree_au); }
+ if(tree == "both") { DrawPlotsToLoad_ExEg_Tracked(tree_both); }
+ } else {
+ if(tree == "cd") { DrawPlotsToLoad_ExEg_AddBack(tree_cd); }
+ if(tree == "au") { DrawPlotsToLoad_ExEg_AddBack(tree_au); }
+ if(tree == "both") { DrawPlotsToLoad_ExEg_AddBack(tree_both); }
+ }
}
+}
+
+//==================================================================================//
+
+//////////////////////////////////////////////////////////////////////////////////////
+void GateGamma_SeeParticle(string treesuffix, string addbackOrTrack, double bing, double binp,
+ double Eg1, double width1)
+{
+ string histname = AddOrTrack(addbackOrTrack);
+ histname = histname + "_" + treesuffix;
+ CheckExEgExists(histname, treesuffix, addbackOrTrack);
//Load in & rebin
TH2F* hCopyMe = (TH2F*) gDirectory->Get(histname.c_str());
@@ -266,15 +329,13 @@ void GateGamma_SeeParticle(string addbackOrTrack, double bing, double binp, doub
}
//////////////////////////////////////////////////////////////////////////////////////
-void GateGamma_SeeParticle(string addbackOrTrack, double bing, double binp, double Eg1, double width1,
+void GateGamma_SeeParticle(string treesuffix, 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();
- }
+ histname = histname + "_" + treesuffix;
+ CheckExEgExists(histname, treesuffix, addbackOrTrack);
//Load in & rebin
TH2F* hCopyMe = (TH2F*) gDirectory->Get(histname.c_str());
@@ -304,16 +365,14 @@ void GateGamma_SeeParticle(string addbackOrTrack, double bing, double binp, doub
}
//////////////////////////////////////////////////////////////////////////////////////
-void GateGamma_SeeParticle(string addbackOrTrack, double bing, double binp, double Eg1, double width1,
+void GateGamma_SeeParticle(string treesuffix, 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();
- }
+ histname = histname + "_" + treesuffix;
+ CheckExEgExists(histname, treesuffix, addbackOrTrack);
//Load in & rebin
TH2F* hCopyMe = (TH2F*) gDirectory->Get(histname.c_str());
@@ -346,17 +405,15 @@ void GateGamma_SeeParticle(string addbackOrTrack, double bing, double binp, doub
}
//////////////////////////////////////////////////////////////////////////////////////
-void GateGamma_SeeParticle(string addbackOrTrack, double bing, double binp, double Eg1, double width1,
+void GateGamma_SeeParticle(string treesuffix, 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();
- }
+ histname = histname + "_" + treesuffix;
+ CheckExEgExists(histname, treesuffix, addbackOrTrack);
//Load in & rebin
TH2F* hCopyMe = (TH2F*) gDirectory->Get(histname.c_str());
@@ -392,18 +449,16 @@ void GateGamma_SeeParticle(string addbackOrTrack, double bing, double binp, doub
}
//////////////////////////////////////////////////////////////////////////////////////
-void GateGamma_SeeParticle(string addbackOrTrack, double bing, double binp, double Eg1, double width1,
+void GateGamma_SeeParticle(string treesuffix, 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();
- }
+ histname = histname + "_" + treesuffix;
+ CheckExEgExists(histname, treesuffix, addbackOrTrack);
//Load in & rebin
TH2F* hCopyMe = (TH2F*) gDirectory->Get(histname.c_str());
@@ -442,14 +497,12 @@ void GateGamma_SeeParticle(string addbackOrTrack, double bing, double binp, doub
}
//////////////////////////////////////////////////////////////////////////////////////
-void GateParticle_SeeGamma(string addbackOrTrack, double bing, double binp, double Ex1, double width1)
+void GateParticle_SeeGamma(string treesuffix, 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();
- }
+ histname = histname + "_" + treesuffix;
+ CheckExEgExists(histname, treesuffix, addbackOrTrack);
//Load in & rebin
TH2F* hCopyMe = (TH2F*) gDirectory->Get(histname.c_str());
@@ -476,15 +529,13 @@ void GateParticle_SeeGamma(string addbackOrTrack, double bing, double binp, doub
}
//////////////////////////////////////////////////////////////////////////////////////
-void GateParticle_SeeGamma(string addbackOrTrack, double bing, double binp, double Ex1, double width1,
+void GateParticle_SeeGamma(string treesuffix, 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();
- }
+ histname = histname + "_" + treesuffix;
+ CheckExEgExists(histname, treesuffix, addbackOrTrack);
//Load in & rebin
TH2F* hCopyMe = (TH2F*) gDirectory->Get(histname.c_str());
@@ -514,16 +565,14 @@ void GateParticle_SeeGamma(string addbackOrTrack, double bing, double binp, doub
}
//////////////////////////////////////////////////////////////////////////////////////
-void GateParticle_SeeGamma(string addbackOrTrack, double bing, double binp, double Ex1, double width1,
+void GateParticle_SeeGamma(string treesuffix, 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();
- }
+ histname = histname + "_" + treesuffix;
+ CheckExEgExists(histname, treesuffix, addbackOrTrack);
//Load in & rebin
TH2F* hCopyMe = (TH2F*) gDirectory->Get(histname.c_str());
@@ -556,17 +605,15 @@ void GateParticle_SeeGamma(string addbackOrTrack, double bing, double binp, doub
}
//////////////////////////////////////////////////////////////////////////////////////
-void GateParticle_SeeGamma(string addbackOrTrack, double bing, double binp, double Ex1, double width1,
+void GateParticle_SeeGamma(string treesuffix, 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();
- }
+ histname = histname + "_" + treesuffix;
+ CheckExEgExists(histname, treesuffix, addbackOrTrack);
//Load in & rebin
TH2F* hCopyMe = (TH2F*) gDirectory->Get(histname.c_str());
@@ -602,18 +649,16 @@ void GateParticle_SeeGamma(string addbackOrTrack, double bing, double binp, doub
}
//////////////////////////////////////////////////////////////////////////////////////
-void GateParticle_SeeGamma(string addbackOrTrack, double bing, double binp, double Ex1, double width1,
+void GateParticle_SeeGamma(string treesuffix, 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();
- }
+ histname = histname + "_" + treesuffix;
+ CheckExEgExists(histname, treesuffix, addbackOrTrack);
//Load in & rebin
TH2F* hCopyMe = (TH2F*) gDirectory->Get(histname.c_str());
@@ -651,18 +696,120 @@ void GateParticle_SeeGamma(string addbackOrTrack, double bing, double binp, doub
projG5->Draw("same");
}
+//==================================================================================//
+
+//////////////////////////////////////////////////////////////////////////////////////
+void GateParticle_CompareTargets(string addbackOrTrack, double bing, double binp,
+ double Ex1, double width1)
+{
+ string histname_cd = AddOrTrack(addbackOrTrack); histname_cd += "_cd";
+ CheckExEgExists(histname_cd, "cd", addbackOrTrack);
+
+ string histname_au = AddOrTrack(addbackOrTrack); histname_au += "_au";
+ CheckExEgExists(histname_au, "au", addbackOrTrack);
+
+ //Load in & rebin CD
+ TH2F* hCopyMe = (TH2F*) gDirectory->Get(histname_cd.c_str());
+ TH2F* hExEg_cd = (TH2F*) hCopyMe->Clone();
+ hExEg_cd->Rebin2D(bing/hExEg_cd->GetXaxis()->GetBinWidth(3),
+ binp/hExEg_cd->GetYaxis()->GetBinWidth(3));
+ //Load in & rebin Au
+ hCopyMe = (TH2F*) gDirectory->Get(histname_au.c_str());
+ TH2F* hExEg_au = (TH2F*) hCopyMe->Clone();
+ hExEg_au->Rebin2D(bing/hExEg_au->GetXaxis()->GetBinWidth(3),
+ binp/hExEg_au->GetYaxis()->GetBinWidth(3));
+ //Set up canvas
+ TCanvas *cGPCT = new TCanvas("cGPCT","cGPCT",1000,750);
+ cGPCT->Divide(1,2);
+ cGPCT->cd(1);
+
+ //Axis being projected
+ TH1F* projP_au = Project(hExEg_au, "Y", kOrange+2);
+ projP_au->Draw();
+ TH1F* projP_cd = Project(hExEg_cd, "Y", kBlack);
+ projP_cd->SetFillColor(kBlack); projP_cd->SetFillStyle(3001);
+ projP_cd->Draw("same");
+
+ //Projection gates
+ Prj_DrawLines(Ex1, width1, kRed);
+
+ //Projected axis
+ cGPCT->cd(2);
+ TH1F* projG1_au = Project(hExEg_au, "X", Ex1, width1, kRed);
+ TH1F* projG1_cd = Project(hExEg_cd, "X", Ex1, width1, kRed+2);
+ projG1_cd->SetFillColor(kRed+2); projG1_cd->SetFillStyle(3001);
+
+ projG1_au->Draw();
+ projG1_cd->Draw("same");
+}
+
+//////////////////////////////////////////////////////////////////////////////////////
+void GateParticle_CompareTargets(string addbackOrTrack, double bing, double binp,
+ double Ex1, double width1,
+ double Ex2, double width2)
+{
+ string histname_cd = AddOrTrack(addbackOrTrack); histname_cd += "_cd";
+ CheckExEgExists(histname_cd, "cd", addbackOrTrack);
+
+ string histname_au = AddOrTrack(addbackOrTrack); histname_au += "_au";
+ CheckExEgExists(histname_au, "au", addbackOrTrack);
+
+ //Load in & rebin CD
+ TH2F* hCopyMe = (TH2F*) gDirectory->Get(histname_cd.c_str());
+ TH2F* hExEg_cd = (TH2F*) hCopyMe->Clone();
+ hExEg_cd->Rebin2D(bing/hExEg_cd->GetXaxis()->GetBinWidth(3),
+ binp/hExEg_cd->GetYaxis()->GetBinWidth(3));
+
+ //Load in & rebin Au
+ hCopyMe = (TH2F*) gDirectory->Get(histname_au.c_str());
+ TH2F* hExEg_au = (TH2F*) hCopyMe->Clone();
+ hExEg_au->Rebin2D(bing/hExEg_au->GetXaxis()->GetBinWidth(3),
+ binp/hExEg_au->GetYaxis()->GetBinWidth(3));
+
+ //Set up canvas
+ TCanvas *cGPCT = new TCanvas("cGPCT","cGPCT",1000,750);
+ cGPCT->Divide(1,2);
+ cGPCT->cd(1);
+
+ //Axis being projected
+ TH1F* projP_au = Project(hExEg_au, "Y", kOrange+2);
+ projP_au->Draw();
+ TH1F* projP_cd = Project(hExEg_cd, "Y", kBlack);
+ projP_cd->SetFillColor(kBlack); projP_cd->SetFillStyle(3001);
+ projP_cd->Draw("same");
+
+ //Projection gates
+ Prj_DrawLines(Ex1, width1, kRed);
+ Prj_DrawLines(Ex2, width2, kBlue);
+
+ //Projected axis
+ cGPCT->cd(2);
+ TH1F* projG1_au = Project(hExEg_au, "X", Ex1, width1, kRed);
+ TH1F* projG1_cd = Project(hExEg_cd, "X", Ex1, width1, kRed+2);
+ projG1_cd->SetFillColor(kRed+2); projG1_cd->SetFillStyle(3001);
+ TH1F* projG2_au = Project(hExEg_au, "X", Ex2, width2, kBlue);
+ TH1F* projG2_cd = Project(hExEg_cd, "X", Ex2, width2, kBlue+2);
+ projG2_cd->SetFillColor(kBlue+2); projG2_cd->SetFillStyle(3001);
+
+ projG1_au->Draw();
+ projG1_cd->Draw("same");
+ projG2_au->Draw("same");
+ projG2_cd->Draw("same");
+}
+
+
+
+//==================================================================================//
//////////////////////////////////////////////////////////////////////////////////////
-void GateGamma_SeeParticle_HighEnergyStaggered(string addbackOrTrack, double bing, double binp)
+void GateGamma_SeeParticle_HighEnergyStaggered(string treesuffix, 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();
- }
+ histname = histname + "_" + treesuffix;
+ CheckExEgExists(histname, treesuffix, addbackOrTrack);
gStyle->SetPalette(kCool);
@@ -703,24 +850,20 @@ void GateGamma_SeeParticle_HighEnergyStaggered(string addbackOrTrack, double bin
}
+//==================================================================================//
-
-
-//=====================================================================//
-//=====================================================================//
-
-////////////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////////////
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");
+ tree->Draw("Ex:ThetaLab>>hGood(26,104,156,150,-2,13)","abs(T_FPMW_HF-13300)<1500 && 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");
+ tree->Draw("Ex:ThetaLab>>hBad(26,104,156,150,-2,13)","abs(T_FPMW_HF-10500)<1500 && Mugast.TelescopeNumber<10","colz");
TH2F* hBad = (TH2F*) gDirectory->Get("hBad");
// ExTL Clean
@@ -811,18 +954,18 @@ void DrawExThetaLab_SubPSp(TChain* tree)
}
-////////////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////////////
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");
+ tree->Draw("Ex>>hGood(2000,-5,15)","abs(T_FPMW_HF-13300)<1500 && 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");
+ tree->Draw("Ex>>hBad(2000,-5,15)","abs(T_FPMW_HF-10500)<1500 && Mugast.TelescopeNumber<10","colz");
TH1F* hBad = (TH1F*) gDirectory->Get("hBad");
// ExTL Clean
@@ -832,7 +975,7 @@ void DrawEx_SubPSp(TChain* tree)
hClean->Add(hBad,-1);
// ExTL PSp
- tree_PSp00->Draw("Ex>>hPSp(300,-2,13)","","colz");
+ tree_PSp00->Draw("Ex>>hPSp(2000,-5,15)","","colz");
TH1F* hPSp = (TH1F*) gDirectory->Get("hPSp");
double scale = 0.134256;
@@ -849,7 +992,7 @@ void DrawEx_SubPSp(TChain* tree)
}
-/////////////////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////////////
void AddKnownStateLines(TH1F* hist, double ymax)
{
hist->Draw();
@@ -874,7 +1017,7 @@ 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",
+ "Mugast.TelescopeNumber<9 && abs(T_FPMW_HF-13300)<1500",
"");
TH1F* hTracked = (TH1F*) gDirectory->Get("hTracked");
hTracked->SetLineColor(kRed);
@@ -882,7 +1025,7 @@ void DrawGammas_AddBackVersusTracked()
tree_both->Draw("AddBack_EDC>>hAddBack(1000,-1,9)",
- "Mugast.TelescopeNumber<9 && abs(T_FPMW_HF-13000)<1000",
+ "Mugast.TelescopeNumber<9 && abs(T_FPMW_HF-13300)<1500",
"");
TH1F* hAddBack = (TH1F*) gDirectory->Get("hAddBack");
hAddBack->SetLineColor(kBlue);
@@ -892,11 +1035,9 @@ void DrawGammas_AddBackVersusTracked()
}
+//==================================================================================//
-//=====================================================================//
-//=====================================================================//
-
-////////////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////////////
void Diagnostics_StripNumber(TTree* tree)
{
TCanvas* c = new TCanvas("c","c",1000,1000);
@@ -942,23 +1083,22 @@ void Diagnostics_StripNumber(TTree* tree)
}
-////////////////////////////////////////////////////////////////////////////////
-void GateThetaLab_MultiWrite(int numGates, double binsize)
+//////////////////////////////////////////////////////////////////////////////////////
+void GateThetaLab_MultiWrite(double startTheta, int numGates, double gatesize, double ExBinsize)
{
- double startTheta= 104.;
- double finishTheta=156.;
+ double finishTheta = startTheta + (numGates * gatesize);
string drawbase = to_string(numGates)
- + ",104,156,"
- + to_string(20.0/binsize)
+ + "," + to_string(startTheta) + "," + to_string(finishTheta) + ","
+ + to_string(20.0/ExBinsize)
+ ",-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");
+ tree_both->Draw(drawgood.c_str(),"abs(T_FPMW_HF-13300)<1500","colz");
TH2F* hGood = (TH2F*) gDirectory->Get("hGood");
- tree_both->Draw(drawbad.c_str(), "abs(T_FPMW_HF-10200)<600","colz");
+ tree_both->Draw(drawbad.c_str(), "abs(T_FPMW_HF-10500)<1500","colz");
TH2F* hBad = (TH2F*) gDirectory->Get("hBad");
TH2F* hClean = (TH2F*) hGood->Clone();
@@ -967,8 +1107,7 @@ void GateThetaLab_MultiWrite(int numGates, double binsize)
hClean->Add(hBad,-1);
hClean->Draw("colz");
- string ytitle = "Counts / " + to_string(binsize) + " MeV";
- double gatesize = (finishTheta-startTheta)/numGates;
+ string ytitle = "Counts / " + to_string(ExBinsize) + " MeV";
TList* list = new TList();
for (int i=0; i<numGates; i++){
@@ -989,24 +1128,78 @@ void GateThetaLab_MultiWrite(int numGates, double binsize)
file->ls();
}
-////////////////////////////////////////////////////////////////////////////////
-void GatePhaseSpace_MultiWrite(int numGates, double binsize)
+//////////////////////////////////////////////////////////////////////////////////////
+void GateThetaLab_MultiWrite_SubPSp(double startTheta, int numGates, double gatesize, double ExBinsize)
{
- double startTheta= 104.;
- double finishTheta=156.;
+ double finishTheta = startTheta + (numGates * gatesize);
+
+ LoadChain_PhaseSpace();
+ TChain* treePSp = tree_PSp00;
+
+ string drawbase = to_string(numGates)
+ + "," + to_string(startTheta) + "," + to_string(finishTheta) + ","
+ + to_string(20.0/ExBinsize)
+ + ",-5,15)";
+ string drawgood = "Ex:ThetaLab>>hGood(" + drawbase;
+ string drawbad = "Ex:ThetaLab>>hBad(" + drawbase;
+ string drawpsp = "Ex:ThetaLab>>hPSp(" + drawbase;
+
+ tree_both->Draw(drawgood.c_str(),"abs(T_FPMW_HF-13300)<1500","colz");
+ TH2F* hGood = (TH2F*) gDirectory->Get("hGood");
+
+ tree_both->Draw(drawbad.c_str(), "abs(T_FPMW_HF-10500)<1500","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");
+
+ treePSp->Draw(drawpsp.c_str(),"","colz");
+ TH2F* hPSp = (TH2F*) gDirectory->Get("hPSp"); hPSp->SetTitle("hPSp"); hPSp->SetName("hPSp");
+ double scale = 0.134256; cout << GREEN << " USING SCALING FROM DrawEx_SubPSp, SCALE = " << scale << RESET << endl;
+
+ hClean->Add(hPSp, -scale);
+
+ string ytitle = "Counts / " + to_string(ExBinsize) + " MeV";
+ 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(double startTheta, int numGates, double gatesize, double ExBinsize)
+{
+ LoadChain_PhaseSpace();
+ double finishTheta = startTheta + (numGates * gatesize);
string draw = "Ex:ThetaLab>>h("
+ to_string(numGates)
- + ",104,156,"
- + to_string(20.0/binsize)
+ + "," + to_string(startTheta) + "," + to_string(finishTheta) + ","
+ + to_string(20.0/ExBinsize)
+ ",-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;
+ string ytitle = "Counts / " + to_string(ExBinsize) + " MeV";
TList* list = new TList();
for (int i=0; i<numGates; i++){
@@ -1027,8 +1220,135 @@ void GatePhaseSpace_MultiWrite(int numGates, double binsize)
file->ls();
}
-//=====================================================================//
-//=====================================================================//
+
+
+//////////////////////////////////////////////////////////////////////////////////////
+void GateThetaLab_MultiWrite_ParallelCuts(double startTheta, int numGates, double gatesize, double ExBinsize)
+{
+ double finishTheta = startTheta + (numGates * gatesize);
+
+ string baseParallelCuts = "(ELab < (0.05458*ThetaLab)-3.456) && (ELab > (0.05458*ThetaLab)-5.696) && Ex < 8.7";
+
+ string drawbase = to_string(numGates)
+ + "," + to_string(startTheta) + "," + to_string(finishTheta) + ","
+ + to_string(20.0/ExBinsize)
+ + ",-5,15)";
+ string drawgood = "Ex:ThetaLab>>hGood(" + drawbase;
+ string drawbad = "Ex:ThetaLab>>hBad(" + drawbase;
+ string drawpsp = "Ex:ThetaLab>>hPSp(" + drawbase;
+
+ string gategood = "abs(T_FPMW_HF-13300)<1500 &&" + baseParallelCuts;
+ string gatebad = "abs(T_FPMW_HF-10500)<1500 &&" + baseParallelCuts;
+
+ tree_both->Draw(drawgood.c_str(), gategood.c_str(), "colz");
+ TH2F* hGood = (TH2F*) gDirectory->Get("hGood");
+
+ tree_both->Draw(drawbad.c_str() , gatebad.c_str() , "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(ExBinsize) + " MeV";
+ 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_Parallel_Experiment.root","RECREATE");
+ list->Write("GateThetaLabHistograms",TObject::kSingleKey);
+ file->ls();
+}
+
+
+
+//////////////////////////////////////////////////////////////////////////////////////
+void GateThetaLab_MultiWrite_ParallelCuts_SubPSp(double startTheta, int numGates, double gatesize, double ExBinsize)
+{
+ double finishTheta = startTheta + (numGates * gatesize);
+
+ string baseParallelCuts = "(ELab < (0.05458*ThetaLab)-3.456) && (ELab > (0.05458*ThetaLab)-5.696) && Ex < 8.7";
+
+ LoadChain_PhaseSpace();
+ TChain* treePSp = tree_PSp00;
+
+ string drawbase = to_string(numGates)
+ + "," + to_string(startTheta) + "," + to_string(finishTheta) + ","
+ + to_string(20.0/ExBinsize)
+ + ",-5,15)";
+ string drawgood = "Ex:ThetaLab>>hGood(" + drawbase;
+ string drawbad = "Ex:ThetaLab>>hBad(" + drawbase;
+ string drawpsp = "Ex:ThetaLab>>hPSp(" + drawbase;
+
+ string gategood = "abs(T_FPMW_HF-13300)<1500 &&" + baseParallelCuts;
+ string gatebad = "abs(T_FPMW_HF-10500)<1500 &&" + baseParallelCuts;
+
+ tree_both->Draw(drawgood.c_str(), gategood.c_str(), "colz");
+ TH2F* hGood = (TH2F*) gDirectory->Get("hGood");
+
+ tree_both->Draw(drawbad.c_str() , gatebad.c_str() , "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");
+
+ treePSp->Draw(drawpsp.c_str(),baseParallelCuts.c_str(),"colz");
+ TH2F* hPSp = (TH2F*) gDirectory->Get("hPSp"); hPSp->SetTitle("hPSp"); hPSp->SetName("hPSp");
+ double scale = 1.0;
+ for(int t=0; t<numGates; t++){
+ for(int b=hClean->GetYaxis()->FindBin(7.7); b<hClean->GetYaxis()->FindBin(8.5);b++){
+ while(hClean->GetBinContent(t,b)>2 && hClean->GetBinContent(t,b)<(hPSp->GetBinContent(t,b)*scale)){
+ scale *= 0.999;
+ }
+ }
+ }
+ cout << "scale = " << scale << endl;
+
+ hClean->Add(hPSp,-scale);
+
+ string ytitle = "Counts / " + to_string(ExBinsize) + " MeV";
+ 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_Parallel_Experiment_SubPSp.root","RECREATE");
+ list->Write("GateThetaLabHistograms",TObject::kSingleKey);
+ file->ls();
+}
+
+
+
+
+
+
+//==================================================================================//
void TWOFNR_ApproachingUnbound()
{
@@ -1040,44 +1360,135 @@ void TWOFNR_ApproachingUnbound()
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");
+ TGraph* s12_7500 = TWOFNR(7.500,2,1,0,0.5,2.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,2.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,2.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,2.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,2.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,2.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");
+ TGraph* p12_7500 = TWOFNR(7.500,2,0,1,0.5,2.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,2.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,2.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,2.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,2.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,2.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");
+ TGraph* d52_7500 = TWOFNR(7.500,2,0,2,2.5,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,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,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,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,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,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");
+ TGraph* f72_7500 = TWOFNR(7.500,3,0,3,3.5,2.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,2.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,2.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,2.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,2.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,2.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");
+}
+void TWOFNR_1MeVJumps()
+{
+ 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(2.000,2,1,0,0.5,2.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(3.000,2,1,0,0.5,2.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(4.000,2,1,0,0.5,2.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(5.000,2,1,0,0.5,2.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(6.000,2,1,0,0.5,2.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.000,2,1,0,0.5,2.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(2.000,2,0,1,0.5,2.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(3.000,2,0,1,0.5,2.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(4.000,2,0,1,0.5,2.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(5.000,2,0,1,0.5,2.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(6.000,2,0,1,0.5,2.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.000,2,0,1,0.5,2.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(2.000,2,0,2,2.5,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(3.000,2,0,2,2.5,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(4.000,2,0,2,2.5,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(5.000,2,0,2,2.5,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(6.000,2,0,2,2.5,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.000,2,0,2,2.5,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(2.000,3,0,3,3.5,2.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(3.000,3,0,3,3.5,2.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(4.000,3,0,3,3.5,2.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(5.000,3,0,3,3.5,2.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(6.000,3,0,3,3.5,2.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.000,3,0,3,3.5,2.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");
+cout << "FROM ENERGY 2 MeV TO 7 MeV" << endl;
+}
+//==================================================================================//
+double AGATAefficiency(double gammaE_MeV){
+ double A = 1.0;
+ double B = 1.0;
+ double C = 0.000; //Fixed to 0, hence not included in IZanon's thesis equation
+ double D = 1.0; //IZ: C
+ double E = 1.0; //IZ: D
+ double F = 1.0; //IZ: E
+ double G = 1.0; //IZ: F
+
+ double constX = log(gammaE_MeV / 0.1); //x = log(Eg/100 keV)
+ double constY = log(gammaE_MeV / 1.0); //x = log(Eg/ 1 MeV)
+
+ double exponent1 = pow(A + (B * constX) + (C * constX * constX) , -G);
+ double exponent2 = pow(D + (E * constY) + (F * constY * constY) , -G);
+
+ double exponentF = pow(exponent1+exponent2, -1/G);
+
+ double effic = exp(exponentF);
+
+
+ cout << " DO NOT HAVE THE VARIABLED FO THE FUNCTIONS! INVALID! " << endl;
+ cout << " DO NOT HAVE THE VARIABLED FO THE FUNCTIONS! INVALID! " << endl;
+ cout << " DO NOT HAVE THE VARIABLED FO THE FUNCTIONS! INVALID! " << endl;
+ cout << " DO NOT HAVE THE VARIABLED FO THE FUNCTIONS! INVALID! " << endl;
+ cout << " DO NOT HAVE THE VARIABLED FO THE FUNCTIONS! INVALID! " << endl;
+ cout << " DO NOT HAVE THE VARIABLED FO THE FUNCTIONS! INVALID! " << endl;
+ cout << " DO NOT HAVE THE VARIABLED FO THE FUNCTIONS! INVALID! " << endl;
+ cout << " DO NOT HAVE THE VARIABLED FO THE FUNCTIONS! INVALID! " << endl;
+ cout << " DO NOT HAVE THE VARIABLED FO THE FUNCTIONS! INVALID! " << endl;
+ cout << " DO NOT HAVE THE VARIABLED FO THE FUNCTIONS! INVALID! " << endl;
+ cout << " DO NOT HAVE THE VARIABLED FO THE FUNCTIONS! INVALID! " << endl;
+ cout << " DO NOT HAVE THE VARIABLED FO THE FUNCTIONS! INVALID! " << endl;
+ cout << " DO NOT HAVE THE VARIABLED FO THE FUNCTIONS! INVALID! " << endl;
+ cout << " DO NOT HAVE THE VARIABLED FO THE FUNCTIONS! INVALID! " << endl;
+ cout << " DO NOT HAVE THE VARIABLED FO THE FUNCTIONS! INVALID! " << endl;
+ cout << " DO NOT HAVE THE VARIABLED FO THE FUNCTIONS! INVALID! " << endl;
+ cout << " DO NOT HAVE THE VARIABLED FO THE FUNCTIONS! INVALID! " << endl;
+ cout << " DO NOT HAVE THE VARIABLED FO THE FUNCTIONS! INVALID! " << endl;
+ cout << " DO NOT HAVE THE VARIABLED FO THE FUNCTIONS! INVALID! " << endl;
+ cout << " DO NOT HAVE THE VARIABLED FO THE FUNCTIONS! INVALID! " << endl;
+ cout << " DO NOT HAVE THE VARIABLED FO THE FUNCTIONS! INVALID! " << endl;
+ cout << " DO NOT HAVE THE VARIABLED FO THE FUNCTIONS! INVALID! " << endl;
+ cout << " DO NOT HAVE THE VARIABLED FO THE FUNCTIONS! INVALID! " << endl;
+ cout << " DO NOT HAVE THE VARIABLED FO THE FUNCTIONS! INVALID! " << endl;
+ cout << " DO NOT HAVE THE VARIABLED FO THE FUNCTIONS! INVALID! " << endl;
+ cout << " DO NOT HAVE THE VARIABLED FO THE FUNCTIONS! INVALID! " << endl;
+ cout << " DO NOT HAVE THE VARIABLED FO THE FUNCTIONS! INVALID! " << endl;
+ cout << " DO NOT HAVE THE VARIABLED FO THE FUNCTIONS! INVALID! " << endl;
+
+
+ return effic;
}
-
diff --git a/Projects/e775s/analysis/Plots_KnownPeakFitter.h b/Projects/e775s/analysis/Plots_KnownPeakFitter.h
index 4ba99d84c44bcc9ef6cc341a1d5a9aea8882945e..fff07484777a9f4215b3abb1011ef316118abb76 100644
--- a/Projects/e775s/analysis/Plots_KnownPeakFitter.h
+++ b/Projects/e775s/analysis/Plots_KnownPeakFitter.h
@@ -5,38 +5,50 @@
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<double> meansGauss = { 0.000,
+ 1.675,
+ 3.572,
+ 4.071,
+ 4.456,
+ 5.004,
+ 5.228,
+ 5.629,
+ 7.252,
+ 7.622,
+ 7.855,
+ 8.313,
+ //8.561, //MOVED TO VOIGT
+ 8.962
+ };
+
+vector<double> meansVoigt = { 8.561
+ };
+
+vector<double> means; //meansGauss + meansVoigt, added together in Plots_Experiment()
+//means.reserve(meansGauss.size() + meansVoigt.size());
+//means.insert(means.end(), meansGauss.begin(),meansGauss.end());
+//means.insert(means.end(), meansVoigt.begin(),meansVoigt.end());
+
+//-------------------------------------------------------
+
+int numGauss;// = 13;
+int numVoigt;// = 1;
+int numPeaks;// = numGauss + numVoigt;
+
+double GetSigmaFromFunction(double Ex){ return (0.000969*Ex*Ex - 0.009116*Ex + 0.128396); }
vector<vector<double>> FitKnownPeaks_RtrnArray(TH1F* hist){
//double minfit = -1.0, maxfit = 7.0;
- double minfit = -1.0, maxfit = 10.0;
+ double minfit = -1.0, maxfit = 8.7;
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];
TF1 **allPeaks = new TF1*[numPeaks];
- for(int i=0; i<numPeaks; i++){
+ for(int i=0; i<numGauss; i++){
string nameHere = nameBase;
nameHere += to_string(i);
allPeaks[i] = new TF1(nameHere.c_str(), function.c_str(), minfit, maxfit);
@@ -44,6 +56,14 @@ vector<vector<double>> FitKnownPeaks_RtrnArray(TH1F* hist){
allPeaks[i]->SetLineStyle(7);
allPeaks[i]->SetNpx(500);
}
+ for(int i=numGauss; i<numPeaks; i++){
+ string nameHere = nameBase;
+ nameHere += to_string(i);
+ allPeaks[i] = new TF1(nameHere.c_str(), "[3] * TMath::Voigt(x-[0],[1],[2])", minfit, maxfit);
+ allPeaks[i]->SetLineColor(kCyan+2);
+ allPeaks[i]->SetLineStyle(7);
+ allPeaks[i]->SetNpx(500);
+ }
// Build background function (SIMPLE FOR NOW!)
TF1 *bg = new TF1("bg","[0]",minfit, maxfit);
@@ -51,40 +71,79 @@ vector<vector<double>> FitKnownPeaks_RtrnArray(TH1F* hist){
bg->SetLineStyle(9);
bg->SetParNames("Background");
- // Build total function
- int bgRefNum = ((numPeaks-1)*3)+3;
+//// Build total function: Original
+//int bgRefNum = ((numGauss-1)*3)+3;
+//string s_full = "gaus(0)";
+//for(int i=1; i<numGauss; 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);
+
+ // Build total function: with Voigt functions
+ int bgRefNum = ((numGauss-1)*3)+3;
string s_full = "gaus(0)";
- for(int i=1; i<numPeaks; i++){
+ for(int i=1; i<numGauss; i++){
s_full += "+ gaus("+ to_string(3*i) + ")";
}
+ for(int i=0; i<numVoigt; i++){
+ s_full += "+ [" + to_string(3*numGauss + i + 3)
+ + "] * TMath::Voigt(x-["
+ + to_string(3*numGauss + i + 0)
+ + "],["
+ + to_string(3*numGauss + i + 1)
+ + "],["
+ + to_string(3*numGauss + i + 2)
+ + "])";
+ }
+ cout << s_full << endl;
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);
+ for(int i=0; i<numGauss; i++){
+ full->FixParameter((i*3)+1, means.at(i));
+ full->FixParameter((i*3)+2, GetSigmaFromFunction(means.at(i)));
full->SetParameter((i*3)+0, 1e1);
full->SetParLimits((i*3)+0, 0.0,1e8);
}
+ for(int i=0; i<numVoigt; i++){
+ full->FixParameter((3*numGauss)+i+0, meansVoigt.at(i));
+ full->FixParameter((3*numGauss)+i+1, GetSigmaFromFunction(meansVoigt.at(i)));
+ full->SetParameter((3*numGauss)+i+2, 0.5);
+ full->SetParLimits((3*numGauss)+i+2, 0.0,1.0);
+ full->SetParameter((3*numGauss)+i+3, 1e2);
+ full->SetParLimits((3*numGauss)+i+3, 0.0,1e8);
+ cout << "fixing " << (3*numGauss)+i+ 0 << " to " << meansVoigt.at(i) << " and " << (3*numGauss)+i+1 << " to " << GetSigmaFromFunction(meansVoigt.at(i)) << endl;
+ }
+
+cout << " BEFORE FIT" << endl;
// Fit full function
hist->Fit(full, "RWQB", "", minfit, maxfit);
+cout << " AFTER FIT" << endl;
// Extract parameters
const Double_t* finalPar = full->GetParameters();
const Double_t* finalErr = full->GetParErrors();
- for(int i=0; i<numPeaks; i++){
+ for(int i=0; i<numGauss; i++){
allPeaks[i]->SetParameters(finalPar[0+(i*3)],
finalPar[1+(i*3)],
finalPar[2+(i*3)]);
allPeaks[i]->SetLineColor(kOrange);
}
-
+cout << " MID EXTRACTION" << endl;
+ for(int i=0; i<numVoigt; i++){
+ allPeaks[numGauss+i]->SetParameters(finalPar[0+i+(numGauss*3)],
+ finalPar[1+i+(numGauss*3)],
+ finalPar[2+i+(numGauss*3)],
+ finalPar[3+i+(numGauss*3)]);
+// allPeaks[i]->SetLineColor(kOrange);
+ }
+cout << " AFTER EXTRACTION" << endl;
+
//Draw all peaks
hist->Draw();
for(int i=0; i<numPeaks; i++){
@@ -97,7 +156,7 @@ vector<vector<double>> FitKnownPeaks_RtrnArray(TH1F* hist){
//Loop over every peak
vector<vector<double>> allpeaks;
- for(int i=0; i<numPeaks; i++){
+ for(int i=0; i<numGauss; i++){
//For AREA = HEIGHT * SIGMA * SQRT(2*PI)
double A = (finalPar[0+(i*3)] * finalPar[2+(i*3)] * sqrt(2*pi)) / binwidth;
@@ -123,6 +182,34 @@ vector<vector<double>> FitKnownPeaks_RtrnArray(TH1F* hist){
allpeaks.push_back(onepeak);
cout << "-------------" << endl;
}
+ for(int i=0; i<numVoigt; i++){
+ //For AREA = ?????????????????
+ //double A =
+
+ cout << "voigt funciton is area of unity" << endl;
+ cout << "therefore area of peaks is multiplier/binwidth = " << finalPar[(numGauss*3)+i+3]/binwidth << endl;
+
+ cout << fixed << setprecision(3)
+ << " #" << numGauss+i << " "
+ << finalPar[(numGauss*3)+i+0] << "\t" << setprecision(0)
+ << finalPar[(numGauss*3)+i+3] / binwidth << "\t+- "
+ << finalErr[(numGauss*3)+i+3] / binwidth << setprecision(3);
+ cout << " SQRT: " << sqrt(finalPar[(numGauss*3)+i+3]) << endl;
+
+ vector<double> onepeak; //energy, area and error for one peak
+ onepeak.push_back(finalPar[(numGauss*3)+i+0]);
+ onepeak.push_back(finalPar[(numGauss*3)+i+3]/binwidth);
+ onepeak.push_back(finalErr[(numGauss*3)+i+3]/binwidth);
+ allpeaks.push_back(onepeak);
+
+ cout << " VOIGT PEAK #" << i << " FUNCTION VALUES: [0] =" << finalPar[(numGauss*3)+i+0] << " , [1] = "
+ << finalPar[(numGauss*3)+i+1] << " , [2] = "
+ << finalPar[(numGauss*3)+i+2] << " , [3] = "
+ << finalPar[(numGauss*3)+i+3] << endl;
+ cout << "-------------" << endl;
+ }
+
+ cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! " << allpeaks.size() << endl;
//if(removing){
// for(int b=1; b<hist->GetNbinsX()-1; b++){
diff --git a/Projects/e775s/analysis/Plots_ParryFunctions.h b/Projects/e775s/analysis/Plots_ParryFunctions.h
index bbcbde7cedbf239f028a7e1bfd463c58e4bd4323..82fa69ea5efaa69b3d0cccb4d6934525667ed317 100644
--- a/Projects/e775s/analysis/Plots_ParryFunctions.h
+++ b/Projects/e775s/analysis/Plots_ParryFunctions.h
@@ -1,20 +1,24 @@
+#include <ctime>
vector<TH1F*> fitEx;
-vector<int> peaksToFit{7622, 7855, 8313, 8554, 8962, 9200, 9400};
+ vector<int> peaksToFit{7252, 7622, 7855, 8313, 8561, 8962, 9200, 9400, 9770};
+//vector<int> peaksToFit{7252, 7622, 7855, 8561, 8962, 9200, 9400, 9770};
+//vector<int> peaksToFit{7622, 7855, 8313, 8554, 8962, 9200, 9400};
+//vector<int> peaksToFit{7622, 7855, 8313, 8554, 8962, 9200, 9770};
//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";
+string basePath = "~/Programs/nptoolv3/Outputs/Analysis/Sim_19Odp_";
+string gateBase = "Mugast.TelescopeNumber<9";
+string gateG = gateBase + " && abs(T_FPMW_HF-13300)<1500";
+string gateB = gateBase + " && abs(T_FPMW_HF-10500)<1500";
-//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 angleMin = 104; int angleMax = 134; double angleBinWidth = 5; double exMin = 6; double exMax = 11; double exBinWidth = 0.05;
+//int angleMin = 104; int angleMax = 134; double angleBinWidth = 30; double exMin = 6; double exMax = 11; double exBinWidth = 0.05;
+//int angleMin = 104; int angleMax = 134; double angleBinWidth = 10; 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) + ","
@@ -113,34 +117,40 @@ vector<TH2F*> Parry_DrawPlots_SimPeaks(vector<int> peaks){
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));
+ TFractionFitter* fracfit = new TFractionFitter(exp, sim,"V");
+ //fracfit->SetRangeX(exp->FindBin(7.4), exp->FindBin(10.5));
+ fracfit->SetRangeX(exp->FindBin(7.0), exp->FindBin(10.5));
+ //fracfit->SetRangeX(exp->FindBin(7.0), exp->FindBin(8.7));
for(int i=0; i<sim->GetEntries(); i++){ fracfit->Constrain(i,0.0,0.4); }
+ //fracfit->ErrorAnalysis(2.0);
+
+////TESTING WAYS TO CHANGE FIT OPTIONS TO FIX FIT FAILURE
+//ROOT::Fit::Fitter* fitoptions = (ROOT::Fit::Fitter*)fracfit->GetFitter();
+//ROOT::Math::MinimizerOptions opt = fitoptions->Config().MinimizerOptions();
+//opt.SetDefaultPrecision(1e-2);
+//opt.SetDefaultTolerance(1e-2);
+//fitoptions->EvalFCN();
+//fitoptions->Config().SetMinimizerOptions(opt);
+
+
//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; }
+
+ //CHECK IF FIT COMPLETED CORRECTLY
+ if(fitstatus!=0){ for(int n=0; n<10; n++){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]); }
+ 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] << " || "; }
@@ -149,6 +159,23 @@ vector<pair<double,double>> Parry_FitSingleAngle(TH1F* exp, TObjArray* sim, int
<< fracfit->GetNDF() << " = "
<< fracfit->GetChisquare() / fracfit->GetNDF()
<< RESET << endl;
+
+ cout << "FIT PROBABILITY = " << fracfit->GetProb() << endl;
+
+ ofstream chi2file;
+ chi2file.open("ParryFit_Chi2NDF.txt",std::ios_base::app);
+ time_t timestamp; time(×tamp);
+ chi2file << ctime(×tamp) << endl;
+ chi2file << " Peaks = ";
+ for(int i=0; i<peaksToFit.size(); i++){ chi2file << peaksToFit[i] << " ";}
+ chi2file << endl;
+ chi2file << " Chi2 / NDF = " << fracfit->GetChisquare() << " / "
+ << fracfit->GetNDF() << " = "
+ << fracfit->GetChisquare() / fracfit->GetNDF()
+ << endl;
+
+ chi2file.close();
+
//SEE OUTPUTS
string canvname = "canv" + to_string(angle);
@@ -181,7 +208,6 @@ vector<pair<double,double>> Parry_FitSingleAngle(TH1F* exp, TObjArray* sim, int
AllPeaksOneAngle.push_back(PointErr);
}
-//return fracfit->GetChisquare() / fracfit->GetNDF();
return AllPeaksOneAngle;
}
@@ -204,6 +230,8 @@ void Parry_RunWholeRange(){
for(int a=0; a<angleBinNumb; a++){
string n = "hAng" + to_string(a);
TH1F* hExp1D = (TH1F*) hExp2D->ProjectionY(n.c_str(), a, a+1);
+
+ //Flatten over-subtracted bins
for(int bin=0; bin<hExp1D->GetNbinsX(); bin++){
if(hExp1D->GetBinContent(bin)<0){
hExp1D->SetBinContent(bin,0);
@@ -211,7 +239,7 @@ void Parry_RunWholeRange(){
}
}
hExp1D->SetLineColor(kBlack); hExp1D->SetLineWidth(2); hExp1D->SetLineStyle(1);
- hExp1D->GetYaxis()->SetRangeUser(0,500);
+ hExp1D->GetYaxis()->SetRangeUser(0,400);
TObjArray *fitSim = new TObjArray(peaksToFit.size());
for(int p=0; p<peaksToFit.size(); p++){
string n2 = n + "_Peak" + to_string(p);
@@ -280,6 +308,3 @@ void Parry_RunWholeRange(){
-
-
-
diff --git a/Projects/e775s/analysis/Plots_TWOFNR.h b/Projects/e775s/analysis/Plots_TWOFNR.h
index 968b70d3b15166c4ad5f6185ec5fd72ce9820bf2..ef8e7d2b6222a3ebc9e433e392a8ea2126aa902f 100644
--- a/Projects/e775s/analysis/Plots_TWOFNR.h
+++ b/Projects/e775s/analysis/Plots_TWOFNR.h
@@ -10,9 +10,9 @@ void MoveDirectoryTo(string newDir){
}
-TGraph* TWOFNR(double Ex, double Jf, double particleN, double particleL, double particleJ){
+TGraph* TWOFNR(double Ex, double Jf, double particleN, double particleL, double particleJ, double Ji){
- /* This funciton moves between directories in roder to execute TWOFNR */
+ /* This funciton moves between directories in order to execute TWOFNR */
/* and as such CANNOT be run on another machine! */
@@ -22,14 +22,19 @@ TGraph* TWOFNR(double Ex, double Jf, double particleN, double particleL, double
cout << "========================================================" << endl;
//double Ex = 0.0;
- double Ji = 2.5;
+ //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
+//int model_out_A = 1; //BechettiGreenlees
+//int model_out_A = 2; //ChapelHill
+//int model_out_A = 3; //Menet
+ int model_out_A = 4; //Perey
+//int model_out_A = 5; //JLM
+//int model_out_A = 6; //KoningDelaroche
string njjj;
@@ -161,10 +166,10 @@ 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);
+ TGraph* s12 = TWOFNR(Ex, 2., 1., 0., 0.5, 2.5);
+ TGraph* p12 = TWOFNR(Ex, 2., 0., 1., 0.5, 2.5);
+ TGraph* d52 = TWOFNR(Ex, 0., 0., 2., 2.5, 2.5);
+ TGraph* f72 = TWOFNR(Ex, 3., 0., 3., 3.5, 2.5);
s12->SetLineColor(kRed); s12->SetLineWidth(2); s12->SetLineStyle(1);
p12->SetLineColor(kGreen); p12->SetLineWidth(1); p12->SetLineStyle(7);
diff --git a/Projects/e775s/analysis/run_autopsp.sh b/Projects/e775s/analysis/run_autopsp.sh
index fd136ee17359fa621cdb9d6f0841411aed08b915..132b7c34dc1ac64464ab2ef22cdb64c954177658 100755
--- a/Projects/e775s/analysis/run_autopsp.sh
+++ b/Projects/e775s/analysis/run_autopsp.sh
@@ -1,3 +1,4 @@
#!/bin/bash
-./run_psp.sh PhaseSpace_InitialTest 0000
-./run_psp.sh PhaseSpace_InitialTest 0096
+#./run_psp.sh PhaseSpace_InitialTest 0000
+#./run_psp.sh PhaseSpace_InitialTest 0096
+./run_psp.sh PhaseSpace_InitialTest 1471
diff --git a/Projects/e775s/analysis/run_autosim.sh b/Projects/e775s/analysis/run_autosim.sh
index 9bc95b15a9683f54e3fe16e574a38911c8dc2436..3f9deef4ad88389c48ccf6345653b7a489936bdd 100755
--- a/Projects/e775s/analysis/run_autosim.sh
+++ b/Projects/e775s/analysis/run_autosim.sh
@@ -4,11 +4,11 @@
#./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 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
@@ -23,8 +23,33 @@
#./run_sim.sh 19Odp 1675
#./run_sim.sh 19Odp 3572
#./run_sim.sh 19Odp 4071
+#./run_sim.sh 19Odp 4456 # Literature
#./run_sim.sh 19Odp 5004
#./run_sim.sh 19Odp 5228
#./run_sim.sh 19Odp 5629
+#./run_sim.sh 19Odp 7252 # Literature
#./run_sim.sh 19Odp 7610
#./run_sim.sh 19Odp 7865
+
+
+#For plotting sigma with energy
+./run_sim.sh 19Odp 0000
+./run_sim.sh 19Odp 0500
+./run_sim.sh 19Odp 1000
+./run_sim.sh 19Odp 1500
+./run_sim.sh 19Odp 2000
+./run_sim.sh 19Odp 2500
+./run_sim.sh 19Odp 3000
+./run_sim.sh 19Odp 3500
+./run_sim.sh 19Odp 4000
+./run_sim.sh 19Odp 4500
+./run_sim.sh 19Odp 5000
+./run_sim.sh 19Odp 5500
+./run_sim.sh 19Odp 6000
+./run_sim.sh 19Odp 6500
+./run_sim.sh 19Odp 7000
+./run_sim.sh 19Odp 7500
+./run_sim.sh 19Odp 8000
+./run_sim.sh 19Odp 8500
+./run_sim.sh 19Odp 9000
+./run_sim.sh 19Odp 9500
diff --git a/Projects/e775s/analysis/run_exp.sh b/Projects/e775s/analysis/run_exp.sh
index 03fbddf8d7dfc0cade9080c3005384571b99a5bb..76324852cf7d2ad715a4ba92d1a30a8c91b07fd8 100755
--- a/Projects/e775s/analysis/run_exp.sh
+++ b/Projects/e775s/analysis/run_exp.sh
@@ -6,10 +6,12 @@ make -j6;
#=====================
#====================================================
#rfile='./InputFiles/e775s.reaction'
- rfile='./InputFiles/e775s_2024-07-30.reaction'
+ #rfile='./InputFiles/e775s_2024-07-30.reaction'
+ rfile='./InputFiles/e775s_2024-09-06.reaction'
#----------------------------------------------------
#dfile='./InputFiles/mugast.detector'
- dfile='./InputFiles/mugast_2024-07-30.detector'
+ #dfile='./InputFiles/mugast_2024-07-30.detector'
+ dfile='./InputFiles/mugast_2024-09-06.detector'
#----------------------------------------------------
cfile='./InputFiles/Calibration.txt'
#====================================================
diff --git a/Projects/e775s/analysis/run_psp.sh b/Projects/e775s/analysis/run_psp.sh
index 764826417408ca72ea6d3c2267a24a2d1ca6ba19..29574797826d50332174e2ebd36909cf45eae0d2 100755
--- a/Projects/e775s/analysis/run_psp.sh
+++ b/Projects/e775s/analysis/run_psp.sh
@@ -5,14 +5,14 @@ rm ./InputFiles/Auto.runtotreat
rm ./InputFiles/AutoPSpace.reaction
sp=" "
-dfile="./InputFiles/mugast_2024-07-30.detector"
-rfile_exp="./InputFiles/e775s_2024-07-30.reaction"
+dfile="./InputFiles/mugast_2024-09-06.detector"
+rfile_exp="./InputFiles/e775s_2024-09-06.reaction"
rfile_sim="./InputFiles/AutoPSpace.reaction"
#Variables
E=$(echo "scale=3; $2/1000" | bc -l)
-X="-3.41"
-Y="-0.20"
+X="-3.8"
+Y="-0.4"
#====================================================
# Symbolic link for mugast map
@@ -65,10 +65,11 @@ 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
+for x in 0 1 2 3 4 5 6 7 8 9
do
outname=$nameconstruct"-"$x
-# npsimulation -D $dfile -E $rfile_sim -B ./InputFiles/runPhaseSpace.mac --random-seed $RANDOM -O $outname;
+ npsimulation -D $dfile -E $rfile_sim -B ./InputFiles/runPhaseSpace.mac --random-seed $RANDOM -O $outname;
+ #npsimulation -D $dfile -E $rfile_sim -B ./InputFiles/runSimSmall.mac --random-seed $RANDOM -O $outname;
done
outfile="PSp_"$nameconstruct
@@ -79,7 +80,7 @@ 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
+for x in 0 1 2 3 4 5 6 7 8 9
do
echo "$sp""$directory""$nameconstruct""-""$x"".root" >> ./InputFiles/Auto.runtotreat
done
diff --git a/Projects/e775s/analysis/run_sim.sh b/Projects/e775s/analysis/run_sim.sh
index 619e9cc6964318ddf85f30067b2ee92b816a5f8f..0e3f21980bed4f63dbd439dc1dd26e8358a8aa68 100755
--- a/Projects/e775s/analysis/run_sim.sh
+++ b/Projects/e775s/analysis/run_sim.sh
@@ -5,13 +5,14 @@ rm ./InputFiles/Auto.runtotreat
rm ./InputFiles/Auto.reaction
sp=" "
-dfile="./InputFiles/mugast_2024-07-30.detector"
+dfile="./InputFiles/mugast_2024-09-06.detector"
rfile="./InputFiles/Auto.reaction"
+#rfile="./InputFiles/SpecificSim.reaction"
#Variables
E=$(echo "scale=3; $2/1000" | bc -l)
-X="-3.41"
-Y="-0.20"
+X="-3.8"
+Y="-0.4"
#====================================================
# Symbolic link for mugast map
@@ -69,11 +70,11 @@ 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
+for x in 0 #1 2 3 4 5 6 7 8 9
do
outname=$nameconstruct"-"$x
- npsimulation -D $dfile -E $rfile -B ./InputFiles/runSimulation.mac --random-seed $RANDOM -O $outname;
+ #npsimulation -D $dfile -E $rfile -B ./InputFiles/runSimulation.mac --random-seed $RANDOM -O $outname;
+ npsimulation -D $dfile -E $rfile -B ./InputFiles/runSimSmall.mac --random-seed $RANDOM -O $outname;
done
outfile="Sim_"$nameconstruct
@@ -87,7 +88,7 @@ defineState='S'$2
echo "Defining state = "$defineState
-for x in 1 2 3 4 5 6 7 8 9 10
+for x in 0 #1 2 3 4 5 6 7 8 9
do
echo "$sp""$directory""$nameconstruct""-""$x"".root" >> ./InputFiles/Auto.runtotreat
done