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 which can be used in standalone cases
  17 // especially as utils for AliTPCCosmicTrackfit
  18 //
  19 // detailed description can be found inside individual function
  20 //
  21 // grep "exitreport" in output log to check abnormal termination
  22 //
  23 //
  24 //  Xianguo Lu
  25 //  lu@physi.uni-heidelberg.de
  26 //  Xianguo.Lu@cern.ch
  27
  28 //
  29
  30 #include <TAxis.h>
  31 #include <TCanvas.h>
  32 #include <TGraph.h>
  33 #include <TTreeStream.h>
  34 #include <TVector3.h>
  35
  36 #include "AliESDtrack.h"
  37 #include "AliESDfriendTrack.h"
  38 #include "AliTPCseed.h"
  39 #include "AliTrackerBase.h"
  40 #include "AliTrackPointArray.h"
  41
  42 #include "AliTPCCosmicUtils.h"
  43
  44 Double_t AliTPCCosmicUtils::Point2LineDist(const TVector3 p0, const TVector3 l1, const TVector3 l2)
  45 {
  46   //
  47   //return distance of p0 to line (l2-l1)
  48   //
  49
  50   const TVector3 va = p0 - l1;
  51   const TVector3 vb = l2 - l1;
  52
  53   const TVector3 dd = va.Cross(vb);
  54
  55   return dd.Mag()/vb.Mag();
  56 }
  57
  58 Double_t AliTPCCosmicUtils::AngleInRange(Double_t phi)
  59 {
  60   //
  61   //get the phi value (in rad) in -pi ~ pi, so that fabs() works naiively
  62   //
  63   const Double_t pmin = -TMath::Pi();
  64   const Double_t pmax = pmin+TMath::TwoPi();
  65
  66   while(phi<pmin)
  67     phi+= TMath::TwoPi();
  68
  69   while(phi>=pmax)
  70     phi-=  TMath::TwoPi();
  71
  72   return phi;
  73 }
  74
  75 Bool_t AliTPCCosmicUtils::RotateSafe(AliExternalTrackParam *trackPar, const Double_t rawalpha)
  76 {
  77   //1. in AliExternalTrackParam::GetXYZ
  78   //r[0]=fX; r[1]=fP[0]; r[2]=fP[1];
  79   //return Local2GlobalPosition(r,fAlpha);
  80   //in AliVParticle::Local2GlobalMomentum
  81   //Double_t cs=TMath::Cos(alpha), sn=TMath::Sin(alpha);
  82   //Double_t r=TMath::Sqrt((1. - p[1])*(1. + p[1]));
  83   //p[0]=pt*(r*cs - p[1]*sn); p[1]=pt*(p[1]*cs + r*sn); p[2]=pt*p[2];
  84   //which is cos(phi_local=Snp+fAlpha) > 0 always.
  85   //2. also in Bool_t AliExternalTrackParam::Rotate(Double_t alpha)
  86   // Double_t sf=fP2, cf=TMath::Sqrt((1.- fP2)*(1.+fP2));
  87   //in AliExternalTrackParam::Set
  88   //mom.RotateZ(-fAlpha);
  89   //fP[2] = TMath::Sin(mom.Phi());
  90   //since only sin is used for mom.Phi(), that is assuming cos(mom.Phi())>0
  91
  92   //the range-changing in AliExternalTrackParam::Rotate not precise enough
  93   const Double_t aa = AngleInRange(rawalpha);
  94
  95   const Double_t a0 = trackPar->GetAlpha();
  96   const Double_t p2 = trackPar->GetParameter()[2];
  97
  98   //copied from AliExternalTrackParam::Rotate
  99   const Double_t ca=TMath::Cos(aa-a0), sa=TMath::Sin(aa-a0);
 100   const Double_t sf=p2, cf=TMath::Sqrt((1.- p2)*(1.+p2));
 101
 102   if((cf*ca+sf*sa)<0) {
 103     return kFALSE;
 104   }
 105
 106   return trackPar->Rotate(aa);
 107
 108   /*
 109   Double_t xyz[3], pxpypz[3];
 110   trackPar->GetXYZ(xyz);
 111   trackPar->GetPxPyPz(pxpypz);
 112   const TVector3 pos0(xyz);
 113   const TVector3 mom0(pxpypz);
 114   TVector3 pos1, mom1, dpos, dmom;
 115   const Double_t eps = 1e-6;
 116   AliExternalTrackParam tmppar = (*trackPar);
 117
 118   if(!tmppar.Rotate(aa)){
 119     return 0;
 120   }
 121
 122   tmppar.GetXYZ(xyz);
 123   tmppar.GetPxPyPz(pxpypz);
 124   pos1.SetXYZ(xyz[0], xyz[1], xyz[2]);
 125   mom1.SetXYZ(pxpypz[0], pxpypz[1], pxpypz[2]);
 126   dpos = pos1-pos0;
 127   dmom = mom1-mom0;
 128   if(dpos.Mag()<eps && dmom.Mag()<eps){
 129     (*trackPar)=tmppar;
 130     return 1;
 131   }
 132   else
 133     return 0;
 134   */
 135
 136   /*
 137   //flip
 138   tmppar = (*trackPar);
 139   if(!tmppar.Rotate(aa+TMath::Pi())){
 140     return 0;
 141   }
 142
 143   tmppar.GetXYZ(xyz);
 144   tmppar.GetPxPyPz(pxpypz);
 145   pos1.SetXYZ(xyz[0], xyz[1], xyz[2]);
 146   mom1.SetXYZ(pxpypz[0], pxpypz[1], pxpypz[2]);
 147   dpos = pos1-pos0;
 148   dmom = mom1-mom0;
 149   if(dpos.Mag()<eps && dmom.Mag()<eps){
 150     (*trackPar)=tmppar;
 151     return -1;
 152   }
 153
 154   return 0;
 155   */
 156 }
 157
 158 void AliTPCCosmicUtils::PrintTrackParam(const Int_t id, const AliExternalTrackParam * trackpar, const char *tag)
 159 {
 160   //
 161   //print out TrackParam
 162   //
 163   Double_t xyz[3]={-999,-999,-999};
 164   trackpar->GetXYZ(xyz);
 165   Double_t pxpypz[3]={-999,-999,-999};
 166   trackpar->GetPxPyPz(pxpypz);
 167
 168   printf("PrintTrackPar %2s %3d : %6.1f %5.1f : %6.6f %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], trackpar->Phi(), pxpypz[0], pxpypz[1], pxpypz[2], trackpar->Charge(), trackpar->Pt(), trackpar->P(), trackpar->GetAlpha()*TMath::RadToDeg());
 169 }
 170
 171 void AliTPCCosmicUtils::DrawTracks(AliESDtrack *esdtrks[], const TString tag, const TString outputformat)
 172 {
 173   //
 174   //draw esdtracks
 175   //
 176   const AliTPCseed *seeds[]={GetTPCseed(esdtrks[0]), GetTPCseed(esdtrks[1])};
 177   DrawSeeds(seeds, tag, outputformat);
 178 }
 179
 180
 181 void AliTPCCosmicUtils::DrawSeeds(const AliTPCseed * seeds[], const TString tag, const TString outputformat)
 182 {
 183   //
 184   //draw seed and output to file
 185   //
 186
 187   if(!seeds[0] || !seeds[1])
 188     return;
 189
 190   TGraph *grsyx[]={new TGraph, new TGraph};
 191   TGraph *grsyz[]={new TGraph, new TGraph};
 192
 193   for(Int_t itrk=0; itrk<2; itrk++){
 194     if(!seeds[itrk])
 195       continue;
 196
 197     Int_t ncl = 0;
 198     for(Int_t irow=0; irow<fgkNRow; irow++){
 199       const AliTPCclusterMI * cls = seeds[itrk]->GetClusterPointer(irow);
 200       if(!cls)
 201         continue;
 202
 203       Float_t xyz[3]={-999,-999,-999};
 204       cls->GetGlobalXYZ(xyz);
 205       //printf("Test %f %f %f \n", xyz[0], xyz[1], xyz[2]);
 206       grsyx[itrk]->SetPoint(ncl, xyz[0], xyz[1]);
 207       grsyz[itrk]->SetPoint(ncl, xyz[2], xyz[1]);
 208       ncl++;
 209     }
 210   }
 211
 212   grsyx[0]->SetTitle(tag);
 213   grsyx[0]->SetMaximum(250);
 214   grsyx[0]->SetMinimum(-250);
 215   grsyx[0]->GetXaxis()->SetLimits(-250,250);
 216   grsyx[0]->GetXaxis()->SetTitle("X (cm)");
 217   grsyx[0]->GetYaxis()->SetTitle("Y (cm)");
 218
 219   grsyx[0]->SetMarkerStyle(20);
 220   grsyx[0]->SetMarkerColor(kRed);
 221   grsyx[1]->SetMarkerStyle(22);
 222   grsyx[1]->SetMarkerColor(kBlue);
 223
 224   grsyz[0]->SetTitle(tag);
 225   grsyz[0]->SetMaximum(250);
 226   grsyz[0]->SetMinimum(-250);
 227   grsyz[0]->GetXaxis()->SetLimits(-250,250);
 228   grsyz[0]->GetXaxis()->SetTitle("Z (cm)");
 229   grsyz[0]->GetYaxis()->SetTitle("Y (cm)");
 230
 231   grsyz[0]->SetMarkerStyle(20);
 232   grsyz[0]->SetMarkerColor(kRed);
 233   grsyz[1]->SetMarkerStyle(22);
 234   grsyz[1]->SetMarkerColor(kBlue);
 235
 236   TCanvas *cc=new TCanvas("cc","",1200,600);
 237   cc->Divide(2);
 238   cc->cd(1);
 239   for(int itrk=0; itrk<2; itrk++){
 240     grsyx[itrk]->SetMarkerSize(1);
 241     grsyx[itrk]->Draw(itrk?"lp same":"alp");
 242   }
 243   cc->cd(2);
 244   for(int itrk=0; itrk<2; itrk++){
 245     grsyz[itrk]->SetMarkerSize(1);
 246     grsyz[itrk]->Draw(itrk?"lp same":"alp");
 247   }
 248
 249   gErrorIgnoreLevel = 1001;
 250   cc->Print(Form("drawTrack%s.%s", tag.Data(), outputformat.Data()));
 251
 252   for(Int_t ii=0; ii<2; ii++){
 253     delete grsyx[ii];
 254     delete grsyz[ii];
 255   }
 256   delete cc;
 257 }
 258
 259 //================================================================================================
 260 AliTPCseed * AliTPCCosmicUtils::GetTPCseed(const AliESDtrack *esdtrack)
 261 {
 262   //
 263   //Get TPC seeds from ESDfriendTrack
 264   //
 265
 266   AliESDfriendTrack * friendtrk = (AliESDfriendTrack *) esdtrack->GetFriendTrack();
 267   if(!friendtrk){
 268     printf("exitreport AliTPCCosmicUtils::GetTPCseed no friend track!\n");
 269     exit(1);
 270   }
 271
 272   TObject *calibObject=0x0;
 273   AliTPCseed *tseed = 0x0;
 274   for(Int_t l=0; (calibObject=friendtrk->GetCalibObject(l)); l++) {
 275     if( (tseed=dynamic_cast<AliTPCseed*>(calibObject)) )
 276       break;
 277   }
 278   return tseed;
 279 }
 280
 281 //================================================================================================
 282 AliExternalTrackParam *AliTPCCosmicUtils::MakeSeed(const AliTPCseed *tseed)
 283 {
 284   //
 285   //make seed for propagation of TrackParam, using np = 3 outer clusters (separated by deltancls clusters) in TPCseed
 286   //
 287
 288   const Int_t rowstart = fgkNRow-1;
 289   const Int_t rowstop = 0;
 290   const Int_t drow = -1;
 291
 292   //---
 293   const Int_t np = 3;
 294   AliTrackPoint *tpos[np];
 295   for(Int_t ii=0; ii<np; ii++)
 296     tpos[ii] = 0x0;
 297
 298   const Float_t cov[6]={0,0,0, 0.01*0.01 ,0, 0.01*0.01};
 299
 300   //---
 301   Int_t npos = 0;
 302   Int_t icl = 0;
 303   const Int_t deltancls = fgkNclsMin/2-1;
 304   Int_t oldcl = -deltancls;
 305
 306   for(Int_t irow=rowstart; drow*irow<=drow*rowstop; irow+=drow){
 307     AliTPCclusterMI *cls = tseed->GetClusterPointer(irow);
 308     if(!cls) {
 309       continue;
 310     }
 311
 312     if( icl == (oldcl+deltancls) ){
 313       Float_t txyz[3];
 314       cls->GetGlobalXYZ(txyz);
 315       tpos[npos++] = new AliTrackPoint(txyz, cov, 0);
 316       //printf("------ %d %f %f %f\n", npos, txyz[0], txyz[1], txyz[2]);
 317
 318       oldcl = icl;
 319       if(npos==np) break;
 320     }
 321     icl++;
 322   }
 323   if(npos!=np){
 324     for(Int_t ii=0; ii<npos; ii++)
 325       delete tpos[ii];
 326
 327     return 0x0;
 328   }
 329
 330   AliExternalTrackParam * trackparam = AliTrackerBase::MakeSeed(*(tpos[0]), *(tpos[1]), *(tpos[2]));
 331   if(!trackparam || trackparam->Pt()==0){
 332     for(Int_t ii=0; ii<npos; ii++)
 333       delete tpos[ii];
 334     delete trackparam;
 335
 336     return 0x0;
 337   }
 338
 339   //------
 340
 341   Double_t sxyz[3]={-999,-999,-999}, spxpypz[3]={-999,-999,-999};
 342   trackparam->GetXYZ(sxyz);
 343   trackparam->GetPxPyPz(spxpypz);
 344   Double_t scov[21];
 345   Int_t sign = trackparam->Charge();
 346
 347   //reset covariance matrix -- necessary, otherwise bad fitting result: trackparam->GetCovariance()[ii] has problem: some are strange values
 348   for(Int_t ii=0; ii<21; ii++) {
 349     scov[ii]=0;
 350   }
 351   trackparam->Set(sxyz, spxpypz, scov, sign);
 352
 353   for(Int_t ii=0; ii<np; ii++)
 354     delete tpos[ii];
 355
 356   return trackparam;
 357 }
 358
 359 //================================================================================================
 360 void AliTPCCosmicUtils::IniCov(AliExternalTrackParam *trackPar, const Double_t ncl)
 361 {
 362   //
 363   //initialize covariance matrix
 364   //
 365
 366   const Double_t ksigma=5.;
 367   Double_t acov[16];
 368   for (Int_t i=0;i<15;i++)
 369     acov[i]=0;
 370
 371   acov[0]=ksigma*ksigma;
 372   acov[2]=ksigma*ksigma;
 373   acov[5]=ksigma*ksigma;
 374   acov[9]=ksigma*ksigma;
 375   acov[14]=0.2*0.2;
 376
 377   acov[5] = ksigma*ksigma/(ncl*ncl);
 378   acov[9] = ksigma*ksigma/(ncl*ncl);
 379
 380   const Double_t resetcov = 4;
 381
 382   trackPar->ResetCovariance(resetcov);
 383   //the following helps a lot!!
 384   trackPar->AddCovariance(acov);
 385 }
 386
 387 void AliTPCCosmicUtils::SingleFit(AliExternalTrackParam * trackInOld, AliExternalTrackParam * trackOutOld, const AliTPCseed *tseed, const Bool_t kinward,  const Int_t rowstartshift, const Int_t rowstep, const Double_t xmin, const Double_t xmax, Int_t &nfit, Int_t &nmiss, Double_t &pchi2, Double_t &lfit, TTreeSRedirector *debugstreamer)
 388 {
 389   //
 390   //fit single track
 391   //
 392   //kinward is the true geometry of the track. Incomming track: 1; outgoing track: 0
 393   const Double_t inde  = kinward? -1 :  1;
 394   const Double_t outde = kinward?  1 : -1;
 395
 396   //PrintTrackParam(9000, trackOutOld);
 397
 398   AliExternalTrackParam trackOut = *trackOutOld;
 399   AliExternalTrackParam trackIn;
 400   Int_t ksite = -999;
 401
 402   //nmiss is from the 2 FitKernel of the last iteration
 403   //nfit, pchi2 is from the last FitKernel of the last iteration
 404
 405   Int_t rowouter = fgkNRow-1-rowstartshift;
 406
 407   //so that when reversed, the same rows are read! important!!!
 408   Int_t rowinner = rowouter - rowouter/rowstep * rowstep;
 409
 410   TVector3 gposStart;
 411   TVector3 gposStop;
 412   TVector3 dpos;
 413   lfit = 0;
 414
 415   for(Int_t ii=0; ii<fgkNiter; ii++){
 416     nmiss = 0;
 417
 418     gposStart.SetXYZ(-999,-999,-999);
 419     ksite = -999;
 420     trackIn  = trackOut;
 421     FitKernel(&trackIn,  tseed, rowouter, rowinner, -rowstep, xmin, xmax, inde , ksite, nfit, nmiss, pchi2, gposStart, gposStop, 0x0, kTRUE);
 422     //PrintTrackParam(9010+ii, &trackIn);
 423
 424     dpos = gposStart-gposStop;
 425     lfit += dpos.Pt();
 426
 427     //---------------------------
 428
 429     nfit = 0;
 430     pchi2 = 0;
 431
 432     gposStart.SetXYZ(-999,-999,-999);
 433     ksite = -999;
 434     trackOut = trackIn;
 435     FitKernel(&trackOut, tseed, rowinner, rowouter,  rowstep, xmin, xmax, outde, ksite, nfit, nmiss, pchi2, gposStart, gposStop, (ii==fgkNiter-1 ? debugstreamer : 0x0), kTRUE);
 436     //PrintTrackParam(90020+ii, &trackOut);
 437
 438     dpos = gposStart-gposStop;
 439     lfit += dpos.Pt();
 440   }
 441
 442   lfit /= 2.*fgkNiter;
 443
 444   if(trackInOld)
 445     (*trackInOld)  = trackIn;
 446
 447   (*trackOutOld) = trackOut;
 448 }
 449
 450 void AliTPCCosmicUtils::CombinedFit(AliExternalTrackParam *trackPars[],  const AliTPCseed *seeds[],  const Int_t rowstartshift, const Int_t rowstep, const Double_t xmin, const Double_t xmax, Int_t &nfit, Int_t &nmiss, Double_t &pchi2, Double_t &lfit, Double_t &vtxD, Double_t &vtxZ, TTreeSRedirector *debugstreamer)
 451 {
 452   //
 453   //combined propagation
 454   //
 455
 456   //seen from 1/pt_fit - 1/pt_gen, for differnt de:
 457   //u+d+: combUp good, combLow bad
 458   //u+d-: combUp good, combLow good
 459   //u-d-: combUp bad,  combLow good
 460   //u-d+: combUp bad,  combLow bad
 461   const Double_t upde  = 1;
 462   const Double_t lowde = -1;
 463
 464   AliExternalTrackParam trackLow = *(trackPars[1]);
 465   AliExternalTrackParam trackUp;
 466
 467   for(Int_t ii=0; ii<fgkNiter; ii++){
 468     lfit = 0;
 469     vtxD = 0;
 470     vtxZ = 0;
 471     Double_t tmpl = -999, tmpd = -999, tmpz = -999;
 472
 473     //lower->upper
 474     trackUp = trackLow;
 475     SubCombined(&trackUp,  seeds, 1, 0, rowstartshift, rowstep, xmin, xmax, upde,  nfit, nmiss, pchi2, tmpl, tmpd, tmpz);
 476     lfit += tmpl;
 477     vtxD += tmpd;
 478     vtxZ += tmpz;
 479
 480     //upper->lower
 481     trackLow = trackUp;
 482     SubCombined(&trackLow, seeds, 0, 1, rowstartshift, rowstep, xmin, xmax, lowde, nfit, nmiss, pchi2, tmpl, tmpd, tmpz, (ii==fgkNiter-1? debugstreamer : 0x0));
 483     lfit += tmpl;
 484     vtxD += tmpd;
 485     vtxZ += tmpz;
 486   }
 487
 488   *(trackPars[0]) = trackUp;
 489   *(trackPars[1]) = trackLow;
 490
 491   //only last iteration used
 492   lfit /= 2;
 493   vtxD /= 2;
 494   vtxZ /= 2;
 495 }
 496
 497 void AliTPCCosmicUtils::SubCombined(AliExternalTrackParam *trackPar, const AliTPCseed *seeds[], const Int_t tk0, const Int_t tk1, const Int_t rowstartshift, const Int_t rowstep, const Double_t xmin, const Double_t xmax, const Double_t eloss, Int_t &nfit, Int_t &nmiss, Double_t &pchi2, Double_t &lfit, Double_t &vtxD, Double_t &vtxZ, TTreeSRedirector *debugstreamer)
 498 {
 499   //
 500   //sub-routine for combined propagation
 501   //
 502
 503   IniCov(trackPar, seeds[0]->GetNumberOfClusters()+seeds[1]->GetNumberOfClusters());
 504
 505   Int_t dtk=1;
 506   if(tk0>tk1)
 507     dtk=-1;
 508
 509   //reset counters
 510   Int_t ksite = -999;
 511   nfit = 0;
 512   nmiss = 0;
 513   pchi2 = 0;
 514   vtxD = -1e10;
 515   vtxZ = -1e10;
 516
 517   //always nrow -> 1 -> nrow
 518   Int_t rowstart= fgkNRow-1-rowstartshift;
 519
 520   //so that when reversed, the same rows are read! important!!!
 521   Int_t rowstop = rowstart - rowstart/rowstep * rowstep;
 522   Int_t drow = -rowstep;
 523
 524   TVector3 gposStart(-999,-999,-999);
 525   TVector3 gposStop(-999,-999,-999);
 526
 527   for(Int_t itrk=tk0; dtk*itrk<=dtk*tk1; itrk+=dtk){
 528     if(itrk==tk1){
 529       Int_t tmprow = rowstart;
 530       rowstart = rowstop;
 531       rowstop = tmprow;
 532       drow *= -1;
 533     }
 534
 535     FitKernel(trackPar, seeds[itrk], rowstart, rowstop, drow, xmin, xmax, eloss, ksite, nfit, nmiss, pchi2, gposStart, gposStop, debugstreamer, kFALSE);
 536
 537     //get the impact parameters at the end of the first propagation (X=0)
 538     if(itrk==tk0){
 539       AliExternalTrackParam vertex(*trackPar);
 540       const Double_t maxStep = 1;
 541       const Bool_t rotateTo = kFALSE;
 542       const Double_t maxSnp = 0.8;
 543       if(AliTrackerBase::PropagateTrackToBxByBz(&vertex, 0, fgkMass, maxStep, rotateTo, maxSnp, eloss)){
 544         vtxD = TMath::Abs(vertex.GetParameter()[0]);
 545         vtxZ = vertex.GetParameter()[1];
 546       }
 547     }
 548   }
 549
 550   TVector3 dpos = gposStart-gposStop;
 551   lfit = dpos.Pt();
 552 }
 553
 554 void AliTPCCosmicUtils::FitKernel(AliExternalTrackParam *trackPar, const AliTPCseed *tseed, const Int_t rowstart, const Int_t rowstop, const Int_t drow, const Double_t xmin, const Double_t xmax, const Double_t eloss, Int_t &ksite, Int_t &nfit, Int_t &nmiss, Double_t &pchi2, TVector3 &gposStart, TVector3 &gposStop, TTreeSRedirector *debugstreamer, const Bool_t kinicov)
 555 {
 556   //
 557   //routine for propagation
 558   //
 559
 560   //only for SingleFit-->
 561   if(kinicov)
 562     IniCov(trackPar, tseed->GetNumberOfClusters());
 563   //<--
 564
 565   Int_t checkstop = -999;
 566
 567   for(Int_t irow=rowstart; drow*irow<=drow*rowstop; irow+=drow){
 568     checkstop = irow;
 569     //printf("test irow %d\n", irow);
 570
 571     AliTPCclusterMI *cl=tseed->GetClusterPointer(irow);
 572     if (!cl) continue;
 573     if (cl->GetX()< fgkXMin) continue;
 574
 575     //cut on cluster X (i.e. r) to specify leverarm
 576     if(cl->GetX()< xmin){
 577       continue;
 578     }
 579     if(cl->GetX()> xmax){
 580       continue;
 581     }
 582     //if propagation not successful, trackPar is not changed
 583     AliExternalTrackParam tmppar(*trackPar);
 584
 585     //--------------------------- to cure large chi2 in z, 2011-05-11. Thandks to M. Ivanov!  ---------------------------
 586     const Int_t tmpsite = ((cl->GetDetector()%36)<18);
 587     if(tmpsite!=ksite && ksite!=-999){
 588       Double_t *clscov=(Double_t*)tmppar.GetCovariance();
 589       clscov[2]+=3*3;
 590     }
 591     ksite = tmpsite;
 592
 593     //--------------------------- rotate cluster position to trackPar position ---------------------------
 594     //DO NOT rotate trackPar, there is "bug" in trackparam::rotate!! stay in the initial one defined as alpha = ATan2(posy0,posx0)
 595     Float_t gxyz[3];
 596     cl->GetGlobalXYZ(gxyz);
 597
 598     const TVector3 gptmp(gxyz[0], gxyz[1], gxyz[2]);
 599
 600     TVector3 ptogo(gxyz);
 601
 602     //printf("test %d : %f %f %f\n", irow, gxyz[0], gxyz[1], gxyz[2]);
 603
 604     const Double_t trkalpha = tmppar.GetAlpha();
 605     ptogo.RotateZ(-trkalpha);
 606     const Double_t xTogo = ptogo.X();
 607     const Double_t yTogo = ptogo.Y();
 608     const Double_t zTogo = ptogo.Z();
 609
 610     //--------------------------- propagate ---------------------------
 611     //AliTrackerBase::PropagateTrackToBxByBz Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
 612     const Double_t maxStep = 1;
 613     const Bool_t rotateTo = kFALSE;
 614     const Double_t maxSnp = 0.8;
 615
 616     //eloss critical only for combine fit!!!
 617     const Bool_t kpro = AliTrackerBase::PropagateTrackToBxByBz(&tmppar, xTogo, fgkMass, maxStep, rotateTo, maxSnp, eloss);
 618     if(!kpro){
 619       nmiss++;
 620       continue;
 621     }
 622
 623     //--------------------------- set error ---------------------------
 624     //necessary, help at mostly low p, and also high p
 625     Double_t cov[3]={0.01, 0., 0.01};
 626     AliTPCseed::GetError(cl, &tmppar, cov[0],cov[2]);
 627     cov[0]*=cov[0];
 628     cov[2]*=cov[2];
 629
 630     if(debugstreamer){
 631       TVector3 tmpgxyz(gxyz);
 632       (*debugstreamer)<<"ErrParam"<<
 633         "Cl.="<<cl<<
 634         "T.="<<&tmppar<<
 635         "covy="<<cov[0]<<
 636         "covz="<<cov[2]<<
 637         "clpos.="<<&ptogo<<
 638         "gpos.="<<&tmpgxyz<<
 639         "\n";
 640     }
 641
 642     //--------------------------- get chi2 and THEN update ---------------------------
 643     Double_t yz[2]={yTogo, zTogo};
 644
 645     pchi2 += tmppar.GetPredictedChi2(yz, cov);
 646
 647     tmppar.Update(yz,cov);
 648
 649     //--------------------------- save change ---------------------------
 650     (*trackPar) = tmppar;
 651
 652     nfit++;
 653
 654     gposStop = gptmp;
 655     if(gposStart.X()<-998)
 656       gposStart = gptmp;
 657   }
 658
 659   //to make sure rowstart and rowstop are the actual ending
 660   if(checkstop != rowstop){
 661     printf("exitreport AliTPCCosmicUtils::FitKernel wrong rowstart, stop, drow!! %d %d %d\n", rowstart, rowstop, drow);
 662     exit(1);
 663   }
 664 }