PWGDQ/dielectron/macros/PlotDataResults.C

   1 #include <TROOT.h>
   2 #include <TStyle.h>
   3 #include <TH1D.h>
   4 #include <TObjArray.h>
   5 #include <TCanvas.h>
   6 #include <TMath.h>
   7 #include <TAxis.h>
   8 #include <TLatex.h>
   9 #include <TFile.h>
  10 #include <TString.h>
  11 #include <TLegend.h>
  12 #include <TList.h>
  13
  14 #include "AliDielectronSignalExt.h"
  15 #include "AliDielectronCFdraw.h"
  16 #include "AliDielectron.h"
  17
  18 AliDielectronSignalBase* GetSignalLS(AliDielectronCFdraw &d, Int_t step,
  19                                      AliDielectronSignalBase::EBackgroundMethod type=AliDielectronSignalBase::kLikeSign);
  20 AliDielectronSignalBase* GetSignalRot(AliDielectronCFdraw &d, Int_t step);
  21 void SetStyle(AliDielectronSignalBase *sig, const char* nameAdd);
  22 void DrawSpectra(AliDielectronSignalBase *sig, const char* cname, TH1  *hEventStat=0x0, Bool_t save=kFALSE);
  23 void FitMCshape(AliDielectronSignalBase *sig);
  24
  25 const char *mcLineShapeFile="$ALICE_ROOT/PWGDQ/dielectron/macros/mcMinv_LHC10f7a.root";
  26 TString addToName="";
  27
  28 //_______________________________________
  29 void PlotDataResults(const char* filenameData, const char* filenameMC="", Bool_t save=kFALSE)
  30 {
  31   if (!addToName.IsNull()) addToName.Prepend("-");
  32   AliDielectronCFdraw d(filenameData);
  33   AliDielectronCFdraw dCorr("corrCont","corrCont");
  34   TString nameCorr(filenameMC);
  35   if (!nameCorr.IsNull()) d.SetCFContainers(nameCorr.Data());
  36   TFile f(filenameData);
  37   TH1 *hStats=(TH1*)f.Get("hEventStat");
  38   if (!f.IsOpen() || f.IsZombie() || !hStats) return;
  39   hStats->SetDirectory(0);
  40   f.Close();
  41
  42   Int_t stepFirst=0, stepAny=1, stepTOFmix=2;
  43
  44   gStyle->SetOptStat(0);
  45   //Set common Ranges
  46   d.SetRangeUser("Leg1_NclsTPC",70.,170.);
  47   d.SetRangeUser("Leg2_NclsTPC",70.,170.);
  48   d.SetRangeUser("Leg1_Pt",1.01,100000);
  49   d.SetRangeUser("Leg2_Pt",1.01,100000);
  50   d.SetRangeUser("Leg1_Eta",-0.899,0.899);
  51   d.SetRangeUser("Leg2_Eta",-0.899,0.899);
  52
  53   d.SetRangeUser("Leg1_TPC_nSigma_Electrons",-3.,2.99);
  54   d.SetRangeUser("Leg2_TPC_nSigma_Electrons",-3.,2.99);
  55   d.SetRangeUser("Leg1_TPC_nSigma_Pions",3.51,20);
  56   d.SetRangeUser("Leg2_TPC_nSigma_Pions",3.51,20);
  57   d.SetRangeUser("Leg1_TPC_nSigma_Protons",3.01,20);
  58   d.SetRangeUser("Leg2_TPC_nSigma_Protons",3.01,20);
  59
  60 //   d.SetRangeUser("Pt",0,1000);
  61
  62   d.SetRangeUser("M",0.5,5.);
  63   //============================
  64   //SPD first
  65   //
  66
  67   //--- Like sign subtraction
  68   AliDielectronSignalBase *sigFirst=GetSignalLS(d,stepFirst);
  69   SetStyle(sigFirst,"ITS First - Like Sign subtraction");
  70   DrawSpectra(sigFirst,"cFirst",hStats,save);
  71   //--- Like sign subtraction Arithmetic mean
  72   AliDielectronSignalBase *sigFirstArith=GetSignalLS(d,stepFirst,AliDielectronSignalBase::kLikeSignArithm);
  73   SetStyle(sigFirstArith,"ITS FirstArith - Like Sign subtraction");
  74   DrawSpectra(sigFirstArith,"cFirstArith",hStats,save);
  75
  76   //============================
  77   //SPD any
  78   //
  79   AliDielectronSignalBase *sigAny=GetSignalLS(d,stepAny);
  80   SetStyle(sigAny,"ITS Any - Like Sign subtraction");
  81   DrawSpectra(sigAny,"cAny",hStats,save);
  82   //--- like sign with arithmetic mean
  83   AliDielectronSignalBase *sigAnyArith=GetSignalLS(d,stepAny,AliDielectronSignalBase::kLikeSignArithm);
  84   SetStyle(sigAnyArith,"ITS Any - Like Sign subtraction (Arithm. mean)");
  85   DrawSpectra(sigAnyArith,"cAnyArith",hStats,save);
  86
  87   if (hStats) delete hStats;
  88 }
  89
  90
  91 //_______________________________________
  92 AliDielectronSignalBase *GetSignalLS(AliDielectronCFdraw &d, Int_t step, AliDielectronSignalBase::EBackgroundMethod type)
  93 {
  94   //
  95   // Get Extracted signal from likesign method
  96   //
  97
  98   TObjArray *arr=new TObjArray;
  99   arr->SetOwner();
 100
 101   for (Int_t iType=0;iType<3;++iType){
 102     d.SetRangeUser("PairType",iType,iType);
 103     arr->AddAt(d.Project("M",step),iType);
 104   }
 105
 106   AliDielectronSignalExt *sig=new AliDielectronSignalExt;
 107   sig->SetScaleRawToBackground(3.2,4.9);
 108   sig->SetIntegralRange(2.92,3.15);
 109   sig->SetMethod(type);
 110   sig->Process(arr);
 111
 112   delete arr;
 113   return sig;
 114 }
 115
 116 //_______________________________________
 117 AliDielectronSignalBase *GetSignalRot(AliDielectronCFdraw &d, Int_t step)
 118 {
 119   //
 120   // Get Extracted signal from likesign method
 121   //
 122
 123   TObjArray *arr=new TObjArray;
 124   arr->SetOwner();
 125
 126   Int_t iType=AliDielectron::kEv1PM;
 127   d.SetRangeUser("PairType",iType,iType);
 128   arr->AddAt(d.Project("M",step),iType);
 129
 130   iType=AliDielectron::kEv1PMRot;
 131   d.SetRangeUser("PairType",iType,iType);
 132   arr->AddAt(d.Project("M",step),iType);
 133
 134   AliDielectronSignalExt *sig=new AliDielectronSignalExt;
 135   sig->SetScaleRawToBackground(3.2,4.9);
 136   sig->SetIntegralRange(2.93,3.15);
 137   sig->SetMethod(AliDielectronSignalBase::kRotation);
 138   sig->Process(arr);
 139
 140
 141   delete arr;
 142   return sig;
 143 }
 144
 145 //_______________________________________
 146 void SetStyle(AliDielectronSignalBase *sig, const char* nameAdd)
 147 {
 148   //
 149   //
 150   //
 151   TH1 *hUS=sig->GetUnlikeSignHistogram();
 152   hUS->SetMarkerStyle(20);
 153   hUS->SetMarkerSize(0.7);
 154   hUS->SetMarkerColor(kRed);
 155   hUS->SetLineColor(kRed);
 156   hUS->SetStats(0);
 157   hUS->SetTitle(Form("Like sign spectrum %s;M inv. ee;counts per %.4g MeV",
 158                      nameAdd,1000*hUS->GetXaxis()->GetBinWidth(1)));
 159
 160   TH1* hBackground=sig->GetBackgroundHistogram();
 161   hBackground->SetMarkerStyle(24);
 162   hBackground->SetMarkerSize(0.7);
 163   hBackground->SetStats(0);
 164   hBackground->SetMarkerColor(kBlue);
 165   hBackground->SetLineColor(kBlue);
 166   hBackground->SetTitle(Form("Like sign spectrum %s;M inv. ee;counts per %.4g MeV",
 167                              nameAdd,1000*hBackground->GetXaxis()->GetBinWidth(1)));
 168
 169   TH1* hSignal=sig->GetSignalHistogram();
 170   hSignal->SetMarkerStyle(20);
 171   hSignal->SetMarkerSize(0.7);
 172   hSignal->SetMarkerColor(kRed);
 173   hSignal->SetLineColor(kRed);
 174   hSignal->SetTitle(Form("Like sign spectrum %s;M inv. ee;counts per %.4g MeV",
 175                          nameAdd,1000*hSignal->GetXaxis()->GetBinWidth(1)));
 176
 177
 178 }
 179
 180 //_______________________________________
 181 void DrawSpectra(AliDielectronSignalBase *sig, const char* cname, TH1  *hEventStat, Bool_t save)
 182 {
 183   //
 184   //
 185   //
 186
 187
 188   FitMCshape(sig);
 189
 190   gStyle->SetOptTitle(0);
 191   TCanvas *c=(TCanvas*)gROOT->FindObject(cname);
 192   if (!c) c=new TCanvas(cname,cname,400*1.3,500*1.3);
 193   c->SetTopMargin(0.04);  c->SetRightMargin(0.04);
 194   c->SetLeftMargin(0.13); c->SetBottomMargin(0.14);
 195   c->Clear();
 196   c->Divide(1,2,0,0);
 197   c->cd(1);
 198   gPad->SetTopMargin(0.01);
 199   gPad->SetRightMargin(0.01);
 200   gPad->SetLeftMargin(0.15);
 201   gPad->SetBottomMargin(0.001);
 202   /*
 203   gPad->SetGridx();
 204   gPad->SetGridy();
 205 */
 206   gPad->SetTickx();
 207   gPad->SetTicky();
 208
 209   TH1 *hUS=sig->GetUnlikeSignHistogram();
 210
 211   TH1* hBackground=sig->GetBackgroundHistogram();
 212
 213   hUS->GetXaxis()->CenterTitle();
 214   hUS->GetXaxis()->SetTitleSize(0.07); hUS->GetXaxis()->SetLabelSize(0.059);
 215 //   hUS->GetXaxis()->SetLabelOffset(1.);
 216   hUS->GetXaxis()->SetTitleOffset(1.1);
 217   hUS->GetXaxis()->SetNdivisions(510);
 218   hUS->GetYaxis()->CenterTitle();
 219   hUS->GetYaxis()->SetTitleSize(0.07); hUS->GetYaxis()->SetLabelSize(0.059);
 220   hUS->GetYaxis()->SetTitleOffset(1.);
 221   hUS->GetXaxis()->SetLabelFont(42); hUS->GetYaxis()->SetLabelFont(42);
 222   hUS->GetXaxis()->SetTitleFont(42); hUS->GetYaxis()->SetTitleFont(42);
 223   hUS->GetYaxis()->SetNdivisions(510);
 224   //hUS->GetXaxis()->SetLimits(mMinPlot,mMmaxPlot);
 225   //hUS->SetLineWidth(.7);
 226   hUS->SetMarkerSize(.8);
 227   hUS->SetMarkerColor(2); hUS->SetMarkerStyle(20); hUS->SetLineColor(2);
 228   hUS->Draw();
 229   hUS->SetMaximum(hUS->GetMaximum()*1.3);
 230
 231   //hBackground->SetLineWidth(.4);
 232   hBackground->SetMarkerSize(.7);
 233   hBackground->SetMarkerColor(4); hBackground->SetMarkerStyle(27); hBackground->SetLineColor(4);
 234
 235   hBackground->Draw("same");
 236
 237   TLatex *lat=new TLatex;
 238   lat->SetNDC(kTRUE);
 239   //  lat->DrawLatex(0.68, 0.67, "ALICE Performance");
 240   lat->SetTextColor(1);lat->SetTextFont(42);lat->SetTextSize(.045);
 241
 242   Double_t sigN=sig->GetSignal();
 243   Double_t sigEr=sig->GetSignalError();
 244   Double_t sigS2B=sig->GetSB();
 245   Double_t sigS2Ber=sig->GetSBError();
 246   Double_t sigSignif= sig->GetSignificance();
 247   Double_t sigSignifEr= sig->GetSignificanceError();
 248   lat->DrawLatex(0.18, 0.92, Form("S: %3d#pm%4.1f, S/B: %3.1f#pm %4.2f, Signif.: %4.1f#pm%4.2f (%4.2f-%4.2f GeV) ",(int)sigN,sigEr,sigS2B,sigS2Ber,sigSignif,sigSignifEr,
 249                        hUS->GetBinLowEdge(hUS->FindBin(sig->GetIntegralMin())),
 250                        hUS->GetBinLowEdge(hUS->FindBin(sig->GetIntegralMax())+1)));
 251
 252   TLegend *leg=new TLegend(0.17,0.72,0.42,0.88);
 253   leg->SetBorderSize(0); leg->SetFillColor(0); leg->SetTextFont(42);
 254   leg->SetFillStyle(0); leg->SetMargin(0.25); //separation symbol-text
 255   leg->SetEntrySeparation(0.15);
 256   leg->AddEntry(hUS,"OS", "p");
 257   leg->AddEntry(hBackground, Form("LS*%4.2f",sig->GetScaleFactor()), "p");
 258   leg->Draw();
 259
 260   c->cd(2);
 261   gPad->SetRightMargin(0.01);
 262   gPad->SetTopMargin(0);
 263   gPad->SetLeftMargin(0.15);
 264   gPad->SetBottomMargin(.15);
 265   /*
 266   gPad->SetGridx();
 267   gPad->SetGridy();
 268   */
 269   gPad->SetTickx();
 270   gPad->SetTicky();
 271
 272   TH1* hSignal=sig->GetSignalHistogram();
 273   hSignal->GetXaxis()->SetTitle("M_{ee} (GeV/c^{2})");
 274   hSignal->GetXaxis()->CenterTitle();
 275   hSignal->GetXaxis()->SetTitleSize(0.06); hSignal->GetXaxis()->SetLabelSize(0.05);
 276   hSignal->GetXaxis()->SetTitleOffset(1.1);
 277   hSignal->GetXaxis()->SetNdivisions(510);
 278   hSignal->GetYaxis()->CenterTitle();
 279   hSignal->GetYaxis()->SetTitleSize(0.06); hSignal->GetYaxis()->SetLabelSize(0.05);
 280   hSignal->GetYaxis()->SetTitleOffset(1.1);
 281   hSignal->GetXaxis()->SetLabelFont(42); hSignal->GetYaxis()->SetLabelFont(42);
 282   hSignal->GetXaxis()->SetTitleFont(42); hSignal->GetYaxis()->SetTitleFont(42);
 283   hSignal->GetYaxis()->SetNdivisions(510);
 284   //hSignal->GetXaxis()->SetLimits(mMinPlot,mMmaxPlot);
 285   //hSignal->SetLineWidth(.7);
 286   hSignal->SetMarkerSize(.8);
 287   hSignal->SetMarkerColor(2); hSignal->SetMarkerStyle(20); hSignal->SetLineColor(2);
 288
 289   hSignal->Draw();
 290
 291   TLegend *leg2=new TLegend(0.16,0.8,0.57,0.93);
 292   leg2->SetBorderSize(0); leg2->SetFillColor(0); leg2->SetTextFont(42);
 293   leg2->SetFillStyle(0); leg2->SetMargin(0.25); //separation symbol-text
 294   leg2->SetEntrySeparation(0.15);
 295   leg2->AddEntry(hSignal, Form("OS-%4.2f*LS",sig->GetScaleFactor()), "p");
 296   TObject *hMmc=hSignal->GetListOfFunctions()->FindObject("mcMinv");
 297   if (hMmc) leg2->AddEntry(hMmc, Form("MC (#chi^{2}/dof=%3.1f)",TString(hMmc->GetTitle()).Atof()), "l");
 298   leg2->Draw();
 299
 300   Double_t beforePhys=0;
 301   Double_t afterPhys=0;
 302   Double_t afterV0and=0;
 303   Double_t afterEventSel=0;
 304   Double_t afterPileupRej=0;
 305
 306   if (hEventStat){
 307     Int_t beforePhysBin=hEventStat->GetXaxis()->FindBin("Before Phys. Sel.");
 308     Int_t afterPhysBin=hEventStat->GetXaxis()->FindBin("After Phys. Sel.");
 309     Int_t afterV0andBin=hEventStat->GetXaxis()->FindBin("V0and triggers");
 310     Int_t afterEventSelBin=hEventStat->GetXaxis()->FindBin("After Event Filter");
 311     Int_t afterPileupRejBin=hEventStat->GetXaxis()->FindBin("After Pileup rejection");
 312
 313     beforePhys=hEventStat->GetBinContent(beforePhysBin);
 314     afterPhys=hEventStat->GetBinContent(afterPhysBin);
 315     afterV0and=hEventStat->GetBinContent(afterV0andBin);
 316     afterEventSel=hEventStat->GetBinContent(afterEventSelBin);
 317     afterPileupRej=hEventStat->GetBinContent(afterPileupRejBin);
 318
 319     printf("Mevents: all: %4.1f, PhysSel: %4.1f, V0AND: %4.1f, EventSel: %4.1f, PileupRej: %4.1f\n",(Float_t)beforePhys/1.e+6,(Float_t)afterPhys/1.e+6,(Float_t)afterV0and/1.e+6,(Float_t)afterEventSel/1.e+6,(Float_t)afterPileupRej/1.e+6);
 320
 321     lat->DrawLatex(0.18, 0.2, Form("%4.1f Mevents",(Float_t)afterPhys/1.e+6));
 322   }
 323
 324   if (save){
 325 //     c->SaveAs(Form("%s%s.eps",cname,addToName.Data()));
 326     c->SaveAs(Form("%s%s.png",cname,addToName.Data()));
 327 /*
 328     FILE *out_file;
 329     if ( (out_file = fopen(Form("sig_%s.txt",cname), "w")) == NULL )
 330     {   fprintf(stderr, "Cannot open file %s\n", Form("sig_%s.txt",cname)); }
 331     fprintf(stdout, "Signal file: %s\n", Form("sig_%s.txt",cname));
 332     fprintf(out_file,"%3d %4.1f  %3.1f %4.2f  %4.1f %4.2f %d\n",(int)sigN,sigEr,sigS2B,sigS2Ber,sigSignif,sigSignifEr,(Int_t)afterPhys);
 333     fclose(out_file);
 334
 335     TFile outMinv(Form("Minv_%s.root",cname), "RECREATE");
 336     hUS->Write();
 337     hBackground->Write();
 338     hSignal->Write();
 339     if (hMmc) hMmc->Write();
 340     outMinv.Close();*/
 341
 342   }
 343
 344   return;
 345 }
 346
 347 //_______________________________________
 348 void FitMCshape(AliDielectronSignalBase *sig)
 349 {
 350   TFile mcFile(mcLineShapeFile);
 351   if (!mcFile.IsOpen()) {
 352     printf("mcMinv_LHC10e2 not found!!!\n");
 353     return;
 354   }
 355   TH1D *hMmc = (TH1D*)mcFile.Get("mcMinv");
 356   if (!hMmc) {
 357     printf("mcMinv not found!!\n");
 358     return;
 359   }
 360   hMmc->SetDirectory(0);
 361   hMmc->SetName("mcMinv");
 362   mcFile.Close();
 363
 364   TH1* hMsub=sig->GetSignalHistogram();
 365   Double_t mb1=sig->GetIntegralMin();
 366   Double_t mb2=sig->GetIntegralMax();
 367
 368   Double_t effInt = 0.;
 369   for(Int_t iBin=hMmc->FindBin(mb1); iBin<hMmc->FindBin(mb2); iBin++) {
 370     effInt += hMmc->GetBinContent(iBin);
 371   }
 372   effInt/=hMmc->Integral();
 373   printf("MC signal fraction in range %4.2f-%4.2f GeV: %5.3f \n",hMmc->GetBinLowEdge(hMmc->FindBin(mb1)),hMmc->GetBinLowEdge(hMmc->FindBin(mb2)+1),effInt);
 374
 375   Float_t mcScale1=(hMsub->GetXaxis()->GetBinWidth(1)/hMmc->GetXaxis()->GetBinWidth(1))*
 376     hMsub->Integral(hMsub->FindBin(mb1),hMsub->FindBin(mb2))/
 377     hMmc->Integral(hMmc->FindBin(mb1),hMmc->FindBin(mb2));
 378
 379   printf("1st guess of MC scale factor: %6.3f\n",mcScale1);
 380
 381   Float_t mcScale=0.;
 382   Float_t chi2_min=100000.;
 383   Int_t iMin=0;
 384   Int_t ndf=0;
 385
 386   for(Int_t i=0; i<20; i++){
 387     Float_t chi2=0.;
 388     Float_t scale=(0.4+0.05*(Float_t)i)*mcScale1;
 389     ndf=0;
 390     for(Int_t ib=1; ib<=hMsub->GetXaxis()->GetNbins(); ib++){
 391       Float_t data=(Float_t)hMsub->GetBinContent(ib);
 392       Float_t err=(Float_t)hMsub->GetBinError(ib);
 393       Float_t mc=scale*((Float_t)hMmc->GetBinContent(hMmc->FindBin(hMsub->GetBinCenter(ib))));
 394       if (err>0) {
 395         chi2 += ((data-mc)*(data-mc))/(err*err);
 396         ndf++;
 397       } else {
 398         //printf("bin %d Err: %6.3f, chi2: %6.1f\n",ib,err,chi2);
 399       }
 400     }
 401       //printf("%d scale factor: %6.3f, chi2: %6.1f\n",i,scale,chi2);
 402     if(chi2 < chi2_min){
 403       chi2_min = chi2;
 404       mcScale = scale;
 405       iMin=i;
 406     }
 407   }
 408   //Float_t chi2dof=chi2_min/(Float_t)(hMinv->GetXaxis()->GetNbins()-1);
 409   Float_t chi2dof=chi2_min/((Float_t)(ndf-1));
 410   printf("MC fit (i=%d): chi2/dof: %6.3f/%d, Scale: %7.4f \n",iMin,chi2_min,(ndf-1),mcScale);
 411   hMmc->SetTitle(Form("%f",chi2dof));
 412
 413   //mcScale=IntData/IntMC;printf("Int Data, MC: %10.1f %10.1f, MC scale: %6.3f\n",IntData,IntMC,mcScale);
 414
 415   hMmc->Scale(mcScale);
 416   hMmc->SetOption("sameHISTC");
 417   hMsub->GetListOfFunctions()->Add(hMmc);
 418 }