changes and additions for the response maker task. fix visualisation (M. Verweij)
[u/mrichter/AliRoot.git] / PWGJE / AliAnaChargedJetResponseMaker.cxx
CommitLineData
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"
12
13#define round(x) ((x)>=0?(long)((x)+0.5):(long)((x)-0.5))
14
c64cb1f6 15using std::cout;
16using std::endl;
17
31d3e4f0 18ClassImp(AliAnaChargedJetResponseMaker)
19
20AliAnaChargedJetResponseMaker::AliAnaChargedJetResponseMaker():
21 fDebug(kFALSE),
22 fResolutionType(kParam),
23 fDeltaPt(0x0),
24 fhDeltaPt(0x0),
25 fDimensions(1),
26 fDimRec(0),
27 fDimGen(1),
28 fPtMin(-999),
29 fPtMax(-999),
30 fNbins(0),
31 fBinArrayPtRec(0x0),
32 fPtMeasured(0x0),
33 fEffFlat(0),
34 fEfficiency(0x0),
35 fEfficiencyFine(0x0),
36 fResponseMatrix(0x0),
37 fResponseMatrixFine(0x0),
38 fPtMinUnfolded(0.),
39 fPtMaxUnfolded(0.),
40 fPtMaxUnfoldedHigh(-1.),
41 fBinWidthFactorUnfolded(2),
42 fSkipBinsUnfolded(0),
43 fExtraBinsUnfolded(5),
44 fbVariableBinning(kFALSE),
45 fPtMaxUnfVarBinning(0),
46 f1MergeFunction(0x0),
47 fFineFrac(10),
48 fbCalcErrors(kFALSE)
49{;}
50
31d3e4f0 51
fa7c34ba 52//--------------------------------------------------------------------------------------------------------------------------------------------------
53AliAnaChargedJetResponseMaker::AliAnaChargedJetResponseMaker(const AliAnaChargedJetResponseMaker& obj):
54 fDebug(obj.fDebug),
55 fResolutionType(obj.fResolutionType),
56 fDeltaPt(obj.fDeltaPt),
57 fhDeltaPt(obj.fhDeltaPt),
58 fDimensions(obj.fDimensions),
59 fDimRec(obj.fDimRec),
60 fDimGen(obj.fDimGen),
61 fPtMin(obj.fPtMin),
62 fPtMax(obj.fPtMax),
63 fNbins(obj.fNbins),
64 fBinArrayPtRec(obj.fBinArrayPtRec),
65 fPtMeasured(obj.fPtMeasured),
66 fEffFlat(obj.fEffFlat),
67 fEfficiency(obj.fEfficiency),
68 fEfficiencyFine(obj.fEfficiencyFine),
69 fResponseMatrix(obj.fResponseMatrix),
70 fResponseMatrixFine(obj.fResponseMatrixFine),
71 fPtMinUnfolded(obj.fPtMinUnfolded),
72 fPtMaxUnfolded(obj.fPtMaxUnfolded),
73 fPtMaxUnfoldedHigh(obj.fPtMaxUnfoldedHigh),
74 fBinWidthFactorUnfolded(obj.fBinWidthFactorUnfolded),
75 fSkipBinsUnfolded(obj.fSkipBinsUnfolded),
76 fExtraBinsUnfolded(obj.fExtraBinsUnfolded),
77 fbVariableBinning(obj.fbVariableBinning),
78 fPtMaxUnfVarBinning(obj.fPtMaxUnfVarBinning),
79 f1MergeFunction(obj.f1MergeFunction),
80 fFineFrac(obj.fFineFrac),
81 fbCalcErrors(obj.fbCalcErrors)
82{;}
83
84//--------------------------------------------------------------------------------------------------------------------------------------------------
85AliAnaChargedJetResponseMaker& AliAnaChargedJetResponseMaker::operator=(const AliAnaChargedJetResponseMaker& other)
86{
87// Assignment
88 if(&other == this) return *this;
89 AliAnaChargedJetResponseMaker::operator=(other);
90 fDebug = other.fDebug;
91 fResolutionType = other.fResolutionType;
92 fDeltaPt = other.fDeltaPt;
93 fhDeltaPt = other.fhDeltaPt;
94 fDimensions = other.fDimensions;
95 fDimRec = other.fDimRec;
96 fDimGen = other.fDimGen;
97 fPtMin = other.fPtMin;
98 fPtMax = other.fPtMax;
99 fNbins = other.fNbins;
100 fBinArrayPtRec = other.fBinArrayPtRec;
101 fPtMeasured = other.fPtMeasured;
102 fEffFlat = other.fEffFlat;
103 fEfficiency = other.fEfficiency;
104 fEfficiencyFine = other.fEfficiencyFine;
105 fResponseMatrix = other.fResponseMatrix;
106 fResponseMatrixFine = other.fResponseMatrixFine;
107 fPtMinUnfolded = other.fPtMinUnfolded;
108 fPtMaxUnfolded = other.fPtMaxUnfolded;
109 fPtMaxUnfoldedHigh = other.fPtMaxUnfoldedHigh;
110 fBinWidthFactorUnfolded = other.fBinWidthFactorUnfolded;
111 fSkipBinsUnfolded = other.fSkipBinsUnfolded;
112 fExtraBinsUnfolded = other.fExtraBinsUnfolded;
113 fbVariableBinning = other.fbVariableBinning;
114 fPtMaxUnfVarBinning = other.fPtMaxUnfVarBinning;
115 f1MergeFunction = other.f1MergeFunction;
116 fFineFrac = other.fFineFrac;
117 fbCalcErrors = other.fbCalcErrors;
118
119 return *this;
31d3e4f0 120}
121
122//--------------------------------------------------------------------------------------------------------------------------------------------------
123void AliAnaChargedJetResponseMaker::SetMeasuredSpectrum(TH1D *hPtMeasured)
124{
125 //
126 // Set measured spectrum in THnSparse format
127 // This defines the minimum and maximum pT on the reconstructed axis of the response matrix
128 //
129 if(fDebug) printf(">>AliAnaChargedJetResponseMaker::SetMeasuredSpectrum \n");
130
131 fNbins = hPtMeasured->GetXaxis()->GetNbins();
132 fPtMin = hPtMeasured->GetXaxis()->GetXmin();
133 fPtMax = hPtMeasured->GetXaxis()->GetXmax();
134
135 if(fDebug) printf("fNbins: %d fPtMin: %f fPtMax: %f \n",fNbins,fPtMin,fPtMax);
136
137 if(fBinArrayPtRec) delete fBinArrayPtRec;
138 fBinArrayPtRec = new Double_t[fNbins+1];
139 for(int j = 0; j<fNbins; j++) {
140 fBinArrayPtRec[j] = hPtMeasured->GetXaxis()->GetBinLowEdge(j+1);
141 }
142 fBinArrayPtRec[fNbins] = hPtMeasured->GetXaxis()->GetBinUpEdge(fNbins);
143
144
145 Int_t nbins[fDimensions];
146 Double_t xmin[fDimensions];
147 Double_t xmax[fDimensions];
148 for(int dim = 0; dim<fDimensions; dim++) {
149 nbins[dim] = fNbins;
150 xmin[dim] = fPtMin;
151 xmax[dim] = fPtMax;
152 }
153
154 if(fPtMeasured) delete fPtMeasured;
155 fPtMeasured = new THnSparseD("fPtMeasured","Measured pT spectrum; p_{T}^{rec}",fDimensions,nbins,xmin,xmax);
156 fPtMeasured->Sumw2();
157 fPtMeasured->GetAxis(0)->SetTitle("p_{T}^{rec}");
158 fPtMeasured->SetBinEdges(0,fBinArrayPtRec);
159
160 //Fill
161 if(fDebug) printf("fill measured THnSparse\n");
162 if(fNbins!=hPtMeasured->GetNbinsX())
163 printf("WARNING: nbins not correct \t %d vs %d !!!\n",fNbins,hPtMeasured->GetNbinsX());
164
165 int bin[1] = {0};
166 bin[0] = 0;
167 for(int i = hPtMeasured->FindBin(fPtMin); i<hPtMeasured->FindBin(fPtMax); i++) {
168 double pT[1];
169 pT[0]= hPtMeasured->GetBinCenter(i);
170 fPtMeasured->SetBinContent(bin,(Double_t)hPtMeasured->GetBinContent(i));
171 fPtMeasured->SetBinError(bin,(Double_t)hPtMeasured->GetBinError(i));
172 bin[0]++;
173 }
174
175 if(fDebug) printf("fPtMeasured->GetNbins(): %d \n",(int)fPtMeasured->GetNbins());
176
177}
178
179//--------------------------------------------------------------------------------------------------------------------------------------------------
180void AliAnaChargedJetResponseMaker::SetFlatEfficiency(Double_t eff) {
181
ef62323a 182 //
183 // Set flat efficiency to value of eff
184 //
185
31d3e4f0 186 fEffFlat = eff;
ef62323a 187 return;
31d3e4f0 188
189 Int_t nbins[fDimensions];
190 Double_t xmin[fDimensions];
191 Double_t xmax[fDimensions];
192 for(int dim = 0; dim<fDimensions; dim++) {
193 nbins[dim] = fNbins;
194 xmin[dim] = fPtMin;
195 xmax[dim] = fPtMax;
196 }
197
198 if(fEfficiency) delete fEfficiency;
199 fEfficiency = new THnSparseD("fEfficiency","Efficiency - no resolution effects; p_{T}^{gen}",fDimensions,nbins,xmin,xmax);
200 fEfficiency->Sumw2();
201 fEfficiency->GetAxis(0)->SetTitle("p_{T}^{gen}");
202 fEfficiency->SetBinEdges(0,fBinArrayPtRec);
203
204 for(int i=0; i<fNbins; i++) {
205 int bin[1];
206 bin[0] = i;
207 fEfficiency->SetBinContent(bin,fEffFlat);
208 fEfficiency->SetBinError(bin,0);
209 }
210
211}
212
213//--------------------------------------------------------------------------------------------------------------------------------------------------
214void AliAnaChargedJetResponseMaker::SetEfficiency(TGraphErrors *grEff)
215{
ef62323a 216 //
217 // Fill fEfficiency (THnSparse) with values from grEff
218 //
219
31d3e4f0 220 Int_t nbins[fDimensions];
221 Double_t xmin[fDimensions];
222 Double_t xmax[fDimensions];
223 for(int dim = 0; dim<fDimensions; dim++) {
224 nbins[dim] = fNbins;
225 xmin[dim] = fPtMin;
226 xmax[dim] = fPtMax;
227 }
228
229 if(fEfficiency) delete fEfficiency;
230 fEfficiency = new THnSparseD("fEfficiency","Efficiency - no resolution effects; p_{T}^{gen}",fDimensions,nbins,xmin,xmax);
231 fEfficiency->Sumw2();
232 fEfficiency->GetAxis(0)->SetTitle("p_{T}^{gen}");
233 fEfficiency->SetBinEdges(0,fBinArrayPtRec);
234
235 double pT[1];
236 double yield = 0.;
237 double error = 0.;
238 double dummy = 0.;
239 int nbinsgr = grEff->GetN();
240
241 for(int i=0; i<nbinsgr; i++) {
242 grEff->GetPoint(i,dummy,yield);
243 pT[0] = dummy;
244 error = grEff->GetErrorY(i);
245
246 fEfficiency->Fill(pT,yield);
247 fEfficiency->SetBinError(i,error);
248
249 }
250
251}
252
253//--------------------------------------------------------------------------------------------------------------------------------------------------
254void AliAnaChargedJetResponseMaker::MakeResponseMatrixJetsFineMerged(Int_t skipBins, Int_t binWidthFactor, Int_t extraBins, Bool_t bVariableBinning, Double_t ptmax)
255{
256 //
257 // Make jet response matrix
258 //
259
ef62323a 260 if(fDebug) printf(">>AliAnaChargedJetResponseMaker::MakeResponseMatrixJetsFineMerged\n");
261
262 if(!fPtMeasured) {
263 if(fDebug) printf("fPtMeasured does not exist. Aborting!!\n");
264 return;
265 }
31d3e4f0 266 if(fResponseMatrix) { delete fResponseMatrix; }
267 if(fResponseMatrixFine) { delete fResponseMatrixFine; }
268
269 SetSkipBinsUnfolded(skipBins);
270 SetBinWidthFactorUnfolded(binWidthFactor);
271 SetExtraBinsUnfolded(extraBins);
272 SetVariableBinning(bVariableBinning,ptmax);
273
274 InitializeResponseMatrix();
275 InitializeEfficiency();
276
277 InitializeResponseMatrixFine();
278 InitializeEfficiencyFine();
279
280 FillResponseMatrixFineAndMerge();
281
282}
283
284//--------------------------------------------------------------------------------------------------------------------------------------------------
285void AliAnaChargedJetResponseMaker::InitializeResponseMatrix() {
286 //
ef62323a 287 //Define bin edges of RM to be used for unfolding
31d3e4f0 288 //
289
290 Int_t nbins[fDimensions*2];
291 nbins[fDimRec] = fNbins;
fa7c34ba 292 nbins[fDimGen] = fNbins;
31d3e4f0 293
294 double binWidthMeas = (double)((fPtMax-fPtMin)/fNbins);
295 double binWidthUnf = (double)fBinWidthFactorUnfolded*binWidthMeas;
296 double binWidthUnfLowPt = -1.;
297 if(fbVariableBinning)
298 binWidthUnfLowPt = binWidthUnf*0.5;
299
300 if(fExtraBinsUnfolded>0) {
301 fPtMaxUnfolded = fPtMax+(double)(fExtraBinsUnfolded)*binWidthUnf;
302 nbins[fDimGen]+=fExtraBinsUnfolded;
303 }
304
305 printf("fPtMinMeas: %f fPtMaxMeas: %f\n",fPtMin,fPtMax);
306 printf("binWidthMeas: %f binWidthUnf: %f fBinWidthFactorUnfolded: %d\n",binWidthMeas,binWidthUnf,fBinWidthFactorUnfolded);
307 printf("binWidthUnfLowPt: %f\n",binWidthUnfLowPt);
308
ef62323a 309 int tmp = round((fPtMaxUnfolded-fPtMinUnfolded)/binWidthUnf); //nbins which fit between 0 and fPtMaxUnfolded
31d3e4f0 310 tmp = tmp - fSkipBinsUnfolded;
311 fPtMinUnfolded = fPtMaxUnfolded-(double)(tmp)*binWidthUnf; //set ptmin unfolded
312 fPtMaxUnfolded = fPtMinUnfolded+(double)(tmp)*binWidthUnf; //set ptmax unfolded
313 nbins[fDimGen] = (int)((fPtMaxUnfolded-fPtMinUnfolded)/binWidthUnf); //adjust nbins to bin width
314
315 if(fbVariableBinning) {
316 tmp = (int)(fPtMaxUnfVarBinning/binWidthUnfLowPt);
317 tmp = tmp - fSkipBinsUnfolded;
318 fPtMinUnfolded = fPtMaxUnfVarBinning-(double)(tmp)*binWidthUnfLowPt;
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 }
322
323 Double_t binWidth[2];
324 binWidth[fDimRec] = (double)((fPtMax-fPtMin)/nbins[fDimRec]);
325 binWidth[fDimGen] = (double)((fPtMaxUnfolded-fPtMinUnfolded)/nbins[fDimGen]);
326
327 Double_t xmin[fDimensions*2];
328 Double_t xmax[fDimensions*2];
329 xmin[fDimRec] = fPtMin;
330 xmax[fDimRec] = fPtMax;
331 xmin[fDimGen] = fPtMinUnfolded;
332 xmax[fDimGen] = fPtMaxUnfolded;
333
334 printf("xmin[fDimRec]: %f xmin[fDimGen]: %f \n",xmin[fDimRec],xmin[fDimGen]);
335 printf("xmax[fDimRec]: %f xmax[fDimGen]: %f \n",xmax[fDimRec],xmax[fDimGen]);
336 printf("nbins[fDimRec]: %d nbins[fDimGen]: %d \n",nbins[fDimRec],nbins[fDimGen]);
337 if(!fbVariableBinning) printf("binWidth[fDimRec]: %f binWidth[fDimGen]: %f \n",binWidth[fDimRec],binWidth[fDimGen]);
338
339 Double_t binArrayPt0[nbins[fDimRec]+1];
340 Double_t binArrayPt1[nbins[fDimGen]+1];
341 Double_t xmaxGen = TMath::Max(xmax[fDimGen],fPtMaxUnfoldedHigh);
342
343 //Define bin limits reconstructed/measured axis
344 for(int i=0; i<nbins[fDimRec]; i++) {
345 binArrayPt0[i] = xmin[fDimRec]+(double)i*binWidth[fDimRec];
346 }
347 binArrayPt0[nbins[fDimRec]]= xmax[fDimRec];
348
349 //Define bin limits generated/unfolded axis
fa7c34ba 350 binArrayPt1[0] = fPtMinUnfolded;
351 for(int i=1; i<nbins[fDimGen]; i++) {
31d3e4f0 352 if(fbVariableBinning) {
353 double test = xmin[fDimGen]+(double)i*binWidthUnfLowPt;
354 if(test<=fPtMaxUnfVarBinning) binArrayPt1[i] = test;
355 else binArrayPt1[i] = binArrayPt1[i-1]+binWidthUnf;
356 }
fa7c34ba 357 else
358 binArrayPt1[i] = xmin[fDimGen]+(double)i*binWidth[fDimGen];
31d3e4f0 359 //printf("RM. i = %d \t binArrayPt[i] = %f \n",i,binArrayPt1[i]);
360 }
361 binArrayPt1[nbins[fDimGen]]= xmaxGen;
362
363
364 // Response matrix : dimensions must be 2N = 2 x (number of variables)
365 // dimensions 0 -> N-1 must be filled with reconstructed values
366 // dimensions N -> 2N-1 must be filled with generated values
ef62323a 367 if(fResponseMatrix) delete fResponseMatrix;
31d3e4f0 368 fResponseMatrix = new THnSparseD("fResponseMatrix","Response Matrix pTMC vs pTrec",fDimensions*2,nbins,xmin,xmax);
369 fResponseMatrix->Sumw2();
370 fResponseMatrix->GetAxis(fDimRec)->SetTitle("p_{T}^{rec}");
371 fResponseMatrix->GetAxis(fDimGen)->SetTitle("p_{T}^{gen}");
372
373 fResponseMatrix->SetBinEdges(fDimRec,binArrayPt0);
374 fResponseMatrix->SetBinEdges(fDimGen,binArrayPt1);
375
ef62323a 376 Int_t bin[2] = {1,1};
377 for(int i=1; i<fResponseMatrix->GetAxis(0)->GetNbins(); i++) {
378 bin[0]=i;
379 for(int j=1; j<fResponseMatrix->GetAxis(1)->GetNbins(); j++) {
380 bin[1]=j;
381 fResponseMatrix->SetBinContent(bin,0.);
382 }
383 }
384
385
31d3e4f0 386
387}
388
389//--------------------------------------------------------------------------------------------------------------------------------------------------
390void AliAnaChargedJetResponseMaker::InitializeEfficiency() {
391 //
392 // Define dimensions of efficiency THnSparse
393 //
394
395 if(!fResponseMatrix) {
396 printf("AliAnaChargedJetResponseMaker::InitializeEfficiency()\n");
397 printf("Can not define dimensions efficiency without fResponseMatrix dimensions. Aborting \n");
398 return;
399 }
400
401 TAxis *genAxis = fResponseMatrix->GetAxis(fDimGen);
402
403 Int_t nbinsEff[fDimensions];
404 Double_t xminEff[fDimensions];
405 Double_t xmaxEff[fDimensions];
406
407 for(int dim = 0; dim<fDimensions; dim++) {
408 nbinsEff[dim] = genAxis->GetNbins();
409 xminEff[dim] = genAxis->GetXmin();
410 xmaxEff[dim] = genAxis->GetXmax();
411 }
412
413 if(fEfficiency) delete fEfficiency;
414 fEfficiency = new THnSparseD("fEfficiency","Efficiency - no resolution effects; p_{T}^{gen}",fDimensions,nbinsEff,xminEff,xmaxEff);
415 fEfficiency->Sumw2();
416 fEfficiency->GetAxis(0)->SetTitle("p_{T}^{gen}");
417
418 const Double_t *binArrayPt = genAxis->GetXbins()->GetArray();
419 fEfficiency->SetBinEdges(0,binArrayPt);
420
421}
422
423//--------------------------------------------------------------------------------------------------------------------------------------------------
424void AliAnaChargedJetResponseMaker::InitializeResponseMatrixFine() {
425 //
426 //Initialize fine response matrix
427 //
428
429 Int_t nbinsFine[fDimensions*2];
fa7c34ba 430 Double_t xminFine[fDimensions*2];
431 Double_t xmaxFine[fDimensions*2];
31d3e4f0 432 Double_t pTarrayFine[fDimensions*2];
433
fa7c34ba 434 nbinsFine[fDimRec] = fResponseMatrix->GetAxis(fDimRec)->GetNbins()*fFineFrac;
435 nbinsFine[fDimGen] = fResponseMatrix->GetAxis(fDimRec)->GetNbins()*fFineFrac;
436 xminFine[fDimRec] = TMath::Min(fPtMin,0.);
437 xminFine[fDimGen] = TMath::Min(fPtMin,0.);
438 xmaxFine[fDimRec] = fResponseMatrix->GetAxis(fDimGen)->GetXmax()+40.;
439 xmaxFine[fDimGen] = fResponseMatrix->GetAxis(fDimGen)->GetXmax()+40.;
440 pTarrayFine[fDimRec] = 0.;
441 pTarrayFine[fDimGen] = 0.;
31d3e4f0 442
443 Double_t binWidth[2];
444 binWidth[fDimRec] = fResponseMatrix->GetAxis(fDimRec)->GetBinWidth(1);
445
446 Double_t binWidthFine[2];
447 binWidthFine[fDimRec] = binWidth[fDimRec]/((double)fFineFrac);
448 binWidthFine[fDimGen] = binWidthFine[fDimRec]*(double)fBinWidthFactorUnfolded;
449
450 nbinsFine[fDimRec] = (int)((xmaxFine[fDimRec]-xminFine[fDimRec])/binWidthFine[fDimRec]); //adjust nbins to bin width
451 nbinsFine[fDimGen] = (int)((xmaxFine[fDimGen]-xminFine[fDimGen])/binWidthFine[fDimGen]); //adjust nbins to bin width
452
453 printf("xminFine[fDimRec]: %f xminFine[fDimGen]: %f \n",xminFine[fDimRec],xminFine[fDimGen]);
454 printf("xmaxFine[fDimRec]: %f xmaxFine[fDimGen]: %f \n",xmaxFine[fDimRec],xmaxFine[fDimGen]);
455 printf("nbinsFine[fDimRec]: %d nbinsFine[fDimGen]: %d \n",nbinsFine[fDimRec],nbinsFine[fDimGen]);
456 printf("binWidthFine[fDimRec]: %f binWidthFine[fDimGen]: %f \n",binWidthFine[fDimRec],binWidthFine[fDimGen]);
457
458
459 Double_t binArrayPt0Fine[nbinsFine[fDimRec]+1];
460 Double_t binArrayPt1Fine[nbinsFine[fDimGen]+1];
461
462 for(int i=0; i<nbinsFine[fDimRec]; i++) {
463 binArrayPt0Fine[i] = xminFine[fDimRec]+(double)i*binWidthFine[fDimRec];
464 // printf("RM. i = %d \t binArrayPtFine[i] = %f \n",i,binArrayPt0Fine[i]);
465 }
466 binArrayPt0Fine[nbinsFine[fDimRec]]= xmaxFine[fDimRec];
467
468 for(int i=0; i<nbinsFine[fDimGen]; i++) {
469 binArrayPt1Fine[i] = xminFine[fDimGen]+(double)i*binWidthFine[fDimGen];
470 // printf("RM. i = %d \t binArrayPtFine[i] = %f \n",i,binArrayPt1Fine[i]);
471 }
472 binArrayPt1Fine[nbinsFine[fDimGen]]= xmaxFine[fDimGen];
473
474 // Response matrix : dimensions must be 2N = 2 x (number of variables)
475 // dimensions 0 -> N-1 must be filled with reconstructed values
476 // dimensions N -> 2N-1 must be filled with generated values
ef62323a 477 if(fResponseMatrixFine) delete fResponseMatrixFine;
31d3e4f0 478 fResponseMatrixFine = new THnSparseD("fResponseMatrixFine","Response Matrix pTMC vs pTrec",fDimensions*2,nbinsFine,xminFine,xmaxFine);
479 fResponseMatrixFine->Sumw2();
480 fResponseMatrixFine->GetAxis(fDimRec)->SetTitle("p_{T}^{rec}");
481 fResponseMatrixFine->GetAxis(fDimGen)->SetTitle("p_{T}^{gen}");
482
483 fResponseMatrixFine->SetBinEdges(fDimRec,binArrayPt0Fine);
484 fResponseMatrixFine->SetBinEdges(fDimGen,binArrayPt1Fine);
485
ef62323a 486 Int_t bin[2] = {1,1};
487 for(int i=1; i<fResponseMatrixFine->GetAxis(0)->GetNbins(); i++) {
488 bin[0]=i;
489 for(int j=1; j<fResponseMatrixFine->GetAxis(1)->GetNbins(); j++) {
490 bin[1]=j;
491 fResponseMatrixFine->SetBinContent(bin,0.);
492 }
493 }
494
31d3e4f0 495}
496
497//--------------------------------------------------------------------------------------------------------------------------------------------------
498void AliAnaChargedJetResponseMaker::InitializeEfficiencyFine() {
499 //
500 // Define dimensions of efficiency THnSparse
501 //
502
503 if(!fResponseMatrixFine) {
504 printf("AliAnaChargedJetResponseMaker::InitializeEfficiencyFine()\n");
505 printf("Can not define dimensions efficiency without fResponseMatrixFine dimensions. Aborting \n");
506 return;
507 }
508
509 TAxis *genAxis = fResponseMatrixFine->GetAxis(fDimGen);
510
511 Int_t nbinsEff[fDimensions];
512 Double_t xminEff[fDimensions];
513 Double_t xmaxEff[fDimensions];
514
515 for(int dim = 0; dim<fDimensions; dim++) {
516 nbinsEff[dim] = genAxis->GetNbins();
517 xminEff[dim] = genAxis->GetXmin();
518 xmaxEff[dim] = genAxis->GetXmax();
519 }
520
521 if(fEfficiencyFine) delete fEfficiencyFine;
522 fEfficiencyFine = new THnSparseD("fEfficiencyFine","EfficiencyFine - no resolution effects; p_{T}^{gen}",fDimensions,nbinsEff,xminEff,xmaxEff);
523 fEfficiencyFine->Sumw2();
524 fEfficiencyFine->GetAxis(0)->SetTitle("p_{T}^{gen}");
525
526 const Double_t *binArrayPt = genAxis->GetXbins()->GetArray();
527 fEfficiencyFine->SetBinEdges(0,binArrayPt);
528
529}
530
531//--------------------------------------------------------------------------------------------------------------------------------------------------
532void AliAnaChargedJetResponseMaker::FillResponseMatrixFineAndMerge() {
533 //
534 // Fill fine response matrix
535 //
536
ef62323a 537 if(!fResponseMatrix) {
538 printf("Dimensions of fResponseMatrix have to be defined first. Aborting!");
539 return;
540 }
541
542 if(!fResponseMatrixFine) {
543 printf("Dimensions of fResponseMatrixFine have to be defined first. Aborting!");
544 return;
545 }
546
31d3e4f0 547 TAxis *genAxis = fResponseMatrixFine->GetAxis(fDimGen);
548 TAxis *recAxis = fResponseMatrixFine->GetAxis(fDimRec);
549
550 Int_t nbinsFine[fDimensions*2];
551 Double_t xminFine[fDimensions*2];
552 Double_t xmaxFine[fDimensions*2];
553 Double_t pTarrayFine[fDimensions*2];
554
555 nbinsFine[fDimGen] = genAxis->GetNbins();
556 nbinsFine[fDimRec] = recAxis->GetNbins();
557 xminFine[fDimGen] = genAxis->GetXmin();
558 xminFine[fDimRec] = recAxis->GetXmin();
559 xmaxFine[fDimGen] = genAxis->GetXmax();
560 xmaxFine[fDimRec] = recAxis->GetXmax();
fa7c34ba 561 pTarrayFine[fDimGen] = 0.;
562 pTarrayFine[fDimRec] = 0.;
31d3e4f0 563
564 double sumyield = 0.;
565 double sumyielderror2 = 0.;
566
567 int ipt[2] = {0.,0.};
568 int iptMerged[2] = {0.,0.};
569 int ierror[2] = {0.,0.};
570
571 Int_t check = 0;
572 Double_t pTgen, pTrec;
573 Double_t yield = 0.;
574 Double_t error = 0.;
575
576 const int nng = fResponseMatrix->GetAxis(fDimGen)->GetNbins();
577 const int nnr = fResponseMatrix->GetAxis(fDimRec)->GetNbins();
578 Double_t errorArray[nng][nnr];
579 for(int iig =0; iig<nng; iig++) {
580 for(int iir =0; iir<nnr; iir++) {
581 errorArray[iig][iir] = 0.;
582 }
583 }
584
585 for(int iy=1; iy<=nbinsFine[fDimGen]; iy++) { //gen
586 pTgen = fResponseMatrixFine->GetAxis(fDimGen)->GetBinCenter(iy);
587 pTarrayFine[fDimGen] = pTgen;
588 ierror[fDimGen]=iy;
589 sumyield = 0.;
590 check = 0;
591
592 for(int ix=1; ix<=nbinsFine[fDimRec]; ix++) { //rec
593 pTrec = fResponseMatrixFine->GetAxis(fDimRec)->GetBinCenter(ix);
594 Double_t width = fResponseMatrixFine->GetAxis(fDimRec)->GetBinWidth(ix);
595 if(fResolutionType==kParam) {
596 yield = fDeltaPt->Eval(pTrec-pTgen);
597 error = 0.;
598 }
599 else if(fResolutionType==kResiduals) {
600 yield = fhDeltaPt->Interpolate(pTrec-pTgen);
601 error = 0.;
602 }
603 else if(fResolutionType==kResidualsErr) {
604 Double_t deltaPt = pTrec-pTgen;
605 int bin = fhDeltaPt->FindBin(deltaPt);
606 yield = fhDeltaPt->GetBinContent(bin);
607 error = fhDeltaPt->GetBinError(bin);
608 }
609 yield=yield*width;
610 error=error*width;
611 //avoid trouble with empty bins in the high pT tail of deltaPt distribution
612 if(check==0 && yield>0. && pTrec>pTgen) check=1;
613 if(check==1 && yield==0.) ix=nbinsFine[fDimRec];
614
615 sumyield+=yield;
616 sumyielderror2 += error*error;
617
618 pTarrayFine[fDimRec] = pTrec;
619 ierror[fDimRec] = ix;
620 fResponseMatrixFine->Fill(pTarrayFine,yield);
621 if(fbCalcErrors) fResponseMatrixFine->SetBinError(ierror,error);
622
623 }// ix (dimRec) loop
624
625 //Normalize to total number of counts =1
626 if(sumyield>1) {
627 ipt[fDimGen]=iy;
628 for(int ix=1; ix<=nbinsFine[fDimRec]; ix++) {
629 ipt[fDimRec]=ix;
630 fResponseMatrixFine->SetBinContent(ipt,fResponseMatrixFine->GetBinContent(ipt)/sumyield);
631 if(fResolutionType==kResidualsErr && fbCalcErrors) {
632 Double_t A = 1./sumyield*fResponseMatrix->GetBinError(ipt);
633 Double_t B = -1.*fResponseMatrix->GetBinContent(ipt)/sumyield/sumyield*TMath::Sqrt(sumyielderror2);
634 Double_t tmp2 = A*A + B*B;
635 fResponseMatrix->SetBinError(ipt,TMath::Sqrt(tmp2));
636 }
637
638 }
639 }
640
641 int bin[1];
642 bin[0] = iy;
643 fEfficiencyFine->SetBinContent(bin,sumyield);
644 if(fbCalcErrors) fEfficiencyFine->SetBinError(bin,TMath::Sqrt(sumyielderror2));
645
646 //fill merged response matrix
647
648 //find bin in fine RM correspoinding to mimimum/maximum bin of merged RM on rec axis
649 int ixMin = fResponseMatrixFine->GetAxis(fDimRec)->FindBin(fResponseMatrix->GetAxis(fDimRec)->GetXmin());
650 int ixMax = fResponseMatrixFine->GetAxis(fDimRec)->FindBin(fResponseMatrix->GetAxis(fDimRec)->GetXmax());
651
652 if(fResponseMatrixFine->GetAxis(fDimGen)->GetBinLowEdge(iy) >= fResponseMatrix->GetAxis(fDimGen)->GetXmin()) {
653 ipt[fDimGen]=iy;
654 iptMerged[fDimGen]=fResponseMatrix->GetAxis(fDimGen)->FindBin(pTgen);
655
656 Double_t weight = 1.;
657 if(f1MergeFunction) {
658 Double_t loEdge = fResponseMatrix->GetAxis(fDimGen)->GetBinLowEdge(iptMerged[fDimGen]);
659 Double_t upEdge = fResponseMatrix->GetAxis(fDimGen)->GetBinUpEdge(iptMerged[fDimGen]);
660 Float_t powInteg = f1MergeFunction->Integral(loEdge,upEdge);
661 //printf("loEdge = %f upEdge = %f powInteg = %f\n",loEdge,upEdge,powInteg);
662 if(powInteg>0.)
663 weight = f1MergeFunction->Integral(fResponseMatrixFine->GetAxis(fDimGen)->GetBinLowEdge(iy),fResponseMatrixFine->GetAxis(fDimGen)->GetBinUpEdge(iy))/powInteg;
664 // printf("weight: %f \n", weight );
665 } else {
666 weight = 1./((double)fFineFrac);
667 }
668
669 for(int ix=ixMin; ix<=ixMax; ix++) { //loop reconstructed axis
670 pTrec = fResponseMatrixFine->GetAxis(fDimRec)->GetBinCenter(ix);
671 ipt[fDimRec]=ix;
672 iptMerged[fDimRec]=fResponseMatrix->GetAxis(fDimRec)->FindBin(pTrec);
673
674 fResponseMatrix->AddBinContent(iptMerged,fResponseMatrixFine->GetBinContent(ipt)*weight);
675 if(fbCalcErrors) errorArray[iptMerged[fDimGen]-1][iptMerged[fDimRec]-1] += fResponseMatrixFine->GetBinError(ipt)*fResponseMatrixFine->GetBinError(ipt)*weight*weight;
676 }
677
678 }//loop gen axis fine response matrix
679
680 } // iy (dimGen) loop
681
682 //Fill Efficiency corresponding to merged response matrix
683 // FillEfficiency(errorArray,fFineFrac);
684
685 for(int iy=1; iy<=fResponseMatrix->GetAxis(fDimGen)->GetNbins(); iy++) { //gen
686 sumyield = 0.;
687 ipt[fDimGen] = iy;
688
689 for(int ix=1; ix<=fResponseMatrix->GetAxis(fDimRec)->GetNbins(); ix++) { //rec
690 ipt[fDimRec] = ix;
691 sumyield += fResponseMatrix->GetBinContent(ipt);
692
693 if(fbCalcErrors) fResponseMatrix->SetBinError(ipt,TMath::Sqrt(errorArray[iy][ix]));
694 }
695 printf("RM: pTgen: %f \t sumyield: %f \n",fResponseMatrix->GetAxis(fDimGen)->GetBinCenter(iy),sumyield);
696 int bin[1];
697 bin[0] = iy;
698 fEfficiency->SetBinContent(bin,sumyield);
699 if(fbCalcErrors) fEfficiency->SetBinError(bin,0);
700 }
701
702 if(fDebug) printf("fResponseMatrixFine->GetNbins(): %d \n",(int)fResponseMatrixFine->GetNbins());
703 if(fDebug) printf("fResponseMatrix->GetNbins(): %d \n",(int)fResponseMatrix->GetNbins());
704
705 if(fDebug) printf("Done constructing response matrix\n");
706
707}
708
ef62323a 709//--------------------------------------------------------------------------------------------------------------------------------------------------
710TH2* AliAnaChargedJetResponseMaker::MakeResponseMatrixRebin(TH2 *hRMFine, TH2 *hRM) {
711
712 //
713 // Rebin matrix hRMFine to dimensions of hRM
714 // function returns matrix in TH2D format with dimensions from hRM
715 //
716
717 TH2 *hRM2 = (TH2*)hRM->Clone("hRM2");
718 hRM2->Reset();
719
720 //Forst normalize columns of input
721 const Int_t nbinsNorm = hRM2->GetNbinsX();
722 cout << "nbinsNorm: " << nbinsNorm << endl;
723
724 TArrayF *normVector = new TArrayF(nbinsNorm);
725
726 for(int ix=1; ix<=hRM2->GetNbinsX(); ix++) {
727 Double_t xLow = hRM2->GetXaxis()->GetBinLowEdge(ix);
728 Double_t xUp = hRM2->GetXaxis()->GetBinUpEdge(ix);
729 //cout << "xLow: " << xLow << " xUp: " << xUp << "\t center: " << hRM2->GetXaxis()->GetBinCenter(ix) << endl;
730 Int_t binxLowFine = hRMFine->GetXaxis()->FindBin(xLow);
731 Int_t binxUpFine = hRMFine->GetXaxis()->FindBin(xUp)-1;
732 //cout << "xLowFine: " << hRMFine->GetXaxis()->GetBinLowEdge(binxLowFine) << "\txUpFine: " << hRMFine->GetXaxis()->GetBinUpEdge(binxUpFine) << endl;
733 normVector->SetAt(hRMFine->Integral(binxLowFine,binxUpFine,1,hRMFine->GetYaxis()->GetNbins()),ix-1);
734 }
735
736 Double_t content, oldcontent = 0.;
737 Int_t ixNew = 0;
738 Int_t iyNew = 0;
739 Double_t xval, yval;
740 Double_t xmin = hRM2->GetXaxis()->GetXmin();
741 Double_t ymin = hRM2->GetYaxis()->GetXmin();
742 Double_t xmax = hRM2->GetXaxis()->GetXmax();
743 Double_t ymax = hRM2->GetYaxis()->GetXmax();
744 for(int ix=1; ix<=hRMFine->GetXaxis()->GetNbins(); ix++) {
745 xval = hRMFine->GetXaxis()->GetBinCenter(ix);
746 if(xval<xmin || xval>xmax) continue;
747 ixNew = hRM2->GetXaxis()->FindBin(xval);
748 for(int iy=1; iy<=hRMFine->GetYaxis()->GetNbins(); iy++) {
749 yval = hRMFine->GetYaxis()->GetBinCenter(iy);
750 if(yval<ymin || yval>ymax) continue;
751 content = hRMFine->GetBinContent(ix,iy);
752 iyNew = hRM2->GetYaxis()->FindBin(yval);
753 oldcontent = hRM2->GetBinContent(ixNew,iyNew);
754
755 // if(fDebug) cout << "ixNew: " << ixNew << " " << xval << " iyNew: " << iyNew << " " << yval << " content: " << content << " oldContent: " << oldcontent << " newContent: " << oldcontent+content << endl;
756 Double_t weight = 1.;
757 if(normVector->At(ixNew-1)>0.) weight = 1./normVector->At(ixNew-1);
758 hRM2->SetBinContent(ixNew,iyNew,oldcontent+content*weight);
759 }
760 }
761
762 if(normVector) delete normVector;
763
764 return hRM2;
765
766}
767
31d3e4f0 768
769//--------------------------------------------------------------------------------------------------------------------------------------------------
770TH1D* AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TH1 *hGen, TH2 *hResponse,TH1 *hEfficiency,Bool_t bDrawSlices) {
771
772 //
773 // Multiply hGen with response matrix to obtain refolded spectrum
ef62323a 774 // Efficiency must be given. In case efficiency is 1 use SetFlatEfficiency(1.) before calling function
31d3e4f0 775 //
776
ef62323a 777 if(!hEfficiency) {
778 printf("Efficiency must be given. In case efficiency is 1 use SetFlatEfficiency(1.) before calling function. Aborting!");
779 return 0;
780 }
781
31d3e4f0 782 //For response
783 //x-axis: generated
784 //y-axis: reconstructed
785 if(fDebug>0) cout << "nbins hGen: " << hGen->GetNbinsX() << "\t nbins hResponseGen: " << hResponse->GetXaxis()->GetNbins() << "\t nbins hResponseRec: " << hResponse->GetYaxis()->GetNbins() << endl;
786
787 TH1D *hRec = hResponse->ProjectionY("hRec");
788 hRec->Sumw2();
789 hRec->Reset();
790 hRec->SetTitle("hRec");
791 hRec->SetName("hRec");
792
793 for(int irec=1; irec<=hRec->GetNbinsX(); irec++)
794 hRec->SetBinContent(irec,0);
795
796 TH1D *hRecGenBin = 0x0;
797 TCanvas *cSlices = 0x0;
798 if(bDrawSlices) {
799 cSlices = new TCanvas("cSlices","cSlices:Slices",400,400);
800 gPad->SetLogy();
801 }
802
803 Double_t yieldMC = 0.;
804 Double_t yieldMCerror = 0.;
805 Double_t sumYield = 0.;
806 const Int_t nbinsRec = hRec->GetNbinsX();
807 Double_t sumError2[nbinsRec+1];
31d3e4f0 808 Double_t eff = 0.;
809
810 for(int igen=1; igen<=hGen->GetNbinsX(); igen++) {
811 //get pTMC
812 sumYield = 0.;
813 if(fEffFlat>0.)
814 eff = fEffFlat;
815 else
816 eff = hEfficiency->GetBinContent(igen);
817 yieldMC = hGen->GetBinContent(igen)*eff;
818 // yieldMCerror = hGen->GetBinError(igen);
819
820 if(bDrawSlices) {
821 hRecGenBin = hResponse->ProjectionY(Form("hRecGenBin%d",igen));
822 hRecGenBin->Sumw2();
823 hRecGenBin->Reset();
824 hRecGenBin->SetTitle(Form("hRecGenBin%d",igen));
825 hRecGenBin->SetName(Form("hRecGenBin%d",igen));
826 }
827
828 for(int irec=1; irec<=hRec->GetNbinsX(); irec++) {
829 hRec->AddBinContent(irec,yieldMC*hResponse->GetBinContent(igen,irec));
830 sumYield+=hResponse->GetBinContent(igen,irec);
831 Double_t A = yieldMC*hResponse->GetBinError(igen,irec);
832 Double_t B = hResponse->GetBinContent(igen,irec)*yieldMCerror;
833 Double_t tmp2 = A*A + B*B;
834 sumError2[irec-1] += tmp2 ;
835
836 if(bDrawSlices)
837 hRecGenBin->SetBinContent(irec,yieldMC*hResponse->GetBinContent(igen,irec));
838
839 }
840 if(bDrawSlices) {
841 cSlices->cd();
842 hRecGenBin->SetLineColor(igen);
843 if(igen==1) hRecGenBin->DrawCopy();
844 else hRecGenBin->DrawCopy("same");
845 }
846
847 if(hRecGenBin) delete hRecGenBin;
848
849 cout << "igen: " << igen << "\tpTMC: " << hGen->GetXaxis()->GetBinCenter(igen) << "\teff:" << eff << "\tsumYield: " << sumYield << endl;
850 }
851
ef62323a 852 for(int i=0; i<=nbinsRec; i++) {
853 if(sumError2[i]>0.)
854 hRec->SetBinError(i+1,TMath::Sqrt(sumError2[i]));
855 }
31d3e4f0 856
857
858 return hRec;
859}
860
861//--------------------------------------------------------------------------------------------------------------------------------------------------
862TH1D* AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TF1 *fGen, TH2 *hResponse,TH1 *hEfficiency) {
863
ef62323a 864 //
865 // Multiply fGen function with response matrix to obtain (re)folded spectrum
866 // Efficiency must be given. In case efficiency is 1 use SetFlatEfficiency(1.) before calling function
867 //
868
31d3e4f0 869 //For response
870 //x-axis: generated
871 //y-axis: reconstructed
872
873 if(fDebug>0) printf(">>AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TF1 *fGen, TH2 *hResponse,TH1 *hEfficiency)");
874
ef62323a 875 if(!hEfficiency) {
876 printf("Efficiency must be given. In case efficiency is 1 use SetFlatEfficiency(1.) before calling function. Aborting!");
877 return 0;
878 }
879
31d3e4f0 880 TH1D *hRec = hResponse->ProjectionY("hRec");
881 hRec->Sumw2();
882 hRec->Reset();
883 hRec->SetTitle("hRec");
884 hRec->SetName("hRec");
885
886 // TH1D *hRec = new TH1D("hRec","hRec",hResponse->GetNbinsY(),hResponse->GetYaxis()->GetXmin(),hResponse->GetYaxis()->GetXmax());
887
888 for(int irec=1; irec<=hRec->GetNbinsX(); irec++)
889 hRec->SetBinContent(irec,0);
890
891 Double_t yieldMC = 0.;
892 Double_t sumYield = 0.;
893 Double_t eff = 0.;
894 for(int igen=1; igen<=hResponse->GetNbinsX(); igen++) {
895 //get pTMC
896 sumYield = 0.;
897 double pTMC = hResponse->GetXaxis()->GetBinCenter(igen);
898 int binEff = hEfficiency->FindBin(pTMC);
899 if(fEffFlat>0.)
900 eff = fEffFlat;
901 else
902 eff = hEfficiency->GetBinContent(binEff);
903 yieldMC = fGen->Eval(pTMC)*eff;
904 for(int irec=1; irec<=hResponse->GetNbinsY(); irec++) {
905 hRec->AddBinContent(irec,yieldMC*hResponse->GetBinContent(igen,irec));
906 sumYield+=hResponse->GetBinContent(igen,irec);
907 }
908 cout << "igen: " << igen << "\tpTMC: " << pTMC << "\tsumYield: " << sumYield << endl;
909 }
910
911 return hRec;
912}
913
914//--------------------------------------------------------------------------------------------------------------------------------------------------
915Double_t AliAnaChargedJetResponseMaker::InterpolateFast(TGraph *gr, Double_t x) {
916
917 Double_t ipmax = gr->GetN()-1.;
918 Double_t x2,y2,xold,yold;
919
920 Double_t xmin,ymin,xmax, ymax;
921 gr->GetPoint(0,xmin,ymin);
922 gr->GetPoint(gr->GetN()-1,xmax,ymax);
923 if(x<xmin || x>xmax) return 0;
924
925 Double_t ip = ipmax/2.;
926 Double_t ipold = 0.;
927 gr->GetPoint((int)(ip),x2,y2);
928
929 Int_t i = 0;
930
931 if(x2>x) {
932 while(x2>x) {
933 if(i==0) ipold=0.;
934 ip -= (ip)/2.;
935 gr->GetPoint((int)(ip),x2,y2);
936 if(x2>x){}
937 else ipold = ip;
938 i++;
939 // cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl;
940 }
941 }
942 else if(x2<x) {
943 while(x2<x) {
944 if(i==0) ipold=ipmax;
945 ip += (ipold-ip)/2.;
946 gr->GetPoint((int)(ip),x2,y2);
947 if(x2>x) ipold = ip;
948 else {}
949 i++;
950 // cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl;
951 }
952 }
953
954 int ip2 = 0;
955 if(x2>x) {
956 ip2 = (int)(ip)-1;
957 gr->GetPoint(ip2,x2,y2);
958 while(x2>x) {
959 ip2--;
960 gr->GetPoint(ip2,x2,y2);
961 }
962 gr->GetPoint(ip2+1,xold,yold);
963 }
964 else {
965 ip2 = (int)(ip)+1;
966 gr->GetPoint(ip2,x2,y2);
967 while(x2<x) {
968 ip2++;
969 gr->GetPoint(ip2,x2,y2);
970 }
971 gr->GetPoint(ip2-1,xold,yold);
972
973 }
974 // cout << "For x=" << x << " interpolate between: " << xold << " and " << x2 << endl;
975 return ((x-xold)*y2 + (x2-x)*yold) / (x2-xold);
976
977}
978
979//--------------------------------------------------------------------------------------------------------------------------------------------------
980Double_t AliAnaChargedJetResponseMaker::InterpolateFast(TH1 *h, Double_t x) {
981
982 // Double_t ipmax = gr->GetN()-1.;
983 Double_t ipmax = (double)h->GetNbinsX();
984 Double_t x2,y2,xold,yold;
985
986 Double_t xmin = h->GetXaxis()->GetXmin();
987 Double_t xmax = h->GetXaxis()->GetXmax();
988 if(x<xmin || x>xmax) return 0;
989
990 Double_t ip = ipmax/2.;
991 Double_t ipold = 0.;
992
993 x2 = h->GetBinCenter((int)ip);
994 y2 = h->GetBinContent((int)ip);
995
996 Int_t i = 0;
997
998 if(x2>x) {
999 while(x2>x) {
1000 if(i==0) ipold=0.;
1001 ip -= (ip)/2.;
1002 x2 = h->GetBinCenter((int)ip);
1003 y2 = h->GetBinContent((int)ip);
1004 if(x2>x) {}
1005 else ipold = ip;
1006 i++;
1007 // cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl;
1008 }
1009 }
1010 else if(x2<x) {
1011 while(x2<x) {
1012 if(i==0) ipold=ipmax;
1013 ip += (ipold-ip)/2.;
1014 x2 = h->GetBinCenter((int)ip);
1015 y2 = h->GetBinContent((int)ip);
1016 if(x2>x) ipold = ip;
1017 else {}
1018 i++;
1019 // cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl;
1020 }
1021 }
1022
1023 int ip2 = 0;
1024 if(x2>x) {
1025 ip2 = (int)(ip)-1;
1026 x2 = h->GetBinCenter(ip2);
1027 y2 = h->GetBinContent(ip2);
1028 while(x2>x) {
1029 ip2--;
1030 x2 = h->GetBinCenter(ip2);
1031 y2 = h->GetBinContent(ip2);
1032 }
1033 xold = h->GetBinCenter(ip2+1);
1034 yold = h->GetBinContent(ip2+1);
1035 }
1036 else {
1037 ip2 = (int)(ip)+1;
1038 x2 = h->GetBinCenter(ip2);
1039 y2 = h->GetBinContent(ip2);
1040 while(x2<x) {
1041 ip2++;
1042 x2 = h->GetBinCenter(ip2);
1043 y2 = h->GetBinContent(ip2);
1044 }
1045 xold = h->GetBinCenter(ip2-1);
1046 yold = h->GetBinContent(ip2-1);
1047
1048 }
1049 // cout << "For x=" << x << " interpolate between: " << xold << " and " << x2 << endl;
1050 return ((x-xold)*y2 + (x2-x)*yold) / (x2-xold);
1051
1052}
1053
1054