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