PWG0/multVScentPbPb/CorrectSpectraMulti.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include "TList.h"
   3 #include "TFile.h"
   4 #include "TStyle.h"
   5 #include "TH1F.h"
   6 #include "TH2F.h"
   7 #include "THnSparse.h"
   8 #include "TLegend.h"
   9 #include "TSystem.h"
  10 #include "TMath.h"
  11 #include "TCanvas.h"
  12 #include "TLegend.h"
  13 #include "TLatex.h"
  14 #include "TF1.h"
  15 #include "TLine.h"
  16 #include "TPaletteAxis.h"
  17 #include "TArrayD.h"
  18 #include "TGraphErrors.h"
  19 //
  20 //
  21 #endif
  22 #include "/home/shahoian/ALICE/mymacros/SaveCanvas.h"
  23
  24
  25
  26 const char kHStatName[]="hStat";
  27 double kEps = 1e-6;
  28 //
  29 double kdPhiBgTailMin = 0.1; // lower limit of dphi tail to use for bg normalization
  30 double kdPhiBgTailMax = 0.3; // upper limit of dphi tail to use for bg normalization
  31 //
  32 double kWDistBgTailMin = 5.; // lower limit of wgh.distance tail to use for bg normalization
  33 double kWDistBgTailMax = 25.; // upper limit of wgh.distance tail to use for bg normalization
  34
  35 double kdPhiSgCut=-1;       // cut in dphi-bent used to extract the signal, extracted from stat histo
  36 double kWDistSgCut=-1;      // cut in w.distance used to extract the signal, extracted from stat histo
  37 //
  38 enum { kNormShapeDist,      // normalize bg tails usig weighted distance shape
  39        kNormShapeDPhi,      // normalize bg tails usig dPhi-bend shape
  40        kNormShapes};
  41
  42 enum { kSclWghMean,         // normalize bg tails to data using weighted mean of bin-by-bin ratios
  43        kSclIntegral,        // normalize bg tails to data using integral
  44        kSclTypes};
  45
  46
  47 const char* figDir = "figMult";
  48 TString  useBgType    = "Inj";
  49 Int_t    useShapeType = kNormShapeDist;    // which distribution to use for bg normalization
  50 Bool_t   useMCLB      = 0;//kFALSE;             // use Comb MC Labels as a template for Bg.
  51 Int_t    useScaleType = kSclIntegral;//kSclWghMean;       // which type of tails normalization to use
  52 const double kEtaFitRange = 0.5;
  53
  54 enum {kBitNormPerEvent=BIT(14)};
  55   // bins for saved parameters in the hStat histo
  56   enum {kDummyBin,
  57         kEvTot0,      // events read
  58         kEvTot,       // events read after vertex quality selection
  59         kOneUnit,     // just 1 to track primate merges
  60         kNWorkers,    // n workers
  61         //
  62         kCentVar,     // cetrality var. used
  63         kDPhi,        // dphi window
  64         kDTht,        // dtheta window
  65         kNStd,        // N.standard deviations to keep
  66         kPhiShift,    // bending shift
  67         kThtS2,       // is dtheta scaled by 1/sin^2
  68         kThtCW,       // on top of w.dist cut cut also on 1 sigma dThetaX
  69         kPhiOvl,      // overlap params
  70         kZEtaOvl,     // overlap params
  71         kNoOvl,       // flag that overlap are suppressed
  72         //
  73         kPhiRot,      // rotation phi
  74         kInjScl,      // injection scaling
  75         kEtaMin,      // eta cut
  76         kEtaMax,      // eta cut
  77         kZVMin,       // min ZVertex to process
  78         kZVMax,       // max ZVertex to process
  79         //
  80         kDPiSCut,     // cut on dphi used to extract signal (when WDist is used in analysis, put it equal to kDPhi
  81         kNStdCut,     // cut on weighted distance (~1) used to extract signal
  82         //
  83         kMCV0Scale,   // scaling value for V0 in MC
  84         //
  85         // here we put entries for each mult.bin
  86         kBinEntries = 50,
  87         kEvProcData,  // events with data mult.object (ESD or reco)
  88         kEvProcInj,   // events Injected, total
  89         kEvProcRot,   // events Rotated
  90         kEvProcMix,   // events Mixed
  91         kEntriesPerBin
  92   };
  93
  94
  95 enum {kSigCorr,kMCPrim,kRawDtCut,kSignalEst,kSignalEstMC,kBgEst,k1MBeta,k1MBetaMC,kAlpha,kAlphaMC,kBgMC,kBgRescFc,kDataDist,kBgDist,kBgMCDist, kMCShift=20, kNHistos=kMCShift+kMCShift};
  96
  97 void    CorrectSpectraMulti(const char* flNameData, const char* flNameMC,const char* unique="");
  98 Bool_t  PrepareHistos(int bin, TList* lst, Bool_t isMC);
  99 void    ProcessHistos(int bin);
 100 TH1*    NormalizeBg(TH1* dataH, TH1* bgH, double &scl, double &scle);
 101 TObject* FindObject(int bin, const char* nameH, const TList* lst, Bool_t normPerEvent=kTRUE);
 102 TList*  LoadList(const char* flName, const char* addPref, const char* nameL="clist");
 103 void    GetRatE(double x,double xe, double y,double ye, double &rat, double &rate);
 104 Int_t   CheckStat(const TList* lst,const char* dtType);
 105 void    Integrate(TH1* hist, double xmn,double xmx, double &val, double& err);
 106 void    CropHisto(TH1* histo, int b00, int b01, int b10=-1, int b11=-1);
 107 void    CropHisto(TH1* histo, double b00, double b01, double b10=-1, double b11=-1);
 108 void    GetRealMinMax(TH1* h, double &vmn, double &vmx);
 109 const char*   HName(const char* prefix,const char* htype);
 110
 111 void PlotResults();
 112 void PlotDNDEta(int bin);
 113 void PlotAlphaBeta(int bin);
 114 void PlotSpecies();
 115
 116 Float_t myMergeFactor = -1; // if the files were manually merged, scale everything except statistics by 1/myMergeFactor
 117 Int_t nCentBins = -1;
 118 TList *listDt=0, *listMC=0;
 119 TObjArray resArr, resDnDeta;
 120 char outStr[1000];
 121 char outTitle[1000];
 122 TString uniqueName="";
 123 //
 124 TArrayD dNdEta,dNdEtaErr;
 125 TCanvas *canvFin=0;
 126 //
 127 Bool_t creatDnDEtaCMacro = kFALSE;
 128 Bool_t creatAlphaBetaCMacro = kFALSE;
 129 Bool_t creatSpeciesCMacro = kFALSE;
 130
 131 void CorrectSpectraMulti(const char* flNameData, const char* flNameMC, const char* uniqueNm)
 132 {
 133   //
 134   uniqueName = uniqueNm;
 135   listDt = LoadList(flNameData,"dt_");
 136   listMC = LoadList(flNameMC,"mc_");
 137   //
 138   resArr.Clear();
 139   //
 140   TH1* hstat = (TH1*)FindObject(-1,kHStatName,listDt,kFALSE);
 141   //
 142   int nbstat = hstat->GetNbinsX();
 143   nCentBins = (nbstat - kBinEntries)/kEntriesPerBin;
 144   printf("%d bins will be processed\n",nCentBins);
 145   if (nCentBins<1) return;
 146   myMergeFactor = hstat->GetBinContent(kOneUnit);
 147   printf("Detected %f mergings\n",myMergeFactor);
 148   //
 149   dNdEta.Set(nCentBins);
 150   dNdEtaErr.Set(nCentBins);
 151   //
 152   kdPhiSgCut  = hstat->GetBinContent(kDPiSCut)/myMergeFactor;
 153   kWDistSgCut = hstat->GetBinContent(kNStdCut)/myMergeFactor;
 154   printf("Signal cuts used: dPhiS: %f WDist:%f\n",kdPhiSgCut,kWDistSgCut);
 155   //
 156   for (int ib=0;ib<nCentBins;ib++) {
 157     if (!PrepareHistos(ib,listMC,kTRUE))  return;
 158     if (!PrepareHistos(ib,listDt,kFALSE)) return;
 159     ProcessHistos(ib);
 160     //
 161   }
 162   //
 163   sprintf(outStr,"CutEta%.1f_%.1f_Zv%.1f_%.1f_bg%s_Shape_%s_mcLB%d_cutSig%.1f_cutBg%.1f",
 164           hstat->GetBinContent(kEtaMin)/myMergeFactor,
 165           hstat->GetBinContent(kEtaMax)/myMergeFactor,
 166           hstat->GetBinContent(kZVMin)/myMergeFactor,
 167           hstat->GetBinContent(kZVMax)/myMergeFactor,
 168           useBgType.Data(),
 169           useShapeType==kNormShapeDist ? "wdst":"dphi",
 170           useMCLB,
 171           useShapeType==kNormShapeDist ? kWDistSgCut:kdPhiSgCut,
 172           useShapeType==kNormShapeDist ? kWDistBgTailMin:kdPhiBgTailMin);
 173   //
 174   PlotResults();
 175   //
 176   printf("Final Results:\n");
 177   printf("dNdEta:    ");
 178   for (int i=nCentBins;i--;) printf("%.2f,",dNdEta[i]); printf("\n");
 179   printf("dNdEtaErr: ");
 180   for (int i=nCentBins;i--;) printf("%.2f,",dNdEtaErr[i]); printf("\n");
 181
 182 }
 183
 184 //_____________________________________________________________________
 185 Bool_t PrepareHistos(int bin, TList* lst, Bool_t isMC)
 186 {
 187   // fill standard histos for given bin
 188   //
 189   char buffn[500];
 190   char bufft[500];
 191   //
 192   double cutBgMin,cutBgMax;
 193   double cutSgMin,cutSgMax;
 194   //
 195   if (useShapeType==kNormShapeDist) {
 196     cutBgMin = kWDistBgTailMin;
 197     cutBgMax = kWDistBgTailMax;
 198     //
 199     cutSgMin = 0;
 200     cutSgMax = kWDistSgCut;
 201   }
 202   else {
 203     cutBgMin =  kdPhiBgTailMin;
 204     cutBgMax =  kdPhiBgTailMax;
 205     //
 206     cutSgMin =  0;
 207     cutSgMax =  kdPhiSgCut;
 208   }
 209   //
 210   const char* zeCut = "ZvEtaCutT";
 211   TObjArray* res = &resArr;
 212   int shift = bin*kNHistos + (isMC ? kMCShift : 0);
 213   //
 214   // histo for "data" Z vs Eta with default signal cut
 215   TH2* hRawDtCut = (TH2*) FindObject(bin,HName("Data",zeCut),lst);
 216   if (!hRawDtCut) return kFALSE;
 217   sprintf(buffn,"bin%d_%s_RawWithCut",bin,isMC ? "mc":"dt");
 218   sprintf(bufft,"bin%d %s Raw Data with cut on tracklets",bin,isMC ? "mc":"dt");
 219   hRawDtCut = (TH2*)hRawDtCut->Clone(buffn);
 220   hRawDtCut->SetTitle(bufft);
 221   res->AddAtAndExpand(hRawDtCut, kRawDtCut+shift);
 222   //
 223   int nbEta = hRawDtCut->GetXaxis()->GetNbins();
 224   int nbZV  = hRawDtCut->GetYaxis()->GetNbins();
 225   //
 226   // "Data - Est.Bg" histo with cut on tails where we look for signal
 227   sprintf(buffn,"bin%d_%s_SignalWithCut",bin,isMC ? "mc":"dt");
 228   sprintf(bufft,"bin%d %s Signal (raw-bg) with cut on tracklets",bin,isMC ? "mc":"dt");
 229   TH2* hSignalEst = (TH2F*)hRawDtCut->Clone(buffn);
 230   hSignalEst->SetTitle(bufft);
 231   res->AddAtAndExpand(hSignalEst, kSignalEst+shift);
 232   //
 233   // "Data - MC.Bg" histo with cut on tails where we look for signal
 234   TH2* hSignalEstMC = 0;
 235   if (isMC) {
 236     sprintf(buffn,"bin%d_%s_SignalWithCut_bgMCLabels",bin,isMC ? "mc":"dt");
 237     sprintf(bufft,"bin%d %s Signal (raw-bg_MCLabels) with cut on tracklets",bin,isMC ? "mc":"dt");
 238     hSignalEstMC = (TH2F*)hRawDtCut->Clone(buffn);
 239     hSignalEstMC->SetTitle(bufft);
 240     res->AddAtAndExpand(hSignalEstMC, kSignalEstMC+shift);
 241   }
 242   //
 243   // Estimated background in the cut range
 244   sprintf(buffn,"bin%d_%s_BgEst",bin,isMC ? "mc":"dt");
 245   sprintf(bufft,"bin%d %s Estimated Bg",bin,isMC ? "mc":"dt");
 246   TH2* hBgEst = (TH2*) FindObject(bin,HName(useBgType.Data(),zeCut),lst);
 247   if (!hBgEst) return kFALSE;
 248   hBgEst = (TH2*)hBgEst->Clone(buffn);
 249   hBgEst->SetTitle(bufft);
 250   res->AddAtAndExpand(hBgEst,kBgEst +shift);
 251   //
 252   // special feature: use MC Labels bg as a shape instead of generated bg
 253   if (useMCLB/* && !isMC*/) {
 254     TString nm  = hBgEst->GetName();   nm  += "_MCLB";
 255     TString tit = hBgEst->GetTitle();  tit += "_MCLB";
 256     TH2* hBMCLB = (TH2*) FindObject(bin,HName("Comb",zeCut),listMC);
 257     if (!hBMCLB) return kFALSE;
 258     hBMCLB = (TH2F*)hBMCLB->Clone(nm.Data());
 259     hBMCLB->SetTitle(tit.Data());
 260     delete hBgEst;
 261     hBgEst = hBMCLB;
 262     res->AddAtAndExpand(hBgEst,kBgEst +shift);
 263   }
 264   //
 265   // 1-beta for "data" = (Data_cut - Bg_cut) / Data_cut
 266   sprintf(buffn,"bin%d_%s_1mBeta",bin,isMC ? "mc":"dt");
 267   sprintf(bufft,"bin%d %s 1-#beta with estimated bg",bin,isMC ? "mc":"dt");
 268   TH2* h1mBeta = (TH2*)hBgEst->Clone(buffn);
 269   h1mBeta->SetTitle(bufft);
 270   h1mBeta->Reset();
 271   res->AddAtAndExpand(h1mBeta, k1MBeta+shift);
 272   //
 273   // If MC labels info is provided, prepare 1-beta with MC bg
 274   TH2* h1mBetaMC = 0;  // 1-beta for MC with bg from labels
 275   TH2* hBgMC = 0;      // bg from MC labels
 276   if (isMC) {
 277     sprintf(buffn,"bin%d_%s_BgMC",bin,isMC ? "mc":"dt");
 278     sprintf(bufft,"bin%d %s Bg from MC labels",bin,isMC ? "mc":"dt");
 279     hBgMC = (TH2*) FindObject(bin,HName("Comb",zeCut),listMC);
 280     if (!hBgMC) return kFALSE;
 281     hBgMC = (TH2F*)hBgMC->Clone(buffn);
 282     hBgMC->SetTitle(bufft);
 283     res->AddAtAndExpand(hBgMC, kBgMC+shift);
 284     //
 285     sprintf(buffn,"bin%d_%s_h1mBetaMC",bin,isMC ? "mc":"dt");
 286     sprintf(bufft,"bin%d %s 1-#beta with bg. from MC labels",bin,isMC ? "mc":"dt");
 287     h1mBetaMC = (TH2F*) hBgMC->Clone(buffn);
 288     h1mBetaMC->SetTitle(bufft);
 289     res->AddAtAndExpand(h1mBetaMC, k1MBetaMC+shift);
 290     h1mBetaMC->Divide(hRawDtCut);
 291     for (int ib0=1;ib0<=nbEta;ib0++)
 292       for (int ib1=1;ib1<=nbZV;ib1++)
 293         h1mBetaMC->SetBinContent(ib0,ib1, 1.- h1mBetaMC->GetBinContent(ib0,ib1));
 294     //
 295     hSignalEstMC->Add(hBgMC,-1);
 296   }
 297   //
 298   // uncut w.distance or dphi distribution for data
 299   TH1* hDstDt = (TH1*) FindObject(bin,HName("Data",useShapeType==kNormShapeDist ? "WDist":"DPhiS"),lst);
 300   if (!hDstDt) return kFALSE;
 301   sprintf(buffn,"bin%d_%s_DistRawData",bin,isMC ? "mc":"dt");
 302   sprintf(bufft,"bin%d %s Raw Distance for Data",bin,isMC ? "mc":"dt");
 303   hDstDt = (TH1*) hDstDt->Clone(buffn);
 304   hDstDt->SetTitle(bufft);
 305   double nrmDst,dumErr = 0;
 306   Integrate(hDstDt, cutBgMin,cutBgMax, nrmDst, dumErr);
 307   hDstDt->Scale(1./nrmDst);
 308   res->AddAtAndExpand(hDstDt, kDataDist+shift);
 309   //
 310   // uncut w.distance or dphi distribution for generated bg
 311   TH1* hDstBg = (TH1*) FindObject(bin,HName(useBgType.Data(),useShapeType==kNormShapeDist ? "WDist":"DPhiS"),lst);
 312   if (!hDstBg) return kFALSE;
 313   sprintf(buffn,"bin%d_%s_DistRawGenBgNorm",bin,isMC ? "mc":"dt");
 314   sprintf(bufft,"bin%d %s Raw Distance for Gen.Bg. Normalized to data",bin,isMC ? "mc":"dt");
 315   hDstBg = (TH1*) hDstBg->Clone(buffn);
 316   hDstBg->SetTitle(bufft);
 317   hDstBg->Scale(1./nrmDst);
 318   //  res->AddAtAndExpand(hDstBg, kBgDist+shift);
 319   //
 320   // uncut w.distance or dphi distribution for comb. MC labels
 321   TH1* hDstBgMC = 0;
 322   if (isMC) {
 323     hDstBgMC = (TH1*) FindObject(bin,HName("Comb",useShapeType==kNormShapeDist ? "WDist":"DPhiS"),lst);
 324     if (!hDstBgMC) return kFALSE;
 325     sprintf(buffn,"bin%d_%s_DistBgMC",bin,isMC ? "mc":"dt");
 326     sprintf(bufft,"bin%d %s Bg. Distance from MC labels",bin,isMC ? "mc":"dt");
 327     hDstBgMC = (TH1*) hDstBgMC->Clone(buffn);
 328     hDstBgMC->SetTitle(bufft);
 329     hDstBgMC->Scale(1./nrmDst);
 330     res->AddAtAndExpand(hDstBgMC, kBgMCDist+shift);
 331   }
 332   //
 333   // fill 1-beta matrix
 334   double scl,sclE;
 335   hDstBg = NormalizeBg(hDstDt,hDstBg,scl,sclE);  // get rescaling factor for bg. from tails comparison
 336   res->AddAtAndExpand(hDstBg, kBgDist+shift);
 337   double bgVal,bgErr;
 338   double dtVal,dtErr;
 339   // integral in the range where we look for signal
 340   Integrate(hDstBg, cutSgMin, cutSgMax, bgVal, bgErr);
 341   Integrate(hDstDt, cutSgMin, cutSgMax, dtVal, dtErr);
 342   double sclb,sclbErr;
 343   GetRatE(bgVal,bgErr, dtVal, dtErr,sclb,sclbErr);
 344   //  hDstBg->Scale(1./nrmDst);
 345   //
 346   // finalize estimated bg and signal matrices
 347   hBgEst->Scale(scl);
 348   //
 349   hSignalEst->Add(hBgEst,-1);
 350   //
 351   // finalize 1-beta
 352   for (int ib0=1;ib0<=nbEta;ib0++) { // eta
 353     for (int ib1=1;ib1<=nbZV;ib1++) { // zv
 354       //      printf("Bin %d %d\n",ib0,ib1);
 355       double bg  = hBgEst->GetBinContent(ib0,ib1);
 356       double bgE = hBgEst->GetBinError(ib0,ib1);
 357       double dt  = hRawDtCut->GetBinContent(ib0,ib1);
 358       double dtE = hRawDtCut->GetBinError(ib0,ib1);
 359       double beta,betaE;
 360       GetRatE(bg,bgE,dt,dtE, beta,betaE );
 361       h1mBeta->SetBinContent(ib0,ib1,1.-beta);
 362       h1mBeta->SetBinError(ib0,ib1,betaE);
 363       //
 364     }
 365   }
 366   //
 367   if (isMC) {
 368     // prepare MC primary signal histo
 369     sprintf(buffn,"bin%d_zvEtaPrimMC",bin);
 370     sprintf(bufft,"bin%d MC True signal Zv vs Eta",bin);
 371     TH2F* mcPrim = (TH2F*)FindObject(bin,"zvEtaPrimMC", lst );
 372     mcPrim = (TH2F*) mcPrim->Clone(buffn);
 373     mcPrim->SetTitle(bufft);
 374     res->AddAtAndExpand(mcPrim, kMCPrim + shift);
 375   }
 376   //
 377   return kTRUE;
 378 }
 379
 380 //_____________________________________________________________________
 381 void ProcessHistos(int bin)
 382 {
 383   //
 384   int shift = bin*kNHistos;
 385   //
 386   TString prefN = "bin"; prefN += bin; prefN+="_";
 387   TString prefT = "bin"; prefT += bin; prefT+=" ";
 388   TObjArray* res = &resArr;
 389   // build alpha matrix
 390   TH2* halp = (TH2*)res->At(shift + kMCShift + kMCPrim);
 391   halp = (TH2*) halp->Clone(prefN+"Alpha");
 392   halp->SetTitle(prefN+"#alpha");
 393   halp->Divide( (TH2*)res->At(shift + kMCShift + k1MBeta) );
 394   halp->Divide( (TH2*)res->At(shift + kMCShift + kRawDtCut) );
 395   res->AddAtAndExpand(halp, shift + kAlpha);
 396   //
 397   // build alpha matrix with MC labels bg
 398   TH2* halpMC = (TH2*)res->At(shift + kMCShift + kMCPrim);
 399   halpMC = (TH2*) halpMC->Clone(prefN + "AlphaMC");
 400   halpMC->SetTitle(prefT + "#alpha MC labels");
 401   halpMC->Divide( (TH2*)res->At(shift + kMCShift + k1MBetaMC) );
 402   halpMC->Divide( (TH2*)res->At(shift + kMCShift + kRawDtCut) );
 403   res->AddAtAndExpand(halpMC, shift + kAlphaMC);
 404   //
 405   // build corrected signal
 406   TH2* hsigCorr = (TH2*)res->At(shift + kSignalEst);
 407   hsigCorr = (TH2*) hsigCorr->Clone(prefN + "SignalEstCorr");
 408   hsigCorr->SetTitle(prefT + "Corrected Signal");
 409   hsigCorr->Multiply( halp );
 410   res->AddAtAndExpand(hsigCorr, shift + kSigCorr);
 411   //
 412   TH1* hsigCorrX = hsigCorr->ProjectionX("DataCorrSignalX");
 413   hsigCorrX->Scale(1./hsigCorr->GetBinWidth(1));
 414   TF1* pl0 = new TF1("pl0","pol0");
 415   pl0->SetParameter(0,hsigCorr->GetMinimum());
 416   hsigCorrX->Fit(pl0,"q0","",-kEtaFitRange,kEtaFitRange);
 417   double fval = pl0->GetParameter(0);
 418   double ferr = pl0->GetParError(0);
 419   delete hsigCorrX;
 420   dNdEta[bin]    = fval;
 421   dNdEtaErr[bin] = ferr;
 422   printf("Bin %d: dN/d#eta_{|#eta|<0.5} = %.2f  %.2f\n",bin, fval,ferr);
 423   //
 424 }
 425
 426 void PlotResults()
 427 {
 428   TString psnm1 = figDir; psnm1 += "/"; psnm1 += uniqueName;
 429   psnm1 += "_"; psnm1 += nCentBins; psnm1+= "bins_";
 430   psnm1 += outStr; psnm1 += ".ps";
 431   TString psnm0 = psnm1.Data();
 432   psnm0 += "[";
 433   TString psnm2 = psnm1.Data();
 434   psnm2 += "]";
 435   //
 436   TH1* hstat = (TH1*)FindObject(-1,kHStatName,listDt,kFALSE);
 437   //
 438   TH1* hbins = (TH1*)FindObject(-1,"EvCentr",listDt,kFALSE);
 439   //
 440   if (!canvFin) canvFin = new TCanvas("canvFin", "canvFin",0,50,700,1000);
 441   canvFin->Clear();
 442   //
 443   canvFin->Print(psnm0.Data());
 444   //
 445   canvFin->Divide(1,2);
 446   canvFin->cd(1);
 447   gPad->SetLeftMargin(0.15);
 448   TGraphErrors* grp = new TGraphErrors(nCentBins);
 449   for (int i=0;i<nCentBins;i++) {
 450     grp->SetPoint(i,hbins->GetBinCenter(i+1),dNdEta[i]);
 451     grp->SetPointError(i,hbins->GetBinWidth(i+1)/2,dNdEtaErr[i]);
 452   }
 453   grp->SetMarkerStyle(20);
 454   grp->SetMarkerColor(kRed);
 455   grp->SetLineColor(kRed);
 456   grp->SetMinimum(1e-6);
 457   grp->Draw("ap");
 458   grp->GetXaxis()->SetTitle("Centrality Variable");
 459   grp->GetYaxis()->SetTitle("dN/d#eta_{|#eta|<0.5}");
 460   grp->GetYaxis()->SetTitleOffset(1.8);
 461   gPad->SetGrid(1,1);
 462   //
 463   canvFin->cd(2);
 464   gPad->SetLeftMargin(0.15);
 465   hbins->Draw();
 466   hbins->SetMinimum(1e-6);
 467   hbins->SetMarkerStyle(20);
 468   hbins->SetMarkerColor(kRed);
 469   hbins->SetLineColor(kRed);
 470   hbins->GetYaxis()->SetTitle("accepted events");
 471   hbins->GetYaxis()->SetTitleOffset(1.8);
 472   gPad->SetGrid(1,1);
 473   //
 474   canvFin->cd(0);
 475   canvFin->Print(psnm1.Data());
 476   //
 477   const TArrayD &binArr = *hbins->GetXaxis()->GetXbins();
 478   //
 479   for (int i=0;i<nCentBins;i++) {
 480     //
 481     sprintf(outTitle,"%s, %d<C_%s<%d, %.1f<#eta<%.1f, %.1f<Z_{V}<%.1f,  Bg.:%s, UseMCLB=%d, CutVar:%s, |sig|<%.2f, %.2f<|bg.nrm|<%.2f",
 482             uniqueName.Data(),
 483             (int)binArr[i],
 484             hstat->GetXaxis()->GetBinLabel(kCentVar),
 485             (int)binArr[i+1],
 486             hstat->GetBinContent(kEtaMin)/myMergeFactor,
 487             hstat->GetBinContent(kEtaMax)/myMergeFactor,
 488             hstat->GetBinContent(kZVMin)/myMergeFactor,
 489             hstat->GetBinContent(kZVMax)/myMergeFactor,
 490             useBgType.Data(),
 491             useMCLB,
 492             useShapeType==kNormShapeDist ? "#Delta":"#Delta#varphi-#delta_{#varphi}",
 493             useShapeType==kNormShapeDist ? kWDistSgCut:kdPhiSgCut,
 494             useShapeType==kNormShapeDist ? kWDistBgTailMin : kdPhiBgTailMin,
 495             useShapeType==kNormShapeDist ? kWDistBgTailMax : kdPhiBgTailMax
 496             );
 497     //
 498     PlotDNDEta(i);
 499     canvFin->Print(psnm1.Data());
 500     PlotAlphaBeta(i);
 501     canvFin->Print(psnm1.Data());
 502   }
 503   PlotSpecies();
 504   canvFin->Print(psnm1.Data());
 505   //
 506   canvFin->Print(psnm2.Data());
 507 }
 508
 509 void PlotDNDEta(int bin)
 510 {
 511   //
 512   TObjArray *res = &resArr;
 513   TString prefN = "bin"; prefN += bin; prefN+="_";
 514   int shift = bin*kNHistos;
 515   //
 516   char buff[1000];
 517   // eta range
 518   gStyle->SetOptFit(0);
 519   gStyle->SetOptStat(0);
 520   gStyle->SetOptTitle(0);
 521   double mn = 1e6,mx = -1e6;
 522   if (!canvFin) canvFin = new TCanvas("canvFin", "canvFin",0,50,700,1000);
 523   canvFin->Clear();
 524   canvFin->Divide(1,2);
 525   canvFin->cd(1);
 526   gPad->SetLeftMargin(0.15);
 527   //
 528   // corrected data
 529   TString nms =  prefN;
 530   nms += "DataCorrSignal";
 531   nms += "_";
 532   nms += uniqueName;
 533   TH1* hsigCorr = ((TH2F*)res->At(shift + kSigCorr))->ProjectionX(nms.Data());
 534   SetHStyle(hsigCorr,kRed,20,1.0);
 535   hsigCorr->Scale(1./hsigCorr->GetBinWidth(1));
 536   hsigCorr->Draw();
 537   mx = TMath::Max(mx, hsigCorr->GetMaximum());
 538   mn = TMath::Min(mn, hsigCorr->GetMinimum());
 539   char ftres[1000];
 540   sprintf(ftres,"dN/d#eta_{|#eta|<%.1f} = %.1f #pm %.1f",kEtaFitRange,dNdEta[bin],dNdEtaErr[bin]);
 541   TLatex *txfit = new TLatex(-0.2,hsigCorr->GetMinimum()*0.9, ftres);
 542   txfit->SetTextSize(0.04);
 543   txfit->Draw();
 544   resDnDeta.AddAtAndExpand( hsigCorr, bin );
 545   hsigCorr->SetDirectory(0);
 546   //
 547   // raw data
 548   TH1* hraw = ((TH2F*)res->At(shift+kRawDtCut))->ProjectionX(prefN+"DataRaw");
 549   SetHStyle(hraw,kRed,21,1.0);
 550   hraw->Scale(1./hraw->GetBinWidth(1));
 551   hraw->Draw("same");
 552   mn = TMath::Min(mn, hraw->GetMinimum());
 553   mx = TMath::Max(mx, hraw->GetMaximum());
 554   //
 555   // raw data bg sub
 556   TH1* hraws = ((TH2F*)res->At(shift+kSignalEst))->ProjectionX(prefN+"DataRawSub");
 557   SetHStyle(hraws,kRed,23,1.0);
 558   hraws->Scale(1./hraw->GetBinWidth(1));
 559   hraws->Draw("same");
 560   mn = TMath::Min(mn, hraw->GetMinimum());
 561   mx = TMath::Max(mx, hraw->GetMaximum());
 562   //
 563   // bg
 564   TH1* hbg = ((TH2F*)res->At(shift+kBgEst))->ProjectionX(prefN+"BgEst");
 565   SetHStyle(hbg,kMagenta,22,1.0);
 566   hbg->Scale(1./hbg->GetBinWidth(1));
 567   hbg->Draw("same");
 568   mn = TMath::Min(mn, hbg->GetMinimum());
 569   mx = TMath::Max(mx, hbg->GetMaximum());
 570   //
 571   // mc part ----------------------------
 572   // raw data
 573   TH1* hrawMC = ((TH2F*)res->At(shift+kRawDtCut+kMCShift))->ProjectionX(prefN+"DataRawMC");
 574   SetHStyle(hrawMC,kBlue,24,1.0);
 575   hrawMC->Scale(1./hrawMC->GetBinWidth(1));
 576   hrawMC->Draw("same");
 577   mn = TMath::Min(mn, hrawMC->GetMinimum());
 578   mx = TMath::Max(mx, hrawMC->GetMaximum());
 579   //
 580   // raw data bg sub
 581   TH1* hrawsMC = ((TH2F*)res->At(shift+kSignalEst+kMCShift))->ProjectionX(prefN+"DataRawSubMC");
 582   SetHStyle(hrawsMC,kBlue,26,1.0);
 583   hrawsMC->Scale(1./hrawMC->GetBinWidth(1));
 584   hrawsMC->Draw("same");
 585   mn = TMath::Min(mn, hrawMC->GetMinimum());
 586   mx = TMath::Max(mx, hrawMC->GetMaximum());
 587   //
 588   // raw data bgMClabels sub
 589   TH1* hrawsMCLB = ((TH2F*)res->At(shift+kSignalEstMC+kMCShift))->ProjectionX(prefN+"DataRawSubMCLB");
 590   SetHStyle(hrawsMCLB,kGreen+2,30,1.0);
 591   hrawsMCLB->Scale(1./hrawsMCLB->GetBinWidth(1));
 592   hrawsMCLB->Draw("same");
 593   mn = TMath::Min(mn, hrawsMCLB->GetMinimum());
 594   mx = TMath::Max(mx, hrawsMCLB->GetMaximum());
 595   //
 596   // bg est
 597   TH1* hbgMCEst = ((TH2F*)res->At(shift+kBgEst+kMCShift))->ProjectionX(prefN+"BgEstMC");
 598   SetHStyle(hbgMCEst,kBlue,26,1.0);
 599   hbgMCEst->Scale(1./hbgMCEst->GetBinWidth(1));
 600   hbgMCEst->Draw("same");
 601   mn = TMath::Min(mn, hbgMCEst->GetMinimum());
 602   mx = TMath::Max(mx, hbgMCEst->GetMaximum());
 603   //
 604   // bg MC
 605   TH1* hbgMC = ((TH2F*)res->At(shift+kBgMC+kMCShift))->ProjectionX(prefN+"BgMC");
 606   SetHStyle(hbgMC,kGreen+2,25,1.0);
 607   hbgMC->Scale(1./hbgMC->GetBinWidth(1));
 608   hbgMC->Draw("same");
 609   mn = TMath::Min(mn, hbgMC->GetMinimum());
 610   mx = TMath::Max(mx, hbgMC->GetMaximum());
 611   //
 612   mn = 0;
 613   hsigCorr->SetMinimum(mn);
 614   hsigCorr->SetMaximum(mx + 0.4*(mx-mn));
 615   gPad->Modified();
 616   //
 617   TLegend *legDnDeta = new TLegend(0.15,0.75, 0.45,0.95);
 618   legDnDeta->SetFillColor(kWhite);
 619   legDnDeta->SetHeader("Data");
 620   legDnDeta->AddEntry(hsigCorr,"Corrected","pl");
 621   legDnDeta->AddEntry(hraw,    "Reconstructed","pl");
 622   sprintf(buff,"Reconstructed - Bckg.%s.",useBgType.Data());
 623   legDnDeta->AddEntry(hraws,   buff,"pl");
 624   sprintf(buff,"Background %s.",useBgType.Data());
 625   legDnDeta->AddEntry(hbg,     buff,"pl");
 626   legDnDeta->Draw();
 627   //
 628   TLegend *legDnDetaMC = new TLegend(0.60,0.72, 0.95,0.95);
 629   legDnDetaMC->SetFillColor(kWhite);
 630   legDnDetaMC->SetHeader("MC");
 631   legDnDetaMC->AddEntry(hrawMC,    "Reconstructed","pl");
 632   sprintf(buff,"Reconstructed - Bckg.%s.",useBgType.Data());
 633   legDnDetaMC->AddEntry(hrawsMC,   buff,"pl");
 634   sprintf(buff,"Reconstructed - Bckg.%s.","MC.Labels");
 635   legDnDetaMC->AddEntry(hrawsMCLB,   buff,"pl");
 636   sprintf(buff,"Background %s.",useBgType.Data());
 637   legDnDetaMC->AddEntry(hbgMCEst,     buff,"pl");
 638   sprintf(buff,"Background %s.","MC Labels");
 639   legDnDetaMC->AddEntry(hbgMC,     buff,"pl");
 640   //
 641   legDnDetaMC->Draw();
 642   //
 643   gPad->SetGrid(1.1);
 644   gPad->Modified();
 645   AddLabel(outTitle,0.1,0.97, kBlack,0.02);
 646   //
 647   canvFin->cd();
 648   //
 649   //---------------- dsitributions
 650   canvFin->cd(2);
 651   //
 652   TH1* mcdst = (TH1*)res->At(shift+kDataDist+kMCShift);
 653   TH1* mcdstbg = (TH1*)res->At(shift+kBgDist+kMCShift);
 654   TH1* mcdstbgLB = (TH1*)res->At(shift+kBgMCDist+kMCShift);
 655   TH1* dtdst = (TH1*)res->At(shift+kDataDist);
 656   TH1* dtdstbg = (TH1*)res->At(shift+kBgDist);
 657   //
 658   TH1* mcDstN     = (TH1*)FindObject(bin,HName("Data", useShapeType==kNormShapeDist ? "WDist":"DPhiS"), listMC );
 659   TH1* mcDstSec   = (TH1*)FindObject(bin,HName("Sec",  useShapeType==kNormShapeDist ? "WDist":"DPhiS"), listMC );
 660   TH1* mcDstCombU = (TH1*)FindObject(bin,HName("CombU",useShapeType==kNormShapeDist ? "WDist":"DPhiS"), listMC );
 661   TH1* mcDstCombC = (TH1*)FindObject(bin,HName("Comb", useShapeType==kNormShapeDist ? "WDist":"DPhiS"), listMC );
 662   //
 663   double scl,sclE;
 664   mcDstN = NormalizeBg(mcdst,mcDstN,scl,sclE);
 665   mcDstSec->Scale(scl);
 666   mcDstCombU->Scale(scl);
 667   mcDstCombC->Scale(scl);
 668   mcDstCombC->Add(mcDstCombU,-1);
 669
 670   dtdst->Draw("");
 671   gPad->Modified();
 672   dtdst->GetXaxis()->SetLabelSize(0.03);
 673   dtdst->GetXaxis()->SetTitleSize(0.03);
 674   dtdst->GetXaxis()->SetTitleOffset(2);
 675   dtdstbg->Draw("same");
 676   //
 677   mcdst->Draw("same");
 678   mcDstSec->Draw("same");
 679   mcdstbgLB->Draw("same");
 680   mcdstbg->Draw("same");
 681   mcDstCombC->Draw("same");
 682   //
 683
 684   SetHStyle(mcdst,kBlue, 25,0.7);
 685   SetHStyle(mcdstbgLB,kGreen, 7/*32*/,0.5);
 686   SetHStyle(mcdstbg,kCyan, 7/*26*/,0.5);
 687   SetHStyle(mcDstCombC,kGreen+2, 21,0.7);
 688   SetHStyle(mcDstSec,kBlue+2, 22,0.7);
 689   //
 690   SetHStyle(dtdst,kRed, 20,0.7);
 691   SetHStyle(dtdstbg,kBlue, 34,0.7);
 692   //
 693   double vmcTot,vmcTotE;
 694   double vmcSec,vmcSecE, ratSec,ratSecE;
 695   double vmcCmbEst,vmcCmbEstE, ratCmbEst,ratCmbEstE;
 696   double vmcCmb,vmcCmbE, ratCmb,ratCmbE;
 697   double vmcCmbC,vmcCmbCE, ratCmbC,ratCmbCE;
 698   double cutSgMin,cutSgMax;
 699   double cutBgMin,cutBgMax;
 700   if (useShapeType==kNormShapeDist) {
 701     cutSgMin = 0;
 702     cutSgMax = kWDistSgCut;
 703     cutBgMin = kWDistBgTailMin;
 704     cutBgMax = kWDistBgTailMax;
 705   }
 706   else {
 707     cutSgMin =  0;
 708     cutSgMax =  kdPhiSgCut;
 709     cutBgMin =  kdPhiBgTailMin;
 710     cutBgMax =  kdPhiBgTailMax;
 711   }
 712   Integrate(mcdst,    cutSgMin,cutSgMax, vmcTot,vmcTotE);
 713   Integrate(mcDstSec, cutSgMin,cutSgMax, vmcSec,vmcSecE);
 714   GetRatE(vmcSec,vmcSecE, vmcTot,vmcTotE, ratSec,ratSecE);
 715   //
 716   Integrate(mcdstbgLB, cutSgMin,cutSgMax, vmcCmb,vmcCmbE);
 717   GetRatE(vmcCmb,vmcCmbE, vmcTot,vmcTotE, ratCmb,ratCmbE);
 718   //
 719   Integrate(mcdstbg,  cutSgMin,cutSgMax, vmcCmbEst,vmcCmbEstE);
 720   GetRatE(vmcCmbEst,vmcCmbEstE, vmcTot,vmcTotE, ratCmbEst,ratCmbEstE);
 721   //
 722   Integrate(mcDstCombC,  cutSgMin,cutSgMax, vmcCmbC,vmcCmbCE);
 723   GetRatE(vmcCmbC,vmcCmbCE, vmcTot,vmcTotE, ratCmbC,ratCmbCE);
 724   //
 725   double vdtTot,vdtTotE;
 726   double vdtBg,vdtBgE, ratdtBg,ratdtBgE;
 727   //
 728   Integrate(dtdst,    cutSgMin,cutSgMax, vdtTot,vdtTotE);
 729   Integrate(dtdstbg,  cutSgMin,cutSgMax, vdtBg,vdtBgE);
 730   GetRatE( vdtBg,vdtBgE,  vdtTot,vdtTotE, ratdtBg,ratdtBgE);
 731   //
 732   //
 733   double dmn = mcdst->GetMinimum();
 734   double dmx = mcdst->GetMaximum();
 735   TLine *ln = new TLine(cutSgMax, dmn, cutSgMax, dmx);
 736   ln->SetLineColor(kBlack);
 737   ln->Draw();
 738   TLine *lnc = new TLine(cutBgMin, dmn, cutBgMin, dmx);
 739   lnc->SetLineColor(kRed);
 740   lnc->Draw();
 741   if (useShapeType==kNormShapeDPhi) {
 742     ln = new TLine(-cutSgMax, dmn, -cutSgMax, dmx);
 743     ln->SetLineColor(kBlack);
 744     ln->Draw();
 745     //
 746     lnc = new TLine(-cutBgMin, dmn, -cutBgMin, dmx);
 747     lnc->SetLineColor(kRed);
 748     lnc->Draw();
 749   }
 750   //
 751   TLegend *legDstMC1 = new TLegend(0.60,0.72, 0.95,0.95);
 752   legDstMC1->SetFillColor(kWhite);
 753
 754   //
 755   legDstMC1->AddEntry(dtdst,    "Data raw","pl");
 756   sprintf(buff,"Data Comb. %s. : %.1f%%",useBgType.Data(),ratdtBg*100);
 757   legDstMC1->AddEntry(dtdstbg,  buff,"pl");
 758   //
 759
 760
 761   legDstMC1->AddEntry(mcdst,    "MC raw","pl");
 762   sprintf(buff,"MC Comb. %s. : %.1f%%",useBgType.Data(),ratCmbEst*100);
 763   legDstMC1->AddEntry(mcdstbg,  buff,"pl");
 764   //
 765   sprintf(buff,"MC Comb. %s. : %.1f%%","MC Lbl.",ratCmb*100);
 766   legDstMC1->AddEntry(mcdstbgLB,  buff,"pl");
 767
 768   sprintf(buff,"MC Comb.Uncorr %s. : %.1f%%","MC Lbl.",ratCmbC*100);
 769   legDstMC1->AddEntry(mcDstCombC,  buff,"pl");
 770
 771   sprintf(buff,"MC Sec.   : %.1f%%",ratSec*100);
 772   legDstMC1->AddEntry(mcDstSec,  buff,"pl");
 773
 774   legDstMC1->Draw();
 775
 776   if (useShapeType==kNormShapeDist) gPad->SetLogx();
 777   gPad->SetLogy();
 778   gPad->SetGrid(1,1);
 779   gPad->Modified();
 780   //
 781   canvFin->cd();
 782   //
 783   if (creatDnDEtaCMacro) {
 784     sprintf(buff,"%s/%s_b%d_dNdEta_%s",figDir,uniqueName.Data(),bin,outStr);
 785     SaveCanvas(canvFin,buff,"cg");
 786   }
 787   //
 788 }
 789 //
 790 void PlotAlphaBeta(int bin)
 791 {
 792   char buff[1000];
 793   int shift = bin*kNHistos;
 794   gStyle->SetOptFit(0);
 795   gStyle->SetOptStat(0);
 796   gStyle->SetOptTitle(0);
 797   TObjArray* res = &resArr;
 798   //------------------------------------------------------
 799   if (!canvFin) canvFin = new TCanvas("canvFin","canvFin",10,10,700,1000);
 800   canvFin->Clear();
 801   canvFin->Divide(2,3,0.01,0.01);
 802   canvFin->cd(1);
 803   TH1* dtBet = (TH1*)res->At(shift + k1MBeta);
 804   TH1* mcBet = (TH1*)res->At(shift + k1MBeta+kMCShift);
 805   TH1* mcBetLB = (TH1*)res->At(shift + k1MBetaMC+kMCShift);
 806   double mn,mx,mnt,mxt;
 807   GetRealMinMax(dtBet,mn,mx);
 808   GetRealMinMax(mcBet,mnt,mxt);
 809   if (mnt<mn) mn = mnt;
 810   if (mxt>mx) mx = mxt;
 811   GetRealMinMax(mcBetLB,mnt,mxt);
 812   if (mnt<mn) mn = mnt;
 813   if (mxt>mx) mx = mxt;
 814   //
 815   dtBet->SetMinimum(mn - 0.05*(mx-mn));
 816   dtBet->SetMaximum(mx + 0.05*(mx-mn));
 817   mcBet->SetMinimum(mn - 0.05*(mx-mn));
 818   mcBet->SetMaximum(mx + 0.05*(mx-mn));
 819   mcBetLB->SetMinimum(mn - 0.05*(mx-mn));
 820   mcBetLB->SetMaximum(mx + 0.05*(mx-mn));
 821   //
 822   canvFin->cd(1);
 823   gPad->SetRightMargin(0.15);
 824   dtBet->Draw("colz");
 825   AddLabel("#beta Data",0.2,0.95,kBlack,0.04);
 826   gPad->Modified();
 827   dtBet->GetYaxis()->SetTitleOffset(1.4);
 828   TPaletteAxis *p = (TPaletteAxis*)dtBet->FindObject("palette");
 829   if (p) p->SetX1NDC(0.85);
 830   canvFin->cd(2);
 831   gPad->SetRightMargin(0.15);
 832   mcBet->Draw("colz");
 833   AddLabel("#beta MC (bckg.estimated)",0.2,0.95,kBlack,0.04);
 834   gPad->Modified();
 835   mcBet->GetYaxis()->SetTitleOffset(1.4);
 836   p = (TPaletteAxis*)mcBet->FindObject("palette");
 837   if (p) p->SetX1NDC(0.85);
 838   canvFin->cd(3);
 839   gPad->SetRightMargin(0.15);
 840   mcBetLB->Draw("colz");
 841   AddLabel("#beta MC (bckg.from MC labels)",0.2,0.95,kBlack,0.04);
 842   gPad->Modified();
 843   mcBetLB->GetYaxis()->SetTitleOffset(1.4);
 844   p = (TPaletteAxis*)mcBetLB->FindObject("palette");
 845   if (p) p->SetX1NDC(0.85);
 846   //
 847   //------------------------------------------------------
 848   TH1* dtAlp = (TH1*)res->At(shift + kAlpha);
 849   TH1* mcAlp = (TH1*)res->At(shift + kAlphaMC);
 850   GetRealMinMax(dtAlp,mn,mx);
 851   GetRealMinMax(mcAlp,mnt,mxt);
 852   if (mnt<mn) mn = mnt;
 853   if (mxt>mx) mx = mxt;
 854   dtAlp->SetMinimum(mn - 0.05*(mx-mn));
 855   dtAlp->SetMaximum(mx + 0.05*(mx-mn));
 856   mcAlp->SetMinimum(mn - 0.05*(mx-mn));
 857   mcAlp->SetMaximum(mx + 0.05*(mx-mn));
 858   //
 859   canvFin->cd(4);
 860   gPad->SetRightMargin(0.15);
 861   dtAlp->Draw("colz");
 862   AddLabel("#alpha (bckg.estimated)",0.2,0.95,kBlack,0.04);
 863   gPad->Modified();
 864   dtAlp->GetYaxis()->SetTitleOffset(1.4);
 865   TPaletteAxis *pa = (TPaletteAxis*)dtBet->FindObject("palette");
 866   if (pa) pa->SetX1NDC(0.85);
 867   canvFin->cd(5);
 868   gPad->SetRightMargin(0.15);
 869   mcAlp->Draw("colz");
 870   AddLabel("#alpha (bckg.from MC labels)",0.2,0.95,kBlack,0.04);
 871   gPad->Modified();
 872   mcAlp->GetYaxis()->SetTitleOffset(1.4);
 873   pa = (TPaletteAxis*)mcBet->FindObject("palette");
 874   if (pa) pa->SetX1NDC(0.85);
 875   gPad->Modified();
 876   canvFin->cd(6);
 877   AddLabel(outTitle,0.1,0.5, kBlack, 0.02);
 878   //
 879   if (creatAlphaBetaCMacro) {
 880     sprintf(buff,"%s/%sAlphaBeta_%s",figDir,uniqueName.Data(),outStr);
 881     SaveCanvas(canvFin,buff,"cg");
 882   }
 883   //
 884 }
 885
 886 void PlotSpecies()
 887 {
 888   char buff[1000];
 889   gStyle->SetOptFit(0);
 890   gStyle->SetOptStat(0);
 891   gStyle->SetOptTitle(0);
 892   //------------------------------------------------------
 893   TH2F* hSpecPrim  = (TH2F*)FindObject(-1, "pdgPrim", listMC,kFALSE);
 894   TH2F* hSpecSec   = (TH2F*)FindObject(-1, "pdgSec", listMC,kFALSE);
 895   TH2F* hSpecPrimP = (TH2F*)FindObject(-1, "pdgPrimPar", listMC,kFALSE);
 896   TH2F* hSpecSecP  = (TH2F*)FindObject(-1, "pdgSecPar", listMC,kFALSE);
 897   int nbd = hSpecPrim->GetXaxis()->GetNbins();
 898   //
 899   TH1* hSpecPrimAll = hSpecPrim->ProjectionX("specPrimAll",1,nbd+1,"e");
 900   hSpecPrimAll->Scale(100./hSpecPrimAll->Integral());
 901   hSpecPrimAll->GetYaxis()->SetTitle("Fraction,%");
 902   hSpecPrimAll->GetXaxis()->SetLabelSize(0.06);
 903   hSpecPrimAll->GetXaxis()->LabelsOption("v");
 904   //
 905   TH1* hSpecSecAll  = hSpecSec->ProjectionX("specSecAll",1,nbd+1,"e");
 906   hSpecSecAll->Scale(100./hSpecSecAll->Integral());
 907   hSpecSecAll->GetYaxis()->SetTitle("Fraction,%");
 908   hSpecSecAll->GetXaxis()->SetLabelSize(0.05);
 909   //
 910   TH1* hSpecPrimPAll = hSpecPrimP->ProjectionX("specPrimPAll",1,nbd+1,"e");
 911   hSpecPrimPAll->Scale(100./hSpecPrimPAll->Integral());
 912   hSpecPrimPAll->GetYaxis()->SetTitle("Fraction,%");
 913   hSpecPrimPAll->GetXaxis()->SetLabelSize(0.06);
 914   hSpecPrimPAll->GetXaxis()->LabelsOption("v");
 915
 916   //
 917   TH1* hSpecSecPAll  = hSpecSecP->ProjectionX("specSecPAll",1,nbd+1,"e");
 918   hSpecSecPAll->Scale(100./hSpecSecPAll->Integral());
 919   hSpecSecPAll->GetYaxis()->SetTitle("Fraction,%");
 920   hSpecSecPAll->GetXaxis()->SetLabelSize(0.05);
 921   //
 922   int binCut = hSpecPrim->GetXaxis()->FindBin(kWDistSgCut-kEps);
 923   TH1* hSpecPrimSel = hSpecPrim->ProjectionX("specPrimSel",1,binCut,"e");
 924   hSpecPrimSel->Scale(100./hSpecPrimSel->Integral());
 925   hSpecPrimSel->GetYaxis()->SetTitle("Fraction,%");
 926   hSpecPrimSel->GetXaxis()->SetLabelSize(0.05);
 927   //
 928   TH1* hSpecSecSel  = hSpecSec->ProjectionX("specSecSel",1,binCut,"e");
 929   hSpecSecSel->Scale(100./hSpecSecSel->Integral());
 930   hSpecSecSel->GetYaxis()->SetTitle("Fraction,%");
 931   hSpecSecSel->GetXaxis()->SetLabelSize(0.05);
 932   //
 933   TH1* hSpecPrimPSel = hSpecPrimP->ProjectionX("specPrimPSel",1,binCut,"e");
 934   hSpecPrimPSel->Scale(100./hSpecPrimPSel->Integral());
 935   hSpecPrimPSel->GetYaxis()->SetTitle("Fraction,%");
 936   hSpecPrimPSel->GetXaxis()->SetLabelSize(0.05);
 937   //
 938   TH1* hSpecSecPSel  = hSpecSecP->ProjectionX("specSecPSel",1,binCut,"e");
 939   hSpecSecPSel->Scale(100./hSpecSecPSel->Integral());
 940   hSpecSecPSel->GetYaxis()->SetTitle("Fraction,%");
 941   hSpecSecPSel->GetXaxis()->SetLabelSize(0.05);
 942   //
 943   if (!canvFin) canvFin = new TCanvas("canvFin","canvFin",10,10,1100,800);
 944   canvFin->Clear();
 945   canvFin->Divide(1,2,0.01,0.01);
 946   canvFin->cd(1);
 947   hSpecPrimAll->Draw();
 948   SetHStyle(hSpecPrimAll,kBlue,25,1.1);
 949   hSpecPrimSel->Draw("same");
 950   SetHStyle(hSpecPrimSel,kRed,20,1);
 951   //
 952   hSpecSecAll->Draw("same");
 953   SetHStyle(hSpecSecAll,kGreen,32,1.1);
 954   hSpecSecSel->Draw("same");
 955   SetHStyle(hSpecSecSel,kBlack,22,1);
 956   //
 957   TLegend *legPart = new TLegend(0.8,0.72, 0.999,0.999);
 958   legPart->SetFillColor(kWhite);
 959   legPart->SetHeader("Tracklet PDG");
 960   //
 961   legPart->AddEntry(hSpecPrimAll,    "Prim., before #Delta cut","pl");
 962   legPart->AddEntry(hSpecPrimSel,    "Prim., after #Delta cut","pl");
 963   legPart->AddEntry(hSpecSecAll,     "Sec., before #Delta cut","pl");
 964   legPart->AddEntry(hSpecSecSel,     "Sec., after #Delta cut","pl");
 965   //
 966   legPart->Draw();
 967   gPad->SetLogy();
 968   gPad->SetGrid(1,1);
 969   gPad->Modified();
 970   //
 971   canvFin->cd(2);
 972   hSpecPrimPAll->Draw();
 973   SetHStyle(hSpecPrimPAll,kBlue,25,1.1);
 974   hSpecPrimPSel->Draw("same");
 975   SetHStyle(hSpecPrimPSel,kRed,20,1);
 976   //
 977   hSpecSecPAll->Draw("same");
 978   SetHStyle(hSpecSecPAll,kGreen,32,1.1);
 979   hSpecSecPSel->Draw("same");
 980   SetHStyle(hSpecSecPSel,kBlack,22,1);
 981   //
 982   TLegend *legPartP = new TLegend(0.8,0.72, 0.999,0.999);
 983   legPartP->SetFillColor(kWhite);
 984   legPartP->SetHeader("Tracklet Parents PDG");
 985   //
 986   legPartP->AddEntry(hSpecPrimPAll,    "Prim., before #Delta cut","pl");
 987   legPartP->AddEntry(hSpecPrimPSel,    "Prim., after #Delta cut","pl");
 988   legPartP->AddEntry(hSpecSecPAll,     "Sec., before #Delta cut","pl");
 989   legPartP->AddEntry(hSpecSecPSel,     "Sec., after #Delta cut","pl");
 990   //
 991   legPartP->Draw();
 992   gPad->SetLogy();
 993   gPad->SetGrid(1,1);
 994   gPad->Modified();
 995   //
 996   canvFin->cd(1);
 997   //  AddLabel(outTitle,0.1,0.97, kBlack, 0.02);
 998   canvFin->cd();
 999   //
1000   if (creatSpeciesCMacro) {
1001     sprintf(buff,"%s/%sSpecies_%s",figDir,uniqueName.Data(),outStr);
1002     SaveCanvas(canvFin,buff,"cg");
1003   }
1004 }
1005
1006 //______________________________________________________________________
1007 void CropHisto(TH1* histo, int bx0, int bx1, int by0, int by1)
1008 {
1009   // fill 0 to all bins outside defined range
1010   TAxis *xax = histo->GetXaxis();
1011   int nbx = xax->GetNbins();
1012   double vmn=1e16,vmx=-1e16;
1013   if (histo->InheritsFrom(TH2::Class())) {
1014     TAxis *yax = histo->GetYaxis();
1015     int nby = yax->GetNbins();
1016     for (int ix=nbx+2;ix--;) {
1017       for (int iy=nby+2;iy--;) {
1018         if ((ix<bx0||ix>bx1)||(iy<by0||iy>by1)) {
1019           histo->SetBinContent(ix,iy,0);
1020           histo->SetBinError(ix,iy,0);
1021         }
1022         else {
1023           double vl = histo->GetBinContent(ix,iy);
1024           if (vl<vmn) vmn = vl;
1025           if (vl>vmx) vmx = vl;
1026         }
1027       }
1028     }
1029   }
1030   else {
1031     for (int ix=nbx+2;ix--;) {
1032       if ((ix<bx0||ix>bx1)) {
1033         histo->SetBinContent(ix,0);
1034         histo->SetBinError(ix,0);
1035       }
1036       else {
1037         double vl = histo->GetBinContent(ix);
1038         if (vl<vmn) vmn = vl;
1039         if (vl>vmx) vmx = vl;
1040       }
1041     }
1042   }
1043   //
1044   if (vmn==vmx) {
1045     vmn = 0.95*vmn;
1046     vmx = 1.05*vmx;
1047   }
1048   histo->SetMaximum(vmx);
1049   histo->SetMinimum(vmn);
1050 }
1051
1052 //______________________________________________________________________
1053 void CropHisto(TH1* histo, double vx0, double vx1, double vy0, double vy1)
1054 {
1055   // fill 0 to all bins outside defined range
1056   TAxis *xax = histo->GetXaxis();
1057   int bx0,bx1,by0=-1,by1=-1;
1058   bx0 = xax->FindBin(vx0+kEps);
1059   bx1 = xax->FindBin(vx1-kEps);
1060   if (histo->InheritsFrom(TH2::Class())) {
1061     TAxis *yax = histo->GetYaxis();
1062     by0 = yax->FindBin(vy0+kEps);
1063     by1 = yax->FindBin(vy1-kEps);
1064   }
1065   CropHisto(histo,bx0,bx1,by0,by1);
1066 }
1067
1068 //______________________________________________________________________
1069 TH1* NormalizeBg(TH1* dataH, TH1* bgH, double &scl, double &sclE)
1070 {
1071   // match generated bg and data tails, calculate normalization, return normalized bg copy
1072   //
1073   TAxis* xax = dataH->GetXaxis();
1074   int nbtot = xax->GetNbins();
1075   int bgBins[2][2] = {{0}}; // limiting bins for tails integration
1076   Int_t ntails; // 0 for dphi plot, 1 for weighted dist plot
1077   if  (useShapeType == kNormShapeDist) { // only positive tail
1078     bgBins[0][0] = xax->FindBin(kWDistBgTailMin+kEps); // positive tail min bin
1079     bgBins[0][1] = xax->FindBin(kWDistBgTailMax-kEps); // positive tail max bin
1080     ntails = 1;
1081   }
1082   else if (useShapeType == kNormShapeDPhi) {         // both tails
1083     bgBins[0][0] = xax->FindBin( kdPhiBgTailMin+kEps); // positive tail min bin
1084     bgBins[0][1] = xax->FindBin( kdPhiBgTailMax-kEps); // positive tail max bin
1085     bgBins[1][0] = xax->FindBin(-kdPhiBgTailMax+kEps); // negative tail min bin
1086     bgBins[1][1] = xax->FindBin(-kdPhiBgTailMin-kEps); // positive tail max bin
1087     ntails = 2;
1088   }
1089   else {printf("NormalizeBg: unknown shape type %d\n",useShapeType);exit(1);}
1090   printf("NormalizeBg: bins for tails: right: %d:%d / left: %d:%d\n",bgBins[0][0],bgBins[0][1],bgBins[1][0],bgBins[1][1]);
1091   //
1092   double meanR=0,meanRE=0,meanRE2=0;
1093   double meanD=0,meanDE2=0;
1094   double meanB=0,meanBE2=0;
1095   double meanRI=0,meanRIE=0;
1096   for (int itp=0;itp<=ntails;itp++) {
1097     for (int ib=bgBins[itp][0];ib<=bgBins[itp][1];ib++) {
1098       if (ib<1||ib>nbtot) continue;
1099       double vD = dataH->GetBinContent(ib);
1100       double vB = bgH->GetBinContent(ib);
1101       double eD = dataH->GetBinError(ib);
1102       double eB = bgH->GetBinError(ib);
1103       meanD += vD; meanDE2 += eD*eD;
1104       meanB += vB; meanBE2 += eB*eB;
1105       if (vD<=0 || vB<=0 || eD<=0 || eB<=0) continue;
1106       double rat = vD/vB;
1107       double ratE2 = rat*rat*(eD*eD/vD/vD + eB*eB/vB/vB);
1108       meanR   += rat/ratE2; meanRE2 += 1.0/ratE2;
1109     }
1110   }
1111   //
1112   if (meanRE2>0) {
1113     meanR  /= meanRE2;
1114     meanRE2 = 1./meanRE2;
1115     meanRE = TMath::Sqrt(meanRE2);
1116   }
1117   if (meanDE2>0 && meanBE2>0) {
1118     meanRI  = meanD/meanB;
1119     meanRIE =  meanRI*TMath::Sqrt(meanDE2/meanD/meanD + meanBE2/meanB/meanB);
1120   }
1121   printf("NormalizeBg: Tails scaling %s wrt %s: Wgh.Mean:%.4f(%.4f) / Integral:%.4f(%.4f)\n",
1122          bgH->GetName(),dataH->GetName(), meanR,meanRE, meanRI,meanRIE);
1123   printf("NormalizeBg: Select scaling type %s\n",useScaleType==kSclWghMean ? "Wgh.Mean":"Integral");
1124   //
1125   scl  = useScaleType==kSclWghMean ? meanR  : meanRI;
1126   sclE = useScaleType==kSclWghMean ? meanRE : meanRIE;
1127   //
1128   // rescaled bg
1129   char buff[1000];
1130   sprintf(buff,"%s_bgNorm",bgH->GetName());
1131   bgH = (TH1*)bgH->Clone(buff);
1132   sprintf(buff,"%s bgNorm%d %.4f+-%.4f",bgH->GetName(),useScaleType,scl,sclE);
1133   TH1* dumH = (TH1*)bgH->Clone("dummySCL$"); dumH->Reset();
1134   for (int i=1;i<=nbtot;i++) {
1135     dumH->SetBinContent(i,scl);
1136     dumH->SetBinError(i,sclE);
1137   }
1138   bgH->Multiply(dumH);
1139   delete dumH;
1140   return bgH;
1141 }
1142
1143 //______________________________________________________________________
1144 TObject* FindObject(int bin, const char* nameH, const TList* lst, Bool_t normPerEvent)
1145 {
1146   // get histo, optionally normalizing it per processed event
1147   if (!lst) {printf("FindObject %s: No list provided\n",nameH); exit(1);}
1148   int nent = lst->GetEntries();
1149   //
1150   char buff[200];
1151   if (bin>=0) sprintf(buff,"b%d_%s",bin,nameH);
1152   else  sprintf(buff,"%s",nameH);
1153   TString nm;
1154   TObject *hst = 0;
1155   for (int i=nent;i--;) {
1156     nm = lst->At(i)->GetName();
1157     if (nm.EndsWith(buff)) {hst = lst->At(i); break;}
1158   }
1159   if (!hst) {printf("FindObject: No bin %d %s histo in list %s\n",bin,nameH,lst->GetName()); exit(1);}
1160   if (!normPerEvent || hst->TestBit(kBitNormPerEvent)) return hst; // already normalized
1161   TString nameHS = nameH;
1162   if (nameHS==kHStatName) return hst;                              // never scale stat. histo
1163   //
1164   TH1* hstat = (TH1*)FindObject(-1,kHStatName,lst,kFALSE);
1165   double nrm = hstat->GetBinContent(kBinEntries + kEvProcInj+bin*kEntriesPerBin);
1166   //  double nrm = hstat->GetBinContent(kBinEntries + kEvProcData+bin*kEntriesPerBin); // HACK
1167   if (nrm<1) {printf("FindObject %s: Anomaluous %d number of events processed in bin %d of list %p\n",
1168                      buff,int(nrm),bin,lst); return 0;}
1169   //
1170   if      (hst->InheritsFrom(TH1::Class())) ((TH1*)hst)->Scale(1./nrm);
1171   else if (hst->InheritsFrom(THnSparse::Class())) {
1172     THnSparse* spr = (THnSparse*) hst;
1173     spr->Sumw2();
1174     int coord[3] = {0,0,0};
1175     for (Long64_t i = 0; i < spr->GetNbins(); ++i) {
1176       // Get the content of the bin from the current histogram
1177       Double_t v = spr->GetBinContent(i, coord);
1178       spr->SetBinContent(coord, v/nrm);
1179       spr->SetBinError(coord,TMath::Sqrt(v)/nrm);
1180     }
1181   }
1182   //
1183   hst->SetBit(kBitNormPerEvent);
1184   return hst;
1185 }
1186
1187 //______________________________________________________________________
1188 TList* LoadList(const char* flName, const char* addPref, const char* nameL)
1189 {
1190   // load list with histos
1191   TString nms = flName;
1192   gSystem->ExpandPathName(nms);
1193   TFile* fl = TFile::Open(nms.Data());
1194   if (!fl) {printf("LoadList: No file %s\n",nms.Data()); exit(1);}
1195   TList* lst = (TList*)fl->Get(nameL);
1196   if (!lst) {printf("LoadList: No list %s in file %s\n",nameL,nms.Data()); exit(1);}
1197   lst->SetName(flName);
1198   //
1199   int nEnt = lst->GetSize();
1200   TString nm;
1201   for (int i=0;i<nEnt;i++) {
1202     TNamed* ob = (TNamed*)lst->At(i);
1203     nm = addPref;
1204     nm += ob->GetName();
1205     ob->SetName(nm.Data());
1206   }
1207   //
1208   return lst;
1209 }
1210
1211 //____________________________________________________________________________
1212 void GetRatE(double x,double xe, double y,double ye, double &rat, double &rate)
1213 {
1214   rat = 0; rate = 0;
1215   if (TMath::Abs(y)<1e-16 || TMath::Abs(x)<1e-16) return;
1216   rat = x/y;
1217   rate = rat*TMath::Sqrt( xe*xe/(x*x) + ye*ye/(y*y));
1218 }
1219
1220 //____________________________________________________________________________
1221 void Integrate(TH1* hist, double xmn,double xmx, double &val, double& err)
1222 {
1223   // integrate 1d histo within given limits
1224   TAxis* xax = hist->GetXaxis();
1225   int bmn = xax->FindBin(xmn+kEps); if (bmn<1) bmn = 0; // include
1226   int bmx = xax->FindBin(xmx-kEps);
1227   val = hist->IntegralAndError(bmn,bmx,err);
1228   // is this histo with symmetric axis ? then integrate also negative half axis
1229   if (TMath::Abs( xax->GetXmin() + xax->GetXmax() )<1e-6) {
1230     bmn = xax->FindBin(-xmx+kEps);
1231     bmx = xax->FindBin(-xmn-kEps);
1232     double errn;
1233     val += hist->IntegralAndError(bmn,bmx,errn);
1234     err = TMath::Sqrt(err*err + errn*errn);
1235   }
1236 }
1237
1238
1239 //____________________________________________________________________________
1240 const char* HName(const char* prefix,const char* htype)
1241 {
1242   // compose the name of histo in the clist
1243   static TString strh;
1244   strh = "Tr"; strh += prefix; strh += "_"; strh += htype;
1245   return strh.Data();
1246 }
1247
1248 //____________________________________________________________________________
1249 Int_t CheckStat(const TList* lst, const char* dtType)
1250 {
1251   // check if needed bg was generated
1252   TH1* hstat = (TH1*)FindObject(-1,kHStatName,lst);
1253   TString dts = dtType;
1254   if (dts=="Data") return int( hstat->GetBinContent(kEvProcData) );
1255   if (dts=="Mix")  return int( hstat->GetBinContent(kEvProcMix) );
1256   if (dts=="Inj")  return int( hstat->GetBinContent(kEvProcInj) );
1257   if (dts=="Rot")  return int( hstat->GetBinContent(kEvProcRot) );
1258   printf("Unknown process %s statistics is checked. Alowed: Data,Mix,Inj,Rot",dtType);
1259   return 0;
1260 }
1261
1262
1263 void GetRealMinMax(TH1* histo, double &vmn, double &vmx)
1264 {
1265   TAxis *xax = histo->GetXaxis();
1266   int nbx = xax->GetNbins();
1267   vmn=1e6, vmx=-1e6;
1268   if (histo->InheritsFrom(TH2::Class())) {
1269     TAxis *yax = histo->GetYaxis();
1270     int nby = yax->GetNbins();
1271     for (int ix=nbx+2;ix--;) {
1272       for (int iy=nby+2;iy--;) {
1273         double vl = histo->GetBinContent(ix,iy);
1274         if (vl<kEps) continue;
1275         if (vl<vmn) vmn = vl;
1276         if (vl>vmx) vmx = vl;
1277       }
1278     }
1279   }
1280   //
1281   else {
1282     for (int ix=nbx+2;ix--;) {
1283       double vl = histo->GetBinContent(ix);
1284       if (vl<vmn) vmn = vl;
1285       if (vl>vmx) vmx = vl;
1286     }
1287   }
1288   //
1289 }