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