[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEpidTPC.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.                  *
**************************************************************************/
//
// Class for TPC PID
// Implements the abstract base class AliHFEpidBase
// 
// Class contains TPC specific cuts and QA histograms
// Two selection strategies are offered: Selection of certain value
// regions in the TPC dE/dx (by IsSelected), and likelihoods
//
// Authors: 
//
//   Markus Fasel <M.Fasel@gsi.de> 
//   Markus Heide <mheide@uni-muenster.de> 
//  
#include <TH2I.h>
#include <TList.h>
#include <TMath.h>
//#include <TParticle.h>

#include "AliAODTrack.h"
#include "AliAODMCParticle.h"
#include "AliESDtrack.h"
#include "AliExternalTrackParam.h"
#include "AliLog.h"
#include "AliMCParticle.h"
#include "AliPID.h"
#include "AliESDpid.h"
//#include "AliVParticle.h"

#include "AliHFEpidTPC.h"


//___________________________________________________________________
AliHFEpidTPC::AliHFEpidTPC(const char* name) :
  // add a list here
  AliHFEpidBase(name)
  , fLineCrossingType(0)
  , fLineCrossingsEnabled(0)
  , fNsigmaTPC(3)
  , fRejectionEnabled(0)
  , fPID(NULL)
  , fESDpid(NULL)
  , fQAList(NULL)
{
  //
  // default  constructor
  // 
  memset(fLineCrossingSigma, 0, sizeof(Double_t) * AliPID::kSPECIES);
  memset(fPAsigCut, 0, sizeof(Float_t) * 2);
  memset(fNAsigmaTPC, 0, sizeof(Float_t) * 2);
  fPID = new AliPID;
  fESDpid = new AliESDpid;
}

//___________________________________________________________________
AliHFEpidTPC::AliHFEpidTPC(const AliHFEpidTPC &ref) :
  AliHFEpidBase("")
  , fLineCrossingType(0)
  , fLineCrossingsEnabled(0)
  , fNsigmaTPC(2)
  , fRejectionEnabled(0)
  , fPID(NULL)
  , fESDpid(NULL)
  , fQAList(NULL)
{
  //
  // Copy constructor
  //
  ref.Copy(*this);
}

//___________________________________________________________________
AliHFEpidTPC &AliHFEpidTPC::operator=(const AliHFEpidTPC &ref){
  //
  // Assignment operator
  //
  if(this != &ref){
    ref.Copy(*this);
  } 
  return *this;
}

//___________________________________________________________________
void AliHFEpidTPC::Copy(TObject &o) const{
  //
  // Copy function 
  // called in copy constructor and assigment operator
  //
  AliHFEpidTPC &target = dynamic_cast<AliHFEpidTPC &>(o);

  target.fLineCrossingsEnabled = fLineCrossingsEnabled;
  target.fNsigmaTPC = fNsigmaTPC;
  target.fRejectionEnabled = fRejectionEnabled;
  target.fPID = new AliPID(*fPID);
  target.fESDpid = new AliESDpid(*fESDpid);
  target.fQAList = dynamic_cast<TList *>(fQAList->Clone());
  memcpy(target.fLineCrossingSigma, fLineCrossingSigma, sizeof(Double_t) * AliPID::kSPECIES);
  memcpy(target.fPAsigCut, fPAsigCut, sizeof(Float_t) * 2);
  memcpy(target.fNAsigmaTPC, fNAsigmaTPC, sizeof(Float_t) * 2);
 
  AliHFEpidBase::Copy(target);
}

//___________________________________________________________________
AliHFEpidTPC::~AliHFEpidTPC(){
  //
  // Destructor
  //
  if(fPID) delete fPID;
  if(fESDpid) delete fESDpid;
  if(fQAList){
    fQAList->Delete();
    delete fQAList;
  }
}

//___________________________________________________________________
Bool_t AliHFEpidTPC::InitializePID(){
  //
  // Add TPC dE/dx Line crossings
  //
  //AddTPCdEdxLineCrossing(AliPID::kKaon, 0.3, 0.018);
  //AddTPCdEdxLineCrossing(AliPID::kProton, 0.9, 0.054);
  return kTRUE;
}

//___________________________________________________________________
Int_t AliHFEpidTPC::IsSelected(AliHFEpidObject *track)
{
  //
  // For the TPC pid we use the 2-sigma band around the bethe bloch curve
  // for electrons
  // exclusion of the crossing points
  //
  if(track->fAnalysisType == AliHFEpidObject::kESDanalysis){
    AliESDtrack *esdTrack = NULL;
    if(!(esdTrack = dynamic_cast<AliESDtrack *>(track->fRecTrack))) return 0;
    AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(track->fMCtrack);
    return MakePIDesd(esdTrack, mctrack);
  }else{
    AliAODTrack *aodtrack = NULL;
    if(!(aodtrack = dynamic_cast<AliAODTrack *>(track->fRecTrack))) return 0;
    AliAODMCParticle *aodmctrack = dynamic_cast<AliAODMCParticle *>(track->fMCtrack);
    return MakePIDaod(aodtrack, aodmctrack);
  }
}
//___________________________________________________________________
Int_t AliHFEpidTPC::MakePIDesd(AliESDtrack *esdTrack, AliMCParticle *mctrack){
  //
  //  Doing TPC PID as explained in IsSelected for ESD tracks
  //
  if(IsQAon()) FillTPChistograms(esdTrack, mctrack);
  Float_t nsigma = fESDpid->NumberOfSigmasTPC(esdTrack, AliPID::kElectron);
  // exclude crossing points:
  // Determine the bethe values for each particle species
  Bool_t isLineCrossing = kFALSE;
  fLineCrossingType = 0;  // default value
  for(Int_t ispecies = 0; ispecies < AliPID::kSPECIES; ispecies++){
    if(ispecies == AliPID::kElectron) continue;
    if(!(fLineCrossingsEnabled & 1 << ispecies)) continue;
    if(TMath::Abs(fESDpid->NumberOfSigmasTPC(esdTrack, (AliPID::EParticleType)ispecies)) < fLineCrossingSigma[ispecies] && TMath::Abs(nsigma) < fNsigmaTPC){
      // Point in a line crossing region, no PID possible, but !PID still possible ;-)
      isLineCrossing = kTRUE;
      fLineCrossingType = ispecies;
      break;
    }
  }
  if(isLineCrossing) return 0;

  // Check particle rejection
  if(HasParticleRejection()){
    Int_t reject = Reject(esdTrack);
    if(reject != 0) return reject;
  }
  // Check whether distance from the electron line is smaller than n-sigma

  // Perform Asymmetric n-sigma cut if required, else perform symmetric TPC sigma cut
  Float_t p = 0.;
  Int_t pdg = 0;
  if(HasAsymmetricSigmaCut() && (p = esdTrack->P()) >= fPAsigCut[0] && p <= fPAsigCut[1]){ 
    if(nsigma >= fNAsigmaTPC[0] && nsigma <= fNAsigmaTPC[1]) pdg = 11; 
  } else {
    if(TMath::Abs(nsigma) < fNsigmaTPC ) pdg = 11;
  }
  if(IsQAon() && pdg != 0) (dynamic_cast<TH2I *>(fQAList->At(kHistTPCselected)))->Fill(esdTrack->GetInnerParam() ? esdTrack->GetInnerParam()->P() : esdTrack->P(), esdTrack->GetTPCsignal());

  return pdg;
}

//___________________________________________________________________
Int_t AliHFEpidTPC::MakePIDaod(AliAODTrack * /*aodTrack*/, AliAODMCParticle * /*mctrack*/){
  AliError("AOD PID not yet implemented");
  return 0;
}

//___________________________________________________________________
Int_t AliHFEpidTPC::Reject(AliESDtrack *track){
  //
  // reject particles based on asymmetric sigma cut
  //
  Int_t pdc[AliPID::kSPECIES] = {11,13,211,321,2212};
  Double_t p = track->GetOuterParam() ? track->GetOuterParam()->P() : track->P();
  for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
    if(!TESTBIT(fRejectionEnabled, ispec)) continue;
    // Particle rejection enabled
    if(p < fRejection[4*ispec] || p > fRejection[4*ispec+2]) continue;
    Double_t sigma = fESDpid->NumberOfSigmasTPC(track, static_cast<AliPID::EParticleType>(ispec));
    if(sigma >= fRejection[4*ispec+1] && sigma <= fRejection[4*ispec+3]) return pdc[ispec] * track->Charge();
  }
  return 0;
}

//___________________________________________________________________
void AliHFEpidTPC::AddTPCdEdxLineCrossing(Int_t species, Double_t sigma){
  //
  // Add exclusion point for the TPC PID where a dEdx line crosses the electron line
  // Stores line center and line sigma
  //
  if(species >= AliPID::kSPECIES){
    AliError("Species doesn't exist");
    return;
  }
  fLineCrossingsEnabled |= 1 << species;
  fLineCrossingSigma[species] = sigma;
}

//___________________________________________________________________
Double_t AliHFEpidTPC::Likelihood(const AliESDtrack *track, Int_t species, Float_t rsig)
{
  //gives probability for track to come from a certain particle species;
  //no priors considered -> probabilities for equal abundances of all species!
  //optional restriction to r-sigma bands of different particle species; default: rsig = 2. (see below)

  //IMPORTANT: Tracks which are judged to be outliers get negative likelihoods -> unusable for combination with further detectors!
  
  if(!track) return -1.;
  Bool_t outlier = kTRUE;
  // Check whether distance from the respective particle line is smaller than r sigma
  for(Int_t hypo = 0; hypo < AliPID::kSPECIES; hypo++){
    if(TMath::Abs(fESDpid->NumberOfSigmasTPC(track, (AliPID::EParticleType)hypo)) > rsig)
      outlier = kTRUE;
    else {
	    outlier = kFALSE;
	    break;
	  }
  }
  if(outlier)
    return -2.;

  Double_t tpcProb[5];

  track->GetTPCpid(tpcProb);

  return tpcProb[species];
}

//___________________________________________________________________
Double_t  AliHFEpidTPC::Suppression(const AliESDtrack *track, Int_t species)
{
  //ratio of likelihoods to be whatever species/to be an electron;
  //as a cross-check for possible particle type suppression compared to electrons
  const Double_t kVerySmall = 1e-10;
  if(!track) return -20;
  if((TMath::Abs(Likelihood(track,species) + 2) < kVerySmall)||(TMath::Abs(Likelihood(track,0) + 2 ) < kVerySmall))
    return -30;
  if(TMath::Abs(Likelihood(track,species)) < kVerySmall)
    return -10;
  if (TMath::Abs(Likelihood(track,0)) < kVerySmall)
    return 10.;
  else
    return TMath::Log10(Likelihood(track,species)/(Likelihood(track,0)));


}

//___________________________________________________________________
void AliHFEpidTPC::FillTPChistograms(const AliESDtrack *track, const AliMCParticle *mctrack){
  // 
  // Fill the QA histogtrams
  //
  if(!track)
    return;
 
  Double_t tpcSignal = track->GetTPCsignal();
  Double_t p = track->GetInnerParam() ? track->GetInnerParam()->P() : track->P();
  if(HasMCData()){
    switch(TMath::Abs(mctrack->Particle()->GetPdgCode())){
      case 11:    (dynamic_cast<TH2I *>(fQAList->At(kHistTPCelectron)))->Fill(p, tpcSignal);
	                (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobEl)))->Fill(p, Likelihood(track, 0));
	                //histograms with ratio of likelihood to be electron/to be other species (a check for quality of likelihood PID);
	                (dynamic_cast<TH2F *>(fQAList->At(kHistTPCenhanceElPi)))->Fill(p, -Suppression(track, 2));
	                (dynamic_cast<TH2F *>(fQAList->At(kHistTPCenhanceElMu)))->Fill(p, -Suppression(track, 1));
	                (dynamic_cast<TH2F *>(fQAList->At(kHistTPCenhanceElKa)))->Fill(p, -Suppression(track, 3));
	                (dynamic_cast<TH2F *>(fQAList->At(kHistTPCenhanceElPro)))->Fill(p, -Suppression(track, 4));
	                //___________________________________________________________________________________________
	                //Likelihoods for electrons to be other particle species
	                (dynamic_cast<TH2F *>(fQAList->At(kHistTPCElprobPi)))->Fill(p, Likelihood(track, 2));
	                (dynamic_cast<TH2F *>(fQAList->At(kHistTPCElprobMu)))->Fill(p, Likelihood(track, 1));
	                (dynamic_cast<TH2F *>(fQAList->At(kHistTPCElprobKa)))->Fill(p, Likelihood(track, 3));
	                (dynamic_cast<TH2F *>(fQAList->At(kHistTPCElprobPro)))->Fill(p, Likelihood(track, 4));
	                break;
	    //___________________________________________________________________________________________
      case 13:    (dynamic_cast<TH2I *>(fQAList->At(kHistTPCmuon)))->Fill(p, tpcSignal);
                  //Likelihood of muon to be an electron
                  (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobMu)))->Fill(p, Likelihood(track, 0));
                  //ratio of likelihood for muon to be a muon/an electron -> indicator for quality of muon suppression
                  //below functions are the same for other species
	                (dynamic_cast<TH2F *>(fQAList->At(kHistTPCsuppressMu)))->Fill(p, Suppression(track, 1));
                  break;
      case 211:   (dynamic_cast<TH2I *>(fQAList->At(kHistTPCpion)))->Fill(p, tpcSignal);
	                (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobPi)))->Fill(p, Likelihood(track, 0));
	                (dynamic_cast<TH2F *>(fQAList->At(kHistTPCsuppressPi)))->Fill(p, Suppression(track, 2));
                  break;
      case 321:   (dynamic_cast<TH2I *>(fQAList->At(kHistTPCkaon)))->Fill(p, tpcSignal);
	                (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobKa)))->Fill(p, Likelihood(track, 0));
	                (dynamic_cast<TH2F *>(fQAList->At(kHistTPCsuppressKa)))->Fill(p, Suppression(track, 3));
                  break;
      case 2212:  (dynamic_cast<TH2I *>(fQAList->At(kHistTPCproton)))->Fill(p, tpcSignal);
	                (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobPro)))->Fill(p, Likelihood(track, 0));
	                (dynamic_cast<TH2F *>(fQAList->At(kHistTPCsuppressPro)))->Fill(p, Suppression(track, 4));
                  break;
      default:    (dynamic_cast<TH2I *>(fQAList->At(kHistTPCothers)))->Fill(p, tpcSignal);
	                (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobOth)))->Fill(p, Likelihood(track, 0));

	                break;
    }
  }
  //TPC signal and Likelihood to be electron for all tracks (independent of MC information)
  (dynamic_cast<TH2I *>(fQAList->At(kHistTPCall)))->Fill(p, tpcSignal);
  (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobAll)))->Fill(p, Likelihood(track, 0));
}

//___________________________________________________________________
void AliHFEpidTPC::AddQAhistograms(TList *qaList){
  //
  // Create QA histograms for TPC PID
  //
  fQAList = new TList;
  fQAList->SetName("fTPCqaHistos");

  fQAList->AddAt(new TH2I("fHistTPCelectron","TPC signal for Electrons", 200, 0, 20, 60, 0, 600), kHistTPCelectron); 
  fQAList->AddAt(new TH2I("fHistTPCmuon","TPC signal for Muons", 200, 0, 20, 60, 0, 600), kHistTPCmuon);
  fQAList->AddAt(new TH2I("fHistTPCpion","TPC signal for Pions", 200, 0, 20, 60, 0, 600), kHistTPCpion);
  fQAList->AddAt(new TH2I("fHistTPCkaon","TPC signal for Kaons", 200, 0, 20, 60, 0, 600), kHistTPCkaon);
  fQAList->AddAt(new TH2I("fHistTPCproton","TPC signal for Protons", 200, 0, 20, 60, 0, 600), kHistTPCproton);
  fQAList->AddAt(new TH2I("fHistTPCothers","TPC signal for other species", 200, 0, 20, 60, 0, 600), kHistTPCothers);
  fQAList->AddAt(new TH2I("fHistTPCall","TPC signal for all species", 200, 0, 20, 60, 0, 600), kHistTPCall);
  fQAList->AddAt(new TH2I("fHistTPCselected","TPC signal for all selected particles", 200, 0, 20, 60, 0, 600), kHistTPCselected);

  fQAList->AddAt(new TH2F("fHistTPCprobEl","TPC likelihood for electrons to be an electron vs. p", 200, 0.,20.,200,0.,1.), kHistTPCprobEl);
  fQAList->AddAt(new TH2F("fHistTPCprobPi","TPC likelihood for pions to be an electron vs. p",  200, 0.,20.,200, 0.,1.), kHistTPCprobPi);
  fQAList->AddAt(new TH2F("fHistTPCprobMu","TPC likelihood for muons to be an electron vs. p",  200, 0.,20.,200, 0.,1.), kHistTPCprobMu);
  fQAList->AddAt(new TH2F("fHistTPCprobKa","TPC likelihood for kaons to be an electron vs. p",  200, 0.,20.,200, 0.,1.), kHistTPCprobKa);
  fQAList->AddAt(new TH2F("fHistTPCprobPro","TPC likelihood for protons to be an electron vs. p",  200, 0.,20.,200, 0.,1.), kHistTPCprobPro);
  fQAList->AddAt(new TH2F("fHistTPCprobOth","TPC likelihood for other particles to be an electron vs. p",  200, 0.,20.,200, 0.,1.), kHistTPCprobOth);
  fQAList->AddAt(new TH2F("fHistTPCprobAll","TPC likelihood for all particles to be an electron vs. p",  200, 0.,20.,200, 0.,1.), kHistTPCprobAll);


 fQAList->AddAt(new TH2F("fHistTPCsuppressPi","log10 of TPC Likelihood(pion)/Likelihood(elec) for pions vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressPi);
 fQAList->AddAt(new TH2F("fHistTPCsuppressMu","log10 of TPC Likelihood(muon)/Likelihood(elec) for muons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressMu);
 fQAList->AddAt(new TH2F("fHistTPCsuppressKa","log10 of TPC Likelihood(kaon)/Likelihood(elec) for kaons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressKa);
 fQAList->AddAt(new TH2F("fHistTPCsuppressPro","log10 of TPC Likelihood(proton)/Likelihood(elec)for protons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressPro);

 fQAList->AddAt(new TH2F("fHistTPCenhanceElPi","log10 of TPC Likelihood(elec)/Likelihood(pion) for electrons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressPi);
 fQAList->AddAt(new TH2F("fHistTPCenhanceElMu","log10 of TPC Likelihood(elec)/Likelihood(muon) for electrons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressMu);
 fQAList->AddAt(new TH2F("fHistTPCenhanceElKa","log10 of TPC Likelihood(elec)/Likelihood(kaon) for electrons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressKa);
 fQAList->AddAt(new TH2F("fHistTPCenhanceElPro","log10 of TPC Likelihood(elec)/Likelihood(proton) for electrons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressPro);

 fQAList->AddAt(new TH2F("fHistTPCElprobPi","TPC likelihood for electrons to be a pion vs. p", 200, 0.,20.,200,0.,1.), kHistTPCElprobPi);
 fQAList->AddAt(new TH2F("fHistTPCElprobMu","TPC likelihood for electrons to be a muon vs. p", 200, 0.,20.,200,0.,1.), kHistTPCElprobMu);
 fQAList->AddAt(new TH2F("fHistTPCElprobKa","TPC likelihood for electrons to be a kaon vs. p", 200, 0.,20.,200,0.,1.), kHistTPCElprobKa);
 fQAList->AddAt(new TH2F("fHistTPCElprobPro","TPC likelihood for electrons to be a proton vs. p", 200, 0.,20.,200,0.,1.), kHistTPCElprobPro);


  qaList->AddLast(fQAList);
}
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	**************************************************************************/
50685501	15	//
	16	// Class for TPC PID
	17	// Implements the abstract base class AliHFEpidBase
	18	//
	19	// Class contains TPC specific cuts and QA histograms
	20	// Two selection strategies are offered: Selection of certain value
	21	// regions in the TPC dE/dx (by IsSelected), and likelihoods
	22	//
	23	// Authors:
	24	//
	25	// Markus Fasel <M.Fasel@gsi.de>
	26	// Markus Heide <mheide@uni-muenster.de>
	27	//
809a4336	28	#include <TH2I.h>
	29	#include <TList.h>
	30	#include <TMath.h>
50685501	31	//#include <TParticle.h>
809a4336	32
722347d8	33	#include "AliAODTrack.h"
722347d8	34	#include "AliAODMCParticle.h"
809a4336	35	#include "AliESDtrack.h"
	36	#include "AliExternalTrackParam.h"
	37	#include "AliLog.h"
722347d8	38	#include "AliMCParticle.h"
809a4336	39	#include "AliPID.h"
10d100d4	40	#include "AliESDpid.h"
50685501	41	//#include "AliVParticle.h"
809a4336	42
	43	#include "AliHFEpidTPC.h"
	44
	45
	46
	47	//___________________________________________________________________
	48	AliHFEpidTPC::AliHFEpidTPC(const char* name) :
	49	// add a list here
9bcfd1ab	50	AliHFEpidBase(name)
9bcfd1ab	51	, fLineCrossingType(0)
809a4336	52	, fLineCrossingsEnabled(0)
75d81601	53	, fNsigmaTPC(3)
0792aa82	54	, fRejectionEnabled(0)
75d81601	55	, fPID(NULL)
10d100d4	56	, fESDpid(NULL)
75d81601	57	, fQAList(NULL)
809a4336	58	{
	59	//
	60	// default constructor
	61	//
809a4336	62	memset(fLineCrossingSigma, 0, sizeof(Double_t) * AliPID::kSPECIES);
722347d8	63	memset(fPAsigCut, 0, sizeof(Float_t) * 2);
722347d8	64	memset(fNAsigmaTPC, 0, sizeof(Float_t) * 2);
809a4336	65	fPID = new AliPID;
10d100d4	66	fESDpid = new AliESDpid;
809a4336	67	}
	68
	69	//___________________________________________________________________
	70	AliHFEpidTPC::AliHFEpidTPC(const AliHFEpidTPC &ref) :
9bcfd1ab	71	AliHFEpidBase("")
9bcfd1ab	72	, fLineCrossingType(0)
809a4336	73	, fLineCrossingsEnabled(0)
809a4336	74	, fNsigmaTPC(2)
0792aa82	75	, fRejectionEnabled(0)
75d81601	76	, fPID(NULL)
10d100d4	77	, fESDpid(NULL)
75d81601	78	, fQAList(NULL)
809a4336	79	{
	80	//
	81	// Copy constructor
	82	//
	83	ref.Copy(*this);
	84	}
	85
	86	//___________________________________________________________________
	87	AliHFEpidTPC &AliHFEpidTPC::operator=(const AliHFEpidTPC &ref){
	88	//
	89	// Assignment operator
	90	//
	91	if(this != &ref){
	92	ref.Copy(*this);
	93	}
	94	return *this;
	95	}
	96
75d81601	97	//___________________________________________________________________
809a4336	98	void AliHFEpidTPC::Copy(TObject &o) const{
	99	//
	100	// Copy function
	101	// called in copy constructor and assigment operator
	102	//
	103	AliHFEpidTPC &target = dynamic_cast<AliHFEpidTPC &>(o);
	104
	105	target.fLineCrossingsEnabled = fLineCrossingsEnabled;
	106	target.fNsigmaTPC = fNsigmaTPC;
0792aa82	107	target.fRejectionEnabled = fRejectionEnabled;
809a4336	108	target.fPID = new AliPID(*fPID);
10d100d4	109	target.fESDpid = new AliESDpid(*fESDpid);
809a4336	110	target.fQAList = dynamic_cast<TList *>(fQAList->Clone());
75d81601	111	memcpy(target.fLineCrossingSigma, fLineCrossingSigma, sizeof(Double_t) * AliPID::kSPECIES);
722347d8	112	memcpy(target.fPAsigCut, fPAsigCut, sizeof(Float_t) * 2);
	113	memcpy(target.fNAsigmaTPC, fNAsigmaTPC, sizeof(Float_t) * 2);
	114
809a4336	115	AliHFEpidBase::Copy(target);
	116	}
	117
	118	//___________________________________________________________________
	119	AliHFEpidTPC::~AliHFEpidTPC(){
	120	//
	121	// Destructor
	122	//
	123	if(fPID) delete fPID;
10d100d4	124	if(fESDpid) delete fESDpid;
809a4336	125	if(fQAList){
	126	fQAList->Delete();
	127	delete fQAList;
	128	}
	129	}
	130
	131	//___________________________________________________________________
	132	Bool_t AliHFEpidTPC::InitializePID(){
	133	//
	134	// Add TPC dE/dx Line crossings
	135	//
75d81601	136	//AddTPCdEdxLineCrossing(AliPID::kKaon, 0.3, 0.018);
75d81601	137	//AddTPCdEdxLineCrossing(AliPID::kProton, 0.9, 0.054);
809a4336	138	return kTRUE;
	139	}
	140
	141	//___________________________________________________________________
722347d8	142	Int_t AliHFEpidTPC::IsSelected(AliHFEpidObject *track)
809a4336	143	{
	144	//
	145	// For the TPC pid we use the 2-sigma band around the bethe bloch curve
	146	// for electrons
	147	// exclusion of the crossing points
	148	//
722347d8	149	if(track->fAnalysisType == AliHFEpidObject::kESDanalysis){
	150	AliESDtrack *esdTrack = NULL;
	151	if(!(esdTrack = dynamic_cast<AliESDtrack *>(track->fRecTrack))) return 0;
	152	AliMCParticle mctrack = dynamic_cast<AliMCParticle >(track->fMCtrack);
	153	return MakePIDesd(esdTrack, mctrack);
	154	}else{
	155	AliAODTrack *aodtrack = NULL;
	156	if(!(aodtrack = dynamic_cast<AliAODTrack *>(track->fRecTrack))) return 0;
	157	AliAODMCParticle aodmctrack = dynamic_cast<AliAODMCParticle >(track->fMCtrack);
	158	return MakePIDaod(aodtrack, aodmctrack);
	159	}
	160	}
722347d8	161	//___________________________________________________________________
	162	Int_t AliHFEpidTPC::MakePIDesd(AliESDtrack esdTrack, AliMCParticle mctrack){
	163	//
	164	// Doing TPC PID as explained in IsSelected for ESD tracks
	165	//
	166	if(IsQAon()) FillTPChistograms(esdTrack, mctrack);
10d100d4	167	Float_t nsigma = fESDpid->NumberOfSigmasTPC(esdTrack, AliPID::kElectron);
809a4336	168	// exclude crossing points:
809a4336	169	// Determine the bethe values for each particle species
809a4336	170	Bool_t isLineCrossing = kFALSE;
9bcfd1ab	171	fLineCrossingType = 0; // default value
809a4336	172	for(Int_t ispecies = 0; ispecies < AliPID::kSPECIES; ispecies++){
75d81601	173	if(ispecies == AliPID::kElectron) continue;
809a4336	174	if(!(fLineCrossingsEnabled & 1 << ispecies)) continue;
10d100d4	175	if(TMath::Abs(fESDpid->NumberOfSigmasTPC(esdTrack, (AliPID::EParticleType)ispecies)) < fLineCrossingSigma[ispecies] && TMath::Abs(nsigma) < fNsigmaTPC){
9bcfd1ab	176	// Point in a line crossing region, no PID possible, but !PID still possible ;-)
10d100d4	177	isLineCrossing = kTRUE;
9bcfd1ab	178	fLineCrossingType = ispecies;
809a4336	179	break;
	180	}
	181	}
	182	if(isLineCrossing) return 0;
0792aa82	183
	184	// Check particle rejection
	185	if(HasParticleRejection()){
	186	Int_t reject = Reject(esdTrack);
0792aa82	187	if(reject != 0) return reject;
0792aa82	188	}
809a4336	189	// Check whether distance from the electron line is smaller than n-sigma
722347d8	190
722347d8	191	// Perform Asymmetric n-sigma cut if required, else perform symmetric TPC sigma cut
722347d8	192	Float_t p = 0.;
0792aa82	193	Int_t pdg = 0;
722347d8	194	if(HasAsymmetricSigmaCut() && (p = esdTrack->P()) >= fPAsigCut[0] && p <= fPAsigCut[1]){
0792aa82	195	if(nsigma >= fNAsigmaTPC[0] && nsigma <= fNAsigmaTPC[1]) pdg = 11;
722347d8	196	} else {
0792aa82	197	if(TMath::Abs(nsigma) < fNsigmaTPC ) pdg = 11;
722347d8	198	}
0792aa82	199	if(IsQAon() && pdg != 0) (dynamic_cast<TH2I *>(fQAList->At(kHistTPCselected)))->Fill(esdTrack->GetInnerParam() ? esdTrack->GetInnerParam()->P() : esdTrack->P(), esdTrack->GetTPCsignal());
9bcfd1ab	200
0792aa82	201	return pdg;
722347d8	202	}
	203
	204	//___________________________________________________________________
	205	Int_t AliHFEpidTPC::MakePIDaod(AliAODTrack * /aodTrack/, AliAODMCParticle * /mctrack/){
	206	AliError("AOD PID not yet implemented");
809a4336	207	return 0;
	208	}
	209
0792aa82	210	//___________________________________________________________________
	211	Int_t AliHFEpidTPC::Reject(AliESDtrack *track){
	212	//
	213	// reject particles based on asymmetric sigma cut
	214	//
	215	Int_t pdc[AliPID::kSPECIES] = {11,13,211,321,2212};
	216	Double_t p = track->GetOuterParam() ? track->GetOuterParam()->P() : track->P();
	217	for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
	218	if(!TESTBIT(fRejectionEnabled, ispec)) continue;
0792aa82	219	// Particle rejection enabled
0792aa82	220	if(p < fRejection[4ispec] \|\| p > fRejection[4ispec+2]) continue;
10d100d4	221	Double_t sigma = fESDpid->NumberOfSigmasTPC(track, static_cast<AliPID::EParticleType>(ispec));
0792aa82	222	if(sigma >= fRejection[4ispec+1] && sigma <= fRejection[4ispec+3]) return pdc[ispec] * track->Charge();
	223	}
	224	return 0;
	225	}
	226
809a4336	227	//___________________________________________________________________
75d81601	228	void AliHFEpidTPC::AddTPCdEdxLineCrossing(Int_t species, Double_t sigma){
809a4336	229	//
	230	// Add exclusion point for the TPC PID where a dEdx line crosses the electron line
	231	// Stores line center and line sigma
	232	//
	233	if(species >= AliPID::kSPECIES){
	234	AliError("Species doesn't exist");
	235	return;
	236	}
	237	fLineCrossingsEnabled \|= 1 << species;
75d81601	238	fLineCrossingSigma[species] = sigma;
809a4336	239	}
	240
	241	//___________________________________________________________________
	242	Double_t AliHFEpidTPC::Likelihood(const AliESDtrack *track, Int_t species, Float_t rsig)
	243	{
	244	//gives probability for track to come from a certain particle species;
	245	//no priors considered -> probabilities for equal abundances of all species!
	246	//optional restriction to r-sigma bands of different particle species; default: rsig = 2. (see below)
	247
	248	//IMPORTANT: Tracks which are judged to be outliers get negative likelihoods -> unusable for combination with further detectors!
	249
	250	if(!track) return -1.;
809a4336	251	Bool_t outlier = kTRUE;
809a4336	252	// Check whether distance from the respective particle line is smaller than r sigma
75d81601	253	for(Int_t hypo = 0; hypo < AliPID::kSPECIES; hypo++){
10d100d4	254	if(TMath::Abs(fESDpid->NumberOfSigmasTPC(track, (AliPID::EParticleType)hypo)) > rsig)
75d81601	255	outlier = kTRUE;
	256	else {
	257	outlier = kFALSE;
	258	break;
	259	}
	260	}
809a4336	261	if(outlier)
	262	return -2.;
	263
75d81601	264	Double_t tpcProb[5];
809a4336	265
75d81601	266	track->GetTPCpid(tpcProb);
809a4336	267
75d81601	268	return tpcProb[species];
809a4336	269	}
809a4336	270
809a4336	271	//___________________________________________________________________
	272	Double_t AliHFEpidTPC::Suppression(const AliESDtrack *track, Int_t species)
	273	{
	274	//ratio of likelihoods to be whatever species/to be an electron;
	275	//as a cross-check for possible particle type suppression compared to electrons
50685501	276	const Double_t kVerySmall = 1e-10;
809a4336	277	if(!track) return -20;
50685501	278	if((TMath::Abs(Likelihood(track,species) + 2) < kVerySmall)\|\|(TMath::Abs(Likelihood(track,0) + 2 ) < kVerySmall))
809a4336	279	return -30;
50685501	280	if(TMath::Abs(Likelihood(track,species)) < kVerySmall)
809a4336	281	return -10;
50685501	282	if (TMath::Abs(Likelihood(track,0)) < kVerySmall)
809a4336	283	return 10.;
	284	else
	285	return TMath::Log10(Likelihood(track,species)/(Likelihood(track,0)));
	286
	287
	288	}
75d81601	289
809a4336	290	//___________________________________________________________________
722347d8	291	void AliHFEpidTPC::FillTPChistograms(const AliESDtrack track, const AliMCParticle mctrack){
50685501	292	//
50685501	293	// Fill the QA histogtrams
809a4336	294	//
	295	if(!track)
	296	return;
	297
75d81601	298	Double_t tpcSignal = track->GetTPCsignal();
809a4336	299	Double_t p = track->GetInnerParam() ? track->GetInnerParam()->P() : track->P();
809a4336	300	if(HasMCData()){
722347d8	301	switch(TMath::Abs(mctrack->Particle()->GetPdgCode())){
	302	case 11: (dynamic_cast<TH2I *>(fQAList->At(kHistTPCelectron)))->Fill(p, tpcSignal);
	303	(dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobEl)))->Fill(p, Likelihood(track, 0));
	304	//histograms with ratio of likelihood to be electron/to be other species (a check for quality of likelihood PID);
	305	(dynamic_cast<TH2F *>(fQAList->At(kHistTPCenhanceElPi)))->Fill(p, -Suppression(track, 2));
	306	(dynamic_cast<TH2F *>(fQAList->At(kHistTPCenhanceElMu)))->Fill(p, -Suppression(track, 1));
	307	(dynamic_cast<TH2F *>(fQAList->At(kHistTPCenhanceElKa)))->Fill(p, -Suppression(track, 3));
	308	(dynamic_cast<TH2F *>(fQAList->At(kHistTPCenhanceElPro)))->Fill(p, -Suppression(track, 4));
	309	//___________________________________________________________________________________________
	310	//Likelihoods for electrons to be other particle species
	311	(dynamic_cast<TH2F *>(fQAList->At(kHistTPCElprobPi)))->Fill(p, Likelihood(track, 2));
	312	(dynamic_cast<TH2F *>(fQAList->At(kHistTPCElprobMu)))->Fill(p, Likelihood(track, 1));
	313	(dynamic_cast<TH2F *>(fQAList->At(kHistTPCElprobKa)))->Fill(p, Likelihood(track, 3));
	314	(dynamic_cast<TH2F *>(fQAList->At(kHistTPCElprobPro)))->Fill(p, Likelihood(track, 4));
	315	break;
	316	//___________________________________________________________________________________________
	317	case 13: (dynamic_cast<TH2I *>(fQAList->At(kHistTPCmuon)))->Fill(p, tpcSignal);
	318	//Likelihood of muon to be an electron
	319	(dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobMu)))->Fill(p, Likelihood(track, 0));
	320	//ratio of likelihood for muon to be a muon/an electron -> indicator for quality of muon suppression
	321	//below functions are the same for other species
	322	(dynamic_cast<TH2F *>(fQAList->At(kHistTPCsuppressMu)))->Fill(p, Suppression(track, 1));
	323	break;
	324	case 211: (dynamic_cast<TH2I *>(fQAList->At(kHistTPCpion)))->Fill(p, tpcSignal);
	325	(dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobPi)))->Fill(p, Likelihood(track, 0));
	326	(dynamic_cast<TH2F *>(fQAList->At(kHistTPCsuppressPi)))->Fill(p, Suppression(track, 2));
	327	break;
	328	case 321: (dynamic_cast<TH2I *>(fQAList->At(kHistTPCkaon)))->Fill(p, tpcSignal);
	329	(dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobKa)))->Fill(p, Likelihood(track, 0));
	330	(dynamic_cast<TH2F *>(fQAList->At(kHistTPCsuppressKa)))->Fill(p, Suppression(track, 3));
	331	break;
	332	case 2212: (dynamic_cast<TH2I *>(fQAList->At(kHistTPCproton)))->Fill(p, tpcSignal);
	333	(dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobPro)))->Fill(p, Likelihood(track, 0));
	334	(dynamic_cast<TH2F *>(fQAList->At(kHistTPCsuppressPro)))->Fill(p, Suppression(track, 4));
	335	break;
	336	default: (dynamic_cast<TH2I *>(fQAList->At(kHistTPCothers)))->Fill(p, tpcSignal);
	337	(dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobOth)))->Fill(p, Likelihood(track, 0));
809a4336	338
722347d8	339	break;
809a4336	340	}
	341	}
	342	//TPC signal and Likelihood to be electron for all tracks (independent of MC information)
75d81601	343	(dynamic_cast<TH2I *>(fQAList->At(kHistTPCall)))->Fill(p, tpcSignal);
722347d8	344	(dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobAll)))->Fill(p, Likelihood(track, 0));
809a4336	345	}
	346
	347	//___________________________________________________________________
	348	void AliHFEpidTPC::AddQAhistograms(TList *qaList){
50685501	349	//
	350	// Create QA histograms for TPC PID
	351	//
809a4336	352	fQAList = new TList;
	353	fQAList->SetName("fTPCqaHistos");
	354
9bcfd1ab	355	fQAList->AddAt(new TH2I("fHistTPCelectron","TPC signal for Electrons", 200, 0, 20, 60, 0, 600), kHistTPCelectron);
	356	fQAList->AddAt(new TH2I("fHistTPCmuon","TPC signal for Muons", 200, 0, 20, 60, 0, 600), kHistTPCmuon);
	357	fQAList->AddAt(new TH2I("fHistTPCpion","TPC signal for Pions", 200, 0, 20, 60, 0, 600), kHistTPCpion);
	358	fQAList->AddAt(new TH2I("fHistTPCkaon","TPC signal for Kaons", 200, 0, 20, 60, 0, 600), kHistTPCkaon);
	359	fQAList->AddAt(new TH2I("fHistTPCproton","TPC signal for Protons", 200, 0, 20, 60, 0, 600), kHistTPCproton);
	360	fQAList->AddAt(new TH2I("fHistTPCothers","TPC signal for other species", 200, 0, 20, 60, 0, 600), kHistTPCothers);
	361	fQAList->AddAt(new TH2I("fHistTPCall","TPC signal for all species", 200, 0, 20, 60, 0, 600), kHistTPCall);
	362	fQAList->AddAt(new TH2I("fHistTPCselected","TPC signal for all selected particles", 200, 0, 20, 60, 0, 600), kHistTPCselected);
809a4336	363
	364	fQAList->AddAt(new TH2F("fHistTPCprobEl","TPC likelihood for electrons to be an electron vs. p", 200, 0.,20.,200,0.,1.), kHistTPCprobEl);
	365	fQAList->AddAt(new TH2F("fHistTPCprobPi","TPC likelihood for pions to be an electron vs. p", 200, 0.,20.,200, 0.,1.), kHistTPCprobPi);
	366	fQAList->AddAt(new TH2F("fHistTPCprobMu","TPC likelihood for muons to be an electron vs. p", 200, 0.,20.,200, 0.,1.), kHistTPCprobMu);
	367	fQAList->AddAt(new TH2F("fHistTPCprobKa","TPC likelihood for kaons to be an electron vs. p", 200, 0.,20.,200, 0.,1.), kHistTPCprobKa);
	368	fQAList->AddAt(new TH2F("fHistTPCprobPro","TPC likelihood for protons to be an electron vs. p", 200, 0.,20.,200, 0.,1.), kHistTPCprobPro);
	369	fQAList->AddAt(new TH2F("fHistTPCprobOth","TPC likelihood for other particles to be an electron vs. p", 200, 0.,20.,200, 0.,1.), kHistTPCprobOth);
	370	fQAList->AddAt(new TH2F("fHistTPCprobAll","TPC likelihood for all particles to be an electron vs. p", 200, 0.,20.,200, 0.,1.), kHistTPCprobAll);
	371
	372
	373	fQAList->AddAt(new TH2F("fHistTPCsuppressPi","log10 of TPC Likelihood(pion)/Likelihood(elec) for pions vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressPi);
	374	fQAList->AddAt(new TH2F("fHistTPCsuppressMu","log10 of TPC Likelihood(muon)/Likelihood(elec) for muons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressMu);
	375	fQAList->AddAt(new TH2F("fHistTPCsuppressKa","log10 of TPC Likelihood(kaon)/Likelihood(elec) for kaons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressKa);
	376	fQAList->AddAt(new TH2F("fHistTPCsuppressPro","log10 of TPC Likelihood(proton)/Likelihood(elec)for protons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressPro);
	377
	378	fQAList->AddAt(new TH2F("fHistTPCenhanceElPi","log10 of TPC Likelihood(elec)/Likelihood(pion) for electrons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressPi);
	379	fQAList->AddAt(new TH2F("fHistTPCenhanceElMu","log10 of TPC Likelihood(elec)/Likelihood(muon) for electrons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressMu);
	380	fQAList->AddAt(new TH2F("fHistTPCenhanceElKa","log10 of TPC Likelihood(elec)/Likelihood(kaon) for electrons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressKa);
	381	fQAList->AddAt(new TH2F("fHistTPCenhanceElPro","log10 of TPC Likelihood(elec)/Likelihood(proton) for electrons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressPro);
	382
	383	fQAList->AddAt(new TH2F("fHistTPCElprobPi","TPC likelihood for electrons to be a pion vs. p", 200, 0.,20.,200,0.,1.), kHistTPCElprobPi);
	384	fQAList->AddAt(new TH2F("fHistTPCElprobMu","TPC likelihood for electrons to be a muon vs. p", 200, 0.,20.,200,0.,1.), kHistTPCElprobMu);
	385	fQAList->AddAt(new TH2F("fHistTPCElprobKa","TPC likelihood for electrons to be a kaon vs. p", 200, 0.,20.,200,0.,1.), kHistTPCElprobKa);
	386	fQAList->AddAt(new TH2F("fHistTPCElprobPro","TPC likelihood for electrons to be a proton vs. p", 200, 0.,20.,200,0.,1.), kHistTPCElprobPro);
	387
	388
	389	qaList->AddLast(fQAList);
	390	}