]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TPC/AliTPCseed.cxx
DQM histograms fixed
[u/mrichter/AliRoot.git] / TPC / AliTPCseed.cxx
index 61543888d4d178147f783c1983e6739061601726..45336b107687a34f26cdab30470f68799f678f47 100644 (file)
 
 
 //-----------------------------------------------------------------
+//
 //           Implementation of the TPC seed class
 //        This class is used by the AliTPCtrackerMI class
 //      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 +35,7 @@
 #include "AliTPCTransform.h"
 #include "AliSplineFit.h"
 #include "AliCDBManager.h"
-
+#include "AliTPCcalibDButil.h"
 
 
 ClassImp(AliTPCseed)
@@ -1035,8 +1037,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){
@@ -1089,9 +1091,9 @@ Float_t  AliTPCseed::CookdEdxNorm(Double_t low, Double_t up, Int_t type, Int_t i
       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);*/ 
       }
     }
   }
@@ -1103,7 +1105,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, Int_t mode){
  
   //
   // calculates dedx using the cluster
@@ -1113,9 +1115,16 @@ Float_t  AliTPCseed::CookdEdxAnalytical(Double_t low, Double_t up, Int_t type, I
   // i1-i2  -  the pad-row range used for calculation
   //           
   // posNorm   - usage of pos normalization 
-  // returnVal - 0 return mean
-  //           - 1 return RMS
-  //           - 2 return number of clusters
+  // returnVal - 0  return mean
+  //           - 1  return RMS
+  //           - 2  return number of clusters
+  //           - 3  ratio
+  //           - 4  mean upper half
+  //           - 5  mean  - lower half
+  //           - 6  third moment
+  // mode      - 0 - linear
+  //           - 1 - logatithmic
+  // rowThres  - number of rows before and after given pad row to check for clusters below threshold
   //           
   // normalization parametrization taken from AliTPCClusterParam
   //
@@ -1129,27 +1138,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
-
   //   
   //
   //
   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;
@@ -1169,14 +1231,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());
@@ -1208,59 +1268,116 @@ 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; 
+  if (type==2) 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 suma3=0;  
+  Float_t sumaS=0;  
   Float_t sumn=0;
-  Int_t icl0=TMath::Nint(ncl*low);
-  Int_t icl1=TMath::Nint(ncl*up);
-  for (Int_t icl=icl0; icl<icl1;icl++){
-    suma+=amp[indexes[icl]];
-    suma2+=amp[indexes[icl]]*amp[indexes[icl]];
-    sumn++;
-  }
-  Float_t mean =suma/sumn;
-  Float_t rms  =TMath::Sqrt(TMath::Abs(suma2/sumn-mean*mean));
+  // upper,and lower part statistic
+  Float_t sumL=0, sumL2=0, sumLN=0;
+  Float_t sumD=0, sumD2=0, sumDN=0;
+
+  Int_t icl0=TMath::Nint((ncl + nclBelowThr)*low);
+  Int_t icl1=TMath::Nint((ncl + nclBelowThr)*up);
+  Int_t iclm=TMath::Nint((ncl + nclBelowThr)*(low +(up+low)*0.5));
   //
-  // 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
+  for (Int_t icl=icl0; icl<icl1;icl++){
+    if (ampWithBelow[icl]<0.1) continue;
+    Double_t camp=ampWithBelow[icl]/corrTimeGain;
+    if (mode==1) camp= TMath::Log(camp);
+    if (icl<icl1){
+      suma+=camp;
+      suma2+=camp*camp;
+      suma3+=camp*camp*camp;
+      sumaS+=TMath::Power(TMath::Abs(camp),1./3.);
+      sumn++;
+    }
+    if (icl>iclm){
+      sumL+=camp;
+      sumL2+=camp*camp;
+      sumLN++;
       }
+    if (icl<=iclm){
+      sumD+=camp;
+      sumD2+=camp*camp;
+      sumDN++;
     }
   }
-  mean /= corrTimeGain;
-  rms /= corrTimeGain;
+  //
+  Float_t mean = 0;
+  Float_t meanL = 0;  
+  Float_t meanD = 0;           // lower half mean
+  if (sumn > 1e-30)   mean =suma/sumn;
+  if (sumLN > 1e-30)  meanL =sumL/sumLN;
+  if (sumDN > 1e-30)  meanD =(sumD/sumDN);
+  /*
+  Float_t mean =suma/sumn;
+  Float_t meanL = sumL/sumLN;  
+  Float_t meanD =(sumD/sumDN);           // lower half mean
+  */
+
+  Float_t rms = 0;
+  Float_t mean2=0;
+  Float_t mean3=0;
+  Float_t meanS=0;
+
+  if(sumn>0){
+    rms = TMath::Sqrt(TMath::Abs(suma2/sumn-mean*mean));
+    mean2=suma2/sumn;
+    mean3=suma3/sumn;
+    meanS=sumaS/sumn;
+  }
+
+  if (mean2>0) mean2=TMath::Power(TMath::Abs(mean2),1./2.);
+  if (mean3>0) mean3=TMath::Power(TMath::Abs(mean3),1./3.);
+  if (meanS>0) meanS=TMath::Power(TMath::Abs(meanS),3.);
+  //
+  if (mode==1) mean=TMath::Exp(mean);
+  if (mode==1) meanL=TMath::Exp(meanL);  // upper truncation
+  if (mode==1) meanD=TMath::Exp(meanD);  // lower truncation
+  //
+  delete [] ampWithBelow;
+  
+
   //
   if (returnVal==1) return rms;
   if (returnVal==2) return ncl;
+  if (returnVal==3) return Double_t(nclBelowThr)/Double_t(nclBelowThr+ncl);
+  if (returnVal==4) return meanL;
+  if (returnVal==5) return meanD;
+  if (returnVal==6) return mean2;
+  if (returnVal==7) return mean3;
+  if (returnVal==8) return meanS;
   return mean;
 }
 
@@ -1411,6 +1528,7 @@ Int_t  AliTPCseed::RefitTrack(AliTPCseed *seed, AliExternalTrackParam * parin, A
   //
   if (parin) (*parin)=paramIn;
   if (parout) (*parout)=paramOut;
+  delete track;
   return ncl;
 }