]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TPC/AliTPCcalibCalib.cxx
Adding Extraction of OCDB entries - spline fits by default
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibCalib.cxx
index d0c60d2682db9e4d4bdbc1093a7d5dccd6de2202..b7ea0f8d72442dec22d4dd739aca956b186441b3 100644 (file)
 //   In reality it overwrites the content of the ESD 
 //    
 
+/*
+
+  gSystem->Load("libANALYSIS");
+  gSystem->Load("libTPCcalib");
+  //
+  gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
+  gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
+  AliXRDPROOFtoolkit tool; 
+  TChain * chainCl = tool.MakeChain("calib.txt","Clusters",0,1);
+  chainCl->Lookup();
+  TChain * chainTr = tool.MakeChain("calib.txt","Tracks",0,1);
+  chainTr->Lookup();
+
+
+
+*/
+
+
+
 //  marian.ivanov@cern.ch
 // 
 #include "AliTPCcalibCalib.h"
 #include "AliESDfriend.h"
 #include "AliESDtrack.h"
 #include "AliTracker.h"
+#include "AliTPCClusterParam.h"
 
 #include "AliTPCcalibDB.h"
 #include "AliTPCTransform.h"
 #include "AliTPCclusterMI.h"
 #include "AliTPCseed.h"
+#include "AliTPCPointCorrection.h"
 
 ClassImp(AliTPCcalibCalib)
 
@@ -98,6 +119,8 @@ void     AliTPCcalibCalib::Process(AliESDEvent *event){
   
   if (GetDebugLevel()>20) printf("Hallo world: Im here\n");
   Int_t ntracks=event->GetNumberOfTracks();   
+  //AliTPCcalibDB::Instance()->SetExBField(fMagF);
+
   //
   //
   //
@@ -126,11 +149,16 @@ Bool_t  AliTPCcalibCalib::RefitTrack(AliESDtrack * track, AliTPCseed *seed){
   // Refit track
   //
 
+  //
+  // 0 - Setup transform object
+  //
+  AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
+  transform->SetCurrentRun(fRun);
+  transform->SetCurrentTimeStamp((UInt_t)fTime);
   //
   // First apply calibration
   //
-
-  AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform();
+  AliTPCPointCorrection * corr =  AliTPCPointCorrection::Instance();
   for (Int_t irow=0;irow<159;irow++) {
     AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
     if (!cluster) continue; 
@@ -138,9 +166,46 @@ Bool_t  AliTPCcalibCalib::RefitTrack(AliESDtrack * track, AliTPCseed *seed){
     Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
     Int_t i[1]={cluster->GetDetector()};
     transform->Transform(x,i,0,1);
-    if (cluster->GetDetector()%36>17){
-      x[1]*=-1;
+    //
+    // get position correction
+    //
+    Int_t ipad=0;
+    if (cluster->GetDetector()>35) ipad=1;
+    Float_t dy =AliTPCClusterParam::SPosCorrection(0,ipad,cluster->GetPad(),cluster->GetTimeBin(),cluster->GetZ(),cluster->GetSigmaY2(),cluster->GetSigmaZ2(),cluster->GetMax());
+    Float_t dz =AliTPCClusterParam::SPosCorrection(1,ipad,cluster->GetPad(),cluster->GetTimeBin(),cluster->GetZ(),cluster->GetSigmaY2(),cluster->GetSigmaZ2(),cluster->GetMax());
+    //x[1]-=dy;
+    //x[2]-=dz;
+    //
+    // Apply sector alignment
+    //
+    Double_t dxq = AliTPCPointCorrection::SGetCorrectionSector(0,cluster->GetDetector()%36,cluster->GetX(),
+                                                              cluster->GetY(),cluster->GetZ());
+    Double_t dyq = AliTPCPointCorrection::SGetCorrectionSector(1,cluster->GetDetector()%36,cluster->GetX(),
+                                                              cluster->GetY(),cluster->GetZ());
+    Double_t dzq = AliTPCPointCorrection::SGetCorrectionSector(2,cluster->GetDetector()%36,cluster->GetX(),
+                                                              cluster->GetY(),cluster->GetZ());
+    if (kTRUE){
+      x[0]-=dxq;
+      x[1]-=dyq;
+      x[2]-=dzq;
+    }
+    //
+    // Apply r-phi correction  - To be done on track level- knowing the track angle !!!
+    //
+    Double_t corrclY =  
+      corr->RPhiCOGCorrection(cluster->GetDetector(),cluster->GetRow(), cluster->GetPad(),
+                                 cluster->GetY(),cluster->GetY(), cluster->GetZ(), 0., cluster->GetMax(),2.5);
+    // R correction
+    Double_t corrR   = corr->CorrectionOutR0(kFALSE,kFALSE,cluster->GetX(),cluster->GetY(),cluster->GetZ(),cluster->GetDetector());
+
+    if (kTRUE){
+      if (cluster->GetY()>0) x[1]+=corrclY;  // rphi correction
+      if (cluster->GetY()<0) x[1]-=corrclY;  // rphi correction
+      x[0]+=corrR;                           // radial correction
     }
+
+    //
+    //
     //
     cluster->SetX(x[0]);
     cluster->SetY(x[1]);
@@ -149,13 +214,37 @@ Bool_t  AliTPCcalibCalib::RefitTrack(AliESDtrack * track, AliTPCseed *seed){
       TTreeSRedirector *cstream = GetDebugStreamer();
       if (cstream){
        (*cstream)<<"Clusters"<<
+         "run="<<fRun<<              //  run number
+         "event="<<fEvent<<          //  event number
+         "time="<<fTime<<            //  time stamp of event
+         "trigger="<<fTrigger<<      //  trigger
+         "triggerClass="<<&fTriggerClass<<      //  trigger      
+         "mag="<<fMagF<<             //  magnetic field
          "cl0.="<<&cl0<<
          "cl.="<<cluster<<
+         "cy="<<dy<<
+         "cz="<<dz<<
+         "cY="<<corrclY<<
+         "cR="<<corrR<<
+         "dxq="<<dxq<<
+         "dyq="<<dyq<<
+         "dzq="<<dzq<<
          "\n";
       }
     }
   }
-  Int_t ncl = seed->GetNumberOfClusters();  
+  //
+  //
+  //
+  Int_t ncl = seed->GetNumberOfClusters();
+  Double_t covar[15];
+  for (Int_t i=0;i<15;i++) covar[i]=0;
+  covar[0]=10.*10.;
+  covar[2]=10.*10.;
+  covar[5]=10.*10./(64.*64.);
+  covar[9]=10.*10./(64.*64.);
+  covar[14]=1*1;
+
   // 
   // And now do refit
   //
@@ -163,14 +252,43 @@ Bool_t  AliTPCcalibCalib::RefitTrack(AliESDtrack * track, AliTPCseed *seed){
   AliExternalTrackParam * trackOutOld = (AliExternalTrackParam*)track->GetOuterParam();
 
   AliExternalTrackParam trackIn  = *trackOutOld;
-  AliExternalTrackParam trackOut = *trackInOld;
-  trackIn.ResetCovariance(10.);
-  trackOut.ResetCovariance(10.);
+  trackIn.ResetCovariance(200.);
+  trackIn.AddCovariance(covar);
   Double_t xyz[3];
   Int_t nclIn=0,nclOut=0;
   //
+  // Refit in
+  //
+
+  for (Int_t irow=159; irow>0; irow--){
+    AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
+    if (!cl) continue;
+    if (cl->GetX()<80) continue;
+    Int_t sector = cl->GetDetector();
+    Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackIn.GetAlpha();
+    if (TMath::Abs(dalpha)>0.01)
+      trackIn.Rotate(TMath::DegToRad()*(sector%18*20.+10.));
+    Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
+    Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
+    AliTPCseed::GetError(cl, &trackIn,cov[0],cov[2]);
+    cov[0]*=cov[0];
+    cov[2]*=cov[2];
+    trackIn.GetXYZ(xyz);
+    Double_t bz = AliTracker::GetBz(xyz);
+
+    if (!trackIn.PropagateTo(r[0],bz)) continue;
+    if (RejectCluster(cl,&trackIn)) continue;
+    nclIn++;
+    trackIn.Update(&r[1],cov);    
+  }
+  //
+  AliExternalTrackParam trackOut = trackIn;
+  trackOut.ResetCovariance(200.);
+  trackOut.AddCovariance(covar);
+  //
   // Refit out
   //
+  //Bool_t lastEdge=kFALSE;
   for (Int_t irow=0; irow<160; irow++){
     AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
     if (!cl) continue;
@@ -180,19 +298,30 @@ Bool_t  AliTPCcalibCalib::RefitTrack(AliESDtrack * track, AliTPCseed *seed){
 
     if (TMath::Abs(dalpha)>0.01)
       trackOut.Rotate(TMath::DegToRad()*(sector%18*20.+10.));
-    //if (RejectCluster(cl,&trackOut)) continue;
     Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
-    Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
+
+    Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation    
+    AliTPCseed::GetError(cl, &trackOut,cov[0],cov[2]);
+    cov[0]*=cov[0];
+    cov[2]*=cov[2];
     trackOut.GetXYZ(xyz);
     Double_t bz = AliTracker::GetBz(xyz);
-    if (trackOut.PropagateTo(r[0],bz)) nclOut++;
-    trackOut.Update(&r[1],cov);    
-    
+    if (!trackOut.PropagateTo(r[0],bz)) continue;
+    if (RejectCluster(cl,&trackOut)) continue;
+    nclOut++;
+    trackOut.Update(&r[1],cov);        
+    //if (cl->GetType()<0) lastEdge=kTRUE;
+    //if (cl->GetType()>=0) lastEdge=kFALSE;    
   }
   //
-  // Refit in
   //
-
+  //
+  nclIn=0;
+  trackIn  = trackOut;
+  trackIn.ResetCovariance(10.);
+  //
+  // Refit in one more time
+  //
   for (Int_t irow=159; irow>0; irow--){
     AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
     if (!cl) continue;
@@ -201,20 +330,43 @@ Bool_t  AliTPCcalibCalib::RefitTrack(AliESDtrack * track, AliTPCseed *seed){
     Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackIn.GetAlpha();
     if (TMath::Abs(dalpha)>0.01)
       trackIn.Rotate(TMath::DegToRad()*(sector%18*20.+10.));
-    //if (RejectCluster(cl,&trackIn)) continue;
     Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
     Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
-    trackOut.GetXYZ(xyz);
+    AliTPCseed::GetError(cl, &trackIn,cov[0],cov[2]);
+    cov[0]*=cov[0];
+    cov[2]*=cov[2];
+    trackIn.GetXYZ(xyz);
     Double_t bz = AliTracker::GetBz(xyz);
 
-    if (trackIn.PropagateTo(r[0],bz)) nclIn++;
+    if (!trackIn.PropagateTo(r[0],bz)) continue;
+    if (RejectCluster(cl,&trackIn)) continue;
+    nclIn++;
     trackIn.Update(&r[1],cov);    
   }
 
+
+  trackIn.Rotate(trackInOld->GetAlpha());
+  trackOut.Rotate(trackOutOld->GetAlpha());
+  //
+  trackInOld->GetXYZ(xyz);
+  Double_t bz = AliTracker::GetBz(xyz);
+  trackIn.PropagateTo(trackInOld->GetX(),bz);
+  //
+  trackOutOld->GetXYZ(xyz);
+  bz = AliTracker::GetBz(xyz);  
+  trackOut.PropagateTo(trackOutOld->GetX(),bz);
+  
+
   if (fStreamLevel>0){
     TTreeSRedirector *cstream = GetDebugStreamer();
     if (cstream){
       (*cstream)<<"Tracks"<<
+       "run="<<fRun<<              //  run number
+       "event="<<fEvent<<          //  event number
+       "time="<<fTime<<            //  time stamp of event
+       "trigger="<<fTrigger<<      //  trigger
+       "triggerClass="<<&fTriggerClass<<      //  trigger
+       "mag="<<fMagF<<             //  magnetic field
        "nclIn="<<nclIn<<
        "nclOut="<<nclOut<<
        "ncl="<<ncl<<
@@ -225,6 +377,17 @@ Bool_t  AliTPCcalibCalib::RefitTrack(AliESDtrack * track, AliTPCseed *seed){
        "\n";
     }
   }
+  //
+  // And now rewrite ESDtrack and TPC seed
+  //
+  (*trackInOld)  = trackIn;
+  (*trackOutOld) = trackOut;
+  AliExternalTrackParam *t = &trackIn;
+  track->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
+  seed->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
+  seed->SetNumberOfClusters((nclIn+nclOut)/2);
+  return kTRUE;
 }
 
 
@@ -234,13 +397,24 @@ Bool_t AliTPCcalibCalib::RejectCluster(AliTPCclusterMI* cl, AliExternalTrackPara
   // check the acceptance of cluster
   // Cut on edge effects
   //
+  Float_t kEdgeCut=2.5;
+  Float_t kSigmaCut=6;
+
   Bool_t isReject = kFALSE;
   Float_t edgeY = cl->GetX()*TMath::Tan(TMath::Pi()/18);
   Float_t dist  = edgeY - TMath::Abs(cl->GetY());
   if (param)  dist  = TMath::Abs(edgeY - TMath::Abs(param->GetY()));
-  if (dist<3) isReject=kTRUE;
-  if (cl->GetType()<0) isReject=kTRUE;
+  if (dist<kEdgeCut) isReject=kTRUE;
+
+  Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation    
+  AliTPCseed::GetError(cl, param,cov[0],cov[2]);
+  Double_t py = (cl->GetY()-param->GetY())/TMath::Sqrt(cov[0]*cov[0]+param->GetSigmaY2());
+  Double_t pz = (cl->GetZ()-param->GetZ())/TMath::Sqrt(cov[2]*cov[2]+param->GetSigmaZ2());
+  //
+  if ((py*py+pz*pz)>kSigmaCut*kSigmaCut) isReject=kTRUE;
+  
   return isReject;
 }
 
 
+