Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:11:15 

0001 #include "sPhenixStyle.C"
0002 #include "sPhenixStyle.h"
0003 
0004 TPad* getPad(int j, int k, int l, float cw, float ch, const int  nEBins, const int nCutBins);
0005 
0006 //void drawCanvas_invMass(TCanvas *c, TH1F *hInvMass, int pad_x, int pad_y, TPad *pad, int isEta);
0007 void drawCanvas_invMass(TCanvas *c, TH1F *hCorr, int pad_x, int pad_y, TPad *pad, int isEta, float peakPos, float peakW);
0008 
0009 void SetHistoStyle(TH1F *histo, int cutBin, int eBin, float low, float hi, int isEta);
0010 
0011 void DrawCutLabel(int x, int y, int isEta);
0012 
0013 void GetSubtractedDist(TH1F *histOrig, TH1F *histSub, TF1 *invMassFit, TF1 *invMassBG);
0014 
0015 //float eBins[] = {1,3,5,7,9,11,13,15,17,20};
0016 float eBins[] = {1,2,3,4,5,6,7,8,9,10,11,12,13};
0017 float eBinsEta[] = {12,14,16,18,20};
0018 float eCuts[] = {0.5,0.6,0.7,0.8,0.9,1.,1.1/*,1.2,1.3,1.4,1.5*/};
0019 float eCutsEta[] = {3,3.5,4,4.5,5,5.5,6};
0020 const int nEtaBins = 5;
0021 void draw_hists_3D(const char *input)
0022 {
0023 
0024   SetsPhenixStyle();
0025   //TH1::SetDefaultSumw2();
0026   //TH2::SetDefaultSumw2();
0027   //TH3::SetDefaultSumw2();
0028   
0029   TFile *fin = new TFile(input);
0030  
0031   string th1List[] = {"photonE_Reco","isoPhotonE_Reco_Tru","truthPi0E","truthDphoE","pi0ETruthReco"};
0032 
0033   string th1xAxis[] = {"E_{#gamma}","E_{#gamma}^{Iso}","E_{#pi^{0}}^{Truth}","E_{#gamma_{Dir}}^{Truth}","E_{#pi^{0}}^{Reco}/E_{#pi^{0}}^{Truth}"};
0034   
0035   string th2List[] = {"clusterChi2","clusterProbPhoton","isoPhotonChi2","isoPhotonProb2","tspE","tspEiso","asymEtruthpi0","deltaREtruthpi0"};
0036   
0037   string th2xAxis[] = {"Cluster #xhi^{2}","Cluster Photon Prob","Iso Cluster #chi^{2}","Iso Cluster Prob","TSP","TSP","#frac{|E_{1} - E_{2}|}{E_{#pi^0}}","#Delta R [rad]"};
0038   
0039   string th2yAxis[] = {"E_{Cluster}","E_{Cluster}","E_{Cluster}","E_{Cluster}","E_{Cluster}","E_{Cluster}","E_{#pi^{0}}","E_{#pi^{0}}"};
0040   
0041   string th3List[] = {"deltaREinvMass","asymEinvMass","invMassEtaE","dphoChi2","dphoProb","pi0Chi2","pi0Prob","etaChi2","etaProb","eChi2","eProb","hChi2","hProb","ePi0InvMassEcut0"};
0042 
0043   string th3xAxis[] = {"#DeltaR","#frac{|E_{1} - E_{2}|}{E_{#pi^0}}","M_{#gamma#gamma}","#chi^{2}","Cluster Photon Prob","#chi^{2}","Cluster Photon Prob","#chi^{2}","Cluster Photon Prob","#chi^{2}","Cluster Photon Prob","#chi^{2}","Cluster Photon Prob","E_{#pi^{0}}^{Reco}"};
0044 
0045   string th3yAxis[] = {"E_{#pi^{0}}^{Reco}","E_{#pi^{0}}^{Reco}","#eta"," E_{Cluster}^{#gamma} [GeV]","E_{Cluster}^{#gamma}[GeV]","E_{Cluster}^{#pi^{0}}[GeV]","E_{Cluster}^{#pi^{0}}[GeV]","E_{Cluster}^{#eta}[GeV]","E_{Cluster}^{#eta}[GeV]","E_{Cluster}^{e^{#pm}}[GeV]","E_{Cluster}^{e^{#pm}}[GeV]","E_{Cluster}^{h^{#pm}}[GeV]","E_{Cluster}^{h^{#pm}}[GeV]","M_{#gamma#gamma}"};
0046 
0047   string th3zAxis[] = {"M_{#gamma#gamma}[GeV/c^{2}]","M_{#gamma#gamma}[GeV/c^{2}]","E^{#gamma}_{Reco}[GeV]","E^{#gamma}_{True}[GeV]","E^{#gamma}_{True}[GeV]","E^{#pi^{0}}_{True}[GeV]","E^{#pi^{0}}_{True}[GeV]","E^{#eta}_{True}[GeV]","E^{#eta}_{True}[GeV]","E^{e^{#pm}}_{True}[GeV]","E_{e^{#pm}}_{True}[GeV]","E_{#pi^{#pm}}_{True}[GeV]","E^{#pi^{#pm}}_{True}[GeV]","E_{Cut}[GeV]"};
0048  
0049   
0050   //loop over and draw the TH1F's
0051   /*TCanvas *c1 = new TCanvas();
0052   const int nTh1 = (int)sizeof(th1List)/sizeof(th1List[0]);
0053   for(int i = 0; i < nTh1; i++)
0054     {
0055       TH1F *hist = (TH1F*)fin -> Get(th1List[i].c_str());
0056       hist -> GetXaxis() -> SetTitle(th1xAxis[i].c_str());
0057       gPad -> SetLogy();
0058       hist -> Draw();
0059       c1 -> SaveAs(Form("plots/%s.pdf",th1List[i].c_str()));
0060     }
0061   gPad -> SetLogy(0);
0062   const int nTh2 = (int)sizeof(th2List)/sizeof(th2List[0]);
0063   for(int i = 0; i < nTh2; i++)
0064     {
0065       gPad -> SetLogz();
0066       TH2F *hist = (TH2F*)fin -> Get(th2List[i].c_str());
0067       hist -> GetYaxis() -> SetTitle(th2yAxis[i].c_str());
0068       hist -> GetXaxis() -> SetTitle(th2xAxis[i].c_str());
0069       hist -> Draw("colz");
0070       c1 -> SaveAs(Form("plots/%s.pdf",th2List[i].c_str()));
0071       }*/
0072   
0073   
0074   const int nTh3 = (int)sizeof(th3List)/sizeof(th3List[0]);
0075   const int nEBins = (int)sizeof(eBins)/sizeof(eBins[0]);
0076   const int nCutBins = (int)sizeof(eCuts)/sizeof(eCuts[0]);
0077   const int nCutBinsEta = (int)sizeof(eCutsEta)/sizeof(eCutsEta[0]);
0078   const int nEBinsEta = (int)sizeof(eBinsEta)/sizeof(eBinsEta[0]);
0079   
0080   TGraphErrors *chi2E[5][2];
0081   //[5]: dir pho, pi0, eta, electron, pi+/-
0082   //[2]: vs. truth energy, vs. cluster energy
0083   TGraphErrors *probE[5][2];
0084 
0085   /*
0086   for(int i = 0; i < nTh3 - 1; i++)
0087     {
0088       gPad -> SetLogy(0);
0089       TH3F *hist = (TH3F*)fin -> Get(th3List[i].c_str());
0090       hist -> GetXaxis() -> SetTitle(th3xAxis[i].c_str());
0091       hist -> GetYaxis() -> SetTitle(th3yAxis[i].c_str());
0092       hist -> GetZaxis() -> SetTitle(th3zAxis[i].c_str());
0093       
0094       TH2F *hist2D = (TH2F*)hist -> Project3D("xy");
0095       c1 -> cd();
0096       hist2D -> Draw("colz");
0097       c1 -> SaveAs(Form("plots/%s_projxy.pdf",th3List[i].c_str()));
0098       
0099       TCanvas *cSec = new TCanvas();
0100 
0101       if(i > 2 && i < nTh3 - 1)
0102     {
0103       for(int j = 0; j < nEBins - 1; j++)
0104         {
0105           TH1F *hist1D = (TH1F*)hist2D->ProjectionY(Form("%s_truthE_%g_%g",th3List[i].c_str(), eBins[j], eBins[j+1]), hist2D -> GetXaxis() -> FindBin(eBins[j]), hist2D -> GetXaxis() -> FindBin(eBins[j+1]));
0106           //hist1D -> Fit("gaus",
0107           //gPad -> WaitPrimitive();
0108           hist1D -> Draw();
0109           TLatex energy;
0110           gPad -> SetLogy();
0111           energy.DrawLatexNDC(0.5,0.7,Form("%g < %s < %g",eBins[j], th3yAxis[i].c_str(), eBins[j+1]));
0112           cSec -> SaveAs(Form("plots/%s_ClusterE_%g_%g.pdf",th3List[i].c_str(), eBins[j], eBins[j+1]));
0113           
0114         }
0115     }
0116 
0117       hist2D = (TH2F*)hist -> Project3D("xz");
0118       c1 -> cd();
0119       hist2D -> Draw("colz");
0120       c1 -> SaveAs(Form("plots/%s_projxz.pdf",th3List[i].c_str()));
0121       
0122        if(i > 2 && i < nTh3 - 1)
0123     {
0124       for(int j = 0; j < nEBins - 1; j++)
0125         {
0126           TH1F *hist1D = (TH1F*)hist2D->ProjectionY(Form("%s_truthE_%g_%g",th3List[i].c_str(), eBins[j], eBins[j+1]), hist2D -> GetXaxis() -> FindBin(eBins[j]), hist2D -> GetXaxis() -> FindBin(eBins[j+1]));
0127           cSec -> cd();
0128           hist1D -> Draw();
0129           TLatex energy;
0130           gPad -> SetLogy();
0131           energy.DrawLatexNDC(0.5,0.7,Form("%g < %s < %g",eBins[j], th3zAxis[i].c_str(), eBins[j+1]));
0132           //hist1D -> Fit("gaus",
0133           //gPad -> WaitPrimitive();
0134           cSec -> SaveAs(Form("plots/%s_TruthE_%g_%g.pdf",th3List[i].c_str(), eBins[j], eBins[j+1]));
0135           
0136         }
0137     }
0138       
0139     }*/
0140   
0141   //going to handle the last TH3 separately since
0142   //it's going on a stupid big canvas. 
0143   
0144   TH3F *hist[nEtaBins]; 
0145   TH3F *histMatched[nEtaBins];
0146   float cw = 2*966, ch = 1.3*637;
0147   
0148   TCanvas *cPi0Cut[nEtaBins];
0149   TPad *pad_invMass[nEBins-1][nCutBins][nEtaBins];
0150 
0151   TCanvas *cSubMass[nEtaBins];
0152   TPad *pad_invMassSub[nEBins-1][nCutBins][nEtaBins];
0153   
0154   TCanvas *cMatchMass[nEtaBins];
0155   TPad *pad_matchMass[nEBins-1][nCutBins][nEtaBins];
0156   
0157 
0158   TGraphErrors *mass[nCutBins][nEtaBins];
0159   TGraphErrors *width[nCutBins][nEtaBins];
0160   TGraphErrors *yields[nCutBins][nEtaBins];
0161   TGraphErrors *yieldsMatched[nCutBins][nEtaBins];
0162   TGraphErrors *matchZeroFrac[nCutBins][nEtaBins];
0163   TF1 *invMassFit, *invMassFitBG;
0164   int colors[] = {1,2,4,kGreen+2, kViolet,kCyan,kOrange+2,kMagenta+2,kAzure-2};
0165   TLegend *leggy = new TLegend(.161, .729, .722, .939 );
0166   leggy -> SetNColumns(3);
0167   leggy -> SetFillStyle(0);
0168   leggy -> SetBorderSize(0);
0169   
0170   float xEndFit = 0.25;
0171   float xStartFit = 0.07;
0172   float xEndHist = 0.25;
0173   float xStartHist = 0.07;
0174   float invMassWindow[2] = {0.112, 0.162};
0175   float peakPos[nEBins - 1][nCutBins][nEtaBins];
0176   float peakW[nEBins - 1][nCutBins][nEtaBins];
0177   float sigma = 3.;
0178   TLatex marks;
0179   
0180   for(int k = 0; k < nEtaBins; k++)
0181     {
0182       cMatchMass[k] = new TCanvas(Form("cMatchMass_Eta%d",k),"Mass from Matched Cluster",cw,ch);
0183       cMatchMass[k] -> Range(0,0,1,1);
0184 
0185       cSubMass[k]= new TCanvas(Form("cMassSubPi0_Eta%d",k),"Subtracted Mass Dist",cw,ch);
0186       cSubMass[k] -> Range(0,0,1,1);
0187 
0188       cPi0Cut[k] = new TCanvas(Form("cMassPi0_Eta%d",k),"Invariant Mass",cw,ch);
0189       cPi0Cut[k] -> Range(0,0,1,1);
0190         
0191       hist[k] = (TH3F*)fin -> Get(Form("ePi0InvMassEcut_0Match_Eta%d",k));
0192       histMatched[k] = (TH3F*)fin -> Get(Form("ePi0InvMassEcut_1Match_Eta%d",k));
0193       
0194       for(int i = 0; i < nEBins - 1; i++)
0195     {
0196       for(int j = 0; j < nCutBins; j++)
0197         {
0198           if(i == 0)
0199         {
0200           mass[j][k] = new TGraphErrors();
0201           mass[j][k] -> SetMarkerColor(colors[j]);
0202 
0203           width[j][k] = new TGraphErrors();
0204           width[j][k] -> SetMarkerColor(colors[j]);
0205 
0206           yields[j][k] = new TGraphErrors();
0207           yields[j][k] -> SetMarkerColor(colors[j]);
0208 
0209           yieldsMatched[j][k] = new TGraphErrors();
0210           yieldsMatched[j][k] -> SetMarkerColor(colors[j]);
0211 
0212           matchZeroFrac[j][k] = new TGraphErrors();
0213           matchZeroFrac[j][k] -> SetMarkerColor(colors[j]);
0214         }
0215      
0216           pad_invMass[i][j][k] = getPad(i,j,k,cw,ch,nEBins , nCutBins);//for unsubtracted invariant mass distributions
0217           pad_invMass[i][j][k] -> cd();
0218 
0219           pad_invMassSub[i][j][k] = getPad(i,j,k,cw,ch,nEBins , nCutBins);//for subtracted mass distributions
0220 
0221           pad_matchMass[i][j][k] = getPad(i,j,k,cw,ch,nEBins, nCutBins);//for distributions composed of clusters matched to truth decay photons
0222 
0223           TH1F *hist1D = (TH1F*)hist[k] -> ProjectionY(Form("invMass_E%d_Cut%d_Eta%d",i,j,k), hist[k] -> GetXaxis() -> FindBin(eBins[i] + 0.001), hist[k] -> GetXaxis() -> FindBin(eBins[i+1] - 0.0001), hist[k] -> GetZaxis() -> FindBin(eCuts[j] + 0.0001), hist[k] -> GetNbinsZ());
0224       
0225           TH1F *hist1DMatched = (TH1F*)histMatched[k] -> ProjectionY(Form("invMass_E%d_Cut%d_eta%d_Matched",i,j,k), histMatched[k] -> GetXaxis() -> FindBin(eBins[i] + 0.001), histMatched[k] -> GetXaxis() -> FindBin(eBins[i+1] - 0.0001), histMatched[k] -> GetZaxis() -> FindBin(eCuts[j] + 0.0001), histMatched[k] -> GetNbinsZ());
0226       
0227       
0228           TH1F *invMassSub = (TH1F*)hist1D -> Clone();
0229           if(eBins[i+1] <= 13) 
0230         {
0231           invMassFit = new TF1(Form("invMassFit_E%d_Cut%d_Eta%d",i,j,k),"[0] + [1]*x + [2]*pow(x,2)+ [3]*pow(x,3)+  [4]*pow(x,4)+ [8]*pow(x,6)  +[5]*exp(-(([7]-x)^2)/(2*[6]^2))");
0232           
0233           
0234           invMassFitBG = new TF1(Form("invMassFit_E%d_Cut%d_Eta%dBG",i,j,k),"[0] + [1]*x + [2]*pow(x,2)+ [3]*pow(x,3)+  [4]*pow(x,4)+ [8]*pow(x,6)  +[5]*exp(-(([7]-x)^2)/(2*[6]^2)) - [5]*exp(-(([7]-x)^2)/(2*[6]^2))"); //hehe
0235           
0236           invMassFit -> FixParameter(8,0);
0237           
0238           invMassFit -> SetLineWidth(1);
0239           invMassFit -> SetParameter(6,0.015);
0240           
0241           invMassFit -> SetLineColor(2);
0242           invMassFit -> SetParameter(7,0.135);
0243           invMassFit -> SetParLimits(7, 0.135 - 0.07, 0.135 + 0.07);
0244           if((eCuts[j] >= 0.9) && (eBins[i] == 3 || eBins[i] == 4)) 
0245             {
0246               invMassFit -> SetParameter(5, hist1D -> GetMaximum()*2);
0247           
0248             }
0249           else invMassFit -> SetParameter(5,hist1D -> GetMaximum());
0250           hist1D -> Fit(invMassFit,"Q","",xStartFit, xEndFit);
0251           GetSubtractedDist(hist1D, invMassSub, invMassFit, invMassFitBG);
0252           invMassFitBG -> SetLineColor(4);
0253          
0254           float peakPosition = invMassFit -> GetParameter(7);
0255           peakPos[i][j][k] = peakPosition;
0256           
0257           float peakWidth = invMassFit -> GetParameter(6);
0258           peakW[i][j][k] = peakWidth;
0259           Double_t yieldErr, yieldErrMatched;
0260       
0261           float yield;
0262           if(hist1D -> Integral() == 0) yield = 0;
0263           else yield = invMassSub -> IntegralAndError(invMassSub -> GetXaxis() -> FindBin((peakPosition - sigma*peakWidth)+0.0001), invMassSub -> GetXaxis() -> FindBin((peakPosition + sigma*peakWidth)+0.0001), yieldErr, "");
0264 
0265           float yieldMatched = hist1DMatched -> IntegralAndError(hist1DMatched -> FindFirstBinAbove(0), hist1DMatched -> GetNbinsX(), yieldErrMatched, "");
0266           if(yield != 0)
0267             {
0268               mass[j][k] -> SetPoint(i, (eBins[i]+eBins[i+1])/2, invMassFit -> GetParameter(7));
0269               mass[j][k] -> SetPointError(i,  0, invMassFit -> GetParError(7));
0270     
0271               width[j][k] -> SetPoint(i, (eBins[i]+eBins[i+1])/2, invMassFit -> GetParameter(6));
0272               width[j][k] -> SetPointError(i, 0, invMassFit -> GetParError(6));
0273     
0274               yields[j][k] -> SetPoint(i, (eBins[i]+eBins[i+1])/2, yield);
0275               yields[j][k] -> SetPointError(i, 0, yieldErr);
0276               
0277              
0278             }
0279           yieldsMatched[j][k] -> SetPoint(i, (eBins[i]+eBins[i+1])/2, yieldMatched);
0280            
0281           yieldsMatched[j][k] -> SetPointError(i, 0, yieldErrMatched);
0282               
0283           float zeroMatchYield = hist1DMatched -> GetBinContent(hist1DMatched -> FindFirstBinAbove(0));
0284           float zeroMatchYieldErr = hist1DMatched -> GetBinError(46);
0285           float matchYieldRatio = zeroMatchYield/yieldMatched;
0286           float zeroMatchYieldRatErr = matchYieldRatio*sqrt(pow(yieldErrMatched/yieldMatched,2) + pow(zeroMatchYieldErr/zeroMatchYieldErr,2));
0287             
0288           matchZeroFrac[j][k] -> SetPoint(i, (eBins[i]+eBins[i+1])/2, matchYieldRatio);
0289           matchZeroFrac[j][k] -> SetPointError(i, 0, zeroMatchYieldRatErr);
0290         }
0291           hist1D -> GetXaxis() -> SetRangeUser(xStartHist,xEndHist);
0292       
0293           TCanvas *cMassFitCheck = new TCanvas();
0294           hist1D -> SetTitle(";M_{#gamma#gamma} [GeV/c^{2}];");
0295           hist1D -> Draw("epx0");
0296           invMassFitBG -> SetLineWidth(1);
0297 
0298           invMassFitBG -> Draw("same");
0299 
0300           marks.DrawLatexNDC(0.5,0.3,Form("#splitline{E_{Min} > %g GeV}{%g > E_{#pi^{0}} > %g}",eCuts[j], eBins[i], eBins[i+1]));
0301           cMassFitCheck -> SaveAs(Form("plots/invMassFitCheck_%dE_%dCut_Sigma%g_Eta%d.pdf",i,j,sigma,k));
0302 
0303 
0304           TCanvas *cMassMatchCheck = new TCanvas();
0305           hist1DMatched -> SetTitle(";M_{#gamma#gamma} [GeV/c^{2}];");
0306           hist1DMatched -> Draw("epx0");
0307           gPad -> SetLogy();
0308           marks.DrawLatexNDC(0.5,0.85,Form("#splitline{E_{Min} > %g GeV}{%g > E_{#pi^{0}} > %g}",eCuts[j], eBins[i], eBins[i+1]));
0309           cMassMatchCheck -> SaveAs(Form("plots/invMassMatchCheck_%dE_%dCut_Sigma%g_Eta%dMatch.pdf",i,j,sigma,k));
0310 
0311           drawCanvas_invMass(cPi0Cut[k], hist1D, j, i, pad_invMass[i][j][k], 0, peakPos[i][j][k], sigma*peakW[i][j][k]);
0312           invMassFitBG -> Draw("same");
0313           pad_invMassSub[i][j][k] -> cd();
0314 
0315           invMassSub -> GetXaxis() -> SetRangeUser(xStartHist, xEndHist);
0316           drawCanvas_invMass(cSubMass[k], invMassSub, j, i, pad_invMassSub[i][j][k], 0, peakPos[i][j][k], sigma*peakW[i][j][k]);
0317 
0318           pad_matchMass[i][j][k] -> cd();
0319           drawCanvas_invMass(cMatchMass[k], hist1DMatched, j, i, pad_matchMass[i][j][k], 1, peakPos[i][j][k], sigma*peakW[i][j][k]);
0320       
0321         }
0322     }
0323     
0324       cPi0Cut[k] -> SaveAs(Form("plots/invMassCuts_Sigma%g_Eta%d.pdf",sigma, k));
0325       cSubMass[k] -> SaveAs(Form("plots/invMassCutsSubtracted_Sigma%g_Eta%d.pdf",sigma, k));
0326       cMatchMass[k] -> SaveAs(Form("plots/invMassMatched_Sigma%g_Eta%d.pdf",sigma,k));
0327     }
0328   
0329   
0330   TCanvas *cMass[nEtaBins];
0331   for(int k = 0; k < nEtaBins; k++)
0332     {
0333       cMass[k] = new TCanvas();
0334       for(int i = 0; i < nCutBins; i++)
0335     {
0336       if(i == 0)
0337         {
0338           mass[i][k] -> Draw("ap");
0339           mass[i][k] -> GetYaxis() -> SetRangeUser(0.130,0.180);
0340           mass[i][k] -> GetXaxis() -> SetLimits(0,13);
0341           mass[i][k] -> GetXaxis() -> SetTitle("E_{#pi^{0}}^{Reco} [GeV]");
0342           mass[i][k] -> GetYaxis() -> SetTitle("Mass Peak Location [GeV/c^{2}]");
0343         }
0344       else mass[i][k] -> Draw("samep");
0345       if(k==0)leggy -> AddEntry(mass[i][k], Form("E_{min} > %g",eCuts[i]),"p");
0346     }
0347  
0348       leggy -> Draw("same");
0349       cMass[k] -> SaveAs(Form("plots/invMassGraphPeakPos_Sigma%g_Eta%d.pdf",sigma,k));
0350     }
0351 
0352   TCanvas *cWidth[nEtaBins];
0353   for(int k = 0; k < nEtaBins; k++)
0354     {
0355       cWidth[k] = new TCanvas();
0356 
0357       for(int i = 0; i < nCutBins; i++)
0358     {
0359       if(i == 0)
0360         {
0361           width[i][k] -> Draw("ap");
0362           width[i][k] -> GetYaxis() -> SetRangeUser(0.013, 0.028);
0363           width[i][k] -> GetXaxis() -> SetLimits(0,13);
0364           width[i][k] -> GetXaxis() -> SetTitle("E_{#pi^{0}}^{Reco} [GeV]");
0365           width[i][k] -> GetYaxis() -> SetTitle("Mass Peak Width [GeV/c^{2}]");
0366         }
0367       else width[i][k] -> Draw("samep");
0368     }
0369  
0370       leggy -> Draw("same");
0371       cWidth[k] -> SaveAs(Form("plots/invMassGraphPeakWidth_Sigma%g_Eta%d.pdf",sigma,k));
0372     }
0373 
0374   TCanvas *cYields[nEtaBins];
0375   for(int k = 0; k < nEtaBins; k++)
0376     {
0377       cYields[k] = new TCanvas();
0378       for(int i = 0; i < nCutBins; i++)
0379     {
0380       if(i == 0)
0381         {
0382           yields[i][k] -> Draw("ap");
0383           yields[i][k] -> GetYaxis() -> SetRangeUser(1e3,1e7); 
0384           yields[i][k] -> GetXaxis() -> SetTitle("E_{#pi^{0}}^{Reco} [GeV]");
0385           yields[i][k] -> GetYaxis() -> SetTitle("N_{#pi^{0}}");
0386         }
0387       else yields[i][k] -> Draw("samep");
0388       gPad -> SetLogy();
0389     }
0390   
0391       leggy -> Draw("same");
0392       cYields[k] -> SaveAs(Form("plots/invMassYield_Sigma%g_Eta%d.pdf",sigma,k));
0393     }
0394 
0395   TCanvas *cYieldsMatch[nEtaBins];
0396   for(int k = 0; k < nEtaBins; k++)
0397     {
0398       cYieldsMatch[k] = new TCanvas();
0399       for(int i = 0; i < nCutBins; i++)
0400     {
0401       if(i == 0)
0402         {
0403           yieldsMatched[i][k] -> Draw("ap");
0404           yieldsMatched[i][k] -> GetYaxis() -> SetRangeUser(1e3,1e7); 
0405           yieldsMatched[i][k] -> GetXaxis() -> SetTitle("E_{#pi^{0}}^{Reco} [GeV]");
0406           yieldsMatched[i][k] -> GetYaxis() -> SetTitle("N_{#pi^{0}}");
0407         }
0408       else yields[i][k] -> Draw("samep");
0409       gPad -> SetLogy();
0410     }
0411   
0412       leggy -> Draw("same");
0413       cYields[k] -> SaveAs(Form("plots/invMassYield_Sigma%g_Eta%d_Matched.pdf",sigma,k));
0414     }
0415   TH1F *pi0TruSpec[nEtaBins];
0416   //for(int i = 0; i < nEtaBins; i++) pi0TruSpec[i] = (TH1F*)fin -> Get(th1List[2].c_str());
0417   for(int i = 0; i < nEtaBins; i++) pi0TruSpec[i] = (TH1F*)fin -> Get(Form("truthPi0E_Eta%d",i));
0418   TGraphErrors *pi0Eff[nCutBins][nEtaBins];
0419   TCanvas *cEff[nEtaBins];
0420   float scaleMax;
0421   float scaleMin;
0422   
0423   for(int k = 0; k < nEtaBins; k++)
0424     {
0425       cEff[k] = new TCanvas();
0426       scaleMax = -9999;
0427       scaleMin = 9999;
0428       for(int i = 0; i < nCutBins ; i++)
0429     {
0430       pi0Eff[i][k] = new TGraphErrors();
0431       for(int j = 0; j < nEBins; j++)
0432         {
0433           Double_t pi0TruthYieldErr;
0434           float pi0TruthYield = pi0TruSpec[k] -> IntegralAndError(pi0TruSpec[k] -> GetXaxis() -> FindBin(eBins[j]+0.001), pi0TruSpec[k] -> GetXaxis() -> FindBin(eBins[j+1]-0.001), pi0TruthYieldErr, "");
0435       
0436       
0437           Double_t x, pi0YieldMeas, pi0YieldErr;
0438           yields[i][k] -> GetPoint(j, x, pi0YieldMeas);
0439           pi0YieldErr = yields[i][k] -> GetErrorY(j);
0440           float eff = pi0YieldMeas/pi0TruthYield;
0441           float err = eff*sqrt(pow(pi0YieldErr/pi0YieldMeas,2) + pow(pi0TruthYieldErr/pi0TruthYield,2));
0442           if(eff <  scaleMin) scaleMin = eff;
0443           if(eff > scaleMax) scaleMax = eff;
0444           pi0Eff[i][k] -> SetPoint(j, (eBins[j]+eBins[j+1])/2, eff);
0445           pi0Eff[i][k] -> SetPointError(j, 0, err);
0446         }
0447       pi0Eff[i][k] -> SetMarkerColor(colors[i]);
0448       if(i == 0)
0449         {
0450           pi0Eff[i][k] -> Draw("ap");
0451           pi0Eff[i][k] -> GetYaxis() -> SetRangeUser(-0.01,scaleMax*1.4);
0452           pi0Eff[i][k] -> GetXaxis() -> SetTitle("E_{#pi^{0}}^{Reco} [GeV]");
0453           pi0Eff[i][k] -> GetYaxis() -> SetTitle("#epsilon");
0454         }
0455       else pi0Eff[i][k] -> Draw("samep");
0456       
0457     }
0458       leggy -> Draw("same");
0459       cEff[k] -> SaveAs(Form("plots/pi0Eff_Sigma%g_Eta%d.pdf",sigma,k));
0460     }
0461 
0462   TGraphErrors *pi0EffMatch[nCutBins][nEtaBins];
0463   TCanvas *cEffMatch[nEtaBins];
0464   
0465   for(int k = 0; k < nEtaBins; k++)
0466     {
0467       scaleMax = -9999;
0468       scaleMin = 9999;
0469       cEffMatch[k] = new TCanvas();
0470       for(int i = 0; i < nCutBins ; i++)
0471     {
0472       pi0EffMatch[i][k] = new TGraphErrors();
0473       for(int j = 0; j < nEBins - 1; j++)
0474         {
0475           Double_t pi0TruthYieldErr;
0476           float pi0TruthYield = pi0TruSpec[k] -> IntegralAndError(pi0TruSpec[k] -> GetXaxis() -> FindBin(eBins[j]+0.001), pi0TruSpec [k]-> GetXaxis() -> FindBin(eBins[j+1]-0.001), pi0TruthYieldErr, "");
0477       
0478           Double_t x, pi0YieldMeas, pi0YieldErr;
0479           yieldsMatched[i][k] -> GetPoint(j, x, pi0YieldMeas);
0480           pi0YieldErr = yieldsMatched[i][k] -> GetErrorY(j);
0481           float eff = pi0YieldMeas/pi0TruthYield;
0482           float err = eff*sqrt(pow(pi0YieldErr/pi0YieldMeas,2) + pow(pi0TruthYieldErr/pi0TruthYield,2));
0483           if(eff <  scaleMin) scaleMin = eff;
0484           if(eff > scaleMax) scaleMax = eff;
0485           pi0EffMatch[i][k] -> SetPoint(j, (eBins[j]+eBins[j+1])/2, eff);
0486           pi0EffMatch[i][k] -> SetPointError(j, 0, err);
0487         }
0488       pi0EffMatch[i][k] -> SetMarkerColor(colors[i]);
0489       if(i == 0)
0490         {
0491           pi0EffMatch[i][k] -> Draw("ap");
0492           pi0EffMatch[i][k] -> GetYaxis() -> SetRangeUser(-0.01,scaleMax*1.25);
0493           pi0EffMatch[i][k] -> GetXaxis() -> SetTitle("E_{#pi^{0}}^{Reco} [GeV]");
0494           pi0EffMatch[i][k] -> GetYaxis() -> SetTitle("#epsilon");
0495         }
0496       else pi0EffMatch[i][k] -> Draw("samep");
0497       
0498     }
0499       leggy -> Draw("same");
0500       cEffMatch[k] -> SaveAs(Form("plots/pi0Eff_Sigma%g_Eta%d_Matched.pdf",sigma,k));
0501     }
0502 
0503    TCanvas *cZeroFrac[nEtaBins];
0504   for(int k = 0; k < nEtaBins; k++)
0505     {
0506       cZeroFrac[k] = new TCanvas();
0507       for(int j = 0; j < nCutBins; j++)
0508     {
0509       if(j == 0)
0510         {
0511           matchZeroFrac[j][k] -> Draw("ap");
0512           //matchZeroFrac[j][k] -> GetYaxis() -> SetRangeUser(-0.1,1.1);
0513           matchZeroFrac[j][k] -> SetTitle(";E_{#pi^{0}}^{Truth} [GeV]; N^{Unmatched}/N^{Matched}");
0514         }
0515       else  matchZeroFrac[j][k] -> Draw("samep");
0516     }
0517       leggy -> Draw("same");
0518       cZeroFrac[k] -> SaveAs(Form("plots/matchedZeroYieldFrac_Sigma%g_Eta%d.pdf",sigma,k));
0519     }
0520           
0521 
0522       
0523   TCanvas *cMiss = new TCanvas();
0524   cMiss -> Divide(4);
0525   cMiss -> cd();
0526   
0527   TH3F *unmatchedLocale = (TH3F*)fin -> Get("unmatchedLocale");
0528  
0529   float eBinsMerge[] = {1,4,5,9,13};
0530   
0531   for(int i = 0; i < 4; i++)
0532     {
0533       unmatchedLocale -> GetZaxis() -> SetRangeUser(eBinsMerge[i], eBinsMerge[i+1]);
0534       TH2F *locationProj = (TH2F*)unmatchedLocale -> Project3D("yx");
0535       locationProj -> SetTitle(Form("%g < E^{#pi^{0}}_{Truth} < %g;#eta;#phi [rad]", eBinsMerge[i], eBinsMerge[i+1]));
0536       locationProj -> Draw("colz");
0537       marks.DrawLatexNDC(0.2,0.5,Form("%g < E^{#pi^{0}}_{Truth} < %g", eBinsMerge[i], eBinsMerge[i+1]));
0538 
0539       cMiss -> SaveAs(Form("plots/missLocationBin_%d.pdf",i));
0540     }
0541 
0542   TCanvas *cMassScale = new TCanvas();
0543   TH3F *invMassPhoPi0 = (TH3F*)fin -> Get("invMassPhoPi0");
0544   
0545   float massZones[] = {0,0.100,0.200,1};
0546   
0547   for(int i = 0; i < 3; i++)
0548     {
0549       TH1F *pi0EScale = (TH1F*)invMassPhoPi0 -> ProjectionZ(Form("pi0EScale_%d",i),invMassPhoPi0 -> GetXaxis() -> FindBin(massZones[i] + 0.0001), invMassPhoPi0 -> GetXaxis() -> FindBin(massZones[i + 1] - 0.0001), 1, invMassPhoPi0 -> GetNbinsY());
0550       pi0EScale -> SetTitle(Form("%g < M_{#gamma#gamma} < %g; E_{Reco}^{#pi^{0}}/E_{True}^{#pi^{0}};",massZones[i], massZones[i+1]));
0551       pi0EScale -> Draw();
0552       marks.DrawLatexNDC(0.7,0.85,Form("%g < M_{#gamma#gamma} < %g",massZones[i], massZones[i+1]));
0553       cMassScale -> SaveAs(Form("plots/pi0EScale_MassBin_%d.pdf",i));
0554       
0555       
0556       TH1F *phoEScale = (TH1F*)invMassPhoPi0 -> ProjectionY(Form("phoEScale_%d",i),invMassPhoPi0 -> GetXaxis() -> FindBin(massZones[i] + 0.0001), invMassPhoPi0 -> GetXaxis() -> FindBin(massZones[i + 1] - 0.0001), 1, invMassPhoPi0 -> GetNbinsZ());
0557       phoEScale -> SetTitle(Form("%g < M_{#gamma#gamma} < %g; E_{Reco}^{#gamma}/E_{True}^{#gamma};",massZones[i], massZones[i+1]));
0558       phoEScale -> Draw();
0559       marks.DrawLatexNDC(0.7,0.85,Form("%g < M_{#gamma#gamma} < %g",massZones[i], massZones[i+1]));
0560 
0561       cMassScale -> SaveAs(Form("plots/phoEScale_MassBin_%d.pdf",i));
0562     }
0563   
0564   TCanvas *cMassPhi = new TCanvas();
0565   TH3F *invMassEPhi = (TH3F*)fin -> Get("invMassEPhi");
0566   
0567   for(int i = 0; i < 4; i++)
0568     {
0569       invMassEPhi -> GetYaxis() -> SetRangeUser(eBinsMerge[i],eBinsMerge[i+1]); 
0570       TH2F *invMassPhi = (TH2F*)invMassEPhi -> Project3D("xz");
0571       
0572       invMassPhi -> SetTitle(Form("%g <  E^{#pi^{0}}_{Truth} < %g;  #phi [rad]; M_{#gamma#gamma};",eBinsMerge[i], eBinsMerge[i+1] ));
0573       invMassPhi -> Draw("colz");
0574       marks.DrawLatexNDC(0.2,0.18,Form("%g <  E^{#pi^{0}}_{Truth} < %g",eBinsMerge[i], eBinsMerge[i+1] ));
0575 
0576       cMassPhi -> SaveAs(Form("plots/invMassPhi_EBin_%d.pdf",i));
0577       }
0578   
0579  
0580       
0581   
0582   /*
0583   TCanvas *cEtaCut = new TCanvas("cMassEta","Invariant Mass Eta Cuts",cw,ch);
0584   cEtaCut -> Range(0,0,1,1);
0585   TPad *pad_invMassEta[nEBinsEta - 1][nCutBinsEta];
0586     
0587   for(int i = 0; i < nEBinsEta - 1; i++)
0588     {
0589       for(int j = 0; j < nCutBinsEta; j++)
0590     {
0591       pad_invMassEta[i][j] = getPad(i,j,cw,ch,nEBinsEta , nCutBinsEta);
0592       pad_invMassEta[i][j] -> cd();
0593       TH1F *hist1D = (TH1F*)hist -> ProjectionY(Form("invMass_EtaE%d_Cut%d",i,j), hist -> GetXaxis() -> FindBin(eBinsEta[i] + 0.001), hist -> GetXaxis() -> FindBin(eBinsEta[i+1] - 0.0001), hist -> GetZaxis() -> FindBin(eCutsEta[j] + 0.0001), hist -> GetNbinsZ());
0594       drawCanvas_invMass(cEtaCut, hist1D, j, i, pad_invMassEta[i][j], 2);
0595     }
0596     }
0597   
0598     cEtaCut -> SaveAs("plots/invMassCutsEta.pdf");*/
0599            
0600 } 
0601   
0602 TPad* getPad(int j, int k, int l, float cw, float ch, const int nEBins, const int nCutBins)
0603 {
0604   float xpos[nCutBins+1], ypos[nEBins];
0605   int nCutOffset = 1;
0606   float marginLeft = 0.07;
0607   float marginRight = 0.035;
0608 
0609   float marginTop = 0.03;
0610   float marginBottom = 0.06;
0611   
0612   float marginColumn = 0.04;
0613   
0614   float w = (1 - marginLeft - marginRight - marginColumn*(nCutBins-2 + nCutOffset))/(nCutBins - 1+ nCutOffset);
0615   float h = (1 - marginTop - marginBottom)/(nEBins - 1);
0616 
0617   for(int m = 0; m < nCutBins + nCutOffset; m++)
0618     {
0619       if(m == 0) xpos[m] = 0;
0620       else if(m == nCutBins - 1+ nCutOffset) xpos[m] = 1.;
0621       else xpos[m] = marginLeft + m*w + (m-1)*marginColumn;
0622     }
0623   
0624   for(int m = 0; m < nEBins; m++)
0625     {
0626       if(m==0) ypos[nEBins -1 -m] = 0;
0627       else if (m == nEBins - 1) ypos[nEBins-1-m] = 1.;
0628       else ypos[nEBins - 1 -m] = marginBottom +m*h;
0629     }
0630   
0631   TPad *myPad = new TPad(Form("pad%d%d%d",j,k,l),Form("pad%d%d%d",j,k,l),xpos[k],ypos[j+1],xpos[k+1],ypos[j]);
0632   
0633   if(k==0) myPad->SetLeftMargin(marginLeft/(xpos[k+1] - xpos[k]));
0634   else myPad -> SetLeftMargin(marginColumn/(xpos[k+1] - xpos[k]));
0635   
0636   if(k == nCutBins - 2+ nCutOffset) myPad -> SetRightMargin(marginRight/(xpos[k+1]-xpos[k]));
0637   else myPad -> SetRightMargin(0.012);
0638   
0639   if(j == 0) myPad -> SetTopMargin(marginTop/(ypos[j] - ypos[j+1]));
0640   else myPad -> SetTopMargin(0);
0641   
0642   if(j == nEBins - 2) myPad -> SetBottomMargin(marginBottom/(ypos[j] - ypos[j+1]));
0643   else myPad -> SetBottomMargin(0);
0644   
0645   return myPad;
0646 }
0647   
0648 void drawCanvas_invMass(TCanvas *c, TH1F *hCorr, int pad_x, int pad_y, TPad *pad, int mode, float peakPos, float peakW)
0649 {
0650   
0651   float range = hCorr -> GetMaximum() - hCorr -> GetMinimum();
0652   SetHistoStyle(hCorr, pad_x, pad_y, hCorr -> GetMinimum(), hCorr->GetMaximum() + 0.2*range, mode);
0653   c -> cd(0);
0654   /*TLine *pi0Mass = new TLine(0.135,0,0.135,hCorr->GetMaximum());
0655   pi0Mass -> SetLineStyle(7);
0656   pi0Mass -> SetLineColor(2);
0657   TLine *etaMass = new TLine(0.548,0,0.548,hCorr->GetMaximum());
0658   etaMass -> SetLineStyle(7);
0659   etaMass -> SetLineColor(2);*/
0660   TLine *massWindow[2];
0661   //massWindow[0] = new TLine(0.112,0,0.112,hCorr->GetMaximum());
0662   pad -> Draw();
0663   pad -> cd();
0664   
0665   hCorr -> Draw("epX0");
0666   DrawCutLabel(pad_x,pad_y,mode);
0667   if(mode == 1) gPad -> SetLogy();
0668   if(peakPos)
0669     {
0670       massWindow[0] = new TLine(peakPos - peakW,0,peakPos - peakW,hCorr->GetMaximum());
0671       massWindow[0] -> SetLineStyle(7);
0672       massWindow[0] -> SetLineColor(2);
0673   
0674       //massWindow[1] = new TLine(0.162,0,0.162,hCorr->GetMaximum());
0675       massWindow[1] = new TLine(peakPos + peakW,0,peakPos + peakW,hCorr->GetMaximum());
0676       massWindow[1] -> SetLineStyle(7);
0677       massWindow[1] -> SetLineColor(2);
0678       massWindow[0] -> Draw("same");
0679       massWindow[1] -> Draw("same");
0680     }
0681   
0682   
0683   //etaMass -> Draw("same");
0684   //pi0Mass -> Draw("same");
0685   
0686 }
0687 
0688 void SetHistoStyle(TH1F *histo, int cutBin, int eBin, float low, float hi, int mode)
0689 {
0690   if(hi > low) histo -> GetYaxis() -> SetRangeUser(low,hi);
0691   if(mode == 1)
0692     {
0693       histo -> GetYaxis() -> SetRangeUser(0.1,10*hi);
0694     }
0695   
0696   histo -> SetTitle(";;");
0697   int xTitleBin[3][2] = {{11,3},{11,3},{3,3}};
0698   int yTitleBin[3][2] = {{5,0},{5,0},{2,0}};
0699 
0700   float yTitleOffset[3] = {0.3, 0.4, 1};
0701   float xTitleOffset[3] = {0.7, 0.7, 0.7};
0702 
0703   float yTitleSize[3] = {0.5, 0.5, 0.15};
0704   float xTitleSize[3] = {0.28, 0.28, 0.13}; 
0705 
0706   if(eBin == xTitleBin[mode][0] && cutBin == xTitleBin[mode][1]) histo -> GetXaxis() -> SetTitle("M_{#gamma#gamma}");
0707   if(eBin == yTitleBin[mode][0] && cutBin == yTitleBin[mode][1]) histo -> GetYaxis() -> SetTitle("N_{#gamma#gamma}");
0708   
0709   histo -> GetXaxis() -> SetTitleSize(xTitleSize[mode]);
0710   histo -> GetXaxis() -> SetTitleOffset(xTitleOffset[mode]);
0711   histo -> GetXaxis() -> CenterTitle();
0712   
0713   histo -> GetXaxis() -> SetLabelFont(43);
0714   histo -> GetXaxis() -> SetLabelSize(18);
0715   histo -> GetXaxis() -> SetLabelOffset(0.0);
0716   histo -> SetTickLength(0.05,"X");
0717   histo -> SetNdivisions(505,"x");
0718 
0719   
0720   histo -> GetYaxis() -> SetNdivisions(505);
0721   histo -> GetYaxis() -> SetTitleSize(yTitleSize[mode]);
0722   histo -> GetYaxis() -> SetTitleOffset(yTitleOffset[mode]);
0723   histo -> GetYaxis() -> CenterTitle();
0724   
0725   histo -> GetYaxis() -> SetLabelFont(43);
0726   histo -> GetYaxis() -> SetLabelSize(19);
0727   histo -> GetYaxis() -> SetLabelOffset(0.003);
0728   histo -> GetYaxis() -> SetNoExponent();
0729   histo -> SetTickLength(0.04,"Y");
0730   
0731   histo -> SetStats(0);
0732   //histo -> SetLineWidth(1.4);
0733   histo -> SetLineColor(1);
0734   histo -> SetMarkerColor(1);
0735   histo -> SetMarkerStyle(20);
0736   histo -> SetMarkerSize(0.4);
0737 }
0738 
0739 void DrawCutLabel(int x, int y, int mode)
0740 {
0741   gPad -> cd();
0742   gPad -> Update();
0743   float Y2 = gPad -> GetUymax();
0744   float Y1 = gPad -> GetUymin();
0745   float rangeY = Y2 - Y1;
0746   float xtop[3] = {0.8, -0.1, 0};
0747   float xbottom[3] = {0.273152, 1.01, 1.1};
0748   double ytop[3] = {Y2+0.01*rangeY, Y2+1*rangeY, Y2+0.01*rangeY};
0749   double ybottom[3] = {  0.98*Y2,  0.98*Y2,  0.9*Y2};
0750   
0751   TText *textCut;
0752   if(y == 0)
0753     {
0754       // if(mode < 2)
0755       //    {
0756       //      if(x <= 4)textCut = new TText(0.08, Y2+0.01*rangeY, Form("Min E > %g GeV", eCuts[x]));
0757       //      else textCut = new TText(0.08, Y2+0.01*rangeY, Form("Min E > %g GeV", eCuts[x]));
0758       //    }
0759       //else textCut = new TText(0, Y2+0.01*rangeY, Form("Min E > %g GeV", eCutsEta[x]));
0760       textCut = new TText(xtop[mode], ytop[mode], Form("Min E > %g GeV", eCuts[x]));
0761       textCut -> SetTextFont(43);
0762       textCut -> SetTextSize(18);
0763       textCut -> Draw();
0764     }
0765   if(x == sizeof(eCuts)/sizeof(eCuts[0]) - 1)
0766     {
0767       //if(mode < 2)textCut = new TText(0.273152, 0.98*Y2, Form("%g-%gGeV",eBins[y], eBins[y+1]));
0768       //else textCut = new TText(1.1, 0.9*Y2, Form("%g-%gGeV",eBinsEta[y], eBinsEta[y+1]));
0769       if(mode < 2)textCut = new TText(xbottom[mode], ybottom[mode], Form("%g-%gGeV",eBins[y], eBins[y+1]));
0770       else textCut = new TText(xbottom[mode], ybottom[mode], Form("%g-%gGeV",eBinsEta[y], eBinsEta[y+1]));
0771       textCut -> SetTextAngle(270);
0772       textCut -> SetTextFont(43);
0773       textCut -> SetTextSize(18);
0774       textCut -> Draw();
0775     }
0776   
0777 }
0778 
0779 void GetSubtractedDist(TH1F *histOrig, TH1F *histSub, TF1 *invMassFit, TF1 *invMassBG)
0780 {
0781   
0782  
0783   for(int i = 0; i < invMassFit -> GetNpar(); i++)
0784     {
0785       invMassBG -> SetParameter(i, invMassFit -> GetParameter(i));
0786     }
0787   
0788   for(int i = 1; i < histSub -> GetNbinsX() + 1; i++)
0789     {
0790       histSub -> SetBinContent(i, histOrig -> GetBinContent(i) - invMassBG -> Eval(histOrig -> GetBinCenter(i)));
0791     }
0792   
0793   
0794 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!