Cleanup the code. Fix memory leak. Now inherit from AliAnalysisTaskSE (Antoine, Phili...
[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**************************************************************************/
27de2dfb 15
16/* $Id$ */
17
50685501 18//
19// Class for TRD PID
20// Implements the abstract base class AliHFEpidbase
21// Make PID does the PID decision
22// Class further contains TRD specific cuts and QA histograms
23//
24// Authors:
25// Markus Fasel <M.Fasel@gsi.de>
26//
809a4336 27#include <TH1F.h>
75d81601 28#include <TH2F.h>
e3fc062d 29#include <THnSparse.h>
50685501 30#include <TList.h>
809a4336 31#include <TString.h>
809a4336 32
3a72645a 33#include "AliAODpidUtil.h"
6555e2ad 34#include "AliAODPid.h"
722347d8 35#include "AliAODTrack.h"
36#include "AliAODMCParticle.h"
809a4336 37#include "AliESDtrack.h"
e3fc062d 38#include "AliESDpid.h"
75d81601 39#include "AliLog.h"
722347d8 40#include "AliMCParticle.h"
809a4336 41#include "AliPID.h"
42
3a72645a 43#include "AliHFEpidQAmanager.h"
809a4336 44#include "AliHFEpidTRD.h"
45
46ClassImp(AliHFEpidTRD)
47
75d81601 48const Double_t AliHFEpidTRD::fgkVerySmall = 1e-12;
49
809a4336 50//___________________________________________________________________
faee3b18 51AliHFEpidTRD::AliHFEpidTRD() :
52 AliHFEpidBase()
e3fc062d 53 , fMinP(1.)
67fe7bd0 54 , fElectronEfficiency(0.91)
faee3b18 55 , fPIDMethod(kNN)
faee3b18 56{
57 //
58 // default constructor
59 //
bf892a6a 60 memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
61 SetUseDefaultParameters();
faee3b18 62}
63
64//___________________________________________________________________
809a4336 65AliHFEpidTRD::AliHFEpidTRD(const char* name) :
66 AliHFEpidBase(name)
e3fc062d 67 , fMinP(1.)
67fe7bd0 68 , fElectronEfficiency(0.91)
809a4336 69 , fPIDMethod(kNN)
809a4336 70{
71 //
72 // default constructor
73 //
dbe3abbe 74 memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
bf892a6a 75 SetUseDefaultParameters();
809a4336 76}
77
78//___________________________________________________________________
79AliHFEpidTRD::AliHFEpidTRD(const AliHFEpidTRD &ref):
80 AliHFEpidBase("")
67fe7bd0 81 , fMinP(ref.fMinP)
82 , fElectronEfficiency(ref.fElectronEfficiency)
83 , fPIDMethod(ref.fPIDMethod)
809a4336 84{
85 //
86 // Copy constructor
87 //
dbe3abbe 88 memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
809a4336 89 ref.Copy(*this);
90}
91
92//___________________________________________________________________
93AliHFEpidTRD &AliHFEpidTRD::operator=(const AliHFEpidTRD &ref){
94 //
95 // Assignment operator
96 //
97 if(this != &ref){
98 ref.Copy(*this);
99 }
100 return *this;
101}
102
103//___________________________________________________________________
104void AliHFEpidTRD::Copy(TObject &ref) const {
105 //
106 // Performs the copying of the object
107 //
108 AliHFEpidTRD &target = dynamic_cast<AliHFEpidTRD &>(ref);
bf892a6a 109
110 Bool_t defaultParameters = UseDefaultParameters();
111 target.SetUseDefaultParameters(defaultParameters);
e3fc062d 112 target.fMinP = fMinP;
809a4336 113 target.fPIDMethod = fPIDMethod;
67fe7bd0 114 target.fElectronEfficiency = fElectronEfficiency;
dbe3abbe 115 memcpy(target.fThreshParams, fThreshParams, sizeof(Double_t) * kThreshParams);
809a4336 116 AliHFEpidBase::Copy(ref);
117}
118
119//___________________________________________________________________
120AliHFEpidTRD::~AliHFEpidTRD(){
121 //
122 // Destructor
123 //
809a4336 124}
125
126//______________________________________________________
127Bool_t AliHFEpidTRD::InitializePID(){
128 //
129 // InitializePID: Load TRD thresholds and create the electron efficiency axis
130 // to navigate
131 //
bf892a6a 132 if(UseDefaultParameters()){
133 if(fPIDMethod == kLQ)
134 InitParameters1DLQ();
135 else
136 InitParameters();
137 }
809a4336 138 return kTRUE;
139}
140
141//______________________________________________________
6555e2ad 142Int_t AliHFEpidTRD::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager *pidqa) const {
809a4336 143 //
144 // Does PID for TRD alone:
145 // PID thresholds based on 90% Electron Efficiency level approximated by a linear
146 // step function
147 //
3a72645a 148 AliDebug(1, "Applying TRD PID");
149 if((!fESDpid && track->IsESDanalysis()) || (!fAODpid && track->IsAODanalysis())){
150 AliDebug(1, "Cannot process track");
151 return 0;
152 }
722347d8 153
3a72645a 154 AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis: AliHFEpidObject::kAODanalysis;
155 Double_t p = GetP(track->GetRecTrack(), anatype);
156 if(p < fMinP){
157 AliDebug(1, Form("Track momentum below %f", fMinP));
158 return 0;
159 }
809a4336 160
3a72645a 161 if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kTRDpid, AliHFEdetPIDqa::kBeforePID);
162 Double_t electronLike = GetElectronLikelihood(track->GetRecTrack(), anatype);
67fe7bd0 163 Double_t threshold = GetTRDthresholds(fElectronEfficiency, p);
0792aa82 164 AliDebug(1, Form("Threshold: %f\n", threshold));
3a72645a 165 if(electronLike > threshold){
166 if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kTRDpid, AliHFEdetPIDqa::kAfterPID);
0792aa82 167 return 11;
168 }
169 return 211;
3a72645a 170
809a4336 171}
172
173//___________________________________________________________________
6555e2ad 174Double_t AliHFEpidTRD::GetTRDthresholds(Double_t electronEff, Double_t p) const {
dbe3abbe 175 //
176 // Return momentum dependent and electron efficiency dependent TRD thresholds
177 //
178 Double_t params[4];
179 GetParameters(electronEff, params);
0792aa82 180 Double_t threshold = 1. - params[0] - params[1] * p - params[2] * TMath::Exp(-params[3] * p);
181 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 182}
183
184//___________________________________________________________________
bf892a6a 185void AliHFEpidTRD::SetThresholdParameters(Double_t electronEff, Double_t *params){
186 //
187 // Set threshold parameters for the given bin
188 //
189 if(electronEff >= 1. || electronEff < 0.7) return;
190 Int_t effbin = static_cast<Int_t>((electronEff - 0.7)/0.05);
191 memcpy(&fThreshParams[effbin * 4], params, sizeof(Double_t) * 4);
192 SetUseDefaultParameters(kFALSE);
193}
194
195//___________________________________________________________________
dbe3abbe 196void AliHFEpidTRD::InitParameters(){
197 //
198 // Fill the Parameters into an array
199 //
200
e3fc062d 201 AliDebug(2, "Loading threshold Parameter");
dbe3abbe 202 // Parameters for 6 Layers
203 fThreshParams[0] = -0.001839; // 0.7 electron eff
204 fThreshParams[1] = 0.000276;
205 fThreshParams[2] = 0.044902;
206 fThreshParams[3] = 1.726751;
207 fThreshParams[4] = -0.002405; // 0.75 electron eff
208 fThreshParams[5] = 0.000372;
209 fThreshParams[6] = 0.061775;
210 fThreshParams[7] = 1.739371;
211 fThreshParams[8] = -0.003178; // 0.8 electron eff
212 fThreshParams[9] = 0.000521;
213 fThreshParams[10] = 0.087585;
214 fThreshParams[11] = 1.749154;
215 fThreshParams[12] = -0.004058; // 0.85 electron eff
216 fThreshParams[13] = 0.000748;
217 fThreshParams[14] = 0.129583;
218 fThreshParams[15] = 1.782323;
219 fThreshParams[16] = -0.004967; // 0.9 electron eff
220 fThreshParams[17] = 0.001216;
221 fThreshParams[18] = 0.210128;
222 fThreshParams[19] = 1.807665;
223 fThreshParams[20] = -0.000996; // 0.95 electron eff
224 fThreshParams[21] = 0.002627;
225 fThreshParams[22] = 0.409099;
226 fThreshParams[23] = 1.787076;
809a4336 227}
228
229//___________________________________________________________________
e3fc062d 230void AliHFEpidTRD::InitParameters1DLQ(){
231 //
232 // Init Parameters for 1DLQ PID (M. Fasel, Sept. 6th, 2010)
233 //
234
235 // Parameters for 6 Layers
236 AliDebug(2, Form("Loading threshold parameter for Method 1DLQ"));
237 fThreshParams[0] = -0.02241; // 0.7 electron eff
238 fThreshParams[1] = 0.05043;
239 fThreshParams[2] = 0.7925;
240 fThreshParams[3] = 2.625;
241 fThreshParams[4] = 0.07438; // 0.75 electron eff
242 fThreshParams[5] = 0.05158;
243 fThreshParams[6] = 2.864;
244 fThreshParams[7] = 4.356;
245 fThreshParams[8] = 0.1977; // 0.8 electron eff
246 fThreshParams[9] = 0.05956;
247 fThreshParams[10] = 2.853;
248 fThreshParams[11] = 3.713;
249 fThreshParams[12] = 0.5206; // 0.85 electron eff
250 fThreshParams[13] = 0.03077;
251 fThreshParams[14] = 2.966;
252 fThreshParams[15] = 4.07;
253 fThreshParams[16] = 0.8808; // 0.9 electron eff
254 fThreshParams[17] = 0.002092;
255 fThreshParams[18] = 1.17;
256 fThreshParams[19] = 4.506;
257 fThreshParams[20] = 1.; // 0.95 electron eff
258 fThreshParams[21] = 0.;
259 fThreshParams[22] = 0.;
260 fThreshParams[23] = 0.;
261
262}
263
264//___________________________________________________________________
6555e2ad 265void AliHFEpidTRD::GetParameters(Double_t electronEff, Double_t *parameters) const {
dbe3abbe 266 //
267 // return parameter set for the given efficiency bin
268 //
0792aa82 269 Int_t effbin = static_cast<Int_t>((electronEff - 0.7)/0.05);
dbe3abbe 270 memcpy(parameters, fThreshParams + effbin * 4, sizeof(Double_t) * 4);
809a4336 271}
75d81601 272
273//___________________________________________________________________
6555e2ad 274Double_t AliHFEpidTRD::GetElectronLikelihood(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const {
3a72645a 275 //
276 // Get TRD likelihoods for ESD respectively AOD tracks
277 //
278 Double_t pidProbs[AliPID::kSPECIES]; memset(pidProbs, 0, sizeof(Double_t) * AliPID::kSPECIES);
279 if(anaType == AliHFEpidObject::kESDanalysis){
280 const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
bf892a6a 281 if(esdtrack) esdtrack->GetTRDpid(pidProbs);
3a72645a 282 } else {
283 const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
bf892a6a 284 if(aodtrack)fAODpid->MakeTRDPID(const_cast<AliAODTrack *>(aodtrack), pidProbs);
3a72645a 285 }
286 return pidProbs[AliPID::kElectron];
287}
288
289//___________________________________________________________________
6555e2ad 290Double_t AliHFEpidTRD::GetP(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const {
3a72645a 291 //
292 // Get the Momentum in the TRD
293 //
294 Double_t p = 0.;
295 if(anaType == AliHFEpidObject::kESDanalysis){
296 const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
bf892a6a 297 if(esdtrack) p = esdtrack->GetOuterParam() ? esdtrack->GetOuterParam()->P() : esdtrack->P();
3a72645a 298 } else {
299 const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
bf892a6a 300 if(aodtrack) p = aodtrack->P();
3a72645a 301 }
302 return p;
303}
304
305//___________________________________________________________________
6555e2ad 306Double_t AliHFEpidTRD::GetChargeLayer(const AliVParticle *track, UInt_t layer, AliHFEpidObject::AnalysisType_t anaType) const {
307 //
308 // Get the Charge in a single TRD layer
309 //
310 if(layer >= 6) return 0.;
311 Double_t charge = 0.;
312 if(anaType == AliHFEpidObject::kESDanalysis){
313 const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
bf892a6a 314 if(esdtrack)
315 for(Int_t islice = 0; islice < esdtrack->GetNumberOfTRDslices(); islice++) charge += esdtrack->GetTRDslice(static_cast<UInt_t>(layer), islice);
6555e2ad 316 } else {
317 const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
bf892a6a 318 AliAODPid *aoddetpid = aodtrack ? aodtrack->GetDetPid() : NULL;
319 if(aoddetpid)
320 for(Int_t islice = 0; islice < aoddetpid->GetTRDnSlices(); islice++) charge += aoddetpid->GetTRDsignal()[layer * aoddetpid->GetTRDnSlices() + islice];
6555e2ad 321 }
322 return charge;
323}
324
325//___________________________________________________________________
bf892a6a 326void AliHFEpidTRD::GetTRDmomenta(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType, Double_t *mom) const {
327 //
328 // Fill Array with momentum information at the TRD tracklet
329 //
330 if(anaType == AliHFEpidObject::kESDanalysis){
331 const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
332 if(esdtrack)
333 for(Int_t itl = 0; itl < 6; itl++)
334 mom[itl] = esdtrack->GetTRDmomentum(itl);
335 } else {
336 const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
337 AliAODPid *aoddetpid = aodtrack ? aodtrack->GetDetPid() : NULL;
338 if(aoddetpid){
339 Float_t *trdmom = aoddetpid->GetTRDmomentum();
340 for(Int_t itl = 0; itl < 6; itl++){
341 mom[itl] = trdmom[itl];
342 }
343 }
344 }
345}
346
347//___________________________________________________________________
348Double_t AliHFEpidTRD::GetTRDSignalV1(const AliESDtrack *track, Float_t truncation) const {
75d81601 349 //
350 // Calculation of the TRD Signal via truncated mean
351 // Method 1: Take all Slices available
352 // cut out 0s
353 // Order them in increasing order
354 // Cut out the upper third
355 // Calculate mean over the last 2/3 slices
356 //
357 const Int_t kNSlices = 48;
bf892a6a 358 const Int_t kSlicePerLayer = 7;
359 // Weight the slice to equalize the MPV of the dQ/dl-distribution per slice to the one in the first slice
360 // Pions are used as reference for the equalization
361 const Double_t kWeightSlice[8] = {1., 2.122, 1.8, 1.635, 1.595, 1.614, 1.16, 7.0};
362 AliDebug(3, Form("Number of Tracklets: %d\n", track->GetTRDntrackletsPID()));
75d81601 363 Double_t trdSlices[kNSlices], tmp[kNSlices];
364 Int_t indices[48];
365 Int_t icnt = 0;
366 for(Int_t idet = 0; idet < 6; idet++)
bf892a6a 367 for(Int_t islice = 0; islice < kSlicePerLayer; islice++){
75d81601 368 AliDebug(2, Form("Chamber[%d], Slice[%d]: TRDSlice = %f", idet, islice, track->GetTRDslice(idet, islice)));
369 if(TMath::Abs(track->GetTRDslice(idet, islice)) < fgkVerySmall) continue;;
bf892a6a 370 trdSlices[icnt++] = track->GetTRDslice(idet, islice) * kWeightSlice[islice];
75d81601 371 }
372 AliDebug(1, Form("Number of Slices: %d\n", icnt));
373 if(icnt < 6) return 0.; // We need at least 6 Slices for the truncated mean
374 TMath::Sort(icnt, trdSlices, indices, kFALSE);
375 memcpy(tmp, trdSlices, sizeof(Double_t) * icnt);
376 for(Int_t ien = 0; ien < icnt; ien++)
377 trdSlices[ien] = tmp[indices[ien]];
bf892a6a 378 Double_t trdSignal = TMath::Mean(static_cast<Int_t>(static_cast<Float_t>(icnt) * truncation), trdSlices);
75d81601 379 Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1;
e3fc062d 380 AliDebug(3, Form("PID Meth. 1: p[%f], TRDSignal[%f]", mom, trdSignal));
75d81601 381 return trdSignal;
382}
383
384//___________________________________________________________________
bf892a6a 385Double_t AliHFEpidTRD::GetTRDSignalV2(const AliESDtrack *track, Float_t truncation) const {
75d81601 386 //
387 // Calculation of the TRD Signal via truncated mean
388 // Method 2: Take only first 5 slices per chamber
389 // Order them in increasing order
390 // Cut out upper half
391 // Now do mean with the reamining 3 slices per chamber
392 //
bf892a6a 393 const Double_t kWeightSlice[8] = {1., 2.122, 1.8, 1.635, 1.595, 1.614, 1.16, 7.0};
394 const Int_t kLayers = 6;
395 const Int_t kSlicesLow = 6;
396 const Int_t kSlicesHigh = 1;
397 Double_t trdSlicesLowTime[kLayers*kSlicesLow], trdSlicesRemaining[kLayers*(kSlicesHigh + kSlicesLow)];
398 Int_t indices[kLayers*kSlicesLow];
75d81601 399 Int_t cntLowTime=0, cntRemaining = 0;
400 for(Int_t idet = 0; idet < 6; idet++)
bf892a6a 401 for(Int_t islice = 0; islice < kSlicesLow+kSlicesHigh; islice++){
75d81601 402 if(TMath::Abs(track->GetTRDslice(idet, islice)) < fgkVerySmall) continue;;
bf892a6a 403 if(islice < kSlicesLow){
e3fc062d 404 AliDebug(3, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
bf892a6a 405 trdSlicesLowTime[cntLowTime++] = track->GetTRDslice(idet, islice) * kWeightSlice[islice];
75d81601 406 } else{
e3fc062d 407 AliDebug(3, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
bf892a6a 408 trdSlicesRemaining[cntRemaining++] = track->GetTRDslice(idet, islice) * kWeightSlice[islice];
75d81601 409 }
410 }
411 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
412 TMath::Sort(cntLowTime, trdSlicesLowTime, indices, kFALSE);
413 // Fill the second array with the lower half of the first time bins
bf892a6a 414 for(Int_t ien = 0; ien < static_cast<Int_t>(static_cast<Float_t>(cntLowTime) * truncation); ien++)
75d81601 415 trdSlicesRemaining[cntRemaining++] = trdSlicesLowTime[indices[ien]];
416 Double_t trdSignal = TMath::Mean(cntRemaining, trdSlicesRemaining);
417 Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1;
e3fc062d 418 AliDebug(3, Form("PID Meth. 2: p[%f], TRDSignal[%f]", mom, trdSignal));
75d81601 419 return trdSignal;
420}