Attached you can find TPC and EMCAL PID response files for the OADB and a larger...
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 16 Apr 2013 06:32:40 +0000 (06:32 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 16 Apr 2013 06:32:40 +0000 (06:32 +0000)
o preparation of multiplicity treatment for TPC PID
o updated TPC splines for (nearly) all periods
o first splines for LHC13b2_fix_1 (no eta maps, yet, to be updated...)
o update in multiplicity treatment for EMCAL PID
Jens Wiechula <Jens.Wiechula@cern.ch>

ANALYSIS/AliAnalysisTaskPIDResponse.cxx
ANALYSIS/AliAnalysisTaskPIDResponse.h
ANALYSIS/macros/AddTaskPIDResponse.C
STEER/AOD/AliAODpidUtil.cxx
STEER/ESD/AliESDpid.cxx
STEER/STEERBase/AliEMCALPIDResponse.cxx
STEER/STEERBase/AliEMCALPIDResponse.h
STEER/STEERBase/AliPIDResponse.cxx
STEER/STEERBase/AliPIDResponse.h
STEER/STEERBase/AliTPCPIDResponse.cxx
STEER/STEERBase/AliTPCPIDResponse.h

index 816e5c4..49af37b 100644 (file)
@@ -46,7 +46,8 @@ fRecoPass(0),
 fIsTunedOnData(kFALSE),
 fTunedOnDataMask(0),
 fRecoPassTuned(0),
-fUseTPCEtaCorrection(kFALSE)//TODO: In future, default kTRUE  
+fUseTPCEtaCorrection(kFALSE),//TODO: In future, default kTRUE 
+fUseTPCMultiplicityCorrection(kFALSE)//TODO: In future, default kTRUE  
 {
   //
   // Dummy constructor
@@ -67,7 +68,8 @@ fRecoPass(0),
 fIsTunedOnData(kFALSE),
 fTunedOnDataMask(0),
 fRecoPassTuned(0),
-fUseTPCEtaCorrection(kFALSE)//TODO: In future, default kTRUE
+fUseTPCEtaCorrection(kFALSE),//TODO: In future, default kTRUE
+fUseTPCMultiplicityCorrection(kFALSE)//TODO: In future, default kTRUE  
 {
   //
   // Default constructor
@@ -120,7 +122,7 @@ void AliAnalysisTaskPIDResponse::UserCreateOutputObjects()
       }
     }
     delete arr;
-  }
+  }  
 }
 
 //______________________________________________________________________________
@@ -135,9 +137,11 @@ void AliAnalysisTaskPIDResponse::UserExec(Option_t */*option*/)
   if (fRun!=fOldRun){
     SetRecoInfo();
     fOldRun=fRun;
+    
+    fPIDResponse->SetUseTPCEtaCorrection(fUseTPCEtaCorrection);
+    fPIDResponse->SetUseTPCMultiplicityCorrection(fUseTPCMultiplicityCorrection);
   }
 
-  fPIDResponse->SetUseTPCEtaCorrection(fUseTPCEtaCorrection);
   fPIDResponse->InitialiseEvent(event,fRecoPass);
   AliESDpid *pidresp = dynamic_cast<AliESDpid*>(fPIDResponse);
   if(pidresp && AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()){
@@ -178,6 +182,8 @@ void AliAnalysisTaskPIDResponse::SetRecoInfo()
     fPIDResponse->SetCurrentFile(fileName.Data());
   }
 
+  fPIDResponse->SetCurrentAliRootRev(prodInfo.GetAlirootSvnVersion());
+  
   if (prodInfo.IsMC() == kTRUE) fIsMC=kTRUE;         // protection if user didn't use macro switch
   if ( (prodInfo.IsMC() == kFALSE) && (fIsMC == kFALSE) ) {      // reco pass is needed only for data
     fRecoPass = prodInfo.GetRecoPass();
index 2e95651..d842ffd 100644 (file)
@@ -47,6 +47,9 @@ public:
   
   void SetUseTPCEtaCorrection(Bool_t useTPCEtaCorrection) { fUseTPCEtaCorrection = useTPCEtaCorrection; };
   Bool_t UseTPCEtaCorrection() const { return fUseTPCEtaCorrection; };
+  
+  void SetUseTPCMultiplicityCorrection(Bool_t useMultiplicityCorrection = kTRUE) { fUseTPCMultiplicityCorrection = useMultiplicityCorrection; };
+  Bool_t UseTPCMultiplicityCorrection() const { return fUseTPCMultiplicityCorrection; };
 
   void SetSpecialDetectorResponse(const char* det) { fSpecialDetResponse=det; }
 
@@ -65,7 +68,8 @@ private:
   Int_t   fTunedOnDataMask;            // mask to activate tuning on data on specific detectors
   Int_t   fRecoPassTuned;              // Reco pass tuned on data for MC
   
-  Bool_t  fUseTPCEtaCorrection;        // Use TPC eta correction
+  Bool_t fUseTPCEtaCorrection;          // Use TPC eta correction
+  Bool_t fUseTPCMultiplicityCorrection; // Use TPC multiplicity correction
   
   //
   void SetRecoInfo();
index a8fb7a9..b407a6f 100644 (file)
@@ -1,7 +1,8 @@
 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 useTPCEtaCorrection = kTRUE)
+                                    Bool_t useTPCEtaCorrection = kTRUE,
+                                    Bool_t useTPCMultiplicityCorrection = kFALSE)
 {
 // Macro to connect a centrality selection task to an existing analysis manager.
   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
@@ -22,10 +23,6 @@ AliAnalysisTask *AddTaskPIDResponse(Bool_t isMC=kFALSE, Bool_t autoMCesd=kTRUE,
     AliPIDResponseInputHandler *pidResponseIH = new AliPIDResponseInputHandler();
     multiInputHandler->AddInputEventHandler(pidResponseIH);
 
-    if (autoMCesd &&
-        multiInputHandler->GetFirstInputEventHandler()->IsA()==AliESDInputHandler::Class() &&
-        multiInputHandler->GetFirstMCEventHandler()
-       ) isMC=kTRUE;
     pidResponseIH->SetIsMC(isMC);
 
     return 0x0;
@@ -36,10 +33,6 @@ AliAnalysisTask *AddTaskPIDResponse(Bool_t isMC=kFALSE, Bool_t autoMCesd=kTRUE,
   printf("PIDResponse: Initialising AliAnalysisTaskPIDResponse\n");
   printf("========================================================================================\n");
   
-  if ( autoMCesd && (inputHandler->IsA() == AliESDInputHandler::Class()) ) {
-    isMC=mgr->GetMCtruthEventHandler()!=0x0;
-  }
-
   AliAnalysisTaskPIDResponse *pidTask = new AliAnalysisTaskPIDResponse("PIDResponseTask");
 //   pidTask->SelectCollisionCandidates(AliVEvent::kMB);
   pidTask->SetIsMC(isMC);
@@ -52,6 +45,7 @@ AliAnalysisTask *AddTaskPIDResponse(Bool_t isMC=kFALSE, Bool_t autoMCesd=kTRUE,
   pidTask->SetCachePID(cachePID);
   pidTask->SetSpecialDetectorResponse(detResponse);
   pidTask->SetUseTPCEtaCorrection(useTPCEtaCorrection);
+  pidTask->SetUseTPCMultiplicityCorrection(useTPCMultiplicityCorrection);
   mgr->AddTask(pidTask);
   
 //   AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("PIDResponseQA",
index f5d31fb..45cda46 100644 (file)
@@ -91,8 +91,10 @@ Float_t AliAODpidUtil::GetTPCsignalTunedOnData(const AliVTrack *t) const {
 
        if(kGood){
            //TODO maybe introduce different dEdxSources?
-        Double_t bethe = fTPCResponse.GetExpectedSignal(track, type, AliTPCPIDResponse::kdEdxDefault, this->UseTPCEtaCorrection());
-        Double_t sigma = fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, this->UseTPCEtaCorrection());
+        Double_t bethe = fTPCResponse.GetExpectedSignal(track, type, AliTPCPIDResponse::kdEdxDefault, this->UseTPCEtaCorrection(),
+                                                        this->UseTPCMultiplicityCorrection());
+        Double_t sigma = fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, this->UseTPCEtaCorrection(),
+                                                       this->UseTPCMultiplicityCorrection());
         dedx = gRandom->Gaus(bethe,sigma);
         
 //         if(iS == AliPID::ParticleCode(AliPID::kHe3) || iS == AliPID::ParticleCode(AliPID::kAlpha)) dedx *= 5;
index 7c39219..7904395 100644 (file)
@@ -123,8 +123,10 @@ Float_t AliESDpid::GetTPCsignalTunedOnData(const AliVTrack *t) const {
 
            if(kGood){
         //TODO maybe introduce different dEdxSources?
-        Double_t bethe = fTPCResponse.GetExpectedSignal(track, type, AliTPCPIDResponse::kdEdxDefault, this->UseTPCEtaCorrection());
-        Double_t sigma = fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, this->UseTPCEtaCorrection());
+        Double_t bethe = fTPCResponse.GetExpectedSignal(track, type, AliTPCPIDResponse::kdEdxDefault, this->UseTPCEtaCorrection(),
+                                                        this->UseTPCMultiplicityCorrection());
+        Double_t sigma = fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, this->UseTPCEtaCorrection(),
+                                                       this->UseTPCMultiplicityCorrection());
                dedx = gRandom->Gaus(bethe,sigma);
 //             if(iS == AliPID::ParticleCode(AliPID::kHe3) || iS == AliPID::ParticleCode(AliPID::kAlpha)) dedx *= 5;
            }
index 7efcb58..91c0618 100644 (file)
@@ -61,6 +61,7 @@ ClassImp(AliEMCALPIDResponse)
 AliEMCALPIDResponse::AliEMCALPIDResponse():
   TObject(),
   fNorm(NULL),
+  fCurrCentrality(-1.),
   fkPIDParams(NULL)
 {
   //
@@ -74,6 +75,7 @@ AliEMCALPIDResponse::AliEMCALPIDResponse():
 AliEMCALPIDResponse::AliEMCALPIDResponse(const AliEMCALPIDResponse &other):
   TObject(other),
   fNorm(other.fNorm),
+  fCurrCentrality(other.fCurrCentrality),
   fkPIDParams(other.fkPIDParams)
 {
   //
@@ -93,6 +95,7 @@ AliEMCALPIDResponse & AliEMCALPIDResponse::operator=( const AliEMCALPIDResponse&
   // Make copy
   TObject::operator=(other);
   fNorm = other.fNorm;
+  fCurrCentrality = other.fCurrCentrality;
   fkPIDParams = other.fkPIDParams;
 
  
@@ -264,6 +267,12 @@ const TVectorD* AliEMCALPIDResponse::GetParams(Int_t nParticle, Float_t fPt, Int
   // 6 = probLow  (not used for electrons)
   // 7 = probHigh (not used for electrons)
   //
+  // for PbPb the parametrization is done centrality dependent (marked by TString "Centrality")
+  // so first the correct centrality bin has to be found
+
+  // **** Centrality bins (hard coded for the moment)
+  const Int_t nCent = 7;
+  Int_t centBins[nCent+1] = {0,10,20,30,40,50,70,90};
 
   if(nParticle > AliPID::kSPECIES || nParticle <0) return NULL;
   if(nParticle == AliPID::kProton && charge == -1) nParticle = AliPID::kSPECIES; // special case for antiprotons
@@ -271,19 +280,53 @@ const TVectorD* AliEMCALPIDResponse::GetParams(Int_t nParticle, Float_t fPt, Int
   TObjArray * particlePar = dynamic_cast<TObjArray *>(fkPIDParams->At(nParticle));
   if(!particlePar) return NULL;
 
-  TIter parIter(particlePar);
   const TVectorD *parameters = NULL;
   Double_t momMin = 0.;
   Double_t momMax = 0.;
 
-  while((parameters = static_cast<const TVectorD *>(parIter()))){
+  // is the centrality dependent parametrization used
+  TString arrayName = particlePar->GetName();
+
+  // centrality dependent parametrization
+  if(arrayName.Contains("Centrality")){
+    
+    for(Int_t iCent = 0; iCent < nCent; iCent++){
 
-    momMin = (*parameters)[0];
-    momMax = (*parameters)[1];
+      if( fCurrCentrality > centBins[iCent] && fCurrCentrality < centBins[iCent+1] ){
+       
+       TObjArray * centPar = dynamic_cast<TObjArray *>(particlePar->At(iCent));
+       if(!centPar) return NULL;
+       
+       TIter centIter(centPar);
+       parameters = NULL;
+       momMin = 0.;
+       momMax = 0.;
+       
+       while((parameters = static_cast<const TVectorD *>(centIter()))){
+         
+         momMin = (*parameters)[0];
+         momMax = (*parameters)[1];
+         
+         if( fPt > momMin && fPt < momMax ) return parameters;
+        
+       } 
+      }
+    }
+  }
 
-    if( fPt > momMin && fPt < momMax ) return parameters;
+  // NO centrality dependent parametrization
+  else{
 
-  }  
+    TIter parIter(particlePar);
+    while((parameters = static_cast<const TVectorD *>(parIter()))){
+      
+      momMin = (*parameters)[0];
+      momMax = (*parameters)[1];
+      
+      if( fPt > momMin && fPt < momMax ) return parameters;
+      
+    }  
+  }
   AliDebug(2, Form("NO params for particle %d and momentum %f \n", nParticle, fPt));
 
   return parameters;
index 57c17fb..ee5b930 100644 (file)
@@ -35,6 +35,7 @@ public :
   
     //Setters
     void   SetPIDParams(const TObjArray * params) { fkPIDParams = params; }
+    void   SetCentrality(Float_t currentCentrality) { fCurrCentrality = currentCentrality;}
     
 
     // EMCAL probability
@@ -46,11 +47,13 @@ private:
 
   TF1 *fNorm;                            // Gauss function for normalizing NON electron probabilities 
 
+  Double_t fCurrCentrality;              // current (in the current event) centrality percentile 
+
   const TObjArray *fkPIDParams;               // PID Params
 
   const TVectorD* GetParams(Int_t nParticle, Float_t fPt, Int_t charge) const; 
 
-  ClassDef(AliEMCALPIDResponse, 1)
+  ClassDef(AliEMCALPIDResponse, 2)
 };
 
 #endif // #ifdef AliEMCALPIDResponse_cxx
index 4806cc9..6875616 100644 (file)
@@ -71,6 +71,7 @@ fLHCperiod(),
 fMCperiodTPC(),
 fMCperiodUser(),
 fCurrentFile(),
+fCurrentAliRootRev(-1),
 fRecoPass(0),
 fRecoPassUser(-1),
 fRun(-1),
@@ -81,7 +82,8 @@ fResT0AC(55.),
 fArrPidResponseMaster(NULL),
 fResolutionCorrection(NULL),
 fOADBvoltageMaps(NULL),
-fUseTPCEtaCorrection(kFALSE),//TODO: In future, default kTRUE
+fUseTPCEtaCorrection(kFALSE),
+fUseTPCMultiplicityCorrection(kFALSE),
 fTRDPIDResponseObject(NULL),
 fTOFtail(1.1),
 fTOFPIDParams(NULL),
@@ -132,6 +134,7 @@ fLHCperiod(),
 fMCperiodTPC(),
 fMCperiodUser(other.fMCperiodUser),
 fCurrentFile(),
+fCurrentAliRootRev(other.fCurrentAliRootRev),
 fRecoPass(0),
 fRecoPassUser(other.fRecoPassUser),
 fRun(-1),
@@ -143,6 +146,7 @@ fArrPidResponseMaster(NULL),
 fResolutionCorrection(NULL),
 fOADBvoltageMaps(NULL),
 fUseTPCEtaCorrection(other.fUseTPCEtaCorrection),
+fUseTPCMultiplicityCorrection(other.fUseTPCMultiplicityCorrection),
 fTRDPIDResponseObject(NULL),
 fTOFtail(1.1),
 fTOFPIDParams(NULL),
@@ -184,6 +188,7 @@ AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other)
     fMCperiodTPC="";
     fMCperiodUser=other.fMCperiodUser;
     fCurrentFile="";
+    fCurrentAliRootRev=other.fCurrentAliRootRev;
     fRecoPass=0;
     fRecoPassUser=other.fRecoPassUser;
     fRun=-1;
@@ -195,12 +200,14 @@ AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other)
     fResolutionCorrection=NULL;
     fOADBvoltageMaps=NULL;
     fUseTPCEtaCorrection=other.fUseTPCEtaCorrection;
+    fUseTPCMultiplicityCorrection=other.fUseTPCMultiplicityCorrection;
     fTRDPIDResponseObject=NULL;
     fEMCALPIDParams=NULL;
     fTOFtail=1.1;
     fTOFPIDParams=NULL;
     fHMPIDPIDParams=NULL;
     fCurrentEvent=other.fCurrentEvent;
+
   }
   return *this;
 }
@@ -271,8 +278,7 @@ Float_t AliPIDResponse::NumberOfSigmasTPC( const AliVParticle *vtrack,
   //get number of sigmas according the selected TPC gain configuration scenario
   const AliVTrack *track=static_cast<const AliVTrack*>(vtrack);
 
-//   return 0.;
-  Float_t nSigma=fTPCResponse.GetNumberOfSigmas(track, type, dedxSource, fUseTPCEtaCorrection);
+  Float_t nSigma=fTPCResponse.GetNumberOfSigmas(track, type, dedxSource, fUseTPCEtaCorrection, fUseTPCMultiplicityCorrection);
 
   return nSigma;
 }
@@ -544,6 +550,20 @@ void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass, Int_t run)
     fTPCResponse.SetSigma(3.79301e-03*corrSigma, 2.21280e+04);
   }
   
+  // Set up TPC multiplicity for PbPb
+  //TODO Will NOT give the desired number for AODs -> Needs new variable/function in future.
+  // Fatal, if AOD event and correction enabled
+  //printf("DETECTED class: %s (%d)\n\n\n\n", event->IsA()->GetName(), fUseTPCMultiplicityCorrection);//TODO
+  if (fUseTPCMultiplicityCorrection && strcmp(event->IsA()->GetName(), "AliESDEvent") != 0) {
+    AliFatal("TPC multiplicity correction is enabled, but will NOT work for AOD events, only for ESD => Disabled multiplicity correction!");
+    fUseTPCMultiplicityCorrection = kFALSE;
+  }
+  
+  if (fUseTPCMultiplicityCorrection)
+    fTPCResponse.SetCurrentEventMultiplicity(event->GetNumberOfTracks());
+  else
+    fTPCResponse.SetCurrentEventMultiplicity(0);
+  
   //TOF resolution
   SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)fTOFPIDParams->GetStartTimeMethod());
 
@@ -556,6 +576,10 @@ void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass, Int_t run)
   else{
     fCurrCentrality = -1;
   }
+
+  // Set centrality percentile for EMCAL
+  fEMCALResponse.SetCentrality(fCurrCentrality);
+
 }
 
 //______________________________________________________________________________
@@ -661,8 +685,18 @@ void AliPIDResponse::SetRecoInfo()
 //   if (fRun >= 188167 && fRun <= 188355 ) { fLHCperiod="LHC12G"; fBeamType="PP"; /*fMCperiodTPC="";*/ }
 //   if (fRun >= 188356 && fRun <= 188503 ) { fLHCperiod="LHC12G"; fBeamType="PPB"; /*fMCperiodTPC="";*/ }
 // for the moment use 12g parametrisation for all full gain runs (LHC12f+)
-  if (fRun >= 186636  && fRun < 194480) { fLHCperiod="LHC12G"; fBeamType="PPB"; fMCperiodTPC="LHC12G"; }
-  if (fRun >= 194480) { fLHCperiod="LHC13B"; fBeamType="PPB"; fMCperiodTPC="LHC12G"; }
+  if (fRun >= 186636 && fRun < 194480) { fLHCperiod="LHC12G"; fBeamType="PPB"; fMCperiodTPC="LHC12G"; }
+
+  // New parametrisation for 2013 pPb runs
+  if (fRun >= 194480) { 
+    fLHCperiod="LHC13B"; 
+    fBeamType="PPB";
+  
+    if (fCurrentAliRootRev >= 61605)
+      fMCperiodTPC="LHC13B2_FIX";
+    else
+      fMCperiodTPC="LHC12G";
+  }
 
   //exception new pp MC productions from 2011
   if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) { fMCperiodTPC="LHC11B2"; fBeamType="PP"; }
@@ -881,13 +915,13 @@ void AliPIDResponse::SetTPCEtaMaps(Double_t refineFactorMapX, Double_t refineFac
   TString period = fLHCperiod.IsNull() ? "No period information" : fLHCperiod;
   
   if (fIsMC)  {
-    if (!fTuneMConData) {
+    if (!(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))) {
       period=fMCperiodTPC;
       dataType="MC";
     }
     fRecoPass = 1;
     
-    if (!fTuneMConData && fMCperiodTPC.IsNull()) {
+    if (!(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC)) && fMCperiodTPC.IsNull()) {
       AliFatal("MC detected, but no MC period set -> Not changing eta maps!");
       return;
     }
@@ -912,13 +946,14 @@ void AliPIDResponse::SetTPCEtaMaps(Double_t refineFactorMapX, Double_t refineFac
                                               Form("TPCetaMaps_%s_pass%d", dataType.Data(), recopass));
   if (statusCont) {
     AliError("Failed initializing TPC eta correction maps from OADB -> Disabled eta correction");
+    fUseTPCEtaCorrection = kFALSE;
   }
   else {
     AliInfo(Form("Loading TPC eta correction map from %s/COMMON/PID/data/TPCetaMaps.root", fOADBPath.Data()));
     
     TH2D* etaMap = 0x0;
     
-    if (fIsMC && !fTuneMConData) {
+    if (fIsMC && !(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))) {
       TString searchMap = Form("TPCetaMaps_%s_%s_pass%d", dataType.Data(), period.Data(), recopass);
       etaMap = dynamic_cast<TH2D *>(etaMapsCont.GetDefaultObject(searchMap.Data()));
       if (!etaMap) {
@@ -933,6 +968,7 @@ void AliPIDResponse::SetTPCEtaMaps(Double_t refineFactorMapX, Double_t refineFac
         
     if (!etaMap) {
       AliError(Form("TPC eta correction map not found for run %d and also no default map found -> Disabled eta correction!!!", fRun));
+      fUseTPCEtaCorrection = kFALSE;
     }
     else {
       TH2D* etaMapRefined = RefineHistoViaLinearInterpolation(etaMap, refineFactorMapX, refineFactorMapY);
@@ -941,6 +977,7 @@ void AliPIDResponse::SetTPCEtaMaps(Double_t refineFactorMapX, Double_t refineFac
         if (!fTPCResponse.SetEtaCorrMap(etaMapRefined)) {
           AliError(Form("Failed to set TPC eta correction map for run %d -> Disabled eta correction!!!", fRun));
           fTPCResponse.SetEtaCorrMap(0x0);
+          fUseTPCEtaCorrection = kFALSE;
         }
         else {
           AliInfo(Form("Loaded TPC eta correction map (refine factors %.2f/%.2f) from %s/COMMON/PID/data/TPCetaMaps.root: %s", 
@@ -951,10 +988,17 @@ void AliPIDResponse::SetTPCEtaMaps(Double_t refineFactorMapX, Double_t refineFac
       }
       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));
+        fUseTPCEtaCorrection = kFALSE;
       }
     }
   }
   
+  // If there was some problem loading the eta maps, it makes no sense to load the sigma maps (that require eta corrected data)
+  if (fUseTPCEtaCorrection == kFALSE) {
+    AliError("Failed to load TPC eta correction map required by sigma maps -> Using old parametrisation for sigma"); 
+    return;
+  }
+  
   // Load the sigma parametrisation (1/dEdx vs tanTheta_local (~eta))
   AliOADBContainer etaSigmaMapsCont(Form("TPCetaSigmaMaps_%s_pass%d", dataType.Data(), recopass)); 
   
@@ -968,7 +1012,7 @@ void AliPIDResponse::SetTPCEtaMaps(Double_t refineFactorMapX, Double_t refineFac
     
     TObjArray* etaSigmaPars = 0x0;
     
-    if (fIsMC && !fTuneMConData) {
+    if (fIsMC && !(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))) {
       TString searchMap = Form("TPCetaSigmaMaps_%s_%s_pass%d", dataType.Data(), period.Data(), recopass);
       etaSigmaPars = dynamic_cast<TObjArray *>(etaSigmaMapsCont.GetDefaultObject(searchMap.Data()));
       if (!etaSigmaPars) {
@@ -1091,16 +1135,16 @@ void AliPIDResponse::SetTPCParametrisation()
   TString datatype="DATA";
   //in case of mc fRecoPass is per default 1
   if (fIsMC) {
-      if(!fTuneMConData) datatype="MC";
+      if(!(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))) datatype="MC";
       fRecoPass=1;
   }
 
   // period
   TString period=fLHCperiod;
-  if (fIsMC && !fTuneMConData) period=fMCperiodTPC;
+  if (fIsMC && !(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))) period=fMCperiodTPC;
 
   Int_t recopass = fRecoPass;
-  if (fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC) ) recopass = fRecoPassUser;
+  if(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC)) recopass = fRecoPassUser;
     
   AliInfo(Form("Searching splines for: %s %s PASS%d %s",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
   Bool_t found=kFALSE;
@@ -1223,8 +1267,86 @@ void AliPIDResponse::SetTPCParametrisation()
     AliError(Form("No splines found for: %s %s PASS%d %s",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
   }
 
+
+  //
+  // Setup multiplicity correction
+  //
+  if (fUseTPCMultiplicityCorrection && (fBeamType.CompareTo("PBPB") == 0 || fBeamType.CompareTo("AA") == 0)) {
+    AliInfo("Multiplicity correction enabled!");
+    
+    //TODO After testing, load parameters from outside       
+    if (period.Contains("LHC11A10"))  {//LHC11A10A
+      AliInfo("Using multiplicity correction parameters for 11a10!");
+      fTPCResponse.SetParameterMultiplicityCorrection(0, 6.90133e-06);
+      fTPCResponse.SetParameterMultiplicityCorrection(1, -1.22123e-03);
+      fTPCResponse.SetParameterMultiplicityCorrection(2, 1.80220e-02);
+      fTPCResponse.SetParameterMultiplicityCorrection(3, 0.1);
+      fTPCResponse.SetParameterMultiplicityCorrection(4, 6.45306e-03);
+      
+      fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(0, -2.85505e-07);
+      fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(1, -1.31911e-06);
+      fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(2, -0.5);
+
+      fTPCResponse.SetParameterMultiplicitySigmaCorrection(0, -4.29665e-05);
+      fTPCResponse.SetParameterMultiplicitySigmaCorrection(1, 1.37023e-02);
+      fTPCResponse.SetParameterMultiplicitySigmaCorrection(2, -6.36337e-01);
+      fTPCResponse.SetParameterMultiplicitySigmaCorrection(3, 1.13479e-02);
+    }
+    else if (period.Contains("LHC10H") && recopass == 2) {    
+      AliInfo("Using multiplicity correction parameters for 10h.pass2!");
+      fTPCResponse.SetParameterMultiplicityCorrection(0, 3.21636e-07);
+      fTPCResponse.SetParameterMultiplicityCorrection(1, -6.65876e-04);
+      fTPCResponse.SetParameterMultiplicityCorrection(2, 1.28786e-03);
+      fTPCResponse.SetParameterMultiplicityCorrection(3, 1.47677e-02);
+      fTPCResponse.SetParameterMultiplicityCorrection(4, 0);
+      
+      fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(0, 7.23591e-08);
+      fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(1, 2.7469e-06);
+      fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(2, -0.5);
+      
+      fTPCResponse.SetParameterMultiplicitySigmaCorrection(0, -1.22590e-05);
+      fTPCResponse.SetParameterMultiplicitySigmaCorrection(1, 6.88888e-03);
+      fTPCResponse.SetParameterMultiplicitySigmaCorrection(2, -3.20788e-01);
+      fTPCResponse.SetParameterMultiplicitySigmaCorrection(3, 1.07345e-02);
+    }
+    else {
+      AliError(Form("Multiplicity correction is enabled, but no multiplicity correction parameters have been found for period %s.pass%d -> Mulitplicity correction DISABLED!", period.Data(), recopass));
+      fUseTPCMultiplicityCorrection = kFALSE;
+      fTPCResponse.ResetMultiplicityCorrectionFunctions();
+    }
+  }
+  else {
+    // Just set parameters such that overall correction factor is 1, i.e. no correction.
+    // This is just a reasonable choice for the parameters for safety reasons. Disabling
+    // the multiplicity correction will anyhow skip the calculation of the corresponding
+    // correction factor inside THIS class. Nevertheless, experts can access the TPCPIDResponse
+    // directly and use it for calculations - which should still give valid results, even if
+    // the multiplicity correction is explicitely enabled in such expert calls.
+    
+    AliInfo(Form("Multiplicity correction %sdisabled (%s)!", fUseTPCMultiplicityCorrection ? "automatically " : "",
+                 fUseTPCMultiplicityCorrection ? "no PbPb or AA" : "requested by user"));
+    
+    fUseTPCMultiplicityCorrection = kFALSE;
+    fTPCResponse.ResetMultiplicityCorrectionFunctions();
+  }
+  
+  /*
+  //TODO NOW start
+  for (Int_t i = 0; i <= 4 + 1; i++) {
+    printf("parMultCorr: %d, %e\n", i, fTPCResponse.GetMultiplicityCorrectionFunction()->GetParameter(i));
+  }
+  for (Int_t j = 0; j <= 2 + 1; j++) {
+    printf("parMultCorrTanTheta: %d, %e\n", j, fTPCResponse.GetMultiplicityCorrectionFunctionTanTheta()->GetParameter(j));
+  }
+  for (Int_t j = 0; j <= 3 + 1; j++) {
+    printf("parMultSigmaCorr: %d, %e\n", j, fTPCResponse.GetMultiplicitySigmaCorrectionFunction()->GetParameter(j));
+  }
+  
+  //TODO NOW end
+  */
+  
   //
-  // Setup resolution parametrisation
+  // Setup old resolution parametrisation
   //
   
   //default
@@ -1805,9 +1927,10 @@ Float_t AliPIDResponse::GetNumberOfSigmasTPC(const AliVParticle *vtrack, AliPID:
   // the following call is needed in order to fill the transient data member
   // fTPCsignalTuned which is used in the TPCPIDResponse to judge
   // if using tuned on data
-  if (fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC) ) this->GetTPCsignalTunedOnData(track);
+  if (fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))
+    this->GetTPCsignalTunedOnData(track);
   
-  return fTPCResponse.GetNumberOfSigmas(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
+  return fTPCResponse.GetNumberOfSigmas(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection, fUseTPCMultiplicityCorrection);
 }
 
 //______________________________________________________________________________
@@ -1886,10 +2009,10 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::GetSignalDeltaTPC(const AliVPartic
   // the following call is needed in order to fill the transient data member
   // fTPCsignalTuned which is used in the TPCPIDResponse to judge
   // if using tuned on data
-  if (fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC) ) 
+  if (fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))
     this->GetTPCsignalTunedOnData(track);
   
-  val=fTPCResponse.GetSignalDelta(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection, ratio);
+  val=fTPCResponse.GetSignalDelta(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection, fUseTPCMultiplicityCorrection, ratio);
   
   return GetTPCPIDStatus(track);
 }
@@ -1902,6 +2025,7 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::GetSignalDeltaTOF(const AliVPartic
   //
   AliVTrack *track=(AliVTrack*)vtrack;
   val=GetSignalDeltaTOFold(track, type, ratio);
+  
   return GetTOFPIDStatus(track);
 }
 
@@ -2005,7 +2129,7 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTPCProbability  (const A
   Double_t dedx=track->GetTPCsignal();
   Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
   
-  if (fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC) ) dedx = this->GetTPCsignalTunedOnData(track);
+  if (fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC)) dedx = this->GetTPCsignalTunedOnData(track);
   
   Double_t bethe = 0.;
   Double_t sigma = 0.;
@@ -2013,8 +2137,8 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTPCProbability  (const A
   for (Int_t j=0; j<nSpecies; j++) {
     AliPID::EParticleType type=AliPID::EParticleType(j);
     
-    bethe=fTPCResponse.GetExpectedSignal(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
-    sigma=fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection);
+    bethe=fTPCResponse.GetExpectedSignal(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection, fUseTPCMultiplicityCorrection);
+    sigma=fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection, fUseTPCMultiplicityCorrection);
     
     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
index 97cbe5d..fa00d13 100644 (file)
@@ -130,6 +130,9 @@ public:
   
   void InitialiseEvent(AliVEvent *event, Int_t pass, Int_t run=-1);
   void SetCurrentFile(const char* file) { fCurrentFile=file; }
+  
+  void SetCurrentAliRootRev(Int_t alirootRev) { fCurrentAliRootRev = alirootRev; }
+  Int_t GetCurrentAliRootRev() const { return fCurrentAliRootRev; }
 
   // cache PID in the track
   void SetCachePID(Bool_t cache)    { fCachePID=cache;  }
@@ -150,6 +153,9 @@ public:
   void SetUseTPCEtaCorrection(Bool_t useEtaCorrection = kTRUE) { fUseTPCEtaCorrection = useEtaCorrection; };
   Bool_t UseTPCEtaCorrection() const { return fUseTPCEtaCorrection; };
   
+  void SetUseTPCMultiplicityCorrection(Bool_t useMultiplicityCorrection = kTRUE) { fUseTPCMultiplicityCorrection = useMultiplicityCorrection; };
+  Bool_t UseTPCMultiplicityCorrection() const { return fUseTPCMultiplicityCorrection; };
+  
   // 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);
@@ -199,6 +205,7 @@ private:
   TString fMCperiodTPC;                //! corresponding MC period to use for the TPC splines
   TString fMCperiodUser;               //  MC prodution requested by the user
   TString fCurrentFile;                //! name of currently processed file
+  Int_t   fCurrentAliRootRev;          //! Aliroot rev. used to reconstruct the data
   Int_t   fRecoPass;                   //! reconstruction pass
   Int_t   fRecoPassUser;               //  reconstruction pass explicitly set by the user
   Int_t   fRun;                        //! current run number
@@ -207,10 +214,11 @@ private:
   Float_t fResT0C;                     //! T0C resolution in current run
   Float_t fResT0AC;                    //! T0A.and.T0C resolution in current run
   
-  TObjArray *fArrPidResponseMaster;    //!  TPC pid splines
-  TF1       *fResolutionCorrection;    //! TPC resolution correction
-  AliOADBContainer* fOADBvoltageMaps;  //! container with the voltage maps
-  Bool_t fUseTPCEtaCorrection;         // Use TPC eta correction
+  TObjArray *fArrPidResponseMaster;     //!  TPC pid splines
+  TF1       *fResolutionCorrection;     //! TPC resolution correction
+  AliOADBContainer* fOADBvoltageMaps;   //! container with the voltage maps
+  Bool_t fUseTPCEtaCorrection;          // Use TPC eta correction
+  Bool_t fUseTPCMultiplicityCorrection; // Use TPC multiplicity correction
 
   AliTRDPIDResponseObject *fTRDPIDResponseObject; //! TRD PID Response Object
 
@@ -305,7 +313,7 @@ private:
   EDetPidStatus GetPHOSPIDStatus(const AliVTrack *track) const;
   EDetPidStatus GetEMCALPIDStatus(const AliVTrack *track) const;
 
-  ClassDef(AliPIDResponse, 11);  //PID response handling
+  ClassDef(AliPIDResponse, 12);  //PID response handling
 };
 
 #endif
index bbf370f..367d8cc 100644 (file)
@@ -71,12 +71,26 @@ AliTPCPIDResponse::AliTPCPIDResponse():
   fMagField(0.),
   fhEtaCorr(0x0),
   fhEtaSigmaPar1(0x0),
-  fSigmaPar0(0.0)
+  fSigmaPar0(0.0),
+  fCurrentEventMultiplicity(0),
+  fCorrFuncMultiplicity(0x0),
+  fCorrFuncMultiplicityTanTheta(0x0),
+  fCorrFuncSigmaMultiplicity(0x0)
 {
   //
   //  The default constructor
   //
   for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=0.07;fResN2[i]=0.0;}
+  
+  fCorrFuncMultiplicity = new TF1("fCorrFuncMultiplicity", 
+                                  "[0] + [1]*TMath::Max([4], TMath::Min(x, [3])) + [2] * TMath::Power(TMath::Max([4], TMath::Min(x, [3])), 2)",
+                                  0., 0.2);
+  fCorrFuncMultiplicityTanTheta = new TF1("fCorrFuncMultiplicityTanTheta", "[0] * (x - [2]) + [1] * (x * x - [2] * [2])", -1.5, 1.5);
+  fCorrFuncSigmaMultiplicity = new TF1("fCorrFuncSigmaMultiplicity",
+                                       "TMath::Max(0, [0] + [1]*TMath::Min(x, [3]) + [2] * TMath::Power(TMath::Min(x, [3]), 2))", 0., 0.2);
+
+  
+  ResetMultiplicityCorrectionFunctions();
 }
 /*TODO remove?
 //_________________________________________________________________________
@@ -122,6 +136,15 @@ AliTPCPIDResponse::~AliTPCPIDResponse()
   
   delete fhEtaSigmaPar1;
   fhEtaSigmaPar1 = 0x0;
+  
+  delete fCorrFuncMultiplicity;
+  fCorrFuncMultiplicity = 0x0;
+  
+  delete fCorrFuncMultiplicityTanTheta;
+  fCorrFuncMultiplicityTanTheta = 0x0;
+  
+  delete fCorrFuncSigmaMultiplicity;
+  fCorrFuncSigmaMultiplicity = 0x0;
 }
 
 
@@ -147,7 +170,11 @@ AliTPCPIDResponse::AliTPCPIDResponse(const AliTPCPIDResponse& that):
   fMagField(that.fMagField),
   fhEtaCorr(0x0),
   fhEtaSigmaPar1(0x0),
-  fSigmaPar0(that.fSigmaPar0)
+  fSigmaPar0(that.fSigmaPar0),
+  fCurrentEventMultiplicity(that.fCurrentEventMultiplicity),
+  fCorrFuncMultiplicity(0x0),
+  fCorrFuncMultiplicityTanTheta(0x0),
+  fCorrFuncSigmaMultiplicity(0x0)
 {
   //copy ctor
   for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=that.fRes0[i];fResN2[i]=that.fResN2[i];}
@@ -162,6 +189,19 @@ AliTPCPIDResponse::AliTPCPIDResponse(const AliTPCPIDResponse& that):
     fhEtaSigmaPar1 = new TH2D(*(that.fhEtaSigmaPar1));
     fhEtaSigmaPar1->SetDirectory(0);
   }
+  
+  // Copy multiplicity correction functions
+  if (that.fCorrFuncMultiplicity) {
+    fCorrFuncMultiplicity = new TF1(*(that.fCorrFuncMultiplicity));
+  }
+  
+  if (that.fCorrFuncMultiplicityTanTheta) {
+    fCorrFuncMultiplicityTanTheta = new TF1(*(that.fCorrFuncMultiplicityTanTheta));
+  }
+  
+  if (that.fCorrFuncSigmaMultiplicity) {
+    fCorrFuncSigmaMultiplicity = new TF1(*(that.fCorrFuncSigmaMultiplicity));
+  }
 }
 
 //_________________________________________________________________________
@@ -185,6 +225,7 @@ AliTPCPIDResponse& AliTPCPIDResponse::operator=(const AliTPCPIDResponse& that)
   fBadOROCthreshhold=that.fBadOROCthreshhold;
   fMaxBadLengthFraction=that.fMaxBadLengthFraction;
   fMagField=that.fMagField;
+  fCurrentEventMultiplicity=that.fCurrentEventMultiplicity;
   for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=that.fRes0[i];fResN2[i]=that.fResN2[i];}
 
   delete fhEtaCorr;
@@ -202,6 +243,24 @@ AliTPCPIDResponse& AliTPCPIDResponse::operator=(const AliTPCPIDResponse& that)
   }
   
   fSigmaPar0 = that.fSigmaPar0;
+  
+  delete fCorrFuncMultiplicity;
+  fCorrFuncMultiplicity = 0x0;
+  if (that.fCorrFuncMultiplicity) {
+    fCorrFuncMultiplicity = new TF1(*(that.fCorrFuncMultiplicity));
+  }
+  
+  delete fCorrFuncMultiplicityTanTheta;
+  fCorrFuncMultiplicityTanTheta = 0x0;
+  if (that.fCorrFuncMultiplicityTanTheta) {
+    fCorrFuncMultiplicityTanTheta = new TF1(*(that.fCorrFuncMultiplicityTanTheta));
+  }
+  
+  delete fCorrFuncSigmaMultiplicity;
+  fCorrFuncSigmaMultiplicity = 0x0;
+  if (that.fCorrFuncSigmaMultiplicity) {
+    fCorrFuncSigmaMultiplicity = new TF1(*(that.fCorrFuncSigmaMultiplicity));
+  }
 
   return *this;
 }
@@ -256,7 +315,7 @@ Double_t AliTPCPIDResponse::GetExpectedSignal(const Float_t mom,
   //
   // Deprecated function (for backward compatibility). Please use 
   // GetExpectedSignal(const AliVTrack* track, AliPID::EParticleType species, ETPCdEdxSource dedxSource,
-  //                   Bool_t correctEta = kTRUE);
+  //                   Bool_t correctEta, Bool_t correctMultiplicity);
   // instead!
   //
   //
@@ -322,14 +381,16 @@ Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track,
                                               AliPID::EParticleType species,
                                               Double_t /*dEdx*/,
                                               const TSpline3* responseFunction,
-                                              Bool_t correctEta) const 
+                                              Bool_t correctEta,
+                                              Bool_t correctMultiplicity) const 
 {
   // Calculates the expected PID signal as the function of 
   // 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 plus, if desired, taking into account the eta dependence. 
+  // Bethe-Bloch formula plus, if desired, taking into account the eta dependence
+  // and the multiplicity dependence (for PbPb). 
   // This can be improved. By taking into account the number of
   // assigned clusters and/or the track dip angle, for example.  
   //
@@ -345,13 +406,23 @@ Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track,
     return Bethe(mom/mass) * chargeFactor;
   
   Double_t dEdxSplines = fMIP*responseFunction->Eval(mom/mass) * chargeFactor;
-  
-  if (!correctEta)
+    
+  if (!correctEta && !correctMultiplicity)
     return dEdxSplines;
   
-  //TODO Alternatively take current track dEdx
-  //return dEdxSplines * GetEtaCorrection(track, dEdx);
-  return dEdxSplines * GetEtaCorrection(track, dEdxSplines);
+  Double_t corrFactorEta = 1.0;
+  Double_t corrFactorMultiplicity = 1.0;
+  
+  if (correctEta) {
+    corrFactorEta = GetEtaCorrection(track, dEdxSplines);
+    //TODO Alternatively take current track dEdx
+    //corrFactorEta = GetEtaCorrection(track, dEdx);
+  }
+  
+  if (correctMultiplicity)
+    corrFactorMultiplicity = GetMultiplicityCorrection(track, dEdxSplines * corrFactorEta, fCurrentEventMultiplicity);
+
+  return dEdxSplines * corrFactorEta * corrFactorMultiplicity;
 }
 
 
@@ -359,13 +430,15 @@ Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track,
 Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track,
                                               AliPID::EParticleType species,
                                               ETPCdEdxSource dedxSource,
-                                              Bool_t correctEta) const
+                                              Bool_t correctEta,
+                                              Bool_t correctMultiplicity) 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 plus, if desired, taking into account the eta dependence. 
+  // Bethe-Bloch formula plus, if desired, taking into account the eta dependence
+  // and the multiplicity dependence (for PbPb). 
   // This can be improved. By taking into account the number of
   // assigned clusters and/or the track dip angle, for example.  
   //
@@ -389,7 +462,7 @@ Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track,
   }
   
   // Charge factor already taken into account inside the following function call
-  return GetExpectedSignal(track, species, dEdx, responseFunction, correctEta);
+  return GetExpectedSignal(track, species, dEdx, responseFunction, correctEta, correctMultiplicity);
 }
   
 //_________________________________________________________________________
@@ -451,7 +524,8 @@ Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track,
                                              Double_t dEdx,
                                              Int_t nPoints,
                                              const TSpline3* responseFunction,
-                                             Bool_t correctEta) const 
+                                             Bool_t correctEta,
+                                             Bool_t correctMultiplicity) const 
 {
   // Calculates the expected sigma of the PID signal as the function of 
   // the information stored in the track and the given parameters,
@@ -460,20 +534,37 @@ Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track,
   
   if (!responseFunction)
     return 999;
-  
+    
+  //TODO Check whether it makes sense to set correctMultiplicity to kTRUE while correctEta might be kFALSE
   // 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) *
+      return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, correctMultiplicity) *
                fRes0[gainScenario] * sqrt(1. + fResN2[gainScenario]/nPoints);
     else
-      return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE)*fRes0[gainScenario];
+      return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, correctMultiplicity)*fRes0[gainScenario];
   }
     
   if (nPoints > 0) {
+    // Use eta correction (+ eta-dependent sigma)
     Double_t sigmaPar1 = GetSigmaPar1(track, species, dEdx, responseFunction);
     
-    return GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE)*sqrt(fSigmaPar0 * fSigmaPar0 + sigmaPar1 * sigmaPar1 / nPoints);
+    if (correctMultiplicity) {
+      // In addition, take into account multiplicity dependence of mean and sigma of dEdx
+      Double_t dEdxExpectedEtaCorrected = GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE);
+      
+      // GetMultiplicityCorrection and GetMultiplicitySigmaCorrection both need the eta corrected dEdxExpected
+      Double_t multiplicityCorrFactor = GetMultiplicityCorrection(track, dEdxExpectedEtaCorrected, fCurrentEventMultiplicity);
+      Double_t multiplicitySigmaCorrFactor = GetMultiplicitySigmaCorrection(dEdxExpectedEtaCorrected, fCurrentEventMultiplicity);
+      
+      // multiplicityCorrFactor to correct dEdxExpected for multiplicity. In addition: Correction factor for sigma
+      return (dEdxExpectedEtaCorrected * multiplicityCorrFactor) 
+              * (sqrt(fSigmaPar0 * fSigmaPar0 + sigmaPar1 * sigmaPar1 / nPoints) * multiplicitySigmaCorrFactor);
+    }
+    else {
+      return GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE)*
+             sqrt(fSigmaPar0 * fSigmaPar0 + sigmaPar1 * sigmaPar1 / nPoints);
+    }
   }
   else { 
     // One should never have/take tracks with 0 dEdx clusters!
@@ -486,7 +577,8 @@ Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track,
 Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track, 
                                              AliPID::EParticleType species,
                                              ETPCdEdxSource dedxSource,
-                                             Bool_t correctEta) const 
+                                             Bool_t correctEta,
+                                             Bool_t correctMultiplicity) const 
 {
   // Calculates the expected sigma of the PID signal as the function of 
   // the information stored in the track, for the specified particle type 
@@ -501,7 +593,7 @@ Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track,
   if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
     return 999; //TODO: better handling!
   
-  return GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta);
+  return GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta, correctMultiplicity);
 }
 
 
@@ -509,7 +601,8 @@ Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track,
 Float_t AliTPCPIDResponse::GetNumberOfSigmas(const AliVTrack* track, 
                              AliPID::EParticleType species,
                              ETPCdEdxSource dedxSource,
-                             Bool_t correctEta) const
+                             Bool_t correctEta,
+                             Bool_t correctMultiplicity) const
 {
   //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
@@ -523,8 +616,8 @@ Float_t AliTPCPIDResponse::GetNumberOfSigmas(const AliVTrack* track,
   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);
+  Double_t bethe = GetExpectedSignal(track, species, dEdx, responseFunction, correctEta, correctMultiplicity);
+  Double_t sigma = GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta, correctMultiplicity);
   // 999 will be returned by GetExpectedSigma e.g. in case of 0 dEdx clusters
   if (sigma >= 998) 
     return -999;
@@ -537,7 +630,8 @@ Float_t AliTPCPIDResponse::GetSignalDelta(const AliVTrack* track,
                                           AliPID::EParticleType species,
                                           ETPCdEdxSource dedxSource,
                                           Bool_t correctEta,
-                                          Bool_t ratio/*=kFALSE*/) const
+                                          Bool_t correctMultiplicity,
+                                          Bool_t ratio/*=kFALSE*/)const
 {
   //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
@@ -551,12 +645,12 @@ Float_t AliTPCPIDResponse::GetSignalDelta(const AliVTrack* track,
   if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
     return -9999.; //TODO: Better handling!
 
-  const Double_t bethe = GetExpectedSignal(track, species, dEdx, responseFunction, correctEta);
+  const Double_t bethe = GetExpectedSignal(track, species, dEdx, responseFunction, correctEta, correctMultiplicity);
 
   Double_t delta=-9999.;
   if (!ratio) delta=dEdx-bethe;
   else if (bethe>1.e-20) delta=dEdx/bethe;
-
+  
   return delta;
 }
 
@@ -677,12 +771,7 @@ Double_t AliTPCPIDResponse::GetEtaCorrection(const AliVTrack *track, Double_t dE
   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);
+  Double_t tanTheta = GetTrackTanTheta(track); 
   Int_t binX = fhEtaCorr->GetXaxis()->FindBin(tanTheta);
   Int_t binY = fhEtaCorr->GetYaxis()->FindBin(1. / tpcSignal);
   
@@ -718,7 +807,8 @@ Double_t AliTPCPIDResponse::GetEtaCorrection(const AliVTrack *track, AliPID::EPa
   if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
     return 1.; 
   
-  Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE);
+  // For the eta correction, do NOT take the multiplicity corrected value of dEdx
+  Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, kFALSE);
   
   //TODO Alternatively take current track dEdx
   //return GetEtaCorrection(track, dEdx);
@@ -756,7 +846,8 @@ Double_t AliTPCPIDResponse::GetEtaCorrectedTrackdEdx(const AliVTrack *track, Ali
   Double_t etaCorr = 0.;
   
   if (species < AliPID::kUnknown) {
-    Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE);
+    // For the eta correction, do NOT take the multiplicity corrected value of dEdx
+    Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, kFALSE);
     etaCorr = GetEtaCorrection(track, dEdxSplines);
   }
   else {
@@ -770,7 +861,6 @@ Double_t AliTPCPIDResponse::GetEtaCorrectedTrackdEdx(const AliVTrack *track, Ali
 }
 
 
-
 //_________________________________________________________________________
 Double_t AliTPCPIDResponse::GetSigmaPar1(const AliVTrack *track, AliPID::EParticleType species, Double_t dEdx, const TSpline3* responseFunction) const
 {
@@ -791,8 +881,9 @@ Double_t AliTPCPIDResponse::GetSigmaPar1(const AliVTrack *track, AliPID::EPartic
   // 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);
+  
+  // For the eta correction, do NOT take the multiplicity corrected value of dEdx
+  Double_t dEdxExpected = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, kFALSE);
   
   if (dEdxExpected < 1.)
     return 999;
@@ -824,7 +915,7 @@ Double_t AliTPCPIDResponse::GetSigmaPar1(const AliVTrack *track, AliPID::EPartic
 Double_t AliTPCPIDResponse::GetSigmaPar1(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const
 {
   //
-  // Get eta correction for the given track.
+  // Get parameter 1 of sigma parametrisation of TPC dEdx from the histogram for the given track.
   //
   
   if (!fhEtaSigmaPar1)
@@ -865,6 +956,7 @@ Bool_t AliTPCPIDResponse::SetEtaCorrMap(TH2D* hMap)
   return kTRUE;
 }
 
+
 //_________________________________________________________________________
 Bool_t AliTPCPIDResponse::SetSigmaParams(TH2D* hSigmaPar1Map, Double_t sigmaPar0)
 {
@@ -894,6 +986,256 @@ Bool_t AliTPCPIDResponse::SetSigmaParams(TH2D* hSigmaPar1Map, Double_t sigmaPar0
 
 
 //_________________________________________________________________________
+Double_t AliTPCPIDResponse::GetTrackTanTheta(const AliVTrack *track) const
+{
+  // Extract the tanTheta from the information available in the AliVTrack
+  
+  // 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.
+  
+  
+  // Alternatively (in average, the effect is found to be negligable!):
+  // Take local tanTheta from TPC inner wall, if available (currently only for ESDs available)
+  /*const AliExternalTrackParam* innerParam = track->GetInnerParam();
+  if (innerParam) {
+    return innerParam->GetTgl();
+  }*/
+  
+  return TMath::Tan(-track->Theta() + TMath::Pi() / 2.0);
+}
+
+
+//_________________________________________________________________________
+Double_t AliTPCPIDResponse::GetMultiplicityCorrection(const AliVTrack *track, const Double_t dEdxExpected, const Int_t multiplicity) const
+{
+  // Calculate the multiplicity correction factor for this track for the given multiplicity.
+  // The parameter dEdxExpected should take into account the eta correction already!
+  
+  // Multiplicity depends on pure dEdx. Therefore, correction factor depends indirectly on eta
+  // => Use eta corrected value for dEdexExpected.
+  
+  if (dEdxExpected <= 0 || multiplicity <= 0)
+    return 1.0;
+  
+  const Double_t dEdxExpectedInv = 1. / dEdxExpected;
+  Double_t relSlope = fCorrFuncMultiplicity->Eval(dEdxExpectedInv);
+  
+  const Double_t tanTheta = GetTrackTanTheta(track);
+  relSlope += fCorrFuncMultiplicityTanTheta->Eval(tanTheta);
+
+  return (1. + relSlope * multiplicity);
+}
+
+
+//_________________________________________________________________________
+Double_t AliTPCPIDResponse::GetMultiplicityCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const
+{
+  //
+  // Get multiplicity correction for the given track (w.r.t. the mulitplicity of the current event)
+  //
+  
+  //TODO Should return error value, if no eta correction is enabled (needs eta correction). OR: Reset correction in case of no eta correction in PIDresponse
+  
+  // No correction in case of multiplicity <= 0
+  if (fCurrentEventMultiplicity <= 0)
+    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.; 
+  
+  //TODO Does it make sense to use the multiplicity correction WITHOUT eta correction?! Are the fit parameters still valid?
+  // To get the expected signal to determine the multiplicity correction, do NOT ask for the multiplicity corrected value (of course)
+  Double_t dEdxExpected = GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE);
+  
+  return GetMultiplicityCorrection(track, dEdxExpected, fCurrentEventMultiplicity);
+}
+
+
+//_________________________________________________________________________
+Double_t AliTPCPIDResponse::GetMultiplicityCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const
+{
+  //
+  // Get multiplicity 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 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!
+  //
+  
+  //TODO Should return error value, if no eta correction is enabled (needs eta correction). OR: Reset correction in case of no eta correction in PIDresponse
+  
+  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.;
+  
+  
+  // No correction in case of multiplicity <= 0
+  if (fCurrentEventMultiplicity <= 0)
+    return dEdx;
+  
+  Double_t multiplicityCorr = 0.;
+  
+  // TODO Normally, one should use the eta corrected values in BOTH of the following cases. Does it make sense to use the uncorrected dEdx values?
+  if (species < AliPID::kUnknown) {
+    // To get the expected signal to determine the multiplicity correction, do NOT ask for the multiplicity corrected value (of course).
+    // However, one needs the eta corrected value!
+    Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE);
+    multiplicityCorr = GetMultiplicityCorrection(track, dEdxSplines, fCurrentEventMultiplicity);
+  }
+  else {
+    // One needs the eta corrected value to determine the multiplicity correction factor!
+    Double_t etaCorr = GetEtaCorrection(track, dEdx);
+    multiplicityCorr = GetMultiplicityCorrection(track, dEdx * etaCorr, fCurrentEventMultiplicity);
+  }
+    
+  if (multiplicityCorr <= 0)
+    return -1.;
+  
+  return dEdx / multiplicityCorr; 
+}
+
+
+//_________________________________________________________________________
+Double_t AliTPCPIDResponse::GetEtaAndMultiplicityCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species,
+                                                                    ETPCdEdxSource dedxSource) const
+{
+  //
+  // Get multiplicity and 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 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 multiplicityCorr = 0.;
+  Double_t etaCorr = 0.;
+    
+  if (species < AliPID::kUnknown) {
+    // To get the expected signal to determine the multiplicity correction, do NOT ask for the multiplicity corrected value (of course)
+    Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, kFALSE);
+    etaCorr = GetEtaCorrection(track, dEdxSplines);
+    multiplicityCorr = GetMultiplicityCorrection(track, dEdxSplines * etaCorr, fCurrentEventMultiplicity);
+  }
+  else {
+    etaCorr = GetEtaCorrection(track, dEdx);
+    multiplicityCorr = GetMultiplicityCorrection(track, dEdx * etaCorr, fCurrentEventMultiplicity);
+  }
+    
+  if (multiplicityCorr <= 0 || etaCorr <= 0)
+    return -1.;
+  
+  return dEdx / multiplicityCorr / etaCorr; 
+}
+
+
+//_________________________________________________________________________
+Double_t AliTPCPIDResponse::GetMultiplicitySigmaCorrection(const Double_t dEdxExpected, const Int_t multiplicity) const
+{
+  // Calculate the multiplicity sigma correction factor for the corresponding expected dEdx and for the given multiplicity.
+  // The parameter dEdxExpected should take into account the eta correction already!
+  
+  // Multiplicity dependence of sigma depends on the real dEdx at zero multiplicity,
+  // i.e. the eta (only) corrected dEdxExpected value has to be used
+  // since all maps etc. have been created for ~zero multiplicity
+  
+  if (dEdxExpected <= 0 || multiplicity <= 0)
+    return 1.0;
+  
+  Double_t relSigmaSlope = fCorrFuncSigmaMultiplicity->Eval(1. / dEdxExpected);
+  return (1. + relSigmaSlope * multiplicity);
+}
+
+
+//_________________________________________________________________________
+Double_t AliTPCPIDResponse::GetMultiplicitySigmaCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const
+{
+  //
+  // Get multiplicity sigma correction for the given track (w.r.t. the mulitplicity of the current event)
+  //
+  
+  //TODO Should return error value, if no eta correction is enabled (needs eta correction). OR: Reset correction in case of no eta correction in PIDresponse
+  
+  // No correction in case of multiplicity <= 0
+  if (fCurrentEventMultiplicity <= 0)
+    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.; 
+  
+  //TODO Does it make sense to use the multiplicity correction WITHOUT eta correction?! Are the fit parameters still valid?
+  // To get the expected signal to determine the multiplicity correction, do NOT ask for the multiplicity corrected value (of course)
+  Double_t dEdxExpected = GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE);
+  
+  return GetMultiplicitySigmaCorrection(dEdxExpected, fCurrentEventMultiplicity);
+}
+
+
+//_________________________________________________________________________
+void AliTPCPIDResponse::ResetMultiplicityCorrectionFunctions()
+{
+  // Default values: No correction, i.e. overall correction factor should be one
+  
+  // This function should always return 0.0, if multcorr disabled
+  fCorrFuncMultiplicity->SetParameter(0, 0.);
+  fCorrFuncMultiplicity->SetParameter(1, 0.);
+  fCorrFuncMultiplicity->SetParameter(2, 0.);
+  fCorrFuncMultiplicity->SetParameter(3, 0.);
+  fCorrFuncMultiplicity->SetParameter(4, 0.);
+    
+  // This function should always return 0., if multCorr disabled
+  fCorrFuncMultiplicityTanTheta->SetParameter(0, 0.);
+  fCorrFuncMultiplicityTanTheta->SetParameter(1, 0.);
+  fCorrFuncMultiplicityTanTheta->SetParameter(2, 0.);
+  
+  // This function should always return 0.0, if mutlcorr disabled
+  fCorrFuncSigmaMultiplicity->SetParameter(0, 0.);
+  fCorrFuncSigmaMultiplicity->SetParameter(1, 0.);
+  fCorrFuncSigmaMultiplicity->SetParameter(2, 0.);
+  fCorrFuncSigmaMultiplicity->SetParameter(3, 0.);
+}
+
+
+//_________________________________________________________________________
 Bool_t AliTPCPIDResponse::sectorNumbersInOut(Double_t* trackPositionInner,
                                              Double_t* trackPositionOuter,
                                              Float_t& inphi,
index 6dc5a98..9bf2639 100644 (file)
@@ -20,6 +20,7 @@
 #include <TNamed.h>
 #include <TVectorF.h>
 #include <TObjArray.h>
+#include <TF1.h>
 
 #include "AliPID.h"
 #include "AliVTrack.h"
@@ -89,8 +90,10 @@ public:
   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 GetTrackTanTheta(const AliVTrack *track) const;
   
+  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; };
@@ -99,25 +102,56 @@ public:
   
   Double_t GetSigmaPar1(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource = kdEdxDefault) const;
 
+  
+  const TF1* GetMultiplicityCorrectionFunction() const  { return fCorrFuncMultiplicity; };
+  void SetParameterMultiplicityCorrection(Int_t parIndex, Double_t parValue)  
+      { if (fCorrFuncMultiplicity) fCorrFuncMultiplicity->SetParameter(parIndex, parValue); };
+  
+  const TF1* GetMultiplicityCorrectionFunctionTanTheta() const  { return fCorrFuncMultiplicityTanTheta; };
+  void SetParameterMultiplicityCorrectionTanTheta(Int_t parIndex, Double_t parValue)  
+      { if (fCorrFuncMultiplicityTanTheta) fCorrFuncMultiplicityTanTheta->SetParameter(parIndex, parValue); };
+
+  const TF1* GetMultiplicitySigmaCorrectionFunction() const  { return fCorrFuncSigmaMultiplicity; };
+  void SetParameterMultiplicitySigmaCorrection(Int_t parIndex, Double_t parValue)  
+      { if (fCorrFuncSigmaMultiplicity) fCorrFuncSigmaMultiplicity->SetParameter(parIndex, parValue); };
+  
+  void ResetMultiplicityCorrectionFunctions(); 
+  
+  void SetCurrentEventMultiplicity(Int_t value) { fCurrentEventMultiplicity = value;  };
+  Int_t GetCurrentEventMultiplicity() const { return fCurrentEventMultiplicity; };
+
+  Double_t GetMultiplicityCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const;
+  
+  Double_t GetMultiplicitySigmaCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const;
+
+  Double_t GetMultiplicityCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource = kdEdxDefault) const;
+  
+  Double_t GetEtaAndMultiplicityCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species,
+                                                   ETPCdEdxSource dedxSource = kdEdxDefault) const;
   //NEW
   void SetSigma(Float_t res0, Float_t resN2, ETPCgainScenario gainScenario );
   Double_t GetExpectedSignal( const AliVTrack* track,
                               AliPID::EParticleType species,
                               ETPCdEdxSource dedxSource = kdEdxDefault,
-                              Bool_t correctEta = kFALSE) const;//TODO: In future, default kTRUE
+                              Bool_t correctEta = kFALSE,
+                              Bool_t correctMultiplicity = kFALSE) const;
   Double_t GetExpectedSigma( const AliVTrack* track, 
                              AliPID::EParticleType species,
                              ETPCdEdxSource dedxSource = kdEdxDefault,
-                             Bool_t correctEta = kFALSE) const;//TODO: In future, default kTRUE
+                             Bool_t correctEta = kFALSE,
+                             Bool_t correctMultiplicity = kFALSE) const;
   Float_t GetNumberOfSigmas( const AliVTrack* track,
                              AliPID::EParticleType species,
                              ETPCdEdxSource dedxSource = kdEdxDefault,
-                             Bool_t correctEta = kFALSE) const;//TODO: In future, default kTRUE
+                             Bool_t correctEta = kFALSE,
+                             Bool_t correctMultiplicity = kFALSE) const;
   
   Float_t GetSignalDelta( const AliVTrack* track,
                           AliPID::EParticleType species,
                           ETPCdEdxSource dedxSource = kdEdxDefault,
-                          Bool_t correctEta = kFALSE, Bool_t ratio=kFALSE) const;
+                          Bool_t correctEta = kFALSE,
+                          Bool_t correctMultiplicity = kFALSE,
+                          Bool_t ratio = kFALSE) const;
   
   void SetResponseFunction(TObject* o,
                            AliPID::EParticleType type,
@@ -155,7 +189,8 @@ public:
                              AliPID::EParticleType n=AliPID::kKaon) const {
     //
     // Deprecated function (for backward compatibility). Please use 
-    // GetNumberOfSigmas(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource )
+    // GetNumberOfSigmas(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource,
+    // Bool_t correctEta, Bool_t correctMultiplicity)
     // instead!TODO
     //
     
@@ -175,7 +210,8 @@ protected:
                              AliPID::EParticleType species,
                              Double_t dEdx,
                              const TSpline3* responseFunction,
-                             Bool_t correctEta) const; 
+                             Bool_t correctEta,
+                             Bool_t correctMultiplicity) const; 
   
   Double_t GetExpectedSigma(const AliVTrack* track, 
                             AliPID::EParticleType species,
@@ -183,10 +219,15 @@ protected:
                             Double_t dEdx,
                             Int_t nPoints,
                             const TSpline3* responseFunction,
-                            Bool_t correctEta) const;
+                            Bool_t correctEta,
+                            Bool_t correctMultiplicity) const;
                              
   Double_t GetEtaCorrection(const AliVTrack *track, Double_t dEdxSplines) const;
   
+  Double_t GetMultiplicityCorrection(const AliVTrack *track, const Double_t dEdxExpected, const Int_t multiplicity) const;
+  
+  Double_t GetMultiplicitySigmaCorrection(const Double_t dEdxExpected, const Int_t multiplicity) const;
+  
   Double_t GetSigmaPar1(const AliVTrack *track, AliPID::EParticleType species,
                         Double_t dEdx, const TSpline3* responseFunction) const;
   
@@ -221,8 +262,13 @@ private:
   TH2D* fhEtaSigmaPar1; //! Map for parameter 1 of the dEdx sigma parametrisation
   
   Double_t fSigmaPar0; // Parameter 0 of the dEdx sigma parametrisation
+  
+  Int_t fCurrentEventMultiplicity; // Multiplicity of the current event
+  TF1* fCorrFuncMultiplicity; //! Function to correct for the multiplicity dependence of the TPC dEdx
+  TF1* fCorrFuncMultiplicityTanTheta; //! Function to correct the additional tanTheta dependence of the multiplicity dependence of the TPC dEdx
+  TF1* fCorrFuncSigmaMultiplicity; //! Function to correct for the multiplicity dependence of the TPC dEdx resolution
 
-  ClassDef(AliTPCPIDResponse,5)   // TPC PID class
+  ClassDef(AliTPCPIDResponse,6)   // TPC PID class
 };
 
 #endif