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