Update of the HFE package
[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 // PID Steering Class 
17 // Interface to the user task
18 // Several strategies for Electron identification implemented.
19 // In addition users can combine different detectors or use
20 // single detector PID
21 // 
22 // Authors:
23 //   Markus Fasel <M.Fasel@gsi.de>
24 //
25 #include <TAxis.h>
26 #include <TClass.h>
27 #include <TF1.h>
28 #include <TIterator.h>
29 #include <TList.h>
30 #include <TObjArray.h>
31 #include <TObjString.h>
32 #include <TString.h>
33
34 #include "AliAODpidUtil.h"
35 #include "AliESDpid.h"
36 #include "AliLog.h"
37 #include "AliPID.h"
38 #include "AliVParticle.h"
39
40 #include "AliHFEcontainer.h"
41 #include "AliHFEpid.h"
42 #include "AliHFEpidQAmanager.h"
43 #include "AliHFEpidITS.h"
44 #include "AliHFEpidTPC.h"
45 #include "AliHFEpidTRD.h"
46 #include "AliHFEpidTOF.h"
47 #include "AliHFEpidEMCAL.h"
48 #include "AliHFEpidMC.h"
49 #include "AliHFEvarManager.h"
50
51 ClassImp(AliHFEpid)
52
53 const Char_t* AliHFEpid::fgkDetectorName[AliHFEpid::kNdetectorPID + 1] = {
54   "MCPID",
55   "ESDPID",
56   "ITSPID",
57   "TPCPID",
58   "TRDPID",
59   "TOFPID",
60   "EMCALPID",
61   "UndefinedPID"
62 };
63
64 //____________________________________________________________
65 AliHFEpid::AliHFEpid():
66   TNamed(),
67   fEnabledDetectors(0),
68   fNPIDdetectors(0),
69   fVarManager(NULL),
70   fCommonObjects(NULL)
71 {
72   //
73   // Default constructor
74   //
75   memset(fDetectorPID, 0, sizeof(AliHFEpidBase *) * kNdetectorPID);
76   memset(fDetectorOrder, kUndefined, sizeof(UInt_t) * kNdetectorPID);
77   memset(fDetectorOrder, 0, sizeof(UInt_t) * kNdetectorPID);
78 }
79
80 //____________________________________________________________
81 AliHFEpid::AliHFEpid(const Char_t *name):
82   TNamed(name, ""),
83   fEnabledDetectors(0),
84   fNPIDdetectors(0),
85   fVarManager(NULL),
86   fCommonObjects(NULL)
87 {
88   //
89   // Default constructor
90   // Create PID objects for all detectors
91   //
92   memset(fDetectorPID, 0, sizeof(AliHFEpidBase *) * kNdetectorPID);
93   memset(fDetectorOrder, kUndefined, sizeof(UInt_t) * kNdetectorPID);
94   memset(fSortedOrder, 0, sizeof(UInt_t) * kNdetectorPID);
95
96   fDetectorPID[kMCpid] = new AliHFEpidMC("MCPID");
97   fDetectorPID[kTPCpid] = new AliHFEpidTPC("TPCPID");
98   fDetectorPID[kTRDpid] = new AliHFEpidTRD("TRDPID");
99   fDetectorPID[kTOFpid] = new AliHFEpidTOF("TOFPID");
100   fDetectorPID[kEMCALpid] = new AliHFEpidEMCAL("EMCALPID");
101
102 }
103
104 //____________________________________________________________
105 AliHFEpid::AliHFEpid(const AliHFEpid &c):
106   TNamed(c),
107   fEnabledDetectors(c.fEnabledDetectors),
108   fNPIDdetectors(c.fNPIDdetectors),
109   fVarManager(c.fVarManager),
110   fCommonObjects(NULL)
111 {
112   //
113   // Copy Constructor
114   //
115   memset(fDetectorPID, 0, sizeof(AliHFEpidBase *) * kNdetectorPID);
116   c.Copy(*this);
117 }
118
119 //____________________________________________________________
120 AliHFEpid& AliHFEpid::operator=(const AliHFEpid &c){
121   //
122   // Assignment operator
123   //
124   if(&c != this) c.Copy(*this);
125   return *this;
126 }
127
128 //____________________________________________________________
129 AliHFEpid::~AliHFEpid(){
130   //
131   // Destructor
132   //
133   for(Int_t idet = 0; idet < kNdetectorPID; idet++)
134     if(fDetectorPID[idet]) delete fDetectorPID[idet];
135   ClearCommonObjects();
136 }
137
138 //____________________________________________________________
139 void AliHFEpid::Copy(TObject &o) const{
140   //
141   // Make copy
142   //
143   
144   TNamed::Copy(o);
145   AliHFEpid &target = dynamic_cast<AliHFEpid &>(o);
146   target.ClearCommonObjects();
147
148   target.fEnabledDetectors = fEnabledDetectors;
149   target.fNPIDdetectors = fNPIDdetectors;
150   target.fVarManager = fVarManager;
151  
152   // Copy detector PIDs
153   for(Int_t idet = 0; idet < kNdetectorPID; idet++){
154     //Cleanup pointers in case of assignment
155     if(target.fDetectorPID[idet])  
156       delete target.fDetectorPID[idet];     
157     if(fDetectorPID[idet]) 
158       target.fDetectorPID[idet] = dynamic_cast<AliHFEpidBase *>(fDetectorPID[idet]->Clone());
159   }
160   memcpy(target.fDetectorOrder, fDetectorOrder, sizeof(UInt_t) * kNdetectorPID);
161   memcpy(target.fSortedOrder, fSortedOrder, sizeof(UInt_t) * kNdetectorPID);
162 }
163
164 //____________________________________________________________
165 void AliHFEpid::AddCommonObject(TObject * const o){
166   //
167   // Add common object to the garbage collection
168   //
169   if(!fCommonObjects) fCommonObjects = new TObjArray;
170   fCommonObjects->Add(o);
171 }
172
173 //____________________________________________________________
174 void AliHFEpid::ClearCommonObjects(){
175   //
176   // empty garbage collection
177   //
178   if(fCommonObjects){
179     fCommonObjects->Delete();
180     delete fCommonObjects;
181     fCommonObjects = NULL;
182   }
183 }
184
185 //____________________________________________________________
186 void AliHFEpid::AddDetector(TString detector, UInt_t position){
187   //
188   // Add Detector in position 
189   //
190   UInt_t detectorID = kUndefined;
191   detector.ToUpper();
192   if(!detector.CompareTo("MC")) detectorID = kMCpid;
193   else if(!detector.CompareTo("TPC")) detectorID = kTPCpid;
194   else if(!detector.CompareTo("TRD")) detectorID = kTRDpid;
195   else if(!detector.CompareTo("TOF")) detectorID = kTOFpid;
196   else if(!detector.CompareTo("EMCAL")) detectorID = kEMCALpid;
197   else AliError("Detector not available");
198
199   if(detectorID == kUndefined) return;
200   if(IsDetectorOn(detectorID)) return;
201   SwitchOnDetector(detectorID);
202   fDetectorOrder[detectorID] = position;
203   fNPIDdetectors++;
204 }
205
206 //____________________________________________________________
207 Bool_t AliHFEpid::InitializePID(){
208   //
209   // Initializes the PID object
210   //
211
212   TMath::Sort(static_cast<UInt_t>(kNdetectorPID), fDetectorOrder, fSortedOrder, kFALSE);
213   // Initialize PID Objects
214   Bool_t status = kTRUE;
215   for(Int_t idet = 0; idet < kNdetectorPID; idet++){
216     if(!IsDetectorOn(idet)) continue;
217     if(fDetectorPID[idet]){ 
218       status &= fDetectorPID[idet]->InitializePID();
219       if(HasMCData() && status) fDetectorPID[idet]->SetHasMCData();
220     }
221   }
222   PrintStatus();
223   return status;
224 }
225
226 //____________________________________________________________
227 Bool_t AliHFEpid::IsSelected(const AliHFEpidObject * const track, AliHFEcontainer *cont, const Char_t *contname, AliHFEpidQAmanager *pidqa){
228   //
229   // Select Tracks
230   //
231   Bool_t isSelected = kTRUE;
232   AliDebug(1, Form("Particle used for PID, QA available: %s", pidqa ? "Yes" : "No"));
233   for(UInt_t idet = 0; idet < fNPIDdetectors; idet++){
234     AliDebug(2, Form("Using Detector %s\n", SortedDetectorName(idet)));
235     if(TMath::Abs(fDetectorPID[fSortedOrder[idet]]->IsSelected(track, pidqa)) != 11){
236       isSelected = kFALSE;
237       break;
238     }
239     AliDebug(2, "Particlae selected by detector");
240     if(fVarManager && cont){
241       TString reccontname = contname; reccontname += "Reco";
242       AliDebug(2, Form("Filling container %s", reccontname.Data()));
243       if(fVarManager->IsSignalTrack())
244         fVarManager->FillContainerStepname(cont, reccontname.Data(), SortedDetectorName(idet));
245       if(HasMCData()){
246         TString mccontname = contname; mccontname += "MC";
247         AliDebug(2, Form("MC Information available, Filling container %s", mccontname.Data()));
248         if(fVarManager->IsSignalTrack()) {
249           fVarManager->FillContainerStepname(cont, mccontname.Data(), SortedDetectorName(idet), kTRUE);
250                 if(cont->GetCorrelationMatrix("correlationstepafterTOF")){
251                   TString tstept("TOFPID"); 
252                   if(!tstept.CompareTo(SortedDetectorName(idet))) {
253                     fVarManager->FillCorrelationMatrix(cont->GetCorrelationMatrix("correlationstepafterTOF"));
254                     //printf("Step %s\n",(const char*) SortedDetectorName(idet));
255                   }
256                 }
257               }
258       }
259       // The PID will NOT fill the double counting information
260     }
261   }
262   return isSelected;
263 }
264
265 //____________________________________________________________
266 void AliHFEpid::SetESDpid(AliESDpid *pid){
267   //
268   // Set ESD PID to the Detector PID objects
269   //
270   for(Int_t idet = 0; idet < kNdetectorPID; idet++){
271     if(fDetectorPID[idet]) fDetectorPID[idet]->SetESDpid(pid);
272   }
273 }
274
275 //____________________________________________________________
276 void AliHFEpid::SetAODpid(AliAODpidUtil *pid){
277   //
278   // Set ESD PID to the Detector PID objects
279   //
280   for(Int_t idet = 0; idet < kNdetectorPID; idet++){
281     if(fDetectorPID[idet]) fDetectorPID[idet]->SetAODpid(pid);
282   }
283 }
284
285 //____________________________________________________________
286 void AliHFEpid::ConfigureTPCasymmetric(Double_t pmin, Double_t pmax, Double_t sigmamin, Double_t sigmamax){
287   //
288   // TPC alone, symmetric 3 sigma cut and asymmetric sigma cut in the momentum region between 2GeV/c and 10 GeV/c and sigma between -1 and 100
289   //
290   AliHFEpidTPC *pid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
291   if(pid){
292     pid->SetTPCnSigma(3);
293     pid->SetAsymmetricTPCsigmaCut(pmin, pmax, sigmamin, sigmamax);
294   }
295 }
296
297 //____________________________________________________________
298 void AliHFEpid::ConfigureTPCrejectionSimple(){
299   //
300   // TPC alone, symmetric 3 sigma cut and 2 - -100 sigma pion rejection
301   //   
302   AliHFEpidTPC *pid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
303   if(pid){
304     pid->SetTPCnSigma(3);
305     pid->SetRejectParticle(AliPID::kPion, 0., -100., 10., 1.);
306   }
307 }
308
309 //____________________________________________________________
310 void AliHFEpid::ConfigureTOF(Float_t TOFCut){
311   //
312   // Set Number of sigmas for TOF PID
313   //
314   AliHFEpidTOF *tofpid = dynamic_cast<AliHFEpidTOF *>(fDetectorPID[kTOFpid]);
315   if(tofpid) tofpid->SetTOFnSigma(TOFCut);
316 }
317
318 //____________________________________________________________
319 void AliHFEpid::ConfigureTPCcentralityCut(Int_t centralityBin, const char *lowerCutParam, const Double_t * const params, Float_t upperTPCCut){
320   //
321   // Cofigure centrality dependent cut function for TPC PID 
322   //
323   ConfigureTPCcut(centralityBin, lowerCutParam, params, upperTPCCut);
324 }
325
326 //____________________________________________________________
327 void AliHFEpid::ConfigureTPCdefaultCut(const char *lowerCutParam, const Double_t * const params, Float_t upperTPCCut){
328   //
329   // Cofigure default cut function for TPC PID 
330   //
331   ConfigureTPCcut(-1, lowerCutParam, params, upperTPCCut);
332 }
333
334 //____________________________________________________________
335 void AliHFEpid::ConfigureTPCcut(Int_t centralityBin, const char *lowerCutParam, const Double_t * const params, Float_t upperTPCCut){
336   //
337   // Cofigure cut function for TPC PID 
338   // if no function parameterizaion is given, then the default one (exponential) is chosen
339   //
340   
341   if(HasMCData()) AliInfo("Configuring TPC for MC\n");
342   AliHFEpidTPC *tpcpid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
343   //TF1 *upperCut = new TF1("upperCut", "[0] * TMath::Exp([1]*x)", 0, 20);
344   TF1 *upperCut = new TF1(Form("upperCut%s", centralityBin <  0 ? "Default" : Form("Bin%d", centralityBin)), "[0]", 0, 20); // Use constant upper cut
345   TF1 *lowerCut = new TF1(Form("lowerCut%s", centralityBin <  0 ? "Default" : Form("Bin%d", centralityBin)), lowerCutParam == NULL ? "[0] * TMath::Exp([1]*x) + [2]": lowerCutParam, 0, 20);
346
347   upperCut->SetParameter(0, upperTPCCut); // pp
348
349   if(params){
350       for(Int_t ipar = 0; ipar < lowerCut->GetNpar(); ipar++)
351       {
352           lowerCut->SetParameter(ipar, params[ipar]);
353         //  printf("printout %i %s %f \n", centralityBin, lowerCutParam, params[ipar]);
354       }
355   } else {
356     // Set default parameterization
357     if(HasMCData()) lowerCut->SetParameter(0, -2.5);
358     else lowerCut->SetParameter(0, -4.03);  //pp
359     lowerCut->SetParameter(1, -0.22); // pp
360     
361     if(HasMCData()) lowerCut->SetParameter(2, -2.2);
362     else lowerCut->SetParameter(2, 0.92); //pp
363   }
364
365
366   if(tpcpid){
367     tpcpid->SetTPCnSigma(2);
368     if(centralityBin < 0){
369       tpcpid->SetUpperSigmaCutDefault(upperCut);
370       tpcpid->SetLowerSigmaCutDefault(lowerCut);
371     } else {
372       tpcpid->SetUpperSigmaCutCentrality(upperCut, centralityBin);
373       tpcpid->SetLowerSigmaCutCentrality(lowerCut, centralityBin);
374     }
375   }
376   AddCommonObject(upperCut);
377   AddCommonObject(lowerCut);
378 }
379
380 //____________________________________________________________
381 void AliHFEpid::PrintStatus() const {
382   //
383   // Print the PID configuration
384   //
385   printf("\n%s: Printing configuration\n", GetName());
386   printf("===============================================\n");
387   printf("PID Detectors: \n");
388   Int_t npid = 0;
389   TString detectors[kNdetectorPID] = {"MC", "ESD", "ITS", "TPC", "TRD", "TOF", "EMCAL"};
390   for(Int_t idet = 0; idet < kNdetectorPID; idet++){
391     if(IsDetectorOn(idet)){
392       printf("\t%s\n", detectors[idet].Data());
393       npid++;
394     }
395   }
396   if(!npid) printf("\tNone\n");
397   printf("\n");
398 }