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