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