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"
32 #include <iostream> // OB TEST!!!
39 ///////////////////////////////////////////////////////////////////////////////
41 // correction class for ouput of AliAnalysisTaskFragmentationFunction //
43 // corrections for: reconstruction efficiency, momentum resolution, //
44 // secondaries, UE / HI background, fluctuations //
45 // back-correction on jet energy on dN/dxi //
47 // read MC ouput and write out efficiency histos, response matrices etc. //
48 // read measured distributions and apply efficiency, response matrices etc. //
50 // contact: o.busch@gsi.de //
52 ///////////////////////////////////////////////////////////////////////////////
55 ClassImp(AliFragmentationFunctionCorrections)
57 //________________________________________________________________________
58 AliFragmentationFunctionCorrections::AliFragmentationFunctionCorrections()
65 ,fNHistoBinsSinglePt(0)
66 ,fHistoBinsSinglePt(0)
73 ,fNCorrectionLevels(0)
75 ,fNCorrectionLevelsBgr(0)
77 ,fNCorrectionLevelsSinglePt(0)
86 ,fh1BbBCorrSinglePt(0)
102 ,fh1FFTrackPtBackFolded(0)
104 ,fh1FFXiBackFolded(0)
105 ,fh1FFRatioTrackPtFolded(0)
106 ,fh1FFRatioZFolded(0)
107 ,fh1FFRatioXiFolded(0)
108 ,fh1FFRatioTrackPtBackFolded(0)
109 ,fh1FFRatioZBackFolded(0)
110 ,fh1FFRatioXiBackFolded(0)
111 ,fh1SingleTrackPtBackFolded(0)
112 ,fh1RatioSingleTrackPtFolded(0)
113 ,fh1RatioSingleTrackPtBackFolded(0)
114 ,fhnResponseSinglePt(0)
118 ,fh1FFTrackPtPrior(0)
121 ,fh1SecCorrSinglePt(0)
123 // default constructor
126 //________________________________________________________________________________________________________________________
127 AliFragmentationFunctionCorrections::AliFragmentationFunctionCorrections(const AliFragmentationFunctionCorrections ©)
130 ,fNJetPtSlices(copy.fNJetPtSlices)
131 ,fJetPtSlices(copy.fJetPtSlices)
133 ,fNJetsBgr(copy.fNJetsBgr)
134 ,fNHistoBinsSinglePt(copy.fNHistoBinsSinglePt)
135 ,fHistoBinsSinglePt(copy.fHistoBinsSinglePt)
136 ,fNHistoBinsPt(copy.fNHistoBinsPt)
137 ,fNHistoBinsZ(copy.fNHistoBinsZ)
138 ,fNHistoBinsXi(copy.fNHistoBinsXi)
139 ,fHistoBinsPt(copy.fHistoBinsPt)
140 ,fHistoBinsZ(copy.fHistoBinsZ)
141 ,fHistoBinsXi(copy.fHistoBinsXi)
142 ,fNCorrectionLevels(copy.fNCorrectionLevels)
143 ,fCorrFF(copy.fCorrFF)
144 ,fNCorrectionLevelsBgr(copy.fNCorrectionLevelsBgr)
145 ,fCorrBgr(copy.fCorrBgr)
146 ,fNCorrectionLevelsSinglePt(copy.fNCorrectionLevelsSinglePt)
147 ,fCorrSinglePt(copy.fCorrSinglePt)
148 ,fh1EffSinglePt(copy.fh1EffSinglePt)
149 ,fh1EffPt(copy.fh1EffPt)
150 ,fh1EffZ(copy.fh1EffZ)
151 ,fh1EffXi(copy.fh1EffXi)
152 ,fh1EffBgrPt(copy.fh1EffBgrPt)
153 ,fh1EffBgrZ(copy.fh1EffBgrZ)
154 ,fh1EffBgrXi(copy.fh1EffBgrXi)
155 ,fh1BbBCorrSinglePt(copy.fh1BbBCorrSinglePt)
156 ,fh1BbBPt(copy.fh1BbBPt)
157 ,fh1BbBZ(copy.fh1BbBZ)
158 ,fh1BbBXi(copy.fh1BbBXi)
159 ,fh1BbBBgrPt(copy.fh1BbBBgrPt)
160 ,fh1BbBBgrZ(copy.fh1BbBBgrZ)
161 ,fh1BbBBgrXi(copy.fh1BbBBgrXi)
162 ,fh1FoldingCorrPt(copy.fh1FoldingCorrPt)
163 ,fh1FoldingCorrZ(copy.fh1FoldingCorrZ)
164 ,fh1FoldingCorrXi(copy.fh1FoldingCorrXi)
165 ,fh1SecCorrPt(copy.fh1SecCorrPt)
166 ,fh1SecCorrZ(copy.fh1SecCorrZ)
167 ,fh1SecCorrXi(copy.fh1SecCorrXi)
168 ,fh1SecCorrBgrPt(copy.fh1SecCorrBgrPt)
169 ,fh1SecCorrBgrZ(copy.fh1SecCorrBgrZ)
170 ,fh1SecCorrBgrXi(copy.fh1SecCorrBgrXi)
171 ,fh1FFTrackPtBackFolded(copy.fh1FFTrackPtBackFolded)
172 ,fh1FFZBackFolded(copy.fh1FFZBackFolded)
173 ,fh1FFXiBackFolded(copy.fh1FFXiBackFolded)
174 ,fh1FFRatioTrackPtFolded(copy.fh1FFRatioTrackPtFolded)
175 ,fh1FFRatioZFolded(copy.fh1FFRatioZFolded)
176 ,fh1FFRatioXiFolded(copy.fh1FFRatioXiFolded)
177 ,fh1FFRatioTrackPtBackFolded(copy.fh1FFRatioTrackPtBackFolded)
178 ,fh1FFRatioZBackFolded(copy.fh1FFRatioZBackFolded)
179 ,fh1FFRatioXiBackFolded(copy.fh1FFRatioXiBackFolded)
180 ,fh1SingleTrackPtBackFolded(copy.fh1SingleTrackPtBackFolded)
181 ,fh1RatioSingleTrackPtFolded(copy.fh1RatioSingleTrackPtFolded)
182 ,fh1RatioSingleTrackPtBackFolded(copy.fh1RatioSingleTrackPtBackFolded)
183 ,fhnResponseSinglePt(copy.fhnResponseSinglePt)
184 ,fhnResponsePt(copy.fhnResponsePt)
185 ,fhnResponseZ(copy.fhnResponseZ)
186 ,fhnResponseXi(copy.fhnResponseXi)
187 ,fh1FFTrackPtPrior(copy.fh1FFTrackPtPrior)
188 ,fh1FFZPrior(copy.fh1FFZPrior)
189 ,fh1FFXiPrior(copy.fh1FFXiPrior)
190 ,fh1SecCorrSinglePt(copy.fh1SecCorrSinglePt)
196 // ______________________________________________________________________________________________________________________________
197 AliFragmentationFunctionCorrections& AliFragmentationFunctionCorrections::operator=(const AliFragmentationFunctionCorrections& o)
202 TObject::operator=(o);
204 fNJetPtSlices = o.fNJetPtSlices;
205 fJetPtSlices = o.fJetPtSlices;
207 fNJetsBgr = o.fNJetsBgr;
208 fNHistoBinsSinglePt = o.fNHistoBinsSinglePt;
209 fHistoBinsSinglePt = o.fHistoBinsSinglePt;
210 fNHistoBinsPt = o.fNHistoBinsPt;
211 fNHistoBinsZ = o.fNHistoBinsZ;
212 fNHistoBinsXi = o.fNHistoBinsXi;
213 fHistoBinsPt = o.fHistoBinsPt;
214 fHistoBinsZ = o.fHistoBinsZ;
215 fHistoBinsXi = o.fHistoBinsXi;
216 fNCorrectionLevels = o.fNCorrectionLevels;
218 fNCorrectionLevelsBgr = o.fNCorrectionLevelsBgr;
219 fCorrBgr = o.fCorrBgr;
220 fNCorrectionLevelsSinglePt = o.fNCorrectionLevelsSinglePt;
221 fCorrSinglePt = o.fCorrSinglePt;
222 fh1EffSinglePt = o.fh1EffSinglePt;
223 fh1EffPt = o.fh1EffPt;
225 fh1EffXi = o.fh1EffXi;
226 fh1EffBgrPt = o.fh1EffBgrPt;
227 fh1EffBgrZ = o.fh1EffBgrZ;
228 fh1EffBgrXi = o.fh1EffBgrXi;
229 fh1BbBCorrSinglePt = o.fh1BbBCorrSinglePt;
230 fh1BbBPt = o.fh1BbBPt;
232 fh1BbBXi = o.fh1BbBXi;
233 fh1BbBBgrPt = o.fh1BbBBgrPt;
234 fh1BbBBgrZ = o.fh1BbBBgrZ;
235 fh1BbBBgrXi = o.fh1BbBBgrXi;
236 fh1FoldingCorrPt = o.fh1FoldingCorrPt;
237 fh1FoldingCorrZ = o.fh1FoldingCorrZ;
238 fh1FoldingCorrXi = o.fh1FoldingCorrXi;
239 fh1SecCorrPt = o.fh1SecCorrPt;
240 fh1SecCorrZ = o.fh1SecCorrZ;
241 fh1SecCorrXi = o.fh1SecCorrXi;
242 fh1SecCorrBgrPt = o.fh1SecCorrBgrPt;
243 fh1SecCorrBgrZ = o.fh1SecCorrBgrZ;
244 fh1SecCorrBgrXi = o.fh1SecCorrBgrXi;
245 fh1FFTrackPtBackFolded = o.fh1FFTrackPtBackFolded;
246 fh1FFZBackFolded = o.fh1FFZBackFolded;
247 fh1FFXiBackFolded = o.fh1FFXiBackFolded;
248 fh1FFRatioTrackPtFolded = o.fh1FFRatioTrackPtFolded;
249 fh1FFRatioZFolded = o.fh1FFRatioZFolded;
250 fh1FFRatioXiFolded = o.fh1FFRatioXiFolded;
251 fh1FFRatioTrackPtBackFolded = o.fh1FFRatioTrackPtBackFolded;
252 fh1FFRatioZBackFolded = o.fh1FFRatioZBackFolded;
253 fh1FFRatioXiBackFolded = o.fh1FFRatioXiBackFolded;
254 fh1SingleTrackPtBackFolded = o.fh1SingleTrackPtBackFolded;
255 fh1RatioSingleTrackPtFolded = o.fh1RatioSingleTrackPtFolded;
256 fh1RatioSingleTrackPtBackFolded = o.fh1RatioSingleTrackPtBackFolded;
257 fhnResponseSinglePt = o.fhnResponseSinglePt;
258 fhnResponsePt = o.fhnResponsePt;
259 fhnResponseZ = o.fhnResponseZ;
260 fhnResponseXi = o.fhnResponseXi;
261 fh1FFTrackPtPrior = o.fh1FFTrackPtPrior;
262 fh1FFZPrior = o.fh1FFZPrior;
263 fh1FFXiPrior = o.fh1FFXiPrior;
264 fh1SecCorrSinglePt = o.fh1SecCorrSinglePt;
270 //_________________________________________________________________________
271 AliFragmentationFunctionCorrections::~AliFragmentationFunctionCorrections()
275 if(fJetPtSlices) delete fJetPtSlices;
276 if(fNJets) delete fNJets;
277 if(fNJetsBgr) delete fNJetsBgr;
279 DeleteHistoArray(fh1EffPt);
280 DeleteHistoArray(fh1EffXi);
281 DeleteHistoArray(fh1EffZ );
283 DeleteHistoArray(fh1EffBgrPt);
284 DeleteHistoArray(fh1EffBgrXi);
285 DeleteHistoArray(fh1EffBgrZ);
288 DeleteHistoArray(fh1BbBPt);
289 DeleteHistoArray(fh1BbBXi);
290 DeleteHistoArray(fh1BbBZ);
292 DeleteHistoArray(fh1BbBBgrPt);
293 DeleteHistoArray(fh1BbBBgrXi);
294 DeleteHistoArray(fh1BbBBgrZ);
296 DeleteHistoArray(fh1FoldingCorrPt);
297 DeleteHistoArray(fh1FoldingCorrXi);
298 DeleteHistoArray(fh1FoldingCorrZ);
301 DeleteHistoArray(fh1SecCorrPt);
302 DeleteHistoArray(fh1SecCorrXi);
303 DeleteHistoArray(fh1SecCorrZ);
305 DeleteHistoArray(fh1SecCorrBgrPt);
306 DeleteHistoArray(fh1SecCorrBgrXi);
307 DeleteHistoArray(fh1SecCorrBgrZ);
311 DeleteHistoArray(fh1FFTrackPtBackFolded);
312 DeleteHistoArray(fh1FFZBackFolded);
313 DeleteHistoArray(fh1FFXiBackFolded);
315 DeleteHistoArray(fh1FFRatioTrackPtFolded);
316 DeleteHistoArray(fh1FFRatioZFolded);
317 DeleteHistoArray(fh1FFRatioXiFolded);
319 DeleteHistoArray(fh1FFRatioTrackPtBackFolded);
320 DeleteHistoArray(fh1FFRatioZBackFolded);
321 DeleteHistoArray(fh1FFRatioXiBackFolded);
323 DeleteTHnSparseArray(fhnResponsePt);
324 DeleteTHnSparseArray(fhnResponseZ);
325 DeleteTHnSparseArray(fhnResponseXi);
327 DeleteHistoArray(fh1FFTrackPtPrior);
328 DeleteHistoArray(fh1FFZPrior);
329 DeleteHistoArray(fh1FFXiPrior);
331 // clean up corrected FF
333 for(Int_t c=0; c<fNCorrectionLevels; c++) delete fCorrFF[c];
338 for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) delete fCorrBgr[c];
341 // clean up inclusive pt
342 for(Int_t c=0; c<fNCorrectionLevelsSinglePt; c++) delete fCorrSinglePt[c];
343 delete[] fCorrSinglePt;
345 delete[] fNHistoBinsPt;
346 delete[] fNHistoBinsZ;
347 delete[] fNHistoBinsXi;
349 // clean up histo bins
351 if(fHistoBinsSinglePt) delete fHistoBinsSinglePt;
353 for(Int_t i=0; i<fNJetPtSlices; i++) delete fHistoBinsPt[i];
354 for(Int_t i=0; i<fNJetPtSlices; i++) delete fHistoBinsZ[i];
355 for(Int_t i=0; i<fNJetPtSlices; i++) delete fHistoBinsXi[i];
357 delete[] fHistoBinsPt;
358 delete[] fHistoBinsZ;
359 delete[] fHistoBinsXi;
362 //_________________________________________________________________________________
363 AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::AliFragFuncCorrHistos()
371 // default constructor
375 //__________________________________________________________________________________________________________________
376 AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::AliFragFuncCorrHistos(const AliFragmentationFunctionCorrections::AliFragFuncCorrHistos& copy)
378 ,fArraySize(copy.fArraySize)
382 ,fCorrLabel(copy.fCorrLabel)
386 fh1CorrFFTrackPt = new TH1F*[copy.fArraySize];
387 fh1CorrFFZ = new TH1F*[copy.fArraySize];
388 fh1CorrFFXi = new TH1F*[copy.fArraySize];
390 for(Int_t i=0; i<copy.fArraySize; i++){
391 fh1CorrFFTrackPt[i] = new TH1F(*(copy.fh1CorrFFTrackPt[i]));
392 fh1CorrFFZ[i] = new TH1F(*(copy.fh1CorrFFZ[i]));
393 fh1CorrFFXi[i] = new TH1F(*(copy.fh1CorrFFXi[i]));
397 //_______________________________________________________________________________________________________________________________________________________________
398 AliFragmentationFunctionCorrections::AliFragFuncCorrHistos& AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::operator=(const AliFragmentationFunctionCorrections::AliFragFuncCorrHistos& o)
403 TObject::operator=(o);
404 Int_t fArraySize_new = o.fArraySize;
405 fCorrLabel = o.fCorrLabel;
407 TH1F** fh1CorrFFTrackPt_new = new TH1F*[fArraySize_new];
408 TH1F** fh1CorrFFZ_new = new TH1F*[fArraySize_new];
409 TH1F** fh1CorrFFXi_new = new TH1F*[fArraySize_new];
411 for(Int_t i=0; i<fArraySize_new; i++){
412 fh1CorrFFTrackPt_new[i] = new TH1F(*(o.fh1CorrFFTrackPt[i]));
413 fh1CorrFFZ_new[i] = new TH1F(*(o.fh1CorrFFZ[i]));
414 fh1CorrFFXi_new[i] = new TH1F(*(o.fh1CorrFFXi[i]));
420 for(int i=0; i<fArraySize; i++) delete fh1CorrFFTrackPt[i];
421 for(int i=0; i<fArraySize; i++) delete fh1CorrFFZ[i];
422 for(int i=0; i<fArraySize; i++) delete fh1CorrFFXi[i];
425 if(fh1CorrFFTrackPt) delete[] fh1CorrFFTrackPt;
426 if(fh1CorrFFZ) delete[] fh1CorrFFZ;
427 if(fh1CorrFFXi) delete[] fh1CorrFFXi;
431 fArraySize = fArraySize_new;
433 fh1CorrFFTrackPt = fh1CorrFFTrackPt_new;
434 fh1CorrFFZ = fh1CorrFFZ_new;
435 fh1CorrFFXi = fh1CorrFFXi_new;
437 for(Int_t i=0; i<fArraySize; i++){
438 fh1CorrFFTrackPt[i] = fh1CorrFFTrackPt_new[i];
439 fh1CorrFFZ[i] = fh1CorrFFZ_new[i];
440 fh1CorrFFXi[i] = fh1CorrFFXi_new[i];
447 //__________________________________________________________________________________
448 AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::~AliFragFuncCorrHistos()
453 for(int i=0; i<fArraySize; i++) delete fh1CorrFFTrackPt[i];
454 for(int i=0; i<fArraySize; i++) delete fh1CorrFFZ[i];
455 for(int i=0; i<fArraySize; i++) delete fh1CorrFFXi[i];
458 if(fh1CorrFFTrackPt) delete[] fh1CorrFFTrackPt;
459 if(fh1CorrFFZ) delete[] fh1CorrFFZ;
460 if(fh1CorrFFXi) delete[] fh1CorrFFXi;
464 //___________________________________________________________________________________________________________________
465 AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::AliFragFuncCorrHistos(const char* label, Int_t arraySize)
475 fArraySize = arraySize;
476 fh1CorrFFTrackPt = new TH1F*[arraySize];
477 fh1CorrFFZ = new TH1F*[arraySize];
478 fh1CorrFFXi = new TH1F*[arraySize];
480 for(int i=0; i<arraySize; i++) fh1CorrFFTrackPt[i] = new TH1F;
481 for(int i=0; i<arraySize; i++) fh1CorrFFZ[i] = new TH1F;
482 for(int i=0; i<arraySize; i++) fh1CorrFFXi[i] = new TH1F;
485 //_______________________________________________________________________________________________________________________________
486 void AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::AddCorrHistos(Int_t slice,TH1F* histPt,TH1F* histZ,TH1F* histXi)
488 // place histo in array, add corrLabel to histo name
490 if(slice>=fArraySize){
491 Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
496 TString name = histPt->GetName();
497 if(fCorrLabel.Length()>0) name += "_"+fCorrLabel;
498 histPt->SetName(name);
500 if(!(histPt->GetSumw2())) histPt->Sumw2();
502 new(fh1CorrFFTrackPt[slice]) TH1F(*histPt);
506 TString name = histZ->GetName();
507 if(fCorrLabel.Length()>0) name += "_"+fCorrLabel;
508 histZ->SetName(name);
510 if(!(histZ->GetSumw2())) histZ->Sumw2();
512 new(fh1CorrFFZ[slice]) TH1F(*histZ);
516 TString name = histXi->GetName();
517 if(fCorrLabel.Length()>0) name += "_"+fCorrLabel;
518 histXi->SetName(name);
520 if(!(histXi->GetSumw2())) histXi->Sumw2();
522 new(fh1CorrFFXi[slice]) TH1F(*histXi);
526 //___________________________________________________________________________________________________________________________________
527 void AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::ReplaceCorrHistos(Int_t slice,TH1F* histPt,TH1F* histZ,TH1F* histXi)
529 // replace histo in array
531 if(slice>=fArraySize){
532 Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
537 if(!(histPt->GetSumw2())) histPt->Sumw2();
539 TString name = histPt->GetName();
540 histPt->SetName(name);
542 delete fh1CorrFFTrackPt[slice];
543 fh1CorrFFTrackPt[slice] = new TH1F;
544 new(fh1CorrFFTrackPt[slice]) TH1F(*histPt);
548 if(!(histZ->GetSumw2())) histZ->Sumw2();
550 TString name = histZ->GetName();
551 histZ->SetName(name);
553 delete fh1CorrFFZ[slice];
554 fh1CorrFFZ[slice] = new TH1F;
555 new(fh1CorrFFZ[slice]) TH1F(*histZ);
559 if(!(histXi->GetSumw2())) histXi->Sumw2();
561 TString name = histXi->GetName();
562 histXi->SetName(name);
564 delete fh1CorrFFXi[slice];
565 fh1CorrFFXi[slice] = new TH1F;
566 new(fh1CorrFFXi[slice]) TH1F(*histXi);
570 // ___________________________________________________________________________________________
571 TH1F* AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::GetTrackPt(const Int_t slice)
575 if(slice>=fArraySize){
576 Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
580 return fh1CorrFFTrackPt[slice];
584 // ______________________________________________________________________________________
585 TH1F* AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::GetZ(const Int_t slice)
589 if(slice>=fArraySize){
590 Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
594 return fh1CorrFFZ[slice];
597 // ________________________________________________________________________________________
598 TH1F* AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::GetXi(const Int_t slice)
602 if(slice>=fArraySize){
603 Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
607 return fh1CorrFFXi[slice];
610 // __________________________________________________________________________
611 void AliFragmentationFunctionCorrections::DeleteHistoArray(TH1F** hist) const
613 // delete array of TH1
616 for(Int_t i=0; i<fNJetPtSlices; i++) delete hist[i];
621 // ____________________________________________________________________________________
622 void AliFragmentationFunctionCorrections::DeleteTHnSparseArray(THnSparse** hist) const
624 // delete array of THnSparse
627 for(Int_t i=0; i<fNJetPtSlices; i++) delete hist[i];
632 // _________________________________________________________
633 TH1F** AliFragmentationFunctionCorrections::BookHistoArray()
638 Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
642 TH1F** hist = new TH1F*[fNJetPtSlices];
643 for(Int_t i=0; i<fNJetPtSlices; i++) hist[i] = new TH1F();
648 //__________________________________________________________________
649 THnSparse** AliFragmentationFunctionCorrections::BookTHnSparseArray()
651 // book array of THnSparse
654 Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
658 THnSparse** hist = new THnSparse*[fNJetPtSlices];
659 for(Int_t i=0; i<fNJetPtSlices; i++) hist[i] = new THnSparseF();
664 //_____________________________________________________________________________
665 void AliFragmentationFunctionCorrections::AddCorrectionLevel(const char* label)
667 // increase corr level
670 Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
674 if(fNCorrectionLevels >= fgMaxNCorrectionLevels){
675 Printf("%s:%d -- max correction level exceeded", (char*)__FILE__,__LINE__);
679 fCorrFF[fNCorrectionLevels] = new AliFragFuncCorrHistos(label,fNJetPtSlices);
680 fNCorrectionLevels++;
684 //________________________________________________________________________________
685 void AliFragmentationFunctionCorrections::AddCorrectionLevelBgr(const char* label)
687 // increase corr level for bgr FF
690 if(fDebug>0) Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
694 if(fNCorrectionLevelsBgr >= fgMaxNCorrectionLevels){
695 Printf("%s:%d -- max correction level exceeded", (char*)__FILE__,__LINE__);
699 fCorrBgr[fNCorrectionLevelsBgr] = new AliFragFuncCorrHistos(label,fNJetPtSlices);
700 fNCorrectionLevelsBgr++;
703 //_____________________________________________________________________________________
704 void AliFragmentationFunctionCorrections::AddCorrectionLevelSinglePt(const char* label)
706 // increase corr level single pt spec
708 Int_t nJetPtSlicesSingle = 1; // pro forma
710 if(fNCorrectionLevelsSinglePt >= fgMaxNCorrectionLevels){
711 Printf("%s:%d -- max correction level exceeded", (char*)__FILE__,__LINE__);
715 fCorrSinglePt[fNCorrectionLevelsSinglePt] = new AliFragFuncCorrHistos(label,nJetPtSlicesSingle);
716 fNCorrectionLevelsSinglePt++;
719 //_____________________________________________________________________________________________
720 void AliFragmentationFunctionCorrections::SetJetPtSlices(Float_t* bins, const Int_t nJetPtSlices)
723 // once slices are known, can book all other histos
725 fNJetPtSlices = nJetPtSlices;
727 const Int_t size = nJetPtSlices+1;
728 fJetPtSlices = new TArrayF(size,bins);
732 fNJets = new TArrayF(size);
735 fNJetsBgr = new TArrayF(size);
740 fNCorrectionLevels = 0;
741 fCorrFF = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels];
742 AddCorrectionLevel(); // first 'correction' level = raw FF
744 fNCorrectionLevelsBgr = 0;
745 fCorrBgr = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels];
746 AddCorrectionLevelBgr(); // first 'correction' level = raw bgr dist
750 fh1EffPt = BookHistoArray();
751 fh1EffXi = BookHistoArray();
752 fh1EffZ = BookHistoArray();
754 fh1EffBgrPt = BookHistoArray();
755 fh1EffBgrXi = BookHistoArray();
756 fh1EffBgrZ = BookHistoArray();
761 fh1BbBPt = BookHistoArray();
762 fh1BbBXi = BookHistoArray();
763 fh1BbBZ = BookHistoArray();
765 fh1BbBBgrPt = BookHistoArray();
766 fh1BbBBgrXi = BookHistoArray();
767 fh1BbBBgrZ = BookHistoArray();
769 fh1FoldingCorrPt = BookHistoArray();
770 fh1FoldingCorrXi = BookHistoArray();
771 fh1FoldingCorrZ = BookHistoArray();
775 fh1SecCorrPt = BookHistoArray();
776 fh1SecCorrXi = BookHistoArray();
777 fh1SecCorrZ = BookHistoArray();
779 fh1SecCorrBgrPt = BookHistoArray();
780 fh1SecCorrBgrXi = BookHistoArray();
781 fh1SecCorrBgrZ = BookHistoArray();
785 fh1FFTrackPtBackFolded = BookHistoArray();
786 fh1FFXiBackFolded = BookHistoArray();
787 fh1FFZBackFolded = BookHistoArray();
789 fh1FFRatioTrackPtFolded = BookHistoArray();
790 fh1FFRatioXiFolded = BookHistoArray();
791 fh1FFRatioZFolded = BookHistoArray();
793 fh1FFRatioTrackPtBackFolded = BookHistoArray();
794 fh1FFRatioXiBackFolded = BookHistoArray();
795 fh1FFRatioZBackFolded = BookHistoArray();
799 fhnResponsePt = BookTHnSparseArray();
800 fhnResponseZ = BookTHnSparseArray();
801 fhnResponseXi = BookTHnSparseArray();
803 fh1FFTrackPtPrior = BookHistoArray();
804 fh1FFZPrior = BookHistoArray();
805 fh1FFXiPrior = BookHistoArray();
809 fNHistoBinsPt = new Int_t[fNJetPtSlices];
810 fNHistoBinsXi = new Int_t[fNJetPtSlices];
811 fNHistoBinsZ = new Int_t[fNJetPtSlices];
813 for(Int_t i=0; i<fNJetPtSlices; i++) fNHistoBinsPt[i] = 0;
814 for(Int_t i=0; i<fNJetPtSlices; i++) fNHistoBinsXi[i] = 0;
815 for(Int_t i=0; i<fNJetPtSlices; i++) fNHistoBinsZ[i] = 0;
817 fHistoBinsPt = new TArrayD*[fNJetPtSlices];
818 fHistoBinsXi = new TArrayD*[fNJetPtSlices];
819 fHistoBinsZ = new TArrayD*[fNJetPtSlices];
821 for(Int_t i=0; i<fNJetPtSlices; i++) fHistoBinsPt[i] = new TArrayD(0);
822 for(Int_t i=0; i<fNJetPtSlices; i++) fHistoBinsXi[i] = new TArrayD(0);
823 for(Int_t i=0; i<fNJetPtSlices; i++) fHistoBinsZ[i] = new TArrayD(0);
826 //_____________________________________________________________________________________________________________________________________
827 void AliFragmentationFunctionCorrections::SetHistoBins(const Int_t jetPtSlice, const Int_t sizeBins, Double_t* bins, const Int_t type)
829 // set histo bins for jet pt slice
830 // if binning undefined for any slice, original binning will be used
833 Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
837 if(jetPtSlice>=fNJetPtSlices){
838 Printf("%s:%d -- jetPtSlice %d exceeds max",(char*)__FILE__,__LINE__,jetPtSlice);
843 fNHistoBinsPt[jetPtSlice] = sizeBins-1;
844 fHistoBinsPt[jetPtSlice]->Set(sizeBins,bins);
846 else if(type==kFlagZ){
847 fNHistoBinsZ[jetPtSlice] = sizeBins-1;
848 fHistoBinsZ[jetPtSlice]->Set(sizeBins,bins);
851 else if(type==kFlagXi){
852 fNHistoBinsXi[jetPtSlice] = sizeBins-1;
853 fHistoBinsXi[jetPtSlice]->Set(sizeBins,bins);
857 //__________________________________________________________________________________________________________________________________________________________________
858 void AliFragmentationFunctionCorrections::SetHistoBins(const Int_t jetPtSlice, const Int_t nBinsLimits, Double_t* binsLimits, Double_t* binsWidth, const Int_t type)
860 // set histo bins for jet pt slice
861 // function takes array of limits and widths (e.g. 1 GeV bins up to 10 GeV, 2 GeV width up to 20 GeV, ...)
862 // array size of binsLimits: nBinsLimits
863 // array size of binsWidth: nBinsLimits-1
864 // binsLimits have to be in increasing order
865 // if binning undefined for any slice, original binning will be used
868 Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
872 if(jetPtSlice>=fNJetPtSlices){
873 Printf("%s:%d -- jetPtSlice %d exceeds max",(char*)__FILE__,__LINE__,jetPtSlice);
878 Double_t binLimitMin = binsLimits[0];
879 Double_t binLimitMax = binsLimits[nBinsLimits-1];
881 Double_t binLimit = binLimitMin; // start value
883 Int_t sizeUpperLim = 10000; //static_cast<Int_t>(binLimitMax/binsWidth[0])+1;
884 TArrayD binsArray(sizeUpperLim);
886 binsArray.SetAt(binLimitMin,nBins++);
888 while(binLimit<0.999999*binLimitMax && nBins<sizeUpperLim){ // 0.999 numerical safety factor
890 Int_t currentSlice = -1;
891 for(Int_t i=0; i<nBinsLimits; i++){
892 if(binLimit >= 0.999999*binsLimits[i]) currentSlice = i;
895 Double_t currentBinWidth = binsWidth[currentSlice];
896 binLimit += currentBinWidth;
898 // if(type == kFlagZ) std::cout<<" SetHistoBins: nBins "<<nBins<<" binLimit "<<binLimit<<" binLimitMax "<<binLimitMax
899 // <<" currentSlice "<<currentSlice<<std::endl;
901 binsArray.SetAt(binLimit,nBins++);
904 Double_t* bins = binsArray.GetArray();
906 SetHistoBins(jetPtSlice,nBins,bins,type);
909 //_____________________________________________________________________________________________________________________________________
910 TArrayD* AliFragmentationFunctionCorrections::GetHistoBins(const Int_t jetPtSlice, const Int_t type)
912 // set histo bins for jet pt slice
913 // if binning undefined for any slice, original binning will be used
916 Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
920 if(jetPtSlice>=fNJetPtSlices){
921 Printf("%s:%d -- jetPtSlice %d exceeds max",(char*)__FILE__,__LINE__,jetPtSlice);
926 return fHistoBinsPt[jetPtSlice];
928 else if(type==kFlagZ){
929 return fHistoBinsZ[jetPtSlice];
932 else if(type==kFlagXi){
933 return fHistoBinsXi[jetPtSlice];
936 Printf("%s%d unknown type",(char*)__FILE__,__LINE__);
941 //__________________________________________________________________________________________________
942 void AliFragmentationFunctionCorrections::SetHistoBinsSinglePt(const Int_t sizeBins, Double_t* bins)
944 // set histo bins for inclusive pt spectra
945 // if binning undefined, original binning will be used
947 fNHistoBinsSinglePt = sizeBins-1;
949 fHistoBinsSinglePt = new TArrayD(sizeBins);
950 fHistoBinsSinglePt->Set(sizeBins,bins);
953 //__________________________________________________________________________________________________________________________________________________________________
954 void AliFragmentationFunctionCorrections::SetHistoBinsSinglePt(const Int_t nBinsLimits, Double_t* binsLimits, Double_t* binsWidth)
956 // set histo bins for inclusive pt spectra
957 // function takes array of limits and widths (e.g. 1 GeV bins up to 10 GeV, 2 GeV width up to 20 GeV, ...)
958 // array size of binsLimits: nBinsLimits
959 // array size of binsWidth: nBinsLimits-1
960 // binsLimits have to be in increasing order
961 // if binning undefined for any slice, original binning will be used
964 Double_t binLimitMin = binsLimits[0];
965 Double_t binLimitMax = binsLimits[nBinsLimits-1];
967 Double_t binLimit = binLimitMin; // start value
969 Int_t sizeUpperLim = 10000; //static_cast<Int_t>(binLimitMax/binsWidth[0])+1;
970 TArrayD binsArray(sizeUpperLim);
972 binsArray.SetAt(binLimitMin,nBins++);
974 while(binLimit<binLimitMax && nBins<sizeUpperLim){
976 Int_t currentSlice = -1;
977 for(Int_t i=0; i<nBinsLimits; i++){
978 if(binLimit >= 0.999*binsLimits[i]) currentSlice = i; // 0.999 numerical saftey factor
981 Double_t currentBinWidth = binsWidth[currentSlice];
982 binLimit += currentBinWidth;
984 binsArray.SetAt(binLimit,nBins++);
987 Double_t* bins = binsArray.GetArray();
989 SetHistoBinsSinglePt(nBins,bins);
992 //____________________________________________________________________________________
993 void AliFragmentationFunctionCorrections::NormalizeTH1(TH1* hist, const Float_t nJets)
995 // FF normalization: divide by bin width and normalize to entries/jets
996 // should also work for TH2, but to be tested !!!
998 if(nJets == 0){ // nothing to do
999 if(fDebug>0) Printf("%s:%d -- normalize: nJets = 0, do nothing", (char*)__FILE__,__LINE__);
1003 Int_t nBins = hist->GetNbinsX()*hist->GetNbinsY()*hist->GetNbinsZ();
1005 for(Int_t bin=0; bin<=nBins; bin++){ // include under-/overflow (?)
1007 Double_t binWidth = hist->GetBinWidth(bin);
1008 Double_t binCont = hist->GetBinContent(bin);
1009 Double_t binErr = hist->GetBinError(bin);
1011 binCont /= binWidth;
1014 hist->SetBinContent(bin,binCont);
1015 hist->SetBinError(bin,binErr);
1018 hist->Scale(1/nJets);
1021 //_____________________________________________________
1022 void AliFragmentationFunctionCorrections::NormalizeFF()
1026 if(fNCorrectionLevels>1){
1027 Printf("%s:%d -- FF normalization should be done BEFORE any correction !!!", (char*)__FILE__,__LINE__);
1030 for(Int_t i=0; i<fNJetPtSlices; i++){
1032 if(fDebug>0) Printf(" normalizeFF: i %d, nJets %f",i,fNJets->At(i));
1034 NormalizeTH1(fCorrFF[0]->GetTrackPt(i),fNJets->At(i)); // always normalize corr level 0 = raw FF
1035 NormalizeTH1(fCorrFF[0]->GetZ(i),fNJets->At(i));
1036 NormalizeTH1(fCorrFF[0]->GetXi(i),fNJets->At(i));
1040 //______________________________________________________
1041 void AliFragmentationFunctionCorrections::NormalizeBgr()
1043 // normalize bgr/UE FF
1045 if(fNCorrectionLevelsBgr>1){
1046 Printf("%s:%d -- FF normalization should be done BEFORE any correction !!!", (char*)__FILE__,__LINE__);
1049 for(Int_t i=0; i<fNJetPtSlices; i++){
1050 NormalizeTH1(fCorrBgr[0]->GetTrackPt(i), fNJetsBgr->At(i)); // always normalize corr level 0 = raw FF
1051 NormalizeTH1(fCorrBgr[0]->GetZ(i), fNJetsBgr->At(i));
1052 NormalizeTH1(fCorrBgr[0]->GetXi(i),fNJetsBgr->At(i));
1057 //__________________________________________________________________________________________________
1058 void AliFragmentationFunctionCorrections::ReadRawFF(TString strfile, TString strID, TString strFFID)
1060 // read raw FF - standard dir/list name
1062 TString strdir = "PWGJE_FragmentationFunction_" + strID;
1063 TString strlist = "fracfunc_" + strID;
1065 ReadRawFF(strfile,strdir,strlist,strFFID);
1068 //____________________________________________________________________________________________________________________
1069 void AliFragmentationFunctionCorrections::ReadRawFF(TString strfile, TString strdir, TString strlist, TString strFFID)
1071 // get raw FF from input file, project in jet pt slice
1072 // normalization done separately
1074 TFile f(strfile,"READ");
1077 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1081 if(fDebug>0) Printf("%s:%d -- read FF from file %s, dir %s ",(char*)__FILE__,__LINE__,strfile.Data(), strdir.Data());
1083 gDirectory->cd(strdir);
1087 if(strlist && strlist.Length()){
1088 if(!(list = (TList*) gDirectory->Get(strlist))){
1089 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1094 TString hnameJetPt(Form("fh1FFJetPt%s",strFFID.Data()));
1095 TString hnameTrackPt(Form("fh2FFTrackPt%s",strFFID.Data()));
1096 TString hnameZ(Form("fh2FFZ%s",strFFID.Data()));
1097 TString hnameXi(Form("fh2FFXi%s",strFFID.Data()));
1099 TH1F* fh1FFJetPt = 0;
1100 TH2F* fh2FFTrackPt = 0;
1105 fh1FFJetPt = (TH1F*) list->FindObject(hnameJetPt);
1106 fh2FFTrackPt = (TH2F*) list->FindObject(hnameTrackPt);
1107 fh2FFZ = (TH2F*) list->FindObject(hnameZ);
1108 fh2FFXi = (TH2F*) list->FindObject(hnameXi);
1111 fh1FFJetPt = (TH1F*) gDirectory->Get(hnameJetPt);
1112 fh2FFTrackPt = (TH2F*) gDirectory->Get(hnameTrackPt);
1113 fh2FFZ = (TH2F*) gDirectory->Get(hnameZ);
1114 fh2FFXi = (TH2F*) gDirectory->Get(hnameXi);
1118 if(!fh1FFJetPt) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameJetPt.Data()); return; }
1119 if(!fh2FFTrackPt){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameTrackPt.Data()); return; }
1120 if(!fh2FFZ) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameZ.Data()); return; }
1121 if(!fh2FFXi) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameXi.Data()); return; }
1123 fh1FFJetPt->SetDirectory(0);
1124 fh2FFTrackPt->SetDirectory(0);
1125 fh2FFZ->SetDirectory(0);
1126 fh2FFXi->SetDirectory(0);
1133 for(Int_t i=0; i<fNJetPtSlices; i++){
1135 Float_t jetPtLoLim = fJetPtSlices->At(i);
1136 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1138 Int_t binLo = static_cast<Int_t>(fh1FFJetPt->FindBin(jetPtLoLim));
1139 Int_t binUp = static_cast<Int_t>(fh1FFJetPt->FindBin(jetPtUpLim)) - 1;
1141 Float_t nJetsBin = fh1FFJetPt->Integral(binLo,binUp);
1143 fNJets->SetAt(nJetsBin,i);
1145 if(fDebug>0) Printf("jet pt %d to %d: nJets %f",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),fNJets->At(i));
1150 for(Int_t i=0; i<fNJetPtSlices; i++){
1152 Float_t jetPtLoLim = fJetPtSlices->At(i);
1153 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1155 Int_t binLo = static_cast<Int_t>(fh2FFTrackPt->GetXaxis()->FindBin(jetPtLoLim));
1156 Int_t binUp = static_cast<Int_t>(fh2FFTrackPt->GetXaxis()->FindBin(jetPtUpLim))-1;
1158 if(binUp > fh2FFTrackPt->GetNbinsX()){
1159 Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
1163 TString strNameFFPt(Form("fh1FFTrackPt%s_%02d_%02d",strFFID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1164 TString strNameFFZ(Form("fh1FFZ%s_%02d_%02d",strFFID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1165 TString strNameFFXi(Form("fh1FFXi%s_%02d_%02d",strFFID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1167 // appendix 'unbinned' to avoid histos with same name after rebinning
1168 TH1F* projPt = (TH1F*) fh2FFTrackPt->ProjectionY(strNameFFPt+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
1169 TH1F* projZ = (TH1F*) fh2FFZ->ProjectionY(strNameFFZ+"_unBinned",binLo,binUp,"o");
1170 TH1F* projXi = (TH1F*) fh2FFXi->ProjectionY(strNameFFXi+"_unBinned",binLo,binUp,"o");
1172 if(fNHistoBinsPt[i]) projPt = (TH1F*) projPt->Rebin(fNHistoBinsPt[i],strNameFFPt,fHistoBinsPt[i]->GetArray());
1173 if(fNHistoBinsZ[i]) projZ = (TH1F*) projZ->Rebin(fNHistoBinsZ[i],strNameFFZ,fHistoBinsZ[i]->GetArray());
1174 if(fNHistoBinsXi[i]) projXi = (TH1F*) projXi->Rebin(fNHistoBinsXi[i],strNameFFXi,fHistoBinsXi[i]->GetArray());
1176 projPt->SetNameTitle(strNameFFPt,"");
1177 projZ->SetNameTitle(strNameFFZ,"");
1178 projXi->SetNameTitle(strNameFFXi,"");
1180 // raw FF = corr level 0
1181 fCorrFF[0]->AddCorrHistos(i,projPt,projZ,projXi);
1186 delete fh2FFTrackPt;
1191 //_____________________________________________________________________________________________________________________
1192 void AliFragmentationFunctionCorrections::ReadRawBgr(TString strfile, TString strID, TString strBgrID, TString strFFID)
1194 // read raw FF - standard dir/list name
1196 TString strdir = "PWGJE_FragmentationFunction_" + strID;
1197 TString strlist = "fracfunc_" + strID;
1199 ReadRawBgr(strfile,strdir,strlist,strBgrID,strFFID);
1202 //_______________________________________________________________________________________________________________________________________
1203 void AliFragmentationFunctionCorrections::ReadRawBgr(TString strfile, TString strdir, TString strlist, TString strBgrID, TString strFFID)
1205 // get raw FF from input file, project in jet pt slice
1206 // use jet dN/dpt corresponding to strFFID, bgr FF to strBgrID+strFFID
1207 // e.g. "fh1FFJetPtRecCuts", "fh2FFXiBgrPerpRecCuts"
1208 // normalization done separately
1210 TString strID = strBgrID + strFFID;
1212 TFile f(strfile,"READ");
1215 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1219 if(fDebug>0) Printf("%s:%d -- read Bgr %s from file %s, dir %s ",(char*)__FILE__,__LINE__,strBgrID.Data(),strfile.Data(),strdir.Data());
1221 gDirectory->cd(strdir);
1225 if(strlist && strlist.Length()){
1226 if(!(list = (TList*) gDirectory->Get(strlist))){
1227 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1232 TString hnameNJets = "fh1nRecJetsCuts";
1233 TString hnameJetPt(Form("fh1FFJetPt%s",strFFID.Data())); // not: strID.Data() !!! would not be proper normalization
1234 TString hnameBgrTrackPt(Form("fh2FFTrackPt%s",strID.Data()));
1235 TString hnameBgrZ(Form("fh2FFZ%s",strID.Data()));
1236 TString hnameBgrXi(Form("fh2FFXi%s",strID.Data()));
1239 TH1F* fh1FFJetPtBgr;
1240 TH2F* fh2FFTrackPtBgr;
1245 fh1NJets = (TH1F*) list->FindObject(hnameNJets); // needed for normalization of bgr out of 2 jets
1246 fh1FFJetPtBgr = (TH1F*) list->FindObject(hnameJetPt);
1247 fh2FFTrackPtBgr = (TH2F*) list->FindObject(hnameBgrTrackPt);
1248 fh2FFZBgr = (TH2F*) list->FindObject(hnameBgrZ);
1249 fh2FFXiBgr = (TH2F*) list->FindObject(hnameBgrXi);
1252 fh1NJets = (TH1F*) gDirectory->Get(hnameNJets); // needed for normalization of bgr out of 2 jets
1253 fh1FFJetPtBgr = (TH1F*) gDirectory->Get(hnameJetPt);
1254 fh2FFTrackPtBgr = (TH2F*) gDirectory->Get(hnameBgrTrackPt);
1255 fh2FFZBgr = (TH2F*) gDirectory->Get(hnameBgrZ);
1256 fh2FFXiBgr = (TH2F*) gDirectory->Get(hnameBgrXi);
1260 if(!fh1FFJetPtBgr) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameJetPt.Data()); return; }
1261 if(!fh1NJets) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameNJets.Data()); }
1262 if(!fh2FFTrackPtBgr){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrTrackPt.Data()); return; }
1263 if(!fh2FFZBgr) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrZ.Data()); return; }
1264 if(!fh2FFXiBgr) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrXi.Data()); return; }
1266 fh1FFJetPtBgr->SetDirectory(0);
1267 fh2FFTrackPtBgr->SetDirectory(0);
1268 fh2FFZBgr->SetDirectory(0);
1269 fh2FFXiBgr->SetDirectory(0);
1270 if(fh1NJets) fh1NJets->SetDirectory(0);
1276 for(Int_t i=0; i<fNJetPtSlices; i++){
1278 Float_t jetPtLoLim = fJetPtSlices->At(i);
1279 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1281 Int_t binLo = static_cast<Int_t>(fh1FFJetPtBgr->FindBin(jetPtLoLim));
1282 Int_t binUp = static_cast<Int_t>(fh1FFJetPtBgr->FindBin(jetPtUpLim)) - 1;
1284 Float_t nJetsBin = fh1FFJetPtBgr->Integral(binLo,binUp);
1285 Double_t scaleF = 1;
1287 //if(strBgrID.Contains("Out2Jets")){ // scale by ratio 2 jets events / all events
1288 // scaleF = fh1NJets->Integral(fh1NJets->FindBin(2),fh1NJets->GetNbinsX())
1289 // / fh1NJets->Integral(fh1NJets->FindBin(1),fh1NJets->GetNbinsX());}
1292 if(strBgrID.Contains("OutAllJets") ){ // scale by ratio >3 jets events / all events
1295 scaleF = fh1NJets->Integral(fh1NJets->FindBin(4),fh1NJets->GetNbinsX())
1296 / fh1NJets->Integral(fh1NJets->FindBin(1),fh1NJets->GetNbinsX());
1299 Printf("%s:%d -- use bgr OutAllJets, but histo fhn1NJets not found",(char*)__FILE__,__LINE__);
1304 fNJetsBgr->SetAt(nJetsBin*scaleF,i);
1306 if(fDebug>0) Printf("bgr jet pt %d to %d: nJets %f, scaleF %.2f",
1307 static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),nJetsBin,scaleF);
1313 for(Int_t i=0; i<fNJetPtSlices; i++){
1315 Float_t jetPtLoLim = fJetPtSlices->At(i);
1316 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1318 Int_t binLo = static_cast<Int_t>(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtLoLim));
1319 Int_t binUp = static_cast<Int_t>(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtUpLim))-1;
1321 if(binUp > fh2FFTrackPtBgr->GetNbinsX()){
1322 Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
1326 TString strNameBgrPt(Form("fh1BgrTrackPt%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1327 TString strNameBgrZ(Form("fh1BgrZ%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1328 TString strNameBgrXi(Form("fh1BgrXi%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1330 // appendix 'unbinned' to avoid histos with same name after rebinning
1331 TH1F* projPt = (TH1F*) fh2FFTrackPtBgr->ProjectionY(strNameBgrPt+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
1332 TH1F* projZ = (TH1F*) fh2FFZBgr->ProjectionY(strNameBgrZ+"_unBinned",binLo,binUp,"o");
1333 TH1F* projXi = (TH1F*) fh2FFXiBgr->ProjectionY(strNameBgrXi+"_unBinned",binLo,binUp,"o");
1335 if(fNHistoBinsPt[i]) projPt = (TH1F*) projPt->Rebin(fNHistoBinsPt[i],strNameBgrPt,fHistoBinsPt[i]->GetArray());
1336 if(fNHistoBinsZ[i]) projZ = (TH1F*) projZ->Rebin(fNHistoBinsZ[i],strNameBgrZ,fHistoBinsZ[i]->GetArray());
1337 if(fNHistoBinsXi[i]) projXi = (TH1F*) projXi->Rebin(fNHistoBinsXi[i],strNameBgrXi,fHistoBinsXi[i]->GetArray());
1339 projPt->SetNameTitle(strNameBgrPt,"");
1340 projZ->SetNameTitle(strNameBgrZ,"");
1341 projXi->SetNameTitle(strNameBgrXi,"");
1343 // raw bgr = corr level 0
1344 fCorrBgr[0]->AddCorrHistos(i,projPt,projZ,projXi);
1347 delete fh1FFJetPtBgr;
1348 if(fh1NJets) delete fh1NJets;
1349 delete fh2FFTrackPtBgr;
1354 //_____________________________________________________________________________________________________________________
1355 void AliFragmentationFunctionCorrections::ReadRawBgrEmbedding(TString strfile, TString strID, TString strFFID)
1357 // read raw FF - standard dir/list name
1359 TString strdir = "PWGJE_FragmentationFunction_" + strID;
1360 TString strlist = "fracfunc_" + strID;
1362 ReadRawBgrEmbedding(strfile,strdir,strlist,strFFID);
1365 //_______________________________________________________________________________________________________________________________________
1366 void AliFragmentationFunctionCorrections::ReadRawBgrEmbedding(TString strfile, TString strdir, TString strlist, TString strFFID)
1368 // get raw FF from input file, project in jet pt slice
1369 // for embedding, the bgr FF are taken from histos "fh1FFJetPtRecCuts", "fh2FFXiRecCuts"
1370 // normalization done separately
1372 TString strBgrID = "BckgEmbed";
1373 TString strID = strBgrID + strFFID;
1375 TFile f(strfile,"READ");
1378 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1382 if(fDebug>0) Printf("%s:%d -- read Bgr %s from file %s ",(char*)__FILE__,__LINE__,strFFID.Data(),strfile.Data());
1384 gDirectory->cd(strdir);
1388 if(!(list = (TList*) gDirectory->Get(strlist))){
1389 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1393 TString hnameNJets = "fh1nRecJetsCuts";
1394 TString hnameJetPt(Form("fh1FFJetPt%s",strFFID.Data()));
1395 TString hnameBgrTrackPt(Form("fh2FFTrackPt%s",strFFID.Data()));
1396 TString hnameBgrZ(Form("fh2FFZ%s",strFFID.Data()));
1397 TString hnameBgrXi(Form("fh2FFXi%s",strFFID.Data()));
1399 TH1F* fh1NJets = (TH1F*) list->FindObject(hnameNJets); // needed for normalization of bgr out of 2 jets
1400 TH1F* fh1FFJetPtBgr = (TH1F*) list->FindObject(hnameJetPt);
1401 TH2F* fh2FFTrackPtBgr = (TH2F*) list->FindObject(hnameBgrTrackPt);
1402 TH2F* fh2FFZBgr = (TH2F*) list->FindObject(hnameBgrZ);
1403 TH2F* fh2FFXiBgr = (TH2F*) list->FindObject(hnameBgrXi);
1405 if(!fh1FFJetPtBgr) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameJetPt.Data()); return; }
1406 if(!fh1NJets) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameNJets.Data()); return; }
1407 if(!fh2FFTrackPtBgr){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrTrackPt.Data()); return; }
1408 if(!fh2FFZBgr) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrZ.Data()); return; }
1409 if(!fh2FFXiBgr) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrXi.Data()); return; }
1411 fh1FFJetPtBgr->SetDirectory(0);
1412 fh1NJets->SetDirectory(0);
1413 fh2FFTrackPtBgr->SetDirectory(0);
1414 fh2FFZBgr->SetDirectory(0);
1415 fh2FFXiBgr->SetDirectory(0);
1421 for(Int_t i=0; i<fNJetPtSlices; i++){
1423 Float_t jetPtLoLim = fJetPtSlices->At(i);
1424 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1426 Int_t binLo = static_cast<Int_t>(fh1FFJetPtBgr->FindBin(jetPtLoLim));
1427 Int_t binUp = static_cast<Int_t>(fh1FFJetPtBgr->FindBin(jetPtUpLim)) - 1;
1429 Float_t nJetsBin = fh1FFJetPtBgr->Integral(binLo,binUp);
1430 Double_t scaleF = 1;
1432 fNJetsBgr->SetAt(nJetsBin*scaleF,i);
1434 if(fDebug>0) Printf("bgr jet pt %d to %d: nJets %f, scaleF %.2f",
1435 static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),nJetsBin,scaleF);
1441 for(Int_t i=0; i<fNJetPtSlices; i++){
1443 Float_t jetPtLoLim = fJetPtSlices->At(i);
1444 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1446 Int_t binLo = static_cast<Int_t>(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtLoLim));
1447 Int_t binUp = static_cast<Int_t>(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtUpLim))-1;
1449 if(binUp > fh2FFTrackPtBgr->GetNbinsX()){
1450 Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
1454 TString strNameBgrPt(Form("fh1BgrTrackPt%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1455 TString strNameBgrZ(Form("fh1BgrZ%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1456 TString strNameBgrXi(Form("fh1BgrXi%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1458 // appendix 'unbinned' to avoid histos with same name after rebinning
1459 TH1F* projPt = (TH1F*) fh2FFTrackPtBgr->ProjectionY(strNameBgrPt+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
1460 TH1F* projZ = (TH1F*) fh2FFZBgr->ProjectionY(strNameBgrZ+"_unBinned",binLo,binUp,"o");
1461 TH1F* projXi = (TH1F*) fh2FFXiBgr->ProjectionY(strNameBgrXi+"_unBinned",binLo,binUp,"o");
1463 if(fNHistoBinsPt[i]) projPt = (TH1F*) projPt->Rebin(fNHistoBinsPt[i],strNameBgrPt,fHistoBinsPt[i]->GetArray());
1464 if(fNHistoBinsZ[i]) projZ = (TH1F*) projZ->Rebin(fNHistoBinsZ[i],strNameBgrZ,fHistoBinsZ[i]->GetArray());
1465 if(fNHistoBinsXi[i]) projXi = (TH1F*) projXi->Rebin(fNHistoBinsXi[i],strNameBgrXi,fHistoBinsXi[i]->GetArray());
1467 projPt->SetNameTitle(strNameBgrPt,"");
1468 projZ->SetNameTitle(strNameBgrZ,"");
1469 projXi->SetNameTitle(strNameBgrXi,"");
1471 // raw bgr = corr level 0
1472 fCorrBgr[0]->AddCorrHistos(i,projPt,projZ,projXi);
1475 delete fh1FFJetPtBgr;
1477 delete fh2FFTrackPtBgr;
1483 //__________________________________________________________________________________________________________
1484 void AliFragmentationFunctionCorrections::WriteOutput(TString strfile, TString strdir, Bool_t updateOutfile)
1486 // write histos to file
1487 // skip histos with 0 entries
1489 TString outfileOption = "RECREATE";
1490 if(updateOutfile) outfileOption = "UPDATE";
1492 TFile f(strfile,outfileOption);
1495 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1499 if(fDebug>0) Printf("%s:%d -- write FF to file %s ",(char*)__FILE__,__LINE__,strfile.Data());
1501 if(strdir && strdir.Length()){
1502 TDirectory* dir = f.mkdir(strdir);
1506 for(Int_t i=0; i<fNJetPtSlices; i++){
1508 for(Int_t c=0; c<fNCorrectionLevels; c++) if(fCorrFF[c]->GetTrackPt(i)->GetEntries()) fCorrFF[c]->GetTrackPt(i)->Write();
1509 for(Int_t c=0; c<fNCorrectionLevels; c++) if(fCorrFF[c]->GetZ(i)->GetEntries()) fCorrFF[c]->GetZ(i)->Write();
1510 for(Int_t c=0; c<fNCorrectionLevels; c++) if(fCorrFF[c]->GetXi(i)->GetEntries()) fCorrFF[c]->GetXi(i)->Write();
1512 for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) if(fCorrBgr[c]->GetTrackPt(i)->GetEntries()) fCorrBgr[c]->GetTrackPt(i)->Write();
1513 for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) if(fCorrBgr[c]->GetZ(i)->GetEntries()) fCorrBgr[c]->GetZ(i)->Write();
1514 for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) if(fCorrBgr[c]->GetXi(i)->GetEntries()) fCorrBgr[c]->GetXi(i)->Write();
1517 if(fh1FFTrackPtBackFolded[i] && fh1FFTrackPtBackFolded[i]->GetEntries()) fh1FFTrackPtBackFolded[i]->Write();
1518 if(fh1FFZBackFolded[i] && fh1FFZBackFolded[i]->GetEntries()) fh1FFZBackFolded[i]->Write();
1519 if(fh1FFXiBackFolded[i] && fh1FFXiBackFolded[i]->GetEntries()) fh1FFXiBackFolded[i]->Write();
1522 if(fh1FFRatioTrackPtFolded[i] && fh1FFRatioTrackPtFolded[i]->GetEntries()) fh1FFRatioTrackPtFolded[i]->Write();
1523 if(fh1FFRatioZFolded[i] && fh1FFRatioZFolded[i]->GetEntries()) fh1FFRatioZFolded[i]->Write();
1524 if(fh1FFRatioXiFolded[i] && fh1FFRatioXiFolded[i]->GetEntries()) fh1FFRatioXiFolded[i]->Write();
1526 if(fh1FFRatioTrackPtBackFolded[i] && fh1FFRatioTrackPtBackFolded[i]->GetEntries()) fh1FFRatioTrackPtBackFolded[i]->Write();
1527 if(fh1FFRatioZBackFolded[i] && fh1FFRatioZBackFolded[i]->GetEntries()) fh1FFRatioZBackFolded[i]->Write();
1528 if(fh1FFRatioXiBackFolded[i] &&fh1FFRatioXiBackFolded[i]->GetEntries()) fh1FFRatioXiBackFolded[i]->Write();
1530 if(fh1FoldingCorrPt[i] && fh1FoldingCorrPt[i]->GetEntries()) fh1FoldingCorrPt[i]->Write(); // TEST!!!
1531 if(fh1FoldingCorrZ[i] && fh1FoldingCorrZ[i]->GetEntries()) fh1FoldingCorrZ[i]->Write(); // TEST!!!
1532 if(fh1FoldingCorrXi[i] && fh1FoldingCorrXi[i]->GetEntries()) fh1FoldingCorrXi[i]->Write(); // TEST!!!
1535 // inclusive track pt
1537 for(Int_t c=0; c<fNCorrectionLevelsSinglePt; c++) if(fCorrSinglePt[c]->GetTrackPt(0)->GetEntries()) fCorrSinglePt[c]->GetTrackPt(0)->Write();
1538 if(fh1SingleTrackPtBackFolded) fh1SingleTrackPtBackFolded->Write();
1539 if(fh1RatioSingleTrackPtFolded) fh1RatioSingleTrackPtFolded->Write();
1540 if(fh1RatioSingleTrackPtBackFolded) fh1RatioSingleTrackPtBackFolded->Write();
1548 //____________________________________________________________________________________________________________________________________
1549 THnSparse* AliFragmentationFunctionCorrections::TH1toSparse(const TH1F* hist, TString strName, TString strTit, const Bool_t fillConst)
1551 // copy 1-dimensional histo to THnSparse
1552 // if fillConst TRUE, create THnSparse with same binning as hist but all entries = 1
1553 // histos with variable bin size are supported
1555 // note: function returns pointer to 'new' THnSparse on heap, object needs to be deleted by user
1557 THnSparseF* fhnHist;
1559 Int_t nBins = hist->GetXaxis()->GetNbins();
1560 Int_t nBinsVec[1] = { nBins };
1562 const Double_t* binsVec = hist->GetXaxis()->GetXbins()->GetArray();
1564 if(binsVec){ // binsVec only neq NULL if histo was rebinned before
1566 fhnHist = new THnSparseF(strName,strTit,1,nBinsVec/*,binMinVec,binMaxVec*/);
1567 fhnHist->SetBinEdges(0,binsVec);
1569 else{ // uniform bin size
1571 Double_t xMin = hist->GetXaxis()->GetXmin();
1572 Double_t xMax = hist->GetXaxis()->GetXmax();
1574 Double_t binMinVec[1] = { xMin };
1575 Double_t binMaxVec[1] = { xMax };
1577 fhnHist = new THnSparseF(strName,strTit,1,nBinsVec,binMinVec,binMaxVec);
1581 for(Int_t bin=1; bin<=nBins; bin++){
1583 Double_t binCenter = fhnHist->GetAxis(0)->GetBinCenter(bin);
1585 Double_t binCoord[] = {binCenter};
1586 fhnHist->Fill(binCoord,1); // initially need to create the bin
1588 Long64_t binIndex = fhnHist->GetBin(binCoord);
1590 Double_t cont = hist->GetBinContent(bin);
1591 Double_t err = hist->GetBinError(bin);
1598 fhnHist->SetBinContent(binIndex,cont);
1599 fhnHist->SetBinError(binIndex,err);
1605 //______________________________________________________________________________________________________________________________________________
1606 TH1F* AliFragmentationFunctionCorrections::Unfold(THnSparse* hnHist, const THnSparse* hnResponse, const THnSparse* hnEff, const Int_t nIter,
1607 const Bool_t useCorrelatedErrors, const THnSparse* hnPrior)
1609 // unfold input histo
1611 AliCFUnfolding unfolding("unfolding","",1,hnResponse,hnEff,hnHist,hnPrior); // arg3: nVar; hnEff required, hnPrior not (defaults to 0x0)
1612 unfolding.SetMaxNumberOfIterations(nIter);
1613 // unfolding.SetMaxChi2PerDOF(1.e-07); // OBSOLETE !!!
1614 // if(useSmoothing) unfolding.UseSmoothing();
1616 if(useCorrelatedErrors) unfolding.SetUseCorrelatedErrors();
1619 THnSparse* unfolded = unfolding.GetUnfolded();
1621 TString hnameUnf = hnHist->GetName();
1622 hnameUnf.Append("_unf");
1623 unfolded->SetNameTitle(hnameUnf,hnHist->GetTitle());
1625 TH1F* h1Unfolded = (TH1F*) unfolded->Projection(0);
1626 h1Unfolded->SetNameTitle(hnHist->GetName(),hnHist->GetTitle());
1631 //___________________________________________________________________________________________________________________________
1632 void AliFragmentationFunctionCorrections::UnfoldHistos(const Int_t nIter, const Bool_t useCorrelatedErrors, const Int_t type)
1634 // loop over jet pt slices and unfold dN/dpt spectra
1636 TString labelCorr = fCorrFF[fNCorrectionLevels-1]->GetLabel();
1637 if(!labelCorr.Contains("unfold")) AddCorrectionLevel("unfold");
1639 for(Int_t i=0; i<fNJetPtSlices; i++){
1642 if(type == kFlagPt) hist = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i); // level -2: before unfolding, level -1: unfolded
1643 else if(type == kFlagZ) hist = fCorrFF[fNCorrectionLevels-2]->GetZ(i); // level -2: before unfolding, level -1: unfolded
1644 else if(type == kFlagXi) hist = fCorrFF[fNCorrectionLevels-2]->GetXi(i); // level -2: before unfolding, level -1: unfolded
1646 Printf("%s%d unknown type",(char*)__FILE__,__LINE__);
1651 Printf("%s%d no histo found ",(char*)__FILE__,__LINE__);
1656 THnSparse* hnResponse = 0;
1657 if(type == kFlagPt) hnResponse = fhnResponsePt[i];
1658 else if(type == kFlagZ) hnResponse = fhnResponseZ[i];
1659 else if(type == kFlagXi) hnResponse = fhnResponseXi[i];
1662 if(type == kFlagPt && fh1FFTrackPtPrior[i] && ((TString(fh1FFTrackPtPrior[i]->GetName())).Length() > 0) ) hPrior = fh1FFTrackPtPrior[i];
1663 else if(type == kFlagZ && fh1FFZPrior[i] && ((TString(fh1FFZPrior[i]->GetName())).Length() > 0) ) hPrior = fh1FFZPrior[i];
1664 else if(type == kFlagXi && fh1FFXiPrior[i] && ((TString(fh1FFXiPrior[i]->GetName())).Length() > 0) ) hPrior = fh1FFXiPrior[i];
1667 TString histNameTHn;
1668 histNameTHn = hist->GetName();
1669 histNameTHn.ReplaceAll("TH1","THn");
1672 TString priorNameTHn;
1674 priorNameTHn = hPrior->GetName();
1675 priorNameTHn.ReplaceAll("TH1","THn");
1678 TString histNameBackFolded;
1679 histNameBackFolded = hist->GetName();
1680 histNameBackFolded.Append("_backfold");
1682 TString histNameRatioFolded;
1683 histNameRatioFolded = hist->GetName();
1684 histNameRatioFolded.ReplaceAll("fh1FF","hRatioFF");
1686 histNameRatioFolded.Append("_unfold");
1688 TString histNameRatioBackFolded;
1689 histNameRatioBackFolded = hist->GetName();
1690 histNameRatioBackFolded.ReplaceAll("fh1FF","hRatioFF");
1691 histNameRatioBackFolded.Append("_backfold");
1693 THnSparse* hnHist = TH1toSparse(hist,histNameTHn,hist->GetTitle());
1694 THnSparse* hnFlatEfficiency = TH1toSparse(hist,"fhnEfficiency","eff",kTRUE); // could optionally also use real eff
1695 THnSparse* hnPrior = 0;
1696 if(hPrior) hnPrior = TH1toSparse(hPrior,priorNameTHn,hPrior->GetTitle());
1698 //THnSparse* hnUnfolded
1699 // = Unfold(hnHist,hnResponse,hnFlatEfficiency,nIter,useCorrelatedErrors,hnPrior);
1701 //TH1F* hUnfolded = (TH1F*) hnUnfolded->Projection(0);
1702 //hUnfolded->SetNameTitle(hist->GetName(),hist->GetTitle());
1705 = Unfold(hnHist,hnResponse,hnFlatEfficiency,nIter,useCorrelatedErrors,hnPrior);
1707 //TH1F* hUnfolded = (TH1F*) hnUnfolded->Projection(0);
1708 //hUnfolded->SetNameTitle(hist->GetName(),hist->GetTitle());
1711 for(Int_t bin=1; bin<=hist->GetNbinsX(); bin++){
1713 Double_t contOrg = hist->GetBinContent(bin);
1714 Double_t errOrg = hist->GetBinError(bin);
1716 Double_t contUnf = hUnfolded->GetBinContent(bin);
1717 //Double_t errUnf = hUnfolded->GetBinError(bin);
1719 Double_t relErrOrg = 0;
1720 if(contOrg) relErrOrg = errOrg/contOrg;
1722 // std::cout<<" cont original "<<contOrg<<" err original "<<errOrg<<std::endl;
1723 // std::cout<<" cont unf "<<contUnf<<" errUnf"<<errUnf<<std::endl;
1724 // std::cout<<" err/cont original "<<relErrOrg<<std::endl;
1725 // if(contUnf) std::cout<<" err/conf unf "<<errUnf/contUnf<<std::endl;
1727 hUnfolded->SetBinError(bin,relErrOrg*contUnf);
1729 // if(hUnfolded->GetBinContent(bin)){
1730 // std::cout<<" modified err unfolded "<<hUnfolded->GetBinError(bin)/hUnfolded->GetBinContent(bin)<<std::endl;
1735 if(type == kFlagPt) fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hUnfolded,0,0);
1736 if(type == kFlagZ) fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,0,hUnfolded,0);
1737 if(type == kFlagXi) fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,0,0,hUnfolded);
1739 // backfolding: apply response matrix to unfolded spectrum
1740 TH1F* hBackFolded = ApplyResponse(hUnfolded,hnResponse);
1741 hBackFolded->SetNameTitle(histNameBackFolded,hist->GetTitle());
1743 if(type == kFlagPt) fh1FFTrackPtBackFolded[i] = hBackFolded;
1744 if(type == kFlagZ) fh1FFZBackFolded[i] = hBackFolded;
1745 if(type == kFlagXi) fh1FFXiBackFolded[i] = hBackFolded;
1747 // ratio unfolded to original histo
1748 TH1F* hRatioUnfolded = (TH1F*) hUnfolded->Clone(histNameRatioFolded);
1749 hRatioUnfolded->Reset();
1750 hRatioUnfolded->Divide(hUnfolded,hist,1,1,"B");
1752 if(type == kFlagPt) fh1FFRatioTrackPtFolded[i] = hRatioUnfolded;
1753 if(type == kFlagZ) fh1FFRatioZFolded[i] = hRatioUnfolded;
1754 if(type == kFlagXi) fh1FFRatioXiFolded[i] = hRatioUnfolded;
1757 // ratio backfolded to original histo
1758 TH1F* hRatioBackFolded = (TH1F*) hBackFolded->Clone(histNameRatioBackFolded);
1759 hRatioBackFolded->Reset();
1760 hRatioBackFolded->Divide(hBackFolded,hist,1,1,"B");
1762 if(type == kFlagPt) fh1FFRatioTrackPtBackFolded[i] = hRatioBackFolded;
1763 if(type == kFlagZ) fh1FFRatioZBackFolded[i] = hRatioBackFolded;
1764 if(type == kFlagXi) fh1FFRatioXiBackFolded[i] = hRatioBackFolded;
1767 delete hnFlatEfficiency;
1771 //_____________________________________________________________________________________________________
1772 void AliFragmentationFunctionCorrections::UnfoldPt(const Int_t nIter, const Bool_t useCorrelatedErrors)
1775 Int_t type = kFlagPt;
1776 UnfoldHistos(nIter, useCorrelatedErrors, type);
1779 //_____________________________________________________________________________________________________
1780 void AliFragmentationFunctionCorrections::UnfoldZ(const Int_t nIter, const Bool_t useCorrelatedErrors)
1783 Int_t type = kFlagZ;
1784 UnfoldHistos(nIter, useCorrelatedErrors, type);
1787 //_____________________________________________________________________________________________________
1788 void AliFragmentationFunctionCorrections::UnfoldXi(const Int_t nIter, const Bool_t useCorrelatedErrors)
1791 Int_t type = kFlagXi;
1792 UnfoldHistos(nIter, useCorrelatedErrors, type);
1796 //______________________________________________________________________________________________
1797 TH1F* AliFragmentationFunctionCorrections::ApplyResponse(const TH1F* hist, THnSparse* hnResponse)
1799 // apply (multiply) response matrix to hist
1801 TH1F* hOut = new TH1F(*hist);
1804 const Int_t axisM = 0;
1805 const Int_t axisT = 1;
1807 Int_t nBinsM = hnResponse->GetAxis(axisM)->GetNbins();
1808 Int_t nBinsT = hnResponse->GetAxis(axisT)->GetNbins();
1810 // response matrix normalization
1811 // do it in this function and not when reading response, since for 'proper' normalization errors are difficult to assign
1812 // so for unfolding proper we leave it to CORRFW ...
1814 Double_t normFacResponse[nBinsT+1];
1816 for(Int_t iT=1; iT<=nBinsT; iT++){
1818 Double_t sumResp = 0;
1820 for(Int_t iM=1; iM<=nBinsM; iM++){
1822 Double_t coordM = hnResponse->GetAxis(axisM)->GetBinCenter(iM);
1823 Double_t coordT = hnResponse->GetAxis(axisT)->GetBinCenter(iT);
1825 Double_t binCoord[] = {coordM,coordT};
1827 Long64_t binIndex = hnResponse->GetBin(binCoord);
1829 Double_t resp = hnResponse->GetBinContent(binIndex);
1834 normFacResponse[iT] = 0;
1835 if(sumResp) normFacResponse[iT] = 1/sumResp;
1840 for(Int_t iM=1; iM<=nBinsM; iM++){
1844 for(Int_t iT=1; iT<=nBinsT; iT++){
1846 Double_t contT = hist->GetBinContent(iT);
1848 Double_t coordM = hnResponse->GetAxis(axisM)->GetBinCenter(iM);
1849 Double_t coordT = hnResponse->GetAxis(axisT)->GetBinCenter(iT);
1851 Double_t binCoord[] = {coordM,coordT};
1853 Long64_t binIndex = hnResponse->GetBin(binCoord);
1855 Double_t resp = hnResponse->GetBinContent(binIndex);
1857 contM += resp*normFacResponse[iT]*contT;
1860 hOut->SetBinContent(iM,contM);
1866 //_______________________________________________________________________________________________________
1867 void AliFragmentationFunctionCorrections::ReadEfficiency(TString strfile, TString strdir, TString strlist)
1869 // read reconstruction efficiency from file
1870 // argument strlist optional - read from directory strdir if not specified
1872 // temporary histos to hold histos from file
1873 TH1F* hEffPt[fNJetPtSlices];
1874 TH1F* hEffZ[fNJetPtSlices];
1875 TH1F* hEffXi[fNJetPtSlices];
1877 for(Int_t i=0; i<fNJetPtSlices; i++) hEffPt[i] = 0;
1878 for(Int_t i=0; i<fNJetPtSlices; i++) hEffZ[i] = 0;
1879 for(Int_t i=0; i<fNJetPtSlices; i++) hEffXi[i] = 0;
1881 TFile f(strfile,"READ");
1884 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1888 if(fDebug>0) Printf("%s:%d -- read efficiencies from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
1890 if(strdir && strdir.Length()) gDirectory->cd(strdir);
1894 if(strlist && strlist.Length()){
1896 if(!(list = (TList*) gDirectory->Get(strlist))){
1897 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1902 for(Int_t i=0; i<fNJetPtSlices; i++){
1904 Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
1905 Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
1907 TString strNameEffPt(Form("hEffPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
1908 TString strNameEffZ(Form("hEffZ_%02d_%02d",jetPtLoLim,jetPtUpLim));
1909 TString strNameEffXi(Form("hEffXi_%02d_%02d",jetPtLoLim,jetPtUpLim));
1913 hEffPt[i] = (TH1F*) list->FindObject(strNameEffPt);
1914 hEffZ[i] = (TH1F*) list->FindObject(strNameEffZ);
1915 hEffXi[i] = (TH1F*) list->FindObject(strNameEffXi);
1918 hEffPt[i] = (TH1F*) gDirectory->Get(strNameEffPt);
1919 hEffZ[i] = (TH1F*) gDirectory->Get(strNameEffZ);
1920 hEffXi[i] = (TH1F*) gDirectory->Get(strNameEffXi);
1924 Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffPt.Data());
1928 Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffZ.Data());
1932 Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffXi.Data());
1936 if(fNHistoBinsPt[i]) hEffPt[i] = (TH1F*) hEffPt[i]->Rebin(fNHistoBinsPt[i],strNameEffPt+"_rebin",fHistoBinsPt[i]->GetArray());
1937 if(fNHistoBinsZ[i]) hEffZ[i] = (TH1F*) hEffZ[i]->Rebin(fNHistoBinsZ[i],strNameEffZ+"_rebin",fHistoBinsZ[i]->GetArray());
1938 if(fNHistoBinsXi[i]) hEffXi[i] = (TH1F*) hEffXi[i]->Rebin(fNHistoBinsXi[i],strNameEffXi+"_rebin",fHistoBinsXi[i]->GetArray());
1940 if(hEffPt[i]) hEffPt[i]->SetDirectory(0);
1941 if(hEffZ[i]) hEffZ[i]->SetDirectory(0);
1942 if(hEffXi[i]) hEffXi[i]->SetDirectory(0);
1944 } // jet slices loop
1948 for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
1949 if(hEffPt[i]) new(fh1EffPt[i]) TH1F(*hEffPt[i]);
1950 if(hEffZ[i]) new(fh1EffZ[i]) TH1F(*hEffZ[i]);
1951 if(hEffXi[i]) new(fh1EffXi[i]) TH1F(*hEffXi[i]);
1955 //___________________________________________________________________________________________________________
1956 void AliFragmentationFunctionCorrections::ReadBgrEfficiency(TString strfile, TString strdir, TString strlist)
1958 // read bgr eff from file
1959 // argument strlist optional - read from directory strdir if not specified
1961 TH1F* hEffPtBgr[fNJetPtSlices];
1962 TH1F* hEffZBgr [fNJetPtSlices];
1963 TH1F* hEffXiBgr[fNJetPtSlices];
1965 for(Int_t i=0; i<fNJetPtSlices; i++) hEffPtBgr[i] = 0;
1966 for(Int_t i=0; i<fNJetPtSlices; i++) hEffZBgr[i] = 0;
1967 for(Int_t i=0; i<fNJetPtSlices; i++) hEffXiBgr[i] = 0;
1970 TFile f(strfile,"READ");
1973 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1977 if(fDebug>0) Printf("%s:%d -- read bgr efficiencies from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
1979 if(strdir && strdir.Length()) gDirectory->cd(strdir);
1983 if(strlist && strlist.Length()){
1985 if(!(list = (TList*) gDirectory->Get(strlist))){
1986 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1991 for(Int_t i=0; i<fNJetPtSlices; i++){
1993 Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
1994 Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
1996 TString strNameEffPtBgr(Form("hEffPtBgr%02dto%02d",jetPtLoLim,jetPtUpLim));
1997 TString strNameEffZBgr(Form("hEffZBgr%02dto%02d",jetPtLoLim,jetPtUpLim));
1998 TString strNameEffXiBgr(Form("hEffXiBgr%02dto%02d",jetPtLoLim,jetPtUpLim));
2002 hEffPtBgr[i] = (TH1F*) list->FindObject(strNameEffPtBgr);
2003 hEffZBgr[i] = (TH1F*) list->FindObject(strNameEffZBgr);
2004 hEffXiBgr[i] = (TH1F*) list->FindObject(strNameEffXiBgr);
2007 hEffPtBgr[i] = (TH1F*) gDirectory->Get(strNameEffPtBgr);
2008 hEffZBgr[i] = (TH1F*) gDirectory->Get(strNameEffZBgr);
2009 hEffXiBgr[i] = (TH1F*) gDirectory->Get(strNameEffXiBgr);
2013 Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffPtBgr.Data());
2017 Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffZBgr.Data());
2021 Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffXiBgr.Data());
2025 if(fNHistoBinsPt[i]) hEffPtBgr[i] = (TH1F*) hEffPtBgr[i]->Rebin(fNHistoBinsPt[i],strNameEffPtBgr+"_rebin",fHistoBinsPt[i]->GetArray());
2026 if(fNHistoBinsZ[i]) hEffZBgr[i] = (TH1F*) hEffZBgr[i]->Rebin(fNHistoBinsZ[i],strNameEffZBgr+"_rebin",fHistoBinsZ[i]->GetArray());
2027 if(fNHistoBinsXi[i]) hEffXiBgr[i] = (TH1F*) hEffXiBgr[i]->Rebin(fNHistoBinsXi[i],strNameEffXiBgr+"_rebin",fHistoBinsXi[i]->GetArray());
2029 if(hEffPtBgr[i]) hEffPtBgr[i]->SetDirectory(0);
2030 if(hEffZBgr[i]) hEffZBgr[i]->SetDirectory(0);
2031 if(hEffXiBgr[i]) hEffXiBgr[i]->SetDirectory(0);
2033 } // jet slices loop
2037 for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
2038 if(hEffPtBgr[i]) new(fh1EffBgrPt[i]) TH1F(*hEffPtBgr[i]);
2039 if(hEffZBgr[i]) new(fh1EffBgrZ[i]) TH1F(*hEffZBgr[i]);
2040 if(hEffXiBgr[i]) new(fh1EffBgrXi[i]) TH1F(*hEffXiBgr[i]);
2044 // ________________________________________________
2045 void AliFragmentationFunctionCorrections::EffCorr()
2047 // apply efficiency correction
2049 AddCorrectionLevel("eff");
2051 for(Int_t i=0; i<fNJetPtSlices; i++){
2053 TH1F* histPt = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i);
2054 TH1F* histZ = fCorrFF[fNCorrectionLevels-2]->GetZ(i);
2055 TH1F* histXi = fCorrFF[fNCorrectionLevels-2]->GetXi(i);
2057 TString histNamePt = histPt->GetName();
2058 TString histNameZ = histZ->GetName();
2059 TString histNameXi = histXi->GetName();
2062 TH1F* hFFTrackPtEffCorr = (TH1F*) histPt->Clone(histNamePt);
2063 hFFTrackPtEffCorr->Divide(histPt,fh1EffPt[i],1,1,"");
2065 TH1F* hFFZEffCorr = (TH1F*) histZ->Clone(histNameZ);
2066 hFFZEffCorr->Divide(histZ,fh1EffZ[i],1,1,"");
2068 TH1F* hFFXiEffCorr = (TH1F*) histXi->Clone(histNameXi);
2069 hFFXiEffCorr->Divide(histXi,fh1EffXi[i],1,1,"");
2071 fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hFFTrackPtEffCorr,hFFZEffCorr,hFFXiEffCorr);
2075 //___________________________________________________
2076 void AliFragmentationFunctionCorrections::EffCorrBgr()
2078 // apply efficiency correction to bgr distributions
2080 AddCorrectionLevelBgr("eff");
2082 Printf("%s:%d -- apply efficiency correction, corrLevel %d",(char*)__FILE__,__LINE__,fNCorrectionLevels-1);
2084 for(Int_t i=0; i<fNJetPtSlices; i++){
2086 TH1F* histPt = fCorrBgr[fNCorrectionLevelsBgr-2]->GetTrackPt(i);
2087 TH1F* histZ = fCorrBgr[fNCorrectionLevelsBgr-2]->GetZ(i);
2088 TH1F* histXi = fCorrBgr[fNCorrectionLevelsBgr-2]->GetXi(i);
2090 TString histNamePt = histPt->GetName();
2091 TString histNameZ = histZ->GetName();
2092 TString histNameXi = histXi->GetName();
2095 TH1F* hFFTrackPtEffCorr = (TH1F*) histPt->Clone(histNamePt);
2096 hFFTrackPtEffCorr->Divide(histPt,fh1EffPt[i],1,1,"");
2098 TH1F* hFFZEffCorr = (TH1F*) histZ->Clone(histNameZ);
2099 hFFZEffCorr->Divide(histZ,fh1EffZ[i],1,1,"");
2101 TH1F* hFFXiEffCorr = (TH1F*) histXi->Clone(histNameXi);
2102 hFFXiEffCorr->Divide(histXi,fh1EffXi[i],1,1,"");
2104 fCorrBgr[fNCorrectionLevelsBgr-1]->AddCorrHistos(i,hFFTrackPtEffCorr,hFFZEffCorr,hFFXiEffCorr);
2108 //_____________________________________________________
2109 void AliFragmentationFunctionCorrections::SubtractBgr(Double_t sysErr)
2111 // subtract bgr distribution from FF
2112 // the current corr level is used for both
2114 AddCorrectionLevel("bgrSub");
2116 fstream ascii_out_pt;
2117 fstream ascii_out_z;
2118 fstream ascii_out_xi;
2120 if(sysErr) ascii_out_pt.open("sysErrUE_pt.txt",std::ios::out);
2121 if(sysErr) ascii_out_z.open("sysErrUE_z.txt",std::ios::out);
2122 if(sysErr) ascii_out_xi.open("sysErrUE_xi.txt",std::ios::out);
2124 for(Int_t i=0; i<fNJetPtSlices; i++){ // jet slices
2126 TH1F* histPt = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i);
2127 TH1F* histZ = fCorrFF[fNCorrectionLevels-2]->GetZ(i);
2128 TH1F* histXi = fCorrFF[fNCorrectionLevels-2]->GetXi(i);
2130 TH1F* histPtBgr = fCorrBgr[fNCorrectionLevelsBgr-1]->GetTrackPt(i);
2131 TH1F* histZBgr = fCorrBgr[fNCorrectionLevelsBgr-1]->GetZ(i);
2132 TH1F* histXiBgr = fCorrBgr[fNCorrectionLevelsBgr-1]->GetXi(i);
2134 TString histNamePt = histPt->GetName();
2135 TString histNameZ = histZ->GetName();
2136 TString histNameXi = histXi->GetName();
2139 TH1F* hFFTrackPtBgrSub = (TH1F*) histPt->Clone(histNamePt);
2140 hFFTrackPtBgrSub->Add(histPtBgr,-1);
2142 TH1F* hFFZBgrSub = (TH1F*) histZ->Clone(histNameZ);
2143 hFFZBgrSub->Add(histZBgr,-1);
2145 TH1F* hFFXiBgrSub = (TH1F*) histXi->Clone(histNameXi);
2146 hFFXiBgrSub->Add(histXiBgr,-1);
2148 fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hFFTrackPtBgrSub,hFFZBgrSub,hFFXiBgrSub);
2152 for(Int_t bin=1; bin<=histPt->GetNbinsX(); bin++){
2154 //Double_t sigPlBgr = histPt->GetBinContent(bin);
2155 Double_t sig = hFFTrackPtBgrSub->GetBinContent(bin);
2156 Double_t bgr = histPtBgr->GetBinContent(bin);
2158 Double_t relErr = 0;
2159 if(sig>0) relErr = sysErr*bgr/sig;
2161 // std::cout<<" pt bin "<<bin<<" mean "<<histPt->GetBinCenter(bin)
2162 // <<" sigPlBgr "<<sigPlBgr<<" sig "<<sig<<" bgr "<<bgr<<" relErr "<<relErr<<std::endl;
2164 ascii_out_pt<<i<<" "<<bin<<" "<<relErr<<std::endl;
2168 for(Int_t bin=1; bin<=histZ->GetNbinsX(); bin++){
2170 //Double_t sigPlBgr = histZ->GetBinContent(bin);
2171 Double_t sig = hFFZBgrSub->GetBinContent(bin);
2172 Double_t bgr = histZBgr->GetBinContent(bin);
2174 Double_t relErr = 0;
2175 if(sig>0) relErr = sysErr*bgr/sig;
2177 std::cout<<" z bin "<<bin<<" mean "<<histZ->GetBinCenter(bin)
2178 <<" sigPlBgr "<<histZ->GetBinContent(bin)<<" sig "<<sig<<" bgr "<<bgr<<" relErr "<<relErr<<std::endl;
2180 ascii_out_z<<i<<" "<<bin<<" "<<relErr<<std::endl;
2184 for(Int_t bin=1; bin<=histXi->GetNbinsX(); bin++){
2186 //Double_t sigPlBgr = histXi->GetBinContent(bin);
2187 Double_t sig = hFFXiBgrSub->GetBinContent(bin);
2188 Double_t bgr = histXiBgr->GetBinContent(bin);
2190 Double_t relErr = 0;
2191 if(sig>0) relErr = sysErr*bgr/sig;
2193 std::cout<<" xi bin "<<bin<<" mean "<<histXi->GetBinCenter(bin)
2194 <<" sigPlBgr "<<histXi->GetBinContent(bin)<<" sig "<<sig<<" bgr "<<bgr<<" relErr "<<relErr<<std::endl;
2196 ascii_out_xi<<i<<" "<<bin<<" "<<relErr<<std::endl;
2201 if(sysErr) ascii_out_pt.close();
2202 if(sysErr) ascii_out_xi.close();
2207 //________________________________________________________________________________________________________________
2208 void AliFragmentationFunctionCorrections::WriteSingleTrackEff(TString strInfile, TString strID, TString strOutfile,
2209 Bool_t updateOutfile, TString strOutDir,TString strPostfix)
2211 // read task ouput from MC and write single track eff - standard dir/list
2213 TString strdir = "PWGJE_FragmentationFunction_" + strID;
2214 TString strlist = "fracfunc_" + strID;
2216 WriteSingleTrackEff(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir,strPostfix);
2219 //___________________________________________________________________________________________________________________________________
2220 void AliFragmentationFunctionCorrections::WriteSingleTrackEff(TString strInfile, TString strdir, TString strlist,
2221 TString strOutfile, Bool_t updateOutfile, TString strOutDir,TString strPostfix)
2223 // read task output from MC and write single track eff as function of pt, eta, phi
2224 // argument strlist optional - read from directory strdir if not specified
2225 // write eff to file stroutfile - by default only 'update' file (don't overwrite)
2228 TH1D* hdNdptTracksMCPrimGen;
2229 TH2D* hdNdetadphiTracksMCPrimGen;
2231 TH1D* hdNdptTracksMCPrimRec;
2232 TH2D* hdNdetadphiTracksMCPrimRec;
2235 TFile f(strInfile,"READ");
2238 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2242 if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: open task ouput file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
2244 if(strdir && strdir.Length()){
2245 if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: change dir to %s",(char*)__FILE__,__LINE__,strdir.Data());
2246 gDirectory->cd(strdir);
2251 if(strlist && strlist.Length()){
2253 if(!(list = (TList*) gDirectory->Get(strlist))){
2254 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2260 TString hnamePtRecEffGen = "fh1TrackQAPtRecEffGen";
2261 if(strPostfix.Length()) hnamePtRecEffGen.Form("fh1TrackQAPtRecEffGen_%s",strPostfix.Data());
2263 TString hnamePtRecEffRec = "fh1TrackQAPtRecEffRec";
2264 if(strPostfix.Length()) hnamePtRecEffRec.Form("fh1TrackQAPtRecEffRec_%s",strPostfix.Data());
2266 TString hnameEtaPhiRecEffGen = "fh2TrackQAEtaPhiRecEffGen";
2267 if(strPostfix.Length()) hnameEtaPhiRecEffGen.Form("fh2TrackQAEtaPhiRecEffGen_%s",strPostfix.Data());
2269 TString hnameEtaPhiRecEffRec = "fh2TrackQAEtaPhiRecEffRec";
2270 if(strPostfix.Length()) hnameEtaPhiRecEffRec.Form("fh2TrackQAEtaPhiRecEffRec_%s",strPostfix.Data());
2274 hdNdptTracksMCPrimGen = (TH1D*) list->FindObject(hnamePtRecEffGen);
2275 hdNdetadphiTracksMCPrimGen = (TH2D*) list->FindObject(hnameEtaPhiRecEffGen);
2277 hdNdptTracksMCPrimRec = (TH1D*) list->FindObject(hnamePtRecEffRec);
2278 hdNdetadphiTracksMCPrimRec = (TH2D*) list->FindObject(hnameEtaPhiRecEffRec);
2281 hdNdptTracksMCPrimGen = (TH1D*) gDirectory->Get(hnamePtRecEffGen);
2282 hdNdetadphiTracksMCPrimGen = (TH2D*) gDirectory->Get(hnameEtaPhiRecEffGen);
2284 hdNdptTracksMCPrimRec = (TH1D*) gDirectory->Get(hnamePtRecEffRec);
2285 hdNdetadphiTracksMCPrimRec = (TH2D*) gDirectory->Get(hnameEtaPhiRecEffRec);
2288 hdNdptTracksMCPrimGen->SetDirectory(0);
2289 hdNdetadphiTracksMCPrimGen->SetDirectory(0);
2290 hdNdptTracksMCPrimRec->SetDirectory(0);
2291 hdNdetadphiTracksMCPrimRec->SetDirectory(0);
2295 // projections: dN/deta, dN/dphi
2297 TH1D* hdNdetaTracksMCPrimGen = (TH1D*) hdNdetadphiTracksMCPrimGen->ProjectionX("hdNdetaTracksMcPrimGen");
2298 TH1D* hdNdphiTracksMCPrimGen = (TH1D*) hdNdetadphiTracksMCPrimGen->ProjectionY("hdNdphiTracksMcPrimGen");
2300 TH1D* hdNdetaTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionX("hdNdetaTracksMcPrimRec");
2301 TH1D* hdNdphiTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionY("hdNdphiTracksMcPrimRec");
2305 TString strNamePtGen = "hTrackPtGenPrim";
2306 TString strNamePtRec = "hTrackPtRecPrim";
2308 if(fNHistoBinsSinglePt) hdNdptTracksMCPrimGen = (TH1D*) hdNdptTracksMCPrimGen->Rebin(fNHistoBinsSinglePt,strNamePtGen,fHistoBinsSinglePt->GetArray());
2309 if(fNHistoBinsSinglePt) hdNdptTracksMCPrimRec = (TH1D*) hdNdptTracksMCPrimRec->Rebin(fNHistoBinsSinglePt,strNamePtRec,fHistoBinsSinglePt->GetArray());
2311 // efficiency: divide
2313 TString hNameTrackEffPt = "hSingleTrackEffPt";
2314 if(strPostfix.Length()) hNameTrackEffPt.Form("hSingleTrackEffPt_%s",strPostfix.Data());
2316 TString hNameTrackEffEta = "hSingleTrackEffEta";
2317 if(strPostfix.Length()) hNameTrackEffEta.Form("hSingleTrackEffEta_%s",strPostfix.Data());
2319 TString hNameTrackEffPhi = "hSingleTrackEffPhi";
2320 if(strPostfix.Length()) hNameTrackEffPhi.Form("hSingleTrackEffPhi_%s",strPostfix.Data());
2323 TH1F* hSingleTrackEffPt = (TH1F*) hdNdptTracksMCPrimRec->Clone(hNameTrackEffPt);
2324 hSingleTrackEffPt->Divide(hdNdptTracksMCPrimRec,hdNdptTracksMCPrimGen,1,1,"B"); // binominal errors
2326 TH1F* hSingleTrackEffEta = (TH1F*) hdNdetaTracksMCPrimRec->Clone(hNameTrackEffEta);
2327 hSingleTrackEffEta->Divide(hdNdetaTracksMCPrimRec,hdNdetaTracksMCPrimGen,1,1,"B"); // binominal errors
2329 TH1F* hSingleTrackEffPhi = (TH1F*) hdNdphiTracksMCPrimRec->Clone(hNameTrackEffPhi);
2330 hSingleTrackEffPhi->Divide(hdNdphiTracksMCPrimRec,hdNdphiTracksMCPrimGen,1,1,"B"); // binominal errors
2333 TString outfileOption = "RECREATE";
2334 if(updateOutfile) outfileOption = "UPDATE";
2336 TFile out(strOutfile,outfileOption);
2339 Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2343 if(fDebug>0) Printf("%s:%d -- write efficiency to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2345 if(strOutDir && strOutDir.Length()){
2348 if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd();
2350 dir = out.mkdir(strOutDir);
2355 hSingleTrackEffPt->Write();
2356 hSingleTrackEffEta->Write();
2357 hSingleTrackEffPhi->Write();
2361 delete hdNdptTracksMCPrimGen;
2362 delete hdNdetadphiTracksMCPrimGen;
2363 delete hdNdptTracksMCPrimRec;
2364 delete hdNdetadphiTracksMCPrimRec;
2367 //________________________________________________________________________________________________________________
2368 void AliFragmentationFunctionCorrections::WriteSingleTrackSecCorr(TString strInfile, TString strID, TString strOutfile,
2369 Bool_t updateOutfile, TString strOutDir)
2371 // read task ouput from MC and write single track eff - standard dir/list
2373 TString strdir = "PWGJE_FragmentationFunction_" + strID;
2374 TString strlist = "fracfunc_" + strID;
2376 WriteSingleTrackSecCorr(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir);
2379 //___________________________________________________________________________________________________________________________________
2380 void AliFragmentationFunctionCorrections::WriteSingleTrackSecCorr(TString strInfile, TString strdir, TString strlist,
2381 TString strOutfile, Bool_t updateOutfile, TString strOutDir)
2383 // read task output from MC and write single track secondaries contamination correction as function of pt, eta, phi
2384 // argument strlist optional - read from directory strdir if not specified
2385 // write corr factor to file stroutfile - by default only 'update' file (don't overwrite)
2387 TH1D* hdNdptTracksMCPrimRec;
2388 TH2D* hdNdetadphiTracksMCPrimRec;
2390 TH1D* hdNdptTracksMCSecRecNS;
2391 TH2D* hdNdetadphiTracksMCSecRecNS;
2393 TH1D* hdNdptTracksMCSecRecSsc;
2394 TH2D* hdNdetadphiTracksMCSecRecSsc;
2396 TFile f(strInfile,"READ");
2399 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2403 if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
2405 if(strdir && strdir.Length()) gDirectory->cd(strdir);
2409 if(strlist && strlist.Length()){
2411 if(!(list = (TList*) gDirectory->Get(strlist))){
2412 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2419 hdNdptTracksMCPrimRec = (TH1D*) list->FindObject("fh1TrackQAPtRecEffGen");
2420 hdNdetadphiTracksMCPrimRec = (TH2D*) list->FindObject("fh2TrackQAEtaPhiRecEffGen");
2422 hdNdptTracksMCSecRecNS = (TH1D*) list->FindObject("fh1TrackQAPtSecRecNS");
2423 hdNdetadphiTracksMCSecRecNS = (TH2D*) list->FindObject("fh2TrackQAEtaPhiSecRecNS");
2425 hdNdptTracksMCSecRecSsc = (TH1D*) list->FindObject("fh1TrackQAPtSecRecSsc");
2426 hdNdetadphiTracksMCSecRecSsc = (TH2D*) list->FindObject("fh2TrackQAEtaPhiSecRecSsc");
2430 hdNdptTracksMCPrimRec = (TH1D*) gDirectory->Get("fh1TrackQAPtRecEffGen");
2431 hdNdetadphiTracksMCPrimRec = (TH2D*) gDirectory->Get("fh2TrackQAEtaPhiRecEffGen");
2433 hdNdptTracksMCSecRecNS = (TH1D*) gDirectory->Get("fh1TrackQAPtSecRecNS");
2434 hdNdetadphiTracksMCSecRecNS = (TH2D*) gDirectory->Get("fh2TrackQAEtaPhiSecRecNS");
2436 hdNdptTracksMCSecRecSsc = (TH1D*) gDirectory->Get("fh1TrackQAPtSecRecSsc");
2437 hdNdetadphiTracksMCSecRecSsc = (TH2D*) gDirectory->Get("fh2TrackQAEtaPhiSecRecSsc");
2440 hdNdptTracksMCPrimRec->SetDirectory(0);
2441 hdNdetadphiTracksMCPrimRec->SetDirectory(0);
2443 hdNdptTracksMCSecRecNS->SetDirectory(0);
2444 hdNdetadphiTracksMCSecRecNS->SetDirectory(0);
2446 hdNdptTracksMCSecRecSsc->SetDirectory(0);
2447 hdNdetadphiTracksMCSecRecSsc->SetDirectory(0);
2451 // projections: dN/deta, dN/dphi
2453 TH1D* hdNdetaTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionX("hdNdetaTracksMcPrimRec");
2454 TH1D* hdNdphiTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionY("hdNdphiTracksMcPrimRec");
2456 TH1D* hdNdetaTracksMCSecRecNS = (TH1D*) hdNdetadphiTracksMCSecRecNS->ProjectionX("hdNdetaTracksMcSecRecNS");
2457 TH1D* hdNdphiTracksMCSecRecNS = (TH1D*) hdNdetadphiTracksMCSecRecNS->ProjectionY("hdNdphiTracksMcSecRecNS");
2459 TH1D* hdNdetaTracksMCSecRecSsc = (TH1D*) hdNdetadphiTracksMCSecRecSsc->ProjectionX("hdNdetaTracksMcSecRecSsc");
2460 TH1D* hdNdphiTracksMCSecRecSsc = (TH1D*) hdNdetadphiTracksMCSecRecSsc->ProjectionY("hdNdphiTracksMcSecRecSsc");
2464 TString strNamePtPrim = "hTrackPtPrim";
2465 TString strNamePtSec = "hTrackPtSec";
2467 if(fNHistoBinsSinglePt) hdNdptTracksMCPrimRec = (TH1D*) hdNdptTracksMCPrimRec->Rebin(fNHistoBinsSinglePt,strNamePtPrim,fHistoBinsSinglePt->GetArray());
2468 if(fNHistoBinsSinglePt) hdNdptTracksMCSecRecNS = (TH1D*) hdNdptTracksMCSecRecNS->Rebin(fNHistoBinsSinglePt,strNamePtSec,fHistoBinsSinglePt->GetArray());
2469 if(fNHistoBinsSinglePt) hdNdptTracksMCSecRecSsc = (TH1D*) hdNdptTracksMCSecRecSsc->Rebin(fNHistoBinsSinglePt,strNamePtSec,fHistoBinsSinglePt->GetArray());
2472 // secondary correction factor: divide prim/(prim+sec)
2474 TH1F* hSingleTrackSecCorrPt = (TH1F*) hdNdptTracksMCPrimRec->Clone("hSingleTrackSecCorrPt");
2475 TH1F* hSumPrimSecPt = (TH1F*) hdNdptTracksMCPrimRec->Clone("hSumPrimSecPt");
2476 hSumPrimSecPt->Add(hdNdptTracksMCSecRecNS);
2477 hSumPrimSecPt->Add(hdNdptTracksMCSecRecSsc);
2478 hSingleTrackSecCorrPt->Divide(hdNdptTracksMCPrimRec,hSumPrimSecPt,1,1,"B"); // binominal errors
2481 TH1F* hSingleTrackSecCorrEta = (TH1F*) hdNdetaTracksMCPrimRec->Clone("hSingleTrackSecCorrEta");
2482 TH1F* hSumPrimSecEta = (TH1F*) hdNdetaTracksMCPrimRec->Clone("hSumPrimSecEta");
2483 hSumPrimSecEta->Add(hdNdetaTracksMCSecRecNS);
2484 hSumPrimSecEta->Add(hdNdetaTracksMCSecRecSsc);
2485 hSingleTrackSecCorrEta->Divide(hdNdetaTracksMCPrimRec,hSumPrimSecEta,1,1,"B"); // binominal errors
2487 TH1F* hSingleTrackSecCorrPhi = (TH1F*) hdNdphiTracksMCPrimRec->Clone("hSingleTrackSecCorrPhi");
2488 TH1F* hSumPrimSecPhi = (TH1F*) hdNdphiTracksMCPrimRec->Clone("hSumPrimSecPhi");
2489 hSumPrimSecPhi->Add(hdNdphiTracksMCSecRecNS);
2490 hSumPrimSecPhi->Add(hdNdphiTracksMCSecRecSsc);
2491 hSingleTrackSecCorrPhi->Divide(hdNdphiTracksMCPrimRec,hSumPrimSecPhi,1,1,"B"); // binominal errors
2495 TString outfileOption = "RECREATE";
2496 if(updateOutfile) outfileOption = "UPDATE";
2498 TFile out(strOutfile,outfileOption);
2501 Printf("%s:%d -- error opening secCorr output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2505 if(fDebug>0) Printf("%s:%d -- write secCorr to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2507 if(strOutDir && strOutDir.Length()){
2510 if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd();
2512 dir = out.mkdir(strOutDir);
2517 hdNdptTracksMCSecRecNS->Write();
2518 hdNdetaTracksMCSecRecNS->Write();
2519 hdNdphiTracksMCSecRecNS->Write();
2521 hdNdptTracksMCSecRecSsc->Write();
2522 hdNdetaTracksMCSecRecSsc->Write();
2523 hdNdphiTracksMCSecRecSsc->Write();
2525 hSingleTrackSecCorrPt->Write();
2526 hSingleTrackSecCorrEta->Write();
2527 hSingleTrackSecCorrPhi->Write();
2531 delete hdNdptTracksMCPrimRec;
2532 delete hdNdetadphiTracksMCPrimRec;
2533 delete hdNdptTracksMCSecRecNS;
2534 delete hdNdetadphiTracksMCSecRecNS;
2535 delete hdNdptTracksMCSecRecSsc;
2536 delete hdNdetadphiTracksMCSecRecSsc;
2539 //________________________________________________________________________________________________________________
2540 void AliFragmentationFunctionCorrections::WriteSingleResponse(TString strInfile, TString strID, TString strOutfile,
2541 Bool_t updateOutfile, TString strOutDir)
2543 // read task ouput from MC and write single track eff - standard dir/list
2545 TString strdir = "PWGJE_FragmentationFunction_" + strID;
2546 TString strlist = "fracfunc_" + strID;
2548 WriteSingleResponse(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir);
2552 //_____________________________________________________________________________________________________________________________________
2553 void AliFragmentationFunctionCorrections::WriteSingleResponse(TString strInfile, TString strdir, TString strlist,
2554 TString strOutfile, Bool_t updateOutfile, TString strOutDir)
2556 // read 2d THnSparse response matrices in pt from file
2558 // write to strOutfile
2560 THnSparse* hnResponseSinglePt;
2561 TH2F* h2ResponseSinglePt;
2563 TFile f(strInfile,"READ");
2566 Printf("%s:%d -- error opening file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2570 if(fDebug>0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
2572 if(strdir && strdir.Length()) gDirectory->cd(strdir);
2576 if(strlist && strlist.Length()){
2578 if(!(list = (TList*) gDirectory->Get(strlist))){
2579 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2584 if(list) hnResponseSinglePt = (THnSparse*) list->FindObject("fhnResponseSinglePt");
2585 else hnResponseSinglePt = (THnSparse*) gDirectory->Get("fhnResponseSinglePt");
2588 if(!hnResponseSinglePt){
2589 Printf("%s:%d -- error retrieving response matrix fhnResponseSinglePt",(char*)__FILE__,__LINE__);
2597 h2ResponseSinglePt = (TH2F*) hnResponseSinglePt->Projection(1,0);// note convention: yDim,xDim
2598 h2ResponseSinglePt->SetNameTitle("h2ResponseSinglePt","");
2603 TString outfileOption = "RECREATE";
2604 if(updateOutfile) outfileOption = "UPDATE";
2606 TFile out(strOutfile,outfileOption);
2609 Printf("%s:%d -- error opening response matrix output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2613 if(fDebug>0) Printf("%s:%d -- write response matrices to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2615 if(strOutDir && strOutDir.Length()){
2618 if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd();
2620 dir = out.mkdir(strOutDir);
2625 hnResponseSinglePt->Write();
2626 h2ResponseSinglePt->Write();
2631 //________________________________________________________________________________________________________________
2632 void AliFragmentationFunctionCorrections::WriteJetTrackEff(TString strInfile, TString strID, TString strOutfile,
2633 Bool_t updateOutfile)
2635 // read task ouput from MC and write single track eff - standard dir/list
2637 TString strdir = "PWGJE_FragmentationFunction_" + strID;
2638 TString strlist = "fracfunc_" + strID;
2640 WriteJetTrackEff(strInfile,strdir,strlist,strOutfile,updateOutfile);
2643 //___________________________________________________________________________________________________________________________________
2644 void AliFragmentationFunctionCorrections::WriteJetTrackEff(TString strInfile, TString strdir, TString strlist,
2645 TString strOutfile, Bool_t updateOutfile)
2647 // read task output from MC and write track eff in jet pt slices as function of pt, z, xi
2648 // argument strlist optional - read from directory strdir if not specified
2649 // write eff to file strOutfile - by default only 'update' file (don't overwrite)
2651 TH1F* hEffPt[fNJetPtSlices];
2652 TH1F* hEffXi[fNJetPtSlices];
2653 TH1F* hEffZ[fNJetPtSlices];
2655 TH1F* hdNdptTracksMCPrimGen[fNJetPtSlices];
2656 TH1F* hdNdxiMCPrimGen[fNJetPtSlices];
2657 TH1F* hdNdzMCPrimGen[fNJetPtSlices];
2659 TH1F* hdNdptTracksMCPrimRec[fNJetPtSlices];
2660 TH1F* hdNdxiMCPrimRec[fNJetPtSlices];
2661 TH1F* hdNdzMCPrimRec[fNJetPtSlices];
2664 TH1F* fh1FFJetPtRecEffGen;
2666 TH2F* fh2FFTrackPtRecEffGen;
2667 TH2F* fh2FFZRecEffGen;
2668 TH2F* fh2FFXiRecEffGen;
2670 TH2F* fh2FFTrackPtRecEffRec;
2671 TH2F* fh2FFZRecEffRec;
2672 TH2F* fh2FFXiRecEffRec;
2675 TFile f(strInfile,"READ");
2678 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2682 if(fDebug>0) Printf("%s:%d -- writeJetTrackEff: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
2684 if(strdir && strdir.Length()) gDirectory->cd(strdir);
2688 if(strlist && strlist.Length()){
2690 if(!(list = (TList*) gDirectory->Get(strlist))){
2691 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2697 fh1FFJetPtRecEffGen = (TH1F*) list->FindObject("fh1FFJetPtRecEffGen");
2699 fh2FFTrackPtRecEffGen = (TH2F*) list->FindObject("fh2FFTrackPtRecEffGen");
2700 fh2FFZRecEffGen = (TH2F*) list->FindObject("fh2FFZRecEffGen");
2701 fh2FFXiRecEffGen = (TH2F*) list->FindObject("fh2FFXiRecEffGen");
2703 fh2FFTrackPtRecEffRec = (TH2F*) list->FindObject("fh2FFTrackPtRecEffRec");
2704 fh2FFZRecEffRec = (TH2F*) list->FindObject("fh2FFZRecEffRec");
2705 fh2FFXiRecEffRec = (TH2F*) list->FindObject("fh2FFXiRecEffRec");
2708 fh1FFJetPtRecEffGen = (TH1F*) gDirectory->Get("fh1FFJetPtRecEffGen");
2710 fh2FFTrackPtRecEffGen = (TH2F*) gDirectory->Get("fh2FFTrackPtRecEffGen");
2711 fh2FFZRecEffGen = (TH2F*) gDirectory->Get("fh2FFZRecEffGen");
2712 fh2FFXiRecEffGen = (TH2F*) gDirectory->Get("fh2FFXiRecEffGen");
2714 fh2FFTrackPtRecEffRec = (TH2F*) gDirectory->Get("fh2FFTrackPtRecEffRec");
2715 fh2FFZRecEffRec = (TH2F*) gDirectory->Get("fh2FFZRecEffRec");
2716 fh2FFXiRecEffRec = (TH2F*) gDirectory->Get("fh2FFXiRecEffRec");
2719 fh1FFJetPtRecEffGen->SetDirectory(0);
2721 fh2FFTrackPtRecEffGen->SetDirectory(0);
2722 fh2FFZRecEffGen->SetDirectory(0);
2723 fh2FFXiRecEffGen->SetDirectory(0);
2725 fh2FFTrackPtRecEffRec->SetDirectory(0);
2726 fh2FFZRecEffRec->SetDirectory(0);
2727 fh2FFXiRecEffRec->SetDirectory(0);
2732 // projections: FF for generated and reconstructed primaries
2734 for(Int_t i=0; i<fNJetPtSlices; i++){
2736 Float_t jetPtLoLim = fJetPtSlices->At(i);
2737 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
2739 Int_t binLo = static_cast<Int_t>(fh2FFTrackPtRecEffGen->GetXaxis()->FindBin(jetPtLoLim));
2740 Int_t binUp = static_cast<Int_t>(fh2FFTrackPtRecEffGen->GetXaxis()->FindBin(jetPtUpLim))-1;
2742 if(binUp > fh2FFTrackPtRecEffGen->GetNbinsX()){
2743 Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
2747 TString strNameFFPtGen(Form("fh1FFTrackPtGenPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2748 TString strNameFFZGen(Form("fh1FFZGenPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2749 TString strNameFFXiGen(Form("fh1FFXiGenPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2751 TString strNameFFPtRec(Form("fh1FFTrackPtRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2752 TString strNameFFZRec(Form("fh1FFZRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2753 TString strNameFFXiRec(Form("fh1FFXiRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2756 // appendix 'unbinned' to avoid histos with same name after rebinning
2758 hdNdptTracksMCPrimGen[i] = (TH1F*) fh2FFTrackPtRecEffGen->ProjectionY(strNameFFPtGen+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
2759 hdNdzMCPrimGen[i] = (TH1F*) fh2FFZRecEffGen->ProjectionY(strNameFFZGen+"_unBinned",binLo,binUp,"o");
2760 hdNdxiMCPrimGen[i] = (TH1F*) fh2FFXiRecEffGen->ProjectionY(strNameFFXiGen+"_unBinned",binLo,binUp,"o");
2762 hdNdptTracksMCPrimRec[i] = (TH1F*) fh2FFTrackPtRecEffRec->ProjectionY(strNameFFPtRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
2763 hdNdzMCPrimRec[i] = (TH1F*) fh2FFZRecEffRec->ProjectionY(strNameFFZRec+"_unBinned",binLo,binUp,"o");
2764 hdNdxiMCPrimRec[i] = (TH1F*) fh2FFXiRecEffRec->ProjectionY(strNameFFXiRec+"_unBinned",binLo,binUp,"o");
2768 if(fNHistoBinsPt[i]) hdNdptTracksMCPrimGen[i] = (TH1F*) hdNdptTracksMCPrimGen[i]->Rebin(fNHistoBinsPt[i],strNameFFPtGen,fHistoBinsPt[i]->GetArray());
2769 if(fNHistoBinsZ[i]) hdNdzMCPrimGen[i] = (TH1F*) hdNdzMCPrimGen[i]->Rebin(fNHistoBinsZ[i],strNameFFZGen,fHistoBinsZ[i]->GetArray());
2770 if(fNHistoBinsXi[i]) hdNdxiMCPrimGen[i] = (TH1F*) hdNdxiMCPrimGen[i]->Rebin(fNHistoBinsXi[i],strNameFFXiGen,fHistoBinsXi[i]->GetArray());
2772 if(fNHistoBinsPt[i]) hdNdptTracksMCPrimRec[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtRec,fHistoBinsPt[i]->GetArray());
2773 if(fNHistoBinsZ[i]) hdNdzMCPrimRec[i] = (TH1F*) hdNdzMCPrimRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZRec,fHistoBinsZ[i]->GetArray());
2774 if(fNHistoBinsXi[i]) hdNdxiMCPrimRec[i] = (TH1F*) hdNdxiMCPrimRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiRec,fHistoBinsXi[i]->GetArray());
2776 hdNdptTracksMCPrimGen[i]->SetNameTitle(strNameFFPtGen,"");
2777 hdNdzMCPrimGen[i]->SetNameTitle(strNameFFZGen,"");
2778 hdNdxiMCPrimGen[i]->SetNameTitle(strNameFFXiGen,"");
2780 hdNdptTracksMCPrimRec[i]->SetNameTitle(strNameFFPtRec,"");
2781 hdNdzMCPrimRec[i]->SetNameTitle(strNameFFZRec,"");
2782 hdNdxiMCPrimRec[i]->SetNameTitle(strNameFFXiRec,"");
2786 Double_t nJetsBin = fh1FFJetPtRecEffGen->Integral(binLo,binUp);
2788 NormalizeTH1(hdNdptTracksMCPrimGen[i],nJetsBin);
2789 NormalizeTH1(hdNdzMCPrimGen[i],nJetsBin);
2790 NormalizeTH1(hdNdxiMCPrimGen[i],nJetsBin);
2792 NormalizeTH1(hdNdptTracksMCPrimRec[i],nJetsBin);
2793 NormalizeTH1(hdNdzMCPrimRec[i],nJetsBin);
2794 NormalizeTH1(hdNdxiMCPrimRec[i],nJetsBin);
2796 // divide rec/gen : efficiency
2798 TString strNameEffPt(Form("hEffPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2799 TString strNameEffZ(Form("hEffZ_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2800 TString strNameEffXi(Form("hEffXi_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2802 hEffPt[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Clone(strNameEffPt);
2803 hEffPt[i]->Divide(hdNdptTracksMCPrimRec[i],hdNdptTracksMCPrimGen[i],1,1,"B"); // binominal errors
2805 hEffXi[i] = (TH1F*) hdNdxiMCPrimRec[i]->Clone(strNameEffXi);
2806 hEffXi[i]->Divide(hdNdxiMCPrimRec[i],hdNdxiMCPrimGen[i],1,1,"B"); // binominal errors
2808 hEffZ[i] = (TH1F*) hdNdzMCPrimRec[i]->Clone(strNameEffZ);
2809 hEffZ[i]->Divide(hdNdzMCPrimRec[i],hdNdzMCPrimGen[i],1,1,"B"); // binominal errors
2814 TString outfileOption = "RECREATE";
2815 if(updateOutfile) outfileOption = "UPDATE";
2817 TFile out(strOutfile,outfileOption);
2820 Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2824 if(fDebug>0) Printf("%s:%d -- write efficiency to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2826 // if(strdir && strdir.Length()){
2827 // TDirectory* dir = out.mkdir(strdir);
2831 for(Int_t i=0; i<fNJetPtSlices; i++){
2833 hdNdptTracksMCPrimGen[i]->Write();
2834 hdNdxiMCPrimGen[i]->Write();
2835 hdNdzMCPrimGen[i]->Write();
2837 hdNdptTracksMCPrimRec[i]->Write();
2838 hdNdxiMCPrimRec[i]->Write();
2839 hdNdzMCPrimRec[i]->Write();
2849 delete fh1FFJetPtRecEffGen;
2851 delete fh2FFTrackPtRecEffGen;
2852 delete fh2FFZRecEffGen;
2853 delete fh2FFXiRecEffGen;
2855 delete fh2FFTrackPtRecEffRec;
2856 delete fh2FFZRecEffRec;
2857 delete fh2FFXiRecEffRec;
2860 //________________________________________________________________________________________________________________
2861 void AliFragmentationFunctionCorrections::WriteJetSecCorr(TString strInfile, TString strID, TString strOutfile,
2862 Bool_t updateOutfile, TString strOutDir)
2864 // read task ouput from MC and write secondary correction - standard dir/list
2866 TString strdir = "PWGJE_FragmentationFunction_" + strID;
2867 TString strlist = "fracfunc_" + strID;
2869 WriteJetSecCorr(strInfile,strdir,strlist,strOutfile,updateOutfile, strOutDir,kFALSE,"");
2872 //________________________________________________________________________________________________________________
2873 void AliFragmentationFunctionCorrections::WriteBgrJetSecCorr(TString strInfile, TString strBgrID, TString strID, TString strOutfile,
2874 Bool_t updateOutfile, TString strOutDir)
2876 // read task ouput from MC and write secondary correction - standard dir/list
2878 TString strdir = "PWGJE_FragmentationFunction_" + strID;
2879 TString strlist = "fracfunc_" + strID;
2881 WriteJetSecCorr(strInfile,strdir,strlist,strOutfile,updateOutfile, strOutDir,kTRUE,strBgrID);
2885 //___________________________________________________________________________________________________________________________________
2886 void AliFragmentationFunctionCorrections::WriteJetSecCorr(TString strInfile, TString strdir, TString strlist,
2887 TString strOutfile, Bool_t updateOutfile, TString strOutDir,
2888 Bool_t writeBgr, TString strBgrID)
2890 // read task output from MC and write secondary correction in jet pt slices as function of pt, z, xi
2891 // argument strlist optional - read from directory strdir if not specified
2892 // write eff to file strOutfile - by default only 'update' file (don't overwrite)
2894 TH1F* hSecCorrPt[fNJetPtSlices];
2895 TH1F* hSecCorrXi[fNJetPtSlices];
2896 TH1F* hSecCorrZ[fNJetPtSlices];
2898 // corr factors using naive Pythia strangeness - for sys uncertainties
2899 TH1F* hSecCorrPt_nonSc[fNJetPtSlices];
2900 TH1F* hSecCorrXi_nonSc[fNJetPtSlices];
2901 TH1F* hSecCorrZ_nonSc[fNJetPtSlices];
2903 TH1F* hdNdptTracksMCPrimRec[fNJetPtSlices];
2904 TH1F* hdNdxiMCPrimRec[fNJetPtSlices];
2905 TH1F* hdNdzMCPrimRec[fNJetPtSlices];
2907 TH1F* hdNdptTracksMCSecRecNS[fNJetPtSlices];
2908 TH1F* hdNdxiMCSecRecNS[fNJetPtSlices];
2909 TH1F* hdNdzMCSecRecNS[fNJetPtSlices];
2911 TH1F* hdNdptTracksMCSecRecS[fNJetPtSlices];
2912 TH1F* hdNdxiMCSecRecS[fNJetPtSlices];
2913 TH1F* hdNdzMCSecRecS[fNJetPtSlices];
2915 TH1F* hdNdptTracksMCSecRecSsc[fNJetPtSlices];
2916 TH1F* hdNdxiMCSecRecSsc[fNJetPtSlices];
2917 TH1F* hdNdzMCSecRecSsc[fNJetPtSlices];
2921 TH1F* fh1FFJetPtRecEffRec;
2923 TH2F* fh2FFTrackPtRecEffRec;
2924 TH2F* fh2FFZRecEffRec;
2925 TH2F* fh2FFXiRecEffRec;
2927 TH2F* fh2FFTrackPtSecRecNS;
2928 TH2F* fh2FFZSecRecNS;
2929 TH2F* fh2FFXiSecRecNS;
2931 TH2F* fh2FFTrackPtSecRecS;
2932 TH2F* fh2FFZSecRecS;
2933 TH2F* fh2FFXiSecRecS;
2935 TH2F* fh2FFTrackPtSecRecSsc;
2936 TH2F* fh2FFZSecRecSsc;
2937 TH2F* fh2FFXiSecRecSsc;
2941 TFile f(strInfile,"READ");
2944 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2948 if(fDebug>0) Printf("%s:%d -- writeJetTrackSecCorr: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
2950 if(strdir && strdir.Length()) gDirectory->cd(strdir);
2954 if(strlist && strlist.Length()){
2956 if(!(list = (TList*) gDirectory->Get(strlist))){
2957 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2963 fh1FFJetPtRecEffRec = (TH1F*) list->FindObject("fh1FFJetPtRecEffRec");
2966 fh2FFTrackPtRecEffRec = (TH2F*) list->FindObject(Form("fh2FFTrackPt%sRecEffRec",strBgrID.Data()));
2967 fh2FFZRecEffRec = (TH2F*) list->FindObject(Form("fh2FFZ%sRecEffRec",strBgrID.Data()));
2968 fh2FFXiRecEffRec = (TH2F*) list->FindObject(Form("fh2FFXi%sRecEffRec",strBgrID.Data()));
2970 fh2FFTrackPtSecRecNS = (TH2F*) list->FindObject(Form("fh2FFTrackPt%sSecRecNS",strBgrID.Data()));
2971 fh2FFZSecRecNS = (TH2F*) list->FindObject(Form("fh2FFZ%sSecRecNS",strBgrID.Data()));
2972 fh2FFXiSecRecNS = (TH2F*) list->FindObject(Form("fh2FFXi%sSecRecNS",strBgrID.Data()));
2974 fh2FFTrackPtSecRecS = (TH2F*) list->FindObject(Form("fh2FFTrackPt%sSecRecS",strBgrID.Data()));
2975 fh2FFZSecRecS = (TH2F*) list->FindObject(Form("fh2FFZ%sSecRecS",strBgrID.Data()));
2976 fh2FFXiSecRecS = (TH2F*) list->FindObject(Form("fh2FFXi%sSecRecS",strBgrID.Data()));
2978 fh2FFTrackPtSecRecSsc = (TH2F*) list->FindObject(Form("fh2FFTrackPt%sSecRecSsc",strBgrID.Data()));
2979 fh2FFZSecRecSsc = (TH2F*) list->FindObject(Form("fh2FFZ%sSecRecSsc",strBgrID.Data()));
2980 fh2FFXiSecRecSsc = (TH2F*) list->FindObject(Form("fh2FFXi%sSecRecSsc",strBgrID.Data()));
2983 fh2FFTrackPtRecEffRec = (TH2F*) list->FindObject("fh2FFTrackPtRecEffRec");
2984 fh2FFZRecEffRec = (TH2F*) list->FindObject("fh2FFZRecEffRec");
2985 fh2FFXiRecEffRec = (TH2F*) list->FindObject("fh2FFXiRecEffRec");
2987 fh2FFTrackPtSecRecNS = (TH2F*) list->FindObject("fh2FFTrackPtSecRecNS");
2988 fh2FFZSecRecNS = (TH2F*) list->FindObject("fh2FFZSecRecNS");
2989 fh2FFXiSecRecNS = (TH2F*) list->FindObject("fh2FFXiSecRecNS");
2991 fh2FFTrackPtSecRecS = (TH2F*) list->FindObject("fh2FFTrackPtSecRecS");
2992 fh2FFZSecRecS = (TH2F*) list->FindObject("fh2FFZSecRecS");
2993 fh2FFXiSecRecS = (TH2F*) list->FindObject("fh2FFXiSecRecS");
2995 fh2FFTrackPtSecRecSsc = (TH2F*) list->FindObject("fh2FFTrackPtSecRecSsc");
2996 fh2FFZSecRecSsc = (TH2F*) list->FindObject("fh2FFZSecRecSsc");
2997 fh2FFXiSecRecSsc = (TH2F*) list->FindObject("fh2FFXiSecRecSsc");
3001 fh1FFJetPtRecEffRec = (TH1F*) gDirectory->Get("fh1FFJetPtRecEffRec");
3004 fh2FFTrackPtRecEffRec = (TH2F*) gDirectory->Get(Form("fh2FFTrackPt%sRecEffRec",strBgrID.Data()));
3005 fh2FFZRecEffRec = (TH2F*) gDirectory->Get(Form("fh2FFZ%sRecEffRec",strBgrID.Data()));
3006 fh2FFXiRecEffRec = (TH2F*) gDirectory->Get(Form("fh2FFXi%sRecEffRec",strBgrID.Data()));
3008 fh2FFTrackPtSecRecNS = (TH2F*) gDirectory->Get(Form("fh2FFTrackPt%sSecRecNS",strBgrID.Data()));
3009 fh2FFZSecRecNS = (TH2F*) gDirectory->Get(Form("fh2FFZ%sSecRecNS",strBgrID.Data()));
3010 fh2FFXiSecRecNS = (TH2F*) gDirectory->Get(Form("fh2FFXi%sSecRecNS",strBgrID.Data()));
3012 fh2FFTrackPtSecRecS = (TH2F*) gDirectory->Get(Form("fh2FFTrackPt%sSecRecS",strBgrID.Data()));
3013 fh2FFZSecRecS = (TH2F*) gDirectory->Get(Form("fh2FFZ%sSecRecS",strBgrID.Data()));
3014 fh2FFXiSecRecS = (TH2F*) gDirectory->Get(Form("fh2FFXi%sSecRecS",strBgrID.Data()));
3016 fh2FFTrackPtSecRecSsc = (TH2F*) gDirectory->Get(Form("fh2FFTrackPt%sSecRecSsc",strBgrID.Data()));
3017 fh2FFZSecRecSsc = (TH2F*) gDirectory->Get(Form("fh2FFZ%sSecRecSsc",strBgrID.Data()));
3018 fh2FFXiSecRecSsc = (TH2F*) gDirectory->Get(Form("fh2FFXi%sSecRecSsc",strBgrID.Data()));
3021 fh2FFTrackPtRecEffRec = (TH2F*) gDirectory->Get("fh2FFTrackPtRecEffRec");
3022 fh2FFZRecEffRec = (TH2F*) gDirectory->Get("fh2FFZRecEffRec");
3023 fh2FFXiRecEffRec = (TH2F*) gDirectory->Get("fh2FFXiRecEffRec");
3025 fh2FFTrackPtSecRecNS = (TH2F*) gDirectory->Get("fh2FFTrackPtSecRecNS");
3026 fh2FFZSecRecNS = (TH2F*) gDirectory->Get("fh2FFZSecRecNS");
3027 fh2FFXiSecRecNS = (TH2F*) gDirectory->Get("fh2FFXiSecRecNS");
3029 fh2FFTrackPtSecRecS = (TH2F*) gDirectory->Get("fh2FFTrackPtSecRecS");
3030 fh2FFZSecRecS = (TH2F*) gDirectory->Get("fh2FFZSecRecS");
3031 fh2FFXiSecRecS = (TH2F*) gDirectory->Get("fh2FFXiSecRecS");
3033 fh2FFTrackPtSecRecSsc = (TH2F*) gDirectory->Get("fh2FFTrackPtSecRecSsc");
3034 fh2FFZSecRecSsc = (TH2F*) gDirectory->Get("fh2FFZSecRecSsc");
3035 fh2FFXiSecRecSsc = (TH2F*) gDirectory->Get("fh2FFXiSecRecSsc");
3039 fh1FFJetPtRecEffRec->SetDirectory(0);
3041 fh2FFTrackPtRecEffRec->SetDirectory(0);
3042 fh2FFZRecEffRec->SetDirectory(0);
3043 fh2FFXiRecEffRec->SetDirectory(0);
3045 fh2FFTrackPtSecRecNS->SetDirectory(0);
3046 fh2FFZSecRecNS->SetDirectory(0);
3047 fh2FFXiSecRecNS->SetDirectory(0);
3049 fh2FFTrackPtSecRecS->SetDirectory(0);
3050 fh2FFZSecRecS->SetDirectory(0);
3051 fh2FFXiSecRecS->SetDirectory(0);
3053 fh2FFTrackPtSecRecSsc->SetDirectory(0);
3054 fh2FFZSecRecSsc->SetDirectory(0);
3055 fh2FFXiSecRecSsc->SetDirectory(0);
3060 // projections: FF for generated and reconstructed primaries
3062 for(Int_t i=0; i<fNJetPtSlices; i++){
3064 Float_t jetPtLoLim = fJetPtSlices->At(i);
3065 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
3067 Int_t binLo = static_cast<Int_t>(fh2FFTrackPtRecEffRec->GetXaxis()->FindBin(jetPtLoLim));
3068 Int_t binUp = static_cast<Int_t>(fh2FFTrackPtRecEffRec->GetXaxis()->FindBin(jetPtUpLim))-1;
3070 if(binUp > fh2FFTrackPtRecEffRec->GetNbinsX()){
3071 Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
3075 TString strNameFFPtPrimRec(Form("fh1FFTrackPtRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3076 TString strNameFFZPrimRec(Form("fh1FFZRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3077 TString strNameFFXiPrimRec(Form("fh1FFXiRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3079 TString strNameFFPtSecRecNS(Form("fh1FFTrackPtRecSecNS_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3080 TString strNameFFZSecRecNS(Form("fh1FFZRecSecNS_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3081 TString strNameFFXiSecRecNS(Form("fh1FFXiRecSecNS_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3083 TString strNameFFPtSecRecS(Form("fh1FFTrackPtRecSecS_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3084 TString strNameFFZSecRecS(Form("fh1FFZRecSecS_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3085 TString strNameFFXiSecRecS(Form("fh1FFXiRecSecS_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3087 TString strNameFFPtSecRecSsc(Form("fh1FFTrackPtRecSecSsc_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3088 TString strNameFFZSecRecSsc(Form("fh1FFZRecSecSsc_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3089 TString strNameFFXiSecRecSsc(Form("fh1FFXiRecSecSsc_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3092 strNameFFPtPrimRec.Form("fh1BgrTrackPt%sRecPrim_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3093 strNameFFZPrimRec.Form("fh1BgrZ%sRecPrim_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3094 strNameFFXiPrimRec.Form("fh1BgrXi%sRecPrim_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3096 strNameFFPtSecRecNS.Form("fh1BgrTrackPt%sRecSecNS_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3097 strNameFFZSecRecNS.Form("fh1BgrZ%sRecSecNS_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3098 strNameFFXiSecRecNS.Form("fh1BgrXi%sRecSecNS_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3100 strNameFFPtSecRecS.Form("fh1BgrTrackPt%sRecSecS_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3101 strNameFFZSecRecS.Form("fh1BgrZ%sRecSecS_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3102 strNameFFXiSecRecS.Form("fh1BgrXi%sRecSecS_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3104 strNameFFPtSecRecSsc.Form("fh1BgrTrackPt%sRecSecSsc_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3105 strNameFFZSecRecSsc.Form("fh1BgrZ%sRecSecSsc_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3106 strNameFFXiSecRecSsc.Form("fh1BgrXi%sRecSecSsc_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3110 // appendix 'unbinned' to avoid histos with same name after rebinning
3112 hdNdptTracksMCPrimRec[i] = (TH1F*) fh2FFTrackPtRecEffRec->ProjectionY(strNameFFPtPrimRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
3113 hdNdzMCPrimRec[i] = (TH1F*) fh2FFZRecEffRec->ProjectionY(strNameFFZPrimRec+"_unBinned",binLo,binUp,"o");
3114 hdNdxiMCPrimRec[i] = (TH1F*) fh2FFXiRecEffRec->ProjectionY(strNameFFXiPrimRec+"_unBinned",binLo,binUp,"o");
3116 hdNdptTracksMCSecRecNS[i] = (TH1F*) fh2FFTrackPtSecRecNS->ProjectionY(strNameFFPtSecRecNS+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
3117 hdNdzMCSecRecNS[i] = (TH1F*) fh2FFZSecRecNS->ProjectionY(strNameFFZSecRecNS+"_unBinned",binLo,binUp,"o");
3118 hdNdxiMCSecRecNS[i] = (TH1F*) fh2FFXiSecRecNS->ProjectionY(strNameFFXiSecRecNS+"_unBinned",binLo,binUp,"o");
3120 hdNdptTracksMCSecRecS[i] = (TH1F*) fh2FFTrackPtSecRecS->ProjectionY(strNameFFPtSecRecS+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
3121 hdNdzMCSecRecS[i] = (TH1F*) fh2FFZSecRecS->ProjectionY(strNameFFZSecRecS+"_unBinned",binLo,binUp,"o");
3122 hdNdxiMCSecRecS[i] = (TH1F*) fh2FFXiSecRecS->ProjectionY(strNameFFXiSecRecS+"_unBinned",binLo,binUp,"o");
3124 hdNdptTracksMCSecRecSsc[i] = (TH1F*) fh2FFTrackPtSecRecSsc->ProjectionY(strNameFFPtSecRecSsc+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
3125 hdNdzMCSecRecSsc[i] = (TH1F*) fh2FFZSecRecSsc->ProjectionY(strNameFFZSecRecSsc+"_unBinned",binLo,binUp,"o");
3126 hdNdxiMCSecRecSsc[i] = (TH1F*) fh2FFXiSecRecSsc->ProjectionY(strNameFFXiSecRecSsc+"_unBinned",binLo,binUp,"o");
3130 if(fNHistoBinsPt[i]) hdNdptTracksMCPrimRec[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtPrimRec,fHistoBinsPt[i]->GetArray());
3131 if(fNHistoBinsZ[i]) hdNdzMCPrimRec[i] = (TH1F*) hdNdzMCPrimRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZPrimRec,fHistoBinsZ[i]->GetArray());
3132 if(fNHistoBinsXi[i]) hdNdxiMCPrimRec[i] = (TH1F*) hdNdxiMCPrimRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiPrimRec,fHistoBinsXi[i]->GetArray());
3134 if(fNHistoBinsPt[i]) hdNdptTracksMCSecRecNS[i] = (TH1F*) hdNdptTracksMCSecRecNS[i]->Rebin(fNHistoBinsPt[i],strNameFFPtSecRecNS,fHistoBinsPt[i]->GetArray());
3135 if(fNHistoBinsZ[i]) hdNdzMCSecRecNS[i] = (TH1F*) hdNdzMCSecRecNS[i]->Rebin(fNHistoBinsZ[i],strNameFFZSecRecNS,fHistoBinsZ[i]->GetArray());
3136 if(fNHistoBinsXi[i]) hdNdxiMCSecRecNS[i] = (TH1F*) hdNdxiMCSecRecNS[i]->Rebin(fNHistoBinsXi[i],strNameFFXiSecRecNS,fHistoBinsXi[i]->GetArray());
3138 if(fNHistoBinsPt[i]) hdNdptTracksMCSecRecS[i] = (TH1F*) hdNdptTracksMCSecRecS[i]->Rebin(fNHistoBinsPt[i],strNameFFPtSecRecS,fHistoBinsPt[i]->GetArray());
3139 if(fNHistoBinsZ[i]) hdNdzMCSecRecS[i] = (TH1F*) hdNdzMCSecRecS[i]->Rebin(fNHistoBinsZ[i],strNameFFZSecRecS,fHistoBinsZ[i]->GetArray());
3140 if(fNHistoBinsXi[i]) hdNdxiMCSecRecS[i] = (TH1F*) hdNdxiMCSecRecS[i]->Rebin(fNHistoBinsXi[i],strNameFFXiSecRecS,fHistoBinsXi[i]->GetArray());
3142 if(fNHistoBinsPt[i]) hdNdptTracksMCSecRecSsc[i] = (TH1F*) hdNdptTracksMCSecRecSsc[i]->Rebin(fNHistoBinsPt[i],strNameFFPtSecRecSsc,fHistoBinsPt[i]->GetArray());
3143 if(fNHistoBinsZ[i]) hdNdzMCSecRecSsc[i] = (TH1F*) hdNdzMCSecRecSsc[i]->Rebin(fNHistoBinsZ[i],strNameFFZSecRecSsc,fHistoBinsZ[i]->GetArray());
3144 if(fNHistoBinsXi[i]) hdNdxiMCSecRecSsc[i] = (TH1F*) hdNdxiMCSecRecSsc[i]->Rebin(fNHistoBinsXi[i],strNameFFXiSecRecSsc,fHistoBinsXi[i]->GetArray());
3147 hdNdptTracksMCPrimRec[i]->SetNameTitle(strNameFFPtPrimRec,"");
3148 hdNdzMCPrimRec[i]->SetNameTitle(strNameFFZPrimRec,"");
3149 hdNdxiMCPrimRec[i]->SetNameTitle(strNameFFXiPrimRec,"");
3151 hdNdptTracksMCSecRecNS[i]->SetNameTitle(strNameFFPtSecRecNS,"");
3152 hdNdzMCSecRecNS[i]->SetNameTitle(strNameFFZSecRecNS,"");
3153 hdNdxiMCSecRecNS[i]->SetNameTitle(strNameFFXiSecRecNS,"");
3155 hdNdptTracksMCSecRecS[i]->SetNameTitle(strNameFFPtSecRecS,"");
3156 hdNdzMCSecRecS[i]->SetNameTitle(strNameFFZSecRecS,"");
3157 hdNdxiMCSecRecS[i]->SetNameTitle(strNameFFXiSecRecS,"");
3159 hdNdptTracksMCSecRecSsc[i]->SetNameTitle(strNameFFPtSecRecSsc,"");
3160 hdNdzMCSecRecSsc[i]->SetNameTitle(strNameFFZSecRecSsc,"");
3161 hdNdxiMCSecRecSsc[i]->SetNameTitle(strNameFFXiSecRecSsc,"");
3164 Double_t nJetsBin = fh1FFJetPtRecEffRec->Integral(binLo,binUp);
3166 NormalizeTH1(hdNdptTracksMCPrimRec[i],nJetsBin);
3167 NormalizeTH1(hdNdzMCPrimRec[i],nJetsBin);
3168 NormalizeTH1(hdNdxiMCPrimRec[i],nJetsBin);
3170 NormalizeTH1(hdNdptTracksMCSecRecNS[i],nJetsBin);
3171 NormalizeTH1(hdNdzMCSecRecNS[i],nJetsBin);
3172 NormalizeTH1(hdNdxiMCSecRecNS[i],nJetsBin);
3174 NormalizeTH1(hdNdptTracksMCSecRecS[i],nJetsBin);
3175 NormalizeTH1(hdNdzMCSecRecS[i],nJetsBin);
3176 NormalizeTH1(hdNdxiMCSecRecS[i],nJetsBin);
3178 NormalizeTH1(hdNdptTracksMCSecRecSsc[i],nJetsBin);
3179 NormalizeTH1(hdNdzMCSecRecSsc[i],nJetsBin);
3180 NormalizeTH1(hdNdxiMCSecRecSsc[i],nJetsBin);
3183 // divide prim / (prim+sec) : corr factor
3184 TString strNameSecCorrPt(Form("hSecCorrPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3185 TString strNameSecCorrZ(Form("hSecCorrZ_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3186 TString strNameSecCorrXi(Form("hSecCorrXi_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3189 strNameSecCorrPt.Form("hSecCorrBgrPt_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
3190 strNameSecCorrZ.Form("hSecCorrBgrZ_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
3191 strNameSecCorrXi.Form("hSecCorrBgrXi_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
3194 hSecCorrPt[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Clone(strNameSecCorrPt);
3195 TH1F* hSumPrimSecPt = (TH1F*) hdNdptTracksMCPrimRec[i]->Clone("hSumPrimSecPt");
3196 hSumPrimSecPt->Add(hdNdptTracksMCSecRecNS[i]);
3197 hSumPrimSecPt->Add(hdNdptTracksMCSecRecSsc[i]);
3198 hSecCorrPt[i]->Divide(hdNdptTracksMCPrimRec[i],hSumPrimSecPt,1,1,"B"); // binominal errors
3200 hSecCorrXi[i] = (TH1F*) hdNdxiMCPrimRec[i]->Clone(strNameSecCorrXi);
3201 TH1F* hSumPrimSecXi = (TH1F*) hdNdxiMCPrimRec[i]->Clone("hSumPrimSecXi");
3202 hSumPrimSecXi->Add(hdNdxiMCSecRecNS[i]);
3203 hSumPrimSecXi->Add(hdNdxiMCSecRecSsc[i]);
3204 hSecCorrXi[i]->Divide(hdNdxiMCPrimRec[i],hSumPrimSecXi,1,1,"B"); // binominal errors
3206 hSecCorrZ[i] = (TH1F*) hdNdzMCPrimRec[i]->Clone(strNameSecCorrZ);
3207 TH1F* hSumPrimSecZ = (TH1F*) hdNdzMCPrimRec[i]->Clone("hSumPrimSecZ");
3208 hSumPrimSecZ->Add(hdNdzMCSecRecNS[i]);
3209 hSumPrimSecZ->Add(hdNdzMCSecRecSsc[i]);
3210 hSecCorrZ[i]->Divide(hdNdzMCPrimRec[i],hSumPrimSecZ,1,1,"B"); // binominal errors
3212 // the same using unscaled strangeness
3213 TString strNameSecCorrPt_nonSc(Form("hSecCorrPt_nonSc_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3214 TString strNameSecCorrZ_nonSc(Form("hSecCorrZ_nonSc_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3215 TString strNameSecCorrXi_nonSc(Form("hSecCorrXi_nonSc_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3218 strNameSecCorrPt_nonSc.Form("hSecCorrBgrPt_nonSc_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
3219 strNameSecCorrZ_nonSc.Form("hSecCorrBgrZ_nonSc_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
3220 strNameSecCorrXi_nonSc.Form("hSecCorrBgrXi_nonSc_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
3223 hSecCorrPt_nonSc[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Clone(strNameSecCorrPt_nonSc);
3224 TH1F* hSumPrimSecPt_nonSc = (TH1F*) hdNdptTracksMCPrimRec[i]->Clone("hSumPrimSecPt_nonSc");
3225 hSumPrimSecPt_nonSc->Add(hdNdptTracksMCSecRecNS[i]);
3226 hSumPrimSecPt_nonSc->Add(hdNdptTracksMCSecRecS[i]); // non-scaled secondaries from strangeness
3227 hSecCorrPt_nonSc[i]->Divide(hdNdptTracksMCPrimRec[i],hSumPrimSecPt_nonSc,1,1,"B"); // binominal errors
3229 hSecCorrZ_nonSc[i] = (TH1F*) hdNdzMCPrimRec[i]->Clone(strNameSecCorrZ_nonSc);
3230 TH1F* hSumPrimSecZ_nonSc = (TH1F*) hdNdzMCPrimRec[i]->Clone("hSumPrimSecZ_nonSc");
3231 hSumPrimSecZ_nonSc->Add(hdNdzMCSecRecNS[i]);
3232 hSumPrimSecZ_nonSc->Add(hdNdzMCSecRecS[i]); // non-scaled secondaries from strangeness
3233 hSecCorrZ_nonSc[i]->Divide(hdNdzMCPrimRec[i],hSumPrimSecZ_nonSc,1,1,"B"); // binominal errors
3235 hSecCorrXi_nonSc[i] = (TH1F*) hdNdxiMCPrimRec[i]->Clone(strNameSecCorrXi_nonSc);
3236 TH1F* hSumPrimSecXi_nonSc = (TH1F*) hdNdxiMCPrimRec[i]->Clone("hSumPrimSecXi_nonSc");
3237 hSumPrimSecXi_nonSc->Add(hdNdxiMCSecRecNS[i]);
3238 hSumPrimSecXi_nonSc->Add(hdNdxiMCSecRecS[i]); // non-scaled secondaries from strangeness
3239 hSecCorrXi_nonSc[i]->Divide(hdNdxiMCPrimRec[i],hSumPrimSecXi_nonSc,1,1,"B"); // binominal errors
3244 TString outfileOption = "RECREATE";
3245 if(updateOutfile) outfileOption = "UPDATE";
3247 TFile out(strOutfile,outfileOption);
3250 Printf("%s:%d -- error opening sec corr output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
3254 if(fDebug>0) Printf("%s:%d -- write jet sec corr to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
3256 if(strOutDir && strOutDir.Length()){
3259 if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd();
3261 dir = out.mkdir(strOutDir);
3267 for(Int_t i=0; i<fNJetPtSlices; i++){
3269 hSecCorrPt[i]->Write();
3270 hSecCorrXi[i]->Write();
3271 hSecCorrZ[i]->Write();
3273 hSecCorrPt_nonSc[i]->Write();
3274 hSecCorrXi_nonSc[i]->Write();
3275 hSecCorrZ_nonSc[i]->Write();
3278 TString strSpectraDir = "spectraSecCorr";
3279 if(writeBgr) strSpectraDir = "spectraBgrSecCorr";
3281 TDirectory *dOut = gDirectory->mkdir(strSpectraDir);
3284 for(Int_t i=0; i<fNJetPtSlices; i++){
3286 hdNdptTracksMCPrimRec[i]->Write();
3287 hdNdzMCPrimRec[i]->Write();
3288 hdNdxiMCPrimRec[i]->Write();
3290 hdNdptTracksMCSecRecNS[i]->Write();
3291 hdNdzMCSecRecNS[i]->Write();
3292 hdNdxiMCSecRecNS[i]->Write();
3294 hdNdptTracksMCSecRecS[i]->Write();
3295 hdNdzMCSecRecS[i]->Write();
3296 hdNdxiMCSecRecS[i]->Write();
3298 hdNdptTracksMCSecRecSsc[i]->Write();
3299 hdNdzMCSecRecSsc[i]->Write();
3300 hdNdxiMCSecRecSsc[i]->Write();
3305 delete fh1FFJetPtRecEffRec;
3307 delete fh2FFTrackPtRecEffRec;
3308 delete fh2FFZRecEffRec;
3309 delete fh2FFXiRecEffRec;
3311 delete fh2FFTrackPtSecRecNS;
3312 delete fh2FFZSecRecNS;
3313 delete fh2FFXiSecRecNS;
3315 delete fh2FFTrackPtSecRecS;
3316 delete fh2FFZSecRecS;
3317 delete fh2FFXiSecRecS;
3319 delete fh2FFTrackPtSecRecSsc;
3320 delete fh2FFZSecRecSsc;
3321 delete fh2FFXiSecRecSsc;
3324 //________________________________________________________________________________________________________________
3325 void AliFragmentationFunctionCorrections::WriteJetResponse(TString strInfile, TString strID, TString strOutfile,
3326 Bool_t updateOutfile, TString strOutDir )
3328 // read task ouput from MC and write single track eff - standard dir/list
3330 TString strdir = "PWGJE_FragmentationFunction_" + strID;
3331 TString strlist = "fracfunc_" + strID;
3333 WriteJetResponse(strInfile,strdir,strlist,strOutfile,updateOutfile, strOutDir);
3336 //_____________________________________________________________________________________________________________________________________
3337 void AliFragmentationFunctionCorrections::WriteJetResponse(TString strInfile, TString strdir, TString strlist,
3338 TString strOutfile, Bool_t updateOutfile, TString strOutDir)
3340 // read 3d THnSparse response matrices in pt,z,xi vs jet pt from file
3341 // project THnSparse + TH2 in jet pt slices
3342 // write to strOutfile
3344 THnSparse* hn3ResponseJetPt;
3345 THnSparse* hn3ResponseJetZ;
3346 THnSparse* hn3ResponseJetXi;
3348 // 2D response matrices
3350 THnSparse* hnResponsePt[fNJetPtSlices];
3351 THnSparse* hnResponseZ[fNJetPtSlices];
3352 THnSparse* hnResponseXi[fNJetPtSlices];
3354 TH2F* h2ResponsePt[fNJetPtSlices];
3355 TH2F* h2ResponseZ[fNJetPtSlices];
3356 TH2F* h2ResponseXi[fNJetPtSlices];
3358 // 1D projections on gen pt / rec pt axes
3360 TH1F* h1FFPtRec[fNJetPtSlices];
3361 TH1F* h1FFZRec[fNJetPtSlices];
3362 TH1F* h1FFXiRec[fNJetPtSlices];
3364 TH1F* h1FFPtGen[fNJetPtSlices];
3365 TH1F* h1FFZGen[fNJetPtSlices];
3366 TH1F* h1FFXiGen[fNJetPtSlices];
3368 TH1F* h1RatioPt[fNJetPtSlices];
3369 TH1F* h1RatioZ[fNJetPtSlices];
3370 TH1F* h1RatioXi[fNJetPtSlices];
3374 TFile f(strInfile,"READ");
3377 Printf("%s:%d -- error opening file %s", (char*)__FILE__,__LINE__,strInfile.Data());
3381 if(fDebug>0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
3383 if(strdir && strdir.Length()) gDirectory->cd(strdir);
3387 if(strlist && strlist.Length()){
3389 if(!(list = (TList*) gDirectory->Get(strlist))){
3390 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
3396 hn3ResponseJetPt = (THnSparse*) list->FindObject("fhnResponseJetTrackPt");
3397 hn3ResponseJetZ = (THnSparse*) list->FindObject("fhnResponseJetZ");
3398 hn3ResponseJetXi = (THnSparse*) list->FindObject("fhnResponseJetXi");
3401 hn3ResponseJetPt = (THnSparse*) gDirectory->Get("fhnResponseJetTrackPt");
3402 hn3ResponseJetZ = (THnSparse*) gDirectory->Get("fhnResponseJetZ");
3403 hn3ResponseJetXi = (THnSparse*) gDirectory->Get("fhnResponseJetXi");
3407 if(!hn3ResponseJetPt){
3408 Printf("%s:%d -- error retrieving response matrix fhnResponseJetTrackPt",(char*)__FILE__,__LINE__);
3412 if(!hn3ResponseJetZ){
3413 Printf("%s:%d -- error retrieving response matrix fhnResponseJetZ",(char*)__FILE__,__LINE__);
3417 if(!hn3ResponseJetXi){
3418 Printf("%s:%d -- error retrieving response matrix fhnResponseJetXi",(char*)__FILE__,__LINE__);
3426 Int_t axisJetPtTHn3 = -1;
3427 Int_t axisGenPtTHn3 = -1;
3428 Int_t axisRecPtTHn3 = -1;
3430 for(Int_t i=0; i<hn3ResponseJetPt->GetNdimensions(); i++){
3432 TString title = hn3ResponseJetPt->GetAxis(i)->GetTitle();
3434 if(title.Contains("jet p_{T}")) axisJetPtTHn3 = i;
3435 if(title.Contains("gen p_{T}")) axisGenPtTHn3 = i;
3436 if(title.Contains("rec p_{T}")) axisRecPtTHn3 = i;
3439 if(axisJetPtTHn3 == -1){
3440 Printf("%s:%d -- error axisJetPtTHn3",(char*)__FILE__,__LINE__);
3444 if(axisGenPtTHn3 == -1){
3445 Printf("%s:%d -- error axisGenPtTHn3",(char*)__FILE__,__LINE__);
3449 if(axisRecPtTHn3 == -1){
3450 Printf("%s:%d -- error axisRecPtTHn3",(char*)__FILE__,__LINE__);
3454 for(Int_t i=0; i<fNJetPtSlices; i++){
3456 Float_t jetPtLoLim = fJetPtSlices->At(i);
3457 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
3459 Int_t binLo = static_cast<Int_t>(hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->FindBin(jetPtLoLim));
3460 Int_t binUp = static_cast<Int_t>(hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->FindBin(jetPtUpLim))-1;
3462 if(binUp > hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->GetNbins()){
3463 Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
3467 TString strNameRespPt(Form("hnResponsePt_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3468 TString strNameRespZ(Form("hnResponseZ_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3469 TString strNameRespXi(Form("hnResponseXi_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3471 TString strNameTH2RespPt(Form("h2ResponsePt_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3472 TString strNameTH2RespZ(Form("h2ResponseZ_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3473 TString strNameTH2RespXi(Form("h2ResponseXi_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3475 TString strNameRecPt(Form("h1FFTrackPtRecPrim_recPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3476 TString strNameRecZ(Form("h1FFZRecPrim_recPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3477 TString strNameRecXi(Form("h1FFXiRecPrim_recPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3479 TString strNameGenPt(Form("h1FFTrackPtRecPrim_genPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3480 TString strNameGenZ(Form("h1FFZRecPrim_genPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3481 TString strNameGenXi(Form("h1FFXiRecPrim_genPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3483 TString strNameRatioPt(Form("h1RatioTrackPtRecPrim_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3484 TString strNameRatioZ(Form("h1RatioZRecPrim_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3485 TString strNameRatioXi(Form("h1RatioXiRecPrim_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3488 // 2D projections in jet pt range
3490 hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->SetRange(binLo,binUp);
3491 hn3ResponseJetZ->GetAxis(axisJetPtTHn3)->SetRange(binLo,binUp);
3492 hn3ResponseJetXi->GetAxis(axisJetPtTHn3)->SetRange(binLo,binUp);
3494 Int_t axesProj[2] = {axisRecPtTHn3,axisGenPtTHn3};
3496 hnResponsePt[i] = hn3ResponseJetPt->Projection(2,axesProj);
3497 hnResponseZ[i] = hn3ResponseJetZ->Projection(2,axesProj);
3498 hnResponseXi[i] = hn3ResponseJetXi->Projection(2,axesProj);
3500 hnResponsePt[i]->SetNameTitle(strNameRespPt,"");
3501 hnResponseZ[i]->SetNameTitle(strNameRespZ,"");
3502 hnResponseXi[i]->SetNameTitle(strNameRespXi,"");
3504 h2ResponsePt[i] = (TH2F*) hnResponsePt[i]->Projection(1,0);// note convention: yDim,xDim
3505 h2ResponseZ[i] = (TH2F*) hnResponseZ[i]->Projection(1,0); // note convention: yDim,xDim
3506 h2ResponseXi[i] = (TH2F*) hnResponseXi[i]->Projection(1,0);// note convention: yDim,xDim
3508 h2ResponsePt[i]->SetNameTitle(strNameTH2RespPt,"");
3509 h2ResponseZ[i]->SetNameTitle(strNameTH2RespZ,"");
3510 h2ResponseXi[i]->SetNameTitle(strNameTH2RespXi,"");
3515 Int_t axisGenPtTHn2 = -1;
3516 Int_t axisRecPtTHn2 = -1;
3518 for(Int_t d=0; d<hnResponsePt[i]->GetNdimensions(); d++){
3520 TString title = hnResponsePt[i]->GetAxis(d)->GetTitle();
3522 if(title.Contains("gen p_{T}")) axisGenPtTHn2 = d;
3523 if(title.Contains("rec p_{T}")) axisRecPtTHn2 = d;
3527 if(axisGenPtTHn2 == -1){
3528 Printf("%s:%d -- error axisGenPtTHn2",(char*)__FILE__,__LINE__);
3532 if(axisRecPtTHn2 == -1){
3533 Printf("%s:%d -- error axisRecPtTHn2",(char*)__FILE__,__LINE__);
3538 h1FFPtRec[i] = (TH1F*) hnResponsePt[i]->Projection(axisRecPtTHn2);// note convention: yDim,xDim
3539 h1FFZRec[i] = (TH1F*) hnResponseZ[i]->Projection(axisRecPtTHn2);// note convention: yDim,xDim
3540 h1FFXiRec[i] = (TH1F*) hnResponseXi[i]->Projection(axisRecPtTHn2);// note convention: yDim,xDim
3542 h1FFPtRec[i]->SetNameTitle(strNameRecPt,"");
3543 h1FFZRec[i]->SetNameTitle(strNameRecZ,"");
3544 h1FFXiRec[i]->SetNameTitle(strNameRecXi,"");
3546 h1FFPtGen[i] = (TH1F*) hnResponsePt[i]->Projection(axisGenPtTHn2);// note convention: yDim,xDim
3547 h1FFZGen[i] = (TH1F*) hnResponseZ[i]->Projection(axisGenPtTHn2);// note convention: yDim,xDim
3548 h1FFXiGen[i] = (TH1F*) hnResponseXi[i]->Projection(axisGenPtTHn2);// note convention: yDim,xDim
3550 h1FFPtGen[i]->SetNameTitle(strNameGenPt,"");
3551 h1FFZGen[i]->SetNameTitle(strNameGenZ,"");
3552 h1FFXiGen[i]->SetNameTitle(strNameGenXi,"");
3554 // normalize 1D projections
3556 if(fNHistoBinsPt[i]) h1FFPtRec[i] = (TH1F*) h1FFPtRec[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray());
3557 if(fNHistoBinsZ[i]) h1FFZRec[i] = (TH1F*) h1FFZRec[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3558 if(fNHistoBinsXi[i]) h1FFXiRec[i] = (TH1F*) h1FFXiRec[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3560 if(fNHistoBinsPt[i]) h1FFPtGen[i] = (TH1F*) h1FFPtGen[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray());
3561 if(fNHistoBinsZ[i]) h1FFZGen[i] = (TH1F*) h1FFZGen[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3562 if(fNHistoBinsXi[i]) h1FFXiGen[i] = (TH1F*) h1FFXiGen[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3564 NormalizeTH1(h1FFPtRec[i],fNJets->At(i));
3565 NormalizeTH1(h1FFZRec[i],fNJets->At(i));
3566 NormalizeTH1(h1FFXiRec[i],fNJets->At(i));
3568 NormalizeTH1(h1FFPtGen[i],fNJets->At(i));
3569 NormalizeTH1(h1FFZGen[i],fNJets->At(i));
3570 NormalizeTH1(h1FFXiGen[i],fNJets->At(i));
3572 // ratios 1D projections
3574 h1RatioPt[i] = (TH1F*) h1FFPtRec[i]->Clone(strNameRatioPt);
3575 h1RatioPt[i]->Reset();
3576 h1RatioPt[i]->Divide(h1FFPtRec[i],h1FFPtGen[i],1,1,"B");
3578 h1RatioZ[i] = (TH1F*) h1FFZRec[i]->Clone(strNameRatioZ);
3579 h1RatioZ[i]->Reset();
3580 h1RatioZ[i]->Divide(h1FFZRec[i],h1FFZGen[i],1,1,"B");
3582 h1RatioXi[i] = (TH1F*) h1FFXiRec[i]->Clone(strNameRatioXi);
3583 h1RatioXi[i]->Reset();
3584 h1RatioXi[i]->Divide(h1FFXiRec[i],h1FFXiGen[i],1,1,"B");
3590 TString outfileOption = "RECREATE";
3591 if(updateOutfile) outfileOption = "UPDATE";
3593 TFile out(strOutfile,outfileOption);
3596 Printf("%s:%d -- error opening response matrix output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
3600 if(fDebug>0) Printf("%s:%d -- write response matrices to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
3602 if(strOutDir && strOutDir.Length()){
3605 if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd();
3607 dir = out.mkdir(strOutDir);
3613 for(Int_t i=0; i<fNJetPtSlices; i++){
3615 hnResponsePt[i]->Write();
3616 hnResponseXi[i]->Write();
3617 hnResponseZ[i]->Write();
3619 h2ResponsePt[i]->Write();
3620 h2ResponseXi[i]->Write();
3621 h2ResponseZ[i]->Write();
3623 h1FFPtRec[i]->Write();
3624 h1FFZRec[i]->Write();
3625 h1FFXiRec[i]->Write();
3627 h1FFPtGen[i]->Write();
3628 h1FFZGen[i]->Write();
3629 h1FFXiGen[i]->Write();
3631 h1RatioPt[i]->Write();
3632 h1RatioZ[i]->Write();
3633 h1RatioXi[i]->Write();
3640 //______________________________________________________________________________________________________
3641 void AliFragmentationFunctionCorrections::ReadResponse(TString strfile, TString strdir, TString strlist)
3643 // read response matrices from file
3644 // argument strlist optional - read from directory strdir if not specified
3645 // note: THnSparse are not rebinned
3647 THnSparse* hRespPt[fNJetPtSlices];
3648 THnSparse* hRespZ[fNJetPtSlices];
3649 THnSparse* hRespXi[fNJetPtSlices];
3651 for(Int_t i=0; i<fNJetPtSlices; i++) hRespPt[i] = 0;
3652 for(Int_t i=0; i<fNJetPtSlices; i++) hRespZ[i] = 0;
3653 for(Int_t i=0; i<fNJetPtSlices; i++) hRespXi[i] = 0;
3655 TFile f(strfile,"READ");
3658 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
3662 if(fDebug>0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3664 if(strdir && strdir.Length()) gDirectory->cd(strdir);
3668 if(strlist && strlist.Length()){
3670 if(!(list = (TList*) gDirectory->Get(strlist))){
3671 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
3676 for(Int_t i=0; i<fNJetPtSlices; i++){
3678 Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
3679 Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
3681 TString strNameRespPt(Form("hnResponsePt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3682 TString strNameRespZ(Form("hnResponseZ_%02d_%02d",jetPtLoLim,jetPtUpLim));
3683 TString strNameRespXi(Form("hnResponseXi_%02d_%02d",jetPtLoLim,jetPtUpLim));
3686 hRespPt[i] = (THnSparse*) list->FindObject(strNameRespPt);
3687 hRespZ[i] = (THnSparse*) list->FindObject(strNameRespZ);
3688 hRespXi[i] = (THnSparse*) list->FindObject(strNameRespXi);
3691 hRespPt[i] = (THnSparse*) gDirectory->Get(strNameRespPt);
3692 hRespZ[i] = (THnSparse*) gDirectory->Get(strNameRespZ);
3693 hRespXi[i] = (THnSparse*) gDirectory->Get(strNameRespXi);
3697 Printf("%s:%d -- error retrieving response %s", (char*)__FILE__,__LINE__,strNameRespPt.Data());
3701 Printf("%s:%d -- error retrieving response %s", (char*)__FILE__,__LINE__,strNameRespZ.Data());
3705 Printf("%s:%d -- error retrieving response %s", (char*)__FILE__,__LINE__,strNameRespXi.Data());
3708 // if(0){ // can't rebin THnSparse ...
3709 // if(fNHistoBinsPt[i]) hRespPt[i]->SetBinEdges(0,fHistoBinsPt[i]->GetArray());
3710 // if(fNHistoBinsPt[i]) hRespPt[i]->SetBinEdges(1,fHistoBinsPt[i]->GetArray());
3712 // if(fNHistoBinsZ[i]) hRespZ[i]->SetBinEdges(0,fHistoBinsZ[i]->GetArray());
3713 // if(fNHistoBinsZ[i]) hRespZ[i]->SetBinEdges(1,fHistoBinsZ[i]->GetArray());
3715 // if(fNHistoBinsXi[i]) hRespXi[i]->SetBinEdges(0,fHistoBinsXi[i]->GetArray());
3716 // if(fNHistoBinsXi[i]) hRespXi[i]->SetBinEdges(1,fHistoBinsXi[i]->GetArray());
3720 } // jet slices loop
3722 f.Close(); // THnSparse pointers still valid even if file closed
3724 // for(Int_t i=0; i<fNJetPtSlices; i++){ // no copy c'tor ...
3725 // if(hRespPt[i]) new(fhnResponsePt[i]) THnSparseF(*hRespPt[i]);
3726 // if(hRespZ[i]) new(fhnResponseZ[i]) THnSparseF(*hRespZ[i]);
3727 // if(hRespXi[i]) new(fhnResponseXi[i]) THnSparseF(*hRespXi[i]);
3730 for(Int_t i=0; i<fNJetPtSlices; i++){
3731 fhnResponsePt[i] = hRespPt[i];
3732 fhnResponseZ[i] = hRespZ[i];
3733 fhnResponseXi[i] = hRespXi[i];
3737 //______________________________________________________________________________________________________________________
3738 void AliFragmentationFunctionCorrections::ReadPriors(TString strfile,const Int_t type)
3740 // read priors from file: rec primaries, gen pt dist
3742 if(fDebug>0) Printf("%s:%d -- read priors from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3744 // temporary histos to store pointers from file
3745 TH1F* hist[fNJetPtSlices];
3747 for(Int_t i=0; i<fNJetPtSlices; i++) hist[i] = 0;
3749 TFile f(strfile,"READ");
3752 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
3756 for(Int_t i=0; i<fNJetPtSlices; i++){
3758 Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
3759 Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
3763 if(type == kFlagPt) strName.Form("h1FFTrackPtRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim);
3764 if(type == kFlagZ) strName.Form("h1FFZRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim);
3765 if(type == kFlagXi) strName.Form("h1FFXiRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim);
3767 hist[i] = (TH1F*) gDirectory->Get(strName);
3770 Printf("%s:%d -- error retrieving prior %s", (char*)__FILE__,__LINE__,strName.Data());
3774 //if(fNHistoBinsPt[i]) hist[i] = (TH1F*) hist[i]->Rebin(fNHistoBinsPt[i],hist[i]->GetName()+"_rebin",fHistoBinsPt[i]->GetArray());
3776 if(hist[i]) hist[i]->SetDirectory(0);
3778 } // jet slices loop
3783 for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
3784 if(hist[i] && type == kFlagPt) new(fh1FFTrackPtPrior[i]) TH1F(*hist[i]);
3785 if(hist[i] && type == kFlagZ) new(fh1FFZPrior[i]) TH1F(*hist[i]);
3786 if(hist[i] && type == kFlagXi) new(fh1FFXiPrior[i]) TH1F(*hist[i]);
3791 //_____________________________________________________
3792 // void AliFragmentationFunctionCorrections::RatioRecGen()
3794 // // create ratio reconstructed over generated FF
3795 // // use current highest corrLevel
3797 // Printf("%s:%d -- build ratio rec.gen, corrLevel %d",(char*)__FILE__,__LINE__,fNCorrectionLevels-1);
3799 // for(Int_t i=0; i<fNJetPtSlices; i++){
3801 // TH1F* histPtRec = fCorrFF[fNCorrectionLevels-1]->GetTrackPt(i); // levels -1: latest corr level
3802 // TH1F* histZRec = fCorrFF[fNCorrectionLevels-1]->GetZ(i); // levels -1: latest corr level
3803 // TH1F* histXiRec = fCorrFF[fNCorrectionLevels-1]->GetXi(i); // levels -1: latest corr level
3805 // TH1F* histPtGen = fh1FFTrackPtGenPrim[i];
3806 // TH1F* histZGen = fh1FFZGenPrim[i];
3807 // TH1F* histXiGen = fh1FFXiGenPrim[i];
3809 // TString histNamePt = histPtRec->GetName();
3810 // TString histNameZ = histZRec->GetName();
3811 // TString histNameXi = histXiRec->GetName();
3813 // histNamePt.ReplaceAll("fh1FF","fh1FFRatioRecGen");
3814 // histNameZ.ReplaceAll("fh1FF","fh1FFRatioRecGen");
3815 // histNameXi.ReplaceAll("fh1FF","fh1FFRatioRecGen");
3818 // TH1F* hRatioRecGenPt = (TH1F*) histPtRec->Clone(histNamePt);
3819 // hRatioRecGenPt->Reset();
3820 // hRatioRecGenPt->Divide(histPtRec,histPtGen,1,1,"B");
3822 // TH1F* hRatioRecGenZ = (TH1F*) histZRec->Clone(histNameZ);
3823 // hRatioRecGenZ->Reset();
3824 // hRatioRecGenZ->Divide(histZRec,histZGen,1,1,"B");
3826 // TH1F* hRatioRecGenXi = (TH1F*) histXiRec->Clone(histNameXi);
3827 // hRatioRecGenXi->Reset();
3828 // hRatioRecGenXi->Divide(histXiRec,histXiGen,1,1,"B");
3830 // new(fh1FFRatioRecGenPt[i]) TH1F(*hRatioRecGenPt);
3831 // new(fh1FFRatioRecGenZ[i]) TH1F(*hRatioRecGenZ);
3832 // new(fh1FFRatioRecGenXi[i]) TH1F(*hRatioRecGenXi);
3836 // //___________________________________________________________
3837 // void AliFragmentationFunctionCorrections::RatioRecPrimaries()
3839 // // create ratio reconstructed tracks over reconstructed primaries
3840 // // use raw FF (corrLevel 0)
3842 // Printf("%s:%d -- build ratio rec tracks /rec primaries",(char*)__FILE__,__LINE__);
3844 // for(Int_t i=0; i<fNJetPtSlices; i++){
3846 // const Int_t corrLevel = 0;
3848 // TH1F* histPtRec = fCorrFF[corrLevel]->GetTrackPt(i); // levels -1: latest corr level
3849 // TH1F* histZRec = fCorrFF[corrLevel]->GetZ(i); // levels -1: latest corr level
3850 // TH1F* histXiRec = fCorrFF[corrLevel]->GetXi(i); // levels -1: latest corr level
3852 // TH1F* histPtRecPrim = fh1FFTrackPtRecPrim[i];
3853 // TH1F* histZRecPrim = fh1FFZRecPrim[i];
3854 // TH1F* histXiRecPrim = fh1FFXiRecPrim[i];
3856 // TString histNamePt = histPtRec->GetName();
3857 // TString histNameZ = histZRec->GetName();
3858 // TString histNameXi = histXiRec->GetName();
3860 // histNamePt.ReplaceAll("fh1FF","fh1FFRatioRecPrim");
3861 // histNameZ.ReplaceAll("fh1FF","fh1FFRatioRecPrim");
3862 // histNameXi.ReplaceAll("fh1FF","fh1FFRatioRecPrim");
3865 // TH1F* hRatioRecPrimPt = (TH1F*) histPtRec->Clone(histNamePt);
3866 // hRatioRecPrimPt->Reset();
3867 // hRatioRecPrimPt->Divide(histPtRec,histPtRecPrim,1,1,"B");
3869 // TH1F* hRatioRecPrimZ = (TH1F*) histZRec->Clone(histNameZ);
3870 // hRatioRecPrimZ->Reset();
3871 // hRatioRecPrimZ->Divide(histZRec,histZRecPrim,1,1,"B");
3873 // TH1F* hRatioRecPrimXi = (TH1F*) histXiRec->Clone(histNameXi);
3874 // hRatioRecPrimXi->Reset();
3875 // hRatioRecPrimXi->Divide(histXiRec,histXiRecPrim,1,1,"B");
3878 // new(fh1FFRatioRecPrimPt[i]) TH1F(*hRatioRecPrimPt);
3879 // new(fh1FFRatioRecPrimZ[i]) TH1F(*hRatioRecPrimZ);
3880 // new(fh1FFRatioRecPrimXi[i]) TH1F(*hRatioRecPrimXi);
3884 // __________________________________________________________________________________
3885 void AliFragmentationFunctionCorrections::ProjectJetResponseMatrices(TString strOutfile)
3888 // project response matrices on both axes:
3889 // FF for rec primaries, in terms of generated and reconstructed momentum
3890 // write FF and ratios to outFile
3892 Printf("%s:%d -- project response matrices, write to %s",(char*)__FILE__,__LINE__,strOutfile.Data());
3894 TH1F* hFFPtRec[fNJetPtSlices];
3895 TH1F* hFFZRec[fNJetPtSlices];
3896 TH1F* hFFXiRec[fNJetPtSlices];
3898 TH1F* hFFPtGen[fNJetPtSlices];
3899 TH1F* hFFZGen[fNJetPtSlices];
3900 TH1F* hFFXiGen[fNJetPtSlices];
3902 TH1F* hRatioPt[fNJetPtSlices];
3903 TH1F* hRatioZ[fNJetPtSlices];
3904 TH1F* hRatioXi[fNJetPtSlices];
3907 Int_t axisGenPt = 1;
3908 Int_t axisRecPt = 0;
3910 for(Int_t i=0; i<fNJetPtSlices; i++){
3912 Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
3913 Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
3915 TString strNameRecPt(Form("h1FFTrackPtRecPrim_recPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3916 TString strNameRecZ(Form("h1FFZRecPrim_recPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3917 TString strNameRecXi(Form("h1FFXiRecPrim_recPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3919 TString strNameGenPt(Form("h1FFTrackPtRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3920 TString strNameGenZ(Form("h1FFZRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3921 TString strNameGenXi(Form("h1FFXiRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3923 TString strNameRatioPt(Form("h1RatioTrackPtRecPrim_%02d_%02d",jetPtLoLim,jetPtUpLim));
3924 TString strNameRatioZ(Form("h1RatioZRecPrim_%02d_%02d",jetPtLoLim,jetPtUpLim));
3925 TString strNameRatioXi(Form("h1RatioXiRecPrim_%02d_%02d",jetPtLoLim,jetPtUpLim));
3928 hFFPtRec[i] = (TH1F*) fhnResponsePt[i]->Projection(axisRecPt);// note convention: yDim,xDim
3929 hFFZRec[i] = (TH1F*) fhnResponseZ[i]->Projection(axisRecPt);// note convention: yDim,xDim
3930 hFFXiRec[i] = (TH1F*) fhnResponseXi[i]->Projection(axisRecPt);// note convention: yDim,xDim
3932 hFFPtRec[i]->SetNameTitle(strNameRecPt,"");
3933 hFFZRec[i]->SetNameTitle(strNameRecZ,"");
3934 hFFXiRec[i]->SetNameTitle(strNameRecXi,"");
3937 hFFPtGen[i] = (TH1F*) fhnResponsePt[i]->Projection(axisGenPt);// note convention: yDim,xDim
3938 hFFZGen[i] = (TH1F*) fhnResponseZ[i]->Projection(axisGenPt);// note convention: yDim,xDim
3939 hFFXiGen[i] = (TH1F*) fhnResponseXi[i]->Projection(axisGenPt);// note convention: yDim,xDim
3941 hFFPtGen[i]->SetNameTitle(strNameGenPt,"");
3942 hFFZGen[i]->SetNameTitle(strNameGenZ,"");
3943 hFFXiGen[i]->SetNameTitle(strNameGenXi,"");
3946 if(fNHistoBinsPt[i]) hFFPtRec[i] = (TH1F*) hFFPtRec[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray());
3947 if(fNHistoBinsZ[i]) hFFZRec[i] = (TH1F*) hFFZRec[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3948 if(fNHistoBinsXi[i]) hFFXiRec[i] = (TH1F*) hFFXiRec[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3950 if(fNHistoBinsPt[i]) hFFPtGen[i] = (TH1F*) hFFPtGen[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray());
3951 if(fNHistoBinsZ[i]) hFFZGen[i] = (TH1F*) hFFZGen[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3952 if(fNHistoBinsXi[i]) hFFXiGen[i] = (TH1F*) hFFXiGen[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3954 NormalizeTH1(hFFPtGen[i],fNJets->At(i));
3955 NormalizeTH1(hFFZGen[i],fNJets->At(i));
3956 NormalizeTH1(hFFXiGen[i],fNJets->At(i));
3958 NormalizeTH1(hFFPtRec[i],fNJets->At(i));
3959 NormalizeTH1(hFFZRec[i],fNJets->At(i));
3960 NormalizeTH1(hFFXiRec[i],fNJets->At(i));
3963 hRatioPt[i] = (TH1F*) hFFPtRec[i]->Clone(strNameRatioPt);
3964 hRatioPt[i]->Reset();
3965 hRatioPt[i]->Divide(hFFPtRec[i],hFFPtGen[i],1,1,"B");
3967 hRatioZ[i] = (TH1F*) hFFZRec[i]->Clone(strNameRatioZ);
3968 hRatioZ[i]->Reset();
3969 hRatioZ[i]->Divide(hFFZRec[i],hFFZGen[i],1,1,"B");
3971 hRatioXi[i] = (TH1F*) hFFXiRec[i]->Clone(strNameRatioXi);
3972 hRatioXi[i]->Reset();
3973 hRatioXi[i]->Divide(hFFXiRec[i],hFFXiGen[i],1,1,"B");
3980 TFile out(strOutfile,"RECREATE");
3983 Printf("%s:%d -- error opening response matrix output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
3987 for(Int_t i=0; i<fNJetPtSlices; i++){
3989 hFFPtRec[i]->Write();
3990 hFFZRec[i]->Write();
3991 hFFXiRec[i]->Write();
3993 hFFPtGen[i]->Write();
3994 hFFZGen[i]->Write();
3995 hFFXiGen[i]->Write();
3997 hRatioPt[i]->Write();
3998 hRatioZ[i]->Write();
3999 hRatioXi[i]->Write();
4005 // ____________________________________________________________________________________________________________________________
4006 void AliFragmentationFunctionCorrections::ProjectSingleResponseMatrix(TString strOutfile, Bool_t updateOutfile, TString strOutDir)
4008 // project response matrix on both axes:
4009 // pt spec for rec primaries, in terms of generated and reconstructed momentum
4010 // write spec and ratios to outFile
4012 Printf("%s:%d -- project single pt response matrix, write to %s",(char*)__FILE__,__LINE__,strOutfile.Data());
4018 Int_t axisGenPt = 1;
4019 Int_t axisRecPt = 0;
4021 TString strNameRecPt = "h1SpecTrackPtRecPrim_recPt";
4022 TString strNameGenPt = "h1SpecTrackPtRecPrim_genPt";
4023 TString strNameRatioPt = "h1RatioTrackPtRecPrim";
4025 hSpecPtRec = (TH1F*) fhnResponseSinglePt->Projection(axisRecPt);// note convention: yDim,xDim
4026 hSpecPtRec->SetNameTitle(strNameRecPt,"");
4028 hSpecPtGen = (TH1F*) fhnResponseSinglePt->Projection(axisGenPt);// note convention: yDim,xDim
4029 hSpecPtGen->SetNameTitle(strNameGenPt,"");
4031 hRatioPt = (TH1F*) hSpecPtRec->Clone(strNameRatioPt);
4033 hRatioPt->Divide(hSpecPtRec,hSpecPtGen,1,1,"B");
4035 TString outfileOption = "RECREATE";
4036 if(updateOutfile) outfileOption = "UPDATE";
4038 TFile out(strOutfile,outfileOption);
4041 Printf("%s:%d -- error opening reponse matrix projections output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
4046 if(strOutDir && strOutDir.Length()){
4049 if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd();
4051 dir = out.mkdir(strOutDir);
4056 hSpecPtRec->Write();
4057 hSpecPtGen->Write();
4064 //__________________________________________________________________________________________________________________________________________________________________
4065 void AliFragmentationFunctionCorrections::RebinHisto(const Int_t jetPtSlice, const Int_t nBinsLimits, Double_t* binsLimits, Double_t* binsWidth, const Int_t type)
4067 // rebin histo, rescale bins according to new width
4068 // only correct for input histos with equal bin size
4070 // args: jetPtSlice, type, use current corr level
4072 // function takes array of limits and widths (e.g. 1 GeV bins up to 10 GeV, 2 GeV width up to 20 GeV, ...)
4073 // array size of binsLimits: nBinsLimits
4074 // array size of binsWidth: nBinsLimits-1
4075 // binsLimits have to be in increasing order
4076 // if binning undefined for any slice, original binning will be kept
4079 Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
4083 if(jetPtSlice>=fNJetPtSlices){
4084 Printf("%s:%d -- jetPtSlice %d exceeds max",(char*)__FILE__,__LINE__,jetPtSlice);
4089 Double_t binLimitMin = binsLimits[0];
4090 Double_t binLimitMax = binsLimits[nBinsLimits-1];
4092 Double_t binLimit = binLimitMin; // start value
4094 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 ...
4095 TArrayD binsArray(sizeUpperLim);
4097 binsArray.SetAt(binLimitMin,nBins++);
4099 while(binLimit<binLimitMax && nBins<sizeUpperLim){
4101 Int_t currentSlice = -1;
4102 for(Int_t i=0; i<nBinsLimits; i++){
4103 if(binLimit >= 0.999*binsLimits[i]) currentSlice = i; // 0.999 numerical saftey factor
4106 Double_t currentBinWidth = binsWidth[currentSlice];
4107 binLimit += currentBinWidth;
4109 binsArray.SetAt(binLimit,nBins++);
4114 if(type == kFlagPt) hist = fCorrFF[fNCorrectionLevels-1]->GetTrackPt(jetPtSlice);
4115 else if(type == kFlagZ) hist = fCorrFF[fNCorrectionLevels-1]->GetZ(jetPtSlice);
4116 else if(type == kFlagXi) hist = fCorrFF[fNCorrectionLevels-1]->GetXi(jetPtSlice);
4117 else if(type == kFlagSinglePt) hist = fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->GetTrackPt(0);
4119 Printf("%s%d unknown type",(char*)__FILE__,__LINE__);
4123 Double_t binWidthNoRebin = hist->GetBinWidth(1);
4125 Double_t* bins = binsArray.GetArray();
4127 hist = (TH1F*) hist->Rebin(nBins-1,"",bins);
4129 for(Int_t bin=0; bin <= hist->GetNbinsX(); bin++){
4131 Double_t binWidthRebin = hist->GetBinWidth(bin);
4132 Double_t scaleF = binWidthNoRebin / binWidthRebin;
4134 Double_t binCont = hist->GetBinContent(bin);
4135 Double_t binErr = hist->GetBinError(bin);
4140 hist->SetBinContent(bin,binCont);
4141 hist->SetBinError(bin,binErr);
4146 TH1F* temp = new TH1F(*hist);
4148 if(type == kFlagPt) fCorrFF[fNCorrectionLevels-1]->ReplaceCorrHistos(jetPtSlice,temp,0,0);
4149 if(type == kFlagZ) fCorrFF[fNCorrectionLevels-1]->ReplaceCorrHistos(jetPtSlice,0,temp,0);
4150 if(type == kFlagXi) fCorrFF[fNCorrectionLevels-1]->ReplaceCorrHistos(jetPtSlice,0,0,temp);
4151 if(type == kFlagSinglePt) fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->ReplaceCorrHistos(0,temp,0,0);
4156 //__________________________________________________________________________________________________________________________________________________________________
4157 void AliFragmentationFunctionCorrections::WriteJetSpecResponse(TString strInfile, TString strdir, TString strlist/*, TString strOutfile*/)
4160 if(fDebug>0) Printf("%s:%d -- read jet spectrum response matrix from file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
4162 if(strdir && strdir.Length()) gDirectory->cd(strdir);
4166 if(strlist && strlist.Length()){
4167 if(!(list = (TList*) gDirectory->Get(strlist))){
4168 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
4172 if(list == 0)return; // catch strlist.Lenght() == 0;
4174 THnSparse* hn6ResponseJetPt = (THnSparse*) list->FindObject("fhnCorrelation");
4176 Int_t axis6RecJetPt = 0;
4177 Int_t axis6GenJetPt = 3;
4179 hn6ResponseJetPt->GetAxis(axis6RecJetPt)->SetTitle("rec jet p_{T} (GeV/c)");
4180 hn6ResponseJetPt->GetAxis(axis6GenJetPt)->SetTitle("gen jet p_{T} (GeV/c)");
4182 Int_t nBinsRecPt = hn6ResponseJetPt->GetAxis(axis6RecJetPt)->GetNbins();
4183 Double_t loLimRecPt = hn6ResponseJetPt->GetAxis(axis6RecJetPt)->GetBinLowEdge(1);
4184 Double_t upLimRecPt = hn6ResponseJetPt->GetAxis(axis6RecJetPt)->GetBinUpEdge(nBinsRecPt);
4186 Int_t nBinsGenPt = hn6ResponseJetPt->GetAxis(axis6GenJetPt)->GetNbins();
4187 Double_t loLimGenPt = hn6ResponseJetPt->GetAxis(axis6GenJetPt)->GetBinLowEdge(1);
4188 Double_t upLimGenPt = hn6ResponseJetPt->GetAxis(axis6GenJetPt)->GetBinUpEdge(nBinsGenPt);
4190 Int_t nBinsTrackPt = 200;
4191 Int_t loLimTrackPt = 0;
4192 Int_t upLimTrackPt = 200;
4195 Int_t nBinsResponse[4] = {nBinsRecPt,nBinsTrackPt,nBinsGenPt,nBinsTrackPt};
4196 Double_t binMinResponse[4] = {loLimRecPt,loLimTrackPt,loLimGenPt,loLimTrackPt};
4197 Double_t binMaxResponse[4] = {upLimRecPt,upLimTrackPt,upLimGenPt,upLimTrackPt};
4199 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)"};
4201 THnSparseD* hn4ResponseTrackPtJetPt = new THnSparseD("hn4ResponseTrackPtJetPt","",4,nBinsResponse,binMinResponse,binMaxResponse);
4203 for(Int_t i=0; i<4; i++){
4204 hn4ResponseTrackPtJetPt->GetAxis(i)->SetTitle(labelsResponseSinglePt[i]);
4213 //_____________________________________________________________________________________________________________________________________
4214 void AliFragmentationFunctionCorrections::ReadSingleTrackEfficiency(TString strfile, TString strdir, TString strlist, TString strname)
4217 ReadSingleTrackCorrection(strfile,strdir,strlist,strname,kFlagEfficiency);
4221 //_____________________________________________________________________________________________________________________________________
4222 void AliFragmentationFunctionCorrections::ReadSingleTrackResponse(TString strfile, TString strdir, TString strlist, TString strname)
4225 ReadSingleTrackCorrection(strfile,strdir,strlist,strname,kFlagResponse);
4229 //_____________________________________________________________________________________________________________________________________
4230 void AliFragmentationFunctionCorrections::ReadSingleTrackSecCorr(TString strfile, TString strdir, TString strlist, TString strname)
4233 ReadSingleTrackCorrection(strfile,strdir,strlist,strname,kFlagSecondaries);
4237 //______________________________________________________________________________________________________________________________________________________
4238 void AliFragmentationFunctionCorrections::ReadSingleTrackCorrection(TString strfile, TString strdir, TString strlist, TString strname, const Int_t type)
4240 // read single track correction (pt) from file
4241 // type: efficiency / response / secondaries correction
4243 if(!((type == kFlagEfficiency) || (type == kFlagResponse) || (type == kFlagSecondaries))){
4244 Printf("%s:%d -- no such correction ",(char*)__FILE__,__LINE__);
4248 TFile f(strfile,"READ");
4251 Printf("%s:%d -- error opening file %s", (char*)__FILE__,__LINE__,strfile.Data());
4255 if(fDebug>0 && type==kFlagEfficiency) Printf("%s:%d -- read single track corr from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
4256 if(fDebug>0 && type==kFlagResponse) Printf("%s:%d -- read single track response from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
4257 if(fDebug>0 && type==kFlagSecondaries) Printf("%s:%d -- read single track secondaries corr from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
4259 if(strdir && strdir.Length()) gDirectory->cd(strdir);
4263 if(strlist && strlist.Length()){
4265 if(!(list = (TList*) gDirectory->Get(strlist))){
4266 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
4271 TH1F* h1CorrHist = 0; // common TObject pointer not possible, need SetDirectory() later
4272 THnSparse* hnCorrHist = 0;
4274 if(type == kFlagEfficiency || type == kFlagSecondaries){
4276 if(list) h1CorrHist = (TH1F*) list->FindObject(strname);
4277 else h1CorrHist = (TH1F*) gDirectory->Get(strname);
4280 Printf("%s:%d -- error retrieving histo %s", (char*)__FILE__,__LINE__,strname.Data());
4285 else if(type == kFlagResponse){
4287 if(list) hnCorrHist = (THnSparse*) list->FindObject(strname);
4288 else hnCorrHist = (THnSparse*) gDirectory->Get(strname);
4291 Printf("%s:%d -- error retrieving histo %s", (char*)__FILE__,__LINE__,strname.Data());
4297 if(h1CorrHist) h1CorrHist->SetDirectory(0);
4298 //if(hnCorrHist) hnCorrHist->SetDirectory(0);
4302 if(type == kFlagEfficiency) fh1EffSinglePt = h1CorrHist;
4303 else if(type == kFlagResponse) fhnResponseSinglePt = hnCorrHist;
4304 else if(type == kFlagSecondaries) fh1SecCorrSinglePt = h1CorrHist;
4308 //________________________________________________________________________________________________________________
4309 void AliFragmentationFunctionCorrections::ReadRawPtSpec(TString strInfile, TString strID)
4311 // read track pt spec from task ouput - standard dir/list
4313 TString strdir = "PWGJE_FragmentationFunction_" + strID;
4314 TString strlist = "fracfunc_" + strID;
4316 ReadRawPtSpec(strInfile,strdir,strlist);
4319 //_______________________________________________________________________________________________________
4320 void AliFragmentationFunctionCorrections::ReadRawPtSpec(TString strfile, TString strdir, TString strlist)
4322 // get raw pt spectra from file
4325 fNCorrectionLevelsSinglePt = 0;
4326 fCorrSinglePt = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels];
4327 AddCorrectionLevelSinglePt(); // first 'correction' level = raw spectrum
4329 // get raw pt spec from input file, normalize
4331 TFile f(strfile,"READ");
4334 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
4338 if(fDebug>0) Printf("%s:%d -- read raw spectra from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
4340 gDirectory->cd(strdir);
4344 if(!(list = (TList*) gDirectory->Get(strlist))){
4345 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
4349 TString hnameTrackPt("fh1TrackQAPtRecCuts");
4350 TString hnameEvtSel("fh1EvtSelection");
4352 TH1F* fh1TrackPt = (TH1F*) list->FindObject(hnameTrackPt);
4353 TH1F* fh1EvtSel = (TH1F*) list->FindObject(hnameEvtSel);
4355 if(!fh1TrackPt){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameTrackPt.Data()); return; }
4356 if(!fh1EvtSel) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameEvtSel.Data()); return; }
4358 Float_t nEvents = fh1EvtSel->GetBinContent(fh1EvtSel->FindBin(0));
4360 fh1TrackPt->SetDirectory(0);
4364 // rebin + normalize
4365 if(fNHistoBinsSinglePt) fh1TrackPt = (TH1F*) fh1TrackPt->Rebin(fNHistoBinsSinglePt,"",fHistoBinsSinglePt->GetArray());
4367 NormalizeTH1(fh1TrackPt,nEvents);
4369 // raw FF = corr level 0
4370 fCorrSinglePt[0]->AddCorrHistos(0,fh1TrackPt);
4374 //_______________________________________________________________________________________________________
4375 void AliFragmentationFunctionCorrections::ReadRawPtSpecQATask(TString strfile, TString strdir, TString strlist)
4377 // get raw pt spectra from file
4378 // for output from Martas QA task
4382 fNCorrectionLevelsSinglePt = 0;
4383 fCorrSinglePt = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels];
4384 AddCorrectionLevelSinglePt(); // first 'correction' level = raw spectrum
4386 // get raw pt spec from input file, normalize
4388 TFile f(strfile,"READ");
4391 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
4395 if(fDebug>0) Printf("%s:%d -- read raw pt spec from QA task output file %s ",(char*)__FILE__,__LINE__,strfile.Data());
4397 gDirectory->cd(strdir);
4401 if(!(list = (TList*) gDirectory->Get(strlist))){
4402 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
4406 TString hnameTrackPt("fPtSel");
4407 TString hnameEvtSel("fNEventAll");
4409 TH1F* fh1TrackPt = (TH1F*) list->FindObject(hnameTrackPt);
4410 TH1F* fh1EvtSel = (TH1F*) list->FindObject(hnameEvtSel);
4412 if(!fh1TrackPt){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameTrackPt.Data()); return; }
4413 if(!fh1EvtSel) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameEvtSel.Data()); return; }
4416 // evts after physics selection
4417 Float_t nEvents = fh1EvtSel->GetEntries();
4419 fh1TrackPt->SetDirectory(0);
4424 NormalizeTH1(fh1TrackPt,nEvents);
4426 // raw FF = corr level 0
4427 fCorrSinglePt[0]->AddCorrHistos(0,fh1TrackPt);
4430 // ________________________________________________________
4431 void AliFragmentationFunctionCorrections::EffCorrSinglePt()
4433 // apply efficiency correction to inclusive track pt spec
4435 AddCorrectionLevelSinglePt("eff");
4437 TH1F* histPt = fCorrSinglePt[fNCorrectionLevelsSinglePt-2]->GetTrackPt(0);
4439 if(histPt->GetNbinsX() != fh1EffSinglePt->GetNbinsX()) Printf("%s:%d: inconsistency pt spec and eff corr bins - rebin effCorr ...", (char*)__FILE__,__LINE__);
4441 TString histNamePt = histPt->GetName();
4442 TH1F* hTrackPtEffCorr = (TH1F*) histPt->Clone(histNamePt);
4443 hTrackPtEffCorr->Divide(histPt,fh1EffSinglePt,1,1,"");
4445 fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->AddCorrHistos(0,hTrackPtEffCorr);
4448 //___________________________________________________________________________________________________________________________
4449 void AliFragmentationFunctionCorrections::UnfoldSinglePt(const Int_t nIter, const Bool_t useCorrelatedErrors)
4451 // unfolde inclusive dN/dpt spectra
4453 AddCorrectionLevelSinglePt("unfold");
4455 TH1F* hist = fCorrSinglePt[fNCorrectionLevelsSinglePt-2]->GetTrackPt(0); // level -2: before unfolding, level -1: unfolded
4456 THnSparse* hnResponse = fhnResponseSinglePt;
4458 TString histNameTHn = hist->GetName();
4459 if(histNameTHn.Contains("TH1")) histNameTHn.ReplaceAll("TH1","THn");
4460 if(histNameTHn.Contains("fPt")) histNameTHn.ReplaceAll("fPt","fhnPt");
4463 TString histNameBackFolded = hist->GetName();
4464 histNameBackFolded.Append("_backfold");
4466 TString histNameRatioFolded = hist->GetName();
4467 if(histNameRatioFolded.Contains("fh1")) histNameRatioFolded.ReplaceAll("fh1","hRatio");
4468 if(histNameRatioFolded.Contains("fPt")) histNameRatioFolded.ReplaceAll("fPt","hRatioPt");
4469 histNameRatioFolded.Append("_unfold");
4471 TString histNameRatioBackFolded = hist->GetName();
4472 if(histNameRatioBackFolded.Contains("fh1")) histNameRatioBackFolded.ReplaceAll("fh1","hRatio");
4473 if(histNameRatioBackFolded.Contains("fPt")) histNameRatioBackFolded.ReplaceAll("fPt","hRatioPt");
4474 histNameRatioBackFolded.Append("_backfold");
4476 THnSparse* hnHist = TH1toSparse(hist,histNameTHn,hist->GetTitle());
4477 THnSparse* hnFlatEfficiency = TH1toSparse(hist,"fhnEfficiency","eff",kTRUE); // could optionally also use real eff
4480 = Unfold(hnHist,hnResponse,hnFlatEfficiency,nIter,useCorrelatedErrors);
4482 //TH1F* hUnfolded = (TH1F*) hnUnfolded->Projection(0);
4483 //hUnfolded->SetNameTitle(hist->GetName(),hist->GetTitle());
4485 fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->AddCorrHistos(0,hUnfolded);
4487 // backfolding: apply response matrix to unfolded spectrum
4488 TH1F* hBackFolded = ApplyResponse(hUnfolded,hnResponse);
4489 hBackFolded->SetNameTitle(histNameBackFolded,hist->GetTitle());
4491 fh1SingleTrackPtBackFolded = hBackFolded;
4494 // ratio unfolded to original histo
4495 TH1F* hRatioUnfolded = (TH1F*) hUnfolded->Clone(histNameRatioFolded);
4496 hRatioUnfolded->Reset();
4497 hRatioUnfolded->Divide(hUnfolded,hist,1,1,"B");
4499 fh1RatioSingleTrackPtFolded = hRatioUnfolded;
4502 // ratio backfolded to original histo
4503 TH1F* hRatioBackFolded = (TH1F*) hBackFolded->Clone(histNameRatioBackFolded);
4504 hRatioBackFolded->Reset();
4505 hRatioBackFolded->Divide(hBackFolded,hist,1,1,"B");
4507 fh1RatioSingleTrackPtBackFolded = hRatioBackFolded;
4510 delete hnFlatEfficiency;
4514 // ________________________________________________________
4515 void AliFragmentationFunctionCorrections::SecCorrSinglePt()
4517 // apply efficiency correction to inclusive track pt spec
4519 AddCorrectionLevelSinglePt("secCorr");
4521 TH1F* histPt = fCorrSinglePt[fNCorrectionLevelsSinglePt-2]->GetTrackPt(0);
4523 if(histPt->GetNbinsX() != fh1SecCorrSinglePt->GetNbinsX())
4524 Printf("%s:%d: inconsistency pt spec and secondaries corr bins - rebin effCorr ...", (char*)__FILE__,__LINE__);
4526 TString histNamePt = histPt->GetName();
4527 TH1F* hTrackPtSecCorr = (TH1F*) histPt->Clone(histNamePt);
4529 hTrackPtSecCorr->Multiply(histPt,fh1SecCorrSinglePt,1,1,"");
4531 fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->AddCorrHistos(0,hTrackPtSecCorr);
4534 //___________________________________________________________________________________________________________________________________
4535 void AliFragmentationFunctionCorrections::dNdz2dNdxi()
4537 // transform dN/dz distribution into dN/dxi
4538 // for current corr level, all jet pt slices
4541 for(Int_t i=0; i<fNJetPtSlices; i++){
4543 TH1F* histZ = fCorrFF[fNCorrectionLevels-1]->GetZ(i);
4544 Int_t nBins = histZ->GetNbinsX();
4546 Double_t* binLims = new Double_t[nBins+1];
4548 for(Int_t bin = 0; bin<nBins; bin++){
4550 Int_t binZ = nBins-bin;
4552 Double_t zLo = histZ->GetXaxis()->GetBinLowEdge(binZ);
4553 Double_t zUp = histZ->GetXaxis()->GetBinUpEdge(binZ);
4555 Double_t xiLo = TMath::Log(1/zUp);
4556 Double_t xiUp = TMath::Log(1/zLo);
4558 if(bin == 0) binLims[0] = xiLo;
4559 binLims[bin+1] = xiUp;
4562 // for(Int_t bin = 0; bin<=nBins; bin++) std::cout<<" bin "<<bin<<" binLims "<<binLims[bin]<<std::endl;
4564 TString strTitle = histZ->GetTitle();
4565 TString strName = histZ->GetName();
4567 strName.ReplaceAll("Z","XiNew");
4568 strTitle.ReplaceAll("Z","XiNew");
4571 TH1F* histXiNew = new TH1F("histXiNew","",nBins,binLims);
4572 histXiNew->SetNameTitle(strName,strTitle);
4575 for(Int_t binZ = 1; binZ<=nBins; binZ++){
4577 Double_t meanZ = histZ->GetBinCenter(binZ);
4578 Double_t cont = histZ->GetBinContent(binZ);
4579 Double_t err = histZ->GetBinError(binZ);
4581 Double_t meanXi = TMath::Log(1/meanZ);
4582 Int_t binXi = histXiNew->FindBin(meanXi);
4584 histXiNew->SetBinContent(binXi,cont);
4585 histXiNew->SetBinError(binXi,err);
4587 //std::cout<<" binZ "<<binZ<<" meanZ "<<meanZ<<" binXi "<<binXi<<" meanXi "<<meanXi<<" cont "<<cont<<" err "<<err<<std::endl;
4590 fCorrFF[fNCorrectionLevels-1]->ReplaceCorrHistos(i,0,0,histXiNew);
4597 //________________________________________________________________________________________________________________
4598 void AliFragmentationFunctionCorrections::WriteBinShiftCorr(TString strInfile, TString strIDGen, TString strIDRec,
4599 TString strOutfile, Bool_t updateOutfile, Bool_t useRecPrim,
4602 TString strdirGen = "PWGJE_FragmentationFunction_" + strIDGen;
4603 TString strlistGen = "fracfunc_" + strIDGen;
4605 TString strdirRec = "PWGJE_FragmentationFunction_" + strIDRec;
4606 TString strlistRec = "fracfunc_" + strIDRec;
4608 WriteBinShiftCorr(strInfile,strdirGen,strlistGen,strdirRec,strlistRec,strOutfile,updateOutfile,useRecPrim, strOutDir, kFALSE, "");
4611 //________________________________________________________________________________________________________________
4612 void AliFragmentationFunctionCorrections::WriteBgrBinShiftCorr(TString strInfile, TString strBgrID, TString strIDGen, TString strIDRec,
4613 TString strOutfile, Bool_t updateOutfile, Bool_t useRecPrim,
4616 TString strdirGen = "PWGJE_FragmentationFunction_" + strIDGen;
4617 TString strlistGen = "fracfunc_" + strIDGen;
4619 TString strdirRec = "PWGJE_FragmentationFunction_" + strIDRec;
4620 TString strlistRec = "fracfunc_" + strIDRec;
4622 WriteBinShiftCorr(strInfile,strdirGen,strlistGen,strdirRec,strlistRec,strOutfile,updateOutfile,useRecPrim, strOutDir, kTRUE, strBgrID);
4625 //___________________________________________________________________________________________________________________________________
4626 void AliFragmentationFunctionCorrections::WriteBinShiftCorr(TString strInfile, TString strdirGen, TString strlistGen,
4627 TString strdirRec, TString strlistRec,
4628 TString strOutfile, Bool_t updateOutfile, Bool_t useRecPrim,
4629 TString strOutDir, Bool_t writeBgr, TString strBgrID)
4632 if((writeBgr && strBgrID.Length() == 0) || (!writeBgr && strBgrID.Length()>0) ){
4633 Printf("%s:%d -- inconsistent arguments to WriteBinShiftCorr FF/UE", (char*)__FILE__,__LINE__);
4637 TH1F* hCorrPt[fNJetPtSlices];
4638 TH1F* hCorrXi[fNJetPtSlices];
4639 TH1F* hCorrZ[fNJetPtSlices];
4641 TH1F* hdNdptTracksMCGen[fNJetPtSlices];
4642 TH1F* hdNdxiMCGen[fNJetPtSlices];
4643 TH1F* hdNdzMCGen[fNJetPtSlices];
4645 TH1F* hdNdptTracksMCRec[fNJetPtSlices];
4646 TH1F* hdNdxiMCRec[fNJetPtSlices];
4647 TH1F* hdNdzMCRec[fNJetPtSlices];
4649 TH1F* fh1FFJetPtMCGen = 0;
4650 TH1F* fh1FFJetPtMCRec = 0;
4652 TH2F* fh2FFTrackPtMCGen = 0;
4653 TH2F* fh2FFZMCGen = 0;
4654 TH2F* fh2FFXiMCGen = 0;
4656 TH2F* fh2FFTrackPtMCRec = 0;
4657 TH2F* fh2FFZMCRec = 0;
4658 TH2F* fh2FFXiMCRec = 0;
4662 TFile f(strInfile,"READ");
4665 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
4669 if(fDebug>0) Printf("%s:%d -- writeBinShiftCorr: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
4671 if(strdirGen && strdirGen.Length()) gDirectory->cd(strdirGen);
4675 if(strlistGen && strlistGen.Length()){
4677 if(!(listGen = (TList*) gDirectory->Get(strlistGen))){
4678 Printf("%s:%d -- error retrieving listGen %s from directory %s", (char*)__FILE__,__LINE__,strlistGen.Data(),strdirGen.Data());
4685 fh1FFJetPtMCGen = (TH1F*) listGen->FindObject("fh1FFJetPtGen");
4688 fh2FFTrackPtMCGen = (TH2F*) listGen->FindObject(Form("fh2FFTrackPt%sGen",strBgrID.Data()));
4689 fh2FFZMCGen = (TH2F*) listGen->FindObject(Form("fh2FFZ%sGen",strBgrID.Data()));
4690 fh2FFXiMCGen = (TH2F*) listGen->FindObject(Form("fh2FFXi%sGen",strBgrID.Data()));
4693 fh2FFTrackPtMCGen = (TH2F*) listGen->FindObject("fh2FFTrackPtGen");
4694 fh2FFZMCGen = (TH2F*) listGen->FindObject("fh2FFZGen");
4695 fh2FFXiMCGen = (TH2F*) listGen->FindObject("fh2FFXiGen");
4699 fh1FFJetPtMCGen = (TH1F*) gDirectory->Get("fh1FFJetPtGen");
4701 fh2FFTrackPtMCGen = (TH2F*) gDirectory->Get(Form("fh2FFTrackPt%sGen",strBgrID.Data()));
4702 fh2FFZMCGen = (TH2F*) gDirectory->Get(Form("fh2FFZ%sGen",strBgrID.Data()));
4703 fh2FFXiMCGen = (TH2F*) gDirectory->Get(Form("fh2FFXi%sGen",strBgrID.Data()));
4706 fh2FFTrackPtMCGen = (TH2F*) gDirectory->Get("fh2FFTrackPtGen");
4707 fh2FFZMCGen = (TH2F*) gDirectory->Get("fh2FFZGen");
4708 fh2FFXiMCGen = (TH2F*) gDirectory->Get("fh2FFXiGen");
4712 fh1FFJetPtMCGen->SetDirectory(0);
4713 fh2FFTrackPtMCGen->SetDirectory(0);
4714 fh2FFZMCGen->SetDirectory(0);
4715 fh2FFXiMCGen->SetDirectory(0);
4721 TFile g(strInfile,"READ");
4724 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
4728 if(fDebug>0) Printf("%s:%d -- writeBinShiftCorr: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
4730 if(strdirRec && strdirRec.Length()) gDirectory->cd(strdirRec);
4734 if(strlistRec && strlistRec.Length()){
4736 if(!(listRec = (TList*) gDirectory->Get(strlistRec))){
4737 Printf("%s:%d -- error retrieving listRec %s from directory %s", (char*)__FILE__,__LINE__,strlistRec.Data(),strdirRec.Data());
4744 fh1FFJetPtMCRec = (TH1F*) listRec->FindObject("fh1FFJetPtRecEffRec");
4747 fh2FFTrackPtMCRec = (TH2F*) listRec->FindObject(Form("fh2FFTrackPt%sRecEffRec",strBgrID.Data()));
4748 fh2FFZMCRec = (TH2F*) listRec->FindObject(Form("fh2FFZ%sRecEffRec",strBgrID.Data()));
4749 fh2FFXiMCRec = (TH2F*) listRec->FindObject(Form("fh2FFXi%sRecEffRec",strBgrID.Data()));
4752 fh2FFTrackPtMCRec = (TH2F*) listRec->FindObject("fh2FFTrackPtRecEffRec");
4753 fh2FFZMCRec = (TH2F*) listRec->FindObject("fh2FFZRecEffRec");
4754 fh2FFXiMCRec = (TH2F*) listRec->FindObject("fh2FFXiRecEffRec");
4758 fh1FFJetPtMCRec = (TH1F*) gDirectory->Get("fh1FFJetPtRecEffRec");
4760 fh2FFTrackPtMCRec = (TH2F*) gDirectory->Get(Form("fh2FFTrackPt%sRecEffRec",strBgrID.Data()));
4761 fh2FFZMCRec = (TH2F*) gDirectory->Get(Form("fh2FFZ%sRecEffRec",strBgrID.Data()));
4762 fh2FFXiMCRec = (TH2F*) gDirectory->Get(Form("fh2FFXi%sRecEffRec",strBgrID.Data()));
4765 fh2FFTrackPtMCRec = (TH2F*) gDirectory->Get("fh2FFTrackPtRecEffRec");
4766 fh2FFZMCRec = (TH2F*) gDirectory->Get("fh2FFZRecEffRec");
4767 fh2FFXiMCRec = (TH2F*) gDirectory->Get("fh2FFXiRecEffRec");
4773 fh1FFJetPtMCRec = (TH1F*) listRec->FindObject("fh1FFJetPtRecCuts");
4775 fh2FFTrackPtMCRec = (TH2F*) listRec->FindObject(Form("fh2FFTrackPt%sRecCuts",strBgrID.Data()));
4776 fh2FFZMCRec = (TH2F*) listRec->FindObject(Form("fh2FFZ%sRecCuts",strBgrID.Data()));
4777 fh2FFXiMCRec = (TH2F*) listRec->FindObject(Form("fh2FFXi%sRecCuts",strBgrID.Data()));
4780 fh2FFTrackPtMCRec = (TH2F*) listRec->FindObject("fh2FFTrackPtRecCuts");
4781 fh2FFZMCRec = (TH2F*) listRec->FindObject("fh2FFZRecCuts");
4782 fh2FFXiMCRec = (TH2F*) listRec->FindObject("fh2FFXiRecCuts");
4786 fh1FFJetPtMCRec = (TH1F*) gDirectory->Get("fh1FFJetPtRecCuts");
4788 fh2FFTrackPtMCRec = (TH2F*) gDirectory->Get(Form("fh2FFTrackPt%sRecCuts",strBgrID.Data()));
4789 fh2FFZMCRec = (TH2F*) gDirectory->Get(Form("fh2FFZ%sRecCuts",strBgrID.Data()));
4790 fh2FFXiMCRec = (TH2F*) gDirectory->Get(Form("fh2FFXi%sRecCuts",strBgrID.Data()));
4793 fh2FFTrackPtMCRec = (TH2F*) gDirectory->Get("fh2FFTrackPtRecCuts");
4794 fh2FFZMCRec = (TH2F*) gDirectory->Get("fh2FFZRecCuts");
4795 fh2FFXiMCRec = (TH2F*) gDirectory->Get("fh2FFXiRecCuts");
4801 fh1FFJetPtMCRec->SetDirectory(0);
4802 fh2FFTrackPtMCRec->SetDirectory(0);
4803 fh2FFZMCRec->SetDirectory(0);
4804 fh2FFXiMCRec->SetDirectory(0);
4808 // projections: FF for generated and reconstructed
4810 for(Int_t i=0; i<fNJetPtSlices; i++){
4812 Float_t jetPtLoLim = fJetPtSlices->At(i);
4813 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
4815 Int_t binLo = static_cast<Int_t>(fh2FFTrackPtMCGen->GetXaxis()->FindBin(jetPtLoLim));
4816 Int_t binUp = static_cast<Int_t>(fh2FFTrackPtMCGen->GetXaxis()->FindBin(jetPtUpLim))-1;
4818 if(binUp > fh2FFTrackPtMCGen->GetNbinsX()){
4819 Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
4823 TString strNameFFPtGen(Form("fh1FFTrackPtGenBinShift_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
4824 TString strNameFFZGen(Form("fh1FFZGenBinShift_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
4825 TString strNameFFXiGen(Form("fh1FFXiGenBinShift_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
4827 TString strNameFFPtRec(Form("fh1FFTrackPtRecBinShift_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
4828 TString strNameFFZRec(Form("fh1FFZRecBinShift_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
4829 TString strNameFFXiRec(Form("fh1FFXiRecBinShift_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
4832 strNameFFPtGen.Form("fh1TrackPtGenBinShiftBgr_%02d_%02d_%s",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim),strBgrID.Data());
4833 strNameFFZGen.Form("fh1ZGenBinShiftBgr_%02d_%02d_%s",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim),strBgrID.Data());
4834 strNameFFXiGen.Form("fh1XiGenBinShiftBgr_%02d_%02d_%s",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim),strBgrID.Data());
4836 strNameFFPtRec.Form("fh1TrackPtRecBinShiftBgr_%02d_%02d_%s",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim),strBgrID.Data());
4837 strNameFFZRec.Form("fh1ZRecBinShiftBgr_%02d_%02d_%s",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim),strBgrID.Data());
4838 strNameFFXiRec.Form("fh1XiRecBinShiftBgr_%02d_%02d_%s",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim),strBgrID.Data());
4842 // appendix 'unbinned' to avoid histos with same name after rebinning
4844 hdNdptTracksMCGen[i] = (TH1F*) fh2FFTrackPtMCGen->ProjectionY(strNameFFPtGen+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
4845 hdNdzMCGen[i] = (TH1F*) fh2FFZMCGen->ProjectionY(strNameFFZGen+"_unBinned",binLo,binUp,"o");
4846 hdNdxiMCGen[i] = (TH1F*) fh2FFXiMCGen->ProjectionY(strNameFFXiGen+"_unBinned",binLo,binUp,"o");
4848 hdNdptTracksMCRec[i] = (TH1F*) fh2FFTrackPtMCRec->ProjectionY(strNameFFPtRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
4849 hdNdzMCRec[i] = (TH1F*) fh2FFZMCRec->ProjectionY(strNameFFZRec+"_unBinned",binLo,binUp,"o");
4850 hdNdxiMCRec[i] = (TH1F*) fh2FFXiMCRec->ProjectionY(strNameFFXiRec+"_unBinned",binLo,binUp,"o");
4855 if(fNHistoBinsPt[i]) hdNdptTracksMCGen[i] = (TH1F*) hdNdptTracksMCGen[i]->Rebin(fNHistoBinsPt[i],strNameFFPtGen,fHistoBinsPt[i]->GetArray());
4856 if(fNHistoBinsZ[i]) hdNdzMCGen[i] = (TH1F*) hdNdzMCGen[i]->Rebin(fNHistoBinsZ[i],strNameFFZGen,fHistoBinsZ[i]->GetArray());
4857 if(fNHistoBinsXi[i]) hdNdxiMCGen[i] = (TH1F*) hdNdxiMCGen[i]->Rebin(fNHistoBinsXi[i],strNameFFXiGen,fHistoBinsXi[i]->GetArray());
4859 if(fNHistoBinsPt[i]) hdNdptTracksMCRec[i] = (TH1F*) hdNdptTracksMCRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtRec,fHistoBinsPt[i]->GetArray());
4860 if(fNHistoBinsZ[i]) hdNdzMCRec[i] = (TH1F*) hdNdzMCRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZRec,fHistoBinsZ[i]->GetArray());
4861 if(fNHistoBinsXi[i]) hdNdxiMCRec[i] = (TH1F*) hdNdxiMCRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiRec,fHistoBinsXi[i]->GetArray());
4864 hdNdptTracksMCGen[i]->SetNameTitle(strNameFFPtGen,"");
4865 hdNdzMCGen[i]->SetNameTitle(strNameFFZGen,"");
4866 hdNdxiMCGen[i]->SetNameTitle(strNameFFXiGen,"");
4868 hdNdptTracksMCRec[i]->SetNameTitle(strNameFFPtRec,"");
4869 hdNdzMCRec[i]->SetNameTitle(strNameFFZRec,"");
4870 hdNdxiMCRec[i]->SetNameTitle(strNameFFXiRec,"");
4874 Double_t nJetsBinGen = fh1FFJetPtMCGen->Integral(binLo,binUp);
4875 Double_t nJetsBinRec = fh1FFJetPtMCRec->Integral(binLo,binUp);
4877 NormalizeTH1(hdNdptTracksMCGen[i],nJetsBinGen);
4878 NormalizeTH1(hdNdzMCGen[i],nJetsBinGen);
4879 NormalizeTH1(hdNdxiMCGen[i],nJetsBinGen);
4881 NormalizeTH1(hdNdptTracksMCRec[i],nJetsBinRec);
4882 NormalizeTH1(hdNdzMCRec[i],nJetsBinRec);
4883 NormalizeTH1(hdNdxiMCRec[i],nJetsBinRec);
4885 // divide gen/rec : corr factor
4887 TString strNameCorrPt(Form("hBbBCorrPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
4888 TString strNameCorrZ(Form("hBbBCorrZ_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
4889 TString strNameCorrXi(Form("hBbBCorrXi_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
4892 strNameCorrPt.Form("hBbBCorrBgrPt_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
4893 strNameCorrZ.Form("hBbBCorrBgrZ_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
4894 strNameCorrXi.Form("hBbBCorrBgrXi_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
4897 hCorrPt[i] = (TH1F*) hdNdptTracksMCGen[i]->Clone(strNameCorrPt);
4898 hCorrPt[i]->Divide(hdNdptTracksMCGen[i],hdNdptTracksMCRec[i],1,1,"B"); // binominal errors
4900 hCorrXi[i] = (TH1F*) hdNdxiMCGen[i]->Clone(strNameCorrXi);
4901 hCorrXi[i]->Divide(hdNdxiMCGen[i],hdNdxiMCRec[i],1,1,"B"); // binominal errors
4903 hCorrZ[i] = (TH1F*) hdNdzMCGen[i]->Clone(strNameCorrZ);
4904 hCorrZ[i]->Divide(hdNdzMCGen[i],hdNdzMCRec[i],1,1,"B"); // binominal errors
4909 TString outfileOption = "RECREATE";
4910 if(updateOutfile) outfileOption = "UPDATE";
4912 TFile out(strOutfile,outfileOption);
4915 Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
4919 if(fDebug>0) Printf("%s:%d -- write bin shift correction to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
4921 if(strOutDir && strOutDir.Length()){
4924 if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd();
4926 dir = out.mkdir(strOutDir);
4931 for(Int_t i=0; i<fNJetPtSlices; i++){
4933 hdNdptTracksMCGen[i]->Write();
4934 hdNdxiMCGen[i]->Write();
4935 hdNdzMCGen[i]->Write();
4937 hdNdptTracksMCRec[i]->Write();
4938 hdNdxiMCRec[i]->Write();
4939 hdNdzMCRec[i]->Write();
4941 hCorrPt[i]->Write();
4942 hCorrXi[i]->Write();
4948 delete fh1FFJetPtMCGen;
4949 delete fh1FFJetPtMCRec;
4951 delete fh2FFTrackPtMCGen;
4953 delete fh2FFXiMCGen;
4955 delete fh2FFTrackPtMCRec;
4957 delete fh2FFXiMCRec;
4960 //________________________________________________________________________________________________________________________________
4961 void AliFragmentationFunctionCorrections::ReadBgrBinShiftCorr(TString strfile, TString strBgrID, TString strdir, TString strlist)
4964 ReadBinShiftCorr(strfile, strdir, strlist, kTRUE, strBgrID);
4967 //___________________________________________________________________________________________________________________________________________
4968 void AliFragmentationFunctionCorrections::ReadBinShiftCorr(TString strfile, TString strdir, TString strlist, Bool_t readBgr, TString strBgrID)
4971 if((readBgr && strBgrID.Length() == 0) || (!readBgr && strBgrID.Length()>0) ){
4972 Printf("%s:%d -- inconsistent arguments to ReadBinShiftCorr FF/UE", (char*)__FILE__,__LINE__);
4976 // temporary histos to hold histos from file
4977 TH1F* hCorrPt[fNJetPtSlices];
4978 TH1F* hCorrZ[fNJetPtSlices];
4979 TH1F* hCorrXi[fNJetPtSlices];
4981 for(Int_t i=0; i<fNJetPtSlices; i++) hCorrPt[i] = 0;
4982 for(Int_t i=0; i<fNJetPtSlices; i++) hCorrZ[i] = 0;
4983 for(Int_t i=0; i<fNJetPtSlices; i++) hCorrXi[i] = 0;
4985 TFile f(strfile,"READ");
4988 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
4992 if(fDebug>0) Printf("%s:%d -- read FF / UE bin shift correction from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
4994 if(strdir && strdir.Length()) gDirectory->cd(strdir);
4998 if(strlist && strlist.Length()){
5000 if(!(list = (TList*) gDirectory->Get(strlist))){
5001 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
5006 for(Int_t i=0; i<fNJetPtSlices; i++){
5008 Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
5009 Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
5011 TString strNameCorrPt(Form("hBbBCorrPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
5012 TString strNameCorrZ(Form("hBbBCorrZ_%02d_%02d",jetPtLoLim,jetPtUpLim));
5013 TString strNameCorrXi(Form("hBbBCorrXi_%02d_%02d",jetPtLoLim,jetPtUpLim));
5016 strNameCorrPt.Form("hBbBCorrBgrPt_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
5017 strNameCorrZ.Form("hBbBCorrBgrZ_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
5018 strNameCorrXi.Form("hBbBCorrBgrXi_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
5023 hCorrPt[i] = (TH1F*) list->FindObject(strNameCorrPt);
5024 hCorrZ[i] = (TH1F*) list->FindObject(strNameCorrZ);
5025 hCorrXi[i] = (TH1F*) list->FindObject(strNameCorrXi);
5028 hCorrPt[i] = (TH1F*) gDirectory->Get(strNameCorrPt);
5029 hCorrZ[i] = (TH1F*) gDirectory->Get(strNameCorrZ);
5030 hCorrXi[i] = (TH1F*) gDirectory->Get(strNameCorrXi);
5034 Printf("%s:%d -- error retrieving bin by bin correction %s", (char*)__FILE__,__LINE__,strNameCorrPt.Data());
5038 Printf("%s:%d -- error retrieving bin by bin correction %s", (char*)__FILE__,__LINE__,strNameCorrZ.Data());
5042 Printf("%s:%d -- error retrieving bin by bin correction %s", (char*)__FILE__,__LINE__,strNameCorrXi.Data());
5046 if(fNHistoBinsPt[i]) hCorrPt[i] = (TH1F*) hCorrPt[i]->Rebin(fNHistoBinsPt[i],strNameCorrPt+"_rebin",fHistoBinsPt[i]->GetArray());
5047 if(fNHistoBinsZ[i]) hCorrZ[i] = (TH1F*) hCorrZ[i]->Rebin(fNHistoBinsZ[i],strNameCorrZ+"_rebin",fHistoBinsZ[i]->GetArray());
5048 if(fNHistoBinsXi[i]) hCorrXi[i] = (TH1F*) hCorrXi[i]->Rebin(fNHistoBinsXi[i],strNameCorrXi+"_rebin",fHistoBinsXi[i]->GetArray());
5050 if(hCorrPt[i]) hCorrPt[i]->SetDirectory(0);
5051 if(hCorrZ[i]) hCorrZ[i]->SetDirectory(0);
5052 if(hCorrXi[i]) hCorrXi[i]->SetDirectory(0);
5054 } // jet slices loop
5059 for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
5060 if(hCorrPt[i]) new(fh1BbBBgrPt[i]) TH1F(*hCorrPt[i]);
5061 if(hCorrZ[i]) new(fh1BbBBgrZ[i]) TH1F(*hCorrZ[i]);
5062 if(hCorrXi[i]) new(fh1BbBBgrXi[i]) TH1F(*hCorrXi[i]);
5066 for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
5067 if(hCorrPt[i]) new(fh1BbBPt[i]) TH1F(*hCorrPt[i]);
5068 if(hCorrZ[i]) new(fh1BbBZ[i]) TH1F(*hCorrZ[i]);
5069 if(hCorrXi[i]) new(fh1BbBXi[i]) TH1F(*hCorrXi[i]);
5074 // ________________________________________________
5075 void AliFragmentationFunctionCorrections::BbBCorr()
5077 // apply bin-by-bin correction
5079 AddCorrectionLevel("BbB");
5081 for(Int_t i=0; i<fNJetPtSlices; i++){
5083 TH1F* histPt = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i);
5084 TH1F* histZ = fCorrFF[fNCorrectionLevels-2]->GetZ(i);
5085 TH1F* histXi = fCorrFF[fNCorrectionLevels-2]->GetXi(i);
5087 TString histNamePt = histPt->GetName();
5088 TString histNameZ = histZ->GetName();
5089 TString histNameXi = histXi->GetName();
5091 TH1F* hFFTrackPtBbBCorr = (TH1F*) histPt->Clone(histNamePt);
5092 hFFTrackPtBbBCorr->Multiply(histPt,fh1BbBPt[i],1,1,"");
5094 TH1F* hFFZBbBCorr = (TH1F*) histZ->Clone(histNameZ);
5095 hFFZBbBCorr->Multiply(histZ,fh1BbBZ[i],1,1,"");
5097 TH1F* hFFXiBbBCorr = (TH1F*) histXi->Clone(histNameXi);
5098 hFFXiBbBCorr->Multiply(histXi,fh1BbBXi[i],1,1,"");
5100 fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hFFTrackPtBbBCorr,hFFZBbBCorr,hFFXiBbBCorr);
5104 // ___________________________________________________
5105 void AliFragmentationFunctionCorrections::BbBCorrBgr()
5107 // apply bin-by-bin correction
5109 AddCorrectionLevelBgr("BbB");
5111 for(Int_t i=0; i<fNJetPtSlices; i++){
5113 TH1F* histPt = fCorrBgr[fNCorrectionLevelsBgr-2]->GetTrackPt(i);
5114 TH1F* histZ = fCorrBgr[fNCorrectionLevelsBgr-2]->GetZ(i);
5115 TH1F* histXi = fCorrBgr[fNCorrectionLevelsBgr-2]->GetXi(i);
5117 TString histNamePt = histPt->GetName();
5118 TString histNameZ = histZ->GetName();
5119 TString histNameXi = histXi->GetName();
5121 TH1F* hBgrTrackPtBbBCorr = (TH1F*) histPt->Clone(histNamePt);
5122 hBgrTrackPtBbBCorr->Multiply(histPt,fh1BbBBgrPt[i],1,1,"");
5124 TH1F* hBgrZBbBCorr = (TH1F*) histZ->Clone(histNameZ);
5125 hBgrZBbBCorr->Multiply(histZ,fh1BbBBgrZ[i],1,1,"");
5127 TH1F* hBgrXiBbBCorr = (TH1F*) histXi->Clone(histNameXi);
5128 hBgrXiBbBCorr->Multiply(histXi,fh1BbBBgrXi[i],1,1,"");
5130 fCorrBgr[fNCorrectionLevelsBgr-1]->AddCorrHistos(i,hBgrTrackPtBbBCorr,hBgrZBbBCorr,hBgrXiBbBCorr);
5134 //_______________________________________________________________________________________________________
5135 void AliFragmentationFunctionCorrections::ReadFoldingCorr(TString strfile, TString strdir, TString strlist)
5138 // temporary histos to hold histos from file
5139 TH1F* hCorrPt[fNJetPtSlices];
5140 TH1F* hCorrZ[fNJetPtSlices];
5141 TH1F* hCorrXi[fNJetPtSlices];
5143 for(Int_t i=0; i<fNJetPtSlices; i++) hCorrPt[i] = 0;
5144 for(Int_t i=0; i<fNJetPtSlices; i++) hCorrZ[i] = 0;
5145 for(Int_t i=0; i<fNJetPtSlices; i++) hCorrXi[i] = 0;
5147 TFile f(strfile,"READ");
5150 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
5154 if(fDebug>0) Printf("%s:%d -- read bin shift correction from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
5156 if(strdir && strdir.Length()) gDirectory->cd(strdir);
5160 if(strlist && strlist.Length()){
5162 if(!(list = (TList*) gDirectory->Get(strlist))){
5163 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
5168 for(Int_t i=0; i<fNJetPtSlices; i++){
5170 Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
5171 Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
5173 TString strNameCorrPt(Form("hCorrFoldingPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
5174 TString strNameCorrZ(Form("hCorrFoldingZ_%02d_%02d",jetPtLoLim,jetPtUpLim));
5175 TString strNameCorrXi(Form("hCorrFoldingXi_%02d_%02d",jetPtLoLim,jetPtUpLim));
5179 hCorrPt[i] = (TH1F*) list->FindObject(strNameCorrPt);
5180 hCorrZ[i] = (TH1F*) list->FindObject(strNameCorrZ);
5181 hCorrXi[i] = (TH1F*) list->FindObject(strNameCorrXi);
5184 hCorrPt[i] = (TH1F*) gDirectory->Get(strNameCorrPt);
5185 hCorrZ[i] = (TH1F*) gDirectory->Get(strNameCorrZ);
5186 hCorrXi[i] = (TH1F*) gDirectory->Get(strNameCorrXi);
5190 //Printf("%s:%d -- error retrieving folding correction %s", (char*)__FILE__,__LINE__,strNameCorrPt.Data());
5194 //Printf("%s:%d -- error retrieving folding correction %s", (char*)__FILE__,__LINE__,strNameCorrZ.Data());
5198 Printf("%s:%d -- error retrieving folding correction %s", (char*)__FILE__,__LINE__,strNameCorrXi.Data());
5201 if(hCorrPt[i]) hCorrPt[i]->SetDirectory(0);
5202 if(hCorrZ[i]) hCorrZ[i]->SetDirectory(0);
5203 if(hCorrXi[i]) hCorrXi[i]->SetDirectory(0);
5205 } // jet slices loop
5209 for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
5210 if(hCorrPt[i]) new(fh1FoldingCorrPt[i]) TH1F(*hCorrPt[i]);
5211 if(hCorrZ[i]) new(fh1FoldingCorrZ[i]) TH1F(*hCorrZ[i]);
5212 if(hCorrXi[i]) new(fh1FoldingCorrXi[i]) TH1F(*hCorrXi[i]);
5216 // ___________________________________________________
5217 void AliFragmentationFunctionCorrections::FoldingCorr()
5219 // apply bin-by-bin correction
5221 AddCorrectionLevel("FoldingCorr");
5223 for(Int_t i=0; i<fNJetPtSlices; i++){
5225 TH1F* histPt = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i);
5226 TH1F* histZ = fCorrFF[fNCorrectionLevels-2]->GetZ(i);
5227 TH1F* histXi = fCorrFF[fNCorrectionLevels-2]->GetXi(i);
5229 TString histNamePt = histPt->GetName();
5230 TString histNameZ = histZ->GetName();
5231 TString histNameXi = histXi->GetName();
5233 std::cout<<" foldingCorr: i "<<i<<" corr pt "<<fh1FoldingCorrPt[i]<<" z "<<fh1FoldingCorrZ[i]<<" xi "
5234 <<fh1FoldingCorrXi[i]<<std::endl;
5236 std::cout<<" foldingCorr: i "<<i<<" mean corr pt "<<fh1FoldingCorrPt[i]->GetMean()<<" z "<<fh1FoldingCorrZ[i]->GetMean()<<" xi "
5237 <<fh1FoldingCorrXi[i]->GetMean()<<std::endl;
5241 TH1F* hFFTrackPtFoldingCorr = (TH1F*) histPt->Clone(histNamePt);
5242 if(fh1FoldingCorrPt[i] && fh1FoldingCorrPt[i]->GetMean()>0) hFFTrackPtFoldingCorr->Multiply(histPt,fh1FoldingCorrPt[i],1,1,"");
5243 else hFFTrackPtFoldingCorr->Reset();
5245 TH1F* hFFZFoldingCorr = (TH1F*) histZ->Clone(histNameZ);
5246 if(fh1FoldingCorrZ[i] && fh1FoldingCorrZ[i]->GetMean()>0) hFFZFoldingCorr->Multiply(histZ,fh1FoldingCorrZ[i],1,1,"");
5247 else hFFZFoldingCorr->Reset();
5249 TH1F* hFFXiFoldingCorr = (TH1F*) histXi->Clone(histNameXi);
5250 if(fh1FoldingCorrXi[i]&& fh1FoldingCorrXi[i]->GetMean()>0) hFFXiFoldingCorr->Multiply(histXi,fh1FoldingCorrXi[i],1,1,"");
5251 else hFFXiFoldingCorr->Reset();
5253 fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hFFTrackPtFoldingCorr,hFFZFoldingCorr,hFFXiFoldingCorr);
5257 //________________________________________________________________________________________________________________________________
5258 void AliFragmentationFunctionCorrections::ReadBgrJetSecCorr(TString strfile, TString strBgrID, TString strdir, TString strlist,
5259 Bool_t useScaledStrangeness)
5262 ReadJetSecCorr(strfile, strdir, strlist, useScaledStrangeness, kTRUE, strBgrID);
5265 //_______________________________________________________________________________________________________
5266 void AliFragmentationFunctionCorrections::ReadJetSecCorr(TString strfile, TString strdir, TString strlist, Bool_t useScaledStrangeness,
5267 Bool_t readBgr, TString strBgrID){
5269 // read reconstruction efficiency from file
5270 // argument strlist optional - read from directory strdir if not specified
5272 // temporary histos to hold histos from file
5273 TH1F* hCorrPt[fNJetPtSlices];
5274 TH1F* hCorrZ[fNJetPtSlices];
5275 TH1F* hCorrXi[fNJetPtSlices];
5277 for(Int_t i=0; i<fNJetPtSlices; i++) hCorrPt[i] = 0;
5278 for(Int_t i=0; i<fNJetPtSlices; i++) hCorrZ[i] = 0;
5279 for(Int_t i=0; i<fNJetPtSlices; i++) hCorrXi[i] = 0;
5281 TFile f(strfile,"READ");
5284 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
5288 if(fDebug>0) Printf("%s:%d -- read secondary correction from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
5290 if(strdir && strdir.Length()) gDirectory->cd(strdir);
5294 if(strlist && strlist.Length()){
5296 if(!(list = (TList*) gDirectory->Get(strlist))){
5297 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
5302 for(Int_t i=0; i<fNJetPtSlices; i++){
5304 Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
5305 Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
5307 TString strNameCorrPt("");
5308 TString strNameCorrZ("");
5309 TString strNameCorrXi("");
5312 strNameCorrPt.Form("hSecCorrBgrPt_%02d_%02d_%s",jetPtLoLim,jetPtUpLim,strBgrID.Data());
5313 strNameCorrZ.Form("hSecCorrBgrZ_%02d_%02d_%s",jetPtLoLim,jetPtUpLim,strBgrID.Data());
5314 strNameCorrXi.Form("hSecCorrBgrXi_%02d_%02d_%s",jetPtLoLim,jetPtUpLim,strBgrID.Data());
5316 if(!useScaledStrangeness){
5317 Printf("%s:%d -- readJetSecCorr bgr: using naive Pythia strangenss ", (char*)__FILE__,__LINE__);
5318 strNameCorrPt.Form("hSecCorrBgrPt_nonSc_%02d_%02d_%s",jetPtLoLim,jetPtUpLim,strBgrID.Data());
5319 strNameCorrZ.Form("hSecCorrBgrZ_nonSc_%02d_%02d_%s",jetPtLoLim,jetPtUpLim,strBgrID.Data());
5320 strNameCorrXi.Form("hSecCorrBgrXi_nonSc_%02d_%02d_%s",jetPtLoLim,jetPtUpLim,strBgrID.Data());
5324 strNameCorrPt.Form("hSecCorrPt_%02d_%02d",jetPtLoLim,jetPtUpLim);
5325 strNameCorrZ.Form("hSecCorrZ_%02d_%02d",jetPtLoLim,jetPtUpLim);
5326 strNameCorrXi.Form("hSecCorrXi_%02d_%02d",jetPtLoLim,jetPtUpLim);
5328 if(!useScaledStrangeness){
5329 Printf("%s:%d -- readJetSecCorr: using naive Pythia strangenss ", (char*)__FILE__,__LINE__);
5330 strNameCorrPt.Form("hSecCorrPt_nonSc_%02d_%02d",jetPtLoLim,jetPtUpLim);
5331 strNameCorrZ.Form("hSecCorrZ_nonSc_%02d_%02d",jetPtLoLim,jetPtUpLim);
5332 strNameCorrXi.Form("hSecCorrXi_nonSc_%02d_%02d",jetPtLoLim,jetPtUpLim);
5337 hCorrPt[i] = (TH1F*) list->FindObject(strNameCorrPt);
5338 hCorrZ[i] = (TH1F*) list->FindObject(strNameCorrZ);
5339 hCorrXi[i] = (TH1F*) list->FindObject(strNameCorrXi);
5342 hCorrPt[i] = (TH1F*) gDirectory->Get(strNameCorrPt);
5343 hCorrZ[i] = (TH1F*) gDirectory->Get(strNameCorrZ);
5344 hCorrXi[i] = (TH1F*) gDirectory->Get(strNameCorrXi);
5348 Printf("%s:%d -- error retrieving secondaries correction %s", (char*)__FILE__,__LINE__,strNameCorrPt.Data());
5352 Printf("%s:%d -- error retrieving secondaries correction %s", (char*)__FILE__,__LINE__,strNameCorrZ.Data());
5356 Printf("%s:%d -- error retrieving secondaries correction %s", (char*)__FILE__,__LINE__,strNameCorrXi.Data());
5359 if(fNHistoBinsPt[i]) hCorrPt[i] = (TH1F*) hCorrPt[i]->Rebin(fNHistoBinsPt[i],strNameCorrPt+"_rebin",fHistoBinsPt[i]->GetArray());
5360 if(fNHistoBinsZ[i]) hCorrZ[i] = (TH1F*) hCorrZ[i]->Rebin(fNHistoBinsZ[i],strNameCorrZ+"_rebin",fHistoBinsZ[i]->GetArray());
5361 if(fNHistoBinsXi[i]) hCorrXi[i] = (TH1F*) hCorrXi[i]->Rebin(fNHistoBinsXi[i],strNameCorrXi+"_rebin",fHistoBinsXi[i]->GetArray());
5363 if(hCorrPt[i]) hCorrPt[i]->SetDirectory(0);
5364 if(hCorrZ[i]) hCorrZ[i]->SetDirectory(0);
5365 if(hCorrXi[i]) hCorrXi[i]->SetDirectory(0);
5367 } // jet slices loop
5371 for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
5374 if(hCorrPt[i]) new(fh1SecCorrBgrPt[i]) TH1F(*hCorrPt[i]);
5375 if(hCorrZ[i]) new(fh1SecCorrBgrZ[i]) TH1F(*hCorrZ[i]);
5376 if(hCorrXi[i]) new(fh1SecCorrBgrXi[i]) TH1F(*hCorrXi[i]);
5379 if(hCorrPt[i]) new(fh1SecCorrPt[i]) TH1F(*hCorrPt[i]);
5380 if(hCorrZ[i]) new(fh1SecCorrZ[i]) TH1F(*hCorrZ[i]);
5381 if(hCorrXi[i]) new(fh1SecCorrXi[i]) TH1F(*hCorrXi[i]);
5386 // ___________________________________________________
5387 void AliFragmentationFunctionCorrections::JetSecCorr()
5389 // apply secondaries correction
5391 AddCorrectionLevel("SecCorr");
5393 Printf("%s:%d -- apply jet secondaries correction", (char*)__FILE__,__LINE__);
5395 for(Int_t i=0; i<fNJetPtSlices; i++){
5397 TH1F* histPt = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i);
5398 TH1F* histZ = fCorrFF[fNCorrectionLevels-2]->GetZ(i);
5399 TH1F* histXi = fCorrFF[fNCorrectionLevels-2]->GetXi(i);
5401 TString histNamePt = histPt->GetName();
5402 TString histNameZ = histZ->GetName();
5403 TString histNameXi = histXi->GetName();
5405 TH1F* hFFTrackPtSecCorr = (TH1F*) histPt->Clone(histNamePt);
5406 hFFTrackPtSecCorr->Multiply(histPt,fh1SecCorrPt[i],1,1,"");
5408 TH1F* hFFZSecCorr = (TH1F*) histZ->Clone(histNameZ);
5409 hFFZSecCorr->Multiply(histZ,fh1SecCorrZ[i],1,1,"");
5411 TH1F* hFFXiSecCorr = (TH1F*) histXi->Clone(histNameXi);
5412 hFFXiSecCorr->Multiply(histXi,fh1SecCorrXi[i],1,1,"");
5414 fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hFFTrackPtSecCorr,hFFZSecCorr,hFFXiSecCorr);
5420 // ___________________________________________________
5421 void AliFragmentationFunctionCorrections::JetSecCorrBgr()
5423 // apply secondaries correction to UE
5425 AddCorrectionLevelBgr("SecCorr");
5427 Printf("%s:%d -- apply jet secondaries correction", (char*)__FILE__,__LINE__);
5429 for(Int_t i=0; i<fNJetPtSlices; i++){
5431 TH1F* histPt = fCorrBgr[fNCorrectionLevelsBgr-2]->GetTrackPt(i);
5432 TH1F* histZ = fCorrBgr[fNCorrectionLevelsBgr-2]->GetZ(i);
5433 TH1F* histXi = fCorrBgr[fNCorrectionLevelsBgr-2]->GetXi(i);
5435 TString histNamePt = histPt->GetName();
5436 TString histNameZ = histZ->GetName();
5437 TString histNameXi = histXi->GetName();
5439 TH1F* hFFTrackPtSecCorr = (TH1F*) histPt->Clone(histNamePt);
5440 hFFTrackPtSecCorr->Multiply(histPt,fh1SecCorrBgrPt[i],1,1,"");
5442 TH1F* hFFZSecCorr = (TH1F*) histZ->Clone(histNameZ);
5443 hFFZSecCorr->Multiply(histZ,fh1SecCorrBgrZ[i],1,1,"");
5445 TH1F* hFFXiSecCorr = (TH1F*) histXi->Clone(histNameXi);
5446 hFFXiSecCorr->Multiply(histXi,fh1SecCorrBgrXi[i],1,1,"");
5448 fCorrBgr[fNCorrectionLevelsBgr-1]->AddCorrHistos(i,hFFTrackPtSecCorr,hFFZSecCorr,hFFXiSecCorr);
5453 //________________________________________________________________________________________________________________
5454 void AliFragmentationFunctionCorrections::WriteBinShiftCorrSinglePt(TString strInfile, TString strIDGen, TString strIDRec,
5455 TString strOutfile, Bool_t updateOutfile, Bool_t useRecPrim,
5458 TString strdirGen = "PWGJE_FragmentationFunction_" + strIDGen;
5459 TString strlistGen = "fracfunc_" + strIDGen;
5461 TString strdirRec = "PWGJE_FragmentationFunction_" + strIDRec;
5462 TString strlistRec = "fracfunc_" + strIDRec;
5464 WriteBinShiftCorrSinglePt(strInfile,strdirGen,strlistGen,strdirRec,strlistRec,strOutfile,updateOutfile,useRecPrim,strOutDir);
5467 //___________________________________________________________________________________________________________________________________
5468 void AliFragmentationFunctionCorrections::WriteBinShiftCorrSinglePt(TString strInfile, TString strdirGen, TString strlistGen,
5469 TString strdirRec, TString strlistRec,
5470 TString strOutfile, Bool_t updateOutfile, Bool_t useRecPrim,
5475 TH1F* hdNdptTracksMCGen = 0;
5476 TH1F* hdNdptTracksMCRec = 0;
5480 TFile f(strInfile,"READ");
5483 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
5487 if(fDebug>0) Printf("%s:%d -- writeBinShiftCorrSinglePt: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
5489 if(strdirGen && strdirGen.Length()) gDirectory->cd(strdirGen);
5493 if(strlistGen && strlistGen.Length()){
5495 if(!(listGen = (TList*) gDirectory->Get(strlistGen))){
5496 Printf("%s:%d -- error retrieving listGen %s from directory %s", (char*)__FILE__,__LINE__,strlistGen.Data(),strdirGen.Data());
5502 hdNdptTracksMCGen = (TH1F*) listGen->FindObject("fh1TrackQAPtGen");
5505 hdNdptTracksMCGen = (TH1F*) gDirectory->Get("fh1TrackQAPtGen");
5508 hdNdptTracksMCGen->SetDirectory(0);
5514 TFile g(strInfile,"READ");
5517 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
5521 if(fDebug>0) Printf("%s:%d -- writeBinShiftCorrSinglePt: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
5523 if(strdirRec && strdirRec.Length()) gDirectory->cd(strdirRec);
5527 if(strlistRec && strlistRec.Length()){
5529 if(!(listRec = (TList*) gDirectory->Get(strlistRec))){
5530 Printf("%s:%d -- error retrieving listRec %s from directory %s", (char*)__FILE__,__LINE__,strlistRec.Data(),strdirRec.Data());
5538 hdNdptTracksMCRec = (TH1F*) listRec->FindObject("fh1TrackQAPtRecEffRec");
5541 hdNdptTracksMCRec = (TH1F*) gDirectory->Get("fh1TrackQAPtRecEffRec");
5546 hdNdptTracksMCRec = (TH1F*) listRec->FindObject("fh1TrackQAPtRecCuts");
5549 hdNdptTracksMCRec = (TH1F*) gDirectory->Get("fh1TrackQAPtRecCuts");
5553 hdNdptTracksMCRec->SetDirectory(0);
5557 TString strNamePtGen = "fh1SinglePtGenBbB";
5558 TString strNamePtRec = "fh1SinglePtRecBbB";
5561 if(fNHistoBinsSinglePt) hdNdptTracksMCGen = (TH1F*) hdNdptTracksMCGen->Rebin(fNHistoBinsSinglePt,strNamePtGen+"_rebin",fHistoBinsSinglePt->GetArray());
5562 if(fNHistoBinsSinglePt) hdNdptTracksMCRec = (TH1F*) hdNdptTracksMCRec->Rebin(fNHistoBinsSinglePt,strNamePtRec+"_rebin",fHistoBinsSinglePt->GetArray());
5564 hdNdptTracksMCGen->SetNameTitle(strNamePtGen,"");
5565 hdNdptTracksMCRec->SetNameTitle(strNamePtRec,"");
5568 TString strTitCorr = "hBbBCorrSinglePt";
5569 if(useRecPrim) strTitCorr = "hBbBCorrRecPrimSinglePt";
5571 hCorrPt = (TH1F*) hdNdptTracksMCGen->Clone(strTitCorr);
5572 hCorrPt->Divide(hdNdptTracksMCGen,hdNdptTracksMCRec,1,1,"B"); // binominal errors
5576 TString outfileOption = "RECREATE";
5577 if(updateOutfile) outfileOption = "UPDATE";
5579 TFile out(strOutfile,outfileOption);
5582 Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
5586 if(fDebug>0) Printf("%s:%d -- write bin shift correction to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
5588 if(strOutDir && strOutDir.Length()){
5591 if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd();
5593 dir = out.mkdir(strOutDir);
5599 hdNdptTracksMCGen->Write();
5600 hdNdptTracksMCRec->Write();
5605 delete hdNdptTracksMCGen;
5606 delete hdNdptTracksMCRec;
5610 //___________________________________________________________________________________________________________________________________
5611 void AliFragmentationFunctionCorrections::ReadBinShiftCorrSinglePt(TString strfile, TString strdir, TString strlist, Bool_t useRecPrim)
5613 // read reconstruction efficiency from file
5614 // argument strlist optional - read from directory strdir if not specified
5616 // temporary histos to hold histos from file
5617 TH1F* hBbBCorrPt = 0;
5619 TFile f(strfile,"READ");
5622 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
5626 if(fDebug>0) Printf("%s:%d -- read BbB corr from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
5628 if(strdir && strdir.Length()) gDirectory->cd(strdir);
5632 if(strlist && strlist.Length()){
5634 if(!(list = (TList*) gDirectory->Get(strlist))){
5635 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
5641 TString strNameBbBCorrPt = "hBbBCorrSinglePt";
5642 if(useRecPrim) strNameBbBCorrPt = "hBbBCorrRecPrimSinglePt";
5645 hBbBCorrPt = (TH1F*) list->FindObject(strNameBbBCorrPt);
5648 hBbBCorrPt = (TH1F*) gDirectory->Get(strNameBbBCorrPt);
5652 Printf("%s:%d -- error retrieving BbB corr single pt %s", (char*)__FILE__,__LINE__,strNameBbBCorrPt.Data());
5656 if(fNHistoBinsPt && hBbBCorrPt)
5657 hBbBCorrPt = (TH1F*) hBbBCorrPt->Rebin(fNHistoBinsSinglePt,strNameBbBCorrPt+"_rebin",fHistoBinsSinglePt->GetArray());
5659 if(hBbBCorrPt) hBbBCorrPt->SetDirectory(0);
5663 fh1BbBCorrSinglePt = hBbBCorrPt;
5666 // ------------------------------------------------------------------
5668 void AliFragmentationFunctionCorrections::BbBCorrSinglePt()
5670 // apply efficiency correction to inclusive track pt spec
5672 AddCorrectionLevelSinglePt("BbB");
5674 TH1F* histPt = fCorrSinglePt[fNCorrectionLevelsSinglePt-2]->GetTrackPt(0);
5676 if(histPt->GetNbinsX() != fh1BbBCorrSinglePt->GetNbinsX()){
5677 Printf("%s:%d: inconsistency pt spec and BbB corr binning ", (char*)__FILE__,__LINE__);
5681 TString histNamePt = histPt->GetName();
5682 TH1F* hTrackPtBbBCorr = (TH1F*) histPt->Clone(histNamePt);
5684 hTrackPtBbBCorr->Multiply(histPt,fh1BbBCorrSinglePt,1,1,"");
5686 fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->AddCorrHistos(0,hTrackPtBbBCorr);