[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEpid.cxx

/**************************************************************************
* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
*                                                                        *
* Author: The ALICE Off-line Project.                                    *
* Contributors are mentioned in the code where appropriate.              *
*                                                                        *
* Permission to use, copy, modify and distribute this software and its   *
* documentation strictly for non-commercial purposes is hereby granted   *
* without fee, provided that the above copyright notice appears in all   *
* copies and that both the copyright notice and this permission notice   *
* appear in the supporting documentation. The authors make no claims     *
* about the suitability of this software for any purpose. It is          *
* provided "as is" without express or implied warranty.                  *
**************************************************************************/
/************************************************************************
 *                                                                      *
 * PID Steering Class                                                   *
 * Interface to the user task                                           *
 *                                                                      *
 * Authors:                                                             *
 *   Markus Fasel <M.Fasel@gsi.de>                                      *
 *                                                                      *
 ************************************************************************/
#include <TClass.h>
#include <THnSparse.h>
#include <TH3F.h>
#include <TIterator.h>
#include <TList.h>
#include <TObjArray.h>
#include <TObjString.h>
#include <TString.h>

#include "AliAODTrack.h"
#include "AliESDtrack.h"
#include "AliLog.h"
#include "AliPID.h"

#include "AliHFEpid.h"
#include "AliHFEpidBase.h"
#include "AliHFEpidITS.h"
#include "AliHFEpidTPC.h"
#include "AliHFEpidTRD.h"
#include "AliHFEpidTOF.h"
#include "AliHFEpidMC.h"

ClassImp(AliHFEpid)

//____________________________________________________________
AliHFEpid::AliHFEpid():
  fEnabledDetectors(0),
  fPIDstrategy(0),
  fQAlist(0x0),
  fDebugLevel(0)
{
  //
  // Default constructor
  //
  memset(fDetectorPID, 0, sizeof(AliHFEpidBase *) * kNdetectorPID);
}

//____________________________________________________________
AliHFEpid::AliHFEpid(const AliHFEpid &c):
  TObject(c),
  fEnabledDetectors(c.fEnabledDetectors),
  fPIDstrategy(0),
  fQAlist(0x0),
  fDebugLevel(c.fDebugLevel)
{
  //
  // Copy Constructor
  //
  memset(fDetectorPID, 0, sizeof(AliHFEpidBase *) * kNdetectorPID);
  if(c.fDetectorPID[kMCpid])
    fDetectorPID[kMCpid] = new AliHFEpidMC(*(dynamic_cast<AliHFEpidMC *>(c.fDetectorPID[kMCpid])));
  if(c.fDetectorPID[kTPCpid])
    fDetectorPID[kTPCpid] = new AliHFEpidTPC(*(dynamic_cast<AliHFEpidTPC *>(c.fDetectorPID[kTPCpid])));
  if(c.fDetectorPID[kTRDpid])
    fDetectorPID[kTRDpid] = new AliHFEpidTRD(*(dynamic_cast<AliHFEpidTRD *>(c.fDetectorPID[kTOFpid])));
  if(c.fDetectorPID[kTOFpid])
    fDetectorPID[kTOFpid] = new AliHFEpidTOF(*(dynamic_cast<AliHFEpidTOF *>(c.fDetectorPID[kTOFpid])));
  if(c.IsQAOn()) SetQAOn();
  if(c.HasMCData()) SetHasMCData(kTRUE);
  for(Int_t idet = 0; idet < kNdetectorPID; idet++){
    if(c.IsQAOn() && fDetectorPID[idet]) fDetectorPID[idet]->SetQAOn(fQAlist);
    if(c.HasMCData() && fDetectorPID[idet]) fDetectorPID[idet]->SetHasMCData(kTRUE);
  }
}

//____________________________________________________________
AliHFEpid& AliHFEpid::operator=(const AliHFEpid &c){
  //
  // Assignment operator
  //
  TObject::operator=(c);

  if(this != &c){
    fEnabledDetectors = c.fEnabledDetectors;
    fPIDstrategy = c.fPIDstrategy;
    fQAlist = 0x0;
    fDebugLevel = c.fDebugLevel;
  
    memset(fDetectorPID, 0, sizeof(AliHFEpidBase *) * kNdetectorPID);
    if(c.fDetectorPID[kMCpid])
      fDetectorPID[kMCpid] = new AliHFEpidMC(*(dynamic_cast<AliHFEpidMC *>(c.fDetectorPID[kMCpid])));
    if(c.fDetectorPID[kTPCpid])
      fDetectorPID[kTPCpid] = new AliHFEpidTPC(*(dynamic_cast<AliHFEpidTPC *>(c.fDetectorPID[kTPCpid])));
    if(c.fDetectorPID[kTRDpid])
      fDetectorPID[kTRDpid] = new AliHFEpidTRD(*(dynamic_cast<AliHFEpidTRD *>(c.fDetectorPID[kTOFpid])));
    if(c.fDetectorPID[kTOFpid])
      fDetectorPID[kTOFpid] = new AliHFEpidTOF(*(dynamic_cast<AliHFEpidTOF *>(c.fDetectorPID[kTOFpid])));
    if(c.IsQAOn()) SetQAOn();
    if(c.HasMCData()) SetHasMCData(kTRUE);
    for(Int_t idet = 0; idet < kNdetectorPID; idet++){
      if(c.IsQAOn() && fDetectorPID[idet]) fDetectorPID[idet]->SetQAOn(fQAlist);
      if(c.HasMCData() && fDetectorPID[idet]) fDetectorPID[idet]->SetHasMCData();
    }
  }
  return *this; 
}

//____________________________________________________________
AliHFEpid::~AliHFEpid(){
  //
  // Destructor
  //
  for(Int_t idet = 0; idet < kNdetectorPID; idet++){
    if(fDetectorPID[idet])
      delete fDetectorPID[idet];
  } 
  if(fQAlist) delete fQAlist; fQAlist = 0x0;  // Each detector has to care about its Histograms
}

//____________________________________________________________
Bool_t AliHFEpid::InitializePID(TString arg){
  //
  // Initializes PID Object:
  // + Defines which detectors to use
  // + Initializes Detector PID objects
  // + Handles QA
  //
  fDetectorPID[kMCpid] = new AliHFEpidMC("Monte Carlo PID"); // Always there
  SETBIT(fEnabledDetectors, kMCpid);
  // Initialize detector PIDs according to PID Strategies
  if(arg.BeginsWith("Strategy")){
    arg.ReplaceAll("Strategy", "");
    fPIDstrategy = arg.Atoi();
    AliInfo(Form("PID Strategy %d enabled", fPIDstrategy));
    Int_t strategyStatus = kTRUE;
    switch(fPIDstrategy){
      case 1: InitStrategy1(); break;
      case 2: InitStrategy2(); break;
      case 3: InitStrategy3(); break;
      case 4: InitStrategy4(); break;
      case 5: InitStrategy5(); break;
      default: strategyStatus = kFALSE;
    }
    return strategyStatus;
  }
  // No Strategy defined, Initialize according to detectors specified
  AliInfo(Form("Doing InitializePID for Detectors %s end", arg.Data()));
  fDetectorPID[kITSpid] = new AliHFEpidITS("ITS development PID");  // Development version of the ITS pid, for the moment always there
  SETBIT(fEnabledDetectors, kITSpid);
  
  TObjArray *detsEnabled = arg.Tokenize(":");
  TIterator *detIterator = detsEnabled->MakeIterator();
  TObjString *det = 0x0;
  while((det = dynamic_cast<TObjString *>(detIterator->Next()))){
    if(det->String().CompareTo("TPC") == 0){
      AliInfo("Doing TPC PID");
      fDetectorPID[kTPCpid] = new AliHFEpidTPC("TPC PID");
      SETBIT(fEnabledDetectors, kTPCpid);
    } else if(det->String().CompareTo("TRD") == 0){
      fDetectorPID[kTRDpid] = new AliHFEpidTRD("TRD PID");
      SETBIT(fEnabledDetectors, kTRDpid);
    } else if(det->String().CompareTo("TOF") == 0){
      fDetectorPID[kTOFpid] = new AliHFEpidTOF("TOF PID");
      SETBIT(fEnabledDetectors, kTOFpid);
    }
    // Here is still space for ESD PID
  }
  // Initialize PID Objects
  Bool_t status = kTRUE;
  for(Int_t idet = 0; idet < kNdetectorPID; idet++){
    if(fDetectorPID[idet]){ 
      status &= fDetectorPID[idet]->InitializePID();
      if(IsQAOn() && status) fDetectorPID[idet]->SetQAOn(fQAlist);
      if(HasMCData() && status) fDetectorPID[idet]->SetHasMCData();
    }
  }
  return status;
}

//____________________________________________________________
Bool_t AliHFEpid::IsSelected(AliHFEpidObject *track){
  //
  // Steers PID decision for single detectors respectively combined
  // PID decision
  //
  if(!track->fRecTrack){
    // MC Event
    return (TMath::Abs(fDetectorPID[kMCpid]->IsSelected(track)) == 11);
  }
  if(fPIDstrategy > 0 && fPIDstrategy < 6){
    Int_t pid = 0;
    switch(fPIDstrategy){
      case 1: pid = IdentifyStrategy1(track); break;
      case 2: pid = IdentifyStrategy2(track); break;
      case 3: pid = IdentifyStrategy3(track); break;
      case 4: pid = IdentifyStrategy4(track); break;
      case 5: pid = IdentifyStrategy5(track); break;
      default: break;
    }
    return pid;
  }
  if(TESTBIT(fEnabledDetectors, kTPCpid)){
    if(IsQAOn() && fDebugLevel > 1){ 
      AliInfo("Filling QA plots");
      MakePlotsItsTpc(track);  // First fill the QA histograms
    }
    if(TESTBIT(fEnabledDetectors, kTOFpid)){
      // case TPC-TOF
      return MakePidTpcTof(track);
    } else if(TESTBIT(fEnabledDetectors, kTRDpid)){
      // case TPC-TRD with low level detector Signals
      return MakePidTpcTrd(track);
    } else
      return (TMath::Abs(fDetectorPID[kTPCpid]->IsSelected(track)) ==11);
  } else if(TESTBIT(fEnabledDetectors, kTRDpid)){
    return (TMath::Abs(fDetectorPID[kTRDpid]->IsSelected(track)) ==11);
  } else if(TESTBIT(fEnabledDetectors, kTOFpid)){
    return (TMath::Abs(fDetectorPID[kTOFpid]->IsSelected(track)) ==11);
  }
  
  return kFALSE;
}

//____________________________________________________________
Bool_t AliHFEpid::MakePidTpcTof(AliHFEpidObject *track){
  //
  // Combines TPC and TOF PID decision
  //
  if(fDetectorPID[kTOFpid]->IsSelected(track)) return fDetectorPID[kTPCpid]->IsSelected(track);
  return kFALSE;
}

//____________________________________________________________
Bool_t AliHFEpid::MakePidTpcTrd(AliHFEpidObject *track){
  //
  // Combination of TPC and TRD PID
  // Fills Histograms TPC Signal vs. TRD signal for different
  // momentum slices
  //
  if(track->fAnalysisType != AliHFEpidObject::kESDanalysis) return kFALSE; //AOD based detector PID combination not yet implemented
  AliDebug(1, "Analysis Type OK, do PID");
  AliESDtrack *esdTrack = dynamic_cast<AliESDtrack *>(track->fRecTrack);
  AliHFEpidTRD *trdPid = dynamic_cast<AliHFEpidTRD *>(fDetectorPID[kTRDpid]);
  Int_t pdg = TMath::Abs(fDetectorPID[kMCpid]->IsSelected(track));
  Int_t pid = -1;
  switch(pdg){
    case 11:    pid = AliPID::kElectron; break;
    case 13:    pid = AliPID::kMuon; break;
    case 211:   pid = AliPID::kPion; break;
    case 321:   pid = AliPID::kKaon; break;
    case 2212:  pid = AliPID::kProton; break;
    default:    pid = -1;
  };
  Double_t content[10];
  content[0] = pid;
  content[1] = esdTrack->P();
  content[2] = esdTrack->GetTPCsignal();
  content[3] = trdPid->GetTRDSignalV1(esdTrack, pid);
  if(IsQAOn() && fDebugLevel > 1){
    content[4] = trdPid->GetTRDSignalV2(esdTrack, pid);
    AliDebug(1, Form("Momentum: %f, TRD Signal: Method 1[%f], Method 2[%f]", content[1], content[3], content[4]));
    (dynamic_cast<THnSparseF *>(fQAlist->At(kTRDSignal)))->Fill(content);
  }
  if(content[1] > 2){ // perform combined
    AliDebug(1, "Momentum bigger 2 GeV/c, doing combined PID"); 
    if(content[2] > 65 && content[3] > 500) return kTRUE;
    else return kFALSE;
  }
  else {
    AliDebug(1, "Momentum smaller 2GeV/c, doing TPC alone PID");
    return fDetectorPID[kTPCpid]->IsSelected(track) == 11;
  }
}

//____________________________________________________________
void AliHFEpid::MakePlotsItsTpc(AliHFEpidObject *track){
  //
  // Make a plot ITS signal - TPC signal for several momentum bins
  //
  if(track->fAnalysisType != AliHFEpidObject::kESDanalysis) return; //AOD based detector PID combination not yet implemented
  AliESDtrack * esdTrack = dynamic_cast<AliESDtrack *>(track->fRecTrack);
   // Fill My Histograms for MC PID
  Int_t pdg = TMath::Abs(fDetectorPID[kMCpid]->IsSelected(track));
  Int_t pid = -1;
  switch(pdg){
    case 11:    pid = AliPID::kElectron; break;
    case 13:    pid = AliPID::kMuon; break;
    case 211:   pid = AliPID::kPion; break;
    case 321:   pid = AliPID::kKaon; break;
    case 2212:  pid = AliPID::kProton; break;
    default:    pid = -1;
  };
  if(IsQAOn() && fDebugLevel > 0){
    Double_t content[10];
    content[0] = pid;
    content[1] = esdTrack->GetTPCInnerParam() ? esdTrack->GetTPCInnerParam()->P() : esdTrack->P();
    content[2] = (dynamic_cast<AliHFEpidITS *>(fDetectorPID[kITSpid]))->GetITSSignalV1(track->fRecTrack, pid);
    content[3] = esdTrack->GetTPCsignal();
    AliDebug(1, Form("Momentum %f, TPC Signal %f, ITS Signal %f", content[1], content[2], content[3]));
    (dynamic_cast<THnSparseF *>(fQAlist->At(kITSSignal)))->Fill(content);
  }
}

//____________________________________________________________
void AliHFEpid::SetQAOn(){
  //
  // Switch on QA
  //
  SetBit(kIsQAOn, kTRUE);
  AliInfo("QA switched on");
  if(fQAlist) return;
  fQAlist = new TList;
  fQAlist->SetName("PIDqa");
  THnSparseF *histo = NULL;

  // Prepare axis for QA histograms
  const Int_t kMomentumBins = 41;
  const Double_t kPtMin = 0.1;
  const Double_t kPtMax = 10.;
  Double_t momentumBins[kMomentumBins];
  for(Int_t ibin = 0; ibin < kMomentumBins; ibin++)
    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)));

  // Add Histogram for combined TPC-TRD PID
  if(fDebugLevel > 1){
    AliDebug(1, "Adding histogram for ITS-TPC investigation");
    const Int_t kDimensionsTRDsig = 5;
    Int_t kNbinsTRDsig[kDimensionsTRDsig] = {AliPID::kSPECIES + 1, kMomentumBins - 1, 200, 3000, 3000};
    Double_t binMinTRDsig[kDimensionsTRDsig] = {-1., 0.1, 0, 0, 0};
    Double_t binMaxTRDsig[kDimensionsTRDsig] = {AliPID::kSPECIES, 10., 200., 3000., 3000.};
    fQAlist->AddAt((histo = new THnSparseF("fCombTPCTRDpid", "Combined TPC-TRD PID", kDimensionsTRDsig, kNbinsTRDsig, binMinTRDsig, binMaxTRDsig)), kTRDSignal);
    histo->GetAxis(1)->Set(kMomentumBins - 1, momentumBins);
    histo->GetAxis(0)->SetTitle("Particle Species");
    histo->GetAxis(1)->SetTitle("p / GeV/c");
    histo->GetAxis(2)->SetTitle("TPC Signal / a.u.");
    histo->GetAxis(3)->SetTitle("TRD Signal / a.u.");
    histo->GetAxis(4)->SetTitle("TRD Signal / a.u.");
  }

  // Add Histogram for combined TPC-ITS PID
  if(fDebugLevel > 0){
    AliDebug(1, "Adding histogram for TPC-TRD investigation");
    const Int_t kDimensionsITSsig = 4;
    Int_t kNbinsITSsig[kDimensionsITSsig] = {AliPID::kSPECIES + 1, kMomentumBins - 1, 300, 3000};
    Double_t binMinITSsig[kDimensionsITSsig] = {-1., 0.1, 0., 0.};
    Double_t binMaxITSsig[kDimensionsITSsig] = {AliPID::kSPECIES, 10., 300., 300.};
    fQAlist->AddAt((histo = new THnSparseF("fCombTPCITSpid", "Combined TPC-ITS PID", kDimensionsITSsig, kNbinsITSsig, binMinITSsig, binMaxITSsig)), kITSSignal);
    histo->GetAxis(1)->Set(kMomentumBins - 1, momentumBins);
    histo->GetAxis(0)->SetTitle("Particle Species");
    histo->GetAxis(1)->SetTitle("p / GeV/c");
    histo->GetAxis(2)->SetTitle("ITS Signal / a.u.");
    histo->GetAxis(3)->SetTitle("TPC Signal / a.u.");
  }
}

//____________________________________________________________
void AliHFEpid::InitStrategy1(){
  //
  // TPC alone, 3-sigma cut
  //
  AliHFEpidTPC *pid = new AliHFEpidTPC("strat1TPCpid");
  pid->SetTPCnSigma(3);
  Bool_t status = pid->InitializePID();
  if(IsQAOn() && status) pid->SetQAOn(fQAlist);
  if(HasMCData() && status) pid->SetHasMCData();
  fDetectorPID[kTPCpid] = pid;
}

//____________________________________________________________
void AliHFEpid::InitStrategy2(){
  //
  // 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
  //
  AliHFEpidTPC *pid = new AliHFEpidTPC("strat2TPCpid");
  pid->SetTPCnSigma(3);
  pid->SetAsymmetricTPCsigmaCut(2., 10., 0., 4.);
  Bool_t status = pid->InitializePID();
  if(IsQAOn() && status) pid->SetQAOn(fQAlist);
  if(HasMCData() && status) pid->SetHasMCData();
  fDetectorPID[kTPCpid] = pid;
}

//____________________________________________________________
void AliHFEpid::InitStrategy3(){
  //
  // TPC alone, symmetric 3 sigma cut and 2 - -100 sigma pion rejection
  //   
  AliHFEpidTPC *pid = new AliHFEpidTPC("strat3TPCpid");
  pid->SetTPCnSigma(3);
  pid->SetRejectParticle(AliPID::kPion, 0., -100., 10., 3.);
  Bool_t status = pid->InitializePID();
  if(IsQAOn() && status) pid->SetQAOn(fQAlist);
  if(HasMCData() && status) pid->SetHasMCData();
  fDetectorPID[kTPCpid] = pid;
}

//____________________________________________________________
void AliHFEpid::InitStrategy4(){
  //
  // TPC and TRD combined, TPC 3 sigma cut and TRD NN 90% el efficiency level above 2 GeV/c
  //
  InitStrategy1();
  AliHFEpidTRD *trdpid = new AliHFEpidTRD("strat4TRDpid");
  Bool_t status = trdpid->InitializePID();
  if(IsQAOn() && status) trdpid->SetQAOn(fQAlist);
  if(HasMCData() && status) trdpid->SetHasMCData();
  fDetectorPID[kTRDpid] = trdpid;
}

//____________________________________________________________
void AliHFEpid::InitStrategy5(){
  //
  // TPC and TRD combined, TPC 3 sigma cut and TRD NN 90% el efficiency level above 2 GeV/c
  //
  InitStrategy1();
  AliHFEpidTRD *trdpid = new AliHFEpidTRD("strat5TRDpid");
  Bool_t status = trdpid->InitializePID();
  if(IsQAOn() && status) trdpid->SetQAOn(fQAlist);
  if(HasMCData() && status) trdpid->SetHasMCData();
  fDetectorPID[kTRDpid] = trdpid;
}

//____________________________________________________________
Bool_t AliHFEpid::IdentifyStrategy1(AliHFEpidObject *track){
  return TMath::Abs(fDetectorPID[kTPCpid]->IsSelected(track)) == 11;
}

//____________________________________________________________
Bool_t AliHFEpid::IdentifyStrategy2(AliHFEpidObject *track){
  return TMath::Abs(fDetectorPID[kTPCpid]->IsSelected(track)) == 11;
}

//____________________________________________________________
Bool_t AliHFEpid::IdentifyStrategy3(AliHFEpidObject *track){
  return TMath::Abs(fDetectorPID[kTPCpid]->IsSelected(track)) == 11;
}

//____________________________________________________________
Bool_t AliHFEpid::IdentifyStrategy4(AliHFEpidObject *track){
  Int_t trdpid = TMath::Abs(fDetectorPID[kTRDpid]->IsSelected(track));
  return (trdpid == 0 || trdpid == 11) && (TMath::Abs(fDetectorPID[kTPCpid]->IsSelected(track)) == 11);
}

//____________________________________________________________
Bool_t AliHFEpid::IdentifyStrategy5(AliHFEpidObject *track){
  return MakePidTpcTrd(track);
}
Commit	Line	Data
809a4336	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>
75d81601	25	#include <THnSparse.h>
75d81601	26	#include <TH3F.h>
809a4336	27	#include <TIterator.h>
	28	#include <TList.h>
	29	#include <TObjArray.h>
	30	#include <TObjString.h>
	31	#include <TString.h>
	32
722347d8	33	#include "AliAODTrack.h"
809a4336	34	#include "AliESDtrack.h"
75d81601	35	#include "AliLog.h"
75d81601	36	#include "AliPID.h"
809a4336	37
	38	#include "AliHFEpid.h"
	39	#include "AliHFEpidBase.h"
75d81601	40	#include "AliHFEpidITS.h"
809a4336	41	#include "AliHFEpidTPC.h"
	42	#include "AliHFEpidTRD.h"
	43	#include "AliHFEpidTOF.h"
	44	#include "AliHFEpidMC.h"
	45
	46	ClassImp(AliHFEpid)
	47
	48	//____________________________________________________________
	49	AliHFEpid::AliHFEpid():
	50	fEnabledDetectors(0),
0792aa82	51	fPIDstrategy(0),
78ea5ef4	52	fQAlist(0x0),
78ea5ef4	53	fDebugLevel(0)
809a4336	54	{
	55	//
	56	// Default constructor
	57	//
	58	memset(fDetectorPID, 0, sizeof(AliHFEpidBase ) kNdetectorPID);
	59	}
	60
	61	//____________________________________________________________
	62	AliHFEpid::AliHFEpid(const AliHFEpid &c):
	63	TObject(c),
	64	fEnabledDetectors(c.fEnabledDetectors),
0792aa82	65	fPIDstrategy(0),
78ea5ef4	66	fQAlist(0x0),
78ea5ef4	67	fDebugLevel(c.fDebugLevel)
809a4336	68	{
	69	//
	70	// Copy Constructor
	71	//
	72	memset(fDetectorPID, 0, sizeof(AliHFEpidBase ) kNdetectorPID);
	73	if(c.fDetectorPID[kMCpid])
	74	fDetectorPID[kMCpid] = new AliHFEpidMC((dynamic_cast<AliHFEpidMC >(c.fDetectorPID[kMCpid])));
	75	if(c.fDetectorPID[kTPCpid])
	76	fDetectorPID[kTPCpid] = new AliHFEpidTPC((dynamic_cast<AliHFEpidTPC >(c.fDetectorPID[kTPCpid])));
	77	if(c.fDetectorPID[kTRDpid])
	78	fDetectorPID[kTRDpid] = new AliHFEpidTRD((dynamic_cast<AliHFEpidTRD >(c.fDetectorPID[kTOFpid])));
	79	if(c.fDetectorPID[kTOFpid])
	80	fDetectorPID[kTOFpid] = new AliHFEpidTOF((dynamic_cast<AliHFEpidTOF >(c.fDetectorPID[kTOFpid])));
	81	if(c.IsQAOn()) SetQAOn();
	82	if(c.HasMCData()) SetHasMCData(kTRUE);
	83	for(Int_t idet = 0; idet < kNdetectorPID; idet++){
	84	if(c.IsQAOn() && fDetectorPID[idet]) fDetectorPID[idet]->SetQAOn(fQAlist);
	85	if(c.HasMCData() && fDetectorPID[idet]) fDetectorPID[idet]->SetHasMCData(kTRUE);
	86	}
	87	}
	88
	89	//____________________________________________________________
	90	AliHFEpid& AliHFEpid::operator=(const AliHFEpid &c){
	91	//
	92	// Assignment operator
	93	//
	94	TObject::operator=(c);
	95
	96	if(this != &c){
	97	fEnabledDetectors = c.fEnabledDetectors;
0792aa82	98	fPIDstrategy = c.fPIDstrategy;
809a4336	99	fQAlist = 0x0;
78ea5ef4	100	fDebugLevel = c.fDebugLevel;
809a4336	101
	102	memset(fDetectorPID, 0, sizeof(AliHFEpidBase ) kNdetectorPID);
	103	if(c.fDetectorPID[kMCpid])
	104	fDetectorPID[kMCpid] = new AliHFEpidMC((dynamic_cast<AliHFEpidMC >(c.fDetectorPID[kMCpid])));
	105	if(c.fDetectorPID[kTPCpid])
	106	fDetectorPID[kTPCpid] = new AliHFEpidTPC((dynamic_cast<AliHFEpidTPC >(c.fDetectorPID[kTPCpid])));
	107	if(c.fDetectorPID[kTRDpid])
	108	fDetectorPID[kTRDpid] = new AliHFEpidTRD((dynamic_cast<AliHFEpidTRD >(c.fDetectorPID[kTOFpid])));
	109	if(c.fDetectorPID[kTOFpid])
	110	fDetectorPID[kTOFpid] = new AliHFEpidTOF((dynamic_cast<AliHFEpidTOF >(c.fDetectorPID[kTOFpid])));
	111	if(c.IsQAOn()) SetQAOn();
	112	if(c.HasMCData()) SetHasMCData(kTRUE);
	113	for(Int_t idet = 0; idet < kNdetectorPID; idet++){
	114	if(c.IsQAOn() && fDetectorPID[idet]) fDetectorPID[idet]->SetQAOn(fQAlist);
	115	if(c.HasMCData() && fDetectorPID[idet]) fDetectorPID[idet]->SetHasMCData();
	116	}
	117	}
	118	return *this;
	119	}
	120
	121	//____________________________________________________________
	122	AliHFEpid::~AliHFEpid(){
	123	//
	124	// Destructor
	125	//
809a4336	126	for(Int_t idet = 0; idet < kNdetectorPID; idet++){
75d81601	127	if(fDetectorPID[idet])
	128	delete fDetectorPID[idet];
	129	}
	130	if(fQAlist) delete fQAlist; fQAlist = 0x0; // Each detector has to care about its Histograms
809a4336	131	}
	132
	133	//____________________________________________________________
0792aa82	134	Bool_t AliHFEpid::InitializePID(TString arg){
809a4336	135	//
	136	// Initializes PID Object:
	137	// + Defines which detectors to use
	138	// + Initializes Detector PID objects
	139	// + Handles QA
	140	//
	141	fDetectorPID[kMCpid] = new AliHFEpidMC("Monte Carlo PID"); // Always there
	142	SETBIT(fEnabledDetectors, kMCpid);
0792aa82	143	// Initialize detector PIDs according to PID Strategies
	144	if(arg.BeginsWith("Strategy")){
	145	arg.ReplaceAll("Strategy", "");
	146	fPIDstrategy = arg.Atoi();
	147	AliInfo(Form("PID Strategy %d enabled", fPIDstrategy));
	148	Int_t strategyStatus = kTRUE;
	149	switch(fPIDstrategy){
	150	case 1: InitStrategy1(); break;
	151	case 2: InitStrategy2(); break;
	152	case 3: InitStrategy3(); break;
	153	case 4: InitStrategy4(); break;
	154	case 5: InitStrategy5(); break;
	155	default: strategyStatus = kFALSE;
	156	}
	157	return strategyStatus;
	158	}
	159	// No Strategy defined, Initialize according to detectors specified
	160	AliInfo(Form("Doing InitializePID for Detectors %s end", arg.Data()));
	161	fDetectorPID[kITSpid] = new AliHFEpidITS("ITS development PID"); // Development version of the ITS pid, for the moment always there
75d81601	162	SETBIT(fEnabledDetectors, kITSpid);
809a4336	163
0792aa82	164	TObjArray *detsEnabled = arg.Tokenize(":");
809a4336	165	TIterator *detIterator = detsEnabled->MakeIterator();
	166	TObjString *det = 0x0;
	167	while((det = dynamic_cast<TObjString *>(detIterator->Next()))){
	168	if(det->String().CompareTo("TPC") == 0){
78ea5ef4	169	AliInfo("Doing TPC PID");
809a4336	170	fDetectorPID[kTPCpid] = new AliHFEpidTPC("TPC PID");
	171	SETBIT(fEnabledDetectors, kTPCpid);
	172	} else if(det->String().CompareTo("TRD") == 0){
	173	fDetectorPID[kTRDpid] = new AliHFEpidTRD("TRD PID");
	174	SETBIT(fEnabledDetectors, kTRDpid);
	175	} else if(det->String().CompareTo("TOF") == 0){
	176	fDetectorPID[kTOFpid] = new AliHFEpidTOF("TOF PID");
	177	SETBIT(fEnabledDetectors, kTOFpid);
	178	}
	179	// Here is still space for ESD PID
	180	}
	181	// Initialize PID Objects
	182	Bool_t status = kTRUE;
	183	for(Int_t idet = 0; idet < kNdetectorPID; idet++){
	184	if(fDetectorPID[idet]){
	185	status &= fDetectorPID[idet]->InitializePID();
	186	if(IsQAOn() && status) fDetectorPID[idet]->SetQAOn(fQAlist);
	187	if(HasMCData() && status) fDetectorPID[idet]->SetHasMCData();
	188	}
	189	}
	190	return status;
	191	}
	192
	193	//____________________________________________________________
722347d8	194	Bool_t AliHFEpid::IsSelected(AliHFEpidObject *track){
809a4336	195	//
	196	// Steers PID decision for single detectors respectively combined
	197	// PID decision
	198	//
722347d8	199	if(!track->fRecTrack){
722347d8	200	// MC Event
809a4336	201	return (TMath::Abs(fDetectorPID[kMCpid]->IsSelected(track)) == 11);
809a4336	202	}
0792aa82	203	if(fPIDstrategy > 0 && fPIDstrategy < 6){
	204	Int_t pid = 0;
	205	switch(fPIDstrategy){
	206	case 1: pid = IdentifyStrategy1(track); break;
	207	case 2: pid = IdentifyStrategy2(track); break;
	208	case 3: pid = IdentifyStrategy3(track); break;
	209	case 4: pid = IdentifyStrategy4(track); break;
	210	case 5: pid = IdentifyStrategy5(track); break;
	211	default: break;
	212	}
	213	return pid;
	214	}
722347d8	215	if(TESTBIT(fEnabledDetectors, kTPCpid)){
	216	if(IsQAOn() && fDebugLevel > 1){
	217	AliInfo("Filling QA plots");
	218	MakePlotsItsTpc(track); // First fill the QA histograms
809a4336	219	}
722347d8	220	if(TESTBIT(fEnabledDetectors, kTOFpid)){
	221	// case TPC-TOF
	222	return MakePidTpcTof(track);
	223	} else if(TESTBIT(fEnabledDetectors, kTRDpid)){
	224	// case TPC-TRD with low level detector Signals
	225	return MakePidTpcTrd(track);
	226	} else
	227	return (TMath::Abs(fDetectorPID[kTPCpid]->IsSelected(track)) ==11);
	228	} else if(TESTBIT(fEnabledDetectors, kTRDpid)){
	229	return (TMath::Abs(fDetectorPID[kTRDpid]->IsSelected(track)) ==11);
	230	} else if(TESTBIT(fEnabledDetectors, kTOFpid)){
	231	return (TMath::Abs(fDetectorPID[kTOFpid]->IsSelected(track)) ==11);
809a4336	232	}
722347d8	233
809a4336	234	return kFALSE;
	235	}
	236
	237	//____________________________________________________________
722347d8	238	Bool_t AliHFEpid::MakePidTpcTof(AliHFEpidObject *track){
809a4336	239	//
	240	// Combines TPC and TOF PID decision
	241	//
	242	if(fDetectorPID[kTOFpid]->IsSelected(track)) return fDetectorPID[kTPCpid]->IsSelected(track);
	243	return kFALSE;
	244	}
	245
75d81601	246	//____________________________________________________________
722347d8	247	Bool_t AliHFEpid::MakePidTpcTrd(AliHFEpidObject *track){
75d81601	248	//
	249	// Combination of TPC and TRD PID
	250	// Fills Histograms TPC Signal vs. TRD signal for different
	251	// momentum slices
	252	//
722347d8	253	if(track->fAnalysisType != AliHFEpidObject::kESDanalysis) return kFALSE; //AOD based detector PID combination not yet implemented
0792aa82	254	AliDebug(1, "Analysis Type OK, do PID");
722347d8	255	AliESDtrack esdTrack = dynamic_cast<AliESDtrack >(track->fRecTrack);
75d81601	256	AliHFEpidTRD trdPid = dynamic_cast<AliHFEpidTRD >(fDetectorPID[kTRDpid]);
722347d8	257	Int_t pdg = TMath::Abs(fDetectorPID[kMCpid]->IsSelected(track));
	258	Int_t pid = -1;
	259	switch(pdg){
	260	case 11: pid = AliPID::kElectron; break;
	261	case 13: pid = AliPID::kMuon; break;
	262	case 211: pid = AliPID::kPion; break;
	263	case 321: pid = AliPID::kKaon; break;
	264	case 2212: pid = AliPID::kProton; break;
	265	default: pid = -1;
	266	};
0792aa82	267	Double_t content[10];
	268	content[0] = pid;
	269	content[1] = esdTrack->P();
	270	content[2] = esdTrack->GetTPCsignal();
	271	content[3] = trdPid->GetTRDSignalV1(esdTrack, pid);
722347d8	272	if(IsQAOn() && fDebugLevel > 1){
722347d8	273	content[4] = trdPid->GetTRDSignalV2(esdTrack, pid);
	274	AliDebug(1, Form("Momentum: %f, TRD Signal: Method 1[%f], Method 2[%f]", content[1], content[3], content[4]));
	275	(dynamic_cast<THnSparseF *>(fQAlist->At(kTRDSignal)))->Fill(content);
75d81601	276	}
0792aa82	277	if(content[1] > 2){ // perform combined
	278	AliDebug(1, "Momentum bigger 2 GeV/c, doing combined PID");
	279	if(content[2] > 65 && content[3] > 500) return kTRUE;
	280	else return kFALSE;
	281	}
	282	else {
	283	AliDebug(1, "Momentum smaller 2GeV/c, doing TPC alone PID");
	284	return fDetectorPID[kTPCpid]->IsSelected(track) == 11;
	285	}
75d81601	286	}
	287
	288	//____________________________________________________________
722347d8	289	void AliHFEpid::MakePlotsItsTpc(AliHFEpidObject *track){
75d81601	290	//
	291	// Make a plot ITS signal - TPC signal for several momentum bins
	292	//
722347d8	293	if(track->fAnalysisType != AliHFEpidObject::kESDanalysis) return; //AOD based detector PID combination not yet implemented
	294	AliESDtrack * esdTrack = dynamic_cast<AliESDtrack *>(track->fRecTrack);
	295	// Fill My Histograms for MC PID
	296	Int_t pdg = TMath::Abs(fDetectorPID[kMCpid]->IsSelected(track));
	297	Int_t pid = -1;
	298	switch(pdg){
	299	case 11: pid = AliPID::kElectron; break;
	300	case 13: pid = AliPID::kMuon; break;
	301	case 211: pid = AliPID::kPion; break;
	302	case 321: pid = AliPID::kKaon; break;
	303	case 2212: pid = AliPID::kProton; break;
	304	default: pid = -1;
	305	};
	306	if(IsQAOn() && fDebugLevel > 0){
	307	Double_t content[10];
75d81601	308	content[0] = pid;
722347d8	309	content[1] = esdTrack->GetTPCInnerParam() ? esdTrack->GetTPCInnerParam()->P() : esdTrack->P();
	310	content[2] = (dynamic_cast<AliHFEpidITS *>(fDetectorPID[kITSpid]))->GetITSSignalV1(track->fRecTrack, pid);
	311	content[3] = esdTrack->GetTPCsignal();
	312	AliDebug(1, Form("Momentum %f, TPC Signal %f, ITS Signal %f", content[1], content[2], content[3]));
	313	(dynamic_cast<THnSparseF *>(fQAlist->At(kITSSignal)))->Fill(content);
75d81601	314	}
75d81601	315	}
75d81601	316
809a4336	317	//____________________________________________________________
	318	void AliHFEpid::SetQAOn(){
	319	//
	320	// Switch on QA
	321	//
75d81601	322	SetBit(kIsQAOn, kTRUE);
722347d8	323	AliInfo("QA switched on");
75d81601	324	if(fQAlist) return;
809a4336	325	fQAlist = new TList;
809a4336	326	fQAlist->SetName("PIDqa");
75d81601	327	THnSparseF *histo = NULL;
	328
	329	// Prepare axis for QA histograms
	330	const Int_t kMomentumBins = 41;
	331	const Double_t kPtMin = 0.1;
	332	const Double_t kPtMax = 10.;
	333	Double_t momentumBins[kMomentumBins];
	334	for(Int_t ibin = 0; ibin < kMomentumBins; ibin++)
	335	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)));
	336
	337	// Add Histogram for combined TPC-TRD PID
78ea5ef4	338	if(fDebugLevel > 1){
722347d8	339	AliDebug(1, "Adding histogram for ITS-TPC investigation");
78ea5ef4	340	const Int_t kDimensionsTRDsig = 5;
	341	Int_t kNbinsTRDsig[kDimensionsTRDsig] = {AliPID::kSPECIES + 1, kMomentumBins - 1, 200, 3000, 3000};
	342	Double_t binMinTRDsig[kDimensionsTRDsig] = {-1., 0.1, 0, 0, 0};
	343	Double_t binMaxTRDsig[kDimensionsTRDsig] = {AliPID::kSPECIES, 10., 200., 3000., 3000.};
	344	fQAlist->AddAt((histo = new THnSparseF("fCombTPCTRDpid", "Combined TPC-TRD PID", kDimensionsTRDsig, kNbinsTRDsig, binMinTRDsig, binMaxTRDsig)), kTRDSignal);
	345	histo->GetAxis(1)->Set(kMomentumBins - 1, momentumBins);
	346	histo->GetAxis(0)->SetTitle("Particle Species");
	347	histo->GetAxis(1)->SetTitle("p / GeV/c");
	348	histo->GetAxis(2)->SetTitle("TPC Signal / a.u.");
	349	histo->GetAxis(3)->SetTitle("TRD Signal / a.u.");
	350	histo->GetAxis(4)->SetTitle("TRD Signal / a.u.");
	351	}
75d81601	352
75d81601	353	// Add Histogram for combined TPC-ITS PID
78ea5ef4	354	if(fDebugLevel > 0){
722347d8	355	AliDebug(1, "Adding histogram for TPC-TRD investigation");
78ea5ef4	356	const Int_t kDimensionsITSsig = 4;
	357	Int_t kNbinsITSsig[kDimensionsITSsig] = {AliPID::kSPECIES + 1, kMomentumBins - 1, 300, 3000};
	358	Double_t binMinITSsig[kDimensionsITSsig] = {-1., 0.1, 0., 0.};
	359	Double_t binMaxITSsig[kDimensionsITSsig] = {AliPID::kSPECIES, 10., 300., 300.};
	360	fQAlist->AddAt((histo = new THnSparseF("fCombTPCITSpid", "Combined TPC-ITS PID", kDimensionsITSsig, kNbinsITSsig, binMinITSsig, binMaxITSsig)), kITSSignal);
	361	histo->GetAxis(1)->Set(kMomentumBins - 1, momentumBins);
	362	histo->GetAxis(0)->SetTitle("Particle Species");
	363	histo->GetAxis(1)->SetTitle("p / GeV/c");
	364	histo->GetAxis(2)->SetTitle("ITS Signal / a.u.");
	365	histo->GetAxis(3)->SetTitle("TPC Signal / a.u.");
	366	}
809a4336	367	}
809a4336	368
0792aa82	369	//____________________________________________________________
	370	void AliHFEpid::InitStrategy1(){
	371	//
	372	// TPC alone, 3-sigma cut
	373	//
	374	AliHFEpidTPC *pid = new AliHFEpidTPC("strat1TPCpid");
	375	pid->SetTPCnSigma(3);
	376	Bool_t status = pid->InitializePID();
	377	if(IsQAOn() && status) pid->SetQAOn(fQAlist);
	378	if(HasMCData() && status) pid->SetHasMCData();
	379	fDetectorPID[kTPCpid] = pid;
	380	}
	381
	382	//____________________________________________________________
	383	void AliHFEpid::InitStrategy2(){
	384	//
	385	// 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
	386	//
	387	AliHFEpidTPC *pid = new AliHFEpidTPC("strat2TPCpid");
	388	pid->SetTPCnSigma(3);
	389	pid->SetAsymmetricTPCsigmaCut(2., 10., 0., 4.);
	390	Bool_t status = pid->InitializePID();
	391	if(IsQAOn() && status) pid->SetQAOn(fQAlist);
	392	if(HasMCData() && status) pid->SetHasMCData();
	393	fDetectorPID[kTPCpid] = pid;
	394	}
	395
	396	//____________________________________________________________
	397	void AliHFEpid::InitStrategy3(){
	398	//
	399	// TPC alone, symmetric 3 sigma cut and 2 - -100 sigma pion rejection
	400	//
	401	AliHFEpidTPC *pid = new AliHFEpidTPC("strat3TPCpid");
	402	pid->SetTPCnSigma(3);
	403	pid->SetRejectParticle(AliPID::kPion, 0., -100., 10., 3.);
	404	Bool_t status = pid->InitializePID();
	405	if(IsQAOn() && status) pid->SetQAOn(fQAlist);
	406	if(HasMCData() && status) pid->SetHasMCData();
	407	fDetectorPID[kTPCpid] = pid;
	408	}
	409
	410	//____________________________________________________________
	411	void AliHFEpid::InitStrategy4(){
	412	//
	413	// TPC and TRD combined, TPC 3 sigma cut and TRD NN 90% el efficiency level above 2 GeV/c
	414	//
	415	InitStrategy1();
	416	AliHFEpidTRD *trdpid = new AliHFEpidTRD("strat4TRDpid");
	417	Bool_t status = trdpid->InitializePID();
	418	if(IsQAOn() && status) trdpid->SetQAOn(fQAlist);
	419	if(HasMCData() && status) trdpid->SetHasMCData();
	420	fDetectorPID[kTRDpid] = trdpid;
	421	}
	422
	423	//____________________________________________________________
	424	void AliHFEpid::InitStrategy5(){
	425	//
	426	// TPC and TRD combined, TPC 3 sigma cut and TRD NN 90% el efficiency level above 2 GeV/c
	427	//
	428	InitStrategy1();
	429	AliHFEpidTRD *trdpid = new AliHFEpidTRD("strat5TRDpid");
	430	Bool_t status = trdpid->InitializePID();
	431	if(IsQAOn() && status) trdpid->SetQAOn(fQAlist);
	432	if(HasMCData() && status) trdpid->SetHasMCData();
433	fDetectorPID[kTRDpid] = trdpid;
434	}
435
436	//____________________________________________________________
437	Bool_t AliHFEpid::IdentifyStrategy1(AliHFEpidObject *track){
438	return TMath::Abs(fDetectorPID[kTPCpid]->IsSelected(track)) == 11;
439	}
440
441	//____________________________________________________________
442	Bool_t AliHFEpid::IdentifyStrategy2(AliHFEpidObject *track){
443	return TMath::Abs(fDetectorPID[kTPCpid]->IsSelected(track)) == 11;
444	}
445
446	//____________________________________________________________
447	Bool_t AliHFEpid::IdentifyStrategy3(AliHFEpidObject *track){
448	return TMath::Abs(fDetectorPID[kTPCpid]->IsSelected(track)) == 11;
449	}
450
451	//____________________________________________________________
452	Bool_t AliHFEpid::IdentifyStrategy4(AliHFEpidObject *track){
453	Int_t trdpid = TMath::Abs(fDetectorPID[kTRDpid]->IsSelected(track));
454	return (trdpid == 0 \|\| trdpid == 11) && (TMath::Abs(fDetectorPID[kTPCpid]->IsSelected(track)) == 11);
455	}
456
457	//____________________________________________________________
458	Bool_t AliHFEpid::IdentifyStrategy5(AliHFEpidObject *track){
459	return MakePidTpcTrd(track);
460	}