]>
Commit | Line | Data |
---|---|---|
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 | ||
26 | const char kHStatName[]="hStat"; | |
27 | double kEps = 1e-6; | |
28 | // | |
29 | double kdPhiBgTailMin = 0.1; // lower limit of dphi tail to use for bg normalization | |
30 | double kdPhiBgTailMax = 0.3; // upper limit of dphi tail to use for bg normalization | |
31 | // | |
32 | double kWDistBgTailMin = 5.; // lower limit of wgh.distance tail to use for bg normalization | |
33 | double kWDistBgTailMax = 25.; // upper limit of wgh.distance tail to use for bg normalization | |
34 | ||
35 | double kdPhiSgCut=-1; // cut in dphi-bent used to extract the signal, extracted from stat histo | |
36 | double kWDistSgCut=-1; // cut in w.distance used to extract the signal, extracted from stat histo | |
37 | // | |
38 | enum { kNormShapeDist, // normalize bg tails usig weighted distance shape | |
39 | kNormShapeDPhi, // normalize bg tails usig dPhi-bend shape | |
40 | kNormShapes}; | |
41 | ||
42 | enum { kSclWghMean, // normalize bg tails to data using weighted mean of bin-by-bin ratios | |
43 | kSclIntegral, // normalize bg tails to data using integral | |
44 | kSclTypes}; | |
45 | ||
46 | ||
47 | const char* figDir = "figMult"; | |
69d0af30 | 48 | const char* resDir = "resMult"; |
a9a39f46 | 49 | TString useBgType = "Inj"; |
50 | Int_t useShapeType = kNormShapeDist; // which distribution to use for bg normalization | |
6c56e693 | 51 | Bool_t useMCLB = 0;//kFALSE; // use Comb MC Labels as a template for Bg. |
a9a39f46 | 52 | Int_t useScaleType = kSclIntegral;//kSclWghMean; // which type of tails normalization to use |
53 | const double kEtaFitRange = 0.5; | |
54 | ||
55 | enum {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 | ||
96 | enum {kSigCorr,kMCPrim,kRawDtCut,kSignalEst,kSignalEstMC,kBgEst,k1MBeta,k1MBetaMC,kAlpha,kAlphaMC,kBgMC,kBgRescFc,kDataDist,kBgDist,kBgMCDist, kMCShift=20, kNHistos=kMCShift+kMCShift}; | |
97 | ||
69d0af30 | 98 | void CorrectSpectraMulti(const char* flNameData, const char* flNameMC,const char* unique="",int maxBins=6); |
a9a39f46 | 99 | Bool_t PrepareHistos(int bin, TList* lst, Bool_t isMC); |
100 | void ProcessHistos(int bin); | |
101 | TH1* NormalizeBg(TH1* dataH, TH1* bgH, double &scl, double &scle); | |
102 | TObject* FindObject(int bin, const char* nameH, const TList* lst, Bool_t normPerEvent=kTRUE); | |
103 | TList* LoadList(const char* flName, const char* addPref, const char* nameL="clist"); | |
104 | void GetRatE(double x,double xe, double y,double ye, double &rat, double &rate); | |
105 | Int_t CheckStat(const TList* lst,const char* dtType); | |
106 | void Integrate(TH1* hist, double xmn,double xmx, double &val, double& err); | |
107 | void CropHisto(TH1* histo, int b00, int b01, int b10=-1, int b11=-1); | |
108 | void CropHisto(TH1* histo, double b00, double b01, double b10=-1, double b11=-1); | |
109 | void GetRealMinMax(TH1* h, double &vmn, double &vmx); | |
110 | const char* HName(const char* prefix,const char* htype); | |
111 | ||
112 | void PlotResults(); | |
113 | void PlotDNDEta(int bin); | |
114 | void PlotAlphaBeta(int bin); | |
115 | void PlotSpecies(); | |
116 | ||
117 | Float_t myMergeFactor = -1; // if the files were manually merged, scale everything except statistics by 1/myMergeFactor | |
118 | Int_t nCentBins = -1; | |
119 | TList *listDt=0, *listMC=0; | |
db1bbbe1 | 120 | TObjArray resArr, resDnDeta; |
a9a39f46 | 121 | char outStr[1000]; |
122 | char outTitle[1000]; | |
123 | TString uniqueName=""; | |
124 | // | |
125 | TArrayD dNdEta,dNdEtaErr; | |
126 | TCanvas *canvFin=0; | |
127 | // | |
128 | Bool_t creatDnDEtaCMacro = kFALSE; | |
129 | Bool_t creatAlphaBetaCMacro = kFALSE; | |
130 | Bool_t creatSpeciesCMacro = kFALSE; | |
131 | ||
69d0af30 | 132 | void 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 | //_____________________________________________________________________ | |
198 | Bool_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 | //_____________________________________________________________________ | |
396 | void 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 | ||
449 | void 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 | ||
534 | void 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 | // | |
815 | void 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 | ||
911 | void 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 | //______________________________________________________________________ | |
1032 | void 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 | //______________________________________________________________________ | |
1078 | void 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 | //______________________________________________________________________ | |
1094 | TH1* 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 | //______________________________________________________________________ | |
1169 | TObject* 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 | //______________________________________________________________________ | |
1213 | TList* 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 | //____________________________________________________________________________ | |
1237 | void 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 | //____________________________________________________________________________ | |
1246 | void 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 | //____________________________________________________________________________ | |
1265 | const 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 | //____________________________________________________________________________ | |
1274 | Int_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 | ||
1288 | void 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 | } |