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