]>
Commit | Line | Data |
---|---|---|
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 | 16 | using std::cout; |
17 | using std::endl; | |
18 | ||
31d3e4f0 | 19 | ClassImp(AliAnaChargedJetResponseMaker) |
20 | ||
21 | AliAnaChargedJetResponseMaker::AliAnaChargedJetResponseMaker(): | |
22 | fDebug(kFALSE), | |
23 | fResolutionType(kParam), | |
24 | fDeltaPt(0x0), | |
25 | fhDeltaPt(0x0), | |
0bf6381d | 26 | fh1MeasuredTruncated(0x0), |
27 | fh2DetectorResponse(0x0), | |
28 | fh2DetectorResponseRebin(0x0), | |
29 | fhEfficiencyDet(0x0), | |
30 | fh2ResponseMatrixCombinedFineFull(0x0), | |
31 | fh2ResponseMatrixCombinedFull(0x0), | |
32 | fh2ResponseMatrixCombined(0x0), | |
33 | fhEfficiencyCombined(0x0), | |
31d3e4f0 | 34 | fDimensions(1), |
35 | fDimRec(0), | |
36 | fDimGen(1), | |
37 | fPtMin(-999), | |
38 | fPtMax(-999), | |
39 | fNbins(0), | |
40 | fBinArrayPtRec(0x0), | |
41 | fPtMeasured(0x0), | |
42 | fEffFlat(0), | |
43 | fEfficiency(0x0), | |
44 | fEfficiencyFine(0x0), | |
45 | fResponseMatrix(0x0), | |
46 | fResponseMatrixFine(0x0), | |
47 | fPtMinUnfolded(0.), | |
48 | fPtMaxUnfolded(0.), | |
49 | fPtMaxUnfoldedHigh(-1.), | |
50 | fBinWidthFactorUnfolded(2), | |
51 | fSkipBinsUnfolded(0), | |
52 | fExtraBinsUnfolded(5), | |
53 | fbVariableBinning(kFALSE), | |
54 | fPtMaxUnfVarBinning(0), | |
55 | f1MergeFunction(0x0), | |
56 | fFineFrac(10), | |
57 | fbCalcErrors(kFALSE) | |
58 | {;} | |
59 | ||
31d3e4f0 | 60 | |
fa7c34ba | 61 | //-------------------------------------------------------------------------------------------------------------------------------------------------- |
62 | AliAnaChargedJetResponseMaker::AliAnaChargedJetResponseMaker(const AliAnaChargedJetResponseMaker& obj): | |
63 | fDebug(obj.fDebug), | |
64 | fResolutionType(obj.fResolutionType), | |
65 | fDeltaPt(obj.fDeltaPt), | |
66 | fhDeltaPt(obj.fhDeltaPt), | |
0bf6381d | 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), | |
fa7c34ba | 75 | fDimensions(obj.fDimensions), |
76 | fDimRec(obj.fDimRec), | |
77 | fDimGen(obj.fDimGen), | |
78 | fPtMin(obj.fPtMin), | |
79 | fPtMax(obj.fPtMax), | |
80 | fNbins(obj.fNbins), | |
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) | |
99 | {;} | |
100 | ||
101 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
102 | AliAnaChargedJetResponseMaker& AliAnaChargedJetResponseMaker::operator=(const AliAnaChargedJetResponseMaker& other) | |
103 | { | |
104 | // Assignment | |
105 | if(&other == this) return *this; | |
106 | AliAnaChargedJetResponseMaker::operator=(other); | |
0bf6381d | 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; | |
fa7c34ba | 143 | |
144 | return *this; | |
31d3e4f0 | 145 | } |
146 | ||
147 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
148 | void AliAnaChargedJetResponseMaker::SetMeasuredSpectrum(TH1D *hPtMeasured) | |
149 | { | |
150 | // | |
151 | // Set measured spectrum in THnSparse format | |
152 | // This defines the minimum and maximum pT on the reconstructed axis of the response matrix | |
153 | // | |
154 | if(fDebug) printf(">>AliAnaChargedJetResponseMaker::SetMeasuredSpectrum \n"); | |
155 | ||
156 | fNbins = hPtMeasured->GetXaxis()->GetNbins(); | |
157 | fPtMin = hPtMeasured->GetXaxis()->GetXmin(); | |
158 | fPtMax = hPtMeasured->GetXaxis()->GetXmax(); | |
159 | ||
160 | if(fDebug) printf("fNbins: %d fPtMin: %f fPtMax: %f \n",fNbins,fPtMin,fPtMax); | |
161 | ||
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); | |
166 | } | |
167 | fBinArrayPtRec[fNbins] = hPtMeasured->GetXaxis()->GetBinUpEdge(fNbins); | |
168 | ||
169 | ||
170 | Int_t nbins[fDimensions]; | |
171 | Double_t xmin[fDimensions]; | |
172 | Double_t xmax[fDimensions]; | |
173 | for(int dim = 0; dim<fDimensions; dim++) { | |
174 | nbins[dim] = fNbins; | |
175 | xmin[dim] = fPtMin; | |
176 | xmax[dim] = fPtMax; | |
177 | } | |
178 | ||
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); | |
184 | ||
185 | //Fill | |
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()); | |
189 | ||
190 | int bin[1] = {0}; | |
191 | bin[0] = 0; | |
192 | for(int i = hPtMeasured->FindBin(fPtMin); i<hPtMeasured->FindBin(fPtMax); i++) { | |
31d3e4f0 | 193 | fPtMeasured->SetBinContent(bin,(Double_t)hPtMeasured->GetBinContent(i)); |
194 | fPtMeasured->SetBinError(bin,(Double_t)hPtMeasured->GetBinError(i)); | |
195 | bin[0]++; | |
196 | } | |
197 | ||
198 | if(fDebug) printf("fPtMeasured->GetNbins(): %d \n",(int)fPtMeasured->GetNbins()); | |
199 | ||
200 | } | |
201 | ||
202 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
203 | void AliAnaChargedJetResponseMaker::SetFlatEfficiency(Double_t eff) { | |
204 | ||
ef62323a | 205 | // |
206 | // Set flat efficiency to value of eff | |
207 | // | |
208 | ||
31d3e4f0 | 209 | fEffFlat = eff; |
ef62323a | 210 | return; |
31d3e4f0 | 211 | |
275d3481 | 212 | /* |
31d3e4f0 | 213 | Int_t nbins[fDimensions]; |
214 | Double_t xmin[fDimensions]; | |
215 | Double_t xmax[fDimensions]; | |
216 | for(int dim = 0; dim<fDimensions; dim++) { | |
217 | nbins[dim] = fNbins; | |
218 | xmin[dim] = fPtMin; | |
219 | xmax[dim] = fPtMax; | |
220 | } | |
221 | ||
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); | |
227 | ||
228 | for(int i=0; i<fNbins; i++) { | |
229 | int bin[1]; | |
230 | bin[0] = i; | |
231 | fEfficiency->SetBinContent(bin,fEffFlat); | |
232 | fEfficiency->SetBinError(bin,0); | |
233 | } | |
5b045191 | 234 | */ |
31d3e4f0 | 235 | } |
236 | ||
237 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
238 | void AliAnaChargedJetResponseMaker::SetEfficiency(TGraphErrors *grEff) | |
239 | { | |
ef62323a | 240 | // |
241 | // Fill fEfficiency (THnSparse) with values from grEff | |
242 | // | |
243 | ||
31d3e4f0 | 244 | Int_t nbins[fDimensions]; |
245 | Double_t xmin[fDimensions]; | |
246 | Double_t xmax[fDimensions]; | |
247 | for(int dim = 0; dim<fDimensions; dim++) { | |
248 | nbins[dim] = fNbins; | |
249 | xmin[dim] = fPtMin; | |
250 | xmax[dim] = fPtMax; | |
251 | } | |
252 | ||
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); | |
258 | ||
259 | double pT[1]; | |
260 | double yield = 0.; | |
261 | double error = 0.; | |
262 | double dummy = 0.; | |
263 | int nbinsgr = grEff->GetN(); | |
264 | ||
265 | for(int i=0; i<nbinsgr; i++) { | |
266 | grEff->GetPoint(i,dummy,yield); | |
267 | pT[0] = dummy; | |
268 | error = grEff->GetErrorY(i); | |
269 | ||
270 | fEfficiency->Fill(pT,yield); | |
271 | fEfficiency->SetBinError(i,error); | |
272 | ||
273 | } | |
274 | ||
275 | } | |
276 | ||
277 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
278 | void AliAnaChargedJetResponseMaker::MakeResponseMatrixJetsFineMerged(Int_t skipBins, Int_t binWidthFactor, Int_t extraBins, Bool_t bVariableBinning, Double_t ptmax) | |
279 | { | |
280 | // | |
281 | // Make jet response matrix | |
282 | // | |
283 | ||
ef62323a | 284 | if(fDebug) printf(">>AliAnaChargedJetResponseMaker::MakeResponseMatrixJetsFineMerged\n"); |
285 | ||
286 | if(!fPtMeasured) { | |
287 | if(fDebug) printf("fPtMeasured does not exist. Aborting!!\n"); | |
288 | return; | |
289 | } | |
31d3e4f0 | 290 | if(fResponseMatrix) { delete fResponseMatrix; } |
291 | if(fResponseMatrixFine) { delete fResponseMatrixFine; } | |
292 | ||
293 | SetSkipBinsUnfolded(skipBins); | |
294 | SetBinWidthFactorUnfolded(binWidthFactor); | |
295 | SetExtraBinsUnfolded(extraBins); | |
296 | SetVariableBinning(bVariableBinning,ptmax); | |
297 | ||
298 | InitializeResponseMatrix(); | |
299 | InitializeEfficiency(); | |
300 | ||
301 | InitializeResponseMatrixFine(); | |
302 | InitializeEfficiencyFine(); | |
303 | ||
304 | FillResponseMatrixFineAndMerge(); | |
305 | ||
306 | } | |
307 | ||
308 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
309 | void AliAnaChargedJetResponseMaker::InitializeResponseMatrix() { | |
310 | // | |
ef62323a | 311 | //Define bin edges of RM to be used for unfolding |
31d3e4f0 | 312 | // |
313 | ||
314 | Int_t nbins[fDimensions*2]; | |
315 | nbins[fDimRec] = fNbins; | |
fa7c34ba | 316 | nbins[fDimGen] = fNbins; |
31d3e4f0 | 317 | |
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; | |
323 | ||
03372fd1 | 324 | fPtMaxUnfolded = fPtMax+(double)(fExtraBinsUnfolded)*binWidthUnf; |
325 | nbins[fDimGen]+=fExtraBinsUnfolded; | |
326 | ||
31d3e4f0 | 327 | |
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); | |
331 | ||
adf3920d | 332 | printf("fPtMinUnfolded: %f fPtMaxUnfolded: %f\n",fPtMinUnfolded,fPtMaxUnfolded); |
333 | ||
31d3e4f0 | 334 | |
335 | if(fbVariableBinning) { | |
adf3920d | 336 | // cout << "fPtMaxUnfVarBinning: " << fPtMaxUnfVarBinning << " \tfPtMinUnfolded: " << fPtMinUnfolded << " binWidthUnfLowPt: " << binWidthUnfLowPt << endl; |
337 | Int_t tmp = (int)((fPtMaxUnfVarBinning-fPtMinUnfolded)/binWidthUnfLowPt); | |
31d3e4f0 | 338 | tmp = tmp - fSkipBinsUnfolded; |
339 | fPtMinUnfolded = fPtMaxUnfVarBinning-(double)(tmp)*binWidthUnfLowPt; | |
adf3920d | 340 | //cout << "tmp = " << tmp << " fSkipBinsUnfolded = " << fSkipBinsUnfolded << " fPtMinUnfolded = " << fPtMinUnfolded << endl; |
31d3e4f0 | 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); | |
343 | } | |
adf3920d | 344 | else { |
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 | |
350 | } | |
351 | ||
352 | printf(" nbins[fDimGen] = %d nbins[fDimRec] = %d\n",nbins[fDimGen],nbins[fDimRec]); | |
31d3e4f0 | 353 | |
354 | Double_t binWidth[2]; | |
355 | binWidth[fDimRec] = (double)((fPtMax-fPtMin)/nbins[fDimRec]); | |
356 | binWidth[fDimGen] = (double)((fPtMaxUnfolded-fPtMinUnfolded)/nbins[fDimGen]); | |
357 | ||
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; | |
364 | ||
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]); | |
369 | ||
370 | Double_t binArrayPt0[nbins[fDimRec]+1]; | |
371 | Double_t binArrayPt1[nbins[fDimGen]+1]; | |
372 | Double_t xmaxGen = TMath::Max(xmax[fDimGen],fPtMaxUnfoldedHigh); | |
373 | ||
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]; | |
377 | } | |
378 | binArrayPt0[nbins[fDimRec]]= xmax[fDimRec]; | |
379 | ||
380 | //Define bin limits generated/unfolded axis | |
fa7c34ba | 381 | binArrayPt1[0] = fPtMinUnfolded; |
382 | for(int i=1; i<nbins[fDimGen]; i++) { | |
31d3e4f0 | 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; | |
387 | } | |
fa7c34ba | 388 | else |
389 | binArrayPt1[i] = xmin[fDimGen]+(double)i*binWidth[fDimGen]; | |
31d3e4f0 | 390 | //printf("RM. i = %d \t binArrayPt[i] = %f \n",i,binArrayPt1[i]); |
391 | } | |
392 | binArrayPt1[nbins[fDimGen]]= xmaxGen; | |
393 | ||
394 | ||
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 | |
ef62323a | 398 | if(fResponseMatrix) delete fResponseMatrix; |
31d3e4f0 | 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}"); | |
403 | ||
404 | fResponseMatrix->SetBinEdges(fDimRec,binArrayPt0); | |
405 | fResponseMatrix->SetBinEdges(fDimGen,binArrayPt1); | |
406 | ||
ef62323a | 407 | Int_t bin[2] = {1,1}; |
408 | for(int i=1; i<fResponseMatrix->GetAxis(0)->GetNbins(); i++) { | |
409 | bin[0]=i; | |
410 | for(int j=1; j<fResponseMatrix->GetAxis(1)->GetNbins(); j++) { | |
411 | bin[1]=j; | |
412 | fResponseMatrix->SetBinContent(bin,0.); | |
413 | } | |
414 | } | |
415 | ||
416 | ||
31d3e4f0 | 417 | |
418 | } | |
419 | ||
420 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
421 | void AliAnaChargedJetResponseMaker::InitializeEfficiency() { | |
422 | // | |
423 | // Define dimensions of efficiency THnSparse | |
424 | // | |
425 | ||
426 | if(!fResponseMatrix) { | |
427 | printf("AliAnaChargedJetResponseMaker::InitializeEfficiency()\n"); | |
428 | printf("Can not define dimensions efficiency without fResponseMatrix dimensions. Aborting \n"); | |
429 | return; | |
430 | } | |
431 | ||
432 | TAxis *genAxis = fResponseMatrix->GetAxis(fDimGen); | |
433 | ||
434 | Int_t nbinsEff[fDimensions]; | |
435 | Double_t xminEff[fDimensions]; | |
436 | Double_t xmaxEff[fDimensions]; | |
437 | ||
438 | for(int dim = 0; dim<fDimensions; dim++) { | |
439 | nbinsEff[dim] = genAxis->GetNbins(); | |
440 | xminEff[dim] = genAxis->GetXmin(); | |
441 | xmaxEff[dim] = genAxis->GetXmax(); | |
442 | } | |
443 | ||
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}"); | |
448 | ||
449 | const Double_t *binArrayPt = genAxis->GetXbins()->GetArray(); | |
450 | fEfficiency->SetBinEdges(0,binArrayPt); | |
451 | ||
452 | } | |
453 | ||
454 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
455 | void AliAnaChargedJetResponseMaker::InitializeResponseMatrixFine() { | |
456 | // | |
457 | //Initialize fine response matrix | |
458 | // | |
459 | ||
460 | Int_t nbinsFine[fDimensions*2]; | |
fa7c34ba | 461 | Double_t xminFine[fDimensions*2]; |
462 | Double_t xmaxFine[fDimensions*2]; | |
0560ee0f | 463 | // Double_t pTarrayFine[fDimensions*2]; |
31d3e4f0 | 464 | |
fa7c34ba | 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.); | |
0bf6381d | 469 | xmaxFine[fDimRec] = fResponseMatrix->GetAxis(fDimGen)->GetXmax();//+40.; |
470 | xmaxFine[fDimGen] = fResponseMatrix->GetAxis(fDimGen)->GetXmax();//+40.; | |
0560ee0f | 471 | // pTarrayFine[fDimRec] = 0.; |
472 | // pTarrayFine[fDimGen] = 0.; | |
31d3e4f0 | 473 | |
474 | Double_t binWidth[2]; | |
475 | binWidth[fDimRec] = fResponseMatrix->GetAxis(fDimRec)->GetBinWidth(1); | |
476 | ||
477 | Double_t binWidthFine[2]; | |
478 | binWidthFine[fDimRec] = binWidth[fDimRec]/((double)fFineFrac); | |
479 | binWidthFine[fDimGen] = binWidthFine[fDimRec]*(double)fBinWidthFactorUnfolded; | |
0bf6381d | 480 | // binWidthFine[fDimRec] = binWidthFine[fDimGen]; |
31d3e4f0 | 481 | |
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 | |
484 | ||
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]); | |
489 | ||
490 | ||
491 | Double_t binArrayPt0Fine[nbinsFine[fDimRec]+1]; | |
492 | Double_t binArrayPt1Fine[nbinsFine[fDimGen]+1]; | |
493 | ||
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]); | |
497 | } | |
498 | binArrayPt0Fine[nbinsFine[fDimRec]]= xmaxFine[fDimRec]; | |
499 | ||
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]); | |
503 | } | |
504 | binArrayPt1Fine[nbinsFine[fDimGen]]= xmaxFine[fDimGen]; | |
505 | ||
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 | |
ef62323a | 509 | if(fResponseMatrixFine) delete fResponseMatrixFine; |
31d3e4f0 | 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}"); | |
514 | ||
515 | fResponseMatrixFine->SetBinEdges(fDimRec,binArrayPt0Fine); | |
516 | fResponseMatrixFine->SetBinEdges(fDimGen,binArrayPt1Fine); | |
517 | ||
ef62323a | 518 | Int_t bin[2] = {1,1}; |
519 | for(int i=1; i<fResponseMatrixFine->GetAxis(0)->GetNbins(); i++) { | |
520 | bin[0]=i; | |
521 | for(int j=1; j<fResponseMatrixFine->GetAxis(1)->GetNbins(); j++) { | |
522 | bin[1]=j; | |
523 | fResponseMatrixFine->SetBinContent(bin,0.); | |
524 | } | |
525 | } | |
526 | ||
31d3e4f0 | 527 | } |
528 | ||
529 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
530 | void AliAnaChargedJetResponseMaker::InitializeEfficiencyFine() { | |
531 | // | |
532 | // Define dimensions of efficiency THnSparse | |
533 | // | |
534 | ||
535 | if(!fResponseMatrixFine) { | |
536 | printf("AliAnaChargedJetResponseMaker::InitializeEfficiencyFine()\n"); | |
537 | printf("Can not define dimensions efficiency without fResponseMatrixFine dimensions. Aborting \n"); | |
538 | return; | |
539 | } | |
540 | ||
541 | TAxis *genAxis = fResponseMatrixFine->GetAxis(fDimGen); | |
542 | ||
543 | Int_t nbinsEff[fDimensions]; | |
544 | Double_t xminEff[fDimensions]; | |
545 | Double_t xmaxEff[fDimensions]; | |
546 | ||
547 | for(int dim = 0; dim<fDimensions; dim++) { | |
548 | nbinsEff[dim] = genAxis->GetNbins(); | |
549 | xminEff[dim] = genAxis->GetXmin(); | |
550 | xmaxEff[dim] = genAxis->GetXmax(); | |
551 | } | |
552 | ||
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}"); | |
557 | ||
558 | const Double_t *binArrayPt = genAxis->GetXbins()->GetArray(); | |
559 | fEfficiencyFine->SetBinEdges(0,binArrayPt); | |
560 | ||
561 | } | |
562 | ||
563 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
564 | void AliAnaChargedJetResponseMaker::FillResponseMatrixFineAndMerge() { | |
565 | // | |
566 | // Fill fine response matrix | |
567 | // | |
568 | ||
ef62323a | 569 | if(!fResponseMatrix) { |
570 | printf("Dimensions of fResponseMatrix have to be defined first. Aborting!"); | |
571 | return; | |
572 | } | |
573 | ||
574 | if(!fResponseMatrixFine) { | |
575 | printf("Dimensions of fResponseMatrixFine have to be defined first. Aborting!"); | |
576 | return; | |
577 | } | |
578 | ||
31d3e4f0 | 579 | TAxis *genAxis = fResponseMatrixFine->GetAxis(fDimGen); |
580 | TAxis *recAxis = fResponseMatrixFine->GetAxis(fDimRec); | |
581 | ||
582 | Int_t nbinsFine[fDimensions*2]; | |
0560ee0f | 583 | //Double_t xminFine[fDimensions*2]; |
584 | //Double_t xmaxFine[fDimensions*2]; | |
31d3e4f0 | 585 | Double_t pTarrayFine[fDimensions*2]; |
586 | ||
587 | nbinsFine[fDimGen] = genAxis->GetNbins(); | |
588 | nbinsFine[fDimRec] = recAxis->GetNbins(); | |
0560ee0f | 589 | // xminFine[fDimGen] = genAxis->GetXmin(); |
590 | // xminFine[fDimRec] = recAxis->GetXmin(); | |
591 | // xmaxFine[fDimGen] = genAxis->GetXmax(); | |
592 | // xmaxFine[fDimRec] = recAxis->GetXmax(); | |
fa7c34ba | 593 | pTarrayFine[fDimGen] = 0.; |
594 | pTarrayFine[fDimRec] = 0.; | |
31d3e4f0 | 595 | |
596 | double sumyield = 0.; | |
597 | double sumyielderror2 = 0.; | |
598 | ||
ada4df5e | 599 | int ipt[2] = {0,0}; |
600 | int iptMerged[2] = {0,0}; | |
601 | int ierror[2] = {0,0}; | |
31d3e4f0 | 602 | |
603 | Int_t check = 0; | |
604 | Double_t pTgen, pTrec; | |
605 | Double_t yield = 0.; | |
606 | Double_t error = 0.; | |
607 | ||
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.; | |
614 | } | |
615 | } | |
616 | ||
617 | for(int iy=1; iy<=nbinsFine[fDimGen]; iy++) { //gen | |
618 | pTgen = fResponseMatrixFine->GetAxis(fDimGen)->GetBinCenter(iy); | |
619 | pTarrayFine[fDimGen] = pTgen; | |
620 | ierror[fDimGen]=iy; | |
621 | sumyield = 0.; | |
622 | check = 0; | |
623 | ||
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); | |
629 | error = 0.; | |
630 | } | |
631 | else if(fResolutionType==kResiduals) { | |
632 | yield = fhDeltaPt->Interpolate(pTrec-pTgen); | |
633 | error = 0.; | |
634 | } | |
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); | |
640 | } | |
641 | yield=yield*width; | |
642 | error=error*width; | |
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]; | |
646 | ||
647 | sumyield+=yield; | |
648 | sumyielderror2 += error*error; | |
649 | ||
650 | pTarrayFine[fDimRec] = pTrec; | |
651 | ierror[fDimRec] = ix; | |
652 | fResponseMatrixFine->Fill(pTarrayFine,yield); | |
653 | if(fbCalcErrors) fResponseMatrixFine->SetBinError(ierror,error); | |
654 | ||
655 | }// ix (dimRec) loop | |
656 | ||
657 | //Normalize to total number of counts =1 | |
658 | if(sumyield>1) { | |
659 | ipt[fDimGen]=iy; | |
660 | for(int ix=1; ix<=nbinsFine[fDimRec]; ix++) { | |
661 | ipt[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)); | |
668 | } | |
669 | ||
670 | } | |
671 | } | |
672 | ||
673 | int bin[1]; | |
674 | bin[0] = iy; | |
675 | fEfficiencyFine->SetBinContent(bin,sumyield); | |
676 | if(fbCalcErrors) fEfficiencyFine->SetBinError(bin,TMath::Sqrt(sumyielderror2)); | |
677 | ||
678 | //fill merged response matrix | |
679 | ||
0bf6381d | 680 | //find bin in fine RM corresponding to mimimum/maximum bin of merged RM on rec axis |
31d3e4f0 | 681 | int ixMin = fResponseMatrixFine->GetAxis(fDimRec)->FindBin(fResponseMatrix->GetAxis(fDimRec)->GetXmin()); |
682 | int ixMax = fResponseMatrixFine->GetAxis(fDimRec)->FindBin(fResponseMatrix->GetAxis(fDimRec)->GetXmax()); | |
683 | ||
684 | if(fResponseMatrixFine->GetAxis(fDimGen)->GetBinLowEdge(iy) >= fResponseMatrix->GetAxis(fDimGen)->GetXmin()) { | |
685 | ipt[fDimGen]=iy; | |
686 | iptMerged[fDimGen]=fResponseMatrix->GetAxis(fDimGen)->FindBin(pTgen); | |
687 | ||
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); | |
694 | if(powInteg>0.) | |
695 | weight = f1MergeFunction->Integral(fResponseMatrixFine->GetAxis(fDimGen)->GetBinLowEdge(iy),fResponseMatrixFine->GetAxis(fDimGen)->GetBinUpEdge(iy))/powInteg; | |
696 | // printf("weight: %f \n", weight ); | |
697 | } else { | |
698 | weight = 1./((double)fFineFrac); | |
adf3920d | 699 | if(fbVariableBinning && pTgen<fPtMaxUnfVarBinning) weight=1./(0.5*(double)fFineFrac); |
31d3e4f0 | 700 | } |
701 | ||
702 | for(int ix=ixMin; ix<=ixMax; ix++) { //loop reconstructed axis | |
703 | pTrec = fResponseMatrixFine->GetAxis(fDimRec)->GetBinCenter(ix); | |
704 | ipt[fDimRec]=ix; | |
705 | iptMerged[fDimRec]=fResponseMatrix->GetAxis(fDimRec)->FindBin(pTrec); | |
706 | ||
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; | |
709 | } | |
710 | ||
711 | }//loop gen axis fine response matrix | |
712 | ||
713 | } // iy (dimGen) loop | |
714 | ||
715 | //Fill Efficiency corresponding to merged response matrix | |
31d3e4f0 | 716 | for(int iy=1; iy<=fResponseMatrix->GetAxis(fDimGen)->GetNbins(); iy++) { //gen |
717 | sumyield = 0.; | |
718 | ipt[fDimGen] = iy; | |
719 | ||
720 | for(int ix=1; ix<=fResponseMatrix->GetAxis(fDimRec)->GetNbins(); ix++) { //rec | |
721 | ipt[fDimRec] = ix; | |
722 | sumyield += fResponseMatrix->GetBinContent(ipt); | |
723 | ||
724 | if(fbCalcErrors) fResponseMatrix->SetBinError(ipt,TMath::Sqrt(errorArray[iy][ix])); | |
725 | } | |
726 | printf("RM: pTgen: %f \t sumyield: %f \n",fResponseMatrix->GetAxis(fDimGen)->GetBinCenter(iy),sumyield); | |
727 | int bin[1]; | |
728 | bin[0] = iy; | |
729 | fEfficiency->SetBinContent(bin,sumyield); | |
730 | if(fbCalcErrors) fEfficiency->SetBinError(bin,0); | |
731 | } | |
732 | ||
733 | if(fDebug) printf("fResponseMatrixFine->GetNbins(): %d \n",(int)fResponseMatrixFine->GetNbins()); | |
734 | if(fDebug) printf("fResponseMatrix->GetNbins(): %d \n",(int)fResponseMatrix->GetNbins()); | |
735 | ||
736 | if(fDebug) printf("Done constructing response matrix\n"); | |
737 | ||
738 | } | |
739 | ||
ef62323a | 740 | //-------------------------------------------------------------------------------------------------------------------------------------------------- |
03372fd1 | 741 | TH2* AliAnaChargedJetResponseMaker::MakeResponseMatrixRebin(TH2 *hRMFine, TH2 *hRM, Bool_t useFunctionWeight) { |
ef62323a | 742 | |
743 | // | |
744 | // Rebin matrix hRMFine to dimensions of hRM | |
03372fd1 | 745 | // function returns matrix in TH2D format (hRM2) with dimensions from hRM |
ef62323a | 746 | // |
747 | ||
748 | TH2 *hRM2 = (TH2*)hRM->Clone("hRM2"); | |
749 | hRM2->Reset(); | |
750 | ||
03372fd1 | 751 | if(useFunctionWeight) cout << "Use function to do weighting" << endl; |
752 | ||
275d3481 | 753 | //First normalize columns of input |
ef62323a | 754 | const Int_t nbinsNorm = hRM2->GetNbinsX(); |
0bf6381d | 755 | if(fDebug) cout << "nbinsNorm: " << nbinsNorm << endl; |
ef62323a | 756 | |
757 | TArrayF *normVector = new TArrayF(nbinsNorm); | |
758 | ||
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; | |
03372fd1 | 766 | if(useFunctionWeight) |
767 | normVector->SetAt(f1MergeFunction->Integral(xLow,xUp),ix-1); | |
768 | else | |
769 | normVector->SetAt(hRMFine->Integral(binxLowFine,binxUpFine,1,hRMFine->GetYaxis()->GetNbins()),ix-1); | |
0bf6381d | 770 | //if(fDebug) cout << "ix norm: " << normVector->At(ix-1) << endl; |
ef62323a | 771 | } |
772 | ||
773 | Double_t content, oldcontent = 0.; | |
774 | Int_t ixNew = 0; | |
775 | Int_t iyNew = 0; | |
275d3481 | 776 | Double_t xvalLo, xvalUp, yvalLo, yvalUp; |
ef62323a | 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++) { | |
275d3481 | 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)); | |
ef62323a | 786 | for(int iy=1; iy<=hRMFine->GetYaxis()->GetNbins(); iy++) { |
275d3481 | 787 | yvalLo = hRMFine->GetYaxis()->GetBinLowEdge(iy); |
788 | yvalUp = hRMFine->GetYaxis()->GetBinUpEdge(iy); | |
789 | if(yvalLo<ymin || yvalUp>ymax) continue; | |
ef62323a | 790 | content = hRMFine->GetBinContent(ix,iy); |
275d3481 | 791 | iyNew = hRM2->GetYaxis()->FindBin(hRMFine->GetYaxis()->GetBinCenter(iy)); |
ef62323a | 792 | oldcontent = hRM2->GetBinContent(ixNew,iyNew); |
793 | ||
275d3481 | 794 | //if(fDebug) cout << "ixNew: " << ixNew << " " << xvalLo << " iyNew: " << iyNew << " " << yvalLo << " content: " << content << " oldContent: " << oldcontent << " newContent: " << oldcontent+content << endl; |
ef62323a | 795 | Double_t weight = 1.; |
03372fd1 | 796 | if(normVector->At(ixNew-1)>0.) { |
797 | if(useFunctionWeight) | |
798 | weight = f1MergeFunction->Integral(xvalLo,xvalUp)/normVector->At(ixNew-1); | |
799 | else | |
800 | weight = 1./normVector->At(ixNew-1); | |
801 | } | |
ef62323a | 802 | hRM2->SetBinContent(ixNew,iyNew,oldcontent+content*weight); |
803 | } | |
804 | } | |
805 | ||
806 | if(normVector) delete normVector; | |
807 | ||
808 | return hRM2; | |
809 | ||
810 | } | |
811 | ||
03372fd1 | 812 | //-------------------------------------------------------------------------------------------------------------------------------------------------- |
813 | TH2* AliAnaChargedJetResponseMaker::CreateTruncated2DHisto(TH2 *h2, Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax) { | |
814 | ||
815 | // | |
816 | // Limit axis range of 2D histogram | |
817 | // | |
818 | ||
819 | Int_t binMinXh2 = h2->GetXaxis()->FindBin(xmin); | |
820 | if(h2->GetXaxis()->GetBinLowEdge(binMinXh2) < xmin ) binMinXh2++; | |
821 | if(h2->GetXaxis()->GetBinLowEdge(binMinXh2) > xmin ) binMinXh2--; | |
822 | ||
823 | Int_t binMinYh2 = h2->GetYaxis()->FindBin(ymin); | |
824 | if(h2->GetYaxis()->GetBinLowEdge(binMinYh2) < ymin ) binMinYh2++; | |
825 | if(h2->GetYaxis()->GetBinLowEdge(binMinYh2) > ymin ) binMinYh2--; | |
826 | ||
827 | Int_t binMaxXh2 = h2->GetXaxis()->FindBin(xmax); | |
828 | if(h2->GetXaxis()->GetBinUpEdge(binMaxXh2) < xmax ) binMaxXh2++; | |
829 | if(h2->GetXaxis()->GetBinUpEdge(binMaxXh2) > xmax ) binMaxXh2--; | |
830 | ||
831 | Int_t binMaxYh2 = h2->GetYaxis()->FindBin(ymax); | |
832 | if(h2->GetYaxis()->GetBinUpEdge(binMaxYh2) < ymax ) binMaxYh2++; | |
833 | if(h2->GetYaxis()->GetBinUpEdge(binMaxYh2) > ymax ) binMaxYh2--; | |
834 | ||
0bf6381d | 835 | Int_t nbinsX = binMaxXh2-binMinXh2+1; |
836 | Int_t nbinsY = binMaxYh2-binMinYh2+1; | |
03372fd1 | 837 | |
838 | Double_t *binsX = new Double_t[nbinsX+1]; | |
839 | Double_t *binsY = new Double_t[nbinsY+1]; | |
840 | ||
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); | |
844 | ||
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); | |
848 | ||
849 | TH2 *h2Lim = new TH2D("h2Lim","h2Lim",nbinsX,binsX,nbinsY,binsY); | |
850 | ||
851 | for(int ix=1; ix<=nbinsX; ix++) { | |
852 | // cout << "ix: " << ix << " " << binsX[ix] << endl; | |
853 | for(int iy=1; iy<=nbinsY; iy++) { | |
0bf6381d | 854 | // cout << "ix: " << ix << " " << binsX[ix] << "\tiy: " << iy << " " << binsY[iy] << endl; |
03372fd1 | 855 | |
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); | |
0bf6381d | 859 | if(fbCalcErrors) h2Lim->SetBinError(ix,iy,error); |
03372fd1 | 860 | |
861 | } | |
862 | } | |
863 | ||
864 | ||
865 | ||
866 | return h2Lim; | |
867 | } | |
868 | ||
869 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
870 | TH2* AliAnaChargedJetResponseMaker::TruncateAxisRangeResponseMatrix(TH2 *hRMOrig, Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax) { | |
871 | ||
872 | // | |
873 | // Limit axis range of response matrix without changing normalization | |
874 | // | |
875 | ||
876 | //TH2 *hRMLimit | |
877 | //TH2 *hRMLimit2 = (TH2*)hRMLimit->Clone("hRMLimit2"); | |
878 | ||
879 | TH2 *hRMLimit2 = CreateTruncated2DHisto(hRMOrig, xmin, xmax, ymin, ymax); | |
880 | hRMLimit2->Reset(); | |
881 | ||
882 | double binCent[2] = {0.,0.}; | |
883 | double content = 0.; | |
884 | double error = 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]); | |
892 | ||
893 | content = hRMOrig->GetBinContent(binOrig[0],binOrig[1]); | |
894 | error = hRMOrig->GetBinError(binOrig[0],binOrig[1]); | |
895 | ||
896 | hRMLimit2->SetBinContent(ix,iy,content); | |
0bf6381d | 897 | if(fbCalcErrors) hRMLimit2->SetBinError(ix,iy,error); |
03372fd1 | 898 | |
899 | } | |
900 | } | |
901 | ||
902 | ||
903 | return hRMLimit2; | |
904 | ||
905 | } | |
906 | ||
0bf6381d | 907 | //-------------------------------------------------------------------------------------------------------------------------------------------------- |
908 | Bool_t AliAnaChargedJetResponseMaker::CheckInputForCombinedResponse() { | |
909 | ||
910 | Bool_t check = kTRUE; | |
911 | ||
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; | |
917 | ||
918 | return check; | |
919 | } | |
920 | ||
921 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
922 | void AliAnaChargedJetResponseMaker::MakeResponseMatrixCombined(Int_t skipBins, Int_t binWidthFactor, Int_t extraBins, Bool_t bVariableBinning, Double_t ptmin) { | |
923 | ||
924 | //Check if all input histos are there otherwise break | |
925 | if(!CheckInputForCombinedResponse()) { | |
926 | printf("Not all input available..... \n"); | |
927 | return; | |
928 | } | |
929 | ||
930 | // Make response bkg fluctuations | |
931 | MakeResponseMatrixJetsFineMerged(skipBins,binWidthFactor,extraBins,bVariableBinning,ptmin); | |
932 | ||
933 | //Get TH2D with dimensions we want final RM | |
934 | TH2D *h2ResponseMatrix = fResponseMatrix->Projection(0,1,"E"); | |
935 | h2ResponseMatrix->Reset(); | |
936 | ||
937 | //Get fine binned response matrix from bkg fluctuations | |
938 | TH2D *h2ResponseMatrixFine = fResponseMatrixFine->Projection(0,1,"E"); | |
939 | ||
940 | // Rebin response detector effects | |
941 | const TArrayD *arrayX = h2ResponseMatrixFine->GetXaxis()->GetXbins(); | |
942 | const Double_t *xbinsArray = arrayX->GetArray(); | |
943 | ||
944 | TH2D *h2ResponseDetTmp = new TH2D("h2ResponseDetTmp","h2ResponseDetTmp",h2ResponseMatrixFine->GetNbinsX(),xbinsArray,h2ResponseMatrixFine->GetNbinsX(),xbinsArray); | |
945 | ||
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"); | |
951 | ||
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"); | |
956 | ||
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})"); | |
963 | ||
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"); | |
967 | ||
968 | fhEfficiencyCombined = (TH1D*)fh2ResponseMatrixCombined->ProjectionX("fhEfficiencyCombined"); | |
969 | fhEfficiencyCombined->Reset(); | |
970 | ||
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); | |
975 | } | |
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.; | |
980 | ||
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); | |
984 | kCounter+=1.; | |
985 | } | |
986 | } | |
987 | kJetEffDet = kJetEffDet/kCounter; | |
988 | ||
989 | if(kJetEffDet==0.) kJetEffDet=1.; | |
990 | ||
991 | fhEfficiencyCombined->SetBinContent(i,sumyield*kJetEffDet); | |
992 | ||
993 | } | |
994 | ||
995 | } | |
996 | ||
275d3481 | 997 | //-------------------------------------------------------------------------------------------------------------------------------------------------- |
998 | TH2* AliAnaChargedJetResponseMaker::MultiplityResponseMatrices(TH2 *h2RMDeltaPt, TH2 *h2RMDetector) { | |
999 | ||
1000 | // | |
1001 | // Make combined response matrix (background flucutuations + detector effects) | |
1002 | // | |
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 | |
1005 | // | |
1006 | // Function assumes that generated/unfolded axis is x-axis and reconstructed is on y-axis on both respone matrices | |
1007 | ||
1008 | ||
03372fd1 | 1009 | TH2D *h2ResponseMatrixCombined = (TH2D*)h2RMDeltaPt->Clone("h2ResponseMatrixCombined"); //h2RMDeltaPt is the bkg fluctuations RM which has the dimensions we want for the combined response matrix |
275d3481 | 1010 | h2ResponseMatrixCombined->SetTitle("h2ResponseMatrixCombined"); |
1011 | h2ResponseMatrixCombined->SetName("h2ResponseMatrixCombined"); | |
1012 | ||
1013 | // M = RM_deltaPt * RM_DetEffects * T (M=measured T=truth) | |
1014 | // Dimensions: | |
1015 | // mx1 = mxd * dxt * tx1 | |
1016 | // m = measured bins | |
1017 | // t = truth 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) | |
1021 | ||
1022 | if(fDebug) { | |
03372fd1 | 1023 | printf("Nt=%d\n",h2ResponseMatrixCombined->GetNbinsX()); |
1024 | printf("Nm=%d\n",h2ResponseMatrixCombined->GetNbinsY()); | |
1025 | printf("Nd=%d\n",h2RMDeltaPt->GetNbinsX()); | |
275d3481 | 1026 | } |
1027 | ||
03372fd1 | 1028 | |
275d3481 | 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); | |
03372fd1 | 1034 | // if(t==10 && m==10) cout << "sum m,d=" << m << "," << d << endl; |
275d3481 | 1035 | }//d-loop |
03372fd1 | 1036 | // if(t==10) cout << "t,m = " << t << "," << m << "\tvalueSum: " << valueSum << endl; |
275d3481 | 1037 | h2ResponseMatrixCombined->SetBinContent(t,m,valueSum); |
1038 | } //m-loop | |
1039 | }//t-loop | |
1040 | ||
1041 | return h2ResponseMatrixCombined; | |
1042 | ||
1043 | } | |
31d3e4f0 | 1044 | |
5d87a047 | 1045 | //-------------------------------------------------------------------------------------------------------------------------------------------------- |
1046 | TH2* AliAnaChargedJetResponseMaker::GetTransposeResponsMatrix(TH2 *h2RM) { | |
1047 | // | |
1048 | // Transpose response matrix | |
1049 | // | |
1050 | ||
e50460de | 1051 | Printf("AliAnaChargedJetResponseMaker::GetTransposeResponsMatrix"); |
1052 | ||
5d87a047 | 1053 | //Initialize transposed response matrix h2RMTranspose |
e50460de | 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(); | |
1057 | ||
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(); | |
1061 | ||
1062 | ||
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)); | |
1067 | } | |
1068 | ybinsArray[h2RM->GetNbinsY()] = h2RM->GetYaxis()->GetBinUpEdge(h2RM->GetNbinsY()); | |
1069 | ||
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()); | |
5d87a047 | 1074 | |
5d87a047 | 1075 | |
1076 | TH2D *h2RMTranspose = new TH2D("h2RMTranspose","h2RMTranspose",h2RM->GetNbinsY(),ybinsArray,h2RM->GetNbinsX(),xbinsArray); | |
1077 | ||
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)); | |
1082 | } | |
1083 | } | |
1084 | ||
1085 | ||
1086 | return h2RMTranspose; | |
1087 | ||
1088 | } | |
1089 | ||
1090 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
1091 | TH2* AliAnaChargedJetResponseMaker::NormalizeResponsMatrixYaxisWithPrior(TH2 *h2RM, TH1 *hPrior) { | |
1092 | // | |
1093 | // Normalize such that the Y projection is the prior | |
1094 | // | |
1095 | ||
03372fd1 | 1096 | // TH1D *hProjRespY = (TH1D*)h2RM->ProjectionY("hProjRespY"); |
5d87a047 | 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); | |
1104 | } | |
1105 | } | |
1106 | ||
1107 | return h2RM; | |
1108 | } | |
1109 | ||
31d3e4f0 | 1110 | //-------------------------------------------------------------------------------------------------------------------------------------------------- |
1111 | TH1D* AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TH1 *hGen, TH2 *hResponse,TH1 *hEfficiency,Bool_t bDrawSlices) { | |
1112 | ||
1113 | // | |
1114 | // Multiply hGen with response matrix to obtain refolded spectrum | |
ef62323a | 1115 | // Efficiency must be given. In case efficiency is 1 use SetFlatEfficiency(1.) before calling function |
31d3e4f0 | 1116 | // |
1117 | ||
ef62323a | 1118 | if(!hEfficiency) { |
5d87a047 | 1119 | // printf("Efficiency must be given. In case efficiency is 1 use SetFlatEfficiency(1.) before calling function. Aborting!"); |
1120 | // return 0; | |
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.); | |
ef62323a | 1125 | } |
1126 | ||
31d3e4f0 | 1127 | //For response |
1128 | //x-axis: generated | |
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; | |
1131 | ||
1132 | TH1D *hRec = hResponse->ProjectionY("hRec"); | |
1133 | hRec->Sumw2(); | |
1134 | hRec->Reset(); | |
1135 | hRec->SetTitle("hRec"); | |
1136 | hRec->SetName("hRec"); | |
1137 | ||
1138 | for(int irec=1; irec<=hRec->GetNbinsX(); irec++) | |
1139 | hRec->SetBinContent(irec,0); | |
1140 | ||
1141 | TH1D *hRecGenBin = 0x0; | |
1142 | TCanvas *cSlices = 0x0; | |
1143 | if(bDrawSlices) { | |
1144 | cSlices = new TCanvas("cSlices","cSlices:Slices",400,400); | |
1145 | gPad->SetLogy(); | |
1146 | } | |
1147 | ||
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]; | |
0bf6381d | 1153 | for(int irec=0; irec<=hRec->GetNbinsX(); irec++) |
1154 | sumError2[irec]=0.; | |
31d3e4f0 | 1155 | Double_t eff = 0.; |
1156 | ||
1157 | for(int igen=1; igen<=hGen->GetNbinsX(); igen++) { | |
1158 | //get pTMC | |
1159 | sumYield = 0.; | |
1160 | if(fEffFlat>0.) | |
1161 | eff = fEffFlat; | |
1162 | else | |
1163 | eff = hEfficiency->GetBinContent(igen); | |
1164 | yieldMC = hGen->GetBinContent(igen)*eff; | |
adf3920d | 1165 | yieldMCerror = hGen->GetBinError(igen)*eff; |
31d3e4f0 | 1166 | |
1167 | if(bDrawSlices) { | |
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)); | |
1173 | } | |
1174 | ||
1175 | for(int irec=1; irec<=hRec->GetNbinsX(); irec++) { | |
1176 | hRec->AddBinContent(irec,yieldMC*hResponse->GetBinContent(igen,irec)); | |
1177 | sumYield+=hResponse->GetBinContent(igen,irec); | |
346cfe4a | 1178 | // Double_t A = yieldMC*hResponse->GetBinError(igen,irec); |
31d3e4f0 | 1179 | Double_t B = hResponse->GetBinContent(igen,irec)*yieldMCerror; |
346cfe4a | 1180 | // Double_t tmp2 = A*A + B*B; |
adf3920d | 1181 | //sumError2[irec-1] += tmp2 ; |
1182 | sumError2[irec-1] += B*B; | |
31d3e4f0 | 1183 | |
1184 | if(bDrawSlices) | |
1185 | hRecGenBin->SetBinContent(irec,yieldMC*hResponse->GetBinContent(igen,irec)); | |
1186 | ||
1187 | } | |
1188 | if(bDrawSlices) { | |
1189 | cSlices->cd(); | |
1190 | hRecGenBin->SetLineColor(igen); | |
1191 | if(igen==1) hRecGenBin->DrawCopy(); | |
1192 | else hRecGenBin->DrawCopy("same"); | |
1193 | } | |
1194 | ||
1195 | if(hRecGenBin) delete hRecGenBin; | |
1196 | ||
1197 | cout << "igen: " << igen << "\tpTMC: " << hGen->GetXaxis()->GetBinCenter(igen) << "\teff:" << eff << "\tsumYield: " << sumYield << endl; | |
1198 | } | |
1199 | ||
ef62323a | 1200 | for(int i=0; i<=nbinsRec; i++) { |
1201 | if(sumError2[i]>0.) | |
1202 | hRec->SetBinError(i+1,TMath::Sqrt(sumError2[i])); | |
1203 | } | |
31d3e4f0 | 1204 | |
1205 | ||
1206 | return hRec; | |
1207 | } | |
1208 | ||
1209 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
1210 | TH1D* AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TF1 *fGen, TH2 *hResponse,TH1 *hEfficiency) { | |
1211 | ||
ef62323a | 1212 | // |
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 | |
1215 | // | |
1216 | ||
31d3e4f0 | 1217 | //For response |
1218 | //x-axis: generated | |
1219 | //y-axis: reconstructed | |
1220 | ||
ada4df5e | 1221 | if(fDebug>0) printf(">>AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TF1 *fGen, TH2 *hResponse,TH1 *hEfficiency)\n"); |
ef62323a | 1222 | |
31d3e4f0 | 1223 | TH1D *hRec = hResponse->ProjectionY("hRec"); |
1224 | hRec->Sumw2(); | |
1225 | hRec->Reset(); | |
1226 | hRec->SetTitle("hRec"); | |
1227 | hRec->SetName("hRec"); | |
1228 | ||
1229 | // TH1D *hRec = new TH1D("hRec","hRec",hResponse->GetNbinsY(),hResponse->GetYaxis()->GetXmin(),hResponse->GetYaxis()->GetXmax()); | |
1230 | ||
1231 | for(int irec=1; irec<=hRec->GetNbinsX(); irec++) | |
1232 | hRec->SetBinContent(irec,0); | |
1233 | ||
1234 | Double_t yieldMC = 0.; | |
1235 | Double_t sumYield = 0.; | |
1236 | Double_t eff = 0.; | |
1237 | for(int igen=1; igen<=hResponse->GetNbinsX(); igen++) { | |
1238 | //get pTMC | |
1239 | sumYield = 0.; | |
1240 | double pTMC = hResponse->GetXaxis()->GetBinCenter(igen); | |
ada4df5e | 1241 | if(hEfficiency) { |
1242 | int binEff = hEfficiency->FindBin(pTMC); | |
31d3e4f0 | 1243 | eff = hEfficiency->GetBinContent(binEff); |
ada4df5e | 1244 | } |
1245 | else eff = 1.; | |
1246 | // yieldMC = fGen->Eval(pTMC)*eff; | |
1247 | yieldMC = fGen->Integral(hResponse->GetXaxis()->GetBinLowEdge(igen),hResponse->GetXaxis()->GetBinUpEdge(igen))*eff; | |
31d3e4f0 | 1248 | for(int irec=1; irec<=hResponse->GetNbinsY(); irec++) { |
1249 | hRec->AddBinContent(irec,yieldMC*hResponse->GetBinContent(igen,irec)); | |
1250 | sumYield+=hResponse->GetBinContent(igen,irec); | |
1251 | } | |
ada4df5e | 1252 | // cout << "igen: " << igen << "\tpTMC: " << pTMC << "\tsumYield: " << sumYield << endl; |
1253 | // cout << "yieldMC: " << yieldMC << endl; | |
1254 | ||
31d3e4f0 | 1255 | } |
1256 | ||
1257 | return hRec; | |
1258 | } | |
1259 | ||
0bf6381d | 1260 | //-------------------------------------------------------------------------------------------------------------------------------------------------- |
1261 | /* static */ Double_t AliAnaChargedJetResponseMaker::GetBetaPerDOFValue(Int_t betaColl, Int_t betaOpt) { | |
1262 | ||
1263 | Float_t betaPerDOFoptions[6] = {0.}; | |
1264 | //define beta per degree of freedom (DOF=nbins in measured) | |
1265 | if(betaColl==0) { | |
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; | |
1272 | } | |
1273 | if(betaColl==1) { | |
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; | |
1280 | } | |
1281 | if(betaColl==2) { | |
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; | |
1288 | } | |
1289 | if(betaColl==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; | |
1296 | } | |
1297 | if(betaColl==4) { | |
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; | |
1304 | } | |
1305 | ||
1306 | return betaPerDOFoptions[betaOpt]; | |
1307 | ||
1308 | } | |
1309 | ||
31d3e4f0 | 1310 | //-------------------------------------------------------------------------------------------------------------------------------------------------- |
1311 | Double_t AliAnaChargedJetResponseMaker::InterpolateFast(TGraph *gr, Double_t x) { | |
1312 | ||
1313 | Double_t ipmax = gr->GetN()-1.; | |
1314 | Double_t x2,y2,xold,yold; | |
1315 | ||
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; | |
1320 | ||
1321 | Double_t ip = ipmax/2.; | |
1322 | Double_t ipold = 0.; | |
1323 | gr->GetPoint((int)(ip),x2,y2); | |
1324 | ||
1325 | Int_t i = 0; | |
1326 | ||
1327 | if(x2>x) { | |
1328 | while(x2>x) { | |
1329 | if(i==0) ipold=0.; | |
1330 | ip -= (ip)/2.; | |
1331 | gr->GetPoint((int)(ip),x2,y2); | |
1332 | if(x2>x){} | |
1333 | else ipold = ip; | |
1334 | i++; | |
1335 | // cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl; | |
1336 | } | |
1337 | } | |
1338 | else if(x2<x) { | |
1339 | while(x2<x) { | |
1340 | if(i==0) ipold=ipmax; | |
1341 | ip += (ipold-ip)/2.; | |
1342 | gr->GetPoint((int)(ip),x2,y2); | |
1343 | if(x2>x) ipold = ip; | |
1344 | else {} | |
1345 | i++; | |
1346 | // cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl; | |
1347 | } | |
1348 | } | |
1349 | ||
1350 | int ip2 = 0; | |
1351 | if(x2>x) { | |
1352 | ip2 = (int)(ip)-1; | |
1353 | gr->GetPoint(ip2,x2,y2); | |
1354 | while(x2>x) { | |
1355 | ip2--; | |
1356 | gr->GetPoint(ip2,x2,y2); | |
1357 | } | |
1358 | gr->GetPoint(ip2+1,xold,yold); | |
1359 | } | |
1360 | else { | |
1361 | ip2 = (int)(ip)+1; | |
1362 | gr->GetPoint(ip2,x2,y2); | |
1363 | while(x2<x) { | |
1364 | ip2++; | |
1365 | gr->GetPoint(ip2,x2,y2); | |
1366 | } | |
1367 | gr->GetPoint(ip2-1,xold,yold); | |
1368 | ||
1369 | } | |
1370 | // cout << "For x=" << x << " interpolate between: " << xold << " and " << x2 << endl; | |
1371 | return ((x-xold)*y2 + (x2-x)*yold) / (x2-xold); | |
1372 | ||
1373 | } | |
1374 | ||
1375 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
1376 | Double_t AliAnaChargedJetResponseMaker::InterpolateFast(TH1 *h, Double_t x) { | |
1377 | ||
1378 | // Double_t ipmax = gr->GetN()-1.; | |
1379 | Double_t ipmax = (double)h->GetNbinsX(); | |
1380 | Double_t x2,y2,xold,yold; | |
1381 | ||
1382 | Double_t xmin = h->GetXaxis()->GetXmin(); | |
1383 | Double_t xmax = h->GetXaxis()->GetXmax(); | |
1384 | if(x<xmin || x>xmax) return 0; | |
1385 | ||
1386 | Double_t ip = ipmax/2.; | |
1387 | Double_t ipold = 0.; | |
1388 | ||
1389 | x2 = h->GetBinCenter((int)ip); | |
1390 | y2 = h->GetBinContent((int)ip); | |
1391 | ||
1392 | Int_t i = 0; | |
1393 | ||
1394 | if(x2>x) { | |
1395 | while(x2>x) { | |
1396 | if(i==0) ipold=0.; | |
1397 | ip -= (ip)/2.; | |
1398 | x2 = h->GetBinCenter((int)ip); | |
1399 | y2 = h->GetBinContent((int)ip); | |
1400 | if(x2>x) {} | |
1401 | else ipold = ip; | |
1402 | i++; | |
1403 | // cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl; | |
1404 | } | |
1405 | } | |
1406 | else if(x2<x) { | |
1407 | while(x2<x) { | |
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; | |
1413 | else {} | |
1414 | i++; | |
1415 | // cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl; | |
1416 | } | |
1417 | } | |
1418 | ||
1419 | int ip2 = 0; | |
1420 | if(x2>x) { | |
1421 | ip2 = (int)(ip)-1; | |
1422 | x2 = h->GetBinCenter(ip2); | |
1423 | y2 = h->GetBinContent(ip2); | |
1424 | while(x2>x) { | |
1425 | ip2--; | |
1426 | x2 = h->GetBinCenter(ip2); | |
1427 | y2 = h->GetBinContent(ip2); | |
1428 | } | |
1429 | xold = h->GetBinCenter(ip2+1); | |
1430 | yold = h->GetBinContent(ip2+1); | |
1431 | } | |
1432 | else { | |
1433 | ip2 = (int)(ip)+1; | |
1434 | x2 = h->GetBinCenter(ip2); | |
1435 | y2 = h->GetBinContent(ip2); | |
1436 | while(x2<x) { | |
1437 | ip2++; | |
1438 | x2 = h->GetBinCenter(ip2); | |
1439 | y2 = h->GetBinContent(ip2); | |
1440 | } | |
1441 | xold = h->GetBinCenter(ip2-1); | |
1442 | yold = h->GetBinContent(ip2-1); | |
1443 | ||
1444 | } | |
1445 | // cout << "For x=" << x << " interpolate between: " << xold << " and " << x2 << endl; | |
1446 | return ((x-xold)*y2 + (x2-x)*yold) / (x2-xold); | |
1447 | ||
1448 | } | |
1449 | ||
1450 |