]>
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++) { | |
193 | double pT[1]; | |
194 | pT[0]= hPtMeasured->GetBinCenter(i); | |
195 | fPtMeasured->SetBinContent(bin,(Double_t)hPtMeasured->GetBinContent(i)); | |
196 | fPtMeasured->SetBinError(bin,(Double_t)hPtMeasured->GetBinError(i)); | |
197 | bin[0]++; | |
198 | } | |
199 | ||
200 | if(fDebug) printf("fPtMeasured->GetNbins(): %d \n",(int)fPtMeasured->GetNbins()); | |
201 | ||
202 | } | |
203 | ||
204 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
205 | void AliAnaChargedJetResponseMaker::SetFlatEfficiency(Double_t eff) { | |
206 | ||
ef62323a | 207 | // |
208 | // Set flat efficiency to value of eff | |
209 | // | |
210 | ||
31d3e4f0 | 211 | fEffFlat = eff; |
ef62323a | 212 | return; |
31d3e4f0 | 213 | |
275d3481 | 214 | /* |
31d3e4f0 | 215 | Int_t nbins[fDimensions]; |
216 | Double_t xmin[fDimensions]; | |
217 | Double_t xmax[fDimensions]; | |
218 | for(int dim = 0; dim<fDimensions; dim++) { | |
219 | nbins[dim] = fNbins; | |
220 | xmin[dim] = fPtMin; | |
221 | xmax[dim] = fPtMax; | |
222 | } | |
223 | ||
224 | if(fEfficiency) delete fEfficiency; | |
225 | fEfficiency = new THnSparseD("fEfficiency","Efficiency - no resolution effects; p_{T}^{gen}",fDimensions,nbins,xmin,xmax); | |
226 | fEfficiency->Sumw2(); | |
227 | fEfficiency->GetAxis(0)->SetTitle("p_{T}^{gen}"); | |
228 | fEfficiency->SetBinEdges(0,fBinArrayPtRec); | |
229 | ||
230 | for(int i=0; i<fNbins; i++) { | |
231 | int bin[1]; | |
232 | bin[0] = i; | |
233 | fEfficiency->SetBinContent(bin,fEffFlat); | |
234 | fEfficiency->SetBinError(bin,0); | |
235 | } | |
5b045191 | 236 | */ |
31d3e4f0 | 237 | } |
238 | ||
239 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
240 | void AliAnaChargedJetResponseMaker::SetEfficiency(TGraphErrors *grEff) | |
241 | { | |
ef62323a | 242 | // |
243 | // Fill fEfficiency (THnSparse) with values from grEff | |
244 | // | |
245 | ||
31d3e4f0 | 246 | Int_t nbins[fDimensions]; |
247 | Double_t xmin[fDimensions]; | |
248 | Double_t xmax[fDimensions]; | |
249 | for(int dim = 0; dim<fDimensions; dim++) { | |
250 | nbins[dim] = fNbins; | |
251 | xmin[dim] = fPtMin; | |
252 | xmax[dim] = fPtMax; | |
253 | } | |
254 | ||
255 | if(fEfficiency) delete fEfficiency; | |
256 | fEfficiency = new THnSparseD("fEfficiency","Efficiency - no resolution effects; p_{T}^{gen}",fDimensions,nbins,xmin,xmax); | |
257 | fEfficiency->Sumw2(); | |
258 | fEfficiency->GetAxis(0)->SetTitle("p_{T}^{gen}"); | |
259 | fEfficiency->SetBinEdges(0,fBinArrayPtRec); | |
260 | ||
261 | double pT[1]; | |
262 | double yield = 0.; | |
263 | double error = 0.; | |
264 | double dummy = 0.; | |
265 | int nbinsgr = grEff->GetN(); | |
266 | ||
267 | for(int i=0; i<nbinsgr; i++) { | |
268 | grEff->GetPoint(i,dummy,yield); | |
269 | pT[0] = dummy; | |
270 | error = grEff->GetErrorY(i); | |
271 | ||
272 | fEfficiency->Fill(pT,yield); | |
273 | fEfficiency->SetBinError(i,error); | |
274 | ||
275 | } | |
276 | ||
277 | } | |
278 | ||
279 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
280 | void AliAnaChargedJetResponseMaker::MakeResponseMatrixJetsFineMerged(Int_t skipBins, Int_t binWidthFactor, Int_t extraBins, Bool_t bVariableBinning, Double_t ptmax) | |
281 | { | |
282 | // | |
283 | // Make jet response matrix | |
284 | // | |
285 | ||
ef62323a | 286 | if(fDebug) printf(">>AliAnaChargedJetResponseMaker::MakeResponseMatrixJetsFineMerged\n"); |
287 | ||
288 | if(!fPtMeasured) { | |
289 | if(fDebug) printf("fPtMeasured does not exist. Aborting!!\n"); | |
290 | return; | |
291 | } | |
31d3e4f0 | 292 | if(fResponseMatrix) { delete fResponseMatrix; } |
293 | if(fResponseMatrixFine) { delete fResponseMatrixFine; } | |
294 | ||
295 | SetSkipBinsUnfolded(skipBins); | |
296 | SetBinWidthFactorUnfolded(binWidthFactor); | |
297 | SetExtraBinsUnfolded(extraBins); | |
298 | SetVariableBinning(bVariableBinning,ptmax); | |
299 | ||
300 | InitializeResponseMatrix(); | |
301 | InitializeEfficiency(); | |
302 | ||
303 | InitializeResponseMatrixFine(); | |
304 | InitializeEfficiencyFine(); | |
305 | ||
306 | FillResponseMatrixFineAndMerge(); | |
307 | ||
308 | } | |
309 | ||
310 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
311 | void AliAnaChargedJetResponseMaker::InitializeResponseMatrix() { | |
312 | // | |
ef62323a | 313 | //Define bin edges of RM to be used for unfolding |
31d3e4f0 | 314 | // |
315 | ||
316 | Int_t nbins[fDimensions*2]; | |
317 | nbins[fDimRec] = fNbins; | |
fa7c34ba | 318 | nbins[fDimGen] = fNbins; |
31d3e4f0 | 319 | |
320 | double binWidthMeas = (double)((fPtMax-fPtMin)/fNbins); | |
321 | double binWidthUnf = (double)fBinWidthFactorUnfolded*binWidthMeas; | |
322 | double binWidthUnfLowPt = -1.; | |
323 | if(fbVariableBinning) | |
324 | binWidthUnfLowPt = binWidthUnf*0.5; | |
325 | ||
03372fd1 | 326 | fPtMaxUnfolded = fPtMax+(double)(fExtraBinsUnfolded)*binWidthUnf; |
327 | nbins[fDimGen]+=fExtraBinsUnfolded; | |
328 | ||
31d3e4f0 | 329 | |
330 | printf("fPtMinMeas: %f fPtMaxMeas: %f\n",fPtMin,fPtMax); | |
331 | printf("binWidthMeas: %f binWidthUnf: %f fBinWidthFactorUnfolded: %d\n",binWidthMeas,binWidthUnf,fBinWidthFactorUnfolded); | |
332 | printf("binWidthUnfLowPt: %f\n",binWidthUnfLowPt); | |
333 | ||
adf3920d | 334 | printf("fPtMinUnfolded: %f fPtMaxUnfolded: %f\n",fPtMinUnfolded,fPtMaxUnfolded); |
335 | ||
31d3e4f0 | 336 | |
337 | if(fbVariableBinning) { | |
adf3920d | 338 | // cout << "fPtMaxUnfVarBinning: " << fPtMaxUnfVarBinning << " \tfPtMinUnfolded: " << fPtMinUnfolded << " binWidthUnfLowPt: " << binWidthUnfLowPt << endl; |
339 | Int_t tmp = (int)((fPtMaxUnfVarBinning-fPtMinUnfolded)/binWidthUnfLowPt); | |
31d3e4f0 | 340 | tmp = tmp - fSkipBinsUnfolded; |
341 | fPtMinUnfolded = fPtMaxUnfVarBinning-(double)(tmp)*binWidthUnfLowPt; | |
adf3920d | 342 | //cout << "tmp = " << tmp << " fSkipBinsUnfolded = " << fSkipBinsUnfolded << " fPtMinUnfolded = " << fPtMinUnfolded << endl; |
31d3e4f0 | 343 | //Redefine also number of bins on generated axis in case of variable binning |
344 | nbins[fDimGen] = (int)((fPtMaxUnfVarBinning-fPtMinUnfolded)/binWidthUnfLowPt)+(int)((fPtMaxUnfolded-fPtMaxUnfVarBinning)/binWidthUnf); | |
345 | } | |
adf3920d | 346 | else { |
347 | int tmp = round((fPtMaxUnfolded-fPtMinUnfolded)/binWidthUnf); //nbins which fit between 0 and fPtMaxUnfolded | |
348 | tmp = tmp - fSkipBinsUnfolded; | |
349 | fPtMinUnfolded = fPtMaxUnfolded-(double)(tmp)*binWidthUnf; //set ptmin unfolded | |
350 | fPtMaxUnfolded = fPtMinUnfolded+(double)(tmp)*binWidthUnf; //set ptmax unfolded | |
351 | nbins[fDimGen] = (int)((fPtMaxUnfolded-fPtMinUnfolded)/binWidthUnf); //adjust nbins to bin width | |
352 | } | |
353 | ||
354 | printf(" nbins[fDimGen] = %d nbins[fDimRec] = %d\n",nbins[fDimGen],nbins[fDimRec]); | |
31d3e4f0 | 355 | |
356 | Double_t binWidth[2]; | |
357 | binWidth[fDimRec] = (double)((fPtMax-fPtMin)/nbins[fDimRec]); | |
358 | binWidth[fDimGen] = (double)((fPtMaxUnfolded-fPtMinUnfolded)/nbins[fDimGen]); | |
359 | ||
360 | Double_t xmin[fDimensions*2]; | |
361 | Double_t xmax[fDimensions*2]; | |
362 | xmin[fDimRec] = fPtMin; | |
363 | xmax[fDimRec] = fPtMax; | |
364 | xmin[fDimGen] = fPtMinUnfolded; | |
365 | xmax[fDimGen] = fPtMaxUnfolded; | |
366 | ||
367 | printf("xmin[fDimRec]: %f xmin[fDimGen]: %f \n",xmin[fDimRec],xmin[fDimGen]); | |
368 | printf("xmax[fDimRec]: %f xmax[fDimGen]: %f \n",xmax[fDimRec],xmax[fDimGen]); | |
369 | printf("nbins[fDimRec]: %d nbins[fDimGen]: %d \n",nbins[fDimRec],nbins[fDimGen]); | |
370 | if(!fbVariableBinning) printf("binWidth[fDimRec]: %f binWidth[fDimGen]: %f \n",binWidth[fDimRec],binWidth[fDimGen]); | |
371 | ||
372 | Double_t binArrayPt0[nbins[fDimRec]+1]; | |
373 | Double_t binArrayPt1[nbins[fDimGen]+1]; | |
374 | Double_t xmaxGen = TMath::Max(xmax[fDimGen],fPtMaxUnfoldedHigh); | |
375 | ||
376 | //Define bin limits reconstructed/measured axis | |
377 | for(int i=0; i<nbins[fDimRec]; i++) { | |
378 | binArrayPt0[i] = xmin[fDimRec]+(double)i*binWidth[fDimRec]; | |
379 | } | |
380 | binArrayPt0[nbins[fDimRec]]= xmax[fDimRec]; | |
381 | ||
382 | //Define bin limits generated/unfolded axis | |
fa7c34ba | 383 | binArrayPt1[0] = fPtMinUnfolded; |
384 | for(int i=1; i<nbins[fDimGen]; i++) { | |
31d3e4f0 | 385 | if(fbVariableBinning) { |
386 | double test = xmin[fDimGen]+(double)i*binWidthUnfLowPt; | |
387 | if(test<=fPtMaxUnfVarBinning) binArrayPt1[i] = test; | |
388 | else binArrayPt1[i] = binArrayPt1[i-1]+binWidthUnf; | |
389 | } | |
fa7c34ba | 390 | else |
391 | binArrayPt1[i] = xmin[fDimGen]+(double)i*binWidth[fDimGen]; | |
31d3e4f0 | 392 | //printf("RM. i = %d \t binArrayPt[i] = %f \n",i,binArrayPt1[i]); |
393 | } | |
394 | binArrayPt1[nbins[fDimGen]]= xmaxGen; | |
395 | ||
396 | ||
397 | // Response matrix : dimensions must be 2N = 2 x (number of variables) | |
398 | // dimensions 0 -> N-1 must be filled with reconstructed values | |
399 | // dimensions N -> 2N-1 must be filled with generated values | |
ef62323a | 400 | if(fResponseMatrix) delete fResponseMatrix; |
31d3e4f0 | 401 | fResponseMatrix = new THnSparseD("fResponseMatrix","Response Matrix pTMC vs pTrec",fDimensions*2,nbins,xmin,xmax); |
402 | fResponseMatrix->Sumw2(); | |
403 | fResponseMatrix->GetAxis(fDimRec)->SetTitle("p_{T}^{rec}"); | |
404 | fResponseMatrix->GetAxis(fDimGen)->SetTitle("p_{T}^{gen}"); | |
405 | ||
406 | fResponseMatrix->SetBinEdges(fDimRec,binArrayPt0); | |
407 | fResponseMatrix->SetBinEdges(fDimGen,binArrayPt1); | |
408 | ||
ef62323a | 409 | Int_t bin[2] = {1,1}; |
410 | for(int i=1; i<fResponseMatrix->GetAxis(0)->GetNbins(); i++) { | |
411 | bin[0]=i; | |
412 | for(int j=1; j<fResponseMatrix->GetAxis(1)->GetNbins(); j++) { | |
413 | bin[1]=j; | |
414 | fResponseMatrix->SetBinContent(bin,0.); | |
415 | } | |
416 | } | |
417 | ||
418 | ||
31d3e4f0 | 419 | |
420 | } | |
421 | ||
422 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
423 | void AliAnaChargedJetResponseMaker::InitializeEfficiency() { | |
424 | // | |
425 | // Define dimensions of efficiency THnSparse | |
426 | // | |
427 | ||
428 | if(!fResponseMatrix) { | |
429 | printf("AliAnaChargedJetResponseMaker::InitializeEfficiency()\n"); | |
430 | printf("Can not define dimensions efficiency without fResponseMatrix dimensions. Aborting \n"); | |
431 | return; | |
432 | } | |
433 | ||
434 | TAxis *genAxis = fResponseMatrix->GetAxis(fDimGen); | |
435 | ||
436 | Int_t nbinsEff[fDimensions]; | |
437 | Double_t xminEff[fDimensions]; | |
438 | Double_t xmaxEff[fDimensions]; | |
439 | ||
440 | for(int dim = 0; dim<fDimensions; dim++) { | |
441 | nbinsEff[dim] = genAxis->GetNbins(); | |
442 | xminEff[dim] = genAxis->GetXmin(); | |
443 | xmaxEff[dim] = genAxis->GetXmax(); | |
444 | } | |
445 | ||
446 | if(fEfficiency) delete fEfficiency; | |
447 | fEfficiency = new THnSparseD("fEfficiency","Efficiency - no resolution effects; p_{T}^{gen}",fDimensions,nbinsEff,xminEff,xmaxEff); | |
448 | fEfficiency->Sumw2(); | |
449 | fEfficiency->GetAxis(0)->SetTitle("p_{T}^{gen}"); | |
450 | ||
451 | const Double_t *binArrayPt = genAxis->GetXbins()->GetArray(); | |
452 | fEfficiency->SetBinEdges(0,binArrayPt); | |
453 | ||
454 | } | |
455 | ||
456 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
457 | void AliAnaChargedJetResponseMaker::InitializeResponseMatrixFine() { | |
458 | // | |
459 | //Initialize fine response matrix | |
460 | // | |
461 | ||
462 | Int_t nbinsFine[fDimensions*2]; | |
fa7c34ba | 463 | Double_t xminFine[fDimensions*2]; |
464 | Double_t xmaxFine[fDimensions*2]; | |
31d3e4f0 | 465 | Double_t pTarrayFine[fDimensions*2]; |
466 | ||
fa7c34ba | 467 | nbinsFine[fDimRec] = fResponseMatrix->GetAxis(fDimRec)->GetNbins()*fFineFrac; |
468 | nbinsFine[fDimGen] = fResponseMatrix->GetAxis(fDimRec)->GetNbins()*fFineFrac; | |
469 | xminFine[fDimRec] = TMath::Min(fPtMin,0.); | |
470 | xminFine[fDimGen] = TMath::Min(fPtMin,0.); | |
0bf6381d | 471 | xmaxFine[fDimRec] = fResponseMatrix->GetAxis(fDimGen)->GetXmax();//+40.; |
472 | xmaxFine[fDimGen] = fResponseMatrix->GetAxis(fDimGen)->GetXmax();//+40.; | |
fa7c34ba | 473 | pTarrayFine[fDimRec] = 0.; |
474 | pTarrayFine[fDimGen] = 0.; | |
31d3e4f0 | 475 | |
476 | Double_t binWidth[2]; | |
477 | binWidth[fDimRec] = fResponseMatrix->GetAxis(fDimRec)->GetBinWidth(1); | |
478 | ||
479 | Double_t binWidthFine[2]; | |
480 | binWidthFine[fDimRec] = binWidth[fDimRec]/((double)fFineFrac); | |
481 | binWidthFine[fDimGen] = binWidthFine[fDimRec]*(double)fBinWidthFactorUnfolded; | |
0bf6381d | 482 | // binWidthFine[fDimRec] = binWidthFine[fDimGen]; |
31d3e4f0 | 483 | |
484 | nbinsFine[fDimRec] = (int)((xmaxFine[fDimRec]-xminFine[fDimRec])/binWidthFine[fDimRec]); //adjust nbins to bin width | |
485 | nbinsFine[fDimGen] = (int)((xmaxFine[fDimGen]-xminFine[fDimGen])/binWidthFine[fDimGen]); //adjust nbins to bin width | |
486 | ||
487 | printf("xminFine[fDimRec]: %f xminFine[fDimGen]: %f \n",xminFine[fDimRec],xminFine[fDimGen]); | |
488 | printf("xmaxFine[fDimRec]: %f xmaxFine[fDimGen]: %f \n",xmaxFine[fDimRec],xmaxFine[fDimGen]); | |
489 | printf("nbinsFine[fDimRec]: %d nbinsFine[fDimGen]: %d \n",nbinsFine[fDimRec],nbinsFine[fDimGen]); | |
490 | printf("binWidthFine[fDimRec]: %f binWidthFine[fDimGen]: %f \n",binWidthFine[fDimRec],binWidthFine[fDimGen]); | |
491 | ||
492 | ||
493 | Double_t binArrayPt0Fine[nbinsFine[fDimRec]+1]; | |
494 | Double_t binArrayPt1Fine[nbinsFine[fDimGen]+1]; | |
495 | ||
496 | for(int i=0; i<nbinsFine[fDimRec]; i++) { | |
497 | binArrayPt0Fine[i] = xminFine[fDimRec]+(double)i*binWidthFine[fDimRec]; | |
498 | // printf("RM. i = %d \t binArrayPtFine[i] = %f \n",i,binArrayPt0Fine[i]); | |
499 | } | |
500 | binArrayPt0Fine[nbinsFine[fDimRec]]= xmaxFine[fDimRec]; | |
501 | ||
502 | for(int i=0; i<nbinsFine[fDimGen]; i++) { | |
503 | binArrayPt1Fine[i] = xminFine[fDimGen]+(double)i*binWidthFine[fDimGen]; | |
504 | // printf("RM. i = %d \t binArrayPtFine[i] = %f \n",i,binArrayPt1Fine[i]); | |
505 | } | |
506 | binArrayPt1Fine[nbinsFine[fDimGen]]= xmaxFine[fDimGen]; | |
507 | ||
508 | // Response matrix : dimensions must be 2N = 2 x (number of variables) | |
509 | // dimensions 0 -> N-1 must be filled with reconstructed values | |
510 | // dimensions N -> 2N-1 must be filled with generated values | |
ef62323a | 511 | if(fResponseMatrixFine) delete fResponseMatrixFine; |
31d3e4f0 | 512 | fResponseMatrixFine = new THnSparseD("fResponseMatrixFine","Response Matrix pTMC vs pTrec",fDimensions*2,nbinsFine,xminFine,xmaxFine); |
513 | fResponseMatrixFine->Sumw2(); | |
514 | fResponseMatrixFine->GetAxis(fDimRec)->SetTitle("p_{T}^{rec}"); | |
515 | fResponseMatrixFine->GetAxis(fDimGen)->SetTitle("p_{T}^{gen}"); | |
516 | ||
517 | fResponseMatrixFine->SetBinEdges(fDimRec,binArrayPt0Fine); | |
518 | fResponseMatrixFine->SetBinEdges(fDimGen,binArrayPt1Fine); | |
519 | ||
ef62323a | 520 | Int_t bin[2] = {1,1}; |
521 | for(int i=1; i<fResponseMatrixFine->GetAxis(0)->GetNbins(); i++) { | |
522 | bin[0]=i; | |
523 | for(int j=1; j<fResponseMatrixFine->GetAxis(1)->GetNbins(); j++) { | |
524 | bin[1]=j; | |
525 | fResponseMatrixFine->SetBinContent(bin,0.); | |
526 | } | |
527 | } | |
528 | ||
31d3e4f0 | 529 | } |
530 | ||
531 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
532 | void AliAnaChargedJetResponseMaker::InitializeEfficiencyFine() { | |
533 | // | |
534 | // Define dimensions of efficiency THnSparse | |
535 | // | |
536 | ||
537 | if(!fResponseMatrixFine) { | |
538 | printf("AliAnaChargedJetResponseMaker::InitializeEfficiencyFine()\n"); | |
539 | printf("Can not define dimensions efficiency without fResponseMatrixFine dimensions. Aborting \n"); | |
540 | return; | |
541 | } | |
542 | ||
543 | TAxis *genAxis = fResponseMatrixFine->GetAxis(fDimGen); | |
544 | ||
545 | Int_t nbinsEff[fDimensions]; | |
546 | Double_t xminEff[fDimensions]; | |
547 | Double_t xmaxEff[fDimensions]; | |
548 | ||
549 | for(int dim = 0; dim<fDimensions; dim++) { | |
550 | nbinsEff[dim] = genAxis->GetNbins(); | |
551 | xminEff[dim] = genAxis->GetXmin(); | |
552 | xmaxEff[dim] = genAxis->GetXmax(); | |
553 | } | |
554 | ||
555 | if(fEfficiencyFine) delete fEfficiencyFine; | |
556 | fEfficiencyFine = new THnSparseD("fEfficiencyFine","EfficiencyFine - no resolution effects; p_{T}^{gen}",fDimensions,nbinsEff,xminEff,xmaxEff); | |
557 | fEfficiencyFine->Sumw2(); | |
558 | fEfficiencyFine->GetAxis(0)->SetTitle("p_{T}^{gen}"); | |
559 | ||
560 | const Double_t *binArrayPt = genAxis->GetXbins()->GetArray(); | |
561 | fEfficiencyFine->SetBinEdges(0,binArrayPt); | |
562 | ||
563 | } | |
564 | ||
565 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
566 | void AliAnaChargedJetResponseMaker::FillResponseMatrixFineAndMerge() { | |
567 | // | |
568 | // Fill fine response matrix | |
569 | // | |
570 | ||
ef62323a | 571 | if(!fResponseMatrix) { |
572 | printf("Dimensions of fResponseMatrix have to be defined first. Aborting!"); | |
573 | return; | |
574 | } | |
575 | ||
576 | if(!fResponseMatrixFine) { | |
577 | printf("Dimensions of fResponseMatrixFine have to be defined first. Aborting!"); | |
578 | return; | |
579 | } | |
580 | ||
31d3e4f0 | 581 | TAxis *genAxis = fResponseMatrixFine->GetAxis(fDimGen); |
582 | TAxis *recAxis = fResponseMatrixFine->GetAxis(fDimRec); | |
583 | ||
584 | Int_t nbinsFine[fDimensions*2]; | |
585 | Double_t xminFine[fDimensions*2]; | |
586 | Double_t xmaxFine[fDimensions*2]; | |
587 | Double_t pTarrayFine[fDimensions*2]; | |
588 | ||
589 | nbinsFine[fDimGen] = genAxis->GetNbins(); | |
590 | nbinsFine[fDimRec] = recAxis->GetNbins(); | |
591 | xminFine[fDimGen] = genAxis->GetXmin(); | |
592 | xminFine[fDimRec] = recAxis->GetXmin(); | |
593 | xmaxFine[fDimGen] = genAxis->GetXmax(); | |
594 | xmaxFine[fDimRec] = recAxis->GetXmax(); | |
fa7c34ba | 595 | pTarrayFine[fDimGen] = 0.; |
596 | pTarrayFine[fDimRec] = 0.; | |
31d3e4f0 | 597 | |
598 | double sumyield = 0.; | |
599 | double sumyielderror2 = 0.; | |
600 | ||
ada4df5e | 601 | int ipt[2] = {0,0}; |
602 | int iptMerged[2] = {0,0}; | |
603 | int ierror[2] = {0,0}; | |
31d3e4f0 | 604 | |
605 | Int_t check = 0; | |
606 | Double_t pTgen, pTrec; | |
607 | Double_t yield = 0.; | |
608 | Double_t error = 0.; | |
609 | ||
610 | const int nng = fResponseMatrix->GetAxis(fDimGen)->GetNbins(); | |
611 | const int nnr = fResponseMatrix->GetAxis(fDimRec)->GetNbins(); | |
612 | Double_t errorArray[nng][nnr]; | |
613 | for(int iig =0; iig<nng; iig++) { | |
614 | for(int iir =0; iir<nnr; iir++) { | |
615 | errorArray[iig][iir] = 0.; | |
616 | } | |
617 | } | |
618 | ||
619 | for(int iy=1; iy<=nbinsFine[fDimGen]; iy++) { //gen | |
620 | pTgen = fResponseMatrixFine->GetAxis(fDimGen)->GetBinCenter(iy); | |
621 | pTarrayFine[fDimGen] = pTgen; | |
622 | ierror[fDimGen]=iy; | |
623 | sumyield = 0.; | |
624 | check = 0; | |
625 | ||
626 | for(int ix=1; ix<=nbinsFine[fDimRec]; ix++) { //rec | |
627 | pTrec = fResponseMatrixFine->GetAxis(fDimRec)->GetBinCenter(ix); | |
628 | Double_t width = fResponseMatrixFine->GetAxis(fDimRec)->GetBinWidth(ix); | |
629 | if(fResolutionType==kParam) { | |
630 | yield = fDeltaPt->Eval(pTrec-pTgen); | |
631 | error = 0.; | |
632 | } | |
633 | else if(fResolutionType==kResiduals) { | |
634 | yield = fhDeltaPt->Interpolate(pTrec-pTgen); | |
635 | error = 0.; | |
636 | } | |
637 | else if(fResolutionType==kResidualsErr) { | |
638 | Double_t deltaPt = pTrec-pTgen; | |
639 | int bin = fhDeltaPt->FindBin(deltaPt); | |
640 | yield = fhDeltaPt->GetBinContent(bin); | |
641 | error = fhDeltaPt->GetBinError(bin); | |
642 | } | |
643 | yield=yield*width; | |
644 | error=error*width; | |
645 | //avoid trouble with empty bins in the high pT tail of deltaPt distribution | |
646 | if(check==0 && yield>0. && pTrec>pTgen) check=1; | |
647 | if(check==1 && yield==0.) ix=nbinsFine[fDimRec]; | |
648 | ||
649 | sumyield+=yield; | |
650 | sumyielderror2 += error*error; | |
651 | ||
652 | pTarrayFine[fDimRec] = pTrec; | |
653 | ierror[fDimRec] = ix; | |
654 | fResponseMatrixFine->Fill(pTarrayFine,yield); | |
655 | if(fbCalcErrors) fResponseMatrixFine->SetBinError(ierror,error); | |
656 | ||
657 | }// ix (dimRec) loop | |
658 | ||
659 | //Normalize to total number of counts =1 | |
660 | if(sumyield>1) { | |
661 | ipt[fDimGen]=iy; | |
662 | for(int ix=1; ix<=nbinsFine[fDimRec]; ix++) { | |
663 | ipt[fDimRec]=ix; | |
664 | fResponseMatrixFine->SetBinContent(ipt,fResponseMatrixFine->GetBinContent(ipt)/sumyield); | |
665 | if(fResolutionType==kResidualsErr && fbCalcErrors) { | |
666 | Double_t A = 1./sumyield*fResponseMatrix->GetBinError(ipt); | |
667 | Double_t B = -1.*fResponseMatrix->GetBinContent(ipt)/sumyield/sumyield*TMath::Sqrt(sumyielderror2); | |
668 | Double_t tmp2 = A*A + B*B; | |
669 | fResponseMatrix->SetBinError(ipt,TMath::Sqrt(tmp2)); | |
670 | } | |
671 | ||
672 | } | |
673 | } | |
674 | ||
675 | int bin[1]; | |
676 | bin[0] = iy; | |
677 | fEfficiencyFine->SetBinContent(bin,sumyield); | |
678 | if(fbCalcErrors) fEfficiencyFine->SetBinError(bin,TMath::Sqrt(sumyielderror2)); | |
679 | ||
680 | //fill merged response matrix | |
681 | ||
0bf6381d | 682 | //find bin in fine RM corresponding to mimimum/maximum bin of merged RM on rec axis |
31d3e4f0 | 683 | int ixMin = fResponseMatrixFine->GetAxis(fDimRec)->FindBin(fResponseMatrix->GetAxis(fDimRec)->GetXmin()); |
684 | int ixMax = fResponseMatrixFine->GetAxis(fDimRec)->FindBin(fResponseMatrix->GetAxis(fDimRec)->GetXmax()); | |
685 | ||
686 | if(fResponseMatrixFine->GetAxis(fDimGen)->GetBinLowEdge(iy) >= fResponseMatrix->GetAxis(fDimGen)->GetXmin()) { | |
687 | ipt[fDimGen]=iy; | |
688 | iptMerged[fDimGen]=fResponseMatrix->GetAxis(fDimGen)->FindBin(pTgen); | |
689 | ||
690 | Double_t weight = 1.; | |
691 | if(f1MergeFunction) { | |
692 | Double_t loEdge = fResponseMatrix->GetAxis(fDimGen)->GetBinLowEdge(iptMerged[fDimGen]); | |
693 | Double_t upEdge = fResponseMatrix->GetAxis(fDimGen)->GetBinUpEdge(iptMerged[fDimGen]); | |
694 | Float_t powInteg = f1MergeFunction->Integral(loEdge,upEdge); | |
695 | //printf("loEdge = %f upEdge = %f powInteg = %f\n",loEdge,upEdge,powInteg); | |
696 | if(powInteg>0.) | |
697 | weight = f1MergeFunction->Integral(fResponseMatrixFine->GetAxis(fDimGen)->GetBinLowEdge(iy),fResponseMatrixFine->GetAxis(fDimGen)->GetBinUpEdge(iy))/powInteg; | |
698 | // printf("weight: %f \n", weight ); | |
699 | } else { | |
700 | weight = 1./((double)fFineFrac); | |
adf3920d | 701 | if(fbVariableBinning && pTgen<fPtMaxUnfVarBinning) weight=1./(0.5*(double)fFineFrac); |
31d3e4f0 | 702 | } |
703 | ||
704 | for(int ix=ixMin; ix<=ixMax; ix++) { //loop reconstructed axis | |
705 | pTrec = fResponseMatrixFine->GetAxis(fDimRec)->GetBinCenter(ix); | |
706 | ipt[fDimRec]=ix; | |
707 | iptMerged[fDimRec]=fResponseMatrix->GetAxis(fDimRec)->FindBin(pTrec); | |
708 | ||
709 | fResponseMatrix->AddBinContent(iptMerged,fResponseMatrixFine->GetBinContent(ipt)*weight); | |
710 | if(fbCalcErrors) errorArray[iptMerged[fDimGen]-1][iptMerged[fDimRec]-1] += fResponseMatrixFine->GetBinError(ipt)*fResponseMatrixFine->GetBinError(ipt)*weight*weight; | |
711 | } | |
712 | ||
713 | }//loop gen axis fine response matrix | |
714 | ||
715 | } // iy (dimGen) loop | |
716 | ||
717 | //Fill Efficiency corresponding to merged response matrix | |
31d3e4f0 | 718 | for(int iy=1; iy<=fResponseMatrix->GetAxis(fDimGen)->GetNbins(); iy++) { //gen |
719 | sumyield = 0.; | |
720 | ipt[fDimGen] = iy; | |
721 | ||
722 | for(int ix=1; ix<=fResponseMatrix->GetAxis(fDimRec)->GetNbins(); ix++) { //rec | |
723 | ipt[fDimRec] = ix; | |
724 | sumyield += fResponseMatrix->GetBinContent(ipt); | |
725 | ||
726 | if(fbCalcErrors) fResponseMatrix->SetBinError(ipt,TMath::Sqrt(errorArray[iy][ix])); | |
727 | } | |
728 | printf("RM: pTgen: %f \t sumyield: %f \n",fResponseMatrix->GetAxis(fDimGen)->GetBinCenter(iy),sumyield); | |
729 | int bin[1]; | |
730 | bin[0] = iy; | |
731 | fEfficiency->SetBinContent(bin,sumyield); | |
732 | if(fbCalcErrors) fEfficiency->SetBinError(bin,0); | |
733 | } | |
734 | ||
735 | if(fDebug) printf("fResponseMatrixFine->GetNbins(): %d \n",(int)fResponseMatrixFine->GetNbins()); | |
736 | if(fDebug) printf("fResponseMatrix->GetNbins(): %d \n",(int)fResponseMatrix->GetNbins()); | |
737 | ||
738 | if(fDebug) printf("Done constructing response matrix\n"); | |
739 | ||
740 | } | |
741 | ||
ef62323a | 742 | //-------------------------------------------------------------------------------------------------------------------------------------------------- |
03372fd1 | 743 | TH2* AliAnaChargedJetResponseMaker::MakeResponseMatrixRebin(TH2 *hRMFine, TH2 *hRM, Bool_t useFunctionWeight) { |
ef62323a | 744 | |
745 | // | |
746 | // Rebin matrix hRMFine to dimensions of hRM | |
03372fd1 | 747 | // function returns matrix in TH2D format (hRM2) with dimensions from hRM |
ef62323a | 748 | // |
749 | ||
750 | TH2 *hRM2 = (TH2*)hRM->Clone("hRM2"); | |
751 | hRM2->Reset(); | |
752 | ||
03372fd1 | 753 | if(useFunctionWeight) cout << "Use function to do weighting" << endl; |
754 | ||
275d3481 | 755 | //First normalize columns of input |
ef62323a | 756 | const Int_t nbinsNorm = hRM2->GetNbinsX(); |
0bf6381d | 757 | if(fDebug) cout << "nbinsNorm: " << nbinsNorm << endl; |
ef62323a | 758 | |
759 | TArrayF *normVector = new TArrayF(nbinsNorm); | |
760 | ||
761 | for(int ix=1; ix<=hRM2->GetNbinsX(); ix++) { | |
762 | Double_t xLow = hRM2->GetXaxis()->GetBinLowEdge(ix); | |
763 | Double_t xUp = hRM2->GetXaxis()->GetBinUpEdge(ix); | |
764 | //cout << "xLow: " << xLow << " xUp: " << xUp << "\t center: " << hRM2->GetXaxis()->GetBinCenter(ix) << endl; | |
765 | Int_t binxLowFine = hRMFine->GetXaxis()->FindBin(xLow); | |
766 | Int_t binxUpFine = hRMFine->GetXaxis()->FindBin(xUp)-1; | |
767 | //cout << "xLowFine: " << hRMFine->GetXaxis()->GetBinLowEdge(binxLowFine) << "\txUpFine: " << hRMFine->GetXaxis()->GetBinUpEdge(binxUpFine) << endl; | |
03372fd1 | 768 | if(useFunctionWeight) |
769 | normVector->SetAt(f1MergeFunction->Integral(xLow,xUp),ix-1); | |
770 | else | |
771 | normVector->SetAt(hRMFine->Integral(binxLowFine,binxUpFine,1,hRMFine->GetYaxis()->GetNbins()),ix-1); | |
0bf6381d | 772 | //if(fDebug) cout << "ix norm: " << normVector->At(ix-1) << endl; |
ef62323a | 773 | } |
774 | ||
775 | Double_t content, oldcontent = 0.; | |
776 | Int_t ixNew = 0; | |
777 | Int_t iyNew = 0; | |
275d3481 | 778 | Double_t xvalLo, xvalUp, yvalLo, yvalUp; |
ef62323a | 779 | Double_t xmin = hRM2->GetXaxis()->GetXmin(); |
780 | Double_t ymin = hRM2->GetYaxis()->GetXmin(); | |
781 | Double_t xmax = hRM2->GetXaxis()->GetXmax(); | |
782 | Double_t ymax = hRM2->GetYaxis()->GetXmax(); | |
783 | for(int ix=1; ix<=hRMFine->GetXaxis()->GetNbins(); ix++) { | |
275d3481 | 784 | xvalLo = hRMFine->GetXaxis()->GetBinLowEdge(ix); |
785 | xvalUp = hRMFine->GetXaxis()->GetBinUpEdge(ix); | |
786 | if(xvalLo<xmin || xvalUp>xmax) continue; | |
787 | ixNew = hRM2->GetXaxis()->FindBin(hRMFine->GetXaxis()->GetBinCenter(ix)); | |
ef62323a | 788 | for(int iy=1; iy<=hRMFine->GetYaxis()->GetNbins(); iy++) { |
275d3481 | 789 | yvalLo = hRMFine->GetYaxis()->GetBinLowEdge(iy); |
790 | yvalUp = hRMFine->GetYaxis()->GetBinUpEdge(iy); | |
791 | if(yvalLo<ymin || yvalUp>ymax) continue; | |
ef62323a | 792 | content = hRMFine->GetBinContent(ix,iy); |
275d3481 | 793 | iyNew = hRM2->GetYaxis()->FindBin(hRMFine->GetYaxis()->GetBinCenter(iy)); |
ef62323a | 794 | oldcontent = hRM2->GetBinContent(ixNew,iyNew); |
795 | ||
275d3481 | 796 | //if(fDebug) cout << "ixNew: " << ixNew << " " << xvalLo << " iyNew: " << iyNew << " " << yvalLo << " content: " << content << " oldContent: " << oldcontent << " newContent: " << oldcontent+content << endl; |
ef62323a | 797 | Double_t weight = 1.; |
03372fd1 | 798 | if(normVector->At(ixNew-1)>0.) { |
799 | if(useFunctionWeight) | |
800 | weight = f1MergeFunction->Integral(xvalLo,xvalUp)/normVector->At(ixNew-1); | |
801 | else | |
802 | weight = 1./normVector->At(ixNew-1); | |
803 | } | |
ef62323a | 804 | hRM2->SetBinContent(ixNew,iyNew,oldcontent+content*weight); |
805 | } | |
806 | } | |
807 | ||
808 | if(normVector) delete normVector; | |
809 | ||
810 | return hRM2; | |
811 | ||
812 | } | |
813 | ||
03372fd1 | 814 | //-------------------------------------------------------------------------------------------------------------------------------------------------- |
815 | TH2* AliAnaChargedJetResponseMaker::CreateTruncated2DHisto(TH2 *h2, Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax) { | |
816 | ||
817 | // | |
818 | // Limit axis range of 2D histogram | |
819 | // | |
820 | ||
821 | Int_t binMinXh2 = h2->GetXaxis()->FindBin(xmin); | |
822 | if(h2->GetXaxis()->GetBinLowEdge(binMinXh2) < xmin ) binMinXh2++; | |
823 | if(h2->GetXaxis()->GetBinLowEdge(binMinXh2) > xmin ) binMinXh2--; | |
824 | ||
825 | Int_t binMinYh2 = h2->GetYaxis()->FindBin(ymin); | |
826 | if(h2->GetYaxis()->GetBinLowEdge(binMinYh2) < ymin ) binMinYh2++; | |
827 | if(h2->GetYaxis()->GetBinLowEdge(binMinYh2) > ymin ) binMinYh2--; | |
828 | ||
829 | Int_t binMaxXh2 = h2->GetXaxis()->FindBin(xmax); | |
830 | if(h2->GetXaxis()->GetBinUpEdge(binMaxXh2) < xmax ) binMaxXh2++; | |
831 | if(h2->GetXaxis()->GetBinUpEdge(binMaxXh2) > xmax ) binMaxXh2--; | |
832 | ||
833 | Int_t binMaxYh2 = h2->GetYaxis()->FindBin(ymax); | |
834 | if(h2->GetYaxis()->GetBinUpEdge(binMaxYh2) < ymax ) binMaxYh2++; | |
835 | if(h2->GetYaxis()->GetBinUpEdge(binMaxYh2) > ymax ) binMaxYh2--; | |
836 | ||
0bf6381d | 837 | Int_t nbinsX = binMaxXh2-binMinXh2+1; |
838 | Int_t nbinsY = binMaxYh2-binMinYh2+1; | |
03372fd1 | 839 | |
840 | Double_t *binsX = new Double_t[nbinsX+1]; | |
841 | Double_t *binsY = new Double_t[nbinsY+1]; | |
842 | ||
843 | for(int ix=1; ix<=nbinsX; ix++) | |
844 | binsX[ix-1] = h2->GetXaxis()->GetBinLowEdge(binMinXh2+ix-1); | |
845 | binsX[nbinsX] = h2->GetXaxis()->GetBinUpEdge(binMaxXh2); | |
846 | ||
847 | for(int iy=1; iy<=nbinsY; iy++) | |
848 | binsY[iy-1] = h2->GetYaxis()->GetBinLowEdge(binMinYh2+iy-1); | |
849 | binsY[nbinsY] = h2->GetYaxis()->GetBinUpEdge(binMaxYh2); | |
850 | ||
851 | TH2 *h2Lim = new TH2D("h2Lim","h2Lim",nbinsX,binsX,nbinsY,binsY); | |
852 | ||
853 | for(int ix=1; ix<=nbinsX; ix++) { | |
854 | // cout << "ix: " << ix << " " << binsX[ix] << endl; | |
855 | for(int iy=1; iy<=nbinsY; iy++) { | |
0bf6381d | 856 | // cout << "ix: " << ix << " " << binsX[ix] << "\tiy: " << iy << " " << binsY[iy] << endl; |
03372fd1 | 857 | |
858 | double content = h2->GetBinContent(binMinXh2+ix-1,binMinYh2+iy-1); | |
859 | double error = h2->GetBinContent(binMinXh2+ix-1,binMinYh2+iy-1); | |
860 | h2Lim->SetBinContent(ix,iy,content); | |
0bf6381d | 861 | if(fbCalcErrors) h2Lim->SetBinError(ix,iy,error); |
03372fd1 | 862 | |
863 | } | |
864 | } | |
865 | ||
866 | ||
867 | ||
868 | return h2Lim; | |
869 | } | |
870 | ||
871 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
872 | TH2* AliAnaChargedJetResponseMaker::TruncateAxisRangeResponseMatrix(TH2 *hRMOrig, Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax) { | |
873 | ||
874 | // | |
875 | // Limit axis range of response matrix without changing normalization | |
876 | // | |
877 | ||
878 | //TH2 *hRMLimit | |
879 | //TH2 *hRMLimit2 = (TH2*)hRMLimit->Clone("hRMLimit2"); | |
880 | ||
881 | TH2 *hRMLimit2 = CreateTruncated2DHisto(hRMOrig, xmin, xmax, ymin, ymax); | |
882 | hRMLimit2->Reset(); | |
883 | ||
884 | double binCent[2] = {0.,0.}; | |
885 | double content = 0.; | |
886 | double error = 0.; | |
887 | Int_t binOrig[2] = {0}; | |
888 | for(int ix=1; ix<=hRMLimit2->GetXaxis()->GetNbins(); ix++) { | |
889 | binCent[0] = hRMLimit2->GetXaxis()->GetBinCenter(ix); | |
890 | binOrig[0] = hRMOrig->GetXaxis()->FindBin(binCent[0]); | |
891 | for(int iy=1; iy<=hRMLimit2->GetYaxis()->GetNbins(); iy++) { | |
892 | binCent[1] = hRMLimit2->GetYaxis()->GetBinCenter(iy); | |
893 | binOrig[1] = hRMOrig->GetYaxis()->FindBin(binCent[1]); | |
894 | ||
895 | content = hRMOrig->GetBinContent(binOrig[0],binOrig[1]); | |
896 | error = hRMOrig->GetBinError(binOrig[0],binOrig[1]); | |
897 | ||
898 | hRMLimit2->SetBinContent(ix,iy,content); | |
0bf6381d | 899 | if(fbCalcErrors) hRMLimit2->SetBinError(ix,iy,error); |
03372fd1 | 900 | |
901 | } | |
902 | } | |
903 | ||
904 | ||
905 | return hRMLimit2; | |
906 | ||
907 | } | |
908 | ||
0bf6381d | 909 | //-------------------------------------------------------------------------------------------------------------------------------------------------- |
910 | Bool_t AliAnaChargedJetResponseMaker::CheckInputForCombinedResponse() { | |
911 | ||
912 | Bool_t check = kTRUE; | |
913 | ||
914 | if(!fhDeltaPt) check = kFALSE; | |
915 | if(!fh2DetectorResponse) check = kFALSE; | |
916 | if(!fhEfficiencyDet) check = kFALSE; | |
917 | if(!fh1MeasuredTruncated) check = kFALSE; | |
918 | if(fPtMin<0. || fPtMax<0.) check = kFALSE; | |
919 | ||
920 | return check; | |
921 | } | |
922 | ||
923 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
924 | void AliAnaChargedJetResponseMaker::MakeResponseMatrixCombined(Int_t skipBins, Int_t binWidthFactor, Int_t extraBins, Bool_t bVariableBinning, Double_t ptmin) { | |
925 | ||
926 | //Check if all input histos are there otherwise break | |
927 | if(!CheckInputForCombinedResponse()) { | |
928 | printf("Not all input available..... \n"); | |
929 | return; | |
930 | } | |
931 | ||
932 | // Make response bkg fluctuations | |
933 | MakeResponseMatrixJetsFineMerged(skipBins,binWidthFactor,extraBins,bVariableBinning,ptmin); | |
934 | ||
935 | //Get TH2D with dimensions we want final RM | |
936 | TH2D *h2ResponseMatrix = fResponseMatrix->Projection(0,1,"E"); | |
937 | h2ResponseMatrix->Reset(); | |
938 | ||
939 | //Get fine binned response matrix from bkg fluctuations | |
940 | TH2D *h2ResponseMatrixFine = fResponseMatrixFine->Projection(0,1,"E"); | |
941 | ||
942 | // Rebin response detector effects | |
943 | const TArrayD *arrayX = h2ResponseMatrixFine->GetXaxis()->GetXbins(); | |
944 | const Double_t *xbinsArray = arrayX->GetArray(); | |
945 | ||
946 | TH2D *h2ResponseDetTmp = new TH2D("h2ResponseDetTmp","h2ResponseDetTmp",h2ResponseMatrixFine->GetNbinsX(),xbinsArray,h2ResponseMatrixFine->GetNbinsX(),xbinsArray); | |
947 | ||
948 | fh2DetectorResponseRebin = (TH2D*)MakeResponseMatrixRebin(fh2DetectorResponse,h2ResponseDetTmp,kFALSE); | |
949 | fh2DetectorResponseRebin->SetName("fh2DetectorResponseRebin"); | |
950 | fh2DetectorResponseRebin->SetTitle("fh2DetectorResponseRebin"); | |
951 | fh2DetectorResponseRebin->SetXTitle("p_{T}^{jet} gen"); | |
952 | fh2DetectorResponseRebin->SetYTitle("p_{T}^{jet} rec"); | |
953 | ||
954 | // Make full combined response matrix fine binning (background flucutuations + detector effects) | |
955 | fh2ResponseMatrixCombinedFineFull = (TH2D*)MultiplityResponseMatrices(h2ResponseMatrixFine,fh2DetectorResponseRebin); | |
956 | fh2ResponseMatrixCombinedFineFull->SetName("fh2ResponseMatrixCombinedFineFull"); | |
957 | fh2ResponseMatrixCombinedFineFull->SetTitle("fh2ResponseMatrixCombinedFineFull"); | |
958 | ||
959 | // Rebin full combined response matrix with weighting | |
960 | fh2ResponseMatrixCombinedFull = (TH2D*)MakeResponseMatrixRebin(fh2ResponseMatrixCombinedFineFull,h2ResponseMatrix,kTRUE); | |
961 | fh2ResponseMatrixCombinedFull->SetName("fh2ResponseMatrixCombinedFull"); | |
962 | fh2ResponseMatrixCombinedFull->SetTitle("fh2ResponseMatrixCombinedFull"); | |
963 | fh2ResponseMatrixCombinedFull->SetXTitle("#it{p}_{T,gen}^{ch jet} (GeV/#it{c})"); | |
964 | fh2ResponseMatrixCombinedFull->SetYTitle("#it{p}_{T,rec}^{ch jet} (GeV/#it{c})"); | |
965 | ||
966 | fh2ResponseMatrixCombined = (TH2D*)TruncateAxisRangeResponseMatrix(fh2ResponseMatrixCombinedFull, h2ResponseMatrix->GetXaxis()->GetXmin(),h2ResponseMatrix->GetXaxis()->GetXmax(),fh1MeasuredTruncated->GetXaxis()->GetXmin(),fh1MeasuredTruncated->GetXaxis()->GetXmax()); | |
967 | fh2ResponseMatrixCombined->SetTitle("fh2ResponseMatrixCombined"); | |
968 | fh2ResponseMatrixCombined->SetName("fh2ResponseMatrixCombined"); | |
969 | ||
970 | fhEfficiencyCombined = (TH1D*)fh2ResponseMatrixCombined->ProjectionX("fhEfficiencyCombined"); | |
971 | fhEfficiencyCombined->Reset(); | |
972 | ||
973 | for(int i=1; i<=fh2ResponseMatrixCombined->GetNbinsX(); i++) { | |
974 | Double_t sumyield = 0.; | |
975 | for(int j=1; j<=fh2ResponseMatrixCombined->GetYaxis()->GetNbins(); j++) { | |
976 | sumyield+=fh2ResponseMatrixCombined->GetBinContent(i,j); | |
977 | } | |
978 | Double_t kCounter = 0.; | |
979 | Double_t kPtLowEdge = fh2ResponseMatrixCombined->GetXaxis()->GetBinLowEdge(i); | |
980 | Double_t kPtUpEdge = fh2ResponseMatrixCombined->GetXaxis()->GetBinUpEdge(i); | |
981 | Double_t kJetEffDet = 0.; | |
982 | ||
983 | for(int k=1; k<=fhEfficiencyDet->GetNbinsX(); k++) { | |
984 | if(fhEfficiencyDet->GetXaxis()->GetBinLowEdge(k)>=kPtLowEdge && fhEfficiencyDet->GetXaxis()->GetBinUpEdge(k)<=kPtUpEdge) { | |
985 | kJetEffDet+=fhEfficiencyDet->GetBinContent(k); | |
986 | kCounter+=1.; | |
987 | } | |
988 | } | |
989 | kJetEffDet = kJetEffDet/kCounter; | |
990 | ||
991 | if(kJetEffDet==0.) kJetEffDet=1.; | |
992 | ||
993 | fhEfficiencyCombined->SetBinContent(i,sumyield*kJetEffDet); | |
994 | ||
995 | } | |
996 | ||
997 | } | |
998 | ||
275d3481 | 999 | //-------------------------------------------------------------------------------------------------------------------------------------------------- |
1000 | TH2* AliAnaChargedJetResponseMaker::MultiplityResponseMatrices(TH2 *h2RMDeltaPt, TH2 *h2RMDetector) { | |
1001 | ||
1002 | // | |
1003 | // Make combined response matrix (background flucutuations + detector effects) | |
1004 | // | |
1005 | // hRMDeltaPt is the response matrix for background fluctuations from \delta(p_t) measurement | |
1006 | // 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 | |
1007 | // | |
1008 | // Function assumes that generated/unfolded axis is x-axis and reconstructed is on y-axis on both respone matrices | |
1009 | ||
1010 | ||
03372fd1 | 1011 | TH2D *h2ResponseMatrixCombined = (TH2D*)h2RMDeltaPt->Clone("h2ResponseMatrixCombined"); //h2RMDeltaPt is the bkg fluctuations RM which has the dimensions we want for the combined response matrix |
275d3481 | 1012 | h2ResponseMatrixCombined->SetTitle("h2ResponseMatrixCombined"); |
1013 | h2ResponseMatrixCombined->SetName("h2ResponseMatrixCombined"); | |
1014 | ||
1015 | // M = RM_deltaPt * RM_DetEffects * T (M=measured T=truth) | |
1016 | // Dimensions: | |
1017 | // mx1 = mxd * dxt * tx1 | |
1018 | // m = measured bins | |
1019 | // t = truth bins | |
1020 | // d = rec for RM_DetEffects and gen for RM_deltaPt | |
1021 | // RM_comb = RM_deltaPt * RM_DetEffects (dimensions mxt) | |
1022 | // RM_comb(m,t) = Sum_d RM_deltaPt(m,d)*RM_DetEffects(d,t) | |
1023 | ||
1024 | if(fDebug) { | |
03372fd1 | 1025 | printf("Nt=%d\n",h2ResponseMatrixCombined->GetNbinsX()); |
1026 | printf("Nm=%d\n",h2ResponseMatrixCombined->GetNbinsY()); | |
1027 | printf("Nd=%d\n",h2RMDeltaPt->GetNbinsX()); | |
275d3481 | 1028 | } |
1029 | ||
03372fd1 | 1030 | |
275d3481 | 1031 | for(Int_t t=1; t<=h2ResponseMatrixCombined->GetNbinsX();t++) { |
1032 | for(Int_t m=1; m<=h2ResponseMatrixCombined->GetNbinsY();m++) { | |
1033 | Double_t valueSum = 0.; | |
1034 | for(Int_t d=1; d<=h2RMDeltaPt->GetNbinsX();d++) { | |
1035 | valueSum += h2RMDeltaPt->GetBinContent(d,m) * h2RMDetector->GetBinContent(t,d); | |
03372fd1 | 1036 | // if(t==10 && m==10) cout << "sum m,d=" << m << "," << d << endl; |
275d3481 | 1037 | }//d-loop |
03372fd1 | 1038 | // if(t==10) cout << "t,m = " << t << "," << m << "\tvalueSum: " << valueSum << endl; |
275d3481 | 1039 | h2ResponseMatrixCombined->SetBinContent(t,m,valueSum); |
1040 | } //m-loop | |
1041 | }//t-loop | |
1042 | ||
1043 | return h2ResponseMatrixCombined; | |
1044 | ||
1045 | } | |
31d3e4f0 | 1046 | |
5d87a047 | 1047 | //-------------------------------------------------------------------------------------------------------------------------------------------------- |
1048 | TH2* AliAnaChargedJetResponseMaker::GetTransposeResponsMatrix(TH2 *h2RM) { | |
1049 | // | |
1050 | // Transpose response matrix | |
1051 | // | |
1052 | ||
1053 | //Initialize transposed response matrix h2RMTranspose | |
1054 | TArrayD *arrayX = (TArrayD*)h2RM->GetXaxis()->GetXbins(); | |
1055 | for(int i=0; i<arrayX->GetSize(); i++) cout << "i: " << arrayX->At(i) << endl; | |
1056 | Double_t *xbinsArray = arrayX->GetArray(); | |
1057 | ||
1058 | TArrayD *arrayY = (TArrayD*)h2RM->GetYaxis()->GetXbins(); | |
1059 | for(int i=0; i<arrayY->GetSize(); i++) cout << "i: " << arrayY->At(i) << endl; | |
1060 | Double_t *ybinsArray = arrayY->GetArray(); | |
1061 | ||
1062 | TH2D *h2RMTranspose = new TH2D("h2RMTranspose","h2RMTranspose",h2RM->GetNbinsY(),ybinsArray,h2RM->GetNbinsX(),xbinsArray); | |
1063 | ||
1064 | //Fill tranposed response matrix | |
1065 | for (Int_t ibin = 1; ibin <= h2RMTranspose->GetNbinsX(); ibin++) { | |
1066 | for (Int_t jbin = 1; jbin <= h2RMTranspose->GetNbinsY(); jbin++) { | |
1067 | h2RMTranspose->SetBinContent(ibin,jbin, h2RM->GetBinContent(jbin,ibin)); | |
1068 | } | |
1069 | } | |
1070 | ||
1071 | ||
1072 | return h2RMTranspose; | |
1073 | ||
1074 | } | |
1075 | ||
1076 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
1077 | TH2* AliAnaChargedJetResponseMaker::NormalizeResponsMatrixYaxisWithPrior(TH2 *h2RM, TH1 *hPrior) { | |
1078 | // | |
1079 | // Normalize such that the Y projection is the prior | |
1080 | // | |
1081 | ||
03372fd1 | 1082 | // TH1D *hProjRespY = (TH1D*)h2RM->ProjectionY("hProjRespY"); |
5d87a047 | 1083 | double intPrior = hPrior->Integral();//"width"); |
1084 | for (Int_t jbin = 1; jbin <= h2RM->GetNbinsY(); jbin++) { | |
1085 | // double corr = hPrior->GetBinContent(jbin)/hProjRespY->GetBinContent(jbin); | |
1086 | for (Int_t ibin = 1; ibin <= h2RM->GetNbinsX(); ibin++) { | |
1087 | double content = h2RM->GetBinContent(ibin,jbin); | |
1088 | // h2RM->SetBinContent(ibin,jbin,content*corr); | |
1089 | h2RM->SetBinContent(ibin,jbin,hPrior->GetBinContent(jbin)/hPrior->GetBinWidth(jbin)/intPrior*content); | |
1090 | } | |
1091 | } | |
1092 | ||
1093 | return h2RM; | |
1094 | } | |
1095 | ||
31d3e4f0 | 1096 | //-------------------------------------------------------------------------------------------------------------------------------------------------- |
1097 | TH1D* AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TH1 *hGen, TH2 *hResponse,TH1 *hEfficiency,Bool_t bDrawSlices) { | |
1098 | ||
1099 | // | |
1100 | // Multiply hGen with response matrix to obtain refolded spectrum | |
ef62323a | 1101 | // Efficiency must be given. In case efficiency is 1 use SetFlatEfficiency(1.) before calling function |
31d3e4f0 | 1102 | // |
1103 | ||
ef62323a | 1104 | if(!hEfficiency) { |
5d87a047 | 1105 | // printf("Efficiency must be given. In case efficiency is 1 use SetFlatEfficiency(1.) before calling function. Aborting!"); |
1106 | // return 0; | |
1107 | printf("Setting efficiency to 1 \n"); | |
1108 | hEfficiency = (TH1D*)hGen->Clone("hEfficiency"); | |
1109 | hEfficiency->Reset(); | |
1110 | for(int i=1; i<=hEfficiency->GetNbinsX(); i++) hEfficiency->SetBinContent(i,1.); | |
ef62323a | 1111 | } |
1112 | ||
31d3e4f0 | 1113 | //For response |
1114 | //x-axis: generated | |
1115 | //y-axis: reconstructed | |
1116 | if(fDebug>0) cout << "nbins hGen: " << hGen->GetNbinsX() << "\t nbins hResponseGen: " << hResponse->GetXaxis()->GetNbins() << "\t nbins hResponseRec: " << hResponse->GetYaxis()->GetNbins() << endl; | |
1117 | ||
1118 | TH1D *hRec = hResponse->ProjectionY("hRec"); | |
1119 | hRec->Sumw2(); | |
1120 | hRec->Reset(); | |
1121 | hRec->SetTitle("hRec"); | |
1122 | hRec->SetName("hRec"); | |
1123 | ||
1124 | for(int irec=1; irec<=hRec->GetNbinsX(); irec++) | |
1125 | hRec->SetBinContent(irec,0); | |
1126 | ||
1127 | TH1D *hRecGenBin = 0x0; | |
1128 | TCanvas *cSlices = 0x0; | |
1129 | if(bDrawSlices) { | |
1130 | cSlices = new TCanvas("cSlices","cSlices:Slices",400,400); | |
1131 | gPad->SetLogy(); | |
1132 | } | |
1133 | ||
1134 | Double_t yieldMC = 0.; | |
1135 | Double_t yieldMCerror = 0.; | |
1136 | Double_t sumYield = 0.; | |
1137 | const Int_t nbinsRec = hRec->GetNbinsX(); | |
1138 | Double_t sumError2[nbinsRec+1]; | |
0bf6381d | 1139 | for(int irec=0; irec<=hRec->GetNbinsX(); irec++) |
1140 | sumError2[irec]=0.; | |
31d3e4f0 | 1141 | Double_t eff = 0.; |
1142 | ||
1143 | for(int igen=1; igen<=hGen->GetNbinsX(); igen++) { | |
1144 | //get pTMC | |
1145 | sumYield = 0.; | |
1146 | if(fEffFlat>0.) | |
1147 | eff = fEffFlat; | |
1148 | else | |
1149 | eff = hEfficiency->GetBinContent(igen); | |
1150 | yieldMC = hGen->GetBinContent(igen)*eff; | |
adf3920d | 1151 | yieldMCerror = hGen->GetBinError(igen)*eff; |
31d3e4f0 | 1152 | |
1153 | if(bDrawSlices) { | |
1154 | hRecGenBin = hResponse->ProjectionY(Form("hRecGenBin%d",igen)); | |
1155 | hRecGenBin->Sumw2(); | |
1156 | hRecGenBin->Reset(); | |
1157 | hRecGenBin->SetTitle(Form("hRecGenBin%d",igen)); | |
1158 | hRecGenBin->SetName(Form("hRecGenBin%d",igen)); | |
1159 | } | |
1160 | ||
1161 | for(int irec=1; irec<=hRec->GetNbinsX(); irec++) { | |
1162 | hRec->AddBinContent(irec,yieldMC*hResponse->GetBinContent(igen,irec)); | |
1163 | sumYield+=hResponse->GetBinContent(igen,irec); | |
346cfe4a | 1164 | // Double_t A = yieldMC*hResponse->GetBinError(igen,irec); |
31d3e4f0 | 1165 | Double_t B = hResponse->GetBinContent(igen,irec)*yieldMCerror; |
346cfe4a | 1166 | // Double_t tmp2 = A*A + B*B; |
adf3920d | 1167 | //sumError2[irec-1] += tmp2 ; |
1168 | sumError2[irec-1] += B*B; | |
31d3e4f0 | 1169 | |
1170 | if(bDrawSlices) | |
1171 | hRecGenBin->SetBinContent(irec,yieldMC*hResponse->GetBinContent(igen,irec)); | |
1172 | ||
1173 | } | |
1174 | if(bDrawSlices) { | |
1175 | cSlices->cd(); | |
1176 | hRecGenBin->SetLineColor(igen); | |
1177 | if(igen==1) hRecGenBin->DrawCopy(); | |
1178 | else hRecGenBin->DrawCopy("same"); | |
1179 | } | |
1180 | ||
1181 | if(hRecGenBin) delete hRecGenBin; | |
1182 | ||
1183 | cout << "igen: " << igen << "\tpTMC: " << hGen->GetXaxis()->GetBinCenter(igen) << "\teff:" << eff << "\tsumYield: " << sumYield << endl; | |
1184 | } | |
1185 | ||
ef62323a | 1186 | for(int i=0; i<=nbinsRec; i++) { |
1187 | if(sumError2[i]>0.) | |
1188 | hRec->SetBinError(i+1,TMath::Sqrt(sumError2[i])); | |
1189 | } | |
31d3e4f0 | 1190 | |
1191 | ||
1192 | return hRec; | |
1193 | } | |
1194 | ||
1195 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
1196 | TH1D* AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TF1 *fGen, TH2 *hResponse,TH1 *hEfficiency) { | |
1197 | ||
ef62323a | 1198 | // |
1199 | // Multiply fGen function with response matrix to obtain (re)folded spectrum | |
1200 | // Efficiency must be given. In case efficiency is 1 use SetFlatEfficiency(1.) before calling function | |
1201 | // | |
1202 | ||
31d3e4f0 | 1203 | //For response |
1204 | //x-axis: generated | |
1205 | //y-axis: reconstructed | |
1206 | ||
ada4df5e | 1207 | if(fDebug>0) printf(">>AliAnaChargedJetResponseMaker::MultiplyResponseGenerated(TF1 *fGen, TH2 *hResponse,TH1 *hEfficiency)\n"); |
ef62323a | 1208 | |
31d3e4f0 | 1209 | TH1D *hRec = hResponse->ProjectionY("hRec"); |
1210 | hRec->Sumw2(); | |
1211 | hRec->Reset(); | |
1212 | hRec->SetTitle("hRec"); | |
1213 | hRec->SetName("hRec"); | |
1214 | ||
1215 | // TH1D *hRec = new TH1D("hRec","hRec",hResponse->GetNbinsY(),hResponse->GetYaxis()->GetXmin(),hResponse->GetYaxis()->GetXmax()); | |
1216 | ||
1217 | for(int irec=1; irec<=hRec->GetNbinsX(); irec++) | |
1218 | hRec->SetBinContent(irec,0); | |
1219 | ||
1220 | Double_t yieldMC = 0.; | |
1221 | Double_t sumYield = 0.; | |
1222 | Double_t eff = 0.; | |
1223 | for(int igen=1; igen<=hResponse->GetNbinsX(); igen++) { | |
1224 | //get pTMC | |
1225 | sumYield = 0.; | |
1226 | double pTMC = hResponse->GetXaxis()->GetBinCenter(igen); | |
ada4df5e | 1227 | if(hEfficiency) { |
1228 | int binEff = hEfficiency->FindBin(pTMC); | |
31d3e4f0 | 1229 | eff = hEfficiency->GetBinContent(binEff); |
ada4df5e | 1230 | } |
1231 | else eff = 1.; | |
1232 | // yieldMC = fGen->Eval(pTMC)*eff; | |
1233 | yieldMC = fGen->Integral(hResponse->GetXaxis()->GetBinLowEdge(igen),hResponse->GetXaxis()->GetBinUpEdge(igen))*eff; | |
31d3e4f0 | 1234 | for(int irec=1; irec<=hResponse->GetNbinsY(); irec++) { |
1235 | hRec->AddBinContent(irec,yieldMC*hResponse->GetBinContent(igen,irec)); | |
1236 | sumYield+=hResponse->GetBinContent(igen,irec); | |
1237 | } | |
ada4df5e | 1238 | // cout << "igen: " << igen << "\tpTMC: " << pTMC << "\tsumYield: " << sumYield << endl; |
1239 | // cout << "yieldMC: " << yieldMC << endl; | |
1240 | ||
31d3e4f0 | 1241 | } |
1242 | ||
1243 | return hRec; | |
1244 | } | |
1245 | ||
0bf6381d | 1246 | //-------------------------------------------------------------------------------------------------------------------------------------------------- |
1247 | /* static */ Double_t AliAnaChargedJetResponseMaker::GetBetaPerDOFValue(Int_t betaColl, Int_t betaOpt) { | |
1248 | ||
1249 | Float_t betaPerDOFoptions[6] = {0.}; | |
1250 | //define beta per degree of freedom (DOF=nbins in measured) | |
1251 | if(betaColl==0) { | |
1252 | betaPerDOFoptions[0] = 9.1e-9; | |
1253 | betaPerDOFoptions[1] = 2.25e-8; | |
1254 | betaPerDOFoptions[2] = 4.55e-8; | |
1255 | betaPerDOFoptions[3] = 9.1e-8; | |
1256 | betaPerDOFoptions[4] = 2.25e-7; | |
1257 | betaPerDOFoptions[5] = 4.55e-7; | |
1258 | } | |
1259 | if(betaColl==1) { | |
1260 | betaPerDOFoptions[0] = 9.1e-7; | |
1261 | betaPerDOFoptions[1] = 2.25e-6; | |
1262 | betaPerDOFoptions[2] = 4.55e-6; | |
1263 | betaPerDOFoptions[3] = 9.1e-6; | |
1264 | betaPerDOFoptions[4] = 2.25e-5; | |
1265 | betaPerDOFoptions[5] = 4.55e-5; | |
1266 | } | |
1267 | if(betaColl==2) { | |
1268 | betaPerDOFoptions[0] = 9.1e-5; | |
1269 | betaPerDOFoptions[1] = 2.25e-4; | |
1270 | betaPerDOFoptions[2] = 4.55e-4; | |
1271 | betaPerDOFoptions[3] = 9.1e-4; | |
1272 | betaPerDOFoptions[4] = 2.25e-3; | |
1273 | betaPerDOFoptions[5] = 4.55e-3; | |
1274 | } | |
1275 | if(betaColl==3) { | |
1276 | betaPerDOFoptions[0] = 9.1e-3; | |
1277 | betaPerDOFoptions[1] = 2.25e-2; | |
1278 | betaPerDOFoptions[2] = 4.55e-2; | |
1279 | betaPerDOFoptions[3] = 9.1e-2; | |
1280 | betaPerDOFoptions[4] = 2.25e-1; | |
1281 | betaPerDOFoptions[5] = 4.55e-1; | |
1282 | } | |
1283 | if(betaColl==4) { | |
1284 | betaPerDOFoptions[0] = 9.1e-1; | |
1285 | betaPerDOFoptions[1] = 2.25; | |
1286 | betaPerDOFoptions[2] = 4.55; | |
1287 | betaPerDOFoptions[3] = 9.1; | |
1288 | betaPerDOFoptions[4] = 22.5; | |
1289 | betaPerDOFoptions[5] = 45.5; | |
1290 | } | |
1291 | ||
1292 | return betaPerDOFoptions[betaOpt]; | |
1293 | ||
1294 | } | |
1295 | ||
31d3e4f0 | 1296 | //-------------------------------------------------------------------------------------------------------------------------------------------------- |
1297 | Double_t AliAnaChargedJetResponseMaker::InterpolateFast(TGraph *gr, Double_t x) { | |
1298 | ||
1299 | Double_t ipmax = gr->GetN()-1.; | |
1300 | Double_t x2,y2,xold,yold; | |
1301 | ||
1302 | Double_t xmin,ymin,xmax, ymax; | |
1303 | gr->GetPoint(0,xmin,ymin); | |
1304 | gr->GetPoint(gr->GetN()-1,xmax,ymax); | |
1305 | if(x<xmin || x>xmax) return 0; | |
1306 | ||
1307 | Double_t ip = ipmax/2.; | |
1308 | Double_t ipold = 0.; | |
1309 | gr->GetPoint((int)(ip),x2,y2); | |
1310 | ||
1311 | Int_t i = 0; | |
1312 | ||
1313 | if(x2>x) { | |
1314 | while(x2>x) { | |
1315 | if(i==0) ipold=0.; | |
1316 | ip -= (ip)/2.; | |
1317 | gr->GetPoint((int)(ip),x2,y2); | |
1318 | if(x2>x){} | |
1319 | else ipold = ip; | |
1320 | i++; | |
1321 | // cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl; | |
1322 | } | |
1323 | } | |
1324 | else if(x2<x) { | |
1325 | while(x2<x) { | |
1326 | if(i==0) ipold=ipmax; | |
1327 | ip += (ipold-ip)/2.; | |
1328 | gr->GetPoint((int)(ip),x2,y2); | |
1329 | if(x2>x) ipold = ip; | |
1330 | else {} | |
1331 | i++; | |
1332 | // cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl; | |
1333 | } | |
1334 | } | |
1335 | ||
1336 | int ip2 = 0; | |
1337 | if(x2>x) { | |
1338 | ip2 = (int)(ip)-1; | |
1339 | gr->GetPoint(ip2,x2,y2); | |
1340 | while(x2>x) { | |
1341 | ip2--; | |
1342 | gr->GetPoint(ip2,x2,y2); | |
1343 | } | |
1344 | gr->GetPoint(ip2+1,xold,yold); | |
1345 | } | |
1346 | else { | |
1347 | ip2 = (int)(ip)+1; | |
1348 | gr->GetPoint(ip2,x2,y2); | |
1349 | while(x2<x) { | |
1350 | ip2++; | |
1351 | gr->GetPoint(ip2,x2,y2); | |
1352 | } | |
1353 | gr->GetPoint(ip2-1,xold,yold); | |
1354 | ||
1355 | } | |
1356 | // cout << "For x=" << x << " interpolate between: " << xold << " and " << x2 << endl; | |
1357 | return ((x-xold)*y2 + (x2-x)*yold) / (x2-xold); | |
1358 | ||
1359 | } | |
1360 | ||
1361 | //-------------------------------------------------------------------------------------------------------------------------------------------------- | |
1362 | Double_t AliAnaChargedJetResponseMaker::InterpolateFast(TH1 *h, Double_t x) { | |
1363 | ||
1364 | // Double_t ipmax = gr->GetN()-1.; | |
1365 | Double_t ipmax = (double)h->GetNbinsX(); | |
1366 | Double_t x2,y2,xold,yold; | |
1367 | ||
1368 | Double_t xmin = h->GetXaxis()->GetXmin(); | |
1369 | Double_t xmax = h->GetXaxis()->GetXmax(); | |
1370 | if(x<xmin || x>xmax) return 0; | |
1371 | ||
1372 | Double_t ip = ipmax/2.; | |
1373 | Double_t ipold = 0.; | |
1374 | ||
1375 | x2 = h->GetBinCenter((int)ip); | |
1376 | y2 = h->GetBinContent((int)ip); | |
1377 | ||
1378 | Int_t i = 0; | |
1379 | ||
1380 | if(x2>x) { | |
1381 | while(x2>x) { | |
1382 | if(i==0) ipold=0.; | |
1383 | ip -= (ip)/2.; | |
1384 | x2 = h->GetBinCenter((int)ip); | |
1385 | y2 = h->GetBinContent((int)ip); | |
1386 | if(x2>x) {} | |
1387 | else ipold = ip; | |
1388 | i++; | |
1389 | // cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl; | |
1390 | } | |
1391 | } | |
1392 | else if(x2<x) { | |
1393 | while(x2<x) { | |
1394 | if(i==0) ipold=ipmax; | |
1395 | ip += (ipold-ip)/2.; | |
1396 | x2 = h->GetBinCenter((int)ip); | |
1397 | y2 = h->GetBinContent((int)ip); | |
1398 | if(x2>x) ipold = ip; | |
1399 | else {} | |
1400 | i++; | |
1401 | // cout << "ipold: " << ipold << "\tip: " << ip << "\tx2: " << x2 << "\ty2: " << y2 << endl; | |
1402 | } | |
1403 | } | |
1404 | ||
1405 | int ip2 = 0; | |
1406 | if(x2>x) { | |
1407 | ip2 = (int)(ip)-1; | |
1408 | x2 = h->GetBinCenter(ip2); | |
1409 | y2 = h->GetBinContent(ip2); | |
1410 | while(x2>x) { | |
1411 | ip2--; | |
1412 | x2 = h->GetBinCenter(ip2); | |
1413 | y2 = h->GetBinContent(ip2); | |
1414 | } | |
1415 | xold = h->GetBinCenter(ip2+1); | |
1416 | yold = h->GetBinContent(ip2+1); | |
1417 | } | |
1418 | else { | |
1419 | ip2 = (int)(ip)+1; | |
1420 | x2 = h->GetBinCenter(ip2); | |
1421 | y2 = h->GetBinContent(ip2); | |
1422 | while(x2<x) { | |
1423 | ip2++; | |
1424 | x2 = h->GetBinCenter(ip2); | |
1425 | y2 = h->GetBinContent(ip2); | |
1426 | } | |
1427 | xold = h->GetBinCenter(ip2-1); | |
1428 | yold = h->GetBinContent(ip2-1); | |
1429 | ||
1430 | } | |
1431 | // cout << "For x=" << x << " interpolate between: " << xold << " and " << x2 << endl; | |
1432 | return ((x-xold)*y2 + (x2-x)*yold) / (x2-xold); | |
1433 | ||
1434 | } | |
1435 | ||
1436 |