During simulation: fill STU region w/ non null time sums
[u/mrichter/AliRoot.git] / 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 }