- Set Sumw2 for all histos in the histo container
[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 #include "TArrayD.h"
13
14 #define round(x) ((x)>=0?(long)((x)+0.5):(long)((x)-0.5))
15
16 using std::cout;
17 using std::endl;
18
19 ClassImp(AliAnaChargedJetResponseMaker)
20
21 AliAnaChargedJetResponseMaker::AliAnaChargedJetResponseMaker(): 
22   fDebug(kFALSE),
23   fResolutionType(kParam),
24   fDeltaPt(0x0),
25   fhDeltaPt(0x0),
26   fh1MeasuredTruncated(0x0),
27   fh2DetectorResponse(0x0),
28   fh2DetectorResponseRebin(0x0),
29   fhEfficiencyDet(0x0),
30   fh2ResponseMatrixCombinedFineFull(0x0),
31   fh2ResponseMatrixCombinedFull(0x0),
32   fh2ResponseMatrixCombined(0x0),
33   fhEfficiencyCombined(0x0),
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
60
61 //--------------------------------------------------------------------------------------------------------------------------------------------------
62 AliAnaChargedJetResponseMaker::AliAnaChargedJetResponseMaker(const AliAnaChargedJetResponseMaker& obj):
63   fDebug(obj.fDebug),
64   fResolutionType(obj.fResolutionType),
65   fDeltaPt(obj.fDeltaPt),
66   fhDeltaPt(obj.fhDeltaPt),
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),
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 //--------------------------------------------------------------------------------------------------------------------------------------------------
102 AliAnaChargedJetResponseMaker& AliAnaChargedJetResponseMaker::operator=(const AliAnaChargedJetResponseMaker& other)
103 {
104 // Assignment
105   if(&other == this) return *this;
106   AliAnaChargedJetResponseMaker::operator=(other);
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;
143
144   return *this;
145 }
146
147 //--------------------------------------------------------------------------------------------------------------------------------------------------
148 void AliAnaChargedJetResponseMaker::SetMeasuredSpectrum(TH1D *hPtMeasured)
149 {
150   //
151   // Set measured spectrum in THnSparse format
152   // This defines the minimum and maximum pT on the reconstructed axis of the response matrix
153   //
154   if(fDebug) printf(">>AliAnaChargedJetResponseMaker::SetMeasuredSpectrum \n");
155
156   fNbins = hPtMeasured->GetXaxis()->GetNbins();
157   fPtMin = hPtMeasured->GetXaxis()->GetXmin();
158   fPtMax = hPtMeasured->GetXaxis()->GetXmax();
159
160   if(fDebug) printf("fNbins: %d  fPtMin: %f  fPtMax: %f \n",fNbins,fPtMin,fPtMax);
161   
162   if(fBinArrayPtRec) delete fBinArrayPtRec;
163   fBinArrayPtRec = new Double_t[fNbins+1];
164   for(int j = 0; j<fNbins; j++) {
165     fBinArrayPtRec[j] = hPtMeasured->GetXaxis()->GetBinLowEdge(j+1);
166   }
167   fBinArrayPtRec[fNbins] = hPtMeasured->GetXaxis()->GetBinUpEdge(fNbins);
168   
169
170   Int_t nbins[fDimensions];
171   Double_t xmin[fDimensions]; 
172   Double_t xmax[fDimensions];
173   for(int dim = 0; dim<fDimensions; dim++) {
174     nbins[dim] = fNbins;
175     xmin[dim]  = fPtMin;
176     xmax[dim]  = fPtMax;
177   }
178
179   if(fPtMeasured) delete fPtMeasured;
180   fPtMeasured = new THnSparseD("fPtMeasured","Measured pT spectrum; p_{T}^{rec}",fDimensions,nbins,xmin,xmax);
181   fPtMeasured->Sumw2();
182   fPtMeasured->GetAxis(0)->SetTitle("p_{T}^{rec}");
183   fPtMeasured->SetBinEdges(0,fBinArrayPtRec);
184
185   //Fill
186   if(fDebug) printf("fill measured THnSparse\n");
187   if(fNbins!=hPtMeasured->GetNbinsX()) 
188     printf("WARNING: nbins not correct \t %d vs %d !!!\n",fNbins,hPtMeasured->GetNbinsX());
189
190   int bin[1] = {0};
191   bin[0] = 0;
192   for(int i = hPtMeasured->FindBin(fPtMin); i<hPtMeasured->FindBin(fPtMax); i++) {
193     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 //--------------------------------------------------------------------------------------------------------------------------------------------------
203 void AliAnaChargedJetResponseMaker::SetFlatEfficiency(Double_t eff) {
204
205   //
206   // Set flat efficiency to value of eff
207   //
208
209   fEffFlat = eff;
210   return;
211
212   /*
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   }
234   */
235 }
236
237 //--------------------------------------------------------------------------------------------------------------------------------------------------
238 void AliAnaChargedJetResponseMaker::SetEfficiency(TGraphErrors *grEff)
239 {
240   //
241   // Fill fEfficiency (THnSparse) with values from grEff
242   //
243
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 //--------------------------------------------------------------------------------------------------------------------------------------------------
278 void 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
284   if(fDebug) printf(">>AliAnaChargedJetResponseMaker::MakeResponseMatrixJetsFineMerged\n");
285
286   if(!fPtMeasured) {
287     if(fDebug) printf("fPtMeasured does not exist. Aborting!!\n");
288     return;
289   }
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 //--------------------------------------------------------------------------------------------------------------------------------------------------
309 void AliAnaChargedJetResponseMaker::InitializeResponseMatrix() {
310   //
311   //Define bin edges of RM to be used for unfolding
312   //
313
314   Int_t nbins[fDimensions*2];
315   nbins[fDimRec] = fNbins;
316   nbins[fDimGen] = fNbins;
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
324   fPtMaxUnfolded = fPtMax+(double)(fExtraBinsUnfolded)*binWidthUnf;
325   nbins[fDimGen]+=fExtraBinsUnfolded;
326
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
332   printf("fPtMinUnfolded: %f  fPtMaxUnfolded: %f\n",fPtMinUnfolded,fPtMaxUnfolded);
333
334
335   if(fbVariableBinning) {
336     // cout << "fPtMaxUnfVarBinning: " << fPtMaxUnfVarBinning << " \tfPtMinUnfolded: " << fPtMinUnfolded << "  binWidthUnfLowPt: " << binWidthUnfLowPt << endl;
337     Int_t tmp = (int)((fPtMaxUnfVarBinning-fPtMinUnfolded)/binWidthUnfLowPt);
338     tmp = tmp - fSkipBinsUnfolded;
339     fPtMinUnfolded = fPtMaxUnfVarBinning-(double)(tmp)*binWidthUnfLowPt;  
340     //cout << "tmp = " << tmp << "  fSkipBinsUnfolded = " << fSkipBinsUnfolded << " fPtMinUnfolded = " << fPtMinUnfolded << endl;
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   }
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]);
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
381   binArrayPt1[0] = fPtMinUnfolded;
382   for(int i=1; i<nbins[fDimGen]; i++) {
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     }
388     else
389       binArrayPt1[i] = xmin[fDimGen]+(double)i*binWidth[fDimGen];
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
398   if(fResponseMatrix) delete fResponseMatrix;
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
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
417
418 }
419
420 //--------------------------------------------------------------------------------------------------------------------------------------------------
421 void 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 //--------------------------------------------------------------------------------------------------------------------------------------------------
455 void AliAnaChargedJetResponseMaker::InitializeResponseMatrixFine() {
456   //
457   //Initialize fine response matrix
458   //
459
460   Int_t nbinsFine[fDimensions*2];
461   Double_t xminFine[fDimensions*2];
462   Double_t xmaxFine[fDimensions*2];
463   // Double_t pTarrayFine[fDimensions*2];
464
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.);
469   xmaxFine[fDimRec] = fResponseMatrix->GetAxis(fDimGen)->GetXmax();//+40.;
470   xmaxFine[fDimGen] = fResponseMatrix->GetAxis(fDimGen)->GetXmax();//+40.;
471   // pTarrayFine[fDimRec] = 0.;
472   // pTarrayFine[fDimGen] = 0.;
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;
480   //  binWidthFine[fDimRec] = binWidthFine[fDimGen];
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
509   if(fResponseMatrixFine) delete fResponseMatrixFine;
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
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
527 }
528
529 //--------------------------------------------------------------------------------------------------------------------------------------------------
530 void 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 //--------------------------------------------------------------------------------------------------------------------------------------------------
564 void AliAnaChargedJetResponseMaker::FillResponseMatrixFineAndMerge() {
565   //
566   // Fill fine response matrix
567   //
568
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
579   TAxis *genAxis = fResponseMatrixFine->GetAxis(fDimGen);
580   TAxis *recAxis = fResponseMatrixFine->GetAxis(fDimRec);
581
582   Int_t nbinsFine[fDimensions*2];
583   //Double_t xminFine[fDimensions*2]; 
584   //Double_t xmaxFine[fDimensions*2];
585   Double_t pTarrayFine[fDimensions*2];
586
587   nbinsFine[fDimGen] = genAxis->GetNbins();
588   nbinsFine[fDimRec] = recAxis->GetNbins();
589   //  xminFine[fDimGen]  = genAxis->GetXmin();
590   //  xminFine[fDimRec]  = recAxis->GetXmin();
591   //  xmaxFine[fDimGen]  = genAxis->GetXmax();
592   //  xmaxFine[fDimRec]  = recAxis->GetXmax();
593   pTarrayFine[fDimGen] = 0.;
594   pTarrayFine[fDimRec] = 0.;
595
596   double sumyield = 0.;
597   double sumyielderror2 = 0.;
598
599   int ipt[2]    = {0,0};
600   int iptMerged[2]    = {0,0};
601   int ierror[2] = {0,0};
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
680     //find bin in fine RM corresponding to mimimum/maximum bin of merged RM on rec axis
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);
699         if(fbVariableBinning && pTgen<fPtMaxUnfVarBinning) weight=1./(0.5*(double)fFineFrac);
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
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
740 //--------------------------------------------------------------------------------------------------------------------------------------------------
741 TH2* AliAnaChargedJetResponseMaker::MakeResponseMatrixRebin(TH2 *hRMFine, TH2 *hRM, Bool_t useFunctionWeight) {
742
743   //
744   // Rebin matrix hRMFine to dimensions of hRM
745   // function returns matrix in TH2D format (hRM2) with dimensions from hRM
746   //
747
748   TH2 *hRM2 = (TH2*)hRM->Clone("hRM2");
749   hRM2->Reset();
750
751   if(useFunctionWeight)  cout << "Use function to do weighting" << endl;
752
753   //First normalize columns of input
754   const Int_t nbinsNorm = hRM2->GetNbinsX();
755   if(fDebug) cout << "nbinsNorm: " << nbinsNorm << endl;
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;
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);
770     //if(fDebug) cout << "ix norm: " << normVector->At(ix-1) << endl;
771   }
772
773   Double_t content, oldcontent = 0.;
774   Int_t ixNew = 0;
775   Int_t iyNew = 0;
776   Double_t xvalLo, xvalUp, yvalLo, yvalUp;
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++) {
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));
786     for(int iy=1; iy<=hRMFine->GetYaxis()->GetNbins(); iy++) {
787       yvalLo = hRMFine->GetYaxis()->GetBinLowEdge(iy);
788       yvalUp = hRMFine->GetYaxis()->GetBinUpEdge(iy);
789       if(yvalLo<ymin || yvalUp>ymax) continue;
790       content = hRMFine->GetBinContent(ix,iy);
791       iyNew = hRM2->GetYaxis()->FindBin(hRMFine->GetYaxis()->GetBinCenter(iy));
792       oldcontent = hRM2->GetBinContent(ixNew,iyNew);
793
794       //if(fDebug) cout << "ixNew: " << ixNew << " " << xvalLo << " iyNew: " << iyNew << " " << yvalLo << " content: " << content << " oldContent: " << oldcontent << " newContent: " << oldcontent+content << endl;
795       Double_t weight = 1.;
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       }
802       hRM2->SetBinContent(ixNew,iyNew,oldcontent+content*weight);
803     }
804   }
805
806   if(normVector) delete normVector;
807   
808   return hRM2;
809
810 }
811
812 //--------------------------------------------------------------------------------------------------------------------------------------------------
813 TH2* 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
835   Int_t nbinsX = binMaxXh2-binMinXh2+1;
836   Int_t nbinsY = binMaxYh2-binMinYh2+1;
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++) {
854       // cout << "ix: " << ix << "  " << binsX[ix] << "\tiy: " << iy << "  " << binsY[iy] << endl;
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);
859       if(fbCalcErrors) h2Lim->SetBinError(ix,iy,error);
860
861     }
862   }
863
864
865
866   return h2Lim;
867 }
868
869 //--------------------------------------------------------------------------------------------------------------------------------------------------
870 TH2* 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);
897       if(fbCalcErrors) hRMLimit2->SetBinError(ix,iy,error);
898
899     }
900   }
901   
902
903   return hRMLimit2;
904
905 }
906
907 //--------------------------------------------------------------------------------------------------------------------------------------------------
908 Bool_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 //--------------------------------------------------------------------------------------------------------------------------------------------------
922 void 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
997 //--------------------------------------------------------------------------------------------------------------------------------------------------
998 TH2* 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
1009   TH2D *h2ResponseMatrixCombined = (TH2D*)h2RMDeltaPt->Clone("h2ResponseMatrixCombined"); //h2RMDeltaPt is the bkg fluctuations RM which has the dimensions we want for the combined response matrix
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) {
1023     printf("Nt=%d\n",h2ResponseMatrixCombined->GetNbinsX());
1024     printf("Nm=%d\n",h2ResponseMatrixCombined->GetNbinsY());
1025     printf("Nd=%d\n",h2RMDeltaPt->GetNbinsX());
1026   }
1027
1028
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);
1034         //      if(t==10 && m==10) cout << "sum m,d=" << m << "," << d << endl;
1035       }//d-loop
1036       //  if(t==10) cout << "t,m = " << t << "," << m << "\tvalueSum: " << valueSum << endl; 
1037       h2ResponseMatrixCombined->SetBinContent(t,m,valueSum);
1038     } //m-loop
1039   }//t-loop
1040
1041   return h2ResponseMatrixCombined;
1042
1043 }
1044
1045 //--------------------------------------------------------------------------------------------------------------------------------------------------
1046 TH2* AliAnaChargedJetResponseMaker::GetTransposeResponsMatrix(TH2 *h2RM) {
1047   //
1048   // Transpose response matrix
1049   //
1050
1051   Printf("AliAnaChargedJetResponseMaker::GetTransposeResponsMatrix");
1052
1053   //Initialize transposed response matrix h2RMTranspose
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());
1074
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 //--------------------------------------------------------------------------------------------------------------------------------------------------
1091 TH2* AliAnaChargedJetResponseMaker::NormalizeResponsMatrixYaxisWithPrior(TH2 *h2RM, TH1 *hPrior) {
1092   //
1093   // Normalize such that the Y projection is the prior
1094   //
1095
1096   //  TH1D *hProjRespY = (TH1D*)h2RM->ProjectionY("hProjRespY");
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
1110 //--------------------------------------------------------------------------------------------------------------------------------------------------
1111 TH1D* AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TH1 *hGen, TH2 *hResponse,TH1 *hEfficiency,Bool_t bDrawSlices) {
1112
1113   //
1114   // Multiply hGen with response matrix to obtain refolded spectrum
1115   // Efficiency must be given. In case efficiency is 1 use SetFlatEfficiency(1.) before calling function
1116   //
1117
1118   if(!hEfficiency) {
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.);
1125   }
1126
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];
1153   for(int irec=0; irec<=hRec->GetNbinsX(); irec++)
1154     sumError2[irec]=0.;
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;
1165     yieldMCerror = hGen->GetBinError(igen)*eff;
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);
1178       //      Double_t A = yieldMC*hResponse->GetBinError(igen,irec);
1179       Double_t B = hResponse->GetBinContent(igen,irec)*yieldMCerror;
1180       //      Double_t tmp2 = A*A + B*B;
1181       //sumError2[irec-1] += tmp2 ;
1182       sumError2[irec-1] += B*B;
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   
1200   for(int i=0; i<=nbinsRec; i++) {
1201     if(sumError2[i]>0.)
1202       hRec->SetBinError(i+1,TMath::Sqrt(sumError2[i]));
1203   }
1204
1205
1206   return hRec;
1207 }
1208
1209 //--------------------------------------------------------------------------------------------------------------------------------------------------
1210 TH1D* AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TF1 *fGen, TH2 *hResponse,TH1 *hEfficiency) {
1211
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
1217   //For response
1218   //x-axis: generated
1219   //y-axis: reconstructed
1220
1221   if(fDebug>0) printf(">>AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TF1 *fGen, TH2 *hResponse,TH1 *hEfficiency)\n");
1222
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);
1241     if(hEfficiency) { 
1242       int binEff = hEfficiency->FindBin(pTMC);
1243       eff = hEfficiency->GetBinContent(binEff);
1244     }
1245     else eff = 1.;
1246     // yieldMC = fGen->Eval(pTMC)*eff;
1247     yieldMC = fGen->Integral(hResponse->GetXaxis()->GetBinLowEdge(igen),hResponse->GetXaxis()->GetBinUpEdge(igen))*eff;
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     }
1252     //    cout << "igen: " << igen << "\tpTMC: " << pTMC << "\tsumYield: " << sumYield << endl;
1253     //    cout << "yieldMC: " << yieldMC << endl;
1254     
1255   }
1256
1257   return hRec;
1258 }
1259
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
1310 //--------------------------------------------------------------------------------------------------------------------------------------------------
1311 Double_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 //--------------------------------------------------------------------------------------------------------------------------------------------------
1376 Double_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