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