#include void plotdedx(char* estim="median", int energymin=500, int energymax=5000, int pmin=900, int pmax=1100, int etamin=0, int etamax=25, int hitmin=3, char* flatvar="Pt", int mode=1, double smear=0, bool anigif=false) { define_style(); gROOT->SetStyle("BABAR"); // allhits("dedx",energymin,energymax,flatvar); // allhits("N",energymin,energymax,flatvar); // allhitsSingle("dedx","mu",energymin,energymax,flatvar); // allhitsSingle("N","mu",energymin,energymax,flatvar); // together(1,estim,energymin,energymax,pmin,pmax,etamin,etamax,hitmin,flatvar,mode,smear); // together(2,estim,energymin,energymax,pmin,pmax,etamin,etamax,hitmin,flatvar,mode,smear); // final(estim,pmin,pmax,etamin,etamax,hitmin,mode,smear); // compareSmear(estim,"mu",energymin,energymax,mode); // together(3,estim,energymin,energymax,0,0,etamin,etamax,hitmin,flatvar,mode,smear,anigif); together(4,estim,energymin,energymax,0,0,etamin,etamax,hitmin,flatvar,mode,smear); /* together(5,estim,energymin,energymax,0,0,etamin,etamax,hitmin,flatvar,mode,smear); together(6,estim,energymin,energymax,0,0,etamin,etamax,hitmin,flatvar,mode,smear); together(7,estim,energymin,energymax,0,0,etamin,etamax,hitmin,flatvar,mode,smear); together(8,estim,energymin,energymax,pmin,pmax,etamin,etamax,hitmin,flatvar,mode,smear,anigif); // NN vs NN profileDeltaP("dp","p",energymin,energymax,etamin,etamax,hitmin,flatvar,mode); profileDeltaP("dp2","p",energymin,energymax,etamin,etamax,hitmin,flatvar,mode); profileDeltaP("dp","pt",energymin,energymax,etamin,etamax,hitmin,flatvar,mode); profileDeltaP("dp2","pt",energymin,energymax,etamin,etamax,hitmin,flatvar,mode); */ // plotSimple("n",3000,5000,"Pt",1); } void allhits(char* var="dedx", int energymin=1000, int energymax=5000, char* flatvar="Pt") { gROOT->Reset(); gStyle->SetOptStat(0000); gStyle->SetCanvasBorderMode(-1); gStyle->SetCanvasBorderSize(1); gStyle->SetCanvasColor(10); TCanvas* myCanvas = new TCanvas("myCanvas","all hits",200,10,900,800); myCanvas->SetFillColor(10); TPad* pad1 = new TPad("pad1","PXB",0.05,0.69,0.45,0.98,21); TPad* pad2 = new TPad("pad2","PXE",0.50,0.69,0.95,0.98,21); TPad* pad3 = new TPad("pad3","TIB",0.05,0.37,0.45,0.66,21); TPad* pad4 = new TPad("pad4","TID",0.50,0.37,0.95,0.66,21); TPad* pad5 = new TPad("pad5","TOB",0.05,0.05,0.45,0.33,21); TPad* pad6 = new TPad("pad6","TEC",0.50,0.05,0.95,0.33,21); // substitute with myCanvas->Divide(2,3) and find out how to set names to pads pad1->Draw(); pad2->Draw(); pad3->Draw(); pad4->Draw(); pad5->Draw(); pad6->Draw(); pad1->cd(); int mode=0; plotHistoHits(var,"pxb","mu",energymin,energymax,flatvar,mode,1); plotHistoHits(var,"pxb","pi",energymin,energymax,flatvar,mode,2,"same"); plotHistoHits(var,"pxb","k",energymin,energymax,flatvar,mode,3,"same"); plotHistoHits(var,"pxb","proton",energymin,energymax,flatvar,mode,4,"same"); plotHistoHits(var,"pxb","e",energymin,energymax,flatvar,mode,6,"same"); pad2->cd(); plotHistoHits(var,"pxe","mu",energymin,energymax,flatvar,mode,1); plotHistoHits(var,"pxe","pi",energymin,energymax,flatvar,mode,2,"same"); plotHistoHits(var,"pxe","k",energymin,energymax,flatvar,mode,3,"same"); plotHistoHits(var,"pxe","proton",energymin,energymax,flatvar,mode,4,"same"); plotHistoHits(var,"pxe","e",energymin,energymax,flatvar,mode,6,"same"); pad3->cd(); plotHistoHits(var,"tib","mu",energymin,energymax,flatvar,mode,1); plotHistoHits(var,"tib","pi",energymin,energymax,flatvar,mode,2,"same"); plotHistoHits(var,"tib","k",energymin,energymax,flatvar,mode,3,"same"); plotHistoHits(var,"tib","proton",energymin,energymax,flatvar,mode,4,"same"); plotHistoHits(var,"tib","e",energymin,energymax,flatvar,mode,6,"same"); pad4->cd(); plotHistoHits(var,"tid","mu",energymin,energymax,flatvar,mode,1); plotHistoHits(var,"tid","pi",energymin,energymax,flatvar,mode,2,"same"); plotHistoHits(var,"tid","k",energymin,energymax,flatvar,mode,3,"same"); plotHistoHits(var,"tid","proton",energymin,energymax,flatvar,mode,4,"same"); plotHistoHits(var,"tid","e",energymin,energymax,flatvar,mode,6,"same"); pad5->cd(); plotHistoHits(var,"tob","mu",energymin,energymax,flatvar,mode,1); plotHistoHits(var,"tob","pi",energymin,energymax,flatvar,mode,2,"same"); plotHistoHits(var,"tob","k",energymin,energymax,flatvar,mode,3,"same"); plotHistoHits(var,"tob","proton",energymin,energymax,flatvar,mode,4,"same"); plotHistoHits(var,"tob","e",energymin,energymax,flatvar,mode,6,"same"); pad6->cd(); plotHistoHits(var,"tec","mu",energymin,energymax,flatvar,mode,1); plotHistoHits(var,"tec","pi",energymin,energymax,flatvar,mode,2,"same"); plotHistoHits(var,"tec","k",energymin,energymax,flatvar,mode,3,"same"); plotHistoHits(var,"tec","proton",energymin,energymax,flatvar,mode,4,"same"); plotHistoHits(var,"tec","e",energymin,energymax,flatvar,mode,6,"same"); TString gif("allhits_"); gif+=var; gif+="_"; gif+=energymin; gif+="_"; gif+=energymax; gif+="mev"; if (flatvar=="E") gif+="_flatE"; if (mode != 0) { gif+="_mode"; gif+=mode; } gif+=".png"; myCanvas->Print(gif); } void allhitsSingle(char* var="dedx", char* type="mu", int energymin=500, int energymax=5000, char* flatvar="Pt") { gROOT->Reset(); gStyle->SetOptStat(0000); define_style(); gROOT->SetStyle("BABAR"); TCanvas* myCanvas = new TCanvas("myCanvas","all hits",200,10,900,800); myCanvas->SetFillColor(10); // myCanvas->Divide(2,3); myCanvas->Divide(2,2); int mode = 1; /* myCanvas->cd(1); plotHistoHits(var,"pxb",type,energymin,energymax,flatvar,mode,1); TPaveLabel *legend1 = new TPaveLabel(0.75,0.80,0.95,0.95,"PXB","NDC"); legend1->Draw(); myCanvas->Update(); myCanvas->cd(2); plotHistoHits(var,"pxe",type,energymin,energymax,flatvar,mode,1); TPaveLabel *legend2 = new TPaveLabel(0.75,0.80,0.95,0.95,"PXF","NDC"); legend2->Draw(); myCanvas->Update(); */ // myCanvas->cd(3); myCanvas->cd(1); plotHistoHits(var,"tib",type,energymin,energymax,flatvar,mode,1); TPaveLabel *legend3 = new TPaveLabel(0.75,0.80,0.95,0.95,"TIB","NDC"); legend3->Draw(); myCanvas->Update(); // myCanvas->cd(4); myCanvas->cd(2); plotHistoHits(var,"tid",type,energymin,energymax,flatvar,mode,1); TPaveLabel *legend4 = new TPaveLabel(0.75,0.80,0.95,0.95,"TID","NDC"); legend4->Draw(); myCanvas->Update(); // myCanvas->cd(5); myCanvas->cd(3); plotHistoHits(var,"tob",type,energymin,energymax,flatvar,mode,1); TPaveLabel *legend5 = new TPaveLabel(0.75,0.80,0.95,0.95,"TOB","NDC"); legend5->Draw(); myCanvas->Update(); // myCanvas->cd(6); myCanvas->cd(4); plotHistoHits(var,"tec",type,energymin,energymax,flatvar,mode,1); TPaveLabel *legend6 = new TPaveLabel(0.75,0.80,0.95,0.95,"TEC","NDC"); legend6->Draw(); myCanvas->Update(); TString gif("allhits_"); gif+=type; gif+="_"; gif+=var; gif+="_"; gif+=energymin; gif+="_"; gif+=energymax; gif+="mev"; if (flatvar=="E") gif+="_flatE"; gif+=".png"; myCanvas->Print(gif); } void together(int switch=1,char* var="mean", int energymin=500, int energymax=5000, int pmin=900, int pmax=1100, int etamin=0, int etamax=25, int hitmin=0, char* flatvar="Pt", int mode=1, double smear=0, bool anigif=false) { gROOT->Reset(); //gStyle->SetOptStat(0000); gStyle->SetCanvasBorderMode(-1); gStyle->SetCanvasBorderSize(1); gStyle->SetCanvasColor(10); TString name1 = filename("mu",energymin,energymax,flatvar,mode,smear); TString name2 = filename("pi",energymin,energymax,flatvar,mode,smear); TString name3 = filename("k",energymin,energymax,flatvar,mode,smear); TString name4 = filename("proton",energymin,energymax,flatvar,mode,smear); TFile f1(name1); TFile f2(name2); TFile f3(name3); TFile f4(name4); TTree *tree1 = (TTree*)f1.Get("T"); TTree *tree2 = (TTree*)f2.Get("T"); TTree *tree3 = (TTree*)f3.Get("T"); TTree *tree4 = (TTree*)f4.Get("T"); bool log=false; if (switch == 2) log=true; bool norma=false; if (switch == 1) norma=true; Int_t nbinsx=50; Double_t xlow,xup; if (log) {xlow=0.;xup=3.;} else {xlow=0.8;xup=8.;} if (var=="nn1" || var=="nn2") {xlow=-1.0;xup=2.0;} TH1F *h1 = new TH1F("hmu", "dE/dx [MeV/cm]", nbinsx, xlow, xup); TH1F *h2 = new TH1F("hpi", "dE/dx [MeV/cm]", nbinsx, xlow, xup); TH1F *h3 = new TH1F("hk", "dE/dx [MeV/cm]", nbinsx, xlow, xup); TH1F *h4 = new TH1F("hp", "dE/dx [MeV/cm]", nbinsx, xlow, xup); TString gif(var); gif+="_"; gif+=energymin; gif+="_"; gif+=energymax; gif+="mev"; if (flatvar=="E") gif+="_flatE"; TCanvas* myCanvas = new TCanvas("myCanvas",var,200,10,900,800); myCanvas->SetFillColor(10); if (switch == 1) { // plot of the 4 particle types, in the same canvas with lin scale myCanvas->cd(); fillHisto(h1,var,tree1,pmin,pmax,etamin,etamax,hitmin,norma,log,"muons",mode,smear,2); fillHisto(h2,var,tree2,pmin,pmax,etamin,etamax,hitmin,norma,log,"pions",mode,smear,1,"same"); fillHisto(h3,var,tree3,pmin,pmax,etamin,etamax,hitmin,norma,log,"kaons",mode,smear,3,"same"); fillHisto(h4,var,tree4,pmin,pmax,etamin,etamax,hitmin,norma,log,"protons",mode,smear,4,"same"); // legend: TLegend *legend = new TLegend(0.75,0.75,1.00,0.95,"","NDC"); legend->AddEntry(h1,"muons","l"); legend->AddEntry(h2,"pions","l"); legend->AddEntry(h3,"kaons","l"); legend->AddEntry(h4,"protons","l"); legend->Draw(); myCanvas->Update(); // png: gif+="_P"; gif+=pmin; gif+="_"; gif+=pmax; if (etamin>0 || etamax<25) { gif+="_eta"; gif+=etamin; gif+="_"; gif+=etamax; } if (hitmin>5) { gif+="_hit"; gif+=hitmin; } if (mode != 0) { gif+="_mode"; gif+=mode; } if (smear != 0) { gif+="_smear"; gif+=(int)smear; } gif+=".png"; } else if (switch == 2) { // plot of the 4 particle types, in different sub-canvases with log scale if (var!="nn1" && var!="nn2") { myCanvas->Divide(2,2); myCanvas->cd(1); fillHisto(h1,var,tree1,pmin,pmax,etamin,etamax,hitmin,norma,log,"muons",mode,smear); myCanvas->cd(2); fillHisto(h2,var,tree2,pmin,pmax,etamin,etamax,hitmin,norma,log,"pions",mode,smear); myCanvas->cd(3); fillHisto(h3,var,tree3,pmin,pmax,etamin,etamax,hitmin,norma,log,"kaons",mode,smear); myCanvas->cd(4); fillHisto(h4,var,tree4,pmin,pmax,etamin,etamax,hitmin,norma,log,"protons",mode,smear); gif+="_P"; gif+=pmin; gif+="_"; gif+=pmax; if (etamin>0 || etamax<25) { gif+="_eta"; gif+=etamin; gif+="_"; gif+=etamax; } if (hitmin>5) { gif+="_hit"; gif+=hitmin; } if (mode != 0) { gif+="_mode"; gif+=mode; } if (smear != 0) { gif+="_smear"; gif+=(int)smear; } gif+="_log.png"; } else { // NN cout << "NO LOG PLOT FOR NN!" << endl; } } else if (switch == 3) { // scatter plot dE/dx vs P myCanvas->cd(); if (etamin>0 || etamax<25) { gif+="_eta"; gif+=etamin; gif+="_"; gif+=etamax; } if (hitmin>5) { gif+="_hit"; gif+=hitmin; } gif+="_vsP"; if (anigif) { TString anigif1(gif); TString anigif2(gif); TString anigif3(gif); TString anigif4(gif); TString command("~/scratch0/gifsicle/bin/gifsicle --delay=100 --loop=5 "); command+=gif; command+="_*.png > "; command+=gif; if (mode != 0) { command+="_mode"; command+=mode; } if (smear != 0) { command+="_smear"; command+=(int)smear; } command+="_anim.gif"; } gif+=".gif"; if (anigif) { anigif1+="_1.gif"; anigif2+="_2.gif"; anigif3+="_3.gif"; anigif4+="_4.gif"; } plotScatter(var,tree1,hitmin,2); if (anigif) myCanvas->Print(anigif1); plotScatter(var,tree2,hitmin,1,"same"); if (anigif) myCanvas->Print(anigif2); plotScatter(var,tree3,hitmin,3,"same"); if (anigif) myCanvas->Print(anigif3); plotScatter(var,tree4,hitmin,4,"same"); if (anigif) myCanvas->Print(anigif4); //make an animated gif file using gifsicle. //for information about gifsicle, see http://www.lcdf.org/~eddietwo/gifsicle/ if (anigif) { cout << "Executing command: " << command << endl; gSystem->Exec(command); gSystem->Exec("rm -f *_vsP_1.gif"); gSystem->Exec("rm -f *_vsP_2.gif"); gSystem->Exec("rm -f *_vsP_3.gif"); gSystem->Exec("rm -f *_vsP_4.gif"); } } else if (switch == 4 || switch == 5 || switch==6 || switch==7) { // profile histo dE/dx vs P/Pt/eta/nhits, and separations vs P/Pt/eta/nhits Int_t nbinsx=30; Double_t xlow; Double_t xup; char* title; char* abscissa; if (switch == 4) {title="dE/dx [MeV/cm] vs P [GeV]"; abscissa="p"; xlow=0.8; xup=4.8;} if (switch == 5) {title="dE/dx [MeV/cm] vs Pt [GeV]"; abscissa="pt"; xlow=0.8; xup=4.8;} if (switch == 6) {title="dE/dx [MeV/cm] vs eta"; abscissa="eta"; xlow=0.; xup=1.3;} if (switch == 7) {title="dE/dx [MeV/cm] vs hits"; abscissa="n"; xlow=5; xup=20; nbins=15;} TProfile *prof1 = new TProfile("pmu", title, nbinsx, xlow, xup); TProfile *prof2 = new TProfile("ppi", title, nbinsx, xlow, xup); TProfile *prof3 = new TProfile("pk", title, nbinsx, xlow, xup); TProfile *prof4 = new TProfile("pp", title, nbinsx, xlow, xup); // myCanvas->Divide(1,2); myCanvas->cd(1); fillProfile(prof4,var,abscissa,tree4,etamin,etamax,hitmin,"protons",4,""); fillProfile(prof3,var,abscissa,tree3,etamin,etamax,hitmin,"kaons",3,"same"); fillProfile(prof1,var,abscissa,tree1,etamin,etamax,hitmin,"muons",2,"same"); fillProfile(prof2,var,abscissa,tree2,etamin,etamax,hitmin,"pions",1,"same"); // legend: TLegend *legend = new TLegend(0.70,0.65,0.95,0.95,"","NDC"); legend->AddEntry(prof1,"muons","l"); legend->AddEntry(prof2,"pions","l"); legend->AddEntry(prof3,"kaons","l"); legend->AddEntry(prof4,"protons","l"); legend->Draw(); myCanvas->Update(); // particle separation: TGraph* gr1 = new TGraph(); /* myCanvas->cd(2); separation(gr1,prof4,prof2,nbinsx,xlow,xup,"proton/pion separation (RMS)",abscissa,23); */ if (etamin>0 || etamax<25) { gif+="_eta"; gif+=etamin; gif+="_"; gif+=etamax; } if (hitmin>5) { gif+="_hit"; gif+=hitmin; } if (mode != 0) { gif+="_mode"; gif+=mode; } if (smear != 0) { gif+="_smear"; gif+=(int)smear; } if (abscissa=="p") {gif+="_prof.png";} else if (abscissa=="pt") {gif+="_profPt.png";} else if (abscissa=="eta") {gif+="_profEta.png";} else if (abscissa=="n") {gif+="_profNhits.png";} else {gif+="else";} delete prof1; delete prof2; delete prof3; delete prof4; delete gr1; } else if (switch == 8) { // scatter plot NN1 vs NN2 myCanvas->cd(); TString gif("nn1_nn2"); gif+="_"; gif+=energymin; gif+="_"; gif+=energymax; gif+="mev"; gif+="_P"; gif+=pmin; gif+="_"; gif+=pmax; if (etamin>0 || etamax<25) { gif+="_eta"; gif+=etamin; gif+="_"; gif+=etamax; } if (hitmin>5) { gif+="_hit"; gif+=hitmin; } if (anigif) { TString anigif1(gif); TString anigif2(gif); TString anigif3(gif); TString anigif4(gif); TString command("~/scratch0/gifsicle/bin/gifsicle --delay=100 --loop=5 "); command+=gif; command+="_*.gif > "; command+=gif; if (mode != 0) { command+="_mode"; command+=mode; } if (smear != 0) { command+="_smear"; command+=(int)smear; } command+="_anim.gif"; } gif+=".gif"; if (anigif) { anigif1+="_1.gif"; anigif2+="_2.gif"; anigif3+="_3.gif"; anigif4+="_4.gif"; } nnScatter(tree1,pmin,pmax,etamin,etamax,hitmin,1); if (anigif) myCanvas->Print(anigif1); nnScatter(tree2,pmin,pmax,etamin,etamax,hitmin,2,"same"); if (anigif) myCanvas->Print(anigif2); nnScatter(tree3,pmin,pmax,etamin,etamax,hitmin,3,"same"); if (anigif) myCanvas->Print(anigif3); nnScatter(tree4,pmin,pmax,etamin,etamax,hitmin,4,"same"); if (anigif) myCanvas->Print(anigif4); //make an animated gif file using gifsicle. //for information about gifsicle, see http://www.lcdf.org/~eddietwo/gifsicle/ if (anigif) { cout << "Executing command: " << command << endl; gSystem->Exec(command); gSystem->Exec("rm -f *_vsP_1.gif"); gSystem->Exec("rm -f *_vsP_2.gif"); gSystem->Exec("rm -f *_vsP_3.gif"); gSystem->Exec("rm -f *_vsP_4.gif"); } } else { cout << "Mode " << switch << " does not exist." << endl; } myCanvas->Print(gif); if (anigif) cout << "See animated figure with the command ~/scratch0/gifsicle/bin/gifview -a" << endl; } void compareSmear(char* var="harm2", char* particle="mu", int energymin=500, int energymax=5000, int mode=1) { gROOT->Reset(); gStyle->SetOptStat(000000); gStyle->SetCanvasBorderMode(-1); gStyle->SetCanvasBorderSize(1); gStyle->SetCanvasColor(10); // no more grey bkg TString name1 = filename(particle,energymin,energymax,"Pt",mode,0.); // TString name2 = filename(particle,energymin,energymax,"Pt",mode,2.); TString name2 = filename(particle,energymin,energymax,"Pt",mode,6.4); // TString name3 = filename(particle,energymin,energymax,"Pt",mode,-0.4); TString name3 = filename(particle,energymin,energymax,"Pt",mode,-3.); TFile f1(name1); TFile f2(name2); TFile f3(name3); // TFile f4(name4); TTree *tree1 = (TTree*)f1.Get("T"); TTree *tree2 = (TTree*)f2.Get("T"); TTree *tree3 = (TTree*)f3.Get("T"); // TTree *tree4 = (TTree*)f4.Get("T"); double variable; tree1->SetBranchAddress(var,&variable); tree2->SetBranchAddress(var,&variable); tree3->SetBranchAddress(var,&variable); // tree4->SetBranchAddress(var,&variable); Int_t nbinsx=30; Double_t xlow,xup; if (var == "like_pi" || var == "like_k" || var == "like_proton") { xlow=0.;xup=1.; } else { xlow=2.4;xup=4.; } TH1F *h1 = new TH1F("h1", "", nbinsx, xlow, xup); TH1F *h2 = new TH1F("h2", "", nbinsx, xlow, xup); TH1F *h3 = new TH1F("h3", "", nbinsx, xlow, xup); // TH1F *h4 = new TH1F("h4", "", nbinsx, xlow, xup); for (int i = 0; i < tree1->GetEntries(); i++) { tree1->GetEntry(i); h1->Fill(variable); } for (int i = 0; i < tree2->GetEntries(); i++) { tree2->GetEntry(i); h2->Fill(variable); } for (int i = 0; i < tree3->GetEntries(); i++) { tree3->GetEntry(i); h3->Fill(variable); } /* for (int i = 0; i < tree4->GetEntries(); i++) { tree4->GetEntry(i); h4->Fill(variable); } */ normalize(h1); normalize(h2); normalize(h3); // normalize(h4); h1->SetLineColor(1); h2->SetLineColor(4); h3->SetLineColor(2); // h4->SetLineColor(4); h1->SetLineWidth(3); h2->SetLineWidth(3); h3->SetLineWidth(3); // h4->SetLineWidth(3); TCanvas* myCanvas = new TCanvas("myCanvas",var,200,10,900,800); myCanvas->cd(); // axis labels: char* xlabel = "dE/dx [MeV/cm]";; char* ylabel = "Arbitrary units"; h1->GetXaxis()->SetTitle(xlabel); h1->GetYaxis()->SetTitle(ylabel); // draw: h1->DrawCopy(); h2->DrawCopy("same"); h3->DrawCopy("same"); // h4->DrawCopy("same"); // legend: TLegend *legend = new TLegend(0.7,0.7,0.95,0.95,"","NDC"); legend->AddEntry(h1,"no smearing","l"); // legend->AddEntry(h2,"2% smearing","l"); legend->AddEntry(h2,"6.4% smearing","l"); // legend->AddEntry(h3,"1 ns async.","l"); legend->AddEntry(h3,"3 ns async.","l"); legend->Draw(); myCanvas->Update(); TString gif(var); gif+="_"; gif+=particle; gif+="_"; gif+=energymin; gif+="_"; gif+=energymax; gif+="mev"; if (mode != 0) { gif+="_mode"; gif+=mode; } gif+="_comparesmearing_pessimistic.png"; myCanvas->Print(gif); // fit (not plotted): h1->Fit("gaus"); h2->Fit("gaus"); h3->Fit("gaus"); TF1 *fit1 = h1->GetFunction("gaus"); TF1 *fit3 = h3->GetFunction("gaus"); Double_t pre=fround(fit1->GetParameter(1),4); Double_t post=fround(fit3->GetParameter(1),4); cout << "Fit: no smearing " << pre << ", async " << post << ", bias: " << 100.*(post-pre)/post << " %" << endl; delete h1; delete h2; delete h3; // delete h4; } void plotSimple(char* var="eta", int energymin=500, int energymax=5000, char* flatvar="Pt", int mode=1) { gROOT->Reset(); gStyle->SetOptStat(000000); gStyle->SetCanvasBorderMode(-1); gStyle->SetCanvasBorderSize(1); gStyle->SetCanvasColor(10); // no more grey bkg gStyle->SetTitleOffset(1.20,"Y"); // move the title // TString name1 = filename("e",energymin,energymax,flatvar,mode); TString name2 = filename("mu",energymin,energymax,flatvar,mode); TString name3 = filename("pi",energymin,energymax,flatvar,mode); TString name4 = filename("k",energymin,energymax,flatvar,mode); TString name5 = filename("proton",energymin,energymax,flatvar,mode); // TFile f1(name1); TFile f2(name2); TFile f3(name3); TFile f4(name4); TFile f5(name5); // TTree *tree1 = (TTree*)f1.Get("T"); TTree *tree2 = (TTree*)f2.Get("T"); TTree *tree3 = (TTree*)f3.Get("T"); TTree *tree4 = (TTree*)f4.Get("T"); TTree *tree5 = (TTree*)f5.Get("T"); double variable; int ivariable; if (var=="n" || var=="landaustatus") { // tree1->SetBranchAddress(var,&ivariable); tree2->SetBranchAddress(var,&ivariable); tree3->SetBranchAddress(var,&ivariable); tree4->SetBranchAddress(var,&ivariable); tree5->SetBranchAddress(var,&ivariable); } else { // tree1->SetBranchAddress(var,&variable); tree2->SetBranchAddress(var,&variable); tree3->SetBranchAddress(var,&variable); tree4->SetBranchAddress(var,&variable); tree5->SetBranchAddress(var,&variable); } Int_t nbinsx=50; if (var=="n") nbinsx=21; Double_t xlow,xup; if (var=="p") {xlow=0.; xup=10.; } if (var=="pt") {xlow=0.; xup=10.; } if (var=="eta") {xlow=-2.5; xup=2.5; } if (var=="n") {xlow=5; xup=26; } // TH1F *h1 = new TH1F("he", "", nbinsx, xlow, xup); TH1F *h2 = new TH1F("hmu", "", nbinsx, xlow, xup); TH1F *h3 = new TH1F("hpi", "", nbinsx, xlow, xup); TH1F *h4 = new TH1F("hk", "", nbinsx, xlow, xup); TH1F *h5 = new TH1F("hp", "", nbinsx, xlow, xup); // for (int i = 0; i < tree1->GetEntries(); i++) { // tree1->GetEntry(i); // if (var=="n") { h1->Fill(ivariable); } else { h1->Fill(variable); } // } for (int i = 0; i < tree2->GetEntries(); i++) { tree2->GetEntry(i); if (var=="n") { h2->Fill(ivariable); } else { h2->Fill(variable); } } for (int i = 0; i < tree3->GetEntries(); i++) { tree3->GetEntry(i); if (var=="n") { h3->Fill(ivariable); } else { h3->Fill(variable); } } for (int i = 0; i < tree4->GetEntries(); i++) { tree4->GetEntry(i); if (var=="n") { h4->Fill(ivariable); } else { h4->Fill(variable); } } for (int i = 0; i < tree5->GetEntries(); i++) { tree5->GetEntry(i); if (var=="n") { h5->Fill(ivariable); } else { h5->Fill(variable); } } // normalize(h1); normalize(h2); normalize(h3); normalize(h4); normalize(h5); // h1->SetLineColor(6); h2->SetLineColor(2); h3->SetLineColor(1); h4->SetLineColor(3); h5->SetLineColor(4); // h1->SetLineWidth(3); h2->SetLineWidth(3); h3->SetLineWidth(3); h4->SetLineWidth(3); h5->SetLineWidth(3); TCanvas* myCanvas = new TCanvas("myCanvas",var,200,10,900,800); myCanvas->cd(); // axis labels: char* xlabel; if (var=="n") xlabel = "Reconstructed hits"; char* ylabel = "Arbitrary units"; h2->GetXaxis()->SetTitle(xlabel); h2->GetYaxis()->SetTitle(ylabel); // draw: h2->DrawCopy(); h3->DrawCopy("same"); h4->DrawCopy("same"); h5->DrawCopy("same"); // h1->DrawCopy("same"); // legend: TLegend *legend = new TLegend(0.7,0.7,0.95,0.95,"particles:","NDC"); legend->AddEntry(h2,"muons","l"); legend->AddEntry(h3,"pions","l"); legend->AddEntry(h4,"kaons","l"); legend->AddEntry(h5,"protons","l"); // legend->AddEntry(h1,"electrons","l"); legend->Draw(); myCanvas->Update(); TString gif(var); gif+="_"; gif+=energymin; gif+="_"; gif+=energymax; gif+="mev"; if (flatvar=="E") gif+="_flatE"; if (mode != 0) { gif+="_mode"; gif+=mode; } gif+=".png"; myCanvas->Print(gif); // delete h1; delete h2; delete h3; delete h4; delete h5; } void profileDeltaP(char* var="dp", char* var2="pt", int energymin=500, int energymax=5000, int etamin=0, int etamax=25, int hitmin=0, char* flatvar="Pt", int mode=0) { gROOT->Reset(); gStyle->SetOptStat(0000); gStyle->SetCanvasBorderMode(-1); gStyle->SetCanvasBorderSize(1); gStyle->SetCanvasColor(10); TString name1 = filename("e",energymin,energymax,flatvar,mode); TString name2 = filename("mu",energymin,energymax,flatvar,mode); TString name3 = filename("pi",energymin,energymax,flatvar,mode); TString name4 = filename("k",energymin,energymax,flatvar,mode); TString name5 = filename("proton",energymin,energymax,flatvar,mode); TFile f1(name1); TFile f2(name2); TFile f3(name3); TFile f4(name4); TFile f5(name5); TTree *tree1 = (TTree*)f1.Get("T"); TTree *tree2 = (TTree*)f2.Get("T"); TTree *tree3 = (TTree*)f3.Get("T"); TTree *tree4 = (TTree*)f4.Get("T"); TTree *tree5 = (TTree*)f5.Get("T"); double variable; double variable2; int n; double p,eta,dr; tree1->SetBranchAddress("n",&n); tree2->SetBranchAddress("n",&n); tree3->SetBranchAddress("n",&n); tree4->SetBranchAddress("n",&n); tree5->SetBranchAddress("n",&n); tree1->SetBranchAddress("eta",&eta); tree2->SetBranchAddress("eta",&eta); tree3->SetBranchAddress("eta",&eta); tree4->SetBranchAddress("eta",&eta); tree5->SetBranchAddress("eta",&eta); tree1->SetBranchAddress("dr",&dr); tree2->SetBranchAddress("dr",&dr); tree3->SetBranchAddress("dr",&dr); tree4->SetBranchAddress("dr",&dr); tree5->SetBranchAddress("dr",&dr); tree1->SetBranchAddress(var,&variable); tree2->SetBranchAddress(var,&variable); tree3->SetBranchAddress(var,&variable); tree4->SetBranchAddress(var,&variable); tree5->SetBranchAddress(var,&variable); tree1->SetBranchAddress(var2,&variable2); tree2->SetBranchAddress(var2,&variable2); tree3->SetBranchAddress(var2,&variable2); tree4->SetBranchAddress(var2,&variable2); tree5->SetBranchAddress(var2,&variable2); Int_t nbinsx=20; Double_t xlow,xup; if (var2=="p") {xlow=0.8; xup=7.8;} if (var2=="pt") {xlow=0.8; xup=4.8;} TProfile *h1 = new TProfile("pe", var, nbinsx, xlow, xup); TProfile *h2 = new TProfile("pmu", var, nbinsx, xlow, xup); TProfile *h3 = new TProfile("ppi", var, nbinsx, xlow, xup); TProfile *h4 = new TProfile("pk", var, nbinsx, xlow, xup); TProfile *h5 = new TProfile("pp", var, nbinsx, xlow, xup); for (int i = 0; i < tree1->GetEntries(); i++) { tree1->GetEntry(i); float v=0.; bool cut=(n>hitmin && fabs(eta)>(double)etamin/10. && fabs(eta)<(double)etamax/10.); if (cut && var=="dp") v=variable/(float)n; if (cut && var=="dp2") v=variable/dr; if (cut) h1->Fill(variable2,v); } for (int i = 0; i < tree2->GetEntries(); i++) { tree2->GetEntry(i); bool cut=(n>hitmin && fabs(eta)>(double)etamin/10. && fabs(eta)<(double)etamax/10.); if (cut && var=="dp") v=variable/(float)n; if (cut && var=="dp2") v=variable/dr; if (cut) h2->Fill(variable2,v); } for (int i = 0; i < tree3->GetEntries(); i++) { tree3->GetEntry(i); bool cut=(n>hitmin && fabs(eta)>(double)etamin/10. && fabs(eta)<(double)etamax/10.); if (cut && var=="dp") v=variable/(float)n; if (cut && var=="dp2") v=variable/dr; if (cut) h3->Fill(variable2,v); } for (int i = 0; i < tree4->GetEntries(); i++) { tree4->GetEntry(i); bool cut=(n>hitmin && fabs(eta)>(double)etamin/10. && fabs(eta)<(double)etamax/10.); if (cut && var=="dp") v=variable/(float)n; if (cut && var=="dp2") v=variable/dr; if (cut) h4->Fill(variable2,v); } for (int i = 0; i < tree5->GetEntries(); i++) { tree5->GetEntry(i); bool cut=(n>hitmin && fabs(eta)>(double)etamin/10. && fabs(eta)<(double)etamax/10.); if (cut && var=="dp") v=variable/(float)n; if (cut && var=="dp2") v=variable/dr; if (cut) h5->Fill(variable2,v); } h1->SetLineColor(6); h2->SetLineColor(2); h3->SetLineColor(1); h4->SetLineColor(3); h5->SetLineColor(4); h1->SetLineWidth(3); h2->SetLineWidth(3); h3->SetLineWidth(3); h4->SetLineWidth(3); h5->SetLineWidth(3); TCanvas* myCanvas = new TCanvas("myCanvas",var,200,10,900,800); myCanvas->Divide(2,2); myCanvas->cd(1); h1->DrawCopy(); h2->DrawCopy("same"); h3->DrawCopy("same"); h4->DrawCopy("same"); h5->DrawCopy("same"); // e/pi separation: myCanvas->cd(2); TGraph* gr1 = new TGraph(); separation(gr1,h1,h3,nbinsx,xlow,xup,"electron/pion separation (RMS)",var2,21); // K/pi separation: myCanvas->cd(3); TGraph* gr2 = new TGraph(); separation(gr2,h4,h3,nbinsx,xlow,xup,"kaon/pion separation (RMS)",var2,21); // proton/pi separation: myCanvas->cd(4); TGraph* gr3 = new TGraph(); separation(gr3,h5,h3,nbinsx,xlow,xup,"proton/pion separation (RMS)",var2,21); TString gif(var); gif+="_vs_"; gif+=var2; gif+="_"; gif+=energymin; gif+="_"; gif+=energymax; gif+="mev"; if (flatvar=="E") gif+="_flatE"; if (etamin>0 || etamax<25) { gif+="_eta"; gif+=etamin; gif+="_"; gif+=etamax; } if (hitmin>5) { gif+="_hit"; gif+=hitmin; } if (mode != 0) { gif+="_mode"; gif+=mode; } gif+=".png"; myCanvas->Print(gif); delete gr1; delete gr2; delete gr3; delete h1; delete h2; delete h3; delete h4; delete h5; } TString filename(char* particle, int energymin, int energymax, char* flatvar, int mode, double smear=0) { TString pippo("trackhisto_"); pippo+=particle; pippo+=energymin; pippo+="_"; pippo+=energymax; pippo+="mev"; if (flatvar=="E") pippo+="_flatE"; pippo+="_mode"; pippo+=mode; if (smear != 0) { pippo+="_smear"; if (smear==-0.4) { pippo+="-0.4"; } else if (smear==6.4) { pippo+="6.4"; } else { pippo+=(int)smear; } } pippo+=".root"; // pippo+="_peak2.root"; // temp return pippo; } void fillProfile(TProfile* prof, char* var="mean", char* var2="p", TTree *tree, int etamin=0, int etamax=25, int hitmin=0, char* label="muons", int color=1, char* superimpose="") { gStyle->SetOptStat(0000); gStyle->SetCanvasBorderMode(-1); gStyle->SetCanvasBorderSize(1); gStyle->SetCanvasColor(10); if (var!="nn1" && var!="nn2") { // normal (one of the dE/dx estimators) double variable; double variable2; double p,pt,eta; int n; tree->SetBranchAddress(var,&variable); if (var2!="n") tree->SetBranchAddress(var2,&variable2); tree->SetBranchAddress("n",&n); if (var2!="eta") tree->SetBranchAddress("eta",&eta); if (var2!="p") tree->SetBranchAddress("p",&p); // int landaustatus=0; // if (var=="landaufit") tree->SetBranchAddress("landaustatus",&landaustatus); int outliers=0; for (int i = 0; i < tree->GetEntries(); i++) { tree->GetEntry(i); if (label=="electrons" && variable>50.) cout << "attenzione! " << variable << endl; if (n>2 && (variable<0 || variable > 50)) { outliers++; cout << "outlier; n = " << n << ", p = " << variable2 << ", eta = " << eta << ", dE/dx = " << variable << endl; } if (var2=="eta") { if (n>hitmin && p<1.5) prof->Fill(fabs(variable2),variable); } else if (var2=="n") { if (n>hitmin && p<1.5) prof->Fill(n,variable); } else { // cout << "landaustatus = " << landaustatus << endl; if (variable<50 && variable>0 && n>hitmin && fabs(eta)>(double)etamin/10. && fabs(eta)<(double)etamax/10.) prof->Fill(variable2,variable); } } cout << " # of outliers: " << outliers << " out of " << tree->GetEntries() << " tracks" << endl; // axis labels: char* xlabel="p [GeV/c]"; // temporary! // if (abscissa=="p") {xlabel="p [GeV/c]";} else if (abscissa=="pt") {xlabel="p_T [GeV/c]";} else if (abscissa=="\eta") {xlabel="\eta";} else if (abscissa=="n") {xlabel="n.hits";} else {continue;} char* ylabel = "dE/dx [MeV/cm]"; prof->GetXaxis()->SetTitle(xlabel); prof->GetYaxis()->SetTitle(ylabel); } else { // NN1 (proton vs pi) or NN2 (e vs pi) if (!gROOT->GetClass("TMultiLayerPerceptron")) { gSystem->Load("libMLP"); cout << endl << "ROOT neural nets library has just been loaded." << endl << endl; } int neuro1,neuro2,ntrain,nAlgo; if (var=="nn1") {// proton vs pi nAlgo=3; neuro1=4; neuro2=0; ntrain=2800; } if (var=="nn2") {// e vs pi nAlgo=3; neuro1=4; neuro2=0; ntrain=3000; } TString structure("n,p,eta,dpdr,dp2,dedx"); // inputs structure+=":"; structure+=neuro1; if (neuro2>0) { // two layers structure+=":"; structure+=neuro2; } structure+=":"; structure+="truth"; //target output; @ to normalize //create new tree with the right format TString outfile("temp_"); outfile+=var; outfile+=".root"; convertTree4NN(tree,outfile); //open the new tree double n,p,eta,dp2,dpdr,dedx,truth; TFile f(outfile); TTree *newtree = (TTree*)gDirectory->Get("T"); newtree->SetBranchAddress("n",&n); newtree->SetBranchAddress("p",&p); newtree->SetBranchAddress("eta",&eta); newtree->SetBranchAddress("dpdr",&dpdr); newtree->SetBranchAddress("dp2",&dp2); newtree->SetBranchAddress("dedx",&dedx); newtree->SetBranchAddress("truth",&truth); TMultiLayerPerceptron *mlp = new TMultiLayerPerceptron(structure,newtree); TString algo; if (nAlgo == 1) algo="Stochastic"; if (nAlgo == 2) algo="Batch"; if (nAlgo == 3) algo="SteepestDescent"; if (nAlgo == 4) algo="RibierePolak"; if (nAlgo == 5) algo="FletcherReeves"; if (nAlgo == 6) algo="BFGS"; cout << "Algorithm: " << algo << endl; TString nomepesi("txt/"); if (var=="nn1") nomepesi+="proton_vs_pi"; if (var=="nn2") nomepesi+="e_vs_pi"; nomepesi+="_"; nomepesi+=algo; nomepesi+="_"; nomepesi+=neuro1; nomepesi+="_"; nomepesi+=neuro2; nomepesi+="_"; nomepesi+=ntrain; nomepesi+=".txt"; mlp->LoadWeights(nomepesi); cout << "File with weights: " << nomepesi << endl; Double_t nn; Double_t inputs[8]; for (int i = 0; i < newtree->GetEntries(); i++) { newtree->GetEntry(i); // evaluate NN output: inputs[0]=n; inputs[1]=p; inputs[2]=eta; inputs[3]=dpdr; inputs[4]=dp2; inputs[5]=dedx; nn = mlp->Evaluate(0,inputs); // apply cuts: vector w; w.clear(); w.push_back(n); w.push_back(eta); w.push_back(p); w.push_back(dpdr); w.push_back(dp2); w.push_back(dedx); vector v=invert(w); if (var2=="eta") { if (v.at(0)>hitmin && v.at(2)<1.5) prof->Fill(fabs(v.at(1)),nn); } else if (var2=="n") { if (v.at(2)<1.5) prof->Fill(v.at(0),nn); } else { Double_t variable2; if (var2=="p") variable2=v.at(2); if (v.at(0)>hitmin && fabs(v.at(1))>(double)etamin/10. && fabs(v.at(1))<(double)etamax/10.) prof->Fill(variable2,nn); } } } prof->SetErrorOption("s"); // to have spread, and not spread/sqrt(N), as error. // prof->SetMarkerStyle(21); prof->SetMarkerColor(color); prof->SetLineColor(color); prof->SetLineWidth(3); prof->DrawCopy(superimpose); } Double_t findmax(TProfile* prof, Int_t nbins) { Double_t max=0; for (unsigned int i=0; iGetBinContent(i); Double_t temp2=prof->GetBinError(i); if (temp1+temp2 > max) max=temp1+temp2; } return max; } Double_t findmin(TProfile* prof, Int_t nbins) { Double_t min=99999.; for (unsigned int i=0; iGetBinContent(i); Double_t temp2=prof->GetBinError(i); if (temp1-temp2 < min) min=temp1-temp2; } return min; } void separation(TGraph* gr, TProfile* prof1, TProfile* prof2, Int_t nbins, Double_t xlow, Double_t xup, char* label="proton/pion separation", char* abscissa="p",int marker=21) { Double_t binwidth=(xup-xlow)/nbins; Double_t* x = NULL; Double_t* y = NULL; x = new Double_t[nbins]; y = new Double_t[nbins]; for (unsigned int i=0; iGetBinContent(i+1); Double_t y2=prof2->GetBinContent(i+1); Double_t e1=prof1->GetBinError(i+1); Double_t e2=prof2->GetBinError(i+1); if (e1>0. && e2>0.) { y[i]= (y1-y2)/sqrt(e1**2+e2**2); // cout << "i = "<< i << ", x = " << x[i] << ", y = " << y[i]<< endl; gr->SetPoint(i,x[i],y[i]); } } gr->SetMarkerStyle(marker); gr->SetTitle(label); if (abscissa=="p") {gr->GetXaxis()->SetTitle("p [GeV/c]");} else if (abscissa=="pt") {gr->GetXaxis()->SetTitle("p_T [GeV/c]");} else if (abscissa=="\eta") {gr->GetXaxis()->SetTitle("\eta");} else if (abscissa=="n") {gr->GetXaxis()->SetTitle("n.hits");} else {continue;} // gr->GetXaxis()->SetTitle("P [GeV]"); gr->GetYaxis()->SetTitle("proton-pion separation (RMS)"); gr->Draw("AP"); delete [] x; delete [] y; } void fillHisto(TH1F* h, char* var="mean", TTree *tree, int pmin, int pmax, int etamin, int etamax, int hitmin, bool norma, bool log, char* label="muons", int mode=1, double smear=0, int color=1, char* superimpose="") { gStyle->SetOptStat(0000); gStyle->SetCanvasBorderMode(-1); gStyle->SetCanvasBorderSize(1); gStyle->SetCanvasColor(10); if (var!="nn1" && var!="nn2") { // normal (one of the dE/dx estimators) double variable; double p,eta; int n; tree->SetBranchAddress("n",&n); tree->SetBranchAddress("p",&p); tree->SetBranchAddress("eta",&eta); tree->SetBranchAddress(var,&variable); for (int i = 0; i < tree->GetEntries(); i++) { tree->GetEntry(i); if (n>hitmin && p > (double)pmin/1000. && p < (double)pmax/1000. && fabs(eta)>(double)etamin/10. && fabs(eta)<(double)etamax/10.) { if (log) { if (variable>0) h->Fill(log(variable)); } else { h->Fill(variable); } } } // axis labels: char* xlabel; if (log) { xlabel = "log(dE/dx [MeV/cm])"; } else { xlabel = "dE/dx [MeV/cm]"; } char* ylabel = "Arbitrary units"; h->GetXaxis()->SetTitle(xlabel); h->GetYaxis()->SetTitle(ylabel); } else { // NN1 (proton vs pi) or NN2 (e vs pi) if (!gROOT->GetClass("TMultiLayerPerceptron")) { gSystem->Load("libMLP"); cout << endl << "ROOT neural nets library has just been loaded." << endl << endl; } int neuro1,neuro2,ntrain,nAlgo; if (var=="nn1") {// proton vs pi nAlgo=3; neuro1=4; neuro2=0; ntrain=2800; } if (var=="nn2") {// e vs pi nAlgo=3; neuro1=4; neuro2=0; ntrain=3000; } TString structure("n,p,eta,dpdr,dp2,dedx"); // inputs structure+=":"; structure+=neuro1; if (neuro2>0) { // two layers structure+=":"; structure+=neuro2; } structure+=":"; structure+="truth"; //target output; @ to normalize //create new tree with the right format TString outfile("temp_"); outfile+=var; outfile+=".root"; convertTree4NN(tree,outfile); //open the new tree double n,p,eta,dp2,dpdr,dedx,truth; TFile f(outfile); TTree *newtree = (TTree*)gDirectory->Get("T"); newtree->SetBranchAddress("n",&n); newtree->SetBranchAddress("p",&p); newtree->SetBranchAddress("eta",&eta); newtree->SetBranchAddress("dpdr",&dpdr); newtree->SetBranchAddress("dp2",&dp2); newtree->SetBranchAddress("dedx",&dedx); newtree->SetBranchAddress("truth",&truth); TMultiLayerPerceptron *mlp = new TMultiLayerPerceptron(structure,newtree); TString algo; if (nAlgo == 1) algo="Stochastic"; if (nAlgo == 2) algo="Batch"; if (nAlgo == 3) algo="SteepestDescent"; if (nAlgo == 4) algo="RibierePolak"; if (nAlgo == 5) algo="FletcherReeves"; if (nAlgo == 6) algo="BFGS"; cout << "Algorithm: " << algo << endl; TString nomepesi("txt/"); if (var=="nn1") nomepesi+="proton_vs_pi"; if (var=="nn2") nomepesi+="e_vs_pi"; nomepesi+="_"; nomepesi+=algo; nomepesi+="_"; nomepesi+=neuro1; nomepesi+="_"; nomepesi+=neuro2; nomepesi+="_"; nomepesi+=ntrain; nomepesi+=".txt"; mlp->LoadWeights(nomepesi); cout << "File with weights: " << nomepesi << endl; Double_t inputs[8]; bool cut=true; for (int i = 0; i < newtree->GetEntries(); i++) { newtree->GetEntry(i); // apply cuts: vector w; w.clear(); w.push_back(n); w.push_back(eta); w.push_back(p); w.push_back(dpdr); w.push_back(dp2); w.push_back(dedx); vector v=invert(w); // cut = (v.at(2)>pmin && v.at(2)etamin && fabs(v.at(1))hitmin && v.at(2) > (double)pmin/1000. && v.at(2) < (double)pmax/1000. && fabs(v.at(1))>(double)etamin/10. && fabs(v.at(1))<(double)etamax/10.); // evaluate NN output: if (cut) { inputs[0]=n; inputs[1]=p; inputs[2]=eta; inputs[3]=dpdr; inputs[4]=dp2; inputs[5]=dedx; h->Fill(mlp->Evaluate(0,inputs)); } } delete newtree; cout << "...deleting file " << outfile << "..." << endl; TString del("rm -f "); del+=outfile; gSystem->Exec(del); } if (norma) normalize(h); h->SetLineColor(color); h->SetLineWidth(3); h->DrawCopy(superimpose); if (log) { // perform a gaussian fit: h->Fit("gaus"); TF1 *fit = h->GetFunction("gaus"); Double_t chi2 = fit->GetChisquare()/(float)fit->GetNDF(); unsigned ndecimals=2; Double_t p1 = fit->GetParameter(1); Double_t e1 = fit->GetParError(1); Double_t p2 = fit->GetParameter(2); Double_t e2 = fit->GetParError(2); cout << "--------------------" << endl; cout << "Result of the fit:" << endl; cout << "--------------------" << endl; cout << "chi2 = " << fround(chi2,1) << endl; cout << "mean = " << fround(p1,ndecimals) << "+-" << fround(e1,ndecimals) << endl; cout << "sigma = " << fround(p2,ndecimals) << "+-" << fround(e2,ndecimals) << endl; cout << "--------------------" << endl; // draw a legend with the fit results: Double_t xmin=1.8; Double_t xmax=3.3; Double_t ymax=h->GetMaximum(); TPaveText *text = new TPaveText(0.65,0.65,.95,0.95,"NDC"); TString t0("chi2 = "); t0+=fround(chi2,1); TString t1("mean = "); t1+=fround(p1,ndecimals); // t1+="+-"; // t1+=e1; TString t2("sigma = "); t2+=fround(p2,ndecimals); // t2+="+-"; // t2+=e2; text->AddText(xmin,0.65*ymax,label); text->AddText(xmin,0.65*ymax,t0); text->AddText(xmin,0.65*ymax,t1); text->AddText(xmin,0.65*ymax,t2); text->SetTextSize(0.05); text->Draw(); // write the fit results in an ascii file TString ascii("txt/"); ascii+=var; ascii+="_"; ascii+=label; ascii+="_P"; ascii+=pmin; ascii+="_"; ascii+=pmax; if (etamin>0 || etamax<25) { ascii+="_eta"; ascii+=etamin; ascii+="_"; ascii+=etamax; } if (hitmin>5) { ascii+="_hit"; ascii+=hitmin; } if (mode != 0) { ascii+="_mode"; ascii+=mode; } if (smear != 0) { ascii+="_smear"; ascii+=(int)smear; } ascii+=".txt"; fstream f1; f1.open(ascii,ios::out); f1 << chi2 << endl; f1 << p1 << endl; f1 << p2 << endl; f1.close(); cout << "Fit results written into: " << ascii << endl; } // delete h; } void final(char* var="mean", int pmin=900, int pmax=1100, int etamin=0, int etamax=25, int hitmin=3, int mode=1, double smear=0) { const int numfiles = 5; Double_t chi2[numfiles], p1[numfiles], p2[numfiles]; TString label[numfiles]={"muons","electrons","pions","kaons","protons"}; Double_t separation[numfiles]; cout << "***************************" << endl; cout << "Estimator: " << var << endl; cout << "***************************" << endl; for (unsigned int i=0; i0 || etamax<25) { ascii+="_eta"; ascii+=etamin; ascii+="_"; ascii+=etamax; } if (hitmin>5) { ascii+="_hit"; ascii+=hitmin; } if (mode != 0) { ascii+="_mode"; ascii+=mode; } if (smear != 0) { ascii+="_smear"; ascii+=(int)smear; } ascii+=".txt"; cout << "Reading file: " << ascii << endl; ifstream f1; f1.open(ascii); if (f1.good()) { f1 >> chi2[i]; f1 >> p1[i]; f1 >> p2[i]; f1.close(); if (i>0) separation[i]=(p1[i]-p1[0])/sqrt(p2[i]**2 + p2[0]**2); } else { cout << ascii << " doesn't exist!" << endl; } } cout << "muon chi2: " << chi2[0] << endl; cout << "muon sigma: " << 100.*p2[0]/p1[0] << " %" << endl; cout << "pion/muon separation: " << separation[2] << " sigmas" << endl; cout << "kaon/muon separation: " << separation[3] << " sigmas" << endl; cout << "proton/muon separation: " << separation[4] << " sigmas" << endl; cout << "kaon/pion separation:" << (p1[3]-p1[2])/sqrt(p2[3]**2 + p2[2]**2) << " sigmas" << endl; cout << "proton/pion separation:" << (p1[4]-p1[2])/sqrt(p2[4]**2 + p2[2]**2) << " sigmas" << endl; cout << "proton/kaon separation:" << (p1[4]-p1[3])/sqrt(p2[4]**2 + p2[3]**2) << " sigmas" << endl; cout << "***************************" << endl; } void nnScatter(TTree *tree, int pmin=900, int pmax=1100, int etamin=0, int etamax=25, int hitmin=0, int color=1, char* superimpose="") { //// WORK IN PROGRESS if (!gROOT->GetClass("TMultiLayerPerceptron")) { gSystem->Load("libMLP"); cout << endl << "ROOT neural nets library has just been loaded." << endl << endl; } // proton vs pi int nAlgo_nn1=3; int neuro1_nn1=4; int neuro2_nn1=0; int ntrain_nn1=2800; TString structure_nn1("n,p,eta,dpdr,dp2,dedx"); // inputs structure_nn1+=":"; structure_nn1+=neuro1_nn1; if (neuro2_nn1>0) { // two layers structure_nn1+=":"; structure_nn1+=neuro2_nn1; } structure_nn1+=":"; structure_nn1+="truth"; //target output; @ to normalize // e vs pi int nAlgo_nn2=3; int neuro1_nn2=4; int neuro2_nn2=0; int ntrain_nn2=3000; TString structure_nn2("n,p,eta,dpdr,dp2,dedx"); // inputs structure_nn2+=":"; structure_nn2+=neuro1_nn2; if (neuro2_nn2>0) { // two layers structure_nn2+=":"; structure_nn2+=neuro2_nn2; } structure_nn2+=":"; structure_nn2+="truth"; //target output; @ to normalize //create new tree with the right format TString outfile("temp_"); outfile+="nn1_nn2"; outfile+=".root"; convertTree4NN(tree,outfile); //open the new tree double n,p,eta,dp2,dpdr,dedx,truth; TFile f(outfile); TTree *newtree = (TTree*)gDirectory->Get("T"); newtree->SetBranchAddress("n",&n); newtree->SetBranchAddress("p",&p); newtree->SetBranchAddress("eta",&eta); newtree->SetBranchAddress("dpdr",&dpdr); newtree->SetBranchAddress("dp2",&dp2); newtree->SetBranchAddress("dedx",&dedx); newtree->SetBranchAddress("truth",&truth); TMultiLayerPerceptron *mlp_nn1 = new TMultiLayerPerceptron(structure_nn1,newtree); TMultiLayerPerceptron *mlp_nn2 = new TMultiLayerPerceptron(structure_nn2,newtree); TString algo_nn1; if (nAlgo_nn1 == 1) algo_nn1="Stochastic"; if (nAlgo_nn1 == 2) algo_nn1="Batch"; if (nAlgo_nn1 == 3) algo_nn1="SteepestDescent"; if (nAlgo_nn1 == 4) algo_nn1="RibierePolak"; if (nAlgo_nn1 == 5) algo_nn1="FletcherReeves"; if (nAlgo_nn1 == 6) algo_nn1="BFGS"; cout << "Algorithm NN1: " << algo_nn1 << endl; TString algo_nn2; if (nAlgo_nn2 == 1) algo_nn2="Stochastic"; if (nAlgo_nn2 == 2) algo_nn2="Batch"; if (nAlgo_nn2 == 3) algo_nn2="SteepestDescent"; if (nAlgo_nn2 == 4) algo_nn2="RibierePolak"; if (nAlgo_nn2 == 5) algo_nn2="FletcherReeves"; if (nAlgo_nn2 == 6) algo_nn2="BFGS"; cout << "Algorithm NN2: " << algo_nn2 << endl; TString nomepesi_nn1("txt/proton_vs_pi"); nomepesi_nn1+="_"; nomepesi_nn1+=algo_nn1; nomepesi_nn1+="_"; nomepesi_nn1+=neuro1_nn1; nomepesi_nn1+="_"; nomepesi_nn1+=neuro2_nn1; nomepesi_nn1+="_"; nomepesi_nn1+=ntrain_nn1; nomepesi_nn1+=".txt"; mlp_nn1->LoadWeights(nomepesi_nn1); cout << "File with weights, NN1: " << nomepesi_nn1 << endl; TString nomepesi_nn2("txt/e_vs_pi"); nomepesi_nn2+="_"; nomepesi_nn2+=algo_nn2; nomepesi_nn2+="_"; nomepesi_nn2+=neuro1_nn2; nomepesi_nn2+="_"; nomepesi_nn2+=neuro2_nn2; nomepesi_nn2+="_"; nomepesi_nn2+=ntrain_nn2; nomepesi_nn2+=".txt"; mlp_nn2->LoadWeights(nomepesi_nn2); cout << "File with weights, NN2: " << nomepesi_nn2 << endl; TH2F *h = new TH2F("h", "NN1 vs NN2",100,-0.5,1.5,100,-0.5,1.5); Double_t nn1; Double_t nn2; Double_t inputs[8]; bool cut=true; for (int i = 0; i < newtree->GetEntries(); i++) { newtree->GetEntry(i); // apply cuts: vector w; w.clear(); w.push_back(n); w.push_back(eta); w.push_back(p); w.push_back(dpdr); w.push_back(dp2); w.push_back(dedx); vector v=invert(w); cut = (v.at(0)>hitmin && v.at(2) > (double)pmin/1000. && v.at(2) < (double)pmax/1000. && fabs(v.at(1))>(double)etamin/10. && fabs(v.at(1))<(double)etamax/10.); // evaluate NN output: if (cut) { inputs[0]=n; inputs[1]=p; inputs[2]=eta; inputs[3]=dpdr; inputs[4]=dp2; inputs[5]=dedx; nn1 = mlp_nn1->Evaluate(0,inputs); nn2 = mlp_nn2->Evaluate(0,inputs); h->Fill(nn1,nn2); } } h->SetMarkerStyle(23); // 21: squares, 20, 22: asymmetric stars, 23: triangles (pointing downwards) h->SetMarkerSize(0.5); h->SetMarkerColor(color); h->DrawCopy(superimpose); delete h; delete mlp_nn1; delete mlp_nn2; delete newtree; cout << "...deleting file " << outfile << "..." << endl; TString del("rm -f "); del+=outfile; gSystem->Exec(del); } void plotScatter(char* var="mean", TTree *tree, int hitmin=3, int color=1, char* superimpose="") { gStyle->SetOptStat(0000); gStyle->SetCanvasBorderMode(-1); gStyle->SetCanvasBorderSize(1); gStyle->SetCanvasColor(10); double variable; double p,eta; int n; tree->SetBranchAddress(var,&variable); tree->SetBranchAddress("p",&p); tree->SetBranchAddress("eta",&eta); tree->SetBranchAddress("n",&n); TH2F *h = new TH2F("h", "dE/dx [MeV/cm] vs P (GeV)",100,0.4,7.0,100,1.5,6.0); for (int i = 0; i < tree->GetEntries(); i++) { tree->GetEntry(i); if (n>hitmin) h->Fill(p,variable); } // axis labels: char* xlabel = "p [GeV/c]"; char* ylabel = "dE/dx [MeV/cm]";; h->GetXaxis()->SetTitle(xlabel); h->GetYaxis()->SetTitle(ylabel); h->SetMarkerStyle(23); // 21: squares, 20, 22: asymmetric stars, 23: triangles (pointing downwards) h->SetMarkerSize(0.5); h->SetMarkerColor(color); h->DrawCopy(superimpose); delete h; } /* round number n to d decimal points */ double fround(double n, unsigned d) { return floor(n * pow(10., d) + .5) / pow(10., d); } void normalize(TH1F* h) { Double_t scale; if (h->Integral()!=0) scale = 1/h->Integral(); h->Scale(scale); } void plotHistoHits(char* var="dedx", char* type="pxb", char* particle="mu", int energymin=1000, int energymax=5000, char* flatvar="Pt", int mode=1, int color=1, char* superimpose="") { gStyle->SetTitleH(0.8); TString name1 = filename(particle,energymin,energymax,flatvar,mode); TFile* f1 = new TFile(name1); TString histo("h"); histo+=var; histo+="_"; histo+=type; TH1F *h = (TH1F*)f1->Get(histo); h->SetLineColor(color); h->SetLineWidth(3); normalize(h); // axis labels: char* xlabel; if (var=="dedx") xlabel = "dE/dx [MeV/cm]"; if (var=="N") xlabel = "n. of hits"; char* ylabel = "Arbitrary units"; h->GetXaxis()->SetLabelSize(0.08); h->GetYaxis()->SetLabelSize(0.08); h->GetXaxis()->SetTitleSize(0.08); h->GetYaxis()->SetTitleSize(0.08); // h->GetXaxis()->SetTitle(xlabel); // h->GetYaxis()->SetTitle(ylabel); h->DrawCopy(superimpose); delete h; } vector regularize(vector v) { vector w; w.clear(); double n=v.at(0); double eta=v.at(1); double p=v.at(2); double dpdr=v.at(3); double dp2=v.at(4); double dedx=v.at(5); w.push_back((1.*n-15.)/8.); w.push_back(eta/2.4); w.push_back((atan(p)-1.2)/.4); w.push_back((dpdr-0.007)/0.005); w.push_back(dp2/0.1); w.push_back(log(dedx)-1.); return w; } vector invert(vector v) { vector w; w.clear(); double n=v.at(0); double eta=v.at(1); double p=v.at(2); double dpdr=v.at(3); double dp2=v.at(4); double dedx=v.at(5); w.push_back(8.*n+15.); w.push_back(eta*2.4); w.push_back(tan(0.4*p+1.2)); w.push_back(dpdr*0.005+0.007); w.push_back(dp2*0.1); w.push_back(exp(dedx+1.)); return w; } void convertTree4NN(TTree* tree, Char_t* outputfile) { cout << "...creating file " << outputfile << "..." << endl; // input tree: Int_t n; Double_t eta; Double_t p; Double_t dp; Double_t dr; Double_t dp2; Double_t harm2; tree->SetBranchAddress("n",&n); tree->SetBranchAddress("eta",&eta); tree->SetBranchAddress("p",&p); tree->SetBranchAddress("dp",&dp); tree->SetBranchAddress("dr",&dr); tree->SetBranchAddress("dp2",&dp2); tree->SetBranchAddress("harm2",&harm2); // new tree: Double_t n_copy;//note int->double conversion! Double_t eta_copy; Double_t p_copy; Double_t dpdr; Double_t dp2_copy; Double_t dedx; Double_t truth; TFile *f = new TFile(outputfile,"RECREATE"); TTree* nntree = new TTree("T","nn inputs"); nntree->Branch("n",&n_copy,"n_copy/D"); nntree->Branch("eta",&eta_copy,"eta_copy/D"); nntree->Branch("p",&p_copy,"p_copy/D"); nntree->Branch("dpdr",&dpdr,"dpdr/D"); nntree->Branch("dp2",&dp2_copy,"dp2_copy/D"); nntree->Branch("dedx",&dedx,"dedx/D"); nntree->Branch("truth",&truth,"truth/D"); // loop on input tree: truth=-1.; // neither signal nor background for (int i = 0; i < tree->GetEntries(); i++) { tree->GetEntry(i); vector v; v.clear(); v.push_back(1.*n); v.push_back(eta); v.push_back(p); v.push_back(dp/dr); v.push_back(dp2); v.push_back(harm2); vector w=regularize(v); n_copy=w.at(0); eta_copy=w.at(1); p_copy=w.at(2); dpdr=w.at(3); dp2_copy=w.at(4); dedx=w.at(5); bool ok=true; for (int k=0;k10.) ok=false;} // exclude outliers! if (ok) nntree->Fill(); } nntree->Write(); f->Write(); f->Close(); } void define_style(){ /* if (getenv("CMSSW_RELEASE_BASE")!='\0'){ printf("CMSSW FWLite loaded\n\n"); gSystem->Load("libFWCoreFWLite"); AutoLibraryLoader::enable(); } cout << "\nFor approved plots use: gROOT->SetStyle(\"BABAR\");" << endl << endl; */ //..BABAR style from RooLogon.C in workdir TStyle *babarStyle= new TStyle("BABAR","BaBar approved plots style"); // use plain black on white colors babarStyle->SetFrameBorderMode(0); babarStyle->SetCanvasBorderMode(0); babarStyle->SetPadBorderMode(0); babarStyle->SetPadColor(0); babarStyle->SetCanvasColor(0); babarStyle->SetStatColor(0); babarStyle->SetFillColor(0); // set the paper & margin sizes babarStyle->SetPaperSize(20,26); babarStyle->SetPadTopMargin(0.05); babarStyle->SetPadRightMargin(0.05); babarStyle->SetPadBottomMargin(0.16); babarStyle->SetPadLeftMargin(0.12); // use large Times-Roman fonts babarStyle->SetTextFont(132); babarStyle->SetTextSize(0.08); babarStyle->SetLabelFont(132,"x"); babarStyle->SetLabelFont(132,"y"); babarStyle->SetLabelFont(132,"z"); babarStyle->SetLabelSize(0.05,"x"); babarStyle->SetTitleSize(0.06,"x"); babarStyle->SetLabelSize(0.05,"y"); babarStyle->SetTitleSize(0.06,"y"); babarStyle->SetLabelSize(0.05,"z"); babarStyle->SetTitleSize(0.06,"z"); // use bold lines and markers babarStyle->SetMarkerStyle(20); babarStyle->SetHistLineWidth(1.85); babarStyle->SetLineStyleString(2,"[12 12]"); // postscript dashes // get rid of X error bars and y error bar caps babarStyle->SetErrorX(0.001); // do not display any of the standard histogram decorations babarStyle->SetOptTitle(0); babarStyle->SetOptStat(0); babarStyle->SetOptFit(0); // put tick marks on top and RHS of plots babarStyle->SetPadTickX(1); babarStyle->SetPadTickY(1); gROOT->SetStyle("Plain"); gStyle->SetOptStat(1111111); gStyle->SetPadTickX(1); gStyle->SetPadTickY(1); }