use Eval const outside of the boundaries
[u/mrichter/AliRoot.git] / TPC / AliTPCseed.cxx
index 0d59602..cbb0c20 100644 (file)
 #include "AliTPCCalROC.h"
 #include "AliTPCcalibDB.h"
 #include "AliTPCParam.h"
-
+#include "AliMathBase.h"
+#include "AliTPCTransform.h"
+#include "AliSplineFit.h"
+#include "AliCDBManager.h"
+#include "AliTPCcalibDButil.h"
 
 
 ClassImp(AliTPCseed)
@@ -575,27 +579,32 @@ Float_t AliTPCseed::CookdEdx(Double_t low, Double_t up,Int_t i1, Int_t i2, Bool_
   //-----------------------------------------------------------------
   // This funtion calculates dE/dX within the "low" and "up" cuts.
   //-----------------------------------------------------------------
+  // CookdEdxAnalytical(Double_t low, Double_t up, Int_t type, Int_t i1, Int_t i2, Int_t returnVal)
   AliTPCParam *param = AliTPCcalibDB::Instance()->GetParameters();
+  
   Int_t row0 = param->GetNRowLow();
   Int_t row1 = row0+param->GetNRowUp1();
   Int_t row2 = row1+param->GetNRowUp2();
+  const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
+  Int_t useTot = 0;
+  if (recoParam) useTot = (recoParam->GetUseTotCharge())? 0:1;
   //
   //
   //
-  fDEDX[0]      = CookdEdxNorm(low,up,0 ,i1  ,i2,  kTRUE,kFALSE,2,0);
-  fDEDX[1]      = CookdEdxNorm(low,up,0 ,0   ,row0,kTRUE,kFALSE,2,0);
-  fDEDX[2]      = CookdEdxNorm(low,up,0 ,row0,row1,kTRUE,kFALSE,2,0);
-  fDEDX[3]      = CookdEdxNorm(low,up,0 ,row1,row2,kTRUE,kFALSE,2,0);
+  fDEDX[0]      = CookdEdxAnalytical(low,up,useTot ,i1  ,i2,   0);
+  fDEDX[1]      = CookdEdxAnalytical(low,up,useTot ,0   ,row0, 0);
+  fDEDX[2]      = CookdEdxAnalytical(low,up,useTot ,row0,row1, 0);
+  fDEDX[3]      = CookdEdxAnalytical(low,up,useTot ,row1,row2, 0);
   //
-  fSDEDX[0]     = CookdEdxNorm(low,up,0 ,i1  ,i2,  kTRUE,kFALSE,2,1);
-  fSDEDX[1]     = CookdEdxNorm(low,up,0 ,0   ,row0,kTRUE,kFALSE,2,1);
-  fSDEDX[2]     = CookdEdxNorm(low,up,0 ,row0,row1,kTRUE,kFALSE,2,1);
-  fSDEDX[3]     = CookdEdxNorm(low,up,0 ,row1,row2,kTRUE,kFALSE,2,1);
+  fSDEDX[0]     = CookdEdxAnalytical(low,up,useTot ,i1  ,i2,   1);
+  fSDEDX[1]     = CookdEdxAnalytical(low,up,useTot ,0   ,row0, 1);
+  fSDEDX[2]     = CookdEdxAnalytical(low,up,useTot ,row0,row1, 1);
+  fSDEDX[3]     = CookdEdxAnalytical(low,up,useTot ,row1,row2, 1);
   //
-  fNCDEDX[0]    = TMath::Nint(CookdEdxNorm(low,up,0 ,i1  ,i2,  kTRUE,kFALSE,2,2));
-  fNCDEDX[1]    = TMath::Nint(CookdEdxNorm(low,up,0 ,0   ,row0,kTRUE,kFALSE,2,2));
-  fNCDEDX[2]    = TMath::Nint(CookdEdxNorm(low,up,0 ,row0,row1,kTRUE,kFALSE,2,2));
-  fNCDEDX[3]    = TMath::Nint(CookdEdxNorm(low,up,0 ,row1,row2,kTRUE,kFALSE,2,2));
+  fNCDEDX[0]    = TMath::Nint(CookdEdxAnalytical(low,up,useTot ,i1  ,i2,   2));
+  fNCDEDX[1]    = TMath::Nint(CookdEdxAnalytical(low,up,useTot ,0   ,row0, 2));
+  fNCDEDX[2]    = TMath::Nint(CookdEdxAnalytical(low,up,useTot ,row0,row1, 2));
+  fNCDEDX[3]    = TMath::Nint(CookdEdxAnalytical(low,up,useTot ,row1,row2, 2));
 
   SetdEdx(fDEDX[0]);
   return fDEDX[0];
@@ -833,18 +842,6 @@ Float_t AliTPCseed::CookdEdx(Double_t low, Double_t up,Int_t i1, Int_t i2, Bool_
 //   SetdEdx(dedx);
 //   return dedx;
 }
-Double_t AliTPCseed::Bethe(Double_t bg){
-  //
-  // This is the Bethe-Bloch function normalised to 1 at the minimum
-  //
-  Double_t bg2=bg*bg;
-  Double_t bethe;
-  if (bg<3.5e1) 
-    bethe=(1.+ bg2)/bg2*(log(5940*bg2) - bg2/(1.+ bg2));
-  else // Density effect ( approximately :) 
-    bethe=1.15*(1.+ bg2)/bg2*(log(3.5*5940*bg) - bg2/(1.+ bg2));
-  return bethe/11.091;
-}
 
 void AliTPCseed::CookPID()
 {
@@ -861,7 +858,7 @@ void AliTPCseed::CookPID()
     Double_t mass=AliPID::ParticleMass(j);
     Double_t mom=GetP();
     Double_t dedx=fdEdx/fMIP;
-    Double_t bethe=Bethe(mom/mass); 
+    Double_t bethe=AliMathBase::BetheBlochAleph(mom/mass); 
     Double_t sigma=fRes*bethe;
     if (sigma>0.001){
       if (TMath::Abs(dedx-bethe) > fRange*sigma) {
@@ -1038,8 +1035,8 @@ Float_t  AliTPCseed::CookdEdxNorm(Double_t low, Double_t up, Int_t type, Int_t i
 
     if (padNorm==1){
       //taken from OCDB
-      if (type==0 && parcl->fQpadTnorm) corrPadType = (*parcl->fQpadTnorm)[ipad];
-      if (type==1 && parcl->fQpadTnorm) corrPadType = (*parcl->fQpadMnorm)[ipad];
+      if (type==0 && parcl->QpadTnorm()) corrPadType = (*parcl->QpadTnorm())[ipad];
+      if (type==1 && parcl->QpadMnorm()) corrPadType = (*parcl->QpadMnorm())[ipad];
 
     }
     if (padNorm==2){
@@ -1077,6 +1074,30 @@ Float_t  AliTPCseed::CookdEdxNorm(Double_t low, Double_t up, Int_t type, Int_t i
   }
   Float_t mean =suma/sumn;
   Float_t rms  =TMath::Sqrt(TMath::Abs(suma2/sumn-mean*mean));
+  //
+  // do time-dependent correction for pressure and temperature variations
+  UInt_t runNumber = 1;
+  Float_t corrTimeGain = 1;
+  AliTPCTransform * trans = AliTPCcalibDB::Instance()->GetTransform();
+  const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
+  if (trans && recoParam->GetUseGainCorrectionTime()>0) {
+    runNumber = trans->GetCurrentRunNumber();
+    //AliTPCcalibDB::Instance()->SetRun(runNumber);
+    TObjArray * timeGainSplines = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(runNumber);
+    if (timeGainSplines) {
+      UInt_t time = trans->GetCurrentTimeStamp();
+      AliSplineFit * fitMIP = (AliSplineFit *) timeGainSplines->At(0);
+      AliSplineFit * fitFPcosmic = (AliSplineFit *) timeGainSplines->At(1);
+      if (fitMIP) {
+       corrTimeGain = AliTPCcalibDButil::EvalGraphConst(fitMIP, time);/*fitMIP->Eval(time);*/
+      } else {
+       if (fitFPcosmic) corrTimeGain = AliTPCcalibDButil::EvalGraphConst(fitFPcosmic, time);/*fitFPcosmic->Eval(time);*/ 
+      }
+    }
+  }
+  mean /= corrTimeGain;
+  rms /= corrTimeGain;
+  //
   if (returnVal==1) return rms;
   if (returnVal==2) return ncl;
   return mean;
@@ -1113,7 +1134,7 @@ Float_t  AliTPCseed::CookdEdxAnalytical(Double_t low, Double_t up, Int_t type, I
   Float_t gainGG      = 1;  // gas gain factor -always enabled
   Float_t gainPad     = 1;  // gain map  - used always
   Float_t corrPos     = 1;  // local position correction - if posNorm enabled
-  Float_t corrNorm    = 1;  // normalization factor - set Q to channel 50
+
   //   
   //
   //
@@ -1214,34 +1235,36 @@ Float_t  AliTPCseed::CookdEdxAnalytical(Double_t low, Double_t up, Int_t type, I
   }
   Float_t mean =suma/sumn;
   Float_t rms  =TMath::Sqrt(TMath::Abs(suma2/sumn-mean*mean));
+  //
+  // do time-dependent correction for pressure and temperature variations
+  UInt_t runNumber = 1;
+  Float_t corrTimeGain = 1;
+  AliTPCTransform * trans = AliTPCcalibDB::Instance()->GetTransform();
+  const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
+  if (trans && recoParam->GetUseGainCorrectionTime()>0) {
+    runNumber = trans->GetCurrentRunNumber();
+    //AliTPCcalibDB::Instance()->SetRun(runNumber);
+    TObjArray * timeGainSplines = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(runNumber);
+    if (timeGainSplines) {
+      UInt_t time = trans->GetCurrentTimeStamp();
+      AliSplineFit * fitMIP = (AliSplineFit *) timeGainSplines->At(0);
+      AliSplineFit * fitFPcosmic = (AliSplineFit *) timeGainSplines->At(1);
+      if (fitMIP) {
+       corrTimeGain = fitMIP->Eval(time);
+      } else {
+       if (fitFPcosmic) corrTimeGain = fitFPcosmic->Eval(time); // This value describes the ratio FP-to-MIP, hardwired for the moment
+      }
+    }
+  }
+  mean /= corrTimeGain;
+  rms /= corrTimeGain;
+  //
   if (returnVal==1) return rms;
   if (returnVal==2) return ncl;
   return mean;
 }
 
 
-Double_t AliTPCseed::BetheMass(Double_t mass){
-  //
-  // return bethe-bloch
-  //
-  Float_t bg= P()/mass; 
-  const Double_t kp1=0.76176e-1;
-  const Double_t kp2=10.632;
-  const Double_t kp3=0.13279e-4;
-  const Double_t kp4=1.8631;
-  const Double_t kp5=1.9479;
-
-  Double_t dbg = (Double_t) bg;
-
-  Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);
-
-  Double_t aa = TMath::Power(beta,kp4);
-  Double_t bb = TMath::Power(1./dbg,kp5);
-
-  bb=TMath::Log(kp3+bb);
-  
-  return ((Float_t)((kp2-aa-bb)*kp1/aa));
-}
 
 
 Float_t  AliTPCseed::CookShape(Int_t type){
@@ -1388,6 +1411,7 @@ Int_t  AliTPCseed::RefitTrack(AliTPCseed *seed, AliExternalTrackParam * parin, A
   //
   if (parin) (*parin)=paramIn;
   if (parout) (*parout)=paramOut;
+  delete track;
   return ncl;
 }