]>
Commit | Line | Data |
---|---|---|
31d3e4f0 | 1 | #include "AliAnaChargedJetResponseMaker.h" |
2 | #include "TGraph.h" | |
3 | #include "TGraphErrors.h" | |
4 | #include "TMath.h" | |
5 | #include "Riostream.h" | |
6 | #include "TH1.h" | |
7 | #include "TRandom.h" | |
8 | #include "TFile.h" | |
9 | #include "TCanvas.h" | |
10 | #include "TF1.h" | |
11 | #include "THnSparse.h" | |
5d87a047 | 12 | #include "TArrayD.h" |
31d3e4f0 | 13 | |
14 | #define round(x) ((x)>=0?(long)((x)+0.5):(long)((x)-0.5)) | |
15 | ||
c64cb1f6 | 16 | using std::cout; |
17 | using std::endl; | |
18 | ||
31d3e4f0 | 19 | ClassImp(AliAnaChargedJetResponseMaker) |
20 | ||
21 | AliAnaChargedJetResponseMaker::AliAnaChargedJetResponseMaker(): | |
22 | fDebug(kFALSE), | |
23 | fResolutionType(kParam), | |
24 | fDeltaPt(0x0), | |
25 | fhDeltaPt(0x0), | |
26 | fDimensions(1), | |
27 | fDimRec(0), | |
28 | fDimGen(1), | |
29 | fPtMin(-999), | |
30 | fPtMax(-999), | |
31 | fNbins(0), | |
32 | fBinArrayPtRec(0x0), | |
33 | fPtMeasured(0x0), | |
34 | fEffFlat(0), | |
35 | fEfficiency(0x0), | |
36 | fEfficiencyFine(0x0), | |
37 | fResponseMatrix(0x0), | |
38 | fResponseMatrixFine(0x0), | |
39 | fPtMinUnfolded(0.), | |
40 | fPtMaxUnfolded(0.), | |
41 | fPtMaxUnfoldedHigh(-1.), | |
42 | fBinWidthFactorUnfolded(2), | |
43 | fSkipBinsUnfolded(0), | |
44 | fExtraBinsUnfolded(5), | |
45 | fbVariableBinning(kFALSE), | |
46 | fPtMaxUnfVarBinning(0), | |
47 | f1MergeFunction(0x0), | |
48 | fFineFrac(10), | |
49 | fbCalcErrors(kFALSE) | |
50 | {;} | |
51 | ||
31d3e4f0 | 52 | |
fa7c34ba | 53 | //-------------------------------------------------------------------------------------------------------------------------------------------------- |
54 | AliAnaChargedJetResponseMaker::AliAnaChargedJetResponseMaker(const AliAnaChargedJetResponseMaker& obj): | |
55 | fDebug(obj.fDebug), | |
56 | fResolutionType(obj.fResolutionType), | |
57 | fDeltaPt(obj.fDeltaPt), | |
58 | fhDeltaPt(obj.fhDeltaPt), | |
59 | fDimensions(obj.fDimensions), | |
60 | fDimRec(obj.fDimRec), | |
61 | fDimGen(obj.fDimGen), | |
62 | fPtMin(obj.fPtMin), | |
63 | fPtMax(obj.fPtMax), | |
64 | fNbins(obj.fNbins), | |
65 | fBinArrayPtRec(obj.fBinArrayPtRec), | |
66 | fPtMeasured(obj.fPtMeasured), | |
67 | fEffFlat(obj.fEffFlat), | |
68 | fEfficiency(obj.fEfficiency), | |
69 | fEfficiencyFine(obj.fEfficiencyFine), | |
70 | fResponseMatrix(obj.fResponseMatrix), | |
71 | fResponseMatrixFine(obj.fResponseMatrixFine), | |
72 | fPtMinUnfolded(obj.fPtMinUnfolded), | |
73 | fPtMaxUnfolded(obj.fPtMaxUnfolded), | |
74 | fPtMaxUnfoldedHigh(obj.fPtMaxUnfoldedHigh), | |
75 | fBinWidthFactorUnfolded(obj.fBinWidthFactorUnfolded), | |
76 | fSkipBinsUnfolded(obj.fSkipBinsUnfolded), | |
77 | fExtraBinsUnfolded(obj.fExtraBinsUnfolded), | |
78 | fbVariableBinning(obj.fbVariableBinning), | |
79 | fPtMaxUnfVarBinning(obj.fPtMaxUnfVarBinning), | |
80 | f1MergeFunction(obj.f1MergeFunction), | |
81 | fFineFrac(obj.fFineFrac), | |
82 | fbCalcErrors(obj.fbCalcErrors) | |
83 | {;} | |
84 | ||
85 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
86 | AliAnaChargedJetResponseMaker& AliAnaChargedJetResponseMaker::operator=(const AliAnaChargedJetResponseMaker& other) | |
87 | { | |
88 | // Assignment | |
89 | if(&other == this) return *this; | |
90 | AliAnaChargedJetResponseMaker::operator=(other); | |
91 | fDebug = other.fDebug; | |
92 | fResolutionType = other.fResolutionType; | |
93 | fDeltaPt = other.fDeltaPt; | |
94 | fhDeltaPt = other.fhDeltaPt; | |
95 | fDimensions = other.fDimensions; | |
96 | fDimRec = other.fDimRec; | |
97 | fDimGen = other.fDimGen; | |
98 | fPtMin = other.fPtMin; | |
99 | fPtMax = other.fPtMax; | |
100 | fNbins = other.fNbins; | |
101 | fBinArrayPtRec = other.fBinArrayPtRec; | |
102 | fPtMeasured = other.fPtMeasured; | |
103 | fEffFlat = other.fEffFlat; | |
104 | fEfficiency = other.fEfficiency; | |
105 | fEfficiencyFine = other.fEfficiencyFine; | |
106 | fResponseMatrix = other.fResponseMatrix; | |
107 | fResponseMatrixFine = other.fResponseMatrixFine; | |
108 | fPtMinUnfolded = other.fPtMinUnfolded; | |
109 | fPtMaxUnfolded = other.fPtMaxUnfolded; | |
110 | fPtMaxUnfoldedHigh = other.fPtMaxUnfoldedHigh; | |
111 | fBinWidthFactorUnfolded = other.fBinWidthFactorUnfolded; | |
112 | fSkipBinsUnfolded = other.fSkipBinsUnfolded; | |
113 | fExtraBinsUnfolded = other.fExtraBinsUnfolded; | |
114 | fbVariableBinning = other.fbVariableBinning; | |
115 | fPtMaxUnfVarBinning = other.fPtMaxUnfVarBinning; | |
116 | f1MergeFunction = other.f1MergeFunction; | |
117 | fFineFrac = other.fFineFrac; | |
118 | fbCalcErrors = other.fbCalcErrors; | |
119 | ||
120 | return *this; | |
31d3e4f0 | 121 | } |
122 | ||
123 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
124 | void AliAnaChargedJetResponseMaker::SetMeasuredSpectrum(TH1D *hPtMeasured) | |
125 | { | |
126 | // | |
127 | // Set measured spectrum in THnSparse format | |
128 | // This defines the minimum and maximum pT on the reconstructed axis of the response matrix | |
129 | // | |
130 | if(fDebug) printf(">>AliAnaChargedJetResponseMaker::SetMeasuredSpectrum \n"); | |
131 | ||
132 | fNbins = hPtMeasured->GetXaxis()->GetNbins(); | |
133 | fPtMin = hPtMeasured->GetXaxis()->GetXmin(); | |
134 | fPtMax = hPtMeasured->GetXaxis()->GetXmax(); | |
135 | ||
136 | if(fDebug) printf("fNbins: %d fPtMin: %f fPtMax: %f \n",fNbins,fPtMin,fPtMax); | |
137 | ||
138 | if(fBinArrayPtRec) delete fBinArrayPtRec; | |
139 | fBinArrayPtRec = new Double_t[fNbins+1]; | |
140 | for(int j = 0; j<fNbins; j++) { | |
141 | fBinArrayPtRec[j] = hPtMeasured->GetXaxis()->GetBinLowEdge(j+1); | |
142 | } | |
143 | fBinArrayPtRec[fNbins] = hPtMeasured->GetXaxis()->GetBinUpEdge(fNbins); | |
144 | ||
145 | ||
146 | Int_t nbins[fDimensions]; | |
147 | Double_t xmin[fDimensions]; | |
148 | Double_t xmax[fDimensions]; | |
149 | for(int dim = 0; dim<fDimensions; dim++) { | |
150 | nbins[dim] = fNbins; | |
151 | xmin[dim] = fPtMin; | |
152 | xmax[dim] = fPtMax; | |
153 | } | |
154 | ||
155 | if(fPtMeasured) delete fPtMeasured; | |
156 | fPtMeasured = new THnSparseD("fPtMeasured","Measured pT spectrum; p_{T}^{rec}",fDimensions,nbins,xmin,xmax); | |
157 | fPtMeasured->Sumw2(); | |
158 | fPtMeasured->GetAxis(0)->SetTitle("p_{T}^{rec}"); | |
159 | fPtMeasured->SetBinEdges(0,fBinArrayPtRec); | |
160 | ||
161 | //Fill | |
162 | if(fDebug) printf("fill measured THnSparse\n"); | |
163 | if(fNbins!=hPtMeasured->GetNbinsX()) | |
164 | printf("WARNING: nbins not correct \t %d vs %d !!!\n",fNbins,hPtMeasured->GetNbinsX()); | |
165 | ||
166 | int bin[1] = {0}; | |
167 | bin[0] = 0; | |
168 | for(int i = hPtMeasured->FindBin(fPtMin); i<hPtMeasured->FindBin(fPtMax); i++) { | |
169 | double pT[1]; | |
170 | pT[0]= hPtMeasured->GetBinCenter(i); | |
171 | fPtMeasured->SetBinContent(bin,(Double_t)hPtMeasured->GetBinContent(i)); | |
172 | fPtMeasured->SetBinError(bin,(Double_t)hPtMeasured->GetBinError(i)); | |
173 | bin[0]++; | |
174 | } | |
175 | ||
176 | if(fDebug) printf("fPtMeasured->GetNbins(): %d \n",(int)fPtMeasured->GetNbins()); | |
177 | ||
178 | } | |
179 | ||
180 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
181 | void AliAnaChargedJetResponseMaker::SetFlatEfficiency(Double_t eff) { | |
182 | ||
ef62323a | 183 | // |
184 | // Set flat efficiency to value of eff | |
185 | // | |
186 | ||
31d3e4f0 | 187 | fEffFlat = eff; |
ef62323a | 188 | return; |
31d3e4f0 | 189 | |
275d3481 | 190 | /* |
31d3e4f0 | 191 | Int_t nbins[fDimensions]; |
192 | Double_t xmin[fDimensions]; | |
193 | Double_t xmax[fDimensions]; | |
194 | for(int dim = 0; dim<fDimensions; dim++) { | |
195 | nbins[dim] = fNbins; | |
196 | xmin[dim] = fPtMin; | |
197 | xmax[dim] = fPtMax; | |
198 | } | |
199 | ||
200 | if(fEfficiency) delete fEfficiency; | |
201 | fEfficiency = new THnSparseD("fEfficiency","Efficiency - no resolution effects; p_{T}^{gen}",fDimensions,nbins,xmin,xmax); | |
202 | fEfficiency->Sumw2(); | |
203 | fEfficiency->GetAxis(0)->SetTitle("p_{T}^{gen}"); | |
204 | fEfficiency->SetBinEdges(0,fBinArrayPtRec); | |
205 | ||
206 | for(int i=0; i<fNbins; i++) { | |
207 | int bin[1]; | |
208 | bin[0] = i; | |
209 | fEfficiency->SetBinContent(bin,fEffFlat); | |
210 | fEfficiency->SetBinError(bin,0); | |
211 | } | |
5b045191 | 212 | */ |
31d3e4f0 | 213 | } |
214 | ||
215 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
216 | void AliAnaChargedJetResponseMaker::SetEfficiency(TGraphErrors *grEff) | |
217 | { | |
ef62323a | 218 | // |
219 | // Fill fEfficiency (THnSparse) with values from grEff | |
220 | // | |
221 | ||
31d3e4f0 | 222 | Int_t nbins[fDimensions]; |
223 | Double_t xmin[fDimensions]; | |
224 | Double_t xmax[fDimensions]; | |
225 | for(int dim = 0; dim<fDimensions; dim++) { | |
226 | nbins[dim] = fNbins; | |
227 | xmin[dim] = fPtMin; | |
228 | xmax[dim] = fPtMax; | |
229 | } | |
230 | ||
231 | if(fEfficiency) delete fEfficiency; | |
232 | fEfficiency = new THnSparseD("fEfficiency","Efficiency - no resolution effects; p_{T}^{gen}",fDimensions,nbins,xmin,xmax); | |
233 | fEfficiency->Sumw2(); | |
234 | fEfficiency->GetAxis(0)->SetTitle("p_{T}^{gen}"); | |
235 | fEfficiency->SetBinEdges(0,fBinArrayPtRec); | |
236 | ||
237 | double pT[1]; | |
238 | double yield = 0.; | |
239 | double error = 0.; | |
240 | double dummy = 0.; | |
241 | int nbinsgr = grEff->GetN(); | |
242 | ||
243 | for(int i=0; i<nbinsgr; i++) { | |
244 | grEff->GetPoint(i,dummy,yield); | |
245 | pT[0] = dummy; | |
246 | error = grEff->GetErrorY(i); | |
247 | ||
248 | fEfficiency->Fill(pT,yield); | |
249 | fEfficiency->SetBinError(i,error); | |
250 | ||
251 | } | |
252 | ||
253 | } | |
254 | ||
255 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
256 | void AliAnaChargedJetResponseMaker::MakeResponseMatrixJetsFineMerged(Int_t skipBins, Int_t binWidthFactor, Int_t extraBins, Bool_t bVariableBinning, Double_t ptmax) | |
257 | { | |
258 | // | |
259 | // Make jet response matrix | |
260 | // | |
261 | ||
ef62323a | 262 | if(fDebug) printf(">>AliAnaChargedJetResponseMaker::MakeResponseMatrixJetsFineMerged\n"); |
263 | ||
264 | if(!fPtMeasured) { | |
265 | if(fDebug) printf("fPtMeasured does not exist. Aborting!!\n"); | |
266 | return; | |
267 | } | |
31d3e4f0 | 268 | if(fResponseMatrix) { delete fResponseMatrix; } |
269 | if(fResponseMatrixFine) { delete fResponseMatrixFine; } | |
270 | ||
271 | SetSkipBinsUnfolded(skipBins); | |
272 | SetBinWidthFactorUnfolded(binWidthFactor); | |
273 | SetExtraBinsUnfolded(extraBins); | |
274 | SetVariableBinning(bVariableBinning,ptmax); | |
275 | ||
276 | InitializeResponseMatrix(); | |
277 | InitializeEfficiency(); | |
278 | ||
279 | InitializeResponseMatrixFine(); | |
280 | InitializeEfficiencyFine(); | |
281 | ||
282 | FillResponseMatrixFineAndMerge(); | |
283 | ||
284 | } | |
285 | ||
286 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
287 | void AliAnaChargedJetResponseMaker::InitializeResponseMatrix() { | |
288 | // | |
ef62323a | 289 | //Define bin edges of RM to be used for unfolding |
31d3e4f0 | 290 | // |
291 | ||
292 | Int_t nbins[fDimensions*2]; | |
293 | nbins[fDimRec] = fNbins; | |
fa7c34ba | 294 | nbins[fDimGen] = fNbins; |
31d3e4f0 | 295 | |
296 | double binWidthMeas = (double)((fPtMax-fPtMin)/fNbins); | |
297 | double binWidthUnf = (double)fBinWidthFactorUnfolded*binWidthMeas; | |
298 | double binWidthUnfLowPt = -1.; | |
299 | if(fbVariableBinning) | |
300 | binWidthUnfLowPt = binWidthUnf*0.5; | |
301 | ||
03372fd1 | 302 | fPtMaxUnfolded = fPtMax+(double)(fExtraBinsUnfolded)*binWidthUnf; |
303 | nbins[fDimGen]+=fExtraBinsUnfolded; | |
304 | ||
31d3e4f0 | 305 | |
306 | printf("fPtMinMeas: %f fPtMaxMeas: %f\n",fPtMin,fPtMax); | |
307 | printf("binWidthMeas: %f binWidthUnf: %f fBinWidthFactorUnfolded: %d\n",binWidthMeas,binWidthUnf,fBinWidthFactorUnfolded); | |
308 | printf("binWidthUnfLowPt: %f\n",binWidthUnfLowPt); | |
309 | ||
adf3920d | 310 | printf("fPtMinUnfolded: %f fPtMaxUnfolded: %f\n",fPtMinUnfolded,fPtMaxUnfolded); |
311 | ||
31d3e4f0 | 312 | |
313 | if(fbVariableBinning) { | |
adf3920d | 314 | // cout << "fPtMaxUnfVarBinning: " << fPtMaxUnfVarBinning << " \tfPtMinUnfolded: " << fPtMinUnfolded << " binWidthUnfLowPt: " << binWidthUnfLowPt << endl; |
315 | Int_t tmp = (int)((fPtMaxUnfVarBinning-fPtMinUnfolded)/binWidthUnfLowPt); | |
31d3e4f0 | 316 | tmp = tmp - fSkipBinsUnfolded; |
317 | fPtMinUnfolded = fPtMaxUnfVarBinning-(double)(tmp)*binWidthUnfLowPt; | |
adf3920d | 318 | //cout << "tmp = " << tmp << " fSkipBinsUnfolded = " << fSkipBinsUnfolded << " fPtMinUnfolded = " << fPtMinUnfolded << endl; |
31d3e4f0 | 319 | //Redefine also number of bins on generated axis in case of variable binning |
320 | nbins[fDimGen] = (int)((fPtMaxUnfVarBinning-fPtMinUnfolded)/binWidthUnfLowPt)+(int)((fPtMaxUnfolded-fPtMaxUnfVarBinning)/binWidthUnf); | |
321 | } | |
adf3920d | 322 | else { |
323 | int tmp = round((fPtMaxUnfolded-fPtMinUnfolded)/binWidthUnf); //nbins which fit between 0 and fPtMaxUnfolded | |
324 | tmp = tmp - fSkipBinsUnfolded; | |
325 | fPtMinUnfolded = fPtMaxUnfolded-(double)(tmp)*binWidthUnf; //set ptmin unfolded | |
326 | fPtMaxUnfolded = fPtMinUnfolded+(double)(tmp)*binWidthUnf; //set ptmax unfolded | |
327 | nbins[fDimGen] = (int)((fPtMaxUnfolded-fPtMinUnfolded)/binWidthUnf); //adjust nbins to bin width | |
328 | } | |
329 | ||
330 | printf(" nbins[fDimGen] = %d nbins[fDimRec] = %d\n",nbins[fDimGen],nbins[fDimRec]); | |
31d3e4f0 | 331 | |
332 | Double_t binWidth[2]; | |
333 | binWidth[fDimRec] = (double)((fPtMax-fPtMin)/nbins[fDimRec]); | |
334 | binWidth[fDimGen] = (double)((fPtMaxUnfolded-fPtMinUnfolded)/nbins[fDimGen]); | |
335 | ||
336 | Double_t xmin[fDimensions*2]; | |
337 | Double_t xmax[fDimensions*2]; | |
338 | xmin[fDimRec] = fPtMin; | |
339 | xmax[fDimRec] = fPtMax; | |
340 | xmin[fDimGen] = fPtMinUnfolded; | |
341 | xmax[fDimGen] = fPtMaxUnfolded; | |
342 | ||
343 | printf("xmin[fDimRec]: %f xmin[fDimGen]: %f \n",xmin[fDimRec],xmin[fDimGen]); | |
344 | printf("xmax[fDimRec]: %f xmax[fDimGen]: %f \n",xmax[fDimRec],xmax[fDimGen]); | |
345 | printf("nbins[fDimRec]: %d nbins[fDimGen]: %d \n",nbins[fDimRec],nbins[fDimGen]); | |
346 | if(!fbVariableBinning) printf("binWidth[fDimRec]: %f binWidth[fDimGen]: %f \n",binWidth[fDimRec],binWidth[fDimGen]); | |
347 | ||
348 | Double_t binArrayPt0[nbins[fDimRec]+1]; | |
349 | Double_t binArrayPt1[nbins[fDimGen]+1]; | |
350 | Double_t xmaxGen = TMath::Max(xmax[fDimGen],fPtMaxUnfoldedHigh); | |
351 | ||
352 | //Define bin limits reconstructed/measured axis | |
353 | for(int i=0; i<nbins[fDimRec]; i++) { | |
354 | binArrayPt0[i] = xmin[fDimRec]+(double)i*binWidth[fDimRec]; | |
355 | } | |
356 | binArrayPt0[nbins[fDimRec]]= xmax[fDimRec]; | |
357 | ||
358 | //Define bin limits generated/unfolded axis | |
fa7c34ba | 359 | binArrayPt1[0] = fPtMinUnfolded; |
360 | for(int i=1; i<nbins[fDimGen]; i++) { | |
31d3e4f0 | 361 | if(fbVariableBinning) { |
362 | double test = xmin[fDimGen]+(double)i*binWidthUnfLowPt; | |
363 | if(test<=fPtMaxUnfVarBinning) binArrayPt1[i] = test; | |
364 | else binArrayPt1[i] = binArrayPt1[i-1]+binWidthUnf; | |
365 | } | |
fa7c34ba | 366 | else |
367 | binArrayPt1[i] = xmin[fDimGen]+(double)i*binWidth[fDimGen]; | |
31d3e4f0 | 368 | //printf("RM. i = %d \t binArrayPt[i] = %f \n",i,binArrayPt1[i]); |
369 | } | |
370 | binArrayPt1[nbins[fDimGen]]= xmaxGen; | |
371 | ||
372 | ||
373 | // Response matrix : dimensions must be 2N = 2 x (number of variables) | |
374 | // dimensions 0 -> N-1 must be filled with reconstructed values | |
375 | // dimensions N -> 2N-1 must be filled with generated values | |
ef62323a | 376 | if(fResponseMatrix) delete fResponseMatrix; |
31d3e4f0 | 377 | fResponseMatrix = new THnSparseD("fResponseMatrix","Response Matrix pTMC vs pTrec",fDimensions*2,nbins,xmin,xmax); |
378 | fResponseMatrix->Sumw2(); | |
379 | fResponseMatrix->GetAxis(fDimRec)->SetTitle("p_{T}^{rec}"); | |
380 | fResponseMatrix->GetAxis(fDimGen)->SetTitle("p_{T}^{gen}"); | |
381 | ||
382 | fResponseMatrix->SetBinEdges(fDimRec,binArrayPt0); | |
383 | fResponseMatrix->SetBinEdges(fDimGen,binArrayPt1); | |
384 | ||
ef62323a | 385 | Int_t bin[2] = {1,1}; |
386 | for(int i=1; i<fResponseMatrix->GetAxis(0)->GetNbins(); i++) { | |
387 | bin[0]=i; | |
388 | for(int j=1; j<fResponseMatrix->GetAxis(1)->GetNbins(); j++) { | |
389 | bin[1]=j; | |
390 | fResponseMatrix->SetBinContent(bin,0.); | |
391 | } | |
392 | } | |
393 | ||
394 | ||
31d3e4f0 | 395 | |
396 | } | |
397 | ||
398 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
399 | void AliAnaChargedJetResponseMaker::InitializeEfficiency() { | |
400 | // | |
401 | // Define dimensions of efficiency THnSparse | |
402 | // | |
403 | ||
404 | if(!fResponseMatrix) { | |
405 | printf("AliAnaChargedJetResponseMaker::InitializeEfficiency()\n"); | |
406 | printf("Can not define dimensions efficiency without fResponseMatrix dimensions. Aborting \n"); | |
407 | return; | |
408 | } | |
409 | ||
410 | TAxis *genAxis = fResponseMatrix->GetAxis(fDimGen); | |
411 | ||
412 | Int_t nbinsEff[fDimensions]; | |
413 | Double_t xminEff[fDimensions]; | |
414 | Double_t xmaxEff[fDimensions]; | |
415 | ||
416 | for(int dim = 0; dim<fDimensions; dim++) { | |
417 | nbinsEff[dim] = genAxis->GetNbins(); | |
418 | xminEff[dim] = genAxis->GetXmin(); | |
419 | xmaxEff[dim] = genAxis->GetXmax(); | |
420 | } | |
421 | ||
422 | if(fEfficiency) delete fEfficiency; | |
423 | fEfficiency = new THnSparseD("fEfficiency","Efficiency - no resolution effects; p_{T}^{gen}",fDimensions,nbinsEff,xminEff,xmaxEff); | |
424 | fEfficiency->Sumw2(); | |
425 | fEfficiency->GetAxis(0)->SetTitle("p_{T}^{gen}"); | |
426 | ||
427 | const Double_t *binArrayPt = genAxis->GetXbins()->GetArray(); | |
428 | fEfficiency->SetBinEdges(0,binArrayPt); | |
429 | ||
430 | } | |
431 | ||
432 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
433 | void AliAnaChargedJetResponseMaker::InitializeResponseMatrixFine() { | |
434 | // | |
435 | //Initialize fine response matrix | |
436 | // | |
437 | ||
438 | Int_t nbinsFine[fDimensions*2]; | |
fa7c34ba | 439 | Double_t xminFine[fDimensions*2]; |
440 | Double_t xmaxFine[fDimensions*2]; | |
31d3e4f0 | 441 | Double_t pTarrayFine[fDimensions*2]; |
442 | ||
fa7c34ba | 443 | nbinsFine[fDimRec] = fResponseMatrix->GetAxis(fDimRec)->GetNbins()*fFineFrac; |
444 | nbinsFine[fDimGen] = fResponseMatrix->GetAxis(fDimRec)->GetNbins()*fFineFrac; | |
445 | xminFine[fDimRec] = TMath::Min(fPtMin,0.); | |
446 | xminFine[fDimGen] = TMath::Min(fPtMin,0.); | |
447 | xmaxFine[fDimRec] = fResponseMatrix->GetAxis(fDimGen)->GetXmax()+40.; | |
448 | xmaxFine[fDimGen] = fResponseMatrix->GetAxis(fDimGen)->GetXmax()+40.; | |
449 | pTarrayFine[fDimRec] = 0.; | |
450 | pTarrayFine[fDimGen] = 0.; | |
31d3e4f0 | 451 | |
452 | Double_t binWidth[2]; | |
453 | binWidth[fDimRec] = fResponseMatrix->GetAxis(fDimRec)->GetBinWidth(1); | |
454 | ||
455 | Double_t binWidthFine[2]; | |
456 | binWidthFine[fDimRec] = binWidth[fDimRec]/((double)fFineFrac); | |
457 | binWidthFine[fDimGen] = binWidthFine[fDimRec]*(double)fBinWidthFactorUnfolded; | |
458 | ||
459 | nbinsFine[fDimRec] = (int)((xmaxFine[fDimRec]-xminFine[fDimRec])/binWidthFine[fDimRec]); //adjust nbins to bin width | |
460 | nbinsFine[fDimGen] = (int)((xmaxFine[fDimGen]-xminFine[fDimGen])/binWidthFine[fDimGen]); //adjust nbins to bin width | |
461 | ||
462 | printf("xminFine[fDimRec]: %f xminFine[fDimGen]: %f \n",xminFine[fDimRec],xminFine[fDimGen]); | |
463 | printf("xmaxFine[fDimRec]: %f xmaxFine[fDimGen]: %f \n",xmaxFine[fDimRec],xmaxFine[fDimGen]); | |
464 | printf("nbinsFine[fDimRec]: %d nbinsFine[fDimGen]: %d \n",nbinsFine[fDimRec],nbinsFine[fDimGen]); | |
465 | printf("binWidthFine[fDimRec]: %f binWidthFine[fDimGen]: %f \n",binWidthFine[fDimRec],binWidthFine[fDimGen]); | |
466 | ||
467 | ||
468 | Double_t binArrayPt0Fine[nbinsFine[fDimRec]+1]; | |
469 | Double_t binArrayPt1Fine[nbinsFine[fDimGen]+1]; | |
470 | ||
471 | for(int i=0; i<nbinsFine[fDimRec]; i++) { | |
472 | binArrayPt0Fine[i] = xminFine[fDimRec]+(double)i*binWidthFine[fDimRec]; | |
473 | // printf("RM. i = %d \t binArrayPtFine[i] = %f \n",i,binArrayPt0Fine[i]); | |
474 | } | |
475 | binArrayPt0Fine[nbinsFine[fDimRec]]= xmaxFine[fDimRec]; | |
476 | ||
477 | for(int i=0; i<nbinsFine[fDimGen]; i++) { | |
478 | binArrayPt1Fine[i] = xminFine[fDimGen]+(double)i*binWidthFine[fDimGen]; | |
479 | // printf("RM. i = %d \t binArrayPtFine[i] = %f \n",i,binArrayPt1Fine[i]); | |
480 | } | |
481 | binArrayPt1Fine[nbinsFine[fDimGen]]= xmaxFine[fDimGen]; | |
482 | ||
483 | // Response matrix : dimensions must be 2N = 2 x (number of variables) | |
484 | // dimensions 0 -> N-1 must be filled with reconstructed values | |
485 | // dimensions N -> 2N-1 must be filled with generated values | |
ef62323a | 486 | if(fResponseMatrixFine) delete fResponseMatrixFine; |
31d3e4f0 | 487 | fResponseMatrixFine = new THnSparseD("fResponseMatrixFine","Response Matrix pTMC vs pTrec",fDimensions*2,nbinsFine,xminFine,xmaxFine); |
488 | fResponseMatrixFine->Sumw2(); | |
489 | fResponseMatrixFine->GetAxis(fDimRec)->SetTitle("p_{T}^{rec}"); | |
490 | fResponseMatrixFine->GetAxis(fDimGen)->SetTitle("p_{T}^{gen}"); | |
491 | ||
492 | fResponseMatrixFine->SetBinEdges(fDimRec,binArrayPt0Fine); | |
493 | fResponseMatrixFine->SetBinEdges(fDimGen,binArrayPt1Fine); | |
494 | ||
ef62323a | 495 | Int_t bin[2] = {1,1}; |
496 | for(int i=1; i<fResponseMatrixFine->GetAxis(0)->GetNbins(); i++) { | |
497 | bin[0]=i; | |
498 | for(int j=1; j<fResponseMatrixFine->GetAxis(1)->GetNbins(); j++) { | |
499 | bin[1]=j; | |
500 | fResponseMatrixFine->SetBinContent(bin,0.); | |
501 | } | |
502 | } | |
503 | ||
31d3e4f0 | 504 | } |
505 | ||
506 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
507 | void AliAnaChargedJetResponseMaker::InitializeEfficiencyFine() { | |
508 | // | |
509 | // Define dimensions of efficiency THnSparse | |
510 | // | |
511 | ||
512 | if(!fResponseMatrixFine) { | |
513 | printf("AliAnaChargedJetResponseMaker::InitializeEfficiencyFine()\n"); | |
514 | printf("Can not define dimensions efficiency without fResponseMatrixFine dimensions. Aborting \n"); | |
515 | return; | |
516 | } | |
517 | ||
518 | TAxis *genAxis = fResponseMatrixFine->GetAxis(fDimGen); | |
519 | ||
520 | Int_t nbinsEff[fDimensions]; | |
521 | Double_t xminEff[fDimensions]; | |
522 | Double_t xmaxEff[fDimensions]; | |
523 | ||
524 | for(int dim = 0; dim<fDimensions; dim++) { | |
525 | nbinsEff[dim] = genAxis->GetNbins(); | |
526 | xminEff[dim] = genAxis->GetXmin(); | |
527 | xmaxEff[dim] = genAxis->GetXmax(); | |
528 | } | |
529 | ||
530 | if(fEfficiencyFine) delete fEfficiencyFine; | |
531 | fEfficiencyFine = new THnSparseD("fEfficiencyFine","EfficiencyFine - no resolution effects; p_{T}^{gen}",fDimensions,nbinsEff,xminEff,xmaxEff); | |
532 | fEfficiencyFine->Sumw2(); | |
533 | fEfficiencyFine->GetAxis(0)->SetTitle("p_{T}^{gen}"); | |
534 | ||
535 | const Double_t *binArrayPt = genAxis->GetXbins()->GetArray(); | |
536 | fEfficiencyFine->SetBinEdges(0,binArrayPt); | |
537 | ||
538 | } | |
539 | ||
540 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
541 | void AliAnaChargedJetResponseMaker::FillResponseMatrixFineAndMerge() { | |
542 | // | |
543 | // Fill fine response matrix | |
544 | // | |
545 | ||
ef62323a | 546 | if(!fResponseMatrix) { |
547 | printf("Dimensions of fResponseMatrix have to be defined first. Aborting!"); | |
548 | return; | |
549 | } | |
550 | ||
551 | if(!fResponseMatrixFine) { | |
552 | printf("Dimensions of fResponseMatrixFine have to be defined first. Aborting!"); | |
553 | return; | |
554 | } | |
555 | ||
31d3e4f0 | 556 | TAxis *genAxis = fResponseMatrixFine->GetAxis(fDimGen); |
557 | TAxis *recAxis = fResponseMatrixFine->GetAxis(fDimRec); | |
558 | ||
559 | Int_t nbinsFine[fDimensions*2]; | |
560 | Double_t xminFine[fDimensions*2]; | |
561 | Double_t xmaxFine[fDimensions*2]; | |
562 | Double_t pTarrayFine[fDimensions*2]; | |
563 | ||
564 | nbinsFine[fDimGen] = genAxis->GetNbins(); | |
565 | nbinsFine[fDimRec] = recAxis->GetNbins(); | |
566 | xminFine[fDimGen] = genAxis->GetXmin(); | |
567 | xminFine[fDimRec] = recAxis->GetXmin(); | |
568 | xmaxFine[fDimGen] = genAxis->GetXmax(); | |
569 | xmaxFine[fDimRec] = recAxis->GetXmax(); | |
fa7c34ba | 570 | pTarrayFine[fDimGen] = 0.; |
571 | pTarrayFine[fDimRec] = 0.; | |
31d3e4f0 | 572 | |
573 | double sumyield = 0.; | |
574 | double sumyielderror2 = 0.; | |
575 | ||
ada4df5e | 576 | int ipt[2] = {0,0}; |
577 | int iptMerged[2] = {0,0}; | |
578 | int ierror[2] = {0,0}; | |
31d3e4f0 | 579 | |
580 | Int_t check = 0; | |
581 | Double_t pTgen, pTrec; | |
582 | Double_t yield = 0.; | |
583 | Double_t error = 0.; | |
584 | ||
585 | const int nng = fResponseMatrix->GetAxis(fDimGen)->GetNbins(); | |
586 | const int nnr = fResponseMatrix->GetAxis(fDimRec)->GetNbins(); | |
587 | Double_t errorArray[nng][nnr]; | |
588 | for(int iig =0; iig<nng; iig++) { | |
589 | for(int iir =0; iir<nnr; iir++) { | |
590 | errorArray[iig][iir] = 0.; | |
591 | } | |
592 | } | |
593 | ||
594 | for(int iy=1; iy<=nbinsFine[fDimGen]; iy++) { //gen | |
595 | pTgen = fResponseMatrixFine->GetAxis(fDimGen)->GetBinCenter(iy); | |
596 | pTarrayFine[fDimGen] = pTgen; | |
597 | ierror[fDimGen]=iy; | |
598 | sumyield = 0.; | |
599 | check = 0; | |
600 | ||
601 | for(int ix=1; ix<=nbinsFine[fDimRec]; ix++) { //rec | |
602 | pTrec = fResponseMatrixFine->GetAxis(fDimRec)->GetBinCenter(ix); | |
603 | Double_t width = fResponseMatrixFine->GetAxis(fDimRec)->GetBinWidth(ix); | |
604 | if(fResolutionType==kParam) { | |
605 | yield = fDeltaPt->Eval(pTrec-pTgen); | |
606 | error = 0.; | |
607 | } | |
608 | else if(fResolutionType==kResiduals) { | |
609 | yield = fhDeltaPt->Interpolate(pTrec-pTgen); | |
610 | error = 0.; | |
611 | } | |
612 | else if(fResolutionType==kResidualsErr) { | |
613 | Double_t deltaPt = pTrec-pTgen; | |
614 | int bin = fhDeltaPt->FindBin(deltaPt); | |
615 | yield = fhDeltaPt->GetBinContent(bin); | |
616 | error = fhDeltaPt->GetBinError(bin); | |
617 | } | |
618 | yield=yield*width; | |
619 | error=error*width; | |
620 | //avoid trouble with empty bins in the high pT tail of deltaPt distribution | |
621 | if(check==0 && yield>0. && pTrec>pTgen) check=1; | |
622 | if(check==1 && yield==0.) ix=nbinsFine[fDimRec]; | |
623 | ||
624 | sumyield+=yield; | |
625 | sumyielderror2 += error*error; | |
626 | ||
627 | pTarrayFine[fDimRec] = pTrec; | |
628 | ierror[fDimRec] = ix; | |
629 | fResponseMatrixFine->Fill(pTarrayFine,yield); | |
630 | if(fbCalcErrors) fResponseMatrixFine->SetBinError(ierror,error); | |
631 | ||
632 | }// ix (dimRec) loop | |
633 | ||
634 | //Normalize to total number of counts =1 | |
635 | if(sumyield>1) { | |
636 | ipt[fDimGen]=iy; | |
637 | for(int ix=1; ix<=nbinsFine[fDimRec]; ix++) { | |
638 | ipt[fDimRec]=ix; | |
639 | fResponseMatrixFine->SetBinContent(ipt,fResponseMatrixFine->GetBinContent(ipt)/sumyield); | |
640 | if(fResolutionType==kResidualsErr && fbCalcErrors) { | |
641 | Double_t A = 1./sumyield*fResponseMatrix->GetBinError(ipt); | |
642 | Double_t B = -1.*fResponseMatrix->GetBinContent(ipt)/sumyield/sumyield*TMath::Sqrt(sumyielderror2); | |
643 | Double_t tmp2 = A*A + B*B; | |
644 | fResponseMatrix->SetBinError(ipt,TMath::Sqrt(tmp2)); | |
645 | } | |
646 | ||
647 | } | |
648 | } | |
649 | ||
650 | int bin[1]; | |
651 | bin[0] = iy; | |
652 | fEfficiencyFine->SetBinContent(bin,sumyield); | |
653 | if(fbCalcErrors) fEfficiencyFine->SetBinError(bin,TMath::Sqrt(sumyielderror2)); | |
654 | ||
655 | //fill merged response matrix | |
656 | ||
657 | //find bin in fine RM correspoinding to mimimum/maximum bin of merged RM on rec axis | |
658 | int ixMin = fResponseMatrixFine->GetAxis(fDimRec)->FindBin(fResponseMatrix->GetAxis(fDimRec)->GetXmin()); | |
659 | int ixMax = fResponseMatrixFine->GetAxis(fDimRec)->FindBin(fResponseMatrix->GetAxis(fDimRec)->GetXmax()); | |
660 | ||
661 | if(fResponseMatrixFine->GetAxis(fDimGen)->GetBinLowEdge(iy) >= fResponseMatrix->GetAxis(fDimGen)->GetXmin()) { | |
662 | ipt[fDimGen]=iy; | |
663 | iptMerged[fDimGen]=fResponseMatrix->GetAxis(fDimGen)->FindBin(pTgen); | |
664 | ||
665 | Double_t weight = 1.; | |
666 | if(f1MergeFunction) { | |
667 | Double_t loEdge = fResponseMatrix->GetAxis(fDimGen)->GetBinLowEdge(iptMerged[fDimGen]); | |
668 | Double_t upEdge = fResponseMatrix->GetAxis(fDimGen)->GetBinUpEdge(iptMerged[fDimGen]); | |
669 | Float_t powInteg = f1MergeFunction->Integral(loEdge,upEdge); | |
670 | //printf("loEdge = %f upEdge = %f powInteg = %f\n",loEdge,upEdge,powInteg); | |
671 | if(powInteg>0.) | |
672 | weight = f1MergeFunction->Integral(fResponseMatrixFine->GetAxis(fDimGen)->GetBinLowEdge(iy),fResponseMatrixFine->GetAxis(fDimGen)->GetBinUpEdge(iy))/powInteg; | |
673 | // printf("weight: %f \n", weight ); | |
674 | } else { | |
675 | weight = 1./((double)fFineFrac); | |
adf3920d | 676 | if(fbVariableBinning && pTgen<fPtMaxUnfVarBinning) weight=1./(0.5*(double)fFineFrac); |
31d3e4f0 | 677 | } |
678 | ||
679 | for(int ix=ixMin; ix<=ixMax; ix++) { //loop reconstructed axis | |
680 | pTrec = fResponseMatrixFine->GetAxis(fDimRec)->GetBinCenter(ix); | |
681 | ipt[fDimRec]=ix; | |
682 | iptMerged[fDimRec]=fResponseMatrix->GetAxis(fDimRec)->FindBin(pTrec); | |
683 | ||
684 | fResponseMatrix->AddBinContent(iptMerged,fResponseMatrixFine->GetBinContent(ipt)*weight); | |
685 | if(fbCalcErrors) errorArray[iptMerged[fDimGen]-1][iptMerged[fDimRec]-1] += fResponseMatrixFine->GetBinError(ipt)*fResponseMatrixFine->GetBinError(ipt)*weight*weight; | |
686 | } | |
687 | ||
688 | }//loop gen axis fine response matrix | |
689 | ||
690 | } // iy (dimGen) loop | |
691 | ||
692 | //Fill Efficiency corresponding to merged response matrix | |
31d3e4f0 | 693 | for(int iy=1; iy<=fResponseMatrix->GetAxis(fDimGen)->GetNbins(); iy++) { //gen |
694 | sumyield = 0.; | |
695 | ipt[fDimGen] = iy; | |
696 | ||
697 | for(int ix=1; ix<=fResponseMatrix->GetAxis(fDimRec)->GetNbins(); ix++) { //rec | |
698 | ipt[fDimRec] = ix; | |
699 | sumyield += fResponseMatrix->GetBinContent(ipt); | |
700 | ||
701 | if(fbCalcErrors) fResponseMatrix->SetBinError(ipt,TMath::Sqrt(errorArray[iy][ix])); | |
702 | } | |
703 | printf("RM: pTgen: %f \t sumyield: %f \n",fResponseMatrix->GetAxis(fDimGen)->GetBinCenter(iy),sumyield); | |
704 | int bin[1]; | |
705 | bin[0] = iy; | |
706 | fEfficiency->SetBinContent(bin,sumyield); | |
707 | if(fbCalcErrors) fEfficiency->SetBinError(bin,0); | |
708 | } | |
709 | ||
710 | if(fDebug) printf("fResponseMatrixFine->GetNbins(): %d \n",(int)fResponseMatrixFine->GetNbins()); | |
711 | if(fDebug) printf("fResponseMatrix->GetNbins(): %d \n",(int)fResponseMatrix->GetNbins()); | |
712 | ||
713 | if(fDebug) printf("Done constructing response matrix\n"); | |
714 | ||
715 | } | |
716 | ||
ef62323a | 717 | //-------------------------------------------------------------------------------------------------------------------------------------------------- |
03372fd1 | 718 | TH2* AliAnaChargedJetResponseMaker::MakeResponseMatrixRebin(TH2 *hRMFine, TH2 *hRM, Bool_t useFunctionWeight) { |
ef62323a | 719 | |
720 | // | |
721 | // Rebin matrix hRMFine to dimensions of hRM | |
03372fd1 | 722 | // function returns matrix in TH2D format (hRM2) with dimensions from hRM |
ef62323a | 723 | // |
724 | ||
725 | TH2 *hRM2 = (TH2*)hRM->Clone("hRM2"); | |
726 | hRM2->Reset(); | |
727 | ||
03372fd1 | 728 | if(useFunctionWeight) cout << "Use function to do weighting" << endl; |
729 | ||
275d3481 | 730 | //First normalize columns of input |
ef62323a | 731 | const Int_t nbinsNorm = hRM2->GetNbinsX(); |
732 | cout << "nbinsNorm: " << nbinsNorm << endl; | |
733 | ||
734 | TArrayF *normVector = new TArrayF(nbinsNorm); | |
735 | ||
736 | for(int ix=1; ix<=hRM2->GetNbinsX(); ix++) { | |
737 | Double_t xLow = hRM2->GetXaxis()->GetBinLowEdge(ix); | |
738 | Double_t xUp = hRM2->GetXaxis()->GetBinUpEdge(ix); | |
739 | //cout << "xLow: " << xLow << " xUp: " << xUp << "\t center: " << hRM2->GetXaxis()->GetBinCenter(ix) << endl; | |
740 | Int_t binxLowFine = hRMFine->GetXaxis()->FindBin(xLow); | |
741 | Int_t binxUpFine = hRMFine->GetXaxis()->FindBin(xUp)-1; | |
742 | //cout << "xLowFine: " << hRMFine->GetXaxis()->GetBinLowEdge(binxLowFine) << "\txUpFine: " << hRMFine->GetXaxis()->GetBinUpEdge(binxUpFine) << endl; | |
03372fd1 | 743 | if(useFunctionWeight) |
744 | normVector->SetAt(f1MergeFunction->Integral(xLow,xUp),ix-1); | |
745 | else | |
746 | normVector->SetAt(hRMFine->Integral(binxLowFine,binxUpFine,1,hRMFine->GetYaxis()->GetNbins()),ix-1); | |
747 | if(fDebug) cout << "ix norm: " << normVector->At(ix-1) << endl; | |
ef62323a | 748 | } |
749 | ||
750 | Double_t content, oldcontent = 0.; | |
751 | Int_t ixNew = 0; | |
752 | Int_t iyNew = 0; | |
275d3481 | 753 | Double_t xvalLo, xvalUp, yvalLo, yvalUp; |
ef62323a | 754 | Double_t xmin = hRM2->GetXaxis()->GetXmin(); |
755 | Double_t ymin = hRM2->GetYaxis()->GetXmin(); | |
756 | Double_t xmax = hRM2->GetXaxis()->GetXmax(); | |
757 | Double_t ymax = hRM2->GetYaxis()->GetXmax(); | |
758 | for(int ix=1; ix<=hRMFine->GetXaxis()->GetNbins(); ix++) { | |
275d3481 | 759 | xvalLo = hRMFine->GetXaxis()->GetBinLowEdge(ix); |
760 | xvalUp = hRMFine->GetXaxis()->GetBinUpEdge(ix); | |
761 | if(xvalLo<xmin || xvalUp>xmax) continue; | |
762 | ixNew = hRM2->GetXaxis()->FindBin(hRMFine->GetXaxis()->GetBinCenter(ix)); | |
ef62323a | 763 | for(int iy=1; iy<=hRMFine->GetYaxis()->GetNbins(); iy++) { |
275d3481 | 764 | yvalLo = hRMFine->GetYaxis()->GetBinLowEdge(iy); |
765 | yvalUp = hRMFine->GetYaxis()->GetBinUpEdge(iy); | |
766 | if(yvalLo<ymin || yvalUp>ymax) continue; | |
ef62323a | 767 | content = hRMFine->GetBinContent(ix,iy); |
275d3481 | 768 | iyNew = hRM2->GetYaxis()->FindBin(hRMFine->GetYaxis()->GetBinCenter(iy)); |
ef62323a | 769 | oldcontent = hRM2->GetBinContent(ixNew,iyNew); |
770 | ||
275d3481 | 771 | //if(fDebug) cout << "ixNew: " << ixNew << " " << xvalLo << " iyNew: " << iyNew << " " << yvalLo << " content: " << content << " oldContent: " << oldcontent << " newContent: " << oldcontent+content << endl; |
ef62323a | 772 | Double_t weight = 1.; |
03372fd1 | 773 | if(normVector->At(ixNew-1)>0.) { |
774 | if(useFunctionWeight) | |
775 | weight = f1MergeFunction->Integral(xvalLo,xvalUp)/normVector->At(ixNew-1); | |
776 | else | |
777 | weight = 1./normVector->At(ixNew-1); | |
778 | } | |
ef62323a | 779 | hRM2->SetBinContent(ixNew,iyNew,oldcontent+content*weight); |
780 | } | |
781 | } | |
782 | ||
783 | if(normVector) delete normVector; | |
784 | ||
785 | return hRM2; | |
786 | ||
787 | } | |
788 | ||
03372fd1 | 789 | //-------------------------------------------------------------------------------------------------------------------------------------------------- |
790 | TH2* AliAnaChargedJetResponseMaker::CreateTruncated2DHisto(TH2 *h2, Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax) { | |
791 | ||
792 | // | |
793 | // Limit axis range of 2D histogram | |
794 | // | |
795 | ||
796 | Int_t binMinXh2 = h2->GetXaxis()->FindBin(xmin); | |
797 | if(h2->GetXaxis()->GetBinLowEdge(binMinXh2) < xmin ) binMinXh2++; | |
798 | if(h2->GetXaxis()->GetBinLowEdge(binMinXh2) > xmin ) binMinXh2--; | |
799 | ||
800 | Int_t binMinYh2 = h2->GetYaxis()->FindBin(ymin); | |
801 | if(h2->GetYaxis()->GetBinLowEdge(binMinYh2) < ymin ) binMinYh2++; | |
802 | if(h2->GetYaxis()->GetBinLowEdge(binMinYh2) > ymin ) binMinYh2--; | |
803 | ||
804 | Int_t binMaxXh2 = h2->GetXaxis()->FindBin(xmax); | |
805 | if(h2->GetXaxis()->GetBinUpEdge(binMaxXh2) < xmax ) binMaxXh2++; | |
806 | if(h2->GetXaxis()->GetBinUpEdge(binMaxXh2) > xmax ) binMaxXh2--; | |
807 | ||
808 | Int_t binMaxYh2 = h2->GetYaxis()->FindBin(ymax); | |
809 | if(h2->GetYaxis()->GetBinUpEdge(binMaxYh2) < ymax ) binMaxYh2++; | |
810 | if(h2->GetYaxis()->GetBinUpEdge(binMaxYh2) > ymax ) binMaxYh2--; | |
811 | ||
812 | Int_t nbinsX = binMaxXh2-binMinXh2; | |
813 | Int_t nbinsY = binMaxYh2-binMinYh2; | |
814 | ||
815 | Double_t *binsX = new Double_t[nbinsX+1]; | |
816 | Double_t *binsY = new Double_t[nbinsY+1]; | |
817 | ||
818 | for(int ix=1; ix<=nbinsX; ix++) | |
819 | binsX[ix-1] = h2->GetXaxis()->GetBinLowEdge(binMinXh2+ix-1); | |
820 | binsX[nbinsX] = h2->GetXaxis()->GetBinUpEdge(binMaxXh2); | |
821 | ||
822 | for(int iy=1; iy<=nbinsY; iy++) | |
823 | binsY[iy-1] = h2->GetYaxis()->GetBinLowEdge(binMinYh2+iy-1); | |
824 | binsY[nbinsY] = h2->GetYaxis()->GetBinUpEdge(binMaxYh2); | |
825 | ||
826 | TH2 *h2Lim = new TH2D("h2Lim","h2Lim",nbinsX,binsX,nbinsY,binsY); | |
827 | ||
828 | for(int ix=1; ix<=nbinsX; ix++) { | |
829 | // cout << "ix: " << ix << " " << binsX[ix] << endl; | |
830 | for(int iy=1; iy<=nbinsY; iy++) { | |
831 | cout << "ix: " << ix << " " << binsX[ix] << "\tiy: " << iy << " " << binsY[iy] << endl; | |
832 | ||
833 | double content = h2->GetBinContent(binMinXh2+ix-1,binMinYh2+iy-1); | |
834 | double error = h2->GetBinContent(binMinXh2+ix-1,binMinYh2+iy-1); | |
835 | h2Lim->SetBinContent(ix,iy,content); | |
836 | h2Lim->SetBinError(ix,iy,error); | |
837 | ||
838 | } | |
839 | } | |
840 | ||
841 | ||
842 | ||
843 | return h2Lim; | |
844 | } | |
845 | ||
846 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
847 | TH2* AliAnaChargedJetResponseMaker::TruncateAxisRangeResponseMatrix(TH2 *hRMOrig, Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax) { | |
848 | ||
849 | // | |
850 | // Limit axis range of response matrix without changing normalization | |
851 | // | |
852 | ||
853 | //TH2 *hRMLimit | |
854 | //TH2 *hRMLimit2 = (TH2*)hRMLimit->Clone("hRMLimit2"); | |
855 | ||
856 | TH2 *hRMLimit2 = CreateTruncated2DHisto(hRMOrig, xmin, xmax, ymin, ymax); | |
857 | hRMLimit2->Reset(); | |
858 | ||
859 | double binCent[2] = {0.,0.}; | |
860 | double content = 0.; | |
861 | double error = 0.; | |
862 | Int_t binOrig[2] = {0}; | |
863 | for(int ix=1; ix<=hRMLimit2->GetXaxis()->GetNbins(); ix++) { | |
864 | binCent[0] = hRMLimit2->GetXaxis()->GetBinCenter(ix); | |
865 | binOrig[0] = hRMOrig->GetXaxis()->FindBin(binCent[0]); | |
866 | for(int iy=1; iy<=hRMLimit2->GetYaxis()->GetNbins(); iy++) { | |
867 | binCent[1] = hRMLimit2->GetYaxis()->GetBinCenter(iy); | |
868 | binOrig[1] = hRMOrig->GetYaxis()->FindBin(binCent[1]); | |
869 | ||
870 | content = hRMOrig->GetBinContent(binOrig[0],binOrig[1]); | |
871 | error = hRMOrig->GetBinError(binOrig[0],binOrig[1]); | |
872 | ||
873 | hRMLimit2->SetBinContent(ix,iy,content); | |
874 | hRMLimit2->SetBinError(ix,iy,error); | |
875 | ||
876 | } | |
877 | } | |
878 | ||
879 | ||
880 | return hRMLimit2; | |
881 | ||
882 | } | |
883 | ||
275d3481 | 884 | //-------------------------------------------------------------------------------------------------------------------------------------------------- |
885 | TH2* AliAnaChargedJetResponseMaker::MultiplityResponseMatrices(TH2 *h2RMDeltaPt, TH2 *h2RMDetector) { | |
886 | ||
887 | // | |
888 | // Make combined response matrix (background flucutuations + detector effects) | |
889 | // | |
890 | // hRMDeltaPt is the response matrix for background fluctuations from \delta(p_t) measurement | |
891 | // hRMDetector is the response matrix for detector effects: needs to be a squared matrix with on each axis the same bins as on the generated axis of the bkg fluctuations response matrix | |
892 | // | |
893 | // Function assumes that generated/unfolded axis is x-axis and reconstructed is on y-axis on both respone matrices | |
894 | ||
895 | ||
03372fd1 | 896 | TH2D *h2ResponseMatrixCombined = (TH2D*)h2RMDeltaPt->Clone("h2ResponseMatrixCombined"); //h2RMDeltaPt is the bkg fluctuations RM which has the dimensions we want for the combined response matrix |
275d3481 | 897 | h2ResponseMatrixCombined->SetTitle("h2ResponseMatrixCombined"); |
898 | h2ResponseMatrixCombined->SetName("h2ResponseMatrixCombined"); | |
899 | ||
900 | // M = RM_deltaPt * RM_DetEffects * T (M=measured T=truth) | |
901 | // Dimensions: | |
902 | // mx1 = mxd * dxt * tx1 | |
903 | // m = measured bins | |
904 | // t = truth bins | |
905 | // d = rec for RM_DetEffects and gen for RM_deltaPt | |
906 | // RM_comb = RM_deltaPt * RM_DetEffects (dimensions mxt) | |
907 | // RM_comb(m,t) = Sum_d RM_deltaPt(m,d)*RM_DetEffects(d,t) | |
908 | ||
909 | if(fDebug) { | |
03372fd1 | 910 | printf("Nt=%d\n",h2ResponseMatrixCombined->GetNbinsX()); |
911 | printf("Nm=%d\n",h2ResponseMatrixCombined->GetNbinsY()); | |
912 | printf("Nd=%d\n",h2RMDeltaPt->GetNbinsX()); | |
275d3481 | 913 | } |
914 | ||
03372fd1 | 915 | |
275d3481 | 916 | for(Int_t t=1; t<=h2ResponseMatrixCombined->GetNbinsX();t++) { |
917 | for(Int_t m=1; m<=h2ResponseMatrixCombined->GetNbinsY();m++) { | |
918 | Double_t valueSum = 0.; | |
919 | for(Int_t d=1; d<=h2RMDeltaPt->GetNbinsX();d++) { | |
920 | valueSum += h2RMDeltaPt->GetBinContent(d,m) * h2RMDetector->GetBinContent(t,d); | |
03372fd1 | 921 | // if(t==10 && m==10) cout << "sum m,d=" << m << "," << d << endl; |
275d3481 | 922 | }//d-loop |
03372fd1 | 923 | // if(t==10) cout << "t,m = " << t << "," << m << "\tvalueSum: " << valueSum << endl; |
275d3481 | 924 | h2ResponseMatrixCombined->SetBinContent(t,m,valueSum); |
925 | } //m-loop | |
926 | }//t-loop | |
927 | ||
928 | return h2ResponseMatrixCombined; | |
929 | ||
930 | } | |
31d3e4f0 | 931 | |
5d87a047 | 932 | //-------------------------------------------------------------------------------------------------------------------------------------------------- |
933 | TH2* AliAnaChargedJetResponseMaker::GetTransposeResponsMatrix(TH2 *h2RM) { | |
934 | // | |
935 | // Transpose response matrix | |
936 | // | |
937 | ||
938 | //Initialize transposed response matrix h2RMTranspose | |
939 | TArrayD *arrayX = (TArrayD*)h2RM->GetXaxis()->GetXbins(); | |
940 | for(int i=0; i<arrayX->GetSize(); i++) cout << "i: " << arrayX->At(i) << endl; | |
941 | Double_t *xbinsArray = arrayX->GetArray(); | |
942 | ||
943 | TArrayD *arrayY = (TArrayD*)h2RM->GetYaxis()->GetXbins(); | |
944 | for(int i=0; i<arrayY->GetSize(); i++) cout << "i: " << arrayY->At(i) << endl; | |
945 | Double_t *ybinsArray = arrayY->GetArray(); | |
946 | ||
947 | TH2D *h2RMTranspose = new TH2D("h2RMTranspose","h2RMTranspose",h2RM->GetNbinsY(),ybinsArray,h2RM->GetNbinsX(),xbinsArray); | |
948 | ||
949 | //Fill tranposed response matrix | |
950 | for (Int_t ibin = 1; ibin <= h2RMTranspose->GetNbinsX(); ibin++) { | |
951 | for (Int_t jbin = 1; jbin <= h2RMTranspose->GetNbinsY(); jbin++) { | |
952 | h2RMTranspose->SetBinContent(ibin,jbin, h2RM->GetBinContent(jbin,ibin)); | |
953 | } | |
954 | } | |
955 | ||
956 | ||
957 | return h2RMTranspose; | |
958 | ||
959 | } | |
960 | ||
961 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
962 | TH2* AliAnaChargedJetResponseMaker::NormalizeResponsMatrixYaxisWithPrior(TH2 *h2RM, TH1 *hPrior) { | |
963 | // | |
964 | // Normalize such that the Y projection is the prior | |
965 | // | |
966 | ||
03372fd1 | 967 | // TH1D *hProjRespY = (TH1D*)h2RM->ProjectionY("hProjRespY"); |
5d87a047 | 968 | double intPrior = hPrior->Integral();//"width"); |
969 | for (Int_t jbin = 1; jbin <= h2RM->GetNbinsY(); jbin++) { | |
970 | // double corr = hPrior->GetBinContent(jbin)/hProjRespY->GetBinContent(jbin); | |
971 | for (Int_t ibin = 1; ibin <= h2RM->GetNbinsX(); ibin++) { | |
972 | double content = h2RM->GetBinContent(ibin,jbin); | |
973 | // h2RM->SetBinContent(ibin,jbin,content*corr); | |
974 | h2RM->SetBinContent(ibin,jbin,hPrior->GetBinContent(jbin)/hPrior->GetBinWidth(jbin)/intPrior*content); | |
975 | } | |
976 | } | |
977 | ||
978 | return h2RM; | |
979 | } | |
980 | ||
31d3e4f0 | 981 | //-------------------------------------------------------------------------------------------------------------------------------------------------- |
982 | TH1D* AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TH1 *hGen, TH2 *hResponse,TH1 *hEfficiency,Bool_t bDrawSlices) { | |
983 | ||
984 | // | |
985 | // Multiply hGen with response matrix to obtain refolded spectrum | |
ef62323a | 986 | // Efficiency must be given. In case efficiency is 1 use SetFlatEfficiency(1.) before calling function |
31d3e4f0 | 987 | // |
988 | ||
ef62323a | 989 | if(!hEfficiency) { |
5d87a047 | 990 | // printf("Efficiency must be given. In case efficiency is 1 use SetFlatEfficiency(1.) before calling function. Aborting!"); |
991 | // return 0; | |
992 | printf("Setting efficiency to 1 \n"); | |
993 | hEfficiency = (TH1D*)hGen->Clone("hEfficiency"); | |
994 | hEfficiency->Reset(); | |
995 | for(int i=1; i<=hEfficiency->GetNbinsX(); i++) hEfficiency->SetBinContent(i,1.); | |
ef62323a | 996 | } |
997 | ||
31d3e4f0 | 998 | //For response |
999 | //x-axis: generated | |
1000 | //y-axis: reconstructed | |
1001 | if(fDebug>0) cout << "nbins hGen: " << hGen->GetNbinsX() << "\t nbins hResponseGen: " << hResponse->GetXaxis()->GetNbins() << "\t nbins hResponseRec: " << hResponse->GetYaxis()->GetNbins() << endl; | |
1002 | ||
1003 | TH1D *hRec = hResponse->ProjectionY("hRec"); | |
1004 | hRec->Sumw2(); | |
1005 | hRec->Reset(); | |
1006 | hRec->SetTitle("hRec"); | |
1007 | hRec->SetName("hRec"); | |
1008 | ||
1009 | for(int irec=1; irec<=hRec->GetNbinsX(); irec++) | |
1010 | hRec->SetBinContent(irec,0); | |
1011 | ||
1012 | TH1D *hRecGenBin = 0x0; | |
1013 | TCanvas *cSlices = 0x0; | |
1014 | if(bDrawSlices) { | |
1015 | cSlices = new TCanvas("cSlices","cSlices:Slices",400,400); | |
1016 | gPad->SetLogy(); | |
1017 | } | |
1018 | ||
1019 | Double_t yieldMC = 0.; | |
1020 | Double_t yieldMCerror = 0.; | |
1021 | Double_t sumYield = 0.; | |
1022 | const Int_t nbinsRec = hRec->GetNbinsX(); | |
1023 | Double_t sumError2[nbinsRec+1]; | |
31d3e4f0 | 1024 | Double_t eff = 0.; |
1025 | ||
1026 | for(int igen=1; igen<=hGen->GetNbinsX(); igen++) { | |
1027 | //get pTMC | |
1028 | sumYield = 0.; | |
1029 | if(fEffFlat>0.) | |
1030 | eff = fEffFlat; | |
1031 | else | |
1032 | eff = hEfficiency->GetBinContent(igen); | |
1033 | yieldMC = hGen->GetBinContent(igen)*eff; | |
adf3920d | 1034 | yieldMCerror = hGen->GetBinError(igen)*eff; |
31d3e4f0 | 1035 | |
1036 | if(bDrawSlices) { | |
1037 | hRecGenBin = hResponse->ProjectionY(Form("hRecGenBin%d",igen)); | |
1038 | hRecGenBin->Sumw2(); | |
1039 | hRecGenBin->Reset(); | |
1040 | hRecGenBin->SetTitle(Form("hRecGenBin%d",igen)); | |
1041 | hRecGenBin->SetName(Form("hRecGenBin%d",igen)); | |
1042 | } | |
1043 | ||
1044 | for(int irec=1; irec<=hRec->GetNbinsX(); irec++) { | |
1045 | hRec->AddBinContent(irec,yieldMC*hResponse->GetBinContent(igen,irec)); | |
1046 | sumYield+=hResponse->GetBinContent(igen,irec); | |
346cfe4a | 1047 | // Double_t A = yieldMC*hResponse->GetBinError(igen,irec); |
31d3e4f0 | 1048 | Double_t B = hResponse->GetBinContent(igen,irec)*yieldMCerror; |
346cfe4a | 1049 | // Double_t tmp2 = A*A + B*B; |
adf3920d | 1050 | //sumError2[irec-1] += tmp2 ; |
1051 | sumError2[irec-1] += B*B; | |
31d3e4f0 | 1052 | |
1053 | if(bDrawSlices) | |
1054 | hRecGenBin->SetBinContent(irec,yieldMC*hResponse->GetBinContent(igen,irec)); | |
1055 | ||
1056 | } | |
1057 | if(bDrawSlices) { | |
1058 | cSlices->cd(); | |
1059 | hRecGenBin->SetLineColor(igen); | |
1060 | if(igen==1) hRecGenBin->DrawCopy(); | |
1061 | else hRecGenBin->DrawCopy("same"); | |
1062 | } | |
1063 | ||
1064 | if(hRecGenBin) delete hRecGenBin; | |
1065 | ||
1066 | cout << "igen: " << igen << "\tpTMC: " << hGen->GetXaxis()->GetBinCenter(igen) << "\teff:" << eff << "\tsumYield: " << sumYield << endl; | |
1067 | } | |
1068 | ||
ef62323a | 1069 | for(int i=0; i<=nbinsRec; i++) { |
1070 | if(sumError2[i]>0.) | |
1071 | hRec->SetBinError(i+1,TMath::Sqrt(sumError2[i])); | |
1072 | } | |
31d3e4f0 | 1073 | |
1074 | ||
1075 | return hRec; | |
1076 | } | |
1077 | ||
1078 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
1079 | TH1D* AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TF1 *fGen, TH2 *hResponse,TH1 *hEfficiency) { | |
1080 | ||
ef62323a | 1081 | // |
1082 | // Multiply fGen function with response matrix to obtain (re)folded spectrum | |
1083 | // Efficiency must be given. In case efficiency is 1 use SetFlatEfficiency(1.) before calling function | |
1084 | // | |
1085 | ||
31d3e4f0 | 1086 | //For response |
1087 | //x-axis: generated | |
1088 | //y-axis: reconstructed | |
1089 | ||
ada4df5e | 1090 | if(fDebug>0) printf(">>AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TF1 *fGen, TH2 *hResponse,TH1 *hEfficiency)\n"); |
ef62323a | 1091 | |
31d3e4f0 | 1092 | TH1D *hRec = hResponse->ProjectionY("hRec"); |
1093 | hRec->Sumw2(); | |
1094 | hRec->Reset(); | |
1095 | hRec->SetTitle("hRec"); | |
1096 | hRec->SetName("hRec"); | |
1097 | ||
1098 | // TH1D *hRec = new TH1D("hRec","hRec",hResponse->GetNbinsY(),hResponse->GetYaxis()->GetXmin(),hResponse->GetYaxis()->GetXmax()); | |
1099 | ||
1100 | for(int irec=1; irec<=hRec->GetNbinsX(); irec++) | |
1101 | hRec->SetBinContent(irec,0); | |
1102 | ||
1103 | Double_t yieldMC = 0.; | |
1104 | Double_t sumYield = 0.; | |
1105 | Double_t eff = 0.; | |
1106 | for(int igen=1; igen<=hResponse->GetNbinsX(); igen++) { | |
1107 | //get pTMC | |
1108 | sumYield = 0.; | |
1109 | double pTMC = hResponse->GetXaxis()->GetBinCenter(igen); | |
ada4df5e | 1110 | if(hEfficiency) { |
1111 | int binEff = hEfficiency->FindBin(pTMC); | |
31d3e4f0 | 1112 | eff = hEfficiency->GetBinContent(binEff); |
ada4df5e | 1113 | } |
1114 | else eff = 1.; | |
1115 | // yieldMC = fGen->Eval(pTMC)*eff; | |
1116 | yieldMC = fGen->Integral(hResponse->GetXaxis()->GetBinLowEdge(igen),hResponse->GetXaxis()->GetBinUpEdge(igen))*eff; | |
31d3e4f0 | 1117 | for(int irec=1; irec<=hResponse->GetNbinsY(); irec++) { |
1118 | hRec->AddBinContent(irec,yieldMC*hResponse->GetBinContent(igen,irec)); | |
1119 | sumYield+=hResponse->GetBinContent(igen,irec); | |
1120 | } | |
ada4df5e | 1121 | // cout << "igen: " << igen << "\tpTMC: " << pTMC << "\tsumYield: " << sumYield << endl; |
1122 | // cout << "yieldMC: " << yieldMC << endl; | |
1123 | ||
31d3e4f0 | 1124 | } |
1125 | ||
1126 | return hRec; | |
1127 | } | |
1128 | ||
1129 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
1130 | Double_t AliAnaChargedJetResponseMaker::InterpolateFast(TGraph *gr, Double_t x) { | |
1131 | ||
1132 | Double_t ipmax = gr->GetN()-1.; | |
1133 | Double_t x2,y2,xold,yold; | |
1134 | ||
1135 | Double_t xmin,ymin,xmax, ymax; | |
1136 | gr->GetPoint(0,xmin,ymin); | |
1137 | gr->GetPoint(gr->GetN()-1,xmax,ymax); | |
1138 | if(x<xmin || x>xmax) return 0; | |
1139 | ||
1140 | Double_t ip = ipmax/2.; | |
1141 | Double_t ipold = 0.; | |
1142 | gr->GetPoint((int)(ip),x2,y2); | |
1143 | ||
1144 | Int_t i = 0; | |
1145 | ||
1146 | if(x2>x) { | |
1147 | while(x2>x) { | |
1148 | if(i==0) ipold=0.; | |
1149 | ip -= (ip)/2.; | |
1150 | gr->GetPoint((int)(ip),x2,y2); | |
1151 | if(x2>x){} | |
1152 | else ipold = ip; | |
1153 | i++; | |
1154 | // cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl; | |
1155 | } | |
1156 | } | |
1157 | else if(x2<x) { | |
1158 | while(x2<x) { | |
1159 | if(i==0) ipold=ipmax; | |
1160 | ip += (ipold-ip)/2.; | |
1161 | gr->GetPoint((int)(ip),x2,y2); | |
1162 | if(x2>x) ipold = ip; | |
1163 | else {} | |
1164 | i++; | |
1165 | // cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl; | |
1166 | } | |
1167 | } | |
1168 | ||
1169 | int ip2 = 0; | |
1170 | if(x2>x) { | |
1171 | ip2 = (int)(ip)-1; | |
1172 | gr->GetPoint(ip2,x2,y2); | |
1173 | while(x2>x) { | |
1174 | ip2--; | |
1175 | gr->GetPoint(ip2,x2,y2); | |
1176 | } | |
1177 | gr->GetPoint(ip2+1,xold,yold); | |
1178 | } | |
1179 | else { | |
1180 | ip2 = (int)(ip)+1; | |
1181 | gr->GetPoint(ip2,x2,y2); | |
1182 | while(x2<x) { | |
1183 | ip2++; | |
1184 | gr->GetPoint(ip2,x2,y2); | |
1185 | } | |
1186 | gr->GetPoint(ip2-1,xold,yold); | |
1187 | ||
1188 | } | |
1189 | // cout << "For x=" << x << " interpolate between: " << xold << " and " << x2 << endl; | |
1190 | return ((x-xold)*y2 + (x2-x)*yold) / (x2-xold); | |
1191 | ||
1192 | } | |
1193 | ||
1194 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
1195 | Double_t AliAnaChargedJetResponseMaker::InterpolateFast(TH1 *h, Double_t x) { | |
1196 | ||
1197 | // Double_t ipmax = gr->GetN()-1.; | |
1198 | Double_t ipmax = (double)h->GetNbinsX(); | |
1199 | Double_t x2,y2,xold,yold; | |
1200 | ||
1201 | Double_t xmin = h->GetXaxis()->GetXmin(); | |
1202 | Double_t xmax = h->GetXaxis()->GetXmax(); | |
1203 | if(x<xmin || x>xmax) return 0; | |
1204 | ||
1205 | Double_t ip = ipmax/2.; | |
1206 | Double_t ipold = 0.; | |
1207 | ||
1208 | x2 = h->GetBinCenter((int)ip); | |
1209 | y2 = h->GetBinContent((int)ip); | |
1210 | ||
1211 | Int_t i = 0; | |
1212 | ||
1213 | if(x2>x) { | |
1214 | while(x2>x) { | |
1215 | if(i==0) ipold=0.; | |
1216 | ip -= (ip)/2.; | |
1217 | x2 = h->GetBinCenter((int)ip); | |
1218 | y2 = h->GetBinContent((int)ip); | |
1219 | if(x2>x) {} | |
1220 | else ipold = ip; | |
1221 | i++; | |
1222 | // cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl; | |
1223 | } | |
1224 | } | |
1225 | else if(x2<x) { | |
1226 | while(x2<x) { | |
1227 | if(i==0) ipold=ipmax; | |
1228 | ip += (ipold-ip)/2.; | |
1229 | x2 = h->GetBinCenter((int)ip); | |
1230 | y2 = h->GetBinContent((int)ip); | |
1231 | if(x2>x) ipold = ip; | |
1232 | else {} | |
1233 | i++; | |
1234 | // cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl; | |
1235 | } | |
1236 | } | |
1237 | ||
1238 | int ip2 = 0; | |
1239 | if(x2>x) { | |
1240 | ip2 = (int)(ip)-1; | |
1241 | x2 = h->GetBinCenter(ip2); | |
1242 | y2 = h->GetBinContent(ip2); | |
1243 | while(x2>x) { | |
1244 | ip2--; | |
1245 | x2 = h->GetBinCenter(ip2); | |
1246 | y2 = h->GetBinContent(ip2); | |
1247 | } | |
1248 | xold = h->GetBinCenter(ip2+1); | |
1249 | yold = h->GetBinContent(ip2+1); | |
1250 | } | |
1251 | else { | |
1252 | ip2 = (int)(ip)+1; | |
1253 | x2 = h->GetBinCenter(ip2); | |
1254 | y2 = h->GetBinContent(ip2); | |
1255 | while(x2<x) { | |
1256 | ip2++; | |
1257 | x2 = h->GetBinCenter(ip2); | |
1258 | y2 = h->GetBinContent(ip2); | |
1259 | } | |
1260 | xold = h->GetBinCenter(ip2-1); | |
1261 | yold = h->GetBinContent(ip2-1); | |
1262 | ||
1263 | } | |
1264 | // cout << "For x=" << x << " interpolate between: " << xold << " and " << x2 << endl; | |
1265 | return ((x-xold)*y2 + (x2-x)*yold) / (x2-xold); | |
1266 | ||
1267 | } | |
1268 | ||
1269 |