Correct initialization of static data members
[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 //_____________________________________________________________________________
57 AliTRDpidESD::AliTRDpidESD(const AliTRDpidESD &p):TObject(p)
58 {
59   //
60   // AliTRDpidESD copy constructor
61   //
62
63   ((AliTRDpidESD &) p).Copy(*this);
64
65 }
66
67 //_____________________________________________________________________________
68 AliTRDpidESD &AliTRDpidESD::operator=(const AliTRDpidESD &p)
69 {
70   //
71   // Assignment operator
72   //
73
74   if (this != &p) ((AliTRDpidESD &) p).Copy(*this);
75   return *this;
76
77 }
78
79 //_____________________________________________________________________________
80 void AliTRDpidESD::Copy(TObject &p) const
81 {
82   //
83   // Copy function
84   //
85
86  ((AliTRDpidESD &) p).fCheckTrackStatus          = fCheckTrackStatus;
87  ((AliTRDpidESD &) p).fCheckKinkStatus           = fCheckKinkStatus;
88  ((AliTRDpidESD &) p).fMinPlane                  = fMinPlane;
89
90 }
91
92 //_____________________________________________________________________________
93 Int_t AliTRDpidESD::MakePID(AliESD *event)
94 {
95   //
96   // This function calculates the PID probabilities based on TRD signals
97   //
98   // So far this method produces probabilities based on the total charge
99   // in each layer and the position of the maximum time bin in each layer.
100   // In a final version this should also exploit the charge measurement in
101   // the different slices of a given layer.
102   //  
103
104   Double_t p[10];
105   Int_t    nSpecies  = AliPID::kSPECIES;
106   Int_t    nPlanePID = 0;
107   Double_t mom       = 0.0;
108   Double_t probTotal = 0.0;
109
110   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
111   if (!calibration) {
112     //AliError("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     //AliError("No access to AliTRDCalPIDLQ\n");
120     return -1;
121   }
122
123   // Loop through all ESD tracks
124   Int_t ntrk = event->GetNumberOfTracks();
125   for (Int_t i = 0; i < ntrk; i++) {
126
127     AliESDtrack *t = event->GetTrack(i);
128
129     // Check the ESD track status
130     if (fCheckTrackStatus) {
131       if (((t->GetStatus() & AliESDtrack::kTRDout  ) == 0) &&
132           ((t->GetStatus() & AliESDtrack::kTRDrefit) == 0)) {
133         continue;
134       }
135     }
136
137     // Check for ESD kink tracks
138     if (fCheckKinkStatus) {
139       if (t->GetKinkIndex(0) != 0) {
140         continue;
141       }
142     }
143
144     // Skip tracks that have no TRD signal at all
145     if (t->GetTRDsignal() == 0) {
146       continue;
147     }
148
149     mom       = t->GetP();
150     probTotal = 0.0;
151     nPlanePID = 0;
152     for (Int_t iSpecies = 0; iSpecies < nSpecies; iSpecies++) {
153       p[iSpecies] = 1.0;
154     }
155
156     // Check the different detector layers
157     for (Int_t iPlan = 0; iPlan < AliTRDgeometry::kNplan; iPlan++) {
158
159       // Use the total charge in a given plane
160       Double_t dedx    = t->GetTRDsignals(iPlan,-1);
161       Int_t    timebin = t->GetTRDTimBin(iPlan);
162       if ((dedx    >  0.0) && 
163           (timebin > -1.0)) {
164
165         nPlanePID++;
166
167         // Get the probabilities for the different particle species
168         for (Int_t iSpecies = 0; iSpecies < nSpecies; iSpecies++) {
169
170           p[iSpecies] *= pd->GetProbability(iSpecies,mom,dedx);
171           p[iSpecies] *= pd->GetProbabilityT(iSpecies,mom,timebin);
172           p[iSpecies] *= 100.0; // ??????????????
173
174           probTotal   += p[iSpecies];
175
176         }
177
178       }
179
180     } 
181
182     for (Int_t iSpecies = 0; iSpecies < nSpecies; iSpecies++) {
183       if ((probTotal >       0.0) &&
184           (nPlanePID > fMinPlane)) {
185         p[iSpecies] /= probTotal;
186       }
187       else {
188         p[iSpecies] = 1.0;
189       }
190     }
191
192     t->SetTRDpid(p);
193
194   }
195
196   return 0;
197
198 }