/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//-----------------------------------------------------------------
//           Implementation of the TPC PID class
// Very naive one... Should be made better by the detector experts...
//      Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
// 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 <TGraph.h>
#include <TObjArray.h>
#include <TSpline.h>
#include <TBits.h>
#include <TMath.h>
#include <TH2D.h>

#include <AliLog.h>
#include "AliExternalTrackParam.h"
#include "AliVTrack.h"
#include "AliTPCPIDResponse.h"
#include "AliTPCdEdxInfo.h"
#include "TFile.h"
#include "TSpline.h"

ClassImp(AliTPCPIDResponse);


AliTPCPIDResponse *AliTPCPIDResponse::fgInstance =0;

const char* AliTPCPIDResponse::fgkGainScenarioName[fgkNumberOfGainScenarios+1]=
{
  "", //default - no name
  "1", //all high
  "2", //OROC only
  "unknown"
};

//_________________________________________________________________________
AliTPCPIDResponse::AliTPCPIDResponse():
  TNamed(),
  fMIP(50.),
  fRes0(),
  fResN2(),
  fKp1(0.0283086),
  fKp2(2.63394e+01),
  fKp3(5.04114e-11),
  fKp4(2.12543),
  fKp5(4.88663),
  fUseDatabase(kFALSE),
  fResponseFunctions(fgkNumberOfParticleSpecies*fgkNumberOfGainScenarios),
  fVoltageMap(72),
  fLowGainIROCthreshold(-40),
  fBadIROCthreshhold(-70),
  fLowGainOROCthreshold(-40),
  fBadOROCthreshhold(-40),
  fMaxBadLengthFraction(0.5),
  fMagField(0.),
  fhEtaCorr(0x0),
  fhEtaSigmaPar1(0x0),
  fSigmaPar0(0.0),
  fCurrentEventMultiplicity(0),
  fCorrFuncMultiplicity(0x0),
  fCorrFuncMultiplicityTanTheta(0x0),
  fCorrFuncSigmaMultiplicity(0x0),
  fSplineArray()
{
  //
  //  The default constructor
  //
  
  AliLog::SetClassDebugLevel("AliTPCPIDResponse", AliLog::kInfo); 
  
  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();
  fgInstance=this;
}
/*TODO remove?
//_________________________________________________________________________
AliTPCPIDResponse::AliTPCPIDResponse(const Double_t *param):
  TNamed(),
  fMIP(param[0]),
  fRes0(),
  fResN2(),
  fKp1(0.0283086),
  fKp2(2.63394e+01),
  fKp3(5.04114e-11),
  fKp4(2.12543),
  fKp5(4.88663),
  fUseDatabase(kFALSE),
  fResponseFunctions(fgkNumberOfParticleSpecies*fgkNumberOfGainScenarios),
  fVoltageMap(72),
  fLowGainIROCthreshold(-40),
  fBadIROCthreshhold(-70),
  fLowGainOROCthreshold(-40),
  fBadOROCthreshhold(-40),
  fMaxBadLengthFraction(0.5),
  fMagField(0.),
  fhEtaCorr(0x0),
  fhEtaSigmaPar1(0x0),
  fSigmaPar0(0.0)
{
  //
  //  The main constructor
  //
  for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=param[1];fResN2[i]=param[2];}
}
*/

//_________________________________________________________________________
AliTPCPIDResponse::~AliTPCPIDResponse()
{
  //
  // Destructor
  //
  
  delete fhEtaCorr;
  fhEtaCorr = 0x0;
  
  delete fhEtaSigmaPar1;
  fhEtaSigmaPar1 = 0x0;
  
  delete fCorrFuncMultiplicity;
  fCorrFuncMultiplicity = 0x0;
  
  delete fCorrFuncMultiplicityTanTheta;
  fCorrFuncMultiplicityTanTheta = 0x0;
  
  delete fCorrFuncSigmaMultiplicity;
  fCorrFuncSigmaMultiplicity = 0x0;
  if (fgInstance==this) fgInstance=0;
}


//_________________________________________________________________________
AliTPCPIDResponse::AliTPCPIDResponse(const AliTPCPIDResponse& that):
  TNamed(that),
  fMIP(that.fMIP),
  fRes0(),
  fResN2(),
  fKp1(that.fKp1),
  fKp2(that.fKp2),
  fKp3(that.fKp3),
  fKp4(that.fKp4),
  fKp5(that.fKp5),
  fUseDatabase(that.fUseDatabase),
  fResponseFunctions(that.fResponseFunctions),
  fVoltageMap(that.fVoltageMap),
  fLowGainIROCthreshold(that.fLowGainIROCthreshold),
  fBadIROCthreshhold(that.fBadIROCthreshhold),
  fLowGainOROCthreshold(that.fLowGainOROCthreshold),
  fBadOROCthreshhold(that.fBadOROCthreshhold),
  fMaxBadLengthFraction(that.fMaxBadLengthFraction),
  fMagField(that.fMagField),
  fhEtaCorr(0x0),
  fhEtaSigmaPar1(0x0),
  fSigmaPar0(that.fSigmaPar0),
  fCurrentEventMultiplicity(that.fCurrentEventMultiplicity),
  fCorrFuncMultiplicity(0x0),
  fCorrFuncMultiplicityTanTheta(0x0),
  fCorrFuncSigmaMultiplicity(0x0),
  fSplineArray()
{
  //copy ctor
  for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=that.fRes0[i];fResN2[i]=that.fResN2[i];}
 
  // Copy eta maps
  if (that.fhEtaCorr) {
    fhEtaCorr = new TH2D(*(that.fhEtaCorr));
    fhEtaCorr->SetDirectory(0);
  }
  
  if (that.fhEtaSigmaPar1) {
    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));
  }
}

//_________________________________________________________________________
AliTPCPIDResponse& AliTPCPIDResponse::operator=(const AliTPCPIDResponse& that)
{
  //assignment
  if (&that==this) return *this;
  TNamed::operator=(that);
  fMIP=that.fMIP;
  fKp1=that.fKp1;
  fKp2=that.fKp2;
  fKp3=that.fKp3;
  fKp4=that.fKp4;
  fKp5=that.fKp5;
  fUseDatabase=that.fUseDatabase;
  fResponseFunctions=that.fResponseFunctions;
  fVoltageMap=that.fVoltageMap;
  fLowGainIROCthreshold=that.fLowGainIROCthreshold;
  fBadIROCthreshhold=that.fBadIROCthreshhold;
  fLowGainOROCthreshold=that.fLowGainOROCthreshold;
  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;
  fhEtaCorr=0x0;
  if (that.fhEtaCorr) {
    fhEtaCorr = new TH2D(*(that.fhEtaCorr));
    fhEtaCorr->SetDirectory(0);
  }
  
  delete fhEtaSigmaPar1;
  fhEtaSigmaPar1=0x0;
  if (that.fhEtaSigmaPar1) {
    fhEtaSigmaPar1 = new TH2D(*(that.fhEtaSigmaPar1));
    fhEtaSigmaPar1->SetDirectory(0);
  }
  
  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;
}

//_________________________________________________________________________
Double_t AliTPCPIDResponse::Bethe(Double_t betaGamma) const {
  //
  // This is the Bethe-Bloch function normalised to 1 at the minimum
  // WARNING
  // Simulated and reconstructed Bethe-Bloch differs
  //           Simulated  curve is the dNprim/dx
  //           Reconstructed is proportianal dNtot/dx
  // Temporary fix for production -  Simple linear correction function
  // Future    2 Bethe Bloch formulas needed
  //           1. for simulation
  //           2. for reconstructed PID
  //
  
//   const Float_t kmeanCorrection =0.1;
  Double_t bb=
    AliExternalTrackParam::BetheBlochAleph(betaGamma,fKp1,fKp2,fKp3,fKp4,fKp5);
  return bb*fMIP;
}

//_________________________________________________________________________
void AliTPCPIDResponse::SetBetheBlochParameters(Double_t kp1,
                             Double_t kp2,
                             Double_t kp3,
                             Double_t kp4,
                             Double_t kp5) {
  //
  // Set the parameters of the ALEPH Bethe-Bloch formula
  //
  fKp1=kp1;
  fKp2=kp2;
  fKp3=kp3;
  fKp4=kp4;
  fKp5=kp5;
}

//_________________________________________________________________________
void AliTPCPIDResponse::SetSigma(Float_t res0, Float_t resN2) {
  //
  // Set the relative resolution  sigma_rel = res0 * sqrt(1+resN2/npoint)
  //
  for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=res0;fResN2[i]=resN2;}
}

//_________________________________________________________________________
Double_t AliTPCPIDResponse::GetExpectedSignal(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, Bool_t correctMultiplicity);
  // instead!
  //
  //
  // 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. 
  // This can be improved. By taking into account the number of
  // 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) * chargeFactor;
  //
  const TSpline3 * responseFunction = (TSpline3 *) fResponseFunctions.UncheckedAt(n);

  if (!responseFunction) return Bethe(mom/mass) * chargeFactor;
  
  return fMIP*responseFunction->Eval(mom/mass)*chargeFactor;

}

//_________________________________________________________________________
Double_t AliTPCPIDResponse::GetExpectedSigma(Float_t mom, 
                                             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 
  //  
  
  if (nPoints != 0) 
    return GetExpectedSignal(mom,n)*fRes0[0]*sqrt(1. + fResN2[0]/nPoints);
  else
    return GetExpectedSignal(mom,n)*fRes0[0];
}

////////////////////////////////////////////////////NEW//////////////////////////////

//_________________________________________________________________________
void AliTPCPIDResponse::SetSigma(Float_t res0, Float_t resN2, ETPCgainScenario gainScenario) {
  //
  // Set the relative resolution  sigma_rel = res0 * sqrt(1+resN2/npoint)
  //
  if ((Int_t)gainScenario>=(Int_t)fgkNumberOfGainScenarios) return; //TODO: better handling!
  fRes0[gainScenario]=res0;
  fResN2[gainScenario]=resN2;
}


//_________________________________________________________________________
Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track,
                                              AliPID::EParticleType species,
                                              Double_t /*dEdx*/,
                                              const TSpline3* responseFunction,
                                              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
  // 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.  
  //
  
  Double_t mom=track->GetTPCmomentum();
  Double_t mass=AliPID::ParticleMassZ(species);
  
  //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;
  
  Double_t dEdxSplines = fMIP*responseFunction->Eval(mom/mass) * chargeFactor;
    
  if (!correctEta && !correctMultiplicity)
    return dEdxSplines;
  
  Double_t corrFactorEta = 1.0;
  Double_t corrFactorMultiplicity = 1.0;
  
  if (correctEta) {
    corrFactorEta = GetEtaCorrectionFast(track, dEdxSplines);
    //TODO Alternatively take current track dEdx
    //corrFactorEta = GetEtaCorrectionFast(track, dEdx);
  }
  
  if (correctMultiplicity)
    corrFactorMultiplicity = GetMultiplicityCorrectionFast(track, dEdxSplines * corrFactorEta, fCurrentEventMultiplicity);

  return dEdxSplines * corrFactorEta * corrFactorMultiplicity;
}


//_________________________________________________________________________
Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track,
                                              AliPID::EParticleType species,
                                              ETPCdEdxSource dedxSource,
                                              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
  // 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.  
  //
  
  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);
  
    return Bethe(track->GetTPCmomentum() / AliPID::ParticleMassZ(species)) * chargeFactor;
  }
  
  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, correctMultiplicity);
}
  
//_________________________________________________________________________
TSpline3* AliTPCPIDResponse::GetResponseFunction( AliPID::EParticleType type,
                                                  AliTPCPIDResponse::ETPCgainScenario gainScenario ) const
{
  //get response function
  return dynamic_cast<TSpline3*>(fResponseFunctions.At(ResponseFunctionIndex(type,gainScenario)));
}

//_________________________________________________________________________
TSpline3* AliTPCPIDResponse::GetResponseFunction( const AliVTrack* track,
                               AliPID::EParticleType species,
                               ETPCdEdxSource dedxSource) const 
{
  //the splines are stored in an array, different scenarios

  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;
}

//_________________________________________________________________________
void AliTPCPIDResponse::ResetSplines()
{
  //reset the array with splines
  for (Int_t i=0;i<fResponseFunctions.GetEntriesFast();i++)
  {
    fResponseFunctions.AddAt(NULL,i);
  }
}
//_________________________________________________________________________
Int_t AliTPCPIDResponse::ResponseFunctionIndex( AliPID::EParticleType species,
                                                ETPCgainScenario gainScenario ) const
{
  //get the index in fResponseFunctions given type and scenario
  return Int_t(species)+Int_t(gainScenario)*fgkNumberOfParticleSpecies;
}

//_________________________________________________________________________
void AliTPCPIDResponse::SetResponseFunction( TObject* o,
                                             AliPID::EParticleType species,
                                             ETPCgainScenario gainScenario )
{
  fResponseFunctions.AddAtAndExpand(o,ResponseFunctionIndex(species,gainScenario));
}


//_________________________________________________________________________
Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track, 
                                             AliPID::EParticleType species,
                                             ETPCgainScenario gainScenario,
                                             Double_t dEdx,
                                             Int_t nPoints,
                                             const TSpline3* responseFunction,
                                             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,
  // for the specified particle type 
  //
  
  if (!responseFunction)
    return 999;
    
  //TODO Check whether it makes sense to set correctMultiplicity to kTRUE while correctEta might be kFALSE
  
  // If eta correction (=> new sigma parametrisation) is requested, but no sigma map is available, print error message
  if (correctEta && !fhEtaSigmaPar1) {
    AliError("New sigma parametrisation requested, but sigma map not initialised (usually via AliPIDResponse). Old sigma parametrisation will be used!");
  }
  
  // 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, correctMultiplicity) *
               fRes0[gainScenario] * sqrt(1. + fResN2[gainScenario]/nPoints);
    else
      return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, correctMultiplicity)*fRes0[gainScenario];
  }
    
  if (nPoints > 0) {
    // Use eta correction (+ eta-dependent sigma)
    Double_t sigmaPar1 = GetSigmaPar1Fast(track, species, dEdx, responseFunction);
    
    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 = GetMultiplicityCorrectionFast(track, dEdxExpectedEtaCorrected, fCurrentEventMultiplicity);
      Double_t multiplicitySigmaCorrFactor = GetMultiplicitySigmaCorrectionFast(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!
    return 999;
  }
}


//_________________________________________________________________________
Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track, 
                                             AliPID::EParticleType species,
                                             ETPCdEdxSource dedxSource,
                                             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 
  // and dedx scenario
  //
  
  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, correctMultiplicity);
}


//_________________________________________________________________________
Float_t AliTPCPIDResponse::GetNumberOfSigmas(const AliVTrack* track, 
                             AliPID::EParticleType species,
                             ETPCdEdxSource dedxSource,
                             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
  //inside the TPC
  
  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, 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;
  else
    return (dEdx-bethe)/sigma;
}

//_________________________________________________________________________
Float_t AliTPCPIDResponse::GetSignalDelta(const AliVTrack* track,
                                          AliPID::EParticleType species,
                                          ETPCdEdxSource dedxSource,
                                          Bool_t correctEta,
                                          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
  //inside the TPC

  Double_t dEdx = -1;
  Int_t nPoints = -1;
  ETPCgainScenario gainScenario = kGainScenarioInvalid;
  TSpline3* responseFunction = 0x0;

  if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
    return -9999.; //TODO: Better handling!

  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;
}

//_________________________________________________________________________
Bool_t AliTPCPIDResponse::ResponseFunctiondEdxN( const AliVTrack* track, 
                                                 AliPID::EParticleType species,
                                                 ETPCdEdxSource dedxSource,
                                                 Double_t& dEdx,
                                                 Int_t& nPoints,
                                                 ETPCgainScenario& gainScenario,
                                                 TSpline3** responseFunction) const 
{
  // Calculates the right parameters for PID
  //   dEdx parametrization for the proper gain scenario, dEdx 
  //   and NPoints used for dEdx
  // based on the track geometry (which chambers it crosses) for the specified particle type 
  // 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
    
    // GetTPCsignalTunedOnData will be non-positive, if it has not been set (i.e. in case of MC NOT tuned to data).
    // If this is the case, just take the normal signal
    dEdx = track->GetTPCsignalTunedOnData();
    if (dEdx <= 0) {
      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;
  }
  
  //TODO Proper handle of tuneMConData for other dEdx sources
  
  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
  const AliTPCdEdxInfo* dEdxInfo = track->GetTPCdEdxInfo();
  
  if (!dEdxInfo && dedxSource!=kdEdxDefault)  //in one case its ok if we dont have the info
  {
    AliError("AliTPCdEdxInfo not available");
    return kFALSE;
  }

  if (dEdxInfo) dEdxInfo->GetTPCSignalRegionInfo(signal,ncl,nrows);

  //check if we cross a bad OROC in which case we reject
  EChamberStatus trackOROCStatus = TrackStatus(track,2);
  if (trackOROCStatus==kChamberOff || trackOROCStatus==kChamberLowGain)
  {
    return kFALSE;
  }

  switch (dedxSource)
  {
    case kdEdxOROC:
      {
        if (trackOROCStatus==kChamberInvalid) return kFALSE; //never reached OROC
        dEdx = signal[3];
        nPoints = ncl[2]+ncl[1];
        gainScenario = kOROChigh;
        break;
      }
    case kdEdxHybrid:
      {
        //if we cross a bad IROC we use OROC dedx, if we dont we use combined
        EChamberStatus status = TrackStatus(track,1);
        if (status!=kChamberHighGain)
        {
          dEdx = signal[3];
          nPoints = ncl[2]+ncl[1];
          gainScenario = kOROChigh;
        }
        else
        {
          dEdx = track->GetTPCsignal();
          nPoints = track->GetTPCsignalN();
          gainScenario = kALLhigh;
        }
        break;
      }
    default:
      {
         dEdx = 0.;
         nPoints = 0;
         gainScenario = kGainScenarioInvalid;
         return kFALSE;
      }
  }
  TObject* obj = fResponseFunctions.UncheckedAt(ResponseFunctionIndex(species,gainScenario));
  *responseFunction = dynamic_cast<TSpline3*>(obj); //TODO:maybe static cast?
  
  return kTRUE;
}


//_________________________________________________________________________
Double_t AliTPCPIDResponse::GetEtaCorrectionFast(const AliVTrack *track, Double_t dEdxSplines) const
{
  // NOTE: For expert use only -> Non-experts are advised to use the function without the "Fast" suffix or stick to AliPIDResponse directly.
  //
  // Get eta correction for the given parameters.
  //
  
  if (!fhEtaCorr) {
    // Calling this function means to request eta correction in some way. Print error message, if no map is available!
    AliError("Eta correction requested, but map not initialised (usually via AliPIDResponse). Returning eta correction factor 1!");
    return 1.;
  }
  
  Double_t tpcSignal = dEdxSplines;
  
  if (tpcSignal < 1.)
    return 1.;
  
  Double_t tanTheta = GetTrackTanTheta(track); 
  Int_t binX = fhEtaCorr->GetXaxis()->FindFixBin(tanTheta);
  Int_t binY = fhEtaCorr->GetYaxis()->FindFixBin(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.; 
  
  // 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 GetEtaCorrectionFast(track, dEdx);
  
  return GetEtaCorrectionFast(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) {
    // For the eta correction, do NOT take the multiplicity corrected value of dEdx
    Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, kFALSE);
    etaCorr = GetEtaCorrectionFast(track, dEdxSplines);
  }
  else {
    etaCorr = GetEtaCorrectionFast(track, dEdx);
  }
    
  if (etaCorr <= 0)
    return -1.;
  
  return dEdx / etaCorr; 
}


//_________________________________________________________________________
Double_t AliTPCPIDResponse::GetSigmaPar1Fast(const AliVTrack *track, AliPID::EParticleType species, Double_t dEdx,
                                             const TSpline3* responseFunction) const
{
  // NOTE: For expert use only -> Non-experts are advised to use the function without the "Fast" suffix or stick to AliPIDResponse directly.
  //
  // Get parameter 1 of sigma parametrisation of TPC dEdx from the histogram for the given track.
  //
  
  if (!fhEtaSigmaPar1) {
    // Calling this function means to request new sigma parametrisation in some way. Print error message, if no map is available!
    AliError("New sigma parametrisation requested, but sigma map not initialised (usually via AliPIDResponse). Returning error value for sigma parameter1 = 999!");
    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
  
  // 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;
  
  Double_t tanTheta = GetTrackTanTheta(track);
  Int_t binX = fhEtaSigmaPar1->GetXaxis()->FindFixBin(tanTheta);
  Int_t binY = fhEtaSigmaPar1->GetYaxis()->FindFixBin(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 parameter 1 of sigma parametrisation of TPC dEdx from the histogram 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 GetSigmaPar1Fast(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;
}


//_________________________________________________________________________
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.
  
  /*
  const AliExternalTrackParam* innerParam = track->GetInnerParam();
  Double_t tanTheta = 0;
  if (innerParam) 
    tanTheta = innerParam->GetTgl();
  else
    tanTheta = TMath::Tan(-track->Theta() + TMath::Pi() / 2.0);
  
  // Constant in formula for B in kGauss (factor 0.1 to convert B from Tesla to kGauss),
  // pT in GeV/c (factor c*1e-9 to arrive at GeV/c) and curvature in 1/cm (factor 0.01 to get from m to cm)
  const Double_t constant = TMath::C()* 1e-9 * 0.1 * 0.01; 
  const Double_t curvature = fMagField * constant / track->Pt(); // in 1./cm
  
  Double_t averageddzdr = 0.;
  Int_t nParts = 0;

  for (Double_t r = 85; r < 245; r++) {
    Double_t sinPhiLocal = TMath::Abs(r*curvature*0.5);
    
    // Cut on |sin(phi)| as during reco
    if (TMath::Abs(sinPhiLocal) <= 0.95) {
      const Double_t phiLocal = TMath::ASin(sinPhiLocal);
      const Double_t tanPhiLocal = TMath::Tan(phiLocal);
      
      averageddzdr += tanTheta * TMath::Sqrt(1. + tanPhiLocal * tanPhiLocal); 
      nParts++;
    }
  }
  
  if (nParts > 0)
    averageddzdr /= nParts; 
  else {
    AliError("Problems during determination of dz/dr. Returning pure tanTheta as best estimate!");
    return tanTheta;
  }
  
  //printf("pT: %f\nFactor/magField(kGs)/curvature^-1: %f / %f /%f\ntanThetaGlobalFromTheta/tanTheta/Averageddzdr: %f / %f / %f\n\n",
  //          track->Pt(), constant, fMagField, 1./curvature, TMath::Tan(-track->Theta() + TMath::Pi() / 2.0), tanTheta, averageddzdr);
  
  return averageddzdr;
  */
  
  
  // 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::GetMultiplicityCorrectionFast(const AliVTrack *track, Double_t dEdxExpected, Int_t multiplicity) const
{
  // NOTE: For expert use only -> Non-experts are advised to use the function without the "Fast" suffix or stick to AliPIDResponse directly.
  //
  // 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 GetMultiplicityCorrectionFast(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 = GetMultiplicityCorrectionFast(track, dEdxSplines, fCurrentEventMultiplicity);
  }
  else {
    // One needs the eta corrected value to determine the multiplicity correction factor!
    Double_t etaCorr = GetEtaCorrectionFast(track, dEdx);
    multiplicityCorr = GetMultiplicityCorrectionFast(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 = GetEtaCorrectionFast(track, dEdxSplines);
    multiplicityCorr = GetMultiplicityCorrectionFast(track, dEdxSplines * etaCorr, fCurrentEventMultiplicity);
  }
  else {
    etaCorr = GetEtaCorrectionFast(track, dEdx);
    multiplicityCorr = GetMultiplicityCorrectionFast(track, dEdx * etaCorr, fCurrentEventMultiplicity);
  }
    
  if (multiplicityCorr <= 0 || etaCorr <= 0)
    return -1.;
  
  return dEdx / multiplicityCorr / etaCorr; 
}


//_________________________________________________________________________
Double_t AliTPCPIDResponse::GetMultiplicitySigmaCorrectionFast(Double_t dEdxExpected, Int_t multiplicity) const
{
  // NOTE: For expert use only -> Non-experts are advised to use the function without the "Fast" suffix or stick to AliPIDResponse directly.
  //
  // 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 GetMultiplicitySigmaCorrectionFast(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,
                                             Float_t& outphi,
                                             Int_t& in,
                                             Int_t& out ) const
{
  //calculate the sector numbers (equivalent to IROC chamber numbers) a track crosses
  //for OROC chamber numbers add 36
  //returned angles are between (0,2pi)

  inphi = TMath::ATan2(trackPositionInner[1],trackPositionInner[0]);
  outphi = TMath::ATan2(trackPositionOuter[1], trackPositionOuter[0]);

  if (inphi<0) {inphi+=TMath::TwoPi();} //because ATan2 gives angles between -pi,pi
  if (outphi<0) {outphi+=TMath::TwoPi();} //because ATan2 gives angles between -pi,pi

  in = sectorNumber(inphi);
  out = sectorNumber(outphi);

  //for the C side (positive z) offset by 18
  if (trackPositionInner[2]>0.0)
  {
    in+=18;
    out+=18;
  }
  return kTRUE;
}


//_____________________________________________________________________________
Int_t AliTPCPIDResponse::sectorNumber(Double_t phi) const
{
  //calculate sector number
  const Float_t width=TMath::TwoPi()/18.;
  return TMath::Floor(phi/width);
}


//_____________________________________________________________________________
void AliTPCPIDResponse::Print(Option_t* /*option*/) const
{
  //Print info
  fResponseFunctions.Print();
}


//_____________________________________________________________________________
AliTPCPIDResponse::EChamberStatus AliTPCPIDResponse::TrackStatus(const AliVTrack* track, Int_t layer) const
{
  //status of the track: if it crosses any bad chambers return kChamberOff;
  //IROC:layer=1, OROC:layer=2
  if (layer<1 || layer>2) layer=1;
  Int_t in=0;
  Int_t out=0;
  Float_t inphi=0.;
  Float_t outphi=0.;
  Float_t innerRadius = (layer==1)?83.0:133.7;
  Float_t outerRadius = (layer==1)?133.5:247.7;

  /////////////////////////////////////////////////////////////////////////////
  //find out where track enters and leaves the layer.
  //
  Double_t trackPositionInner[3];
  Double_t trackPositionOuter[3];
 
  //if there is no inner param this could mean we're using the AOD track,
  //we still can extrapolate from the vertex - so use those params.
  const AliExternalTrackParam* ip = track->GetInnerParam();
  if (ip) track=dynamic_cast<const AliVTrack*>(ip);

  Bool_t trackAtInner = track->GetXYZAt(innerRadius, fMagField, trackPositionInner);
  Bool_t trackAtOuter = track->GetXYZAt(outerRadius, fMagField, trackPositionOuter);

  if (!trackAtInner)
  {
    //if we dont even enter inner radius we do nothing and return invalid
    inphi=0.0;
    outphi=0.0;
    in=0;
    out=0;
    return kChamberInvalid;
  }

  if (!trackAtOuter)
  {
    //if we don't reach the outer radius check that the apex is indeed within the outer radius and use apex position
    Bool_t haveApex = TrackApex(track, fMagField, trackPositionOuter);
    Float_t apexRadius = TMath::Sqrt(trackPositionOuter[0]*trackPositionOuter[0]+trackPositionOuter[1]*trackPositionOuter[1]);
    if ( haveApex && apexRadius<=outerRadius && apexRadius>innerRadius)
    {
      //printf("pt: %.2f, apexRadius: %.2f(%s), x: %.2f, y: %.2f\n",track->Pt(),apexRadius,(haveApex)?"OK":"BAD",trackPositionOuter[0],trackPositionOuter[1]);
    }
    else
    {
      inphi=0.0;
      outphi=0.0;
      in=0;
      out=0;
      return kChamberInvalid;
    }
  }


  if (!sectorNumbersInOut(trackPositionInner,
                          trackPositionOuter,
                          inphi,
                          outphi,
                          in,
                          out)) return kChamberInvalid;

  /////////////////////////////////////////////////////////////////////////////
  //now we have the location of the track we can check
  //if it is in a good/bad chamber
  //
  Bool_t sideA = kTRUE;
 
  if (((in/18)%2==1) && ((out/18)%2==1)) sideA=kFALSE;

  in=in%18;
  out=out%18;

  if (TMath::Abs(in-out)>9)
  {
    if (TMath::Max(in,out)==out)
    {
      Int_t tmp=out;
      out = in;
      in = tmp;
      Float_t tmpphi=outphi;
      outphi=inphi;
      inphi=tmpphi;
    }
    in-=18;
    inphi-=TMath::TwoPi();
  }
  else
  {
    if (TMath::Max(in,out)==in)
    {
      Int_t tmp=out;
      out=in;
      in=tmp;
      Float_t tmpphi=outphi;
      outphi=inphi;
      inphi=tmpphi;
    }
  }

  Float_t trackLengthInBad = 0.;
  Float_t trackLengthInLowGain = 0.;
  Float_t trackLengthTotal = TMath::Abs(outphi-inphi);
  Float_t lengthFractionInBadSectors = 0.;

  const Float_t sectorWidth = TMath::TwoPi()/18.;  
 
  for (Int_t i=in; i<=out; i++)
  {
    int j=i;
    if (i<0) j+=18;    //correct for the negative values
    if (!sideA) j+=18; //move to the correct side
   
    Float_t deltaPhi = 0.;
    Float_t phiEdge=sectorWidth*i;
    if (inphi>phiEdge) {deltaPhi=phiEdge+sectorWidth-inphi;}
    else if ((outphi>=phiEdge) && (outphi<(phiEdge+sectorWidth))) {deltaPhi=outphi-phiEdge;}
    else {deltaPhi=sectorWidth;}
   
    Float_t v = fVoltageMap[(layer==1)?(j):(j+36)];
    if (v<=fBadIROCthreshhold)
    {
      trackLengthInBad+=deltaPhi;
      lengthFractionInBadSectors=1.;
    }
    if (v<=fLowGainIROCthreshold && v>fBadIROCthreshhold)
    {
      trackLengthInLowGain+=deltaPhi;
      lengthFractionInBadSectors=1.;
    }
  }

  //for now low gain and bad (off) chambers are treated equally
  if (trackLengthTotal>0)
    lengthFractionInBadSectors = (trackLengthInLowGain+trackLengthInBad)/trackLengthTotal;

  //printf("### side: %s, pt: %.2f, pz: %.2f, in: %i, out: %i, phiIN: %.2f, phiOUT: %.2f, rIN: %.2f, rOUT: %.2f\n",(sideA)?"A":"C",track->Pt(),track->Pz(),in,out,inphi,outphi,innerRadius,outerRadius);
 
  if (lengthFractionInBadSectors>fMaxBadLengthFraction)
  {
    //printf("%%%%%%%% %s kChamberLowGain\n",(layer==1)?"IROC":"OROC");
    return kChamberLowGain;
  }
 
  return kChamberHighGain;
}


//_____________________________________________________________________________
Float_t AliTPCPIDResponse::MaxClusterRadius(const AliVTrack* track) const
{
  //return the radius of the outermost padrow containing a cluster in TPC
  //for the track
  const TBits* clusterMap=track->GetTPCClusterMapPtr();
  if (!clusterMap) return 0.;

  //from AliTPCParam, radius of first IROC row
  const Float_t rfirstIROC = 8.52249984741210938e+01;
  const Float_t padrowHeightIROC = 0.75;
  const Float_t rfirstOROC0 = 1.35100006103515625e+02;
  const Float_t padrowHeightOROC0 = 1.0;
  const Float_t rfirstOROC1 = 1.99350006103515625e+02;
  const Float_t padrowHeightOROC1 = 1.5;

  Int_t maxPadRow=160;
  while ((--maxPadRow)>0 && !clusterMap->TestBitNumber(maxPadRow)){}
  if (maxPadRow>126) return rfirstOROC1+(maxPadRow-126-1)*padrowHeightOROC1;
  if (maxPadRow>62) return rfirstOROC0+(maxPadRow-62-1)*padrowHeightOROC0;
  if (maxPadRow>0) return rfirstIROC+(maxPadRow-1)*padrowHeightIROC;
  return 0.0;
}


//_____________________________________________________________________________
Bool_t AliTPCPIDResponse::TrackApex(const AliVTrack* track, Float_t magField, Double_t position[3]) const
{
  //calculate the coordinates of the apex of the track
  Double_t x[3];
  track->GetXYZ(x);
  Double_t p[3];
  track->GetPxPyPz(p);
  Double_t r = 1./track->OneOverPt()/0.0299792458/magField; //signed - will determine the direction of b
  //printf("b: %.2f, x:%.2f, y:%.2f, pt: %.2f, px:%.2f, py%.2f, r: %.2f\n",magField, x[0],x[1],track->Pt(), p[0],p[1],r);
  //find orthogonal vector (Gram-Schmidt)
  Double_t alpha = (p[0]*x[0] + p[1]*x[1])/(p[0]*p[0] + p[1]*p[1]);
  Double_t b[2];
  b[0]=x[0]-alpha*p[0];
  b[1]=x[1]-alpha*p[1];
 
  Double_t norm = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]);
  if (TMath::AreEqualAbs(norm,0.0,1e-10)) return kFALSE;
  b[0]/=norm;
  b[1]/=norm;
  b[0]*=r;
  b[1]*=r;
  b[0]+=x[0];
  b[1]+=x[1];
  //printf("center: x:%.2f, y:%.2f\n",b[0],b[1]);
 
  norm = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]);
  if (TMath::AreEqualAbs(norm,0.0,1e-10)) return kFALSE;
 
  position[0]=b[0]+b[0]*TMath::Abs(r)/norm;
  position[1]=b[1]+b[1]*TMath::Abs(r)/norm;
  position[2]=0.;
  return kTRUE;
}

Double_t AliTPCPIDResponse::EvaldEdxSpline(Double_t bg,Int_t entry){
  //
  // Evaluate the dEdx response for given entry
  //
  TSpline * spline = (TSpline*)fSplineArray.At(entry);
  if (spline) return spline->Eval(bg);
  return 0;
}


Bool_t   AliTPCPIDResponse::RegisterSpline(const char * name, Int_t index){
  //
  // register spline to be used for drawing comparisons
  // 
  TFile * fTPCBB = TFile::Open("$ALICE_ROOT/OADB/COMMON/PID/data/TPCPIDResponse.root");
  TObjArray  *arrayTPCPID= (TObjArray*)  fTPCBB->Get("TPCPIDResponse");
  if (fSplineArray.GetEntriesFast()<index) fSplineArray.Expand(index*2);
  TSpline3 *spline=0;
  if (arrayTPCPID){
    spline = (TSpline3*)arrayTPCPID->FindObject(name);
    if (spline) fSplineArray.AddAt(spline->Clone(),index);    
  }
  delete arrayTPCPID;
  delete fTPCBB;
  return (spline!=0);
}