X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=STEER%2FAOD%2FAliAODpidUtil.cxx;h=4c61bb46bd55150f54d0fa1ae53f012a766dfd38;hb=7e68d49f13a55ff01ccf9bb45e8ba4cab4402592;hp=a005b67289dd10f1b26a9e96392bfe867c472426;hpb=7330f0e52a55c9762ac260712cb54cd9446b876d;p=u%2Fmrichter%2FAliRoot.git diff --git a/STEER/AOD/AliAODpidUtil.cxx b/STEER/AOD/AliAODpidUtil.cxx index a005b67289d..4c61bb46bd5 100644 --- a/STEER/AOD/AliAODpidUtil.cxx +++ b/STEER/AOD/AliAODpidUtil.cxx @@ -22,6 +22,7 @@ // Origin: Rosa Romita, GSI, r.romita@gsi.de //----------------------------------------------------------------- +#include "TRandom.h" #include "AliLog.h" #include "AliPID.h" #include "AliAODpidUtil.h" @@ -30,200 +31,173 @@ #include "AliAODPid.h" #include "AliTRDPIDResponse.h" #include "AliESDtrack.h" +#include "AliAODMCHeader.h" +#include "AliAODMCParticle.h" +#include "AliTOFPIDParams.h" + +#include ClassImp(AliAODpidUtil) - Int_t AliAODpidUtil::MakePID(const AliAODTrack *track,Double_t *p) const { - // - // Calculate probabilities for all detectors, except if TPConly==kTRUE - // and combine PID - // - // Option TPConly==kTRUE is used during reconstruction, - // because ITS tracking uses TPC pid - // HMPID and TRD pid are done in detector reconstructors - // +//_________________________________________________________________________ +Float_t AliAODpidUtil::GetTPCsignalTunedOnData(const AliVTrack *t) const { + AliAODTrack *track = (AliAODTrack *) t; + Float_t dedx = track->GetTPCsignalTunedOnData(); + if(dedx > 0) return dedx; + + dedx = t->GetTPCsignal(); + track->SetTPCsignalTunedOnData(dedx); + + if(dedx < 20) return dedx; + + + AliPID::EParticleType type = AliPID::kPion; + + AliAODMCHeader *mcHeader = dynamic_cast(track->GetAODEvent()->GetList()->FindObject(AliAODMCHeader::StdBranchName())); + if (mcHeader) { + TClonesArray *mcArray = (TClonesArray*)track->GetAODEvent()->GetList()->FindObject(AliAODMCParticle::StdBranchName()); + + Bool_t kGood = kTRUE; + + Int_t iS = TMath::Abs(((AliAODMCParticle*)mcArray->At(TMath::Abs(t->GetLabel())))->GetPdgCode()); + if(iS==AliPID::ParticleCode(AliPID::kElectron)){ + type = AliPID::kElectron; + } + else if(iS==AliPID::ParticleCode(AliPID::kMuon)){ + type = AliPID::kMuon; + } + else if(iS==AliPID::ParticleCode(AliPID::kPion)){ + type = AliPID::kPion; + } + else if(iS==AliPID::ParticleCode(AliPID::kKaon)){ + type = AliPID::kKaon; + } + else if(iS==AliPID::ParticleCode(AliPID::kProton)){ + type = AliPID::kProton; + } + else if(iS==AliPID::ParticleCode(AliPID::kDeuteron)){ // d + type = AliPID::kDeuteron; + } + else if(iS==AliPID::ParticleCode(AliPID::kTriton)){ // t + type = AliPID::kTriton; + } + else if(iS==AliPID::ParticleCode(AliPID::kHe3)){ // 3He + type = AliPID::kHe3; + } + else if(iS==AliPID::ParticleCode(AliPID::kAlpha)){ // 4He + type = AliPID::kAlpha; + } + else + kGood = kFALSE; + + if(kGood){ + //TODO maybe introduce different dEdxSources? + Double_t bethe = fTPCResponse.GetExpectedSignal(track, type, AliTPCPIDResponse::kdEdxDefault, this->UseTPCEtaCorrection(), + this->UseTPCMultiplicityCorrection()); + Double_t sigma = fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, this->UseTPCEtaCorrection(), + this->UseTPCMultiplicityCorrection()); + dedx = gRandom->Gaus(bethe,sigma); + +// if(iS == AliPID::ParticleCode(AliPID::kHe3) || iS == AliPID::ParticleCode(AliPID::kAlpha)) dedx *= 5; + } - /* - Float_t TimeZeroTOF = 0; - if (subtractT0) - TimeZeroTOF = event->GetT0(); - */ - Int_t ns=AliPID::kSPECIES; - Double_t tpcPid[AliPID::kSPECIES]; - MakeTPCPID(track,tpcPid); - Double_t itsPid[AliPID::kSPECIES]; - Double_t tofPid[AliPID::kSPECIES]; - Double_t trdPid[AliPID::kSPECIES]; - MakeITSPID(track,itsPid); - MakeTOFPID(track,tofPid); - //MakeHMPIDPID(track); - MakeTRDPID(track,trdPid); - for (Int_t j=0; jSetTPCsignalTunedOnData(dedx); + return dedx; } //_________________________________________________________________________ -void AliAODpidUtil::MakeTPCPID(const AliAODTrack *track,Double_t *p) const -{ - // - // TPC pid using bethe-bloch and gaussian response - // - - if ((track->GetStatus()&AliESDtrack::kTPCin )==0) return; - - AliAODPid *pidObj = track->GetDetPid(); - Double_t mom = track->P(); - Double_t dedx = 0.; - UShort_t nTPCClus = 0; - if (pidObj) { - nTPCClus = pidObj->GetTPCsignalN(); - dedx = pidObj->GetTPCsignal(); - mom = pidObj->GetTPCmomentum(); - } - - Bool_t mismatch=kTRUE; - - for (Int_t j=0; j fRange*sigma) { - p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma; - } else { - p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma; - mismatch=kFALSE; +Float_t AliAODpidUtil::GetTOFsignalTunedOnData(const AliVTrack *t) const { + AliAODTrack *track = (AliAODTrack *) t; + Double_t tofSignal = track->GetTOFsignalTunedOnData(); + + if(tofSignal < 99999) return (Float_t)tofSignal; // it has been already set + + // read additional mismatch fraction + Float_t addmism = GetTOFPIDParams()->GetTOFadditionalMismForMC(); + if(addmism > 1.){ + Float_t centr = GetCurrentCentrality(); + if(centr > 50) addmism *= 0.1667; + else if(centr > 20) addmism *= 0.33; } - } - - if (mismatch) - for (Int_t j=0; jGetDetPid(); + tofSignal = pidObj->GetTOFsignal() + fTOFResponse.GetTailRandomValue(t->Pt(),t->Eta(),pidObj->GetTOFsignal(),addmism); + track->SetTOFsignalTunedOnData(tofSignal); + return (Float_t)tofSignal; } + //_________________________________________________________________________ -void AliAODpidUtil::MakeITSPID(const AliAODTrack *track,Double_t *p) const +Float_t AliAODpidUtil::GetSignalDeltaTOFold(const AliVParticle *vtrack, AliPID::EParticleType type, Bool_t ratio/*=kFALSE*/) const { // - // ITS PID - // 1) Truncated mean method + // Number of sigma implementation for the TOF // - - - if ((track->GetStatus()&AliESDtrack::kITSin)==0) return; - UChar_t clumap=track->GetITSClusterMap(); - Int_t nPointsForPid=0; - for(Int_t i=2; i<6; i++){ - if(clumap&(1<P(); - AliAODPid *pidObj = track->GetDetPid(); - - Double_t dedx = 0.; - if (pidObj) { - dedx = pidObj->GetITSsignal(); - } - Bool_t mismatch = kTRUE; - Bool_t isSA = kTRUE; - if(track->GetStatus() & AliESDtrack::kTPCin){ - isSA = kFALSE; - if (pidObj) - mom = pidObj->GetTPCmomentum(); - } - for (Int_t j=0; j fRange*sigma) { - p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma; - } else { - p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma; - mismatch=kFALSE; + AliAODTrack *track=(AliAODTrack*)vtrack; + AliAODPid *pidObj = track->GetDetPid(); + if (!pidObj) return -9999.; + Double_t tofTime = 99999; + if (fTuneMConData && ((fTuneMConDataMask & kDetTOF) == kDetTOF) ) tofTime = (Double_t)this->GetTOFsignalTunedOnData((AliVTrack*)vtrack); + else tofTime=pidObj->GetTOFsignal(); + const Double_t expTime=fTOFResponse.GetExpectedSignal((AliVTrack*)vtrack,type); + Double_t sigmaTOFPid[AliPID::kSPECIES]; + pidObj->GetTOFpidResolution(sigmaTOFPid); + AliAODEvent *event=(AliAODEvent*)track->GetAODEvent(); + if (event) { // protection if the user didn't call GetTrack, which sets the internal pointer + AliTOFHeader* tofH=(AliTOFHeader*)event->GetTOFHeader(); + if (tofH && (TMath::Abs(sigmaTOFPid[0]) <= 1.E-16) ) { // new AOD + tofTime -= fTOFResponse.GetStartTime(vtrack->P()); } - - // Check for particles heavier than (AliPID::kSPECIES - 1) - + } else { + AliError("pointer to AliAODEvent not found, please call GetTrack to set it"); + return -9999.; } - if (mismatch) - for (Int_t j=0; j1.e-20) delta=tofTime/expTime; + + return delta; } + //_________________________________________________________________________ -void AliAODpidUtil::MakeTOFPID(const AliAODTrack *track, Double_t *p) const +Float_t AliAODpidUtil::GetNumberOfSigmasTOFold(const AliVParticle *vtrack, AliPID::EParticleType type) const { // - // TOF PID using gaussian response + // Number of sigma implementation for the TOF // - if ((track->GetStatus()&AliESDtrack::kTOFout )==0) return; - if ((track->GetStatus()&AliESDtrack::kTIME )==0) return; - if ((track->GetStatus()&AliESDtrack::kTOFpid )==0) return; + + AliAODTrack *track=(AliAODTrack*)vtrack; - Double_t time[AliPID::kSPECIESN]; - Double_t sigma[AliPID::kSPECIESN]; + Bool_t oldAod=kTRUE; + Double_t sigTOF=0.; AliAODPid *pidObj = track->GetDetPid(); - pidObj->GetIntegratedTimes(time); - - for (Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++) { - sigma[iPart] = fTOFResponse.GetExpectedSigma(track->P(),time[iPart],AliPID::ParticleMass(iPart)); + if (!pidObj) return -999.; + Double_t tofTime = 99999; + if (fTuneMConData && ((fTuneMConDataMask & kDetTOF) == kDetTOF) ) tofTime = (Double_t)this->GetTOFsignalTunedOnData((AliVTrack*)vtrack); + else tofTime=pidObj->GetTOFsignal(); + Double_t expTime=fTOFResponse.GetExpectedSignal((AliVTrack*)vtrack,type); + Double_t sigmaTOFPid[AliPID::kSPECIES]; + pidObj->GetTOFpidResolution(sigmaTOFPid); + AliAODEvent *event=(AliAODEvent*)track->GetAODEvent(); + if (event) { // protection if the user didn't call GetTrack, which sets the internal pointer + AliTOFHeader* tofH=(AliTOFHeader*)event->GetTOFHeader(); + if (tofH && (TMath::Abs(sigmaTOFPid[0]) <= 1.E-16) ) { // new AOD + sigTOF=fTOFResponse.GetExpectedSigma(track->P(),expTime,AliPID::ParticleMassZ(type)); //fTOFResponse is set in InitialiseEvent + tofTime -= fTOFResponse.GetStartTime(vtrack->P()); + oldAod=kFALSE; + } + } else { + AliError("pointer to AliAODEvent not found, please call GetTrack to set it"); + return -996.; } - - AliDebugGeneral("AliESDpid::MakeTOFPID",2, - Form("Expected TOF signals [ps]: %f %f %f %f %f", - time[AliPID::kElectron], - time[AliPID::kMuon], - time[AliPID::kPion], - time[AliPID::kKaon], - time[AliPID::kProton])); - - AliDebugGeneral("AliESDpid::MakeTOFPID",2, - Form("Expected TOF std deviations [ps]: %f %f %f %f %f", - sigma[AliPID::kElectron], - sigma[AliPID::kMuon], - sigma[AliPID::kPion], - sigma[AliPID::kKaon], - sigma[AliPID::kProton] - )); - - Double_t tof = pidObj->GetTOFsignal(); - - Bool_t mismatch = kTRUE; - for (Int_t j=0; j fRange*sig) { - p[j] = TMath::Exp(-0.5*fRange*fRange)/sig; - } else - p[j] = TMath::Exp(-0.5*(tof-time[j])*(tof-time[j])/(sig*sig))/sig; - - // Check the mismatching - Double_t mass = AliPID::ParticleMass(j); - Double_t pm = fTOFResponse.GetMismatchProbability(track->P(),mass); - if (p[j]>pm) mismatch = kFALSE; - - // Check for particles heavier than (AliPID::kSPECIES - 1) - + if (oldAod) { // old AOD + if (type <= AliPID::kProton) { + sigTOF=sigmaTOFPid[type]; + } else return -998.; // light nuclei cannot be supported on old AOD because we don't have timeZero resolution } - - if (mismatch) - for (Int_t j=0; j0) return (tofTime - expTime)/sigTOF; + else return -997.; } -//_________________________________________________________________________ -void AliAODpidUtil::MakeTRDPID(const AliAODTrack *track,Double_t *p) const -{ - ComputeTRDProbability(track, AliPID::kSPECIES, p); - return; -} -