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