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 "AliHFEpidBase.h"
  43 #include "AliHFEpidQAmanager.h"
  44 #include "AliHFEpidITS.h"
  45 #include "AliHFEpidTPC.h"
  46 #include "AliHFEpidTRD.h"
  47 #include "AliHFEpidTOF.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   "UndefinedPID"
  61 };
  62
  63 //____________________________________________________________
  64 AliHFEpid::AliHFEpid():
  65   TNamed(),
  66   fEnabledDetectors(0),
  67   fNPIDdetectors(0),
  68   fVarManager(NULL),
  69   fCommonObjects(NULL)
  70 {
  71   //
  72   // Default constructor
  73   //
  74   memset(fDetectorPID, 0, sizeof(AliHFEpidBase *) * kNdetectorPID);
  75   memset(fDetectorOrder, kUndefined, sizeof(UInt_t) * kNdetectorPID);
  76   memset(fDetectorOrder, 0, sizeof(UInt_t) * kNdetectorPID);
  77 }
  78
  79 //____________________________________________________________
  80 AliHFEpid::AliHFEpid(const Char_t *name):
  81   TNamed(name, ""),
  82   fEnabledDetectors(0),
  83   fNPIDdetectors(0),
  84   fVarManager(NULL),
  85   fCommonObjects(NULL)
  86 {
  87   //
  88   // Default constructor
  89   // Create PID objects for all detectors
  90   //
  91   memset(fDetectorPID, 0, sizeof(AliHFEpidBase *) * kNdetectorPID);
  92   memset(fDetectorOrder, kUndefined, sizeof(UInt_t) * kNdetectorPID);
  93   memset(fSortedOrder, 0, sizeof(UInt_t) * kNdetectorPID);
  94
  95   fDetectorPID[kMCpid] = new AliHFEpidMC("MCPID");
  96   fDetectorPID[kTPCpid] = new AliHFEpidTPC("TPCPID");
  97   fDetectorPID[kTRDpid] = new AliHFEpidTRD("TRDPID");
  98   fDetectorPID[kTOFpid] = new AliHFEpidTOF("TOFPID");
  99
 100 }
 101
 102 //____________________________________________________________
 103 AliHFEpid::AliHFEpid(const AliHFEpid &c):
 104   TNamed(c),
 105   fEnabledDetectors(c.fEnabledDetectors),
 106   fNPIDdetectors(c.fNPIDdetectors),
 107   fVarManager(c.fVarManager),
 108   fCommonObjects(NULL)
 109 {
 110   //
 111   // Copy Constructor
 112   //
 113   memset(fDetectorPID, 0, sizeof(AliHFEpidBase *) * kNdetectorPID);
 114   c.Copy(*this);
 115 }
 116
 117 //____________________________________________________________
 118 AliHFEpid& AliHFEpid::operator=(const AliHFEpid &c){
 119   //
 120   // Assignment operator
 121   //
 122   if(&c != this) c.Copy(*this);
 123   return *this;
 124 }
 125
 126 //____________________________________________________________
 127 AliHFEpid::~AliHFEpid(){
 128   //
 129   // Destructor
 130   //
 131   for(Int_t idet = 0; idet < kNdetectorPID; idet++)
 132     if(fDetectorPID[idet]) delete fDetectorPID[idet];
 133   ClearCommonObjects();
 134 }
 135
 136 //____________________________________________________________
 137 void AliHFEpid::Copy(TObject &o) const{
 138   //
 139   // Make copy
 140   //
 141
 142   TNamed::Copy(o);
 143   AliHFEpid &target = dynamic_cast<AliHFEpid &>(o);
 144   target.ClearCommonObjects();
 145
 146   target.fEnabledDetectors = fEnabledDetectors;
 147   target.fNPIDdetectors = fNPIDdetectors;
 148   target.fVarManager = fVarManager;
 149
 150   // Copy detector PIDs
 151   for(Int_t idet = 0; idet < kNdetectorPID; idet++){
 152     //Cleanup pointers in case of assignment
 153     if(target.fDetectorPID[idet])
 154       delete target.fDetectorPID[idet];
 155     if(fDetectorPID[idet])
 156       target.fDetectorPID[idet] = dynamic_cast<AliHFEpidBase *>(fDetectorPID[idet]->Clone());
 157   }
 158   memcpy(target.fDetectorOrder, fDetectorOrder, sizeof(UInt_t) * kNdetectorPID);
 159   memcpy(target.fSortedOrder, fSortedOrder, sizeof(UInt_t) * kNdetectorPID);
 160 }
 161
 162 //____________________________________________________________
 163 void AliHFEpid::AddCommonObject(TObject * const o){
 164   //
 165   // Add common object to the garbage collection
 166   //
 167   if(!fCommonObjects) fCommonObjects = new TObjArray;
 168   fCommonObjects->Add(o);
 169 }
 170
 171 //____________________________________________________________
 172 void AliHFEpid::ClearCommonObjects(){
 173   //
 174   // empty garbage collection
 175   //
 176   if(fCommonObjects){
 177     fCommonObjects->Delete();
 178     delete fCommonObjects;
 179     fCommonObjects = NULL;
 180   }
 181 }
 182
 183 //____________________________________________________________
 184 void AliHFEpid::AddDetector(TString detector, UInt_t position){
 185   //
 186   // Add Detector in position
 187   //
 188   UInt_t detectorID = kUndefined;
 189   detector.ToUpper();
 190   if(!detector.CompareTo("MC")) detectorID = kMCpid;
 191   else if(!detector.CompareTo("TPC")) detectorID = kTPCpid;
 192   else if(!detector.CompareTo("TRD")) detectorID = kTRDpid;
 193   else if(!detector.CompareTo("TOF")) detectorID = kTOFpid;
 194   else AliError("Detector not available");
 195
 196   if(detectorID == kUndefined) return;
 197   if(IsDetectorOn(detectorID)) return;
 198   SwitchOnDetector(detectorID);
 199   fDetectorOrder[detectorID] = position;
 200   fNPIDdetectors++;
 201 }
 202
 203 //____________________________________________________________
 204 Bool_t AliHFEpid::InitializePID(){
 205   //
 206   // Initializes the PID object
 207   //
 208
 209   TMath::Sort(static_cast<UInt_t>(kNdetectorPID), fDetectorOrder, fSortedOrder, kFALSE);
 210   // Initialize PID Objects
 211   Bool_t status = kTRUE;
 212   for(Int_t idet = 0; idet < kNdetectorPID; idet++){
 213     if(!IsDetectorOn(idet)) continue;
 214     if(fDetectorPID[idet]){
 215       status &= fDetectorPID[idet]->InitializePID();
 216       if(HasMCData() && status) fDetectorPID[idet]->SetHasMCData();
 217     }
 218   }
 219   PrintStatus();
 220   return status;
 221 }
 222
 223 //____________________________________________________________
 224 Bool_t AliHFEpid::IsSelected(AliHFEpidObject *track, AliHFEcontainer *cont, const Char_t *contname, AliHFEpidQAmanager *pidqa){
 225   //
 226   // Select Tracks
 227   //
 228   Bool_t isSelected = kTRUE;
 229   AliDebug(1, Form("Particle used for PID, QA available: %s", pidqa ? "Yes" : "No"));
 230   for(UInt_t idet = 0; idet < fNPIDdetectors; idet++){
 231     AliDebug(2, Form("Using Detector %s\n", SortedDetectorName(idet)));
 232     if(TMath::Abs(fDetectorPID[fSortedOrder[idet]]->IsSelected(track, pidqa)) != 11){
 233       isSelected = kFALSE;
 234       break;
 235     }
 236     AliDebug(2, "Particlae selected by detector");
 237     if(fVarManager && cont){
 238       Char_t reccontname[256];
 239       sprintf(reccontname, "%sReco", contname);
 240       AliDebug(2, Form("Filling container %s", reccontname));
 241       if(fVarManager->IsSignalTrack())
 242         fVarManager->FillContainerStepname(cont, reccontname, SortedDetectorName(idet));
 243       if(HasMCData()){
 244         Char_t mccontname[256];
 245         sprintf(mccontname, "%sMC", contname);
 246         AliDebug(2, Form("MC Information available, Filling container %s", mccontname));
 247         if(fVarManager->IsSignalTrack())
 248           fVarManager->FillContainerStepname(cont, mccontname, SortedDetectorName(idet), kTRUE);
 249       }
 250       // The PID will NOT fill the double counting information
 251     }
 252   }
 253   return isSelected;
 254 }
 255
 256 //____________________________________________________________
 257 void AliHFEpid::SetESDpid(AliESDpid *pid){
 258   //
 259   // Set ESD PID to the Detector PID objects
 260   //
 261   for(Int_t idet = 0; idet < kNdetectorPID; idet++){
 262     if(fDetectorPID[idet]) fDetectorPID[idet]->SetESDpid(pid);
 263   }
 264 }
 265
 266 //____________________________________________________________
 267 void AliHFEpid::SetAODpid(AliAODpidUtil *pid){
 268   //
 269   // Set ESD PID to the Detector PID objects
 270   //
 271   for(Int_t idet = 0; idet < kNdetectorPID; idet++){
 272     if(fDetectorPID[idet]) fDetectorPID[idet]->SetAODpid(pid);
 273   }
 274 }
 275
 276 //____________________________________________________________
 277 void AliHFEpid::ConfigureTPCasymmetric(Double_t pmin, Double_t pmax, Double_t sigmamin, Double_t sigmamax){
 278   //
 279   // 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
 280   //
 281   AliHFEpidTPC *pid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
 282   pid->SetTPCnSigma(3);
 283   pid->SetAsymmetricTPCsigmaCut(pmin, pmax, sigmamin, sigmamax);
 284 }
 285
 286 //____________________________________________________________
 287 void AliHFEpid::ConfigureTPCrejectionSimple(){
 288   //
 289   // TPC alone, symmetric 3 sigma cut and 2 - -100 sigma pion rejection
 290   //
 291   AliHFEpidTPC *pid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
 292   pid->SetTPCnSigma(3);
 293   pid->SetRejectParticle(AliPID::kPion, 0., -100., 10., 1.);
 294 }
 295
 296 //____________________________________________________________
 297 void AliHFEpid::ConfigureTPCrejection(){
 298   //
 299   // Combined TPC-TOF PID, combination is discribed in the funtion MakePidTpcTof
 300   //
 301   if(HasMCData()) printf("Configuring TPC for MC\n");
 302   AliHFEpidTPC *tpcpid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
 303   AliHFEpidTOF *tofpid = dynamic_cast<AliHFEpidTOF *>(fDetectorPID[kTOFpid]);
 304   tpcpid->SetTPCnSigma(2);
 305   tofpid->SetTOFnSigma(3);
 306
 307   //TF1 *upperCut = new TF1("upperCut", "[0] * TMath::Exp([1]*x)", 0, 20);
 308   TF1 *upperCut = new TF1("upperCut", "[0]", 0, 20); // Use constant upper cut
 309   TF1 *lowerCut = new TF1("lowerCut", "[0] * TMath::Exp([1]*x) + [2]", 0, 20);
 310   upperCut->SetParameter(0, 3.);
 311   //upperCut->SetParameter(0, 2.7);
 312   //upperCut->SetParameter(1, -0.4357);
 313
 314   if(HasMCData()) lowerCut->SetParameter(0, -2.5);
 315   else lowerCut->SetParameter(0, -3.7);
 316
 317   lowerCut->SetParameter(1, -0.8);
 318
 319   if(HasMCData()) lowerCut->SetParameter(2, -2.2);
 320   else lowerCut->SetParameter(2, -0.35);
 321
 322   tpcpid->SetUpperSigmaCut(upperCut);
 323   tpcpid->SetLowerSigmaCut(lowerCut);
 324   AddCommonObject(upperCut);
 325   AddCommonObject(lowerCut);
 326 }
 327
 328 //____________________________________________________________
 329 void AliHFEpid::ConfigureTPCstrategyParis(){
 330   //
 331   // TPC alone, symmetric 3 sigma cut and 2 - -100 sigma pion rejection
 332   //
 333   AliHFEpidTPC *pid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
 334   pid->SetTPCnSigma(2);
 335   pid->SetRejectParticle(AliPID::kProton, 0., -3., 10., 3.);
 336   pid->SetRejectParticle(AliPID::kKaon, 0., -3., 10., 3.);
 337 }
 338
 339 //____________________________________________________________
 340 void AliHFEpid::PrintStatus() const {
 341   //
 342   // Print the PID configuration
 343   //
 344   printf("\n%s: Printing configuration\n", GetName());
 345   printf("===============================================\n");
 346   printf("PID Detectors: \n");
 347   Int_t npid = 0;
 348   TString detectors[kNdetectorPID] = {"MC", "ESD", "ITS", "TPC", "TRD", "TOF"};
 349   for(Int_t idet = 0; idet < kNdetectorPID; idet++){
 350     if(IsDetectorOn(idet)){
 351       printf("\t%s\n", detectors[idet].Data());
 352       npid++;
 353     }
 354   }
 355   if(!npid) printf("\tNone\n");
 356   printf("\n");
 357 }