]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/SPECTRA/PiKaPr/TOF/pPb502/macros/TOFpid.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 / TOFpid.C
CommitLineData
59e49925 1#include "CommonDefs.C"
2
3#define DISPLAY 0
4
5/* SIGNAL SHAPE */
6Bool_t GAUSSIAN_SIGNAL = kFALSE;
7Bool_t GAUSSIANTAIL_SIGNAL = kFALSE;
8Bool_t GAUSSIANTAIL2_SIGNAL = kFALSE;
9Bool_t GAUSSIANPLUSGAUSSIANTAIL_SIGNAL = kFALSE;
10Bool_t GAUSSIANPLUSEXPONENTIAL_SIGNAL = kFALSE;
11Bool_t EXPECTED_SIGNAL_TAIL = kTRUE;
12Bool_t EXPECTED_SIGNAL_TEMPLATE = kTRUE;
13Bool_t EXPECTED_SIGNAL_FIT = kFALSE;
14/* SIGNAL PARAMETERS */
15Bool_t FIX_SIGNAL_MEAN = kFALSE;
16Bool_t FIX_SIGNAL_SIGMA = kFALSE;
17Bool_t FIX_SIGNAL_TAIL = kFALSE;
18Float_t SCALE_SIGNAL_SIGMA = 1.;
19Float_t SCALE_SIGNAL_TAIL = 1.;
20/* OTHER STUFF */
21Char_t *SIGNAL_PARAM_FILE = NULL;//"signalParamFile.root";
22Float_t DEFAULT_SIGNAL_MEAN = 0.;
23Float_t MIN_SIGNAL_MEAN = -0.2;
24Float_t MAX_SIGNAL_MEAN = 0.2.;
25Float_t DEFAULT_SIGNAL_SIGMA = 1.;
26Float_t MIN_SIGNAL_SIGMA = 0.8;
27Float_t MAX_SIGNAL_SIGMA = 1.2;
28Float_t DEFAULT_SIGNAL_TAIL = 1.;
29Float_t MIN_SIGNAL_TAIL = 0.5;
30Float_t MAX_SIGNAL_TAIL = 1.5;
31/* BACKGROUND */
32Bool_t EXPECTED_BACKGROUND_TAIL = kTRUE;
33Bool_t EXPECTED_BACKGROUND_TEMPLATE = kTRUE;
34Bool_t EXPECTED_BACKGROUND_FIT = kFALSE;
35Bool_t GAUSSIAN_BACKGROUND = kFALSE;
36Bool_t USE_ELECTRON_BACKGROUND = kTRUE;
37/* BACKGROUND PARAMETERS */
38Bool_t FIX_BACKGROUND_MEAN = kTRUE;
39Bool_t FIX_BACKGROUND_SIGMA = kTRUE;
40Bool_t FIX_BACKGROUND_TAIL = kTRUE;
41Float_t SCALE_BACKGROUND_SIGMA = 1.;
42Float_t SCALE_BACKGROUND_TAIL = 1.;
43/* MISMATCH */
44Bool_t NO_MISMATCH = kFALSE;
45Bool_t EXPECTED_MISMATCH = kTRUE;
46Bool_t EXPONENTIAL_MISMATCH = kFALSE;
47Bool_t UNIFORM_MISMATCH = kFALSE;
48Bool_t DOUBLEEXPONENTIAL_MISMATCH = kFALSE;
49Bool_t UNIFORMPLUSEXPONENTIAL_MISMATCH = kFALSE;
50/* AUTOMATIC CONSTRAINS */
51Float_t FIT_ELECTRONS_PT_MIN = 0.;
52Float_t FIT_ELECTRONS_PT_MAX = 1.0;
53Float_t FIT_MUONS_PT_MIN = 0.7;
54Float_t FIT_MUONS_PT_MAX = 0.;
55Float_t FIT_PIONS_PT_MIN = 0.;
56Float_t FIT_PIONS_PT_MAX = 5.;
57
58Float_t CONSTRAINSIGNAL_LIMIT = 0.;
59Float_t CONSTRAINTAIL_LIMIT = 0.;
60
61enum EFitParams_t {
62 kMean,
63 kSigma,
64 kTail,
65 kTotalCounts,
66 kIntegralCounts,
67 kSignalCounts,
68 kBkg1Counts,
69 kBkg2Counts,
70 kBkg3Counts,
71 kBkg4Counts,
72 kMismatchCounts,
73 kNFitParams
74};
75
76/* fit output params name */
77const Char_t *fitParamName[kNFitParams] = {
78 "Mean",
79 "Sigma",
80 "Tail",
81 "TotalCounts",
82 "IntegralCounts",
83 "SignalCounts",
84 "Bkg1Counts",
85 "Bkg2Counts",
86 "Bkg3Counts",
87 "Bkg4Counts",
88 "MismatchCounts"
89};
90
91/* fit output params title */
92const Char_t *fitParamTitle[kNFitParams] = {
93 "Signal mean;p_{T} (GeV/c);#mu (ps)",
94 "Signal sigma;p_{T} (GeV/c);#sigma (ps)",
95 "Signal tail;p_{T} (GeV/c);#sigma_{tail} (ps)",
96 "Total counts;p_{T} (GeV/c);counts",
97 "Total counts within fit range;p_{T} (GeV/c);counts",
98 "Signal counts;p_{T} (GeV/c);counts",
99 "Background-1 counts;p_{T} (GeV/c);counts",
100 "Background-2 counts;p_{T} (GeV/c);counts",
101 "Background-3 counts;p_{T} (GeV/c);counts",
102 "Background-4 counts;p_{T} (GeV/c);counts",
103 "Mismatch counts within fit range;p_{T} (GeV/c);counts"
104};
105
106/**************************************************************/
107/*** GENERATION OF TEMPLATE HISTOS ****************************/
108/**************************************************************/
109
110const Int_t NmismatchTrials = 1;
111const Int_t NexpectedTrials = 1;
112
113/**************************************************************/
114/*** HISTOS AND BINNING ***************************************/
115/**************************************************************/
116
117/**************************************************************/
118/**************************************************************/
119enum EHistoParam_t {
120 kCentrality,
121 kTOFsigma,
122 kPt,
123 kTPCsigma,
124 kNHistoParams
125};
126/**************************************************************/
127/**************************************************************/
128const Int_t NtofsigmaBins = 1750;
129Double_t tofsigmaBin[NtofsigmaBins + 1];
130Double_t tofsigmaMin = -100., tofsigmaMax = 250., tofsigmaStep = (tofsigmaMax - tofsigmaMin) / NtofsigmaBins;
131/**************************************************************/
132const Int_t NtofsignalBins = 2000;
133Double_t tofsignalBin[NtofsignalBins + 1];
134Double_t tofsignalMin = -24400., tofsignalMax = 24400., tofsignalStep = (tofsignalMax - tofsignalMin) / NtofsignalBins;
135/**************************************************************/
136const Int_t NtofresoBins = 500;
137Double_t tofresoBin[NtofsignalBins + 1];
138Double_t tofresoMin = 0., tofresoMax = 500., tofresoStep = (tofresoMax - tofresoMin) / NtofresoBins;
139/**************************************************************/
140/**************************************************************/
141const Int_t NtpcsigmaBins = 10;
142Double_t tpcsigmaBin[NtpcsigmaBins + 1];
143Double_t tpcsigmaMin = -5., tpcsigmaMax = 5., tpcsigmaStep = (tpcsigmaMax - tpcsigmaMin) / NtpcsigmaBins;
144/**************************************************************/
145Int_t NparamsBins[kNHistoParams] = {NcentralityBins, NtofsigmaBins, NptBins, NtpcsigmaBins};
146Double_t *paramBin[kNHistoParams] = {centralityBin, tofsigmaBin, ptBin, tpcsigmaBin};
147/**************************************************************/
148
149/**************************************************************/
150/**************************************************************/
151/**************************************************************/
152
153/**************************************************************/
154
155AliTOFGeometry tofGeo;
156Float_t c = TMath::C() * 1.e2 / 1.e12; /* cm/ps */
157Float_t c_1 = 1. / c;
158
159Double_t
160GenerateRandomHit(TH1F *hT0Fill, Double_t t0fill, Int_t index)
161{
162
163 Int_t det[5];
164 Float_t length, timeexp, pos[3];
165
166 /* compute length and expected time */
167 tofGeo.GetVolumeIndices(index, det);
168 tofGeo.GetPosPar(det, pos);
169 length = 0.;
170 for (Int_t i = 0; i < 3; i++) length += pos[i] * pos[i];
171 length = TMath::Sqrt(length);
172 timeexp = length * c_1;
173
174 Double_t hittime = hT0Fill->GetRandom() - t0fill + timeexp;
175 return hittime;
176
177}
178
179/**************************************************************/
180
181TOFpid(const Char_t *filename, Int_t ipart, Int_t icharge, Int_t iipart, Bool_t dorapidityCut = kTRUE, Bool_t electronCut = kFALSE, Bool_t cutOnTPC = kFALSE, Float_t tpcsignalMin = -2., Float_t tpcsignalMax = 2., Int_t evMax = kMaxInt, Int_t startEv = 0, Bool_t mcFlag = kFALSE)
182{
183
184 printf("****************************************\n");
185 printf("RUNNING TOF PID:\n");
186 if (dorapidityCut)
187 printf("RAPIDITY-CUT: %s\n", AliPID::ParticleName(ipart));
188 else
189 printf("NO RAPIDITY-CUT\n");
190 printf("CHARGE: %s\n", chargeName[icharge]);
191 printf("PARTICLE-ID: %s\n", AliPID::ParticleName(iipart));
192 if (electronCut)
193 printf("-> ELECTRON CUT REQUESTED\n");
194 if (cutOnTPC)
195 printf(" -> CUT-ON-TPC [%3.1f,%3.1f]\n", tpcsignalMin, tpcsignalMax);
196 printf("****************************************\n");
197
198 /* include path for ACLic */
199 gSystem->AddIncludePath("-I$ALICE_ROOT/include");
200 gSystem->AddIncludePath("-I$ALICE_ROOT/TOF");
201 /* load libraries */
202 gSystem->Load("libANALYSIS");
203 gSystem->Load("libANALYSISalice");
204 /* build analysis task class */
205 gROOT->LoadMacro("AliAnalysisParticle.cxx+g");
206 gROOT->LoadMacro("AliAnalysisEvent.cxx+g");
207 gROOT->LoadMacro("AliAnalysisTrack.cxx+g");
208
209 /* create TOF response with tail */
210 // gROOT->LoadMacro("~/ALICE.2011/ANALYSIS/TOFSpectraPbPb/macros/TOFsignal.C");
211 gROOT->LoadMacro("TOFsignal.C");
212 TF1 *fTOFsignal = new TF1("fTOFsignal", TOFsignal, -2440., 2440., 4);
213 fTOFsignal->SetParameter(0, 1.);
214 fTOFsignal->SetParameter(1, 0.);
215 fTOFsignal->SetParameter(2, tofReso);
216 fTOFsignal->SetParameter(3, tofTail);
217
218 /* open file, get tree and connect */
219 TFile *filein = TFile::Open(filename);
220 TTree *treein = (TTree *)filein->Get("aodTree");
221 printf("got \"aodTree\": %d entries\n", treein->GetEntries());
222 AliAnalysisEvent *analysisEvent = new AliAnalysisEvent();
223 TClonesArray *analysisTrackArray = new TClonesArray("AliAnalysisTrack");
224 AliAnalysisTrack *analysisTrack = NULL;
225 treein->SetBranchAddress("AnalysisEvent", &analysisEvent);
226 treein->SetBranchAddress("AnalysisTrack", &analysisTrackArray);
227
228 /* open hT0fill for mismatch */
229 TFile *filein_T0Fill = TFile::Open(t0FillOnlineFileName);
230 TH1F *hT0Fill = (TH1F *)filein_T0Fill->Get("hT0Fill");
231 Double_t t0fill = t0Fill_offset;
232
233 /* open enabled flag map */
234 TH1F *hEnabledFlag = NULL;
235 if (enabledChannelsFileName) {
236 TFile *enabledfile = TFile::Open(enabledChannelsFileName);
237 hEnabledFlag = (TH1F *)enabledfile->Get("hEnabledFlag");
238 }
239
240 /**************************************************************/
241 /*** HISTOS ***************************************************/
242 /**************************************************************/
243
244 /* run-time binning */
245 for (Int_t ibin = 0; ibin < NtofsigmaBins + 1; ibin++)
246 tofsigmaBin[ibin] = tofsigmaMin + ibin * tofsigmaStep;
247 for (Int_t ibin = 0; ibin < NtofsignalBins + 1; ibin++)
248 tofsignalBin[ibin] = tofsignalMin + ibin * tofsignalStep;
249 for (Int_t ibin = 0; ibin < NtpcsigmaBins + 1; ibin++)
250 tpcsigmaBin[ibin] = tpcsigmaMin + ibin * tpcsigmaStep;
251 for (Int_t ibin = 0; ibin < NtofresoBins + 1; ibin++)
252 tofresoBin[ibin] = tofresoMin + ibin * tofresoStep;
253
254 /* histos */
255 TH1F *hAcceptedEvents = new TH1F(Form("hAcceptedEvents_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)), "", NcentralityBins, centralityBin);
256 TH2I *hAcceptedTracks = new TH2I(Form("hAcceptedTracks_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)), "", NcentralityBins, centralityBin, NptBins, ptBin);
257
258 TH3I *hTOFreso = new TH3I(Form("hTOFreso_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)), "", NcentralityBins, centralityBin, NptBins, ptBin, NtofresoBins, tofresoBin);
259
260 TH3I *hTOFpid = new TH3I(Form("hTOFpid_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)), "", NcentralityBins, centralityBin, NptBins, ptBin, NtofsigmaBins, tofsigmaBin);
261 TH3I *hTOFpid_withTZERO = new TH3I(Form("hTOFpid_withTZERO_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)), "", NcentralityBins, centralityBin, NptBins, ptBin, NtofsigmaBins, tofsigmaBin);
262 TH3I *hTOFmismatch = new TH3I(Form("hTOFmismatch_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)), "", NcentralityBins, centralityBin, NptBins, ptBin, NtofsigmaBins, tofsigmaBin);
263 TH3I *hTOFexpected[AliPID::kSPECIES];
264 for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) {
265 hTOFexpected[iiipart] = new TH3I(Form("hTOFexpected_%s_%s_%sID_%sBKG", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart), AliPID::ParticleName(iiipart)), "", NcentralityBins, centralityBin, NptBins, ptBin, NtofsigmaBins, tofsigmaBin);
266 }
267
268 TH3I *hTOFpid_delta = new TH3I(Form("hTOFpid_delta_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)), "", NcentralityBins, centralityBin, NptBins, ptBin, NtofsignalBins, tofsignalBin);
269 TH3I *hTOFmismatch_delta = new TH3I(Form("hTOFmismatch_delta_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)), "", NcentralityBins, centralityBin, NptBins, ptBin, NtofsignalBins, tofsignalBin);
270 TH3I *hTOFexpected_delta[AliPID::kSPECIES];
271 for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) {
272 hTOFexpected_delta[iiipart] = new TH3I(Form("hTOFexpected_delta_%s_%s_%sID_%sBKG", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart), AliPID::ParticleName(iiipart)), "", NcentralityBins, centralityBin, NptBins, ptBin, NtofsignalBins, tofsignalBin);
273 }
274
275 /**************************************************************/
276 /**************************************************************/
277 /**************************************************************/
278
279 /* TOF PID response */
280 AliTOFPIDResponse tofResponse;
281 tofResponse.SetTimeResolution(tofReso);
282 /* TPC PID response */
283 AliTPCPIDResponse *tpcResponse = AliAnalysisTrack::GetTPCResponse();
284
285 /* start stopwatch */
286 TStopwatch timer;
287 timer.Start();
288
289 /* verbose */
290 printf("***** RUNNING for %s %s *****\n", chargeName[icharge], AliPID::ParticleName(ipart));
291 if (!dorapidityCut)
292 printf("***** NO RAPIDITY CUT REQUESTED *****\n");
293 if (cutOnTPC) {
294 printf("***** CUT-ON-TPC requested %3.1f-%3.1f *****\n", tpcsignalMin, tpcsignalMax);
295 }
296
297 /* loop over events */
298 Bool_t hastofpid;
299 Int_t charge, index;
300 UShort_t dedxN;
301 Double_t cent, p, pt, mt, tofsignal, tpcsignal;
302 Double_t dedx, bethe, deltadedx, dedx_sigma, ptpc;
303 Double_t time, time_sigma, timezero, timezeroTOF, timezeroA, timezeroC, timezero_sigmaTOF, timezero_sigma, tof, tof_sigma, texp, texp_sigma, deltat, deltat_sigma, tof_rnd, tof_th, signal_smear, timezero_smear, texp_smear;
304
305 /* open TZERO calibfile */
306 TH1 *hCentrality_TZEROA_mean = NULL;
307 TH1 *hCentrality_TZEROA_sigma = NULL;
308 TH1 *hCentrality_TZEROC_mean = NULL;
309 TH1 *hCentrality_TZEROC_sigma = NULL;
310 TH1 *hCentrality_TOF_mean = NULL;
311 TH1 *hCentrality_TOF_TZEROA_mean = NULL;
312 TH1 *hCentrality_TOF_TZEROC_mean = NULL;
313 TH1 *hCentrality_TOF_TZEROTOF_mean = NULL;
314 TH1 *hCentrality_TZEROA_reso = NULL;
315 TH1 *hCentrality_TZEROC_reso = NULL;
316 TH1 *hCentrality_TZEROAND_reso = NULL;
317 if (1) {
318 TFile *calibfile = TFile::Open("TZEROcalibration.root");
319 hCentrality_TZEROA_mean = (TH1 *)calibfile->Get("hCentrality_TZEROA_mean");
320 hCentrality_TZEROA_sigma = (TH1 *)calibfile->Get("hCentrality_TZEROA_sigma");
321 hCentrality_TZEROC_mean = (TH1 *)calibfile->Get("hCentrality_TZEROC_calib");
322 hCentrality_TZEROC_sigma = (TH1 *)calibfile->Get("hCentrality_TZEROC_sigma");
323 hCentrality_TOF_mean = (TH1 *)calibfile->Get("hCentrality_TOF_mean");
324 hCentrality_TOF_TZEROA_mean = (TH1 *)calibfile->Get("hCentrality_TOF_TZEROA_mean");
325 hCentrality_TOF_TZEROC_mean = (TH1 *)calibfile->Get("hCentrality_TOF_TZEROC_mean");
326 hCentrality_TOF_TZEROTOF_mean = (TH1 *)calibfile->Get("hCentrality_TOF_TZEROTOF_mean");
327
328 TFile *resofile = TFile::Open("TZEROresolution.root");
329 hCentrality_TZEROA_reso = (TH1 *)resofile->Get("hTZEROA_reso");
330 hCentrality_TZEROC_reso = (TH1 *)resofile->Get("hTZEROC_reso");
331 hCentrality_TZEROAND_reso = (TH1 *)resofile->Get("hTZEROAND_reso");
332 }
333 Double_t TZEROA_mean;
334 Double_t TZEROA_sigma;
335 Double_t TZEROC_mean;
336 Double_t TZEROC_sigma;
337 Double_t TOF_mean;
338 Double_t TOF_TZEROA_mean;
339 Double_t TOF_TZEROC_mean;
340 Double_t TOF_TZEROTOF_mean;
341 Double_t TZEROA;
342 Double_t TZEROA_reso;
343 Bool_t hasTZEROA;
344 Double_t TZEROC;
345 Double_t TZEROC_reso;
346 Bool_t hasTZEROC;
347 Double_t TZEROAND;
348 Double_t TZEROAND_reso;
349 Bool_t hasTZEROAND;
350 Double_t TZEROTOF;
351 Double_t TZEROTOF_reso;
352 Bool_t hasTZEROTOF;
353 Double_t TZEROMEAN;
354 Double_t TZEROMEAN_weight;
355 Double_t TZEROBEST;
356 Double_t TZEROBEST_reso;
357
358 Double_t param[kNHistoParams];
359 for (Int_t iev = startEv; iev < treein->GetEntries() && iev < evMax; iev++) {
360 /* get event */
361 treein->GetEvent(iev);
362 if (iev % 10000 == 0) printf("iev = %d\n", iev);
363 /* check event */
364 if (!analysisEvent->AcceptEvent(acceptEventType)) continue;
365
366 /*** ACCEPTED EVENT ***/
367
368 /* apply time-zero TOF correction */
369 // analysisEvent->ApplyTimeZeroTOFCorrection();
370
371 /* get centrality */
372 cent = analysisEvent->GetCentralityPercentile(centralityEstimator);
373
374 /* vertex */
375 Double_t vertexz = analysisEvent->GetVertexZ();
376
377 /* fill histos */
378 hAcceptedEvents->Fill(cent);
379
380 /* TZERO corrections */
381 Int_t icent;
382 for (icent = 0; icent < NcentralityBins; icent++)
383 if (cent < centralityBin[icent + 1])
384 break;
385 TZEROA_mean = hCentrality_TZEROA_mean ? hCentrality_TZEROA_mean->GetBinContent(icent + 1) : 0.;
386 TZEROA_sigma = hCentrality_TZEROA_sigma ? hCentrality_TZEROA_sigma->GetBinContent(icent + 1) : 1000.;
387 TZEROC_mean = hCentrality_TZEROC_mean ? hCentrality_TZEROC_mean->GetBinContent(icent + 1) : 0.;
388 TZEROC_sigma = hCentrality_TZEROC_sigma ? hCentrality_TZEROC_sigma->GetBinContent(icent + 1) : 1000.;
389
390 TOF_mean = hCentrality_TOF_mean ? hCentrality_TOF_mean->GetBinContent(icent + 1) : 0.;
391 TOF_TZEROA_mean = hCentrality_TOF_TZEROA_mean ? hCentrality_TOF_TZEROA_mean->GetBinContent(icent + 1) : 0.;
392 TOF_TZEROC_mean = hCentrality_TOF_TZEROC_mean ? hCentrality_TOF_TZEROC_mean->GetBinContent(icent + 1) : 0.;
393 TOF_TZEROTOF_mean = hCentrality_TOF_TZEROTOF_mean ? hCentrality_TOF_TZEROTOF_mean->GetBinContent(icent + 1) : 0.;
394
395 TZEROA_reso = hCentrality_TZEROA_reso ? hCentrality_TZEROA_reso->GetBinContent(icent + 1) : 70.;
396 TZEROC_reso = hCentrality_TZEROC_reso ? hCentrality_TZEROC_reso->GetBinContent(icent + 1) : 70.;
397 TZEROAND_reso = hCentrality_TZEROAND_reso ? hCentrality_TZEROAND_reso->GetBinContent(icent + 1) : 50.;
398
399 /* simulate inefficient TZERO */
400 Bool_t resetTZERO = kFALSE;
401 if (forcetimezeroineff > 0.)
402 if (gRandom->Uniform() < forcetimezeroineff)
403 resetTZERO = kTRUE;
404
405 /* loop over tracks */
406 for (Int_t itrk = 0; itrk < analysisTrackArray->GetEntries(); itrk++) {
407 /* get track */
408 analysisTrack = (AliAnalysisTrack *)analysisTrackArray->At(itrk);
409 if (!analysisTrack) continue;
410 /* check accepted track */
411 if (!analysisTrack->AcceptTrack()) continue;
412 /* check rapidity */
413 if (dorapidityCut)
414 if (analysisTrack->GetY(AliPID::ParticleMass(ipart)) - rapidityShift > rapidityMaxCut ||
415 analysisTrack->GetY(AliPID::ParticleMass(ipart)) - rapidityShift < rapidityMinCut) continue;
416
417 /* check charge */
418 charge = analysisTrack->GetSign() > 0. ? kPositive : kNegative;
419 if (charge != icharge) continue;
420
421 /*** ACCEPTED TRACK ***/
422
423 /* get track info */
424 p = analysisTrack->GetP();
425 pt = analysisTrack->GetPt();
426
427 /* compute track mt */
428 mt = TMath::Sqrt(pt * pt + AliPID::ParticleMass(ipart) * AliPID::ParticleMass(ipart));
429
430 /* get TPC info */
431 dedx = analysisTrack->GetTPCdEdx();
432 dedxN = analysisTrack->GetTPCdEdxN();
433 ptpc = analysisTrack->GetTPCmomentum();
434
435 /* TPC signal */
436 bethe = tpcResponse->GetExpectedSignal(ptpc, iipart);
437 /* fix electron signal */
438 if (iipart == AliPID::kElectron)
439 bethe += 23.;
440 deltadedx = dedx - bethe;
441 dedx_sigma = 0.07 * bethe;
442 tpcsignal = deltadedx / dedx_sigma;
443
444 /* electronCut requested, remove electrons */
445 if (electronCut) {
446 /* TPC signal */
447 bethe = tpcResponse->GetExpectedSignal(ptpc, AliPID::kElectron);
448 /* fix electron signal */
449 bethe += 23.;
450 deltadedx = dedx - bethe;
451 dedx_sigma = 0.07 * bethe;
452 tpcsignal = deltadedx / dedx_sigma;
453 if (TMath::Abs(tpcsignal) < 1.5) continue;
454 }
455
456 /* cut on TPC signal if requested */
457 if (cutOnTPC && (tpcsignal < tpcsignalMin || tpcsignal > tpcsignalMax))
458 continue;
459
460 /* fill histos */
461 hAcceptedTracks->Fill(cent, pt);
462
463 /* set TOF pid flag */
464 hastofpid = analysisTrack->HasTOFPID(hEnabledFlag);
465
466 /* check TOF pid */
467 if (!hastofpid)
468 continue;
469
470 /*** ACCEPTED TRACK WITH TOF PID ***/
471
472 /* apply expected time correction */
473 // analysisTrack->ApplyTOFExpectedTimeCorrection();
474
475 /* get TOF info */
476 index = analysisTrack->GetTOFIndex();
477 time = analysisTrack->GetTOFTime() - TOF_mean;
478 time_sigma = tofReso;
479 /* TZEROTOF */
480 TZEROTOF = analysisEvent->GetTimeZeroTOF(analysisTrack->GetP());
481 TZEROTOF_reso = analysisEvent->GetTimeZeroTOFSigma(analysisTrack->GetP());
482 hasTZEROTOF = TZEROTOF_reso < 150.;
483 if (hasTZEROTOF) {
484 TZEROTOF += TOF_TZEROTOF_mean - TOF_mean;
485 TZEROTOF_reso *= TZEROTOF_resoScaleFactor;
486 }
487 /* TZEROA */
488 TZEROA = analysisEvent->GetTimeZeroT0(1) - TZEROA_shift;
489 // TZEROA_reso = TZEROA_resolution[icent];
490 hasTZEROA = TMath::Abs(TZEROA - TZEROA_mean) < 3. * TZEROA_sigma;
491 TZEROA += TOF_TZEROA_mean - TOF_mean;
492 /* TZEROC */
493 TZEROC = analysisEvent->GetTimeZeroT0(2) - TZEROC_shift;
494 // TZEROC_reso = TZEROC_resolution[icent];
495 hasTZEROC = TMath::Abs(TZEROC - TZEROC_mean) < 3. * TZEROC_sigma;
496 TZEROC += TOF_TZEROC_mean - TOF_mean;
497 /* TZEROAND */
498 TZEROAND = (TZEROA + TZEROC) * 0.5;
499 // TZEROAND_reso = TZEROAND_resolution[icent];
500 hasTZEROAND = hasTZEROA && hasTZEROC;
501 /* TZEROMEAN */
502 TZEROMEAN = TZEROTOF / TZEROTOF_reso / TZEROTOF_reso;
503 TZEROMEAN_weight = 1. / TZEROTOF_reso / TZEROTOF_reso;
504 if (hasTZEROAND) {
505 // printf("TZEROAND\n");
506 TZEROMEAN += TZEROAND / TZEROAND_reso / TZEROAND_reso;
507 TZEROMEAN_weight = 1. / TZEROAND_reso / TZEROAND_reso;
508 }
509 else if (hasTZEROA) {
510 // printf("TZEROA\n");
511 TZEROMEAN += TZEROA / TZEROA_reso / TZEROA_reso;
512 TZEROMEAN_weight = 1. / TZEROA_reso / TZEROA_reso;
513 }
514 else if (hasTZEROC) {
515 // printf("TZEROC\n");
516 TZEROMEAN += TZEROC / TZEROC_reso / TZEROC_reso;
517 TZEROMEAN_weight = 1. / TZEROC_reso / TZEROC_reso;
518 }
519 timezero = TZEROMEAN / TZEROMEAN_weight;
520 timezero_sigma = TMath::Sqrt(1. / TZEROMEAN_weight);
521 /* TZEROBEST */
522 TZEROBEST = TZEROTOF;
523 TZEROBEST_reso = TZEROTOF_reso;
524 if (hasTZEROAND && TZEROAND_reso < TZEROBEST_reso) {
525 TZEROBEST = TZEROAND;
526 TZEROBEST_reso = TZEROAND_reso;
527 }
528 else if (hasTZEROA && TZEROA_reso < TZEROBEST_reso) {
529 TZEROBEST = TZEROA;
530 TZEROBEST_reso = TZEROA_reso;
531 }
532 if (hasTZEROC && TZEROC_reso < TZEROBEST_reso) {
533 TZEROBEST = TZEROC;
534 TZEROBEST_reso = TZEROC_reso;
535 }
536 timezero = TZEROBEST;
537 timezero_sigma = TZEROBEST_reso;
538
539 /* DEBUG */
540 // timezero = 0.;//TZEROTOF;
541 // timezero_sigma = 203.854691;//TZEROTOF_reso;
542
543 // if (timezero == 0.)
544 // printf("%f %f\n", timezero, timezero_sigma);
545
546 timezero_sigma *= scaletimezerosigma;
547
548 if (resetTZERO) {
549 timezero = 0.;
550 timezero_sigma = timezero_spread;
551 }
552
553
554 tof = time - timezero;
555 tof_sigma = TMath::Sqrt(time_sigma * time_sigma + timezero_sigma * timezero_sigma);
556
557 /* TOF expected time */
558 texp = analysisTrack->GetTOFExpTime(iipart);
559 texp_sigma = analysisTrack->GetTOFExpTimeSigma(iipart) * scaletexpreso[iipart];
560
561 /* TOF signal */
562 deltat = tof - texp;
563 deltat_sigma = TMath::Sqrt(tof_sigma * tof_sigma + texp_sigma * texp_sigma);
564 tofsignal = deltat / deltat_sigma;
565
566 /* fill histo */
567 hTOFpid->Fill(cent, pt, tofsignal);
568 hTOFpid_delta->Fill(cent, p, deltat);
569 hTOFreso->Fill(cent, pt, deltat_sigma);
570 if (hasTZEROTOF || hasTZEROA || hasTZEROC || hasTZEROAND)
571 hTOFpid_withTZERO->Fill(cent, pt, tofsignal);
572
573
574 /*** EXPECTED MISMATCH ***/
575
576 /* loop to generated random hits */
577 for (Int_t irnd = 0; irnd < NmismatchTrials; irnd++) {
578
579 /* generate ramdom tof values */
580 tof_rnd = GenerateRandomHit(hT0Fill, t0fill, index);
581
582 /* TOF signal */
583 deltat = tof_rnd - texp;
584 tofsignal = deltat / deltat_sigma;
585
586 /* fill histo */
587 hTOFmismatch->Fill(cent, pt, tofsignal);
588 hTOFmismatch_delta->Fill(cent, p, deltat);
589
590 } /* end of loop over generated random hits */
591
592 /*** EXPECTED SIGNALS ***/
593
594 /* loop over other particles */
595 for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) {
596
597 /* generate expected tof value */
598 tof_th = analysisTrack->GetTOFExpTime(iiipart);
599 texp_sigma = analysisTrack->GetTOFExpTimeSigma(iiipart) * scaletexpreso[iiipart];
600
601 /* loop over many trials */
602 for (Int_t irnd = 0; irnd < NexpectedTrials; irnd++) {
603
604 /* tof response smearing */
605 signal_smear = fTOFsignal->GetRandom();
606 /* timezero resolution smearing */
607 timezero_smear = gRandom->Gaus(0., timezero_sigma);
608 /* texp resolution smearing */
609 texp_smear = gRandom->Gaus(0., texp_sigma);
610
611 /* deltat and sigma */
612 deltat = tof_th - texp + signal_smear + timezero_smear + texp_smear;
613 tofsignal = deltat / deltat_sigma;
614
615 /* fill histo */
616 hTOFexpected[iiipart]->Fill(cent, pt, tofsignal);
617 hTOFexpected_delta[iiipart]->Fill(cent, p, deltat);
618
619 } /* end of loop over many trials */
620 } /* end of loop over other particle */
621 } /* end of loop over tracks */
622 } /* end of loop over events */
623
624 /* stop stopwatch */
625 timer.Stop();
626 timer.Print();
627
628 /* output */
629 TString outputstring = "TOFpid";
630 if (!dorapidityCut)
631 outputstring += "_noRapidityCut";
632 if (electronCut)
633 outputstring += "_electronCut";
634 if (cutOnTPC)
635 outputstring += Form("_cutOnTPC[%3.1f,%3.1f]", , tpcsignalMin, tpcsignalMax);
636 outputstring += Form("_%s_%s_%sID.%s", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart), filename);
637 TFile *fileout = TFile::Open(outputstring.Data(), "RECREATE");
638 hAcceptedEvents->Write();
639 hAcceptedTracks->Write();
640 hTOFpid->Write();
641 hTOFreso->Write();
642 hTOFpid_withTZERO->Write();
643 hTOFmismatch->Write();
644 for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++)
645 hTOFexpected[iiipart]->Write();
646 hTOFpid_delta->Write();
647 hTOFmismatch_delta->Write();
648 for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++)
649 hTOFexpected_delta[iiipart]->Write();
650
651 fileout->Close();
652
653}
654
655//___________________________________________________________________________________
656
657TOFspectra_defaultFit(const Char_t *filename)
658{
659
660Bool_t EXPECTED_SIGNAL_TEMPLATE = kTRUE;
661Bool_t EXPECTED_SIGNAL_FIT = kFALSE;
662Bool_t EXPECTED_BACKGROUND_TEMPLATE = kFALSE;
663Bool_t EXPECTED_BACKGROUND_FIT = kTRUE;
664
665}
666
667//___________________________________________________________________________________
668
669TOFspectra_defaultFit_fitElectrons(const Char_t *filename, Float_t electronLimit = 5.)
670{
671
672
673 TOFspectra(filename, electronLimit);
674}
675
676//___________________________________________________________________________________
677
678TOFspectra_signalFit(Bool_t fixParams = kTRUE, Float_t scaleSigma = 1., Float_t scaleTail = 1.)
679{
680
681 EXPECTED_SIGNAL_TEMPLATE = kFALSE;
682 EXPECTED_SIGNAL_FIT = kTRUE;
683 FIX_SIGNAL_MEAN = fixParams;
684 FIX_SIGNAL_SIGMA = fixParams;
685 FIX_SIGNAL_TAIL = fixParams;
686 SCALE_SIGNAL_SIGMA = scaleSigma;
687 SCALE_SIGNAL_TAIL = scaleTail;
688
689}
690
691//___________________________________________________________________________________
692
693TOFspectra_bkgFit(Bool_t fixParams = kTRUE, Float_t scaleSigma = 1., Float_t scaleTail = 1.)
694{
695
696 EXPECTED_BACKGROUND_TEMPLATE = kFALSE;
697 EXPECTED_BACKGROUND_FIT = kTRUE;
698 FIX_BACKGROUND_MEAN = fixParams;
699 FIX_BACKGROUND_SIGMA = fixParams;
700 FIX_BACKGROUND_TAIL = fixParams;
701 SCALE_BACKGROUND_SIGMA = scaleSigma;
702 SCALE_BACKGROUND_TAIL = scaleTail;
703
704 TOFspectra(filename);
705
706}
707
708//___________________________________________________________________________________
709
710void
711TOFspectra(const Char_t *filename, Float_t electronLimit = 0.)
712{
713
714 for (Int_t icent = 0; icent < NcentralityBins; icent++)
715 for (Int_t icharge = 0; icharge < kNCharges; icharge++)
716 for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++)
717 TOFspectrum(filename, ipart, icharge, ipart, icent, -1., -1., electronLimit);
718
719}
720
721//___________________________________________________________________________________
722
723#define USENSIGMA 1
724
725/* fit ranges */
726Double_t fitPtMin[AliPID::kSPECIES] = {0.5, 0.5, 0.3, 0.4, 0.5};
727Double_t fitPtMax[AliPID::kSPECIES] = {3.0, 3.0, 3.0, 3.0, 5.0};
728#if USENSIGMA
729Double_t fitSigmaMin[AliPID::kSPECIES] = {tofsigmaMin, tofsigmaMin, -25., -75., -65.};
730Double_t fitSigmaMax[AliPID::kSPECIES] = {tofsigmaMax, tofsigmaMax, 225., 200., 100.};
731#else
732Double_t fitSigmaMin[AliPID::kSPECIES] = {-24400., -24400., -24400., -24400., -24400.};
733Double_t fitSigmaMax[AliPID::kSPECIES] = {24400., 24400., 24400., 24400., 24400.};
734#endif
735
736void
737TOFspectrum(const Char_t *filename, Int_t ipart, Int_t icharge, Int_t iipart, Int_t icent, Float_t ptMin = -1., Float_t ptMax = -1., Bool_t checkHistoFlag = kFALSE)
738{
739
740 printf("****************************************\n");
741 printf("RUNNING SPECTRA FIT:\n");
742 printf("RAPIDITY-CUT: %s\n", AliPID::ParticleName(ipart));
743 printf("RAPIDITY-CUT: %s\n", AliPID::ParticleName(ipart));
744 printf("CHARGE: %s\n", chargeName[icharge]);
745 printf("PARTICLE: %s\n", AliPID::ParticleName(iipart));
746 printf("CENTRALITY BIN: %d\n", icent);
747 printf("****************************************\n");
748
749 /* mismatch balance functions */
750 TF1 *fMismatchBalanceFunction[5][2];
751 Double_t mismatchBalanceParam0[5][2] = {
752 0., 0.,
753 0., 0.,
754 0.669488, 0.599374,
755 -1.2582, -2.53613,
756 5.48850e-01, 10.1622
757 };
758 Double_t mismatchBalanceParam1[5][2] = {
759 0., 0.,
760 0., 0.,
761 -1.64894, -0.825764,
762 -1.63683, -2.55156,
763 -1.64894, -7.09531
764 };
765 for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) {
766 for (Int_t iiicharge = 0; iiicharge < 2; iiicharge++) {
767 fMismatchBalanceFunction[iiipart][iiicharge] = new TF1(Form(""), "1. - [0]*TMath::Exp(x*[1])", 0., 5.);
768 fMismatchBalanceFunction[iiipart][iiicharge]->SetParameter(0, mismatchBalanceParam0[iiipart][iiicharge]);
769 fMismatchBalanceFunction[iiipart][iiicharge]->SetParameter(1, mismatchBalanceParam1[iiipart][iiicharge]);
770 }
771 }
772
773 /* open data */
774 TFile *filein = TFile::Open(filename);
775
776 /* get number of events */
777 TH1F *hAcceptedEvents = (TH1F *)filein->Get(Form("hAcceptedEvents_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)));
778 if (!hAcceptedEvents) {
779 printf("cannot find %s\n", Form("hAcceptedEvents_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)));
780 return;
781 }
782 Double_t nevents;
783 if (icent < 0 || icent >= NcentralityBins)
784 nevents = hAcceptedEvents->GetEntries();
785 else
786 nevents = hAcceptedEvents->Integral(icent + 1, icent + 1);
787 printf("N_EVENTS : %d\n", nevents);
788 printf("****************************************\n");
789
790 /* get histos */
791 TH2I *hAcceptedTracks = (TH2I *)filein->Get(Form("hAcceptedTracks_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)));
792 if (!hAcceptedTracks) {
793 printf("cannot find %s\n", Form("hAcceptedTracks_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)));
794 // return;
795 }
796#if USENSIGMA
797 TH3I *hTOFpid = (TH3I *)filein->Get(Form("hTOFpid_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)));
798#else
799 TH3I *hTOFpid = (TH3I *)filein->Get(Form("hTOFpid_delta_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)));
800#endif
801 if (!hTOFpid) {
802 printf("cannot find %s\n", Form("hTOFpid_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)));
803 return;
804 }
805#if USENSIGMA
806 TH3I *hTOFmismatch = (TH3I *)filein->Get(Form("hTOFmismatch_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)));
807#else
808 TH3I *hTOFmismatch = (TH3I *)filein->Get(Form("hTOFmismatch_delta_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)));
809#endif
810 if (!hTOFmismatch) {
811 printf("cannot find %s\n", Form("hTOFpid_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)));
812 return;
813 }
814 TH3I *hTOFexpected[AliPID::kSPECIES];
815 for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) {
816#if USENSIGMA
817 hTOFexpected[iiipart] = (TH3I *)filein->Get(Form("hTOFexpected_%s_%s_%sID_%sBKG", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart), AliPID::ParticleName(iiipart)));
818#else
819 hTOFexpected[iiipart] = (TH3I *)filein->Get(Form("hTOFexpected_delta_%s_%s_%sID_%sBKG", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart), AliPID::ParticleName(iiipart)));
820#endif
821 if (!hTOFexpected[iiipart]) {
822 printf("cannot find %s\n", Form("hTOFexpected_%s_%s_%sID_%sBKG", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart), AliPID::ParticleName(iiipart)));
823 return;
824 }
825 }
826
827 /* setup centrality range */
828 if (icent < 0 || icent >= NcentralityBins) {
829 printf("WARNING: undefined centrality -> using 00-90\% range\n");
830 if (hAcceptedTracks) hAcceptedTracks->GetXaxis()->SetRange(1, NcentralityBins);
831 hTOFpid->GetXaxis()->SetRange(1, NcentralityBins);
832 hTOFmismatch->GetXaxis()->SetRange(1, NcentralityBins);
833 for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++)
834 hTOFexpected[iiipart]->GetXaxis()->SetRange(1, NcentralityBins);
835 }
836 else {
837 printf("***** FITTING CENTRALITY-BIN [%02d, %02d] %% *****\n", centralityBin[icent], centralityBin[icent + 1]);
838 if (hAcceptedTracks) hAcceptedTracks->GetXaxis()->SetRange(icent + 1, icent + 1);
839 hTOFpid->GetXaxis()->SetRange(icent + 1, icent + 1);
840 hTOFmismatch->GetXaxis()->SetRange(icent + 1, icent + 1);
841 for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++)
842 hTOFexpected[iiipart]->GetXaxis()->SetRange(icent + 1, icent + 1);
843 }
844
845 /* init flags */
846 Bool_t requestedRange = kFALSE;
847 Bool_t fitElectrons = kTRUE;
848 Bool_t fitMuons = kTRUE;
849 Bool_t fitPions = kTRUE;
850
851 /* setup pt range if requested */
852 if (ptMin > -0.001 && ptMax > -0.001 && ptMax > ptMin) {
853 printf("***** FITTING PT-BIN [%f, %f] GeV/c *****\n", ptMin, ptMax);
854 requestedRange = kTRUE;
855
856 /* check electron-fit is allowed */
857 fitElectrons = kTRUE;
858 if ((ptMin + 0.001) < FIT_ELECTRONS_PT_MIN || (ptMax - 0.001) > FIT_ELECTRONS_PT_MAX)
859 fitElectrons = kFALSE;
860 /* check muon-fit is allowed */
861 fitMuons = kTRUE;
862 if ((ptMin + 0.001) < FIT_MUONS_PT_MIN || (ptMax - 0.001) > FIT_MUONS_PT_MAX)
863 fitMuons = kFALSE;
864 /* check pion-fit is allowed */
865 fitPions = kTRUE;
866 if ((ptMin + 0.001) < FIT_PIONS_PT_MIN || (ptMax - 0.001) > FIT_PIONS_PT_MAX)
867 fitPions = kFALSE;
868
869
870 hTOFpid->GetYaxis()->SetRangeUser(ptMin + 0.001, ptMax - 0.001);
871 hTOFmismatch->GetYaxis()->SetRangeUser(ptMin + 0.001, ptMax - 0.001);
872 for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++)
873 hTOFexpected[iiipart]->GetYaxis()->SetRangeUser(ptMin + 0.001, ptMax - 0.001);
874 }
875
876 /* output */
877 Char_t outfilename[1024];
878 if (icent < 0 || icent >= NcentralityBins)
879 sprintf(outfilename, "TOFspectrum_cent0090_%s_%s_%sID.root", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart));
880 else {
881 sprintf(outfilename, "TOFspectrum_cent%02d%02d_%s_%s_%sID.root", centralityBin[icent], centralityBin[icent + 1], AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart));
882 }
883 TFile *fileout = TFile::Open(outfilename, "RECREATE");
884 TDirectory *fitDir = fileout->mkdir("FitParams");
885 /* canvas */
886 TCanvas *canvas = new TCanvas("canvas");
887 canvas->SetLogy();
888 /* histo */
889 TH1D *hFitParamHisto[kNFitParams];
890 for (Int_t iparam = 0; iparam < kNFitParams; iparam++)
891 hFitParamHisto[iparam] = new TH1D(Form("h%s", fitParamName[iparam]), fitParamTitle[iparam], NptBins, ptBin);
892
893 /* loop over ptBins */
894 for (Int_t ipt = 0; ipt < NptBins; ipt++) {
895
896 if (!requestedRange) {
897 if ((ptBin[ipt] + 0.001) < fitPtMin[ipart] || (ptBin[ipt + 1] - 0.001) > fitPtMax[ipart]) continue;
898 printf("***** FITTING PT-BIN [%f, %f] GeV/c *****\n", ptBin[ipt], ptBin[ipt + 1]);
899
900 /* check electron-fit is allowed */
901 fitElectrons = kTRUE;
902 if ((ptBin[ipt] + 0.001) < FIT_ELECTRONS_PT_MIN || (ptBin[ipt + 1] - 0.001) > FIT_ELECTRONS_PT_MAX)
903 fitElectrons = kFALSE;
904 /* check muon-fit is allowed */
905 fitMuons = kTRUE;
906 if ((ptBin[ipt] + 0.001) < FIT_MUONS_PT_MIN || (ptBin[ipt + 1] - 0.001) > FIT_MUONS_PT_MAX)
907 fitMuons = kFALSE;
908 /* check pion-fit is allowed */
909 fitPions = kTRUE;
910 if ((ptBin[ipt] + 0.001) < FIT_PIONS_PT_MIN || (ptBin[ipt + 1] - 0.001) > FIT_PIONS_PT_MAX)
911 fitPions = kFALSE;
912
913 hTOFpid->GetYaxis()->SetRange(ipt + 1, ipt + 1);
914 hTOFmismatch->GetYaxis()->SetRange(ipt + 1, ipt + 1);
915 for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++)
916 hTOFexpected[iiipart]->GetYaxis()->SetRange(ipt + 1, ipt + 1);
917 }
918
919 /* nsigma projections */
920 TH1 *hSignal_py = hTOFpid->Project3D("z");
921 TH1 *hMismatch_py = hTOFmismatch->Project3D("z");
922 TH1 *hSignalExp_py[AliPID::kSPECIES];
923 TH1 *hSignalExpTail_py[AliPID::kSPECIES];
924 for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) {
925 hSignalExp_py[iiipart] = hTOFexpected[iiipart]->Project3D("z");
926 hSignalExpTail_py[iiipart] = hTOFexpected[iiipart]->Project3D("z");
927 }
928
929 /* prepare histos for the fitter */
930 Int_t partbkg1[AliPID::kSPECIES] = {AliPID::kKaon, AliPID::kKaon, AliPID::kKaon, AliPID::kPion, AliPID::kPion};
931 Int_t partbkg2[AliPID::kSPECIES] = {AliPID::kProton, AliPID::kProton, AliPID::kProton, AliPID::kProton, AliPID::kKaon};
932 Int_t partbkg3[AliPID::kSPECIES] = {AliPID::kPion, AliPID::kPion, AliPID::kElectron, AliPID::kElectron, AliPID::kElectron};
933 Int_t partbkg4[AliPID::kSPECIES] = {AliPID::kMuon, AliPID::kElectron, AliPID::kMuon, AliPID::kMuon, AliPID::kMuon};
934 TH1 *hSigExp_py, *hBkgExp1_py, *hBkgExp2_py, *hBkgExp3_py, *hBkgExp4_py;
935 hSigExp_py = EXPECTED_SIGNAL_TAIL ? hSignalExpTail_py[iipart] : hSignalExp_py[iipart];
936 hBkgExp1_py = EXPECTED_BACKGROUND_TAIL ? hSignalExpTail_py[partbkg1[iipart]] : hSignalExp_py[partbkg1[iipart]];
937 hBkgExp2_py = EXPECTED_BACKGROUND_TAIL ? hSignalExpTail_py[partbkg2[iipart]] : hSignalExp_py[partbkg2[iipart]];
938 hBkgExp3_py = EXPECTED_BACKGROUND_TAIL ? hSignalExpTail_py[partbkg3[iipart]] : hSignalExp_py[partbkg3[iipart]];
939 hBkgExp4_py = EXPECTED_BACKGROUND_TAIL ? hSignalExpTail_py[partbkg4[iipart]] : hSignalExp_py[partbkg4[iipart]];
940
941 /* check histos if requested */
942 if (checkHistoFlag) {
943 TCanvas *cCheckHisto = new TCanvas("cCheckHisto");
944 cCheckHisto->Divide(2, 3);
945 cCheckHisto->cd(1);
946 hSignal_py->Draw();
947 cCheckHisto->cd(2);
948 hSigExp_py->Draw();
949 cCheckHisto->cd(3);
950 hBkgExp1_py->Draw();
951 cCheckHisto->cd(4);
952 hBkgExp2_py->Draw();
953 cCheckHisto->cd(5);
954 hBkgExp3_py->Draw();
955 cCheckHisto->cd(6);
956 hMismatch_py->Draw();
957 return;
958 }
959
960 Double_t rangeMin = fitSigmaMin[iipart], rangeMax = fitSigmaMax[iipart];
961 Bool_t constrainSignal = kFALSE;
962 Bool_t constrainBkg1 = kFALSE;
963 Bool_t constrainBkg2 = kFALSE;
964 Bool_t forceGaussianSignal = kFALSE;
965 Bool_t fitBkg1 = kTRUE, fitBkg2 = kTRUE, fitBkg3 = kTRUE, fitBkg4 = kTRUE;
966
967 /* check whether we can fit electrons */
968 if (!fitElectrons) {
969 printf("INHIBIT FIT ELECTRONS\n");
970 if (partbkg1[iipart] == AliPID::kElectron) fitBkg1 = kFALSE;
971 if (partbkg2[iipart] == AliPID::kElectron) fitBkg2 = kFALSE;
972 if (partbkg3[iipart] == AliPID::kElectron) fitBkg3 = kFALSE;
973 if (partbkg4[iipart] == AliPID::kElectron) fitBkg4 = kFALSE;
974 }
975 /* check whether we can fit muons */
976 if (!fitMuons) {
977 printf("INHIBIT FIT MUONS\n");
978 if (partbkg1[iipart] == AliPID::kMuon) fitBkg1 = kFALSE;
979 if (partbkg2[iipart] == AliPID::kMuon) fitBkg2 = kFALSE;
980 if (partbkg3[iipart] == AliPID::kMuon) fitBkg3 = kFALSE;
981 if (partbkg4[iipart] == AliPID::kMuon) fitBkg4 = kFALSE;
982 }
983 /* check whether we can fit pions */
984 if (!fitPions) {
985 printf("INHIBIT FIT PIONS\n");
986 if (partbkg1[iipart] == AliPID::kPion) fitBkg1 = kFALSE;
987 if (partbkg2[iipart] == AliPID::kPion) fitBkg2 = kFALSE;
988 if (partbkg3[iipart] == AliPID::kPion) fitBkg3 = kFALSE;
989 if (partbkg4[iipart] == AliPID::kPion) fitBkg4 = kFALSE;
990 }
991
992
993 /* fit */
994 Double_t param[kNFitParams];
995 Double_t param_err[kNFitParams];
996 TOFpid_fit(hSignal_py, hSigExp_py, hBkgExp1_py, hBkgExp2_py, hBkgExp3_py, hBkgExp4_py, hMismatch_py, rangeMin, rangeMax, fitBkg1, fitBkg2, fitBkg3, fitBkg4, constrainSignal, constrainBkg1, constrainBkg2, forceGaussianSignal, param, param_err, canvas);
997
998
999 /* if fitting pions add electrons and to the signal */
1000 if (iipart == 2) {
1001 param[kSignalCounts] += param[kBkg3Counts] + param[kBkg4Counts];
1002 param_err[kSignalCounts] = TMath::Sqrt(param_err[kSignalCounts] * param_err[kSignalCounts] + param_err[kBkg3Counts] * param_err[kBkg3Counts] + param_err[kBkg4Counts] * param_err[kBkg3Counts]);
1003 }
1004 /* else add to bkg1 (pions) */
1005 else {
1006 param[kBkg1Counts] += param[kBkg3Counts] + param[kBkg4Counts];
1007 param_err[kBkg1Counts] = TMath::Sqrt(param_err[kBkg1Counts] * param_err[kBkg1Counts] + param_err[kBkg3Counts] * param_err[kBkg3Counts] + param_err[kBkg4Counts] * param_err[kBkg3Counts]);
1008 }
1009
1010 /* check requested pt-range */
1011 if (requestedRange)
1012 return;
1013
1014 /* write canvas */
1015 fitDir->cd();
1016 canvas->Write(Form("fitDisplay_ptBin_%3.2f_%3.2f", ptBin[ipt], ptBin[ipt + 1]));
1017
1018 /* set histo */
1019 for (Int_t iparam = 0; iparam < kNFitParams; iparam++) {
1020 hFitParamHisto[iparam]->SetBinContent(ipt + 1, param[iparam]);
1021 hFitParamHisto[iparam]->SetBinError(ipt + 1, param_err[iparam]);
1022 }
1023
1024 /* delete */
1025 delete hSignal_py;
1026 delete hMismatch_py;
1027 for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) {
1028 delete hSignalExp_py[iiipart];
1029 // delete hSignalExpTail_py[iiipart];
1030 }
1031 }
1032
1033 /* check requested pt-range */
1034 if (requestedRange)
1035 return;
1036
1037 /*** POST-ANALYSIS ***/
1038
1039 TDirectory *postDir = fileout->mkdir("PostAnalysis");
1040
1041 /* compute overflows */
1042 TH1D *hOverflowCounts = new TH1D(*hFitParamHisto[kTotalCounts]);
1043 hOverflowCounts->SetNameTitle("hOverflowCounts", "Overflow counts: TotalCounts - IntegralCounts;p_{T} (GeV/c);counts");
1044 hOverflowCounts->Add(hFitParamHisto[kIntegralCounts], -1.);
1045 /* compute total mismatch counts */
1046 TH1D *hTotalMismatchCounts = new TH1D(*hFitParamHisto[kMismatchCounts]);
1047 hTotalMismatchCounts->SetNameTitle("hTotalMismatchCounts", "Total mismatch counts: MismatchCounts + OverflowCounts;p_{T} (GeV/c);counts");
1048 hTotalMismatchCounts->Add(hOverflowCounts);
1049 /* computed mismatch fraction */
1050 TH1D *hTotalMismatchFraction = new TH1D(*hTotalMismatchCounts);
1051 hTotalMismatchFraction->SetNameTitle("hTotalMismatchFraction", "Total mismatch fraction: TotalMismatchCounts / TotalCounts;p_{T} (GeV/c);");
1052 hTotalMismatchFraction->Divide(hFitParamHisto[kTotalCounts]);
1053 /* compute identified counts */
1054 TH1D *hIdentifiedCounts = new TH1D(*hFitParamHisto[kSignalCounts]);
1055 hIdentifiedCounts->SetNameTitle("hIdentifiedCounts", "Identified counts: SignalCounts + sum(BkgCounts);p_{T} (GeV/c);counts");
1056 hIdentifiedCounts->Add(hFitParamHisto[kBkg1Counts]);
1057 hIdentifiedCounts->Add(hFitParamHisto[kBkg2Counts]);
1058 // hIdentifiedCounts->Add(hFitParamHisto[kBkg3Counts]);
1059 // hIdentifiedCounts->Add(hFitParamHisto[kBkg4Counts]);
1060 /* compute signal fraction */
1061 TH1D *hSignalFraction = new TH1D(*hFitParamHisto[kSignalCounts]);
1062 hSignalFraction->SetNameTitle("hSignalFraction", "Signal fraction: SignalCounts / IdentifiedCounts;p_{T} (GeV/c);");
1063 hSignalFraction->Divide(hSignalFraction, hIdentifiedCounts, 1., 1., "B");
1064 /* compute bkg1 fraction */
1065 TH1D *hBkg1Fraction = new TH1D(*hFitParamHisto[kBkg1Counts]);
1066 hBkg1Fraction->SetNameTitle("hBkg1Fraction", "Bkg1 fraction: Bkg1Counts / IdentifiedCounts;p_{T} (GeV/c);");
1067 hBkg1Fraction->Divide(hBkg1Fraction, hIdentifiedCounts, 1., 1., "B");
1068 /* compute bkg2 fraction */
1069 TH1D *hBkg2Fraction = new TH1D(*hFitParamHisto[kBkg2Counts]);
1070 hBkg2Fraction->SetNameTitle("hBkg2Fraction", "Bkg2 fraction: Bkg2Counts / IdentifiedCounts;p_{T} (GeV/c);");
1071 hBkg2Fraction->Divide(hBkg2Fraction, hIdentifiedCounts, 1., 1., "B");
1072 /* compute bkg3 fraction */
1073 TH1D *hBkg3Fraction = new TH1D(*hFitParamHisto[kBkg3Counts]);
1074 hBkg3Fraction->SetNameTitle("hBkg3Fraction", "Bkg3 fraction: Bkg3Counts / IdentifiedCounts;p_{T} (GeV/c);");
1075 hBkg3Fraction->Divide(hBkg3Fraction, hIdentifiedCounts, 1., 1., "B");
1076 /* compute bkg4 fraction */
1077 TH1D *hBkg4Fraction = new TH1D(*hFitParamHisto[kBkg4Counts]);
1078 hBkg4Fraction->SetNameTitle("hBkg4Fraction", "Bkg4 fraction: Bkg4Counts / IdentifiedCounts;p_{T} (GeV/c);");
1079 hBkg4Fraction->Divide(hBkg4Fraction, hIdentifiedCounts, 1., 1., "B");
1080 /* rework */
1081 TH1D *hSignalFractionReworked = new TH1D(*hFitParamHisto[kSignalCounts]);
1082 hSignalFractionReworked->SetNameTitle("hSignalFractionReworked", "Signal fraction reworked: SignalCounts / IdentifiedCounts;p_{T} (GeV/c);");
1083 hSignalFractionReworked->Multiply(fMismatchBalanceFunction[iipart][icharge]);
1084 TH1D *hBkg1FractionReworked = new TH1D(*hFitParamHisto[kBkg1Counts]);
1085 hBkg1FractionReworked->SetNameTitle("hBkg1FractionReworked", "Bkg1 fraction reworked: SignalCounts / IdentifiedCounts;p_{T} (GeV/c);");
1086 hBkg1FractionReworked->Multiply(fMismatchBalanceFunction[partbkg1[iipart]][icharge]);
1087 TH1D *hBkg2FractionReworked = new TH1D(*hFitParamHisto[kBkg2Counts]);
1088 hBkg2FractionReworked->SetNameTitle("hBkg2FractionReworked", "Bkg1 fraction reworked: SignalCounts / IdentifiedCounts;p_{T} (GeV/c);");
1089 hBkg2FractionReworked->Multiply(fMismatchBalanceFunction[partbkg2[iipart]][icharge]);
1090 TH1D *hIdentifiedCountsReworked = new TH1D(*hSignalFractionReworked);
1091 hIdentifiedCountsReworked->SetNameTitle("hIdentifiedCountsReworked", "Identified counts reworked: SignalCounts / IdentifiedCounts;p_{T} (GeV/c);");
1092 hIdentifiedCountsReworked->Add(hBkg1FractionReworked);
1093 hIdentifiedCountsReworked->Add(hBkg2FractionReworked);
1094
1095 hSignalFractionReworked->Divide(hSignalFractionReworked, hIdentifiedCountsReworked, 1., 1., "B");
1096
1097 /* compute mismatch-correction counts */
1098 TH1D *hMismatchCorrectionCounts = new TH1D(*hTotalMismatchCounts);
1099 hMismatchCorrectionCounts->SetNameTitle("hMismatchCorrectionCounts", "Mismatch-correction counts: TotalMismatchCounts * SignalFraction;p_{T} (GeV/c);counts");
1100 hMismatchCorrectionCounts->Multiply(hSignalFractionReworked);
1101 /* compute mismatch-corrected signal counts */
1102 TH1D *hMismatchCorrectedSignalCounts = new TH1D(*hFitParamHisto[kSignalCounts]);
1103 hMismatchCorrectedSignalCounts->SetNameTitle("hMismatchCorrectedSignalCounts", "Mismatch-corrected signal counts: SignalCounts + MismatchCorrectionCounts;p_{T} (GeV/c);counts");
1104 hMismatchCorrectedSignalCounts->Add(hMismatchCorrectionCounts);
1105
1106
1107 /* accepted tracks histo */
1108 if (hAcceptedTracks) {
1109 TH1D *hAcceptedTracks_py = hAcceptedTracks->ProjectionY();
1110 hAcceptedTracks_py->SetNameTitle("hAcceptedTracks", "Accepted tracks;p_{T} (GeV/c);");
1111 hAcceptedTracks_py->Sumw2();
1112 }
1113
1114 /*** RAW SPECTRA ***/
1115
1116 TDirectory *rawDir = fileout->mkdir("RawSpectra");
1117
1118 /* compute normalized raw yield */
1119 TH1D *hNormalizedRawYield = new TH1D(*hFitParamHisto[kSignalCounts]);
1120 hNormalizedRawYield->SetNameTitle("hNormalizedRawYield", "Raw yield;p_{T} (GeV/c);#frac{1}{N_{ev}} #frac{d^{2}N}{dy dp_{T}}");
1121 TOFpid_normalize(hNormalizedRawYield, nevents);
1122 /* compute normalized mismatch-corrected raw yield */
1123 TH1D *hNormalizedMismatchCorrectedRawYield = new TH1D(*hMismatchCorrectedSignalCounts);
1124 hNormalizedMismatchCorrectedRawYield->SetNameTitle("hNormalizedMismatchCorrectedRawYield", "Mismatch-corrected raw yield;p_{T} (GeV/c);#frac{1}{N_{ev}} #frac{d^{2}N}{dy dp_{T}}");
1125 TOFpid_normalize(hNormalizedMismatchCorrectedRawYield, nevents);
1126
1127 /*** OUTPUT ***/
1128
1129 /* write fir params histos */
1130 fitDir->cd();
1131 for (Int_t iparam = 0; iparam < kNFitParams; iparam++)
1132 hFitParamHisto[iparam]->Write();
1133 /* write post-analysis histos */
1134 postDir->cd();
1135 hOverflowCounts->Write();
1136 hTotalMismatchCounts->Write();
1137 hTotalMismatchFraction->Write();
1138 hIdentifiedCounts->Write();
1139 hSignalFraction->Write();
1140 hSignalFractionReworked->Write();
1141 hBkg1Fraction->Write();
1142 hBkg2Fraction->Write();
1143 hBkg3Fraction->Write();
1144 hBkg4Fraction->Write();
1145 hMismatchCorrectionCounts->Write();
1146 hMismatchCorrectedSignalCounts->Write();
1147 if (hAcceptedTracks) hAcceptedTracks_py->Write();
1148 /* write raw spectra histos */
1149 rawDir->cd();
1150 hNormalizedRawYield->Write();
1151 hNormalizedMismatchCorrectedRawYield->Write();
1152
1153 /* clean up */
1154 delete canvas;
1155 for (Int_t iparam = 0; iparam < kNFitParams; iparam++)
1156 delete hFitParamHisto[iparam];
1157 delete hOverflowCounts;
1158 delete hTotalMismatchCounts;
1159 delete hTotalMismatchFraction;
1160 delete hIdentifiedCounts;
1161 delete hSignalFraction;
1162 delete hBkg1Fraction;
1163 delete hBkg2Fraction;
1164 delete hBkg3Fraction;
1165 delete hBkg4Fraction;
1166 delete hMismatchCorrectionCounts;
1167 delete hMismatchCorrectedSignalCounts;
1168 if (hAcceptedTracks) {
1169 delete hAcceptedTracks_py;
1170 }
1171 delete hNormalizedRawYield;
1172 delete hNormalizedMismatchCorrectedRawYield;
1173
1174 /* close file */
1175 fileout->Close();
1176}
1177
1178//___________________________________________________________________________________
1179
1180Float_t
1181TOFpid_histomin(TH1 *h)
1182{
1183
1184 for (Int_t ibin = 0; ibin < h->GetNbinsX(); ibin++)
1185 if (h->GetBinContent(ibin + 1) > 0.)
1186 return h->GetXaxis()->GetBinCenter(ibin + 1);
1187 return kMaxInt;
1188}
1189
1190//___________________________________________________________________________________
1191
1192Float_t
1193TOFpid_histomax(TH1 *h)
1194{
1195
1196 for (Int_t ibin = h->GetNbinsX(); ibin > 0; ibin--)
1197 if (h->GetBinContent(ibin) > 0.)
1198 return h->GetXaxis()->GetBinCenter(ibin);
1199 return -kMaxInt;
1200}
1201
1202//___________________________________________________________________________________
1203
1204TOFpid_checkneg(TH1 *h)
1205{
1206
1207 for (Int_t ibin = 0; ibin < h->GetNbinsX(); ibin++)
1208 if (h->GetBinContent(ibin + 1) <= 0.) {
1209 h->SetBinContent(ibin + 1, 1.e-300);
1210 // h->SetBinError(ibin + 1, 0.1);
1211 }
1212}
1213
1214//___________________________________________________________________________________
1215
1216Double_t
1217TOFpid_fit(TH1 *hSignal, TH1 *hSigExp, TH1 *hBkgExp1, TH1 *hBkgExp2, TH1 *hBkgExp3, TH1 *hBkgExp4, TH1 *hMismatch, Double_t rangeMin, Double_t rangeMax, Bool_t fitBkg1, Bool_t fitBkg2, Bool_t fitBkg3, Bool_t fitBkg4, Bool_t constrainSignal, Bool_t constrainBkg1, Bool_t constrainBkg2, Bool_t forceGaussianSignal, Double_t *param = NULL, Double_t *param_err = NULL, TCanvas *canvas = NULL)
1218{
1219
1220 /** ROOFIT ***/
1221 gSystem->Load("libRooFit");
1222 using namespace RooFit;
1223 /*** LOAD GAUSSIANTAIL CLASS ***/
b0635849 1224 gSystem->Load("RooFermiCutoff_cxx");
1225 gSystem->Load("RooGaussianTail_cxx");
59e49925 1226
1227 /*** DEFINE FIT RANGE ***/
1228
1229 printf("***** FIT RANGE DEFINITION *****\n");
1230
1231 /* check mismatch histogram to define min/max fit range */
1232 // rangeMin = TMath::Max(rangeMin, TOFpid_histomin(hMismatch));
1233 // rangeMax = TMath::Min(rangeMax, TOFpid_histomax(hMismatch));
1234 /* fix zeroes */
1235 TOFpid_checkneg(hMismatch);
1236
1237 /* define range */
1238 RooRealVar x("x", "n_{#sigma}", 0., rangeMin, rangeMax, "");
1239 printf("FIT RANGE DEFINED: %f -> %f\n", rangeMin, rangeMax);
1240 printf("********************************\n");
1241
1242 /*** DEFINE HISTOGRAM DATA ***/
1243
1244 /* define data to fit and background from input histogram */
1245 RooDataHist hdata("hdata", "hdata", x, hSignal);
1246 RooDataHist hsig("hsig", "hsig", x, hSigExp);
1247 RooDataHist hbkg1("hbkg1", "hbkg1", x, hBkgExp1);
1248 RooDataHist hbkg2("hbkg2", "hbkg2", x, hBkgExp2);
1249 RooDataHist hbkg3("hbkg3", "hbkg3", x, hBkgExp3);
1250 RooDataHist hbkg4("hbkg4", "hbkg4", x, hBkgExp4);
1251 RooDataHist hmismatch("hmismatch", "hmismatch", x, hMismatch);
1252
1253 /*** DEFINE SIGNAL SHAPE ***/
1254
1255 printf("***** SIGNAL SHAPE DEFINITION *****\n");
1256
1257 /* variables */
1258 if (FIX_SIGNAL_MEAN) {
1259 RooRealVar mean("mean", "mean", DEFAULT_SIGNAL_MEAN, "");
1260 printf("FIXED SIGNAL_MEAN = %f\n", mean.getVal());
1261 }
1262 else {
1263 RooRealVar mean("mean", "mean", DEFAULT_SIGNAL_MEAN, MIN_SIGNAL_MEAN, MAX_SIGNAL_MEAN, "");
1264 printf("FREE SIGNAL_MEAN = %f [%f, %f]\n", mean.getVal(), MIN_SIGNAL_MEAN, MAX_SIGNAL_MEAN);
1265 }
1266 if (FIX_SIGNAL_SIGMA) {
1267 RooRealVar sigma("sigma", "sigma", DEFAULT_SIGNAL_SIGMA, "");
1268 printf("FIXED SIGNAL_SIGMA = %f\n", sigma.getVal());
1269 }
1270 else {
1271 RooRealVar sigma("sigma", "sigma", DEFAULT_SIGNAL_SIGMA, MIN_SIGNAL_SIGMA, MAX_SIGNAL_SIGMA, "");
1272 printf("FREE SIGNAL_SIGMA = %f\n", sigma.getVal());
1273 }
1274 if (FIX_SIGNAL_TAIL) {
1275 RooRealVar tail("tail", "tail", DEFAULT_SIGNAL_TAIL, "");
1276 printf("FIXED SIGNAL_TAIL = %f\n", tail.getVal());
1277 }
1278 else {
1279 RooRealVar tail("tail", "tail", DEFAULT_SIGNAL_TAIL, MIN_SIGNAL_TAIL, MAX_SIGNAL_TAIL, "");
1280 printf("FREE SIGNAL_TAIL = %f\n", tail.getVal());
1281 }
1282 RooRealVar gaussianfrac("gaussianfrac", "gaussianfrac", 1., 0., 1., "");
1283 RooRealVar sigalpha("sigalpha", "sigalpha", 0., -10., 0.);
1284
1285 /* shapes */
1286 if (GAUSSIAN_SIGNAL || forceGaussianSignal) {
1287 printf("USING GAUSSIAN_SIGNAL SHAPE\n");
1288 RooGaussian signal("signal", "signal", x, mean, sigma);
1289 }
1290 else if (GAUSSIANTAIL_SIGNAL) {
1291 printf("USING GAUSSIANTAIL_SIGNAL SHAPE\n");
1292 RooGaussianTail signal("signal", "signal", x, mean, sigma, tail);
1293 }
1294 else if (GAUSSIANTAIL2_SIGNAL) {
1295 printf("USING GAUSSIANTAIL2_SIGNAL SHAPE\n");
1296 RooGaussianTail signal("signal", "signal", x, mean, sigma, sigma);
1297 }
1298 else if (GAUSSIANPLUSGAUSSIANTAIL_SIGNAL) {
1299 printf("USING GAUSSIANPLUSGAUSSIANTAIL_SIGNAL SHAPE\n");
1300 RooGaussian gaussian("gaussian", "gaussian", x, mean, sigma);
1301 RooGaussianTail gaussiantail("gaussiantail", "gaussiantail", x, mean, sigma, tail);
1302 RooAddPdf signal("signal", "signal", RooArgList(gaussian, gaussiantail), gaussianfrac, kTRUE);
1303
1304 }
1305 else if (GAUSSIANPLUSEXPONENTIAL_SIGNAL) {
1306 printf("USING GAUSSIANPLUSEXPONENTIAL_SIGNAL SHAPE\n");
1307 RooGaussian gaussian("gaussian", "gaussian", x, mean, sigma);
1308 RooExponential sigexpo("sigexpo", "sigexpo", x, sigalpha);
1309 RooRealVar sigcutoff("sigcutoff", "sigcutoff", 0.);
1310 RooRealVar sigpower("sigpower", "sigpower", 0.1);
1311 RooFermiCutoff sigfermi("sigfermi", "sigfermi", x, sigcutoff, sigpower);
1312 RooProdPdf exposignal("exposignal", "exposignal", sigfermi, sigexpo, -100.);
1313 RooRealVar gaussianfrac("gaussianfrac", "gaussianfrac", 0.9, 0.7, 1., "");
1314 RooAddPdf signal("signal", "signal", RooArgList(gaussian, exposignal), gaussianfrac, kTRUE);
1315 }
1316 else if (EXPECTED_SIGNAL_TEMPLATE) {
1317 printf("SHAPE OF SIGNAL FROM TEMPLATE HISTOGRAM\n");
1318 RooHistPdf signal("signal", "signal", x, hsig);
1319 }
1320 else if (EXPECTED_SIGNAL_FIT) {
1321 /* fitting bkg1 */
1322 TF1 *fGaus = (TF1 *)gROOT->GetFunction("gaus");
1323 hSigExp->Fit(fGaus, "0q");
1324 Double_t sig_mean = fGaus->GetParameter(1);
1325 Double_t sig_sigma = fGaus->GetParameter(2);
1326 mean.setConstant(kFALSE);
1327 mean.setVal(sig_mean);
1328 mean.setRange(sig_mean - 10., sig_mean + 10.);
1329 mean.setConstant(kFALSE);
1330 sigma.setConstant(kFALSE);
1331 sigma.setVal(sig_sigma);
1332 sigma.setRange(sig_sigma * 0.5, sig_sigma * 1.5);
1333 tail.setConstant(kFALSE);
1334 tail.setVal(1.);
1335 tail.setRange(0.5, 5.);
1336 RooGaussianTail signal("signal", "signal", x, mean, sigma, tail);
1337 signal.fitTo(hsig, Range(sig_mean - 5. * sig_sigma, sig_mean + 10. * sig_sigma), SumW2Error(kFALSE), Verbose(kFALSE), PrintEvalErrors(10));
1338#if DISPLAY
1339 TCanvas *cSig_fit = new TCanvas("cSig_fit");
1340 RooPlot *sig_frame = x.frame();
1341 hsig.plotOn(sig_frame);
1342 signal.plotOn(sig_frame, LineColor(kRed));
1343 sig_frame->Draw();
1344 cSig_fit->Update();
1345#endif
1346 printf("SIGNAL PARAMETERS AFTER FIT OF EXPECTED SIGNAL\n");
1347 printf("mean = %f +- %f\n", mean.getVal(), mean.getError());
1348 printf("sigma = %f +- %f\n", sigma.getVal(), sigma.getError());
1349 printf("tail = %f +- %f\n", tail.getVal(), tail.getError());
1350 /* scale parameters if requested */
1351 if (SCALE_SIGNAL_SIGMA != 1.) {
1352 printf("SCALE FITTED SIGNAL SIGMA BY %f\n", SCALE_SIGNAL_SIGMA);
1353 sigma.setVal(sigma.getVal() * SCALE_SIGNAL_SIGMA);
1354 }
1355 if (SCALE_SIGNAL_TAIL != 1.) {
1356 printf("SCALE FITTED SIGNAL TAIL BY %f\n", SCALE_SIGNAL_TAIL);
1357 tail.setVal(tail.getVal() * SCALE_SIGNAL_TAIL);
1358 }
1359 /* fix/release parameters if requested */
1360 if (FIX_SIGNAL_MEAN) {
1361 printf("SETTING FITTED SIGNAL MEAN AS CONSTANTS\n");
1362 mean.setConstant(kTRUE);
1363 }
1364 else {
1365 printf("SETTING FITTED SIGNAL MEAN AS FREE\n");
1366 // mean.setRange(mean.getVal() - 0.25 * TMath::Abs(mean.getVal()), mean.getVal() + 0.25 * TMath::Abs(mean.getVal()));
1367 // mean.setRange(mean.getVal() - 0.5 * TMath::Abs(mean.getVal()), mean.getVal() + 0.5 * TMath::Abs(mean.getVal()));
1368 mean.setRange(mean.getVal() - 0.2, mean.getVal() + 0.2);
1369 }
1370 if (FIX_SIGNAL_SIGMA) {
1371 printf("SETTING FITTED SIGNAL SIGMA AS CONSTANTS\n");
1372 sigma.setConstant(kTRUE);
1373 }
1374 else {
1375 printf("SETTING FITTED SIGNAL SIGMA AS FREE\n");
1376 sigma.setRange(sigma.getVal() * 0.8, sigma.getVal() * 1.2);
1377 // sigma.setRange(sigma.getVal() - 0.5, sigma.getVal() + 0.5);
1378 }
1379 if (FIX_SIGNAL_TAIL) {
1380 printf("SETTING FITTED SIGNAL TAIL AS CONSTANTS\n");
1381 tail.setConstant(kTRUE);
1382 }
1383 else {
1384 printf("SETTING FITTED SIGNAL TAIL AS FREE\n");
1385 // tail.setRange(0.5, 5.0);
1386 tail.setRange(tail.getVal() * 0.8, tail.getVal() * 1.2);
1387 // tail.setRange(tail.getVal() - 0.5, tail.getVal() + 0.5);
1388 }
1389 }
1390 else {
1391 printf("SHAPE OF SIGNAL NOT DEFINE: using GAUSSIAN_SIGNAL\n");
1392 RooGaussian signal("signal", "signal", x, mean, sigma);
1393 }
1394
1395
1396#if 0
1397 if (constrainSignal) {
1398#if 0
1399 /* fix expected signal and constrain parameters if requested */
1400 signal.fitTo(hsig);
1401#if 0
1402 TCanvas *cConstrainSignal = new TCanvas("cConstrainSignal");
1403 RooPlot *xfr = x.frame();
1404 hsig.plotOn(xfr);
1405 signal.plotOn(xfr, LineColor(kRed));
1406 xfr->Draw();
1407 cConstrainSignal->Update();
1408#endif
1409 printf("SIGNAL PARAMETERS AFTER FIT OF EXPECTED SIGNAL\n");
1410 printf("mean = %f +- %f\n", mean.getVal(), mean.getError());
1411 printf("sigma = %f +- %f\n", sigma.getVal(), sigma.getError());
1412 printf("tail = %f +- %f\n", tail.getVal(), tail.getError());
1413 if (constrainSignal) {
1414 mean.setConstant(kTRUE);
1415 sigma.setConstant(kTRUE);
1416 tail.setConstant(kTRUE);
1417 printf("SIGNAL PARAMETERS CONSTRAINED AFTER FIT OF EXPECTED SIGNAL\n");
1418 }
1419#endif
1420 }
1421#endif
1422
1423 if (constrainSignal) {
1424 mean.setConstant(kTRUE);
1425 sigma.setConstant(kTRUE);
1426 tail.setConstant(kTRUE);
1427 // mean.setRange(-0.1, 0.1);
1428 // sigma.setRange(0.95, 1.05);
1429 // tail.setRange(0.95, 1.25);
1430 }
1431
1432 printf("***********************************\n");
1433
1434 /*** DEFINE IDENTIFIED BACKGROUND SHAPES ***/
1435
1436 printf("***** IDENTIFIED BACKGROUND SHAPE DEFINITION *****\n");
1437
1438 /* shapes */
1439if (EXPECTED_BACKGROUND_TEMPLATE) {
1440 printf("USING EXPECTED BACKGROUND TEMPLATE SHAPES FROM HISTOGRAMS\n");
1441 RooHistPdf bkg1("bkg1", "bkg1", x, hbkg1, 0);
1442 RooHistPdf bkg2("bkg2", "bkg2", x, hbkg2, 0);
1443 RooHistPdf bkg3("bkg3", "bkg3", x, hbkg3, 0);
1444 RooHistPdf bkg4("bkg4", "bkg4", x, hbkg4, 0);
1445 }
1446 else if (EXPECTED_BACKGROUND_FIT) {
1447 printf("USING EXPECTED BACKGROUND FITTED SHAPES FROM HISTOGRAMS\n");
1448 /* fitting bkg1 */
1449 TF1 *fGaus = (TF1 *)gROOT->GetFunction("gaus");
1450 hBkgExp1->Fit(fGaus, "0q");
1451 Double_t bkgexp1_mean = fGaus->GetParameter(1);
1452 Double_t bkgexp1_sigma = fGaus->GetParameter(2);
1453 RooRealVar mean_bkg1("mean_bkg1", "mean_bkg1", bkgexp1_mean, bkgexp1_mean - 10., bkgexp1_mean + 10., "");
1454 RooRealVar sigma_bkg1("sigma_bkg1", "sigma_bkg1", bkgexp1_sigma, bkgexp1_sigma * 0.5, bkgexp1_sigma * 1.5, "");
1455 RooRealVar tail_bkg1("tail_bkg1", "tail_bkg1", 1.0, 0.5, 5., "");
1456 RooGaussianTail bkg1("bkg1", "bkg1", x, mean_bkg1, sigma_bkg1, tail_bkg1);
1457 bkg1.fitTo(hbkg1, Range(bkgexp1_mean - 5. * bkgexp1_sigma, bkgexp1_mean + 10. * bkgexp1_sigma), SumW2Error(kFALSE), Verbose(kFALSE), PrintEvalErrors(10));
1458#if DISPLAY
1459 TCanvas *cBkg1_fit = new TCanvas("cBkg1_fit");
1460 RooPlot *bkg1_frame = x.frame();
1461 hbkg1.plotOn(bkg1_frame);
1462 bkg1.plotOn(bkg1_frame, LineColor(kCyan+1));
1463 bkg1_frame->Draw();
1464 cBkg1_fit->Update();
1465#endif
1466 printf("BACKGROUND-1 PARAMETERS AFTER FIT OF EXPECTED BACKGROUND-1\n");
1467 printf("mean_bkg1 = %f +- %f\n", mean_bkg1.getVal(), mean_bkg1.getError());
1468 printf("sigma_bkg1 = %f +- %f\n", sigma_bkg1.getVal(), sigma_bkg1.getError());
1469 printf("tail_bkg1 = %f +- %f\n", tail_bkg1.getVal(), tail_bkg1.getError());
1470 /* scale parameters if requested */
1471 if (SCALE_BACKGROUND_SIGMA != 1.) {
1472 printf("SCALE FITTED BACKGROUND-1 SIGMA BY %f\n", SCALE_BACKGROUND_SIGMA);
1473 sigma_bkg1.setVal(sigma_bkg1.getVal() * SCALE_BACKGROUND_SIGMA);
1474 }
1475 if (SCALE_BACKGROUND_TAIL != 1.) {
1476 printf("SCALE FITTED BACKGROUND-1 TAIL BY %f\n", SCALE_BACKGROUND_TAIL);
1477 tail_bkg1.setVal(tail_bkg1.getVal() * SCALE_BACKGROUND_TAIL);
1478 }
1479 /* fix/release parameters if requested */
1480 if (FIX_BACKGROUND_MEAN) {
1481 printf("SETTING BACKGROUND-1 FITTED MEAN AS CONSTANTS\n");
1482 mean_bkg1.setConstant(kTRUE);
1483 }
1484 else {
1485 printf("SETTING BACKGROUND-1 FITTED MEAN AS FREE\n");
1486 mean_bkg1.setRange(mean_bkg1.getVal() - 0.25 * TMath::Abs(mean_bkg1.getVal()), mean_bkg1.getVal() + 0.25 * TMath::Abs(mean_bkg1.getVal()));
1487 }
1488 if (FIX_BACKGROUND_SIGMA) {
1489 printf("SETTING BACKGROUND-1 FITTED SIGMA AS CONSTANTS\n");
1490 sigma_bkg1.setConstant(kTRUE);
1491 }
1492 else {
1493 printf("SETTING BACKGROUND-1 FITTED SIGMA AS FREE\n");
1494 sigma_bkg1.setRange(sigma_bkg1.getVal() * 0.75, sigma_bkg1.getVal() * 1.25);
1495 }
1496 if (FIX_BACKGROUND_TAIL) {
1497 printf("SETTING BACKGROUND-1 FITTED TAIL AS CONSTANTS\n");
1498 tail_bkg1.setConstant(kTRUE);
1499 }
1500 else {
1501 printf("SETTING BACKGROUND-1 FITTED TAIL AS FREE\n");
1502 tail_bkg1.setRange(tail_bkg1.getVal() * 0.75, tail_bkg1.getVal() * 1.25);
1503 }
1504 /* fitting bkg2 */
1505 TF1 *fGaus = (TF1 *)gROOT->GetFunction("gaus");
1506 hBkgExp2->Fit(fGaus, "0q");
1507 Double_t bkgexp2_mean = fGaus->GetParameter(1);
1508 Double_t bkgexp2_sigma = fGaus->GetParameter(2);
1509 RooRealVar mean_bkg2("mean_bkg2", "mean_bkg2", bkgexp2_mean, bkgexp2_mean - 10., bkgexp2_mean + 10., "");
1510 RooRealVar sigma_bkg2("sigma_bkg2", "sigma_bkg2", bkgexp2_sigma, bkgexp2_sigma * 0.5, bkgexp2_sigma * 1.5, "");
1511 RooRealVar tail_bkg2("tail_bkg2", "tail_bkg2", 1.0, 0.5, 5., "");
1512 RooGaussianTail bkg2("bkg2", "bkg2", x, mean_bkg2, sigma_bkg2, tail_bkg2);
1513 bkg2.fitTo(hbkg2, Range(bkgexp2_mean - 5. * bkgexp2_sigma, bkgexp2_mean + 10. * bkgexp2_sigma), SumW2Error(kFALSE), Verbose(kFALSE), PrintEvalErrors(10));
1514#if DISPLAY
1515 TCanvas *cBkg2_fit = new TCanvas("cBkg2_fit");
1516 RooPlot *bkg2_frame = x.frame();
1517 hbkg2.plotOn(bkg2_frame);
1518 bkg2.plotOn(bkg2_frame, LineColor(kGreen+1));
1519 bkg2_frame->Draw();
1520 cBkg2_fit->Update();
1521#endif
1522 printf("BACKGROUND-2 PARAMETERS AFTER FIT OF EXPECTED BACKGROUND-2\n");
1523 printf("mean_bkg2 = %f +- %f\n", mean_bkg2.getVal(), mean_bkg2.getError());
1524 printf("sigma_bkg2 = %f +- %f\n", sigma_bkg2.getVal(), sigma_bkg2.getError());
1525 printf("tail_bkg2 = %f +- %f\n", tail_bkg2.getVal(), tail_bkg2.getError());
1526 /* scale parameters if requested */
1527 if (SCALE_BACKGROUND_SIGMA != 1.) {
1528 printf("SCALE FITTED BACKGROUND-2 SIGMA BY %f\n", SCALE_BACKGROUND_SIGMA);
1529 sigma_bkg2.setVal(sigma_bkg2.getVal() * SCALE_BACKGROUND_SIGMA);
1530 }
1531 if (SCALE_BACKGROUND_TAIL != 1.) {
1532 printf("SCALE FITTED BACKGROUND-2 TAIL BY %f\n", SCALE_BACKGROUND_TAIL);
1533 tail_bkg2.setVal(tail_bkg2.getVal() * SCALE_BACKGROUND_TAIL);
1534 }
1535 /* fix/release parameters if requested */
1536 if (FIX_BACKGROUND_MEAN) {
1537 printf("SETTING BACKGROUND-2 FITTED MEAN AS CONSTANTS\n");
1538 mean_bkg2.setConstant(kTRUE);
1539 }
1540 else {
1541 printf("SETTING BACKGROUND-2 FITTED MEAN AS FREE\n");
1542 mean_bkg2.setRange(mean_bkg2.getVal() - 0.25 * TMath::Abs(mean_bkg2.getVal()), mean_bkg2.getVal() + 0.25 * TMath::Abs(mean_bkg2.getVal()));
1543 }
1544 if (FIX_BACKGROUND_SIGMA) {
1545 printf("SETTING BACKGROUND-2 FITTED SIGMA AS CONSTANTS\n");
1546 sigma_bkg2.setConstant(kTRUE);
1547 }
1548 else {
1549 printf("SETTING BACKGROUND-2 FITTED SIGMA AS FREE\n");
1550 sigma_bkg2.setRange(sigma_bkg2.getVal() * 0.75, sigma_bkg2.getVal() * 1.25);
1551 }
1552 if (FIX_BACKGROUND_TAIL) {
1553 printf("SETTING BACKGROUND-2 FITTED TAIL AS CONSTANTS\n");
1554 tail_bkg2.setConstant(kTRUE);
1555 }
1556 else {
1557 printf("SETTING BACKGROUND-2 FITTED TAIL AS FREE\n");
1558 tail_bkg2.setRange(tail_bkg2.getVal() * 0.75, tail_bkg2.getVal() * 1.25);
1559 }
1560 /* fitting bkg3 */
1561 TF1 *fGaus = (TF1 *)gROOT->GetFunction("gaus");
1562 hBkgExp3->Fit(fGaus, "0q");
1563 Double_t bkgexp3_mean = fGaus->GetParameter(1);
1564 Double_t bkgexp3_sigma = fGaus->GetParameter(2);
1565 RooRealVar mean_bkg3("mean_bkg3", "mean_bkg3", bkgexp3_mean, bkgexp3_mean - 10., bkgexp3_mean + 10., "");
1566 RooRealVar sigma_bkg3("sigma_bkg3", "sigma_bkg3", bkgexp3_sigma, bkgexp3_sigma * 0.5, bkgexp3_sigma * 1.5, "");
1567 RooRealVar tail_bkg3("tail_bkg3", "tail_bkg3", 1., 0.5, 5., "");
1568 RooGaussianTail bkg3("bkg3", "bkg3", x, mean_bkg3, sigma_bkg3, tail_bkg3);
1569 bkg3.fitTo(hbkg3, Range(bkgexp3_mean - 5. * bkgexp3_sigma, bkgexp3_mean + 10. * bkgexp3_sigma), SumW2Error(kFALSE), Verbose(kFALSE), PrintEvalErrors(10));
1570#if DISPLAY
1571 TCanvas *cBkg3_fit = new TCanvas("cBkg3_fit");
1572 RooPlot *bkg3_frame = x.frame();
1573 hbkg3.plotOn(bkg3_frame);
1574 bkg3.plotOn(bkg3_frame, LineColor(kYellow+1));
1575 bkg3_frame->Draw();
1576 cBkg3_fit->Update();
1577#endif
1578 printf("BACKGROUND-3 PARAMETERS AFTER FIT OF EXPECTED BACKGROUND-3\n");
1579 printf("mean_bkg3 = %f +- %f\n", mean_bkg3.getVal(), mean_bkg3.getError());
1580 printf("sigma_bkg3 = %f +- %f\n", sigma_bkg3.getVal(), sigma_bkg3.getError());
1581 printf("tail_bkg3 = %f +- %f\n", tail_bkg3.getVal(), tail_bkg3.getError());
1582 /* scale parameters if requested */
1583 if (SCALE_BACKGROUND_SIGMA != 1.) {
1584 printf("SCALE FITTED BACKGROUND-3 SIGMA BY %f\n", SCALE_BACKGROUND_SIGMA);
1585 sigma_bkg3.setVal(sigma_bkg3.getVal() * SCALE_BACKGROUND_SIGMA);
1586 }
1587 if (SCALE_BACKGROUND_TAIL != 1.) {
1588 printf("SCALE FITTED BACKGROUND-3 TAIL BY %f\n", SCALE_BACKGROUND_TAIL);
1589 tail_bkg3.setVal(tail_bkg3.getVal() * SCALE_BACKGROUND_TAIL);
1590 }
1591 /* fix/release parameters if requested */
1592 if (FIX_BACKGROUND_MEAN) {
1593 printf("SETTING BACKGROUND-3 FITTED MEAN AS CONSTANTS\n");
1594 mean_bkg3.setConstant(kTRUE);
1595 }
1596 else {
1597 printf("SETTING BACKGROUND-3 FITTED MEAN AS FREE\n");
1598 mean_bkg3.setRange(mean_bkg3.getVal() - 0.25 * TMath::Abs(mean_bkg3.getVal()), mean_bkg3.getVal() + 0.25 * TMath::Abs(mean_bkg3.getVal()));
1599 }
1600 if (FIX_BACKGROUND_SIGMA) {
1601 printf("SETTING BACKGROUND-3 FITTED SIGMA AS CONSTANTS\n");
1602 sigma_bkg3.setConstant(kTRUE);
1603 }
1604 else {
1605 printf("SETTING BACKGROUND-3 FITTED SIGMA AS FREE\n");
1606 sigma_bkg3.setRange(sigma_bkg3.getVal() * 0.75, sigma_bkg3.getVal() * 1.25);
1607 }
1608 if (FIX_BACKGROUND_TAIL) {
1609 printf("SETTING BACKGROUND-3 FITTED TAIL AS CONSTANTS\n");
1610 tail_bkg3.setConstant(kTRUE);
1611 }
1612 else {
1613 printf("SETTING BACKGROUND-3 FITTED TAIL AS FREE\n");
1614 tail_bkg3.setRange(tail_bkg3.getVal() * 0.75, tail_bkg3.getVal() * 1.25);
1615 }
1616 /* fitting bkg4 */
1617 TF1 *fGaus = (TF1 *)gROOT->GetFunction("gaus");
1618 hBkgExp4->Fit(fGaus, "0q");
1619 Double_t bkgexp4_mean = fGaus->GetParameter(1);
1620 Double_t bkgexp4_sigma = fGaus->GetParameter(2);
1621 RooRealVar mean_bkg4("mean_bkg4", "mean_bkg4", bkgexp4_mean, bkgexp4_mean - 10., bkgexp4_mean + 10., "");
1622 RooRealVar sigma_bkg4("sigma_bkg4", "sigma_bkg4", bkgexp4_sigma, bkgexp4_sigma * 0.5, bkgexp4_sigma * 1.5, "");
1623 RooRealVar tail_bkg4("tail_bkg4", "tail_bkg4", 1., 0.5, 5., "");
1624 RooGaussianTail bkg4("bkg4", "bkg4", x, mean_bkg4, sigma_bkg4, tail_bkg4);
1625 bkg4.fitTo(hbkg4, Range(bkgexp4_mean - 5. * bkgexp4_sigma, bkgexp4_mean + 10. * bkgexp4_sigma), SumW2Error(kFALSE), Verbose(kFALSE), PrintEvalErrors(10));
1626#if DISPLAY
1627 TCanvas *cBkg4_fit = new TCanvas("cBkg4_fit");
1628 RooPlot *bkg4_frame = x.frame();
1629 hbkg4.plotOn(bkg4_frame);
1630 bkg4.plotOn(bkg4_frame, LineColor(kYellow+2));
1631 bkg4_frame->Draw();
1632 cBkg4_fit->Update();
1633#endif
1634 printf("BACKGROUND-4 PARAMETERS AFTER FIT OF EXPECTED BACKGROUND-4\n");
1635 printf("mean_bkg4 = %f +- %f\n", mean_bkg4.getVal(), mean_bkg4.getError());
1636 printf("sigma_bkg4 = %f +- %f\n", sigma_bkg4.getVal(), sigma_bkg4.getError());
1637 printf("tail_bkg4 = %f +- %f\n", tail_bkg4.getVal(), tail_bkg4.getError());
1638 /* scale parameters if requested */
1639 if (SCALE_BACKGROUND_SIGMA != 1.) {
1640 printf("SCALE FITTED BACKGROUND-4 SIGMA BY %f\n", SCALE_BACKGROUND_SIGMA);
1641 sigma_bkg4.setVal(sigma_bkg4.getVal() * SCALE_BACKGROUND_SIGMA);
1642 }
1643 if (SCALE_BACKGROUND_TAIL != 1.) {
1644 printf("SCALE FITTED BACKGROUND-4 TAIL BY %f\n", SCALE_BACKGROUND_TAIL);
1645 tail_bkg4.setVal(tail_bkg4.getVal() * SCALE_BACKGROUND_TAIL);
1646 }
1647 /* fix/release parameters if requested */
1648 if (FIX_BACKGROUND_MEAN) {
1649 printf("SETTING BACKGROUND-4 FITTED MEAN AS CONSTANTS\n");
1650 mean_bkg4.setConstant(kTRUE);
1651 }
1652 else {
1653 printf("SETTING BACKGROUND-4 FITTED MEAN AS FREE\n");
1654 mean_bkg4.setRange(mean_bkg4.getVal() - 0.25 * TMath::Abs(mean_bkg4.getVal()), mean_bkg4.getVal() + 0.25 * TMath::Abs(mean_bkg4.getVal()));
1655 }
1656 if (FIX_BACKGROUND_SIGMA) {
1657 printf("SETTING BACKGROUND-4 FITTED SIGMA AS CONSTANTS\n");
1658 sigma_bkg4.setConstant(kTRUE);
1659 }
1660 else {
1661 printf("SETTING BACKGROUND-4 FITTED SIGMA AS FREE\n");
1662 sigma_bkg4.setRange(sigma_bkg4.getVal() * 0.75, sigma_bkg4.getVal() * 1.25);
1663 }
1664 if (FIX_BACKGROUND_TAIL) {
1665 printf("SETTING BACKGROUND-4 FITTED TAIL AS CONSTANTS\n");
1666 tail_bkg4.setConstant(kTRUE);
1667 }
1668 else {
1669 printf("SETTING BACKGROUND-4 FITTED TAIL AS FREE\n");
1670 tail_bkg4.setRange(tail_bkg4.getVal() * 0.75, tail_bkg4.getVal() * 1.25);
1671 }
1672 }
1673 else if (GAUSSIAN_BACKGROUND) {
1674 printf("USING GAUSSIAN BACKGROUND SHAPES (not reccomended)\n");
1675 RooRealVar mean1("mean1", "mean1", hBkgExp1->GetMean(), hBkgExp1->GetMean() * 0.95, hBkgExp1->GetMean() * 1.05, "");
1676 RooRealVar sigma1("sigma1", "sigma1", hBkgExp1->GetRMS(), hBkgExp1->GetRMS() * 0.5, hBkgExp1->GetRMS() * 1.5, "");
1677 RooGaussian bkg1("bkg1", "bkg1", x, mean1, sigma1);
1678
1679 RooRealVar mean2("mean2", "mean2", hBkgExp2->GetMean(), hBkgExp2->GetMean() * 0.95, hBkgExp2->GetMean() * 1.05, "");
1680 RooRealVar sigma2("sigma2", "sigma2", hBkgExp2->GetRMS(), hBkgExp2->GetRMS() * 0.5, hBkgExp2->GetRMS() * 1.5, "");
1681 RooGaussian bkg2("bkg2", "bkg2", x, mean2, sigma2);
1682
1683 RooRealVar mean3("mean3", "mean3", hBkgExp3->GetMean(), hBkgExp3->GetMean() * 0.95, hBkgExp3->GetMean() * 1.05, "");
1684 RooRealVar sigma3("sigma3", "sigma3", hBkgExp3->GetRMS(), hBkgExp3->GetRMS() * 0.5, hBkgExp3->GetRMS() * 1.5, "");
1685 RooGaussian bkg3("bkg3", "bkg3", x, mean3, sigma3);
1686 }
1687 else {
1688 printf("SHAPE OF BACKGROUND NOT DEFINE: using EXPECTED_BACKGROUND\n");
1689 RooHistPdf bkg1("bkg1", "bkg1", x, hbkg1, 0);
1690 RooHistPdf bkg2("bkg2", "bkg2", x, hbkg2, 0);
1691 RooHistPdf bkg3("bkg3", "bkg3", x, hbkg3, 0);
1692 RooHistPdf bkg4("bkg4", "bkg4", x, hbkg4, 0);
1693 }
1694 printf("**************************************************\n");
1695
1696 /*** DEFINE MISMATCH BACKGROUND SHAPE ***/
1697
1698 printf("***** MISMATCH BACKGROUND SHAPE DEFINITION *****\n");
1699
1700 /* variables and generic shapes */
1701 Double_t expectedCutoff = hBkgExp3->GetMean();
1702 Double_t expectedCutoffRMS = hBkgExp3->GetRMS();
1703 // RooRealVar cutoff("cutoff", "cutoff", -30., -50., 0., "");
1704 RooRealVar cutoff("cutoff", "cutoff", expectedCutoff, expectedCutoff - 3. * expectedCutoffRMS, expectedCutoff, "");
1705 // RooRealVar cutoff("cutoff", "cutoff", expectedCutoff, "");
1706 // RooRealVar power("power", "power", 1., 0.5, 5.0, "");
1707 RooRealVar power("power", "power", 1., "");
1708 RooFermiCutoff fermi("fermi", "fermi", x, cutoff, power);
1709 RooRealVar alpha("alpha", "alpha", 0., -10., 0., "");
1710 RooExponential expo("expo", "expo", x, alpha);
1711 RooUniform uniform("uniform", "uniform", x);
1712
1713 /* mismatch shape */
1714 if (EXPECTED_MISMATCH) {
1715 printf("USING EXPECTED MISMATCH SHAPE FROM HISTOGRAMS\n");
1716 RooHistPdf mismatch("mismatch", "mismatch", x, hmismatch, 0);
1717 }
1718 else if (UNIFORM_MISMATCH) {
1719 printf("USING UNIFORM BACKGROUND SHAPE WITH CUTOFF\n");
1720 RooProdPdf mismatch("mismatch", "mismatch", fermi, uniform, -100.);
1721 }
1722 else if (EXPONENTIAL_MISMATCH) {
1723 printf("USING EXPONENTIAL BACKGROUND SHAPE WITH CUTOFF\n");
1724 RooProdPdf mismatch("mismatch", "mismatch", fermi, expo, -100.);
1725 }
1726 else if (DOUBLEEXPONENTIAL_MISMATCH) {
1727 printf("USING DOUBLE-EXPONENTIAL BACKGROUND SHAPE WITH CUTOFF\n");
1728 RooRealVar alpha1("alpha1", "alpha1", 0., -10., 0., "");
1729 RooRealVar alpha2("alpha2", "alpha2", 0., -10., 0., "");
1730 RooRealVar frac("frac", "frac", 0.5, 0., 1., "");
1731 RooGenericPdf doubleexpo("doubleexpo", "TMath::Exp(alpha1 * x) + frac * TMath::Exp(alpha2 * x)", RooArgSet(x, alpha1, alpha2, frac));
1732 RooProdPdf mismatch("mismatch", "mismatch", fermi, doubleexpo, -100.);
1733 }
1734 else if (UNIFORMPLUSEXPONENTIAL_MISMATCH) {
1735 printf("USING UNIFORM-PLUS-EXPONENTIAL BACKGROUND SHAPE WITH CUTOFF\n");
1736 RooRealVar funiform("funiform", "funiform", 100., 0., 100000., "");
1737 RooRealVar fexpo("fexpo", "fexpo", 100., 0., 100000., "");
1738 RooAddPdf uniformplusexpo("uniformplusexpo", "uniformplusexpo", RooArgList(uniform, expo), RooArgList(funiform, fexpo), kFALSE);
1739 RooProdPdf mismatch("mismatch", "mismatch", fermi, uniformplusexpo, -100.);
1740 }
1741 else {
1742 printf("SHAPE OF MISMATCH NOT DEFINE: using EXPECTED_MISMATCH\n");
1743 RooHistPdf mismatch("mismatch", "mismatch", x, hmismatch, 0);
1744 }
1745 printf("************************************************\n");
1746
1747 /*** DEFINE THE MODEL ***/
1748
1749 printf("***** MODEL DEFINITION *****\n");
1750
1751 /* variables */
1752 Double_t integral = hdata.sumEntries();
1753 RooRealVar nsignal("nsignal", "nsignal", 0.3 * integral, 0., integral);
1754 RooRealVar nbkg1("nbkg1", "nbkg1", 0.3 * integral, 0., integral);
1755 RooRealVar nbkg2("nbkg2", "nbkg2", 0.3 * integral, 0., integral);
1756 RooRealVar nbkg3("nbkg3", "nbkg3", 0.3 * integral, 0., integral);
1757 RooRealVar nbkg4("nbkg4", "nbkg4", 0.3 * integral, 0., integral);
1758 RooRealVar nmismatch("nmismatch", "nmismatch", 0.1 * integral, 0., integral);
1759
1760if (!fitBkg1) {
1761 nbkg1.setVal(0.);
1762 nbkg1.setConstant(kTRUE);
1763 }
1764if (!fitBkg2) {
1765 nbkg2.setVal(0.);
1766 nbkg2.setConstant(kTRUE);
1767 }
1768if (!fitBkg3) {
1769 nbkg3.setVal(0.);
1770 nbkg3.setConstant(kTRUE);
1771 }
1772if (!fitBkg4) {
1773 nbkg4.setVal(0.);
1774 nbkg4.setConstant(kTRUE);
1775 }
1776
1777 RooAddPdf model("model", "model p.d.f.", RooArgList(signal, bkg1, bkg2, bkg3, bkg4, mismatch), RooArgList(nsignal, nbkg1, nbkg2, nbkg3, nbkg4, nmismatch));
1778
1779#if 0
1780 /* the model */
1781 if (USE_ELECTRON_BACKGROUND && fitElectrons && fitPions) {
1782 printf("USING ELECTRON BACKGROUND\n");
1783 if (NO_MISMATCH) {
1784 printf("NOT USING MISMATCH BACKGROUND\n");
1785 nmismatch.setVal(0.);
1786 RooAddPdf model("model", "model p.d.f.", RooArgList(signal, bkg1, bkg2, bkg3/*, bkg4*/), RooArgList(nsignal, nbkg1, nbkg2, nbkg3/*, nbkg4*/));
1787 }
1788 else {
1789 printf("USING MISMATCH BACKGROUND\n");
1790 RooAddPdf model("model", "model p.d.f.", RooArgList(signal, bkg1, bkg2, bkg3/*, bkg4*/, mismatch), RooArgList(nsignal, nbkg1, nbkg2, nbkg3/*, nbkg4*/, nmismatch));
1791 }
1792 }
1793 else if (!USE_ELECTRON_BACKGROUND || !fitElectrons) {
1794 printf("NOT USING ELECTRON BACKGROUND\n");
1795 nbkg3.setVal(0.);
1796 nbkg4.setVal(0.);
1797 if (NO_MISMATCH) {
1798 printf("NOT USING MISMATCH BACKGROUND\n");
1799 nmismatch.setVal(0.);
1800 RooAddPdf model("model", "model p.d.f.", RooArgList(signal, bkg1, bkg2), RooArgList(nsignal, nbkg1, nbkg2));
1801 }
1802 else {
1803 printf("USING MISMATCH BACKGROUND\n");
1804 RooAddPdf model("model", "model p.d.f.", RooArgList(signal, bkg1, bkg2, mismatch), RooArgList(nsignal, nbkg1, nbkg2, nmismatch));
1805 }
1806 }
1807 printf("****************************\n");
1808#endif
1809
1810
1811
1812
1813 /*** FIT ***/
1814
1815 printf("***** FIT *****\n");
1816
1817 printf("SIGNAL PARAMETERS BEFORE GLOBAL FIT\n");
1818 printf("mean = %f +- %f\n", mean.getVal(), mean.getError());
1819 printf("sigma = %f +- %f\n", sigma.getVal(), sigma.getError());
1820 printf("tail = %f +- %f\n", tail.getVal(), tail.getError());
1821
1822 /* fit and draw */
1823 if (canvas) canvas->cd();
1824 // model.chi2FitTo(hdata, Extended(kTRUE), Verbose(kFALSE), SumW2Error(kFALSE), Range(-40., 140.), Binned(kTRUE));
1825// model.fitTo(hdata, Extended(kTRUE), SumW2Error(kFALSE), Verbose(kFALSE), PrintEvalErrors(10), Range(-10., 10.));
1826model.fitTo(hdata, Range(rangeMin, rangeMax), Extended(kTRUE), SumW2Error(kFALSE), Verbose(kFALSE), PrintEvalErrors(10));
1827
1828 printf("***************\n");
1829
1830 /*** DRAW ***/
1831#if DISPLAY
1832 RooPlot *xframe = x.frame();
1833hdata.plotOn(xframe, XErrorSize(0), DrawOption("PZ"));
1834 printf("plotting model...\n");
1835model.plotOn(xframe, LineWidth(2)/*, Precision(1.e-4)*/);
1836 printf("plotting signal...\n");
1837model.plotOn(xframe, Components(signal), LineWidth(2), LineColor(kRed)/*, Precision(1.e-4)*/);
1838 printf("plotting bkg1...\n");
1839 model.plotOn(xframe, Components(bkg1), LineWidth(2), LineStyle(kDashed), LineColor(kCyan+1));
1840 printf("plotting bkg2...\n");
1841 model.plotOn(xframe, Components(bkg2), LineWidth(2), LineStyle(kDashed), LineColor(kGreen+1));
1842 if (USE_ELECTRON_BACKGROUND) {
1843 model.plotOn(xframe, Components(bkg3), LineWidth(2), LineStyle(kDashed), LineColor(kYellow+1));
1844 model.plotOn(xframe, Components(bkg4), LineWidth(2), LineStyle(kDashed), LineColor(kYellow+2));
1845 }
1846 if (!NO_MISMATCH)
1847 model.plotOn(xframe, Components(mismatch), LineWidth(2), LineStyle(kDashed), LineColor(kGray+1));
1848 hSignal->SetFillColor(kYellow);
1849 hSignal->SetLineWidth(0);
1850 hSignal->SetFillStyle(0);
1851 hSignal->SetMinimum(0.1);
1852 hSignal->GetXaxis()->SetRangeUser(rangeMin, rangeMax);
1853// hSignal->Draw();
1854xframe->SetMinimum(0.1);
1855 xframe->Draw();
1856#endif
1857 if (canvas) canvas->Update();
1858
1859 /*** COMPUTE CHI2 ***/
1860 Double_t datax, datapoint, datapoint_err, modelpoint;
1861 Double_t chi2 = 0.;
1862 Int_t ndf = 0;
1863 for (Int_t ibin = 0; ibin < hSignal->GetNbinsX(); ibin++) {
1864 datax = hSignal->GetBinCenter(ibin + 1);
1865 if (datax < rangeMin || datax > rangeMax) continue;
1866 datapoint = hSignal->GetBinContent(ibin + 1);
1867 datapoint_err = hSignal->GetBinError(ibin + 1);
1868 if (datapoint_err == 0.) continue;
1869 x.setVal(datax);
1870 modelpoint = model.getVal();
1871 chi2 += (datapoint - modelpoint) * (datapoint - modelpoint) / (datapoint_err * datapoint_err);
1872 ndf++;
1873 }
1874
1875
1876/*** PRINT FIT OUTPUT ***/
1877printf("***** CHI-SQUARE *****\n");
1878 printf("chi-square = %f\n", chi2);
1879 printf("NDF = %d\n", ndf);
1880 printf("chi2/NDF = %f\n", ndf > 0 ? chi2 / ndf : 0.);
1881 printf("***** SIGNAL-SHAPE INFO *****\n");
1882 printf("mean = %f +- %f\n", mean.getVal(), mean.getError());
1883 printf("sigma = %f +- %f\n", sigma.getVal(), sigma.getError());
1884 printf("tail = %f +- %f\n", tail.getVal(), tail.getError());
1885 printf("pure gaus = %f +- %f\n", gaussianfrac.getVal(), gaussianfrac.getError());
1886printf("*****************************\n");
1887printf("***** COUNTS *****\n");
1888printf("total = %f\n", hSignal->Integral(-1, -1));
1889 printf("integral = %f\n", hdata.sumEntries());
1890 printf("nsignal = %f +- %f\n", nsignal.getVal(), nsignal.getError());
1891 printf("nbkg1 = %f +- %f\n", nbkg1.getVal(), nbkg1.getError());
1892 printf("nbkg2 = %f +- %f\n", nbkg2.getVal(), nbkg2.getError());
1893 printf("nbkg3 = %f +- %f\n", nbkg3.getVal(), nbkg3.getError());
1894 printf("nbkg4 = %f +- %f\n", nbkg4.getVal(), nbkg4.getError());
1895 printf("nmismatch = %f +- %f\n", nmismatch.getVal(), nmismatch.getError());
1896 printf("******************\n");
1897
1898 /*** OUTPUT FIT PARAMS ***/
1899
1900 param[kMean] = mean.getVal();
1901 param_err[kMean] = mean.getError();
1902 param[kSigma] = sigma.getVal();
1903 param_err[kSigma] = sigma.getError();
1904 param[kTail] = tail.getVal();
1905 param_err[kTail] = tail.getError();
1906 param[kTotalCounts] = hSignal->GetEntries();
1907 param_err[kTotalCounts] = sqrt(hSignal->GetEntries());
1908 param[kIntegralCounts] = hdata.sumEntries();
1909 param_err[kIntegralCounts] = sqrt(hdata.sumEntries());
1910 param[kSignalCounts] = nsignal.getVal();
1911 param_err[kSignalCounts] = nsignal.getError();
1912 param[kBkg1Counts] = nbkg1.getVal();
1913 param_err[kBkg1Counts] = nbkg1.getError();
1914 param[kBkg2Counts] = nbkg2.getVal();
1915 param_err[kBkg2Counts] = nbkg2.getError();
1916 param[kBkg3Counts] = nbkg3.getVal();
1917 param_err[kBkg3Counts] = nbkg3.getError();
1918 param[kBkg4Counts] = nbkg4.getVal();
1919 param_err[kBkg4Counts] = nbkg4.getError();
1920 param[kMismatchCounts] = nmismatch.getVal();
1921 param_err[kMismatchCounts] = nmismatch.getError();
1922
1923 return;
1924}
1925
1926//___________________________________________________________________________________
1927
1928TOFpid_efficiency(Int_t ipart, Int_t icharge, Int_t icent, Int_t useTPCcut = -1, Double_t minsignalFrac = 0., Int_t nTrials = 1000)
1929{
1930
1931 /* prepare centrality name */
1932 Char_t centName[1024];
1933 if (icent < 0 || icent >= NcentralityBins)
1934 sprintf(centName, "cent0090");
1935 else
1936 sprintf(centName, "cent%02d%02d", centralityBin[icent], centralityBin[icent + 1]);
1937
1938 /* fraction names */
1939 const Char_t *fractionName[AliPID::kSPECIES][AliPID::kSPECIES] = {
1940 "hSignalFraction", "hBkg4Fraction", "hBkg3Fraction", "hBkg1Fraction", "hBkg2Fraction",
1941 "hBkg4Fraction", "hSignalFraction", "hBkg3Fraction", "hBkg1Fraction", "hBkg2Fraction",
1942 "hBkg3Fraction", "hBkg4Fraction", "hSignalFraction", "hBkg1Fraction", "hBkg2Fraction",
1943 "hBkg3Fraction", "hBkg4Fraction", "hBkg1Fraction", "hSignalFraction", "hBkg2Fraction",
1944 "hBkg3Fraction", "hBkg4Fraction", "hBkg1Fraction", "hBkg2Fraction", "hSignalFraction"
1945 };
1946
1947 enum ETPCcut_t {
1948 kCurrentDir,
1949 k11Cut,
1950 k12Cut,
1951 k21Cut,
1952 k22Cut,
1953 k23Cut,
1954 k32Cut,
1955 k33Cut,
1956 kNTPCcuts
1957 };
1958 const Char_t *tpccutdir[8] = {
1959 ".",
1960 "TOFpid_cutOnTPC[-1.0,1.0]",
1961 "TOFpid_cutOnTPC[1.0,2.0]",
1962 "TOFpid_cutOnTPC[-2.0,-1.0]",
1963 "TOFpid_cutOnTPC[-2.0,2.0]",
1964 "TOFpid_cutOnTPC[2.0,3.0]",
1965 "TOFpid_cutOnTPC[-3.0,-2.0]",
1966 "TOFpid_cutOnTPC[-3.0,3.0]"
1967 };
1968
1969 /* get data */
1970 Char_t filename[1024];
1971 TH1D *hAccepted[AliPID::kSPECIES][kNTPCcuts];
1972 TH1D *hIdentified[AliPID::kSPECIES][kNTPCcuts];
1973 TH1D *hEfficiencyIn[AliPID::kSPECIES][kNTPCcuts];
1974 TH1D *hFraction[AliPID::kSPECIES][AliPID::kSPECIES][kNTPCcuts];
1975 for (Int_t itpccut = 0; itpccut < kNTPCcuts; itpccut++) {
1976
1977 /* check whether we use this cut */
1978 if (useTPCcut >= 0 && itpccut != useTPCcut)
1979 continue;
1980
1981 for (Int_t iipart = 0; iipart < AliPID::kSPECIES; iipart++) {
1982 /* skip electrons and muons */
1983 if (iipart == AliPID::kMuon)
1984 continue;
1985
1986 /* open file */
1987 sprintf(filename, "%s/TOFspectrum_%s_%s_%s_%sID.root", tpccutdir[itpccut], centName, AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart));
1988 TFile *filein = TFile::Open(filename);
1989 if (!filein || !filein->IsOpen()) {
1990 printf("cannot open %s\n", filename);
1991 return;
1992 }
1993
1994 /* get accepted tracks */
1995 hAccepted[iipart][itpccut] = (TH1D *)filein->Get("PostAnalysis/hAcceptedTracks");
1996 if (!hAccepted[iipart][itpccut]) {
1997 printf("cannot find PostAnalysis/hAcceptedTracks in %s\n", filename);
1998 return;
1999 }
2000
2001 /* get identified tracks */
2002 hIdentified[iipart][itpccut] = (TH1D *)filein->Get("PostAnalysis/hIdentifiedCounts");
2003 if (!hIdentified[iipart][itpccut]) {
2004 printf("cannot find PostAnalysis/hIdentifiedCounts in %s\n", filename);
2005 return;
2006 }
2007
2008 /* compute efficiency */
2009 hEfficiencyIn[iipart][itpccut] = new TH1D(*hIdentified[iipart][itpccut]);
2010 hEfficiencyIn[iipart][itpccut]->Divide(hEfficiencyIn[iipart][itpccut], hAccepted[iipart][itpccut], 1., 1., "B");
2011
2012 /* get fractions */
2013 for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) {
2014 /* skip muons */
2015 if (iiipart == AliPID::kMuon)
2016 continue;
2017
2018 hFraction[iipart][iiipart][itpccut] = (TH1D *)filein->Get(Form("PostAnalysis/%s", fractionName[iipart][iiipart]));
2019 if (!hFraction[iipart][iiipart][itpccut]) {
2020 printf("cannot find PostAnalysis/%s in %s\n", fractionName[iipart][iiipart], filename);
2021 return;
2022 }
2023 }
2024 }
2025 }
2026
2027 /* prepare output efficiency histos */
2028 TH1D *hEfficiencyOut[AliPID::kSPECIES];
2029 for (Int_t iipart = 0; iipart < AliPID::kSPECIES; iipart++) {
2030 hEfficiencyOut[iipart] = new TH1D(Form("hPIDEff_%d_%s_%s_%sID", icent, AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)), "", NptBins, ptBin);
2031 }
2032
2033
2034 /*** MATRIX ***/
2035
2036 const Int_t nPart = 4;
2037 Int_t partIndex[4] = {AliPID::kElectron, AliPID::kPion, AliPID::kKaon, AliPID::kProton};
2038
2039 TMatrixD Meffin(nPart, 1);
2040 TMatrixD Meffout(nPart, 1);
2041 TMatrixD Mfrac(nPart, nPart);
2042
2043 Double_t eff[4], effe[4], frac[4][4], frace[4][4], geneff, genfrac;
2044 Bool_t gotbestcut[4];
2045
2046 TH1F *hEffTemp[4];
2047 for (Int_t iipart = 0; iipart < nPart; iipart++)
2048 hEffTemp[iipart] = new TH1F(Form("hEffTemp_%d", iipart), "", 20000, -10., 10.);
2049
2050 /* loop over pt-bins */
2051 for (Int_t ibin = 0; ibin < NptBins; ibin++) {
2052
2053 /* reset temp histos */
2054 for (Int_t iipart = 0; iipart < nPart; iipart++)
2055 hEffTemp[iipart]->Reset();
2056
2057 /* get measured data */
2058 for (Int_t iipart = 0; iipart < nPart; iipart++) {
2059
2060 /* select the best set of cuts */
2061 Double_t signalFrac, maxsignalFrac = minsignalFrac;
2062 Int_t bestcut = -1;
2063 gotbestcut[iipart] = kFALSE;
2064 for (Int_t itpccut = 0; itpccut < kNTPCcuts; itpccut++) {
2065
2066 /* check whether we use this cut */
2067 if (useTPCcut >= 0 && itpccut != useTPCcut)
2068 continue;
2069
2070 /* maximize the signal fraction */
2071 signalFrac = hFraction[partIndex[iipart]][partIndex[iipart]][itpccut]->GetBinContent(ibin + 1);
2072 if (signalFrac > maxsignalFrac) {
2073 maxsignalFrac = signalFrac;
2074 bestcut = itpccut;
2075 gotbestcut[iipart] = kTRUE;
2076 }
2077 }
2078 if (bestcut < 0)
2079 continue;
2080
2081 eff[iipart] = hEfficiencyIn[partIndex[iipart]][bestcut]->GetBinContent(ibin + 1);
2082 effe[iipart] = hEfficiencyIn[partIndex[iipart]][bestcut]->GetBinError(ibin + 1);
2083 for (Int_t iiipart = 0; iiipart < nPart; iiipart++) {
2084 frac[iipart][iiipart] = hFraction[partIndex[iipart]][partIndex[iiipart]][bestcut]->GetBinContent(ibin + 1);
2085 frace[iipart][iiipart] = hFraction[partIndex[iipart]][partIndex[iiipart]][bestcut]->GetBinError(ibin + 1);
2086 }
2087 }
2088
2089 /* check best cuts */
2090 Bool_t skip = kFALSE;
2091 for (Int_t iipart = 0; iipart < nPart; iipart++)
2092 if (!gotbestcut[iipart])
2093 skip = kTRUE;
2094 if (skip) continue;
2095
2096 /* loop over trials */
2097 for (Int_t itry = 0; itry < nTrials; itry++) {
2098
2099 /* setup matrix */
2100 for (Int_t iipart = 0; iipart < nPart; iipart++) {
2101 geneff = gRandom->Gaus(eff[iipart], effe[iipart]);
2102 if (geneff < 0.) geneff = 0.;
2103 if (geneff > 1.) geneff = 1.;
2104 Meffin[iipart] = geneff != 0. ? 1. / geneff : 0.;
2105 for (Int_t iiipart = 0; iiipart < nPart; iiipart++) {
2106 genfrac = gRandom->Gaus(frac[iipart][iiipart], frace[iipart][iiipart]);
2107 if (genfrac < 0.) genfrac = 0.;
2108 if (genfrac > 1.) genfrac = 1.;
2109 Mfrac[iipart][iiipart] = genfrac;
2110 }
2111 }
2112
2113 /* invert matrix */
2114 TDecompLU lu(Mfrac);
2115 TMatrixD Minv(nPart, nPart);
2116 if (!lu.Decompose()) continue;
2117 lu.Invert(Minv);
2118 Meffout = Minv * Meffin;
2119
2120 /* fill histos */
2121 for (Int_t iipart = 0; iipart < nPart; iipart++) {
2122 if (Meffout.GetMatrixArray()[iipart] > 0.)
2123 hEffTemp[iipart]->Fill(1. / Meffout.GetMatrixArray()[iipart]);
2124 }
2125
2126 } /* end of loop over trials */
2127
2128
2129 /* get average and RMS */
2130 for (Int_t iipart = 0; iipart < nPart; iipart++) {
2131 hEfficiencyOut[partIndex[iipart]]->SetBinContent(ibin + 1, hEffTemp[iipart]->GetMean());
2132 hEfficiencyOut[partIndex[iipart]]->SetBinError(ibin + 1, hEffTemp[iipart]->GetRMS());
2133 }
2134
2135 } /* end of loop over pt-bins */
2136
2137
2138 hEfficiencyOut[AliPID::kPion]->SetMarkerStyle(20);
2139 hEfficiencyOut[AliPID::kKaon]->SetMarkerStyle(20);
2140 hEfficiencyOut[AliPID::kProton]->SetMarkerStyle(20);
2141 hEfficiencyOut[AliPID::kElectron]->SetMarkerStyle(20);
2142
2143 hEfficiencyOut[AliPID::kPion]->SetMarkerColor(4);
2144 hEfficiencyOut[AliPID::kKaon]->SetMarkerColor(8);
2145 hEfficiencyOut[AliPID::kProton]->SetMarkerColor(2);
2146 hEfficiencyOut[AliPID::kElectron]->SetMarkerColor(1);
2147
2148 hEfficiencyOut[AliPID::kPion]->Draw();
2149 hEfficiencyOut[AliPID::kKaon]->Draw("same");
2150 hEfficiencyOut[AliPID::kProton]->Draw("same");
2151 hEfficiencyOut[AliPID::kElectron]->Draw("same");
2152
2153
2154 /* output */
2155 TFile *fileout = TFile::Open(Form("TOFpid_efficiency_cent%d_%s_%s.root", icent, AliPID::ParticleName(ipart), chargeName[icharge]), "RECREATE");
2156 hEfficiencyOut[AliPID::kPion]->Write();
2157 hEfficiencyOut[AliPID::kKaon]->Write();
2158 hEfficiencyOut[AliPID::kProton]->Write();
2159 hEfficiencyOut[AliPID::kElectron]->Write();
2160 fileout->Close();
2161
2162}
2163
2164TOFpid_normalize(TH1D *h, Double_t nevents = 8.42446600000000000e+06)
2165{
2166
2167 Double_t counts, counts_err;
2168 for (Int_t ibin = 0; ibin < h->GetNbinsX(); ibin++) {
2169 counts = h->GetBinContent(ibin + 1);
2170 counts_err = h->GetBinError(ibin + 1);
2171 counts /= h->GetBinWidth(ibin + 1);
2172 counts_err /= h->GetBinWidth(ibin + 1);
2173 h->SetBinContent(ibin + 1, counts);
2174 h->SetBinError(ibin + 1, counts_err);
2175 }
2176
2177 h->Scale(1. / nevents);
2178
2179}
2180
2181TOFpid_normalizeAndwrite(const Char_t *fileoutname, const Char_t *corrstring = "")
2182{
2183
2184 TFile *fpiplus = TFile::Open("TOFpid_spectrum_pion_plus.root");
2185 TFile *fpiminus = TFile::Open("TOFpid_spectrum_pion_minus.root");
2186 TFile *fkaplus = TFile::Open("TOFpid_spectrum_kaon_plus.root");
2187 TFile *fkaminus = TFile::Open("TOFpid_spectrum_kaon_minus.root");
2188 TFile *fprplus = TFile::Open("TOFpid_spectrum_proton_plus.root");
2189 TFile *fprminus = TFile::Open("TOFpid_spectrum_proton_minus.root");
2190
2191 hpiplus = (TH1F *)fpiplus->Get(Form("hSignal%sCounts", corrstring));
2192 hpiminus = (TH1F *)fpiminus->Get(Form("hSignal%sCounts", corrstring));
2193 hkaplus = (TH1F *)fkaplus->Get(Form("hSignal%sCounts", corrstring));
2194 hkaminus = (TH1F *)fkaminus->Get(Form("hSignal%sCounts", corrstring));
2195 hprplus = (TH1F *)fprplus->Get(Form("hSignal%sCounts", corrstring));
2196 hprminus = (TH1F *)fprminus->Get(Form("hSignal%sCounts", corrstring));
2197
2198 hpiplus->SetName("hpiplus");
2199 hpiminus->SetName("hpiminus");
2200 hkaplus->SetName("hkaplus");
2201 hkaminus->SetName("hkaminus");
2202 hprplus->SetName("hprplus");
2203 hprminus->SetName("hprminus");
2204
2205 TOFpid_normalize(hpiplus);
2206 TOFpid_normalize(hpiminus);
2207 TOFpid_normalize(hkaplus);
2208 TOFpid_normalize(hkaminus);
2209 TOFpid_normalize(hprplus);
2210 TOFpid_normalize(hprminus);
2211
2212 TFile *fileout = TFile::Open(fileoutname, "RECREATE");
2213 hpiplus->Write("hpiplus");
2214 hpiminus->Write("hpiminus");
2215 hkaplus->Write("hkaplus");
2216 hkaminus->Write("hkaminus");
2217 hprplus->Write("hprplus");
2218 hprminus->Write("hprminus");
2219 fileout->Close();
2220
2221}
2222
2223/**************************************************************/
2224
2225TOFpid_rawSpectra(const Char_t *destdir = ".")
2226{
2227
2228 Int_t marker[2] = {20, 21};
2229 Int_t color[AliPID::kSPECIES] = {1, 1, 4, 8, 2};
2230 Char_t *partLatex[AliPID::kSPECIES][2] = {
2231 "", "", "", "", "#pi^{+}", "#pi^{-}", "K^{+}", "K^{-}", "p", "#bar{p}"
2232 };
2233
2234 Float_t deltaRapidity = rapidityMaxCut - rapidityMinCut;
2235
2236 TFile *fileout = TFile::Open("TOF_rawSpectra.root", "RECREATE");
2237 TH1D *hRaw;
2238 Char_t title[1024];
2239 for (Int_t icent = -1; icent < NcentralityBins; icent++)
2240 for (Int_t icharge = 0; icharge < kNCharges; icharge++)
2241 for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++) {
2242 hRaw = TOFpid_rawSpectra(destdir, ipart, icharge, icent);
2243 if (!hRaw) continue;
2244 hRaw->Scale(1. / deltaRapidity);
2245 hRaw->SetMarkerStyle(marker[icharge]);
2246 hRaw->SetMarkerColor(color[ipart]);
2247 hRaw->SetLineColor(1);
2248 hRaw->SetLineWidth(1);
2249 if (icent == -1)
2250 sprintf(title, "%s (MB);p_{T} (GeV/c);#frac{d^{2}N}{dy dp_{T}} (c/GeV);", partLatex[ipart][icharge]);
2251 else
2252 sprintf(title, "%s (%d-%d%%);p_{T} (GeV/c);#frac{d^{2}N}{dy dp_{T}} (c/GeV);", partLatex[ipart][icharge], (Int_t)centralityBin[icent], (Int_t)centralityBin[icent + 1]);
2253 hRaw->SetTitle(title);
2254 if (icent == -1)
2255 hRaw->SetName(Form("hRaw_MB_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]));
2256 else
2257 hRaw->SetName(Form("hRaw_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
2258 fileout->cd();
2259 hRaw->Write();
2260 }
2261
2262 fileout->Close();
2263}
2264
2265TH1D *
2266TOFpid_rawSpectra(const Char_t *dirname, Int_t ipart, Int_t icharge, Int_t icent)
2267{
2268
2269 /* open data */
2270 Char_t outfilename[1024];
2271 if (icent < 0 || icent >= NcentralityBins) {
2272 sprintf(outfilename, "%s/TOFspectrum_cent0090_%s_%s_%sID.root", dirname, AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(ipart));
2273 }
2274 else {
2275 sprintf(outfilename, "%s/TOFspectrum_cent%02d%02d_%s_%s_%sID.root", dirname, centralityBin[icent], centralityBin[icent + 1], AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(ipart));
2276 }
2277 TFile *filein = TFile::Open(outfilename);
2278 if (!filein || !filein->IsOpen()) {
2279 printf("cannot open %s\n", outfilename);
2280 return NULL;
2281 }
2282 /* get data */
2283 TH1D *h = (TH1D *)filein->Get("RawSpectra/hNormalizedRawYield");
2284 if (!h) {
2285 printf("cannot get RawSpectra/hNormalizedRawYield from %s\n", outfilename);
2286 return NULL;
2287 }
2288
2289 return h;
2290
2291}
2292
2293TOFpid_rawSpectra_mismatchCorrected(const Char_t *destdir = ".")
2294{
2295
2296 Int_t marker[2] = {20, 21};
2297 Int_t color[AliPID::kSPECIES] = {1, 1, 4, 8, 2};
2298 Char_t *partLatex[AliPID::kSPECIES][2] = {
2299 "", "", "", "", "#pi^{+}", "#pi^{-}", "K^{+}", "K^{-}", "p", "#bar{p}"
2300 };
2301
2302 Float_t deltaRapidity = rapidityMaxCut - rapidityMinCut;
2303
2304 TFile *fileout = TFile::Open("TOF_rawSpectra_mismatchCorrected.root", "RECREATE");
2305 TH1D *hRaw;
2306 Char_t title[1024];
2307 for (Int_t icent = -1; icent < NcentralityBins; icent++)
2308 for (Int_t icharge = 0; icharge < kNCharges; icharge++)
2309 for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++) {
2310 hRaw = TOFpid_rawSpectra_mismatchCorrected(destdir, ipart, icharge, icent);
2311 if (!hRaw) continue;
2312 hRaw->Scale(1. / deltaRapidity);
2313 hRaw->SetMarkerStyle(marker[icharge]);
2314 hRaw->SetMarkerColor(color[ipart]);
2315 hRaw->SetLineColor(1);
2316 hRaw->SetLineWidth(1);
2317 if (icent == -1)
2318 sprintf(title, "%s (MB);p_{T} (GeV/c);#frac{d^{2}N}{dy dp_{T}} (c/GeV);", partLatex[ipart][icharge]);
2319 else
2320 sprintf(title, "%s (%d-%d%%);p_{T} (GeV/c);#frac{d^{2}N}{dy dp_{T}} (c/GeV);", partLatex[ipart][icharge], (Int_t)centralityBin[icent], (Int_t)centralityBin[icent + 1]);
2321 hRaw->SetTitle(title);
2322 if (icent == -1)
2323 hRaw->SetName(Form("hRaw_MB_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]));
2324 else
2325 hRaw->SetName(Form("hRaw_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
2326 fileout->cd();
2327 hRaw->Write();
2328 }
2329
2330 fileout->Close();
2331}
2332
2333TH1D *
2334TOFpid_rawSpectra_mismatchCorrected(const Char_t *dirname, Int_t ipart, Int_t icharge, Int_t icent)
2335{
2336
2337 /* open data */
2338 Char_t outfilename[1024];
2339 if (icent < 0 || icent >= NcentralityBins) {
2340 sprintf(outfilename, "%s/TOFspectrum_cent0090_%s_%s_%sID.root", dirname, AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(ipart));
2341 }
2342 else {
2343 sprintf(outfilename, "%s/TOFspectrum_cent%02d%02d_%s_%s_%sID.root", dirname, centralityBin[icent], centralityBin[icent + 1], AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(ipart));
2344 }
2345 TFile *filein = TFile::Open(outfilename);
2346 if (!filein || !filein->IsOpen()) {
2347 printf("cannot open %s\n", outfilename);
2348 return NULL;
2349 }
2350 /* get data */
2351 TH1D *h = (TH1D *)filein->Get("RawSpectra/hNormalizedMismatchCorrectedRawYield");
2352 if (!h) {
2353 printf("cannot get RawSpectra/hNormalizedRawYield from %s\n", outfilename);
2354 return NULL;
2355 }
2356
2357 return h;
2358
2359}
2360
2361/**************************************************************/
2362
2363TOFpid_electronCorrection()
2364{
2365
2366 Int_t marker[2] = {20, 21};
2367 Int_t color[AliPID::kSPECIES] = {1, 1, 4, 8, 2};
2368 Char_t *partLatex[AliPID::kSPECIES][2] = {
2369 "", "", "", "", "#pi^{+}", "#pi^{-}", "K^{+}", "K^{-}", "p", "#bar{p}"
2370 };
2371
2372
2373 TFile *fileout = TFile::Open("TOF_electronCorrection.root", "RECREATE");
2374 TH1D *hCorr;
2375 TProfile *pCorrAv = new TProfile("pCorrAv", "", NptBins, ptBin, "s");
2376 Char_t title[1024];
2377 for (Int_t icent = 0; icent < NcentralityBins; icent++)
2378 for (Int_t icharge = 0; icharge < kNCharges; icharge++)
2379 for (Int_t ipart = 2; ipart < 3; ipart++) {
2380 hCorr = TOFpid_electronCorrection(ipart, icharge, icent);
2381 hCorr->SetMarkerStyle(marker[icharge]);
2382 hCorr->SetMarkerColor(color[ipart]);
2383 hCorr->SetLineColor(1);
2384 hCorr->SetLineWidth(1);
2385 sprintf(title, "%s (%d-%d%%);p_{T} (GeV/c);electron-subtracted pion fraction;", partLatex[ipart][icharge], (Int_t)centralityBin[icent], (Int_t)centralityBin[icent + 1]);
2386 hCorr->SetTitle(title);
2387 hCorr->SetName(Form("hElectronCorr_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
2388 fileout->cd();
2389 hCorr->Write();
2390 /* fill for average correction */
2391 for (Int_t ipt = 0; ipt < NptBins; ipt++) {
2392 pCorrAv->Fill(hCorr->GetBinCenter(ipt + 1), hCorr->GetBinContent(ipt + 1));
2393 }
2394 }
2395 hCorr = new TH1D("hElectronCorr_average", ";p_{T} (GeV/c);electron-subtracted pion fraction;", NptBins, ptBin);
2396 for (Int_t ipt = 0; ipt < NptBins; ipt++) {
2397 hCorr->SetBinContent(ipt + 1, pCorrAv->GetBinContent(ipt + 1));
2398 hCorr->SetBinError(ipt + 1, pCorrAv->GetBinError(ipt + 1));
2399 }
2400 hCorr->Write();
2401 fileout->Close();
2402}
2403
2404TH1D *
2405TOFpid_electronCorrection(Int_t ipart, Int_t icharge, Int_t icent)
2406{
2407
2408 Int_t marker[2] = {20, 21};
2409 Int_t color[AliPID::kSPECIES] = {1, 1, 4, 8, 2};
2410 Char_t *partLatex[AliPID::kSPECIES][2] = {
2411 "", "", "", "", "#pi^{+}", "#pi^{-}", "K^{+}", "K^{-}", "p", "#bar{p}"
2412 };
2413
2414 TH1D *hr = TOFpid_systematics_ratio("defaultFit_electronCut", "defaultFit", ipart, icharge, icent, "electronCorrection", "", marker[icharge], color[ipart], kTRUE);
2415 TF1 fOne("fOne", "1.", fitPtMin[ipart], fitPtMax[ipart]);
2416 hr->Add(&fOne, -1.);
2417 hr->Scale(1. / 0.866);
2418 hr->Add(&fOne, 1.);
2419
2420 return hr;
2421}
2422
2423/**************************************************************/
2424
2425SystematicCheck(const Char_t *defaultname = "defaultFit")
2426{
2427
2428 const Int_t ndata = 6;
2429 const Char_t *name[ndata] = {
2430 "signalFit_fixed_scaleSigma_09",
2431 "signalFit_fixed_scaleSigma_11",
2432 "signalFit_fixed_scaleTail_09",
2433 "signalFit_fixed_scaleTail_11",
2434 "signalFit_fixed_scaleSigma_09_scaleTail_11",
2435 "signalFit_fixed_scaleSigma_11_scaleTail_09"
2436 };
2437 for (Int_t idata = 0; idata < ndata; idata++)
2438 SystematicCheck(name[idata], defaultname);
2439
2440}
2441
2442SystematicCheck(const Char_t *checkname, const Char_t *defaultname = "defaultFit")
2443{
2444
2445 gROOT->LoadMacro("HistoUtils.C");
2446
2447 Char_t filename1[1024];
2448 Char_t filename2[1024];
2449
2450 Int_t marker[2] = {20, 25};
2451 Int_t color[AliPID::kSPECIES] = {1, 1, 4, 8, 2};
2452 Char_t *partLatex[AliPID::kSPECIES][2] = {
2453 "", "", "", "", "#pi^{+}", "#pi^{-}", "K^{+}", "K^{-}", "p", "#bar{p}"
2454 };
2455
2456
2457 TFile *fileout = TFile::Open(Form("SystematicCheck_%s.root", checkname), "RECREATE");
2458 Char_t title[1024];
2459 TH1 *hd;
2460 for (Int_t icent = 0; icent < NcentralityBins; icent++)
2461 for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++)
2462 for (Int_t icharge = 0; icharge < kNCharges; icharge++) {
2463
2464 if (icent < 0 || icent >= NcentralityBins) {
2465 sprintf(filename1, "%s/TOFspectrum_cent0090_%s_%s_%sID.root", checkname, AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(ipart));
2466 sprintf(filename2, "%s/TOFspectrum_cent0090_%s_%s_%sID.root", defaultname, AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(ipart));
2467 }
2468 else {
2469 sprintf(filename1, "%s/TOFspectrum_cent%02d%02d_%s_%s_%sID.root", checkname, centralityBin[icent], centralityBin[icent + 1], AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(ipart));
2470 sprintf(filename2, "%s/TOFspectrum_cent%02d%02d_%s_%s_%sID.root", defaultname, centralityBin[icent], centralityBin[icent + 1], AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(ipart));
2471 }
2472
2473 hd = HistoUtils_systematics(filename1, filename2, "RawSpectra/hNormalizedRawYield", NULL);
2474 hd->SetName(Form("hRaw_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
2475 sprintf(title, "%s (%d-%d%%);p_{T} (GeV/c);#frac{d^{2}N}{dy dp_{T}} (c/GeV);", partLatex[ipart][icharge], centralityBin[icent], centralityBin[icent + 1]);
2476 hd->SetTitle(title);
2477 hd->SetLineColor(1);
2478 hd->SetLineWidth(1);
2479 hd->SetMarkerStyle(marker[icharge]);
2480 hd->SetMarkerColor(color[ipart]);
2481 fileout->cd();
2482 hd->Write();
2483 delete hd;
2484 }
2485 fileout->Close();
2486
2487}
2488
2489/**************************************************************/
2490
2491TOFpid_systematics()
2492{
2493
2494 TFile *fileout = TFile::Open("TOFpid_systematics.root", "RECREATE");
2495 TH1D *hSys;
2496 for (Int_t icent = 0; icent < NcentralityBins; icent++)
2497 for (Int_t icharge = 0; icharge < kNCharges; icharge++)
2498 for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++) {
2499 hSys = TOFpid_systematics(ipart, icharge, icent);
2500 fileout->cd();
2501 hSys->Write(Form("hRawSys_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
2502 delete hSys;
2503 }
2504
2505 fileout->Close();
2506}
2507
2508TH1D *
2509TOFpid_systematics(Int_t ipart, Int_t icharge, Int_t icent)
2510{
2511
2512 TH1D *hSignalFit = TOFpid_systematics_signalFit(ipart, icharge, icent);
2513 TH1D *hBkgFit = TOFpid_systematics_bkgFit(ipart, icharge, icent);
2514 TH1D *hFitRange = TOFpid_systematics_fitRange(ipart, icharge, icent);
2515
2516 TH1D *hSys = new TH1D("hSys", "", NptBins, ptBin);
2517 Double_t sigsys, bkgsys, rangesys, totsys = 0.;
2518 for (Int_t ipt = 0; ipt < NptBins; ipt++) {
2519 sigsys = hSignalFit->GetBinError(ipt + 1);
2520 bkgsys = hBkgFit->GetBinError(ipt + 1);
2521 rangesys = hFitRange->GetBinError(ipt + 1);
2522 totsys = TMath::Sqrt(sigsys * sigsys + bkgsys * bkgsys + rangesys * rangesys);
2523 hSys->SetBinContent(ipt + 1, totsys);
2524 }
2525
2526 hSys->Draw();
2527
2528 delete hSignalFit;
2529 delete hBkgFit;
2530 delete hFitRange;
2531
2532 return hSys;
2533}
2534
2535TH1D *
2536TOFpid_systematics_fitRange(Int_t ipart, Int_t icharge, Int_t icent)
2537{
2538
2539 TH1D *hArea = new TH1D("hArea", "", NptBins, ptBin);
2540 hArea->SetMinimum(0.8);
2541 hArea->SetMaximum(1.2);
2542 hArea->Draw();
2543
2544 TH1D *hr = TOFpid_systematics_ratio("default_wideRange", "default", ipart, icharge, icent, "wide", "wide-range fit;p_{T} (GeV/c);ratio wrt. default;", 20, 4);
2545 hr->Draw("same");
2546
2547 TH1D *hSys = new TH1D("hSys", "", NptBins, ptBin);
2548 hSys->SetFillStyle(0);
2549 hSys->SetMarkerSize(0);
2550 for (Int_t ipt = 0; ipt < NptBins; ipt++) {
2551 Double_t val, max = 0.;
2552 if (hr->GetBinContent(ipt + 1) == 0.) continue;
2553 val = TMath::Abs(hr->GetBinContent(ipt + 1) - 1.);
2554 hSys->SetBinContent(ipt + 1, 1.);
2555 hSys->SetBinError(ipt + 1, val);
2556 }
2557 hSys->Draw("same, E2");
2558
2559 delete hr;
2560
2561 return hSys;
2562}
2563
2564TH1D *
2565TOFpid_systematics_signalFit(Int_t ipart, Int_t icharge, Int_t icent)
2566{
2567
2568 const Int_t ndata = 6;
2569 const Char_t *name[ndata] = {
2570 "signalFit_fixed_scaleSigma_09",
2571 "signalFit_fixed_scaleSigma_11",
2572 "signalFit_fixed_scaleTail_09",
2573 "signalFit_fixed_scaleTail_11",
2574 "signalFit_fixed_scaleSigma_09_scaleTail_11",
2575 "signalFit_fixed_scaleSigma_11_scaleTail_09"
2576 };
2577 const Char_t *title[ndata] = {
2578 "-10% #sigma;p_{T} (GeV/c); ratio",
2579 "+10% #sigma;p_{T} (GeV/c); ratio",
2580 "-10% #tau;p_{T} (GeV/c); ratio",
2581 "+10% #tau;p_{T} (GeV/c); ratio",
2582 "-10% #sigma, +10% #tau;p_{T} (GeV/c); ratio",
2583 "+10% #sigma, -10% #tau;p_{T} (GeV/c); ratio"
2584 };
2585 Int_t marker[ndata] = {22, 28, 22, 28, 22, 28};
2586 Int_t color[ndata] = {2, 2, 8, 8, 4, 4};
2587
2588 TH1D *hArea = new TH1D("hArea", "", NptBins, ptBin);
2589 hArea->SetMinimum(0.8);
2590 hArea->SetMaximum(1.2);
2591 hArea->Draw();
2592
2593 TH1D *hr[ndata];
2594 for (Int_t idata = 0; idata < ndata; idata++) {
2595 hr[idata] = TOFpid_systematics_ratio(name[idata], "defaultFit", ipart, icharge, icent, name[idata], title[idata], marker[idata], color[idata]);
2596 hr[idata]->Draw("same");
2597 }
2598
2599 TH1D *hSys = new TH1D("hSys", "", NptBins, ptBin);
2600 hSys->SetFillStyle(0);
2601 hSys->SetMarkerSize(0);
2602 for (Int_t ipt = 0; ipt < NptBins; ipt++) {
2603 Double_t val, max = 0.;
2604 for (Int_t idata = 0; idata < ndata; idata++) {
2605 if (hr[idata]->GetBinContent(ipt + 1) == 0.) continue;
2606 val = TMath::Abs(hr[idata]->GetBinContent(ipt + 1) - 1.);
2607 if (val > 0.9) continue;
2608 if (val > max)
2609 max = val;
2610 }
2611 hSys->SetBinContent(ipt + 1, 1.);
2612 hSys->SetBinError(ipt + 1, max);
2613 }
2614 hSys->Draw("same, E2");
2615
2616 // delete hArea;
2617 // for (Int_t idata = 0; idata < ndata; idata++)
2618 // delete hr[idata];
2619
2620 return hSys;
2621}
2622
2623TH1D *
2624TOFpid_systematics_bkgFit(Int_t ipart, Int_t icharge, Int_t icent)
2625{
2626
2627 const Int_t ndata = 6;
2628 const Char_t *name[ndata] = {
2629 "bkgFit_fixed_scaleSigma_09",
2630 "bkgFit_fixed_scaleSigma_11",
2631 "bkgFit_fixed_scaleTail_09",
2632 "bkgFit_fixed_scaleTail_11",
2633 "bkgFit_fixed_scaleSigma_09_scaleTail_11",
2634 "bkgFit_fixed_scaleSigma_11_scaleTail_09"
2635 };
2636 const Char_t *title[ndata] = {
2637 "-10% #sigma;p_{T} (GeV/c); ratio",
2638 "+10% #sigma;p_{T} (GeV/c); ratio",
2639 "-10% #tau;p_{T} (GeV/c); ratio",
2640 "+10% #tau;p_{T} (GeV/c); ratio",
2641 "-10% #sigma, +10% #tau;p_{T} (GeV/c); ratio",
2642 "+10% #sigma, -10% #tau;p_{T} (GeV/c); ratio"
2643 };
2644 Int_t marker[ndata] = {22, 28, 22, 28, 22, 28};
2645 Int_t color[ndata] = {2, 2, 8, 8, 4, 4};
2646
2647 TH1D *hArea = new TH1D("hArea", "", NptBins, ptBin);
2648 hArea->SetMinimum(0.5);
2649 hArea->SetMaximum(1.5);
2650 hArea->Draw();
2651
2652 TH1D *hr[ndata];
2653 for (Int_t idata = 0; idata < ndata; idata++) {
2654 hr[idata] = TOFpid_systematics_ratio(name[idata], "defaultFit", ipart, icharge, icent, name[idata], title[idata], marker[idata], color[idata]);
2655 hr[idata]->Draw("same");
2656 }
2657
2658 TH1D *hSys = new TH1D("hSys", "", NptBins, ptBin);
2659 hSys->SetFillStyle(0);
2660 hSys->SetMarkerSize(0);
2661 for (Int_t ipt = 0; ipt < NptBins; ipt++) {
2662 Double_t val, max = 0.;
2663 for (Int_t idata = 0; idata < ndata; idata++) {
2664 if (hr[idata]->GetBinContent(ipt + 1) == 0.) continue;
2665 val = TMath::Abs(hr[idata]->GetBinContent(ipt + 1) - 1.);
2666 if (val > 0.9) continue;
2667 if (val > max)
2668 max = val;
2669 }
2670 hSys->SetBinContent(ipt + 1, 1.);
2671 hSys->SetBinError(ipt + 1, max);
2672 }
2673 hSys->Draw("same, E2");
2674
2675 // delete hArea;
2676 // for (Int_t idata = 0; idata < ndata; idata++)
2677 // delete hr[idata];
2678
2679 return hSys;
2680}
2681
2682TH1D *
2683TOFpid_systematics_ratio(const Char_t *dirname1, const Char_t *dirname2, Int_t ipart, Int_t icharge, Int_t icent, const Char_t *name = "rawRatio", const Char_t *title = ";p_{T} (GeV/c);raw yield ratio;", Int_t marker = 20, Int_t color = 2, Bool_t correlated = kFALSE)
2684{
2685
2686 TH1D *hr = new TH1D("hr", "", NptBins, ptBin);
2687
2688 /* open data */
2689 Char_t outfilename1[1024];
2690 Char_t outfilename2[1024];
2691 if (icent < 0 || icent >= NcentralityBins) {
2692 sprintf(outfilename1, "%s/TOFspectrum_cent0090_%s_%s_%sID.root", dirname1, AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(ipart));
2693 sprintf(outfilename2, "%s/TOFspectrum_cent0090_%s_%s_%sID.root", dirname2, AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(ipart));
2694 }
2695 else {
2696 sprintf(outfilename1, "%s/TOFspectrum_cent%02d%02d_%s_%s_%sID.root", dirname1, centralityBin[icent], centralityBin[icent + 1], AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(ipart));
2697 sprintf(outfilename2, "%s/TOFspectrum_cent%02d%02d_%s_%s_%sID.root", dirname2, centralityBin[icent], centralityBin[icent + 1], AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(ipart));
2698 }
2699 TFile *filein1 = TFile::Open(outfilename1);
2700 if (!filein1 || !filein1->IsOpen()) {
2701 printf("cannot open %s\n", outfilename1);
2702 return;
2703 }
2704 TFile *filein2 = TFile::Open(outfilename2);
2705 if (!filein2 || !filein2->IsOpen()) {
2706 printf("cannot open %s\n", outfilename2);
2707 return;
2708 }
2709 /* get data */
2710 TH1D *h1 = (TH1D *)filein1->Get("RawSpectra/hNormalizedRawYield");
2711 if (!h1) {
2712 printf("cannot get RawSpectra/hNormalizedRawYield from %s\n", outfilename1);
2713 return;
2714 }
2715 TH1D *h2 = (TH1D *)filein2->Get("RawSpectra/hNormalizedRawYield");
2716 if (!h2) {
2717 printf("cannot get RawSpectra/hNormalizedRawYield from %s\n", outfilename2);
2718 return;
2719 }
2720 /* ratio */
2721 if (correlated) hr->Divide(h1, h2, 1., 1., "B");
2722 else hr->Divide(h1, h2);
2723 hr->SetNameTitle(name, title);
2724 hr->SetMarkerStyle(marker);
2725 hr->SetMarkerColor(color);
2726
2727 filein1->Close();
2728 filein2->Close();
2729
2730 return hr;
2731
2732}
2733