1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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
23 // Markus Fasel <M.Fasel@gsi.de>
28 #include <TIterator.h>
31 #include <TObjArray.h>
32 #include <TObjString.h>
37 #include "AliVParticle.h"
39 #include "AliHFEcontainer.h"
40 #include "AliHFEpid.h"
41 #include "AliHFEpidQAmanager.h"
42 #include "AliHFEpidITS.h"
43 #include "AliHFEpidTPC.h"
44 #include "AliHFEpidTRD.h"
45 #include "AliHFEpidTOF.h"
46 #include "AliHFEpidEMCAL.h"
47 #include "AliHFEpidMC.h"
48 #include "AliHFEpidBayes.h"
49 #include "AliHFEvarManager.h"
53 const Char_t* AliHFEpid::fgkDetectorName[AliHFEpid::kNdetectorPID + 1] = {
64 //____________________________________________________________
65 AliHFEpid::AliHFEpid():
73 // Default constructor
75 memset(fDetectorPID, 0, sizeof(AliHFEpidBase *) * kNdetectorPID);
76 memset(fDetectorOrder, kUndefined, sizeof(UInt_t) * kNdetectorPID);
77 memset(fSortedOrder, 0, sizeof(UInt_t) * kNdetectorPID);
80 //____________________________________________________________
81 AliHFEpid::AliHFEpid(const Char_t *name):
89 // Default constructor
90 // Create PID objects for all detectors
92 memset(fDetectorPID, 0, sizeof(AliHFEpidBase *) * kNdetectorPID);
93 memset(fDetectorOrder, kUndefined, sizeof(UInt_t) * kNdetectorPID);
94 memset(fSortedOrder, 0, sizeof(UInt_t) * kNdetectorPID);
96 fDetectorPID[kMCpid] = new AliHFEpidMC("MCPID");
97 fDetectorPID[kBAYESpid] = new AliHFEpidBayes("BAYESPID");
98 fDetectorPID[kTPCpid] = new AliHFEpidTPC("TPCPID");
99 fDetectorPID[kTRDpid] = new AliHFEpidTRD("TRDPID");
100 fDetectorPID[kTOFpid] = new AliHFEpidTOF("TOFPID");
101 fDetectorPID[kEMCALpid] = new AliHFEpidEMCAL("EMCALPID");
105 //____________________________________________________________
106 AliHFEpid::AliHFEpid(const AliHFEpid &c):
108 fEnabledDetectors(c.fEnabledDetectors),
109 fNPIDdetectors(c.fNPIDdetectors),
110 fVarManager(c.fVarManager),
116 memset(fDetectorPID, 0, sizeof(AliHFEpidBase *) * kNdetectorPID);
120 //____________________________________________________________
121 AliHFEpid& AliHFEpid::operator=(const AliHFEpid &c){
123 // Assignment operator
125 if(&c != this) c.Copy(*this);
129 //____________________________________________________________
130 AliHFEpid::~AliHFEpid(){
134 for(Int_t idet = 0; idet < kNdetectorPID; idet++)
135 if(fDetectorPID[idet]) delete fDetectorPID[idet];
136 ClearCommonObjects();
139 //____________________________________________________________
140 void AliHFEpid::Copy(TObject &o) const{
146 AliHFEpid &target = dynamic_cast<AliHFEpid &>(o);
147 target.ClearCommonObjects();
149 target.fEnabledDetectors = fEnabledDetectors;
150 target.fNPIDdetectors = fNPIDdetectors;
151 target.fVarManager = fVarManager;
153 // Copy detector PIDs
154 for(Int_t idet = 0; idet < kNdetectorPID; idet++){
155 //Cleanup pointers in case of assignment
156 if(target.fDetectorPID[idet])
157 delete target.fDetectorPID[idet];
158 if(fDetectorPID[idet])
159 target.fDetectorPID[idet] = dynamic_cast<AliHFEpidBase *>(fDetectorPID[idet]->Clone());
161 memcpy(target.fDetectorOrder, fDetectorOrder, sizeof(UInt_t) * kNdetectorPID);
162 memcpy(target.fSortedOrder, fSortedOrder, sizeof(UInt_t) * kNdetectorPID);
165 //____________________________________________________________
166 void AliHFEpid::AddCommonObject(TObject * const o){
168 // Add common object to the garbage collection
170 if(!fCommonObjects) fCommonObjects = new TObjArray;
171 fCommonObjects->Add(o);
174 //____________________________________________________________
175 void AliHFEpid::ClearCommonObjects(){
177 // empty garbage collection
180 fCommonObjects->Delete();
181 delete fCommonObjects;
182 fCommonObjects = NULL;
186 //____________________________________________________________
187 void AliHFEpid::SetDetectorsForAnalysis(TString detectors){
189 // Set detectors used in Analysis to the position corresponding to their
190 // position in the string
191 // Detectors are separated by ","
193 TObjArray *detarray = detectors.Tokenize(",");
194 TObjString *detString(NULL);
196 for(int idet = 0; idet < detarray->GetEntries(); idet++){
197 detString = dynamic_cast<TObjString *>(detarray->At(idet));
198 if(detString) AddDetector(detString->String(), ndet++);
200 AliDebug(1, Form("%d detectors used in Analysis", ndet));
203 //____________________________________________________________
204 void AliHFEpid::AddDetector(TString detector, UInt_t position){
206 // Add Detector in position
208 UInt_t detectorID = kUndefined;
210 if(!detector.CompareTo("MC")) detectorID = kMCpid;
211 else if(!detector.CompareTo("BAYES")) detectorID = kBAYESpid;
212 else if(!detector.CompareTo("TPC")) detectorID = kTPCpid;
213 else if(!detector.CompareTo("TRD")) detectorID = kTRDpid;
214 else if(!detector.CompareTo("TOF")) detectorID = kTOFpid;
215 else if(!detector.CompareTo("EMCAL")) detectorID = kEMCALpid;
216 else AliError("Detector not available");
218 if(detectorID == kUndefined) return;
219 if(IsDetectorOn(detectorID)) return;
220 SwitchOnDetector(detectorID);
221 fDetectorOrder[detectorID] = position;
225 //____________________________________________________________
226 Bool_t AliHFEpid::InitializePID(Int_t run){
228 // Initializes the PID object
231 if(!TestBit(kDetectorsSorted)) SortDetectors();
232 Bool_t status = kTRUE;
233 for(Int_t idet = 0; idet < kNdetectorPID; idet++){
234 if(!IsDetectorOn(idet)) continue;
235 if(fDetectorPID[idet]){
236 status &= fDetectorPID[idet]->InitializePID(run);
237 if(HasMCData() && status) fDetectorPID[idet]->SetHasMCData();
245 //____________________________________________________________
246 Bool_t AliHFEpid::IsSelected(const AliHFEpidObject * const track, AliHFEcontainer *cont, const Char_t *contname, AliHFEpidQAmanager *pidqa){
250 Bool_t isSelected = kTRUE;
251 AliDebug(1, Form("Particle used for PID, QA available: %s", pidqa ? "Yes" : "No"));
252 for(UInt_t idet = 0; idet < fNPIDdetectors; idet++){
253 AliDebug(2, Form("Using Detector %s\n", SortedDetectorName(idet)));
254 if(TMath::Abs(fDetectorPID[fSortedOrder[idet]]->IsSelected(track, pidqa)) != 11){
258 AliDebug(2, "Particlae selected by detector");
259 if(fVarManager && cont){
260 TString reccontname = contname; reccontname += "Reco";
261 AliDebug(2, Form("Filling container %s", reccontname.Data()));
262 if(fVarManager->IsSignalTrack())
263 fVarManager->FillContainerStepname(cont, reccontname.Data(), SortedDetectorName(idet));
265 TString mccontname = contname; mccontname += "MC";
266 AliDebug(2, Form("MC Information available, Filling container %s", mccontname.Data()));
267 if(fVarManager->IsSignalTrack()) {
268 fVarManager->FillContainerStepname(cont, mccontname.Data(), SortedDetectorName(idet), kTRUE);
269 if(cont->GetCorrelationMatrix("correlationstepafterTOF")){
270 TString tstept("TOFPID");
271 if(!tstept.CompareTo(SortedDetectorName(idet))) {
272 fVarManager->FillCorrelationMatrix(cont->GetCorrelationMatrix("correlationstepafterTOF"));
273 //printf("Step %s\n",(const char*) SortedDetectorName(idet));
278 // The PID will NOT fill the double counting information
284 //____________________________________________________________
285 void AliHFEpid::SortDetectors(){
287 // Make sorted list of detectors
289 if(TestBit(kDetectorsSorted)) return; // Don't sort detectors when they are already sorted
290 TMath::Sort(static_cast<UInt_t>(kNdetectorPID), fDetectorOrder, fSortedOrder, kFALSE);
291 SetBit(kDetectorsSorted);
294 //____________________________________________________________
295 void AliHFEpid::SetPIDResponse(const AliPIDResponse * const pid){
297 // Set ESD PID to the Detector PID objects
299 for(Int_t idet = 0; idet < kNdetectorPID; idet++){
300 if(fDetectorPID[idet]) fDetectorPID[idet]->SetPIDResponse(pid);
304 //____________________________________________________________
305 const AliPIDResponse *AliHFEpid::GetPIDResponse() const {
307 // Return PID response function
309 const AliPIDResponse *response = NULL;
310 for(Int_t idet = 0; idet < kNdetectorPID; idet++){
311 if(fDetectorPID[idet]){
312 response = fDetectorPID[idet]->GetPIDResponse();
319 //____________________________________________________________
320 void AliHFEpid::ConfigureTPCasymmetric(Double_t pmin, Double_t pmax, Double_t sigmamin, Double_t sigmamax){
322 // 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
324 AliHFEpidTPC *pid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
326 pid->SetTPCnSigma(3);
327 pid->SetAsymmetricTPCsigmaCut(pmin, pmax, sigmamin, sigmamax);
331 //____________________________________________________________
332 void AliHFEpid::ConfigureTPCrejectionSimple(){
334 // TPC alone, symmetric 3 sigma cut and 2 - -100 sigma pion rejection
336 AliHFEpidTPC *pid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
338 pid->SetTPCnSigma(3);
339 pid->SetRejectParticle(AliPID::kPion, 0., -100., 10., 1.);
343 //____________________________________________________________
344 void AliHFEpid::ConfigureTOF(Float_t TOFCut){
346 // Set Number of sigmas for TOF PID
348 AliHFEpidTOF *tofpid = dynamic_cast<AliHFEpidTOF *>(fDetectorPID[kTOFpid]);
349 if(tofpid) tofpid->SetTOFnSigma(TOFCut);
352 //____________________________________________________________
353 void AliHFEpid::ConfigureTPCcentralityCut(Int_t centralityBin, const char *lowerCutParam, const Double_t * const params, Float_t upperTPCCut){
355 // Cofigure centrality dependent cut function for TPC PID
357 ConfigureTPCcut(centralityBin, lowerCutParam, params, upperTPCCut);
360 //____________________________________________________________
361 void AliHFEpid::ConfigureTPCdefaultCut(const char *lowerCutParam, const Double_t * const params, Float_t upperTPCCut){
363 // Cofigure default cut function for TPC PID
365 ConfigureTPCcut(-1, lowerCutParam, params, upperTPCCut);
368 //____________________________________________________________
369 void AliHFEpid::ConfigureTPCcut(Int_t centralityBin, const char *lowerCutParam, const Double_t * const params, Float_t upperTPCCut){
371 // Cofigure cut function for TPC PID
372 // if no function parameterizaion is given, then the default one (exponential) is chosen
375 if(HasMCData()) AliInfo("Configuring TPC for MC\n");
376 AliHFEpidTPC *tpcpid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
377 //TF1 *upperCut = new TF1("upperCut", "[0] * TMath::Exp([1]*x)", 0, 20);
378 TF1 *upperCut = new TF1(Form("upperCut%s", centralityBin < 0 ? "Default" : Form("Bin%d", centralityBin)), "[0]", 0, 20); // Use constant upper cut
379 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);
381 upperCut->SetParameter(0, upperTPCCut); // pp
384 for(Int_t ipar = 0; ipar < lowerCut->GetNpar(); ipar++)
386 lowerCut->SetParameter(ipar, params[ipar]);
387 // printf("printout %i %s %f \n", centralityBin, lowerCutParam, params[ipar]);
390 // Set default parameterization
391 if(HasMCData()) lowerCut->SetParameter(0, -2.5);
392 else lowerCut->SetParameter(0, -4.03); //pp
393 lowerCut->SetParameter(1, -0.22); // pp
395 if(HasMCData()) lowerCut->SetParameter(2, -2.2);
396 else lowerCut->SetParameter(2, 0.92); //pp
401 tpcpid->SetTPCnSigma(2);
402 if(centralityBin < 0){
403 tpcpid->SetUpperSigmaCutDefault(upperCut);
404 tpcpid->SetLowerSigmaCutDefault(lowerCut);
406 tpcpid->SetUpperSigmaCutCentrality(upperCut, centralityBin);
407 tpcpid->SetLowerSigmaCutCentrality(lowerCut, centralityBin);
410 AddCommonObject(upperCut);
411 AddCommonObject(lowerCut);
414 //____________________________________________________________
415 void AliHFEpid::ConfigureBayesDetectorMask(Int_t detmask){
417 // Configure detector mask for Bayes PID
418 // if no detector mask is set the default mask is chosen
421 if(HasMCData()) AliInfo("Configuring Bayes for MC\n");
422 AliHFEpidBayes *bayespid = dynamic_cast<AliHFEpidBayes *>(fDetectorPID[kBAYESpid]);
426 bayespid->SetBayesDetectorMask(detmask);
427 printf("detector mask in pid class %i \n",detmask);
432 //____________________________________________________________
433 void AliHFEpid::ConfigureBayesPIDThreshold(Float_t pidthres){
435 // Configure pid threshold for Bayes PID
436 // if no threshold is set the default threshold is chosen
439 if(HasMCData()) AliInfo("Configuring Bayes for MC\n");
440 AliHFEpidBayes *bayespid = dynamic_cast<AliHFEpidBayes *>(fDetectorPID[kBAYESpid]);
444 bayespid->SetBayesPIDThreshold(pidthres);
445 printf("combined pidthreshold %f \n",pidthres);
450 //____________________________________________________________
451 void AliHFEpid::PrintStatus() const {
453 // Print the PID configuration
455 printf("\n%s: Printing configuration\n", GetName());
456 printf("===============================================\n");
457 printf("PID Detectors: \n");
459 TString detectors[kNdetectorPID] = {"MC", "BAYES", "ITS", "TPC", "TRD", "TOF", "EMCAL"};
460 for(Int_t idet = 0; idet < kNdetectorPID; idet++){
461 if(IsDetectorOn(idet)){
462 printf("\t%s\n", detectors[idet].Data());
466 if(!npid) printf("\tNone\n");