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