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