- introduction of gain scenarios (e.g. OROC only - for homogeneous gain in 11h)
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 3 Sep 2012 13:57:07 +0000 (13:57 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 3 Sep 2012 13:57:07 +0000 (13:57 +0000)
- splines for gain scenarios in 11h
- splines for light nuclei in 10h
Jens Wiechula <Jens.Wiechula@cern.ch>
mikolaj <m.krzewicki@gsi.de>

ANALYSIS/AliAnalysisTaskPIDqa.cxx
OADB/COMMON/PID/data/TPCPIDResponse.root
OADB/COMMON/PID/data/TPCvoltageSettings.root [new file with mode: 0644]
STEER/STEERBase/AliPIDResponse.cxx
STEER/STEERBase/AliPIDResponse.h
STEER/STEERBase/AliTPCPIDResponse.cxx
STEER/STEERBase/AliTPCPIDResponse.h

index 54b8ff1..5e1d623 100644 (file)
@@ -428,15 +428,35 @@ void AliAnalysisTaskPIDqa::FillTPCqa()
     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;\r
     \r
     Double_t mom=track->GetTPCmomentum();\r
-    \r
+    // the default scenario\r
     for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){\r
       TH2 *h=(TH2*)fListQAtpc->At(ispecie);\r
       if (!h) continue;\r
       Double_t nSigma=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie);\r
       h->Fill(mom,nSigma);\r
     }\r
+    // the "hybrid" scenario\r
+    for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){\r
+      TH2 *h=(TH2*)fListQAtpc->At(ispecie+AliPID::kSPECIESC);\r
+      if (!h) continue;\r
+      Double_t nSigma=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie, AliTPCPIDResponse::kdEdxHybrid);\r
+      h->Fill(mom,nSigma);\r
+    }\r
+    \r
+    // the "OROC" scenario\r
+    for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){\r
+      TH2 *h=(TH2*)fListQAtpc->At(ispecie+2*AliPID::kSPECIESC);\r
+      if (!h) continue;\r
+      Double_t nSigma=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie, AliTPCPIDResponse::kdEdxOROC);\r
+      //TSpline3* spline = fPIDResponse->GetTPCResponse().GetCurrentResponseFunction();\r
+      //std::cout<<ispecie<<" "<<nSigma<<" phi:"<<track->Phi()<<". "<<std::endl;\r
+      //if (spline) {cout<<spline->GetName()<<endl;}\r
+      //else {cout<<"NULL spline"<<endl;}\r
+      h->Fill(mom,nSigma);\r
+    }\r
     \r
-    TH2 *h=(TH2*)fListQAtpc->At(AliPID::kSPECIESC);\r
+    TH2 *h=(TH2*)fListQAtpc->At(3*AliPID::kSPECIESC);\r
+\r
     if (h) {\r
       Double_t sig=track->GetTPCsignal();\r
       h->Fill(mom,sig);\r
@@ -989,13 +1009,32 @@ void AliAnalysisTaskPIDqa::SetupTPCqa()
                               200,-10,10);\r
     fListQAtpc->Add(hNsigmaP);\r
   }\r
+\r
+  // the "hybrid" scenario\r
+  for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){\r
+    TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TPC_%s_Hybrid",AliPID::ParticleName(ispecie)),\r
+                              Form("TPC n#sigma %s vs. p (Hybrid gain scenario);p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),\r
+                              vX->GetNrows()-1,vX->GetMatrixArray(),\r
+                              200,-10,10);\r
+    fListQAtpc->Add(hNsigmaP);\r
+  }\r
+   \r
+  // the "OROC high" scenario\r
+  for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){\r
+    TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TPC_%s_OROChigh",AliPID::ParticleName(ispecie)),\r
+                              Form("TPC n#sigma %s vs. p (OROChigh gain scenario);p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),\r
+                              vX->GetNrows()-1,vX->GetMatrixArray(),\r
+                              200,-10,10);\r
+    fListQAtpc->Add(hNsigmaP);\r
+  }\r
+  \r
   \r
   \r
   TH2F *hSig = new TH2F("hSigP_TPC",\r
                         "TPC signal vs. p;p [GeV]; TPC signal [arb. units]",\r
                         vX->GetNrows()-1,vX->GetMatrixArray(),\r
                         300,0,300);\r
-  fListQAtpc->Add(hSig);\r
+  fListQAtpc->Add(hSig); //3*AliPID::kSPECIESC\r
 \r
   delete vX;  \r
 }\r
index 630263b..ee658bd 100644 (file)
Binary files a/OADB/COMMON/PID/data/TPCPIDResponse.root and b/OADB/COMMON/PID/data/TPCPIDResponse.root differ
diff --git a/OADB/COMMON/PID/data/TPCvoltageSettings.root b/OADB/COMMON/PID/data/TPCvoltageSettings.root
new file mode 100644 (file)
index 0000000..e49618f
Binary files /dev/null and b/OADB/COMMON/PID/data/TPCvoltageSettings.root differ
index 08f33d2..138c29b 100644 (file)
@@ -68,13 +68,14 @@ fRecoPass(0),
 fRecoPassUser(-1),
 fRun(0),
 fOldRun(0),
-fArrPidResponseMaster(0x0),
-fResolutionCorrection(0x0),
-fTRDPIDResponseObject(0x0),
+fArrPidResponseMaster(NULL),
+fResolutionCorrection(NULL),
+fOADBvoltageMaps(NULL),
+fTRDPIDResponseObject(NULL),
 fTOFtail(1.1),
-fTOFPIDParams(0x0),
-fEMCALPIDParams(0x0),
-fCurrentEvent(0x0),
+fTOFPIDParams(NULL),
+fEMCALPIDParams(NULL),
+fCurrentEvent(NULL),
 fCurrCentrality(0.0),
 fTuneMConData(kFALSE)
 {
@@ -121,13 +122,14 @@ fRecoPass(0),
 fRecoPassUser(other.fRecoPassUser),
 fRun(0),
 fOldRun(0),
-fArrPidResponseMaster(0x0),
-fResolutionCorrection(0x0),
-fTRDPIDResponseObject(0x0),
+fArrPidResponseMaster(NULL),
+fResolutionCorrection(NULL),
+fOADBvoltageMaps(NULL),
+fTRDPIDResponseObject(NULL),
 fTOFtail(1.1),
-fTOFPIDParams(0x0),
-fEMCALPIDParams(0x0),
-fCurrentEvent(0x0),
+fTOFPIDParams(NULL),
+fEMCALPIDParams(NULL),
+fCurrentEvent(NULL),
 fCurrCentrality(0.0),
 fTuneMConData(kFALSE)
 {
@@ -165,13 +167,14 @@ AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other)
     fRecoPassUser=other.fRecoPassUser;
     fRun=0;
     fOldRun=0;
-    fArrPidResponseMaster=0x0;
-    fResolutionCorrection=0x0;
-    fTRDPIDResponseObject=0x0;
-    fEMCALPIDParams=0x0;
+    fArrPidResponseMaster=NULL;
+    fResolutionCorrection=NULL;
+    fOADBvoltageMaps=NULL;
+    fTRDPIDResponseObject=NULL;
+    fEMCALPIDParams=NULL;
     memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
     fTOFtail=1.1;
-    fTOFPIDParams=0x0;
+    fTOFPIDParams=NULL;
     fCurrentEvent=other.fCurrentEvent;
   }
   return *this;
@@ -263,6 +266,19 @@ Float_t AliPIDResponse::NumberOfSigmasTPC(const AliVParticle *vtrack, AliPID::EP
 }
 
 //______________________________________________________________________________
+Float_t AliPIDResponse::NumberOfSigmasTPC( const AliVParticle *vtrack, 
+                                           AliPID::EParticleType type,
+                                           AliTPCPIDResponse::ETPCdEdxSource dedxSource) 
+{
+  //get number of sigmas according the selected TPC gain configuration scenario
+  const AliVTrack *track=static_cast<const AliVTrack*>(vtrack);
+
+  Float_t nSigma=fTPCResponse.GetNumberOfSigmas(track, type, dedxSource);
+
+  return nSigma;
+}
+
+//______________________________________________________________________________
 Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type) const
 {
   //
@@ -696,7 +712,7 @@ void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass, Int_t run)
   //
   fRecoPass=pass;
   
-  fCurrentEvent=0x0;
+  fCurrentEvent=NULL;
   if (!event) return;
   fCurrentEvent=event;
   if (run>0) fRun=run;
@@ -748,6 +764,8 @@ void AliPIDResponse::ExecNewRun()
   
   SetTOFPidResponseMaster();
   InitializeTOFResponse();
+
+  if (fCurrentEvent) fTPCResponse.SetMagField(fCurrentEvent->GetMagneticField());
 }
 
 //_____________________________________________________
@@ -783,7 +801,6 @@ void AliPIDResponse::SetRecoInfo()
     
   fBeamType="PP";
   
-
   TPRegexp reg(".*(LHC1[1-2][a-z]+[0-9]+[a-z_]*)/.*");
   TPRegexp reg12a17(".*(LHC12a17[a-z]+)/.*");
 
@@ -835,6 +852,7 @@ void AliPIDResponse::SetTPCPidResponseMaster()
 {
   //
   // Load the TPC pid response functions from the OADB
+  // Load the TPC voltage maps from OADB
   //
   //don't load twice for the moment
    if (fArrPidResponseMaster) return;
@@ -842,22 +860,37 @@ void AliPIDResponse::SetTPCPidResponseMaster()
 
   //reset the PID response functions
   delete fArrPidResponseMaster;
-  fArrPidResponseMaster=0x0;
+  fArrPidResponseMaster=NULL;
   
   TString fileName(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data()));
+  TFile *f=NULL;
   if (!fCustomTPCpidResponse.IsNull()) fileName=fCustomTPCpidResponse;
   
-  TFile *f=TFile::Open(fileName.Data());
+  TString fileNamePIDresponse(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data()));
+  f=TFile::Open(fileNamePIDresponse.Data());
   if (f && f->IsOpen() && !f->IsZombie()){
     fArrPidResponseMaster=dynamic_cast<TObjArray*>(f->Get("TPCPIDResponse"));
   }
   delete f;
+
+  TString fileNameVoltageMaps(Form("%s/COMMON/PID/data/TPCvoltageSettings.root", fOADBPath.Data()));
+  f=TFile::Open(fileNameVoltageMaps.Data());
+  if (f && f->IsOpen() && !f->IsZombie()){
+    fOADBvoltageMaps=dynamic_cast<AliOADBContainer*>(f->Get("TPCvoltageSettings"));
+  }
+  delete f;
   
   if (!fArrPidResponseMaster){
-    AliFatal(Form("Could not retrieve the TPC pid response from: %s",fileName.Data()));
+    AliFatal(Form("Could not retrieve the TPC pid response from: %s",fileNamePIDresponse.Data()));
     return;
   }
   fArrPidResponseMaster->SetOwner();
+
+  if (!fOADBvoltageMaps)
+  {
+    AliFatal(Form("Could not retrieve the TPC voltage maps from: %s",fileNameVoltageMaps.Data()));
+  }
+  fArrPidResponseMaster->SetOwner();
 }
 
 //______________________________________________________________________________
@@ -887,9 +920,7 @@ void AliPIDResponse::SetTPCParametrisation()
   //
   //reset old splines
   //
-  for (Int_t ispec=0; ispec<AliPID::kSPECIESC; ++ispec){
-    fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,0x0);
-  }
+  fTPCResponse.ResetSplines();
 
   // period
   TString period=fLHCperiod;
@@ -903,59 +934,82 @@ void AliPIDResponse::SetTPCParametrisation()
   if (fArrPidResponseMaster){
     Int_t recopass = fRecoPass;
     if(fTuneMConData) recopass = fRecoPassUser;
-    TObject *grAll=0x0;
     //for MC don't use period information
-//     if (fIsMC) period="[A-Z0-9]*";
+    //if (fIsMC) period="[A-Z0-9]*";
     //for MC use MC period information
-//pattern for the default entry (valid for all particles)
-    TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
-   
-    //loop over entries and filter them
-    for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp){
-      TObject *responseFunction=fArrPidResponseMaster->At(iresp);
-      if (responseFunction==0x0) continue;
-      TString responseName=responseFunction->GetName();
-      
-      if (!reg.MatchB(responseName)) continue;
-      
-      TObjArray *arr=reg.MatchS(responseName);
-      TString particleName=arr->At(1)->GetName();
-      delete arr;
-      if (particleName.IsNull()) continue;
-      if (!grAll&&particleName=="ALL") grAll=responseFunction;
-      else {
-        //find particle id
-        for (Int_t ispec=0; ispec<AliPID::kSPECIESC; ++ispec){
-          TString particle=AliPID::ParticleName(ispec);
-          particle.ToUpper();
-          if ( particle == particleName ){
-            fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,responseFunction);
-            fTPCResponse.SetUseDatabase(kTRUE);
-            AliInfo(Form("Adding graph: %d - %s",ispec,responseFunction->GetName()));
-            // overwrite default with proton spline (for light nuclei)
-            if (ispec==AliPID::kProton) grAll=responseFunction;
-            found=kTRUE;
-            break;
+    //pattern for the default entry (valid for all particles)
+    TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN(_*)([A-Z1-9]*)",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
+
+    //find particle id ang gain scenario
+    for (Int_t igainScenario=0; igainScenario<AliTPCPIDResponse::fgkNumberOfGainScenarios; igainScenario++)
+    {
+      TObject *grAll=NULL;
+      TString gainScenario = AliTPCPIDResponse::GainScenarioName(igainScenario);
+      gainScenario.ToUpper();
+      //loop over entries and filter them
+      for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp)
+      {
+        TObject *responseFunction=fArrPidResponseMaster->At(iresp);
+        if (responseFunction==NULL) continue;
+        TString responseName=responseFunction->GetName();
+         
+        if (!reg.MatchB(responseName)) continue;
+
+        TObjArray *arr=reg.MatchS(responseName); if (!arr) continue;
+        TObject* tmp=NULL;
+        tmp=arr->At(1); if (!tmp) continue;
+        TString particleName=tmp->GetName();
+        tmp=arr->At(3); if (!tmp) continue;
+        TString gainScenarioName=tmp->GetName();
+        delete arr;
+        if (particleName.IsNull()) continue;
+        if (!grAll && particleName=="ALL" && gainScenarioName==gainScenario) grAll=responseFunction;
+        else 
+        {
+          for (Int_t ispec=0; ispec<(AliTPCPIDResponse::fgkNumberOfParticleSpecies); ++ispec)
+          {
+            TString particle=AliPID::ParticleName(ispec);
+            particle.ToUpper();
+            //std::cout<<responseName<<" "<<particle<<" "<<particleName<<" "<<gainScenario<<" "<<gainScenarioName<<std::endl;
+            if ( particle == particleName && gainScenario == gainScenarioName )
+            {
+              fTPCResponse.SetResponseFunction( responseFunction,
+                                                (AliPID::EParticleType)ispec,
+                                                (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
+              fTPCResponse.SetUseDatabase(kTRUE);
+              AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,responseFunction->GetName()));
+              found=kTRUE;
+              // overwrite default with proton spline (for light nuclei)
+              if (ispec==AliPID::kProton) grAll=responseFunction;
+              break;
+            }
           }
         }
       }
-    }
-    
-    //set default response function to all particles which don't have a specific one
-    if (grAll){
-      for (Int_t ispec=0; ispec<AliPID::kSPECIESC; ++ispec){
-        if (!fTPCResponse.GetResponseFunction((AliPID::EParticleType)ispec)){
-          fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,grAll);
-          AliInfo(Form("Adding graph: %d - %s",ispec,grAll->GetName()));
-          found=kTRUE;
+      if (grAll)
+      {
+        for (Int_t ispec=0; ispec<(AliTPCPIDResponse::fgkNumberOfParticleSpecies); ++ispec)
+        {
+          if (!fTPCResponse.GetResponseFunction( (AliPID::EParticleType)ispec,
+                                                 (AliTPCPIDResponse::ETPCgainScenario)igainScenario))
+          {
+              fTPCResponse.SetResponseFunction( grAll,
+                                                (AliPID::EParticleType)ispec,
+                                                (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
+              fTPCResponse.SetUseDatabase(kTRUE);
+              AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,grAll->GetName()));
+              found=kTRUE;
+          }
         }
       }
     }
   }
+  else AliInfo("no fArrPidResponseMaster");
 
   if (!found){
     AliError(Form("No splines found for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
   }
+
   //
   // Setup resolution parametrisation
   //
@@ -970,6 +1024,24 @@ void AliPIDResponse::SetTPCParametrisation()
   fResolutionCorrection=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_SIGMA",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
   
   if (fResolutionCorrection) AliInfo(Form("Setting multiplicity correction function: %s",fResolutionCorrection->GetName()));
+
+  //read in the voltage map
+  TVectorF* gsm = dynamic_cast<TVectorF*>(fOADBvoltageMaps->GetObject(fRun));
+  if (gsm) 
+  {
+    fTPCResponse.SetVoltageMap(*gsm);
+    TString vals;
+    AliInfo(Form("Reading the voltage map for run %d\n",fRun));
+    vals="IROC A: "; for (Int_t i=0; i<18; i++){vals+=Form("%.2f ",(*gsm)[i]);}
+    AliInfo(vals.Data());
+    vals="IROC C: "; for (Int_t i=18; i<36; i++){vals+=Form("%.2f ",(*gsm)[i]);}
+    AliInfo(vals.Data());
+    vals="OROC A: "; for (Int_t i=36; i<54; i++){vals+=Form("%.2f ",(*gsm)[i]);}
+    AliInfo(vals.Data());
+    vals="OROC C: "; for (Int_t i=54; i<72; i++){vals+=Form("%.2f ",(*gsm)[i]);}
+    AliInfo(vals.Data());
+  }
+  else AliInfo("no voltage map, ideal default assumed");
 }
 
 //______________________________________________________________________________
@@ -1044,7 +1116,7 @@ void AliPIDResponse::SetTOFPidResponseMaster()
   //
 
   if (fTOFPIDParams) delete fTOFPIDParams;
-  fTOFPIDParams=0x0;
+  fTOFPIDParams=NULL;
 
   TFile *oadbf = new TFile(Form("%s/COMMON/PID/data/TOFPIDParams.root",fOADBPath.Data()));
   if (oadbf && oadbf->IsOpen()) {
index 307d5b4..82f158b 100644 (file)
@@ -24,6 +24,7 @@
 #include "TNamed.h"
 
 class TF1;
+class AliOADBContainer;
 class TObjArray;
 
 class AliVEvent;
@@ -65,7 +66,7 @@ public:
     kDetPidOk=1,
     kDetMismatch=2
   };
-  
+
   AliITSPIDResponse &GetITSResponse() {return fITSResponse;}
   AliTPCPIDResponse &GetTPCResponse() {return fTPCResponse;}
   AliTOFPIDResponse &GetTOFResponse() {return fTOFResponse;}
@@ -77,6 +78,7 @@ public:
   
   virtual Float_t NumberOfSigmasITS  (const AliVParticle *track, AliPID::EParticleType type) const;
   virtual Float_t NumberOfSigmasTPC  (const AliVParticle *track, AliPID::EParticleType type) const;
+  virtual Float_t NumberOfSigmasTPC  (const AliVParticle *track, AliPID::EParticleType type, AliTPCPIDResponse::ETPCdEdxSource dedxSource);
   virtual Float_t NumberOfSigmasEMCAL(const AliVParticle *track, AliPID::EParticleType type, Double_t &eop, Double_t showershape[4]) const;
   virtual Float_t NumberOfSigmasTOF  (const AliVParticle *track, AliPID::EParticleType type) const = 0;
   virtual Float_t NumberOfSigmasEMCAL(const AliVParticle *track, AliPID::EParticleType type) const;
@@ -160,6 +162,7 @@ private:
   
   TObjArray *fArrPidResponseMaster;    //!  TPC pid splines
   TF1       *fResolutionCorrection;    //! TPC resolution correction
+  AliOADBContainer* fOADBvoltageMaps;  //! container with the voltage maps
 
   AliTRDPIDResponseObject *fTRDPIDResponseObject; //! TRD PID Response Object
   UInt_t fTRDslicesForPID[2];           //! TRD PID slices
index b0bf074..9a9def6 100644 (file)
 // With many additions and modifications suggested by
 //      Alexander Kalweit, GSI, alexander.philipp.kalweit@cern.ch
 //      Dariusz Miskowiec, GSI, D.Miskowiec@gsi.de
+// ...and some modifications by
+//      Mikolaj Krzewicki, GSI, mikolaj.krzewicki@cern.ch
 //-----------------------------------------------------------------
 
 #include <TGraph.h>
 #include <TObjArray.h>
 #include <TSpline.h>
+#include <TBits.h>
 #include <TMath.h>
 
+#include <AliLog.h>
 #include "AliExternalTrackParam.h"
-
+#include "AliVTrack.h"
 #include "AliTPCPIDResponse.h"
+#include "AliTPCdEdxInfo.h"
+
+ClassImp(AliTPCPIDResponse);
 
-ClassImp(AliTPCPIDResponse)
+const char* AliTPCPIDResponse::fgkGainScenarioName[fgkNumberOfGainScenarios+1]=
+{
+  "", //default - no name
+  "1", //all high
+  "2", //OROC only
+  "unknown"
+};
 
 //_________________________________________________________________________
 AliTPCPIDResponse::AliTPCPIDResponse():
+  TNamed(),
   fMIP(50.),
-  fRes0(0.07),
-  fResN2(0.),
+  fRes0(),
+  fResN2(),
   fKp1(0.0283086),
   fKp2(2.63394e+01),
   fKp3(5.04114e-11),
   fKp4(2.12543),
   fKp5(4.88663),
   fUseDatabase(kFALSE),
-  fResponseFunctions(AliPID::kUnknown+1)
+  fResponseFunctions(fgkNumberOfParticleSpecies*fgkNumberOfGainScenarios),
+  fVoltageMap(72),
+  fLowGainIROCthreshold(-40),
+  fBadIROCthreshhold(-70),
+  fLowGainOROCthreshold(-40),
+  fBadOROCthreshhold(-40),
+  fMaxBadLengthFraction(0.5),
+  fCurrentResponseFunction(NULL),
+  fCurrentdEdx(0.0),
+  fCurrentNPoints(0),
+  fCurrentGainScenario(kGainScenarioInvalid),
+  fMagField(0.)
 {
   //
   //  The default constructor
   //
+  for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=0.07;fResN2[i]=0.0;}
 }
 
 //_________________________________________________________________________
 AliTPCPIDResponse::AliTPCPIDResponse(const Double_t *param):
+  TNamed(),
   fMIP(param[0]),
-  fRes0(param[1]),
-  fResN2(param[2]),
+  fRes0(),
+  fResN2(),
   fKp1(0.0283086),
   fKp2(2.63394e+01),
   fKp3(5.04114e-11),
   fKp4(2.12543),
   fKp5(4.88663),
   fUseDatabase(kFALSE),
-  fResponseFunctions(AliPID::kUnknown+1)
+  fResponseFunctions(fgkNumberOfParticleSpecies*fgkNumberOfGainScenarios),
+  fVoltageMap(72),
+  fLowGainIROCthreshold(-40),
+  fBadIROCthreshhold(-70),
+  fLowGainOROCthreshold(-40),
+  fBadOROCthreshhold(-40),
+  fMaxBadLengthFraction(0.5),
+  fCurrentResponseFunction(NULL),
+  fCurrentdEdx(0.0),
+  fCurrentNPoints(0),
+  fCurrentGainScenario(kGainScenarioInvalid),
+  fMagField(0.)
 {
   //
   //  The main constructor
   //
+  for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=param[1];fResN2[i]=param[2];}
+}
+
+//_________________________________________________________________________
+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),
+  fCurrentResponseFunction(NULL),
+  fCurrentdEdx(0.0),
+  fCurrentNPoints(0),
+  fCurrentGainScenario(kGainScenarioInvalid),
+  fMagField(that.fMagField)
+{
+  //copy ctor
+  for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=that.fRes0[i];fResN2[i]=that.fResN2[i];}
+}
+
+//_________________________________________________________________________
+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;
+  for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=that.fRes0[i];fResN2[i]=that.fResN2[i];}
+  return *this;
 }
 
 //_________________________________________________________________________
@@ -104,13 +197,13 @@ void AliTPCPIDResponse::SetBetheBlochParameters(Double_t kp1,
   fKp4=kp4;
   fKp5=kp5;
 }
+
 //_________________________________________________________________________
 void AliTPCPIDResponse::SetSigma(Float_t res0, Float_t resN2) {
   //
   // Set the relative resolution  sigma_rel = res0 * sqrt(1+resN2/npoint)
   //
-  fRes0 = res0;
-  fResN2 = resN2;
+  for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=res0;fResN2[i]=resN2;}
 }
 
 //_________________________________________________________________________
@@ -127,9 +220,9 @@ Double_t AliTPCPIDResponse::GetExpectedSignal(const Float_t mom,
   //
   
   Double_t mass=AliPID::ParticleMassZ(n);
-  if (!fUseDatabase||fResponseFunctions.GetEntriesFast()>AliPID::kUnknown) return Bethe(mom/mass);
+  if (!fUseDatabase) return Bethe(mom/mass);
   //
-  TSpline3 * responseFunction = (TSpline3 *) fResponseFunctions.UncheckedAt(n);
+  const TSpline3 * responseFunction = (TSpline3 *) fResponseFunctions.UncheckedAt(n);
 
   if (!responseFunction) return Bethe(mom/mass);
   //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...)
@@ -141,16 +234,521 @@ Double_t AliTPCPIDResponse::GetExpectedSignal(const Float_t mom,
 
 //_________________________________________________________________________
 Double_t AliTPCPIDResponse::GetExpectedSigma(const Float_t mom, 
-                                            const Int_t nPoints,
-                                         AliPID::EParticleType n) const {
+                                            const Int_t nPoints,
+                                             AliPID::EParticleType n) const {
   //
   // 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(Double_t mom,
+                                              AliPID::EParticleType species,
+                                              const TSpline3* responseFunction) 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. 
+  // This can be improved. By taking into account the number of
+  // assigned clusters and/or the track dip angle, for example.  
   //
   
+  
+  Double_t mass=AliPID::ParticleMass(species);
+  
+  if (!fUseDatabase) return Bethe(mom/mass);
+  
+  if (!responseFunction) return Bethe(mom/mass);
+  return fMIP*responseFunction->Eval(mom/mass);
+}
+
+//_________________________________________________________________________
+Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track,
+                                              AliPID::EParticleType species,
+                                              ETPCdEdxSource dedxSource) 
+{
+  // 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.  
+  //
+  
+  Double_t mom = track->GetTPCmomentum();
+  Double_t mass=AliPID::ParticleMass(species);
+  
+  if (!fUseDatabase) return Bethe(mom/mass);
+  
+  const TSpline3* responseFunction = GetResponseFunction(track,species,dedxSource);
+  if (!responseFunction) return Bethe(mom/mass);
+  return fMIP*responseFunction->Eval(mom/mass);
+
+}
+  
+//_________________________________________________________________________
+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 ) 
+{
+  //the splines are stored in an array, different scenarios
+
+  if (ResponseFunctiondEdxN(track,
+                            species,
+                            dedxSource))
+    return fCurrentResponseFunction;
+  
+  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( Double_t mom,
+                                              Int_t nPoints,
+                                              AliPID::EParticleType species,
+                                              ETPCgainScenario gainScenario,
+                                              const TSpline3* responseFunction) const
+{
+  // Calculates the expected sigma of the PID signal as the function of 
+  // the information stored in the track, for the specified particle type 
+  // and dedx scenario
+  //
+  
+  if (!responseFunction) return 999;
   if (nPoints != 0) 
-    return GetExpectedSignal(mom,n)*fRes0*sqrt(1. + fResN2/nPoints);
+    return GetExpectedSignal(mom,species,responseFunction) *
+           fRes0[gainScenario] * sqrt(1. + fResN2[gainScenario]/nPoints);
+  else
+    return GetExpectedSignal(mom,species,responseFunction)*fRes0[gainScenario];
+}
+
+//_________________________________________________________________________
+Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track, 
+                                             AliPID::EParticleType species,
+                                             ETPCdEdxSource dedxSource) 
+{
+  // Calculates the expected sigma of the PID signal as the function of 
+  // the information stored in the track, for the specified particle type 
+  // and dedx scenario
+  //
+  
+  if (!ResponseFunctiondEdxN(track,
+                             species,
+                             dedxSource)) return 998; //TODO: better handling!
+
+  return GetExpectedSigma( track->GetTPCmomentum(),
+                           fCurrentNPoints,
+                           species,
+                           fCurrentGainScenario,
+                           fCurrentResponseFunction );
+}
+
+//_________________________________________________________________________
+Float_t AliTPCPIDResponse::GetNumberOfSigmas(const AliVTrack* track, 
+                             AliPID::EParticleType species,
+                             ETPCdEdxSource dedxSource) 
+{
+  //calculates the number of sigmas of the PID signal from the expected value
+  //for a gicven particle species in the presence of multiple gain scenarios
+  //inside the TPC
+
+  if (!ResponseFunctiondEdxN(track,
+                             species,
+                             dedxSource)) return 997; //TODO: better handling!
+
+  Double_t mom = track->GetTPCmomentum();
+  Double_t bethe=GetExpectedSignal(mom,species,fCurrentResponseFunction);
+  Double_t sigma=GetExpectedSigma( mom,
+                                   fCurrentNPoints,
+                                   species,
+                                   fCurrentGainScenario,
+                                   fCurrentResponseFunction );
+  return (fCurrentdEdx-bethe)/sigma;
+}
+
+//_________________________________________________________________________
+Bool_t AliTPCPIDResponse::ResponseFunctiondEdxN( const AliVTrack* track, 
+                                                 AliPID::EParticleType species,
+                                                 ETPCdEdxSource dedxSource ) 
+{
+  // 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
+  
+  Double32_t signal[4]; //0: IROC, 1: OROC medium, 2:OROC long, 3: OROC all (def. truncation used)
+  Char_t ncl[3];        //same
+  Char_t nrows[3];      //same
+  AliTPCdEdxInfo* dEdxInfo = track->GetTPCdEdxInfo();
+  
+  if (!dEdxInfo && dedxSource!=kdEdxDefault)  //in one case its ok if we dont have the info
+  {
+    AliError("AliTPCdEdxInfo not available");
+    InvalidateCurrentValues();
+    return kFALSE;
+  }
+
+  if (dEdxInfo) dEdxInfo->GetTPCSignalRegionInfo(signal,ncl,nrows);
+
+  //check if we cross a bad OROC in which case we reject
+  EChamberStatus trackStatus = TrackStatus(track,2);
+  if (trackStatus==kChamberOff || trackStatus==kChamberLowGain)
+  {
+    InvalidateCurrentValues();
+    return kFALSE;
+  }
+
+  switch (dedxSource)
+  {
+    case kdEdxDefault:
+      {
+        fCurrentdEdx = track->GetTPCsignal();
+        fCurrentNPoints = track->GetTPCsignalN();
+        fCurrentGainScenario = kDefault;
+        break;
+      }
+    case kdEdxOROC:
+      {
+        fCurrentdEdx = signal[3];
+        fCurrentNPoints = ncl[2]+ncl[1];
+        fCurrentGainScenario = kOROChigh;
+        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)
+        {
+          fCurrentdEdx = signal[3];
+          fCurrentNPoints = ncl[2]+ncl[1];
+          fCurrentGainScenario = kOROChigh;
+        }
+        else
+        {
+          fCurrentdEdx = track->GetTPCsignal();
+          fCurrentNPoints = track->GetTPCsignalN();
+          fCurrentGainScenario = kALLhigh;
+        }
+        break;
+      }
+    default:
+      {
+         fCurrentdEdx = 0.;
+         fCurrentNPoints = 0;
+         fCurrentGainScenario = kGainScenarioInvalid;
+         return kFALSE;
+      }
+  }
+  TObject* obj = fResponseFunctions.UncheckedAt(ResponseFunctionIndex(species,fCurrentGainScenario));
+  fCurrentResponseFunction = dynamic_cast<TSpline3*>(obj); //TODO:maybe static cast?
+  return kTRUE;
+}
+
+//_________________________________________________________________________
+Bool_t AliTPCPIDResponse::sectorNumbersInOut(const AliVTrack* track,
+                                             Double_t innerRadius,
+                                             Double_t outerRadius,
+                                             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)
+  
+  Double_t trackPositionInner[3]; 
+  Double_t trackPositionOuter[3]; 
+
+  Bool_t trackAtInner=kTRUE;
+  Bool_t trackAtOuter=kTRUE;
+  const AliExternalTrackParam* ip = track->GetInnerParam();
+  
+  //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.
+  if (ip) track=ip;
+
+  if (!track->GetXYZAt(innerRadius, fMagField, trackPositionInner))
+    trackAtInner=kFALSE;
+
+  if (!track->GetXYZAt(outerRadius, fMagField, trackPositionOuter))
+    trackAtOuter=kFALSE;
+
+  if (!trackAtInner)
+  {
+    //if we dont even enter inner radius we do nothing
+    inphi=0.0;
+    outphi=0.0;
+    in=0;
+    out=0;
+    return kFALSE;
+  }
+
+  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 kFALSE;
+    }
+  }
+
+  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;
+  if (!sectorNumbersInOut(track, 
+                          innerRadius, 
+                          outerRadius, 
+                          inphi, 
+                          outphi, 
+                          in, 
+                          out)) return kChamberOff;
+
+  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
-    return GetExpectedSignal(mom,n)*fRes0;
+  {
+    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);
+
+  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;
+    if (v<=fLowGainIROCthreshold && v>fBadIROCthreshhold) trackLengthInLowGain+=deltaPhi;
+  }
+
+  //for now low gain and bad (off) chambers are treated equally
+  Float_t 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;
 }
+
+//_____________________________________________________________________________
+void AliTPCPIDResponse::InvalidateCurrentValues()
+{
+  //make the "current" stored values from the last processed track invalid
+  fCurrentResponseFunction=NULL;
+  fCurrentdEdx=0.;
+  fCurrentNPoints=0;
+  fCurrentGainScenario=kGainScenarioInvalid;
+}
+
+//_____________________________________________________________________________
+Bool_t AliTPCPIDResponse::TrackApex(const AliVTrack* track, Float_t magField, Double_t position[3]) const
+{
+  //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;
+}
+
index 94673a0..eb6a3e3 100644 (file)
 //      Dariusz Miskowiec, GSI, D.Miskowiec@gsi.de
 //-------------------------------------------------------
 #include <Rtypes.h>
+
+#include <TNamed.h>
+#include <TVectorF.h>
 #include <TObjArray.h>
 
 #include "AliPID.h"
 
-class AliTPCPIDResponse {
+class AliVTrack;
+class TSpline3;
+
+class AliTPCPIDResponse: public TNamed {
 public:
   AliTPCPIDResponse();
   AliTPCPIDResponse(const Double_t *param);
+  AliTPCPIDResponse(const AliTPCPIDResponse&);
+  AliTPCPIDResponse& operator=(const AliTPCPIDResponse&);
   virtual ~AliTPCPIDResponse() {}
+
+  enum EChamberStatus {
+    kChamberOff=0,
+    kChamberHighGain=1,
+    kChamberLowGain=2,
+    kChamberInvalid=3
+  };
+  
+  enum ETPCgainScenario {
+    kDefault= 0,
+    kALLhigh = 1,
+    kOROChigh = 2,
+    kGainScenarioInvalid = 3
+  };
+
+  static const Int_t fgkNumberOfParticleSpecies=AliPID::kSPECIESC;
+  static const Int_t fgkNumberOfGainScenarios=3;
+  static const Int_t fgkNumberOfdEdxSourceScenarios=3;
+
+  enum ETPCdEdxSource {
+    kdEdxDefault=0,        // use combined dEdx from IROC+OROC (assumes ideal detector)
+    kdEdxOROC=1,       // use only OROC
+    kdEdxHybrid=2,   // Use IROC+OROC dEdx only where IROCS are good (high gain), otherwise fall back to OROC only
+    kdEdxInvalid=3     //invalid
+  };
+
   void SetSigma(Float_t res0, Float_t resN2);
   void SetBetheBlochParameters(Double_t kp1,
                                Double_t kp2,
@@ -31,17 +65,80 @@ public:
   void SetMip(Float_t mip) { fMIP = mip; } // Set overall normalisation; mean dE/dx for MIP
   Double_t Bethe(Double_t bg) const;
   void SetUseDatabase(Bool_t useDatabase) { fUseDatabase = useDatabase;}
+  Bool_t GetUseDatabase() const { return fUseDatabase;}
   
   void SetResponseFunction(AliPID::EParticleType type, TObject * const o) { fResponseFunctions.AddAt(o,(Int_t)type); }
   const TObject * GetResponseFunction(AliPID::EParticleType type) { return fResponseFunctions.At((Int_t)type); }
+  void SetVoltage(Int_t n, Float_t v) {fVoltageMap[n]=v;}
+  void SetVoltageMap(const TVectorF& a) {fVoltageMap=a;} //resets ownership, ~ will not delete contents
+  Float_t GetVoltage(Int_t n) const {return fVoltageMap[n];}
+  void SetLowGainIROCthreshold(Float_t v) {fLowGainIROCthreshold=v;}
+  void SetBadIROCthreshold(Float_t v) {fBadIROCthreshhold=v;}
+  void SetLowGainOROCthreshold(Float_t v) {fLowGainOROCthreshold=v;}
+  void SetBadOROCthreshold(Float_t v) {fBadOROCthreshhold=v;}
+  void SetMaxBadLengthFraction(Float_t f) {fMaxBadLengthFraction=f;}
+
+  void SetMagField(Double_t mf) { fMagField=mf; }
   
+  //NEW
+  void SetSigma(Float_t res0, Float_t resN2, ETPCgainScenario gainScenario );
+  Double_t GetExpectedSignal( Double_t momentum,
+                              AliPID::EParticleType species,
+                              const TSpline3* responseFunction ) const;
+  Double_t GetExpectedSignal( const AliVTrack* track,
+                              AliPID::EParticleType species,
+                              ETPCdEdxSource dedxSource );
+  Double_t GetExpectedSigma( const AliVTrack* track, 
+                             AliPID::EParticleType species,
+                             ETPCdEdxSource dedxSource );
+  Double_t GetExpectedSigma( Double_t mom,
+                             Int_t nPoints,
+                             AliPID::EParticleType species,
+                             ETPCgainScenario gainScenario,
+                             const TSpline3* responseFunction) const;
+  Float_t GetNumberOfSigmas( const AliVTrack* track,
+                             AliPID::EParticleType species,
+                             ETPCdEdxSource dedxSource );
+
+  void SetResponseFunction(TObject* o,
+                           AliPID::EParticleType type,
+                           ETPCgainScenario gainScenario);
+  void Print(Option_t* option="") const;
+  TSpline3* GetResponseFunction( AliPID::EParticleType species,
+                                 ETPCgainScenario gainScenario ) const;
+  TSpline3* GetResponseFunction( const AliVTrack* track,
+                                 AliPID::EParticleType species,
+                                 ETPCdEdxSource dedxSource );
+  Bool_t ResponseFunctiondEdxN(const AliVTrack* track, 
+                               AliPID::EParticleType species,
+                               ETPCdEdxSource dedxSource);
+  Bool_t sectorNumbersInOut(const AliVTrack* track, 
+                            Double_t innerRadius, Double_t outerRadius, 
+                            Float_t& phiIn, Float_t& phiOut, 
+                            Int_t& in, Int_t& out ) const;
+  AliTPCPIDResponse::EChamberStatus TrackStatus(const AliVTrack* track, Int_t layer) const;
+  Float_t MaxClusterRadius(const AliVTrack* track) const;
+  Bool_t TrackApex(const AliVTrack* track, Float_t magField, Double_t position[3]) const;
+  static const char* GainScenarioName(Int_t n) {return fgkGainScenarioName[(n>fgkNumberOfdEdxSourceScenarios)?fgkNumberOfdEdxSourceScenarios+1:n];}
+  Int_t ResponseFunctionIndex( AliPID::EParticleType species,
+                               ETPCgainScenario gainScenario ) const;
+  void ResetSplines();
+
+  void InvalidateCurrentValues();
+  TSpline3* GetCurrentResponseFunction() const {return fCurrentResponseFunction;}
+  Double_t GetCurrentdEdx() const {return fCurrentdEdx;}
+  Int_t GetCurrentNPoints() const {return fCurrentNPoints;}
+  ETPCgainScenario GetCurrentGainScenario() const {return fCurrentGainScenario;}
+
+  //OLD
   Double_t GetExpectedSignal(const Float_t mom,
                      AliPID::EParticleType n=AliPID::kKaon) const;
   Double_t GetExpectedSigma(const Float_t mom, const Int_t nPoints,
-                     AliPID::EParticleType n=AliPID::kKaon) const;
-  Float_t  GetNumberOfSigmas(const Float_t mom, const Float_t dEdx, 
-                            const Int_t nPoints,
-                     AliPID::EParticleType n=AliPID::kKaon) const {
+                            AliPID::EParticleType n=AliPID::kKaon) const;
+  Float_t  GetNumberOfSigmas(const Float_t mom, 
+                             const Float_t dEdx, 
+                                              const Int_t nPoints,
+                             AliPID::EParticleType n=AliPID::kKaon) const {
 
     Double_t bethe=GetExpectedSignal(mom,n);
     Double_t sigma=GetExpectedSigma(mom,nPoints,n);
@@ -49,13 +146,15 @@ public:
   }
 
   Double_t GetMIP() const { return fMIP;} 
-  Float_t  GetRes0()  const { return fRes0;  }
-  Float_t  GetResN2() const { return fResN2; }
+  Float_t  GetRes0()  const { return fRes0[0];  }
+  Float_t  GetResN2() const { return fResN2[0]; }
+  Float_t  GetRes0(ETPCgainScenario s)  const { return fRes0[s];  }
+  Float_t  GetResN2(ETPCgainScenario s) const { return fResN2[s]; }
 
 private:
   Float_t fMIP;          // dEdx for MIP
-  Float_t fRes0;         // relative dEdx resolution  rel sigma = fRes0*sqrt(1+fResN2/npoint)
-  Float_t fResN2;        // relative Npoint dependence rel  sigma = fRes0*sqrt(1+fResN2/npoint)
+  Float_t fRes0[fgkNumberOfGainScenarios];  // relative dEdx resolution  rel sigma = fRes0*sqrt(1+fResN2/npoint)
+  Float_t fResN2[fgkNumberOfGainScenarios]; // relative Npoint dependence rel  sigma = fRes0*sqrt(1+fResN2/npoint)
 
   Double_t fKp1;   // Parameters
   Double_t fKp2;   //    of
@@ -64,9 +163,26 @@ private:
   Double_t fKp5;   // formula
 
   Bool_t fUseDatabase; // flag if fine-tuned database-response or simple ALEPH BB should be used
+  
   TObjArray fResponseFunctions; //! ObjArray of response functions individually for each particle
+  TVectorF fVoltageMap; //!stores a map of voltages wrt nominal for all chambers
+  Float_t fLowGainIROCthreshold;  //voltage threshold below which the IROC is considered low gain
+  Float_t fBadIROCthreshhold;     //voltage threshold for bad IROCS
+  Float_t fLowGainOROCthreshold;  //voltage threshold below which the OROC is considered low gain
+  Float_t fBadOROCthreshhold;     //voltage threshold for bad OROCS
+  Float_t fMaxBadLengthFraction;  //the maximum allowed fraction of track length in a bad sector.
+
+  TSpline3* fCurrentResponseFunction;      //!response function for current track
+  Double_t fCurrentdEdx;                  //!dEdx for currently processed track
+  Int_t fCurrentNPoints;                  //!number of points used for dEdx calculation for current track
+  ETPCgainScenario fCurrentGainScenario;  //!gain scenario used for current track
+  Int_t sectorNumber(Double_t phi) const;
+
+  Double_t fMagField;  //! Magnetic field
+
+  static const char* fgkGainScenarioName[fgkNumberOfGainScenarios+1];
 
-  ClassDef(AliTPCPIDResponse,3)   // TPC PID class
+  ClassDef(AliTPCPIDResponse,4)   // TPC PID class
 };
 
 #endif