PWGLF/totEt/macros/emEt/hadCorr/PlotBackgroundClusters.C

   1 void SetStyles(TH1 *histo,int marker, int color,char *xtitle, char *ytitle);
   2 TH1* bayneseffdiv(TH1* numerator, TH1* denominator,Char_t* name) ;
   3 void PlotBackgroundClusters(TString filename="rootFiles/LHC11a10a_bis/Et.ESD.simPbPb.EMCal.LHC11a10a_bis.Run139465.root", Bool_t peripheral = kFALSE, int bin = 10, int binLast = 10, TString det = "EMCal"){
   4
   5   gStyle->SetOptTitle(0);
   6   gStyle->SetOptStat(0);
   7   gStyle->SetOptFit(0);
   8   TString outname = "";
   9   TString outnamebin = Form("%iTo%i",bin,binLast);
  10
  11   TString suffix = "";
  12   if(peripheral) suffix ="Peripheral";
  13
  14   TFile *f = TFile::Open(filename, "READ");
  15
  16   TList *l = (TList*)f->Get("out1");
  17   TH1F  *fHistChargedTracksCut = l->FindObject(("fHistChargedTracksCut"+suffix).Data());
  18   TH1F  *fHistChargedTracksAccepted = l->FindObject(("fHistChargedTracksAccepted"+suffix).Data());
  19   TH1F  *fHistGammasCut = l->FindObject(("fHistGammasTracksCut"+suffix).Data());
  20   TH1F  *fHistGammasAccepted = l->FindObject(("fHistGammasTracksAccepted"+suffix).Data());
  21
  22   int rebin = 1;
  23   fHistChargedTracksCut->Rebin(rebin);
  24   fHistChargedTracksAccepted->Rebin(rebin);
  25   fHistGammasCut->Rebin(rebin);
  26   fHistGammasAccepted->Rebin(rebin);
  27
  28   int bin1 = 1;
  29   int bin2 = fHistChargedTracksCut->GetXaxis()->FindBin(0.5);
  30   int bin3 = fHistChargedTracksCut->GetNbinsX();
  31   float nGammasCutBelow = fHistGammasCut->Integral(bin1,bin2);
  32   float nGammasCutAbove = fHistGammasCut->Integral(bin2+1,bin3);
  33   float nGammasCut = fHistGammasCut->Integral();
  34   float nTracksCutBelow = fHistChargedTracksCut->Integral(bin1,bin2);
  35   float nTracksCutAbove = fHistChargedTracksCut->Integral(bin2+1,bin3);
  36   float nTracksCut = fHistChargedTracksCut->Integral();
  37   cout<<"Gamma Fraction above: "<<nGammasCutAbove/nGammasCut<<" below: "<<nGammasCutBelow/nGammasCut<<endl;
  38   cout<<"Track Fraction above: "<<nTracksCutAbove/nTracksCut<<" below: "<<nTracksCutBelow/nTracksCut<<endl;
  39
  40   TH1F *hTotalCut = fHistChargedTracksCut->Clone("hTotalCut");
  41   hTotalCut->Add(fHistGammasCut);
  42   TH1F *hTotalAccepted = fHistChargedTracksAccepted->Clone("hTotalAccepted");
  43   hTotalAccepted->Add(fHistGammasAccepted);
  44   TH1F *hTotalClusters = fHistChargedTracksCut->Clone("hTotalClusters");
  45   hTotalClusters->Add(fHistChargedTracksAccepted);
  46   hTotalClusters->Add(fHistGammasCut);
  47   hTotalClusters->Add(fHistGammasAccepted);
  48   TH1F *hTotalCharged = fHistChargedTracksCut->Clone("hTotalCharged");
  49   hTotalCharged->Add(fHistChargedTracksAccepted);
  50   TH1F *hTotalGamma = fHistGammasCut->Clone("hTotalGamma");
  51   hTotalGamma->Add(fHistGammasAccepted);
  52
  53   TH1F *hSignalCut = bayneseffdiv(fHistGammasCut,hTotalCut,"hSignalCut");
  54   TH1F *hBkgdCut = bayneseffdiv(fHistChargedTracksCut,hTotalCut,"hBkgdCut");
  55   TH1F *hSignalAcc = bayneseffdiv(fHistGammasAccepted,hTotalAccepted,"hSignalAcc");
  56   TH1F *hBkgdAcc = bayneseffdiv(fHistChargedTracksAccepted,hTotalAccepted,"hBkgdAcc");
  57   TH1F *hFracCharged = bayneseffdiv(hTotalCharged,hTotalClusters,"hFracCharged");
  58   TH1F *hFracGamma = bayneseffdiv(hTotalGamma,hTotalClusters,"hFracGamma");
  59   SetStyles(hSignalCut,21,2,"E_{T}","fraction cut");
  60   SetStyles(hSignalAcc,25,2,"E_{T}","fraction cut");
  61   SetStyles(hBkgdCut,20,4,"E_{T}","fraction cut");
  62   SetStyles(hBkgdAcc,24,4,"E_{T}","fraction cut");
  63   SetStyles(hFracCharged,24,4,"E_{T}","fraction cut");
  64   SetStyles(hFracGamma,25,2,"E_{T}","fraction cut");
  65
  66   float ECut = 0.5;
  67   TLine *lineEDep = new TLine(ECut,0,ECut,1);
  68   lineEDep->SetLineWidth(2);
  69   TLine *lineEDep2D = new TLine(ECut,0,ECut,3);
  70   lineEDep2D->SetLineWidth(2);
  71
  72   TLegend *leg = new TLegend(0.357383,0.47043,0.458054,0.75);
  73   leg->SetFillStyle(0);
  74   leg->SetFillColor(0);
  75   leg->SetBorderSize(0);
  76   leg->SetTextSize(0.03);
  77   leg->AddEntry(hSignalCut,"Fraction of cut particles that are signals");
  78   leg->AddEntry(hSignalAcc,"Fraction of accepted particles that are signals");
  79   leg->AddEntry(hBkgdCut,"Fraction of cut particles that are bkgd");
  80   leg->AddEntry(hBkgdAcc,"Fraction of accepted particles that are background");
  81   leg->SetTextSize(0.061828);
  82
  83   TCanvas *c1 = new TCanvas("c1","c1",600,400);
  84   c1->SetTopMargin(0.02);
  85   c1->SetRightMargin(0.02);
  86   c1->SetBorderSize(0);
  87   c1->SetFillColor(0);
  88   c1->SetFillColor(0);
  89   c1->SetBorderMode(0);
  90   c1->SetFrameFillColor(0);
  91   c1->SetFrameBorderMode(0);
  92   hSignalCut->SetMaximum(1.0);
  93   hSignalCut->Draw();
  94   hSignalAcc->Draw("same");
  95   hBkgdCut->Draw("same");
  96   hBkgdAcc->Draw("same");
  97   lineEDep->Draw();
  98   leg->Draw();
  99
 100   c1->SaveAs("/tmp/SignalBkgdEmcal.png");
 101
 102   TLegend *leg1 = new TLegend(0.357383,0.47043,0.458054,0.75);
 103   leg1->SetFillStyle(0);
 104   leg1->SetFillColor(0);
 105   leg1->SetBorderSize(0);
 106   leg1->SetTextSize(0.03);
 107   leg1->AddEntry(hFracCharged,"Fraction of clusters from charged particles");
 108   leg1->AddEntry(hFracGamma,"Fraction of clusters from gammas");
 109   leg1->SetTextSize(0.061828);
 110
 111   TCanvas *c2 = new TCanvas("c2","c2",600,400);
 112   c2->SetTopMargin(0.02);
 113   c2->SetRightMargin(0.02);
 114   c2->SetBorderSize(0);
 115   c2->SetFillColor(0);
 116   c2->SetFillColor(0);
 117   c2->SetBorderMode(0);
 118   c2->SetFrameFillColor(0);
 119   c2->SetFrameBorderMode(0);
 120 //   hFracCharged->SetMaximum(1.0);
 121 //   hFracCharged->Draw();
 122 //   hFracGamma->Draw("same");
 123 //   leg1->Draw();
 124
 125 //   c2->SaveAs("/tmp/FractionsEmcal.png");
 126
 127
 128   TH2F  *fHistMatchedTracksEvspTBkgd = l->FindObject(("fHistMatchedTracksEvspTBkgd"+suffix).Data());
 129   TProfile * profBkgd = fHistMatchedTracksEvspTBkgd->ProfileX();
 130   TH2F  *fHistMatchedTracksEvspTSignal = l->FindObject(("fHistMatchedTracksEvspTSignal"+suffix).Data());
 131   TProfile * profSignal = fHistMatchedTracksEvspTSignal->ProfileX();
 132   TF1 *func = new TF1("func","x",0,3);
 133   fHistMatchedTracksEvspTBkgd->GetXaxis()->SetTitle("p_{T}");
 134   fHistMatchedTracksEvspTBkgd->GetYaxis()->SetTitle("E_{cluster}");
 135   fHistMatchedTracksEvspTSignal->GetXaxis()->SetTitle("p_{T}");
 136   fHistMatchedTracksEvspTSignal->GetYaxis()->SetTitle("E_{cluster}");
 137
 138
 139   TF1 *funcAvgSig = new TF1("funcAvgSig","[0]*x+[1]",0,3);
 140   funcAvgSig->SetParameter(0,1);
 141   funcAvgSig->SetParameter(1,0.3);
 142
 143   TF1 *funcCut = new TF1("funcCut","[0]*x+[1]",0,3);
 144   funcCut->SetParameter(0,9.54979e-01);
 145   funcCut->SetParameter(1,-2.47491e-01);
 146
 147
 148   TH3F  *fHistMatchedTracksEvspTBkgdvsMult = l->FindObject("fHistMatchedTracksEvspTBkgdvsCent");
 149   TH3F  *fHistMatchedTracksEvspTSignalvsMult = l->FindObject("fHistMatchedTracksEvspTSignalvsCent");
 150   //DoProjectProfile2D(const char* name, const char* title, TAxis* projX, TAxis* projY, bool originalRange, bool useUF, bool useOF) const
 151   fHistMatchedTracksEvspTBkgdvsMult->GetZaxis()->SetRange(bin,binLast);
 152   TH2D *hBkgd2D = (TProfile2D*) fHistMatchedTracksEvspTBkgdvsMult->Project3D("yx");
 153   TProfile * profBkgd2D = hBkgd2D->ProfileX();
 154   fHistMatchedTracksEvspTSignalvsMult->GetZaxis()->SetRange(bin,binLast);
 155   TH2D *hSignal2D = (TProfile2D*) fHistMatchedTracksEvspTSignalvsMult->Project3D("yx");
 156   TProfile * profSignal2D = hSignal2D->ProfileX();
 157
 158   TH2F  *fHistChargedTracksCutMult = l->FindObject("fHistChargedTracksCutMult");
 159   TH2F  *fHistChargedTracksAcceptedMult = l->FindObject("fHistChargedTracksAcceptedMult");
 160   TH2F  *fHistGammasCutMult = l->FindObject("fHistGammasTracksCutMult");
 161   TH2F  *fHistGammasAcceptedMult = l->FindObject("fHistGammasTracksAcceptedMult");
 162
 163   TH1D  *fHistChargedTracksCutPeri = fHistChargedTracksCutMult->ProjectionX("fHistChargedTracksCutPeri",bin,binLast);
 164   TH1D  *fHistChargedTracksAcceptedPeri = fHistChargedTracksAcceptedMult->ProjectionX("fHistChargedTracksAcceptedPeri",bin,binLast);
 165   TH1D  *fHistGammasCutPeri = fHistGammasCutMult->ProjectionX("fHistGammasTracksCutPeri",bin,binLast);
 166   TH1D  *fHistGammasAcceptedPeri = fHistGammasAcceptedMult->ProjectionX("fHistGammasTracksAcceptedPeri",bin,binLast);
 167
 168   int rebin = 1;
 169   fHistChargedTracksCutPeri->Rebin(rebin);
 170   fHistChargedTracksAcceptedPeri->Rebin(rebin);
 171   fHistGammasCutPeri->Rebin(rebin);
 172   fHistGammasAcceptedPeri->Rebin(rebin);
 173
 174   TH1D *hTotalCutPeri = fHistChargedTracksCutPeri->Clone("hTotalCutPeri");
 175   hTotalCutPeri->Add(fHistGammasCutPeri);
 176   TH1D *hTotalAcceptedPeri = fHistChargedTracksAcceptedPeri->Clone("hTotalAcceptedPeri");
 177   hTotalAcceptedPeri->Add(fHistGammasAcceptedPeri);
 178   TH1D *hTotalClustersPeri = fHistChargedTracksCutPeri->Clone("hTotalClustersPeri");
 179   hTotalClustersPeri->Add(fHistChargedTracksAcceptedPeri);
 180   hTotalClustersPeri->Add(fHistGammasCutPeri);
 181   hTotalClustersPeri->Add(fHistGammasAcceptedPeri);
 182   TH1D *hTotalChargedPeri = fHistChargedTracksCutPeri->Clone("hTotalChargedPeri");
 183   hTotalChargedPeri->Add(fHistChargedTracksAcceptedPeri);
 184   TH1D *hTotalGammaPeri = fHistGammasCutPeri->Clone("hTotalGammaPeri");
 185   hTotalGammaPeri->Add(fHistGammasAcceptedPeri);
 186
 187   TH1D *hSignalCutPeri = bayneseffdiv(fHistGammasCutPeri,hTotalCutPeri,"hSignalCutPeri");
 188   TH1D *hBkgdCutPeri = bayneseffdiv(fHistChargedTracksCutPeri,hTotalCutPeri,"hBkgdCutPeri");
 189   TH1D *hSignalAccPeri = bayneseffdiv(fHistGammasAcceptedPeri,hTotalAcceptedPeri,"hSignalAccPeri");
 190   TH1D *hBkgdAccPeri = bayneseffdiv(fHistChargedTracksAcceptedPeri,hTotalAcceptedPeri,"hBkgdAccPeri");
 191   TH1D *hFracChargedPeri = bayneseffdiv(hTotalChargedPeri,hTotalClustersPeri,"hFracChargedPeri");
 192   TH1D *hFracGammaPeri = bayneseffdiv(hTotalGammaPeri,hTotalClustersPeri,"hFracGammaPeri");
 193   SetStyles(hSignalCutPeri,21,2,"E_{T}","fraction cut");
 194   SetStyles(hSignalAccPeri,25,2,"E_{T}","fraction cut");
 195   SetStyles(hBkgdCutPeri,20,4,"E_{T}","fraction cut");
 196   SetStyles(hBkgdAccPeri,24,4,"E_{T}","fraction cut");
 197   SetStyles(hFracChargedPeri,24,4,"E_{T}","fraction cut");
 198   SetStyles(hFracGammaPeri,25,2,"E_{T}","fraction cut");
 199
 200   TLegend *legPeri = new TLegend(0.357383,0.47043,0.458054,0.75);
 201   legPeri->SetFillStyle(0);
 202   legPeri->SetFillColor(0);
 203   legPeri->SetBorderSize(0);
 204   legPeri->SetTextSize(0.03);
 205   legPeri->AddEntry(hSignalCutPeri,"Fraction of cut particles that are signals");
 206   legPeri->AddEntry(hSignalAccPeri,"Fraction of accepted particles that are signals");
 207   legPeri->AddEntry(hBkgdCutPeri,"Fraction of cut particles that are bkgd");
 208   legPeri->AddEntry(hBkgdAccPeri,"Fraction of accepted particles that are background");
 209   legPeri->SetTextSize(0.061828);
 210
 211   //profSignal->Fit(funcAvgSig,"","",ECut,3);
 212   funcAvgSig->SetLineStyle(2);
 213   funcAvgSig->SetLineColor(1);
 214   funcCut->SetLineStyle(2);
 215   funcCut->SetLineColor(1);
 216
 217   TCanvas *c3 = new TCanvas("c3","Signal All",600,400);
 218   c3->SetTopMargin(0.02);
 219   c3->SetRightMargin(0.02);
 220   c3->SetBorderSize(0);
 221   c3->SetFillColor(0);
 222   c3->SetFillColor(0);
 223   c3->SetBorderMode(0);
 224   c3->SetFrameFillColor(0);
 225   c3->SetFrameBorderMode(0);
 226   fHistMatchedTracksEvspTSignal->Draw("colz");
 227   profSignal->Draw("esame");
 228   profBkgd->Draw("esame");
 229   func->Draw("same");
 230   //funcCut->Draw("same");
 231   lineEDep2D->Draw();
 232
 233
 234   TCanvas *c10 = new TCanvas("c10","c10",600,400);
 235   c10->SetTopMargin(0.02);
 236   c10->SetRightMargin(0.02);
 237   c10->SetBorderSize(0);
 238   c10->SetFillColor(0);
 239   c10->SetFillColor(0);
 240   c10->SetBorderMode(0);
 241   c10->SetFrameFillColor(0);
 242   c10->SetFrameBorderMode(0);
 243   int firstbin = fHistMatchedTracksEvspTBkgd->GetXaxis()->FindBin(0.75);
 244   TProfile * profE = fHistMatchedTracksEvspTBkgd->ProfileY("Test",firstbin,firstbin+1);
 245   profE->Draw();
 246
 247   TCanvas *c4 = new TCanvas("c4","Background All",600,400);
 248   c4->SetTopMargin(0.02);
 249   c4->SetRightMargin(0.02);
 250   c4->SetBorderSize(0);
 251   c4->SetFillColor(0);
 252   c4->SetFillColor(0);
 253   c4->SetBorderMode(0);
 254   c4->SetFrameFillColor(0);
 255   c4->SetFrameBorderMode(0);
 256   fHistMatchedTracksEvspTBkgd->Draw("colz");
 257   profSignal->Draw("esame");
 258   profBkgd->Draw("esame");
 259   func->Draw("same");
 260   //funcCut->Draw("same");
 261   //funcAvgSig->Draw("same");
 262   lineEDep2D->Draw();
 263
 264   TCanvas *c7 = new TCanvas("c7","Background Binned",600,400);
 265   c7->SetTopMargin(0.02);
 266   c7->SetRightMargin(0.02);
 267   c7->SetBorderSize(0);
 268   c7->SetFillColor(0);
 269   c7->SetFillColor(0);
 270   c7->SetBorderMode(0);
 271   c7->SetFrameFillColor(0);
 272   c7->SetFrameBorderMode(0);
 273   hBkgd2D->Draw("colz");
 274   profBkgd2D->Draw("same");
 275   func->Draw("same");
 276   lineEDep2D->Draw();
 277   //funcCut->Draw("same");
 278   //funcAvgSig->Draw("same");
 279   //c7->SaveAs(Form("/tmp/Bkgd%i.png",bin));
 280   outname = "/tmp/TrackMatchingBkgd2D"+det+outnamebin+".png";
 281   c7->SaveAs(outname.Data());
 282
 283   TCanvas *c8 = new TCanvas("c8","Signal Binned",600,400);
 284   c8->SetTopMargin(0.02);
 285   c8->SetRightMargin(0.02);
 286   c8->SetBorderSize(0);
 287   c8->SetFillColor(0);
 288   c8->SetFillColor(0);
 289   c8->SetBorderMode(0);
 290   c8->SetFrameFillColor(0);
 291   c8->SetFrameBorderMode(0);
 292   hSignal2D->Draw("colz");
 293   profSignal2D->Draw("same");
 294   func->Draw("same");
 295   lineEDep2D->Draw();
 296   //funcCut->Draw("same");
 297   //funcAvgSig->Draw("same");
 298   //c8->SaveAs(Form("/tmp/Signal%i.png",bin));
 299   outname = "/tmp/TrackMatchingSignal2D"+det+outnamebin+".png";
 300   c8->SaveAs(outname.Data());
 301
 302
 303   TCanvas *c9 = new TCanvas("c9","c9",600,400);
 304   c9->SetTopMargin(0.02);
 305   c9->SetRightMargin(0.02);
 306   c9->SetBorderSize(0);
 307   c9->SetFillColor(0);
 308   c9->SetFillColor(0);
 309   c9->SetBorderMode(0);
 310   c9->SetFrameFillColor(0);
 311   c9->SetFrameBorderMode(0);
 312   hSignalCutPeri->SetMaximum(1.0);
 313   hSignalCutPeri->Draw();
 314   hSignalAccPeri->Draw("same");
 315   hBkgdCutPeri->Draw("same");
 316   hBkgdAccPeri->Draw("same");
 317   lineEDep->Draw();
 318   legPeri->Draw();
 319   outname = "/tmp/TrackMatchingForCuts"+det+outnamebin+".png";
 320   c9->SaveAs(outname.Data());
 321   return;
 322
 323   TCanvas *c5 = new TCanvas("c5","Signal Binned",600,400);
 324   c5->SetTopMargin(0.02);
 325   c5->SetRightMargin(0.02);
 326   c5->SetBorderSize(0);
 327   c5->SetFillColor(0);
 328   c5->SetFillColor(0);
 329   c5->SetBorderMode(0);
 330   c5->SetFrameFillColor(0);
 331   c5->SetFrameBorderMode(0);
 332   TH2F *hSignalOverBkgd = (TH2F*) fHistMatchedTracksEvspTSignal->Clone("SignalOverBackground");
 333   hSignalOverBkgd->Divide(fHistMatchedTracksEvspTBkgd);
 334   hSignalOverBkgd->Draw("colz");
 335   func->Draw("same");
 336   profSignal->Draw("esame");
 337   profBkgd->Draw("esame");
 338
 339   TCanvas *c6 = new TCanvas("c6","c6",600,400);
 340   c6->SetTopMargin(0.02);
 341   c6->SetRightMargin(0.02);
 342   c6->SetBorderSize(0);
 343   c6->SetFillColor(0);
 344   c6->SetFillColor(0);
 345   c6->SetBorderMode(0);
 346   c6->SetFrameFillColor(0);
 347   c6->SetFrameBorderMode(0);
 348   TH2F *hBkgdOverSignal = (TH2F*) fHistMatchedTracksEvspTBkgd->Clone("BackgroundOverSignal");
 349   hBkgdOverSignal->Divide(fHistMatchedTracksEvspTSignal);
 350   hBkgdOverSignal->Draw("colz");
 351   func->Draw("same");
 352   profSignal->Draw("esame");
 353   profBkgd->Draw("esame");
 354
 355   //c3->SaveAs("/tmp/SignalBkgdEmcal.png");
 356
 357 }
 358
 359 void SetStyles(TH1 *histo,int marker, int color,char *xtitle, char *ytitle){
 360   histo->SetMarkerStyle(marker);
 361   histo->SetMarkerColor(color);
 362   histo->SetLineColor(color);
 363   histo->GetXaxis()->SetTitle(xtitle);
 364   histo->GetYaxis()->SetTitle(ytitle);
 365 }
 366
 367
 368 TH1* bayneseffdiv(TH1* numerator, TH1* denominator,Char_t* name)
 369 {
 370     if(!numerator){
 371       cerr<<"Error:  numerator does not exist!"<<endl;
 372       return NULL;
 373     }
 374     if(!denominator){
 375       cerr<<"Error:  denominator does not exist!"<<endl;
 376       return NULL;
 377     }
 378     TH1* result = (TH1*)numerator->Clone(name);
 379     Int_t nbins = numerator->GetNbinsX();
 380     for (Int_t ibin=0; ibin<= nbins+1; ++ibin) {
 381       Double_t numeratorVal = numerator->GetBinContent(ibin);
 382       Double_t denominatorVal = denominator->GetBinContent(ibin);
 383       // Check if the errors are right or the thing is scaled
 384       Double_t numeratorValErr = numerator->GetBinError(ibin);
 385       if (!(numeratorValErr==0. || numeratorVal ==0.) ) {
 386         Double_t rescale = numeratorValErr*numeratorValErr/numeratorVal;
 387         numeratorVal /= rescale;
 388       }
 389       Double_t denominatorValErr = denominator->GetBinError(ibin);
 390       if (!(denominatorValErr==0. || denominatorVal==0. )) {
 391         Double_t rescale = denominatorValErr*denominatorValErr/denominatorVal;
 392         denominatorVal /= rescale;
 393       }
 394       Double_t quotient = 0.;
 395       if (denominatorVal!=0.) {
 396         quotient = numeratorVal/denominatorVal;
 397       }
 398       Double_t quotientErr=0;
 399       quotientErr = TMath::Sqrt(
 400                                 (numeratorVal+1.0)/(denominatorVal+2.0)*
 401                                 ((numeratorVal+2.0)/(denominatorVal+3.0)-(numeratorVal+1.0)/(denominatorVal+2.0)));
 402       result->SetBinContent(ibin,quotient);
 403       result->SetBinError(ibin,quotientErr);
 404       //cout<<"Setting bin "<<ibin<<" to "<<quotient<<" "<<numeratorVal<<"/"<<denominatorVal<<endl;
 405     }
 406     return result;
 407 }