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