Using TMath::Abs instead of fabs
[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 /* $Id$ */
17
18 ////////////////////////////////////////////////////////////////////////////
19 //                                                                        //
20 // Implementation of the TRD PID class                                    //
21 //                                                                        //
22 // Assigns the electron and pion likelihoods to each ESD track.           //
23 // The function MakePID(AliESD *event) calculates the probability         //
24 // of having dedx and a maximum timbin at a given                         //
25 // momentum (mom) and particle type k                                     //
26 // from the precalculated distributions.                                  //
27 //                                                                        //
28 // Authors :                                                              //
29 // Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de> (Original version)//
30 // Alex Bercuci (a.bercuci@gsi.de)                                        //
31 //                                                                        //
32 ////////////////////////////////////////////////////////////////////////////
33
34 #include "AliLog.h"
35 #include "AliESD.h"
36 #include "AliESDtrack.h"
37
38 #include "AliTRDpidESD.h"
39 #include "AliTRDgeometry.h"
40 #include "AliTRDcalibDB.h"
41 #include "AliRun.h"
42 #include "AliTRDtrack.h"
43 #include "Cal/AliTRDCalPIDLQ.h"
44
45 ClassImp(AliTRDpidESD)
46
47   Bool_t AliTRDpidESD::fCheckTrackStatus = kTRUE;
48   Bool_t AliTRDpidESD::fCheckKinkStatus  = kFALSE;
49   Int_t AliTRDpidESD::fMinPlane          = 0;
50
51 //_____________________________________________________________________________
52 AliTRDpidESD::AliTRDpidESD()
53   :TObject()
54 {
55   //
56   // Default constructor
57   //
58
59 }
60
61 //_____________________________________________________________________________
62 AliTRDpidESD::AliTRDpidESD(const AliTRDpidESD &p)
63   :TObject(p)
64 {
65   //
66   // AliTRDpidESD copy constructor
67   //
68
69   ((AliTRDpidESD &) p).Copy(*this);
70
71 }
72
73 //_____________________________________________________________________________
74 AliTRDpidESD &AliTRDpidESD::operator=(const AliTRDpidESD &p)
75 {
76   //
77   // Assignment operator
78   //
79
80   if (this != &p) ((AliTRDpidESD &) p).Copy(*this);
81   return *this;
82
83 }
84
85 //_____________________________________________________________________________
86 void AliTRDpidESD::Copy(TObject &p) const
87 {
88   //
89   // Copy function
90   //
91
92   ((AliTRDpidESD &) p).fCheckTrackStatus          = fCheckTrackStatus;
93   ((AliTRDpidESD &) p).fCheckKinkStatus           = fCheckKinkStatus;
94   ((AliTRDpidESD &) p).fMinPlane                  = fMinPlane;
95
96 }
97
98 //_____________________________________________________________________________
99 Int_t AliTRDpidESD::MakePID(AliESD *event)
100 {
101   //
102   // This function calculates the PID probabilities based on TRD signals
103   //
104   // The method produces probabilities based on the charge
105   // and the position of the maximum time bin in each layer.
106   // The dE/dx information can be used as global charge or 2 to 3
107   // slices. Check AliTRDCalPIDLQ and AliTRDCalPIDLQRef for the actual
108   // implementation.
109   //
110   // Author
111   // Alex Bercuci (A.Bercuci@gsi.de) 2nd May 2007
112
113         AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
114         if (!calibration) {
115                 AliErrorGeneral("AliTRDpidESD::MakePID()"
116                         ,"No access to calibration data");
117                 return -1;
118         }
119         
120         // Retrieve the CDB container class with the probability distributions
121         const AliTRDCalPIDLQ *pd = calibration->GetPIDLQObject();
122         if (!pd) {
123                 AliErrorGeneral("AliTRDpidESD::MakePID()"
124                         ,"No access to AliTRDCalPIDLQ");
125                 return -1;
126         }
127
128
129         // Loop through all ESD tracks
130         Double_t p[10];
131         AliESDtrack *t = 0x0;
132         Double_t dedx[AliTRDtrack::kNslice], dEdx;
133         Int_t    timebin;
134         Float_t mom, length, probTotal;
135         Int_t nPlanePID;
136         for (Int_t i=0; i<event->GetNumberOfTracks(); i++) {
137                 t = event->GetTrack(i);
138
139                 // Check track
140                 if(!CheckTrack(t)) continue;
141                                                 
142                 // Skip tracks which have no TRD signal at all
143                 if (t->GetTRDsignal() == 0.) continue;
144         
145                 // Loop over detector layers
146                 mom       = 0.; //t->GetP();
147                 length    = 0.;
148                 probTotal = 0.;
149                 nPlanePID    = 0;
150                 for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) p[iSpecies] = 1.;
151                 for (Int_t iPlan = 0; iPlan < AliTRDgeometry::kNplan; iPlan++) {
152                         // read data for track segment
153                         for(int iSlice=0; iSlice<AliTRDtrack::kNslice; iSlice++)
154                                 dedx[iSlice] = t->GetTRDsignals(iPlan, iSlice);
155                         dEdx    = t->GetTRDsignals(iPlan, -1);
156                         timebin = t->GetTRDTimBin(iPlan);
157
158                         // check data
159                         if ((dEdx <=  0.) || (timebin <= -1.)) continue;
160
161                         // retrive kinematic info for this track segment
162                         if(!GetTrackSegmentKine(t, iPlan, mom, length)) continue;
163                         
164                         // this track segment has fulfilled all requierments
165                         nPlanePID++;
166                         
167                         // Get the probabilities for the different particle species
168                         for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) {
169                                 p[iSpecies] *= pd->GetProbability(iSpecies, mom, dedx, length);
170                                 p[iSpecies] *= pd->GetProbabilityT(iSpecies, mom, timebin);
171                                 probTotal   += p[iSpecies];
172                         }
173                 }
174
175                 // normalize probabilities
176                 if(probTotal > 0.)
177                         for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
178                                 if(nPlanePID > fMinPlane) p[iSpecies] /= probTotal;
179                                 else p[iSpecies] = 1.0;
180
181
182                 // book PID to the track
183                 t->SetTRDpid(p);
184         }
185         
186         return 0;
187 }
188
189 //_____________________________________________________________________________
190 Bool_t AliTRDpidESD::CheckTrack(AliESDtrack *t)
191 {
192   //
193   // Check if track is eligible for PID calculations
194   //
195         
196         // Check the ESD track status
197         if (fCheckTrackStatus) {
198                 if (((t->GetStatus() & AliESDtrack::kTRDout  ) == 0) &&
199                         ((t->GetStatus() & AliESDtrack::kTRDrefit) == 0)) return kFALSE;
200         }
201
202         // Check for ESD kink tracks
203         if (fCheckKinkStatus && (t->GetKinkIndex(0) != 0)) return kFALSE;
204
205         return kTRUE;
206 }
207
208 //_____________________________________________________________________________
209 Bool_t AliTRDpidESD::GetTrackSegmentKine(AliESDtrack *t, Int_t plan, Float_t &mom, Float_t &length)
210 {
211   //
212   // Retrive momentum "mom" and track "length" in TRD chamber from plane
213   // "plan" according to information stored in AliESDtrack "t".
214   // 
215
216         if(!gAlice){
217                 AliErrorGeneral("AliTRDpidESD::GetTrackSegmentKine()"
218                 ,"No gAlice object to retrive TRDgeometry and Magnetic fied  - this has to be removed in the future.");
219                 return kFALSE;
220         }
221         
222         // Retrieve TRD geometry -> Maybe there is a better way to do this
223         Bool_t kSelfGeom = kFALSE;
224         AliTRDgeometry *TRDgeom =0x0;
225         if(gAlice) TRDgeom = AliTRDgeometry::GetGeometry(gAlice->GetRunLoader());
226         if(!TRDgeom){
227                 AliWarningGeneral("AliTRDpidESD::GetTrackSegmentKine()", "Cannot load TRD geometry from gAlice! Build a new one.\n");
228                 TRDgeom = new AliTRDgeometry();
229                 kSelfGeom = kTRUE;
230         }
231         const Float_t kAmHalfWidth = AliTRDgeometry::AmThick() / 2.;
232   const Float_t kDrWidth = AliTRDgeometry::DrThick();
233         
234
235         // retrive the magnetic field
236         Double_t xyz0[3] = { 0., 0., 0.}, xyz1[3];
237         Double_t b[3], alpha;
238         gAlice->Field(xyz0,b);      // b[] is in kilo Gauss
239         Float_t field = b[2] * 0.1; // Tesla
240
241                 
242         // find momentum at chamber entrance and track length in chamber
243         AliExternalTrackParam *param = (plan<3) ? new AliExternalTrackParam(*t->GetInnerParam()) : new AliExternalTrackParam(*t->GetOuterParam());
244
245         param->PropagateTo(TRDgeom->GetTime0(plan)+kAmHalfWidth, field);
246         param->GetXYZ(xyz0);
247         alpha = param->GetAlpha();
248         param->PropagateTo(TRDgeom->GetTime0(plan)-kAmHalfWidth-kDrWidth, field);
249         // eliminate track segments which are crossing SM boundaries along chamber
250         if(TMath::Abs(alpha-param->GetAlpha())>.01){
251                 delete param;
252                 if(kSelfGeom) delete TRDgeom;
253                 return kFALSE;
254         }
255         param->GetXYZ(xyz1);
256         length = sqrt(
257                 (xyz1[0]-xyz0[0])*(xyz1[0]-xyz0[0])+
258                 (xyz1[1]-xyz0[1])*(xyz1[1]-xyz0[1])+
259                 (xyz1[2]-xyz0[2])*(xyz1[2]-xyz0[2])
260         );
261         param->GetPxPyPz(xyz1);
262         mom = sqrt(xyz1[0]*xyz1[0] + xyz1[1]*xyz1[1] + xyz1[2]*xyz1[2]);
263         delete param;
264         if(kSelfGeom) delete TRDgeom;
265
266         return kTRUE;
267 }
268