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