]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/AliFragmentationFunctionCorrections.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Coverity fixes.
[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"
32
33///////////////////////////////////////////////////////////////////////////////
34// //
35// correction class for ouput of AliAnalysisTaskFragmentationFunction //
36// //
37// corrections for: reconstruction efficiency, momentum resolution, //
38// secondaries, UE / HI background, fluctuations //
39// back-correction on jet energy on dN/dxi //
40// //
41// read MC ouput and write out efficiency histos, response matrices etc. //
42// read measured distributions and apply efficiency, response matrices etc. //
43// //
44// contact: o.busch@gsi.de //
45// //
46///////////////////////////////////////////////////////////////////////////////
47
48
49ClassImp(AliFragmentationFunctionCorrections)
50
51//________________________________________________________________________
52AliFragmentationFunctionCorrections::AliFragmentationFunctionCorrections()
53 : TObject()
54 ,fDebug(0)
55 ,fNJetPtSlices(0)
56 ,fJetPtSlices(0)
57 ,fNJets(0)
58 ,fNJetsBgr(0)
59 ,fNHistoBinsSinglePt(0)
60 ,fHistoBinsSinglePt(0)
61 ,fNHistoBinsPt(0)
62 ,fNHistoBinsZ(0)
63 ,fNHistoBinsXi(0)
64 ,fHistoBinsPt(0)
65 ,fHistoBinsZ(0)
66 ,fHistoBinsXi(0)
67 ,fNCorrectionLevels(0)
68 ,fCorrFF(0)
69 ,fNCorrectionLevelsBgr(0)
70 ,fCorrBgr(0)
71 ,fNCorrectionLevelsSinglePt(0)
72 ,fCorrSinglePt(0)
73 ,fh1FFXiShift(0)
74 ,fh1EffSinglePt(0)
75 ,fh1EffPt(0)
76 ,fh1EffZ(0)
77 ,fh1EffXi(0)
78 ,fh1EffBgrPt(0)
79 ,fh1EffBgrZ(0)
80 ,fh1EffBgrXi(0)
81 ,fh1FFTrackPtBackFolded(0)
82 ,fh1FFZBackFolded(0)
83 ,fh1FFXiBackFolded(0)
84 ,fh1FFRatioTrackPtFolded(0)
85 ,fh1FFRatioZFolded(0)
86 ,fh1FFRatioXiFolded(0)
87 ,fh1FFRatioTrackPtBackFolded(0)
88 ,fh1FFRatioZBackFolded(0)
89 ,fh1FFRatioXiBackFolded(0)
90 ,fh1SingleTrackPtBackFolded(0)
91 ,fh1RatioSingleTrackPtFolded(0)
92 ,fh1RatioSingleTrackPtBackFolded(0)
93 ,fhnResponseSinglePt(0)
94 ,fhnResponsePt(0)
95 ,fhnResponseZ(0)
96 ,fhnResponseXi(0)
97 ,fh1FFTrackPtPrior(0)
98 ,fh1FFZPrior(0)
99 ,fh1FFXiPrior(0)
100 ,fh1SecCorrSinglePt(0)
101{
102 // default constructor
103}
104
105//________________________________________________________________________________________________________________________
106AliFragmentationFunctionCorrections::AliFragmentationFunctionCorrections(const AliFragmentationFunctionCorrections &copy)
107 : TObject()
108 ,fDebug(copy.fDebug)
109 ,fNJetPtSlices(copy.fNJetPtSlices)
110 ,fJetPtSlices(copy.fJetPtSlices)
111 ,fNJets(copy.fNJets)
112 ,fNJetsBgr(copy.fNJetsBgr)
1aa4f09f 113 ,fNHistoBinsSinglePt(copy.fNHistoBinsSinglePt)
114 ,fHistoBinsSinglePt(copy.fHistoBinsSinglePt)
39e2b057 115 ,fNHistoBinsPt(copy.fNHistoBinsPt)
116 ,fNHistoBinsZ(copy.fNHistoBinsZ)
117 ,fNHistoBinsXi(copy.fNHistoBinsXi)
118 ,fHistoBinsPt(copy.fHistoBinsPt)
119 ,fHistoBinsZ(copy.fHistoBinsZ)
120 ,fHistoBinsXi(copy.fHistoBinsXi)
121 ,fNCorrectionLevels(copy.fNCorrectionLevels)
122 ,fCorrFF(copy.fCorrFF)
123 ,fNCorrectionLevelsBgr(copy.fNCorrectionLevelsBgr)
124 ,fCorrBgr(copy.fCorrBgr)
125 ,fNCorrectionLevelsSinglePt(copy.fNCorrectionLevelsSinglePt)
126 ,fCorrSinglePt(copy.fCorrSinglePt)
127 ,fh1FFXiShift(copy.fh1FFXiShift)
1aa4f09f 128 ,fh1EffSinglePt(copy.fh1EffSinglePt)
39e2b057 129 ,fh1EffPt(copy.fh1EffPt)
130 ,fh1EffZ(copy.fh1EffZ)
131 ,fh1EffXi(copy.fh1EffXi)
132 ,fh1EffBgrPt(copy.fh1EffBgrPt)
133 ,fh1EffBgrZ(copy.fh1EffBgrZ)
134 ,fh1EffBgrXi(copy.fh1EffBgrXi)
135 ,fh1FFTrackPtBackFolded(copy.fh1FFTrackPtBackFolded)
136 ,fh1FFZBackFolded(copy.fh1FFZBackFolded)
137 ,fh1FFXiBackFolded(copy.fh1FFXiBackFolded)
138 ,fh1FFRatioTrackPtFolded(copy.fh1FFRatioTrackPtFolded)
139 ,fh1FFRatioZFolded(copy.fh1FFRatioZFolded)
140 ,fh1FFRatioXiFolded(copy.fh1FFRatioXiFolded)
141 ,fh1FFRatioTrackPtBackFolded(copy.fh1FFRatioTrackPtBackFolded)
142 ,fh1FFRatioZBackFolded(copy.fh1FFRatioZBackFolded)
143 ,fh1FFRatioXiBackFolded(copy.fh1FFRatioXiBackFolded)
144 ,fh1SingleTrackPtBackFolded(copy.fh1SingleTrackPtBackFolded)
145 ,fh1RatioSingleTrackPtFolded(copy.fh1RatioSingleTrackPtFolded)
146 ,fh1RatioSingleTrackPtBackFolded(copy.fh1RatioSingleTrackPtBackFolded)
147 ,fhnResponseSinglePt(copy.fhnResponseSinglePt)
148 ,fhnResponsePt(copy.fhnResponsePt)
149 ,fhnResponseZ(copy.fhnResponseZ)
150 ,fhnResponseXi(copy.fhnResponseXi)
151 ,fh1FFTrackPtPrior(copy.fh1FFTrackPtPrior)
152 ,fh1FFZPrior(copy.fh1FFZPrior)
153 ,fh1FFXiPrior(copy.fh1FFXiPrior)
154 ,fh1SecCorrSinglePt(copy.fh1SecCorrSinglePt)
155{
156 // copy constructor
157
158}
159
160// ______________________________________________________________________________________________________________________________
161AliFragmentationFunctionCorrections& AliFragmentationFunctionCorrections::operator=(const AliFragmentationFunctionCorrections& o)
162{
163 // assignment
164
165 if(this!=&o){
166 TObject::operator=(o);
167 fDebug = o.fDebug;
168 fNJetPtSlices = o.fNJetPtSlices;
169 fJetPtSlices = o.fJetPtSlices;
170 fNJets = o.fNJets;
171 fNJetsBgr = o.fNJetsBgr;
172 fNHistoBinsSinglePt = o.fNHistoBinsSinglePt;
173 fHistoBinsSinglePt = o.fHistoBinsSinglePt;
174 fNHistoBinsPt = o.fNHistoBinsPt;
175 fNHistoBinsZ = o.fNHistoBinsZ;
176 fNHistoBinsXi = o.fNHistoBinsXi;
177 fHistoBinsPt = o.fHistoBinsPt;
178 fHistoBinsZ = o.fHistoBinsZ;
179 fHistoBinsXi = o.fHistoBinsXi;
180 fNCorrectionLevels = o.fNCorrectionLevels;
181 fCorrFF = o.fCorrFF;
182 fNCorrectionLevelsBgr = o.fNCorrectionLevelsBgr;
183 fCorrBgr = o.fCorrBgr;
184 fNCorrectionLevelsSinglePt = o.fNCorrectionLevelsSinglePt;
185 fCorrSinglePt = o.fCorrSinglePt;
186 fh1FFXiShift = o.fh1FFXiShift;
187 fh1EffSinglePt = o.fh1EffSinglePt;
188 fh1EffPt = o.fh1EffPt;
189 fh1EffZ = o.fh1EffZ;
190 fh1EffXi = o.fh1EffXi;
191 fh1EffBgrPt = o.fh1EffBgrPt;
192 fh1EffBgrZ = o.fh1EffBgrZ;
193 fh1EffBgrXi = o.fh1EffBgrXi;
194 fh1FFTrackPtBackFolded = o.fh1FFTrackPtBackFolded;
195 fh1FFZBackFolded = o.fh1FFZBackFolded;
196 fh1FFXiBackFolded = o.fh1FFXiBackFolded;
197 fh1FFRatioTrackPtFolded = o.fh1FFRatioTrackPtFolded;
198 fh1FFRatioZFolded = o.fh1FFRatioZFolded;
199 fh1FFRatioXiFolded = o.fh1FFRatioXiFolded;
200 fh1FFRatioTrackPtBackFolded = o.fh1FFRatioTrackPtBackFolded;
201 fh1FFRatioZBackFolded = o.fh1FFRatioZBackFolded;
202 fh1FFRatioXiBackFolded = o.fh1FFRatioXiBackFolded;
203 fh1SingleTrackPtBackFolded = o.fh1SingleTrackPtBackFolded;
204 fh1RatioSingleTrackPtFolded = o.fh1RatioSingleTrackPtFolded;
205 fh1RatioSingleTrackPtBackFolded = o.fh1RatioSingleTrackPtBackFolded;
206 fhnResponseSinglePt = o.fhnResponseSinglePt;
207 fhnResponsePt = o.fhnResponsePt;
208 fhnResponseZ = o.fhnResponseZ;
209 fhnResponseXi = o.fhnResponseXi;
210 fh1FFTrackPtPrior = o.fh1FFTrackPtPrior;
211 fh1FFZPrior = o.fh1FFZPrior;
212 fh1FFXiPrior = o.fh1FFXiPrior;
213 fh1SecCorrSinglePt = o.fh1SecCorrSinglePt;
214 }
215
216 return *this;
217}
218
219//_________________________________________________________________________
220AliFragmentationFunctionCorrections::~AliFragmentationFunctionCorrections()
221{
222 // destructor
223
224 if(fJetPtSlices) delete fJetPtSlices;
225 if(fNJets) delete fNJets;
226 if(fNJetsBgr) delete fNJetsBgr;
227
228 DeleteHistoArray(fh1FFXiShift);
229
230 DeleteHistoArray(fh1EffPt);
231 DeleteHistoArray(fh1EffXi);
232 DeleteHistoArray(fh1EffZ );
233
234 DeleteHistoArray(fh1EffBgrPt);
235 DeleteHistoArray(fh1EffBgrXi);
236 DeleteHistoArray(fh1EffBgrZ);
237
238 // unfolding
239
240 DeleteHistoArray(fh1FFTrackPtBackFolded);
241 DeleteHistoArray(fh1FFZBackFolded);
242 DeleteHistoArray(fh1FFXiBackFolded);
243
244 DeleteHistoArray(fh1FFRatioTrackPtFolded);
245 DeleteHistoArray(fh1FFRatioZFolded);
246 DeleteHistoArray(fh1FFRatioXiFolded);
247
248 DeleteHistoArray(fh1FFRatioTrackPtBackFolded);
249 DeleteHistoArray(fh1FFRatioZBackFolded);
250 DeleteHistoArray(fh1FFRatioXiBackFolded);
251
252 DeleteTHnSparseArray(fhnResponsePt);
253 DeleteTHnSparseArray(fhnResponseZ);
254 DeleteTHnSparseArray(fhnResponseXi);
255
256 DeleteHistoArray(fh1FFTrackPtPrior);
257 DeleteHistoArray(fh1FFZPrior);
258 DeleteHistoArray(fh1FFXiPrior);
259
260
261 // clean up corrected FF
262
263 for(Int_t c=0; c<fNCorrectionLevels; c++) delete fCorrFF[c];
264 delete[] fCorrFF;
265
266 // clean up bgr
267
268 for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) delete fCorrBgr[c];
269 delete[] fCorrBgr;
270
271 // clean up inclusive pt
272 for(Int_t c=0; c<fNCorrectionLevelsSinglePt; c++) delete fCorrSinglePt[c];
273 delete[] fCorrSinglePt;
274
275 delete[] fNHistoBinsPt;
276 delete[] fNHistoBinsZ;
277 delete[] fNHistoBinsXi;
278
279 // clean up histo bins
280
281 if(fHistoBinsSinglePt) delete fHistoBinsSinglePt;
282
283 for(Int_t i=0; i<fNJetPtSlices; i++) delete fHistoBinsPt[i];
284 for(Int_t i=0; i<fNJetPtSlices; i++) delete fHistoBinsZ[i];
285 for(Int_t i=0; i<fNJetPtSlices; i++) delete fHistoBinsXi[i];
286
287 delete[] fHistoBinsPt;
288 delete[] fHistoBinsZ;
289 delete[] fHistoBinsXi;
290}
291
292//_________________________________________________________________________________
293AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::AliFragFuncCorrHistos()
294 : TObject()
295 ,fArraySize(0)
296 ,fh1CorrFFTrackPt(0)
297 ,fh1CorrFFZ(0)
298 ,fh1CorrFFXi(0)
299 ,fCorrLabel(0)
300{
301 // default constructor
302
303}
304
305//__________________________________________________________________________________________________________________
6daac008 306AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::AliFragFuncCorrHistos(const AliFragmentationFunctionCorrections::AliFragFuncCorrHistos& copy)
39e2b057 307 : TObject()
308 ,fArraySize(copy.fArraySize)
6daac008 309 ,fh1CorrFFTrackPt(0)
310 ,fh1CorrFFZ(0)
311 ,fh1CorrFFXi(0)
39e2b057 312 ,fCorrLabel(copy.fCorrLabel)
313{
314 // copy constructor
6daac008 315
316 fh1CorrFFTrackPt = new TH1F*[copy.fArraySize];
317 fh1CorrFFZ = new TH1F*[copy.fArraySize];
318 fh1CorrFFXi = new TH1F*[copy.fArraySize];
319
1aa4f09f 320 for(Int_t i=0; i<copy.fArraySize; i++){
6daac008 321 fh1CorrFFTrackPt[i] = new TH1F(*(copy.fh1CorrFFTrackPt[i]));
322 fh1CorrFFZ[i] = new TH1F(*(copy.fh1CorrFFZ[i]));
323 fh1CorrFFXi[i] = new TH1F(*(copy.fh1CorrFFXi[i]));
1aa4f09f 324 }
39e2b057 325}
326
39e2b057 327//_______________________________________________________________________________________________________________________________________________________________
328AliFragmentationFunctionCorrections::AliFragFuncCorrHistos& AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::operator=(const AliFragmentationFunctionCorrections::AliFragFuncCorrHistos& o)
329{
330 // assignment
331
332 if(this!=&o){
333 TObject::operator=(o);
6daac008 334 Int_t fArraySize_new = o.fArraySize;
39e2b057 335 fCorrLabel = o.fCorrLabel;
1aa4f09f 336
6daac008 337 TH1F** fh1CorrFFTrackPt_new = new TH1F*[fArraySize_new];
338 TH1F** fh1CorrFFZ_new = new TH1F*[fArraySize_new];
339 TH1F** fh1CorrFFXi_new = new TH1F*[fArraySize_new];
340
341 for(Int_t i=0; i<fArraySize_new; i++){
342 fh1CorrFFTrackPt_new[i] = new TH1F(*(o.fh1CorrFFTrackPt[i]));
343 fh1CorrFFZ_new[i] = new TH1F(*(o.fh1CorrFFZ[i]));
344 fh1CorrFFXi_new[i] = new TH1F(*(o.fh1CorrFFXi[i]));
1aa4f09f 345 }
39e2b057 346
6daac008 347 // ---
348
349 if(fArraySize){
350 for(int i=0; i<fArraySize; i++) delete fh1CorrFFTrackPt[i];
351 for(int i=0; i<fArraySize; i++) delete fh1CorrFFZ[i];
352 for(int i=0; i<fArraySize; i++) delete fh1CorrFFXi[i];
353 }
354
355 if(fh1CorrFFTrackPt) delete[] fh1CorrFFTrackPt;
356 if(fh1CorrFFZ) delete[] fh1CorrFFZ;
357 if(fh1CorrFFXi) delete[] fh1CorrFFXi;
358
359 // ---
360
361 fArraySize = fArraySize_new;
362
363 fh1CorrFFTrackPt = fh1CorrFFTrackPt_new;
364 fh1CorrFFZ = fh1CorrFFZ_new;
365 fh1CorrFFXi = fh1CorrFFXi_new;
366
367 for(Int_t i=0; i<fArraySize; i++){
368 fh1CorrFFTrackPt[i] = fh1CorrFFTrackPt_new[i];
369 fh1CorrFFZ[i] = fh1CorrFFZ_new[i];
370 fh1CorrFFXi[i] = fh1CorrFFXi_new[i];
371 }
372 }
373
39e2b057 374 return *this;
375}
376
377//__________________________________________________________________________________
378AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::~AliFragFuncCorrHistos()
379{
380 // destructor
381
382 if(fArraySize){
383 for(int i=0; i<fArraySize; i++) delete fh1CorrFFTrackPt[i];
384 for(int i=0; i<fArraySize; i++) delete fh1CorrFFZ[i];
385 for(int i=0; i<fArraySize; i++) delete fh1CorrFFXi[i];
386 }
387
388 if(fh1CorrFFTrackPt) delete[] fh1CorrFFTrackPt;
389 if(fh1CorrFFZ) delete[] fh1CorrFFZ;
390 if(fh1CorrFFXi) delete[] fh1CorrFFXi;
391
392}
393
394//___________________________________________________________________________________________________________________
395AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::AliFragFuncCorrHistos(const char* label, Int_t arraySize)
396 : TObject()
397 ,fArraySize(0)
398 ,fh1CorrFFTrackPt(0)
399 ,fh1CorrFFZ(0)
400 ,fh1CorrFFXi(0)
401 ,fCorrLabel(label)
402{
403 // constructor
404
405 fArraySize = arraySize;
406 fh1CorrFFTrackPt = new TH1F*[arraySize];
407 fh1CorrFFZ = new TH1F*[arraySize];
408 fh1CorrFFXi = new TH1F*[arraySize];
409
410 for(int i=0; i<arraySize; i++) fh1CorrFFTrackPt[i] = new TH1F;
411 for(int i=0; i<arraySize; i++) fh1CorrFFZ[i] = new TH1F;
412 for(int i=0; i<arraySize; i++) fh1CorrFFXi[i] = new TH1F;
413}
414
415//_______________________________________________________________________________________________________________________________
416void AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::AddCorrHistos(Int_t slice,TH1F* histPt,TH1F* histZ,TH1F* histXi)
417{
418 // place histo in array, add corrLabel to histo name
419
420 if(slice>=fArraySize){
421 Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
422 return;
423 }
424
425 if(histPt){
426 TString name = histPt->GetName();
427 if(fCorrLabel.Length()>0) name += "_"+fCorrLabel;
428 histPt->SetName(name);
429
430 if(!(histPt->GetSumw2())) histPt->Sumw2();
431
432 new(fh1CorrFFTrackPt[slice]) TH1F(*histPt);
433 }
434
435 if(histZ){
436 TString name = histZ->GetName();
437 if(fCorrLabel.Length()>0) name += "_"+fCorrLabel;
438 histZ->SetName(name);
439
440 if(!(histZ->GetSumw2())) histZ->Sumw2();
441
442 new(fh1CorrFFZ[slice]) TH1F(*histZ);
443 }
444
445 if(histXi){
446 TString name = histXi->GetName();
447 if(fCorrLabel.Length()>0) name += "_"+fCorrLabel;
448 histXi->SetName(name);
449
450 if(!(histXi->GetSumw2())) histXi->Sumw2();
451
452 new(fh1CorrFFXi[slice]) TH1F(*histXi);
453 }
454}
455
456//___________________________________________________________________________________________________________________________________
457void AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::ReplaceCorrHistos(Int_t slice,TH1F* histPt,TH1F* histZ,TH1F* histXi)
458{
459 // replace histo in array
460
461 if(slice>=fArraySize){
462 Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
463 return;
464 }
465
466 if(histPt){
467 if(!(histPt->GetSumw2())) histPt->Sumw2();
468
469 TString name = histPt->GetName();
470 histPt->SetName(name);
471
472 delete fh1CorrFFTrackPt[slice];
473 fh1CorrFFTrackPt[slice] = new TH1F;
474 new(fh1CorrFFTrackPt[slice]) TH1F(*histPt);
475 }
476
477 if(histZ){
478 if(!(histZ->GetSumw2())) histZ->Sumw2();
479
480 TString name = histZ->GetName();
481 histZ->SetName(name);
482
483 delete fh1CorrFFZ[slice];
484 fh1CorrFFZ[slice] = new TH1F;
485 new(fh1CorrFFZ[slice]) TH1F(*histZ);
486 }
487
488 if(histXi){
489 if(!(histXi->GetSumw2())) histXi->Sumw2();
490
491 TString name = histXi->GetName();
492 histXi->SetName(name);
493
494 delete fh1CorrFFXi[slice];
495 fh1CorrFFXi[slice] = new TH1F;
496 new(fh1CorrFFXi[slice]) TH1F(*histXi);
497 }
498}
499
500// ___________________________________________________________________________________________
501TH1F* AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::GetTrackPt(const Int_t slice)
502{
503 // return pt histo
504
505 if(slice>=fArraySize){
506 Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
507 return 0;
508 }
509
510 return fh1CorrFFTrackPt[slice];
511
512}
513
514// ______________________________________________________________________________________
515TH1F* AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::GetZ(const Int_t slice)
516{
517 // return z histo
518
519 if(slice>=fArraySize){
520 Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
521 return 0;
522 }
523
524 return fh1CorrFFZ[slice];
525}
526
527// ________________________________________________________________________________________
528TH1F* AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::GetXi(const Int_t slice)
529{
530 // return xi histo
531
532 if(slice>=fArraySize){
533 Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
534 return 0;
535 }
536
537 return fh1CorrFFXi[slice];
538}
539
540// __________________________________________________________________________
541void AliFragmentationFunctionCorrections::DeleteHistoArray(TH1F** hist) const
542{
543 // delete array of TH1
544
545 if(hist){
546 for(Int_t i=0; i<fNJetPtSlices; i++) delete hist[i];
547 delete[] hist;
548 }
549}
550
551// ____________________________________________________________________________________
552void AliFragmentationFunctionCorrections::DeleteTHnSparseArray(THnSparse** hist) const
553{
554 // delete array of THnSparse
555
556 if(hist){
557 for(Int_t i=0; i<fNJetPtSlices; i++) delete hist[i];
558 delete[] hist;
559 }
560}
561
562// _________________________________________________________
563TH1F** AliFragmentationFunctionCorrections::BookHistoArray()
564{
565 // book array of TH1
566
567 if(!fNJetPtSlices){
568 Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
569 return 0;
570 }
571
572 TH1F** hist = new TH1F*[fNJetPtSlices];
573 for(Int_t i=0; i<fNJetPtSlices; i++) hist[i] = new TH1F();
574
575 return hist;
576}
577
578//__________________________________________________________________
579THnSparse** AliFragmentationFunctionCorrections::BookTHnSparseArray()
580{
581 // book array of THnSparse
582
583 if(!fNJetPtSlices){
584 Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
585 return 0;
586 }
587
588 THnSparse** hist = new THnSparse*[fNJetPtSlices];
589 for(Int_t i=0; i<fNJetPtSlices; i++) hist[i] = new THnSparseF();
590
591 return hist;
592}
593
594//_____________________________________________________________________________
595void AliFragmentationFunctionCorrections::AddCorrectionLevel(const char* label)
596{
597 // increase corr level
598
599 if(!fNJetPtSlices){
600 Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
601 return;
602 }
603
604 if(fNCorrectionLevels >= fgMaxNCorrectionLevels){
605 Printf("%s:%d -- max correction level exceeded", (char*)__FILE__,__LINE__);
606 return;
607 }
608
609 fCorrFF[fNCorrectionLevels] = new AliFragFuncCorrHistos(label,fNJetPtSlices);
610 fNCorrectionLevels++;
611}
612
613
614//________________________________________________________________________________
615void AliFragmentationFunctionCorrections::AddCorrectionLevelBgr(const char* label)
616{
617 // increase corr level for bgr FF
618
619 if(!fNJetPtSlices){
620 if(fDebug>0) Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
621 return;
622 }
623
624 if(fNCorrectionLevelsBgr >= fgMaxNCorrectionLevels){
625 Printf("%s:%d -- max correction level exceeded", (char*)__FILE__,__LINE__);
626 return;
627 }
628
629 fCorrBgr[fNCorrectionLevelsBgr] = new AliFragFuncCorrHistos(label,fNJetPtSlices);
630 fNCorrectionLevelsBgr++;
631}
632
633//_____________________________________________________________________________________
634void AliFragmentationFunctionCorrections::AddCorrectionLevelSinglePt(const char* label)
635{
636 // increase corr level single pt spec
637
638 Int_t nJetPtSlicesSingle = 1; // pro forma
639
640 if(fNCorrectionLevelsSinglePt >= fgMaxNCorrectionLevels){
641 Printf("%s:%d -- max correction level exceeded", (char*)__FILE__,__LINE__);
642 return;
643 }
644
645 fCorrSinglePt[fNCorrectionLevelsSinglePt] = new AliFragFuncCorrHistos(label,nJetPtSlicesSingle);
646 fNCorrectionLevelsSinglePt++;
647}
648
649//_____________________________________________________________________________________________
650void AliFragmentationFunctionCorrections::SetJetPtSlices(Float_t* bins, const Int_t nJetPtSlices)
651{
652 // set jet pt slices
653 // once slices are known, can book all other histos
654
655 fNJetPtSlices = nJetPtSlices;
656
657 const Int_t size = nJetPtSlices+1;
658 fJetPtSlices = new TArrayF(size,bins);
659
660 // nJets array
661
662 fNJets = new TArrayF(size);
663 fNJets->Reset(0);
664
665 fNJetsBgr = new TArrayF(size);
666 fNJetsBgr->Reset(0);
667
668 // histos
669
670 fNCorrectionLevels = 0;
671 fCorrFF = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels];
672 AddCorrectionLevel(); // first 'correction' level = raw FF
673
674 fNCorrectionLevelsBgr = 0;
675 fCorrBgr = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels];
676 AddCorrectionLevelBgr(); // first 'correction' level = raw bgr dist
677
678 fh1FFXiShift = BookHistoArray();
679
680 // eff histos
681
682 fh1EffPt = BookHistoArray();
683 fh1EffXi = BookHistoArray();
684 fh1EffZ = BookHistoArray();
685
686 fh1EffBgrPt = BookHistoArray();
687 fh1EffBgrXi = BookHistoArray();
688 fh1EffBgrZ = BookHistoArray();
689
690
691 // unfolding
692
693 fh1FFTrackPtBackFolded = BookHistoArray();
694 fh1FFXiBackFolded = BookHistoArray();
695 fh1FFZBackFolded = BookHistoArray();
696
697 fh1FFRatioTrackPtFolded = BookHistoArray();
698 fh1FFRatioXiFolded = BookHistoArray();
699 fh1FFRatioZFolded = BookHistoArray();
700
701 fh1FFRatioTrackPtBackFolded = BookHistoArray();
702 fh1FFRatioXiBackFolded = BookHistoArray();
703 fh1FFRatioZBackFolded = BookHistoArray();
704
705 //
706
707 fhnResponsePt = BookTHnSparseArray();
708 fhnResponseZ = BookTHnSparseArray();
709 fhnResponseXi = BookTHnSparseArray();
710
711 fh1FFTrackPtPrior = BookHistoArray();
712 fh1FFZPrior = BookHistoArray();
713 fh1FFXiPrior = BookHistoArray();
714
715
716 // histos bins
717
718 fNHistoBinsPt = new Int_t[fNJetPtSlices];
719 fNHistoBinsXi = new Int_t[fNJetPtSlices];
720 fNHistoBinsZ = new Int_t[fNJetPtSlices];
721
722 for(Int_t i=0; i<fNJetPtSlices; i++) fNHistoBinsPt[i] = 0;
723 for(Int_t i=0; i<fNJetPtSlices; i++) fNHistoBinsXi[i] = 0;
724 for(Int_t i=0; i<fNJetPtSlices; i++) fNHistoBinsZ[i] = 0;
725
726 fHistoBinsPt = new TArrayD*[fNJetPtSlices];
727 fHistoBinsXi = new TArrayD*[fNJetPtSlices];
728 fHistoBinsZ = new TArrayD*[fNJetPtSlices];
729
730 for(Int_t i=0; i<fNJetPtSlices; i++) fHistoBinsPt[i] = new TArrayD(0);
731 for(Int_t i=0; i<fNJetPtSlices; i++) fHistoBinsXi[i] = new TArrayD(0);
732 for(Int_t i=0; i<fNJetPtSlices; i++) fHistoBinsZ[i] = new TArrayD(0);
733}
734
735//_____________________________________________________________________________________________________________________________________
736void AliFragmentationFunctionCorrections::SetHistoBins(const Int_t jetPtSlice, const Int_t sizeBins, Double_t* bins, const Int_t type)
737{
738 // set histo bins for jet pt slice
739 // if binning undefined for any slice, original binning will be used
740
741 if(!fNJetPtSlices){
742 Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
743 return;
744 }
745
746 if(jetPtSlice>=fNJetPtSlices){
747 Printf("%s:%d -- jetPtSlice %d exceeds max",(char*)__FILE__,__LINE__,jetPtSlice);
748 return;
749 }
750
751 if(type == kFlagPt){
752 fNHistoBinsPt[jetPtSlice] = sizeBins-1;
753 fHistoBinsPt[jetPtSlice]->Set(sizeBins,bins);
754 }
755 else if(type==kFlagZ){
756 fNHistoBinsZ[jetPtSlice] = sizeBins-1;
757 fHistoBinsZ[jetPtSlice]->Set(sizeBins,bins);
758 }
759
760 else if(type==kFlagXi){
761 fNHistoBinsXi[jetPtSlice] = sizeBins-1;
762 fHistoBinsXi[jetPtSlice]->Set(sizeBins,bins);
763 }
764}
765
766//__________________________________________________________________________________________________________________________________________________________________
767void AliFragmentationFunctionCorrections::SetHistoBins(const Int_t jetPtSlice, const Int_t nBinsLimits, Double_t* binsLimits, Double_t* binsWidth, const Int_t type)
768{
769 // set histo bins for jet pt slice
770 // function takes array of limits and widths (e.g. 1 GeV bins up to 10 GeV, 2 GeV width up to 20 GeV, ...)
771 // array size of binsLimits: nBinsLimits
772 // array size of binsWidth: nBinsLimits-1
773 // binsLimits have to be in increasing order
774 // if binning undefined for any slice, original binning will be used
775
776 if(!fNJetPtSlices){
777 Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
778 return;
779 }
780
781 if(jetPtSlice>=fNJetPtSlices){
782 Printf("%s:%d -- jetPtSlice %d exceeds max",(char*)__FILE__,__LINE__,jetPtSlice);
783 return;
784 }
785
786
787 Double_t binLimitMin = binsLimits[0];
788 Double_t binLimitMax = binsLimits[nBinsLimits-1];
789
790 Double_t binLimit = binLimitMin; // start value
791
792 Int_t sizeUpperLim = 10000; //static_cast<Int_t>(binLimitMax/binsWidth[0])+1;
793 TArrayD binsArray(sizeUpperLim);
794 Int_t nBins = 0;
795 binsArray.SetAt(binLimitMin,nBins++);
796
797 while(binLimit<binLimitMax && nBins<sizeUpperLim){
798
799 Int_t currentSlice = -1;
800 for(Int_t i=0; i<nBinsLimits; i++){
801 if(binLimit >= 0.999*binsLimits[i]) currentSlice = i; // 0.999 numerical saftey factor
802 }
803
804 Double_t currentBinWidth = binsWidth[currentSlice];
805 binLimit += currentBinWidth;
806
807 binsArray.SetAt(binLimit,nBins++);
808 }
809
810 Double_t* bins = binsArray.GetArray();
811
812 SetHistoBins(jetPtSlice,nBins,bins,type);
813}
814
815//__________________________________________________________________________________________________
816void AliFragmentationFunctionCorrections::SetHistoBinsSinglePt(const Int_t sizeBins, Double_t* bins)
817{
818 // set histo bins for inclusive pt spectra
819 // if binning undefined, original binning will be used
820
821 fNHistoBinsSinglePt = sizeBins-1;
822
823 fHistoBinsSinglePt = new TArrayD(sizeBins);
824 fHistoBinsSinglePt->Set(sizeBins,bins);
825}
826
827//__________________________________________________________________________________________________________________________________________________________________
828void AliFragmentationFunctionCorrections::SetHistoBinsSinglePt(const Int_t nBinsLimits, Double_t* binsLimits, Double_t* binsWidth)
829{
830 // set histo bins for inclusive pt spectra
831 // function takes array of limits and widths (e.g. 1 GeV bins up to 10 GeV, 2 GeV width up to 20 GeV, ...)
832 // array size of binsLimits: nBinsLimits
833 // array size of binsWidth: nBinsLimits-1
834 // binsLimits have to be in increasing order
835 // if binning undefined for any slice, original binning will be used
836
837
838 Double_t binLimitMin = binsLimits[0];
839 Double_t binLimitMax = binsLimits[nBinsLimits-1];
840
841 Double_t binLimit = binLimitMin; // start value
842
843 Int_t sizeUpperLim = 10000; //static_cast<Int_t>(binLimitMax/binsWidth[0])+1;
844 TArrayD binsArray(sizeUpperLim);
845 Int_t nBins = 0;
846 binsArray.SetAt(binLimitMin,nBins++);
847
848 while(binLimit<binLimitMax && nBins<sizeUpperLim){
849
850 Int_t currentSlice = -1;
851 for(Int_t i=0; i<nBinsLimits; i++){
852 if(binLimit >= 0.999*binsLimits[i]) currentSlice = i; // 0.999 numerical saftey factor
853 }
854
855 Double_t currentBinWidth = binsWidth[currentSlice];
856 binLimit += currentBinWidth;
857
858 binsArray.SetAt(binLimit,nBins++);
859 }
860
861 Double_t* bins = binsArray.GetArray();
862
863 SetHistoBinsSinglePt(nBins,bins);
864}
865
866//____________________________________________________________________________________
867void AliFragmentationFunctionCorrections::NormalizeTH1(TH1* hist, const Float_t nJets)
868{
869 // FF normalization: divide by bin width and normalize to entries/jets
870 // should also work for TH2, but to be tested !!!
871
872 if(nJets == 0){ // nothing to do
873 if(fDebug>0) Printf("%s:%d -- normalize: nJets = 0, do nothing", (char*)__FILE__,__LINE__);
874 return;
875 }
876
877 Int_t nBins = hist->GetNbinsX()*hist->GetNbinsY()*hist->GetNbinsZ();
878
879 for(Int_t bin=0; bin<=nBins; bin++){ // include under-/overflow (?)
880
881 Double_t binWidth = hist->GetBinWidth(bin);
882 Double_t binCont = hist->GetBinContent(bin);
883 Double_t binErr = hist->GetBinError(bin);
884
885 binCont /= binWidth;
886 binErr /= binWidth;
887
888 hist->SetBinContent(bin,binCont);
889 hist->SetBinError(bin,binErr);
890 }
891
892 hist->Scale(1/nJets);
893}
894
895//_____________________________________________________
896void AliFragmentationFunctionCorrections::NormalizeFF()
897{
898 // normalize FF
899
900 if(fNCorrectionLevels>1){
901 Printf("%s:%d -- FF normalization should be done BEFORE any correction !!!", (char*)__FILE__,__LINE__);
902 }
903
904 for(Int_t i=0; i<fNJetPtSlices; i++){
905
906 if(fDebug>0) Printf(" normalizeFF: i %d, nJets %f",i,fNJets->At(i));
907
908 NormalizeTH1(fCorrFF[0]->GetTrackPt(i),fNJets->At(i)); // always normalize corr level 0 = raw FF
909 NormalizeTH1(fCorrFF[0]->GetZ(i),fNJets->At(i));
910 NormalizeTH1(fCorrFF[0]->GetXi(i),fNJets->At(i));
911 }
912}
913
914//______________________________________________________
915void AliFragmentationFunctionCorrections::NormalizeBgr()
916{
917 // normalize bgr/UE FF
918
919 if(fNCorrectionLevelsBgr>1){
920 Printf("%s:%d -- FF normalization should be done BEFORE any correction !!!", (char*)__FILE__,__LINE__);
921 }
922
923 for(Int_t i=0; i<fNJetPtSlices; i++){
924 NormalizeTH1(fCorrBgr[0]->GetTrackPt(i), fNJetsBgr->At(i)); // always normalize corr level 0 = raw FF
925 NormalizeTH1(fCorrBgr[0]->GetZ(i), fNJetsBgr->At(i));
926 NormalizeTH1(fCorrBgr[0]->GetXi(i),fNJetsBgr->At(i));
927 }
928
929}
930
931//__________________________________________________________________________________________________
932void AliFragmentationFunctionCorrections::ReadRawFF(TString strfile, TString strID, TString strFFID)
933{
934 // read raw FF - standard dir/list name
935
936 TString strdir = "PWG4_FragmentationFunction_" + strID;
937 TString strlist = "fracfunc_" + strID;
938
939 ReadRawFF(strfile,strdir,strlist,strFFID);
940}
941
942//____________________________________________________________________________________________________________________
943void AliFragmentationFunctionCorrections::ReadRawFF(TString strfile, TString strdir, TString strlist, TString strFFID)
944{
945 // get raw FF from input file, project in jet pt slice
946 // normalization done separately
947
948 TFile f(strfile,"READ");
949
950 if(!f.IsOpen()){
951 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
952 return;
953 }
954
955 if(fDebug>0) Printf("%s:%d -- read FF from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
956
957 gDirectory->cd(strdir);
958
959 TList* list = 0;
960
961 if(!(list = (TList*) gDirectory->Get(strlist))){
962 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
963 return;
964 }
965
966 TString hnameJetPt(Form("fh1FFJetPt%s",strFFID.Data()));
967 TString hnameTrackPt(Form("fh2FFTrackPt%s",strFFID.Data()));
968 TString hnameZ(Form("fh2FFZ%s",strFFID.Data()));
969 TString hnameXi(Form("fh2FFXi%s",strFFID.Data()));
970
971 TH1F* fh1FFJetPt = (TH1F*) list->FindObject(hnameJetPt);
972 TH2F* fh2FFTrackPt = (TH2F*) list->FindObject(hnameTrackPt);
973 TH2F* fh2FFZ = (TH2F*) list->FindObject(hnameZ);
974 TH2F* fh2FFXi = (TH2F*) list->FindObject(hnameXi);
975
976 if(!fh1FFJetPt) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameJetPt.Data()); return; }
977 if(!fh2FFTrackPt){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameTrackPt.Data()); return; }
978 if(!fh2FFZ) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameZ.Data()); return; }
979 if(!fh2FFXi) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameXi.Data()); return; }
980
981 fh1FFJetPt->SetDirectory(0);
982 fh2FFTrackPt->SetDirectory(0);
983 fh2FFZ->SetDirectory(0);
984 fh2FFXi->SetDirectory(0);
985
986 f.Close();
987
988
989 // nJets per bin
990
991 for(Int_t i=0; i<fNJetPtSlices; i++){
992
993 Float_t jetPtLoLim = fJetPtSlices->At(i);
994 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
995
996 Int_t binLo = static_cast<Int_t>(fh1FFJetPt->FindBin(jetPtLoLim));
997 Int_t binUp = static_cast<Int_t>(fh1FFJetPt->FindBin(jetPtUpLim)) - 1;
998
999 Float_t nJetsBin = fh1FFJetPt->Integral(binLo,binUp);
1000
1001 fNJets->SetAt(nJetsBin,i);
1002
1003 if(fDebug>0) Printf("jet pt %d to %d: nJets %f",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),fNJets->At(i));
1004 }
1005
1006 // projections: FF
1007
1008 for(Int_t i=0; i<fNJetPtSlices; i++){
1009
1010 Float_t jetPtLoLim = fJetPtSlices->At(i);
1011 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1012
1013 Int_t binLo = static_cast<Int_t>(fh2FFTrackPt->GetXaxis()->FindBin(jetPtLoLim));
1014 Int_t binUp = static_cast<Int_t>(fh2FFTrackPt->GetXaxis()->FindBin(jetPtUpLim))-1;
1015
1016 if(binUp > fh2FFTrackPt->GetNbinsX()){
1017 Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
1018 return;
1019 }
1020
1021 TString strNameFFPt(Form("fh1FFTrackPt%s_%02d_%02d",strFFID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1022 TString strNameFFZ(Form("fh1FFZ%s_%02d_%02d",strFFID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1023 TString strNameFFXi(Form("fh1FFXi%s_%02d_%02d",strFFID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1024
1025 // appendix 'unbinned' to avoid histos with same name after rebinning
1026 TH1F* projPt = (TH1F*) fh2FFTrackPt->ProjectionY(strNameFFPt+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
1027 TH1F* projZ = (TH1F*) fh2FFZ->ProjectionY(strNameFFZ+"_unBinned",binLo,binUp,"o");
1028 TH1F* projXi = (TH1F*) fh2FFXi->ProjectionY(strNameFFXi+"_unBinned",binLo,binUp,"o");
1029
1030 if(fNHistoBinsPt[i]) projPt = (TH1F*) projPt->Rebin(fNHistoBinsPt[i],strNameFFPt,fHistoBinsPt[i]->GetArray());
1031 if(fNHistoBinsZ[i]) projZ = (TH1F*) projZ->Rebin(fNHistoBinsZ[i],strNameFFZ,fHistoBinsZ[i]->GetArray());
1032 if(fNHistoBinsXi[i]) projXi = (TH1F*) projXi->Rebin(fNHistoBinsXi[i],strNameFFXi,fHistoBinsXi[i]->GetArray());
1033
1034 projPt->SetNameTitle(strNameFFPt,"");
1035 projZ->SetNameTitle(strNameFFZ,"");
1036 projXi->SetNameTitle(strNameFFXi,"");
1037
1038 // raw FF = corr level 0
1039 fCorrFF[0]->AddCorrHistos(i,projPt,projZ,projXi);
1040 }
1041}
1042
1043//_____________________________________________________________________________________________________________________
1044void AliFragmentationFunctionCorrections::ReadRawBgr(TString strfile, TString strID, TString strBgrID, TString strFFID)
1045{
1046 // read raw FF - standard dir/list name
1047
1048 TString strdir = "PWG4_FragmentationFunction_" + strID;
1049 TString strlist = "fracfunc_" + strID;
1050
1051 ReadRawBgr(strfile,strdir,strlist,strBgrID,strFFID);
1052}
1053
1054//_______________________________________________________________________________________________________________________________________
1055void AliFragmentationFunctionCorrections::ReadRawBgr(TString strfile, TString strdir, TString strlist, TString strBgrID, TString strFFID)
1056{
1057 // get raw FF from input file, project in jet pt slice
1058 // use jet dN/dpt corresponding to strFFID, bgr FF to strBgrID+strFFID
1059 // e.g. "fh1FFJetPtRecCuts", "fh2FFXiBgrPerpRecCuts"
1060 // normalization done separately
1061
1062 TString strID = strBgrID + strFFID;
1063
1064 TFile f(strfile,"READ");
1065
1066 if(!f.IsOpen()){
1067 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1068 return;
1069 }
1070
1071 if(fDebug>0) Printf("%s:%d -- read Bgr %s from file %s ",(char*)__FILE__,__LINE__,strBgrID.Data(),strfile.Data());
1072
1073 gDirectory->cd(strdir);
1074
1075 TList* list = 0;
1076
1077 if(!(list = (TList*) gDirectory->Get(strlist))){
1078 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1079 return;
1080 }
1081
1082 TString hnameNJets = "fh1nRecJetsCuts";
1083 TString hnameJetPt(Form("fh1FFJetPt%s",strFFID.Data())); // not: strID.Data() !!! would not be proper normalization
1084 TString hnameBgrTrackPt(Form("fh2FFTrackPt%s",strID.Data()));
1085 TString hnameBgrZ(Form("fh2FFZ%s",strID.Data()));
1086 TString hnameBgrXi(Form("fh2FFXi%s",strID.Data()));
1087
1088 TH1F* fh1NJets = (TH1F*) list->FindObject(hnameNJets); // needed for normalization of bgr out of 2 jets
1089 TH1F* fh1FFJetPtBgr = (TH1F*) list->FindObject(hnameJetPt);
1090 TH2F* fh2FFTrackPtBgr = (TH2F*) list->FindObject(hnameBgrTrackPt);
1091 TH2F* fh2FFZBgr = (TH2F*) list->FindObject(hnameBgrZ);
1092 TH2F* fh2FFXiBgr = (TH2F*) list->FindObject(hnameBgrXi);
1093
1094 if(!fh1FFJetPtBgr) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameJetPt.Data()); return; }
1095 if(!fh1NJets) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameNJets.Data()); return; }
1096 if(!fh2FFTrackPtBgr){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrTrackPt.Data()); return; }
1097 if(!fh2FFZBgr) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrZ.Data()); return; }
1098 if(!fh2FFXiBgr) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrXi.Data()); return; }
1099
1100 fh1FFJetPtBgr->SetDirectory(0);
1101 fh1NJets->SetDirectory(0);
1102 fh2FFTrackPtBgr->SetDirectory(0);
1103 fh2FFZBgr->SetDirectory(0);
1104 fh2FFXiBgr->SetDirectory(0);
1105
1106 f.Close();
1107
1108 // nJets per bin
1109
1110 for(Int_t i=0; i<fNJetPtSlices; i++){
1111
1112 Float_t jetPtLoLim = fJetPtSlices->At(i);
1113 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1114
1115 Int_t binLo = static_cast<Int_t>(fh1FFJetPtBgr->FindBin(jetPtLoLim));
1116 Int_t binUp = static_cast<Int_t>(fh1FFJetPtBgr->FindBin(jetPtUpLim)) - 1;
1117
1118 Float_t nJetsBin = fh1FFJetPtBgr->Integral(binLo,binUp);
1119 Double_t scaleF = 1;
1120
1121 //if(strBgrID.Contains("Out2Jets")){ // scale by ratio 2 jets events / all events
1122 // scaleF = fh1NJets->Integral(fh1NJets->FindBin(2),fh1NJets->GetNbinsX())
1123 // / fh1NJets->Integral(fh1NJets->FindBin(1),fh1NJets->GetNbinsX());}
1124
1125
1126 if(strBgrID.Contains("OutAllJets")){ // scale by ratio >3 jets events / all events
1127 scaleF = fh1NJets->Integral(fh1NJets->FindBin(4),fh1NJets->GetNbinsX())
1128 / fh1NJets->Integral(fh1NJets->FindBin(1),fh1NJets->GetNbinsX());
1129 }
1130
1131 fNJetsBgr->SetAt(nJetsBin*scaleF,i);
1132
1133 if(fDebug>0) Printf("bgr jet pt %d to %d: nJets %f, scaleF %.2f",
1134 static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),nJetsBin,scaleF);
1135
1136 }
1137
1138 // projections: FF
1139
1140 for(Int_t i=0; i<fNJetPtSlices; i++){
1141
1142 Float_t jetPtLoLim = fJetPtSlices->At(i);
1143 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1144
1145 Int_t binLo = static_cast<Int_t>(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtLoLim));
1146 Int_t binUp = static_cast<Int_t>(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtUpLim))-1;
1147
1148 if(binUp > fh2FFTrackPtBgr->GetNbinsX()){
1149 Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
1150 return;
1151 }
1152
1153 TString strNameBgrPt(Form("fh1BgrTrackPt%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1154 TString strNameBgrZ(Form("fh1BgrZ%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1155 TString strNameBgrXi(Form("fh1BgrXi%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1156
1157 // appendix 'unbinned' to avoid histos with same name after rebinning
1158 TH1F* projPt = (TH1F*) fh2FFTrackPtBgr->ProjectionY(strNameBgrPt+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
1159 TH1F* projZ = (TH1F*) fh2FFZBgr->ProjectionY(strNameBgrZ+"_unBinned",binLo,binUp,"o");
1160 TH1F* projXi = (TH1F*) fh2FFXiBgr->ProjectionY(strNameBgrXi+"_unBinned",binLo,binUp,"o");
1161
1162 if(fNHistoBinsPt[i]) projPt = (TH1F*) projPt->Rebin(fNHistoBinsPt[i],strNameBgrPt,fHistoBinsPt[i]->GetArray());
1163 if(fNHistoBinsZ[i]) projZ = (TH1F*) projZ->Rebin(fNHistoBinsZ[i],strNameBgrZ,fHistoBinsZ[i]->GetArray());
1164 if(fNHistoBinsXi[i]) projXi = (TH1F*) projXi->Rebin(fNHistoBinsXi[i],strNameBgrXi,fHistoBinsXi[i]->GetArray());
1165
1166 projPt->SetNameTitle(strNameBgrPt,"");
1167 projZ->SetNameTitle(strNameBgrZ,"");
1168 projXi->SetNameTitle(strNameBgrXi,"");
1169
1170 // raw bgr = corr level 0
1171 fCorrBgr[0]->AddCorrHistos(i,projPt,projZ,projXi);
1172 }
1173}
1174
1175//_____________________________________________________________________________________________________________________
1176void AliFragmentationFunctionCorrections::ReadRawBgrEmbedding(TString strfile, TString strID, TString strFFID)
1177{
1178 // read raw FF - standard dir/list name
1179
1180 TString strdir = "PWG4_FragmentationFunction_" + strID;
1181 TString strlist = "fracfunc_" + strID;
1182
1183 ReadRawBgrEmbedding(strfile,strdir,strlist,strFFID);
1184}
1185
1186//_______________________________________________________________________________________________________________________________________
1187void AliFragmentationFunctionCorrections::ReadRawBgrEmbedding(TString strfile, TString strdir, TString strlist, TString strFFID)
1188{
1189 // get raw FF from input file, project in jet pt slice
1190 // for embedding, the bgr FF are taken from histos "fh1FFJetPtRecCuts", "fh2FFXiRecCuts"
1191 // normalization done separately
1192
1193 TString strBgrID = "BckgEmbed";
1194 TString strID = strBgrID + strFFID;
1195
1196 TFile f(strfile,"READ");
1197
1198 if(!f.IsOpen()){
1199 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1200 return;
1201 }
1202
1203 if(fDebug>0) Printf("%s:%d -- read Bgr %s from file %s ",(char*)__FILE__,__LINE__,strFFID.Data(),strfile.Data());
1204
1205 gDirectory->cd(strdir);
1206
1207 TList* list = 0;
1208
1209 if(!(list = (TList*) gDirectory->Get(strlist))){
1210 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1211 return;
1212 }
1213
1214 TString hnameNJets = "fh1nRecJetsCuts";
1215 TString hnameJetPt(Form("fh1FFJetPt%s",strFFID.Data()));
1216 TString hnameBgrTrackPt(Form("fh2FFTrackPt%s",strFFID.Data()));
1217 TString hnameBgrZ(Form("fh2FFZ%s",strFFID.Data()));
1218 TString hnameBgrXi(Form("fh2FFXi%s",strFFID.Data()));
1219
1220 TH1F* fh1NJets = (TH1F*) list->FindObject(hnameNJets); // needed for normalization of bgr out of 2 jets
1221 TH1F* fh1FFJetPtBgr = (TH1F*) list->FindObject(hnameJetPt);
1222 TH2F* fh2FFTrackPtBgr = (TH2F*) list->FindObject(hnameBgrTrackPt);
1223 TH2F* fh2FFZBgr = (TH2F*) list->FindObject(hnameBgrZ);
1224 TH2F* fh2FFXiBgr = (TH2F*) list->FindObject(hnameBgrXi);
1225
1226 if(!fh1FFJetPtBgr) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameJetPt.Data()); return; }
1227 if(!fh1NJets) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameNJets.Data()); return; }
1228 if(!fh2FFTrackPtBgr){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrTrackPt.Data()); return; }
1229 if(!fh2FFZBgr) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrZ.Data()); return; }
1230 if(!fh2FFXiBgr) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrXi.Data()); return; }
1231
1232 fh1FFJetPtBgr->SetDirectory(0);
1233 fh1NJets->SetDirectory(0);
1234 fh2FFTrackPtBgr->SetDirectory(0);
1235 fh2FFZBgr->SetDirectory(0);
1236 fh2FFXiBgr->SetDirectory(0);
1237
1238 f.Close();
1239
1240 // nJets per bin
1241
1242 for(Int_t i=0; i<fNJetPtSlices; i++){
1243
1244 Float_t jetPtLoLim = fJetPtSlices->At(i);
1245 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1246
1247 Int_t binLo = static_cast<Int_t>(fh1FFJetPtBgr->FindBin(jetPtLoLim));
1248 Int_t binUp = static_cast<Int_t>(fh1FFJetPtBgr->FindBin(jetPtUpLim)) - 1;
1249
1250 Float_t nJetsBin = fh1FFJetPtBgr->Integral(binLo,binUp);
1251 Double_t scaleF = 1;
1252
1253 fNJetsBgr->SetAt(nJetsBin*scaleF,i);
1254
1255 if(fDebug>0) Printf("bgr jet pt %d to %d: nJets %f, scaleF %.2f",
1256 static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),nJetsBin,scaleF);
1257
1258 }
1259
1260 // projections: FF
1261
1262 for(Int_t i=0; i<fNJetPtSlices; i++){
1263
1264 Float_t jetPtLoLim = fJetPtSlices->At(i);
1265 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1266
1267 Int_t binLo = static_cast<Int_t>(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtLoLim));
1268 Int_t binUp = static_cast<Int_t>(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtUpLim))-1;
1269
1270 if(binUp > fh2FFTrackPtBgr->GetNbinsX()){
1271 Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
1272 return;
1273 }
1274
1275 TString strNameBgrPt(Form("fh1BgrTrackPt%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1276 TString strNameBgrZ(Form("fh1BgrZ%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1277 TString strNameBgrXi(Form("fh1BgrXi%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1278
1279 // appendix 'unbinned' to avoid histos with same name after rebinning
1280 TH1F* projPt = (TH1F*) fh2FFTrackPtBgr->ProjectionY(strNameBgrPt+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
1281 TH1F* projZ = (TH1F*) fh2FFZBgr->ProjectionY(strNameBgrZ+"_unBinned",binLo,binUp,"o");
1282 TH1F* projXi = (TH1F*) fh2FFXiBgr->ProjectionY(strNameBgrXi+"_unBinned",binLo,binUp,"o");
1283
1284 if(fNHistoBinsPt[i]) projPt = (TH1F*) projPt->Rebin(fNHistoBinsPt[i],strNameBgrPt,fHistoBinsPt[i]->GetArray());
1285 if(fNHistoBinsZ[i]) projZ = (TH1F*) projZ->Rebin(fNHistoBinsZ[i],strNameBgrZ,fHistoBinsZ[i]->GetArray());
1286 if(fNHistoBinsXi[i]) projXi = (TH1F*) projXi->Rebin(fNHistoBinsXi[i],strNameBgrXi,fHistoBinsXi[i]->GetArray());
1287
1288 projPt->SetNameTitle(strNameBgrPt,"");
1289 projZ->SetNameTitle(strNameBgrZ,"");
1290 projXi->SetNameTitle(strNameBgrXi,"");
1291
1292 // raw bgr = corr level 0
1293 fCorrBgr[0]->AddCorrHistos(i,projPt,projZ,projXi);
1294 }
1295}
1296
1297
1298//__________________________________________________________________________________________________________
1299void AliFragmentationFunctionCorrections::WriteOutput(TString strfile, TString strdir, Bool_t updateOutfile)
1300{
1301 // write histos to file
1302 // skip histos with 0 entries
1303
1304 TString outfileOption = "RECREATE";
1305 if(updateOutfile) outfileOption = "UPDATE";
1306
1307 TFile f(strfile,outfileOption);
1308
1309 if(!f.IsOpen()){
1310 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1311 return;
1312 }
1313
1314 if(fDebug>0) Printf("%s:%d -- write FF to file %s ",(char*)__FILE__,__LINE__,strfile.Data());
1315
1316 if(strdir && strdir.Length()){
1317 TDirectory* dir = f.mkdir(strdir);
1318 dir->cd();
1319 }
1320
1321 for(Int_t i=0; i<fNJetPtSlices; i++){
1322
1323 for(Int_t c=0; c<fNCorrectionLevels; c++) if(fCorrFF[c]->GetTrackPt(i)->GetEntries()) fCorrFF[c]->GetTrackPt(i)->Write();
1324 for(Int_t c=0; c<fNCorrectionLevels; c++) if(fCorrFF[c]->GetZ(i)->GetEntries()) fCorrFF[c]->GetZ(i)->Write();
1325 for(Int_t c=0; c<fNCorrectionLevels; c++) if(fCorrFF[c]->GetXi(i)->GetEntries()) fCorrFF[c]->GetXi(i)->Write();
1326
1327 if(fh1FFXiShift[i]->GetEntries()) fh1FFXiShift[i]->Write();
1328
1329 for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) if(fCorrBgr[c]->GetTrackPt(i)->GetEntries()) fCorrBgr[c]->GetTrackPt(i)->Write();
1330 for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) if(fCorrBgr[c]->GetZ(i)->GetEntries()) fCorrBgr[c]->GetZ(i)->Write();
1331 for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) if(fCorrBgr[c]->GetXi(i)->GetEntries()) fCorrBgr[c]->GetXi(i)->Write();
1332
1333
1334 if(fh1FFTrackPtBackFolded[i] && fh1FFTrackPtBackFolded[i]->GetEntries()) fh1FFTrackPtBackFolded[i]->Write();
1335 if(fh1FFZBackFolded[i] && fh1FFZBackFolded[i]->GetEntries()) fh1FFZBackFolded[i]->Write();
1336 if(fh1FFXiBackFolded[i] && fh1FFXiBackFolded[i]->GetEntries()) fh1FFXiBackFolded[i]->Write();
1337
1338
1339 if(fh1FFRatioTrackPtFolded[i] && fh1FFRatioTrackPtFolded[i]->GetEntries()) fh1FFRatioTrackPtFolded[i]->Write();
1340 if(fh1FFRatioZFolded[i] && fh1FFRatioZFolded[i]->GetEntries()) fh1FFRatioZFolded[i]->Write();
1341 if(fh1FFRatioXiFolded[i] && fh1FFRatioXiFolded[i]->GetEntries()) fh1FFRatioXiFolded[i]->Write();
1342
1343 if(fh1FFRatioTrackPtBackFolded[i] && fh1FFRatioTrackPtBackFolded[i]->GetEntries()) fh1FFRatioTrackPtBackFolded[i]->Write();
1344 if(fh1FFRatioZBackFolded[i] && fh1FFRatioZBackFolded[i]->GetEntries()) fh1FFRatioZBackFolded[i]->Write();
1345 if(fh1FFRatioXiBackFolded[i] &&fh1FFRatioXiBackFolded[i]->GetEntries()) fh1FFRatioXiBackFolded[i]->Write();
1346
1347 }
1348
1349 // inclusive track pt
1350
1351 for(Int_t c=0; c<fNCorrectionLevelsSinglePt; c++) if(fCorrSinglePt[c]->GetTrackPt(0)->GetEntries()) fCorrSinglePt[c]->GetTrackPt(0)->Write();
1352 if(fh1SingleTrackPtBackFolded) fh1SingleTrackPtBackFolded->Write();
1353 if(fh1RatioSingleTrackPtFolded) fh1RatioSingleTrackPtFolded->Write();
1354 if(fh1RatioSingleTrackPtBackFolded) fh1RatioSingleTrackPtBackFolded->Write();
1355
1356 f.Close();
1357}
1358
1359//____________________________________________________________________________________________________________________________________
1360THnSparse* AliFragmentationFunctionCorrections::TH1toSparse(const TH1F* hist, TString strName, TString strTit, const Bool_t fillConst)
1361{
1362 // copy 1-dimensional histo to THnSparse
1363 // if fillConst TRUE, create THnSparse with same binning as hist but all entries = 1
1364 // histos with variable bin size are supported
1365
1366 // note: function returns pointer to 'new' THnSparse on heap, object needs to be deleted by user
1367
1368 THnSparseF* fhnHist;
1369
1370 Int_t nBins = hist->GetXaxis()->GetNbins();
1371 Int_t nBinsVec[1] = { nBins };
1372
1373 const Double_t* binsVec = hist->GetXaxis()->GetXbins()->GetArray();
1374
1375 if(binsVec){ // binsVec only neq NULL if histo was rebinned before
1376
1377 fhnHist = new THnSparseF(strName,strTit,1,nBinsVec/*,binMinVec,binMaxVec*/);
1378 fhnHist->SetBinEdges(0,binsVec);
1379 }
1380 else{ // uniform bin size
1381
1382 Double_t xMin = hist->GetXaxis()->GetXmin();
1383 Double_t xMax = hist->GetXaxis()->GetXmax();
1384
1385 Double_t binMinVec[1] = { xMin };
1386 Double_t binMaxVec[1] = { xMax };
1387
1388 fhnHist = new THnSparseF(strName,strTit,1,nBinsVec,binMinVec,binMaxVec);
1389 }
1390
1391
1392 for(Int_t bin=1; bin<=nBins; bin++){
1393
1394 Double_t binCenter = fhnHist->GetAxis(0)->GetBinCenter(bin);
1395
1396 Double_t binCoord[] = {binCenter};
1397 fhnHist->Fill(binCoord,1); // initially need to create the bin
1398
1399 Long64_t binIndex = fhnHist->GetBin(binCoord);
1400
1401 Double_t cont = hist->GetBinContent(bin);
1402 Double_t err = hist->GetBinError(bin);
1403
1404 if(fillConst){
1405 cont = 1;
1406 err = 0;
1407 }
1408
1409 fhnHist->SetBinContent(binIndex,cont);
1410 fhnHist->SetBinError(binIndex,err);
1411 }
1412
1413 return fhnHist;
1414}
1415
1416//______________________________________________________________________________________________________________________________________________
1417THnSparse* AliFragmentationFunctionCorrections::Unfold(THnSparse* hnHist, const THnSparse* hnResponse, const THnSparse* hnEff, const Int_t nIter,
1418 const Bool_t useCorrelatedErrors, const THnSparse* hnPrior)
1419{
1420 // unfold input histo
1421
1422 AliCFUnfolding unfolding("unfolding","",1,hnResponse,hnEff,hnHist,hnPrior); // arg3: nVar; hnEff required, hnPrior not (defaults to 0x0)
1423 unfolding.SetMaxNumberOfIterations(nIter);
1424 // unfolding.SetMaxChi2PerDOF(1.e-07); // OBSOLETE !!!
1425 // if(useSmoothing) unfolding.UseSmoothing();
1426
1427 if(useCorrelatedErrors) unfolding.SetUseCorrelatedErrors();
1428 unfolding.Unfold();
1429
1430 THnSparse* unfolded = unfolding.GetUnfolded();
1431
1432 TString hnameUnf = hnHist->GetName();
1433 hnameUnf.Append("_unf");
1434 unfolded->SetNameTitle(hnameUnf,hnHist->GetTitle());
1435
1436 return unfolded;
1437}
1438
1439//___________________________________________________________________________________________________________________________
1440void AliFragmentationFunctionCorrections::UnfoldHistos(const Int_t nIter, const Bool_t useCorrelatedErrors, const Int_t type)
1441{
1442 // loop over jet pt slices and unfold dN/dpt spectra
1443
1444 TString labelCorr = fCorrFF[fNCorrectionLevels-1]->GetLabel();
1445 if(!labelCorr.Contains("unfold")) AddCorrectionLevel("unfold");
1446
1447 for(Int_t i=0; i<fNJetPtSlices; i++){
1448
1449 TH1F* hist = 0;
1aa4f09f 1450 if(type == kFlagPt) hist = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i); // level -2: before unfolding, level -1: unfolded
1451 else if(type == kFlagZ) hist = fCorrFF[fNCorrectionLevels-2]->GetZ(i); // level -2: before unfolding, level -1: unfolded
1452 else if(type == kFlagXi) hist = fCorrFF[fNCorrectionLevels-2]->GetXi(i); // level -2: before unfolding, level -1: unfolded
6daac008 1453 else{
1454 Printf("%s%d unknown type",(char*)__FILE__,__LINE__);
1455 return;
1456 }
39e2b057 1457
1458 THnSparse* hnResponse = 0;
6daac008 1459 if(type == kFlagPt) hnResponse = fhnResponsePt[i];
1aa4f09f 1460 else if(type == kFlagZ) hnResponse = fhnResponseZ[i];
1461 else if(type == kFlagXi) hnResponse = fhnResponseXi[i];
6daac008 1462 else{
1463 Printf("%s%d unknown type",(char*)__FILE__,__LINE__);
1464 return;
1465 }
39e2b057 1466
1467 TH1F* hPrior = 0;
1468 if(type == kFlagPt && fh1FFTrackPtPrior[i] && ((TString(fh1FFTrackPtPrior[i]->GetName())).Length() > 0) ) hPrior = fh1FFTrackPtPrior[i];
a5592cfa 1469 else if(type == kFlagZ && fh1FFZPrior[i] && ((TString(fh1FFZPrior[i]->GetName())).Length() > 0) ) hPrior = fh1FFZPrior[i];
1470 else if(type == kFlagXi && fh1FFXiPrior[i] && ((TString(fh1FFXiPrior[i]->GetName())).Length() > 0) ) hPrior = fh1FFXiPrior[i];
1aa4f09f 1471
39e2b057 1472
1473 TString histNameTHn = hist->GetName();
1474 histNameTHn.ReplaceAll("TH1","THn");
1475
1476 TString priorNameTHn;
1477 if(hPrior){
1478 priorNameTHn = hPrior->GetName();
1479 priorNameTHn.ReplaceAll("TH1","THn");
1480 }
1481
1482 TString histNameBackFolded = hist->GetName();
1483 histNameBackFolded.Append("_backfold");
1484
1485 TString histNameRatioFolded = hist->GetName();
1486 histNameRatioFolded.ReplaceAll("fh1FF","hRatioFF");
1487 histNameRatioFolded.Append("_unfold");
1488
1489 TString histNameRatioBackFolded = hist->GetName();
1490 histNameRatioBackFolded.ReplaceAll("fh1FF","hRatioFF");
1491 histNameRatioBackFolded.Append("_backfold");
1492
1493 THnSparse* hnHist = TH1toSparse(hist,histNameTHn,hist->GetTitle());
1494 THnSparse* hnFlatEfficiency = TH1toSparse(hist,"fhnEfficiency","eff",kTRUE); // could optionally also use real eff
1495 THnSparse* hnPrior = 0;
1496 if(hPrior) hnPrior = TH1toSparse(hPrior,priorNameTHn,hPrior->GetTitle());
1497
1498 THnSparse* hnUnfolded
1499 = Unfold(hnHist,hnResponse,hnFlatEfficiency,nIter,useCorrelatedErrors,hnPrior);
1500
6daac008 1501 TH1F* hUnfolded = (TH1F*) hnUnfolded->Projection(0);
1502 if(hist)
47f544f7 1503 hUnfolded->SetNameTitle(hist->GetName(),hist->GetTitle());
39e2b057 1504
1505 if(type == kFlagPt) fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hUnfolded,0,0);
1506 if(type == kFlagZ) fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,0,hUnfolded,0);
1507 if(type == kFlagXi) fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,0,0,hUnfolded);
1508
1509 // backfolding: apply response matrix to unfolded spectrum
6daac008 1510 TH1F* hBackFolded = ApplyResponse(hUnfolded,hnResponse);
1511 if(hist)
47f544f7 1512 hBackFolded->SetNameTitle(histNameBackFolded,hist->GetTitle());
39e2b057 1513
1514 if(type == kFlagPt) fh1FFTrackPtBackFolded[i] = hBackFolded;
1515 if(type == kFlagZ) fh1FFZBackFolded[i] = hBackFolded;
1516 if(type == kFlagXi) fh1FFXiBackFolded[i] = hBackFolded;
1517
1518 // ratio unfolded to original histo
1519 TH1F* hRatioUnfolded = (TH1F*) hUnfolded->Clone(histNameRatioFolded);
1520 hRatioUnfolded->Reset();
6daac008 1521 if (hist)
47f544f7 1522 hRatioUnfolded->Divide(hUnfolded,hist,1,1,"B");
39e2b057 1523
1524 if(type == kFlagPt) fh1FFRatioTrackPtFolded[i] = hRatioUnfolded;
1525 if(type == kFlagZ) fh1FFRatioZFolded[i] = hRatioUnfolded;
1526 if(type == kFlagXi) fh1FFRatioXiFolded[i] = hRatioUnfolded;
1527
1528
1529 // ratio backfolded to original histo
1530 TH1F* hRatioBackFolded = (TH1F*) hBackFolded->Clone(histNameRatioBackFolded);
1531 hRatioBackFolded->Reset();
1532 hRatioBackFolded->Divide(hBackFolded,hist,1,1,"B");
1533
1534 if(type == kFlagPt) fh1FFRatioTrackPtBackFolded[i] = hRatioBackFolded;
1535 if(type == kFlagZ) fh1FFRatioZBackFolded[i] = hRatioBackFolded;
1536 if(type == kFlagXi) fh1FFRatioXiBackFolded[i] = hRatioBackFolded;
1537
1538 delete hnHist;
1539 delete hnFlatEfficiency;
1540
1541 }
1542}
1543
1544//_____________________________________________________________________________________________________
1545void AliFragmentationFunctionCorrections::UnfoldPt(const Int_t nIter, const Bool_t useCorrelatedErrors)
1546{
1547
1548 Int_t type = kFlagPt;
1549 UnfoldHistos(nIter, useCorrelatedErrors, type);
1550}
1551
1552//_____________________________________________________________________________________________________
1553void AliFragmentationFunctionCorrections::UnfoldZ(const Int_t nIter, const Bool_t useCorrelatedErrors)
1554{
1555
1556 Int_t type = kFlagZ;
1557 UnfoldHistos(nIter, useCorrelatedErrors, type);
1558}
1559
1560//_____________________________________________________________________________________________________
1561void AliFragmentationFunctionCorrections::UnfoldXi(const Int_t nIter, const Bool_t useCorrelatedErrors)
1562{
1563
1564 Int_t type = kFlagXi;
1565 UnfoldHistos(nIter, useCorrelatedErrors, type);
1566}
1567
1568//______________________________________________________________________________________________
1569TH1F* AliFragmentationFunctionCorrections::ApplyResponse(const TH1F* hist, THnSparse* hnResponse)
1570{
1571 // apply (multiply) response matrix to hist
1572
1573 TH1F* hOut = new TH1F(*hist);
1574 hOut->Reset();
1575
1576 const Int_t axisM = 0;
1577 const Int_t axisT = 1;
1578
1579 Int_t nBinsM = hnResponse->GetAxis(axisM)->GetNbins();
1580 Int_t nBinsT = hnResponse->GetAxis(axisT)->GetNbins();
1581
1582 // response matrix normalization
1583 // do it in this function and not when reading response, since for 'proper' normalization errors are difficult to assign
1584 // so for unfolding proper we leave it to CORRFW ...
1585
1586 Double_t normFacResponse[nBinsT];
1587
1588 for(Int_t iT=1; iT<=nBinsT; iT++){
1589
1590 Double_t sumResp = 0;
1591
1592 for(Int_t iM=1; iM<=nBinsM; iM++){
1593
1594 Double_t coordM = hnResponse->GetAxis(axisM)->GetBinCenter(iM);
1595 Double_t coordT = hnResponse->GetAxis(axisT)->GetBinCenter(iT);
1596
1597 Double_t binCoord[] = {coordM,coordT};
1598
1599 Long64_t binIndex = hnResponse->GetBin(binCoord);
1600
1601 Double_t resp = hnResponse->GetBinContent(binIndex);
1602
1603 sumResp += resp;
1604 }
1605
1606 normFacResponse[iT] = 0;
1607 if(sumResp) normFacResponse[iT] = 1/sumResp;
1608 }
1609
1610
1611
1612 for(Int_t iM=1; iM<=nBinsM; iM++){
1613
1614 Double_t contM = 0;
1615
1616 for(Int_t iT=1; iT<=nBinsT; iT++){
1617
1618 Double_t contT = hist->GetBinContent(iT);
1619
1620 Double_t coordM = hnResponse->GetAxis(axisM)->GetBinCenter(iM);
1621 Double_t coordT = hnResponse->GetAxis(axisT)->GetBinCenter(iT);
1622
1623 Double_t binCoord[] = {coordM,coordT};
1624
1625 Long64_t binIndex = hnResponse->GetBin(binCoord);
1626
1627 Double_t resp = hnResponse->GetBinContent(binIndex);
1628
1629 contM += resp*normFacResponse[iT]*contT;
1630 }
1631
1632 hOut->SetBinContent(iM,contM);
1633 }
1634
1635 return hOut;
1636}
1637
1638//_______________________________________________________________________________________________________
1639void AliFragmentationFunctionCorrections::ReadEfficiency(TString strfile, TString strdir, TString strlist)
1640{
1641 // read reconstruction efficiency from file
1642 // argument strlist optional - read from directory strdir if not specified
1643
1644 // temporary histos to hold histos from file
1645 TH1F* hEffPt[fNJetPtSlices];
1646 TH1F* hEffZ[fNJetPtSlices];
1647 TH1F* hEffXi[fNJetPtSlices];
1648
1649 for(Int_t i=0; i<fNJetPtSlices; i++) hEffPt[i] = 0;
1650 for(Int_t i=0; i<fNJetPtSlices; i++) hEffZ[i] = 0;
1651 for(Int_t i=0; i<fNJetPtSlices; i++) hEffXi[i] = 0;
1652
1653 TFile f(strfile,"READ");
1654
1655 if(!f.IsOpen()){
1656 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1657 return;
1658 }
1659
1660 if(fDebug>0) Printf("%s:%d -- read efficiencies from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
1661
1662 if(strdir && strdir.Length()) gDirectory->cd(strdir);
1663
1664 TList* list = 0;
1665
1666 if(strlist && strlist.Length()){
1667
1668 if(!(list = (TList*) gDirectory->Get(strlist))){
1669 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1670 return;
1671 }
1672 }
1673
1674 for(Int_t i=0; i<fNJetPtSlices; i++){
1675
1676 Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
1677 Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
1678
1679 TString strNameEffPt(Form("hEffPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
1680 TString strNameEffZ(Form("hEffZ_%02d_%02d",jetPtLoLim,jetPtUpLim));
1681 TString strNameEffXi(Form("hEffXi_%02d_%02d",jetPtLoLim,jetPtUpLim));
1682
1683
1684 if(list){
1685 hEffPt[i] = (TH1F*) list->FindObject(strNameEffPt);
1686 hEffZ[i] = (TH1F*) list->FindObject(strNameEffZ);
1687 hEffXi[i] = (TH1F*) list->FindObject(strNameEffXi);
1688 }
1689 else{
1690 hEffPt[i] = (TH1F*) gDirectory->Get(strNameEffPt);
1691 hEffZ[i] = (TH1F*) gDirectory->Get(strNameEffZ);
1692 hEffXi[i] = (TH1F*) gDirectory->Get(strNameEffXi);
1693 }
1694
1695 if(!hEffPt[i]){
1696 Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffPt.Data());
1697 }
1698
1699 if(!hEffZ[i]){
1700 Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffZ.Data());
1701 }
1702
1703 if(!hEffXi[i]){
1704 Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffXi.Data());
1705 }
1706
1707
1708 if(fNHistoBinsPt[i]) hEffPt[i] = (TH1F*) hEffPt[i]->Rebin(fNHistoBinsPt[i],strNameEffPt+"_rebin",fHistoBinsPt[i]->GetArray());
1709 if(fNHistoBinsZ[i]) hEffZ[i] = (TH1F*) hEffZ[i]->Rebin(fNHistoBinsZ[i],strNameEffZ+"_rebin",fHistoBinsZ[i]->GetArray());
1710 if(fNHistoBinsXi[i]) hEffXi[i] = (TH1F*) hEffXi[i]->Rebin(fNHistoBinsXi[i],strNameEffXi+"_rebin",fHistoBinsXi[i]->GetArray());
1711
1712 if(hEffPt[i]) hEffPt[i]->SetDirectory(0);
1713 if(hEffZ[i]) hEffZ[i]->SetDirectory(0);
1714 if(hEffXi[i]) hEffXi[i]->SetDirectory(0);
1715
1716 } // jet slices loop
1717
1718 f.Close();
1719
1720 for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
1721 if(hEffPt[i]) new(fh1EffPt[i]) TH1F(*hEffPt[i]);
1722 if(hEffZ[i]) new(fh1EffZ[i]) TH1F(*hEffZ[i]);
1723 if(hEffXi[i]) new(fh1EffXi[i]) TH1F(*hEffXi[i]);
1724 }
1725}
1726
1727//___________________________________________________________________________________________________________
1728void AliFragmentationFunctionCorrections::ReadBgrEfficiency(TString strfile, TString strdir, TString strlist)
1729{
1730 // read bgr eff from file
1731 // argument strlist optional - read from directory strdir if not specified
1732
1733 TH1F* hEffPtBgr[fNJetPtSlices];
1734 TH1F* hEffZBgr [fNJetPtSlices];
1735 TH1F* hEffXiBgr[fNJetPtSlices];
1736
1737 for(Int_t i=0; i<fNJetPtSlices; i++) hEffPtBgr[i] = 0;
1738 for(Int_t i=0; i<fNJetPtSlices; i++) hEffZBgr[i] = 0;
1739 for(Int_t i=0; i<fNJetPtSlices; i++) hEffXiBgr[i] = 0;
1740
1741
1742 TFile f(strfile,"READ");
1743
1744 if(!f.IsOpen()){
1745 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1746 return;
1747 }
1748
1749 if(fDebug>0) Printf("%s:%d -- read bgr efficiencies from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
1750
1751 if(strdir && strdir.Length()) gDirectory->cd(strdir);
1752
1753 TList* list = 0;
1754
1755 if(strlist && strlist.Length()){
1756
1757 if(!(list = (TList*) gDirectory->Get(strlist))){
1758 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1759 return;
1760 }
1761 }
1762
1763 for(Int_t i=0; i<fNJetPtSlices; i++){
1764
1765 Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
1766 Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
1767
1768 TString strNameEffPtBgr(Form("hEffPtBgr%02dto%02d",jetPtLoLim,jetPtUpLim));
1769 TString strNameEffZBgr(Form("hEffZBgr%02dto%02d",jetPtLoLim,jetPtUpLim));
1770 TString strNameEffXiBgr(Form("hEffXiBgr%02dto%02d",jetPtLoLim,jetPtUpLim));
1771
1772
1773 if(list){
1774 hEffPtBgr[i] = (TH1F*) list->FindObject(strNameEffPtBgr);
1775 hEffZBgr[i] = (TH1F*) list->FindObject(strNameEffZBgr);
1776 hEffXiBgr[i] = (TH1F*) list->FindObject(strNameEffXiBgr);
1777 }
1778 else{
1779 hEffPtBgr[i] = (TH1F*) gDirectory->Get(strNameEffPtBgr);
1780 hEffZBgr[i] = (TH1F*) gDirectory->Get(strNameEffZBgr);
1781 hEffXiBgr[i] = (TH1F*) gDirectory->Get(strNameEffXiBgr);
1782 }
1783
1784 if(!hEffPtBgr[i]){
1785 Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffPtBgr.Data());
1786 }
1787
1788 if(!hEffZBgr[i]){
1789 Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffZBgr.Data());
1790 }
1791
1792 if(!hEffXiBgr[i]){
1793 Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffXiBgr.Data());
1794 }
1795
1796
1797 if(fNHistoBinsPt[i]) hEffPtBgr[i] = (TH1F*) hEffPtBgr[i]->Rebin(fNHistoBinsPt[i],strNameEffPtBgr+"_rebin",fHistoBinsPt[i]->GetArray());
1798 if(fNHistoBinsZ[i]) hEffZBgr[i] = (TH1F*) hEffZBgr[i]->Rebin(fNHistoBinsZ[i],strNameEffZBgr+"_rebin",fHistoBinsZ[i]->GetArray());
1799 if(fNHistoBinsXi[i]) hEffXiBgr[i] = (TH1F*) hEffXiBgr[i]->Rebin(fNHistoBinsXi[i],strNameEffXiBgr+"_rebin",fHistoBinsXi[i]->GetArray());
1800
1801 if(hEffPtBgr[i]) hEffPtBgr[i]->SetDirectory(0);
1802 if(hEffZBgr[i]) hEffZBgr[i]->SetDirectory(0);
1803 if(hEffXiBgr[i]) hEffXiBgr[i]->SetDirectory(0);
1804
1805 } // jet slices loop
1806
1807 f.Close();
1808
1809 for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
1810 if(hEffPtBgr[i]) new(fh1EffBgrPt[i]) TH1F(*hEffPtBgr[i]);
1811 if(hEffZBgr[i]) new(fh1EffBgrZ[i]) TH1F(*hEffZBgr[i]);
1812 if(hEffXiBgr[i]) new(fh1EffBgrXi[i]) TH1F(*hEffXiBgr[i]);
1813 }
1814}
1815
1816// ________________________________________________
1817void AliFragmentationFunctionCorrections::EffCorr()
1818{
1819 // apply efficiency correction
1820
1821 AddCorrectionLevel("eff");
1822
1823 for(Int_t i=0; i<fNJetPtSlices; i++){
1824
1825 TH1F* histPt = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i);
1826 TH1F* histZ = fCorrFF[fNCorrectionLevels-2]->GetZ(i);
1827 TH1F* histXi = fCorrFF[fNCorrectionLevels-2]->GetXi(i);
1828
1829 TString histNamePt = histPt->GetName();
1830 TString histNameZ = histZ->GetName();
1831 TString histNameXi = histXi->GetName();
1832
1833
1834 TH1F* hFFTrackPtEffCorr = (TH1F*) histPt->Clone(histNamePt);
1835 hFFTrackPtEffCorr->Divide(histPt,fh1EffPt[i],1,1,"");
1836
1837 TH1F* hFFZEffCorr = (TH1F*) histZ->Clone(histNameZ);
1838 hFFZEffCorr->Divide(histZ,fh1EffZ[i],1,1,"");
1839
1840 TH1F* hFFXiEffCorr = (TH1F*) histXi->Clone(histNameXi);
1841 hFFXiEffCorr->Divide(histXi,fh1EffXi[i],1,1,"");
1842
1843 fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hFFTrackPtEffCorr,hFFZEffCorr,hFFXiEffCorr);
1844 }
1845}
1846
1847//___________________________________________________
1848void AliFragmentationFunctionCorrections::EffCorrBgr()
1849{
1850 // apply efficiency correction to bgr distributions
1851
1852 AddCorrectionLevelBgr("eff");
1853
1854 Printf("%s:%d -- apply efficiency correction, corrLevel %d",(char*)__FILE__,__LINE__,fNCorrectionLevels-1);
1855
1856 for(Int_t i=0; i<fNJetPtSlices; i++){
1857
1858 TH1F* histPt = fCorrBgr[fNCorrectionLevelsBgr-2]->GetTrackPt(i);
1859 TH1F* histZ = fCorrBgr[fNCorrectionLevelsBgr-2]->GetZ(i);
1860 TH1F* histXi = fCorrBgr[fNCorrectionLevelsBgr-2]->GetXi(i);
1861
1862 TString histNamePt = histPt->GetName();
1863 TString histNameZ = histZ->GetName();
1864 TString histNameXi = histXi->GetName();
1865
1866
1867 TH1F* hFFTrackPtEffCorr = (TH1F*) histPt->Clone(histNamePt);
1868 hFFTrackPtEffCorr->Divide(histPt,fh1EffPt[i],1,1,"");
1869
1870 TH1F* hFFZEffCorr = (TH1F*) histZ->Clone(histNameZ);
1871 hFFZEffCorr->Divide(histZ,fh1EffZ[i],1,1,"");
1872
1873 TH1F* hFFXiEffCorr = (TH1F*) histXi->Clone(histNameXi);
1874 hFFXiEffCorr->Divide(histXi,fh1EffXi[i],1,1,"");
1875
1876 fCorrBgr[fNCorrectionLevelsBgr-1]->AddCorrHistos(i,hFFTrackPtEffCorr,hFFZEffCorr,hFFXiEffCorr);
1877 }
1878}
1879
1880//______________________________________________________________________
1881void AliFragmentationFunctionCorrections::XiShift(const Int_t corrLevel)
1882{
1883 // re-evaluate jet energy after FF corrections from dN/dpt distribution
1884 // apply correction (shift) to dN/dxi distribution: xi = ln (pt/E) -> xi' = ln (pt/E') = ln (pt/E x E/E') = xi + ln E/E'
1885 // argument corrlevel: which level of correction to be corrected/shifted to
1886
1887 if(corrLevel>=fNCorrectionLevels){
1888 Printf(" calc xi shift: corrLevel exceeded - do nothing");
1889 return;
1890 }
1891
1892 Double_t* jetPtUncorr = new Double_t[fNJetPtSlices];
1893 Double_t* jetPtCorr = new Double_t[fNJetPtSlices];
1894 Double_t* deltaXi = new Double_t[fNJetPtSlices];
1895
1896 for(Int_t i=0; i<fNJetPtSlices; i++){
1897
1898 TH1F* histPtRaw = fCorrFF[0]->GetTrackPt(i);
1899 TH1F* histPt = fCorrFF[corrLevel]->GetTrackPt(i);
1900
1901 Double_t ptUncorr = 0;
1902 Double_t ptCorr = 0;
1903
1904 for(Int_t bin = 1; bin<=histPtRaw->GetNbinsX(); bin++){
1905
1906 Double_t cont = histPtRaw->GetBinContent(bin);
1907 Double_t width = histPtRaw->GetBinWidth(bin);
1908 Double_t meanPt = histPtRaw->GetBinCenter(bin);
1909
1910 ptUncorr += meanPt*cont*width;
1911 }
1912
1913 for(Int_t bin = 1; bin<=histPt->GetNbinsX(); bin++){
1914
1915 Double_t cont = histPt->GetBinContent(bin);
1916 Double_t width = histPt->GetBinWidth(bin);
1917 Double_t meanPt = histPt->GetBinCenter(bin);
1918
1919 ptCorr += meanPt*cont*width;
1920 }
1921
1922 jetPtUncorr[i] = ptUncorr;
1923 jetPtCorr[i] = ptCorr;
1924 }
1925
1926 // calc dXi from dN/dpt distribution :
1927 // sum over track pt for raw and corrected FF is equivalent to raw/corrected jet pt
1928
1929 for(Int_t i=0; i<fNJetPtSlices; i++){
1930
1931 Float_t jetPtLoLim = fJetPtSlices->At(i);
1932 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1933
1934 Double_t meanJetPt = 0.5*(jetPtUpLim+jetPtLoLim);
1935
1936 Double_t ptUncorr = jetPtUncorr[i];
1937 Double_t ptCorr = jetPtCorr[i];
1938
1939 Double_t dXi = TMath::Log(ptCorr/ptUncorr);
1940
1941 Printf(" calc xi shift: jet pt slice %d, mean jet pt %f, ptUncorr %f, ptCorr %f, ratio corr/uncorr %f, dXi %f "
1942 ,i,meanJetPt,ptUncorr,ptCorr,ptCorr/ptUncorr,dXi);
1943
1944 deltaXi[i] = dXi;
1945 }
1946
1947 // book & fill new dN/dxi histos
1948
1949 for(Int_t i=0; i<fNJetPtSlices; i++){
1950
1951 TH1F* histXi = fCorrFF[corrLevel]->GetXi(i);
1952
1953 Double_t dXi = deltaXi[i];
1954
1955 Int_t nBins = histXi->GetNbinsX();
1956 const Double_t* binsVec = histXi->GetXaxis()->GetXbins()->GetArray();
1957 Float_t binsVecNew[nBins+1];
1958
1959 TString strName = histXi->GetName();
1960 strName.Append("_shift");
1961 TString strTit = histXi->GetTitle();
1962
1963 TH1F* hXiCorr;
1964
1965 // create shifted histo ...
1966
1967 if(binsVec){ // binsVec only neq NULL if histo was rebinned before
1968
1969 for(Int_t bin=0; bin<nBins+1; bin++) binsVecNew[bin] = binsVec[bin] + dXi;
1970 hXiCorr = new TH1F(strName,strTit,nBins,binsVecNew);
1971 }
1972 else{ // uniform bin size
1973
1974 Double_t xMin = histXi->GetXaxis()->GetXmin();
1975 Double_t xMax = histXi->GetXaxis()->GetXmax();
1976
1977 xMin += dXi;
1978 xMax += dXi;
1979
1980 hXiCorr = new TH1F(strName,strTit,nBins,xMin,xMax);
1981 }
1982
1983 // ... and fill
1984
1985 for(Int_t bin=1; bin<nBins+1; bin++){
1986 Double_t cont = histXi->GetBinContent(bin);
1987 Double_t err = histXi->GetBinError(bin);
1988
1989 hXiCorr->SetBinContent(bin,cont);
1990 hXiCorr->SetBinError(bin,err);
1991 }
1992
1993 new(fh1FFXiShift[i]) TH1F(*hXiCorr);
1994 delete hXiCorr;
1995 }
1996
1997 delete[] jetPtUncorr;
1998 delete[] jetPtCorr;
1999 delete[] deltaXi;
2000}
2001
2002//_____________________________________________________
2003void AliFragmentationFunctionCorrections::SubtractBgr()
2004{
2005 // subtract bgr distribution from FF
2006 // the current corr level is used for both
2007
2008 AddCorrectionLevel("bgrSub");
2009
2010 for(Int_t i=0; i<fNJetPtSlices; i++){
2011
2012 TH1F* histPt = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i);
2013 TH1F* histZ = fCorrFF[fNCorrectionLevels-2]->GetZ(i);
2014 TH1F* histXi = fCorrFF[fNCorrectionLevels-2]->GetXi(i);
2015
2016 TH1F* histPtBgr = fCorrBgr[fNCorrectionLevelsBgr-1]->GetTrackPt(i);
2017 TH1F* histZBgr = fCorrBgr[fNCorrectionLevelsBgr-1]->GetZ(i);
2018 TH1F* histXiBgr = fCorrBgr[fNCorrectionLevelsBgr-1]->GetXi(i);
2019
2020 TString histNamePt = histPt->GetName();
2021 TString histNameZ = histZ->GetName();
2022 TString histNameXi = histXi->GetName();
2023
2024
2025 TH1F* hFFTrackPtBgrSub = (TH1F*) histPt->Clone(histNamePt);
2026 hFFTrackPtBgrSub->Add(histPtBgr,-1);
2027
2028 TH1F* hFFZBgrSub = (TH1F*) histZ->Clone(histNameZ);
2029 hFFZBgrSub->Add(histZBgr,-1);
2030
2031 TH1F* hFFXiBgrSub = (TH1F*) histXi->Clone(histNameXi);
2032 hFFXiBgrSub->Add(histXiBgr,-1);
2033
2034 fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hFFTrackPtBgrSub,hFFZBgrSub,hFFXiBgrSub);
2035 }
2036}
2037
2038//________________________________________________________________________________________________________________
2039void AliFragmentationFunctionCorrections::WriteSingleTrackEff(TString strInfile, TString strID, TString strOutfile,
2040 Bool_t updateOutfile, TString strOutDir,TString strPostfix)
2041{
2042 // read task ouput from MC and write single track eff - standard dir/list
2043
2044 TString strdir = "PWG4_FragmentationFunction_" + strID;
2045 TString strlist = "fracfunc_" + strID;
2046
2047 WriteSingleTrackEff(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir,strPostfix);
2048}
2049
2050//___________________________________________________________________________________________________________________________________
2051void AliFragmentationFunctionCorrections::WriteSingleTrackEff(TString strInfile, TString strdir, TString strlist,
2052 TString strOutfile, Bool_t updateOutfile, TString strOutDir,TString strPostfix)
2053{
2054 // read task output from MC and write single track eff as function of pt, eta, phi
2055 // argument strlist optional - read from directory strdir if not specified
2056 // write eff to file stroutfile - by default only 'update' file (don't overwrite)
2057
2058
2059 TH1D* hdNdptTracksMCPrimGen;
2060 TH2D* hdNdetadphiTracksMCPrimGen;
2061
2062 TH1D* hdNdptTracksMCPrimRec;
2063 TH2D* hdNdetadphiTracksMCPrimRec;
2064
2065
2066 TFile f(strInfile,"READ");
2067
2068 if(!f.IsOpen()){
2069 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2070 return;
2071 }
2072
2073 if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: open task ouput file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
2074
2075 if(strdir && strdir.Length()){
2076 if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: change dir to %s",(char*)__FILE__,__LINE__,strdir.Data());
2077 gDirectory->cd(strdir);
2078 }
2079
2080 TList* list = 0;
2081
2082 if(strlist && strlist.Length()){
2083
2084 if(!(list = (TList*) gDirectory->Get(strlist))){
2085 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2086 return;
2087 }
2088 }
2089
2090
2091 TString hnamePtRecEffGen = "fh1TrackQAPtRecEffGen";
2092 if(strPostfix.Length()) hnamePtRecEffGen.Form("fh1TrackQAPtRecEffGen_%s",strPostfix.Data());
2093
2094 TString hnamePtRecEffRec = "fh1TrackQAPtRecEffRec";
2095 if(strPostfix.Length()) hnamePtRecEffRec.Form("fh1TrackQAPtRecEffRec_%s",strPostfix.Data());
2096
2097 TString hnameEtaPhiRecEffGen = "fh2TrackQAEtaPhiRecEffGen";
2098 if(strPostfix.Length()) hnameEtaPhiRecEffGen.Form("fh2TrackQAEtaPhiRecEffGen_%s",strPostfix.Data());
2099
2100 TString hnameEtaPhiRecEffRec = "fh2TrackQAEtaPhiRecEffRec";
2101 if(strPostfix.Length()) hnameEtaPhiRecEffRec.Form("fh2TrackQAEtaPhiRecEffRec_%s",strPostfix.Data());
2102
2103
2104 if(list){
2105 hdNdptTracksMCPrimGen = (TH1D*) list->FindObject(hnamePtRecEffGen);
2106 hdNdetadphiTracksMCPrimGen = (TH2D*) list->FindObject(hnameEtaPhiRecEffGen);
2107
2108 hdNdptTracksMCPrimRec = (TH1D*) list->FindObject(hnamePtRecEffRec);
2109 hdNdetadphiTracksMCPrimRec = (TH2D*) list->FindObject(hnameEtaPhiRecEffRec);
2110 }
2111 else{
2112 hdNdptTracksMCPrimGen = (TH1D*) gDirectory->Get(hnamePtRecEffGen);
2113 hdNdetadphiTracksMCPrimGen = (TH2D*) gDirectory->Get(hnameEtaPhiRecEffGen);
2114
2115 hdNdptTracksMCPrimRec = (TH1D*) gDirectory->Get(hnamePtRecEffRec);
2116 hdNdetadphiTracksMCPrimRec = (TH2D*) gDirectory->Get(hnameEtaPhiRecEffRec);
2117 }
2118
2119 hdNdptTracksMCPrimGen->SetDirectory(0);
2120 hdNdetadphiTracksMCPrimGen->SetDirectory(0);
2121 hdNdptTracksMCPrimRec->SetDirectory(0);
2122 hdNdetadphiTracksMCPrimRec->SetDirectory(0);
2123
2124 f.Close();
2125
2126 // projections: dN/deta, dN/dphi
2127
2128 TH1D* hdNdetaTracksMCPrimGen = (TH1D*) hdNdetadphiTracksMCPrimGen->ProjectionX("hdNdetaTracksMcPrimGen");
2129 TH1D* hdNdphiTracksMCPrimGen = (TH1D*) hdNdetadphiTracksMCPrimGen->ProjectionY("hdNdphiTracksMcPrimGen");
2130
2131 TH1D* hdNdetaTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionX("hdNdetaTracksMcPrimRec");
2132 TH1D* hdNdphiTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionY("hdNdphiTracksMcPrimRec");
2133
2134 // rebin
2135
2136 TString strNamePtGen = "hTrackPtGenPrim";
2137 TString strNamePtRec = "hTrackPtRecPrim";
2138
2139 if(fNHistoBinsSinglePt) hdNdptTracksMCPrimGen = (TH1D*) hdNdptTracksMCPrimGen->Rebin(fNHistoBinsSinglePt,strNamePtGen,fHistoBinsSinglePt->GetArray());
2140 if(fNHistoBinsSinglePt) hdNdptTracksMCPrimRec = (TH1D*) hdNdptTracksMCPrimRec->Rebin(fNHistoBinsSinglePt,strNamePtRec,fHistoBinsSinglePt->GetArray());
2141
2142 // efficiency: divide
2143
2144 TString hNameTrackEffPt = "hSingleTrackEffPt";
2145 if(strPostfix.Length()) hNameTrackEffPt.Form("hSingleTrackEffPt_%s",strPostfix.Data());
2146
2147 TString hNameTrackEffEta = "hSingleTrackEffEta";
2148 if(strPostfix.Length()) hNameTrackEffEta.Form("hSingleTrackEffEta_%s",strPostfix.Data());
2149
2150 TString hNameTrackEffPhi = "hSingleTrackEffPhi";
2151 if(strPostfix.Length()) hNameTrackEffPhi.Form("hSingleTrackEffPhi_%s",strPostfix.Data());
2152
2153
2154 TH1F* hSingleTrackEffPt = (TH1F*) hdNdptTracksMCPrimRec->Clone(hNameTrackEffPt);
2155 hSingleTrackEffPt->Divide(hdNdptTracksMCPrimRec,hdNdptTracksMCPrimGen,1,1,"B"); // binominal errors
2156
2157 TH1F* hSingleTrackEffEta = (TH1F*) hdNdetaTracksMCPrimRec->Clone(hNameTrackEffEta);
2158 hSingleTrackEffEta->Divide(hdNdetaTracksMCPrimRec,hdNdetaTracksMCPrimGen,1,1,"B"); // binominal errors
2159
2160 TH1F* hSingleTrackEffPhi = (TH1F*) hdNdphiTracksMCPrimRec->Clone(hNameTrackEffPhi);
2161 hSingleTrackEffPhi->Divide(hdNdphiTracksMCPrimRec,hdNdphiTracksMCPrimGen,1,1,"B"); // binominal errors
2162
2163
2164 TString outfileOption = "RECREATE";
2165 if(updateOutfile) outfileOption = "UPDATE";
2166
2167 TFile out(strOutfile,outfileOption);
2168
2169 if(!out.IsOpen()){
2170 Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2171 return;
2172 }
2173
2174 if(fDebug>0) Printf("%s:%d -- write efficiency to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2175
2176 if(strOutDir && strOutDir.Length()){
2177
2178 TDirectory* dir;
2179 if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd();
2180 else{
2181 dir = out.mkdir(strOutDir);
2182 dir->cd();
2183 }
2184 }
2185
2186 hSingleTrackEffPt->Write();
2187 hSingleTrackEffEta->Write();
2188 hSingleTrackEffPhi->Write();
2189
2190 out.Close();
2191}
2192
2193//________________________________________________________________________________________________________________
2194void AliFragmentationFunctionCorrections::WriteSingleTrackSecCorr(TString strInfile, TString strID, TString strOutfile,
2195 Bool_t updateOutfile, TString strOutDir)
2196{
2197 // read task ouput from MC and write single track eff - standard dir/list
2198
2199 TString strdir = "PWG4_FragmentationFunction_" + strID;
2200 TString strlist = "fracfunc_" + strID;
2201
2202 WriteSingleTrackSecCorr(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir);
2203}
2204
2205//___________________________________________________________________________________________________________________________________
2206void AliFragmentationFunctionCorrections::WriteSingleTrackSecCorr(TString strInfile, TString strdir, TString strlist,
2207 TString strOutfile, Bool_t updateOutfile, TString strOutDir)
2208{
2209 // read task output from MC and write single track secondaries contamination correction as function of pt, eta, phi
2210 // argument strlist optional - read from directory strdir if not specified
2211 // write corr factor to file stroutfile - by default only 'update' file (don't overwrite)
2212
2213 TH1D* hdNdptTracksMCPrimRec;
2214 TH2D* hdNdetadphiTracksMCPrimRec;
2215
2216 TH1D* hdNdptTracksMCSecRec;
2217 TH2D* hdNdetadphiTracksMCSecRec;
2218
2219 TFile f(strInfile,"READ");
2220
2221 if(!f.IsOpen()){
2222 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2223 return;
2224 }
2225
2226 if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
2227
2228 if(strdir && strdir.Length()) gDirectory->cd(strdir);
2229
2230 TList* list = 0;
2231
2232 if(strlist && strlist.Length()){
2233
2234 if(!(list = (TList*) gDirectory->Get(strlist))){
2235 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2236 return;
2237 }
2238 }
2239
2240
2241 if(list){
2242 hdNdptTracksMCPrimRec = (TH1D*) list->FindObject("fh1TrackQAPtRecEffGen");
2243 hdNdetadphiTracksMCPrimRec = (TH2D*) list->FindObject("fh2TrackQAEtaPhiRecEffGen");
2244
2245 hdNdptTracksMCSecRec = (TH1D*) list->FindObject("fh1TrackQAPtSecRec");
2246 hdNdetadphiTracksMCSecRec = (TH2D*) list->FindObject("fh2TrackQAEtaPhiSecRec");
2247 }
2248 else{
2249 hdNdptTracksMCPrimRec = (TH1D*) gDirectory->Get("fh1TrackQAPtRecEffGen");
2250 hdNdetadphiTracksMCPrimRec = (TH2D*) gDirectory->Get("fh2TrackQAEtaPhiRecEffGen");
2251
2252 hdNdptTracksMCSecRec = (TH1D*) gDirectory->Get("fh1TrackQAPtSecRec");
2253 hdNdetadphiTracksMCSecRec = (TH2D*) gDirectory->Get("fh2TrackQAEtaPhiSecRec");
2254 }
2255
2256 hdNdptTracksMCPrimRec->SetDirectory(0);
2257 hdNdetadphiTracksMCPrimRec->SetDirectory(0);
2258 hdNdptTracksMCSecRec->SetDirectory(0);
2259 hdNdetadphiTracksMCSecRec->SetDirectory(0);
2260
2261 f.Close();
2262
2263 // projections: dN/deta, dN/dphi
2264
2265 TH1D* hdNdetaTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionX("hdNdetaTracksMcPrimRec");
2266 TH1D* hdNdphiTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionY("hdNdphiTracksMcPrimRec");
2267
2268 TH1D* hdNdetaTracksMCSecRec = (TH1D*) hdNdetadphiTracksMCSecRec->ProjectionX("hdNdetaTracksMcSecRec");
2269 TH1D* hdNdphiTracksMCSecRec = (TH1D*) hdNdetadphiTracksMCSecRec->ProjectionY("hdNdphiTracksMcSecRec");
2270
2271
2272 // rebin
2273
2274 TString strNamePtPrim = "hTrackPtPrim";
2275 TString strNamePtSec = "hTrackPtSec";
2276
2277 if(fNHistoBinsSinglePt) hdNdptTracksMCPrimRec = (TH1D*) hdNdptTracksMCPrimRec->Rebin(fNHistoBinsSinglePt,strNamePtPrim,fHistoBinsSinglePt->GetArray());
2278 if(fNHistoBinsSinglePt) hdNdptTracksMCSecRec = (TH1D*) hdNdptTracksMCSecRec->Rebin(fNHistoBinsSinglePt,strNamePtSec,fHistoBinsSinglePt->GetArray());
2279
2280
2281 // secondary correction factor: divide prim/(prim+sec)
2282
2283 TH1F* hSingleTrackSecCorrPt = (TH1F*) hdNdptTracksMCSecRec->Clone("hSingleTrackSecCorrPt");
2284 TH1F* hSumPrimSecPt = (TH1F*) hdNdptTracksMCSecRec->Clone("hSumPrimSecPt");
2285 hSumPrimSecPt->Add(hdNdptTracksMCPrimRec);
2286 hSingleTrackSecCorrPt->Divide(hdNdptTracksMCPrimRec,hSumPrimSecPt,1,1,"B"); // binominal errors
2287
2288 TH1F* hSingleTrackSecCorrEta = (TH1F*) hdNdetaTracksMCSecRec->Clone("hSingleTrackSecCorrEta");
2289 TH1F* hSumPrimSecEta = (TH1F*) hdNdetaTracksMCSecRec->Clone("hSumPrimSecEta");
2290 hSumPrimSecEta->Add(hdNdetaTracksMCPrimRec);
2291 hSingleTrackSecCorrEta->Divide(hdNdetaTracksMCPrimRec,hSumPrimSecEta,1,1,"B"); // binominal errors
2292
2293 TH1F* hSingleTrackSecCorrPhi = (TH1F*) hdNdphiTracksMCSecRec->Clone("hSingleTrackSecCorrPhi");
2294 TH1F* hSumPrimSecPhi = (TH1F*) hdNdphiTracksMCSecRec->Clone("hSumPrimSecPhi");
2295 hSumPrimSecPhi->Add(hdNdphiTracksMCPrimRec);
2296 hSingleTrackSecCorrPhi->Divide(hdNdphiTracksMCPrimRec,hSumPrimSecPhi,1,1,"B"); // binominal errors
2297
2298
2299 TString outfileOption = "RECREATE";
2300 if(updateOutfile) outfileOption = "UPDATE";
2301
2302 TFile out(strOutfile,outfileOption);
2303
2304 if(!out.IsOpen()){
2305 Printf("%s:%d -- error opening secCorr output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2306 return;
2307 }
2308
2309 if(fDebug>0) Printf("%s:%d -- write secCorr to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2310
2311 if(strOutDir && strOutDir.Length()){
2312
2313 TDirectory* dir;
2314 if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd();
2315 else{
2316 dir = out.mkdir(strOutDir);
2317 dir->cd();
2318 }
2319 }
2320
2321 hdNdptTracksMCSecRec->Write();
2322 hdNdetaTracksMCSecRec->Write();
2323 hdNdphiTracksMCSecRec->Write();
2324
2325 hSingleTrackSecCorrPt->Write();
2326 hSingleTrackSecCorrEta->Write();
2327 hSingleTrackSecCorrPhi->Write();
2328
2329 out.Close();
2330}
2331
2332//________________________________________________________________________________________________________________
2333void AliFragmentationFunctionCorrections::WriteSingleResponse(TString strInfile, TString strID, TString strOutfile,
2334 Bool_t updateOutfile, TString strOutDir)
2335{
2336 // read task ouput from MC and write single track eff - standard dir/list
2337
2338 TString strdir = "PWG4_FragmentationFunction_" + strID;
2339 TString strlist = "fracfunc_" + strID;
2340
2341 WriteSingleResponse(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir);
2342}
2343
2344
2345//_____________________________________________________________________________________________________________________________________
2346void AliFragmentationFunctionCorrections::WriteSingleResponse(TString strInfile, TString strdir, TString strlist,
2347 TString strOutfile, Bool_t updateOutfile, TString strOutDir)
2348{
2349 // read 2d THnSparse response matrices in pt from file
2350 // project TH2
2351 // write to strOutfile
2352
2353 THnSparse* hnResponseSinglePt;
2354 TH2F* h2ResponseSinglePt;
2355
2356 TFile f(strInfile,"READ");
2357
2358 if(!f.IsOpen()){
2359 Printf("%s:%d -- error opening file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2360 return;
2361 }
2362
2363 if(fDebug>0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
2364
2365 if(strdir && strdir.Length()) gDirectory->cd(strdir);
2366
2367 TList* list = 0;
2368
2369 if(strlist && strlist.Length()){
2370
2371 if(!(list = (TList*) gDirectory->Get(strlist))){
2372 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2373 return;
2374 }
2375 }
2376
2377 if(list) hnResponseSinglePt = (THnSparse*) list->FindObject("fhnResponseSinglePt");
2378 else hnResponseSinglePt = (THnSparse*) gDirectory->Get("fhnResponseSinglePt");
2379
2380
2381 if(!hnResponseSinglePt){
2382 Printf("%s:%d -- error retrieving response matrix fhnResponseSinglePt",(char*)__FILE__,__LINE__);
2383 return;
2384 }
2385
2386 f.Close();
2387
2388
2389 // project 2d histo
2390 h2ResponseSinglePt = (TH2F*) hnResponseSinglePt->Projection(1,0);// note convention: yDim,xDim
2391 h2ResponseSinglePt->SetNameTitle("h2ResponseSinglePt","");
2392
2393
2394 // write
2395
2396 TString outfileOption = "RECREATE";
2397 if(updateOutfile) outfileOption = "UPDATE";
2398
2399 TFile out(strOutfile,outfileOption);
2400
2401 if(!out.IsOpen()){
2402 Printf("%s:%d -- error opening response matrix output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2403 return;
2404 }
2405
2406 if(fDebug>0) Printf("%s:%d -- write response matrices to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2407
2408 if(strOutDir && strOutDir.Length()){
2409
2410 TDirectory* dir;
2411 if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd();
2412 else{
2413 dir = out.mkdir(strOutDir);
2414 dir->cd();
2415 }
2416 }
2417
2418 hnResponseSinglePt->Write();
2419 h2ResponseSinglePt->Write();
2420
2421 out.Close();
2422}
2423
2424//________________________________________________________________________________________________________________
2425void AliFragmentationFunctionCorrections::WriteJetTrackEff(TString strInfile, TString strID, TString strOutfile,
2426 Bool_t updateOutfile)
2427{
2428 // read task ouput from MC and write single track eff - standard dir/list
2429
2430 TString strdir = "PWG4_FragmentationFunction_" + strID;
2431 TString strlist = "fracfunc_" + strID;
2432
2433 WriteJetTrackEff(strInfile,strdir,strlist,strOutfile,updateOutfile);
2434}
2435
2436//___________________________________________________________________________________________________________________________________
2437void AliFragmentationFunctionCorrections::WriteJetTrackEff(TString strInfile, TString strdir, TString strlist,
2438 TString strOutfile, Bool_t updateOutfile)
2439{
2440 // read task output from MC and write track eff in jet pt slices as function of pt, z, xi
2441 // argument strlist optional - read from directory strdir if not specified
2442 // write eff to file strOutfile - by default only 'update' file (don't overwrite)
2443
2444 TH1F* hEffPt[fNJetPtSlices];
2445 TH1F* hEffXi[fNJetPtSlices];
2446 TH1F* hEffZ[fNJetPtSlices];
2447
2448 TH1F* hdNdptTracksMCPrimGen[fNJetPtSlices];
2449 TH1F* hdNdxiMCPrimGen[fNJetPtSlices];
2450 TH1F* hdNdzMCPrimGen[fNJetPtSlices];
2451
2452 TH1F* hdNdptTracksMCPrimRec[fNJetPtSlices];
2453 TH1F* hdNdxiMCPrimRec[fNJetPtSlices];
2454 TH1F* hdNdzMCPrimRec[fNJetPtSlices];
2455
2456
2457 TH1F* fh1FFJetPtRecEffGen;
2458
2459 TH2F* fh2FFTrackPtRecEffGen;
2460 TH2F* fh2FFZRecEffGen;
2461 TH2F* fh2FFXiRecEffGen;
2462
2463 TH2F* fh2FFTrackPtRecEffRec;
2464 TH2F* fh2FFZRecEffRec;
2465 TH2F* fh2FFXiRecEffRec;
2466
2467
2468 TFile f(strInfile,"READ");
2469
2470 if(!f.IsOpen()){
2471 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2472 return;
2473 }
2474
2475 if(fDebug>0) Printf("%s:%d -- writeJetTrackEff: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
2476
2477 if(strdir && strdir.Length()) gDirectory->cd(strdir);
2478
2479 TList* list = 0;
2480
2481 if(strlist && strlist.Length()){
2482
2483 if(!(list = (TList*) gDirectory->Get(strlist))){
2484 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2485 return;
2486 }
2487 }
2488
2489 if(list){
2490 fh1FFJetPtRecEffGen = (TH1F*) list->FindObject("fh1FFJetPtRecEffGen");
2491
2492 fh2FFTrackPtRecEffGen = (TH2F*) list->FindObject("fh2FFTrackPtRecEffGen");
2493 fh2FFZRecEffGen = (TH2F*) list->FindObject("fh2FFZRecEffGen");
2494 fh2FFXiRecEffGen = (TH2F*) list->FindObject("fh2FFXiRecEffGen");
2495
2496 fh2FFTrackPtRecEffRec = (TH2F*) list->FindObject("fh2FFTrackPtRecEffRec");
2497 fh2FFZRecEffRec = (TH2F*) list->FindObject("fh2FFZRecEffRec");
2498 fh2FFXiRecEffRec = (TH2F*) list->FindObject("fh2FFXiRecEffRec");
2499 }
2500 else{
2501 fh1FFJetPtRecEffGen = (TH1F*) gDirectory->FindObject("fh1FFJetPtRecEffGen");
2502
2503 fh2FFTrackPtRecEffGen = (TH2F*) gDirectory->FindObject("fh2FFTrackPtRecEffGen");
2504 fh2FFZRecEffGen = (TH2F*) gDirectory->FindObject("fh2FFZRecEffGen");
2505 fh2FFXiRecEffGen = (TH2F*) gDirectory->FindObject("fh2FFXiRecEffGen");
2506
2507 fh2FFTrackPtRecEffRec = (TH2F*) gDirectory->FindObject("fh2FFTrackPtRecEffRec");
2508 fh2FFZRecEffRec = (TH2F*) gDirectory->FindObject("fh2FFZRecEffRec");
2509 fh2FFXiRecEffRec = (TH2F*) gDirectory->FindObject("fh2FFXiRecEffRec");
2510 }
2511
2512 fh1FFJetPtRecEffGen->SetDirectory(0);
2513
2514 fh2FFTrackPtRecEffGen->SetDirectory(0);
2515 fh2FFZRecEffGen->SetDirectory(0);
2516 fh2FFXiRecEffGen->SetDirectory(0);
2517
2518 fh2FFTrackPtRecEffRec->SetDirectory(0);
2519 fh2FFZRecEffRec->SetDirectory(0);
2520 fh2FFXiRecEffRec->SetDirectory(0);
2521
2522 f.Close();
2523
2524
2525 // projections: FF for generated and reconstructed primaries
2526
2527 for(Int_t i=0; i<fNJetPtSlices; i++){
2528
2529 Float_t jetPtLoLim = fJetPtSlices->At(i);
2530 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
2531
2532 Int_t binLo = static_cast<Int_t>(fh2FFTrackPtRecEffGen->GetXaxis()->FindBin(jetPtLoLim));
2533 Int_t binUp = static_cast<Int_t>(fh2FFTrackPtRecEffGen->GetXaxis()->FindBin(jetPtUpLim))-1;
2534
2535 if(binUp > fh2FFTrackPtRecEffGen->GetNbinsX()){
2536 Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
2537 return;
2538 }
2539
2540 TString strNameFFPtGen(Form("fh1FFTrackPtGenPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2541 TString strNameFFZGen(Form("fh1FFZGenPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2542 TString strNameFFXiGen(Form("fh1FFXiGenPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2543
2544 TString strNameFFPtRec(Form("fh1FFTrackPtRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2545 TString strNameFFZRec(Form("fh1FFZRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2546 TString strNameFFXiRec(Form("fh1FFXiRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2547
2548 // project
2549 // appendix 'unbinned' to avoid histos with same name after rebinning
2550
2551 hdNdptTracksMCPrimGen[i] = (TH1F*) fh2FFTrackPtRecEffGen->ProjectionY(strNameFFPtGen+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
2552 hdNdzMCPrimGen[i] = (TH1F*) fh2FFZRecEffGen->ProjectionY(strNameFFZGen+"_unBinned",binLo,binUp,"o");
2553 hdNdxiMCPrimGen[i] = (TH1F*) fh2FFXiRecEffGen->ProjectionY(strNameFFXiGen+"_unBinned",binLo,binUp,"o");
2554
2555 hdNdptTracksMCPrimRec[i] = (TH1F*) fh2FFTrackPtRecEffRec->ProjectionY(strNameFFPtRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
2556 hdNdzMCPrimRec[i] = (TH1F*) fh2FFZRecEffRec->ProjectionY(strNameFFZRec+"_unBinned",binLo,binUp,"o");
2557 hdNdxiMCPrimRec[i] = (TH1F*) fh2FFXiRecEffRec->ProjectionY(strNameFFXiRec+"_unBinned",binLo,binUp,"o");
2558
2559 // rebin
2560
2561 if(fNHistoBinsPt[i]) hdNdptTracksMCPrimGen[i] = (TH1F*) hdNdptTracksMCPrimGen[i]->Rebin(fNHistoBinsPt[i],strNameFFPtGen,fHistoBinsPt[i]->GetArray());
2562 if(fNHistoBinsZ[i]) hdNdzMCPrimGen[i] = (TH1F*) hdNdzMCPrimGen[i]->Rebin(fNHistoBinsZ[i],strNameFFZGen,fHistoBinsZ[i]->GetArray());
2563 if(fNHistoBinsXi[i]) hdNdxiMCPrimGen[i] = (TH1F*) hdNdxiMCPrimGen[i]->Rebin(fNHistoBinsXi[i],strNameFFXiGen,fHistoBinsXi[i]->GetArray());
2564
2565 if(fNHistoBinsPt[i]) hdNdptTracksMCPrimRec[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtRec,fHistoBinsPt[i]->GetArray());
2566 if(fNHistoBinsZ[i]) hdNdzMCPrimRec[i] = (TH1F*) hdNdzMCPrimRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZRec,fHistoBinsZ[i]->GetArray());
2567 if(fNHistoBinsXi[i]) hdNdxiMCPrimRec[i] = (TH1F*) hdNdxiMCPrimRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiRec,fHistoBinsXi[i]->GetArray());
2568
2569 hdNdptTracksMCPrimGen[i]->SetNameTitle(strNameFFPtGen,"");
2570 hdNdzMCPrimGen[i]->SetNameTitle(strNameFFZGen,"");
2571 hdNdxiMCPrimGen[i]->SetNameTitle(strNameFFXiGen,"");
2572
2573 hdNdptTracksMCPrimRec[i]->SetNameTitle(strNameFFPtRec,"");
2574 hdNdzMCPrimRec[i]->SetNameTitle(strNameFFZRec,"");
2575 hdNdxiMCPrimRec[i]->SetNameTitle(strNameFFXiRec,"");
2576
2577 // normalize
2578
2579 Double_t nJetsBin = fh1FFJetPtRecEffGen->Integral(binLo,binUp);
2580
2581 NormalizeTH1(hdNdptTracksMCPrimGen[i],nJetsBin);
2582 NormalizeTH1(hdNdzMCPrimGen[i],nJetsBin);
2583 NormalizeTH1(hdNdxiMCPrimGen[i],nJetsBin);
2584
2585 NormalizeTH1(hdNdptTracksMCPrimRec[i],nJetsBin);
2586 NormalizeTH1(hdNdzMCPrimRec[i],nJetsBin);
2587 NormalizeTH1(hdNdxiMCPrimRec[i],nJetsBin);
2588
2589 // divide rec/gen : efficiency
2590
2591 TString strNameEffPt(Form("hEffPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2592 TString strNameEffZ(Form("hEffZ_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2593 TString strNameEffXi(Form("hEffXi_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2594
2595 hEffPt[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Clone(strNameEffPt);
2596 hEffPt[i]->Divide(hdNdptTracksMCPrimRec[i],hdNdptTracksMCPrimGen[i],1,1,"B"); // binominal errors
2597
2598 hEffXi[i] = (TH1F*) hdNdxiMCPrimRec[i]->Clone(strNameEffXi);
2599 hEffXi[i]->Divide(hdNdxiMCPrimRec[i],hdNdxiMCPrimGen[i],1,1,"B"); // binominal errors
2600
2601 hEffZ[i] = (TH1F*) hdNdzMCPrimRec[i]->Clone(strNameEffZ);
2602 hEffZ[i]->Divide(hdNdzMCPrimRec[i],hdNdzMCPrimGen[i],1,1,"B"); // binominal errors
2603 }
2604
2605 // write
2606
2607 TString outfileOption = "RECREATE";
2608 if(updateOutfile) outfileOption = "UPDATE";
2609
2610 TFile out(strOutfile,outfileOption);
2611
2612 if(!out.IsOpen()){
2613 Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2614 return;
2615 }
2616
2617 if(fDebug>0) Printf("%s:%d -- write efficiency to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2618
2619// if(strdir && strdir.Length()){
2620// TDirectory* dir = out.mkdir(strdir);
2621// dir->cd();
2622// }
2623
2624 for(Int_t i=0; i<fNJetPtSlices; i++){
2625
2626 hdNdptTracksMCPrimGen[i]->Write();
2627 hdNdxiMCPrimGen[i]->Write();
2628 hdNdzMCPrimGen[i]->Write();
2629
2630 hdNdptTracksMCPrimRec[i]->Write();
2631 hdNdxiMCPrimRec[i]->Write();
2632 hdNdzMCPrimRec[i]->Write();
2633
2634 hEffPt[i]->Write();
2635 hEffXi[i]->Write();
2636 hEffZ[i]->Write();
2637 }
2638
2639 out.Close();
2640
2641}
2642
2643//________________________________________________________________________________________________________________
2644void AliFragmentationFunctionCorrections::WriteJetSecCorr(TString strInfile, TString strID, TString strOutfile,
2645 Bool_t updateOutfile)
2646{
2647 // read task ouput from MC and write secondary correction - standard dir/list
2648
2649 TString strdir = "PWG4_FragmentationFunction_" + strID;
2650 TString strlist = "fracfunc_" + strID;
2651
2652 WriteJetSecCorr(strInfile,strdir,strlist,strOutfile,updateOutfile);
2653}
2654
2655//___________________________________________________________________________________________________________________________________
2656void AliFragmentationFunctionCorrections::WriteJetSecCorr(TString strInfile, TString strdir, TString strlist,
2657 TString strOutfile, Bool_t updateOutfile)
2658{
2659 // read task output from MC and write secondary correction in jet pt slices as function of pt, z, xi
2660 // argument strlist optional - read from directory strdir if not specified
2661 // write eff to file strOutfile - by default only 'update' file (don't overwrite)
2662
2663 TH1F* hSecCorrPt[fNJetPtSlices];
2664 TH1F* hSecCorrXi[fNJetPtSlices];
2665 TH1F* hSecCorrZ[fNJetPtSlices];
2666
2667 TH1F* hdNdptTracksMCPrimRec[fNJetPtSlices];
2668 TH1F* hdNdxiMCPrimRec[fNJetPtSlices];
2669 TH1F* hdNdzMCPrimRec[fNJetPtSlices];
2670
2671 TH1F* hdNdptTracksMCSecRec[fNJetPtSlices];
2672 TH1F* hdNdxiMCSecRec[fNJetPtSlices];
2673 TH1F* hdNdzMCSecRec[fNJetPtSlices];
2674
2675 TH1F* fh1FFJetPtRecEffGen;
2676
2677 TH2F* fh2FFTrackPtRecEffRec;
2678 TH2F* fh2FFZRecEffRec;
2679 TH2F* fh2FFXiRecEffRec;
2680
2681 TH2F* fh2FFTrackPtSecRec;
2682 TH2F* fh2FFZSecRec;
2683 TH2F* fh2FFXiSecRec;
2684
2685 TFile f(strInfile,"READ");
2686
2687 if(!f.IsOpen()){
2688 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2689 return;
2690 }
2691
2692 if(fDebug>0) Printf("%s:%d -- writeJetTrackEff: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
2693
2694 if(strdir && strdir.Length()) gDirectory->cd(strdir);
2695
2696 TList* list = 0;
2697
2698 if(strlist && strlist.Length()){
2699
2700 if(!(list = (TList*) gDirectory->Get(strlist))){
2701 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2702 return;
2703 }
2704 }
2705
2706 if(list){
2707 fh1FFJetPtRecEffGen = (TH1F*) list->FindObject("fh1FFJetPtRecEffGen");
2708
2709 fh2FFTrackPtRecEffRec = (TH2F*) list->FindObject("fh2FFTrackPtRecEffRec");
2710 fh2FFZRecEffRec = (TH2F*) list->FindObject("fh2FFZRecEffRec");
2711 fh2FFXiRecEffRec = (TH2F*) list->FindObject("fh2FFXiRecEffRec");
2712
2713 fh2FFTrackPtSecRec = (TH2F*) list->FindObject("fh2FFTrackPtSecRec");
2714 fh2FFZSecRec = (TH2F*) list->FindObject("fh2FFZSecRec");
2715 fh2FFXiSecRec = (TH2F*) list->FindObject("fh2FFXiSecRec");
2716 }
2717 else{
2718 fh1FFJetPtRecEffGen = (TH1F*) gDirectory->FindObject("fh1FFJetPtRecEffGen");
2719
2720 fh2FFTrackPtRecEffRec = (TH2F*) gDirectory->FindObject("fh2FFTrackPtRecEffRec");
2721 fh2FFZRecEffRec = (TH2F*) gDirectory->FindObject("fh2FFZRecEffRec");
2722 fh2FFXiRecEffRec = (TH2F*) gDirectory->FindObject("fh2FFXiRecEffRec");
2723
2724 fh2FFTrackPtSecRec = (TH2F*) gDirectory->FindObject("fh2FFTrackPtSecRec");
2725 fh2FFZSecRec = (TH2F*) gDirectory->FindObject("fh2FFZSecRec");
2726 fh2FFXiSecRec = (TH2F*) gDirectory->FindObject("fh2FFXiSecRec");
2727 }
2728
2729 fh1FFJetPtRecEffGen->SetDirectory(0);
2730
2731 fh2FFTrackPtRecEffRec->SetDirectory(0);
2732 fh2FFZRecEffRec->SetDirectory(0);
2733 fh2FFXiRecEffRec->SetDirectory(0);
2734
2735 fh2FFTrackPtSecRec->SetDirectory(0);
2736 fh2FFZSecRec->SetDirectory(0);
2737 fh2FFXiSecRec->SetDirectory(0);
2738
2739 f.Close();
2740
2741
2742 // projections: FF for generated and reconstructed primaries
2743
2744 for(Int_t i=0; i<fNJetPtSlices; i++){
2745
2746 Float_t jetPtLoLim = fJetPtSlices->At(i);
2747 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
2748
2749 Int_t binLo = static_cast<Int_t>(fh2FFTrackPtRecEffRec->GetXaxis()->FindBin(jetPtLoLim));
2750 Int_t binUp = static_cast<Int_t>(fh2FFTrackPtRecEffRec->GetXaxis()->FindBin(jetPtUpLim))-1;
2751
2752 if(binUp > fh2FFTrackPtRecEffRec->GetNbinsX()){
2753 Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
2754 return;
2755 }
2756
2757 TString strNameFFPtPrimRec(Form("fh1FFTrackPtRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2758 TString strNameFFZPrimRec(Form("fh1FFZRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2759 TString strNameFFXiPrimRec(Form("fh1FFXiRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2760
2761 TString strNameFFPtSecRec(Form("fh1FFTrackPtRecSec_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2762 TString strNameFFZSecRec(Form("fh1FFZRecSec_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2763 TString strNameFFXiSecRec(Form("fh1FFXiRecSec_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2764
2765 // project
2766 // appendix 'unbinned' to avoid histos with same name after rebinning
2767
2768 hdNdptTracksMCPrimRec[i] = (TH1F*) fh2FFTrackPtRecEffRec->ProjectionY(strNameFFPtPrimRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
2769 hdNdzMCPrimRec[i] = (TH1F*) fh2FFZRecEffRec->ProjectionY(strNameFFZPrimRec+"_unBinned",binLo,binUp,"o");
2770 hdNdxiMCPrimRec[i] = (TH1F*) fh2FFXiRecEffRec->ProjectionY(strNameFFXiPrimRec+"_unBinned",binLo,binUp,"o");
2771
2772 hdNdptTracksMCSecRec[i] = (TH1F*) fh2FFTrackPtSecRec->ProjectionY(strNameFFPtSecRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
2773 hdNdzMCSecRec[i] = (TH1F*) fh2FFZSecRec->ProjectionY(strNameFFZSecRec+"_unBinned",binLo,binUp,"o");
2774 hdNdxiMCSecRec[i] = (TH1F*) fh2FFXiSecRec->ProjectionY(strNameFFXiSecRec+"_unBinned",binLo,binUp,"o");
2775
2776 // rebin
2777
2778 if(fNHistoBinsPt[i]) hdNdptTracksMCPrimRec[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtPrimRec,fHistoBinsPt[i]->GetArray());
2779 if(fNHistoBinsZ[i]) hdNdzMCPrimRec[i] = (TH1F*) hdNdzMCPrimRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZPrimRec,fHistoBinsZ[i]->GetArray());
2780 if(fNHistoBinsXi[i]) hdNdxiMCPrimRec[i] = (TH1F*) hdNdxiMCPrimRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiPrimRec,fHistoBinsXi[i]->GetArray());
2781
2782 if(fNHistoBinsPt[i]) hdNdptTracksMCSecRec[i] = (TH1F*) hdNdptTracksMCSecRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtSecRec,fHistoBinsPt[i]->GetArray());
2783 if(fNHistoBinsZ[i]) hdNdzMCSecRec[i] = (TH1F*) hdNdzMCSecRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZSecRec,fHistoBinsZ[i]->GetArray());
2784 if(fNHistoBinsXi[i]) hdNdxiMCSecRec[i] = (TH1F*) hdNdxiMCSecRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiSecRec,fHistoBinsXi[i]->GetArray());
2785
2786 hdNdptTracksMCPrimRec[i]->SetNameTitle(strNameFFPtPrimRec,"");
2787 hdNdzMCPrimRec[i]->SetNameTitle(strNameFFZPrimRec,"");
2788 hdNdxiMCPrimRec[i]->SetNameTitle(strNameFFXiPrimRec,"");
2789
2790 hdNdptTracksMCSecRec[i]->SetNameTitle(strNameFFPtSecRec,"");
2791 hdNdzMCSecRec[i]->SetNameTitle(strNameFFZSecRec,"");
2792 hdNdxiMCSecRec[i]->SetNameTitle(strNameFFXiSecRec,"");
2793
2794 // normalize
2795
2796 Double_t nJetsBin = fh1FFJetPtRecEffGen->Integral(binLo,binUp);
2797
2798 NormalizeTH1(hdNdptTracksMCPrimRec[i],nJetsBin);
2799 NormalizeTH1(hdNdzMCPrimRec[i],nJetsBin);
2800 NormalizeTH1(hdNdxiMCPrimRec[i],nJetsBin);
2801
2802 NormalizeTH1(hdNdptTracksMCSecRec[i],nJetsBin);
2803 NormalizeTH1(hdNdzMCSecRec[i],nJetsBin);
2804 NormalizeTH1(hdNdxiMCSecRec[i],nJetsBin);
2805
2806 // divide rec/gen : efficiency
2807
2808 TString strNameSecCorrPt(Form("hSecCorrPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2809 TString strNameSecCorrZ(Form("hSecCorrZ_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2810 TString strNameSecCorrXi(Form("hSecCorrXi_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2811
2812 hSecCorrPt[i] = (TH1F*) hdNdptTracksMCSecRec[i]->Clone(strNameSecCorrPt);
2813 TH1F* hSumPrimSecPt = (TH1F*) hdNdptTracksMCSecRec[i]->Clone("hSumPrimSecPt");
2814 hSumPrimSecPt->Add(hdNdptTracksMCPrimRec[i]);
2815 hSecCorrPt[i]->Divide(hdNdptTracksMCPrimRec[i],hSumPrimSecPt,1,1,"B"); // binominal errors
2816
2817 hSecCorrXi[i] = (TH1F*) hdNdxiMCSecRec[i]->Clone(strNameSecCorrXi);
2818 TH1F* hSumPrimSecXi = (TH1F*) hdNdxiMCSecRec[i]->Clone("hSumPrimSecXi");
2819 hSumPrimSecXi->Add(hdNdxiMCPrimRec[i]);
2820 hSecCorrXi[i]->Divide(hdNdxiMCPrimRec[i],hSumPrimSecXi,1,1,"B"); // binominal errors
2821
2822 hSecCorrZ[i] = (TH1F*) hdNdzMCSecRec[i]->Clone(strNameSecCorrZ);
2823 TH1F* hSumPrimSecZ = (TH1F*) hdNdzMCSecRec[i]->Clone("hSumPrimSecZ");
2824 hSumPrimSecZ->Add(hdNdzMCPrimRec[i]);
2825 hSecCorrZ[i]->Divide(hdNdzMCPrimRec[i],hSumPrimSecZ,1,1,"B"); // binominal errors
2826 }
2827
2828 // write
2829
2830 TString outfileOption = "RECREATE";
2831 if(updateOutfile) outfileOption = "UPDATE";
2832
2833 TFile out(strOutfile,outfileOption);
2834
2835 if(!out.IsOpen()){
2836 Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2837 return;
2838 }
2839
2840 if(fDebug>0) Printf("%s:%d -- write efficiency to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2841
2842 for(Int_t i=0; i<fNJetPtSlices; i++){
2843
2844 // hdNdptTracksMCSecRec[i]->Write();
2845 // hdNdxiMCSecRec[i]->Write();
2846 // hdNdzMCSecRec[i]->Write();
2847
2848 hSecCorrPt[i]->Write();
2849 hSecCorrXi[i]->Write();
2850 hSecCorrZ[i]->Write();
2851 }
2852
2853 out.Close();
2854}
2855
2856//________________________________________________________________________________________________________________
2857void AliFragmentationFunctionCorrections::WriteJetResponse(TString strInfile, TString strID, TString strOutfile,
2858 Bool_t updateOutfile)
2859{
2860 // read task ouput from MC and write single track eff - standard dir/list
2861
2862 TString strdir = "PWG4_FragmentationFunction_" + strID;
2863 TString strlist = "fracfunc_" + strID;
2864
2865 WriteJetResponse(strInfile,strdir,strlist,strOutfile,updateOutfile);
2866}
2867
2868//_____________________________________________________________________________________________________________________________________
2869void AliFragmentationFunctionCorrections::WriteJetResponse(TString strInfile, TString strdir, TString strlist,
2870 TString strOutfile, Bool_t updateOutfile)
2871{
2872 // read 3d THnSparse response matrices in pt,z,xi vs jet pt from file
2873 // project THnSparse + TH2 in jet pt slices
2874 // write to strOutfile
2875
2876 THnSparse* hn3ResponseJetPt;
2877 THnSparse* hn3ResponseJetZ;
2878 THnSparse* hn3ResponseJetXi;
2879
2880 // 2D response matrices
2881
2882 THnSparse* hnResponsePt[fNJetPtSlices];
2883 THnSparse* hnResponseZ[fNJetPtSlices];
2884 THnSparse* hnResponseXi[fNJetPtSlices];
2885
2886 TH2F* h2ResponsePt[fNJetPtSlices];
2887 TH2F* h2ResponseZ[fNJetPtSlices];
2888 TH2F* h2ResponseXi[fNJetPtSlices];
2889
2890 // 1D projections on gen pt / rec pt axes
2891
2892 TH1F* h1FFPtRec[fNJetPtSlices];
2893 TH1F* h1FFZRec[fNJetPtSlices];
2894 TH1F* h1FFXiRec[fNJetPtSlices];
2895
2896 TH1F* h1FFPtGen[fNJetPtSlices];
2897 TH1F* h1FFZGen[fNJetPtSlices];
2898 TH1F* h1FFXiGen[fNJetPtSlices];
2899
2900 TH1F* h1RatioPt[fNJetPtSlices];
2901 TH1F* h1RatioZ[fNJetPtSlices];
2902 TH1F* h1RatioXi[fNJetPtSlices];
2903
2904
2905
2906 TFile f(strInfile,"READ");
2907
2908 if(!f.IsOpen()){
2909 Printf("%s:%d -- error opening file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2910 return;
2911 }
2912
2913 if(fDebug>0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
2914
2915 if(strdir && strdir.Length()) gDirectory->cd(strdir);
2916
2917 TList* list = 0;
2918
2919 if(strlist && strlist.Length()){
2920
2921 if(!(list = (TList*) gDirectory->Get(strlist))){
2922 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2923 return;
2924 }
2925 }
2926
2927 if(list){
2928 hn3ResponseJetPt = (THnSparse*) list->FindObject("fhnResponseJetTrackPt");
2929 hn3ResponseJetZ = (THnSparse*) list->FindObject("fhnResponseJetZ");
2930 hn3ResponseJetXi = (THnSparse*) list->FindObject("fhnResponseJetXi");
2931 }
2932 else{
2933 hn3ResponseJetPt = (THnSparse*) gDirectory->Get("fhnResponseJetTrackPt");
2934 hn3ResponseJetZ = (THnSparse*) gDirectory->Get("fhnResponseJetZ");
2935 hn3ResponseJetXi = (THnSparse*) gDirectory->Get("fhnResponseJetXi");
2936 }
2937
2938
2939 if(!hn3ResponseJetPt){
2940 Printf("%s:%d -- error retrieving response matrix fhnResponseJetTrackPt",(char*)__FILE__,__LINE__);
2941 return;
2942 }
2943
2944 if(!hn3ResponseJetZ){
2945 Printf("%s:%d -- error retrieving response matrix fhnResponseJetZ",(char*)__FILE__,__LINE__);
2946 return;
2947 }
2948
2949 if(!hn3ResponseJetXi){
2950 Printf("%s:%d -- error retrieving response matrix fhnResponseJetXi",(char*)__FILE__,__LINE__);
2951 return;
2952 }
2953
2954 f.Close();
2955
2956 // axes
2957
2958 Int_t axisJetPtTHn3 = -1;
2959 Int_t axisGenPtTHn3 = -1;
2960 Int_t axisRecPtTHn3 = -1;
2961
2962 for(Int_t i=0; i<hn3ResponseJetPt->GetNdimensions(); i++){
2963
2964 TString title = hn3ResponseJetPt->GetAxis(i)->GetTitle();
2965
2966 if(title.Contains("jet p_{T}")) axisJetPtTHn3 = i;
2967 if(title.Contains("gen p_{T}")) axisGenPtTHn3 = i;
2968 if(title.Contains("rec p_{T}")) axisRecPtTHn3 = i;
2969 }
2970
2971 if(axisJetPtTHn3 == -1){
2972 Printf("%s:%d -- error axisJetPtTHn3",(char*)__FILE__,__LINE__);
2973 return;
2974 }
2975
2976 if(axisGenPtTHn3 == -1){
2977 Printf("%s:%d -- error axisGenPtTHn3",(char*)__FILE__,__LINE__);
2978 return;
2979 }
2980
2981 if(axisRecPtTHn3 == -1){
2982 Printf("%s:%d -- error axisRecPtTHn3",(char*)__FILE__,__LINE__);
2983 return;
2984 }
2985
2986 for(Int_t i=0; i<fNJetPtSlices; i++){
2987
2988 Float_t jetPtLoLim = fJetPtSlices->At(i);
2989 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
2990
2991 Int_t binLo = static_cast<Int_t>(hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->FindBin(jetPtLoLim));
2992 Int_t binUp = static_cast<Int_t>(hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->FindBin(jetPtUpLim))-1;
2993
2994 if(binUp > hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->GetNbins()){
2995 Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
2996 return;
2997 }
2998
2999 TString strNameRespPt(Form("hnResponsePt_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3000 TString strNameRespZ(Form("hnResponseZ_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3001 TString strNameRespXi(Form("hnResponseXi_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3002
3003 TString strNameTH2RespPt(Form("h2ResponsePt_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3004 TString strNameTH2RespZ(Form("h2ResponseZ_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3005 TString strNameTH2RespXi(Form("h2ResponseXi_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3006
3007 TString strNameRecPt(Form("h1FFTrackPtRecPrim_recPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3008 TString strNameRecZ(Form("h1FFZRecPrim_recPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3009 TString strNameRecXi(Form("h1FFXiRecPrim_recPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3010
3011 TString strNameGenPt(Form("h1FFTrackPtRecPrim_genPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3012 TString strNameGenZ(Form("h1FFZRecPrim_genPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3013 TString strNameGenXi(Form("h1FFXiRecPrim_genPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3014
3015 TString strNameRatioPt(Form("h1RatioTrackPtRecPrim_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3016 TString strNameRatioZ(Form("h1RatioZRecPrim_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3017 TString strNameRatioXi(Form("h1RatioXiRecPrim_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3018
3019
3020 // 2D projections in jet pt range
3021
3022 hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->SetRange(binLo,binUp);
3023 hn3ResponseJetZ->GetAxis(axisJetPtTHn3)->SetRange(binLo,binUp);
3024 hn3ResponseJetXi->GetAxis(axisJetPtTHn3)->SetRange(binLo,binUp);
3025
3026 Int_t axesProj[2] = {axisRecPtTHn3,axisGenPtTHn3};
3027
3028 hnResponsePt[i] = hn3ResponseJetPt->Projection(2,axesProj);
3029 hnResponseZ[i] = hn3ResponseJetZ->Projection(2,axesProj);
3030 hnResponseXi[i] = hn3ResponseJetXi->Projection(2,axesProj);
3031
3032 hnResponsePt[i]->SetNameTitle(strNameRespPt,"");
3033 hnResponseZ[i]->SetNameTitle(strNameRespZ,"");
3034 hnResponseXi[i]->SetNameTitle(strNameRespXi,"");
3035
3036 h2ResponsePt[i] = (TH2F*) hnResponsePt[i]->Projection(1,0);// note convention: yDim,xDim
3037 h2ResponseZ[i] = (TH2F*) hnResponseZ[i]->Projection(1,0); // note convention: yDim,xDim
3038 h2ResponseXi[i] = (TH2F*) hnResponseXi[i]->Projection(1,0);// note convention: yDim,xDim
3039
3040 h2ResponsePt[i]->SetNameTitle(strNameTH2RespPt,"");
3041 h2ResponseZ[i]->SetNameTitle(strNameTH2RespZ,"");
3042 h2ResponseXi[i]->SetNameTitle(strNameTH2RespXi,"");
3043
3044
3045 // 1D projections
3046
3047 Int_t axisGenPtTHn2 = -1;
3048 Int_t axisRecPtTHn2 = -1;
3049
3050 for(Int_t d=0; d<hnResponsePt[i]->GetNdimensions(); d++){
3051
3052 TString title = hnResponsePt[i]->GetAxis(d)->GetTitle();
3053
3054 if(title.Contains("gen p_{T}")) axisGenPtTHn2 = d;
3055 if(title.Contains("rec p_{T}")) axisRecPtTHn2 = d;
3056 }
3057
3058
3059 if(axisGenPtTHn2 == -1){
3060 Printf("%s:%d -- error axisGenPtTHn2",(char*)__FILE__,__LINE__);
3061 return;
3062 }
3063
3064 if(axisRecPtTHn2 == -1){
3065 Printf("%s:%d -- error axisRecPtTHn2",(char*)__FILE__,__LINE__);
3066 return;
3067 }
3068
3069
3070 h1FFPtRec[i] = (TH1F*) hnResponsePt[i]->Projection(axisRecPtTHn2);// note convention: yDim,xDim
3071 h1FFZRec[i] = (TH1F*) hnResponseZ[i]->Projection(axisRecPtTHn2);// note convention: yDim,xDim
3072 h1FFXiRec[i] = (TH1F*) hnResponseXi[i]->Projection(axisRecPtTHn2);// note convention: yDim,xDim
3073
3074 h1FFPtRec[i]->SetNameTitle(strNameRecPt,"");
3075 h1FFZRec[i]->SetNameTitle(strNameRecZ,"");
3076 h1FFXiRec[i]->SetNameTitle(strNameRecXi,"");
3077
3078 h1FFPtGen[i] = (TH1F*) hnResponsePt[i]->Projection(axisGenPtTHn2);// note convention: yDim,xDim
3079 h1FFZGen[i] = (TH1F*) hnResponseZ[i]->Projection(axisGenPtTHn2);// note convention: yDim,xDim
3080 h1FFXiGen[i] = (TH1F*) hnResponseXi[i]->Projection(axisGenPtTHn2);// note convention: yDim,xDim
3081
3082 h1FFPtGen[i]->SetNameTitle(strNameGenPt,"");
3083 h1FFZGen[i]->SetNameTitle(strNameGenZ,"");
3084 h1FFXiGen[i]->SetNameTitle(strNameGenXi,"");
3085
3086 // normalize 1D projections
3087
3088 if(fNHistoBinsPt[i]) h1FFPtRec[i] = (TH1F*) h1FFPtRec[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray());
3089 if(fNHistoBinsZ[i]) h1FFZRec[i] = (TH1F*) h1FFZRec[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3090 if(fNHistoBinsXi[i]) h1FFXiRec[i] = (TH1F*) h1FFXiRec[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3091
3092 if(fNHistoBinsPt[i]) h1FFPtGen[i] = (TH1F*) h1FFPtGen[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray());
3093 if(fNHistoBinsZ[i]) h1FFZGen[i] = (TH1F*) h1FFZGen[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3094 if(fNHistoBinsXi[i]) h1FFXiGen[i] = (TH1F*) h1FFXiGen[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3095
3096 NormalizeTH1(h1FFPtRec[i],fNJets->At(i));
3097 NormalizeTH1(h1FFZRec[i],fNJets->At(i));
3098 NormalizeTH1(h1FFXiRec[i],fNJets->At(i));
3099
3100 NormalizeTH1(h1FFPtGen[i],fNJets->At(i));
3101 NormalizeTH1(h1FFZGen[i],fNJets->At(i));
3102 NormalizeTH1(h1FFXiGen[i],fNJets->At(i));
3103
3104 // ratios 1D projections
3105
3106 h1RatioPt[i] = (TH1F*) h1FFPtRec[i]->Clone(strNameRatioPt);
3107 h1RatioPt[i]->Reset();
3108 h1RatioPt[i]->Divide(h1FFPtRec[i],h1FFPtGen[i],1,1,"B");
3109
3110 h1RatioZ[i] = (TH1F*) h1FFZRec[i]->Clone(strNameRatioZ);
3111 h1RatioZ[i]->Reset();
3112 h1RatioZ[i]->Divide(h1FFZRec[i],h1FFZGen[i],1,1,"B");
3113
3114 h1RatioXi[i] = (TH1F*) h1FFXiRec[i]->Clone(strNameRatioXi);
3115 h1RatioXi[i]->Reset();
3116 h1RatioXi[i]->Divide(h1FFXiRec[i],h1FFXiGen[i],1,1,"B");
3117 }
3118
3119
3120 // write
3121
3122 TString outfileOption = "RECREATE";
3123 if(updateOutfile) outfileOption = "UPDATE";
3124
3125 TFile out(strOutfile,outfileOption);
3126
3127 if(!out.IsOpen()){
3128 Printf("%s:%d -- error opening response matrix output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
3129 return;
3130 }
3131
3132 if(fDebug>0) Printf("%s:%d -- write response matrices to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
3133
3134 //if(strdir && strdir.Length()){
3135 // TDirectory* dir = out.mkdir(strdir);
3136 // dir->cd();
3137 //}
3138
3139 for(Int_t i=0; i<fNJetPtSlices; i++){
3140
3141 hnResponsePt[i]->Write();
3142 hnResponseXi[i]->Write();
3143 hnResponseZ[i]->Write();
3144
3145 h2ResponsePt[i]->Write();
3146 h2ResponseXi[i]->Write();
3147 h2ResponseZ[i]->Write();
3148
3149 h1FFPtRec[i]->Write();
3150 h1FFZRec[i]->Write();
3151 h1FFXiRec[i]->Write();
3152
3153 h1FFPtGen[i]->Write();
3154 h1FFZGen[i]->Write();
3155 h1FFXiGen[i]->Write();
3156
3157 h1RatioPt[i]->Write();
3158 h1RatioZ[i]->Write();
3159 h1RatioXi[i]->Write();
3160
3161 }
3162
3163 out.Close();
3164}
3165
3166//______________________________________________________________________________________________________
3167void AliFragmentationFunctionCorrections::ReadResponse(TString strfile, TString strdir, TString strlist)
3168{
3169 // read response matrices from file
3170 // argument strlist optional - read from directory strdir if not specified
3171 // note: THnSparse are not rebinned
3172
3173 THnSparse* hRespPt[fNJetPtSlices];
3174 THnSparse* hRespZ[fNJetPtSlices];
3175 THnSparse* hRespXi[fNJetPtSlices];
3176
3177 for(Int_t i=0; i<fNJetPtSlices; i++) hRespPt[i] = 0;
3178 for(Int_t i=0; i<fNJetPtSlices; i++) hRespZ[i] = 0;
3179 for(Int_t i=0; i<fNJetPtSlices; i++) hRespXi[i] = 0;
3180
3181 TFile f(strfile,"READ");
3182
3183 if(!f.IsOpen()){
3184 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
3185 return;
3186 }
3187
3188 if(fDebug>0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3189
3190 if(strdir && strdir.Length()) gDirectory->cd(strdir);
3191
3192 TList* list = 0;
3193
3194 if(strlist && strlist.Length()){
3195
3196 if(!(list = (TList*) gDirectory->Get(strlist))){
3197 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
3198 return;
3199 }
3200 }
3201
3202 for(Int_t i=0; i<fNJetPtSlices; i++){
3203
3204 Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
3205 Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
3206
3207 TString strNameRespPt(Form("hnResponsePt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3208 TString strNameRespZ(Form("hnResponseZ_%02d_%02d",jetPtLoLim,jetPtUpLim));
3209 TString strNameRespXi(Form("hnResponseXi_%02d_%02d",jetPtLoLim,jetPtUpLim));
3210
3211 if(list){
3212 hRespPt[i] = (THnSparse*) list->FindObject(strNameRespPt);
3213 hRespZ[i] = (THnSparse*) list->FindObject(strNameRespZ);
3214 hRespXi[i] = (THnSparse*) list->FindObject(strNameRespXi);
3215 }
3216 else{
3217 hRespPt[i] = (THnSparse*) gDirectory->Get(strNameRespPt);
3218 hRespZ[i] = (THnSparse*) gDirectory->Get(strNameRespZ);
3219 hRespXi[i] = (THnSparse*) gDirectory->Get(strNameRespXi);
3220 }
3221
3222 if(!hRespPt[i]){
3223 Printf("%s:%d -- error retrieving response %s", (char*)__FILE__,__LINE__,strNameRespPt.Data());
3224 }
3225
3226 if(!hRespZ[i]){
3227 Printf("%s:%d -- error retrieving response %s", (char*)__FILE__,__LINE__,strNameRespZ.Data());
3228 }
3229
3230 if(!hRespXi[i]){
3231 Printf("%s:%d -- error retrieving response %s", (char*)__FILE__,__LINE__,strNameRespXi.Data());
3232 }
3233
3234 // if(0){ // can't rebin THnSparse ...
3235 // if(fNHistoBinsPt[i]) hRespPt[i]->SetBinEdges(0,fHistoBinsPt[i]->GetArray());
3236 // if(fNHistoBinsPt[i]) hRespPt[i]->SetBinEdges(1,fHistoBinsPt[i]->GetArray());
3237
3238 // if(fNHistoBinsZ[i]) hRespZ[i]->SetBinEdges(0,fHistoBinsZ[i]->GetArray());
3239 // if(fNHistoBinsZ[i]) hRespZ[i]->SetBinEdges(1,fHistoBinsZ[i]->GetArray());
3240
3241 // if(fNHistoBinsXi[i]) hRespXi[i]->SetBinEdges(0,fHistoBinsXi[i]->GetArray());
3242 // if(fNHistoBinsXi[i]) hRespXi[i]->SetBinEdges(1,fHistoBinsXi[i]->GetArray());
3243 // }
3244
3245
3246 } // jet slices loop
3247
3248 f.Close(); // THnSparse pointers still valid even if file closed
3249
3250// for(Int_t i=0; i<fNJetPtSlices; i++){ // no copy c'tor ...
3251// if(hRespPt[i]) new(fhnResponsePt[i]) THnSparseF(*hRespPt[i]);
3252// if(hRespZ[i]) new(fhnResponseZ[i]) THnSparseF(*hRespZ[i]);
3253// if(hRespXi[i]) new(fhnResponseXi[i]) THnSparseF(*hRespXi[i]);
3254// }
3255
3256 for(Int_t i=0; i<fNJetPtSlices; i++){
3257 fhnResponsePt[i] = hRespPt[i];
3258 fhnResponseZ[i] = hRespZ[i];
3259 fhnResponseXi[i] = hRespXi[i];
3260 }
3261}
3262
3263//______________________________________________________________________________________________________________________
3264void AliFragmentationFunctionCorrections::ReadPriors(TString strfile,const Int_t type)
3265{
3266 // read priors from file: rec primaries, gen pt dist
3267
3268 if(fDebug>0) Printf("%s:%d -- read priors from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3269
3270 // temporary histos to store pointers from file
3271 TH1F* hist[fNJetPtSlices];
3272
3273 for(Int_t i=0; i<fNJetPtSlices; i++) hist[i] = 0;
3274
3275 TFile f(strfile,"READ");
3276
3277 if(!f.IsOpen()){
3278 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
3279 return;
3280 }
3281
3282 for(Int_t i=0; i<fNJetPtSlices; i++){
3283
3284 Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
3285 Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
3286
3287 TString strName;
3288
3289 if(type == kFlagPt) strName.Form("h1FFTrackPtRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim);
3290 if(type == kFlagZ) strName.Form("h1FFZRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim);
3291 if(type == kFlagXi) strName.Form("h1FFXiRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim);
3292
3293 hist[i] = (TH1F*) gDirectory->Get(strName);
3294
3295 if(!hist[i]){
3296 Printf("%s:%d -- error retrieving prior %s", (char*)__FILE__,__LINE__,strName.Data());
3297 }
3298
3299
3300 //if(fNHistoBinsPt[i]) hist[i] = (TH1F*) hist[i]->Rebin(fNHistoBinsPt[i],hist[i]->GetName()+"_rebin",fHistoBinsPt[i]->GetArray());
3301
3302 if(hist[i]) hist[i]->SetDirectory(0);
3303
3304 } // jet slices loop
3305
3306 f.Close();
3307
3308
3309 for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
3310 if(hist[i] && type == kFlagPt) new(fh1FFTrackPtPrior[i]) TH1F(*hist[i]);
3311 if(hist[i] && type == kFlagZ) new(fh1FFZPrior[i]) TH1F(*hist[i]);
3312 if(hist[i] && type == kFlagXi) new(fh1FFXiPrior[i]) TH1F(*hist[i]);
3313 }
3314}
3315
3316//_____________________________________________________
3317// void AliFragmentationFunctionCorrections::RatioRecGen()
3318// {
3319// // create ratio reconstructed over generated FF
3320// // use current highest corrLevel
3321
3322// Printf("%s:%d -- build ratio rec.gen, corrLevel %d",(char*)__FILE__,__LINE__,fNCorrectionLevels-1);
3323
3324// for(Int_t i=0; i<fNJetPtSlices; i++){
3325
3326// TH1F* histPtRec = fCorrFF[fNCorrectionLevels-1]->GetTrackPt(i); // levels -1: latest corr level
3327// TH1F* histZRec = fCorrFF[fNCorrectionLevels-1]->GetZ(i); // levels -1: latest corr level
3328// TH1F* histXiRec = fCorrFF[fNCorrectionLevels-1]->GetXi(i); // levels -1: latest corr level
3329
3330// TH1F* histPtGen = fh1FFTrackPtGenPrim[i];
3331// TH1F* histZGen = fh1FFZGenPrim[i];
3332// TH1F* histXiGen = fh1FFXiGenPrim[i];
3333
3334// TString histNamePt = histPtRec->GetName();
3335// TString histNameZ = histZRec->GetName();
3336// TString histNameXi = histXiRec->GetName();
3337
3338// histNamePt.ReplaceAll("fh1FF","fh1FFRatioRecGen");
3339// histNameZ.ReplaceAll("fh1FF","fh1FFRatioRecGen");
3340// histNameXi.ReplaceAll("fh1FF","fh1FFRatioRecGen");
3341
3342// // ratio
3343// TH1F* hRatioRecGenPt = (TH1F*) histPtRec->Clone(histNamePt);
3344// hRatioRecGenPt->Reset();
3345// hRatioRecGenPt->Divide(histPtRec,histPtGen,1,1,"B");
3346
3347// TH1F* hRatioRecGenZ = (TH1F*) histZRec->Clone(histNameZ);
3348// hRatioRecGenZ->Reset();
3349// hRatioRecGenZ->Divide(histZRec,histZGen,1,1,"B");
3350
3351// TH1F* hRatioRecGenXi = (TH1F*) histXiRec->Clone(histNameXi);
3352// hRatioRecGenXi->Reset();
3353// hRatioRecGenXi->Divide(histXiRec,histXiGen,1,1,"B");
3354
3355// new(fh1FFRatioRecGenPt[i]) TH1F(*hRatioRecGenPt);
3356// new(fh1FFRatioRecGenZ[i]) TH1F(*hRatioRecGenZ);
3357// new(fh1FFRatioRecGenXi[i]) TH1F(*hRatioRecGenXi);
3358// }
3359// }
3360
3361// //___________________________________________________________
3362// void AliFragmentationFunctionCorrections::RatioRecPrimaries()
3363// {
3364// // create ratio reconstructed tracks over reconstructed primaries
3365// // use raw FF (corrLevel 0)
3366
3367// Printf("%s:%d -- build ratio rec tracks /rec primaries",(char*)__FILE__,__LINE__);
3368
3369// for(Int_t i=0; i<fNJetPtSlices; i++){
3370
3371// const Int_t corrLevel = 0;
3372
3373// TH1F* histPtRec = fCorrFF[corrLevel]->GetTrackPt(i); // levels -1: latest corr level
3374// TH1F* histZRec = fCorrFF[corrLevel]->GetZ(i); // levels -1: latest corr level
3375// TH1F* histXiRec = fCorrFF[corrLevel]->GetXi(i); // levels -1: latest corr level
3376
3377// TH1F* histPtRecPrim = fh1FFTrackPtRecPrim[i];
3378// TH1F* histZRecPrim = fh1FFZRecPrim[i];
3379// TH1F* histXiRecPrim = fh1FFXiRecPrim[i];
3380
3381// TString histNamePt = histPtRec->GetName();
3382// TString histNameZ = histZRec->GetName();
3383// TString histNameXi = histXiRec->GetName();
3384
3385// histNamePt.ReplaceAll("fh1FF","fh1FFRatioRecPrim");
3386// histNameZ.ReplaceAll("fh1FF","fh1FFRatioRecPrim");
3387// histNameXi.ReplaceAll("fh1FF","fh1FFRatioRecPrim");
3388
3389// // ratio
3390// TH1F* hRatioRecPrimPt = (TH1F*) histPtRec->Clone(histNamePt);
3391// hRatioRecPrimPt->Reset();
3392// hRatioRecPrimPt->Divide(histPtRec,histPtRecPrim,1,1,"B");
3393
3394// TH1F* hRatioRecPrimZ = (TH1F*) histZRec->Clone(histNameZ);
3395// hRatioRecPrimZ->Reset();
3396// hRatioRecPrimZ->Divide(histZRec,histZRecPrim,1,1,"B");
3397
3398// TH1F* hRatioRecPrimXi = (TH1F*) histXiRec->Clone(histNameXi);
3399// hRatioRecPrimXi->Reset();
3400// hRatioRecPrimXi->Divide(histXiRec,histXiRecPrim,1,1,"B");
3401
3402
3403// new(fh1FFRatioRecPrimPt[i]) TH1F(*hRatioRecPrimPt);
3404// new(fh1FFRatioRecPrimZ[i]) TH1F(*hRatioRecPrimZ);
3405// new(fh1FFRatioRecPrimXi[i]) TH1F(*hRatioRecPrimXi);
3406// }
3407// }
3408
3409// __________________________________________________________________________________
3410void AliFragmentationFunctionCorrections::ProjectJetResponseMatrices(TString strOutfile)
3411{
3412
3413 // project response matrices on both axes:
3414 // FF for rec primaries, in terms of generated and reconstructed momentum
3415 // write FF and ratios to outFile
3416
3417 Printf("%s:%d -- project response matrices, write to %s",(char*)__FILE__,__LINE__,strOutfile.Data());
3418
3419 TH1F* hFFPtRec[fNJetPtSlices];
3420 TH1F* hFFZRec[fNJetPtSlices];
3421 TH1F* hFFXiRec[fNJetPtSlices];
3422
3423 TH1F* hFFPtGen[fNJetPtSlices];
3424 TH1F* hFFZGen[fNJetPtSlices];
3425 TH1F* hFFXiGen[fNJetPtSlices];
3426
3427 TH1F* hRatioPt[fNJetPtSlices];
3428 TH1F* hRatioZ[fNJetPtSlices];
3429 TH1F* hRatioXi[fNJetPtSlices];
3430
3431
3432 Int_t axisGenPt = 1;
3433 Int_t axisRecPt = 0;
3434
3435 for(Int_t i=0; i<fNJetPtSlices; i++){
3436
3437 Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
3438 Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
3439
3440 TString strNameRecPt(Form("h1FFTrackPtRecPrim_recPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3441 TString strNameRecZ(Form("h1FFZRecPrim_recPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3442 TString strNameRecXi(Form("h1FFXiRecPrim_recPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3443
3444 TString strNameGenPt(Form("h1FFTrackPtRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3445 TString strNameGenZ(Form("h1FFZRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3446 TString strNameGenXi(Form("h1FFXiRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3447
3448 TString strNameRatioPt(Form("h1RatioTrackPtRecPrim_%02d_%02d",jetPtLoLim,jetPtUpLim));
3449 TString strNameRatioZ(Form("h1RatioZRecPrim_%02d_%02d",jetPtLoLim,jetPtUpLim));
3450 TString strNameRatioXi(Form("h1RatioXiRecPrim_%02d_%02d",jetPtLoLim,jetPtUpLim));
3451
3452
3453 hFFPtRec[i] = (TH1F*) fhnResponsePt[i]->Projection(axisRecPt);// note convention: yDim,xDim
3454 hFFZRec[i] = (TH1F*) fhnResponseZ[i]->Projection(axisRecPt);// note convention: yDim,xDim
3455 hFFXiRec[i] = (TH1F*) fhnResponseXi[i]->Projection(axisRecPt);// note convention: yDim,xDim
3456
3457 hFFPtRec[i]->SetNameTitle(strNameRecPt,"");
3458 hFFZRec[i]->SetNameTitle(strNameRecZ,"");
3459 hFFXiRec[i]->SetNameTitle(strNameRecXi,"");
3460
3461
3462 hFFPtGen[i] = (TH1F*) fhnResponsePt[i]->Projection(axisGenPt);// note convention: yDim,xDim
3463 hFFZGen[i] = (TH1F*) fhnResponseZ[i]->Projection(axisGenPt);// note convention: yDim,xDim
3464 hFFXiGen[i] = (TH1F*) fhnResponseXi[i]->Projection(axisGenPt);// note convention: yDim,xDim
3465
3466 hFFPtGen[i]->SetNameTitle(strNameGenPt,"");
3467 hFFZGen[i]->SetNameTitle(strNameGenZ,"");
3468 hFFXiGen[i]->SetNameTitle(strNameGenXi,"");
3469
3470
3471 if(fNHistoBinsPt[i]) hFFPtRec[i] = (TH1F*) hFFPtRec[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray());
3472 if(fNHistoBinsZ[i]) hFFZRec[i] = (TH1F*) hFFZRec[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3473 if(fNHistoBinsXi[i]) hFFXiRec[i] = (TH1F*) hFFXiRec[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3474
3475 if(fNHistoBinsPt[i]) hFFPtGen[i] = (TH1F*) hFFPtGen[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray());
3476 if(fNHistoBinsZ[i]) hFFZGen[i] = (TH1F*) hFFZGen[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3477 if(fNHistoBinsXi[i]) hFFXiGen[i] = (TH1F*) hFFXiGen[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3478
3479 NormalizeTH1(hFFPtGen[i],fNJets->At(i));
3480 NormalizeTH1(hFFZGen[i],fNJets->At(i));
3481 NormalizeTH1(hFFXiGen[i],fNJets->At(i));
3482
3483 NormalizeTH1(hFFPtRec[i],fNJets->At(i));
3484 NormalizeTH1(hFFZRec[i],fNJets->At(i));
3485 NormalizeTH1(hFFXiRec[i],fNJets->At(i));
3486
3487
3488 hRatioPt[i] = (TH1F*) hFFPtRec[i]->Clone(strNameRatioPt);
3489 hRatioPt[i]->Reset();
3490 hRatioPt[i]->Divide(hFFPtRec[i],hFFPtGen[i],1,1,"B");
3491
3492 hRatioZ[i] = (TH1F*) hFFZRec[i]->Clone(strNameRatioZ);
3493 hRatioZ[i]->Reset();
3494 hRatioZ[i]->Divide(hFFZRec[i],hFFZGen[i],1,1,"B");
3495
3496 hRatioXi[i] = (TH1F*) hFFXiRec[i]->Clone(strNameRatioXi);
3497 hRatioXi[i]->Reset();
3498 hRatioXi[i]->Divide(hFFXiRec[i],hFFXiGen[i],1,1,"B");
3499 }
3500
3501
3502
3503 // write
3504
3505 TFile out(strOutfile,"RECREATE");
3506
3507 if(!out.IsOpen()){
3508 Printf("%s:%d -- error opening response matrix output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
3509 return;
3510 }
3511
3512 for(Int_t i=0; i<fNJetPtSlices; i++){
3513
3514 hFFPtRec[i]->Write();
3515 hFFZRec[i]->Write();
3516 hFFXiRec[i]->Write();
3517
3518 hFFPtGen[i]->Write();
3519 hFFZGen[i]->Write();
3520 hFFXiGen[i]->Write();
3521
3522 hRatioPt[i]->Write();
3523 hRatioZ[i]->Write();
3524 hRatioXi[i]->Write();
3525 }
3526
3527 out.Close();
3528}
3529
3530// ____________________________________________________________________________________________________________________________
3531void AliFragmentationFunctionCorrections::ProjectSingleResponseMatrix(TString strOutfile, Bool_t updateOutfile, TString strOutDir )
3532{
3533 // project response matrix on both axes:
3534 // pt spec for rec primaries, in terms of generated and reconstructed momentum
3535 // write spec and ratios to outFile
3536
3537 Printf("%s:%d -- project single pt response matrix, write to %s",(char*)__FILE__,__LINE__,strOutfile.Data());
3538
3539 TH1F* hSpecPtRec;
3540 TH1F* hSpecPtGen;
3541 TH1F* hRatioPt;
3542
3543 Int_t axisGenPt = 1;
3544 Int_t axisRecPt = 0;
3545
3546 TString strNameRecPt = "h1SpecTrackPtRecPrim_recPt";
3547 TString strNameGenPt = "h1SpecTrackPtRecPrim_genPt";
3548 TString strNameRatioPt = "h1RatioTrackPtRecPrim";
3549
3550 hSpecPtRec = (TH1F*) fhnResponseSinglePt->Projection(axisRecPt);// note convention: yDim,xDim
3551 hSpecPtRec->SetNameTitle(strNameRecPt,"");
3552
3553 hSpecPtGen = (TH1F*) fhnResponseSinglePt->Projection(axisGenPt);// note convention: yDim,xDim
3554 hSpecPtGen->SetNameTitle(strNameGenPt,"");
3555
3556 hRatioPt = (TH1F*) hSpecPtRec->Clone(strNameRatioPt);
3557 hRatioPt->Reset();
3558 hRatioPt->Divide(hSpecPtRec,hSpecPtGen,1,1,"B");
3559
3560 TString outfileOption = "RECREATE";
3561 if(updateOutfile) outfileOption = "UPDATE";
3562
3563 TFile out(strOutfile,outfileOption);
3564
3565 if(!out.IsOpen()){
3566 Printf("%s:%d -- error opening reponse matrix projections output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
3567 return;
3568 }
3569
3570
3571 if(strOutDir && strOutDir.Length()){
3572
3573 TDirectory* dir;
3574 if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd();
3575 else{
3576 dir = out.mkdir(strOutDir);
3577 dir->cd();
3578 }
3579 }
3580
3581 hSpecPtRec->Write();
3582 hSpecPtGen->Write();
3583 hRatioPt->Write();
3584
3585 out.Close();
3586}
3587
3588
3589//__________________________________________________________________________________________________________________________________________________________________
3590void AliFragmentationFunctionCorrections::RebinHisto(const Int_t jetPtSlice, const Int_t nBinsLimits, Double_t* binsLimits, Double_t* binsWidth, const Int_t type)
3591{
3592 // rebin histo, rescale bins according to new width
3593 // only correct for input histos with equal bin size
3594
3595 // args: jetPtSlice, type, use current corr level
3596
3597 // function takes array of limits and widths (e.g. 1 GeV bins up to 10 GeV, 2 GeV width up to 20 GeV, ...)
3598 // array size of binsLimits: nBinsLimits
3599 // array size of binsWidth: nBinsLimits-1
3600 // binsLimits have to be in increasing order
3601 // if binning undefined for any slice, original binning will be kept
3602
3603 if(!fNJetPtSlices){
3604 Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
3605 return;
3606 }
3607
3608 if(jetPtSlice>=fNJetPtSlices){
3609 Printf("%s:%d -- jetPtSlice %d exceeds max",(char*)__FILE__,__LINE__,jetPtSlice);
3610 return;
3611 }
3612
3613
3614 Double_t binLimitMin = binsLimits[0];
3615 Double_t binLimitMax = binsLimits[nBinsLimits-1];
3616
3617 Double_t binLimit = binLimitMin; // start value
3618
3619 Int_t sizeUpperLim = 1000; //static_cast<Int_t>(binLimitMax/binsWidth[0])+1; - only if first bin has smallest width, but not case for dN/dxi ...
3620 TArrayD binsArray(sizeUpperLim);
3621 Int_t nBins = 0;
3622 binsArray.SetAt(binLimitMin,nBins++);
3623
3624 while(binLimit<binLimitMax && nBins<sizeUpperLim){
3625
3626 Int_t currentSlice = -1;
3627 for(Int_t i=0; i<nBinsLimits; i++){
3628 if(binLimit >= 0.999*binsLimits[i]) currentSlice = i; // 0.999 numerical saftey factor
3629 }
3630
3631 Double_t currentBinWidth = binsWidth[currentSlice];
3632 binLimit += currentBinWidth;
3633
3634 binsArray.SetAt(binLimit,nBins++);
3635 }
3636
3637
3638 TH1F* hist = 0;
a5592cfa 3639 if(type == kFlagPt) hist = fCorrFF[fNCorrectionLevels-1]->GetTrackPt(jetPtSlice);
1aa4f09f 3640 else if(type == kFlagZ) hist = fCorrFF[fNCorrectionLevels-1]->GetZ(jetPtSlice);
3641 else if(type == kFlagXi) hist = fCorrFF[fNCorrectionLevels-1]->GetXi(jetPtSlice);
3642 else if(type == kFlagSinglePt) hist = fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->GetTrackPt(0);
3643 else{
3644 Printf("%s%d unknown type",(char*)__FILE__,__LINE__);
3645 return;
3646 }
39e2b057 3647
3648 Double_t binWidthNoRebin = hist->GetBinWidth(1);
3649
3650 Double_t* bins = binsArray.GetArray();
3651
3652 hist = (TH1F*) hist->Rebin(nBins-1,"",bins);
3653
3654 for(Int_t bin=0; bin <= hist->GetNbinsX(); bin++){
3655
3656 Double_t binWidthRebin = hist->GetBinWidth(bin);
3657 Double_t scaleF = binWidthNoRebin / binWidthRebin;
3658
3659 Double_t binCont = hist->GetBinContent(bin);
3660 Double_t binErr = hist->GetBinError(bin);
3661
3662 binCont *= scaleF;
3663 binErr *= scaleF;
3664
3665 hist->SetBinContent(bin,binCont);
3666 hist->SetBinError(bin,binErr);
3667 }
3668
3669
3670
3671 TH1F* temp = new TH1F(*hist);
3672
3673 if(type == kFlagPt) fCorrFF[fNCorrectionLevels-1]->ReplaceCorrHistos(jetPtSlice,temp,0,0);
3674 if(type == kFlagZ) fCorrFF[fNCorrectionLevels-1]->ReplaceCorrHistos(jetPtSlice,0,temp,0);
3675 if(type == kFlagXi) fCorrFF[fNCorrectionLevels-1]->ReplaceCorrHistos(jetPtSlice,0,0,temp);
3676 if(type == kFlagSinglePt) fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->ReplaceCorrHistos(0,temp,0,0);
3677
3678
3679 delete temp;
3680}
3681//__________________________________________________________________________________________________________________________________________________________________
633ed180 3682void AliFragmentationFunctionCorrections::WriteJetSpecResponse(TString strInfile, TString strdir, TString strlist/*, TString strOutfile*/)
39e2b057 3683{
3684
3685 if(fDebug>0) Printf("%s:%d -- read jet spectrum response matrix from file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
3686
3687 if(strdir && strdir.Length()) gDirectory->cd(strdir);
3688
3689 TList* list = 0;
3690
3691 if(strlist && strlist.Length()){
3692
3693 if(!(list = (TList*) gDirectory->Get(strlist))){
3694 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
3695 return;
3696 }
3697 }
3698
a5592cfa 3699 THnSparse* hn6ResponseJetPt = (THnSparse*) list->FindObject("fhnCorrelation");
3700
39e2b057 3701 Int_t axis6RecJetPt = 0;
3702 Int_t axis6GenJetPt = 3;
3703
3704 hn6ResponseJetPt->GetAxis(axis6RecJetPt)->SetTitle("rec jet p_{T} (GeV/c)");
3705 hn6ResponseJetPt->GetAxis(axis6GenJetPt)->SetTitle("gen jet p_{T} (GeV/c)");
3706
3707 Int_t nBinsRecPt = hn6ResponseJetPt->GetAxis(axis6RecJetPt)->GetNbins();
3708 Double_t loLimRecPt = hn6ResponseJetPt->GetAxis(axis6RecJetPt)->GetBinLowEdge(1);
3709 Double_t upLimRecPt = hn6ResponseJetPt->GetAxis(axis6RecJetPt)->GetBinUpEdge(nBinsRecPt);
3710
3711 Int_t nBinsGenPt = hn6ResponseJetPt->GetAxis(axis6GenJetPt)->GetNbins();
3712 Double_t loLimGenPt = hn6ResponseJetPt->GetAxis(axis6GenJetPt)->GetBinLowEdge(1);
3713 Double_t upLimGenPt = hn6ResponseJetPt->GetAxis(axis6GenJetPt)->GetBinUpEdge(nBinsGenPt);
3714
3715 Int_t nBinsTrackPt = 200;
3716 Int_t loLimTrackPt = 0;
3717 Int_t upLimTrackPt = 200;
3718
3719
3720 Int_t nBinsResponse[4] = {nBinsRecPt,nBinsTrackPt,nBinsGenPt,nBinsTrackPt};
3721 Double_t binMinResponse[4] = {loLimRecPt,loLimTrackPt,loLimGenPt,loLimTrackPt};
3722 Double_t binMaxResponse[4] = {upLimRecPt,upLimTrackPt,upLimGenPt,upLimTrackPt};
3723
3724 const char* labelsResponseSinglePt[4] = {"rec jet p_{T} (GeV/c)", "rec track p_{T} (GeV/c)", "gen jet p_{T} (GeV/c)", "gen track p_{T} (GeV/c)"};
3725
3726 THnSparseD* hn4ResponseTrackPtJetPt = new THnSparseD("hn4ResponseTrackPtJetPt","",4,nBinsResponse,binMinResponse,binMaxResponse);
3727
3728 for(Int_t i=0; i<4; i++){
3729 hn4ResponseTrackPtJetPt->GetAxis(i)->SetTitle(labelsResponseSinglePt[i]);
3730 }
3731
3732
3733 // fill
3734
3735
3736}
3737
3738//_____________________________________________________________________________________________________________________________________
3739void AliFragmentationFunctionCorrections::ReadSingleTrackEfficiency(TString strfile, TString strdir, TString strlist, TString strname)
3740{
3741
3742 ReadSingleTrackCorrection(strfile,strdir,strlist,strname,kFlagEfficiency);
3743
3744}
3745
3746//_____________________________________________________________________________________________________________________________________
3747void AliFragmentationFunctionCorrections::ReadSingleTrackResponse(TString strfile, TString strdir, TString strlist, TString strname)
3748{
3749
3750 ReadSingleTrackCorrection(strfile,strdir,strlist,strname,kFlagResponse);
3751
3752}
3753
3754//_____________________________________________________________________________________________________________________________________
3755void AliFragmentationFunctionCorrections::ReadSingleTrackSecCorr(TString strfile, TString strdir, TString strlist, TString strname)
3756{
3757
3758 ReadSingleTrackCorrection(strfile,strdir,strlist,strname,kFlagSecondaries);
3759
3760}
3761
3762//______________________________________________________________________________________________________________________________________________________
3763void AliFragmentationFunctionCorrections::ReadSingleTrackCorrection(TString strfile, TString strdir, TString strlist, TString strname, const Int_t type)
3764{
3765 // read single track correction (pt) from file
3766 // type: efficiency / response / secondaries correction
3767
3768 if(!((type == kFlagEfficiency) || (type == kFlagResponse) || (type == kFlagSecondaries))){
3769 Printf("%s:%d -- no such correction ",(char*)__FILE__,__LINE__);
3770 return;
3771 }
3772
3773 TFile f(strfile,"READ");
3774
3775 if(!f.IsOpen()){
3776 Printf("%s:%d -- error opening file %s", (char*)__FILE__,__LINE__,strfile.Data());
3777 return;
3778 }
3779
3780 if(fDebug>0 && type==kFlagEfficiency) Printf("%s:%d -- read single track corr from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3781 if(fDebug>0 && type==kFlagResponse) Printf("%s:%d -- read single track response from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3782 if(fDebug>0 && type==kFlagSecondaries) Printf("%s:%d -- read single track secondaries corr from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3783
3784 if(strdir && strdir.Length()) gDirectory->cd(strdir);
3785
3786 TList* list = 0;
3787
3788 if(strlist && strlist.Length()){
3789
3790 if(!(list = (TList*) gDirectory->Get(strlist))){
3791 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
3792 return;
3793 }
3794 }
3795
3796 TH1F* h1CorrHist = 0; // common TObject pointer not possible, need SetDirectory() later
3797 THnSparse* hnCorrHist = 0;
3798
3799 if(type == kFlagEfficiency || type == kFlagSecondaries){
3800
3801 if(list) h1CorrHist = (TH1F*) list->FindObject(strname);
3802 else h1CorrHist = (TH1F*) gDirectory->Get(strname);
3803
3804 if(!h1CorrHist){
3805 Printf("%s:%d -- error retrieving histo %s", (char*)__FILE__,__LINE__,strname.Data());
3806 return;
3807 }
3808
3809 }
3810 else if(type == kFlagResponse){
3811
3812 if(list) hnCorrHist = (THnSparse*) list->FindObject(strname);
3813 else hnCorrHist = (THnSparse*) gDirectory->Get(strname);
3814
3815 if(!hnCorrHist){
3816 Printf("%s:%d -- error retrieving histo %s", (char*)__FILE__,__LINE__,strname.Data());
3817 return;
3818 }
3819
3820 }
3821
3822 if(h1CorrHist) h1CorrHist->SetDirectory(0);
3823 //if(hnCorrHist) hnCorrHist->SetDirectory(0);
3824
3825 f.Close();
3826
3827 if(type == kFlagEfficiency) fh1EffSinglePt = h1CorrHist;
3828 else if(type == kFlagResponse) fhnResponseSinglePt = hnCorrHist;
3829 else if(type == kFlagSecondaries) fh1SecCorrSinglePt = h1CorrHist;
3830
3831}
3832
3833//________________________________________________________________________________________________________________
3834void AliFragmentationFunctionCorrections::ReadRawPtSpec(TString strInfile, TString strID)
3835{
3836 // read track pt spec from task ouput - standard dir/list
3837
3838 TString strdir = "PWG4_FragmentationFunction_" + strID;
3839 TString strlist = "fracfunc_" + strID;
3840
3841 ReadRawPtSpec(strInfile,strdir,strlist);
3842}
3843
3844//_______________________________________________________________________________________________________
3845void AliFragmentationFunctionCorrections::ReadRawPtSpec(TString strfile, TString strdir, TString strlist)
3846{
3847 // get raw pt spectra from file
3848
3849
3850 // book histos
3851 fNCorrectionLevelsSinglePt = 0;
3852 fCorrSinglePt = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels];
3853 AddCorrectionLevelSinglePt(); // first 'correction' level = raw spectrum
3854
3855 // get raw pt spec from input file, normalize
3856
3857 TFile f(strfile,"READ");
3858
3859 if(!f.IsOpen()){
3860 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
3861 return;
3862 }
3863
3864 if(fDebug>0) Printf("%s:%d -- read raw spectra from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3865
3866 gDirectory->cd(strdir);
3867
3868 TList* list = 0;
3869
3870 if(!(list = (TList*) gDirectory->Get(strlist))){
3871 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
3872 return;
3873 }
3874
3875 TString hnameTrackPt("fh1TrackQAPtRecCuts");
3876 TString hnameEvtSel("fh1EvtSelection");
3877
3878 TH1F* fh1TrackPt = (TH1F*) list->FindObject(hnameTrackPt);
3879 TH1F* fh1EvtSel = (TH1F*) list->FindObject(hnameEvtSel);
3880
3881 if(!fh1TrackPt){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameTrackPt.Data()); return; }
3882 if(!fh1EvtSel) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameEvtSel.Data()); return; }
3883
3884
3885 //Float_t nEvents = fh1EvtSel->GetBinContent(fh1EvtSel->FindBin(0));
3886
3887 // evts after physics selection
3888 Float_t nEvents = fh1EvtSel->GetEntries()
3889 - fh1EvtSel->GetBinContent(fh1EvtSel->FindBin(1)) // evts rejected by trigger selection ( = PhysicsSelection
3890 - fh1EvtSel->GetBinContent(fh1EvtSel->FindBin(2)); // evts with wrong centrality class
3891
3892
3893 fh1TrackPt->SetDirectory(0);
3894
3895 f.Close();
3896
3897
3898 NormalizeTH1(fh1TrackPt,nEvents);
3899
3900 // raw FF = corr level 0
3901 fCorrSinglePt[0]->AddCorrHistos(0,fh1TrackPt);
3902}
3903
3904
3905//_______________________________________________________________________________________________________
3906void AliFragmentationFunctionCorrections::ReadRawPtSpecQATask(TString strfile, TString strdir, TString strlist)
3907{
3908 // get raw pt spectra from file
3909 // for output from Martas QA task
3910
3911
3912 // book histos
3913 fNCorrectionLevelsSinglePt = 0;
3914 fCorrSinglePt = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels];
3915 AddCorrectionLevelSinglePt(); // first 'correction' level = raw spectrum
3916
3917 // get raw pt spec from input file, normalize
3918
3919 TFile f(strfile,"READ");
3920
3921 if(!f.IsOpen()){
3922 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
3923 return;
3924 }
3925
3926 if(fDebug>0) Printf("%s:%d -- read raw pt spec from QA task output file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3927
3928 gDirectory->cd(strdir);
3929
3930 TList* list = 0;
3931
3932 if(!(list = (TList*) gDirectory->Get(strlist))){
3933 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
3934 return;
3935 }
3936
3937 TString hnameTrackPt("fPtSel");
3938 TString hnameEvtSel("fNEventAll");
3939
3940 TH1F* fh1TrackPt = (TH1F*) list->FindObject(hnameTrackPt);
3941 TH1F* fh1EvtSel = (TH1F*) list->FindObject(hnameEvtSel);
3942
3943 if(!fh1TrackPt){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameTrackPt.Data()); return; }
3944 if(!fh1EvtSel) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameEvtSel.Data()); return; }
3945
3946
3947 // evts after physics selection
3948 Float_t nEvents = fh1EvtSel->GetEntries();
3949
3950 fh1TrackPt->SetDirectory(0);
3951
3952 f.Close();
3953
3954
3955 NormalizeTH1(fh1TrackPt,nEvents);
3956
3957 // raw FF = corr level 0
3958 fCorrSinglePt[0]->AddCorrHistos(0,fh1TrackPt);
3959}
3960
3961// ________________________________________________________
3962void AliFragmentationFunctionCorrections::EffCorrSinglePt()
3963{
3964 // apply efficiency correction to inclusive track pt spec
3965
3966 AddCorrectionLevelSinglePt("eff");
3967
3968 TH1F* histPt = fCorrSinglePt[fNCorrectionLevelsSinglePt-2]->GetTrackPt(0);
3969
3970 if(histPt->GetNbinsX() != fh1EffSinglePt->GetNbinsX()) Printf("%s:%d: inconsistency pt spec and eff corr bins - rebin effCorr ...", (char*)__FILE__,__LINE__);
3971
3972 TString histNamePt = histPt->GetName();
3973 TH1F* hTrackPtEffCorr = (TH1F*) histPt->Clone(histNamePt);
3974 hTrackPtEffCorr->Divide(histPt,fh1EffSinglePt,1,1,"");
3975
3976 fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->AddCorrHistos(0,hTrackPtEffCorr);
3977}
3978
3979//___________________________________________________________________________________________________________________________
3980void AliFragmentationFunctionCorrections::UnfoldSinglePt(const Int_t nIter, const Bool_t useCorrelatedErrors)
3981{
3982 // unfolde inclusive dN/dpt spectra
3983
3984 AddCorrectionLevelSinglePt("unfold");
3985
3986 TH1F* hist = fCorrSinglePt[fNCorrectionLevelsSinglePt-2]->GetTrackPt(0); // level -2: before unfolding, level -1: unfolded
3987 THnSparse* hnResponse = fhnResponseSinglePt;
3988
3989 TString histNameTHn = hist->GetName();
3990 if(histNameTHn.Contains("TH1")) histNameTHn.ReplaceAll("TH1","THn");
3991 if(histNameTHn.Contains("fPt")) histNameTHn.ReplaceAll("fPt","fhnPt");
3992
3993
3994 TString histNameBackFolded = hist->GetName();
3995 histNameBackFolded.Append("_backfold");
3996
3997 TString histNameRatioFolded = hist->GetName();
3998 if(histNameRatioFolded.Contains("fh1")) histNameRatioFolded.ReplaceAll("fh1","hRatio");
3999 if(histNameRatioFolded.Contains("fPt")) histNameRatioFolded.ReplaceAll("fPt","hRatioPt");
4000 histNameRatioFolded.Append("_unfold");
4001
4002 TString histNameRatioBackFolded = hist->GetName();
4003 if(histNameRatioBackFolded.Contains("fh1")) histNameRatioBackFolded.ReplaceAll("fh1","hRatio");
4004 if(histNameRatioBackFolded.Contains("fPt")) histNameRatioBackFolded.ReplaceAll("fPt","hRatioPt");
4005 histNameRatioBackFolded.Append("_backfold");
4006
4007 THnSparse* hnHist = TH1toSparse(hist,histNameTHn,hist->GetTitle());
4008 THnSparse* hnFlatEfficiency = TH1toSparse(hist,"fhnEfficiency","eff",kTRUE); // could optionally also use real eff
4009
4010 THnSparse* hnUnfolded
4011 = Unfold(hnHist,hnResponse,hnFlatEfficiency,nIter,useCorrelatedErrors);
4012
4013 TH1F* hUnfolded = (TH1F*) hnUnfolded->Projection(0);
4014 hUnfolded->SetNameTitle(hist->GetName(),hist->GetTitle());
4015
4016 fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->AddCorrHistos(0,hUnfolded);
4017
4018 // backfolding: apply response matrix to unfolded spectrum
4019 TH1F* hBackFolded = ApplyResponse(hUnfolded,hnResponse);
4020 hBackFolded->SetNameTitle(histNameBackFolded,hist->GetTitle());
4021
4022 fh1SingleTrackPtBackFolded = hBackFolded;
4023
4024
4025 // ratio unfolded to original histo
4026 TH1F* hRatioUnfolded = (TH1F*) hUnfolded->Clone(histNameRatioFolded);
4027 hRatioUnfolded->Reset();
4028 hRatioUnfolded->Divide(hUnfolded,hist,1,1,"B");
4029
4030 fh1RatioSingleTrackPtFolded = hRatioUnfolded;
4031
4032
4033 // ratio backfolded to original histo
4034 TH1F* hRatioBackFolded = (TH1F*) hBackFolded->Clone(histNameRatioBackFolded);
4035 hRatioBackFolded->Reset();
4036 hRatioBackFolded->Divide(hBackFolded,hist,1,1,"B");
4037
4038 fh1RatioSingleTrackPtBackFolded = hRatioBackFolded;
4039
4040 delete hnHist;
4041 delete hnFlatEfficiency;
4042
4043}
4044
4045// ________________________________________________________
4046void AliFragmentationFunctionCorrections::SecCorrSinglePt()
4047{
4048 // apply efficiency correction to inclusive track pt spec
4049
4050 AddCorrectionLevelSinglePt("secCorr");
4051
4052 TH1F* histPt = fCorrSinglePt[fNCorrectionLevelsSinglePt-2]->GetTrackPt(0);
4053
4054 if(histPt->GetNbinsX() != fh1SecCorrSinglePt->GetNbinsX())
4055 Printf("%s:%d: inconsistency pt spec and secondaries corr bins - rebin effCorr ...", (char*)__FILE__,__LINE__);
4056
4057 TString histNamePt = histPt->GetName();
4058 TH1F* hTrackPtSecCorr = (TH1F*) histPt->Clone(histNamePt);
4059
4060 hTrackPtSecCorr->Multiply(histPt,fh1SecCorrSinglePt,1,1,"");
4061
4062 fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->AddCorrHistos(0,hTrackPtSecCorr);
4063}