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