Since it contains fixes of coding rule violations, all classes are involved. Further...
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEpid.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  * PID Steering Class                                                   *
18  * Interface to the user task                                           *
19  *                                                                      *
20  * Authors:                                                             *
21  *   Markus Fasel <M.Fasel@gsi.de>                                      *
22  *                                                                      *
23  ************************************************************************/
24 #include <TClass.h>
25 #include <THnSparse.h>
26 #include <TH3F.h>
27 #include <TIterator.h>
28 #include <TList.h>
29 #include <TObjArray.h>
30 #include <TObjString.h>
31 #include <TString.h>
32
33 #include "AliESDtrack.h"
34 #include "AliLog.h"
35 #include "AliPID.h"
36
37 #include "AliHFEpid.h"
38 #include "AliHFEpidBase.h"
39 #include "AliHFEpidITS.h"
40 #include "AliHFEpidTPC.h"
41 #include "AliHFEpidTRD.h"
42 #include "AliHFEpidTOF.h"
43 #include "AliHFEpidMC.h"
44
45 ClassImp(AliHFEpid)
46
47 //____________________________________________________________
48 AliHFEpid::AliHFEpid():
49   fEnabledDetectors(0),
50   fQAlist(0x0)
51 {
52   //
53   // Default constructor
54   //
55   memset(fDetectorPID, 0, sizeof(AliHFEpidBase *) * kNdetectorPID);
56 }
57
58 //____________________________________________________________
59 AliHFEpid::AliHFEpid(const AliHFEpid &c):
60   TObject(c),
61   fEnabledDetectors(c.fEnabledDetectors),
62   fQAlist(0x0)
63 {
64   //
65   // Copy Constructor
66   //
67   memset(fDetectorPID, 0, sizeof(AliHFEpidBase *) * kNdetectorPID);
68   if(c.fDetectorPID[kMCpid])
69     fDetectorPID[kMCpid] = new AliHFEpidMC(*(dynamic_cast<AliHFEpidMC *>(c.fDetectorPID[kMCpid])));
70   if(c.fDetectorPID[kTPCpid])
71     fDetectorPID[kTPCpid] = new AliHFEpidTPC(*(dynamic_cast<AliHFEpidTPC *>(c.fDetectorPID[kTPCpid])));
72   if(c.fDetectorPID[kTRDpid])
73     fDetectorPID[kTRDpid] = new AliHFEpidTRD(*(dynamic_cast<AliHFEpidTRD *>(c.fDetectorPID[kTOFpid])));
74   if(c.fDetectorPID[kTOFpid])
75     fDetectorPID[kTOFpid] = new AliHFEpidTOF(*(dynamic_cast<AliHFEpidTOF *>(c.fDetectorPID[kTOFpid])));
76   if(c.IsQAOn()) SetQAOn();
77   if(c.HasMCData()) SetHasMCData(kTRUE);
78   for(Int_t idet = 0; idet < kNdetectorPID; idet++){
79     if(c.IsQAOn() && fDetectorPID[idet]) fDetectorPID[idet]->SetQAOn(fQAlist);
80     if(c.HasMCData() && fDetectorPID[idet]) fDetectorPID[idet]->SetHasMCData(kTRUE);
81   }
82 }
83
84 //____________________________________________________________
85 AliHFEpid& AliHFEpid::operator=(const AliHFEpid &c){
86   //
87   // Assignment operator
88   //
89   TObject::operator=(c);
90
91   if(this != &c){
92     fEnabledDetectors = c.fEnabledDetectors;
93     fQAlist = 0x0;
94   
95     memset(fDetectorPID, 0, sizeof(AliHFEpidBase *) * kNdetectorPID);
96     if(c.fDetectorPID[kMCpid])
97       fDetectorPID[kMCpid] = new AliHFEpidMC(*(dynamic_cast<AliHFEpidMC *>(c.fDetectorPID[kMCpid])));
98     if(c.fDetectorPID[kTPCpid])
99       fDetectorPID[kTPCpid] = new AliHFEpidTPC(*(dynamic_cast<AliHFEpidTPC *>(c.fDetectorPID[kTPCpid])));
100     if(c.fDetectorPID[kTRDpid])
101       fDetectorPID[kTRDpid] = new AliHFEpidTRD(*(dynamic_cast<AliHFEpidTRD *>(c.fDetectorPID[kTOFpid])));
102     if(c.fDetectorPID[kTOFpid])
103       fDetectorPID[kTOFpid] = new AliHFEpidTOF(*(dynamic_cast<AliHFEpidTOF *>(c.fDetectorPID[kTOFpid])));
104     if(c.IsQAOn()) SetQAOn();
105     if(c.HasMCData()) SetHasMCData(kTRUE);
106     for(Int_t idet = 0; idet < kNdetectorPID; idet++){
107       if(c.IsQAOn() && fDetectorPID[idet]) fDetectorPID[idet]->SetQAOn(fQAlist);
108       if(c.HasMCData() && fDetectorPID[idet]) fDetectorPID[idet]->SetHasMCData();
109     }
110   }
111   return *this; 
112 }
113
114 //____________________________________________________________
115 AliHFEpid::~AliHFEpid(){
116   //
117   // Destructor
118   //
119   for(Int_t idet = 0; idet < kNdetectorPID; idet++){
120     if(fDetectorPID[idet])
121       delete fDetectorPID[idet];
122   } 
123   if(fQAlist) delete fQAlist; fQAlist = 0x0;  // Each detector has to care about its Histograms
124 }
125
126 //____________________________________________________________
127 Bool_t AliHFEpid::InitializePID(TString detectors){
128   //
129   // Initializes PID Object:
130   // + Defines which detectors to use
131   // + Initializes Detector PID objects
132   // + Handles QA
133   //
134   AliInfo(Form("Doing InitializePID for Detectors %s", detectors.Data()));
135   fDetectorPID[kMCpid] = new AliHFEpidMC("Monte Carlo PID"); // Always there
136   fDetectorPID[kITSpid] = new AliHFEpidITS("ITS development PID");  // Development version of the ITS pid, for the moment always there
137   SETBIT(fEnabledDetectors, kMCpid);
138   SETBIT(fEnabledDetectors, kITSpid);
139   
140   TObjArray *detsEnabled = detectors.Tokenize(":");
141   TIterator *detIterator = detsEnabled->MakeIterator();
142   TObjString *det = 0x0;
143   while((det = dynamic_cast<TObjString *>(detIterator->Next()))){
144     if(det->String().CompareTo("TPC") == 0){
145       fDetectorPID[kTPCpid] = new AliHFEpidTPC("TPC PID");
146       SETBIT(fEnabledDetectors, kTPCpid);
147     } else if(det->String().CompareTo("TRD") == 0){
148       fDetectorPID[kTRDpid] = new AliHFEpidTRD("TRD PID");
149       SETBIT(fEnabledDetectors, kTRDpid);
150     } else if(det->String().CompareTo("TOF") == 0){
151       fDetectorPID[kTOFpid] = new AliHFEpidTOF("TOF PID");
152       SETBIT(fEnabledDetectors, kTOFpid);
153     }
154     // Here is still space for ESD PID
155   }
156   // Initialize PID Objects
157   Bool_t status = kTRUE;
158   for(Int_t idet = 0; idet < kNdetectorPID; idet++){
159     if(fDetectorPID[idet]){ 
160       status &= fDetectorPID[idet]->InitializePID();
161       if(IsQAOn() && status) fDetectorPID[idet]->SetQAOn(fQAlist);
162       if(HasMCData() && status) fDetectorPID[idet]->SetHasMCData();
163     }
164   }
165   return status;
166 }
167
168 //____________________________________________________________
169 Bool_t AliHFEpid::IsSelected(AliVParticle *track){
170   //
171   // Steers PID decision for single detectors respectively combined
172   // PID decision
173   //
174   if(TString(track->IsA()->GetName()).CompareTo("AliMCparticle") == 0){
175     return (TMath::Abs(fDetectorPID[kMCpid]->IsSelected(track)) == 11);
176   }
177   if(TString(track->IsA()->GetName()).CompareTo("AliESDtrack") == 0){
178     if(TESTBIT(fEnabledDetectors, kTPCpid)){
179       if(IsQAOn()) 
180         MakePlotsItsTpc(dynamic_cast<AliESDtrack *>(track));  // First fill the QA histograms
181       if(TESTBIT(fEnabledDetectors, kTOFpid)){
182         // case TPC-TOF
183         return MakePidTpcTof(dynamic_cast<AliESDtrack *>(track));
184       } else if(TESTBIT(fEnabledDetectors, kTRDpid)){
185         // case TPC-TRD with low level detector Signals
186         return MakePidTpcTrd(dynamic_cast<AliESDtrack *>(track));
187       } else
188         return (TMath::Abs(fDetectorPID[kTPCpid]->IsSelected(track)) ==11);
189     } else if(TESTBIT(fEnabledDetectors, kTRDpid)){
190       return (TMath::Abs(fDetectorPID[kTRDpid]->IsSelected(track)) ==11);
191     } else if(TESTBIT(fEnabledDetectors, kTOFpid)){
192       return (TMath::Abs(fDetectorPID[kTOFpid]->IsSelected(track)) ==11);
193     }
194     
195   }
196   return kFALSE;
197 }
198
199 //____________________________________________________________
200 Bool_t AliHFEpid::MakePidTpcTof(AliESDtrack *track){
201   //
202   // Combines TPC and TOF PID decision
203   //
204   if(fDetectorPID[kTOFpid]->IsSelected(track)) return fDetectorPID[kTPCpid]->IsSelected(track);
205   return kFALSE;
206 }
207
208 //____________________________________________________________
209 Bool_t AliHFEpid::MakePidTpcTrd(AliESDtrack *track){
210   //
211   // Combination of TPC and TRD PID
212   // Fills Histograms TPC Signal vs. TRD signal for different
213   // momentum slices
214   //
215   Double_t content[10];
216   content[0] = -1;
217   content[1] = track->P();
218   content[2] = track->GetTPCsignal();
219   AliHFEpidTRD *trdPid = dynamic_cast<AliHFEpidTRD *>(fDetectorPID[kTRDpid]);
220   content[3] = trdPid->GetTRDSignalV1(track);
221   content[4] = trdPid->GetTRDSignalV2(track);
222   AliDebug(1, Form("Momentum: %f, TRD Signal: Method 1[%f], Method 2[%f]", content[1], content[3], content[4]));
223   if(IsQAOn()){
224     if(HasMCData()){
225       // Fill My Histograms for MC PID
226       Int_t pdg = TMath::Abs(fDetectorPID[kMCpid]->IsSelected(track));
227       Int_t pid = -1;
228       switch(pdg){
229         case 11:    pid = AliPID::kElectron; break;
230         case 13:    pid = AliPID::kMuon; break;
231         case 211:   pid = AliPID::kPion; break;
232         case 321:   pid = AliPID::kKaon; break;
233         case 2212:  pid = AliPID::kProton; break;
234         default:    pid = -1;
235       };
236       content[0] = pid;
237     }
238     (dynamic_cast<THnSparse *>(fQAlist->At(kTRDSignal)))->Fill(content);
239   }
240   return trdPid->IsSelected(track);
241 }
242
243 //____________________________________________________________
244 void AliHFEpid::MakePlotsItsTpc(AliESDtrack *track){
245   //
246   // Make a plot ITS signal - TPC signal for several momentum bins
247   //
248   Double_t content[10];
249   content[0] = -1;
250   content[1] = track->GetTPCInnerParam() ? track->GetTPCInnerParam()->P() : track->P();
251   content[2] = (dynamic_cast<AliHFEpidITS *>(fDetectorPID[kITSpid]))->GetITSSignalV1(track);
252   content[3] = track->GetTPCsignal();
253   AliDebug(1, Form("Momentum %f, TPC Signal %f, ITS Signal %f", content[1], content[2], content[3]));
254   if(HasMCData()){
255     // Fill My Histograms for MC PID
256     Int_t pdg = TMath::Abs(fDetectorPID[kMCpid]->IsSelected(track));
257     Int_t pid = -1;
258     switch(pdg){
259       case 11:    pid = AliPID::kElectron; break;
260       case 13:    pid = AliPID::kMuon; break;
261       case 211:   pid = AliPID::kPion; break;
262       case 321:   pid = AliPID::kKaon; break;
263       case 2212:  pid = AliPID::kProton; break;
264       default:    pid = -1;
265     };
266     content[0] = pid;
267   }
268   (dynamic_cast<THnSparse *>(fQAlist->At(kITSSignal)))->Fill(content);
269 }
270
271 //____________________________________________________________
272 void AliHFEpid::SetQAOn(){
273   //
274   // Switch on QA
275   //
276   SetBit(kIsQAOn, kTRUE);
277   if(fQAlist) return;
278   fQAlist = new TList;
279   fQAlist->SetName("PIDqa");
280   THnSparseF *histo = NULL;
281
282   // Prepare axis for QA histograms
283   const Int_t kMomentumBins = 41;
284   const Double_t kPtMin = 0.1;
285   const Double_t kPtMax = 10.;
286   Double_t momentumBins[kMomentumBins];
287   for(Int_t ibin = 0; ibin < kMomentumBins; ibin++)
288     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)));
289
290   // Add Histogram for combined TPC-TRD PID
291   const Int_t kDimensionsTRDsig = 5;
292   Int_t kNbinsTRDsig[kDimensionsTRDsig] = {AliPID::kSPECIES + 1, kMomentumBins - 1, 200, 3000, 3000};
293   Double_t binMinTRDsig[kDimensionsTRDsig] = {-1., 0.1, 0, 0, 0};
294   Double_t binMaxTRDsig[kDimensionsTRDsig] = {AliPID::kSPECIES, 10., 200., 3000., 3000.};
295   fQAlist->AddAt((histo = new THnSparseF("fCombTPCTRDpid", "Combined TPC-TRD PID", kDimensionsTRDsig, kNbinsTRDsig, binMinTRDsig, binMaxTRDsig)), kTRDSignal);
296   histo->GetAxis(1)->Set(kMomentumBins - 1, momentumBins);
297   histo->GetAxis(0)->SetTitle("Particle Species");
298   histo->GetAxis(1)->SetTitle("p / GeV/c");
299   histo->GetAxis(2)->SetTitle("TPC Signal / a.u.");
300   histo->GetAxis(3)->SetTitle("TRD Signal / a.u.");
301   histo->GetAxis(4)->SetTitle("TRD Signal / a.u.");
302
303   // Add Histogram for combined TPC-ITS PID
304   const Int_t kDimensionsITSsig = 4;
305   Int_t kNbinsITSsig[kDimensionsITSsig] = {AliPID::kSPECIES + 1, kMomentumBins - 1, 200, 1000};
306   Double_t binMinITSsig[kDimensionsITSsig] = {-1., 0.1, 0., 0.};
307   Double_t binMaxITSsig[kDimensionsITSsig] = {AliPID::kSPECIES, 10., 200., 300.};
308   fQAlist->AddAt((histo = new THnSparseF("fCombTPCITSpid", "Combined TPC-ITS PID", kDimensionsITSsig, kNbinsITSsig, binMinITSsig, binMaxITSsig)), kITSSignal);
309   histo->GetAxis(0)->Set(kMomentumBins - 1, momentumBins);
310   histo->GetAxis(0)->SetTitle("Particle Species");
311   histo->GetAxis(1)->SetTitle("p / GeV/c");
312   histo->GetAxis(2)->SetTitle("TPC Signal / a.u.");
313   histo->GetAxis(3)->SetTitle("ITS Signal / a.u.");
314 }
315
316 //____________________________________________________________
317 void AliHFEpid::SetMCEvent(AliMCEvent *event){
318   //
319   // Set MC Event for the detector PID classes
320   //
321   for(Int_t idet = 0; idet < kNdetectorPID; idet++)
322     if(fDetectorPID[idet]) fDetectorPID[idet]->SetMCEvent(event);
323 }
324