1 #include "AliAnaChargedJetResponseMaker.h"
3 #include "TGraphErrors.h"
11 #include "THnSparse.h"
14 #define round(x) ((x)>=0?(long)((x)+0.5):(long)((x)-0.5))
19 ClassImp(AliAnaChargedJetResponseMaker)
21 AliAnaChargedJetResponseMaker::AliAnaChargedJetResponseMaker():
23 fResolutionType(kParam),
38 fResponseMatrixFine(0x0),
41 fPtMaxUnfoldedHigh(-1.),
42 fBinWidthFactorUnfolded(2),
44 fExtraBinsUnfolded(5),
45 fbVariableBinning(kFALSE),
46 fPtMaxUnfVarBinning(0),
53 //--------------------------------------------------------------------------------------------------------------------------------------------------
54 AliAnaChargedJetResponseMaker::AliAnaChargedJetResponseMaker(const AliAnaChargedJetResponseMaker& obj):
56 fResolutionType(obj.fResolutionType),
57 fDeltaPt(obj.fDeltaPt),
58 fhDeltaPt(obj.fhDeltaPt),
59 fDimensions(obj.fDimensions),
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)
85 //--------------------------------------------------------------------------------------------------------------------------------------------------
86 AliAnaChargedJetResponseMaker& AliAnaChargedJetResponseMaker::operator=(const AliAnaChargedJetResponseMaker& other)
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;
123 //--------------------------------------------------------------------------------------------------------------------------------------------------
124 void AliAnaChargedJetResponseMaker::SetMeasuredSpectrum(TH1D *hPtMeasured)
127 // Set measured spectrum in THnSparse format
128 // This defines the minimum and maximum pT on the reconstructed axis of the response matrix
130 if(fDebug) printf(">>AliAnaChargedJetResponseMaker::SetMeasuredSpectrum \n");
132 fNbins = hPtMeasured->GetXaxis()->GetNbins();
133 fPtMin = hPtMeasured->GetXaxis()->GetXmin();
134 fPtMax = hPtMeasured->GetXaxis()->GetXmax();
136 if(fDebug) printf("fNbins: %d fPtMin: %f fPtMax: %f \n",fNbins,fPtMin,fPtMax);
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);
143 fBinArrayPtRec[fNbins] = hPtMeasured->GetXaxis()->GetBinUpEdge(fNbins);
146 Int_t nbins[fDimensions];
147 Double_t xmin[fDimensions];
148 Double_t xmax[fDimensions];
149 for(int dim = 0; dim<fDimensions; dim++) {
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);
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());
168 for(int i = hPtMeasured->FindBin(fPtMin); i<hPtMeasured->FindBin(fPtMax); i++) {
170 pT[0]= hPtMeasured->GetBinCenter(i);
171 fPtMeasured->SetBinContent(bin,(Double_t)hPtMeasured->GetBinContent(i));
172 fPtMeasured->SetBinError(bin,(Double_t)hPtMeasured->GetBinError(i));
176 if(fDebug) printf("fPtMeasured->GetNbins(): %d \n",(int)fPtMeasured->GetNbins());
180 //--------------------------------------------------------------------------------------------------------------------------------------------------
181 void AliAnaChargedJetResponseMaker::SetFlatEfficiency(Double_t eff) {
184 // Set flat efficiency to value of eff
191 Int_t nbins[fDimensions];
192 Double_t xmin[fDimensions];
193 Double_t xmax[fDimensions];
194 for(int dim = 0; dim<fDimensions; dim++) {
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);
206 for(int i=0; i<fNbins; i++) {
209 fEfficiency->SetBinContent(bin,fEffFlat);
210 fEfficiency->SetBinError(bin,0);
215 //--------------------------------------------------------------------------------------------------------------------------------------------------
216 void AliAnaChargedJetResponseMaker::SetEfficiency(TGraphErrors *grEff)
219 // Fill fEfficiency (THnSparse) with values from grEff
222 Int_t nbins[fDimensions];
223 Double_t xmin[fDimensions];
224 Double_t xmax[fDimensions];
225 for(int dim = 0; dim<fDimensions; dim++) {
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);
241 int nbinsgr = grEff->GetN();
243 for(int i=0; i<nbinsgr; i++) {
244 grEff->GetPoint(i,dummy,yield);
246 error = grEff->GetErrorY(i);
248 fEfficiency->Fill(pT,yield);
249 fEfficiency->SetBinError(i,error);
255 //--------------------------------------------------------------------------------------------------------------------------------------------------
256 void AliAnaChargedJetResponseMaker::MakeResponseMatrixJetsFineMerged(Int_t skipBins, Int_t binWidthFactor, Int_t extraBins, Bool_t bVariableBinning, Double_t ptmax)
259 // Make jet response matrix
262 if(fDebug) printf(">>AliAnaChargedJetResponseMaker::MakeResponseMatrixJetsFineMerged\n");
265 if(fDebug) printf("fPtMeasured does not exist. Aborting!!\n");
268 if(fResponseMatrix) { delete fResponseMatrix; }
269 if(fResponseMatrixFine) { delete fResponseMatrixFine; }
271 SetSkipBinsUnfolded(skipBins);
272 SetBinWidthFactorUnfolded(binWidthFactor);
273 SetExtraBinsUnfolded(extraBins);
274 SetVariableBinning(bVariableBinning,ptmax);
276 InitializeResponseMatrix();
277 InitializeEfficiency();
279 InitializeResponseMatrixFine();
280 InitializeEfficiencyFine();
282 FillResponseMatrixFineAndMerge();
286 //--------------------------------------------------------------------------------------------------------------------------------------------------
287 void AliAnaChargedJetResponseMaker::InitializeResponseMatrix() {
289 //Define bin edges of RM to be used for unfolding
292 Int_t nbins[fDimensions*2];
293 nbins[fDimRec] = fNbins;
294 nbins[fDimGen] = fNbins;
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;
302 fPtMaxUnfolded = fPtMax+(double)(fExtraBinsUnfolded)*binWidthUnf;
303 nbins[fDimGen]+=fExtraBinsUnfolded;
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);
310 printf("fPtMinUnfolded: %f fPtMaxUnfolded: %f\n",fPtMinUnfolded,fPtMaxUnfolded);
313 if(fbVariableBinning) {
314 // cout << "fPtMaxUnfVarBinning: " << fPtMaxUnfVarBinning << " \tfPtMinUnfolded: " << fPtMinUnfolded << " binWidthUnfLowPt: " << binWidthUnfLowPt << endl;
315 Int_t tmp = (int)((fPtMaxUnfVarBinning-fPtMinUnfolded)/binWidthUnfLowPt);
316 tmp = tmp - fSkipBinsUnfolded;
317 fPtMinUnfolded = fPtMaxUnfVarBinning-(double)(tmp)*binWidthUnfLowPt;
318 //cout << "tmp = " << tmp << " fSkipBinsUnfolded = " << fSkipBinsUnfolded << " fPtMinUnfolded = " << fPtMinUnfolded << endl;
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);
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
330 printf(" nbins[fDimGen] = %d nbins[fDimRec] = %d\n",nbins[fDimGen],nbins[fDimRec]);
332 Double_t binWidth[2];
333 binWidth[fDimRec] = (double)((fPtMax-fPtMin)/nbins[fDimRec]);
334 binWidth[fDimGen] = (double)((fPtMaxUnfolded-fPtMinUnfolded)/nbins[fDimGen]);
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;
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]);
348 Double_t binArrayPt0[nbins[fDimRec]+1];
349 Double_t binArrayPt1[nbins[fDimGen]+1];
350 Double_t xmaxGen = TMath::Max(xmax[fDimGen],fPtMaxUnfoldedHigh);
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];
356 binArrayPt0[nbins[fDimRec]]= xmax[fDimRec];
358 //Define bin limits generated/unfolded axis
359 binArrayPt1[0] = fPtMinUnfolded;
360 for(int i=1; i<nbins[fDimGen]; i++) {
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;
367 binArrayPt1[i] = xmin[fDimGen]+(double)i*binWidth[fDimGen];
368 //printf("RM. i = %d \t binArrayPt[i] = %f \n",i,binArrayPt1[i]);
370 binArrayPt1[nbins[fDimGen]]= xmaxGen;
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
376 if(fResponseMatrix) delete fResponseMatrix;
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}");
382 fResponseMatrix->SetBinEdges(fDimRec,binArrayPt0);
383 fResponseMatrix->SetBinEdges(fDimGen,binArrayPt1);
385 Int_t bin[2] = {1,1};
386 for(int i=1; i<fResponseMatrix->GetAxis(0)->GetNbins(); i++) {
388 for(int j=1; j<fResponseMatrix->GetAxis(1)->GetNbins(); j++) {
390 fResponseMatrix->SetBinContent(bin,0.);
398 //--------------------------------------------------------------------------------------------------------------------------------------------------
399 void AliAnaChargedJetResponseMaker::InitializeEfficiency() {
401 // Define dimensions of efficiency THnSparse
404 if(!fResponseMatrix) {
405 printf("AliAnaChargedJetResponseMaker::InitializeEfficiency()\n");
406 printf("Can not define dimensions efficiency without fResponseMatrix dimensions. Aborting \n");
410 TAxis *genAxis = fResponseMatrix->GetAxis(fDimGen);
412 Int_t nbinsEff[fDimensions];
413 Double_t xminEff[fDimensions];
414 Double_t xmaxEff[fDimensions];
416 for(int dim = 0; dim<fDimensions; dim++) {
417 nbinsEff[dim] = genAxis->GetNbins();
418 xminEff[dim] = genAxis->GetXmin();
419 xmaxEff[dim] = genAxis->GetXmax();
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}");
427 const Double_t *binArrayPt = genAxis->GetXbins()->GetArray();
428 fEfficiency->SetBinEdges(0,binArrayPt);
432 //--------------------------------------------------------------------------------------------------------------------------------------------------
433 void AliAnaChargedJetResponseMaker::InitializeResponseMatrixFine() {
435 //Initialize fine response matrix
438 Int_t nbinsFine[fDimensions*2];
439 Double_t xminFine[fDimensions*2];
440 Double_t xmaxFine[fDimensions*2];
441 Double_t pTarrayFine[fDimensions*2];
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.;
452 Double_t binWidth[2];
453 binWidth[fDimRec] = fResponseMatrix->GetAxis(fDimRec)->GetBinWidth(1);
455 Double_t binWidthFine[2];
456 binWidthFine[fDimRec] = binWidth[fDimRec]/((double)fFineFrac);
457 binWidthFine[fDimGen] = binWidthFine[fDimRec]*(double)fBinWidthFactorUnfolded;
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
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]);
468 Double_t binArrayPt0Fine[nbinsFine[fDimRec]+1];
469 Double_t binArrayPt1Fine[nbinsFine[fDimGen]+1];
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]);
475 binArrayPt0Fine[nbinsFine[fDimRec]]= xmaxFine[fDimRec];
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]);
481 binArrayPt1Fine[nbinsFine[fDimGen]]= xmaxFine[fDimGen];
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
486 if(fResponseMatrixFine) delete fResponseMatrixFine;
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}");
492 fResponseMatrixFine->SetBinEdges(fDimRec,binArrayPt0Fine);
493 fResponseMatrixFine->SetBinEdges(fDimGen,binArrayPt1Fine);
495 Int_t bin[2] = {1,1};
496 for(int i=1; i<fResponseMatrixFine->GetAxis(0)->GetNbins(); i++) {
498 for(int j=1; j<fResponseMatrixFine->GetAxis(1)->GetNbins(); j++) {
500 fResponseMatrixFine->SetBinContent(bin,0.);
506 //--------------------------------------------------------------------------------------------------------------------------------------------------
507 void AliAnaChargedJetResponseMaker::InitializeEfficiencyFine() {
509 // Define dimensions of efficiency THnSparse
512 if(!fResponseMatrixFine) {
513 printf("AliAnaChargedJetResponseMaker::InitializeEfficiencyFine()\n");
514 printf("Can not define dimensions efficiency without fResponseMatrixFine dimensions. Aborting \n");
518 TAxis *genAxis = fResponseMatrixFine->GetAxis(fDimGen);
520 Int_t nbinsEff[fDimensions];
521 Double_t xminEff[fDimensions];
522 Double_t xmaxEff[fDimensions];
524 for(int dim = 0; dim<fDimensions; dim++) {
525 nbinsEff[dim] = genAxis->GetNbins();
526 xminEff[dim] = genAxis->GetXmin();
527 xmaxEff[dim] = genAxis->GetXmax();
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}");
535 const Double_t *binArrayPt = genAxis->GetXbins()->GetArray();
536 fEfficiencyFine->SetBinEdges(0,binArrayPt);
540 //--------------------------------------------------------------------------------------------------------------------------------------------------
541 void AliAnaChargedJetResponseMaker::FillResponseMatrixFineAndMerge() {
543 // Fill fine response matrix
546 if(!fResponseMatrix) {
547 printf("Dimensions of fResponseMatrix have to be defined first. Aborting!");
551 if(!fResponseMatrixFine) {
552 printf("Dimensions of fResponseMatrixFine have to be defined first. Aborting!");
556 TAxis *genAxis = fResponseMatrixFine->GetAxis(fDimGen);
557 TAxis *recAxis = fResponseMatrixFine->GetAxis(fDimRec);
559 Int_t nbinsFine[fDimensions*2];
560 Double_t xminFine[fDimensions*2];
561 Double_t xmaxFine[fDimensions*2];
562 Double_t pTarrayFine[fDimensions*2];
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();
570 pTarrayFine[fDimGen] = 0.;
571 pTarrayFine[fDimRec] = 0.;
573 double sumyield = 0.;
574 double sumyielderror2 = 0.;
577 int iptMerged[2] = {0,0};
578 int ierror[2] = {0,0};
581 Double_t pTgen, pTrec;
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.;
594 for(int iy=1; iy<=nbinsFine[fDimGen]; iy++) { //gen
595 pTgen = fResponseMatrixFine->GetAxis(fDimGen)->GetBinCenter(iy);
596 pTarrayFine[fDimGen] = pTgen;
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);
608 else if(fResolutionType==kResiduals) {
609 yield = fhDeltaPt->Interpolate(pTrec-pTgen);
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);
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];
625 sumyielderror2 += error*error;
627 pTarrayFine[fDimRec] = pTrec;
628 ierror[fDimRec] = ix;
629 fResponseMatrixFine->Fill(pTarrayFine,yield);
630 if(fbCalcErrors) fResponseMatrixFine->SetBinError(ierror,error);
634 //Normalize to total number of counts =1
637 for(int ix=1; ix<=nbinsFine[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));
652 fEfficiencyFine->SetBinContent(bin,sumyield);
653 if(fbCalcErrors) fEfficiencyFine->SetBinError(bin,TMath::Sqrt(sumyielderror2));
655 //fill merged response matrix
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());
661 if(fResponseMatrixFine->GetAxis(fDimGen)->GetBinLowEdge(iy) >= fResponseMatrix->GetAxis(fDimGen)->GetXmin()) {
663 iptMerged[fDimGen]=fResponseMatrix->GetAxis(fDimGen)->FindBin(pTgen);
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);
672 weight = f1MergeFunction->Integral(fResponseMatrixFine->GetAxis(fDimGen)->GetBinLowEdge(iy),fResponseMatrixFine->GetAxis(fDimGen)->GetBinUpEdge(iy))/powInteg;
673 // printf("weight: %f \n", weight );
675 weight = 1./((double)fFineFrac);
676 if(fbVariableBinning && pTgen<fPtMaxUnfVarBinning) weight=1./(0.5*(double)fFineFrac);
679 for(int ix=ixMin; ix<=ixMax; ix++) { //loop reconstructed axis
680 pTrec = fResponseMatrixFine->GetAxis(fDimRec)->GetBinCenter(ix);
682 iptMerged[fDimRec]=fResponseMatrix->GetAxis(fDimRec)->FindBin(pTrec);
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;
688 }//loop gen axis fine response matrix
690 } // iy (dimGen) loop
692 //Fill Efficiency corresponding to merged response matrix
693 for(int iy=1; iy<=fResponseMatrix->GetAxis(fDimGen)->GetNbins(); iy++) { //gen
697 for(int ix=1; ix<=fResponseMatrix->GetAxis(fDimRec)->GetNbins(); ix++) { //rec
699 sumyield += fResponseMatrix->GetBinContent(ipt);
701 if(fbCalcErrors) fResponseMatrix->SetBinError(ipt,TMath::Sqrt(errorArray[iy][ix]));
703 printf("RM: pTgen: %f \t sumyield: %f \n",fResponseMatrix->GetAxis(fDimGen)->GetBinCenter(iy),sumyield);
706 fEfficiency->SetBinContent(bin,sumyield);
707 if(fbCalcErrors) fEfficiency->SetBinError(bin,0);
710 if(fDebug) printf("fResponseMatrixFine->GetNbins(): %d \n",(int)fResponseMatrixFine->GetNbins());
711 if(fDebug) printf("fResponseMatrix->GetNbins(): %d \n",(int)fResponseMatrix->GetNbins());
713 if(fDebug) printf("Done constructing response matrix\n");
717 //--------------------------------------------------------------------------------------------------------------------------------------------------
718 TH2* AliAnaChargedJetResponseMaker::MakeResponseMatrixRebin(TH2 *hRMFine, TH2 *hRM, Bool_t useFunctionWeight) {
721 // Rebin matrix hRMFine to dimensions of hRM
722 // function returns matrix in TH2D format (hRM2) with dimensions from hRM
725 TH2 *hRM2 = (TH2*)hRM->Clone("hRM2");
728 if(useFunctionWeight) cout << "Use function to do weighting" << endl;
730 //First normalize columns of input
731 const Int_t nbinsNorm = hRM2->GetNbinsX();
732 cout << "nbinsNorm: " << nbinsNorm << endl;
734 TArrayF *normVector = new TArrayF(nbinsNorm);
736 for(int ix=1; ix<=hRM2->GetNbinsX(); ix++) {
737 Double_t xLow = hRM2->GetXaxis()->GetBinLowEdge(ix);
738 Double_t xUp = hRM2->GetXaxis()->GetBinUpEdge(ix);
739 //cout << "xLow: " << xLow << " xUp: " << xUp << "\t center: " << hRM2->GetXaxis()->GetBinCenter(ix) << endl;
740 Int_t binxLowFine = hRMFine->GetXaxis()->FindBin(xLow);
741 Int_t binxUpFine = hRMFine->GetXaxis()->FindBin(xUp)-1;
742 //cout << "xLowFine: " << hRMFine->GetXaxis()->GetBinLowEdge(binxLowFine) << "\txUpFine: " << hRMFine->GetXaxis()->GetBinUpEdge(binxUpFine) << endl;
743 if(useFunctionWeight)
744 normVector->SetAt(f1MergeFunction->Integral(xLow,xUp),ix-1);
746 normVector->SetAt(hRMFine->Integral(binxLowFine,binxUpFine,1,hRMFine->GetYaxis()->GetNbins()),ix-1);
747 if(fDebug) cout << "ix norm: " << normVector->At(ix-1) << endl;
750 Double_t content, oldcontent = 0.;
753 Double_t xvalLo, xvalUp, yvalLo, yvalUp;
754 Double_t xmin = hRM2->GetXaxis()->GetXmin();
755 Double_t ymin = hRM2->GetYaxis()->GetXmin();
756 Double_t xmax = hRM2->GetXaxis()->GetXmax();
757 Double_t ymax = hRM2->GetYaxis()->GetXmax();
758 for(int ix=1; ix<=hRMFine->GetXaxis()->GetNbins(); ix++) {
759 xvalLo = hRMFine->GetXaxis()->GetBinLowEdge(ix);
760 xvalUp = hRMFine->GetXaxis()->GetBinUpEdge(ix);
761 if(xvalLo<xmin || xvalUp>xmax) continue;
762 ixNew = hRM2->GetXaxis()->FindBin(hRMFine->GetXaxis()->GetBinCenter(ix));
763 for(int iy=1; iy<=hRMFine->GetYaxis()->GetNbins(); iy++) {
764 yvalLo = hRMFine->GetYaxis()->GetBinLowEdge(iy);
765 yvalUp = hRMFine->GetYaxis()->GetBinUpEdge(iy);
766 if(yvalLo<ymin || yvalUp>ymax) continue;
767 content = hRMFine->GetBinContent(ix,iy);
768 iyNew = hRM2->GetYaxis()->FindBin(hRMFine->GetYaxis()->GetBinCenter(iy));
769 oldcontent = hRM2->GetBinContent(ixNew,iyNew);
771 //if(fDebug) cout << "ixNew: " << ixNew << " " << xvalLo << " iyNew: " << iyNew << " " << yvalLo << " content: " << content << " oldContent: " << oldcontent << " newContent: " << oldcontent+content << endl;
772 Double_t weight = 1.;
773 if(normVector->At(ixNew-1)>0.) {
774 if(useFunctionWeight)
775 weight = f1MergeFunction->Integral(xvalLo,xvalUp)/normVector->At(ixNew-1);
777 weight = 1./normVector->At(ixNew-1);
779 hRM2->SetBinContent(ixNew,iyNew,oldcontent+content*weight);
783 if(normVector) delete normVector;
789 //--------------------------------------------------------------------------------------------------------------------------------------------------
790 TH2* AliAnaChargedJetResponseMaker::CreateTruncated2DHisto(TH2 *h2, Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax) {
793 // Limit axis range of 2D histogram
796 Int_t binMinXh2 = h2->GetXaxis()->FindBin(xmin);
797 if(h2->GetXaxis()->GetBinLowEdge(binMinXh2) < xmin ) binMinXh2++;
798 if(h2->GetXaxis()->GetBinLowEdge(binMinXh2) > xmin ) binMinXh2--;
800 Int_t binMinYh2 = h2->GetYaxis()->FindBin(ymin);
801 if(h2->GetYaxis()->GetBinLowEdge(binMinYh2) < ymin ) binMinYh2++;
802 if(h2->GetYaxis()->GetBinLowEdge(binMinYh2) > ymin ) binMinYh2--;
804 Int_t binMaxXh2 = h2->GetXaxis()->FindBin(xmax);
805 if(h2->GetXaxis()->GetBinUpEdge(binMaxXh2) < xmax ) binMaxXh2++;
806 if(h2->GetXaxis()->GetBinUpEdge(binMaxXh2) > xmax ) binMaxXh2--;
808 Int_t binMaxYh2 = h2->GetYaxis()->FindBin(ymax);
809 if(h2->GetYaxis()->GetBinUpEdge(binMaxYh2) < ymax ) binMaxYh2++;
810 if(h2->GetYaxis()->GetBinUpEdge(binMaxYh2) > ymax ) binMaxYh2--;
812 Int_t nbinsX = binMaxXh2-binMinXh2;
813 Int_t nbinsY = binMaxYh2-binMinYh2;
815 Double_t *binsX = new Double_t[nbinsX+1];
816 Double_t *binsY = new Double_t[nbinsY+1];
818 for(int ix=1; ix<=nbinsX; ix++)
819 binsX[ix-1] = h2->GetXaxis()->GetBinLowEdge(binMinXh2+ix-1);
820 binsX[nbinsX] = h2->GetXaxis()->GetBinUpEdge(binMaxXh2);
822 for(int iy=1; iy<=nbinsY; iy++)
823 binsY[iy-1] = h2->GetYaxis()->GetBinLowEdge(binMinYh2+iy-1);
824 binsY[nbinsY] = h2->GetYaxis()->GetBinUpEdge(binMaxYh2);
826 TH2 *h2Lim = new TH2D("h2Lim","h2Lim",nbinsX,binsX,nbinsY,binsY);
828 for(int ix=1; ix<=nbinsX; ix++) {
829 // cout << "ix: " << ix << " " << binsX[ix] << endl;
830 for(int iy=1; iy<=nbinsY; iy++) {
831 cout << "ix: " << ix << " " << binsX[ix] << "\tiy: " << iy << " " << binsY[iy] << endl;
833 double content = h2->GetBinContent(binMinXh2+ix-1,binMinYh2+iy-1);
834 double error = h2->GetBinContent(binMinXh2+ix-1,binMinYh2+iy-1);
835 h2Lim->SetBinContent(ix,iy,content);
836 h2Lim->SetBinError(ix,iy,error);
846 //--------------------------------------------------------------------------------------------------------------------------------------------------
847 TH2* AliAnaChargedJetResponseMaker::TruncateAxisRangeResponseMatrix(TH2 *hRMOrig, Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax) {
850 // Limit axis range of response matrix without changing normalization
854 //TH2 *hRMLimit2 = (TH2*)hRMLimit->Clone("hRMLimit2");
856 TH2 *hRMLimit2 = CreateTruncated2DHisto(hRMOrig, xmin, xmax, ymin, ymax);
859 double binCent[2] = {0.,0.};
862 Int_t binOrig[2] = {0};
863 for(int ix=1; ix<=hRMLimit2->GetXaxis()->GetNbins(); ix++) {
864 binCent[0] = hRMLimit2->GetXaxis()->GetBinCenter(ix);
865 binOrig[0] = hRMOrig->GetXaxis()->FindBin(binCent[0]);
866 for(int iy=1; iy<=hRMLimit2->GetYaxis()->GetNbins(); iy++) {
867 binCent[1] = hRMLimit2->GetYaxis()->GetBinCenter(iy);
868 binOrig[1] = hRMOrig->GetYaxis()->FindBin(binCent[1]);
870 content = hRMOrig->GetBinContent(binOrig[0],binOrig[1]);
871 error = hRMOrig->GetBinError(binOrig[0],binOrig[1]);
873 hRMLimit2->SetBinContent(ix,iy,content);
874 hRMLimit2->SetBinError(ix,iy,error);
884 //--------------------------------------------------------------------------------------------------------------------------------------------------
885 TH2* AliAnaChargedJetResponseMaker::MultiplityResponseMatrices(TH2 *h2RMDeltaPt, TH2 *h2RMDetector) {
888 // Make combined response matrix (background flucutuations + detector effects)
890 // hRMDeltaPt is the response matrix for background fluctuations from \delta(p_t) measurement
891 // 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
893 // Function assumes that generated/unfolded axis is x-axis and reconstructed is on y-axis on both respone matrices
896 TH2D *h2ResponseMatrixCombined = (TH2D*)h2RMDeltaPt->Clone("h2ResponseMatrixCombined"); //h2RMDeltaPt is the bkg fluctuations RM which has the dimensions we want for the combined response matrix
897 h2ResponseMatrixCombined->SetTitle("h2ResponseMatrixCombined");
898 h2ResponseMatrixCombined->SetName("h2ResponseMatrixCombined");
900 // M = RM_deltaPt * RM_DetEffects * T (M=measured T=truth)
902 // mx1 = mxd * dxt * tx1
905 // d = rec for RM_DetEffects and gen for RM_deltaPt
906 // RM_comb = RM_deltaPt * RM_DetEffects (dimensions mxt)
907 // RM_comb(m,t) = Sum_d RM_deltaPt(m,d)*RM_DetEffects(d,t)
910 printf("Nt=%d\n",h2ResponseMatrixCombined->GetNbinsX());
911 printf("Nm=%d\n",h2ResponseMatrixCombined->GetNbinsY());
912 printf("Nd=%d\n",h2RMDeltaPt->GetNbinsX());
916 for(Int_t t=1; t<=h2ResponseMatrixCombined->GetNbinsX();t++) {
917 for(Int_t m=1; m<=h2ResponseMatrixCombined->GetNbinsY();m++) {
918 Double_t valueSum = 0.;
919 for(Int_t d=1; d<=h2RMDeltaPt->GetNbinsX();d++) {
920 valueSum += h2RMDeltaPt->GetBinContent(d,m) * h2RMDetector->GetBinContent(t,d);
921 // if(t==10 && m==10) cout << "sum m,d=" << m << "," << d << endl;
923 // if(t==10) cout << "t,m = " << t << "," << m << "\tvalueSum: " << valueSum << endl;
924 h2ResponseMatrixCombined->SetBinContent(t,m,valueSum);
928 return h2ResponseMatrixCombined;
932 //--------------------------------------------------------------------------------------------------------------------------------------------------
933 TH2* AliAnaChargedJetResponseMaker::GetTransposeResponsMatrix(TH2 *h2RM) {
935 // Transpose response matrix
938 //Initialize transposed response matrix h2RMTranspose
939 TArrayD *arrayX = (TArrayD*)h2RM->GetXaxis()->GetXbins();
940 for(int i=0; i<arrayX->GetSize(); i++) cout << "i: " << arrayX->At(i) << endl;
941 Double_t *xbinsArray = arrayX->GetArray();
943 TArrayD *arrayY = (TArrayD*)h2RM->GetYaxis()->GetXbins();
944 for(int i=0; i<arrayY->GetSize(); i++) cout << "i: " << arrayY->At(i) << endl;
945 Double_t *ybinsArray = arrayY->GetArray();
947 TH2D *h2RMTranspose = new TH2D("h2RMTranspose","h2RMTranspose",h2RM->GetNbinsY(),ybinsArray,h2RM->GetNbinsX(),xbinsArray);
949 //Fill tranposed response matrix
950 for (Int_t ibin = 1; ibin <= h2RMTranspose->GetNbinsX(); ibin++) {
951 for (Int_t jbin = 1; jbin <= h2RMTranspose->GetNbinsY(); jbin++) {
952 h2RMTranspose->SetBinContent(ibin,jbin, h2RM->GetBinContent(jbin,ibin));
957 return h2RMTranspose;
961 //--------------------------------------------------------------------------------------------------------------------------------------------------
962 TH2* AliAnaChargedJetResponseMaker::NormalizeResponsMatrixYaxisWithPrior(TH2 *h2RM, TH1 *hPrior) {
964 // Normalize such that the Y projection is the prior
967 // TH1D *hProjRespY = (TH1D*)h2RM->ProjectionY("hProjRespY");
968 double intPrior = hPrior->Integral();//"width");
969 for (Int_t jbin = 1; jbin <= h2RM->GetNbinsY(); jbin++) {
970 // double corr = hPrior->GetBinContent(jbin)/hProjRespY->GetBinContent(jbin);
971 for (Int_t ibin = 1; ibin <= h2RM->GetNbinsX(); ibin++) {
972 double content = h2RM->GetBinContent(ibin,jbin);
973 // h2RM->SetBinContent(ibin,jbin,content*corr);
974 h2RM->SetBinContent(ibin,jbin,hPrior->GetBinContent(jbin)/hPrior->GetBinWidth(jbin)/intPrior*content);
981 //--------------------------------------------------------------------------------------------------------------------------------------------------
982 TH1D* AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TH1 *hGen, TH2 *hResponse,TH1 *hEfficiency,Bool_t bDrawSlices) {
985 // Multiply hGen with response matrix to obtain refolded spectrum
986 // Efficiency must be given. In case efficiency is 1 use SetFlatEfficiency(1.) before calling function
990 // printf("Efficiency must be given. In case efficiency is 1 use SetFlatEfficiency(1.) before calling function. Aborting!");
992 printf("Setting efficiency to 1 \n");
993 hEfficiency = (TH1D*)hGen->Clone("hEfficiency");
994 hEfficiency->Reset();
995 for(int i=1; i<=hEfficiency->GetNbinsX(); i++) hEfficiency->SetBinContent(i,1.);
1000 //y-axis: reconstructed
1001 if(fDebug>0) cout << "nbins hGen: " << hGen->GetNbinsX() << "\t nbins hResponseGen: " << hResponse->GetXaxis()->GetNbins() << "\t nbins hResponseRec: " << hResponse->GetYaxis()->GetNbins() << endl;
1003 TH1D *hRec = hResponse->ProjectionY("hRec");
1006 hRec->SetTitle("hRec");
1007 hRec->SetName("hRec");
1009 for(int irec=1; irec<=hRec->GetNbinsX(); irec++)
1010 hRec->SetBinContent(irec,0);
1012 TH1D *hRecGenBin = 0x0;
1013 TCanvas *cSlices = 0x0;
1015 cSlices = new TCanvas("cSlices","cSlices:Slices",400,400);
1019 Double_t yieldMC = 0.;
1020 Double_t yieldMCerror = 0.;
1021 Double_t sumYield = 0.;
1022 const Int_t nbinsRec = hRec->GetNbinsX();
1023 Double_t sumError2[nbinsRec+1];
1026 for(int igen=1; igen<=hGen->GetNbinsX(); igen++) {
1032 eff = hEfficiency->GetBinContent(igen);
1033 yieldMC = hGen->GetBinContent(igen)*eff;
1034 yieldMCerror = hGen->GetBinError(igen)*eff;
1037 hRecGenBin = hResponse->ProjectionY(Form("hRecGenBin%d",igen));
1038 hRecGenBin->Sumw2();
1039 hRecGenBin->Reset();
1040 hRecGenBin->SetTitle(Form("hRecGenBin%d",igen));
1041 hRecGenBin->SetName(Form("hRecGenBin%d",igen));
1044 for(int irec=1; irec<=hRec->GetNbinsX(); irec++) {
1045 hRec->AddBinContent(irec,yieldMC*hResponse->GetBinContent(igen,irec));
1046 sumYield+=hResponse->GetBinContent(igen,irec);
1047 // Double_t A = yieldMC*hResponse->GetBinError(igen,irec);
1048 Double_t B = hResponse->GetBinContent(igen,irec)*yieldMCerror;
1049 // Double_t tmp2 = A*A + B*B;
1050 //sumError2[irec-1] += tmp2 ;
1051 sumError2[irec-1] += B*B;
1054 hRecGenBin->SetBinContent(irec,yieldMC*hResponse->GetBinContent(igen,irec));
1059 hRecGenBin->SetLineColor(igen);
1060 if(igen==1) hRecGenBin->DrawCopy();
1061 else hRecGenBin->DrawCopy("same");
1064 if(hRecGenBin) delete hRecGenBin;
1066 cout << "igen: " << igen << "\tpTMC: " << hGen->GetXaxis()->GetBinCenter(igen) << "\teff:" << eff << "\tsumYield: " << sumYield << endl;
1069 for(int i=0; i<=nbinsRec; i++) {
1071 hRec->SetBinError(i+1,TMath::Sqrt(sumError2[i]));
1078 //--------------------------------------------------------------------------------------------------------------------------------------------------
1079 TH1D* AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TF1 *fGen, TH2 *hResponse,TH1 *hEfficiency) {
1082 // Multiply fGen function with response matrix to obtain (re)folded spectrum
1083 // Efficiency must be given. In case efficiency is 1 use SetFlatEfficiency(1.) before calling function
1088 //y-axis: reconstructed
1090 if(fDebug>0) printf(">>AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TF1 *fGen, TH2 *hResponse,TH1 *hEfficiency)\n");
1092 TH1D *hRec = hResponse->ProjectionY("hRec");
1095 hRec->SetTitle("hRec");
1096 hRec->SetName("hRec");
1098 // TH1D *hRec = new TH1D("hRec","hRec",hResponse->GetNbinsY(),hResponse->GetYaxis()->GetXmin(),hResponse->GetYaxis()->GetXmax());
1100 for(int irec=1; irec<=hRec->GetNbinsX(); irec++)
1101 hRec->SetBinContent(irec,0);
1103 Double_t yieldMC = 0.;
1104 Double_t sumYield = 0.;
1106 for(int igen=1; igen<=hResponse->GetNbinsX(); igen++) {
1109 double pTMC = hResponse->GetXaxis()->GetBinCenter(igen);
1111 int binEff = hEfficiency->FindBin(pTMC);
1112 eff = hEfficiency->GetBinContent(binEff);
1115 // yieldMC = fGen->Eval(pTMC)*eff;
1116 yieldMC = fGen->Integral(hResponse->GetXaxis()->GetBinLowEdge(igen),hResponse->GetXaxis()->GetBinUpEdge(igen))*eff;
1117 for(int irec=1; irec<=hResponse->GetNbinsY(); irec++) {
1118 hRec->AddBinContent(irec,yieldMC*hResponse->GetBinContent(igen,irec));
1119 sumYield+=hResponse->GetBinContent(igen,irec);
1121 // cout << "igen: " << igen << "\tpTMC: " << pTMC << "\tsumYield: " << sumYield << endl;
1122 // cout << "yieldMC: " << yieldMC << endl;
1129 //--------------------------------------------------------------------------------------------------------------------------------------------------
1130 Double_t AliAnaChargedJetResponseMaker::InterpolateFast(TGraph *gr, Double_t x) {
1132 Double_t ipmax = gr->GetN()-1.;
1133 Double_t x2,y2,xold,yold;
1135 Double_t xmin,ymin,xmax, ymax;
1136 gr->GetPoint(0,xmin,ymin);
1137 gr->GetPoint(gr->GetN()-1,xmax,ymax);
1138 if(x<xmin || x>xmax) return 0;
1140 Double_t ip = ipmax/2.;
1141 Double_t ipold = 0.;
1142 gr->GetPoint((int)(ip),x2,y2);
1150 gr->GetPoint((int)(ip),x2,y2);
1154 // cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl;
1159 if(i==0) ipold=ipmax;
1160 ip += (ipold-ip)/2.;
1161 gr->GetPoint((int)(ip),x2,y2);
1162 if(x2>x) ipold = ip;
1165 // cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl;
1172 gr->GetPoint(ip2,x2,y2);
1175 gr->GetPoint(ip2,x2,y2);
1177 gr->GetPoint(ip2+1,xold,yold);
1181 gr->GetPoint(ip2,x2,y2);
1184 gr->GetPoint(ip2,x2,y2);
1186 gr->GetPoint(ip2-1,xold,yold);
1189 // cout << "For x=" << x << " interpolate between: " << xold << " and " << x2 << endl;
1190 return ((x-xold)*y2 + (x2-x)*yold) / (x2-xold);
1194 //--------------------------------------------------------------------------------------------------------------------------------------------------
1195 Double_t AliAnaChargedJetResponseMaker::InterpolateFast(TH1 *h, Double_t x) {
1197 // Double_t ipmax = gr->GetN()-1.;
1198 Double_t ipmax = (double)h->GetNbinsX();
1199 Double_t x2,y2,xold,yold;
1201 Double_t xmin = h->GetXaxis()->GetXmin();
1202 Double_t xmax = h->GetXaxis()->GetXmax();
1203 if(x<xmin || x>xmax) return 0;
1205 Double_t ip = ipmax/2.;
1206 Double_t ipold = 0.;
1208 x2 = h->GetBinCenter((int)ip);
1209 y2 = h->GetBinContent((int)ip);
1217 x2 = h->GetBinCenter((int)ip);
1218 y2 = h->GetBinContent((int)ip);
1222 // cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl;
1227 if(i==0) ipold=ipmax;
1228 ip += (ipold-ip)/2.;
1229 x2 = h->GetBinCenter((int)ip);
1230 y2 = h->GetBinContent((int)ip);
1231 if(x2>x) ipold = ip;
1234 // cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl;
1241 x2 = h->GetBinCenter(ip2);
1242 y2 = h->GetBinContent(ip2);
1245 x2 = h->GetBinCenter(ip2);
1246 y2 = h->GetBinContent(ip2);
1248 xold = h->GetBinCenter(ip2+1);
1249 yold = h->GetBinContent(ip2+1);
1253 x2 = h->GetBinCenter(ip2);
1254 y2 = h->GetBinContent(ip2);
1257 x2 = h->GetBinCenter(ip2);
1258 y2 = h->GetBinContent(ip2);
1260 xold = h->GetBinCenter(ip2-1);
1261 yold = h->GetBinContent(ip2-1);
1264 // cout << "For x=" << x << " interpolate between: " << xold << " and " << x2 << endl;
1265 return ((x-xold)*y2 + (x2-x)*yold) / (x2-xold);