PWGHF/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 <TMath.h>
  31 #include <TObjArray.h>
  32 #include <TObjString.h>
  33 #include <TString.h>
  34
  35 #include "AliLog.h"
  36 #include "AliPID.h"
  37 #include "AliVParticle.h"
  38
  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"
  50
  51 ClassImp(AliHFEpid)
  52
  53 const Char_t* AliHFEpid::fgkDetectorName[AliHFEpid::kNdetectorPID + 1] = {
  54   "MCPID",
  55   "BAYESPID",
  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(fSortedOrder, 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[kBAYESpid] = new AliHFEpidBayes("BAYESPID");
  98   fDetectorPID[kTPCpid] = new AliHFEpidTPC("TPCPID");
  99   fDetectorPID[kTRDpid] = new AliHFEpidTRD("TRDPID");
 100   fDetectorPID[kITSpid] = new AliHFEpidITS("ITSPID");
 101   fDetectorPID[kTOFpid] = new AliHFEpidTOF("TOFPID");
 102   fDetectorPID[kEMCALpid] = new AliHFEpidEMCAL("EMCALPID");
 103
 104 }
 105
 106 //____________________________________________________________
 107 AliHFEpid::AliHFEpid(const AliHFEpid &c):
 108   TNamed(c),
 109   fEnabledDetectors(c.fEnabledDetectors),
 110   fNPIDdetectors(c.fNPIDdetectors),
 111   fVarManager(c.fVarManager),
 112   fCommonObjects(NULL)
 113 {
 114   //
 115   // Copy Constructor
 116   //
 117   memset(fDetectorPID, 0, sizeof(AliHFEpidBase *) * kNdetectorPID);
 118   c.Copy(*this);
 119 }
 120
 121 //____________________________________________________________
 122 AliHFEpid& AliHFEpid::operator=(const AliHFEpid &c){
 123   //
 124   // Assignment operator
 125   //
 126   if(&c != this) c.Copy(*this);
 127   return *this;
 128 }
 129
 130 //____________________________________________________________
 131 AliHFEpid::~AliHFEpid(){
 132   //
 133   // Destructor
 134   //
 135   for(Int_t idet = 0; idet < kNdetectorPID; idet++)
 136     if(fDetectorPID[idet]) delete fDetectorPID[idet];
 137   ClearCommonObjects();
 138 }
 139
 140 //____________________________________________________________
 141 void AliHFEpid::Copy(TObject &o) const{
 142   //
 143   // Make copy
 144   //
 145
 146   TNamed::Copy(o);
 147   AliHFEpid &target = dynamic_cast<AliHFEpid &>(o);
 148   target.ClearCommonObjects();
 149
 150   target.fEnabledDetectors = fEnabledDetectors;
 151   target.fNPIDdetectors = fNPIDdetectors;
 152   target.fVarManager = fVarManager;
 153
 154   // Copy detector PIDs
 155   for(Int_t idet = 0; idet < kNdetectorPID; idet++){
 156     //Cleanup pointers in case of assignment
 157     if(target.fDetectorPID[idet])
 158       delete target.fDetectorPID[idet];
 159     if(fDetectorPID[idet])
 160       target.fDetectorPID[idet] = dynamic_cast<AliHFEpidBase *>(fDetectorPID[idet]->Clone());
 161   }
 162   memcpy(target.fDetectorOrder, fDetectorOrder, sizeof(UInt_t) * kNdetectorPID);
 163   memcpy(target.fSortedOrder, fSortedOrder, sizeof(UInt_t) * kNdetectorPID);
 164 }
 165
 166 //____________________________________________________________
 167 void AliHFEpid::AddCommonObject(TObject * const o){
 168   //
 169   // Add common object to the garbage collection
 170   //
 171   if(!fCommonObjects) fCommonObjects = new TObjArray;
 172   fCommonObjects->Add(o);
 173 }
 174
 175 //____________________________________________________________
 176 void AliHFEpid::ClearCommonObjects(){
 177   //
 178   // empty garbage collection
 179   //
 180   if(fCommonObjects){
 181     fCommonObjects->Delete();
 182     delete fCommonObjects;
 183     fCommonObjects = NULL;
 184   }
 185 }
 186
 187 //____________________________________________________________
 188 void AliHFEpid::SetDetectorsForAnalysis(TString detectors){
 189   //
 190   // Set detectors used in Analysis to the position corresponding to their
 191   // position in the string
 192   // Detectors are separated by ","
 193   //
 194   TObjArray *detarray = detectors.Tokenize(",");
 195   TObjString *detString(NULL);
 196   int ndet(0);
 197   for(int idet = 0; idet < detarray->GetEntries(); idet++){
 198     detString = dynamic_cast<TObjString *>(detarray->At(idet));
 199     if(detString) AddDetector(detString->String(), ndet++);
 200   }
 201   AliDebug(1, Form("%d detectors used in Analysis", ndet));
 202 }
 203
 204 //____________________________________________________________
 205 void AliHFEpid::AddDetector(TString detector, UInt_t position){
 206   //
 207   // Add Detector in position
 208   //
 209   UInt_t detectorID = kUndefined;
 210   detector.ToUpper();
 211   if(!detector.CompareTo("MC")) detectorID = kMCpid;
 212   else if(!detector.CompareTo("BAYES")) detectorID = kBAYESpid;
 213   else if(!detector.CompareTo("TPC")) detectorID = kTPCpid;
 214   else if(!detector.CompareTo("TRD")) detectorID = kTRDpid;
 215   else if(!detector.CompareTo("TOF")) detectorID = kTOFpid;
 216   else if(!detector.CompareTo("ITS")) detectorID = kITSpid;
 217   else if(!detector.CompareTo("EMCAL")) detectorID = kEMCALpid;
 218   else AliError("Detector not available");
 219
 220   if(detectorID == kUndefined) return;
 221   if(IsDetectorOn(detectorID)) return;
 222   SwitchOnDetector(detectorID);
 223   fDetectorOrder[detectorID] = position;
 224   fNPIDdetectors++;
 225 }
 226
 227 //____________________________________________________________
 228 Bool_t AliHFEpid::InitializePID(Int_t run){
 229   //
 230   // Initializes the PID object
 231   //
 232
 233   if(!TestBit(kDetectorsSorted)) SortDetectors();
 234   Bool_t status = kTRUE;
 235   for(Int_t idet = 0; idet < kNdetectorPID; idet++){
 236     if(!IsDetectorOn(idet)) continue;
 237     if(fDetectorPID[idet]){
 238       status &= fDetectorPID[idet]->InitializePID(run);
 239       if(HasMCData() && status) fDetectorPID[idet]->SetHasMCData();
 240     }
 241   }
 242   SetBit(kIsInit);
 243   PrintStatus();
 244   return status;
 245 }
 246
 247 //____________________________________________________________
 248 Bool_t AliHFEpid::IsSelected(const AliHFEpidObject * const track, AliHFEcontainer *cont, const Char_t *contname, AliHFEpidQAmanager *pidqa){
 249   //
 250   // Select Tracks
 251   //
 252   Bool_t isSelected = kTRUE;
 253   AliDebug(1, Form("Particle used for PID, QA available: %s", pidqa ? "Yes" : "No"));
 254   for(UInt_t idet = 0; idet < fNPIDdetectors; idet++){
 255     AliDebug(2, Form("Using Detector %s\n", SortedDetectorName(idet)));
 256     if(TMath::Abs(fDetectorPID[fSortedOrder[idet]]->IsSelected(track, pidqa)) != 11){
 257       isSelected = kFALSE;
 258       break;
 259     }
 260     AliDebug(2, "Particlae selected by detector");
 261     if(fVarManager && cont){
 262       TString reccontname = contname; reccontname += "Reco";
 263       AliDebug(2, Form("Filling container %s", reccontname.Data()));
 264       if(fVarManager->IsSignalTrack())
 265         fVarManager->FillContainerStepname(cont, reccontname.Data(), SortedDetectorName(idet));
 266       if(HasMCData()){
 267         TString mccontname = contname; mccontname += "MC";
 268         AliDebug(2, Form("MC Information available, Filling container %s", mccontname.Data()));
 269         if(fVarManager->IsSignalTrack()) {
 270           fVarManager->FillContainerStepname(cont, mccontname.Data(), SortedDetectorName(idet), kTRUE);
 271                 if(cont->GetCorrelationMatrix("correlationstepafterTOF")){
 272                   TString tstept("TOFPID");
 273                   if(!tstept.CompareTo(SortedDetectorName(idet))) {
 274                     fVarManager->FillCorrelationMatrix(cont->GetCorrelationMatrix("correlationstepafterTOF"));
 275                     //printf("Step %s\n",(const char*) SortedDetectorName(idet));
 276                   }
 277                 }
 278               }
 279       }
 280       // The PID will NOT fill the double counting information
 281     }
 282   }
 283   return isSelected;
 284 }
 285
 286 //____________________________________________________________
 287 void AliHFEpid::SortDetectors(){
 288   //
 289   // Make sorted list of detectors
 290   //
 291   if(TestBit(kDetectorsSorted)) return; // Don't sort detectors when they are already sorted
 292   TMath::Sort(static_cast<UInt_t>(kNdetectorPID), fDetectorOrder, fSortedOrder, kFALSE);
 293   SetBit(kDetectorsSorted);
 294 }
 295
 296 //____________________________________________________________
 297 void AliHFEpid::SetPIDResponse(const AliPIDResponse * const pid){
 298   //
 299   // Set ESD PID to the Detector PID objects
 300   //
 301   for(Int_t idet = 0; idet < kNdetectorPID; idet++){
 302     if(fDetectorPID[idet]) fDetectorPID[idet]->SetPIDResponse(pid);
 303   }
 304 }
 305
 306 //____________________________________________________________
 307 const AliPIDResponse *AliHFEpid::GetPIDResponse() const {
 308   //
 309   // Return PID response function
 310   //
 311   const AliPIDResponse *response = NULL;
 312   for(Int_t idet = 0; idet < kNdetectorPID; idet++){
 313     if(fDetectorPID[idet]){
 314       response = fDetectorPID[idet]->GetPIDResponse();
 315       break;
 316     }
 317   }
 318   return response;
 319 }
 320
 321 //____________________________________________________________
 322 void AliHFEpid::ConfigureTPCasymmetric(Double_t pmin, Double_t pmax, Double_t sigmamin, Double_t sigmamax){
 323   //
 324   // 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
 325   //
 326   AliHFEpidTPC *pid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
 327   if(pid){
 328     pid->SetTPCnSigma(3);
 329     pid->SetAsymmetricTPCsigmaCut(pmin, pmax, sigmamin, sigmamax);
 330   }
 331 }
 332
 333 //____________________________________________________________
 334 void AliHFEpid::ConfigureTPCrejectionSimple(){
 335   //
 336   // TPC alone, symmetric 3 sigma cut and 2 - -100 sigma pion rejection
 337   //
 338   AliHFEpidTPC *pid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
 339   if(pid){
 340     pid->SetTPCnSigma(3);
 341     pid->SetRejectParticle(AliPID::kPion, 0., -100., 10., 1.);
 342   }
 343 }
 344
 345 //____________________________________________________________
 346 void AliHFEpid::ConfigureTOF(Float_t TOFCut){
 347   //
 348   // Set Number of sigmas for TOF PID
 349   //
 350   AliHFEpidTOF *tofpid = dynamic_cast<AliHFEpidTOF *>(fDetectorPID[kTOFpid]);
 351   if(tofpid) tofpid->SetTOFnSigma(TOFCut);
 352 }
 353
 354 //____________________________________________________________
 355 void AliHFEpid::ConfigureTPCcentralityCut(Int_t centralityBin, const char *lowerCutParam, const Double_t * const params, Float_t upperTPCCut){
 356   //
 357   // Cofigure centrality dependent cut function for TPC PID
 358   //
 359   ConfigureTPCcut(centralityBin, lowerCutParam, params, upperTPCCut);
 360 }
 361
 362 //____________________________________________________________
 363 void AliHFEpid::ConfigureTPCdefaultCut(const char *lowerCutParam, const Double_t * const params, Float_t upperTPCCut){
 364   //
 365   // Cofigure default cut function for TPC PID
 366   //
 367   ConfigureTPCcut(-1, lowerCutParam, params, upperTPCCut);
 368 }
 369
 370 //____________________________________________________________
 371 void AliHFEpid::ConfigureTPCcut(Int_t centralityBin, const char *lowerCutParam, const Double_t * const params, Float_t upperTPCCut){
 372   //
 373   // Cofigure cut function for TPC PID
 374   // if no function parameterizaion is given, then the default one (exponential) is chosen
 375   //
 376
 377   if(HasMCData()) AliInfo("Configuring TPC for MC\n");
 378   AliHFEpidTPC *tpcpid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
 379   //TF1 *upperCut = new TF1("upperCut", "[0] * TMath::Exp([1]*x)", 0, 20);
 380   TF1 *upperCut = new TF1(Form("upperCut%s", centralityBin <  0 ? "Default" : Form("Bin%d", centralityBin)), "[0]", 0, 20); // Use constant upper cut
 381   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);
 382
 383   upperCut->SetParameter(0, upperTPCCut); // pp
 384
 385   if(params){
 386       for(Int_t ipar = 0; ipar < lowerCut->GetNpar(); ipar++)
 387       {
 388           lowerCut->SetParameter(ipar, params[ipar]);
 389         //  printf("printout %i %s %f \n", centralityBin, lowerCutParam, params[ipar]);
 390       }
 391   } else {
 392     // Set default parameterization
 393     if(HasMCData()) lowerCut->SetParameter(0, -2.5);
 394     else lowerCut->SetParameter(0, -4.03);  //pp
 395     lowerCut->SetParameter(1, -0.22); // pp
 396
 397     if(HasMCData()) lowerCut->SetParameter(2, -2.2);
 398     else lowerCut->SetParameter(2, 0.92); //pp
 399   }
 400
 401
 402   if(tpcpid){
 403     tpcpid->SetTPCnSigma(2);
 404     if(centralityBin < 0){
 405       tpcpid->SetUpperSigmaCutDefault(upperCut);
 406       tpcpid->SetLowerSigmaCutDefault(lowerCut);
 407     } else {
 408       tpcpid->SetUpperSigmaCutCentrality(upperCut, centralityBin);
 409       tpcpid->SetLowerSigmaCutCentrality(lowerCut, centralityBin);
 410     }
 411   }
 412   AddCommonObject(upperCut);
 413   AddCommonObject(lowerCut);
 414 }
 415
 416 //____________________________________________________________
 417 void AliHFEpid::ConfigureBayesDetectorMask(Int_t detmask){
 418   //
 419   // Configure detector mask for Bayes PID
 420   // if no detector mask is set the default mask is chosen
 421   //
 422
 423   if(HasMCData()) AliInfo("Configuring Bayes for MC\n");
 424   AliHFEpidBayes *bayespid = dynamic_cast<AliHFEpidBayes *>(fDetectorPID[kBAYESpid]);
 425
 426   if(bayespid)
 427   {
 428       bayespid->SetBayesDetectorMask(detmask);
 429       printf("detector mask in pid class %i \n",detmask);
 430   }
 431
 432 }
 433
 434 //____________________________________________________________
 435 void AliHFEpid::ConfigureBayesPIDThreshold(Float_t pidthres){
 436   //
 437   // Configure pid threshold for Bayes PID
 438   // if no threshold is set the default threshold is chosen
 439   //
 440
 441   if(HasMCData()) AliInfo("Configuring Bayes for MC\n");
 442   AliHFEpidBayes *bayespid = dynamic_cast<AliHFEpidBayes *>(fDetectorPID[kBAYESpid]);
 443
 444   if(bayespid)
 445   {
 446       bayespid->SetBayesPIDThreshold(pidthres);
 447       printf("combined pidthreshold %f \n",pidthres);
 448   }
 449
 450 }
 451
 452 //____________________________________________________________
 453 void AliHFEpid::PrintStatus() const {
 454   //
 455   // Print the PID configuration
 456   //
 457   printf("\n%s: Printing configuration\n", GetName());
 458   printf("===============================================\n");
 459   printf("PID Detectors: \n");
 460   Int_t npid = 0;
 461   TString detectors[kNdetectorPID] = {"MC", "BAYES", "ITS", "TPC", "TRD", "TOF", "EMCAL"};
 462   for(Int_t idet = 0; idet < kNdetectorPID; idet++){
 463     if(IsDetectorOn(idet)){
 464       printf("\t%s\n", detectors[idet].Data());
 465       npid++;
 466     }
 467   }
 468   if(!npid) printf("\tNone\n");
 469   printf("\n");
 470 }