]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TPC/AliTPCseed.cxx
Patch for the tracker
[u/mrichter/AliRoot.git] / TPC / AliTPCseed.cxx
index 3b784d70d2b146931774116c362cc2865328dc25..c05f1f82c3e87ef25a1c15ac7fd931ef895a4627 100644 (file)
@@ -22,6 +22,7 @@
 //      Origin: Marian Ivanov, CERN, Marian.Ivanov@cern.ch
 //-----------------------------------------------------------------
 #include "TClonesArray.h"
+#include "TGraphErrors.h"
 #include "AliTPCseed.h"
 #include "AliTPCReconstructor.h"
 #include "AliTPCClusterParam.h"
@@ -33,7 +34,7 @@
 #include "AliTPCTransform.h"
 #include "AliSplineFit.h"
 #include "AliCDBManager.h"
-
+#include "AliTPCcalibDButil.h"
 
 
 ClassImp(AliTPCseed)
@@ -579,27 +580,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];
@@ -1030,8 +1036,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){
@@ -1074,18 +1080,19 @@ Float_t  AliTPCseed::CookdEdxNorm(Double_t low, Double_t up, Int_t type, Int_t i
   UInt_t runNumber = 1;
   Float_t corrTimeGain = 1;
   AliTPCTransform * trans = AliTPCcalibDB::Instance()->GetTransform();
-  if (trans) {
+  const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
+  if (trans && recoParam->GetUseGainCorrectionTime()>0) {
     runNumber = trans->GetCurrentRunNumber();
-               //AliTPCcalibDB::Instance()->SetRun(runNumber);
+    //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);
+       corrTimeGain = AliTPCcalibDButil::EvalGraphConst(fitMIP, time);/*fitMIP->Eval(time);*/
       } else {
-       if (fitFPcosmic) corrTimeGain = fitFPcosmic->Eval(time); // This value describes the ratio FP-to-MIP, hardwired for the moment
+       if (fitFPcosmic) corrTimeGain = AliTPCcalibDButil::EvalGraphConst(fitFPcosmic, time);/*fitFPcosmic->Eval(time);*/ 
       }
     }
   }
@@ -1097,7 +1104,7 @@ Float_t  AliTPCseed::CookdEdxNorm(Double_t low, Double_t up, Int_t type, Int_t i
   return mean;
 }
 
-Float_t  AliTPCseed::CookdEdxAnalytical(Double_t low, Double_t up, Int_t type, Int_t i1, Int_t i2, Int_t returnVal){
+Float_t  AliTPCseed::CookdEdxAnalytical(Double_t low, Double_t up, Int_t type, Int_t i1, Int_t i2, Int_t returnVal, Int_t rowThres){
  
   //
   // calculates dedx using the cluster
@@ -1110,6 +1117,8 @@ Float_t  AliTPCseed::CookdEdxAnalytical(Double_t low, Double_t up, Int_t type, I
   // returnVal - 0 return mean
   //           - 1 return RMS
   //           - 2 return number of clusters
+  //
+  // rowThres  - number of rows before and after given pad row to check for clusters below threshold
   //           
   // normalization parametrization taken from AliTPCClusterParam
   //
@@ -1123,27 +1132,80 @@ Float_t  AliTPCseed::CookdEdxAnalytical(Double_t low, Double_t up, Int_t type, I
   Float_t amp[160];
   Int_t   indexes[160];
   Int_t   ncl=0;
+  Int_t   nclBelowThr = 0; // counts number of clusters below threshold
   //
   //
   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
   //   
   //
   //
   if (AliTPCcalibDB::Instance()->GetParameters()){
     gainGG= AliTPCcalibDB::Instance()->GetParameters()->GetGasGain()/20000;  //relative gas gain
   }
-
+  //
+  // extract time-dependent correction for pressure and temperature variations
+  //
+  UInt_t runNumber = 1;
+  Float_t corrTimeGain = 1;
+  TObjArray * timeGainSplines = 0x0;
+  TGraphErrors * grPadEqual = 0x0;
+  //
+  AliTPCTransform * trans = AliTPCcalibDB::Instance()->GetTransform();
+  const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
+  //
+  if (recoParam->GetNeighborRowsDedx() == 0) rowThres = 0;
+  //
+  if (trans) {
+      runNumber = trans->GetCurrentRunNumber();
+      //AliTPCcalibDB::Instance()->SetRun(runNumber);
+      timeGainSplines = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(runNumber);
+      if (timeGainSplines && recoParam->GetUseGainCorrectionTime()>0) {
+       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); */
+       }
+       //
+       if (type==1) grPadEqual = (TGraphErrors * ) timeGainSplines->FindObject("TGRAPHERRORS_MEANQMAX_PADREGIONGAIN_BEAM_ALL");
+       if (type==0) grPadEqual = (TGraphErrors * ) timeGainSplines->FindObject("TGRAPHERRORS_MEANQTOT_PADREGIONGAIN_BEAM_ALL");
+      }
+  }   
+  
+  const Float_t kClusterShapeCut = 1.5; // IMPPRTANT TO DO: move value to AliTPCRecoParam
   const Float_t ktany = TMath::Tan(TMath::DegToRad()*10);
   const Float_t kedgey =3.;
   //
   //
   for (Int_t irow=i1; irow<i2; irow++){
     AliTPCclusterMI* cluster = GetClusterPointer(irow);
+    if (!cluster && irow > 1 && irow < 157) {
+      Bool_t isClBefore = kFALSE;
+      Bool_t isClAfter  = kFALSE;
+      for(Int_t ithres = 1; ithres <= rowThres; ithres++) {
+       AliTPCclusterMI * clusterBefore = GetClusterPointer(irow - ithres);
+       if (clusterBefore) isClBefore = kTRUE;
+       AliTPCclusterMI * clusterAfter  = GetClusterPointer(irow + ithres);
+       if (clusterAfter) isClAfter = kTRUE;
+      }
+      if (isClBefore && isClAfter) nclBelowThr++;
+    }
     if (!cluster) continue;
+    //
+    //
     if (TMath::Abs(cluster->GetY())>cluster->GetX()*ktany-kedgey) continue; // edge cluster
+    //
+    AliTPCTrackerPoint * point = GetTrackPoint(irow);
+    if (point==0) continue;    
+    Float_t rsigmay = TMath::Sqrt(point->GetSigmaY());
+    if (rsigmay > kClusterShapeCut) continue;
+    //
+    if (cluster->IsUsed(11)) continue; // remove shared clusters for PbPb
+    //
     Float_t charge= (type%2)? cluster->GetMax():cluster->GetQ();
     Int_t  ipad= 0;
     if (irow>=row0) ipad=1;
@@ -1163,14 +1225,12 @@ Float_t  AliTPCseed::CookdEdxAnalytical(Double_t low, Double_t up, Int_t type, I
       } else {         // OROC
        factor = roc->GetValue(irow - row0, TMath::Nint(cluster->GetPad()));
       }
-      if (factor>0.5) gainPad=factor;
+      if (factor>0.3) gainPad=factor;
     }
-    
     //
     // Do position normalization - relative distance to 
     // center of pad- time bin
     
-    AliTPCTrackerPoint * point = GetTrackPoint(irow);
     Float_t              ty = TMath::Abs(point->GetAngleY());
     Float_t              tz = TMath::Abs(point->GetAngleZ()*TMath::Sqrt(1+ty*ty));
     Float_t yres0 = parcl->GetRMS0(0,ipad,0,0)/param->GetPadPitchWidth(cluster->GetDetector());
@@ -1202,53 +1262,49 @@ Float_t  AliTPCseed::CookdEdxAnalytical(Double_t low, Double_t up, Int_t type, I
       corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,2)*tz*tz);
       //
     }
-
+    //
+    // pad region equalization outside of cluster param
+    //
+    Float_t gainEqualPadRegion = 1;
+    if (grPadEqual) gainEqualPadRegion = grPadEqual->Eval(ipad);
     //
     amp[ncl]=charge;
     amp[ncl]/=gainGG;
     amp[ncl]/=gainPad;
     amp[ncl]/=corrPos;
+    amp[ncl]/=gainEqualPadRegion;
     //
     ncl++;
   }
 
   if (type>3) return ncl; 
   TMath::Sort(ncl,amp, indexes, kFALSE);
-
+  //
   if (ncl<10) return 0;
-  
+  //
+  Double_t * ampWithBelow = new Double_t[ncl + nclBelowThr];
+  for(Int_t iCl = 0; iCl < ncl + nclBelowThr; iCl++) {
+    if (iCl < nclBelowThr) {
+      ampWithBelow[iCl] = amp[indexes[0]];
+    } else {
+      ampWithBelow[iCl] = amp[indexes[iCl - nclBelowThr]];
+    }
+  }
+  //printf("DEBUG: %i shit %f", nclBelowThr, amp[indexes[0]]);
+  //
   Float_t suma=0;
   Float_t suma2=0;  
   Float_t sumn=0;
-  Int_t icl0=TMath::Nint(ncl*low);
-  Int_t icl1=TMath::Nint(ncl*up);
+  Int_t icl0=TMath::Nint((ncl + nclBelowThr)*low);
+  Int_t icl1=TMath::Nint((ncl + nclBelowThr)*up);
   for (Int_t icl=icl0; icl<icl1;icl++){
-    suma+=amp[indexes[icl]];
-    suma2+=amp[indexes[icl]]*amp[indexes[icl]];
+    suma+=ampWithBelow[icl];
+    suma2+=ampWithBelow[icl]*ampWithBelow[icl];
     sumn++;
   }
   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();
-  if (trans) {
-    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;
   //
@@ -1404,6 +1460,7 @@ Int_t  AliTPCseed::RefitTrack(AliTPCseed *seed, AliExternalTrackParam * parin, A
   //
   if (parin) (*parin)=paramIn;
   if (parout) (*parout)=paramOut;
+  delete track;
   return ncl;
 }