TPC/AliTPCcalibCalib.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17 ///////////////////////////////////////////////////////////////////////////////
  18 //                                                                           //
  19 //   Component for redoing the reconstruction from the clusters and tracks
  20 //
  21 //   The new calibration data used
  22 //
  23 //   In reality it overwrites the content of the ESD
  24 //
  25
  26 /*
  27
  28   gSystem->Load("libANALYSIS");
  29   gSystem->Load("libTPCcalib");
  30   //
  31   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
  32   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
  33   AliXRDPROOFtoolkit tool;
  34   TChain * chainCl = tool.MakeChain("calib.txt","Clusters",0,1);
  35   chainCl->Lookup();
  36   TChain * chainTr = tool.MakeChain("calib.txt","Tracks",0,1);
  37   chainTr->Lookup();
  38
  39
  40
  41 */
  42
  43
  44
  45 //  marian.ivanov@cern.ch
  46 //
  47 #include "AliTPCcalibCalib.h"
  48 #include "TSystem.h"
  49 #include "TFile.h"
  50 #include "TTreeStream.h"
  51 #include "AliLog.h"
  52 #include "TTimeStamp.h"
  53 #include "AliESDEvent.h"
  54 #include "AliESDfriend.h"
  55 #include "AliESDtrack.h"
  56 #include "AliTracker.h"
  57 #include "AliTPCClusterParam.h"
  58
  59 #include "AliTPCcalibDB.h"
  60 #include "AliTPCTransform.h"
  61 #include "AliTPCRecoParam.h"
  62 #include "AliTPCclusterMI.h"
  63 #include "AliTPCseed.h"
  64 #include "AliTPCPointCorrection.h"
  65
  66 ClassImp(AliTPCcalibCalib)
  67
  68 AliTPCcalibCalib::AliTPCcalibCalib():
  69 AliTPCcalibBase(),
  70   fApplyExBCorrection(1),      // apply ExB correction
  71   fApplyTOFCorrection(1),      // apply TOF correction
  72   fApplyPositionCorrection(1), // apply position correction
  73   fApplySectorAlignment(1),    // apply sector alignment
  74   fApplyRPhiCorrection(1),     // apply R-Phi correction
  75   fApplyRCorrection(1)         // apply Radial correction
  76
  77 {
  78   //
  79   // Constructor
  80   //
  81 }
  82
  83
  84 AliTPCcalibCalib::AliTPCcalibCalib(const Text_t *name, const Text_t *title)
  85   :AliTPCcalibBase(),
  86   fApplyExBCorrection(1),      // apply ExB correction
  87   fApplyTOFCorrection(1),      // apply TOF correction
  88   fApplyPositionCorrection(1), // apply position correction
  89   fApplySectorAlignment(1),    // apply sector alignment
  90   fApplyRPhiCorrection(1),     // apply R-Phi correction
  91   fApplyRCorrection(1)         // apply Radial correction
  92 {
  93   SetName(name);
  94   SetTitle(title);
  95 }
  96
  97
  98 AliTPCcalibCalib::AliTPCcalibCalib(const AliTPCcalibCalib&calib):
  99   AliTPCcalibBase(calib),
 100   fApplyExBCorrection(calib.GetApplyExBCorrection()),
 101   fApplyTOFCorrection(calib.GetApplyTOFCorrection()),
 102   fApplyPositionCorrection(calib.GetApplyPositionCorrection()),
 103   fApplySectorAlignment(calib.GetApplySectorAlignment()),
 104   fApplyRPhiCorrection(calib.GetApplyRPhiCorrection()),
 105   fApplyRCorrection(calib.GetApplyRCorrection())
 106
 107 {
 108   //
 109   // copy constructor
 110   //
 111 }
 112
 113 AliTPCcalibCalib &AliTPCcalibCalib::operator=(const AliTPCcalibCalib&calib){
 114   //
 115   //
 116   //
 117   ((AliTPCcalibBase *)this)->operator=(calib);
 118   return *this;
 119 }
 120
 121
 122 AliTPCcalibCalib::~AliTPCcalibCalib() {
 123   //
 124   // destructor
 125   //
 126 }
 127
 128
 129 void     AliTPCcalibCalib::Process(AliESDEvent *event){
 130   //
 131   //
 132   //
 133   if (!event) {
 134     return;
 135   }
 136   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
 137   if (!ESDfriend) {
 138    return;
 139   }
 140
 141   if (GetDebugLevel()>20) printf("Hallo world: Im here\n");
 142   Int_t ntracks=event->GetNumberOfTracks();
 143   //AliTPCcalibDB::Instance()->SetExBField(fMagF);
 144
 145   //
 146   //
 147   //
 148
 149   for (Int_t i=0;i<ntracks;++i) {
 150     AliESDtrack *track = event->GetTrack(i);
 151     AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
 152
 153     const AliExternalTrackParam * trackIn  = track->GetInnerParam();
 154     const AliExternalTrackParam * trackOut = track->GetOuterParam();
 155     AliExternalTrackParam * tpcOut   = (AliExternalTrackParam *)friendTrack->GetTPCOut();
 156     if (!trackIn) continue;
 157     if (!trackOut) continue;
 158     if (!tpcOut) continue;
 159     TObject *calibObject;
 160     AliTPCseed *seed = 0;
 161     for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
 162       if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
 163     }
 164     if (!seed) continue;
 165     RefitTrack(track, seed,event->GetMagneticField());
 166     (*tpcOut)=*(track->GetOuterParam());
 167   }
 168   return;
 169 }
 170
 171 Bool_t  AliTPCcalibCalib::RefitTrack(AliESDtrack * track, AliTPCseed *seed, Float_t magesd){
 172   //
 173   // Refit track
 174   // if magesd==0 forget the curvature
 175
 176   //
 177   // 0 - Setup transform object
 178   //
 179   AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
 180   AliTPCParam     *param     = AliTPCcalibDB::Instance()->GetParameters();
 181   transform->SetCurrentRun(fRun);
 182   transform->SetCurrentTimeStamp((UInt_t)fTime);
 183   if(!fApplyExBCorrection) { // disable ExB correction in transform
 184     if(transform->GetCurrentRecoParam())
 185       transform->GetCurrentRecoParamNonConst()->SetUseExBCorrection(0);
 186   }
 187   if(!fApplyTOFCorrection) { // disable TOF correction in transform
 188     if(transform->GetCurrentRecoParam())
 189       transform->GetCurrentRecoParamNonConst()->SetUseTOFCorrection(kFALSE);
 190   }
 191
 192   //
 193   // First apply calibration
 194   //
 195   AliTPCPointCorrection * corr =  AliTPCPointCorrection::Instance();
 196   for (Int_t irow=0;irow<159;irow++) {
 197     AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
 198     if (!cluster) continue;
 199     AliTPCclusterMI cl0(*cluster);
 200     Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
 201     Int_t i[1]={cluster->GetDetector()};
 202
 203     transform->Transform(x,i,0,1);
 204     //
 205     // get position correction
 206     //
 207     Int_t ipad=0;
 208     if (cluster->GetDetector()>35) ipad=1;
 209     Float_t dy =AliTPCClusterParam::SPosCorrection(0,ipad,cluster->GetPad(),cluster->GetTimeBin(),cluster->GetZ(),cluster->GetSigmaY2(),cluster->GetSigmaZ2(),cluster->GetMax());
 210     Float_t dz =AliTPCClusterParam::SPosCorrection(1,ipad,cluster->GetPad(),cluster->GetTimeBin(),cluster->GetZ(),cluster->GetSigmaY2(),cluster->GetSigmaZ2(),cluster->GetMax());
 211     // if(fApplyPositionCorrection) {
 212     //x[1]-=dy;
 213     //x[2]-=dz;
 214     // }
 215     //
 216     // Apply sector alignment
 217     //
 218     Double_t dxq = AliTPCPointCorrection::SGetCorrectionSector(0,cluster->GetDetector()%36,cluster->GetX(),
 219                                                                cluster->GetY(),cluster->GetZ());
 220     Double_t dyq = AliTPCPointCorrection::SGetCorrectionSector(1,cluster->GetDetector()%36,cluster->GetX(),
 221                                                                cluster->GetY(),cluster->GetZ());
 222     Double_t dzq = AliTPCPointCorrection::SGetCorrectionSector(2,cluster->GetDetector()%36,cluster->GetX(),
 223                                                                cluster->GetY(),cluster->GetZ());
 224     if (fApplySectorAlignment){
 225       x[0]-=dxq;
 226       x[1]-=dyq;
 227       x[2]-=dzq;
 228     }
 229     //
 230     // Apply r-phi correction  - To be done on track level- knowing the track angle !!!
 231     //
 232     Double_t corrclY =
 233       corr->RPhiCOGCorrection(cluster->GetDetector(),cluster->GetRow(), cluster->GetPad(),
 234                                   cluster->GetY(),cluster->GetY(), cluster->GetZ(), 0., cluster->GetMax(),2.5);
 235     // R correction
 236     Double_t corrR   = corr->CorrectionOutR0(kFALSE,kFALSE,cluster->GetX(),cluster->GetY(),cluster->GetZ(),cluster->GetDetector());
 237
 238     if (fApplyRPhiCorrection){
 239       if (cluster->GetY()>0) x[1]+=corrclY;  // rphi correction
 240       if (cluster->GetY()<0) x[1]-=corrclY;  // rphi correction
 241     }
 242     if (fApplyRCorrection){
 243       x[0]+=corrR;                           // radial correction
 244     }
 245
 246     //
 247     //
 248     //
 249     cluster->SetX(x[0]);
 250     cluster->SetY(x[1]);
 251     cluster->SetZ(x[2]);
 252
 253     //
 254     // Apply alignemnt
 255     //
 256
 257 //     if (!param->IsGeoRead()) param->ReadGeoMatrices();
 258 //     TGeoHMatrix  *mat = param->GetClusterMatrix(cluster->GetDetector());
 259 //     //TGeoHMatrix  mat;
 260 //     Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
 261 //     Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
 262 //     if (mat) mat->LocalToMaster(pos,posC);
 263 //     else{
 264 //       // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
 265 //     }
 266 //     cluster->SetX(posC[0]);
 267 //     cluster->SetY(posC[1]);
 268 //     cluster->SetZ(posC[2]);
 269
 270
 271
 272
 273     if (fStreamLevel>2){
 274       TTreeSRedirector *cstream = GetDebugStreamer();
 275       if (cstream){
 276         (*cstream)<<"Clusters"<<
 277           "run="<<fRun<<              //  run number
 278           "event="<<fEvent<<          //  event number
 279           "time="<<fTime<<            //  time stamp of event
 280           "trigger="<<fTrigger<<      //  trigger
 281           "triggerClass="<<&fTriggerClass<<      //  trigger
 282           "mag="<<fMagF<<             //  magnetic field
 283           "cl0.="<<&cl0<<
 284           "cl.="<<cluster<<
 285           "cy="<<dy<<
 286           "cz="<<dz<<
 287           "cY="<<corrclY<<
 288           "cR="<<corrR<<
 289           "dxq="<<dxq<<
 290           "dyq="<<dyq<<
 291           "dzq="<<dzq<<
 292           "\n";
 293       }
 294     }
 295   }
 296   //
 297   //
 298   //
 299   Int_t ncl = seed->GetNumberOfClusters();
 300   Double_t covar[15];
 301   for (Int_t i=0;i<15;i++) covar[i]=0;
 302   covar[0]=10.*10.;
 303   covar[2]=10.*10.;
 304   covar[5]=10.*10./(64.*64.);
 305   covar[9]=10.*10./(64.*64.);
 306   covar[14]=0.3*0.3;
 307   if (TMath::Abs(magesd)<0.05) {
 308      covar[14]=0.025*0.025;
 309   }
 310   //
 311   // And now do refit
 312   //
 313   AliExternalTrackParam * trackInOld  = (AliExternalTrackParam*)track->GetInnerParam();
 314   AliExternalTrackParam * trackOutOld = (AliExternalTrackParam*)track->GetOuterParam();
 315
 316
 317   AliExternalTrackParam trackIn  = *trackOutOld;
 318   if (TMath::Abs(magesd)<0.05) {
 319     ((Double_t&)(trackIn.GetParameter()[4]))=0.000000001;
 320   }
 321   trackIn.ResetCovariance(20.);
 322   trackIn.AddCovariance(covar);
 323   Double_t xyz[3];
 324   Int_t nclIn=0,nclOut=0;
 325   //
 326   // Refit in
 327   //
 328
 329   for (Int_t irow=159; irow>0; irow--){
 330     AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
 331     if (!cl) continue;
 332     if (cl->GetX()<80) continue;
 333     Int_t sector = cl->GetDetector();
 334     Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackIn.GetAlpha();
 335     if (TMath::Abs(dalpha)>0.01){
 336       if (!trackIn.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
 337     }
 338     Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
 339     Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
 340     AliTPCseed::GetError(cl, &trackIn,cov[0],cov[2]);
 341     cov[0]*=cov[0];
 342     cov[2]*=cov[2];
 343     trackIn.GetXYZ(xyz);
 344     Double_t bz = AliTracker::GetBz(xyz);
 345
 346     if (!trackIn.PropagateTo(r[0],bz)) continue;
 347     if (RejectCluster(cl,&trackIn)) continue;
 348     nclIn++;
 349     trackIn.Update(&r[1],cov);
 350   }
 351   //
 352   AliExternalTrackParam trackOut = trackIn;
 353   if (TMath::Abs(magesd)<0.05) {
 354     ((Double_t&)(trackOut.GetParameter()[4]))=0.000000001;
 355   }
 356   trackOut.ResetCovariance(20.);
 357   trackOut.AddCovariance(covar);
 358   //
 359   // Refit out
 360   //
 361   //Bool_t lastEdge=kFALSE;
 362   for (Int_t irow=0; irow<160; irow++){
 363     AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
 364     if (!cl) continue;
 365     if (cl->GetX()<80) continue;
 366     Int_t sector = cl->GetDetector();
 367     Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackOut.GetAlpha();
 368
 369     if (TMath::Abs(dalpha)>0.01){
 370       if (!trackOut.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
 371     }
 372     Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
 373
 374     Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
 375     AliTPCseed::GetError(cl, &trackOut,cov[0],cov[2]);
 376     cov[0]*=cov[0];
 377     cov[2]*=cov[2];
 378     trackOut.GetXYZ(xyz);
 379     Double_t bz = AliTracker::GetBz(xyz);
 380     if (!trackOut.PropagateTo(r[0],bz)) continue;
 381     if (RejectCluster(cl,&trackOut)) continue;
 382     nclOut++;
 383     trackOut.Update(&r[1],cov);
 384     //if (cl->GetType()<0) lastEdge=kTRUE;
 385     //if (cl->GetType()>=0) lastEdge=kFALSE;
 386   }
 387   //
 388   //
 389   //
 390   nclIn=0;
 391   trackIn  = trackOut;
 392   trackIn.ResetCovariance(10.);
 393   //
 394   // Refit in one more time
 395   //
 396   for (Int_t irow=159; irow>0; irow--){
 397     AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
 398     if (!cl) continue;
 399     if (cl->GetX()<80) continue;
 400     Int_t sector = cl->GetDetector();
 401     Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackIn.GetAlpha();
 402     if (TMath::Abs(dalpha)>0.01){
 403       if (!trackIn.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
 404     }
 405     Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
 406     Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
 407     AliTPCseed::GetError(cl, &trackIn,cov[0],cov[2]);
 408     cov[0]*=cov[0];
 409     cov[2]*=cov[2];
 410     trackIn.GetXYZ(xyz);
 411     Double_t bz = AliTracker::GetBz(xyz);
 412
 413     if (!trackIn.PropagateTo(r[0],bz)) continue;
 414     if (RejectCluster(cl,&trackIn)) continue;
 415     nclIn++;
 416     trackIn.Update(&r[1],cov);
 417   }
 418
 419
 420   trackIn.Rotate(trackInOld->GetAlpha());
 421   trackOut.Rotate(trackOutOld->GetAlpha());
 422   //
 423   trackInOld->GetXYZ(xyz);
 424   Double_t bz = AliTracker::GetBz(xyz);
 425   trackIn.PropagateTo(trackInOld->GetX(),bz);
 426   //
 427   trackOutOld->GetXYZ(xyz);
 428   bz = AliTracker::GetBz(xyz);
 429   trackOut.PropagateTo(trackOutOld->GetX(),bz);
 430
 431
 432   if (fStreamLevel>0){
 433     TTreeSRedirector *cstream = GetDebugStreamer();
 434     if (cstream){
 435       (*cstream)<<"Tracks"<<
 436         "run="<<fRun<<              //  run number
 437         "event="<<fEvent<<          //  event number
 438         "time="<<fTime<<            //  time stamp of event
 439         "trigger="<<fTrigger<<      //  trigger
 440         "triggerClass="<<&fTriggerClass<<      //  trigger
 441         "mag="<<fMagF<<             //  magnetic field
 442         "nclIn="<<nclIn<<
 443         "nclOut="<<nclOut<<
 444         "ncl="<<ncl<<
 445         "TrIn0.="<<trackInOld<<
 446         "TrOut0.="<<trackOutOld<<
 447         "TrIn1.="<<&trackIn<<
 448         "TrOut1.="<<&trackOut<<
 449         "\n";
 450     }
 451   }
 452   //
 453   // And now rewrite ESDtrack and TPC seed
 454   //
 455
 456   (*trackInOld)  = trackIn;
 457   (*trackOutOld) = trackOut;
 458   AliExternalTrackParam *t = &trackIn;
 459   track->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
 460   seed->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
 461   seed->SetNumberOfClusters((nclIn+nclOut)/2);
 462   return kTRUE;
 463 }
 464
 465
 466
 467 Bool_t AliTPCcalibCalib::RejectCluster(AliTPCclusterMI* cl, AliExternalTrackParam * param){
 468   //
 469   // check the acceptance of cluster
 470   // Cut on edge effects
 471   //
 472   Float_t kEdgeCut=2.5;
 473   Float_t kSigmaCut=6;
 474
 475   Bool_t isReject = kFALSE;
 476   Float_t edgeY = cl->GetX()*TMath::Tan(TMath::Pi()/18);
 477   Float_t dist  = edgeY - TMath::Abs(cl->GetY());
 478   if (param)  dist  = TMath::Abs(edgeY - TMath::Abs(param->GetY()));
 479   if (dist<kEdgeCut) isReject=kTRUE;
 480
 481   Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
 482   AliTPCseed::GetError(cl, param,cov[0],cov[2]);
 483   if (param->GetSigmaY2()<0 || param->GetSigmaZ2()<0){
 484     AliError("Wrong parameters");
 485     return kFALSE;
 486   }
 487   Double_t py = (cl->GetY()-param->GetY())/TMath::Sqrt(cov[0]*cov[0]+param->GetSigmaY2());
 488   Double_t pz = (cl->GetZ()-param->GetZ())/TMath::Sqrt(cov[2]*cov[2]+param->GetSigmaZ2());
 489   //
 490   if ((py*py+pz*pz)>kSigmaCut*kSigmaCut) isReject=kTRUE;
 491
 492   return isReject;
 493 }
 494
 495
 496