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