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,Double_t scaleFacBgrRec)
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,scaleFacBgrRec);
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,Double_t scaleFacBgrRec)
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 // scale fac for perp2 bgr
3167 if(writeBgr && scaleFacBgrRec && (scaleFacBgrRec != 1)) nJetsBin /= scaleFacBgrRec;
3169 NormalizeTH1(hdNdptTracksMCPrimRec[i],nJetsBin);
3170 NormalizeTH1(hdNdzMCPrimRec[i],nJetsBin);
3171 NormalizeTH1(hdNdxiMCPrimRec[i],nJetsBin);
3173 NormalizeTH1(hdNdptTracksMCSecRecNS[i],nJetsBin);
3174 NormalizeTH1(hdNdzMCSecRecNS[i],nJetsBin);
3175 NormalizeTH1(hdNdxiMCSecRecNS[i],nJetsBin);
3177 NormalizeTH1(hdNdptTracksMCSecRecS[i],nJetsBin);
3178 NormalizeTH1(hdNdzMCSecRecS[i],nJetsBin);
3179 NormalizeTH1(hdNdxiMCSecRecS[i],nJetsBin);
3181 NormalizeTH1(hdNdptTracksMCSecRecSsc[i],nJetsBin);
3182 NormalizeTH1(hdNdzMCSecRecSsc[i],nJetsBin);
3183 NormalizeTH1(hdNdxiMCSecRecSsc[i],nJetsBin);
3186 // divide prim / (prim+sec) : corr factor
3187 TString strNameSecCorrPt(Form("hSecCorrPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3188 TString strNameSecCorrZ(Form("hSecCorrZ_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3189 TString strNameSecCorrXi(Form("hSecCorrXi_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3192 strNameSecCorrPt.Form("hSecCorrBgrPt_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
3193 strNameSecCorrZ.Form("hSecCorrBgrZ_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
3194 strNameSecCorrXi.Form("hSecCorrBgrXi_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
3197 hSecCorrPt[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Clone(strNameSecCorrPt);
3198 TH1F* hSumPrimSecPt = (TH1F*) hdNdptTracksMCPrimRec[i]->Clone("hSumPrimSecPt");
3199 hSumPrimSecPt->Add(hdNdptTracksMCSecRecNS[i]);
3200 hSumPrimSecPt->Add(hdNdptTracksMCSecRecSsc[i]);
3201 hSecCorrPt[i]->Divide(hdNdptTracksMCPrimRec[i],hSumPrimSecPt,1,1,"B"); // binominal errors
3203 hSecCorrXi[i] = (TH1F*) hdNdxiMCPrimRec[i]->Clone(strNameSecCorrXi);
3204 TH1F* hSumPrimSecXi = (TH1F*) hdNdxiMCPrimRec[i]->Clone("hSumPrimSecXi");
3205 hSumPrimSecXi->Add(hdNdxiMCSecRecNS[i]);
3206 hSumPrimSecXi->Add(hdNdxiMCSecRecSsc[i]);
3207 hSecCorrXi[i]->Divide(hdNdxiMCPrimRec[i],hSumPrimSecXi,1,1,"B"); // binominal errors
3209 hSecCorrZ[i] = (TH1F*) hdNdzMCPrimRec[i]->Clone(strNameSecCorrZ);
3210 TH1F* hSumPrimSecZ = (TH1F*) hdNdzMCPrimRec[i]->Clone("hSumPrimSecZ");
3211 hSumPrimSecZ->Add(hdNdzMCSecRecNS[i]);
3212 hSumPrimSecZ->Add(hdNdzMCSecRecSsc[i]);
3213 hSecCorrZ[i]->Divide(hdNdzMCPrimRec[i],hSumPrimSecZ,1,1,"B"); // binominal errors
3215 // the same using unscaled strangeness
3216 TString strNameSecCorrPt_nonSc(Form("hSecCorrPt_nonSc_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3217 TString strNameSecCorrZ_nonSc(Form("hSecCorrZ_nonSc_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3218 TString strNameSecCorrXi_nonSc(Form("hSecCorrXi_nonSc_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3221 strNameSecCorrPt_nonSc.Form("hSecCorrBgrPt_nonSc_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
3222 strNameSecCorrZ_nonSc.Form("hSecCorrBgrZ_nonSc_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
3223 strNameSecCorrXi_nonSc.Form("hSecCorrBgrXi_nonSc_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
3226 hSecCorrPt_nonSc[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Clone(strNameSecCorrPt_nonSc);
3227 TH1F* hSumPrimSecPt_nonSc = (TH1F*) hdNdptTracksMCPrimRec[i]->Clone("hSumPrimSecPt_nonSc");
3228 hSumPrimSecPt_nonSc->Add(hdNdptTracksMCSecRecNS[i]);
3229 hSumPrimSecPt_nonSc->Add(hdNdptTracksMCSecRecS[i]); // non-scaled secondaries from strangeness
3230 hSecCorrPt_nonSc[i]->Divide(hdNdptTracksMCPrimRec[i],hSumPrimSecPt_nonSc,1,1,"B"); // binominal errors
3232 hSecCorrZ_nonSc[i] = (TH1F*) hdNdzMCPrimRec[i]->Clone(strNameSecCorrZ_nonSc);
3233 TH1F* hSumPrimSecZ_nonSc = (TH1F*) hdNdzMCPrimRec[i]->Clone("hSumPrimSecZ_nonSc");
3234 hSumPrimSecZ_nonSc->Add(hdNdzMCSecRecNS[i]);
3235 hSumPrimSecZ_nonSc->Add(hdNdzMCSecRecS[i]); // non-scaled secondaries from strangeness
3236 hSecCorrZ_nonSc[i]->Divide(hdNdzMCPrimRec[i],hSumPrimSecZ_nonSc,1,1,"B"); // binominal errors
3238 hSecCorrXi_nonSc[i] = (TH1F*) hdNdxiMCPrimRec[i]->Clone(strNameSecCorrXi_nonSc);
3239 TH1F* hSumPrimSecXi_nonSc = (TH1F*) hdNdxiMCPrimRec[i]->Clone("hSumPrimSecXi_nonSc");
3240 hSumPrimSecXi_nonSc->Add(hdNdxiMCSecRecNS[i]);
3241 hSumPrimSecXi_nonSc->Add(hdNdxiMCSecRecS[i]); // non-scaled secondaries from strangeness
3242 hSecCorrXi_nonSc[i]->Divide(hdNdxiMCPrimRec[i],hSumPrimSecXi_nonSc,1,1,"B"); // binominal errors
3247 TString outfileOption = "RECREATE";
3248 if(updateOutfile) outfileOption = "UPDATE";
3250 TFile out(strOutfile,outfileOption);
3253 Printf("%s:%d -- error opening sec corr output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
3257 if(fDebug>0) Printf("%s:%d -- write jet sec corr to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
3259 if(strOutDir && strOutDir.Length()){
3262 if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd();
3264 dir = out.mkdir(strOutDir);
3270 for(Int_t i=0; i<fNJetPtSlices; i++){
3272 hSecCorrPt[i]->Write();
3273 hSecCorrXi[i]->Write();
3274 hSecCorrZ[i]->Write();
3276 hSecCorrPt_nonSc[i]->Write();
3277 hSecCorrXi_nonSc[i]->Write();
3278 hSecCorrZ_nonSc[i]->Write();
3281 TString strSpectraDir = "spectraSecCorr";
3282 if(writeBgr) strSpectraDir = "spectraBgrSecCorr";
3284 TDirectory *dOut = gDirectory->mkdir(strSpectraDir);
3287 for(Int_t i=0; i<fNJetPtSlices; i++){
3289 hdNdptTracksMCPrimRec[i]->Write();
3290 hdNdzMCPrimRec[i]->Write();
3291 hdNdxiMCPrimRec[i]->Write();
3293 hdNdptTracksMCSecRecNS[i]->Write();
3294 hdNdzMCSecRecNS[i]->Write();
3295 hdNdxiMCSecRecNS[i]->Write();
3297 hdNdptTracksMCSecRecS[i]->Write();
3298 hdNdzMCSecRecS[i]->Write();
3299 hdNdxiMCSecRecS[i]->Write();
3301 hdNdptTracksMCSecRecSsc[i]->Write();
3302 hdNdzMCSecRecSsc[i]->Write();
3303 hdNdxiMCSecRecSsc[i]->Write();
3308 delete fh1FFJetPtRecEffRec;
3310 delete fh2FFTrackPtRecEffRec;
3311 delete fh2FFZRecEffRec;
3312 delete fh2FFXiRecEffRec;
3314 delete fh2FFTrackPtSecRecNS;
3315 delete fh2FFZSecRecNS;
3316 delete fh2FFXiSecRecNS;
3318 delete fh2FFTrackPtSecRecS;
3319 delete fh2FFZSecRecS;
3320 delete fh2FFXiSecRecS;
3322 delete fh2FFTrackPtSecRecSsc;
3323 delete fh2FFZSecRecSsc;
3324 delete fh2FFXiSecRecSsc;
3327 //________________________________________________________________________________________________________________
3328 void AliFragmentationFunctionCorrections::WriteJetResponse(TString strInfile, TString strID, TString strOutfile,
3329 Bool_t updateOutfile, TString strOutDir )
3331 // read task ouput from MC and write single track eff - standard dir/list
3333 TString strdir = "PWGJE_FragmentationFunction_" + strID;
3334 TString strlist = "fracfunc_" + strID;
3336 WriteJetResponse(strInfile,strdir,strlist,strOutfile,updateOutfile, strOutDir);
3339 //_____________________________________________________________________________________________________________________________________
3340 void AliFragmentationFunctionCorrections::WriteJetResponse(TString strInfile, TString strdir, TString strlist,
3341 TString strOutfile, Bool_t updateOutfile, TString strOutDir)
3343 // read 3d THnSparse response matrices in pt,z,xi vs jet pt from file
3344 // project THnSparse + TH2 in jet pt slices
3345 // write to strOutfile
3347 THnSparse* hn3ResponseJetPt;
3348 THnSparse* hn3ResponseJetZ;
3349 THnSparse* hn3ResponseJetXi;
3351 // 2D response matrices
3353 THnSparse* hnResponsePt[fNJetPtSlices];
3354 THnSparse* hnResponseZ[fNJetPtSlices];
3355 THnSparse* hnResponseXi[fNJetPtSlices];
3357 TH2F* h2ResponsePt[fNJetPtSlices];
3358 TH2F* h2ResponseZ[fNJetPtSlices];
3359 TH2F* h2ResponseXi[fNJetPtSlices];
3361 // 1D projections on gen pt / rec pt axes
3363 TH1F* h1FFPtRec[fNJetPtSlices];
3364 TH1F* h1FFZRec[fNJetPtSlices];
3365 TH1F* h1FFXiRec[fNJetPtSlices];
3367 TH1F* h1FFPtGen[fNJetPtSlices];
3368 TH1F* h1FFZGen[fNJetPtSlices];
3369 TH1F* h1FFXiGen[fNJetPtSlices];
3371 TH1F* h1RatioPt[fNJetPtSlices];
3372 TH1F* h1RatioZ[fNJetPtSlices];
3373 TH1F* h1RatioXi[fNJetPtSlices];
3377 TFile f(strInfile,"READ");
3380 Printf("%s:%d -- error opening file %s", (char*)__FILE__,__LINE__,strInfile.Data());
3384 if(fDebug>0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
3386 if(strdir && strdir.Length()) gDirectory->cd(strdir);
3390 if(strlist && strlist.Length()){
3392 if(!(list = (TList*) gDirectory->Get(strlist))){
3393 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
3399 hn3ResponseJetPt = (THnSparse*) list->FindObject("fhnResponseJetTrackPt");
3400 hn3ResponseJetZ = (THnSparse*) list->FindObject("fhnResponseJetZ");
3401 hn3ResponseJetXi = (THnSparse*) list->FindObject("fhnResponseJetXi");
3404 hn3ResponseJetPt = (THnSparse*) gDirectory->Get("fhnResponseJetTrackPt");
3405 hn3ResponseJetZ = (THnSparse*) gDirectory->Get("fhnResponseJetZ");
3406 hn3ResponseJetXi = (THnSparse*) gDirectory->Get("fhnResponseJetXi");
3410 if(!hn3ResponseJetPt){
3411 Printf("%s:%d -- error retrieving response matrix fhnResponseJetTrackPt",(char*)__FILE__,__LINE__);
3415 if(!hn3ResponseJetZ){
3416 Printf("%s:%d -- error retrieving response matrix fhnResponseJetZ",(char*)__FILE__,__LINE__);
3420 if(!hn3ResponseJetXi){
3421 Printf("%s:%d -- error retrieving response matrix fhnResponseJetXi",(char*)__FILE__,__LINE__);
3429 Int_t axisJetPtTHn3 = -1;
3430 Int_t axisGenPtTHn3 = -1;
3431 Int_t axisRecPtTHn3 = -1;
3433 for(Int_t i=0; i<hn3ResponseJetPt->GetNdimensions(); i++){
3435 TString title = hn3ResponseJetPt->GetAxis(i)->GetTitle();
3437 if(title.Contains("jet p_{T}")) axisJetPtTHn3 = i;
3438 if(title.Contains("gen p_{T}")) axisGenPtTHn3 = i;
3439 if(title.Contains("rec p_{T}")) axisRecPtTHn3 = i;
3442 if(axisJetPtTHn3 == -1){
3443 Printf("%s:%d -- error axisJetPtTHn3",(char*)__FILE__,__LINE__);
3447 if(axisGenPtTHn3 == -1){
3448 Printf("%s:%d -- error axisGenPtTHn3",(char*)__FILE__,__LINE__);
3452 if(axisRecPtTHn3 == -1){
3453 Printf("%s:%d -- error axisRecPtTHn3",(char*)__FILE__,__LINE__);
3457 for(Int_t i=0; i<fNJetPtSlices; i++){
3459 Float_t jetPtLoLim = fJetPtSlices->At(i);
3460 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
3462 Int_t binLo = static_cast<Int_t>(hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->FindBin(jetPtLoLim));
3463 Int_t binUp = static_cast<Int_t>(hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->FindBin(jetPtUpLim))-1;
3465 if(binUp > hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->GetNbins()){
3466 Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
3470 TString strNameRespPt(Form("hnResponsePt_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3471 TString strNameRespZ(Form("hnResponseZ_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3472 TString strNameRespXi(Form("hnResponseXi_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3474 TString strNameTH2RespPt(Form("h2ResponsePt_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3475 TString strNameTH2RespZ(Form("h2ResponseZ_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3476 TString strNameTH2RespXi(Form("h2ResponseXi_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3478 TString strNameRecPt(Form("h1FFTrackPtRecPrim_recPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3479 TString strNameRecZ(Form("h1FFZRecPrim_recPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3480 TString strNameRecXi(Form("h1FFXiRecPrim_recPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3482 TString strNameGenPt(Form("h1FFTrackPtRecPrim_genPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3483 TString strNameGenZ(Form("h1FFZRecPrim_genPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3484 TString strNameGenXi(Form("h1FFXiRecPrim_genPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3486 TString strNameRatioPt(Form("h1RatioTrackPtRecPrim_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3487 TString strNameRatioZ(Form("h1RatioZRecPrim_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3488 TString strNameRatioXi(Form("h1RatioXiRecPrim_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3491 // 2D projections in jet pt range
3493 hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->SetRange(binLo,binUp);
3494 hn3ResponseJetZ->GetAxis(axisJetPtTHn3)->SetRange(binLo,binUp);
3495 hn3ResponseJetXi->GetAxis(axisJetPtTHn3)->SetRange(binLo,binUp);
3497 Int_t axesProj[2] = {axisRecPtTHn3,axisGenPtTHn3};
3499 hnResponsePt[i] = hn3ResponseJetPt->Projection(2,axesProj);
3500 hnResponseZ[i] = hn3ResponseJetZ->Projection(2,axesProj);
3501 hnResponseXi[i] = hn3ResponseJetXi->Projection(2,axesProj);
3503 hnResponsePt[i]->SetNameTitle(strNameRespPt,"");
3504 hnResponseZ[i]->SetNameTitle(strNameRespZ,"");
3505 hnResponseXi[i]->SetNameTitle(strNameRespXi,"");
3507 h2ResponsePt[i] = (TH2F*) hnResponsePt[i]->Projection(1,0);// note convention: yDim,xDim
3508 h2ResponseZ[i] = (TH2F*) hnResponseZ[i]->Projection(1,0); // note convention: yDim,xDim
3509 h2ResponseXi[i] = (TH2F*) hnResponseXi[i]->Projection(1,0);// note convention: yDim,xDim
3511 h2ResponsePt[i]->SetNameTitle(strNameTH2RespPt,"");
3512 h2ResponseZ[i]->SetNameTitle(strNameTH2RespZ,"");
3513 h2ResponseXi[i]->SetNameTitle(strNameTH2RespXi,"");
3518 Int_t axisGenPtTHn2 = -1;
3519 Int_t axisRecPtTHn2 = -1;
3521 for(Int_t d=0; d<hnResponsePt[i]->GetNdimensions(); d++){
3523 TString title = hnResponsePt[i]->GetAxis(d)->GetTitle();
3525 if(title.Contains("gen p_{T}")) axisGenPtTHn2 = d;
3526 if(title.Contains("rec p_{T}")) axisRecPtTHn2 = d;
3530 if(axisGenPtTHn2 == -1){
3531 Printf("%s:%d -- error axisGenPtTHn2",(char*)__FILE__,__LINE__);
3535 if(axisRecPtTHn2 == -1){
3536 Printf("%s:%d -- error axisRecPtTHn2",(char*)__FILE__,__LINE__);
3541 h1FFPtRec[i] = (TH1F*) hnResponsePt[i]->Projection(axisRecPtTHn2);// note convention: yDim,xDim
3542 h1FFZRec[i] = (TH1F*) hnResponseZ[i]->Projection(axisRecPtTHn2);// note convention: yDim,xDim
3543 h1FFXiRec[i] = (TH1F*) hnResponseXi[i]->Projection(axisRecPtTHn2);// note convention: yDim,xDim
3545 h1FFPtRec[i]->SetNameTitle(strNameRecPt,"");
3546 h1FFZRec[i]->SetNameTitle(strNameRecZ,"");
3547 h1FFXiRec[i]->SetNameTitle(strNameRecXi,"");
3549 h1FFPtGen[i] = (TH1F*) hnResponsePt[i]->Projection(axisGenPtTHn2);// note convention: yDim,xDim
3550 h1FFZGen[i] = (TH1F*) hnResponseZ[i]->Projection(axisGenPtTHn2);// note convention: yDim,xDim
3551 h1FFXiGen[i] = (TH1F*) hnResponseXi[i]->Projection(axisGenPtTHn2);// note convention: yDim,xDim
3553 h1FFPtGen[i]->SetNameTitle(strNameGenPt,"");
3554 h1FFZGen[i]->SetNameTitle(strNameGenZ,"");
3555 h1FFXiGen[i]->SetNameTitle(strNameGenXi,"");
3557 // normalize 1D projections
3559 if(fNHistoBinsPt[i]) h1FFPtRec[i] = (TH1F*) h1FFPtRec[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray());
3560 if(fNHistoBinsZ[i]) h1FFZRec[i] = (TH1F*) h1FFZRec[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3561 if(fNHistoBinsXi[i]) h1FFXiRec[i] = (TH1F*) h1FFXiRec[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3563 if(fNHistoBinsPt[i]) h1FFPtGen[i] = (TH1F*) h1FFPtGen[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray());
3564 if(fNHistoBinsZ[i]) h1FFZGen[i] = (TH1F*) h1FFZGen[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3565 if(fNHistoBinsXi[i]) h1FFXiGen[i] = (TH1F*) h1FFXiGen[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3567 NormalizeTH1(h1FFPtRec[i],fNJets->At(i));
3568 NormalizeTH1(h1FFZRec[i],fNJets->At(i));
3569 NormalizeTH1(h1FFXiRec[i],fNJets->At(i));
3571 NormalizeTH1(h1FFPtGen[i],fNJets->At(i));
3572 NormalizeTH1(h1FFZGen[i],fNJets->At(i));
3573 NormalizeTH1(h1FFXiGen[i],fNJets->At(i));
3575 // ratios 1D projections
3577 h1RatioPt[i] = (TH1F*) h1FFPtRec[i]->Clone(strNameRatioPt);
3578 h1RatioPt[i]->Reset();
3579 h1RatioPt[i]->Divide(h1FFPtRec[i],h1FFPtGen[i],1,1,"B");
3581 h1RatioZ[i] = (TH1F*) h1FFZRec[i]->Clone(strNameRatioZ);
3582 h1RatioZ[i]->Reset();
3583 h1RatioZ[i]->Divide(h1FFZRec[i],h1FFZGen[i],1,1,"B");
3585 h1RatioXi[i] = (TH1F*) h1FFXiRec[i]->Clone(strNameRatioXi);
3586 h1RatioXi[i]->Reset();
3587 h1RatioXi[i]->Divide(h1FFXiRec[i],h1FFXiGen[i],1,1,"B");
3593 TString outfileOption = "RECREATE";
3594 if(updateOutfile) outfileOption = "UPDATE";
3596 TFile out(strOutfile,outfileOption);
3599 Printf("%s:%d -- error opening response matrix output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
3603 if(fDebug>0) Printf("%s:%d -- write response matrices to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
3605 if(strOutDir && strOutDir.Length()){
3608 if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd();
3610 dir = out.mkdir(strOutDir);
3616 for(Int_t i=0; i<fNJetPtSlices; i++){
3618 hnResponsePt[i]->Write();
3619 hnResponseXi[i]->Write();
3620 hnResponseZ[i]->Write();
3622 h2ResponsePt[i]->Write();
3623 h2ResponseXi[i]->Write();
3624 h2ResponseZ[i]->Write();
3626 h1FFPtRec[i]->Write();
3627 h1FFZRec[i]->Write();
3628 h1FFXiRec[i]->Write();
3630 h1FFPtGen[i]->Write();
3631 h1FFZGen[i]->Write();
3632 h1FFXiGen[i]->Write();
3634 h1RatioPt[i]->Write();
3635 h1RatioZ[i]->Write();
3636 h1RatioXi[i]->Write();
3643 //______________________________________________________________________________________________________
3644 void AliFragmentationFunctionCorrections::ReadResponse(TString strfile, TString strdir, TString strlist)
3646 // read response matrices from file
3647 // argument strlist optional - read from directory strdir if not specified
3648 // note: THnSparse are not rebinned
3650 THnSparse* hRespPt[fNJetPtSlices];
3651 THnSparse* hRespZ[fNJetPtSlices];
3652 THnSparse* hRespXi[fNJetPtSlices];
3654 for(Int_t i=0; i<fNJetPtSlices; i++) hRespPt[i] = 0;
3655 for(Int_t i=0; i<fNJetPtSlices; i++) hRespZ[i] = 0;
3656 for(Int_t i=0; i<fNJetPtSlices; i++) hRespXi[i] = 0;
3658 TFile f(strfile,"READ");
3661 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
3665 if(fDebug>0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3667 if(strdir && strdir.Length()) gDirectory->cd(strdir);
3671 if(strlist && strlist.Length()){
3673 if(!(list = (TList*) gDirectory->Get(strlist))){
3674 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
3679 for(Int_t i=0; i<fNJetPtSlices; i++){
3681 Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
3682 Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
3684 TString strNameRespPt(Form("hnResponsePt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3685 TString strNameRespZ(Form("hnResponseZ_%02d_%02d",jetPtLoLim,jetPtUpLim));
3686 TString strNameRespXi(Form("hnResponseXi_%02d_%02d",jetPtLoLim,jetPtUpLim));
3689 hRespPt[i] = (THnSparse*) list->FindObject(strNameRespPt);
3690 hRespZ[i] = (THnSparse*) list->FindObject(strNameRespZ);
3691 hRespXi[i] = (THnSparse*) list->FindObject(strNameRespXi);
3694 hRespPt[i] = (THnSparse*) gDirectory->Get(strNameRespPt);
3695 hRespZ[i] = (THnSparse*) gDirectory->Get(strNameRespZ);
3696 hRespXi[i] = (THnSparse*) gDirectory->Get(strNameRespXi);
3700 Printf("%s:%d -- error retrieving response %s", (char*)__FILE__,__LINE__,strNameRespPt.Data());
3704 Printf("%s:%d -- error retrieving response %s", (char*)__FILE__,__LINE__,strNameRespZ.Data());
3708 Printf("%s:%d -- error retrieving response %s", (char*)__FILE__,__LINE__,strNameRespXi.Data());
3711 // if(0){ // can't rebin THnSparse ...
3712 // if(fNHistoBinsPt[i]) hRespPt[i]->SetBinEdges(0,fHistoBinsPt[i]->GetArray());
3713 // if(fNHistoBinsPt[i]) hRespPt[i]->SetBinEdges(1,fHistoBinsPt[i]->GetArray());
3715 // if(fNHistoBinsZ[i]) hRespZ[i]->SetBinEdges(0,fHistoBinsZ[i]->GetArray());
3716 // if(fNHistoBinsZ[i]) hRespZ[i]->SetBinEdges(1,fHistoBinsZ[i]->GetArray());
3718 // if(fNHistoBinsXi[i]) hRespXi[i]->SetBinEdges(0,fHistoBinsXi[i]->GetArray());
3719 // if(fNHistoBinsXi[i]) hRespXi[i]->SetBinEdges(1,fHistoBinsXi[i]->GetArray());
3723 } // jet slices loop
3725 f.Close(); // THnSparse pointers still valid even if file closed
3727 // for(Int_t i=0; i<fNJetPtSlices; i++){ // no copy c'tor ...
3728 // if(hRespPt[i]) new(fhnResponsePt[i]) THnSparseF(*hRespPt[i]);
3729 // if(hRespZ[i]) new(fhnResponseZ[i]) THnSparseF(*hRespZ[i]);
3730 // if(hRespXi[i]) new(fhnResponseXi[i]) THnSparseF(*hRespXi[i]);
3733 for(Int_t i=0; i<fNJetPtSlices; i++){
3734 fhnResponsePt[i] = hRespPt[i];
3735 fhnResponseZ[i] = hRespZ[i];
3736 fhnResponseXi[i] = hRespXi[i];
3740 //______________________________________________________________________________________________________________________
3741 void AliFragmentationFunctionCorrections::ReadPriors(TString strfile,const Int_t type)
3743 // read priors from file: rec primaries, gen pt dist
3745 if(fDebug>0) Printf("%s:%d -- read priors from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3747 // temporary histos to store pointers from file
3748 TH1F* hist[fNJetPtSlices];
3750 for(Int_t i=0; i<fNJetPtSlices; i++) hist[i] = 0;
3752 TFile f(strfile,"READ");
3755 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
3759 for(Int_t i=0; i<fNJetPtSlices; i++){
3761 Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
3762 Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
3766 if(type == kFlagPt) strName.Form("h1FFTrackPtRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim);
3767 if(type == kFlagZ) strName.Form("h1FFZRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim);
3768 if(type == kFlagXi) strName.Form("h1FFXiRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim);
3770 hist[i] = (TH1F*) gDirectory->Get(strName);
3773 Printf("%s:%d -- error retrieving prior %s", (char*)__FILE__,__LINE__,strName.Data());
3777 //if(fNHistoBinsPt[i]) hist[i] = (TH1F*) hist[i]->Rebin(fNHistoBinsPt[i],hist[i]->GetName()+"_rebin",fHistoBinsPt[i]->GetArray());
3779 if(hist[i]) hist[i]->SetDirectory(0);
3781 } // jet slices loop
3786 for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
3787 if(hist[i] && type == kFlagPt) new(fh1FFTrackPtPrior[i]) TH1F(*hist[i]);
3788 if(hist[i] && type == kFlagZ) new(fh1FFZPrior[i]) TH1F(*hist[i]);
3789 if(hist[i] && type == kFlagXi) new(fh1FFXiPrior[i]) TH1F(*hist[i]);
3794 //_____________________________________________________
3795 // void AliFragmentationFunctionCorrections::RatioRecGen()
3797 // // create ratio reconstructed over generated FF
3798 // // use current highest corrLevel
3800 // Printf("%s:%d -- build ratio rec.gen, corrLevel %d",(char*)__FILE__,__LINE__,fNCorrectionLevels-1);
3802 // for(Int_t i=0; i<fNJetPtSlices; i++){
3804 // TH1F* histPtRec = fCorrFF[fNCorrectionLevels-1]->GetTrackPt(i); // levels -1: latest corr level
3805 // TH1F* histZRec = fCorrFF[fNCorrectionLevels-1]->GetZ(i); // levels -1: latest corr level
3806 // TH1F* histXiRec = fCorrFF[fNCorrectionLevels-1]->GetXi(i); // levels -1: latest corr level
3808 // TH1F* histPtGen = fh1FFTrackPtGenPrim[i];
3809 // TH1F* histZGen = fh1FFZGenPrim[i];
3810 // TH1F* histXiGen = fh1FFXiGenPrim[i];
3812 // TString histNamePt = histPtRec->GetName();
3813 // TString histNameZ = histZRec->GetName();
3814 // TString histNameXi = histXiRec->GetName();
3816 // histNamePt.ReplaceAll("fh1FF","fh1FFRatioRecGen");
3817 // histNameZ.ReplaceAll("fh1FF","fh1FFRatioRecGen");
3818 // histNameXi.ReplaceAll("fh1FF","fh1FFRatioRecGen");
3821 // TH1F* hRatioRecGenPt = (TH1F*) histPtRec->Clone(histNamePt);
3822 // hRatioRecGenPt->Reset();
3823 // hRatioRecGenPt->Divide(histPtRec,histPtGen,1,1,"B");
3825 // TH1F* hRatioRecGenZ = (TH1F*) histZRec->Clone(histNameZ);
3826 // hRatioRecGenZ->Reset();
3827 // hRatioRecGenZ->Divide(histZRec,histZGen,1,1,"B");
3829 // TH1F* hRatioRecGenXi = (TH1F*) histXiRec->Clone(histNameXi);
3830 // hRatioRecGenXi->Reset();
3831 // hRatioRecGenXi->Divide(histXiRec,histXiGen,1,1,"B");
3833 // new(fh1FFRatioRecGenPt[i]) TH1F(*hRatioRecGenPt);
3834 // new(fh1FFRatioRecGenZ[i]) TH1F(*hRatioRecGenZ);
3835 // new(fh1FFRatioRecGenXi[i]) TH1F(*hRatioRecGenXi);
3839 // //___________________________________________________________
3840 // void AliFragmentationFunctionCorrections::RatioRecPrimaries()
3842 // // create ratio reconstructed tracks over reconstructed primaries
3843 // // use raw FF (corrLevel 0)
3845 // Printf("%s:%d -- build ratio rec tracks /rec primaries",(char*)__FILE__,__LINE__);
3847 // for(Int_t i=0; i<fNJetPtSlices; i++){
3849 // const Int_t corrLevel = 0;
3851 // TH1F* histPtRec = fCorrFF[corrLevel]->GetTrackPt(i); // levels -1: latest corr level
3852 // TH1F* histZRec = fCorrFF[corrLevel]->GetZ(i); // levels -1: latest corr level
3853 // TH1F* histXiRec = fCorrFF[corrLevel]->GetXi(i); // levels -1: latest corr level
3855 // TH1F* histPtRecPrim = fh1FFTrackPtRecPrim[i];
3856 // TH1F* histZRecPrim = fh1FFZRecPrim[i];
3857 // TH1F* histXiRecPrim = fh1FFXiRecPrim[i];
3859 // TString histNamePt = histPtRec->GetName();
3860 // TString histNameZ = histZRec->GetName();
3861 // TString histNameXi = histXiRec->GetName();
3863 // histNamePt.ReplaceAll("fh1FF","fh1FFRatioRecPrim");
3864 // histNameZ.ReplaceAll("fh1FF","fh1FFRatioRecPrim");
3865 // histNameXi.ReplaceAll("fh1FF","fh1FFRatioRecPrim");
3868 // TH1F* hRatioRecPrimPt = (TH1F*) histPtRec->Clone(histNamePt);
3869 // hRatioRecPrimPt->Reset();
3870 // hRatioRecPrimPt->Divide(histPtRec,histPtRecPrim,1,1,"B");
3872 // TH1F* hRatioRecPrimZ = (TH1F*) histZRec->Clone(histNameZ);
3873 // hRatioRecPrimZ->Reset();
3874 // hRatioRecPrimZ->Divide(histZRec,histZRecPrim,1,1,"B");
3876 // TH1F* hRatioRecPrimXi = (TH1F*) histXiRec->Clone(histNameXi);
3877 // hRatioRecPrimXi->Reset();
3878 // hRatioRecPrimXi->Divide(histXiRec,histXiRecPrim,1,1,"B");
3881 // new(fh1FFRatioRecPrimPt[i]) TH1F(*hRatioRecPrimPt);
3882 // new(fh1FFRatioRecPrimZ[i]) TH1F(*hRatioRecPrimZ);
3883 // new(fh1FFRatioRecPrimXi[i]) TH1F(*hRatioRecPrimXi);
3887 // __________________________________________________________________________________
3888 void AliFragmentationFunctionCorrections::ProjectJetResponseMatrices(TString strOutfile)
3891 // project response matrices on both axes:
3892 // FF for rec primaries, in terms of generated and reconstructed momentum
3893 // write FF and ratios to outFile
3895 Printf("%s:%d -- project response matrices, write to %s",(char*)__FILE__,__LINE__,strOutfile.Data());
3897 TH1F* hFFPtRec[fNJetPtSlices];
3898 TH1F* hFFZRec[fNJetPtSlices];
3899 TH1F* hFFXiRec[fNJetPtSlices];
3901 TH1F* hFFPtGen[fNJetPtSlices];
3902 TH1F* hFFZGen[fNJetPtSlices];
3903 TH1F* hFFXiGen[fNJetPtSlices];
3905 TH1F* hRatioPt[fNJetPtSlices];
3906 TH1F* hRatioZ[fNJetPtSlices];
3907 TH1F* hRatioXi[fNJetPtSlices];
3910 Int_t axisGenPt = 1;
3911 Int_t axisRecPt = 0;
3913 for(Int_t i=0; i<fNJetPtSlices; i++){
3915 Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
3916 Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
3918 TString strNameRecPt(Form("h1FFTrackPtRecPrim_recPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3919 TString strNameRecZ(Form("h1FFZRecPrim_recPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3920 TString strNameRecXi(Form("h1FFXiRecPrim_recPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3922 TString strNameGenPt(Form("h1FFTrackPtRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3923 TString strNameGenZ(Form("h1FFZRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3924 TString strNameGenXi(Form("h1FFXiRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3926 TString strNameRatioPt(Form("h1RatioTrackPtRecPrim_%02d_%02d",jetPtLoLim,jetPtUpLim));
3927 TString strNameRatioZ(Form("h1RatioZRecPrim_%02d_%02d",jetPtLoLim,jetPtUpLim));
3928 TString strNameRatioXi(Form("h1RatioXiRecPrim_%02d_%02d",jetPtLoLim,jetPtUpLim));
3931 hFFPtRec[i] = (TH1F*) fhnResponsePt[i]->Projection(axisRecPt);// note convention: yDim,xDim
3932 hFFZRec[i] = (TH1F*) fhnResponseZ[i]->Projection(axisRecPt);// note convention: yDim,xDim
3933 hFFXiRec[i] = (TH1F*) fhnResponseXi[i]->Projection(axisRecPt);// note convention: yDim,xDim
3935 hFFPtRec[i]->SetNameTitle(strNameRecPt,"");
3936 hFFZRec[i]->SetNameTitle(strNameRecZ,"");
3937 hFFXiRec[i]->SetNameTitle(strNameRecXi,"");
3940 hFFPtGen[i] = (TH1F*) fhnResponsePt[i]->Projection(axisGenPt);// note convention: yDim,xDim
3941 hFFZGen[i] = (TH1F*) fhnResponseZ[i]->Projection(axisGenPt);// note convention: yDim,xDim
3942 hFFXiGen[i] = (TH1F*) fhnResponseXi[i]->Projection(axisGenPt);// note convention: yDim,xDim
3944 hFFPtGen[i]->SetNameTitle(strNameGenPt,"");
3945 hFFZGen[i]->SetNameTitle(strNameGenZ,"");
3946 hFFXiGen[i]->SetNameTitle(strNameGenXi,"");
3949 if(fNHistoBinsPt[i]) hFFPtRec[i] = (TH1F*) hFFPtRec[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray());
3950 if(fNHistoBinsZ[i]) hFFZRec[i] = (TH1F*) hFFZRec[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3951 if(fNHistoBinsXi[i]) hFFXiRec[i] = (TH1F*) hFFXiRec[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3953 if(fNHistoBinsPt[i]) hFFPtGen[i] = (TH1F*) hFFPtGen[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray());
3954 if(fNHistoBinsZ[i]) hFFZGen[i] = (TH1F*) hFFZGen[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3955 if(fNHistoBinsXi[i]) hFFXiGen[i] = (TH1F*) hFFXiGen[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3957 NormalizeTH1(hFFPtGen[i],fNJets->At(i));
3958 NormalizeTH1(hFFZGen[i],fNJets->At(i));
3959 NormalizeTH1(hFFXiGen[i],fNJets->At(i));
3961 NormalizeTH1(hFFPtRec[i],fNJets->At(i));
3962 NormalizeTH1(hFFZRec[i],fNJets->At(i));
3963 NormalizeTH1(hFFXiRec[i],fNJets->At(i));
3966 hRatioPt[i] = (TH1F*) hFFPtRec[i]->Clone(strNameRatioPt);
3967 hRatioPt[i]->Reset();
3968 hRatioPt[i]->Divide(hFFPtRec[i],hFFPtGen[i],1,1,"B");
3970 hRatioZ[i] = (TH1F*) hFFZRec[i]->Clone(strNameRatioZ);
3971 hRatioZ[i]->Reset();
3972 hRatioZ[i]->Divide(hFFZRec[i],hFFZGen[i],1,1,"B");
3974 hRatioXi[i] = (TH1F*) hFFXiRec[i]->Clone(strNameRatioXi);
3975 hRatioXi[i]->Reset();
3976 hRatioXi[i]->Divide(hFFXiRec[i],hFFXiGen[i],1,1,"B");
3983 TFile out(strOutfile,"RECREATE");
3986 Printf("%s:%d -- error opening response matrix output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
3990 for(Int_t i=0; i<fNJetPtSlices; i++){
3992 hFFPtRec[i]->Write();
3993 hFFZRec[i]->Write();
3994 hFFXiRec[i]->Write();
3996 hFFPtGen[i]->Write();
3997 hFFZGen[i]->Write();
3998 hFFXiGen[i]->Write();
4000 hRatioPt[i]->Write();
4001 hRatioZ[i]->Write();
4002 hRatioXi[i]->Write();
4008 // ____________________________________________________________________________________________________________________________
4009 void AliFragmentationFunctionCorrections::ProjectSingleResponseMatrix(TString strOutfile, Bool_t updateOutfile, TString strOutDir)
4011 // project response matrix on both axes:
4012 // pt spec for rec primaries, in terms of generated and reconstructed momentum
4013 // write spec and ratios to outFile
4015 Printf("%s:%d -- project single pt response matrix, write to %s",(char*)__FILE__,__LINE__,strOutfile.Data());
4021 Int_t axisGenPt = 1;
4022 Int_t axisRecPt = 0;
4024 TString strNameRecPt = "h1SpecTrackPtRecPrim_recPt";
4025 TString strNameGenPt = "h1SpecTrackPtRecPrim_genPt";
4026 TString strNameRatioPt = "h1RatioTrackPtRecPrim";
4028 hSpecPtRec = (TH1F*) fhnResponseSinglePt->Projection(axisRecPt);// note convention: yDim,xDim
4029 hSpecPtRec->SetNameTitle(strNameRecPt,"");
4031 hSpecPtGen = (TH1F*) fhnResponseSinglePt->Projection(axisGenPt);// note convention: yDim,xDim
4032 hSpecPtGen->SetNameTitle(strNameGenPt,"");
4034 hRatioPt = (TH1F*) hSpecPtRec->Clone(strNameRatioPt);
4036 hRatioPt->Divide(hSpecPtRec,hSpecPtGen,1,1,"B");
4038 TString outfileOption = "RECREATE";
4039 if(updateOutfile) outfileOption = "UPDATE";
4041 TFile out(strOutfile,outfileOption);
4044 Printf("%s:%d -- error opening reponse matrix projections output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
4049 if(strOutDir && strOutDir.Length()){
4052 if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd();
4054 dir = out.mkdir(strOutDir);
4059 hSpecPtRec->Write();
4060 hSpecPtGen->Write();
4067 //__________________________________________________________________________________________________________________________________________________________________
4068 void AliFragmentationFunctionCorrections::RebinHisto(const Int_t jetPtSlice, const Int_t nBinsLimits, Double_t* binsLimits, Double_t* binsWidth, const Int_t type)
4070 // rebin histo, rescale bins according to new width
4071 // only correct for input histos with equal bin size
4073 // args: jetPtSlice, type, use current corr level
4075 // function takes array of limits and widths (e.g. 1 GeV bins up to 10 GeV, 2 GeV width up to 20 GeV, ...)
4076 // array size of binsLimits: nBinsLimits
4077 // array size of binsWidth: nBinsLimits-1
4078 // binsLimits have to be in increasing order
4079 // if binning undefined for any slice, original binning will be kept
4082 Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
4086 if(jetPtSlice>=fNJetPtSlices){
4087 Printf("%s:%d -- jetPtSlice %d exceeds max",(char*)__FILE__,__LINE__,jetPtSlice);
4092 Double_t binLimitMin = binsLimits[0];
4093 Double_t binLimitMax = binsLimits[nBinsLimits-1];
4095 Double_t binLimit = binLimitMin; // start value
4097 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 ...
4098 TArrayD binsArray(sizeUpperLim);
4100 binsArray.SetAt(binLimitMin,nBins++);
4102 while(binLimit<binLimitMax && nBins<sizeUpperLim){
4104 Int_t currentSlice = -1;
4105 for(Int_t i=0; i<nBinsLimits; i++){
4106 if(binLimit >= 0.999*binsLimits[i]) currentSlice = i; // 0.999 numerical saftey factor
4109 Double_t currentBinWidth = binsWidth[currentSlice];
4110 binLimit += currentBinWidth;
4112 binsArray.SetAt(binLimit,nBins++);
4117 if(type == kFlagPt) hist = fCorrFF[fNCorrectionLevels-1]->GetTrackPt(jetPtSlice);
4118 else if(type == kFlagZ) hist = fCorrFF[fNCorrectionLevels-1]->GetZ(jetPtSlice);
4119 else if(type == kFlagXi) hist = fCorrFF[fNCorrectionLevels-1]->GetXi(jetPtSlice);
4120 else if(type == kFlagSinglePt) hist = fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->GetTrackPt(0);
4122 Printf("%s%d unknown type",(char*)__FILE__,__LINE__);
4126 Double_t binWidthNoRebin = hist->GetBinWidth(1);
4128 Double_t* bins = binsArray.GetArray();
4130 hist = (TH1F*) hist->Rebin(nBins-1,"",bins);
4132 for(Int_t bin=0; bin <= hist->GetNbinsX(); bin++){
4134 Double_t binWidthRebin = hist->GetBinWidth(bin);
4135 Double_t scaleF = binWidthNoRebin / binWidthRebin;
4137 Double_t binCont = hist->GetBinContent(bin);
4138 Double_t binErr = hist->GetBinError(bin);
4143 hist->SetBinContent(bin,binCont);
4144 hist->SetBinError(bin,binErr);
4149 TH1F* temp = new TH1F(*hist);
4151 if(type == kFlagPt) fCorrFF[fNCorrectionLevels-1]->ReplaceCorrHistos(jetPtSlice,temp,0,0);
4152 if(type == kFlagZ) fCorrFF[fNCorrectionLevels-1]->ReplaceCorrHistos(jetPtSlice,0,temp,0);
4153 if(type == kFlagXi) fCorrFF[fNCorrectionLevels-1]->ReplaceCorrHistos(jetPtSlice,0,0,temp);
4154 if(type == kFlagSinglePt) fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->ReplaceCorrHistos(0,temp,0,0);
4159 //__________________________________________________________________________________________________________________________________________________________________
4160 void AliFragmentationFunctionCorrections::WriteJetSpecResponse(TString strInfile, TString strdir, TString strlist/*, TString strOutfile*/)
4163 if(fDebug>0) Printf("%s:%d -- read jet spectrum response matrix from file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
4165 if(strdir && strdir.Length()) gDirectory->cd(strdir);
4169 if(strlist && strlist.Length()){
4170 if(!(list = (TList*) gDirectory->Get(strlist))){
4171 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
4175 if(list == 0)return; // catch strlist.Lenght() == 0;
4177 THnSparse* hn6ResponseJetPt = (THnSparse*) list->FindObject("fhnCorrelation");
4179 Int_t axis6RecJetPt = 0;
4180 Int_t axis6GenJetPt = 3;
4182 hn6ResponseJetPt->GetAxis(axis6RecJetPt)->SetTitle("rec jet p_{T} (GeV/c)");
4183 hn6ResponseJetPt->GetAxis(axis6GenJetPt)->SetTitle("gen jet p_{T} (GeV/c)");
4185 Int_t nBinsRecPt = hn6ResponseJetPt->GetAxis(axis6RecJetPt)->GetNbins();
4186 Double_t loLimRecPt = hn6ResponseJetPt->GetAxis(axis6RecJetPt)->GetBinLowEdge(1);
4187 Double_t upLimRecPt = hn6ResponseJetPt->GetAxis(axis6RecJetPt)->GetBinUpEdge(nBinsRecPt);
4189 Int_t nBinsGenPt = hn6ResponseJetPt->GetAxis(axis6GenJetPt)->GetNbins();
4190 Double_t loLimGenPt = hn6ResponseJetPt->GetAxis(axis6GenJetPt)->GetBinLowEdge(1);
4191 Double_t upLimGenPt = hn6ResponseJetPt->GetAxis(axis6GenJetPt)->GetBinUpEdge(nBinsGenPt);
4193 Int_t nBinsTrackPt = 200;
4194 Int_t loLimTrackPt = 0;
4195 Int_t upLimTrackPt = 200;
4198 Int_t nBinsResponse[4] = {nBinsRecPt,nBinsTrackPt,nBinsGenPt,nBinsTrackPt};
4199 Double_t binMinResponse[4] = {loLimRecPt,loLimTrackPt,loLimGenPt,loLimTrackPt};
4200 Double_t binMaxResponse[4] = {upLimRecPt,upLimTrackPt,upLimGenPt,upLimTrackPt};
4202 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)"};
4204 THnSparseD* hn4ResponseTrackPtJetPt = new THnSparseD("hn4ResponseTrackPtJetPt","",4,nBinsResponse,binMinResponse,binMaxResponse);
4206 for(Int_t i=0; i<4; i++){
4207 hn4ResponseTrackPtJetPt->GetAxis(i)->SetTitle(labelsResponseSinglePt[i]);
4216 //_____________________________________________________________________________________________________________________________________
4217 void AliFragmentationFunctionCorrections::ReadSingleTrackEfficiency(TString strfile, TString strdir, TString strlist, TString strname)
4220 ReadSingleTrackCorrection(strfile,strdir,strlist,strname,kFlagEfficiency);
4224 //_____________________________________________________________________________________________________________________________________
4225 void AliFragmentationFunctionCorrections::ReadSingleTrackResponse(TString strfile, TString strdir, TString strlist, TString strname)
4228 ReadSingleTrackCorrection(strfile,strdir,strlist,strname,kFlagResponse);
4232 //_____________________________________________________________________________________________________________________________________
4233 void AliFragmentationFunctionCorrections::ReadSingleTrackSecCorr(TString strfile, TString strdir, TString strlist, TString strname)
4236 ReadSingleTrackCorrection(strfile,strdir,strlist,strname,kFlagSecondaries);
4240 //______________________________________________________________________________________________________________________________________________________
4241 void AliFragmentationFunctionCorrections::ReadSingleTrackCorrection(TString strfile, TString strdir, TString strlist, TString strname, const Int_t type)
4243 // read single track correction (pt) from file
4244 // type: efficiency / response / secondaries correction
4246 if(!((type == kFlagEfficiency) || (type == kFlagResponse) || (type == kFlagSecondaries))){
4247 Printf("%s:%d -- no such correction ",(char*)__FILE__,__LINE__);
4251 TFile f(strfile,"READ");
4254 Printf("%s:%d -- error opening file %s", (char*)__FILE__,__LINE__,strfile.Data());
4258 if(fDebug>0 && type==kFlagEfficiency) Printf("%s:%d -- read single track corr from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
4259 if(fDebug>0 && type==kFlagResponse) Printf("%s:%d -- read single track response from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
4260 if(fDebug>0 && type==kFlagSecondaries) Printf("%s:%d -- read single track secondaries corr from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
4262 if(strdir && strdir.Length()) gDirectory->cd(strdir);
4266 if(strlist && strlist.Length()){
4268 if(!(list = (TList*) gDirectory->Get(strlist))){
4269 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
4274 TH1F* h1CorrHist = 0; // common TObject pointer not possible, need SetDirectory() later
4275 THnSparse* hnCorrHist = 0;
4277 if(type == kFlagEfficiency || type == kFlagSecondaries){
4279 if(list) h1CorrHist = (TH1F*) list->FindObject(strname);
4280 else h1CorrHist = (TH1F*) gDirectory->Get(strname);
4283 Printf("%s:%d -- error retrieving histo %s", (char*)__FILE__,__LINE__,strname.Data());
4288 else if(type == kFlagResponse){
4290 if(list) hnCorrHist = (THnSparse*) list->FindObject(strname);
4291 else hnCorrHist = (THnSparse*) gDirectory->Get(strname);
4294 Printf("%s:%d -- error retrieving histo %s", (char*)__FILE__,__LINE__,strname.Data());
4300 if(h1CorrHist) h1CorrHist->SetDirectory(0);
4301 //if(hnCorrHist) hnCorrHist->SetDirectory(0);
4305 if(type == kFlagEfficiency) fh1EffSinglePt = h1CorrHist;
4306 else if(type == kFlagResponse) fhnResponseSinglePt = hnCorrHist;
4307 else if(type == kFlagSecondaries) fh1SecCorrSinglePt = h1CorrHist;
4311 //________________________________________________________________________________________________________________
4312 void AliFragmentationFunctionCorrections::ReadRawPtSpec(TString strInfile, TString strID)
4314 // read track pt spec from task ouput - standard dir/list
4316 TString strdir = "PWGJE_FragmentationFunction_" + strID;
4317 TString strlist = "fracfunc_" + strID;
4319 ReadRawPtSpec(strInfile,strdir,strlist);
4322 //_______________________________________________________________________________________________________
4323 void AliFragmentationFunctionCorrections::ReadRawPtSpec(TString strfile, TString strdir, TString strlist)
4325 // get raw pt spectra from file
4328 fNCorrectionLevelsSinglePt = 0;
4329 fCorrSinglePt = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels];
4330 AddCorrectionLevelSinglePt(); // first 'correction' level = raw spectrum
4332 // get raw pt spec from input file, normalize
4334 TFile f(strfile,"READ");
4337 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
4341 if(fDebug>0) Printf("%s:%d -- read raw spectra from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
4343 gDirectory->cd(strdir);
4347 if(!(list = (TList*) gDirectory->Get(strlist))){
4348 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
4352 TString hnameTrackPt("fh1TrackQAPtRecCuts");
4353 TString hnameEvtSel("fh1EvtSelection");
4355 TH1F* fh1TrackPt = (TH1F*) list->FindObject(hnameTrackPt);
4356 TH1F* fh1EvtSel = (TH1F*) list->FindObject(hnameEvtSel);
4358 if(!fh1TrackPt){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameTrackPt.Data()); return; }
4359 if(!fh1EvtSel) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameEvtSel.Data()); return; }
4361 Float_t nEvents = fh1EvtSel->GetBinContent(fh1EvtSel->FindBin(0));
4363 fh1TrackPt->SetDirectory(0);
4367 // rebin + normalize
4368 if(fNHistoBinsSinglePt) fh1TrackPt = (TH1F*) fh1TrackPt->Rebin(fNHistoBinsSinglePt,"",fHistoBinsSinglePt->GetArray());
4370 NormalizeTH1(fh1TrackPt,nEvents);
4372 // raw FF = corr level 0
4373 fCorrSinglePt[0]->AddCorrHistos(0,fh1TrackPt);
4377 //_______________________________________________________________________________________________________
4378 void AliFragmentationFunctionCorrections::ReadRawPtSpecQATask(TString strfile, TString strdir, TString strlist)
4380 // get raw pt spectra from file
4381 // for output from Martas QA task
4385 fNCorrectionLevelsSinglePt = 0;
4386 fCorrSinglePt = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels];
4387 AddCorrectionLevelSinglePt(); // first 'correction' level = raw spectrum
4389 // get raw pt spec from input file, normalize
4391 TFile f(strfile,"READ");
4394 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
4398 if(fDebug>0) Printf("%s:%d -- read raw pt spec from QA task output file %s ",(char*)__FILE__,__LINE__,strfile.Data());
4400 gDirectory->cd(strdir);
4404 if(!(list = (TList*) gDirectory->Get(strlist))){
4405 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
4409 TString hnameTrackPt("fPtSel");
4410 TString hnameEvtSel("fNEventAll");
4412 TH1F* fh1TrackPt = (TH1F*) list->FindObject(hnameTrackPt);
4413 TH1F* fh1EvtSel = (TH1F*) list->FindObject(hnameEvtSel);
4415 if(!fh1TrackPt){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameTrackPt.Data()); return; }
4416 if(!fh1EvtSel) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameEvtSel.Data()); return; }
4419 // evts after physics selection
4420 Float_t nEvents = fh1EvtSel->GetEntries();
4422 fh1TrackPt->SetDirectory(0);
4427 NormalizeTH1(fh1TrackPt,nEvents);
4429 // raw FF = corr level 0
4430 fCorrSinglePt[0]->AddCorrHistos(0,fh1TrackPt);
4433 // ________________________________________________________
4434 void AliFragmentationFunctionCorrections::EffCorrSinglePt()
4436 // apply efficiency correction to inclusive track pt spec
4438 AddCorrectionLevelSinglePt("eff");
4440 TH1F* histPt = fCorrSinglePt[fNCorrectionLevelsSinglePt-2]->GetTrackPt(0);
4442 if(histPt->GetNbinsX() != fh1EffSinglePt->GetNbinsX()) Printf("%s:%d: inconsistency pt spec and eff corr bins - rebin effCorr ...", (char*)__FILE__,__LINE__);
4444 TString histNamePt = histPt->GetName();
4445 TH1F* hTrackPtEffCorr = (TH1F*) histPt->Clone(histNamePt);
4446 hTrackPtEffCorr->Divide(histPt,fh1EffSinglePt,1,1,"");
4448 fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->AddCorrHistos(0,hTrackPtEffCorr);
4451 //___________________________________________________________________________________________________________________________
4452 void AliFragmentationFunctionCorrections::UnfoldSinglePt(const Int_t nIter, const Bool_t useCorrelatedErrors)
4454 // unfolde inclusive dN/dpt spectra
4456 AddCorrectionLevelSinglePt("unfold");
4458 TH1F* hist = fCorrSinglePt[fNCorrectionLevelsSinglePt-2]->GetTrackPt(0); // level -2: before unfolding, level -1: unfolded
4459 THnSparse* hnResponse = fhnResponseSinglePt;
4461 TString histNameTHn = hist->GetName();
4462 if(histNameTHn.Contains("TH1")) histNameTHn.ReplaceAll("TH1","THn");
4463 if(histNameTHn.Contains("fPt")) histNameTHn.ReplaceAll("fPt","fhnPt");
4466 TString histNameBackFolded = hist->GetName();
4467 histNameBackFolded.Append("_backfold");
4469 TString histNameRatioFolded = hist->GetName();
4470 if(histNameRatioFolded.Contains("fh1")) histNameRatioFolded.ReplaceAll("fh1","hRatio");
4471 if(histNameRatioFolded.Contains("fPt")) histNameRatioFolded.ReplaceAll("fPt","hRatioPt");
4472 histNameRatioFolded.Append("_unfold");
4474 TString histNameRatioBackFolded = hist->GetName();
4475 if(histNameRatioBackFolded.Contains("fh1")) histNameRatioBackFolded.ReplaceAll("fh1","hRatio");
4476 if(histNameRatioBackFolded.Contains("fPt")) histNameRatioBackFolded.ReplaceAll("fPt","hRatioPt");
4477 histNameRatioBackFolded.Append("_backfold");
4479 THnSparse* hnHist = TH1toSparse(hist,histNameTHn,hist->GetTitle());
4480 THnSparse* hnFlatEfficiency = TH1toSparse(hist,"fhnEfficiency","eff",kTRUE); // could optionally also use real eff
4483 = Unfold(hnHist,hnResponse,hnFlatEfficiency,nIter,useCorrelatedErrors);
4485 //TH1F* hUnfolded = (TH1F*) hnUnfolded->Projection(0);
4486 //hUnfolded->SetNameTitle(hist->GetName(),hist->GetTitle());
4488 fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->AddCorrHistos(0,hUnfolded);
4490 // backfolding: apply response matrix to unfolded spectrum
4491 TH1F* hBackFolded = ApplyResponse(hUnfolded,hnResponse);
4492 hBackFolded->SetNameTitle(histNameBackFolded,hist->GetTitle());
4494 fh1SingleTrackPtBackFolded = hBackFolded;
4497 // ratio unfolded to original histo
4498 TH1F* hRatioUnfolded = (TH1F*) hUnfolded->Clone(histNameRatioFolded);
4499 hRatioUnfolded->Reset();
4500 hRatioUnfolded->Divide(hUnfolded,hist,1,1,"B");
4502 fh1RatioSingleTrackPtFolded = hRatioUnfolded;
4505 // ratio backfolded to original histo
4506 TH1F* hRatioBackFolded = (TH1F*) hBackFolded->Clone(histNameRatioBackFolded);
4507 hRatioBackFolded->Reset();
4508 hRatioBackFolded->Divide(hBackFolded,hist,1,1,"B");
4510 fh1RatioSingleTrackPtBackFolded = hRatioBackFolded;
4513 delete hnFlatEfficiency;
4517 // ________________________________________________________
4518 void AliFragmentationFunctionCorrections::SecCorrSinglePt()
4520 // apply efficiency correction to inclusive track pt spec
4522 AddCorrectionLevelSinglePt("secCorr");
4524 TH1F* histPt = fCorrSinglePt[fNCorrectionLevelsSinglePt-2]->GetTrackPt(0);
4526 if(histPt->GetNbinsX() != fh1SecCorrSinglePt->GetNbinsX())
4527 Printf("%s:%d: inconsistency pt spec and secondaries corr bins - rebin effCorr ...", (char*)__FILE__,__LINE__);
4529 TString histNamePt = histPt->GetName();
4530 TH1F* hTrackPtSecCorr = (TH1F*) histPt->Clone(histNamePt);
4532 hTrackPtSecCorr->Multiply(histPt,fh1SecCorrSinglePt,1,1,"");
4534 fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->AddCorrHistos(0,hTrackPtSecCorr);
4537 //___________________________________________________________________________________________________________________________________
4538 void AliFragmentationFunctionCorrections::dNdz2dNdxi()
4540 // transform dN/dz distribution into dN/dxi
4541 // for current corr level, all jet pt slices
4544 for(Int_t i=0; i<fNJetPtSlices; i++){
4546 TH1F* histZ = fCorrFF[fNCorrectionLevels-1]->GetZ(i);
4547 Int_t nBins = histZ->GetNbinsX();
4549 Double_t* binLims = new Double_t[nBins+1];
4551 for(Int_t bin = 0; bin<nBins; bin++){
4553 Int_t binZ = nBins-bin;
4555 Double_t zLo = histZ->GetXaxis()->GetBinLowEdge(binZ);
4556 Double_t zUp = histZ->GetXaxis()->GetBinUpEdge(binZ);
4558 Double_t xiLo = TMath::Log(1/zUp);
4559 Double_t xiUp = TMath::Log(1/zLo);
4561 if(bin == 0) binLims[0] = xiLo;
4562 binLims[bin+1] = xiUp;
4565 // for(Int_t bin = 0; bin<=nBins; bin++) std::cout<<" bin "<<bin<<" binLims "<<binLims[bin]<<std::endl;
4567 TString strTitle = histZ->GetTitle();
4568 TString strName = histZ->GetName();
4570 strName.ReplaceAll("Z","XiNew");
4571 strTitle.ReplaceAll("Z","XiNew");
4574 TH1F* histXiNew = new TH1F("histXiNew","",nBins,binLims);
4575 histXiNew->SetNameTitle(strName,strTitle);
4578 for(Int_t binZ = 1; binZ<=nBins; binZ++){
4580 Double_t meanZ = histZ->GetBinCenter(binZ);
4581 Double_t cont = histZ->GetBinContent(binZ);
4582 Double_t err = histZ->GetBinError(binZ);
4584 Double_t meanXi = TMath::Log(1/meanZ);
4585 Int_t binXi = histXiNew->FindBin(meanXi);
4587 histXiNew->SetBinContent(binXi,cont);
4588 histXiNew->SetBinError(binXi,err);
4590 //std::cout<<" binZ "<<binZ<<" meanZ "<<meanZ<<" binXi "<<binXi<<" meanXi "<<meanXi<<" cont "<<cont<<" err "<<err<<std::endl;
4593 fCorrFF[fNCorrectionLevels-1]->ReplaceCorrHistos(i,0,0,histXiNew);
4600 //________________________________________________________________________________________________________________
4601 void AliFragmentationFunctionCorrections::WriteBinShiftCorr(TString strInfile, TString strIDGen, TString strIDRec,
4602 TString strOutfile, Bool_t updateOutfile, Bool_t useRecPrim,
4603 TString strOutDir,Double_t scaleFacBgrRec)
4605 TString strdirGen = "PWGJE_FragmentationFunction_" + strIDGen;
4606 TString strlistGen = "fracfunc_" + strIDGen;
4608 TString strdirRec = "PWGJE_FragmentationFunction_" + strIDRec;
4609 TString strlistRec = "fracfunc_" + strIDRec;
4611 WriteBinShiftCorr(strInfile,strdirGen,strlistGen,strdirRec,strlistRec,strOutfile,updateOutfile,useRecPrim, strOutDir, kFALSE, "",scaleFacBgrRec);
4614 //________________________________________________________________________________________________________________
4615 void AliFragmentationFunctionCorrections::WriteBgrBinShiftCorr(TString strInfile, TString strBgrID, TString strIDGen, TString strIDRec,
4616 TString strOutfile, Bool_t updateOutfile, Bool_t useRecPrim,
4617 TString strOutDir,Double_t scaleFacBgrRec)
4619 TString strdirGen = "PWGJE_FragmentationFunction_" + strIDGen;
4620 TString strlistGen = "fracfunc_" + strIDGen;
4622 TString strdirRec = "PWGJE_FragmentationFunction_" + strIDRec;
4623 TString strlistRec = "fracfunc_" + strIDRec;
4625 WriteBinShiftCorr(strInfile,strdirGen,strlistGen,strdirRec,strlistRec,strOutfile,updateOutfile,useRecPrim, strOutDir, kTRUE, strBgrID,scaleFacBgrRec);
4628 //___________________________________________________________________________________________________________________________________
4629 void AliFragmentationFunctionCorrections::WriteBinShiftCorr(TString strInfile, TString strdirGen, TString strlistGen,
4630 TString strdirRec, TString strlistRec,
4631 TString strOutfile, Bool_t updateOutfile, Bool_t useRecPrim,
4632 TString strOutDir, Bool_t writeBgr, TString strBgrID, Double_t scaleFacBgrRec)
4635 if((writeBgr && strBgrID.Length() == 0) || (!writeBgr && strBgrID.Length()>0) ){
4636 Printf("%s:%d -- inconsistent arguments to WriteBinShiftCorr FF/UE", (char*)__FILE__,__LINE__);
4640 TH1F* hCorrPt[fNJetPtSlices];
4641 TH1F* hCorrXi[fNJetPtSlices];
4642 TH1F* hCorrZ[fNJetPtSlices];
4644 TH1F* hdNdptTracksMCGen[fNJetPtSlices];
4645 TH1F* hdNdxiMCGen[fNJetPtSlices];
4646 TH1F* hdNdzMCGen[fNJetPtSlices];
4648 TH1F* hdNdptTracksMCRec[fNJetPtSlices];
4649 TH1F* hdNdxiMCRec[fNJetPtSlices];
4650 TH1F* hdNdzMCRec[fNJetPtSlices];
4652 TH1F* fh1FFJetPtMCGen = 0;
4653 TH1F* fh1FFJetPtMCRec = 0;
4655 TH2F* fh2FFTrackPtMCGen = 0;
4656 TH2F* fh2FFZMCGen = 0;
4657 TH2F* fh2FFXiMCGen = 0;
4659 TH2F* fh2FFTrackPtMCRec = 0;
4660 TH2F* fh2FFZMCRec = 0;
4661 TH2F* fh2FFXiMCRec = 0;
4665 TFile f(strInfile,"READ");
4668 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
4672 if(fDebug>0) Printf("%s:%d -- writeBinShiftCorr: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
4674 if(strdirGen && strdirGen.Length()) gDirectory->cd(strdirGen);
4678 if(strlistGen && strlistGen.Length()){
4680 if(!(listGen = (TList*) gDirectory->Get(strlistGen))){
4681 Printf("%s:%d -- error retrieving listGen %s from directory %s", (char*)__FILE__,__LINE__,strlistGen.Data(),strdirGen.Data());
4688 fh1FFJetPtMCGen = (TH1F*) listGen->FindObject("fh1FFJetPtGen");
4691 fh2FFTrackPtMCGen = (TH2F*) listGen->FindObject(Form("fh2FFTrackPt%sGen",strBgrID.Data()));
4692 fh2FFZMCGen = (TH2F*) listGen->FindObject(Form("fh2FFZ%sGen",strBgrID.Data()));
4693 fh2FFXiMCGen = (TH2F*) listGen->FindObject(Form("fh2FFXi%sGen",strBgrID.Data()));
4696 fh2FFTrackPtMCGen = (TH2F*) listGen->FindObject("fh2FFTrackPtGen");
4697 fh2FFZMCGen = (TH2F*) listGen->FindObject("fh2FFZGen");
4698 fh2FFXiMCGen = (TH2F*) listGen->FindObject("fh2FFXiGen");
4702 fh1FFJetPtMCGen = (TH1F*) gDirectory->Get("fh1FFJetPtGen");
4704 fh2FFTrackPtMCGen = (TH2F*) gDirectory->Get(Form("fh2FFTrackPt%sGen",strBgrID.Data()));
4705 fh2FFZMCGen = (TH2F*) gDirectory->Get(Form("fh2FFZ%sGen",strBgrID.Data()));
4706 fh2FFXiMCGen = (TH2F*) gDirectory->Get(Form("fh2FFXi%sGen",strBgrID.Data()));
4709 fh2FFTrackPtMCGen = (TH2F*) gDirectory->Get("fh2FFTrackPtGen");
4710 fh2FFZMCGen = (TH2F*) gDirectory->Get("fh2FFZGen");
4711 fh2FFXiMCGen = (TH2F*) gDirectory->Get("fh2FFXiGen");
4715 fh1FFJetPtMCGen->SetDirectory(0);
4716 fh2FFTrackPtMCGen->SetDirectory(0);
4717 fh2FFZMCGen->SetDirectory(0);
4718 fh2FFXiMCGen->SetDirectory(0);
4724 TFile g(strInfile,"READ");
4727 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
4731 if(fDebug>0) Printf("%s:%d -- writeBinShiftCorr: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
4733 if(strdirRec && strdirRec.Length()) gDirectory->cd(strdirRec);
4737 if(strlistRec && strlistRec.Length()){
4739 if(!(listRec = (TList*) gDirectory->Get(strlistRec))){
4740 Printf("%s:%d -- error retrieving listRec %s from directory %s", (char*)__FILE__,__LINE__,strlistRec.Data(),strdirRec.Data());
4747 fh1FFJetPtMCRec = (TH1F*) listRec->FindObject("fh1FFJetPtRecEffRec");
4750 fh2FFTrackPtMCRec = (TH2F*) listRec->FindObject(Form("fh2FFTrackPt%sRecEffRec",strBgrID.Data()));
4751 fh2FFZMCRec = (TH2F*) listRec->FindObject(Form("fh2FFZ%sRecEffRec",strBgrID.Data()));
4752 fh2FFXiMCRec = (TH2F*) listRec->FindObject(Form("fh2FFXi%sRecEffRec",strBgrID.Data()));
4755 fh2FFTrackPtMCRec = (TH2F*) listRec->FindObject("fh2FFTrackPtRecEffRec");
4756 fh2FFZMCRec = (TH2F*) listRec->FindObject("fh2FFZRecEffRec");
4757 fh2FFXiMCRec = (TH2F*) listRec->FindObject("fh2FFXiRecEffRec");
4761 fh1FFJetPtMCRec = (TH1F*) gDirectory->Get("fh1FFJetPtRecEffRec");
4763 fh2FFTrackPtMCRec = (TH2F*) gDirectory->Get(Form("fh2FFTrackPt%sRecEffRec",strBgrID.Data()));
4764 fh2FFZMCRec = (TH2F*) gDirectory->Get(Form("fh2FFZ%sRecEffRec",strBgrID.Data()));
4765 fh2FFXiMCRec = (TH2F*) gDirectory->Get(Form("fh2FFXi%sRecEffRec",strBgrID.Data()));
4768 fh2FFTrackPtMCRec = (TH2F*) gDirectory->Get("fh2FFTrackPtRecEffRec");
4769 fh2FFZMCRec = (TH2F*) gDirectory->Get("fh2FFZRecEffRec");
4770 fh2FFXiMCRec = (TH2F*) gDirectory->Get("fh2FFXiRecEffRec");
4776 fh1FFJetPtMCRec = (TH1F*) listRec->FindObject("fh1FFJetPtRecCuts");
4778 fh2FFTrackPtMCRec = (TH2F*) listRec->FindObject(Form("fh2FFTrackPt%sRecCuts",strBgrID.Data()));
4779 fh2FFZMCRec = (TH2F*) listRec->FindObject(Form("fh2FFZ%sRecCuts",strBgrID.Data()));
4780 fh2FFXiMCRec = (TH2F*) listRec->FindObject(Form("fh2FFXi%sRecCuts",strBgrID.Data()));
4783 fh2FFTrackPtMCRec = (TH2F*) listRec->FindObject("fh2FFTrackPtRecCuts");
4784 fh2FFZMCRec = (TH2F*) listRec->FindObject("fh2FFZRecCuts");
4785 fh2FFXiMCRec = (TH2F*) listRec->FindObject("fh2FFXiRecCuts");
4789 fh1FFJetPtMCRec = (TH1F*) gDirectory->Get("fh1FFJetPtRecCuts");
4791 fh2FFTrackPtMCRec = (TH2F*) gDirectory->Get(Form("fh2FFTrackPt%sRecCuts",strBgrID.Data()));
4792 fh2FFZMCRec = (TH2F*) gDirectory->Get(Form("fh2FFZ%sRecCuts",strBgrID.Data()));
4793 fh2FFXiMCRec = (TH2F*) gDirectory->Get(Form("fh2FFXi%sRecCuts",strBgrID.Data()));
4796 fh2FFTrackPtMCRec = (TH2F*) gDirectory->Get("fh2FFTrackPtRecCuts");
4797 fh2FFZMCRec = (TH2F*) gDirectory->Get("fh2FFZRecCuts");
4798 fh2FFXiMCRec = (TH2F*) gDirectory->Get("fh2FFXiRecCuts");
4804 fh1FFJetPtMCRec->SetDirectory(0);
4805 fh2FFTrackPtMCRec->SetDirectory(0);
4806 fh2FFZMCRec->SetDirectory(0);
4807 fh2FFXiMCRec->SetDirectory(0);
4811 // projections: FF for generated and reconstructed
4813 for(Int_t i=0; i<fNJetPtSlices; i++){
4815 Float_t jetPtLoLim = fJetPtSlices->At(i);
4816 Float_t jetPtUpLim = fJetPtSlices->At(i+1);
4818 Int_t binLo = static_cast<Int_t>(fh2FFTrackPtMCGen->GetXaxis()->FindBin(jetPtLoLim));
4819 Int_t binUp = static_cast<Int_t>(fh2FFTrackPtMCGen->GetXaxis()->FindBin(jetPtUpLim))-1;
4821 if(binUp > fh2FFTrackPtMCGen->GetNbinsX()){
4822 Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim);
4826 TString strNameFFPtGen(Form("fh1FFTrackPtGenBinShift_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
4827 TString strNameFFZGen(Form("fh1FFZGenBinShift_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
4828 TString strNameFFXiGen(Form("fh1FFXiGenBinShift_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
4830 TString strNameFFPtRec(Form("fh1FFTrackPtRecBinShift_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
4831 TString strNameFFZRec(Form("fh1FFZRecBinShift_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
4832 TString strNameFFXiRec(Form("fh1FFXiRecBinShift_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
4835 strNameFFPtGen.Form("fh1TrackPtGenBinShiftBgr_%02d_%02d_%s",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim),strBgrID.Data());
4836 strNameFFZGen.Form("fh1ZGenBinShiftBgr_%02d_%02d_%s",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim),strBgrID.Data());
4837 strNameFFXiGen.Form("fh1XiGenBinShiftBgr_%02d_%02d_%s",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim),strBgrID.Data());
4839 strNameFFPtRec.Form("fh1TrackPtRecBinShiftBgr_%02d_%02d_%s",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim),strBgrID.Data());
4840 strNameFFZRec.Form("fh1ZRecBinShiftBgr_%02d_%02d_%s",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim),strBgrID.Data());
4841 strNameFFXiRec.Form("fh1XiRecBinShiftBgr_%02d_%02d_%s",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim),strBgrID.Data());
4845 // appendix 'unbinned' to avoid histos with same name after rebinning
4847 hdNdptTracksMCGen[i] = (TH1F*) fh2FFTrackPtMCGen->ProjectionY(strNameFFPtGen+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
4848 hdNdzMCGen[i] = (TH1F*) fh2FFZMCGen->ProjectionY(strNameFFZGen+"_unBinned",binLo,binUp,"o");
4849 hdNdxiMCGen[i] = (TH1F*) fh2FFXiMCGen->ProjectionY(strNameFFXiGen+"_unBinned",binLo,binUp,"o");
4851 hdNdptTracksMCRec[i] = (TH1F*) fh2FFTrackPtMCRec->ProjectionY(strNameFFPtRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range
4852 hdNdzMCRec[i] = (TH1F*) fh2FFZMCRec->ProjectionY(strNameFFZRec+"_unBinned",binLo,binUp,"o");
4853 hdNdxiMCRec[i] = (TH1F*) fh2FFXiMCRec->ProjectionY(strNameFFXiRec+"_unBinned",binLo,binUp,"o");
4858 if(fNHistoBinsPt[i]) hdNdptTracksMCGen[i] = (TH1F*) hdNdptTracksMCGen[i]->Rebin(fNHistoBinsPt[i],strNameFFPtGen,fHistoBinsPt[i]->GetArray());
4859 if(fNHistoBinsZ[i]) hdNdzMCGen[i] = (TH1F*) hdNdzMCGen[i]->Rebin(fNHistoBinsZ[i],strNameFFZGen,fHistoBinsZ[i]->GetArray());
4860 if(fNHistoBinsXi[i]) hdNdxiMCGen[i] = (TH1F*) hdNdxiMCGen[i]->Rebin(fNHistoBinsXi[i],strNameFFXiGen,fHistoBinsXi[i]->GetArray());
4862 if(fNHistoBinsPt[i]) hdNdptTracksMCRec[i] = (TH1F*) hdNdptTracksMCRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtRec,fHistoBinsPt[i]->GetArray());
4863 if(fNHistoBinsZ[i]) hdNdzMCRec[i] = (TH1F*) hdNdzMCRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZRec,fHistoBinsZ[i]->GetArray());
4864 if(fNHistoBinsXi[i]) hdNdxiMCRec[i] = (TH1F*) hdNdxiMCRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiRec,fHistoBinsXi[i]->GetArray());
4867 hdNdptTracksMCGen[i]->SetNameTitle(strNameFFPtGen,"");
4868 hdNdzMCGen[i]->SetNameTitle(strNameFFZGen,"");
4869 hdNdxiMCGen[i]->SetNameTitle(strNameFFXiGen,"");
4871 hdNdptTracksMCRec[i]->SetNameTitle(strNameFFPtRec,"");
4872 hdNdzMCRec[i]->SetNameTitle(strNameFFZRec,"");
4873 hdNdxiMCRec[i]->SetNameTitle(strNameFFXiRec,"");
4877 Double_t nJetsBinGen = fh1FFJetPtMCGen->Integral(binLo,binUp);
4878 Double_t nJetsBinRec = fh1FFJetPtMCRec->Integral(binLo,binUp);
4880 // scale fac for perp2 bgr
4881 if(useRecPrim && writeBgr && scaleFacBgrRec && (scaleFacBgrRec != 1)) nJetsBinRec /= scaleFacBgrRec;
4883 NormalizeTH1(hdNdptTracksMCGen[i],nJetsBinGen);
4884 NormalizeTH1(hdNdzMCGen[i],nJetsBinGen);
4885 NormalizeTH1(hdNdxiMCGen[i],nJetsBinGen);
4887 NormalizeTH1(hdNdptTracksMCRec[i],nJetsBinRec);
4888 NormalizeTH1(hdNdzMCRec[i],nJetsBinRec);
4889 NormalizeTH1(hdNdxiMCRec[i],nJetsBinRec);
4891 // divide gen/rec : corr factor
4893 TString strNameCorrPt(Form("hBbBCorrPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
4894 TString strNameCorrZ(Form("hBbBCorrZ_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
4895 TString strNameCorrXi(Form("hBbBCorrXi_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
4898 strNameCorrPt.Form("hBbBCorrBgrPt_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
4899 strNameCorrZ.Form("hBbBCorrBgrZ_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
4900 strNameCorrXi.Form("hBbBCorrBgrXi_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
4903 hCorrPt[i] = (TH1F*) hdNdptTracksMCGen[i]->Clone(strNameCorrPt);
4904 hCorrPt[i]->Divide(hdNdptTracksMCGen[i],hdNdptTracksMCRec[i],1,1,"B"); // binominal errors
4906 hCorrXi[i] = (TH1F*) hdNdxiMCGen[i]->Clone(strNameCorrXi);
4907 hCorrXi[i]->Divide(hdNdxiMCGen[i],hdNdxiMCRec[i],1,1,"B"); // binominal errors
4909 hCorrZ[i] = (TH1F*) hdNdzMCGen[i]->Clone(strNameCorrZ);
4910 hCorrZ[i]->Divide(hdNdzMCGen[i],hdNdzMCRec[i],1,1,"B"); // binominal errors
4915 TString outfileOption = "RECREATE";
4916 if(updateOutfile) outfileOption = "UPDATE";
4918 TFile out(strOutfile,outfileOption);
4921 Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
4925 if(fDebug>0) Printf("%s:%d -- write bin shift correction to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
4927 if(strOutDir && strOutDir.Length()){
4930 if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd();
4932 dir = out.mkdir(strOutDir);
4937 for(Int_t i=0; i<fNJetPtSlices; i++){
4939 hdNdptTracksMCGen[i]->Write();
4940 hdNdxiMCGen[i]->Write();
4941 hdNdzMCGen[i]->Write();
4943 hdNdptTracksMCRec[i]->Write();
4944 hdNdxiMCRec[i]->Write();
4945 hdNdzMCRec[i]->Write();
4947 hCorrPt[i]->Write();
4948 hCorrXi[i]->Write();
4954 delete fh1FFJetPtMCGen;
4955 delete fh1FFJetPtMCRec;
4957 delete fh2FFTrackPtMCGen;
4959 delete fh2FFXiMCGen;
4961 delete fh2FFTrackPtMCRec;
4963 delete fh2FFXiMCRec;
4966 //________________________________________________________________________________________________________________________________
4967 void AliFragmentationFunctionCorrections::ReadBgrBinShiftCorr(TString strfile, TString strBgrID, TString strdir, TString strlist)
4970 ReadBinShiftCorr(strfile, strdir, strlist, kTRUE, strBgrID);
4973 //___________________________________________________________________________________________________________________________________________
4974 void AliFragmentationFunctionCorrections::ReadBinShiftCorr(TString strfile, TString strdir, TString strlist, Bool_t readBgr, TString strBgrID)
4977 if((readBgr && strBgrID.Length() == 0) || (!readBgr && strBgrID.Length()>0) ){
4978 Printf("%s:%d -- inconsistent arguments to ReadBinShiftCorr FF/UE", (char*)__FILE__,__LINE__);
4982 // temporary histos to hold histos from file
4983 TH1F* hCorrPt[fNJetPtSlices];
4984 TH1F* hCorrZ[fNJetPtSlices];
4985 TH1F* hCorrXi[fNJetPtSlices];
4987 for(Int_t i=0; i<fNJetPtSlices; i++) hCorrPt[i] = 0;
4988 for(Int_t i=0; i<fNJetPtSlices; i++) hCorrZ[i] = 0;
4989 for(Int_t i=0; i<fNJetPtSlices; i++) hCorrXi[i] = 0;
4991 TFile f(strfile,"READ");
4994 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
4998 if(fDebug>0) Printf("%s:%d -- read FF / UE bin shift correction from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
5000 if(strdir && strdir.Length()) gDirectory->cd(strdir);
5004 if(strlist && strlist.Length()){
5006 if(!(list = (TList*) gDirectory->Get(strlist))){
5007 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
5012 for(Int_t i=0; i<fNJetPtSlices; i++){
5014 Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
5015 Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
5017 TString strNameCorrPt(Form("hBbBCorrPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
5018 TString strNameCorrZ(Form("hBbBCorrZ_%02d_%02d",jetPtLoLim,jetPtUpLim));
5019 TString strNameCorrXi(Form("hBbBCorrXi_%02d_%02d",jetPtLoLim,jetPtUpLim));
5022 strNameCorrPt.Form("hBbBCorrBgrPt_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
5023 strNameCorrZ.Form("hBbBCorrBgrZ_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
5024 strNameCorrXi.Form("hBbBCorrBgrXi_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
5029 hCorrPt[i] = (TH1F*) list->FindObject(strNameCorrPt);
5030 hCorrZ[i] = (TH1F*) list->FindObject(strNameCorrZ);
5031 hCorrXi[i] = (TH1F*) list->FindObject(strNameCorrXi);
5034 hCorrPt[i] = (TH1F*) gDirectory->Get(strNameCorrPt);
5035 hCorrZ[i] = (TH1F*) gDirectory->Get(strNameCorrZ);
5036 hCorrXi[i] = (TH1F*) gDirectory->Get(strNameCorrXi);
5040 Printf("%s:%d -- error retrieving bin by bin correction %s", (char*)__FILE__,__LINE__,strNameCorrPt.Data());
5044 Printf("%s:%d -- error retrieving bin by bin correction %s", (char*)__FILE__,__LINE__,strNameCorrZ.Data());
5048 Printf("%s:%d -- error retrieving bin by bin correction %s", (char*)__FILE__,__LINE__,strNameCorrXi.Data());
5052 if(fNHistoBinsPt[i]) hCorrPt[i] = (TH1F*) hCorrPt[i]->Rebin(fNHistoBinsPt[i],strNameCorrPt+"_rebin",fHistoBinsPt[i]->GetArray());
5053 if(fNHistoBinsZ[i]) hCorrZ[i] = (TH1F*) hCorrZ[i]->Rebin(fNHistoBinsZ[i],strNameCorrZ+"_rebin",fHistoBinsZ[i]->GetArray());
5054 if(fNHistoBinsXi[i]) hCorrXi[i] = (TH1F*) hCorrXi[i]->Rebin(fNHistoBinsXi[i],strNameCorrXi+"_rebin",fHistoBinsXi[i]->GetArray());
5056 if(hCorrPt[i]) hCorrPt[i]->SetDirectory(0);
5057 if(hCorrZ[i]) hCorrZ[i]->SetDirectory(0);
5058 if(hCorrXi[i]) hCorrXi[i]->SetDirectory(0);
5060 } // jet slices loop
5065 for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
5066 if(hCorrPt[i]) new(fh1BbBBgrPt[i]) TH1F(*hCorrPt[i]);
5067 if(hCorrZ[i]) new(fh1BbBBgrZ[i]) TH1F(*hCorrZ[i]);
5068 if(hCorrXi[i]) new(fh1BbBBgrXi[i]) TH1F(*hCorrXi[i]);
5072 for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
5073 if(hCorrPt[i]) new(fh1BbBPt[i]) TH1F(*hCorrPt[i]);
5074 if(hCorrZ[i]) new(fh1BbBZ[i]) TH1F(*hCorrZ[i]);
5075 if(hCorrXi[i]) new(fh1BbBXi[i]) TH1F(*hCorrXi[i]);
5080 // ________________________________________________
5081 void AliFragmentationFunctionCorrections::BbBCorr()
5083 // apply bin-by-bin correction
5085 AddCorrectionLevel("BbB");
5087 for(Int_t i=0; i<fNJetPtSlices; i++){
5089 TH1F* histPt = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i);
5090 TH1F* histZ = fCorrFF[fNCorrectionLevels-2]->GetZ(i);
5091 TH1F* histXi = fCorrFF[fNCorrectionLevels-2]->GetXi(i);
5093 TString histNamePt = histPt->GetName();
5094 TString histNameZ = histZ->GetName();
5095 TString histNameXi = histXi->GetName();
5097 TH1F* hFFTrackPtBbBCorr = (TH1F*) histPt->Clone(histNamePt);
5098 hFFTrackPtBbBCorr->Multiply(histPt,fh1BbBPt[i],1,1,"");
5100 TH1F* hFFZBbBCorr = (TH1F*) histZ->Clone(histNameZ);
5101 hFFZBbBCorr->Multiply(histZ,fh1BbBZ[i],1,1,"");
5103 TH1F* hFFXiBbBCorr = (TH1F*) histXi->Clone(histNameXi);
5104 hFFXiBbBCorr->Multiply(histXi,fh1BbBXi[i],1,1,"");
5106 fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hFFTrackPtBbBCorr,hFFZBbBCorr,hFFXiBbBCorr);
5110 // ___________________________________________________
5111 void AliFragmentationFunctionCorrections::BbBCorrBgr()
5113 // apply bin-by-bin correction
5115 AddCorrectionLevelBgr("BbB");
5117 for(Int_t i=0; i<fNJetPtSlices; i++){
5119 TH1F* histPt = fCorrBgr[fNCorrectionLevelsBgr-2]->GetTrackPt(i);
5120 TH1F* histZ = fCorrBgr[fNCorrectionLevelsBgr-2]->GetZ(i);
5121 TH1F* histXi = fCorrBgr[fNCorrectionLevelsBgr-2]->GetXi(i);
5123 TString histNamePt = histPt->GetName();
5124 TString histNameZ = histZ->GetName();
5125 TString histNameXi = histXi->GetName();
5127 TH1F* hBgrTrackPtBbBCorr = (TH1F*) histPt->Clone(histNamePt);
5128 hBgrTrackPtBbBCorr->Multiply(histPt,fh1BbBBgrPt[i],1,1,"");
5130 TH1F* hBgrZBbBCorr = (TH1F*) histZ->Clone(histNameZ);
5131 hBgrZBbBCorr->Multiply(histZ,fh1BbBBgrZ[i],1,1,"");
5133 TH1F* hBgrXiBbBCorr = (TH1F*) histXi->Clone(histNameXi);
5134 hBgrXiBbBCorr->Multiply(histXi,fh1BbBBgrXi[i],1,1,"");
5136 fCorrBgr[fNCorrectionLevelsBgr-1]->AddCorrHistos(i,hBgrTrackPtBbBCorr,hBgrZBbBCorr,hBgrXiBbBCorr);
5140 //_______________________________________________________________________________________________________
5141 void AliFragmentationFunctionCorrections::ReadFoldingCorr(TString strfile, TString strdir, TString strlist)
5144 // temporary histos to hold histos from file
5145 TH1F* hCorrPt[fNJetPtSlices];
5146 TH1F* hCorrZ[fNJetPtSlices];
5147 TH1F* hCorrXi[fNJetPtSlices];
5149 for(Int_t i=0; i<fNJetPtSlices; i++) hCorrPt[i] = 0;
5150 for(Int_t i=0; i<fNJetPtSlices; i++) hCorrZ[i] = 0;
5151 for(Int_t i=0; i<fNJetPtSlices; i++) hCorrXi[i] = 0;
5153 TFile f(strfile,"READ");
5156 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
5160 if(fDebug>0) Printf("%s:%d -- read bin shift correction from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
5162 if(strdir && strdir.Length()) gDirectory->cd(strdir);
5166 if(strlist && strlist.Length()){
5168 if(!(list = (TList*) gDirectory->Get(strlist))){
5169 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
5174 for(Int_t i=0; i<fNJetPtSlices; i++){
5176 Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
5177 Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
5179 TString strNameCorrPt(Form("hCorrFoldingPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
5180 TString strNameCorrZ(Form("hCorrFoldingZ_%02d_%02d",jetPtLoLim,jetPtUpLim));
5181 TString strNameCorrXi(Form("hCorrFoldingXi_%02d_%02d",jetPtLoLim,jetPtUpLim));
5185 hCorrPt[i] = (TH1F*) list->FindObject(strNameCorrPt);
5186 hCorrZ[i] = (TH1F*) list->FindObject(strNameCorrZ);
5187 hCorrXi[i] = (TH1F*) list->FindObject(strNameCorrXi);
5190 hCorrPt[i] = (TH1F*) gDirectory->Get(strNameCorrPt);
5191 hCorrZ[i] = (TH1F*) gDirectory->Get(strNameCorrZ);
5192 hCorrXi[i] = (TH1F*) gDirectory->Get(strNameCorrXi);
5196 //Printf("%s:%d -- error retrieving folding correction %s", (char*)__FILE__,__LINE__,strNameCorrPt.Data());
5200 //Printf("%s:%d -- error retrieving folding correction %s", (char*)__FILE__,__LINE__,strNameCorrZ.Data());
5204 Printf("%s:%d -- error retrieving folding correction %s", (char*)__FILE__,__LINE__,strNameCorrXi.Data());
5207 if(hCorrPt[i]) hCorrPt[i]->SetDirectory(0);
5208 if(hCorrZ[i]) hCorrZ[i]->SetDirectory(0);
5209 if(hCorrXi[i]) hCorrXi[i]->SetDirectory(0);
5211 } // jet slices loop
5215 for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
5216 if(hCorrPt[i]) new(fh1FoldingCorrPt[i]) TH1F(*hCorrPt[i]);
5217 if(hCorrZ[i]) new(fh1FoldingCorrZ[i]) TH1F(*hCorrZ[i]);
5218 if(hCorrXi[i]) new(fh1FoldingCorrXi[i]) TH1F(*hCorrXi[i]);
5222 // ___________________________________________________
5223 void AliFragmentationFunctionCorrections::FoldingCorr()
5225 // apply bin-by-bin correction
5227 AddCorrectionLevel("FoldingCorr");
5229 for(Int_t i=0; i<fNJetPtSlices; i++){
5231 TH1F* histPt = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i);
5232 TH1F* histZ = fCorrFF[fNCorrectionLevels-2]->GetZ(i);
5233 TH1F* histXi = fCorrFF[fNCorrectionLevels-2]->GetXi(i);
5235 TString histNamePt = histPt->GetName();
5236 TString histNameZ = histZ->GetName();
5237 TString histNameXi = histXi->GetName();
5239 std::cout<<" foldingCorr: i "<<i<<" corr pt "<<fh1FoldingCorrPt[i]<<" z "<<fh1FoldingCorrZ[i]<<" xi "
5240 <<fh1FoldingCorrXi[i]<<std::endl;
5242 std::cout<<" foldingCorr: i "<<i<<" mean corr pt "<<fh1FoldingCorrPt[i]->GetMean()<<" z "<<fh1FoldingCorrZ[i]->GetMean()<<" xi "
5243 <<fh1FoldingCorrXi[i]->GetMean()<<std::endl;
5247 TH1F* hFFTrackPtFoldingCorr = (TH1F*) histPt->Clone(histNamePt);
5248 if(fh1FoldingCorrPt[i] && fh1FoldingCorrPt[i]->GetMean()>0) hFFTrackPtFoldingCorr->Multiply(histPt,fh1FoldingCorrPt[i],1,1,"");
5249 else hFFTrackPtFoldingCorr->Reset();
5251 TH1F* hFFZFoldingCorr = (TH1F*) histZ->Clone(histNameZ);
5252 if(fh1FoldingCorrZ[i] && fh1FoldingCorrZ[i]->GetMean()>0) hFFZFoldingCorr->Multiply(histZ,fh1FoldingCorrZ[i],1,1,"");
5253 else hFFZFoldingCorr->Reset();
5255 TH1F* hFFXiFoldingCorr = (TH1F*) histXi->Clone(histNameXi);
5256 if(fh1FoldingCorrXi[i]&& fh1FoldingCorrXi[i]->GetMean()>0) hFFXiFoldingCorr->Multiply(histXi,fh1FoldingCorrXi[i],1,1,"");
5257 else hFFXiFoldingCorr->Reset();
5259 fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hFFTrackPtFoldingCorr,hFFZFoldingCorr,hFFXiFoldingCorr);
5263 //________________________________________________________________________________________________________________________________
5264 void AliFragmentationFunctionCorrections::ReadBgrJetSecCorr(TString strfile, TString strBgrID, TString strdir, TString strlist,
5265 Bool_t useScaledStrangeness)
5268 ReadJetSecCorr(strfile, strdir, strlist, useScaledStrangeness, kTRUE, strBgrID);
5271 //_______________________________________________________________________________________________________
5272 void AliFragmentationFunctionCorrections::ReadJetSecCorr(TString strfile, TString strdir, TString strlist, Bool_t useScaledStrangeness,
5273 Bool_t readBgr, TString strBgrID){
5275 // read reconstruction efficiency from file
5276 // argument strlist optional - read from directory strdir if not specified
5278 // temporary histos to hold histos from file
5279 TH1F* hCorrPt[fNJetPtSlices];
5280 TH1F* hCorrZ[fNJetPtSlices];
5281 TH1F* hCorrXi[fNJetPtSlices];
5283 for(Int_t i=0; i<fNJetPtSlices; i++) hCorrPt[i] = 0;
5284 for(Int_t i=0; i<fNJetPtSlices; i++) hCorrZ[i] = 0;
5285 for(Int_t i=0; i<fNJetPtSlices; i++) hCorrXi[i] = 0;
5287 TFile f(strfile,"READ");
5290 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
5294 if(fDebug>0) Printf("%s:%d -- read secondary correction from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
5296 if(strdir && strdir.Length()) gDirectory->cd(strdir);
5300 if(strlist && strlist.Length()){
5302 if(!(list = (TList*) gDirectory->Get(strlist))){
5303 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
5308 for(Int_t i=0; i<fNJetPtSlices; i++){
5310 Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
5311 Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
5313 TString strNameCorrPt("");
5314 TString strNameCorrZ("");
5315 TString strNameCorrXi("");
5318 strNameCorrPt.Form("hSecCorrBgrPt_%02d_%02d_%s",jetPtLoLim,jetPtUpLim,strBgrID.Data());
5319 strNameCorrZ.Form("hSecCorrBgrZ_%02d_%02d_%s",jetPtLoLim,jetPtUpLim,strBgrID.Data());
5320 strNameCorrXi.Form("hSecCorrBgrXi_%02d_%02d_%s",jetPtLoLim,jetPtUpLim,strBgrID.Data());
5322 if(!useScaledStrangeness){
5323 Printf("%s:%d -- readJetSecCorr bgr: using naive Pythia strangenss ", (char*)__FILE__,__LINE__);
5324 strNameCorrPt.Form("hSecCorrBgrPt_nonSc_%02d_%02d_%s",jetPtLoLim,jetPtUpLim,strBgrID.Data());
5325 strNameCorrZ.Form("hSecCorrBgrZ_nonSc_%02d_%02d_%s",jetPtLoLim,jetPtUpLim,strBgrID.Data());
5326 strNameCorrXi.Form("hSecCorrBgrXi_nonSc_%02d_%02d_%s",jetPtLoLim,jetPtUpLim,strBgrID.Data());
5330 strNameCorrPt.Form("hSecCorrPt_%02d_%02d",jetPtLoLim,jetPtUpLim);
5331 strNameCorrZ.Form("hSecCorrZ_%02d_%02d",jetPtLoLim,jetPtUpLim);
5332 strNameCorrXi.Form("hSecCorrXi_%02d_%02d",jetPtLoLim,jetPtUpLim);
5334 if(!useScaledStrangeness){
5335 Printf("%s:%d -- readJetSecCorr: using naive Pythia strangenss ", (char*)__FILE__,__LINE__);
5336 strNameCorrPt.Form("hSecCorrPt_nonSc_%02d_%02d",jetPtLoLim,jetPtUpLim);
5337 strNameCorrZ.Form("hSecCorrZ_nonSc_%02d_%02d",jetPtLoLim,jetPtUpLim);
5338 strNameCorrXi.Form("hSecCorrXi_nonSc_%02d_%02d",jetPtLoLim,jetPtUpLim);
5343 hCorrPt[i] = (TH1F*) list->FindObject(strNameCorrPt);
5344 hCorrZ[i] = (TH1F*) list->FindObject(strNameCorrZ);
5345 hCorrXi[i] = (TH1F*) list->FindObject(strNameCorrXi);
5348 hCorrPt[i] = (TH1F*) gDirectory->Get(strNameCorrPt);
5349 hCorrZ[i] = (TH1F*) gDirectory->Get(strNameCorrZ);
5350 hCorrXi[i] = (TH1F*) gDirectory->Get(strNameCorrXi);
5354 Printf("%s:%d -- error retrieving secondaries correction %s", (char*)__FILE__,__LINE__,strNameCorrPt.Data());
5358 Printf("%s:%d -- error retrieving secondaries correction %s", (char*)__FILE__,__LINE__,strNameCorrZ.Data());
5362 Printf("%s:%d -- error retrieving secondaries correction %s", (char*)__FILE__,__LINE__,strNameCorrXi.Data());
5365 if(fNHistoBinsPt[i]) hCorrPt[i] = (TH1F*) hCorrPt[i]->Rebin(fNHistoBinsPt[i],strNameCorrPt+"_rebin",fHistoBinsPt[i]->GetArray());
5366 if(fNHistoBinsZ[i]) hCorrZ[i] = (TH1F*) hCorrZ[i]->Rebin(fNHistoBinsZ[i],strNameCorrZ+"_rebin",fHistoBinsZ[i]->GetArray());
5367 if(fNHistoBinsXi[i]) hCorrXi[i] = (TH1F*) hCorrXi[i]->Rebin(fNHistoBinsXi[i],strNameCorrXi+"_rebin",fHistoBinsXi[i]->GetArray());
5369 if(hCorrPt[i]) hCorrPt[i]->SetDirectory(0);
5370 if(hCorrZ[i]) hCorrZ[i]->SetDirectory(0);
5371 if(hCorrXi[i]) hCorrXi[i]->SetDirectory(0);
5373 } // jet slices loop
5377 for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
5380 if(hCorrPt[i]) new(fh1SecCorrBgrPt[i]) TH1F(*hCorrPt[i]);
5381 if(hCorrZ[i]) new(fh1SecCorrBgrZ[i]) TH1F(*hCorrZ[i]);
5382 if(hCorrXi[i]) new(fh1SecCorrBgrXi[i]) TH1F(*hCorrXi[i]);
5385 if(hCorrPt[i]) new(fh1SecCorrPt[i]) TH1F(*hCorrPt[i]);
5386 if(hCorrZ[i]) new(fh1SecCorrZ[i]) TH1F(*hCorrZ[i]);
5387 if(hCorrXi[i]) new(fh1SecCorrXi[i]) TH1F(*hCorrXi[i]);
5392 // ___________________________________________________
5393 void AliFragmentationFunctionCorrections::JetSecCorr()
5395 // apply secondaries correction
5397 AddCorrectionLevel("SecCorr");
5399 Printf("%s:%d -- apply jet secondaries correction", (char*)__FILE__,__LINE__);
5401 for(Int_t i=0; i<fNJetPtSlices; i++){
5403 TH1F* histPt = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i);
5404 TH1F* histZ = fCorrFF[fNCorrectionLevels-2]->GetZ(i);
5405 TH1F* histXi = fCorrFF[fNCorrectionLevels-2]->GetXi(i);
5407 TString histNamePt = histPt->GetName();
5408 TString histNameZ = histZ->GetName();
5409 TString histNameXi = histXi->GetName();
5411 TH1F* hFFTrackPtSecCorr = (TH1F*) histPt->Clone(histNamePt);
5412 hFFTrackPtSecCorr->Multiply(histPt,fh1SecCorrPt[i],1,1,"");
5414 TH1F* hFFZSecCorr = (TH1F*) histZ->Clone(histNameZ);
5415 hFFZSecCorr->Multiply(histZ,fh1SecCorrZ[i],1,1,"");
5417 TH1F* hFFXiSecCorr = (TH1F*) histXi->Clone(histNameXi);
5418 hFFXiSecCorr->Multiply(histXi,fh1SecCorrXi[i],1,1,"");
5420 fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hFFTrackPtSecCorr,hFFZSecCorr,hFFXiSecCorr);
5426 // ___________________________________________________
5427 void AliFragmentationFunctionCorrections::JetSecCorrBgr()
5429 // apply secondaries correction to UE
5431 AddCorrectionLevelBgr("SecCorr");
5433 Printf("%s:%d -- apply jet secondaries correction", (char*)__FILE__,__LINE__);
5435 for(Int_t i=0; i<fNJetPtSlices; i++){
5437 TH1F* histPt = fCorrBgr[fNCorrectionLevelsBgr-2]->GetTrackPt(i);
5438 TH1F* histZ = fCorrBgr[fNCorrectionLevelsBgr-2]->GetZ(i);
5439 TH1F* histXi = fCorrBgr[fNCorrectionLevelsBgr-2]->GetXi(i);
5441 TString histNamePt = histPt->GetName();
5442 TString histNameZ = histZ->GetName();
5443 TString histNameXi = histXi->GetName();
5445 TH1F* hFFTrackPtSecCorr = (TH1F*) histPt->Clone(histNamePt);
5446 hFFTrackPtSecCorr->Multiply(histPt,fh1SecCorrBgrPt[i],1,1,"");
5448 TH1F* hFFZSecCorr = (TH1F*) histZ->Clone(histNameZ);
5449 hFFZSecCorr->Multiply(histZ,fh1SecCorrBgrZ[i],1,1,"");
5451 TH1F* hFFXiSecCorr = (TH1F*) histXi->Clone(histNameXi);
5452 hFFXiSecCorr->Multiply(histXi,fh1SecCorrBgrXi[i],1,1,"");
5454 fCorrBgr[fNCorrectionLevelsBgr-1]->AddCorrHistos(i,hFFTrackPtSecCorr,hFFZSecCorr,hFFXiSecCorr);
5459 //________________________________________________________________________________________________________________
5460 void AliFragmentationFunctionCorrections::WriteBinShiftCorrSinglePt(TString strInfile, TString strIDGen, TString strIDRec,
5461 TString strOutfile, Bool_t updateOutfile, Bool_t useRecPrim,
5464 TString strdirGen = "PWGJE_FragmentationFunction_" + strIDGen;
5465 TString strlistGen = "fracfunc_" + strIDGen;
5467 TString strdirRec = "PWGJE_FragmentationFunction_" + strIDRec;
5468 TString strlistRec = "fracfunc_" + strIDRec;
5470 WriteBinShiftCorrSinglePt(strInfile,strdirGen,strlistGen,strdirRec,strlistRec,strOutfile,updateOutfile,useRecPrim,strOutDir);
5473 //___________________________________________________________________________________________________________________________________
5474 void AliFragmentationFunctionCorrections::WriteBinShiftCorrSinglePt(TString strInfile, TString strdirGen, TString strlistGen,
5475 TString strdirRec, TString strlistRec,
5476 TString strOutfile, Bool_t updateOutfile, Bool_t useRecPrim,
5481 TH1F* hdNdptTracksMCGen = 0;
5482 TH1F* hdNdptTracksMCRec = 0;
5486 TFile f(strInfile,"READ");
5489 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
5493 if(fDebug>0) Printf("%s:%d -- writeBinShiftCorrSinglePt: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
5495 if(strdirGen && strdirGen.Length()) gDirectory->cd(strdirGen);
5499 if(strlistGen && strlistGen.Length()){
5501 if(!(listGen = (TList*) gDirectory->Get(strlistGen))){
5502 Printf("%s:%d -- error retrieving listGen %s from directory %s", (char*)__FILE__,__LINE__,strlistGen.Data(),strdirGen.Data());
5508 hdNdptTracksMCGen = (TH1F*) listGen->FindObject("fh1TrackQAPtGen");
5511 hdNdptTracksMCGen = (TH1F*) gDirectory->Get("fh1TrackQAPtGen");
5514 hdNdptTracksMCGen->SetDirectory(0);
5520 TFile g(strInfile,"READ");
5523 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
5527 if(fDebug>0) Printf("%s:%d -- writeBinShiftCorrSinglePt: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
5529 if(strdirRec && strdirRec.Length()) gDirectory->cd(strdirRec);
5533 if(strlistRec && strlistRec.Length()){
5535 if(!(listRec = (TList*) gDirectory->Get(strlistRec))){
5536 Printf("%s:%d -- error retrieving listRec %s from directory %s", (char*)__FILE__,__LINE__,strlistRec.Data(),strdirRec.Data());
5544 hdNdptTracksMCRec = (TH1F*) listRec->FindObject("fh1TrackQAPtRecEffRec");
5547 hdNdptTracksMCRec = (TH1F*) gDirectory->Get("fh1TrackQAPtRecEffRec");
5552 hdNdptTracksMCRec = (TH1F*) listRec->FindObject("fh1TrackQAPtRecCuts");
5555 hdNdptTracksMCRec = (TH1F*) gDirectory->Get("fh1TrackQAPtRecCuts");
5559 hdNdptTracksMCRec->SetDirectory(0);
5563 TString strNamePtGen = "fh1SinglePtGenBbB";
5564 TString strNamePtRec = "fh1SinglePtRecBbB";
5567 if(fNHistoBinsSinglePt) hdNdptTracksMCGen = (TH1F*) hdNdptTracksMCGen->Rebin(fNHistoBinsSinglePt,strNamePtGen+"_rebin",fHistoBinsSinglePt->GetArray());
5568 if(fNHistoBinsSinglePt) hdNdptTracksMCRec = (TH1F*) hdNdptTracksMCRec->Rebin(fNHistoBinsSinglePt,strNamePtRec+"_rebin",fHistoBinsSinglePt->GetArray());
5570 hdNdptTracksMCGen->SetNameTitle(strNamePtGen,"");
5571 hdNdptTracksMCRec->SetNameTitle(strNamePtRec,"");
5574 TString strTitCorr = "hBbBCorrSinglePt";
5575 if(useRecPrim) strTitCorr = "hBbBCorrRecPrimSinglePt";
5577 hCorrPt = (TH1F*) hdNdptTracksMCGen->Clone(strTitCorr);
5578 hCorrPt->Divide(hdNdptTracksMCGen,hdNdptTracksMCRec,1,1,"B"); // binominal errors
5582 TString outfileOption = "RECREATE";
5583 if(updateOutfile) outfileOption = "UPDATE";
5585 TFile out(strOutfile,outfileOption);
5588 Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
5592 if(fDebug>0) Printf("%s:%d -- write bin shift correction to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
5594 if(strOutDir && strOutDir.Length()){
5597 if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd();
5599 dir = out.mkdir(strOutDir);
5605 hdNdptTracksMCGen->Write();
5606 hdNdptTracksMCRec->Write();
5611 delete hdNdptTracksMCGen;
5612 delete hdNdptTracksMCRec;
5616 //___________________________________________________________________________________________________________________________________
5617 void AliFragmentationFunctionCorrections::ReadBinShiftCorrSinglePt(TString strfile, TString strdir, TString strlist, Bool_t useRecPrim)
5619 // read reconstruction efficiency from file
5620 // argument strlist optional - read from directory strdir if not specified
5622 // temporary histos to hold histos from file
5623 TH1F* hBbBCorrPt = 0;
5625 TFile f(strfile,"READ");
5628 Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
5632 if(fDebug>0) Printf("%s:%d -- read BbB corr from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
5634 if(strdir && strdir.Length()) gDirectory->cd(strdir);
5638 if(strlist && strlist.Length()){
5640 if(!(list = (TList*) gDirectory->Get(strlist))){
5641 Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
5647 TString strNameBbBCorrPt = "hBbBCorrSinglePt";
5648 if(useRecPrim) strNameBbBCorrPt = "hBbBCorrRecPrimSinglePt";
5651 hBbBCorrPt = (TH1F*) list->FindObject(strNameBbBCorrPt);
5654 hBbBCorrPt = (TH1F*) gDirectory->Get(strNameBbBCorrPt);
5658 Printf("%s:%d -- error retrieving BbB corr single pt %s", (char*)__FILE__,__LINE__,strNameBbBCorrPt.Data());
5662 if(fNHistoBinsPt && hBbBCorrPt)
5663 hBbBCorrPt = (TH1F*) hBbBCorrPt->Rebin(fNHistoBinsSinglePt,strNameBbBCorrPt+"_rebin",fHistoBinsSinglePt->GetArray());
5665 if(hBbBCorrPt) hBbBCorrPt->SetDirectory(0);
5669 fh1BbBCorrSinglePt = hBbBCorrPt;
5672 // ------------------------------------------------------------------
5674 void AliFragmentationFunctionCorrections::BbBCorrSinglePt()
5676 // apply efficiency correction to inclusive track pt spec
5678 AddCorrectionLevelSinglePt("BbB");
5680 TH1F* histPt = fCorrSinglePt[fNCorrectionLevelsSinglePt-2]->GetTrackPt(0);
5682 if(histPt->GetNbinsX() != fh1BbBCorrSinglePt->GetNbinsX()){
5683 Printf("%s:%d: inconsistency pt spec and BbB corr binning ", (char*)__FILE__,__LINE__);
5687 TString histNamePt = histPt->GetName();
5688 TH1F* hTrackPtBbBCorr = (TH1F*) histPt->Clone(histNamePt);
5690 hTrackPtBbBCorr->Multiply(histPt,fh1BbBCorrSinglePt,1,1,"");
5692 fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->AddCorrHistos(0,hTrackPtBbBCorr);