#include <TF1.h>
#include <TSpline.h>
#include <TFile.h>
+#include <TArrayI.h>
+#include <TArrayF.h>
#include <AliVEvent.h>
#include <AliVTrack.h>
#include <AliLog.h>
#include <AliPID.h>
#include <AliOADBContainer.h>
-#include <AliTRDPIDParams.h>
-#include <AliTRDPIDReference.h>
+#include <AliTRDPIDResponseObject.h>
#include <AliTOFPIDParams.h>
#include "AliPIDResponse.h"
+#include "AliDetectorPID.h"
#include "AliCentrality.h"
fITSPIDmethod(kITSTruncMean),
fIsMC(isMC),
fOADBPath(),
+fCustomTPCpidResponse(),
fBeamType("PP"),
fLHCperiod(),
fMCperiodTPC(),
fRecoPassUser(-1),
fRun(0),
fOldRun(0),
-fArrPidResponseMaster(0x0),
-fResolutionCorrection(0x0),
-fTRDPIDParams(0x0),
-fTRDPIDReference(0x0),
+fArrPidResponseMaster(NULL),
+fResolutionCorrection(NULL),
+fOADBvoltageMaps(NULL),
+fTRDPIDResponseObject(NULL),
fTOFtail(1.1),
-fTOFPIDParams(0x0),
-fEMCALPIDParams(0x0),
-fCurrentEvent(0x0),
-fCurrCentrality(0.0)
+fTOFPIDParams(NULL),
+fEMCALPIDParams(NULL),
+fCurrentEvent(NULL),
+fCurrCentrality(0.0),
+fTuneMConData(kFALSE)
{
//
// default ctor
//
- AliLog::SetClassDebugLevel("AliPIDResponse",10);
- AliLog::SetClassDebugLevel("AliESDpid",10);
- AliLog::SetClassDebugLevel("AliAODpidUtil",10);
+ AliLog::SetClassDebugLevel("AliPIDResponse",0);
+ AliLog::SetClassDebugLevel("AliESDpid",0);
+ AliLog::SetClassDebugLevel("AliAODpidUtil",0);
memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
}
// dtor
//
delete fArrPidResponseMaster;
- delete fTRDPIDParams;
- delete fTRDPIDReference;
- if (fTOFPIDParams) delete fTOFPIDParams;
+ delete fTRDPIDResponseObject;
+ delete fTOFPIDParams;
}
//______________________________________________________________________________
fITSPIDmethod(other.fITSPIDmethod),
fIsMC(other.fIsMC),
fOADBPath(other.fOADBPath),
+fCustomTPCpidResponse(other.fCustomTPCpidResponse),
fBeamType("PP"),
fLHCperiod(),
fMCperiodTPC(),
fRecoPassUser(other.fRecoPassUser),
fRun(0),
fOldRun(0),
-fArrPidResponseMaster(0x0),
-fResolutionCorrection(0x0),
-fTRDPIDParams(0x0),
-fTRDPIDReference(0x0),
+fArrPidResponseMaster(NULL),
+fResolutionCorrection(NULL),
+fOADBvoltageMaps(NULL),
+fTRDPIDResponseObject(NULL),
fTOFtail(1.1),
-fTOFPIDParams(0x0),
-fEMCALPIDParams(0x0),
-fCurrentEvent(0x0),
-fCurrCentrality(0.0)
+fTOFPIDParams(NULL),
+fEMCALPIDParams(NULL),
+fCurrentEvent(NULL),
+fCurrCentrality(0.0),
+fTuneMConData(kFALSE)
{
//
// copy ctor
fRange=other.fRange;
fITSPIDmethod=other.fITSPIDmethod;
fOADBPath=other.fOADBPath;
+ fCustomTPCpidResponse=other.fCustomTPCpidResponse;
fIsMC=other.fIsMC;
fBeamType="PP";
fLHCperiod="";
fRecoPassUser=other.fRecoPassUser;
fRun=0;
fOldRun=0;
- fArrPidResponseMaster=0x0;
- fResolutionCorrection=0x0;
- fTRDPIDParams=0x0;
- fTRDPIDReference=0x0;
- fEMCALPIDParams=0x0;
+ fArrPidResponseMaster=NULL;
+ fResolutionCorrection=NULL;
+ fOADBvoltageMaps=NULL;
+ fTRDPIDResponseObject=NULL;
+ fEMCALPIDParams=NULL;
memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
fTOFtail=1.1;
- fTOFPIDParams=0x0;
+ fTOFPIDParams=NULL;
fCurrentEvent=other.fCurrentEvent;
}
return *this;
case kDetITS: return NumberOfSigmasITS(track, type); break;
case kDetTPC: return NumberOfSigmasTPC(track, type); break;
case kDetTOF: return NumberOfSigmasTOF(track, type); break;
-// case kDetTRD: return ComputeTRDProbability(track, type); break;
-// case kDetPHOS: return ComputePHOSProbability(track, type); break;
-// case kDetEMCAL: return NumberOfSigmasEMCAL(track, type); break;
-// case kDetHMPID: return ComputeHMPIDProbability(track, type); break;
+ case kDetEMCAL: return NumberOfSigmasEMCAL(track, type); break;
default: return -999.;
}
}
//______________________________________________________________________________
-Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVTrack *track, AliPID::EParticleType type) const {
+Float_t AliPIDResponse::NumberOfSigmas(EDetector detCode, const AliVParticle *track, AliPID::EParticleType type) const
+{
+ //
+ // NumberOfSigmas for 'detCode'
+ //
+ return NumberOfSigmas((EDetCode)(1<<detCode), track, type);
+}
+
+//______________________________________________________________________________
+Float_t AliPIDResponse::NumberOfSigmasITS(const AliVParticle *vtrack, AliPID::EParticleType type) const
+{
+ //
+ // Calculate the number of sigmas in the ITS
+ //
+
+ AliVTrack *track=(AliVTrack*)vtrack;
+
+ // look for cached value first
+ // only the non SA tracks are cached
+ if ( track->GetDetectorPID() ){
+ return track->GetDetectorPID()->GetNumberOfSigmas(kITS, type);
+ }
+
+ Float_t dEdx=track->GetITSsignal();
+ if (dEdx<=0) return -999.;
+
+ UChar_t clumap=track->GetITSClusterMap();
+ Int_t nPointsForPid=0;
+ for(Int_t i=2; i<6; i++){
+ if(clumap&(1<<i)) ++nPointsForPid;
+ }
+ Float_t mom=track->P();
+
+ //check for ITS standalone tracks
+ Bool_t isSA=kTRUE;
+ if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
+
+ //TODO: in case of the electron, use the SA parametrisation,
+ // this needs to be changed if ITS provides a parametrisation
+ // for electrons also for ITS+TPC tracks
+ return fITSResponse.GetNumberOfSigmas(mom,dEdx,type,nPointsForPid,isSA || (type==AliPID::kElectron));
+}
+
+//______________________________________________________________________________
+Float_t AliPIDResponse::NumberOfSigmasTPC(const AliVParticle *vtrack, AliPID::EParticleType type) const
+{
+ //
+ // Calculate the number of sigmas in the TPC
+ //
+
+ AliVTrack *track=(AliVTrack*)vtrack;
+
+ // look for cached value first
+ if (track->GetDetectorPID()){
+ return track->GetDetectorPID()->GetNumberOfSigmas(kTPC, type);
+ }
+
+ Double_t mom = track->GetTPCmomentum();
+ Double_t sig = track->GetTPCsignal();
+ UInt_t sigN = track->GetTPCsignalN();
+
+ Double_t nSigma = -999.;
+ if (sigN>0) nSigma=fTPCResponse.GetNumberOfSigmas(mom,sig,sigN,type);
+
+ return nSigma;
+}
+
+//______________________________________________________________________________
+Float_t AliPIDResponse::NumberOfSigmasTPC( const AliVParticle *vtrack,
+ AliPID::EParticleType type,
+ AliTPCPIDResponse::ETPCdEdxSource dedxSource)
+{
+ //get number of sigmas according the selected TPC gain configuration scenario
+ const AliVTrack *track=static_cast<const AliVTrack*>(vtrack);
+
+ Float_t nSigma=fTPCResponse.GetNumberOfSigmas(track, type, dedxSource);
+
+ return nSigma;
+}
+
+//______________________________________________________________________________
+Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type) const
+{
+ //
+ // Calculate the number of sigmas in the EMCAL
+ //
+
+ AliVTrack *track=(AliVTrack*)vtrack;
+ // look for cached value first
+ if (track->GetDetectorPID()){
+ return track->GetDetectorPID()->GetNumberOfSigmas(kEMCAL, type);
+ }
+
AliVCluster *matchedClus = NULL;
Double_t mom = -1.;
}
//______________________________________________________________________________
-Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVTrack *track, AliPID::EParticleType type, Double_t &eop, Double_t showershape[4]) const {
+Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &eop, Double_t showershape[4]) const {
+ AliVTrack *track=(AliVTrack*)vtrack;
+
AliVCluster *matchedClus = NULL;
Double_t mom = -1.;
Double_t pt = -1.;
Double_t EovP = -1.;
Double_t fClsE = -1.;
+
+ // initialize eop and shower shape parameters
+ eop = -1.;
+ for(Int_t i = 0; i < 4; i++){
+ showershape[i] = -1.;
+ }
Int_t nMatchClus = -1;
Int_t charge = 0;
switch (detCode){
case kDetITS: return ComputeITSProbability(track, nSpecies, p); break;
case kDetTPC: return ComputeTPCProbability(track, nSpecies, p); break;
- case kDetTOF: return ComputeTOFProbability(track, nSpecies, p); break;
case kDetTRD: return ComputeTRDProbability(track, nSpecies, p); break;
+ case kDetTOF: return ComputeTOFProbability(track, nSpecies, p); break;
case kDetPHOS: return ComputePHOSProbability(track, nSpecies, p); break;
case kDetEMCAL: return ComputeEMCALProbability(track, nSpecies, p); break;
case kDetHMPID: return ComputeHMPIDProbability(track, nSpecies, p); break;
}
}
+//______________________________________________________________________________
+AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability (EDetector detCode, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
+{
+ //
+ // Compute PID response of 'detCode'
+ //
+
+ return ComputePIDProbability((EDetCode)(1<<detCode),track,nSpecies,p);
+}
+
//______________________________________________________________________________
AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeITSProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
{
// Compute PID response for the ITS
//
+ // look for cached value first
+ // only the non SA tracks are cached
+ if (track->GetDetectorPID()){
+ return track->GetDetectorPID()->GetRawProbability(kITS, p, nSpecies);
+ }
+
// set flat distribution (no decision)
for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
if ((track->GetStatus()&AliVTrack::kITSin)==0 &&
(track->GetStatus()&AliVTrack::kITSout)==0) return kDetNoSignal;
+
+ //check for ITS standalone tracks
+ Bool_t isSA=kTRUE;
+ if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
Double_t mom=track->P();
Double_t dedx=track->GetITSsignal();
- Bool_t isSA=kTRUE;
Double_t momITS=mom;
- ULong_t trStatus=track->GetStatus();
- if(trStatus&AliVTrack::kTPCin) isSA=kFALSE;
UChar_t clumap=track->GetITSClusterMap();
Int_t nPointsForPid=0;
for(Int_t i=2; i<6; i++){
Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
for (Int_t j=0; j<AliPID::kSPECIES; j++) {
- Double_t mass=AliPID::ParticleMass(j);//GeV/c^2
- Double_t bethe=fITSResponse.Bethe(momITS,mass);
- Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA);
+ Double_t mass=AliPID::ParticleMassZ(j);//GeV/c^2
+ const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(j),2.);
+ Double_t bethe=fITSResponse.Bethe(momITS,mass)*chargeFactor;
+ //TODO: in case of the electron, use the SA parametrisation,
+ // this needs to be changed if ITS provides a parametrisation
+ // for electrons also for ITS+TPC tracks
+ Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA || (j==(Int_t)AliPID::kElectron));
if (TMath::Abs(dedx-bethe) > fRange*sigma) {
p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
} else {
//
// Compute PID response for the TPC
//
-
+
+ // look for cached value first
+ if (track->GetDetectorPID()){
+ return track->GetDetectorPID()->GetRawProbability(kTPC, p, nSpecies);
+ }
+
// set flat distribution (no decision)
for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
Double_t dedx=track->GetTPCsignal();
Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
- for (Int_t j=0; j<AliPID::kSPECIES; j++) {
+ if(fTuneMConData) dedx = this->GetTPCsignalTunedOnData(track);
+
+ for (Int_t j=0; j<AliPID::kSPECIESC; j++) {
AliPID::EParticleType type=AliPID::EParticleType(j);
Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
Double_t sigma=fTPCResponse.GetExpectedSigma(mom,track->GetTPCsignalN(),type);
p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
mismatch=kFALSE;
}
-
- // TODO: Light nuclei, also in TPC pid response
-
- // Check for particles heavier than (AliPID::kSPECIES - 1)
-// if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
-
}
if (mismatch){
//
// Compute PID response for the
//
-
+
+ // look for cached value first
+ if (track->GetDetectorPID()){
+ return track->GetDetectorPID()->GetRawProbability(kTOF, p, nSpecies);
+ }
+
Double_t meanCorrFactor = 0.11/fTOFtail; // Correction factor on the mean because of the tail (should be ~ 0.1 with tail = 1.1)
// set flat distribution (no decision)
if ((track->GetStatus()&AliVTrack::kTOFout)==0) return kDetNoSignal;
if ((track->GetStatus()&AliVTrack::kTIME)==0) return kDetNoSignal;
- Double_t time[AliPID::kSPECIESN];
- track->GetIntegratedTimes(time);
-
- Double_t sigma[AliPID::kSPECIES];
- for (Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++) {
- sigma[iPart] = fTOFResponse.GetExpectedSigma(track->P(),time[iPart],AliPID::ParticleMass(iPart));
- }
-
Bool_t mismatch = kTRUE/*, heavy = kTRUE*/;
- for (Int_t j=0; j<AliPID::kSPECIES; j++) {
+ for (Int_t j=0; j<AliPID::kSPECIESC; j++) {
AliPID::EParticleType type=AliPID::EParticleType(j);
Double_t nsigmas=NumberOfSigmasTOF(track,type) + meanCorrFactor;
- Double_t sig = sigma[j];
+ Double_t expTime = fTOFResponse.GetExpectedSignal(track,type);
+ Double_t sig = fTOFResponse.GetExpectedSigma(track->P(),expTime,AliPID::ParticleMassZ(type));
if (TMath::Abs(nsigmas) > (fRange+2)) {
if(nsigmas < fTOFtail)
p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
p[j] = TMath::Exp(-(nsigmas - fTOFtail*0.5)*fTOFtail)/sig;
}
- /* OLD Gaussian shape
- if (TMath::Abs(nsigmas) > (fRange+2)) {
- p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
- } else
- p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
- */
-
if (TMath::Abs(nsigmas)<5.){
Double_t nsigmasTPC=NumberOfSigmasTPC(track,type);
if (TMath::Abs(nsigmasTPC)<5.) mismatch=kFALSE;
//
// Compute PID response for the
//
-
+
+ // look for cached value first
+ if (track->GetDetectorPID()){
+ return track->GetDetectorPID()->GetRawProbability(kTRD, p, nSpecies);
+ }
+
// set flat distribution (no decision)
for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
if((track->GetStatus()&AliVTrack::kTRDout)==0) return kDetNoSignal;
- Float_t mom[6];
- Double_t dedx[48]; // Allocate space for the maximum number of TRD slices
+ Float_t mom[6]={0.};
+ Double_t dedx[48]={0.}; // Allocate space for the maximum number of TRD slices
Int_t nslices = fTRDslicesForPID[1] - fTRDslicesForPID[0] + 1;
AliDebug(1, Form("First Slice: %d, Last Slice: %d, Number of slices: %d", fTRDslicesForPID[0], fTRDslicesForPID[1], nslices));
for(UInt_t ilayer = 0; ilayer < 6; ilayer++){
//
// Compute PID response for the EMCAL
//
+
+ // look for cached value first
+ if (track->GetDetectorPID()){
+ return track->GetDetectorPID()->GetRawProbability(kEMCAL, p, nSpecies);
+ }
+
+ for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
AliVCluster *matchedClus = NULL;
nMatchClus = track->GetEMCALcluster();
if(nMatchClus > -1){
-
+
mom = track->P();
pt = track->Pt();
charge = track->Charge();
matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
if(matchedClus){
-
- // matched cluster is EMCAL
- if(matchedClus->IsEMCAL()){
-
+
+ // matched cluster is EMCAL
+ if(matchedClus->IsEMCAL()){
+
fClsE = matchedClus->E();
EovP = fClsE/mom;
// compute the probabilities
- if( 999 != fEMCALResponse.ComputeEMCALProbability(pt,EovP,charge,p)){
-
- // in case everything is OK
- return kDetPidOk;
-
+ if(fEMCALResponse.ComputeEMCALProbability(nSpecies,pt,EovP,charge,p)){
+
+ // in case everything is OK
+ return kDetPidOk;
+ }
}
}
}
- }
// in all other cases set flat distribution (no decision)
- for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
+ for (Int_t j=0; j<nSpecies; j++) p[j] = 1./nSpecies;
return kDetNoSignal;
}
//
// Compute PID response for the PHOS
//
-
+
+ // look for cached value first
+// if (track->GetDetectorPID()){
+// return track->GetDetectorPID()->GetRawProbability(kPHOS, p, nSpecies);
+// }
+
// set flat distribution (no decision)
for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
return kDetNoSignal;
// Compute PID response for the HMPID
//
+
+ // look for cached value first
+ if (track->GetDetectorPID()){
+ return track->GetDetectorPID()->GetRawProbability(kHMPID, p, nSpecies);
+ }
+
// set flat distribution (no decision)
for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
if((track->GetStatus()&AliVTrack::kHMPIDpid)==0) return kDetNoSignal;
}
//______________________________________________________________________________
-void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass)
+void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass, Int_t run)
{
//
// Apply settings for the current event
//
fRecoPass=pass;
- fCurrentEvent=0x0;
+ fCurrentEvent=NULL;
if (!event) return;
fCurrentEvent=event;
- fRun=event->GetRunNumber();
+ if (run>0) fRun=run;
+ else fRun=event->GetRunNumber();
if (fRun!=fOldRun){
ExecNewRun();
SetTOFPidResponseMaster();
InitializeTOFResponse();
+
+ if (fCurrentEvent) fTPCResponse.SetMagField(fCurrentEvent->GetMagneticField());
}
//_____________________________________________________
fBeamType="PP";
- TPRegexp reg(".*(LHC11[a-z]+[0-9]+[a-z_]*)/.*");
+ TPRegexp reg(".*(LHC1[1-2][a-z]+[0-9]+[a-z_]*)/.*");
+ TPRegexp reg12a17(".*(LHC12a17[a-z]+)/.*");
+
//find the period by run number (UGLY, but not stored in ESD and AOD... )
if (fRun>=114737&&fRun<=117223) { fLHCperiod="LHC10B"; fMCperiodTPC="LHC10D1"; }
else if (fRun>=118503&&fRun<=121040) { fLHCperiod="LHC10C"; fMCperiodTPC="LHC10D1"; }
//TODO: periods 11B, 11C are not yet treated assume 11d for the moment
else if (fRun>=148531&&fRun<=155384) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
else if (fRun>=156477&&fRun<=159635) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
+ // also for 11e,f use 11d
+ else if (fRun>=160676&&fRun<=162740) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
+ else if (fRun>=162933&&fRun<=165746) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
+
else if (fRun>=166529) {
fLHCperiod="LHC11H";
fMCperiodTPC="LHC11A10";
fBeamType="PBPB";
+ if (reg12a17.MatchB(fCurrentFile)) fMCperiodTPC="LHC12A17";
}
//exception new pp MC productions from 2011
if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11B2";
+ // exception for 11f1
+ if (fCurrentFile.Contains("LHC11f1/")) fMCperiodTPC="LHC11F1";
}
//______________________________________________________________________________
{
//
// Load the TPC pid response functions from the OADB
+ // Load the TPC voltage maps from OADB
//
//don't load twice for the moment
if (fArrPidResponseMaster) return;
//reset the PID response functions
delete fArrPidResponseMaster;
- fArrPidResponseMaster=0x0;
+ fArrPidResponseMaster=NULL;
TString fileName(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data()));
+ TFile *f=NULL;
+ if (!fCustomTPCpidResponse.IsNull()) fileName=fCustomTPCpidResponse;
- TFile *f=TFile::Open(fileName.Data());
+ TString fileNamePIDresponse(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data()));
+ f=TFile::Open(fileNamePIDresponse.Data());
if (f && f->IsOpen() && !f->IsZombie()){
fArrPidResponseMaster=dynamic_cast<TObjArray*>(f->Get("TPCPIDResponse"));
}
delete f;
+
+ TString fileNameVoltageMaps(Form("%s/COMMON/PID/data/TPCvoltageSettings.root", fOADBPath.Data()));
+ f=TFile::Open(fileNameVoltageMaps.Data());
+ if (f && f->IsOpen() && !f->IsZombie()){
+ fOADBvoltageMaps=dynamic_cast<AliOADBContainer*>(f->Get("TPCvoltageSettings"));
+ }
+ delete f;
if (!fArrPidResponseMaster){
- AliFatal(Form("Could not retrieve the TPC pid response from: %s",fileName.Data()));
+ AliFatal(Form("Could not retrieve the TPC pid response from: %s",fileNamePIDresponse.Data()));
return;
}
fArrPidResponseMaster->SetOwner();
+
+ if (!fOADBvoltageMaps)
+ {
+ AliFatal(Form("Could not retrieve the TPC voltage maps from: %s",fileNameVoltageMaps.Data()));
+ }
+ fArrPidResponseMaster->SetOwner();
}
//______________________________________________________________________________
TString datatype="DATA";
//in case of mc fRecoPass is per default 1
if (fIsMC) {
- datatype="MC";
- fRecoPass=1;
+ if(!fTuneMConData) datatype="MC";
+ fRecoPass=1;
}
//
//reset old splines
//
- for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
- fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,0x0);
- }
-
+ fTPCResponse.ResetSplines();
+
+ // period
+ TString period=fLHCperiod;
+ if (fIsMC && !fTuneMConData) period=fMCperiodTPC;
+
+ AliInfo(Form("Searching splines for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
+ Bool_t found=kFALSE;
//
//set the new PID splines
//
- TString period=fLHCperiod;
if (fArrPidResponseMaster){
- TObject *grAll=0x0;
+ Int_t recopass = fRecoPass;
+ if(fTuneMConData) recopass = fRecoPassUser;
//for MC don't use period information
-// if (fIsMC) period="[A-Z0-9]*";
+ //if (fIsMC) period="[A-Z0-9]*";
//for MC use MC period information
- if (fIsMC) period=fMCperiodTPC;
-//pattern for the default entry (valid for all particles)
- TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
-
- //loop over entries and filter them
- for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp){
- TObject *responseFunction=fArrPidResponseMaster->At(iresp);
- if (responseFunction==0x0) continue;
- TString responseName=responseFunction->GetName();
-
- if (!reg.MatchB(responseName)) continue;
-
- TObjArray *arr=reg.MatchS(responseName);
- TString particleName=arr->At(1)->GetName();
- delete arr;
- if (particleName.IsNull()) continue;
- if (particleName=="ALL") grAll=responseFunction;
- else {
- //find particle id
- for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
- TString particle=AliPID::ParticleName(ispec);
- particle.ToUpper();
- if ( particle == particleName ){
- fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,responseFunction);
- fTPCResponse.SetUseDatabase(kTRUE);
- AliInfo(Form("Adding graph: %d - %s",ispec,responseFunction->GetName()));
- break;
+ //pattern for the default entry (valid for all particles)
+ TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN(_*)([A-Z1-9]*)",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
+
+ //find particle id ang gain scenario
+ for (Int_t igainScenario=0; igainScenario<AliTPCPIDResponse::fgkNumberOfGainScenarios; igainScenario++)
+ {
+ TObject *grAll=NULL;
+ TString gainScenario = AliTPCPIDResponse::GainScenarioName(igainScenario);
+ gainScenario.ToUpper();
+ //loop over entries and filter them
+ for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp)
+ {
+ TObject *responseFunction=fArrPidResponseMaster->At(iresp);
+ if (responseFunction==NULL) continue;
+ TString responseName=responseFunction->GetName();
+
+ if (!reg.MatchB(responseName)) continue;
+
+ TObjArray *arr=reg.MatchS(responseName); if (!arr) continue;
+ TObject* tmp=NULL;
+ tmp=arr->At(1); if (!tmp) continue;
+ TString particleName=tmp->GetName();
+ tmp=arr->At(3); if (!tmp) continue;
+ TString gainScenarioName=tmp->GetName();
+ delete arr;
+ if (particleName.IsNull()) continue;
+ if (!grAll && particleName=="ALL" && gainScenarioName==gainScenario) grAll=responseFunction;
+ else
+ {
+ for (Int_t ispec=0; ispec<(AliTPCPIDResponse::fgkNumberOfParticleSpecies); ++ispec)
+ {
+ TString particle=AliPID::ParticleName(ispec);
+ particle.ToUpper();
+ //std::cout<<responseName<<" "<<particle<<" "<<particleName<<" "<<gainScenario<<" "<<gainScenarioName<<std::endl;
+ if ( particle == particleName && gainScenario == gainScenarioName )
+ {
+ fTPCResponse.SetResponseFunction( responseFunction,
+ (AliPID::EParticleType)ispec,
+ (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
+ fTPCResponse.SetUseDatabase(kTRUE);
+ AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,responseFunction->GetName()));
+ found=kTRUE;
+ // overwrite default with proton spline (for light nuclei)
+ if (ispec==AliPID::kProton) grAll=responseFunction;
+ break;
+ }
}
}
}
- }
-
- //set default response function to all particles which don't have a specific one
- if (grAll){
- for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
- if (!fTPCResponse.GetResponseFunction((AliPID::EParticleType)ispec)){
- fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,grAll);
- AliInfo(Form("Adding graph: %d - %s",ispec,grAll->GetName()));
+ if (grAll)
+ {
+ for (Int_t ispec=0; ispec<(AliTPCPIDResponse::fgkNumberOfParticleSpecies); ++ispec)
+ {
+ if (!fTPCResponse.GetResponseFunction( (AliPID::EParticleType)ispec,
+ (AliTPCPIDResponse::ETPCgainScenario)igainScenario))
+ {
+ fTPCResponse.SetResponseFunction( grAll,
+ (AliPID::EParticleType)ispec,
+ (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
+ fTPCResponse.SetUseDatabase(kTRUE);
+ AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,grAll->GetName()));
+ found=kTRUE;
+ }
}
}
}
}
-
+ else AliInfo("no fArrPidResponseMaster");
+
+ if (!found){
+ AliError(Form("No splines found for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
+ }
+
//
// Setup resolution parametrisation
//
fResolutionCorrection=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_SIGMA",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
if (fResolutionCorrection) AliInfo(Form("Setting multiplicity correction function: %s",fResolutionCorrection->GetName()));
+
+ //read in the voltage map
+ TVectorF* gsm = dynamic_cast<TVectorF*>(fOADBvoltageMaps->GetObject(fRun));
+ if (gsm)
+ {
+ fTPCResponse.SetVoltageMap(*gsm);
+ TString vals;
+ AliInfo(Form("Reading the voltage map for run %d\n",fRun));
+ vals="IROC A: "; for (Int_t i=0; i<18; i++){vals+=Form("%.2f ",(*gsm)[i]);}
+ AliInfo(vals.Data());
+ vals="IROC C: "; for (Int_t i=18; i<36; i++){vals+=Form("%.2f ",(*gsm)[i]);}
+ AliInfo(vals.Data());
+ vals="OROC A: "; for (Int_t i=36; i<54; i++){vals+=Form("%.2f ",(*gsm)[i]);}
+ AliInfo(vals.Data());
+ vals="OROC C: "; for (Int_t i=54; i<72; i++){vals+=Form("%.2f ",(*gsm)[i]);}
+ AliInfo(vals.Data());
+ }
+ else AliInfo("no voltage map, ideal default assumed");
}
//______________________________________________________________________________
//
// Load the TRD pid params and references from the OADB
//
- if(fTRDPIDParams) return;
+ if(fTRDPIDResponseObject) return;
AliOADBContainer contParams("contParams");
- Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDParams.root", fOADBPath.Data()), "AliTRDPIDParams");
- if(statusPars){
- AliError("Failed initializing PID Params from OADB");
+ Int_t statusResponse = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()), "AliTRDPIDResponseObject");
+ if(statusResponse){
+ AliError("Failed initializing PID Response Object from OADB");
} else {
- AliInfo(Form("Loading TRD Params from %s/COMMON/PID/data/TRDPIDParams.root", fOADBPath.Data()));
- fTRDPIDParams = dynamic_cast<AliTRDPIDParams *>(contParams.GetObject(fRun));
- if(!fTRDPIDParams){
- AliError(Form("TRD Params not found in run %d", fRun));
+ AliInfo(Form("Loading TRD Response from %s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()));
+ fTRDPIDResponseObject = dynamic_cast<AliTRDPIDResponseObject *>(contParams.GetObject(fRun));
+ if(!fTRDPIDResponseObject){
+ AliError(Form("TRD Response not found in run %d", fRun));
}
}
-
+ /*
AliOADBContainer contRefs("contRefs");
Int_t statusRefs = contRefs.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDReferenceLQ1D.root", fOADBPath.Data()), "AliTRDPIDReference");
if(statusRefs){
if(!fTRDPIDReference){
AliError(Form("TRD References not found in OADB Container for run %d", fRun));
}
- }
+ }
+ */
}
//______________________________________________________________________________
//
// Set PID Params and references to the TRD PID response
//
- fTRDResponse.SetPIDParams(fTRDPIDParams);
- fTRDResponse.Load(fTRDPIDReference);
- if(fLHCperiod == "LHC10b" || fLHCperiod == "LHC10c" || fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
+ fTRDResponse.SetPIDResponseObject(fTRDPIDResponseObject);
+ SetTRDPIDmethod();
+}
+
+void AliPIDResponse::SetTRDPIDmethod(AliTRDPIDResponse::ETRDPIDMethod method){
+
+ fTRDResponse.SetPIDmethod(method);
+ if(fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
+ // backward compatibility for setting with 8 slices
fTRDslicesForPID[0] = 0;
fTRDslicesForPID[1] = 7;
}
+ else{
+ if(method==AliTRDPIDResponse::kLQ1D){
+ fTRDslicesForPID[0] = 0; // first Slice contains normalized dEdx
+ fTRDslicesForPID[1] = 0;
+ }
+ if(method==AliTRDPIDResponse::kLQ2D){
+ fTRDslicesForPID[0] = 1;
+ fTRDslicesForPID[1] = 7;
+ }
+ }
+ AliDebug(1,Form("Slice Range set to %d - %d",fTRDslicesForPID[0],fTRDslicesForPID[1]));
}
//______________________________________________________________________________
//
// Load the TOF pid params from the OADB
//
+
+ if (fTOFPIDParams) delete fTOFPIDParams;
+ fTOFPIDParams=NULL;
+
TFile *oadbf = new TFile(Form("%s/COMMON/PID/data/TOFPIDParams.root",fOADBPath.Data()));
- if (oadbf->IsOpen()) {
+ if (oadbf && oadbf->IsOpen()) {
AliInfo(Form("Loading TOF Params from %s/COMMON/PID/data/TOFPIDParams.root", fOADBPath.Data()));
AliOADBContainer *oadbc = (AliOADBContainer *)oadbf->Get("TOFoadb");
- if (fTOFPIDParams) delete fTOFPIDParams;
- fTOFPIDParams = dynamic_cast<AliTOFPIDParams *>(oadbc->GetObject(fRun));
+ if (oadbc) fTOFPIDParams = dynamic_cast<AliTOFPIDParams *>(oadbc->GetObject(fRun,"TOFparams"));
oadbf->Close();
delete oadbc;
- } else {
- AliError(Form("TOFPIDParams.root not found in %s/COMMON/PID/data !!",fOADBPath.Data()));
}
delete oadbf;
- }
+ if (!fTOFPIDParams) AliFatal("TOFPIDParams could not be retrieved");
+}
//______________________________________________________________________________
void AliPIDResponse::InitializeTOFResponse(){
//
// Set PID Params to the TOF PID response
- //
- TString stMethod[4]={"kFILL_T0","kTOF_T0","kT0_T0","kBest_T0"};
+ //
+
+ AliInfo("TOF PID Params loaded from OADB");
+ AliInfo(Form(" TOF resolution %5.2f [ps]",fTOFPIDParams->GetTOFresolution()));
+ AliInfo(Form(" StartTime method %d",fTOFPIDParams->GetStartTimeMethod()));
+ AliInfo(Form(" TOF res. mom. params: %5.2f %5.2f %5.2f %5.2f",
+ fTOFPIDParams->GetSigParams(0),fTOFPIDParams->GetSigParams(1),fTOFPIDParams->GetSigParams(2),fTOFPIDParams->GetSigParams(3)));
+
for (Int_t i=0;i<4;i++) {
fTOFResponse.SetTrackParameter(i,fTOFPIDParams->GetSigParams(i));
}
fEMCALResponse.SetPIDParams(fEMCALPIDParams);
}
+
+//_____________________________________________________
+void AliPIDResponse::FillTrackDetectorPID()
+{
+ //
+ // create detector PID information and setup the transient pointer in the track
+ //
+
+ if (!fCurrentEvent) return;
+
+ //TODO: which particles to include? See also the loops below...
+ Double_t values[AliPID::kSPECIESC]={0};
+
+ for (Int_t itrack=0; itrack<fCurrentEvent->GetNumberOfTracks(); ++itrack){
+ AliVTrack *track=dynamic_cast<AliVTrack*>(fCurrentEvent->GetTrack(itrack));
+ if (!track) continue;
+
+ AliDetectorPID *detPID=new AliDetectorPID;
+ for (Int_t idet=0; idet<kNdetectors; ++idet){
+
+ //nsigmas
+ for (Int_t ipart=0; ipart<AliPID::kSPECIESC; ++ipart)
+ values[ipart]=NumberOfSigmas((EDetector)idet,track,(AliPID::EParticleType)ipart);
+ detPID->SetNumberOfSigmas((EDetector)idet, values, (Int_t)AliPID::kSPECIESC);
+
+ //probabilities
+ EDetPidStatus status=ComputePIDProbability((EDetector)idet,track,AliPID::kSPECIESC,values);
+ detPID->SetRawProbability((EDetector)idet, values, (Int_t)AliPID::kSPECIESC, status);
+ }
+
+ track->SetDetectorPID(detPID);
+ }
+}
+
//_________________________________________________________________________
void AliPIDResponse::SetTOFResponse(AliVEvent *vevent,EStartTimeType_t option){
//
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff