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  *                                                                      *
  17  * PID Steering Class                                                   *
  18  * Interface to the user task                                           *
  19  *                                                                      *
  20  * Authors:                                                             *
  21  *   Markus Fasel <M.Fasel@gsi.de>                                      *
  22  *                                                                      *
  23  ************************************************************************/
  24 #include <TClass.h>
  25 #include <THnSparse.h>
  26 #include <TH3F.h>
  27 #include <TIterator.h>
  28 #include <TList.h>
  29 #include <TObjArray.h>
  30 #include <TObjString.h>
  31 #include <TString.h>
  32
  33 #include "AliESDtrack.h"
  34 #include "AliLog.h"
  35 #include "AliPID.h"
  36
  37 #include "AliHFEpid.h"
  38 #include "AliHFEpidBase.h"
  39 #include "AliHFEpidITS.h"
  40 #include "AliHFEpidTPC.h"
  41 #include "AliHFEpidTRD.h"
  42 #include "AliHFEpidTOF.h"
  43 #include "AliHFEpidMC.h"
  44
  45 ClassImp(AliHFEpid)
  46
  47 //____________________________________________________________
  48 AliHFEpid::AliHFEpid():
  49   fEnabledDetectors(0),
  50   fQAlist(0x0)
  51 {
  52   //
  53   // Default constructor
  54   //
  55   memset(fDetectorPID, 0, sizeof(AliHFEpidBase *) * kNdetectorPID);
  56 }
  57
  58 //____________________________________________________________
  59 AliHFEpid::AliHFEpid(const AliHFEpid &c):
  60   TObject(c),
  61   fEnabledDetectors(c.fEnabledDetectors),
  62   fQAlist(0x0)
  63 {
  64   //
  65   // Copy Constructor
  66   //
  67   memset(fDetectorPID, 0, sizeof(AliHFEpidBase *) * kNdetectorPID);
  68   if(c.fDetectorPID[kMCpid])
  69     fDetectorPID[kMCpid] = new AliHFEpidMC(*(dynamic_cast<AliHFEpidMC *>(c.fDetectorPID[kMCpid])));
  70   if(c.fDetectorPID[kTPCpid])
  71     fDetectorPID[kTPCpid] = new AliHFEpidTPC(*(dynamic_cast<AliHFEpidTPC *>(c.fDetectorPID[kTPCpid])));
  72   if(c.fDetectorPID[kTRDpid])
  73     fDetectorPID[kTRDpid] = new AliHFEpidTRD(*(dynamic_cast<AliHFEpidTRD *>(c.fDetectorPID[kTOFpid])));
  74   if(c.fDetectorPID[kTOFpid])
  75     fDetectorPID[kTOFpid] = new AliHFEpidTOF(*(dynamic_cast<AliHFEpidTOF *>(c.fDetectorPID[kTOFpid])));
  76   if(c.IsQAOn()) SetQAOn();
  77   if(c.HasMCData()) SetHasMCData(kTRUE);
  78   for(Int_t idet = 0; idet < kNdetectorPID; idet++){
  79     if(c.IsQAOn() && fDetectorPID[idet]) fDetectorPID[idet]->SetQAOn(fQAlist);
  80     if(c.HasMCData() && fDetectorPID[idet]) fDetectorPID[idet]->SetHasMCData(kTRUE);
  81   }
  82 }
  83
  84 //____________________________________________________________
  85 AliHFEpid& AliHFEpid::operator=(const AliHFEpid &c){
  86   //
  87   // Assignment operator
  88   //
  89   TObject::operator=(c);
  90
  91   if(this != &c){
  92     fEnabledDetectors = c.fEnabledDetectors;
  93     fQAlist = 0x0;
  94
  95     memset(fDetectorPID, 0, sizeof(AliHFEpidBase *) * kNdetectorPID);
  96     if(c.fDetectorPID[kMCpid])
  97       fDetectorPID[kMCpid] = new AliHFEpidMC(*(dynamic_cast<AliHFEpidMC *>(c.fDetectorPID[kMCpid])));
  98     if(c.fDetectorPID[kTPCpid])
  99       fDetectorPID[kTPCpid] = new AliHFEpidTPC(*(dynamic_cast<AliHFEpidTPC *>(c.fDetectorPID[kTPCpid])));
 100     if(c.fDetectorPID[kTRDpid])
 101       fDetectorPID[kTRDpid] = new AliHFEpidTRD(*(dynamic_cast<AliHFEpidTRD *>(c.fDetectorPID[kTOFpid])));
 102     if(c.fDetectorPID[kTOFpid])
 103       fDetectorPID[kTOFpid] = new AliHFEpidTOF(*(dynamic_cast<AliHFEpidTOF *>(c.fDetectorPID[kTOFpid])));
 104     if(c.IsQAOn()) SetQAOn();
 105     if(c.HasMCData()) SetHasMCData(kTRUE);
 106     for(Int_t idet = 0; idet < kNdetectorPID; idet++){
 107       if(c.IsQAOn() && fDetectorPID[idet]) fDetectorPID[idet]->SetQAOn(fQAlist);
 108       if(c.HasMCData() && fDetectorPID[idet]) fDetectorPID[idet]->SetHasMCData();
 109     }
 110   }
 111   return *this;
 112 }
 113
 114 //____________________________________________________________
 115 AliHFEpid::~AliHFEpid(){
 116   //
 117   // Destructor
 118   //
 119   for(Int_t idet = 0; idet < kNdetectorPID; idet++){
 120     if(fDetectorPID[idet])
 121       delete fDetectorPID[idet];
 122   }
 123   if(fQAlist) delete fQAlist; fQAlist = 0x0;  // Each detector has to care about its Histograms
 124 }
 125
 126 //____________________________________________________________
 127 Bool_t AliHFEpid::InitializePID(TString detectors){
 128   //
 129   // Initializes PID Object:
 130   // + Defines which detectors to use
 131   // + Initializes Detector PID objects
 132   // + Handles QA
 133   //
 134   AliInfo(Form("Doing InitializePID for Detectors %s", detectors.Data()));
 135   fDetectorPID[kMCpid] = new AliHFEpidMC("Monte Carlo PID"); // Always there
 136   fDetectorPID[kITSpid] = new AliHFEpidITS("ITS development PID");  // Development version of the ITS pid, for the moment always there
 137   SETBIT(fEnabledDetectors, kMCpid);
 138   SETBIT(fEnabledDetectors, kITSpid);
 139
 140   TObjArray *detsEnabled = detectors.Tokenize(":");
 141   TIterator *detIterator = detsEnabled->MakeIterator();
 142   TObjString *det = 0x0;
 143   while((det = dynamic_cast<TObjString *>(detIterator->Next()))){
 144     if(det->String().CompareTo("TPC") == 0){
 145       fDetectorPID[kTPCpid] = new AliHFEpidTPC("TPC PID");
 146       SETBIT(fEnabledDetectors, kTPCpid);
 147     } else if(det->String().CompareTo("TRD") == 0){
 148       fDetectorPID[kTRDpid] = new AliHFEpidTRD("TRD PID");
 149       SETBIT(fEnabledDetectors, kTRDpid);
 150     } else if(det->String().CompareTo("TOF") == 0){
 151       fDetectorPID[kTOFpid] = new AliHFEpidTOF("TOF PID");
 152       SETBIT(fEnabledDetectors, kTOFpid);
 153     }
 154     // Here is still space for ESD PID
 155   }
 156   // Initialize PID Objects
 157   Bool_t status = kTRUE;
 158   for(Int_t idet = 0; idet < kNdetectorPID; idet++){
 159     if(fDetectorPID[idet]){
 160       status &= fDetectorPID[idet]->InitializePID();
 161       if(IsQAOn() && status) fDetectorPID[idet]->SetQAOn(fQAlist);
 162       if(HasMCData() && status) fDetectorPID[idet]->SetHasMCData();
 163     }
 164   }
 165   return status;
 166 }
 167
 168 //____________________________________________________________
 169 Bool_t AliHFEpid::IsSelected(AliVParticle *track){
 170   //
 171   // Steers PID decision for single detectors respectively combined
 172   // PID decision
 173   //
 174   if(TString(track->IsA()->GetName()).CompareTo("AliMCparticle") == 0){
 175     return (TMath::Abs(fDetectorPID[kMCpid]->IsSelected(track)) == 11);
 176   }
 177   if(TString(track->IsA()->GetName()).CompareTo("AliESDtrack") == 0){
 178     if(TESTBIT(fEnabledDetectors, kTPCpid)){
 179       if(IsQAOn())
 180         MakePlotsItsTpc(dynamic_cast<AliESDtrack *>(track));  // First fill the QA histograms
 181       if(TESTBIT(fEnabledDetectors, kTOFpid)){
 182         // case TPC-TOF
 183         return MakePidTpcTof(dynamic_cast<AliESDtrack *>(track));
 184       } else if(TESTBIT(fEnabledDetectors, kTRDpid)){
 185         // case TPC-TRD with low level detector Signals
 186         return MakePidTpcTrd(dynamic_cast<AliESDtrack *>(track));
 187       } else
 188         return (TMath::Abs(fDetectorPID[kTPCpid]->IsSelected(track)) ==11);
 189     } else if(TESTBIT(fEnabledDetectors, kTRDpid)){
 190       return (TMath::Abs(fDetectorPID[kTRDpid]->IsSelected(track)) ==11);
 191     } else if(TESTBIT(fEnabledDetectors, kTOFpid)){
 192       return (TMath::Abs(fDetectorPID[kTOFpid]->IsSelected(track)) ==11);
 193     }
 194
 195   }
 196   return kFALSE;
 197 }
 198
 199 //____________________________________________________________
 200 Bool_t AliHFEpid::MakePidTpcTof(AliESDtrack *track){
 201   //
 202   // Combines TPC and TOF PID decision
 203   //
 204   if(fDetectorPID[kTOFpid]->IsSelected(track)) return fDetectorPID[kTPCpid]->IsSelected(track);
 205   return kFALSE;
 206 }
 207
 208 //____________________________________________________________
 209 Bool_t AliHFEpid::MakePidTpcTrd(AliESDtrack *track){
 210   //
 211   // Combination of TPC and TRD PID
 212   // Fills Histograms TPC Signal vs. TRD signal for different
 213   // momentum slices
 214   //
 215   Double_t content[10];
 216   content[0] = -1;
 217   content[1] = track->P();
 218   content[2] = track->GetTPCsignal();
 219   AliHFEpidTRD *trdPid = dynamic_cast<AliHFEpidTRD *>(fDetectorPID[kTRDpid]);
 220   content[3] = trdPid->GetTRDSignalV1(track);
 221   content[4] = trdPid->GetTRDSignalV2(track);
 222   AliDebug(1, Form("Momentum: %f, TRD Signal: Method 1[%f], Method 2[%f]", content[1], content[3], content[4]));
 223   if(IsQAOn()){
 224     if(HasMCData()){
 225       // Fill My Histograms for MC PID
 226       Int_t pdg = TMath::Abs(fDetectorPID[kMCpid]->IsSelected(track));
 227       Int_t pid = -1;
 228       switch(pdg){
 229         case 11:    pid = AliPID::kElectron; break;
 230         case 13:    pid = AliPID::kMuon; break;
 231         case 211:   pid = AliPID::kPion; break;
 232         case 321:   pid = AliPID::kKaon; break;
 233         case 2212:  pid = AliPID::kProton; break;
 234         default:    pid = -1;
 235       };
 236       content[0] = pid;
 237     }
 238     (dynamic_cast<THnSparse *>(fQAlist->At(kTRDSignal)))->Fill(content);
 239   }
 240   return trdPid->IsSelected(track);
 241 }
 242
 243 //____________________________________________________________
 244 void AliHFEpid::MakePlotsItsTpc(AliESDtrack *track){
 245   //
 246   // Make a plot ITS signal - TPC signal for several momentum bins
 247   //
 248   Double_t content[10];
 249   content[0] = -1;
 250   content[1] = track->GetTPCInnerParam() ? track->GetTPCInnerParam()->P() : track->P();
 251   content[2] = (dynamic_cast<AliHFEpidITS *>(fDetectorPID[kITSpid]))->GetITSSignalV1(track);
 252   content[3] = track->GetTPCsignal();
 253   AliDebug(1, Form("Momentum %f, TPC Signal %f, ITS Signal %f", content[1], content[2], content[3]));
 254   if(HasMCData()){
 255     // Fill My Histograms for MC PID
 256     Int_t pdg = TMath::Abs(fDetectorPID[kMCpid]->IsSelected(track));
 257     Int_t pid = -1;
 258     switch(pdg){
 259       case 11:    pid = AliPID::kElectron; break;
 260       case 13:    pid = AliPID::kMuon; break;
 261       case 211:   pid = AliPID::kPion; break;
 262       case 321:   pid = AliPID::kKaon; break;
 263       case 2212:  pid = AliPID::kProton; break;
 264       default:    pid = -1;
 265     };
 266     content[0] = pid;
 267   }
 268   (dynamic_cast<THnSparse *>(fQAlist->At(kITSSignal)))->Fill(content);
 269 }
 270
 271 //____________________________________________________________
 272 void AliHFEpid::SetQAOn(){
 273   //
 274   // Switch on QA
 275   //
 276   SetBit(kIsQAOn, kTRUE);
 277   if(fQAlist) return;
 278   fQAlist = new TList;
 279   fQAlist->SetName("PIDqa");
 280   THnSparseF *histo = NULL;
 281
 282   // Prepare axis for QA histograms
 283   const Int_t kMomentumBins = 41;
 284   const Double_t kPtMin = 0.1;
 285   const Double_t kPtMax = 10.;
 286   Double_t momentumBins[kMomentumBins];
 287   for(Int_t ibin = 0; ibin < kMomentumBins; ibin++)
 288     momentumBins[ibin] = static_cast<Double_t>(TMath::Power(10,TMath::Log10(kPtMin) + (TMath::Log10(kPtMax)-TMath::Log10(kPtMin))/(kMomentumBins-1)*static_cast<Double_t>(ibin)));
 289
 290   // Add Histogram for combined TPC-TRD PID
 291   const Int_t kDimensionsTRDsig = 5;
 292   Int_t kNbinsTRDsig[kDimensionsTRDsig] = {AliPID::kSPECIES + 1, kMomentumBins - 1, 200, 3000, 3000};
 293   Double_t binMinTRDsig[kDimensionsTRDsig] = {-1., 0.1, 0, 0, 0};
 294   Double_t binMaxTRDsig[kDimensionsTRDsig] = {AliPID::kSPECIES, 10., 200., 3000., 3000.};
 295   fQAlist->AddAt((histo = new THnSparseF("fCombTPCTRDpid", "Combined TPC-TRD PID", kDimensionsTRDsig, kNbinsTRDsig, binMinTRDsig, binMaxTRDsig)), kTRDSignal);
 296   histo->GetAxis(1)->Set(kMomentumBins - 1, momentumBins);
 297   histo->GetAxis(0)->SetTitle("Particle Species");
 298   histo->GetAxis(1)->SetTitle("p / GeV/c");
 299   histo->GetAxis(2)->SetTitle("TPC Signal / a.u.");
 300   histo->GetAxis(3)->SetTitle("TRD Signal / a.u.");
 301   histo->GetAxis(4)->SetTitle("TRD Signal / a.u.");
 302
 303   // Add Histogram for combined TPC-ITS PID
 304   const Int_t kDimensionsITSsig = 4;
 305   Int_t kNbinsITSsig[kDimensionsITSsig] = {AliPID::kSPECIES + 1, kMomentumBins - 1, 200, 1000};
 306   Double_t binMinITSsig[kDimensionsITSsig] = {-1., 0.1, 0., 0.};
 307   Double_t binMaxITSsig[kDimensionsITSsig] = {AliPID::kSPECIES, 10., 200., 300.};
 308   fQAlist->AddAt((histo = new THnSparseF("fCombTPCITSpid", "Combined TPC-ITS PID", kDimensionsITSsig, kNbinsITSsig, binMinITSsig, binMaxITSsig)), kITSSignal);
 309   histo->GetAxis(0)->Set(kMomentumBins - 1, momentumBins);
 310   histo->GetAxis(0)->SetTitle("Particle Species");
 311   histo->GetAxis(1)->SetTitle("p / GeV/c");
 312   histo->GetAxis(2)->SetTitle("TPC Signal / a.u.");
 313   histo->GetAxis(3)->SetTitle("ITS Signal / a.u.");
 314 }
 315
 316 //____________________________________________________________
 317 void AliHFEpid::SetMCEvent(AliMCEvent *event){
 318   //
 319   // Set MC Event for the detector PID classes
 320   //
 321   for(Int_t idet = 0; idet < kNdetectorPID; idet++)
 322     if(fDetectorPID[idet]) fDetectorPID[idet]->SetMCEvent(event);
 323 }
 324