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