DA rpms
[u/mrichter/AliRoot.git] / TOF / TestTOFPID.C
1 TestTOFPID(const Char_t *filename, Bool_t calibrateESD = kTRUE, Bool_t correctTExp = kTRUE, Bool_t useT0TOF = kTRUE, Double_t timeResolution = 100., Bool_t tuneTOFMC = kFALSE)
2 {
3   /* PID analysis */
4
5   /* check MC flag */
6   if (tuneTOFMC) calibrateESD = kFALSE;
7
8   /* init ESD */
9   TFile *filein = TFile::Open(filename);
10   TTree *treein = (TTree *)filein->Get("esdTree");
11   AliESDEvent *event = new AliESDEvent();
12   event->ReadFromTree(treein);
13   /* init OCDB */
14   treein->GetEvent(0);
15   Int_t run = event->GetRunNumber();
16   AliCDBManager *cdb = AliCDBManager::Instance();
17   cdb->SetDefaultStorage("raw://");
18   cdb->SetRun(run);
19   /* init TOF calibration */
20   AliTOFcalib *tofCalib = new AliTOFcalib();
21   if (correctTExp)
22     tofCalib->SetCorrectTExp(kTRUE);
23   tofCalib->Init();
24   /* init TOF T0-maker */
25   AliESDpid *fPIDesd = new AliESDpid();
26   AliTOFT0maker *t0maker = new AliTOFT0maker(fPIDesd, tofCalib); 
27   t0maker->SetTimeResolution(timeResolution);
28
29   /* pid histos */
30   TH2F *hTOFpid[AliPID::kSPECIES];
31   TH2F *hTOFpidSigma[AliPID::kSPECIES];
32   for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) {
33     hTOFpid[ipart] = new TH2F(Form("hTOFpid_%s", AliPID::ParticleName(ipart)), Form("%s-ID;p (GeV/c);t - t^{%s}_{exp} (ps)", AliPID::ParticleName(ipart), AliPID::ParticleLatexName(ipart)), 200, 0., 10., 2000, -24400., 24400.);
34     hTOFpidSigma[ipart] = new TH2F(Form("hTOFpidSigma_%s", AliPID::ParticleName(ipart)), Form("%s-ID;p (GeV/c);(t - t^{%s}_{exp}) / #sigma^{%s}_{exp} (ps)", AliPID::ParticleName(ipart), AliPID::ParticleLatexName(ipart), AliPID::ParticleLatexName(ipart)), 200, 0., 10., 2000, -10., 10.);
35   }
36
37   AliESDtrack *track;
38   Double_t p, time, timei[AliPID::kSPECIESC], sigma[AliPID::kSPECIES];
39   /* loop over events */
40   for (Int_t iev = 0; iev < treein->GetEntries(); iev++) {
41     
42     /* get event */
43     treein->GetEvent(iev);
44
45     /* calibrate ESD */
46     if (calibrateESD)
47       tofCalib->CalibrateESD(event);
48
49     /* tune TOF if requested for MC */
50     if (tuneTOFMC) {
51       t0maker->TuneForMC(event);
52     }
53
54     /* compute and apply T0-TOF */
55     if (useT0TOF) {
56       t0maker->ComputeT0TOF(event);
57       t0maker->ApplyT0TOF(event);
58       fPIDesd->MakePID(event,kFALSE,0.);
59     }
60
61     /* loop over tracks */
62     for (Int_t itrk = 0; itrk < event->GetNumberOfTracks(); itrk++) {
63       track = event->GetTrack(itrk);
64       /* check TOF match */
65       if (!track || 
66           !(track->GetStatus() & AliESDtrack::kTOFout) ||
67           !(track->GetStatus() & AliESDtrack::kTIME)) continue;
68
69       /* get track momentum */
70       p = track->P();
71       /* get TOF time */
72       time = track->GetTOFsignal();
73       /* get expected times */
74       track->GetIntegratedTimes(timei);
75
76       /* fill PID histos */
77       for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++){
78         sigma[ipart] = t0maker->GetExpectedSigma(p, timei[ipart], AliPID::ParticleMass(ipart));
79         hTOFpid[ipart]->Fill(p, (time - timei[ipart])); 
80         hTOFpidSigma[ipart]->Fill(p, (time - timei[ipart]) / sigma[ipart]); 
81       }
82     }
83   }
84   
85   /* write output */
86   TFile *fileout = TFile::Open("testPID.root", "RECREATE");
87   for (Int_t ipart = 0; ipart < 5; ipart++){
88     hTOFpid[ipart]->Write();
89     hTOFpidSigma[ipart]->Write();
90   }
91   fileout->Close();
92 }