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 "AliTPCclusterMI.h"
  62 #include "AliTPCseed.h"
  63 #include "AliTPCPointCorrection.h"
  64
  65 ClassImp(AliTPCcalibCalib)
  66
  67 AliTPCcalibCalib::AliTPCcalibCalib():
  68     AliTPCcalibBase()
  69 {
  70   //
  71   // Constructor
  72   //
  73 }
  74
  75
  76 AliTPCcalibCalib::AliTPCcalibCalib(const Text_t *name, const Text_t *title)
  77   :AliTPCcalibBase()
  78 {
  79   SetName(name);
  80   SetTitle(title);
  81 }
  82
  83
  84 AliTPCcalibCalib::AliTPCcalibCalib(const AliTPCcalibCalib&calib):
  85   AliTPCcalibBase(calib)
  86 {
  87   //
  88   // copy constructor
  89   //
  90 }
  91
  92 AliTPCcalibCalib &AliTPCcalibCalib::operator=(const AliTPCcalibCalib&calib){
  93   //
  94   //
  95   //
  96   ((AliTPCcalibBase *)this)->operator=(calib);
  97   return *this;
  98 }
  99
 100
 101 AliTPCcalibCalib::~AliTPCcalibCalib() {
 102   //
 103   // destructor
 104   //
 105 }
 106
 107
 108 void     AliTPCcalibCalib::Process(AliESDEvent *event){
 109   //
 110   //
 111   //
 112   if (!event) {
 113     return;
 114   }
 115   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
 116   if (!ESDfriend) {
 117    return;
 118   }
 119
 120   if (GetDebugLevel()>20) printf("Hallo world: Im here\n");
 121   Int_t ntracks=event->GetNumberOfTracks();
 122   //AliTPCcalibDB::Instance()->SetExBField(fMagF);
 123
 124   //
 125   //
 126   //
 127
 128   for (Int_t i=0;i<ntracks;++i) {
 129     AliESDtrack *track = event->GetTrack(i);
 130     const AliExternalTrackParam * trackIn = track->GetInnerParam();
 131     const AliExternalTrackParam * trackOut = track->GetOuterParam();
 132     if (!trackIn) continue;
 133     if (!trackOut) continue;
 134
 135     AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
 136     TObject *calibObject;
 137     AliTPCseed *seed = 0;
 138     for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
 139       if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
 140     }
 141     if (!seed) continue;
 142     RefitTrack(track, seed);
 143   }
 144   return;
 145 }
 146
 147 Bool_t  AliTPCcalibCalib::RefitTrack(AliESDtrack * track, AliTPCseed *seed){
 148   //
 149   // Refit track
 150   //
 151
 152   //
 153   // 0 - Setup transform object
 154   //
 155   AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
 156   transform->SetCurrentRun(fRun);
 157   transform->SetCurrentTimeStamp((UInt_t)fTime);
 158   //
 159   // First apply calibration
 160   //
 161   AliTPCPointCorrection * corr =  AliTPCPointCorrection::Instance();
 162   for (Int_t irow=0;irow<159;irow++) {
 163     AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
 164     if (!cluster) continue;
 165     AliTPCclusterMI cl0(*cluster);
 166     Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
 167     Int_t i[1]={cluster->GetDetector()};
 168     transform->Transform(x,i,0,1);
 169     //
 170     // get position correction
 171     //
 172     Int_t ipad=0;
 173     if (cluster->GetDetector()>35) ipad=1;
 174     Float_t dy =AliTPCClusterParam::SPosCorrection(0,ipad,cluster->GetPad(),cluster->GetTimeBin(),cluster->GetZ(),cluster->GetSigmaY2(),cluster->GetSigmaZ2(),cluster->GetMax());
 175     Float_t dz =AliTPCClusterParam::SPosCorrection(1,ipad,cluster->GetPad(),cluster->GetTimeBin(),cluster->GetZ(),cluster->GetSigmaY2(),cluster->GetSigmaZ2(),cluster->GetMax());
 176     //x[1]-=dy;
 177     //x[2]-=dz;
 178     //
 179     // Apply sector alignment
 180     //
 181     Double_t dxq = AliTPCPointCorrection::SGetCorrectionSector(0,cluster->GetDetector()%36,cluster->GetX(),
 182                                                                cluster->GetY(),cluster->GetZ());
 183     Double_t dyq = AliTPCPointCorrection::SGetCorrectionSector(1,cluster->GetDetector()%36,cluster->GetX(),
 184                                                                cluster->GetY(),cluster->GetZ());
 185     Double_t dzq = AliTPCPointCorrection::SGetCorrectionSector(2,cluster->GetDetector()%36,cluster->GetX(),
 186                                                                cluster->GetY(),cluster->GetZ());
 187     if (kTRUE){
 188       x[0]-=dxq;
 189       x[1]-=dyq;
 190       x[2]-=dzq;
 191     }
 192     //
 193     // Apply r-phi correction  - To be done on track level- knowing the track angle !!!
 194     //
 195     Double_t corrclY =
 196       corr->RPhiCOGCorrection(cluster->GetDetector(),cluster->GetRow(), cluster->GetPad(),
 197                                   cluster->GetY(),cluster->GetY(), cluster->GetZ(), 0., cluster->GetMax(),2.5);
 198     // R correction
 199     Double_t corrR   = corr->CorrectionOutR0(kFALSE,kFALSE,cluster->GetX(),cluster->GetY(),cluster->GetZ(),cluster->GetDetector());
 200
 201     if (kTRUE){
 202       if (cluster->GetY()>0) x[1]+=corrclY;  // rphi correction
 203       if (cluster->GetY()<0) x[1]-=corrclY;  // rphi correction
 204       x[0]+=corrR;                           // radial correction
 205     }
 206
 207     //
 208     //
 209     //
 210     cluster->SetX(x[0]);
 211     cluster->SetY(x[1]);
 212     cluster->SetZ(x[2]);
 213     if (fStreamLevel>2){
 214       TTreeSRedirector *cstream = GetDebugStreamer();
 215       if (cstream){
 216         (*cstream)<<"Clusters"<<
 217           "run="<<fRun<<              //  run number
 218           "event="<<fEvent<<          //  event number
 219           "time="<<fTime<<            //  time stamp of event
 220           "trigger="<<fTrigger<<      //  trigger
 221           "triggerClass="<<&fTriggerClass<<      //  trigger
 222           "mag="<<fMagF<<             //  magnetic field
 223           "cl0.="<<&cl0<<
 224           "cl.="<<cluster<<
 225           "cy="<<dy<<
 226           "cz="<<dz<<
 227           "cY="<<corrclY<<
 228           "cR="<<corrR<<
 229           "dxq="<<dxq<<
 230           "dyq="<<dyq<<
 231           "dzq="<<dzq<<
 232           "\n";
 233       }
 234     }
 235   }
 236   //
 237   //
 238   //
 239   Int_t ncl = seed->GetNumberOfClusters();
 240   Double_t covar[15];
 241   for (Int_t i=0;i<15;i++) covar[i]=0;
 242   covar[0]=10.*10.;
 243   covar[2]=10.*10.;
 244   covar[5]=10.*10./(64.*64.);
 245   covar[9]=10.*10./(64.*64.);
 246   covar[14]=1*1;
 247
 248   //
 249   // And now do refit
 250   //
 251   AliExternalTrackParam * trackInOld  = (AliExternalTrackParam*)track->GetInnerParam();
 252   AliExternalTrackParam * trackOutOld = (AliExternalTrackParam*)track->GetOuterParam();
 253
 254   AliExternalTrackParam trackIn  = *trackOutOld;
 255   trackIn.ResetCovariance(200.);
 256   trackIn.AddCovariance(covar);
 257   Double_t xyz[3];
 258   Int_t nclIn=0,nclOut=0;
 259   //
 260   // Refit in
 261   //
 262
 263   for (Int_t irow=159; irow>0; irow--){
 264     AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
 265     if (!cl) continue;
 266     if (cl->GetX()<80) continue;
 267     Int_t sector = cl->GetDetector();
 268     Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackIn.GetAlpha();
 269     if (TMath::Abs(dalpha)>0.01)
 270       trackIn.Rotate(TMath::DegToRad()*(sector%18*20.+10.));
 271     Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
 272     Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
 273     AliTPCseed::GetError(cl, &trackIn,cov[0],cov[2]);
 274     cov[0]*=cov[0];
 275     cov[2]*=cov[2];
 276     trackIn.GetXYZ(xyz);
 277     Double_t bz = AliTracker::GetBz(xyz);
 278
 279     if (!trackIn.PropagateTo(r[0],bz)) continue;
 280     if (RejectCluster(cl,&trackIn)) continue;
 281     nclIn++;
 282     trackIn.Update(&r[1],cov);
 283   }
 284   //
 285   AliExternalTrackParam trackOut = trackIn;
 286   trackOut.ResetCovariance(200.);
 287   trackOut.AddCovariance(covar);
 288   //
 289   // Refit out
 290   //
 291   //Bool_t lastEdge=kFALSE;
 292   for (Int_t irow=0; irow<160; irow++){
 293     AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
 294     if (!cl) continue;
 295     if (cl->GetX()<80) continue;
 296     Int_t sector = cl->GetDetector();
 297     Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackOut.GetAlpha();
 298
 299     if (TMath::Abs(dalpha)>0.01)
 300       trackOut.Rotate(TMath::DegToRad()*(sector%18*20.+10.));
 301     Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
 302
 303     Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
 304     AliTPCseed::GetError(cl, &trackOut,cov[0],cov[2]);
 305     cov[0]*=cov[0];
 306     cov[2]*=cov[2];
 307     trackOut.GetXYZ(xyz);
 308     Double_t bz = AliTracker::GetBz(xyz);
 309     if (!trackOut.PropagateTo(r[0],bz)) continue;
 310     if (RejectCluster(cl,&trackOut)) continue;
 311     nclOut++;
 312     trackOut.Update(&r[1],cov);
 313     //if (cl->GetType()<0) lastEdge=kTRUE;
 314     //if (cl->GetType()>=0) lastEdge=kFALSE;
 315   }
 316   //
 317   //
 318   //
 319   nclIn=0;
 320   trackIn  = trackOut;
 321   trackIn.ResetCovariance(10.);
 322   //
 323   // Refit in one more time
 324   //
 325   for (Int_t irow=159; irow>0; irow--){
 326     AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
 327     if (!cl) continue;
 328     if (cl->GetX()<80) continue;
 329     Int_t sector = cl->GetDetector();
 330     Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackIn.GetAlpha();
 331     if (TMath::Abs(dalpha)>0.01)
 332       trackIn.Rotate(TMath::DegToRad()*(sector%18*20.+10.));
 333     Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
 334     Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
 335     AliTPCseed::GetError(cl, &trackIn,cov[0],cov[2]);
 336     cov[0]*=cov[0];
 337     cov[2]*=cov[2];
 338     trackIn.GetXYZ(xyz);
 339     Double_t bz = AliTracker::GetBz(xyz);
 340
 341     if (!trackIn.PropagateTo(r[0],bz)) continue;
 342     if (RejectCluster(cl,&trackIn)) continue;
 343     nclIn++;
 344     trackIn.Update(&r[1],cov);
 345   }
 346
 347
 348   trackIn.Rotate(trackInOld->GetAlpha());
 349   trackOut.Rotate(trackOutOld->GetAlpha());
 350   //
 351   trackInOld->GetXYZ(xyz);
 352   Double_t bz = AliTracker::GetBz(xyz);
 353   trackIn.PropagateTo(trackInOld->GetX(),bz);
 354   //
 355   trackOutOld->GetXYZ(xyz);
 356   bz = AliTracker::GetBz(xyz);
 357   trackOut.PropagateTo(trackOutOld->GetX(),bz);
 358
 359
 360   if (fStreamLevel>0){
 361     TTreeSRedirector *cstream = GetDebugStreamer();
 362     if (cstream){
 363       (*cstream)<<"Tracks"<<
 364         "run="<<fRun<<              //  run number
 365         "event="<<fEvent<<          //  event number
 366         "time="<<fTime<<            //  time stamp of event
 367         "trigger="<<fTrigger<<      //  trigger
 368         "triggerClass="<<&fTriggerClass<<      //  trigger
 369         "mag="<<fMagF<<             //  magnetic field
 370         "nclIn="<<nclIn<<
 371         "nclOut="<<nclOut<<
 372         "ncl="<<ncl<<
 373         "TrIn0.="<<trackInOld<<
 374         "TrOut0.="<<trackOutOld<<
 375         "TrIn1.="<<&trackIn<<
 376         "TrOut1.="<<&trackOut<<
 377         "\n";
 378     }
 379   }
 380   //
 381   // And now rewrite ESDtrack and TPC seed
 382   //
 383
 384   (*trackInOld)  = trackIn;
 385   (*trackOutOld) = trackOut;
 386   AliExternalTrackParam *t = &trackIn;
 387   track->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
 388   seed->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
 389   seed->SetNumberOfClusters((nclIn+nclOut)/2);
 390   return kTRUE;
 391 }
 392
 393
 394
 395 Bool_t AliTPCcalibCalib::RejectCluster(AliTPCclusterMI* cl, AliExternalTrackParam * param){
 396   //
 397   // check the acceptance of cluster
 398   // Cut on edge effects
 399   //
 400   Float_t kEdgeCut=2.5;
 401   Float_t kSigmaCut=6;
 402
 403   Bool_t isReject = kFALSE;
 404   Float_t edgeY = cl->GetX()*TMath::Tan(TMath::Pi()/18);
 405   Float_t dist  = edgeY - TMath::Abs(cl->GetY());
 406   if (param)  dist  = TMath::Abs(edgeY - TMath::Abs(param->GetY()));
 407   if (dist<kEdgeCut) isReject=kTRUE;
 408
 409   Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
 410   AliTPCseed::GetError(cl, param,cov[0],cov[2]);
 411   Double_t py = (cl->GetY()-param->GetY())/TMath::Sqrt(cov[0]*cov[0]+param->GetSigmaY2());
 412   Double_t pz = (cl->GetZ()-param->GetZ())/TMath::Sqrt(cov[2]*cov[2]+param->GetSigmaZ2());
 413   //
 414   if ((py*py+pz*pz)>kSigmaCut*kSigmaCut) isReject=kTRUE;
 415
 416   return isReject;
 417 }
 418
 419
 420