From 65cbe1faad57566ca6ca7471311b2d5127dc7554 Mon Sep 17 00:00:00 2001 From: Charlie Paxman <cp00474@surrey.ac.uk> Date: Fri, 25 Jun 2021 14:04:57 +0100 Subject: [PATCH] * e793s plotting updates + Far more useful plotting/gating functions --- Projects/e793s/macro/DrawPlots.C | 237 ++-------------------- Projects/e793s/macro/DrawPlots.h | 336 +++++++++++++++++++++++++++++++ 2 files changed, 357 insertions(+), 216 deletions(-) create mode 100644 Projects/e793s/macro/DrawPlots.h diff --git a/Projects/e793s/macro/DrawPlots.C b/Projects/e793s/macro/DrawPlots.C index e076d73bc..38c11d1fe 100755 --- a/Projects/e793s/macro/DrawPlots.C +++ b/Projects/e793s/macro/DrawPlots.C @@ -1,196 +1,11 @@ -//void AddTiStates(double E); -#include "NPReaction.h" -#include <string> -#include <sstream> +#include "DrawPlots.h" -using namespace std; -//TCutG* ETOF=NULL; -//TCutG* EDE=NULL; -TChain* chain=NULL ; +/* USE THIS SPACE TO TEST NEW FEATURES */ +/* ONCE HAPPY, MOVE TO HEADER */ -char cond[1000]; - -//////////////////////////////////////////////////////////////////////////////// -TChain* Chain(std::string TreeName, std::vector<std::string>& file, bool EventList){ - TChain* chain = new TChain(TreeName.c_str()); - unsigned int size =file.size(); - for(unsigned int i = 0 ; i < size ; i++){ - cout << "Adding " << file[i] << endl; - chain->Add(file[i].c_str()); - /* TFile* file = new TFile(file[i].c_str(),"READ"); - if(EventList){ - TEventList* current = (TEventList*) file->Get("GoodEntries"); - if(!EL) - EL=current; - else{ - - EL->Add(current); - } - }*/ - } - //chain->SetEntryListFile("$.root/"); - return chain; -} -//////////////////////////////////////////////////////////////////////////////// -void LoadChainNP(){ - vector<string> files; - files.push_back("../../../Outputs/Analysis/47K_Full_09June_MG3_Target.root"); - chain = Chain("PhysicsTree",files,true); -} - -void plot_kine(NPL::Reaction r, double Ex,Color_t c,int w, int s){ - r.SetExcitation4(Ex); - TGraph* g= r.GetKinematicLine3(); - g->SetLineColor(c) ; - g->SetLineStyle(s) ; - g->SetLineWidth(w) ; - g->Draw("c"); -} - -void plot_state(double Ex,double max,Color_t c,int w, int s){ - TLine* line = new TLine(Ex,0,Ex,max) ; - line->SetLineColor(c) ; - line->SetLineStyle(s) ; - line->SetLineWidth(w) ; - line->Draw(); -} - -void AddTiStates(double E){ - NPL::Reaction Ti("47Ti(d,p)48Ti@362"); - Ti.SetExcitationHeavy(E); - auto g = Ti.GetKinematicLine3(); - g->SetLineWidth(1); - g->SetLineStyle(2); - g->Draw("c"); -} -//////////////////////////////////////////////////////////////////////////////// -//////////////////////////////////////////////////////////////////////////////// -//////////////////////////////////////////////////////////////////////////////// -//////////////////////////////////////////////////////////////////////////////// -//////////////////////////////////////////////////////////////////////////////// -void Draw_1DGamma(){ - TCanvas *cEg = new TCanvas("cEg","cEg",1000,1000); - chain->Draw("AddBack_EDC>>Eg(2500,0,5)","abs(T_MUGAST_VAMOS-2777)<600"); - TH1F* Eg = (TH1F*) gDirectory->Get("Eg"); - Eg->SetTitle("Egamma (using 09Jun Full run)"); - Eg->GetXaxis()->SetTitle("Eg [MeV]"); -} -///////////////////////////////////// -void Draw_1DParticle(){ - TCanvas *cEx = new TCanvas("cEx","cEx",1000,1000); - chain->Draw("Ex>>Ep(220,-1,10)","abs(T_MUGAST_VAMOS-2777)<600"); - TH1F* Ep = (TH1F*) gDirectory->Get("Ep"); - Ep->SetTitle("Ex (using 09Jun Full run)"); - Ep->GetXaxis()->SetTitle("Ex [MeV]"); -} -///////////////////////////////////// -void Draw_2DParticleGamma(){ - TCanvas *cExEg = new TCanvas("cExEg","cExEg",1000,1000); - chain->Draw("AddBack_EDC:Ex>>ExEg(100,-1,5,2500,0,5)","abs(T_MUGAST_VAMOS-2777)<600","colz"); - TH1F* ExEg = (TH1F*) gDirectory->Get("ExEg"); - ExEg->SetTitle("Ex-Egamma (using 09Jun Full run)"); - ExEg->GetXaxis()->SetTitle("Ex [MeV]"); - ExEg->GetYaxis()->SetTitle("Eg [MeV]"); - TLine *XeqY = new TLine(0,0,9,9); - XeqY->SetLineColor(kRed); - XeqY->Draw(); -} -/////////////////////////////////// -void GateOnGamma(double gamma, double width){ - string gating = "abs(T_MUGAST_VAMOS-2777)<600 && abs(AddBack_EDC-" - + to_string(gamma) - + ")<" - + to_string(width); - - TCanvas *cEx_Gate = new TCanvas("cEx_Gate","cEx_Gate",1000,1000); - chain->Draw("Ex>>ExGate(220,-1,10)",gating.c_str(),"colz"); - TH1F* ExGate = (TH1F*) gDirectory->Get("ExGate"); - //ExGate->SetTitle("Ex gated on (using 09Jun Full run)"); - ExGate->GetXaxis()->SetTitle("Ex [MeV]"); - ExGate->GetYaxis()->SetTitle("Counts / 0.05 MeV"); -} -/////////////////////////////////// -void GateOnParticle(double particle, double width){ - string gating = "abs(T_MUGAST_VAMOS-2777)<600 && abs(Ex-" - + to_string(particle) - + ")<" - + to_string(width); - - TCanvas *cEg_Gate = new TCanvas("cEg_Gate","cEg_Gate",1000,1000); - chain->Draw("AddBack_EDC>>EgGate(10000,0,10)",gating.c_str(),"colz"); - TH1F* EgGate = (TH1F*) gDirectory->Get("EgGate"); - //ExGate->SetTitle("Ex gated on (using 09Jun Full run)"); - EgGate->GetXaxis()->SetTitle("Eg [MeV]"); - EgGate->GetYaxis()->SetTitle("Counts / 1 keV"); -} -void CompareExsAt4MeV(){ - TCanvas *cExCompare = new TCanvas("cExCompare","cExCompare",1000,1000); - chain->Draw("AddBack_EDC>>gate3p0(1000,0,10)", - "abs(T_MUGAST_VAMOS-2777)<600 && abs(Ex-3.0)<0.1","same"); - chain->Draw("AddBack_EDC>>gate3p5(1000,0,10)", - "abs(T_MUGAST_VAMOS-2777)<600 && abs(Ex-3.5)<0.1","same"); - chain->Draw("AddBack_EDC>>gate3p9(1000,0,10)", - "abs(T_MUGAST_VAMOS-2777)<600 && abs(Ex-3.9)<0.1","same"); - chain->Draw("AddBack_EDC>>gate4p3(1000,0,10)", - "abs(T_MUGAST_VAMOS-2777)<600 && abs(Ex-4.3)<0.1","same"); - - //cExCompare->Divide(4,1); - - //cExCompare->cd(1); - TH1F* gate3p0 = (TH1F*) gDirectory->Get("gate3p0"); - gate3p0->GetXaxis()->SetTitle("Egamma [MeV]"); - gate3p0->GetYaxis()->SetTitle("Counts"); - gate3p0->SetLineColor(kRed); - - //cExCompare->cd(2); - TH1F* gate3p5 = (TH1F*) gDirectory->Get("gate3p5"); - gate3p5->SetLineColor(kBlue); - - //cExCompare->cd(3); - TH1F* gate3p9 = (TH1F*) gDirectory->Get("gate3p9"); - gate3p9->SetLineColor(kGreen); - - //cExCompare->cd(4); - TH1F* gate4p3 = (TH1F*) gDirectory->Get("gate4p3"); - gate4p3->SetLineColor(kViolet); - gate4p3->GetXaxis()->SetRangeUser(0.,5.); - - cout << " 3.0 - Red\n" - << " 3.5 - Blue\n" - << " 3.9 - Green\n" - << " 4.3 - Violet" - << endl; -} - - -void CompareSimExp(){ - - TCanvas *cSimExp = new TCanvas("cSimExp","cSimExp",1000,1000); - gStyle->SetOptStat(0); - - chain->Draw("Ex>>hexp(70,-1,6)","abs(T_MUGAST_VAMOS-2777)<600",""); - TH1F* hexp = (TH1F*) gDirectory->Get("hexp"); - hexp->SetTitle("Comparing Simulation to Experiment"); - hexp->GetXaxis()->SetTitle("Ex [MeV]"); - hexp->GetYaxis()->SetTitle("Counts / 0.1 MeV"); - hexp->SetLineColor(kRed); - - TFile* simfile = new TFile("../../../Outputs/Analysis/SimTest_18June_30mill.root", "READ"); - TTree* simtree = (TTree*) simfile->FindObjectAny("PhysicsTree"); - - simtree->Draw("Ex>>hsim(70,-1,6)","Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8","same"); - TH1F* hsim = (TH1F*) gDirectory->Get("hsim"); - hsim->SetLineColor(kBlue); - - auto legend = new TLegend(0.7,0.8,0.9,0.9); - legend->AddEntry(hexp,"Experiment","l"); - legend->AddEntry(hsim,"Simulation","l"); - legend->Draw(); -} - -//////////////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////// @@ -200,27 +15,27 @@ void DrawPlots(){ gStyle->SetOptStat("nei"); LoadChainNP(); - NPL::Reaction Kdp("47K(d,p)48K@362"); - NPL::Reaction Kdt("47K(d,t)46K@362"); - NPL::Reaction Kdd("47K(d,d)47K@362"); - NPL::Reaction Kpp("47K(p,p)47K@362"); - NPL::Reaction K12C12C("47K(12C,12C)47K@362"); - NPL::Reaction Tidp("47Ti(d,p)48Ti@362"); - NPL::Reaction Tidt("47Ti(d,t)46Ti@362"); - NPL::Reaction Tidd("47Ti(d,d)47Ti@362"); - NPL::Reaction Ti12C12C("47Ti(12C,12C)47Ti@362"); - - - cout << "==========================================" << endl; - cout << " AVAILABLE FUNCTIONS:: " << endl; + cout << "========== AVAILABLE FUNCTIONS ===========" << endl; + cout << " 2D Matrices " << endl; cout << "\t- Draw_2DParticleGamma() "<< endl; - cout << "\t- Draw_1DGamma() "<< endl; + cout << "\t- Draw_2DGammaGamma() "<< endl; + cout << ""<< endl; + cout << " Ungated histograms " << endl; cout << "\t- Draw_1DParticle() "<< endl; + cout << "\t- Draw_1DGamma() "<< endl; cout << ""<< endl; - cout << "\t- GateOnGamma(gamma, width) "<< endl; - cout << "\t- GateOnParticle(particle, width) "<< endl; + cout << " Gated histograms " << endl; + cout << "\t- GateParticle_SeeGamma(particle, width) "<< endl; + cout << "\t- GateGamma_SeeParticle(gamma, width) "<< endl; + cout << "\t- GateGamma_SeeGamma(gamma, width) "<< endl; cout << ""<< endl; + cout << " Gated histograms w/ background selection" << endl; + cout << "\t- GateParticle_SeeGamma_WithBG(particle, width, bgrnd) "<< endl; + cout << "\t- GateGamma_SeeParticle_WithBG(gamma, width, bgrnd) "<< endl; + cout << "\t- GateGamma_SeeGamma_WithBG(gamma, width, bgrnd) "<< endl; + cout << ""<< endl; + cout << " Specific functions" << endl; cout << "\t- CompareExsAt4MeV() "<< endl; cout << "\t- CompareSimExp() "<< endl; cout << "==========================================" << endl; @@ -228,17 +43,7 @@ void DrawPlots(){ - /*** Functions for drawing ***/ - - //Draw_1DGamma(); - //Draw_2DParticleGamma(); - - //GateOnGamma(0.1426, 0.0015); - //GateOnGamma(1.2680, 0.0193); - //GateOnGamma(1.8382, 0.0062); - //GateOnGamma(3.5179, 0.0109); - - /*****************************/ + /************** OLD STUFF! Make these into functions, eventually ***************/ /* TCanvas *cF = new TCanvas("cF","cF",1000,1000); diff --git a/Projects/e793s/macro/DrawPlots.h b/Projects/e793s/macro/DrawPlots.h new file mode 100644 index 000000000..94ec4cd06 --- /dev/null +++ b/Projects/e793s/macro/DrawPlots.h @@ -0,0 +1,336 @@ +#include "NPReaction.h" +#include <string> +#include <sstream> + +using namespace std; + +TChain* chain=NULL ; +char cond[1000]; + +NPL::Reaction Kdp("47K(d,p)48K@362"); +NPL::Reaction Kdt("47K(d,t)46K@362"); +NPL::Reaction Kdd("47K(d,d)47K@362"); +NPL::Reaction Kpp("47K(p,p)47K@362"); +NPL::Reaction K12C12C("47K(12C,12C)47K@362"); +NPL::Reaction Tidp("47Ti(d,p)48Ti@362"); +NPL::Reaction Tidt("47Ti(d,t)46Ti@362"); +NPL::Reaction Tidd("47Ti(d,d)47Ti@362"); +NPL::Reaction Ti12C12C("47Ti(12C,12C)47Ti@362"); + +//////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////// +TChain* Chain(std::string TreeName, std::vector<std::string>& file, bool EventList){ + TChain* chain = new TChain(TreeName.c_str()); + unsigned int size =file.size(); + for(unsigned int i = 0 ; i < size ; i++){ + cout << "Adding " << file[i] << endl; + chain->Add(file[i].c_str()); + } + return chain; +} +///////////////////////////////////// +void LoadChainNP(){ + vector<string> files; + //files.push_back("../../../Outputs/Analysis/47K_Full_09June_MG3_Target.root"); + files.push_back("../../../Outputs/Analysis/47K_Full_25Jun_TestingAddBack2.root"); + chain = Chain("PhysicsTree",files,true); +} +//////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////// +void plot_kine(NPL::Reaction r, double Ex,Color_t c,int w, int s){ + r.SetExcitation4(Ex); + TGraph* g= r.GetKinematicLine3(); + g->SetLineColor(c) ; + g->SetLineStyle(s) ; + g->SetLineWidth(w) ; + g->Draw("c"); +} +///////////////////////////////////// +void plot_state(double Ex,double max,Color_t c,int w, int s){ + TLine* line = new TLine(Ex,0,Ex,max) ; + line->SetLineColor(c) ; + line->SetLineStyle(s) ; + line->SetLineWidth(w) ; + line->Draw(); +} +///////////////////////////////////// +void AddTiStates(double E){ + NPL::Reaction Ti("47Ti(d,p)48Ti@362"); + Ti.SetExcitationHeavy(E); + auto g = Ti.GetKinematicLine3(); + g->SetLineWidth(1); + g->SetLineStyle(2); + g->Draw("c"); +} +//////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////// +void Draw_1DGamma(){ + TCanvas *cEg = new TCanvas("cEg","cEg",1000,1000); + chain->Draw("AddBack_EDC>>Eg(2500,0,5)","abs(T_MUGAST_VAMOS-2777)<600"); + TH1F* Eg = (TH1F*) gDirectory->Get("Eg"); + Eg->SetTitle("Egamma (using 25Jun Full run)"); + Eg->GetXaxis()->SetTitle("Eg [MeV]"); +} +///////////////////////////////////// +void Draw_1DParticle(){ + TCanvas *cEx = new TCanvas("cEx","cEx",1000,1000); + chain->Draw("Ex>>Ep(220,-1,10)","abs(T_MUGAST_VAMOS-2777)<600"); + TH1F* Ep = (TH1F*) gDirectory->Get("Ep"); + Ep->SetTitle("Ex (using 25Jun Full run)"); + Ep->GetXaxis()->SetTitle("Ex [MeV]"); +} +///////////////////////////////////// +void Draw_2DParticleGamma(){ + TCanvas *cExEg = new TCanvas("cExEg","cExEg",1000,1000); + chain->Draw("AddBack_EDC:Ex>>ExEg(100,-1,5,2500,0,5)", + "abs(T_MUGAST_VAMOS-2777)<600","colz"); + TH1F* ExEg = (TH1F*) gDirectory->Get("ExEg"); + ExEg->SetTitle("Ex-Egamma (using 25Jun Full run)"); + ExEg->GetXaxis()->SetTitle("Ex [MeV]"); + ExEg->GetYaxis()->SetTitle("Eg [MeV]"); + TLine *XeqY = new TLine(0,0,9,9); + XeqY->SetLineColor(kRed); + XeqY->Draw(); +} +///////////////////////////////////// +void Draw_2DGammaGamma(){ + TCanvas *cEgEg = new TCanvas("cEgEg","cEgEg",1000,1000); + chain->Draw("AddBack_EDC:AddBack_EDC2>>EgEg(590,0.05,3,590,0.05,3)","","colz"); + TH1F* EgEg = (TH1F*) gDirectory->Get("EgEg"); + EgEg->SetTitle("Egamma-Egamma (25Jun Full run)"); + EgEg->GetXaxis()->SetTitle("Eg [MeV]"); + EgEg->GetYaxis()->SetTitle("Eg [MeV]"); + //TLine *XeqY = new TLine(0,0,3,3); + //XeqY->SetLineColor(kRed); + //XeqY->Draw(); +} +//////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////// +void GateGamma_SeeParticle(double gamma, double width){ + string gating = "abs(T_MUGAST_VAMOS-2777)<600 && abs(AddBack_EDC-" + + to_string(gamma) + + ")<" + + to_string(width); + + TCanvas *cEx_Gate = new TCanvas("cEx_Gate","cEx_Gate",1000,1000); + chain->Draw("Ex>>ExGate(220,-1,10)",gating.c_str(),"colz"); + TH1F* ExGate = (TH1F*) gDirectory->Get("ExGate"); + //ExGate->SetTitle("Ex gated on (using 09Jun Full run)"); + ExGate->GetXaxis()->SetTitle("Ex [MeV]"); + ExGate->GetYaxis()->SetTitle("Counts / 0.05 MeV"); +} +/////////////////////////////////// +void GateGamma_SeeParticle_WithBG(double gamma, double width, double bg){ + string gating = "abs(T_MUGAST_VAMOS-2777)<600 && abs(AddBack_EDC-" + + to_string(gamma) + + ")<" + + to_string(width); + string bggate = "abs(T_MUGAST_VAMOS-2777)<600 && abs(AddBack_EDC-" + + to_string(bg) + + ")<" + + to_string(width); + + TCanvas *cEx_Gate = new TCanvas("cEx_Gate","cEx_Gate",1000,1000); + chain->Draw("Ex>>ExGate(220,-1,10)",gating.c_str(),""); + TH1F* ExGate = (TH1F*) gDirectory->Get("ExGate"); + ExGate->GetXaxis()->SetTitle("Ex [MeV]"); + ExGate->GetYaxis()->SetTitle("Counts / 0.05 MeV"); + ExGate->SetLineColor(kGreen); + ExGate->SetFillColor(kGreen); + ExGate->SetFillStyle(3154); + + chain->Draw("Ex>>ExBG(220,-1,10)",bggate.c_str(),"same"); + TH1F* ExBG = (TH1F*) gDirectory->Get("ExBG"); + ExBG->SetLineColor(kRed); + ExBG->SetFillColor(kRed); + ExBG->SetFillStyle(3345); +} +/////////////////////////////////// +void GateParticle_SeeGamma(double particle, double width){ + string gating = "abs(T_MUGAST_VAMOS-2777)<600 && abs(Ex-" + + to_string(particle) + + ")<" + + to_string(width); + + TCanvas *cEg_Gate = new TCanvas("cEg_Gate","cEg_Gate",1000,1000); + chain->Draw("AddBack_EDC>>EgGate(1000,0,10)",gating.c_str(),"colz"); + TH1F* EgGate = (TH1F*) gDirectory->Get("EgGate"); + //ExGate->SetTitle("Ex gated on (using 09Jun Full run)"); + EgGate->GetXaxis()->SetTitle("Eg [MeV]"); + EgGate->GetYaxis()->SetTitle("Counts / 10 keV"); +} +/////////////////////////////////// +void GateParticle_SeeGamma_WithBG(double particle, double width, double bg){ + string gating = "abs(T_MUGAST_VAMOS-2777)<600 && abs(Ex-" + + to_string(particle) + + ")<" + + to_string(width); + string bggate = "abs(T_MUGAST_VAMOS-2777)<600 && abs(Ex-" + + to_string(bg) + + ")<" + + to_string(width); + + TCanvas *cEg_Gate = new TCanvas("cEg_Gate","cEg_Gate",1000,1000); + chain->Draw("AddBack_EDC>>EgGate(1000,0,10)",gating.c_str(),""); + TH1F* EgGate = (TH1F*) gDirectory->Get("EgGate"); + EgGate->GetXaxis()->SetTitle("Eg [MeV]"); + EgGate->GetYaxis()->SetTitle("Counts / 10 keV"); + EgGate->SetLineColor(kGreen); + EgGate->SetFillColor(kGreen); + EgGate->SetFillStyle(3154); + + chain->Draw("AddBack_EDC>>EgBG(1000,0,10)",bggate.c_str(),"same"); + TH1F* EgBG = (TH1F*) gDirectory->Get("EgBG"); + EgBG->SetLineColor(kRed); + EgBG->SetFillColor(kRed); + EgBG->SetFillStyle(3345); +} +/////////////////////////////////// +void GateGamma_SeeGamma(double gamma, double width){ + string gating = "abs(AddBack_EDC2-" + + to_string(gamma) + + ")<" + + to_string(width); + + TCanvas *cEx_Gate = new TCanvas("cggGate","cggGate",1000,1000); + chain->Draw("AddBack_EDC>>ggGate(590,0.05,3)",gating.c_str(),""); + TH1F* ggGate = (TH1F*) gDirectory->Get("ggGate"); + //ExGate->SetTitle("Ex gated on (using 09Jun Full run)"); + ggGate->GetXaxis()->SetTitle("Eg [MeV]"); + ggGate->GetYaxis()->SetTitle("Counts / 0.05 MeV"); +} +/////////////////////////////////// +void GateGamma_SeeGamma_WithBG(double gamma, double width, double bg){ + string gating = "abs(AddBack_EDC2-" + + to_string(gamma) + + ")<" + + to_string(width); + string bggate = "abs(AddBack_EDC2-" + + to_string(bg) + + ")<" + + to_string(width); + + TCanvas *cggGate = new TCanvas("cggGate","cggGate",1000,1000); + chain->Draw("AddBack_EDC>>ggGate(590,0.05,3)",gating.c_str(),""); + TH1F* ggGate = (TH1F*) gDirectory->Get("ggGate"); + ggGate->GetXaxis()->SetTitle("Eg [MeV]"); + ggGate->GetYaxis()->SetTitle("Counts / 0.05 MeV"); + ggGate->SetLineColor(kGreen); + ggGate->SetFillColor(kGreen); + ggGate->SetFillStyle(3154); + + chain->Draw("AddBack_EDC>>ggBG(590,0.05,3)",bggate.c_str(),"same"); + TH1F* ggBG = (TH1F*) gDirectory->Get("ggBG"); + ggBG->SetLineColor(kRed); + ggBG->SetFillColor(kRed); + ggBG->SetFillStyle(3345); +} +//////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////// +void CompareExsAt4MeV(){ + TCanvas *cExCompare = new TCanvas("cExCompare","cExCompare",1000,1000); + chain->Draw("AddBack_EDC>>gate3p0(1000,0,10)", + "abs(T_MUGAST_VAMOS-2777)<600 && abs(Ex-3.0)<0.1","same"); + chain->Draw("AddBack_EDC>>gate3p5(1000,0,10)", + "abs(T_MUGAST_VAMOS-2777)<600 && abs(Ex-3.5)<0.1","same"); + chain->Draw("AddBack_EDC>>gate3p9(1000,0,10)", + "abs(T_MUGAST_VAMOS-2777)<600 && abs(Ex-3.9)<0.1","same"); + chain->Draw("AddBack_EDC>>gate4p3(1000,0,10)", + "abs(T_MUGAST_VAMOS-2777)<600 && abs(Ex-4.3)<0.1","same"); + + TH1F* gate3p0 = (TH1F*) gDirectory->Get("gate3p0"); + gate3p0->GetXaxis()->SetTitle("Egamma [MeV]"); + gate3p0->GetYaxis()->SetTitle("Counts"); + gate3p0->SetLineColor(kRed); + + TH1F* gate3p5 = (TH1F*) gDirectory->Get("gate3p5"); + gate3p5->SetLineColor(kBlue); + + TH1F* gate3p9 = (TH1F*) gDirectory->Get("gate3p9"); + gate3p9->SetLineColor(kGreen); + + TH1F* gate4p3 = (TH1F*) gDirectory->Get("gate4p3"); + gate4p3->SetLineColor(kViolet); + gate4p3->GetXaxis()->SetRangeUser(0.,5.); + + cout << " 3.0 - Red\n" + << " 3.5 - Blue\n" + << " 3.9 - Green\n" + << " 4.3 - Violet" + << endl; +} +/////////////////////////////////// +void CompareSimExp(){ + TCanvas *cSimExp = new TCanvas("cSimExp","cSimExp",1000,1000); + gStyle->SetOptStat(0); + + chain->Draw("Ex>>hexp(70,-1,6)","abs(T_MUGAST_VAMOS-2777)<600",""); + TH1F* hexp = (TH1F*) gDirectory->Get("hexp"); + hexp->SetTitle("Comparing Simulation to Experiment"); + hexp->GetXaxis()->SetTitle("Ex [MeV]"); + hexp->GetYaxis()->SetTitle("Counts / 0.1 MeV"); + hexp->SetLineColor(kRed); + + TFile* simfile = new TFile("../../../Outputs/Analysis/SimTest_Jun22_TWOFNR_p32.root", + "READ"); + TTree* simtree = (TTree*) simfile->FindObjectAny("PhysicsTree"); + simtree->Draw("Ex>>hsimMGp32(70,-1,6)", + "Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8","same"); + TH1F* hsimMGp32 = (TH1F*) gDirectory->Get("hsimMGp32"); + hsimMGp32->SetLineColor(kBlue); + //simtree->Draw("Ex>>hsimALLp32(70,-1,6)","","same"); + // TH1F* hsimALLp32 = (TH1F*) gDirectory->Get("hsimALLp32"); + // hsimALLp32->SetLineColor(kBlue); + // hsimALLp32->SetLineStyle(kDashed); + + TFile* simfile3 = new TFile("../../../Outputs/Analysis/SimTest_Jun22_TWOFNR_p12.root", + "READ"); + TTree* simtree3 = (TTree*) simfile3->FindObjectAny("PhysicsTree"); + simtree3->Draw("Ex>>hsimMGp12(70,-1,6)", + "Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8","same"); + TH1F* hsimMGp12 = (TH1F*) gDirectory->Get("hsimMGp12"); + hsimMGp12->SetLineColor(kOrange); + //simtree3->Draw("Ex>>hsimALLp12(70,-1,6)","","same"); + // TH1F* hsimALLp12 = (TH1F*) gDirectory->Get("hsimALLp12"); + // hsimALLp12->SetLineColor(kOrange); + // hsimALLp12->SetLineStyle(kDashed); + + TFile* simfile2 = new TFile("../../../Outputs/Analysis/SimTest_Jun22_TWOFNR_f72.root", + "READ"); + TTree* simtree2 = (TTree*) simfile2->FindObjectAny("PhysicsTree"); + simtree2->Draw("Ex>>hsimMGf72(70,-1,6)", + "Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8","same"); + TH1F* hsimMGf72 = (TH1F*) gDirectory->Get("hsimMGf72"); + hsimMGf72->SetLineColor(kGreen); + //simtree2->Draw("Ex>>hsimALLf72(70,-1,6)","","same"); + // TH1F* hsimALLf72 = (TH1F*) gDirectory->Get("hsimALLf72"); + // hsimALLf72->SetLineColor(kGreen); + // hsimALLf72->SetLineStyle(kDashed); + + TFile* simfile4 = new TFile("../../../Outputs/Analysis/SimTest_Jun22_TWOFNR_f52.root", + "READ"); + TTree* simtree4 = (TTree*) simfile4->FindObjectAny("PhysicsTree"); + simtree4->Draw("Ex>>hsimMGf52(70,-1,6)", + "Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8","same"); + TH1F* hsimMGf52 = (TH1F*) gDirectory->Get("hsimMGf52"); + hsimMGf52->SetLineColor(kViolet); + //simtree4->Draw("Ex>>hsimALLf72(70,-1,6)","","same"); + // TH1F* hsimALLf52 = (TH1F*) gDirectory->Get("hsimALLf52"); + // hsimALLf52->SetLineColor(kViolet); + // hsimALLf52->SetLineStyle(kDashed); + + + auto legend = new TLegend(0.7,0.8,0.9,0.9); + legend->AddEntry(hexp, "Experiment","l"); + legend->AddEntry(hsimMGp32, "Simulation p3/2, MG","l"); + legend->AddEntry(hsimMGp12, "Simulation p1/2, MG","l"); + //legend->AddEntry(hsimALLp32, "Simulation p3/2, MG+MM","l"); + //legend->AddEntry(hsimALLp12, "Simulation p1/2, MG+MM","l"); + legend->AddEntry(hsimMGf72, "Simulation f7/2, MG","l"); + legend->AddEntry(hsimMGf52, "Simulation f5/2, MG","l"); + //legend->AddEntry(hsimALLf72, "Simulation f7/2, MG+MM","l"); + //legend->AddEntry(hsimALLf52, "Simulation f5/2, MG+MM","l"); + legend->Draw(); +} +//////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////// -- GitLab