TPC/Upgrade/macros/testRec.C

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