Coding rules (Markus)
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEpidTRD.cxx
CommitLineData
809a4336 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**************************************************************************/
50685501 15//
16// Class for TRD PID
17// Implements the abstract base class AliHFEpidbase
18// Make PID does the PID decision
19// Class further contains TRD specific cuts and QA histograms
20//
21// Authors:
22// Markus Fasel <M.Fasel@gsi.de>
23//
809a4336 24#include <TH1F.h>
75d81601 25#include <TH2F.h>
50685501 26#include <TList.h>
809a4336 27#include <TString.h>
809a4336 28
722347d8 29#include "AliAODTrack.h"
30#include "AliAODMCParticle.h"
809a4336 31#include "AliESDtrack.h"
75d81601 32#include "AliLog.h"
722347d8 33#include "AliMCParticle.h"
809a4336 34#include "AliPID.h"
35
36#include "AliHFEpidTRD.h"
37
38ClassImp(AliHFEpidTRD)
39
75d81601 40const Double_t AliHFEpidTRD::fgkVerySmall = 1e-12;
41
809a4336 42//___________________________________________________________________
43AliHFEpidTRD::AliHFEpidTRD(const char* name) :
44 AliHFEpidBase(name)
809a4336 45 , fPIDMethod(kNN)
75d81601 46 , fContainer(0x0)
809a4336 47{
48 //
49 // default constructor
50 //
dbe3abbe 51 memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
809a4336 52}
53
54//___________________________________________________________________
55AliHFEpidTRD::AliHFEpidTRD(const AliHFEpidTRD &ref):
56 AliHFEpidBase("")
809a4336 57 , fPIDMethod(kLQ)
75d81601 58 , fContainer(0x0)
809a4336 59{
60 //
61 // Copy constructor
62 //
dbe3abbe 63 memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
809a4336 64 ref.Copy(*this);
65}
66
67//___________________________________________________________________
68AliHFEpidTRD &AliHFEpidTRD::operator=(const AliHFEpidTRD &ref){
69 //
70 // Assignment operator
71 //
72 if(this != &ref){
73 ref.Copy(*this);
74 }
75 return *this;
76}
77
78//___________________________________________________________________
79void AliHFEpidTRD::Copy(TObject &ref) const {
80 //
81 // Performs the copying of the object
82 //
83 AliHFEpidTRD &target = dynamic_cast<AliHFEpidTRD &>(ref);
84
809a4336 85 target.fPIDMethod = fPIDMethod;
dbe3abbe 86 memcpy(target.fThreshParams, fThreshParams, sizeof(Double_t) * kThreshParams);
75d81601 87 if(fContainer) target.fContainer = dynamic_cast<TList *>(fContainer->Clone());
809a4336 88 AliHFEpidBase::Copy(ref);
89}
90
91//___________________________________________________________________
92AliHFEpidTRD::~AliHFEpidTRD(){
93 //
94 // Destructor
95 //
75d81601 96 if(fContainer){
97 fContainer->Clear();
98 delete fContainer;
99 }
809a4336 100}
101
102//______________________________________________________
103Bool_t AliHFEpidTRD::InitializePID(){
104 //
105 // InitializePID: Load TRD thresholds and create the electron efficiency axis
106 // to navigate
107 //
dbe3abbe 108 InitParameters();
809a4336 109 return kTRUE;
110}
111
112//______________________________________________________
722347d8 113Int_t AliHFEpidTRD::IsSelected(AliHFEpidObject *track){
809a4336 114 //
115 // Does PID for TRD alone:
116 // PID thresholds based on 90% Electron Efficiency level approximated by a linear
117 // step function
118 //
722347d8 119 if(track->fAnalysisType == AliHFEpidObject::kESDanalysis){
120 AliESDtrack *esdTrack = dynamic_cast<AliESDtrack *>(track->fRecTrack);
121 if(!esdTrack) return 0;
122 AliMCParticle *esdmc = dynamic_cast<AliMCParticle *>(track->fMCtrack);
123 return MakePIDesd(esdTrack, esdmc);
124 } else {
125 AliAODTrack *aodTrack = dynamic_cast<AliAODTrack *>(track->fRecTrack);
126 if(!aodTrack) return 0;
127 AliAODMCParticle *aodmc = dynamic_cast<AliAODMCParticle *>(track->fMCtrack);
128 return MakePIDaod(aodTrack, aodmc);
129 }
130}
131
132//______________________________________________________
133Int_t AliHFEpidTRD::MakePIDaod(AliAODTrack * /*aodTrack*/, AliAODMCParticle * /*aodmc*/){
50685501 134 //
135 // Does PID decision for AOD tracks as discribed above
136 //
722347d8 137 AliError("AOD PID not yet implemented");
138 return 0;
139}
140
141//______________________________________________________
142Int_t AliHFEpidTRD::MakePIDesd(AliESDtrack *esdTrack, AliMCParticle * /*mcTrack*/){
50685501 143 //
144 // Does PID decision for ESD tracks as discribed above
145 //
809a4336 146 Double_t p = esdTrack->GetOuterParam() ? esdTrack->GetOuterParam()->P() : esdTrack->P();
722347d8 147 if(p < 2.0) return 0;
809a4336 148
809a4336 149 Double_t pidProbs[AliPID::kSPECIES];
150 esdTrack->GetTRDpid(pidProbs);
0792aa82 151 if(IsQAon())FillHistogramsLikelihood(0, p, pidProbs[AliPID::kElectron]);
152 Double_t threshold = GetTRDthresholds(0.91, p);
153 AliDebug(1, Form("Threshold: %f\n", threshold));
154 if(IsQAon()) (dynamic_cast<TH2F *>(fContainer->At(kHistTRDthresholds)))->Fill(p, threshold);
155 if(pidProbs[AliPID::kElectron] > threshold){
156 if(IsQAon()) FillHistogramsLikelihood(1, p, pidProbs[AliPID::kElectron]);
157 return 11;
158 }
159 return 211;
809a4336 160}
161
162//___________________________________________________________________
dbe3abbe 163Double_t AliHFEpidTRD::GetTRDthresholds(Double_t electronEff, Double_t p){
164 //
165 // Return momentum dependent and electron efficiency dependent TRD thresholds
166 //
167 Double_t params[4];
168 GetParameters(electronEff, params);
0792aa82 169 Double_t threshold = 1. - params[0] - params[1] * p - params[2] * TMath::Exp(-params[3] * p);
170 return TMath::Max(TMath::Min(threshold, 0.99), 0.2); // truncate the threshold upperwards to 0.999 and lowerwards to 0.2 and exclude unphysical values
809a4336 171}
172
173//___________________________________________________________________
dbe3abbe 174void AliHFEpidTRD::InitParameters(){
175 //
176 // Fill the Parameters into an array
177 //
178
179 // Parameters for 6 Layers
180 fThreshParams[0] = -0.001839; // 0.7 electron eff
181 fThreshParams[1] = 0.000276;
182 fThreshParams[2] = 0.044902;
183 fThreshParams[3] = 1.726751;
184 fThreshParams[4] = -0.002405; // 0.75 electron eff
185 fThreshParams[5] = 0.000372;
186 fThreshParams[6] = 0.061775;
187 fThreshParams[7] = 1.739371;
188 fThreshParams[8] = -0.003178; // 0.8 electron eff
189 fThreshParams[9] = 0.000521;
190 fThreshParams[10] = 0.087585;
191 fThreshParams[11] = 1.749154;
192 fThreshParams[12] = -0.004058; // 0.85 electron eff
193 fThreshParams[13] = 0.000748;
194 fThreshParams[14] = 0.129583;
195 fThreshParams[15] = 1.782323;
196 fThreshParams[16] = -0.004967; // 0.9 electron eff
197 fThreshParams[17] = 0.001216;
198 fThreshParams[18] = 0.210128;
199 fThreshParams[19] = 1.807665;
200 fThreshParams[20] = -0.000996; // 0.95 electron eff
201 fThreshParams[21] = 0.002627;
202 fThreshParams[22] = 0.409099;
203 fThreshParams[23] = 1.787076;
809a4336 204}
205
206//___________________________________________________________________
dbe3abbe 207void AliHFEpidTRD::GetParameters(Double_t electronEff, Double_t *parameters){
208 //
209 // return parameter set for the given efficiency bin
210 //
0792aa82 211 Int_t effbin = static_cast<Int_t>((electronEff - 0.7)/0.05);
dbe3abbe 212 memcpy(parameters, fThreshParams + effbin * 4, sizeof(Double_t) * 4);
809a4336 213}
75d81601 214
215//___________________________________________________________________
722347d8 216Double_t AliHFEpidTRD::GetTRDSignalV1(AliESDtrack *track, Int_t mcPID){
75d81601 217 //
218 // Calculation of the TRD Signal via truncated mean
219 // Method 1: Take all Slices available
220 // cut out 0s
221 // Order them in increasing order
222 // Cut out the upper third
223 // Calculate mean over the last 2/3 slices
224 //
225 const Int_t kNSlices = 48;
226 AliDebug(2, Form("Number of Tracklets: %d\n", track->GetTRDpidQuality()));
227 Double_t trdSlices[kNSlices], tmp[kNSlices];
228 Int_t indices[48];
229 Int_t icnt = 0;
230 for(Int_t idet = 0; idet < 6; idet++)
231 for(Int_t islice = 0; islice < 8; islice++){
232 AliDebug(2, Form("Chamber[%d], Slice[%d]: TRDSlice = %f", idet, islice, track->GetTRDslice(idet, islice)));
233 if(TMath::Abs(track->GetTRDslice(idet, islice)) < fgkVerySmall) continue;;
234 trdSlices[icnt++] = track->GetTRDslice(idet, islice);
235 }
236 AliDebug(1, Form("Number of Slices: %d\n", icnt));
237 if(icnt < 6) return 0.; // We need at least 6 Slices for the truncated mean
238 TMath::Sort(icnt, trdSlices, indices, kFALSE);
239 memcpy(tmp, trdSlices, sizeof(Double_t) * icnt);
240 for(Int_t ien = 0; ien < icnt; ien++)
241 trdSlices[ien] = tmp[indices[ien]];
242 Double_t trdSignal = TMath::Mean(static_cast<Int_t>(static_cast<Double_t>(icnt) * 2./3.), trdSlices);
243 Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1;
244 AliDebug(1, Form("PID Meth. 1: p[%f], TRDSignal[%f]", mom, trdSignal));
722347d8 245 if(IsQAon() && mom > 0.) FillHistogramsTRDSignalV1(trdSignal, mom, mcPID);
75d81601 246 return trdSignal;
247}
248
249//___________________________________________________________________
722347d8 250Double_t AliHFEpidTRD::GetTRDSignalV2(AliESDtrack *track, Int_t mcPID){
75d81601 251 //
252 // Calculation of the TRD Signal via truncated mean
253 // Method 2: Take only first 5 slices per chamber
254 // Order them in increasing order
255 // Cut out upper half
256 // Now do mean with the reamining 3 slices per chamber
257 //
258 Double_t trdSlicesLowTime[30], trdSlicesRemaining[32];
259 Int_t indices[30];
260 Int_t cntLowTime=0, cntRemaining = 0;
261 for(Int_t idet = 0; idet < 6; idet++)
262 for(Int_t islice = 0; islice < 8; islice++){
263 if(TMath::Abs(track->GetTRDslice(idet, islice)) < fgkVerySmall) continue;;
264 if(islice < 5){
265 AliDebug(2, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
266 trdSlicesLowTime[cntLowTime++] = track->GetTRDslice(idet, islice);
267 } else{
268 AliDebug(2, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
269 trdSlicesRemaining[cntRemaining++] = track->GetTRDslice(idet, islice);
270 }
271 }
272 if(cntLowTime < 4 || cntRemaining < 2) return 0.; // Min. Number of Slices at high time is 2 (matches with 1 layer), for the truncated mean we need at least 4 Slices
273 TMath::Sort(cntLowTime, trdSlicesLowTime, indices, kFALSE);
274 // Fill the second array with the lower half of the first time bins
275 for(Int_t ien = 0; ien < static_cast<Int_t>(static_cast<Double_t>(cntLowTime) * 0.5); ien++)
276 trdSlicesRemaining[cntRemaining++] = trdSlicesLowTime[indices[ien]];
277 Double_t trdSignal = TMath::Mean(cntRemaining, trdSlicesRemaining);
278 Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1;
279 AliDebug(1, Form("PID Meth. 2: p[%f], TRDSignal[%f]", mom, trdSignal));
722347d8 280 if(IsQAon() && mom > 0.) FillHistogramsTRDSignalV2(trdSignal, mom, mcPID);
75d81601 281 return trdSignal;
282}
283
284//___________________________________________________________________
285void AliHFEpidTRD::FillHistogramsTRDSignalV1(Double_t signal, Double_t mom, Int_t species){
286 //
287 // Fill histograms for TRD Signal from Method 1 vs. p for different particle species
288 //
289 (dynamic_cast<TH2F *>(fContainer->At(kHistTRDSigV1)))->Fill(mom, signal);
290 if(species < 0 || species >= AliPID::kSPECIES) return;
291 (dynamic_cast<TH2F *>(fContainer->At(kHistOverallSpecies + species)))->Fill(mom, signal);
292}
293
294//___________________________________________________________________
295void AliHFEpidTRD::FillHistogramsTRDSignalV2(Double_t signal, Double_t mom, Int_t species){
296 //
297 // Fill histograms for TRD Signal from Method 2 vs. p for different particle species
298 //
299 (dynamic_cast<TH2F *>(fContainer->At(kHistTRDSigV2)))->Fill(mom, signal);
300 if(species < 0 || species >= AliPID::kSPECIES) return;
301 (dynamic_cast<TH2F *>(fContainer->At(kHistOverallSpecies + AliPID::kSPECIES + species)))->Fill(mom, signal);
302}
303
304//___________________________________________________________________
75d81601 305void AliHFEpidTRD::AddQAhistograms(TList *l){
306 //
307 // Adding QA histograms for the TRD PID
308 // QA histograms are:
309 // + TRD Signal from Meth. 1 vs p for all species
310 // + TRD Signal from Meth. 2 vs p for all species
311 // + For each species
312 // - TRD Signal from Meth. 1 vs p
313 // - TRD Signal from Meth. 2 vs p
314 //
315 const Int_t kMomentumBins = 41;
316 const Double_t kPtMin = 0.1;
317 const Double_t kPtMax = 10.;
318 const Int_t kSigBinsMeth1 = 100;
319 const Int_t kSigBinsMeth2 = 100;
320 const Double_t kMinSig = 0.;
321 const Double_t kMaxSigMeth1 = 10000.;
322 const Double_t kMaxSigMeth2 = 10000.;
323
324 if(!fContainer) fContainer = new TList;
325 fContainer->SetName("fTRDqaHistograms");
326
327 Double_t momentumBins[kMomentumBins];
328 for(Int_t ibin = 0; ibin < kMomentumBins; ibin++)
329 momentumBins[ibin] = static_cast<Double_t>(TMath::Power(10,TMath::Log10(kPtMin) + (TMath::Log10(kPtMax)-TMath::Log10(kPtMin))/(kMomentumBins-1)*static_cast<Double_t>(ibin)));
0792aa82 330 // Likelihood Histograms
331 fContainer->AddAt(new TH2F("fTRDlikeBefore", "TRD Electron Likelihood before cut", kMomentumBins - 1, momentumBins, 1000, 0, 1), kHistTRDlikeBefore);
332 fContainer->AddAt(new TH2F("fTRDlikeAfter", "TRD Electron Likelihood after cut", kMomentumBins - 1, momentumBins, 1000, 0, 1), kHistTRDlikeAfter);
333 fContainer->AddAt(new TH2F("fTRDthesholds", "TRD Electron thresholds", kMomentumBins - 1, momentumBins, 1000, 0, 1), kHistTRDthresholds);
334 // Signal Histograms
75d81601 335 fContainer->AddAt(new TH2F("fTRDSigV1all", "TRD Signal (all particles, Method 1)", kMomentumBins - 1, momentumBins, kSigBinsMeth1, kMinSig, kMaxSigMeth1), kHistTRDSigV1);
336 fContainer->AddAt(new TH2F("fTRDSigV2all", "TRD Signal (all particles, Method 2)", kMomentumBins - 1, momentumBins, kSigBinsMeth2, kMinSig, kMaxSigMeth2), kHistTRDSigV2);
337 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
338 fContainer->AddAt(new TH2F(Form("fTRDSigV1%s", AliPID::ParticleName(ispec)), Form("TRD Signal (%s, Method 1)", AliPID::ParticleName(ispec)), kMomentumBins - 1, momentumBins, kSigBinsMeth1, kMinSig, kMaxSigMeth1), kHistOverallSpecies + ispec);
339 fContainer->AddAt(new TH2F(Form("fTRDSigV2%s", AliPID::ParticleName(ispec)), Form("TRD Signal (%s, Method 2)", AliPID::ParticleName(ispec)), kMomentumBins - 1, momentumBins, kSigBinsMeth1, kMinSig, kMaxSigMeth2), kHistOverallSpecies + AliPID::kSPECIES + ispec);
340 }
341 l->AddLast(fContainer);
342}
343
0792aa82 344//___________________________________________________________________
345void AliHFEpidTRD::FillHistogramsLikelihood(Int_t whenFilled, Float_t p, Float_t elProb){
346 //
347 // Fill Likelihood Histogram before respectively after decision
348 //
349 TH2F *histo = NULL;
350 if(whenFilled)
351 histo = dynamic_cast<TH2F *>(fContainer->At(kHistTRDlikeAfter));
352 else
353 histo = dynamic_cast<TH2F *>(fContainer->At(kHistTRDlikeBefore));
354 if(!histo){
355 AliError("QA histograms not found");
356 return;
357 }
358 histo->Fill(p, elProb);
359}