TPC/Upgrade/macros/testRec.C

   1 Int_t GetTimeAtVertex(Float_t &tVtx, Float_t &x, AliToyMCTrack *tr, Int_t clsType=0, Int_t seedRow=140, Int_t seedDist=10, Int_t correctionType=0);
   2 void SetTrackPointFromCluster(AliTPCclusterMI *cl, AliTrackPoint &p);
   3 void ClusterToSpacePoint(AliTPCclusterMI *cl, Float_t xyz[3]);
   4 void InitSpaceCharge(TTree *t=0x0);
   5 /*
   6
   7 root.exe -l $ALICE_ROOT/TPC/Upgrade/macros/{loadlibs.C,ConfigOCDB.C}
   8 .x $ALICE_ROOT/TPC/Upgrade/macros/testRec.C
   9
  10 */
  11
  12 AliTPCParam *fTPCParam=AliTPCcalibDB::Instance()->GetParameters();
  13 AliTPCSpaceCharge3D *fSpaceCharge=0x0;
  14 TTreeSRedirector *fStreamer=0x0;
  15
  16 Int_t fEvent=-1;
  17 Int_t fTrack=-1;
  18 Float_t fT0event=-1;
  19 Float_t fZevent=-1;
  20
  21 void testRec(const char* filename="toyMC.root", Int_t nmaxEv=-1)
  22 {
  23   //
  24   //
  25   //
  26
  27   TFile f(filename);
  28   if (!f.IsOpen() || f.IsZombie()) {
  29     printf("ERROR: couldn't open the file '%s'\n", filename);
  30     return;
  31   }
  32
  33   TTree *t=(TTree*)f.Get("toyMCtree");
  34   if (!t) {
  35     printf("ERROR: couldn't read the 'toyMCtree' from file '%s'\n", filename);
  36     return;
  37   }
  38
  39   AliToyMCEvent *ev=0x0;
  40   t->SetBranchAddress("event",&ev);
  41
  42   // read spacecharge from the Userinfo ot the tree
  43   InitSpaceCharge(t);
  44
  45   TString debugName=filename;
  46   debugName.ReplaceAll(".root","");
  47   debugName.Append(".debug.root");
  48
  49   gSystem->Exec(Form("test -f %s && rm %s", debugName.Data(), debugName.Data()));
  50   if (!fStreamer) fStreamer=new TTreeSRedirector(debugName.Data());
  51
  52   gROOT->cd();
  53
  54 //   const Double_t kDriftVel = fTPCParam->GetDriftV()/1000000;
  55   const Double_t kDriftVel = fTPCParam->GetDriftV();
  56   const Double_t kMaxZ0=fTPCParam->GetZLength();
  57
  58   TH1F *h0=new TH1F("h0","h0",1000,0,0);
  59   TH1F *hX=new TH1F("hX","hX",1000,0,0);
  60   TH1F *h1=new TH1F("h1","h1",1000,0,0);
  61   TH1I *hcount0=new TH1I("count0","Failed extrapolation1",5,0,5);
  62   TH1I *hcount1=new TH1I("count1","Failed extrapolation2",5,0,5);
  63
  64   Int_t maxev=t->GetEntries();
  65   if (nmaxEv>0&&nmaxEv<maxev) maxev=nmaxEv;
  66
  67   for (Int_t iev=0; iev<maxev; ++iev){
  68     t->GetEvent(iev);
  69     fEvent=iev;
  70     for (Int_t itr=0; itr<ev->GetNumberOfTracks(); ++itr){
  71       printf("==============  Processing Track %6d\n",itr);
  72       fTrack=itr;
  73       fT0event = ev->GetT0();
  74       fZevent  = ev->GetZ();
  75
  76       //Float_t &tVtx, Float_t &x, AliToyMCTrack *tr,
  77       //  Int_t clsType=0, Int_t seedRow=140, Int_t seedDist=10, Int_t correctionType=0
  78       // correctionType: 0 none, 1 center, 2 mean tan,
  79       //                 3 full from seed (iterative), 4 ideal (real z-Position)
  80       AliToyMCTrack *tr=ev->GetTrack(itr);
  81       tr->SetUniqueID(itr);
  82       Float_t tVtx0=0;
  83       Float_t xmin=0;
  84       Int_t ret0=GetTimeAtVertex(tVtx0,xmin,tr);
  85       hX->Fill(xmin);
  86       //fully distorted
  87       Float_t tVtx1=0;
  88       Int_t ret1=GetTimeAtVertex(tVtx1,xmin,tr,1);// seeding at the outside
  89       GetTimeAtVertex(tVtx1,xmin,tr,1,70); // seeding in the center
  90       GetTimeAtVertex(tVtx1,xmin,tr,1,0); // seeding at the inside
  91       //correction at tpc center
  92       GetTimeAtVertex(tVtx1,xmin,tr,1,140, 10, 1);
  93       //correction with mean tan theta
  94       GetTimeAtVertex(tVtx1,xmin,tr,1,140, 10, 2);
  95       //correction with ideal z
  96       GetTimeAtVertex(tVtx1,xmin,tr,1,140, 10, 4);
  97
  98       hcount0->Fill(ret0);
  99       hcount1->Fill(ret1);
 100       if (ret0==0) {
 101         h0->Fill(tVtx0);
 102       }
 103
 104       if (ret1==0) {
 105         h1->Fill(tVtx1);
 106       }
 107 //       printf("TVtx: %f, %f\n",tVtx0,0);
 108     }
 109   }
 110
 111   TCanvas *c=(TCanvas*)gROOT->GetListOfCanvases()->FindObject("cOutput");
 112   if (!c) c=new TCanvas("cOutput","Results");
 113   c->Clear();
 114   c->Divide(2,2);
 115   c->cd(1);
 116   h0->Draw();
 117   h1->SetLineColor(kRed);
 118   h1->Draw("same");
 119   c->cd(2);
 120   hcount0->Draw();
 121   hcount1->SetLineColor(kRed);
 122   hcount1->Draw("same");
 123   c->cd(3);
 124   hX->Draw();
 125
 126   delete fStreamer;
 127   fStreamer=0x0;
 128 }
 129
 130 //____________________________________________________________________________
 131 Float_t GetTimeAtVertex(Float_t &tVtx,  Float_t &x, AliToyMCTrack *tr, Int_t clsType, Int_t seedRow, Int_t seedDist, Int_t correctionType)
 132 {
 133   //
 134   // clsType:    0 undistorted; 1 distorted
 135   // seedRow:    seeding row
 136   // seedDist:   distance of seeding points
 137   // correctionType: 0 none, 1 center, 2 mean tan,
 138   //                 3 full from seed (iterative), 4 ideal (real z-Position)
 139   //
 140
 141   // seed point informaion
 142   AliTrackPoint    seedPoint[3];
 143   AliTPCclusterMI *seedCluster[3]={0x0,0x0,0x0};
 144
 145   // number of clusters to loop over
 146   const Int_t ncls=(clsType==0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
 147
 148   UChar_t nextSeedRow=seedRow;
 149   Int_t   seed=0;
 150
 151   //assumes sorted clusters
 152   for (Int_t icl=0;icl<ncls;++icl) {
 153     AliTPCclusterMI *cl=tr->GetSpacePoint(icl);
 154     if (clsType==1) cl=tr->GetDistortedSpacePoint(icl);
 155     if (!cl) continue;
 156     // use row in sector
 157     const UChar_t row=cl->GetRow() + 63*(cl->GetDetector()>35);
 158     // skip clusters without proper pad row
 159     if (row>200) continue;
 160
 161     //check seeding row
 162     // if we are in the last row and still miss a seed we use the last row
 163     //   even if the row spacing will not be equal
 164     if (row>=nextSeedRow || icl==ncls-1){
 165       seedCluster[seed]=cl;
 166       SetTrackPointFromCluster(cl, seedPoint[seed]);
 167 //       printf("\nSeed point %d: %d, %d, %.2f, %.2f, %.2f, %.2f, %.2f\n",seed, cl->GetDetector(), row, seedPoint[seed].GetX(),seedPoint[seed].GetY(),seedPoint[seed].GetZ(), seedPoint[seed].GetAngle(), ((cl->GetDetector()%18)*20.+10.)/180.*TMath::Pi());
 168       ++seed;
 169       nextSeedRow+=seedDist;
 170
 171       if (seed==3) break;
 172     }
 173   }
 174
 175   // check we really have 3 seeds
 176   if (seed!=3) {
 177
 178     // debug output for failed seeding
 179     (*fStreamer) << "TracksFailSeed" <<
 180       "iev="      << fEvent <<
 181       "t0="       << fT0event <<
 182       "z0="       << fZevent <<
 183       "itrack="   << fTrack <<
 184       "clsType="  << clsType <<
 185       "seedRow="  << seedRow <<
 186       "seedDist=" << seedDist <<
 187       "corrType=" << correctionType <<
 188       "track.="   << tr    <<
 189       "\n";
 190
 191     printf("Seeding failed for parameters %d, %d\n",seedRow,seedDist, seed);
 192     return 1;
 193   }
 194
 195   // do cluster correction and
 196   // assign the cluster abs time as z component to all seeds
 197   for (Int_t iseed=0; iseed<3; ++iseed) {
 198     Float_t xyz[3]={0,0,0};
 199     seedPoint[iseed].GetXYZ(xyz);
 200
 201     Int_t sector=seedCluster[iseed]->GetDetector();
 202     Int_t sign=1-2*((sector/18)%2);
 203
 204     if (clsType && correctionType) {
 205       if (correctionType==1) xyz[2]=125.;
 206       if (correctionType==2) xyz[2]=TMath::Tan(45./2.*TMath::DegToRad())*xyz[1]*sign;
 207 //       if (correctionType==3) xyz[2]=125.;
 208       if (correctionType==4) xyz[2]=seedCluster[iseed]->GetZ();
 209
 210       fSpaceCharge->CorrectPoint(xyz, seedCluster[iseed]->GetDetector());
 211     }
 212
 213     // set different sign for c-Side (only for testing: makes half of the times negative)
 214     // xyz[2]=seedCluster[iseed]->GetTimeBin() * sign;
 215
 216     // correct with track z (only for testing: not possible in exp?)
 217     // xyz[2]=seedCluster[iseed]->GetTimeBin() + sign * tr->GetZ()/(fTPCParam->GetDriftV());
 218
 219     // no correction (default)
 220     xyz[2]=seedCluster[iseed]->GetTimeBin();
 221
 222     seedPoint[iseed].SetXYZ(xyz);
 223   }
 224
 225   // create seed and Propagate to r=0;
 226
 227   const Double_t kMaxSnp = 0.85;
 228   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 229
 230   AliExternalTrackParam *track = 0x0;
 231   track = AliTrackerBase::MakeSeed(seedPoint[0], seedPoint[1], seedPoint[2]);
 232   track->ResetCovariance(10);
 233
 234   printf("orig:  %.2f, %.2f, %.2f, %.2f, %.2f, %.2f\n",
 235          tr->GetParameter()[0],tr->GetParameter()[1],tr->GetParameter()[2],
 236          tr->GetParameter()[3],tr->GetParameter()[4],tr->GetAlpha());
 237
 238   printf("seed:  %.2f, %.2f, %.2f, %.2f, %.2f, %.2f\n",
 239          track->GetParameter()[0],track->GetParameter()[1],track->GetParameter()[2],
 240          track->GetParameter()[3],track->GetParameter()[4],track->GetAlpha());
 241
 242   //   printf("Track: %.2f, %.2f, %.2f, %.2f, %.2f\n",track->GetX(),track->GetY(),track->GetZ(), track->GetAlpha(),track->Phi());
 243   AliExternalTrackParam pInit(*track);
 244
 245   // NOTE:
 246   // when propagating with the time binwe need to switch off the material correction
 247   // otherwise we will be quite off ...
 248   //
 249   AliTrackerBase::PropagateTrackTo(track,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
 250   if (TMath::Abs(track->GetX())>3) {
 251     printf("Could not propagate track to 0, %.2f, %.2f, %.2f\n",track->GetX(),track->GetAlpha(),track->Pt());
 252 //     return 2;
 253   }
 254 //   printf("Track2: %.2f, %.2f, %.2f, %.2f\n",track->GetX(),track->GetY(),track->GetZ(), track->GetAlpha());
 255
 256   // simple linear fit
 257   TGraph gr;
 258   for (Int_t i=0; i<3; ++i)
 259     gr.SetPoint(gr.GetN(),seedPoint[i].GetX(),seedPoint[i].GetZ());
 260 //   gr.Print();
 261   TF1 fpol1("fpol1","pol1");
 262   gr.Fit(&fpol1,"QN");
 263   Float_t fitT0=fpol1.Eval(0);
 264 //   fpol1.Print();
 265   AliExternalTrackParam pOrig(*tr);
 266   (*fStreamer) << "Tracks" <<
 267     "iev="      << fEvent <<
 268     "t0="       << fT0event <<
 269     "z0="       << fZevent <<
 270     "itrack="   << fTrack <<
 271     "clsType="  << clsType <<
 272     "seedRow="  << seedRow <<
 273     "seedDist=" << seedDist <<
 274     "corrType=" << correctionType <<
 275     "track.="   << tr    <<
 276     "seed.="    << track <<
 277     "seedI.="   << &pInit <<
 278 //     "seedcl0.=" << seedCluster[0] <<
 279 //     "seedcl1.=" << seedCluster[1] <<
 280 //     "seedcl2.=" << seedCluster[2] <<
 281 //     "seedp0.="  << &seedPoint[0] <<
 282 //     "seedp1.="  << &seedPoint[1] <<
 283 //     "seedp2.="  << &seedPoint[2] <<
 284     "fitT0="    << fitT0 <<
 285     "\n";
 286
 287   tVtx=track->GetZ();
 288   x=track->GetX();
 289   delete track;
 290   return 0;
 291 }
 292
 293 //____________________________________________________________________________
 294 void SetTrackPointFromCluster(AliTPCclusterMI *cl, AliTrackPoint &p ) {
 295   //
 296   // make AliTrackPoint out of AliTPCclusterMI
 297   //
 298
 299   if (!cl) return;
 300 //   Float_t xyz[3]={0.,0.,0.};
 301 //   ClusterToSpacePoint(cl,xyz);
 302 //   cl->GetGlobalCov(cov);
 303   //TODO: what to do with the covariance matrix???
 304   //TODO: the problem is that it is used in GetAngle in AliTrackPoint
 305   //TODO: which is used by AliTrackerBase::MakeSeed to get alpha correct ...
 306   //TODO: for the moment simply assign 1 permill squared
 307   // in AliTrackPoint the cov is xx, xy, xz, yy, yz, zz
 308 //   Float_t cov[6]={xyz[0]*xyz[0]*1e-6,xyz[0]*xyz[1]*1e-6,xyz[0]*xyz[2]*1e-6,
 309 //                   xyz[1]*xyz[1]*1e-6,xyz[1]*xyz[2]*1e-6,xyz[2]*xyz[2]*1e-6};
 310 //   cl->GetGlobalXYZ(xyz);
 311 //   cl->GetGlobalCov(cov);
 312   // voluem ID to add later ....
 313 //   p.SetXYZ(xyz);
 314 //   p.SetCov(cov);
 315   AliTrackPoint *tp=cl->MakePoint();
 316   p=*tp;
 317   delete tp;
 318 //   cl->Print();
 319 //   p.Print();
 320   p.SetVolumeID(cl->GetDetector());
 321 //   p.Rotate(p.GetAngle()).Print();
 322 }
 323
 324 //____________________________________________________________________________
 325 void ClusterToSpacePoint(AliTPCclusterMI *cl, Float_t xyz[3])
 326 {
 327   //
 328   // convert the cluster to a space point in global coordinates
 329   //
 330   if (!cl) return;
 331   xyz[0]=cl->GetRow();
 332   xyz[1]=cl->GetPad();
 333   xyz[2]=cl->GetTimeBin(); // this will not be correct at all
 334   Int_t i[3]={0,cl->GetDetector(),cl->GetRow()};
 335 //   printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
 336   fTPCParam->Transform8to4(xyz,i);
 337 //   printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
 338   fTPCParam->Transform4to3(xyz,i);
 339 //   printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
 340   fTPCParam->Transform2to1(xyz,i);
 341 //   printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
 342 }
 343
 344 //____________________________________________________________________________
 345 void InitSpaceCharge(TTree *t)
 346 {
 347   //
 348   // Init the space charge map
 349   //
 350
 351   TString filename="$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps5_50kHz_precal.root";
 352   if (t) {
 353     TList *l=t->GetUserInfo();
 354     for (Int_t i=0; i<l->GetEntries(); ++i) {
 355       TString s(l->At(i)->GetName());
 356       if (s.Contains("SC_")) {
 357         filename=s;
 358         break;
 359       }
 360     }
 361   }
 362
 363   printf("Initialising the space charge map using the file: '%s'\n",filename.Data());
 364   TFile f(filename.Data());
 365   fSpaceCharge=(AliTPCSpaceCharge3D*)f.Get("map");
 366
 367 //   fSpaceCharge = new AliTPCSpaceCharge3D();
 368 //   fSpaceCharge->SetSCDataFileName("$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps10_50kHz.root");
 369 //   fSpaceCharge->SetOmegaTauT1T2(0.325,1,1); // Ne CO2
 370 // //   fSpaceCharge->SetOmegaTauT1T2(0.41,1,1.05); // Ar CO2
 371 //   fSpaceCharge->InitSpaceCharge3DDistortion();
 372
 373 }
 374
 375 //____________________________________________________________________________
 376
 377
 378 AliExternalTrackParam* GetFullTrack(AliToyMCTrack *tr, Int_t clsType=0, Int_t corrType=0, Bool_t useMaterial=kFALSE)
 379 {
 380   //
 381   // clsType:  0=undistorted clusters; 1: distorted clusters
 382   // corrType: 0=none; 1: ideal
 383   //
 384
 385   // no correction for undistorted clusters
 386   if (clsType==0) corrType=0;
 387
 388   AliTPCROC * roc = AliTPCROC::Instance();
 389   const Int_t    npoints0=roc->GetNRows(0)+roc->GetNRows(36);
 390   const Double_t kRTPC0  =roc->GetPadRowRadii(0,0);
 391   const Double_t kRTPC1  =roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
 392   const Double_t kMaxSnp = 0.85;
 393   const Double_t kSigmaY=0.1;
 394   const Double_t kSigmaZ=0.1;
 395   const Double_t kMaxR=500;
 396   const Double_t kMaxZ=500;
 397
 398   const Double_t kMaxZ0=220;
 399   const Double_t kZcut=3;
 400   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 401
 402   Int_t ncls=(clsType==0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
 403
 404   Int_t dir = -1;
 405   Double_t refX = 0.;
 406
 407   // get points for the seed
 408   Int_t seed=0;
 409   AliTrackPoint    seedPoint[3];
 410   for (Int_t ipoint=ncls-1; ipoint>=ncls-3; --ipoint){
 411     AliTPCclusterMI *cl=tr->GetSpacePoint(ipoint);
 412     if (clsType==1) cl=tr->GetDistortedSpacePoint(ipoint);
 413
 414     SetTrackPointFromCluster(cl, seedPoint[seed]);
 415
 416     if (corrType==1){
 417       Float_t xyz[3]={0,0,0};
 418       seedPoint[seed].GetXYZ(xyz);
 419       fSpaceCharge->CorrectPoint(xyz, cl->GetDetector());
 420       seedPoint[seed].SetXYZ(xyz);
 421     }
 422     seedPoint[seed].Print();
 423     ++seed;
 424   }
 425
 426   AliExternalTrackParam *track = 0x0;
 427   track = AliTrackerBase::MakeSeed(seedPoint[2], seedPoint[1], seedPoint[0]);
 428   track->ResetCovariance(10);
 429
 430 //   printf("============================================\n");
 431 //   printf("orig:  %.2f, %.2f, %.2f, %.2f, %.2f, %.2f\n",
 432 //          tr->GetParameter()[0],tr->GetParameter()[1],tr->GetParameter()[2],
 433 //          tr->GetParameter()[3],tr->GetParameter()[4],tr->GetAlpha());
 434
 435 //   printf("seed:  %.2f, %.2f, %.2f, %.2f, %.2f, %.2f\n",
 436 //          track->GetParameter()[0],track->GetParameter()[1],track->GetParameter()[2],
 437 //          track->GetParameter()[3],track->GetParameter()[4],track->GetAlpha());
 438
 439   // loop over all other points and add to the track
 440   for (Int_t ipoint=ncls-4; ipoint>=0; --ipoint){
 441     AliTrackPoint pIn;
 442     AliTPCclusterMI *cl=tr->GetSpacePoint(ipoint);
 443     if (clsType==1) cl=tr->GetDistortedSpacePoint(ipoint);
 444     SetTrackPointFromCluster(cl, pIn);
 445     if (corrType==1){
 446       Float_t xyz[3]={0,0,0};
 447       pIn.GetXYZ(xyz);
 448       fSpaceCharge->CorrectPoint(xyz, cl->GetDetector());
 449       pIn.SetXYZ(xyz);
 450     }
 451     // rotate the cluster to the local detector frame
 452     track->Rotate(((cl->GetDetector()%18)*20+10)*TMath::DegToRad());
 453     AliTrackPoint prot = pIn.Rotate(track->GetAlpha());   // rotate to the local frame - non distoted  point
 454     if (TMath::Abs(prot.GetX())<kRTPC0) continue;
 455     if (TMath::Abs(prot.GetX())>kRTPC1) continue;
 456     //
 457     Int_t ret=0;
 458 //     printf("before: %.2f, %.2f, %.2f, %.2f, %.2f, %.2f\n",
 459 //            track->GetParameter()[0],track->GetParameter()[1],track->GetParameter()[2],
 460 //            track->GetParameter()[3],track->GetParameter()[4],track->GetAlpha());
 461     if (useMaterial) ret=AliTrackerBase::PropagateTrackTo2(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp);
 462     else ret=AliTrackerBase::PropagateTrackTo2(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp,0,kFALSE,kFALSE);
 463 //     printf("after: %.2f, %.2f, %.2f, %.2f, %.2f, %.2f\n",
 464 //            track->GetParameter()[0],track->GetParameter()[1],track->GetParameter()[2],
 465 //            track->GetParameter()[3],track->GetParameter()[4],track->GetAlpha());
 466     if (ret<0) {
 467       (*fStreamer) << "np" <<
 468       "iev="      << fEvent <<
 469       "itr="      << fTrack <<
 470       "track.="   << tr     <<
 471       "seed.="    << track  <<
 472       "clsType="  << clsType <<
 473       "corrType=" << corrType <<
 474       "\n";
 475       printf("Could not propagate track: %d\n",ret);
 476       break;
 477     }
 478 //     printf("\n=========\n%d:\n",ipoint);
 479 //     printf("%.2f, %.2f, %.2f - %d, %d, %.2f, %.2g\n",cl->GetX(),cl->GetY(),cl->GetZ(),cl->GetDetector(),cl->GetRow(),cl->GetPad(),cl->GetTimeBin());
 480 //     printf("%.2f, %.2f, %.2f - %.2f\n", prot.GetX(),prot.GetY(),prot.GetZ(), prot.GetAngle());
 481 //     printf("%.2f, %.2f, %.2f - %.2f\n", track->GetX(),track->GetY(),track->GetZ(), track->GetAlpha());
 482
 483     if (TMath::Abs(track->GetZ())>kMaxZ) break;
 484     if (TMath::Abs(track->GetX())>kMaxR) break;
 485     if (dir>0 && track->GetX()>refX) continue;
 486     if (dir<0 && track->GetX()<refX) continue;
 487     if (TMath::Abs(track->GetZ())<kZcut)continue;
 488     //
 489     Double_t pointPos[2]={0,0};
 490     Double_t pointCov[3]={0,0,0};
 491     pointPos[0]=prot.GetY();//local y
 492     pointPos[1]=prot.GetZ();//local z
 493     pointCov[0]=prot.GetCov()[3];//simay^2
 494     pointCov[1]=prot.GetCov()[4];//sigmayz
 495     pointCov[2]=prot.GetCov()[5];//sigmaz^2
 496     if (!track->Update(pointPos,pointCov)) {printf("no update\n"); break;}
 497   }
 498 //   printf(">>> before2: %.2f, %.2f, %.2f, %.2f, %.2f, %.2f\n",
 499 //          track->GetParameter()[0],track->GetParameter()[1],track->GetParameter()[2],
 500 //          track->GetParameter()[3],track->GetParameter()[4],track->GetAlpha());
 501   if (useMaterial) AliTrackerBase::PropagateTrackTo2(track,refX,kMass,5.,kTRUE,kMaxSnp);
 502   else AliTrackerBase::PropagateTrackTo2(track,refX,kMass,5.,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
 503 //   printf(">>> after2: %.2f, %.2f, %.2f, %.2f, %.2f, %.2f\n",
 504 //          track->GetParameter()[0],track->GetParameter()[1],track->GetParameter()[2],
 505 //          track->GetParameter()[3],track->GetParameter()[4],track->GetAlpha());
 506   if (useMaterial) ret=AliTrackerBase::PropagateTrackTo2(track,refX,kMass,1.,kTRUE,kMaxSnp);
 507   else ret=AliTrackerBase::PropagateTrackTo2(track,refX,kMass,1.,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
 508 //   printf(">>> after2.2: %.2f, %.2f, %.2f, %.2f, %.2f, %.2f\n",
 509 //          track->GetParameter()[0],track->GetParameter()[1],track->GetParameter()[2],
 510 //          track->GetParameter()[3],track->GetParameter()[4],track->GetAlpha());
 511   printf("Propagation to 0 stopped at %.2f with %d\n",track->GetX(),ret);
 512   track->Rotate(tr->GetAlpha());
 513   return track;
 514 }
 515
 516 //____________________________________________________________________________
 517 void testResolution(const char* filename, Int_t nmaxEv=-1, Bool_t useMaterial=kFALSE)
 518 {
 519   TFile f(filename);
 520   if (!f.IsOpen() || f.IsZombie()) {
 521     printf("ERROR: couldn't open the file '%s'\n", filename);
 522     return;
 523   }
 524
 525   TTree *t=(TTree*)f.Get("toyMCtree");
 526   if (!t) {
 527     printf("ERROR: couldn't read the 'toyMCtree' from file '%s'\n", filename);
 528     return;
 529   }
 530
 531   AliToyMCEvent *ev=0x0;
 532   t->SetBranchAddress("event",&ev);
 533
 534   TString debugName=filename;
 535   debugName.ReplaceAll(".root","");
 536   if (useMaterial) debugName.Append(".Mat");
 537   else debugName.Append(".noMat");
 538   debugName.Append(".testRes.root");
 539
 540   gSystem->Exec(Form("test -f %s && rm %s", debugName.Data(), debugName.Data()));
 541   if (!fStreamer) fStreamer=new TTreeSRedirector(debugName.Data());
 542
 543   gROOT->cd();
 544
 545   // read spacecharge from the Userinfo ot the tree
 546   InitSpaceCharge(t);
 547
 548   Int_t maxev=t->GetEntries();
 549   if (nmaxEv>0&&nmaxEv<maxev) maxev=nmaxEv;
 550
 551   for (Int_t iev=0; iev<maxev; ++iev){
 552     t->GetEvent(iev);
 553     fEvent=iev;
 554     printf("========== Processing event %3d =============\n",iev);
 555     for (Int_t itr=0; itr<ev->GetNumberOfTracks(); ++itr){
 556       fTrack=itr;
 557       printf("   ======= Processing track %3d ==========\n",itr);
 558       AliToyMCTrack *tr=ev->GetTrack(itr);
 559       AliExternalTrackParam tOrig(*tr);
 560       AliExternalTrackParam *tIdeal    = GetFullTrack(tr,0,0,useMaterial);
 561       AliExternalTrackParam *tDist     = GetFullTrack(tr,1);
 562       AliExternalTrackParam *tDistCorr = GetFullTrack(tr,1,1);
 563
 564       (*fStreamer) << "res" <<
 565       "iev="        << iev <<
 566       "itr="        << itr <<
 567       "tOrig.="     << &tOrig <<
 568       "tIdeal.="    << tIdeal <<
 569       "tDist.="     << tDist  <<
 570       "tDistCorr.=" << tDistCorr <<
 571       "\n";
 572
 573       delete tIdeal;
 574       delete tDist;
 575       delete tDistCorr;
 576     }
 577   }
 578
 579   delete fStreamer;
 580   fStreamer=0x0;
 581 }
 582