Go back to previous version of simulation
[u/mrichter/AliRoot.git] / TRD / AliTRDpidESD.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
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  **************************************************************************/
15
16 ////////////////////////////////////////////////////////////////////////////
17 //                                                                        //
18 // Implementation of the TRD PID class                                    //
19 //                                                                        //
20 // Assigns the electron and pion likelihoods to each ESD track.           //
21 // The function MakePID(AliESD *event) calculates the probability         //
22 // of having dedx and a maximum timbin at a given                         //
23 // momentum (mom) and particle type k                                     //
24 // from the precalculated distributions.                                  //
25 //                                                                        //
26 // Original version:                                                      //
27 // Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de>                   //
28 //                                                                        //
29 ////////////////////////////////////////////////////////////////////////////
30
31 #include "AliLog.h"
32 #include "AliESD.h"
33 #include "AliESDtrack.h"
34
35 #include "AliTRDpidESD.h"
36 #include "AliTRDgeometry.h"
37 #include "AliTRDcalibDB.h"
38 #include "Cal/AliTRDCalPIDLQ.h"
39
40 ClassImp(AliTRDpidESD)
41
42   Bool_t AliTRDpidESD::fCheckTrackStatus = kTRUE;
43   Bool_t AliTRDpidESD::fCheckKinkStatus  = kFALSE;
44   Int_t AliTRDpidESD::fMinPlane          = 0;
45
46 //_____________________________________________________________________________
47 AliTRDpidESD::AliTRDpidESD():TObject()
48 {
49   //
50   // Default constructor
51   //
52
53 }
54
55 //_____________________________________________________________________________
56 AliTRDpidESD::AliTRDpidESD(const AliTRDpidESD &p):TObject(p)
57 {
58   //
59   // AliTRDpidESD copy constructor
60   //
61
62   ((AliTRDpidESD &) p).Copy(*this);
63
64 }
65
66 //_____________________________________________________________________________
67 AliTRDpidESD &AliTRDpidESD::operator=(const AliTRDpidESD &p)
68 {
69   //
70   // Assignment operator
71   //
72
73   if (this != &p) ((AliTRDpidESD &) p).Copy(*this);
74   return *this;
75
76 }
77
78 //_____________________________________________________________________________
79 void AliTRDpidESD::Copy(TObject &p) const
80 {
81   //
82   // Copy function
83   //
84
85   ((AliTRDpidESD &) p).fCheckTrackStatus          = fCheckTrackStatus;
86   ((AliTRDpidESD &) p).fCheckKinkStatus           = fCheckKinkStatus;
87   ((AliTRDpidESD &) p).fMinPlane                  = fMinPlane;
88
89 }
90
91 //_____________________________________________________________________________
92 Int_t AliTRDpidESD::MakePID(AliESD *event)
93 {
94   //
95   // This function calculates the PID probabilities based on TRD signals
96   //
97   // So far this method produces probabilities based on the total charge
98   // in each layer and the position of the maximum time bin in each layer.
99   // In a final version this should also exploit the charge measurement in
100   // the different slices of a given layer.
101   //  
102
103   Double_t p[10];
104   Int_t    nSpecies  = AliPID::kSPECIES;
105   Int_t    nPlanePID = 0;
106   Double_t mom       = 0.0;
107   Double_t probTotal = 0.0;
108
109   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
110   if (!calibration) {
111     AliErrorGeneral("AliTRDpidESD::MakePID"
112                    ,"No access to calibration data\n");
113     return -1;
114   }  
115
116   // Retrieve the CDB container class with the probability distributions
117   const AliTRDCalPIDLQ *pd = calibration->GetPIDLQObject();
118   if (!pd) {
119     AliErrorGeneral("AliTRDpidESD::MakePID"
120                    ,"No access to AliTRDCalPIDLQ\n");
121     return -1;
122   }
123
124   // Loop through all ESD tracks
125   Int_t ntrk = event->GetNumberOfTracks();
126   for (Int_t i = 0; i < ntrk; i++) {
127
128     AliESDtrack *t = event->GetTrack(i);
129
130     // Check the ESD track status
131     if (fCheckTrackStatus) {
132       if (((t->GetStatus() & AliESDtrack::kTRDout  ) == 0) &&
133           ((t->GetStatus() & AliESDtrack::kTRDrefit) == 0)) {
134         continue;
135       }
136     }
137
138     // Check for ESD kink tracks
139     if (fCheckKinkStatus) {
140       if (t->GetKinkIndex(0) != 0) {
141         continue;
142       }
143     }
144
145     // Skip tracks that have no TRD signal at all
146     if (t->GetTRDsignal() == 0) {
147       continue;
148     }
149
150     mom       = t->GetP();
151     probTotal = 0.0;
152     nPlanePID = 0;
153     for (Int_t iSpecies = 0; iSpecies < nSpecies; iSpecies++) {
154       p[iSpecies] = 1.0;
155     }
156
157     // Check the different detector layers
158     for (Int_t iPlan = 0; iPlan < AliTRDgeometry::kNplan; iPlan++) {
159
160       // Use the total charge in a given plane
161       Double_t dedx    = t->GetTRDsignals(iPlan,-1);
162       Int_t    timebin = t->GetTRDTimBin(iPlan);
163       if ((dedx    >  0.0) && 
164           (timebin > -1.0)) {
165
166         nPlanePID++;
167
168         // Get the probabilities for the different particle species
169         for (Int_t iSpecies = 0; iSpecies < nSpecies; iSpecies++) {
170
171           p[iSpecies] *= pd->GetProbability(iSpecies,mom,dedx);
172           p[iSpecies] *= pd->GetProbabilityT(iSpecies,mom,timebin);
173           p[iSpecies] *= 100.0; // ??????????????
174
175         }
176
177       }
178
179     }
180
181     for (Int_t iSpecies = 0; iSpecies < nSpecies; iSpecies++) {
182       probTotal += p[iSpecies];
183     }
184
185     for (Int_t iSpecies = 0; iSpecies < nSpecies; iSpecies++) {
186       if ((probTotal >       0.0) &&
187           (nPlanePID > fMinPlane)) {
188         p[iSpecies] /= probTotal;
189       }
190       else {
191         p[iSpecies] = 1.0;
192       }
193     }
194
195     t->SetTRDpid(p);
196
197   }
198
199   return 0;
200
201 }