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