TPC/AliTPCCosmicUtils.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 // Static function member
  17 // as utils for
  18 // AliTPCCombinedTrackfit
  19 //
  20 // grep "exitreport" in output log to check abnormal termination
  21 // comment tbc
  22 //
  23 //  Xianguo Lu <lu@physi.uni-heidelberg.de>
  24 //
  25
  26 #include <TAxis.h>
  27 #include <TCanvas.h>
  28 #include <TGraph.h>
  29 #include <TTreeStream.h>
  30 #include <TVector3.h>
  31
  32 #include "AliESDtrack.h"
  33 #include "AliESDfriendTrack.h"
  34 #include "AliTPCseed.h"
  35 #include "AliTrackerBase.h"
  36 #include "AliTrackPointArray.h"
  37
  38 #include "AliTPCCosmicUtils.h"
  39
  40 Int_t AliTPCCosmicUtils::RotateSafe(AliExternalTrackParam *trackPar, const Double_t aa)
  41 {
  42   //1. in AliExternalTrackParam::GetXYZ
  43   //r[0]=fX; r[1]=fP[0]; r[2]=fP[1];
  44   //return Local2GlobalPosition(r,fAlpha);
  45   //in AliVParticle::Local2GlobalMomentum
  46   //Double_t cs=TMath::Cos(alpha), sn=TMath::Sin(alpha);
  47   //Double_t r=TMath::Sqrt((1. - p[1])*(1. + p[1]));
  48   //p[0]=pt*(r*cs - p[1]*sn); p[1]=pt*(p[1]*cs + r*sn); p[2]=pt*p[2];
  49   //which is cos(phi_local=Snp+fAlpha) > 0 always.
  50   //2. also in Bool_t AliExternalTrackParam::Rotate(Double_t alpha)
  51   // Double_t sf=fP2, cf=TMath::Sqrt((1.- fP2)*(1.+fP2));
  52   //in AliExternalTrackParam::Set
  53   //mom.RotateZ(-fAlpha);
  54   //fP[2] = TMath::Sin(mom.Phi());
  55   //since only sin is used for mom.Phi(), that is assuming cos(mom.Phi())>0
  56
  57   Double_t xyz[3], pxpypz[3];
  58   trackPar->GetXYZ(xyz);
  59   trackPar->GetPxPyPz(pxpypz);
  60   const TVector3 pos0(xyz);
  61   const TVector3 mom0(pxpypz);
  62   TVector3 pos1, mom1, dpos, dmom;
  63   const Double_t eps = 1e-6;
  64   AliExternalTrackParam tmppar = (*trackPar);
  65
  66   if(!tmppar.Rotate(aa)){
  67     return 0;
  68   }
  69
  70   tmppar.GetXYZ(xyz);
  71   tmppar.GetPxPyPz(pxpypz);
  72   pos1.SetXYZ(xyz[0], xyz[1], xyz[2]);
  73   mom1.SetXYZ(pxpypz[0], pxpypz[1], pxpypz[2]);
  74   dpos = pos1-pos0;
  75   dmom = mom1-mom0;
  76   if(dpos.Mag()<eps && dmom.Mag()<eps){
  77     (*trackPar)=tmppar;
  78     return 1;
  79   }
  80   else
  81     return 0;
  82
  83   /*
  84   //flip
  85   tmppar = (*trackPar);
  86   if(!tmppar.Rotate(aa+TMath::Pi())){
  87     return 0;
  88   }
  89
  90   tmppar.GetXYZ(xyz);
  91   tmppar.GetPxPyPz(pxpypz);
  92   pos1.SetXYZ(xyz[0], xyz[1], xyz[2]);
  93   mom1.SetXYZ(pxpypz[0], pxpypz[1], pxpypz[2]);
  94   dpos = pos1-pos0;
  95   dmom = mom1-mom0;
  96   if(dpos.Mag()<eps && dmom.Mag()<eps){
  97     (*trackPar)=tmppar;
  98     return -1;
  99   }
 100
 101   return 0;
 102   */
 103 }
 104
 105 void AliTPCCosmicUtils::PrintTrackParam(const Int_t id, const AliExternalTrackParam * trackpar, const char *tag)
 106 {
 107   //
 108   //print out TrackParam
 109   //
 110   Double_t xyz[3]={-999,-999,-999};
 111   trackpar->GetXYZ(xyz);
 112   Double_t pxpypz[3]={-999,-999,-999};
 113   trackpar->GetPxPyPz(pxpypz);
 114
 115   printf("PrintTrackPar %2s %3d : %6.1f %5.1f : %6.6f %6.6f %6.6f : %8.3f %8.3f %8.3f : (%2d) %8.3f %11.6f %.2f\n", tag, id, trackpar->GetX(), sqrt(pow(xyz[0],2)+pow(xyz[1],2)), xyz[0], xyz[1], xyz[2], pxpypz[0], pxpypz[1], pxpypz[2], trackpar->Charge(), trackpar->Pt(), trackpar->P(), trackpar->GetAlpha()*TMath::RadToDeg());
 116 }
 117
 118 void AliTPCCosmicUtils::DrawSeeds(const AliTPCseed * seeds[], const TString tag, const TString outputformat)
 119 {
 120   //
 121   //draw seed and output to file
 122   //
 123   TGraph *grsyx[]={new TGraph, new TGraph};
 124   TGraph *grsyz[]={new TGraph, new TGraph};
 125
 126   for(Int_t itrk=0; itrk<2; itrk++){
 127     if(!seeds[itrk])
 128       continue;
 129
 130     Int_t ncl = 0;
 131     for(Int_t irow=0; irow<fgkNRow; irow++){
 132       const AliTPCclusterMI * cls = seeds[itrk]->GetClusterPointer(irow);
 133       if(!cls)
 134         continue;
 135
 136       Float_t xyz[3]={-999,-999,-999};
 137       cls->GetGlobalXYZ(xyz);
 138       //printf("Test %f %f %f \n", xyz[0], xyz[1], xyz[2]);
 139       grsyx[itrk]->SetPoint(ncl, xyz[0], xyz[1]);
 140       grsyz[itrk]->SetPoint(ncl, xyz[2], xyz[1]);
 141       ncl++;
 142     }
 143   }
 144
 145   grsyx[0]->SetTitle(tag);
 146   grsyx[0]->SetMaximum(250);
 147   grsyx[0]->SetMinimum(-250);
 148   grsyx[0]->GetXaxis()->SetLimits(-250,250);
 149   grsyx[0]->GetXaxis()->SetTitle("X (cm)");
 150   grsyx[0]->GetYaxis()->SetTitle("Y (cm)");
 151
 152   grsyx[0]->SetMarkerStyle(20);
 153   grsyx[0]->SetMarkerColor(kRed);
 154   grsyx[1]->SetMarkerStyle(22);
 155   grsyx[1]->SetMarkerColor(kBlue);
 156
 157   grsyz[0]->SetTitle(tag);
 158   grsyz[0]->SetMaximum(250);
 159   grsyz[0]->SetMinimum(-250);
 160   grsyz[0]->GetXaxis()->SetLimits(-250,250);
 161   grsyz[0]->GetXaxis()->SetTitle("Z (cm)");
 162   grsyz[0]->GetYaxis()->SetTitle("Y (cm)");
 163
 164   grsyz[0]->SetMarkerStyle(20);
 165   grsyz[0]->SetMarkerColor(kRed);
 166   grsyz[1]->SetMarkerStyle(22);
 167   grsyz[1]->SetMarkerColor(kBlue);
 168
 169   TCanvas *cc=new TCanvas("cc","",1200,600);
 170   cc->Divide(2);
 171   cc->cd(1);
 172   for(int itrk=0; itrk<2; itrk++){
 173     grsyx[itrk]->SetMarkerSize(1);
 174     grsyx[itrk]->Draw(itrk?"lp same":"alp");
 175   }
 176   cc->cd(2);
 177   for(int itrk=0; itrk<2; itrk++){
 178     grsyz[itrk]->SetMarkerSize(1);
 179     grsyz[itrk]->Draw(itrk?"lp same":"alp");
 180   }
 181   cc->Print(Form("drawTrack%s.%s", tag.Data(), outputformat.Data()));
 182
 183   for(Int_t ii=0; ii<2; ii++){
 184     delete grsyx[ii];
 185     delete grsyz[ii];
 186   }
 187   delete cc;
 188 }
 189
 190 //================================================================================================
 191 AliTPCseed * AliTPCCosmicUtils::GetTPCseed(const AliESDtrack *esdtrack)
 192 {
 193   //
 194   //Get TPC seeds from ESDfriendTrack
 195   //
 196
 197   AliESDfriendTrack * friendtrk = (AliESDfriendTrack *) esdtrack->GetFriendTrack();
 198   if(!friendtrk){
 199     printf("exitreport AliTPCCosmicUtils::GetTPCseed no friend track!\n");
 200     exit(1);
 201   }
 202
 203   TObject *calibObject=0x0;
 204   AliTPCseed *tseed = 0x0;
 205   for(Int_t l=0; (calibObject=friendtrk->GetCalibObject(l)); l++) {
 206     if( (tseed=dynamic_cast<AliTPCseed*>(calibObject)) )
 207       break;
 208   }
 209   return tseed;
 210 }
 211
 212 //================================================================================================
 213 AliExternalTrackParam *AliTPCCosmicUtils::MakeSeed(const AliTPCseed *tseed)
 214 {
 215   //
 216   //make seed for propagation of TrackParam, using np = 3 outer clusters (separated by deltancls clusters) in TPCseed
 217   //
 218
 219   const Int_t rowstart = fgkNRow-1;
 220   const Int_t rowstop = 0;
 221   const Int_t drow = -1;
 222
 223   //---
 224   const Int_t np = 3;
 225   AliTrackPoint *tpos[np];
 226   for(Int_t ii=0; ii<np; ii++)
 227     tpos[ii] = 0x0;
 228
 229   const Float_t cov[6]={0,0,0, 0.01*0.01 ,0, 0.01*0.01};
 230
 231   //---
 232   Int_t npos = 0;
 233   Int_t icl = 0;
 234   const Int_t deltancls = fgkNclsMin/2-1;
 235   Int_t oldcl = -deltancls;
 236
 237   for(Int_t irow=rowstart; drow*irow<=drow*rowstop; irow+=drow){
 238     AliTPCclusterMI *cls = tseed->GetClusterPointer(irow);
 239     if(!cls) {
 240       continue;
 241     }
 242
 243     if( icl == (oldcl+deltancls) ){
 244       Float_t txyz[3];
 245       cls->GetGlobalXYZ(txyz);
 246       tpos[npos++] = new AliTrackPoint(txyz, cov, 0);
 247       //printf("------ %d %f %f %f\n", npos, txyz[0], txyz[1], txyz[2]);
 248
 249       oldcl = icl;
 250       if(npos==np) break;
 251     }
 252     icl++;
 253   }
 254   if(npos!=np){
 255     for(Int_t ii=0; ii<npos; ii++)
 256       delete tpos[ii];
 257
 258     return 0x0;
 259   }
 260
 261   AliExternalTrackParam * trackparam = AliTrackerBase::MakeSeed(*(tpos[0]), *(tpos[1]), *(tpos[2]));
 262   if(!trackparam || trackparam->Pt()==0){
 263     for(Int_t ii=0; ii<npos; ii++)
 264       delete tpos[ii];
 265     delete trackparam;
 266
 267     return 0x0;
 268   }
 269
 270   //------
 271
 272   Double_t sxyz[3]={-999,-999,-999}, spxpypz[3]={-999,-999,-999};
 273   trackparam->GetXYZ(sxyz);
 274   trackparam->GetPxPyPz(spxpypz);
 275   Double_t scov[21];
 276   Int_t sign = trackparam->Charge();
 277
 278   //reset covariance matrix -- necessary, otherwise bad fitting result: trackparam->GetCovariance()[ii] has problem: some are strange values
 279   for(Int_t ii=0; ii<21; ii++) {
 280     scov[ii]=0;
 281   }
 282   trackparam->Set(sxyz, spxpypz, scov, sign);
 283
 284   for(Int_t ii=0; ii<np; ii++)
 285     delete tpos[ii];
 286
 287   return trackparam;
 288 }
 289
 290 //================================================================================================
 291 void AliTPCCosmicUtils::IniCov(AliExternalTrackParam *trackPar, const Double_t ncl)
 292 {
 293   //
 294   //initialize covariance matrix
 295   //
 296
 297   const Double_t ksigma=5.;
 298   Double_t acov[16];
 299   for (Int_t i=0;i<15;i++)
 300     acov[i]=0;
 301
 302   acov[0]=ksigma*ksigma;
 303   acov[2]=ksigma*ksigma;
 304   acov[5]=ksigma*ksigma;
 305   acov[9]=ksigma*ksigma;
 306   acov[14]=0.2*0.2;
 307
 308   acov[5] = ksigma*ksigma/(ncl*ncl);
 309   acov[9] = ksigma*ksigma/(ncl*ncl);
 310
 311   const Double_t resetcov = 4;
 312
 313   trackPar->ResetCovariance(resetcov);
 314   //the following helps a lot!!
 315   trackPar->AddCovariance(acov);
 316 }
 317
 318 void AliTPCCosmicUtils::SingleFit(AliExternalTrackParam * trackInOld, AliExternalTrackParam * trackOutOld, const AliTPCseed *tseed, const Bool_t kinward,  Int_t &nfit, Int_t &nmiss, Double_t &pchi2, TTreeSRedirector *debugstreamer)
 319 {
 320   //
 321   //fit single track
 322   //
 323   //kinward is the true geometry of the track. Incomming track: 1; outgoing track: 0
 324   const Double_t inde  = kinward? -1 :  1;
 325   const Double_t outde = kinward?  1 : -1;
 326
 327   //PrintTrackParam(9000, trackOutOld);
 328
 329   AliExternalTrackParam trackOut = *trackOutOld;
 330   AliExternalTrackParam trackIn;
 331   Int_t ksite = -999;
 332
 333   //nmiss is from the 2 FitKernel of the last iteration
 334   //nfit, pchi2 is from the last FitKernel of the last iteration
 335   for(Int_t ii=0; ii<fgkNiter; ii++){
 336     nmiss = 0;
 337
 338     ksite = -999;
 339     trackIn  = trackOut;
 340     FitKernel(&trackIn,  tseed, fgkNRow-1, 0, inde , ksite, nfit, nmiss, pchi2, 0x0, kTRUE);
 341     //PrintTrackParam(9010+ii, &trackIn);
 342
 343     nfit = 0;
 344     pchi2 = 0;
 345
 346     ksite = -999;
 347     trackOut = trackIn;
 348     FitKernel(&trackOut, tseed, 0, fgkNRow-1, outde, ksite, nfit, nmiss, pchi2, (ii==fgkNiter-1 ? debugstreamer : 0x0), kTRUE);
 349     //PrintTrackParam(90020+ii, &trackOut);
 350   }
 351
 352   if(trackInOld)
 353     (*trackInOld)  = trackIn;
 354
 355   (*trackOutOld) = trackOut;
 356 }
 357
 358 void AliTPCCosmicUtils::CombinedFit(AliExternalTrackParam *trackPars[],  const AliTPCseed *seeds[],  Int_t &nfit, Int_t &nmiss, Double_t &pchi2, TTreeSRedirector *debugstreamer)
 359 {
 360   //
 361   //combined propagation
 362   //
 363
 364   //seen from 1/pt_fit - 1/pt_gen, for differnt de:
 365   //u+d+: combUp good, combLow bad
 366   //u+d-: combUp good, combLow good
 367   //u-d-: combUp bad,  combLow good
 368   //u-d+: combUp bad,  combLow bad
 369   const Double_t upde  = 1;
 370   const Double_t lowde = -1;
 371
 372   AliExternalTrackParam trackLow = *(trackPars[1]);
 373   AliExternalTrackParam trackUp;
 374
 375   for(Int_t ii=0; ii<fgkNiter; ii++){
 376     nmiss = 0;
 377
 378     //lower->upper
 379     trackUp = trackLow;
 380     SubCombined(&trackUp,  seeds, 1, 0, upde,  nfit, nmiss, pchi2);
 381
 382     nfit = 0;
 383     pchi2 = 0;
 384
 385     //upper->lower
 386     trackLow = trackUp;
 387     SubCombined(&trackLow, seeds, 0, 1, lowde, nfit, nmiss, pchi2, (ii==fgkNiter-1? debugstreamer : 0x0));
 388   }
 389
 390   *(trackPars[0]) = trackUp;
 391   *(trackPars[1]) = trackLow;
 392 }
 393
 394 void AliTPCCosmicUtils::SubCombined(AliExternalTrackParam *trackPar, const AliTPCseed *seeds[], const Int_t tk0, const Int_t tk1, const Double_t eloss, Int_t &nfit, Int_t &nmiss, Double_t &pchi2, TTreeSRedirector *debugstreamer)
 395 {
 396   //
 397   //sub-routine for combined propagation
 398   //
 399
 400   IniCov(trackPar, seeds[0]->GetNumberOfClusters()+seeds[1]->GetNumberOfClusters());
 401
 402   Int_t dtk=1;
 403   if(tk0>tk1)
 404     dtk=-1;
 405
 406   //reset counters
 407   Int_t ksite = -999;
 408   nfit = 0;
 409   nmiss = 0;
 410   pchi2 = 0;
 411
 412   //always nrow -> 1 -> nrow
 413   Int_t rowstart=fgkNRow-1;
 414   Int_t rowstop=0;
 415   for(Int_t itrk=tk0; dtk*itrk<=dtk*tk1; itrk+=dtk){
 416     if(itrk==tk1){
 417       Int_t tmprow = rowstart;
 418       rowstart = rowstop;
 419       rowstop = tmprow;
 420     }
 421     FitKernel(trackPar, seeds[itrk], rowstart, rowstop, eloss, ksite, nfit, nmiss, pchi2, debugstreamer, kFALSE);
 422   }
 423 }
 424
 425 void AliTPCCosmicUtils::FitKernel(AliExternalTrackParam *trackPar, const AliTPCseed *tseed, const Int_t rowstart, const Int_t rowstop, const Double_t eloss, Int_t &ksite, Int_t &nfit, Int_t &nmiss, Double_t &pchi2, TTreeSRedirector *debugstreamer, const Bool_t kinicov)
 426 {
 427   //
 428   //routine for propagation
 429   //
 430
 431   //only for SingleFit-->
 432   if(kinicov)
 433     IniCov(trackPar, tseed->GetNumberOfClusters());
 434   //<--
 435
 436   Int_t drow=1;
 437   if(rowstart>rowstop)
 438     drow = -1;
 439
 440   for(Int_t irow=rowstart; drow*irow<=drow*rowstop; irow+=drow){
 441     if(0){
 442       //only for momentum resolution investigation, dumm for usual usage
 443       if(eloss>0){
 444         if(irow%2==0)
 445           continue;
 446       }
 447       else{
 448         if(irow%2!=0)
 449           continue;
 450       }
 451     }
 452
 453     AliTPCclusterMI *cl=tseed->GetClusterPointer(irow);
 454     if (!cl) continue;
 455     if (cl->GetX()< fgkXMin) continue;
 456
 457     //if propagation not successful, trackPar is not changed
 458     AliExternalTrackParam tmppar(*trackPar);
 459
 460     //--------------------------- to cure large chi2 in z, 2011-05-11. Thandks to M. Ivanov!  ---------------------------
 461     const Int_t tmpsite = ((cl->GetDetector()%36)<18);
 462     if(tmpsite!=ksite && ksite!=-999){
 463       Double_t *clscov=(Double_t*)tmppar.GetCovariance();
 464       clscov[2]+=3*3;
 465     }
 466     ksite = tmpsite;
 467
 468     //--------------------------- rotate cluster position to trackPar position ---------------------------
 469     //DO NOT rotate trackPar, there is "bug" in trackparam::rotate!! stay in the initial one defined as alpha = ATan2(posy0,posx0)
 470     Float_t gxyz[3];
 471     cl->GetGlobalXYZ(gxyz);
 472     TVector3 ptogo(gxyz);
 473
 474     //printf("test %d : %f %f %f\n", irow, gxyz[0], gxyz[1], gxyz[2]);
 475
 476     const Double_t trkalpha = tmppar.GetAlpha();
 477     ptogo.RotateZ(-trkalpha);
 478     const Double_t xTogo = ptogo.X();
 479     const Double_t yTogo = ptogo.Y();
 480     const Double_t zTogo = ptogo.Z();
 481
 482     //--------------------------- propagate ---------------------------
 483     //AliTrackerBase::PropagateTrackToBxByBz Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
 484     const Double_t maxStep = 1;
 485     const Bool_t rotateTo = kFALSE;
 486     const Double_t maxSnp = 0.8;
 487
 488     //eloss critical only for combine fit!!!
 489     const Bool_t kpro = AliTrackerBase::PropagateTrackToBxByBz(&tmppar, xTogo, fgkMass, maxStep, rotateTo, maxSnp, eloss);
 490     if(!kpro){
 491       nmiss++;
 492       continue;
 493     }
 494
 495     //--------------------------- set error ---------------------------
 496     //necessary, help at mostly low p, and also high p
 497     Double_t cov[3]={0.01, 0., 0.01};
 498     AliTPCseed::GetError(cl, &tmppar, cov[0],cov[2]);
 499     cov[0]*=cov[0];
 500     cov[2]*=cov[2];
 501
 502     if(debugstreamer){
 503       TVector3 tmpgxyz(gxyz);
 504       (*debugstreamer)<<"ErrParam"<<
 505         "Cl.="<<cl<<
 506         "T.="<<&tmppar<<
 507         "covy="<<cov[0]<<
 508         "covz="<<cov[2]<<
 509         "clpos.="<<&ptogo<<
 510         "gpos.="<<&tmpgxyz<<
 511         "\n";
 512     }
 513
 514     //--------------------------- get chi2 and THEN update ---------------------------
 515     Double_t yz[2]={yTogo, zTogo};
 516
 517     pchi2 += tmppar.GetPredictedChi2(yz, cov);
 518
 519     tmppar.Update(yz,cov);
 520
 521     //--------------------------- save change ---------------------------
 522     (*trackPar) = tmppar;
 523
 524     nfit++;
 525   }
 526 }