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