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