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(AliHFEpidObject *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::ConfigureTPCrejection(const char *lowerCutParam, Double_t * const params, Float_t upperTPCCut, Float_t TOFCut){
 311   //
 312   // Combined TPC-TOF PID
 313   // if no function parameterizaion is given, then the default one (exponential) is chosen
 314   //
 315   if(HasMCData()) AliInfo("Configuring TPC for MC\n");
 316   AliHFEpidTPC *tpcpid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
 317   AliHFEpidTOF *tofpid = dynamic_cast<AliHFEpidTOF *>(fDetectorPID[kTOFpid]);
 318   if(tofpid) tofpid->SetTOFnSigma(TOFCut);
 319
 320   //TF1 *upperCut = new TF1("upperCut", "[0] * TMath::Exp([1]*x)", 0, 20);
 321   TF1 *upperCut = new TF1("upperCut", "[0]", 0, 20); // Use constant upper cut
 322   TF1 *lowerCut = new TF1("lowerCut", lowerCutParam == NULL ? "[0] * TMath::Exp([1]*x) + [2]": lowerCutParam, 0, 20);
 323
 324   upperCut->SetParameter(0, upperTPCCut); // pp
 325
 326 //  upperCut->SetParameter(0, 3.5); //PbPb
 327 //  printf("upper %f \n",upperCut->GetParameter(0));
 328   //upperCut->SetParameter(0, 2.7);
 329   //upperCut->SetParameter(1, -0.4357);
 330   if(params){
 331     for(Int_t ipar = 0; ipar < lowerCut->GetNpar(); ipar++) lowerCut->SetParameter(ipar, params[ipar]);
 332   } else {
 333     // Set default parameterization
 334       if(HasMCData()) lowerCut->SetParameter(0, -2.5);
 335       else lowerCut->SetParameter(0, -4.03);  //pp
 336 //      else lowerCut->SetParameter(0, -3.83);  //PbPb
 337
 338
 339     //else lowerCut->SetParameter(0, -3.71769);
 340     //else lowerCut->SetParameter(0, -3.7);
 341
 342     lowerCut->SetParameter(1, -0.22); // pp
 343 //     lowerCut->SetParameter(1,-0.36 ); // PbPb
 344     //lowerCut->SetParameter(1, -0.40263);
 345     //lowerCut->SetParameter(1, -0.8);
 346
 347     if(HasMCData()) lowerCut->SetParameter(2, -2.2);
 348     else lowerCut->SetParameter(2, 0.92); //pp
 349 //     else lowerCut->SetParameter(2, 0.27); //PbPb
 350    // else lowerCut->SetParameter(2, 0.92); //pp
 351 //     printf("lower0 %f \n",lowerCut->GetParameter(0));
 352 //     printf("lower1 %f \n",lowerCut->GetParameter(1));
 353 //     printf("lower2 %f \n",lowerCut->GetParameter(2));
 354     //else lowerCut->SetParameter(2, 0.267857);
 355     //else lowerCut->SetParameter(2, -0.35);
 356   }
 357
 358
 359   if(tpcpid){
 360     tpcpid->SetTPCnSigma(2);
 361     tpcpid->SetUpperSigmaCut(upperCut);
 362     tpcpid->SetLowerSigmaCut(lowerCut);
 363   }
 364   AddCommonObject(upperCut);
 365   AddCommonObject(lowerCut);
 366 }
 367
 368 //____________________________________________________________
 369 void AliHFEpid::ConfigureTPCstrategyParis(){
 370   //
 371   // TPC alone, symmetric 3 sigma cut and 2 - -100 sigma pion rejection
 372   //
 373   AliHFEpidTPC *pid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
 374   if(pid){
 375     pid->SetTPCnSigma(2);
 376     pid->SetRejectParticle(AliPID::kProton, 0., -3., 10., 3.);
 377     pid->SetRejectParticle(AliPID::kKaon, 0., -3., 10., 3.);
 378   }
 379 }
 380
 381 //____________________________________________________________
 382 void AliHFEpid::PrintStatus() const {
 383   //
 384   // Print the PID configuration
 385   //
 386   printf("\n%s: Printing configuration\n", GetName());
 387   printf("===============================================\n");
 388   printf("PID Detectors: \n");
 389   Int_t npid = 0;
 390   TString detectors[kNdetectorPID] = {"MC", "ESD", "ITS", "TPC", "TRD", "TOF", "EMCAL"};
 391   for(Int_t idet = 0; idet < kNdetectorPID; idet++){
 392     if(IsDetectorOn(idet)){
 393       printf("\t%s\n", detectors[idet].Data());
 394       npid++;
 395     }
 396   }
 397   if(!npid) printf("\tNone\n");
 398   printf("\n");
 399 }