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),
26 fh1MeasuredTruncated(0x0),
27 fh2DetectorResponse(0x0),
28 fh2DetectorResponseRebin(0x0),
30 fh2ResponseMatrixCombinedFineFull(0x0),
31 fh2ResponseMatrixCombinedFull(0x0),
32 fh2ResponseMatrixCombined(0x0),
33 fhEfficiencyCombined(0x0),
46 fResponseMatrixFine(0x0),
49 fPtMaxUnfoldedHigh(-1.),
50 fBinWidthFactorUnfolded(2),
52 fExtraBinsUnfolded(5),
53 fbVariableBinning(kFALSE),
54 fPtMaxUnfVarBinning(0),
61 //--------------------------------------------------------------------------------------------------------------------------------------------------
62 AliAnaChargedJetResponseMaker::AliAnaChargedJetResponseMaker(const AliAnaChargedJetResponseMaker& obj):
64 fResolutionType(obj.fResolutionType),
65 fDeltaPt(obj.fDeltaPt),
66 fhDeltaPt(obj.fhDeltaPt),
67 fh1MeasuredTruncated(obj.fh1MeasuredTruncated),
68 fh2DetectorResponse(obj.fh2DetectorResponse),
69 fh2DetectorResponseRebin(obj.fh2DetectorResponseRebin),
70 fhEfficiencyDet(obj.fhEfficiencyDet),
71 fh2ResponseMatrixCombinedFineFull(obj.fh2ResponseMatrixCombinedFineFull),
72 fh2ResponseMatrixCombinedFull(obj.fh2ResponseMatrixCombinedFull),
73 fh2ResponseMatrixCombined(obj.fh2ResponseMatrixCombined),
74 fhEfficiencyCombined(obj.fhEfficiencyCombined),
75 fDimensions(obj.fDimensions),
81 fBinArrayPtRec(obj.fBinArrayPtRec),
82 fPtMeasured(obj.fPtMeasured),
83 fEffFlat(obj.fEffFlat),
84 fEfficiency(obj.fEfficiency),
85 fEfficiencyFine(obj.fEfficiencyFine),
86 fResponseMatrix(obj.fResponseMatrix),
87 fResponseMatrixFine(obj.fResponseMatrixFine),
88 fPtMinUnfolded(obj.fPtMinUnfolded),
89 fPtMaxUnfolded(obj.fPtMaxUnfolded),
90 fPtMaxUnfoldedHigh(obj.fPtMaxUnfoldedHigh),
91 fBinWidthFactorUnfolded(obj.fBinWidthFactorUnfolded),
92 fSkipBinsUnfolded(obj.fSkipBinsUnfolded),
93 fExtraBinsUnfolded(obj.fExtraBinsUnfolded),
94 fbVariableBinning(obj.fbVariableBinning),
95 fPtMaxUnfVarBinning(obj.fPtMaxUnfVarBinning),
96 f1MergeFunction(obj.f1MergeFunction),
97 fFineFrac(obj.fFineFrac),
98 fbCalcErrors(obj.fbCalcErrors)
101 //--------------------------------------------------------------------------------------------------------------------------------------------------
102 AliAnaChargedJetResponseMaker& AliAnaChargedJetResponseMaker::operator=(const AliAnaChargedJetResponseMaker& other)
105 if(&other == this) return *this;
106 AliAnaChargedJetResponseMaker::operator=(other);
107 fDebug = other.fDebug;
108 fResolutionType = other.fResolutionType;
109 fDeltaPt = other.fDeltaPt;
110 fhDeltaPt = other.fhDeltaPt;
111 fh1MeasuredTruncated = other.fh1MeasuredTruncated;
112 fh2DetectorResponse = other.fh2DetectorResponse;
113 fh2DetectorResponseRebin = other.fh2DetectorResponseRebin;
114 fhEfficiencyDet = other.fhEfficiencyDet;
115 fh2ResponseMatrixCombinedFineFull = other.fh2ResponseMatrixCombinedFineFull;
116 fh2ResponseMatrixCombinedFull = other.fh2ResponseMatrixCombinedFull;
117 fh2ResponseMatrixCombined = other.fh2ResponseMatrixCombined;
118 fhEfficiencyCombined = other.fhEfficiencyCombined;
119 fDimensions = other.fDimensions;
120 fDimRec = other.fDimRec;
121 fDimGen = other.fDimGen;
122 fPtMin = other.fPtMin;
123 fPtMax = other.fPtMax;
124 fNbins = other.fNbins;
125 fBinArrayPtRec = other.fBinArrayPtRec;
126 fPtMeasured = other.fPtMeasured;
127 fEffFlat = other.fEffFlat;
128 fEfficiency = other.fEfficiency;
129 fEfficiencyFine = other.fEfficiencyFine;
130 fResponseMatrix = other.fResponseMatrix;
131 fResponseMatrixFine = other.fResponseMatrixFine;
132 fPtMinUnfolded = other.fPtMinUnfolded;
133 fPtMaxUnfolded = other.fPtMaxUnfolded;
134 fPtMaxUnfoldedHigh = other.fPtMaxUnfoldedHigh;
135 fBinWidthFactorUnfolded = other.fBinWidthFactorUnfolded;
136 fSkipBinsUnfolded = other.fSkipBinsUnfolded;
137 fExtraBinsUnfolded = other.fExtraBinsUnfolded;
138 fbVariableBinning = other.fbVariableBinning;
139 fPtMaxUnfVarBinning = other.fPtMaxUnfVarBinning;
140 f1MergeFunction = other.f1MergeFunction;
141 fFineFrac = other.fFineFrac;
142 fbCalcErrors = other.fbCalcErrors;
147 //--------------------------------------------------------------------------------------------------------------------------------------------------
148 void AliAnaChargedJetResponseMaker::SetMeasuredSpectrum(TH1D *hPtMeasured)
151 // Set measured spectrum in THnSparse format
152 // This defines the minimum and maximum pT on the reconstructed axis of the response matrix
154 if(fDebug) printf(">>AliAnaChargedJetResponseMaker::SetMeasuredSpectrum \n");
156 fNbins = hPtMeasured->GetXaxis()->GetNbins();
157 fPtMin = hPtMeasured->GetXaxis()->GetXmin();
158 fPtMax = hPtMeasured->GetXaxis()->GetXmax();
160 if(fDebug) printf("fNbins: %d fPtMin: %f fPtMax: %f \n",fNbins,fPtMin,fPtMax);
162 if(fBinArrayPtRec) delete fBinArrayPtRec;
163 fBinArrayPtRec = new Double_t[fNbins+1];
164 for(int j = 0; j<fNbins; j++) {
165 fBinArrayPtRec[j] = hPtMeasured->GetXaxis()->GetBinLowEdge(j+1);
167 fBinArrayPtRec[fNbins] = hPtMeasured->GetXaxis()->GetBinUpEdge(fNbins);
170 Int_t nbins[fDimensions];
171 Double_t xmin[fDimensions];
172 Double_t xmax[fDimensions];
173 for(int dim = 0; dim<fDimensions; dim++) {
179 if(fPtMeasured) delete fPtMeasured;
180 fPtMeasured = new THnSparseD("fPtMeasured","Measured pT spectrum; p_{T}^{rec}",fDimensions,nbins,xmin,xmax);
181 fPtMeasured->Sumw2();
182 fPtMeasured->GetAxis(0)->SetTitle("p_{T}^{rec}");
183 fPtMeasured->SetBinEdges(0,fBinArrayPtRec);
186 if(fDebug) printf("fill measured THnSparse\n");
187 if(fNbins!=hPtMeasured->GetNbinsX())
188 printf("WARNING: nbins not correct \t %d vs %d !!!\n",fNbins,hPtMeasured->GetNbinsX());
192 for(int i = hPtMeasured->FindBin(fPtMin); i<hPtMeasured->FindBin(fPtMax); i++) {
193 fPtMeasured->SetBinContent(bin,(Double_t)hPtMeasured->GetBinContent(i));
194 fPtMeasured->SetBinError(bin,(Double_t)hPtMeasured->GetBinError(i));
198 if(fDebug) printf("fPtMeasured->GetNbins(): %d \n",(int)fPtMeasured->GetNbins());
202 //--------------------------------------------------------------------------------------------------------------------------------------------------
203 void AliAnaChargedJetResponseMaker::SetFlatEfficiency(Double_t eff) {
206 // Set flat efficiency to value of eff
213 Int_t nbins[fDimensions];
214 Double_t xmin[fDimensions];
215 Double_t xmax[fDimensions];
216 for(int dim = 0; dim<fDimensions; dim++) {
222 if(fEfficiency) delete fEfficiency;
223 fEfficiency = new THnSparseD("fEfficiency","Efficiency - no resolution effects; p_{T}^{gen}",fDimensions,nbins,xmin,xmax);
224 fEfficiency->Sumw2();
225 fEfficiency->GetAxis(0)->SetTitle("p_{T}^{gen}");
226 fEfficiency->SetBinEdges(0,fBinArrayPtRec);
228 for(int i=0; i<fNbins; i++) {
231 fEfficiency->SetBinContent(bin,fEffFlat);
232 fEfficiency->SetBinError(bin,0);
237 //--------------------------------------------------------------------------------------------------------------------------------------------------
238 void AliAnaChargedJetResponseMaker::SetEfficiency(TGraphErrors *grEff)
241 // Fill fEfficiency (THnSparse) with values from grEff
244 Int_t nbins[fDimensions];
245 Double_t xmin[fDimensions];
246 Double_t xmax[fDimensions];
247 for(int dim = 0; dim<fDimensions; dim++) {
253 if(fEfficiency) delete fEfficiency;
254 fEfficiency = new THnSparseD("fEfficiency","Efficiency - no resolution effects; p_{T}^{gen}",fDimensions,nbins,xmin,xmax);
255 fEfficiency->Sumw2();
256 fEfficiency->GetAxis(0)->SetTitle("p_{T}^{gen}");
257 fEfficiency->SetBinEdges(0,fBinArrayPtRec);
263 int nbinsgr = grEff->GetN();
265 for(int i=0; i<nbinsgr; i++) {
266 grEff->GetPoint(i,dummy,yield);
268 error = grEff->GetErrorY(i);
270 fEfficiency->Fill(pT,yield);
271 fEfficiency->SetBinError(i,error);
277 //--------------------------------------------------------------------------------------------------------------------------------------------------
278 void AliAnaChargedJetResponseMaker::MakeResponseMatrixJetsFineMerged(Int_t skipBins, Int_t binWidthFactor, Int_t extraBins, Bool_t bVariableBinning, Double_t ptmax)
281 // Make jet response matrix
284 if(fDebug) printf(">>AliAnaChargedJetResponseMaker::MakeResponseMatrixJetsFineMerged\n");
287 if(fDebug) printf("fPtMeasured does not exist. Aborting!!\n");
290 if(fResponseMatrix) { delete fResponseMatrix; }
291 if(fResponseMatrixFine) { delete fResponseMatrixFine; }
293 SetSkipBinsUnfolded(skipBins);
294 SetBinWidthFactorUnfolded(binWidthFactor);
295 SetExtraBinsUnfolded(extraBins);
296 SetVariableBinning(bVariableBinning,ptmax);
298 InitializeResponseMatrix();
299 InitializeEfficiency();
301 InitializeResponseMatrixFine();
302 InitializeEfficiencyFine();
304 FillResponseMatrixFineAndMerge();
308 //--------------------------------------------------------------------------------------------------------------------------------------------------
309 void AliAnaChargedJetResponseMaker::InitializeResponseMatrix() {
311 //Define bin edges of RM to be used for unfolding
314 Int_t nbins[fDimensions*2];
315 nbins[fDimRec] = fNbins;
316 nbins[fDimGen] = fNbins;
318 double binWidthMeas = (double)((fPtMax-fPtMin)/fNbins);
319 double binWidthUnf = (double)fBinWidthFactorUnfolded*binWidthMeas;
320 double binWidthUnfLowPt = -1.;
321 if(fbVariableBinning)
322 binWidthUnfLowPt = binWidthUnf*0.5;
324 fPtMaxUnfolded = fPtMax+(double)(fExtraBinsUnfolded)*binWidthUnf;
325 nbins[fDimGen]+=fExtraBinsUnfolded;
328 printf("fPtMinMeas: %f fPtMaxMeas: %f\n",fPtMin,fPtMax);
329 printf("binWidthMeas: %f binWidthUnf: %f fBinWidthFactorUnfolded: %d\n",binWidthMeas,binWidthUnf,fBinWidthFactorUnfolded);
330 printf("binWidthUnfLowPt: %f\n",binWidthUnfLowPt);
332 printf("fPtMinUnfolded: %f fPtMaxUnfolded: %f\n",fPtMinUnfolded,fPtMaxUnfolded);
335 if(fbVariableBinning) {
336 // cout << "fPtMaxUnfVarBinning: " << fPtMaxUnfVarBinning << " \tfPtMinUnfolded: " << fPtMinUnfolded << " binWidthUnfLowPt: " << binWidthUnfLowPt << endl;
337 Int_t tmp = (int)((fPtMaxUnfVarBinning-fPtMinUnfolded)/binWidthUnfLowPt);
338 tmp = tmp - fSkipBinsUnfolded;
339 fPtMinUnfolded = fPtMaxUnfVarBinning-(double)(tmp)*binWidthUnfLowPt;
340 //cout << "tmp = " << tmp << " fSkipBinsUnfolded = " << fSkipBinsUnfolded << " fPtMinUnfolded = " << fPtMinUnfolded << endl;
341 //Redefine also number of bins on generated axis in case of variable binning
342 nbins[fDimGen] = (int)((fPtMaxUnfVarBinning-fPtMinUnfolded)/binWidthUnfLowPt)+(int)((fPtMaxUnfolded-fPtMaxUnfVarBinning)/binWidthUnf);
345 int tmp = round((fPtMaxUnfolded-fPtMinUnfolded)/binWidthUnf); //nbins which fit between 0 and fPtMaxUnfolded
346 tmp = tmp - fSkipBinsUnfolded;
347 fPtMinUnfolded = fPtMaxUnfolded-(double)(tmp)*binWidthUnf; //set ptmin unfolded
348 fPtMaxUnfolded = fPtMinUnfolded+(double)(tmp)*binWidthUnf; //set ptmax unfolded
349 nbins[fDimGen] = (int)((fPtMaxUnfolded-fPtMinUnfolded)/binWidthUnf); //adjust nbins to bin width
352 printf(" nbins[fDimGen] = %d nbins[fDimRec] = %d\n",nbins[fDimGen],nbins[fDimRec]);
354 Double_t binWidth[2];
355 binWidth[fDimRec] = (double)((fPtMax-fPtMin)/nbins[fDimRec]);
356 binWidth[fDimGen] = (double)((fPtMaxUnfolded-fPtMinUnfolded)/nbins[fDimGen]);
358 Double_t xmin[fDimensions*2];
359 Double_t xmax[fDimensions*2];
360 xmin[fDimRec] = fPtMin;
361 xmax[fDimRec] = fPtMax;
362 xmin[fDimGen] = fPtMinUnfolded;
363 xmax[fDimGen] = fPtMaxUnfolded;
365 printf("xmin[fDimRec]: %f xmin[fDimGen]: %f \n",xmin[fDimRec],xmin[fDimGen]);
366 printf("xmax[fDimRec]: %f xmax[fDimGen]: %f \n",xmax[fDimRec],xmax[fDimGen]);
367 printf("nbins[fDimRec]: %d nbins[fDimGen]: %d \n",nbins[fDimRec],nbins[fDimGen]);
368 if(!fbVariableBinning) printf("binWidth[fDimRec]: %f binWidth[fDimGen]: %f \n",binWidth[fDimRec],binWidth[fDimGen]);
370 Double_t binArrayPt0[nbins[fDimRec]+1];
371 Double_t binArrayPt1[nbins[fDimGen]+1];
372 Double_t xmaxGen = TMath::Max(xmax[fDimGen],fPtMaxUnfoldedHigh);
374 //Define bin limits reconstructed/measured axis
375 for(int i=0; i<nbins[fDimRec]; i++) {
376 binArrayPt0[i] = xmin[fDimRec]+(double)i*binWidth[fDimRec];
378 binArrayPt0[nbins[fDimRec]]= xmax[fDimRec];
380 //Define bin limits generated/unfolded axis
381 binArrayPt1[0] = fPtMinUnfolded;
382 for(int i=1; i<nbins[fDimGen]; i++) {
383 if(fbVariableBinning) {
384 double test = xmin[fDimGen]+(double)i*binWidthUnfLowPt;
385 if(test<=fPtMaxUnfVarBinning) binArrayPt1[i] = test;
386 else binArrayPt1[i] = binArrayPt1[i-1]+binWidthUnf;
389 binArrayPt1[i] = xmin[fDimGen]+(double)i*binWidth[fDimGen];
390 //printf("RM. i = %d \t binArrayPt[i] = %f \n",i,binArrayPt1[i]);
392 binArrayPt1[nbins[fDimGen]]= xmaxGen;
395 // Response matrix : dimensions must be 2N = 2 x (number of variables)
396 // dimensions 0 -> N-1 must be filled with reconstructed values
397 // dimensions N -> 2N-1 must be filled with generated values
398 if(fResponseMatrix) delete fResponseMatrix;
399 fResponseMatrix = new THnSparseD("fResponseMatrix","Response Matrix pTMC vs pTrec",fDimensions*2,nbins,xmin,xmax);
400 fResponseMatrix->Sumw2();
401 fResponseMatrix->GetAxis(fDimRec)->SetTitle("p_{T}^{rec}");
402 fResponseMatrix->GetAxis(fDimGen)->SetTitle("p_{T}^{gen}");
404 fResponseMatrix->SetBinEdges(fDimRec,binArrayPt0);
405 fResponseMatrix->SetBinEdges(fDimGen,binArrayPt1);
407 Int_t bin[2] = {1,1};
408 for(int i=1; i<fResponseMatrix->GetAxis(0)->GetNbins(); i++) {
410 for(int j=1; j<fResponseMatrix->GetAxis(1)->GetNbins(); j++) {
412 fResponseMatrix->SetBinContent(bin,0.);
420 //--------------------------------------------------------------------------------------------------------------------------------------------------
421 void AliAnaChargedJetResponseMaker::InitializeEfficiency() {
423 // Define dimensions of efficiency THnSparse
426 if(!fResponseMatrix) {
427 printf("AliAnaChargedJetResponseMaker::InitializeEfficiency()\n");
428 printf("Can not define dimensions efficiency without fResponseMatrix dimensions. Aborting \n");
432 TAxis *genAxis = fResponseMatrix->GetAxis(fDimGen);
434 Int_t nbinsEff[fDimensions];
435 Double_t xminEff[fDimensions];
436 Double_t xmaxEff[fDimensions];
438 for(int dim = 0; dim<fDimensions; dim++) {
439 nbinsEff[dim] = genAxis->GetNbins();
440 xminEff[dim] = genAxis->GetXmin();
441 xmaxEff[dim] = genAxis->GetXmax();
444 if(fEfficiency) delete fEfficiency;
445 fEfficiency = new THnSparseD("fEfficiency","Efficiency - no resolution effects; p_{T}^{gen}",fDimensions,nbinsEff,xminEff,xmaxEff);
446 fEfficiency->Sumw2();
447 fEfficiency->GetAxis(0)->SetTitle("p_{T}^{gen}");
449 const Double_t *binArrayPt = genAxis->GetXbins()->GetArray();
450 fEfficiency->SetBinEdges(0,binArrayPt);
454 //--------------------------------------------------------------------------------------------------------------------------------------------------
455 void AliAnaChargedJetResponseMaker::InitializeResponseMatrixFine() {
457 //Initialize fine response matrix
460 Int_t nbinsFine[fDimensions*2];
461 Double_t xminFine[fDimensions*2];
462 Double_t xmaxFine[fDimensions*2];
463 // Double_t pTarrayFine[fDimensions*2];
465 nbinsFine[fDimRec] = fResponseMatrix->GetAxis(fDimRec)->GetNbins()*fFineFrac;
466 nbinsFine[fDimGen] = fResponseMatrix->GetAxis(fDimRec)->GetNbins()*fFineFrac;
467 xminFine[fDimRec] = TMath::Min(fPtMin,0.);
468 xminFine[fDimGen] = TMath::Min(fPtMin,0.);
469 xmaxFine[fDimRec] = fResponseMatrix->GetAxis(fDimGen)->GetXmax();//+40.;
470 xmaxFine[fDimGen] = fResponseMatrix->GetAxis(fDimGen)->GetXmax();//+40.;
471 // pTarrayFine[fDimRec] = 0.;
472 // pTarrayFine[fDimGen] = 0.;
474 Double_t binWidth[2];
475 binWidth[fDimRec] = fResponseMatrix->GetAxis(fDimRec)->GetBinWidth(1);
477 Double_t binWidthFine[2];
478 binWidthFine[fDimRec] = binWidth[fDimRec]/((double)fFineFrac);
479 binWidthFine[fDimGen] = binWidthFine[fDimRec]*(double)fBinWidthFactorUnfolded;
480 // binWidthFine[fDimRec] = binWidthFine[fDimGen];
482 nbinsFine[fDimRec] = (int)((xmaxFine[fDimRec]-xminFine[fDimRec])/binWidthFine[fDimRec]); //adjust nbins to bin width
483 nbinsFine[fDimGen] = (int)((xmaxFine[fDimGen]-xminFine[fDimGen])/binWidthFine[fDimGen]); //adjust nbins to bin width
485 printf("xminFine[fDimRec]: %f xminFine[fDimGen]: %f \n",xminFine[fDimRec],xminFine[fDimGen]);
486 printf("xmaxFine[fDimRec]: %f xmaxFine[fDimGen]: %f \n",xmaxFine[fDimRec],xmaxFine[fDimGen]);
487 printf("nbinsFine[fDimRec]: %d nbinsFine[fDimGen]: %d \n",nbinsFine[fDimRec],nbinsFine[fDimGen]);
488 printf("binWidthFine[fDimRec]: %f binWidthFine[fDimGen]: %f \n",binWidthFine[fDimRec],binWidthFine[fDimGen]);
491 Double_t binArrayPt0Fine[nbinsFine[fDimRec]+1];
492 Double_t binArrayPt1Fine[nbinsFine[fDimGen]+1];
494 for(int i=0; i<nbinsFine[fDimRec]; i++) {
495 binArrayPt0Fine[i] = xminFine[fDimRec]+(double)i*binWidthFine[fDimRec];
496 // printf("RM. i = %d \t binArrayPtFine[i] = %f \n",i,binArrayPt0Fine[i]);
498 binArrayPt0Fine[nbinsFine[fDimRec]]= xmaxFine[fDimRec];
500 for(int i=0; i<nbinsFine[fDimGen]; i++) {
501 binArrayPt1Fine[i] = xminFine[fDimGen]+(double)i*binWidthFine[fDimGen];
502 // printf("RM. i = %d \t binArrayPtFine[i] = %f \n",i,binArrayPt1Fine[i]);
504 binArrayPt1Fine[nbinsFine[fDimGen]]= xmaxFine[fDimGen];
506 // Response matrix : dimensions must be 2N = 2 x (number of variables)
507 // dimensions 0 -> N-1 must be filled with reconstructed values
508 // dimensions N -> 2N-1 must be filled with generated values
509 if(fResponseMatrixFine) delete fResponseMatrixFine;
510 fResponseMatrixFine = new THnSparseD("fResponseMatrixFine","Response Matrix pTMC vs pTrec",fDimensions*2,nbinsFine,xminFine,xmaxFine);
511 fResponseMatrixFine->Sumw2();
512 fResponseMatrixFine->GetAxis(fDimRec)->SetTitle("p_{T}^{rec}");
513 fResponseMatrixFine->GetAxis(fDimGen)->SetTitle("p_{T}^{gen}");
515 fResponseMatrixFine->SetBinEdges(fDimRec,binArrayPt0Fine);
516 fResponseMatrixFine->SetBinEdges(fDimGen,binArrayPt1Fine);
518 Int_t bin[2] = {1,1};
519 for(int i=1; i<fResponseMatrixFine->GetAxis(0)->GetNbins(); i++) {
521 for(int j=1; j<fResponseMatrixFine->GetAxis(1)->GetNbins(); j++) {
523 fResponseMatrixFine->SetBinContent(bin,0.);
529 //--------------------------------------------------------------------------------------------------------------------------------------------------
530 void AliAnaChargedJetResponseMaker::InitializeEfficiencyFine() {
532 // Define dimensions of efficiency THnSparse
535 if(!fResponseMatrixFine) {
536 printf("AliAnaChargedJetResponseMaker::InitializeEfficiencyFine()\n");
537 printf("Can not define dimensions efficiency without fResponseMatrixFine dimensions. Aborting \n");
541 TAxis *genAxis = fResponseMatrixFine->GetAxis(fDimGen);
543 Int_t nbinsEff[fDimensions];
544 Double_t xminEff[fDimensions];
545 Double_t xmaxEff[fDimensions];
547 for(int dim = 0; dim<fDimensions; dim++) {
548 nbinsEff[dim] = genAxis->GetNbins();
549 xminEff[dim] = genAxis->GetXmin();
550 xmaxEff[dim] = genAxis->GetXmax();
553 if(fEfficiencyFine) delete fEfficiencyFine;
554 fEfficiencyFine = new THnSparseD("fEfficiencyFine","EfficiencyFine - no resolution effects; p_{T}^{gen}",fDimensions,nbinsEff,xminEff,xmaxEff);
555 fEfficiencyFine->Sumw2();
556 fEfficiencyFine->GetAxis(0)->SetTitle("p_{T}^{gen}");
558 const Double_t *binArrayPt = genAxis->GetXbins()->GetArray();
559 fEfficiencyFine->SetBinEdges(0,binArrayPt);
563 //--------------------------------------------------------------------------------------------------------------------------------------------------
564 void AliAnaChargedJetResponseMaker::FillResponseMatrixFineAndMerge() {
566 // Fill fine response matrix
569 if(!fResponseMatrix) {
570 printf("Dimensions of fResponseMatrix have to be defined first. Aborting!");
574 if(!fResponseMatrixFine) {
575 printf("Dimensions of fResponseMatrixFine have to be defined first. Aborting!");
579 TAxis *genAxis = fResponseMatrixFine->GetAxis(fDimGen);
580 TAxis *recAxis = fResponseMatrixFine->GetAxis(fDimRec);
582 Int_t nbinsFine[fDimensions*2];
583 //Double_t xminFine[fDimensions*2];
584 //Double_t xmaxFine[fDimensions*2];
585 Double_t pTarrayFine[fDimensions*2];
587 nbinsFine[fDimGen] = genAxis->GetNbins();
588 nbinsFine[fDimRec] = recAxis->GetNbins();
589 // xminFine[fDimGen] = genAxis->GetXmin();
590 // xminFine[fDimRec] = recAxis->GetXmin();
591 // xmaxFine[fDimGen] = genAxis->GetXmax();
592 // xmaxFine[fDimRec] = recAxis->GetXmax();
593 pTarrayFine[fDimGen] = 0.;
594 pTarrayFine[fDimRec] = 0.;
596 double sumyield = 0.;
597 double sumyielderror2 = 0.;
600 int iptMerged[2] = {0,0};
601 int ierror[2] = {0,0};
604 Double_t pTgen, pTrec;
608 const int nng = fResponseMatrix->GetAxis(fDimGen)->GetNbins();
609 const int nnr = fResponseMatrix->GetAxis(fDimRec)->GetNbins();
610 Double_t errorArray[nng][nnr];
611 for(int iig =0; iig<nng; iig++) {
612 for(int iir =0; iir<nnr; iir++) {
613 errorArray[iig][iir] = 0.;
617 for(int iy=1; iy<=nbinsFine[fDimGen]; iy++) { //gen
618 pTgen = fResponseMatrixFine->GetAxis(fDimGen)->GetBinCenter(iy);
619 pTarrayFine[fDimGen] = pTgen;
624 for(int ix=1; ix<=nbinsFine[fDimRec]; ix++) { //rec
625 pTrec = fResponseMatrixFine->GetAxis(fDimRec)->GetBinCenter(ix);
626 Double_t width = fResponseMatrixFine->GetAxis(fDimRec)->GetBinWidth(ix);
627 if(fResolutionType==kParam) {
628 yield = fDeltaPt->Eval(pTrec-pTgen);
631 else if(fResolutionType==kResiduals) {
632 yield = fhDeltaPt->Interpolate(pTrec-pTgen);
635 else if(fResolutionType==kResidualsErr) {
636 Double_t deltaPt = pTrec-pTgen;
637 int bin = fhDeltaPt->FindBin(deltaPt);
638 yield = fhDeltaPt->GetBinContent(bin);
639 error = fhDeltaPt->GetBinError(bin);
643 //avoid trouble with empty bins in the high pT tail of deltaPt distribution
644 if(check==0 && yield>0. && pTrec>pTgen) check=1;
645 if(check==1 && yield==0.) ix=nbinsFine[fDimRec];
648 sumyielderror2 += error*error;
650 pTarrayFine[fDimRec] = pTrec;
651 ierror[fDimRec] = ix;
652 fResponseMatrixFine->Fill(pTarrayFine,yield);
653 if(fbCalcErrors) fResponseMatrixFine->SetBinError(ierror,error);
657 //Normalize to total number of counts =1
660 for(int ix=1; ix<=nbinsFine[fDimRec]; ix++) {
662 fResponseMatrixFine->SetBinContent(ipt,fResponseMatrixFine->GetBinContent(ipt)/sumyield);
663 if(fResolutionType==kResidualsErr && fbCalcErrors) {
664 Double_t A = 1./sumyield*fResponseMatrix->GetBinError(ipt);
665 Double_t B = -1.*fResponseMatrix->GetBinContent(ipt)/sumyield/sumyield*TMath::Sqrt(sumyielderror2);
666 Double_t tmp2 = A*A + B*B;
667 fResponseMatrix->SetBinError(ipt,TMath::Sqrt(tmp2));
675 fEfficiencyFine->SetBinContent(bin,sumyield);
676 if(fbCalcErrors) fEfficiencyFine->SetBinError(bin,TMath::Sqrt(sumyielderror2));
678 //fill merged response matrix
680 //find bin in fine RM corresponding to mimimum/maximum bin of merged RM on rec axis
681 int ixMin = fResponseMatrixFine->GetAxis(fDimRec)->FindBin(fResponseMatrix->GetAxis(fDimRec)->GetXmin());
682 int ixMax = fResponseMatrixFine->GetAxis(fDimRec)->FindBin(fResponseMatrix->GetAxis(fDimRec)->GetXmax());
684 if(fResponseMatrixFine->GetAxis(fDimGen)->GetBinLowEdge(iy) >= fResponseMatrix->GetAxis(fDimGen)->GetXmin()) {
686 iptMerged[fDimGen]=fResponseMatrix->GetAxis(fDimGen)->FindBin(pTgen);
688 Double_t weight = 1.;
689 if(f1MergeFunction) {
690 Double_t loEdge = fResponseMatrix->GetAxis(fDimGen)->GetBinLowEdge(iptMerged[fDimGen]);
691 Double_t upEdge = fResponseMatrix->GetAxis(fDimGen)->GetBinUpEdge(iptMerged[fDimGen]);
692 Float_t powInteg = f1MergeFunction->Integral(loEdge,upEdge);
693 //printf("loEdge = %f upEdge = %f powInteg = %f\n",loEdge,upEdge,powInteg);
695 weight = f1MergeFunction->Integral(fResponseMatrixFine->GetAxis(fDimGen)->GetBinLowEdge(iy),fResponseMatrixFine->GetAxis(fDimGen)->GetBinUpEdge(iy))/powInteg;
696 // printf("weight: %f \n", weight );
698 weight = 1./((double)fFineFrac);
699 if(fbVariableBinning && pTgen<fPtMaxUnfVarBinning) weight=1./(0.5*(double)fFineFrac);
702 for(int ix=ixMin; ix<=ixMax; ix++) { //loop reconstructed axis
703 pTrec = fResponseMatrixFine->GetAxis(fDimRec)->GetBinCenter(ix);
705 iptMerged[fDimRec]=fResponseMatrix->GetAxis(fDimRec)->FindBin(pTrec);
707 fResponseMatrix->AddBinContent(iptMerged,fResponseMatrixFine->GetBinContent(ipt)*weight);
708 if(fbCalcErrors) errorArray[iptMerged[fDimGen]-1][iptMerged[fDimRec]-1] += fResponseMatrixFine->GetBinError(ipt)*fResponseMatrixFine->GetBinError(ipt)*weight*weight;
711 }//loop gen axis fine response matrix
713 } // iy (dimGen) loop
715 //Fill Efficiency corresponding to merged response matrix
716 for(int iy=1; iy<=fResponseMatrix->GetAxis(fDimGen)->GetNbins(); iy++) { //gen
720 for(int ix=1; ix<=fResponseMatrix->GetAxis(fDimRec)->GetNbins(); ix++) { //rec
722 sumyield += fResponseMatrix->GetBinContent(ipt);
724 if(fbCalcErrors) fResponseMatrix->SetBinError(ipt,TMath::Sqrt(errorArray[iy][ix]));
726 printf("RM: pTgen: %f \t sumyield: %f \n",fResponseMatrix->GetAxis(fDimGen)->GetBinCenter(iy),sumyield);
729 fEfficiency->SetBinContent(bin,sumyield);
730 if(fbCalcErrors) fEfficiency->SetBinError(bin,0);
733 if(fDebug) printf("fResponseMatrixFine->GetNbins(): %d \n",(int)fResponseMatrixFine->GetNbins());
734 if(fDebug) printf("fResponseMatrix->GetNbins(): %d \n",(int)fResponseMatrix->GetNbins());
736 if(fDebug) printf("Done constructing response matrix\n");
740 //--------------------------------------------------------------------------------------------------------------------------------------------------
741 TH2* AliAnaChargedJetResponseMaker::MakeResponseMatrixRebin(TH2 *hRMFine, TH2 *hRM, Bool_t useFunctionWeight) {
744 // Rebin matrix hRMFine to dimensions of hRM
745 // function returns matrix in TH2D format (hRM2) with dimensions from hRM
748 TH2 *hRM2 = (TH2*)hRM->Clone("hRM2");
751 if(useFunctionWeight) cout << "Use function to do weighting" << endl;
753 //First normalize columns of input
754 const Int_t nbinsNorm = hRM2->GetNbinsX();
755 if(fDebug) cout << "nbinsNorm: " << nbinsNorm << endl;
757 TArrayF *normVector = new TArrayF(nbinsNorm);
759 for(int ix=1; ix<=hRM2->GetNbinsX(); ix++) {
760 Double_t xLow = hRM2->GetXaxis()->GetBinLowEdge(ix);
761 Double_t xUp = hRM2->GetXaxis()->GetBinUpEdge(ix);
762 //cout << "xLow: " << xLow << " xUp: " << xUp << "\t center: " << hRM2->GetXaxis()->GetBinCenter(ix) << endl;
763 Int_t binxLowFine = hRMFine->GetXaxis()->FindBin(xLow);
764 Int_t binxUpFine = hRMFine->GetXaxis()->FindBin(xUp)-1;
765 //cout << "xLowFine: " << hRMFine->GetXaxis()->GetBinLowEdge(binxLowFine) << "\txUpFine: " << hRMFine->GetXaxis()->GetBinUpEdge(binxUpFine) << endl;
766 if(useFunctionWeight)
767 normVector->SetAt(f1MergeFunction->Integral(xLow,xUp),ix-1);
769 normVector->SetAt(hRMFine->Integral(binxLowFine,binxUpFine,1,hRMFine->GetYaxis()->GetNbins()),ix-1);
770 //if(fDebug) cout << "ix norm: " << normVector->At(ix-1) << endl;
773 Double_t content, oldcontent = 0.;
776 Double_t xvalLo, xvalUp, yvalLo, yvalUp;
777 Double_t xmin = hRM2->GetXaxis()->GetXmin();
778 Double_t ymin = hRM2->GetYaxis()->GetXmin();
779 Double_t xmax = hRM2->GetXaxis()->GetXmax();
780 Double_t ymax = hRM2->GetYaxis()->GetXmax();
781 for(int ix=1; ix<=hRMFine->GetXaxis()->GetNbins(); ix++) {
782 xvalLo = hRMFine->GetXaxis()->GetBinLowEdge(ix);
783 xvalUp = hRMFine->GetXaxis()->GetBinUpEdge(ix);
784 if(xvalLo<xmin || xvalUp>xmax) continue;
785 ixNew = hRM2->GetXaxis()->FindBin(hRMFine->GetXaxis()->GetBinCenter(ix));
786 for(int iy=1; iy<=hRMFine->GetYaxis()->GetNbins(); iy++) {
787 yvalLo = hRMFine->GetYaxis()->GetBinLowEdge(iy);
788 yvalUp = hRMFine->GetYaxis()->GetBinUpEdge(iy);
789 if(yvalLo<ymin || yvalUp>ymax) continue;
790 content = hRMFine->GetBinContent(ix,iy);
791 iyNew = hRM2->GetYaxis()->FindBin(hRMFine->GetYaxis()->GetBinCenter(iy));
792 oldcontent = hRM2->GetBinContent(ixNew,iyNew);
794 //if(fDebug) cout << "ixNew: " << ixNew << " " << xvalLo << " iyNew: " << iyNew << " " << yvalLo << " content: " << content << " oldContent: " << oldcontent << " newContent: " << oldcontent+content << endl;
795 Double_t weight = 1.;
796 if(normVector->At(ixNew-1)>0.) {
797 if(useFunctionWeight)
798 weight = f1MergeFunction->Integral(xvalLo,xvalUp)/normVector->At(ixNew-1);
800 weight = 1./normVector->At(ixNew-1);
802 hRM2->SetBinContent(ixNew,iyNew,oldcontent+content*weight);
806 if(normVector) delete normVector;
812 //--------------------------------------------------------------------------------------------------------------------------------------------------
813 TH2* AliAnaChargedJetResponseMaker::CreateTruncated2DHisto(TH2 *h2, Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax) {
816 // Limit axis range of 2D histogram
819 Int_t binMinXh2 = h2->GetXaxis()->FindBin(xmin);
820 if(h2->GetXaxis()->GetBinLowEdge(binMinXh2) < xmin ) binMinXh2++;
821 if(h2->GetXaxis()->GetBinLowEdge(binMinXh2) > xmin ) binMinXh2--;
823 Int_t binMinYh2 = h2->GetYaxis()->FindBin(ymin);
824 if(h2->GetYaxis()->GetBinLowEdge(binMinYh2) < ymin ) binMinYh2++;
825 if(h2->GetYaxis()->GetBinLowEdge(binMinYh2) > ymin ) binMinYh2--;
827 Int_t binMaxXh2 = h2->GetXaxis()->FindBin(xmax);
828 if(h2->GetXaxis()->GetBinUpEdge(binMaxXh2) < xmax ) binMaxXh2++;
829 if(h2->GetXaxis()->GetBinUpEdge(binMaxXh2) > xmax ) binMaxXh2--;
831 Int_t binMaxYh2 = h2->GetYaxis()->FindBin(ymax);
832 if(h2->GetYaxis()->GetBinUpEdge(binMaxYh2) < ymax ) binMaxYh2++;
833 if(h2->GetYaxis()->GetBinUpEdge(binMaxYh2) > ymax ) binMaxYh2--;
835 Int_t nbinsX = binMaxXh2-binMinXh2+1;
836 Int_t nbinsY = binMaxYh2-binMinYh2+1;
838 Double_t *binsX = new Double_t[nbinsX+1];
839 Double_t *binsY = new Double_t[nbinsY+1];
841 for(int ix=1; ix<=nbinsX; ix++)
842 binsX[ix-1] = h2->GetXaxis()->GetBinLowEdge(binMinXh2+ix-1);
843 binsX[nbinsX] = h2->GetXaxis()->GetBinUpEdge(binMaxXh2);
845 for(int iy=1; iy<=nbinsY; iy++)
846 binsY[iy-1] = h2->GetYaxis()->GetBinLowEdge(binMinYh2+iy-1);
847 binsY[nbinsY] = h2->GetYaxis()->GetBinUpEdge(binMaxYh2);
849 TH2 *h2Lim = new TH2D("h2Lim","h2Lim",nbinsX,binsX,nbinsY,binsY);
851 for(int ix=1; ix<=nbinsX; ix++) {
852 // cout << "ix: " << ix << " " << binsX[ix] << endl;
853 for(int iy=1; iy<=nbinsY; iy++) {
854 // cout << "ix: " << ix << " " << binsX[ix] << "\tiy: " << iy << " " << binsY[iy] << endl;
856 double content = h2->GetBinContent(binMinXh2+ix-1,binMinYh2+iy-1);
857 double error = h2->GetBinContent(binMinXh2+ix-1,binMinYh2+iy-1);
858 h2Lim->SetBinContent(ix,iy,content);
859 if(fbCalcErrors) h2Lim->SetBinError(ix,iy,error);
869 //--------------------------------------------------------------------------------------------------------------------------------------------------
870 TH2* AliAnaChargedJetResponseMaker::TruncateAxisRangeResponseMatrix(TH2 *hRMOrig, Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax) {
873 // Limit axis range of response matrix without changing normalization
877 //TH2 *hRMLimit2 = (TH2*)hRMLimit->Clone("hRMLimit2");
879 TH2 *hRMLimit2 = CreateTruncated2DHisto(hRMOrig, xmin, xmax, ymin, ymax);
882 double binCent[2] = {0.,0.};
885 Int_t binOrig[2] = {0};
886 for(int ix=1; ix<=hRMLimit2->GetXaxis()->GetNbins(); ix++) {
887 binCent[0] = hRMLimit2->GetXaxis()->GetBinCenter(ix);
888 binOrig[0] = hRMOrig->GetXaxis()->FindBin(binCent[0]);
889 for(int iy=1; iy<=hRMLimit2->GetYaxis()->GetNbins(); iy++) {
890 binCent[1] = hRMLimit2->GetYaxis()->GetBinCenter(iy);
891 binOrig[1] = hRMOrig->GetYaxis()->FindBin(binCent[1]);
893 content = hRMOrig->GetBinContent(binOrig[0],binOrig[1]);
894 error = hRMOrig->GetBinError(binOrig[0],binOrig[1]);
896 hRMLimit2->SetBinContent(ix,iy,content);
897 if(fbCalcErrors) hRMLimit2->SetBinError(ix,iy,error);
907 //--------------------------------------------------------------------------------------------------------------------------------------------------
908 Bool_t AliAnaChargedJetResponseMaker::CheckInputForCombinedResponse() {
910 Bool_t check = kTRUE;
912 if(!fhDeltaPt) check = kFALSE;
913 if(!fh2DetectorResponse) check = kFALSE;
914 if(!fhEfficiencyDet) check = kFALSE;
915 if(!fh1MeasuredTruncated) check = kFALSE;
916 if(fPtMin<0. || fPtMax<0.) check = kFALSE;
921 //--------------------------------------------------------------------------------------------------------------------------------------------------
922 void AliAnaChargedJetResponseMaker::MakeResponseMatrixCombined(Int_t skipBins, Int_t binWidthFactor, Int_t extraBins, Bool_t bVariableBinning, Double_t ptmin) {
924 //Check if all input histos are there otherwise break
925 if(!CheckInputForCombinedResponse()) {
926 printf("Not all input available..... \n");
930 // Make response bkg fluctuations
931 MakeResponseMatrixJetsFineMerged(skipBins,binWidthFactor,extraBins,bVariableBinning,ptmin);
933 //Get TH2D with dimensions we want final RM
934 TH2D *h2ResponseMatrix = fResponseMatrix->Projection(0,1,"E");
935 h2ResponseMatrix->Reset();
937 //Get fine binned response matrix from bkg fluctuations
938 TH2D *h2ResponseMatrixFine = fResponseMatrixFine->Projection(0,1,"E");
940 // Rebin response detector effects
941 const TArrayD *arrayX = h2ResponseMatrixFine->GetXaxis()->GetXbins();
942 const Double_t *xbinsArray = arrayX->GetArray();
944 TH2D *h2ResponseDetTmp = new TH2D("h2ResponseDetTmp","h2ResponseDetTmp",h2ResponseMatrixFine->GetNbinsX(),xbinsArray,h2ResponseMatrixFine->GetNbinsX(),xbinsArray);
946 fh2DetectorResponseRebin = (TH2D*)MakeResponseMatrixRebin(fh2DetectorResponse,h2ResponseDetTmp,kFALSE);
947 fh2DetectorResponseRebin->SetName("fh2DetectorResponseRebin");
948 fh2DetectorResponseRebin->SetTitle("fh2DetectorResponseRebin");
949 fh2DetectorResponseRebin->SetXTitle("p_{T}^{jet} gen");
950 fh2DetectorResponseRebin->SetYTitle("p_{T}^{jet} rec");
952 // Make full combined response matrix fine binning (background flucutuations + detector effects)
953 fh2ResponseMatrixCombinedFineFull = (TH2D*)MultiplityResponseMatrices(h2ResponseMatrixFine,fh2DetectorResponseRebin);
954 fh2ResponseMatrixCombinedFineFull->SetName("fh2ResponseMatrixCombinedFineFull");
955 fh2ResponseMatrixCombinedFineFull->SetTitle("fh2ResponseMatrixCombinedFineFull");
957 // Rebin full combined response matrix with weighting
958 fh2ResponseMatrixCombinedFull = (TH2D*)MakeResponseMatrixRebin(fh2ResponseMatrixCombinedFineFull,h2ResponseMatrix,kTRUE);
959 fh2ResponseMatrixCombinedFull->SetName("fh2ResponseMatrixCombinedFull");
960 fh2ResponseMatrixCombinedFull->SetTitle("fh2ResponseMatrixCombinedFull");
961 fh2ResponseMatrixCombinedFull->SetXTitle("#it{p}_{T,gen}^{ch jet} (GeV/#it{c})");
962 fh2ResponseMatrixCombinedFull->SetYTitle("#it{p}_{T,rec}^{ch jet} (GeV/#it{c})");
964 fh2ResponseMatrixCombined = (TH2D*)TruncateAxisRangeResponseMatrix(fh2ResponseMatrixCombinedFull, h2ResponseMatrix->GetXaxis()->GetXmin(),h2ResponseMatrix->GetXaxis()->GetXmax(),fh1MeasuredTruncated->GetXaxis()->GetXmin(),fh1MeasuredTruncated->GetXaxis()->GetXmax());
965 fh2ResponseMatrixCombined->SetTitle("fh2ResponseMatrixCombined");
966 fh2ResponseMatrixCombined->SetName("fh2ResponseMatrixCombined");
968 fhEfficiencyCombined = (TH1D*)fh2ResponseMatrixCombined->ProjectionX("fhEfficiencyCombined");
969 fhEfficiencyCombined->Reset();
971 for(int i=1; i<=fh2ResponseMatrixCombined->GetNbinsX(); i++) {
972 Double_t sumyield = 0.;
973 for(int j=1; j<=fh2ResponseMatrixCombined->GetYaxis()->GetNbins(); j++) {
974 sumyield+=fh2ResponseMatrixCombined->GetBinContent(i,j);
976 Double_t kCounter = 0.;
977 Double_t kPtLowEdge = fh2ResponseMatrixCombined->GetXaxis()->GetBinLowEdge(i);
978 Double_t kPtUpEdge = fh2ResponseMatrixCombined->GetXaxis()->GetBinUpEdge(i);
979 Double_t kJetEffDet = 0.;
981 for(int k=1; k<=fhEfficiencyDet->GetNbinsX(); k++) {
982 if(fhEfficiencyDet->GetXaxis()->GetBinLowEdge(k)>=kPtLowEdge && fhEfficiencyDet->GetXaxis()->GetBinUpEdge(k)<=kPtUpEdge) {
983 kJetEffDet+=fhEfficiencyDet->GetBinContent(k);
987 kJetEffDet = kJetEffDet/kCounter;
989 if(kJetEffDet==0.) kJetEffDet=1.;
991 fhEfficiencyCombined->SetBinContent(i,sumyield*kJetEffDet);
997 //--------------------------------------------------------------------------------------------------------------------------------------------------
998 TH2* AliAnaChargedJetResponseMaker::MultiplityResponseMatrices(TH2 *h2RMDeltaPt, TH2 *h2RMDetector) {
1001 // Make combined response matrix (background flucutuations + detector effects)
1003 // hRMDeltaPt is the response matrix for background fluctuations from \delta(p_t) measurement
1004 // 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
1006 // Function assumes that generated/unfolded axis is x-axis and reconstructed is on y-axis on both respone matrices
1009 TH2D *h2ResponseMatrixCombined = (TH2D*)h2RMDeltaPt->Clone("h2ResponseMatrixCombined"); //h2RMDeltaPt is the bkg fluctuations RM which has the dimensions we want for the combined response matrix
1010 h2ResponseMatrixCombined->SetTitle("h2ResponseMatrixCombined");
1011 h2ResponseMatrixCombined->SetName("h2ResponseMatrixCombined");
1013 // M = RM_deltaPt * RM_DetEffects * T (M=measured T=truth)
1015 // mx1 = mxd * dxt * tx1
1016 // m = measured bins
1018 // d = rec for RM_DetEffects and gen for RM_deltaPt
1019 // RM_comb = RM_deltaPt * RM_DetEffects (dimensions mxt)
1020 // RM_comb(m,t) = Sum_d RM_deltaPt(m,d)*RM_DetEffects(d,t)
1023 printf("Nt=%d\n",h2ResponseMatrixCombined->GetNbinsX());
1024 printf("Nm=%d\n",h2ResponseMatrixCombined->GetNbinsY());
1025 printf("Nd=%d\n",h2RMDeltaPt->GetNbinsX());
1029 for(Int_t t=1; t<=h2ResponseMatrixCombined->GetNbinsX();t++) {
1030 for(Int_t m=1; m<=h2ResponseMatrixCombined->GetNbinsY();m++) {
1031 Double_t valueSum = 0.;
1032 for(Int_t d=1; d<=h2RMDeltaPt->GetNbinsX();d++) {
1033 valueSum += h2RMDeltaPt->GetBinContent(d,m) * h2RMDetector->GetBinContent(t,d);
1034 // if(t==10 && m==10) cout << "sum m,d=" << m << "," << d << endl;
1036 // if(t==10) cout << "t,m = " << t << "," << m << "\tvalueSum: " << valueSum << endl;
1037 h2ResponseMatrixCombined->SetBinContent(t,m,valueSum);
1041 return h2ResponseMatrixCombined;
1045 //--------------------------------------------------------------------------------------------------------------------------------------------------
1046 TH2* AliAnaChargedJetResponseMaker::GetTransposeResponsMatrix(TH2 *h2RM) {
1048 // Transpose response matrix
1051 Printf("AliAnaChargedJetResponseMaker::GetTransposeResponsMatrix");
1053 //Initialize transposed response matrix h2RMTranspose
1054 // TArrayD *arrayX = (TArrayD*)h2RM->GetXaxis()->GetXbins();
1055 // for(int i=0; i<arrayX->GetSize(); i++) Printf("i: %d %f", i,arrayX->At(i));
1056 // Double_t *xbinsArray = arrayX->GetArray();
1058 // TArrayD *arrayY = (TArrayD*)h2RM->GetYaxis()->GetXbins();
1059 // for(int i=0; i<arrayY->GetSize(); i++) Printf("i: %d %f", i,arrayY->At(i));
1060 // Double_t *ybinsArray = arrayY->GetArray();
1063 Double_t *ybinsArray = new Double_t[h2RM->GetNbinsY()+1];
1064 for(int i=1; i<=h2RM->GetNbinsY(); i++) {
1065 ybinsArray[i-1] = h2RM->GetYaxis()->GetBinLowEdge(i);
1066 Printf("i=%d %f %f",i,ybinsArray[i-1],h2RM->GetYaxis()->GetBinLowEdge(i));
1068 ybinsArray[h2RM->GetNbinsY()] = h2RM->GetYaxis()->GetBinUpEdge(h2RM->GetNbinsY());
1070 Double_t *xbinsArray = new Double_t[h2RM->GetNbinsX()+1];
1071 for(int i=1; i<=h2RM->GetNbinsX(); i++)
1072 xbinsArray[i-1] = h2RM->GetXaxis()->GetBinLowEdge(i);
1073 xbinsArray[h2RM->GetNbinsX()] = h2RM->GetXaxis()->GetBinUpEdge(h2RM->GetNbinsX());
1076 TH2D *h2RMTranspose = new TH2D("h2RMTranspose","h2RMTranspose",h2RM->GetNbinsY(),ybinsArray,h2RM->GetNbinsX(),xbinsArray);
1078 //Fill tranposed response matrix
1079 for (Int_t ibin = 1; ibin <= h2RMTranspose->GetNbinsX(); ibin++) {
1080 for (Int_t jbin = 1; jbin <= h2RMTranspose->GetNbinsY(); jbin++) {
1081 h2RMTranspose->SetBinContent(ibin,jbin, h2RM->GetBinContent(jbin,ibin));
1086 return h2RMTranspose;
1090 //--------------------------------------------------------------------------------------------------------------------------------------------------
1091 TH2* AliAnaChargedJetResponseMaker::NormalizeResponsMatrixYaxisWithPrior(TH2 *h2RM, TH1 *hPrior) {
1093 // Normalize such that the Y projection is the prior
1096 // TH1D *hProjRespY = (TH1D*)h2RM->ProjectionY("hProjRespY");
1097 double intPrior = hPrior->Integral();//"width");
1098 for (Int_t jbin = 1; jbin <= h2RM->GetNbinsY(); jbin++) {
1099 // double corr = hPrior->GetBinContent(jbin)/hProjRespY->GetBinContent(jbin);
1100 for (Int_t ibin = 1; ibin <= h2RM->GetNbinsX(); ibin++) {
1101 double content = h2RM->GetBinContent(ibin,jbin);
1102 // h2RM->SetBinContent(ibin,jbin,content*corr);
1103 h2RM->SetBinContent(ibin,jbin,hPrior->GetBinContent(jbin)/hPrior->GetBinWidth(jbin)/intPrior*content);
1110 //--------------------------------------------------------------------------------------------------------------------------------------------------
1111 TH1D* AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TH1 *hGen, TH2 *hResponse,TH1 *hEfficiency,Bool_t bDrawSlices) {
1114 // Multiply hGen with response matrix to obtain refolded spectrum
1115 // Efficiency must be given. In case efficiency is 1 use SetFlatEfficiency(1.) before calling function
1119 // printf("Efficiency must be given. In case efficiency is 1 use SetFlatEfficiency(1.) before calling function. Aborting!");
1121 printf("Setting efficiency to 1 \n");
1122 hEfficiency = (TH1D*)hGen->Clone("hEfficiency");
1123 hEfficiency->Reset();
1124 for(int i=1; i<=hEfficiency->GetNbinsX(); i++) hEfficiency->SetBinContent(i,1.);
1129 //y-axis: reconstructed
1130 if(fDebug>0) cout << "nbins hGen: " << hGen->GetNbinsX() << "\t nbins hResponseGen: " << hResponse->GetXaxis()->GetNbins() << "\t nbins hResponseRec: " << hResponse->GetYaxis()->GetNbins() << endl;
1132 TH1D *hRec = hResponse->ProjectionY("hRec");
1135 hRec->SetTitle("hRec");
1136 hRec->SetName("hRec");
1138 for(int irec=1; irec<=hRec->GetNbinsX(); irec++)
1139 hRec->SetBinContent(irec,0);
1141 TH1D *hRecGenBin = 0x0;
1142 TCanvas *cSlices = 0x0;
1144 cSlices = new TCanvas("cSlices","cSlices:Slices",400,400);
1148 Double_t yieldMC = 0.;
1149 Double_t yieldMCerror = 0.;
1150 Double_t sumYield = 0.;
1151 const Int_t nbinsRec = hRec->GetNbinsX();
1152 Double_t sumError2[nbinsRec+1];
1153 for(int irec=0; irec<=hRec->GetNbinsX(); irec++)
1157 for(int igen=1; igen<=hGen->GetNbinsX(); igen++) {
1163 eff = hEfficiency->GetBinContent(igen);
1164 yieldMC = hGen->GetBinContent(igen)*eff;
1165 yieldMCerror = hGen->GetBinError(igen)*eff;
1168 hRecGenBin = hResponse->ProjectionY(Form("hRecGenBin%d",igen));
1169 hRecGenBin->Sumw2();
1170 hRecGenBin->Reset();
1171 hRecGenBin->SetTitle(Form("hRecGenBin%d",igen));
1172 hRecGenBin->SetName(Form("hRecGenBin%d",igen));
1175 for(int irec=1; irec<=hRec->GetNbinsX(); irec++) {
1176 hRec->AddBinContent(irec,yieldMC*hResponse->GetBinContent(igen,irec));
1177 sumYield+=hResponse->GetBinContent(igen,irec);
1178 // Double_t A = yieldMC*hResponse->GetBinError(igen,irec);
1179 Double_t B = hResponse->GetBinContent(igen,irec)*yieldMCerror;
1180 // Double_t tmp2 = A*A + B*B;
1181 //sumError2[irec-1] += tmp2 ;
1182 sumError2[irec-1] += B*B;
1185 hRecGenBin->SetBinContent(irec,yieldMC*hResponse->GetBinContent(igen,irec));
1190 hRecGenBin->SetLineColor(igen);
1191 if(igen==1) hRecGenBin->DrawCopy();
1192 else hRecGenBin->DrawCopy("same");
1195 if(hRecGenBin) delete hRecGenBin;
1197 cout << "igen: " << igen << "\tpTMC: " << hGen->GetXaxis()->GetBinCenter(igen) << "\teff:" << eff << "\tsumYield: " << sumYield << endl;
1200 for(int i=0; i<=nbinsRec; i++) {
1202 hRec->SetBinError(i+1,TMath::Sqrt(sumError2[i]));
1209 //--------------------------------------------------------------------------------------------------------------------------------------------------
1210 TH1D* AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TF1 *fGen, TH2 *hResponse,TH1 *hEfficiency) {
1213 // Multiply fGen function with response matrix to obtain (re)folded spectrum
1214 // Efficiency must be given. In case efficiency is 1 use SetFlatEfficiency(1.) before calling function
1219 //y-axis: reconstructed
1221 if(fDebug>0) printf(">>AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TF1 *fGen, TH2 *hResponse,TH1 *hEfficiency)\n");
1223 TH1D *hRec = hResponse->ProjectionY("hRec");
1226 hRec->SetTitle("hRec");
1227 hRec->SetName("hRec");
1229 // TH1D *hRec = new TH1D("hRec","hRec",hResponse->GetNbinsY(),hResponse->GetYaxis()->GetXmin(),hResponse->GetYaxis()->GetXmax());
1231 for(int irec=1; irec<=hRec->GetNbinsX(); irec++)
1232 hRec->SetBinContent(irec,0);
1234 Double_t yieldMC = 0.;
1235 Double_t sumYield = 0.;
1237 for(int igen=1; igen<=hResponse->GetNbinsX(); igen++) {
1240 double pTMC = hResponse->GetXaxis()->GetBinCenter(igen);
1242 int binEff = hEfficiency->FindBin(pTMC);
1243 eff = hEfficiency->GetBinContent(binEff);
1246 // yieldMC = fGen->Eval(pTMC)*eff;
1247 yieldMC = fGen->Integral(hResponse->GetXaxis()->GetBinLowEdge(igen),hResponse->GetXaxis()->GetBinUpEdge(igen))*eff;
1248 for(int irec=1; irec<=hResponse->GetNbinsY(); irec++) {
1249 hRec->AddBinContent(irec,yieldMC*hResponse->GetBinContent(igen,irec));
1250 sumYield+=hResponse->GetBinContent(igen,irec);
1252 // cout << "igen: " << igen << "\tpTMC: " << pTMC << "\tsumYield: " << sumYield << endl;
1253 // cout << "yieldMC: " << yieldMC << endl;
1260 //--------------------------------------------------------------------------------------------------------------------------------------------------
1261 /* static */ Double_t AliAnaChargedJetResponseMaker::GetBetaPerDOFValue(Int_t betaColl, Int_t betaOpt) {
1263 Float_t betaPerDOFoptions[6] = {0.};
1264 //define beta per degree of freedom (DOF=nbins in measured)
1266 betaPerDOFoptions[0] = 9.1e-9;
1267 betaPerDOFoptions[1] = 2.25e-8;
1268 betaPerDOFoptions[2] = 4.55e-8;
1269 betaPerDOFoptions[3] = 9.1e-8;
1270 betaPerDOFoptions[4] = 2.25e-7;
1271 betaPerDOFoptions[5] = 4.55e-7;
1274 betaPerDOFoptions[0] = 9.1e-7;
1275 betaPerDOFoptions[1] = 2.25e-6;
1276 betaPerDOFoptions[2] = 4.55e-6;
1277 betaPerDOFoptions[3] = 9.1e-6;
1278 betaPerDOFoptions[4] = 2.25e-5;
1279 betaPerDOFoptions[5] = 4.55e-5;
1282 betaPerDOFoptions[0] = 9.1e-5;
1283 betaPerDOFoptions[1] = 2.25e-4;
1284 betaPerDOFoptions[2] = 4.55e-4;
1285 betaPerDOFoptions[3] = 9.1e-4;
1286 betaPerDOFoptions[4] = 2.25e-3;
1287 betaPerDOFoptions[5] = 4.55e-3;
1290 betaPerDOFoptions[0] = 9.1e-3;
1291 betaPerDOFoptions[1] = 2.25e-2;
1292 betaPerDOFoptions[2] = 4.55e-2;
1293 betaPerDOFoptions[3] = 9.1e-2;
1294 betaPerDOFoptions[4] = 2.25e-1;
1295 betaPerDOFoptions[5] = 4.55e-1;
1298 betaPerDOFoptions[0] = 9.1e-1;
1299 betaPerDOFoptions[1] = 2.25;
1300 betaPerDOFoptions[2] = 4.55;
1301 betaPerDOFoptions[3] = 9.1;
1302 betaPerDOFoptions[4] = 22.5;
1303 betaPerDOFoptions[5] = 45.5;
1306 return betaPerDOFoptions[betaOpt];
1310 //--------------------------------------------------------------------------------------------------------------------------------------------------
1311 Double_t AliAnaChargedJetResponseMaker::InterpolateFast(TGraph *gr, Double_t x) {
1313 Double_t ipmax = gr->GetN()-1.;
1314 Double_t x2,y2,xold,yold;
1316 Double_t xmin,ymin,xmax, ymax;
1317 gr->GetPoint(0,xmin,ymin);
1318 gr->GetPoint(gr->GetN()-1,xmax,ymax);
1319 if(x<xmin || x>xmax) return 0;
1321 Double_t ip = ipmax/2.;
1322 Double_t ipold = 0.;
1323 gr->GetPoint((int)(ip),x2,y2);
1331 gr->GetPoint((int)(ip),x2,y2);
1335 // cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl;
1340 if(i==0) ipold=ipmax;
1341 ip += (ipold-ip)/2.;
1342 gr->GetPoint((int)(ip),x2,y2);
1343 if(x2>x) ipold = ip;
1346 // cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl;
1353 gr->GetPoint(ip2,x2,y2);
1356 gr->GetPoint(ip2,x2,y2);
1358 gr->GetPoint(ip2+1,xold,yold);
1362 gr->GetPoint(ip2,x2,y2);
1365 gr->GetPoint(ip2,x2,y2);
1367 gr->GetPoint(ip2-1,xold,yold);
1370 // cout << "For x=" << x << " interpolate between: " << xold << " and " << x2 << endl;
1371 return ((x-xold)*y2 + (x2-x)*yold) / (x2-xold);
1375 //--------------------------------------------------------------------------------------------------------------------------------------------------
1376 Double_t AliAnaChargedJetResponseMaker::InterpolateFast(TH1 *h, Double_t x) {
1378 // Double_t ipmax = gr->GetN()-1.;
1379 Double_t ipmax = (double)h->GetNbinsX();
1380 Double_t x2,y2,xold,yold;
1382 Double_t xmin = h->GetXaxis()->GetXmin();
1383 Double_t xmax = h->GetXaxis()->GetXmax();
1384 if(x<xmin || x>xmax) return 0;
1386 Double_t ip = ipmax/2.;
1387 Double_t ipold = 0.;
1389 x2 = h->GetBinCenter((int)ip);
1390 y2 = h->GetBinContent((int)ip);
1398 x2 = h->GetBinCenter((int)ip);
1399 y2 = h->GetBinContent((int)ip);
1403 // cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl;
1408 if(i==0) ipold=ipmax;
1409 ip += (ipold-ip)/2.;
1410 x2 = h->GetBinCenter((int)ip);
1411 y2 = h->GetBinContent((int)ip);
1412 if(x2>x) ipold = ip;
1415 // cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl;
1422 x2 = h->GetBinCenter(ip2);
1423 y2 = h->GetBinContent(ip2);
1426 x2 = h->GetBinCenter(ip2);
1427 y2 = h->GetBinContent(ip2);
1429 xold = h->GetBinCenter(ip2+1);
1430 yold = h->GetBinContent(ip2+1);
1434 x2 = h->GetBinCenter(ip2);
1435 y2 = h->GetBinContent(ip2);
1438 x2 = h->GetBinCenter(ip2);
1439 y2 = h->GetBinContent(ip2);
1441 xold = h->GetBinCenter(ip2-1);
1442 yold = h->GetBinContent(ip2-1);
1445 // cout << "For x=" << x << " interpolate between: " << xold << " and " << x2 << endl;
1446 return ((x-xold)*y2 + (x2-x)*yold) / (x2-xold);