1 // *************************************************************************
3 // * corrections to Fragmentation Functions *
5 // *************************************************************************
8 /**************************************************************************
9 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
11 * Author: The ALICE Off-line Project. *
12 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
27 #include "THnSparse.h"
29 #include "TDirectory.h"
30 #include "AliCFUnfolding.h"
31 #include "AliFragmentationFunctionCorrections.h"
33 ///////////////////////////////////////////////////////////////////////////////
35 // correction class for ouput of AliAnalysisTaskFragmentationFunction //
37 // corrections for: reconstruction efficiency, momentum resolution, //
38 // secondaries, UE / HI background, fluctuations //
39 // back-correction on jet energy on dN/dxi //
41 // read MC ouput and write out efficiency histos, response matrices etc. //
42 // read measured distributions and apply efficiency, response matrices etc. //
44 // contact: o.busch@gsi.de //
46 ///////////////////////////////////////////////////////////////////////////////
49 ClassImp(AliFragmentationFunctionCorrections)
51 //________________________________________________________________________
52 AliFragmentationFunctionCorrections::AliFragmentationFunctionCorrections()
59 ,fNHistoBinsSinglePt(0)
60 ,fHistoBinsSinglePt(0)
67 ,fNCorrectionLevels(0)
69 ,fNCorrectionLevelsBgr(0)
71 ,fNCorrectionLevelsSinglePt(0)
81 ,fh1FFTrackPtBackFolded(0)
84 ,fh1FFRatioTrackPtFolded(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)
100 ,fh1SecCorrSinglePt(0)
102 // default constructor
105 //________________________________________________________________________________________________________________________
106 AliFragmentationFunctionCorrections::AliFragmentationFunctionCorrections(const AliFragmentationFunctionCorrections ©)
109 ,fNJetPtSlices(copy.fNJetPtSlices)
110 ,fJetPtSlices(copy.fJetPtSlices)
112 ,fNJetsBgr(copy.fNJetsBgr)
113 ,fNHistoBinsSinglePt(copy.fNHistoBinsSinglePt)
114 ,fHistoBinsSinglePt(copy.fHistoBinsSinglePt)
115 ,fNHistoBinsPt(copy.fNHistoBinsPt)
116 ,fNHistoBinsZ(copy.fNHistoBinsZ)
117 ,fNHistoBinsXi(copy.fNHistoBinsXi)
118 ,fHistoBinsPt(copy.fHistoBinsPt)
119 ,fHistoBinsZ(copy.fHistoBinsZ)
120 ,fHistoBinsXi(copy.fHistoBinsXi)
121 ,fNCorrectionLevels(copy.fNCorrectionLevels)
122 ,fCorrFF(copy.fCorrFF)
123 ,fNCorrectionLevelsBgr(copy.fNCorrectionLevelsBgr)
124 ,fCorrBgr(copy.fCorrBgr)
125 ,fNCorrectionLevelsSinglePt(copy.fNCorrectionLevelsSinglePt)
126 ,fCorrSinglePt(copy.fCorrSinglePt)
127 ,fh1FFXiShift(copy.fh1FFXiShift)
128 ,fh1EffSinglePt(copy.fh1EffSinglePt)
129 ,fh1EffPt(copy.fh1EffPt)
130 ,fh1EffZ(copy.fh1EffZ)
131 ,fh1EffXi(copy.fh1EffXi)
132 ,fh1EffBgrPt(copy.fh1EffBgrPt)
133 ,fh1EffBgrZ(copy.fh1EffBgrZ)
134 ,fh1EffBgrXi(copy.fh1EffBgrXi)
135 ,fh1FFTrackPtBackFolded(copy.fh1FFTrackPtBackFolded)
136 ,fh1FFZBackFolded(copy.fh1FFZBackFolded)
137 ,fh1FFXiBackFolded(copy.fh1FFXiBackFolded)
138 ,fh1FFRatioTrackPtFolded(copy.fh1FFRatioTrackPtFolded)
139 ,fh1FFRatioZFolded(copy.fh1FFRatioZFolded)
140 ,fh1FFRatioXiFolded(copy.fh1FFRatioXiFolded)
141 ,fh1FFRatioTrackPtBackFolded(copy.fh1FFRatioTrackPtBackFolded)
142 ,fh1FFRatioZBackFolded(copy.fh1FFRatioZBackFolded)
143 ,fh1FFRatioXiBackFolded(copy.fh1FFRatioXiBackFolded)
144 ,fh1SingleTrackPtBackFolded(copy.fh1SingleTrackPtBackFolded)
145 ,fh1RatioSingleTrackPtFolded(copy.fh1RatioSingleTrackPtFolded)
146 ,fh1RatioSingleTrackPtBackFolded(copy.fh1RatioSingleTrackPtBackFolded)
147 ,fhnResponseSinglePt(copy.fhnResponseSinglePt)
148 ,fhnResponsePt(copy.fhnResponsePt)
149 ,fhnResponseZ(copy.fhnResponseZ)
150 ,fhnResponseXi(copy.fhnResponseXi)
151 ,fh1FFTrackPtPrior(copy.fh1FFTrackPtPrior)
152 ,fh1FFZPrior(copy.fh1FFZPrior)
153 ,fh1FFXiPrior(copy.fh1FFXiPrior)
154 ,fh1SecCorrSinglePt(copy.fh1SecCorrSinglePt)
160 // ______________________________________________________________________________________________________________________________
161 AliFragmentationFunctionCorrections& AliFragmentationFunctionCorrections::operator=(const AliFragmentationFunctionCorrections& o)
166 TObject::operator=(o);
168 fNJetPtSlices = o.fNJetPtSlices;
169 fJetPtSlices = o.fJetPtSlices;
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;
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;
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;
219 //_________________________________________________________________________
220 AliFragmentationFunctionCorrections::~AliFragmentationFunctionCorrections()
224 if(fJetPtSlices) delete fJetPtSlices;
225 if(fNJets) delete fNJets;
226 if(fNJetsBgr) delete fNJetsBgr;
228 DeleteHistoArray(fh1FFXiShift);
230 DeleteHistoArray(fh1EffPt);
231 DeleteHistoArray(fh1EffXi);
232 DeleteHistoArray(fh1EffZ );
234 DeleteHistoArray(fh1EffBgrPt);
235 DeleteHistoArray(fh1EffBgrXi);
236 DeleteHistoArray(fh1EffBgrZ);
240 DeleteHistoArray(fh1FFTrackPtBackFolded);
241 DeleteHistoArray(fh1FFZBackFolded);
242 DeleteHistoArray(fh1FFXiBackFolded);
244 DeleteHistoArray(fh1FFRatioTrackPtFolded);
245 DeleteHistoArray(fh1FFRatioZFolded);
246 DeleteHistoArray(fh1FFRatioXiFolded);
248 DeleteHistoArray(fh1FFRatioTrackPtBackFolded);
249 DeleteHistoArray(fh1FFRatioZBackFolded);
250 DeleteHistoArray(fh1FFRatioXiBackFolded);
252 DeleteTHnSparseArray(fhnResponsePt);
253 DeleteTHnSparseArray(fhnResponseZ);
254 DeleteTHnSparseArray(fhnResponseXi);
256 DeleteHistoArray(fh1FFTrackPtPrior);
257 DeleteHistoArray(fh1FFZPrior);
258 DeleteHistoArray(fh1FFXiPrior);
261 // clean up corrected FF
263 for(Int_t c=0; c<fNCorrectionLevels; c++) delete fCorrFF[c];
268 for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) delete fCorrBgr[c];
271 // clean up inclusive pt
272 for(Int_t c=0; c<fNCorrectionLevelsSinglePt; c++) delete fCorrSinglePt[c];
273 delete[] fCorrSinglePt;
275 delete[] fNHistoBinsPt;
276 delete[] fNHistoBinsZ;
277 delete[] fNHistoBinsXi;
279 // clean up histo bins
281 if(fHistoBinsSinglePt) delete fHistoBinsSinglePt;
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];
287 delete[] fHistoBinsPt;
288 delete[] fHistoBinsZ;
289 delete[] fHistoBinsXi;
292 //_________________________________________________________________________________
293 AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::AliFragFuncCorrHistos()
301 // default constructor
305 //__________________________________________________________________________________________________________________
306 AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::AliFragFuncCorrHistos(const AliFragmentationFunctionCorrections::AliFragFuncCorrHistos& copy)
308 ,fArraySize(copy.fArraySize)
312 ,fCorrLabel(copy.fCorrLabel)
316 fh1CorrFFTrackPt = new TH1F*[copy.fArraySize];
317 fh1CorrFFZ = new TH1F*[copy.fArraySize];
318 fh1CorrFFXi = new TH1F*[copy.fArraySize];
320 for(Int_t i=0; i<copy.fArraySize; i++){
321 fh1CorrFFTrackPt[i] = new TH1F(*(copy.fh1CorrFFTrackPt[i]));
322 fh1CorrFFZ[i] = new TH1F(*(copy.fh1CorrFFZ[i]));
323 fh1CorrFFXi[i] = new TH1F(*(copy.fh1CorrFFXi[i]));
327 //_______________________________________________________________________________________________________________________________________________________________
328 AliFragmentationFunctionCorrections::AliFragFuncCorrHistos& AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::operator=(const AliFragmentationFunctionCorrections::AliFragFuncCorrHistos& o)
333 TObject::operator=(o);
334 Int_t fArraySize_new = o.fArraySize;
335 fCorrLabel = o.fCorrLabel;
337 TH1F** fh1CorrFFTrackPt_new = new TH1F*[fArraySize_new];
338 TH1F** fh1CorrFFZ_new = new TH1F*[fArraySize_new];
339 TH1F** fh1CorrFFXi_new = new TH1F*[fArraySize_new];
341 for(Int_t i=0; i<fArraySize_new; i++){
342 fh1CorrFFTrackPt_new[i] = new TH1F(*(o.fh1CorrFFTrackPt[i]));
343 fh1CorrFFZ_new[i] = new TH1F(*(o.fh1CorrFFZ[i]));
344 fh1CorrFFXi_new[i] = new TH1F(*(o.fh1CorrFFXi[i]));
350 for(int i=0; i<fArraySize; i++) delete fh1CorrFFTrackPt[i];
351 for(int i=0; i<fArraySize; i++) delete fh1CorrFFZ[i];
352 for(int i=0; i<fArraySize; i++) delete fh1CorrFFXi[i];
355 if(fh1CorrFFTrackPt) delete[] fh1CorrFFTrackPt;
356 if(fh1CorrFFZ) delete[] fh1CorrFFZ;
357 if(fh1CorrFFXi) delete[] fh1CorrFFXi;
361 fArraySize = fArraySize_new;
363 fh1CorrFFTrackPt = fh1CorrFFTrackPt_new;
364 fh1CorrFFZ = fh1CorrFFZ_new;
365 fh1CorrFFXi = fh1CorrFFXi_new;
367 for(Int_t i=0; i<fArraySize; i++){
368 fh1CorrFFTrackPt[i] = fh1CorrFFTrackPt_new[i];
369 fh1CorrFFZ[i] = fh1CorrFFZ_new[i];
370 fh1CorrFFXi[i] = fh1CorrFFXi_new[i];
377 //__________________________________________________________________________________
378 AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::~AliFragFuncCorrHistos()
383 for(int i=0; i<fArraySize; i++) delete fh1CorrFFTrackPt[i];
384 for(int i=0; i<fArraySize; i++) delete fh1CorrFFZ[i];
385 for(int i=0; i<fArraySize; i++) delete fh1CorrFFXi[i];
388 if(fh1CorrFFTrackPt) delete[] fh1CorrFFTrackPt;
389 if(fh1CorrFFZ) delete[] fh1CorrFFZ;
390 if(fh1CorrFFXi) delete[] fh1CorrFFXi;
394 //___________________________________________________________________________________________________________________
395 AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::AliFragFuncCorrHistos(const char* label, Int_t arraySize)
405 fArraySize = arraySize;
406 fh1CorrFFTrackPt = new TH1F*[arraySize];
407 fh1CorrFFZ = new TH1F*[arraySize];
408 fh1CorrFFXi = new TH1F*[arraySize];
410 for(int i=0; i<arraySize; i++) fh1CorrFFTrackPt[i] = new TH1F;
411 for(int i=0; i<arraySize; i++) fh1CorrFFZ[i] = new TH1F;
412 for(int i=0; i<arraySize; i++) fh1CorrFFXi[i] = new TH1F;
415 //_______________________________________________________________________________________________________________________________
416 void AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::AddCorrHistos(Int_t slice,TH1F* histPt,TH1F* histZ,TH1F* histXi)
418 // place histo in array, add corrLabel to histo name
420 if(slice>=fArraySize){
421 Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
426 TString name = histPt->GetName();
427 if(fCorrLabel.Length()>0) name += "_"+fCorrLabel;
428 histPt->SetName(name);
430 if(!(histPt->GetSumw2())) histPt->Sumw2();
432 new(fh1CorrFFTrackPt[slice]) TH1F(*histPt);
436 TString name = histZ->GetName();
437 if(fCorrLabel.Length()>0) name += "_"+fCorrLabel;
438 histZ->SetName(name);
440 if(!(histZ->GetSumw2())) histZ->Sumw2();
442 new(fh1CorrFFZ[slice]) TH1F(*histZ);
446 TString name = histXi->GetName();
447 if(fCorrLabel.Length()>0) name += "_"+fCorrLabel;
448 histXi->SetName(name);
450 if(!(histXi->GetSumw2())) histXi->Sumw2();
452 new(fh1CorrFFXi[slice]) TH1F(*histXi);
456 //___________________________________________________________________________________________________________________________________
457 void AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::ReplaceCorrHistos(Int_t slice,TH1F* histPt,TH1F* histZ,TH1F* histXi)
459 // replace histo in array
461 if(slice>=fArraySize){
462 Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
467 if(!(histPt->GetSumw2())) histPt->Sumw2();
469 TString name = histPt->GetName();
470 histPt->SetName(name);
472 delete fh1CorrFFTrackPt[slice];
473 fh1CorrFFTrackPt[slice] = new TH1F;
474 new(fh1CorrFFTrackPt[slice]) TH1F(*histPt);
478 if(!(histZ->GetSumw2())) histZ->Sumw2();
480 TString name = histZ->GetName();
481 histZ->SetName(name);
483 delete fh1CorrFFZ[slice];
484 fh1CorrFFZ[slice] = new TH1F;
485 new(fh1CorrFFZ[slice]) TH1F(*histZ);
489 if(!(histXi->GetSumw2())) histXi->Sumw2();
491 TString name = histXi->GetName();
492 histXi->SetName(name);
494 delete fh1CorrFFXi[slice];
495 fh1CorrFFXi[slice] = new TH1F;
496 new(fh1CorrFFXi[slice]) TH1F(*histXi);
500 // ___________________________________________________________________________________________
501 TH1F* AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::GetTrackPt(const Int_t slice)
505 if(slice>=fArraySize){
506 Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
510 return fh1CorrFFTrackPt[slice];
514 // ______________________________________________________________________________________
515 TH1F* AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::GetZ(const Int_t slice)
519 if(slice>=fArraySize){
520 Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
524 return fh1CorrFFZ[slice];
527 // ________________________________________________________________________________________
528 TH1F* AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::GetXi(const Int_t slice)
532 if(slice>=fArraySize){
533 Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
537 return fh1CorrFFXi[slice];
540 // __________________________________________________________________________
541 void AliFragmentationFunctionCorrections::DeleteHistoArray(TH1F** hist) const
543 // delete array of TH1
546 for(Int_t i=0; i<fNJetPtSlices; i++) delete hist[i];
551 // ____________________________________________________________________________________
552 void AliFragmentationFunctionCorrections::DeleteTHnSparseArray(THnSparse** hist) const
554 // delete array of THnSparse
557 for(Int_t i=0; i<fNJetPtSlices; i++) delete hist[i];
562 // _________________________________________________________
563 TH1F** AliFragmentationFunctionCorrections::BookHistoArray()
568 Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
572 TH1F** hist = new TH1F*[fNJetPtSlices];
573 for(Int_t i=0; i<fNJetPtSlices; i++) hist[i] = new TH1F();
578 //__________________________________________________________________
579 THnSparse** AliFragmentationFunctionCorrections::BookTHnSparseArray()
581 // book array of THnSparse
584 Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
588 THnSparse** hist = new THnSparse*[fNJetPtSlices];
589 for(Int_t i=0; i<fNJetPtSlices; i++) hist[i] = new THnSparseF();
594 //_____________________________________________________________________________
595 void AliFragmentationFunctionCorrections::AddCorrectionLevel(const char* label)
597 // increase corr level
600 Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
604 if(fNCorrectionLevels >= fgMaxNCorrectionLevels){
605 Printf("%s:%d -- max correction level exceeded", (char*)__FILE__,__LINE__);
609 fCorrFF[fNCorrectionLevels] = new AliFragFuncCorrHistos(label,fNJetPtSlices);
610 fNCorrectionLevels++;
614 //________________________________________________________________________________
615 void AliFragmentationFunctionCorrections::AddCorrectionLevelBgr(const char* label)
617 // increase corr level for bgr FF
620 if(fDebug>0) Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
624 if(fNCorrectionLevelsBgr >= fgMaxNCorrectionLevels){
625 Printf("%s:%d -- max correction level exceeded", (char*)__FILE__,__LINE__);
629 fCorrBgr[fNCorrectionLevelsBgr] = new AliFragFuncCorrHistos(label,fNJetPtSlices);
630 fNCorrectionLevelsBgr++;
633 //_____________________________________________________________________________________
634 void AliFragmentationFunctionCorrections::AddCorrectionLevelSinglePt(const char* label)
636 // increase corr level single pt spec
638 Int_t nJetPtSlicesSingle = 1; // pro forma
640 if(fNCorrectionLevelsSinglePt >= fgMaxNCorrectionLevels){
641 Printf("%s:%d -- max correction level exceeded", (char*)__FILE__,__LINE__);
645 fCorrSinglePt[fNCorrectionLevelsSinglePt] = new AliFragFuncCorrHistos(label,nJetPtSlicesSingle);
646 fNCorrectionLevelsSinglePt++;
649 //_____________________________________________________________________________________________
650 void AliFragmentationFunctionCorrections::SetJetPtSlices(Float_t* bins, const Int_t nJetPtSlices)
653 // once slices are known, can book all other histos
655 fNJetPtSlices = nJetPtSlices;
657 const Int_t size = nJetPtSlices+1;
658 fJetPtSlices = new TArrayF(size,bins);
662 fNJets = new TArrayF(size);
665 fNJetsBgr = new TArrayF(size);
670 fNCorrectionLevels = 0;
671 fCorrFF = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels];
672 AddCorrectionLevel(); // first 'correction' level = raw FF
674 fNCorrectionLevelsBgr = 0;
675 fCorrBgr = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels];
676 AddCorrectionLevelBgr(); // first 'correction' level = raw bgr dist
678 fh1FFXiShift = BookHistoArray();
682 fh1EffPt = BookHistoArray();
683 fh1EffXi = BookHistoArray();
684 fh1EffZ = BookHistoArray();
686 fh1EffBgrPt = BookHistoArray();
687 fh1EffBgrXi = BookHistoArray();
688 fh1EffBgrZ = BookHistoArray();
693 fh1FFTrackPtBackFolded = BookHistoArray();
694 fh1FFXiBackFolded = BookHistoArray();
695 fh1FFZBackFolded = BookHistoArray();
697 fh1FFRatioTrackPtFolded = BookHistoArray();
698 fh1FFRatioXiFolded = BookHistoArray();
699 fh1FFRatioZFolded = BookHistoArray();
701 fh1FFRatioTrackPtBackFolded = BookHistoArray();
702 fh1FFRatioXiBackFolded = BookHistoArray();
703 fh1FFRatioZBackFolded = BookHistoArray();
707 fhnResponsePt = BookTHnSparseArray();
708 fhnResponseZ = BookTHnSparseArray();
709 fhnResponseXi = BookTHnSparseArray();
711 fh1FFTrackPtPrior = BookHistoArray();
712 fh1FFZPrior = BookHistoArray();
713 fh1FFXiPrior = BookHistoArray();
718 fNHistoBinsPt = new Int_t[fNJetPtSlices];
719 fNHistoBinsXi = new Int_t[fNJetPtSlices];
720 fNHistoBinsZ = new Int_t[fNJetPtSlices];
722 for(Int_t i=0; i<fNJetPtSlices; i++) fNHistoBinsPt[i] = 0;
723 for(Int_t i=0; i<fNJetPtSlices; i++) fNHistoBinsXi[i] = 0;
724 for(Int_t i=0; i<fNJetPtSlices; i++) fNHistoBinsZ[i] = 0;
726 fHistoBinsPt = new TArrayD*[fNJetPtSlices];
727 fHistoBinsXi = new TArrayD*[fNJetPtSlices];
728 fHistoBinsZ = new TArrayD*[fNJetPtSlices];
730 for(Int_t i=0; i<fNJetPtSlices; i++) fHistoBinsPt[i] = new TArrayD(0);
731 for(Int_t i=0; i<fNJetPtSlices; i++) fHistoBinsXi[i] = new TArrayD(0);
732 for(Int_t i=0; i<fNJetPtSlices; i++) fHistoBinsZ[i] = new TArrayD(0);
735 //_____________________________________________________________________________________________________________________________________
736 void AliFragmentationFunctionCorrections::SetHistoBins(const Int_t jetPtSlice, const Int_t sizeBins, Double_t* bins, const Int_t type)
738 // set histo bins for jet pt slice
739 // if binning undefined for any slice, original binning will be used
742 Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
746 if(jetPtSlice>=fNJetPtSlices){
747 Printf("%s:%d -- jetPtSlice %d exceeds max",(char*)__FILE__,__LINE__,jetPtSlice);
752 fNHistoBinsPt[jetPtSlice] = sizeBins-1;
753 fHistoBinsPt[jetPtSlice]->Set(sizeBins,bins);
755 else if(type==kFlagZ){
756 fNHistoBinsZ[jetPtSlice] = sizeBins-1;
757 fHistoBinsZ[jetPtSlice]->Set(sizeBins,bins);
760 else if(type==kFlagXi){
761 fNHistoBinsXi[jetPtSlice] = sizeBins-1;
762 fHistoBinsXi[jetPtSlice]->Set(sizeBins,bins);
766 //__________________________________________________________________________________________________________________________________________________________________
767 void AliFragmentationFunctionCorrections::SetHistoBins(const Int_t jetPtSlice, const Int_t nBinsLimits, Double_t* binsLimits, Double_t* binsWidth, const Int_t type)
769 // set histo bins for jet pt slice
770 // function takes array of limits and widths (e.g. 1 GeV bins up to 10 GeV, 2 GeV width up to 20 GeV, ...)
771 // array size of binsLimits: nBinsLimits
772 // array size of binsWidth: nBinsLimits-1
773 // binsLimits have to be in increasing order
774 // if binning undefined for any slice, original binning will be used
777 Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
781 if(jetPtSlice>=fNJetPtSlices){
782 Printf("%s:%d -- jetPtSlice %d exceeds max",(char*)__FILE__,__LINE__,jetPtSlice);
787 Double_t binLimitMin = binsLimits[0];
788 Double_t binLimitMax = binsLimits[nBinsLimits-1];
790 Double_t binLimit = binLimitMin; // start value
792 Int_t sizeUpperLim = 10000; //static_cast<Int_t>(binLimitMax/binsWidth[0])+1;
793 TArrayD binsArray(sizeUpperLim);
795 binsArray.SetAt(binLimitMin,nBins++);
797 while(binLimit<binLimitMax && nBins<sizeUpperLim){
799 Int_t currentSlice = -1;
800 for(Int_t i=0; i<nBinsLimits; i++){
801 if(binLimit >= 0.999*binsLimits[i]) currentSlice = i; // 0.999 numerical saftey factor
804 Double_t currentBinWidth = binsWidth[currentSlice];
805 binLimit += currentBinWidth;
807 binsArray.SetAt(binLimit,nBins++);
810 Double_t* bins = binsArray.GetArray();
812 SetHistoBins(jetPtSlice,nBins,bins,type);
815 //__________________________________________________________________________________________________
816 void AliFragmentationFunctionCorrections::SetHistoBinsSinglePt(const Int_t sizeBins, Double_t* bins)
818 // set histo bins for inclusive pt spectra
819 // if binning undefined, original binning will be used
821 fNHistoBinsSinglePt = sizeBins-1;
823 fHistoBinsSinglePt = new TArrayD(sizeBins);
824 fHistoBinsSinglePt->Set(sizeBins,bins);
827 //__________________________________________________________________________________________________________________________________________________________________
828 void AliFragmentationFunctionCorrections::SetHistoBinsSinglePt(const Int_t nBinsLimits, Double_t* binsLimits, Double_t* binsWidth)
830 // set histo bins for inclusive pt spectra
831 // function takes array of limits and widths (e.g. 1 GeV bins up to 10 GeV, 2 GeV width up to 20 GeV, ...)
832 // array size of binsLimits: nBinsLimits
833 // array size of binsWidth: nBinsLimits-1
834 // binsLimits have to be in increasing order
835 // if binning undefined for any slice, original binning will be used
838 Double_t binLimitMin = binsLimits[0];
839 Double_t binLimitMax = binsLimits[nBinsLimits-1];
841 Double_t binLimit = binLimitMin; // start value
843 Int_t sizeUpperLim = 10000; //static_cast<Int_t>(binLimitMax/binsWidth[0])+1;
844 TArrayD binsArray(sizeUpperLim);
846 binsArray.SetAt(binLimitMin,nBins++);
848 while(binLimit<binLimitMax && nBins<sizeUpperLim){
850 Int_t currentSlice = -1;
851 for(Int_t i=0; i<nBinsLimits; i++){
852 if(binLimit >= 0.999*binsLimits[i]) currentSlice = i; // 0.999 numerical saftey factor
855 Double_t currentBinWidth = binsWidth[currentSlice];
856 binLimit += currentBinWidth;
858 binsArray.SetAt(binLimit,nBins++);
861 Double_t* bins = binsArray.GetArray();
863 SetHistoBinsSinglePt(nBins,bins);
866 //____________________________________________________________________________________
867 void AliFragmentationFunctionCorrections::NormalizeTH1(TH1* hist, const Float_t nJets)
869 // FF normalization: divide by bin width and normalize to entries/jets
870 // should also work for TH2, but to be tested !!!
872 if(nJets == 0){ // nothing to do
873 if(fDebug>0) Printf("%s:%d -- normalize: nJets = 0, do nothing", (char*)__FILE__,__LINE__);
877 Int_t nBins = hist->GetNbinsX()*hist->GetNbinsY()*hist->GetNbinsZ();
879 for(Int_t bin=0; bin<=nBins; bin++){ // include under-/overflow (?)
881 Double_t binWidth = hist->GetBinWidth(bin);
882 Double_t binCont = hist->GetBinContent(bin);
883 Double_t binErr = hist->GetBinError(bin);
888 hist->SetBinContent(bin,binCont);
889 hist->SetBinError(bin,binErr);
892 hist->Scale(1/nJets);
895 //_____________________________________________________
896 void AliFragmentationFunctionCorrections::NormalizeFF()
900 if(fNCorrectionLevels>1){
901 Printf("%s:%d -- FF normalization should be done BEFORE any correction !!!", (char*)__FILE__,__LINE__);
904 for(Int_t i=0; i<fNJetPtSlices; i++){
906 if(fDebug>0) Printf(" normalizeFF: i %d, nJets %f",i,fNJets->At(i));
908 NormalizeTH1(fCorrFF[0]->GetTrackPt(i),fNJets->At(i)); // always normalize corr level 0 = raw FF
909 NormalizeTH1(fCorrFF[0]->GetZ(i),fNJets->At(i));
910 NormalizeTH1(fCorrFF[0]->GetXi(i),fNJets->At(i));
914 //______________________________________________________
915 void AliFragmentationFunctionCorrections::NormalizeBgr()
917 // normalize bgr/UE FF
919 if(fNCorrectionLevelsBgr>1){
920 Printf("%s:%d -- FF normalization should be done BEFORE any correction !!!", (char*)__FILE__,__LINE__);
923 for(Int_t i=0; i<fNJetPtSlices; i++){
924 NormalizeTH1(fCorrBgr[0]->GetTrackPt(i), fNJetsBgr->At(i)); // always normalize corr level 0 = raw FF
925 NormalizeTH1(fCorrBgr[0]->GetZ(i), fNJetsBgr->At(i));
926 NormalizeTH1(fCorrBgr[0]->GetXi(i),fNJetsBgr->At(i));
931 //__________________________________________________________________________________________________
932 void AliFragmentationFunctionCorrections::ReadRawFF(TString strfile, TString strID, TString strFFID)
934 // read raw FF - standard dir/list name
936 TString strdir = "PWGJE_FragmentationFunction_" + strID;
937 TString strlist = "fracfunc_" + strID;
939 ReadRawFF(strfile,strdir,strlist,strFFID);
942 //____________________________________________________________________________________________________________________
943 void AliFragmentationFunctionCorrections::ReadRawFF(TString strfile, TString strdir, TString strlist, TString strFFID)
945 // get raw FF from input file, project in jet pt slice
946 // normalization done separately
948 TFile f(strfile,"READ");
951 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
955 if(fDebug>0) Printf("%s:%d -- read FF from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
957 gDirectory->cd(strdir);
961 if(!(list = (TList*) gDirectory->Get(strlist))){
962 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
966 TString hnameJetPt(Form("fh1FFJetPt%s",strFFID.Data()));
967 TString hnameTrackPt(Form("fh2FFTrackPt%s",strFFID.Data()));
968 TString hnameZ(Form("fh2FFZ%s",strFFID.Data()));
969 TString hnameXi(Form("fh2FFXi%s",strFFID.Data()));
971 TH1F* fh1FFJetPt = (TH1F*) list->FindObject(hnameJetPt);
972 TH2F* fh2FFTrackPt = (TH2F*) list->FindObject(hnameTrackPt);
973 TH2F* fh2FFZ = (TH2F*) list->FindObject(hnameZ);
974 TH2F* fh2FFXi = (TH2F*) list->FindObject(hnameXi);
976 if(!fh1FFJetPt) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameJetPt.Data()); return; }
977 if(!fh2FFTrackPt){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameTrackPt.Data()); return; }
978 if(!fh2FFZ) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameZ.Data()); return; }
979 if(!fh2FFXi) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameXi.Data()); return; }
981 fh1FFJetPt->SetDirectory(0);
982 fh2FFTrackPt->SetDirectory(0);
983 fh2FFZ->SetDirectory(0);
984 fh2FFXi->SetDirectory(0);
991 for(Int_t i=0; i<fNJetPtSlices; i++){
993 Float_t jetPtLoLim = fJetPtSlices->At(i);
994 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
996 Int_t binLo = static_cast<Int_t>(fh1FFJetPt->FindBin(jetPtLoLim));
997 Int_t binUp = static_cast<Int_t>(fh1FFJetPt->FindBin(jetPtUpLim)) - 1;
999 Float_t nJetsBin = fh1FFJetPt->Integral(binLo,binUp);
1001 fNJets->SetAt(nJetsBin,i);
1003 if(fDebug>0) Printf("jet pt %d to %d: nJets %f",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),fNJets->At(i));
1008 for(Int_t i=0; i<fNJetPtSlices; i++){
1010 Float_t jetPtLoLim = fJetPtSlices->At(i);
1011 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1013 Int_t binLo = static_cast<Int_t>(fh2FFTrackPt->GetXaxis()->FindBin(jetPtLoLim));
1014 Int_t binUp = static_cast<Int_t>(fh2FFTrackPt->GetXaxis()->FindBin(jetPtUpLim))-1;
1016 if(binUp > fh2FFTrackPt->GetNbinsX()){
1017 Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
1021 TString strNameFFPt(Form("fh1FFTrackPt%s_%02d_%02d",strFFID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1022 TString strNameFFZ(Form("fh1FFZ%s_%02d_%02d",strFFID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1023 TString strNameFFXi(Form("fh1FFXi%s_%02d_%02d",strFFID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1025 // appendix 'unbinned' to avoid histos with same name after rebinning
1026 TH1F* projPt = (TH1F*) fh2FFTrackPt->ProjectionY(strNameFFPt+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
1027 TH1F* projZ = (TH1F*) fh2FFZ->ProjectionY(strNameFFZ+"_unBinned",binLo,binUp,"o");
1028 TH1F* projXi = (TH1F*) fh2FFXi->ProjectionY(strNameFFXi+"_unBinned",binLo,binUp,"o");
1030 if(fNHistoBinsPt[i]) projPt = (TH1F*) projPt->Rebin(fNHistoBinsPt[i],strNameFFPt,fHistoBinsPt[i]->GetArray());
1031 if(fNHistoBinsZ[i]) projZ = (TH1F*) projZ->Rebin(fNHistoBinsZ[i],strNameFFZ,fHistoBinsZ[i]->GetArray());
1032 if(fNHistoBinsXi[i]) projXi = (TH1F*) projXi->Rebin(fNHistoBinsXi[i],strNameFFXi,fHistoBinsXi[i]->GetArray());
1034 projPt->SetNameTitle(strNameFFPt,"");
1035 projZ->SetNameTitle(strNameFFZ,"");
1036 projXi->SetNameTitle(strNameFFXi,"");
1038 // raw FF = corr level 0
1039 fCorrFF[0]->AddCorrHistos(i,projPt,projZ,projXi);
1043 //_____________________________________________________________________________________________________________________
1044 void AliFragmentationFunctionCorrections::ReadRawBgr(TString strfile, TString strID, TString strBgrID, TString strFFID)
1046 // read raw FF - standard dir/list name
1048 TString strdir = "PWGJE_FragmentationFunction_" + strID;
1049 TString strlist = "fracfunc_" + strID;
1051 ReadRawBgr(strfile,strdir,strlist,strBgrID,strFFID);
1054 //_______________________________________________________________________________________________________________________________________
1055 void AliFragmentationFunctionCorrections::ReadRawBgr(TString strfile, TString strdir, TString strlist, TString strBgrID, TString strFFID)
1057 // get raw FF from input file, project in jet pt slice
1058 // use jet dN/dpt corresponding to strFFID, bgr FF to strBgrID+strFFID
1059 // e.g. "fh1FFJetPtRecCuts", "fh2FFXiBgrPerpRecCuts"
1060 // normalization done separately
1062 TString strID = strBgrID + strFFID;
1064 TFile f(strfile,"READ");
1067 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1071 if(fDebug>0) Printf("%s:%d -- read Bgr %s from file %s ",(char*)__FILE__,__LINE__,strBgrID.Data(),strfile.Data());
1073 gDirectory->cd(strdir);
1077 if(!(list = (TList*) gDirectory->Get(strlist))){
1078 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1082 TString hnameNJets = "fh1nRecJetsCuts";
1083 TString hnameJetPt(Form("fh1FFJetPt%s",strFFID.Data())); // not: strID.Data() !!! would not be proper normalization
1084 TString hnameBgrTrackPt(Form("fh2FFTrackPt%s",strID.Data()));
1085 TString hnameBgrZ(Form("fh2FFZ%s",strID.Data()));
1086 TString hnameBgrXi(Form("fh2FFXi%s",strID.Data()));
1088 TH1F* fh1NJets = (TH1F*) list->FindObject(hnameNJets); // needed for normalization of bgr out of 2 jets
1089 TH1F* fh1FFJetPtBgr = (TH1F*) list->FindObject(hnameJetPt);
1090 TH2F* fh2FFTrackPtBgr = (TH2F*) list->FindObject(hnameBgrTrackPt);
1091 TH2F* fh2FFZBgr = (TH2F*) list->FindObject(hnameBgrZ);
1092 TH2F* fh2FFXiBgr = (TH2F*) list->FindObject(hnameBgrXi);
1094 if(!fh1FFJetPtBgr) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameJetPt.Data()); return; }
1095 if(!fh1NJets) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameNJets.Data()); return; }
1096 if(!fh2FFTrackPtBgr){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrTrackPt.Data()); return; }
1097 if(!fh2FFZBgr) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrZ.Data()); return; }
1098 if(!fh2FFXiBgr) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrXi.Data()); return; }
1100 fh1FFJetPtBgr->SetDirectory(0);
1101 fh1NJets->SetDirectory(0);
1102 fh2FFTrackPtBgr->SetDirectory(0);
1103 fh2FFZBgr->SetDirectory(0);
1104 fh2FFXiBgr->SetDirectory(0);
1110 for(Int_t i=0; i<fNJetPtSlices; i++){
1112 Float_t jetPtLoLim = fJetPtSlices->At(i);
1113 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1115 Int_t binLo = static_cast<Int_t>(fh1FFJetPtBgr->FindBin(jetPtLoLim));
1116 Int_t binUp = static_cast<Int_t>(fh1FFJetPtBgr->FindBin(jetPtUpLim)) - 1;
1118 Float_t nJetsBin = fh1FFJetPtBgr->Integral(binLo,binUp);
1119 Double_t scaleF = 1;
1121 //if(strBgrID.Contains("Out2Jets")){ // scale by ratio 2 jets events / all events
1122 // scaleF = fh1NJets->Integral(fh1NJets->FindBin(2),fh1NJets->GetNbinsX())
1123 // / fh1NJets->Integral(fh1NJets->FindBin(1),fh1NJets->GetNbinsX());}
1126 if(strBgrID.Contains("OutAllJets")){ // scale by ratio >3 jets events / all events
1127 scaleF = fh1NJets->Integral(fh1NJets->FindBin(4),fh1NJets->GetNbinsX())
1128 / fh1NJets->Integral(fh1NJets->FindBin(1),fh1NJets->GetNbinsX());
1131 fNJetsBgr->SetAt(nJetsBin*scaleF,i);
1133 if(fDebug>0) Printf("bgr jet pt %d to %d: nJets %f, scaleF %.2f",
1134 static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),nJetsBin,scaleF);
1140 for(Int_t i=0; i<fNJetPtSlices; i++){
1142 Float_t jetPtLoLim = fJetPtSlices->At(i);
1143 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1145 Int_t binLo = static_cast<Int_t>(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtLoLim));
1146 Int_t binUp = static_cast<Int_t>(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtUpLim))-1;
1148 if(binUp > fh2FFTrackPtBgr->GetNbinsX()){
1149 Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
1153 TString strNameBgrPt(Form("fh1BgrTrackPt%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1154 TString strNameBgrZ(Form("fh1BgrZ%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1155 TString strNameBgrXi(Form("fh1BgrXi%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1157 // appendix 'unbinned' to avoid histos with same name after rebinning
1158 TH1F* projPt = (TH1F*) fh2FFTrackPtBgr->ProjectionY(strNameBgrPt+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
1159 TH1F* projZ = (TH1F*) fh2FFZBgr->ProjectionY(strNameBgrZ+"_unBinned",binLo,binUp,"o");
1160 TH1F* projXi = (TH1F*) fh2FFXiBgr->ProjectionY(strNameBgrXi+"_unBinned",binLo,binUp,"o");
1162 if(fNHistoBinsPt[i]) projPt = (TH1F*) projPt->Rebin(fNHistoBinsPt[i],strNameBgrPt,fHistoBinsPt[i]->GetArray());
1163 if(fNHistoBinsZ[i]) projZ = (TH1F*) projZ->Rebin(fNHistoBinsZ[i],strNameBgrZ,fHistoBinsZ[i]->GetArray());
1164 if(fNHistoBinsXi[i]) projXi = (TH1F*) projXi->Rebin(fNHistoBinsXi[i],strNameBgrXi,fHistoBinsXi[i]->GetArray());
1166 projPt->SetNameTitle(strNameBgrPt,"");
1167 projZ->SetNameTitle(strNameBgrZ,"");
1168 projXi->SetNameTitle(strNameBgrXi,"");
1170 // raw bgr = corr level 0
1171 fCorrBgr[0]->AddCorrHistos(i,projPt,projZ,projXi);
1175 //_____________________________________________________________________________________________________________________
1176 void AliFragmentationFunctionCorrections::ReadRawBgrEmbedding(TString strfile, TString strID, TString strFFID)
1178 // read raw FF - standard dir/list name
1180 TString strdir = "PWGJE_FragmentationFunction_" + strID;
1181 TString strlist = "fracfunc_" + strID;
1183 ReadRawBgrEmbedding(strfile,strdir,strlist,strFFID);
1186 //_______________________________________________________________________________________________________________________________________
1187 void AliFragmentationFunctionCorrections::ReadRawBgrEmbedding(TString strfile, TString strdir, TString strlist, TString strFFID)
1189 // get raw FF from input file, project in jet pt slice
1190 // for embedding, the bgr FF are taken from histos "fh1FFJetPtRecCuts", "fh2FFXiRecCuts"
1191 // normalization done separately
1193 TString strBgrID = "BckgEmbed";
1194 TString strID = strBgrID + strFFID;
1196 TFile f(strfile,"READ");
1199 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1203 if(fDebug>0) Printf("%s:%d -- read Bgr %s from file %s ",(char*)__FILE__,__LINE__,strFFID.Data(),strfile.Data());
1205 gDirectory->cd(strdir);
1209 if(!(list = (TList*) gDirectory->Get(strlist))){
1210 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1214 TString hnameNJets = "fh1nRecJetsCuts";
1215 TString hnameJetPt(Form("fh1FFJetPt%s",strFFID.Data()));
1216 TString hnameBgrTrackPt(Form("fh2FFTrackPt%s",strFFID.Data()));
1217 TString hnameBgrZ(Form("fh2FFZ%s",strFFID.Data()));
1218 TString hnameBgrXi(Form("fh2FFXi%s",strFFID.Data()));
1220 TH1F* fh1NJets = (TH1F*) list->FindObject(hnameNJets); // needed for normalization of bgr out of 2 jets
1221 TH1F* fh1FFJetPtBgr = (TH1F*) list->FindObject(hnameJetPt);
1222 TH2F* fh2FFTrackPtBgr = (TH2F*) list->FindObject(hnameBgrTrackPt);
1223 TH2F* fh2FFZBgr = (TH2F*) list->FindObject(hnameBgrZ);
1224 TH2F* fh2FFXiBgr = (TH2F*) list->FindObject(hnameBgrXi);
1226 if(!fh1FFJetPtBgr) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameJetPt.Data()); return; }
1227 if(!fh1NJets) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameNJets.Data()); return; }
1228 if(!fh2FFTrackPtBgr){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrTrackPt.Data()); return; }
1229 if(!fh2FFZBgr) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrZ.Data()); return; }
1230 if(!fh2FFXiBgr) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrXi.Data()); return; }
1232 fh1FFJetPtBgr->SetDirectory(0);
1233 fh1NJets->SetDirectory(0);
1234 fh2FFTrackPtBgr->SetDirectory(0);
1235 fh2FFZBgr->SetDirectory(0);
1236 fh2FFXiBgr->SetDirectory(0);
1242 for(Int_t i=0; i<fNJetPtSlices; i++){
1244 Float_t jetPtLoLim = fJetPtSlices->At(i);
1245 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1247 Int_t binLo = static_cast<Int_t>(fh1FFJetPtBgr->FindBin(jetPtLoLim));
1248 Int_t binUp = static_cast<Int_t>(fh1FFJetPtBgr->FindBin(jetPtUpLim)) - 1;
1250 Float_t nJetsBin = fh1FFJetPtBgr->Integral(binLo,binUp);
1251 Double_t scaleF = 1;
1253 fNJetsBgr->SetAt(nJetsBin*scaleF,i);
1255 if(fDebug>0) Printf("bgr jet pt %d to %d: nJets %f, scaleF %.2f",
1256 static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),nJetsBin,scaleF);
1262 for(Int_t i=0; i<fNJetPtSlices; i++){
1264 Float_t jetPtLoLim = fJetPtSlices->At(i);
1265 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1267 Int_t binLo = static_cast<Int_t>(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtLoLim));
1268 Int_t binUp = static_cast<Int_t>(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtUpLim))-1;
1270 if(binUp > fh2FFTrackPtBgr->GetNbinsX()){
1271 Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
1275 TString strNameBgrPt(Form("fh1BgrTrackPt%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1276 TString strNameBgrZ(Form("fh1BgrZ%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1277 TString strNameBgrXi(Form("fh1BgrXi%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1279 // appendix 'unbinned' to avoid histos with same name after rebinning
1280 TH1F* projPt = (TH1F*) fh2FFTrackPtBgr->ProjectionY(strNameBgrPt+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
1281 TH1F* projZ = (TH1F*) fh2FFZBgr->ProjectionY(strNameBgrZ+"_unBinned",binLo,binUp,"o");
1282 TH1F* projXi = (TH1F*) fh2FFXiBgr->ProjectionY(strNameBgrXi+"_unBinned",binLo,binUp,"o");
1284 if(fNHistoBinsPt[i]) projPt = (TH1F*) projPt->Rebin(fNHistoBinsPt[i],strNameBgrPt,fHistoBinsPt[i]->GetArray());
1285 if(fNHistoBinsZ[i]) projZ = (TH1F*) projZ->Rebin(fNHistoBinsZ[i],strNameBgrZ,fHistoBinsZ[i]->GetArray());
1286 if(fNHistoBinsXi[i]) projXi = (TH1F*) projXi->Rebin(fNHistoBinsXi[i],strNameBgrXi,fHistoBinsXi[i]->GetArray());
1288 projPt->SetNameTitle(strNameBgrPt,"");
1289 projZ->SetNameTitle(strNameBgrZ,"");
1290 projXi->SetNameTitle(strNameBgrXi,"");
1292 // raw bgr = corr level 0
1293 fCorrBgr[0]->AddCorrHistos(i,projPt,projZ,projXi);
1298 //__________________________________________________________________________________________________________
1299 void AliFragmentationFunctionCorrections::WriteOutput(TString strfile, TString strdir, Bool_t updateOutfile)
1301 // write histos to file
1302 // skip histos with 0 entries
1304 TString outfileOption = "RECREATE";
1305 if(updateOutfile) outfileOption = "UPDATE";
1307 TFile f(strfile,outfileOption);
1310 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1314 if(fDebug>0) Printf("%s:%d -- write FF to file %s ",(char*)__FILE__,__LINE__,strfile.Data());
1316 if(strdir && strdir.Length()){
1317 TDirectory* dir = f.mkdir(strdir);
1321 for(Int_t i=0; i<fNJetPtSlices; i++){
1323 for(Int_t c=0; c<fNCorrectionLevels; c++) if(fCorrFF[c]->GetTrackPt(i)->GetEntries()) fCorrFF[c]->GetTrackPt(i)->Write();
1324 for(Int_t c=0; c<fNCorrectionLevels; c++) if(fCorrFF[c]->GetZ(i)->GetEntries()) fCorrFF[c]->GetZ(i)->Write();
1325 for(Int_t c=0; c<fNCorrectionLevels; c++) if(fCorrFF[c]->GetXi(i)->GetEntries()) fCorrFF[c]->GetXi(i)->Write();
1327 if(fh1FFXiShift[i]->GetEntries()) fh1FFXiShift[i]->Write();
1329 for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) if(fCorrBgr[c]->GetTrackPt(i)->GetEntries()) fCorrBgr[c]->GetTrackPt(i)->Write();
1330 for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) if(fCorrBgr[c]->GetZ(i)->GetEntries()) fCorrBgr[c]->GetZ(i)->Write();
1331 for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) if(fCorrBgr[c]->GetXi(i)->GetEntries()) fCorrBgr[c]->GetXi(i)->Write();
1334 if(fh1FFTrackPtBackFolded[i] && fh1FFTrackPtBackFolded[i]->GetEntries()) fh1FFTrackPtBackFolded[i]->Write();
1335 if(fh1FFZBackFolded[i] && fh1FFZBackFolded[i]->GetEntries()) fh1FFZBackFolded[i]->Write();
1336 if(fh1FFXiBackFolded[i] && fh1FFXiBackFolded[i]->GetEntries()) fh1FFXiBackFolded[i]->Write();
1339 if(fh1FFRatioTrackPtFolded[i] && fh1FFRatioTrackPtFolded[i]->GetEntries()) fh1FFRatioTrackPtFolded[i]->Write();
1340 if(fh1FFRatioZFolded[i] && fh1FFRatioZFolded[i]->GetEntries()) fh1FFRatioZFolded[i]->Write();
1341 if(fh1FFRatioXiFolded[i] && fh1FFRatioXiFolded[i]->GetEntries()) fh1FFRatioXiFolded[i]->Write();
1343 if(fh1FFRatioTrackPtBackFolded[i] && fh1FFRatioTrackPtBackFolded[i]->GetEntries()) fh1FFRatioTrackPtBackFolded[i]->Write();
1344 if(fh1FFRatioZBackFolded[i] && fh1FFRatioZBackFolded[i]->GetEntries()) fh1FFRatioZBackFolded[i]->Write();
1345 if(fh1FFRatioXiBackFolded[i] &&fh1FFRatioXiBackFolded[i]->GetEntries()) fh1FFRatioXiBackFolded[i]->Write();
1349 // inclusive track pt
1351 for(Int_t c=0; c<fNCorrectionLevelsSinglePt; c++) if(fCorrSinglePt[c]->GetTrackPt(0)->GetEntries()) fCorrSinglePt[c]->GetTrackPt(0)->Write();
1352 if(fh1SingleTrackPtBackFolded) fh1SingleTrackPtBackFolded->Write();
1353 if(fh1RatioSingleTrackPtFolded) fh1RatioSingleTrackPtFolded->Write();
1354 if(fh1RatioSingleTrackPtBackFolded) fh1RatioSingleTrackPtBackFolded->Write();
1359 //____________________________________________________________________________________________________________________________________
1360 THnSparse* AliFragmentationFunctionCorrections::TH1toSparse(const TH1F* hist, TString strName, TString strTit, const Bool_t fillConst)
1362 // copy 1-dimensional histo to THnSparse
1363 // if fillConst TRUE, create THnSparse with same binning as hist but all entries = 1
1364 // histos with variable bin size are supported
1366 // note: function returns pointer to 'new' THnSparse on heap, object needs to be deleted by user
1368 THnSparseF* fhnHist;
1370 Int_t nBins = hist->GetXaxis()->GetNbins();
1371 Int_t nBinsVec[1] = { nBins };
1373 const Double_t* binsVec = hist->GetXaxis()->GetXbins()->GetArray();
1375 if(binsVec){ // binsVec only neq NULL if histo was rebinned before
1377 fhnHist = new THnSparseF(strName,strTit,1,nBinsVec/*,binMinVec,binMaxVec*/);
1378 fhnHist->SetBinEdges(0,binsVec);
1380 else{ // uniform bin size
1382 Double_t xMin = hist->GetXaxis()->GetXmin();
1383 Double_t xMax = hist->GetXaxis()->GetXmax();
1385 Double_t binMinVec[1] = { xMin };
1386 Double_t binMaxVec[1] = { xMax };
1388 fhnHist = new THnSparseF(strName,strTit,1,nBinsVec,binMinVec,binMaxVec);
1392 for(Int_t bin=1; bin<=nBins; bin++){
1394 Double_t binCenter = fhnHist->GetAxis(0)->GetBinCenter(bin);
1396 Double_t binCoord[] = {binCenter};
1397 fhnHist->Fill(binCoord,1); // initially need to create the bin
1399 Long64_t binIndex = fhnHist->GetBin(binCoord);
1401 Double_t cont = hist->GetBinContent(bin);
1402 Double_t err = hist->GetBinError(bin);
1409 fhnHist->SetBinContent(binIndex,cont);
1410 fhnHist->SetBinError(binIndex,err);
1416 //______________________________________________________________________________________________________________________________________________
1417 THnSparse* AliFragmentationFunctionCorrections::Unfold(THnSparse* hnHist, const THnSparse* hnResponse, const THnSparse* hnEff, const Int_t nIter,
1418 const Bool_t useCorrelatedErrors, const THnSparse* hnPrior)
1420 // unfold input histo
1422 AliCFUnfolding unfolding("unfolding","",1,hnResponse,hnEff,hnHist,hnPrior); // arg3: nVar; hnEff required, hnPrior not (defaults to 0x0)
1423 unfolding.SetMaxNumberOfIterations(nIter);
1424 // unfolding.SetMaxChi2PerDOF(1.e-07); // OBSOLETE !!!
1425 // if(useSmoothing) unfolding.UseSmoothing();
1427 if(useCorrelatedErrors) unfolding.SetUseCorrelatedErrors();
1430 THnSparse* unfolded = unfolding.GetUnfolded();
1432 TString hnameUnf = hnHist->GetName();
1433 hnameUnf.Append("_unf");
1434 unfolded->SetNameTitle(hnameUnf,hnHist->GetTitle());
1439 //___________________________________________________________________________________________________________________________
1440 void AliFragmentationFunctionCorrections::UnfoldHistos(const Int_t nIter, const Bool_t useCorrelatedErrors, const Int_t type)
1442 // loop over jet pt slices and unfold dN/dpt spectra
1444 TString labelCorr = fCorrFF[fNCorrectionLevels-1]->GetLabel();
1445 if(!labelCorr.Contains("unfold")) AddCorrectionLevel("unfold");
1447 for(Int_t i=0; i<fNJetPtSlices; i++){
1450 if(type == kFlagPt) hist = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i); // level -2: before unfolding, level -1: unfolded
1451 else if(type == kFlagZ) hist = fCorrFF[fNCorrectionLevels-2]->GetZ(i); // level -2: before unfolding, level -1: unfolded
1452 else if(type == kFlagXi) hist = fCorrFF[fNCorrectionLevels-2]->GetXi(i); // level -2: before unfolding, level -1: unfolded
1454 Printf("%s%d unknown type",(char*)__FILE__,__LINE__);
1459 Printf("%s%d no histo found ",(char*)__FILE__,__LINE__);
1464 THnSparse* hnResponse = 0;
1465 if(type == kFlagPt) hnResponse = fhnResponsePt[i];
1466 else if(type == kFlagZ) hnResponse = fhnResponseZ[i];
1467 else if(type == kFlagXi) hnResponse = fhnResponseXi[i];
1470 if(type == kFlagPt && fh1FFTrackPtPrior[i] && ((TString(fh1FFTrackPtPrior[i]->GetName())).Length() > 0) ) hPrior = fh1FFTrackPtPrior[i];
1471 else if(type == kFlagZ && fh1FFZPrior[i] && ((TString(fh1FFZPrior[i]->GetName())).Length() > 0) ) hPrior = fh1FFZPrior[i];
1472 else if(type == kFlagXi && fh1FFXiPrior[i] && ((TString(fh1FFXiPrior[i]->GetName())).Length() > 0) ) hPrior = fh1FFXiPrior[i];
1475 TString histNameTHn;
1476 histNameTHn = hist->GetName();
1477 histNameTHn.ReplaceAll("TH1","THn");
1480 TString priorNameTHn;
1482 priorNameTHn = hPrior->GetName();
1483 priorNameTHn.ReplaceAll("TH1","THn");
1486 TString histNameBackFolded;
1487 histNameBackFolded = hist->GetName();
1488 histNameBackFolded.Append("_backfold");
1490 TString histNameRatioFolded;
1491 histNameRatioFolded = hist->GetName();
1492 histNameRatioFolded.ReplaceAll("fh1FF","hRatioFF");
1494 histNameRatioFolded.Append("_unfold");
1496 TString histNameRatioBackFolded;
1497 histNameRatioBackFolded = hist->GetName();
1498 histNameRatioBackFolded.ReplaceAll("fh1FF","hRatioFF");
1499 histNameRatioBackFolded.Append("_backfold");
1501 THnSparse* hnHist = TH1toSparse(hist,histNameTHn,hist->GetTitle());
1502 THnSparse* hnFlatEfficiency = TH1toSparse(hist,"fhnEfficiency","eff",kTRUE); // could optionally also use real eff
1503 THnSparse* hnPrior = 0;
1504 if(hPrior) hnPrior = TH1toSparse(hPrior,priorNameTHn,hPrior->GetTitle());
1506 THnSparse* hnUnfolded
1507 = Unfold(hnHist,hnResponse,hnFlatEfficiency,nIter,useCorrelatedErrors,hnPrior);
1509 TH1F* hUnfolded = (TH1F*) hnUnfolded->Projection(0);
1510 hUnfolded->SetNameTitle(hist->GetName(),hist->GetTitle());
1512 if(type == kFlagPt) fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hUnfolded,0,0);
1513 if(type == kFlagZ) fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,0,hUnfolded,0);
1514 if(type == kFlagXi) fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,0,0,hUnfolded);
1516 // backfolding: apply response matrix to unfolded spectrum
1517 TH1F* hBackFolded = ApplyResponse(hUnfolded,hnResponse);
1518 hBackFolded->SetNameTitle(histNameBackFolded,hist->GetTitle());
1520 if(type == kFlagPt) fh1FFTrackPtBackFolded[i] = hBackFolded;
1521 if(type == kFlagZ) fh1FFZBackFolded[i] = hBackFolded;
1522 if(type == kFlagXi) fh1FFXiBackFolded[i] = hBackFolded;
1524 // ratio unfolded to original histo
1525 TH1F* hRatioUnfolded = (TH1F*) hUnfolded->Clone(histNameRatioFolded);
1526 hRatioUnfolded->Reset();
1527 hRatioUnfolded->Divide(hUnfolded,hist,1,1,"B");
1529 if(type == kFlagPt) fh1FFRatioTrackPtFolded[i] = hRatioUnfolded;
1530 if(type == kFlagZ) fh1FFRatioZFolded[i] = hRatioUnfolded;
1531 if(type == kFlagXi) fh1FFRatioXiFolded[i] = hRatioUnfolded;
1534 // ratio backfolded to original histo
1535 TH1F* hRatioBackFolded = (TH1F*) hBackFolded->Clone(histNameRatioBackFolded);
1536 hRatioBackFolded->Reset();
1537 hRatioBackFolded->Divide(hBackFolded,hist,1,1,"B");
1539 if(type == kFlagPt) fh1FFRatioTrackPtBackFolded[i] = hRatioBackFolded;
1540 if(type == kFlagZ) fh1FFRatioZBackFolded[i] = hRatioBackFolded;
1541 if(type == kFlagXi) fh1FFRatioXiBackFolded[i] = hRatioBackFolded;
1544 delete hnFlatEfficiency;
1549 //_____________________________________________________________________________________________________
1550 void AliFragmentationFunctionCorrections::UnfoldPt(const Int_t nIter, const Bool_t useCorrelatedErrors)
1553 Int_t type = kFlagPt;
1554 UnfoldHistos(nIter, useCorrelatedErrors, type);
1557 //_____________________________________________________________________________________________________
1558 void AliFragmentationFunctionCorrections::UnfoldZ(const Int_t nIter, const Bool_t useCorrelatedErrors)
1561 Int_t type = kFlagZ;
1562 UnfoldHistos(nIter, useCorrelatedErrors, type);
1565 //_____________________________________________________________________________________________________
1566 void AliFragmentationFunctionCorrections::UnfoldXi(const Int_t nIter, const Bool_t useCorrelatedErrors)
1569 Int_t type = kFlagXi;
1570 UnfoldHistos(nIter, useCorrelatedErrors, type);
1573 //______________________________________________________________________________________________
1574 TH1F* AliFragmentationFunctionCorrections::ApplyResponse(const TH1F* hist, THnSparse* hnResponse)
1576 // apply (multiply) response matrix to hist
1578 TH1F* hOut = new TH1F(*hist);
1581 const Int_t axisM = 0;
1582 const Int_t axisT = 1;
1584 Int_t nBinsM = hnResponse->GetAxis(axisM)->GetNbins();
1585 Int_t nBinsT = hnResponse->GetAxis(axisT)->GetNbins();
1587 // response matrix normalization
1588 // do it in this function and not when reading response, since for 'proper' normalization errors are difficult to assign
1589 // so for unfolding proper we leave it to CORRFW ...
1591 Double_t normFacResponse[nBinsT];
1593 for(Int_t iT=1; iT<=nBinsT; iT++){
1595 Double_t sumResp = 0;
1597 for(Int_t iM=1; iM<=nBinsM; iM++){
1599 Double_t coordM = hnResponse->GetAxis(axisM)->GetBinCenter(iM);
1600 Double_t coordT = hnResponse->GetAxis(axisT)->GetBinCenter(iT);
1602 Double_t binCoord[] = {coordM,coordT};
1604 Long64_t binIndex = hnResponse->GetBin(binCoord);
1606 Double_t resp = hnResponse->GetBinContent(binIndex);
1611 normFacResponse[iT] = 0;
1612 if(sumResp) normFacResponse[iT] = 1/sumResp;
1617 for(Int_t iM=1; iM<=nBinsM; iM++){
1621 for(Int_t iT=1; iT<=nBinsT; iT++){
1623 Double_t contT = hist->GetBinContent(iT);
1625 Double_t coordM = hnResponse->GetAxis(axisM)->GetBinCenter(iM);
1626 Double_t coordT = hnResponse->GetAxis(axisT)->GetBinCenter(iT);
1628 Double_t binCoord[] = {coordM,coordT};
1630 Long64_t binIndex = hnResponse->GetBin(binCoord);
1632 Double_t resp = hnResponse->GetBinContent(binIndex);
1634 contM += resp*normFacResponse[iT]*contT;
1637 hOut->SetBinContent(iM,contM);
1643 //_______________________________________________________________________________________________________
1644 void AliFragmentationFunctionCorrections::ReadEfficiency(TString strfile, TString strdir, TString strlist)
1646 // read reconstruction efficiency from file
1647 // argument strlist optional - read from directory strdir if not specified
1649 // temporary histos to hold histos from file
1650 TH1F* hEffPt[fNJetPtSlices];
1651 TH1F* hEffZ[fNJetPtSlices];
1652 TH1F* hEffXi[fNJetPtSlices];
1654 for(Int_t i=0; i<fNJetPtSlices; i++) hEffPt[i] = 0;
1655 for(Int_t i=0; i<fNJetPtSlices; i++) hEffZ[i] = 0;
1656 for(Int_t i=0; i<fNJetPtSlices; i++) hEffXi[i] = 0;
1658 TFile f(strfile,"READ");
1661 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1665 if(fDebug>0) Printf("%s:%d -- read efficiencies from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
1667 if(strdir && strdir.Length()) gDirectory->cd(strdir);
1671 if(strlist && strlist.Length()){
1673 if(!(list = (TList*) gDirectory->Get(strlist))){
1674 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1679 for(Int_t i=0; i<fNJetPtSlices; i++){
1681 Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
1682 Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
1684 TString strNameEffPt(Form("hEffPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
1685 TString strNameEffZ(Form("hEffZ_%02d_%02d",jetPtLoLim,jetPtUpLim));
1686 TString strNameEffXi(Form("hEffXi_%02d_%02d",jetPtLoLim,jetPtUpLim));
1690 hEffPt[i] = (TH1F*) list->FindObject(strNameEffPt);
1691 hEffZ[i] = (TH1F*) list->FindObject(strNameEffZ);
1692 hEffXi[i] = (TH1F*) list->FindObject(strNameEffXi);
1695 hEffPt[i] = (TH1F*) gDirectory->Get(strNameEffPt);
1696 hEffZ[i] = (TH1F*) gDirectory->Get(strNameEffZ);
1697 hEffXi[i] = (TH1F*) gDirectory->Get(strNameEffXi);
1701 Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffPt.Data());
1705 Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffZ.Data());
1709 Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffXi.Data());
1713 if(fNHistoBinsPt[i]) hEffPt[i] = (TH1F*) hEffPt[i]->Rebin(fNHistoBinsPt[i],strNameEffPt+"_rebin",fHistoBinsPt[i]->GetArray());
1714 if(fNHistoBinsZ[i]) hEffZ[i] = (TH1F*) hEffZ[i]->Rebin(fNHistoBinsZ[i],strNameEffZ+"_rebin",fHistoBinsZ[i]->GetArray());
1715 if(fNHistoBinsXi[i]) hEffXi[i] = (TH1F*) hEffXi[i]->Rebin(fNHistoBinsXi[i],strNameEffXi+"_rebin",fHistoBinsXi[i]->GetArray());
1717 if(hEffPt[i]) hEffPt[i]->SetDirectory(0);
1718 if(hEffZ[i]) hEffZ[i]->SetDirectory(0);
1719 if(hEffXi[i]) hEffXi[i]->SetDirectory(0);
1721 } // jet slices loop
1725 for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
1726 if(hEffPt[i]) new(fh1EffPt[i]) TH1F(*hEffPt[i]);
1727 if(hEffZ[i]) new(fh1EffZ[i]) TH1F(*hEffZ[i]);
1728 if(hEffXi[i]) new(fh1EffXi[i]) TH1F(*hEffXi[i]);
1732 //___________________________________________________________________________________________________________
1733 void AliFragmentationFunctionCorrections::ReadBgrEfficiency(TString strfile, TString strdir, TString strlist)
1735 // read bgr eff from file
1736 // argument strlist optional - read from directory strdir if not specified
1738 TH1F* hEffPtBgr[fNJetPtSlices];
1739 TH1F* hEffZBgr [fNJetPtSlices];
1740 TH1F* hEffXiBgr[fNJetPtSlices];
1742 for(Int_t i=0; i<fNJetPtSlices; i++) hEffPtBgr[i] = 0;
1743 for(Int_t i=0; i<fNJetPtSlices; i++) hEffZBgr[i] = 0;
1744 for(Int_t i=0; i<fNJetPtSlices; i++) hEffXiBgr[i] = 0;
1747 TFile f(strfile,"READ");
1750 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1754 if(fDebug>0) Printf("%s:%d -- read bgr efficiencies from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
1756 if(strdir && strdir.Length()) gDirectory->cd(strdir);
1760 if(strlist && strlist.Length()){
1762 if(!(list = (TList*) gDirectory->Get(strlist))){
1763 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1768 for(Int_t i=0; i<fNJetPtSlices; i++){
1770 Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
1771 Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
1773 TString strNameEffPtBgr(Form("hEffPtBgr%02dto%02d",jetPtLoLim,jetPtUpLim));
1774 TString strNameEffZBgr(Form("hEffZBgr%02dto%02d",jetPtLoLim,jetPtUpLim));
1775 TString strNameEffXiBgr(Form("hEffXiBgr%02dto%02d",jetPtLoLim,jetPtUpLim));
1779 hEffPtBgr[i] = (TH1F*) list->FindObject(strNameEffPtBgr);
1780 hEffZBgr[i] = (TH1F*) list->FindObject(strNameEffZBgr);
1781 hEffXiBgr[i] = (TH1F*) list->FindObject(strNameEffXiBgr);
1784 hEffPtBgr[i] = (TH1F*) gDirectory->Get(strNameEffPtBgr);
1785 hEffZBgr[i] = (TH1F*) gDirectory->Get(strNameEffZBgr);
1786 hEffXiBgr[i] = (TH1F*) gDirectory->Get(strNameEffXiBgr);
1790 Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffPtBgr.Data());
1794 Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffZBgr.Data());
1798 Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffXiBgr.Data());
1802 if(fNHistoBinsPt[i]) hEffPtBgr[i] = (TH1F*) hEffPtBgr[i]->Rebin(fNHistoBinsPt[i],strNameEffPtBgr+"_rebin",fHistoBinsPt[i]->GetArray());
1803 if(fNHistoBinsZ[i]) hEffZBgr[i] = (TH1F*) hEffZBgr[i]->Rebin(fNHistoBinsZ[i],strNameEffZBgr+"_rebin",fHistoBinsZ[i]->GetArray());
1804 if(fNHistoBinsXi[i]) hEffXiBgr[i] = (TH1F*) hEffXiBgr[i]->Rebin(fNHistoBinsXi[i],strNameEffXiBgr+"_rebin",fHistoBinsXi[i]->GetArray());
1806 if(hEffPtBgr[i]) hEffPtBgr[i]->SetDirectory(0);
1807 if(hEffZBgr[i]) hEffZBgr[i]->SetDirectory(0);
1808 if(hEffXiBgr[i]) hEffXiBgr[i]->SetDirectory(0);
1810 } // jet slices loop
1814 for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
1815 if(hEffPtBgr[i]) new(fh1EffBgrPt[i]) TH1F(*hEffPtBgr[i]);
1816 if(hEffZBgr[i]) new(fh1EffBgrZ[i]) TH1F(*hEffZBgr[i]);
1817 if(hEffXiBgr[i]) new(fh1EffBgrXi[i]) TH1F(*hEffXiBgr[i]);
1821 // ________________________________________________
1822 void AliFragmentationFunctionCorrections::EffCorr()
1824 // apply efficiency correction
1826 AddCorrectionLevel("eff");
1828 for(Int_t i=0; i<fNJetPtSlices; i++){
1830 TH1F* histPt = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i);
1831 TH1F* histZ = fCorrFF[fNCorrectionLevels-2]->GetZ(i);
1832 TH1F* histXi = fCorrFF[fNCorrectionLevels-2]->GetXi(i);
1834 TString histNamePt = histPt->GetName();
1835 TString histNameZ = histZ->GetName();
1836 TString histNameXi = histXi->GetName();
1839 TH1F* hFFTrackPtEffCorr = (TH1F*) histPt->Clone(histNamePt);
1840 hFFTrackPtEffCorr->Divide(histPt,fh1EffPt[i],1,1,"");
1842 TH1F* hFFZEffCorr = (TH1F*) histZ->Clone(histNameZ);
1843 hFFZEffCorr->Divide(histZ,fh1EffZ[i],1,1,"");
1845 TH1F* hFFXiEffCorr = (TH1F*) histXi->Clone(histNameXi);
1846 hFFXiEffCorr->Divide(histXi,fh1EffXi[i],1,1,"");
1848 fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hFFTrackPtEffCorr,hFFZEffCorr,hFFXiEffCorr);
1852 //___________________________________________________
1853 void AliFragmentationFunctionCorrections::EffCorrBgr()
1855 // apply efficiency correction to bgr distributions
1857 AddCorrectionLevelBgr("eff");
1859 Printf("%s:%d -- apply efficiency correction, corrLevel %d",(char*)__FILE__,__LINE__,fNCorrectionLevels-1);
1861 for(Int_t i=0; i<fNJetPtSlices; i++){
1863 TH1F* histPt = fCorrBgr[fNCorrectionLevelsBgr-2]->GetTrackPt(i);
1864 TH1F* histZ = fCorrBgr[fNCorrectionLevelsBgr-2]->GetZ(i);
1865 TH1F* histXi = fCorrBgr[fNCorrectionLevelsBgr-2]->GetXi(i);
1867 TString histNamePt = histPt->GetName();
1868 TString histNameZ = histZ->GetName();
1869 TString histNameXi = histXi->GetName();
1872 TH1F* hFFTrackPtEffCorr = (TH1F*) histPt->Clone(histNamePt);
1873 hFFTrackPtEffCorr->Divide(histPt,fh1EffPt[i],1,1,"");
1875 TH1F* hFFZEffCorr = (TH1F*) histZ->Clone(histNameZ);
1876 hFFZEffCorr->Divide(histZ,fh1EffZ[i],1,1,"");
1878 TH1F* hFFXiEffCorr = (TH1F*) histXi->Clone(histNameXi);
1879 hFFXiEffCorr->Divide(histXi,fh1EffXi[i],1,1,"");
1881 fCorrBgr[fNCorrectionLevelsBgr-1]->AddCorrHistos(i,hFFTrackPtEffCorr,hFFZEffCorr,hFFXiEffCorr);
1885 //______________________________________________________________________
1886 void AliFragmentationFunctionCorrections::XiShift(const Int_t corrLevel)
1888 // re-evaluate jet energy after FF corrections from dN/dpt distribution
1889 // 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'
1890 // argument corrlevel: which level of correction to be corrected/shifted to
1892 if(corrLevel>=fNCorrectionLevels){
1893 Printf(" calc xi shift: corrLevel exceeded - do nothing");
1897 Double_t* jetPtUncorr = new Double_t[fNJetPtSlices];
1898 Double_t* jetPtCorr = new Double_t[fNJetPtSlices];
1899 Double_t* deltaXi = new Double_t[fNJetPtSlices];
1901 for(Int_t i=0; i<fNJetPtSlices; i++){
1903 TH1F* histPtRaw = fCorrFF[0]->GetTrackPt(i);
1904 TH1F* histPt = fCorrFF[corrLevel]->GetTrackPt(i);
1906 Double_t ptUncorr = 0;
1907 Double_t ptCorr = 0;
1909 for(Int_t bin = 1; bin<=histPtRaw->GetNbinsX(); bin++){
1911 Double_t cont = histPtRaw->GetBinContent(bin);
1912 Double_t width = histPtRaw->GetBinWidth(bin);
1913 Double_t meanPt = histPtRaw->GetBinCenter(bin);
1915 ptUncorr += meanPt*cont*width;
1918 for(Int_t bin = 1; bin<=histPt->GetNbinsX(); bin++){
1920 Double_t cont = histPt->GetBinContent(bin);
1921 Double_t width = histPt->GetBinWidth(bin);
1922 Double_t meanPt = histPt->GetBinCenter(bin);
1924 ptCorr += meanPt*cont*width;
1927 jetPtUncorr[i] = ptUncorr;
1928 jetPtCorr[i] = ptCorr;
1931 // calc dXi from dN/dpt distribution :
1932 // sum over track pt for raw and corrected FF is equivalent to raw/corrected jet pt
1934 for(Int_t i=0; i<fNJetPtSlices; i++){
1936 Float_t jetPtLoLim = fJetPtSlices->At(i);
1937 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1939 Double_t meanJetPt = 0.5*(jetPtUpLim+jetPtLoLim);
1941 Double_t ptUncorr = jetPtUncorr[i];
1942 Double_t ptCorr = jetPtCorr[i];
1944 Double_t dXi = TMath::Log(ptCorr/ptUncorr);
1946 Printf(" calc xi shift: jet pt slice %d, mean jet pt %f, ptUncorr %f, ptCorr %f, ratio corr/uncorr %f, dXi %f "
1947 ,i,meanJetPt,ptUncorr,ptCorr,ptCorr/ptUncorr,dXi);
1952 // book & fill new dN/dxi histos
1954 for(Int_t i=0; i<fNJetPtSlices; i++){
1956 TH1F* histXi = fCorrFF[corrLevel]->GetXi(i);
1958 Double_t dXi = deltaXi[i];
1960 Int_t nBins = histXi->GetNbinsX();
1961 const Double_t* binsVec = histXi->GetXaxis()->GetXbins()->GetArray();
1962 Float_t binsVecNew[nBins+1];
1964 TString strName = histXi->GetName();
1965 strName.Append("_shift");
1966 TString strTit = histXi->GetTitle();
1970 // create shifted histo ...
1972 if(binsVec){ // binsVec only neq NULL if histo was rebinned before
1974 for(Int_t bin=0; bin<nBins+1; bin++) binsVecNew[bin] = binsVec[bin] + dXi;
1975 hXiCorr = new TH1F(strName,strTit,nBins,binsVecNew);
1977 else{ // uniform bin size
1979 Double_t xMin = histXi->GetXaxis()->GetXmin();
1980 Double_t xMax = histXi->GetXaxis()->GetXmax();
1985 hXiCorr = new TH1F(strName,strTit,nBins,xMin,xMax);
1990 for(Int_t bin=1; bin<nBins+1; bin++){
1991 Double_t cont = histXi->GetBinContent(bin);
1992 Double_t err = histXi->GetBinError(bin);
1994 hXiCorr->SetBinContent(bin,cont);
1995 hXiCorr->SetBinError(bin,err);
1998 new(fh1FFXiShift[i]) TH1F(*hXiCorr);
2002 delete[] jetPtUncorr;
2007 //_____________________________________________________
2008 void AliFragmentationFunctionCorrections::SubtractBgr()
2010 // subtract bgr distribution from FF
2011 // the current corr level is used for both
2013 AddCorrectionLevel("bgrSub");
2015 for(Int_t i=0; i<fNJetPtSlices; i++){
2017 TH1F* histPt = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i);
2018 TH1F* histZ = fCorrFF[fNCorrectionLevels-2]->GetZ(i);
2019 TH1F* histXi = fCorrFF[fNCorrectionLevels-2]->GetXi(i);
2021 TH1F* histPtBgr = fCorrBgr[fNCorrectionLevelsBgr-1]->GetTrackPt(i);
2022 TH1F* histZBgr = fCorrBgr[fNCorrectionLevelsBgr-1]->GetZ(i);
2023 TH1F* histXiBgr = fCorrBgr[fNCorrectionLevelsBgr-1]->GetXi(i);
2025 TString histNamePt = histPt->GetName();
2026 TString histNameZ = histZ->GetName();
2027 TString histNameXi = histXi->GetName();
2030 TH1F* hFFTrackPtBgrSub = (TH1F*) histPt->Clone(histNamePt);
2031 hFFTrackPtBgrSub->Add(histPtBgr,-1);
2033 TH1F* hFFZBgrSub = (TH1F*) histZ->Clone(histNameZ);
2034 hFFZBgrSub->Add(histZBgr,-1);
2036 TH1F* hFFXiBgrSub = (TH1F*) histXi->Clone(histNameXi);
2037 hFFXiBgrSub->Add(histXiBgr,-1);
2039 fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hFFTrackPtBgrSub,hFFZBgrSub,hFFXiBgrSub);
2043 //________________________________________________________________________________________________________________
2044 void AliFragmentationFunctionCorrections::WriteSingleTrackEff(TString strInfile, TString strID, TString strOutfile,
2045 Bool_t updateOutfile, TString strOutDir,TString strPostfix)
2047 // read task ouput from MC and write single track eff - standard dir/list
2049 TString strdir = "PWGJE_FragmentationFunction_" + strID;
2050 TString strlist = "fracfunc_" + strID;
2052 WriteSingleTrackEff(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir,strPostfix);
2055 //___________________________________________________________________________________________________________________________________
2056 void AliFragmentationFunctionCorrections::WriteSingleTrackEff(TString strInfile, TString strdir, TString strlist,
2057 TString strOutfile, Bool_t updateOutfile, TString strOutDir,TString strPostfix)
2059 // read task output from MC and write single track eff as function of pt, eta, phi
2060 // argument strlist optional - read from directory strdir if not specified
2061 // write eff to file stroutfile - by default only 'update' file (don't overwrite)
2064 TH1D* hdNdptTracksMCPrimGen;
2065 TH2D* hdNdetadphiTracksMCPrimGen;
2067 TH1D* hdNdptTracksMCPrimRec;
2068 TH2D* hdNdetadphiTracksMCPrimRec;
2071 TFile f(strInfile,"READ");
2074 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2078 if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: open task ouput file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
2080 if(strdir && strdir.Length()){
2081 if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: change dir to %s",(char*)__FILE__,__LINE__,strdir.Data());
2082 gDirectory->cd(strdir);
2087 if(strlist && strlist.Length()){
2089 if(!(list = (TList*) gDirectory->Get(strlist))){
2090 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2096 TString hnamePtRecEffGen = "fh1TrackQAPtRecEffGen";
2097 if(strPostfix.Length()) hnamePtRecEffGen.Form("fh1TrackQAPtRecEffGen_%s",strPostfix.Data());
2099 TString hnamePtRecEffRec = "fh1TrackQAPtRecEffRec";
2100 if(strPostfix.Length()) hnamePtRecEffRec.Form("fh1TrackQAPtRecEffRec_%s",strPostfix.Data());
2102 TString hnameEtaPhiRecEffGen = "fh2TrackQAEtaPhiRecEffGen";
2103 if(strPostfix.Length()) hnameEtaPhiRecEffGen.Form("fh2TrackQAEtaPhiRecEffGen_%s",strPostfix.Data());
2105 TString hnameEtaPhiRecEffRec = "fh2TrackQAEtaPhiRecEffRec";
2106 if(strPostfix.Length()) hnameEtaPhiRecEffRec.Form("fh2TrackQAEtaPhiRecEffRec_%s",strPostfix.Data());
2110 hdNdptTracksMCPrimGen = (TH1D*) list->FindObject(hnamePtRecEffGen);
2111 hdNdetadphiTracksMCPrimGen = (TH2D*) list->FindObject(hnameEtaPhiRecEffGen);
2113 hdNdptTracksMCPrimRec = (TH1D*) list->FindObject(hnamePtRecEffRec);
2114 hdNdetadphiTracksMCPrimRec = (TH2D*) list->FindObject(hnameEtaPhiRecEffRec);
2117 hdNdptTracksMCPrimGen = (TH1D*) gDirectory->Get(hnamePtRecEffGen);
2118 hdNdetadphiTracksMCPrimGen = (TH2D*) gDirectory->Get(hnameEtaPhiRecEffGen);
2120 hdNdptTracksMCPrimRec = (TH1D*) gDirectory->Get(hnamePtRecEffRec);
2121 hdNdetadphiTracksMCPrimRec = (TH2D*) gDirectory->Get(hnameEtaPhiRecEffRec);
2124 hdNdptTracksMCPrimGen->SetDirectory(0);
2125 hdNdetadphiTracksMCPrimGen->SetDirectory(0);
2126 hdNdptTracksMCPrimRec->SetDirectory(0);
2127 hdNdetadphiTracksMCPrimRec->SetDirectory(0);
2131 // projections: dN/deta, dN/dphi
2133 TH1D* hdNdetaTracksMCPrimGen = (TH1D*) hdNdetadphiTracksMCPrimGen->ProjectionX("hdNdetaTracksMcPrimGen");
2134 TH1D* hdNdphiTracksMCPrimGen = (TH1D*) hdNdetadphiTracksMCPrimGen->ProjectionY("hdNdphiTracksMcPrimGen");
2136 TH1D* hdNdetaTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionX("hdNdetaTracksMcPrimRec");
2137 TH1D* hdNdphiTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionY("hdNdphiTracksMcPrimRec");
2141 TString strNamePtGen = "hTrackPtGenPrim";
2142 TString strNamePtRec = "hTrackPtRecPrim";
2144 if(fNHistoBinsSinglePt) hdNdptTracksMCPrimGen = (TH1D*) hdNdptTracksMCPrimGen->Rebin(fNHistoBinsSinglePt,strNamePtGen,fHistoBinsSinglePt->GetArray());
2145 if(fNHistoBinsSinglePt) hdNdptTracksMCPrimRec = (TH1D*) hdNdptTracksMCPrimRec->Rebin(fNHistoBinsSinglePt,strNamePtRec,fHistoBinsSinglePt->GetArray());
2147 // efficiency: divide
2149 TString hNameTrackEffPt = "hSingleTrackEffPt";
2150 if(strPostfix.Length()) hNameTrackEffPt.Form("hSingleTrackEffPt_%s",strPostfix.Data());
2152 TString hNameTrackEffEta = "hSingleTrackEffEta";
2153 if(strPostfix.Length()) hNameTrackEffEta.Form("hSingleTrackEffEta_%s",strPostfix.Data());
2155 TString hNameTrackEffPhi = "hSingleTrackEffPhi";
2156 if(strPostfix.Length()) hNameTrackEffPhi.Form("hSingleTrackEffPhi_%s",strPostfix.Data());
2159 TH1F* hSingleTrackEffPt = (TH1F*) hdNdptTracksMCPrimRec->Clone(hNameTrackEffPt);
2160 hSingleTrackEffPt->Divide(hdNdptTracksMCPrimRec,hdNdptTracksMCPrimGen,1,1,"B"); // binominal errors
2162 TH1F* hSingleTrackEffEta = (TH1F*) hdNdetaTracksMCPrimRec->Clone(hNameTrackEffEta);
2163 hSingleTrackEffEta->Divide(hdNdetaTracksMCPrimRec,hdNdetaTracksMCPrimGen,1,1,"B"); // binominal errors
2165 TH1F* hSingleTrackEffPhi = (TH1F*) hdNdphiTracksMCPrimRec->Clone(hNameTrackEffPhi);
2166 hSingleTrackEffPhi->Divide(hdNdphiTracksMCPrimRec,hdNdphiTracksMCPrimGen,1,1,"B"); // binominal errors
2169 TString outfileOption = "RECREATE";
2170 if(updateOutfile) outfileOption = "UPDATE";
2172 TFile out(strOutfile,outfileOption);
2175 Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2179 if(fDebug>0) Printf("%s:%d -- write efficiency to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2181 if(strOutDir && strOutDir.Length()){
2184 if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd();
2186 dir = out.mkdir(strOutDir);
2191 hSingleTrackEffPt->Write();
2192 hSingleTrackEffEta->Write();
2193 hSingleTrackEffPhi->Write();
2198 //________________________________________________________________________________________________________________
2199 void AliFragmentationFunctionCorrections::WriteSingleTrackSecCorr(TString strInfile, TString strID, TString strOutfile,
2200 Bool_t updateOutfile, TString strOutDir)
2202 // read task ouput from MC and write single track eff - standard dir/list
2204 TString strdir = "PWGJE_FragmentationFunction_" + strID;
2205 TString strlist = "fracfunc_" + strID;
2207 WriteSingleTrackSecCorr(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir);
2210 //___________________________________________________________________________________________________________________________________
2211 void AliFragmentationFunctionCorrections::WriteSingleTrackSecCorr(TString strInfile, TString strdir, TString strlist,
2212 TString strOutfile, Bool_t updateOutfile, TString strOutDir)
2214 // read task output from MC and write single track secondaries contamination correction as function of pt, eta, phi
2215 // argument strlist optional - read from directory strdir if not specified
2216 // write corr factor to file stroutfile - by default only 'update' file (don't overwrite)
2218 TH1D* hdNdptTracksMCPrimRec;
2219 TH2D* hdNdetadphiTracksMCPrimRec;
2221 TH1D* hdNdptTracksMCSecRec;
2222 TH2D* hdNdetadphiTracksMCSecRec;
2224 TFile f(strInfile,"READ");
2227 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2231 if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
2233 if(strdir && strdir.Length()) gDirectory->cd(strdir);
2237 if(strlist && strlist.Length()){
2239 if(!(list = (TList*) gDirectory->Get(strlist))){
2240 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2247 hdNdptTracksMCPrimRec = (TH1D*) list->FindObject("fh1TrackQAPtRecEffGen");
2248 hdNdetadphiTracksMCPrimRec = (TH2D*) list->FindObject("fh2TrackQAEtaPhiRecEffGen");
2250 hdNdptTracksMCSecRec = (TH1D*) list->FindObject("fh1TrackQAPtSecRec");
2251 hdNdetadphiTracksMCSecRec = (TH2D*) list->FindObject("fh2TrackQAEtaPhiSecRec");
2254 hdNdptTracksMCPrimRec = (TH1D*) gDirectory->Get("fh1TrackQAPtRecEffGen");
2255 hdNdetadphiTracksMCPrimRec = (TH2D*) gDirectory->Get("fh2TrackQAEtaPhiRecEffGen");
2257 hdNdptTracksMCSecRec = (TH1D*) gDirectory->Get("fh1TrackQAPtSecRec");
2258 hdNdetadphiTracksMCSecRec = (TH2D*) gDirectory->Get("fh2TrackQAEtaPhiSecRec");
2261 hdNdptTracksMCPrimRec->SetDirectory(0);
2262 hdNdetadphiTracksMCPrimRec->SetDirectory(0);
2263 hdNdptTracksMCSecRec->SetDirectory(0);
2264 hdNdetadphiTracksMCSecRec->SetDirectory(0);
2268 // projections: dN/deta, dN/dphi
2270 TH1D* hdNdetaTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionX("hdNdetaTracksMcPrimRec");
2271 TH1D* hdNdphiTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionY("hdNdphiTracksMcPrimRec");
2273 TH1D* hdNdetaTracksMCSecRec = (TH1D*) hdNdetadphiTracksMCSecRec->ProjectionX("hdNdetaTracksMcSecRec");
2274 TH1D* hdNdphiTracksMCSecRec = (TH1D*) hdNdetadphiTracksMCSecRec->ProjectionY("hdNdphiTracksMcSecRec");
2279 TString strNamePtPrim = "hTrackPtPrim";
2280 TString strNamePtSec = "hTrackPtSec";
2282 if(fNHistoBinsSinglePt) hdNdptTracksMCPrimRec = (TH1D*) hdNdptTracksMCPrimRec->Rebin(fNHistoBinsSinglePt,strNamePtPrim,fHistoBinsSinglePt->GetArray());
2283 if(fNHistoBinsSinglePt) hdNdptTracksMCSecRec = (TH1D*) hdNdptTracksMCSecRec->Rebin(fNHistoBinsSinglePt,strNamePtSec,fHistoBinsSinglePt->GetArray());
2286 // secondary correction factor: divide prim/(prim+sec)
2288 TH1F* hSingleTrackSecCorrPt = (TH1F*) hdNdptTracksMCSecRec->Clone("hSingleTrackSecCorrPt");
2289 TH1F* hSumPrimSecPt = (TH1F*) hdNdptTracksMCSecRec->Clone("hSumPrimSecPt");
2290 hSumPrimSecPt->Add(hdNdptTracksMCPrimRec);
2291 hSingleTrackSecCorrPt->Divide(hdNdptTracksMCPrimRec,hSumPrimSecPt,1,1,"B"); // binominal errors
2293 TH1F* hSingleTrackSecCorrEta = (TH1F*) hdNdetaTracksMCSecRec->Clone("hSingleTrackSecCorrEta");
2294 TH1F* hSumPrimSecEta = (TH1F*) hdNdetaTracksMCSecRec->Clone("hSumPrimSecEta");
2295 hSumPrimSecEta->Add(hdNdetaTracksMCPrimRec);
2296 hSingleTrackSecCorrEta->Divide(hdNdetaTracksMCPrimRec,hSumPrimSecEta,1,1,"B"); // binominal errors
2298 TH1F* hSingleTrackSecCorrPhi = (TH1F*) hdNdphiTracksMCSecRec->Clone("hSingleTrackSecCorrPhi");
2299 TH1F* hSumPrimSecPhi = (TH1F*) hdNdphiTracksMCSecRec->Clone("hSumPrimSecPhi");
2300 hSumPrimSecPhi->Add(hdNdphiTracksMCPrimRec);
2301 hSingleTrackSecCorrPhi->Divide(hdNdphiTracksMCPrimRec,hSumPrimSecPhi,1,1,"B"); // binominal errors
2304 TString outfileOption = "RECREATE";
2305 if(updateOutfile) outfileOption = "UPDATE";
2307 TFile out(strOutfile,outfileOption);
2310 Printf("%s:%d -- error opening secCorr output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2314 if(fDebug>0) Printf("%s:%d -- write secCorr to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2316 if(strOutDir && strOutDir.Length()){
2319 if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd();
2321 dir = out.mkdir(strOutDir);
2326 hdNdptTracksMCSecRec->Write();
2327 hdNdetaTracksMCSecRec->Write();
2328 hdNdphiTracksMCSecRec->Write();
2330 hSingleTrackSecCorrPt->Write();
2331 hSingleTrackSecCorrEta->Write();
2332 hSingleTrackSecCorrPhi->Write();
2337 //________________________________________________________________________________________________________________
2338 void AliFragmentationFunctionCorrections::WriteSingleResponse(TString strInfile, TString strID, TString strOutfile,
2339 Bool_t updateOutfile, TString strOutDir)
2341 // read task ouput from MC and write single track eff - standard dir/list
2343 TString strdir = "PWGJE_FragmentationFunction_" + strID;
2344 TString strlist = "fracfunc_" + strID;
2346 WriteSingleResponse(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir);
2350 //_____________________________________________________________________________________________________________________________________
2351 void AliFragmentationFunctionCorrections::WriteSingleResponse(TString strInfile, TString strdir, TString strlist,
2352 TString strOutfile, Bool_t updateOutfile, TString strOutDir)
2354 // read 2d THnSparse response matrices in pt from file
2356 // write to strOutfile
2358 THnSparse* hnResponseSinglePt;
2359 TH2F* h2ResponseSinglePt;
2361 TFile f(strInfile,"READ");
2364 Printf("%s:%d -- error opening file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2368 if(fDebug>0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
2370 if(strdir && strdir.Length()) gDirectory->cd(strdir);
2374 if(strlist && strlist.Length()){
2376 if(!(list = (TList*) gDirectory->Get(strlist))){
2377 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2382 if(list) hnResponseSinglePt = (THnSparse*) list->FindObject("fhnResponseSinglePt");
2383 else hnResponseSinglePt = (THnSparse*) gDirectory->Get("fhnResponseSinglePt");
2386 if(!hnResponseSinglePt){
2387 Printf("%s:%d -- error retrieving response matrix fhnResponseSinglePt",(char*)__FILE__,__LINE__);
2395 h2ResponseSinglePt = (TH2F*) hnResponseSinglePt->Projection(1,0);// note convention: yDim,xDim
2396 h2ResponseSinglePt->SetNameTitle("h2ResponseSinglePt","");
2401 TString outfileOption = "RECREATE";
2402 if(updateOutfile) outfileOption = "UPDATE";
2404 TFile out(strOutfile,outfileOption);
2407 Printf("%s:%d -- error opening response matrix output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2411 if(fDebug>0) Printf("%s:%d -- write response matrices to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2413 if(strOutDir && strOutDir.Length()){
2416 if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd();
2418 dir = out.mkdir(strOutDir);
2423 hnResponseSinglePt->Write();
2424 h2ResponseSinglePt->Write();
2429 //________________________________________________________________________________________________________________
2430 void AliFragmentationFunctionCorrections::WriteJetTrackEff(TString strInfile, TString strID, TString strOutfile,
2431 Bool_t updateOutfile)
2433 // read task ouput from MC and write single track eff - standard dir/list
2435 TString strdir = "PWGJE_FragmentationFunction_" + strID;
2436 TString strlist = "fracfunc_" + strID;
2438 WriteJetTrackEff(strInfile,strdir,strlist,strOutfile,updateOutfile);
2441 //___________________________________________________________________________________________________________________________________
2442 void AliFragmentationFunctionCorrections::WriteJetTrackEff(TString strInfile, TString strdir, TString strlist,
2443 TString strOutfile, Bool_t updateOutfile)
2445 // read task output from MC and write track eff in jet pt slices as function of pt, z, xi
2446 // argument strlist optional - read from directory strdir if not specified
2447 // write eff to file strOutfile - by default only 'update' file (don't overwrite)
2449 TH1F* hEffPt[fNJetPtSlices];
2450 TH1F* hEffXi[fNJetPtSlices];
2451 TH1F* hEffZ[fNJetPtSlices];
2453 TH1F* hdNdptTracksMCPrimGen[fNJetPtSlices];
2454 TH1F* hdNdxiMCPrimGen[fNJetPtSlices];
2455 TH1F* hdNdzMCPrimGen[fNJetPtSlices];
2457 TH1F* hdNdptTracksMCPrimRec[fNJetPtSlices];
2458 TH1F* hdNdxiMCPrimRec[fNJetPtSlices];
2459 TH1F* hdNdzMCPrimRec[fNJetPtSlices];
2462 TH1F* fh1FFJetPtRecEffGen;
2464 TH2F* fh2FFTrackPtRecEffGen;
2465 TH2F* fh2FFZRecEffGen;
2466 TH2F* fh2FFXiRecEffGen;
2468 TH2F* fh2FFTrackPtRecEffRec;
2469 TH2F* fh2FFZRecEffRec;
2470 TH2F* fh2FFXiRecEffRec;
2473 TFile f(strInfile,"READ");
2476 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2480 if(fDebug>0) Printf("%s:%d -- writeJetTrackEff: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
2482 if(strdir && strdir.Length()) gDirectory->cd(strdir);
2486 if(strlist && strlist.Length()){
2488 if(!(list = (TList*) gDirectory->Get(strlist))){
2489 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2495 fh1FFJetPtRecEffGen = (TH1F*) list->FindObject("fh1FFJetPtRecEffGen");
2497 fh2FFTrackPtRecEffGen = (TH2F*) list->FindObject("fh2FFTrackPtRecEffGen");
2498 fh2FFZRecEffGen = (TH2F*) list->FindObject("fh2FFZRecEffGen");
2499 fh2FFXiRecEffGen = (TH2F*) list->FindObject("fh2FFXiRecEffGen");
2501 fh2FFTrackPtRecEffRec = (TH2F*) list->FindObject("fh2FFTrackPtRecEffRec");
2502 fh2FFZRecEffRec = (TH2F*) list->FindObject("fh2FFZRecEffRec");
2503 fh2FFXiRecEffRec = (TH2F*) list->FindObject("fh2FFXiRecEffRec");
2506 fh1FFJetPtRecEffGen = (TH1F*) gDirectory->FindObject("fh1FFJetPtRecEffGen");
2508 fh2FFTrackPtRecEffGen = (TH2F*) gDirectory->FindObject("fh2FFTrackPtRecEffGen");
2509 fh2FFZRecEffGen = (TH2F*) gDirectory->FindObject("fh2FFZRecEffGen");
2510 fh2FFXiRecEffGen = (TH2F*) gDirectory->FindObject("fh2FFXiRecEffGen");
2512 fh2FFTrackPtRecEffRec = (TH2F*) gDirectory->FindObject("fh2FFTrackPtRecEffRec");
2513 fh2FFZRecEffRec = (TH2F*) gDirectory->FindObject("fh2FFZRecEffRec");
2514 fh2FFXiRecEffRec = (TH2F*) gDirectory->FindObject("fh2FFXiRecEffRec");
2517 fh1FFJetPtRecEffGen->SetDirectory(0);
2519 fh2FFTrackPtRecEffGen->SetDirectory(0);
2520 fh2FFZRecEffGen->SetDirectory(0);
2521 fh2FFXiRecEffGen->SetDirectory(0);
2523 fh2FFTrackPtRecEffRec->SetDirectory(0);
2524 fh2FFZRecEffRec->SetDirectory(0);
2525 fh2FFXiRecEffRec->SetDirectory(0);
2530 // projections: FF for generated and reconstructed primaries
2532 for(Int_t i=0; i<fNJetPtSlices; i++){
2534 Float_t jetPtLoLim = fJetPtSlices->At(i);
2535 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
2537 Int_t binLo = static_cast<Int_t>(fh2FFTrackPtRecEffGen->GetXaxis()->FindBin(jetPtLoLim));
2538 Int_t binUp = static_cast<Int_t>(fh2FFTrackPtRecEffGen->GetXaxis()->FindBin(jetPtUpLim))-1;
2540 if(binUp > fh2FFTrackPtRecEffGen->GetNbinsX()){
2541 Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
2545 TString strNameFFPtGen(Form("fh1FFTrackPtGenPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2546 TString strNameFFZGen(Form("fh1FFZGenPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2547 TString strNameFFXiGen(Form("fh1FFXiGenPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2549 TString strNameFFPtRec(Form("fh1FFTrackPtRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2550 TString strNameFFZRec(Form("fh1FFZRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2551 TString strNameFFXiRec(Form("fh1FFXiRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2554 // appendix 'unbinned' to avoid histos with same name after rebinning
2556 hdNdptTracksMCPrimGen[i] = (TH1F*) fh2FFTrackPtRecEffGen->ProjectionY(strNameFFPtGen+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
2557 hdNdzMCPrimGen[i] = (TH1F*) fh2FFZRecEffGen->ProjectionY(strNameFFZGen+"_unBinned",binLo,binUp,"o");
2558 hdNdxiMCPrimGen[i] = (TH1F*) fh2FFXiRecEffGen->ProjectionY(strNameFFXiGen+"_unBinned",binLo,binUp,"o");
2560 hdNdptTracksMCPrimRec[i] = (TH1F*) fh2FFTrackPtRecEffRec->ProjectionY(strNameFFPtRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
2561 hdNdzMCPrimRec[i] = (TH1F*) fh2FFZRecEffRec->ProjectionY(strNameFFZRec+"_unBinned",binLo,binUp,"o");
2562 hdNdxiMCPrimRec[i] = (TH1F*) fh2FFXiRecEffRec->ProjectionY(strNameFFXiRec+"_unBinned",binLo,binUp,"o");
2566 if(fNHistoBinsPt[i]) hdNdptTracksMCPrimGen[i] = (TH1F*) hdNdptTracksMCPrimGen[i]->Rebin(fNHistoBinsPt[i],strNameFFPtGen,fHistoBinsPt[i]->GetArray());
2567 if(fNHistoBinsZ[i]) hdNdzMCPrimGen[i] = (TH1F*) hdNdzMCPrimGen[i]->Rebin(fNHistoBinsZ[i],strNameFFZGen,fHistoBinsZ[i]->GetArray());
2568 if(fNHistoBinsXi[i]) hdNdxiMCPrimGen[i] = (TH1F*) hdNdxiMCPrimGen[i]->Rebin(fNHistoBinsXi[i],strNameFFXiGen,fHistoBinsXi[i]->GetArray());
2570 if(fNHistoBinsPt[i]) hdNdptTracksMCPrimRec[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtRec,fHistoBinsPt[i]->GetArray());
2571 if(fNHistoBinsZ[i]) hdNdzMCPrimRec[i] = (TH1F*) hdNdzMCPrimRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZRec,fHistoBinsZ[i]->GetArray());
2572 if(fNHistoBinsXi[i]) hdNdxiMCPrimRec[i] = (TH1F*) hdNdxiMCPrimRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiRec,fHistoBinsXi[i]->GetArray());
2574 hdNdptTracksMCPrimGen[i]->SetNameTitle(strNameFFPtGen,"");
2575 hdNdzMCPrimGen[i]->SetNameTitle(strNameFFZGen,"");
2576 hdNdxiMCPrimGen[i]->SetNameTitle(strNameFFXiGen,"");
2578 hdNdptTracksMCPrimRec[i]->SetNameTitle(strNameFFPtRec,"");
2579 hdNdzMCPrimRec[i]->SetNameTitle(strNameFFZRec,"");
2580 hdNdxiMCPrimRec[i]->SetNameTitle(strNameFFXiRec,"");
2584 Double_t nJetsBin = fh1FFJetPtRecEffGen->Integral(binLo,binUp);
2586 NormalizeTH1(hdNdptTracksMCPrimGen[i],nJetsBin);
2587 NormalizeTH1(hdNdzMCPrimGen[i],nJetsBin);
2588 NormalizeTH1(hdNdxiMCPrimGen[i],nJetsBin);
2590 NormalizeTH1(hdNdptTracksMCPrimRec[i],nJetsBin);
2591 NormalizeTH1(hdNdzMCPrimRec[i],nJetsBin);
2592 NormalizeTH1(hdNdxiMCPrimRec[i],nJetsBin);
2594 // divide rec/gen : efficiency
2596 TString strNameEffPt(Form("hEffPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2597 TString strNameEffZ(Form("hEffZ_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2598 TString strNameEffXi(Form("hEffXi_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2600 hEffPt[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Clone(strNameEffPt);
2601 hEffPt[i]->Divide(hdNdptTracksMCPrimRec[i],hdNdptTracksMCPrimGen[i],1,1,"B"); // binominal errors
2603 hEffXi[i] = (TH1F*) hdNdxiMCPrimRec[i]->Clone(strNameEffXi);
2604 hEffXi[i]->Divide(hdNdxiMCPrimRec[i],hdNdxiMCPrimGen[i],1,1,"B"); // binominal errors
2606 hEffZ[i] = (TH1F*) hdNdzMCPrimRec[i]->Clone(strNameEffZ);
2607 hEffZ[i]->Divide(hdNdzMCPrimRec[i],hdNdzMCPrimGen[i],1,1,"B"); // binominal errors
2612 TString outfileOption = "RECREATE";
2613 if(updateOutfile) outfileOption = "UPDATE";
2615 TFile out(strOutfile,outfileOption);
2618 Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2622 if(fDebug>0) Printf("%s:%d -- write efficiency to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2624 // if(strdir && strdir.Length()){
2625 // TDirectory* dir = out.mkdir(strdir);
2629 for(Int_t i=0; i<fNJetPtSlices; i++){
2631 hdNdptTracksMCPrimGen[i]->Write();
2632 hdNdxiMCPrimGen[i]->Write();
2633 hdNdzMCPrimGen[i]->Write();
2635 hdNdptTracksMCPrimRec[i]->Write();
2636 hdNdxiMCPrimRec[i]->Write();
2637 hdNdzMCPrimRec[i]->Write();
2648 //________________________________________________________________________________________________________________
2649 void AliFragmentationFunctionCorrections::WriteJetSecCorr(TString strInfile, TString strID, TString strOutfile,
2650 Bool_t updateOutfile)
2652 // read task ouput from MC and write secondary correction - standard dir/list
2654 TString strdir = "PWGJE_FragmentationFunction_" + strID;
2655 TString strlist = "fracfunc_" + strID;
2657 WriteJetSecCorr(strInfile,strdir,strlist,strOutfile,updateOutfile);
2660 //___________________________________________________________________________________________________________________________________
2661 void AliFragmentationFunctionCorrections::WriteJetSecCorr(TString strInfile, TString strdir, TString strlist,
2662 TString strOutfile, Bool_t updateOutfile)
2664 // read task output from MC and write secondary correction in jet pt slices as function of pt, z, xi
2665 // argument strlist optional - read from directory strdir if not specified
2666 // write eff to file strOutfile - by default only 'update' file (don't overwrite)
2668 TH1F* hSecCorrPt[fNJetPtSlices];
2669 TH1F* hSecCorrXi[fNJetPtSlices];
2670 TH1F* hSecCorrZ[fNJetPtSlices];
2672 TH1F* hdNdptTracksMCPrimRec[fNJetPtSlices];
2673 TH1F* hdNdxiMCPrimRec[fNJetPtSlices];
2674 TH1F* hdNdzMCPrimRec[fNJetPtSlices];
2676 TH1F* hdNdptTracksMCSecRec[fNJetPtSlices];
2677 TH1F* hdNdxiMCSecRec[fNJetPtSlices];
2678 TH1F* hdNdzMCSecRec[fNJetPtSlices];
2680 TH1F* fh1FFJetPtRecEffGen;
2682 TH2F* fh2FFTrackPtRecEffRec;
2683 TH2F* fh2FFZRecEffRec;
2684 TH2F* fh2FFXiRecEffRec;
2686 TH2F* fh2FFTrackPtSecRec;
2688 TH2F* fh2FFXiSecRec;
2690 TFile f(strInfile,"READ");
2693 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2697 if(fDebug>0) Printf("%s:%d -- writeJetTrackEff: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
2699 if(strdir && strdir.Length()) gDirectory->cd(strdir);
2703 if(strlist && strlist.Length()){
2705 if(!(list = (TList*) gDirectory->Get(strlist))){
2706 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2712 fh1FFJetPtRecEffGen = (TH1F*) list->FindObject("fh1FFJetPtRecEffGen");
2714 fh2FFTrackPtRecEffRec = (TH2F*) list->FindObject("fh2FFTrackPtRecEffRec");
2715 fh2FFZRecEffRec = (TH2F*) list->FindObject("fh2FFZRecEffRec");
2716 fh2FFXiRecEffRec = (TH2F*) list->FindObject("fh2FFXiRecEffRec");
2718 fh2FFTrackPtSecRec = (TH2F*) list->FindObject("fh2FFTrackPtSecRec");
2719 fh2FFZSecRec = (TH2F*) list->FindObject("fh2FFZSecRec");
2720 fh2FFXiSecRec = (TH2F*) list->FindObject("fh2FFXiSecRec");
2723 fh1FFJetPtRecEffGen = (TH1F*) gDirectory->FindObject("fh1FFJetPtRecEffGen");
2725 fh2FFTrackPtRecEffRec = (TH2F*) gDirectory->FindObject("fh2FFTrackPtRecEffRec");
2726 fh2FFZRecEffRec = (TH2F*) gDirectory->FindObject("fh2FFZRecEffRec");
2727 fh2FFXiRecEffRec = (TH2F*) gDirectory->FindObject("fh2FFXiRecEffRec");
2729 fh2FFTrackPtSecRec = (TH2F*) gDirectory->FindObject("fh2FFTrackPtSecRec");
2730 fh2FFZSecRec = (TH2F*) gDirectory->FindObject("fh2FFZSecRec");
2731 fh2FFXiSecRec = (TH2F*) gDirectory->FindObject("fh2FFXiSecRec");
2734 fh1FFJetPtRecEffGen->SetDirectory(0);
2736 fh2FFTrackPtRecEffRec->SetDirectory(0);
2737 fh2FFZRecEffRec->SetDirectory(0);
2738 fh2FFXiRecEffRec->SetDirectory(0);
2740 fh2FFTrackPtSecRec->SetDirectory(0);
2741 fh2FFZSecRec->SetDirectory(0);
2742 fh2FFXiSecRec->SetDirectory(0);
2747 // projections: FF for generated and reconstructed primaries
2749 for(Int_t i=0; i<fNJetPtSlices; i++){
2751 Float_t jetPtLoLim = fJetPtSlices->At(i);
2752 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
2754 Int_t binLo = static_cast<Int_t>(fh2FFTrackPtRecEffRec->GetXaxis()->FindBin(jetPtLoLim));
2755 Int_t binUp = static_cast<Int_t>(fh2FFTrackPtRecEffRec->GetXaxis()->FindBin(jetPtUpLim))-1;
2757 if(binUp > fh2FFTrackPtRecEffRec->GetNbinsX()){
2758 Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
2762 TString strNameFFPtPrimRec(Form("fh1FFTrackPtRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2763 TString strNameFFZPrimRec(Form("fh1FFZRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2764 TString strNameFFXiPrimRec(Form("fh1FFXiRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2766 TString strNameFFPtSecRec(Form("fh1FFTrackPtRecSec_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2767 TString strNameFFZSecRec(Form("fh1FFZRecSec_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2768 TString strNameFFXiSecRec(Form("fh1FFXiRecSec_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2771 // appendix 'unbinned' to avoid histos with same name after rebinning
2773 hdNdptTracksMCPrimRec[i] = (TH1F*) fh2FFTrackPtRecEffRec->ProjectionY(strNameFFPtPrimRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
2774 hdNdzMCPrimRec[i] = (TH1F*) fh2FFZRecEffRec->ProjectionY(strNameFFZPrimRec+"_unBinned",binLo,binUp,"o");
2775 hdNdxiMCPrimRec[i] = (TH1F*) fh2FFXiRecEffRec->ProjectionY(strNameFFXiPrimRec+"_unBinned",binLo,binUp,"o");
2777 hdNdptTracksMCSecRec[i] = (TH1F*) fh2FFTrackPtSecRec->ProjectionY(strNameFFPtSecRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
2778 hdNdzMCSecRec[i] = (TH1F*) fh2FFZSecRec->ProjectionY(strNameFFZSecRec+"_unBinned",binLo,binUp,"o");
2779 hdNdxiMCSecRec[i] = (TH1F*) fh2FFXiSecRec->ProjectionY(strNameFFXiSecRec+"_unBinned",binLo,binUp,"o");
2783 if(fNHistoBinsPt[i]) hdNdptTracksMCPrimRec[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtPrimRec,fHistoBinsPt[i]->GetArray());
2784 if(fNHistoBinsZ[i]) hdNdzMCPrimRec[i] = (TH1F*) hdNdzMCPrimRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZPrimRec,fHistoBinsZ[i]->GetArray());
2785 if(fNHistoBinsXi[i]) hdNdxiMCPrimRec[i] = (TH1F*) hdNdxiMCPrimRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiPrimRec,fHistoBinsXi[i]->GetArray());
2787 if(fNHistoBinsPt[i]) hdNdptTracksMCSecRec[i] = (TH1F*) hdNdptTracksMCSecRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtSecRec,fHistoBinsPt[i]->GetArray());
2788 if(fNHistoBinsZ[i]) hdNdzMCSecRec[i] = (TH1F*) hdNdzMCSecRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZSecRec,fHistoBinsZ[i]->GetArray());
2789 if(fNHistoBinsXi[i]) hdNdxiMCSecRec[i] = (TH1F*) hdNdxiMCSecRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiSecRec,fHistoBinsXi[i]->GetArray());
2791 hdNdptTracksMCPrimRec[i]->SetNameTitle(strNameFFPtPrimRec,"");
2792 hdNdzMCPrimRec[i]->SetNameTitle(strNameFFZPrimRec,"");
2793 hdNdxiMCPrimRec[i]->SetNameTitle(strNameFFXiPrimRec,"");
2795 hdNdptTracksMCSecRec[i]->SetNameTitle(strNameFFPtSecRec,"");
2796 hdNdzMCSecRec[i]->SetNameTitle(strNameFFZSecRec,"");
2797 hdNdxiMCSecRec[i]->SetNameTitle(strNameFFXiSecRec,"");
2801 Double_t nJetsBin = fh1FFJetPtRecEffGen->Integral(binLo,binUp);
2803 NormalizeTH1(hdNdptTracksMCPrimRec[i],nJetsBin);
2804 NormalizeTH1(hdNdzMCPrimRec[i],nJetsBin);
2805 NormalizeTH1(hdNdxiMCPrimRec[i],nJetsBin);
2807 NormalizeTH1(hdNdptTracksMCSecRec[i],nJetsBin);
2808 NormalizeTH1(hdNdzMCSecRec[i],nJetsBin);
2809 NormalizeTH1(hdNdxiMCSecRec[i],nJetsBin);
2811 // divide rec/gen : efficiency
2813 TString strNameSecCorrPt(Form("hSecCorrPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2814 TString strNameSecCorrZ(Form("hSecCorrZ_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2815 TString strNameSecCorrXi(Form("hSecCorrXi_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2817 hSecCorrPt[i] = (TH1F*) hdNdptTracksMCSecRec[i]->Clone(strNameSecCorrPt);
2818 TH1F* hSumPrimSecPt = (TH1F*) hdNdptTracksMCSecRec[i]->Clone("hSumPrimSecPt");
2819 hSumPrimSecPt->Add(hdNdptTracksMCPrimRec[i]);
2820 hSecCorrPt[i]->Divide(hdNdptTracksMCPrimRec[i],hSumPrimSecPt,1,1,"B"); // binominal errors
2822 hSecCorrXi[i] = (TH1F*) hdNdxiMCSecRec[i]->Clone(strNameSecCorrXi);
2823 TH1F* hSumPrimSecXi = (TH1F*) hdNdxiMCSecRec[i]->Clone("hSumPrimSecXi");
2824 hSumPrimSecXi->Add(hdNdxiMCPrimRec[i]);
2825 hSecCorrXi[i]->Divide(hdNdxiMCPrimRec[i],hSumPrimSecXi,1,1,"B"); // binominal errors
2827 hSecCorrZ[i] = (TH1F*) hdNdzMCSecRec[i]->Clone(strNameSecCorrZ);
2828 TH1F* hSumPrimSecZ = (TH1F*) hdNdzMCSecRec[i]->Clone("hSumPrimSecZ");
2829 hSumPrimSecZ->Add(hdNdzMCPrimRec[i]);
2830 hSecCorrZ[i]->Divide(hdNdzMCPrimRec[i],hSumPrimSecZ,1,1,"B"); // binominal errors
2835 TString outfileOption = "RECREATE";
2836 if(updateOutfile) outfileOption = "UPDATE";
2838 TFile out(strOutfile,outfileOption);
2841 Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2845 if(fDebug>0) Printf("%s:%d -- write efficiency to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2847 for(Int_t i=0; i<fNJetPtSlices; i++){
2849 // hdNdptTracksMCSecRec[i]->Write();
2850 // hdNdxiMCSecRec[i]->Write();
2851 // hdNdzMCSecRec[i]->Write();
2853 hSecCorrPt[i]->Write();
2854 hSecCorrXi[i]->Write();
2855 hSecCorrZ[i]->Write();
2861 //________________________________________________________________________________________________________________
2862 void AliFragmentationFunctionCorrections::WriteJetResponse(TString strInfile, TString strID, TString strOutfile,
2863 Bool_t updateOutfile)
2865 // read task ouput from MC and write single track eff - standard dir/list
2867 TString strdir = "PWGJE_FragmentationFunction_" + strID;
2868 TString strlist = "fracfunc_" + strID;
2870 WriteJetResponse(strInfile,strdir,strlist,strOutfile,updateOutfile);
2873 //_____________________________________________________________________________________________________________________________________
2874 void AliFragmentationFunctionCorrections::WriteJetResponse(TString strInfile, TString strdir, TString strlist,
2875 TString strOutfile, Bool_t updateOutfile)
2877 // read 3d THnSparse response matrices in pt,z,xi vs jet pt from file
2878 // project THnSparse + TH2 in jet pt slices
2879 // write to strOutfile
2881 THnSparse* hn3ResponseJetPt;
2882 THnSparse* hn3ResponseJetZ;
2883 THnSparse* hn3ResponseJetXi;
2885 // 2D response matrices
2887 THnSparse* hnResponsePt[fNJetPtSlices];
2888 THnSparse* hnResponseZ[fNJetPtSlices];
2889 THnSparse* hnResponseXi[fNJetPtSlices];
2891 TH2F* h2ResponsePt[fNJetPtSlices];
2892 TH2F* h2ResponseZ[fNJetPtSlices];
2893 TH2F* h2ResponseXi[fNJetPtSlices];
2895 // 1D projections on gen pt / rec pt axes
2897 TH1F* h1FFPtRec[fNJetPtSlices];
2898 TH1F* h1FFZRec[fNJetPtSlices];
2899 TH1F* h1FFXiRec[fNJetPtSlices];
2901 TH1F* h1FFPtGen[fNJetPtSlices];
2902 TH1F* h1FFZGen[fNJetPtSlices];
2903 TH1F* h1FFXiGen[fNJetPtSlices];
2905 TH1F* h1RatioPt[fNJetPtSlices];
2906 TH1F* h1RatioZ[fNJetPtSlices];
2907 TH1F* h1RatioXi[fNJetPtSlices];
2911 TFile f(strInfile,"READ");
2914 Printf("%s:%d -- error opening file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2918 if(fDebug>0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
2920 if(strdir && strdir.Length()) gDirectory->cd(strdir);
2924 if(strlist && strlist.Length()){
2926 if(!(list = (TList*) gDirectory->Get(strlist))){
2927 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2933 hn3ResponseJetPt = (THnSparse*) list->FindObject("fhnResponseJetTrackPt");
2934 hn3ResponseJetZ = (THnSparse*) list->FindObject("fhnResponseJetZ");
2935 hn3ResponseJetXi = (THnSparse*) list->FindObject("fhnResponseJetXi");
2938 hn3ResponseJetPt = (THnSparse*) gDirectory->Get("fhnResponseJetTrackPt");
2939 hn3ResponseJetZ = (THnSparse*) gDirectory->Get("fhnResponseJetZ");
2940 hn3ResponseJetXi = (THnSparse*) gDirectory->Get("fhnResponseJetXi");
2944 if(!hn3ResponseJetPt){
2945 Printf("%s:%d -- error retrieving response matrix fhnResponseJetTrackPt",(char*)__FILE__,__LINE__);
2949 if(!hn3ResponseJetZ){
2950 Printf("%s:%d -- error retrieving response matrix fhnResponseJetZ",(char*)__FILE__,__LINE__);
2954 if(!hn3ResponseJetXi){
2955 Printf("%s:%d -- error retrieving response matrix fhnResponseJetXi",(char*)__FILE__,__LINE__);
2963 Int_t axisJetPtTHn3 = -1;
2964 Int_t axisGenPtTHn3 = -1;
2965 Int_t axisRecPtTHn3 = -1;
2967 for(Int_t i=0; i<hn3ResponseJetPt->GetNdimensions(); i++){
2969 TString title = hn3ResponseJetPt->GetAxis(i)->GetTitle();
2971 if(title.Contains("jet p_{T}")) axisJetPtTHn3 = i;
2972 if(title.Contains("gen p_{T}")) axisGenPtTHn3 = i;
2973 if(title.Contains("rec p_{T}")) axisRecPtTHn3 = i;
2976 if(axisJetPtTHn3 == -1){
2977 Printf("%s:%d -- error axisJetPtTHn3",(char*)__FILE__,__LINE__);
2981 if(axisGenPtTHn3 == -1){
2982 Printf("%s:%d -- error axisGenPtTHn3",(char*)__FILE__,__LINE__);
2986 if(axisRecPtTHn3 == -1){
2987 Printf("%s:%d -- error axisRecPtTHn3",(char*)__FILE__,__LINE__);
2991 for(Int_t i=0; i<fNJetPtSlices; i++){
2993 Float_t jetPtLoLim = fJetPtSlices->At(i);
2994 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
2996 Int_t binLo = static_cast<Int_t>(hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->FindBin(jetPtLoLim));
2997 Int_t binUp = static_cast<Int_t>(hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->FindBin(jetPtUpLim))-1;
2999 if(binUp > hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->GetNbins()){
3000 Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
3004 TString strNameRespPt(Form("hnResponsePt_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3005 TString strNameRespZ(Form("hnResponseZ_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3006 TString strNameRespXi(Form("hnResponseXi_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3008 TString strNameTH2RespPt(Form("h2ResponsePt_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3009 TString strNameTH2RespZ(Form("h2ResponseZ_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3010 TString strNameTH2RespXi(Form("h2ResponseXi_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3012 TString strNameRecPt(Form("h1FFTrackPtRecPrim_recPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3013 TString strNameRecZ(Form("h1FFZRecPrim_recPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3014 TString strNameRecXi(Form("h1FFXiRecPrim_recPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3016 TString strNameGenPt(Form("h1FFTrackPtRecPrim_genPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3017 TString strNameGenZ(Form("h1FFZRecPrim_genPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3018 TString strNameGenXi(Form("h1FFXiRecPrim_genPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3020 TString strNameRatioPt(Form("h1RatioTrackPtRecPrim_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3021 TString strNameRatioZ(Form("h1RatioZRecPrim_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3022 TString strNameRatioXi(Form("h1RatioXiRecPrim_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3025 // 2D projections in jet pt range
3027 hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->SetRange(binLo,binUp);
3028 hn3ResponseJetZ->GetAxis(axisJetPtTHn3)->SetRange(binLo,binUp);
3029 hn3ResponseJetXi->GetAxis(axisJetPtTHn3)->SetRange(binLo,binUp);
3031 Int_t axesProj[2] = {axisRecPtTHn3,axisGenPtTHn3};
3033 hnResponsePt[i] = hn3ResponseJetPt->Projection(2,axesProj);
3034 hnResponseZ[i] = hn3ResponseJetZ->Projection(2,axesProj);
3035 hnResponseXi[i] = hn3ResponseJetXi->Projection(2,axesProj);
3037 hnResponsePt[i]->SetNameTitle(strNameRespPt,"");
3038 hnResponseZ[i]->SetNameTitle(strNameRespZ,"");
3039 hnResponseXi[i]->SetNameTitle(strNameRespXi,"");
3041 h2ResponsePt[i] = (TH2F*) hnResponsePt[i]->Projection(1,0);// note convention: yDim,xDim
3042 h2ResponseZ[i] = (TH2F*) hnResponseZ[i]->Projection(1,0); // note convention: yDim,xDim
3043 h2ResponseXi[i] = (TH2F*) hnResponseXi[i]->Projection(1,0);// note convention: yDim,xDim
3045 h2ResponsePt[i]->SetNameTitle(strNameTH2RespPt,"");
3046 h2ResponseZ[i]->SetNameTitle(strNameTH2RespZ,"");
3047 h2ResponseXi[i]->SetNameTitle(strNameTH2RespXi,"");
3052 Int_t axisGenPtTHn2 = -1;
3053 Int_t axisRecPtTHn2 = -1;
3055 for(Int_t d=0; d<hnResponsePt[i]->GetNdimensions(); d++){
3057 TString title = hnResponsePt[i]->GetAxis(d)->GetTitle();
3059 if(title.Contains("gen p_{T}")) axisGenPtTHn2 = d;
3060 if(title.Contains("rec p_{T}")) axisRecPtTHn2 = d;
3064 if(axisGenPtTHn2 == -1){
3065 Printf("%s:%d -- error axisGenPtTHn2",(char*)__FILE__,__LINE__);
3069 if(axisRecPtTHn2 == -1){
3070 Printf("%s:%d -- error axisRecPtTHn2",(char*)__FILE__,__LINE__);
3075 h1FFPtRec[i] = (TH1F*) hnResponsePt[i]->Projection(axisRecPtTHn2);// note convention: yDim,xDim
3076 h1FFZRec[i] = (TH1F*) hnResponseZ[i]->Projection(axisRecPtTHn2);// note convention: yDim,xDim
3077 h1FFXiRec[i] = (TH1F*) hnResponseXi[i]->Projection(axisRecPtTHn2);// note convention: yDim,xDim
3079 h1FFPtRec[i]->SetNameTitle(strNameRecPt,"");
3080 h1FFZRec[i]->SetNameTitle(strNameRecZ,"");
3081 h1FFXiRec[i]->SetNameTitle(strNameRecXi,"");
3083 h1FFPtGen[i] = (TH1F*) hnResponsePt[i]->Projection(axisGenPtTHn2);// note convention: yDim,xDim
3084 h1FFZGen[i] = (TH1F*) hnResponseZ[i]->Projection(axisGenPtTHn2);// note convention: yDim,xDim
3085 h1FFXiGen[i] = (TH1F*) hnResponseXi[i]->Projection(axisGenPtTHn2);// note convention: yDim,xDim
3087 h1FFPtGen[i]->SetNameTitle(strNameGenPt,"");
3088 h1FFZGen[i]->SetNameTitle(strNameGenZ,"");
3089 h1FFXiGen[i]->SetNameTitle(strNameGenXi,"");
3091 // normalize 1D projections
3093 if(fNHistoBinsPt[i]) h1FFPtRec[i] = (TH1F*) h1FFPtRec[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray());
3094 if(fNHistoBinsZ[i]) h1FFZRec[i] = (TH1F*) h1FFZRec[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3095 if(fNHistoBinsXi[i]) h1FFXiRec[i] = (TH1F*) h1FFXiRec[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3097 if(fNHistoBinsPt[i]) h1FFPtGen[i] = (TH1F*) h1FFPtGen[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray());
3098 if(fNHistoBinsZ[i]) h1FFZGen[i] = (TH1F*) h1FFZGen[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3099 if(fNHistoBinsXi[i]) h1FFXiGen[i] = (TH1F*) h1FFXiGen[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3101 NormalizeTH1(h1FFPtRec[i],fNJets->At(i));
3102 NormalizeTH1(h1FFZRec[i],fNJets->At(i));
3103 NormalizeTH1(h1FFXiRec[i],fNJets->At(i));
3105 NormalizeTH1(h1FFPtGen[i],fNJets->At(i));
3106 NormalizeTH1(h1FFZGen[i],fNJets->At(i));
3107 NormalizeTH1(h1FFXiGen[i],fNJets->At(i));
3109 // ratios 1D projections
3111 h1RatioPt[i] = (TH1F*) h1FFPtRec[i]->Clone(strNameRatioPt);
3112 h1RatioPt[i]->Reset();
3113 h1RatioPt[i]->Divide(h1FFPtRec[i],h1FFPtGen[i],1,1,"B");
3115 h1RatioZ[i] = (TH1F*) h1FFZRec[i]->Clone(strNameRatioZ);
3116 h1RatioZ[i]->Reset();
3117 h1RatioZ[i]->Divide(h1FFZRec[i],h1FFZGen[i],1,1,"B");
3119 h1RatioXi[i] = (TH1F*) h1FFXiRec[i]->Clone(strNameRatioXi);
3120 h1RatioXi[i]->Reset();
3121 h1RatioXi[i]->Divide(h1FFXiRec[i],h1FFXiGen[i],1,1,"B");
3127 TString outfileOption = "RECREATE";
3128 if(updateOutfile) outfileOption = "UPDATE";
3130 TFile out(strOutfile,outfileOption);
3133 Printf("%s:%d -- error opening response matrix output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
3137 if(fDebug>0) Printf("%s:%d -- write response matrices to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
3139 //if(strdir && strdir.Length()){
3140 // TDirectory* dir = out.mkdir(strdir);
3144 for(Int_t i=0; i<fNJetPtSlices; i++){
3146 hnResponsePt[i]->Write();
3147 hnResponseXi[i]->Write();
3148 hnResponseZ[i]->Write();
3150 h2ResponsePt[i]->Write();
3151 h2ResponseXi[i]->Write();
3152 h2ResponseZ[i]->Write();
3154 h1FFPtRec[i]->Write();
3155 h1FFZRec[i]->Write();
3156 h1FFXiRec[i]->Write();
3158 h1FFPtGen[i]->Write();
3159 h1FFZGen[i]->Write();
3160 h1FFXiGen[i]->Write();
3162 h1RatioPt[i]->Write();
3163 h1RatioZ[i]->Write();
3164 h1RatioXi[i]->Write();
3171 //______________________________________________________________________________________________________
3172 void AliFragmentationFunctionCorrections::ReadResponse(TString strfile, TString strdir, TString strlist)
3174 // read response matrices from file
3175 // argument strlist optional - read from directory strdir if not specified
3176 // note: THnSparse are not rebinned
3178 THnSparse* hRespPt[fNJetPtSlices];
3179 THnSparse* hRespZ[fNJetPtSlices];
3180 THnSparse* hRespXi[fNJetPtSlices];
3182 for(Int_t i=0; i<fNJetPtSlices; i++) hRespPt[i] = 0;
3183 for(Int_t i=0; i<fNJetPtSlices; i++) hRespZ[i] = 0;
3184 for(Int_t i=0; i<fNJetPtSlices; i++) hRespXi[i] = 0;
3186 TFile f(strfile,"READ");
3189 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
3193 if(fDebug>0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3195 if(strdir && strdir.Length()) gDirectory->cd(strdir);
3199 if(strlist && strlist.Length()){
3201 if(!(list = (TList*) gDirectory->Get(strlist))){
3202 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
3207 for(Int_t i=0; i<fNJetPtSlices; i++){
3209 Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
3210 Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
3212 TString strNameRespPt(Form("hnResponsePt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3213 TString strNameRespZ(Form("hnResponseZ_%02d_%02d",jetPtLoLim,jetPtUpLim));
3214 TString strNameRespXi(Form("hnResponseXi_%02d_%02d",jetPtLoLim,jetPtUpLim));
3217 hRespPt[i] = (THnSparse*) list->FindObject(strNameRespPt);
3218 hRespZ[i] = (THnSparse*) list->FindObject(strNameRespZ);
3219 hRespXi[i] = (THnSparse*) list->FindObject(strNameRespXi);
3222 hRespPt[i] = (THnSparse*) gDirectory->Get(strNameRespPt);
3223 hRespZ[i] = (THnSparse*) gDirectory->Get(strNameRespZ);
3224 hRespXi[i] = (THnSparse*) gDirectory->Get(strNameRespXi);
3228 Printf("%s:%d -- error retrieving response %s", (char*)__FILE__,__LINE__,strNameRespPt.Data());
3232 Printf("%s:%d -- error retrieving response %s", (char*)__FILE__,__LINE__,strNameRespZ.Data());
3236 Printf("%s:%d -- error retrieving response %s", (char*)__FILE__,__LINE__,strNameRespXi.Data());
3239 // if(0){ // can't rebin THnSparse ...
3240 // if(fNHistoBinsPt[i]) hRespPt[i]->SetBinEdges(0,fHistoBinsPt[i]->GetArray());
3241 // if(fNHistoBinsPt[i]) hRespPt[i]->SetBinEdges(1,fHistoBinsPt[i]->GetArray());
3243 // if(fNHistoBinsZ[i]) hRespZ[i]->SetBinEdges(0,fHistoBinsZ[i]->GetArray());
3244 // if(fNHistoBinsZ[i]) hRespZ[i]->SetBinEdges(1,fHistoBinsZ[i]->GetArray());
3246 // if(fNHistoBinsXi[i]) hRespXi[i]->SetBinEdges(0,fHistoBinsXi[i]->GetArray());
3247 // if(fNHistoBinsXi[i]) hRespXi[i]->SetBinEdges(1,fHistoBinsXi[i]->GetArray());
3251 } // jet slices loop
3253 f.Close(); // THnSparse pointers still valid even if file closed
3255 // for(Int_t i=0; i<fNJetPtSlices; i++){ // no copy c'tor ...
3256 // if(hRespPt[i]) new(fhnResponsePt[i]) THnSparseF(*hRespPt[i]);
3257 // if(hRespZ[i]) new(fhnResponseZ[i]) THnSparseF(*hRespZ[i]);
3258 // if(hRespXi[i]) new(fhnResponseXi[i]) THnSparseF(*hRespXi[i]);
3261 for(Int_t i=0; i<fNJetPtSlices; i++){
3262 fhnResponsePt[i] = hRespPt[i];
3263 fhnResponseZ[i] = hRespZ[i];
3264 fhnResponseXi[i] = hRespXi[i];
3268 //______________________________________________________________________________________________________________________
3269 void AliFragmentationFunctionCorrections::ReadPriors(TString strfile,const Int_t type)
3271 // read priors from file: rec primaries, gen pt dist
3273 if(fDebug>0) Printf("%s:%d -- read priors from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3275 // temporary histos to store pointers from file
3276 TH1F* hist[fNJetPtSlices];
3278 for(Int_t i=0; i<fNJetPtSlices; i++) hist[i] = 0;
3280 TFile f(strfile,"READ");
3283 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
3287 for(Int_t i=0; i<fNJetPtSlices; i++){
3289 Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
3290 Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
3294 if(type == kFlagPt) strName.Form("h1FFTrackPtRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim);
3295 if(type == kFlagZ) strName.Form("h1FFZRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim);
3296 if(type == kFlagXi) strName.Form("h1FFXiRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim);
3298 hist[i] = (TH1F*) gDirectory->Get(strName);
3301 Printf("%s:%d -- error retrieving prior %s", (char*)__FILE__,__LINE__,strName.Data());
3305 //if(fNHistoBinsPt[i]) hist[i] = (TH1F*) hist[i]->Rebin(fNHistoBinsPt[i],hist[i]->GetName()+"_rebin",fHistoBinsPt[i]->GetArray());
3307 if(hist[i]) hist[i]->SetDirectory(0);
3309 } // jet slices loop
3314 for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
3315 if(hist[i] && type == kFlagPt) new(fh1FFTrackPtPrior[i]) TH1F(*hist[i]);
3316 if(hist[i] && type == kFlagZ) new(fh1FFZPrior[i]) TH1F(*hist[i]);
3317 if(hist[i] && type == kFlagXi) new(fh1FFXiPrior[i]) TH1F(*hist[i]);
3321 //_____________________________________________________
3322 // void AliFragmentationFunctionCorrections::RatioRecGen()
3324 // // create ratio reconstructed over generated FF
3325 // // use current highest corrLevel
3327 // Printf("%s:%d -- build ratio rec.gen, corrLevel %d",(char*)__FILE__,__LINE__,fNCorrectionLevels-1);
3329 // for(Int_t i=0; i<fNJetPtSlices; i++){
3331 // TH1F* histPtRec = fCorrFF[fNCorrectionLevels-1]->GetTrackPt(i); // levels -1: latest corr level
3332 // TH1F* histZRec = fCorrFF[fNCorrectionLevels-1]->GetZ(i); // levels -1: latest corr level
3333 // TH1F* histXiRec = fCorrFF[fNCorrectionLevels-1]->GetXi(i); // levels -1: latest corr level
3335 // TH1F* histPtGen = fh1FFTrackPtGenPrim[i];
3336 // TH1F* histZGen = fh1FFZGenPrim[i];
3337 // TH1F* histXiGen = fh1FFXiGenPrim[i];
3339 // TString histNamePt = histPtRec->GetName();
3340 // TString histNameZ = histZRec->GetName();
3341 // TString histNameXi = histXiRec->GetName();
3343 // histNamePt.ReplaceAll("fh1FF","fh1FFRatioRecGen");
3344 // histNameZ.ReplaceAll("fh1FF","fh1FFRatioRecGen");
3345 // histNameXi.ReplaceAll("fh1FF","fh1FFRatioRecGen");
3348 // TH1F* hRatioRecGenPt = (TH1F*) histPtRec->Clone(histNamePt);
3349 // hRatioRecGenPt->Reset();
3350 // hRatioRecGenPt->Divide(histPtRec,histPtGen,1,1,"B");
3352 // TH1F* hRatioRecGenZ = (TH1F*) histZRec->Clone(histNameZ);
3353 // hRatioRecGenZ->Reset();
3354 // hRatioRecGenZ->Divide(histZRec,histZGen,1,1,"B");
3356 // TH1F* hRatioRecGenXi = (TH1F*) histXiRec->Clone(histNameXi);
3357 // hRatioRecGenXi->Reset();
3358 // hRatioRecGenXi->Divide(histXiRec,histXiGen,1,1,"B");
3360 // new(fh1FFRatioRecGenPt[i]) TH1F(*hRatioRecGenPt);
3361 // new(fh1FFRatioRecGenZ[i]) TH1F(*hRatioRecGenZ);
3362 // new(fh1FFRatioRecGenXi[i]) TH1F(*hRatioRecGenXi);
3366 // //___________________________________________________________
3367 // void AliFragmentationFunctionCorrections::RatioRecPrimaries()
3369 // // create ratio reconstructed tracks over reconstructed primaries
3370 // // use raw FF (corrLevel 0)
3372 // Printf("%s:%d -- build ratio rec tracks /rec primaries",(char*)__FILE__,__LINE__);
3374 // for(Int_t i=0; i<fNJetPtSlices; i++){
3376 // const Int_t corrLevel = 0;
3378 // TH1F* histPtRec = fCorrFF[corrLevel]->GetTrackPt(i); // levels -1: latest corr level
3379 // TH1F* histZRec = fCorrFF[corrLevel]->GetZ(i); // levels -1: latest corr level
3380 // TH1F* histXiRec = fCorrFF[corrLevel]->GetXi(i); // levels -1: latest corr level
3382 // TH1F* histPtRecPrim = fh1FFTrackPtRecPrim[i];
3383 // TH1F* histZRecPrim = fh1FFZRecPrim[i];
3384 // TH1F* histXiRecPrim = fh1FFXiRecPrim[i];
3386 // TString histNamePt = histPtRec->GetName();
3387 // TString histNameZ = histZRec->GetName();
3388 // TString histNameXi = histXiRec->GetName();
3390 // histNamePt.ReplaceAll("fh1FF","fh1FFRatioRecPrim");
3391 // histNameZ.ReplaceAll("fh1FF","fh1FFRatioRecPrim");
3392 // histNameXi.ReplaceAll("fh1FF","fh1FFRatioRecPrim");
3395 // TH1F* hRatioRecPrimPt = (TH1F*) histPtRec->Clone(histNamePt);
3396 // hRatioRecPrimPt->Reset();
3397 // hRatioRecPrimPt->Divide(histPtRec,histPtRecPrim,1,1,"B");
3399 // TH1F* hRatioRecPrimZ = (TH1F*) histZRec->Clone(histNameZ);
3400 // hRatioRecPrimZ->Reset();
3401 // hRatioRecPrimZ->Divide(histZRec,histZRecPrim,1,1,"B");
3403 // TH1F* hRatioRecPrimXi = (TH1F*) histXiRec->Clone(histNameXi);
3404 // hRatioRecPrimXi->Reset();
3405 // hRatioRecPrimXi->Divide(histXiRec,histXiRecPrim,1,1,"B");
3408 // new(fh1FFRatioRecPrimPt[i]) TH1F(*hRatioRecPrimPt);
3409 // new(fh1FFRatioRecPrimZ[i]) TH1F(*hRatioRecPrimZ);
3410 // new(fh1FFRatioRecPrimXi[i]) TH1F(*hRatioRecPrimXi);
3414 // __________________________________________________________________________________
3415 void AliFragmentationFunctionCorrections::ProjectJetResponseMatrices(TString strOutfile)
3418 // project response matrices on both axes:
3419 // FF for rec primaries, in terms of generated and reconstructed momentum
3420 // write FF and ratios to outFile
3422 Printf("%s:%d -- project response matrices, write to %s",(char*)__FILE__,__LINE__,strOutfile.Data());
3424 TH1F* hFFPtRec[fNJetPtSlices];
3425 TH1F* hFFZRec[fNJetPtSlices];
3426 TH1F* hFFXiRec[fNJetPtSlices];
3428 TH1F* hFFPtGen[fNJetPtSlices];
3429 TH1F* hFFZGen[fNJetPtSlices];
3430 TH1F* hFFXiGen[fNJetPtSlices];
3432 TH1F* hRatioPt[fNJetPtSlices];
3433 TH1F* hRatioZ[fNJetPtSlices];
3434 TH1F* hRatioXi[fNJetPtSlices];
3437 Int_t axisGenPt = 1;
3438 Int_t axisRecPt = 0;
3440 for(Int_t i=0; i<fNJetPtSlices; i++){
3442 Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
3443 Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
3445 TString strNameRecPt(Form("h1FFTrackPtRecPrim_recPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3446 TString strNameRecZ(Form("h1FFZRecPrim_recPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3447 TString strNameRecXi(Form("h1FFXiRecPrim_recPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3449 TString strNameGenPt(Form("h1FFTrackPtRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3450 TString strNameGenZ(Form("h1FFZRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3451 TString strNameGenXi(Form("h1FFXiRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3453 TString strNameRatioPt(Form("h1RatioTrackPtRecPrim_%02d_%02d",jetPtLoLim,jetPtUpLim));
3454 TString strNameRatioZ(Form("h1RatioZRecPrim_%02d_%02d",jetPtLoLim,jetPtUpLim));
3455 TString strNameRatioXi(Form("h1RatioXiRecPrim_%02d_%02d",jetPtLoLim,jetPtUpLim));
3458 hFFPtRec[i] = (TH1F*) fhnResponsePt[i]->Projection(axisRecPt);// note convention: yDim,xDim
3459 hFFZRec[i] = (TH1F*) fhnResponseZ[i]->Projection(axisRecPt);// note convention: yDim,xDim
3460 hFFXiRec[i] = (TH1F*) fhnResponseXi[i]->Projection(axisRecPt);// note convention: yDim,xDim
3462 hFFPtRec[i]->SetNameTitle(strNameRecPt,"");
3463 hFFZRec[i]->SetNameTitle(strNameRecZ,"");
3464 hFFXiRec[i]->SetNameTitle(strNameRecXi,"");
3467 hFFPtGen[i] = (TH1F*) fhnResponsePt[i]->Projection(axisGenPt);// note convention: yDim,xDim
3468 hFFZGen[i] = (TH1F*) fhnResponseZ[i]->Projection(axisGenPt);// note convention: yDim,xDim
3469 hFFXiGen[i] = (TH1F*) fhnResponseXi[i]->Projection(axisGenPt);// note convention: yDim,xDim
3471 hFFPtGen[i]->SetNameTitle(strNameGenPt,"");
3472 hFFZGen[i]->SetNameTitle(strNameGenZ,"");
3473 hFFXiGen[i]->SetNameTitle(strNameGenXi,"");
3476 if(fNHistoBinsPt[i]) hFFPtRec[i] = (TH1F*) hFFPtRec[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray());
3477 if(fNHistoBinsZ[i]) hFFZRec[i] = (TH1F*) hFFZRec[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3478 if(fNHistoBinsXi[i]) hFFXiRec[i] = (TH1F*) hFFXiRec[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3480 if(fNHistoBinsPt[i]) hFFPtGen[i] = (TH1F*) hFFPtGen[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray());
3481 if(fNHistoBinsZ[i]) hFFZGen[i] = (TH1F*) hFFZGen[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3482 if(fNHistoBinsXi[i]) hFFXiGen[i] = (TH1F*) hFFXiGen[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3484 NormalizeTH1(hFFPtGen[i],fNJets->At(i));
3485 NormalizeTH1(hFFZGen[i],fNJets->At(i));
3486 NormalizeTH1(hFFXiGen[i],fNJets->At(i));
3488 NormalizeTH1(hFFPtRec[i],fNJets->At(i));
3489 NormalizeTH1(hFFZRec[i],fNJets->At(i));
3490 NormalizeTH1(hFFXiRec[i],fNJets->At(i));
3493 hRatioPt[i] = (TH1F*) hFFPtRec[i]->Clone(strNameRatioPt);
3494 hRatioPt[i]->Reset();
3495 hRatioPt[i]->Divide(hFFPtRec[i],hFFPtGen[i],1,1,"B");
3497 hRatioZ[i] = (TH1F*) hFFZRec[i]->Clone(strNameRatioZ);
3498 hRatioZ[i]->Reset();
3499 hRatioZ[i]->Divide(hFFZRec[i],hFFZGen[i],1,1,"B");
3501 hRatioXi[i] = (TH1F*) hFFXiRec[i]->Clone(strNameRatioXi);
3502 hRatioXi[i]->Reset();
3503 hRatioXi[i]->Divide(hFFXiRec[i],hFFXiGen[i],1,1,"B");
3510 TFile out(strOutfile,"RECREATE");
3513 Printf("%s:%d -- error opening response matrix output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
3517 for(Int_t i=0; i<fNJetPtSlices; i++){
3519 hFFPtRec[i]->Write();
3520 hFFZRec[i]->Write();
3521 hFFXiRec[i]->Write();
3523 hFFPtGen[i]->Write();
3524 hFFZGen[i]->Write();
3525 hFFXiGen[i]->Write();
3527 hRatioPt[i]->Write();
3528 hRatioZ[i]->Write();
3529 hRatioXi[i]->Write();
3535 // ____________________________________________________________________________________________________________________________
3536 void AliFragmentationFunctionCorrections::ProjectSingleResponseMatrix(TString strOutfile, Bool_t updateOutfile, TString strOutDir )
3538 // project response matrix on both axes:
3539 // pt spec for rec primaries, in terms of generated and reconstructed momentum
3540 // write spec and ratios to outFile
3542 Printf("%s:%d -- project single pt response matrix, write to %s",(char*)__FILE__,__LINE__,strOutfile.Data());
3548 Int_t axisGenPt = 1;
3549 Int_t axisRecPt = 0;
3551 TString strNameRecPt = "h1SpecTrackPtRecPrim_recPt";
3552 TString strNameGenPt = "h1SpecTrackPtRecPrim_genPt";
3553 TString strNameRatioPt = "h1RatioTrackPtRecPrim";
3555 hSpecPtRec = (TH1F*) fhnResponseSinglePt->Projection(axisRecPt);// note convention: yDim,xDim
3556 hSpecPtRec->SetNameTitle(strNameRecPt,"");
3558 hSpecPtGen = (TH1F*) fhnResponseSinglePt->Projection(axisGenPt);// note convention: yDim,xDim
3559 hSpecPtGen->SetNameTitle(strNameGenPt,"");
3561 hRatioPt = (TH1F*) hSpecPtRec->Clone(strNameRatioPt);
3563 hRatioPt->Divide(hSpecPtRec,hSpecPtGen,1,1,"B");
3565 TString outfileOption = "RECREATE";
3566 if(updateOutfile) outfileOption = "UPDATE";
3568 TFile out(strOutfile,outfileOption);
3571 Printf("%s:%d -- error opening reponse matrix projections output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
3576 if(strOutDir && strOutDir.Length()){
3579 if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd();
3581 dir = out.mkdir(strOutDir);
3586 hSpecPtRec->Write();
3587 hSpecPtGen->Write();
3594 //__________________________________________________________________________________________________________________________________________________________________
3595 void AliFragmentationFunctionCorrections::RebinHisto(const Int_t jetPtSlice, const Int_t nBinsLimits, Double_t* binsLimits, Double_t* binsWidth, const Int_t type)
3597 // rebin histo, rescale bins according to new width
3598 // only correct for input histos with equal bin size
3600 // args: jetPtSlice, type, use current corr level
3602 // function takes array of limits and widths (e.g. 1 GeV bins up to 10 GeV, 2 GeV width up to 20 GeV, ...)
3603 // array size of binsLimits: nBinsLimits
3604 // array size of binsWidth: nBinsLimits-1
3605 // binsLimits have to be in increasing order
3606 // if binning undefined for any slice, original binning will be kept
3609 Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
3613 if(jetPtSlice>=fNJetPtSlices){
3614 Printf("%s:%d -- jetPtSlice %d exceeds max",(char*)__FILE__,__LINE__,jetPtSlice);
3619 Double_t binLimitMin = binsLimits[0];
3620 Double_t binLimitMax = binsLimits[nBinsLimits-1];
3622 Double_t binLimit = binLimitMin; // start value
3624 Int_t sizeUpperLim = 1000; //static_cast<Int_t>(binLimitMax/binsWidth[0])+1; - only if first bin has smallest width, but not case for dN/dxi ...
3625 TArrayD binsArray(sizeUpperLim);
3627 binsArray.SetAt(binLimitMin,nBins++);
3629 while(binLimit<binLimitMax && nBins<sizeUpperLim){
3631 Int_t currentSlice = -1;
3632 for(Int_t i=0; i<nBinsLimits; i++){
3633 if(binLimit >= 0.999*binsLimits[i]) currentSlice = i; // 0.999 numerical saftey factor
3636 Double_t currentBinWidth = binsWidth[currentSlice];
3637 binLimit += currentBinWidth;
3639 binsArray.SetAt(binLimit,nBins++);
3644 if(type == kFlagPt) hist = fCorrFF[fNCorrectionLevels-1]->GetTrackPt(jetPtSlice);
3645 else if(type == kFlagZ) hist = fCorrFF[fNCorrectionLevels-1]->GetZ(jetPtSlice);
3646 else if(type == kFlagXi) hist = fCorrFF[fNCorrectionLevels-1]->GetXi(jetPtSlice);
3647 else if(type == kFlagSinglePt) hist = fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->GetTrackPt(0);
3649 Printf("%s%d unknown type",(char*)__FILE__,__LINE__);
3653 Double_t binWidthNoRebin = hist->GetBinWidth(1);
3655 Double_t* bins = binsArray.GetArray();
3657 hist = (TH1F*) hist->Rebin(nBins-1,"",bins);
3659 for(Int_t bin=0; bin <= hist->GetNbinsX(); bin++){
3661 Double_t binWidthRebin = hist->GetBinWidth(bin);
3662 Double_t scaleF = binWidthNoRebin / binWidthRebin;
3664 Double_t binCont = hist->GetBinContent(bin);
3665 Double_t binErr = hist->GetBinError(bin);
3670 hist->SetBinContent(bin,binCont);
3671 hist->SetBinError(bin,binErr);
3676 TH1F* temp = new TH1F(*hist);
3678 if(type == kFlagPt) fCorrFF[fNCorrectionLevels-1]->ReplaceCorrHistos(jetPtSlice,temp,0,0);
3679 if(type == kFlagZ) fCorrFF[fNCorrectionLevels-1]->ReplaceCorrHistos(jetPtSlice,0,temp,0);
3680 if(type == kFlagXi) fCorrFF[fNCorrectionLevels-1]->ReplaceCorrHistos(jetPtSlice,0,0,temp);
3681 if(type == kFlagSinglePt) fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->ReplaceCorrHistos(0,temp,0,0);
3686 //__________________________________________________________________________________________________________________________________________________________________
3687 void AliFragmentationFunctionCorrections::WriteJetSpecResponse(TString strInfile, TString strdir, TString strlist/*, TString strOutfile*/)
3690 if(fDebug>0) Printf("%s:%d -- read jet spectrum response matrix from file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
3692 if(strdir && strdir.Length()) gDirectory->cd(strdir);
3696 if(strlist && strlist.Length()){
3697 if(!(list = (TList*) gDirectory->Get(strlist))){
3698 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
3702 if(list == 0)return; // catch strlist.Lenght() == 0;
3704 THnSparse* hn6ResponseJetPt = (THnSparse*) list->FindObject("fhnCorrelation");
3706 Int_t axis6RecJetPt = 0;
3707 Int_t axis6GenJetPt = 3;
3709 hn6ResponseJetPt->GetAxis(axis6RecJetPt)->SetTitle("rec jet p_{T} (GeV/c)");
3710 hn6ResponseJetPt->GetAxis(axis6GenJetPt)->SetTitle("gen jet p_{T} (GeV/c)");
3712 Int_t nBinsRecPt = hn6ResponseJetPt->GetAxis(axis6RecJetPt)->GetNbins();
3713 Double_t loLimRecPt = hn6ResponseJetPt->GetAxis(axis6RecJetPt)->GetBinLowEdge(1);
3714 Double_t upLimRecPt = hn6ResponseJetPt->GetAxis(axis6RecJetPt)->GetBinUpEdge(nBinsRecPt);
3716 Int_t nBinsGenPt = hn6ResponseJetPt->GetAxis(axis6GenJetPt)->GetNbins();
3717 Double_t loLimGenPt = hn6ResponseJetPt->GetAxis(axis6GenJetPt)->GetBinLowEdge(1);
3718 Double_t upLimGenPt = hn6ResponseJetPt->GetAxis(axis6GenJetPt)->GetBinUpEdge(nBinsGenPt);
3720 Int_t nBinsTrackPt = 200;
3721 Int_t loLimTrackPt = 0;
3722 Int_t upLimTrackPt = 200;
3725 Int_t nBinsResponse[4] = {nBinsRecPt,nBinsTrackPt,nBinsGenPt,nBinsTrackPt};
3726 Double_t binMinResponse[4] = {loLimRecPt,loLimTrackPt,loLimGenPt,loLimTrackPt};
3727 Double_t binMaxResponse[4] = {upLimRecPt,upLimTrackPt,upLimGenPt,upLimTrackPt};
3729 const char* labelsResponseSinglePt[4] = {"rec jet p_{T} (GeV/c)", "rec track p_{T} (GeV/c)", "gen jet p_{T} (GeV/c)", "gen track p_{T} (GeV/c)"};
3731 THnSparseD* hn4ResponseTrackPtJetPt = new THnSparseD("hn4ResponseTrackPtJetPt","",4,nBinsResponse,binMinResponse,binMaxResponse);
3733 for(Int_t i=0; i<4; i++){
3734 hn4ResponseTrackPtJetPt->GetAxis(i)->SetTitle(labelsResponseSinglePt[i]);
3743 //_____________________________________________________________________________________________________________________________________
3744 void AliFragmentationFunctionCorrections::ReadSingleTrackEfficiency(TString strfile, TString strdir, TString strlist, TString strname)
3747 ReadSingleTrackCorrection(strfile,strdir,strlist,strname,kFlagEfficiency);
3751 //_____________________________________________________________________________________________________________________________________
3752 void AliFragmentationFunctionCorrections::ReadSingleTrackResponse(TString strfile, TString strdir, TString strlist, TString strname)
3755 ReadSingleTrackCorrection(strfile,strdir,strlist,strname,kFlagResponse);
3759 //_____________________________________________________________________________________________________________________________________
3760 void AliFragmentationFunctionCorrections::ReadSingleTrackSecCorr(TString strfile, TString strdir, TString strlist, TString strname)
3763 ReadSingleTrackCorrection(strfile,strdir,strlist,strname,kFlagSecondaries);
3767 //______________________________________________________________________________________________________________________________________________________
3768 void AliFragmentationFunctionCorrections::ReadSingleTrackCorrection(TString strfile, TString strdir, TString strlist, TString strname, const Int_t type)
3770 // read single track correction (pt) from file
3771 // type: efficiency / response / secondaries correction
3773 if(!((type == kFlagEfficiency) || (type == kFlagResponse) || (type == kFlagSecondaries))){
3774 Printf("%s:%d -- no such correction ",(char*)__FILE__,__LINE__);
3778 TFile f(strfile,"READ");
3781 Printf("%s:%d -- error opening file %s", (char*)__FILE__,__LINE__,strfile.Data());
3785 if(fDebug>0 && type==kFlagEfficiency) Printf("%s:%d -- read single track corr from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3786 if(fDebug>0 && type==kFlagResponse) Printf("%s:%d -- read single track response from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3787 if(fDebug>0 && type==kFlagSecondaries) Printf("%s:%d -- read single track secondaries corr from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3789 if(strdir && strdir.Length()) gDirectory->cd(strdir);
3793 if(strlist && strlist.Length()){
3795 if(!(list = (TList*) gDirectory->Get(strlist))){
3796 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
3801 TH1F* h1CorrHist = 0; // common TObject pointer not possible, need SetDirectory() later
3802 THnSparse* hnCorrHist = 0;
3804 if(type == kFlagEfficiency || type == kFlagSecondaries){
3806 if(list) h1CorrHist = (TH1F*) list->FindObject(strname);
3807 else h1CorrHist = (TH1F*) gDirectory->Get(strname);
3810 Printf("%s:%d -- error retrieving histo %s", (char*)__FILE__,__LINE__,strname.Data());
3815 else if(type == kFlagResponse){
3817 if(list) hnCorrHist = (THnSparse*) list->FindObject(strname);
3818 else hnCorrHist = (THnSparse*) gDirectory->Get(strname);
3821 Printf("%s:%d -- error retrieving histo %s", (char*)__FILE__,__LINE__,strname.Data());
3827 if(h1CorrHist) h1CorrHist->SetDirectory(0);
3828 //if(hnCorrHist) hnCorrHist->SetDirectory(0);
3832 if(type == kFlagEfficiency) fh1EffSinglePt = h1CorrHist;
3833 else if(type == kFlagResponse) fhnResponseSinglePt = hnCorrHist;
3834 else if(type == kFlagSecondaries) fh1SecCorrSinglePt = h1CorrHist;
3838 //________________________________________________________________________________________________________________
3839 void AliFragmentationFunctionCorrections::ReadRawPtSpec(TString strInfile, TString strID)
3841 // read track pt spec from task ouput - standard dir/list
3843 TString strdir = "PWGJE_FragmentationFunction_" + strID;
3844 TString strlist = "fracfunc_" + strID;
3846 ReadRawPtSpec(strInfile,strdir,strlist);
3849 //_______________________________________________________________________________________________________
3850 void AliFragmentationFunctionCorrections::ReadRawPtSpec(TString strfile, TString strdir, TString strlist)
3852 // get raw pt spectra from file
3856 fNCorrectionLevelsSinglePt = 0;
3857 fCorrSinglePt = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels];
3858 AddCorrectionLevelSinglePt(); // first 'correction' level = raw spectrum
3860 // get raw pt spec from input file, normalize
3862 TFile f(strfile,"READ");
3865 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
3869 if(fDebug>0) Printf("%s:%d -- read raw spectra from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3871 gDirectory->cd(strdir);
3875 if(!(list = (TList*) gDirectory->Get(strlist))){
3876 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
3880 TString hnameTrackPt("fh1TrackQAPtRecCuts");
3881 TString hnameEvtSel("fh1EvtSelection");
3883 TH1F* fh1TrackPt = (TH1F*) list->FindObject(hnameTrackPt);
3884 TH1F* fh1EvtSel = (TH1F*) list->FindObject(hnameEvtSel);
3886 if(!fh1TrackPt){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameTrackPt.Data()); return; }
3887 if(!fh1EvtSel) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameEvtSel.Data()); return; }
3890 //Float_t nEvents = fh1EvtSel->GetBinContent(fh1EvtSel->FindBin(0));
3892 // evts after physics selection
3893 Float_t nEvents = fh1EvtSel->GetEntries()
3894 - fh1EvtSel->GetBinContent(fh1EvtSel->FindBin(1)) // evts rejected by trigger selection ( = PhysicsSelection
3895 - fh1EvtSel->GetBinContent(fh1EvtSel->FindBin(2)); // evts with wrong centrality class
3898 fh1TrackPt->SetDirectory(0);
3903 NormalizeTH1(fh1TrackPt,nEvents);
3905 // raw FF = corr level 0
3906 fCorrSinglePt[0]->AddCorrHistos(0,fh1TrackPt);
3910 //_______________________________________________________________________________________________________
3911 void AliFragmentationFunctionCorrections::ReadRawPtSpecQATask(TString strfile, TString strdir, TString strlist)
3913 // get raw pt spectra from file
3914 // for output from Martas QA task
3918 fNCorrectionLevelsSinglePt = 0;
3919 fCorrSinglePt = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels];
3920 AddCorrectionLevelSinglePt(); // first 'correction' level = raw spectrum
3922 // get raw pt spec from input file, normalize
3924 TFile f(strfile,"READ");
3927 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
3931 if(fDebug>0) Printf("%s:%d -- read raw pt spec from QA task output file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3933 gDirectory->cd(strdir);
3937 if(!(list = (TList*) gDirectory->Get(strlist))){
3938 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
3942 TString hnameTrackPt("fPtSel");
3943 TString hnameEvtSel("fNEventAll");
3945 TH1F* fh1TrackPt = (TH1F*) list->FindObject(hnameTrackPt);
3946 TH1F* fh1EvtSel = (TH1F*) list->FindObject(hnameEvtSel);
3948 if(!fh1TrackPt){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameTrackPt.Data()); return; }
3949 if(!fh1EvtSel) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameEvtSel.Data()); return; }
3952 // evts after physics selection
3953 Float_t nEvents = fh1EvtSel->GetEntries();
3955 fh1TrackPt->SetDirectory(0);
3960 NormalizeTH1(fh1TrackPt,nEvents);
3962 // raw FF = corr level 0
3963 fCorrSinglePt[0]->AddCorrHistos(0,fh1TrackPt);
3966 // ________________________________________________________
3967 void AliFragmentationFunctionCorrections::EffCorrSinglePt()
3969 // apply efficiency correction to inclusive track pt spec
3971 AddCorrectionLevelSinglePt("eff");
3973 TH1F* histPt = fCorrSinglePt[fNCorrectionLevelsSinglePt-2]->GetTrackPt(0);
3975 if(histPt->GetNbinsX() != fh1EffSinglePt->GetNbinsX()) Printf("%s:%d: inconsistency pt spec and eff corr bins - rebin effCorr ...", (char*)__FILE__,__LINE__);
3977 TString histNamePt = histPt->GetName();
3978 TH1F* hTrackPtEffCorr = (TH1F*) histPt->Clone(histNamePt);
3979 hTrackPtEffCorr->Divide(histPt,fh1EffSinglePt,1,1,"");
3981 fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->AddCorrHistos(0,hTrackPtEffCorr);
3984 //___________________________________________________________________________________________________________________________
3985 void AliFragmentationFunctionCorrections::UnfoldSinglePt(const Int_t nIter, const Bool_t useCorrelatedErrors)
3987 // unfolde inclusive dN/dpt spectra
3989 AddCorrectionLevelSinglePt("unfold");
3991 TH1F* hist = fCorrSinglePt[fNCorrectionLevelsSinglePt-2]->GetTrackPt(0); // level -2: before unfolding, level -1: unfolded
3992 THnSparse* hnResponse = fhnResponseSinglePt;
3994 TString histNameTHn = hist->GetName();
3995 if(histNameTHn.Contains("TH1")) histNameTHn.ReplaceAll("TH1","THn");
3996 if(histNameTHn.Contains("fPt")) histNameTHn.ReplaceAll("fPt","fhnPt");
3999 TString histNameBackFolded = hist->GetName();
4000 histNameBackFolded.Append("_backfold");
4002 TString histNameRatioFolded = hist->GetName();
4003 if(histNameRatioFolded.Contains("fh1")) histNameRatioFolded.ReplaceAll("fh1","hRatio");
4004 if(histNameRatioFolded.Contains("fPt")) histNameRatioFolded.ReplaceAll("fPt","hRatioPt");
4005 histNameRatioFolded.Append("_unfold");
4007 TString histNameRatioBackFolded = hist->GetName();
4008 if(histNameRatioBackFolded.Contains("fh1")) histNameRatioBackFolded.ReplaceAll("fh1","hRatio");
4009 if(histNameRatioBackFolded.Contains("fPt")) histNameRatioBackFolded.ReplaceAll("fPt","hRatioPt");
4010 histNameRatioBackFolded.Append("_backfold");
4012 THnSparse* hnHist = TH1toSparse(hist,histNameTHn,hist->GetTitle());
4013 THnSparse* hnFlatEfficiency = TH1toSparse(hist,"fhnEfficiency","eff",kTRUE); // could optionally also use real eff
4015 THnSparse* hnUnfolded
4016 = Unfold(hnHist,hnResponse,hnFlatEfficiency,nIter,useCorrelatedErrors);
4018 TH1F* hUnfolded = (TH1F*) hnUnfolded->Projection(0);
4019 hUnfolded->SetNameTitle(hist->GetName(),hist->GetTitle());
4021 fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->AddCorrHistos(0,hUnfolded);
4023 // backfolding: apply response matrix to unfolded spectrum
4024 TH1F* hBackFolded = ApplyResponse(hUnfolded,hnResponse);
4025 hBackFolded->SetNameTitle(histNameBackFolded,hist->GetTitle());
4027 fh1SingleTrackPtBackFolded = hBackFolded;
4030 // ratio unfolded to original histo
4031 TH1F* hRatioUnfolded = (TH1F*) hUnfolded->Clone(histNameRatioFolded);
4032 hRatioUnfolded->Reset();
4033 hRatioUnfolded->Divide(hUnfolded,hist,1,1,"B");
4035 fh1RatioSingleTrackPtFolded = hRatioUnfolded;
4038 // ratio backfolded to original histo
4039 TH1F* hRatioBackFolded = (TH1F*) hBackFolded->Clone(histNameRatioBackFolded);
4040 hRatioBackFolded->Reset();
4041 hRatioBackFolded->Divide(hBackFolded,hist,1,1,"B");
4043 fh1RatioSingleTrackPtBackFolded = hRatioBackFolded;
4046 delete hnFlatEfficiency;
4050 // ________________________________________________________
4051 void AliFragmentationFunctionCorrections::SecCorrSinglePt()
4053 // apply efficiency correction to inclusive track pt spec
4055 AddCorrectionLevelSinglePt("secCorr");
4057 TH1F* histPt = fCorrSinglePt[fNCorrectionLevelsSinglePt-2]->GetTrackPt(0);
4059 if(histPt->GetNbinsX() != fh1SecCorrSinglePt->GetNbinsX())
4060 Printf("%s:%d: inconsistency pt spec and secondaries corr bins - rebin effCorr ...", (char*)__FILE__,__LINE__);
4062 TString histNamePt = histPt->GetName();
4063 TH1F* hTrackPtSecCorr = (TH1F*) histPt->Clone(histNamePt);
4065 hTrackPtSecCorr->Multiply(histPt,fh1SecCorrSinglePt,1,1,"");
4067 fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->AddCorrHistos(0,hTrackPtSecCorr);