Fixes for coverity.
[u/mrichter/AliRoot.git] / PWGJE / AliFragmentationFunctionCorrections.cxx
CommitLineData
39e2b057 1// *************************************************************************
2// * *
3// * corrections to Fragmentation Functions *
4// * *
5// *************************************************************************
6
7
8/**************************************************************************
9 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
10 * *
11 * Author: The ALICE Off-line Project. *
12 * Contributors are mentioned in the code where appropriate. *
13 * *
14 * Permission to use, copy, modify and distribute this software and its *
15 * documentation strictly for non-commercial purposes is hereby granted *
16 * without fee, provided that the above copyright notice appears in all *
17 * copies and that both the copyright notice and this permission notice *
18 * appear in the supporting documentation. The authors make no claims *
19 * about the suitability of this software for any purpose. It is *
20 * provided "as is" without express or implied warranty. *
21 **************************************************************************/
22
23/* $Id: */
24
25#include "TMath.h"
26#include "TH2F.h"
27#include "THnSparse.h"
28#include "TFile.h"
29#include "TDirectory.h"
30#include "AliCFUnfolding.h"
31#include "AliFragmentationFunctionCorrections.h"
ce55b926 32#include <iostream> // OB TEST!!!
33#include <fstream>
c2aad3ae 34using std::cout;
35using std::endl;
36using std::cerr;
37using std::fstream;
39e2b057 38
39///////////////////////////////////////////////////////////////////////////////
40// //
41// correction class for ouput of AliAnalysisTaskFragmentationFunction //
42// //
43// corrections for: reconstruction efficiency, momentum resolution, //
44// secondaries, UE / HI background, fluctuations //
45// back-correction on jet energy on dN/dxi //
46// //
47// read MC ouput and write out efficiency histos, response matrices etc. //
48// read measured distributions and apply efficiency, response matrices etc. //
49// //
50// contact: o.busch@gsi.de //
51// //
52///////////////////////////////////////////////////////////////////////////////
53
54
55ClassImp(AliFragmentationFunctionCorrections)
56
57//________________________________________________________________________
58AliFragmentationFunctionCorrections::AliFragmentationFunctionCorrections()
59 : TObject()
60 ,fDebug(0)
61 ,fNJetPtSlices(0)
62 ,fJetPtSlices(0)
63 ,fNJets(0)
64 ,fNJetsBgr(0)
65 ,fNHistoBinsSinglePt(0)
66 ,fHistoBinsSinglePt(0)
67 ,fNHistoBinsPt(0)
68 ,fNHistoBinsZ(0)
69 ,fNHistoBinsXi(0)
70 ,fHistoBinsPt(0)
71 ,fHistoBinsZ(0)
72 ,fHistoBinsXi(0)
73 ,fNCorrectionLevels(0)
74 ,fCorrFF(0)
75 ,fNCorrectionLevelsBgr(0)
76 ,fCorrBgr(0)
77 ,fNCorrectionLevelsSinglePt(0)
78 ,fCorrSinglePt(0)
39e2b057 79 ,fh1EffSinglePt(0)
80 ,fh1EffPt(0)
81 ,fh1EffZ(0)
82 ,fh1EffXi(0)
83 ,fh1EffBgrPt(0)
84 ,fh1EffBgrZ(0)
85 ,fh1EffBgrXi(0)
ce55b926 86 ,fh1BbBCorrSinglePt(0)
87 ,fh1BbBPt(0)
88 ,fh1BbBZ(0)
89 ,fh1BbBXi(0)
90 ,fh1BbBBgrPt(0)
91 ,fh1BbBBgrZ(0)
92 ,fh1BbBBgrXi(0)
93 ,fh1FoldingCorrPt(0)
94 ,fh1FoldingCorrZ(0)
95 ,fh1FoldingCorrXi(0)
96 ,fh1SecCorrPt(0)
97 ,fh1SecCorrZ(0)
98 ,fh1SecCorrXi(0)
99 ,fh1SecCorrBgrPt(0)
100 ,fh1SecCorrBgrZ(0)
101 ,fh1SecCorrBgrXi(0)
39e2b057 102 ,fh1FFTrackPtBackFolded(0)
103 ,fh1FFZBackFolded(0)
104 ,fh1FFXiBackFolded(0)
105 ,fh1FFRatioTrackPtFolded(0)
106 ,fh1FFRatioZFolded(0)
107 ,fh1FFRatioXiFolded(0)
108 ,fh1FFRatioTrackPtBackFolded(0)
109 ,fh1FFRatioZBackFolded(0)
110 ,fh1FFRatioXiBackFolded(0)
111 ,fh1SingleTrackPtBackFolded(0)
112 ,fh1RatioSingleTrackPtFolded(0)
113 ,fh1RatioSingleTrackPtBackFolded(0)
114 ,fhnResponseSinglePt(0)
115 ,fhnResponsePt(0)
116 ,fhnResponseZ(0)
117 ,fhnResponseXi(0)
118 ,fh1FFTrackPtPrior(0)
119 ,fh1FFZPrior(0)
120 ,fh1FFXiPrior(0)
121 ,fh1SecCorrSinglePt(0)
122{
123 // default constructor
124}
125
126//________________________________________________________________________________________________________________________
127AliFragmentationFunctionCorrections::AliFragmentationFunctionCorrections(const AliFragmentationFunctionCorrections &copy)
128 : TObject()
129 ,fDebug(copy.fDebug)
130 ,fNJetPtSlices(copy.fNJetPtSlices)
131 ,fJetPtSlices(copy.fJetPtSlices)
132 ,fNJets(copy.fNJets)
133 ,fNJetsBgr(copy.fNJetsBgr)
1aa4f09f 134 ,fNHistoBinsSinglePt(copy.fNHistoBinsSinglePt)
135 ,fHistoBinsSinglePt(copy.fHistoBinsSinglePt)
39e2b057 136 ,fNHistoBinsPt(copy.fNHistoBinsPt)
137 ,fNHistoBinsZ(copy.fNHistoBinsZ)
138 ,fNHistoBinsXi(copy.fNHistoBinsXi)
139 ,fHistoBinsPt(copy.fHistoBinsPt)
140 ,fHistoBinsZ(copy.fHistoBinsZ)
141 ,fHistoBinsXi(copy.fHistoBinsXi)
142 ,fNCorrectionLevels(copy.fNCorrectionLevels)
143 ,fCorrFF(copy.fCorrFF)
144 ,fNCorrectionLevelsBgr(copy.fNCorrectionLevelsBgr)
145 ,fCorrBgr(copy.fCorrBgr)
146 ,fNCorrectionLevelsSinglePt(copy.fNCorrectionLevelsSinglePt)
147 ,fCorrSinglePt(copy.fCorrSinglePt)
1aa4f09f 148 ,fh1EffSinglePt(copy.fh1EffSinglePt)
39e2b057 149 ,fh1EffPt(copy.fh1EffPt)
150 ,fh1EffZ(copy.fh1EffZ)
151 ,fh1EffXi(copy.fh1EffXi)
152 ,fh1EffBgrPt(copy.fh1EffBgrPt)
153 ,fh1EffBgrZ(copy.fh1EffBgrZ)
154 ,fh1EffBgrXi(copy.fh1EffBgrXi)
ce55b926 155 ,fh1BbBCorrSinglePt(copy.fh1BbBCorrSinglePt)
156 ,fh1BbBPt(copy.fh1BbBPt)
157 ,fh1BbBZ(copy.fh1BbBZ)
158 ,fh1BbBXi(copy.fh1BbBXi)
159 ,fh1BbBBgrPt(copy.fh1BbBBgrPt)
160 ,fh1BbBBgrZ(copy.fh1BbBBgrZ)
161 ,fh1BbBBgrXi(copy.fh1BbBBgrXi)
162 ,fh1FoldingCorrPt(copy.fh1FoldingCorrPt)
163 ,fh1FoldingCorrZ(copy.fh1FoldingCorrZ)
164 ,fh1FoldingCorrXi(copy.fh1FoldingCorrXi)
165 ,fh1SecCorrPt(copy.fh1SecCorrPt)
166 ,fh1SecCorrZ(copy.fh1SecCorrZ)
167 ,fh1SecCorrXi(copy.fh1SecCorrXi)
168 ,fh1SecCorrBgrPt(copy.fh1SecCorrBgrPt)
169 ,fh1SecCorrBgrZ(copy.fh1SecCorrBgrZ)
170 ,fh1SecCorrBgrXi(copy.fh1SecCorrBgrXi)
39e2b057 171 ,fh1FFTrackPtBackFolded(copy.fh1FFTrackPtBackFolded)
172 ,fh1FFZBackFolded(copy.fh1FFZBackFolded)
173 ,fh1FFXiBackFolded(copy.fh1FFXiBackFolded)
174 ,fh1FFRatioTrackPtFolded(copy.fh1FFRatioTrackPtFolded)
175 ,fh1FFRatioZFolded(copy.fh1FFRatioZFolded)
176 ,fh1FFRatioXiFolded(copy.fh1FFRatioXiFolded)
177 ,fh1FFRatioTrackPtBackFolded(copy.fh1FFRatioTrackPtBackFolded)
178 ,fh1FFRatioZBackFolded(copy.fh1FFRatioZBackFolded)
179 ,fh1FFRatioXiBackFolded(copy.fh1FFRatioXiBackFolded)
180 ,fh1SingleTrackPtBackFolded(copy.fh1SingleTrackPtBackFolded)
181 ,fh1RatioSingleTrackPtFolded(copy.fh1RatioSingleTrackPtFolded)
182 ,fh1RatioSingleTrackPtBackFolded(copy.fh1RatioSingleTrackPtBackFolded)
183 ,fhnResponseSinglePt(copy.fhnResponseSinglePt)
184 ,fhnResponsePt(copy.fhnResponsePt)
185 ,fhnResponseZ(copy.fhnResponseZ)
186 ,fhnResponseXi(copy.fhnResponseXi)
187 ,fh1FFTrackPtPrior(copy.fh1FFTrackPtPrior)
188 ,fh1FFZPrior(copy.fh1FFZPrior)
189 ,fh1FFXiPrior(copy.fh1FFXiPrior)
190 ,fh1SecCorrSinglePt(copy.fh1SecCorrSinglePt)
191{
192 // copy constructor
193
194}
195
196// ______________________________________________________________________________________________________________________________
197AliFragmentationFunctionCorrections& AliFragmentationFunctionCorrections::operator=(const AliFragmentationFunctionCorrections& o)
198{
199 // assignment
200
201 if(this!=&o){
202 TObject::operator=(o);
203 fDebug = o.fDebug;
204 fNJetPtSlices = o.fNJetPtSlices;
205 fJetPtSlices = o.fJetPtSlices;
206 fNJets = o.fNJets;
207 fNJetsBgr = o.fNJetsBgr;
208 fNHistoBinsSinglePt = o.fNHistoBinsSinglePt;
209 fHistoBinsSinglePt = o.fHistoBinsSinglePt;
210 fNHistoBinsPt = o.fNHistoBinsPt;
211 fNHistoBinsZ = o.fNHistoBinsZ;
212 fNHistoBinsXi = o.fNHistoBinsXi;
213 fHistoBinsPt = o.fHistoBinsPt;
214 fHistoBinsZ = o.fHistoBinsZ;
215 fHistoBinsXi = o.fHistoBinsXi;
216 fNCorrectionLevels = o.fNCorrectionLevels;
217 fCorrFF = o.fCorrFF;
218 fNCorrectionLevelsBgr = o.fNCorrectionLevelsBgr;
219 fCorrBgr = o.fCorrBgr;
220 fNCorrectionLevelsSinglePt = o.fNCorrectionLevelsSinglePt;
221 fCorrSinglePt = o.fCorrSinglePt;
39e2b057 222 fh1EffSinglePt = o.fh1EffSinglePt;
223 fh1EffPt = o.fh1EffPt;
224 fh1EffZ = o.fh1EffZ;
225 fh1EffXi = o.fh1EffXi;
226 fh1EffBgrPt = o.fh1EffBgrPt;
227 fh1EffBgrZ = o.fh1EffBgrZ;
228 fh1EffBgrXi = o.fh1EffBgrXi;
ce55b926 229 fh1BbBCorrSinglePt = o.fh1BbBCorrSinglePt;
230 fh1BbBPt = o.fh1BbBPt;
231 fh1BbBZ = o.fh1BbBZ;
232 fh1BbBXi = o.fh1BbBXi;
233 fh1BbBBgrPt = o.fh1BbBBgrPt;
234 fh1BbBBgrZ = o.fh1BbBBgrZ;
235 fh1BbBBgrXi = o.fh1BbBBgrXi;
236 fh1FoldingCorrPt = o.fh1FoldingCorrPt;
237 fh1FoldingCorrZ = o.fh1FoldingCorrZ;
238 fh1FoldingCorrXi = o.fh1FoldingCorrXi;
239 fh1SecCorrPt = o.fh1SecCorrPt;
240 fh1SecCorrZ = o.fh1SecCorrZ;
241 fh1SecCorrXi = o.fh1SecCorrXi;
242 fh1SecCorrBgrPt = o.fh1SecCorrBgrPt;
243 fh1SecCorrBgrZ = o.fh1SecCorrBgrZ;
244 fh1SecCorrBgrXi = o.fh1SecCorrBgrXi;
39e2b057 245 fh1FFTrackPtBackFolded = o.fh1FFTrackPtBackFolded;
246 fh1FFZBackFolded = o.fh1FFZBackFolded;
247 fh1FFXiBackFolded = o.fh1FFXiBackFolded;
248 fh1FFRatioTrackPtFolded = o.fh1FFRatioTrackPtFolded;
249 fh1FFRatioZFolded = o.fh1FFRatioZFolded;
250 fh1FFRatioXiFolded = o.fh1FFRatioXiFolded;
251 fh1FFRatioTrackPtBackFolded = o.fh1FFRatioTrackPtBackFolded;
252 fh1FFRatioZBackFolded = o.fh1FFRatioZBackFolded;
253 fh1FFRatioXiBackFolded = o.fh1FFRatioXiBackFolded;
254 fh1SingleTrackPtBackFolded = o.fh1SingleTrackPtBackFolded;
255 fh1RatioSingleTrackPtFolded = o.fh1RatioSingleTrackPtFolded;
256 fh1RatioSingleTrackPtBackFolded = o.fh1RatioSingleTrackPtBackFolded;
257 fhnResponseSinglePt = o.fhnResponseSinglePt;
258 fhnResponsePt = o.fhnResponsePt;
259 fhnResponseZ = o.fhnResponseZ;
260 fhnResponseXi = o.fhnResponseXi;
261 fh1FFTrackPtPrior = o.fh1FFTrackPtPrior;
262 fh1FFZPrior = o.fh1FFZPrior;
263 fh1FFXiPrior = o.fh1FFXiPrior;
264 fh1SecCorrSinglePt = o.fh1SecCorrSinglePt;
265 }
266
267 return *this;
268}
269
270//_________________________________________________________________________
271AliFragmentationFunctionCorrections::~AliFragmentationFunctionCorrections()
272{
273 // destructor
274
275 if(fJetPtSlices) delete fJetPtSlices;
276 if(fNJets) delete fNJets;
277 if(fNJetsBgr) delete fNJetsBgr;
278
39e2b057 279 DeleteHistoArray(fh1EffPt);
280 DeleteHistoArray(fh1EffXi);
281 DeleteHistoArray(fh1EffZ );
282
283 DeleteHistoArray(fh1EffBgrPt);
284 DeleteHistoArray(fh1EffBgrXi);
285 DeleteHistoArray(fh1EffBgrZ);
286
ce55b926 287 // BbB
288 DeleteHistoArray(fh1BbBPt);
289 DeleteHistoArray(fh1BbBXi);
290 DeleteHistoArray(fh1BbBZ);
291
292 DeleteHistoArray(fh1BbBBgrPt);
293 DeleteHistoArray(fh1BbBBgrXi);
294 DeleteHistoArray(fh1BbBBgrZ);
295
296 DeleteHistoArray(fh1FoldingCorrPt);
297 DeleteHistoArray(fh1FoldingCorrXi);
298 DeleteHistoArray(fh1FoldingCorrZ);
299
300 // sec corr
301 DeleteHistoArray(fh1SecCorrPt);
302 DeleteHistoArray(fh1SecCorrXi);
303 DeleteHistoArray(fh1SecCorrZ);
304
305 DeleteHistoArray(fh1SecCorrBgrPt);
306 DeleteHistoArray(fh1SecCorrBgrXi);
307 DeleteHistoArray(fh1SecCorrBgrZ);
308
39e2b057 309 // unfolding
310
311 DeleteHistoArray(fh1FFTrackPtBackFolded);
312 DeleteHistoArray(fh1FFZBackFolded);
313 DeleteHistoArray(fh1FFXiBackFolded);
314
315 DeleteHistoArray(fh1FFRatioTrackPtFolded);
316 DeleteHistoArray(fh1FFRatioZFolded);
317 DeleteHistoArray(fh1FFRatioXiFolded);
318
319 DeleteHistoArray(fh1FFRatioTrackPtBackFolded);
320 DeleteHistoArray(fh1FFRatioZBackFolded);
321 DeleteHistoArray(fh1FFRatioXiBackFolded);
322
323 DeleteTHnSparseArray(fhnResponsePt);
324 DeleteTHnSparseArray(fhnResponseZ);
325 DeleteTHnSparseArray(fhnResponseXi);
326
327 DeleteHistoArray(fh1FFTrackPtPrior);
328 DeleteHistoArray(fh1FFZPrior);
329 DeleteHistoArray(fh1FFXiPrior);
330
39e2b057 331 // clean up corrected FF
332
333 for(Int_t c=0; c<fNCorrectionLevels; c++) delete fCorrFF[c];
334 delete[] fCorrFF;
335
336 // clean up bgr
337
338 for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) delete fCorrBgr[c];
339 delete[] fCorrBgr;
340
341 // clean up inclusive pt
342 for(Int_t c=0; c<fNCorrectionLevelsSinglePt; c++) delete fCorrSinglePt[c];
343 delete[] fCorrSinglePt;
344
345 delete[] fNHistoBinsPt;
346 delete[] fNHistoBinsZ;
347 delete[] fNHistoBinsXi;
348
349 // clean up histo bins
350
351 if(fHistoBinsSinglePt) delete fHistoBinsSinglePt;
352
353 for(Int_t i=0; i<fNJetPtSlices; i++) delete fHistoBinsPt[i];
354 for(Int_t i=0; i<fNJetPtSlices; i++) delete fHistoBinsZ[i];
355 for(Int_t i=0; i<fNJetPtSlices; i++) delete fHistoBinsXi[i];
356
357 delete[] fHistoBinsPt;
358 delete[] fHistoBinsZ;
359 delete[] fHistoBinsXi;
360}
361
362//_________________________________________________________________________________
363AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::AliFragFuncCorrHistos()
364 : TObject()
365 ,fArraySize(0)
366 ,fh1CorrFFTrackPt(0)
367 ,fh1CorrFFZ(0)
368 ,fh1CorrFFXi(0)
369 ,fCorrLabel(0)
370{
371 // default constructor
372
373}
374
375//__________________________________________________________________________________________________________________
6daac008 376AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::AliFragFuncCorrHistos(const AliFragmentationFunctionCorrections::AliFragFuncCorrHistos& copy)
39e2b057 377 : TObject()
378 ,fArraySize(copy.fArraySize)
6daac008 379 ,fh1CorrFFTrackPt(0)
380 ,fh1CorrFFZ(0)
381 ,fh1CorrFFXi(0)
39e2b057 382 ,fCorrLabel(copy.fCorrLabel)
383{
384 // copy constructor
6daac008 385
386 fh1CorrFFTrackPt = new TH1F*[copy.fArraySize];
387 fh1CorrFFZ = new TH1F*[copy.fArraySize];
388 fh1CorrFFXi = new TH1F*[copy.fArraySize];
389
1aa4f09f 390 for(Int_t i=0; i<copy.fArraySize; i++){
6daac008 391 fh1CorrFFTrackPt[i] = new TH1F(*(copy.fh1CorrFFTrackPt[i]));
392 fh1CorrFFZ[i] = new TH1F(*(copy.fh1CorrFFZ[i]));
393 fh1CorrFFXi[i] = new TH1F(*(copy.fh1CorrFFXi[i]));
1aa4f09f 394 }
39e2b057 395}
396
39e2b057 397//_______________________________________________________________________________________________________________________________________________________________
398AliFragmentationFunctionCorrections::AliFragFuncCorrHistos& AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::operator=(const AliFragmentationFunctionCorrections::AliFragFuncCorrHistos& o)
399{
400 // assignment
401
402 if(this!=&o){
403 TObject::operator=(o);
6daac008 404 Int_t fArraySize_new = o.fArraySize;
39e2b057 405 fCorrLabel = o.fCorrLabel;
1aa4f09f 406
6daac008 407 TH1F** fh1CorrFFTrackPt_new = new TH1F*[fArraySize_new];
408 TH1F** fh1CorrFFZ_new = new TH1F*[fArraySize_new];
409 TH1F** fh1CorrFFXi_new = new TH1F*[fArraySize_new];
410
411 for(Int_t i=0; i<fArraySize_new; i++){
412 fh1CorrFFTrackPt_new[i] = new TH1F(*(o.fh1CorrFFTrackPt[i]));
413 fh1CorrFFZ_new[i] = new TH1F(*(o.fh1CorrFFZ[i]));
414 fh1CorrFFXi_new[i] = new TH1F(*(o.fh1CorrFFXi[i]));
1aa4f09f 415 }
39e2b057 416
6daac008 417 // ---
418
419 if(fArraySize){
420 for(int i=0; i<fArraySize; i++) delete fh1CorrFFTrackPt[i];
421 for(int i=0; i<fArraySize; i++) delete fh1CorrFFZ[i];
422 for(int i=0; i<fArraySize; i++) delete fh1CorrFFXi[i];
423 }
424
425 if(fh1CorrFFTrackPt) delete[] fh1CorrFFTrackPt;
426 if(fh1CorrFFZ) delete[] fh1CorrFFZ;
427 if(fh1CorrFFXi) delete[] fh1CorrFFXi;
428
429 // ---
430
431 fArraySize = fArraySize_new;
432
433 fh1CorrFFTrackPt = fh1CorrFFTrackPt_new;
434 fh1CorrFFZ = fh1CorrFFZ_new;
435 fh1CorrFFXi = fh1CorrFFXi_new;
436
437 for(Int_t i=0; i<fArraySize; i++){
438 fh1CorrFFTrackPt[i] = fh1CorrFFTrackPt_new[i];
439 fh1CorrFFZ[i] = fh1CorrFFZ_new[i];
440 fh1CorrFFXi[i] = fh1CorrFFXi_new[i];
441 }
442 }
443
39e2b057 444 return *this;
445}
446
447//__________________________________________________________________________________
448AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::~AliFragFuncCorrHistos()
449{
450 // destructor
451
452 if(fArraySize){
453 for(int i=0; i<fArraySize; i++) delete fh1CorrFFTrackPt[i];
454 for(int i=0; i<fArraySize; i++) delete fh1CorrFFZ[i];
455 for(int i=0; i<fArraySize; i++) delete fh1CorrFFXi[i];
456 }
457
458 if(fh1CorrFFTrackPt) delete[] fh1CorrFFTrackPt;
459 if(fh1CorrFFZ) delete[] fh1CorrFFZ;
460 if(fh1CorrFFXi) delete[] fh1CorrFFXi;
461
462}
463
464//___________________________________________________________________________________________________________________
465AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::AliFragFuncCorrHistos(const char* label, Int_t arraySize)
466 : TObject()
467 ,fArraySize(0)
468 ,fh1CorrFFTrackPt(0)
469 ,fh1CorrFFZ(0)
470 ,fh1CorrFFXi(0)
471 ,fCorrLabel(label)
472{
473 // constructor
474
475 fArraySize = arraySize;
476 fh1CorrFFTrackPt = new TH1F*[arraySize];
477 fh1CorrFFZ = new TH1F*[arraySize];
478 fh1CorrFFXi = new TH1F*[arraySize];
479
480 for(int i=0; i<arraySize; i++) fh1CorrFFTrackPt[i] = new TH1F;
481 for(int i=0; i<arraySize; i++) fh1CorrFFZ[i] = new TH1F;
482 for(int i=0; i<arraySize; i++) fh1CorrFFXi[i] = new TH1F;
483}
484
485//_______________________________________________________________________________________________________________________________
486void AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::AddCorrHistos(Int_t slice,TH1F* histPt,TH1F* histZ,TH1F* histXi)
487{
488 // place histo in array, add corrLabel to histo name
489
490 if(slice>=fArraySize){
491 Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
492 return;
493 }
494
495 if(histPt){
496 TString name = histPt->GetName();
497 if(fCorrLabel.Length()>0) name += "_"+fCorrLabel;
498 histPt->SetName(name);
499
500 if(!(histPt->GetSumw2())) histPt->Sumw2();
501
502 new(fh1CorrFFTrackPt[slice]) TH1F(*histPt);
503 }
504
505 if(histZ){
506 TString name = histZ->GetName();
507 if(fCorrLabel.Length()>0) name += "_"+fCorrLabel;
508 histZ->SetName(name);
509
510 if(!(histZ->GetSumw2())) histZ->Sumw2();
511
512 new(fh1CorrFFZ[slice]) TH1F(*histZ);
513 }
514
515 if(histXi){
516 TString name = histXi->GetName();
517 if(fCorrLabel.Length()>0) name += "_"+fCorrLabel;
518 histXi->SetName(name);
519
520 if(!(histXi->GetSumw2())) histXi->Sumw2();
521
522 new(fh1CorrFFXi[slice]) TH1F(*histXi);
523 }
524}
525
526//___________________________________________________________________________________________________________________________________
527void AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::ReplaceCorrHistos(Int_t slice,TH1F* histPt,TH1F* histZ,TH1F* histXi)
528{
529 // replace histo in array
530
531 if(slice>=fArraySize){
532 Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
533 return;
534 }
535
536 if(histPt){
537 if(!(histPt->GetSumw2())) histPt->Sumw2();
538
539 TString name = histPt->GetName();
540 histPt->SetName(name);
541
542 delete fh1CorrFFTrackPt[slice];
543 fh1CorrFFTrackPt[slice] = new TH1F;
544 new(fh1CorrFFTrackPt[slice]) TH1F(*histPt);
545 }
546
547 if(histZ){
548 if(!(histZ->GetSumw2())) histZ->Sumw2();
549
550 TString name = histZ->GetName();
551 histZ->SetName(name);
552
553 delete fh1CorrFFZ[slice];
554 fh1CorrFFZ[slice] = new TH1F;
555 new(fh1CorrFFZ[slice]) TH1F(*histZ);
556 }
557
558 if(histXi){
559 if(!(histXi->GetSumw2())) histXi->Sumw2();
560
561 TString name = histXi->GetName();
562 histXi->SetName(name);
563
564 delete fh1CorrFFXi[slice];
565 fh1CorrFFXi[slice] = new TH1F;
566 new(fh1CorrFFXi[slice]) TH1F(*histXi);
567 }
568}
569
570// ___________________________________________________________________________________________
571TH1F* AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::GetTrackPt(const Int_t slice)
572{
573 // return pt histo
574
575 if(slice>=fArraySize){
576 Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
577 return 0;
578 }
579
580 return fh1CorrFFTrackPt[slice];
581
582}
583
584// ______________________________________________________________________________________
585TH1F* AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::GetZ(const Int_t slice)
586{
587 // return z histo
588
589 if(slice>=fArraySize){
590 Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
591 return 0;
592 }
593
594 return fh1CorrFFZ[slice];
595}
596
597// ________________________________________________________________________________________
598TH1F* AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::GetXi(const Int_t slice)
599{
600 // return xi histo
601
602 if(slice>=fArraySize){
603 Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
604 return 0;
605 }
606
607 return fh1CorrFFXi[slice];
608}
609
610// __________________________________________________________________________
611void AliFragmentationFunctionCorrections::DeleteHistoArray(TH1F** hist) const
612{
613 // delete array of TH1
614
615 if(hist){
616 for(Int_t i=0; i<fNJetPtSlices; i++) delete hist[i];
617 delete[] hist;
618 }
619}
620
621// ____________________________________________________________________________________
622void AliFragmentationFunctionCorrections::DeleteTHnSparseArray(THnSparse** hist) const
623{
624 // delete array of THnSparse
625
626 if(hist){
627 for(Int_t i=0; i<fNJetPtSlices; i++) delete hist[i];
628 delete[] hist;
629 }
630}
631
632// _________________________________________________________
633TH1F** AliFragmentationFunctionCorrections::BookHistoArray()
634{
635 // book array of TH1
636
637 if(!fNJetPtSlices){
638 Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
639 return 0;
640 }
641
642 TH1F** hist = new TH1F*[fNJetPtSlices];
643 for(Int_t i=0; i<fNJetPtSlices; i++) hist[i] = new TH1F();
644
645 return hist;
646}
647
648//__________________________________________________________________
649THnSparse** AliFragmentationFunctionCorrections::BookTHnSparseArray()
650{
651 // book array of THnSparse
652
653 if(!fNJetPtSlices){
654 Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
655 return 0;
656 }
657
658 THnSparse** hist = new THnSparse*[fNJetPtSlices];
659 for(Int_t i=0; i<fNJetPtSlices; i++) hist[i] = new THnSparseF();
660
661 return hist;
662}
663
664//_____________________________________________________________________________
665void AliFragmentationFunctionCorrections::AddCorrectionLevel(const char* label)
666{
667 // increase corr level
668
669 if(!fNJetPtSlices){
670 Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
671 return;
672 }
673
674 if(fNCorrectionLevels >= fgMaxNCorrectionLevels){
675 Printf("%s:%d -- max correction level exceeded", (char*)__FILE__,__LINE__);
676 return;
677 }
678
679 fCorrFF[fNCorrectionLevels] = new AliFragFuncCorrHistos(label,fNJetPtSlices);
680 fNCorrectionLevels++;
681}
682
683
684//________________________________________________________________________________
685void AliFragmentationFunctionCorrections::AddCorrectionLevelBgr(const char* label)
686{
687 // increase corr level for bgr FF
688
689 if(!fNJetPtSlices){
690 if(fDebug>0) Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
691 return;
692 }
693
694 if(fNCorrectionLevelsBgr >= fgMaxNCorrectionLevels){
695 Printf("%s:%d -- max correction level exceeded", (char*)__FILE__,__LINE__);
696 return;
697 }
698
699 fCorrBgr[fNCorrectionLevelsBgr] = new AliFragFuncCorrHistos(label,fNJetPtSlices);
700 fNCorrectionLevelsBgr++;
701}
702
703//_____________________________________________________________________________________
704void AliFragmentationFunctionCorrections::AddCorrectionLevelSinglePt(const char* label)
705{
706 // increase corr level single pt spec
707
708 Int_t nJetPtSlicesSingle = 1; // pro forma
709
710 if(fNCorrectionLevelsSinglePt >= fgMaxNCorrectionLevels){
711 Printf("%s:%d -- max correction level exceeded", (char*)__FILE__,__LINE__);
712 return;
713 }
714
715 fCorrSinglePt[fNCorrectionLevelsSinglePt] = new AliFragFuncCorrHistos(label,nJetPtSlicesSingle);
716 fNCorrectionLevelsSinglePt++;
717}
718
719//_____________________________________________________________________________________________
720void AliFragmentationFunctionCorrections::SetJetPtSlices(Float_t* bins, const Int_t nJetPtSlices)
721{
722 // set jet pt slices
723 // once slices are known, can book all other histos
724
725 fNJetPtSlices = nJetPtSlices;
726
727 const Int_t size = nJetPtSlices+1;
728 fJetPtSlices = new TArrayF(size,bins);
729
730 // nJets array
731
732 fNJets = new TArrayF(size);
733 fNJets->Reset(0);
734
735 fNJetsBgr = new TArrayF(size);
736 fNJetsBgr->Reset(0);
737
738 // histos
739
740 fNCorrectionLevels = 0;
741 fCorrFF = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels];
742 AddCorrectionLevel(); // first 'correction' level = raw FF
743
744 fNCorrectionLevelsBgr = 0;
745 fCorrBgr = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels];
746 AddCorrectionLevelBgr(); // first 'correction' level = raw bgr dist
747
39e2b057 748 // eff histos
749
750 fh1EffPt = BookHistoArray();
751 fh1EffXi = BookHistoArray();
752 fh1EffZ = BookHistoArray();
753
754 fh1EffBgrPt = BookHistoArray();
755 fh1EffBgrXi = BookHistoArray();
756 fh1EffBgrZ = BookHistoArray();
757
758
ce55b926 759 // BbB
760
761 fh1BbBPt = BookHistoArray();
762 fh1BbBXi = BookHistoArray();
763 fh1BbBZ = BookHistoArray();
764
765 fh1BbBBgrPt = BookHistoArray();
766 fh1BbBBgrXi = BookHistoArray();
767 fh1BbBBgrZ = BookHistoArray();
768
769 fh1FoldingCorrPt = BookHistoArray();
770 fh1FoldingCorrXi = BookHistoArray();
771 fh1FoldingCorrZ = BookHistoArray();
772
773 // SecCorr
774
775 fh1SecCorrPt = BookHistoArray();
776 fh1SecCorrXi = BookHistoArray();
777 fh1SecCorrZ = BookHistoArray();
778
779 fh1SecCorrBgrPt = BookHistoArray();
780 fh1SecCorrBgrXi = BookHistoArray();
781 fh1SecCorrBgrZ = BookHistoArray();
782
39e2b057 783 // unfolding
784
785 fh1FFTrackPtBackFolded = BookHistoArray();
786 fh1FFXiBackFolded = BookHistoArray();
787 fh1FFZBackFolded = BookHistoArray();
788
789 fh1FFRatioTrackPtFolded = BookHistoArray();
790 fh1FFRatioXiFolded = BookHistoArray();
791 fh1FFRatioZFolded = BookHistoArray();
792
793 fh1FFRatioTrackPtBackFolded = BookHistoArray();
794 fh1FFRatioXiBackFolded = BookHistoArray();
795 fh1FFRatioZBackFolded = BookHistoArray();
796
797 //
798
799 fhnResponsePt = BookTHnSparseArray();
800 fhnResponseZ = BookTHnSparseArray();
801 fhnResponseXi = BookTHnSparseArray();
802
803 fh1FFTrackPtPrior = BookHistoArray();
804 fh1FFZPrior = BookHistoArray();
805 fh1FFXiPrior = BookHistoArray();
806
39e2b057 807 // histos bins
808
809 fNHistoBinsPt = new Int_t[fNJetPtSlices];
810 fNHistoBinsXi = new Int_t[fNJetPtSlices];
811 fNHistoBinsZ = new Int_t[fNJetPtSlices];
812
813 for(Int_t i=0; i<fNJetPtSlices; i++) fNHistoBinsPt[i] = 0;
814 for(Int_t i=0; i<fNJetPtSlices; i++) fNHistoBinsXi[i] = 0;
815 for(Int_t i=0; i<fNJetPtSlices; i++) fNHistoBinsZ[i] = 0;
816
817 fHistoBinsPt = new TArrayD*[fNJetPtSlices];
818 fHistoBinsXi = new TArrayD*[fNJetPtSlices];
819 fHistoBinsZ = new TArrayD*[fNJetPtSlices];
820
821 for(Int_t i=0; i<fNJetPtSlices; i++) fHistoBinsPt[i] = new TArrayD(0);
822 for(Int_t i=0; i<fNJetPtSlices; i++) fHistoBinsXi[i] = new TArrayD(0);
823 for(Int_t i=0; i<fNJetPtSlices; i++) fHistoBinsZ[i] = new TArrayD(0);
824}
825
826//_____________________________________________________________________________________________________________________________________
827void AliFragmentationFunctionCorrections::SetHistoBins(const Int_t jetPtSlice, const Int_t sizeBins, Double_t* bins, const Int_t type)
828{
829 // set histo bins for jet pt slice
830 // if binning undefined for any slice, original binning will be used
831
832 if(!fNJetPtSlices){
833 Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
834 return;
835 }
836
837 if(jetPtSlice>=fNJetPtSlices){
838 Printf("%s:%d -- jetPtSlice %d exceeds max",(char*)__FILE__,__LINE__,jetPtSlice);
839 return;
840 }
841
842 if(type == kFlagPt){
843 fNHistoBinsPt[jetPtSlice] = sizeBins-1;
844 fHistoBinsPt[jetPtSlice]->Set(sizeBins,bins);
845 }
846 else if(type==kFlagZ){
847 fNHistoBinsZ[jetPtSlice] = sizeBins-1;
848 fHistoBinsZ[jetPtSlice]->Set(sizeBins,bins);
849 }
850
851 else if(type==kFlagXi){
852 fNHistoBinsXi[jetPtSlice] = sizeBins-1;
853 fHistoBinsXi[jetPtSlice]->Set(sizeBins,bins);
854 }
855}
856
857//__________________________________________________________________________________________________________________________________________________________________
858void AliFragmentationFunctionCorrections::SetHistoBins(const Int_t jetPtSlice, const Int_t nBinsLimits, Double_t* binsLimits, Double_t* binsWidth, const Int_t type)
859{
860 // set histo bins for jet pt slice
861 // function takes array of limits and widths (e.g. 1 GeV bins up to 10 GeV, 2 GeV width up to 20 GeV, ...)
862 // array size of binsLimits: nBinsLimits
863 // array size of binsWidth: nBinsLimits-1
864 // binsLimits have to be in increasing order
865 // if binning undefined for any slice, original binning will be used
866
867 if(!fNJetPtSlices){
868 Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
869 return;
870 }
871
872 if(jetPtSlice>=fNJetPtSlices){
873 Printf("%s:%d -- jetPtSlice %d exceeds max",(char*)__FILE__,__LINE__,jetPtSlice);
874 return;
875 }
876
877
878 Double_t binLimitMin = binsLimits[0];
879 Double_t binLimitMax = binsLimits[nBinsLimits-1];
880
881 Double_t binLimit = binLimitMin; // start value
882
883 Int_t sizeUpperLim = 10000; //static_cast<Int_t>(binLimitMax/binsWidth[0])+1;
884 TArrayD binsArray(sizeUpperLim);
885 Int_t nBins = 0;
886 binsArray.SetAt(binLimitMin,nBins++);
887
ce55b926 888 while(binLimit<0.999999*binLimitMax && nBins<sizeUpperLim){ // 0.999 numerical safety factor
39e2b057 889
890 Int_t currentSlice = -1;
891 for(Int_t i=0; i<nBinsLimits; i++){
ce55b926 892 if(binLimit >= 0.999999*binsLimits[i]) currentSlice = i;
39e2b057 893 }
894
895 Double_t currentBinWidth = binsWidth[currentSlice];
896 binLimit += currentBinWidth;
897
ce55b926 898 // if(type == kFlagZ) std::cout<<" SetHistoBins: nBins "<<nBins<<" binLimit "<<binLimit<<" binLimitMax "<<binLimitMax
899 // <<" currentSlice "<<currentSlice<<std::endl;
900
39e2b057 901 binsArray.SetAt(binLimit,nBins++);
902 }
903
904 Double_t* bins = binsArray.GetArray();
905
906 SetHistoBins(jetPtSlice,nBins,bins,type);
907}
908
ce55b926 909//_____________________________________________________________________________________________________________________________________
80acde99 910TArrayD* AliFragmentationFunctionCorrections::GetHistoBins(const Int_t& jetPtSlice, const Int_t& type)
ce55b926 911{
912 // set histo bins for jet pt slice
913 // if binning undefined for any slice, original binning will be used
914
915 if(!fNJetPtSlices){
916 Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
917 return 0;
918 }
919
920 if(jetPtSlice>=fNJetPtSlices){
921 Printf("%s:%d -- jetPtSlice %d exceeds max",(char*)__FILE__,__LINE__,jetPtSlice);
922 return 0;
923 }
924
925 if(type == kFlagPt){
926 return fHistoBinsPt[jetPtSlice];
927 }
928 else if(type==kFlagZ){
929 return fHistoBinsZ[jetPtSlice];
930 }
931
932 else if(type==kFlagXi){
933 return fHistoBinsXi[jetPtSlice];
934 }
935 else{
936 Printf("%s%d unknown type",(char*)__FILE__,__LINE__);
937 return 0;
938 }
939}
940
39e2b057 941//__________________________________________________________________________________________________
942void AliFragmentationFunctionCorrections::SetHistoBinsSinglePt(const Int_t sizeBins, Double_t* bins)
943{
944 // set histo bins for inclusive pt spectra
945 // if binning undefined, original binning will be used
946
947 fNHistoBinsSinglePt = sizeBins-1;
948
949 fHistoBinsSinglePt = new TArrayD(sizeBins);
950 fHistoBinsSinglePt->Set(sizeBins,bins);
951}
952
953//__________________________________________________________________________________________________________________________________________________________________
954void AliFragmentationFunctionCorrections::SetHistoBinsSinglePt(const Int_t nBinsLimits, Double_t* binsLimits, Double_t* binsWidth)
955{
956 // set histo bins for inclusive pt spectra
957 // function takes array of limits and widths (e.g. 1 GeV bins up to 10 GeV, 2 GeV width up to 20 GeV, ...)
958 // array size of binsLimits: nBinsLimits
959 // array size of binsWidth: nBinsLimits-1
960 // binsLimits have to be in increasing order
961 // if binning undefined for any slice, original binning will be used
962
963
964 Double_t binLimitMin = binsLimits[0];
965 Double_t binLimitMax = binsLimits[nBinsLimits-1];
966
967 Double_t binLimit = binLimitMin; // start value
968
969 Int_t sizeUpperLim = 10000; //static_cast<Int_t>(binLimitMax/binsWidth[0])+1;
970 TArrayD binsArray(sizeUpperLim);
971 Int_t nBins = 0;
972 binsArray.SetAt(binLimitMin,nBins++);
973
974 while(binLimit<binLimitMax && nBins<sizeUpperLim){
975
976 Int_t currentSlice = -1;
977 for(Int_t i=0; i<nBinsLimits; i++){
978 if(binLimit >= 0.999*binsLimits[i]) currentSlice = i; // 0.999 numerical saftey factor
979 }
980
981 Double_t currentBinWidth = binsWidth[currentSlice];
982 binLimit += currentBinWidth;
983
984 binsArray.SetAt(binLimit,nBins++);
985 }
986
987 Double_t* bins = binsArray.GetArray();
988
989 SetHistoBinsSinglePt(nBins,bins);
990}
991
992//____________________________________________________________________________________
993void AliFragmentationFunctionCorrections::NormalizeTH1(TH1* hist, const Float_t nJets)
994{
995 // FF normalization: divide by bin width and normalize to entries/jets
996 // should also work for TH2, but to be tested !!!
997
998 if(nJets == 0){ // nothing to do
999 if(fDebug>0) Printf("%s:%d -- normalize: nJets = 0, do nothing", (char*)__FILE__,__LINE__);
1000 return;
1001 }
1002
1003 Int_t nBins = hist->GetNbinsX()*hist->GetNbinsY()*hist->GetNbinsZ();
1004
1005 for(Int_t bin=0; bin<=nBins; bin++){ // include under-/overflow (?)
1006
1007 Double_t binWidth = hist->GetBinWidth(bin);
1008 Double_t binCont = hist->GetBinContent(bin);
1009 Double_t binErr = hist->GetBinError(bin);
1010
1011 binCont /= binWidth;
1012 binErr /= binWidth;
1013
1014 hist->SetBinContent(bin,binCont);
1015 hist->SetBinError(bin,binErr);
1016 }
1017
1018 hist->Scale(1/nJets);
1019}
1020
1021//_____________________________________________________
1022void AliFragmentationFunctionCorrections::NormalizeFF()
1023{
1024 // normalize FF
1025
1026 if(fNCorrectionLevels>1){
1027 Printf("%s:%d -- FF normalization should be done BEFORE any correction !!!", (char*)__FILE__,__LINE__);
1028 }
1029
1030 for(Int_t i=0; i<fNJetPtSlices; i++){
1031
1032 if(fDebug>0) Printf(" normalizeFF: i %d, nJets %f",i,fNJets->At(i));
1033
1034 NormalizeTH1(fCorrFF[0]->GetTrackPt(i),fNJets->At(i)); // always normalize corr level 0 = raw FF
1035 NormalizeTH1(fCorrFF[0]->GetZ(i),fNJets->At(i));
1036 NormalizeTH1(fCorrFF[0]->GetXi(i),fNJets->At(i));
1037 }
1038}
1039
1040//______________________________________________________
1041void AliFragmentationFunctionCorrections::NormalizeBgr()
1042{
1043 // normalize bgr/UE FF
1044
1045 if(fNCorrectionLevelsBgr>1){
1046 Printf("%s:%d -- FF normalization should be done BEFORE any correction !!!", (char*)__FILE__,__LINE__);
1047 }
1048
1049 for(Int_t i=0; i<fNJetPtSlices; i++){
1050 NormalizeTH1(fCorrBgr[0]->GetTrackPt(i), fNJetsBgr->At(i)); // always normalize corr level 0 = raw FF
1051 NormalizeTH1(fCorrBgr[0]->GetZ(i), fNJetsBgr->At(i));
1052 NormalizeTH1(fCorrBgr[0]->GetXi(i),fNJetsBgr->At(i));
1053 }
1054
1055}
1056
1057//__________________________________________________________________________________________________
1058void AliFragmentationFunctionCorrections::ReadRawFF(TString strfile, TString strID, TString strFFID)
1059{
1060 // read raw FF - standard dir/list name
1061
b541fbca 1062 TString strdir = "PWGJE_FragmentationFunction_" + strID;
39e2b057 1063 TString strlist = "fracfunc_" + strID;
1064
1065 ReadRawFF(strfile,strdir,strlist,strFFID);
1066}
1067
1068//____________________________________________________________________________________________________________________
1069void AliFragmentationFunctionCorrections::ReadRawFF(TString strfile, TString strdir, TString strlist, TString strFFID)
1070{
1071 // get raw FF from input file, project in jet pt slice
1072 // normalization done separately
1073
1074 TFile f(strfile,"READ");
1075
1076 if(!f.IsOpen()){
1077 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1078 return;
1079 }
1080
ce55b926 1081 if(fDebug>0) Printf("%s:%d -- read FF from file %s, dir %s ",(char*)__FILE__,__LINE__,strfile.Data(), strdir.Data());
39e2b057 1082
1083 gDirectory->cd(strdir);
1084
1085 TList* list = 0;
1086
ce55b926 1087 if(strlist && strlist.Length()){
1088 if(!(list = (TList*) gDirectory->Get(strlist))){
1089 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1090 return;
1091 }
39e2b057 1092 }
1093
1094 TString hnameJetPt(Form("fh1FFJetPt%s",strFFID.Data()));
1095 TString hnameTrackPt(Form("fh2FFTrackPt%s",strFFID.Data()));
1096 TString hnameZ(Form("fh2FFZ%s",strFFID.Data()));
1097 TString hnameXi(Form("fh2FFXi%s",strFFID.Data()));
1098
ce55b926 1099 TH1F* fh1FFJetPt = 0;
1100 TH2F* fh2FFTrackPt = 0;
1101 TH2F* fh2FFZ = 0;
1102 TH2F* fh2FFXi = 0;
1103
1104 if(list){
1105 fh1FFJetPt = (TH1F*) list->FindObject(hnameJetPt);
1106 fh2FFTrackPt = (TH2F*) list->FindObject(hnameTrackPt);
1107 fh2FFZ = (TH2F*) list->FindObject(hnameZ);
1108 fh2FFXi = (TH2F*) list->FindObject(hnameXi);
1109 }
1110 else{
1111 fh1FFJetPt = (TH1F*) gDirectory->Get(hnameJetPt);
1112 fh2FFTrackPt = (TH2F*) gDirectory->Get(hnameTrackPt);
1113 fh2FFZ = (TH2F*) gDirectory->Get(hnameZ);
1114 fh2FFXi = (TH2F*) gDirectory->Get(hnameXi);
1115 }
1116
39e2b057 1117
1118 if(!fh1FFJetPt) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameJetPt.Data()); return; }
1119 if(!fh2FFTrackPt){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameTrackPt.Data()); return; }
1120 if(!fh2FFZ) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameZ.Data()); return; }
1121 if(!fh2FFXi) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameXi.Data()); return; }
1122
1123 fh1FFJetPt->SetDirectory(0);
1124 fh2FFTrackPt->SetDirectory(0);
1125 fh2FFZ->SetDirectory(0);
1126 fh2FFXi->SetDirectory(0);
1127
1128 f.Close();
1129
1130
1131 // nJets per bin
1132
1133 for(Int_t i=0; i<fNJetPtSlices; i++){
1134
1135 Float_t jetPtLoLim = fJetPtSlices->At(i);
1136 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1137
1138 Int_t binLo = static_cast<Int_t>(fh1FFJetPt->FindBin(jetPtLoLim));
1139 Int_t binUp = static_cast<Int_t>(fh1FFJetPt->FindBin(jetPtUpLim)) - 1;
1140
1141 Float_t nJetsBin = fh1FFJetPt->Integral(binLo,binUp);
1142
1143 fNJets->SetAt(nJetsBin,i);
1144
1145 if(fDebug>0) Printf("jet pt %d to %d: nJets %f",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),fNJets->At(i));
1146 }
1147
1148 // projections: FF
1149
1150 for(Int_t i=0; i<fNJetPtSlices; i++){
1151
1152 Float_t jetPtLoLim = fJetPtSlices->At(i);
1153 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1154
1155 Int_t binLo = static_cast<Int_t>(fh2FFTrackPt->GetXaxis()->FindBin(jetPtLoLim));
1156 Int_t binUp = static_cast<Int_t>(fh2FFTrackPt->GetXaxis()->FindBin(jetPtUpLim))-1;
1157
1158 if(binUp > fh2FFTrackPt->GetNbinsX()){
1159 Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
1160 return;
1161 }
1162
1163 TString strNameFFPt(Form("fh1FFTrackPt%s_%02d_%02d",strFFID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1164 TString strNameFFZ(Form("fh1FFZ%s_%02d_%02d",strFFID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1165 TString strNameFFXi(Form("fh1FFXi%s_%02d_%02d",strFFID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1166
1167 // appendix 'unbinned' to avoid histos with same name after rebinning
1168 TH1F* projPt = (TH1F*) fh2FFTrackPt->ProjectionY(strNameFFPt+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
1169 TH1F* projZ = (TH1F*) fh2FFZ->ProjectionY(strNameFFZ+"_unBinned",binLo,binUp,"o");
1170 TH1F* projXi = (TH1F*) fh2FFXi->ProjectionY(strNameFFXi+"_unBinned",binLo,binUp,"o");
1171
1172 if(fNHistoBinsPt[i]) projPt = (TH1F*) projPt->Rebin(fNHistoBinsPt[i],strNameFFPt,fHistoBinsPt[i]->GetArray());
1173 if(fNHistoBinsZ[i]) projZ = (TH1F*) projZ->Rebin(fNHistoBinsZ[i],strNameFFZ,fHistoBinsZ[i]->GetArray());
1174 if(fNHistoBinsXi[i]) projXi = (TH1F*) projXi->Rebin(fNHistoBinsXi[i],strNameFFXi,fHistoBinsXi[i]->GetArray());
1175
1176 projPt->SetNameTitle(strNameFFPt,"");
1177 projZ->SetNameTitle(strNameFFZ,"");
1178 projXi->SetNameTitle(strNameFFXi,"");
1179
1180 // raw FF = corr level 0
1181 fCorrFF[0]->AddCorrHistos(i,projPt,projZ,projXi);
ce55b926 1182 }
1183
1184
1185 delete fh1FFJetPt;
1186 delete fh2FFTrackPt;
1187 delete fh2FFZ;
1188 delete fh2FFXi;
39e2b057 1189}
1190
1191//_____________________________________________________________________________________________________________________
1192void AliFragmentationFunctionCorrections::ReadRawBgr(TString strfile, TString strID, TString strBgrID, TString strFFID)
1193{
1194 // read raw FF - standard dir/list name
1195
b541fbca 1196 TString strdir = "PWGJE_FragmentationFunction_" + strID;
39e2b057 1197 TString strlist = "fracfunc_" + strID;
1198
1199 ReadRawBgr(strfile,strdir,strlist,strBgrID,strFFID);
1200}
1201
1202//_______________________________________________________________________________________________________________________________________
1203void AliFragmentationFunctionCorrections::ReadRawBgr(TString strfile, TString strdir, TString strlist, TString strBgrID, TString strFFID)
1204{
1205 // get raw FF from input file, project in jet pt slice
1206 // use jet dN/dpt corresponding to strFFID, bgr FF to strBgrID+strFFID
1207 // e.g. "fh1FFJetPtRecCuts", "fh2FFXiBgrPerpRecCuts"
1208 // normalization done separately
1209
1210 TString strID = strBgrID + strFFID;
1211
1212 TFile f(strfile,"READ");
1213
1214 if(!f.IsOpen()){
1215 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1216 return;
1217 }
1218
ce55b926 1219 if(fDebug>0) Printf("%s:%d -- read Bgr %s from file %s, dir %s ",(char*)__FILE__,__LINE__,strBgrID.Data(),strfile.Data(),strdir.Data());
39e2b057 1220
1221 gDirectory->cd(strdir);
1222
1223 TList* list = 0;
1224
ce55b926 1225 if(strlist && strlist.Length()){
1226 if(!(list = (TList*) gDirectory->Get(strlist))){
1227 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1228 return;
1229 }
39e2b057 1230 }
1231
1232 TString hnameNJets = "fh1nRecJetsCuts";
1233 TString hnameJetPt(Form("fh1FFJetPt%s",strFFID.Data())); // not: strID.Data() !!! would not be proper normalization
1234 TString hnameBgrTrackPt(Form("fh2FFTrackPt%s",strID.Data()));
1235 TString hnameBgrZ(Form("fh2FFZ%s",strID.Data()));
1236 TString hnameBgrXi(Form("fh2FFXi%s",strID.Data()));
1237
ce55b926 1238 TH1F* fh1NJets;
1239 TH1F* fh1FFJetPtBgr;
1240 TH2F* fh2FFTrackPtBgr;
1241 TH2F* fh2FFZBgr;
1242 TH2F* fh2FFXiBgr;
1243
1244 if(list){
1245 fh1NJets = (TH1F*) list->FindObject(hnameNJets); // needed for normalization of bgr out of 2 jets
1246 fh1FFJetPtBgr = (TH1F*) list->FindObject(hnameJetPt);
1247 fh2FFTrackPtBgr = (TH2F*) list->FindObject(hnameBgrTrackPt);
1248 fh2FFZBgr = (TH2F*) list->FindObject(hnameBgrZ);
1249 fh2FFXiBgr = (TH2F*) list->FindObject(hnameBgrXi);
1250 }
1251 else{
1252 fh1NJets = (TH1F*) gDirectory->Get(hnameNJets); // needed for normalization of bgr out of 2 jets
1253 fh1FFJetPtBgr = (TH1F*) gDirectory->Get(hnameJetPt);
1254 fh2FFTrackPtBgr = (TH2F*) gDirectory->Get(hnameBgrTrackPt);
1255 fh2FFZBgr = (TH2F*) gDirectory->Get(hnameBgrZ);
1256 fh2FFXiBgr = (TH2F*) gDirectory->Get(hnameBgrXi);
1257 }
1258
39e2b057 1259
1260 if(!fh1FFJetPtBgr) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameJetPt.Data()); return; }
ce55b926 1261 if(!fh1NJets) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameNJets.Data()); }
39e2b057 1262 if(!fh2FFTrackPtBgr){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrTrackPt.Data()); return; }
1263 if(!fh2FFZBgr) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrZ.Data()); return; }
1264 if(!fh2FFXiBgr) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrXi.Data()); return; }
1265
1266 fh1FFJetPtBgr->SetDirectory(0);
39e2b057 1267 fh2FFTrackPtBgr->SetDirectory(0);
1268 fh2FFZBgr->SetDirectory(0);
1269 fh2FFXiBgr->SetDirectory(0);
ce55b926 1270 if(fh1NJets) fh1NJets->SetDirectory(0);
39e2b057 1271
1272 f.Close();
1273
1274 // nJets per bin
1275
1276 for(Int_t i=0; i<fNJetPtSlices; i++){
1277
1278 Float_t jetPtLoLim = fJetPtSlices->At(i);
1279 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1280
1281 Int_t binLo = static_cast<Int_t>(fh1FFJetPtBgr->FindBin(jetPtLoLim));
1282 Int_t binUp = static_cast<Int_t>(fh1FFJetPtBgr->FindBin(jetPtUpLim)) - 1;
1283
1284 Float_t nJetsBin = fh1FFJetPtBgr->Integral(binLo,binUp);
1285 Double_t scaleF = 1;
1286
1287 //if(strBgrID.Contains("Out2Jets")){ // scale by ratio 2 jets events / all events
1288 // scaleF = fh1NJets->Integral(fh1NJets->FindBin(2),fh1NJets->GetNbinsX())
1289 // / fh1NJets->Integral(fh1NJets->FindBin(1),fh1NJets->GetNbinsX());}
1290
1291
ce55b926 1292 if(strBgrID.Contains("OutAllJets") ){ // scale by ratio >3 jets events / all events
1293
1294 if(fh1NJets){
1295 scaleF = fh1NJets->Integral(fh1NJets->FindBin(4),fh1NJets->GetNbinsX())
1296 / fh1NJets->Integral(fh1NJets->FindBin(1),fh1NJets->GetNbinsX());
1297 }
1298 else{
1299 Printf("%s:%d -- use bgr OutAllJets, but histo fhn1NJets not found",(char*)__FILE__,__LINE__);
1300 }
39e2b057 1301 }
ce55b926 1302
1303
39e2b057 1304 fNJetsBgr->SetAt(nJetsBin*scaleF,i);
1305
1306 if(fDebug>0) Printf("bgr jet pt %d to %d: nJets %f, scaleF %.2f",
1307 static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),nJetsBin,scaleF);
1308
1309 }
1310
1311 // projections: FF
1312
1313 for(Int_t i=0; i<fNJetPtSlices; i++){
1314
1315 Float_t jetPtLoLim = fJetPtSlices->At(i);
1316 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1317
1318 Int_t binLo = static_cast<Int_t>(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtLoLim));
1319 Int_t binUp = static_cast<Int_t>(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtUpLim))-1;
1320
1321 if(binUp > fh2FFTrackPtBgr->GetNbinsX()){
1322 Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
1323 return;
1324 }
1325
1326 TString strNameBgrPt(Form("fh1BgrTrackPt%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1327 TString strNameBgrZ(Form("fh1BgrZ%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1328 TString strNameBgrXi(Form("fh1BgrXi%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1329
1330 // appendix 'unbinned' to avoid histos with same name after rebinning
1331 TH1F* projPt = (TH1F*) fh2FFTrackPtBgr->ProjectionY(strNameBgrPt+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
1332 TH1F* projZ = (TH1F*) fh2FFZBgr->ProjectionY(strNameBgrZ+"_unBinned",binLo,binUp,"o");
1333 TH1F* projXi = (TH1F*) fh2FFXiBgr->ProjectionY(strNameBgrXi+"_unBinned",binLo,binUp,"o");
1334
1335 if(fNHistoBinsPt[i]) projPt = (TH1F*) projPt->Rebin(fNHistoBinsPt[i],strNameBgrPt,fHistoBinsPt[i]->GetArray());
1336 if(fNHistoBinsZ[i]) projZ = (TH1F*) projZ->Rebin(fNHistoBinsZ[i],strNameBgrZ,fHistoBinsZ[i]->GetArray());
1337 if(fNHistoBinsXi[i]) projXi = (TH1F*) projXi->Rebin(fNHistoBinsXi[i],strNameBgrXi,fHistoBinsXi[i]->GetArray());
1338
1339 projPt->SetNameTitle(strNameBgrPt,"");
1340 projZ->SetNameTitle(strNameBgrZ,"");
1341 projXi->SetNameTitle(strNameBgrXi,"");
1342
1343 // raw bgr = corr level 0
1344 fCorrBgr[0]->AddCorrHistos(i,projPt,projZ,projXi);
ce55b926 1345 }
1346
1347 delete fh1FFJetPtBgr;
1348 if(fh1NJets) delete fh1NJets;
1349 delete fh2FFTrackPtBgr;
1350 delete fh2FFZBgr;
1351 delete fh2FFXiBgr;
39e2b057 1352}
1353
1354//_____________________________________________________________________________________________________________________
1355void AliFragmentationFunctionCorrections::ReadRawBgrEmbedding(TString strfile, TString strID, TString strFFID)
1356{
1357 // read raw FF - standard dir/list name
1358
b541fbca 1359 TString strdir = "PWGJE_FragmentationFunction_" + strID;
39e2b057 1360 TString strlist = "fracfunc_" + strID;
1361
1362 ReadRawBgrEmbedding(strfile,strdir,strlist,strFFID);
1363}
1364
1365//_______________________________________________________________________________________________________________________________________
1366void AliFragmentationFunctionCorrections::ReadRawBgrEmbedding(TString strfile, TString strdir, TString strlist, TString strFFID)
1367{
1368 // get raw FF from input file, project in jet pt slice
1369 // for embedding, the bgr FF are taken from histos "fh1FFJetPtRecCuts", "fh2FFXiRecCuts"
1370 // normalization done separately
1371
1372 TString strBgrID = "BckgEmbed";
1373 TString strID = strBgrID + strFFID;
1374
1375 TFile f(strfile,"READ");
1376
1377 if(!f.IsOpen()){
1378 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1379 return;
1380 }
1381
1382 if(fDebug>0) Printf("%s:%d -- read Bgr %s from file %s ",(char*)__FILE__,__LINE__,strFFID.Data(),strfile.Data());
1383
1384 gDirectory->cd(strdir);
1385
1386 TList* list = 0;
1387
1388 if(!(list = (TList*) gDirectory->Get(strlist))){
1389 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1390 return;
1391 }
1392
1393 TString hnameNJets = "fh1nRecJetsCuts";
1394 TString hnameJetPt(Form("fh1FFJetPt%s",strFFID.Data()));
1395 TString hnameBgrTrackPt(Form("fh2FFTrackPt%s",strFFID.Data()));
1396 TString hnameBgrZ(Form("fh2FFZ%s",strFFID.Data()));
1397 TString hnameBgrXi(Form("fh2FFXi%s",strFFID.Data()));
1398
1399 TH1F* fh1NJets = (TH1F*) list->FindObject(hnameNJets); // needed for normalization of bgr out of 2 jets
1400 TH1F* fh1FFJetPtBgr = (TH1F*) list->FindObject(hnameJetPt);
1401 TH2F* fh2FFTrackPtBgr = (TH2F*) list->FindObject(hnameBgrTrackPt);
1402 TH2F* fh2FFZBgr = (TH2F*) list->FindObject(hnameBgrZ);
1403 TH2F* fh2FFXiBgr = (TH2F*) list->FindObject(hnameBgrXi);
1404
1405 if(!fh1FFJetPtBgr) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameJetPt.Data()); return; }
1406 if(!fh1NJets) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameNJets.Data()); return; }
1407 if(!fh2FFTrackPtBgr){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrTrackPt.Data()); return; }
1408 if(!fh2FFZBgr) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrZ.Data()); return; }
1409 if(!fh2FFXiBgr) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrXi.Data()); return; }
1410
1411 fh1FFJetPtBgr->SetDirectory(0);
1412 fh1NJets->SetDirectory(0);
1413 fh2FFTrackPtBgr->SetDirectory(0);
1414 fh2FFZBgr->SetDirectory(0);
1415 fh2FFXiBgr->SetDirectory(0);
1416
1417 f.Close();
1418
1419 // nJets per bin
1420
1421 for(Int_t i=0; i<fNJetPtSlices; i++){
1422
1423 Float_t jetPtLoLim = fJetPtSlices->At(i);
1424 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1425
1426 Int_t binLo = static_cast<Int_t>(fh1FFJetPtBgr->FindBin(jetPtLoLim));
1427 Int_t binUp = static_cast<Int_t>(fh1FFJetPtBgr->FindBin(jetPtUpLim)) - 1;
1428
1429 Float_t nJetsBin = fh1FFJetPtBgr->Integral(binLo,binUp);
1430 Double_t scaleF = 1;
1431
1432 fNJetsBgr->SetAt(nJetsBin*scaleF,i);
1433
1434 if(fDebug>0) Printf("bgr jet pt %d to %d: nJets %f, scaleF %.2f",
1435 static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),nJetsBin,scaleF);
1436
1437 }
1438
1439 // projections: FF
1440
1441 for(Int_t i=0; i<fNJetPtSlices; i++){
1442
1443 Float_t jetPtLoLim = fJetPtSlices->At(i);
1444 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1445
1446 Int_t binLo = static_cast<Int_t>(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtLoLim));
1447 Int_t binUp = static_cast<Int_t>(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtUpLim))-1;
1448
1449 if(binUp > fh2FFTrackPtBgr->GetNbinsX()){
1450 Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
1451 return;
1452 }
1453
1454 TString strNameBgrPt(Form("fh1BgrTrackPt%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1455 TString strNameBgrZ(Form("fh1BgrZ%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1456 TString strNameBgrXi(Form("fh1BgrXi%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1457
1458 // appendix 'unbinned' to avoid histos with same name after rebinning
1459 TH1F* projPt = (TH1F*) fh2FFTrackPtBgr->ProjectionY(strNameBgrPt+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
1460 TH1F* projZ = (TH1F*) fh2FFZBgr->ProjectionY(strNameBgrZ+"_unBinned",binLo,binUp,"o");
1461 TH1F* projXi = (TH1F*) fh2FFXiBgr->ProjectionY(strNameBgrXi+"_unBinned",binLo,binUp,"o");
1462
1463 if(fNHistoBinsPt[i]) projPt = (TH1F*) projPt->Rebin(fNHistoBinsPt[i],strNameBgrPt,fHistoBinsPt[i]->GetArray());
1464 if(fNHistoBinsZ[i]) projZ = (TH1F*) projZ->Rebin(fNHistoBinsZ[i],strNameBgrZ,fHistoBinsZ[i]->GetArray());
1465 if(fNHistoBinsXi[i]) projXi = (TH1F*) projXi->Rebin(fNHistoBinsXi[i],strNameBgrXi,fHistoBinsXi[i]->GetArray());
1466
1467 projPt->SetNameTitle(strNameBgrPt,"");
1468 projZ->SetNameTitle(strNameBgrZ,"");
1469 projXi->SetNameTitle(strNameBgrXi,"");
1470
1471 // raw bgr = corr level 0
1472 fCorrBgr[0]->AddCorrHistos(i,projPt,projZ,projXi);
ce55b926 1473 }
1474
1475 delete fh1FFJetPtBgr;
1476 delete fh1NJets;
1477 delete fh2FFTrackPtBgr;
1478 delete fh2FFZBgr;
1479 delete fh2FFXiBgr;
39e2b057 1480}
1481
1482
1483//__________________________________________________________________________________________________________
1484void AliFragmentationFunctionCorrections::WriteOutput(TString strfile, TString strdir, Bool_t updateOutfile)
1485{
1486 // write histos to file
1487 // skip histos with 0 entries
1488
1489 TString outfileOption = "RECREATE";
1490 if(updateOutfile) outfileOption = "UPDATE";
1491
1492 TFile f(strfile,outfileOption);
1493
1494 if(!f.IsOpen()){
1495 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1496 return;
1497 }
1498
1499 if(fDebug>0) Printf("%s:%d -- write FF to file %s ",(char*)__FILE__,__LINE__,strfile.Data());
1500
1501 if(strdir && strdir.Length()){
1502 TDirectory* dir = f.mkdir(strdir);
1503 dir->cd();
1504 }
1505
1506 for(Int_t i=0; i<fNJetPtSlices; i++){
1507
1508 for(Int_t c=0; c<fNCorrectionLevels; c++) if(fCorrFF[c]->GetTrackPt(i)->GetEntries()) fCorrFF[c]->GetTrackPt(i)->Write();
1509 for(Int_t c=0; c<fNCorrectionLevels; c++) if(fCorrFF[c]->GetZ(i)->GetEntries()) fCorrFF[c]->GetZ(i)->Write();
1510 for(Int_t c=0; c<fNCorrectionLevels; c++) if(fCorrFF[c]->GetXi(i)->GetEntries()) fCorrFF[c]->GetXi(i)->Write();
1511
39e2b057 1512 for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) if(fCorrBgr[c]->GetTrackPt(i)->GetEntries()) fCorrBgr[c]->GetTrackPt(i)->Write();
1513 for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) if(fCorrBgr[c]->GetZ(i)->GetEntries()) fCorrBgr[c]->GetZ(i)->Write();
1514 for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) if(fCorrBgr[c]->GetXi(i)->GetEntries()) fCorrBgr[c]->GetXi(i)->Write();
1515
1516
1517 if(fh1FFTrackPtBackFolded[i] && fh1FFTrackPtBackFolded[i]->GetEntries()) fh1FFTrackPtBackFolded[i]->Write();
1518 if(fh1FFZBackFolded[i] && fh1FFZBackFolded[i]->GetEntries()) fh1FFZBackFolded[i]->Write();
1519 if(fh1FFXiBackFolded[i] && fh1FFXiBackFolded[i]->GetEntries()) fh1FFXiBackFolded[i]->Write();
1520
1521
1522 if(fh1FFRatioTrackPtFolded[i] && fh1FFRatioTrackPtFolded[i]->GetEntries()) fh1FFRatioTrackPtFolded[i]->Write();
1523 if(fh1FFRatioZFolded[i] && fh1FFRatioZFolded[i]->GetEntries()) fh1FFRatioZFolded[i]->Write();
1524 if(fh1FFRatioXiFolded[i] && fh1FFRatioXiFolded[i]->GetEntries()) fh1FFRatioXiFolded[i]->Write();
1525
1526 if(fh1FFRatioTrackPtBackFolded[i] && fh1FFRatioTrackPtBackFolded[i]->GetEntries()) fh1FFRatioTrackPtBackFolded[i]->Write();
1527 if(fh1FFRatioZBackFolded[i] && fh1FFRatioZBackFolded[i]->GetEntries()) fh1FFRatioZBackFolded[i]->Write();
1528 if(fh1FFRatioXiBackFolded[i] &&fh1FFRatioXiBackFolded[i]->GetEntries()) fh1FFRatioXiBackFolded[i]->Write();
1529
ce55b926 1530 if(fh1FoldingCorrPt[i] && fh1FoldingCorrPt[i]->GetEntries()) fh1FoldingCorrPt[i]->Write(); // TEST!!!
1531 if(fh1FoldingCorrZ[i] && fh1FoldingCorrZ[i]->GetEntries()) fh1FoldingCorrZ[i]->Write(); // TEST!!!
1532 if(fh1FoldingCorrXi[i] && fh1FoldingCorrXi[i]->GetEntries()) fh1FoldingCorrXi[i]->Write(); // TEST!!!
39e2b057 1533 }
1534
1535 // inclusive track pt
1536
1537 for(Int_t c=0; c<fNCorrectionLevelsSinglePt; c++) if(fCorrSinglePt[c]->GetTrackPt(0)->GetEntries()) fCorrSinglePt[c]->GetTrackPt(0)->Write();
1538 if(fh1SingleTrackPtBackFolded) fh1SingleTrackPtBackFolded->Write();
1539 if(fh1RatioSingleTrackPtFolded) fh1RatioSingleTrackPtFolded->Write();
1540 if(fh1RatioSingleTrackPtBackFolded) fh1RatioSingleTrackPtBackFolded->Write();
1541
ce55b926 1542
1543
1544
39e2b057 1545 f.Close();
1546}
1547
1548//____________________________________________________________________________________________________________________________________
1549THnSparse* AliFragmentationFunctionCorrections::TH1toSparse(const TH1F* hist, TString strName, TString strTit, const Bool_t fillConst)
1550{
1551 // copy 1-dimensional histo to THnSparse
1552 // if fillConst TRUE, create THnSparse with same binning as hist but all entries = 1
1553 // histos with variable bin size are supported
1554
1555 // note: function returns pointer to 'new' THnSparse on heap, object needs to be deleted by user
1556
1557 THnSparseF* fhnHist;
1558
1559 Int_t nBins = hist->GetXaxis()->GetNbins();
1560 Int_t nBinsVec[1] = { nBins };
1561
1562 const Double_t* binsVec = hist->GetXaxis()->GetXbins()->GetArray();
1563
1564 if(binsVec){ // binsVec only neq NULL if histo was rebinned before
1565
1566 fhnHist = new THnSparseF(strName,strTit,1,nBinsVec/*,binMinVec,binMaxVec*/);
1567 fhnHist->SetBinEdges(0,binsVec);
1568 }
1569 else{ // uniform bin size
1570
1571 Double_t xMin = hist->GetXaxis()->GetXmin();
1572 Double_t xMax = hist->GetXaxis()->GetXmax();
1573
1574 Double_t binMinVec[1] = { xMin };
1575 Double_t binMaxVec[1] = { xMax };
1576
1577 fhnHist = new THnSparseF(strName,strTit,1,nBinsVec,binMinVec,binMaxVec);
1578 }
1579
1580
1581 for(Int_t bin=1; bin<=nBins; bin++){
1582
1583 Double_t binCenter = fhnHist->GetAxis(0)->GetBinCenter(bin);
1584
1585 Double_t binCoord[] = {binCenter};
1586 fhnHist->Fill(binCoord,1); // initially need to create the bin
1587
1588 Long64_t binIndex = fhnHist->GetBin(binCoord);
1589
1590 Double_t cont = hist->GetBinContent(bin);
1591 Double_t err = hist->GetBinError(bin);
1592
1593 if(fillConst){
1594 cont = 1;
1595 err = 0;
1596 }
1597
1598 fhnHist->SetBinContent(binIndex,cont);
1599 fhnHist->SetBinError(binIndex,err);
1600 }
1601
1602 return fhnHist;
1603}
1604
1605//______________________________________________________________________________________________________________________________________________
ce55b926 1606TH1F* AliFragmentationFunctionCorrections::Unfold(THnSparse* hnHist, const THnSparse* hnResponse, const THnSparse* hnEff, const Int_t nIter,
1607 const Bool_t useCorrelatedErrors, const THnSparse* hnPrior)
39e2b057 1608{
1609 // unfold input histo
1610
1611 AliCFUnfolding unfolding("unfolding","",1,hnResponse,hnEff,hnHist,hnPrior); // arg3: nVar; hnEff required, hnPrior not (defaults to 0x0)
1612 unfolding.SetMaxNumberOfIterations(nIter);
1613 // unfolding.SetMaxChi2PerDOF(1.e-07); // OBSOLETE !!!
1614 // if(useSmoothing) unfolding.UseSmoothing();
1615
1616 if(useCorrelatedErrors) unfolding.SetUseCorrelatedErrors();
1617 unfolding.Unfold();
1618
1619 THnSparse* unfolded = unfolding.GetUnfolded();
1620
1621 TString hnameUnf = hnHist->GetName();
1622 hnameUnf.Append("_unf");
1623 unfolded->SetNameTitle(hnameUnf,hnHist->GetTitle());
ce55b926 1624
1625 TH1F* h1Unfolded = (TH1F*) unfolded->Projection(0);
1626 h1Unfolded->SetNameTitle(hnHist->GetName(),hnHist->GetTitle());
1627
1628 return h1Unfolded;
39e2b057 1629}
1630
1631//___________________________________________________________________________________________________________________________
1632void AliFragmentationFunctionCorrections::UnfoldHistos(const Int_t nIter, const Bool_t useCorrelatedErrors, const Int_t type)
1633{
1634 // loop over jet pt slices and unfold dN/dpt spectra
1635
1636 TString labelCorr = fCorrFF[fNCorrectionLevels-1]->GetLabel();
1637 if(!labelCorr.Contains("unfold")) AddCorrectionLevel("unfold");
1638
1639 for(Int_t i=0; i<fNJetPtSlices; i++){
1640
1641 TH1F* hist = 0;
1aa4f09f 1642 if(type == kFlagPt) hist = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i); // level -2: before unfolding, level -1: unfolded
1643 else if(type == kFlagZ) hist = fCorrFF[fNCorrectionLevels-2]->GetZ(i); // level -2: before unfolding, level -1: unfolded
1644 else if(type == kFlagXi) hist = fCorrFF[fNCorrectionLevels-2]->GetXi(i); // level -2: before unfolding, level -1: unfolded
6daac008 1645 else{
1646 Printf("%s%d unknown type",(char*)__FILE__,__LINE__);
1647 return;
1648 }
b6c0e285 1649
1650 if(!hist){
1651 Printf("%s%d no histo found ",(char*)__FILE__,__LINE__);
1652 return;
1653 }
1654
39e2b057 1655
1656 THnSparse* hnResponse = 0;
6daac008 1657 if(type == kFlagPt) hnResponse = fhnResponsePt[i];
1aa4f09f 1658 else if(type == kFlagZ) hnResponse = fhnResponseZ[i];
1659 else if(type == kFlagXi) hnResponse = fhnResponseXi[i];
39e2b057 1660
1661 TH1F* hPrior = 0;
1662 if(type == kFlagPt && fh1FFTrackPtPrior[i] && ((TString(fh1FFTrackPtPrior[i]->GetName())).Length() > 0) ) hPrior = fh1FFTrackPtPrior[i];
a5592cfa 1663 else if(type == kFlagZ && fh1FFZPrior[i] && ((TString(fh1FFZPrior[i]->GetName())).Length() > 0) ) hPrior = fh1FFZPrior[i];
1664 else if(type == kFlagXi && fh1FFXiPrior[i] && ((TString(fh1FFXiPrior[i]->GetName())).Length() > 0) ) hPrior = fh1FFXiPrior[i];
1aa4f09f 1665
39e2b057 1666
191e5919 1667 TString histNameTHn;
b6c0e285 1668 histNameTHn = hist->GetName();
1669 histNameTHn.ReplaceAll("TH1","THn");
1670
39e2b057 1671
1672 TString priorNameTHn;
1673 if(hPrior){
1674 priorNameTHn = hPrior->GetName();
1675 priorNameTHn.ReplaceAll("TH1","THn");
1676 }
1677
191e5919 1678 TString histNameBackFolded;
b6c0e285 1679 histNameBackFolded = hist->GetName();
39e2b057 1680 histNameBackFolded.Append("_backfold");
b6c0e285 1681
191e5919 1682 TString histNameRatioFolded;
b6c0e285 1683 histNameRatioFolded = hist->GetName();
1684 histNameRatioFolded.ReplaceAll("fh1FF","hRatioFF");
1685
39e2b057 1686 histNameRatioFolded.Append("_unfold");
1687
191e5919 1688 TString histNameRatioBackFolded;
b6c0e285 1689 histNameRatioBackFolded = hist->GetName();
1690 histNameRatioBackFolded.ReplaceAll("fh1FF","hRatioFF");
39e2b057 1691 histNameRatioBackFolded.Append("_backfold");
1692
1693 THnSparse* hnHist = TH1toSparse(hist,histNameTHn,hist->GetTitle());
1694 THnSparse* hnFlatEfficiency = TH1toSparse(hist,"fhnEfficiency","eff",kTRUE); // could optionally also use real eff
1695 THnSparse* hnPrior = 0;
1696 if(hPrior) hnPrior = TH1toSparse(hPrior,priorNameTHn,hPrior->GetTitle());
b6c0e285 1697
ce55b926 1698 //THnSparse* hnUnfolded
1699 // = Unfold(hnHist,hnResponse,hnFlatEfficiency,nIter,useCorrelatedErrors,hnPrior);
1700
1701 //TH1F* hUnfolded = (TH1F*) hnUnfolded->Projection(0);
1702 //hUnfolded->SetNameTitle(hist->GetName(),hist->GetTitle());
1703
1704 TH1F* hUnfolded
39e2b057 1705 = Unfold(hnHist,hnResponse,hnFlatEfficiency,nIter,useCorrelatedErrors,hnPrior);
1706
ce55b926 1707 //TH1F* hUnfolded = (TH1F*) hnUnfolded->Projection(0);
1708 //hUnfolded->SetNameTitle(hist->GetName(),hist->GetTitle());
1709
1710 // errors
1711 for(Int_t bin=1; bin<=hist->GetNbinsX(); bin++){
1712
1713 Double_t contOrg = hist->GetBinContent(bin);
1714 Double_t errOrg = hist->GetBinError(bin);
1715
1716 Double_t contUnf = hUnfolded->GetBinContent(bin);
1717 //Double_t errUnf = hUnfolded->GetBinError(bin);
1718
1719 Double_t relErrOrg = 0;
1720 if(contOrg) relErrOrg = errOrg/contOrg;
1721
1722 // std::cout<<" cont original "<<contOrg<<" err original "<<errOrg<<std::endl;
1723 // std::cout<<" cont unf "<<contUnf<<" errUnf"<<errUnf<<std::endl;
1724 // std::cout<<" err/cont original "<<relErrOrg<<std::endl;
1725 // if(contUnf) std::cout<<" err/conf unf "<<errUnf/contUnf<<std::endl;
1726
1727 hUnfolded->SetBinError(bin,relErrOrg*contUnf);
1728
1729 // if(hUnfolded->GetBinContent(bin)){
1730 // std::cout<<" modified err unfolded "<<hUnfolded->GetBinError(bin)/hUnfolded->GetBinContent(bin)<<std::endl;
1731 // }
1732 }
39e2b057 1733
ce55b926 1734
39e2b057 1735 if(type == kFlagPt) fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hUnfolded,0,0);
1736 if(type == kFlagZ) fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,0,hUnfolded,0);
1737 if(type == kFlagXi) fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,0,0,hUnfolded);
1738
1739 // backfolding: apply response matrix to unfolded spectrum
6daac008 1740 TH1F* hBackFolded = ApplyResponse(hUnfolded,hnResponse);
b6c0e285 1741 hBackFolded->SetNameTitle(histNameBackFolded,hist->GetTitle());
39e2b057 1742
1743 if(type == kFlagPt) fh1FFTrackPtBackFolded[i] = hBackFolded;
1744 if(type == kFlagZ) fh1FFZBackFolded[i] = hBackFolded;
1745 if(type == kFlagXi) fh1FFXiBackFolded[i] = hBackFolded;
1746
1747 // ratio unfolded to original histo
1748 TH1F* hRatioUnfolded = (TH1F*) hUnfolded->Clone(histNameRatioFolded);
1749 hRatioUnfolded->Reset();
b6c0e285 1750 hRatioUnfolded->Divide(hUnfolded,hist,1,1,"B");
39e2b057 1751
1752 if(type == kFlagPt) fh1FFRatioTrackPtFolded[i] = hRatioUnfolded;
1753 if(type == kFlagZ) fh1FFRatioZFolded[i] = hRatioUnfolded;
1754 if(type == kFlagXi) fh1FFRatioXiFolded[i] = hRatioUnfolded;
1755
1756
1757 // ratio backfolded to original histo
1758 TH1F* hRatioBackFolded = (TH1F*) hBackFolded->Clone(histNameRatioBackFolded);
1759 hRatioBackFolded->Reset();
1760 hRatioBackFolded->Divide(hBackFolded,hist,1,1,"B");
1761
1762 if(type == kFlagPt) fh1FFRatioTrackPtBackFolded[i] = hRatioBackFolded;
1763 if(type == kFlagZ) fh1FFRatioZBackFolded[i] = hRatioBackFolded;
1764 if(type == kFlagXi) fh1FFRatioXiBackFolded[i] = hRatioBackFolded;
1765
1766 delete hnHist;
1767 delete hnFlatEfficiency;
39e2b057 1768 }
1769}
1770
1771//_____________________________________________________________________________________________________
1772void AliFragmentationFunctionCorrections::UnfoldPt(const Int_t nIter, const Bool_t useCorrelatedErrors)
1773{
1774
1775 Int_t type = kFlagPt;
1776 UnfoldHistos(nIter, useCorrelatedErrors, type);
1777}
1778
1779//_____________________________________________________________________________________________________
1780void AliFragmentationFunctionCorrections::UnfoldZ(const Int_t nIter, const Bool_t useCorrelatedErrors)
1781{
1782
1783 Int_t type = kFlagZ;
1784 UnfoldHistos(nIter, useCorrelatedErrors, type);
1785}
1786
1787//_____________________________________________________________________________________________________
1788void AliFragmentationFunctionCorrections::UnfoldXi(const Int_t nIter, const Bool_t useCorrelatedErrors)
1789{
1790
1791 Int_t type = kFlagXi;
1792 UnfoldHistos(nIter, useCorrelatedErrors, type);
1793}
1794
ce55b926 1795
39e2b057 1796//______________________________________________________________________________________________
1797TH1F* AliFragmentationFunctionCorrections::ApplyResponse(const TH1F* hist, THnSparse* hnResponse)
1798{
1799 // apply (multiply) response matrix to hist
1800
1801 TH1F* hOut = new TH1F(*hist);
1802 hOut->Reset();
1803
1804 const Int_t axisM = 0;
1805 const Int_t axisT = 1;
1806
1807 Int_t nBinsM = hnResponse->GetAxis(axisM)->GetNbins();
1808 Int_t nBinsT = hnResponse->GetAxis(axisT)->GetNbins();
1809
1810 // response matrix normalization
1811 // do it in this function and not when reading response, since for 'proper' normalization errors are difficult to assign
1812 // so for unfolding proper we leave it to CORRFW ...
1813
ce55b926 1814 Double_t normFacResponse[nBinsT+1];
39e2b057 1815
1816 for(Int_t iT=1; iT<=nBinsT; iT++){
1817
1818 Double_t sumResp = 0;
1819
1820 for(Int_t iM=1; iM<=nBinsM; iM++){
1821
1822 Double_t coordM = hnResponse->GetAxis(axisM)->GetBinCenter(iM);
1823 Double_t coordT = hnResponse->GetAxis(axisT)->GetBinCenter(iT);
1824
1825 Double_t binCoord[] = {coordM,coordT};
1826
1827 Long64_t binIndex = hnResponse->GetBin(binCoord);
1828
1829 Double_t resp = hnResponse->GetBinContent(binIndex);
1830
1831 sumResp += resp;
1832 }
1833
1834 normFacResponse[iT] = 0;
1835 if(sumResp) normFacResponse[iT] = 1/sumResp;
1836 }
1837
1838
1839
1840 for(Int_t iM=1; iM<=nBinsM; iM++){
1841
1842 Double_t contM = 0;
1843
1844 for(Int_t iT=1; iT<=nBinsT; iT++){
1845
1846 Double_t contT = hist->GetBinContent(iT);
1847
1848 Double_t coordM = hnResponse->GetAxis(axisM)->GetBinCenter(iM);
1849 Double_t coordT = hnResponse->GetAxis(axisT)->GetBinCenter(iT);
1850
1851 Double_t binCoord[] = {coordM,coordT};
1852
1853 Long64_t binIndex = hnResponse->GetBin(binCoord);
1854
1855 Double_t resp = hnResponse->GetBinContent(binIndex);
1856
1857 contM += resp*normFacResponse[iT]*contT;
1858 }
1859
1860 hOut->SetBinContent(iM,contM);
1861 }
1862
1863 return hOut;
1864}
1865
1866//_______________________________________________________________________________________________________
1867void AliFragmentationFunctionCorrections::ReadEfficiency(TString strfile, TString strdir, TString strlist)
1868{
1869 // read reconstruction efficiency from file
1870 // argument strlist optional - read from directory strdir if not specified
1871
1872 // temporary histos to hold histos from file
1873 TH1F* hEffPt[fNJetPtSlices];
1874 TH1F* hEffZ[fNJetPtSlices];
1875 TH1F* hEffXi[fNJetPtSlices];
1876
1877 for(Int_t i=0; i<fNJetPtSlices; i++) hEffPt[i] = 0;
1878 for(Int_t i=0; i<fNJetPtSlices; i++) hEffZ[i] = 0;
1879 for(Int_t i=0; i<fNJetPtSlices; i++) hEffXi[i] = 0;
1880
1881 TFile f(strfile,"READ");
1882
1883 if(!f.IsOpen()){
1884 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1885 return;
1886 }
1887
1888 if(fDebug>0) Printf("%s:%d -- read efficiencies from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
1889
1890 if(strdir && strdir.Length()) gDirectory->cd(strdir);
1891
1892 TList* list = 0;
1893
1894 if(strlist && strlist.Length()){
1895
1896 if(!(list = (TList*) gDirectory->Get(strlist))){
1897 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1898 return;
1899 }
1900 }
1901
1902 for(Int_t i=0; i<fNJetPtSlices; i++){
1903
1904 Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
1905 Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
1906
1907 TString strNameEffPt(Form("hEffPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
1908 TString strNameEffZ(Form("hEffZ_%02d_%02d",jetPtLoLim,jetPtUpLim));
1909 TString strNameEffXi(Form("hEffXi_%02d_%02d",jetPtLoLim,jetPtUpLim));
1910
1911
1912 if(list){
1913 hEffPt[i] = (TH1F*) list->FindObject(strNameEffPt);
1914 hEffZ[i] = (TH1F*) list->FindObject(strNameEffZ);
1915 hEffXi[i] = (TH1F*) list->FindObject(strNameEffXi);
1916 }
1917 else{
1918 hEffPt[i] = (TH1F*) gDirectory->Get(strNameEffPt);
1919 hEffZ[i] = (TH1F*) gDirectory->Get(strNameEffZ);
1920 hEffXi[i] = (TH1F*) gDirectory->Get(strNameEffXi);
1921 }
1922
1923 if(!hEffPt[i]){
1924 Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffPt.Data());
1925 }
1926
1927 if(!hEffZ[i]){
1928 Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffZ.Data());
1929 }
1930
1931 if(!hEffXi[i]){
1932 Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffXi.Data());
1933 }
1934
1935
1936 if(fNHistoBinsPt[i]) hEffPt[i] = (TH1F*) hEffPt[i]->Rebin(fNHistoBinsPt[i],strNameEffPt+"_rebin",fHistoBinsPt[i]->GetArray());
1937 if(fNHistoBinsZ[i]) hEffZ[i] = (TH1F*) hEffZ[i]->Rebin(fNHistoBinsZ[i],strNameEffZ+"_rebin",fHistoBinsZ[i]->GetArray());
1938 if(fNHistoBinsXi[i]) hEffXi[i] = (TH1F*) hEffXi[i]->Rebin(fNHistoBinsXi[i],strNameEffXi+"_rebin",fHistoBinsXi[i]->GetArray());
1939
1940 if(hEffPt[i]) hEffPt[i]->SetDirectory(0);
1941 if(hEffZ[i]) hEffZ[i]->SetDirectory(0);
1942 if(hEffXi[i]) hEffXi[i]->SetDirectory(0);
1943
1944 } // jet slices loop
1945
1946 f.Close();
1947
1948 for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
1949 if(hEffPt[i]) new(fh1EffPt[i]) TH1F(*hEffPt[i]);
1950 if(hEffZ[i]) new(fh1EffZ[i]) TH1F(*hEffZ[i]);
1951 if(hEffXi[i]) new(fh1EffXi[i]) TH1F(*hEffXi[i]);
1952 }
1953}
1954
1955//___________________________________________________________________________________________________________
1956void AliFragmentationFunctionCorrections::ReadBgrEfficiency(TString strfile, TString strdir, TString strlist)
1957{
1958 // read bgr eff from file
1959 // argument strlist optional - read from directory strdir if not specified
1960
1961 TH1F* hEffPtBgr[fNJetPtSlices];
1962 TH1F* hEffZBgr [fNJetPtSlices];
1963 TH1F* hEffXiBgr[fNJetPtSlices];
1964
1965 for(Int_t i=0; i<fNJetPtSlices; i++) hEffPtBgr[i] = 0;
1966 for(Int_t i=0; i<fNJetPtSlices; i++) hEffZBgr[i] = 0;
1967 for(Int_t i=0; i<fNJetPtSlices; i++) hEffXiBgr[i] = 0;
1968
1969
1970 TFile f(strfile,"READ");
1971
1972 if(!f.IsOpen()){
1973 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1974 return;
1975 }
1976
1977 if(fDebug>0) Printf("%s:%d -- read bgr efficiencies from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
1978
1979 if(strdir && strdir.Length()) gDirectory->cd(strdir);
1980
1981 TList* list = 0;
1982
1983 if(strlist && strlist.Length()){
1984
1985 if(!(list = (TList*) gDirectory->Get(strlist))){
1986 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1987 return;
1988 }
1989 }
1990
1991 for(Int_t i=0; i<fNJetPtSlices; i++){
1992
1993 Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
1994 Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
1995
1996 TString strNameEffPtBgr(Form("hEffPtBgr%02dto%02d",jetPtLoLim,jetPtUpLim));
1997 TString strNameEffZBgr(Form("hEffZBgr%02dto%02d",jetPtLoLim,jetPtUpLim));
1998 TString strNameEffXiBgr(Form("hEffXiBgr%02dto%02d",jetPtLoLim,jetPtUpLim));
1999
2000
2001 if(list){
2002 hEffPtBgr[i] = (TH1F*) list->FindObject(strNameEffPtBgr);
2003 hEffZBgr[i] = (TH1F*) list->FindObject(strNameEffZBgr);
2004 hEffXiBgr[i] = (TH1F*) list->FindObject(strNameEffXiBgr);
2005 }
2006 else{
2007 hEffPtBgr[i] = (TH1F*) gDirectory->Get(strNameEffPtBgr);
2008 hEffZBgr[i] = (TH1F*) gDirectory->Get(strNameEffZBgr);
2009 hEffXiBgr[i] = (TH1F*) gDirectory->Get(strNameEffXiBgr);
2010 }
2011
2012 if(!hEffPtBgr[i]){
2013 Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffPtBgr.Data());
2014 }
2015
2016 if(!hEffZBgr[i]){
2017 Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffZBgr.Data());
2018 }
2019
2020 if(!hEffXiBgr[i]){
2021 Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffXiBgr.Data());
2022 }
2023
2024
2025 if(fNHistoBinsPt[i]) hEffPtBgr[i] = (TH1F*) hEffPtBgr[i]->Rebin(fNHistoBinsPt[i],strNameEffPtBgr+"_rebin",fHistoBinsPt[i]->GetArray());
2026 if(fNHistoBinsZ[i]) hEffZBgr[i] = (TH1F*) hEffZBgr[i]->Rebin(fNHistoBinsZ[i],strNameEffZBgr+"_rebin",fHistoBinsZ[i]->GetArray());
2027 if(fNHistoBinsXi[i]) hEffXiBgr[i] = (TH1F*) hEffXiBgr[i]->Rebin(fNHistoBinsXi[i],strNameEffXiBgr+"_rebin",fHistoBinsXi[i]->GetArray());
2028
2029 if(hEffPtBgr[i]) hEffPtBgr[i]->SetDirectory(0);
2030 if(hEffZBgr[i]) hEffZBgr[i]->SetDirectory(0);
2031 if(hEffXiBgr[i]) hEffXiBgr[i]->SetDirectory(0);
2032
2033 } // jet slices loop
2034
2035 f.Close();
2036
2037 for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
2038 if(hEffPtBgr[i]) new(fh1EffBgrPt[i]) TH1F(*hEffPtBgr[i]);
2039 if(hEffZBgr[i]) new(fh1EffBgrZ[i]) TH1F(*hEffZBgr[i]);
2040 if(hEffXiBgr[i]) new(fh1EffBgrXi[i]) TH1F(*hEffXiBgr[i]);
2041 }
2042}
2043
2044// ________________________________________________
2045void AliFragmentationFunctionCorrections::EffCorr()
2046{
2047 // apply efficiency correction
2048
2049 AddCorrectionLevel("eff");
2050
2051 for(Int_t i=0; i<fNJetPtSlices; i++){
2052
2053 TH1F* histPt = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i);
2054 TH1F* histZ = fCorrFF[fNCorrectionLevels-2]->GetZ(i);
2055 TH1F* histXi = fCorrFF[fNCorrectionLevels-2]->GetXi(i);
2056
2057 TString histNamePt = histPt->GetName();
2058 TString histNameZ = histZ->GetName();
2059 TString histNameXi = histXi->GetName();
2060
2061
2062 TH1F* hFFTrackPtEffCorr = (TH1F*) histPt->Clone(histNamePt);
2063 hFFTrackPtEffCorr->Divide(histPt,fh1EffPt[i],1,1,"");
2064
2065 TH1F* hFFZEffCorr = (TH1F*) histZ->Clone(histNameZ);
2066 hFFZEffCorr->Divide(histZ,fh1EffZ[i],1,1,"");
2067
2068 TH1F* hFFXiEffCorr = (TH1F*) histXi->Clone(histNameXi);
2069 hFFXiEffCorr->Divide(histXi,fh1EffXi[i],1,1,"");
2070
2071 fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hFFTrackPtEffCorr,hFFZEffCorr,hFFXiEffCorr);
2072 }
2073}
2074
2075//___________________________________________________
2076void AliFragmentationFunctionCorrections::EffCorrBgr()
2077{
2078 // apply efficiency correction to bgr distributions
2079
2080 AddCorrectionLevelBgr("eff");
2081
2082 Printf("%s:%d -- apply efficiency correction, corrLevel %d",(char*)__FILE__,__LINE__,fNCorrectionLevels-1);
2083
2084 for(Int_t i=0; i<fNJetPtSlices; i++){
2085
2086 TH1F* histPt = fCorrBgr[fNCorrectionLevelsBgr-2]->GetTrackPt(i);
2087 TH1F* histZ = fCorrBgr[fNCorrectionLevelsBgr-2]->GetZ(i);
2088 TH1F* histXi = fCorrBgr[fNCorrectionLevelsBgr-2]->GetXi(i);
2089
2090 TString histNamePt = histPt->GetName();
2091 TString histNameZ = histZ->GetName();
2092 TString histNameXi = histXi->GetName();
2093
2094
2095 TH1F* hFFTrackPtEffCorr = (TH1F*) histPt->Clone(histNamePt);
2096 hFFTrackPtEffCorr->Divide(histPt,fh1EffPt[i],1,1,"");
2097
2098 TH1F* hFFZEffCorr = (TH1F*) histZ->Clone(histNameZ);
2099 hFFZEffCorr->Divide(histZ,fh1EffZ[i],1,1,"");
2100
2101 TH1F* hFFXiEffCorr = (TH1F*) histXi->Clone(histNameXi);
2102 hFFXiEffCorr->Divide(histXi,fh1EffXi[i],1,1,"");
2103
2104 fCorrBgr[fNCorrectionLevelsBgr-1]->AddCorrHistos(i,hFFTrackPtEffCorr,hFFZEffCorr,hFFXiEffCorr);
2105 }
2106}
2107
39e2b057 2108//_____________________________________________________
ce55b926 2109void AliFragmentationFunctionCorrections::SubtractBgr(Double_t sysErr)
39e2b057 2110{
2111 // subtract bgr distribution from FF
2112 // the current corr level is used for both
2113
2114 AddCorrectionLevel("bgrSub");
2115
ce55b926 2116 fstream ascii_out_pt;
2117 fstream ascii_out_z;
2118 fstream ascii_out_xi;
2119
2120 if(sysErr) ascii_out_pt.open("sysErrUE_pt.txt",std::ios::out);
2121 if(sysErr) ascii_out_z.open("sysErrUE_z.txt",std::ios::out);
2122 if(sysErr) ascii_out_xi.open("sysErrUE_xi.txt",std::ios::out);
2123
2124 for(Int_t i=0; i<fNJetPtSlices; i++){ // jet slices
39e2b057 2125
2126 TH1F* histPt = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i);
2127 TH1F* histZ = fCorrFF[fNCorrectionLevels-2]->GetZ(i);
2128 TH1F* histXi = fCorrFF[fNCorrectionLevels-2]->GetXi(i);
2129
2130 TH1F* histPtBgr = fCorrBgr[fNCorrectionLevelsBgr-1]->GetTrackPt(i);
2131 TH1F* histZBgr = fCorrBgr[fNCorrectionLevelsBgr-1]->GetZ(i);
2132 TH1F* histXiBgr = fCorrBgr[fNCorrectionLevelsBgr-1]->GetXi(i);
2133
2134 TString histNamePt = histPt->GetName();
2135 TString histNameZ = histZ->GetName();
2136 TString histNameXi = histXi->GetName();
2137
2138
2139 TH1F* hFFTrackPtBgrSub = (TH1F*) histPt->Clone(histNamePt);
2140 hFFTrackPtBgrSub->Add(histPtBgr,-1);
2141
2142 TH1F* hFFZBgrSub = (TH1F*) histZ->Clone(histNameZ);
2143 hFFZBgrSub->Add(histZBgr,-1);
2144
2145 TH1F* hFFXiBgrSub = (TH1F*) histXi->Clone(histNameXi);
2146 hFFXiBgrSub->Add(histXiBgr,-1);
2147
2148 fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hFFTrackPtBgrSub,hFFZBgrSub,hFFXiBgrSub);
ce55b926 2149
2150 if(sysErr){
2151
2152 for(Int_t bin=1; bin<=histPt->GetNbinsX(); bin++){
2153
2154 //Double_t sigPlBgr = histPt->GetBinContent(bin);
2155 Double_t sig = hFFTrackPtBgrSub->GetBinContent(bin);
2156 Double_t bgr = histPtBgr->GetBinContent(bin);
2157
2158 Double_t relErr = 0;
2159 if(sig>0) relErr = sysErr*bgr/sig;
2160
2161// std::cout<<" pt bin "<<bin<<" mean "<<histPt->GetBinCenter(bin)
2162// <<" sigPlBgr "<<sigPlBgr<<" sig "<<sig<<" bgr "<<bgr<<" relErr "<<relErr<<std::endl;
2163
2164 ascii_out_pt<<i<<" "<<bin<<" "<<relErr<<std::endl;
2165 }
2166
2167
2168 for(Int_t bin=1; bin<=histZ->GetNbinsX(); bin++){
2169
2170 //Double_t sigPlBgr = histZ->GetBinContent(bin);
2171 Double_t sig = hFFZBgrSub->GetBinContent(bin);
2172 Double_t bgr = histZBgr->GetBinContent(bin);
2173
2174 Double_t relErr = 0;
2175 if(sig>0) relErr = sysErr*bgr/sig;
2176
2177 std::cout<<" z bin "<<bin<<" mean "<<histZ->GetBinCenter(bin)
2178 <<" sigPlBgr "<<histZ->GetBinContent(bin)<<" sig "<<sig<<" bgr "<<bgr<<" relErr "<<relErr<<std::endl;
2179
2180 ascii_out_z<<i<<" "<<bin<<" "<<relErr<<std::endl;
2181 }
2182
2183
2184 for(Int_t bin=1; bin<=histXi->GetNbinsX(); bin++){
2185
2186 //Double_t sigPlBgr = histXi->GetBinContent(bin);
2187 Double_t sig = hFFXiBgrSub->GetBinContent(bin);
2188 Double_t bgr = histXiBgr->GetBinContent(bin);
2189
2190 Double_t relErr = 0;
2191 if(sig>0) relErr = sysErr*bgr/sig;
2192
2193 std::cout<<" xi bin "<<bin<<" mean "<<histXi->GetBinCenter(bin)
2194 <<" sigPlBgr "<<histXi->GetBinContent(bin)<<" sig "<<sig<<" bgr "<<bgr<<" relErr "<<relErr<<std::endl;
2195
2196 ascii_out_xi<<i<<" "<<bin<<" "<<relErr<<std::endl;
2197 }
2198 }
39e2b057 2199 }
ce55b926 2200
2201 if(sysErr) ascii_out_pt.close();
2202 if(sysErr) ascii_out_xi.close();
2203
2204
39e2b057 2205}
2206
2207//________________________________________________________________________________________________________________
2208void AliFragmentationFunctionCorrections::WriteSingleTrackEff(TString strInfile, TString strID, TString strOutfile,
2209 Bool_t updateOutfile, TString strOutDir,TString strPostfix)
2210{
2211 // read task ouput from MC and write single track eff - standard dir/list
2212
b541fbca 2213 TString strdir = "PWGJE_FragmentationFunction_" + strID;
39e2b057 2214 TString strlist = "fracfunc_" + strID;
2215
2216 WriteSingleTrackEff(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir,strPostfix);
2217}
2218
2219//___________________________________________________________________________________________________________________________________
2220void AliFragmentationFunctionCorrections::WriteSingleTrackEff(TString strInfile, TString strdir, TString strlist,
2221 TString strOutfile, Bool_t updateOutfile, TString strOutDir,TString strPostfix)
2222{
2223 // read task output from MC and write single track eff as function of pt, eta, phi
2224 // argument strlist optional - read from directory strdir if not specified
2225 // write eff to file stroutfile - by default only 'update' file (don't overwrite)
2226
2227
2228 TH1D* hdNdptTracksMCPrimGen;
2229 TH2D* hdNdetadphiTracksMCPrimGen;
2230
2231 TH1D* hdNdptTracksMCPrimRec;
2232 TH2D* hdNdetadphiTracksMCPrimRec;
2233
2234
2235 TFile f(strInfile,"READ");
2236
2237 if(!f.IsOpen()){
2238 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2239 return;
2240 }
2241
2242 if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: open task ouput file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
2243
2244 if(strdir && strdir.Length()){
2245 if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: change dir to %s",(char*)__FILE__,__LINE__,strdir.Data());
2246 gDirectory->cd(strdir);
2247 }
2248
2249 TList* list = 0;
2250
2251 if(strlist && strlist.Length()){
2252
2253 if(!(list = (TList*) gDirectory->Get(strlist))){
2254 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2255 return;
2256 }
2257 }
2258
2259
2260 TString hnamePtRecEffGen = "fh1TrackQAPtRecEffGen";
2261 if(strPostfix.Length()) hnamePtRecEffGen.Form("fh1TrackQAPtRecEffGen_%s",strPostfix.Data());
2262
2263 TString hnamePtRecEffRec = "fh1TrackQAPtRecEffRec";
2264 if(strPostfix.Length()) hnamePtRecEffRec.Form("fh1TrackQAPtRecEffRec_%s",strPostfix.Data());
2265
2266 TString hnameEtaPhiRecEffGen = "fh2TrackQAEtaPhiRecEffGen";
2267 if(strPostfix.Length()) hnameEtaPhiRecEffGen.Form("fh2TrackQAEtaPhiRecEffGen_%s",strPostfix.Data());
2268
2269 TString hnameEtaPhiRecEffRec = "fh2TrackQAEtaPhiRecEffRec";
2270 if(strPostfix.Length()) hnameEtaPhiRecEffRec.Form("fh2TrackQAEtaPhiRecEffRec_%s",strPostfix.Data());
2271
2272
2273 if(list){
2274 hdNdptTracksMCPrimGen = (TH1D*) list->FindObject(hnamePtRecEffGen);
2275 hdNdetadphiTracksMCPrimGen = (TH2D*) list->FindObject(hnameEtaPhiRecEffGen);
2276
2277 hdNdptTracksMCPrimRec = (TH1D*) list->FindObject(hnamePtRecEffRec);
2278 hdNdetadphiTracksMCPrimRec = (TH2D*) list->FindObject(hnameEtaPhiRecEffRec);
2279 }
2280 else{
2281 hdNdptTracksMCPrimGen = (TH1D*) gDirectory->Get(hnamePtRecEffGen);
2282 hdNdetadphiTracksMCPrimGen = (TH2D*) gDirectory->Get(hnameEtaPhiRecEffGen);
2283
2284 hdNdptTracksMCPrimRec = (TH1D*) gDirectory->Get(hnamePtRecEffRec);
2285 hdNdetadphiTracksMCPrimRec = (TH2D*) gDirectory->Get(hnameEtaPhiRecEffRec);
2286 }
2287
2288 hdNdptTracksMCPrimGen->SetDirectory(0);
2289 hdNdetadphiTracksMCPrimGen->SetDirectory(0);
2290 hdNdptTracksMCPrimRec->SetDirectory(0);
2291 hdNdetadphiTracksMCPrimRec->SetDirectory(0);
2292
2293 f.Close();
2294
2295 // projections: dN/deta, dN/dphi
2296
2297 TH1D* hdNdetaTracksMCPrimGen = (TH1D*) hdNdetadphiTracksMCPrimGen->ProjectionX("hdNdetaTracksMcPrimGen");
2298 TH1D* hdNdphiTracksMCPrimGen = (TH1D*) hdNdetadphiTracksMCPrimGen->ProjectionY("hdNdphiTracksMcPrimGen");
2299
2300 TH1D* hdNdetaTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionX("hdNdetaTracksMcPrimRec");
2301 TH1D* hdNdphiTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionY("hdNdphiTracksMcPrimRec");
2302
2303 // rebin
2304
2305 TString strNamePtGen = "hTrackPtGenPrim";
2306 TString strNamePtRec = "hTrackPtRecPrim";
2307
2308 if(fNHistoBinsSinglePt) hdNdptTracksMCPrimGen = (TH1D*) hdNdptTracksMCPrimGen->Rebin(fNHistoBinsSinglePt,strNamePtGen,fHistoBinsSinglePt->GetArray());
2309 if(fNHistoBinsSinglePt) hdNdptTracksMCPrimRec = (TH1D*) hdNdptTracksMCPrimRec->Rebin(fNHistoBinsSinglePt,strNamePtRec,fHistoBinsSinglePt->GetArray());
2310
2311 // efficiency: divide
2312
2313 TString hNameTrackEffPt = "hSingleTrackEffPt";
2314 if(strPostfix.Length()) hNameTrackEffPt.Form("hSingleTrackEffPt_%s",strPostfix.Data());
2315
2316 TString hNameTrackEffEta = "hSingleTrackEffEta";
2317 if(strPostfix.Length()) hNameTrackEffEta.Form("hSingleTrackEffEta_%s",strPostfix.Data());
2318
2319 TString hNameTrackEffPhi = "hSingleTrackEffPhi";
2320 if(strPostfix.Length()) hNameTrackEffPhi.Form("hSingleTrackEffPhi_%s",strPostfix.Data());
2321
2322
2323 TH1F* hSingleTrackEffPt = (TH1F*) hdNdptTracksMCPrimRec->Clone(hNameTrackEffPt);
2324 hSingleTrackEffPt->Divide(hdNdptTracksMCPrimRec,hdNdptTracksMCPrimGen,1,1,"B"); // binominal errors
2325
2326 TH1F* hSingleTrackEffEta = (TH1F*) hdNdetaTracksMCPrimRec->Clone(hNameTrackEffEta);
2327 hSingleTrackEffEta->Divide(hdNdetaTracksMCPrimRec,hdNdetaTracksMCPrimGen,1,1,"B"); // binominal errors
2328
2329 TH1F* hSingleTrackEffPhi = (TH1F*) hdNdphiTracksMCPrimRec->Clone(hNameTrackEffPhi);
2330 hSingleTrackEffPhi->Divide(hdNdphiTracksMCPrimRec,hdNdphiTracksMCPrimGen,1,1,"B"); // binominal errors
2331
2332
2333 TString outfileOption = "RECREATE";
2334 if(updateOutfile) outfileOption = "UPDATE";
2335
2336 TFile out(strOutfile,outfileOption);
2337
2338 if(!out.IsOpen()){
2339 Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2340 return;
2341 }
2342
2343 if(fDebug>0) Printf("%s:%d -- write efficiency to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2344
2345 if(strOutDir && strOutDir.Length()){
2346
2347 TDirectory* dir;
2348 if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd();
2349 else{
2350 dir = out.mkdir(strOutDir);
2351 dir->cd();
2352 }
2353 }
2354
2355 hSingleTrackEffPt->Write();
2356 hSingleTrackEffEta->Write();
2357 hSingleTrackEffPhi->Write();
2358
2359 out.Close();
ce55b926 2360
2361 delete hdNdptTracksMCPrimGen;
2362 delete hdNdetadphiTracksMCPrimGen;
2363 delete hdNdptTracksMCPrimRec;
2364 delete hdNdetadphiTracksMCPrimRec;
39e2b057 2365}
2366
2367//________________________________________________________________________________________________________________
2368void AliFragmentationFunctionCorrections::WriteSingleTrackSecCorr(TString strInfile, TString strID, TString strOutfile,
2369 Bool_t updateOutfile, TString strOutDir)
2370{
2371 // read task ouput from MC and write single track eff - standard dir/list
2372
b541fbca 2373 TString strdir = "PWGJE_FragmentationFunction_" + strID;
39e2b057 2374 TString strlist = "fracfunc_" + strID;
2375
2376 WriteSingleTrackSecCorr(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir);
2377}
2378
2379//___________________________________________________________________________________________________________________________________
2380void AliFragmentationFunctionCorrections::WriteSingleTrackSecCorr(TString strInfile, TString strdir, TString strlist,
2381 TString strOutfile, Bool_t updateOutfile, TString strOutDir)
2382{
2383 // read task output from MC and write single track secondaries contamination correction as function of pt, eta, phi
2384 // argument strlist optional - read from directory strdir if not specified
2385 // write corr factor to file stroutfile - by default only 'update' file (don't overwrite)
2386
2387 TH1D* hdNdptTracksMCPrimRec;
2388 TH2D* hdNdetadphiTracksMCPrimRec;
2389
ce55b926 2390 TH1D* hdNdptTracksMCSecRecNS;
2391 TH2D* hdNdetadphiTracksMCSecRecNS;
2392
2393 TH1D* hdNdptTracksMCSecRecSsc;
2394 TH2D* hdNdetadphiTracksMCSecRecSsc;
39e2b057 2395
2396 TFile f(strInfile,"READ");
2397
2398 if(!f.IsOpen()){
2399 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2400 return;
2401 }
2402
2403 if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
2404
2405 if(strdir && strdir.Length()) gDirectory->cd(strdir);
2406
2407 TList* list = 0;
2408
2409 if(strlist && strlist.Length()){
2410
2411 if(!(list = (TList*) gDirectory->Get(strlist))){
2412 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2413 return;
2414 }
2415 }
2416
2417
2418 if(list){
ce55b926 2419 hdNdptTracksMCPrimRec = (TH1D*) list->FindObject("fh1TrackQAPtRecEffGen");
2420 hdNdetadphiTracksMCPrimRec = (TH2D*) list->FindObject("fh2TrackQAEtaPhiRecEffGen");
39e2b057 2421
ce55b926 2422 hdNdptTracksMCSecRecNS = (TH1D*) list->FindObject("fh1TrackQAPtSecRecNS");
2423 hdNdetadphiTracksMCSecRecNS = (TH2D*) list->FindObject("fh2TrackQAEtaPhiSecRecNS");
2424
2425 hdNdptTracksMCSecRecSsc = (TH1D*) list->FindObject("fh1TrackQAPtSecRecSsc");
2426 hdNdetadphiTracksMCSecRecSsc = (TH2D*) list->FindObject("fh2TrackQAEtaPhiSecRecSsc");
2427
39e2b057 2428 }
2429 else{
ce55b926 2430 hdNdptTracksMCPrimRec = (TH1D*) gDirectory->Get("fh1TrackQAPtRecEffGen");
2431 hdNdetadphiTracksMCPrimRec = (TH2D*) gDirectory->Get("fh2TrackQAEtaPhiRecEffGen");
2432
2433 hdNdptTracksMCSecRecNS = (TH1D*) gDirectory->Get("fh1TrackQAPtSecRecNS");
2434 hdNdetadphiTracksMCSecRecNS = (TH2D*) gDirectory->Get("fh2TrackQAEtaPhiSecRecNS");
39e2b057 2435
ce55b926 2436 hdNdptTracksMCSecRecSsc = (TH1D*) gDirectory->Get("fh1TrackQAPtSecRecSsc");
2437 hdNdetadphiTracksMCSecRecSsc = (TH2D*) gDirectory->Get("fh2TrackQAEtaPhiSecRecSsc");
39e2b057 2438 }
2439
2440 hdNdptTracksMCPrimRec->SetDirectory(0);
2441 hdNdetadphiTracksMCPrimRec->SetDirectory(0);
ce55b926 2442
2443 hdNdptTracksMCSecRecNS->SetDirectory(0);
2444 hdNdetadphiTracksMCSecRecNS->SetDirectory(0);
2445
2446 hdNdptTracksMCSecRecSsc->SetDirectory(0);
2447 hdNdetadphiTracksMCSecRecSsc->SetDirectory(0);
39e2b057 2448
2449 f.Close();
2450
2451 // projections: dN/deta, dN/dphi
2452
ce55b926 2453 TH1D* hdNdetaTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionX("hdNdetaTracksMcPrimRec");
2454 TH1D* hdNdphiTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionY("hdNdphiTracksMcPrimRec");
39e2b057 2455
ce55b926 2456 TH1D* hdNdetaTracksMCSecRecNS = (TH1D*) hdNdetadphiTracksMCSecRecNS->ProjectionX("hdNdetaTracksMcSecRecNS");
2457 TH1D* hdNdphiTracksMCSecRecNS = (TH1D*) hdNdetadphiTracksMCSecRecNS->ProjectionY("hdNdphiTracksMcSecRecNS");
39e2b057 2458
ce55b926 2459 TH1D* hdNdetaTracksMCSecRecSsc = (TH1D*) hdNdetadphiTracksMCSecRecSsc->ProjectionX("hdNdetaTracksMcSecRecSsc");
2460 TH1D* hdNdphiTracksMCSecRecSsc = (TH1D*) hdNdetadphiTracksMCSecRecSsc->ProjectionY("hdNdphiTracksMcSecRecSsc");
39e2b057 2461
2462 // rebin
2463
2464 TString strNamePtPrim = "hTrackPtPrim";
2465 TString strNamePtSec = "hTrackPtSec";
2466
ce55b926 2467 if(fNHistoBinsSinglePt) hdNdptTracksMCPrimRec = (TH1D*) hdNdptTracksMCPrimRec->Rebin(fNHistoBinsSinglePt,strNamePtPrim,fHistoBinsSinglePt->GetArray());
2468 if(fNHistoBinsSinglePt) hdNdptTracksMCSecRecNS = (TH1D*) hdNdptTracksMCSecRecNS->Rebin(fNHistoBinsSinglePt,strNamePtSec,fHistoBinsSinglePt->GetArray());
2469 if(fNHistoBinsSinglePt) hdNdptTracksMCSecRecSsc = (TH1D*) hdNdptTracksMCSecRecSsc->Rebin(fNHistoBinsSinglePt,strNamePtSec,fHistoBinsSinglePt->GetArray());
39e2b057 2470
2471
2472 // secondary correction factor: divide prim/(prim+sec)
2473
ce55b926 2474 TH1F* hSingleTrackSecCorrPt = (TH1F*) hdNdptTracksMCPrimRec->Clone("hSingleTrackSecCorrPt");
2475 TH1F* hSumPrimSecPt = (TH1F*) hdNdptTracksMCPrimRec->Clone("hSumPrimSecPt");
2476 hSumPrimSecPt->Add(hdNdptTracksMCSecRecNS);
2477 hSumPrimSecPt->Add(hdNdptTracksMCSecRecSsc);
39e2b057 2478 hSingleTrackSecCorrPt->Divide(hdNdptTracksMCPrimRec,hSumPrimSecPt,1,1,"B"); // binominal errors
2479
ce55b926 2480
2481 TH1F* hSingleTrackSecCorrEta = (TH1F*) hdNdetaTracksMCPrimRec->Clone("hSingleTrackSecCorrEta");
2482 TH1F* hSumPrimSecEta = (TH1F*) hdNdetaTracksMCPrimRec->Clone("hSumPrimSecEta");
2483 hSumPrimSecEta->Add(hdNdetaTracksMCSecRecNS);
2484 hSumPrimSecEta->Add(hdNdetaTracksMCSecRecSsc);
39e2b057 2485 hSingleTrackSecCorrEta->Divide(hdNdetaTracksMCPrimRec,hSumPrimSecEta,1,1,"B"); // binominal errors
2486
ce55b926 2487 TH1F* hSingleTrackSecCorrPhi = (TH1F*) hdNdphiTracksMCPrimRec->Clone("hSingleTrackSecCorrPhi");
2488 TH1F* hSumPrimSecPhi = (TH1F*) hdNdphiTracksMCPrimRec->Clone("hSumPrimSecPhi");
2489 hSumPrimSecPhi->Add(hdNdphiTracksMCSecRecNS);
2490 hSumPrimSecPhi->Add(hdNdphiTracksMCSecRecSsc);
39e2b057 2491 hSingleTrackSecCorrPhi->Divide(hdNdphiTracksMCPrimRec,hSumPrimSecPhi,1,1,"B"); // binominal errors
ce55b926 2492
2493 //
39e2b057 2494
2495 TString outfileOption = "RECREATE";
2496 if(updateOutfile) outfileOption = "UPDATE";
2497
2498 TFile out(strOutfile,outfileOption);
2499
2500 if(!out.IsOpen()){
2501 Printf("%s:%d -- error opening secCorr output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2502 return;
2503 }
2504
2505 if(fDebug>0) Printf("%s:%d -- write secCorr to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2506
2507 if(strOutDir && strOutDir.Length()){
2508
2509 TDirectory* dir;
2510 if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd();
2511 else{
2512 dir = out.mkdir(strOutDir);
2513 dir->cd();
2514 }
2515 }
2516
ce55b926 2517 hdNdptTracksMCSecRecNS->Write();
2518 hdNdetaTracksMCSecRecNS->Write();
2519 hdNdphiTracksMCSecRecNS->Write();
2520
2521 hdNdptTracksMCSecRecSsc->Write();
2522 hdNdetaTracksMCSecRecSsc->Write();
2523 hdNdphiTracksMCSecRecSsc->Write();
39e2b057 2524
2525 hSingleTrackSecCorrPt->Write();
2526 hSingleTrackSecCorrEta->Write();
2527 hSingleTrackSecCorrPhi->Write();
2528
2529 out.Close();
ce55b926 2530
2531 delete hdNdptTracksMCPrimRec;
2532 delete hdNdetadphiTracksMCPrimRec;
2533 delete hdNdptTracksMCSecRecNS;
2534 delete hdNdetadphiTracksMCSecRecNS;
2535 delete hdNdptTracksMCSecRecSsc;
2536 delete hdNdetadphiTracksMCSecRecSsc;
39e2b057 2537}
2538
2539//________________________________________________________________________________________________________________
2540void AliFragmentationFunctionCorrections::WriteSingleResponse(TString strInfile, TString strID, TString strOutfile,
2541 Bool_t updateOutfile, TString strOutDir)
2542{
2543 // read task ouput from MC and write single track eff - standard dir/list
2544
b541fbca 2545 TString strdir = "PWGJE_FragmentationFunction_" + strID;
39e2b057 2546 TString strlist = "fracfunc_" + strID;
2547
2548 WriteSingleResponse(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir);
2549}
2550
2551
2552//_____________________________________________________________________________________________________________________________________
2553void AliFragmentationFunctionCorrections::WriteSingleResponse(TString strInfile, TString strdir, TString strlist,
2554 TString strOutfile, Bool_t updateOutfile, TString strOutDir)
2555{
2556 // read 2d THnSparse response matrices in pt from file
2557 // project TH2
2558 // write to strOutfile
2559
2560 THnSparse* hnResponseSinglePt;
2561 TH2F* h2ResponseSinglePt;
2562
2563 TFile f(strInfile,"READ");
2564
2565 if(!f.IsOpen()){
2566 Printf("%s:%d -- error opening file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2567 return;
2568 }
2569
2570 if(fDebug>0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
2571
2572 if(strdir && strdir.Length()) gDirectory->cd(strdir);
2573
2574 TList* list = 0;
2575
2576 if(strlist && strlist.Length()){
2577
2578 if(!(list = (TList*) gDirectory->Get(strlist))){
2579 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2580 return;
2581 }
2582 }
2583
2584 if(list) hnResponseSinglePt = (THnSparse*) list->FindObject("fhnResponseSinglePt");
2585 else hnResponseSinglePt = (THnSparse*) gDirectory->Get("fhnResponseSinglePt");
2586
2587
2588 if(!hnResponseSinglePt){
2589 Printf("%s:%d -- error retrieving response matrix fhnResponseSinglePt",(char*)__FILE__,__LINE__);
2590 return;
2591 }
2592
2593 f.Close();
2594
2595
2596 // project 2d histo
2597 h2ResponseSinglePt = (TH2F*) hnResponseSinglePt->Projection(1,0);// note convention: yDim,xDim
2598 h2ResponseSinglePt->SetNameTitle("h2ResponseSinglePt","");
2599
2600
2601 // write
2602
2603 TString outfileOption = "RECREATE";
2604 if(updateOutfile) outfileOption = "UPDATE";
2605
2606 TFile out(strOutfile,outfileOption);
2607
2608 if(!out.IsOpen()){
2609 Printf("%s:%d -- error opening response matrix output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2610 return;
2611 }
2612
2613 if(fDebug>0) Printf("%s:%d -- write response matrices to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2614
2615 if(strOutDir && strOutDir.Length()){
2616
2617 TDirectory* dir;
2618 if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd();
2619 else{
2620 dir = out.mkdir(strOutDir);
2621 dir->cd();
2622 }
2623 }
2624
2625 hnResponseSinglePt->Write();
2626 h2ResponseSinglePt->Write();
2627
2628 out.Close();
2629}
2630
2631//________________________________________________________________________________________________________________
2632void AliFragmentationFunctionCorrections::WriteJetTrackEff(TString strInfile, TString strID, TString strOutfile,
2633 Bool_t updateOutfile)
2634{
2635 // read task ouput from MC and write single track eff - standard dir/list
2636
b541fbca 2637 TString strdir = "PWGJE_FragmentationFunction_" + strID;
39e2b057 2638 TString strlist = "fracfunc_" + strID;
2639
2640 WriteJetTrackEff(strInfile,strdir,strlist,strOutfile,updateOutfile);
2641}
2642
2643//___________________________________________________________________________________________________________________________________
2644void AliFragmentationFunctionCorrections::WriteJetTrackEff(TString strInfile, TString strdir, TString strlist,
2645 TString strOutfile, Bool_t updateOutfile)
2646{
2647 // read task output from MC and write track eff in jet pt slices as function of pt, z, xi
2648 // argument strlist optional - read from directory strdir if not specified
2649 // write eff to file strOutfile - by default only 'update' file (don't overwrite)
2650
2651 TH1F* hEffPt[fNJetPtSlices];
2652 TH1F* hEffXi[fNJetPtSlices];
2653 TH1F* hEffZ[fNJetPtSlices];
2654
2655 TH1F* hdNdptTracksMCPrimGen[fNJetPtSlices];
2656 TH1F* hdNdxiMCPrimGen[fNJetPtSlices];
2657 TH1F* hdNdzMCPrimGen[fNJetPtSlices];
2658
2659 TH1F* hdNdptTracksMCPrimRec[fNJetPtSlices];
2660 TH1F* hdNdxiMCPrimRec[fNJetPtSlices];
2661 TH1F* hdNdzMCPrimRec[fNJetPtSlices];
2662
2663
2664 TH1F* fh1FFJetPtRecEffGen;
2665
2666 TH2F* fh2FFTrackPtRecEffGen;
2667 TH2F* fh2FFZRecEffGen;
2668 TH2F* fh2FFXiRecEffGen;
2669
2670 TH2F* fh2FFTrackPtRecEffRec;
2671 TH2F* fh2FFZRecEffRec;
2672 TH2F* fh2FFXiRecEffRec;
2673
2674
2675 TFile f(strInfile,"READ");
2676
2677 if(!f.IsOpen()){
2678 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2679 return;
2680 }
2681
2682 if(fDebug>0) Printf("%s:%d -- writeJetTrackEff: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
2683
2684 if(strdir && strdir.Length()) gDirectory->cd(strdir);
2685
2686 TList* list = 0;
2687
2688 if(strlist && strlist.Length()){
2689
2690 if(!(list = (TList*) gDirectory->Get(strlist))){
2691 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2692 return;
2693 }
2694 }
2695
2696 if(list){
2697 fh1FFJetPtRecEffGen = (TH1F*) list->FindObject("fh1FFJetPtRecEffGen");
2698
2699 fh2FFTrackPtRecEffGen = (TH2F*) list->FindObject("fh2FFTrackPtRecEffGen");
2700 fh2FFZRecEffGen = (TH2F*) list->FindObject("fh2FFZRecEffGen");
2701 fh2FFXiRecEffGen = (TH2F*) list->FindObject("fh2FFXiRecEffGen");
2702
2703 fh2FFTrackPtRecEffRec = (TH2F*) list->FindObject("fh2FFTrackPtRecEffRec");
2704 fh2FFZRecEffRec = (TH2F*) list->FindObject("fh2FFZRecEffRec");
2705 fh2FFXiRecEffRec = (TH2F*) list->FindObject("fh2FFXiRecEffRec");
2706 }
2707 else{
ce55b926 2708 fh1FFJetPtRecEffGen = (TH1F*) gDirectory->Get("fh1FFJetPtRecEffGen");
39e2b057 2709
ce55b926 2710 fh2FFTrackPtRecEffGen = (TH2F*) gDirectory->Get("fh2FFTrackPtRecEffGen");
2711 fh2FFZRecEffGen = (TH2F*) gDirectory->Get("fh2FFZRecEffGen");
2712 fh2FFXiRecEffGen = (TH2F*) gDirectory->Get("fh2FFXiRecEffGen");
39e2b057 2713
ce55b926 2714 fh2FFTrackPtRecEffRec = (TH2F*) gDirectory->Get("fh2FFTrackPtRecEffRec");
2715 fh2FFZRecEffRec = (TH2F*) gDirectory->Get("fh2FFZRecEffRec");
2716 fh2FFXiRecEffRec = (TH2F*) gDirectory->Get("fh2FFXiRecEffRec");
39e2b057 2717 }
2718
2719 fh1FFJetPtRecEffGen->SetDirectory(0);
2720
2721 fh2FFTrackPtRecEffGen->SetDirectory(0);
2722 fh2FFZRecEffGen->SetDirectory(0);
2723 fh2FFXiRecEffGen->SetDirectory(0);
2724
2725 fh2FFTrackPtRecEffRec->SetDirectory(0);
2726 fh2FFZRecEffRec->SetDirectory(0);
2727 fh2FFXiRecEffRec->SetDirectory(0);
2728
2729 f.Close();
2730
2731
2732 // projections: FF for generated and reconstructed primaries
2733
2734 for(Int_t i=0; i<fNJetPtSlices; i++){
2735
2736 Float_t jetPtLoLim = fJetPtSlices->At(i);
2737 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
2738
2739 Int_t binLo = static_cast<Int_t>(fh2FFTrackPtRecEffGen->GetXaxis()->FindBin(jetPtLoLim));
2740 Int_t binUp = static_cast<Int_t>(fh2FFTrackPtRecEffGen->GetXaxis()->FindBin(jetPtUpLim))-1;
2741
2742 if(binUp > fh2FFTrackPtRecEffGen->GetNbinsX()){
2743 Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
2744 return;
2745 }
2746
2747 TString strNameFFPtGen(Form("fh1FFTrackPtGenPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2748 TString strNameFFZGen(Form("fh1FFZGenPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2749 TString strNameFFXiGen(Form("fh1FFXiGenPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2750
2751 TString strNameFFPtRec(Form("fh1FFTrackPtRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2752 TString strNameFFZRec(Form("fh1FFZRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2753 TString strNameFFXiRec(Form("fh1FFXiRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2754
2755 // project
2756 // appendix 'unbinned' to avoid histos with same name after rebinning
2757
2758 hdNdptTracksMCPrimGen[i] = (TH1F*) fh2FFTrackPtRecEffGen->ProjectionY(strNameFFPtGen+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
2759 hdNdzMCPrimGen[i] = (TH1F*) fh2FFZRecEffGen->ProjectionY(strNameFFZGen+"_unBinned",binLo,binUp,"o");
2760 hdNdxiMCPrimGen[i] = (TH1F*) fh2FFXiRecEffGen->ProjectionY(strNameFFXiGen+"_unBinned",binLo,binUp,"o");
2761
2762 hdNdptTracksMCPrimRec[i] = (TH1F*) fh2FFTrackPtRecEffRec->ProjectionY(strNameFFPtRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
2763 hdNdzMCPrimRec[i] = (TH1F*) fh2FFZRecEffRec->ProjectionY(strNameFFZRec+"_unBinned",binLo,binUp,"o");
2764 hdNdxiMCPrimRec[i] = (TH1F*) fh2FFXiRecEffRec->ProjectionY(strNameFFXiRec+"_unBinned",binLo,binUp,"o");
2765
2766 // rebin
2767
2768 if(fNHistoBinsPt[i]) hdNdptTracksMCPrimGen[i] = (TH1F*) hdNdptTracksMCPrimGen[i]->Rebin(fNHistoBinsPt[i],strNameFFPtGen,fHistoBinsPt[i]->GetArray());
2769 if(fNHistoBinsZ[i]) hdNdzMCPrimGen[i] = (TH1F*) hdNdzMCPrimGen[i]->Rebin(fNHistoBinsZ[i],strNameFFZGen,fHistoBinsZ[i]->GetArray());
2770 if(fNHistoBinsXi[i]) hdNdxiMCPrimGen[i] = (TH1F*) hdNdxiMCPrimGen[i]->Rebin(fNHistoBinsXi[i],strNameFFXiGen,fHistoBinsXi[i]->GetArray());
2771
2772 if(fNHistoBinsPt[i]) hdNdptTracksMCPrimRec[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtRec,fHistoBinsPt[i]->GetArray());
2773 if(fNHistoBinsZ[i]) hdNdzMCPrimRec[i] = (TH1F*) hdNdzMCPrimRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZRec,fHistoBinsZ[i]->GetArray());
2774 if(fNHistoBinsXi[i]) hdNdxiMCPrimRec[i] = (TH1F*) hdNdxiMCPrimRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiRec,fHistoBinsXi[i]->GetArray());
2775
2776 hdNdptTracksMCPrimGen[i]->SetNameTitle(strNameFFPtGen,"");
2777 hdNdzMCPrimGen[i]->SetNameTitle(strNameFFZGen,"");
2778 hdNdxiMCPrimGen[i]->SetNameTitle(strNameFFXiGen,"");
2779
2780 hdNdptTracksMCPrimRec[i]->SetNameTitle(strNameFFPtRec,"");
2781 hdNdzMCPrimRec[i]->SetNameTitle(strNameFFZRec,"");
2782 hdNdxiMCPrimRec[i]->SetNameTitle(strNameFFXiRec,"");
2783
2784 // normalize
2785
2786 Double_t nJetsBin = fh1FFJetPtRecEffGen->Integral(binLo,binUp);
2787
2788 NormalizeTH1(hdNdptTracksMCPrimGen[i],nJetsBin);
2789 NormalizeTH1(hdNdzMCPrimGen[i],nJetsBin);
2790 NormalizeTH1(hdNdxiMCPrimGen[i],nJetsBin);
2791
2792 NormalizeTH1(hdNdptTracksMCPrimRec[i],nJetsBin);
2793 NormalizeTH1(hdNdzMCPrimRec[i],nJetsBin);
2794 NormalizeTH1(hdNdxiMCPrimRec[i],nJetsBin);
2795
2796 // divide rec/gen : efficiency
2797
2798 TString strNameEffPt(Form("hEffPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2799 TString strNameEffZ(Form("hEffZ_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2800 TString strNameEffXi(Form("hEffXi_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2801
2802 hEffPt[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Clone(strNameEffPt);
2803 hEffPt[i]->Divide(hdNdptTracksMCPrimRec[i],hdNdptTracksMCPrimGen[i],1,1,"B"); // binominal errors
2804
2805 hEffXi[i] = (TH1F*) hdNdxiMCPrimRec[i]->Clone(strNameEffXi);
2806 hEffXi[i]->Divide(hdNdxiMCPrimRec[i],hdNdxiMCPrimGen[i],1,1,"B"); // binominal errors
2807
2808 hEffZ[i] = (TH1F*) hdNdzMCPrimRec[i]->Clone(strNameEffZ);
2809 hEffZ[i]->Divide(hdNdzMCPrimRec[i],hdNdzMCPrimGen[i],1,1,"B"); // binominal errors
2810 }
2811
2812 // write
2813
2814 TString outfileOption = "RECREATE";
2815 if(updateOutfile) outfileOption = "UPDATE";
2816
2817 TFile out(strOutfile,outfileOption);
2818
2819 if(!out.IsOpen()){
2820 Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2821 return;
2822 }
2823
2824 if(fDebug>0) Printf("%s:%d -- write efficiency to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2825
2826// if(strdir && strdir.Length()){
2827// TDirectory* dir = out.mkdir(strdir);
2828// dir->cd();
2829// }
2830
2831 for(Int_t i=0; i<fNJetPtSlices; i++){
2832
2833 hdNdptTracksMCPrimGen[i]->Write();
2834 hdNdxiMCPrimGen[i]->Write();
2835 hdNdzMCPrimGen[i]->Write();
2836
2837 hdNdptTracksMCPrimRec[i]->Write();
2838 hdNdxiMCPrimRec[i]->Write();
2839 hdNdzMCPrimRec[i]->Write();
2840
2841 hEffPt[i]->Write();
2842 hEffXi[i]->Write();
2843 hEffZ[i]->Write();
2844 }
2845
2846 out.Close();
2847
ce55b926 2848
2849 delete fh1FFJetPtRecEffGen;
2850
2851 delete fh2FFTrackPtRecEffGen;
2852 delete fh2FFZRecEffGen;
2853 delete fh2FFXiRecEffGen;
2854
2855 delete fh2FFTrackPtRecEffRec;
2856 delete fh2FFZRecEffRec;
2857 delete fh2FFXiRecEffRec;
39e2b057 2858}
2859
2860//________________________________________________________________________________________________________________
2861void AliFragmentationFunctionCorrections::WriteJetSecCorr(TString strInfile, TString strID, TString strOutfile,
ce55b926 2862 Bool_t updateOutfile, TString strOutDir)
2863{
2864 // read task ouput from MC and write secondary correction - standard dir/list
2865
2866 TString strdir = "PWGJE_FragmentationFunction_" + strID;
2867 TString strlist = "fracfunc_" + strID;
2868
2869 WriteJetSecCorr(strInfile,strdir,strlist,strOutfile,updateOutfile, strOutDir,kFALSE,"");
2870}
2871
2872//________________________________________________________________________________________________________________
2873void AliFragmentationFunctionCorrections::WriteBgrJetSecCorr(TString strInfile, TString strBgrID, TString strID, TString strOutfile,
37f35567 2874 Bool_t updateOutfile, TString strOutDir,Double_t scaleFacBgrRec)
39e2b057 2875{
2876 // read task ouput from MC and write secondary correction - standard dir/list
2877
b541fbca 2878 TString strdir = "PWGJE_FragmentationFunction_" + strID;
39e2b057 2879 TString strlist = "fracfunc_" + strID;
2880
37f35567 2881 WriteJetSecCorr(strInfile,strdir,strlist,strOutfile,updateOutfile, strOutDir,kTRUE,strBgrID,scaleFacBgrRec);
39e2b057 2882}
2883
ce55b926 2884
39e2b057 2885//___________________________________________________________________________________________________________________________________
2886void AliFragmentationFunctionCorrections::WriteJetSecCorr(TString strInfile, TString strdir, TString strlist,
ce55b926 2887 TString strOutfile, Bool_t updateOutfile, TString strOutDir,
37f35567 2888 Bool_t writeBgr, TString strBgrID,Double_t scaleFacBgrRec)
39e2b057 2889{
2890 // read task output from MC and write secondary correction in jet pt slices as function of pt, z, xi
2891 // argument strlist optional - read from directory strdir if not specified
2892 // write eff to file strOutfile - by default only 'update' file (don't overwrite)
2893
2894 TH1F* hSecCorrPt[fNJetPtSlices];
2895 TH1F* hSecCorrXi[fNJetPtSlices];
2896 TH1F* hSecCorrZ[fNJetPtSlices];
2897
ce55b926 2898 // corr factors using naive Pythia strangeness - for sys uncertainties
2899 TH1F* hSecCorrPt_nonSc[fNJetPtSlices];
2900 TH1F* hSecCorrXi_nonSc[fNJetPtSlices];
2901 TH1F* hSecCorrZ_nonSc[fNJetPtSlices];
2902
39e2b057 2903 TH1F* hdNdptTracksMCPrimRec[fNJetPtSlices];
2904 TH1F* hdNdxiMCPrimRec[fNJetPtSlices];
2905 TH1F* hdNdzMCPrimRec[fNJetPtSlices];
2906
ce55b926 2907 TH1F* hdNdptTracksMCSecRecNS[fNJetPtSlices];
2908 TH1F* hdNdxiMCSecRecNS[fNJetPtSlices];
2909 TH1F* hdNdzMCSecRecNS[fNJetPtSlices];
39e2b057 2910
ce55b926 2911 TH1F* hdNdptTracksMCSecRecS[fNJetPtSlices];
2912 TH1F* hdNdxiMCSecRecS[fNJetPtSlices];
2913 TH1F* hdNdzMCSecRecS[fNJetPtSlices];
2914
2915 TH1F* hdNdptTracksMCSecRecSsc[fNJetPtSlices];
2916 TH1F* hdNdxiMCSecRecSsc[fNJetPtSlices];
2917 TH1F* hdNdzMCSecRecSsc[fNJetPtSlices];
2918
2919 // ---
2920
2921 TH1F* fh1FFJetPtRecEffRec;
39e2b057 2922
2923 TH2F* fh2FFTrackPtRecEffRec;
2924 TH2F* fh2FFZRecEffRec;
2925 TH2F* fh2FFXiRecEffRec;
2926
ce55b926 2927 TH2F* fh2FFTrackPtSecRecNS;
2928 TH2F* fh2FFZSecRecNS;
2929 TH2F* fh2FFXiSecRecNS;
2930
2931 TH2F* fh2FFTrackPtSecRecS;
2932 TH2F* fh2FFZSecRecS;
2933 TH2F* fh2FFXiSecRecS;
2934
2935 TH2F* fh2FFTrackPtSecRecSsc;
2936 TH2F* fh2FFZSecRecSsc;
2937 TH2F* fh2FFXiSecRecSsc;
2938
2939 // ---
39e2b057 2940
2941 TFile f(strInfile,"READ");
2942
2943 if(!f.IsOpen()){
2944 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2945 return;
2946 }
2947
ce55b926 2948 if(fDebug>0) Printf("%s:%d -- writeJetTrackSecCorr: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
39e2b057 2949
2950 if(strdir && strdir.Length()) gDirectory->cd(strdir);
2951
2952 TList* list = 0;
2953
2954 if(strlist && strlist.Length()){
2955
2956 if(!(list = (TList*) gDirectory->Get(strlist))){
2957 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2958 return;
2959 }
2960 }
2961
2962 if(list){
ce55b926 2963 fh1FFJetPtRecEffRec = (TH1F*) list->FindObject("fh1FFJetPtRecEffRec");
39e2b057 2964
ce55b926 2965 if(writeBgr){
2966 fh2FFTrackPtRecEffRec = (TH2F*) list->FindObject(Form("fh2FFTrackPt%sRecEffRec",strBgrID.Data()));
2967 fh2FFZRecEffRec = (TH2F*) list->FindObject(Form("fh2FFZ%sRecEffRec",strBgrID.Data()));
2968 fh2FFXiRecEffRec = (TH2F*) list->FindObject(Form("fh2FFXi%sRecEffRec",strBgrID.Data()));
2969
2970 fh2FFTrackPtSecRecNS = (TH2F*) list->FindObject(Form("fh2FFTrackPt%sSecRecNS",strBgrID.Data()));
2971 fh2FFZSecRecNS = (TH2F*) list->FindObject(Form("fh2FFZ%sSecRecNS",strBgrID.Data()));
2972 fh2FFXiSecRecNS = (TH2F*) list->FindObject(Form("fh2FFXi%sSecRecNS",strBgrID.Data()));
2973
2974 fh2FFTrackPtSecRecS = (TH2F*) list->FindObject(Form("fh2FFTrackPt%sSecRecS",strBgrID.Data()));
2975 fh2FFZSecRecS = (TH2F*) list->FindObject(Form("fh2FFZ%sSecRecS",strBgrID.Data()));
2976 fh2FFXiSecRecS = (TH2F*) list->FindObject(Form("fh2FFXi%sSecRecS",strBgrID.Data()));
2977
2978 fh2FFTrackPtSecRecSsc = (TH2F*) list->FindObject(Form("fh2FFTrackPt%sSecRecSsc",strBgrID.Data()));
2979 fh2FFZSecRecSsc = (TH2F*) list->FindObject(Form("fh2FFZ%sSecRecSsc",strBgrID.Data()));
2980 fh2FFXiSecRecSsc = (TH2F*) list->FindObject(Form("fh2FFXi%sSecRecSsc",strBgrID.Data()));
2981 }
2982 else{
2983 fh2FFTrackPtRecEffRec = (TH2F*) list->FindObject("fh2FFTrackPtRecEffRec");
2984 fh2FFZRecEffRec = (TH2F*) list->FindObject("fh2FFZRecEffRec");
2985 fh2FFXiRecEffRec = (TH2F*) list->FindObject("fh2FFXiRecEffRec");
2986
2987 fh2FFTrackPtSecRecNS = (TH2F*) list->FindObject("fh2FFTrackPtSecRecNS");
2988 fh2FFZSecRecNS = (TH2F*) list->FindObject("fh2FFZSecRecNS");
2989 fh2FFXiSecRecNS = (TH2F*) list->FindObject("fh2FFXiSecRecNS");
2990
2991 fh2FFTrackPtSecRecS = (TH2F*) list->FindObject("fh2FFTrackPtSecRecS");
2992 fh2FFZSecRecS = (TH2F*) list->FindObject("fh2FFZSecRecS");
2993 fh2FFXiSecRecS = (TH2F*) list->FindObject("fh2FFXiSecRecS");
2994
2995 fh2FFTrackPtSecRecSsc = (TH2F*) list->FindObject("fh2FFTrackPtSecRecSsc");
2996 fh2FFZSecRecSsc = (TH2F*) list->FindObject("fh2FFZSecRecSsc");
2997 fh2FFXiSecRecSsc = (TH2F*) list->FindObject("fh2FFXiSecRecSsc");
2998 }
39e2b057 2999 }
3000 else{
ce55b926 3001 fh1FFJetPtRecEffRec = (TH1F*) gDirectory->Get("fh1FFJetPtRecEffRec");
39e2b057 3002
ce55b926 3003 if(writeBgr){
3004 fh2FFTrackPtRecEffRec = (TH2F*) gDirectory->Get(Form("fh2FFTrackPt%sRecEffRec",strBgrID.Data()));
3005 fh2FFZRecEffRec = (TH2F*) gDirectory->Get(Form("fh2FFZ%sRecEffRec",strBgrID.Data()));
3006 fh2FFXiRecEffRec = (TH2F*) gDirectory->Get(Form("fh2FFXi%sRecEffRec",strBgrID.Data()));
3007
3008 fh2FFTrackPtSecRecNS = (TH2F*) gDirectory->Get(Form("fh2FFTrackPt%sSecRecNS",strBgrID.Data()));
3009 fh2FFZSecRecNS = (TH2F*) gDirectory->Get(Form("fh2FFZ%sSecRecNS",strBgrID.Data()));
3010 fh2FFXiSecRecNS = (TH2F*) gDirectory->Get(Form("fh2FFXi%sSecRecNS",strBgrID.Data()));
3011
3012 fh2FFTrackPtSecRecS = (TH2F*) gDirectory->Get(Form("fh2FFTrackPt%sSecRecS",strBgrID.Data()));
3013 fh2FFZSecRecS = (TH2F*) gDirectory->Get(Form("fh2FFZ%sSecRecS",strBgrID.Data()));
3014 fh2FFXiSecRecS = (TH2F*) gDirectory->Get(Form("fh2FFXi%sSecRecS",strBgrID.Data()));
3015
3016 fh2FFTrackPtSecRecSsc = (TH2F*) gDirectory->Get(Form("fh2FFTrackPt%sSecRecSsc",strBgrID.Data()));
3017 fh2FFZSecRecSsc = (TH2F*) gDirectory->Get(Form("fh2FFZ%sSecRecSsc",strBgrID.Data()));
3018 fh2FFXiSecRecSsc = (TH2F*) gDirectory->Get(Form("fh2FFXi%sSecRecSsc",strBgrID.Data()));
3019 }
3020 else{
3021 fh2FFTrackPtRecEffRec = (TH2F*) gDirectory->Get("fh2FFTrackPtRecEffRec");
3022 fh2FFZRecEffRec = (TH2F*) gDirectory->Get("fh2FFZRecEffRec");
3023 fh2FFXiRecEffRec = (TH2F*) gDirectory->Get("fh2FFXiRecEffRec");
39e2b057 3024
ce55b926 3025 fh2FFTrackPtSecRecNS = (TH2F*) gDirectory->Get("fh2FFTrackPtSecRecNS");
3026 fh2FFZSecRecNS = (TH2F*) gDirectory->Get("fh2FFZSecRecNS");
3027 fh2FFXiSecRecNS = (TH2F*) gDirectory->Get("fh2FFXiSecRecNS");
3028
3029 fh2FFTrackPtSecRecS = (TH2F*) gDirectory->Get("fh2FFTrackPtSecRecS");
3030 fh2FFZSecRecS = (TH2F*) gDirectory->Get("fh2FFZSecRecS");
3031 fh2FFXiSecRecS = (TH2F*) gDirectory->Get("fh2FFXiSecRecS");
3032
3033 fh2FFTrackPtSecRecSsc = (TH2F*) gDirectory->Get("fh2FFTrackPtSecRecSsc");
3034 fh2FFZSecRecSsc = (TH2F*) gDirectory->Get("fh2FFZSecRecSsc");
3035 fh2FFXiSecRecSsc = (TH2F*) gDirectory->Get("fh2FFXiSecRecSsc");
3036 }
39e2b057 3037 }
3038
ce55b926 3039 fh1FFJetPtRecEffRec->SetDirectory(0);
3040
39e2b057 3041 fh2FFTrackPtRecEffRec->SetDirectory(0);
3042 fh2FFZRecEffRec->SetDirectory(0);
3043 fh2FFXiRecEffRec->SetDirectory(0);
3044
ce55b926 3045 fh2FFTrackPtSecRecNS->SetDirectory(0);
3046 fh2FFZSecRecNS->SetDirectory(0);
3047 fh2FFXiSecRecNS->SetDirectory(0);
39e2b057 3048
ce55b926 3049 fh2FFTrackPtSecRecS->SetDirectory(0);
3050 fh2FFZSecRecS->SetDirectory(0);
3051 fh2FFXiSecRecS->SetDirectory(0);
3052
3053 fh2FFTrackPtSecRecSsc->SetDirectory(0);
3054 fh2FFZSecRecSsc->SetDirectory(0);
3055 fh2FFXiSecRecSsc->SetDirectory(0);
3056
39e2b057 3057 f.Close();
3058
3059
3060 // projections: FF for generated and reconstructed primaries
3061
3062 for(Int_t i=0; i<fNJetPtSlices; i++){
3063
3064 Float_t jetPtLoLim = fJetPtSlices->At(i);
3065 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
3066
3067 Int_t binLo = static_cast<Int_t>(fh2FFTrackPtRecEffRec->GetXaxis()->FindBin(jetPtLoLim));
3068 Int_t binUp = static_cast<Int_t>(fh2FFTrackPtRecEffRec->GetXaxis()->FindBin(jetPtUpLim))-1;
3069
3070 if(binUp > fh2FFTrackPtRecEffRec->GetNbinsX()){
3071 Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
3072 return;
3073 }
3074
3075 TString strNameFFPtPrimRec(Form("fh1FFTrackPtRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3076 TString strNameFFZPrimRec(Form("fh1FFZRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3077 TString strNameFFXiPrimRec(Form("fh1FFXiRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3078
ce55b926 3079 TString strNameFFPtSecRecNS(Form("fh1FFTrackPtRecSecNS_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3080 TString strNameFFZSecRecNS(Form("fh1FFZRecSecNS_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3081 TString strNameFFXiSecRecNS(Form("fh1FFXiRecSecNS_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3082
3083 TString strNameFFPtSecRecS(Form("fh1FFTrackPtRecSecS_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3084 TString strNameFFZSecRecS(Form("fh1FFZRecSecS_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3085 TString strNameFFXiSecRecS(Form("fh1FFXiRecSecS_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3086
3087 TString strNameFFPtSecRecSsc(Form("fh1FFTrackPtRecSecSsc_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3088 TString strNameFFZSecRecSsc(Form("fh1FFZRecSecSsc_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3089 TString strNameFFXiSecRecSsc(Form("fh1FFXiRecSecSsc_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3090
3091 if(writeBgr){
3092 strNameFFPtPrimRec.Form("fh1BgrTrackPt%sRecPrim_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3093 strNameFFZPrimRec.Form("fh1BgrZ%sRecPrim_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3094 strNameFFXiPrimRec.Form("fh1BgrXi%sRecPrim_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3095
3096 strNameFFPtSecRecNS.Form("fh1BgrTrackPt%sRecSecNS_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3097 strNameFFZSecRecNS.Form("fh1BgrZ%sRecSecNS_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3098 strNameFFXiSecRecNS.Form("fh1BgrXi%sRecSecNS_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3099
3100 strNameFFPtSecRecS.Form("fh1BgrTrackPt%sRecSecS_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3101 strNameFFZSecRecS.Form("fh1BgrZ%sRecSecS_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3102 strNameFFXiSecRecS.Form("fh1BgrXi%sRecSecS_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3103
3104 strNameFFPtSecRecSsc.Form("fh1BgrTrackPt%sRecSecSsc_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3105 strNameFFZSecRecSsc.Form("fh1BgrZ%sRecSecSsc_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3106 strNameFFXiSecRecSsc.Form("fh1BgrXi%sRecSecSsc_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3107 }
39e2b057 3108
3109 // project
3110 // appendix 'unbinned' to avoid histos with same name after rebinning
3111
3112 hdNdptTracksMCPrimRec[i] = (TH1F*) fh2FFTrackPtRecEffRec->ProjectionY(strNameFFPtPrimRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
3113 hdNdzMCPrimRec[i] = (TH1F*) fh2FFZRecEffRec->ProjectionY(strNameFFZPrimRec+"_unBinned",binLo,binUp,"o");
3114 hdNdxiMCPrimRec[i] = (TH1F*) fh2FFXiRecEffRec->ProjectionY(strNameFFXiPrimRec+"_unBinned",binLo,binUp,"o");
3115
ce55b926 3116 hdNdptTracksMCSecRecNS[i] = (TH1F*) fh2FFTrackPtSecRecNS->ProjectionY(strNameFFPtSecRecNS+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
3117 hdNdzMCSecRecNS[i] = (TH1F*) fh2FFZSecRecNS->ProjectionY(strNameFFZSecRecNS+"_unBinned",binLo,binUp,"o");
3118 hdNdxiMCSecRecNS[i] = (TH1F*) fh2FFXiSecRecNS->ProjectionY(strNameFFXiSecRecNS+"_unBinned",binLo,binUp,"o");
39e2b057 3119
ce55b926 3120 hdNdptTracksMCSecRecS[i] = (TH1F*) fh2FFTrackPtSecRecS->ProjectionY(strNameFFPtSecRecS+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
3121 hdNdzMCSecRecS[i] = (TH1F*) fh2FFZSecRecS->ProjectionY(strNameFFZSecRecS+"_unBinned",binLo,binUp,"o");
3122 hdNdxiMCSecRecS[i] = (TH1F*) fh2FFXiSecRecS->ProjectionY(strNameFFXiSecRecS+"_unBinned",binLo,binUp,"o");
3123
3124 hdNdptTracksMCSecRecSsc[i] = (TH1F*) fh2FFTrackPtSecRecSsc->ProjectionY(strNameFFPtSecRecSsc+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
3125 hdNdzMCSecRecSsc[i] = (TH1F*) fh2FFZSecRecSsc->ProjectionY(strNameFFZSecRecSsc+"_unBinned",binLo,binUp,"o");
3126 hdNdxiMCSecRecSsc[i] = (TH1F*) fh2FFXiSecRecSsc->ProjectionY(strNameFFXiSecRecSsc+"_unBinned",binLo,binUp,"o");
39e2b057 3127
ce55b926 3128
3129 // rebin
39e2b057 3130 if(fNHistoBinsPt[i]) hdNdptTracksMCPrimRec[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtPrimRec,fHistoBinsPt[i]->GetArray());
3131 if(fNHistoBinsZ[i]) hdNdzMCPrimRec[i] = (TH1F*) hdNdzMCPrimRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZPrimRec,fHistoBinsZ[i]->GetArray());
3132 if(fNHistoBinsXi[i]) hdNdxiMCPrimRec[i] = (TH1F*) hdNdxiMCPrimRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiPrimRec,fHistoBinsXi[i]->GetArray());
3133
ce55b926 3134 if(fNHistoBinsPt[i]) hdNdptTracksMCSecRecNS[i] = (TH1F*) hdNdptTracksMCSecRecNS[i]->Rebin(fNHistoBinsPt[i],strNameFFPtSecRecNS,fHistoBinsPt[i]->GetArray());
3135 if(fNHistoBinsZ[i]) hdNdzMCSecRecNS[i] = (TH1F*) hdNdzMCSecRecNS[i]->Rebin(fNHistoBinsZ[i],strNameFFZSecRecNS,fHistoBinsZ[i]->GetArray());
3136 if(fNHistoBinsXi[i]) hdNdxiMCSecRecNS[i] = (TH1F*) hdNdxiMCSecRecNS[i]->Rebin(fNHistoBinsXi[i],strNameFFXiSecRecNS,fHistoBinsXi[i]->GetArray());
3137
3138 if(fNHistoBinsPt[i]) hdNdptTracksMCSecRecS[i] = (TH1F*) hdNdptTracksMCSecRecS[i]->Rebin(fNHistoBinsPt[i],strNameFFPtSecRecS,fHistoBinsPt[i]->GetArray());
3139 if(fNHistoBinsZ[i]) hdNdzMCSecRecS[i] = (TH1F*) hdNdzMCSecRecS[i]->Rebin(fNHistoBinsZ[i],strNameFFZSecRecS,fHistoBinsZ[i]->GetArray());
3140 if(fNHistoBinsXi[i]) hdNdxiMCSecRecS[i] = (TH1F*) hdNdxiMCSecRecS[i]->Rebin(fNHistoBinsXi[i],strNameFFXiSecRecS,fHistoBinsXi[i]->GetArray());
3141
3142 if(fNHistoBinsPt[i]) hdNdptTracksMCSecRecSsc[i] = (TH1F*) hdNdptTracksMCSecRecSsc[i]->Rebin(fNHistoBinsPt[i],strNameFFPtSecRecSsc,fHistoBinsPt[i]->GetArray());
3143 if(fNHistoBinsZ[i]) hdNdzMCSecRecSsc[i] = (TH1F*) hdNdzMCSecRecSsc[i]->Rebin(fNHistoBinsZ[i],strNameFFZSecRecSsc,fHistoBinsZ[i]->GetArray());
3144 if(fNHistoBinsXi[i]) hdNdxiMCSecRecSsc[i] = (TH1F*) hdNdxiMCSecRecSsc[i]->Rebin(fNHistoBinsXi[i],strNameFFXiSecRecSsc,fHistoBinsXi[i]->GetArray());
3145
39e2b057 3146
3147 hdNdptTracksMCPrimRec[i]->SetNameTitle(strNameFFPtPrimRec,"");
3148 hdNdzMCPrimRec[i]->SetNameTitle(strNameFFZPrimRec,"");
3149 hdNdxiMCPrimRec[i]->SetNameTitle(strNameFFXiPrimRec,"");
3150
ce55b926 3151 hdNdptTracksMCSecRecNS[i]->SetNameTitle(strNameFFPtSecRecNS,"");
3152 hdNdzMCSecRecNS[i]->SetNameTitle(strNameFFZSecRecNS,"");
3153 hdNdxiMCSecRecNS[i]->SetNameTitle(strNameFFXiSecRecNS,"");
39e2b057 3154
ce55b926 3155 hdNdptTracksMCSecRecS[i]->SetNameTitle(strNameFFPtSecRecS,"");
3156 hdNdzMCSecRecS[i]->SetNameTitle(strNameFFZSecRecS,"");
3157 hdNdxiMCSecRecS[i]->SetNameTitle(strNameFFXiSecRecS,"");
3158
3159 hdNdptTracksMCSecRecSsc[i]->SetNameTitle(strNameFFPtSecRecSsc,"");
3160 hdNdzMCSecRecSsc[i]->SetNameTitle(strNameFFZSecRecSsc,"");
3161 hdNdxiMCSecRecSsc[i]->SetNameTitle(strNameFFXiSecRecSsc,"");
39e2b057 3162
37f35567 3163 // normalize
ce55b926 3164 Double_t nJetsBin = fh1FFJetPtRecEffRec->Integral(binLo,binUp);
39e2b057 3165
37f35567 3166 // scale fac for perp2 bgr
3167 if(writeBgr && scaleFacBgrRec && (scaleFacBgrRec != 1)) nJetsBin /= scaleFacBgrRec;
3168
39e2b057 3169 NormalizeTH1(hdNdptTracksMCPrimRec[i],nJetsBin);
3170 NormalizeTH1(hdNdzMCPrimRec[i],nJetsBin);
3171 NormalizeTH1(hdNdxiMCPrimRec[i],nJetsBin);
3172
ce55b926 3173 NormalizeTH1(hdNdptTracksMCSecRecNS[i],nJetsBin);
3174 NormalizeTH1(hdNdzMCSecRecNS[i],nJetsBin);
3175 NormalizeTH1(hdNdxiMCSecRecNS[i],nJetsBin);
39e2b057 3176
ce55b926 3177 NormalizeTH1(hdNdptTracksMCSecRecS[i],nJetsBin);
3178 NormalizeTH1(hdNdzMCSecRecS[i],nJetsBin);
3179 NormalizeTH1(hdNdxiMCSecRecS[i],nJetsBin);
3180
3181 NormalizeTH1(hdNdptTracksMCSecRecSsc[i],nJetsBin);
3182 NormalizeTH1(hdNdzMCSecRecSsc[i],nJetsBin);
3183 NormalizeTH1(hdNdxiMCSecRecSsc[i],nJetsBin);
3184
3185
3186 // divide prim / (prim+sec) : corr factor
39e2b057 3187 TString strNameSecCorrPt(Form("hSecCorrPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3188 TString strNameSecCorrZ(Form("hSecCorrZ_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3189 TString strNameSecCorrXi(Form("hSecCorrXi_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
ce55b926 3190
3191 if(writeBgr){
3192 strNameSecCorrPt.Form("hSecCorrBgrPt_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
3193 strNameSecCorrZ.Form("hSecCorrBgrZ_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
3194 strNameSecCorrXi.Form("hSecCorrBgrXi_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
3195 }
3196
3197 hSecCorrPt[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Clone(strNameSecCorrPt);
3198 TH1F* hSumPrimSecPt = (TH1F*) hdNdptTracksMCPrimRec[i]->Clone("hSumPrimSecPt");
3199 hSumPrimSecPt->Add(hdNdptTracksMCSecRecNS[i]);
3200 hSumPrimSecPt->Add(hdNdptTracksMCSecRecSsc[i]);
39e2b057 3201 hSecCorrPt[i]->Divide(hdNdptTracksMCPrimRec[i],hSumPrimSecPt,1,1,"B"); // binominal errors
3202
ce55b926 3203 hSecCorrXi[i] = (TH1F*) hdNdxiMCPrimRec[i]->Clone(strNameSecCorrXi);
3204 TH1F* hSumPrimSecXi = (TH1F*) hdNdxiMCPrimRec[i]->Clone("hSumPrimSecXi");
3205 hSumPrimSecXi->Add(hdNdxiMCSecRecNS[i]);
3206 hSumPrimSecXi->Add(hdNdxiMCSecRecSsc[i]);
39e2b057 3207 hSecCorrXi[i]->Divide(hdNdxiMCPrimRec[i],hSumPrimSecXi,1,1,"B"); // binominal errors
3208
ce55b926 3209 hSecCorrZ[i] = (TH1F*) hdNdzMCPrimRec[i]->Clone(strNameSecCorrZ);