1 #include "CommonDefs.C"
6 Bool_t GAUSSIAN_SIGNAL = kFALSE;
7 Bool_t GAUSSIANTAIL_SIGNAL = kFALSE;
8 Bool_t GAUSSIANTAIL2_SIGNAL = kFALSE;
9 Bool_t GAUSSIANPLUSGAUSSIANTAIL_SIGNAL = kFALSE;
10 Bool_t GAUSSIANPLUSEXPONENTIAL_SIGNAL = kFALSE;
11 Bool_t EXPECTED_SIGNAL_TAIL = kTRUE;
12 Bool_t EXPECTED_SIGNAL_TEMPLATE = kTRUE;
13 Bool_t EXPECTED_SIGNAL_FIT = kFALSE;
14 /* SIGNAL PARAMETERS */
15 Bool_t FIX_SIGNAL_MEAN = kFALSE;
16 Bool_t FIX_SIGNAL_SIGMA = kFALSE;
17 Bool_t FIX_SIGNAL_TAIL = kFALSE;
18 Float_t SCALE_SIGNAL_SIGMA = 1.;
19 Float_t SCALE_SIGNAL_TAIL = 1.;
21 Char_t *SIGNAL_PARAM_FILE = NULL;//"signalParamFile.root";
22 Float_t DEFAULT_SIGNAL_MEAN = 0.;
23 Float_t MIN_SIGNAL_MEAN = -0.2;
24 Float_t MAX_SIGNAL_MEAN = 0.2.;
25 Float_t DEFAULT_SIGNAL_SIGMA = 1.;
26 Float_t MIN_SIGNAL_SIGMA = 0.8;
27 Float_t MAX_SIGNAL_SIGMA = 1.2;
28 Float_t DEFAULT_SIGNAL_TAIL = 1.;
29 Float_t MIN_SIGNAL_TAIL = 0.5;
30 Float_t MAX_SIGNAL_TAIL = 1.5;
32 Bool_t EXPECTED_BACKGROUND_TAIL = kTRUE;
33 Bool_t EXPECTED_BACKGROUND_TEMPLATE = kTRUE;
34 Bool_t EXPECTED_BACKGROUND_FIT = kFALSE;
35 Bool_t GAUSSIAN_BACKGROUND = kFALSE;
36 Bool_t USE_ELECTRON_BACKGROUND = kTRUE;
37 /* BACKGROUND PARAMETERS */
38 Bool_t FIX_BACKGROUND_MEAN = kTRUE;
39 Bool_t FIX_BACKGROUND_SIGMA = kTRUE;
40 Bool_t FIX_BACKGROUND_TAIL = kTRUE;
41 Float_t SCALE_BACKGROUND_SIGMA = 1.;
42 Float_t SCALE_BACKGROUND_TAIL = 1.;
44 Bool_t NO_MISMATCH = kFALSE;
45 Bool_t EXPECTED_MISMATCH = kTRUE;
46 Bool_t EXPONENTIAL_MISMATCH = kFALSE;
47 Bool_t UNIFORM_MISMATCH = kFALSE;
48 Bool_t DOUBLEEXPONENTIAL_MISMATCH = kFALSE;
49 Bool_t UNIFORMPLUSEXPONENTIAL_MISMATCH = kFALSE;
50 /* AUTOMATIC CONSTRAINS */
51 Float_t FIT_ELECTRONS_PT_MIN = 0.;
52 Float_t FIT_ELECTRONS_PT_MAX = 1.0;
53 Float_t FIT_MUONS_PT_MIN = 0.7;
54 Float_t FIT_MUONS_PT_MAX = 0.;
55 Float_t FIT_PIONS_PT_MIN = 0.;
56 Float_t FIT_PIONS_PT_MAX = 5.;
58 Float_t CONSTRAINSIGNAL_LIMIT = 0.;
59 Float_t CONSTRAINTAIL_LIMIT = 0.;
76 /* fit output params name */
77 const Char_t *fitParamName[kNFitParams] = {
91 /* fit output params title */
92 const 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"
106 /**************************************************************/
107 /*** GENERATION OF TEMPLATE HISTOS ****************************/
108 /**************************************************************/
110 const Int_t NmismatchTrials = 1;
111 const Int_t NexpectedTrials = 1;
113 /**************************************************************/
114 /*** HISTOS AND BINNING ***************************************/
115 /**************************************************************/
117 /**************************************************************/
118 /**************************************************************/
126 /**************************************************************/
127 /**************************************************************/
128 const Int_t NtofsigmaBins = 1750;
129 Double_t tofsigmaBin[NtofsigmaBins + 1];
130 Double_t tofsigmaMin = -100., tofsigmaMax = 250., tofsigmaStep = (tofsigmaMax - tofsigmaMin) / NtofsigmaBins;
131 /**************************************************************/
132 const Int_t NtofsignalBins = 2000;
133 Double_t tofsignalBin[NtofsignalBins + 1];
134 Double_t tofsignalMin = -24400., tofsignalMax = 24400., tofsignalStep = (tofsignalMax - tofsignalMin) / NtofsignalBins;
135 /**************************************************************/
136 const Int_t NtofresoBins = 500;
137 Double_t tofresoBin[NtofsignalBins + 1];
138 Double_t tofresoMin = 0., tofresoMax = 500., tofresoStep = (tofresoMax - tofresoMin) / NtofresoBins;
139 /**************************************************************/
140 /**************************************************************/
141 const Int_t NtpcsigmaBins = 10;
142 Double_t tpcsigmaBin[NtpcsigmaBins + 1];
143 Double_t tpcsigmaMin = -5., tpcsigmaMax = 5., tpcsigmaStep = (tpcsigmaMax - tpcsigmaMin) / NtpcsigmaBins;
144 /**************************************************************/
145 Int_t NparamsBins[kNHistoParams] = {NcentralityBins, NtofsigmaBins, NptBins, NtpcsigmaBins};
146 Double_t *paramBin[kNHistoParams] = {centralityBin, tofsigmaBin, ptBin, tpcsigmaBin};
147 /**************************************************************/
149 /**************************************************************/
150 /**************************************************************/
151 /**************************************************************/
153 /**************************************************************/
155 AliTOFGeometry tofGeo;
156 Float_t c = TMath::C() * 1.e2 / 1.e12; /* cm/ps */
157 Float_t c_1 = 1. / c;
160 GenerateRandomHit(TH1F *hT0Fill, Double_t t0fill, Int_t index)
164 Float_t length, timeexp, pos[3];
166 /* compute length and expected time */
167 tofGeo.GetVolumeIndices(index, det);
168 tofGeo.GetPosPar(det, pos);
170 for (Int_t i = 0; i < 3; i++) length += pos[i] * pos[i];
171 length = TMath::Sqrt(length);
172 timeexp = length * c_1;
174 Double_t hittime = hT0Fill->GetRandom() - t0fill + timeexp;
179 /**************************************************************/
181 TOFpid(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)
184 printf("****************************************\n");
185 printf("RUNNING TOF PID:\n");
187 printf("RAPIDITY-CUT: %s\n", AliPID::ParticleName(ipart));
189 printf("NO RAPIDITY-CUT\n");
190 printf("CHARGE: %s\n", chargeName[icharge]);
191 printf("PARTICLE-ID: %s\n", AliPID::ParticleName(iipart));
193 printf("-> ELECTRON CUT REQUESTED\n");
195 printf(" -> CUT-ON-TPC [%3.1f,%3.1f]\n", tpcsignalMin, tpcsignalMax);
196 printf("****************************************\n");
198 /* include path for ACLic */
199 gSystem->AddIncludePath("-I$ALICE_ROOT/include");
200 gSystem->AddIncludePath("-I$ALICE_ROOT/TOF");
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");
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);
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);
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;
233 /* open enabled flag map */
234 TH1F *hEnabledFlag = NULL;
235 if (enabledChannelsFileName) {
236 TFile *enabledfile = TFile::Open(enabledChannelsFileName);
237 hEnabledFlag = (TH1F *)enabledfile->Get("hEnabledFlag");
240 /**************************************************************/
241 /*** HISTOS ***************************************************/
242 /**************************************************************/
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;
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);
258 TH3I *hTOFreso = new TH3I(Form("hTOFreso_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)), "", NcentralityBins, centralityBin, NptBins, ptBin, NtofresoBins, tofresoBin);
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);
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);
275 /**************************************************************/
276 /**************************************************************/
277 /**************************************************************/
279 /* TOF PID response */
280 AliTOFPIDResponse tofResponse;
281 tofResponse.SetTimeResolution(tofReso);
282 /* TPC PID response */
283 AliTPCPIDResponse *tpcResponse = AliAnalysisTrack::GetTPCResponse();
285 /* start stopwatch */
290 printf("***** RUNNING for %s %s *****\n", chargeName[icharge], AliPID::ParticleName(ipart));
292 printf("***** NO RAPIDITY CUT REQUESTED *****\n");
294 printf("***** CUT-ON-TPC requested %3.1f-%3.1f *****\n", tpcsignalMin, tpcsignalMax);
297 /* loop over events */
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;
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;
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");
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");
333 Double_t TZEROA_mean;
334 Double_t TZEROA_sigma;
335 Double_t TZEROC_mean;
336 Double_t TZEROC_sigma;
338 Double_t TOF_TZEROA_mean;
339 Double_t TOF_TZEROC_mean;
340 Double_t TOF_TZEROTOF_mean;
342 Double_t TZEROA_reso;
345 Double_t TZEROC_reso;
348 Double_t TZEROAND_reso;
351 Double_t TZEROTOF_reso;
354 Double_t TZEROMEAN_weight;
356 Double_t TZEROBEST_reso;
358 Double_t param[kNHistoParams];
359 for (Int_t iev = startEv; iev < treein->GetEntries() && iev < evMax; iev++) {
361 treein->GetEvent(iev);
362 if (iev % 10000 == 0) printf("iev = %d\n", iev);
364 if (!analysisEvent->AcceptEvent(acceptEventType)) continue;
366 /*** ACCEPTED EVENT ***/
368 /* apply time-zero TOF correction */
369 // analysisEvent->ApplyTimeZeroTOFCorrection();
372 cent = analysisEvent->GetCentralityPercentile(centralityEstimator);
375 Double_t vertexz = analysisEvent->GetVertexZ();
378 hAcceptedEvents->Fill(cent);
380 /* TZERO corrections */
382 for (icent = 0; icent < NcentralityBins; icent++)
383 if (cent < centralityBin[icent + 1])
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.;
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.;
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.;
399 /* simulate inefficient TZERO */
400 Bool_t resetTZERO = kFALSE;
401 if (forcetimezeroineff > 0.)
402 if (gRandom->Uniform() < forcetimezeroineff)
405 /* loop over tracks */
406 for (Int_t itrk = 0; itrk < analysisTrackArray->GetEntries(); itrk++) {
408 analysisTrack = (AliAnalysisTrack *)analysisTrackArray->At(itrk);
409 if (!analysisTrack) continue;
410 /* check accepted track */
411 if (!analysisTrack->AcceptTrack()) continue;
414 if (analysisTrack->GetY(AliPID::ParticleMass(ipart)) - rapidityShift > rapidityMaxCut ||
415 analysisTrack->GetY(AliPID::ParticleMass(ipart)) - rapidityShift < rapidityMinCut) continue;
418 charge = analysisTrack->GetSign() > 0. ? kPositive : kNegative;
419 if (charge != icharge) continue;
421 /*** ACCEPTED TRACK ***/
424 p = analysisTrack->GetP();
425 pt = analysisTrack->GetPt();
427 /* compute track mt */
428 mt = TMath::Sqrt(pt * pt + AliPID::ParticleMass(ipart) * AliPID::ParticleMass(ipart));
431 dedx = analysisTrack->GetTPCdEdx();
432 dedxN = analysisTrack->GetTPCdEdxN();
433 ptpc = analysisTrack->GetTPCmomentum();
436 bethe = tpcResponse->GetExpectedSignal(ptpc, iipart);
437 /* fix electron signal */
438 if (iipart == AliPID::kElectron)
440 deltadedx = dedx - bethe;
441 dedx_sigma = 0.07 * bethe;
442 tpcsignal = deltadedx / dedx_sigma;
444 /* electronCut requested, remove electrons */
447 bethe = tpcResponse->GetExpectedSignal(ptpc, AliPID::kElectron);
448 /* fix electron signal */
450 deltadedx = dedx - bethe;
451 dedx_sigma = 0.07 * bethe;
452 tpcsignal = deltadedx / dedx_sigma;
453 if (TMath::Abs(tpcsignal) < 1.5) continue;
456 /* cut on TPC signal if requested */
457 if (cutOnTPC && (tpcsignal < tpcsignalMin || tpcsignal > tpcsignalMax))
461 hAcceptedTracks->Fill(cent, pt);
463 /* set TOF pid flag */
464 hastofpid = analysisTrack->HasTOFPID(hEnabledFlag);
470 /*** ACCEPTED TRACK WITH TOF PID ***/
472 /* apply expected time correction */
473 // analysisTrack->ApplyTOFExpectedTimeCorrection();
476 index = analysisTrack->GetTOFIndex();
477 time = analysisTrack->GetTOFTime() - TOF_mean;
478 time_sigma = tofReso;
480 TZEROTOF = analysisEvent->GetTimeZeroTOF(analysisTrack->GetP());
481 TZEROTOF_reso = analysisEvent->GetTimeZeroTOFSigma(analysisTrack->GetP());
482 hasTZEROTOF = TZEROTOF_reso < 150.;
484 TZEROTOF += TOF_TZEROTOF_mean - TOF_mean;
485 TZEROTOF_reso *= TZEROTOF_resoScaleFactor;
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;
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;
498 TZEROAND = (TZEROA + TZEROC) * 0.5;
499 // TZEROAND_reso = TZEROAND_resolution[icent];
500 hasTZEROAND = hasTZEROA && hasTZEROC;
502 TZEROMEAN = TZEROTOF / TZEROTOF_reso / TZEROTOF_reso;
503 TZEROMEAN_weight = 1. / TZEROTOF_reso / TZEROTOF_reso;
505 // printf("TZEROAND\n");
506 TZEROMEAN += TZEROAND / TZEROAND_reso / TZEROAND_reso;
507 TZEROMEAN_weight = 1. / TZEROAND_reso / TZEROAND_reso;
509 else if (hasTZEROA) {
510 // printf("TZEROA\n");
511 TZEROMEAN += TZEROA / TZEROA_reso / TZEROA_reso;
512 TZEROMEAN_weight = 1. / TZEROA_reso / TZEROA_reso;
514 else if (hasTZEROC) {
515 // printf("TZEROC\n");
516 TZEROMEAN += TZEROC / TZEROC_reso / TZEROC_reso;
517 TZEROMEAN_weight = 1. / TZEROC_reso / TZEROC_reso;
519 timezero = TZEROMEAN / TZEROMEAN_weight;
520 timezero_sigma = TMath::Sqrt(1. / TZEROMEAN_weight);
522 TZEROBEST = TZEROTOF;
523 TZEROBEST_reso = TZEROTOF_reso;
524 if (hasTZEROAND && TZEROAND_reso < TZEROBEST_reso) {
525 TZEROBEST = TZEROAND;
526 TZEROBEST_reso = TZEROAND_reso;
528 else if (hasTZEROA && TZEROA_reso < TZEROBEST_reso) {
530 TZEROBEST_reso = TZEROA_reso;
532 if (hasTZEROC && TZEROC_reso < TZEROBEST_reso) {
534 TZEROBEST_reso = TZEROC_reso;
536 timezero = TZEROBEST;
537 timezero_sigma = TZEROBEST_reso;
540 // timezero = 0.;//TZEROTOF;
541 // timezero_sigma = 203.854691;//TZEROTOF_reso;
543 // if (timezero == 0.)
544 // printf("%f %f\n", timezero, timezero_sigma);
546 timezero_sigma *= scaletimezerosigma;
550 timezero_sigma = timezero_spread;
554 tof = time - timezero;
555 tof_sigma = TMath::Sqrt(time_sigma * time_sigma + timezero_sigma * timezero_sigma);
557 /* TOF expected time */
558 texp = analysisTrack->GetTOFExpTime(iipart);
559 texp_sigma = analysisTrack->GetTOFExpTimeSigma(iipart) * scaletexpreso[iipart];
563 deltat_sigma = TMath::Sqrt(tof_sigma * tof_sigma + texp_sigma * texp_sigma);
564 tofsignal = deltat / deltat_sigma;
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);
574 /*** EXPECTED MISMATCH ***/
576 /* loop to generated random hits */
577 for (Int_t irnd = 0; irnd < NmismatchTrials; irnd++) {
579 /* generate ramdom tof values */
580 tof_rnd = GenerateRandomHit(hT0Fill, t0fill, index);
583 deltat = tof_rnd - texp;
584 tofsignal = deltat / deltat_sigma;
587 hTOFmismatch->Fill(cent, pt, tofsignal);
588 hTOFmismatch_delta->Fill(cent, p, deltat);
590 } /* end of loop over generated random hits */
592 /*** EXPECTED SIGNALS ***/
594 /* loop over other particles */
595 for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) {
597 /* generate expected tof value */
598 tof_th = analysisTrack->GetTOFExpTime(iiipart);
599 texp_sigma = analysisTrack->GetTOFExpTimeSigma(iiipart) * scaletexpreso[iiipart];
601 /* loop over many trials */
602 for (Int_t irnd = 0; irnd < NexpectedTrials; irnd++) {
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);
611 /* deltat and sigma */
612 deltat = tof_th - texp + signal_smear + timezero_smear + texp_smear;
613 tofsignal = deltat / deltat_sigma;
616 hTOFexpected[iiipart]->Fill(cent, pt, tofsignal);
617 hTOFexpected_delta[iiipart]->Fill(cent, p, deltat);
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 */
629 TString outputstring = "TOFpid";
631 outputstring += "_noRapidityCut";
633 outputstring += "_electronCut";
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();
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();
655 //___________________________________________________________________________________
657 TOFspectra_defaultFit(const Char_t *filename)
660 Bool_t EXPECTED_SIGNAL_TEMPLATE = kTRUE;
661 Bool_t EXPECTED_SIGNAL_FIT = kFALSE;
662 Bool_t EXPECTED_BACKGROUND_TEMPLATE = kFALSE;
663 Bool_t EXPECTED_BACKGROUND_FIT = kTRUE;
667 //___________________________________________________________________________________
669 TOFspectra_defaultFit_fitElectrons(const Char_t *filename, Float_t electronLimit = 5.)
673 TOFspectra(filename, electronLimit);
676 //___________________________________________________________________________________
678 TOFspectra_signalFit(Bool_t fixParams = kTRUE, Float_t scaleSigma = 1., Float_t scaleTail = 1.)
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;
691 //___________________________________________________________________________________
693 TOFspectra_bkgFit(Bool_t fixParams = kTRUE, Float_t scaleSigma = 1., Float_t scaleTail = 1.)
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;
704 TOFspectra(filename);
708 //___________________________________________________________________________________
711 TOFspectra(const Char_t *filename, Float_t electronLimit = 0.)
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);
721 //___________________________________________________________________________________
726 Double_t fitPtMin[AliPID::kSPECIES] = {0.5, 0.5, 0.3, 0.4, 0.5};
727 Double_t fitPtMax[AliPID::kSPECIES] = {3.0, 3.0, 3.0, 3.0, 5.0};
729 Double_t fitSigmaMin[AliPID::kSPECIES] = {tofsigmaMin, tofsigmaMin, -25., -75., -65.};
730 Double_t fitSigmaMax[AliPID::kSPECIES] = {tofsigmaMax, tofsigmaMax, 225., 200., 100.};
732 Double_t fitSigmaMin[AliPID::kSPECIES] = {-24400., -24400., -24400., -24400., -24400.};
733 Double_t fitSigmaMax[AliPID::kSPECIES] = {24400., 24400., 24400., 24400., 24400.};
737 TOFspectrum(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)
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");
749 /* mismatch balance functions */
750 TF1 *fMismatchBalanceFunction[5][2];
751 Double_t mismatchBalanceParam0[5][2] = {
758 Double_t mismatchBalanceParam1[5][2] = {
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]);
774 TFile *filein = TFile::Open(filename);
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)));
783 if (icent < 0 || icent >= NcentralityBins)
784 nevents = hAcceptedEvents->GetEntries();
786 nevents = hAcceptedEvents->Integral(icent + 1, icent + 1);
787 printf("N_EVENTS : %d\n", nevents);
788 printf("****************************************\n");
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)));
797 TH3I *hTOFpid = (TH3I *)filein->Get(Form("hTOFpid_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)));
799 TH3I *hTOFpid = (TH3I *)filein->Get(Form("hTOFpid_delta_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)));
802 printf("cannot find %s\n", Form("hTOFpid_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)));
806 TH3I *hTOFmismatch = (TH3I *)filein->Get(Form("hTOFmismatch_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)));
808 TH3I *hTOFmismatch = (TH3I *)filein->Get(Form("hTOFmismatch_delta_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)));
811 printf("cannot find %s\n", Form("hTOFpid_%s_%s_%sID", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart)));
814 TH3I *hTOFexpected[AliPID::kSPECIES];
815 for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) {
817 hTOFexpected[iiipart] = (TH3I *)filein->Get(Form("hTOFexpected_%s_%s_%sID_%sBKG", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart), AliPID::ParticleName(iiipart)));
819 hTOFexpected[iiipart] = (TH3I *)filein->Get(Form("hTOFexpected_delta_%s_%s_%sID_%sBKG", AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart), AliPID::ParticleName(iiipart)));
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)));
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);
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);
846 Bool_t requestedRange = kFALSE;
847 Bool_t fitElectrons = kTRUE;
848 Bool_t fitMuons = kTRUE;
849 Bool_t fitPions = kTRUE;
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;
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 */
862 if ((ptMin + 0.001) < FIT_MUONS_PT_MIN || (ptMax - 0.001) > FIT_MUONS_PT_MAX)
864 /* check pion-fit is allowed */
866 if ((ptMin + 0.001) < FIT_PIONS_PT_MIN || (ptMax - 0.001) > FIT_PIONS_PT_MAX)
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);
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));
881 sprintf(outfilename, "TOFspectrum_cent%02d%02d_%s_%s_%sID.root", centralityBin[icent], centralityBin[icent + 1], AliPID::ParticleName(ipart), chargeName[icharge], AliPID::ParticleName(iipart));
883 TFile *fileout = TFile::Open(outfilename, "RECREATE");
884 TDirectory *fitDir = fileout->mkdir("FitParams");
886 TCanvas *canvas = new TCanvas("canvas");
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);
893 /* loop over ptBins */
894 for (Int_t ipt = 0; ipt < NptBins; ipt++) {
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]);
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 */
906 if ((ptBin[ipt] + 0.001) < FIT_MUONS_PT_MIN || (ptBin[ipt + 1] - 0.001) > FIT_MUONS_PT_MAX)
908 /* check pion-fit is allowed */
910 if ((ptBin[ipt] + 0.001) < FIT_PIONS_PT_MIN || (ptBin[ipt + 1] - 0.001) > FIT_PIONS_PT_MAX)
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);
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");
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]];
941 /* check histos if requested */
942 if (checkHistoFlag) {
943 TCanvas *cCheckHisto = new TCanvas("cCheckHisto");
944 cCheckHisto->Divide(2, 3);
956 hMismatch_py->Draw();
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;
967 /* check whether we can fit electrons */
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;
975 /* check whether we can fit muons */
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;
983 /* check whether we can fit pions */
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;
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);
999 /* if fitting pions add electrons and to the signal */
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]);
1004 /* else add to bkg1 (pions) */
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]);
1010 /* check requested pt-range */
1016 canvas->Write(Form("fitDisplay_ptBin_%3.2f_%3.2f", ptBin[ipt], ptBin[ipt + 1]));
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]);
1026 delete hMismatch_py;
1027 for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) {
1028 delete hSignalExp_py[iiipart];
1029 // delete hSignalExpTail_py[iiipart];
1033 /* check requested pt-range */
1037 /*** POST-ANALYSIS ***/
1039 TDirectory *postDir = fileout->mkdir("PostAnalysis");
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");
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);
1095 hSignalFractionReworked->Divide(hSignalFractionReworked, hIdentifiedCountsReworked, 1., 1., "B");
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);
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();
1114 /*** RAW SPECTRA ***/
1116 TDirectory *rawDir = fileout->mkdir("RawSpectra");
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);
1129 /* write fir params histos */
1131 for (Int_t iparam = 0; iparam < kNFitParams; iparam++)
1132 hFitParamHisto[iparam]->Write();
1133 /* write post-analysis histos */
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 */
1150 hNormalizedRawYield->Write();
1151 hNormalizedMismatchCorrectedRawYield->Write();
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;
1171 delete hNormalizedRawYield;
1172 delete hNormalizedMismatchCorrectedRawYield;
1178 //___________________________________________________________________________________
1181 TOFpid_histomin(TH1 *h)
1184 for (Int_t ibin = 0; ibin < h->GetNbinsX(); ibin++)
1185 if (h->GetBinContent(ibin + 1) > 0.)
1186 return h->GetXaxis()->GetBinCenter(ibin + 1);
1190 //___________________________________________________________________________________
1193 TOFpid_histomax(TH1 *h)
1196 for (Int_t ibin = h->GetNbinsX(); ibin > 0; ibin--)
1197 if (h->GetBinContent(ibin) > 0.)
1198 return h->GetXaxis()->GetBinCenter(ibin);
1202 //___________________________________________________________________________________
1204 TOFpid_checkneg(TH1 *h)
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);
1214 //___________________________________________________________________________________
1217 TOFpid_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)
1221 gSystem->Load("libRooFit");
1222 using namespace RooFit;
1223 /*** LOAD GAUSSIANTAIL CLASS ***/
1224 gSystem->Load("RooFermiCutoff_cxx");
1225 gSystem->Load("RooGaussianTail_cxx");
1227 /*** DEFINE FIT RANGE ***/
1229 printf("***** FIT RANGE DEFINITION *****\n");
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));
1235 TOFpid_checkneg(hMismatch);
1238 RooRealVar x("x", "n_{#sigma}", 0., rangeMin, rangeMax, "");
1239 printf("FIT RANGE DEFINED: %f -> %f\n", rangeMin, rangeMax);
1240 printf("********************************\n");
1242 /*** DEFINE HISTOGRAM DATA ***/
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);
1253 /*** DEFINE SIGNAL SHAPE ***/
1255 printf("***** SIGNAL SHAPE DEFINITION *****\n");
1258 if (FIX_SIGNAL_MEAN) {
1259 RooRealVar mean("mean", "mean", DEFAULT_SIGNAL_MEAN, "");
1260 printf("FIXED SIGNAL_MEAN = %f\n", mean.getVal());
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);
1266 if (FIX_SIGNAL_SIGMA) {
1267 RooRealVar sigma("sigma", "sigma", DEFAULT_SIGNAL_SIGMA, "");
1268 printf("FIXED SIGNAL_SIGMA = %f\n", sigma.getVal());
1271 RooRealVar sigma("sigma", "sigma", DEFAULT_SIGNAL_SIGMA, MIN_SIGNAL_SIGMA, MAX_SIGNAL_SIGMA, "");
1272 printf("FREE SIGNAL_SIGMA = %f\n", sigma.getVal());
1274 if (FIX_SIGNAL_TAIL) {
1275 RooRealVar tail("tail", "tail", DEFAULT_SIGNAL_TAIL, "");
1276 printf("FIXED SIGNAL_TAIL = %f\n", tail.getVal());
1279 RooRealVar tail("tail", "tail", DEFAULT_SIGNAL_TAIL, MIN_SIGNAL_TAIL, MAX_SIGNAL_TAIL, "");
1280 printf("FREE SIGNAL_TAIL = %f\n", tail.getVal());
1282 RooRealVar gaussianfrac("gaussianfrac", "gaussianfrac", 1., 0., 1., "");
1283 RooRealVar sigalpha("sigalpha", "sigalpha", 0., -10., 0.);
1286 if (GAUSSIAN_SIGNAL || forceGaussianSignal) {
1287 printf("USING GAUSSIAN_SIGNAL SHAPE\n");
1288 RooGaussian signal("signal", "signal", x, mean, sigma);
1290 else if (GAUSSIANTAIL_SIGNAL) {
1291 printf("USING GAUSSIANTAIL_SIGNAL SHAPE\n");
1292 RooGaussianTail signal("signal", "signal", x, mean, sigma, tail);
1294 else if (GAUSSIANTAIL2_SIGNAL) {
1295 printf("USING GAUSSIANTAIL2_SIGNAL SHAPE\n");
1296 RooGaussianTail signal("signal", "signal", x, mean, sigma, sigma);
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);
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);
1316 else if (EXPECTED_SIGNAL_TEMPLATE) {
1317 printf("SHAPE OF SIGNAL FROM TEMPLATE HISTOGRAM\n");
1318 RooHistPdf signal("signal", "signal", x, hsig);
1320 else if (EXPECTED_SIGNAL_FIT) {
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);
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));
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));
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);
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);
1359 /* fix/release parameters if requested */
1360 if (FIX_SIGNAL_MEAN) {
1361 printf("SETTING FITTED SIGNAL MEAN AS CONSTANTS\n");
1362 mean.setConstant(kTRUE);
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);
1370 if (FIX_SIGNAL_SIGMA) {
1371 printf("SETTING FITTED SIGNAL SIGMA AS CONSTANTS\n");
1372 sigma.setConstant(kTRUE);
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);
1379 if (FIX_SIGNAL_TAIL) {
1380 printf("SETTING FITTED SIGNAL TAIL AS CONSTANTS\n");
1381 tail.setConstant(kTRUE);
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);
1391 printf("SHAPE OF SIGNAL NOT DEFINE: using GAUSSIAN_SIGNAL\n");
1392 RooGaussian signal("signal", "signal", x, mean, sigma);
1397 if (constrainSignal) {
1399 /* fix expected signal and constrain parameters if requested */
1402 TCanvas *cConstrainSignal = new TCanvas("cConstrainSignal");
1403 RooPlot *xfr = x.frame();
1405 signal.plotOn(xfr, LineColor(kRed));
1407 cConstrainSignal->Update();
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");
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);
1432 printf("***********************************\n");
1434 /*** DEFINE IDENTIFIED BACKGROUND SHAPES ***/
1436 printf("***** IDENTIFIED BACKGROUND SHAPE DEFINITION *****\n");
1439 if (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);
1446 else if (EXPECTED_BACKGROUND_FIT) {
1447 printf("USING EXPECTED BACKGROUND FITTED SHAPES FROM HISTOGRAMS\n");
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));
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));
1464 cBkg1_fit->Update();
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);
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);
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);
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()));
1488 if (FIX_BACKGROUND_SIGMA) {
1489 printf("SETTING BACKGROUND-1 FITTED SIGMA AS CONSTANTS\n");
1490 sigma_bkg1.setConstant(kTRUE);
1493 printf("SETTING BACKGROUND-1 FITTED SIGMA AS FREE\n");
1494 sigma_bkg1.setRange(sigma_bkg1.getVal() * 0.75, sigma_bkg1.getVal() * 1.25);
1496 if (FIX_BACKGROUND_TAIL) {
1497 printf("SETTING BACKGROUND-1 FITTED TAIL AS CONSTANTS\n");
1498 tail_bkg1.setConstant(kTRUE);
1501 printf("SETTING BACKGROUND-1 FITTED TAIL AS FREE\n");
1502 tail_bkg1.setRange(tail_bkg1.getVal() * 0.75, tail_bkg1.getVal() * 1.25);
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));
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));
1520 cBkg2_fit->Update();
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);
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);
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);
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()));
1544 if (FIX_BACKGROUND_SIGMA) {
1545 printf("SETTING BACKGROUND-2 FITTED SIGMA AS CONSTANTS\n");
1546 sigma_bkg2.setConstant(kTRUE);
1549 printf("SETTING BACKGROUND-2 FITTED SIGMA AS FREE\n");
1550 sigma_bkg2.setRange(sigma_bkg2.getVal() * 0.75, sigma_bkg2.getVal() * 1.25);
1552 if (FIX_BACKGROUND_TAIL) {
1553 printf("SETTING BACKGROUND-2 FITTED TAIL AS CONSTANTS\n");
1554 tail_bkg2.setConstant(kTRUE);
1557 printf("SETTING BACKGROUND-2 FITTED TAIL AS FREE\n");
1558 tail_bkg2.setRange(tail_bkg2.getVal() * 0.75, tail_bkg2.getVal() * 1.25);
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));
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));
1576 cBkg3_fit->Update();
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);
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);
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);
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()));
1600 if (FIX_BACKGROUND_SIGMA) {
1601 printf("SETTING BACKGROUND-3 FITTED SIGMA AS CONSTANTS\n");
1602 sigma_bkg3.setConstant(kTRUE);
1605 printf("SETTING BACKGROUND-3 FITTED SIGMA AS FREE\n");
1606 sigma_bkg3.setRange(sigma_bkg3.getVal() * 0.75, sigma_bkg3.getVal() * 1.25);
1608 if (FIX_BACKGROUND_TAIL) {
1609 printf("SETTING BACKGROUND-3 FITTED TAIL AS CONSTANTS\n");
1610 tail_bkg3.setConstant(kTRUE);
1613 printf("SETTING BACKGROUND-3 FITTED TAIL AS FREE\n");
1614 tail_bkg3.setRange(tail_bkg3.getVal() * 0.75, tail_bkg3.getVal() * 1.25);
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));
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));
1632 cBkg4_fit->Update();
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);
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);
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);
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()));
1656 if (FIX_BACKGROUND_SIGMA) {
1657 printf("SETTING BACKGROUND-4 FITTED SIGMA AS CONSTANTS\n");
1658 sigma_bkg4.setConstant(kTRUE);
1661 printf("SETTING BACKGROUND-4 FITTED SIGMA AS FREE\n");
1662 sigma_bkg4.setRange(sigma_bkg4.getVal() * 0.75, sigma_bkg4.getVal() * 1.25);
1664 if (FIX_BACKGROUND_TAIL) {
1665 printf("SETTING BACKGROUND-4 FITTED TAIL AS CONSTANTS\n");
1666 tail_bkg4.setConstant(kTRUE);
1669 printf("SETTING BACKGROUND-4 FITTED TAIL AS FREE\n");
1670 tail_bkg4.setRange(tail_bkg4.getVal() * 0.75, tail_bkg4.getVal() * 1.25);
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);
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);
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);
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);
1694 printf("**************************************************\n");
1696 /*** DEFINE MISMATCH BACKGROUND SHAPE ***/
1698 printf("***** MISMATCH BACKGROUND SHAPE DEFINITION *****\n");
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);
1713 /* mismatch shape */
1714 if (EXPECTED_MISMATCH) {
1715 printf("USING EXPECTED MISMATCH SHAPE FROM HISTOGRAMS\n");
1716 RooHistPdf mismatch("mismatch", "mismatch", x, hmismatch, 0);
1718 else if (UNIFORM_MISMATCH) {
1719 printf("USING UNIFORM BACKGROUND SHAPE WITH CUTOFF\n");
1720 RooProdPdf mismatch("mismatch", "mismatch", fermi, uniform, -100.);
1722 else if (EXPONENTIAL_MISMATCH) {
1723 printf("USING EXPONENTIAL BACKGROUND SHAPE WITH CUTOFF\n");
1724 RooProdPdf mismatch("mismatch", "mismatch", fermi, expo, -100.);
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.);
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.);
1742 printf("SHAPE OF MISMATCH NOT DEFINE: using EXPECTED_MISMATCH\n");
1743 RooHistPdf mismatch("mismatch", "mismatch", x, hmismatch, 0);
1745 printf("************************************************\n");
1747 /*** DEFINE THE MODEL ***/
1749 printf("***** MODEL DEFINITION *****\n");
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);
1762 nbkg1.setConstant(kTRUE);
1766 nbkg2.setConstant(kTRUE);
1770 nbkg3.setConstant(kTRUE);
1774 nbkg4.setConstant(kTRUE);
1777 RooAddPdf model("model", "model p.d.f.", RooArgList(signal, bkg1, bkg2, bkg3, bkg4, mismatch), RooArgList(nsignal, nbkg1, nbkg2, nbkg3, nbkg4, nmismatch));
1781 if (USE_ELECTRON_BACKGROUND && fitElectrons && fitPions) {
1782 printf("USING ELECTRON BACKGROUND\n");
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*/));
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));
1793 else if (!USE_ELECTRON_BACKGROUND || !fitElectrons) {
1794 printf("NOT USING ELECTRON BACKGROUND\n");
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));
1803 printf("USING MISMATCH BACKGROUND\n");
1804 RooAddPdf model("model", "model p.d.f.", RooArgList(signal, bkg1, bkg2, mismatch), RooArgList(nsignal, nbkg1, nbkg2, nmismatch));
1807 printf("****************************\n");
1815 printf("***** FIT *****\n");
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());
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.));
1826 model.fitTo(hdata, Range(rangeMin, rangeMax), Extended(kTRUE), SumW2Error(kFALSE), Verbose(kFALSE), PrintEvalErrors(10));
1828 printf("***************\n");
1832 RooPlot *xframe = x.frame();
1833 hdata.plotOn(xframe, XErrorSize(0), DrawOption("PZ"));
1834 printf("plotting model...\n");
1835 model.plotOn(xframe, LineWidth(2)/*, Precision(1.e-4)*/);
1836 printf("plotting signal...\n");
1837 model.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));
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);
1854 xframe->SetMinimum(0.1);
1857 if (canvas) canvas->Update();
1859 /*** COMPUTE CHI2 ***/
1860 Double_t datax, datapoint, datapoint_err, modelpoint;
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;
1870 modelpoint = model.getVal();
1871 chi2 += (datapoint - modelpoint) * (datapoint - modelpoint) / (datapoint_err * datapoint_err);
1876 /*** PRINT FIT OUTPUT ***/
1877 printf("***** 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());
1886 printf("*****************************\n");
1887 printf("***** COUNTS *****\n");
1888 printf("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");
1898 /*** OUTPUT FIT PARAMS ***/
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();
1926 //___________________________________________________________________________________
1928 TOFpid_efficiency(Int_t ipart, Int_t icharge, Int_t icent, Int_t useTPCcut = -1, Double_t minsignalFrac = 0., Int_t nTrials = 1000)
1931 /* prepare centrality name */
1932 Char_t centName[1024];
1933 if (icent < 0 || icent >= NcentralityBins)
1934 sprintf(centName, "cent0090");
1936 sprintf(centName, "cent%02d%02d", centralityBin[icent], centralityBin[icent + 1]);
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"
1958 const Char_t *tpccutdir[8] = {
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]"
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++) {
1977 /* check whether we use this cut */
1978 if (useTPCcut >= 0 && itpccut != useTPCcut)
1981 for (Int_t iipart = 0; iipart < AliPID::kSPECIES; iipart++) {
1982 /* skip electrons and muons */
1983 if (iipart == AliPID::kMuon)
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);
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);
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);
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");
2013 for (Int_t iiipart = 0; iiipart < AliPID::kSPECIES; iiipart++) {
2015 if (iiipart == AliPID::kMuon)
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);
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);
2036 const Int_t nPart = 4;
2037 Int_t partIndex[4] = {AliPID::kElectron, AliPID::kPion, AliPID::kKaon, AliPID::kProton};
2039 TMatrixD Meffin(nPart, 1);
2040 TMatrixD Meffout(nPart, 1);
2041 TMatrixD Mfrac(nPart, nPart);
2043 Double_t eff[4], effe[4], frac[4][4], frace[4][4], geneff, genfrac;
2044 Bool_t gotbestcut[4];
2047 for (Int_t iipart = 0; iipart < nPart; iipart++)
2048 hEffTemp[iipart] = new TH1F(Form("hEffTemp_%d", iipart), "", 20000, -10., 10.);
2050 /* loop over pt-bins */
2051 for (Int_t ibin = 0; ibin < NptBins; ibin++) {
2053 /* reset temp histos */
2054 for (Int_t iipart = 0; iipart < nPart; iipart++)
2055 hEffTemp[iipart]->Reset();
2057 /* get measured data */
2058 for (Int_t iipart = 0; iipart < nPart; iipart++) {
2060 /* select the best set of cuts */
2061 Double_t signalFrac, maxsignalFrac = minsignalFrac;
2063 gotbestcut[iipart] = kFALSE;
2064 for (Int_t itpccut = 0; itpccut < kNTPCcuts; itpccut++) {
2066 /* check whether we use this cut */
2067 if (useTPCcut >= 0 && itpccut != useTPCcut)
2070 /* maximize the signal fraction */
2071 signalFrac = hFraction[partIndex[iipart]][partIndex[iipart]][itpccut]->GetBinContent(ibin + 1);
2072 if (signalFrac > maxsignalFrac) {
2073 maxsignalFrac = signalFrac;
2075 gotbestcut[iipart] = kTRUE;
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);
2089 /* check best cuts */
2090 Bool_t skip = kFALSE;
2091 for (Int_t iipart = 0; iipart < nPart; iipart++)
2092 if (!gotbestcut[iipart])
2096 /* loop over trials */
2097 for (Int_t itry = 0; itry < nTrials; itry++) {
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;
2114 TDecompLU lu(Mfrac);
2115 TMatrixD Minv(nPart, nPart);
2116 if (!lu.Decompose()) continue;
2118 Meffout = Minv * Meffin;
2121 for (Int_t iipart = 0; iipart < nPart; iipart++) {
2122 if (Meffout.GetMatrixArray()[iipart] > 0.)
2123 hEffTemp[iipart]->Fill(1. / Meffout.GetMatrixArray()[iipart]);
2126 } /* end of loop over trials */
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());
2135 } /* end of loop over pt-bins */
2138 hEfficiencyOut[AliPID::kPion]->SetMarkerStyle(20);
2139 hEfficiencyOut[AliPID::kKaon]->SetMarkerStyle(20);
2140 hEfficiencyOut[AliPID::kProton]->SetMarkerStyle(20);
2141 hEfficiencyOut[AliPID::kElectron]->SetMarkerStyle(20);
2143 hEfficiencyOut[AliPID::kPion]->SetMarkerColor(4);
2144 hEfficiencyOut[AliPID::kKaon]->SetMarkerColor(8);
2145 hEfficiencyOut[AliPID::kProton]->SetMarkerColor(2);
2146 hEfficiencyOut[AliPID::kElectron]->SetMarkerColor(1);
2148 hEfficiencyOut[AliPID::kPion]->Draw();
2149 hEfficiencyOut[AliPID::kKaon]->Draw("same");
2150 hEfficiencyOut[AliPID::kProton]->Draw("same");
2151 hEfficiencyOut[AliPID::kElectron]->Draw("same");
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();
2164 TOFpid_normalize(TH1D *h, Double_t nevents = 8.42446600000000000e+06)
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);
2177 h->Scale(1. / nevents);
2181 TOFpid_normalizeAndwrite(const Char_t *fileoutname, const Char_t *corrstring = "")
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");
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));
2198 hpiplus->SetName("hpiplus");
2199 hpiminus->SetName("hpiminus");
2200 hkaplus->SetName("hkaplus");
2201 hkaminus->SetName("hkaminus");
2202 hprplus->SetName("hprplus");
2203 hprminus->SetName("hprminus");
2205 TOFpid_normalize(hpiplus);
2206 TOFpid_normalize(hpiminus);
2207 TOFpid_normalize(hkaplus);
2208 TOFpid_normalize(hkaminus);
2209 TOFpid_normalize(hprplus);
2210 TOFpid_normalize(hprminus);
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");
2223 /**************************************************************/
2225 TOFpid_rawSpectra(const Char_t *destdir = ".")
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}"
2234 Float_t deltaRapidity = rapidityMaxCut - rapidityMinCut;
2236 TFile *fileout = TFile::Open("TOF_rawSpectra.root", "RECREATE");
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);
2250 sprintf(title, "%s (MB);p_{T} (GeV/c);#frac{d^{2}N}{dy dp_{T}} (c/GeV);", partLatex[ipart][icharge]);
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);
2255 hRaw->SetName(Form("hRaw_MB_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]));
2257 hRaw->SetName(Form("hRaw_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
2266 TOFpid_rawSpectra(const Char_t *dirname, Int_t ipart, Int_t icharge, Int_t icent)
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));
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));
2277 TFile *filein = TFile::Open(outfilename);
2278 if (!filein || !filein->IsOpen()) {
2279 printf("cannot open %s\n", outfilename);
2283 TH1D *h = (TH1D *)filein->Get("RawSpectra/hNormalizedRawYield");
2285 printf("cannot get RawSpectra/hNormalizedRawYield from %s\n", outfilename);
2293 TOFpid_rawSpectra_mismatchCorrected(const Char_t *destdir = ".")
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}"
2302 Float_t deltaRapidity = rapidityMaxCut - rapidityMinCut;
2304 TFile *fileout = TFile::Open("TOF_rawSpectra_mismatchCorrected.root", "RECREATE");
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);
2318 sprintf(title, "%s (MB);p_{T} (GeV/c);#frac{d^{2}N}{dy dp_{T}} (c/GeV);", partLatex[ipart][icharge]);
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);
2323 hRaw->SetName(Form("hRaw_MB_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]));
2325 hRaw->SetName(Form("hRaw_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
2334 TOFpid_rawSpectra_mismatchCorrected(const Char_t *dirname, Int_t ipart, Int_t icharge, Int_t icent)
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));
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));
2345 TFile *filein = TFile::Open(outfilename);
2346 if (!filein || !filein->IsOpen()) {
2347 printf("cannot open %s\n", outfilename);
2351 TH1D *h = (TH1D *)filein->Get("RawSpectra/hNormalizedMismatchCorrectedRawYield");
2353 printf("cannot get RawSpectra/hNormalizedRawYield from %s\n", outfilename);
2361 /**************************************************************/
2363 TOFpid_electronCorrection()
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}"
2373 TFile *fileout = TFile::Open("TOF_electronCorrection.root", "RECREATE");
2375 TProfile *pCorrAv = new TProfile("pCorrAv", "", NptBins, ptBin, "s");
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]));
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));
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));
2405 TOFpid_electronCorrection(Int_t ipart, Int_t icharge, Int_t icent)
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}"
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);
2423 /**************************************************************/
2425 SystematicCheck(const Char_t *defaultname = "defaultFit")
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"
2437 for (Int_t idata = 0; idata < ndata; idata++)
2438 SystematicCheck(name[idata], defaultname);
2442 SystematicCheck(const Char_t *checkname, const Char_t *defaultname = "defaultFit")
2445 gROOT->LoadMacro("HistoUtils.C");
2447 Char_t filename1[1024];
2448 Char_t filename2[1024];
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}"
2457 TFile *fileout = TFile::Open(Form("SystematicCheck_%s.root", checkname), "RECREATE");
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++) {
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));
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));
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]);
2489 /**************************************************************/
2491 TOFpid_systematics()
2494 TFile *fileout = TFile::Open("TOFpid_systematics.root", "RECREATE");
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);
2501 hSys->Write(Form("hRawSys_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
2509 TOFpid_systematics(Int_t ipart, Int_t icharge, Int_t icent)
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);
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);
2536 TOFpid_systematics_fitRange(Int_t ipart, Int_t icharge, Int_t icent)
2539 TH1D *hArea = new TH1D("hArea", "", NptBins, ptBin);
2540 hArea->SetMinimum(0.8);
2541 hArea->SetMaximum(1.2);
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);
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);
2557 hSys->Draw("same, E2");
2565 TOFpid_systematics_signalFit(Int_t ipart, Int_t icharge, Int_t icent)
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"
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"
2585 Int_t marker[ndata] = {22, 28, 22, 28, 22, 28};
2586 Int_t color[ndata] = {2, 2, 8, 8, 4, 4};
2588 TH1D *hArea = new TH1D("hArea", "", NptBins, ptBin);
2589 hArea->SetMinimum(0.8);
2590 hArea->SetMaximum(1.2);
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");
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;
2611 hSys->SetBinContent(ipt + 1, 1.);
2612 hSys->SetBinError(ipt + 1, max);
2614 hSys->Draw("same, E2");
2617 // for (Int_t idata = 0; idata < ndata; idata++)
2618 // delete hr[idata];
2624 TOFpid_systematics_bkgFit(Int_t ipart, Int_t icharge, Int_t icent)
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"
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"
2644 Int_t marker[ndata] = {22, 28, 22, 28, 22, 28};
2645 Int_t color[ndata] = {2, 2, 8, 8, 4, 4};
2647 TH1D *hArea = new TH1D("hArea", "", NptBins, ptBin);
2648 hArea->SetMinimum(0.5);
2649 hArea->SetMaximum(1.5);
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");
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;
2670 hSys->SetBinContent(ipt + 1, 1.);
2671 hSys->SetBinError(ipt + 1, max);
2673 hSys->Draw("same, E2");
2676 // for (Int_t idata = 0; idata < ndata; idata++)
2677 // delete hr[idata];
2683 TOFpid_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)
2686 TH1D *hr = new TH1D("hr", "", NptBins, ptBin);
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));
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));
2699 TFile *filein1 = TFile::Open(outfilename1);
2700 if (!filein1 || !filein1->IsOpen()) {
2701 printf("cannot open %s\n", outfilename1);
2704 TFile *filein2 = TFile::Open(outfilename2);
2705 if (!filein2 || !filein2->IsOpen()) {
2706 printf("cannot open %s\n", outfilename2);
2710 TH1D *h1 = (TH1D *)filein1->Get("RawSpectra/hNormalizedRawYield");
2712 printf("cannot get RawSpectra/hNormalizedRawYield from %s\n", outfilename1);
2715 TH1D *h2 = (TH1D *)filein2->Get("RawSpectra/hNormalizedRawYield");
2717 printf("cannot get RawSpectra/hNormalizedRawYield from %s\n", outfilename2);
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);