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 #include "AliTPCParam.h"
  59
  60 #include "AliTPCcalibDB.h"
  61 #include "AliTPCTransform.h"
  62 #include "AliTPCRecoParam.h"
  63 #include "AliTPCclusterMI.h"
  64 #include "AliTPCseed.h"
  65 #include "AliTPCPointCorrection.h"
  66 #include <TGeoManager.h>
  67 #include <TGeoPhysicalNode.h>
  68 #include "TDatabasePDG.h"
  69 ClassImp(AliTPCcalibCalib)
  70
  71 AliTPCcalibCalib::AliTPCcalibCalib():
  72 AliTPCcalibBase(),
  73   fApplyExBCorrection(1),      // apply ExB correction
  74   fApplyTOFCorrection(1),      // apply TOF correction
  75   fApplyPositionCorrection(1), // apply position correction
  76   fApplySectorAlignment(1),    // apply sector alignment
  77   fApplyRPhiCorrection(1),     // apply R-Phi correction
  78   fApplyRCorrection(1)         // apply Radial correction
  79
  80 {
  81   //
  82   // Constructor
  83   //
  84 }
  85
  86
  87 AliTPCcalibCalib::AliTPCcalibCalib(const Text_t *name, const Text_t *title)
  88   :AliTPCcalibBase(),
  89   fApplyExBCorrection(1),      // apply ExB correction
  90   fApplyTOFCorrection(1),      // apply TOF correction
  91   fApplyPositionCorrection(1), // apply position correction
  92   fApplySectorAlignment(1),    // apply sector alignment
  93   fApplyRPhiCorrection(1),     // apply R-Phi correction
  94   fApplyRCorrection(1)         // apply Radial correction
  95 {
  96   SetName(name);
  97   SetTitle(title);
  98 }
  99
 100
 101 AliTPCcalibCalib::AliTPCcalibCalib(const AliTPCcalibCalib&calib):
 102   AliTPCcalibBase(calib),
 103   fApplyExBCorrection(calib.GetApplyExBCorrection()),
 104   fApplyTOFCorrection(calib.GetApplyTOFCorrection()),
 105   fApplyPositionCorrection(calib.GetApplyPositionCorrection()),
 106   fApplySectorAlignment(calib.GetApplySectorAlignment()),
 107   fApplyRPhiCorrection(calib.GetApplyRPhiCorrection()),
 108   fApplyRCorrection(calib.GetApplyRCorrection())
 109
 110 {
 111   //
 112   // copy constructor
 113   //
 114 }
 115
 116 AliTPCcalibCalib &AliTPCcalibCalib::operator=(const AliTPCcalibCalib&calib){
 117   //
 118   //
 119   //
 120   ((AliTPCcalibBase *)this)->operator=(calib);
 121   return *this;
 122 }
 123
 124
 125 AliTPCcalibCalib::~AliTPCcalibCalib() {
 126   //
 127   // destructor
 128   //
 129 }
 130
 131
 132 void     AliTPCcalibCalib::Process(AliESDEvent *event){
 133   //
 134   //
 135   //
 136   if (!event) {
 137     return;
 138   }
 139   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
 140   if (!ESDfriend) {
 141    return;
 142   }
 143
 144   if (GetDebugLevel()>20) printf("Hallo world: Im here\n");
 145   Int_t ntracks=event->GetNumberOfTracks();
 146   //AliTPCcalibDB::Instance()->SetExBField(fMagF);
 147
 148   //
 149   //
 150   //
 151
 152   for (Int_t i=0;i<ntracks;++i) {
 153     AliESDtrack *track = event->GetTrack(i);
 154     AliESDfriendTrack *friendTrack = (AliESDfriendTrack*) ESDfriend->GetTrack(i);
 155     if (!friendTrack) continue;
 156     track->SetFriendTrack(friendTrack);
 157     const AliExternalTrackParam * trackIn  = track->GetInnerParam();
 158     const AliExternalTrackParam * trackOut = track->GetOuterParam();
 159     AliExternalTrackParam * tpcOut   = (AliExternalTrackParam *)friendTrack->GetTPCOut();
 160     if (!trackIn) continue;
 161     if (!trackOut) continue;
 162     if (!tpcOut) continue;
 163     TObject *calibObject;
 164     AliTPCseed *seed = 0;
 165     for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
 166       if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
 167     }
 168     if (!seed) continue;
 169     RefitTrack(track, seed, event->GetMagneticField());
 170     (*tpcOut)=*(track->GetOuterParam());
 171   }
 172   return;
 173 }
 174
 175 Bool_t  AliTPCcalibCalib::RefitTrack(AliESDtrack * track, AliTPCseed *seed, Float_t magesd){
 176   //
 177   // Refit track
 178   // if magesd==0 forget the curvature
 179
 180   //
 181   // 0 - Setup transform object
 182   //
 183   AliESDfriendTrack *friendTrack = (AliESDfriendTrack *)track->GetFriendTrack();
 184
 185   AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
 186   AliTPCParam     *param     = AliTPCcalibDB::Instance()->GetParameters();
 187   transform->SetCurrentRun(fRun);
 188   transform->SetCurrentTimeStamp((UInt_t)fTime);
 189   if(!fApplyExBCorrection) { // disable ExB correction in transform
 190     if(transform->GetCurrentRecoParam())
 191       transform->GetCurrentRecoParamNonConst()->SetUseExBCorrection(0);
 192   }
 193   if(!fApplyTOFCorrection) { // disable TOF correction in transform
 194     if(transform->GetCurrentRecoParam())
 195       transform->GetCurrentRecoParamNonConst()->SetUseTOFCorrection(kFALSE);
 196   }
 197
 198   //
 199   // First apply calibration
 200   //
 201   //  AliTPCPointCorrection * corr =  AliTPCPointCorrection::Instance();
 202   for (Int_t irow=0;irow<159;irow++) {
 203     AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
 204     if (!cluster) continue;
 205     AliTPCclusterMI cl0(*cluster);
 206     Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
 207     Int_t i[1]={cluster->GetDetector()};
 208
 209     transform->Transform(x,i,0,1);
 210     //
 211     // get position correction
 212     //
 213     Int_t ipad=0;
 214     if (cluster->GetDetector()>35) ipad=1;
 215     Float_t dy =0;//AliTPCClusterParam::SPosCorrection(0,ipad,cluster->GetPad(),cluster->GetTimeBin(),cluster->GetZ(),cluster->GetSigmaY2(),cluster->GetSigmaZ2(),cluster->GetMax());
 216     Float_t dz =0;//AliTPCClusterParam::SPosCorrection(1,ipad,cluster->GetPad(),cluster->GetTimeBin(),cluster->GetZ(),cluster->GetSigmaY2(),cluster->GetSigmaZ2(),cluster->GetMax());
 217     // if(fApplyPositionCorrection) {
 218     //x[1]-=dy;
 219     //x[2]-=dz;
 220     // }
 221     //
 222     // Apply sector alignment
 223     //
 224     Double_t dxq = AliTPCPointCorrection::SGetCorrectionSector(0,cluster->GetDetector()%36,cluster->GetX(),
 225                                                                cluster->GetY(),cluster->GetZ());
 226     Double_t dyq = AliTPCPointCorrection::SGetCorrectionSector(1,cluster->GetDetector()%36,cluster->GetX(),
 227                                                                cluster->GetY(),cluster->GetZ());
 228     Double_t dzq = AliTPCPointCorrection::SGetCorrectionSector(2,cluster->GetDetector()%36,cluster->GetX(),
 229                                                                cluster->GetY(),cluster->GetZ());
 230     if (0&fApplySectorAlignment){
 231       x[0]-=dxq;
 232       x[1]-=dyq;
 233       x[2]-=dzq;
 234     }
 235 //     //
 236 //     // Apply r-phi correction  - To be done on track level- knowing the track angle !!!
 237 //     //
 238 //     Double_t corrclY =
 239 //       corr->RPhiCOGCorrection(cluster->GetDetector(),cluster->GetRow(), cluster->GetPad(),
 240 //                                cluster->GetY(),cluster->GetY(), cluster->GetZ(), 0., cluster->GetMax(),2.5);
 241 //     // R correction
 242 //     Double_t corrR   = corr->CorrectionOutR0(kFALSE,kFALSE,cluster->GetX(),cluster->GetY(),cluster->GetZ(),cluster->GetDetector());
 243
 244 //     if (0&fApplyRPhiCorrection){
 245 //       if (cluster->GetY()>0) x[1]+=corrclY;  // rphi correction
 246 //       if (cluster->GetY()<0) x[1]-=corrclY;  // rphi correction
 247 //     }
 248 //     if (0&fApplyRCorrection){
 249 //       x[0]+=corrR;                           // radial correction
 250 //     }
 251
 252     //
 253     //
 254     //
 255     cluster->SetX(x[0]);
 256     cluster->SetY(x[1]);
 257     cluster->SetZ(x[2]);
 258
 259     //
 260     // Apply alignemnt
 261     //
 262     if (1){
 263       if (!param->IsGeoRead()) param->ReadGeoMatrices();
 264       TGeoHMatrix  *mat = param->GetClusterMatrix(cluster->GetDetector());
 265       //TGeoHMatrix  mat;
 266       Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
 267       Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
 268       if (mat) mat->LocalToMaster(pos,posC);
 269       else{
 270         // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
 271       }
 272       cluster->SetX(posC[0]);
 273       cluster->SetY(posC[1]);
 274       cluster->SetZ(posC[2]);
 275     }
 276
 277
 278
 279     if (fStreamLevel>2){
 280       TTreeSRedirector *cstream = GetDebugStreamer();
 281       if (cstream){
 282         (*cstream)<<"Clusters"<<
 283           "run="<<fRun<<              //  run number
 284           "event="<<fEvent<<          //  event number
 285           "time="<<fTime<<            //  time stamp of event
 286           "trigger="<<fTrigger<<      //  trigger
 287           "triggerClass="<<&fTriggerClass<<      //  trigger
 288           "mag="<<fMagF<<             //  magnetic field
 289           "cl0.="<<&cl0<<
 290           "cl.="<<cluster<<
 291           "cy="<<dy<<
 292           "cz="<<dz<<
 293           //      "cY="<<corrclY<<
 294           //      "cR="<<corrR<<
 295           "dxq="<<dxq<<
 296           "dyq="<<dyq<<
 297           "dzq="<<dzq<<
 298           "\n";
 299       }
 300     }
 301   }
 302   //
 303   //
 304   //
 305   Int_t ncl = seed->GetNumberOfClusters();
 306   const Double_t kResetCov=4.;
 307   const Double_t kSigma=5.;
 308   Double_t covar[15];
 309   for (Int_t i=0;i<15;i++) covar[i]=0;
 310   covar[0]=kSigma*kSigma;
 311   covar[2]=kSigma*kSigma;
 312   covar[5]=kSigma*kSigma/Float_t(ncl*ncl);
 313   covar[9]=kSigma*kSigma/Float_t(ncl*ncl);
 314   covar[14]=0.2*0.2;
 315   if (TMath::Abs(magesd)<0.05) {
 316      covar[14]=0.025*0.025;
 317   }
 318   //
 319   // And now do refit
 320   //
 321   AliExternalTrackParam * trackInOld  = (AliExternalTrackParam*)track->GetInnerParam();
 322   AliExternalTrackParam * trackOuter = (AliExternalTrackParam*)track->GetOuterParam();
 323   AliExternalTrackParam * trackOutOld   = (AliExternalTrackParam *)friendTrack->GetTPCOut();
 324   Double_t mass =    TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 325
 326
 327
 328   AliExternalTrackParam trackIn  = *trackOutOld;
 329   trackIn.ResetCovariance(kResetCov);
 330   trackIn.AddCovariance(covar);
 331   if (TMath::Abs(magesd)<0.05) {
 332     ((Double_t&)(trackIn.GetParameter()[4]))=0.000000001;
 333     ((Double_t&)(trackIn.GetCovariance()[14]))=covar[14];  // fix the line
 334   }
 335
 336   Double_t xyz[3];
 337   Int_t nclIn=0,nclOut=0;
 338   //
 339   // Refit in
 340   //
 341   for (Int_t irow=159; irow>0; irow--){
 342     AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
 343     if (!cl) continue;
 344     if (cl->GetX()<80) continue;
 345     Int_t sector = cl->GetDetector();
 346     Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackIn.GetAlpha();
 347     if (TMath::Abs(dalpha)>0.01){
 348       if (!trackIn.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
 349     }
 350     Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
 351     Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
 352     AliTPCseed::GetError(cl, &trackIn,cov[0],cov[2]);
 353     cov[0]*=cov[0];
 354     cov[2]*=cov[2];
 355     trackIn.GetXYZ(xyz);
 356     //    Double_t bz = AliTracker::GetBz(xyz);
 357
 358     //    if (!trackIn.PropagateTo(r[0],bz)) continue;
 359     if (!AliTracker::PropagateTrackToBxByBz(&trackIn, r[0],mass,0.1,kFALSE)) continue;
 360
 361     if (RejectCluster(cl,&trackIn)) continue;
 362     nclIn++;
 363     trackIn.Update(&r[1],cov);
 364   }
 365   //
 366   AliExternalTrackParam trackOut = trackIn;
 367   trackOut.ResetCovariance(kResetCov);
 368   trackOut.AddCovariance(covar);
 369   if (TMath::Abs(magesd)<0.05) {
 370     ((Double_t&)(trackOut.GetParameter()[4]))=0.000000001;
 371     ((Double_t&)(trackOut.GetCovariance()[14]))=covar[14];  // fix the line
 372   }
 373
 374   //
 375   // Refit out
 376   //
 377   //Bool_t lastEdge=kFALSE;
 378   for (Int_t irow=0; irow<160; irow++){
 379     AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
 380     if (!cl) continue;
 381     if (cl->GetX()<80) continue;
 382     Int_t sector = cl->GetDetector();
 383     Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackOut.GetAlpha();
 384
 385     if (TMath::Abs(dalpha)>0.01){
 386       if (!trackOut.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
 387     }
 388     Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
 389
 390     Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
 391     AliTPCseed::GetError(cl, &trackOut,cov[0],cov[2]);
 392     cov[0]*=cov[0];
 393     cov[2]*=cov[2];
 394     trackOut.GetXYZ(xyz);
 395     //Double_t bz = AliTracker::GetBz(xyz);
 396     //    if (!trackOut.PropagateTo(r[0],bz)) continue;
 397     if (!AliTracker::PropagateTrackToBxByBz(&trackOut, r[0],mass,0.1,kFALSE)) continue;
 398
 399     if (RejectCluster(cl,&trackOut)) continue;
 400     nclOut++;
 401     trackOut.Update(&r[1],cov);
 402     //if (cl->GetType()<0) lastEdge=kTRUE;
 403     //if (cl->GetType()>=0) lastEdge=kFALSE;
 404   }
 405   //
 406   //
 407   //
 408   nclIn=0;
 409   trackIn  = trackOut;
 410   trackIn.ResetCovariance(kResetCov);
 411   if (TMath::Abs(magesd)<0.05) {
 412     ((Double_t&)(trackIn.GetParameter()[4]))=0.000000001;
 413     ((Double_t&)(trackIn.GetCovariance()[14]))=covar[14];  // fix the line
 414   }
 415   //
 416   // Refit in one more time
 417   //
 418   for (Int_t irow=159; irow>0; irow--){
 419     AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
 420     if (!cl) continue;
 421     if (cl->GetX()<80) continue;
 422     Int_t sector = cl->GetDetector();
 423     Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackIn.GetAlpha();
 424     if (TMath::Abs(dalpha)>0.01){
 425       if (!trackIn.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
 426     }
 427     Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
 428     Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
 429     AliTPCseed::GetError(cl, &trackIn,cov[0],cov[2]);
 430     cov[0]*=cov[0];
 431     cov[2]*=cov[2];
 432     trackIn.GetXYZ(xyz);
 433     //Double_t bz = AliTracker::GetBz(xyz);
 434
 435     //    if (!trackIn.PropagateTo(r[0],bz)) continue;
 436     if (!AliTracker::PropagateTrackToBxByBz(&trackIn, r[0],mass,0.1,kFALSE)) continue;
 437
 438     if (RejectCluster(cl,&trackIn)) continue;
 439     nclIn++;
 440     trackIn.Update(&r[1],cov);
 441   }
 442
 443
 444   trackIn.Rotate(trackInOld->GetAlpha());
 445   trackOut.Rotate(trackOutOld->GetAlpha());
 446   //
 447   trackInOld->GetXYZ(xyz);
 448   Double_t bz = AliTracker::GetBz(xyz);
 449   trackIn.PropagateTo(trackInOld->GetX(),bz);
 450   //
 451   trackOutOld->GetXYZ(xyz);
 452   bz = AliTracker::GetBz(xyz);
 453   trackOut.PropagateTo(trackOutOld->GetX(),bz);
 454
 455
 456   if (fStreamLevel>0){
 457     TTreeSRedirector *cstream = GetDebugStreamer();
 458     if (cstream){
 459       (*cstream)<<"Tracks"<<
 460         "run="<<fRun<<              //  run number
 461         "event="<<fEvent<<          //  event number
 462         "time="<<fTime<<            //  time stamp of event
 463         "trigger="<<fTrigger<<      //  trigger
 464         "triggerClass="<<&fTriggerClass<<      //  trigger
 465         "mag="<<fMagF<<             //  magnetic field
 466         "nclIn="<<nclIn<<
 467         "nclOut="<<nclOut<<
 468         "ncl="<<ncl<<
 469         "TrIn0.="<<trackInOld<<
 470         "TrOut0.="<<trackOutOld<<
 471         "TrIn1.="<<&trackIn<<
 472         "TrOut1.="<<&trackOut<<
 473         "\n";
 474     }
 475   }
 476   //
 477   // And now rewrite ESDtrack and TPC seed
 478   //
 479
 480   (*trackInOld)  = trackIn;
 481   (*trackOutOld) = trackOut;
 482   (*trackOuter) = trackOut;
 483   AliExternalTrackParam *t = &trackIn;
 484   track->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
 485   seed->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
 486   seed->SetNumberOfClusters((nclIn+nclOut)/2);
 487   return kTRUE;
 488 }
 489
 490
 491
 492 Bool_t AliTPCcalibCalib::RejectCluster(AliTPCclusterMI* cl, AliExternalTrackParam * param){
 493   //
 494   // check the acceptance of cluster
 495   // Cut on edge effects
 496   //
 497   Float_t kEdgeCut=2.5;
 498   Float_t kSigmaCut=6;
 499
 500   Bool_t isReject = kFALSE;
 501   Float_t edgeY = cl->GetX()*TMath::Tan(TMath::Pi()/18);
 502   Float_t dist  = edgeY - TMath::Abs(cl->GetY());
 503   if (param)  dist  = TMath::Abs(edgeY - TMath::Abs(param->GetY()));
 504   if (dist<kEdgeCut) isReject=kTRUE;
 505
 506   Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
 507   AliTPCseed::GetError(cl, param,cov[0],cov[2]);
 508   if (param->GetSigmaY2()<0 || param->GetSigmaZ2()<0){
 509     AliError("Wrong parameters");
 510     return kFALSE;
 511   }
 512   Double_t py = (cl->GetY()-param->GetY())/TMath::Sqrt(cov[0]*cov[0]+param->GetSigmaY2());
 513   Double_t pz = (cl->GetZ()-param->GetZ())/TMath::Sqrt(cov[2]*cov[2]+param->GetSigmaZ2());
 514   //
 515   if ((py*py+pz*pz)>kSigmaCut*kSigmaCut) isReject=kTRUE;
 516
 517   return isReject;
 518 }
 519
 520
 521