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 #include <TF1.h>
  14
  15 #include "AliDielectronSignalExt.h"
  16 #include "AliDielectronCFdraw.h"
  17 #include "AliDielectron.h"
  18 #include "AliDielectronHelper.h"
  19
  20 AliDielectronSignalBase* GetSignalLS(AliDielectronCFdraw &d, Int_t step,
  21                                      AliDielectronSignalBase::EBackgroundMethod type=AliDielectronSignalBase::kLikeSign);
  22 AliDielectronSignalBase* GetSignalLSY(AliDielectronCFdraw &d, Int_t step,
  23                                      AliDielectronSignalBase::EBackgroundMethod type=AliDielectronSignalBase::kLikeSign);
  24 AliDielectronSignalBase* GetSignalRot(AliDielectronCFdraw &d, Int_t step);
  25 void SetStyle(AliDielectronSignalBase *sig, const char* nameAdd);
  26 void DrawSpectra(AliDielectronSignalBase *sig, const char* cname, TH1  *hEventStat=0x0, Float_t fNevSel=1., Bool_t save=kTRUE);
  27 Float_t FitMCshape(AliDielectronSignalBase *sig);
  28
  29 const char *mcLineShapeFile="$ALICE_ROOT/PWGDQ/dielectron/macros/mcMinv_LHC10f7a.root";
  30
  31 //_______________________________________
  32 void PlotDataResults(const char* filenameData, const char* filenameMC="", Bool_t save=kTRUE, Bool_t DoPt=kTRUE )
  33 {
  34   AliDielectronCFdraw d(filenameData);
  35   AliDielectronCFdraw dCorr("corrCont","corrCont");
  36   TString nameCorr(filenameMC);
  37   if (!nameCorr.IsNull()) d.SetCFContainers(nameCorr.Data());
  38
  39   TFile f(filenameData);
  40   TH1 *hStats=(TH1*)f.Get("hEventStat");
  41   hStats->SetDirectory(0);
  42   TList* listData = (TList*)f.Get("jpsi_QA");
  43   //THashList *a = (THashList*)listData->FindObject("basicQ+SPDfirst+pt>.6+PID");
  44   THashList *a = (THashList*)listData->FindObject("basicQ+SPDfirst+pt>.8+PID"); //01.02.11
  45   //THashList *a = (THashList*)listData->FindObject("FullCuts+SPDany+Pi3.0+P3.0"); //Ionut
  46   THashList *b = (THashList*)a->FindObject("Event");
  47   TH1F *hZvtx= (TH1F*) b->FindObject("VtxZ");
  48   hZvtx->SetDirectory(0);
  49   f.Close();
  50
  51   hZvtx->GetXaxis()->SetTitleSize(.05); hZvtx->GetXaxis()->SetLabelSize(.038);
  52   hZvtx->GetXaxis()->SetTitleOffset(1.1);
  53   hZvtx->GetXaxis()->SetNdivisions(510);
  54   //hZvtx->GetYaxis()->CenterTitle();  hZvtx->GetXaxis()->CenterTitle();
  55   hZvtx->GetYaxis()->SetTitleSize(.05); hZvtx->GetYaxis()->SetLabelSize(.038);
  56   hZvtx->GetYaxis()->SetTitleOffset(1.3); //1.0
  57   hZvtx->GetXaxis()->SetLabelFont(42); hZvtx->GetYaxis()->SetLabelFont(42);
  58   hZvtx->GetXaxis()->SetTitleFont(42); hZvtx->GetYaxis()->SetTitleFont(42);
  59   hZvtx->GetYaxis()->SetNdivisions(510);
  60   hZvtx->SetLineWidth(2.5);
  61   hZvtx->SetMarkerSize(.8);
  62   hZvtx->SetMarkerColor(2); hZvtx->SetMarkerStyle(20); hZvtx->SetLineColor(2);
  63   hZvtx->GetXaxis()->SetTitle("Zvtx (cm)");
  64   hZvtx->GetYaxis()->SetTitle("Events");
  65
  66   TCanvas *c0=new TCanvas("c0","c0",10,0,450,450);
  67   c0->SetTopMargin(0.04);  c0->SetRightMargin(0.03);
  68   c0->SetLeftMargin(0.13); c0->SetBottomMargin(0.14);
  69
  70   TF1* gausFit=new TF1("gausFit", "gaus",-9.,9.);
  71   gausFit->SetLineColor(4); gausFit->SetLineWidth(.7);// gausFit->Draw("same");
  72   hZvtx->Fit(gausFit, "QM", "", -9.,9.);
  73   gStyle->SetOptFit(0);
  74   hZvtx->Draw();
  75
  76   //Float_t IntGaus=gausFit->Integral(gausFit->GetParameter(1)-2.0*gausFit->GetParameter(2),gausFit->GetParameter(1)+2.0*gausFit->GetParameter(2))/h->GetBinWidth(1);
  77   Float_t NevSel=hZvtx->Integral();
  78   Float_t IntGaus=gausFit->Integral(-25.,25.)/hZvtx->GetBinWidth(1);
  79   printf("NevSel: %6.1f, IntGauss: %6.1f \n",NevSel/1.e+6,IntGaus/1.e+6);
  80   Float_t fNevSel=NevSel/IntGaus;
  81
  82   TLatex *lat=new TLatex();
  83   lat->SetNDC(kTRUE);
  84   lat->SetTextColor(2);lat->SetTextFont(42);lat->SetTextSize(.035);
  85   lat->DrawLatex(0.16, 0.9, Form("NeventsSel: %5.1f M",NevSel/1.e+6));
  86
  87   TLatex *lat2=new TLatex();
  88   lat2->SetNDC(kTRUE);
  89   lat2->SetTextColor(4);lat2->SetTextFont(42);lat2->SetTextSize(.035);
  90   lat2->DrawLatex(0.7, 0.9, "Gauss fit:");
  91   lat2->DrawLatex(0.7, 0.86, Form("Nevents: %5.1f M",IntGaus/1.e+6));
  92   lat2->DrawLatex(0.7, 0.82, Form("Mean: %5.3f cm",gausFit->GetParameter(1)));
  93   lat2->DrawLatex(0.7, 0.78, Form("Sigma: %4.2f cm",gausFit->GetParameter(2)));
  94
  95   Int_t afterPhysBin=hStats->GetXaxis()->FindBin("After Phys. Sel.");
  96   Double_t afterPhys=hStats->GetBinContent(afterPhysBin);
  97   lat2->DrawLatex(0.4, 0.18, Form("NevPhysSel: %5.1f",(Float_t)afterPhys/1.e+6));
  98
  99   c0->SaveAs("Zvtx_cSel.eps");
 100
 101   //if (hZvtx) delete hZvtx;
 102
 103   Int_t stepFirst=0, stepAny=1, stepAnyV0=2, stepAnyV0T=3; //Data ...from 05.02.11 on
 104   //Int_t stepFirst=0, stepAny=1, stepTOFmix=2; //Data
 105   //Int_t stepFirst=2, stepAny=4, stepTOFmix=6; d.SetRangeUser("PairType",1,1); //MC only
 106   //Int_t stepFirst=1, stepAny=3, stepAnyV0T=5; //MC truth
 107   //Int_t stepFirst=2, stepAny=4, stepAnyV0T=5; //MC, all pairs ...or vice-versa
 108
 109   gStyle->SetOptStat(0);
 110   //Set common Ranges
 111   d.SetRangeUser("Leg1_NclsTPC",70.,170.); //was 90
 112   d.SetRangeUser("Leg2_NclsTPC",70.,170.);
 113   d.SetRangeUser("Leg1_Pt",1.01,100000);
 114   d.SetRangeUser("Leg2_Pt",1.01,100000);
 115   d.SetRangeUser("Leg1_Eta",-0.899,0.899);
 116   d.SetRangeUser("Leg2_Eta",-0.899,0.899);
 117
 118 //   d.SetRangeUser("Leg1_TPC_nSigma_Electrons",-3,3);
 119 //   d.SetRangeUser("Leg2_TPC_nSigma_Electrons",-3,3);
 120   d.SetRangeUser("Leg1_TPC_nSigma_Pions",3.51,20); // 2011-02-09_0027.5195 only
 121   d.SetRangeUser("Leg2_TPC_nSigma_Pions",3.51,20);
 122   d.SetRangeUser("Leg1_TPC_nSigma_Protons",3.01,20);
 123   d.SetRangeUser("Leg2_TPC_nSigma_Protons",3.01,20);
 124   d.SetRangeUser("M",0.5,5.);
 125   //============================
 126   //SPD first
 127   //
 128
 129   //--- Like sign subtraction
 130   AliDielectronSignalBase *sigFirst=GetSignalLS(d,stepFirst);
 131   SetStyle(sigFirst,"SPDfirst - Like Sign subtraction");
 132   DrawSpectra(sigFirst,"cFirst",hStats,fNevSel,save);
 133   //--- Like sign subtraction Arithmetic mean
 134   AliDielectronSignalBase *sigFirstArith=GetSignalLS(d,stepFirst,AliDielectronSignalBase::kLikeSignArithm);
 135   SetStyle(sigFirstArith,"SPDfirstArith - Like Sign subtraction");
 136   DrawSpectra(sigFirstArith,"cFirstArith",hStats,fNevSel,save);
 137   //--- Rotation subtraction
 138   AliDielectronSignalBase *sigFirstRot=GetSignalRot(d,stepFirst);
 139   SetStyle(sigFirstRot,"SPDfirst - Track rotation subtraction");
 140   DrawSpectra(sigFirstRot,"cFirstRot",hStats,save);
 141
 142   //============================
 143   //SPD any
 144   //
 145   AliDielectronSignalBase *sigAny=GetSignalLS(d,stepAny);
 146   SetStyle(sigAny,"ITS Any - Like Sign subtraction");
 147   DrawSpectra(sigAny,"cAny",hStats,fNevSel,save);
 148   //--- like sign with arithmetic mean
 149   AliDielectronSignalBase *sigAnyArith=GetSignalLS(d,stepAny,AliDielectronSignalBase::kLikeSignArithm);
 150   SetStyle(sigAnyArith,"SPDany - Like Sign subtraction (Arithm. mean)");
 151   DrawSpectra(sigAnyArith,"cAnyArith",hStats,fNevSel,save);
 152
 153   //--- Rotation subtraction
 154   AliDielectronSignalBase *sigAnyRot=GetSignalRot(d,stepAny);
 155   SetStyle(sigAnyRot,"SPDany - Track rotation subtraction");
 156   DrawSpectra(sigAnyRot,"cAnyRot",hStats,fNevSel,save);
 157
 158   //AliDielectronSignalBase *sigAnyRot2=GetSignalRot(d,stepAnyV0); //diff. rot. angle
 159   //SetStyle(sigAnyRot2,"ITS First - Track rotation subtraction");
 160   //DrawSpectra(sigAnyRot2,"cAnyRot2",hStats,fNevSel,save);
 161
 162   //AliDielectronSignalBase *sigAnyRot3=GetSignalRot(d,stepAnyV0T);  //diff. rot. angle
 163   //SetStyle(sigAnyRot3,"ITS First - Track rotation subtraction");
 164   //DrawSpectra(sigAnyRot3,"cAnyRot3",hStats,fNevSel,save);
 165
 166   //--- like sign with arithmetic mean, V0 (tender) conversions excl.
 167   //AliDielectronSignalBase *sigAnyV0tArith=GetSignalLS(d,stepAnyV0T,AliDielectronSignalBase::kLikeSignArithm);
 168   //SetStyle(sigAnyV0tArith,"ITS Any V0 (Tender) excl. - LS-Arithm.");
 169   //DrawSpectra(sigAnyV0tArith,"cAnyV0tArith",hStats,fNevSel,save);
 170
 171   //=============================
 172   //TOF up to 1.2, parametrisation in TPC ele
 173   //
 174 /*
 175   AliDielectronSignalBase *sigTOFmix=GetSignalLS(d,stepTOFmix);
 176   SetStyle(sigTOFmix,"TOF + TPC - Like Sign subtraction");
 177   DrawSpectra(sigTOFmix,"cTOFTPC",hStats,fNevSel,save);
 178   //--- Rotation subtraction
 179   AliDielectronSignalBase *sigTOFmixRot=GetSignalRot(d,stepTOFmix);
 180   SetStyle(sigTOFmixRot,"TOF + TPC - Track rotation subtraction");
 181   DrawSpectra(sigTOFmixRot,"cTOFTPCrot",hStats,fNevSel,save);
 182 */
 183
 184  //================================
 185  // y bins
 186   AliDielectronSignalBase *sigY_03_09=GetSignalLSY(d,stepAny,AliDielectronSignalBase::kLikeSignArithm);
 187   DrawSpectra(sigY_03_09,"cAny_Y03_09",hStats,fNevSel,save);
 188
 189   d.SetRangeUser("Y",-.299,+.299);
 190   AliDielectronSignalBase *sigY_03=GetSignalLS(d,stepAny,AliDielectronSignalBase::kLikeSignArithm);
 191   DrawSpectra(sigY_03,"cAny_Y03",hStats,fNevSel,save);
 192
 193 //  simulator.AliModule::SetDensityFactor(1.07);
 194
 195   if (DoPt){
 196
 197   //===============================
 198   // Pt bins
 199   //   "0.0, 0.8, 1.4, 2.8, 5., 9.9"
 200   d.SetRangeUser("Y",-.899,+.899);
 201   //Char_t *metPt="cAnyArith"; //method for background for Minv in pt bins
 202   Char_t *metPt="cAnyRot";
 203
 204
 205   FILE *out_file;
 206   if (save){
 207       if ( (out_file = fopen(Form("sigPt_%s.txt",metPt), "w")) == NULL )
 208       { fprintf(stderr, "Cannot open file %s\n", Form("sigPt_%s.txt",metPt)); }
 209   }
 210   const Int_t NptBins=6;
 211   //Float_t PtC[5]={0.5,1.25,2.25,4.,7.5}; //pt bins...need clever way
 212   //Float_t PtW[5]={0.5,0.25,0.75,1.0,2.5};
 213   //Float_t PtC[5]={0.5,1.5,2.5,4.,6.5}; //pt bins...23.01.11
 214   //Float_t PtW[5]={0.5,0.5,0.5,1.0,1.5};
 215   Float_t PtC[NptBins]={0.5,1.5,2.5,4.,6.,8.5}; //pt bins...01.02.11
 216   Float_t PtW[NptBins]={0.5,0.5,0.5,1.,1.,1.5};
 217
 218
 219   //TVectorD *vBins=AliDielectronHelper::MakeArbitraryBinning("0.0, 1, 1.5, 3, 5., 9.9");
 220   //TVectorD *vBins=AliDielectronHelper::MakeArbitraryBinning("0.0, 1, 2, 3, 5, 8"); //26.01
 221   TVectorD *vBins=AliDielectronHelper::MakeArbitraryBinning("0.0,1,2,3,5,7,10"); //01.02
 222
 223
 224   //AliDielectronSignalBase *sig_Pt[5];
 225   AliDielectronSignalBase *sig_Pt[NptBins];
 226
 227   //for (Int_t i=0; i<5; ++i){
 228   for (Int_t i=0; i<NptBins; ++i){
 229     d.SetRangeUser("Pt",(*vBins)[i]+.001, (*vBins)[i+1]-0.001);
 230     if ("cAnyArith" == metPt){
 231         sig_Pt[i]=GetSignalLS(d,stepAny,AliDielectronSignalBase::kLikeSignArithm);
 232     }
 233     else if ("cAnyRot" == metPt){
 234         sig_Pt[i]=GetSignalRot(d,stepAny); //rotation
 235     }
 236     else {
 237         cout << " ??? No such method ??? " <<metPt <<endl;
 238     }
 239     DrawSpectra(sig_Pt[i],Form("%s_Pt%d",metPt,i),hStats,fNevSel,save);
 240   }
 241
 242   TH1D *hSigPt=new TH1D("hSigPt","hSigPt",NptBins,vBins->GetMatrixArray());
 243   TH1D *hSigSign=new TH1D("hSigSign","hSigSign",NptBins,vBins->GetMatrixArray());
 244   TH1D *hSigSOB=new TH1D("hSigSOB","hSigSOB",NptBins,vBins->GetMatrixArray());
 245
 246   for (Int_t i=0; i<NptBins; ++i){
 247       hSigPt->SetBinContent(i+1,sig_Pt[i]->GetSignal());
 248       hSigPt->SetBinError(i+1,sig_Pt[i]->GetSignalError());
 249       hSigSign->SetBinContent(i+1,sig_Pt[i]->GetSignificance());
 250       hSigSign->SetBinError(i+1,sig_Pt[i]->GetSignificanceError());
 251       hSigSOB->SetBinContent(i+1,sig_Pt[i]->GetSB());
 252       hSigSOB->SetBinError(i+1,sig_Pt[i]->GetSBError());
 253       if (save){
 254           fprintf(out_file,"%4.1f %5.2f %5.1f  %4.1f  %3.1f %4.2f  %4.1f %4.2f \n",PtC[i],PtW[i],hSigPt->GetBinContent(i+1),hSigPt->GetBinError(i+1),hSigSign->GetBinContent(i+1),hSigSign->GetBinError(i+1),hSigSOB->GetBinContent(i+1),hSigSOB->GetBinError(i+1));
 255       }
 256
 257   }
 258
 259   if (save){
 260       fclose(out_file);
 261       fprintf(stdout, " *** Signal file vs. pt: %s\n", Form("sigPt_%s.txt",metPt));
 262   }
 263
 264   const Int_t kMarkTyp=20; //marker type
 265   const Int_t kMarkCol=2; //...and color
 266   const Float_t kTitSize=0.08; //axis title size
 267   const Float_t kAsize=0.85*kTitSize; //...and label size
 268   const Float_t kToffset=0.8;
 269   const Int_t kFont=42;
 270
 271   TCanvas *cSigPt=(TCanvas*)gROOT->FindObject("cSigPt");
 272
 273   if (!cSigPt) cSigPt=new TCanvas("cSigPt","cSigPt",5,5,450,700);
 274   cSigPt->Clear();
 275   //cSigPt->Divide(2,2);
 276   cSigPt->SetTopMargin(0.04);  cSigPt->SetRightMargin(0.04);
 277   cSigPt->SetLeftMargin(0.13);  cSigPt->SetBottomMargin(0.2);
 278   cSigPt->Divide(1,3,0,0);
 279   cSigPt->cd(1);
 280
 281   hSigPt->SetMarkerStyle(kMarkTyp);
 282   hSigPt->SetMarkerColor(kMarkCol); hSigPt->SetLineColor(kMarkCol);
 283   hSigPt->GetYaxis()->SetTitleOffset(kToffset);
 284   hSigPt->SetTitleSize(kTitSize,"XY");
 285   hSigPt->SetTitleFont(kFont,"XY");
 286   hSigPt->SetLabelSize(kAsize,"YX");
 287   hSigPt->SetLabelFont(kFont,"XY");
 288   //gS->SetAxisRange(0.,Pt2[Npt-1],"X");
 289   //hSigPt->GSetAxisRange(Y1m,Y1M,"Y");
 290   hSigPt->GetYaxis()->SetTitle("N_{J/#psi}");
 291   hSigPt->Draw();
 292
 293
 294   cSigPt->cd(2);
 295   hSigSign->SetMarkerStyle(kMarkTyp);
 296   hSigSign->SetMarkerColor(kMarkCol); hSigSign->SetLineColor(kMarkCol);
 297   hSigSign->GetYaxis()->SetTitleOffset(kToffset);
 298   hSigSign->SetTitleSize(kTitSize,"XY");
 299   hSigSign->SetTitleFont(kFont,"XY");
 300   hSigSign->SetLabelSize(kAsize,"YX");
 301   hSigSign->SetLabelFont(kFont,"XY");
 302   hSigSign->GetYaxis()->SetTitle("Significance");
 303   hSigSign->Draw();
 304
 305   cSigPt->cd(3);
 306   hSigSOB->SetMaximum(4.);
 307   hSigSOB->SetMinimum(0.);
 308   hSigSOB->SetMarkerStyle(kMarkTyp);
 309   hSigSOB->SetMarkerColor(kMarkCol); hSigSOB->SetLineColor(kMarkCol);
 310   hSigSOB->GetYaxis()->SetTitleOffset(kToffset);
 311   hSigSOB->SetTitleSize(kTitSize,"XY");
 312   hSigSOB->SetTitleFont(kFont,"XY");
 313   hSigSOB->SetLabelSize(kAsize,"YX");
 314   hSigSOB->SetLabelFont(kFont,"XY");
 315   hSigSOB->GetXaxis()->SetTitle("p_{t}^{J/#psi} (GeV/c");
 316   hSigSOB->GetYaxis()->SetTitle("S/B");
 317   hSigSOB->Draw();
 318
 319   cSigPt->SaveAs(Form("SigPt_%s.eps",metPt));
 320
 321   } //DoPt
 322
 323   if (hStats) delete hStats;
 324
 325
 326 }
 327
 328
 329 //_______________________________________
 330 AliDielectronSignalBase *GetSignalLS(AliDielectronCFdraw &d, Int_t step, AliDielectronSignalBase::EBackgroundMethod type)
 331 {
 332   //
 333   // Get Extracted signal from likesign method
 334   //
 335
 336   TObjArray *arr=new TObjArray;
 337   arr->SetOwner();
 338
 339   for (Int_t iType=0;iType<3;++iType){
 340     d.SetRangeUser("PairType",iType,iType);
 341     TH1D* hh = (TH1D*)d.Project("M",step);
 342     if(TMath::Abs(hh->GetBinWidth(1)-0.02)<0.0001) {
 343         //hh->Rebin(2);
 344     }
 345     arr->AddAt(hh,iType);
 346   }
 347
 348   AliDielectronSignalExt *sig=new AliDielectronSignalExt;
 349   sig->SetScaleRawToBackground(3.2,4.9);
 350   //sig->SetScaleRawToBackground(3.2,5.0);
 351   sig->SetIntegralRange(2.92,3.15);
 352   sig->SetMethod(type);
 353   sig->Process(arr);
 354
 355   delete arr;
 356   return sig;
 357 }
 358
 359
 360 //_______________________________________
 361 AliDielectronSignalBase *GetSignalLSY(AliDielectronCFdraw &d, Int_t step, AliDielectronSignalBase::EBackgroundMethod type)
 362 {
 363   //
 364   // Get Extracted signal from likesign method
 365   //
 366
 367   TObjArray *arr=new TObjArray;
 368   arr->SetOwner();
 369
 370   d.SetRangeUser("Y",-.899,-.301);
 371
 372   for (Int_t iType=0;iType<3;++iType){
 373     d.SetRangeUser("PairType",iType,iType);
 374     arr->AddAt(d.Project("M",step),iType);
 375   }
 376
 377   d.SetRangeUser("Y",+.301,+.899);
 378
 379   for (Int_t iType=0;iType<3;++iType){
 380     d.SetRangeUser("PairType",iType,iType);
 381     TH1 *h=d.Project("M",step);
 382
 383     TH1 *h2=(TH1*)arr->At(iType);
 384     h2->Add(h);
 385     delete h;
 386   }
 387
 388   AliDielectronSignalExt *sig=new AliDielectronSignalExt;
 389   sig->SetScaleRawToBackground(3.2,4.9);
 390   sig->SetIntegralRange(2.93,3.15);
 391   sig->SetMethod(type);
 392   sig->Process(arr);
 393
 394   delete arr;
 395   return sig;
 396 }
 397
 398
 399 //_______________________________________
 400 AliDielectronSignalBase *GetSignalRot(AliDielectronCFdraw &d, Int_t step)
 401 {
 402   //
 403   // Get Extracted signal from likesign method
 404   //
 405
 406   TObjArray *arr=new TObjArray;
 407   arr->SetOwner();
 408
 409   Int_t iType=AliDielectron::kEv1PM;
 410   d.SetRangeUser("PairType",iType,iType);
 411   arr->AddAt(d.Project("M",step),iType);
 412
 413   iType=AliDielectron::kEv1PMRot;
 414   d.SetRangeUser("PairType",iType,iType);
 415   arr->AddAt(d.Project("M",step),iType);
 416
 417   AliDielectronSignalExt *sig=new AliDielectronSignalExt;
 418   sig->SetScaleRawToBackground(3.2,4.9);
 419   sig->SetIntegralRange(2.93,3.15);
 420   sig->SetMethod(AliDielectronSignalBase::kRotation);
 421   sig->Process(arr);
 422
 423
 424   delete arr;
 425   return sig;
 426 }
 427
 428 //_______________________________________
 429 void SetStyle(AliDielectronSignalBase *sig, const char* nameAdd)
 430 {
 431   //
 432   //
 433   //
 434   TH1 *hUS=sig->GetUnlikeSignHistogram();
 435   hUS->SetMarkerStyle(20);
 436   hUS->SetMarkerSize(0.7);
 437   hUS->SetMarkerColor(kRed);
 438   hUS->SetLineColor(kRed);
 439   hUS->SetStats(0);
 440   hUS->SetTitle(Form("Like sign spectrum %s;M inv. ee;counts per %.4g MeV",
 441                      nameAdd,1000*hUS->GetXaxis()->GetBinWidth(1)));
 442
 443   TH1* hBackground=sig->GetBackgroundHistogram();
 444   hBackground->SetMarkerStyle(24);
 445   hBackground->SetMarkerSize(0.7);
 446   hBackground->SetStats(0);
 447   hBackground->SetMarkerColor(kBlue);
 448   hBackground->SetLineColor(kBlue);
 449   hBackground->SetTitle(Form("Like sign spectrum %s;M inv. ee;counts per %.4g MeV",
 450                              nameAdd,1000*hBackground->GetXaxis()->GetBinWidth(1)));
 451
 452   TH1* hSignal=sig->GetSignalHistogram();
 453   hSignal->SetMarkerStyle(20);
 454   hSignal->SetMarkerSize(0.7);
 455   hSignal->SetMarkerColor(kRed);
 456   hSignal->SetLineColor(kRed);
 457   hSignal->SetTitle(Form("Like sign spectrum %s;M inv. ee;counts per %.4g MeV",
 458                          nameAdd,1000*hSignal->GetXaxis()->GetBinWidth(1)));
 459
 460
 461 }
 462
 463 //_______________________________________
 464 void DrawSpectra(AliDielectronSignalBase *sig, const char* cname, TH1  *hEventStat, Float_t fNevSel, Bool_t save)
 465 {
 466   //
 467   //
 468   //
 469   Float_t effInt=FitMCshape(sig);
 470
 471   gStyle->SetOptTitle(0);
 472   TCanvas *c=(TCanvas*)gROOT->FindObject(cname);
 473   if (!c) c=new TCanvas(cname,cname,400,500);
 474   c->SetTopMargin(0.04);  c->SetRightMargin(0.04);
 475   c->SetLeftMargin(0.13); c->SetBottomMargin(0.14);
 476   c->Clear();
 477   c->Divide(1,2,0,0);
 478   c->cd(1);
 479   gPad->SetTopMargin(0.01);
 480   gPad->SetRightMargin(0.01);
 481   gPad->SetLeftMargin(0.15);
 482   gPad->SetBottomMargin(0.001);
 483   /*
 484   gPad->SetGridx();
 485   gPad->SetGridy();
 486 */
 487   gPad->SetTickx();
 488   gPad->SetTicky();
 489
 490   TH1 *hUS=sig->GetUnlikeSignHistogram();
 491
 492   TH1* hBackground=sig->GetBackgroundHistogram();
 493
 494   hUS->GetXaxis()->CenterTitle();
 495   hUS->GetXaxis()->SetTitleSize(0.07); hUS->GetXaxis()->SetLabelSize(0.059);
 496 //   hUS->GetXaxis()->SetLabelOffset(1.);
 497   hUS->GetXaxis()->SetTitleOffset(1.1);
 498   hUS->GetXaxis()->SetNdivisions(510);
 499   hUS->GetYaxis()->CenterTitle();
 500   hUS->GetYaxis()->SetTitleSize(0.07); hUS->GetYaxis()->SetLabelSize(0.059);
 501   hUS->GetYaxis()->SetTitleOffset(1.);
 502   hUS->GetXaxis()->SetLabelFont(42); hUS->GetYaxis()->SetLabelFont(42);
 503   hUS->GetXaxis()->SetTitleFont(42); hUS->GetYaxis()->SetTitleFont(42);
 504   hUS->GetYaxis()->SetNdivisions(510);
 505   //hUS->GetXaxis()->SetLimits(mMinPlot,mMmaxPlot);
 506   //hUS->SetLineWidth(.7);
 507   hUS->SetMarkerSize(.8);
 508   hUS->SetMarkerColor(2); hUS->SetMarkerStyle(20); hUS->SetLineColor(2);
 509   hUS->Draw();
 510   hUS->SetMaximum(hUS->GetMaximum()*1.3);
 511
 512   //hBackground->SetLineWidth(.4);
 513   hBackground->SetMarkerSize(.7);
 514   hBackground->SetMarkerColor(4); hBackground->SetMarkerStyle(27); hBackground->SetLineColor(4);
 515
 516   hBackground->Draw("same");
 517
 518   TLatex *lat=new TLatex;
 519   lat->SetNDC(kTRUE);
 520   //  lat->DrawLatex(0.68, 0.67, "ALICE Performance");
 521   lat->SetTextColor(1);lat->SetTextFont(42);lat->SetTextSize(.045);
 522
 523   Double_t sigN=sig->GetSignal();
 524   Double_t sigEr=sig->GetSignalError();
 525   Double_t sigS2B=sig->GetSB();
 526   Double_t sigS2Ber=sig->GetSBError();
 527   Double_t sigSignif= sig->GetSignificance();
 528   Double_t sigSignifEr= sig->GetSignificanceError();
 529   lat->DrawLatex(0.18, 0.92, Form("S: %5.1f#pm%4.1f, S/B: %3.1f#pm %4.2f, Signif.: %4.1f#pm%4.2f (%4.2f-%4.2f GeV) ",sigN,sigEr,sigS2B,sigS2Ber,sigSignif,sigSignifEr,
 530                        hUS->GetBinLowEdge(hUS->FindBin(sig->GetIntegralMin())),
 531                        hUS->GetBinLowEdge(hUS->FindBin(sig->GetIntegralMax())+1)));
 532
 533   TLegend *leg=new TLegend(0.17,0.72,0.42,0.88);
 534   leg->SetBorderSize(0); leg->SetFillColor(0); leg->SetTextFont(42);
 535   leg->SetFillStyle(0); leg->SetMargin(0.25); //separation symbol-text
 536   leg->SetEntrySeparation(0.15);
 537   leg->AddEntry(hUS,"OS", "p");
 538   leg->AddEntry(hBackground, Form("LS*%4.2f",sig->GetScaleFactor()), "p");
 539   leg->Draw();
 540
 541   c->cd(2);
 542   gPad->SetRightMargin(0.01);
 543   gPad->SetTopMargin(0);
 544   gPad->SetLeftMargin(0.15);
 545   gPad->SetBottomMargin(.15);
 546   /*
 547   gPad->SetGridx();
 548   gPad->SetGridy();
 549   */
 550   gPad->SetTickx();
 551   gPad->SetTicky();
 552
 553   TH1* hSignal=sig->GetSignalHistogram();
 554   hSignal->GetXaxis()->SetTitle("M_{ee} (GeV/c^{2})");
 555   hSignal->GetXaxis()->CenterTitle();
 556   hSignal->GetXaxis()->SetTitleSize(0.06); hSignal->GetXaxis()->SetLabelSize(0.05);
 557   hSignal->GetXaxis()->SetTitleOffset(1.1);
 558   hSignal->GetXaxis()->SetNdivisions(510);
 559   hSignal->GetYaxis()->CenterTitle();
 560   hSignal->GetYaxis()->SetTitleSize(0.06); hSignal->GetYaxis()->SetLabelSize(0.05);
 561   hSignal->GetYaxis()->SetTitleOffset(1.1);
 562   hSignal->GetXaxis()->SetLabelFont(42); hSignal->GetYaxis()->SetLabelFont(42);
 563   hSignal->GetXaxis()->SetTitleFont(42); hSignal->GetYaxis()->SetTitleFont(42);
 564   hSignal->GetYaxis()->SetNdivisions(510);
 565   //hSignal->GetXaxis()->SetLimits(mMinPlot,mMmaxPlot);
 566   //hSignal->SetLineWidth(.7);
 567   hSignal->SetMarkerSize(.8);
 568   hSignal->SetMarkerColor(2); hSignal->SetMarkerStyle(20); hSignal->SetLineColor(2);
 569
 570   hSignal->Draw();
 571
 572   TLegend *leg2=new TLegend(0.16,0.8,0.57,0.93);
 573   leg2->SetBorderSize(0); leg2->SetFillColor(0); leg2->SetTextFont(42);
 574   leg2->SetFillStyle(0); leg2->SetMargin(0.25); //separation symbol-text
 575   leg2->SetEntrySeparation(0.15);
 576   leg2->AddEntry(hSignal, Form("OS-%4.2f*LS",sig->GetScaleFactor()), "p");
 577   TObject *hMmc=hSignal->GetListOfFunctions()->FindObject("mcMinv");
 578   if (hMmc) leg2->AddEntry(hMmc, Form("MC (#chi^{2}/dof=%3.1f)",TString(hMmc->GetTitle()).Atof()), "l");
 579   leg2->Draw();
 580
 581   Double_t beforePhys=0;
 582   Double_t afterPhys=0;
 583   Double_t afterV0and=0;
 584   Double_t afterEventSel=0;
 585   Double_t afterPileupRej=0;
 586
 587   if (hEventStat){
 588     Int_t beforePhysBin=hEventStat->GetXaxis()->FindBin("Before Phys. Sel.");
 589     Int_t afterPhysBin=hEventStat->GetXaxis()->FindBin("After Phys. Sel.");
 590     Int_t afterV0andBin=hEventStat->GetXaxis()->FindBin("V0and triggers");
 591     Int_t afterEventSelBin=hEventStat->GetXaxis()->FindBin("After Event Filter");
 592     Int_t afterPileupRejBin=hEventStat->GetXaxis()->FindBin("After Pileup rejection");
 593
 594     beforePhys=hEventStat->GetBinContent(beforePhysBin);
 595     afterPhys=hEventStat->GetBinContent(afterPhysBin);
 596     afterV0and=hEventStat->GetBinContent(afterV0andBin);
 597     afterEventSel=hEventStat->GetBinContent(afterEventSelBin);
 598     afterPileupRej=hEventStat->GetBinContent(afterPileupRejBin);
 599
 600     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);
 601
 602     afterPhys*=fNevSel; //fraction after Zvtx cut
 603
 604     lat->DrawLatex(0.18, 0.2, Form("%4.1f Mevents",(Float_t)afterPhys/1.e+6));
 605   }
 606
 607   if (save){
 608
 609     FILE *out_file;
 610     if ( (out_file = fopen(Form("sig_%s.txt",cname), "w")) == NULL )
 611     {   fprintf(stderr, "Cannot open file %s\n", Form("sig_%s.txt",cname)); }
 612     fprintf(stdout, "Signal file: %s\n", Form("sig_%s.txt",cname));
 613     fprintf(out_file,"%3d %4.1f  %3.1f %4.2f  %4.1f %4.2f  %d  %4.2f  %6.3f\n",(int)sigN,sigEr,sigS2B,sigS2Ber,sigSignif,sigSignifEr,(Int_t)afterPhys,sig->GetScaleFactor(),effInt);
 614     fclose(out_file);
 615
 616     TFile outMinv(Form("Minv_%s.root",cname), "RECREATE");
 617     hUS->Write();
 618     hBackground->Write();
 619     hSignal->Write();
 620     if (hMmc) hMmc->Write();
 621     outMinv.Close();
 622
 623     c->SaveAs(Form("Minv_%s.eps",cname));
 624     //c->SaveAs(Form("%s.png",cname));
 625
 626   }
 627
 628   return;
 629 }
 630
 631 //_______________________________________
 632 Float_t FitMCshape(AliDielectronSignalBase *sig)
 633 {
 634   TFile mcFile(mcLineShapeFile);
 635   if (!mcFile.IsOpen()) {
 636     printf("MC Minv file not found!!!\n");
 637     return 0;
 638   }
 639   //TH1D *hMmc = (TH1D*)mcFile.Get("mcMinv");
 640   TH1D *hMmc = (TH1D*)mcFile.Get("HistSignal");
 641   if (!hMmc) {
 642     printf("mcMinv not found!!\n");
 643     return 0;
 644   }
 645   hMmc->SetDirectory(0);
 646   hMmc->SetName("mcMinv");
 647   mcFile.Close();
 648
 649   TH1* hMsub=sig->GetSignalHistogram();
 650   Double_t mb1=sig->GetIntegralMin();
 651   Double_t mb2=sig->GetIntegralMax();
 652
 653   Double_t effInt = 0.;
 654   for(Int_t iBin=hMmc->FindBin(mb1); iBin<=hMmc->FindBin(mb2); iBin++) {
 655     effInt += hMmc->GetBinContent(iBin);
 656   }
 657   effInt/=hMmc->Integral();
 658   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);
 659
 660   Float_t MminFit=1.2; //min., Minv for fit (and chi2 calc.)
 661   Float_t mcScale1=(hMsub->GetXaxis()->GetBinWidth(1)/hMmc->GetXaxis()->GetBinWidth(1))*
 662     hMsub->Integral(hMsub->FindBin(mb1),hMsub->FindBin(mb2))/
 663     hMmc->Integral(hMmc->FindBin(mb1),hMmc->FindBin(mb2));
 664
 665   printf("1st guess of MC scale factor: %6.3f ...chi2 for Minv>%4.1f\n",mcScale1,MminFit);
 666
 667   Float_t mcScale=0.;
 668   Float_t chi2_min=100000.;
 669   Int_t iMin=0;
 670   Int_t ndf=0;
 671
 672   for(Int_t i=0; i<20; i++){
 673     Float_t chi2=0.;
 674     Float_t scale=(0.4+0.05*(Float_t)i)*mcScale1;
 675     ndf=0;
 676     for(Int_t ib=1; ib<=hMsub->GetXaxis()->GetNbins(); ib++){
 677         if (hMsub->GetBinCenter(ib) > MminFit) {
 678             Float_t data=(Float_t)hMsub->GetBinContent(ib);
 679             Float_t err=(Float_t)hMsub->GetBinError(ib);
 680             Float_t mc=scale*((Float_t)hMmc->GetBinContent(hMmc->FindBin(hMsub->GetBinCenter(ib))));
 681             if (err>0) {
 682                 chi2 += ((data-mc)*(data-mc))/(err*err);
 683                 ndf++;
 684             } else {
 685                 printf("bin %d Err: %6.3f, chi2: %6.1f\n",ib,err,chi2);
 686             }
 687         }
 688     }
 689       //printf("%d scale factor: %6.3f, chi2: %6.1f\n",i,scale,chi2);
 690     if(chi2 < chi2_min){
 691       chi2_min = chi2;
 692       mcScale = scale;
 693       iMin=i;
 694     }
 695   }
 696   //Float_t chi2dof=chi2_min/(Float_t)(hMinv->GetXaxis()->GetNbins()-1);
 697   Float_t chi2dof=chi2_min/((Float_t)(ndf-1));
 698   printf("MC fit (i=%d): chi2/dof: %6.3f/%d, Scale: %7.4f \n",iMin,chi2_min,(ndf-1),mcScale);
 699   hMmc->SetTitle(Form("%f",chi2dof));
 700
 701   //mcScale=IntData/IntMC;printf("Int Data, MC: %10.1f %10.1f, MC scale: %6.3f\n",IntData,IntMC,mcScale);
 702
 703   hMmc->Scale(mcScale);
 704   hMmc->SetOption("sameHISTC");
 705   hMmc->SetLineColor(1); hMmc->SetMarkerColor(1);
 706   hMsub->GetListOfFunctions()->Add(hMmc);
 707
 708   return effInt;
 709
 710 }