TPC PID update:
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 12 Nov 2012 12:23:58 +0000 (12:23 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 12 Nov 2012 12:23:58 +0000 (12:23 +0000)
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)

ANALYSIS/AliAnalysisTaskPIDResponse.cxx
ANALYSIS/AliAnalysisTaskPIDResponse.h
ANALYSIS/macros/AddTaskPIDResponse.C
OADB/COMMON/PID/data/TPCPIDResponse.root
OADB/COMMON/PID/data/TPCetaMaps.root [new file with mode: 0644]
STEER/STEERBase/AliPIDResponse.cxx
STEER/STEERBase/AliPIDResponse.h
STEER/STEERBase/AliTPCPIDResponse.cxx
STEER/STEERBase/AliTPCPIDResponse.h

index ac57341e0dfb06279a837de7897cf0542996d110..c112754b11a944b569cb0e9971c0c6e4b2f48228 100644 (file)
@@ -43,7 +43,8 @@ fRun(0),
 fOldRun(0),
 fRecoPass(0),
 fIsTunedOnData(kFALSE),
-fRecoPassTuned(0)  
+fRecoPassTuned(0),
+fUseTPCEtaCorrection(kFALSE)//TODO: In future, default kTRUE  
 {
   //
   // Dummy constructor
@@ -62,7 +63,8 @@ fRun(0),
 fOldRun(0),
 fRecoPass(0),
 fIsTunedOnData(kFALSE),
-fRecoPassTuned(0)
+fRecoPassTuned(0),
+fUseTPCEtaCorrection(kFALSE)//TODO: In future, default kTRUE
 {
   //
   // Default constructor
@@ -131,6 +133,7 @@ void AliAnalysisTaskPIDResponse::UserExec(Option_t */*option*/)
     fOldRun=fRun;
   }
 
+  fPIDResponse->SetUseTPCEtaCorrection(fUseTPCEtaCorrection);
   fPIDResponse->InitialiseEvent(event,fRecoPass);
   AliESDpid *pidresp = dynamic_cast<AliESDpid*>(fPIDResponse);
   if(pidresp && AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()){
index 59805b1a4742ad10ad1eeb53780f534aeb103bfb..ce2a0b38ed34f2e762bc76eb41a0074aa6414a71 100644 (file)
@@ -43,6 +43,9 @@ public:
   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; }
 
@@ -60,12 +63,14 @@ private:
   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
index 1fcaef1ed7a5915b5b981fd8c05e832d985e23cd..fdc40cd12d51ae07ad063e7832e428a2f2647038 100644 (file)
@@ -1,6 +1,7 @@
 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();
@@ -45,6 +46,7 @@ AliAnalysisTask *AddTaskPIDResponse(Bool_t isMC=kFALSE, Bool_t autoMCesd=kTRUE,
   if(isMC&&tuneOnData) pidTask->SetTuneOnData(kTRUE,recoPass);
   pidTask->SetCachePID(cachePID);
   pidTask->SetSpecialDetectorResponse(detResponse);
+  pidTask->SetUseTPCEtaCorrection(useTPCEtaCorrection);
   mgr->AddTask(pidTask);
   
 //   AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("PIDResponseQA",
index 11dfe07e976e50a29470b578fff84a365913b972..7a642ce0dedebcf5e8bbd8b1bafc0d43141198c0 100644 (file)
Binary files a/OADB/COMMON/PID/data/TPCPIDResponse.root and b/OADB/COMMON/PID/data/TPCPIDResponse.root differ
diff --git a/OADB/COMMON/PID/data/TPCetaMaps.root b/OADB/COMMON/PID/data/TPCetaMaps.root
new file mode 100644 (file)
index 0000000..9187100
Binary files /dev/null and b/OADB/COMMON/PID/data/TPCetaMaps.root differ
index d89780b0a3745db90d62ce9d2db62ad115a2b492..3dcfb1ce649fda36bc0bfcaca832f10c4a2f41da 100644 (file)
 #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>
@@ -75,6 +77,7 @@ fResT0AC(55.),
 fArrPidResponseMaster(NULL),
 fResolutionCorrection(NULL),
 fOADBvoltageMaps(NULL),
+fUseTPCEtaCorrection(kFALSE),//TODO: In future, default kTRUE
 fTRDPIDResponseObject(NULL),
 fTOFtail(1.1),
 fTOFPIDParams(NULL),
@@ -132,6 +135,7 @@ fResT0AC(55.),
 fArrPidResponseMaster(NULL),
 fResolutionCorrection(NULL),
 fOADBvoltageMaps(NULL),
+fUseTPCEtaCorrection(other.fUseTPCEtaCorrection),
 fTRDPIDResponseObject(NULL),
 fTOFtail(1.1),
 fTOFPIDParams(NULL),
@@ -180,6 +184,7 @@ AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other)
     fArrPidResponseMaster=NULL;
     fResolutionCorrection=NULL;
     fOADBvoltageMaps=NULL;
+       fUseTPCEtaCorrection=other.fUseTPCEtaCorrection;
     fTRDPIDResponseObject=NULL;
     fEMCALPIDParams=NULL;
     fTOFtail=1.1;
@@ -271,12 +276,12 @@ Float_t AliPIDResponse::NumberOfSigmasTPC(const AliVParticle *vtrack, AliPID::EP
 //______________________________________________________________________________
 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;
 }
@@ -400,7 +405,7 @@ Float_t  AliPIDResponse::NumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID:
 
 
 //______________________________________________________________________________
-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'
@@ -622,6 +627,7 @@ void AliPIDResponse::ExecNewRun()
   
   SetTPCPidResponseMaster();
   SetTPCParametrisation();
+  SetTPCEtaMaps();
 
   SetTRDPidResponseMaster(); 
   InitializeTRDResponse();
@@ -725,6 +731,292 @@ void AliPIDResponse::SetITSParametrisation()
   //
 }
 
+//______________________________________________________________________________
+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()
 {
@@ -799,7 +1091,7 @@ void AliPIDResponse::SetTPCParametrisation()
       if(!fTuneMConData) datatype="MC";
       fRecoPass=1;
   }
-  
+
   // period
   TString period=fLHCperiod;
   if (fIsMC && !fTuneMConData) period=fMCperiodTPC;
@@ -1261,10 +1553,10 @@ void AliPIDResponse::SetTOFResponse(AliVEvent *vevent,EStartTimeType_t option){
        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;
@@ -1336,10 +1628,10 @@ void AliPIDResponse::SetTOFResponse(AliVEvent *vevent,EStartTimeType_t option){
        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){
@@ -1442,13 +1734,20 @@ Float_t AliPIDResponse::GetNumberOfSigmasTPC(const AliVParticle *vtrack, AliPID:
   
   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;
 }
@@ -1600,17 +1899,27 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTPCProbability  (const A
   // 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 {
index 55441e6a7bb2bcdadcdbb3c048c2dd5c6034bcc7..2cdf5f83a14087dc7323d376c31a45757b97d852 100644 (file)
@@ -25,6 +25,7 @@
 
 class TF1;
 class TObjArray;
+class TLinearFitter;
 
 class AliVEvent;
 class AliTRDPIDResponseObject;
@@ -66,20 +67,20 @@ public:
     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;
@@ -123,6 +124,10 @@ public:
   // 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);
@@ -172,7 +177,8 @@ private:
   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
@@ -196,9 +202,15 @@ private:
   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();
@@ -231,7 +243,7 @@ private:
   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
index 9a9def6a4b88ea2be032f7e0c5aaaee296c24b6f..4ec93ca1a671092af5a9699673b72558f4b78026 100644 (file)
@@ -22,6 +22,8 @@
 //      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>
@@ -29,6 +31,7 @@
 #include <TSpline.h>
 #include <TBits.h>
 #include <TMath.h>
+#include <TH2D.h>
 
 #include <AliLog.h>
 #include "AliExternalTrackParam.h"
@@ -65,11 +68,10 @@ AliTPCPIDResponse::AliTPCPIDResponse():
   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
@@ -96,11 +98,10 @@ AliTPCPIDResponse::AliTPCPIDResponse(const Double_t *param):
   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
@@ -108,6 +109,22 @@ AliTPCPIDResponse::AliTPCPIDResponse(const Double_t *param):
   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),
@@ -127,14 +144,20 @@ AliTPCPIDResponse::AliTPCPIDResponse(const AliTPCPIDResponse& 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);
 }
 
 //_________________________________________________________________________
@@ -159,6 +182,17 @@ AliTPCPIDResponse& AliTPCPIDResponse::operator=(const AliTPCPIDResponse& that)
   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;
 }
 
@@ -210,6 +244,12 @@ void AliTPCPIDResponse::SetSigma(Float_t res0, Float_t resN2) {
 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 
   //  
@@ -219,24 +259,31 @@ Double_t AliTPCPIDResponse::GetExpectedSignal(const Float_t mom,
   // 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 
   //  
@@ -259,52 +306,80 @@ void AliTPCPIDResponse::SetSigma(Float_t res0, Float_t resN2, ETPCgainScenario g
   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);
 }
   
 //_________________________________________________________________________
@@ -318,14 +393,17 @@ TSpline3* AliTPCPIDResponse::GetResponseFunction( AliPID::EParticleType type,
 //_________________________________________________________________________
 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;
 }
@@ -355,74 +433,100 @@ void AliTPCPIDResponse::SetResponseFunction( TObject* o,
   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 
@@ -431,15 +535,29 @@ Bool_t AliTPCPIDResponse::ResponseFunctiondEdxN( const AliVTrack* track,
   // 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;
   }
 
@@ -449,24 +567,16 @@ Bool_t AliTPCPIDResponse::ResponseFunctiondEdxN( const AliVTrack* track,
   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:
@@ -475,31 +585,264 @@ Bool_t AliTPCPIDResponse::ResponseFunctiondEdxN( const AliVTrack* track,
         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,
@@ -707,16 +1050,6 @@ Float_t AliTPCPIDResponse::MaxClusterRadius(const AliVTrack* track) const
   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
 {
@@ -750,5 +1083,4 @@ Bool_t AliTPCPIDResponse::TrackApex(const AliVTrack* track, Float_t magField, Do
   position[1]=b[1]+b[1]*TMath::Abs(r)/norm;
   position[2]=0.;
   return kTRUE;
-}
-
+}
\ No newline at end of file
index eb6a3e383f6ec710c3647f1546b48f3b191495b9..ce09635572d7fc65282bafc2e0b0ac69369885f8 100644 (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>
 
@@ -18,8 +22,9 @@
 #include <TObjArray.h>
 
 #include "AliPID.h"
+#include "AliVTrack.h"
 
-class AliVTrack;
+class TH2D;
 class TSpline3;
 
 class AliTPCPIDResponse: public TNamed {
@@ -28,7 +33,7 @@ public:
   AliTPCPIDResponse(const Double_t *param);
   AliTPCPIDResponse(const AliTPCPIDResponse&);
   AliTPCPIDResponse& operator=(const AliTPCPIDResponse&);
-  virtual ~AliTPCPIDResponse() {}
+  virtual ~AliTPCPIDResponse();
 
   enum EChamberStatus {
     kChamberOff=0,
@@ -80,25 +85,33 @@ public:
 
   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,
@@ -108,10 +121,11 @@ public:
                                  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, 
@@ -124,12 +138,6 @@ public:
                                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;
@@ -139,7 +147,12 @@ public:
                              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;
@@ -151,6 +164,26 @@ public:
   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)
@@ -172,17 +205,18 @@ private:
   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