Compatibility with the Root trunk
[u/mrichter/AliRoot.git] / PWGJE / AliAnaChargedJetResponseMaker.cxx
CommitLineData
31d3e4f0 1#include "AliAnaChargedJetResponseMaker.h"
2#include "TGraph.h"
3#include "TGraphErrors.h"
4#include "TMath.h"
5#include "Riostream.h"
6#include "TH1.h"
7#include "TRandom.h"
8#include "TFile.h"
9#include "TCanvas.h"
10#include "TF1.h"
11#include "THnSparse.h"
12
13#define round(x) ((x)>=0?(long)((x)+0.5):(long)((x)-0.5))
14
c64cb1f6 15using std::cout;
16using std::endl;
17
31d3e4f0 18ClassImp(AliAnaChargedJetResponseMaker)
19
20AliAnaChargedJetResponseMaker::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
31d3e4f0 51
fa7c34ba 52//--------------------------------------------------------------------------------------------------------------------------------------------------
53AliAnaChargedJetResponseMaker::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//--------------------------------------------------------------------------------------------------------------------------------------------------
85AliAnaChargedJetResponseMaker& 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;
31d3e4f0 120}
121
122//--------------------------------------------------------------------------------------------------------------------------------------------------
123void 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//--------------------------------------------------------------------------------------------------------------------------------------------------
180void 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//--------------------------------------------------------------------------------------------------------------------------------------------------
209void 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//--------------------------------------------------------------------------------------------------------------------------------------------------
245void 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//--------------------------------------------------------------------------------------------------------------------------------------------------
271void 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;
fa7c34ba 278 nbins[fDimGen] = fNbins;
31d3e4f0 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
fa7c34ba 336 binArrayPt1[0] = fPtMinUnfolded;
337 for(int i=1; i<nbins[fDimGen]; i++) {
31d3e4f0 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 }
fa7c34ba 343 else
344 binArrayPt1[i] = xmin[fDimGen]+(double)i*binWidth[fDimGen];
31d3e4f0 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//--------------------------------------------------------------------------------------------------------------------------------------------------
365void 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//--------------------------------------------------------------------------------------------------------------------------------------------------
399void AliAnaChargedJetResponseMaker::InitializeResponseMatrixFine() {
400 //
401 //Initialize fine response matrix
402 //
403
404 Int_t nbinsFine[fDimensions*2];
fa7c34ba 405 Double_t xminFine[fDimensions*2];
406 Double_t xmaxFine[fDimensions*2];
31d3e4f0 407 Double_t pTarrayFine[fDimensions*2];
408
fa7c34ba 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.;
31d3e4f0 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//--------------------------------------------------------------------------------------------------------------------------------------------------
463void 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//--------------------------------------------------------------------------------------------------------------------------------------------------
497void 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();
fa7c34ba 516 pTarrayFine[fDimGen] = 0.;
517 pTarrayFine[fDimRec] = 0.;
31d3e4f0 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//--------------------------------------------------------------------------------------------------------------------------------------------------
666TH1D* 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];
31d3e4f0 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//--------------------------------------------------------------------------------------------------------------------------------------------------
749TH1D* 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//--------------------------------------------------------------------------------------------------------------------------------------------------
792Double_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//--------------------------------------------------------------------------------------------------------------------------------------------------
857Double_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