1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 // Implements the abstract base class AliHFEpidbase
18 // Make PID does the PID decision
19 // Class further contains TRD specific cuts and QA histograms
22 // Markus Fasel <M.Fasel@gsi.de>
26 #include <THnSparse.h>
30 #include "AliAODPid.h"
31 #include "AliAODTrack.h"
32 #include "AliAODMCParticle.h"
33 #include "AliESDtrack.h"
35 #include "AliMCParticle.h"
36 #include "AliOADBContainer.h"
38 #include "AliPIDResponse.h"
40 #include "AliHFEOADBThresholdsTRD.h"
41 #include "AliHFEpidQAmanager.h"
42 #include "AliHFEpidTRD.h"
44 ClassImp(AliHFEpidTRD)
46 //___________________________________________________________________
47 AliHFEpidTRD::AliHFEpidTRD() :
49 , fOADBThresholds(NULL)
53 , fElectronEfficiency(0.90)
54 , fTotalChargeInSlice0(kFALSE)
57 // default constructor
59 memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
62 //___________________________________________________________________
63 AliHFEpidTRD::AliHFEpidTRD(const char* name) :
65 , fOADBThresholds(NULL)
69 , fElectronEfficiency(0.91)
70 , fTotalChargeInSlice0(kFALSE)
73 // default constructor
75 memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
78 //___________________________________________________________________
79 AliHFEpidTRD::AliHFEpidTRD(const AliHFEpidTRD &ref):
81 , fOADBThresholds(NULL)
83 , fNTracklets(ref.fNTracklets)
84 , fRunNumber(ref.fRunNumber)
85 , fElectronEfficiency(ref.fElectronEfficiency)
86 , fTotalChargeInSlice0(ref.fTotalChargeInSlice0)
91 memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
95 //___________________________________________________________________
96 AliHFEpidTRD &AliHFEpidTRD::operator=(const AliHFEpidTRD &ref){
98 // Assignment operator
106 //___________________________________________________________________
107 void AliHFEpidTRD::Copy(TObject &ref) const {
109 // Performs the copying of the object
111 AliHFEpidTRD &target = dynamic_cast<AliHFEpidTRD &>(ref);
113 target.fMinP = fMinP;
114 target.fNTracklets = fNTracklets;
115 target.fRunNumber = fRunNumber;
116 target.fTotalChargeInSlice0 = fTotalChargeInSlice0;
117 target.fElectronEfficiency = fElectronEfficiency;
118 memcpy(target.fThreshParams, fThreshParams, sizeof(Double_t) * kThreshParams);
119 AliHFEpidBase::Copy(ref);
122 //___________________________________________________________________
123 AliHFEpidTRD::~AliHFEpidTRD(){
129 //______________________________________________________
130 Bool_t AliHFEpidTRD::InitializePID(Int_t run){
132 // InitializePID: Load TRD thresholds and create the electron efficiency axis
135 AliDebug(1, Form("Initializing TRD PID for run %d", run));
136 if(InitParamsFromOADB(run)){
137 SetBit(kThresholdsInitialized);
140 AliDebug(1, Form("Threshold Parameters for %d tracklets and an electron efficiency %f loaded:", fNTracklets, fElectronEfficiency));
141 AliDebug(1, Form("Params: [%f|%f|%f|%f]", fThreshParams[0], fThreshParams[1], fThreshParams[2], fThreshParams[3]));
146 //______________________________________________________
147 Int_t AliHFEpidTRD::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager *pidqa) const {
149 // Does PID for TRD alone:
150 // PID thresholds based on 90% Electron Efficiency level approximated by a linear
153 if(!TestBit(kThresholdsInitialized)) {
154 AliDebug(1,"Threshold Parameters not available");
157 AliDebug(2, "Applying TRD PID");
159 AliDebug(2, "Cannot process track");
164 const AliESDtrack *esdt = dynamic_cast<const AliESDtrack *>(track->GetRecTrack());
165 printf("checking IdentifiedAsElectronTRD, number of Tracklets: %d\n", esdt->GetTRDntrackletsPID());
166 if(fkPIDResponse->IdentifiedAsElectronTRD(dynamic_cast<const AliVTrack *>(track->GetRecTrack()), 0.8)) printf("Track identified as electron\n");
167 else printf("Track rejected\n");
169 AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis: AliHFEpidObject::kAODanalysis;
170 Double_t p = GetP(track->GetRecTrack(), anatype);
172 AliDebug(2, Form("Track momentum below %f", fMinP));
176 if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kTRDpid, AliHFEdetPIDqa::kBeforePID);
177 Double_t electronLike = GetElectronLikelihood(static_cast<const AliVTrack *>(track->GetRecTrack()), anatype);
179 if(TestBit(kSelectCutOnTheFly)){
180 threshold = GetTRDthresholds(p, (static_cast<const AliVTrack *>(track->GetRecTrack()))->GetTRDntrackletsPID());
182 threshold = GetTRDthresholds(p);
184 AliDebug(2, Form("Threshold: %f\n", threshold));
185 if(electronLike > threshold){
186 if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kTRDpid, AliHFEdetPIDqa::kAfterPID);
193 //___________________________________________________________________
194 Double_t AliHFEpidTRD::GetTRDthresholds(Double_t p, UInt_t nTracklets) const {
196 // Return momentum dependent and electron efficiency dependent TRD thresholds
197 // Determine threshold based on the number of tracklets on the fly, electron efficiency not modified
199 Double_t threshParams[4];
200 // Get threshold paramters for the given number of tracklets from OADB container
201 AliHFEOADBThresholdsTRD *thresholds = dynamic_cast<AliHFEOADBThresholdsTRD *>(fOADBThresholds->GetObject(fRunNumber));
203 AliDebug(1, Form("Thresholds for run %d not in the OADB", fRunNumber));
206 if(!thresholds->GetThresholdParameters(nTracklets, fElectronEfficiency, threshParams)) return 0.;
208 Double_t threshold = 1. - threshParams[0] - threshParams[1] * p - threshParams[2] * TMath::Exp(-threshParams[3] * p);
209 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
212 //___________________________________________________________________
213 Double_t AliHFEpidTRD::GetTRDthresholds(Double_t p) const {
215 // Return momentum dependent and electron efficiency dependent TRD thresholds
217 Double_t threshold = 1. - fThreshParams[0] - fThreshParams[1] * p - fThreshParams[2] * TMath::Exp(-fThreshParams[3] * p);
218 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
222 //___________________________________________________________________
223 Bool_t AliHFEpidTRD::InitParamsFromOADB(Int_t run){
225 // The name of the function says it all
227 AliHFEOADBThresholdsTRD *thresholds = dynamic_cast<AliHFEOADBThresholdsTRD *>(fOADBThresholds->GetObject(run));
229 AliDebug(1, Form("Thresholds for run %d not in the OADB", run));
232 return thresholds->GetThresholdParameters(fNTracklets, fElectronEfficiency, fThreshParams);
235 //___________________________________________________________________
236 void AliHFEpidTRD::RenormalizeElPi(const Double_t * const likein, Double_t * const likeout) const {
238 // Renormalize likelihoods for electrons and pions neglecting the
239 // likelihoods for protons, kaons and muons
241 memset(likeout, 0, sizeof(Double_t) * AliPID::kSPECIES);
242 Double_t norm = likein[AliPID::kElectron] + likein[AliPID::kPion];
243 if(norm == 0.) norm = 1.; // Safety
244 likeout[AliPID::kElectron] = likein[AliPID::kElectron] / norm;
245 likeout[AliPID::kPion] = likein[AliPID::kPion] / norm;
248 //___________________________________________________________________
249 Double_t AliHFEpidTRD::GetElectronLikelihood(const AliVTrack *track, AliHFEpidObject::AnalysisType_t anaType) const {
251 // Get TRD likelihoods for ESD respectively AOD tracks
253 Double_t pidProbs[AliPID::kSPECIES]; memset(pidProbs, 0, sizeof(Double_t) * AliPID::kSPECIES);
254 if(anaType == AliHFEpidObject::kESDanalysis){
255 const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
256 if(esdtrack) esdtrack->GetTRDpid(pidProbs);
258 fkPIDResponse->ComputeTRDProbability(track, AliPID::kSPECIES, pidProbs);
260 if(!IsRenormalizeElPi()) return pidProbs[AliPID::kElectron];
261 Double_t probsNew[AliPID::kSPECIES];
262 RenormalizeElPi(pidProbs, probsNew);
263 return probsNew[AliPID::kElectron];
266 //___________________________________________________________________
267 Double_t AliHFEpidTRD::GetP(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const {
269 // Get the Momentum in the TRD
272 if(anaType == AliHFEpidObject::kESDanalysis){
273 const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
274 if(esdtrack) p = esdtrack->GetOuterParam() ? esdtrack->GetOuterParam()->P() : esdtrack->P();
276 const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
277 if(aodtrack) p = aodtrack->P();
282 //___________________________________________________________________
283 Double_t AliHFEpidTRD::GetChargeLayer(const AliVParticle *track, UInt_t layer, AliHFEpidObject::AnalysisType_t anaType) const {
285 // Get the Charge in a single TRD layer
287 if(layer >= 6) return 0.;
288 Double_t charge = 0.;
289 if(anaType == AliHFEpidObject::kESDanalysis){
290 const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
292 // Distinction between old and new reconstruction: in the new reconstruction, the total charge is stored in slice 0, slices 1 to 8 are used for the slices for
293 // the neural network.
294 if(fTotalChargeInSlice0)
295 charge = esdtrack->GetTRDslice(static_cast<UInt_t>(layer), 0);
297 for(Int_t islice = 0; islice < esdtrack->GetNumberOfTRDslices(); islice++) charge += esdtrack->GetTRDslice(static_cast<UInt_t>(layer), islice);
300 const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
301 AliAODPid *aoddetpid = aodtrack ? aodtrack->GetDetPid() : NULL;
303 if(fTotalChargeInSlice0)
304 charge = aoddetpid->GetTRDsignal()[layer * aoddetpid->GetTRDnSlices()];
306 for(Int_t islice = 0; islice < aoddetpid->GetTRDnSlices(); islice++) charge += aoddetpid->GetTRDsignal()[layer * aoddetpid->GetTRDnSlices() + islice];
312 //___________________________________________________________________
313 void AliHFEpidTRD::GetTRDmomenta(const AliVTrack *track, Double_t *mom) const {
315 // Fill Array with momentum information at the TRD tracklet
317 for(Int_t itl = 0; itl < 6; itl++)
318 mom[itl] = track->GetTRDmomentum(itl);
321 //___________________________________________________________________
322 Double_t AliHFEpidTRD::GetTRDSignalV1(const AliESDtrack *track, Float_t truncation) const {
324 // Calculation of the TRD Signal via truncated mean
325 // Method 1: Take all Slices available
327 // Order them in increasing order
328 // Cut out the upper third
329 // Calculate mean over the last 2/3 slices
331 const Int_t kNSlices = 48;
332 const Int_t kSlicePerLayer = 7;
333 // Weight the slice to equalize the MPV of the dQ/dl-distribution per slice to the one in the first slice
334 // Pions are used as reference for the equalization
335 const Double_t kWeightSlice[8] = {1., 2.122, 1.8, 1.635, 1.595, 1.614, 1.16, 7.0};
336 AliDebug(3, Form("Number of Tracklets: %d\n", track->GetTRDntrackletsPID()));
337 Double_t trdSlices[kNSlices], tmp[kNSlices];
340 for(Int_t idet = 0; idet < 6; idet++)
341 for(Int_t islice = fTotalChargeInSlice0 ? 1 : 0 ; islice < kSlicePerLayer; islice++){
342 AliDebug(2, Form("Chamber[%d], Slice[%d]: TRDSlice = %f", idet, islice, track->GetTRDslice(idet, islice)));
343 if(TMath::Abs(track->GetTRDslice(idet, islice)) < fgkVerySmall) continue;;
344 trdSlices[icnt++] = track->GetTRDslice(idet, islice) * kWeightSlice[islice];
346 AliDebug(1, Form("Number of Slices: %d\n", icnt));
347 if(icnt < 6) return 0.; // We need at least 6 Slices for the truncated mean
348 TMath::Sort(icnt, trdSlices, indices, kFALSE);
349 memcpy(tmp, trdSlices, sizeof(Double_t) * icnt);
350 for(Int_t ien = 0; ien < icnt; ien++)
351 trdSlices[ien] = tmp[indices[ien]];
352 Double_t trdSignal = TMath::Mean(static_cast<Int_t>(static_cast<Float_t>(icnt) * truncation), trdSlices);
353 Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1;
354 AliDebug(3, Form("PID Meth. 1: p[%f], TRDSignal[%f]", mom, trdSignal));
358 //___________________________________________________________________
359 Double_t AliHFEpidTRD::GetTRDSignalV2(const AliESDtrack *track, Float_t truncation) const {
361 // Calculation of the TRD Signal via truncated mean
362 // Method 2: Take only first 5 slices per chamber
363 // Order them in increasing order
364 // Cut out upper half
365 // Now do mean with the reamining 3 slices per chamber
367 const Double_t kWeightSlice[8] = {1., 2.122, 1.8, 1.635, 1.595, 1.614, 1.16, 7.0};
368 const Int_t kLayers = 6;
369 const Int_t kSlicesLow = 6;
370 const Int_t kSlicesHigh = 1;
371 Double_t trdSlicesLowTime[kLayers*kSlicesLow], trdSlicesRemaining[kLayers*(kSlicesHigh + kSlicesLow)];
372 Int_t indices[kLayers*kSlicesLow];
373 Int_t cntLowTime=0, cntRemaining = 0;
374 for(Int_t idet = 0; idet < 6; idet++)
375 for(Int_t islice = fTotalChargeInSlice0 ? 1 : 0; islice < kSlicesLow+kSlicesHigh; islice++){
376 if(TMath::Abs(track->GetTRDslice(idet, islice)) < fgkVerySmall) continue;;
377 if(islice < kSlicesLow){
378 AliDebug(3, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
379 trdSlicesLowTime[cntLowTime++] = track->GetTRDslice(idet, islice) * kWeightSlice[islice];
381 AliDebug(3, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
382 trdSlicesRemaining[cntRemaining++] = track->GetTRDslice(idet, islice) * kWeightSlice[islice];
385 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
386 TMath::Sort(cntLowTime, trdSlicesLowTime, indices, kFALSE);
387 // Fill the second array with the lower half of the first time bins
388 for(Int_t ien = 0; ien < static_cast<Int_t>(static_cast<Float_t>(cntLowTime) * truncation); ien++)
389 trdSlicesRemaining[cntRemaining++] = trdSlicesLowTime[indices[ien]];
390 Double_t trdSignal = TMath::Mean(cntRemaining, trdSlicesRemaining);
391 Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1;
392 AliDebug(3, Form("PID Meth. 2: p[%f], TRDSignal[%f]", mom, trdSignal));