AliPhysics  9538fdd (9538fdd)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CorrectSpectraMultiMCBG.C
Go to the documentation of this file.
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 "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 = "corrFig";//"figMult";
48 const char* resDir = "corrRes";//"resMult";
49 Bool_t use1mBeta = true;
50 TString useBgType = "Comb"; // "Inj"; // "Inj"; "Comb";
51 Int_t useShapeType = kNormShapeDist; // which distribution to use for bg normalization
52 // If set 1 one, then
53 Int_t useMCLB = 2; // use Comb MC Labels as a template for Bg.
54 Double_t scaleBG = 1.3; // apply scaling to extracted bg for data
55 Int_t useScaleType = kSclIntegral;//kSclWghMean; // which type of tails normalization to use
56 Bool_t useZbinWAv = kFALSE; // kFALSE;
57 //Bool_t normToMB = kFALSE; // normalize with MB MC truth, rather than for selected events
58 Bool_t normToMB = kTRUE; // normalize with MB MC truth, rather than for selected events
59 const double kEtaFitRange = 0.5;
60 
61 Double_t minAlpha = 0.;
62 Double_t maxAlpha = 2.5;
63 
64 enum {kBitNormPerEvent=1<<14};
65 // bins for saved parameters in the hStat histo
66 enum {kDummyBin,
67  kEvTot0, // events read
68  kEvTot, // events read after vertex quality selection
69  kOneUnit, // just 1 to track primate merges
70  kNWorkers, // n workers
71  //
72  kCentVar, // cetrality var. used
73  kDPhi, // dphi window
74  kDTht, // dtheta window
75  kNStd, // N.standard deviations to keep
76  kPhiShift, // bending shift
77  kThtS2, // is dtheta scaled by 1/sin^2
78  kThtCW, // on top of w.dist cut cut also on 1 sigma dThetaX
79  kPhiOvl, // overlap params
80  kZEtaOvl, // overlap params
81  kNoOvl, // flag that overlap are suppressed
82  //
83  kPhiRot, // rotation phi
84  kInjScl, // injection scaling
85  kEtaMin, // eta cut
86  kEtaMax, // eta cut
87  kZVMin, // min ZVertex to process
88  kZVMax, // max ZVertex to process
89  //
90  kDPiSCut, // cut on dphi used to extract signal (when WDist is used in analysis, put it equal to kDPhi
91  kNStdCut, // cut on weighted distance (~1) used to extract signal
92  //
93  kMCV0Scale, // scaling value for V0 in MC
94  //
95  // here we put entries for each mult.bin
96  kBinEntries = 50,
97  kEvInMltBin = 0,
98  kEvProcData, // events with data mult.object (ESD or reco) passing all selections
99  kEvProcInj, // events Injected, total
100  kEvProcRot, // events Rotated
101  kEvProcMix, // events Mixed
102  kEntriesPerBin = 10
103 };
104 
105 
106 enum {kSigCorr,kMCPrim,kRawDtCut,kSignalEst,kSignalEstMC,kBgEst,k1MBeta,k1MBetaMC,k1MBetaMCscl,kAlpha,kAlphaMC,kBgMC,kBgRescFc,
107  kDataDist,kBgDist,kBgMCDist,kZvDist,kZvMCDistNS,kZvDistCorr,kZvEff,kMCShift=20, kNHistos=kMCShift+kMCShift};
108 
109 #ifndef __CINT__
110 #include <TROOT.h>
111 UShort_t fgDebug = 1;
112 void _MyPrint(UShort_t lvl, const char* fmt, ...)
113 {
114  if (lvl > fgDebug) return;
115  char buf[512];
116  va_list ap;
117  va_start(ap, fmt);
118  vsnprintf(buf, 511, fmt, ap);
119  buf[511] = '\0';
120  va_end(ap);
121  gROOT->IndentLevel();
122  Printf("%s",buf);
123 }
124 void Incr() { gROOT->IncreaseDirLevel(); }
125 void Decr() { gROOT->DecreaseDirLevel(); }
126 struct _MyGuard
127 {
128  _MyGuard(UShort_t lvl, const char* fmt, ...)
129  : ok(lvl <= fgDebug)
130  {
131  if (!ok) return;
132  char buf[512];
133  va_list ap;
134  va_start(ap, fmt);
135  vsnprintf(buf, 511, fmt, ap);
136  buf[511] = '\0';
137  va_end(ap);
138  gROOT->IndentLevel();
139  Printf("%s",buf);
140  Incr();
141  }
142  ~_MyGuard() { if (ok) Decr(); }
143  Bool_t ok;
144 };
145 #define MyPrint(L,F,...) _MyPrint(L,F, ## __VA_ARGS__)
146 #define MyGuard(L,F,...) _MyGuard _guard(L,F, ## __VA_ARGS__)
147 #endif
148 
149 void PrintAndPause(TCanvas* c, const TString& what, Bool_t wait)
150 {
151  Printf("Wat is %d", wait);
152  c->Modified();
153  c->Update();
154  c->cd();
155  c->Print(what);
156  if (wait) c->WaitPrimitive();
157 }
158 
159 void CorrectSpectraMultiMCBG(const char* flNameData, const char* flNameMC,const char* unique="",int maxBins=10, Bool_t waitForUser=false, const char* bgType="Comb");
160 Bool_t PrepareHistos(int bin, TList* lst, TList* lisMC, Bool_t isMC);
161 void ProcessHistos(int bin);
162 TH1* NormalizeBg(TH1* dataH, TH1* bgH, double &scl, double &scle);
163 TObject* FindObject(int bin, const char* nameH, const TList* lst, Bool_t normPerEvent=kTRUE);
164 TList* LoadList(const char* flName, const char* addPref, const char* nameL="clist");
165 void GetRatE(double x,double xe, double y,double ye, double &rat, double &rate);
166 Int_t CheckStat(const TList* lst,const char* dtType);
167 void Integrate(TH1* hist, double xmn,double xmx, double &val, double& err);
168 void CropHisto(TH1* histo, int b00, int b01, int b10=-1, int b11=-1);
169 void CropHisto(TH1* histo, double b00, double b01, double b10=-1, double b11=-1);
170 void GetRealMinMax(TH1* h, double &vmn, double &vmx);
171 const char* HName(const char* prefix,const char* htype);
172 TH1* ProjNorm(TH2* hEtaZ, TH1* hZv,const char* name="_px", Int_t firstbin=0, Int_t lastbin=-1);
173 TH1* ProjectWghMean(TH2* hEtaZ, const char* name = "_px", Int_t firstbin = 0, Int_t lastbin = -1, double rejOutliers=6.);
174 void CorrectForZV(TH2* hEtaZ, TH1* hZv);
175 void KillBadBins(TH2* histo, double mn=-1e50,double mx=1e50);
176 
177 void PlotResults(Bool_t waitForUser);
178 void PlotDNDEta(int bin);
179 void PlotAlphaBeta(int bin);
180 void PlotSpecies();
181 void PrintH(TH2* h, Int_t prec=2);
182 void PrintH(TH1* h, Int_t prec=2);
183 
184 Int_t nCentBins = -1;
185 TList *listDt=0, *listMC=0;
186 TObjArray resArr, resDnDeta;
187 char outStr[1000];
188 char outTitle[1000];
189 TString uniqueName="";
190 //
191 TArrayD dNdEta,dNdEtaErr;
192 TCanvas *canvFin=0;
193 //
194 Bool_t creatDnDEtaCMacro = kTRUE;
195 Bool_t creatAlphaBetaCMacro = kTRUE;
196 Bool_t creatSpeciesCMacro = kTRUE;
197 
198 void CorrectSpectraMultiMCBG(const char* flNameData, const char* flNameMC, const char* uniqueNm,int maxBin,Bool_t waitForUser,const char* bgType)
199 {
200  useBgType = bgType;
201  //
202  uniqueName = uniqueNm;
203  Printf("WaitForUser is %d", waitForUser);
204  MyPrint(1,"Loading lists from %s and %s", flNameData, flNameMC);
205  {
206  MyGuard(1,"Loading data lists");
207  listDt = LoadList(flNameData,"dt_");
208  }
209  {
210  MyGuard(1,"loading MC lists");
211  listMC = LoadList(flNameMC,"mc_");
212  }
213  //
214  resArr.Clear();
215  //
216  TH1* hstat = (TH1*)FindObject(-1,kHStatName,listDt,kFALSE);
217  MyPrint(2,"Got statistics histogram %p", hstat);
218 
219  //TH1* hstat = (TH1*)FindObject(-1,kHStatName,listMC,kFALSE);
220  //
221  int nbstat = hstat->GetNbinsX();
222  nCentBins = (nbstat - kBinEntries)/kEntriesPerBin;
223  nCentBins = nCentBins>maxBin ? maxBin : nCentBins;
224  MyPrint(2,"Will process %d centrality bins", nCentBins);
225 
226  if (nCentBins<1) return;
227  double myMergeFactor = hstat->GetBinContent(kOneUnit);
228  if (myMergeFactor<1) myMergeFactor = 1.;
229  MyPrint(2,"Output was merged from %5.1f jobs", myMergeFactor);
230 
231  //
232  dNdEta.Set(nCentBins);
233  dNdEtaErr.Set(nCentBins);
234  //
235  kdPhiSgCut = hstat->GetBinContent(kDPiSCut)/myMergeFactor;
236  kWDistSgCut = hstat->GetBinContent(kNStdCut)/myMergeFactor;
237  MyPrint(2,"Signal cuts used: dPhiS: %f WDist:%f",
238  kdPhiSgCut,kWDistSgCut);
239 
240  {
241  MyGuard(1,"Now looping over bins");
242  for (int ib=0;ib<nCentBins;ib++) {
243  {
244  MyGuard(1,"Extracting from MC");
245  if (!PrepareHistos(ib,listMC,listMC,kTRUE)) continue;//return;
246  }
247  {
248  MyGuard(1,"Extracting from Data");
249  if (!PrepareHistos(ib,listDt,listMC,kFALSE)) continue;//return;
250  }
251  {
252  MyGuard(1,"Processing information");
253  ProcessHistos(ib);
254  }
255  }
256  }
257 
258  //
259  sprintf(outStr,"CutEta%.1f_%.1f_Zv%.1f_%.1f_bg%s_Shape_%s_%s_cutSig%.1f_cutBg%.1f",
260  hstat->GetBinContent(kEtaMin)/myMergeFactor,
261  hstat->GetBinContent(kEtaMax)/myMergeFactor,
262  hstat->GetBinContent(kZVMin)/myMergeFactor,
263  hstat->GetBinContent(kZVMax)/myMergeFactor,
264  useBgType.Data(),
265  useShapeType==kNormShapeDist ? "wdst":"dphi",
266  (use1mBeta ? "MCBG" : "DTBG"),
267  useShapeType==kNormShapeDist ? kWDistSgCut:kdPhiSgCut,
268  useShapeType==kNormShapeDist ? kWDistBgTailMin:kdPhiBgTailMin);
269  //
270  PlotResults(waitForUser);
271  //
272  printf("Final Results:\n");
273  printf("dNdEta: ");
274  for (int i=nCentBins;i--;) printf("%.2f,",dNdEta[i]); printf("\n");
275  printf("dNdEtaErr: ");
276  for (int i=nCentBins;i--;) printf("%.2f,",dNdEtaErr[i]); printf("\n");
277  //
278  TString rtnm1 = resDir; rtnm1 += "/"; rtnm1 += uniqueName;
279  rtnm1 += "_"; rtnm1 += nCentBins; rtnm1+= "bins_";
280  rtnm1 += outStr; rtnm1 += ".root";
281  TFile* flRes = TFile::Open(rtnm1.Data(),"recreate");
282  flRes->WriteObject(&resDnDeta, "TObjArray", "kSingleKey");
283  flRes->WriteObject(&resArr, "TObjArrayAux", "kSingleKey");
284  flRes->Close();
285  delete flRes;
286  printf("Stored result in %s\n",rtnm1.Data());
287  //
288 }
289 
290 //_____________________________________________________________________
291 TH1* CopyAdd(TH1* h,
292  const char* name,
293  const char* title,
294  TObjArray* col,
295  Int_t location,
296  Int_t shift)
297 {
298  TH1* ret = static_cast<TH1*>(h->Clone(name));
299  ret->SetTitle(title);
300  if (col) col->AddAtAndExpand(ret, location+shift);
301  return ret;
302 }
303 
304 //_____________________________________________________________________
305 Bool_t PrepareHistos(int bin, TList* lst, TList* lstMC, Bool_t isMC)
306 {
307  if (!lst) {
308  Warning("PrepareHistos", "No list for %s", isMC ? "simulation" : "data");
309  return false;
310  }
311  // fill standard histos for given bin
312  //
313  char buffn[500];
314  char bufft[500];
315  //
316  double cutBgMin,cutBgMax;
317  double cutSgMin,cutSgMax;
318  //
319  if (useShapeType==kNormShapeDist) {
320  cutBgMin = kWDistBgTailMin;
321  cutBgMax = kWDistBgTailMax;
322  //
323  cutSgMin = 0;
324  cutSgMax = kWDistSgCut;
325  }
326  else {
327  cutBgMin = kdPhiBgTailMin;
328  cutBgMax = kdPhiBgTailMax;
329  //
330  cutSgMin = 0;
331  cutSgMax = kdPhiSgCut;
332  }
333  //
334  TH1* hstat = (TH1*)FindObject(-1,kHStatName,lst,kFALSE);
335  MyPrint(2,"Got statistics histogram %p", hstat);
336 
337  double nrmBin = hstat->GetBinContent(kBinEntries + kEvProcData+bin*kEntriesPerBin);
338  if (nrmBin<1) {
339  Warning("PrepareHistos", "Norm is 0 for bin %d",bin);
340  return kFALSE;
341  }
342  MyPrint(0,"Normalization for bin %d: %f", bin, nrmBin);
343 
344  const char* zeCut = "ZvEtaCutT";
345  TObjArray* res = &resArr;
346  int shift = bin*kNHistos + (isMC ? kMCShift : 0);
347  //
348  // histo for "data" Z vs Eta with default signal cut - scaled by IPz dist
349  TH2* tmp2 = (TH2*)FindObject(bin,HName("Data",zeCut),lst);
350  if (!tmp2) {
351  Warning("PrepareHistos", "Didn't find %s", HName("Data",zeCut));
352  return kFALSE;
353  }
354  TH2* hRawDtCut = (TH2*)CopyAdd(tmp2,
355  Form("bin%d_%s_RawWithCut",bin,isMC?"mc":"dt"),
356  Form("bin%d %s Raw Data with cut on tracklets",
357  bin,isMC ? "mc":"dt"),
358  res,kRawDtCut,shift);
359  MyPrint(1,"Measurements %s (clone of %s)", hRawDtCut->GetName(),
360  tmp2->GetName());
361 
362  // Number of eta, IPz bins
363  int nbEta = hRawDtCut->GetXaxis()->GetNbins();
364  int nbZV = hRawDtCut->GetYaxis()->GetNbins();
365  //
366  // Zv nonuniformity histo - IPz vs Centrality
367  TH2* hzv2 = (TH2*) FindObject(-1,"zv",lst, kFALSE);
368  TH1* hz = hzv2->ProjectionX(Form("zv%d_%s",bin,isMC ? "mc":"dt"),
369  bin+1,bin+1,"e");
370  MyPrint(1,"Get histogram of vertex disp for bin %d and scale by %f",
371  bin, nrmBin);
372  hz->Scale(1./nrmBin);
373  res->AddAtAndExpand(hz, kZvDist+shift);
374 
375  //
376  // "Data - Est.Bg" histo with cut on tails where we look for signal
377  TH2* hSignalEst =
378  (TH2*)CopyAdd(hRawDtCut,
379  Form("bin%d_%s_SignalWithCut",bin,isMC ? "mc":"dt"),
380  Form("bin%d %s Signal (raw-bg) with cut on tracklets",
381  bin,isMC ? "mc":"dt"),
382  res,kSignalEst,shift);
383  MyPrint(1,"Signal estimate: %s (clone of %s)", hSignalEst->GetName(),
384  hRawDtCut->GetName());
385 
386  //
387  // "Data - MC.Bg" histo with cut on tails where we look for signal
388  TH2* hSignalEstMC = 0;
389  if (isMC) {
390  MyGuard(1, "MC specific code - MC signal estimate");
391  hSignalEstMC =
392  (TH2*)CopyAdd(hRawDtCut,
393  Form("bin%d_%s_SignalWithCut_bgMCLabels",
394  bin,isMC ? "mc":"dt"),
395  Form("bin%d %s Signal (raw-bg_MCLabels) w/cut on tracklets",
396  bin,isMC ? "mc":"dt"),
397  res,kSignalEstMC,shift);
398  MyPrint(1,"MC signal estimate %s (copy of %s)",hSignalEstMC->GetName(),
399  hRawDtCut->GetName());
400  // also currently a copy of eta vs IPz
401  }
402  //
403 
404  // Estimated background in the cut range
405  TList* localLst = (useBgType.EqualTo("Comb") ? lstMC : lst);
406  tmp2 = (TH2*)FindObject(bin,HName(useBgType.Data(),zeCut),localLst);
407  TH2* hBgEst = (TH2*)CopyAdd(tmp2,
408  Form("bin%d_%s_BgEst",bin,isMC ? "mc":"dt"),
409  Form("bin%d %s Estimated Bg",bin,isMC?"mc":"dt"),
410  res,kBgEst,shift);
411  MyPrint(1,"Bg estimate %s (clone of %s)",hBgEst->GetName(),tmp2->GetName());
412 
413  TH2* h1mBeta = (TH2*)CopyAdd(hBgEst,
414  Form("bin%d_%s_1mBeta",bin,isMC ? "mc":"dt"),
415  Form("bin%d %s 1-#beta with estimated bg",
416  bin,isMC ? "mc":"dt"),
417  res, k1MBeta, shift);
418  h1mBeta->Reset();
419  MyPrint(1, "1-beta %s (clone of %s)", h1mBeta->GetName(), hBgEst->GetName());
420 
421  //
422  // If MC labels info is provided, prepare 1-beta with MC bg
423  TH2* h1mBetaMC = 0; // 1-beta for MC with bg from labels
424  TH2* h1mBetaMCscl = 0; // scaled by ad hoc factor
425  TH2* hBgMC = 0; // bg from MC labels
426  if (isMC) {
427  MyGuard(1, "MC specific code - prepare 1-beta");
428  tmp2 = (TH2*)FindObject(bin,HName("Comb",zeCut),lstMC);
429  MyPrint(1,"Get MC background from %s", HName("Comb",zeCut));
430  if (!tmp2) return kFALSE;
431  hBgMC = (TH2*)CopyAdd(tmp2,
432  Form("bin%d_%s_BgMC",bin,isMC ? "mc":"dt"),
433  Form("bin%d %s Bg from MC labels",
434  bin,isMC ? "mc":"dt"),
435  res,kBgMC,shift);
436  MyPrint(1,"Set MC background %s (clone of %s)",hBgMC->GetName(),
437  tmp2->GetName());
438 
439  // printf(">>Here1\n");
440  // Now scale by measured eta vs IPz
441  h1mBetaMC = (TH2*)CopyAdd(hBgMC,
442  Form("bin%d_%s_h1mBetaMC",bin,isMC ? "mc":"dt"),
443  Form("bin%d %s 1-#beta with bg. from MC labels",
444  bin,isMC ? "mc":"dt"),
445  res, k1MBetaMC,shift);
446  h1mBetaMC->Divide(hRawDtCut);
447  MyPrint(1, "1-beta (MC) %s (ratio of %s to %s)",
448  h1mBetaMC->GetName(),
449  hBgMC->GetName(), hRawDtCut->GetName());
450 
451  h1mBetaMCscl = (TH2*)CopyAdd(h1mBetaMC,
452  Form("%s_scl",h1mBetaMC->GetName()),
453  h1mBetaMC->GetTitle(),
454  res,k1MBetaMCscl,shift);
455  h1mBetaMCscl->Scale(scaleBG);
456  MyPrint(1, "1-beta (MC) %s (clone of %s)", h1mBetaMCscl->GetName(),
457  h1mBetaMC->GetName());
458  // printf("<<Here1\n");
459  for (int ib0=1;ib0<=nbEta;ib0++)
460  for (int ib1=1;ib1<=nbZV;ib1++) {
461  // Do the "1-" part
462  Double_t oneMBeta = h1mBetaMC->GetBinContent(ib0,ib1);
463  Double_t oneMBetaScl = h1mBetaMCscl->GetBinContent(ib0,ib1);
464  h1mBetaMCscl->SetBinContent(ib0,ib1, 1.- oneMBetaScl);
465  h1mBetaMC->SetBinContent(ib0,ib1, 1.- oneMBeta);
466  }
467  MyPrint(1,"Did 1-X on %s and %s (%f)", h1mBetaMC->GetName(),
468  h1mBetaMCscl->GetName(), scaleBG);
469  //
470  // Now we have (1-beta)*meassured=no_combinatorics_signal
471  hSignalEstMC->Multiply(h1mBetaMC);
472  MyPrint(1,"Multiply %s onto %s", h1mBetaMC->GetName(),
473  hSignalEstMC->GetName());
474  }
475  //
476  // uncut w.distance or dphi distribution for data
477  TH1* tmp1 = (TH1*)FindObject(bin,HName("Data",useShapeType==kNormShapeDist?
478  "WDist":"DPhiS"),lst);
479  if (!tmp1) return kFALSE;
480  TH1* hDstDt = CopyAdd(tmp1,
481  Form("bin%d_%s_DistRawData",bin,isMC ? "mc":"dt"),
482  Form("bin%d %s Raw Distance for Data",
483  bin,isMC?"mc":"dt"),
484  res,kDataDist,shift);
485  MyPrint(1,"Delta dist %s (clone of %s)", hDstDt->GetName(), tmp1->GetName());
486 
487  // Integrate tail of measured Delta dist
488  double nrmDst,dumErr = 0;
489  Integrate(hDstDt, cutBgMin,cutBgMax, nrmDst, dumErr);
490  MyPrint(1,"Integral of %s is %f +/- %f", hDstDt->GetName(), nrmDst, dumErr);
491 
492  if (nrmDst<1e-10) nrmDst = 1.; // RS resque
493  // Scale measured Delta distribution by integral
494  hDstDt->Scale(1./nrmDst);
495  MyPrint(1,"Scaled %s and by integral %f", hDstDt->GetName(),nrmDst);
496 
497  //
498  // uncut w.distance or dphi distribution for generated bg
499  TH1* hDstBg = 0;
500  if (!useBgType.IsNull()) {
501  tmp1 = (TH1*)FindObject(bin,HName(useBgType.Data(),
502  useShapeType==kNormShapeDist ?
503  "WDist":"DPhiS"),localLst);
504  if (!tmp1) return false;
505  hDstBg = CopyAdd(tmp1,
506  Form("bin%d_%s_DistRawGenBgNorm",bin,isMC ? "mc":"dt"),
507  Form("bin%d %s Raw Dist. for Gen.Bg. Normalized to data",
508  bin,isMC ? "mc":"dt"),
509  res, kBgDist, shift);
510  // Scale bg Delta to mesaured Delta integral
511  hDstBg->Scale(1./nrmDst);
512  MyPrint(1,"%s (copy of %s) scaled by %f", hDstBg->GetName(),
513  tmp1->GetName(), nrmDst);
514  }
515  //
516  // fill 1-beta matrix
517  double scl,sclE;
518  // get rescaling factor for bg. from tails comparison
519  hDstBg = NormalizeBg(hDstDt,hDstBg,scl,sclE);
520  MyPrint(1, "Scaling between %s and %s: %f +/- %f",
521  hDstDt->GetName(), hDstBg->GetName(), scl, sclE);
522  res->AddAtAndExpand(hDstBg, kBgDist+shift);
523  double bgVal,bgErr;
524  double dtVal,dtErr;
525  // integral in the range where we look for signal
526  Integrate(hDstBg, cutSgMin, cutSgMax, bgVal, bgErr);
527  Integrate(hDstDt, cutSgMin, cutSgMax, dtVal, dtErr);
528 
529  double sclb,sclbErr;
530  GetRatE(bgVal,bgErr, dtVal, dtErr,sclb,sclbErr);
531  MyPrint(1,"Ratio of (%s) %f +/- %f over (%s) %f +/- %f -> %f +/- %f",
532  hDstBg->GetName(), bgVal, bgErr,
533  hDstDt->GetName(), dtVal, dtErr,
534  sclb, sclbErr);
535 
536  // hDstBg->Scale(1./nrmDst);
537  //
538  // finalize estimated bg and signal matrices
539  hBgEst->Scale(scl);
540  MyPrint(1, "Scaling %s by %f", hBgEst->GetName(), scl);
541  //
542  if (use1mBeta) {
543  MyGuard(1,"Using MC Combinatorials fraction for real data bg.estimate");
544  TH2* bet1m = (TH2*)res->At((isMC ? k1MBetaMC : k1MBetaMCscl)
545  + shift + (isMC ? 0 :kMCShift) );
546  MyPrint(1,"Multiply %s with %s", hSignalEst->GetName(), bet1m->GetName());
547  hSignalEst->Multiply(bet1m);
548  // if (isMC) PrintH(hSignalEst);
549  }
550  else {
551  MyPrint(1,"Subtracting %s from %s",hBgEst->GetName(),hSignalEst->GetName());
552  hSignalEst->Add(hBgEst,-1);
553  }
554  //
555  // finalize 1-beta
556  {
557  MyGuard(1,"Finalize 1-beta");
558  for (int ib0=1;ib0<=nbEta;ib0++) { // eta
559  for (int ib1=1;ib1<=nbZV;ib1++) { // zv
560  // printf("Bin %d %d\n",ib0,ib1);
561  double bg = hBgEst->GetBinContent(ib0,ib1);
562  double bgE = hBgEst->GetBinError(ib0,ib1);
563  double dt = hRawDtCut->GetBinContent(ib0,ib1);
564  double dtE = hRawDtCut->GetBinError(ib0,ib1);
565  double beta,betaE;
566  GetRatE(bg,bgE,dt,dtE, beta,betaE );
567  h1mBeta->SetBinContent(ib0,ib1,1.-beta);
568  h1mBeta->SetBinError(ib0,ib1,betaE);
569  //
570  }
571  }
572  MyPrint(1, "%s calculated as 1 minus ratio of %s to %s",
573  h1mBeta->GetName(), hBgEst->GetName(), hRawDtCut->GetName());
574  }
575  //
576  if (isMC) {
577  MyGuard(1,"MC specific stuff");
578  // Zv nonuniformity histo: MCGen before cuts
579  TH2* hzv2ns = (TH2*) FindObject(-1,"zvMCNoPS",lst, kFALSE);
580  MyPrint(1, "MC truth vertex distribution %s", hzv2ns->GetName());
581  TH1* hzns = hzv2ns->ProjectionX(Form("zvMCNS%d",bin),bin+1,bin+1,"e");
582  MyPrint(1, "Projection %s", hzns->GetName());
583  //
584  double myMergeFactor = hstat->GetBinContent(kOneUnit);
585  double zmn = hstat->GetBinContent(kZVMin)/myMergeFactor;
586  double zmx = hstat->GetBinContent(kZVMax)/myMergeFactor;
587  int zbmn = hzns->FindBin(zmn+1e-3);
588  int zbmx = hzns->FindBin(zmx-1e-3);
589  double nEvTot = hzns->Integral(zbmn,zbmx);
590  MyPrint(1,"Scale MC Truth By %9.1f",nEvTot);
591  if (nEvTot<1) exit(1);
592  hzns->Scale(1./nEvTot);
593  /*
594  // here there could be some over/under flows, add them to
595  // first/last used bins correspondingly
596  //
597  double under = hzns->Integral(0,zbmn-1);
598  double over = hzns->Integral(zbmx+1,hzns->GetNbinsX()+1);
599  hzns->SetBinContent(zbmn, hzns->GetBinContent(zbmn)+under);
600  hzns->SetBinContent(zbmx, hzns->GetBinContent(zbmx)+over);
601  printf("Add %e and %e of under/overflows to extreme bins\n",under,over);
602  */
603  res->AddAtAndExpand(hzns, kZvMCDistNS+shift);
604  //
605  TH2* mcPrim = 0;
606  // prepare MC primary signal histo
607  if (normToMB) {
608  // don't normalise to NSel
609  tmp2 = (TH2*)FindObject(bin,"zvEtaPrimMC", lst, kFALSE );
610  mcPrim = (TH2*)CopyAdd(tmp2,
611  Form("bin%d_zvEtaPrimMC",bin),
612  Form("bin%d MC True signal Zv vs Eta",bin),
613  res, kMCPrim, shift);
614  mcPrim->Scale(1./nEvTot);
615  MyPrint(1,"%s (copy of %s) scaled by %f", mcPrim->GetName(),
616  tmp2->GetName(), nEvTot);
617  }
618  else {
619  // don't normalise to NSel
620  tmp2 = (TH2*)FindObject(bin,"zvrecEtaPrimMCSel", lst, kTRUE );
621  mcPrim = (TH2*)CopyAdd(tmp2,
622  Form("bin%d_zvEtaPrimMCSel", bin),
623  Form("bin%d MC True signal Zv vs Eta (sel)",bin),
624  res, kMCPrim, shift);
625  MyPrint(1,"%s (copy of %s)", mcPrim->GetName(), tmp2->GetName());
626  }
627  }
628  //
629  return kTRUE;
630 }
631 
632 //_____________________________________________________________________
633 void ProcessHistos(int bin)
634 {
635  //
636  int shift = bin*kNHistos;
637  MyPrint(2,"Shift for bin %d is %d", bin, shift);
638  //
639  TString prefN = "bin"; prefN += bin; prefN+="_";
640  TString prefT = "bin"; prefT += bin; prefT+=" ";
641  TObjArray* res = &resArr;
642  //
643  // // build MC truth
644  // TH2* mctru2 = (TH2*)res->At(shift + kMCShift + kMCPrim);
645  // mctru2 = (TH2*) mctru2->Clone(prefN+"MCtruth2D");
646  // mctru2->SetTitle(prefN+"MCtruth2D");
647  // TH1F* mctru = mctru2->ProjectionX(prefN+"MCtruth");
648  // mctru->SetTitle(prefN+"MCtruth");
649  // double binNorm = mctru2->GetXaxis()->GetBinWidth(1)*mctru2->GetYaxis()->GetBinWidth(1);
650  // mctru->Scal(1./binNorm);
651  // res->AddAtAndExpand(mctru, shift + kMCTruth);
652  //
653  TH1* hZV = (TH1*)res->At(shift + kZvDist);
654 
655  TH1* hZVmcRec = (TH1*)res->At(shift + kMCShift + kZvDist);
656  TH1* hZVmcGen = (TH1*)res->At(shift + kMCShift + kZvMCDistNS);
657  TH1* mcVtxEff = (TH1*)hZVmcRec->Clone(Form("Eff_%s",hZVmcRec->GetName()));
658  MyPrint(1,"Got vertex distributions %s (data) %s %s (MC)",
659  hZV->GetName(), hZVmcRec->GetName(), hZVmcGen->GetName());
660 
661  mcVtxEff->Divide(hZVmcRec,hZVmcGen,1,1,"B"); // vertex rec eff
662  MyPrint(1,"Calculate vertex efficiency as %s over %s -> %s",
663  hZVmcRec->GetName(), hZVmcGen->GetName(), mcVtxEff->GetName());
664 
665  TH1* hZVcorr = (TH1*)hZV->Clone(Form("corr_%s",hZV->GetName()));
666  hZVcorr->Divide(mcVtxEff);
667  MyPrint(1,"Calculate corrected vertex distribution %s", hZVcorr->GetName());
668  //
669  res->AddAt(hZVcorr, shift+kZvDistCorr);
670  res->AddAt(mcVtxEff, shift+kMCShift+kZvEff);
671  MyPrint(2,"Adding vtx corr (%s) and eff (%s) at %d and %d",
672  hZVcorr->GetName(), mcVtxEff->GetName(), shift+kZvDistCorr,
673  shift+kMCShift+kZvEff);
674  //
675  // build alpha matrix
676  TH2* halp = (TH2*)res->At(shift + kMCShift + kMCPrim);
677  MyPrint(1,"Start of alpha: %s", halp->GetName());
678  halp = (TH2*) halp->Clone(prefN+"Alpha");
679  MyPrint(1,"Cloned as %s", halp->GetName());
680  halp->SetTitle(prefN+"#alpha");
681 
682  TH2* omBeta = (TH2*)res->At(shift + kMCShift + k1MBeta);
683  TH2* mcRaw = (TH2*)res->At(shift + kMCShift + kRawDtCut);
684  halp->Divide( omBeta );
685  halp->Divide( mcRaw );
686  MyPrint(1,"%s scaled by %s and %s", halp->GetName(), omBeta->GetName(),
687  mcRaw->GetName());
688  res->AddAtAndExpand(halp, shift + kAlpha);
689  MyPrint(2,"Adding alpha %s at %d", halp->GetName(), shift+kAlpha);
690 
691  //
692  KillBadBins(halp,minAlpha,maxAlpha);
693  //
694  // build alpha matrix with MC labels bg
695  TH2* halpMC = (TH2*)res->At(shift + kMCShift + kMCPrim);
696  // PrintH(halpMC);
697  MyPrint(1,"Start of alpha MC: %s", halpMC->GetName());
698  halpMC = (TH2*) halpMC->Clone(prefN + "AlphaMC");
699  MyPrint(1,"Cloned as %s", halpMC->GetName());
700  halpMC->SetTitle(prefT + "#alpha MC labels");
701 
702  TH2* omBetaMC = (TH2*)res->At(shift + kMCShift + k1MBetaMC);
703  halpMC->Divide( omBetaMC );
704  halpMC->Divide( mcRaw );
705  MyPrint(1, "%s scaled by %s and %s", halpMC->GetName(), omBetaMC->GetName(),
706  mcRaw->GetName());
707  res->AddAtAndExpand(halpMC, shift + kAlphaMC);
708  //
709  KillBadBins(halpMC,minAlpha,maxAlpha);
710  // PrintH(halpMC);
711  //
712  // build corrected signal
713  TH2* hsigCorr = (TH2*)res->At(shift + kSignalEst);
714  MyPrint(1,"Start of corrected signal: %s", hsigCorr->GetName());
715  hsigCorr = (TH2*) hsigCorr->Clone(prefN + "SignalEstCorr");
716  MyPrint(1,"Cloned as %s", hsigCorr->GetName());
717  hsigCorr->SetTitle(prefT + "Corrected Signal");
718 
719  hsigCorr->Multiply(use1mBeta ? halpMC : halp );
720  // PrintH(hsigCorr);
721  MyPrint(1,"Multiply by %s", use1mBeta ? halpMC->GetName() : halp->GetName());
722  res->AddAtAndExpand(hsigCorr, shift + kSigCorr);
723  MyPrint(2,"Add corrected signal %s at %d", hsigCorr->GetName(), shift+kSigCorr);
724  //
725 
726  TH2* tmpSC = (TH2*)hsigCorr->Clone("tmpSC");
727  TH1* useZV = normToMB ? hZVcorr : hZV;
728  TH1* hsigCorrX = 0;
729  if (useZbinWAv) {
730  // CorrectForZV(tmpSC,hZV);
731  CorrectForZV(tmpSC,useZV);
732  hsigCorrX = ProjectWghMean(tmpSC,"DataCorrSignalX");
733  MyPrint(1,"Corrected X signal from weighted mean: %s",
734  hsigCorrX->GetName());
735  }
736  else {
737  // hsigCorrX = ProjNorm(tmpSC,hZV,"DataCorrSignalX");
738  hsigCorrX = ProjNorm(tmpSC,useZV,"DataCorrSignalX");
739  MyPrint(1,"Corrected X signal from normal mean %s", hsigCorrX->GetName());
740  }
741  //
742  delete tmpSC;
743  hsigCorrX->Scale(1./hsigCorrX->GetBinWidth(1));
744  MyPrint(1,"Scaling corrected X signal %s with bin width",
745  hsigCorrX->GetName());
746 
747  TF1* pl0 = new TF1("pl0","pol0");
748  pl0->SetParameter(0,hsigCorr->GetMinimum());
749  hsigCorrX->Fit(pl0,"q0","",-kEtaFitRange,kEtaFitRange);
750  MyPrint(1,"Fit straight line to X signal %s", hsigCorrX->GetName());
751  double fval = pl0->GetParameter(0);
752  double ferr = pl0->GetParError(0);
753  delete hsigCorrX;
754  dNdEta[bin] = fval;
755  dNdEtaErr[bin] = ferr;
756  printf("Bin %d: dN/d#eta_{|#eta|<0.5} = %.2f %.2f\n",bin, fval,ferr);
757  //
758 }
759 
760 void PlotResults(Bool_t waitForUser)
761 {
762  MyGuard(1,"Plottinmg results");
763  TString psnm1 = figDir; psnm1 += "/"; psnm1 += uniqueName;
764  psnm1 += "_"; psnm1 += nCentBins; psnm1+= "bins_";
765  psnm1 += outStr; psnm1 += ".pdf";
766  TString psnm0 = psnm1.Data();
767  psnm0 += "[";
768  TString psnm2 = psnm1.Data();
769  psnm2 += "]";
770  //
771  TH1* hstat = (TH1*)FindObject(-1,kHStatName,listDt,kFALSE);
772  double myMergeFactor = hstat->GetBinContent(kOneUnit);
773  if (myMergeFactor<1) myMergeFactor = 1.;
774  //
775  TH1* hbins = (TH1*)FindObject(-1,"EvCentr",listDt,kFALSE);
776  //
777  if (!canvFin) canvFin = new TCanvas("canvFin", "canvFin",0,50,700,1000);
778  canvFin->Clear();
779  //
780  PrintAndPause(canvFin,psnm0.Data(),false);
781  //
782  canvFin->Divide(1,2);
783  canvFin->cd(1);
784  gPad->SetLeftMargin(0.15);
785  TGraphErrors* grp = new TGraphErrors(nCentBins);
786  for (int i=0;i<nCentBins;i++) {
787  grp->SetPoint(i,hbins->GetBinCenter(i+1),dNdEta[i]);
788  grp->SetPointError(i,hbins->GetBinWidth(i+1)/2,dNdEtaErr[i]);
789  }
790  grp->SetMarkerStyle(20);
791  grp->SetMarkerColor(kRed);
792  grp->SetLineColor(kRed);
793  grp->SetMinimum(1e-6);
794  grp->Draw("ap");
795  grp->GetXaxis()->SetTitle("Centrality Variable");
796  grp->GetYaxis()->SetTitle("dN/d#eta_{|#eta|<0.5}");
797  grp->GetYaxis()->SetTitleOffset(1.8);
798  MyPrint(2,"Created graph of central dNch/deta %s", grp->GetName());
799  gPad->SetGrid(1,1);
800  //
801  canvFin->cd(2);
802  gPad->SetLeftMargin(0.15);
803  hbins->Draw();
804  hbins->SetMinimum(1e-6);
805  hbins->SetMarkerStyle(20);
806  hbins->SetMarkerColor(kRed);
807  hbins->SetLineColor(kRed);
808  hbins->GetYaxis()->SetTitle("accepted events");
809  hbins->GetYaxis()->SetTitleOffset(1.8);
810  gPad->SetGrid(1,1);
811  MyPrint(2,"Drew %s", hbins->GetName());
812  //
813  canvFin->cd(0);
814  PrintAndPause(canvFin,psnm1.Data(),waitForUser);
815  //
816  const TArrayD &binArr = *hbins->GetXaxis()->GetXbins();
817  //
818  {
819  MyGuard(1,"Looping over centrality bins");
820  for (int i=0;i<nCentBins;i++) {
821  //
822  sprintf(outTitle,"%s, %d<C_%s<%d, %.1f<#eta<%.1f, %.1f<Z_{V}<%.1f, Bg.:%s, %s, CutVar:%s, |sig|<%.2f, %.2f<|bg.nrm|<%.2f",
823  uniqueName.Data(),
824  (int)binArr[i],
825  hstat->GetXaxis()->GetBinLabel(kCentVar),
826  (int)binArr[i+1],
827  hstat->GetBinContent(kEtaMin)/myMergeFactor,
828  hstat->GetBinContent(kEtaMax)/myMergeFactor,
829  hstat->GetBinContent(kZVMin)/myMergeFactor,
830  hstat->GetBinContent(kZVMax)/myMergeFactor,
831  useBgType.Data(),
832  (use1mBeta ? "MC-BG" : "DT-BG"),
833  useShapeType==kNormShapeDist ? "#Delta":"#Delta#varphi-#delta_{#varphi}",
834  useShapeType==kNormShapeDist ? kWDistSgCut:kdPhiSgCut,
835  useShapeType==kNormShapeDist ? kWDistBgTailMin : kdPhiBgTailMin,
836  useShapeType==kNormShapeDist ? kWDistBgTailMax : kdPhiBgTailMax
837  );
838  //
839  MyPrint(2,"Title: %s", outTitle);
840  PlotDNDEta(i);
841  PrintAndPause(canvFin,psnm1.Data(),waitForUser);
842  PlotAlphaBeta(i);
843  PrintAndPause(canvFin,psnm1.Data(),waitForUser);
844  }
845  }
846  MyPrint(2,"Now plotting species");
847  PlotSpecies();
848  PrintAndPause(canvFin,psnm1.Data(),waitForUser);
849  //
850  PrintAndPause(canvFin,psnm2.Data(),false);
851 }
852 
853 void PlotDNDEta(int bin)
854 {
855  MyGuard(1,"Plotting dN/deta in bin %d", bin);
856  //
857  TObjArray *res = &resArr;
858  TString prefN = "bin"; prefN += bin; prefN+="_";
859  int shift = bin*kNHistos;
860  //
861  char buff[1000];
862  // eta range
863  gStyle->SetOptFit(0);
864  gStyle->SetOptStat(0);
865  gStyle->SetOptTitle(0);
866  double mn = 1e6,mx = -1e6;
867  if (!canvFin) canvFin = new TCanvas("canvFin", "canvFin",0,50,700,1000);
868  canvFin->Clear();
869  canvFin->Divide(1,2);
870  canvFin->cd(1);
871  gPad->SetLeftMargin(0.15);
872  //
873  TH1* hZV = (TH1*)res->At(shift + kZvDist);
874  TH1* hZVcorr = (TH1*)res->At(shift + kZvDistCorr);
875  if (!res->At(shift + kSigCorr)) {
876  MyPrint(2,"No histogram at %d", shift+kSigCorr);
877  return;
878  }
879  MyPrint(2,"Got Z vertex and correction histograms: %s %s", hZV->GetName(), hZVcorr->GetName());
880 
881  // corrected data
882  TString nms = prefN;
883  nms += "DataCorrSignal";
884  nms += "_";
885  nms += uniqueName;
886  TH2* tmpSCorr = (TH2F*)res->At(shift + kSigCorr)->Clone(Form("%stmpSCorrBin%d",prefN.Data(),bin));
887  TH1* hsigCorr = 0;
888  TH1* useZV = normToMB ? hZVcorr : hZV;
889  MyPrint(2,"Got corrrected signal at %d (%s): %s", shift + kSigCorr,
890  res->At(shift + kSigCorr)->GetName(), tmpSCorr->GetName());
891  if (useZbinWAv) {
892  // CorrectForZV(tmpSCorr,hZV);
893  CorrectForZV(tmpSCorr,useZV);
894  hsigCorr = ProjectWghMean( tmpSCorr, nms.Data());
895  }
896  else {
897  // hsigCorr = ProjNorm(tmpSCorr,hZV,nms.Data());
898  hsigCorr = ProjNorm(tmpSCorr,useZV,nms.Data());
899  MyPrint(2,"Corrected X signal from normal mean %s", hsigCorr->GetName());
900  }
901  // delete tmpSCorr;
902  SetHStyle(hsigCorr,kRed,20,1.0);
903  hsigCorr->Scale(1./hsigCorr->GetBinWidth(1));
904  hsigCorr->Draw();
905  MyPrint(2,"Scaled %s by bin width, and drew", hsigCorr->GetName());
906  mx = TMath::Max(mx, hsigCorr->GetMaximum());
907  mn = TMath::Min(mn, hsigCorr->GetMinimum());
908  resDnDeta.AddAtAndExpand( hsigCorr, bin );
909  MyPrint(2,"Added %s at %d", hsigCorr->GetName(), bin);
910  hsigCorr->SetDirectory(0);
911  resDnDeta.AddAtAndExpand( res->At(shift + kSigCorr), 100+bin );
912  MyPrint(2,"Added unnormalized %s at %d", hsigCorr->GetName(), bin);
913  resDnDeta.AddAtAndExpand( tmpSCorr, 200+bin );
914  tmpSCorr->SetDirectory(0);
915  resDnDeta.AddAtAndExpand(hZV , 250+bin );
916  MyPrint(2,"Added vertex %s at %d", hZV->GetName(), 250+bin);
917  hZV->SetDirectory(0);
918 
919  //
920  // raw data
921  TH2* tmpRaw = (TH2F*)res->At(shift+kRawDtCut)->Clone("tmpRaw");
922  TH1* hraw = tmpRaw->ProjectionX(prefN+"DataRaw");
923  delete tmpRaw;
924  MyPrint(2,"Got the raw data %s and made projection %s",
925  res->At(shift+kRawDtCut)->GetName(), hraw->GetName());
926  SetHStyle(hraw,kRed,21,1.0);
927  hraw->Scale(1./hraw->GetBinWidth(1));
928  hraw->Draw("same");
929  mn = TMath::Min(mn, hraw->GetMinimum());
930  mx = TMath::Max(mx, hraw->GetMaximum());
931  //
932  // raw data bg sub
933  TH2* tmpRawB = (TH2F*)res->At(shift+kSignalEst)->Clone("tmpRawB");
934  TH1* hraws = tmpRawB->ProjectionX(prefN+"DataRawSub");
935  delete tmpRawB;
936  SetHStyle(hraws,kRed,23,1.0);
937  MyPrint(2,"Got the raw data %s and made projection %s",
938  res->At(shift+kSignalEst)->GetName(), hraws->GetName());
939  hraws->Scale(1./hraw->GetBinWidth(1));
940  hraws->Draw("same");
941  mn = TMath::Min(mn, hraw->GetMinimum());
942  mx = TMath::Max(mx, hraw->GetMaximum());
943  //
944  // bg
945  TH2* tmpBg = (TH2F*)res->At(shift+kBgEst)->Clone("tmpBg");
946  TH1* hbg = tmpBg->ProjectionX(prefN+"BgEst");
947  delete tmpBg;
948  MyPrint(2,"Got the bg estimate %s and made projection %s",
949  res->At(shift+kBgEst)->GetName(), hbg->GetName());
950  SetHStyle(hbg,kMagenta,22,1.0);
951  hbg->Scale(1./hbg->GetBinWidth(1));
952  hbg->Draw("same");
953  mn = TMath::Min(mn, hbg->GetMinimum());
954  mx = TMath::Max(mx, hbg->GetMaximum());
955  //
956  // mc part ----------------------------
957  TH1* hZVMC = (TH1*)res->At(shift + kZvDist + kMCShift);
958  TH1* hZVMCNS = (TH1*)res->At(shift + kZvMCDistNS + kMCShift);
959  //
960  // raw data
961  TH2* tmpRawMC = (TH2F*)res->At(shift+kRawDtCut+kMCShift)->Clone("tmpRawMC");
962  TH1* hrawMC = tmpRawMC->ProjectionX(prefN+"DataRawMC");
963  delete tmpRawMC;
964  MyPrint(2,"Got the MC raw %s and made projection %s",
965  res->At(shift+kRawDtCut+kMCShift)->GetName(), hrawMC->GetName());
966  SetHStyle(hrawMC,kBlue,24,1.0);
967  hrawMC->Scale(1./hrawMC->GetBinWidth(1));
968  hrawMC->Draw("same");
969  mn = TMath::Min(mn, hrawMC->GetMinimum());
970  mx = TMath::Max(mx, hrawMC->GetMaximum());
971  //
972  // raw data bg sub
973  TH2* tmpRawBMC = (TH2F*)res->At(shift+kSignalEst+kMCShift)->Clone("tmpRawBMC");
974  TH1* hrawsMC = tmpRawBMC->ProjectionX(prefN+"DataRawSubMC");
975  delete tmpRawBMC;
976  MyPrint(2,"Got the MC signal estimate %s and made projection %s",
977  res->At(shift+kSignalEst+kMCShift)->GetName(), hrawsMC->GetName());
978  SetHStyle(hrawsMC,kBlue,26,1.0);
979  hrawsMC->Scale(1./hrawMC->GetBinWidth(1));
980  hrawsMC->Draw("same");
981  mn = TMath::Min(mn, hrawMC->GetMinimum());
982  mx = TMath::Max(mx, hrawMC->GetMaximum());
983  //
984  // raw data bgMClabels sub
985  TH2* tmpRawLBMC = (TH2F*)res->At(shift+kSignalEstMC+kMCShift)->Clone("tmpRawLBMC");
986  TH1* hrawsMCLB = tmpRawLBMC->ProjectionX(prefN+"DataRawSubMCLB");
987  delete tmpRawLBMC;
988  MyPrint(2,"Got the MC signal estimate %s from labels and made projection %s",
989  res->At(shift+kSignalEstMC+kMCShift)->GetName(), hrawsMCLB->GetName());
990  SetHStyle(hrawsMCLB,kGreen+2,30,1.0);
991  hrawsMCLB->Scale(1./hrawsMCLB->GetBinWidth(1));
992  hrawsMCLB->Draw("same");
993  mn = TMath::Min(mn, hrawsMCLB->GetMinimum());
994  mx = TMath::Max(mx, hrawsMCLB->GetMaximum());
995  //
996  // bg est
997  TH2* tmpBgMCEst = (TH2F*)res->At(shift+kBgEst+kMCShift)->Clone("tmpRawBbMCEst");
998  TH1* hbgMCEst = tmpBgMCEst->ProjectionX(prefN+"BgEstMC");
999  MyPrint(2,"Got the MC bg estimate %s and made projection %s",
1000  res->At(shift+kBgEst+kMCShift)->GetName(), hbgMCEst->GetName());
1001  delete tmpBgMCEst;
1002  SetHStyle(hbgMCEst,kBlue,26,1.0);
1003  hbgMCEst->Scale(1./hbgMCEst->GetBinWidth(1));
1004  hbgMCEst->Draw("same");
1005  mn = TMath::Min(mn, hbgMCEst->GetMinimum());
1006  mx = TMath::Max(mx, hbgMCEst->GetMaximum());
1007  //
1008  // bg MC
1009  TH2* tmpBgMC = (TH2F*)res->At(shift+kBgMC+kMCShift)->Clone("tmpBgMC");
1010  TH1* hbgMC = tmpBgMC->ProjectionX(prefN+"BgMC");
1011  delete tmpBgMC;
1012  MyPrint(2,"Got the MC bg %s and made projection %s",
1013  res->At(shift+kBgMC+kMCShift)->GetName(), hbgMC->GetName());
1014  SetHStyle(hbgMC,kGreen+2,25,1.0);
1015  hbgMC->Scale(1./hbgMC->GetBinWidth(1));
1016  hbgMC->Draw("same");
1017  mn = TMath::Min(mn, hbgMC->GetMinimum());
1018  mx = TMath::Max(mx, hbgMC->GetMaximum());
1019  //
1020  // mc truth
1021  TH2* tmpMCTrue = (TH2F*)res->At(shift+kMCPrim+kMCShift)->Clone(Form("%stmpMCTrueBin%d",prefN.Data(),bin));
1022  TH1* hMCtrue = 0;
1023  TH1* hzMC = normToMB ? hZVMCNS : hZVMC;
1024  if (useZbinWAv) {
1025  CorrectForZV(tmpMCTrue,hzMC);
1026  hMCtrue = ProjectWghMean(tmpMCTrue, prefN+"MCTruth");
1027  }
1028  else {
1029  hMCtrue = ProjNorm(tmpMCTrue,hzMC,prefN+"MCTruth");
1030  // hMCtrue = ProjNorm(tmpMCTrue,hZVMCNS,prefN+"MCTruth");
1031  }
1032  // delete tmpMCTrue;
1033  MyPrint(2,"Got the MC truth %s and made projection %s",
1034  res->At(shift+kMCPrim+kMCShift)->GetName(), hMCtrue->GetName());
1035  SetHStyle(hMCtrue,kCyan,20,0.8);
1036  hMCtrue->Scale(1./hMCtrue->GetBinWidth(1));
1037  hMCtrue->SetLineStyle(1);
1038  hMCtrue->Draw("same");
1039  mn = TMath::Min(mn, hMCtrue->GetMinimum());
1040  mx = TMath::Max(mx, hMCtrue->GetMaximum());
1041  //
1042  hMCtrue->SetDirectory(0);
1043  tmpMCTrue->SetDirectory(0);
1044  resDnDeta.AddAtAndExpand( hMCtrue, 300+bin);
1045  resDnDeta.AddAtAndExpand( res->At(shift+kMCPrim+kMCShift), 400+bin );
1046  resDnDeta.AddAtAndExpand( tmpMCTrue, 500+bin );
1047  resDnDeta.AddAtAndExpand( hZVMCNS, 550+bin );
1048  //
1049 
1050  mn = 0;
1051  hsigCorr->SetMinimum(mn);
1052  hsigCorr->SetMaximum(mx + 0.4*(mx-mn));
1053  //
1054  char ftres[1000];
1055  sprintf(ftres,"dN/d#eta_{|#eta|<%.2f} = %.2f #pm %.2f",kEtaFitRange,dNdEta[bin],dNdEtaErr[bin]);
1056  double xtxt = (hsigCorr->GetXaxis()->GetXmin()+hsigCorr->GetXaxis()->GetXmax())/2;
1057  TLatex *txfit = new TLatex(xtxt,mn+0.25*(mx-mn), ftres);
1058  txfit->SetTextSize(0.04);
1059  txfit->Draw();
1060 
1061 
1062  gPad->Modified();
1063  //
1064  MyPrint(2,"Making the legend");
1065  TLegend *legDnDeta = new TLegend(0.15,0.75, 0.45,0.95);
1066  legDnDeta->SetFillColor(kWhite);
1067  legDnDeta->SetHeader("Data");
1068  legDnDeta->AddEntry(hsigCorr,"Corrected","pl");
1069  legDnDeta->AddEntry(hraw, "Reconstructed","pl");
1070  sprintf(buff,"Reconstructed - Bckg.%s.",useBgType.Data());
1071  legDnDeta->AddEntry(hraws, buff,"pl");
1072  sprintf(buff,"Background %s.",useBgType.Data());
1073  legDnDeta->AddEntry(hbg, buff,"pl");
1074  legDnDeta->Draw();
1075  //
1076  TLegend *legDnDetaMC = new TLegend(0.60,0.72, 0.95,0.95);
1077  legDnDetaMC->SetFillColor(kWhite);
1078  legDnDetaMC->SetHeader("MC");
1079  legDnDetaMC->AddEntry(hrawMC, "Reconstructed","pl");
1080  sprintf(buff,"Reconstructed - Bckg.%s.",useBgType.Data());
1081  legDnDetaMC->AddEntry(hrawsMC, buff,"pl");
1082  sprintf(buff,"Reconstructed - Bckg.%s.","MC.Labels");
1083  legDnDetaMC->AddEntry(hrawsMCLB, buff,"pl");
1084  sprintf(buff,"Background %s.",useBgType.Data());
1085  legDnDetaMC->AddEntry(hbgMCEst, buff,"pl");
1086  sprintf(buff,"Background %s.","MC Labels");
1087  legDnDetaMC->AddEntry(hbgMC, buff,"pl");
1088  //
1089  legDnDetaMC->Draw();
1090  //
1091  gPad->SetGrid(1.1);
1092  gPad->Modified();
1093  AddLabel(outTitle,0.1,0.97, kBlack,0.02);
1094  //
1095  canvFin->cd();
1096  //
1097  //---------------- dsitributions
1098  canvFin->cd(2);
1099  //
1100  MyPrint(2,"Drawing AUX distributions");
1101  TH1* mcdst = (TH1*)res->At(shift+kDataDist+kMCShift);
1102  TH1* mcdstbg = (TH1*)res->At(shift+kBgDist+kMCShift);
1103  TH1* mcdstbgLB = (TH1*)res->At(shift+kBgMCDist+kMCShift);
1104  TH1* dtdst = (TH1*)res->At(shift+kDataDist);
1105  TH1* dtdstbg = (TH1*)res->At(shift+kBgDist);
1106  //
1107  TH1* mcDstN = (TH1*)FindObject(bin,HName("Data", useShapeType==kNormShapeDist ? "WDist":"DPhiS"), listMC );
1108  TH1* mcDstSec = (TH1*)FindObject(bin,HName("Sec", useShapeType==kNormShapeDist ? "WDist":"DPhiS"), listMC );
1109  TH1* mcDstCombU = (TH1*)FindObject(bin,HName("CombU",useShapeType==kNormShapeDist ? "WDist":"DPhiS"), listMC );
1110  TH1* mcDstCombC = (TH1*)FindObject(bin,HName("Comb", useShapeType==kNormShapeDist ? "WDist":"DPhiS"), listMC );
1111  //
1112  double scl,sclE;
1113  mcDstN = NormalizeBg(mcdst,mcDstN,scl,sclE);
1114  mcDstSec->Scale(scl);
1115  mcDstCombU->Scale(scl);
1116  mcDstCombC->Scale(scl);
1117  mcDstCombC->Add(mcDstCombU,-1);
1118 
1119  dtdst->Draw("");
1120  gPad->Modified();
1121  dtdst->GetXaxis()->SetLabelSize(0.03);
1122  dtdst->GetXaxis()->SetTitleSize(0.03);
1123  dtdst->GetXaxis()->SetTitleOffset(2);
1124  dtdstbg->Draw("same");
1125  //
1126  mcdst->Draw("same");
1127  mcDstSec->Draw("same");
1128  if (mcdstbgLB) mcdstbgLB->Draw("same");
1129  mcdstbg->Draw("same");
1130  mcDstCombC->Draw("same");
1131  //
1132 
1133  SetHStyle(mcdst,kBlue, 25,0.7);
1134  if (mcdstbgLB) SetHStyle(mcdstbgLB,kGreen, 7/*32*/,0.5);
1135  SetHStyle(mcdstbg,kCyan, 7/*26*/,0.5);
1136  SetHStyle(mcDstCombC,kGreen+2, 21,0.7);
1137  SetHStyle(mcDstSec,kBlue+2, 22,0.7);
1138  //
1139  SetHStyle(dtdst,kRed, 20,0.7);
1140  SetHStyle(dtdstbg,kBlue, 34,0.7);
1141  //
1142  double vmcTot,vmcTotE;
1143  double vmcSec,vmcSecE, ratSec,ratSecE;
1144  double vmcCmbEst,vmcCmbEstE, ratCmbEst,ratCmbEstE;
1145  double vmcCmb,vmcCmbE, ratCmb,ratCmbE;
1146  double vmcCmbC,vmcCmbCE, ratCmbC,ratCmbCE;
1147  double cutSgMin,cutSgMax;
1148  double cutBgMin,cutBgMax;
1149  if (useShapeType==kNormShapeDist) {
1150  cutSgMin = 0;
1151  cutSgMax = kWDistSgCut;
1152  cutBgMin = kWDistBgTailMin;
1153  cutBgMax = kWDistBgTailMax;
1154  }
1155  else {
1156  cutSgMin = 0;
1157  cutSgMax = kdPhiSgCut;
1158  cutBgMin = kdPhiBgTailMin;
1159  cutBgMax = kdPhiBgTailMax;
1160  }
1161  Integrate(mcdst, cutSgMin,cutSgMax, vmcTot,vmcTotE);
1162  Integrate(mcDstSec, cutSgMin,cutSgMax, vmcSec,vmcSecE);
1163  GetRatE(vmcSec,vmcSecE, vmcTot,vmcTotE, ratSec,ratSecE);
1164  //
1165  if (mcdstbgLB) {
1166  Integrate(mcdstbgLB, cutSgMin,cutSgMax, vmcCmb,vmcCmbE);
1167  GetRatE(vmcCmb,vmcCmbE, vmcTot,vmcTotE, ratCmb,ratCmbE);
1168  }
1169  //
1170  Integrate(mcdstbg, cutSgMin,cutSgMax, vmcCmbEst,vmcCmbEstE);
1171  GetRatE(vmcCmbEst,vmcCmbEstE, vmcTot,vmcTotE, ratCmbEst,ratCmbEstE);
1172  //
1173  Integrate(mcDstCombC, cutSgMin,cutSgMax, vmcCmbC,vmcCmbCE);
1174  GetRatE(vmcCmbC,vmcCmbCE, vmcTot,vmcTotE, ratCmbC,ratCmbCE);
1175  //
1176  double vdtTot,vdtTotE;
1177  double vdtBg,vdtBgE, ratdtBg,ratdtBgE;
1178  //
1179  Integrate(dtdst, cutSgMin,cutSgMax, vdtTot,vdtTotE);
1180  Integrate(dtdstbg, cutSgMin,cutSgMax, vdtBg,vdtBgE);
1181  GetRatE( vdtBg,vdtBgE, vdtTot,vdtTotE, ratdtBg,ratdtBgE);
1182  //
1183  //
1184  double dmn = mcdst->GetMinimum();
1185  double dmx = mcdst->GetMaximum();
1186  TLine *ln = new TLine(cutSgMax, dmn, cutSgMax, dmx);
1187  ln->SetLineColor(kBlack);
1188  ln->Draw();
1189  TLine *lnc = new TLine(cutBgMin, dmn, cutBgMin, dmx);
1190  lnc->SetLineColor(kRed);
1191  lnc->Draw();
1192  if (useShapeType==kNormShapeDPhi) {
1193  ln = new TLine(-cutSgMax, dmn, -cutSgMax, dmx);
1194  ln->SetLineColor(kBlack);
1195  ln->Draw();
1196  //
1197  lnc = new TLine(-cutBgMin, dmn, -cutBgMin, dmx);
1198  lnc->SetLineColor(kRed);
1199  lnc->Draw();
1200  }
1201  //
1202  TLegend *legDstMC1 = new TLegend(0.60,0.72, 0.95,0.95);
1203  legDstMC1->SetFillColor(kWhite);
1204 
1205  //
1206  legDstMC1->AddEntry(dtdst, "Data raw","pl");
1207  sprintf(buff,"Data Comb. %s. : %.1f%%",useBgType.Data(),ratdtBg*100);
1208  legDstMC1->AddEntry(dtdstbg, buff,"pl");
1209  //
1210 
1211 
1212  legDstMC1->AddEntry(mcdst, "MC raw","pl");
1213  sprintf(buff,"MC Comb. %s. : %.1f%%",useBgType.Data(),ratCmbEst*100);
1214  legDstMC1->AddEntry(mcdstbg, buff,"pl");
1215  //
1216  sprintf(buff,"MC Comb. %s. : %.1f%%","MC Lbl.",ratCmb*100);
1217  legDstMC1->AddEntry(mcdstbgLB, buff,"pl");
1218 
1219  sprintf(buff,"MC Comb.Uncorr %s. : %.1f%%","MC Lbl.",ratCmbC*100);
1220  legDstMC1->AddEntry(mcDstCombC, buff,"pl");
1221 
1222  sprintf(buff,"MC Sec. : %.1f%%",ratSec*100);
1223  legDstMC1->AddEntry(mcDstSec, buff,"pl");
1224 
1225  legDstMC1->Draw();
1226 
1227  if (useShapeType==kNormShapeDist) gPad->SetLogx();
1228  gPad->SetLogy();
1229  gPad->SetGrid(1,1);
1230  gPad->Modified();
1231  //
1232  canvFin->cd();
1233  //
1234  if (creatDnDEtaCMacro) {
1235  sprintf(buff,"%s/%s_b%d_dNdEta_%s",figDir,uniqueName.Data(),bin,outStr);
1236  SaveCanvas(canvFin,buff,"cg");
1237  }
1238  //
1239 }
1240 //
1241 void PlotAlphaBeta(int bin)
1242 {
1243  MyGuard(2,"Draw alpha and beta");
1244  char buff[1000];
1245  int shift = bin*kNHistos;
1246  gStyle->SetOptFit(0);
1247  gStyle->SetOptStat(0);
1248  gStyle->SetOptTitle(0);
1249  TObjArray* res = &resArr;
1250  //------------------------------------------------------
1251  if (!canvFin) canvFin = new TCanvas("canvFin","canvFin",10,10,700,1000);
1252  canvFin->Clear();
1253  canvFin->Divide(2,3,0.01,0.01);
1254  canvFin->cd(1);
1255  TH1* dtBet = (TH1*)res->At(shift + k1MBeta);
1256  TH1* mcBet = (TH1*)res->At(shift + k1MBeta+kMCShift);
1257  TH1* mcBetLB = (TH1*)res->At(shift + k1MBetaMC+kMCShift);
1258  if (!dtBet) return;
1259  double mn,mx,mnt,mxt;
1260  GetRealMinMax(dtBet,mn,mx);
1261  GetRealMinMax(mcBet,mnt,mxt);
1262  if (mnt<mn) mn = mnt;
1263  if (mxt>mx) mx = mxt;
1264  GetRealMinMax(mcBetLB,mnt,mxt);
1265  if (mnt<mn) mn = mnt;
1266  if (mxt>mx) mx = mxt;
1267  //
1268  dtBet->SetMinimum(mn - 0.05*(mx-mn));
1269  dtBet->SetMaximum(mx + 0.05*(mx-mn));
1270  mcBet->SetMinimum(mn - 0.05*(mx-mn));
1271  mcBet->SetMaximum(mx + 0.05*(mx-mn));
1272  mcBetLB->SetMinimum(mn - 0.05*(mx-mn));
1273  mcBetLB->SetMaximum(mx + 0.05*(mx-mn));
1274  //
1275  canvFin->cd(1);
1276  gPad->SetRightMargin(0.15);
1277  dtBet->Draw("colz");
1278  AddLabel("#beta Data",0.2,0.95,kBlack,0.04);
1279  gPad->Modified();
1280  dtBet->GetYaxis()->SetTitleOffset(1.4);
1281 
1282  // TPaletteAxis *p = (TPaletteAxis*)dtBet->FindObject("palette");
1283  // if (p) p->SetX1NDC(0.85);
1284  canvFin->cd(2);
1285  gPad->SetRightMargin(0.15);
1286  mcBet->Draw("colz");
1287  AddLabel("#beta MC (bckg.estimated)",0.2,0.95,kBlack,0.04);
1288  gPad->Modified();
1289  mcBet->GetYaxis()->SetTitleOffset(1.4);
1290  // p = (TPaletteAxis*)mcBet->FindObject("palette");
1291  // if (p) p->SetX1NDC(0.85);
1292  canvFin->cd(3);
1293  gPad->SetRightMargin(0.15);
1294  mcBetLB->Draw("colz");
1295  AddLabel("#beta MC (bckg.from MC labels)",0.2,0.95,kBlack,0.04);
1296  gPad->Modified();
1297  mcBetLB->GetYaxis()->SetTitleOffset(1.4);
1298  // p = (TPaletteAxis*)mcBetLB->FindObject("palette");
1300  //
1301  //------------------------------------------------------
1302  TH1* dtAlp = (TH1*)res->At(shift + kAlpha);
1303  TH1* mcAlp = (TH1*)res->At(shift + kAlphaMC);
1304  GetRealMinMax(dtAlp,mn,mx);
1305  GetRealMinMax(mcAlp,mnt,mxt);
1306  if (mnt<mn) mn = mnt;
1307  if (mxt>mx) mx = mxt;
1308  dtAlp->SetMinimum(mn - 0.05*(mx-mn));
1309  dtAlp->SetMaximum(mx + 0.05*(mx-mn));
1310  mcAlp->SetMinimum(mn - 0.05*(mx-mn));
1311  mcAlp->SetMaximum(mx + 0.05*(mx-mn));
1312  //
1313  canvFin->cd(4);
1314  gPad->SetRightMargin(0.15);
1315  dtAlp->Draw("colz");
1316  AddLabel("#alpha (bckg.estimated)",0.2,0.95,kBlack,0.04);
1317  gPad->Modified();
1318  dtAlp->GetYaxis()->SetTitleOffset(1.4);
1319  // TPaletteAxis *pa = (TPaletteAxis*)dtBet->FindObject("palette");
1320  // if (pa) pa->SetX1NDC(0.85);
1321  canvFin->cd(5);
1322  gPad->SetRightMargin(0.15);
1323  mcAlp->Draw("colz");
1324  AddLabel("#alpha (bckg.from MC labels)",0.2,0.95,kBlack,0.04);
1325  gPad->Modified();
1326  mcAlp->GetYaxis()->SetTitleOffset(1.4);
1327  // pa = (TPaletteAxis*)mcBet->FindObject("palette");
1328  // if (pa) pa->SetX1NDC(0.85);
1329  gPad->Modified();
1330  canvFin->cd(6);
1331  AddLabel(outTitle,0.1,0.5, kBlack, 0.02);
1332  //
1333  if (creatAlphaBetaCMacro) {
1334  sprintf(buff,"%s/%sAlphaBeta_%s",figDir,uniqueName.Data(),outStr);
1335  SaveCanvas(canvFin,buff,"cg");
1336  }
1337  //
1338 }
1339 
1340 void PlotSpecies()
1341 {
1342  MyGuard(2,"Plotting species");
1343  char buff[1000];
1344  gStyle->SetOptFit(0);
1345  gStyle->SetOptStat(0);
1346  gStyle->SetOptTitle(0);
1347  //------------------------------------------------------
1348  TH2F* hSpecPrim = (TH2F*)FindObject(-1, "pdgPrim", listMC,kFALSE);
1349  TH2F* hSpecSec = (TH2F*)FindObject(-1, "pdgSec", listMC,kFALSE);
1350  TH2F* hSpecPrimP = (TH2F*)FindObject(-1, "pdgPrimPar", listMC,kFALSE);
1351  TH2F* hSpecSecP = (TH2F*)FindObject(-1, "pdgSecPar", listMC,kFALSE);
1352  int nbd = hSpecPrim->GetXaxis()->GetNbins();
1353  //
1354  TH1* hSpecPrimAll = hSpecPrim->ProjectionX("specPrimAll",1,nbd+1,"e");
1355  hSpecPrimAll->Scale(100./hSpecPrimAll->Integral());
1356  hSpecPrimAll->GetYaxis()->SetTitle("Fraction,%");
1357  hSpecPrimAll->GetXaxis()->SetLabelSize(0.06);
1358  hSpecPrimAll->GetXaxis()->LabelsOption("v");
1359  //
1360  TH1* hSpecSecAll = hSpecSec->ProjectionX("specSecAll",1,nbd+1,"e");
1361  hSpecSecAll->Scale(100./hSpecSecAll->Integral());
1362  hSpecSecAll->GetYaxis()->SetTitle("Fraction,%");
1363  hSpecSecAll->GetXaxis()->SetLabelSize(0.05);
1364  //
1365  TH1* hSpecPrimPAll = hSpecPrimP->ProjectionX("specPrimPAll",1,nbd+1,"e");
1366  hSpecPrimPAll->Scale(100./hSpecPrimPAll->Integral());
1367  hSpecPrimPAll->GetYaxis()->SetTitle("Fraction,%");
1368  hSpecPrimPAll->GetXaxis()->SetLabelSize(0.06);
1369  hSpecPrimPAll->GetXaxis()->LabelsOption("v");
1370 
1371  //
1372  TH1* hSpecSecPAll = hSpecSecP->ProjectionX("specSecPAll",1,nbd+1,"e");
1373  hSpecSecPAll->Scale(100./hSpecSecPAll->Integral());
1374  hSpecSecPAll->GetYaxis()->SetTitle("Fraction,%");
1375  hSpecSecPAll->GetXaxis()->SetLabelSize(0.05);
1376  //
1377  int binCut = hSpecPrim->GetXaxis()->FindBin(kWDistSgCut-kEps);
1378  TH1* hSpecPrimSel = hSpecPrim->ProjectionX("specPrimSel",1,binCut,"e");
1379  if (hSpecPrimSel->Integral()<1) return;
1380  hSpecPrimSel->Scale(100./hSpecPrimSel->Integral());
1381  hSpecPrimSel->GetYaxis()->SetTitle("Fraction,%");
1382  hSpecPrimSel->GetXaxis()->SetLabelSize(0.05);
1383  //
1384  TH1* hSpecSecSel = hSpecSec->ProjectionX("specSecSel",1,binCut,"e");
1385  hSpecSecSel->Scale(100./hSpecSecSel->Integral());
1386  hSpecSecSel->GetYaxis()->SetTitle("Fraction,%");
1387  hSpecSecSel->GetXaxis()->SetLabelSize(0.05);
1388  //
1389  TH1* hSpecPrimPSel = hSpecPrimP->ProjectionX("specPrimPSel",1,binCut,"e");
1390  hSpecPrimPSel->Scale(100./hSpecPrimPSel->Integral());
1391  hSpecPrimPSel->GetYaxis()->SetTitle("Fraction,%");
1392  hSpecPrimPSel->GetXaxis()->SetLabelSize(0.05);
1393  //
1394  TH1* hSpecSecPSel = hSpecSecP->ProjectionX("specSecPSel",1,binCut,"e");
1395  hSpecSecPSel->Scale(100./hSpecSecPSel->Integral());
1396  hSpecSecPSel->GetYaxis()->SetTitle("Fraction,%");
1397  hSpecSecPSel->GetXaxis()->SetLabelSize(0.05);
1398  //
1399  if (!canvFin) canvFin = new TCanvas("canvFin","canvFin",10,10,1100,800);
1400  canvFin->Clear();
1401  canvFin->Divide(1,2,0.01,0.01);
1402  canvFin->cd(1);
1403  hSpecPrimAll->Draw();
1404  SetHStyle(hSpecPrimAll,kBlue,25,1.1);
1405  hSpecPrimSel->Draw("same");
1406  SetHStyle(hSpecPrimSel,kRed,20,1);
1407  //
1408  hSpecSecAll->Draw("same");
1409  SetHStyle(hSpecSecAll,kGreen,32,1.1);
1410  hSpecSecSel->Draw("same");
1411  SetHStyle(hSpecSecSel,kBlack,22,1);
1412  //
1413  TLegend *legPart = new TLegend(0.8,0.72, 0.999,0.999);
1414  legPart->SetFillColor(kWhite);
1415  legPart->SetHeader("Tracklet PDG");
1416  //
1417  legPart->AddEntry(hSpecPrimAll, "Prim., before #Delta cut","pl");
1418  legPart->AddEntry(hSpecPrimSel, "Prim., after #Delta cut","pl");
1419  legPart->AddEntry(hSpecSecAll, "Sec., before #Delta cut","pl");
1420  legPart->AddEntry(hSpecSecSel, "Sec., after #Delta cut","pl");
1421  //
1422  legPart->Draw();
1423  gPad->SetLogy();
1424  gPad->SetGrid(1,1);
1425  gPad->Modified();
1426  //
1427  canvFin->cd(2);
1428  hSpecPrimPAll->Draw();
1429  SetHStyle(hSpecPrimPAll,kBlue,25,1.1);
1430  hSpecPrimPSel->Draw("same");
1431  SetHStyle(hSpecPrimPSel,kRed,20,1);
1432  //
1433  hSpecSecPAll->Draw("same");
1434  SetHStyle(hSpecSecPAll,kGreen,32,1.1);
1435  hSpecSecPSel->Draw("same");
1436  SetHStyle(hSpecSecPSel,kBlack,22,1);
1437  //
1438  TLegend *legPartP = new TLegend(0.8,0.72, 0.999,0.999);
1439  legPartP->SetFillColor(kWhite);
1440  legPartP->SetHeader("Tracklet Parents PDG");
1441  //
1442  legPartP->AddEntry(hSpecPrimPAll, "Prim., before #Delta cut","pl");
1443  legPartP->AddEntry(hSpecPrimPSel, "Prim., after #Delta cut","pl");
1444  legPartP->AddEntry(hSpecSecPAll, "Sec., before #Delta cut","pl");
1445  legPartP->AddEntry(hSpecSecPSel, "Sec., after #Delta cut","pl");
1446  //
1447  legPartP->Draw();
1448  gPad->SetLogy();
1449  gPad->SetGrid(1,1);
1450  gPad->Modified();
1451  //
1452  canvFin->cd(1);
1453  // AddLabel(outTitle,0.1,0.97, kBlack, 0.02);
1454  canvFin->cd();
1455  //
1456  if (creatSpeciesCMacro) {
1457  sprintf(buff,"%s/%sSpecies_%s",figDir,uniqueName.Data(),outStr);
1458  SaveCanvas(canvFin,buff,"cg");
1459  }
1460 }
1461 
1462 //______________________________________________________________________
1463 void CropHisto(TH1* histo, int bx0, int bx1, int by0, int by1)
1464 {
1465  MyGuard(3,"Cropping histogram %s between bins %d,%d x %d,%d", bx0, bx1, by0, by1);
1466  // fill 0 to all bins outside defined range
1467  TAxis *xax = histo->GetXaxis();
1468  int nbx = xax->GetNbins();
1469  double vmn=1e16,vmx=-1e16;
1470  if (histo->InheritsFrom(TH2::Class())) {
1471  TAxis *yax = histo->GetYaxis();
1472  int nby = yax->GetNbins();
1473  for (int ix=nbx+2;ix--;) {
1474  for (int iy=nby+2;iy--;) {
1475  if ((ix<bx0||ix>bx1)||(iy<by0||iy>by1)) {
1476  histo->SetBinContent(ix,iy,0);
1477  histo->SetBinError(ix,iy,0);
1478  }
1479  else {
1480  double vl = histo->GetBinContent(ix,iy);
1481  if (vl<vmn) vmn = vl;
1482  if (vl>vmx) vmx = vl;
1483  }
1484  }
1485  }
1486  }
1487  else {
1488  for (int ix=nbx+2;ix--;) {
1489  if ((ix<bx0||ix>bx1)) {
1490  histo->SetBinContent(ix,0);
1491  histo->SetBinError(ix,0);
1492  }
1493  else {
1494  double vl = histo->GetBinContent(ix);
1495  if (vl<vmn) vmn = vl;
1496  if (vl>vmx) vmx = vl;
1497  }
1498  }
1499  }
1500  //
1501  if (vmn==vmx) {
1502  vmn = 0.95*vmn;
1503  vmx = 1.05*vmx;
1504  }
1505  histo->SetMaximum(vmx);
1506  histo->SetMinimum(vmn);
1507 }
1508 
1509 //______________________________________________________________________
1510 void CropHisto(TH1* histo, double vx0, double vx1, double vy0, double vy1)
1511 {
1512  MyGuard(3,"Cropping histogram %s between values %f,%f x %f,%f", vx0, vx1, vy0, vy1);
1513  // fill 0 to all bins outside defined range
1514  TAxis *xax = histo->GetXaxis();
1515  int bx0,bx1,by0=-1,by1=-1;
1516  bx0 = xax->FindBin(vx0+kEps);
1517  bx1 = xax->FindBin(vx1-kEps);
1518  if (histo->InheritsFrom(TH2::Class())) {
1519  TAxis *yax = histo->GetYaxis();
1520  by0 = yax->FindBin(vy0+kEps);
1521  by1 = yax->FindBin(vy1-kEps);
1522  }
1523  CropHisto(histo,bx0,bx1,by0,by1);
1524 }
1525 
1526 //______________________________________________________________________
1527 TH1* NormalizeBg(TH1* dataH, TH1* bgH, double &scl, double &sclE)
1528 {
1529  // match generated bg and data tails, calculate normalization,
1530  // return normalized bg copy
1531  //
1532  MyGuard(1,"Normalize background from %s and %s",
1533  dataH->GetName(), bgH->GetName());
1534 
1535  TAxis* xax = dataH->GetXaxis();
1536  int nbtot = xax->GetNbins();
1537  int bgBins[2][2] = {{0}}; // limiting bins for tails integration
1538  Int_t ntails; // 0 for dphi plot, 1 for weighted dist plot
1539  if (useShapeType == kNormShapeDist) { // only positive tail
1540  bgBins[0][0] = xax->FindBin(kWDistBgTailMin+kEps); // positive tail min bin
1541  bgBins[0][1] = xax->FindBin(kWDistBgTailMax-kEps); // positive tail max bin
1542  ntails = 1;
1543  }
1544  else if (useShapeType == kNormShapeDPhi) { // both tails
1545  bgBins[0][0] = xax->FindBin( kdPhiBgTailMin+kEps); // positive tail min bin
1546  bgBins[0][1] = xax->FindBin( kdPhiBgTailMax-kEps); // positive tail max bin
1547  bgBins[1][0] = xax->FindBin(-kdPhiBgTailMax+kEps); // negative tail min bin
1548  bgBins[1][1] = xax->FindBin(-kdPhiBgTailMin-kEps); // positive tail max bin
1549  ntails = 2;
1550  }
1551  else {
1552  MyPrint(2,"NormalizeBg: unknown shape type %d",useShapeType);
1553  return 0;
1554  }
1555  MyPrint(2,"NormalizeBg: bins for tails: right: %d:%d / left: %d:%d",
1556  bgBins[0][0],bgBins[0][1],bgBins[1][0],bgBins[1][1]);
1557  //
1558  double meanR=0,meanRE=0,meanRE2=0;
1559  double meanD=0,meanDE2=0;
1560  double meanB=0,meanBE2=0;
1561  double meanRI=0,meanRIE=0;
1562  for (int itp=0;itp<=ntails;itp++) {
1563  for (int ib=bgBins[itp][0];ib<=bgBins[itp][1];ib++) {
1564  if (ib<1||ib>nbtot) continue;
1565  double vD = dataH->GetBinContent(ib);
1566  double vB = bgH->GetBinContent(ib);
1567  double eD = dataH->GetBinError(ib);
1568  double eB = bgH->GetBinError(ib);
1569  meanD += vD; meanDE2 += eD*eD;
1570  meanB += vB; meanBE2 += eB*eB;
1571  if (vD<=0 || vB<=0 || eD<=0 || eB<=0) continue;
1572  double rat = vD/vB;
1573  double ratE2 = rat*rat*(eD*eD/vD/vD + eB*eB/vB/vB);
1574  meanR += rat/ratE2; meanRE2 += 1.0/ratE2;
1575  }
1576  }
1577  //
1578  if (meanRE2>0) {
1579  meanR /= meanRE2;
1580  meanRE2 = 1./meanRE2;
1581  meanRE = TMath::Sqrt(meanRE2);
1582  }
1583  if (meanDE2>0 && meanBE2>0) {
1584  meanRI = meanD/meanB;
1585  meanRIE = meanRI*TMath::Sqrt(meanDE2/meanD/meanD + meanBE2/meanB/meanB);
1586  }
1587  MyPrint(1,"NormalizeBg: Tails scaling %s wrt %s: "
1588  "Wgh.Mean:%.4f(%.4f) / Integral:%.4f(%.4f)",
1589  bgH->GetName(),dataH->GetName(), meanR,meanRE, meanRI,meanRIE);
1590  MyPrint(2,"NormalizeBg: Select scaling type %s",
1591  useScaleType==kSclWghMean ? "Wgh.Mean":"Integral");
1592  //
1593  scl = useScaleType==kSclWghMean ? meanR : meanRI;
1594  sclE = useScaleType==kSclWghMean ? meanRE : meanRIE;
1595  //
1596  // rescaled bg
1597  char buff[1000];
1598  sprintf(buff,"%s_bgNorm",bgH->GetName());
1599  TH1* tmp = bgH;
1600  bgH = (TH1*)tmp->Clone(buff);
1601  sprintf(buff,"%s bgNorm%d %.4f+-%.4f",bgH->GetName(),useScaleType,scl,sclE);
1602  TH1* dumH = (TH1*)bgH->Clone("dummySCL$"); dumH->Reset();
1603  for (int i=1;i<=nbtot;i++) {
1604  dumH->SetBinContent(i,scl);
1605  dumH->SetBinError(i,sclE);
1606  }
1607  bgH->Multiply(dumH);
1608  MyPrint(1,"Returning normal backg delta distribution %s (copy of %s)",
1609  bgH->GetName(), tmp->GetName());
1610  delete dumH;
1611  return bgH;
1612 }
1613 
1614 //______________________________________________________________________
1615 TObject* FindObject(int bin, const char* nameH, const TList* lst, Bool_t normPerEvent)
1616 {
1617  // get histo, optionally normalizing it per processed event
1618  if (!lst) {
1619  Warning("FindObject","%s: No list provided",nameH);
1620  return 0;
1621  }
1622 
1623  int nent = lst->GetEntries();
1624  char buff[200];
1625  if (bin>=0) sprintf(buff,"b%d_%s",bin,nameH);
1626  else sprintf(buff,"%s",nameH);
1627 
1628  MyGuard(1,"Looking for bin %d of %s: %s", bin, nameH, buff);
1629 
1630  TString nm;
1631  TObject *hst = 0;
1632  for (int i=nent;i--;) {
1633  nm = lst->At(i)->GetName();
1634  if (nm.EndsWith(buff)) {hst = lst->At(i); break;}
1635  }
1636  if (!hst) {
1637  Warning("FindObject",
1638  "No bin %d %s histo in list %s\n",bin,nameH,lst->GetName());
1639  lst->Print();
1640  return 0;
1641  }
1642  // already normalized
1643  if (!normPerEvent || hst->TestBit(kBitNormPerEvent)) return hst;
1644 
1645  TString nameHS = nameH;
1646  // never scale stat. histo
1647  if (nameHS==kHStatName) return hst;
1648 
1649  //
1650  TH1* hstat = (TH1*)FindObject(-1,kHStatName,lst,kFALSE);
1651  // MyPrint(2,"Find statistics histogram: %p", hstat);
1652 
1653  //double nrm = hstat->GetBinContent(kBinEntries +
1654  // kEvProcInj+bin*kEntriesPerBin);
1655  double nrm = hstat->GetBinContent(kBinEntries + kEvProcData+
1656  bin*kEntriesPerBin); // HACK
1657  if (nrm<1) {
1658  Warning("FindObject",
1659  "%s: Anomaluous %d number of events processed in bin %d of list %p",
1660  buff,int(nrm),bin,lst);
1661  return 0;
1662  }
1663  MyPrint(2,"Normalization: %f", nrm);
1664 
1665  //
1666  if (hst->InheritsFrom(TH1::Class())) {
1667  MyPrint(0,"Scaling histogram %s by %f", hst->GetName(), nrm);
1668  ((TH1*)hst)->Scale(1./nrm);
1669  }
1670  else if (hst->InheritsFrom(THnSparse::Class())) {
1671  THnSparse* spr = (THnSparse*) hst;
1672  spr->Sumw2();
1673  int coord[3] = {0,0,0};
1674  for (Long64_t i = 0; i < spr->GetNbins(); ++i) {
1675  // Get the content of the bin from the current histogram
1676  Double_t v = spr->GetBinContent(i, coord);
1677  spr->SetBinContent(coord, v/nrm);
1678  spr->SetBinError(coord,TMath::Sqrt(v)/nrm);
1679  }
1680  }
1681  //
1682  MyPrint(2,"Flagging %s a scaled", hst->GetName());
1683  hst->SetBit(kBitNormPerEvent);
1684  return hst;
1685 }
1686 
1687 //______________________________________________________________________
1688 TList* LoadList(const char* flName, const char* addPref, const char* nameL)
1689 {
1690  MyPrint(2,"Loading lists from %s (prefix=%s, name=%s)", flName, addPref, nameL);
1691 
1692  // load list with histos
1693  TString nms = flName;
1694  gSystem->ExpandPathName(nms);
1695 
1696  TFile* fl = TFile::Open(nms.Data());
1697  if (!fl) {
1698  Error("LoadList","No file %s\n",nms.Data());
1699  return 0;
1700  }
1701 
1702  TList* lst = (TList*)fl->Get(nameL);
1703  if (!lst) {
1704  Error("LoadList","No list %s in file %s\n",nameL,nms.Data());
1705  return 0;
1706  }
1707 
1708  MyPrint(2,"Setting name to %s", flName);
1709  lst->SetName(flName);
1710  int nEnt = lst->GetSize();
1711  TString nm;
1712  for (int i=0;i<nEnt;i++) {
1713  TNamed* ob = (TNamed*)lst->At(i);
1714  nm = addPref;
1715  nm += ob->GetName();
1716  MyPrint(2," Renaming %40s to %s", ob->GetName(), nm.Data());
1717  ob->SetName(nm.Data());
1718  }
1719  //
1720  return lst;
1721 }
1722 
1723 //____________________________________________________________________________
1724 void GetRatE(double x,double xe, double y,double ye, double &rat, double &rate)
1725 {
1726  MyGuard(10,"Calculating ratio of %f +/- %f over %f +/- %f", x, xe, y, ye);
1727  rat = 0; rate = 0;
1728  if (TMath::Abs(y)<1e-16 || TMath::Abs(x)<1e-16) return;
1729  rat = x/y;
1730  rate = rat*TMath::Sqrt( xe*xe/(x*x) + ye*ye/(y*y));
1731  MyPrint(10,"Result: %f +/- %f", rat, rate);
1732 }
1733 
1734 //____________________________________________________________________________
1735 void Integrate(TH1* hist, double xmn,double xmx, double &val, double& err)
1736 {
1737  MyGuard(1,"Integrate %s over [%f,%f]", hist->GetName(), xmn, xmx);
1738  // integrate 1d histo within given limits
1739  TAxis* xax = hist->GetXaxis();
1740  int bmn = xax->FindBin(xmn+kEps); if (bmn<1) bmn = 0; // include
1741  int bmx = xax->FindBin(xmx-kEps);
1742  val = hist->IntegralAndError(bmn,bmx,err);
1743  // is this histo with symmetric axis ? then integrate also negative half axis
1744  if (TMath::Abs( xax->GetXmin() + xax->GetXmax() )<1e-6) {
1745  bmn = xax->FindBin(-xmx+kEps);
1746  bmx = xax->FindBin(-xmn-kEps);
1747  double errn;
1748  val += hist->IntegralAndError(bmn,bmx,errn);
1749  err = TMath::Sqrt(err*err + errn*errn);
1750  }
1751  MyPrint(1,"Integral: %f +/- %f", val, err);
1752 }
1753 
1754 
1755 //____________________________________________________________________________
1756 const char* HName(const char* prefix,const char* htype)
1757 {
1758  // compose the name of histo in the clist
1759  static TString strh;
1760  strh = "Tr"; strh += prefix; strh += "_"; strh += htype;
1761  return strh.Data();
1762 }
1763 
1764 //____________________________________________________________________________
1765 Int_t CheckStat(const TList* lst, const char* dtType)
1766 {
1767  MyGuard(3,"Check input list %s for %s", lst->GetName(), dtType);
1768  // check if needed bg was generated
1769  TH1* hstat = (TH1*)FindObject(-1,kHStatName,lst);
1770  TString dts = dtType;
1771  if (dts=="Data") return int( hstat->GetBinContent(kEvProcData) );
1772  if (dts=="Mix") return int( hstat->GetBinContent(kEvProcMix) );
1773  if (dts=="Inj") return int( hstat->GetBinContent(kEvProcInj) );
1774  if (dts=="Rot") return int( hstat->GetBinContent(kEvProcRot) );
1775  MyPrint(3,"Unknown process %s statistics is checked. Alowed: Data,Mix,Inj,Rot",dtType);
1776  return 0;
1777 }
1778 
1779 
1780 void GetRealMinMax(TH1* histo, double &vmn, double &vmx)
1781 {
1782  MyGuard(3,"Get real min and max of %s", histo->GetName());
1783  TAxis *xax = histo->GetXaxis();
1784  int nbx = xax->GetNbins();
1785  vmn=1e6, vmx=-1e6;
1786  if (histo->InheritsFrom(TH2::Class())) {
1787  TAxis *yax = histo->GetYaxis();
1788  int nby = yax->GetNbins();
1789  for (int ix=nbx+2;ix--;) {
1790  for (int iy=nby+2;iy--;) {
1791  double vl = histo->GetBinContent(ix,iy);
1792  if (vl<kEps) continue;
1793  if (vl<vmn) vmn = vl;
1794  if (vl>vmx) vmx = vl;
1795  }
1796  }
1797  }
1798  //
1799  else {
1800  for (int ix=nbx+2;ix--;) {
1801  double vl = histo->GetBinContent(ix);
1802  if (vl<vmn) vmn = vl;
1803  if (vl>vmx) vmx = vl;
1804  }
1805  }
1806  //
1807 }
1808 
1809 //____________________________________________________________________________
1810 void CorrectForZV(TH2* hEtaZ, TH1* hZv)
1811 {
1812  MyGuard(2,"Correcting %s for vertex efficiency using %s",
1813  hEtaZ->GetName(), hZv->GetName());
1814  // correct for the nonuniformity of Zvertices along Z
1815  int nbEta = hEtaZ->GetNbinsX();
1816  int nbZv = hEtaZ->GetNbinsY();
1817  //
1818  // printf("CorrectForZV: %s %s\n",hEtaZ->GetName(), hZv->GetName());
1819  TAxis* zax = hEtaZ->GetYaxis();
1820  for (int ibz=1;ibz<=nbZv;ibz++) {
1821  double z = zax->GetBinCenter(ibz);
1822  int izbh = hZv->FindBin(z);
1823  double scl = hZv->GetBinContent(izbh);
1824  double scle = hZv->GetBinError(izbh);
1825  MyPrint(2,"Vertex correction for z=%f: %f +/- %f", z, scl, scle);
1826  //scle = 0;
1827  if (scl<1e-9) {
1828  MyPrint(2,"Did not find ZV weight for Z=%f, stop",z);
1829  continue;
1830  }
1831  for (int ibe=1;ibe<=nbEta;ibe++) {
1832  double val = hEtaZ->GetBinContent(ibe,ibz);
1833  double err = hEtaZ->GetBinError(ibe,ibz);
1834  // err = 0;
1835  double rat = val/scl;
1836  double ratE = (val>0 && scl>0) ? rat*TMath::Sqrt((err*err)/(val*val) + (scle*scle)/(scl*scl) ) : 0.0;
1837  // double eta = hEtaZ->GetXaxis()->GetBinCenter(ibe);
1838  // printf("Zc: %+.5e Eta: %+.5e Val: %+.5e Scl:%+.5e ->%.5e\n",z,eta,val,scl,val/scl);
1839  hEtaZ->SetBinContent(ibe,ibz, rat);
1840  hEtaZ->SetBinError(ibe,ibz, ratE);
1841  }
1842  }
1843 }
1844 
1845 //____________________________________________________________________________
1846 TH1* ProjNorm(TH2* hEtaZ, TH1* hZv,const char* name, Int_t firstbin, Int_t lastbin)
1847 {
1848  MyGuard(1,"Projection of %s - %s", hEtaZ->GetName(), hEtaZ->GetTitle());
1849  // correct for the nonuniformity of Zvertices along Z
1850  int nbEta = hEtaZ->GetNbinsX();
1851  int nbZv = hEtaZ->GetNbinsY();
1852  if (firstbin<1) firstbin = 1;
1853  if (lastbin<firstbin) lastbin = nbZv;
1854  //
1855  TH1* prj = hEtaZ->ProjectionX(name,firstbin,lastbin,"e");
1856  prj->Reset();
1857  MyPrint(2,"Made projection and reset it");
1858  double *arrV = new double[nbZv];
1859  double *arrE = new double[nbZv];
1860  double *arZV = new double[nbZv];
1861  double *arZE = new double[nbZv];
1862  double *relE = new double[nbZv];
1863  int *ind = new int[nbZv];
1864  //
1865  TAxis* zax = hEtaZ->GetYaxis();
1866  for (int ibe=1;ibe<=nbEta;ibe++) {
1867  //
1868  int cnt = 0;
1869  for (int ibz=firstbin;ibz<=lastbin;ibz++) {
1870  double val = hEtaZ->GetBinContent(ibe,ibz);
1871  double err = hEtaZ->GetBinError(ibe,ibz);
1872  if (val<1e-9) continue;
1873  //
1874  double z = zax->GetBinCenter(ibz);
1875  int izbh = hZv->FindBin(z);
1876  //
1877  arrV[cnt] = val;
1878  arrE[cnt] = err;
1879  relE[cnt] = err/val;
1880  arZV[cnt] = hZv->GetBinContent(izbh);
1881  arZE[cnt] = hZv->GetBinError(izbh);
1882  cnt++;
1883  }
1884  // suppress bins increasing the error
1885  TMath::Sort(cnt,relE,ind,kFALSE);
1886  double av=0,ave=0,vsm=0,vsme=0,zsm=0,zsme=0;
1887  double rat = 1e99;
1888  int iu=0;
1889  for (iu=0;iu<cnt;iu++) {
1890  int j = ind[iu];
1891  double rt = TMath::Sqrt(vsme+arrE[j]*arrE[j])/(vsm+ arrV[j]);
1892  // printf("%2d(%2d/%2d) New:%.e Old:%.e | %e %e %e\n",ibe,j,iu,rt,rat,arrV[j],arrE[j],arrE[j]/arrV[j]);
1893  MyPrint(10," %10s - bin %3d,%3d %9.3f+/-%9.3f (%5.2f%%) [%9.3f+/-%9.3f] "
1894  "-> %9.3f (%9.3g)",
1895  hEtaZ->GetName(), ibe, j, arrV[j], arrE[j],
1896  100*relE[j], arZV[j], arZE[j], rt, rat);
1897  if (rt>rat) continue;//break;
1898  rat = rt;
1899  vsm += arrV[j];
1900  vsme+= arrE[j]*arrE[j];
1901  zsm += arZV[j];
1902  zsme+= arZE[j]*arZE[j];
1903  // printf("%2d(%2d/%2d) New:%.e Old:%.e | %e %e %e | %e %e | %e %e\n",ibe,j,iu,rt,rat,arrV[j],arrE[j],arrE[j]/arrV[j],
1904  // vsm,TMath::Sqrt(vsme), zsm, TMath::Sqrt(zsme) );
1905  //
1906  }
1907  if (iu==0) continue; // no contribution
1908  av = vsm/zsm;
1909  zsme /= zsm*zsm;
1910  vsme /= vsm*vsm;
1911  ave = av*TMath::Sqrt( vsme + zsme );
1912  MyPrint(10,
1913  "%10s - bin %3d (%+6.3f) count=%2d n=%9.3f+/-%9.3f d=%9.3f+/-%9.3f"
1914  " -> %9.3f +/- %9.3f", hEtaZ->GetName(), ibe,
1915  hEtaZ->GetXaxis()->GetBinCenter(ibe), iu, vsm, TMath::Sqrt(vsme),
1916  zsm, TMath::Sqrt(zsme), av, ave);
1917  // printf("->%d %.4e(%.4e) %2d bins out of %2d\n",ibe, av,ave,iu,cnt);
1918  prj->SetBinContent(ibe, av);
1919  prj->SetBinError(ibe, ave);
1920  }
1921  //
1922  delete[] arrV;
1923  delete[] arrE;
1924  delete[] arZV;
1925  delete[] arZE;
1926  delete[] relE;
1927  delete[] ind;
1928  //
1929  return prj;
1930 }
1931 
1932 //____________________________________________________________________________
1933 TH1* ProjectWghMean(TH2* hEtaZ, const char* name, Int_t firstbin, Int_t lastbin, double rejOutliers)
1934 {
1935  MyGuard(2,"Calculating weighted mean projection of %s from %d to %d",
1936  hEtaZ, firstbin, lastbin);
1937  // for each X bin calculated the weighted average over Y bins
1938  int nbEta = hEtaZ->GetNbinsX();
1939  int nbZv = hEtaZ->GetNbinsY();
1940  if (firstbin<1) firstbin = 1;
1941  if (lastbin<firstbin) lastbin = nbZv;
1942  //
1943  MyPrint(2,"HISTO: %s %s %p | bins: %d %d",hEtaZ->GetName(),hEtaZ->GetTitle(), hEtaZ, firstbin, lastbin);
1944  TArrayD val(nbZv),err(nbZv),prc(nbZv);
1945  TArrayI ind(nbZv);
1946  //
1947  TH1* prj = hEtaZ->ProjectionX(name,firstbin,lastbin,"e");
1948  prj->Reset();
1949  MyPrint(2,"Made projection, and reset it");
1950  for (int ibe=1;ibe<=nbEta;ibe++) {
1951  int cnt = 0;
1952  double valAv = 0;
1953  double errAv = 0;
1954  //
1955  for (int ibz=firstbin;ibz<=lastbin;ibz++) {
1956  val[cnt] = hEtaZ->GetBinContent(ibe,ibz);
1957  err[cnt] = hEtaZ->GetBinError(ibe,ibz);
1958  if (!err[cnt]>0) continue;
1959  valAv += val[cnt]/(err[cnt]*err[cnt]);
1960  errAv += 1./(err[cnt]*err[cnt]);
1961  // if (errAv>0) printf("%2d| (%+e %+e) %+e %+e\n",cnt,val[cnt],err[cnt],valAv/errAv,1./sqrt(errAv));
1962  cnt++;
1963  }
1964  if (cnt<1) continue;
1965  valAv /= errAv;
1966  errAv = 1./TMath::Sqrt(errAv);
1967  // printf("#%3d Full W.Mean for %3d bins: %+e +- %e\n",ibe,cnt,valAv,errAv);
1968  //
1969  int cntSkip=-1;
1970  int useCnt = 0;
1971  double valFin=0,errFin=0;
1972  while(cntSkip<cnt-1) {
1973  // get truncated mean
1974  TMath::Sort(cnt,val.GetArray(),ind.GetArray(), kFALSE); // sort in increasing value
1975  if (cntSkip>-1) err[ind[cntSkip]] = -1;
1976  valAv = 0;
1977  errAv = 0;
1978  for (int ic=0;ic<cnt;ic++) {
1979  double erb = err[ind[ic]];
1980  double vlb = val[ind[ic]];
1981  if (!(erb>0)) continue;
1982  valAv += vlb/(erb*erb);
1983  errAv += 1./(erb*erb);
1984  useCnt++;
1985  }
1986  if (useCnt<1 || !(errAv>0)) break;
1987  valAv /= errAv;
1988  errAv = 1./TMath::Sqrt(errAv);
1989  // printf("#%3d Truncated0 W.Mean for %d out of %d bins: %+e +- %e\n",ibe,useCnt,cnt,valAv,errAv);
1990  //
1991  for (int ic=0;ic<cnt;ic++) prc[ic] = err[ic]>0 ? TMath::Abs(valAv-val[ic])/err[ic] : 1e9;
1992  // take most precize values
1993  TMath::Sort(cnt,prc.GetArray(),ind.GetArray(), kFALSE); // sort in increasing error
1994  useCnt = cnt>10 ? cnt/2 : cnt*0.7;
1995  if (useCnt<1) useCnt = 1;
1996  valAv = 0;
1997  errAv = 0;
1998  for (int ic=0;ic<useCnt;ic++) {
1999  double erb = err[ind[ic]];
2000  double vlb = val[ind[ic]];
2001  if (!(erb>0)) continue;
2002  valAv += vlb/(erb*erb);
2003  errAv += 1./(erb*erb);
2004  }
2005  if (!(errAv>0)) break;
2006  valAv /= errAv;
2007  errAv = 1./TMath::Sqrt(errAv);
2008  // printf("#%3d Truncated W.Mean for %d out of %d bins: %+e +- %e\n",ibe,useCnt,cnt,valAv,errAv);
2009  valFin=0;
2010  errFin=0;
2011  useCnt = 0;
2012  for (int ic=0;ic<cnt;ic++) {
2013  double erb = err[ic];
2014  double vlb = val[ic];
2015  if (!(erb>0)) continue;
2016  double dev = TMath::Abs(valAv-vlb)/erb;
2017  if (rejOutliers>0 && dev>rejOutliers) continue;
2018  //
2019  valFin += vlb/(erb*erb);
2020  errFin += 1./(erb*erb);
2021  useCnt++;
2022  }
2023  if (useCnt<1 || !(errFin)>0) {
2024  cntSkip++;
2025  continue;
2026  }
2027  valFin /= errFin;
2028  errFin = 1./TMath::Sqrt(errFin);
2029  break;
2030  }
2031  // estimate error
2032  if (!(errFin>0)) {
2033  printf("failed to evaluate eta bin %d\n",ibe);
2034  continue;
2035  }
2036  //
2037  // printf("#%3d Truncated W.Mean for %d out of %d selected bins: %+e +- %e\n",ibe,useCnt,cnt,valFin,errFin);
2038  //
2039  prj->SetBinContent(ibe,valFin);
2040  prj->SetBinError(ibe,errFin);
2041  }
2042  //
2043  return prj;
2044 }
2045 
2046 void KillBadBins(TH2* histo, double mn,double mx)
2047 {
2048  MyGuard(2,"Killing bad (<%f or >%f) bins of %s",
2049  mn, mx, histo->GetName());
2050  int nbx = histo->GetNbinsX();
2051  int nby = histo->GetNbinsY();
2052  for (int ix=1;ix<=nbx;ix++) {
2053  for (int iy=1;iy<=nby;iy++) {
2054  double vl = histo->GetBinContent(ix,iy);
2055  if (vl>mx || vl<mn) {
2056  histo->SetBinContent(ix,iy,0);
2057  histo->SetBinError(ix,iy,0);
2058  }
2059  }
2060  }
2061  histo->SetMinimum(mn);
2062  histo->SetMaximum(mx);
2063 }
2064 
2065 void PrintH(TH2* h, Int_t prec)
2066 {
2067  Printf("Content of %s - %s", h->GetName(), h->GetTitle());
2068  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
2069  printf("%3d: ", i);
2070  for (Int_t j = 1; j <= h->GetNbinsY(); j++) {
2071  printf("%.*f+/-%.*f ",
2072  prec, h->GetBinContent(i,j),
2073  prec, h->GetBinError(i,j));
2074  }
2075  printf("\n");
2076  }
2077 }
2078 //____________________________________________________________________
2079 void PrintH(TH1* h, Int_t prec)
2080 {
2081  Printf("Content of %s - %s", h->GetName(), h->GetTitle());
2082  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
2083  if (h->GetBinContent(i) <= 1e-6) continue;
2084  printf("%3d (%+6.3f): %.*f+/-%.*f\n", i,
2085  h->GetXaxis()->GetBinCenter(i),
2086  prec, h->GetBinContent(i),
2087  prec, h->GetBinError(i));
2088  }
2089 }
kSclIntegral
Definition: CorrectSpectraMultiMCBG.C:43
kDTht
Definition: CorrectSpectraMultiMCBG.C:74
kSclTypes
Definition: CorrectSpectraMultiMCBG.C:44
kNormShapeDPhi
Definition: CorrectSpectraMultiMCBG.C:39
kHStatName
const char kHStatName[]
Definition: CorrectSpectraMultiMCBG.C:26
Decr
void Decr()
Definition: CorrectSpectraMultiMCBG.C:125
kWDistSgCut
double kWDistSgCut
Definition: CorrectSpectraMultiMCBG.C:36
outTitle
char outTitle[1000]
Definition: CorrectSpectraMultiMCBG.C:188
Double_t
double Double_t
Definition: External.C:58
HName
const char * HName(const char *prefix, const char *htype)
Definition: CorrectSpectraMultiMCBG.C:1756
MyGuard
#define MyGuard(L, F,...)
Definition: CorrectSpectraMultiMCBG.C:146
kMCV0Scale
Definition: CorrectSpectraMultiMCBG.C:93
AddLabel
TLatex * AddLabel(const char *txt, float x=0.1, float y=0.9, int color=kBlack, float size=0.04)
Definition: SaveCanvas.C:61
TH2F
Definition: External.C:236
creatSpeciesCMacro
Bool_t creatSpeciesCMacro
Definition: CorrectSpectraMultiMCBG.C:196
title
const char * title
Definition: MakeQAPdf.C:26
SaveCanvas
void SaveCanvas(TCanvas *canv, const char *path="canv", const Option_t *option="ecn")
Definition: SaveCanvas.C:28
kBinEntries
Definition: CorrectSpectraMultiMCBG.C:96
kNWorkers
Definition: CorrectSpectraMultiMCBG.C:70
kPhiRot
Definition: CorrectSpectraMultiMCBG.C:83
k1MBetaMC
Definition: CorrectSpectraMultiMCBG.C:106
useScaleType
Int_t useScaleType
Definition: CorrectSpectraMultiMCBG.C:55
maxAlpha
Double_t maxAlpha
Definition: CorrectSpectraMultiMCBG.C:62
CropHisto
void CropHisto(TH1 *histo, int b00, int b01, int b10=-1, int b11=-1)
Definition: CorrectSpectraMultiMCBG.C:1463
kNormShapes
Definition: CorrectSpectraMultiMCBG.C:40
kZEtaOvl
Definition: CorrectSpectraMultiMCBG.C:80
PrepareHistos
Bool_t PrepareHistos(int bin, TList *lst, TList *lisMC, Bool_t isMC)
Definition: CorrectSpectraMultiMCBG.C:305
kNHistos
Definition: CorrectSpectraMultiMCBG.C:107
kDummyBin
Definition: CorrectSpectraMultiMCBG.C:66
Long64_t
long long Long64_t
Definition: External.C:43
fmt
const char * fmt
Definition: BackOfTheEnvelope.C:41
GetRealMinMax
void GetRealMinMax(TH1 *h, double &vmn, double &vmx)
Definition: CorrectSpectraMultiMCBG.C:1780
kThtCW
Definition: CorrectSpectraMultiMCBG.C:78
kZvDistCorr
Definition: CorrectSpectraMultiMCBG.C:107
gSystem
TSystem * gSystem
SaveCanvas.h
use1mBeta
Bool_t use1mBeta
Definition: CorrectSpectraMultiMCBG.C:49
kDataDist
Definition: CorrectSpectraMultiMCBG.C:107
kSignalEst
Definition: CorrectSpectraMultiMCBG.C:106
kSigCorr
Definition: CorrectSpectraMultiMCBG.C:106
kBgRescFc
Definition: CorrectSpectraMultiMCBG.C:106
kEvTot0
Definition: CorrectSpectraMultiMCBG.C:67
c
TCanvas * c
Definition: TestFitELoss.C:172
kdPhiBgTailMin
double kdPhiBgTailMin
Definition: CorrectSpectraMultiMCBG.C:29
KillBadBins
void KillBadBins(TH2 *histo, double mn=-1e50, double mx=1e50)
Definition: CorrectSpectraMultiMCBG.C:2046
nCentBins
Int_t nCentBins
Definition: CorrectSpectraMultiMCBG.C:184
kThtS2
Definition: CorrectSpectraMultiMCBG.C:77
useShapeType
Int_t useShapeType
Definition: CorrectSpectraMultiMCBG.C:51
kEps
double kEps
Definition: CorrectSpectraMultiMCBG.C:27
kZVMin
Definition: CorrectSpectraMultiMCBG.C:87
PrintH
void PrintH(TH2 *h, Int_t prec=2)
Definition: CorrectSpectraMultiMCBG.C:2065
ProjectWghMean
TH1 * ProjectWghMean(TH2 *hEtaZ, const char *name="_px", Int_t firstbin=0, Int_t lastbin=-1, double rejOutliers=6.)
Definition: CorrectSpectraMultiMCBG.C:1933
figDir
const char * figDir
Definition: CorrectSpectraMultiMCBG.C:47
PrintAndPause
void PrintAndPause(TCanvas *c, const TString &what, Bool_t wait)
Definition: CorrectSpectraMultiMCBG.C:149
kZvDist
Definition: CorrectSpectraMultiMCBG.C:107
_MyGuard::_MyGuard
_MyGuard(UShort_t lvl, const char *fmt,...)
Definition: CorrectSpectraMultiMCBG.C:128
kPhiShift
Definition: CorrectSpectraMultiMCBG.C:76
dNdEta
TArrayD dNdEta
Definition: CorrectSpectraMultiMCBG.C:191
kDPiSCut
Definition: CorrectSpectraMultiMCBG.C:90
kCentVar
Definition: CorrectSpectraMultiMCBG.C:72
kEntriesPerBin
Definition: CorrectSpectraMultiMCBG.C:102
Int_t
int Int_t
Definition: External.C:63
kEvInMltBin
Definition: CorrectSpectraMultiMCBG.C:97
_MyGuard::ok
Bool_t ok
Definition: CorrectSpectraMultiMCBG.C:143
CopyAdd
TH1 * CopyAdd(TH1 *h, const char *name, const char *title, TObjArray *col, Int_t location, Int_t shift)
Definition: CorrectSpectraMultiMCBG.C:291
kInjScl
Definition: CorrectSpectraMultiMCBG.C:84
kWDistBgTailMin
double kWDistBgTailMin
Definition: CorrectSpectraMultiMCBG.C:32
kBitNormPerEvent
Definition: CorrectSpectraMultiMCBG.C:64
kSclWghMean
Definition: CorrectSpectraMultiMCBG.C:42
kAlphaMC
Definition: CorrectSpectraMultiMCBG.C:106
_MyPrint
void _MyPrint(UShort_t lvl, const char *fmt,...)
Definition: CorrectSpectraMultiMCBG.C:112
NormalizeBg
TH1 * NormalizeBg(TH1 *dataH, TH1 *bgH, double &scl, double &scle)
Definition: CorrectSpectraMultiMCBG.C:1527
kMCShift
Definition: CorrectSpectraMultiMCBG.C:107
kPhiOvl
Definition: CorrectSpectraMultiMCBG.C:79
kOneUnit
Definition: CorrectSpectraMultiMCBG.C:69
kBgEst
Definition: CorrectSpectraMultiMCBG.C:106
minAlpha
Double_t minAlpha
Definition: CorrectSpectraMultiMCBG.C:61
kNStd
Definition: CorrectSpectraMultiMCBG.C:75
CorrectSpectraMultiMCBG
void CorrectSpectraMultiMCBG(const char *flNameData, const char *flNameMC, const char *unique="", int maxBins=10, Bool_t waitForUser=false, const char *bgType="Comb")
Definition: CorrectSpectraMultiMCBG.C:198
kBgDist
Definition: CorrectSpectraMultiMCBG.C:107
_MyGuard::~_MyGuard
~_MyGuard()
Definition: CorrectSpectraMultiMCBG.C:142
PlotResults
void PlotResults(Bool_t waitForUser)
Definition: CorrectSpectraMultiMCBG.C:760
kNoOvl
Definition: CorrectSpectraMultiMCBG.C:81
kEvProcInj
Definition: CorrectSpectraMultiMCBG.C:99
kNStdCut
Definition: CorrectSpectraMultiMCBG.C:91
kBgMC
Definition: CorrectSpectraMultiMCBG.C:106
kEvProcMix
Definition: CorrectSpectraMultiMCBG.C:101
Incr
void Incr()
Definition: CorrectSpectraMultiMCBG.C:124
resArr
TObjArray resArr
Definition: CorrectSpectraMultiMCBG.C:186
kEvTot
Definition: CorrectSpectraMultiMCBG.C:68
kZvMCDistNS
Definition: CorrectSpectraMultiMCBG.C:107
Integrate
void Integrate(TH1 *hist, double xmn, double xmx, double &val, double &err)
Definition: CorrectSpectraMultiMCBG.C:1735
ProjNorm
TH1 * ProjNorm(TH2 *hEtaZ, TH1 *hZv, const char *name="_px", Int_t firstbin=0, Int_t lastbin=-1)
Definition: CorrectSpectraMultiMCBG.C:1846
PlotDNDEta
void PlotDNDEta(int bin)
Definition: CorrectSpectraMultiMCBG.C:853
PlotSpecies
void PlotSpecies()
Definition: CorrectSpectraMultiMCBG.C:1340
PlotAlphaBeta
void PlotAlphaBeta(int bin)
Definition: CorrectSpectraMultiMCBG.C:1241
kWDistBgTailMax
double kWDistBgTailMax
Definition: CorrectSpectraMultiMCBG.C:33
creatAlphaBetaCMacro
Bool_t creatAlphaBetaCMacro
Definition: CorrectSpectraMultiMCBG.C:195
FindObject
TObject * FindObject(int bin, const char *nameH, const TList *lst, Bool_t normPerEvent=kTRUE)
Definition: CorrectSpectraMultiMCBG.C:1615
outStr
char outStr[1000]
Definition: CorrectSpectraMultiMCBG.C:187
normToMB
Bool_t normToMB
Definition: CorrectSpectraMultiMCBG.C:58
resDir
const char * resDir
Definition: CorrectSpectraMultiMCBG.C:48
kSignalEstMC
Definition: CorrectSpectraMultiMCBG.C:106
kRawDtCut
Definition: CorrectSpectraMultiMCBG.C:106
dNdEtaErr
TArrayD dNdEtaErr
Definition: CorrectSpectraMultiMCBG.C:191
resDnDeta
TObjArray resDnDeta
Definition: CorrectSpectraMultiMCBG.C:186
kEtaMax
Definition: CorrectSpectraMultiMCBG.C:86
LoadList
TList * LoadList(const char *flName, const char *addPref, const char *nameL="clist")
Definition: CorrectSpectraMultiMCBG.C:1688
creatDnDEtaCMacro
Bool_t creatDnDEtaCMacro
Definition: CorrectSpectraMultiMCBG.C:194
grp
GRPData * grp
Definition: GRP.C:361
isMC
Bool_t isMC
Definition: ReadDvsMultiplicity.C:60
TH2
Definition: External.C:220
scaleBG
Double_t scaleBG
Definition: CorrectSpectraMultiMCBG.C:54
CorrectForZV
void CorrectForZV(TH2 *hEtaZ, TH1 *hZv)
Definition: CorrectSpectraMultiMCBG.C:1810
canvFin
TCanvas * canvFin
Definition: CorrectSpectraMultiMCBG.C:192
kNormShapeDist
Definition: CorrectSpectraMultiMCBG.C:38
kAlpha
Definition: CorrectSpectraMultiMCBG.C:106
kMCPrim
Definition: CorrectSpectraMultiMCBG.C:106
kZvEff
Definition: CorrectSpectraMultiMCBG.C:107
k1MBeta
Definition: CorrectSpectraMultiMCBG.C:106
MyPrint
#define MyPrint(L, F,...)
Definition: CorrectSpectraMultiMCBG.C:145
UShort_t
unsigned short UShort_t
Definition: External.C:28
uniqueName
TString uniqueName
Definition: CorrectSpectraMultiMCBG.C:189
k1MBetaMCscl
Definition: CorrectSpectraMultiMCBG.C:106
useZbinWAv
Bool_t useZbinWAv
Definition: CorrectSpectraMultiMCBG.C:56
useMCLB
Int_t useMCLB
Definition: CorrectSpectraMultiMCBG.C:53
listMC
TList * listMC
Definition: CorrectSpectraMultiMCBG.C:185
useBgType
TString useBgType
Definition: CorrectSpectraMultiMCBG.C:50
kDPhi
Definition: CorrectSpectraMultiMCBG.C:73
kEvProcData
Definition: CorrectSpectraMultiMCBG.C:98
fgDebug
UShort_t fgDebug
Definition: CorrectSpectraMultiMCBG.C:111
_MyGuard
Definition: CorrectSpectraMultiMCBG.C:126
Bool_t
bool Bool_t
Definition: External.C:53
kEvProcRot
Definition: CorrectSpectraMultiMCBG.C:100
kEtaFitRange
const double kEtaFitRange
Definition: CorrectSpectraMultiMCBG.C:59
kZVMax
Definition: CorrectSpectraMultiMCBG.C:88
kdPhiSgCut
double kdPhiSgCut
Definition: CorrectSpectraMultiMCBG.C:35
kdPhiBgTailMax
double kdPhiBgTailMax
Definition: CorrectSpectraMultiMCBG.C:30
GetRatE
void GetRatE(double x, double xe, double y, double ye, double &rat, double &rate)
Definition: CorrectSpectraMultiMCBG.C:1724
CheckStat
Int_t CheckStat(const TList *lst, const char *dtType)
Definition: CorrectSpectraMultiMCBG.C:1765
kBgMCDist
Definition: CorrectSpectraMultiMCBG.C:107
TH1
Definition: External.C:196
kEtaMin
Definition: CorrectSpectraMultiMCBG.C:85
SetHStyle
void SetHStyle(TH1 *hst, int col=kRed, int mark=20, float mrsize=0.7)
Definition: SaveCanvas.C:130
listDt
TList * listDt
Definition: CorrectSpectraMultiMCBG.C:185
ProcessHistos
void ProcessHistos(int bin)
Definition: CorrectSpectraMultiMCBG.C:633