]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/SPECTRA/PiKaPr/TOF/pPb502/macros/DCA.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TOF / pPb502 / macros / DCA.C
CommitLineData
59e49925 1#include "CommonDefs.C"
2
3enum EMCHisto_t {
4 kPrimary,
5 kWeakDecay,
6 kMaterial,
7 kNMCHistos
8};
9const Char_t *mchistoName[kNMCHistos] = {
10 "primary",
11 "weakdecay",
12 "material"
13};
14
15DCAdata(const Char_t *filename, Int_t evMax = kMaxInt, Int_t startEv = 0)
16{
17
18 /* include path for ACLic */
19 gSystem->AddIncludePath("-I$ALICE_ROOT/include");
20 gSystem->AddIncludePath("-I$ALICE_ROOT/TOF");
21 /* load libraries */
22 gSystem->Load("libANALYSIS");
23 gSystem->Load("libANALYSISalice");
24 /* build analysis task class */
25 gROOT->LoadMacro("AliAnalysisParticle.cxx+g");
26 gROOT->LoadMacro("AliAnalysisEvent.cxx+g");
27 gROOT->LoadMacro("AliAnalysisTrack.cxx+g");
28
29 /* open file, get tree and connect */
30 TFile *filein = TFile::Open(filename);
31 TTree *treein = (TTree *)filein->Get("aodTree");
32 printf("got \"aodTree\": %d entries\n", treein->GetEntries());
33 AliAnalysisEvent *analysisEvent = new AliAnalysisEvent();
34 TClonesArray *analysisTrackArray = new TClonesArray("AliAnalysisTrack");
35 AliAnalysisTrack *analysisTrack = NULL;
36 treein->SetBranchAddress("AnalysisEvent", &analysisEvent);
37 treein->SetBranchAddress("AnalysisTrack", &analysisTrackArray);
38
39 /* open enabled flag map */
40 TH1F *hEnabledFlag = NULL;
41 if (enabledChannelsFileName) {
42 TFile *enabledfile = TFile::Open(enabledChannelsFileName);
43 hEnabledFlag = (TH1F *)enabledfile->Get("hEnabledFlag");
44 }
45
46 /**************************************************************/
47 /*** HISTOS ***************************************************/
48 /**************************************************************/
49
50 /* run-time binning */
51 for (Int_t ibin = 0; ibin < NdcaBins + 1; ibin++)
52 dcaBin[ibin] = dcaMin + ibin * dcaStep;
53
54 /* histos */
55 TH3I *hDCAopen[AliPID::kSPECIES][kNCharges];
56 TH3I *hDCAcut[AliPID::kSPECIES][kNCharges];
57 for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) {
58 for (Int_t icharge = 0; icharge < kNCharges; icharge++) {
59 hDCAopen[ipart][icharge] = new TH3I(Form("hDCAopen_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", NcentralityBins, centralityBin, NptBins, ptBin, NdcaBins, dcaBin);
60 hDCAcut[ipart][icharge] = new TH3I(Form("hDCAcut_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", NcentralityBins, centralityBin, NptBins, ptBin, NdcaBins, dcaBin);
61 }
62 }
63
64 /**************************************************************/
65 /**************************************************************/
66 /**************************************************************/
67
68 /* TOF PID response */
69 AliTOFPIDResponse tofResponse;
70 tofResponse.SetTimeResolution(tofReso);
71 /* TPC PID response */
72 AliTPCPIDResponse *tpcResponse = AliAnalysisTrack::GetTPCResponse();
73
74 /* start stopwatch */
75 TStopwatch timer;
76 timer.Start();
77
78 /* loop over events */
79 Int_t charge, index;
80 UShort_t dedxN;
81 Double_t cent, p, pt, mt, dca, tofsignal, tpcsignal, tpctofsignal;
82 Double_t dedx, bethe, deltadedx, dedx_sigma, ptpc;
83 Double_t time, time_sigma, timezero, timezero_sigma, tof, tof_sigma, texp, texp_sigma, deltat, deltat_sigma;
84
85 /* open TZERO calibfile */
86 TH1 *hCentrality_TZEROA_mean = NULL;
87 TH1 *hCentrality_TZEROA_sigma = NULL;
88 TH1 *hCentrality_TZEROC_mean = NULL;
89 TH1 *hCentrality_TZEROC_sigma = NULL;
90 TH1 *hCentrality_TOF_mean = NULL;
91 TH1 *hCentrality_TOF_TZEROA_mean = NULL;
92 TH1 *hCentrality_TOF_TZEROC_mean = NULL;
93 TH1 *hCentrality_TOF_TZEROTOF_mean = NULL;
94 TH1 *hCentrality_TZEROA_reso = NULL;
95 TH1 *hCentrality_TZEROC_reso = NULL;
96 TH1 *hCentrality_TZEROAND_reso = NULL;
97 if (1) {
98 TFile *calibfile = TFile::Open("TZEROcalibration.root");
99 hCentrality_TZEROA_mean = (TH1 *)calibfile->Get("hCentrality_TZEROA_mean");
100 hCentrality_TZEROA_sigma = (TH1 *)calibfile->Get("hCentrality_TZEROA_sigma");
101 hCentrality_TZEROC_mean = (TH1 *)calibfile->Get("hCentrality_TZEROC_calib");
102 hCentrality_TZEROC_sigma = (TH1 *)calibfile->Get("hCentrality_TZEROC_sigma");
103 hCentrality_TOF_mean = (TH1 *)calibfile->Get("hCentrality_TOF_mean");
104 hCentrality_TOF_TZEROA_mean = (TH1 *)calibfile->Get("hCentrality_TOF_TZEROA_mean");
105 hCentrality_TOF_TZEROC_mean = (TH1 *)calibfile->Get("hCentrality_TOF_TZEROC_mean");
106 hCentrality_TOF_TZEROTOF_mean = (TH1 *)calibfile->Get("hCentrality_TOF_TZEROTOF_mean");
107
108 TFile *resofile = TFile::Open("TZEROresolution.root");
109 hCentrality_TZEROA_reso = (TH1 *)resofile->Get("hTZEROA_reso");
110 hCentrality_TZEROC_reso = (TH1 *)resofile->Get("hTZEROC_reso");
111 hCentrality_TZEROAND_reso = (TH1 *)resofile->Get("hTZEROAND_reso");
112 }
113 Double_t TZEROA_mean;
114 Double_t TZEROA_sigma;
115 Double_t TZEROC_mean;
116 Double_t TZEROC_sigma;
117 Double_t TOF_mean;
118 Double_t TOF_TZEROA_mean;
119 Double_t TOF_TZEROC_mean;
120 Double_t TOF_TZEROTOF_mean;
121 Double_t TZEROA;
122 Double_t TZEROA_reso;
123 Bool_t hasTZEROA;
124 Double_t TZEROC;
125 Double_t TZEROC_reso;
126 Bool_t hasTZEROC;
127 Double_t TZEROAND;
128 Double_t TZEROAND_reso;
129 Bool_t hasTZEROAND;
130 Double_t TZEROTOF;
131 Double_t TZEROTOF_reso;
132 Bool_t hasTZEROTOF;
133 Double_t TZEROMEAN;
134 Double_t TZEROMEAN_weight;
135 Double_t TZEROBEST;
136 Double_t TZEROBEST_reso;
137
138
139 for (Int_t iev = startEv; iev < treein->GetEntries() && iev < evMax; iev++) {
140 /* get event */
141 treein->GetEvent(iev);
142 if (iev % 100000 == 0) printf("iev = %d\n", iev);
143 /* check event */
144 if (!analysisEvent->AcceptEvent(acceptEventType)) continue;
145
146 /*** ACCEPTED EVENT ***/
147
148 /* apply time-zero TOF correction */
149 analysisEvent->ApplyTimeZeroTOFCorrection();
150
151 /* get centrality */
152 cent = analysisEvent->GetCentralityPercentile(centralityEstimator);
153
154 /* TZERO corrections */
155 Int_t icent;
156 for (icent = 0; icent < NcentralityBins; icent++)
157 if (cent < centralityBin[icent + 1])
158 break;
159 TZEROA_mean = hCentrality_TZEROA_mean ? hCentrality_TZEROA_mean->GetBinContent(icent + 1) : 0.;
160 TZEROA_sigma = hCentrality_TZEROA_sigma ? hCentrality_TZEROA_sigma->GetBinContent(icent + 1) : 1000.;
161 TZEROC_mean = hCentrality_TZEROC_mean ? hCentrality_TZEROC_mean->GetBinContent(icent + 1) : 0.;
162 TZEROC_sigma = hCentrality_TZEROC_sigma ? hCentrality_TZEROC_sigma->GetBinContent(icent + 1) : 1000.;
163
164 TOF_mean = hCentrality_TOF_mean ? hCentrality_TOF_mean->GetBinContent(icent + 1) : 0.;
165 TOF_TZEROA_mean = hCentrality_TOF_TZEROA_mean ? hCentrality_TOF_TZEROA_mean->GetBinContent(icent + 1) : 0.;
166 TOF_TZEROC_mean = hCentrality_TOF_TZEROC_mean ? hCentrality_TOF_TZEROC_mean->GetBinContent(icent + 1) : 0.;
167 TOF_TZEROTOF_mean = hCentrality_TOF_TZEROTOF_mean ? hCentrality_TOF_TZEROTOF_mean->GetBinContent(icent + 1) : 0.;
168
169 TZEROA_reso = hCentrality_TZEROA_reso ? hCentrality_TZEROA_reso->GetBinContent(icent + 1) : 70.;
170 TZEROC_reso = hCentrality_TZEROC_reso ? hCentrality_TZEROC_reso->GetBinContent(icent + 1) : 70.;
171 TZEROAND_reso = hCentrality_TZEROAND_reso ? hCentrality_TZEROAND_reso->GetBinContent(icent + 1) : 50.;
172
173 /* loop over tracks */
174 for (Int_t itrk = 0; itrk < analysisTrackArray->GetEntries(); itrk++) {
175 /* get track */
176 analysisTrack = (AliAnalysisTrack *)analysisTrackArray->At(itrk);
177 if (!analysisTrack) continue;
178 /* check accepted track (no DCA cut) */
179 if (!analysisTrack->AcceptTrack(kFALSE)) continue;
180 /* get charge */
181 charge = analysisTrack->GetSign() > 0. ? kPositive : kNegative;
182 /* check TPC pid */
183 if (!analysisTrack->HasTPCPID()) continue;
184 /* check TOF pid */
185 if (!analysisTrack->HasTOFPID(hEnabledFlag)) continue;
186
187 /*** ACCEPTED TRACK WITH TPC+TOF PID ***/
188
189 /* get track info */
190 p = analysisTrack->GetP();
191 pt = analysisTrack->GetPt();
192 dca = analysisTrack->GetImpactParameter(0);
193
194 /* get TPC info */
195 dedx = analysisTrack->GetTPCdEdx();
196 dedxN = analysisTrack->GetTPCdEdxN();
197 ptpc = analysisTrack->GetTPCmomentum();
198
199 /* apply expected time correction */
200 // analysisTrack->ApplyTOFExpectedTimeCorrection();
201
202 /* get TOF info */
203 time = analysisTrack->GetTOFTime() - TOF_mean;
204 time_sigma = tofReso;
205 /* TZEROTOF */
206 TZEROTOF = analysisEvent->GetTimeZeroTOF(analysisTrack->GetP());
207 TZEROTOF_reso = analysisEvent->GetTimeZeroTOFSigma(analysisTrack->GetP());
208 hasTZEROTOF = TZEROTOF_reso < 150.;
209 if (hasTZEROTOF) {
210 TZEROTOF += TOF_TZEROTOF_mean - TOF_mean;
211 TZEROTOF_reso *= TZEROTOF_resoScaleFactor;
212 }
213 /* TZEROA */
214 TZEROA = analysisEvent->GetTimeZeroT0(1) - TZEROA_shift;
215 // TZEROA_reso = TZEROA_resolution[icent];
216 hasTZEROA = TMath::Abs(TZEROA - TZEROA_mean) < 3. * TZEROA_sigma;
217 TZEROA += TOF_TZEROA_mean - TOF_mean;
218 /* TZEROC */
219 TZEROC = analysisEvent->GetTimeZeroT0(2) - TZEROC_shift;
220 // TZEROC_reso = TZEROC_resolution[icent];
221 hasTZEROC = TMath::Abs(TZEROC - TZEROC_mean) < 3. * TZEROC_sigma;
222 TZEROC += TOF_TZEROC_mean - TOF_mean;
223 /* TZEROAND */
224 TZEROAND = (TZEROA + TZEROC) * 0.5;
225 // TZEROAND_reso = TZEROAND_resolution[icent];
226 hasTZEROAND = hasTZEROA && hasTZEROC;
227 /* TZEROMEAN */
228 TZEROMEAN = TZEROTOF / TZEROTOF_reso / TZEROTOF_reso;
229 TZEROMEAN_weight = 1. / TZEROTOF_reso / TZEROTOF_reso;
230 if (hasTZEROAND) {
231 // printf("TZEROAND\n");
232 TZEROMEAN += TZEROAND / TZEROAND_reso / TZEROAND_reso;
233 TZEROMEAN_weight = 1. / TZEROAND_reso / TZEROAND_reso;
234 }
235 else if (hasTZEROA) {
236 // printf("TZEROA\n");
237 TZEROMEAN += TZEROA / TZEROA_reso / TZEROA_reso;
238 TZEROMEAN_weight = 1. / TZEROA_reso / TZEROA_reso;
239 }
240 else if (hasTZEROC) {
241 // printf("TZEROC\n");
242 TZEROMEAN += TZEROC / TZEROC_reso / TZEROC_reso;
243 TZEROMEAN_weight = 1. / TZEROC_reso / TZEROC_reso;
244 }
245 timezero = TZEROMEAN / TZEROMEAN_weight;
246 timezero_sigma = TMath::Sqrt(1. / TZEROMEAN_weight);
247 /* TZEROBEST */
248 TZEROBEST = TZEROTOF;
249 TZEROBEST_reso = TZEROTOF_reso;
250 if (hasTZEROAND && TZEROAND_reso < TZEROBEST_reso) {
251 TZEROBEST = TZEROAND;
252 TZEROBEST_reso = TZEROAND_reso;
253 }
254 else if (hasTZEROA && TZEROA_reso < TZEROBEST_reso) {
255 TZEROBEST = TZEROA;
256 TZEROBEST_reso = TZEROA_reso;
257 }
258 if (hasTZEROC && TZEROC_reso < TZEROBEST_reso) {
259 TZEROBEST = TZEROC;
260 TZEROBEST_reso = TZEROC_reso;
261 }
262 timezero = TZEROBEST;
263 timezero_sigma = TZEROBEST_reso;
264
265 tof = time - timezero;
266 tof_sigma = TMath::Sqrt(time_sigma * time_sigma + timezero_sigma * timezero_sigma);
267
268 /* loop over particle IDs */
269 for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) {
270
271 /* check rapidity */
272 if (analysisTrack->GetY(AliPID::ParticleMass(ipart)) - rapidityShift > rapidityMaxCut ||
273 analysisTrack->GetY(AliPID::ParticleMass(ipart)) - rapidityShift < rapidityMinCut) continue;
274
275 /*** ACCEPTED TRACK WITHIN CORRECT RAPIDTY ***/
276
277 /* TPC signal */
278 bethe = tpcResponse->GetExpectedSignal(ptpc, ipart);
279 deltadedx = dedx - bethe;
280 dedx_sigma = tpcResponse->GetExpectedSigma(ptpc, dedxN, ipart);
281 tpcsignal = deltadedx / dedx_sigma;
282
283 /* TOF expected time */
284 texp = analysisTrack->GetTOFExpTime(ipart);
285 texp_sigma = analysisTrack->GetTOFExpTimeSigma(ipart);
286
287 /* TOF signal */
288 deltat = tof - texp;
289 deltat_sigma = TMath::Sqrt(tof_sigma * tof_sigma + texp_sigma * texp_sigma);
290 tofsignal = deltat / deltat_sigma;
291
292 /* TPC+TOF signal */
293 tpctofsignal = TMath::Sqrt(tpcsignal * tpcsignal + tofsignal * tofsignal);
294
295 /* check PID cuts */
296 if (tpctofsignal > 2.) continue;
297
298 /* fill histo */
299 hDCAopen[ipart][charge]->Fill(cent, pt, dca);
300
301 /* check accepted track (with DCA cut) */
302 if (!analysisTrack->AcceptTrack(kTRUE)) continue;
303
304 /* fill histo */
305 hDCAcut[ipart][charge]->Fill(cent, pt, dca);
306
307 } /* end of loop over particle IDs */
308 } /* end of loop over tracks */
309 } /* end of loop over events */
310
311 /* stop stopwatch */
312 timer.Stop();
313 timer.Print();
314
315 /* output */
316 TFile *fileout = TFile::Open(Form("DCAdata.%s", filename), "RECREATE");
317 for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) {
318 for (Int_t icharge = 0; icharge < kNCharges; icharge++) {
319 hDCAopen[ipart][icharge]->Write();
320 hDCAcut[ipart][icharge]->Write();
321 }
322 }
323
324 fileout->Close();
325
326}
327
328/**************************************************************/
329
330DCAmc(const Char_t *filename, Int_t evMax = kMaxInt, Int_t startEv = 0)
331{
332
333 /* include path for ACLic */
334 gSystem->AddIncludePath("-I$ALICE_ROOT/include");
335 gSystem->AddIncludePath("-I$ALICE_ROOT/TOF");
336 /* load libraries */
337 gSystem->Load("libANALYSIS");
338 gSystem->Load("libANALYSISalice");
339 /* build analysis task class */
340 gROOT->LoadMacro("AliAnalysisParticle.cxx+g");
341 gROOT->LoadMacro("AliAnalysisEvent.cxx+g");
342 gROOT->LoadMacro("AliAnalysisTrack.cxx+g");
343
344 /* open file, get tree and connect */
345 TFile *filein = TFile::Open(filename);
346 TTree *treein = (TTree *)filein->Get("aodTree");
347 printf("got \"aodTree\": %d entries\n", treein->GetEntries());
348 AliAnalysisEvent *analysisEvent = new AliAnalysisEvent();
349 TClonesArray *analysisTrackArray = new TClonesArray("AliAnalysisTrack");
350 AliAnalysisTrack *analysisTrack = NULL;
351 treein->SetBranchAddress("AnalysisEvent", &analysisEvent);
352 treein->SetBranchAddress("AnalysisTrack", &analysisTrackArray);
353
354 /**************************************************************/
355 /*** HISTOS ***************************************************/
356 /**************************************************************/
357
358 /* run-time binning */
359 for (Int_t ibin = 0; ibin < NdcaBins + 1; ibin++)
360 dcaBin[ibin] = dcaMin + ibin * dcaStep;
361
362 /* histos */
363 TH3I *hDCAopen[kNMCHistos][AliPID::kSPECIES][kNCharges];
364 TH3I *hDCAcut[kNMCHistos][AliPID::kSPECIES][kNCharges];
365 for (Int_t ihisto = 0; ihisto < kNMCHistos; ihisto++ ) {
366 for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) {
367 for (Int_t icharge = 0; icharge < kNCharges; icharge++) {
368 hDCAopen[ihisto][ipart][icharge] = new TH3I(Form("hDCAopen_%s_%s_%s", mchistoName[ihisto], AliPID::ParticleName(ipart), chargeName[icharge]), "", NcentralityBins, centralityBin, NptBins, ptBin, NdcaBins, dcaBin);
369 hDCAcut[ihisto][ipart][icharge] = new TH3I(Form("hDCAcut_%s_%s_%s", mchistoName[ihisto], AliPID::ParticleName(ipart), chargeName[icharge]), "", NcentralityBins, centralityBin, NptBins, ptBin, NdcaBins, dcaBin);
370 }
371 }
372 }
373
374 /**************************************************************/
375 /**************************************************************/
376 /**************************************************************/
377
378 /* start stopwatch */
379 TStopwatch timer;
380 timer.Start();
381
382 /* loop over events */
383 Int_t part, charge, type;
384 Double_t cent, pt, mt, dca;
385
386 for (Int_t iev = startEv; iev < treein->GetEntries() && iev < evMax; iev++) {
387 /* get event */
388 treein->GetEvent(iev);
389 if (iev % 100000 == 0) printf("iev = %d\n", iev);
390 /* check event */
391 if (!analysisEvent->AcceptEvent(acceptEventType)) continue;
392
393 /*** ACCEPTED EVENT ***/
394
395 /* get centrality */
396 cent = analysisEvent->GetCentralityPercentile(centralityEstimator);
397
398 /* loop over tracks */
399 for (Int_t itrk = 0; itrk < analysisTrackArray->GetEntries(); itrk++) {
400 /* get track */
401 analysisTrack = (AliAnalysisTrack *)analysisTrackArray->At(itrk);
402 if (!analysisTrack) continue;
403 /* check defined PID */
404 part = analysisTrack->GetMCPID();
405 if (part < 0 || analysisTrack->GetSign() == 0.) continue;
406 /* check accepted track (no DCA cut) */
407 if (!analysisTrack->AcceptTrack(kFALSE)) continue;
408 /* check rapidity */
409 if (analysisTrack->GetY(AliPID::ParticleMass(part)) - rapidityShift > rapidityMaxCut ||
410 analysisTrack->GetY(AliPID::ParticleMass(part)) - rapidityShift < rapidityMinCut) continue;
411
412 /*** ACCEPTED TRACK WITHIN CORRECT RAPIDTY ***/
413
414 /* get track info */
415 pt = analysisTrack->GetPt();
416 dca = analysisTrack->GetImpactParameter(0);
417 charge = analysisTrack->GetMCCharge() > 0. ? kPositive : kNegative;
418 if (analysisTrack->IsMCPrimary())
419 type = kPrimary;
420 else if (analysisTrack->IsMCSecondaryWeakDecay())
421 type = kWeakDecay;
422 else if (analysisTrack->IsMCSecondaryMaterial())
423 type = kMaterial;
424 else
425 continue;
426
427 /* fill histo */
428 hDCAopen[type][part][charge]->Fill(cent, pt, dca);
429
430 /* check accepted track (with DCA cut) */
431 if (!analysisTrack->AcceptTrack(kTRUE)) continue;
432
433 /* fill histo */
434 hDCAcut[type][part][charge]->Fill(cent, pt, dca);
435
436 } /* end of loop over tracks */
437 } /* end of loop over events */
438
439 /* stop stopwatch */
440 timer.Stop();
441 timer.Print();
442
443 /* output */
444 TFile *fileout = TFile::Open(Form("DCAmc.%s", filename), "RECREATE");
445 for (Int_t ihisto = 0; ihisto < kNMCHistos; ihisto++) {
446 for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) {
447 for (Int_t icharge = 0; icharge < kNCharges; icharge++) {
448 hDCAopen[ihisto][ipart][icharge]->Write();
449 hDCAcut[ihisto][ipart][icharge]->Write();
450 }
451 }
452 }
453
454 fileout->Close();
455
456}
457
458/**************************************************************/
459
460enum EFitParams_t {
461 kIntegralCounts,
462 kPrimaryCounts,
463 kWeakDecayCounts,
464 kMaterialCounts,
465 kPrimaryIntegral,
466 kWeakDecayIntegral,
467 kMaterialIntegral,
468 kDataCutIntegral,
469 kPrimaryCutIntegral,
470 kWeakDecayCutIntegral,
471 kMaterialCutIntegral,
472 kNFitParams
473};
474
475/* fit output params name */
476const Char_t *fitParamName[kNFitParams] = {
477 "IntegralCounts",
478 "PrimaryCounts",
479 "WeakDecayCounts",
480 "MaterialCounts",
481 "PrimaryIntegral",
482 "WeakDecayIntegral",
483 "MaterialIntegral",
484 "DataCutIntegral",
485 "PrimaryCutIntegral",
486 "WeakDecayCutIntegral",
487 "MaterialCutIntegral"
488};
489
490/* fit output params title */
491const Char_t *fitParamTitle[kNFitParams] = {
492 "Integral counts;p_{T} (GeV/c);",
493 "Primary counts;p_{T} (GeV/c);",
494 "Weak-decay counts;p_{T} (GeV/c);",
495 "Material counts;p_{T} (GeV/c);",
496 "Primary integral;p_{T} (GeV/c);",
497 "Weak-decay integral;p_{T} (GeV/c);",
498 "Material integral;p_{T} (GeV/c);",
499 "Data cut integral;p_{T} (GeV/c);",
500 "Primary cut integral;p_{T} (GeV/c);",
501 "Weak-decay cut integral;p_{T} (GeV/c);",
502 "Material cut integral;p_{T} (GeV/c);"
503};
504
505/* fit ranges */
506Double_t fitPtMin[AliPID::kSPECIES] = {0.5, 0.5, 0.3, 0.4, 0.5};
507Double_t fitPtMax[AliPID::kSPECIES] = {2.0, 2.0, 2.0, 2.0, 3.0};
508
509/* rebin DCA */
510Int_t rebindca = 10;
511
512DCAdeltafeed(const Char_t *datafilename, const Char_t *mcfilename)
513{
514 for (Int_t icharge = 0; icharge < kNCharges; icharge++) {
515 DCAdeltafeed(datafilename, mcfilename, 2, icharge, -1);
516 DCAdeltafeed(datafilename, mcfilename, 4, icharge, -1);
517 for (Int_t icent = 0; icent < NcentralityBins; icent++) {
518 DCAdeltafeed(datafilename, mcfilename, 2, icharge, icent);
519 DCAdeltafeed(datafilename, mcfilename, 4, icharge, icent);
520 }
521 }
522}
523
524DCAdeltafeed(const Char_t *datafilename, const Char_t *mcfilename, Int_t ipart, Int_t icharge, Int_t icent, Float_t ptMin = -1., Float_t ptMax = -1., Bool_t checkHistoFlag = kFALSE)
525{
526
527 printf("****************************************\n");
528 printf("RUNNING DCA FIT:\n");
529 printf("RAPIDITY-CUT: %s\n", AliPID::ParticleName(ipart));
530 printf("CHARGE: %s\n", chargeName[icharge]);
531 printf("PARTICLE: %s\n", AliPID::ParticleName(ipart));
532 printf("CENTRALITY BIN: %d\n", icent);
533 printf("****************************************\n");
534
535 /* open data */
536 TFile *datafilein = TFile::Open(datafilename);
537 TFile *mcfilein = TFile::Open(mcfilename);
538
539 /* get histos */
540 TH3I *hDCAopen = (TH3I *)datafilein->Get(Form("hDCAopen_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]));
541 TH3I *hDCAopen_primary = (TH3I *)mcfilein->Get(Form("hDCAopen_primary_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]));
542 TH3I *hDCAopen_weakdecay = (TH3I *)mcfilein->Get(Form("hDCAopen_weakdecay_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]));
543 TH3I *hDCAopen_material = (TH3I *)mcfilein->Get(Form("hDCAopen_material_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]));
544 TH3I *hDCAcut = (TH3I *)datafilein->Get(Form("hDCAcut_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]));
545 TH3I *hDCAcut_primary = (TH3I *)mcfilein->Get(Form("hDCAcut_primary_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]));
546 TH3I *hDCAcut_weakdecay = (TH3I *)mcfilein->Get(Form("hDCAcut_weakdecay_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]));
547 TH3I *hDCAcut_material = (TH3I *)mcfilein->Get(Form("hDCAcut_material_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]));
548
549 /* setup centrality range */
550 if (icent < 0 || icent >= NcentralityBins) {
551 printf("WARNING: undefined centrality -> using 00-90\% range\n");
552 hDCAopen->GetXaxis()->SetRange(1, NcentralityBins);
553 hDCAopen_primary->GetXaxis()->SetRange(1, NcentralityBins);
554 hDCAopen_weakdecay->GetXaxis()->SetRange(1, NcentralityBins);
555 hDCAopen_material->GetXaxis()->SetRange(1, NcentralityBins);
556 hDCAcut->GetXaxis()->SetRange(1, NcentralityBins);
557 hDCAcut_primary->GetXaxis()->SetRange(1, NcentralityBins);
558 hDCAcut_weakdecay->GetXaxis()->SetRange(1, NcentralityBins);
559 hDCAcut_material->GetXaxis()->SetRange(1, NcentralityBins);
560 }
561 else {
562 printf("***** FITTING CENTRALITY-BIN [%02d, %02d] %% *****\n", centralityBin[icent], centralityBin[icent + 1]);
563 hDCAopen->GetXaxis()->SetRange(icent + 1, icent + 1);
564#if 0
565 hDCAopen_primary->GetXaxis()->SetRange(icent + 1, icent + 1);
566 hDCAopen_weakdecay->GetXaxis()->SetRange(icent + 1, icent + 1);
567 hDCAopen_material->GetXaxis()->SetRange(icent + 1, icent + 1);
568#else
569 hDCAopen_primary->GetXaxis()->SetRange(1, NcentralityBins);
570 hDCAopen_weakdecay->GetXaxis()->SetRange(1, NcentralityBins);
571 hDCAopen_material->GetXaxis()->SetRange(1, NcentralityBins);
572#endif
573 hDCAcut->GetXaxis()->SetRange(icent + 1, icent + 1);
574#if 0
575 hDCAcut_primary->GetXaxis()->SetRange(icent + 1, icent + 1);
576 hDCAcut_weakdecay->GetXaxis()->SetRange(icent + 1, icent + 1);
577 hDCAcut_material->GetXaxis()->SetRange(icent + 1, icent + 1);
578#else
579 hDCAcut_primary->GetXaxis()->SetRange(1, NcentralityBins);
580 hDCAcut_weakdecay->GetXaxis()->SetRange(1, NcentralityBins);
581 hDCAcut_material->GetXaxis()->SetRange(1, NcentralityBins);
582#endif
583 }
584
585 /* setup pt range */
586 Bool_t requestedRange = kFALSE;
587 if (ptMin > -0.001 && ptMax > -0.001 && ptMax > ptMin) {
588 printf("***** FITTING PT-BIN [%f, %f] GeV/c *****\n", ptMin, ptMax);
589 requestedRange = kTRUE;
590 hDCAopen->GetYaxis()->SetRangeUser(ptMin + 0.001, ptMax - 0.001);
591 hDCAopen_primary->GetYaxis()->SetRangeUser(ptMin + 0.001, ptMax - 0.001);
592 hDCAopen_weakdecay->GetYaxis()->SetRangeUser(ptMin + 0.001, ptMax - 0.001);
593 hDCAopen_material->GetYaxis()->SetRangeUser(ptMin + 0.001, ptMax - 0.001);
594 hDCAcut->GetYaxis()->SetRangeUser(ptMin + 0.001, ptMax - 0.001);
595 hDCAcut_primary->GetYaxis()->SetRangeUser(ptMin + 0.001, ptMax - 0.001);
596 hDCAcut_weakdecay->GetYaxis()->SetRangeUser(ptMin + 0.001, ptMax - 0.001);
597 hDCAcut_material->GetYaxis()->SetRangeUser(ptMin + 0.001, ptMax - 0.001);
598 }
599
600 /* output */
601 Char_t outfilename[1024];
602 if (icent < 0 || icent >= NcentralityBins)
603 sprintf(outfilename, "DCAdeltafeed_cent0090_%s_%s.root", AliPID::ParticleName(ipart), chargeName[icharge]);
604 else {
605 sprintf(outfilename, "DCAdeltafeed_cent%02d%02d_%s_%s.root", centralityBin[icent], centralityBin[icent + 1], AliPID::ParticleName(ipart), chargeName[icharge]);
606 }
607 TFile *fileout = TFile::Open(outfilename, "RECREATE");
608 TDirectory *fitDir = fileout->mkdir("FitParams");
609 /* canvas */
610 TCanvas *canvas = new TCanvas("canvas");
611 canvas->SetLogy();
612 /* histo */
613 TH1D *hFitParamHisto[kNFitParams];
614 for (Int_t iparam = 0; iparam < kNFitParams; iparam++)
615 hFitParamHisto[iparam] = new TH1D(Form("h%s", fitParamName[iparam]), fitParamTitle[iparam], NptBins, ptBin);
616
617
618 /* loop over ptBins */
619 for (Int_t ipt = 0; ipt < NptBins; ipt++) {
620
621 if (!requestedRange) {
622 if ((ptBin[ipt] + 0.001) < fitPtMin[ipart] || (ptBin[ipt + 1] - 0.001) > fitPtMax[ipart]) continue;
623 printf("***** FITTING PT-BIN [%f, %f] GeV/c *****\n", ptBin[ipt], ptBin[ipt + 1]);
624 hDCAopen->GetYaxis()->SetRange(ipt + 1, ipt + 1);
625 hDCAopen_primary->GetYaxis()->SetRange(ipt + 1, ipt + 1);
626 hDCAopen_weakdecay->GetYaxis()->SetRange(ipt + 1, ipt + 1);
627 hDCAopen_material->GetYaxis()->SetRange(ipt + 1, ipt + 1);
628 hDCAcut->GetYaxis()->SetRange(ipt + 1, ipt + 1);
629 hDCAcut_primary->GetYaxis()->SetRange(ipt + 1, ipt + 1);
630 hDCAcut_weakdecay->GetYaxis()->SetRange(ipt + 1, ipt + 1);
631 hDCAcut_material->GetYaxis()->SetRange(ipt + 1, ipt + 1);
632 }
633
634 /* dca projections */
635 TH1 *hDCAopen_py = hDCAopen->Project3D("z");
636 TH1 *hDCAopen_primary_py = hDCAopen_primary->Project3D("z");
637 TH1 *hDCAopen_weakdecay_py = hDCAopen_weakdecay->Project3D("z");
638 TH1 *hDCAopen_material_py = hDCAopen_material->Project3D("z");
639 TH1 *hDCAcut_py = hDCAcut->Project3D("z");
640 TH1 *hDCAcut_primary_py = hDCAcut_primary->Project3D("z");
641 TH1 *hDCAcut_weakdecay_py = hDCAcut_weakdecay->Project3D("z");
642 TH1 *hDCAcut_material_py = hDCAcut_material->Project3D("z");
643
644 /* rebin */
645 hDCAopen_py->Rebin(rebindca);
646 hDCAopen_primary_py->Rebin(rebindca);
647 hDCAopen_weakdecay_py->Rebin(rebindca);
648 hDCAopen_material_py->Rebin(rebindca);
649 hDCAcut_py->Rebin(rebindca);
650 hDCAcut_primary_py->Rebin(rebindca);
651 hDCAcut_weakdecay_py->Rebin(rebindca);
652 hDCAcut_material_py->Rebin(rebindca);
653
654 /* check histos if requested */
655 if (checkHistoFlag) {
656 TCanvas *cCheckHisto = new TCanvas("cCheckHisto");
657 cCheckHisto->Divide(2, 4);
658 cCheckHisto->cd(1);
659 hDCAopen_py->Draw();
660 cCheckHisto->cd(2);
661 hDCAcut_py->Draw();
662 cCheckHisto->cd(3);
663 hDCAopen_primary_py->Draw();
664 cCheckHisto->cd(4);
665 hDCAcut_primary_py->Draw();
666 cCheckHisto->cd(5);
667 hDCAopen_weakdecay_py->Draw();
668 cCheckHisto->cd(6);
669 hDCAcut_weakdecay_py->Draw();
670 cCheckHisto->cd(7);
671 hDCAopen_material_py->Draw();
672 cCheckHisto->cd(8);
673 hDCAcut_material_py->Draw();
674 return;
675 }
676
677 Double_t param[kNFitParams];
678 Double_t param_err[kNFitParams];
679
680 param[kPrimaryIntegral] = hDCAopen_primary_py->Integral();
681 param_err[kPrimaryIntegral] = TMath::Sqrt(hDCAopen_primary_py->Integral());
682 param[kWeakDecayIntegral] = hDCAopen_weakdecay_py->Integral();
683 param_err[kWeakDecayIntegral] = TMath::Sqrt(hDCAopen_weakdecay_py->Integral());
684 param[kMaterialIntegral] = hDCAopen_material_py->Integral();
685 param_err[kMaterialIntegral] = TMath::Sqrt(hDCAopen_material_py->Integral());
686
687 param[kDataCutIntegral] = hDCAcut_py->Integral();
688 param_err[kDataCutIntegral] = TMath::Sqrt(hDCAcut_py->Integral());
689 param[kPrimaryCutIntegral] = hDCAcut_primary_py->Integral();
690 param_err[kPrimaryCutIntegral] = TMath::Sqrt(hDCAcut_primary_py->Integral());
691 param[kWeakDecayCutIntegral] = hDCAcut_weakdecay_py->Integral();
692 param_err[kWeakDecayCutIntegral] = TMath::Sqrt(hDCAcut_weakdecay_py->Integral());
693 param[kMaterialCutIntegral] = hDCAcut_material_py->Integral();
694 param_err[kMaterialCutIntegral] = TMath::Sqrt(hDCAcut_material_py->Integral());
695
696 /* fit */
697 DCAfit(hDCAopen_py, hDCAopen_primary_py, hDCAopen_weakdecay_py, hDCAopen_material_py, param, param_err, canvas);
698
699
700 /* check requested pt-range */
701 if (requestedRange)
702 return;
703
704 /* write canvas */
705 fitDir->cd();
706 canvas->Write(Form("fitDisplay_ptBin_%3.2f_%3.2f", ptBin[ipt], ptBin[ipt + 1]));
707
708 /* set histo */
709 for (Int_t iparam = 0; iparam < kNFitParams; iparam++) {
710 hFitParamHisto[iparam]->SetBinContent(ipt + 1, param[iparam]);
711 hFitParamHisto[iparam]->SetBinError(ipt + 1, param_err[iparam]);
712 }
713
714 /* delete */
715 delete hDCAopen_py;
716 delete hDCAopen_primary_py;
717 delete hDCAopen_weakdecay_py;
718 delete hDCAopen_material_py;
719 delete hDCAcut_py;
720 delete hDCAcut_primary_py;
721 delete hDCAcut_weakdecay_py;
722 delete hDCAcut_material_py;
723
724 }
725
726 /* check requested pt-range */
727 if (requestedRange)
728 return;
729
730 /*** POST-ANALYSIS ***/
731
732 TDirectory *postDir = fileout->mkdir("PostAnalysis");
733
734 /* compute fractions */
735 TH1D *hPrimaryFraction = new TH1D(*hFitParamHisto[kPrimaryCounts]);
736 hPrimaryFraction->Divide(hFitParamHisto[kPrimaryCounts], hFitParamHisto[kIntegralCounts], 1., 1., "B");
737 TH1D *hWeakDecayFraction = new TH1D(*hFitParamHisto[kWeakDecayCounts]);
738 hWeakDecayFraction->Divide(hFitParamHisto[kWeakDecayCounts], hFitParamHisto[kIntegralCounts], 1., 1., "B");
739 TH1D *hMaterialFraction = new TH1D(*hFitParamHisto[kMaterialCounts]);
740 hMaterialFraction->Divide(hFitParamHisto[kMaterialCounts], hFitParamHisto[kIntegralCounts], 1., 1., "B");
741
742 /* compute scale factors */
743 TH1D *hPrimaryScale = new TH1D(*hFitParamHisto[kPrimaryCounts]);
744 hPrimaryScale->Divide(hFitParamHisto[kPrimaryCounts], hFitParamHisto[kPrimaryIntegral], 1., 1., "B");
745 TH1D *hWeakDecayScale = new TH1D(*hFitParamHisto[kWeakDecayCounts]);
746 hWeakDecayScale->Divide(hFitParamHisto[kWeakDecayCounts], hFitParamHisto[kWeakDecayIntegral], 1., 1., "B");
747 TH1D *hMaterialScale = new TH1D(*hFitParamHisto[kMaterialCounts]);
748 hMaterialScale->Divide(hFitParamHisto[kMaterialCounts], hFitParamHisto[kMaterialIntegral], 1., 1., "B");
749
750 /* compute cut integrals */
751 TH1D *hPrimaryCutIntegral = new TH1D(*hFitParamHisto[kPrimaryCutIntegral]);
752 hPrimaryCutIntegral->Multiply(hPrimaryScale);
753 TH1D *hWeakDecayCutIntegral = new TH1D(*hFitParamHisto[kWeakDecayCutIntegral]);
754 hWeakDecayCutIntegral->Multiply(hWeakDecayScale);
755 TH1D *hMaterialCutIntegral = new TH1D(*hFitParamHisto[kMaterialCutIntegral]);
756 hMaterialCutIntegral->Multiply(hMaterialScale);
757
758 /* compute cut fractions */
759 TH1D *hPrimaryCutFraction = new TH1D(*hPrimaryCutIntegral);
760 hPrimaryCutFraction->Divide(hPrimaryCutIntegral, hFitParamHisto[kDataCutIntegral], 1., 1., "B");
761 TH1D *hWeakDecayCutFraction = new TH1D(*hWeakDecayCutIntegral);
762 hWeakDecayCutFraction->Divide(hWeakDecayCutIntegral, hFitParamHisto[kDataCutIntegral], 1., 1., "B");
763 TH1D *hMaterialCutFraction = new TH1D(*hMaterialCutIntegral);
764 hMaterialCutFraction->Divide(hMaterialCutIntegral, hFitParamHisto[kDataCutIntegral], 1., 1., "B");
765
766
767 /*** OUTPUT ***/
768
769 /* write fir params histos */
770 fitDir->cd();
771 for (Int_t iparam = 0; iparam < kNFitParams; iparam++)
772 hFitParamHisto[iparam]->Write();
773 /* write post-analysis histos */
774 postDir->cd();
775 hPrimaryFraction->Write("hPrimaryFraction");
776 hWeakDecayFraction->Write("hWeakDecayFraction");
777 hMaterialFraction->Write("hMaterialFraction");
778 hPrimaryCutFraction->Write("hPrimaryCutFraction");
779 hWeakDecayCutFraction->Write("hWeakDecayCutFraction");
780 hMaterialCutFraction->Write("hMaterialCutFraction");
781
782 /* cleanup */
783 delete hPrimaryFraction;
784 delete hWeakDecayFraction;
785 delete hMaterialFraction;
786 delete hPrimaryScale;
787 delete hWeakDecayScale;
788 delete hMaterialScale;
789 delete hPrimaryCutIntegral;
790 delete hWeakDecayCutIntegral;
791 delete hMaterialCutIntegral;
792 delete hPrimaryCutFraction;
793 delete hWeakDecayCutFraction;
794 delete hMaterialCutFraction;
795 delete canvas;
796 for (Int_t iparam = 0; iparam < kNFitParams; iparam++)
797 delete hFitParamHisto[iparam];
798
799 /* close file */
800 fileout->Close();
801
802
803}
804
805/**************************************************************/
806
807//___________________________________________________________________________________
808
809Float_t
810TOFpid_histomin(TH1 *h)
811{
812
813 for (Int_t ibin = 0; ibin < h->GetNbinsX(); ibin++)
814 if (h->GetBinContent(ibin + 1) > 0.)
815 return h->GetXaxis()->GetBinCenter(ibin + 1);
816 return kMaxInt;
817}
818
819//___________________________________________________________________________________
820
821Float_t
822TOFpid_histomax(TH1 *h)
823{
824
825 for (Int_t ibin = h->GetNbinsX(); ibin > 0; ibin--)
826 if (h->GetBinContent(ibin) > 0.)
827 return h->GetXaxis()->GetBinCenter(ibin);
828 return -kMaxInt;
829}
830
831//___________________________________________________________________________________
832
833TOFpid_checkneg(TH1 *h)
834{
835
836 for (Int_t ibin = 0; ibin < h->GetNbinsX(); ibin++)
837 if (h->GetBinContent(ibin + 1) <= 0.) {
838 h->SetBinContent(ibin + 1, 1.e-300);
839 // h->SetBinError(ibin + 1, 0.1);
840 }
841}
842
843
844Double_t
845DCAfit(TH1 *hData, TH1 *hPrimary, TH1 *hWeakDecay, TH1 *hMaterial, Double_t *param = NULL, Double_t *param_err = NULL, TCanvas *canvas = NULL)
846{
847
848 /** ROOFIT ***/
849 gSystem->Load("libRooFit");
850 using namespace RooFit;
851
852 /*** DEFINE FIT RANGE ***/
853
854 printf("***** FIT RANGE DEFINITION *****\n");
855
856 /* check material histogram to define min/max fit range */
857 Double_t rangeMin = TMath::Max(dcaMin, TOFpid_histomin(hMaterial));
858 Double_t rangeMax = TMath::Min(dcaMax, TOFpid_histomax(hMaterial));
859 /* fix zeroes */
860 TOFpid_checkneg(hMaterial);
861
862 /* define range */
863 RooRealVar x("x", "DCA_{xy}", 0., -5., 5., "");
864 printf("FIT RANGE DEFINED: %f -> %f\n", dcaMin, dcaMax);
865 printf("********************************\n");
866
867 /*** DEFINE HISTOGRAM DATA ***/
868
869 /* define data to fit and background from input histogram */
870 RooDataHist hdata("hdata", "hdata", x, hData);
871 RooDataHist hprimary("hprimary", "hprimary", x, hPrimary);
872 RooDataHist hweakdecay("hweakdecay", "hweakdecay", x, hWeakDecay);
873 RooDataHist hmaterial("hmaterial", "hmaterial", x, hMaterial);
874
875 /*** DEFINE SHAPES ***/
876
877 RooHistPdf primary("primary", "primary", x, hprimary);
878 RooHistPdf weakdecay("weakdecay", "weakdecay", x, hweakdecay);
879 RooHistPdf material("material", "material", x, hmaterial);
880
881 /*** DEFINE MODEL ***/
882
883 Double_t integral = hdata.sumEntries();
884 RooRealVar nprimary("nprimary", "nprimary", 0.80 * integral, 0., integral);
885 RooRealVar nweakdecay("nweakdecay", "nweakdecay", 0.1 * integral, 0., integral);
886 RooRealVar nmaterial("nmaterial", "nmaterial", 0.1 * integral, 0., integral);
887
888 RooAddPdf model("model", "model p.d.f.", RooArgList(primary, weakdecay, material), RooArgList(nprimary, nweakdecay, nmaterial));
889
890 /*** FIT ***/
891
892 if (canvas) canvas->cd();
893 model.fitTo(hdata, Extended(kTRUE), SumW2Error(kFALSE));
894
895 /*** DRAW ***/
896
897 RooPlot *xframe = x.frame();
898 hdata.plotOn(xframe, XErrorSize(0), DrawOption("PZ"));
899 model.plotOn(xframe, LineWidth(2), Precision(1.e-4));
900 model.plotOn(xframe, Components(primary), LineWidth(2), LineColor(kRed), Precision(1.e-4));
901 model.plotOn(xframe, Components(weakdecay), LineWidth(2), LineStyle(kDashed), LineColor(kCyan+1));
902 model.plotOn(xframe, Components(material), LineWidth(2), LineStyle(kDashed), LineColor(kGreen+1));
903 xframe->SetMinimum(0.1);
904 xframe->Draw();
905 if (canvas) canvas->Update();
906
907 printf("*****************************\n");
908 printf("***** FRACTIONS *****\n");
909 printf("primary = %f +- %f\n", nprimary.getVal(), nprimary.getError());
910 printf("weak-decay = %f +- %f\n", nweakdecay.getVal(), nweakdecay.getError());
911 printf("material = %f +- %f\n", nmaterial.getVal(), nmaterial.getError());
912 printf("******************\n");
913
914 /*** OUTPUT FIT PARAMS ***/
915
916 param[kIntegralCounts] = integral;
917 param_err[kIntegralCounts] = sqrt(integral);
918 param[kPrimaryCounts] = nprimary.getVal();
919 param_err[kPrimaryCounts] = nprimary.getError();
920 param[kWeakDecayCounts] = nweakdecay.getVal();
921 param_err[kWeakDecayCounts] = nweakdecay.getError();
922 param[kMaterialCounts] = nmaterial.getVal();
923 param_err[kMaterialCounts] = nmaterial.getError();
924
925 return;
926
927}
928
929enum EFitFunc_t {
930 kExpo,
931 kInverseExpo,
932 kPowerLaw,
933 kInversePowerLaw,
934 kExpoPowerLaw
935};
936
937DCAdeltafeed_param()
938{
939#if 0
940 DCAdeltafeed_param(2, 0, "WeakDecayCutFraction", kExpo);
941 DCAdeltafeed_param(2, 1, "WeakDecayCutFraction", kExpo);
942 DCAdeltafeed_param(4, 0, "WeakDecayCutFraction", kExpo);
943 DCAdeltafeed_param(4, 1, "WeakDecayCutFraction", kExpo);
944 DCAdeltafeed_param(2, 0, "MaterialCutFraction", kPowerLaw);
945 DCAdeltafeed_param(2, 1, "MaterialCutFraction", kPowerLaw);
946 DCAdeltafeed_param(4, 0, "MaterialCutFraction", kPowerLaw);
947 DCAdeltafeed_param(4, 1, "MaterialCutFraction", kPowerLaw);
948#endif
949 DCAdeltafeed_param_bothcharges(2, 0, "PrimaryCutFraction", kInverseExpo);
950 DCAdeltafeed_param_bothcharges(2, 1, "PrimaryCutFraction", kInverseExpo);
951 DCAdeltafeed_param_bothcharges(4, 0, "PrimaryCutFraction", kInverseExpo);
952 DCAdeltafeed_param_bothcharges(4, 1, "PrimaryCutFraction", kInverseExpo);
953}
954
955DCAdeltafeed_param_bothcharges(Int_t ipart, Int_t icharge, const Char_t *name, Int_t fitFunc)
956{
957
958 TVirtualFitter::SetMaxIterations(1000000);
959
960 /* load HistoUtils */
961 gROOT->LoadMacro("HistoUtils.C");
962
963 Char_t outfilename[1024];
964 TH1D *hDeltaFeed_mb[kNCharges];
965 TH1D *hDeltaFeed_cent[kNCharges][NcentralityBins];
966 for (Int_t iicharge = 0; iicharge < kNCharges; iicharge++) {
967
968 /* get minimum-bias results */
969 sprintf(outfilename, "DCAdeltafeed_cent0090_%s_%s.root", AliPID::ParticleName(ipart), chargeName[iicharge]);
970 TFile *filein = TFile::Open(outfilename);
971 hDeltaFeed_mb[iicharge] = (TH1D *)filein->Get(Form("PostAnalysis/h%s", name));
972 if (!hDeltaFeed_mb[iicharge]) {
973 printf("cannot find PostAnalysis/h%s in %s\n", name, outfilename);
974 return;
975 }
976
977 /* get centrality-binned results */
978 for (Int_t icent = 0; icent < NcentralityBins; icent++) {
979 sprintf(outfilename, "DCAdeltafeed_cent%02d%02d_%s_%s.root", centralityBin[icent], centralityBin[icent + 1], AliPID::ParticleName(ipart), chargeName[iicharge]);
980 filein = TFile::Open(outfilename);
981 if (!filein || !filein->IsOpen()) continue;
982 hDeltaFeed_cent[iicharge][icent] = (TH1D *)filein->Get(Form("PostAnalysis/h%s", name));
983 if (!hDeltaFeed_cent[iicharge][icent]) {
984 printf("cannot find PostAnalysis/h%s in %s\n", name, outfilename);
985 return;
986 }
987 }
988 }
989
990 /* define delta feed-down function */
991 switch (fitFunc) {
992 case kExpo:
993 printf("expo function selected\n");
994 TF1 *fDeltaFeed = new TF1("fDeltaFeed_expo", "[0] * ([1] + [2] * TMath::Exp([3] * x))", 0.5, 5.);
995 fDeltaFeed->SetParameter(1, 0.01);
996 fDeltaFeed->SetParLimits(1, 0., 0.1);
997 fDeltaFeed->SetParameter(2, 0.1);
998 fDeltaFeed->SetParameter(3, -1.);
999 break;
1000 case kInverseExpo:
1001 printf("inverse-expo function selected\n");
1002 TF1 *fDeltaFeed = new TF1("fDeltaFeed_expo", "[0] * ([1] + [2] * TMath::Exp([3] * x))", 0.5, 5.);
1003 fDeltaFeed->SetParameter(1, 1.);
1004 fDeltaFeed->SetParLimits(1, 0.9, 1.0);
1005 fDeltaFeed->SetParameter(2, -0.1);
1006 fDeltaFeed->SetParLimits(2, -1., 0.);
1007 fDeltaFeed->SetParameter(3, -1.);
1008 break;
1009 case kPowerLaw:
1010 printf("powerlaw function selected\n");
1011 TF1 *fDeltaFeed = new TF1("fDeltaFeed_powerlaw", "[0] * ([1] + [2] * TMath::Power(x, [3]))", 0.5, 5.);
1012 fDeltaFeed->SetParameter(1, 0.01);
1013 fDeltaFeed->SetParLimits(1, 0., 0.1);
1014 fDeltaFeed->SetParameter(2, 0.1);
1015 fDeltaFeed->SetParameter(3, -1.);
1016 break;
1017 case kInversePowerLaw:
1018 printf("inverse-powerlaw function selected\n");
1019 TF1 *fDeltaFeed = new TF1("fDeltaFeed_powerlaw", "[0] * ([1] + [2] * TMath::Power(x, [3]) + [4] * x)", 0.5, 5.);
1020 fDeltaFeed->SetParameter(1, 1.);
1021 fDeltaFeed->SetParLimits(1, 0.9, 1.);
1022 fDeltaFeed->SetParameter(2, -0.1);
1023 fDeltaFeed->SetParameter(3, -1.);
1024 fDeltaFeed->SetParameter(4, 0.);
1025 fDeltaFeed->SetParLimits(4, 0., 1.)
1026 break;
1027 case kExpoPowerLaw:
1028 printf("expo+powerlaw function selected\n");
1029 TF1 *fDeltaFeed = new TF1("fDeltaFeed_expopowerlaw", "[0] * ([1] + [2] * TMath::Exp([3] * x) + [4] * TMath::Exp([5] * x))", 0.5, 5.);
1030
1031 fDeltaFeed->SetParameter(1, 0.95);
1032 fDeltaFeed->SetParLimits(1, 0.9, 1.);
1033 fDeltaFeed->SetParameter(2, -0.1);
1034 fDeltaFeed->SetParLimits(2, -1000., 0.);
1035 fDeltaFeed->SetParameter(3, -1.);
1036 fDeltaFeed->SetParameter(4, -100.);
1037 fDeltaFeed->SetParLimits(4, -100., 0.);
1038 fDeltaFeed->SetParameter(5, -20.);
1039 fDeltaFeed->SetParLimits(5, -100., 10.);
1040 break;
1041 default:
1042 printf("fit function not defined\n");
1043 return;
1044 break;
1045 }
1046 fDeltaFeed->FixParameter(0, 1.);
1047
1048 /* output */
1049 TFile *fileout = TFile::Open(Form("DCAdeltafeed_param_%s_%s_%s.root", name, AliPID::ParticleName(ipart), chargeName[icharge]), "RECREATE");
1050
1051 /* extra-material function */
1052 TF1 *fExtraMaterial = new TF1("fExtraMaterial", "1. + [0] * TMath::Exp([1] * x)", 0., 5.);
1053 fExtraMaterial->SetParameter(0, 0.1);
1054 fExtraMaterial->SetParameter(1, -0.5);
1055
1056 /* loop over centrality bins */
1057 for (Int_t icent = -1; icent < NcentralityBins; icent++) {
1058
1059 printf(" icent = %d\n", icent);
1060
1061 /* fit material from proton-antiproton ratio and remove it */
1062 if (ipart == AliPID::kProton) {
1063 printf("fitting extra-material for positive protons\n");
1064 if (icent == -1) {
1065 TH1D *hr = new TH1D(*hDeltaFeed_mb[kPositive]);
1066 hr->Divide(hDeltaFeed_mb[kNegative]);
1067 hr->Fit(fExtraMaterial, "R", "", fitPtMin[ipart], fitPtMax[ipart]);
1068 hDeltaFeed_mb[kPositive]->Divide(fExtraMaterial);
1069 }
1070 else {
1071 TH1D *hr = new TH1D(*hDeltaFeed_cent[kPositive][icent]);
1072 hr->Divide(hDeltaFeed_cent[kNegative][icent]);
1073 hr->Fit(fExtraMaterial, "R", "", fitPtMin[ipart], fitPtMax[ipart]);
1074 hDeltaFeed_cent[kPositive][icent]->Divide(fExtraMaterial);
1075 }
1076 }
1077
1078 /* build graph with both charges */
1079 TGraphErrors *gDeltaFeed = new TGraphErrors();
1080 Int_t npoints = 0;
1081 Double_t pt, pte, val, vale;
1082 for (Int_t iicharge = 0; iicharge < kNCharges; iicharge++) {
1083 for (Int_t ipt = 0; ipt < NptBins; ipt++) {
1084 if (icent == -1) {
1085 if (hDeltaFeed_mb[iicharge]->GetBinError(ipt + 1) == 0.)
1086 continue;
1087 pt = hDeltaFeed_mb[iicharge]->GetBinCenter(ipt + 1);
1088 pte = 0.5 * hDeltaFeed_mb[iicharge]->GetBinWidth(ipt + 1);
1089 val = hDeltaFeed_mb[iicharge]->GetBinContent(ipt + 1);
1090 vale = hDeltaFeed_mb[iicharge]->GetBinError(ipt + 1);
1091 gDeltaFeed->SetPoint(npoints, pt, val);
1092 gDeltaFeed->SetPointError(npoints, pte, vale);
1093 npoints++;
1094 }
1095 else {
1096 if (hDeltaFeed_cent[iicharge][icent]->GetBinError(ipt + 1) == 0.)
1097 continue;
1098 pt = hDeltaFeed_cent[iicharge][icent]->GetBinCenter(ipt + 1);
1099 pte = 0.5 * hDeltaFeed_cent[iicharge][icent]->GetBinWidth(ipt + 1);
1100 val = hDeltaFeed_cent[iicharge][icent]->GetBinContent(ipt + 1);
1101 vale = hDeltaFeed_cent[iicharge][icent]->GetBinError(ipt + 1);
1102 gDeltaFeed->SetPoint(npoints, pt, val);
1103 gDeltaFeed->SetPointError(npoints, pte, vale);
1104 npoints++;
1105 }
1106 }
1107 }
1108
1109 /* fit graph */
1110 Int_t fitres = -1;
1111 fitres = gDeltaFeed->Fit(fDeltaFeed, "R", "", 0.5, 2.0);
1112 if (fitres != 0) {
1113 fDeltaFeed->SetParLimits(1, 0.9, 1.0);
1114 fDeltaFeed->SetParameter(2, -0.1);
1115 fDeltaFeed->SetParLimits(2, -1., 0.);
1116 fDeltaFeed->SetParameter(3, -1.);
1117 }
1118 fitres = -1;
1119 while(fitres != 0)
1120 fitres = gDeltaFeed->Fit(fDeltaFeed, "R", "", fitPtMin[ipart], fitPtMax[ipart]);
1121
1122 gDeltaFeed->Draw("ap*");
1123 fDeltaFeed->Draw("same");
1124 gPad->Update();
1125
1126 printf("done part = %d, charge = %d, cent = %d\n", ipart, icharge, icent);
1127
1128 /* build fit profile */
1129 TProfile *pDeltaFeed = new TProfile("pDeltaFeed", "", NptBins, ptBin);
1130 HistoUtils_Function2Profile(fDeltaFeed, pDeltaFeed);
1131
1132 /* put material back */
1133 if (ipart == AliPID::kProton && icharge == kPositive) {
1134 if (icent == -1) {
1135 hDeltaFeed_mb[kPositive]->Multiply(fExtraMaterial);
1136 pDeltaFeed->Multiply(fExtraMaterial);
1137 }
1138 else {
1139 hDeltaFeed_cent[kPositive][icent]->Multiply(fExtraMaterial);
1140 pDeltaFeed->Multiply(fExtraMaterial);
1141 }
1142 }
1143
1144 /* write */
1145 fileout->cd();
1146 if (icent == -1) {
1147 hDeltaFeed_mb[icharge]->Write(Form("hDeltaFeed_MB"));
1148 fDeltaFeed->Write(Form("fDeltaFeed_MB"));
1149 pDeltaFeed->Write(Form("pDeltaFeed_MB"));
1150 fExtraMaterial->Write(Form("fExtraMaterial_MB"));
1151 }
1152 else {
1153 hDeltaFeed_cent[icharge][icent]->Write(Form("hDeltaFeed_cent%d", icent));
1154 fDeltaFeed->Write(Form("fDeltaFeed_cent%d", icent));
1155 pDeltaFeed->Write(Form("pDeltaFeed_cent%d", icent));
1156 fExtraMaterial->Write(Form("fExtraMaterial_cent%d", icent));
1157 }
1158 }
1159
1160 fileout->Close();
1161}
1162
1163DCAdeltafeed_param(Int_t ipart, Int_t icharge, const Char_t *name, Int_t fitFunc)
1164{
1165
1166 TVirtualFitter::SetMaxIterations(1000000);
1167
1168 /* load HistoUtils */
1169 gROOT->LoadMacro("HistoUtils.C");
1170
1171 /* get minimum-bias results */
1172 Char_t outfilename[1024];
1173 sprintf(outfilename, "DCAdeltafeed_cent0090_%s_%s.root", AliPID::ParticleName(ipart), chargeName[icharge]);
1174 TFile *filein = TFile::Open(outfilename);
1175 TH1D *hDeltaFeed_mb = (TH1D *)filein->Get(Form("PostAnalysis/h%s", name));
1176 if (!hDeltaFeed_mb) {
1177 printf("cannot find PostAnalysis/h%s in %s\n", name, outfilename);
1178 return;
1179 }
1180
1181 /* get centrality-binned results */
1182 TH1D *hDeltaFeed_cent[NcentralityBins];
1183 for (Int_t icent = 0; icent < NcentralityBins; icent++) {
1184 sprintf(outfilename, "DCAdeltafeed_cent%02d%02d_%s_%s.root", centralityBin[icent], centralityBin[icent + 1], AliPID::ParticleName(ipart), chargeName[icharge]);
1185 filein = TFile::Open(outfilename);
1186 if (!filein || !filein->IsOpen()) continue;
1187 hDeltaFeed_cent[icent] = (TH1D *)filein->Get(Form("PostAnalysis/h%s", name));
1188 if (!hDeltaFeed_cent[icent]) {
1189 printf("cannot find PostAnalysis/h%s in %s\n", name, outfilename);
1190 return;
1191 }
1192 }
1193
1194 /* define delta feed-down function */
1195 switch (fitFunc) {
1196 case kExpo:
1197 printf("expo function selected\n");
1198 TF1 *fDeltaFeed = new TF1("fDeltaFeed_expo", "[0] * ([1] + [2] * TMath::Exp([3] * x))", 0.5, 5.);
1199 fDeltaFeed->SetParameter(1, 0.01);
1200 fDeltaFeed->SetParLimits(1, 0., 0.1);
1201 fDeltaFeed->SetParameter(2, 0.1);
1202 fDeltaFeed->SetParameter(3, -1.);
1203 break;
1204 case kInverseExpo:
1205 printf("inverse-expo function selected\n");
1206 TF1 *fDeltaFeed = new TF1("fDeltaFeed_expo", "[0] * ([1] + [2] * TMath::Exp([3] * x))", 0.5, 5.);
1207 fDeltaFeed->SetParameter(1, 1.);
1208 fDeltaFeed->SetParLimits(1, 0.9, 1.0);
1209 fDeltaFeed->SetParameter(2, -0.1);
1210 fDeltaFeed->SetParameter(3, -1.);
1211 break;
1212 case kPowerLaw:
1213 printf("powerlaw function selected\n");
1214 TF1 *fDeltaFeed = new TF1("fDeltaFeed_powerlaw", "[0] * ([1] + [2] * TMath::Power(x, [3]))", 0.5, 5.);
1215 fDeltaFeed->SetParameter(1, 0.01);
1216 fDeltaFeed->SetParLimits(1, 0., 0.1);
1217 fDeltaFeed->SetParameter(2, 0.1);
1218 fDeltaFeed->SetParameter(3, -1.);
1219 break;
1220 case kInversePowerLaw:
1221 printf("inverse-powerlaw function selected\n");
1222 TF1 *fDeltaFeed = new TF1("fDeltaFeed_powerlaw", "[0] * ([1] + [2] * TMath::Power(x, [3]) + [4] * x)", 0.5, 5.);
1223 fDeltaFeed->SetParameter(1, 1.);
1224 fDeltaFeed->SetParLimits(1, 0.9, 1.);
1225 fDeltaFeed->SetParameter(2, -0.1);
1226 fDeltaFeed->SetParameter(3, -1.);
1227 fDeltaFeed->SetParameter(4, 0.);
1228 fDeltaFeed->SetParLimits(4, 0., 1.)
1229 break;
1230 case kExpoPowerLaw:
1231 printf("expo+powerlaw function selected\n");
1232 // TF1 *fDeltaFeed = new TF1("fDeltaFeed_expopowerlaw", "[0] * ([1] + [2] * TMath::Exp([3] * x) + [4] * TMath::Exp([5] * x))", 0.5, 5.);
1233 TF1 *fDeltaFeed = new TF1("fDeltaFeed_expopowerlaw", "[0] * ([1] + [2] * TMath::Exp([3] * x) * TMath::Exp([5] / x))", 0.5, 5.);
1234
1235 fDeltaFeed->SetParameter(1, 0.95);
1236 fDeltaFeed->SetParLimits(1, 0.9, 1.);
1237 fDeltaFeed->SetParameter(2, -0.1);
1238 fDeltaFeed->SetParLimits(2, -1000., 0.);
1239 fDeltaFeed->SetParameter(3, -1.);
1240 fDeltaFeed->SetParameter(4, -100.);
1241 fDeltaFeed->SetParLimits(4, -100., 0.);
1242 fDeltaFeed->SetParameter(5, -20.);
1243 fDeltaFeed->SetParLimits(5, -100., 10.);
1244 break;
1245 default:
1246 printf("fit function not defined\n");
1247 return;
1248 break;
1249 }
1250 fDeltaFeed->FixParameter(0, 1.);
1251
1252 /* output */
1253 TFile *fileout = TFile::Open(Form("DCAdeltafeed_param_%s_%s_%s.root", name, AliPID::ParticleName(ipart), chargeName[icharge]), "RECREATE");
1254
1255 /* fit minimum-bias */
1256 TCanvas *cMinimumBias = new TCanvas("cMinimumBias");
1257 if (fitFunc == kExpoPowerLaw) {
1258 fDeltaFeed->FixParameter(4, 0.);
1259 hDeltaFeed_mb->Fit(fDeltaFeed, "IMRE", "", 1.0, 3.0);
1260 fDeltaFeed->FixParameter(1, fDeltaFeed->GetParameter(1));
1261 fDeltaFeed->FixParameter(2, fDeltaFeed->GetParameter(2));
1262 fDeltaFeed->FixParameter(3, fDeltaFeed->GetParameter(3));
1263 fDeltaFeed->ReleaseParameter(4);
1264 hDeltaFeed_mb->Fit(fDeltaFeed, "IMRE", "", 0.5, 3.0);
1265 }
1266 else {
1267 hDeltaFeed_mb->Fit(fDeltaFeed, "IMRE", "", 0.5, 3.0);
1268 }
1269 cMinimumBias->SetLogy();
1270
1271 /* build fit profile */
1272 TProfile *pDeltaFeed = new TProfile("pDeltaFeed", "", NptBins, ptBin);
1273 HistoUtils_Function2Profile(fDeltaFeed, pDeltaFeed);
1274
1275 /* save MB params */
1276 Double_t parMB[100];
1277 for (Int_t ipar = 0; ipar < fDeltaFeed->GetNpar(); ipar++)
1278 parMB[ipar] = fDeltaFeed->GetParameter(ipar);
1279
1280 /* write */
1281 fileout->cd();
1282 hDeltaFeed_mb->Write("hDeltaFeed_mb");
1283 fDeltaFeed->Write("fDeltaFeed_mb");
1284 pDeltaFeed->Write("pDeltaFeed_mb");
1285
1286 /* fix pt-depencence and release scale factor
1287 to fit centrality bins */
1288 TCanvas *cCentralityDependence = new TCanvas("cCentralityDependence");
1289 TH1D *hCentDep = new TH1D("hCentDep", "", NcentralityBins, centralityBin);
1290 //fDeltaFeed->ReleaseParameter(0);
1291 // for (Int_t ipar = 1; ipar < fDeltaFeed->GetNpar(); ipar++)
1292 // fDeltaFeed->FixParameter(ipar, fDeltaFeed->GetParameter(ipar));
1293 //fDeltaFeed->FixParameter(1, fDeltaFeed->GetParameter(1));
1294 TProfile *pDeltaFeed_cent[NcentralityBins]; = new TProfile("pDeltaFeed", "", NptBins, ptBin);
1295 for (Int_t icent = 0; icent < NcentralityBins; icent++) {
1296 if (!hDeltaFeed_cent[icent]) continue;
1297
1298 /* restore MB parameters */
1299 for (Int_t ipar = 0; ipar < fDeltaFeed->GetNpar(); ipar++)
1300 fDeltaFeed->SetParameter(ipar, parMB[ipar]);
1301
1302 if (fitFunc == kExpoPowerLaw) {
1303 fDeltaFeed->FixParameter(4, 0.);
1304 hDeltaFeed_mb->Fit(fDeltaFeed, "IMRE", "", 1.0, 3.0);
1305 fDeltaFeed->FixParameter(1, fDeltaFeed->GetParameter(1));
1306 fDeltaFeed->FixParameter(2, fDeltaFeed->GetParameter(2));
1307 fDeltaFeed->FixParameter(3, fDeltaFeed->GetParameter(3));
1308 fDeltaFeed->ReleaseParameter(4);
1309 hDeltaFeed_cent[icent]->Fit(fDeltaFeed, "IMRE", "", 0.5, 3.0);
1310 }
1311 else {
1312 hDeltaFeed_cent[icent]->Fit(fDeltaFeed, "IMRE", "", 0.5, 3.0);
1313 }
1314
1315 pDeltaFeed_cent[icent] = new TProfile(Form("pDeltaFeed_cent%d", icent), "", NptBins, ptBin);
1316 HistoUtils_Function2Profile(fDeltaFeed, pDeltaFeed_cent[icent]);
1317 hCentDep->SetBinContent(icent + 1, fDeltaFeed->GetParameter(0.));
1318 hCentDep->SetBinError(icent + 1, fDeltaFeed->GetParError(0.));
1319
1320 /* write */
1321 fileout->cd();
1322 hDeltaFeed_cent[icent]->Write(Form("hDeltaFeed_cent%d", icent));
1323 fDeltaFeed->Write(Form("fDeltaFeed_cent%d", icent));
1324 pDeltaFeed_cent[icent]->Write(Form("pDeltaFeed_cent%d", icent));
1325
1326 }
1327
1328 /* fit centrality dependence */
1329 TF1 *fCentDep = (TF1 *)gROOT->GetFunction("pol3");
1330 // hCentDep->Fit(fCentDep);
1331
1332 /* build fit profile */
1333 TProfile *pCentDep = new TProfile("pCentDep", "", NcentralityBins, centralityBin);
1334 HistoUtils_Function2Profile(fCentDep, pCentDep);
1335
1336 /* write */
1337 fileout->cd();
1338 hCentDep->Write("hCentDep");
1339 fCentDep->Write("fCentDep");
1340 pCentDep->Write("pCentDep");
1341 fileout->Close();
1342
1343}
1344
1345DCA_primaryFraction(const Char_t *destdir)
1346{
1347
1348 Int_t marker[2] = {20, 21};
1349 Int_t color[AliPID::kSPECIES] = {1, 1, 4, 8, 2};
1350 Char_t *partLatex[AliPID::kSPECIES][2] = {
1351 "", "", "", "", "#pi^{+}", "#pi^{-}", "K^{+}", "K^{-}", "p", "#bar{p}"
1352 };
1353
1354 TFile *fileout = TFile::Open("TOF_primaryFraction.root", "RECREATE");
1355
1356 TProfile *pFrac;
1357 TH1D *hFrac = new TH1D("hFrac", "", NptBins, ptBin);
1358 Char_t title[1024];
1359 for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++) {
1360 if (ipart == 3) continue;
1361 for (Int_t icharge = 0; icharge < kNCharges; icharge++) {
1362 TFile *filein = TFile::Open(Form("%s/DCAdeltafeed_param_PrimaryCutFraction_%s_%s.root", destdir, AliPID::ParticleName(ipart), chargeName[icharge]));
1363 for (Int_t icent = -1; icent < NcentralityBins; icent++) {
1364 if (icent == -1)
1365 pFrac = (TProfile *)filein->Get(Form("pDeltaFeed_MB"));
1366 else
1367 pFrac = (TProfile *)filein->Get(Form("pDeltaFeed_cent%d", icent));
1368 /* copy profile in TH1D */
1369 for (Int_t ipt = 0; ipt < NptBins; ipt++) {
1370 hFrac->SetBinContent(ipt + 1, pFrac->GetBinContent(ipt + 1));
1371 hFrac->SetBinError(ipt + 1, pFrac->GetBinError(ipt + 1));
1372 }
1373 if (icent == -1)
1374 sprintf(title, "%s (MB);p_{T} (GeV/c);fraction of primary particle;", partLatex[ipart][icharge]);
1375 else
1376 sprintf(title, "%s (%d-%d%%);p_{T} (GeV/c);fraction of primary particle;", partLatex[ipart][icharge], (Int_t)centralityBin[icent], (Int_t)centralityBin[icent + 1]);
1377 hFrac->SetMarkerStyle(marker[icharge]);
1378 hFrac->SetMarkerColor(color[ipart]);
1379 hFrac->SetTitle(title);
1380 if (icent == -1)
1381 hFrac->SetName(Form("hPrimaryFrac_MB_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]));
1382 else
1383 hFrac->SetName(Form("hPrimaryFrac_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
1384 fileout->cd();
1385 hFrac->Write();
1386 }
1387 }
1388 }
1389
1390 fileout->Close();
1391
1392}