1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
17 // Implements the abstract base class AliHFEpidBase
18 // IsInitialized() does the PID decision
21 // Markus Fasel <M.Fasel@gsi.de>
22 // Matus Kalisky <matus.kalisky@cern.ch> (contact)
27 #include <THnSparse.h>
28 #include <TDatabasePDG.h>
30 #include "AliAODTrack.h"
31 #include "AliAODMCParticle.h"
32 #include "AliESDtrack.h"
33 #include "AliMCParticle.h"
35 #include "AliESDpid.h"
37 #include "AliHFEcollection.h"
38 #include "AliHFEpidTOF.h"
39 #include "AliHFEpidBase.h"
42 ClassImp(AliHFEpidTOF)
44 //___________________________________________________________________
45 AliHFEpidTOF::AliHFEpidTOF(const Char_t *name):
55 //___________________________________________________________________
56 AliHFEpidTOF::AliHFEpidTOF(const AliHFEpidTOF &c):
68 //___________________________________________________________________
69 AliHFEpidTOF &AliHFEpidTOF::operator=(const AliHFEpidTOF &ref){
71 // Assignment operator
80 //___________________________________________________________________
81 AliHFEpidTOF::~AliHFEpidTOF(){
91 //___________________________________________________________________
92 void AliHFEpidTOF::Copy(TObject &ref) const {
94 // Performs the copying of the object
96 AliHFEpidTOF &target = dynamic_cast<AliHFEpidTOF &>(ref);
99 target.fQAList = fQAList;
100 target.fNsigmaTOF = fNsigmaTOF;
102 AliHFEpidBase::Copy(ref);
104 //___________________________________________________________________
105 Bool_t AliHFEpidTOF::InitializePID(){
107 // InitializePID: TOF experts have to implement code here
112 //___________________________________________________________________
113 Int_t AliHFEpidTOF::IsSelected(AliHFEpidObject *vtrack)
117 // as of 22/05/2006 :
118 // returns AliPID based on the ESD TOF PID decision
119 // the ESD PID will be checked and if necessary improved
120 // in the mean time. Best of luck
122 // returns 10 (== kUnknown)if PID can not be assigned for whatever reason
125 if(vtrack->fAnalysisType == AliHFEpidObject::kESDanalysis){
126 AliESDtrack *esdTrack = dynamic_cast<AliESDtrack *>(vtrack->fRecTrack);
127 if(!esdTrack) return 0;
128 AliMCParticle *mcTrack = dynamic_cast<AliMCParticle *>(vtrack->fMCtrack);
129 return MakePIDesdV3(esdTrack, mcTrack);
131 AliAODTrack *aodTrack = dynamic_cast<AliAODTrack *>(vtrack->fRecTrack);
132 if(!aodTrack) return 0;
133 AliAODMCParticle *aodmc = dynamic_cast<AliAODMCParticle *>(vtrack->fMCtrack);
134 return MakePIDaod(aodTrack, aodmc);
138 //___________________________________________________________________
139 Int_t AliHFEpidTOF::MakePIDesd(AliESDtrack *track, AliMCParticle * /*mcTrack*/){
141 // Does particle identification as discribed in IsSelected
144 AliError("No ESD PID object available");
148 status = track->GetStatus();
150 if(!(status & AliESDtrack::kTOFout)) return 0;
152 if(IsQAon()) fQAList->Fill("hTOF_flags", 0.);
154 Double_t tItrackL = track->GetIntegratedLength();
155 Double_t tTOFsignal = track->GetTOFsignal();
159 fQAList->Fill("hTOF_flags", 1.);
162 fQAList->Fill("hTOF_flags", 2.);
166 if(tItrackL <=0 || tTOFsignal <=0) return 0;
169 fQAList->Fill("hTOF_flags", 3.);
170 fQAList->Fill("hTOF_signal", tTOFsignal/1000.);
171 fQAList->Fill("hTOF_length", tItrackL);
173 // get the TOF pid probabilities
174 Double_t tESDpid[5] = {0., 0., 0., 0., 0.};
175 Float_t tTOFpidSum = 0.;
176 // find the largest PID probability
177 track->GetTOFpid(tESDpid);
178 Double_t tMAXpid = 0.;
179 Int_t tMAXindex = -1;
180 for(Int_t i=0; i<5; ++i){
181 tTOFpidSum += tESDpid[i];
182 if(tESDpid[i] > tMAXpid){
183 tMAXpid = tESDpid[i];
192 case 0: pdg = 11; specname = "electron"; break;
193 case 1: pdg = 13; specname = "muon"; break;
194 case 2: pdg = 211; specname = "pion"; break;
195 case 3: pdg = 321; specname = "kaon"; break;
196 case 4: pdg = 2212; specname = "proton"; break;
201 Double_t p = track->GetOuterParam()->P();
202 Double_t beta = (tItrackL/100.)/(TMath::C()*(tTOFsignal/1e12));
204 // sanity check, should not be necessary
205 if(TMath::Abs(tTOFpidSum - 1) > 0.01) return 0;
207 Double_t nSigmas = fESDpid->NumberOfSigmasTOF(track, (AliPID::EParticleType)tMAXindex, 0.);
208 if(TMath::Abs(nSigmas) > fNsigmaTOF) return 0;
211 // should be the same as AliPID flags
214 TString histname = "hTOFpid_" + specname;
215 fQAList->Fill(histname.Data(), beta, p);
216 fQAList->Fill("fTOFbeta_v_P_no", beta, p);
222 //__________________________________________________________________
223 Int_t AliHFEpidTOF::MakePIDesdV2(AliESDtrack *track, AliMCParticle * /*mcTrack*/){
225 // Computes the PID response based on TOF & T0 signal
229 AliError("No ESD PID object available");
233 status = track->GetStatus();
234 if(!(status & AliESDtrack::kTOFpid)) return 0;
236 Double_t p = track->GetOuterParam()->P();
237 // track integrated times for 5 different hypothesis and T0 time
239 track->GetIntegratedTimes(times);
240 Double_t tItrackL = track->GetIntegratedLength();
241 Double_t tTOFsignal = track->GetTOFsignal();
242 Double_t t0 = fESDpid->GetTOFResponse().GetTimeZero();
243 //printf("-D: tof: %f, T0: %f\n", tTOFsignal, t0);
244 // suppress missing or wrong T0 information
245 if(t0 > 999990.0) return 0;
246 Double_t tof = tTOFsignal - t0;
247 Double_t beta = (tItrackL/100.)/((TMath::C()*tof)/1e12);
248 //if(IsQAon())fQAList->Fill("hTOFbetaV2all", p, beta);
250 const Int_t pdg[5] = {11, 13, 211, 321, 2212};
251 const Double_t invMass[5] = {TDatabasePDG::Instance()->GetParticle(11)->Mass(),
252 TDatabasePDG::Instance()->GetParticle(13)->Mass(),
253 TDatabasePDG::Instance()->GetParticle(211)->Mass(),
254 TDatabasePDG::Instance()->GetParticle(321)->Mass(),
255 TDatabasePDG::Instance()->GetParticle(2212)->Mass()};
258 // accepted beta bands as function of momentum - parameters
259 // line: par[0]/p + par[1] + expected_tof
260 const Double_t bMin[5][2] = {{0., -0.03}, {-0.005, -0.02}, {-0.005, -0.02}, {-0.02, -0.006}, {-0.03, -0.005}};
261 const Double_t bMax[5][2] = {{0., 0.03}, {0.005, 0.02}, {0.005, 0.02}, {0.02, 0.006}, {0.03, 0.005}};
265 for(Int_t i=0; i<5; ++i){
266 Double_t d = (TMath::Abs(times[i] - tof))/times[i];
273 // stupid and unnecessary complicated - to be improved soon
274 Double_t a = p/(invMass[index]);
276 Double_t betaMatch = TMath::Sqrt(a/(1+a));
278 // check wheter the most probable match is within allowed region of beta for given momentum and species
279 Double_t min = bMin[index][0]/p + bMin[index][1] + betaMatch;
280 Double_t max = bMax[index][0]/p + bMax[index][1] + betaMatch;
283 //printf("-D: p: %f, beta: %f, pdg: %i, min: %f, max: %f\n", p, beta, pdg[index], min, max);
286 // PID decision - can be simpler than the QA histograms above could indicate !
290 if(beta < 0.2) return 0;
292 // 1) Simple version - protect electrons
293 if(beta > (1+bMin[0][1]) && beta < (1+bMax[0][1])){
294 //if(IsQAon())fQAList->Fill("hTOFbetaV2electron", p, beta);
300 // NOT ACTIVE when (1) activated
301 // 2) more complex version - return true PID of the particle based on the best TOF estimate
302 // under development - still keep protecting electrons
303 if(beta > (1+bMin[0][1]) && beta < (1+bMax[0][1])){
304 if(IsQAon())fQAList->Fill("hTOFbetaV2_electron", p, beta);
307 // above 3 GeV/c the supression gets weak
308 if(p > 3.0) return 0;
309 if(beta > min && beta < max) {
310 if(IsQAon())fQAList->Fill("hTOFbetaV2selected", p, beta);
318 //___________________________________________________________________
319 Int_t AliHFEpidTOF::MakePIDesdV3(AliESDtrack *track, AliMCParticle * /*mctrack*/){
321 // TOF PID based on n-Sigma cut
322 // Selects Protons and Kaons via n-sigma cut up to 3 GeV/c
323 // In addition histos for n-sigma before (all species) and after (only closest species) are filled
326 AliError("No ESD pid Object available. Return");
329 if(!(track->GetStatus() & AliESDtrack::kTOFpid)) return 0;
331 Double_t t0 = fESDpid->GetTOFResponse().GetTimeZero();
332 Double_t p = track->GetOuterParam() ? track->GetOuterParam()->P() : track->P();
334 // Fill before selection
335 Double_t sigEle = fESDpid->NumberOfSigmasTOF(track, AliPID::kElectron, t0);
336 //printf("-D: p: %f, t0: %f, nSigma: %f\n", p, t0, sigEle);
338 if(TMath::Abs(sigEle) < fNsigmaTOF)
342 Double_t nsigTPC = fESDpid->NumberOfSigmasTPC(track, AliPID::kElectron);
343 if(nsigTPC > -3 && nsigTPC < 5) fQAList->Fill("hTOFsigmaTPCcleanup", p, sigEle);
345 Double_t hcontent[3] = {p, sigEle, 0};
346 THnSparseF * hptr = dynamic_cast<THnSparseF *>(fQAList->Get("hTOFsigmaElectron"));
347 hptr->Fill(hcontent);
350 hptr->Fill(hcontent);
356 //___________________________________________________________________
357 Double_t AliHFEpidTOF::Likelihood(const AliESDtrack *track, Int_t species, Float_t rsig){
359 //gives probability for track to come from a certain particle species;
360 //no priors considered -> probabilities for equal abundances of all species!
361 //optional restriction to r-sigma bands of different particle species; default: rsig = 2. (see below)
363 //IMPORTANT: Tracks which are judged to be outliers get negative likelihoods -> unusable for combination with further detectors!
365 if(!track) return -1.;
366 Bool_t outlier = kTRUE;
367 // Check whether distance from the respective particle line is smaller than r sigma
368 for(Int_t hypo = 0; hypo < AliPID::kSPECIES; hypo++){
369 if(TMath::Abs(fESDpid->NumberOfSigmasTOF(track, (AliPID::EParticleType)hypo, 0.)) > rsig)
381 track->GetTOFpid(tofProb);
383 return tofProb[species];
385 //___________________________________________________________________
386 Int_t AliHFEpidTOF::MakePIDaod(AliAODTrack * /*aodTrack*/, AliAODMCParticle * /*mcTrack*/){
387 AliError("AOD PID not yet implemented");
391 //___________________________________________________________________
392 void AliHFEpidTOF::AddQAhistograms(TList *qaList){
394 // Create QA histograms for TOF PID
397 fQAList = new AliHFEcollection("TOFqaHistos", "Collection for TOF PID histograms");
398 //fQAList->SetName("fTOFqaHistos");
399 fQAList->CreateTH1F("hTOF_flags", "TOF flags;flags (see code for info);counts", 10, -0.25, 4.75);
400 fQAList->CreateTH2F("fTOFbeta_v_P_no","beta -v- P; beta;momentum [GeV/c]", 120, 0, 1.2, 200, 0, 20);
401 fQAList->CreateTH1F("hTOF_signal", "TOF signal; TOF signal [ns];counts", 1000, 12, 50);
402 fQAList->CreateTH1F("hTOF_length", "TOF track length; length [cm];counts", 400, 300, 700);
403 fQAList->CreateTH2F("hTOFpid_electron", "TOF reco electron; beta ; momentum [GeV/c]", 120, 0, 1.2, 200, 0, 5);
404 fQAList->CreateTH2F("hTOFpid_muon", "TOF reco muon; beta ; momentum [GeV/c]", 120, 0, 1.2, 200, 0, 5);
405 fQAList->CreateTH2F("hTOFpid_pion", "TOF reco pion; beta ; momentum [GeV/c]", 120, 0, 1.2, 200, 0, 5);
406 fQAList->CreateTH2F("hTOFpid_kaon", "TOF reco kaon; beta ; momentum [GeV/c]", 120, 0, 1.2, 200, 0, 5);
407 fQAList->CreateTH2F("hTOFpid_proton", "TOF reco proton; beta ; momentum [GeV/c]", 120, 0, 1.2, 200, 0, 5);
408 //fQAList->CreateTH2F("hTOFbetaV2all", "TOF #beta vs p for all tracks; momentum [GeV/c]; beta ", 400, 0.1, 10., 1200, 0, 1.2, 0);
409 //fQAList->CreateTH2F("hTOFbetaV2electron", "TOF #beta vs p for selected electron tracks; momentum [GeV/c]; beta", 400, 0.1, 10., 1200, 0, 1.2, 0);
411 // histograms for sigma cut
412 const Int_t kNdim= 3;
413 Int_t nBins[kNdim] = {1000, 1400, 2};
414 Double_t binMin[kNdim] = {0.1, -12., 0};
415 Double_t binMax[kNdim] = {20, 12., 2};
416 fQAList->CreateTHnSparse("hTOFsigmaElectron", "TOF N#sigma around the Electron Line for all tracks; p (GeV/c); N sigma; Selection Step", kNdim, nBins, binMin, binMax);
417 fQAList->CreateTH2F("hTOFsigmaTPCcleanup", "TOF N#sigma around the Electron Line for all tracks after TPC cleanup; p (GeV/c); N sigma;", 100, 0.1, 20, 140, -12, 12);
419 qaList->AddLast(fQAList->GetList());