o Add possibility to correct for the eta dependence of the TPC signal when calculating nSigma
o Add Eta correction maps
o Updated TPC splines: optimised for eta corrected signal, better performance at low p
o Currently supported periods 10b,c,d,e pass2 11a pass1,2
(Benjamin Hess, Jens Wiechula)
fOldRun(0),
fRecoPass(0),
fIsTunedOnData(kFALSE),
-fRecoPassTuned(0)
+fRecoPassTuned(0),
+fUseTPCEtaCorrection(kFALSE)//TODO: In future, default kTRUE
{
//
// Dummy constructor
fOldRun(0),
fRecoPass(0),
fIsTunedOnData(kFALSE),
-fRecoPassTuned(0)
+fRecoPassTuned(0),
+fUseTPCEtaCorrection(kFALSE)//TODO: In future, default kTRUE
{
//
// Default constructor
fOldRun=fRun;
}
+ fPIDResponse->SetUseTPCEtaCorrection(fUseTPCEtaCorrection);
fPIDResponse->InitialiseEvent(event,fRecoPass);
AliESDpid *pidresp = dynamic_cast<AliESDpid*>(fPIDResponse);
if(pidresp && AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()){
void SetOADBPath(const char* path) {fOADBPath=path;}
const char* GetOADBPath() const { return fOADBPath.Data(); }
void SetTuneOnData(Bool_t flag,Int_t recopass){fIsTunedOnData=flag;fRecoPassTuned=recopass;};
+
+ void SetUseTPCEtaCorrection(Bool_t useTPCEtaCorrection) { fUseTPCEtaCorrection = useTPCEtaCorrection; };
+ Bool_t UseTPCEtaCorrection() const { return fUseTPCEtaCorrection; };
void SetSpecialDetectorResponse(const char* det) { fSpecialDetResponse=det; }
Bool_t fIsTunedOnData; // flag to tune MC on data (dE/dx)
Int_t fRecoPassTuned; // Reco pass tuned on data for MC
+ Bool_t fUseTPCEtaCorrection; //! Use TPC eta correction
+
//
void SetRecoInfo();
AliAnalysisTaskPIDResponse(const AliAnalysisTaskPIDResponse &other);
AliAnalysisTaskPIDResponse& operator=(const AliAnalysisTaskPIDResponse &other);
- ClassDef(AliAnalysisTaskPIDResponse,3) // Task to properly set the PID response functions of all detectors
+ ClassDef(AliAnalysisTaskPIDResponse,4) // Task to properly set the PID response functions of all detectors
};
#endif
AliAnalysisTask *AddTaskPIDResponse(Bool_t isMC=kFALSE, Bool_t autoMCesd=kTRUE,
Bool_t tuneOnData=kFALSE, Int_t recoPass=2,
- Bool_t cachePID=kFALSE, TString detResponse="")
+ Bool_t cachePID=kFALSE, TString detResponse="",
+ Bool_t useTPCEtaCorrection = kFALSE)
{
// Macro to connect a centrality selection task to an existing analysis manager.
AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
if(isMC&&tuneOnData) pidTask->SetTuneOnData(kTRUE,recoPass);
pidTask->SetCachePID(cachePID);
pidTask->SetSpecialDetectorResponse(detResponse);
+ pidTask->SetUseTPCEtaCorrection(useTPCEtaCorrection);
mgr->AddTask(pidTask);
// AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("PIDResponseQA",
#include <TObjArray.h>
#include <TPRegexp.h>
#include <TF1.h>
+#include <TH2D.h>
#include <TSpline.h>
#include <TFile.h>
#include <TArrayI.h>
#include <TArrayF.h>
+#include <TLinearFitter.h>
#include <AliVEvent.h>
#include <AliVTrack.h>
fArrPidResponseMaster(NULL),
fResolutionCorrection(NULL),
fOADBvoltageMaps(NULL),
+fUseTPCEtaCorrection(kFALSE),//TODO: In future, default kTRUE
fTRDPIDResponseObject(NULL),
fTOFtail(1.1),
fTOFPIDParams(NULL),
fArrPidResponseMaster(NULL),
fResolutionCorrection(NULL),
fOADBvoltageMaps(NULL),
+fUseTPCEtaCorrection(other.fUseTPCEtaCorrection),
fTRDPIDResponseObject(NULL),
fTOFtail(1.1),
fTOFPIDParams(NULL),
fArrPidResponseMaster=NULL;
fResolutionCorrection=NULL;
fOADBvoltageMaps=NULL;
+ fUseTPCEtaCorrection=other.fUseTPCEtaCorrection;
fTRDPIDResponseObject=NULL;
fEMCALPIDParams=NULL;
fTOFtail=1.1;
//______________________________________________________________________________
Float_t AliPIDResponse::NumberOfSigmasTPC( const AliVParticle *vtrack,
AliPID::EParticleType type,
- AliTPCPIDResponse::ETPCdEdxSource dedxSource)
+ AliTPCPIDResponse::ETPCdEdxSource dedxSource) const
{
//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);
+ Float_t nSigma=fTPCResponse.GetNumberOfSigmas(track, type, dedxSource, fUseTPCEtaCorrection);
return nSigma;
}
//______________________________________________________________________________
-AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability (EDetCode detCode, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
+AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability (EDetCode detCode, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
{
//
// Compute PID response of 'detCode'
SetTPCPidResponseMaster();
SetTPCParametrisation();
+ SetTPCEtaMaps();
SetTRDPidResponseMaster();
InitializeTRDResponse();
//
}
+
+//______________________________________________________________________________
+void AliPIDResponse::AddPointToHyperplane(TH2D* h, TLinearFitter* linExtrapolation, Int_t binX, Int_t binY)
+{
+ if (h->GetBinContent(binX, binY) <= 1e-4)
+ return; // Reject bins without content (within some numerical precision) or with strange content
+
+ Double_t coord[2] = {0, 0};
+ coord[0] = h->GetXaxis()->GetBinCenter(binX);
+ coord[1] = h->GetYaxis()->GetBinCenter(binY);
+ Double_t binError = h->GetBinError(binX, binY);
+ if (binError <= 0) {
+ binError = 1000; // Should not happen because bins without content are rejected for the map (TH2D* h)
+ printf("ERROR: This should never happen: Trying to add bin in addPointToHyperplane with error not set....\n");
+ }
+ linExtrapolation->AddPoint(coord, h->GetBinContent(binX, binY, binError));
+}
+
+
+//______________________________________________________________________________
+TH2D* AliPIDResponse::RefineHistoViaLinearInterpolation(TH2D* h, Double_t refineFactorX, Double_t refineFactorY)
+{
+ if (!h)
+ return 0x0;
+
+ // Interpolate to finer map
+ TLinearFitter* linExtrapolation = new TLinearFitter(2, "hyp2", "");
+
+ Double_t upperMapBoundY = h->GetYaxis()->GetBinUpEdge(h->GetYaxis()->GetNbins());
+ Double_t lowerMapBoundY = h->GetYaxis()->GetBinLowEdge(1);
+ Int_t nBinsX = 20;
+ // Binning was find to yield good results, if 40 bins are chosen for the range 0.0016 to 0.02. For the new variable range,
+ // scale the number of bins correspondingly
+ Int_t nBinsY = TMath::Nint((upperMapBoundY - lowerMapBoundY) / (0.02 - lowerMapBoundY) * 40);
+ Int_t nBinsXrefined = nBinsX * refineFactorX;
+ Int_t nBinsYrefined = nBinsY * refineFactorY;
+
+ TH2D* hRefined = new TH2D(Form("%s_refined", h->GetName()), Form("%s (refined)", h->GetTitle()),
+ nBinsXrefined, h->GetXaxis()->GetBinLowEdge(1), h->GetXaxis()->GetBinUpEdge(h->GetXaxis()->GetNbins()),
+ nBinsYrefined, lowerMapBoundY, upperMapBoundY);
+
+ for (Int_t binX = 1; binX <= nBinsXrefined; binX++) {
+ for (Int_t binY = 1; binY <= nBinsYrefined; binY++) {
+
+ hRefined->SetBinContent(binX, binY, 1); // Default value is 1
+
+ Double_t centerX = hRefined->GetXaxis()->GetBinCenter(binX);
+ Double_t centerY = hRefined->GetYaxis()->GetBinCenter(binY);
+
+ /*TODO NOW NOW
+ linExtrapolation->ClearPoints();
+
+ // For interpolation: Just take the corresponding bin from the old histo.
+ // For extrapolation: take the last available bin from the old histo.
+ // If the boundaries are to be skipped, also skip the corresponding bins
+ Int_t oldBinX = h->GetXaxis()->FindBin(centerX);
+ if (oldBinX < 1)
+ oldBinX = 1;
+ if (oldBinX > nBinsX)
+ oldBinX = nBinsX;
+
+ Int_t oldBinY = h->GetYaxis()->FindBin(centerY);
+ if (oldBinY < 1)
+ oldBinY = 1;
+ if (oldBinY > nBinsY)
+ oldBinY = nBinsY;
+
+ // Neighbours left column
+ if (oldBinX >= 2) {
+ if (oldBinY >= 2) {
+ AddPointToHyperplane(h, linExtrapolation, oldBinX - 1, oldBinY - 1);
+ }
+
+ AddPointToHyperplane(h, linExtrapolation, oldBinX - 1, oldBinY);
+
+ if (oldBinY < nBinsY) {
+ AddPointToHyperplane(h, linExtrapolation, oldBinX - 1, oldBinY + 1);
+ }
+ }
+
+ // Neighbours (and point itself) same column
+ if (oldBinY >= 2) {
+ AddPointToHyperplane(h, linExtrapolation, oldBinX, oldBinY - 1);
+ }
+
+ AddPointToHyperplane(h, linExtrapolation, oldBinX, oldBinY);
+
+ if (oldBinY < nBinsY) {
+ AddPointToHyperplane(h, linExtrapolation, oldBinX, oldBinY + 1);
+ }
+
+ // Neighbours right column
+ if (oldBinX < nBinsX) {
+ if (oldBinY >= 2) {
+ AddPointToHyperplane(h, linExtrapolation, oldBinX + 1, oldBinY - 1);
+ }
+
+ AddPointToHyperplane(h, linExtrapolation, oldBinX + 1, oldBinY);
+
+ if (oldBinY < nBinsY) {
+ AddPointToHyperplane(h, linExtrapolation, oldBinX + 1, oldBinY + 1);
+ }
+ }
+
+
+ // Fit 2D-hyperplane
+ if (linExtrapolation->GetNpoints() <= 0)
+ continue;
+
+ if (linExtrapolation->Eval() != 0)// EvalRobust -> Takes much, much, [...], much more time (~hours instead of seconds)
+ continue;
+
+ // Fill the bin of the refined histogram with the extrapolated value
+ Double_t interpolatedValue = linExtrapolation->GetParameter(0) + linExtrapolation->GetParameter(1) * centerX
+ + linExtrapolation->GetParameter(2) * centerY;
+ */
+ Double_t interpolatedValue = h->Interpolate(centerX, centerY) ;//TODO NOW NOW NOW
+ hRefined->SetBinContent(binX, binY, interpolatedValue);
+ }
+ }
+
+ delete linExtrapolation;
+
+ return hRefined;
+}
+
+//______________________________________________________________________________
+void AliPIDResponse::SetTPCEtaMaps(Double_t refineFactorMapX, Double_t refineFactorMapY,
+ Double_t refineFactorSigmaMapX, Double_t refineFactorSigmaMapY)
+{
+ //
+ // Load the TPC eta correction maps from the OADB
+ //
+
+ TString dataType = "DATA";
+ TString period = fLHCperiod.IsNull() ? "No period information" : fLHCperiod;
+
+ if (fIsMC) {
+ dataType = "MC";
+ fRecoPass = 1;
+
+ if (fMCperiodTPC.IsNull()) {
+ AliFatal("MC detected, but no MC period set -> Not changing eta maps!");
+ return;
+ }
+
+ period = fMCperiodTPC;
+ }
+
+ TString defaultObj = Form("Default_%s_pass%d", dataType.Data(), fRecoPass);
+
+ AliInfo(Form("Current period and reco pass: %s.pass%d", period.Data(), fRecoPass));
+ if (fUseTPCEtaCorrection == kFALSE) {
+ // Disable eta correction via setting no maps
+ if (!fTPCResponse.SetEtaCorrMap(0x0))
+ AliInfo("Request to disable TPC eta correction -> Eta correction has been disabled");
+ else
+ AliError("Request to disable TPC eta correction -> Some error occured when unloading the correction maps");
+
+ if (!fTPCResponse.SetSigmaParams(0x0, 0))
+ AliInfo("Request to disable TPC eta correction -> Using old parametrisation for sigma");
+ else
+ AliError("Request to disable TPC eta correction -> Some error occured when unloading the sigma maps");
+
+ return;
+ }
+
+ // Invalidate old maps
+ fTPCResponse.SetEtaCorrMap(0x0);
+ fTPCResponse.SetSigmaParams(0x0, 0);
+
+ // Load the eta correction maps
+ AliOADBContainer etaMapsCont(Form("TPCetaMaps_%s_pass%d", dataType.Data(), fRecoPass));
+
+ Int_t statusCont = etaMapsCont.InitFromFile(Form("%s/COMMON/PID/data/TPCetaMaps.root", fOADBPath.Data()),
+ Form("TPCetaMaps_%s_pass%d", dataType.Data(), fRecoPass));
+ if (statusCont) {
+ AliError("Failed initializing TPC eta correction maps from OADB -> Disabled eta correction");
+ }
+ else {
+ AliInfo(Form("Loading TPC eta correction map from %s/COMMON/PID/data/TPCetaMaps.root", fOADBPath.Data()));
+
+ TH2D* etaMap = 0x0;
+
+ if (fIsMC) {
+ TString searchMap = Form("TPCetaMaps_%s_%s_pass%d", dataType.Data(), period.Data(), fRecoPass);
+ etaMap = dynamic_cast<TH2D *>(etaMapsCont.GetDefaultObject(searchMap.Data()));
+ if (!etaMap) {
+ // Try default object
+ etaMap = dynamic_cast<TH2D *>(etaMapsCont.GetDefaultObject(defaultObj.Data()));
+ }
+ }
+ else {
+ etaMap = dynamic_cast<TH2D *>(etaMapsCont.GetObject(fRun, defaultObj.Data()));
+ }
+
+
+ if (!etaMap) {
+ AliError(Form("TPC eta correction map not found for run %d and also no default map found -> Disabled eta correction!!!", fRun));
+ }
+ else {
+ TH2D* etaMapRefined = RefineHistoViaLinearInterpolation(etaMap, refineFactorMapX, refineFactorMapY);
+
+ if (etaMapRefined) {
+ if (!fTPCResponse.SetEtaCorrMap(etaMapRefined)) {
+ AliError(Form("Failed to set TPC eta correction map for run %d -> Disabled eta correction!!!", fRun));
+ fTPCResponse.SetEtaCorrMap(0x0);
+ }
+ else {
+ AliInfo(Form("Loaded TPC eta correction map (refine factors %.2f/%.2f) from %s/COMMON/PID/data/TPCetaMaps.root: %s",
+ refineFactorMapX, refineFactorMapY, fOADBPath.Data(), fTPCResponse.GetEtaCorrMap()->GetTitle()));
+ }
+
+ delete etaMapRefined;
+ }
+ else {
+ AliError(Form("Failed to set TPC eta correction map for run %d (map was loaded, but couldn't be refined) -> Disabled eta correction!!!", fRun));
+ }
+ }
+ }
+
+ // Load the sigma parametrisation (1/dEdx vs tanTheta_local (~eta))
+ AliOADBContainer etaSigmaMapsCont(Form("TPCetaSigmaMaps_%s_pass%d", dataType.Data(), fRecoPass));
+
+ statusCont = etaSigmaMapsCont.InitFromFile(Form("%s/COMMON/PID/data/TPCetaMaps.root", fOADBPath.Data()),
+ Form("TPCetaSigmaMaps_%s_pass%d", dataType.Data(), fRecoPass));
+ if (statusCont) {
+ AliError("Failed initializing TPC eta sigma maps from OADB -> Using old sigma parametrisation");
+ }
+ else {
+ AliInfo(Form("Loading TPC eta sigma map from %s/COMMON/PID/data/TPCetaMaps.root", fOADBPath.Data()));
+
+ TObjArray* etaSigmaPars = 0x0;
+
+ if (fIsMC) {
+ TString searchMap = Form("TPCetaSigmaMaps_%s_%s_pass%d", dataType.Data(), period.Data(), fRecoPass);
+ etaSigmaPars = dynamic_cast<TObjArray *>(etaSigmaMapsCont.GetDefaultObject(searchMap.Data()));
+ if (!etaSigmaPars) {
+ // Try default object
+ etaSigmaPars = dynamic_cast<TObjArray *>(etaSigmaMapsCont.GetDefaultObject(defaultObj.Data()));
+ }
+ }
+ else {
+ etaSigmaPars = dynamic_cast<TObjArray *>(etaSigmaMapsCont.GetObject(fRun, defaultObj.Data()));
+ }
+
+ if (!etaSigmaPars) {
+ AliError(Form("TPC eta sigma parametrisation not found for run %d -> Using old sigma parametrisation!!!", fRun));
+ }
+ else {
+ TH2D* etaSigmaPar1Map = dynamic_cast<TH2D *>(etaSigmaPars->FindObject("sigmaPar1Map"));
+ TNamed* sigmaPar0Info = dynamic_cast<TNamed *>(etaSigmaPars->FindObject("sigmaPar0"));
+ Double_t sigmaPar0 = 0.0;
+
+ if (sigmaPar0Info) {
+ TString sigmaPar0String = sigmaPar0Info->GetTitle();
+ sigmaPar0 = sigmaPar0String.Atof();
+ }
+ else {
+ // Something is weired because the object for parameter 0 could not be loaded -> New sigma parametrisation can not be used!
+ etaSigmaPar1Map = 0x0;
+ }
+
+ TH2D* etaSigmaPar1MapRefined = RefineHistoViaLinearInterpolation(etaSigmaPar1Map, refineFactorSigmaMapX, refineFactorSigmaMapY);
+
+
+ if (etaSigmaPar1MapRefined) {
+ if (!fTPCResponse.SetSigmaParams(etaSigmaPar1MapRefined, sigmaPar0)) {
+ AliError(Form("Failed to set TPC eta sigma map for run %d -> Using old sigma parametrisation!!!", fRun));
+ fTPCResponse.SetSigmaParams(0x0, 0);
+ }
+ else {
+ AliInfo(Form("Loaded TPC sigma correction map (refine factors %.2f/%.2f) from %s/COMMON/PID/data/TPCetaMaps.root: %s",
+ refineFactorSigmaMapX, refineFactorSigmaMapY, fOADBPath.Data(), fTPCResponse.GetSigmaPar1Map()->GetTitle()));
+ }
+
+ delete etaSigmaPar1MapRefined;
+ }
+ else {
+ AliError(Form("Failed to set TPC eta sigma map for run %d (map was loaded, but couldn't be refined) -> Using old sigma parametrisation!!!",
+ fRun));
+ }
+ }
+ }
+}
+
//______________________________________________________________________________
void AliPIDResponse::SetTPCPidResponseMaster()
{
if(!fTuneMConData) datatype="MC";
fRecoPass=1;
}
-
+
// period
TString period=fLHCperiod;
if (fIsMC && !fTuneMConData) period=fMCperiodTPC;
if(flagT0T0){
t0A= vevent->GetT0TOF()[1];
t0C= vevent->GetT0TOF()[2];
- // t0AC= vevent->GetT0TOF()[0];
- t0AC= t0A/resT0A/resT0A + t0C/resT0C/resT0C;
- resT0AC= TMath::Sqrt(1./resT0A/resT0A + 1./resT0C/resT0C);
- t0AC /= resT0AC*resT0AC;
+ // t0AC= vevent->GetT0TOF()[0];
+ t0AC= t0A/resT0A/resT0A + t0C/resT0C/resT0C;
+ resT0AC= TMath::Sqrt(1./resT0A/resT0A + 1./resT0C/resT0C);
+ t0AC /= resT0AC*resT0AC;
}
Float_t t0t0Best = 0;
if(flagT0T0){
t0A= vevent->GetT0TOF()[1];
t0C= vevent->GetT0TOF()[2];
- // t0AC= vevent->GetT0TOF()[0];
- t0AC= t0A/resT0A/resT0A + t0C/resT0C/resT0C;
- resT0AC= TMath::Sqrt(1./resT0A/resT0A + 1./resT0C/resT0C);
- t0AC /= resT0AC*resT0AC;
+ // t0AC= vevent->GetT0TOF()[0];
+ t0AC= t0A/resT0A/resT0A + t0C/resT0C/resT0C;
+ resT0AC= TMath::Sqrt(1./resT0A/resT0A + 1./resT0C/resT0C);
+ t0AC /= resT0AC*resT0AC;
}
if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
AliVTrack *track=(AliVTrack*)vtrack;
- Double_t mom = track->GetTPCmomentum();
- Double_t sig = track->GetTPCsignal();
- if(fTuneMConData) sig = this->GetTPCsignalTunedOnData(track);
- UInt_t sigN = track->GetTPCsignalN();
-
Double_t nSigma = -999.;
- if (sigN>0) nSigma=fTPCResponse.GetNumberOfSigmas(mom,sig,sigN,type);
+
+ //TODO: TPCsignalTunedOnData not taken into account for the new TPCPIDresponse functions, so take the old functions
+ if (fTuneMConData) {
+ Double_t mom = track->GetTPCmomentum();
+ Double_t sig = this->GetTPCsignalTunedOnData(track);
+ UInt_t sigN = track->GetTPCsignalN();
+
+ if (sigN>0)
+ nSigma = fTPCResponse.GetNumberOfSigmas(mom, sig, sigN, type);
+ }
+ else {
+ nSigma = fTPCResponse.GetNumberOfSigmas(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
+ }
return nSigma;
}
// check quality of the track
if ( (track->GetStatus()&AliVTrack::kTPCin )==0 && (track->GetStatus()&AliVTrack::kTPCout)==0 ) return kDetNoSignal;
- Double_t mom = track->GetTPCmomentum();
-
Double_t dedx=track->GetTPCsignal();
Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
if(fTuneMConData) dedx = this->GetTPCsignalTunedOnData(track);
+ Double_t bethe = 0.;
+ Double_t sigma = 0.;
+
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);
+
+ //TODO: TPCsignalTunedOnData not taken into account for the new TPCPIDresponse functions, so take the old functions
+ if (fTuneMConData) {
+ Double_t mom = track->GetTPCmomentum();
+ bethe=fTPCResponse.GetExpectedSignal(mom,type);
+ sigma=fTPCResponse.GetExpectedSigma(mom,track->GetTPCsignalN(),type);
+ }
+ else {
+ bethe=fTPCResponse.GetExpectedSignal(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
+ sigma=fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
+ }
if (TMath::Abs(dedx-bethe) > fRange*sigma) {
p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
} else {
class TF1;
class TObjArray;
+class TLinearFitter;
class AliVEvent;
class AliTRDPIDResponseObject;
kDetPidOk=1,
kDetMismatch=2
};
-
+
AliITSPIDResponse &GetITSResponse() {return fITSResponse;}
AliTPCPIDResponse &GetTPCResponse() {return fTPCResponse;}
AliTOFPIDResponse &GetTOFResponse() {return fTOFResponse;}
AliTRDPIDResponse &GetTRDResponse() {return fTRDResponse;}
AliEMCALPIDResponse &GetEMCALResponse() {return fEMCALResponse;}
-
+
//buffered PID calculation
Float_t NumberOfSigmas(EDetector detCode, const AliVParticle *track, AliPID::EParticleType type) const;
Float_t NumberOfSigmas(EDetCode detCode, const AliVParticle *track, AliPID::EParticleType type) const;
virtual Float_t NumberOfSigmasITS (const AliVParticle *track, AliPID::EParticleType type) const;
virtual Float_t NumberOfSigmasTPC (const AliVParticle *track, AliPID::EParticleType type) const;
- virtual Float_t NumberOfSigmasTPC (const AliVParticle *track, AliPID::EParticleType type, AliTPCPIDResponse::ETPCdEdxSource dedxSource);
+ virtual Float_t NumberOfSigmasTPC (const AliVParticle *track, AliPID::EParticleType type, AliTPCPIDResponse::ETPCdEdxSource dedxSource) const;
virtual Float_t NumberOfSigmasEMCAL(const AliVParticle *track, AliPID::EParticleType type, Double_t &eop, Double_t showershape[4]) const;
virtual Float_t NumberOfSigmasTOF (const AliVParticle *track, AliPID::EParticleType type) const;
virtual Float_t NumberOfSigmasEMCAL(const AliVParticle *track, AliPID::EParticleType type) const;
// event info
Float_t GetCurrentCentrality() const {return fCurrCentrality;};
+ // TPC setting
+ void SetUseTPCEtaCorrection(Bool_t useEtaCorrection = kTRUE) { fUseTPCEtaCorrection = useEtaCorrection; };
+ Bool_t UseTPCEtaCorrection() const { return fUseTPCEtaCorrection; };
+
// TOF setting
void SetTOFtail(Float_t tail=1.1){if(tail > 0) fTOFtail=tail; else printf("TOF tail should be greater than 0 (nothing done)\n");};
void SetTOFResponse(AliVEvent *vevent,EStartTimeType_t option);
TObjArray *fArrPidResponseMaster; //! TPC pid splines
TF1 *fResolutionCorrection; //! TPC resolution correction
AliOADBContainer* fOADBvoltageMaps; //! container with the voltage maps
-
+ Bool_t fUseTPCEtaCorrection; // Use TPC eta correction
+
AliTRDPIDResponseObject *fTRDPIDResponseObject; //! TRD PID Response Object
Float_t fTOFtail; //! TOF tail effect used in TOF probability
void SetITSParametrisation();
//TPC
+ void SetTPCEtaMaps(Double_t refineFactorMapX = 6.0, Double_t refineFactorMapY = 6.0, Double_t refineFactorSigmaMapX = 6.0,
+ Double_t refineFactorSigmaMapY = 6.0);
void SetTPCPidResponseMaster();
void SetTPCParametrisation();
Double_t GetTPCMultiplicityBin(const AliVEvent * const event);
+
+ // TPC helpers for the eta maps
+ void AddPointToHyperplane(TH2D* h, TLinearFitter* linExtrapolation, Int_t binX, Int_t binY);
+ TH2D* RefineHistoViaLinearInterpolation(TH2D* h, Double_t refineFactorX = 6.0, Double_t refineFactorY = 6.0);
//TRD
void SetTRDPidResponseMaster();
EDetPidStatus GetComputePHOSProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const;
EDetPidStatus GetComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const;
- ClassDef(AliPIDResponse, 10); //PID response handling
+ ClassDef(AliPIDResponse, 11); //PID response handling
};
#endif
// Dariusz Miskowiec, GSI, D.Miskowiec@gsi.de
// ...and some modifications by
// Mikolaj Krzewicki, GSI, mikolaj.krzewicki@cern.ch
+// ...and some modifications plus eta correction functions by
+// Benjamin Hess, University of Tuebingen, bhess@cern.ch
//-----------------------------------------------------------------
#include <TGraph.h>
#include <TSpline.h>
#include <TBits.h>
#include <TMath.h>
+#include <TH2D.h>
#include <AliLog.h>
#include "AliExternalTrackParam.h"
fLowGainOROCthreshold(-40),
fBadOROCthreshhold(-40),
fMaxBadLengthFraction(0.5),
- fCurrentResponseFunction(NULL),
- fCurrentdEdx(0.0),
- fCurrentNPoints(0),
- fCurrentGainScenario(kGainScenarioInvalid),
- fMagField(0.)
+ fMagField(0.),
+ fhEtaCorr(0x0),
+ fhEtaSigmaPar1(0x0),
+ fSigmaPar0(0.0)
{
//
// The default constructor
fLowGainOROCthreshold(-40),
fBadOROCthreshhold(-40),
fMaxBadLengthFraction(0.5),
- fCurrentResponseFunction(NULL),
- fCurrentdEdx(0.0),
- fCurrentNPoints(0),
- fCurrentGainScenario(kGainScenarioInvalid),
- fMagField(0.)
+ fMagField(0.),
+ fhEtaCorr(0x0),
+ fhEtaSigmaPar1(0x0),
+ fSigmaPar0(0.0)
{
//
// The main constructor
for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=param[1];fResN2[i]=param[2];}
}
+
+//_________________________________________________________________________
+AliTPCPIDResponse::~AliTPCPIDResponse()
+{
+ //
+ // Destructor
+ //
+
+ delete fhEtaCorr;
+ fhEtaCorr = 0x0;
+
+ delete fhEtaSigmaPar1;
+ fhEtaSigmaPar1 = 0x0;
+}
+
+
//_________________________________________________________________________
AliTPCPIDResponse::AliTPCPIDResponse(const AliTPCPIDResponse& that):
TNamed(that),
fLowGainOROCthreshold(that.fLowGainOROCthreshold),
fBadOROCthreshhold(that.fBadOROCthreshhold),
fMaxBadLengthFraction(that.fMaxBadLengthFraction),
- fCurrentResponseFunction(NULL),
- fCurrentdEdx(0.0),
- fCurrentNPoints(0),
- fCurrentGainScenario(kGainScenarioInvalid),
- fMagField(that.fMagField)
+ fMagField(that.fMagField),
+ fhEtaCorr(0x0),
+ fhEtaSigmaPar1(0x0),
+ fSigmaPar0(that.fSigmaPar0)
{
//copy ctor
for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=that.fRes0[i];fResN2[i]=that.fResN2[i];}
+
+ // Copy eta maps
+ fhEtaCorr = new TH2D(*(that.fhEtaCorr));
+ fhEtaCorr->SetDirectory(0);
+
+ fhEtaSigmaPar1 = new TH2D(*(that.fhEtaSigmaPar1));
+ fhEtaSigmaPar1->SetDirectory(0);
}
//_________________________________________________________________________
fMaxBadLengthFraction=that.fMaxBadLengthFraction;
fMagField=that.fMagField;
for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=that.fRes0[i];fResN2[i]=that.fResN2[i];}
+
+ delete fhEtaCorr;
+ fhEtaCorr = new TH2D(*(that.fhEtaCorr));
+ fhEtaCorr->SetDirectory(0);
+
+ delete fhEtaSigmaPar1;
+ fhEtaSigmaPar1 = new TH2D(*(that.fhEtaSigmaPar1));
+ fhEtaSigmaPar1->SetDirectory(0);
+
+ fSigmaPar0 = that.fSigmaPar0;
+
return *this;
}
Double_t AliTPCPIDResponse::GetExpectedSignal(const Float_t mom,
AliPID::EParticleType n) const {
//
+ // Deprecated function (for backward compatibility). Please use
+ // GetExpectedSignal(const AliVTrack* track, AliPID::EParticleType species, ETPCdEdxSource dedxSource,
+ // Bool_t correctEta = kTRUE);
+ // instead!
+ //
+ //
// Calculates the expected PID signal as the function of
// the information stored in the track, for the specified particle type
//
// assigned clusters and/or the track dip angle, for example.
//
+ //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...)
+ // !!! Splines for light nuclei need to be normalised to this factor !!!
+ const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(n),2.3);
+
Double_t mass=AliPID::ParticleMassZ(n);
- if (!fUseDatabase) return Bethe(mom/mass);
+ if (!fUseDatabase) return Bethe(mom/mass) * chargeFactor;
//
const TSpline3 * responseFunction = (TSpline3 *) fResponseFunctions.UncheckedAt(n);
- if (!responseFunction) return Bethe(mom/mass);
- //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...)
- // !!! Splines for light nuclei need to be normalised to this factor !!!
- const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(n),2.3);
+ if (!responseFunction) return Bethe(mom/mass) * chargeFactor;
+
return fMIP*responseFunction->Eval(mom/mass)*chargeFactor;
}
//_________________________________________________________________________
Double_t AliTPCPIDResponse::GetExpectedSigma(const Float_t mom,
- const Int_t nPoints,
+ const Int_t nPoints,
AliPID::EParticleType n) const {
//
+ // Deprecated function (for backward compatibility). Please use
+ // GetExpectedSigma(onst AliVTrack* track, AliPID::EParticleType species,
+ // ETPCdEdxSource dedxSource, Bool_t correctEta) instead!
+ //
+ //
// Calculates the expected sigma of the PID signal as the function of
// the information stored in the track, for the specified particle type
//
fResN2[gainScenario]=resN2;
}
+
//_________________________________________________________________________
-Double_t AliTPCPIDResponse::GetExpectedSignal(Double_t mom,
+Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track,
AliPID::EParticleType species,
- const TSpline3* responseFunction) const
+ Double_t /*dEdx*/,
+ const TSpline3* responseFunction,
+ Bool_t correctEta) const
{
// Calculates the expected PID signal as the function of
- // the information stored in the track, for the specified particle type
+ // the information stored in the track and the given parameters,
+ // for the specified particle type
//
// At the moment, these signals are just the results of calling the
- // Bethe-Bloch formula.
+ // Bethe-Bloch formula plus, if desired, taking into account the eta dependence.
// This can be improved. By taking into account the number of
// assigned clusters and/or the track dip angle, for example.
//
-
+ Double_t mom=track->GetTPCmomentum();
Double_t mass=AliPID::ParticleMass(species);
- if (!fUseDatabase) return Bethe(mom/mass);
+ //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...)
+ // !!! Splines for light nuclei need to be normalised to this factor !!!
+ const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(species),2.3);
+
+ if (!responseFunction)
+ return Bethe(mom/mass) * chargeFactor;
- if (!responseFunction) return Bethe(mom/mass);
- return fMIP*responseFunction->Eval(mom/mass);
+ Double_t dEdxSplines = fMIP*responseFunction->Eval(mom/mass) * chargeFactor;
+
+ if (!correctEta)
+ return dEdxSplines;
+
+ //TODO Alternatively take current track dEdx
+ //return dEdxSplines * GetEtaCorrection(track, dEdx);
+ return dEdxSplines * GetEtaCorrection(track, dEdxSplines);
}
+
//_________________________________________________________________________
Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track,
AliPID::EParticleType species,
- ETPCdEdxSource dedxSource)
+ ETPCdEdxSource dedxSource,
+ Bool_t correctEta) const
{
// Calculates the expected PID signal as the function of
// the information stored in the track, for the specified particle type
//
// At the moment, these signals are just the results of calling the
- // Bethe-Bloch formula.
+ // Bethe-Bloch formula plus, if desired, taking into account the eta dependence.
// This can be improved. By taking into account the number of
// assigned clusters and/or the track dip angle, for example.
//
- Double_t mom = track->GetTPCmomentum();
- Double_t mass=AliPID::ParticleMass(species);
+ if (!fUseDatabase) {
+ //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...)
+ // !!! Splines for light nuclei need to be normalised to this factor !!!
+ const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(species),2.3);
- if (!fUseDatabase) return Bethe(mom/mass);
+ return Bethe(track->GetTPCmomentum() / AliPID::ParticleMass(species)) * chargeFactor;
+ }
- const TSpline3* responseFunction = GetResponseFunction(track,species,dedxSource);
- if (!responseFunction) return Bethe(mom/mass);
- return fMIP*responseFunction->Eval(mom/mass);
-
+ Double_t dEdx = -1;
+ Int_t nPoints = -1;
+ ETPCgainScenario gainScenario = kGainScenarioInvalid;
+ TSpline3* responseFunction = 0x0;
+
+ if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) {
+ // Something is wrong with the track -> Return obviously invalid value
+ return -999;
+ }
+
+ // Charge factor already taken into account inside the following function call
+ return GetExpectedSignal(track, species, dEdx, responseFunction, correctEta);
}
//_________________________________________________________________________
//_________________________________________________________________________
TSpline3* AliTPCPIDResponse::GetResponseFunction( const AliVTrack* track,
AliPID::EParticleType species,
- ETPCdEdxSource dedxSource )
+ ETPCdEdxSource dedxSource) const
{
//the splines are stored in an array, different scenarios
- if (ResponseFunctiondEdxN(track,
- species,
- dedxSource))
- return fCurrentResponseFunction;
+ Double_t dEdx = -1;
+ Int_t nPoints = -1;
+ ETPCgainScenario gainScenario = kGainScenarioInvalid;
+ TSpline3* responseFunction = 0x0;
+
+ if (ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
+ return responseFunction;
return NULL;
}
fResponseFunctions.AddAtAndExpand(o,ResponseFunctionIndex(species,gainScenario));
}
+
//_________________________________________________________________________
-Double_t AliTPCPIDResponse::GetExpectedSigma( Double_t mom,
- Int_t nPoints,
- AliPID::EParticleType species,
- ETPCgainScenario gainScenario,
- const TSpline3* responseFunction) const
+Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track,
+ AliPID::EParticleType species,
+ ETPCgainScenario gainScenario,
+ Double_t dEdx,
+ Int_t nPoints,
+ const TSpline3* responseFunction,
+ Bool_t correctEta) const
{
// Calculates the expected sigma of the PID signal as the function of
- // the information stored in the track, for the specified particle type
- // and dedx scenario
+ // the information stored in the track and the given parameters,
+ // for the specified particle type
//
- if (!responseFunction) return 999;
- if (nPoints != 0)
- return GetExpectedSignal(mom,species,responseFunction) *
- fRes0[gainScenario] * sqrt(1. + fResN2[gainScenario]/nPoints);
- else
- return GetExpectedSignal(mom,species,responseFunction)*fRes0[gainScenario];
+ if (!responseFunction)
+ return 999;
+
+ // If no sigma map is available or if no eta correction is requested (sigma maps only for corrected eta!), use the old parametrisation
+ if (!fhEtaSigmaPar1 || !correctEta) {
+ if (nPoints != 0)
+ return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE) *
+ fRes0[gainScenario] * sqrt(1. + fResN2[gainScenario]/nPoints);
+ else
+ return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE)*fRes0[gainScenario];
+ }
+
+ if (nPoints > 0) {
+ Double_t sigmaPar1 = GetSigmaPar1(track, species, dEdx, responseFunction);
+
+ return GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE)*sqrt(fSigmaPar0 * fSigmaPar0 + sigmaPar1 * sigmaPar1 / nPoints);
+ }
+ else {
+ // One should never have/take tracks with 0 dEdx clusters!
+ return 999;
+ }
}
+
//_________________________________________________________________________
Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track,
AliPID::EParticleType species,
- ETPCdEdxSource dedxSource)
+ ETPCdEdxSource dedxSource,
+ Bool_t correctEta) const
{
// Calculates the expected sigma of the PID signal as the function of
// the information stored in the track, for the specified particle type
// and dedx scenario
//
- if (!ResponseFunctiondEdxN(track,
- species,
- dedxSource)) return 998; //TODO: better handling!
-
- return GetExpectedSigma( track->GetTPCmomentum(),
- fCurrentNPoints,
- species,
- fCurrentGainScenario,
- fCurrentResponseFunction );
+ Double_t dEdx = -1;
+ Int_t nPoints = -1;
+ ETPCgainScenario gainScenario = kGainScenarioInvalid;
+ TSpline3* responseFunction = 0x0;
+
+ if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
+ return 999; //TODO: better handling!
+
+ return GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta);
}
+
//_________________________________________________________________________
Float_t AliTPCPIDResponse::GetNumberOfSigmas(const AliVTrack* track,
AliPID::EParticleType species,
- ETPCdEdxSource dedxSource)
+ ETPCdEdxSource dedxSource,
+ Bool_t correctEta) const
{
- //calculates the number of sigmas of the PID signal from the expected value
- //for a gicven particle species in the presence of multiple gain scenarios
+ //Calculates the number of sigmas of the PID signal from the expected value
+ //for a given particle species in the presence of multiple gain scenarios
//inside the TPC
-
- if (!ResponseFunctiondEdxN(track,
- species,
- dedxSource)) return 997; //TODO: better handling!
-
- Double_t mom = track->GetTPCmomentum();
- Double_t bethe=GetExpectedSignal(mom,species,fCurrentResponseFunction);
- Double_t sigma=GetExpectedSigma( mom,
- fCurrentNPoints,
- species,
- fCurrentGainScenario,
- fCurrentResponseFunction );
- return (fCurrentdEdx-bethe)/sigma;
+
+ Double_t dEdx = -1;
+ Int_t nPoints = -1;
+ ETPCgainScenario gainScenario = kGainScenarioInvalid;
+ TSpline3* responseFunction = 0x0;
+
+ if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
+ return -999; //TODO: Better handling!
+
+ Double_t bethe = GetExpectedSignal(track, species, dEdx, responseFunction, correctEta);
+ Double_t sigma = GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta);
+ // 999 will be returned by GetExpectedSigma e.g. in case of 0 dEdx clusters
+ if (sigma >= 998)
+ return -999;
+ else
+ return (dEdx-bethe)/sigma;
}
//_________________________________________________________________________
Bool_t AliTPCPIDResponse::ResponseFunctiondEdxN( const AliVTrack* track,
AliPID::EParticleType species,
- ETPCdEdxSource dedxSource )
+ ETPCdEdxSource dedxSource,
+ Double_t& dEdx, Int_t& nPoints, ETPCgainScenario& gainScenario, TSpline3** responseFunction) const
{
// Calculates the right parameters for PID
// dEdx parametrization for the proper gain scenario, dEdx
// and preferred source of dedx.
// returns true on success
+
+ if (dedxSource == kdEdxDefault) {
+ // Fast handling for default case. In addition: Keep it simple (don't call additional functions) to
+ // avoid possible bugs
+ dEdx = track->GetTPCsignal();
+ nPoints = track->GetTPCsignalN();
+ gainScenario = kDefault;
+
+ TObject* obj = fResponseFunctions.UncheckedAt(ResponseFunctionIndex(species,gainScenario));
+ *responseFunction = dynamic_cast<TSpline3*>(obj); //TODO:maybe static cast?
+
+ return kTRUE;
+ }
+
+
Double32_t signal[4]; //0: IROC, 1: OROC medium, 2:OROC long, 3: OROC all (def. truncation used)
Char_t ncl[3]; //same
Char_t nrows[3]; //same
- AliTPCdEdxInfo* dEdxInfo = track->GetTPCdEdxInfo();
+ const AliTPCdEdxInfo* dEdxInfo = track->GetTPCdEdxInfo();
if (!dEdxInfo && dedxSource!=kdEdxDefault) //in one case its ok if we dont have the info
{
AliError("AliTPCdEdxInfo not available");
- InvalidateCurrentValues();
return kFALSE;
}
EChamberStatus trackStatus = TrackStatus(track,2);
if (trackStatus==kChamberOff || trackStatus==kChamberLowGain)
{
- InvalidateCurrentValues();
return kFALSE;
}
switch (dedxSource)
{
- case kdEdxDefault:
- {
- fCurrentdEdx = track->GetTPCsignal();
- fCurrentNPoints = track->GetTPCsignalN();
- fCurrentGainScenario = kDefault;
- break;
- }
case kdEdxOROC:
{
- fCurrentdEdx = signal[3];
- fCurrentNPoints = ncl[2]+ncl[1];
- fCurrentGainScenario = kOROChigh;
+ dEdx = signal[3];
+ nPoints = ncl[2]+ncl[1];
+ gainScenario = kOROChigh;
break;
}
case kdEdxHybrid:
EChamberStatus status = TrackStatus(track,1);
if (status!=kChamberHighGain)
{
- fCurrentdEdx = signal[3];
- fCurrentNPoints = ncl[2]+ncl[1];
- fCurrentGainScenario = kOROChigh;
+ dEdx = signal[3];
+ nPoints = ncl[2]+ncl[1];
+ gainScenario = kOROChigh;
}
else
{
- fCurrentdEdx = track->GetTPCsignal();
- fCurrentNPoints = track->GetTPCsignalN();
- fCurrentGainScenario = kALLhigh;
+ dEdx = track->GetTPCsignal();
+ nPoints = track->GetTPCsignalN();
+ gainScenario = kALLhigh;
}
break;
}
default:
{
- fCurrentdEdx = 0.;
- fCurrentNPoints = 0;
- fCurrentGainScenario = kGainScenarioInvalid;
+ dEdx = 0.;
+ nPoints = 0;
+ gainScenario = kGainScenarioInvalid;
return kFALSE;
}
}
- TObject* obj = fResponseFunctions.UncheckedAt(ResponseFunctionIndex(species,fCurrentGainScenario));
- fCurrentResponseFunction = dynamic_cast<TSpline3*>(obj); //TODO:maybe static cast?
+ TObject* obj = fResponseFunctions.UncheckedAt(ResponseFunctionIndex(species,gainScenario));
+ *responseFunction = dynamic_cast<TSpline3*>(obj); //TODO:maybe static cast?
+
return kTRUE;
}
+
+//_________________________________________________________________________
+Double_t AliTPCPIDResponse::GetEtaCorrection(const AliVTrack *track, Double_t dEdxSplines) const
+{
+ //
+ // Get eta correction for the given parameters.
+ //
+
+ if (!fhEtaCorr)
+ return 1.;
+
+ Double_t tpcSignal = dEdxSplines;
+
+ if (tpcSignal < 1.)
+ return 1.;
+
+ // For ESD tracks, the local tanTheta could be used (esdTrack->GetInnerParam()->GetTgl()).
+ // However, this value is not available for AODs and, thus, not for AliVTrack.
+ // Fortunately, the following formula allows to approximate the local tanTheta with the
+ // global theta angle -> This is for by far most of the tracks the same, but gives at
+ // maybe the percent level differences within +- 0.2 in tanTheta -> Which is still ok.
+ Double_t tanTheta = TMath::Tan(-track->Theta() + TMath::Pi() / 2.0);
+ Int_t binX = fhEtaCorr->GetXaxis()->FindBin(tanTheta);
+ Int_t binY = fhEtaCorr->GetYaxis()->FindBin(1. / tpcSignal);
+
+ if (binX == 0)
+ binX = 1;
+ if (binX > fhEtaCorr->GetXaxis()->GetNbins())
+ binX = fhEtaCorr->GetXaxis()->GetNbins();
+
+ if (binY == 0)
+ binY = 1;
+ if (binY > fhEtaCorr->GetYaxis()->GetNbins())
+ binY = fhEtaCorr->GetYaxis()->GetNbins();
+
+ return fhEtaCorr->GetBinContent(binX, binY);
+}
+
+
+//_________________________________________________________________________
+Double_t AliTPCPIDResponse::GetEtaCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const
+{
+ //
+ // Get eta correction for the given track.
+ //
+
+ if (!fhEtaCorr)
+ return 1.;
+
+ Double_t dEdx = -1;
+ Int_t nPoints = -1;
+ ETPCgainScenario gainScenario = kGainScenarioInvalid;
+ TSpline3* responseFunction = 0x0;
+
+ if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
+ return 1.;
+
+ Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE);
+
+ //TODO Alternatively take current track dEdx
+ //return GetEtaCorrection(track, dEdx);
+
+ return GetEtaCorrection(track, dEdxSplines);
+}
+
+
+//_________________________________________________________________________
+Double_t AliTPCPIDResponse::GetEtaCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const
+{
+ //
+ // Get eta corrected dEdx for the given track. For the correction, the expected dEdx of
+ // the specified species will be used. If the species is set to AliPID::kUnknown, the
+ // dEdx of the track is used instead.
+ // WARNING: In the latter case, the eta correction might not be as good as if the
+ // expected dEdx is used, which is the way the correction factor is designed
+ // for.
+ // In any case, one has to decide carefully to which expected signal one wants to
+ // compare the corrected value - to the corrected or uncorrected.
+ // Anyhow, a safer way of looking e.g. at the n-sigma is to call
+ // the corresponding function GetNumberOfSigmas!
+ //
+
+ Double_t dEdx = -1;
+ Int_t nPoints = -1;
+ ETPCgainScenario gainScenario = kGainScenarioInvalid;
+ TSpline3* responseFunction = 0x0;
+
+ // Note: In case of species == AliPID::kUnknown, the responseFunction might not be set. However, in this case
+ // it is not used anyway, so this causes no trouble.
+ if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
+ return -1.;
+
+ Double_t etaCorr = 0.;
+
+ if (species < AliPID::kUnknown) {
+ Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE);
+ etaCorr = GetEtaCorrection(track, dEdxSplines);
+ }
+ else {
+ etaCorr = GetEtaCorrection(track, dEdx);
+ }
+
+ if (etaCorr <= 0)
+ return -1.;
+
+ return dEdx / etaCorr;
+}
+
+
+
+//_________________________________________________________________________
+Double_t AliTPCPIDResponse::GetSigmaPar1(const AliVTrack *track, AliPID::EParticleType species, Double_t dEdx, const TSpline3* responseFunction) const
+{
+ //
+ // Get parameter 1 of sigma parametrisation of TPC dEdx from the histogram for the given track.
+ //
+
+ if (!fhEtaSigmaPar1)
+ return 999;
+
+ // The sigma maps are created with data sets that are already eta corrected and for which the
+ // splines have been re-created. Therefore, the value for the lookup needs to be the value of
+ // the splines without any additional eta correction.
+ // NOTE: This is due to the method the maps are created. The track dEdx (not the expected one!)
+ // is corrected to uniquely related a momemtum bin with an expected dEdx, where the expected dEdx
+ // equals the track dEdx for all eta (thanks to the correction and the re-fit of the splines).
+ // Consequently, looking up the uncorrected expected dEdx at a given tanTheta yields the correct
+ // sigma parameter!
+ // Also: It has to be the spline dEdx, since one wants to get the sigma for the assumption(!)
+ // of such a particle, which by assumption then has this dEdx value
+
+ Double_t dEdxExpected = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE);
+
+ if (dEdxExpected < 1.)
+ return 999;
+
+ // For ESD tracks, the local tanTheta could be used (esdTrack->GetInnerParam()->GetTgl()).
+ // However, this value is not available for AODs and, thus, not or AliVTrack.
+ // Fortunately, the following formula allows to approximate the local tanTheta with the
+ // global theta angle -> This is for by far most of the tracks the same, but gives at
+ // maybe the percent level differences within +- 0.2 in tanTheta -> Which is still ok.
+ Double_t tanTheta = TMath::Tan(-track->Theta() + TMath::Pi() / 2.0);
+ Int_t binX = fhEtaSigmaPar1->GetXaxis()->FindBin(tanTheta);
+ Int_t binY = fhEtaSigmaPar1->GetYaxis()->FindBin(1. / dEdxExpected);
+
+ if (binX == 0)
+ binX = 1;
+ if (binX > fhEtaSigmaPar1->GetXaxis()->GetNbins())
+ binX = fhEtaSigmaPar1->GetXaxis()->GetNbins();
+
+ if (binY == 0)
+ binY = 1;
+ if (binY > fhEtaSigmaPar1->GetYaxis()->GetNbins())
+ binY = fhEtaSigmaPar1->GetYaxis()->GetNbins();
+
+ return fhEtaSigmaPar1->GetBinContent(binX, binY);
+}
+
+
+//_________________________________________________________________________
+Double_t AliTPCPIDResponse::GetSigmaPar1(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const
+{
+ //
+ // Get eta correction for the given track.
+ //
+
+ if (!fhEtaSigmaPar1)
+ return 999;
+
+ Double_t dEdx = -1;
+ Int_t nPoints = -1;
+ ETPCgainScenario gainScenario = kGainScenarioInvalid;
+ TSpline3* responseFunction = 0x0;
+
+ if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
+ return 999;
+
+ return GetSigmaPar1(track, species, dEdx, responseFunction);
+}
+
+
+//_________________________________________________________________________
+Bool_t AliTPCPIDResponse::SetEtaCorrMap(TH2D* hMap)
+{
+ //
+ // Load map for TPC eta correction (a copy is stored and will be deleted automatically).
+ // If hMap is 0x0,the eta correction will be disabled and kFALSE is returned.
+ // If the map can be set, kTRUE is returned.
+ //
+
+ delete fhEtaCorr;
+
+ if (!hMap) {
+ fhEtaCorr = 0x0;
+
+ return kFALSE;
+ }
+
+ fhEtaCorr = (TH2D*)(hMap->Clone());
+ fhEtaCorr->SetDirectory(0);
+
+ return kTRUE;
+}
+
+//_________________________________________________________________________
+Bool_t AliTPCPIDResponse::SetSigmaParams(TH2D* hSigmaPar1Map, Double_t sigmaPar0)
+{
+ //
+ // Load map for TPC sigma map (a copy is stored and will be deleted automatically):
+ // Parameter 1 is stored as a 2D map (1/dEdx vs. tanTheta_local) and parameter 0 is
+ // a constant. If hSigmaPar1Map is 0x0, the old sigma parametrisation will be used
+ // (and sigmaPar0 is ignored!) and kFALSE is returned.
+ // If the map can be set, sigmaPar0 is also set and kTRUE will be returned.
+ //
+
+ delete fhEtaSigmaPar1;
+
+ if (!hSigmaPar1Map) {
+ fhEtaSigmaPar1 = 0x0;
+ fSigmaPar0 = 0.0;
+
+ return kFALSE;
+ }
+
+ fhEtaSigmaPar1 = (TH2D*)(hSigmaPar1Map->Clone());
+ fhEtaSigmaPar1->SetDirectory(0);
+ fSigmaPar0 = sigmaPar0;
+
+ return kTRUE;
+}
+
+
//_________________________________________________________________________
Bool_t AliTPCPIDResponse::sectorNumbersInOut(const AliVTrack* track,
Double_t innerRadius,
return 0.0;
}
-//_____________________________________________________________________________
-void AliTPCPIDResponse::InvalidateCurrentValues()
-{
- //make the "current" stored values from the last processed track invalid
- fCurrentResponseFunction=NULL;
- fCurrentdEdx=0.;
- fCurrentNPoints=0;
- fCurrentGainScenario=kGainScenarioInvalid;
-}
-
//_____________________________________________________________________________
Bool_t AliTPCPIDResponse::TrackApex(const AliVTrack* track, Float_t magField, Double_t position[3]) const
{
position[1]=b[1]+b[1]*TMath::Abs(r)/norm;
position[2]=0.;
return kTRUE;
-}
-
+}
\ No newline at end of file
// With many additions and modifications suggested by
// Alexander Kalweit, GSI, alexander.philipp.kalweit@cern.ch
// Dariusz Miskowiec, GSI, D.Miskowiec@gsi.de
+// ...and some modifications by
+// Mikolaj Krzewicki, GSI, mikolaj.krzewicki@cern.ch
+// ...and some modifications plus eta correction functions by
+// Benjamin Hess, University of Tuebingen, bhess@cern.ch
//-------------------------------------------------------
#include <Rtypes.h>
#include <TObjArray.h>
#include "AliPID.h"
+#include "AliVTrack.h"
-class AliVTrack;
+class TH2D;
class TSpline3;
class AliTPCPIDResponse: public TNamed {
AliTPCPIDResponse(const Double_t *param);
AliTPCPIDResponse(const AliTPCPIDResponse&);
AliTPCPIDResponse& operator=(const AliTPCPIDResponse&);
- virtual ~AliTPCPIDResponse() {}
+ virtual ~AliTPCPIDResponse();
enum EChamberStatus {
kChamberOff=0,
void SetMagField(Double_t mf) { fMagField=mf; }
+ const TH2D* GetEtaCorrMap() const { return fhEtaCorr; };
+ Bool_t SetEtaCorrMap(TH2D* hMap);
+
+ Double_t GetEtaCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource = kdEdxDefault) const;
+
+ Double_t GetEtaCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource = kdEdxDefault) const;
+
+ const TH2D* GetSigmaPar1Map() const { return fhEtaSigmaPar1; };
+ Double_t GetSigmaPar0() const { return fSigmaPar0; };
+ Bool_t SetSigmaParams(TH2D* hSigmaPar1Map, Double_t sigmaPar0);
+
+ Double_t GetSigmaPar1(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource = kdEdxDefault) const;
+
//NEW
void SetSigma(Float_t res0, Float_t resN2, ETPCgainScenario gainScenario );
- Double_t GetExpectedSignal( Double_t momentum,
- AliPID::EParticleType species,
- const TSpline3* responseFunction ) const;
Double_t GetExpectedSignal( const AliVTrack* track,
AliPID::EParticleType species,
- ETPCdEdxSource dedxSource );
+ ETPCdEdxSource dedxSource = kdEdxDefault,
+ Bool_t correctEta = kFALSE) const;//TODO: In future, default kTRUE
Double_t GetExpectedSigma( const AliVTrack* track,
AliPID::EParticleType species,
- ETPCdEdxSource dedxSource );
- Double_t GetExpectedSigma( Double_t mom,
- Int_t nPoints,
- AliPID::EParticleType species,
- ETPCgainScenario gainScenario,
- const TSpline3* responseFunction) const;
+ ETPCdEdxSource dedxSource = kdEdxDefault,
+ Bool_t correctEta = kFALSE) const;//TODO: In future, default kTRUE
Float_t GetNumberOfSigmas( const AliVTrack* track,
AliPID::EParticleType species,
- ETPCdEdxSource dedxSource );
+ ETPCdEdxSource dedxSource = kdEdxDefault,
+ Bool_t correctEta = kFALSE) const;//TODO: In future, default kTRUE
void SetResponseFunction(TObject* o,
AliPID::EParticleType type,
ETPCgainScenario gainScenario ) const;
TSpline3* GetResponseFunction( const AliVTrack* track,
AliPID::EParticleType species,
- ETPCdEdxSource dedxSource );
+ ETPCdEdxSource dedxSource = kdEdxDefault) const;
Bool_t ResponseFunctiondEdxN(const AliVTrack* track,
AliPID::EParticleType species,
- ETPCdEdxSource dedxSource);
+ ETPCdEdxSource dedxSource,
+ Double_t& dEdx, Int_t& nPoints, ETPCgainScenario& gainScenario, TSpline3** responseFunction) const;
Bool_t sectorNumbersInOut(const AliVTrack* track,
Double_t innerRadius, Double_t outerRadius,
Float_t& phiIn, Float_t& phiOut,
ETPCgainScenario gainScenario ) const;
void ResetSplines();
- void InvalidateCurrentValues();
- TSpline3* GetCurrentResponseFunction() const {return fCurrentResponseFunction;}
- Double_t GetCurrentdEdx() const {return fCurrentdEdx;}
- Int_t GetCurrentNPoints() const {return fCurrentNPoints;}
- ETPCgainScenario GetCurrentGainScenario() const {return fCurrentGainScenario;}
-
//OLD
Double_t GetExpectedSignal(const Float_t mom,
AliPID::EParticleType n=AliPID::kKaon) const;
const Float_t dEdx,
const Int_t nPoints,
AliPID::EParticleType n=AliPID::kKaon) const {
-
+ //
+ // Deprecated function (for backward compatibility). Please use
+ // GetNumberOfSigmas(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource )
+ // instead!TODO
+ //
+
Double_t bethe=GetExpectedSignal(mom,n);
Double_t sigma=GetExpectedSigma(mom,nPoints,n);
return (dEdx-bethe)/sigma;
Float_t GetRes0(ETPCgainScenario s) const { return fRes0[s]; }
Float_t GetResN2(ETPCgainScenario s) const { return fResN2[s]; }
+protected:
+ Double_t GetExpectedSignal(const AliVTrack* track,
+ AliPID::EParticleType species,
+ Double_t dEdx,
+ const TSpline3* responseFunction,
+ Bool_t correctEta) const;
+
+ Double_t GetExpectedSigma(const AliVTrack* track,
+ AliPID::EParticleType species,
+ ETPCgainScenario gainScenario,
+ Double_t dEdx,
+ Int_t nPoints,
+ const TSpline3* responseFunction,
+ Bool_t correctEta) const;
+
+ Double_t GetEtaCorrection(const AliVTrack *track, Double_t dEdxSplines) const;
+
+ Double_t GetSigmaPar1(const AliVTrack *track, AliPID::EParticleType species,
+ Double_t dEdx, const TSpline3* responseFunction) const;
+
private:
Float_t fMIP; // dEdx for MIP
Float_t fRes0[fgkNumberOfGainScenarios]; // relative dEdx resolution rel sigma = fRes0*sqrt(1+fResN2/npoint)
Float_t fBadOROCthreshhold; //voltage threshold for bad OROCS
Float_t fMaxBadLengthFraction; //the maximum allowed fraction of track length in a bad sector.
- TSpline3* fCurrentResponseFunction; //!response function for current track
- Double_t fCurrentdEdx; //!dEdx for currently processed track
- Int_t fCurrentNPoints; //!number of points used for dEdx calculation for current track
- ETPCgainScenario fCurrentGainScenario; //!gain scenario used for current track
Int_t sectorNumber(Double_t phi) const;
Double_t fMagField; //! Magnetic field
static const char* fgkGainScenarioName[fgkNumberOfGainScenarios+1];
- ClassDef(AliTPCPIDResponse,4) // TPC PID class
+ TH2D* fhEtaCorr; //! Map for TPC eta correction
+ TH2D* fhEtaSigmaPar1; //! Map for parameter 1 of the dEdx sigma parametrisation
+
+ Double_t fSigmaPar0; // Parameter 0 of the dEdx sigma parametrisation
+
+ ClassDef(AliTPCPIDResponse,5) // TPC PID class
};
#endif
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff