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