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