Since it contains fixes of coding rule violations, all classes are involved. Further...
[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  *                                                                      *
17  * Class for TRD PID                                                    *
18  * Implements the abstract base class AliHFEpidbase                     *
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>
29 #include <TH2F.h>
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"
39 #include "AliLog.h"
40 #include "AliPID.h"
41
42 #include "AliHFEpidTRD.h"
43
44 ClassImp(AliHFEpidTRD)
45
46 const Double_t AliHFEpidTRD::fgkVerySmall = 1e-12;
47
48 //___________________________________________________________________
49 AliHFEpidTRD::AliHFEpidTRD(const char* name) :
50     AliHFEpidBase(name)
51   , fPIDMethod(kNN)
52   , fContainer(0x0)
53 {
54   //
55   // default  constructor
56   // 
57   memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
58 }
59
60 //___________________________________________________________________
61 AliHFEpidTRD::AliHFEpidTRD(const AliHFEpidTRD &ref):
62     AliHFEpidBase("")
63   , fPIDMethod(kLQ)
64   , fContainer(0x0)
65 {
66   //
67   // Copy constructor
68   //
69   memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
70   ref.Copy(*this);
71 }
72
73 //___________________________________________________________________
74 AliHFEpidTRD &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 //___________________________________________________________________
85 void AliHFEpidTRD::Copy(TObject &ref) const {
86   //
87   // Performs the copying of the object
88   //
89   AliHFEpidTRD &target = dynamic_cast<AliHFEpidTRD &>(ref);
90
91   target.fPIDMethod = fPIDMethod;
92   memcpy(target.fThreshParams, fThreshParams, sizeof(Double_t) * kThreshParams);
93   if(fContainer) target.fContainer = dynamic_cast<TList *>(fContainer->Clone());
94   AliHFEpidBase::Copy(ref);
95 }
96
97 //___________________________________________________________________
98 AliHFEpidTRD::~AliHFEpidTRD(){
99   //
100   // Destructor
101   //
102   if(fContainer){
103     fContainer->Clear();
104     delete fContainer;
105   }
106 }
107
108 //______________________________________________________
109 Bool_t AliHFEpidTRD::InitializePID(){
110   //
111   // InitializePID: Load TRD thresholds and create the electron efficiency axis
112   // to navigate 
113   //
114   InitParameters();
115   return kTRUE;
116 }
117
118 //______________________________________________________
119 Int_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
130   Double_t pidProbs[AliPID::kSPECIES];
131   esdTrack->GetTRDpid(pidProbs);
132   if(pidProbs[AliPID::kElectron] > GetTRDthresholds(0.91, p)) return 11;
133   return 0;
134 }
135
136 //___________________________________________________________________
137 Double_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];
144 }
145
146 //___________________________________________________________________
147 void 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;
177 }
178
179 //___________________________________________________________________
180 void 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);
186 }
187
188 //___________________________________________________________________
189 Double_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 //___________________________________________________________________
223 Double_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 //___________________________________________________________________
258 void 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 //___________________________________________________________________
268 void 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 //___________________________________________________________________
278 Int_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 //___________________________________________________________________
296 void 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