Now dumping all data from simulation via MCDataInterface.
[u/mrichter/AliRoot.git] / TRD / AliTRDpidESD.cxx
CommitLineData
b0f03c34 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
dab811d0 16/* $Id$ */
17
b52cb1b6 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// //
dab811d0 28// Authors : //
29// Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de> (Original version)//
30// Alex Bercuci (a.bercuci@gsi.de) //
b52cb1b6 31// //
32////////////////////////////////////////////////////////////////////////////
b0f03c34 33
b52cb1b6 34#include "AliLog.h"
b0f03c34 35#include "AliESD.h"
36#include "AliESDtrack.h"
b52cb1b6 37
38#include "AliTRDpidESD.h"
39#include "AliTRDgeometry.h"
cc7cef99 40#include "AliTRDcalibDB.h"
dab811d0 41#include "AliRun.h"
42#include "AliTRDtrack.h"
7754cd1f 43#include "Cal/AliTRDCalPIDLQ.h"
b0f03c34 44
45ClassImp(AliTRDpidESD)
46
df34b21e 47 Bool_t AliTRDpidESD::fCheckTrackStatus = kTRUE;
48 Bool_t AliTRDpidESD::fCheckKinkStatus = kFALSE;
85815482 49 Int_t AliTRDpidESD::fMinPlane = 0;
df34b21e 50
b52cb1b6 51//_____________________________________________________________________________
dab811d0 52AliTRDpidESD::AliTRDpidESD()
53 :TObject()
b0f03c34 54{
55 //
b52cb1b6 56 // Default constructor
b0f03c34 57 //
b52cb1b6 58
b0f03c34 59}
60
b52cb1b6 61//_____________________________________________________________________________
dab811d0 62AliTRDpidESD::AliTRDpidESD(const AliTRDpidESD &p)
63 :TObject(p)
eab5961e 64{
b0f03c34 65 //
b52cb1b6 66 // AliTRDpidESD copy constructor
b0f03c34 67 //
eab5961e 68
b52cb1b6 69 ((AliTRDpidESD &) p).Copy(*this);
70
71}
72
73//_____________________________________________________________________________
74AliTRDpidESD &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//_____________________________________________________________________________
86void AliTRDpidESD::Copy(TObject &p) const
87{
88 //
89 // Copy function
90 //
91
25ca55ce 92 ((AliTRDpidESD &) p).fCheckTrackStatus = fCheckTrackStatus;
93 ((AliTRDpidESD &) p).fCheckKinkStatus = fCheckKinkStatus;
94 ((AliTRDpidESD &) p).fMinPlane = fMinPlane;
eab5961e 95
b0f03c34 96}
97
b52cb1b6 98//_____________________________________________________________________________
b0f03c34 99Int_t AliTRDpidESD::MakePID(AliESD *event)
100{
101 //
b52cb1b6 102 // This function calculates the PID probabilities based on TRD signals
b0f03c34 103 //
dab811d0 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;
b52cb1b6 126 }
127
b52cb1b6 128
dab811d0 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}
b52cb1b6 188
dab811d0 189//_____________________________________________________________________________
190Bool_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 }
b52cb1b6 201
dab811d0 202 // Check for ESD kink tracks
203 if (fCheckKinkStatus && (t->GetKinkIndex(0) != 0)) return kFALSE;
b52cb1b6 204
dab811d0 205 return kTRUE;
206}
b52cb1b6 207
dab811d0 208//_____________________________________________________________________________
209Bool_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
f162af62 223 AliTRDgeometry trdGeom;
dab811d0 224 const Float_t kAmHalfWidth = AliTRDgeometry::AmThick() / 2.;
f162af62 225 const Float_t kDrWidth = AliTRDgeometry::DrThick();
dab811d0 226
227
228 // retrive the magnetic field
229 Double_t xyz0[3] = { 0., 0., 0.}, xyz1[3];
230 Double_t b[3], alpha;
231 gAlice->Field(xyz0,b); // b[] is in kilo Gauss
232 Float_t field = b[2] * 0.1; // Tesla
dab811d0 233
234 // find momentum at chamber entrance and track length in chamber
235 AliExternalTrackParam *param = (plan<3) ? new AliExternalTrackParam(*t->GetInnerParam()) : new AliExternalTrackParam(*t->GetOuterParam());
236
f162af62 237 param->PropagateTo(trdGeom.GetTime0(plan)+kAmHalfWidth, field);
dab811d0 238 param->GetXYZ(xyz0);
239 alpha = param->GetAlpha();
f162af62 240 param->PropagateTo(trdGeom.GetTime0(plan)-kAmHalfWidth-kDrWidth, field);
dab811d0 241 // eliminate track segments which are crossing SM boundaries along chamber
aa612059 242 if(TMath::Abs(alpha-param->GetAlpha())>.01){
dab811d0 243 delete param;
dab811d0 244 return kFALSE;
245 }
246 param->GetXYZ(xyz1);
247 length = sqrt(
248 (xyz1[0]-xyz0[0])*(xyz1[0]-xyz0[0])+
249 (xyz1[1]-xyz0[1])*(xyz1[1]-xyz0[1])+
250 (xyz1[2]-xyz0[2])*(xyz1[2]-xyz0[2])
251 );
252 param->GetPxPyPz(xyz1);
253 mom = sqrt(xyz1[0]*xyz1[0] + xyz1[1]*xyz1[1] + xyz1[2]*xyz1[2]);
254 delete param;
dab811d0 255
256 return kTRUE;
b0f03c34 257}
dab811d0 258