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 /* $Id$ */
  17
  18 //
  19 // PID Steering Class
  20 // Interface to the user task
  21 // Several strategies for Electron identification implemented.
  22 // In addition users can combine different detectors or use
  23 // single detector PID
  24 //
  25 // Authors:
  26 //   Markus Fasel <M.Fasel@gsi.de>
  27 //
  28 #include <TAxis.h>
  29 #include <TClass.h>
  30 #include <TF1.h>
  31 #include <TIterator.h>
  32 #include <TList.h>
  33 #include <TObjArray.h>
  34 #include <TObjString.h>
  35 #include <TString.h>
  36
  37 #include "AliAODpidUtil.h"
  38 #include "AliESDpid.h"
  39 #include "AliLog.h"
  40 #include "AliPID.h"
  41 #include "AliVParticle.h"
  42
  43 #include "AliHFEcontainer.h"
  44 #include "AliHFEpid.h"
  45 #include "AliHFEpidBase.h"
  46 #include "AliHFEpidQAmanager.h"
  47 #include "AliHFEpidITS.h"
  48 #include "AliHFEpidTPC.h"
  49 #include "AliHFEpidTRD.h"
  50 #include "AliHFEpidTOF.h"
  51 #include "AliHFEpidMC.h"
  52 #include "AliHFEvarManager.h"
  53
  54 ClassImp(AliHFEpid)
  55
  56 const Char_t* AliHFEpid::fgkDetectorName[AliHFEpid::kNdetectorPID + 1] = {
  57   "MCPID",
  58   "ESDPID",
  59   "ITSPID",
  60   "TPCPID",
  61   "TRDPID",
  62   "TOFPID",
  63   "EMCALPID",
  64   "UndefinedPID"
  65 };
  66
  67 //____________________________________________________________
  68 AliHFEpid::AliHFEpid():
  69   TNamed(),
  70   fEnabledDetectors(0),
  71   fNPIDdetectors(0),
  72   fVarManager(NULL),
  73   fCommonObjects(NULL)
  74 {
  75   //
  76   // Default constructor
  77   //
  78   memset(fDetectorPID, 0, sizeof(AliHFEpidBase *) * kNdetectorPID);
  79   memset(fDetectorOrder, kUndefined, sizeof(UInt_t) * kNdetectorPID);
  80   memset(fDetectorOrder, 0, sizeof(UInt_t) * kNdetectorPID);
  81 }
  82
  83 //____________________________________________________________
  84 AliHFEpid::AliHFEpid(const Char_t *name):
  85   TNamed(name, ""),
  86   fEnabledDetectors(0),
  87   fNPIDdetectors(0),
  88   fVarManager(NULL),
  89   fCommonObjects(NULL)
  90 {
  91   //
  92   // Default constructor
  93   // Create PID objects for all detectors
  94   //
  95   memset(fDetectorPID, 0, sizeof(AliHFEpidBase *) * kNdetectorPID);
  96   memset(fDetectorOrder, kUndefined, sizeof(UInt_t) * kNdetectorPID);
  97   memset(fSortedOrder, 0, sizeof(UInt_t) * kNdetectorPID);
  98
  99   fDetectorPID[kMCpid] = new AliHFEpidMC("MCPID");
 100   fDetectorPID[kTPCpid] = new AliHFEpidTPC("TPCPID");
 101   fDetectorPID[kTRDpid] = new AliHFEpidTRD("TRDPID");
 102   fDetectorPID[kTOFpid] = new AliHFEpidTOF("TOFPID");
 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::AddDetector(TString detector, UInt_t position){
 189   //
 190   // Add Detector in position
 191   //
 192   UInt_t detectorID = kUndefined;
 193   detector.ToUpper();
 194   if(!detector.CompareTo("MC")) detectorID = kMCpid;
 195   else if(!detector.CompareTo("TPC")) detectorID = kTPCpid;
 196   else if(!detector.CompareTo("TRD")) detectorID = kTRDpid;
 197   else if(!detector.CompareTo("TOF")) detectorID = kTOFpid;
 198   else if(!detector.CompareTo("EMCAL")) detectorID = kEMCALpid;
 199   else AliError("Detector not available");
 200
 201   if(detectorID == kUndefined) return;
 202   if(IsDetectorOn(detectorID)) return;
 203   SwitchOnDetector(detectorID);
 204   fDetectorOrder[detectorID] = position;
 205   fNPIDdetectors++;
 206 }
 207
 208 //____________________________________________________________
 209 Bool_t AliHFEpid::InitializePID(){
 210   //
 211   // Initializes the PID object
 212   //
 213
 214   TMath::Sort(static_cast<UInt_t>(kNdetectorPID), fDetectorOrder, fSortedOrder, kFALSE);
 215   // Initialize PID Objects
 216   Bool_t status = kTRUE;
 217   for(Int_t idet = 0; idet < kNdetectorPID; idet++){
 218     if(!IsDetectorOn(idet)) continue;
 219     if(fDetectorPID[idet]){
 220       status &= fDetectorPID[idet]->InitializePID();
 221       if(HasMCData() && status) fDetectorPID[idet]->SetHasMCData();
 222     }
 223   }
 224   PrintStatus();
 225   return status;
 226 }
 227
 228 //____________________________________________________________
 229 Bool_t AliHFEpid::IsSelected(AliHFEpidObject *track, AliHFEcontainer *cont, const Char_t *contname, AliHFEpidQAmanager *pidqa){
 230   //
 231   // Select Tracks
 232   //
 233   Bool_t isSelected = kTRUE;
 234   AliDebug(1, Form("Particle used for PID, QA available: %s", pidqa ? "Yes" : "No"));
 235   for(UInt_t idet = 0; idet < fNPIDdetectors; idet++){
 236     AliDebug(2, Form("Using Detector %s\n", SortedDetectorName(idet)));
 237     if(TMath::Abs(fDetectorPID[fSortedOrder[idet]]->IsSelected(track, pidqa)) != 11){
 238       isSelected = kFALSE;
 239       break;
 240     }
 241     AliDebug(2, "Particlae selected by detector");
 242     if(fVarManager && cont){
 243       TString reccontname = contname; reccontname += "Reco";
 244       AliDebug(2, Form("Filling container %s", reccontname.Data()));
 245       if(fVarManager->IsSignalTrack())
 246         fVarManager->FillContainerStepname(cont, reccontname.Data(), SortedDetectorName(idet));
 247       if(HasMCData()){
 248         TString mccontname = contname; mccontname += "MC";
 249         AliDebug(2, Form("MC Information available, Filling container %s", mccontname.Data()));
 250         if(fVarManager->IsSignalTrack()) {
 251           fVarManager->FillContainerStepname(cont, mccontname.Data(), SortedDetectorName(idet), kTRUE);
 252                 if(cont->GetCorrelationMatrix("correlationstepafterTOF")){
 253                   TString tstept("TOFPID");
 254                   if(!tstept.CompareTo(SortedDetectorName(idet))) {
 255                     fVarManager->FillCorrelationMatrix(cont->GetCorrelationMatrix("correlationstepafterTOF"));
 256                     //printf("Step %s\n",(const char*) SortedDetectorName(idet));
 257                   }
 258                 }
 259               }
 260       }
 261       // The PID will NOT fill the double counting information
 262     }
 263   }
 264   return isSelected;
 265 }
 266
 267 //____________________________________________________________
 268 void AliHFEpid::SetESDpid(AliESDpid *pid){
 269   //
 270   // Set ESD PID to the Detector PID objects
 271   //
 272   for(Int_t idet = 0; idet < kNdetectorPID; idet++){
 273     if(fDetectorPID[idet]) fDetectorPID[idet]->SetESDpid(pid);
 274   }
 275 }
 276
 277 //____________________________________________________________
 278 void AliHFEpid::SetAODpid(AliAODpidUtil *pid){
 279   //
 280   // Set ESD PID to the Detector PID objects
 281   //
 282   for(Int_t idet = 0; idet < kNdetectorPID; idet++){
 283     if(fDetectorPID[idet]) fDetectorPID[idet]->SetAODpid(pid);
 284   }
 285 }
 286
 287 //____________________________________________________________
 288 void AliHFEpid::ConfigureTPCasymmetric(Double_t pmin, Double_t pmax, Double_t sigmamin, Double_t sigmamax){
 289   //
 290   // 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
 291   //
 292   AliHFEpidTPC *pid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
 293   if(pid){
 294     pid->SetTPCnSigma(3);
 295     pid->SetAsymmetricTPCsigmaCut(pmin, pmax, sigmamin, sigmamax);
 296   }
 297 }
 298
 299 //____________________________________________________________
 300 void AliHFEpid::ConfigureTPCrejectionSimple(){
 301   //
 302   // TPC alone, symmetric 3 sigma cut and 2 - -100 sigma pion rejection
 303   //
 304   AliHFEpidTPC *pid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
 305   if(pid){
 306     pid->SetTPCnSigma(3);
 307     pid->SetRejectParticle(AliPID::kPion, 0., -100., 10., 1.);
 308   }
 309 }
 310
 311 //____________________________________________________________
 312 void AliHFEpid::ConfigureTPCrejection(const char *lowerCutParam, Double_t *params){
 313   //
 314   // Combined TPC-TOF PID
 315   // if no function parameterizaion is given, then the default one (exponential) is chosen
 316   //
 317   if(HasMCData()) AliInfo("Configuring TPC for MC\n");
 318   AliHFEpidTPC *tpcpid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
 319   AliHFEpidTOF *tofpid = dynamic_cast<AliHFEpidTOF *>(fDetectorPID[kTOFpid]);
 320   if(tofpid) tofpid->SetTOFnSigma(3);
 321
 322   //TF1 *upperCut = new TF1("upperCut", "[0] * TMath::Exp([1]*x)", 0, 20);
 323   TF1 *upperCut = new TF1("upperCut", "[0]", 0, 20); // Use constant upper cut
 324   TF1 *lowerCut = new TF1("lowerCut", lowerCutParam == NULL ? "[0] * TMath::Exp([1]*x) + [2]": lowerCutParam, 0, 20);
 325   upperCut->SetParameter(0, 3.);
 326   //upperCut->SetParameter(0, 2.7);
 327   //upperCut->SetParameter(1, -0.4357);
 328   if(params){
 329     for(Int_t ipar = 0; ipar < lowerCut->GetNpar(); ipar++) lowerCut->SetParameter(ipar, params[ipar]);
 330   } else {
 331     // Set default parameterization
 332     if(HasMCData()) lowerCut->SetParameter(0, -2.5);
 333     else lowerCut->SetParameter(0, -3.71769);
 334     //else lowerCut->SetParameter(0, -3.7);
 335
 336     lowerCut->SetParameter(1, -0.40263);
 337     //lowerCut->SetParameter(1, -0.8);
 338
 339     if(HasMCData()) lowerCut->SetParameter(2, -2.2);
 340     else lowerCut->SetParameter(2, 0.267857);
 341     //else lowerCut->SetParameter(2, -0.35);
 342   }
 343
 344   if(tpcpid){
 345     tpcpid->SetTPCnSigma(2);
 346     tpcpid->SetUpperSigmaCut(upperCut);
 347     tpcpid->SetLowerSigmaCut(lowerCut);
 348   }
 349   AddCommonObject(upperCut);
 350   AddCommonObject(lowerCut);
 351 }
 352
 353 //____________________________________________________________
 354 void AliHFEpid::ConfigureTPCstrategyParis(){
 355   //
 356   // TPC alone, symmetric 3 sigma cut and 2 - -100 sigma pion rejection
 357   //
 358   AliHFEpidTPC *pid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
 359   if(pid){
 360     pid->SetTPCnSigma(2);
 361     pid->SetRejectParticle(AliPID::kProton, 0., -3., 10., 3.);
 362     pid->SetRejectParticle(AliPID::kKaon, 0., -3., 10., 3.);
 363   }
 364 }
 365
 366 //____________________________________________________________
 367 void AliHFEpid::PrintStatus() const {
 368   //
 369   // Print the PID configuration
 370   //
 371   printf("\n%s: Printing configuration\n", GetName());
 372   printf("===============================================\n");
 373   printf("PID Detectors: \n");
 374   Int_t npid = 0;
 375   TString detectors[kNdetectorPID] = {"MC", "ESD", "ITS", "TPC", "TRD", "TOF"};
 376   for(Int_t idet = 0; idet < kNdetectorPID; idet++){
 377     if(IsDetectorOn(idet)){
 378       printf("\t%s\n", detectors[idet].Data());
 379       npid++;
 380     }
 381   }
 382   if(!npid) printf("\tNone\n");
 383   printf("\n");
 384 }