New package for heavy flavour electrons analysis (M.Fasel)
[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 <TIterator.h>
30 #include <TKey.h>
31 #include <TMap.h>
32 #include <TObjArray.h>
33 #include <TObjString.h>
34 #include <TString.h>
35 #include <TROOT.h>
36
37 #include "AliESDtrack.h"
38 #include "AliPID.h"
39
40 #include "AliHFEpidTRD.h"
41
42 ClassImp(AliHFEpidTRD)
43
44 //___________________________________________________________________
45 AliHFEpidTRD::AliHFEpidTRD(const char* name) :
46     AliHFEpidBase(name)
47   , fThresholdFile("TRD.PIDthresholds.root")
48   , fPIDMethod(kNN)
49   , fTRDthresholds(0x0)
50   , fTRDelectronEfficiencies(0x0)
51 {
52   //
53   // default  constructor
54   // 
55   fTRDthresholds = new TMap();
56 }
57
58 //___________________________________________________________________
59 AliHFEpidTRD::AliHFEpidTRD(const AliHFEpidTRD &ref):
60     AliHFEpidBase("")
61   , fThresholdFile("")
62   , fPIDMethod(kLQ)
63   , fTRDthresholds(0x0)
64   , fTRDelectronEfficiencies(0x0)
65 {
66   //
67   // Copy constructor
68   //
69   ref.Copy(*this);
70 }
71
72 //___________________________________________________________________
73 AliHFEpidTRD &AliHFEpidTRD::operator=(const AliHFEpidTRD &ref){
74   //
75   // Assignment operator
76   //
77   if(this != &ref){
78     ref.Copy(*this);
79   }
80   return *this;
81 }
82
83 //___________________________________________________________________
84 void AliHFEpidTRD::Copy(TObject &ref) const {
85   //
86   // Performs the copying of the object
87   //
88   AliHFEpidTRD &target = dynamic_cast<AliHFEpidTRD &>(ref);
89
90   target.fThresholdFile = fThresholdFile;
91   target.fPIDMethod = fPIDMethod;
92   if(fTRDthresholds){
93     target.fTRDthresholds = dynamic_cast<TMap *>(fTRDthresholds->Clone());
94   }
95   if(fTRDelectronEfficiencies){
96     target.fTRDelectronEfficiencies = new TAxis(*fTRDelectronEfficiencies);
97   }
98   AliHFEpidBase::Copy(ref);
99 }
100
101 //___________________________________________________________________
102 AliHFEpidTRD::~AliHFEpidTRD(){
103   //
104   // Destructor
105   //
106   if(fTRDthresholds){
107     fTRDthresholds->Delete();
108     delete fTRDthresholds;
109   }
110   if(fTRDelectronEfficiencies) delete fTRDelectronEfficiencies;
111 }
112
113 //______________________________________________________
114 Bool_t AliHFEpidTRD::InitializePID(){
115   //
116   // InitializePID: Load TRD thresholds and create the electron efficiency axis
117   // to navigate 
118   //
119   LoadTRDthresholds();
120   // Fill the electron efficiencies axis for the TRD alone PID
121   Int_t nEffs = fTRDthresholds->GetEntries() + 1;
122   TIterator* thres =fTRDthresholds->MakeIterator();
123   TObjString *key = 0x0;
124   Double_t *tmp = new Double_t[nEffs], *titer = tmp;
125   while((key = dynamic_cast<TObjString *>((*thres)()))){
126     (*titer++) = static_cast<Double_t>(key->String().Atoi())/100.;
127   }
128   delete thres;
129   *titer = 1.;
130   // Sort the electron efficiencies and put them into the TAxis for later navigation
131   Int_t *ind = new Int_t[nEffs], *iiter = ind;
132   TMath::Sort(nEffs, tmp, ind, kFALSE);
133   Double_t *eleffs = new Double_t[nEffs], *eiter = eleffs;
134   while(eiter < &eleffs[nEffs]) *(eiter++) = tmp[*(iiter++)];
135   // print the content
136   Int_t cnt = 0;
137   if(GetDebugLevel() > 1){
138     printf("Printing electron efficiency bins:\n");
139     eiter = eleffs;
140     thres = fTRDthresholds->MakeIterator();
141     TObject *object = 0x0;
142     while(eiter < &eleffs[nEffs - 1]){
143       printf("eleffs[%d] = %f", cnt++, *(eiter++));
144       key = dynamic_cast<TObjString *>(thres->Next());
145       object = fTRDthresholds->GetValue(key->String().Data());
146       printf(", Content: %p\n", (void *)object);
147     }
148   }
149   delete[] tmp; delete[] ind;
150   fTRDelectronEfficiencies = new TAxis(nEffs - 1, eleffs);
151   if(GetDebugLevel() > 1){
152     printf("Printing axis content:\n");
153     for(Int_t ibin = fTRDelectronEfficiencies->GetFirst(); ibin <= fTRDelectronEfficiencies->GetLast(); ibin++)
154       printf("%d.) minimum: %f, maximum? %f\n", ibin, fTRDelectronEfficiencies->GetBinLowEdge(ibin), fTRDelectronEfficiencies->GetBinUpEdge(ibin));
155   }
156   delete[] eleffs;
157   return kTRUE;
158 }
159
160 //______________________________________________________
161 Int_t AliHFEpidTRD::IsSelected(AliVParticle *track){
162   //
163   // Does PID for TRD alone:
164   // PID thresholds based on 90% Electron Efficiency level approximated by a linear 
165   // step function
166   //
167   AliESDtrack *esdTrack = 0x0;
168   if(!(esdTrack = dynamic_cast<AliESDtrack *>(track))) return kFALSE;
169   Double_t p = esdTrack->GetOuterParam() ? esdTrack->GetOuterParam()->P() : esdTrack->P();
170   if(p < 0.6) return 0;
171
172   // Get the Histograms
173   TH1 *threshist = GetTRDthresholds(0.91);
174   Int_t bin = 0;
175   if(p > threshist->GetXaxis()->GetXmax()) 
176     bin = threshist->GetXaxis()->GetLast();
177   else if(p < threshist->GetXaxis()->GetXmin()) 
178     bin = threshist->GetXaxis()->GetFirst();
179   else
180     bin = threshist->GetXaxis()->FindBin(p);
181
182   Double_t pidProbs[AliPID::kSPECIES];
183   esdTrack->GetTRDpid(pidProbs);
184   if(pidProbs[AliPID::kElectron] > threshist->GetBinContent(bin)) return 11;
185   return 0;
186 }
187
188 //___________________________________________________________________
189 void AliHFEpidTRD::LoadTRDthresholds(){
190   //
191   // Load TRD threshold histograms from File
192   //
193   TFile *mythresholds = TFile::Open(fThresholdFile);
194   TKey *object = 0x0;
195   TString electron_eff;
196   TObjArray *histos = 0x0;
197   Float_t eff;
198   TH1F *refhist = 0x0;
199   TIterator *keyIterator = mythresholds->GetListOfKeys()->MakeIterator();
200   TString histnames[2] = {"fHistThreshLQ", "fHistThreshNN"};
201   gROOT->cd();
202   while((object = dynamic_cast<TKey *>((*keyIterator)()))){
203     // Get the electron efficiency bin this histogram was taken with
204     electron_eff = object->GetName();
205     electron_eff = electron_eff.Remove(0,3);
206     eff = static_cast<Float_t>(electron_eff.Atoi())/100.;
207
208     // Get the threshold according to the selected 
209     histos = dynamic_cast<TObjArray *>(object->ReadObj());
210     refhist = dynamic_cast<TH1F *>(histos->FindObject(histnames[fPIDMethod].Data()));
211     SetTRDthresholds(refhist, eff);
212     histos->Delete();
213     delete histos;
214   }
215   delete keyIterator;
216   mythresholds->Close();
217   delete mythresholds;
218 }
219
220 //___________________________________________________________________
221 void AliHFEpidTRD::SetTRDthresholds(TH1F *thresholds, Float_t electronEff){
222   //
223   // Set the threshold histogram for the TRD pid
224   //
225   fTRDthresholds->Add(new TObjString(Form("%d", TMath::Nint(electronEff * 100.))), new TH1F(*thresholds));
226 }
227
228 //___________________________________________________________________
229 TH1F *AliHFEpidTRD::GetTRDthresholds(Float_t electronEff){ 
230   Int_t bin = 0;
231   if(electronEff < fTRDelectronEfficiencies->GetXmin()) bin = fTRDelectronEfficiencies->GetFirst();
232   else if(electronEff > fTRDelectronEfficiencies->GetXmax()) bin = fTRDelectronEfficiencies->GetLast();
233   else bin = fTRDelectronEfficiencies->FindBin(electronEff);
234  TObjString keyname = Form("%d", TMath::Nint(fTRDelectronEfficiencies->GetBinLowEdge(bin)* 100.));
235 /*  printf("Key: %s\n", keyname.String().Data());*/
236   TH1F *thresholds = dynamic_cast<TH1F *>((dynamic_cast<TPair *>(fTRDthresholds->FindObject(&keyname)))->Value());
237 /*  printf("thresholds: %p\n", thresholds);*/
238   return thresholds;
239 }