]> 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
standalone jet finder class (Ruediger)
[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
1053 //Initialize transposed response matrix h2RMTranspose
1054 TArrayD *arrayX = (TArrayD*)h2RM->GetXaxis()->GetXbins();
1055 for(int i=0; i<arrayX->GetSize(); i++) cout << "i: " << arrayX->At(i) << endl;
1056 Double_t *xbinsArray = arrayX->GetArray();
1057
1058 TArrayD *arrayY = (TArrayD*)h2RM->GetYaxis()->GetXbins();
1059 for(int i=0; i<arrayY->GetSize(); i++) cout << "i: " << arrayY->At(i) << endl;
1060 Double_t *ybinsArray = arrayY->GetArray();
1061
1062 TH2D *h2RMTranspose = new TH2D("h2RMTranspose","h2RMTranspose",h2RM->GetNbinsY(),ybinsArray,h2RM->GetNbinsX(),xbinsArray);
1063
1064 //Fill tranposed response matrix
1065 for (Int_t ibin = 1; ibin <= h2RMTranspose->GetNbinsX(); ibin++) {
1066 for (Int_t jbin = 1; jbin <= h2RMTranspose->GetNbinsY(); jbin++) {
1067 h2RMTranspose->SetBinContent(ibin,jbin, h2RM->GetBinContent(jbin,ibin));
1068 }
1069 }
1070
1071
1072 return h2RMTranspose;
1073
1074}
1075
1076//--------------------------------------------------------------------------------------------------------------------------------------------------
1077TH2* AliAnaChargedJetResponseMaker::NormalizeResponsMatrixYaxisWithPrior(TH2 *h2RM, TH1 *hPrior) {
1078 //
1079 // Normalize such that the Y projection is the prior
1080 //
1081
03372fd1 1082 // TH1D *hProjRespY = (TH1D*)h2RM->ProjectionY("hProjRespY");
5d87a047 1083 double intPrior = hPrior->Integral();//"width");
1084 for (Int_t jbin = 1; jbin <= h2RM->GetNbinsY(); jbin++) {
1085 // double corr = hPrior->GetBinContent(jbin)/hProjRespY->GetBinContent(jbin);
1086 for (Int_t ibin = 1; ibin <= h2RM->GetNbinsX(); ibin++) {
1087 double content = h2RM->GetBinContent(ibin,jbin);
1088 // h2RM->SetBinContent(ibin,jbin,content*corr);
1089 h2RM->SetBinContent(ibin,jbin,hPrior->GetBinContent(jbin)/hPrior->GetBinWidth(jbin)/intPrior*content);
1090 }
1091 }
1092
1093 return h2RM;
1094}
1095
31d3e4f0 1096//--------------------------------------------------------------------------------------------------------------------------------------------------
1097TH1D* AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TH1 *hGen, TH2 *hResponse,TH1 *hEfficiency,Bool_t bDrawSlices) {
1098
1099 //
1100 // Multiply hGen with response matrix to obtain refolded spectrum
ef62323a 1101 // Efficiency must be given. In case efficiency is 1 use SetFlatEfficiency(1.) before calling function
31d3e4f0 1102 //
1103
ef62323a 1104 if(!hEfficiency) {
5d87a047 1105 // printf("Efficiency must be given. In case efficiency is 1 use SetFlatEfficiency(1.) before calling function. Aborting!");
1106 // return 0;
1107 printf("Setting efficiency to 1 \n");
1108 hEfficiency = (TH1D*)hGen->Clone("hEfficiency");
1109 hEfficiency->Reset();
1110 for(int i=1; i<=hEfficiency->GetNbinsX(); i++) hEfficiency->SetBinContent(i,1.);
ef62323a 1111 }
1112
31d3e4f0 1113 //For response
1114 //x-axis: generated
1115 //y-axis: reconstructed
1116 if(fDebug>0) cout << "nbins hGen: " << hGen->GetNbinsX() << "\t nbins hResponseGen: " << hResponse->GetXaxis()->GetNbins() << "\t nbins hResponseRec: " << hResponse->GetYaxis()->GetNbins() << endl;
1117
1118 TH1D *hRec = hResponse->ProjectionY("hRec");
1119 hRec->Sumw2();
1120 hRec->Reset();
1121 hRec->SetTitle("hRec");
1122 hRec->SetName("hRec");
1123
1124 for(int irec=1; irec<=hRec->GetNbinsX(); irec++)
1125 hRec->SetBinContent(irec,0);
1126
1127 TH1D *hRecGenBin = 0x0;
1128 TCanvas *cSlices = 0x0;
1129 if(bDrawSlices) {
1130 cSlices = new TCanvas("cSlices","cSlices:Slices",400,400);
1131 gPad->SetLogy();
1132 }
1133
1134 Double_t yieldMC = 0.;
1135 Double_t yieldMCerror = 0.;
1136 Double_t sumYield = 0.;
1137 const Int_t nbinsRec = hRec->GetNbinsX();
1138 Double_t sumError2[nbinsRec+1];
0bf6381d 1139 for(int irec=0; irec<=hRec->GetNbinsX(); irec++)
1140 sumError2[irec]=0.;
31d3e4f0 1141 Double_t eff = 0.;
1142
1143 for(int igen=1; igen<=hGen->GetNbinsX(); igen++) {
1144 //get pTMC
1145 sumYield = 0.;
1146 if(fEffFlat>0.)
1147 eff = fEffFlat;
1148 else
1149 eff = hEfficiency->GetBinContent(igen);
1150 yieldMC = hGen->GetBinContent(igen)*eff;
adf3920d 1151 yieldMCerror = hGen->GetBinError(igen)*eff;
31d3e4f0 1152
1153 if(bDrawSlices) {
1154 hRecGenBin = hResponse->ProjectionY(Form("hRecGenBin%d",igen));
1155 hRecGenBin->Sumw2();
1156 hRecGenBin->Reset();
1157 hRecGenBin->SetTitle(Form("hRecGenBin%d",igen));
1158 hRecGenBin->SetName(Form("hRecGenBin%d",igen));
1159 }
1160
1161 for(int irec=1; irec<=hRec->GetNbinsX(); irec++) {
1162 hRec->AddBinContent(irec,yieldMC*hResponse->GetBinContent(igen,irec));
1163 sumYield+=hResponse->GetBinContent(igen,irec);
346cfe4a 1164 // Double_t A = yieldMC*hResponse->GetBinError(igen,irec);
31d3e4f0 1165 Double_t B = hResponse->GetBinContent(igen,irec)*yieldMCerror;
346cfe4a 1166 // Double_t tmp2 = A*A + B*B;
adf3920d 1167 //sumError2[irec-1] += tmp2 ;
1168 sumError2[irec-1] += B*B;
31d3e4f0 1169
1170 if(bDrawSlices)
1171 hRecGenBin->SetBinContent(irec,yieldMC*hResponse->GetBinContent(igen,irec));
1172
1173 }
1174 if(bDrawSlices) {
1175 cSlices->cd();
1176 hRecGenBin->SetLineColor(igen);
1177 if(igen==1) hRecGenBin->DrawCopy();
1178 else hRecGenBin->DrawCopy("same");
1179 }
1180
1181 if(hRecGenBin) delete hRecGenBin;
1182
1183 cout << "igen: " << igen << "\tpTMC: " << hGen->GetXaxis()->GetBinCenter(igen) << "\teff:" << eff << "\tsumYield: " << sumYield << endl;
1184 }
1185
ef62323a 1186 for(int i=0; i<=nbinsRec; i++) {
1187 if(sumError2[i]>0.)
1188 hRec->SetBinError(i+1,TMath::Sqrt(sumError2[i]));
1189 }
31d3e4f0 1190
1191
1192 return hRec;
1193}
1194
1195//--------------------------------------------------------------------------------------------------------------------------------------------------
1196TH1D* AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TF1 *fGen, TH2 *hResponse,TH1 *hEfficiency) {
1197
ef62323a 1198 //
1199 // Multiply fGen function with response matrix to obtain (re)folded spectrum
1200 // Efficiency must be given. In case efficiency is 1 use SetFlatEfficiency(1.) before calling function
1201 //
1202
31d3e4f0 1203 //For response
1204 //x-axis: generated
1205 //y-axis: reconstructed
1206
ada4df5e 1207 if(fDebug>0) printf(">>AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TF1 *fGen, TH2 *hResponse,TH1 *hEfficiency)\n");
ef62323a 1208
31d3e4f0 1209 TH1D *hRec = hResponse->ProjectionY("hRec");
1210 hRec->Sumw2();
1211 hRec->Reset();
1212 hRec->SetTitle("hRec");
1213 hRec->SetName("hRec");
1214
1215 // TH1D *hRec = new TH1D("hRec","hRec",hResponse->GetNbinsY(),hResponse->GetYaxis()->GetXmin(),hResponse->GetYaxis()->GetXmax());
1216
1217 for(int irec=1; irec<=hRec->GetNbinsX(); irec++)
1218 hRec->SetBinContent(irec,0);
1219
1220 Double_t yieldMC = 0.;
1221 Double_t sumYield = 0.;
1222 Double_t eff = 0.;
1223 for(int igen=1; igen<=hResponse->GetNbinsX(); igen++) {
1224 //get pTMC
1225 sumYield = 0.;
1226 double pTMC = hResponse->GetXaxis()->GetBinCenter(igen);
ada4df5e 1227 if(hEfficiency) {
1228 int binEff = hEfficiency->FindBin(pTMC);
31d3e4f0 1229 eff = hEfficiency->GetBinContent(binEff);
ada4df5e 1230 }
1231 else eff = 1.;
1232 // yieldMC = fGen->Eval(pTMC)*eff;
1233 yieldMC = fGen->Integral(hResponse->GetXaxis()->GetBinLowEdge(igen),hResponse->GetXaxis()->GetBinUpEdge(igen))*eff;
31d3e4f0 1234 for(int irec=1; irec<=hResponse->GetNbinsY(); irec++) {
1235 hRec->AddBinContent(irec,yieldMC*hResponse->GetBinContent(igen,irec));
1236 sumYield+=hResponse->GetBinContent(igen,irec);
1237 }
ada4df5e 1238 // cout << "igen: " << igen << "\tpTMC: " << pTMC << "\tsumYield: " << sumYield << endl;
1239 // cout << "yieldMC: " << yieldMC << endl;
1240
31d3e4f0 1241 }
1242
1243 return hRec;
1244}
1245
0bf6381d 1246//--------------------------------------------------------------------------------------------------------------------------------------------------
1247/* static */ Double_t AliAnaChargedJetResponseMaker::GetBetaPerDOFValue(Int_t betaColl, Int_t betaOpt) {
1248
1249 Float_t betaPerDOFoptions[6] = {0.};
1250 //define beta per degree of freedom (DOF=nbins in measured)
1251 if(betaColl==0) {
1252 betaPerDOFoptions[0] = 9.1e-9;
1253 betaPerDOFoptions[1] = 2.25e-8;
1254 betaPerDOFoptions[2] = 4.55e-8;
1255 betaPerDOFoptions[3] = 9.1e-8;
1256 betaPerDOFoptions[4] = 2.25e-7;
1257 betaPerDOFoptions[5] = 4.55e-7;
1258 }
1259 if(betaColl==1) {
1260 betaPerDOFoptions[0] = 9.1e-7;
1261 betaPerDOFoptions[1] = 2.25e-6;
1262 betaPerDOFoptions[2] = 4.55e-6;
1263 betaPerDOFoptions[3] = 9.1e-6;
1264 betaPerDOFoptions[4] = 2.25e-5;
1265 betaPerDOFoptions[5] = 4.55e-5;
1266 }
1267 if(betaColl==2) {
1268 betaPerDOFoptions[0] = 9.1e-5;
1269 betaPerDOFoptions[1] = 2.25e-4;
1270 betaPerDOFoptions[2] = 4.55e-4;
1271 betaPerDOFoptions[3] = 9.1e-4;
1272 betaPerDOFoptions[4] = 2.25e-3;
1273 betaPerDOFoptions[5] = 4.55e-3;
1274 }
1275 if(betaColl==3) {
1276 betaPerDOFoptions[0] = 9.1e-3;
1277 betaPerDOFoptions[1] = 2.25e-2;
1278 betaPerDOFoptions[2] = 4.55e-2;
1279 betaPerDOFoptions[3] = 9.1e-2;
1280 betaPerDOFoptions[4] = 2.25e-1;
1281 betaPerDOFoptions[5] = 4.55e-1;
1282 }
1283 if(betaColl==4) {
1284 betaPerDOFoptions[0] = 9.1e-1;
1285 betaPerDOFoptions[1] = 2.25;
1286 betaPerDOFoptions[2] = 4.55;
1287 betaPerDOFoptions[3] = 9.1;
1288 betaPerDOFoptions[4] = 22.5;
1289 betaPerDOFoptions[5] = 45.5;
1290 }
1291
1292 return betaPerDOFoptions[betaOpt];
1293
1294}
1295
31d3e4f0 1296//--------------------------------------------------------------------------------------------------------------------------------------------------
1297Double_t AliAnaChargedJetResponseMaker::InterpolateFast(TGraph *gr, Double_t x) {
1298
1299 Double_t ipmax = gr->GetN()-1.;
1300 Double_t x2,y2,xold,yold;
1301
1302 Double_t xmin,ymin,xmax, ymax;
1303 gr->GetPoint(0,xmin,ymin);
1304 gr->GetPoint(gr->GetN()-1,xmax,ymax);
1305 if(x<xmin || x>xmax) return 0;
1306
1307 Double_t ip = ipmax/2.;
1308 Double_t ipold = 0.;
1309 gr->GetPoint((int)(ip),x2,y2);
1310
1311 Int_t i = 0;
1312
1313 if(x2>x) {
1314 while(x2>x) {
1315 if(i==0) ipold=0.;
1316 ip -= (ip)/2.;
1317 gr->GetPoint((int)(ip),x2,y2);
1318 if(x2>x){}
1319 else ipold = ip;
1320 i++;
1321 // cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl;
1322 }
1323 }
1324 else if(x2<x) {
1325 while(x2<x) {
1326 if(i==0) ipold=ipmax;
1327 ip += (ipold-ip)/2.;
1328 gr->GetPoint((int)(ip),x2,y2);
1329 if(x2>x) ipold = ip;
1330 else {}
1331 i++;
1332 // cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl;
1333 }
1334 }
1335
1336 int ip2 = 0;
1337 if(x2>x) {
1338 ip2 = (int)(ip)-1;
1339 gr->GetPoint(ip2,x2,y2);
1340 while(x2>x) {
1341 ip2--;
1342 gr->GetPoint(ip2,x2,y2);
1343 }
1344 gr->GetPoint(ip2+1,xold,yold);
1345 }
1346 else {
1347 ip2 = (int)(ip)+1;
1348 gr->GetPoint(ip2,x2,y2);
1349 while(x2<x) {
1350 ip2++;
1351 gr->GetPoint(ip2,x2,y2);
1352 }
1353 gr->GetPoint(ip2-1,xold,yold);
1354
1355 }
1356 // cout << "For x=" << x << " interpolate between: " << xold << " and " << x2 << endl;
1357 return ((x-xold)*y2 + (x2-x)*yold) / (x2-xold);
1358
1359}
1360
1361//--------------------------------------------------------------------------------------------------------------------------------------------------
1362Double_t AliAnaChargedJetResponseMaker::InterpolateFast(TH1 *h, Double_t x) {
1363
1364 // Double_t ipmax = gr->GetN()-1.;
1365 Double_t ipmax = (double)h->GetNbinsX();
1366 Double_t x2,y2,xold,yold;
1367
1368 Double_t xmin = h->GetXaxis()->GetXmin();
1369 Double_t xmax = h->GetXaxis()->GetXmax();
1370 if(x<xmin || x>xmax) return 0;
1371
1372 Double_t ip = ipmax/2.;
1373 Double_t ipold = 0.;
1374
1375 x2 = h->GetBinCenter((int)ip);
1376 y2 = h->GetBinContent((int)ip);
1377
1378 Int_t i = 0;
1379
1380 if(x2>x) {
1381 while(x2>x) {
1382 if(i==0) ipold=0.;
1383 ip -= (ip)/2.;
1384 x2 = h->GetBinCenter((int)ip);
1385 y2 = h->GetBinContent((int)ip);
1386 if(x2>x) {}
1387 else ipold = ip;
1388 i++;
1389 // cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl;
1390 }
1391 }
1392 else if(x2<x) {
1393 while(x2<x) {
1394 if(i==0) ipold=ipmax;
1395 ip += (ipold-ip)/2.;
1396 x2 = h->GetBinCenter((int)ip);
1397 y2 = h->GetBinContent((int)ip);
1398 if(x2>x) ipold = ip;
1399 else {}
1400 i++;
1401 // cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl;
1402 }
1403 }
1404
1405 int ip2 = 0;
1406 if(x2>x) {
1407 ip2 = (int)(ip)-1;
1408 x2 = h->GetBinCenter(ip2);
1409 y2 = h->GetBinContent(ip2);
1410 while(x2>x) {
1411 ip2--;
1412 x2 = h->GetBinCenter(ip2);
1413 y2 = h->GetBinContent(ip2);
1414 }
1415 xold = h->GetBinCenter(ip2+1);
1416 yold = h->GetBinContent(ip2+1);
1417 }
1418 else {
1419 ip2 = (int)(ip)+1;
1420 x2 = h->GetBinCenter(ip2);
1421 y2 = h->GetBinContent(ip2);
1422 while(x2<x) {
1423 ip2++;
1424 x2 = h->GetBinCenter(ip2);
1425 y2 = h->GetBinContent(ip2);
1426 }
1427 xold = h->GetBinCenter(ip2-1);
1428 yold = h->GetBinContent(ip2-1);
1429
1430 }
1431 // cout << "For x=" << x << " interpolate between: " << xold << " and " << x2 << endl;
1432 return ((x-xold)*y2 + (x2-x)*yold) / (x2-xold);
1433
1434}
1435
1436