diff --git a/Projects/e793s/Analysis.cxx b/Projects/e793s/Analysis.cxx old mode 100755 new mode 100644 index 23d8d12000e961c2a14fd4ed2e95242c5ed32a2e..10aee429561c4b3161b00fea9f20b25c17480c76 --- a/Projects/e793s/Analysis.cxx +++ b/Projects/e793s/Analysis.cxx @@ -62,10 +62,10 @@ void Analysis::Init() { } agata_zShift=51*mm; - //BrhoRef=0.65; + + //For calculating solid angle when processing simulated data if(isSim && !isPhaseSpace){ -//cout << "here_InIsSimLoop" << endl; Initial = new TInitialConditions(); ReactionConditions = new TReactionConditions(); RootInput::getInstance()->GetChain()->SetBranchAddress("InitialConditions",&Initial); @@ -73,24 +73,48 @@ void Analysis::Init() { // Initilize MGX histograms, because must be inside function string base1 = "ThetaCM_detected_MG"; string base2 = "ThetaLab_detected_MG"; + string base1v = "ThetaCM_detected_MG_wVAMOS"; + string base2v = "ThetaLab_detected_MG_wVAMOS"; + string base1unb = "ThetaCM_detected_MG_Unb"; + string base2unb = "ThetaLab_detected_MG_Unb"; for(int i=0; i<6; i++){ int j=i+1; if (i==6){int j=7;} string name1 = base1 + to_string(j); string name2 = base2 + to_string(j); + string name1v = base1v + to_string(j); + string name2v = base2v + 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); //900 bins for 0.2 angular bin width ThetaLab_detected_MGX[i] = new TH1F(name2.c_str(),name2.c_str(),NumThetaAngleBins,0,180); + ThetaCM_detected_MGX_wV[i] = new TH1F(name1v.c_str(),name1v.c_str(),NumThetaAngleBins,0,180); //900 bins for 0.2 angular bin width + ThetaLab_detected_MGX_wV[i] = new TH1F(name2v.c_str(),name2v.c_str(),NumThetaAngleBins,0,180); + ThetaCM_detected_MGX_Unb[i] = new TH1F(name1unb.c_str(),name1unb.c_str(),NumThetaAngleBins,0,180); //900 bins for 0.2 angular bin width + ThetaLab_detected_MGX_Unb[i] = new TH1F(name2unb.c_str(),name2unb.c_str(),NumThetaAngleBins,0,180); } // Initilize MMX histograms string base3 = "ThetaCM_detected_MM"; string base4 = "ThetaLab_detected_MM"; + string base3v = "ThetaCM_detected_MM_wVAMOS"; + string base4v = "ThetaLab_detected_MM_wVAMOS"; + string base3unb = "ThetaCM_detected_MM_Unb"; + string base4unb = "ThetaLab_detected_MM_Unb"; for(int i=0; i<5; i++){ int j=i+1; string name1 = base3 + to_string(j); string name2 = base4 + to_string(j); + string name1v = base3v + to_string(j); + string name2v = base4v + to_string(j); + string name1unb = base3unb + to_string(j); + string name2unb = base4unb + to_string(j); ThetaCM_detected_MMX[i] = new TH1F(name1.c_str(),name1.c_str(),NumThetaAngleBins,0,180); ThetaLab_detected_MMX[i] = new TH1F(name2.c_str(),name2.c_str(),NumThetaAngleBins,0,180); + ThetaCM_detected_MMX_wV[i] = new TH1F(name1v.c_str(),name1v.c_str(),NumThetaAngleBins,0,180); + ThetaLab_detected_MMX_wV[i] = new TH1F(name2v.c_str(),name2v.c_str(),NumThetaAngleBins,0,180); + ThetaCM_detected_MMX_Unb[i] = new TH1F(name1unb.c_str(),name1unb.c_str(),NumThetaAngleBins,0,180); + ThetaLab_detected_MMX_Unb[i] = new TH1F(name2unb.c_str(),name2unb.c_str(),NumThetaAngleBins,0,180); } } @@ -101,7 +125,6 @@ void Analysis::Init() { M2 = (TMust2Physics*) m_DetectorManager -> GetDetector("M2Telescope"); 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()); @@ -133,7 +156,7 @@ void Analysis::Init() { LightAl = NPL::EnergyLoss(light+"_Al.G4table" ,"G4Table",100); LightSi = NPL::EnergyLoss(light+"_Si.G4table" ,"G4Table",100); BeamTargetELoss = NPL::EnergyLoss(beam+"_"+TargetMaterial+".G4table","G4Table",100); - //HeavyTargetELoss = NPL::EnergyLoss(heavy+"_"+TargetMaterial+".G4table","G4Table",100); +//HeavyTargetELoss = NPL::EnergyLoss(heavy+"_"+TargetMaterial+".G4table","G4Table",100); FinalBeamEnergy = BeamTargetELoss.Slow(OriginalBeamEnergy, 0.5*TargetThickness, 0); reaction.SetBeamEnergy(FinalBeamEnergy); @@ -167,28 +190,19 @@ void Analysis::Init() { Z.clear(); dE = 0; BeamDirection = TVector3(0,0,1); - //nbTrack=0; - //nbHits=0; - //count=0; AHeavy=reaction.GetParticle4()->GetA(); ALight=reaction.GetParticle3()->GetA(); MHeavy=reaction.GetParticle4()->Mass(); MLight=reaction.GetParticle3()->Mass(); - //bool writetoscreen=true; for(int i=0;i<GATCONF_SIZE;i++){ // loop over the bits GATCONF_Counter[i] = 0 ; } - -// ThetaCM_detected->Sumw2(); -// ThetaLab_detected->Sumw2(); -//cout << "here_endInit" << endl; } //////////////////////////////////////////////////////////////////////////////// void Analysis::TreatEvent(){ // Reinitiate calculated variable -//cout << "here_TreatEvent" << endl; ReInitValue(); if(isSim && !isPhaseSpace){ ThetaCM_emmitted->Fill(ReactionConditions->GetThetaCM()); @@ -199,7 +213,6 @@ void Analysis::TreatEvent(){ for(int i=0;i<GATCONF_SIZE;i++){ // loop over the bits if(GATCONF_MASTER & (unsigned int)pow(2,i)){ // test if ith bit is on GATCONF_Counter[i]++ ; // increment the array - // cout << std::bitset<16>(GATCONF_MASTER) << " " << i+1 << endl; } } @@ -212,47 +225,17 @@ void Analysis::TreatEvent(){ TVector3 BeamDirection(0.,0.,1.); BeamImpact = TVector3(XBeam,YBeam,m_DetectorManager->GetTargetZ()); - //ParticleMult=M2->Si_E.size();////+MG->DSSD_E.size(); ParticleMult=M2->Si_E.size()+MG->DSSD_E.size(); -//cout << "here_BeforeMustLoop" << endl; //////////////////////////////////////////////////////////////////////////// //////////////////////////////// LOOP on MUST2 //////////////////////////// //////////////////////////////////////////////////////////////////////////// unsigned int sizeM2 = M2->Si_E.size(); for(unsigned int countMust2 = 0; countMust2 < sizeM2; countMust2++){ -//cout << "here_InMustLoop" << endl; - /************************************************/ // Part 0 : Get the useful Data // MUST2 int TelescopeNumber = M2->TelescopeNumber[countMust2]; - /**/ - if(isSim && !isPhaseSpace){ - if(M2->TelescopeNumber[countMust2]<5){ - if(M2->Si_E[countMust2]>0 && // DSSD count - M2->CsI_E[countMust2]<=0 && // No CsI count - M2->Si_T[countMust2]<460 // Triton kinematic line, not punch through - ){ - ThetaCM_detected_MM->Fill(ReactionConditions->GetThetaCM()); - ThetaLab_detected_MM->Fill(ReactionConditions->GetTheta(0)); - - int MMX = TelescopeNumber-1; - ThetaCM_detected_MMX[MMX]->Fill(ReactionConditions->GetThetaCM()); - ThetaLab_detected_MMX[MMX]->Fill(ReactionConditions->GetTheta(0)); - } - } else { - //No triton requirement for MM5 - ThetaCM_detected_MM->Fill(ReactionConditions->GetThetaCM()); - ThetaLab_detected_MM->Fill(ReactionConditions->GetTheta(0)); - - int MMX = TelescopeNumber-1; - ThetaCM_detected_MMX[MMX]->Fill(ReactionConditions->GetThetaCM()); - ThetaLab_detected_MMX[MMX]->Fill(ReactionConditions->GetTheta(0)); - } - } - /**/ - /************************************************/ // Part 1 : Impact Angle ThetaM2Surface = 0; @@ -279,78 +262,123 @@ void Analysis::TreatEvent(){ // if CsI if(CsI_E_M2>0 ){ - // The energy in CsI is calculate form dE/dx Table because Energy = CsI_E_M2; - //if(!isSim){ Energy = LightAl.EvaluateInitialEnergy(Energy, 0.4*micrometer, ThetaM2Surface); - //} Energy+=Si_E_M2; - } - else + }else{ Energy = Si_E_M2; - + } + RawEnergy.push_back(Energy); -// if(!isSim){ - // Evaluate energy using the thickness - elab_tmp = LightAl.EvaluateInitialEnergy( - Energy, - 0.4*micrometer, - ThetaM2Surface); - ELoss_Al.push_back(Energy-elab_tmp); - double elab_tmp2 = elab_tmp; - // Target Correction - elab_tmp = LightTarget.EvaluateInitialEnergy( - elab_tmp, - 0.5*TargetThickness, - ThetaNormalTarget); - ELoss_Target.push_back(elab_tmp2-elab_tmp); - ELoss.push_back(Energy-elab_tmp); -// } else {elab_tmp = Energy;} + // Evaluate energy using the thickness + elab_tmp = LightAl.EvaluateInitialEnergy( + Energy, + 0.4*micrometer, + ThetaM2Surface); + ELoss_Al.push_back(Energy-elab_tmp); + double elab_tmp2 = elab_tmp; + + // Target Correction + elab_tmp = LightTarget.EvaluateInitialEnergy( + elab_tmp, + 0.5*TargetThickness, + ThetaNormalTarget); + ELoss_Target.push_back(elab_tmp2-elab_tmp); + ELoss.push_back(Energy-elab_tmp); ELab.push_back(elab_tmp); - /************************************************/ /************************************************/ // Part 3 : Excitation Energy Calculation - //Ex.push_back(reaction.ReconstructRelativistic(elab_tmp,thetalab_tmp)); + //Ex.push_back(reaction.ReconstructRelativistic(elab_tmp,thetalab_tmp)); Ex.push_back(reaction.ReconstructRelativistic(elab_tmp,thetalab_tmp, philab_tmp)); Ecm.push_back(Energy*(AHeavy+ALight)/(4*AHeavy*cos(thetalab_tmp)*cos(thetalab_tmp))); /************************************************/ /************************************************/ // Part 4 : Theta CM Calculation - ThetaCM.push_back(reaction.EnergyLabToThetaCM(elab_tmp, thetalab_tmp)/deg); /************************************************/ ThetaLab.push_back(thetalab_tmp/deg); PhiLab.push_back(philab_tmp/deg); - //if(isSim && !isPhaseSpace){ - // ThetaCM_detected_MM->Fill(reaction.EnergyLabToThetaCM(elab_tmp, thetalab_tmp)/deg); - // ThetaLab_detected_MM->Fill(thetalab_tmp/deg); - - // int MMX = TelescopeNumber-1; - // ThetaCM_detected_MMX[MMX]->Fill(reaction.EnergyLabToThetaCM(elab_tmp, thetalab_tmp)/deg); - // ThetaLab_detected_MMX[MMX]->Fill(thetalab_tmp/deg); - //} - + if(isSim && !isPhaseSpace){ + if(M2->TelescopeNumber[countMust2]<5){ + if(M2->Si_E[countMust2]>0 && // DSSD count + M2->CsI_E[countMust2]<=0 && // No CsI count + M2->Si_T[countMust2]<460 // Triton kinematic line, not punch through + ){ + ThetaCM_detected_MM->Fill(reaction.EnergyLabToThetaCM(elab_tmp, thetalab_tmp) + /deg); + ThetaLab_detected_MM->Fill(thetalab_tmp/deg); + + int MMX = TelescopeNumber-1; + ThetaCM_detected_MMX[MMX]->Fill(reaction.EnergyLabToThetaCM(elab_tmp, + thetalab_tmp)/deg); + ThetaLab_detected_MMX[MMX]->Fill(thetalab_tmp/deg); + } + } else { + //No triton requirement for MM5 + ThetaCM_detected_MM->Fill(ReactionConditions->GetThetaCM()); + ThetaLab_detected_MM->Fill(ReactionConditions->GetTheta(0)); + + int MMX = TelescopeNumber-1; + ThetaCM_detected_MMX[MMX]->Fill(ReactionConditions->GetThetaCM()); + ThetaLab_detected_MMX[MMX]->Fill(ReactionConditions->GetTheta(0)); + } + if(M2->TelescopeNumber[countMust2]<5){ + if(M2->Si_E[countMust2]>0 && // DSSD count + M2->CsI_E[countMust2]<=0 && // No CsI count + M2->Si_T[countMust2]<460 // Triton kinematic line, not punch through + ){ + ThetaCM_detected_MM_wV->Fill(reaction.EnergyLabToThetaCM(elab_tmp, thetalab_tmp) + /deg); + ThetaLab_detected_MM_wV->Fill(thetalab_tmp/deg); + + int MMX = TelescopeNumber-1; + ThetaCM_detected_MMX_wV[MMX]->Fill(reaction.EnergyLabToThetaCM(elab_tmp, + thetalab_tmp)/deg); + ThetaLab_detected_MMX_wV[MMX]->Fill(thetalab_tmp/deg); + } + } else { + //No triton requirement for MM5 + ThetaCM_detected_MM_wV->Fill(ReactionConditions->GetThetaCM()); + ThetaLab_detected_MM_wV->Fill(ReactionConditions->GetTheta(0)); + + int MMX = TelescopeNumber-1; + ThetaCM_detected_MMX_wV[MMX]->Fill(ReactionConditions->GetThetaCM()); + ThetaLab_detected_MMX_wV[MMX]->Fill(ReactionConditions->GetTheta(0)); + } + } + double HeavyThetaLab_tmp = reaction.GetEnergyImpulsionLab_4().Theta(); + double HeavyPhiLab_tmp = reaction.GetEnergyImpulsionLab_4().Phi(); + HeavyThetaLab.push_back(HeavyThetaLab_tmp/deg); + HeavyPhiLab.push_back(HeavyPhiLab_tmp/deg); + + auto tmpHeavy = new TVector3(1,1,1); + tmpHeavy->SetTheta(HeavyThetaLab_tmp); + tmpHeavy->SetPhi(HeavyPhiLab_tmp); + tmpHeavy->SetMag(WindowDistance/cos(HeavyThetaLab_tmp)); + 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 << "here_EndMustLoop" << endl; } -//cout << "here_BeforeMugastLoop" << endl; //////////////////////////////////////////////////////////////////////////// //////////////////////////////// LOOP on MUGAST //////////////////////////// //////////////////////////////////////////////////////////////////////////// unsigned int sizeMG = MG->DSSD_E.size(); for(unsigned int countMugast = 0; countMugast<sizeMG; countMugast++){ - if(isSim && !isPhaseSpace){ int MGX = MG->TelescopeNumber[0]; MGX = MGX-1; @@ -369,14 +397,6 @@ void Analysis::TreatEvent(){ thetalab_tmp = 0; philab_tmp = 0; TVector3 HitDirection = MG->GetPositionOfInteraction(countMugast) - BeamImpact; - //TVector3 tempVector; - //if(MG->TelescopeNumber[0]==3){ - // tempVector = {MG->GetPositionOfInteraction(countMugast).X(), - // MG->GetPositionOfInteraction(countMugast).Y(), - // MG->GetPositionOfInteraction(countMugast).Z()+1.0}; //add 1mm to MG3 - // HitDirection = tempVector - BeamImpact; - // if(warning){cout << "!!! EDITING MG3 !!!" << endl; warning=false;} - //} thetalab_tmp = HitDirection.Angle(BeamDirection); philab_tmp = HitDirection.Phi(); @@ -390,6 +410,7 @@ void Analysis::TreatEvent(){ // Part 2 : Impact Energy Energy = elab_tmp = 0; Energy = MG->GetEnergyDeposit(countMugast); + double rawe = Energy; //for parallel line unbound state calc later, new addiiton RawEnergy.push_back(Energy); if(!isSim){ @@ -407,14 +428,14 @@ void Analysis::TreatEvent(){ ELoss.push_back(Energy-elab_tmp); } else { //TESTING DIFFERENT ENERGY LOSSES IN SIMULATION elab_tmp = Energy; //so I can add and remove sections - //elab_tmp = LightSi.EvaluateInitialEnergy( - // elab_tmp, //particle energy after Si - // 0.5*500.*micrometer, //thickness of Si - // ThetaMGSurface); //angle of impingement - //elab_tmp = LightAl.EvaluateInitialEnergy( - // elab_tmp, //particle energy after Al - // 0.4*micrometer, //thickness of Al - // ThetaMGSurface); //angle of impingement + //elab_tmp = LightSi.EvaluateInitialEnergy( + // elab_tmp, //particle energy after Si + // 0.5*500.*micrometer, //thickness of Si + // ThetaMGSurface); //angle of impingement + //elab_tmp = LightAl.EvaluateInitialEnergy( + // elab_tmp, //particle energy after Al + // 0.4*micrometer, //thickness of Al + // ThetaMGSurface); //angle of impingement elab_tmp = LightTarget.EvaluateInitialEnergy( elab_tmp, //particle energy after leaving target TargetThickness*0.5, //distance passed through target @@ -423,20 +444,9 @@ void Analysis::TreatEvent(){ ELab.push_back(elab_tmp); - //cout << "===============" << endl; - //cout << "RawE:\t" << RawEnergy.back() << endl; - //cout << "ELAl:\t" << ELoss_Al.back() << endl; - //cout << "ELCD:\t" << ELoss_Target.back() << endl; - //cout << "ELTt:\t" << ELoss.back() << endl; - //cout << "ELab:\t" << ELab.back() << endl; - - // Part 3 : Excitation Energy Calculation - //if(!isSim){ //TESTING!!!! - //Ex.push_back(reaction.ReconstructRelativistic(elab_tmp,thetalab_tmp)); - Ex.push_back(reaction.ReconstructRelativistic(elab_tmp,thetalab_tmp, philab_tmp)); - Ecm.push_back(elab_tmp*(AHeavy+ALight)/(4*AHeavy*cos(thetalab_tmp)*cos(thetalab_tmp))); - //} + Ex.push_back(reaction.ReconstructRelativistic(elab_tmp,thetalab_tmp, philab_tmp)); + Ecm.push_back(elab_tmp*(AHeavy+ALight)/(4*AHeavy*cos(thetalab_tmp)*cos(thetalab_tmp))); // Part 4 : Theta CM Calculation ThetaLab.push_back(thetalab_tmp/deg); @@ -451,9 +461,35 @@ void Analysis::TreatEvent(){ MG_D = MG->TelescopeNumber[0]; } - }//end loop Mugast + if(isSim && !isPhaseSpace){ + if(reaction.ReconstructRelativistic(elab_tmp,thetalab_tmp, philab_tmp)<5.9 + && elab_tmp>1.5 + && thetalab_tmp/deg>104. + && elab_tmp<((0.07*(thetalab_tmp/deg))-5.39) + && elab_tmp>((0.07*(thetalab_tmp/deg))-6.18)) + { + ThetaCM_detected_MG_Unb->Fill(ReactionConditions->GetThetaCM()); + ThetaLab_detected_MG_Unb->Fill(ReactionConditions->GetTheta(0)); + } + } + double HeavyThetaLab_tmp = reaction.GetEnergyImpulsionLab_4().Theta(); + double HeavyPhiLab_tmp = reaction.GetEnergyImpulsionLab_4().Phi(); + HeavyThetaLab.push_back(HeavyThetaLab_tmp/deg); + HeavyPhiLab.push_back(HeavyPhiLab_tmp/deg); + + auto tmpHeavy = new TVector3(1,1,1); + tmpHeavy->SetTheta(HeavyThetaLab_tmp*(180./M_PI)); + tmpHeavy->SetPhi(HeavyPhiLab_tmp*(180./M_PI)); + tmpHeavy->SetMag(WindowDistance/cos(HeavyThetaLab_tmp));//cos=A/H, using average of MUST2 corners as 'Z' + HeavyPosXAtWindow.push_back(tmpHeavy->X()); + HeavyPosYAtWindow.push_back(tmpHeavy->Y()); + HeavyAngX.push_back(atan(tmpHeavy->X()/WindowDistance)); + HeavyAngY.push_back(atan(tmpHeavy->Y()/WindowDistance)); + HeavyPosXAtMUST2.push_back(MUST2InnerDistance*tan( atan(tmpHeavy->X()/WindowDistance) )); + HeavyPosYAtMUST2.push_back(MUST2InnerDistance*tan( atan(tmpHeavy->Y()/WindowDistance) )); + }//end loop Mugast //////////////////////////////////////////////////////////////////////////// @@ -502,20 +538,6 @@ void Analysis::TreatEvent(){ TVector3 GammaHit(AddX[i],AddY[i],AddZ[i]+agata_zShift); // TVector3 GammaHit(trackX1[0],trackY1[0],trackZ1[0]); - -// cout << "gDir " << GammaHit.X() << " " -// << GammaHit.Y() << " " -// << GammaHit.Z() << " "; - /* TEST adding constant to Phi */ - //double phitemp = GammaHit.Phi(); - //double phiconst = -2.0; - //GammaHit.SetPhi(phitemp+phiconst); -// GammaHit.RotateZ(0.0); - -// cout << "gDir " << GammaHit.X() << " " -// << GammaHit.Y() << " " -// << GammaHit.Z() << endl; - // Gamma Direction TVector3 GammaDirection = GammaHit-BeamImpact; GammaDirection = GammaDirection.Unit(); @@ -523,21 +545,14 @@ void Analysis::TreatEvent(){ /* Beta from Two body kinematic */ TVector3 beta = reaction.GetEnergyImpulsionLab_4().BoostVector(); - //TVector3 beta = reaction->GetEnergyImpulsionLab_4().BoostVector(); - -// cout << "bDir " << beta.X() << " " -// << beta.Y() << " " -// << beta.Z() ; -// beta.RotateX(0.002847); beta.RotateY(3.144869); beta.RotateZ(0.095923); beta.RotateY(M_PI); -// cout << "bDir " << beta.X() << " " -// << beta.Y() << " " -// << beta.Z() << endl; - + /* Original beta */ -// TVector3 beta(0,0,-0.1257); + //TVector3 beta(0,0,-0.1257); // For beta rotation minimization + AGATA_Theta.push_back(GammaDirection.Theta()); + AGATA_Phi.push_back(GammaDirection.Phi()); AGATA_GammaPx.push_back(Egamma*GammaDirection.X()); AGATA_GammaPy.push_back(Egamma*GammaDirection.Y()); AGATA_GammaPz.push_back(Egamma*GammaDirection.Z()); @@ -579,7 +594,6 @@ void Analysis::TreatEvent(){ MWT[i] = (MW_T[i]-offset[j])/slope[j]; } - } //////////////////////////////////////////////////////////////////////////////// @@ -612,17 +626,6 @@ void FillSolidAngles(TH1F* hSA, TH1F* hDet, TH1F* hEmm){ } } - -//void DivideByCline(TH1F* histo){ -// -// for(int b=0; b<NumThetaAngleBins; b++){ -// int store = histo->GetBinContent(b); -// histo->SetBinContent(b,(int) store/clineVal.at(b)); -/// } -//} - - - //////////////////////////////////////////////////////////////////////////////// void Analysis::End(){ cout << endl << "\e[1;32m GATCONF Statistics " << endl ; @@ -653,23 +656,56 @@ void Analysis::End(){ 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); - //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); - + 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); + FillSolidAngles(SolidAngle_CM_MM, ThetaCM_detected_MM, ThetaCM_emmitted); FillSolidAngles(SolidAngle_Lab_MM, ThetaLab_detected_MM, 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); - + 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); + HistList->Add(ThetaCM_emmitted); HistList->Add(ThetaLab_emmitted); HistList->Add(ThetaCM_detected_MM); @@ -678,9 +714,18 @@ void Analysis::End(){ HistList->Add(Efficiency_Lab_MM); HistList->Add(SolidAngle_CM_MM); HistList->Add(SolidAngle_Lab_MM); - //HistList->Add(Cline_MM); - //HistList->Add(Cline); - + 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); // MUGAST auto Efficiency_CM_MG = new TH1F(*ThetaCM_detected_MG); @@ -699,20 +744,54 @@ 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); - //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); + 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); + /* Testing method for better errors in SolidAngle histograms */ FillSolidAngles(SolidAngle_CM_MG, ThetaCM_detected_MG, ThetaCM_emmitted); - //SolidAngle_CM_MG->Divide(Cline_MG,1); - //SolidAngle_CM_MG->Divide(Cline); + FillSolidAngles(SolidAngle_CM_MG_wV, ThetaCM_detected_MG_wV, ThetaCM_emmitted); + FillSolidAngles(SolidAngle_CM_MG_Unb, ThetaCM_detected_MG_Unb, ThetaCM_emmitted); FillSolidAngles(SolidAngle_Lab_MG, ThetaLab_detected_MG, ThetaLab_emmitted); - //SolidAngle_Lab_MG->Divide(Cline_MG,1); - //SolidAngle_Lab_MG->Divide(Cline); + FillSolidAngles(SolidAngle_Lab_MG_wV, ThetaLab_detected_MG_wV, ThetaLab_emmitted); + FillSolidAngles(SolidAngle_Lab_MG_Unb, ThetaLab_detected_MG_Unb, ThetaLab_emmitted); HistList->Add(ThetaCM_detected_MG); HistList->Add(ThetaLab_detected_MG); @@ -720,7 +799,18 @@ void Analysis::End(){ HistList->Add(Efficiency_Lab_MG); HistList->Add(SolidAngle_CM_MG); HistList->Add(SolidAngle_Lab_MG); - //HistList->Add(Cline_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); // MUGAST INDIVIDUAL TH1F *SolidAngle_CM_MGX[6]; @@ -733,8 +823,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[6]; @@ -747,8 +835,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); } // MUST2 INDIVIDUAL @@ -761,8 +847,6 @@ void Analysis::End(){ SolidAngle_CM_MMX[i]->SetName(name.c_str()); SolidAngle_CM_MMX[i]->SetTitle(name.c_str()); FillSolidAngles(SolidAngle_CM_MMX[i], ThetaCM_detected_MMX[i], ThetaCM_emmitted); - //SolidAngle_CM_MMX[i]->Divide(Cline_MM,1); - //SolidAngle_CM_MMX[i]->Divide(Cline); } TH1F *SolidAngle_Lab_MMX[6]; @@ -774,8 +858,6 @@ void Analysis::End(){ SolidAngle_Lab_MMX[i]->SetName(name.c_str()); SolidAngle_Lab_MMX[i]->SetTitle(name.c_str()); FillSolidAngles(SolidAngle_Lab_MMX[i], ThetaLab_detected_MMX[i], ThetaLab_emmitted); - //SolidAngle_Lab_MMX[i]->Divide(Cline_MM,1); - //SolidAngle_Lab_MMX[i]->Divide(Cline); } for(int i=0; i<6; i++){HistList->Add(ThetaCM_detected_MGX[i]);} @@ -800,21 +882,19 @@ void Analysis::End(){ HistList->Write(); HistoFile->Close(); - cout << "!!! MAKE SURE YOU RENAME THE SOLID ANGLE FILE! !!!" << endl; } } //////////////////////////////////////////////////////////////////////////////// void Analysis::InitOutputBranch(){ -//cout << "here_InitOutput" << endl; - //RootOutput::getInstance()->GetTree()->Branch("Ex",&Ex,"Ex/D"); RootOutput::getInstance()->GetTree()->Branch("Ex",&Ex); - //RootOutput::getInstance()->GetTree()->Branch("EDC",&EDC,"EDC/D"); RootOutput::getInstance()->GetTree()->Branch("EDC",&EDC); RootOutput::getInstance()->GetTree()->Branch("AddBack_EDC",&AddBack_EDC); RootOutput::getInstance()->GetTree()->Branch("AddBack_EDC_Event1",&AddBack_EDC_Event1); RootOutput::getInstance()->GetTree()->Branch("AddBack_EDC_Event2",&AddBack_EDC_Event2); + RootOutput::getInstance()->GetTree()->Branch("AGATA_Theta",&AGATA_Theta); + RootOutput::getInstance()->GetTree()->Branch("AGATA_Phi",&AGATA_Phi); RootOutput::getInstance()->GetTree()->Branch("AGATA_GammaPx",&AGATA_GammaPx); RootOutput::getInstance()->GetTree()->Branch("AGATA_GammaPy",&AGATA_GammaPy); RootOutput::getInstance()->GetTree()->Branch("AGATA_GammaPz",&AGATA_GammaPz); @@ -832,6 +912,14 @@ void Analysis::InitOutputBranch(){ RootOutput::getInstance()->GetTree()->Branch("ELoss",&ELoss); RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab); RootOutput::getInstance()->GetTree()->Branch("PhiLab",&PhiLab); + RootOutput::getInstance()->GetTree()->Branch("HeavyThetaLab",&HeavyThetaLab); + RootOutput::getInstance()->GetTree()->Branch("HeavyPhiLab",&HeavyPhiLab); + RootOutput::getInstance()->GetTree()->Branch("HeavyPosXAtWindow",&HeavyPosXAtWindow); + RootOutput::getInstance()->GetTree()->Branch("HeavyPosYAtWindow",&HeavyPosYAtWindow); + RootOutput::getInstance()->GetTree()->Branch("HeavyPosXAtMUST2",&HeavyPosXAtMUST2); + RootOutput::getInstance()->GetTree()->Branch("HeavyPosYAtMUST2",&HeavyPosYAtMUST2); + RootOutput::getInstance()->GetTree()->Branch("HeavyAngX",&HeavyAngX); + RootOutput::getInstance()->GetTree()->Branch("HeavyAngY",&HeavyAngY); RootOutput::getInstance()->GetTree()->Branch("ThetaCM",&ThetaCM); RootOutput::getInstance()->GetTree()->Branch("Run",&Run,"Run/I"); RootOutput::getInstance()->GetTree()->Branch("X",&X); @@ -939,11 +1027,7 @@ void Analysis::InitOutputBranch(){ //////////////////////////////////////////////////////////////////////////////// void Analysis::InitInputBranch(){ - -//cout << "here_InitInput" << endl; SetBranchStatus(); - // RootInput:: getInstance()->GetChain()->SetBranchAddress("GATCONF",&vGATCONF); - // // Vamos RootInput::getInstance()->GetChain()->SetBranchAddress("LTS",<S); /* RootInput::getInstance()->GetChain()->SetBranchAddress("T_FPMW_CATS1",&T_FPMW_CATS1); @@ -996,14 +1080,10 @@ void Analysis::InitInputBranch(){ } if(isSim && !isPhaseSpace){ - //RootInput:: getInstance()->GetChain()->SetBranchStatus("InitialConditions",true ); - //RootInput:: getInstance()->GetChain()->SetBranchStatus("fIC_*",true ); RootInput::getInstance()->GetChain()->SetBranchAddress("InitialConditions", &Initial); RootInput::getInstance()->GetChain()->SetBranchAddress("InteractionCoordinates", &Interaction); - //RootInput:: getInstance()->GetChain()->SetBranchStatus("ReactionConditions",true ); - //RootInput:: getInstance()->GetChain()->SetBranchStatus("fRC_*",true ); RootInput::getInstance()->GetChain()->SetBranchAddress("ReactionConditions", &ReactionConditions); } @@ -1011,7 +1091,6 @@ void Analysis::InitInputBranch(){ //////////////////////////////////////////////////////////////////////////////// void Analysis::SetBranchStatus(){ -//cout << "here_SetBranchStatus" << endl; // Set Branch status RootInput::getInstance()->GetChain()->SetBranchStatus("LTS",true); /* RootInput::getInstance()->GetChain()->SetBranchStatus("T_FPMW_CATS1",true); @@ -1075,12 +1154,13 @@ void Analysis::SetBranchStatus(){ //////////////////////////////////////////////////////////////////////////////// void Analysis::ReInitValue(){ -//cout << "here_ReInit" << endl; Ex.clear(); Ecm.clear(); AddBack_EDC.clear(); AddBack_EDC_Event1.clear(); AddBack_EDC_Event2.clear(); + AGATA_Theta.clear(); + AGATA_Phi.clear(); AGATA_GammaPx.clear(); AGATA_GammaPy.clear(); AGATA_GammaPz.clear(); @@ -1098,8 +1178,14 @@ void Analysis::ReInitValue(){ ThetaLab.clear(); PhiLab.clear(); ThetaCM.clear(); - //relative_angle = -1000; - //sum_angle = -1000; + HeavyThetaLab.clear(); + HeavyPhiLab.clear(); + HeavyPosXAtWindow.clear(); + HeavyPosYAtWindow.clear(); + HeavyPosXAtMUST2.clear(); + HeavyPosYAtMUST2.clear(); + HeavyAngX.clear(); + HeavyAngY.clear(); X.clear(); Y.clear(); Z.clear(); @@ -1127,7 +1213,6 @@ void Analysis::ReInitValue(){ MG_Y=-1000; MG_D=-1000; - //EDC= -1000; EDC.clear(); } @@ -1152,3 +1237,4 @@ extern "C"{ proxy_analysis p_analysis; } + diff --git a/Projects/e793s/Analysis.h b/Projects/e793s/Analysis.h old mode 100755 new mode 100644 index 8b15a8c64a671bb7fa2a92c476bcf827322921b0..c022c97fa6d906cde853453c66d51ef0e0356476 --- a/Projects/e793s/Analysis.h +++ b/Projects/e793s/Analysis.h @@ -35,7 +35,6 @@ #include "TMust2Physics.h" #include "TMugastPhysics.h" #include "TCATSPhysics.h" -//#include "../../NPLib/Detectors/ModularLeaf/TModularLeafPhysics.h" #include "TModularLeafPhysics.h" #include <TRandom3.h> #include <TVector3.h> @@ -47,12 +46,18 @@ 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_MM = new TH1F("ThetaCM_detected_MM","ThetaCM_detected_MM",NumThetaAngleBins,0,180); auto ThetaCM_detected_MG = new TH1F("ThetaCM_detected_MG","ThetaCM_detected_MG",NumThetaAngleBins,0,180); +auto ThetaCM_detected_MM_wV = new TH1F("ThetaCM_detected_MM_wV","ThetaCM_detected_MM_wV",NumThetaAngleBins,0,180); +auto ThetaCM_detected_MG_wV = new TH1F("ThetaCM_detected_MG_wV","ThetaCM_detected_MG_wV",NumThetaAngleBins,0,180); +auto ThetaCM_detected_MM_Unb = new TH1F("ThetaCM_detected_MM_Unb","ThetaCM_detected_MM_Unb",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_MM = new TH1F("ThetaLab_detected_MM","ThetaLab_detected_MM",NumThetaAngleBins,0,180); auto ThetaLab_detected_MG = new TH1F("ThetaLab_detected_MG","ThetaLab_detected_MG",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); +auto ThetaLab_detected_MM_wV = new TH1F("ThetaLab_detected_MM_wV","ThetaLab_detected_MM_wV",NumThetaAngleBins,0,180); +auto ThetaLab_detected_MG_wV = new TH1F("ThetaLab_detected_MG_wV","ThetaLab_detected_MG_wV",NumThetaAngleBins,0,180); +auto ThetaLab_detected_MM_Unb = new TH1F("ThetaLab_detected_MM_Unb","ThetaLab_detected_MM_Unb",NumThetaAngleBins,0,180); +auto ThetaLab_detected_MG_Unb = new TH1F("ThetaLab_detected_MG_Unb","ThetaLab_detected_MG_Unb",NumThetaAngleBins,0,180); double degtorad = M_PI/180.; vector<double> clineVal, clineX; @@ -63,6 +68,16 @@ TH1F *ThetaCM_detected_MMX[5]; TH1F *ThetaLab_detected_MGX[6]; TH1F *ThetaLab_detected_MMX[5]; +TH1F *ThetaCM_detected_MGX_wV[6]; +TH1F *ThetaCM_detected_MMX_wV[5]; +TH1F *ThetaLab_detected_MGX_wV[6]; +TH1F *ThetaLab_detected_MMX_wV[5]; + +TH1F *ThetaCM_detected_MGX_Unb[6]; +TH1F *ThetaCM_detected_MMX_Unb[5]; +TH1F *ThetaLab_detected_MGX_Unb[6]; +TH1F *ThetaLab_detected_MMX_Unb[5]; + class Analysis: public NPL::VAnalysis{ public: Analysis(); @@ -88,6 +103,8 @@ class Analysis: public NPL::VAnalysis{ vector<double> AddBack_EDC; vector<double> AddBack_EDC_Event1; vector<double> AddBack_EDC_Event2; + vector<double> AGATA_Theta; + vector<double> AGATA_Phi; vector<double> AGATA_GammaPx; vector<double> AGATA_GammaPy; vector<double> AGATA_GammaPz; @@ -107,6 +124,17 @@ class Analysis: public NPL::VAnalysis{ std::vector<double> ThetaLab; std::vector<double> PhiLab; std::vector<double> ThetaCM; + std::vector<double> HeavyThetaLab; + std::vector<double> HeavyPhiLab; + std::vector<double> HeavyPosXAtWindow; + std::vector<double> HeavyPosYAtWindow; + std::vector<double> HeavyPosXAtMUST2; + std::vector<double> HeavyPosYAtMUST2; + std::vector<double> HeavyAngX; + std::vector<double> HeavyAngY; + + double WindowDistance = 681.33; + double MUST2InnerDistance = 391.4; //For use in simulations double OriginalELab; @@ -114,8 +142,8 @@ class Analysis: public NPL::VAnalysis{ double BeamEnergy; NPL::Reaction reaction; - //NPL::Reaction *reaction; - // Energy loss table: the G4Table are generated by the simulation + + //Enery loss tables are generated by G4 simulation NPL::EnergyLoss LightTarget; NPL::EnergyLoss LightAl; NPL::EnergyLoss LightSi; @@ -124,7 +152,6 @@ class Analysis: public NPL::VAnalysis{ double TargetThickness ; double WindowsThickness; - // Beam Energy double OriginalBeamEnergy ; // AMEV double FinalBeamEnergy; @@ -147,7 +174,7 @@ class Analysis: public NPL::VAnalysis{ double Si_E_M2; double CsI_E_M2; double Energy; -#define GATCONF_SIZE 16 + #define GATCONF_SIZE 16 unsigned int GATCONF_Counter[GATCONF_SIZE]; string GATCONF_Label[GATCONF_SIZE] = { "MUVI1", @@ -261,7 +288,6 @@ class Analysis: public NPL::VAnalysis{ // Branches and detectors TMust2Physics* M2; TMugastPhysics* MG; - //TCATSPhysics* CATS; TModularLeafPhysics* ML; bool warning=true; @@ -275,6 +301,7 @@ class Analysis: public NPL::VAnalysis{ TInitialConditions* Initial; TInteractionCoordinates* Interaction; TReactionConditions* ReactionConditions; + bool hitVAMOS = false; // Beam object NPL::Beam* Beam; @@ -285,3 +312,4 @@ class Analysis: public NPL::VAnalysis{ }; #endif + diff --git a/Projects/e793s/CMakeLists.txt b/Projects/e793s/CMakeLists.txt index 18d9ab375d9946eb246d297aacc7aee6813ff8a6..08d7404639f00cc7be6175ea8a25e9d7089a32d4 100755 --- a/Projects/e793s/CMakeLists.txt +++ b/Projects/e793s/CMakeLists.txt @@ -1,4 +1,4 @@ -cmake_minimum_required (VERSION 2.8) +cmake_minimum_required (VERSION 2.8...3.23) project(e793s) # Setting the policy to match Cmake version cmake_policy(VERSION ${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION}) diff --git a/Projects/e793s/RunToTreat.txt b/Projects/e793s/RunToTreat.txt index ec67499e059bacbdc96a5ddd41351eb0352de992..1b5db402bf6ebfa5ef639fa9ed7c26466863ac77 100755 --- a/Projects/e793s/RunToTreat.txt +++ b/Projects/e793s/RunToTreat.txt @@ -1,14 +1,14 @@ TTreeName TreeMaster RootFileName - /home/charlottepaxman/StoreData/K48/run_0051/Tree_0000.root - /home/charlottepaxman/StoreData/K48/run_0052/Tree_0000.root - /home/charlottepaxman/StoreData/K48/run_0053/Tree_0000.root - /home/charlottepaxman/StoreData/K48/run_0054/Tree_0000.root - /home/charlottepaxman/StoreData/K48/run_0055/Tree_0000.root - /home/charlottepaxman/StoreData/K48/run_0056/Tree_0000.root - /home/charlottepaxman/StoreData/K48/run_0060/Tree_0000.root - /home/charlottepaxman/StoreData/K48/run_0061/Tree_0000.root - /home/charlottepaxman/StoreData/K48/run_0062/Tree_0000.root - /home/charlottepaxman/StoreData/K48/run_0063/Tree_0000.root - /home/charlottepaxman/StoreData/K48/run_0064/Tree_0000.root + ~/StoreData/K48/run_0051/Tree_0000.root + ~/StoreData/K48/run_0052/Tree_0000.root + ~/StoreData/K48/run_0053/Tree_0000.root + ~/StoreData/K48/run_0054/Tree_0000.root + ~/StoreData/K48/run_0055/Tree_0000.root + ~/StoreData/K48/run_0056/Tree_0000.root + ~/StoreData/K48/run_0060/Tree_0000.root + ~/StoreData/K48/run_0061/Tree_0000.root + ~/StoreData/K48/run_0062/Tree_0000.root + ~/StoreData/K48/run_0063/Tree_0000.root + ~/StoreData/K48/run_0064/Tree_0000.root diff --git a/Projects/e793s/RunToTreat_PartII.txt b/Projects/e793s/RunToTreat_PartII.txt new file mode 100755 index 0000000000000000000000000000000000000000..10492073d0a1a01edc2113073760789a312201dc --- /dev/null +++ b/Projects/e793s/RunToTreat_PartII.txt @@ -0,0 +1,5 @@ +TTreeName + TreeMaster +RootFileName + ~/StoreData/K48/run_0063/Tree_0000.root + ~/StoreData/K48/run_0064/Tree_0000.root diff --git a/Projects/e793s/RunToTreat_PartI_2.txt b/Projects/e793s/RunToTreat_PartI_2.txt new file mode 100644 index 0000000000000000000000000000000000000000..5ee524b2fcb9147ee4b7ff8dc04f68db6a06bb6b --- /dev/null +++ b/Projects/e793s/RunToTreat_PartI_2.txt @@ -0,0 +1,10 @@ +TTreeName + TreeMaster +RootFileName + ~/StoreData/K48/run_0053/Tree_0000.root + ~/StoreData/K48/run_0054/Tree_0000.root + ~/StoreData/K48/run_0055/Tree_0000.root + ~/StoreData/K48/run_0056/Tree_0000.root + ~/StoreData/K48/run_0060/Tree_0000.root + ~/StoreData/K48/run_0061/Tree_0000.root + ~/StoreData/K48/run_0062/Tree_0000.root diff --git a/Projects/e793s/autophs.sh b/Projects/e793s/autophs.sh new file mode 100755 index 0000000000000000000000000000000000000000..d2d9fc4682be07dd0ddea2b37e221148ed3ea899 --- /dev/null +++ b/Projects/e793s/autophs.sh @@ -0,0 +1,8 @@ +#!/bin/bash + +#./phs.sh PhaseSpace_Thresh1p2 0000 +#./phs.sh PhaseSpace_Thresh1p4 0000 + + +#./phs.sh PhaseSpace_MegaNoThresh 0000 + diff --git a/Projects/e793s/autosim.sh b/Projects/e793s/autosim.sh index e077c83d1850d1f7580b3be0c2e146e327485397..c91bd794303c38c395beb2de8ea3ad5064579036 100755 --- a/Projects/e793s/autosim.sh +++ b/Projects/e793s/autosim.sh @@ -1,23 +1,67 @@ #!/bin/bash -#./sim.sh 47Kdp 0000 -#./sim.sh 47Kdp 0143 -#./sim.sh 47Kdp 1981 -#./sim.sh 47Kdp 1410 -#./sim.sh 47Kdp 0279 -#./sim.sh 47Kdp 0728 -#./sim.sh 47Kdp 0968 -#./sim.sh 47Kdp 2412 -#./sim.sh 47Kdp 2910 + ./sim.sh 47Kdp 0000 + ./sim.sh 47Kdp 0143 + ./sim.sh 47Kdp 0279 + ./sim.sh 47Kdp 0728 + ./sim.sh 47Kdp 0967 + ./sim.sh 47Kdp 1409 + ./sim.sh 47Kdp 1978 + ./sim.sh 47Kdp 2407 + ./sim.sh 47Kdp 2908 + ./sim.sh 47Kdp 3254 + ./sim.sh 47Kdp 3601 + ./sim.sh 47Kdp 3830 + ./sim.sh 47Kdp 4363 -#./sim.sh 47Kdt 1945 -#./sim.sh 47Kdt 3344 -#./sim.sh 47Kdt 2730 +#./sim.sh 47Kdp 4860 +#./sim.sh 47Kdp 5280 +#./sim.sh 47Kdp 5860 +###################### +#Sims for unbound fits +#./sim.sh 47Kdp 4820 +#./sim.sh 47Kdp 4840 +#./sim.sh 47Kdp 4860 +#./sim.sh 47Kdp 4880 + #around 4900 +#./sim.sh 47Kdp 4920 +#./sim.sh 47Kdp 4940 +#./sim.sh 47Kdp 4960 +#./sim.sh 47Kdp 4980 +#--------------------# +#./sim.sh 47Kdp 5220 +#./sim.sh 47Kdp 5240 +#./sim.sh 47Kdp 5260 +#./sim.sh 47Kdp 5280 + #around 5300 +#./sim.sh 47Kdp 5320 +#./sim.sh 47Kdp 5340 +#./sim.sh 47Kdp 5360 +#./sim.sh 47Kdp 5380 +#--------------------# +#./sim.sh 47Kdp 5820 +#./sim.sh 47Kdp 5840 +#./sim.sh 47Kdp 5860 +#./sim.sh 47Kdp 5880 + #around 5900 +#./sim.sh 47Kdp 5920 +#./sim.sh 47Kdp 5940 +#./sim.sh 47Kdp 5960 +#./sim.sh 47Kdp 5980 +###################### -./sim.sh 47Kdt 0000 -./sim.sh 47Kdt 1000 -./sim.sh 47Kdt 2000 -./sim.sh 47Kdt 3000 -./sim.sh 47Kdt 4000 -./sim.sh 47Kdt 5000 -./sim.sh 47Kdt 6000 +# ./sim.sh 47Kdt 0000 +### ./sim.sh 47Kdt 0587 +### ./sim.sh 47Kdt 0691 +# ./sim.sh 47Kdt 0639 +###./sim.sh 47Kdt 0886 +## ./sim.sh 47Kdt 1370 +# ./sim.sh 47Kdt 1554 +## ./sim.sh 47Kdt 1738 +# ./sim.sh 47Kdt 1944 +# ./sim.sh 47Kdt 2233 +# ./sim.sh 47Kdt 2732 +### ./sim.sh 47Kdt 4150 +# ./sim.sh 47Kdt 3377 +## ./sim.sh 47Kdt 4297 +# ./sim.sh 47Kdt 4150 diff --git a/Projects/e793s/configs/ConfigMugast.dat b/Projects/e793s/configs/ConfigMugast.dat index 4f23355a7664046a71c839adaddffd9a1b4b11dc..c14b30b918f6882a31ac8be943b63fc4ffd14656 100755 --- a/Projects/e793s/configs/ConfigMugast.dat +++ b/Projects/e793s/configs/ConfigMugast.dat @@ -1,102 +1,196 @@ ConfigMugast TAKE_E_Y= 1 - TAKE_T_Y= 1 - + MG1_TAKE_T_Y= 0 + MG2_TAKE_T_Y= 0 + MG3_TAKE_T_Y= 0 + MG4_TAKE_T_Y= 0 + MG5_TAKE_T_Y= 1 + MG7_TAKE_T_Y= 0 +ConfigMugast DISABLE_CHANNEL_X= 2 4 +ConfigMugast DISABLE_CHANNEL_X= 3 7 +ConfigMugast DISABLE_CHANNEL_X= 3 5 +ConfigMugast DISABLE_CHANNEL_X= 4 31 +ConfigMugast DISABLE_CHANNEL_X= 5 61 +ConfigMugast DISABLE_CHANNEL_X= 5 24 +ConfigMugast DISABLE_CHANNEL_X= 5 32 +ConfigMugast DISABLE_CHANNEL_X= 5 125 +ConfigMugast DISABLE_CHANNEL_X= 5 124 +ConfigMugast DISABLE_CHANNEL_X= 5 109 +ConfigMugast DISABLE_CHANNEL_X= 5 106 +ConfigMugast DISABLE_CHANNEL_X= 5 104 +ConfigMugast DISABLE_CHANNEL_X= 5 94 +ConfigMugast DISABLE_CHANNEL_X= 5 76 +ConfigMugast DISABLE_CHANNEL_X= 5 73 +ConfigMugast DISABLE_CHANNEL_X= 5 71 +ConfigMugast DISABLE_CHANNEL_X= 5 72 +ConfigMugast DISABLE_CHANNEL_X= 7 9 +ConfigMugast DISABLE_CHANNEL_X= 7 77 +ConfigMugast DISABLE_CHANNEL_X= 7 75 +ConfigMugast DISABLE_CHANNEL_X= 7 73 +ConfigMugast DISABLE_CHANNEL_X= 7 74 +ConfigMugast DISABLE_CHANNEL_X= 7 71 +ConfigMugast DISABLE_CHANNEL_X= 7 72 +ConfigMugast DISABLE_CHANNEL_X= 7 69 +ConfigMugast DISABLE_CHANNEL_X= 7 70 +ConfigMugast DISABLE_CHANNEL_X= 7 67 +ConfigMugast DISABLE_CHANNEL_X= 7 68 +ConfigMugast DISABLE_CHANNEL_X= 7 65 - +ConfigMugast DISABLE_CHANNEL_Y= 1 85 +ConfigMugast DISABLE_CHANNEL_Y= 1 31 +ConfigMugast DISABLE_CHANNEL_Y= 1 36 +ConfigMugast DISABLE_CHANNEL_Y= 1 27 +ConfigMugast DISABLE_CHANNEL_Y= 1 35 +ConfigMugast DISABLE_CHANNEL_Y= 1 24 +ConfigMugast + DISABLE_CHANNEL_Y= 1 94 +ConfigMugast DISABLE_CHANNEL_Y= 1 33 +ConfigMugast DISABLE_CHANNEL_Y= 1 39 +ConfigMugast DISABLE_CHANNEL_Y= 1 37 +ConfigMugast DISABLE_CHANNEL_Y= 1 22 +ConfigMugast DISABLE_CHANNEL_Y= 1 46 +ConfigMugast DISABLE_CHANNEL_Y= 1 26 +ConfigMugast DISABLE_CHANNEL_Y= 1 13 +ConfigMugast DISABLE_CHANNEL_Y= 1 15 +ConfigMugast DISABLE_CHANNEL_Y= 2 39 +ConfigMugast DISABLE_CHANNEL_Y= 2 16 +ConfigMugast DISABLE_CHANNEL_Y= 2 62 +ConfigMugast DISABLE_CHANNEL_Y= 3 11 +ConfigMugast DISABLE_CHANNEL_Y= 4 100 +ConfigMugast DISABLE_CHANNEL_Y= 4 90 +ConfigMugast DISABLE_CHANNEL_Y= 4 102 +ConfigMugast DISABLE_CHANNEL_Y= 4 117 +ConfigMugast DISABLE_CHANNEL_Y= 4 105 +ConfigMugast DISABLE_CHANNEL_Y= 4 38 +ConfigMugast DISABLE_CHANNEL_Y= 4 36 +ConfigMugast DISABLE_CHANNEL_Y= 4 22 +ConfigMugast DISABLE_CHANNEL_Y= 4 50 +ConfigMugast DISABLE_CHANNEL_Y= 4 26 +ConfigMugast DISABLE_CHANNEL_Y= 4 23 +ConfigMugast DISABLE_CHANNEL_Y= 4 64 +ConfigMugast DISABLE_CHANNEL_Y= 4 4 +ConfigMugast DISABLE_CHANNEL_Y= 4 2 +ConfigMugast DISABLE_CHANNEL_Y= 4 55 +ConfigMugast DISABLE_CHANNEL_Y= 5 120 +ConfigMugast DISABLE_CHANNEL_Y= 5 122 +ConfigMugast DISABLE_CHANNEL_Y= 5 113 +ConfigMugast DISABLE_CHANNEL_Y= 5 111 +ConfigMugast DISABLE_CHANNEL_Y= 5 109 +ConfigMugast DISABLE_CHANNEL_Y= 5 93 +ConfigMugast DISABLE_CHANNEL_Y= 5 107 +ConfigMugast DISABLE_CHANNEL_Y= 5 105 +ConfigMugast DISABLE_CHANNEL_Y= 5 38 +ConfigMugast DISABLE_CHANNEL_Y= 5 36 +ConfigMugast DISABLE_CHANNEL_Y= 5 4 +ConfigMugast DISABLE_CHANNEL_Y= 7 126 +ConfigMugast DISABLE_CHANNEL_Y= 7 103 +ConfigMugast DISABLE_CHANNEL_Y= 7 128 +ConfigMugast DISABLE_CHANNEL_Y= 7 123 +ConfigMugast DISABLE_CHANNEL_Y= 7 121 +ConfigMugast DISABLE_CHANNEL_Y= 7 105 +ConfigMugast DISABLE_CHANNEL_Y= 7 38 +ConfigMugast DISABLE_CHANNEL_Y= 7 36 +ConfigMugast DISABLE_CHANNEL_Y= 7 35 +ConfigMugast DISABLE_CHANNEL_Y= 7 40 +ConfigMugast DISABLE_CHANNEL_Y= 7 43 +ConfigMugast DISABLE_CHANNEL_Y= 7 46 +ConfigMugast DISABLE_CHANNEL_Y= 7 54 +ConfigMugast DISABLE_CHANNEL_Y= 7 60 +ConfigMugast DISABLE_CHANNEL_Y= 7 63 +ConfigMugast DISABLE_CHANNEL_Y= 7 57 - +ConfigMugast MAX_STRIP_MULTIPLICITY= 10 STRIP_ENERGY_MATCHING= 0.06 MeV DSSD_X_E_RAW_THRESHOLD= 8250 DSSD_Y_E_RAW_THRESHOLD= 8100 - DSSD_X_E_THRESHOLD= 1 - DSSD_Y_E_THRESHOLD= 1 + DSSD_X_E_THRESHOLD= 0.01 MeV + DSSD_Y_E_THRESHOLD= 0.01 MeV diff --git a/Projects/e793s/configs/ConfigMust2.dat b/Projects/e793s/configs/ConfigMust2.dat index beee75509e8a157ba1c76de44ccb71e3bb664191..15bbb53435c46ee4311c60ae054cbeffaea3acd7 100755 --- a/Projects/e793s/configs/ConfigMust2.dat +++ b/Projects/e793s/configs/ConfigMust2.dat @@ -149,7 +149,34 @@ ConfigMust2 DISABLE_CHANNEL MM5STRX101 DISABLE_CHANNEL MM5STRX102 DISABLE_CHANNEL MM5STRX103 - + %CHECKING HOW REMOVING THE REST WORKS!!!! + DISABLE_CHANNEL MM5STRX104 + DISABLE_CHANNEL MM5STRX105 + DISABLE_CHANNEL MM5STRX106 + DISABLE_CHANNEL MM5STRX107 + DISABLE_CHANNEL MM5STRX108 + DISABLE_CHANNEL MM5STRX109 + DISABLE_CHANNEL MM5STRX110 + DISABLE_CHANNEL MM5STRX111 + DISABLE_CHANNEL MM5STRX112 + DISABLE_CHANNEL MM5STRX113 + DISABLE_CHANNEL MM5STRX114 + DISABLE_CHANNEL MM5STRX115 + DISABLE_CHANNEL MM5STRX116 + DISABLE_CHANNEL MM5STRX117 + DISABLE_CHANNEL MM5STRX118 + DISABLE_CHANNEL MM5STRX119 + DISABLE_CHANNEL MM5STRX120 + DISABLE_CHANNEL MM5STRX121 + DISABLE_CHANNEL MM5STRX122 + DISABLE_CHANNEL MM5STRX123 + DISABLE_CHANNEL MM5STRX124 + DISABLE_CHANNEL MM5STRX125 + DISABLE_CHANNEL MM5STRX126 + DISABLE_CHANNEL MM5STRX127 + DISABLE_CHANNEL MM5STRX128 + + %MM5 Y DISABLE_CHANNEL MM5STRY0 DISABLE_CHANNEL MM5STRY1 @@ -171,6 +198,7 @@ ConfigMust2 DISABLE_CHANNEL MM5STRY113 DISABLE_CHANNEL MM5STRY122 DISABLE_CHANNEL MM5STRY127 + DISABLE_CHANNEL MM5STRY128 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% diff --git a/Projects/e793s/exp.sh b/Projects/e793s/exp.sh index e60f417cb1bb417a13fa7dd466d5f9373902d4e7..4fd3810d1dccde0f526b00bbe5a72ee7f37a8445 100755 --- a/Projects/e793s/exp.sh +++ b/Projects/e793s/exp.sh @@ -4,98 +4,20 @@ cmake ./; make -j6; #==================================================== -#rfile='Reaction/47Kdp_08Nov.reaction' -#rfile='Reaction/47Kdd_08Nov.reaction' -#rfile='Reaction/47Kpp_08Nov.reaction' -#rfile='Reaction/47Kdt_08Nov.reaction' -#---------------------------------------------------- -#rfile='Reaction/47Kdp_11Jul22.reaction' -#rfile='Reaction/47Kdt_11Jul22.reaction' -#rfile='Reaction/47Kdd_11Jul22.reaction' -#rfile='Reaction/47Kpp_11Jul22.reaction' -#---------------------------------------------------- -#rfile='Reaction/47Kdp_30Aug22_2p06um.reaction' -#---------------------------------------------------- -#rfile='Reaction/47Kdp_14Oct22.reaction' -#rfile='Reaction/47Kdp_14Oct22_2.reaction' -#---------------------------------------------------- -rfile='Reaction/47Kdp_18Oct22.reaction' -#rfile='Reaction/47Kdt_18Oct22.reaction' +##rfile='./Reaction/47Kdp_zerozero.reaction' + rfile='./Reaction/47Kdp_18Oct22.reaction' +##rfile='./Reaction/47Kpp_18Oct22.reaction' +##rfile='./Reaction/47Kdd_18Oct22.reaction' +##rfile='./Reaction/47Kdt_18Oct22.reaction' +##-------------------------- +##rfile='./Reaction/47Kdp_25Jul23.reaction' +##rfile='./Reaction/47Kdd_25Jul23.reaction' #==================================================== -#dfile='Detector/mugast_08Nov.detector' -#---------------------------------------------------- -#dfile='Detector/mugast_11Jul22.detector' -#dfile='Detector/mugast_Test2um.detector' -#dfile='Detector/mugast_17Jul22_MM+200mm.detector' -#dfile='Detector/mugast_17Jul22_MM+203mm.detector' -#dfile='Detector/mugast_17Jul22_MM+210mm.detector' -#dfile='Detector/mugast_17Jul22_MM+190mm.detector' #DO ME!!! -#---------------------------------------------------- -#dfile='Detector/mugast_01Sep22.detector' -#---------------------------------------------------- -#dfile='Detector/mugast_14Oct22.detector' -#dfile='Detector/mugast_14Oct22_2.detector' -dfile='Detector/mugast_18Oct22.detector' +##dfile='./Detector/mugast_zerozero.detector' + dfile='./Detector/mugast_18Oct22.detector' +##-------------------------- +##dfile='./Detector/mugast_25Jul23.detector' #==================================================== - -smallTest='_Test' -ptI='_PartI' -ptII='_PartII' -ptA='_PartA' -ptB='_PartB' -ptC='_PartC' - -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' - -#npanalysis --definition Exp -R RunToTreat_PartI.txt -E $rfile -D $dfile -C Calibration.txt -O $1$ptI; -#npanalysis --definition Exp -R RunToTreat_PartI_2.txt -E $rfile -D $dfile -C Calibration.txt -O $1$ptI; -#npanalysis --definition Exp -R RunToTreat_PartII.txt -E $rfile -D $dfile -C Calibration.txt -O $1$ptII; -npanalysis --definition Exp -R RunToTreat_PartA.txt -E $rfile -D $dfile -C Calibration.txt -O $1$ptA; -npanalysis --definition Exp -R RunToTreat_PartB.txt -E $rfile -D $dfile -C Calibration.txt -O $1$ptB; -npanalysis --definition Exp -R RunToTreat_PartC.txt -E $rfile -D $dfile -C Calibration.txt -O $1$ptC; - -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -echo 'CURRENTLY EXCLUDING RUNS 51 AND 52!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' -#npanalysis --definition Exp --definition ExcludeThePoor -R RunToTreat_PartI.txt -E $rfile -D $dfile -C Calibration.txt -O $1$ptI; -#npanalysis --definition Exp --definition ExcludeThePoor -R RunToTreat_PartII.txt -E $rfile -D $dfile -C Calibration.txt -O $1$ptII; + npanalysis --definition Exp -R RunToTreat_PartI_2.txt -E $rfile -D $dfile -C Calibration.txt -O $1$ptI; + npanalysis --definition Exp -R RunToTreat_PartII.txt -E $rfile -D $dfile -C Calibration.txt -O $1$ptII; diff --git a/Projects/e793s/macro/BeamSpot/EventReader.C b/Projects/e793s/macro/BeamSpot/EventReader.C index 227da085011d85ff770880f3b38105eb5477bc3c..4dbfbf4a26b22a4ed0f11606c5d7794ec13b12dc 100755 --- a/Projects/e793s/macro/BeamSpot/EventReader.C +++ b/Projects/e793s/macro/BeamSpot/EventReader.C @@ -8,27 +8,26 @@ #include <fstream> #include <vector> #include <cmath> -//Root -//#include <TVector3.h> -//NPTool -//#include "NPEnergyLoss.h" -//#include "NPReaction.h" -//#include "NPSystemOfUnits.h" using namespace std; //------------------------------- void EventReader(){ // Read ROOT file and pull the tree - //auto DataFile = new TFile("../../../../Outputs/Analysis/47K_Full_08July.root", "READ"); - auto DataFile = new TFile("../../../../Outputs/Analysis/47Kdp_11Apr22_PartII.root", "READ"); + int breakEvent; + string outname; + auto DataFile = new TFile("../../../../Outputs/Analysis/47Kdp_25Jul23_5p5um_PartI.root", "READ"); breakEvent = 240178; outname = "./XYZE_5p5um_09Apr24_PartI.txt"; + //auto DataFile = new TFile("../../../../Outputs/Analysis/47Kdp_25Jul23_5p5um_PartII.root", "READ"); breakEvent = 1000000000000; outname = "./XYZE_5p5um_09Apr24_PartII.txt"; + auto PhysTree = (TTree*) DataFile->FindObjectAny("PhysicsTree"); // Initilise the Mugast branch auto Mugast = new TMugastPhysics(); // Initilise access variables for PhysTree + static vector<double> *MWT; static double T_MUGAST_VAMOS; + static double T_MUGAST_CATS2; static vector<double> *X, *Y, *Z, *RawEnergy, *AddBack_EDC, *ThetaLab; // Pull PhysTree branches @@ -37,7 +36,9 @@ void EventReader(){ auto X_Branch = PhysTree->GetBranch("X"); auto Y_Branch = PhysTree->GetBranch("Y"); auto Z_Branch = PhysTree->GetBranch("Z"); + auto MWT_Branch = PhysTree->GetBranch("MWT"); auto MugVam_Branch = PhysTree->GetBranch("T_MUGAST_VAMOS"); + auto MugCat_Branch = PhysTree->GetBranch("T_MUGAST_CATS2"); auto ThetaLab_Branch = PhysTree->GetBranch("ThetaLab"); // Set Mugast branch address @@ -49,7 +50,9 @@ void EventReader(){ X_Branch->SetAddress(&X); Y_Branch->SetAddress(&Y); Z_Branch->SetAddress(&Z); + MWT_Branch->SetAddress(&MWT); MugVam_Branch->SetAddress(&T_MUGAST_VAMOS); + MugCat_Branch->SetAddress(&T_MUGAST_CATS2); ThetaLab_Branch->SetAddress(&ThetaLab); // Build loop variables @@ -58,41 +61,31 @@ void EventReader(){ // Open output file ofstream outfile; - //outfile.open("./XYZE_ShiftMG3_Sep08_GateOn0143.txt"); - outfile.open("./XYZE_GateOn1838_11Apr22_PtII.txt"); + outfile.open(outname.c_str()); // Loop on entries - for(unsigned int i=0; i<numEntries; i++){ + for(unsigned int i=0; i<breakEvent; i++){ PhysTree->GetEntry(i); multiplicity = Mugast->TelescopeNumber.size(); // Loop on MUGAST multiplicity for(int m=0; m<multiplicity; m++){ - // Gate on Timing - if(abs(T_MUGAST_VAMOS-2700)<400){ + //if(abs(T_MUGAST_VAMOS-2700)<400){ int gammaMultip = AddBack_EDC->size(); -// cout << "GAMMA MULT " << gammaMultip << endl; - //if(ThetaLab->at(m) <= 130. && ThetaLab->at(m) >= 100.){ - - // gamma stuff (fails otherwise??) - //if(gammaMultip!=0){ if(gammaMultip>=1){ - //for(int g; g<gammaMultip; g++){ - // Gate on E_Gamma - if(abs(AddBack_EDC->at(0)-1.838) < 0.005){ + if(abs(AddBack_EDC->at(0)-0.143) < 0.005){ outfile << Mugast->TelescopeNumber.at(m) << " " << X->at(m) << " " << Y->at(m) << " " << Z->at(m) << " " << RawEnergy->at(m) << endl; }//if 0.143 - //}//for g - }//if gamma + }//if gammamultip //}//if theta (new!! testing!!! - }//if timing + //}//if timing }//for m }//for i outfile.close(); diff --git a/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.cxx b/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.cxx index 35521abf090b546c16505c657753b93c9e7cc6e3..e2984121dcc2bad15a27c53bb17c2c5664ea1203 100755 --- a/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.cxx +++ b/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.cxx @@ -34,6 +34,13 @@ double devE(const double * parameter) { h5->SetAxisRange(-1.0, +2.6, "X"); h7->SetAxisRange(-1.0, +2.6, "X"); } + //Set final beam energy (added 05July) + double FinalBeamEnergy = BeamTarget.Slow( + 354.75, //post-CATS beam energy + 0.5*parameter[3]*micrometer, //thickness + 0); //angle, probably + reaction.SetBeamEnergy(FinalBeamEnergy); + //Initilize results array // 7 => Sum in 0 and them MG's in 1-6 @@ -46,14 +53,6 @@ double devE(const double * parameter) { //Particle path vector dir = * (pos[i]) - offset; - //Set final beam energy (added 05July) - double FinalBeamEnergy = BeamTarget.Slow( -// reaction.GetBeamEnergy(), - 354.75, //post-CATS beam energy - 0.5*parameter[3]*micrometer, //thickness - 0); //angle, probably - reaction.SetBeamEnergy(FinalBeamEnergy); - //Define normal vector for the MG# of detection DetectorSwitch(detnum[i], MugastNormal); @@ -61,7 +60,7 @@ double devE(const double * parameter) { double ThetaTarget = dir.Angle(TVector3(0.0, 0.0, 1.0)); double ThetaMugast = dir.Angle(MugastNormal); double Energy = energy[i]; -//cout << "@ angle " << ThetaTarget/deg << " " << Energy << " -> "; + //Energy loss in Al Energy = Al.EvaluateInitialEnergy( Energy, //energy Al @@ -69,7 +68,6 @@ double devE(const double * parameter) { ThetaMugast //angle impinging on MUGAST ); -//cout << Energy << " -> "; //Energy loss in target Energy = CD2.EvaluateInitialEnergy( Energy, //energy after leaving target @@ -77,22 +75,15 @@ double devE(const double * parameter) { ThetaTarget //angle leaving target ); -//cout << Energy << endl; //Final value of Ex double Ex = reaction.ReconstructRelativistic(Energy, ThetaTarget); //Fill histograms with -//if(ThetaTarget/deg<130.){ // TESTING -//if(ThetaTarget/deg>130.){ // TESTING - if(allButMG3){ - if(detnum[i]!=3){ - h -> Fill(Ex); - cout << "removed MG3 from sum!!" << endl; - } - } else { - h -> Fill(Ex); + if(ThetaTarget/deg<120.){ + if(allButMG3){ if(detnum[i]!=3){ h -> Fill(Ex); } } + else { h -> Fill(Ex); } } -//} + DetectorSwitch(detnum[i], Ex); hT -> Fill(ThetaTarget/deg,Ex); hEL -> Fill(ThetaTarget/deg,Energy-energy[i]);//ELoss in target and Al @@ -104,18 +95,28 @@ double devE(const double * parameter) { //Write vals to screen if (flagDraw) {cout << "==================================================" << endl;} - /*** Minimize by one peak ***/ + /*** Minimize by total peak ***/ /**/ - double multiplier = 0.05; //0.08; - double metric = abs(FitResultMatrix[mgSelect][0]-refE) + abs(multiplier*FitResultMatrix[mgSelect][2]); - //double metric = abs(FitResultMatrix[mgSelect][0]-0.143) + abs(multiplier*FitResultMatrix[mgSelect][2]) + abs(FitResultMatrix[mgSelect][5]-1.410) + abs(multiplier*FitResultMatrix[mgSelect][6]); - //double metric = abs(FitResultMatrix[mgSelect][0]-0.143) + abs(multiplier*FitResultMatrix[mgSelect][2]) + abs(FitResultMatrix[mgSelect][5]-1.410) + abs(multiplier*FitResultMatrix[mgSelect][6])+ abs(FitResultMatrix[mgSelect][7]-1.980) + abs(multiplier*FitResultMatrix[mgSelect][8]) ; - //double metric = abs(FitResultMatrix[mgSelect][7]-1.981) + abs(multiplier*FitResultMatrix[mgSelect][8]); ; + //double multiplier = 0.08; //0.08; + double multiplier = 0.10; //0.08; + double metric = 0.0 + + 4.0 * abs(FitResultMatrix[mgSelect][0]-0.143) + + abs(multiplier*FitResultMatrix[mgSelect][2]) + + 3.0 * abs(FitResultMatrix[mgSelect][5]-1.410) + + abs(multiplier*FitResultMatrix[mgSelect][6]) + + 1.0 * abs(FitResultMatrix[mgSelect][7]-1.978) + + abs(multiplier*FitResultMatrix[mgSelect][8]) + ; /**/ +/*** + // Minimize by making distance from 0.143 to 1.981 = 1.838 // + double metric = 10 * abs(abs(FitResultMatrix[mgSelect][0]-FitResultMatrix[mgSelect][7]) - 1.838); +***/ + /*** Minimize by all peaks ***/ -/** +/*** //double multiplier = 0.125; double multiplier = 0.05; double metric = @@ -125,26 +126,33 @@ double devE(const double * parameter) { + (1.0/6.0)*abs(FitResultMatrix[4][0]-refE) + (1.0/6.0)*abs(FitResultMatrix[5][0]-refE) + (1.0/6.0)*abs(FitResultMatrix[6][0]-refE) - //+ 1.0*abs(FitResultMatrix[0][0]-refE) - + multiplier*FitResultMatrix[1][2] - + multiplier*FitResultMatrix[2][2] - + multiplier*FitResultMatrix[3][2] - + multiplier*FitResultMatrix[4][2] - + multiplier*FitResultMatrix[5][2] - + multiplier*FitResultMatrix[6][2] + + 1.0*abs(FitResultMatrix[0][0]-refE) + //+ multiplier*FitResultMatrix[1][2] + //+ multiplier*FitResultMatrix[2][2] + //+ multiplier*FitResultMatrix[3][2] + //+ multiplier*FitResultMatrix[4][2] + //+ multiplier*FitResultMatrix[5][2] + //+ multiplier*FitResultMatrix[6][2] ; -**/ +***/ /*** Minimize by variation in peaks ***/ -/** +/*** vector<double> r = {FitResultMatrix[1][0], FitResultMatrix[2][0], FitResultMatrix[3][0], FitResultMatrix[4][0], FitResultMatrix[5][0], - FitResultMatrix[6][0]}; + FitResultMatrix[6][0], + FitResultMatrix[1][7], + FitResultMatrix[2][7], + FitResultMatrix[3][7], + FitResultMatrix[4][7], + FitResultMatrix[5][7], + FitResultMatrix[6][7]}; + double metric = (*max_element(r.begin(),r.end()) - *min_element(r.begin(),r.end())); -**/ +***/ WriteToCout(FitResultMatrix[mgSelect], parameter[3], metric); if (flagDraw) { WriteToFile(FitResultMatrix[mgSelect], parameter, metric); } @@ -154,15 +162,13 @@ double devE(const double * parameter) { //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// void MinimizeBeamSpot() { - //filename = "XYZE_GateOn1838_11Apr22.txt"; refE = 1.981; // the energy of the selected states - filename = "XYZE_GateOn0143_11Apr22.txt"; refE = 0.143; // the energy of the selected states - + filename = "XYZE_GateOn0143_25Jul23.txt"; refE = 0.143; // the energy of the selected states //Read data in LoadFile(); //Output formatting - cout << fixed << showpoint << setprecision(6) << showpos; + cout << fixed << showpoint << setprecision(4) << showpos; //Read in cout << "==================================================" << endl; @@ -190,7 +196,7 @@ void MinimizeBeamSpot() { //Open output file file.open("gridMinResults.txt", ios::out); file << "minimX \tminimY \tminimZ \tminimT \tMetric \tMean \tMeanErr \tStdDev \tStdDvErr\tChi2/NDF" << endl; - file << setprecision(6) << fixed; + file << setprecision(4) << fixed; //Start timer auto start = high_resolution_clock::now(); @@ -198,35 +204,8 @@ void MinimizeBeamSpot() { //Start with beam (0,0,0) and 4.76um 0.5mg/c2 target double parameter[4] = { - //0.0, 0.0, 0.0, 4.76 - //-3.9164, +0.0550, 1.3958, 1.3008 - - //-3.9164, +0.0550, +1.0, +2.4 - //-3.9164, +0.0550, +0.681611, +2.256769 - //-3.9164, +0.0550, +0.786721, +2.113878 - - //-3.9164, 0.0550, 0.5, 2.6 - //-3.9164, 0.0550, +0.146952, +2.517784 - - //-4.675715, +0.143686, +0.115426, +2.586484 - //-4.636827, -0.206620, +0.115426, +2.586484 - - - //-4., +1.0, +0.115426, +2.586484 - //-5.088554, +1.476001, +0.342552, +2.586484 - //-5.088554, +1.476001, +0.236014, +2.523850 - //-3.989139, 0.99116, 0.115426, 2.586484 - //-3.989139, +0.991160, +0.587660, +2.586484 - - //-3.989139, +0.991160, +0.587660, +2.586484 - //-3.9164, +0.0550, 0.0, 2.6 - //-3.9164, +0.0550, 0.0, +2.838327 - //-3.912110, +0.083580, +0.000000, +2.838327 - //-3.912110, +0.083580, +0.000000, +2.6 - //-3.9164, +0.0550, +0.0, 2.3 - //-3.9164, +0.0550, +0.0, +2.591512 - -5.334864, -0.698900, +0.0, +2.591512 - + //+0.0, +0.0, +0.0, 3.00 + -4.0, +0.5, +0.0, 3.00 }; //Don't draw iterations of minimizer @@ -250,16 +229,15 @@ void MinimizeBeamSpot() { minim -> SetFunction(func); //Assign variable limits - //minim -> SetLimitedVariable(0, "X", parameter[0], 0.10, -8.0, -0.0); minim -> SetFixedVariable(0, "X", parameter[0]); - //minim -> SetLimitedVariable(1, "Y", parameter[1], 0.10, -6.0, +6.0); minim -> SetFixedVariable(1, "Y", parameter[1]); - //minim -> SetLimitedVariable(2, "Z", parameter[2], 0.05, 0.90, +1.10); - //minim -> SetLimitedVariable(2, "Z", parameter[2], 0.10, -2.00, +2.00); - minim -> SetFixedVariable(2, "Z", parameter[2]); - //minim -> SetLimitedVariable(3, "T", parameter[3], 0.05, +0.9, +1.5); // ELASTICS, 1.2(3) - //minim -> SetLimitedVariable(3, "T", parameter[3], 0.05, +2.05, +3.25); // ELASTICS, 2.65(59) - minim -> SetFixedVariable(3, "T", parameter[3]); + minim -> SetLimitedVariable(2, "Z", parameter[2], 0.01, -2.0, +2.0); + minim -> SetLimitedVariable(3, "T", parameter[3], 0.01, +2.7, +3.1); // ELASTICS + +//minim -> SetLimitedVariable(0, "X", parameter[0], 0.02, -5.5, -3.5); +//minim -> SetLimitedVariable(1, "Y", parameter[1], 0.02, -2.0, +2.0); +//minim -> SetFixedVariable(2, "Z", parameter[2]); +//minim -> SetFixedVariable(3, "T", parameter[3]); //Don't draw iterations of minimizer diff --git a/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.h b/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.h index 4fef196796bed609af46879d7d0807ba12522d75..d11f73582e72905a296bc0cebac2f908cc22df3a 100755 --- a/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.h +++ b/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.h @@ -6,6 +6,7 @@ #include "TError.h" #include <iostream> #include <chrono> +#include "../DefineColours.h" using namespace std; using namespace std::chrono; @@ -13,7 +14,6 @@ using namespace std::chrono; //////////////////////////////////////////////////////////////////////////////// /* Global */ //Various numbers and objects - vector<TVector3*> pos; vector<double> energy; vector<int> detnum; @@ -32,7 +32,7 @@ int writeCount = 0; //Histograms string filename; double refE; // the energy of the selected states -static auto h = new TH1D("h","All MG#'s", 200,-1.0,3.0); +static auto h = new TH1D("h","All MG#'s", 80,-1.0,3.0); static auto h1 = new TH1D("h1","Individual MG#'s", 80,-1.0,3.0); static auto h2 = new TH1D("h2","h2", 80,-1.0,3.0); static auto h3 = new TH1D("h3","h3", 80,-1.0,3.0); @@ -42,14 +42,10 @@ static auto h7 = new TH1D("h7","h7", 80,-1.0,3.0); static auto hT = new TH2F("hT","hT", 60,100.,160.,80,-1.0,3.0); static auto hEL = new TH2F("hEL","hEL", 60,100.,160.,2000,0.0,0.1); - - - +//////////////////////////////////////////////////////////////////////////////// Double_t f_full(Double_t *x, Double_t *par) { float xx = x[0]; double result, norm; - // Flat background - //result = par[0]; result = 0; // Add N peaks for(int pk=0; pk<3; pk++){ @@ -58,7 +54,7 @@ Double_t f_full(Double_t *x, Double_t *par) { } return result; } - +//////////////////////////////////////////////////////////////////////////////// Double_t f_one(Double_t *x, Double_t *par) { float xx = x[0]; double result, norm; @@ -66,8 +62,6 @@ Double_t f_one(Double_t *x, Double_t *par) { * exp(-0.5*pow((xx-par[2])/par[1],2)); return result; } - -//static auto hT = new TH2F("hT","hT", 20,100.,160.,20,-1.0,1.0); //////////////////////////////////////////////////////////////////////////////// void LoadFile(){ // Open XYZE gamma-gated file @@ -84,10 +78,6 @@ void LoadFile(){ double x,y,z,e; while(file >> mg >> x >> y >> z >> e ){ auto p = new TVector3(x,y,z); -// if(mg==3){ -// p->SetZ(z+1.0); //Edit MG3 position directly -// cout << "EDITING MG3 POSITION" << endl; -// } detnum.push_back(mg); pos.push_back(p); energy.push_back(e); @@ -171,29 +161,34 @@ void DetectorSwitch(int MG, double Ex){ //////////////////////////////////////////////////////////////////////////////// void WriteToCout(double* result, double thick, double metric){ cout - << "Mean: " + //<< " Thick: " + //<< thick + //<< " um " + << " Chi2/NDF = " + << result[4] + << " Metric: " + << metric + << " |" + << RED + << " Mean1: " << result[0] << " +/- " << result[1] - << " StdDev: " + << " StdDev1: " << result[2] << " +/- " << result[3] - << " Thick: " - << thick - << " um" - << " Chi2/NDF = " - << result[4] - << " Metric: " - << metric - << " Mean2: " + << BLUE + << " Mean2: " << result[5] - << " StdDev2: " + << " StdDev2: " << result[6] - << " Mean3: " + << GREEN + << " Mean3: " << result[7] - << " StdDev3: " + << " StdDev3: " << result[8] + << RESET << endl; } //////////////////////////////////////////////////////////////////////////////// @@ -236,31 +231,13 @@ void DrawOneHistogram(TH1D* hist, int mg, int colour, int fill, double *FitResul const char* settings; if (drawFit){ - settings = "RBWQS"; + //settings = "RBWQS"; + settings = "RBQS"; } else { - settings = "RBWQSN"; + //settings = "RBWQSN"; + settings = "RBQSN"; } -/* - //TFitResultPtr fit; - //if(abs(refE-0.143)<0.05){ - TF1 *full = new TF1("fitThreePeaks", f_full, -1.0, +2.6, (int) 1+(3*3)); - for(int i=0; i<3; i++) { - full->SetParameter((i*3)+1,0.14); - full->SetParameter((i*3)+3,1e3); - } - full->SetParameter((0*3)+2,0.143); - full->SetParameter((1*3)+2,1.410); - full->SetParameter((2*3)+2,1.981); - //fit = hist->Fit(full, settings, "", -1.0, +2.6); - TFitResultPtr fit = hist->Fit(full, settings, "", -1.0, +2.6); - //} else { - // TF1 *one = new TF1("fitOne", f_one, -1.0, +2.6, (int) 1+(3*3)); - // one->SetParameter(1,0.14); - // one->SetParameter(3,1e3); - // one->SetParameter(2,refE); - // fit = hist->Fit(one, settings, "", -1.0, +2.6); - //} -*/ + TF1 *full = new TF1("fitThreePeaks", f_full, -1.0, +2.6, (int) 1+(3*3)); for(int i=0; i<3; i++) { full->SetParameter((i*3)+1,0.14); @@ -279,8 +256,6 @@ void InitiliseCanvas(double FitResultMatrix[7][9]){ //Canvas setup TCanvas *canv = new TCanvas("canv","Ex Histograms",20,20,1600,800); gStyle->SetOptStat(0); - //canv->Divide(2,1,0.005,0.005,0); - //canv->Divide(3,1,0.005,0.005,0); canv->Divide(2,2,0.005,0.005,0); canv->cd(1)->SetLeftMargin(0.15); canv->cd(1)->SetBottomMargin(0.15); @@ -318,10 +293,18 @@ void InitiliseCanvas(double FitResultMatrix[7][9]){ DrawOneHistogram(h7, 6, 880, 0, FitResultMatrix[6], 0); canv->Update(); - TLine *line=new TLine(refE,0.0,refE,150.0); + TLine *line=new TLine(0.143,0.0,0.143,150.0); + TLine *linex=new TLine(1.981,0.0,1.981,150.0); + TLine *linea=new TLine(1.410,0.0,1.410,150.0); line->SetLineColor(kBlack); + linex->SetLineColor(kBlack); + linea->SetLineColor(kBlack); line->SetLineStyle(7); + linex->SetLineStyle(7); + linea->SetLineStyle(7); line->Draw(); + linex->Draw(); + linea->Draw(); //Format legend auto legend = new TLegend(0.15,0.7,0.45,0.9); @@ -331,7 +314,6 @@ void InitiliseCanvas(double FitResultMatrix[7][9]){ legend->AddEntry(h4,"MG #4","f"); legend->AddEntry(h5,"MG #5","f"); legend->AddEntry(h7,"MG #7","f"); - //legend->SetTextSize(20); legend->Draw(); //Histogram setup - Sum @@ -343,10 +325,18 @@ void InitiliseCanvas(double FitResultMatrix[7][9]){ DrawOneHistogram(h, 0, 1, 0, FitResultMatrix[0],1); canv->Update(); - TLine *line2=new TLine(refE,0.0,refE,h->GetMaximum()+10.); + TLine *line2=new TLine(0.143,0.0,0.143,h->GetMaximum()+10.); + TLine *liney=new TLine(1.981,0.0,1.981,h->GetMaximum()+10.); + TLine *lineb=new TLine(1.410,0.0,1.410,h->GetMaximum()+10.); line2->SetLineColor(kBlack); + liney->SetLineColor(kBlack); + lineb->SetLineColor(kBlack); line2->SetLineStyle(7); + liney->SetLineStyle(7); + lineb->SetLineStyle(7); line2->Draw(); + liney->Draw(); + lineb->Draw(); canv->cd(3); hT->GetXaxis()->SetTitle("Theta (degrees)"); diff --git a/Projects/e793s/macro/CS2_master.h b/Projects/e793s/macro/CS2_master.h index de44b9d098bc4e68f31847d55a4252f64ed1a4af..f4769642caebae4349b46d9b6d92e80dcb793a1b 100644 --- a/Projects/e793s/macro/CS2_master.h +++ b/Projects/e793s/macro/CS2_master.h @@ -49,6 +49,14 @@ string statename; string inputdate; +/////////////////////////// +double testGrad = 0.0109; +double testCons = -0.8263; +/////////////////////////// + +string chooseVAMOSoption = "original"; +//string chooseVAMOSoption = "newgating"; + //////////////////////////////////////////////////////////////////////////////// void WriteToFile(TGraphErrors* gr, string filename){ int npoints = gr->GetN(); @@ -104,7 +112,6 @@ TGraphErrors* GetExpDCS_ForFig(string filename){ return g; } -//////////////////////////////////////////////////////////////////////////////// TGraph* GetSimDCS_ForFig(string filename){ vector<double> x, dx, y, dy; ifstream f(filename.c_str()); @@ -127,6 +134,8 @@ TGraph* GetSimDCS_ForFig(string filename){ return g; } +//////////////////////////////////////////////////////////////////////////////// + void CanvasPartition(TCanvas *C,const Int_t Nx = 2,const Int_t Ny = 2, Float_t lMargin = 0.15, Float_t rMargin = 0.05, Float_t bMargin = 0.15, Float_t tMargin = 0.05); @@ -501,6 +510,122 @@ void FigureTest(){ C->SaveAs("FigureCS2s.svg"); } +void FigureTest_dt(string state){ + gStyle->SetOptStat(0); + +// auto C = (TCanvas*) gROOT->FindObject("C"); +// if (C) delete C; +// auto C = new TCanvas("C","canvas",700,700); +// C->SetFillStyle(4000); + + TCanvas *canv = new TCanvas("cEg","cEg",700,2000); + canv->Divide(2,3); + +// // Margins +// Float_t lMargin = 0.15; +// Float_t rMargin = 0.01; +// Float_t bMargin = 0.15; +// Float_t tMargin = 0.01; + + // Dummy histogram. + string baseEx = "./CS2_TextFiles/22Sep23_Exp_dt_"; + string baseTs = "./CS2_TextFiles/22Sep23_ThryS_dt_"; + string baseTp = "./CS2_TextFiles/22Sep23_ThryP_dt_"; + string baseTd = "./CS2_TextFiles/22Sep23_ThryD_dt_"; + string baseTf = "./CS2_TextFiles/22Sep23_ThryF_dt_"; + string baseTm = "./CS2_TextFiles/22Sep23_ThryM_dt_"; + + cout << "here!!!!!!!!!!!!!!!!!!!!!!!!!" << endl; + // Get the pads previously created. + //TPad* pad = (TPad*) C->FindObject(TString::Format("pad_%d_%d",1,1).Data()); + //TPad* pad = new TPad(); + TPad* pad = (TPad*) canv->GetPad(1); + pad->Draw(); + pad->SetFillStyle(4000); + pad->SetFrameFillStyle(4000); + pad->SetLogy(); + pad->cd(); + + // Size factors + Float_t xFactor = pad->GetAbsWNDC()/pad->GetAbsWNDC(); + Float_t yFactor = pad->GetAbsHNDC()/pad->GetAbsHNDC(); + + vector<double> xTp, xTf, yTp, yTf; + + string sEx = baseEx + state + ".txt"; + string sTs = baseTs + state + ".txt"; + string sTp = baseTp + state + ".txt"; + string sTd = baseTd + state + ".txt"; + string sTf = baseTf + state + ".txt"; + string sTm = baseTm + state + ".txt"; + + cout << sEx << endl; + TGraphErrors* gEx = GetExpDCS_ForFig(sEx); + gEx->SetMarkerStyle(7); + gEx->SetTitle(""); + + cout << sTs << endl; + TGraphErrors* gTs = GetExpDCS_ForFig(sTs); + gTs->SetLineColor(kMagenta); + gTs->SetTitle(""); + + cout << sTp << endl; + TGraphErrors* gTp = GetExpDCS_ForFig(sTp); + gTp->SetLineColor(kRed); + gTp->SetTitle(""); + + cout << sTd << endl; + TGraphErrors* gTd = GetExpDCS_ForFig(sTd); + gTd->SetLineColor(kGreen); + gTd->SetTitle(""); + + cout << sTf << endl; + TGraphErrors* gTf = GetExpDCS_ForFig(sTf); + gTf->SetLineColor(kBlue); + gTf->SetTitle(""); + + TGraphErrors* gTm; + if(state == "3377"){// || state == "4150"){ + cout << sTm << endl; + gTm = GetExpDCS_ForFig(sTm); + gTm->SetLineColor(kCyan); + gTm->SetTitle(""); + } + + // Format for y axis + gTs->GetYaxis()->SetLabelFont(43); + gTs->GetYaxis()->SetLabelSize(16); + gTs->GetYaxis()->SetLabelOffset(0.02); + gTs->GetYaxis()->SetTitle(""); + gTs->GetYaxis()->SetTitleFont(43); + gTs->GetYaxis()->SetTitleSize(25); + gTs->GetYaxis()->SetTitleOffset(6); + gTs->GetYaxis()->CenterTitle(); + gTs->GetYaxis()->SetRangeUser(5e-5,2e-2); + gTs->GetYaxis()->SetTitle("d#sigma/d#Omega [mb/msr]"); + + // Format for x axis + gTs->GetXaxis()->SetLabelFont(43); + gTs->GetXaxis()->SetLabelSize(16); + gTs->GetXaxis()->SetLabelOffset(0.02); + gTs->GetXaxis()->SetTitleFont(43); + gTs->GetXaxis()->SetTitleSize(25); + gTs->GetXaxis()->SetTitleOffset(5); + gTs->GetXaxis()->CenterTitle(); + gTs->GetXaxis()->SetRangeUser(2.,15.); + gTs->GetXaxis()->SetTitle("#theta_{CM} [deg]"); + + gTs->Draw("ac"); + gTp->Draw("same c"); + gTd->Draw("same c"); + gTf->Draw("same c"); + if(state == "3377"){// || state == "4150"){ + gTm->Draw("same c"); + } + gEx->Draw("same P"); + +} + void FigureTest_dt(){ gStyle->SetOptStat(0); @@ -510,9 +635,9 @@ void FigureTest_dt(){ C->SetFillStyle(4000); // Number of PADS - const Int_t Nx = 2; + const Int_t Nx = 2;//2; //const Int_t Ny = 5; - const Int_t Ny = 3; + const Int_t Ny = 3;//3; // Margins Float_t lMargin = 0.15; @@ -524,16 +649,21 @@ void FigureTest_dt(){ CanvasPartition(C,Nx,Ny,lMargin,rMargin,bMargin,tMargin); // Dummy histogram. - string baseEx = "./Exp_dt"; - string baseTp = "./ThryS_dt"; - string baseTf = "./ThryD_dt"; - string baseTm = "./ThryF_dt"; + //string baseEx = "./Exp_dt"; + //string baseTp = "./ThryS_dt"; + //string baseTf = "./ThryD_dt"; + //string baseTm = "./ThryF_dt"; //string baseTp = "./CS2_TextFiles/24Mar22_TheoryP_"; //string baseTf = "./CS2_TextFiles/24Mar22_TheoryF_"; //string baseTm = "./CS2_TextFiles/24Mar22_TheoryMixed_"; + string baseEx = "./CS2_TextFiles/22Sep23_Exp_dt"; + string baseTs = "./CS2_TextFiles/22Sep23_ThryS_dt"; + string baseTp = "./CS2_TextFiles/22Sep23_ThryP_dt"; + string baseTd = "./CS2_TextFiles/22Sep23_ThryD_dt"; + string baseTf = "./CS2_TextFiles/22Sep23_ThryF_dt"; vector<string> states = - {"3344", "3344", "3344"}; + {"1945", "1945", "1945"}; vector<int> l = //{ 3, 3, 1, 2, 1, 1, 1, 1}; //{1, 3}; @@ -568,12 +698,6 @@ void FigureTest_dt(){ Float_t yFactor = pad[0][0]->GetAbsHNDC()/pad[i][j]->GetAbsHNDC(); - -// int i,j; -// if(p%2==0){ i=0; } -// else { i=1; } -// j=(p-i)/2; - int p = (2*j)+i; cout << p << " " << i << " " << j << endl; @@ -581,9 +705,10 @@ void FigureTest_dt(){ vector<double> xTp, xTf, yTp, yTf; string sEx = baseEx + states[p] + ".txt"; + string sTs = baseTs + states[p] + ".txt"; string sTp = baseTp + states[p] + ".txt"; + string sTd = baseTd + states[p] + ".txt"; string sTf = baseTf + states[p] + ".txt"; - string sTm = baseTm + states[p] + ".txt"; cout << sEx << endl; @@ -591,6 +716,11 @@ void FigureTest_dt(){ gEx->SetMarkerStyle(7); gEx->SetTitle(""); + cout << sTs << endl; + TGraphErrors* gTs = GetExpDCS_ForFig(sTs); + gTs->SetLineColor(kMagenta); + gTs->SetTitle(""); + cout << sTp << endl; TGraphErrors* gTp = GetExpDCS_ForFig(sTp); gTp->SetLineColor(kRed); @@ -607,30 +737,35 @@ void FigureTest_dt(){ // gTp->GetYaxis()->SetTitle("d#sigma/d#Omega [mb/msr]"); // gTp->GetYaxis()->CenterTitle(true); + cout << sTd << endl; + TGraphErrors* gTd = GetExpDCS_ForFig(sTd); + gTd->SetLineColor(kGreen); + gTd->SetTitle(""); + cout << sTf << endl; TGraphErrors* gTf = GetExpDCS_ForFig(sTf); gTf->SetLineColor(kBlue); gTf->SetTitle(""); - TGraphErrors* gTm; - - gTp->SetLineWidth(2); - gTf->SetLineWidth(2); - if(l[p] == 1){ -// gTf->SetLineStyle(7); - } else if (l[p] == 3){ -// gTp->SetLineStyle(7); - } else if (l[p] == 2){ -// gTp->SetLineStyle(7); -// gTf->SetLineStyle(7); - - cout << sTm << endl; - gTm = GetExpDCS_ForFig(sTm); - gTm->SetLineColor(kViolet); - gTm->SetLineWidth(2); - gTm->SetTitle(""); - - } +// TGraphErrors* gTm; +// +// gTp->SetLineWidth(2); +// gTf->SetLineWidth(2); +// if(l[p] == 1){ +//// gTf->SetLineStyle(7); +// } else if (l[p] == 3){ +//// gTp->SetLineStyle(7); +// } else if (l[p] == 2){ +//// gTp->SetLineStyle(7); +//// gTf->SetLineStyle(7); +// +// cout << sTm << endl; +// gTm = GetExpDCS_ForFig(sTm); +// gTm->SetLineColor(kViolet); +// gTm->SetLineWidth(2); +// gTm->SetTitle(""); +// +// } @@ -640,34 +775,34 @@ void FigureTest_dt(){ //gTp->SetMaximum(1.2*h->GetMaximum()); // Format for y axis - gTp->GetYaxis()->SetLabelFont(43); - gTp->GetYaxis()->SetLabelSize(16); - gTp->GetYaxis()->SetLabelOffset(0.02); - gTp->GetYaxis()->SetTitle(""); - gTp->GetYaxis()->SetTitleFont(43); - gTp->GetYaxis()->SetTitleSize(25); - gTp->GetYaxis()->SetTitleOffset(6); - gTp->GetYaxis()->CenterTitle(); + gTs->GetYaxis()->SetLabelFont(43); + gTs->GetYaxis()->SetLabelSize(16); + gTs->GetYaxis()->SetLabelOffset(0.02); + gTs->GetYaxis()->SetTitle(""); + gTs->GetYaxis()->SetTitleFont(43); + gTs->GetYaxis()->SetTitleSize(25); + gTs->GetYaxis()->SetTitleOffset(6); + gTs->GetYaxis()->CenterTitle(); ////gTp->GetYaxis()->SetRangeUser(0.000006,0.015); ////gTp->GetYaxis()->SetRangeUser(0.00009,0.006); ////gTp->GetYaxis()->SetRangeUser(0.00002,0.0011); //gTp->GetYaxis()->SetRangeUser(0.0002,0.02); - gTp->GetYaxis()->SetTitle("d#sigma/d#Omega [mb/msr]"); + gTs->GetYaxis()->SetTitle("d#sigma/d#Omega [mb/msr]"); // TICKS Y Axis //gTp->GetYaxis()->SetTickLength(xFactor*0.04/yFactor); // Format for x axis - gTp->GetXaxis()->SetLabelFont(43); - gTp->GetXaxis()->SetLabelSize(16); - gTp->GetXaxis()->SetLabelOffset(0.02); - gTp->GetXaxis()->SetTitleFont(43); - gTp->GetXaxis()->SetTitleSize(25); - gTp->GetXaxis()->SetTitleOffset(5); - gTp->GetXaxis()->CenterTitle(); + gTs->GetXaxis()->SetLabelFont(43); + gTs->GetXaxis()->SetLabelSize(16); + gTs->GetXaxis()->SetLabelOffset(0.02); + gTs->GetXaxis()->SetTitleFont(43); + gTs->GetXaxis()->SetTitleSize(25); + gTs->GetXaxis()->SetTitleOffset(5); + gTs->GetXaxis()->CenterTitle(); //gTp->GetXaxis()->SetRangeUser(102.,158.); - gTp->GetXaxis()->SetRangeUser(3.,15.); - gTp->GetXaxis()->SetTitle("#theta_{CM} [deg]"); + gTs->GetXaxis()->SetRangeUser(3.,15.); + gTs->GetXaxis()->SetTitle("#theta_{CM} [deg]"); // TICKS X Axis //gTp->GetXaxis()->SetTickLength(yFactor*0.06/xFactor); @@ -686,111 +821,113 @@ void FigureTest_dt(){ - gTp->Draw("ac"); + gTs->Draw("ac"); + gTp->Draw("same c"); + gTd->Draw("same c"); gTf->Draw("same c"); - if(l[p]==2){gTm->Draw("same c");} + //if(l[p]==2){gTm->Draw("same c");} gEx->Draw("same P"); - - double e = (double)stoi(states[p])/1000.; - - - cout << e << endl; - - -// string titletext = to_string(e) + " MeV, S = "+-" -// + to_string(dS[p]); - - std::ostringstream out; - - TText text; - TText text2; - if(e==3.83){ - text.SetTextAlign(21); - text.SetTextFont(43); - text.SetTextSize(20); - text2.SetTextAlign(21); - text2.SetTextFont(43); - text2.SetTextSize(20); - - text.DrawTextNDC(XtoPad(0.5), YtoPad(0.88), "3.792 MeV, S = 0.18(1)"); - text2.DrawTextNDC(XtoPad(0.5), YtoPad(0.78), "3.868 MeV, S = 0.16(1)"); - - } else if (e==2.909){ - - text.SetTextAlign(21); - text.SetTextFont(43); - text.SetTextSize(20); - text2.SetTextAlign(21); - text2.SetTextFont(43); - text2.SetTextSize(20); - - text.DrawTextNDC(XtoPad(0.5), YtoPad(0.88), "2.908 MeV"); - text2.DrawTextNDC(XtoPad(0.5), YtoPad(0.78), "L=1 S=0.02(1), L=3 S=0.07(3)"); - - } else { - out.precision(3); - out << std::fixed << e << " MeV, S = " + S[p]; -// out.precision(2); -// out << S[p] << " +- " -// << dS[p]; - - text.SetTextAlign(21); - text.SetTextFont(43); - text.SetTextSize(20); - text.DrawTextNDC(XtoPad(0.5), YtoPad(0.88), (out.str()).c_str()); - - } - +// +// double e = (double)stoi(states[p])/1000.; +// // +// cout << e << endl; +// +// +//// string titletext = to_string(e) + " MeV, S = "+-" +//// + to_string(dS[p]); +// +// std::ostringstream out; // -// TH1F *hFrame = (TH1F*) h->Clone(TString::Format("h_%d_%d",i,j).Data()); -// -// // y axis range -// hFrame->SetMinimum(0.0001); // do not show 0 -// hFrame->SetMaximum(1.2*h->GetMaximum()); -// -// // Format for y axis -// hFrame->GetYaxis()->SetLabelFont(43); -// hFrame->GetYaxis()->SetLabelSize(16); -// hFrame->GetYaxis()->SetLabelOffset(0.02); -// hFrame->GetYaxis()->SetTitleFont(43); -// hFrame->GetYaxis()->SetTitleSize(16); -// hFrame->GetYaxis()->SetTitleOffset(2); -// -// hFrame->GetYaxis()->CenterTitle(); -// hFrame->GetYaxis()->SetNdivisions(505); -// -// // TICKS Y Axis -// hFrame->GetYaxis()->SetTickLength(xFactor*0.04/yFactor); -// -// // Format for x axis -// hFrame->GetXaxis()->SetLabelFont(43); -// hFrame->GetXaxis()->SetLabelSize(16); -// hFrame->GetXaxis()->SetLabelOffset(0.02); -// hFrame->GetXaxis()->SetTitleFont(43); -// hFrame->GetXaxis()->SetTitleSize(16); -// hFrame->GetXaxis()->SetTitleOffset(1); -// hFrame->GetXaxis()->CenterTitle(); -// hFrame->GetXaxis()->SetNdivisions(505); -// -// // TICKS X Axis -// hFrame->GetXaxis()->SetTickLength(yFactor*0.06/xFactor); -// -// // Draw cloned histogram with individual settings -// hFrame->Draw(); -// // TText text; -// text.SetTextAlign(31); +// TText text2; +// if(e==3.83){ +// text.SetTextAlign(21); // text.SetTextFont(43); -// text.SetTextSize(10); -// text.DrawTextNDC(XtoPad(0.9), YtoPad(0.8), gPad->GetName()); +// text.SetTextSize(20); +// text2.SetTextAlign(21); +// text2.SetTextFont(43); +// text2.SetTextSize(20); +// +// text.DrawTextNDC(XtoPad(0.5), YtoPad(0.88), "3.792 MeV, S = 0.18(1)"); +// text2.DrawTextNDC(XtoPad(0.5), YtoPad(0.78), "3.868 MeV, S = 0.16(1)"); +// +// } else if (e==2.909){ +// +// text.SetTextAlign(21); +// text.SetTextFont(43); +// text.SetTextSize(20); +// text2.SetTextAlign(21); +// text2.SetTextFont(43); +// text2.SetTextSize(20); +// +// text.DrawTextNDC(XtoPad(0.5), YtoPad(0.88), "2.908 MeV"); +// text2.DrawTextNDC(XtoPad(0.5), YtoPad(0.78), "L=1 S=0.02(1), L=3 S=0.07(3)"); +// +// } else { +// out.precision(3); +// out << std::fixed << e << " MeV, S = " + S[p]; +//// out.precision(2); +//// out << S[p] << " +- " +//// << dS[p]; +// +// text.SetTextAlign(21); +// text.SetTextFont(43); +// text.SetTextSize(20); +// text.DrawTextNDC(XtoPad(0.5), YtoPad(0.88), (out.str()).c_str()); +// +// } +// +//// +//// +//// TH1F *hFrame = (TH1F*) h->Clone(TString::Format("h_%d_%d",i,j).Data()); +//// +//// // y axis range +//// hFrame->SetMinimum(0.0001); // do not show 0 +//// hFrame->SetMaximum(1.2*h->GetMaximum()); +//// +//// // Format for y axis +//// hFrame->GetYaxis()->SetLabelFont(43); +//// hFrame->GetYaxis()->SetLabelSize(16); +//// hFrame->GetYaxis()->SetLabelOffset(0.02); +//// hFrame->GetYaxis()->SetTitleFont(43); +//// hFrame->GetYaxis()->SetTitleSize(16); +//// hFrame->GetYaxis()->SetTitleOffset(2); +//// +//// hFrame->GetYaxis()->CenterTitle(); +//// hFrame->GetYaxis()->SetNdivisions(505); +//// +//// // TICKS Y Axis +//// hFrame->GetYaxis()->SetTickLength(xFactor*0.04/yFactor); +//// +//// // Format for x axis +//// hFrame->GetXaxis()->SetLabelFont(43); +//// hFrame->GetXaxis()->SetLabelSize(16); +//// hFrame->GetXaxis()->SetLabelOffset(0.02); +//// hFrame->GetXaxis()->SetTitleFont(43); +//// hFrame->GetXaxis()->SetTitleSize(16); +//// hFrame->GetXaxis()->SetTitleOffset(1); +//// hFrame->GetXaxis()->CenterTitle(); +//// hFrame->GetXaxis()->SetNdivisions(505); +//// +//// // TICKS X Axis +//// hFrame->GetXaxis()->SetTickLength(yFactor*0.06/xFactor); +//// +//// // Draw cloned histogram with individual settings +//// hFrame->Draw(); +//// +//// TText text; +//// text.SetTextAlign(31); +//// text.SetTextFont(43); +//// text.SetTextSize(10); +//// text.DrawTextNDC(XtoPad(0.9), YtoPad(0.8), gPad->GetName()); } } - C->cd(); - -// C->SaveAs("FigureCS2s.pdf"); -// C->SaveAs("FigureCS2s.svg"); +// C->cd(); +// +//// C->SaveAs("FigureCS2s.pdf"); +//// C->SaveAs("FigureCS2s.svg"); } double XtoPad(double x){ @@ -873,7 +1010,8 @@ void Figure_AllCSFigures(){ gTp->GetXaxis()->SetRangeUser(102.,158.); gTp->GetYaxis()->SetRangeUser(0.00002,0.018); } else if (reactionName=="47K(d,t)") { - gTp->GetXaxis()->SetRangeUser(2.,12.); + //gTp->GetXaxis()->SetRangeUser(2.,12.); + gTp->GetXaxis()->SetRangeUser(2.,16.); gTp->GetYaxis()->SetRangeUser(0.002,0.18); } ////gTp->GetYaxis()->SetRangeUser(0.0002,0.008); @@ -1128,7 +1266,7 @@ void CS(double Energy, double Spin, double spdf, double angmom, double binning, // p3/2 -> spdf = 1, angmom = 1.5 // J0 is incident spin, which is 47K g.s. therefore J0 = 1/2 double J0 = 0.5; - double ElasticNorm = 0.000234, ElasticNormErr = 0.000; //18Oct22 + double ElasticNorm = 0.00023, ElasticNormErr = 0.00000; //18Oct22 inputdate = "18Oct22"; string backupFileName = "SolidAngle_HistFiles/"; @@ -1137,8 +1275,9 @@ void CS(double Energy, double Spin, double spdf, double angmom, double binning, if(reactionName=="47K(d,p)"){ // Which angular bin set to pull from - //numAngleBins=20; numAngleBinsBase=20; widthAngleBins=2.5; firstAngle=105.0; - numAngleBins=26; numAngleBinsBase=26; widthAngleBins=2.0; firstAngle=104.0; + numAngleBins=24; numAngleBinsBase=24; widthAngleBins=2.0; firstAngle=106.0; + //numAngleBins=26; numAngleBinsBase=26; widthAngleBins=2.0; firstAngle=104.0; + ////numAngleBins=15; numAngleBinsBase=15; widthAngleBins=2.0; firstAngle=126.0; // Selecting (d,p) global values numPeaks_CS = numPeaks; @@ -1153,7 +1292,8 @@ void CS(double Energy, double Spin, double spdf, double angmom, double binning, CSvertH = 5e-3; } else if(reactionName=="47K(d,t)"){ // Which angular bin set to pull from - numAngleBins=18; numAngleBinsBase=18; widthAngleBins=1.0; firstAngle=2.0; + //numAngleBins=18; numAngleBinsBase=18; widthAngleBins=1.0; firstAngle=2.0; + numAngleBins=8; numAngleBinsBase=8; widthAngleBins=2.0; firstAngle=2.0; // Selecting (d,p) global values numPeaks_CS = numPeaks_dt; @@ -1313,6 +1453,17 @@ void CS(double Energy, double Spin, double spdf, double angmom, double binning, * sqrt( pow(areaArray[i][2]/areaArray[i][1],2) + pow(SAerr/SA,2) ); + + //cout << BOLDRED << "!!!!\n!!!!!!!\n!!!!!!!!\n!!!!!!!!!" << endl; + //cout << "!! APPLYING ROUGH THETALAB EFFICIENCY 'CORRECTION' !! " << endl; + //double angle = ((areaArray[i][4]-areaArray[i][3])/2.)+areaArray[i][3]; + //tempvalue *= ((-2.27e-4 * angle * angle) + (6.98e-2 * angle) -4.60e0 ); + ////values here from fitting the trend of ratio clean/dirty against thetalab, + ////from the Test_RemoveBG_2D() funciton! + //cout << "!!!!!!!!!!\n!!!!!!!!\n!!!!!!!\n!!!!" << RESET <<endl; + + + //if(temperror<tempvalue){ // exclude very large error bars expDCS.push_back(tempvalue); expDCSerr.push_back(temperror); @@ -1374,6 +1525,13 @@ void CS(double Energy, double Spin, double spdf, double angmom, double binning, TCanvas* c_TWOFNR = new TCanvas("c_TWOFNR","c_TWOFNR",1000,1000); c_TWOFNR->SetLogy(); TGraph* TheoryDiffCross = TWOFNR(means_CS.at(indexE), J0, Spin, nodes, spdf, angmom); + cout << "TWOFNR(" + << means_CS.at(indexE) << " , " + << J0 << " , " + << 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(); @@ -1432,6 +1590,8 @@ void CS(double Energy, double Spin, double spdf, double angmom, double binning, } else { Final = FindNormalisation(TheoryDiffCross,gdSdO); } + + cout << "TEMP!!! INTEGR 4 to 14 = " << setprecision(8) << gdSdO2->Integral(4,14) << endl; Final->SetName("Final"); Final->SetTitle("Final"); @@ -1569,7 +1729,7 @@ vector<vector<double>> GetExpDiffCross(double Energy){ if(means_CS.at(indexE)>3.5){ //Find max angle, then round down to nearest integer double thetaMax = -14.36*means_CS.at(indexE) + 205.50; - thetaMax = floor(thetaMax-5.0);//TEMP REMOVE MORE!!!! + //thetaMax = floor(thetaMax-5.0);//TEMP REMOVE MORE!!!! if(loud){cout << RED << thetaMax << " rounds by " << widthAngleBins << " to ";} //Round down to nearest angular bin (i.e. 2.5deg sections) @@ -1589,16 +1749,39 @@ vector<vector<double>> GetExpDiffCross(double Energy){ //thetaMin = ceil(thetaMin); double thetaMin = 2; - if(Energy>3.0){thetaMin=3;} + //if(abs(Energy-3.3)<0.2){thetaMin=3;} + ////if(abs(Energy-1.7)<0.2){thetaMin=3;} + //if(Energy>4.0){thetaMin=4;} + if(abs(Energy-1.37)<0.2){thetaMin=4;} + if(abs(Energy-3.3)<0.2){thetaMin=4;} + //if(abs(Energy-1.7)<0.2){thetaMin=3;} if(Energy>4.0){thetaMin=4;} - double thetaMax = 11; - if(Energy<0.5){thetaMax=7;} - else if(Energy<2.0){thetaMax=8;} - else if(Energy<3.0){thetaMax=9;} - else if(Energy<4.0){thetaMax=10;} + //double thetaMax = 11; + //if(Energy<0.5){thetaMax=7;} + //else if(Energy<2.0){thetaMax=8;} + //else if(Energy<3.0){thetaMax=9;} + //else if(Energy<4.0){thetaMax=10;} + + //double thetaMax = 15; + //if(Energy<0.5){thetaMax=8;} + //else if(Energy<1.0){thetaMax=9;} + ////else if(Energy<1.8){thetaMax=5;} + //else if(Energy<2.0){thetaMax=12;} + //else if(Energy<3.0){thetaMax=13;} + //else if(Energy<4.0){thetaMax=14;} + //if(abs(Energy-1.7)<0.2){thetaMax=10;} + + double thetaMax = 16; + if(Energy<1.0){thetaMax=6;} + if(Energy<1.8){thetaMax=10;} + else if(Energy<2.8){thetaMax=12;} + else if(Energy<4.0){thetaMax=14;} + + - numAngleBins = (int)thetaMax - (int)thetaMin; + + numAngleBins = ((int)thetaMax - (int)thetaMin)/widthAngleBins; firstAngle = (int)thetaMin; cout << RED << "ANGLES:: " << firstAngle << " -> " << thetaMax << RESET << endl; @@ -1608,7 +1791,7 @@ vector<vector<double>> GetExpDiffCross(double Energy){ /* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */ - //for(int i=0; i<1; i++){ + //for(int i=0; i<5; i++){ for(int i=0; i<numAngleBins;i++){ double min = firstAngle + (i*widthAngleBins); double max = min + widthAngleBins; @@ -1668,7 +1851,13 @@ vector<vector<double>> GetExpDiffCross(double Energy){ // Retrieve array containing all fits, for one angle gate. // // Specific peak of interest selected from the vector by // // global variable indexE // - AllPeaks_OneGate = FitKnownPeaks_RtrnArry(gate, 0.0); + + if(chooseVAMOSoption == "original"){ + AllPeaks_OneGate = FitKnownPeaks_RtrnArry(gate, 0.0); + }else if(chooseVAMOSoption == "newgating"){ + //AllPeaks_OneGate = FitKnownPeaks_Pol2BG(gate); + AllPeaks_OneGate = FitKnownPeaks_SkewedGaussBG(gate); + } } else if (reactionName=="47K(d,t)"){ gate = PullThetaCMHist(i,firstAngle,widthAngleBins); @@ -1696,6 +1885,27 @@ vector<vector<double>> GetExpDiffCross(double Energy){ << endl; } + + /////// TESTING SCALE BY UNTIMEGATED RATIO /////// + //cout << "TESTING SCALE BY UNTIMEGATED RATIO" << endl; + //cout << "TESTING SCALE BY UNTIMEGATED RATIO" << endl; + //cout << "TESTING SCALE BY UNTIMEGATED RATIO" << endl; + //cout << "TESTING SCALE BY UNTIMEGATED RATIO" << endl; + //cout << "TESTING SCALE BY UNTIMEGATED RATIO" << endl; + //cout << "TESTING SCALE BY UNTIMEGATED RATIO" << endl; + //cout << "TESTING SCALE BY UNTIMEGATED RATIO" << endl; + //cout << "TESTING SCALE BY UNTIMEGATED RATIO" << endl; + //cout << "TESTING SCALE BY UNTIMEGATED RATIO" << endl; + //cout << "TESTING SCALE BY UNTIMEGATED RATIO" << endl; + //cout << "TESTING SCALE BY UNTIMEGATED RATIO" << endl; + //cout << "TESTING SCALE BY UNTIMEGATED RATIO" << endl; + //cout << "TESTING SCALE BY UNTIMEGATED RATIO" << endl; + //double ang = (min + max)/2.; + //AllPeaks_OneGate[indexE][1] = AllPeaks_OneGate[indexE][1]/(testGrad*ang+testCons); + //AllPeaks_OneGate[indexE][2] = AllPeaks_OneGate[indexE][2]/(testGrad*ang+testCons); + /////// TESTING SCALE BY UNTIMEGATED RATIO /////// + + /* Add min and max angle to end of relevant OneGate vector */ AllPeaks_OneGate[indexE].push_back(min); AllPeaks_OneGate[indexE].push_back(max); @@ -1715,7 +1925,18 @@ vector<vector<double>> GetExpDiffCross(double Energy){ //////////////////////////////////////////////////////////////////////////////// TH1F* PullThetaLabHist(int i, double minTheta, double gatesize){ //string name = "GateThetaLabHistograms_18Oct22_47Kdp.root"; - string name = "GateThetaLabHistograms_18Oct22_47Kdp_26angles.root"; + //string name = "GateThetaLabHistograms_11Aug23NewGating.root"; + + //string name = "GateThetaLabHistograms_18Oct22_47Kdp_26angles.root"; + //string name = "GateThetaLabHistograms_18Oct22_NewBroadTimeGates.root"; + + string name; + if(chooseVAMOSoption == "original"){ + name = "GateThetaLabHistograms_18Oct22_47Kdp_26angles.root"; + }else if(chooseVAMOSoption == "newgating"){ + name = "GateThetaLabHistograms_18Oct22_NewBroadTimeGates.root"; + } + TFile* file = new TFile(name.c_str(),"READ"); string histname = "cThetaLabGate_" + to_string((int)(minTheta+(i*gatesize))) + "-" @@ -1737,7 +1958,8 @@ TH1F* PullThetaLabHist(int i, double minTheta, double gatesize){ 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"; + //string name = "GateThetaCMHistograms_47Kdt_09Jun23"; + string name = "GateThetaCMHistograms_47Kdt_29Jan24"; if(!GammaGate){ //name = name + "bin0p1.root"; name = name + ".root"; @@ -1950,6 +2172,8 @@ TGraph* TWOFNR(double E, double J0, double J, double n, double l, double j){ system("pwd"); cout << "========================================================" << endl; + cout << "TEMP!!! INTEGR 4 to 14 = " << setprecision(8) << CS->Integral(4,14) << endl; + return CS; } @@ -2033,7 +2257,7 @@ TGraph* FindNormalisation(TGraph* theory, TGraphErrors* experiment){ int SpecFactorUpperLimit = 2; if(reactionName=="47K(d,t)"){ - SpecFactorUpperLimit = 12; + SpecFactorUpperLimit = 40; } // Set range of parameter(??) @@ -2168,3 +2392,4 @@ void CS(double Energy, double Spin, double spdf, double angmom, double binning, CS(Energy, Spin, spdf, angmom, binning, ""); } + diff --git a/Projects/e793s/macro/DrawPlots.h b/Projects/e793s/macro/DrawPlots.h index 995f67f9458dc94d53ab68bc19ef8876162a328f..7aeb8e22be75cc985e448ad68bd277f23348b4a4 100755 --- a/Projects/e793s/macro/DrawPlots.h +++ b/Projects/e793s/macro/DrawPlots.h @@ -2,10 +2,12 @@ #include "NPReaction.h" #include <string> #include <sstream> +#include "TGraphErrors.h" using namespace std; TChain* chain=NULL ; TChain* chainps=NULL ; +TChain* chainpsOld=NULL ; char cond[1000]; NPL::Reaction Cadp("47Ca(d,p)48Ca@355"); @@ -39,11 +41,18 @@ string exclBmDcy = "abs(AGATA_GammaE-2.013)>0.004 && abs(AGATA_GammaE-0.511)>0.0 /* BASE FUNCTIONS */ TF1* f_efficAGATA(){ - TF1 *f_E = new TF1("fit_1","TMath::Exp([0]+[1]*TMath::Log(x)+[2]*pow(TMath::Log(x),2.0)+[3]*pow(TMath::Log(x),3.0)+[4]*pow(TMath::Log(x),4.0))",10,6000); + TF1 *f_E = new TF1("f_E","TMath::Exp([0]+[1]*TMath::Log(x)+[2]*pow(TMath::Log(x),2.0)+[3]*pow(TMath::Log(x),3.0)+[4]*pow(TMath::Log(x),4.0))",10,6000); f_E->SetParameters(-6.34543e+01, +4.24746e+01, -1.00304e+01, +1.03468e+00, -3.97076e-02); return f_E; } +TF1* f_badEfficAGATA(){ + //before high energy correction + TF1 *f_Ebad = new TF1("f_Ebad","TMath::Exp([0]+[1]*TMath::Log(x)+[2]*pow(TMath::Log(x),2.0)+[3]*pow(TMath::Log(x),3.0)+[4]*pow(TMath::Log(x),4.0))",10,6000); + f_Ebad->SetParameters(-9.87216e+01, 6.98517e+01, -1.78362e+01, 2.00526e+00, -8.41808e-02); + return f_Ebad; +} + TChain* Chain(std::string TreeName, std::vector<std::string>& file, bool EventList){ TChain* chain = new TChain(TreeName.c_str()); unsigned int size =file.size(); @@ -55,12 +64,24 @@ TChain* Chain(std::string TreeName, std::vector<std::string>& file, bool EventLi } void LoadChainPhaseSpace(){ + vector<string> filesOld; vector<string> files; - for(int i=1; i<47; i++){ - string base = "../../../Outputs/Analysis/Sim_18Oct22_PhaseSpace_RandomSeed-"; - files.push_back((base + to_string(i) + ".root").c_str()); - } +// for(int i=1; i<47; i++){ +// //for(int i=1; i<5; i++){ +// string base = "../../../Outputs/Analysis/Sim_18Oct22_PhaseSpace_RandomSeed-"; +// filesOld.push_back((base + to_string(i) + ".root").c_str()); +// } +// +// chainpsOld = Chain("PhysicsTree",files,true); +// +// for(int i=1; i<25; i++){ +// string base = "../../../Outputs/Analysis/Sim_18Oct22_PhaseSpace_Thresh1p4_17Nov23-"; +// files.push_back((base + to_string(i) + ".root").c_str()); +// } +// + + files.push_back("../../../Outputs/Analysis/Sim_18Oct22_PhaseSpace_MegaNoThresh_01Dec23.root"); chainps = Chain("PhysicsTree",files,true); } @@ -100,8 +121,52 @@ void LoadChain47Kdp(){ /*************************/ - files.push_back("../../../Outputs/Analysis/47Kdp_18Oct22_10Jul23_PartI.root"); - files.push_back("../../../Outputs/Analysis/47Kdp_18Oct22_10Jul23_PartII.root"); + //files.push_back("../../../Outputs/Analysis/47Kdp_18Oct22_10Jul23_PartI.root"); + //files.push_back("../../../Outputs/Analysis/47Kdp_18Oct22_10Jul23_PartII.root"); + + //files.push_back("../../../Outputs/Analysis/47Kdp_25Jul23_PartI.root"); + //files.push_back("../../../Outputs/Analysis/47Kdp_25Jul23_PartII.root"); + + files.push_back("../../../Outputs/Analysis/47Kdp_08Aug23_18Oct22_TimeFromX_PartI.root"); + files.push_back("../../../Outputs/Analysis/47Kdp_08Aug23_18Oct22_TimeFromX_PartII.root"); + + /*************************/ + + //files.push_back("../../../Outputs/Analysis/47Kdp_UncorrectedXYZT_Run53.root"); + //files.push_back("../../../Outputs/Analysis/47Kdp_UncorrectedXYZT_Run54.root"); + //files.push_back("../../../Outputs/Analysis/47Kdp_UncorrectedXYZT_Run55.root"); + //files.push_back("../../../Outputs/Analysis/47Kdp_UncorrectedXYZT_Run56.root"); + //files.push_back("../../../Outputs/Analysis/47Kdp_UncorrectedXYZT_Run60.root"); + //files.push_back("../../../Outputs/Analysis/47Kdp_UncorrectedXYZT_Run61.root"); + //files.push_back("../../../Outputs/Analysis/47Kdp_UncorrectedXYZT_Run62.root"); + //files.push_back("../../../Outputs/Analysis/47Kdp_UncorrectedXYZT_Run63.root"); + //files.push_back("../../../Outputs/Analysis/47Kdp_UncorrectedXYZT_Run64.root"); + + + + //files.push_back("../../../Outputs/Analysis/47Kdp_10Aug23_18Oct22_IndividTimingMGX_PartI.root"); + //files.push_back("../../../Outputs/Analysis/47Kdp_10Aug23_18Oct22_IndividTimingMGX_PartII.root"); + + //files.push_back("../../../Outputs/Analysis/47Kdp_18Oct22_3p0um_PartI.root"); files.push_back("../../../Outputs/Analysis/47Kdp_18Oct22_3p0um_PartII.root"); + //files.push_back("../../../Outputs/Analysis/47Kdp_22Aug23_CalcVal_1p7um_PartI.root"); + //files.push_back("../../../Outputs/Analysis/47Kdp_25Jul23_5p5um_PartI.root"); files.push_back("../../../Outputs/Analysis/47Kdp_25Jul23_5p5um_PartII.root"); + //files.push_back("../../../Outputs/Analysis/47Kdp_22Aug23_CalcVal_1p7um_PartI.root"); files.push_back("../../../Outputs/Analysis/47Kdp_22Aug23_CalcVal_1p7um_PartII.root"); + //files.push_back("../../../Outputs/Analysis/47Kdp_18Oct23_HeavyXY_PartI.root"); files.push_back("../../../Outputs/Analysis/47Kdp_18Oct23_HeavyXY_PartII.root"); + + //files.push_back("../../../Outputs/Analysis/47Kdp_04Dec23_NewNoThreshold_PartI.root"); files.push_back("../../../Outputs/Analysis/47Kdp_04Dec23_NewNoThreshold_PartII.root"); + + + + //files.push_back("../../../Outputs/Analysis/47Kdp_08Jan24_OriginalValues_PartI.root"); files.push_back("../../../Outputs/Analysis/47Kdp_08Jan24_OriginalValues_PartII.root"); + + //files.push_back("../../../Outputs/Analysis/47Kdp_05Dec23_NoThreshold_AngleAGATA_PartI.root"); files.push_back("../../../Outputs/Analysis/47Kdp_05Dec23_NoThreshold_AngleAGATA_PartII.root"); + + + //////files.push_back("../../../Outputs/Analysis/47Kdp_18Oct23_10Jan24_HeavyTwoDownstreamPositions_PartI.root"); files.push_back("../../../Outputs/Analysis/47Kdp_18Oct23_10Jan24_HeavyTwoDownstreamPositions_PartII.root"); + + + + chain = Chain("PhysicsTree",files,true); } @@ -121,9 +186,14 @@ void LoadChain47Kdt(){ //files.push_back("../../../Outputs/Analysis/47Kdt_24May23_PostPhiFix_PartI.root"); //files.push_back("../../../Outputs/Analysis/47Kdt_24May23_PostPhiFix_PartII.root"); - files.push_back("../../../Outputs/Analysis/47Kdt_30May23_PhiFix_BetaMag_PartI.root"); - files.push_back("../../../Outputs/Analysis/47Kdt_30May23_PhiFix_BetaMag_PartII.root"); + //files.push_back("../../../Outputs/Analysis/47Kdt_30May23_PhiFix_BetaMag_PartI.root"); + //files.push_back("../../../Outputs/Analysis/47Kdt_30May23_PhiFix_BetaMag_PartII.root"); + + //files.push_back("../../../Outputs/Analysis/47Kdt_18Oct22_3p0um_PartI.root"); + //files.push_back("../../../Outputs/Analysis/47Kdt_18Oct22_3p0um_PartII.root"); + //files.push_back("../../../Outputs/Analysis/47Kdt_18Oct23_OldPositionMM_PartI.root"); files.push_back("../../../Outputs/Analysis/47Kdt_18Oct23_OldPositionMM_PartII.root"); + files.push_back("../../../Outputs/Analysis/47Kdt_18Oct23_HeavyXY_PartI.root"); files.push_back("../../../Outputs/Analysis/47Kdt_18Oct23_HeavyXY_PartII.root"); chain = Chain("PhysicsTree",files,true); } @@ -141,10 +211,19 @@ void LoadChain47Kdd(){ //files.push_back("../../../Outputs/Analysis/47Kdd_01Sep22_RemoveMoreMM5_PartII.root"); /* Now without runs 51&52, due to CATS not being in the trigger for these runs */ - files.push_back("../../../Outputs/Analysis/47Kdd_21Sep22_RmvMM5_NoRun51-52_PartI.root"); - files.push_back("../../../Outputs/Analysis/47Kdd_21Sep22_RmvMM5_NoRun51-52_PartII.root"); + //files.push_back("../../../Outputs/Analysis/47Kdd_21Sep22_RmvMM5_NoRun51-52_PartI.root"); + //files.push_back("../../../Outputs/Analysis/47Kdd_21Sep22_RmvMM5_NoRun51-52_PartII.root"); //cout << RED << " ONLY USING ONE PART OF THE SORT" << RESET << endl; +// files.push_back("../../../Outputs/Analysis/47Kdd_18Oct22_20Jul23_PartI.root"); +// files.push_back("../../../Outputs/Analysis/47Kdd_18Oct22_20Jul23_PartII.root"); + + + //files.push_back("../../../Outputs/Analysis/47Kdd_08Aug23_5umTarget_PartI.root"); + //files.push_back("../../../Outputs/Analysis/47Kdd_08Aug23_5umTarget_PartII.root"); + + files.push_back("../../../Outputs/Analysis/47Kdd_18Oct22_16Aug23_PartI.root"); + files.push_back("../../../Outputs/Analysis/47Kdd_18Oct22_16Aug23_PartII.root"); chain = Chain("PhysicsTree",files,true); } @@ -160,9 +239,14 @@ void LoadChain47Kpp(){ //files.push_back("../../../Outputs/Analysis/24Oct22_47Kpp_PartI.root"); //files.push_back("../../../Outputs/Analysis/24Oct22_47Kpp_PartII.root"); - files.push_back("../../../Outputs/Analysis/25Oct22_47Kpp_PartI.root"); - files.push_back("../../../Outputs/Analysis/25Oct22_47Kpp_PartII.root"); + //files.push_back("../../../Outputs/Analysis/25Oct22_47Kpp_PartI.root"); + //files.push_back("../../../Outputs/Analysis/25Oct22_47Kpp_PartII.root"); + //files.push_back("../../../Outputs/Analysis/47Kpp_18Oct22_20Jul23_PartI.root"); + //files.push_back("../../../Outputs/Analysis/47Kpp_18Oct22_20Jul23_PartII.root"); + + files.push_back("../../../Outputs/Analysis/47Kpp_18Oct22_16Aug23_PartI.root"); + files.push_back("../../../Outputs/Analysis/47Kpp_18Oct22_16Aug23_PartII.root"); chain = Chain("PhysicsTree",files,true); } @@ -354,8 +438,14 @@ void Load_1DGamma_MG(){ void Load_1DParticle(){ TCanvas *cEg = new TCanvas("cEg","cEg",1000,1000); TH1F *hEx = new TH1F("hEx","Loaded 1D Particle Spectrum",600,-15,15); - TFile *file = new TFile("LoadHistograms_48K/Load_1DParticle.root","READ"); - hEx = (TH1F*)file->Get("Ep"); + TFile *file; + if(reactionName=="47K(d,p)"){ + file = new TFile("LoadHistograms_48K/Load_1DParticle.root","READ"); + hEx = (TH1F*)file->Get("Ep"); + } else if(reactionName=="47K(d,t)"){ + file = new TFile("LoadHistograms_46K/Load_1DParticle.root","READ"); + hEx = (TH1F*)file->Get("hEx"); + } hEx->Draw(); } @@ -470,7 +560,11 @@ void Draw_1DParticle(){ void Load_2DParticleGamma(){ TCanvas *cExEg = new TCanvas("cExEg","cExEg",1000,1000); TH2F *hExEg = new TH2F("hExEg","Loaded 2D Particle-Gamma",600,-15,15,2500,0,5); - TFile *file = new TFile("LoadHistograms_15Feb23/Load_2DParticleGamma.root","READ"); + //TFile *file = new TFile("LoadHistograms_15Feb23/Load_2DParticleGamma.root","READ"); + //TFile *file = new TFile("LoadHistograms_48K/16Aug24_ExEg.root","READ"); + //TFile *file = new TFile("LoadHistograms_48K/16Aug24_ExEg_NoGradient.root","READ"); + //TFile *file = new TFile("LoadHistograms_48K/16Aug24_ExEg_AntiGradient.root","READ"); + TFile *file = new TFile("LoadHistograms_48K/16Aug24_ExEg_WithGradient.root","READ"); hExEg = (TH2F*)file->Get("ExEg"); hExEg->Draw("colz"); } @@ -486,10 +580,12 @@ void Draw_2DParticleGamma(){ TCanvas *cExEg = new TCanvas("cExEg","cExEg",1000,1000); gStyle->SetOptStat(0); //chain->Draw("AddBack_EDC:Ex>>ExEg(600,-15,15,2500,0,5)", gate.c_str(), "colz"); - chain->Draw("Ex:AddBack_EDC>>ExEg(10000,0,5,6000,-15,15)", gate.c_str(), ""); + chain->Draw("(Ex)/0.989:AddBack_EDC>>ExEg(10000,0,5,6000,-15,15)", gate.c_str(), ""); + //chain->Draw("Ex:AddBack_EDC>>ExEg(10000,0,5,6000,-15,15)", gate.c_str(), ""); + //chain->Draw("Ex:AddBack_EDC>>ExEg(10000,0,5,6000,-15,15)", gate.c_str(), ""); TH1F* ExEg = (TH1F*) gDirectory->Get("ExEg"); ExEg->SetTitle(""); - ExEg->GetYaxis()->SetTitle("Ex [MeV]"); + ExEg->GetYaxis()->SetTitle("E_{x} [MeV]"); ExEg->GetXaxis()->SetTitle("E_{#gamma} [MeV]"); ExEg->GetYaxis()->SetRangeUser(-1.0,7.0); ExEg->Draw(); @@ -1161,12 +1257,12 @@ void ForPoster_DiffCrossSec(){ void MugastMisalignment(){ TCanvas *cMisaligned = new TCanvas("cMisaligned","cMisaligned",1000,1000); gStyle->SetOptStat(0); - chain->Draw("Ex>>hcG_T1(120,-1,5)","abs(T_MUGAST_VAMOS-2700)<400 && Mugast.TelescopeNumber==1",""); - chain->Draw("Ex>>hcG_T2(120,-1,5)","abs(T_MUGAST_VAMOS-2700)<400 && Mugast.TelescopeNumber==2","same"); - chain->Draw("Ex>>hcG_T3(120,-1,5)","abs(T_MUGAST_VAMOS-2700)<400 && Mugast.TelescopeNumber==3","same"); - chain->Draw("Ex>>hcG_T4(120,-1,5)","abs(T_MUGAST_VAMOS-2700)<400 && Mugast.TelescopeNumber==4","same"); - chain->Draw("Ex>>hcG_T5(120,-1,5)","abs(T_MUGAST_VAMOS-2700)<400 && Mugast.TelescopeNumber==5","same"); - chain->Draw("Ex>>hcG_T7(120,-1,5)","abs(T_MUGAST_VAMOS-2700)<400 && Mugast.TelescopeNumber==7","same"); + chain->Draw("Ex>>hcG_T1(120,-1,5)","abs(T_MUGAST_VAMOS-270)<40 && Mugast.TelescopeNumber==1",""); + chain->Draw("Ex>>hcG_T2(120,-1,5)","abs(T_MUGAST_VAMOS-270)<40 && Mugast.TelescopeNumber==2","same"); + chain->Draw("Ex>>hcG_T3(120,-1,5)","abs(T_MUGAST_VAMOS-270)<40 && Mugast.TelescopeNumber==3","same"); + chain->Draw("Ex>>hcG_T4(120,-1,5)","abs(T_MUGAST_VAMOS-270)<40 && Mugast.TelescopeNumber==4","same"); + chain->Draw("Ex>>hcG_T5(120,-1,5)","abs(T_MUGAST_VAMOS-270)<40 && Mugast.TelescopeNumber==5","same"); + chain->Draw("Ex>>hcG_T7(120,-1,5)","abs(T_MUGAST_VAMOS-270)<40 && Mugast.TelescopeNumber==7","same"); TH1F* hcG_T1 = (TH1F*) gDirectory->Get("hcG_T1"); hcG_T1->SetTitle("Misalignment of MUGAST telescopes"); hcG_T1->GetXaxis()->SetTitle("E_{x} [MeV]"); @@ -1210,7 +1306,7 @@ void MugastMisalignment(){ legend->AddEntry(hcG_T7,"MUGAST 7","f"); legend->Draw(); - DrawParticleStates(cMisaligned); +// DrawParticleStates(cMisaligned); } void MugastMisalignment(double gamma, double width){ @@ -1282,7 +1378,7 @@ void ExPhiLab(){ TCanvas *diagnosePhi = new TCanvas("diagnosePhi","diagnosePhi",1000,1000); chain->Draw( "Ex:PhiLab>>phiHist(180,-180,180,120,-1,5)", - "abs(T_MUGAST_VAMOS-2700)<400 && Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8", + "abs(T_MUGAST_VAMOS-270)<40 && Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8", "colz"); TH1F* phiHist = (TH1F*) gDirectory->Get("phiHist"); phiHist->GetXaxis()->SetTitle("Phi (degrees)"); @@ -1965,6 +2061,581 @@ double expPyErr2[15]={0.5, canvThick->BuildLegend(); +} + +void FillVector_TWOFNR_Figure(string filename, + vector<double> &vecX, vector<double> &vecY){ + double x, y; + std::ifstream infile(filename.c_str()); + while (infile >> x >> y){ + if(y>0){ + vecX.push_back(x); + vecY.push_back(y); + } + } +} + +void ApplySettings_TWOFNR_Figure(TMultiGraph* mg){ + mg->GetXaxis()->SetTickLength(0.05); + mg->GetYaxis()->SetTickLength(0.05); + //mg->GetXaxis()->SetTitle("#theta_{lab}"); + mg->GetXaxis()->SetTitle("#theta_{CM}"); + mg->GetYaxis()->SetTitle("Differnetial cross section [mb/msr]"); + mg->GetXaxis()->SetTitleSize(0.06); + mg->GetYaxis()->SetTitleSize(0.06); + mg->GetXaxis()->SetLabelSize(0.05); + mg->GetYaxis()->SetLabelSize(0.05); + mg->GetXaxis()->CenterTitle(); + mg->GetYaxis()->CenterTitle(); + //mg->GetXaxis()->SetRangeUser(0.,180.); + //mg->GetXaxis()->SetRangeUser(91.,179.); + mg->GetXaxis()->SetRangeUser(1.,89.); + //mg->GetYaxis()->SetRangeUser(1e-5,1e-2); + //mg->GetYaxis()->SetRangeUser(1.5e-4,7e-3); + mg->GetYaxis()->SetRangeUser(2e-6,5e-2); +} + +void Figure_Transfer_dp(){ + + // THEORY ------------------------------------- + double blank, x, y; + string blankstr; + vector<double> vx, + vp12_sp0_e0, vp12_sp0_e2, vp12_sp0_e4, + vp12_sp1_e0, vp12_sp1_e2, vp12_sp1_e4, + vp32_sp1_e0, vp32_sp1_e2, vp32_sp1_e4, + vp32_sp2_e0, vp32_sp2_e2, vp32_sp2_e4, + vf52_sp2_e0, vf52_sp2_e2, vf52_sp2_e4, + vf52_sp3_e0, vf52_sp3_e2, vf52_sp3_e4, + vBlank; + + FillVector_TWOFNR_Figure("Transfer_dp_p12_sp0_e0.txt", vx, vp12_sp0_e0); + FillVector_TWOFNR_Figure("Transfer_dp_p12_sp1_e0.txt", vBlank, vp12_sp1_e0); + FillVector_TWOFNR_Figure("Transfer_dp_p32_sp1_e0.txt", vBlank, vp32_sp1_e0); + FillVector_TWOFNR_Figure("Transfer_dp_p32_sp2_e0.txt", vBlank, vp32_sp2_e0); + FillVector_TWOFNR_Figure("Transfer_dp_f52_sp2_e0.txt", vBlank, vf52_sp2_e0); + FillVector_TWOFNR_Figure("Transfer_dp_f52_sp3_e0.txt", vBlank, vf52_sp3_e0); + + FillVector_TWOFNR_Figure("Transfer_dp_p12_sp0_e2.txt", vBlank, vp12_sp0_e2); + FillVector_TWOFNR_Figure("Transfer_dp_p12_sp1_e2.txt", vBlank, vp12_sp1_e2); + FillVector_TWOFNR_Figure("Transfer_dp_p32_sp1_e2.txt", vBlank, vp32_sp1_e2); + FillVector_TWOFNR_Figure("Transfer_dp_p32_sp2_e2.txt", vBlank, vp32_sp2_e2); + FillVector_TWOFNR_Figure("Transfer_dp_f52_sp2_e2.txt", vBlank, vf52_sp2_e2); + FillVector_TWOFNR_Figure("Transfer_dp_f52_sp3_e2.txt", vBlank, vf52_sp3_e2); + + FillVector_TWOFNR_Figure("Transfer_dp_p12_sp0_e4.txt", vBlank, vp12_sp0_e4); + FillVector_TWOFNR_Figure("Transfer_dp_p12_sp1_e4.txt", vBlank, vp12_sp1_e4); + FillVector_TWOFNR_Figure("Transfer_dp_p32_sp1_e4.txt", vBlank, vp32_sp1_e4); + FillVector_TWOFNR_Figure("Transfer_dp_p32_sp2_e4.txt", vBlank, vp32_sp2_e4); + FillVector_TWOFNR_Figure("Transfer_dp_f52_sp2_e4.txt", vBlank, vf52_sp2_e4); + FillVector_TWOFNR_Figure("Transfer_dp_f52_sp3_e4.txt", vBlank, vf52_sp3_e4); + + TCanvas *canv = new TCanvas("canv","canv",600,700); + canv->Divide(2,3, 0.00001, 0.000001); + + TGraph* gp12_sp0_e0 = new TGraph(vx.size(), &vx[0], &vp12_sp0_e0[0]); + gp12_sp0_e0->SetLineColor(kOrange); + gp12_sp0_e0->SetLineStyle(1); + gp12_sp0_e0->SetLineWidth(2); + gp12_sp0_e0->SetTitle("p_{1/2}, spin = 0, E = 0 MeV"); + + TGraph* gp12_sp1_e0 = new TGraph(vx.size(), &vx[0], &vp12_sp1_e0[0]); + gp12_sp1_e0->SetLineColor(kOrange); + gp12_sp1_e0->SetLineStyle(1); + gp12_sp1_e0->SetLineWidth(2); + gp12_sp1_e0->SetTitle("p_{1/2}, spin = 1, E = 0 MeV"); + + TGraph* gp32_sp1_e0 = new TGraph(vx.size(), &vx[0], &vp32_sp1_e0[0]); + gp32_sp1_e0->SetLineColor(kRed); + gp32_sp1_e0->SetLineStyle(1); + gp32_sp1_e0->SetLineWidth(2); + gp32_sp1_e0->SetTitle("p_{3/2}, spin = 1, E = 0 MeV"); + + TGraph* gp32_sp2_e0 = new TGraph(vx.size(), &vx[0], &vp32_sp2_e0[0]); + gp32_sp2_e0->SetLineColor(kRed); + gp32_sp2_e0->SetLineStyle(1); + gp32_sp2_e0->SetLineWidth(2); + gp32_sp2_e0->SetTitle("p_{3/2}, spin = 2, E = 0 MeV"); + + TGraph* gf52_sp2_e0 = new TGraph(vx.size(), &vx[0], &vf52_sp2_e0[0]); + gf52_sp2_e0->SetLineColor(kBlue); + gf52_sp2_e0->SetLineStyle(1); + gf52_sp2_e0->SetLineWidth(2); + gf52_sp2_e0->SetTitle("f_{5/2}, spin = 2, E = 0 MeV"); + + TGraph* gf52_sp3_e0 = new TGraph(vx.size(), &vx[0], &vf52_sp3_e0[0]); + gf52_sp3_e0->SetLineColor(kBlue); + gf52_sp3_e0->SetLineStyle(1); + gf52_sp3_e0->SetLineWidth(2); + gf52_sp3_e0->SetTitle("f_{5/2}, spin = 3, E = 0 MeV"); + + TGraph* gp12_sp0_e2 = new TGraph(vx.size(), &vx[0], &vp12_sp0_e2[0]); + gp12_sp0_e2->SetLineColor(kOrange); + gp12_sp0_e2->SetLineStyle(9); + gp12_sp0_e2->SetLineWidth(2); + gp12_sp0_e2->SetTitle("p_{1/2}, spin = 0, E = 2 MeV"); + + TGraph* gp12_sp1_e2 = new TGraph(vx.size(), &vx[0], &vp12_sp1_e2[0]); + gp12_sp1_e2->SetLineColor(kOrange); + gp12_sp1_e2->SetLineStyle(9); + gp12_sp1_e2->SetLineWidth(2); + gp12_sp1_e2->SetTitle("p_{1/2}, spin = 1, E = 2 MeV"); + + TGraph* gp32_sp1_e2 = new TGraph(vx.size(), &vx[0], &vp32_sp1_e2[0]); + gp32_sp1_e2->SetLineColor(kRed); + gp32_sp1_e2->SetLineStyle(9); + gp32_sp1_e2->SetLineWidth(2); + gp32_sp1_e2->SetTitle("p_{3/2}, spin = 1, E = 2 MeV"); + + TGraph* gp32_sp2_e2 = new TGraph(vx.size(), &vx[0], &vp32_sp2_e2[0]); + gp32_sp2_e2->SetLineColor(kRed); + gp32_sp2_e2->SetLineStyle(9); + gp32_sp2_e2->SetLineWidth(2); + gp32_sp2_e2->SetTitle("p_{3/2}, spin = 2, E = 2 MeV"); + + TGraph* gf52_sp2_e2 = new TGraph(vx.size(), &vx[0], &vf52_sp2_e2[0]); + gf52_sp2_e2->SetLineColor(kBlue); + gf52_sp2_e2->SetLineStyle(9); + gf52_sp2_e2->SetLineWidth(2); + gf52_sp2_e2->SetTitle("f_{5/2}, spin = 2, E = 2 MeV"); + + TGraph* gf52_sp3_e2 = new TGraph(vx.size(), &vx[0], &vf52_sp3_e2[0]); + gf52_sp3_e2->SetLineColor(kBlue); + gf52_sp3_e2->SetLineStyle(9); + gf52_sp3_e2->SetLineWidth(2); + gf52_sp3_e2->SetTitle("f_{5/2}, spin = 3, E = 2 MeV"); + + TGraph* gp12_sp0_e4 = new TGraph(vx.size(), &vx[0], &vp12_sp0_e4[0]); + gp12_sp0_e4->SetLineColor(kOrange); + gp12_sp0_e4->SetLineStyle(2); + gp12_sp0_e4->SetLineWidth(2); + gp12_sp0_e4->SetTitle("p_{1/2}, spin = 0, E = 4 MeV"); + + TGraph* gp12_sp1_e4 = new TGraph(vx.size(), &vx[0], &vp12_sp1_e4[0]); + gp12_sp1_e4->SetLineColor(kOrange); + gp12_sp1_e4->SetLineStyle(2); + gp12_sp1_e4->SetLineWidth(2); + gp12_sp1_e4->SetTitle("p_{1/2}, spin = 1, E = 4 MeV"); + + TGraph* gp32_sp1_e4 = new TGraph(vx.size(), &vx[0], &vp32_sp1_e4[0]); + gp32_sp1_e4->SetLineColor(kRed); + gp32_sp1_e4->SetLineStyle(2); + gp32_sp1_e4->SetLineWidth(2); + gp32_sp1_e4->SetTitle("p_{3/2}, spin = 1, E = 4 MeV"); + + TGraph* gp32_sp2_e4 = new TGraph(vx.size(), &vx[0], &vp32_sp2_e4[0]); + gp32_sp2_e4->SetLineColor(kRed); + gp32_sp2_e4->SetLineStyle(2); + gp32_sp2_e4->SetLineWidth(2); + gp32_sp2_e4->SetTitle("p_{3/2}, spin = 2, E = 4 MeV"); + + TGraph* gf52_sp2_e4 = new TGraph(vx.size(), &vx[0], &vf52_sp2_e4[0]); + gf52_sp2_e4->SetLineColor(kBlue); + gf52_sp2_e4->SetLineStyle(2); + gf52_sp2_e4->SetLineWidth(2); + gf52_sp2_e4->SetTitle("f_{5/2}, spin = 2, E = 4 MeV"); + + TGraph* gf52_sp3_e4 = new TGraph(vx.size(), &vx[0], &vf52_sp3_e4[0]); + gf52_sp3_e4->SetLineColor(kBlue); + gf52_sp3_e4->SetLineStyle(2); + gf52_sp3_e4->SetLineWidth(2); + gf52_sp3_e4->SetTitle("f_{5/2}, spin = 3, E = 4 MeV"); + + canv->cd(1); + canv->cd(1)->SetLogy(); + canv->cd(1)->SetTopMargin( 0.00001); + canv->cd(1)->SetBottomMargin(0.15); + canv->cd(1)->SetLeftMargin( 0.25); + canv->cd(1)->SetRightMargin( 0.01); + TMultiGraph* mg_p12_sp0 = new TMultiGraph(); + mg_p12_sp0->Add(gp12_sp0_e0); + mg_p12_sp0->Add(gp12_sp0_e2); + mg_p12_sp0->Add(gp12_sp0_e4); + ApplySettings_TWOFNR_Figure(mg_p12_sp0); + mg_p12_sp0->GetYaxis()->SetRangeUser(1.5e-4,5e-3); + mg_p12_sp0->Draw("AC"); + + canv->cd(2); + canv->cd(2)->SetLogy(); + canv->cd(2)->SetTopMargin( 0.00001); + canv->cd(2)->SetBottomMargin(0.15); + canv->cd(2)->SetLeftMargin( 0.25); + canv->cd(2)->SetRightMargin( 0.01); + TMultiGraph* mg_p12_sp1 = new TMultiGraph(); + mg_p12_sp1->Add(gp12_sp1_e0); + mg_p12_sp1->Add(gp12_sp1_e2); + mg_p12_sp1->Add(gp12_sp1_e4); + ApplySettings_TWOFNR_Figure(mg_p12_sp1); + mg_p12_sp1->GetYaxis()->SetRangeUser(1.5e-4,5e-3); + mg_p12_sp1->Draw("AC"); + + canv->cd(3); + canv->cd(3)->SetLogy(); + canv->cd(3)->SetTopMargin( 0.00001); + canv->cd(3)->SetBottomMargin(0.15); + canv->cd(3)->SetLeftMargin( 0.25); + canv->cd(3)->SetRightMargin( 0.01); + TMultiGraph* mg_p32_sp1 = new TMultiGraph(); + mg_p32_sp1->Add(gp32_sp1_e0); + mg_p32_sp1->Add(gp32_sp1_e2); + mg_p32_sp1->Add(gp32_sp1_e4); + ApplySettings_TWOFNR_Figure(mg_p32_sp1); + mg_p32_sp1->GetYaxis()->SetRangeUser(5e-4,7e-3); + mg_p32_sp1->Draw("AC"); + + canv->cd(4); + canv->cd(4)->SetLogy(); + canv->cd(4)->SetTopMargin( 0.00001); + canv->cd(4)->SetBottomMargin(0.15); + canv->cd(4)->SetLeftMargin( 0.25); + canv->cd(4)->SetRightMargin( 0.01); + TMultiGraph* mg_p32_sp2 = new TMultiGraph(); + mg_p32_sp2->Add(gp32_sp2_e0); + mg_p32_sp2->Add(gp32_sp2_e2); + mg_p32_sp2->Add(gp32_sp2_e4); + ApplySettings_TWOFNR_Figure(mg_p32_sp2); + mg_p32_sp2->GetYaxis()->SetRangeUser(5e-4,7e-3); + mg_p32_sp2->Draw("AC"); + + canv->cd(5); + canv->cd(5)->SetLogy(); + canv->cd(5)->SetTopMargin( 0.00001); + canv->cd(5)->SetBottomMargin(0.15); + canv->cd(5)->SetLeftMargin( 0.25); + canv->cd(5)->SetRightMargin( 0.01); + TMultiGraph* mg_f52_sp2 = new TMultiGraph(); + mg_f52_sp2->Add(gf52_sp2_e0); + mg_f52_sp2->Add(gf52_sp2_e2); + mg_f52_sp2->Add(gf52_sp2_e4); + ApplySettings_TWOFNR_Figure(mg_f52_sp2); + mg_f52_sp2->GetYaxis()->SetRangeUser(1.5e-4,4e-3); + mg_f52_sp2->Draw("AC"); + + canv->cd(6); + canv->cd(6)->SetLogy(); + canv->cd(6)->SetTopMargin( 0.00001); + canv->cd(6)->SetBottomMargin(0.15); + canv->cd(6)->SetLeftMargin( 0.25); + canv->cd(6)->SetRightMargin( 0.01); + TMultiGraph* mg_f52_sp3 = new TMultiGraph(); + mg_f52_sp3->Add(gf52_sp3_e0); + mg_f52_sp3->Add(gf52_sp3_e2); + mg_f52_sp3->Add(gf52_sp3_e4); + ApplySettings_TWOFNR_Figure(mg_f52_sp3); + mg_f52_sp3->GetYaxis()->SetRangeUser(1.5e-4,4e-3); + mg_f52_sp3->Draw("AC"); + + + +} + +void Figure_Transfer_dt(){ + + // THEORY ------------------------------------- + double blank, x, y; + string blankstr; + vector<double> vx, + vs12_sp0_e0, vs12_sp0_e2, vs12_sp0_e4, + vs12_sp1_e0, vs12_sp1_e2, vs12_sp1_e4, + vp32_sp1_e0, vp32_sp1_e2, vp32_sp1_e4, + vp32_sp2_e0, vp32_sp2_e2, vp32_sp2_e4, + vd32_sp1_e0, vd32_sp1_e2, vd32_sp1_e4, + vd32_sp2_e0, vd32_sp2_e2, vd32_sp2_e4, + vf72_sp3_e0, vf72_sp3_e2, vf72_sp3_e4, + vf72_sp4_e0, vf72_sp4_e2, vf72_sp4_e4, + vBlank; + + FillVector_TWOFNR_Figure("Transfer_dt_s12_sp0_e0.txt", vx, vs12_sp0_e0); + FillVector_TWOFNR_Figure("Transfer_dt_s12_sp1_e0.txt", vBlank, vs12_sp1_e0); + FillVector_TWOFNR_Figure("Transfer_dt_p32_sp1_e0.txt", vBlank, vp32_sp1_e0); + FillVector_TWOFNR_Figure("Transfer_dt_p32_sp2_e0.txt", vBlank, vp32_sp2_e0); + FillVector_TWOFNR_Figure("Transfer_dt_d32_sp1_e0.txt", vBlank, vd32_sp1_e0); + FillVector_TWOFNR_Figure("Transfer_dt_d32_sp2_e0.txt", vBlank, vd32_sp2_e0); + FillVector_TWOFNR_Figure("Transfer_dt_f72_sp3_e0.txt", vBlank, vf72_sp3_e0); + FillVector_TWOFNR_Figure("Transfer_dt_f72_sp4_e0.txt", vBlank, vf72_sp4_e0); + + FillVector_TWOFNR_Figure("Transfer_dt_s12_sp0_e2.txt", vBlank, vs12_sp0_e2); + FillVector_TWOFNR_Figure("Transfer_dt_s12_sp1_e2.txt", vBlank, vs12_sp1_e2); + FillVector_TWOFNR_Figure("Transfer_dt_p32_sp1_e2.txt", vBlank, vp32_sp1_e2); + FillVector_TWOFNR_Figure("Transfer_dt_p32_sp2_e2.txt", vBlank, vp32_sp2_e2); + FillVector_TWOFNR_Figure("Transfer_dt_d32_sp1_e2.txt", vBlank, vd32_sp1_e2); + FillVector_TWOFNR_Figure("Transfer_dt_d32_sp2_e2.txt", vBlank, vd32_sp2_e2); + FillVector_TWOFNR_Figure("Transfer_dt_f72_sp3_e2.txt", vBlank, vf72_sp3_e2); + FillVector_TWOFNR_Figure("Transfer_dt_f72_sp4_e2.txt", vBlank, vf72_sp4_e2); + + FillVector_TWOFNR_Figure("Transfer_dt_s12_sp0_e4.txt", vBlank, vs12_sp0_e4); + FillVector_TWOFNR_Figure("Transfer_dt_s12_sp1_e4.txt", vBlank, vs12_sp1_e4); + FillVector_TWOFNR_Figure("Transfer_dt_p32_sp1_e4.txt", vBlank, vp32_sp1_e4); + FillVector_TWOFNR_Figure("Transfer_dt_p32_sp2_e4.txt", vBlank, vp32_sp2_e4); + FillVector_TWOFNR_Figure("Transfer_dt_d32_sp1_e4.txt", vBlank, vd32_sp1_e4); + FillVector_TWOFNR_Figure("Transfer_dt_d32_sp2_e4.txt", vBlank, vd32_sp2_e4); + FillVector_TWOFNR_Figure("Transfer_dt_f72_sp3_e4.txt", vBlank, vf72_sp3_e4); + FillVector_TWOFNR_Figure("Transfer_dt_f72_sp4_e4.txt", vBlank, vf72_sp4_e4); + + TCanvas *canv = new TCanvas("canv","canv",600,933); + canv->Divide(2,4, 0.00001, 0.000001); + + TGraph* gs12_sp0_e0 = new TGraph(vx.size(), &vx[0], &vs12_sp0_e0[0]); + gs12_sp0_e0->SetLineColor(kViolet); + gs12_sp0_e0->SetLineStyle(1); + gs12_sp0_e0->SetLineWidth(2); + gs12_sp0_e0->SetTitle("s_{1/2}, spin = 0, E = 0 MeV"); + + TGraph* gs12_sp1_e0 = new TGraph(vx.size(), &vx[0], &vs12_sp1_e0[0]); + gs12_sp1_e0->SetLineColor(kViolet); + gs12_sp1_e0->SetLineStyle(1); + gs12_sp1_e0->SetLineWidth(2); + gs12_sp1_e0->SetTitle("s_{1/2}, spin = 1, E = 0 MeV"); + + TGraph* gp32_sp1_e0 = new TGraph(vx.size(), &vx[0], &vp32_sp1_e0[0]); + gp32_sp1_e0->SetLineColor(kRed); + gp32_sp1_e0->SetLineStyle(1); + gp32_sp1_e0->SetLineWidth(2); + gp32_sp1_e0->SetTitle("p_{3/2}, spin = 1, E = 0 MeV"); + + TGraph* gp32_sp2_e0 = new TGraph(vx.size(), &vx[0], &vp32_sp2_e0[0]); + gp32_sp2_e0->SetLineColor(kRed); + gp32_sp2_e0->SetLineStyle(1); + gp32_sp2_e0->SetLineWidth(2); + gp32_sp2_e0->SetTitle("p_{3/2}, spin = 2, E = 0 MeV"); + + TGraph* gd32_sp1_e0 = new TGraph(vx.size(), &vx[0], &vd32_sp1_e0[0]); + gd32_sp1_e0->SetLineColor(kGreen); + gd32_sp1_e0->SetLineStyle(1); + gd32_sp1_e0->SetLineWidth(2); + gd32_sp1_e0->SetTitle("f_{5/2}, spin = 2, E = 0 MeV"); + + TGraph* gd32_sp2_e0 = new TGraph(vx.size(), &vx[0], &vd32_sp2_e0[0]); + gd32_sp2_e0->SetLineColor(kGreen); + gd32_sp2_e0->SetLineStyle(1); + gd32_sp2_e0->SetLineWidth(2); + gd32_sp2_e0->SetTitle("f_{5/2}, spin = 3, E = 0 MeV"); + + TGraph* gf72_sp3_e0 = new TGraph(vx.size(), &vx[0], &vf72_sp3_e0[0]); + gf72_sp3_e0->SetLineColor(kBlue); + gf72_sp3_e0->SetLineStyle(1); + gf72_sp3_e0->SetLineWidth(2); + gf72_sp3_e0->SetTitle("f_{5/2}, spin = 2, E = 0 MeV"); + + TGraph* gf72_sp4_e0 = new TGraph(vx.size(), &vx[0], &vf72_sp4_e0[0]); + gf72_sp4_e0->SetLineColor(kBlue); + gf72_sp4_e0->SetLineStyle(1); + gf72_sp4_e0->SetLineWidth(2); + gf72_sp4_e0->SetTitle("f_{5/2}, spin = 3, E = 0 MeV"); + + + TGraph* gs12_sp0_e2 = new TGraph(vx.size(), &vx[0], &vs12_sp0_e2[0]); + gs12_sp0_e2->SetLineColor(kViolet); + gs12_sp0_e2->SetLineStyle(9); + gs12_sp0_e2->SetLineWidth(2); + gs12_sp0_e2->SetTitle("s_{1/2}, spin = 0, E = 2 MeV"); + + TGraph* gs12_sp1_e2 = new TGraph(vx.size(), &vx[0], &vs12_sp1_e2[0]); + gs12_sp1_e2->SetLineColor(kViolet); + gs12_sp1_e2->SetLineStyle(9); + gs12_sp1_e2->SetLineWidth(2); + gs12_sp1_e2->SetTitle("s_{1/2}, spin = 1, E = 2 MeV"); + + TGraph* gp32_sp1_e2 = new TGraph(vx.size(), &vx[0], &vp32_sp1_e2[0]); + gp32_sp1_e2->SetLineColor(kRed); + gp32_sp1_e2->SetLineStyle(9); + gp32_sp1_e2->SetLineWidth(2); + gp32_sp1_e2->SetTitle("p_{3/2}, spin = 1, E = 2 MeV"); + + TGraph* gp32_sp2_e2 = new TGraph(vx.size(), &vx[0], &vp32_sp2_e2[0]); + gp32_sp2_e2->SetLineColor(kRed); + gp32_sp2_e2->SetLineStyle(9); + gp32_sp2_e2->SetLineWidth(2); + gp32_sp2_e2->SetTitle("p_{3/2}, spin = 2, E = 2 MeV"); + + TGraph* gd32_sp1_e2 = new TGraph(vx.size(), &vx[0], &vd32_sp1_e2[0]); + gd32_sp1_e2->SetLineColor(kGreen); + gd32_sp1_e2->SetLineStyle(9); + gd32_sp1_e2->SetLineWidth(2); + gd32_sp1_e2->SetTitle("f_{5/2}, spin = 2, E = 2 MeV"); + + TGraph* gd32_sp2_e2 = new TGraph(vx.size(), &vx[0], &vd32_sp2_e2[0]); + gd32_sp2_e2->SetLineColor(kGreen); + gd32_sp2_e2->SetLineStyle(9); + gd32_sp2_e2->SetLineWidth(2); + gd32_sp2_e2->SetTitle("f_{5/2}, spin = 3, E = 2 MeV"); + + TGraph* gf72_sp3_e2 = new TGraph(vx.size(), &vx[0], &vf72_sp3_e2[0]); + gf72_sp3_e2->SetLineColor(kBlue); + gf72_sp3_e2->SetLineStyle(9); + gf72_sp3_e2->SetLineWidth(2); + gf72_sp3_e2->SetTitle("f_{5/2}, spin = 2, E = 2 MeV"); + + TGraph* gf72_sp4_e2 = new TGraph(vx.size(), &vx[0], &vf72_sp4_e2[0]); + gf72_sp4_e2->SetLineColor(kBlue); + gf72_sp4_e2->SetLineStyle(9); + gf72_sp4_e2->SetLineWidth(2); + gf72_sp4_e2->SetTitle("f_{5/2}, spin = 3, E = 2 MeV"); + + TGraph* gs12_sp0_e4 = new TGraph(vx.size(), &vx[0], &vs12_sp0_e4[0]); + gs12_sp0_e4->SetLineColor(kViolet); + gs12_sp0_e4->SetLineStyle(2); + gs12_sp0_e4->SetLineWidth(2); + gs12_sp0_e4->SetTitle("s_{1/2}, spin = 0, E = 4 MeV"); + + TGraph* gs12_sp1_e4 = new TGraph(vx.size(), &vx[0], &vs12_sp1_e4[0]); + gs12_sp1_e4->SetLineColor(kViolet); + gs12_sp1_e4->SetLineStyle(2); + gs12_sp1_e4->SetLineWidth(2); + gs12_sp1_e4->SetTitle("s_{1/2}, spin = 1, E = 4 MeV"); + + TGraph* gp32_sp1_e4 = new TGraph(vx.size(), &vx[0], &vp32_sp1_e4[0]); + gp32_sp1_e4->SetLineColor(kRed); + gp32_sp1_e4->SetLineStyle(2); + gp32_sp1_e4->SetLineWidth(2); + gp32_sp1_e4->SetTitle("p_{3/2}, spin = 1, E = 4 MeV"); + + TGraph* gp32_sp2_e4 = new TGraph(vx.size(), &vx[0], &vp32_sp2_e4[0]); + gp32_sp2_e4->SetLineColor(kRed); + gp32_sp2_e4->SetLineStyle(2); + gp32_sp2_e4->SetLineWidth(2); + gp32_sp2_e4->SetTitle("p_{3/2}, spin = 2, E = 4 MeV"); + + TGraph* gd32_sp1_e4 = new TGraph(vx.size(), &vx[0], &vd32_sp1_e4[0]); + gd32_sp1_e4->SetLineColor(kGreen); + gd32_sp1_e4->SetLineStyle(2); + gd32_sp1_e4->SetLineWidth(2); + gd32_sp1_e4->SetTitle("f_{5/2}, spin = 2, E = 4 MeV"); + + TGraph* gd32_sp2_e4 = new TGraph(vx.size(), &vx[0], &vd32_sp2_e4[0]); + gd32_sp2_e4->SetLineColor(kGreen); + gd32_sp2_e4->SetLineStyle(2); + gd32_sp2_e4->SetLineWidth(2); + gd32_sp2_e4->SetTitle("f_{5/2}, spin = 3, E = 4 MeV"); + + TGraph* gf72_sp3_e4 = new TGraph(vx.size(), &vx[0], &vf72_sp3_e4[0]); + gf72_sp3_e4->SetLineColor(kBlue); + gf72_sp3_e4->SetLineStyle(2); + gf72_sp3_e4->SetLineWidth(2); + gf72_sp3_e4->SetTitle("f_{5/2}, spin = 2, E = 4 MeV"); + + TGraph* gf72_sp4_e4 = new TGraph(vx.size(), &vx[0], &vf72_sp4_e4[0]); + gf72_sp4_e4->SetLineColor(kBlue); + gf72_sp4_e4->SetLineStyle(2); + gf72_sp4_e4->SetLineWidth(2); + gf72_sp4_e4->SetTitle("f_{5/2}, spin = 3, E = 4 MeV"); + + canv->cd(1); + canv->cd(1)->SetLogy(); + canv->cd(1)->SetTopMargin( 0.00001); + canv->cd(1)->SetBottomMargin(0.15); + canv->cd(1)->SetLeftMargin( 0.25); + canv->cd(1)->SetRightMargin( 0.01); + TMultiGraph* mg_s12_sp0 = new TMultiGraph(); + mg_s12_sp0->Add(gs12_sp0_e0); + mg_s12_sp0->Add(gs12_sp0_e2); + mg_s12_sp0->Add(gs12_sp0_e4); + ApplySettings_TWOFNR_Figure(mg_s12_sp0); + mg_s12_sp0->GetYaxis()->SetRangeUser(8e-6,5e-2); + mg_s12_sp0->Draw("AC"); + +// canv->cd(2); +// canv->cd(2)->SetLogy(); +// canv->cd(2)->SetTopMargin( 0.00001); +// canv->cd(2)->SetBottomMargin(0.15); +// canv->cd(2)->SetLeftMargin( 0.25); +// canv->cd(2)->SetRightMargin( 0.01); +// TMultiGraph* mg_s12_sp1 = new TMultiGraph(); +// mg_s12_sp1->Add(gs12_sp1_e0); +// mg_s12_sp1->Add(gs12_sp1_e2); +// mg_s12_sp1->Add(gs12_sp1_e4); +// ApplySettings_TWOFNR_Figure(mg_s12_sp1); +// //mg_s12_sp1->GetYaxis()->SetRangeUser(1.5e-4,5e-3); +// mg_s12_sp1->Draw("AC"); + + //orig cd3 + canv->cd(2); + canv->cd(2)->SetLogy(); + canv->cd(2)->SetTopMargin( 0.00001); + canv->cd(2)->SetBottomMargin(0.15); + canv->cd(2)->SetLeftMargin( 0.25); + canv->cd(2)->SetRightMargin( 0.01); + TMultiGraph* mg_p32_sp1 = new TMultiGraph(); + mg_p32_sp1->Add(gp32_sp1_e0); + mg_p32_sp1->Add(gp32_sp1_e2); + mg_p32_sp1->Add(gp32_sp1_e4); + ApplySettings_TWOFNR_Figure(mg_p32_sp1); + mg_p32_sp1->GetYaxis()->SetRangeUser(8e-6,5e-2); + mg_p32_sp1->Draw("AC"); + +// canv->cd(4); +// canv->cd(4)->SetLogy(); +// canv->cd(4)->SetTopMargin( 0.00001); +// canv->cd(4)->SetBottomMargin(0.15); +// canv->cd(4)->SetLeftMargin( 0.25); +// canv->cd(4)->SetRightMargin( 0.01); +// TMultiGraph* mg_p32_sp2 = new TMultiGraph(); +// mg_p32_sp2->Add(gp32_sp2_e0); +// mg_p32_sp2->Add(gp32_sp2_e2); +// mg_p32_sp2->Add(gp32_sp2_e4); +// ApplySettings_TWOFNR_Figure(mg_p32_sp2); +// //mg_p32_sp2->GetYaxis()->SetRangeUser(5e-4,7e-3); +// mg_p32_sp2->Draw("AC"); + + //orig cd5 + canv->cd(3); + canv->cd(3)->SetLogy(); + canv->cd(3)->SetTopMargin( 0.00001); + canv->cd(3)->SetBottomMargin(0.15); + canv->cd(3)->SetLeftMargin( 0.25); + canv->cd(3)->SetRightMargin( 0.01); + TMultiGraph* mg_d32_sp1 = new TMultiGraph(); + mg_d32_sp1->Add(gd32_sp1_e0); + mg_d32_sp1->Add(gd32_sp1_e2); + mg_d32_sp1->Add(gd32_sp1_e4); + ApplySettings_TWOFNR_Figure(mg_d32_sp1); + mg_d32_sp1->GetYaxis()->SetRangeUser(2e-6,2e-3); + mg_d32_sp1->Draw("AC"); + +// canv->cd(6); +// canv->cd(6)->SetLogy(); +// canv->cd(6)->SetTopMargin( 0.00001); +// canv->cd(6)->SetBottomMargin(0.15); +// canv->cd(6)->SetLeftMargin( 0.25); +// canv->cd(6)->SetRightMargin( 0.01); +// TMultiGraph* mg_d32_sp2 = new TMultiGraph(); +// mg_d32_sp2->Add(gd32_sp2_e0); +// mg_d32_sp2->Add(gd32_sp2_e2); +// mg_d32_sp2->Add(gd32_sp2_e4); +// ApplySettings_TWOFNR_Figure(mg_d32_sp2); +// //mg_f72_sp4->GetYaxis()->SetRangeUser(1.5e-4,4e-3); +// mg_d32_sp2->Draw("AC"); + + //orig cd7 + canv->cd(4); + canv->cd(4)->SetLogy(); + canv->cd(4)->SetTopMargin( 0.00001); + canv->cd(4)->SetBottomMargin(0.15); + canv->cd(4)->SetLeftMargin( 0.25); + canv->cd(4)->SetRightMargin( 0.01); + TMultiGraph* mg_f72_sp3 = new TMultiGraph(); + mg_f72_sp3->Add(gf72_sp3_e0); + mg_f72_sp3->Add(gf72_sp3_e2); + mg_f72_sp3->Add(gf72_sp3_e4); + ApplySettings_TWOFNR_Figure(mg_f72_sp3); + mg_f72_sp3->GetYaxis()->SetRangeUser(2e-6,2e-3); + mg_f72_sp3->Draw("AC"); + +// canv->cd(8); +// canv->cd(8)->SetLogy(); +// canv->cd(8)->SetTopMargin( 0.00001); +// canv->cd(8)->SetBottomMargin(0.15); +// canv->cd(8)->SetLeftMargin( 0.25); +// canv->cd(8)->SetRightMargin( 0.01); +// TMultiGraph* mg_f72_sp4 = new TMultiGraph(); +// mg_f72_sp4->Add(gf72_sp4_e0); +// mg_f72_sp4->Add(gf72_sp4_e2); +// mg_f72_sp4->Add(gf72_sp4_e4); +// ApplySettings_TWOFNR_Figure(mg_f72_sp4); +// //mg_f72_sp4->GetYaxis()->SetRangeUser(1.5e-4,4e-3); +// mg_f72_sp4->Draw("AC"); + + } /* @@ -2165,13 +2836,24 @@ void AGATA_efficiency(double Energy_keV){ void ElasticsGate(double EMin, double EMax){ double width = EMax-EMin; double centre = EMin+(0.5*width); -// string gates = "abs(T_MUGAST_VAMOS-2700)<400 && MUST2.TelescopeNumber==5 && abs(ELab - " - string gate = timegate + "&&" + det_gate +// string gate = timegate + "&&" + det_gate +// + "&& abs(ELab - " +// + to_string(centre) +// + ")< " +// + to_string(width); + + //As of 21Jul23 - VAMOS timing investigation + //string gate = "T_MUGAST_VAMOS>10 && T_MUGAST_CATS2>115 &&" + string gate = " " + //string gate = "abs(MWT-35)>20 &&" + //string gate = "abs(T_MUGAST_CATS2-230)<200 && " + + det_gate + "&& abs(ELab - " + to_string(centre) + ")< " + to_string(width); + chain->Draw("ThetaLab>>hist(80,50,90)", gate.c_str(), ""); } @@ -2361,12 +3043,52 @@ void GateThetaLab_MultiWrite(double startTheta, double finishTheta, int numGates file->ls(); } +void GateThetaCM_MultiWrite(double startTheta, double finishTheta, int numGates, double binsize, double gateGammaE, double gateGammaWdth){ + string core = timegate + + " && " + det_gate + + " && abs(AddBack_EDC-" + to_string(gateGammaE) + ")<" + to_string(gateGammaWdth); + if(reactionName=="47K(d,t)"){ + core += " && cutTritons && cutTime"; + } + core += " && Ex@.size()==1 && ThetaLab > "; + string ytitle = "Counts / " + to_string(binsize) + " MeV"; + double gatesize = (finishTheta-startTheta)/numGates; + TList* list = new TList(); + + for (int i=0; i<numGates; i++){ + double minTheta = startTheta + (i * gatesize); + string title = to_string((int) minTheta)+" < ThetaLab < "+to_string((int) (minTheta+gatesize)); + string gating = core + + to_string(minTheta) + + " && ThetaCM < " + + to_string(minTheta+gatesize); + string histname = "cThetaCMGate_" + to_string((int) minTheta) + "-" + to_string((int) (minTheta+gatesize)); + string draw = "Ex>>" + histname + "(" + to_string(30.0/binsize) + ",-15,15)"; + + TCanvas *cEx_ThetaCMGate = new TCanvas(histname.c_str(),histname.c_str(),1000,1000); + chain->Draw(draw.c_str(),gating.c_str(),"colz"); + TH1F* Ex_ThetaCMGate = (TH1F*) gDirectory->Get(histname.c_str()); + Ex_ThetaCMGate->GetXaxis()->SetTitle("Ex [MeV]"); + Ex_ThetaCMGate->GetYaxis()->SetTitle(ytitle.c_str()); + Ex_ThetaCMGate->Sumw2(); + Ex_ThetaCMGate->SetTitle(title.c_str()); + list->Add(Ex_ThetaCMGate); + delete cEx_ThetaCMGate; + } + + TFile* file = new TFile("GateThetaCMHistograms_GammaGate.root","RECREATE"); + list->Write("GateThetaCMHistograms",TObject::kSingleKey); + file->ls(); +} + void GateThetaLab_MultiWrite(double startTheta, double finishTheta, int numGates, double binsize, double gateGammaE, double gateGammaWdth){ string core = timegate + " && " + det_gate - + " && abs(AddBack_EDC-" + to_string(gateGammaE) + ")<" + to_string(gateGammaWdth) - + " && Ex@.size()==1 && ThetaLab > "; -// string core = "abs(T_MUGAST_VAMOS-2700)<400 && Mugast.TelescopeNumber>0 && abs(AddBack_EDC-" + to_string(gateGammaE) + ")<" + to_string(gateGammaWdth) + " && Mugast.TelescopeNumber<8 && ThetaLab > "; + + " && abs(AddBack_EDC-" + to_string(gateGammaE) + ")<" + to_string(gateGammaWdth); + if(reactionName=="47K(d,t)"){ + core += " && cutTritons && cutTime"; + } + core += " && Ex@.size()==1 && ThetaLab > "; string ytitle = "Counts / " + to_string(binsize) + " MeV"; double gatesize = (finishTheta-startTheta)/numGates; TList* list = new TList(); @@ -3430,7 +4152,8 @@ void Figure_SolidAngle_dt(){ TH2F* Func_LoadIn2DPartGamma(){ string loadname = "./"; if(reactionName=="47K(d,p)"){ - loadname.append("LoadHistograms_48K/LoadMe_2DParticleGamma.root"); + //loadname.append("LoadHistograms_48K/LoadMe_2DParticleGamma.root"); + loadname.append("LoadHistograms_48K/LoadMe_15Feb24_2DParticleGamma.root"); } else if (reactionName=="47K(d,t)"){ loadname.append("LoadHistograms_46K/LoadMe_2DParticleGamma.root"); } @@ -3536,14 +4259,17 @@ void Func_AddGatingLines(double centre, double width, double height, EColor col) lHig->Draw("SAME"); } +//================================================================== + void GateGamma_SeeParticle_LoadFromFile(double gamma, double width, double gammaBin, double particleBin){ TH2F* h2D = (TH2F*)Func_LoadIn2DPartGamma(); gStyle->SetOptStat(0); h2D = Func_Rebin2DPartGamma(h2D, gammaBin, particleBin); - - TCanvas *cGGSP = new TCanvas("cGGSP","cGGSP",1000,1000); + stThesisWideStacked->cd(); + TCanvas *cGGSP = new TCanvas(); + //TCanvas *cGGSP = new TCanvas("cGGSP","cGGSP",1000,1000); cGGSP->Divide(1,2); double zoom = 0.2; @@ -3560,7 +4286,7 @@ void GateGamma_SeeParticle_LoadFromFile(double gamma, double width, double gamma cGGSP->cd(2); TH1F* hEx1 = Func_Gater(h2D, "Y", gamma, width, kBlack); - + //TH1F* hEx = (TH1F*)h2D->ProjectionY("hEx", // h2D->GetXaxis()->FindBin(gamma-width), // h2D->GetXaxis()->FindBin(gamma+width) @@ -3579,11 +4305,13 @@ void GateGamma_SeeParticle_LoadFromFile(double gamma1, double width1, double gam gStyle->SetPadLeftMargin(0.09); gStyle->SetOptStat(0); - TH1F* hExTotal = Return_1DParticle(); + //TH1F* hExTotal = Return_1DParticle(); TH2F* h2D = (TH2F*)Func_LoadIn2DPartGamma(); h2D = Func_Rebin2DPartGamma(h2D, gammaBin, particleBin); - TCanvas *cGGSP = new TCanvas("cGGSP","cGGSP",1000,600); + stThesisWideStacked->cd(); + TCanvas *cGGSP = new TCanvas(); + //TCanvas *cGGSP = new TCanvas("cGGSP","cGGSP",1000,600); cGGSP->Divide(1,2); double zoom = 0.2; @@ -3595,7 +4323,7 @@ void GateGamma_SeeParticle_LoadFromFile(double gamma1, double width1, double gam hEg->SetTitle(""); - hEg->GetXaxis()->SetTitle("Ex [MeV]"); + hEg->GetXaxis()->SetTitle("E_{#gamma} [MeV]"); hEg->GetXaxis()->CenterTitle(); hEg->GetXaxis()->SetTitleFont(133); hEg->GetXaxis()->SetTitleSize(22); @@ -3603,9 +4331,9 @@ void GateGamma_SeeParticle_LoadFromFile(double gamma1, double width1, double gam hEg->GetXaxis()->SetLabelSize(20); hEg->GetXaxis()->SetNdivisions(511); - string exytit = "Counts / " + to_string(particleBin) + " MeV"; + string egytit = "Counts / " + to_string(gammaBin) + " MeV"; - hEg->GetYaxis()->SetTitle(exytit.c_str()); + hEg->GetYaxis()->SetTitle(egytit.c_str()); hEg->GetYaxis()->CenterTitle(); hEg->GetYaxis()->SetTitleFont(133); hEg->GetYaxis()->SetTitleSize(22); @@ -3625,26 +4353,23 @@ void GateGamma_SeeParticle_LoadFromFile(double gamma1, double width1, double gam //hEx1->Draw(); //hEx2->Draw("same"); - string egytit = "Counts / " + to_string(particleBin) + " MeV"; - - ////// Un-comment for preseentable figure ////// - //hExTotal->SetTitle(""); - //hExTotal->GetXaxis()->CenterTitle(); - //hExTotal->GetXaxis()->SetTitleFont(133); - //hExTotal->GetXaxis()->SetTitleSize(22); - //hExTotal->GetXaxis()->SetLabelFont(133); - //hExTotal->GetXaxis()->SetLabelSize(20); - //hExTotal->GetXaxis()->SetNdivisions(511); - //hExTotal->GetYaxis()->CenterTitle(); - //hExTotal->GetYaxis()->SetTitleFont(133); - //hExTotal->GetYaxis()->SetTitleSize(22); - //hExTotal->GetYaxis()->SetLabelFont(133); - //hExTotal->GetYaxis()->SetLabelSize(20); - //hExTotal->GetYaxis()->SetTitleOffset(0.45); - //hExTotal->SetLineColor(kBlack); - //hExTotal->Scale(0.02); - //hExTotal->Draw(); - //////////////////////////////////////////////// + //// Un-comment for preseentable figure ////// + hEg->SetTitle(""); + hEg->GetXaxis()->CenterTitle(); + hEg->GetXaxis()->SetTitleFont(133); + hEg->GetXaxis()->SetTitleSize(30); + hEg->GetXaxis()->SetLabelFont(133); + hEg->GetXaxis()->SetLabelSize(30); + hEg->GetXaxis()->SetNdivisions(511); + hEg->GetYaxis()->CenterTitle(); + hEg->GetYaxis()->SetTitleFont(133); + hEg->GetYaxis()->SetTitleSize(30); + hEg->GetYaxis()->SetLabelFont(133); + hEg->GetYaxis()->SetLabelSize(30); + hEg->GetYaxis()->SetTitleOffset(0.45); + hEg->SetLineColor(kBlack); + hEg->Draw(); + ////////////////////////////////////////////// TH1F* hEx = (TH1F*)h2D->ProjectionY("hEx", h2D->GetXaxis()->FindBin(gamma1-width1), @@ -3652,8 +4377,48 @@ void GateGamma_SeeParticle_LoadFromFile(double gamma1, double width1, double gam ); hEx->GetXaxis()->SetRangeUser(-1,7); hEx->SetLineColor(kRed); + hEx->SetTitle(""); + hEx->GetXaxis()->SetTitle("E_{x} [MeV]"); + hEx->GetXaxis()->CenterTitle(); + hEx->GetXaxis()->SetTitleFont(133); + hEx->GetXaxis()->SetTitleSize(22); + hEx->GetXaxis()->SetLabelFont(133); + hEx->GetXaxis()->SetLabelSize(20); + hEx->GetXaxis()->SetNdivisions(511); + string exytit = "Counts / " + to_string(particleBin) + " MeV"; + hEx->GetYaxis()->SetTitle(exytit.c_str()); + hEx->GetYaxis()->CenterTitle(); + hEx->GetYaxis()->SetTitleFont(133); + hEx->GetYaxis()->SetTitleSize(22); + hEx->GetYaxis()->SetLabelFont(133); + hEx->GetYaxis()->SetLabelSize(20); + hEx->GetYaxis()->SetTitleOffset(0.45); + + + + //// Un-comment for preseentable figure ////// + hEx->SetTitle(""); + hEx->GetXaxis()->CenterTitle(); + hEx->GetXaxis()->SetTitleFont(133); + hEx->GetXaxis()->SetTitleSize(30); + hEx->GetXaxis()->SetLabelFont(133); + hEx->GetXaxis()->SetLabelSize(30); + hEx->GetXaxis()->SetNdivisions(511); + hEx->GetYaxis()->CenterTitle(); + hEx->GetYaxis()->SetTitleFont(133); + hEx->GetYaxis()->SetTitleSize(30); + hEx->GetYaxis()->SetLabelFont(133); + hEx->GetYaxis()->SetLabelSize(30); + hEx->GetYaxis()->SetTitleOffset(0.45); + hEx->Draw(); + ////////////////////////////////////////////// + + + + + + hEx->Draw(); - //hEx->Draw("same"); TH1F* hEx2 = (TH1F*)h2D->ProjectionY("hEx2", h2D->GetXaxis()->FindBin(gamma2-width2), @@ -3677,7 +4442,9 @@ void GateGamma_SeeParticle_LoadFromFile(double gamma, double width, double gamma gStyle->SetOptStat(0); h2D = Func_Rebin2DPartGamma(h2D, gammaBin, particleBin); - TCanvas *cGGSP = new TCanvas("cGGSP","cGGSP",1000,600); + stThesisWideStacked->cd(); + TCanvas *cGGSP = new TCanvas(); + //TCanvas *cGGSP = new TCanvas("cGGSP","cGGSP",1000,600); cGGSP->Divide(1,2); double zoom = 0.2; @@ -3784,7 +4551,9 @@ void GateParticle_SeeGamma_LoadFromFile(double gammaBin, double particleBin){ gStyle->SetOptStat(0); h2D = Func_Rebin2DPartGamma(h2D, gammaBin, particleBin); - TCanvas *cGPSG = new TCanvas("cGPSG","cGPSG",1000,1000); + stThesisWideStacked->cd(); + TCanvas *cGPSG = new TCanvas(); + //TCanvas *cGPSG = new TCanvas("cGPSG","cGPSG",1000,1000); cGPSG->Divide(1,2); cGPSG->cd(1); @@ -3831,7 +4600,10 @@ void GateParticle_SeeGamma_LoadFromFile(double particle, double width, double ga gStyle->SetOptStat(0); h2D = Func_Rebin2DPartGamma(h2D, gammaBin, particleBin); - TCanvas *cGPSG = new TCanvas("cGPSG","cGPSG",1000,1000); + stThesisWideStacked->cd(); + TCanvas *cGPSG = new TCanvas(); + //TCanvas *cGPSG = new TCanvas("cGPSG","cGPSG",1000,1000); + //TCanvas *cGPSG = new TCanvas("cGPSG","cGPSG",1000,600); cGPSG->Divide(1,2); double zoom = 0.2; @@ -3862,7 +4634,9 @@ void GateParticle_SeeGamma_LoadFromFile(double particle, double width, double pa gStyle->SetOptStat(0); h2D = Func_Rebin2DPartGamma(h2D, gammaBin, particleBin); - TCanvas *cGPSG = new TCanvas("cGPSG","cGPSG",1000,1000); + stThesisWideStacked->cd(); + TCanvas *cGPSG = new TCanvas(); + //TCanvas *cGPSG = new TCanvas("cGPSG","cGPSG",1000,600); cGPSG->Divide(1,2); cGPSG->cd(1); @@ -3898,7 +4672,9 @@ void GateParticle_SeeGamma_LoadFromFile(double particle, double width, double pa gStyle->SetOptStat(0); h2D = Func_Rebin2DPartGamma(h2D, gammaBin, particleBin); - TCanvas *cGPSG = new TCanvas("cGPSG","cGPSG",1000,1000); + stThesisWideStacked->cd(); + TCanvas *cGPSG = new TCanvas(); + //TCanvas *cGPSG = new TCanvas("cGPSG","cGPSG",1000,1000); cGPSG->Divide(1,2); cGPSG->cd(1); @@ -3934,6 +4710,8 @@ void GateParticle_SeeGamma_LoadFromFile(double particle, double width, double pa hEg3->Draw("same"); } +//================================================================== + void GateGamma_SeeGamma_LoadFromFile(double gamma, double width, double gammaBin1, double gammaBin2){ TH2F* h2D = (TH2F*)Func_LoadIn2DGammaGamma(); @@ -4277,10 +5055,10 @@ void Fig_Ex_Eg_SepCanvases_SameScale(){ } //chain->Draw("Ex>>Ep(22,-0.5,5)", gate.c_str(),""); - chain->Draw("Ex>>Ep(220,-0.5,5)", gate.c_str(),""); + chain->Draw("Ex>>Ep(160,-1,7)", gate.c_str(),""); TH1F* Ep = (TH1F*) gDirectory->Get("Ep"); Ep->GetXaxis()->SetTitle("Ex [MeV]"); - Ep->GetYaxis()->SetTitle("Counts / 0.025 MeV"); + Ep->GetYaxis()->SetTitle("Counts / 0.05 MeV"); Ep->SetTitle(""); Ep->GetXaxis()->CenterTitle(); Ep->GetXaxis()->SetTitleFont(133); @@ -4341,7 +5119,7 @@ void Fig_Ex_Eg_SepCanvases_SameScale(){ cEg->SetLogy(); gate = gate + " && " + exclBmDcy; - chain->Draw("AddBack_EDC>>Eg(1100,-0.5,5)",gate.c_str(),""); + chain->Draw("AddBack_EDC>>Eg(1600,-1,7)",gate.c_str(),""); //chain->Draw("AddBack_EDC>>Eg(110,-0.5,5)",gate.c_str(),""); TH1F* Eg = (TH1F*) gDirectory->Get("Eg"); Eg->GetXaxis()->SetTitle("E_{#gamma} [MeV]"); @@ -4389,55 +5167,6 @@ void Fig_Ex_Eg_SepCanvases_SameScale(){ EgT->Draw("same"); -} - -void Testing_BetaCorrection_dt(){ - -TCanvas* canv = new TCanvas("canv","canv",1000,1000); - -chain->Draw("AGATA_OrigBetaMag:ThetaLab>>hAll3(60,0,30,250,0.125,0.127)", - "abs(T_MUGAST_VAMOS-280)<30 && MUST2.TelescopeNumber<5 && cutTritons && cutTime && (abs(Ex-2.0)<0.05 || abs(Ex-3.4)<0.05 || abs(Ex-0.6)<0.05)", - //"abs(Ex-2.0)<0.05 || abs(Ex-3.4)<0.05 || abs(Ex-0.6)<0.05", - ""); - -TH1F* hAll3 = (TH1F*) gDirectory->Get("hAll3"); -hAll3->Draw(); - -auto t0p6 = TMarker(5.0,0.12658,22); -t0p6.SetMarkerColor(kRed); -t0p6.SetMarkerSize(3); -t0p6.SetMarkerStyle(22); - -t0p6.DrawMarker(5.0 ,12.6576); -t0p6.DrawMarker(10.0,12.6567); -t0p6.DrawMarker(15.0,12.6551); -t0p6.DrawMarker(20.0,12.6526); - - -auto t2p0 = TMarker(5.0,0.12658,22); -t2p0.SetMarkerColor(kOrange); -t2p0.SetMarkerSize(3); -t2p0.SetMarkerStyle(23); - -t2p0.DrawMarker(5.0 ,12.6224); -t2p0.DrawMarker(10.0,12.6211); -t2p0.DrawMarker(15.0,12.6187); -t2p0.DrawMarker(20.0,12.6150); - - - -auto t3p4 = TMarker(5.0,0.12658,22); -t3p4.SetMarkerColor(kViolet+3); -t3p4.SetMarkerSize(3); -t3p4.SetMarkerStyle(20); - -t3p4.DrawMarker(5.0 ,12.5846); -t3p4.DrawMarker(10.0,12.5828); -t3p4.DrawMarker(15.0,12.5793); -t3p4.DrawMarker(20.0,12.5738); - - - } double sq(double val){ @@ -4541,8 +5270,7 @@ void Func_SetXYAxisSize(TGraph* hist, int labSize, int titSize, double labOff, d hist->GetYaxis()->CenterTitle(); } -void Func_RotateFigure(TH1F* hist){ -//TGraph* Func_RotateFigure(TH1F* hist){ +pair<TGraph*,TString> Func_RotateFigure(TH1F* hist){ /////////////////// ROTATING Ex FIGURE ///////////////////// // From root-forum.cern.ch/t/rotation-of-1d-histogram/758 // Double_t x, y, dx2, dy2, err; @@ -4552,110 +5280,120 @@ void Func_RotateFigure(TH1F* hist){ Bool_t logZscale = kFALSE; Bool_t logXscale = kFALSE; Bool_t logYscale = kFALSE; - TString option_2dim("COLZ"); - TString option_1dim("HIST Y+"); // options for the 1-dim hists + TString option_1dim("HIST Y+"); // options for the 1-dim hists TString option_1dim_y(option_1dim.Data()); // y-projection is a TGraph TString fill_option(""); - Int_t fill_style = 0; - Int_t fill_color = 0; - - Int_t nbins = hist->GetNbinsX(); - Double_t ymin = hist->GetXaxis()->GetXmin(); // nota bene; x <-> y - Double_t ymax = hist->GetXaxis()->GetXmax(); - Double_t xmin = hist->GetMinimum(); - Double_t xmax = hist->GetMaximum(); - Double_t logxmin = xmax; - xmax += 0.1 * xmax; - - Double_t * xv = new Double_t[nbins]; - Double_t * yv = new Double_t[nbins]; - Double_t * xe = new Double_t[nbins]; - Double_t * ye = new Double_t[nbins]; - Double_t logymin = ymax; - - for(Int_t i=1; i <= nbins; i++){ - y = hist->GetBinLowEdge(i); - dy2 = 0.5 * hist->GetBinWidth(i); - x = hist->GetBinContent(i); - if(x > 0 && x < logxmin) logxmin = x; - err = hist->GetBinError(i); - - xv[i-1] = x; - yv[i-1] = y+dy2; - xe[i-1] = err; - ye[i-1] = dy2; - } - if(logZscale)xmin = logxmin; - else if (xmin>0) xmin =0; - - TH2F * hy = new TH2F("hy", "", 10, xmin, xmax, 10, ymin, ymax); - hy->Draw(); - hy->GetXaxis()->SetNdivisions(505); - if(option_1dim_y.Contains("E", TString::kIgnoreCase)){ - if(!option_1dim_y.Contains("1", TString::kIgnoreCase)){ - option_1dim_y += "z"; // switch off the little lines at error bars - } - ge = new TGraphErrors(nbins,xv,yv, xe, ye); - ge->Draw(option_1dim_y.Data()); - } else { - if(option_1dim_y.Contains("H", TString::kIgnoreCase)){ - yv[0] -= dy2; - yv[1] = yv[nbins-1] + dy2; - option_1dim_y += "R"; // rotate - g = new TGraph(nbins, xv, yv); - if(fill_option.Contains("F", TString::kIgnoreCase)){ - option_1dim_y += "F"; - g->SetFillStyle(fill_style); - g->SetFillColor(fill_color); - } - } else { - g = new TGraph(nbins, xv, yv); - cout <<xv[nbins-1] << " " << yv[nbins-1] << endl; - } + Int_t fill_style = 4000; + + Int_t trans; + new TColor(trans, 0.5, 0.5, 0.5, "transparent", 0.0); + + Int_t nbins = hist->GetNbinsX(); + Double_t ymin = hist->GetXaxis()->GetXmin(); // nota bene; x <-> y + Double_t ymax = hist->GetXaxis()->GetXmax(); + Double_t xmin = hist->GetMinimum(); + Double_t xmax = hist->GetMaximum(); + Double_t logxmin = xmax; + xmax += 0.1 * xmax; + + Double_t * xv = new Double_t[nbins]; + Double_t * yv = new Double_t[nbins]; + Double_t * xe = new Double_t[nbins]; + Double_t * ye = new Double_t[nbins]; + Double_t logymin = ymax; + + for(Int_t i=1; i <= nbins; i++){ + y = hist->GetBinLowEdge(i); + dy2 = 0.5 * hist->GetBinWidth(i); + x = hist->GetBinContent(i); + if(x > 0 && x < logxmin) logxmin = x; + err = hist->GetBinError(i); + + xv[i-1] = x; + yv[i-1] = y+dy2; + xe[i-1] = err; + ye[i-1] = dy2; + } - int labSize = 20; - int titSize = 20; - double labOff = 0.005; - double titOff = 1.0; - - g->GetXaxis()->SetLabelFont(133); - g->GetXaxis()->SetLabelSize(labSize); //in pixels - g->GetYaxis()->SetLabelFont(133); - g->GetYaxis()->SetLabelSize(labSize); //in pixels - g->GetXaxis()->SetLabelOffset(labOff); - g->GetYaxis()->SetLabelOffset(labOff); - - g->GetXaxis()->SetTitleFont(133); - g->GetXaxis()->SetTitleSize(titSize); //in pixels - g->GetYaxis()->SetTitleFont(133); - g->GetYaxis()->SetTitleSize(titSize); //in pixels - g->GetXaxis()->SetTitleOffset(titOff); - g->GetYaxis()->SetTitleOffset(titOff); - - g->GetXaxis()->CenterTitle(); - g->GetYaxis()->CenterTitle(); - - g->Draw(option_1dim_y.Data()); - - } - delete [] xv; - delete [] yv; - if(xe) {delete [] xe; xe = 0;}; - if(ye) {delete [] ye; ye = 0;}; - - //return g; + if(logZscale) xmin = logxmin; + else if (xmin>0) xmin = 0; + + TH2F * hy = new TH2F("hy", "", 10, xmin, xmax, 10, ymin, ymax); + hy->Draw(); + hy->GetXaxis()->SetNdivisions(505); + if(option_1dim_y.Contains("E", TString::kIgnoreCase)){ + if(!option_1dim_y.Contains("1", TString::kIgnoreCase)){ + option_1dim_y += "z"; // switch off the little lines at error bars + } + ge = new TGraphErrors(nbins,xv,yv, xe, ye); + ge->Draw(option_1dim_y.Data()); + } else { + if(option_1dim_y.Contains("H", TString::kIgnoreCase)){ + yv[0] -= dy2; + yv[1] = yv[nbins-1] + dy2; + option_1dim_y += "R"; // rotate + g = new TGraph(nbins, xv, yv); + if(fill_option.Contains("F", TString::kIgnoreCase)){ + option_1dim_y += "F"; + g->SetFillStyle(fill_style); + g->SetFillColor(trans); + } + } else { + g = new TGraph(nbins, xv, yv); + cout <<xv[nbins-1] << " " << yv[nbins-1] << endl; + } + + int labSize = 20; + int titSize = 20; + double labOff = 0.001; + double titOff = 0.3; + + g->GetXaxis()->SetLabelFont(133); + g->GetYaxis()->SetLabelFont(133); + g->GetXaxis()->SetLabelSize(labSize); //in pixels + g->GetYaxis()->SetLabelSize(labSize); //in pixels + g->GetXaxis()->SetLabelOffset(labOff); + g->GetYaxis()->SetLabelOffset(labOff); + + g->GetXaxis()->SetTitle("Counts"); + g->GetYaxis()->SetTitle("Ex []"); + g->GetXaxis()->CenterTitle(); + g->GetYaxis()->CenterTitle(); + g->GetXaxis()->SetTitleFont(133); + g->GetYaxis()->SetTitleFont(133); + g->GetXaxis()->SetTitleSize(titSize); //in pixels + g->GetYaxis()->SetTitleSize(titSize); //in pixels + g->GetXaxis()->SetTitleOffset(titOff); + g->GetYaxis()->SetTitleOffset(titOff); + + g->GetXaxis()->SetRangeUser(0,500); + + //g->Draw(option_1dim_y.Data()); + + } + delete [] xv; + delete [] yv; + if(xe) {delete [] xe; xe = 0;}; + if(ye) {delete [] ye; ye = 0;}; + + + pair<TGraph*,TString> output; + output.first = g; + output.second = option_1dim_y; + + //return g; + return output; //////////////////////////////////////////////////////////// } void Figure_ExEg(){ - int xwidth = 1000; + int xwidth = 1000;//1000; int ywidth = 1000; int border = 5; int single = 300; - TH2F* h2D = (TH2F*)Func_LoadIn2DPartGamma(); @@ -4665,7 +5403,7 @@ void Figure_ExEg(){ TCanvas *cExEg = new TCanvas("cExEg","cExEg",xwidth,ywidth); gStyle->SetOptStat(0); - gStyle->SetCanvasColor(0); + gStyle->SetCanvasColor(4000); //////////////////////////////// // ---------- ExEg ---------- // @@ -4682,21 +5420,22 @@ void Figure_ExEg(){ pExEg->SetBottomMargin(0.11); pExEg->SetLeftMargin( 0.11); pExEg->SetRightMargin( 0.00001); - pExEg->SetFillStyle(0); + pExEg->SetFillStyle(4000); + pExEg->SetFillColor(kWhite); pExEg->cd(); h2D->GetXaxis()->SetRangeUser(0.0,4.75); h2D->GetYaxis()->SetRangeUser(-1.0,7.5); Func_SetXYAxisSize(h2D, 20, 22, 0.005, 1.2); -// h2D->GetXaxis()->SetLabelFont(133); -// h2D->GetXaxis()->SetLabelSize(16); //in pixels -// h2D->GetYaxis()->SetLabelFont(133); -// h2D->GetYaxis()->SetLabelSize(16); //in pixels -// h2D->GetXaxis()->SetTitleFont(133); -// h2D->GetXaxis()->SetTitleSize(20); //in pixels -// h2D->GetXaxis()->CenterTitle(); -// h2D->GetYaxis()->SetTitleFont(133); -// h2D->GetYaxis()->SetTitleSize(20); //in pixels -// h2D->GetYaxis()->CenterTitle(); + h2D->GetXaxis()->SetLabelFont(133); + h2D->GetXaxis()->SetLabelSize(20); //in pixels + h2D->GetYaxis()->SetLabelFont(133); + h2D->GetYaxis()->SetLabelSize(20); //in pixels + h2D->GetXaxis()->SetTitleFont(133); + h2D->GetXaxis()->SetTitleSize(25); //in pixels + h2D->GetXaxis()->CenterTitle(); + h2D->GetYaxis()->SetTitleFont(133); + h2D->GetYaxis()->SetTitleSize(55); //in pixels + h2D->GetYaxis()->CenterTitle(); h2D->Draw(); @@ -4715,11 +5454,13 @@ void Figure_ExEg(){ pEg->SetBottomMargin(0.00001); pEg->SetLeftMargin( 0.11); pEg->SetRightMargin( 0.00001); - pEg->SetFillStyle(0); + pEg->SetFillStyle(4000); + pEg->SetFillColor(kWhite); pEg->SetBorderMode(0); pEg->SetLogy(); pEg->SetTickx(); pEg->SetTicky(); + pEg->Draw(); pEg->cd(); TH1F* hEg = (TH1F*)h2D->ProjectionX("hEg",0,10000); @@ -4743,49 +5484,521 @@ void Figure_ExEg(){ pEx->SetBottomMargin(0.11); pEx->SetLeftMargin( 0.00001); pEx->SetRightMargin( 0.11); - pEx->SetFillStyle(0); + pEx->SetFillStyle(4000); + pEx->SetFillColor(kWhite); pEx->SetBorderMode(0); pEx->SetTickx(); pEg->SetTicky(); + pEx->Draw(); pEx->cd(); TH1F* hEx = (TH1F*)h2D->ProjectionY("hEx",0,10000); hEx->GetXaxis()->SetRangeUser(-1.0,7.5); Func_SetXYAxisSize(hEx, 20, 22, 0.005, 1.2); + hEx->GetYaxis()->SetTitle("TESTME!"); hEx->Draw(); - Func_RotateFigure(hEx); - //auto gEx = Func_RotateFigure(hEx); + auto pairEx = Func_RotateFigure(hEx); + pairEx.first->Draw(pairEx.second.Data()); - //TH1F* hExG1 = (TH1F*)h2D->ProjectionY("hExG1", - // h2D->GetXaxis()->FindBin(0.140), - // h2D->GetXaxis()->FindBin(0.145) - // ); - //hExG1->SetLineColor(kRed); + double gamma = 1.838; + double width = 0.005; +// +// cExEg->cd(); +// TPad *pExR = new TPad("pExR","pExR",0.3,0.3,0.9,0.9); +// pExR->SetPad("pExR","pExR", +// cExEg->AbsPixeltoX(xwidth - border - single), +// cExEg->AbsPixeltoY(single + border), +// cExEg->AbsPixeltoX(xwidth - border), +// cExEg->AbsPixeltoY(ywidth - border) +// ); +// pExR->Draw(); +// pExR->SetTopMargin( 0.00001); +// pExR->SetBottomMargin(0.11); +// pExR->SetLeftMargin( 0.00001); +// pExR->SetRightMargin( 0.11); +// pExR->SetFillStyle(4000); +// pExR->SetFillColor(kWhite); +// pExR->SetBorderMode(0); +// pExR->SetTickx(); +// pExR->SetTicky(); +// pExR->Draw(); +// pExR->cd(); +// +// +// TH1F* hExRed = (TH1F*)h2D->ProjectionY("hExRed",h2D->GetXaxis()->FindBin(gamma-width),h2D->GetXaxis()->FindBin(gamma+width)); +// hExRed->SetLineColor(kRed); +// hExRed->Draw("hist"); +// +// auto pairExR = Func_RotateFigure(hExRed); +// pairExR.first->Draw(pairExR.second.Data()); - //auto gExG1 = Func_RotateFigure(hExG1); - //gEx->Draw("hist"); -// TMultiGraph* mg = new TMultiGraph(); -// mg->Add(gEx, "hist"); -// mg->Add(gExG1, "hist"); -// mg->Draw("a"); +} +void Figure_ExEg_Vers2(){ + + int border = 5; + int single = 270; + int doubel = 500; + int xwidth = doubel + doubel + 3*border; + int ywidth = doubel + single + 2*border; + double margbig = 0.21; + double margsml = 0.00001; + + int textfont = 133; + int titlesize = 28; double titleoff1 = 1.9; double titleoff2 = titleoff1;//+0.4; + int labelsize = 28;// double labeloffs = 0.005; + + + + //Int_t col1; col1 = TColor::GetColor("#20e0eb"); + //Int_t col2; col2 = TColor::GetColor("#535acf"); + //Int_t col3; col3 = TColor::GetColor("#e100e6"); + + Int_t col1; col1 = TColor::GetColor("#882E72"); + Int_t col2; col2 = TColor::GetColor("#4EB265"); + Int_t col3; col3 = TColor::GetColor("#E8601C"); + Int_t col4; col4 = TColor::GetColor("#DC050C"); + + + + + TH2F* hExEg = (TH2F*)Func_LoadIn2DPartGamma(); + + hExEg->RebinX(10); + hExEg->RebinY(16); + hExEg->GetXaxis()->SetRangeUser(-0.00,4.90); + hExEg->GetYaxis()->SetRangeUser(-0.30,4.64);//7.50); + + TCanvas *cExEg = new TCanvas("cExEg","cExEg",xwidth,ywidth); + gStyle->SetOptStat(0); + gStyle->SetCanvasColor(4000); + + //////////////////////////////// + // ---------- ExEg ---------- // + cExEg->cd(); + TPad *pExEg = new TPad("pExEg","pExEg",0.1,0.1,0.6,0.6); + pExEg->SetPad("pExEg","pExEg", + cExEg->AbsPixeltoX(border ), + cExEg->AbsPixeltoY(border + single ), + cExEg->AbsPixeltoX(border + doubel ), + cExEg->AbsPixeltoY(border + single + doubel) + ); + pExEg->SetTopMargin( margsml); + pExEg->SetBottomMargin(margbig); + pExEg->SetLeftMargin( margbig); + pExEg->SetRightMargin( margsml); + pExEg->SetFillStyle(4000); + pExEg->SetFillColor(kWhite); + pExEg->SetBorderMode(0); + pExEg->SetTickx(); + pExEg->SetTicky(); + pExEg->Draw(); + pExEg->cd(); + hExEg->GetXaxis()->SetNdivisions(508); + hExEg->GetXaxis()->SetTitleFont(textfont); + hExEg->GetXaxis()->SetTitleSize(titlesize); + hExEg->GetXaxis()->SetTitleOffset(titleoff1); + hExEg->GetXaxis()->SetLabelFont(textfont); +// hExEg->GetXaxis()->SetLabelOffset(labeloffs); + hExEg->GetXaxis()->SetLabelSize(labelsize); + hExEg->GetXaxis()->CenterTitle(); + hExEg->GetXaxis()->SetDecimals(1); + hExEg->GetYaxis()->SetNdivisions(508); + hExEg->GetYaxis()->SetTitleFont(textfont); + hExEg->GetYaxis()->SetTitleSize(titlesize); + hExEg->GetYaxis()->SetTitleOffset(titleoff1); + hExEg->GetYaxis()->SetLabelFont(textfont); + hExEg->GetYaxis()->SetLabelSize(labelsize); +// hExEg->GetYaxis()->SetLabelOffset(labeloffs); + hExEg->GetYaxis()->CenterTitle(); + hExEg->GetYaxis()->SetDecimals(1); + hExEg->Draw(); + + TLine *lExEg = new TLine(0.000,0.000,4.644,4.644); + lExEg->SetLineColor(col4); + lExEg->SetLineStyle(2); + lExEg->SetLineWidth(2); + lExEg->Draw(); +// TLine *Sn = new TLine(0.000,4.644,4.75,4.644); +// Sn->SetLineColor(kViolet); +// Sn->SetLineStyle(7); +// Sn->Draw(); + + + //////////////////////////////// + // ----------- Eg ----------- // + cExEg->cd(); + TPad *p__Eg = new TPad("pEg","pEg",0.3,0.3,0.9,0.9); + p__Eg->SetPad("pEg","pEg", + cExEg->AbsPixeltoX(border ), + cExEg->AbsPixeltoY(border ), + cExEg->AbsPixeltoX(border + doubel ), + cExEg->AbsPixeltoY(border + single ) + ); + p__Eg->SetTopMargin( margbig); + p__Eg->SetBottomMargin(margsml); + p__Eg->SetLeftMargin( margbig); + p__Eg->SetRightMargin( margsml); + p__Eg->SetFillStyle(4000); + p__Eg->SetFillColor(kWhite); + p__Eg->SetBorderMode(0); + p__Eg->SetLogy(); + p__Eg->SetTickx(); + p__Eg->SetTicky(); + p__Eg->Draw(); + p__Eg->cd(); + + TH1F* h__Eg = (TH1F*)hExEg->ProjectionX("h__Eg",0,10000); + h__Eg->GetXaxis()->SetTitleFont(textfont); + h__Eg->GetXaxis()->SetTitleSize(titlesize); + h__Eg->GetXaxis()->SetTitleOffset(titleoff2); + h__Eg->GetXaxis()->SetLabelFont(textfont); + h__Eg->GetXaxis()->SetLabelSize(labelsize); +// h__Eg->GetXaxis()->SetLabelOffset(labeloffs); + h__Eg->GetXaxis()->CenterTitle(); + h__Eg->GetYaxis()->SetTitleFont(textfont); + h__Eg->GetYaxis()->SetTitleSize(titlesize); + h__Eg->GetYaxis()->SetTitleOffset(titleoff2); + h__Eg->GetYaxis()->SetLabelFont(textfont); + h__Eg->GetYaxis()->SetLabelSize(labelsize); +// h__Eg->GetYaxis()->SetLabelOffset(labeloffs); + h__Eg->GetYaxis()->CenterTitle(); + h__Eg->GetYaxis()->SetRangeUser(2e-1,7e3); + + h__Eg->GetXaxis()->SetTickSize(0.03*(doubel/single)); + + h__Eg->SetLineColor(kBlack); + h__Eg->GetYaxis()->SetTitle("Counts / 0.005 MeV"); + h__Eg->Draw(); + + //////////////////////////////// + // ----------- Ex ----------- // + cExEg->cd(); + TPad *pEx__ = new TPad("pEx","pEx",0.1,0.1,0.7,0.7); + pEx__->SetPad("pEx","pEx", + cExEg->AbsPixeltoX(border + doubel + border), + cExEg->AbsPixeltoY(border ), + cExEg->AbsPixeltoX(border + doubel + border + doubel), + cExEg->AbsPixeltoY(border + single ) + ); + pEx__->SetTopMargin( margbig); + pEx__->SetBottomMargin(margsml); + pEx__->SetLeftMargin( margbig); + pEx__->SetRightMargin( margsml); + pEx__->SetFillStyle(4000); + pEx__->SetFillColor(kWhite); + pEx__->SetBorderMode(0); + pEx__->SetTickx(); + pEx__->SetTicky(); + pEx__->Draw(); + pEx__->cd(); + + TH1F* hEx__ = (TH1F*)hExEg->ProjectionY("hEx__",0,10000); + hEx__->GetXaxis()->SetTitleFont(textfont); + hEx__->GetXaxis()->SetTitleSize(titlesize); + hEx__->GetXaxis()->SetTitleOffset(titleoff2); + hEx__->GetXaxis()->SetLabelFont(textfont); + hEx__->GetXaxis()->SetLabelSize(labelsize); +// hEx__->GetXaxis()->SetLabelOffset(labeloffs); + hEx__->GetXaxis()->CenterTitle(); + hEx__->GetYaxis()->SetTitleFont(textfont); + hEx__->GetYaxis()->SetTitleSize(titlesize); + hEx__->GetYaxis()->SetTitleOffset(titleoff2); + hEx__->GetYaxis()->SetLabelFont(textfont); + hEx__->GetYaxis()->SetLabelSize(labelsize); +// hEx__->GetYaxis()->SetLabelOffset(labeloffs); + hEx__->GetYaxis()->CenterTitle(); + hEx__->GetYaxis()->SetNdivisions(506); + hEx__->GetYaxis()->SetRangeUser(1,900.); + + hEx__->GetXaxis()->SetTickSize(0.03*(doubel/single)); + + hEx__->SetLineColor(kBlack); + //hEx__->GetYaxis()->SetRangeUser(5.,550.); + //hEx__->GetYaxis()->SetTitle("#splitline{Counts}{/0.0XX MeV}"); + hEx__->GetYaxis()->SetTitle("Counts / 0.08 MeV"); + hEx__->SetFillColorAlpha(kBlack,0.15); + hEx__->SetFillStyle(1001); + + // ----------- Ex ----------- // + double widthA = 0.005; + double gammaR = 1.838; + TH1F* hEx18 = (TH1F*)hExEg->ProjectionY( + "hEx18", + hExEg->GetXaxis()->FindBin(gammaR-widthA), + hExEg->GetXaxis()->FindBin(gammaR+widthA) + ); + hEx18->SetLineColor(col1); + //hEx18->SetLineWidth(0); + hEx18->SetFillColorAlpha(col1,0.6); + hEx18->SetFillStyle(1001); + + // ----------- Ex ----------- // + double gammaB = 2.872; + TH1F* hEx36 = (TH1F*)hExEg->ProjectionY( + "hEx36", + hExEg->GetXaxis()->FindBin(gammaB-widthA), + hExEg->GetXaxis()->FindBin(gammaB+widthA) + ); + hEx36->SetLineColor(col2); + //hEx36->SetLineWidth(0); + hEx36->SetFillColorAlpha(col2,0.6); + hEx36->SetFillStyle(1001); + + // ----------- Ex ----------- // + double gammaG = 3.515; + TH1F* hEx38 = (TH1F*)hExEg->ProjectionY( + "hEx38", + hExEg->GetXaxis()->FindBin(gammaG-widthA), + hExEg->GetXaxis()->FindBin(gammaG+widthA) + ); + hEx38->SetLineColor(col3); + //hEx38->SetLineWidth(0); + hEx38->SetFillColorAlpha(col3,0.6); + hEx38->SetFillStyle(1001); + +// hEx18->Rebin(2); hEx36->Rebin(2); hEx38->Rebin(2); + +// hEx18->SetFillColor(col1); +// hEx36->SetFillColor(col2); +// hEx38->SetFillColor(col3); + + + + hEx__->Draw("hist RX"); + hEx18->Scale( 7.14); hEx18->Draw("hist same"); + hEx36->Scale(40.39); hEx36->Draw("hist same"); + hEx38->Scale(12.05); hEx38->Draw("hist same"); + //hEx18->Scale( 8.); hEx18->Draw("hist same"); + //hEx36->Scale(30.); hEx36->Draw("hist same"); + //hEx38->Scale(26.); hEx38->Draw("hist same"); + + pExEg->cd(); + TLine *lineR = new TLine(gammaR, -0.35, gammaR, 7.5); + lineR->SetLineColorAlpha(col1,0.4); + lineR->SetLineWidth(5); + lineR->Draw(); + TLine *lineB = new TLine(gammaB, -0.35, gammaB, 7.5); + lineB->SetLineColorAlpha(col2,0.4); + lineB->SetLineWidth(5); + lineB->Draw(); + TLine *lineG = new TLine(gammaG, -0.35, gammaG, 7.5); + lineG->SetLineColorAlpha(col3,0.4); + lineG->SetLineWidth(5); + lineG->Draw(); + TLatex *Texeg = new TLatex(.5,.5,"Ex = E_{#gamma}"); + Texeg->SetTextColor(kRed); + Texeg->SetTextSize(0.05); + Texeg->SetX(4.05); + Texeg->SetY(3.60); + Texeg->SetTextAngle(44); + Texeg->Draw("same"); + + + p__Eg->cd(); + TLine *liner = new TLine(gammaR, 1e-9, gammaR, 7e3); + liner->SetLineColorAlpha(col1,0.4); + liner->SetLineWidth(5); + liner->Draw(); + TLine *lineb = new TLine(gammaB, 1e-9, gammaB, 7e3); + lineb->SetLineColorAlpha(col2,0.4); + lineb->SetLineWidth(5); + lineb->Draw(); + TLine *lineg = new TLine(gammaG, 1e-9, gammaG, 7e3); + lineg->SetLineColorAlpha(col3,0.4); + lineg->SetLineWidth(5); + lineg->Draw(); + +} + +void Testing_BGremovalELTL(){ + //chain->Draw("ELab:ThetaLab>>hELTL(26,104,156,200,1,11)","!cutBG_MGTvE && Mugast.DSSD_T<770 && abs(T_MUGAST_CATS2-230)<200","colz"); + chain->Draw("Ex:ThetaLab>>hELTL(26,104,156,500,-15,10)","!cutBG_MGTvE && Mugast.DSSD_T<770 && abs(T_MUGAST_CATS2-230)<200","colz"); + TH2F* hELTL = (TH2F*) gDirectory->Get("hELTL"); + +// plot_kine(Kdp, 0.000, kBlack, 2, 1); + + + +// Kdp.SetExcitation4(-1.0); +// TGraph* gKine= Kdp.GetKinematicLine3(); +// //gKine->Draw("c same"); + + TList* listFull = new TList(); + TList* listNoBG = new TList(); + + for(int i=0; i<26; i++){ + string nameFull = "GateFull_"+to_string(104+(i*2))+"-"+to_string(106+(i*2)); + string nameNoBG = "GateNoBG_"+to_string(104+(i*2))+"-"+to_string(106+(i*2)); + + //double min = gKine->Eval(105+(2*i)); + //cout << "MIN = " << min << endl; + + TH1F* hProj = (TH1F*) hELTL->ProjectionY("hProj",i,i+1); + //hProj->GetXaxis()->SetRangeUser(min,11.); + hProj->GetXaxis()->SetRangeUser(-10.,-0.7); + hProj->Fit("expo","s"); + hProj->GetXaxis()->SetRangeUser(-15.,10.); + //fitPtr->Draw("same"); + + + //TF1* fitFunc = fitPtr->FittedFunction(); + TF1* fitFunc = (TF1*) hProj->GetListOfFunctions()->FindObject("expo"); + + hProj->SetName(nameFull.c_str()); + hProj->SetTitle(nameFull.c_str()); + + listFull->Add(hProj); + + TH1F* hNoBG = (TH1F*) hProj->Clone(); + + for(int j=0; j<hProj->GetNbinsX(); j++){ + hNoBG->SetBinContent(j, hProj->GetBinContent(j) - fitFunc->Eval(105+(i*2))); + hNoBG->SetName(nameNoBG.c_str()); + hNoBG->SetTitle(nameNoBG.c_str()); + + cout << j << " " + << hProj->GetBinContent(j) << " - " + << fitFunc->Eval(105+(i*2)) << " = " + << hNoBG->GetBinContent(j) << endl; + } + + listNoBG->Add(hProj); + + } + + TFile* file = new TFile("Testing_SubtractBackgroundELTL.root","RECREATE"); + listFull->Write("BeforeBGSubtraction",TObject::kSingleKey); + listNoBG->Write("AfterBGSubtraction",TObject::kSingleKey); + file->ls(); + + +} + +TChain* ch3p0=NULL ; +TChain* ch5p5=NULL ; + +void Testing_CompareExForDiffThicknesses(){ + + /////////////////////////////////////////////////////////// + vector<string> files; + string gate = "Mugast.TelescopeNumber>0 && Ex@.size()==1 && abs(T_MUGAST_VAMOS-270)<40"; + + /////////////////////////////////////////////////////////// + files.push_back("../../../Outputs/Analysis/47Kdp_18Oct22_3p0um_PartI.root"); + files.push_back("../../../Outputs/Analysis/47Kdp_18Oct22_3p0um_PartII.root"); + ch3p0 = Chain("PhysicsTree",files,true); + files.clear(); + + ch3p0->Draw("Ex>>h3p0(600,-1,5)", gate.c_str(), ""); + TH1F* h3p0 = (TH1F*) gDirectory->Get("h3p0"); + h3p0->SetLineColor(kBlack); + + /////////////////////////////////////////////////////////// + files.push_back("../../../Outputs/Analysis/47Kdp_25Jul23_5p5um_PartI.root"); + files.push_back("../../../Outputs/Analysis/47Kdp_25Jul23_5p5um_PartII.root"); + ch5p5 = Chain("PhysicsTree",files,true); + files.clear(); + + ch5p5->Draw("Ex>>h5p5(600,-1,5)", gate.c_str(), ""); + TH1F* h5p5 = (TH1F*) gDirectory->Get("h5p5"); + h5p5->SetLineColor(kRed); + + /////////////////////////////////////////////////////////// + + + h3p0->Draw(); + h5p5->Draw("same"); } +void Test_RemoveBG_2D(){ + vector<TH1F*> vDirty; + vector<TH1F*> vClean; + vector<TH1F*> vBackG; + + string dirtygate = "Mugast.TelescopeNumber>0 && Mugast.DSSD_T<770 && abs(T_MUGAST_CATS2-230)<200 && Ex@.size()==1"; + + string cleangate = "Mugast.TelescopeNumber>0 && abs(T_MUGAST_VAMOS-270)<40 && Ex@.size()==1"; + chain->Draw("Ex:ThetaLab>>hDirty2D(26,104,156,500,-15,10)", + dirtygate.c_str(), + ""); + TH2F* hDirty2D = (TH2F*) gDirectory->Get("hDirty2D"); + chain->Draw("Ex:ThetaLab>>hClean2D(26,104,156,500,-15,10)", + cleangate.c_str(), + ""); + TH2F* hClean2D = (TH2F*) gDirectory->Get("hClean2D"); + + for(int i=0; i<26; i++){ + string nameD = "hDirty_" + to_string(104+(i*2)) + "to" + to_string(106+(i*2)); + string nameC = "hClean_" + to_string(104+(i*2)) + "to" + to_string(106+(i*2)); + string nameB = "hBackG_" + to_string(104+(i*2)) + "to" + to_string(106+(i*2)); + + //take cut for each angle bin + TH1F* hDirty = (TH1F*) hDirty2D->ProjectionY(nameD.c_str(), i, i); + TH1F* hClean = (TH1F*) hClean2D->ProjectionY(nameC.c_str(), i, i); + + hDirty->SetLineColor(kBlack); hDirty->SetLineWidth(2); + hClean->SetLineColor(kGreen); hClean->SetLineWidth(2); + + //focus on 0MeV peak, fit straight line bg to dirty to find "real" peak height + hDirty->GetXaxis()->SetRangeUser(-4.0,0.7); + hClean->GetXaxis()->SetRangeUser(-4.0,0.7); + + TF1* fitGausLinear = new TF1("fitGausLinear","[0] + [1]*x + gaus(2)",-4.0,0.7); + fitGausLinear->SetParName(0,"Intercept"); + fitGausLinear->SetParName(1,"Gradient"); + fitGausLinear->SetParName(2,"Constant"); fitGausLinear->SetParameter(2, 1e2); + fitGausLinear->SetParName(3,"Mean"); fitGausLinear->SetParameter(3, 0.1); + fitGausLinear->SetParName(4,"Sigma"); fitGausLinear->SetParameter(4,0.15); + TFitResultPtr fitDirty = hDirty->Fit(fitGausLinear, "sq"); + + TF1* fitGaus = new TF1("fitGaus","gaus(0)",-4.0,0.7); + fitGaus->SetParName(0,"Constant"); fitGaus->SetParameter(0, 1e2); + fitGaus->SetParName(1,"Mean"); fitGaus->SetParameter(1, 0.1); + fitGaus->SetParName(2,"Sigma"); fitGaus->SetParameter(2,0.15); + TFitResultPtr fitClean = hClean->Fit(fitGaus, "sq"); + + vDirty.push_back(hDirty); + vClean.push_back(hClean); + + //maybe align largest bins? + //scale clean to dirty and subtract + double heightD = fitDirty->Parameter(2); + double heightC = fitClean->Parameter(0); + double ratioCD = fitClean->Parameter(0) / fitDirty->Parameter(2); + cout << "Height Dirty = " << heightD + << "\tHeight Clean = " << heightC + << "\tRatio = " << ratioCD + << endl; + TH1F* hBackG = (TH1F*) hDirty->Clone(); + hBackG->SetTitle(nameB.c_str()); + hBackG->SetName(nameB.c_str()); + hBackG->SetLineColor(kRed); + TH1F* hBGrmv = (TH1F*) hClean->Clone(); + hBGrmv->SetTitle("xx"); + hBGrmv->SetName("xx"); + hBGrmv->Scale(ratioCD); + hBackG->Add(hBGrmv,-1); + vBackG.push_back(hBackG); + } + + +} diff --git a/Projects/e793s/macro/ElasticsFitELabGates.h b/Projects/e793s/macro/ElasticsFitELabGates.h new file mode 100644 index 0000000000000000000000000000000000000000..ca21dad73e43710c92927709e2c4fb403f40967f --- /dev/null +++ b/Projects/e793s/macro/ElasticsFitELabGates.h @@ -0,0 +1,93 @@ +void ElasticsFitAngleGates(string ddpp, double MinThetaCM, double MaxThetaCM, double StepWidthThetaCM){ + double Amp, PShft; + int waitpoint; + TList* list = new TList(); + + auto t = time(nullptr); + auto tm = *localtime(&t); + + cout << " Elastics Equation:: " << endl; + cout << " Amplitude / PShift? " << endl; + + if(ddpp == "dd"){ + cout << GREEN << " (d,d) T = 11Jul22" << RESET << endl; + Amp = 27.817908; + PShft = 27.817762; + } else if (ddpp == "pp"){ + cout << GREEN << " (p,p) T = 11Jul22" << RESET << endl; + Amp = 14.947797; + PShft = 14.947895; + } + + int nIter = (int) (MaxThetaCM-MinThetaCM)/StepWidthThetaCM; + cout << "Number of steps: " << nIter << endl; + + ostringstream ossText; + ossText << "ElasticsFitGates/outputElasticsAreas_" + << ddpp << "_" + << put_time(&tm, "%Y-%m-%d_%H%M") + << ".root"; + auto strText = ossText.str(); + + ofstream outfile; + outfile.open (strText.c_str()); + + for(int i=0; i<nIter; i++){ + double MinELab = (Amp * cos((180.-(MinThetaCM+(i*StepWidthThetaCM))) * (M_PI/180.))) + PShft; + double MaxELab = (Amp * cos((180.-(MinThetaCM+((i+1)*StepWidthThetaCM))) * (M_PI/180.))) + PShft; + + cout << RESET; + cout << "==========================================================" << endl; + cout << MinELab << " " << MaxELab << endl; + cout << "----------------------------------------------------------" << endl; + cout << "Gate: " << MinThetaCM+(i*StepWidthThetaCM) + << " to " << MinThetaCM+((i+1)*StepWidthThetaCM) + << endl; + cout << "----------------------------------------------------------" << endl; + cout << "Constructing histogram..." << endl; + + string histname = "Gate_" + to_string(MinThetaCM+(i*StepWidthThetaCM)) + "-" + to_string(MinThetaCM+((i+1)*StepWidthThetaCM)); + + ElasticsGate(MinELab, MaxELab); + TH1F* hist = (TH1F*) gDirectory->Get("hist"); + hist->Draw(); + hist->SetName(histname.c_str()); + gPad->Update(); + list->Add(hist); + cout << "----------------------------------------------------------" << endl; + + double ddLine, ppLine; + if(ddpp == "dd"){ + ddLine = 81.-5.5*log(MaxELab); + ppLine = 79.-7.5*log(MaxELab); + }else if (ddpp == "pp"){ + ddLine = 83.-7.5*log(MaxELab); + ppLine = 82.-10.5*log(MaxELab); + } + + cout << "Fitting with dd = " << ddLine << " and pp = " << ppLine << endl; + vector<double> areas = DoubleGausNumbs(hist, 52., 87., ppLine, ddLine, true); + + outfile << areas.at(0) << "\t" + << areas.at(1) << "\t" + << areas.at(2) << "\t" + << areas.at(3) << endl; + + gPad->Update(); + + } + + outfile.close(); + + ostringstream oss; + oss << "ElasticsFitGates/ElasticsFitGates_" + << ddpp << "_" + << put_time(&tm, "%Y-%m-%d_%H%M") + << ".root"; + auto str = oss.str(); + + TFile* file = new TFile(str.c_str(),"RECREATE"); + list->Write("ElasticsGateFits",TObject::kSingleKey); + file->ls(); + +} diff --git a/Projects/e793s/macro/KnownPeakFitter.h b/Projects/e793s/macro/KnownPeakFitter.h index c7e685674d12003fd7416e75a64d08c0ee90a9e5..bb56484def43e6afceb1ceff0c9a657a8f93497d 100644 --- a/Projects/e793s/macro/KnownPeakFitter.h +++ b/Projects/e793s/macro/KnownPeakFitter.h @@ -4,12 +4,14 @@ #include "stdlib.h" bool removing = false; +bool bgForCS2 = true; vector<double> fitContByBin; //each entry is a bin, containing eval of fit function in that bin for removal function +int numPeaks_temp = 0; vector<double> means = { 0.000, 0.143, - 0.279, - 0.728, + //0.279, + //0.728, 0.967, 1.409, 1.978, @@ -25,24 +27,28 @@ vector<double> means = { 0.000, }; const int numPeaks = means.size(); - -vector<double> means_dt = { 0.000, - 0.587, - 0.691, +vector<double> means_dt = { + 0.000, + 0.639,//FIT AS DOUBLET + //0.587, + //0.691, //0.886, + //1.370, + //1.554,//TESTING HALFWAY BETWEEN 1.37 AND 1.74 //1.738, 1.944, 2.233, 2.732, 3.377,//DOUBLET + //3.330,//3.377,//DOUBLET //3.344, - //3.410, - 4.297 + //3.409, + 4.15 + //4.15//4.297 }; const int numPeaks_dt = means_dt.size(); - array<double,27> knowngammas = { 0.143, 0.279, 0.449, @@ -93,61 +99,124 @@ Double_t f_peak(Double_t *x, Double_t *par){ */ string f_full_string(){ + string result = "gaus(0)"; + //for(int pk=1; pk<numPeaks; pk++){ + for(int pk=1; pk<numPeaks_temp; pk++){ + result += "+gaus(" + to_string(3*pk) + ")"; + } + return result; +} + +string f_full_pol2bg_string(){ string result = "gaus(0)"; for(int pk=1; pk<numPeaks; pk++){ result += "+gaus(" + to_string(3*pk) + ")"; } + result += " + [" + to_string(3*numPeaks+0) + + "] + [" + to_string(3*numPeaks+1) + + "]*x + [" + to_string(3*numPeaks+2) + + "]*x*x"; return result; } -Double_t f_full(Double_t *x, Double_t *par) { - float xx = x[0]; - double result, norm; - result = 0; +string f_full_pol4bg_string(){ + string result = "gaus(0)"; + for(int pk=1; pk<numPeaks; pk++){ + result += "+gaus(" + to_string(3*pk) + ")"; + } +//[0] * ( [1] + [2]*x + [3]*x*x + [4]*x*x*x + [5]*x*x*x*x) + result += " + ([" + to_string(3*numPeaks+1) + "] " + + " + [" + to_string(3*numPeaks+2) + "]*x " + + " + [" + to_string(3*numPeaks+3) + "]*x*x " + + " + [" + to_string(3*numPeaks+4) + "]*x*x*x " + + " + [" + to_string(3*numPeaks+5) + "]*x*x*x*x " + + " ) * " + to_string(3*numPeaks+0) + "]"; - // Add N peaks - for(int pk=0; pk<numPeaks; pk++){ - result += (par[1+(pk*3)]*exp(-0.5*pow(((xx-par[2+(pk*3)]-par[0])/par[3+(pk*3)]),2))); //REGULAR GAUS + + + +// result += " + [" + to_string(3*numPeaks+0) +// + "] *([" + to_string(3*numPeaks+1) +// + "] + [" + to_string(3*numPeaks+1) +// + "]*x + [" + to_string(3*numPeaks+2) +// + "]*x*x"; + return result; +} + +string f_full_skewgaussbg_string(){ + string result = "gaus(0)"; + for(int pk=1; pk<numPeaks; pk++){ + result += "+gaus(" + to_string(3*pk) + ")"; } +//"2.*gaus(x-[3],[0],0,[1])*ROOT::Math::normal_cdf([2]*(x-[3]),1,0)" + result += " + 2.*gaus(x-[" + + to_string(3*numPeaks+3) + "], [" + + to_string(3*numPeaks+0) + "], 0 ,[" + + to_string(3*numPeaks+1) + "])*ROOT::Math::normal_cdf([" + + to_string(3*numPeaks+2) + "]*(x-[" + + to_string(3*numPeaks+3) + "]),1,0)"; + return result; } -Double_t f_full_dt(Double_t *x, Double_t *par) { +Double_t f_full(Double_t *x, Double_t *par) { float xx = x[0]; double result, norm; - // Flat background -// result = par[0]; result = 0; + // Add N peaks - for(int pk=0; pk<numPeaks_dt; pk++){ - result += (par[3+(pk*3)]/(par[1+(pk*3)]*sqrt(2*pi))) - //* exp(-0.5*pow((xx-par[2+(pk*3)])/par[1+(pk*3)],2)); - * exp(-0.5*pow((xx-par[2+(pk*3)]-par[0])/par[1+(pk*3)],2)); //added par 0 as shift in energy + for(int pk=0; pk<numPeaks; pk++){ + result += (par[1+(pk*3)]*exp(-0.5*pow(((xx-par[2+(pk*3)]-par[0])/par[3+(pk*3)]),2))); //REGULAR GAUS } return result; } +//Double_t f_full_dt(Double_t *x, Double_t *par) { +// float xx = x[0]; +// double result, norm; +// // Flat background +//// result = par[0]; +// result = 0; +// // Add N peaks +// for(int pk=0; pk<numPeaks_dt; pk++){ +// result += (par[3+(pk*3)]/(par[1+(pk*3)]*sqrt(2*pi))) +// //* exp(-0.5*pow((xx-par[2+(pk*3)])/par[1+(pk*3)],2)); +// * exp(-0.5*pow((xx-par[2+(pk*3)]-par[0])/par[1+(pk*3)],2)); //added par 0 as shift in energy +// } +// return result; +//} + vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist, double slideshift){ - double minFit=-1.0, maxFit=4.644; + double minFit, maxFit; double binWidth = hist->GetXaxis()->GetBinWidth(3); //double sigma = 0.14; double sigma = 0.0;// = 0.14; - int numPeaks_temp = 0; + numPeaks_temp = 0; vector<double> means_temp; //Gradient on mean value double gradient; + double intercept; if(reactionName=="47K(d,p)"){ sigma = 0.14; - gradient = 0.989; + //gradient = 0.989; + gradient = 1.0;//0.989; + intercept = 0.000; numPeaks_temp = numPeaks; means_temp = means; + minFit = -1.0; + maxFit = 4.644; } else if (reactionName=="47K(d,t)"){ - sigma = 0.18; - gradient = 0.986; + sigma=0.16;//sigma = 0.18; + //gradient = 0.986; + gradient = 1.002; //0.995;//1.005; + //gradient = 1.002; //0.995;//1.005; + intercept = -0.060; //-0.055;//-0.061; numPeaks_temp = numPeaks_dt; means_temp = means_dt; + minFit = -1.; + maxFit = +10.; } cout << BOLDGREEN @@ -166,35 +235,100 @@ vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist, double slideshift){ string nameHere = nameBase; nameHere +=to_string(i); allPeaks[i] = new TF1(nameHere.c_str(), function.c_str(), minFit, maxFit); - allPeaks[i]->SetLineColor(kBlack); + allPeaks[i]->SetLineColor(4001); allPeaks[i]->SetLineStyle(7); allPeaks[i]->SetNpx(500); } //Subtract flat background equal to smallest bin in range if(!removing){ - hist->GetXaxis()->SetRange(hist->FindBin(-0.8), hist->FindBin(-0.0)); - double bgmin = hist->GetBinContent(hist->GetMinimumBin()); - hist->GetXaxis()->UnZoom(); - cout << "Subtracting background of " << bgmin << endl; - for(int b=1; b<hist->GetNbinsX() ; b++){ - hist->SetBinContent(b,hist->GetBinContent(b)-bgmin); + TF1 *flat = new TF1 ("flat","[0]",5,8); + cout << RED << "REMOVING FLAT BACKGROUND!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << RESET << endl; + if(reactionName=="47K(d,p)"){ + hist->GetXaxis()->SetRange(hist->FindBin(-0.8), hist->FindBin(-0.0)); + double bgmin = hist->GetBinContent(hist->GetMinimumBin()); + hist->GetXaxis()->UnZoom(); + cout << "Subtracting background of " << bgmin << endl; + for(int b=1; b<hist->GetNbinsX() ; b++){ + hist->SetBinContent(b,hist->GetBinContent(b)-bgmin); + } + } else if(reactionName=="47K(d,t)"){ +// int count = 0; +// int sum = 0; +// cout << "here" << endl; +// for(int b=hist->FindBin(5.0); b<hist->FindBin(7.0) ; b++){ +// sum += hist->GetBinContent(b); +// count++; +// } +// double avg = (double)count / (double)sum; +// cout << "Subtracting background of " << avg << endl; +// for(int b=1; b<hist->GetNbinsX() ; b++){ +// hist->SetBinContent(b,hist->GetBinContent(b)-avg); +// if(hist->GetBinContent(b)<0){hist->SetBinContent(b,0);} +// } + } + } //Build background function - TF1 *bg = new TF1 ("bg","[0]",minFit, maxFit); + string backend_bgfuncstring; + //if(bgForCS2){ + backend_bgfuncstring = "]"; //cout << "fitting with FLAT bg because makes no diff for CS2" << endl; + //} else { + //backend_bgfuncstring = "]-0.478*(x-5.0)^2"; //pol2 + //backend_bgfuncstring = "]*exp(-(x-4.5)^2/(2*(2.9^2)))"; //gaus + //backend_bgfuncstring = "]/4.5 * sqrt((4.5)^2 - (x-4.5)^2)"; //semicirc + //} + string bgfuncstring; + double bgminFit, bgmaxFit; + if(reactionName=="47K(d,p)"){ + bgfuncstring = "[0]"; + bgminFit = minFit; + bgmaxFit = maxFit; + } else if(reactionName=="47K(d,t)"){ + bgfuncstring = "[0"; + if(!bgForCS2){bgfuncstring += backend_bgfuncstring;} + else{bgfuncstring += "]";} + bgminFit = 5.; + bgmaxFit = 7.; + } + TF1 *bg = new TF1 ("bg",bgfuncstring.c_str(),bgminFit, bgmaxFit); bg->SetLineColor(kGreen); bg->SetLineStyle(9); bg->SetParNames("Background"); - + //bg->SetParameter(0,10); + //bg->SetParameter(0,-10); + if(bgForCS2){bg->FixParameter(0,0);} + hist->Fit(bg, "RBL", "", bgminFit, bgmaxFit); + bg->SetRange(minFit,maxFit); + cout << bg->GetParameter(0) << endl; + //Build total function + int bgRefNum = ((numPeaks_temp-1)*3)+3; string s_full = f_full_string(); + if(reactionName=="47K(d,p)"){ + //background? + } else if(reactionName=="47K(d,t)"){ +// if(!bgForCS2){ +// string s_full = f_full_string(); + s_full.append(" + (["); + s_full+=to_string(bgRefNum); + //s_full.append("]/4.5 * sqrt(25.0 - (x-4.5)^2))"); +// s_full+=backend_bgfuncstring; + s_full.append("])"); +// } + } + + cout << "BG REF NUMBER " << bgRefNum << endl; + cout << s_full << endl; + TF1 *full = new TF1("fitAllPeaks", s_full.c_str(), minFit, maxFit); - full->SetLineColor(kRed); + full->SetLineColor(4000); const int numParams = (numPeaks_temp*3); for(int i=0; i<numPeaks_temp; i++) { - full->FixParameter((i*3)+1,means_temp.at(i)*gradient); + //full->FixParameter((i*3)+1,means_temp.at(i)*gradient); + full->FixParameter((i*3)+1,(means_temp.at(i)*gradient) + intercept); full->FixParameter((i*3)+2,sigma); full->SetParameter((i*3)+0,1e1); full->SetParLimits((i*3)+0,0.0,1e5); @@ -203,18 +337,57 @@ vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist, double slideshift){ //Set and fix specific parameters for (int i=0; i<numPeaks; i++){ - if(means[i]>4.0){ - full->SetParameter((i*3)+2,0.23); - } - if(abs(means[i]-0.279)<0.05){ - full->FixParameter((i*3)+0,0.0); + if(reactionName=="47K(d,p)"){ + if(means[i]>4.0){ + full->SetParameter((i*3)+2,0.23); + } + if(abs(means[i]-0.279)<0.05){ + full->FixParameter((i*3)+0,0.0); + } + if(abs(means[i]-0.728)<0.05){ + full->FixParameter((i*3)+0,0.0); + } + } else if(reactionName=="47K(d,t)"){ + full->SetParLimits(bgRefNum,0,30); +// full->FixParameter(bgRefNum,0); + //if(abs(means[i]-1.55)<0.3){ + // cout << "TESTING FORCE FOR FWAVE STATES!!!" << endl; + // full->SetParameter((i*3)+0,10.); + // full->SetParLimits((i*3)+0,0.0,50.); + //} + if(means[i]>4.0){ + cout << "HERE!! " << endl; + full->SetParameter((i*3)+0,10.); + full->SetParLimits((i*3)+0,0.0,50.); + } +// if(bgForCS2){full->FixParameter(bgRefNum,0);} +// else { +// full->SetParLimits(bgRefNum,bg->GetParameter(0),bg->GetParameter(0)); +// full->FixParameter(bgRefNum,bg->GetParameter(0)); +// cout << "fixing to ...... " << bg->GetParameter(0) << endl; +// } + //cout << "HERE!! IN DT FIT!!! " << bg->GetParameter(0) << endl; + //cout << "HERE!! IN DT FIT!!! " << bg->GetParameter(0) << endl; + //cout << "HERE!! IN DT FIT!!! " << bg->GetParameter(0) << endl; + //cout << "HERE!! IN DT FIT!!! " << bg->GetParameter(0) << endl; + //cout << "HERE!! IN DT FIT!!! " << bg->GetParameter(0) << endl; + //cout << "HERE!! IN DT FIT!!! " << bg->GetParameter(0) << endl; + //cout << "HERE!! IN DT FIT!!! " << bg->GetParameter(0) << endl; + //cout << "HERE!! IN DT FIT!!! " << bg->GetParameter(0) << endl; + //cout << "HERE!! IN DT FIT!!! " << bg->GetParameter(0) << endl; + + //full->SetParLimits(bgRefNum,8.,8.); + //full->FixParameter(bgRefNum,8.); + //cout << "fixing to ...... " << 8. << endl; } } //Fit full function to histogram - //hist->Fit(full, "RWQB", "", minFit, maxFit); - hist->Fit(full, "RWQBL", "", minFit, maxFit); cout << GREEN << "FITTING WITH MAX LIKELYHOOD METHOD" << RESET << endl; - + if(reactionName=="47K(d,p)"){ + hist->Fit(full, "RWQBL", "", minFit, maxFit); cout << GREEN << "FITTING WITH MAX LIKELYHOOD METHOD" << RESET << endl; + } else if(reactionName=="47K(d,t)"){ + hist->Fit(full, "RWBL", "", minFit, maxFit); cout << GREEN << "FITTING WITH MAX LIKELYHOOD METHOD" << RESET << endl; + } hist->Draw(); //Extract fitted variables, assign them to individual fits, and draw them @@ -230,27 +403,36 @@ vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist, double slideshift){ //Make irresolvable doublet clear if(reactionName=="47K(d,p)"){ for (int i=0; i<numPeaks_temp; i++){ - if(abs(means[i]-3.83)<0.1){ + if(abs(means_temp[i]-3.83)<0.1){ cout << BOLDCYAN << "PEAK #" << i << " IS FITTING IRRESOLVABLE DOUBLET!" << RESET << endl; - allPeaks[i]->SetLineColor(kCyan+1); + allPeaks[i]->SetLineColor(4003); } - if(means[i]>4.0){ + if(means_temp[i]>4.0){ cout << BOLDMAGENTA << "PEAK #" << i << " IS FITTING UNKNOWN ZONE WITH SIGMA = " << full->GetParameter((i*3)+2) << RESET << endl; - allPeaks[i]->SetLineColor(kMagenta); + allPeaks[i]->SetLineColor(4004); } } } else if(reactionName=="47K(d,t)"){ for (int i=0; i<numPeaks_temp; i++){ - if(abs(means[i]-3.37)<0.1){ + cout << "here" << endl; + if(abs(means_temp[i]-3.377)<0.10){ + cout << BOLDCYAN + << "PEAK #" << i + << " IS FITTING IRRESOLVABLE DOUBLET!" + << RESET + << endl; + allPeaks[i]->SetLineColor(kCyan+1); + } + if(abs(means_temp[i]-0.639)<0.10){ cout << BOLDCYAN << "PEAK #" << i << " IS FITTING IRRESOLVABLE DOUBLET!" @@ -258,6 +440,28 @@ vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist, double slideshift){ << endl; allPeaks[i]->SetLineColor(kCyan+1); } +// +// +// +// +// +// if(abs(means[i]-3.37)<0.1){ +// cout << BOLDCYAN +// << "PEAK #" << i +// << " IS FITTING IRRESOLVABLE DOUBLET!" +// << RESET +// << endl; +// allPeaks[i]->SetLineColor(kCyan+1); +// } +// if(abs(means[i]-0.639)<0.1){ +// cout << BOLDCYAN +// << "PEAK #" << i +// << " IS FITTING IRRESOLVABLE DOUBLET!" +// << RESET +// << endl; +// allPeaks[i]->SetLineColor(kCyan+1); +// } +// } } @@ -304,10 +508,8 @@ vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist, double slideshift){ } if(removing){ - - for(int b=1; b<hist->GetNbinsX()-1; b++){ - cout << b << " " << hist->GetBinContent(b) << " " << full->Eval(hist->GetBinCenter(b)) << endl; + //cout << b << " " << hist->GetBinContent(b) << " " << full->Eval(hist->GetBinCenter(b)) << endl; fitContByBin.push_back(full->Eval(hist->GetBinCenter(b))); } } @@ -320,7 +522,8 @@ void FitKnownPeaks(TH1F* hist){ vector<vector<double>> shell = FitKnownPeaks_RtrnArry(hist,0.0); } -void FitKnownPeaks_dt(TH1F* hist){ +void FitKnownPeaks_dt(TH1F* hist, bool forCS2){ + bgForCS2 = forCS2; //Shell function to call Rtrn_Arry without writing vector<vector<double>> to screen //vector<vector<double>> shell = FitKnownPeaks_dt_RtrnArry(hist,0.0); vector<vector<double>> shell = FitKnownPeaks_RtrnArry(hist,0.0); @@ -347,3 +550,653 @@ TH1F* RemoveKnownPeaks(TH1F* hist){ return hist; } +vector<vector<double>> FitKnownPeaks_Pol2BG(TH1F* hist){ + //double minFit=-14.0, maxFit=4.644; + double minFit, maxFit; + double bgminfit, bgmaxfit; + double binWidth = hist->GetXaxis()->GetBinWidth(3); + + double sigma = 0.0;// = 0.14; + int numPeaks_temp = 0; + vector<double> means_temp; + + //Gradient on mean value + double gradient; + if(reactionName=="47K(d,p)"){ + sigma = 0.14; + gradient = 0.989; + //gradient = 0.968; + numPeaks_temp = numPeaks; + means_temp = means; + minFit=-10.0; + maxFit=4.644; + bgminfit=minFit; + bgmaxfit=-0.7; + } else if (reactionName=="47K(d,t)"){ + sigma = 0.18; + gradient = 0.986; + numPeaks_temp = numPeaks_dt; + means_temp = means_dt; + minFit=-1.0; + maxFit=10.0; + bgminfit=5.0; + bgmaxfit=maxFit; + } + + string histnamepull = hist->GetName(); + cout << histnamepull << " -> " << histnamepull.substr(14,3) << endl; + int anglenumber = stoi(histnamepull.substr(14,3)); + + + vector<double> means_GradientScaled; + for(int i=0; i<means_temp.size(); i++){ + //if(anglenumber < 120){ + //means_GradientScaled.push_back((means_temp.at(i) * 0.987) - 0.039); + //} else { + // means_GradientScaled.push_back(((means_temp.at(i) * 0.987) - 0.039) * 0.0015*(anglenumber) + 1.154); + //} + means_GradientScaled.push_back(means_temp.at(i) * gradient); + } + + cout << BOLDGREEN + << "FITTING WITH A REAL:OBSERVED Ex GRADIENT GRADIENT OF " + //<< gradient + << " CUSTOM EQUATION!!!!!!" + << RESET + << endl; + + //Build individual peak fit functions + string nameBase = "Peak "; + string function = "gaus";//"[0]*exp(-0.5*pow(((x-[1])/[2]),2))";//REGULAR GAUS + TF1 **allPeaks = new TF1*[numPeaks_temp]; + for(int i=0; i<numPeaks_temp; i++) { + string nameHere = nameBase; + nameHere +=to_string(i); + allPeaks[i] = new TF1(nameHere.c_str(), function.c_str(), minFit, maxFit); + allPeaks[i]->SetLineColor(kBlack); + allPeaks[i]->SetLineStyle(7); + allPeaks[i]->SetNpx(500); + } + + //Build background function + TF1 *bg = new TF1 ("bg","[0] + [1]*x + abs([2])*x*x", bgminfit, bgmaxfit);//-1.5); + bg->SetParLimits(2,0,1e2); + bg->SetLineColor(kGreen); + //bg->SetLineStyle(9); + bg->SetParNames("Background"); + hist->Fit(bg, "RWQBL", "", bgminfit, bgmaxfit); + cout << bg->GetParameter(0) << " + " + << bg->GetParameter(1) << "x + " + << bg->GetParameter(2) << "x*x" + << endl; + bg->SetRange(minFit,maxFit); + + //Build total function + string s_full = f_full_pol2bg_string(); + TF1 *full = new TF1("fitAllPeaks", s_full.c_str(), minFit, maxFit); + full->SetLineColor(kRed); + const int numParams = (numPeaks_temp*3); + for(int i=0; i<numPeaks_temp; i++) { + //full->FixParameter((i*3)+1,means_temp.at(i)*gradient); + full->FixParameter((i*3)+1,means_GradientScaled.at(i)); + full->FixParameter((i*3)+2,sigma); + full->SetParameter((i*3)+0,1e1); + full->SetParLimits((i*3)+0,0.0,1e5); + } + full->FixParameter(3*numPeaks_temp+0,bg->GetParameter(0)); + full->FixParameter(3*numPeaks_temp+1,bg->GetParameter(1)); + full->FixParameter(3*numPeaks_temp+2,bg->GetParameter(2)); + full->SetNpx(500); + + //Set and fix specific parameters + for (int i=0; i<numPeaks_temp; i++){ + if(reactionName=="47K(d,p)"){ + if(means[i]>4.0){ + full->SetParameter((i*3)+2,0.23); + } + if(abs(means[i]-0.279)<0.05){ + full->FixParameter((i*3)+0,0.0); + } + } else if(reactionName=="47K(d,t)"){ + if(abs(means[i]-1.55)<0.3){ + full->SetParameter((i*3)+2,0.23); + } + } + } + + //Fit full function to histogram + hist->Fit(full, "RWQBL", "", minFit, maxFit); + cout << GREEN << "FITTING WITH MAX LIKELYHOOD METHOD" << RESET << endl; + hist->Draw(); + + //Extract fitted variables, assign them to individual fits, and draw them + const Double_t* finalPar = full->GetParameters(); + const Double_t* finalErr = full->GetParErrors(); + for (int i=0; i<numPeaks_temp; i++){ + allPeaks[i]->SetParameters(finalPar[0+(i*3)], + finalPar[1+(i*3)], + finalPar[2+(i*3)]); + } + bg->Draw("SAME"); + full->Draw("SAME"); + + //Make irresolvable doublet clear + if(reactionName=="47K(d,p)"){ + for (int i=0; i<numPeaks_temp; i++){ + if(abs(means[i]-3.83)<0.1){ + cout << BOLDCYAN + << "PEAK #" << i + << " IS FITTING IRRESOLVABLE DOUBLET!" + << RESET + << endl; + allPeaks[i]->SetLineColor(kCyan+1); + } + if(means[i]>4.0){ + cout << BOLDMAGENTA + << "PEAK #" << i + << " IS FITTING UNKNOWN ZONE WITH SIGMA = " + << full->GetParameter((i*3)+2) + << RESET + << endl; + allPeaks[i]->SetLineColor(kMagenta); + } + } + } else if(reactionName=="47K(d,t)"){ + for (int i=0; i<numPeaks_temp; i++){ + if(abs(means[i]-3.37)<0.1){ + cout << BOLDCYAN + << "PEAK #" << i + << " IS FITTING IRRESOLVABLE DOUBLET!" + << RESET + << endl; + allPeaks[i]->SetLineColor(kCyan+1); + } + } + } + + //Draw all peaks + for (int i=0; i<numPeaks_temp; i++){ + allPeaks[i]->Draw("SAME"); + } + + //Write to screen + cout << "===========================" << endl; + cout << "== PEAK =========== AREA ==" << endl; + + //Loop over every peak + vector<vector<double>> allpeaks; + for(int i=0; i<numPeaks_temp; i++){ + //For AREA = HEIGHT * SIGMA * SQRT(2*PI) + double A = (finalPar[0+(i*3)] * finalPar[2+(i*3)] * sqrt(2*pi)) / binWidth; + + double deltaA = A * sqrt( pow( finalErr[0+(i*3)]/finalPar[0+(i*3)] ,2) + + pow( finalErr[2+(i*3)]/finalPar[2+(i*3)] ,2) ); + +// cout << "DELTAAREA = " +// << finalErr[0+(i*3)] << " / " +// << finalPar[0+(i*3)] << " and " +// << finalErr[2+(i*3)] << " / " +// << finalPar[2+(i*3)] << " = " +// << deltaA << endl; + + cout << fixed << setprecision(3) + << " #" << i << " " + << finalPar[(i*3)+1] << "\t" << setprecision(0) + << A << "\t+- " + << deltaA << setprecision(3); + cout << " SQRT: " << sqrt(A) << endl; + + vector<double> onepeak; //energy, area and error for one peak + onepeak.push_back(finalPar[(i*3)+1]); + onepeak.push_back(A); + onepeak.push_back(deltaA); + allpeaks.push_back(onepeak); + cout << "-------------" << endl; + } + + if(removing){ + for(int b=1; b<hist->GetNbinsX()-1; b++){ + //cout << b << " " << hist->GetBinContent(b) << " " << full->Eval(hist->GetBinCenter(b)) << endl; + fitContByBin.push_back(full->Eval(hist->GetBinCenter(b))); + } + } + + return allpeaks; +} + +vector<vector<double>> FitKnownPeaks_Pol4BG(TH1F* hist){ + //double minFit=-14.0, maxFit=4.644; + double minFit=-10.0, maxFit=4.644; + double binWidth = hist->GetXaxis()->GetBinWidth(3); + + double sigma = 0.0;// = 0.14; + int numPeaks_temp = 0; + vector<double> means_temp; + + //FIX SHAPE OF POL4 FROM FITTING BG GATE, ONLY ALLOW VERTICAL SCALING + double fix1 = 136.101; + double fix2 = 15.4108; + double fix3 = -0.906266; + double fix4 = -0.0737336; + double fix5 = 0.00402861; + + + //Gradient on mean value + double gradient; + if(reactionName=="47K(d,p)"){ + sigma = 0.14; + gradient = 0.989; + //gradient = 0.968; + numPeaks_temp = numPeaks; + means_temp = means; + } else if (reactionName=="47K(d,t)"){ + sigma = 0.18; + gradient = 0.986; + numPeaks_temp = numPeaks_dt; + means_temp = means_dt; + } + + string histnamepull = hist->GetName(); + cout << histnamepull << " -> " << histnamepull.substr(14,3) << endl; + int anglenumber = stoi(histnamepull.substr(14,3)); + + + vector<double> means_GradientScaled; + for(int i=0; i<means_temp.size(); i++){ + //if(anglenumber < 120){ + //means_GradientScaled.push_back((means_temp.at(i) * 0.987) - 0.039); + //} else { + // means_GradientScaled.push_back(((means_temp.at(i) * 0.987) - 0.039) * 0.0015*(anglenumber) + 1.154); + //} + means_GradientScaled.push_back(means_temp.at(i) * gradient); + } + + cout << BOLDGREEN + << "FITTING WITH A REAL:OBSERVED Ex GRADIENT GRADIENT OF " + //<< gradient + << " CUSTOM EQUATION!!!!!!" + << RESET + << endl; + + //Build individual peak fit functions + string nameBase = "Peak "; + string function = "gaus";//"[0]*exp(-0.5*pow(((x-[1])/[2]),2))";//REGULAR GAUS + TF1 **allPeaks = new TF1*[numPeaks_temp]; + for(int i=0; i<numPeaks_temp; i++) { + string nameHere = nameBase; + nameHere +=to_string(i); + allPeaks[i] = new TF1(nameHere.c_str(), function.c_str(), minFit, maxFit); + allPeaks[i]->SetLineColor(kBlack); + allPeaks[i]->SetLineStyle(7); + allPeaks[i]->SetNpx(500); + } + + //Build background function + TF1 *bg = new TF1 ("bg","[0] * ( [1] + [2]*x + [3]*x*x + [4]*x*x*x + [5]*x*x*x*x)", -8., -1.);//-1.5); + bg->FixParameter(1,fix1); + bg->FixParameter(2,fix2); + bg->FixParameter(3,fix3); + bg->FixParameter(4,fix4); + bg->FixParameter(5,fix5); + bg->SetLineColor(kGreen); + //bg->SetLineStyle(9); + bg->SetParNames("Background"); + hist->Fit(bg, "RWQBL", "", minFit, -1.5); + bg->SetRange(minFit,maxFit); + + //Build total function + string s_full = f_full_pol2bg_string(); + TF1 *full = new TF1("fitAllPeaks", s_full.c_str(), minFit, maxFit); + full->SetLineColor(kRed); + const int numParams = (numPeaks_temp*3); + for(int i=0; i<numPeaks_temp; i++) { + //full->FixParameter((i*3)+1,means_temp.at(i)*gradient); + full->FixParameter((i*3)+1,means_GradientScaled.at(i)); + full->FixParameter((i*3)+2,sigma); + full->SetParameter((i*3)+0,1e1); + full->SetParLimits((i*3)+0,0.0,1e5); + } + full->FixParameter(3*numPeaks_temp+0,bg->GetParameter(0)); + full->FixParameter(3*numPeaks_temp+1,fix1); + full->FixParameter(3*numPeaks_temp+2,fix2); + full->FixParameter(3*numPeaks_temp+3,fix3); + full->FixParameter(3*numPeaks_temp+4,fix4); + full->FixParameter(3*numPeaks_temp+5,fix5); + + full->SetNpx(500); + + //Set and fix specific parameters + for (int i=0; i<numPeaks; i++){ + if(means[i]>4.0){ + full->SetParameter((i*3)+2,0.23); + } + if(abs(means[i]-0.279)<0.05){ + full->FixParameter((i*3)+0,0.0); + } + } + + //Fit full function to histogram + hist->Fit(full, "RWQBL", "", minFit, maxFit); + cout << GREEN << "FITTING WITH MAX LIKELYHOOD METHOD" << RESET << endl; + hist->Draw(); + + //Extract fitted variables, assign them to individual fits, and draw them + const Double_t* finalPar = full->GetParameters(); + const Double_t* finalErr = full->GetParErrors(); + for (int i=0; i<numPeaks_temp; i++){ + allPeaks[i]->SetParameters(finalPar[0+(i*3)], + finalPar[1+(i*3)], + finalPar[2+(i*3)]); + } + bg->Draw("SAME"); + full->Draw("SAME"); + + //Make irresolvable doublet clear + if(reactionName=="47K(d,p)"){ + for (int i=0; i<numPeaks_temp; i++){ + if(abs(means[i]-3.83)<0.1){ + cout << BOLDCYAN + << "PEAK #" << i + << " IS FITTING IRRESOLVABLE DOUBLET!" + << RESET + << endl; + allPeaks[i]->SetLineColor(kCyan+1); + } + if(means[i]>4.0){ + cout << BOLDMAGENTA + << "PEAK #" << i + << " IS FITTING UNKNOWN ZONE WITH SIGMA = " + << full->GetParameter((i*3)+2) + << RESET + << endl; + allPeaks[i]->SetLineColor(kMagenta); + } + } + } else if(reactionName=="47K(d,t)"){ + for (int i=0; i<numPeaks_temp; i++){ + cout << "here" << i << endl; + if(abs(means[i]-3.377)<0.05){ + cout << BOLDCYAN + << "PEAK #" << i + << " IS FITTING IRRESOLVABLE DOUBLET!" + << RESET + << endl; + allPeaks[i]->SetLineColor(kCyan+1); + } + if(abs(means[i]-0.639)<0.05){ + cout << BOLDCYAN + << "PEAK #" << i + << " IS FITTING IRRESOLVABLE DOUBLET!" + << RESET + << endl; + allPeaks[i]->SetLineColor(kCyan+1); + } + + } + } + + //Draw all peaks + for (int i=0; i<numPeaks_temp; i++){ + allPeaks[i]->Draw("SAME"); + } + + //Write to screen + cout << "===========================" << endl; + cout << "== PEAK =========== AREA ==" << endl; + + //Loop over every peak + vector<vector<double>> allpeaks; + for(int i=0; i<numPeaks_temp; i++){ + //For AREA = HEIGHT * SIGMA * SQRT(2*PI) + double A = (finalPar[0+(i*3)] * finalPar[2+(i*3)] * sqrt(2*pi)) / binWidth; + + double deltaA = A * sqrt( pow( finalErr[0+(i*3)]/finalPar[0+(i*3)] ,2) + + pow( finalErr[2+(i*3)]/finalPar[2+(i*3)] ,2) ); + +// cout << "DELTAAREA = " +// << finalErr[0+(i*3)] << " / " +// << finalPar[0+(i*3)] << " and " +// << finalErr[2+(i*3)] << " / " +// << finalPar[2+(i*3)] << " = " +// << deltaA << endl; + + cout << fixed << setprecision(3) + << " #" << i << " " + << finalPar[(i*3)+1] << "\t" << setprecision(0) + << A << "\t+- " + << deltaA << setprecision(3); + cout << " SQRT: " << sqrt(A) << endl; + + vector<double> onepeak; //energy, area and error for one peak + onepeak.push_back(finalPar[(i*3)+1]); + onepeak.push_back(A); + onepeak.push_back(deltaA); + allpeaks.push_back(onepeak); + cout << "-------------" << endl; + } + + if(removing){ + for(int b=1; b<hist->GetNbinsX()-1; b++){ + //cout << b << " " << hist->GetBinContent(b) << " " << full->Eval(hist->GetBinCenter(b)) << endl; + fitContByBin.push_back(full->Eval(hist->GetBinCenter(b))); + } + } + + return allpeaks; +} + +vector<vector<double>> FitKnownPeaks_SkewedGaussBG(TH1F* hist){ + // + //24Aug23 - BG gound from fitting skewed gaussian to 24Aug23_BackgroundShape.root, + //which was itself determined by subtraction old low-bg gate from new high-bg gate. + // + + //double minFit=-14.0, maxFit=4.644; + double minFit=-10.0, maxFit=4.644; + double binWidth = hist->GetXaxis()->GetBinWidth(3); + + double sigma = 0.0;// = 0.14; + int numPeaks_temp = 0; + vector<double> means_temp; + + //FIX SHAPE OF SKEW FROM FITTING BG GATE, ONLY ALLOW VERTICAL SCALING + double fix1 = +6.54282e+00; + double fix2 = -3.17214e-01; + double fix3 = +5.64278e+00; + + + //Gradient on mean value + double gradient; + if(reactionName=="47K(d,p)"){ + sigma = 0.14; + gradient = 0.989; + //gradient = 0.968; + numPeaks_temp = numPeaks; + means_temp = means; + } else if (reactionName=="47K(d,t)"){ + sigma = 0.18; + gradient = 0.986; + numPeaks_temp = numPeaks_dt; + means_temp = means_dt; + } + + string histnamepull = hist->GetName(); + cout << histnamepull << " -> " << histnamepull.substr(14,3) << endl; + int anglenumber = stoi(histnamepull.substr(14,3)); +// int anglenumber = 1; + + + vector<double> means_GradientScaled; + for(int i=0; i<means_temp.size(); i++){ + //if(anglenumber < 120){ + //means_GradientScaled.push_back((means_temp.at(i) * 0.987) - 0.039); + //} else { + // means_GradientScaled.push_back(((means_temp.at(i) * 0.987) - 0.039) * 0.0015*(anglenumber) + 1.154); + //} + means_GradientScaled.push_back(means_temp.at(i) * gradient); + } + + cout << BOLDGREEN + << "FITTING WITH A REAL:OBSERVED Ex GRADIENT GRADIENT OF " + << gradient + //<< " CUSTOM EQUATION!!!!!!" + << RESET + << endl; + + //Build individual peak fit functions + string nameBase = "Peak "; + string function = "gaus";//"[0]*exp(-0.5*pow(((x-[1])/[2]),2))";//REGULAR GAUS + TF1 **allPeaks = new TF1*[numPeaks_temp]; + for(int i=0; i<numPeaks_temp; i++) { + string nameHere = nameBase; + nameHere +=to_string(i); + allPeaks[i] = new TF1(nameHere.c_str(), function.c_str(), minFit, maxFit); + allPeaks[i]->SetLineColor(kBlack); + allPeaks[i]->SetLineStyle(7); + allPeaks[i]->SetNpx(500); + } + + //Build background function + TF1 *bg = new TF1 ("bg","2.*gaus(x-[3],[0],0,[1])*ROOT::Math::normal_cdf([2]*(x-[3]),1,0)", + -15., 10.); + bg->FixParameter(1,fix1); + bg->FixParameter(2,fix2); + bg->FixParameter(3,fix3); + bg->SetLineColor(kGreen); + //bg->SetLineStyle(9); + bg->SetParNames("Background"); + hist->Fit(bg, "RWQBL", "", minFit, -3.0);// -1.5); + bg->SetRange(minFit,maxFit); + + //Build total function + string s_full = f_full_skewgaussbg_string(); + TF1 *full = new TF1("fitAllPeaks", s_full.c_str(), minFit, maxFit); + full->SetLineColor(kRed); + const int numParams = (numPeaks_temp*3); + for(int i=0; i<numPeaks_temp; i++) { + //full->FixParameter((i*3)+1,means_temp.at(i)*gradient); + full->FixParameter((i*3)+1,means_GradientScaled.at(i)); + full->FixParameter((i*3)+2,sigma); + full->SetParameter((i*3)+0,1e1); + full->SetParLimits((i*3)+0,0.0,1e5); + } + full->FixParameter(3*numPeaks_temp+0,bg->GetParameter(0)); + full->FixParameter(3*numPeaks_temp+1,fix1); + full->FixParameter(3*numPeaks_temp+2,fix2); + full->FixParameter(3*numPeaks_temp+3,fix3); + + full->SetNpx(500); + + //Set and fix specific parameters + for (int i=0; i<numPeaks; i++){ + if(means[i]>4.0){ + full->SetParameter((i*3)+2,0.23); + } + if(abs(means[i]-0.279)<0.05){ + full->FixParameter((i*3)+0,0.0); + } + } + + //Fit full function to histogram + hist->Fit(full, "RWQBL", "", minFit, maxFit); + cout << GREEN << "FITTING WITH MAX LIKELYHOOD METHOD" << RESET << endl; + hist->Draw(); + + //Extract fitted variables, assign them to individual fits, and draw them + const Double_t* finalPar = full->GetParameters(); + const Double_t* finalErr = full->GetParErrors(); + for (int i=0; i<numPeaks_temp; i++){ + allPeaks[i]->SetParameters(finalPar[0+(i*3)], + finalPar[1+(i*3)], + finalPar[2+(i*3)]); + } + bg->Draw("SAME"); + full->Draw("SAME"); + + //Make irresolvable doublet clear + if(reactionName=="47K(d,p)"){ + for (int i=0; i<numPeaks_temp; i++){ + if(abs(means[i]-3.83)<0.1){ + cout << BOLDCYAN + << "PEAK #" << i + << " IS FITTING IRRESOLVABLE DOUBLET!" + << RESET + << endl; + allPeaks[i]->SetLineColor(kCyan+1); + } + if(means[i]>4.0){ + cout << BOLDMAGENTA + << "PEAK #" << i + << " IS FITTING UNKNOWN ZONE WITH SIGMA = " + << full->GetParameter((i*3)+2) + << RESET + << endl; + allPeaks[i]->SetLineColor(kMagenta); + } + } + } else if(reactionName=="47K(d,t)"){ + for (int i=0; i<numPeaks_temp; i++){ + if(abs(means[i]-3.37)<0.1){ + cout << BOLDCYAN + << "PEAK #" << i + << " IS FITTING IRRESOLVABLE DOUBLET!" + << RESET + << endl; + allPeaks[i]->SetLineColor(kCyan+1); + } + } + } + + //Draw all peaks + for (int i=0; i<numPeaks_temp; i++){ + allPeaks[i]->Draw("SAME"); + } + + //Write to screen + cout << "===========================" << endl; + cout << "== PEAK =========== AREA ==" << endl; + + //Loop over every peak + vector<vector<double>> allpeaks; + for(int i=0; i<numPeaks_temp; i++){ + //For AREA = HEIGHT * SIGMA * SQRT(2*PI) + double A = (finalPar[0+(i*3)] * finalPar[2+(i*3)] * sqrt(2*pi)) / binWidth; + + double deltaA = A * sqrt( pow( finalErr[0+(i*3)]/finalPar[0+(i*3)] ,2) + + pow( finalErr[2+(i*3)]/finalPar[2+(i*3)] ,2) ); + +// cout << "DELTAAREA = " +// << finalErr[0+(i*3)] << " / " +// << finalPar[0+(i*3)] << " and " +// << finalErr[2+(i*3)] << " / " +// << finalPar[2+(i*3)] << " = " +// << deltaA << endl; + + cout << fixed << setprecision(3) + << " #" << i << " " + << finalPar[(i*3)+1] << "\t" << setprecision(0) + << A << "\t+- " + << deltaA << setprecision(3); + cout << " SQRT: " << sqrt(A) << endl; + + vector<double> onepeak; //energy, area and error for one peak + onepeak.push_back(finalPar[(i*3)+1]); + onepeak.push_back(A); + onepeak.push_back(deltaA); + allpeaks.push_back(onepeak); + cout << "-------------" << endl; + } + + if(removing){ + for(int b=1; b<hist->GetNbinsX()-1; b++){ + //cout << b << " " << hist->GetBinContent(b) << " " << full->Eval(hist->GetBinCenter(b)) << endl; + fitContByBin.push_back(full->Eval(hist->GetBinCenter(b))); + } + } + + return allpeaks; +} + + + + diff --git a/Projects/e793s/macro/Plots_47Kdp.C b/Projects/e793s/macro/Plots_47Kdp.C index 82ab351325bae99004e9fe953734c8420900923e..cfbb1297c7e9fcdfd93b471f2b51a6985c8691bd 100644 --- a/Projects/e793s/macro/Plots_47Kdp.C +++ b/Projects/e793s/macro/Plots_47Kdp.C @@ -3,7 +3,8 @@ string reactionName; /* defined by choice of dp or dt */ #include "GausFit.h" #include "KnownPeakFitter.h" #include "DrawPlots.h" -#include "SimHistFitter.h" +//#include "SimHistFitter.h" +#include "SimHistFitter_EditFromBackup21Dec23.h" #include "CS2_master.h" //#include "CS2.h" @@ -11,6 +12,7 @@ string reactionName; /* defined by choice of dp or dt */ #include "ThreeBodyBreakup.h" #include "ThreeBodyBreakup_FitPhaseSpace.h" //#include "20Oct22_CompareYield.h" +#include "ThesisFigurePlotting.h" void AddGammaLines(TH1F* hist, double particle, double ymax){ @@ -130,7 +132,7 @@ void CS(){ cout << "---- 2.908, p3/2 = CS(2.908, 2, 1, 1.5, 0.05, \"mixed\") " << endl; cout << "---- 2.908, f5/2 = CS(2.908, 2, 3, 2.5, 0.05, \"mixed\") " << endl; cout << "---- 3.254, f5/2 = CS(3.254, 3, 3, 2.5, 0.05, \"\") " << endl; - cout << "---- 3.601, f5/2 = CS(3.601, 3, 3, 2.5, 0.05, \"\") " << endl; + cout << "---- 3.601, f5/2 = CS(3.601, 2, 3, 2.5, 0.05, \"\") " << endl; cout << "---- 3.792 (x-?) " << endl; cout << " FIT TOGETHER = CS(3.830, 2, 3, 2.5, 0.05, \"\")" << endl; cout << "---- 3.868 (2-?) " << endl; @@ -207,19 +209,22 @@ void CompareGammas_ParticleRegions_48K(double binwidth, double minEx, double max void Plots_47Kdp(){ -cout << "test" << endl; - LoadChain47Kdp(); gStyle->SetOptStat("nemMrRi"); -cout << "test" << endl; + + TFile cutIn("cutBG_MGTvE.root"); + TCutG* cutBG_MGTvE = (TCutG*) cutIn.FindObjectAny("cutBG_MGTvE"); + cutBG_MGTvE->SetName("cutBG_MGTvE"); + + tCentre = 270; tRange = 40; //tCentre = 2700; tRange = 400; timegate = "abs(T_MUGAST_VAMOS-" + to_string(tCentre) + ")<" + to_string(tRange); + //timegate = "!cutBG_MGTvE && Mugast.DSSD_T<770 && abs(T_MUGAST_CATS2-230)<200"; det_gate = "Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8"; reactionName = "47K(d,p)"; -cout << "test" << endl; cout << "==============================================" << endl; cout << "=============== (d,p) reaction ===============" << endl; cout << "==============================================" << endl; diff --git a/Projects/e793s/macro/Plots_47Kdt.C b/Projects/e793s/macro/Plots_47Kdt.C index 6a940f252fe5ed5811103e32e574d869ed61ce40..d1d22cd4eecd6db5910f2d3191c086f39387acc6 100644 --- a/Projects/e793s/macro/Plots_47Kdt.C +++ b/Projects/e793s/macro/Plots_47Kdt.C @@ -9,20 +9,21 @@ string reactionName; /* defined by choice of dp or dt */ //#include "ThreeBodyBreakup.h" //#include "ThreeBodyBreakup_FitPhaseSpace.h" +#include "ThesisFigurePlotting.h" void AddGammaLines(TH1F* hist, double particle, double ymax){ -// string base = "sub "; -// -// for(int i=1; i<means.size();i++){ -// string name = base + to_string(means.at(i)); -// TLine *line = new TLine(particle-means.at(i), 0.0, particle-means.at(i), ymax); -// line->SetLineColor(kBlack); line->SetLineStyle(kDotted); -// line->Draw(); -// TText *text = new TText((1.-(means.at(i)/particle))*particle,0.8*ymax,name.c_str()); -// text->SetTextAngle(90); -// //text->SetTextSize(40); -// text->Draw(); -// } + string base = "sub "; + + for(int i=1; i<means_dt.size();i++){ + string name = base + to_string(means_dt.at(i)); + TLine *line = new TLine(particle-means_dt.at(i), 0.0, particle-means_dt.at(i), ymax); + line->SetLineColor(kBlack); line->SetLineStyle(kDotted); + line->Draw(); + TText *text = new TText((1.-(means_dt.at(i)/particle))*particle,0.8*ymax,name.c_str()); + text->SetTextAngle(90); + //text->SetTextSize(40); + text->Draw(); + } } void AddPlacedGammas(TH1F* hist, double ymax){ @@ -74,16 +75,18 @@ void CS_Diagnosis(){ void CS(){ /* Overload function */ cout << "- CS(stateE, stateSp, orb_l, orb_j, nodes) "<< endl; - cout << "---- 0.000, f7/2 = CS(0.000, 4, 3, 3.5, 0.10, \"\") "<< endl; - cout << "---- 0.579, f7/2 = CS(0.579, 4, 3, 3.5, 0.10, \"\") "<< endl; - cout << "---- 0.691, f7/2 = CS(0.691, 4, 3, 3.5, 0.10, \"\") "<< endl; - cout << "---- 1.945, s1/2 = CS(1.945, 1, 0, 0.5, 0.10, \"\") "<< endl; - cout << "---- 2.233, s1/2 = CS(2.233, 1, 0, 0.5, 0.10, \"\") "<< endl; - cout << "---- 2.732, dX/2 = CS(2.732, 2, 2, 1.5, 0.10, \"\") "<< endl; - cout << "---- 3.344 (x-?) " << endl; - cout << " FIT TOGETHER = CS(3.377, 1, 0, 0.5, 0.10, \"\")" << endl; - cout << "---- 3.410 (?-?) " << endl; - cout << "---- 4.297, d5/2 = CS(4.297, 3, 2, 2.5, 0.10, \"\") "<< endl; + cout << "---- 0.000, (2-), p3/2 = CS(0.000, 2, 1, 1.5, 0.10, \"\")"<< endl; + cout << "---- 0.579, 3- " << endl; + cout << " FIT TOGETHER = CS(0.639, 3, 3, 3.5, 0.10, \"\") "<< endl; + cout << "---- 0.691, (4-) " << endl; + cout << "---- 1.738, 3- , f7/2 = CS(1.738, 3, 3, 3.5, 0.10, \"\") "<< endl; + cout << "---- 1.945, 1+ , s1/2 = CS(1.945, 1, 0, 0.5, 0.10, \"\") "<< endl; + cout << "---- 2.233, (0+), s1/2 = CS(2.233, 0, 0, 0.5, 0.10, \"\") "<< endl; + cout << "---- 2.732, (XX), d3/2 = CS(2.732, 2, 2, 1.5, 0.10, \"\") "<< endl; + cout << "---- 3.344, (XX) " << endl; + cout << " FIT TOGETHER = CS(3.377, 1, 0, 0.5, 0.10, \"mixed\") : 2 1.5"<< endl; + cout << "---- 3.410, (XX) " << endl; + cout << "---- 4.297, (3-), d3/2 = CS(4.20, 1, 0, 0.5, 0.10, \"mixed\") : 2 1.5 "<< endl; cout << "NOT SURE OF THE L TRANSFERS OR SPINS OF THESE STATES!!! CORRECT AS I GO ALONG!!!!" << endl; } @@ -124,17 +127,14 @@ void MM_Timing_Comparison(){ void Plots_47Kdt(){ /* Load graphical cut */ - //TFile gcIn("GraphicalCut_22Jun22.root"); - //TCutG* cutTritons = (TCutG*) gcIn.FindObjectAny("cutTritons"); - - //TFile gcIn("cutTritonsWide.root"); - TFile gcIn("cutTritons_26Aug22Long.root"); - //TCutG* cutTritons = (TCutG*) gcIn.FindObjectAny("cutTritonsWide"); + //TFile gcIn("cutTritons_04Jan24.root"); + TFile gcIn("cutTritons_Tight29Jan24.root"); TCutG* cutTritons = (TCutG*) gcIn.FindObjectAny("cutTritons"); cutTritons->SetName("cutTritons"); //TFile gcIn2("cutTime.root"); - TFile gcIn2("cutTimeNew.root"); + //TFile gcIn2("cutTimeNew.root"); + TFile gcIn2("cutTime_05Jan24.root"); TCutG* cutTime = (TCutG*) gcIn2.FindObjectAny("cutTime"); cutTime->SetName("cutTime"); diff --git a/Projects/e793s/phs.sh b/Projects/e793s/phs.sh new file mode 100755 index 0000000000000000000000000000000000000000..f9b835fccb85db957fc9245b2eaf0cdac281a5df --- /dev/null +++ b/Projects/e793s/phs.sh @@ -0,0 +1,97 @@ +#!/bin/bash + +rm RunToTreat_AutoGenerated.txt +rm ./Reaction/Auto.reaction +rm ./Reaction/AutoPhaseSpace.reaction + +XYZTdate="18Oct22" +sp=" " + +rfile_phspc='Reaction/AutoPhaseSpace.reaction' +rfile_npana='Reaction/Auto.reaction' +dfile="./Detector/mugast_"$XYZTdate".detector" + +E=$(echo "scale=3; $2/1000" | bc -l) + +#============================================================================================== +#Boilerplate +echo -e "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" >> ./Reaction/Auto.reaction +echo -e "Beam" >> ./Reaction/Auto.reaction +echo -e "$sp""Particle=""$sp""47K" >> ./Reaction/Auto.reaction +echo -e "$sp""ExcitationEnergy=""$sp""0""$sp""MeV" >> ./Reaction/Auto.reaction +echo -e "$sp""Energy=""$sp""354.75""$sp""MeV" >> ./Reaction/Auto.reaction +echo -e "$sp""SigmaEnergy=""$sp""0.""$sp"" MeV" >> ./Reaction/Auto.reaction +echo -e "$sp""SigmaThetaX=""$sp""0.0""$sp""deg" >> ./Reaction/Auto.reaction +echo -e "$sp""SigmaPhiY=""$sp""0.0""$sp""deg" >> ./Reaction/Auto.reaction +echo -e "$sp""SigmaX=""$sp""2.0""$sp""millimeter" >> ./Reaction/Auto.reaction +echo -e "$sp""SigmaY=""$sp""2.0""$sp""millimeter" >> ./Reaction/Auto.reaction +echo -e "$sp""MeanThetaX=""$sp""0""$sp""deg" >> ./Reaction/Auto.reaction +echo -e "$sp""MeanPhiY=""$sp""0""$sp""deg" >> ./Reaction/Auto.reaction +echo -e "$sp""MeanX=""$sp""-4.1561""$sp""millimeter" >> ./Reaction/Auto.reaction +echo -e "$sp""MeanY=""$sp""+0.4701""$sp""millimeter" >> ./Reaction/Auto.reaction +echo -e "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" >> ./Reaction/Auto.reaction +echo -e "TwoBodyReaction" >> ./Reaction/Auto.reaction +echo -e "$sp""Beam=""$sp""47K" >> ./Reaction/Auto.reaction +echo -e "$sp""Target=""$sp""2H" >> ./Reaction/Auto.reaction +echo -e "$sp""Light=""$sp""1H" >> ./Reaction/Auto.reaction +echo -e "$sp""Heavy=""$sp""48K" >> ./Reaction/Auto.reaction +echo -e "$sp""ExcitationEnergyLight=""$sp""0.0""$sp""MeV" >> ./Reaction/Auto.reaction +echo -e "$sp""ExcitationEnergyHeavy=""$sp""$E""$sp""MeV" >> ./Reaction/Auto.reaction +echo -e "$sp""CrossSectionPath=""$sp""flat.txt""$sp""CSR" >> ./Reaction/Auto.reaction +echo -e "$sp""ShootLight=""$sp""1" >> ./Reaction/Auto.reaction +echo -e "$sp""ShootHeavy=""$sp""0" >> ./Reaction/Auto.reaction + +#============================================================================================== + +echo -e "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" >> ./Reaction/AutoPhaseSpace.reaction +echo -e "Beam" >> ./Reaction/AutoPhaseSpace.reaction +echo -e "$sp""Particle=""$sp""47K" >> ./Reaction/AutoPhaseSpace.reaction +echo -e "$sp""ExcitationEnergy=""$sp""0""$sp""MeV" >> ./Reaction/AutoPhaseSpace.reaction +echo -e "$sp""Energy=""$sp""354.75""$sp""MeV" >> ./Reaction/AutoPhaseSpace.reaction +echo -e "$sp""SigmaEnergy=""$sp""0.""$sp"" MeV" >> ./Reaction/AutoPhaseSpace.reaction +echo -e "$sp""SigmaThetaX=""$sp""0.0""$sp""deg" >> ./Reaction/AutoPhaseSpace.reaction +echo -e "$sp""SigmaPhiY=""$sp""0.0""$sp""deg" >> ./Reaction/AutoPhaseSpace.reaction +echo -e "$sp""SigmaX=""$sp""2.0""$sp""millimeter" >> ./Reaction/AutoPhaseSpace.reaction +echo -e "$sp""SigmaY=""$sp""2.0""$sp""millimeter" >> ./Reaction/AutoPhaseSpace.reaction +echo -e "$sp""MeanThetaX=""$sp""0""$sp""deg" >> ./Reaction/AutoPhaseSpace.reaction +echo -e "$sp""MeanPhiY=""$sp""0""$sp""deg" >> ./Reaction/AutoPhaseSpace.reaction +echo -e "$sp""MeanX=""$sp""-4.1561""$sp""millimeter" >> ./Reaction/AutoPhaseSpace.reaction +echo -e "$sp""MeanY=""$sp""+0.4701""$sp""millimeter" >> ./Reaction/AutoPhaseSpace.reaction +echo -e "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" >> ./Reaction/AutoPhaseSpace.reaction +echo -e "PhaseSpace" >> ./Reaction/AutoPhaseSpace.reaction +echo -e "$sp""Beam=""$sp""47K" >> ./Reaction/AutoPhaseSpace.reaction +echo -e "$sp""Target=""$sp""2H" >> ./Reaction/AutoPhaseSpace.reaction +echo -e "$sp""Daughters=""$sp""1H""$sp""47K""$sp""n" >> ./Reaction/AutoPhaseSpace.reaction +echo -e "$sp""ExcitationEnergies=""$sp""0.0""$sp""$E""$sp""0.0""$sp""MeV" >> ./Reaction/AutoPhaseSpace.reaction +echo -e "$sp""Fermi=""$sp""1" >> ./Reaction/AutoPhaseSpace.reaction + +cd ~/Programs/nptool/Projects/e793s; +cmake ./; +make -j6; + +directory=' /home/charlottepaxman/Programs/nptool/Outputs/Simulation/' +dotroot='.root' + +nameconstruct=$XYZTdate"_"$1"_01Dec23" +echo "Name: $nameconstruct" +echo "Reaction file: $rfile_phspc" +echo "Detector file: $dfile" + +mvnamesim_part="/media/charlottepaxman/Elements1/ActiveWork_47K/Simulation/""$nameconstruct" +sp=" " + +echo "TTreeName" > RunToTreat_AutoGenerated.txt +echo " SimulatedTree" >> RunToTreat_AutoGenerated.txt +echo "RootFileName" >> RunToTreat_AutoGenerated.txt + +for x in {1..64} +do + outname=$nameconstruct"-"$x + npsimulation -D $dfile -E $rfile_phspc -B ./runsimulation.mac --random-seed $RANDOM -O $outname; + mv "../../Outputs/Simulation/""$outname"".root" "$mvnamesim_part""-""$x"".root" + echo "$sp""$mvnamesim_part""-""$x"".root" >> RunToTreat_AutoGenerated.txt +done + +outfile="Sim_"$nameconstruct + +npanalysis --definition PhS -R RunToTreat_AutoGenerated.txt -E $rfile_npana -D $dfile -O $outfile; diff --git a/Projects/e793s/sim.sh b/Projects/e793s/sim.sh index ea06cb9e117e453520867988559921065f2db363..9d330da3d4f615a8698fac26aa37568060e8fa2b 100755 --- a/Projects/e793s/sim.sh +++ b/Projects/e793s/sim.sh @@ -4,9 +4,12 @@ rm RunToTreat_AutoGenerated.txt rm ./Reaction/Auto.reaction XYZTdate="18Oct22" +#XYZTdate="25Jul23" sp=" " +dfile="./Detector/mugast_"$XYZTdate".detector" rfile="./Reaction/Auto.reaction" + E=$(echo "scale=3; $2/1000" | bc -l) #Boilerplate @@ -56,13 +59,10 @@ fi #Define heavy recoil excitation & finish echo -e "$sp""ExcitationEnergyLight=""$sp""0.0""$sp""MeV" >> ./Reaction/Auto.reaction echo -e "$sp""ExcitationEnergyHeavy=""$sp""$E""$sp""MeV" >> ./Reaction/Auto.reaction -echo -e "$sp""CrossSectionPath=""$sp""flat.txt""$sp""CSR" >> ./Reaction/Auto.reaction +echo -e "$sp""LabCrossSectionPath=""$sp""flat.txt""$sp""CSR" >> ./Reaction/Auto.reaction echo -e "$sp""ShootLight=""$sp""1" >> ./Reaction/Auto.reaction echo -e "$sp""ShootHeavy=""$sp""0" >> ./Reaction/Auto.reaction - -#rfile="./Reaction/18Oct22_47Kdp_"$2".reaction" -dfile="./Detector/mugast_"$XYZTdate".detector" #==================================================== cd ~/Programs/nptool/Projects/e793s; @@ -76,12 +76,10 @@ 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 1 2 3 4 5 #6 7 8 9 10 do outname=$nameconstruct"-"$x - npsimulation -D $dfile -E $rfile -B ./runsimulation.mac -O $outname; - #npsimulation -D $dfile -E $rfile -B ./runsimulation_small.mac -O $outname; - #filename=$directory$nameconstruct"-"$x".root" + npsimulation -D $dfile -E $rfile -B ./runsimulation.mac --random-seed $RANDOM -O $outname; done outfile="Sim_"$nameconstruct @@ -90,13 +88,12 @@ echo "TTreeName" >> RunToTreat_AutoGenerated.txt echo " SimulatedTree" >> RunToTreat_AutoGenerated.txt echo "RootFileName" >> RunToTreat_AutoGenerated.txt -for x in 1 #2 3 4 #5 6 7 8 9 10 +for x in 1 2 3 4 5 #6 7 8 9 10 do echo "$sp""$directory""$nameconstruct""-""$x"".root" >> RunToTreat_AutoGenerated.txt done npanalysis --definition Sim -R RunToTreat_AutoGenerated.txt -E $rfile -D $dfile -O $outfile; -#####npanalysis --definition Sim --definition ExcludeThePoor -R RunToTreat_AutoGenerated.txt -E $rfile -D $dfile -O $outfile; mvname="./macro/SolidAngle_HistFiles/SAHF_"$nameconstruct".root" mv SolidAngle_HistFile_New.root $mvname