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