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 <AliTPCCorrection.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 <AliTPCReconstructor.h>
  20 #include <AliTPCTransform.h>
  21 #include <AliTPCseed.h>
  22 #include <AliTPCtracker.h>
  23 #include <AliTPCtrackerSector.h>
  24
  25 #include "AliToyMCTrack.h"
  26 #include "AliToyMCEvent.h"
  27
  28 #include "AliToyMCReconstruction.h"
  29
  30 /*
  31
  32
  33
  34 */
  35
  36 //____________________________________________________________________________________
  37 AliToyMCReconstruction::AliToyMCReconstruction() : TObject()
  38 , fSeedingRow(140)
  39 , fSeedingDist(10)
  40 , fClusterType(0)
  41 , fCorrectionType(kNoCorrection)
  42 , fDoTrackFit(kTRUE)
  43 , fUseMaterial(kFALSE)
  44 , fIdealTracking(kFALSE)
  45 , fTime0(-1)
  46 , fCreateT0seed(kFALSE)
  47 , fStreamer(0x0)
  48 , fTree(0x0)
  49 , fEvent(0x0)
  50 , fTPCParam(0x0)
  51 , fTPCCorrection(0x0)
  52 , fkNSectorInner(18) // hard-coded to avoid loading the parameters before
  53 , fInnerSectorArray(0x0)
  54 , fkNSectorOuter(18) // hard-coded to avoid loading the parameters before
  55 , fOuterSectorArray(0x0)
  56 {
  57   //
  58   //  ctor
  59   //
  60   fTPCParam=AliTPCcalibDB::Instance()->GetParameters();
  61
  62 }
  63
  64 //____________________________________________________________________________________
  65 AliToyMCReconstruction::~AliToyMCReconstruction()
  66 {
  67   //
  68   //  dtor
  69   //
  70
  71   delete fStreamer;
  72 //   delete fTree;
  73 }
  74
  75 //____________________________________________________________________________________
  76 void AliToyMCReconstruction::RunReco(const char* file, Int_t nmaxEv)
  77 {
  78   //
  79   // Recostruction from associated clusters
  80   //
  81
  82   TFile f(file);
  83   if (!f.IsOpen() || f.IsZombie()) {
  84     printf("ERROR: couldn't open the file '%s'\n", file);
  85     return;
  86   }
  87
  88  fTree=(TTree*)f.Get("toyMCtree");
  89   if (!fTree) {
  90     printf("ERROR: couldn't read the 'toyMCtree' from file '%s'\n", file);
  91     return;
  92   }
  93
  94   fEvent=0x0;
  95   fTree->SetBranchAddress("event",&fEvent);
  96
  97   // read spacecharge from the Userinfo ot the tree
  98   InitSpaceCharge();
  99
 100   TString debugName=file;
 101   debugName.ReplaceAll(".root","");
 102   debugName.Append(Form(".%1d.%1d_%1d_%1d_%03d_%02d",
 103                         fUseMaterial,fIdealTracking,fClusterType,
 104                         Int_t(fCorrectionType),fSeedingRow,fSeedingDist));
 105   debugName.Append(".debug.root");
 106
 107   gSystem->Exec(Form("test -f %s && rm %s", debugName.Data(), debugName.Data()));
 108   if (!fStreamer) fStreamer=new TTreeSRedirector(debugName.Data());
 109
 110   gROOT->cd();
 111
 112   static AliExternalTrackParam dummySeedT0;
 113   static AliExternalTrackParam dummySeed;
 114   static AliExternalTrackParam dummyTrack;
 115
 116   AliExternalTrackParam t0seed;
 117   AliExternalTrackParam seed;
 118   AliExternalTrackParam track;
 119   AliExternalTrackParam tOrig;
 120
 121   AliExternalTrackParam *dummy;
 122
 123   Int_t maxev=fTree->GetEntries();
 124   if (nmaxEv>0&&nmaxEv<maxev) maxev=nmaxEv;
 125
 126   for (Int_t iev=0; iev<maxev; ++iev){
 127     printf("==============  Processing Event %6d =================\n",iev);
 128     fTree->GetEvent(iev);
 129     for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
 130       printf(" > ======  Processing Track %6d ========  \n",itr);
 131       const AliToyMCTrack *tr=fEvent->GetTrack(itr);
 132       tOrig = *tr;
 133
 134
 135       // set dummy
 136       t0seed    = dummySeedT0;
 137       seed      = dummySeed;
 138       track     = dummyTrack;
 139
 140       Float_t z0=fEvent->GetZ();
 141       Float_t t0=fEvent->GetT0();
 142
 143       Float_t vdrift=GetVDrift();
 144       Float_t zLength=GetZLength(0);
 145
 146       // crate time0 seed, steered by fCreateT0seed
 147 //       printf("t0 seed\n");
 148       fTime0=-1.;
 149       fCreateT0seed=kTRUE;
 150       dummy = GetSeedFromTrack(tr);
 151
 152       if (dummy) {
 153         t0seed = *dummy;
 154         delete dummy;
 155
 156         // crate real seed using the time 0 from the first seed
 157         // set fCreateT0seed now to false to get the seed in z coordinates
 158         fTime0 = t0seed.GetZ()-zLength/vdrift;
 159         fCreateT0seed = kFALSE;
 160 //         printf("seed (%.2g)\n",fTime0);
 161         dummy  = GetSeedFromTrack(tr);
 162         if (dummy) {
 163           seed = *dummy;
 164           delete dummy;
 165
 166           // create fitted track
 167           if (fDoTrackFit){
 168 //             printf("track\n");
 169             dummy = GetFittedTrackFromSeed(tr, &seed);
 170             track = *dummy;
 171             delete dummy;
 172           }
 173
 174           // propagate seed to 0
 175           const Double_t kMaxSnp = 0.85;
 176           const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 177           AliTrackerBase::PropagateTrackTo(&seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
 178
 179         }
 180       }
 181
 182       Int_t ctype(fCorrectionType);
 183
 184       if (fStreamer) {
 185         (*fStreamer) << "Tracks" <<
 186         "iev="         << iev             <<
 187         "z0="          << z0              <<
 188         "t0="          << t0              <<
 189         "fTime0="      << fTime0          <<
 190         "itr="         << itr             <<
 191         "clsType="     << fClusterType    <<
 192         "corrType="    << ctype           <<
 193         "seedRow="     << fSeedingRow     <<
 194         "seedDist="    << fSeedingDist    <<
 195         "vdrift="      << vdrift          <<
 196         "zLength="     << zLength         <<
 197         "t0seed.="     << &t0seed         <<
 198         "seed.="       << &seed           <<
 199         "track.="      << &track          <<
 200         "tOrig.="      << &tOrig          <<
 201         "\n";
 202       }
 203
 204
 205     }
 206   }
 207
 208   delete fStreamer;
 209   fStreamer=0x0;
 210
 211   delete fEvent;
 212   fEvent = 0x0;
 213
 214   delete fTree;
 215   fTree=0x0;
 216   f.Close();
 217 }
 218
 219
 220 //____________________________________________________________________________________
 221 void AliToyMCReconstruction::RunRecoAllClusters(const char* file, Int_t nmaxEv)
 222 {
 223   //
 224   // Reconstruction for seed from associated clusters, but array of clusters:
 225   // Step 1) Filling of cluster arrays
 226   // Step 2) Seeding from clusters associated to tracks
 227   // Step 3) Free track reconstruction using all clusters
 228   //
 229
 230   TFile f(file);
 231   if (!f.IsOpen() || f.IsZombie()) {
 232     printf("ERROR: couldn't open the file '%s'\n", file);
 233     return;
 234   }
 235
 236  fTree=(TTree*)f.Get("toyMCtree");
 237   if (!fTree) {
 238     printf("ERROR: couldn't read the 'toyMCtree' from file '%s'\n", file);
 239     return;
 240   }
 241
 242   fEvent=0x0;
 243   fTree->SetBranchAddress("event",&fEvent);
 244
 245   // read spacecharge from the Userinfo ot the tree
 246   InitSpaceCharge();
 247
 248   TString debugName=file;
 249   debugName.ReplaceAll(".root","");
 250   debugName.Append(Form(".%1d.%1d_%1d_%1d_%03d_%02d",
 251                         fUseMaterial,fIdealTracking,fClusterType,
 252                         Int_t(fCorrectionType),fSeedingRow,fSeedingDist));
 253   debugName.Append(".allClusters.debug.root");
 254
 255   gSystem->Exec(Form("test -f %s && rm %s", debugName.Data(), debugName.Data()));
 256   if (!fStreamer) fStreamer=new TTreeSRedirector(debugName.Data());
 257
 258   gROOT->cd();
 259
 260   static AliExternalTrackParam dummySeedT0;
 261   static AliExternalTrackParam dummySeed;
 262   static AliExternalTrackParam dummyTrack;
 263
 264   AliExternalTrackParam t0seed;
 265   AliExternalTrackParam seed;
 266   AliExternalTrackParam track;
 267   AliExternalTrackParam tOrig;
 268
 269   AliExternalTrackParam *dummy;
 270
 271   Int_t maxev=fTree->GetEntries();
 272   if (nmaxEv>0&&nmaxEv<maxev) maxev=nmaxEv;
 273
 274   // ===========================================================================================
 275   // Loop 1: Fill AliTPCtrackerSector structure
 276   // ===========================================================================================
 277   FillSectorStructure(maxev);
 278
 279   // settings (TODO: find the correct settings)
 280   AliTPCRecoParam *tpcRecoParam = new AliTPCRecoParam();
 281   tpcRecoParam->SetDoKinks(kFALSE);
 282   AliTPCcalibDB::Instance()->GetTransform()->SetCurrentRecoParam(tpcRecoParam);
 283   //tpcRecoParam->Print();
 284
 285   // need AliTPCReconstructor for parameter settings in AliTPCtracker
 286   AliTPCReconstructor *tpcRec   = new AliTPCReconstructor();
 287   tpcRec->SetRecoParam(tpcRecoParam);
 288
 289
 290   // ===========================================================================================
 291   // Loop 2: Seeding from clusters associated to tracks
 292   // TODO: Implement tracking from given seed!
 293   // ===========================================================================================
 294   for (Int_t iev=0; iev<maxev; ++iev){
 295     printf("==============  Processing Event %6d =================\n",iev);
 296     fTree->GetEvent(iev);
 297     for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
 298       printf(" > ======  Processing Track %6d ========  \n",itr);
 299       const AliToyMCTrack *tr=fEvent->GetTrack(itr);
 300       tOrig = *tr;
 301
 302
 303       // set dummy
 304       t0seed    = dummySeedT0;
 305       seed      = dummySeed;
 306       track     = dummyTrack;
 307
 308       Float_t z0=fEvent->GetZ();
 309       Float_t t0=fEvent->GetT0();
 310
 311       Float_t vdrift=GetVDrift();
 312       Float_t zLength=GetZLength(0);
 313
 314       Int_t nClus = 0;
 315
 316       // crate time0 seed, steered by fCreateT0seed
 317       printf("t0 seed\n");
 318       fTime0=-1.;
 319       fCreateT0seed=kTRUE;
 320       dummy = GetSeedFromTrack(tr);
 321
 322       if (dummy) {
 323         t0seed = *dummy;
 324         delete dummy;
 325
 326         // crate real seed using the time 0 from the first seed
 327         // set fCreateT0seed now to false to get the seed in z coordinates
 328         fTime0 = t0seed.GetZ()-zLength/vdrift;
 329         fCreateT0seed = kFALSE;
 330         printf("seed (%.2g)\n",fTime0);
 331         dummy  = GetSeedFromTrack(tr);
 332         if (dummy) {
 333           seed = *dummy;
 334           delete dummy;
 335
 336           // create fitted track
 337           if (fDoTrackFit){
 338             printf("track\n");
 339             dummy = GetFittedTrackFromSeedAllClusters(tr, &seed,nClus);
 340             track = *dummy;
 341             delete dummy;
 342           }
 343
 344           // propagate seed to 0
 345           const Double_t kMaxSnp = 0.85;
 346           const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 347           AliTrackerBase::PropagateTrackTo(&seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
 348
 349         }
 350       }
 351
 352       Int_t ctype(fCorrectionType);
 353
 354       if (fStreamer) {
 355         (*fStreamer) << "Tracks" <<
 356         "iev="         << iev             <<
 357         "z0="          << z0              <<
 358         "t0="          << t0              <<
 359         "fTime0="      << fTime0          <<
 360         "itr="         << itr             <<
 361         "clsType="     << fClusterType    <<
 362         "corrType="    << ctype           <<
 363         "seedRow="     << fSeedingRow     <<
 364         "seedDist="    << fSeedingDist    <<
 365         "vdrift="      << vdrift          <<
 366         "zLength="     << zLength         <<
 367         "nClus="       << nClus           <<
 368         "t0seed.="     << &t0seed         <<
 369         "seed.="       << &seed           <<
 370         "track.="      << &track          <<
 371         "tOrig.="      << &tOrig          <<
 372         "\n";
 373       }
 374
 375
 376     }
 377   }
 378
 379
 380   delete fStreamer;
 381   fStreamer=0x0;
 382
 383   delete fEvent;
 384   fEvent = 0x0;
 385
 386   delete fTree;
 387   fTree=0x0;
 388   f.Close();
 389 }
 390
 391 //____________________________________________________________________________________
 392 void AliToyMCReconstruction::RunRecoAllClustersStandardTracking(const char* file, Int_t nmaxEv)
 393 {
 394   //
 395   // Reconstruction for seed from associated clusters, but array of clusters
 396   // Step 1) Filling of cluster arrays
 397   // Step 2) Use the standard tracking: AliTPCtracker::Clusters2Tracks();
 398   //
 399
 400   TFile f(file);
 401   if (!f.IsOpen() || f.IsZombie()) {
 402     printf("ERROR: couldn't open the file '%s'\n", file);
 403     return;
 404   }
 405
 406  fTree=(TTree*)f.Get("toyMCtree");
 407   if (!fTree) {
 408     printf("ERROR: couldn't read the 'toyMCtree' from file '%s'\n", file);
 409     return;
 410   }
 411
 412   fEvent=0x0;
 413   fTree->SetBranchAddress("event",&fEvent);
 414
 415   // read spacecharge from the Userinfo ot the tree
 416   InitSpaceCharge();
 417
 418   TString debugName=file;
 419   debugName.ReplaceAll(".root","");
 420   debugName.Append(Form(".%1d.%1d_%1d_%1d_%03d_%02d",
 421                         fUseMaterial,fIdealTracking,fClusterType,
 422                         Int_t(fCorrectionType),fSeedingRow,fSeedingDist));
 423   debugName.Append(".allClusters.debug.root");
 424
 425   gSystem->Exec(Form("test -f %s && rm %s", debugName.Data(), debugName.Data()));
 426   if (!fStreamer) fStreamer=new TTreeSRedirector(debugName.Data());
 427
 428   gROOT->cd();
 429
 430   static AliExternalTrackParam dummySeedT0;
 431   static AliExternalTrackParam dummySeed;
 432   static AliExternalTrackParam dummyTrack;
 433
 434   AliExternalTrackParam t0seed;
 435   AliExternalTrackParam seed;
 436   AliExternalTrackParam track;
 437   AliExternalTrackParam tOrig;
 438
 439   AliExternalTrackParam *dummy;
 440   AliTPCseed            *seedBest;
 441   AliTPCseed            *seedTmp;
 442   AliTPCclusterMI       *cluster;
 443
 444   Int_t maxev=fTree->GetEntries();
 445   if (nmaxEv>0&&nmaxEv<maxev) maxev=nmaxEv;
 446
 447
 448   // ===========================================================================================
 449   // Loop 1: Fill AliTPCtrackerSector structure
 450   // ===========================================================================================
 451   FillSectorStructure(maxev);
 452
 453   // ===========================================================================================
 454   // Loop 2: Use the TPC tracker for seeding (MakeSeeds2)
 455   // TODO: - check tracking configuration
 456   //       - add clusters and original tracks to output (how?)
 457   // ===========================================================================================
 458
 459   // settings (TODO: find the correct settings)
 460   AliTPCRecoParam *tpcRecoParam = new AliTPCRecoParam();
 461   tpcRecoParam->SetDoKinks(kFALSE);
 462   AliTPCcalibDB::Instance()->GetTransform()->SetCurrentRecoParam(tpcRecoParam);
 463   //tpcRecoParam->Print();
 464
 465   // need AliTPCReconstructor for parameter settings in AliTPCtracker
 466   AliTPCReconstructor *tpcRec   = new AliTPCReconstructor();
 467   tpcRec->SetRecoParam(tpcRecoParam);
 468
 469   // AliTPCtracker
 470   AliTPCtracker *tpcTracker = new AliTPCtracker(fTPCParam);
 471   tpcTracker->SetDebug(10);
 472
 473   // set sector arrays
 474   tpcTracker->SetTPCtrackerSectors(fInnerSectorArray,fOuterSectorArray);
 475   tpcTracker->LoadInnerSectors();
 476   tpcTracker->LoadOuterSectors();
 477
 478   // seeding
 479   static TObjArray arrTracks;
 480   TObjArray * arr    = &arrTracks;
 481   TObjArray * seeds  = new TObjArray;
 482
 483   // cuts for seeding
 484   Float_t cuts[4];
 485   cuts[0]=0.0070;  // cuts[0]   - fP4 cut
 486   cuts[1] = 1.5;   // cuts[1]   - tan(phi) cut
 487   cuts[2] = 3.;    // cuts[2]   - zvertex cut
 488   cuts[3] = 3.;    // cuts[3]   - fP3 cut
 489
 490   // rows for seeding
 491   Int_t lowerRow = fSeedingRow;
 492   Int_t upperRow = fSeedingRow+2*fSeedingDist;
 493   const AliTPCROC * roc      = AliTPCROC::Instance();
 494   const Int_t kNRowsInnerTPC = roc->GetNRows(0);
 495   const Int_t kNRowsTPC      = kNRowsInnerTPC + roc->GetNRows(36);
 496   if(lowerRow < kNRowsInnerTPC){
 497     Printf("Seeding row requested (%d) is lower than kNRowsInnerTPC --> use %d",lowerRow,kNRowsInnerTPC);
 498     lowerRow = kNRowsInnerTPC;
 499     upperRow = lowerRow + 20;
 500   }
 501   if(upperRow >= kNRowsTPC){
 502     Printf("Seeding row requested (%d) is larger than kNRowsTPC --> use %d",upperRow,kNRowsTPC-1);
 503     upperRow = kNRowsTPC-1;
 504     lowerRow = upperRow-20;
 505   }
 506
 507   // do the seeding
 508   for (Int_t sec=0;sec<fkNSectorOuter;sec++){
 509     //
 510     tpcTracker->MakeSeeds3(arr, sec,upperRow,lowerRow,cuts,-1,1);
 511     tpcTracker->SumTracks(seeds,arr);
 512     tpcTracker->SignClusters(seeds,3.0,3.0);
 513   }
 514
 515   Printf("After seeding we have %d tracks",seeds->GetEntriesFast());
 516
 517   // Loop over all input tracks and connect to found seeds
 518   for (Int_t iev=0; iev<maxev; ++iev){
 519     printf("==============  Fill Tracks: Processing Event %6d  =================\n",iev);
 520     fTree->GetEvent(iev);
 521     for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
 522       printf(" > ======  Fill Tracks: Processing Track %6d  ========  \n",itr);
 523       const AliToyMCTrack *tr=fEvent->GetTrack(itr);
 524       tOrig = *tr;
 525
 526       Int_t nClus=0;
 527       Float_t z0=fEvent->GetZ();
 528       Float_t t0=fEvent->GetT0();
 529       Float_t vdrift=GetVDrift();
 530       Float_t zLength=GetZLength(0);
 531
 532       // find the corresponding seed (and track)
 533       Int_t trackID          = tr->GetUniqueID();
 534       Int_t idxSeed          = 0;
 535       Int_t nSeeds           = 0;
 536       Int_t nSeedClusters    = 0;
 537       Int_t nSeedClustersTmp = 0;
 538       for(Int_t iSeeds = 0; iSeeds < seeds->GetEntriesFast(); ++iSeeds){
 539
 540         seedTmp = (AliTPCseed*)(seeds->At(iSeeds));
 541         nSeedClustersTmp = 0;
 542         for(Int_t iRow = seedTmp->GetRow(); iRow < seedTmp->GetRow() + seedTmp->GetNumberOfClustersIndices(); iRow++ ){
 543
 544           cluster = seedTmp->GetClusterFast(iRow);
 545           if(cluster){
 546             if(cluster->GetLabel(0)==trackID){
 547               nSeedClustersTmp++;
 548             }
 549           }
 550         }
 551
 552         // if number of corresponding clusters > 0,
 553         // increase nSeeds
 554         if(nSeedClustersTmp > 0){
 555           nSeeds++;
 556         }
 557
 558         // if number of corresponding clusters bigger than current nSeedClusters,
 559         // take this seed as the main one
 560         if(nSeedClustersTmp > nSeedClusters){
 561           idxSeed  = iSeeds;
 562           seedBest = seedTmp;
 563           nSeedClusters = nSeedClustersTmp;
 564         }
 565
 566       }
 567
 568       // cluster to track association
 569       if (nSeeds>0&&nSeedClusters>0) {
 570
 571         t0seed = (AliExternalTrackParam)*seedBest;
 572         fTime0 = t0seed.GetZ()-zLength/vdrift;
 573         printf("seed (%.2g): %d seeds with %d clusters\n",fTime0,nSeeds,nSeedClusters);
 574
 575         // cluster to track association for all good seeds
 576         // set fCreateT0seed to true to get the seed in time coordinates
 577         fCreateT0seed = kTRUE;
 578         dummy = ClusterToTrackAssociation(seedBest,trackID,nClus);
 579         track = *dummy;
 580         delete dummy;
 581
 582          // propagate track to 0
 583         const Double_t kMaxSnp = 0.85;
 584         const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 585         AliTrackerBase::PropagateTrackTo(&track,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
 586
 587       }
 588
 589       Int_t ctype(fCorrectionType);
 590
 591       if (fStreamer) {
 592         (*fStreamer) << "Tracks" <<
 593           "iev="         << iev             <<
 594           "z0="          << z0              <<
 595           "t0="          << t0              <<
 596           "fTime0="      << fTime0          <<
 597           "itr="         << itr             <<
 598           "clsType="     << fClusterType    <<
 599           "corrType="    << ctype           <<
 600           "seedRow="     << fSeedingRow     <<
 601           "seedDist="    << fSeedingDist    <<
 602           "vdrift="      << vdrift          <<
 603           "zLength="     << zLength         <<
 604           "nSeeds="      << nSeeds          <<
 605           "nSeedClusters="<<nSeedClusters   <<
 606           "nClus="       << nClus           <<
 607           "t0seed.="     << &t0seed         <<
 608           "track.="      << &track          <<
 609           "tOrig.="      << &tOrig          <<
 610           "\n";
 611       }
 612     }
 613   }
 614
 615
 616   delete fStreamer;
 617   fStreamer=0x0;
 618
 619   delete fEvent;
 620   fEvent = 0x0;
 621
 622   delete fTree;
 623   fTree=0x0;
 624   f.Close();
 625 }
 626
 627
 628 //____________________________________________________________________________________
 629 AliExternalTrackParam* AliToyMCReconstruction::GetSeedFromTrack(const AliToyMCTrack * const tr)
 630 {
 631   //
 632   // if we don't have a valid time0 informaion (fTime0) available yet
 633   // assume we create a seed for the time0 estimate
 634   //
 635
 636   // seed point informaion
 637   AliTrackPoint    seedPoint[3];
 638   const AliTPCclusterMI *seedCluster[3]={0x0,0x0,0x0};
 639
 640   // number of clusters to loop over
 641   const Int_t ncls=(fClusterType==0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
 642
 643   UChar_t nextSeedRow=fSeedingRow;
 644   Int_t   nseeds=0;
 645
 646   //assumes sorted clusters
 647   for (Int_t icl=0;icl<ncls;++icl) {
 648     const AliTPCclusterMI *cl=tr->GetSpacePoint(icl);
 649     if (fClusterType==1) cl=tr->GetDistortedSpacePoint(icl);
 650     if (!cl) continue;
 651     // use row in sector
 652     const UChar_t row=cl->GetRow() + 63*(cl->GetDetector()>35);
 653     // skip clusters without proper pad row
 654     if (row>200) continue;
 655
 656     //check seeding row
 657     // if we are in the last row and still miss a seed we use the last row
 658     //   even if the row spacing will not be equal
 659     if (row>=nextSeedRow || icl==ncls-1){
 660       seedCluster[nseeds]=cl;
 661       SetTrackPointFromCluster(cl, seedPoint[nseeds]);
 662       ++nseeds;
 663       nextSeedRow+=fSeedingDist;
 664
 665       if (nseeds==3) break;
 666     }
 667   }
 668
 669   // check we really have 3 seeds
 670   if (nseeds!=3) {
 671     AliError(Form("Seeding failed for parameters %d, %d\n",fSeedingDist,fSeedingRow));
 672     return 0x0;
 673   }
 674
 675   // do cluster correction for fCorrectionType:
 676   //   0 - no correction
 677   //   1 - TPC center
 678   //   2 - average eta
 679   //   3 - ideal
 680   // assign the cluster abs time as z component to all seeds
 681   for (Int_t iseed=0; iseed<3; ++iseed) {
 682     Float_t xyz[3]={0,0,0};
 683     seedPoint[iseed].GetXYZ(xyz);
 684     const Float_t r=TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
 685
 686     const Int_t sector=seedCluster[iseed]->GetDetector();
 687     const Int_t sign=1-2*((sector/18)%2);
 688
 689     if ( (fClusterType == 1) && (fCorrectionType != kNoCorrection) ) {
 690 //       printf("correction type: %d\n",(Int_t)fCorrectionType);
 691
 692       // the settings below are for the T0 seed
 693       // for known T0 the z position is already calculated in SetTrackPointFromCluster
 694       if ( fCreateT0seed ){
 695         if ( fCorrectionType == kTPCCenter  ) xyz[2] = 125.*sign;
 696         //!!! TODO: is this the correct association?
 697         if ( fCorrectionType == kAverageEta ) xyz[2] = TMath::Tan(45./2.*TMath::DegToRad())*r*sign;
 698       }
 699
 700       if ( fCorrectionType == kIdeal      ) xyz[2] = seedCluster[iseed]->GetZ();
 701
 702       //!!! TODO: to be replaced with the proper correction
 703       fTPCCorrection->CorrectPoint(xyz, seedCluster[iseed]->GetDetector());
 704     }
 705
 706     // after the correction set the time bin as z-Position in case of a T0 seed
 707     if ( fCreateT0seed )
 708       xyz[2]=seedCluster[iseed]->GetTimeBin();
 709
 710     seedPoint[iseed].SetXYZ(xyz);
 711   }
 712
 713   const Double_t kMaxSnp = 0.85;
 714   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 715
 716   AliExternalTrackParam *seed = AliTrackerBase::MakeSeed(seedPoint[0], seedPoint[1], seedPoint[2]);
 717   seed->ResetCovariance(10);
 718
 719   if (fCreateT0seed){
 720     // if fTime0 < 0 we assume that we create a seed for the T0 estimate
 721     AliTrackerBase::PropagateTrackTo(seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
 722     if (TMath::Abs(seed->GetX())>3) {
 723 //       printf("Could not propagate track to 0, %.2f, %.2f, %.2f\n",seed->GetX(),seed->GetAlpha(),seed->Pt());
 724     }
 725   }
 726
 727   return seed;
 728
 729 }
 730
 731
 732 //____________________________________________________________________________________
 733 void AliToyMCReconstruction::SetTrackPointFromCluster(const AliTPCclusterMI *cl, AliTrackPoint &p )
 734 {
 735   //
 736   // make AliTrackPoint out of AliTPCclusterMI
 737   //
 738
 739   if (!cl) return;
 740     Float_t xyz[3]={0.,0.,0.};
 741   //   ClusterToSpacePoint(cl,xyz);
 742   //   cl->GetGlobalCov(cov);
 743   //TODO: what to do with the covariance matrix???
 744   //TODO: the problem is that it is used in GetAngle in AliTrackPoint
 745   //TODO: which is used by AliTrackerBase::MakeSeed to get alpha correct ...
 746   //TODO: for the moment simply assign 1 permill squared
 747   // in AliTrackPoint the cov is xx, xy, xz, yy, yz, zz
 748   //   Float_t cov[6]={xyz[0]*xyz[0]*1e-6,xyz[0]*xyz[1]*1e-6,xyz[0]*xyz[2]*1e-6,
 749   //                   xyz[1]*xyz[1]*1e-6,xyz[1]*xyz[2]*1e-6,xyz[2]*xyz[2]*1e-6};
 750   //   cl->GetGlobalXYZ(xyz);
 751   //   cl->GetGlobalCov(cov);
 752   // voluem ID to add later ....
 753   //   p.SetXYZ(xyz);
 754   //   p.SetCov(cov);
 755   AliTrackPoint *tp=const_cast<AliTPCclusterMI*>(cl)->MakePoint();
 756   p=*tp;
 757   delete tp;
 758   //   cl->Print();
 759   //   p.Print();
 760   p.SetVolumeID(cl->GetDetector());
 761
 762
 763   if ( !fCreateT0seed && !fIdealTracking ) {
 764     p.GetXYZ(xyz);
 765     const Int_t sector=cl->GetDetector();
 766     const Int_t sign=1-2*((sector/18)%2);
 767     const Float_t zT0=( GetZLength(sector) - (cl->GetTimeBin()-fTime0)*GetVDrift() )*sign;
 768 //     printf(" z:  %.2f  %.2f\n",xyz[2],zT0);
 769     xyz[2]=zT0;
 770     p.SetXYZ(xyz);
 771   }
 772
 773
 774   //   p.Rotate(p.GetAngle()).Print();
 775 }
 776
 777 //____________________________________________________________________________________
 778 void AliToyMCReconstruction::ClusterToSpacePoint(const AliTPCclusterMI *cl, Float_t xyz[3])
 779 {
 780   //
 781   // convert the cluster to a space point in global coordinates
 782   //
 783   if (!cl) return;
 784   xyz[0]=cl->GetRow();
 785   xyz[1]=cl->GetPad();
 786   xyz[2]=cl->GetTimeBin(); // this will not be correct at all
 787   Int_t i[3]={0,cl->GetDetector(),cl->GetRow()};
 788   //   printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
 789   fTPCParam->Transform8to4(xyz,i);
 790   //   printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
 791   fTPCParam->Transform4to3(xyz,i);
 792   //   printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
 793   fTPCParam->Transform2to1(xyz,i);
 794   //   printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
 795 }
 796
 797 //____________________________________________________________________________________
 798 AliExternalTrackParam* AliToyMCReconstruction::GetFittedTrackFromSeed(const AliToyMCTrack *tr, const AliExternalTrackParam *seed)
 799 {
 800   //
 801   //
 802   //
 803
 804   // create track
 805   AliExternalTrackParam *track = new AliExternalTrackParam(*seed);
 806
 807   Int_t ncls=(fClusterType == 0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
 808
 809   const AliTPCROC * roc = AliTPCROC::Instance();
 810
 811   const Double_t kRTPC0  = roc->GetPadRowRadii(0,0);
 812   const Double_t kRTPC1  = roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
 813   const Double_t kMaxSnp = 0.85;
 814   const Double_t kMaxR   = 500.;
 815   const Double_t kMaxZ   = 500.;
 816
 817   //   const Double_t kMaxZ0=220;
 818 //   const Double_t kZcut=3;
 819
 820   const Double_t refX = tr->GetX();
 821
 822   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 823
 824   // loop over all other points and add to the track
 825   for (Int_t ipoint=ncls-1; ipoint>=0; --ipoint){
 826     AliTrackPoint pIn;
 827     const AliTPCclusterMI *cl=tr->GetSpacePoint(ipoint);
 828     if (fClusterType == 1) cl=tr->GetDistortedSpacePoint(ipoint);
 829     SetTrackPointFromCluster(cl, pIn);
 830     if (fCorrectionType != kNoCorrection){
 831       Float_t xyz[3]={0,0,0};
 832       pIn.GetXYZ(xyz);
 833 //       if ( fCorrectionType == kIdeal ) xyz[2] = cl->GetZ();
 834       fTPCCorrection->CorrectPoint(xyz, cl->GetDetector());
 835       pIn.SetXYZ(xyz);
 836     }
 837     // rotate the cluster to the local detector frame
 838     track->Rotate(((cl->GetDetector()%18)*20+10)*TMath::DegToRad());
 839     AliTrackPoint prot = pIn.Rotate(track->GetAlpha());   // rotate to the local frame - non distoted  point
 840     if (TMath::Abs(prot.GetX())<kRTPC0) continue;
 841     if (TMath::Abs(prot.GetX())>kRTPC1) continue;
 842     //
 843     Bool_t res=AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
 844
 845     if (!res) break;
 846
 847     if (TMath::Abs(track->GetZ())>kMaxZ) break;
 848     if (TMath::Abs(track->GetX())>kMaxR) break;
 849 //     if (TMath::Abs(track->GetZ())<kZcut)continue;
 850     //
 851     Double_t pointPos[2]={0,0};
 852     Double_t pointCov[3]={0,0,0};
 853     pointPos[0]=prot.GetY();//local y
 854     pointPos[1]=prot.GetZ();//local z
 855     pointCov[0]=prot.GetCov()[3];//simay^2
 856     pointCov[1]=prot.GetCov()[4];//sigmayz
 857     pointCov[2]=prot.GetCov()[5];//sigmaz^2
 858
 859     if (!track->Update(pointPos,pointCov)) {printf("no update\n"); break;}
 860   }
 861
 862   AliTrackerBase::PropagateTrackTo2(track,refX,kMass,5.,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
 863
 864   // rotate fittet track to the frame of the original track and propagate to same reference
 865   track->Rotate(tr->GetAlpha());
 866
 867   AliTrackerBase::PropagateTrackTo2(track,refX,kMass,1.,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
 868
 869   return track;
 870 }
 871
 872
 873 //____________________________________________________________________________________
 874 AliExternalTrackParam* AliToyMCReconstruction::GetFittedTrackFromSeedAllClusters(const AliToyMCTrack *tr, const AliExternalTrackParam *seed, Int_t &nClus)
 875 {
 876   //
 877   // Tracking for given seed on an array of clusters
 878   //
 879
 880   // create track
 881   AliExternalTrackParam *track = new AliExternalTrackParam(*seed);
 882
 883   const AliTPCROC * roc = AliTPCROC::Instance();
 884
 885   const Double_t kRTPC0    = roc->GetPadRowRadii(0,0);
 886   const Double_t kRTPC1    = roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
 887   const Int_t kNRowsTPC    = roc->GetNRows(0) + roc->GetNRows(36) - 1;
 888   const Int_t kIRowsTPC    = roc->GetNRows(0) - 1;
 889   const Double_t kMaxSnp   = 0.85;
 890   const Double_t kMaxR     = 500.;
 891   const Double_t kMaxZ     = 500.;
 892   const Double_t roady     = 100.;
 893   const Double_t roadz     = 100.;
 894
 895   const Double_t refX = tr->GetX();
 896
 897   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 898
 899   Int_t  secCur   = -1;
 900   UInt_t indexCur = 0;
 901   Double_t xCur, yCur, zCur = 0.;
 902
 903   Float_t vDrift = GetVDrift();
 904
 905   // first propagate seed to outermost row
 906   AliTrackerBase::PropagateTrackTo(track,kRTPC1,kMass,5,kFALSE,kMaxSnp);
 907
 908   // Loop over rows and find the cluster candidates
 909   for( Int_t iRow = kNRowsTPC; iRow >= 0; --iRow ){
 910
 911     // inner or outer sector
 912     Bool_t bInnerSector = kTRUE;
 913     if(iRow > kIRowsTPC) bInnerSector = kFALSE;
 914
 915     // nearest track point/cluster (to be found)
 916     AliTrackPoint nearestPoint;
 917     AliTPCclusterMI *nearestCluster = NULL;
 918
 919     // Inner Sector
 920     if(bInnerSector){
 921
 922       // Propagate to center of pad row and extract parameters
 923       AliTrackerBase::PropagateTrackTo(track,roc->GetPadRowRadii(0,iRow),kMass,5,kFALSE,kMaxSnp);
 924       xCur   = track->GetX();
 925       yCur   = track->GetY();
 926       zCur   = track->GetZ();
 927       if ( !fIdealTracking ) {
 928         zCur = zCur/vDrift + fTime0; // Look at time, not at z!
 929       }
 930       secCur = GetSector(track);
 931
 932       // Find the nearest cluster (TODO: correct road settings!)
 933       Printf("inner tracking here: x = %.2f, y = %.2f, z = %.2f (Row %d Sector %d)",xCur,yCur,zCur,iRow,secCur);
 934       nearestCluster = fInnerSectorArray[secCur%fkNSectorInner][iRow].FindNearest2(yCur,zCur,roady,roadz,indexCur);
 935
 936       // Move to next row if now cluster found
 937       if(!nearestCluster) continue;
 938       //Printf("Nearest Clusters = %d (of %d) ",indexCur,fInnerSectorArray[secCur%fkNSectorInner][iRow].GetN());
 939
 940     }
 941
 942     // Outer sector
 943     else{
 944
 945       // Propagate to center of pad row and extract parameters
 946       AliTrackerBase::PropagateTrackTo(track,roc->GetPadRowRadii(36,iRow-kIRowsTPC-1),kMass,5,kFALSE,kMaxSnp);
 947       xCur   = track->GetX();
 948       yCur   = track->GetY();
 949       zCur   = track->GetZ();
 950       if ( !fIdealTracking ) {
 951         zCur = zCur/vDrift + fTime0; // Look at time, not at z!
 952       }
 953       secCur = GetSector(track);
 954
 955       // Find the nearest cluster (TODO: correct road settings!)
 956       Printf("outer tracking here: x = %.2f, y = %.2f, z = %.2f (Row %d Sector %d)",xCur,yCur,zCur,iRow,secCur);
 957       nearestCluster = fOuterSectorArray[(secCur-fkNSectorInner*2)%fkNSectorOuter][iRow-kIRowsTPC-1].FindNearest2(yCur,zCur,roady,roadz,indexCur);
 958
 959       // Move to next row if now cluster found
 960       if(!nearestCluster) continue;
 961       //Printf("Nearest Clusters = %d (of %d)",indexCur,fOuterSectorArray[(secCur-fkNSectorInner*2)%fkNSectorOuter][iRow-kIRowsTPC-1].GetN());
 962
 963     }
 964
 965     // create track point from cluster
 966     SetTrackPointFromCluster(nearestCluster,nearestPoint);
 967
 968     //Printf("Track point = %.2f %.2f %.2f",nearestPoint.GetX(),nearestPoint.GetY(),nearestPoint.GetZ());
 969
 970     // correction
 971     // TODO: also correction when looking for the next cluster?
 972     if (fCorrectionType != kNoCorrection){
 973       Float_t xyz[3]={0,0,0};
 974       nearestPoint.GetXYZ(xyz);
 975       fTPCCorrection->CorrectPoint(xyz, nearestCluster->GetDetector());
 976       nearestPoint.SetXYZ(xyz);
 977     }
 978
 979     // rotate the cluster to the local detector frame
 980     track->Rotate(((nearestCluster->GetDetector()%18)*20+10)*TMath::DegToRad());
 981     AliTrackPoint prot = nearestPoint.Rotate(track->GetAlpha());   // rotate to the local frame - non distoted  point
 982     if (TMath::Abs(prot.GetX())<kRTPC0) continue;
 983     if (TMath::Abs(prot.GetX())>kRTPC1) continue;
 984
 985
 986     //Printf("Rotated Track point = %.2f %.2f %.2f",prot.GetX(),prot.GetY(),prot.GetZ());
 987
 988     // update track with the nearest track point
 989     Bool_t res=AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
 990
 991     if (!res) break;
 992
 993     if (TMath::Abs(track->GetZ())>kMaxZ) break;
 994     if (TMath::Abs(track->GetX())>kMaxR) break;
 995     //if (TMath::Abs(track->GetZ())<kZcut)continue;
 996
 997       Double_t pointPos[2]={0,0};
 998       Double_t pointCov[3]={0,0,0};
 999       pointPos[0]=prot.GetY();//local y
1000       pointPos[1]=prot.GetZ();//local z
1001       pointCov[0]=prot.GetCov()[3];//simay^2
1002       pointCov[1]=prot.GetCov()[4];//sigmayz
1003       pointCov[2]=prot.GetCov()[5];//sigmaz^2
1004
1005       if (!track->Update(pointPos,pointCov)) {printf("no update\n");}
1006
1007       ++nClus;
1008   }
1009
1010
1011   // propagation to refX
1012   AliTrackerBase::PropagateTrackTo(track,refX,kMass,5.,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
1013
1014   // rotate fittet track to the frame of the original track and propagate to same reference
1015   track->Rotate(tr->GetAlpha());
1016
1017   AliTrackerBase::PropagateTrackTo(track,refX,kMass,1.,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
1018
1019   Printf("We have %d clusters in this track!",nClus);
1020
1021   return track;
1022 }
1023
1024 //____________________________________________________________________________________
1025 AliExternalTrackParam* AliToyMCReconstruction::ClusterToTrackAssociation(const AliTPCseed *seed, Int_t trackID, Int_t &nClus)
1026 {
1027   //
1028   // Cluster to track association for given seed on an array of clusters
1029   //
1030
1031   // create track
1032   AliExternalTrackParam *track = new AliExternalTrackParam(*seed);
1033
1034   const AliTPCROC * roc = AliTPCROC::Instance();
1035
1036   const Double_t kRTPC0    = roc->GetPadRowRadii(0,0);
1037   const Double_t kRTPC1    = roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
1038   const Int_t kNRowsTPC    = roc->GetNRows(0) + roc->GetNRows(36) - 1;
1039   const Int_t kIRowsTPC    = roc->GetNRows(0) - 1;
1040   const Double_t kMaxSnp   = 0.85;
1041   const Double_t kMaxR     = 500.;
1042   const Double_t kMaxZ     = 500.;
1043   const Double_t roady     = 0.1;
1044   const Double_t roadz     = 0.01;
1045
1046   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1047
1048   Int_t  secCur, secOld   = -1;
1049   UInt_t indexCur = 0;
1050   Double_t xCur, yCur, zCur = 0.;
1051
1052   Float_t vDrift = GetVDrift();
1053
1054   // first propagate seed to outermost row
1055   AliTrackerBase::PropagateTrackTo(track,kRTPC1,kMass,5,kFALSE,kMaxSnp);
1056
1057   // Loop over rows and find the cluster candidates
1058   for( Int_t iRow = kNRowsTPC; iRow >= 0; --iRow ){
1059
1060     // inner or outer sector
1061     Bool_t bInnerSector = kTRUE;
1062     if(iRow > kIRowsTPC) bInnerSector = kFALSE;
1063
1064     // nearest track point/cluster (to be found)
1065     AliTrackPoint nearestPoint;
1066     AliTPCclusterMI *nearestCluster = NULL;
1067
1068     // Inner Sector
1069     if(bInnerSector){
1070
1071       // Propagate to center of pad row and extract parameters
1072       AliTrackerBase::PropagateTrackTo(track,roc->GetPadRowRadii(0,iRow),kMass,5,kFALSE,kMaxSnp);
1073       xCur   = track->GetX();
1074       yCur   = track->GetY();
1075       zCur   = track->GetZ();
1076       secCur = GetSector(track);
1077
1078       // Find the nearest cluster (TODO: correct road settings!)
1079       //Printf("inner tracking here: x = %.2f, y = %.2f, z = %.6f (Row %d Sector %d)",xCur,yCur,zCur,iRow,secCur);
1080       nearestCluster = fInnerSectorArray[secCur%fkNSectorInner][iRow].FindNearest2(yCur,zCur,roady,roadz,indexCur);
1081
1082       // Look again at -y if nothing found here and the sector changed with respect to the last found cluster
1083       // Increase also the road in this case
1084       if(!nearestCluster && secCur != secOld){
1085         //Printf("inner tracking here 2: x = %.2f, y = %.2f, z = %.6f (Row %d Sector %d)",xCur,-yCur,zCur,iRow,secCur);
1086         nearestCluster = fInnerSectorArray[secCur%fkNSectorInner][iRow].FindNearest2(-yCur,zCur,roady*100,roadz,indexCur);
1087       }
1088
1089       // Move to next row if now cluster found
1090       if(!nearestCluster) continue;
1091       //Printf("Nearest Clusters = %d (of %d) ",indexCur,fInnerSectorArray[secCur%fkNSectorInner][iRow].GetN());
1092
1093     }
1094
1095     // Outer sector
1096     else{
1097
1098       // Propagate to center of pad row and extract parameters
1099       AliTrackerBase::PropagateTrackTo(track,roc->GetPadRowRadii(36,iRow-kIRowsTPC-1),kMass,5,kFALSE,kMaxSnp);
1100       xCur   = track->GetX();
1101       yCur   = track->GetY();
1102       zCur   = track->GetZ();
1103       secCur = GetSector(track);
1104
1105       // Find the nearest cluster (TODO: correct road settings!)
1106       //Printf("outer tracking here: x = %.2f, y = %.2f, z = %.6f (Row %d Sector %d)",xCur,yCur,zCur,iRow,secCur);
1107       nearestCluster = fOuterSectorArray[(secCur-fkNSectorInner*2)%fkNSectorOuter][iRow-kIRowsTPC-1].FindNearest2(yCur,zCur,roady,roadz,indexCur);
1108
1109       // Look again at -y if nothing found here and the sector changed with respect to the last found cluster
1110       // Increase also the road in this case
1111       if(!nearestCluster && secCur != secOld){
1112         //Printf("outer tracking here 2: x = %.2f, y = %.2f, z = %.6f (Row %d Sector %d)",xCur,-yCur,zCur,iRow,secCur);
1113         nearestCluster = fOuterSectorArray[(secCur-fkNSectorInner*2)%fkNSectorOuter][iRow-kIRowsTPC-1].FindNearest2(-yCur,zCur,roady*100,roadz,indexCur);
1114       }
1115
1116
1117       // Move to next row if now cluster found
1118       if(!nearestCluster) continue;
1119       //Printf("Nearest Clusters = %d (of %d)",indexCur,fOuterSectorArray[(secCur-fkNSectorInner*2)%fkNSectorOuter][iRow-kIRowsTPC-1].GetN());
1120
1121     }
1122
1123     // create track point from cluster
1124     SetTrackPointFromCluster(nearestCluster,nearestPoint);
1125
1126     //Printf("Track point = %.2f %.2f %.2f",nearestPoint.GetX(),nearestPoint.GetY(),nearestPoint.GetZ());
1127
1128     // rotate the cluster to the local detector frame
1129     track->Rotate(((nearestCluster->GetDetector()%18)*20+10)*TMath::DegToRad());
1130     AliTrackPoint prot = nearestPoint.Rotate(track->GetAlpha());   // rotate to the local frame - non distoted  point
1131     if (TMath::Abs(prot.GetX())<kRTPC0) continue;
1132     if (TMath::Abs(prot.GetX())>kRTPC1) continue;
1133
1134
1135     //Printf("Rotated Track point = %.2f %.2f %.2f",prot.GetX(),prot.GetY(),prot.GetZ());
1136
1137     // update track with the nearest track point
1138     Bool_t res=AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
1139
1140     if (!res) break;
1141
1142     if (TMath::Abs(track->GetZ())>kMaxZ) break;
1143     if (TMath::Abs(track->GetX())>kMaxR) break;
1144     //if (TMath::Abs(track->GetZ())<kZcut)continue;
1145
1146       Double_t pointPos[2]={0,0};
1147       Double_t pointCov[3]={0,0,0};
1148       pointPos[0]=prot.GetY();//local y
1149       pointPos[1]=prot.GetZ();//local z
1150       pointCov[0]=prot.GetCov()[3];//simay^2
1151       pointCov[1]=prot.GetCov()[4];//sigmayz
1152       pointCov[2]=prot.GetCov()[5];//sigmaz^2
1153
1154       if (!track->Update(pointPos,pointCov)) {printf("no update\n");}
1155       secOld = secCur;
1156
1157       //Printf("Cluster belongs to track = %d",nearestCluster->GetLabel(0));
1158
1159       // only count as associate cluster if it belongs to correct track!
1160       if(nearestCluster->GetLabel(0) == trackID)
1161         ++nClus;
1162   }
1163
1164   Printf("We have %d clusters in this track!",nClus);
1165
1166   return track;
1167 }
1168
1169
1170 //____________________________________________________________________________________
1171 void AliToyMCReconstruction::InitSpaceCharge()
1172 {
1173   //
1174   // Init the space charge map
1175   //
1176
1177   TString filename="$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps5_50kHz_precal.root";
1178   if (fTree) {
1179     TList *l=fTree->GetUserInfo();
1180     for (Int_t i=0; i<l->GetEntries(); ++i) {
1181       TString s(l->At(i)->GetName());
1182       if (s.Contains("SC_")) {
1183         filename=s;
1184         break;
1185       }
1186     }
1187   }
1188
1189   printf("Initialising the space charge map using the file: '%s'\n",filename.Data());
1190   TFile f(filename.Data());
1191   fTPCCorrection=(AliTPCCorrection*)f.Get("map");
1192
1193   //   fTPCCorrection = new AliTPCSpaceCharge3D();
1194   //   fTPCCorrection->SetSCDataFileName("$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps10_50kHz.root");
1195   //   fTPCCorrection->SetOmegaTauT1T2(0.325,1,1); // Ne CO2
1196   // //   fTPCCorrection->SetOmegaTauT1T2(0.41,1,1.05); // Ar CO2
1197   //   fTPCCorrection->InitSpaceCharge3DDistortion();
1198
1199 }
1200
1201 //____________________________________________________________________________________
1202 Double_t AliToyMCReconstruction::GetVDrift() const
1203 {
1204   //
1205   //
1206   //
1207   return fTPCParam->GetDriftV();
1208 }
1209
1210 //____________________________________________________________________________________
1211 Double_t AliToyMCReconstruction::GetZLength(Int_t roc) const
1212 {
1213   //
1214   //
1215   //
1216   if (roc<0 || roc>71) return -1;
1217   return fTPCParam->GetZLength(roc);
1218 }
1219
1220 //____________________________________________________________________________________
1221 TTree* AliToyMCReconstruction::ConnectTrees (const char* files) {
1222   TString s=gSystem->GetFromPipe(Form("ls %s",files));
1223
1224   TTree *tFirst=0x0;
1225   TObjArray *arrFiles=s.Tokenize("\n");
1226
1227   for (Int_t ifile=0; ifile<arrFiles->GetEntriesFast(); ++ifile){
1228     TString name(arrFiles->At(ifile)->GetName());
1229
1230     TPRegexp reg(".*([0-9]_[0-9]_[0-9]_[0-9]{3}_[0-9]{2}).debug.root");
1231     TObjArray *arrMatch=0x0;
1232     arrMatch=reg.MatchS(name);
1233     TString matchName;
1234     if (arrMatch && arrMatch->At(1)) matchName=arrMatch->At(1)->GetName();
1235     else matchName=Form("%02d",ifile);
1236     delete arrMatch;
1237
1238     if (!tFirst) {
1239       TFile *f=TFile::Open(name.Data());
1240       if (!f) continue;
1241       TTree *t=(TTree*)f->Get("Tracks");
1242       if (!t) {
1243         delete f;
1244         continue;
1245       }
1246
1247       t->SetName(matchName.Data());
1248       tFirst=t;
1249     } else {
1250       tFirst->AddFriend(Form("t%s=Tracks",matchName.Data()), name.Data());
1251 //       tFirst->AddFriend(Form("t%d=Tracks",ifile), name.Data());
1252     }
1253   }
1254
1255   tFirst->GetListOfFriends()->Print();
1256   return tFirst;
1257 }
1258
1259 //_____________________________________________________________________________
1260 Int_t AliToyMCReconstruction::LoadOuterSectors() {
1261   //-----------------------------------------------------------------
1262   // This function fills outer TPC sectors with clusters.
1263   // Copy and paste from AliTPCtracker
1264   //-----------------------------------------------------------------
1265   Int_t nrows = fOuterSectorArray->GetNRows();
1266   UInt_t index=0;
1267   for (Int_t sec = 0;sec<fkNSectorOuter;sec++)
1268     for (Int_t row = 0;row<nrows;row++){
1269       AliTPCtrackerRow*  tpcrow = &(fOuterSectorArray[sec%fkNSectorOuter][row]);
1270       Int_t sec2 = sec+2*fkNSectorInner;
1271       //left
1272       Int_t ncl = tpcrow->GetN1();
1273       while (ncl--) {
1274         AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1275         index=(((sec2<<8)+row)<<16)+ncl;
1276         tpcrow->InsertCluster(c,index);
1277       }
1278       //right
1279       ncl = tpcrow->GetN2();
1280       while (ncl--) {
1281         AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1282         index=((((sec2+fkNSectorOuter)<<8)+row)<<16)+ncl;
1283         tpcrow->InsertCluster(c,index);
1284       }
1285       //
1286       // write indexes for fast acces
1287       //
1288       for (Int_t i=0;i<510;i++)
1289         tpcrow->SetFastCluster(i,-1);
1290       for (Int_t i=0;i<tpcrow->GetN();i++){
1291         Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1292         tpcrow->SetFastCluster(zi,i);  // write index
1293       }
1294       Int_t last = 0;
1295       for (Int_t i=0;i<510;i++){
1296         if (tpcrow->GetFastCluster(i)<0)
1297           tpcrow->SetFastCluster(i,last);
1298         else
1299           last = tpcrow->GetFastCluster(i);
1300       }
1301     }
1302   return 0;
1303 }
1304
1305
1306 //_____________________________________________________________________________
1307 Int_t  AliToyMCReconstruction::LoadInnerSectors() {
1308   //-----------------------------------------------------------------
1309   // This function fills inner TPC sectors with clusters.
1310   // Copy and paste from AliTPCtracker
1311   //-----------------------------------------------------------------
1312   Int_t nrows = fInnerSectorArray->GetNRows();
1313   UInt_t index=0;
1314   for (Int_t sec = 0;sec<fkNSectorInner;sec++)
1315     for (Int_t row = 0;row<nrows;row++){
1316       AliTPCtrackerRow*  tpcrow = &(fInnerSectorArray[sec%fkNSectorInner][row]);
1317       //
1318       //left
1319       Int_t ncl = tpcrow->GetN1();
1320       while (ncl--) {
1321         AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1322         index=(((sec<<8)+row)<<16)+ncl;
1323         tpcrow->InsertCluster(c,index);
1324       }
1325       //right
1326       ncl = tpcrow->GetN2();
1327       while (ncl--) {
1328         AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1329         index=((((sec+fkNSectorInner)<<8)+row)<<16)+ncl;
1330         tpcrow->InsertCluster(c,index);
1331       }
1332       //
1333       // write indexes for fast acces
1334       //
1335       for (Int_t i=0;i<510;i++)
1336         tpcrow->SetFastCluster(i,-1);
1337       for (Int_t i=0;i<tpcrow->GetN();i++){
1338         Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1339         tpcrow->SetFastCluster(zi,i);  // write index
1340       }
1341       Int_t last = 0;
1342       for (Int_t i=0;i<510;i++){
1343         if (tpcrow->GetFastCluster(i)<0)
1344           tpcrow->SetFastCluster(i,last);
1345         else
1346           last = tpcrow->GetFastCluster(i);
1347       }
1348
1349     }
1350   return 0;
1351 }
1352
1353 //_____________________________________________________________________________
1354 Int_t  AliToyMCReconstruction::GetSector(AliExternalTrackParam *track) {
1355   //-----------------------------------------------------------------
1356   // This function returns the sector number for a given track
1357   //-----------------------------------------------------------------
1358
1359   Int_t sector = -1;
1360
1361   // get the sector number
1362   // rotate point to global system (track is already global!)
1363   Double_t xd[3];
1364   track->GetXYZ(xd);
1365   //track->Local2GlobalPosition(xd,track->GetAlpha());
1366
1367   // use TPCParams to get the sector number
1368   Float_t xyz[3] = {xd[0],xd[1],xd[2]};
1369   Int_t   i[3]   = {0,0,0};
1370   if(fTPCParam){
1371     sector  = fTPCParam->Transform0to1(xyz,i);
1372   }
1373
1374   return sector;
1375 }
1376
1377 //_____________________________________________________________________________
1378 void  AliToyMCReconstruction::FillSectorStructure(Int_t maxev) {
1379   //-----------------------------------------------------------------
1380   // This function fills the sector structure of AliToyMCReconstruction
1381   //-----------------------------------------------------------------
1382
1383   // cluster array of all sectors
1384   fInnerSectorArray = new AliTPCtrackerSector[fkNSectorInner];
1385   fOuterSectorArray = new AliTPCtrackerSector[fkNSectorOuter];
1386
1387   for (Int_t i=0; i<fkNSectorInner; ++i) fInnerSectorArray[i].Setup(fTPCParam,0);
1388   for (Int_t i=0; i<fkNSectorOuter; ++i) fOuterSectorArray[i].Setup(fTPCParam,1);
1389
1390   Int_t count[72][96] = { {0} , {0} };
1391
1392   for (Int_t iev=0; iev<maxev; ++iev){
1393     printf("==============  Fill Clusters: Processing Event %6d  =================\n",iev);
1394     fTree->GetEvent(iev);
1395     for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
1396       printf(" > ======  Fill Clusters: Processing Track %6d ========  \n",itr);
1397       const AliToyMCTrack *tr=fEvent->GetTrack(itr);
1398
1399       // number of clusters to loop over
1400       const Int_t ncls=(fClusterType==0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
1401
1402       for(Int_t icl=0; icl<ncls; ++icl){
1403
1404         AliTPCclusterMI *cl=const_cast<AliTPCclusterMI *>(tr->GetSpacePoint(icl));
1405         if (fClusterType==1) cl=const_cast<AliTPCclusterMI *>(tr->GetDistortedSpacePoint(icl));
1406         if (!cl) continue;
1407
1408         Int_t sec = cl->GetDetector();
1409         Int_t row = cl->GetRow();
1410
1411         // set cluster time to cluster Z (if not ideal tracking)
1412         if ( !fIdealTracking ) {
1413           cl->SetZ(cl->GetTimeBin());
1414         }
1415         //Printf("Fill clusters (sector %d row %d): %.2f %.2f %.2f %.2f",sec,row,cl->GetX(),cl->GetY(),cl->GetZ(),cl->GetTimeBin());
1416
1417         // fill arrays for inner and outer sectors (A/C side handled internally)
1418         if (sec<fkNSectorInner*2){
1419           fInnerSectorArray[sec%fkNSectorInner].InsertCluster(cl, count[sec][row], fTPCParam);
1420         }
1421         else{
1422           fOuterSectorArray[(sec-fkNSectorInner*2)%fkNSectorOuter].InsertCluster(cl, count[sec][row], fTPCParam);
1423         }
1424
1425         ++count[sec][row];
1426       }
1427     }
1428   }
1429
1430   // fill the arrays completely
1431   // LoadOuterSectors();
1432   // LoadInnerSectors();
1433
1434   // // check the arrays
1435   // for (Int_t i=0; i<fkNSectorInner; ++i){
1436   //   for (Int_t j=0; j<fInnerSectorArray[i].GetNRows(); ++j){
1437   //     if(fInnerSectorArray[i][j].GetN()>0){
1438   //    Printf("Inner: Sector %d Row %d : %d",i,j,fInnerSectorArray[i][j].GetN());
1439   //     }
1440   //   }
1441   // }
1442   // for (Int_t i=0; i<fkNSectorInner; ++i){
1443   //   for (Int_t j=0; j<fOuterSectorArray[i].GetNRows(); ++j){
1444   //     if(fOuterSectorArray[i][j].GetN()>0){
1445   //    Printf("Outer: Sector %d Row %d : %d",i,j,fOuterSectorArray[i][j].GetN());
1446   //     }
1447   //   }
1448   // }
1449 }