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