TPC/Upgrade/AliToyMCReconstruction.cxx

   1
   2 #include <TDatabasePDG.h>
   3 #include <TString.h>
   4 #include <TSystem.h>
   5 #include <TROOT.h>
   6 #include <TFile.h>
   7 #include <TPRegexp.h>
   8
   9 #include <AliExternalTrackParam.h>
  10 #include <AliTPCcalibDB.h>
  11 #include <AliTPCclusterMI.h>
  12 #include <AliTPCSpaceCharge3D.h>
  13 #include <AliTrackerBase.h>
  14 #include <AliTrackPointArray.h>
  15 #include <AliLog.h>
  16 #include <AliTPCParam.h>
  17 #include <AliTPCROC.h>
  18 #include <TTreeStream.h>
  19 #include <AliTPCtracker.h>
  20 #include <AliTPCtrackerSector.h>
  21
  22 #include "AliToyMCTrack.h"
  23 #include "AliToyMCEvent.h"
  24
  25 #include "AliToyMCReconstruction.h"
  26
  27
  28 //____________________________________________________________________________________
  29 AliToyMCReconstruction::AliToyMCReconstruction() : TObject()
  30 , fSeedingRow(140)
  31 , fSeedingDist(10)
  32 , fClusterType(0)
  33 , fCorrectionType(kNoCorrection)
  34 , fDoTrackFit(kTRUE)
  35 , fUseMaterial(kFALSE)
  36 , fIdealTracking(kFALSE)
  37 , fTime0(-1)
  38 , fCreateT0seed(kFALSE)
  39 , fStreamer(0x0)
  40 , fTree(0x0)
  41 , fEvent(0x0)
  42 , fTPCParam(0x0)
  43 , fSpaceCharge(0x0)
  44 {
  45   //
  46   //  ctor
  47   //
  48   fTPCParam=AliTPCcalibDB::Instance()->GetParameters();
  49 }
  50
  51 //____________________________________________________________________________________
  52 AliToyMCReconstruction::~AliToyMCReconstruction()
  53 {
  54   //
  55   //  dtor
  56   //
  57
  58   delete fStreamer;
  59 //   delete fTree;
  60 }
  61
  62 //____________________________________________________________________________________
  63 void AliToyMCReconstruction::RunReco(const char* file, Int_t nmaxEv)
  64 {
  65   //
  66   // Recostruction from associated clusters
  67   //
  68
  69   TFile f(file);
  70   if (!f.IsOpen() || f.IsZombie()) {
  71     printf("ERROR: couldn't open the file '%s'\n", file);
  72     return;
  73   }
  74
  75  fTree=(TTree*)f.Get("toyMCtree");
  76   if (!fTree) {
  77     printf("ERROR: couldn't read the 'toyMCtree' from file '%s'\n", file);
  78     return;
  79   }
  80
  81   fEvent=0x0;
  82   fTree->SetBranchAddress("event",&fEvent);
  83
  84   // read spacecharge from the Userinfo ot the tree
  85   InitSpaceCharge();
  86
  87   TString debugName=file;
  88   debugName.ReplaceAll(".root","");
  89   debugName.Append(Form(".%1d.%1d_%1d_%1d_%03d_%02d",
  90                         fUseMaterial,fIdealTracking,fClusterType,
  91                         Int_t(fCorrectionType),fSeedingRow,fSeedingDist));
  92   debugName.Append(".debug.root");
  93
  94   gSystem->Exec(Form("test -f %s && rm %s", debugName.Data(), debugName.Data()));
  95   if (!fStreamer) fStreamer=new TTreeSRedirector(debugName.Data());
  96
  97   gROOT->cd();
  98
  99   static AliExternalTrackParam dummySeedT0;
 100   static AliExternalTrackParam dummySeed;
 101   static AliExternalTrackParam dummyTrack;
 102
 103   AliExternalTrackParam t0seed;
 104   AliExternalTrackParam seed;
 105   AliExternalTrackParam track;
 106   AliExternalTrackParam tOrig;
 107
 108   AliExternalTrackParam *dummy;
 109
 110   Int_t maxev=fTree->GetEntries();
 111   if (nmaxEv>0&&nmaxEv<maxev) maxev=nmaxEv;
 112
 113   for (Int_t iev=0; iev<maxev; ++iev){
 114     printf("==============  Processing Event %6d =================\n",iev);
 115     fTree->GetEvent(iev);
 116     for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
 117       printf(" > ======  Processing Track %6d ========  \n",itr);
 118       const AliToyMCTrack *tr=fEvent->GetTrack(itr);
 119       tOrig = *tr;
 120
 121
 122       // set dummy
 123       t0seed    = dummySeedT0;
 124       seed      = dummySeed;
 125       track     = dummyTrack;
 126
 127       Float_t z0=fEvent->GetZ();
 128       Float_t t0=fEvent->GetT0();
 129
 130       Float_t vdrift=GetVDrift();
 131       Float_t zLength=GetZLength(0);
 132
 133       // crate time0 seed, steered by fCreateT0seed
 134       printf("t0 seed\n");
 135       fTime0=-1.;
 136       fCreateT0seed=kTRUE;
 137       dummy = GetSeedFromTrack(tr);
 138
 139       if (dummy) {
 140         t0seed = *dummy;
 141         delete dummy;
 142
 143         // crate real seed using the time 0 from the first seed
 144         // set fCreateT0seed now to false to get the seed in z coordinates
 145         fTime0 = t0seed.GetZ()-zLength/vdrift;
 146         fCreateT0seed = kFALSE;
 147         printf("seed (%.2g)\n",fTime0);
 148         dummy  = GetSeedFromTrack(tr);
 149         if (dummy) {
 150           seed = *dummy;
 151           delete dummy;
 152
 153           // create fitted track
 154           if (fDoTrackFit){
 155             printf("track\n");
 156             dummy = GetFittedTrackFromSeed(tr, &seed);
 157             track = *dummy;
 158             delete dummy;
 159           }
 160
 161           // propagate seed to 0
 162           const Double_t kMaxSnp = 0.85;
 163           const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 164 //           AliTrackerBase::PropagateTrackTo(&seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
 165
 166         }
 167       }
 168
 169       Int_t ctype(fCorrectionType);
 170
 171       if (fStreamer) {
 172         (*fStreamer) << "Tracks" <<
 173         "iev="         << iev             <<
 174         "z0="          << z0              <<
 175         "t0="          << t0              <<
 176         "fTime0="      << fTime0          <<
 177         "itr="         << itr             <<
 178         "clsType="     << fClusterType    <<
 179         "corrType="    << ctype           <<
 180         "seedRow="     << fSeedingRow     <<
 181         "seedDist="    << fSeedingDist    <<
 182         "vdrift="      << vdrift          <<
 183         "zLength="     << zLength         <<
 184         "t0seed.="     << &t0seed         <<
 185         "seed.="       << &seed           <<
 186         "track.="      << &track          <<
 187         "tOrig.="      << &tOrig          <<
 188         "\n";
 189       }
 190
 191
 192     }
 193   }
 194
 195   delete fStreamer;
 196   fStreamer=0x0;
 197
 198   delete fEvent;
 199   fEvent = 0x0;
 200
 201   delete fTree;
 202   fTree=0x0;
 203   f.Close();
 204 }
 205
 206
 207 //____________________________________________________________________________________
 208 void AliToyMCReconstruction::RunRecoAllClusters(const char* file, Int_t nmaxEv)
 209 {
 210   //
 211   // Reconstruction for seed from associated clusters, but array of clusters
 212   //
 213
 214   TFile f(file);
 215   if (!f.IsOpen() || f.IsZombie()) {
 216     printf("ERROR: couldn't open the file '%s'\n", file);
 217     return;
 218   }
 219
 220  fTree=(TTree*)f.Get("toyMCtree");
 221   if (!fTree) {
 222     printf("ERROR: couldn't read the 'toyMCtree' from file '%s'\n", file);
 223     return;
 224   }
 225
 226   fEvent=0x0;
 227   fTree->SetBranchAddress("event",&fEvent);
 228
 229   // read spacecharge from the Userinfo ot the tree
 230   InitSpaceCharge();
 231
 232   TString debugName=file;
 233   debugName.ReplaceAll(".root","");
 234   debugName.Append(Form(".%1d.%1d_%1d_%1d_%03d_%02d",
 235                         fUseMaterial,fIdealTracking,fClusterType,
 236                         Int_t(fCorrectionType),fSeedingRow,fSeedingDist));
 237   debugName.Append(".allClusters.debug.root");
 238
 239   gSystem->Exec(Form("test -f %s && rm %s", debugName.Data(), debugName.Data()));
 240   if (!fStreamer) fStreamer=new TTreeSRedirector(debugName.Data());
 241
 242   gROOT->cd();
 243
 244   static AliExternalTrackParam dummySeedT0;
 245   static AliExternalTrackParam dummySeed;
 246   static AliExternalTrackParam dummyTrack;
 247
 248   AliExternalTrackParam t0seed;
 249   AliExternalTrackParam seed;
 250   AliExternalTrackParam track;
 251   AliExternalTrackParam tOrig;
 252
 253   // cluster array of all sectors
 254   const Int_t fkNSectorInner = fTPCParam->GetNInnerSector()/2;
 255   const Int_t fkNSectorOuter = fTPCParam->GetNOuterSector()/2;
 256
 257   AliTPCtrackerSector *fInnerSectorArray = new AliTPCtrackerSector[fkNSectorInner];
 258   AliTPCtrackerSector *fOuterSectorArray = new AliTPCtrackerSector[fkNSectorOuter];
 259
 260   for (Int_t i=0; i<fkNSectorInner; ++i) fInnerSectorArray[i].Setup(fTPCParam,0);
 261   for (Int_t i=0; i<fkNSectorOuter; ++i) fOuterSectorArray[i].Setup(fTPCParam,1);
 262
 263   Int_t count[72][96] = { {0} , {0} };
 264
 265
 266
 267   AliExternalTrackParam *dummy;
 268
 269   Int_t maxev=fTree->GetEntries();
 270   if (nmaxEv>0&&nmaxEv<maxev) maxev=nmaxEv;
 271
 272
 273   // ===========================================================================================
 274   // Loop 1: Fill AliTPCtrackerSector structure
 275   // ===========================================================================================
 276   for (Int_t iev=0; iev<maxev; ++iev){
 277     printf("==============  Fill Clusters: Processing Event %6d  =================\n",iev);
 278     fTree->GetEvent(iev);
 279     for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
 280       printf(" > ======  Fill Clusters: Processing Track %6d ========  \n",itr);
 281       const AliToyMCTrack *tr=fEvent->GetTrack(itr);
 282
 283       // number of clusters to loop over
 284       const Int_t ncls=(fClusterType==0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
 285
 286       for(Int_t icl=0; icl<ncls; ++icl){
 287
 288         const AliTPCclusterMI *cl=tr->GetSpacePoint(icl);
 289         if (fClusterType==1) cl=tr->GetDistortedSpacePoint(icl);
 290         if (!cl) continue;
 291
 292         Int_t sec = cl->GetDetector();
 293         Int_t row = cl->GetRow();
 294
 295         // fill arrays for inner and outer sectors (A/C side handled internally)
 296         if (sec<fkNSectorInner*2){
 297           fInnerSectorArray[sec%fkNSectorInner].InsertCluster(const_cast<AliTPCclusterMI*>(cl), count[sec][row], fTPCParam);
 298         }
 299         else{
 300           fOuterSectorArray[(sec-fkNSectorInner*2)%fkNSectorOuter].InsertCluster(const_cast<AliTPCclusterMI*>(cl), count[sec][row], fTPCParam);
 301         }
 302
 303         ++count[sec][row];
 304       }
 305     }
 306   }
 307
 308
 309   // ===========================================================================================
 310   // Loop 2a: Use the full TPC tracker
 311   // TODO: how to use???
 312   // ===========================================================================================
 313   AliTPCtracker *tpcTracker = new AliTPCtracker;
 314
 315   // ===========================================================================================
 316   // Loop 2b: Seeding from clusters associated to tracks
 317   // TODO: Implement tracking from given seed!
 318   // ===========================================================================================
 319   for (Int_t iev=0; iev<maxev; ++iev){
 320     printf("==============  Processing Event %6d =================\n",iev);
 321     fTree->GetEvent(iev);
 322     for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
 323       printf(" > ======  Processing Track %6d ========  \n",itr);
 324       const AliToyMCTrack *tr=fEvent->GetTrack(itr);
 325       tOrig = *tr;
 326
 327
 328       // set dummy
 329       t0seed    = dummySeedT0;
 330       seed      = dummySeed;
 331       track     = dummyTrack;
 332
 333       Float_t z0=fEvent->GetZ();
 334       Float_t t0=fEvent->GetT0();
 335
 336       Float_t vdrift=GetVDrift();
 337       Float_t zLength=GetZLength(0);
 338
 339       // crate time0 seed, steered by fCreateT0seed
 340       printf("t0 seed\n");
 341       fTime0=-1.;
 342       fCreateT0seed=kTRUE;
 343       dummy = GetSeedFromTrack(tr);
 344
 345       if (dummy) {
 346         t0seed = *dummy;
 347         delete dummy;
 348
 349         // crate real seed using the time 0 from the first seed
 350         // set fCreateT0seed now to false to get the seed in z coordinates
 351         fTime0 = t0seed.GetZ()-zLength/vdrift;
 352         fCreateT0seed = kFALSE;
 353         printf("seed (%.2g)\n",fTime0);
 354         dummy  = GetSeedFromTrack(tr);
 355         if (dummy) {
 356           seed = *dummy;
 357           delete dummy;
 358
 359           // create fitted track
 360           if (fDoTrackFit){
 361             printf("track\n");
 362             dummy = GetFittedTrackFromSeedAllClusters(tr, &seed, fInnerSectorArray, fOuterSectorArray);
 363             track = *dummy;
 364             delete dummy;
 365           }
 366
 367           // propagate seed to 0
 368           const Double_t kMaxSnp = 0.85;
 369           const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 370           AliTrackerBase::PropagateTrackTo(&seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
 371
 372         }
 373       }
 374
 375       Int_t ctype(fCorrectionType);
 376
 377       if (fStreamer) {
 378         (*fStreamer) << "Tracks" <<
 379         "iev="         << iev             <<
 380         "z0="          << z0              <<
 381         "t0="          << t0              <<
 382         "fTime0="      << fTime0          <<
 383         "itr="         << itr             <<
 384         "clsType="     << fClusterType    <<
 385         "corrType="    << ctype           <<
 386         "seedRow="     << fSeedingRow     <<
 387         "seedDist="    << fSeedingDist    <<
 388         "vdrift="      << vdrift          <<
 389         "zLength="     << zLength         <<
 390         "t0seed.="     << &t0seed         <<
 391         "seed.="       << &seed           <<
 392         "track.="      << &track          <<
 393         "tOrig.="      << &tOrig          <<
 394         "\n";
 395       }
 396
 397
 398     }
 399   }
 400
 401
 402   delete fStreamer;
 403   fStreamer=0x0;
 404
 405   delete fEvent;
 406   fEvent = 0x0;
 407
 408   delete fTree;
 409   fTree=0x0;
 410   f.Close();
 411 }
 412
 413
 414 //____________________________________________________________________________________
 415 AliExternalTrackParam* AliToyMCReconstruction::GetSeedFromTrack(const AliToyMCTrack * const tr)
 416 {
 417   //
 418   // if we don't have a valid time0 informaion (fTime0) available yet
 419   // assume we create a seed for the time0 estimate
 420   //
 421
 422   // seed point informaion
 423   AliTrackPoint    seedPoint[3];
 424   const AliTPCclusterMI *seedCluster[3]={0x0,0x0,0x0};
 425
 426   // number of clusters to loop over
 427   const Int_t ncls=(fClusterType==0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
 428
 429   UChar_t nextSeedRow=fSeedingRow;
 430   Int_t   nseeds=0;
 431
 432   //assumes sorted clusters
 433   for (Int_t icl=0;icl<ncls;++icl) {
 434     const AliTPCclusterMI *cl=tr->GetSpacePoint(icl);
 435     if (fClusterType==1) cl=tr->GetDistortedSpacePoint(icl);
 436     if (!cl) continue;
 437     // use row in sector
 438     const UChar_t row=cl->GetRow() + 63*(cl->GetDetector()>35);
 439     // skip clusters without proper pad row
 440     if (row>200) continue;
 441
 442     //check seeding row
 443     // if we are in the last row and still miss a seed we use the last row
 444     //   even if the row spacing will not be equal
 445     if (row>=nextSeedRow || icl==ncls-1){
 446       seedCluster[nseeds]=cl;
 447       SetTrackPointFromCluster(cl, seedPoint[nseeds]);
 448       ++nseeds;
 449       nextSeedRow+=fSeedingDist;
 450
 451       if (nseeds==3) break;
 452     }
 453   }
 454
 455   // check we really have 3 seeds
 456   if (nseeds!=3) {
 457     AliError(Form("Seeding failed for parameters %d, %d\n",fSeedingDist,fSeedingRow));
 458     return 0x0;
 459   }
 460
 461   // do cluster correction for fCorrectionType:
 462   //   0 - no correction
 463   //   1 - TPC center
 464   //   2 - average eta
 465   //   3 - ideal
 466   // assign the cluster abs time as z component to all seeds
 467   for (Int_t iseed=0; iseed<3; ++iseed) {
 468     Float_t xyz[3]={0,0,0};
 469     seedPoint[iseed].GetXYZ(xyz);
 470     const Float_t r=TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
 471
 472     const Int_t sector=seedCluster[iseed]->GetDetector();
 473     const Int_t sign=1-2*((sector/18)%2);
 474
 475     if ( (fClusterType == 1) && (fCorrectionType != kNoCorrection) ) {
 476       printf("correction type: %d\n",(Int_t)fCorrectionType);
 477
 478       // the settings below are for the T0 seed
 479       // for known T0 the z position is already calculated in SetTrackPointFromCluster
 480       if ( fCreateT0seed ){
 481         if ( fCorrectionType == kTPCCenter  ) xyz[2] = 125.*sign;
 482         //!!! TODO: is this the correct association?
 483         if ( fCorrectionType == kAverageEta ) xyz[2] = TMath::Tan(45./2.*TMath::DegToRad())*r*sign;
 484       }
 485
 486       if ( fCorrectionType == kIdeal      ) xyz[2] = seedCluster[iseed]->GetZ();
 487
 488       //!!! TODO: to be replaced with the proper correction
 489       fSpaceCharge->CorrectPoint(xyz, seedCluster[iseed]->GetDetector());
 490     }
 491
 492     // after the correction set the time bin as z-Position in case of a T0 seed
 493     if ( fCreateT0seed )
 494       xyz[2]=seedCluster[iseed]->GetTimeBin();
 495
 496     seedPoint[iseed].SetXYZ(xyz);
 497   }
 498
 499   const Double_t kMaxSnp = 0.85;
 500   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 501
 502   AliExternalTrackParam *seed = AliTrackerBase::MakeSeed(seedPoint[0], seedPoint[1], seedPoint[2]);
 503   seed->ResetCovariance(10);
 504
 505   if (fCreateT0seed){
 506     // if fTime0 < 0 we assume that we create a seed for the T0 estimate
 507     AliTrackerBase::PropagateTrackTo(seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
 508     if (TMath::Abs(seed->GetX())>3) {
 509       printf("Could not propagate track to 0, %.2f, %.2f, %.2f\n",seed->GetX(),seed->GetAlpha(),seed->Pt());
 510     }
 511   }
 512
 513   return seed;
 514
 515 }
 516
 517 //____________________________________________________________________________________
 518 void AliToyMCReconstruction::SetTrackPointFromCluster(const AliTPCclusterMI *cl, AliTrackPoint &p )
 519 {
 520   //
 521   // make AliTrackPoint out of AliTPCclusterMI
 522   //
 523
 524   if (!cl) return;
 525     Float_t xyz[3]={0.,0.,0.};
 526   //   ClusterToSpacePoint(cl,xyz);
 527   //   cl->GetGlobalCov(cov);
 528   //TODO: what to do with the covariance matrix???
 529   //TODO: the problem is that it is used in GetAngle in AliTrackPoint
 530   //TODO: which is used by AliTrackerBase::MakeSeed to get alpha correct ...
 531   //TODO: for the moment simply assign 1 permill squared
 532   // in AliTrackPoint the cov is xx, xy, xz, yy, yz, zz
 533   //   Float_t cov[6]={xyz[0]*xyz[0]*1e-6,xyz[0]*xyz[1]*1e-6,xyz[0]*xyz[2]*1e-6,
 534   //                   xyz[1]*xyz[1]*1e-6,xyz[1]*xyz[2]*1e-6,xyz[2]*xyz[2]*1e-6};
 535   //   cl->GetGlobalXYZ(xyz);
 536   //   cl->GetGlobalCov(cov);
 537   // voluem ID to add later ....
 538   //   p.SetXYZ(xyz);
 539   //   p.SetCov(cov);
 540   AliTrackPoint *tp=const_cast<AliTPCclusterMI*>(cl)->MakePoint();
 541   p=*tp;
 542   delete tp;
 543   //   cl->Print();
 544   //   p.Print();
 545   p.SetVolumeID(cl->GetDetector());
 546
 547
 548   if ( !fCreateT0seed && !fIdealTracking ) {
 549     p.GetXYZ(xyz);
 550     const Int_t sector=cl->GetDetector();
 551     const Int_t sign=1-2*((sector/18)%2);
 552     const Float_t zT0=( GetZLength(sector) - (cl->GetTimeBin()-fTime0)*GetVDrift() )*sign;
 553     printf(" z:  %.2f  %.2f\n",xyz[2],zT0);
 554     xyz[2]=zT0;
 555     p.SetXYZ(xyz);
 556   }
 557
 558
 559   //   p.Rotate(p.GetAngle()).Print();
 560 }
 561
 562 //____________________________________________________________________________________
 563 void AliToyMCReconstruction::ClusterToSpacePoint(const AliTPCclusterMI *cl, Float_t xyz[3])
 564 {
 565   //
 566   // convert the cluster to a space point in global coordinates
 567   //
 568   if (!cl) return;
 569   xyz[0]=cl->GetRow();
 570   xyz[1]=cl->GetPad();
 571   xyz[2]=cl->GetTimeBin(); // this will not be correct at all
 572   Int_t i[3]={0,cl->GetDetector(),cl->GetRow()};
 573   //   printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
 574   fTPCParam->Transform8to4(xyz,i);
 575   //   printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
 576   fTPCParam->Transform4to3(xyz,i);
 577   //   printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
 578   fTPCParam->Transform2to1(xyz,i);
 579   //   printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
 580 }
 581
 582 //____________________________________________________________________________________
 583 AliExternalTrackParam* AliToyMCReconstruction::GetFittedTrackFromSeed(const AliToyMCTrack *tr, const AliExternalTrackParam *seed)
 584 {
 585   //
 586   //
 587   //
 588
 589   // create track
 590   AliExternalTrackParam *track = new AliExternalTrackParam(*seed);
 591
 592   Int_t ncls=(fClusterType == 0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
 593
 594   const AliTPCROC * roc = AliTPCROC::Instance();
 595
 596   const Double_t kRTPC0  = roc->GetPadRowRadii(0,0);
 597   const Double_t kRTPC1  = roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
 598   const Double_t kMaxSnp = 0.85;
 599   const Double_t kMaxR   = 500.;
 600   const Double_t kMaxZ   = 500.;
 601
 602   //   const Double_t kMaxZ0=220;
 603 //   const Double_t kZcut=3;
 604
 605   const Double_t refX = tr->GetX();
 606
 607   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 608
 609   // loop over all other points and add to the track
 610   for (Int_t ipoint=ncls-1; ipoint>=0; --ipoint){
 611     AliTrackPoint pIn;
 612     const AliTPCclusterMI *cl=tr->GetSpacePoint(ipoint);
 613     if (fClusterType == 1) cl=tr->GetDistortedSpacePoint(ipoint);
 614     SetTrackPointFromCluster(cl, pIn);
 615     if (fCorrectionType != kNoCorrection){
 616       Float_t xyz[3]={0,0,0};
 617       pIn.GetXYZ(xyz);
 618 //       if ( fCorrectionType == kIdeal ) xyz[2] = cl->GetZ();
 619       fSpaceCharge->CorrectPoint(xyz, cl->GetDetector());
 620       pIn.SetXYZ(xyz);
 621     }
 622     // rotate the cluster to the local detector frame
 623     track->Rotate(((cl->GetDetector()%18)*20+10)*TMath::DegToRad());
 624     AliTrackPoint prot = pIn.Rotate(track->GetAlpha());   // rotate to the local frame - non distoted  point
 625     if (TMath::Abs(prot.GetX())<kRTPC0) continue;
 626     if (TMath::Abs(prot.GetX())>kRTPC1) continue;
 627     //
 628     Bool_t res=kTRUE;
 629     if (fUseMaterial) res=AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp);
 630     else res=AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp,0,kFALSE,kFALSE);
 631
 632     if (!res) break;
 633
 634     if (TMath::Abs(track->GetZ())>kMaxZ) break;
 635     if (TMath::Abs(track->GetX())>kMaxR) break;
 636 //     if (TMath::Abs(track->GetZ())<kZcut)continue;
 637     //
 638     Double_t pointPos[2]={0,0};
 639     Double_t pointCov[3]={0,0,0};
 640     pointPos[0]=prot.GetY();//local y
 641     pointPos[1]=prot.GetZ();//local z
 642     pointCov[0]=prot.GetCov()[3];//simay^2
 643     pointCov[1]=prot.GetCov()[4];//sigmayz
 644     pointCov[2]=prot.GetCov()[5];//sigmaz^2
 645
 646     if (!track->Update(pointPos,pointCov)) {printf("no update\n"); break;}
 647   }
 648
 649   if (fUseMaterial) AliTrackerBase::PropagateTrackTo2(track,refX,kMass,5.,kTRUE,kMaxSnp);
 650   else AliTrackerBase::PropagateTrackTo2(track,refX,kMass,5.,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
 651
 652   // rotate fittet track to the frame of the original track and propagate to same reference
 653   track->Rotate(tr->GetAlpha());
 654
 655   if (fUseMaterial) AliTrackerBase::PropagateTrackTo2(track,refX,kMass,1.,kFALSE,kMaxSnp);
 656   else AliTrackerBase::PropagateTrackTo2(track,refX,kMass,1.,kFALSE,kMaxSnp,0,kFALSE,kFALSE);
 657
 658   return track;
 659 }
 660
 661
 662 //____________________________________________________________________________________
 663 AliExternalTrackParam* AliToyMCReconstruction::GetFittedTrackFromSeedAllClusters(const AliToyMCTrack *tr, const AliExternalTrackParam *seed, AliTPCtrackerSector *fInnerSectorArray, AliTPCtrackerSector *fOuterSectorArray)
 664 {
 665   //
 666   // Tracking for given seed on an array of clusters
 667   //
 668
 669   // create track
 670   AliExternalTrackParam *track = new AliExternalTrackParam(*seed);
 671
 672   Printf("Tracking NOT YET DONE");
 673
 674   return track;
 675 }
 676
 677
 678 //____________________________________________________________________________________
 679 void AliToyMCReconstruction::InitSpaceCharge()
 680 {
 681   //
 682   // Init the space charge map
 683   //
 684
 685   TString filename="$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps5_50kHz_precal.root";
 686   if (fTree) {
 687     TList *l=fTree->GetUserInfo();
 688     for (Int_t i=0; i<l->GetEntries(); ++i) {
 689       TString s(l->At(i)->GetName());
 690       if (s.Contains("SC_")) {
 691         filename=s;
 692         break;
 693       }
 694     }
 695   }
 696
 697   printf("Initialising the space charge map using the file: '%s'\n",filename.Data());
 698   TFile f(filename.Data());
 699   fSpaceCharge=(AliTPCSpaceCharge3D*)f.Get("map");
 700
 701   //   fSpaceCharge = new AliTPCSpaceCharge3D();
 702   //   fSpaceCharge->SetSCDataFileName("$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps10_50kHz.root");
 703   //   fSpaceCharge->SetOmegaTauT1T2(0.325,1,1); // Ne CO2
 704   // //   fSpaceCharge->SetOmegaTauT1T2(0.41,1,1.05); // Ar CO2
 705   //   fSpaceCharge->InitSpaceCharge3DDistortion();
 706
 707 }
 708
 709 //____________________________________________________________________________________
 710 Double_t AliToyMCReconstruction::GetVDrift() const
 711 {
 712   //
 713   //
 714   //
 715   return fTPCParam->GetDriftV();
 716 }
 717
 718 //____________________________________________________________________________________
 719 Double_t AliToyMCReconstruction::GetZLength(Int_t roc) const
 720 {
 721   //
 722   //
 723   //
 724   if (roc<0 || roc>71) return -1;
 725   return fTPCParam->GetZLength(roc);
 726 }
 727
 728 //____________________________________________________________________________________
 729 TTree* AliToyMCReconstruction::ConnectTrees (const char* files) {
 730   TString s=gSystem->GetFromPipe(Form("ls %s",files));
 731
 732   TTree *tFirst=0x0;
 733   TObjArray *arrFiles=s.Tokenize("\n");
 734
 735   for (Int_t ifile=0; ifile<arrFiles->GetEntriesFast(); ++ifile){
 736     TString name(arrFiles->At(ifile)->GetName());
 737
 738     TPRegexp reg(".*([0-9]_[0-9]_[0-9]_[0-9]{3}_[0-9]{2}).debug.root");
 739     TObjArray *arrMatch=0x0;
 740     arrMatch=reg.MatchS(name);
 741
 742     if (!tFirst) {
 743       TFile *f=TFile::Open(name.Data());
 744       if (!f) continue;
 745       TTree *t=(TTree*)f->Get("Tracks");
 746       if (!t) {
 747         delete f;
 748         continue;
 749       }
 750
 751       t->SetName(arrMatch->At(1)->GetName());
 752       tFirst=t;
 753     } else {
 754       tFirst->AddFriend(Form("t%s=Tracks",arrMatch->At(1)->GetName()), name.Data());
 755 //       tFirst->AddFriend(Form("t%d=Tracks",ifile), name.Data());
 756     }
 757   }
 758
 759   tFirst->GetListOfFriends()->Print();
 760   return tFirst;
 761 }
 762