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 #include <AliRieman.h>
  25
  26 #include "AliToyMCTrack.h"
  27 #include "AliToyMCEvent.h"
  28
  29 #include "AliToyMCReconstruction.h"
  30
  31 /*
  32
  33
  34
  35 */
  36
  37 //____________________________________________________________________________________
  38 AliToyMCReconstruction::AliToyMCReconstruction() : TObject()
  39 , fSeedingRow(140)
  40 , fSeedingDist(10)
  41 , fClusterType(0)
  42 , fCorrectionType(kNoCorrection)
  43 , fDoTrackFit(kTRUE)
  44 , fUseMaterial(kFALSE)
  45 , fIdealTracking(kFALSE)
  46 , fTime0(-1)
  47 , fCreateT0seed(kFALSE)
  48 , fStreamer(0x0)
  49 , fInputFile(0x0)
  50 , fTree(0x0)
  51 , fEvent(0x0)
  52 , fTPCParam(0x0)
  53 , fTPCCorrection(0x0)
  54 , fkNSectorInner(18) // hard-coded to avoid loading the parameters before
  55 , fInnerSectorArray(0x0)
  56 , fkNSectorOuter(18) // hard-coded to avoid loading the parameters before
  57 , fOuterSectorArray(0x0)
  58 , fAllClusters("AliTPCclusterMI",10000)
  59 , fMapTrackEvent(10000)
  60 , fMapTrackTrackInEvent(10000)
  61 , fIsAC(kFALSE)
  62 {
  63   //
  64   //  ctor
  65   //
  66   fTPCParam=AliTPCcalibDB::Instance()->GetParameters();
  67
  68 }
  69
  70 //____________________________________________________________________________________
  71 AliToyMCReconstruction::~AliToyMCReconstruction()
  72 {
  73   //
  74   //  dtor
  75   //
  76
  77   Cleanup();
  78 }
  79
  80 //____________________________________________________________________________________
  81 void AliToyMCReconstruction::RunReco(const char* file, Int_t nmaxEv)
  82 {
  83   //
  84   // Recostruction from associated clusters
  85   //
  86
  87   ConnectInputFile(file);
  88   if (!fTree) return;
  89
  90   // read spacecharge from the Userinfo ot the tree
  91   InitSpaceCharge();
  92
  93   TString debugName=file;
  94   debugName.ReplaceAll(".root","");
  95   debugName.Append(Form(".%1d.%1d_%1d_%1d_%03d_%02d",
  96                         fUseMaterial,fIdealTracking,fClusterType,
  97                         Int_t(fCorrectionType),fSeedingRow,fSeedingDist));
  98   debugName.Append(".debug.root");
  99
 100   gSystem->Exec(Form("test -f %s && rm %s", debugName.Data(), debugName.Data()));
 101   if (!fStreamer) fStreamer=new TTreeSRedirector(debugName.Data());
 102
 103   gROOT->cd();
 104
 105   static AliExternalTrackParam dummySeedT0;
 106   static AliExternalTrackParam dummySeed;
 107   static AliExternalTrackParam dummyTrack;
 108
 109   AliExternalTrackParam t0seed;
 110   AliExternalTrackParam seed;
 111   AliExternalTrackParam track;
 112   AliExternalTrackParam tOrig;
 113
 114   AliExternalTrackParam *dummy;
 115
 116   Int_t maxev=fTree->GetEntries();
 117   if (nmaxEv>0&&nmaxEv<maxev) maxev=nmaxEv;
 118
 119   for (Int_t iev=0; iev<maxev; ++iev){
 120     printf("==============  Processing Event %6d =================\n",iev);
 121     fTree->GetEvent(iev);
 122     for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
 123 //       printf(" > ======  Processing Track %6d ========  \n",itr);
 124       const AliToyMCTrack *tr=fEvent->GetTrack(itr);
 125       tOrig = *tr;
 126
 127
 128       // set dummy
 129       t0seed    = dummySeedT0;
 130       seed      = dummySeed;
 131       track     = dummyTrack;
 132
 133       Float_t z0=fEvent->GetZ();
 134       Float_t t0=fEvent->GetT0();
 135
 136       Float_t vDrift=GetVDrift();
 137       Float_t zLength=GetZLength(0);
 138
 139       // crate time0 seed, steered by fCreateT0seed
 140 //       printf("t0 seed\n");
 141       fTime0=-1.;
 142       fCreateT0seed=kTRUE;
 143       dummy = GetSeedFromTrack(tr);
 144
 145       if (dummy) {
 146         t0seed = *dummy;
 147         delete dummy;
 148
 149         // crate real seed using the time 0 from the first seed
 150         // set fCreateT0seed now to false to get the seed in z coordinates
 151         fTime0 = t0seed.GetZ()-zLength/vDrift;
 152         fCreateT0seed = kFALSE;
 153 //         printf("seed (%.2g)\n",fTime0);
 154         dummy  = GetSeedFromTrack(tr);
 155         if (dummy) {
 156           seed = *dummy;
 157           delete dummy;
 158
 159           // create fitted track
 160           if (fDoTrackFit){
 161 //             printf("track\n");
 162             dummy = GetFittedTrackFromSeed(tr, &seed);
 163             track = *dummy;
 164             delete dummy;
 165           }
 166
 167           // propagate seed to 0
 168           const Double_t kMaxSnp = 0.85;
 169           const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 170           AliTrackerBase::PropagateTrackTo(&seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
 171
 172         }
 173       }
 174
 175       Int_t ctype(fCorrectionType);
 176
 177       if (fStreamer) {
 178         (*fStreamer) << "Tracks" <<
 179         "iev="         << iev             <<
 180         "z0="          << z0              <<
 181         "t0="          << t0              <<
 182         "fTime0="      << fTime0          <<
 183         "itr="         << itr             <<
 184         "clsType="     << fClusterType    <<
 185         "corrType="    << ctype           <<
 186         "seedRow="     << fSeedingRow     <<
 187         "seedDist="    << fSeedingDist    <<
 188         "vDrift="      << vDrift          <<
 189         "zLength="     << zLength         <<
 190         "t0seed.="     << &t0seed         <<
 191         "seed.="       << &seed           <<
 192         "track.="      << &track          <<
 193         "tOrig.="      << &tOrig          <<
 194         "\n";
 195       }
 196
 197
 198     }
 199   }
 200
 201   Cleanup();
 202 }
 203
 204
 205 //____________________________________________________________________________________
 206 void AliToyMCReconstruction::RunRecoAllClusters(const char* file, Int_t nmaxEv)
 207 {
 208   //
 209   // Reconstruction for seed from associated clusters, but array of clusters:
 210   // Step 1) Filling of cluster arrays
 211   // Step 2) Seeding from clusters associated to tracks
 212   // Step 3) Free track reconstruction using all clusters
 213   //
 214
 215   TFile f(file);
 216   if (!f.IsOpen() || f.IsZombie()) {
 217     printf("ERROR: couldn't open the file '%s'\n", file);
 218     return;
 219   }
 220
 221  fTree=(TTree*)f.Get("toyMCtree");
 222   if (!fTree) {
 223     printf("ERROR: couldn't read the 'toyMCtree' from file '%s'\n", file);
 224     return;
 225   }
 226
 227   fEvent=0x0;
 228   fTree->SetBranchAddress("event",&fEvent);
 229
 230   // read spacecharge from the Userinfo ot the tree
 231   InitSpaceCharge();
 232
 233   TString debugName=file;
 234   debugName.ReplaceAll(".root","");
 235   debugName.Append(Form(".%1d.%1d_%1d_%1d_%03d_%02d",
 236                         fUseMaterial,fIdealTracking,fClusterType,
 237                         Int_t(fCorrectionType),fSeedingRow,fSeedingDist));
 238   debugName.Append(".allClusters.debug.root");
 239
 240   gSystem->Exec(Form("test -f %s && rm %s", debugName.Data(), debugName.Data()));
 241   if (!fStreamer) fStreamer=new TTreeSRedirector(debugName.Data());
 242
 243   gROOT->cd();
 244
 245   static AliExternalTrackParam dummySeedT0;
 246   static AliExternalTrackParam dummySeed;
 247   static AliExternalTrackParam dummyTrack;
 248
 249   AliExternalTrackParam t0seed;
 250   AliExternalTrackParam seed;
 251   AliExternalTrackParam track;
 252   AliExternalTrackParam tOrig;
 253
 254   AliExternalTrackParam *dummy;
 255
 256   Int_t maxev=fTree->GetEntries();
 257   if (nmaxEv>0&&nmaxEv<maxev) maxev=nmaxEv;
 258
 259   // ===========================================================================================
 260   // Loop 1: Fill AliTPCtrackerSector structure
 261   // ===========================================================================================
 262   FillSectorStructure(maxev);
 263
 264   // settings (TODO: find the correct settings)
 265   AliTPCRecoParam *tpcRecoParam = new AliTPCRecoParam();
 266   tpcRecoParam->SetDoKinks(kFALSE);
 267   AliTPCcalibDB::Instance()->GetTransform()->SetCurrentRecoParam(tpcRecoParam);
 268   //tpcRecoParam->Print();
 269
 270   // need AliTPCReconstructor for parameter settings in AliTPCtracker
 271   AliTPCReconstructor *tpcRec   = new AliTPCReconstructor();
 272   tpcRec->SetRecoParam(tpcRecoParam);
 273
 274
 275   // ===========================================================================================
 276   // Loop 2: Seeding from clusters associated to tracks
 277   // TODO: Implement tracking from given seed!
 278   // ===========================================================================================
 279   for (Int_t iev=0; iev<maxev; ++iev){
 280     printf("==============  Processing Event %6d =================\n",iev);
 281     fTree->GetEvent(iev);
 282     for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
 283       printf(" > ======  Processing Track %6d ========  \n",itr);
 284       const AliToyMCTrack *tr=fEvent->GetTrack(itr);
 285       tOrig = *tr;
 286
 287
 288       // set dummy
 289       t0seed    = dummySeedT0;
 290       seed      = dummySeed;
 291       track     = dummyTrack;
 292
 293       Float_t z0=fEvent->GetZ();
 294       Float_t t0=fEvent->GetT0();
 295
 296       Float_t vDrift=GetVDrift();
 297       Float_t zLength=GetZLength(0);
 298
 299       Int_t nClus = 0;
 300
 301       // crate time0 seed, steered by fCreateT0seed
 302       printf("t0 seed\n");
 303       fTime0=-1.;
 304       fCreateT0seed=kTRUE;
 305       dummy = GetSeedFromTrack(tr);
 306
 307       if (dummy) {
 308         t0seed = *dummy;
 309         delete dummy;
 310
 311         // crate real seed using the time 0 from the first seed
 312         // set fCreateT0seed now to false to get the seed in z coordinates
 313         fTime0 = t0seed.GetZ()-zLength/vDrift;
 314         fCreateT0seed = kFALSE;
 315         printf("seed (%.2g)\n",fTime0);
 316         dummy  = GetSeedFromTrack(tr);
 317         if (dummy) {
 318           seed = *dummy;
 319           delete dummy;
 320
 321           // create fitted track
 322           if (fDoTrackFit){
 323             printf("track\n");
 324             dummy = GetFittedTrackFromSeedAllClusters(tr, &seed,nClus);
 325             track = *dummy;
 326             delete dummy;
 327           }
 328
 329           // propagate seed to 0
 330           const Double_t kMaxSnp = 0.85;
 331           const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 332           AliTrackerBase::PropagateTrackTo(&seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
 333
 334         }
 335       }
 336
 337       Int_t ctype(fCorrectionType);
 338
 339       if (fStreamer) {
 340         (*fStreamer) << "Tracks" <<
 341         "iev="         << iev             <<
 342         "z0="          << z0              <<
 343         "t0="          << t0              <<
 344         "fTime0="      << fTime0          <<
 345         "itr="         << itr             <<
 346         "clsType="     << fClusterType    <<
 347         "corrType="    << ctype           <<
 348         "seedRow="     << fSeedingRow     <<
 349         "seedDist="    << fSeedingDist    <<
 350         "vDrift="      << vDrift          <<
 351         "zLength="     << zLength         <<
 352         "nClus="       << nClus           <<
 353         "t0seed.="     << &t0seed         <<
 354         "seed.="       << &seed           <<
 355         "track.="      << &track          <<
 356         "tOrig.="      << &tOrig          <<
 357         "\n";
 358       }
 359
 360
 361     }
 362   }
 363
 364
 365   delete fStreamer;
 366   fStreamer=0x0;
 367
 368   delete fEvent;
 369   fEvent = 0x0;
 370
 371   delete fTree;
 372   fTree=0x0;
 373   f.Close();
 374 }
 375
 376 //____________________________________________________________________________________
 377 void AliToyMCReconstruction::RunRecoAllClustersStandardTracking(const char* file, Int_t nmaxEv)
 378 {
 379   //
 380   // Reconstruction for seed from associated clusters, but array of clusters
 381   // Step 1) Filling of cluster arrays
 382   // Step 2) Use the standard tracking: AliTPCtracker::Clusters2Tracks();
 383   //
 384
 385   TFile f(file);
 386   if (!f.IsOpen() || f.IsZombie()) {
 387     printf("ERROR: couldn't open the file '%s'\n", file);
 388     return;
 389   }
 390
 391  fTree=(TTree*)f.Get("toyMCtree");
 392   if (!fTree) {
 393     printf("ERROR: couldn't read the 'toyMCtree' from file '%s'\n", file);
 394     return;
 395   }
 396
 397   fEvent=0x0;
 398   fTree->SetBranchAddress("event",&fEvent);
 399
 400   // read spacecharge from the Userinfo ot the tree
 401   InitSpaceCharge();
 402
 403   TString debugName=file;
 404   debugName.ReplaceAll(".root","");
 405   debugName.Append(Form(".%1d.%1d_%1d_%1d_%03d_%02d",
 406                         fUseMaterial,fIdealTracking,fClusterType,
 407                         Int_t(fCorrectionType),fSeedingRow,fSeedingDist));
 408   debugName.Append(".allClusters.debug.root");
 409
 410   gSystem->Exec(Form("test -f %s && rm %s", debugName.Data(), debugName.Data()));
 411   if (!fStreamer) fStreamer=new TTreeSRedirector(debugName.Data());
 412
 413   gROOT->cd();
 414
 415   AliExternalTrackParam t0seed;
 416   AliExternalTrackParam seed;
 417   AliExternalTrackParam track;
 418   AliExternalTrackParam tOrig;
 419   AliToyMCTrack tOrigToy;
 420
 421   AliExternalTrackParam *dummy;
 422   AliTPCseed            *seedBest;
 423   AliTPCseed            *seedTmp;
 424   AliTPCclusterMI       *cluster;
 425
 426   Int_t maxev=fTree->GetEntries();
 427   if (nmaxEv>0&&nmaxEv<maxev) maxev=nmaxEv;
 428
 429
 430   // ===========================================================================================
 431   // Loop 1: Fill AliTPCtrackerSector structure
 432   // ===========================================================================================
 433   FillSectorStructure(maxev);
 434
 435   // ===========================================================================================
 436   // Loop 2: Use the TPC tracker for seeding (MakeSeeds3)
 437   // TODO: - check tracking configuration
 438   //       - add clusters and original tracks to output (how?)
 439   // ===========================================================================================
 440
 441   // settings (TODO: find the correct settings)
 442   AliTPCRecoParam *tpcRecoParam = new AliTPCRecoParam();
 443   tpcRecoParam->SetDoKinks(kFALSE);
 444   AliTPCcalibDB::Instance()->GetTransform()->SetCurrentRecoParam(tpcRecoParam);
 445   //tpcRecoParam->Print();
 446
 447   // need AliTPCReconstructor for parameter settings in AliTPCtracker
 448   AliTPCReconstructor *tpcRec   = new AliTPCReconstructor();
 449   tpcRec->SetRecoParam(tpcRecoParam);
 450
 451   // AliTPCtracker
 452   AliTPCtracker *tpcTracker = new AliTPCtracker(fTPCParam);
 453   tpcTracker->SetDebug(10);
 454
 455   // set sector arrays
 456   tpcTracker->SetTPCtrackerSectors(fInnerSectorArray,fOuterSectorArray);
 457   tpcTracker->LoadInnerSectors();
 458   tpcTracker->LoadOuterSectors();
 459
 460   // seeding
 461   static TObjArray arrTracks;
 462   TObjArray * arr    = &arrTracks;
 463   TObjArray * seeds  = new TObjArray;
 464
 465   // cuts for seeding
 466 //   Float_t cuts[4];
 467 //   cuts[0]=0.0070;  // cuts[0]   - fP4 cut
 468 //   cuts[1] = 1.5;   // cuts[1]   - tan(phi) cut
 469 //   cuts[2] = 3.;    // cuts[2]   - zvertex cut
 470 //   cuts[3] = 3.;    // cuts[3]   - fP3 cut
 471
 472   // rows for seeding
 473   Int_t lowerRow = fSeedingRow;
 474   Int_t upperRow = fSeedingRow+2*fSeedingDist;
 475   const AliTPCROC * roc      = AliTPCROC::Instance();
 476   const Int_t kNRowsInnerTPC = roc->GetNRows(0);
 477   const Int_t kNRowsTPC      = kNRowsInnerTPC + roc->GetNRows(36);
 478   if(lowerRow < kNRowsInnerTPC){
 479     Printf("Seeding row requested (%d) is lower than kNRowsInnerTPC --> use %d",lowerRow,kNRowsInnerTPC);
 480     lowerRow = kNRowsInnerTPC;
 481     upperRow = lowerRow + 20;
 482   }
 483   if(upperRow >= kNRowsTPC){
 484     Printf("Seeding row requested (%d) is larger than kNRowsTPC --> use %d",upperRow,kNRowsTPC-1);
 485     upperRow = kNRowsTPC-1;
 486     lowerRow = upperRow-20;
 487   }
 488
 489   // do the seeding
 490   for (Int_t sec=0;sec<fkNSectorOuter;sec++){
 491     //
 492     //tpcTracker->MakeSeeds3(arr, sec,upperRow,lowerRow,cuts,-1,1);
 493 //     MakeSeeds(arr, sec,upperRow,lowerRow); // own function (based on TLinearFitter)
 494     MakeSeeds2(arr, sec,upperRow,lowerRow); // own function (based on TLinearFitter)
 495     //tpcTracker->SumTracks(seeds,arr);
 496     //tpcTracker->SignClusters(seeds,3.0,3.0);
 497
 498   }
 499
 500   Printf("After seeding we have %d tracks",seeds->GetEntriesFast());
 501
 502   // Standard tracking
 503   tpcTracker->SetSeeds(seeds);
 504   tpcTracker->PropagateForward();
 505   Printf("After trackinging we have %d tracks",seeds->GetEntriesFast());
 506 return;
 507
 508   // Loop over all input tracks and connect to found seeds
 509   for (Int_t iev=0; iev<maxev; ++iev){
 510     printf("==============  Fill Tracks: Processing Event %6d  =================\n",iev);
 511     fTree->GetEvent(iev);
 512     for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
 513       printf(" > ======  Fill Tracks: Processing Track %6d  ========  \n",itr);
 514       const AliToyMCTrack *tr=fEvent->GetTrack(itr);
 515       tOrig = *tr;
 516       tOrigToy = *tr;
 517
 518       Float_t z0=fEvent->GetZ();
 519       Float_t t0=fEvent->GetT0();
 520       Float_t vDrift=GetVDrift();
 521       Float_t zLength=GetZLength(0);
 522
 523       // find the corresponding seed (and track)
 524       Int_t trackID            = tr->GetUniqueID();
 525       Int_t nClustersMC        = tr->GetNumberOfSpacePoints();      // number of findable clusters (ideal)
 526       if(fClusterType==1)
 527             nClustersMC        = tr->GetNumberOfDistSpacePoints();  // number of findable clusters (distorted)
 528 //       Int_t idxSeed            = 0; // index of best seed (best is with maximum number of clusters with correct ID)
 529       Int_t nSeeds             = 0; // number of seeds for MC track
 530       Int_t nSeedClusters      = 0; // number of clusters for best seed
 531       Int_t nSeedClustersTmp   = 0; // number of clusters for current seed
 532       Int_t nSeedClustersID    = 0; // number of clusters with correct ID for best seed
 533       Int_t nSeedClustersIDTmp = 0; // number of clusters with correct ID for current seed
 534       for(Int_t iSeeds = 0; iSeeds < seeds->GetEntriesFast(); ++iSeeds){
 535
 536         // set current seed and reset counters
 537         seedTmp = (AliTPCseed*)(seeds->At(iSeeds));
 538         nSeedClustersTmp   = 0;
 539         nSeedClustersIDTmp = 0;
 540
 541         if(!seedTmp) continue;
 542
 543         // loop over all rows
 544         for(Int_t iRow = seedTmp->GetRow(); iRow < seedTmp->GetRow() + seedTmp->GetNumberOfClustersIndices(); iRow++ ){
 545
 546           // get cluster and increment counters
 547           cluster = seedTmp->GetClusterFast(iRow);
 548           if(cluster){
 549             nSeedClustersTmp++;
 550             if(cluster->GetLabel(0)==trackID){
 551               nSeedClustersIDTmp++;
 552             }
 553           }
 554         }
 555
 556         // if number of corresponding clusters > 0,
 557         // increase nSeeds
 558         if(nSeedClustersTmp > 0){
 559           nSeeds++;
 560         }
 561
 562         // if number of corresponding clusters bigger than current nSeedClusters,
 563         // take this seed as the best one
 564         if(nSeedClustersIDTmp > nSeedClustersID){
 565 //        idxSeed  = iSeeds;
 566           seedBest = seedTmp;
 567           nSeedClusters   = nSeedClustersTmp;   // number of correctly assigned clusters
 568           nSeedClustersID = nSeedClustersIDTmp; // number of all clusters
 569         }
 570
 571       }
 572
 573       // cluster to track association (commented out, when used standard tracking)
 574       if (nSeeds>0&&nSeedClusters>0) {
 575         t0seed = (AliExternalTrackParam)*seedBest;
 576 //              fTime0 = t0seed.GetZ()-zLength/vDrift;
 577         // get the refitted track from the seed
 578         // this will also set the fTime0 from the seed extrapolation
 579         dummy=GetRefittedTrack(*seedBest);
 580         track=*dummy;
 581         delete dummy;
 582         //printf("seed (%.2g): %d seeds with %d clusters\n",fTime0,nSeeds,nSeedClusters);
 583
 584         //      // cluster to track association for all good seeds
 585         //      // set fCreateT0seed to true to get the seed in time coordinates
 586         //      fCreateT0seed = kTRUE;
 587         //      dummy = ClusterToTrackAssociation(seedBest,trackID,nClus);
 588
 589         //// propagate track to 0
 590         //const Double_t kMaxSnp = 0.85;
 591         //const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 592         //AliTrackerBase::PropagateTrackTo(&track,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
 593
 594       }
 595
 596       Int_t ctype(fCorrectionType);
 597
 598       if (fStreamer) {
 599         (*fStreamer) << "Tracks" <<
 600           "iev="             << iev             <<
 601           "z0="              << z0              <<
 602           "t0="              << t0              <<
 603           "fTime0="          << fTime0          <<
 604           "itr="             << itr             <<
 605           "clsType="         << fClusterType    <<
 606           "corrType="        << ctype           <<
 607           "seedRow="         << fSeedingRow     <<
 608           "seedDist="        << fSeedingDist    <<
 609           "vDrift="          << vDrift          <<
 610           "zLength="         << zLength         <<
 611           "nClustersMC="     << nClustersMC     <<
 612           "nSeeds="          << nSeeds          <<
 613           "nSeedClusters="   << nSeedClusters   <<
 614           "nSeedClustersID=" << nSeedClustersID <<
 615           "t0seed.="         << &t0seed         <<
 616           "track.="          << &track          <<
 617           "tOrig.="          << &tOrig          <<
 618           "tOrigToy.="       << &tOrigToy       <<
 619           "\n";
 620       }
 621     }
 622   }
 623
 624
 625   delete fStreamer;
 626   fStreamer=0x0;
 627
 628   delete fEvent;
 629   fEvent = 0x0;
 630
 631   delete fTree;
 632   fTree=0x0;
 633   f.Close();
 634 }
 635
 636
 637 //____________________________________________________________________________________
 638 AliExternalTrackParam* AliToyMCReconstruction::GetSeedFromTrack(const AliToyMCTrack * const tr)
 639 {
 640   //
 641   // if we don't have a valid time0 informaion (fTime0) available yet
 642   // assume we create a seed for the time0 estimate
 643   //
 644
 645   // seed point informaion
 646   AliTrackPoint    seedPoint[3];
 647   const AliTPCclusterMI *seedCluster[3]={0x0,0x0,0x0};
 648
 649   // number of clusters to loop over
 650   const Int_t ncls=(fClusterType==0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
 651
 652   UChar_t nextSeedRow=fSeedingRow;
 653   Int_t   nseeds=0;
 654
 655   //assumes sorted clusters
 656   for (Int_t icl=0;icl<ncls;++icl) {
 657     const AliTPCclusterMI *cl=tr->GetSpacePoint(icl);
 658     if (fClusterType==1) cl=tr->GetDistortedSpacePoint(icl);
 659     if (!cl) continue;
 660     // use row in sector
 661     const UChar_t row=cl->GetRow() + 63*(cl->GetDetector()>35);
 662     // skip clusters without proper pad row
 663     if (row>200) continue;
 664
 665     //check seeding row
 666     // if we are in the last row and still miss a seed we use the last row
 667     //   even if the row spacing will not be equal
 668     if (row>=nextSeedRow || icl==ncls-1){
 669       seedCluster[nseeds]=cl;
 670       SetTrackPointFromCluster(cl, seedPoint[nseeds]);
 671       ++nseeds;
 672       nextSeedRow+=fSeedingDist;
 673
 674       if (nseeds==3) break;
 675     }
 676   }
 677
 678   // check we really have 3 seeds
 679   if (nseeds!=3) {
 680     AliError(Form("Seeding failed for parameters %d, %d\n",fSeedingDist,fSeedingRow));
 681     return 0x0;
 682   }
 683
 684   // do cluster correction for fCorrectionType:
 685   //   0 - no correction
 686   //   1 - TPC center
 687   //   2 - average eta
 688   //   3 - ideal
 689   // assign the cluster abs time as z component to all seeds
 690   for (Int_t iseed=0; iseed<3; ++iseed) {
 691     Float_t xyz[3]={0,0,0};
 692     seedPoint[iseed].GetXYZ(xyz);
 693     const Float_t r=TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
 694
 695     const Int_t sector=seedCluster[iseed]->GetDetector();
 696     const Int_t sign=1-2*((sector/18)%2);
 697
 698     if ( (fClusterType == 1) && (fCorrectionType != kNoCorrection) ) {
 699 //       printf("correction type: %d\n",(Int_t)fCorrectionType);
 700
 701       // the settings below are for the T0 seed
 702       // for known T0 the z position is already calculated in SetTrackPointFromCluster
 703       if ( fCreateT0seed ){
 704         if ( fCorrectionType == kTPCCenter  ) xyz[2] = 125.*sign;
 705         //!!! TODO: is this the correct association?
 706         if ( fCorrectionType == kAverageEta ) xyz[2] = TMath::Tan(45./2.*TMath::DegToRad())*r*sign;
 707       }
 708
 709       if ( fCorrectionType == kIdeal      ) xyz[2] = seedCluster[iseed]->GetZ();
 710
 711       //!!! TODO: to be replaced with the proper correction
 712       fTPCCorrection->CorrectPoint(xyz, seedCluster[iseed]->GetDetector());
 713     }
 714
 715     // after the correction set the time bin as z-Position in case of a T0 seed
 716     if ( fCreateT0seed )
 717       xyz[2]=seedCluster[iseed]->GetTimeBin();
 718
 719     seedPoint[iseed].SetXYZ(xyz);
 720   }
 721
 722   const Double_t kMaxSnp = 0.85;
 723   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 724
 725   AliExternalTrackParam *seed = AliTrackerBase::MakeSeed(seedPoint[0], seedPoint[1], seedPoint[2]);
 726   seed->ResetCovariance(10);
 727
 728   if (fCreateT0seed){
 729     // if fTime0 < 0 we assume that we create a seed for the T0 estimate
 730     AliTrackerBase::PropagateTrackTo(seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
 731     if (TMath::Abs(seed->GetX())>3) {
 732 //       printf("Could not propagate track to 0, %.2f, %.2f, %.2f\n",seed->GetX(),seed->GetAlpha(),seed->Pt());
 733     }
 734   }
 735
 736   return seed;
 737
 738 }
 739
 740
 741 //____________________________________________________________________________________
 742 void AliToyMCReconstruction::SetTrackPointFromCluster(const AliTPCclusterMI *cl, AliTrackPoint &p )
 743 {
 744   //
 745   // make AliTrackPoint out of AliTPCclusterMI
 746   //
 747
 748   if (!cl) return;
 749     Float_t xyz[3]={0.,0.,0.};
 750   //   ClusterToSpacePoint(cl,xyz);
 751   //   cl->GetGlobalCov(cov);
 752   //TODO: what to do with the covariance matrix???
 753   //TODO: the problem is that it is used in GetAngle in AliTrackPoint
 754   //TODO: which is used by AliTrackerBase::MakeSeed to get alpha correct ...
 755   //TODO: for the moment simply assign 1 permill squared
 756   // in AliTrackPoint the cov is xx, xy, xz, yy, yz, zz
 757   //   Float_t cov[6]={xyz[0]*xyz[0]*1e-6,xyz[0]*xyz[1]*1e-6,xyz[0]*xyz[2]*1e-6,
 758   //                   xyz[1]*xyz[1]*1e-6,xyz[1]*xyz[2]*1e-6,xyz[2]*xyz[2]*1e-6};
 759   //   cl->GetGlobalXYZ(xyz);
 760   //   cl->GetGlobalCov(cov);
 761   // voluem ID to add later ....
 762   //   p.SetXYZ(xyz);
 763   //   p.SetCov(cov);
 764 //   AliTrackPoint *tp=const_cast<AliTPCclusterMI*>(cl)->MakePoint(p);
 765 //   p=*tp;
 766 //   delete tp;
 767   const_cast<AliTPCclusterMI*>(cl)->MakePoint(p);
 768   //   cl->Print();
 769   //   p.Print();
 770   p.SetVolumeID(cl->GetDetector());
 771
 772
 773   if ( !fCreateT0seed && !fIdealTracking ) {
 774     p.GetXYZ(xyz);
 775     const Int_t sector=cl->GetDetector();
 776     const Int_t sign=1-2*((sector/18)%2);
 777     const Float_t zT0=( GetZLength(sector) - (cl->GetTimeBin()-fTime0)*GetVDrift() )*sign;
 778 //     printf(" z:  %.2f  %.2f\n",xyz[2],zT0);
 779     xyz[2]=zT0;
 780     p.SetXYZ(xyz);
 781   }
 782
 783
 784   //   p.Rotate(p.GetAngle()).Print();
 785 }
 786
 787 //____________________________________________________________________________________
 788 void AliToyMCReconstruction::ClusterToSpacePoint(const AliTPCclusterMI *cl, Float_t xyz[3])
 789 {
 790   //
 791   // convert the cluster to a space point in global coordinates
 792   //
 793   if (!cl) return;
 794   xyz[0]=cl->GetRow();
 795   xyz[1]=cl->GetPad();
 796   xyz[2]=cl->GetTimeBin(); // this will not be correct at all
 797   Int_t i[3]={0,cl->GetDetector(),cl->GetRow()};
 798   //   printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
 799   fTPCParam->Transform8to4(xyz,i);
 800   //   printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
 801   fTPCParam->Transform4to3(xyz,i);
 802   //   printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
 803   fTPCParam->Transform2to1(xyz,i);
 804   //   printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
 805 }
 806
 807 //____________________________________________________________________________________
 808 AliExternalTrackParam* AliToyMCReconstruction::GetFittedTrackFromSeed(const AliToyMCTrack *tr, const AliExternalTrackParam *seed)
 809 {
 810   //
 811   //
 812   //
 813
 814   // create track
 815   AliExternalTrackParam *track = new AliExternalTrackParam(*seed);
 816
 817   Int_t ncls=(fClusterType == 0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
 818
 819   const AliTPCROC * roc = AliTPCROC::Instance();
 820
 821   const Double_t kRTPC0  = roc->GetPadRowRadii(0,0);
 822   const Double_t kRTPC1  = roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
 823   const Double_t kMaxSnp = 0.85;
 824   const Double_t kMaxR   = 500.;
 825   const Double_t kMaxZ   = 500.;
 826
 827   //   const Double_t kMaxZ0=220;
 828 //   const Double_t kZcut=3;
 829
 830   const Double_t refX = tr->GetX();
 831
 832   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 833
 834   // loop over all other points and add to the track
 835   for (Int_t ipoint=ncls-1; ipoint>=0; --ipoint){
 836     AliTrackPoint pIn;
 837     const AliTPCclusterMI *cl=tr->GetSpacePoint(ipoint);
 838     if (fClusterType == 1) cl=tr->GetDistortedSpacePoint(ipoint);
 839     SetTrackPointFromCluster(cl, pIn);
 840     if (fCorrectionType != kNoCorrection){
 841       Float_t xyz[3]={0,0,0};
 842       pIn.GetXYZ(xyz);
 843 //       if ( fCorrectionType == kIdeal ) xyz[2] = cl->GetZ();
 844       fTPCCorrection->CorrectPoint(xyz, cl->GetDetector());
 845       pIn.SetXYZ(xyz);
 846     }
 847     // rotate the cluster to the local detector frame
 848     track->Rotate(((cl->GetDetector()%18)*20+10)*TMath::DegToRad());
 849     AliTrackPoint prot = pIn.Rotate(track->GetAlpha());   // rotate to the local frame - non distoted  point
 850     if (TMath::Abs(prot.GetX())<kRTPC0) continue;
 851     if (TMath::Abs(prot.GetX())>kRTPC1) continue;
 852     //
 853     Bool_t res=AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
 854
 855     if (!res) break;
 856
 857     if (TMath::Abs(track->GetZ())>kMaxZ) break;
 858     if (TMath::Abs(track->GetX())>kMaxR) break;
 859 //     if (TMath::Abs(track->GetZ())<kZcut)continue;
 860     //
 861     Double_t pointPos[2]={0,0};
 862     Double_t pointCov[3]={0,0,0};
 863     pointPos[0]=prot.GetY();//local y
 864     pointPos[1]=prot.GetZ();//local z
 865     pointCov[0]=prot.GetCov()[3];//simay^2
 866     pointCov[1]=prot.GetCov()[4];//sigmayz
 867     pointCov[2]=prot.GetCov()[5];//sigmaz^2
 868
 869     if (!track->Update(pointPos,pointCov)) {printf("no update\n"); break;}
 870   }
 871
 872   AliTrackerBase::PropagateTrackTo2(track,refX,kMass,5.,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
 873
 874   // rotate fittet track to the frame of the original track and propagate to same reference
 875   track->Rotate(tr->GetAlpha());
 876
 877   AliTrackerBase::PropagateTrackTo2(track,refX,kMass,1.,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
 878
 879   return track;
 880 }
 881
 882
 883 //____________________________________________________________________________________
 884 AliExternalTrackParam* AliToyMCReconstruction::GetFittedTrackFromSeedAllClusters(const AliToyMCTrack *tr, const AliExternalTrackParam *seed, Int_t &nClus)
 885 {
 886   //
 887   // Tracking for given seed on an array of clusters
 888   //
 889
 890   // create track
 891   AliExternalTrackParam *track = new AliExternalTrackParam(*seed);
 892
 893   const AliTPCROC * roc = AliTPCROC::Instance();
 894
 895   const Double_t kRTPC0    = roc->GetPadRowRadii(0,0);
 896   const Double_t kRTPC1    = roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
 897   const Int_t kNRowsTPC    = roc->GetNRows(0) + roc->GetNRows(36) - 1;
 898   const Int_t kIRowsTPC    = roc->GetNRows(0) - 1;
 899   const Double_t kMaxSnp   = 0.85;
 900   const Double_t kMaxR     = 500.;
 901   const Double_t kMaxZ     = 500.;
 902   const Double_t roady     = 100.;
 903   const Double_t roadz     = 100.;
 904
 905   const Double_t refX = tr->GetX();
 906
 907   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 908
 909   Int_t  secCur   = -1;
 910   UInt_t indexCur = 0;
 911   Double_t xCur, yCur, zCur = 0.;
 912
 913   Float_t vDrift = GetVDrift();
 914
 915   // first propagate seed to outermost row
 916   AliTrackerBase::PropagateTrackTo(track,kRTPC1,kMass,5,kFALSE,kMaxSnp);
 917
 918   // Loop over rows and find the cluster candidates
 919   for( Int_t iRow = kNRowsTPC; iRow >= 0; --iRow ){
 920
 921     // inner or outer sector
 922     Bool_t bInnerSector = kTRUE;
 923     if(iRow > kIRowsTPC) bInnerSector = kFALSE;
 924
 925     // nearest track point/cluster (to be found)
 926     AliTrackPoint nearestPoint;
 927     AliTPCclusterMI *nearestCluster = NULL;
 928
 929     // Inner Sector
 930     if(bInnerSector){
 931
 932       // Propagate to center of pad row and extract parameters
 933       AliTrackerBase::PropagateTrackTo(track,roc->GetPadRowRadii(0,iRow),kMass,5,kFALSE,kMaxSnp);
 934       xCur   = track->GetX();
 935       yCur   = track->GetY();
 936       zCur   = track->GetZ();
 937       if ( !fIdealTracking ) {
 938         zCur = zCur/vDrift + fTime0; // Look at time, not at z!
 939       }
 940       secCur = GetSector(track);
 941
 942       // Find the nearest cluster (TODO: correct road settings!)
 943       Printf("inner tracking here: x = %.2f, y = %.2f, z = %.2f (Row %d Sector %d)",xCur,yCur,zCur,iRow,secCur);
 944       nearestCluster = fInnerSectorArray[secCur%fkNSectorInner][iRow].FindNearest2(yCur,zCur,roady,roadz,indexCur);
 945
 946       // Move to next row if now cluster found
 947       if(!nearestCluster) continue;
 948       //Printf("Nearest Clusters = %d (of %d) ",indexCur,fInnerSectorArray[secCur%fkNSectorInner][iRow].GetN());
 949
 950     }
 951
 952     // Outer sector
 953     else{
 954
 955       // Propagate to center of pad row and extract parameters
 956       AliTrackerBase::PropagateTrackTo(track,roc->GetPadRowRadii(36,iRow-kIRowsTPC-1),kMass,5,kFALSE,kMaxSnp);
 957       xCur   = track->GetX();
 958       yCur   = track->GetY();
 959       zCur   = track->GetZ();
 960       if ( !fIdealTracking ) {
 961         zCur = zCur/vDrift + fTime0; // Look at time, not at z!
 962       }
 963       secCur = GetSector(track);
 964
 965       // Find the nearest cluster (TODO: correct road settings!)
 966       Printf("outer tracking here: x = %.2f, y = %.2f, z = %.2f (Row %d Sector %d)",xCur,yCur,zCur,iRow,secCur);
 967       nearestCluster = fOuterSectorArray[(secCur-fkNSectorInner*2)%fkNSectorOuter][iRow-kIRowsTPC-1].FindNearest2(yCur,zCur,roady,roadz,indexCur);
 968
 969       // Move to next row if now cluster found
 970       if(!nearestCluster) continue;
 971       //Printf("Nearest Clusters = %d (of %d)",indexCur,fOuterSectorArray[(secCur-fkNSectorInner*2)%fkNSectorOuter][iRow-kIRowsTPC-1].GetN());
 972
 973     }
 974
 975     // create track point from cluster
 976     SetTrackPointFromCluster(nearestCluster,nearestPoint);
 977
 978     //Printf("Track point = %.2f %.2f %.2f",nearestPoint.GetX(),nearestPoint.GetY(),nearestPoint.GetZ());
 979
 980     // correction
 981     // TODO: also correction when looking for the next cluster?
 982     if (fCorrectionType != kNoCorrection){
 983       Float_t xyz[3]={0,0,0};
 984       nearestPoint.GetXYZ(xyz);
 985       fTPCCorrection->CorrectPoint(xyz, nearestCluster->GetDetector());
 986       nearestPoint.SetXYZ(xyz);
 987     }
 988
 989     // rotate the cluster to the local detector frame
 990     track->Rotate(((nearestCluster->GetDetector()%18)*20+10)*TMath::DegToRad());
 991     AliTrackPoint prot = nearestPoint.Rotate(track->GetAlpha());   // rotate to the local frame - non distoted  point
 992     if (TMath::Abs(prot.GetX())<kRTPC0) continue;
 993     if (TMath::Abs(prot.GetX())>kRTPC1) continue;
 994
 995
 996     //Printf("Rotated Track point = %.2f %.2f %.2f",prot.GetX(),prot.GetY(),prot.GetZ());
 997
 998     // update track with the nearest track point
 999     Bool_t res=AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
1000
1001     if (!res) break;
1002
1003     if (TMath::Abs(track->GetZ())>kMaxZ) break;
1004     if (TMath::Abs(track->GetX())>kMaxR) break;
1005     //if (TMath::Abs(track->GetZ())<kZcut)continue;
1006
1007       Double_t pointPos[2]={0,0};
1008       Double_t pointCov[3]={0,0,0};
1009       pointPos[0]=prot.GetY();//local y
1010       pointPos[1]=prot.GetZ();//local z
1011       pointCov[0]=prot.GetCov()[3];//simay^2
1012       pointCov[1]=prot.GetCov()[4];//sigmayz
1013       pointCov[2]=prot.GetCov()[5];//sigmaz^2
1014
1015       if (!track->Update(pointPos,pointCov)) {printf("no update\n");}
1016
1017       ++nClus;
1018   }
1019
1020
1021   // propagation to refX
1022   AliTrackerBase::PropagateTrackTo(track,refX,kMass,5.,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
1023
1024   // rotate fittet track to the frame of the original track and propagate to same reference
1025   track->Rotate(tr->GetAlpha());
1026
1027   AliTrackerBase::PropagateTrackTo(track,refX,kMass,1.,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
1028
1029   Printf("We have %d clusters in this track!",nClus);
1030
1031   return track;
1032 }
1033
1034 //____________________________________________________________________________________
1035 AliExternalTrackParam* AliToyMCReconstruction::ClusterToTrackAssociation(const AliTPCseed *seed, Int_t trackID, Int_t &nClus)
1036 {
1037   //
1038   // Cluster to track association for given seed on an array of clusters
1039   //
1040
1041   // create track
1042   AliExternalTrackParam *track = new AliExternalTrackParam(*seed);
1043
1044   const AliTPCROC * roc = AliTPCROC::Instance();
1045
1046   const Double_t kRTPC0    = roc->GetPadRowRadii(0,0);
1047   const Double_t kRTPC1    = roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
1048   const Int_t kNRowsTPC    = roc->GetNRows(0) + roc->GetNRows(36) - 1;
1049   const Int_t kIRowsTPC    = roc->GetNRows(0) - 1;
1050   const Double_t kMaxSnp   = 0.85;
1051   const Double_t kMaxR     = 500.;
1052   const Double_t kMaxZ     = 500.;
1053   const Double_t roady     = 0.1;
1054   const Double_t roadz     = 0.01;
1055
1056   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1057
1058   Int_t  secCur, secOld   = -1;
1059   UInt_t indexCur = 0;
1060   Double_t xCur, yCur, zCur = 0.;
1061
1062 //   Float_t vDrift = GetVDrift();
1063
1064   // first propagate seed to outermost row
1065   Bool_t res0=AliTrackerBase::PropagateTrackTo(track,kRTPC1,kMass,5,kFALSE,kMaxSnp);
1066
1067   // Loop over rows and find the cluster candidates
1068   for( Int_t iRow = kNRowsTPC; iRow >= 0; --iRow ){
1069
1070     // inner or outer sector
1071     Bool_t bInnerSector = kTRUE;
1072     if(iRow > kIRowsTPC) bInnerSector = kFALSE;
1073
1074     // nearest track point/cluster (to be found)
1075     AliTrackPoint nearestPoint;
1076     AliTPCclusterMI *nearestCluster = NULL;
1077
1078     // Inner Sector
1079     if(bInnerSector){
1080
1081       // Propagate to center of pad row and extract parameters
1082       AliTrackerBase::PropagateTrackTo(track,roc->GetPadRowRadii(0,iRow),kMass,5,kFALSE,kMaxSnp);
1083       xCur   = track->GetX();
1084       yCur   = track->GetY();
1085       zCur   = track->GetZ();
1086       secCur = GetSector(track);
1087
1088       // Find the nearest cluster (TODO: correct road settings!)
1089       //Printf("inner tracking here: x = %.2f, y = %.2f, z = %.6f (Row %d Sector %d)",xCur,yCur,zCur,iRow,secCur);
1090       nearestCluster = fInnerSectorArray[secCur%fkNSectorInner][iRow].FindNearest2(yCur,zCur,roady,roadz,indexCur);
1091
1092       // Look again at -y if nothing found here and the sector changed with respect to the last found cluster
1093       // Increase also the road in this case
1094       if(!nearestCluster && secCur != secOld && secOld > -1){
1095         //Printf("inner tracking here 2: x = %.2f, y = %.2f, z = %.6f (Row %d Sector %d)",xCur,-yCur,zCur,iRow,secCur);
1096         nearestCluster = fInnerSectorArray[secCur%fkNSectorInner][iRow].FindNearest2(-yCur,zCur,roady*100,roadz,indexCur);
1097       }
1098
1099       // Move to next row if now cluster found
1100       if(!nearestCluster) continue;
1101       //Printf("Nearest Clusters = %d (of %d) ",indexCur,fInnerSectorArray[secCur%fkNSectorInner][iRow].GetN());
1102
1103     }
1104
1105     // Outer sector
1106     else{
1107
1108       // Propagate to center of pad row and extract parameters
1109       AliTrackerBase::PropagateTrackTo(track,roc->GetPadRowRadii(36,iRow-kIRowsTPC-1),kMass,5,kFALSE,kMaxSnp);
1110       xCur   = track->GetX();
1111       yCur   = track->GetY();
1112       zCur   = track->GetZ();
1113       secCur = GetSector(track);
1114
1115       // Find the nearest cluster (TODO: correct road settings!)
1116       Printf("res0 = %d, outer tracking here: x = %.2f, y = %.2f, z = %.6f (Row %d Sector %d)",res0,xCur,yCur,zCur,iRow,secCur);
1117       nearestCluster = fOuterSectorArray[(secCur-fkNSectorInner*2)%fkNSectorOuter][iRow-kIRowsTPC-1].FindNearest2(yCur,zCur,roady,roadz,indexCur);
1118
1119       // Look again at -y if nothing found here and the sector changed with respect to the last found cluster
1120       // Increase also the road in this case
1121       if(!nearestCluster && secCur != secOld && secOld > -1){
1122         Printf("outer tracking here 2: x = %.2f, y = %.2f, z = %.6f (Row %d Sector %d)",xCur,-yCur,zCur,iRow,secCur);
1123         nearestCluster = fOuterSectorArray[(secCur-fkNSectorInner*2)%fkNSectorOuter][iRow-kIRowsTPC-1].FindNearest2(-yCur,zCur,roady*100,roadz,indexCur);
1124       }
1125
1126
1127       // Move to next row if now cluster found
1128       if(!nearestCluster) continue;
1129       Printf("Nearest Clusters = %d (of %d)",indexCur,fOuterSectorArray[(secCur-fkNSectorInner*2)%fkNSectorOuter][iRow-kIRowsTPC-1].GetN());
1130
1131     }
1132
1133     // create track point from cluster
1134     SetTrackPointFromCluster(nearestCluster,nearestPoint);
1135
1136     //Printf("Track point = %.2f %.2f %.2f",nearestPoint.GetX(),nearestPoint.GetY(),nearestPoint.GetZ());
1137
1138     // rotate the cluster to the local detector frame
1139     track->Rotate(((nearestCluster->GetDetector()%18)*20+10)*TMath::DegToRad());
1140     AliTrackPoint prot = nearestPoint.Rotate(track->GetAlpha());   // rotate to the local frame - non distoted  point
1141     if (TMath::Abs(prot.GetX())<kRTPC0) continue;
1142     if (TMath::Abs(prot.GetX())>kRTPC1) continue;
1143
1144
1145     //Printf("Rotated Track point = %.2f %.2f %.2f",prot.GetX(),prot.GetY(),prot.GetZ());
1146
1147     // update track with the nearest track point
1148     Bool_t res=AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
1149
1150     if (!res) break;
1151
1152     if (TMath::Abs(track->GetZ())>kMaxZ) break;
1153     if (TMath::Abs(track->GetX())>kMaxR) break;
1154     //if (TMath::Abs(track->GetZ())<kZcut)continue;
1155
1156       Double_t pointPos[2]={0,0};
1157       Double_t pointCov[3]={0,0,0};
1158       pointPos[0]=prot.GetY();//local y
1159       pointPos[1]=prot.GetZ();//local z
1160       pointCov[0]=prot.GetCov()[3];//simay^2
1161       pointCov[1]=prot.GetCov()[4];//sigmayz
1162       pointCov[2]=prot.GetCov()[5];//sigmaz^2
1163
1164       if (!track->Update(pointPos,pointCov)) {printf("no update\n");}
1165       secOld = secCur;
1166
1167       //Printf("Cluster belongs to track = %d",nearestCluster->GetLabel(0));
1168
1169       // only count as associate cluster if it belongs to correct track!
1170       if(nearestCluster->GetLabel(0) == trackID)
1171         ++nClus;
1172   }
1173
1174   Printf("We have %d clusters in this track!",nClus);
1175
1176   return track;
1177 }
1178
1179
1180 //____________________________________________________________________________________
1181 void AliToyMCReconstruction::InitSpaceCharge()
1182 {
1183   //
1184   // Init the space charge map
1185   //
1186
1187   TString filename="$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps5_50kHz_precal.root";
1188   if (fTree) {
1189     TList *l=fTree->GetUserInfo();
1190     for (Int_t i=0; i<l->GetEntries(); ++i) {
1191       TString s(l->At(i)->GetName());
1192       if (s.Contains("SC_")) {
1193         filename=s;
1194         break;
1195       }
1196     }
1197   }
1198
1199   printf("Initialising the space charge map using the file: '%s'\n",filename.Data());
1200   TFile f(filename.Data());
1201   fTPCCorrection=(AliTPCCorrection*)f.Get("map");
1202
1203   //   fTPCCorrection = new AliTPCSpaceCharge3D();
1204   //   fTPCCorrection->SetSCDataFileName("$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps10_50kHz.root");
1205   //   fTPCCorrection->SetOmegaTauT1T2(0.325,1,1); // Ne CO2
1206   // //   fTPCCorrection->SetOmegaTauT1T2(0.41,1,1.05); // Ar CO2
1207   //   fTPCCorrection->InitSpaceCharge3DDistortion();
1208
1209 }
1210
1211 //____________________________________________________________________________________
1212 Double_t AliToyMCReconstruction::GetVDrift() const
1213 {
1214   //
1215   //
1216   //
1217   return fTPCParam->GetDriftV();
1218 }
1219
1220 //____________________________________________________________________________________
1221 Double_t AliToyMCReconstruction::GetZLength(Int_t roc) const
1222 {
1223   //
1224   //
1225   //
1226   if (roc<0 || roc>71) return -1;
1227   return fTPCParam->GetZLength(roc);
1228 }
1229
1230 //____________________________________________________________________________________
1231 TTree* AliToyMCReconstruction::ConnectTrees (const char* files) {
1232   TString s=gSystem->GetFromPipe(Form("ls %s",files));
1233
1234   TTree *tFirst=0x0;
1235   TObjArray *arrFiles=s.Tokenize("\n");
1236
1237   for (Int_t ifile=0; ifile<arrFiles->GetEntriesFast(); ++ifile){
1238     TString name(arrFiles->At(ifile)->GetName());
1239
1240     TPRegexp reg(".*([0-9]_[0-9]_[0-9]_[0-9]{3}_[0-9]{2}).debug.root");
1241     TObjArray *arrMatch=0x0;
1242     arrMatch=reg.MatchS(name);
1243     TString matchName;
1244     if (arrMatch && arrMatch->At(1)) matchName=arrMatch->At(1)->GetName();
1245     else matchName=Form("%02d",ifile);
1246     delete arrMatch;
1247
1248     if (!tFirst) {
1249       TFile *f=TFile::Open(name.Data());
1250       if (!f) continue;
1251       TTree *t=(TTree*)f->Get("Tracks");
1252       if (!t) {
1253         delete f;
1254         continue;
1255       }
1256
1257       t->SetName(matchName.Data());
1258       tFirst=t;
1259     } else {
1260       tFirst->AddFriend(Form("t%s=Tracks",matchName.Data()), name.Data());
1261 //       tFirst->AddFriend(Form("t%d=Tracks",ifile), name.Data());
1262     }
1263   }
1264
1265   tFirst->GetListOfFriends()->Print();
1266   return tFirst;
1267 }
1268
1269 //____________________________________________________________________________________
1270 Int_t AliToyMCReconstruction::LoadOuterSectors() {
1271   //-----------------------------------------------------------------
1272   // This function fills outer TPC sectors with clusters.
1273   // Copy and paste from AliTPCtracker
1274   //-----------------------------------------------------------------
1275   Int_t nrows = fOuterSectorArray->GetNRows();
1276   UInt_t index=0;
1277   for (Int_t sec = 0;sec<fkNSectorOuter;sec++)
1278     for (Int_t row = 0;row<nrows;row++){
1279       AliTPCtrackerRow*  tpcrow = &(fOuterSectorArray[sec%fkNSectorOuter][row]);
1280       Int_t sec2 = sec+2*fkNSectorInner;
1281       //left
1282       Int_t ncl = tpcrow->GetN1();
1283       while (ncl--) {
1284         AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1285         index=(((sec2<<8)+row)<<16)+ncl;
1286         tpcrow->InsertCluster(c,index);
1287       }
1288       //right
1289       ncl = tpcrow->GetN2();
1290       while (ncl--) {
1291         AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1292         index=((((sec2+fkNSectorOuter)<<8)+row)<<16)+ncl;
1293         tpcrow->InsertCluster(c,index);
1294       }
1295       //
1296       // write indexes for fast acces
1297       //
1298       for (Int_t i=0;i<510;i++)
1299         tpcrow->SetFastCluster(i,-1);
1300       for (Int_t i=0;i<tpcrow->GetN();i++){
1301         Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1302         tpcrow->SetFastCluster(zi,i);  // write index
1303       }
1304       Int_t last = 0;
1305       for (Int_t i=0;i<510;i++){
1306         if (tpcrow->GetFastCluster(i)<0)
1307           tpcrow->SetFastCluster(i,last);
1308         else
1309           last = tpcrow->GetFastCluster(i);
1310       }
1311     }
1312   return 0;
1313 }
1314
1315
1316 //____________________________________________________________________________________
1317 Int_t  AliToyMCReconstruction::LoadInnerSectors() {
1318   //-----------------------------------------------------------------
1319   // This function fills inner TPC sectors with clusters.
1320   // Copy and paste from AliTPCtracker
1321   //-----------------------------------------------------------------
1322   Int_t nrows = fInnerSectorArray->GetNRows();
1323   UInt_t index=0;
1324   for (Int_t sec = 0;sec<fkNSectorInner;sec++)
1325     for (Int_t row = 0;row<nrows;row++){
1326       AliTPCtrackerRow*  tpcrow = &(fInnerSectorArray[sec%fkNSectorInner][row]);
1327       //
1328       //left
1329       Int_t ncl = tpcrow->GetN1();
1330       while (ncl--) {
1331         AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1332         index=(((sec<<8)+row)<<16)+ncl;
1333         tpcrow->InsertCluster(c,index);
1334       }
1335       //right
1336       ncl = tpcrow->GetN2();
1337       while (ncl--) {
1338         AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1339         index=((((sec+fkNSectorInner)<<8)+row)<<16)+ncl;
1340         tpcrow->InsertCluster(c,index);
1341       }
1342       //
1343       // write indexes for fast acces
1344       //
1345       for (Int_t i=0;i<510;i++)
1346         tpcrow->SetFastCluster(i,-1);
1347       for (Int_t i=0;i<tpcrow->GetN();i++){
1348         Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1349         tpcrow->SetFastCluster(zi,i);  // write index
1350       }
1351       Int_t last = 0;
1352       for (Int_t i=0;i<510;i++){
1353         if (tpcrow->GetFastCluster(i)<0)
1354           tpcrow->SetFastCluster(i,last);
1355         else
1356           last = tpcrow->GetFastCluster(i);
1357       }
1358
1359     }
1360   return 0;
1361 }
1362
1363 //____________________________________________________________________________________
1364 Int_t  AliToyMCReconstruction::GetSector(AliExternalTrackParam *track) {
1365   //-----------------------------------------------------------------
1366   // This function returns the sector number for a given track
1367   //-----------------------------------------------------------------
1368
1369   Int_t sector = -1;
1370
1371   // get the sector number
1372   // rotate point to global system (track is already global!)
1373   Double_t xd[3];
1374   track->GetXYZ(xd);
1375   //track->Local2GlobalPosition(xd,track->GetAlpha());
1376
1377   // use TPCParams to get the sector number
1378   Float_t xyz[3] = {xd[0],xd[1],xd[2]};
1379   Int_t   i[3]   = {0,0,0};
1380   if(fTPCParam){
1381     sector  = fTPCParam->Transform0to1(xyz,i);
1382   }
1383
1384   return sector;
1385 }
1386
1387 //____________________________________________________________________________________
1388 void  AliToyMCReconstruction::FillSectorStructure(Int_t maxev) {
1389   //-----------------------------------------------------------------
1390   // This function fills the sector structure of AliToyMCReconstruction
1391   //-----------------------------------------------------------------
1392
1393   // cluster array of all sectors
1394   fInnerSectorArray = new AliTPCtrackerSector[fkNSectorInner];
1395   fOuterSectorArray = new AliTPCtrackerSector[fkNSectorOuter];
1396
1397   for (Int_t i=0; i<fkNSectorInner; ++i) fInnerSectorArray[i].Setup(fTPCParam,0);
1398   for (Int_t i=0; i<fkNSectorOuter; ++i) fOuterSectorArray[i].Setup(fTPCParam,1);
1399
1400   Int_t count[72][96] = { {0} , {0} };
1401
1402   for (Int_t iev=0; iev<maxev; ++iev){
1403     printf("==============  Fill Clusters: Processing Event %6d  =================\n",iev);
1404     fTree->GetEvent(iev);
1405     for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
1406       printf(" > ======  Fill Clusters: Processing Track %6d ========  \n",itr);
1407       const AliToyMCTrack *tr=fEvent->GetTrack(itr);
1408
1409       // number of clusters to loop over
1410       const Int_t ncls=(fClusterType==0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
1411
1412       for(Int_t icl=0; icl<ncls; ++icl){
1413
1414         AliTPCclusterMI *cl=const_cast<AliTPCclusterMI *>(tr->GetSpacePoint(icl));
1415         if (fClusterType==1) cl=const_cast<AliTPCclusterMI *>(tr->GetDistortedSpacePoint(icl));
1416         if (!cl) continue;
1417
1418         Int_t sec = cl->GetDetector();
1419         Int_t row = cl->GetRow();
1420
1421         // set Q of the cluster to 1, Q==0 does not work for the seeding
1422         cl->SetQ(1);
1423
1424         // set cluster time to cluster Z (if not ideal tracking)
1425         if ( !fIdealTracking ) {
1426           // a 'valid' position in z is needed for the seeding procedure
1427 //           cl->SetZ(cl->GetTimeBin()*GetVDrift());
1428           cl->SetZ(cl->GetTimeBin());
1429         }
1430         //Printf("Fill clusters (sector %d row %d): %.2f %.2f %.2f %.2f",sec,row,cl->GetX(),cl->GetY(),cl->GetZ(),cl->GetTimeBin());
1431
1432         // fill arrays for inner and outer sectors (A/C side handled internally)
1433         if (sec<fkNSectorInner*2){
1434           fInnerSectorArray[sec%fkNSectorInner].InsertCluster(cl, count[sec][row], fTPCParam);
1435         }
1436         else{
1437           fOuterSectorArray[(sec-fkNSectorInner*2)%fkNSectorOuter].InsertCluster(cl, count[sec][row], fTPCParam);
1438         }
1439
1440         ++count[sec][row];
1441       }
1442     }
1443   }
1444
1445   // fill the arrays completely
1446   // LoadOuterSectors();
1447   // LoadInnerSectors();
1448
1449   // // check the arrays
1450   // for (Int_t i=0; i<fkNSectorInner; ++i){
1451   //   for (Int_t j=0; j<fInnerSectorArray[i].GetNRows(); ++j){
1452   //     if(fInnerSectorArray[i][j].GetN()>0){
1453   //    Printf("Inner: Sector %d Row %d : %d",i,j,fInnerSectorArray[i][j].GetN());
1454   //     }
1455   //   }
1456   // }
1457   // for (Int_t i=0; i<fkNSectorInner; ++i){
1458   //   for (Int_t j=0; j<fOuterSectorArray[i].GetNRows(); ++j){
1459   //     if(fOuterSectorArray[i][j].GetN()>0){
1460   //    Printf("Outer: Sector %d Row %d : %d",i,j,fOuterSectorArray[i][j].GetN());
1461   //     }
1462   //   }
1463   // }
1464 }
1465
1466 //____________________________________________________________________________________
1467 void  AliToyMCReconstruction::FillSectorStructureAC(Int_t maxev) {
1468   //-----------------------------------------------------------------
1469   // This function fills the sector structure of AliToyMCReconstruction
1470   //-----------------------------------------------------------------
1471
1472   /*
1473    my god is the AliTPCtrackerSector stuff complicated!!!
1474    Ok, so here we will not fill the fClusters1 and fClusters2 of AliTPCtrackerRow,
1475    using InsertCluster of AliTPCtrackerSector, but only the fClusters via InsertCluster
1476    of AliTPCtrackerRow itself which then will not be owner, but we create an array in
1477    here (fAllClusters) which owns all clusters ...
1478   */
1479
1480   fIsAC=kTRUE;
1481   // cluster array of all sectors
1482   fInnerSectorArray = new AliTPCtrackerSector[2*fkNSectorInner];
1483   fOuterSectorArray = new AliTPCtrackerSector[2*fkNSectorOuter];
1484
1485   for (Int_t i=0; i<2*fkNSectorInner; ++i) {
1486     fInnerSectorArray[i].Setup(fTPCParam,0);
1487   }
1488
1489   for (Int_t i=0; i<2*fkNSectorOuter; ++i) {
1490     fOuterSectorArray[i].Setup(fTPCParam,1);
1491   }
1492
1493   Int_t count[72][96] = { {0} , {0} };
1494
1495   for (Int_t iev=0; iev<maxev; ++iev){
1496     printf("==============  Fill Clusters: Processing Event %6d  =================\n",iev);
1497     fTree->GetEvent(iev);
1498     for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
1499       printf(" > ======  Fill Clusters: Processing Track %6d ========  \n",itr);
1500       const AliToyMCTrack *tr=fEvent->GetTrack(itr);
1501
1502       // number of clusters to loop over
1503       const Int_t ncls=(fClusterType==0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
1504
1505       // check if expansion of the cluster arrays is needed.
1506       if (fAllClusters.GetEntriesFast()+ncls>=fAllClusters.Capacity()) fAllClusters.Expand(2*fAllClusters.Capacity());
1507       for(Int_t icl=0; icl<ncls; ++icl){
1508
1509         AliTPCclusterMI *cl=const_cast<AliTPCclusterMI *>(tr->GetSpacePoint(icl));
1510         if (fClusterType==1) cl=const_cast<AliTPCclusterMI *>(tr->GetDistortedSpacePoint(icl));
1511         if (!cl) continue;
1512
1513         // register copy to the cluster array
1514         cl = new(fAllClusters[fAllClusters.GetEntriesFast()]) AliTPCclusterMI(*cl);
1515
1516         Int_t sec = cl->GetDetector();
1517         Int_t row = cl->GetRow();
1518
1519         // set Q of the cluster to 1, Q==0 does not work for the seeding
1520         cl->SetQ(1);
1521
1522         // set cluster time to cluster Z (if not ideal tracking)
1523         if ( !fIdealTracking ) {
1524           // a 'valid' position in z is needed for the seeding procedure
1525           cl->SetZ(cl->GetTimeBin()*GetVDrift());
1526 //           cl->SetZ(cl->GetTimeBin());
1527         }
1528         //Printf("Fill clusters (sector %d row %d): %.2f %.2f %.2f %.2f",sec,row,cl->GetX(),cl->GetY(),cl->GetZ(),cl->GetTimeBin());
1529
1530         // fill arrays for inner and outer sectors (A/C side handled internally)
1531         if (sec<fkNSectorInner*2){
1532           fInnerSectorArray[sec][row].InsertCluster(cl, count[sec][row]);
1533         }
1534         else{
1535           fOuterSectorArray[(sec-fkNSectorInner*2)][row].InsertCluster(cl, count[sec][row]);
1536         }
1537
1538         ++count[sec][row];
1539       }
1540     }
1541   }
1542
1543 }
1544
1545 //____________________________________________________________________________________
1546 AliToyMCTrack *AliToyMCReconstruction::ConvertTPCSeedToToyMCTrack(const AliTPCseed &seed)
1547 {
1548   //
1549   //
1550   //
1551
1552
1553   AliToyMCTrack *tToy = new AliToyMCTrack(seed);
1554
1555   for (Int_t icl=0; icl<159; ++icl){
1556     const AliTPCclusterMI *cl=seed.GetClusterFast(icl);
1557     if (!cl) continue;
1558     if (fClusterType==0){
1559       tToy->AddSpacePoint(*cl);
1560     } else {
1561       tToy->AddDistortedSpacePoint(*cl);
1562     }
1563   }
1564
1565   return tToy;
1566 }
1567
1568 //____________________________________________________________________________________
1569 AliExternalTrackParam* AliToyMCReconstruction::GetRefittedTrack(const AliTPCseed &seed)
1570 {
1571   //
1572   //
1573   //
1574
1575   AliExternalTrackParam *track=0x0;
1576
1577   const Float_t vDrift=GetVDrift();
1578   const Float_t zLength=GetZLength(0);
1579   const Double_t kMaxSnp = 0.85;
1580   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1581
1582   // for propagation use a copy
1583   AliExternalTrackParam t0seed(seed);
1584   AliTrackerBase::PropagateTrackTo(&t0seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
1585   fTime0 = t0seed.GetZ()-zLength/vDrift;
1586   fCreateT0seed = kFALSE;
1587
1588   AliToyMCTrack *toyTrack = ConvertTPCSeedToToyMCTrack(seed);
1589
1590   fTime0 = t0seed.GetZ()-zLength/vDrift;
1591   fCreateT0seed = kFALSE;
1592   //         printf("seed (%.2g)\n",fTime0);
1593   AliExternalTrackParam *dummy = GetSeedFromTrack(toyTrack);
1594   if (dummy) {
1595     track = GetFittedTrackFromSeed(toyTrack, dummy);
1596     delete dummy;
1597     // propagate seed to 0
1598     AliTrackerBase::PropagateTrackTo(track,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
1599
1600   }
1601
1602   return track;
1603 }
1604
1605 //____________________________________________________________________________________
1606 AliTPCclusterMI* AliToyMCReconstruction::FindMiddleCluster(const AliTPCclusterMI *clInner, const AliTPCclusterMI *clOuter,
1607                                                            Double_t roady, Double_t roadz,
1608                                                            AliRieman &seedFit)
1609 {
1610   //
1611   //
1612   //
1613
1614   const Int_t rocInner = clInner->GetDetector();
1615   const Int_t rocOuter = clOuter->GetDetector();
1616
1617   if ( (rocInner%18) != (rocOuter%18) ){
1618     AliError("Currently only same Sector implemented");
1619     return 0x0;
1620   }
1621
1622   const Int_t innerDet=clInner->GetDetector();
1623   const Int_t outerDet=clOuter->GetDetector();
1624   const Int_t globalRowInner  = clInner->GetRow() +(innerDet >35)*63;
1625   const Int_t globalRowOuter  = clOuter->GetRow() +(outerDet >35)*63;
1626
1627   Int_t iMiddle  = (globalRowInner+globalRowOuter)/2;
1628   Int_t roc = innerDet;
1629   if (iMiddle>62){
1630     iMiddle-=63;
1631     roc=outerDet;
1632   }
1633
1634   const AliTPCtrackerRow& krMiddle = fOuterSectorArray[roc%36][iMiddle]; // middle
1635   // initial guess use simple linear interpolation
1636   Double_t y=(clInner->GetY()+clOuter->GetY())/2;
1637   Double_t z=(clInner->GetZ()+clOuter->GetZ())/2;
1638   if (seedFit.IsValid()) {
1639     // update values once the fit was performed
1640     y=seedFit.GetYat(krMiddle.GetX());
1641     z=seedFit.GetZat(krMiddle.GetX());
1642   }
1643
1644   AliTPCclusterMI *n=krMiddle.FindNearest(y,z,roady,roadz);
1645 //   if (n)
1646 //     printf("      Nearest cluster (%.2f, %.2f, %.2f) = m(%.2f, %.2f, %.2f : %d) i(%.2f, %.2f , %.2f : %d) o(%.2f, %.2f, %.2f : %d)\n",krMiddle.GetX(),y,z,
1647 //          n->GetX(), n->GetY(),n->GetZ(),n->GetLabel(0),
1648 //            clInner->GetX(), clInner->GetY(),clInner->GetZ(),clInner->GetLabel(0),
1649 //            clOuter->GetX(), clOuter->GetY(),clOuter->GetZ(),clOuter->GetLabel(0)
1650 //         );
1651   return n;
1652 }
1653
1654 //____________________________________________________________________________________
1655 void AliToyMCReconstruction::AddMiddleClusters(AliTPCseed *seed,
1656                                                const AliTPCclusterMI *clInner, const AliTPCclusterMI *clOuter,
1657                                                Double_t roady, Double_t roadz,
1658                                                Int_t &nTotalClusters, AliRieman &seedFit)
1659 {
1660   //
1661   // Iteratively add the middle clusters
1662   //
1663
1664   // update rieman fit with every second point added
1665   AliTPCclusterMI *clMiddle=FindMiddleCluster(clInner,clOuter,roady,roadz,seedFit);
1666
1667   // break iterative process
1668   if (!clMiddle) return;
1669
1670   const Int_t globalRowInner  = clInner->GetRow() +(clInner->GetDetector() >35)*63;
1671   const Int_t globalRowMiddle = clMiddle->GetRow()+(clMiddle->GetDetector()>35)*63;
1672   const Int_t globalRowOuter  = clOuter->GetRow() +(clOuter->GetDetector() >35)*63;
1673
1674   seed->SetClusterPointer(globalRowMiddle,clMiddle);
1675   ++nTotalClusters;
1676 //   printf("    > Total clusters: %d\n",nTotalClusters);
1677   seedFit.AddPoint(clMiddle->GetX(), clMiddle->GetY(), clMiddle->GetZ(),
1678                    TMath::Sqrt(clMiddle->GetSigmaY2()), TMath::Sqrt(clMiddle->GetSigmaZ2()));
1679
1680   if (!(seedFit.GetN()%3)) {
1681 //     printf("      call update: %d (%d)\n",seedFit.GetN(),nTotalClusters);
1682 //     printf("      Riemann results: valid=%d, Chi2=%.2f, Chi2Y=%.2f, Chi2Z=%.2f\n",
1683 //            seedFit.IsValid(), seedFit.GetChi2(), seedFit.GetChi2Y(), seedFit.GetChi2Z());
1684     seedFit.Update();
1685   }
1686   if ( seedFit.IsValid() && seedFit.GetChi2()>1000 ) return;
1687
1688   // Add middle clusters towards outer radius
1689   if (globalRowMiddle+1<globalRowOuter) AddMiddleClusters(seed,clMiddle,clOuter,roady,roadz,nTotalClusters,seedFit);
1690
1691   // Add middle clusters towards innter radius
1692   if (globalRowInner+1<globalRowMiddle) AddMiddleClusters(seed,clInner,clMiddle,roady,roadz,nTotalClusters,seedFit);
1693
1694   return;
1695 }
1696
1697 //____________________________________________________________________________________
1698 void AliToyMCReconstruction::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t iRowInner, Int_t iRowOuter)
1699 {
1700   //
1701   // find seeds in a sector, requires iRowInner < iRowOuter
1702   // iRowXXX runs from 0-159, so over IROC and OROC
1703   //
1704
1705   if (!fIsAC) {
1706     AliError("This function requires the sector arrays filled for AC tracking");
1707     return;
1708   }
1709
1710   // swap rows in case they are in the wrong order
1711   if (iRowInner>iRowOuter) {
1712     Int_t tmp=iRowInner;
1713     iRowInner=iRowOuter;
1714     iRowOuter=tmp;
1715   }
1716
1717   // only for CookLabel
1718   AliTPCtracker tpcTracker(fTPCParam);
1719
1720   // Define the search roads:
1721   // timeRoadCombinatorics - the maximum time difference used for the
1722   //    combinatorics. Since we don't have a 'z-Vertex' estimate this will
1723   //    reduce the combinatorics significantly when iterating on one TF
1724   //    use a little more than one full drift length of the TPC to allow for
1725   //    distortions
1726   //
1727   // timeRoad - the time difference allowed when associating the cluster
1728   //    in the middle of the other two used for the initial search
1729   //
1730   // padRoad  - the local y difference allowed when associating the middle cluster
1731   Float_t vDrift=GetVDrift();
1732   Float_t zLength=GetZLength(0);
1733
1734 //   Double_t timeRoadCombinatorics = 270./vDrift;
1735 //   Double_t timeRoad = 20./vDrift;
1736   Double_t timeRoadCombinatorics = 270.;
1737   Double_t timeRoad = 20.;
1738   Double_t  padRoad = 10.;
1739
1740
1741   // fOuterSectorArray runs from 0-95, so from iRowXXX which runs from 0-159
1742   //   the number of rows in the IROC has to be subtracted
1743   const Int_t innerRows=fInnerSectorArray->GetNRows();
1744   const AliTPCtrackerRow& krOuter  = fOuterSectorArray[sec][iRowOuter   - innerRows];   // down
1745   const AliTPCtrackerRow& krInner  = fOuterSectorArray[sec][iRowInner   - innerRows];   // up
1746
1747   AliTPCseed *seed = new AliTPCseed;
1748
1749   const Int_t nMaxClusters=iRowOuter-iRowInner+1;
1750 //   Int_t nScannedClusters = 0;
1751
1752   // loop over all points in the firstand last search row
1753   for (Int_t iOuter=0; iOuter < krOuter; iOuter++) {
1754     const AliTPCclusterMI *clOuter = krOuter[iOuter];
1755     for (Int_t iInner=0; iInner < krInner; iInner++) {
1756       const AliTPCclusterMI *clInner = krInner[iInner];
1757 // printf("\n\n Check combination %d (%d), %d (%d)\n",iOuter, iInner, clOuter->GetLabel(0), clInner->GetLabel(0));
1758       // check maximum distance for combinatorics
1759       if (TMath::Abs(clOuter->GetZ()-clInner->GetZ())>timeRoadCombinatorics) continue;
1760 // printf("  Is inside one drift\n");
1761
1762       // use rieman fit for seed description
1763       AliRieman seedFit(159);
1764       // add first two points
1765       seedFit.AddPoint(clInner->GetX(), clInner->GetY(), clInner->GetZ(),
1766                        TMath::Sqrt(clInner->GetSigmaY2()), TMath::Sqrt(clInner->GetSigmaZ2()));
1767       seedFit.AddPoint(clOuter->GetX(), clOuter->GetY(), clOuter->GetZ(),
1768                        TMath::Sqrt(clOuter->GetSigmaY2()), TMath::Sqrt(clOuter->GetSigmaZ2()));
1769
1770       // Iteratively add all clusters in the respective middle
1771       Int_t nFoundClusters=2;
1772       AddMiddleClusters(seed,clInner,clOuter,padRoad,timeRoad,nFoundClusters,seedFit);
1773 //       printf("  Clusters attached: %d\n",nFoundClusters);
1774       seedFit.Update();
1775 //       printf("  Riemann results: valid=%d, Chi2=%.2f, Chi2Y=%.2f, Chi2Z=%.2f\n",
1776 //              seedFit.IsValid(), seedFit.GetChi2(), seedFit.GetChi2Y(), seedFit.GetChi2Z());
1777
1778       // check for minimum number of assigned clusters and a decent chi2
1779       if ( nFoundClusters<0.5*nMaxClusters || seedFit.GetChi2()>1000 ){
1780         seed->Reset();
1781         continue;
1782       }
1783 //       printf("  >>> Seed accepted\n");
1784       // add original outer clusters
1785       const Int_t globalRowInner  = clInner->GetRow() +(clInner->GetDetector() >35)*63;
1786       const Int_t globalRowOuter  = clOuter->GetRow() +(clOuter->GetDetector() >35)*63;
1787
1788       seed->SetClusterPointer(globalRowInner,const_cast<AliTPCclusterMI*>(clInner));
1789       seed->SetClusterPointer(globalRowOuter,const_cast<AliTPCclusterMI*>(clOuter));
1790
1791       // set parameters retrieved from AliRieman
1792       Double_t params[5]={0};
1793       Double_t covar[15]={0};
1794       const Double_t alpha=TMath::DegToRad()*(clInner->GetDetector()%18*20.+10.);
1795       const Double_t x=clInner->GetX();
1796       seedFit.GetExternalParameters(x,params,covar);
1797
1798       seed->Set(x,alpha,params,covar);
1799
1800       // set label of the seed. At least 60% of the clusters need the correct label
1801       tpcTracker.CookLabel(seed,.6);
1802
1803       arr->Add(seed);
1804       seed=new AliTPCseed;
1805     }
1806   }
1807   //delete surplus seed
1808   delete seed;
1809   seed=0x0;
1810
1811   // for debugging
1812   if (fStreamer && fTree) {
1813     //loop over all events and tracks and try to associate the seed to the track
1814
1815     const Double_t kMaxSnp = 0.85;
1816     const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1817
1818     for (Int_t iseed=0; iseed<arr->GetEntriesFast(); ++iseed){
1819       seed = static_cast<AliTPCseed*>(arr->UncheckedAt(iseed));
1820       AliExternalTrackParam extSeed(*seed);
1821       AliTrackerBase::PropagateTrackTo(&extSeed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
1822       Int_t seedLabel=seed->GetLabel();
1823
1824       // get original track
1825       Int_t iev=fMapTrackEvent.GetValue(TMath::Abs(seedLabel));
1826       Int_t itr=fMapTrackTrackInEvent.GetValue(TMath::Abs(seedLabel));
1827
1828       fTree->GetEvent(iev);
1829       const AliToyMCTrack *toyTrack = fEvent->GetTrack(itr);
1830
1831       AliExternalTrackParam extTrack(*toyTrack);
1832
1833       //propagate to same reference frame
1834       extSeed.Rotate(extTrack.GetAlpha());
1835       AliTrackerBase::PropagateTrackTo(&extSeed,0,kMass,1,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
1836
1837       fTime0=extSeed.GetZ()/vDrift;
1838
1839       Float_t z0=fEvent->GetZ();
1840       Float_t t0=fEvent->GetT0();
1841
1842       Int_t ctype(fCorrectionType);
1843
1844       (*fStreamer) << "Seeds" <<
1845       "iev="             << iev             <<
1846       "z0="              << z0              <<
1847       "t0="              << t0              <<
1848       "fTime0="          << fTime0          <<
1849       "itr="             << itr             <<
1850       "clsType="         << fClusterType    <<
1851       "corrType="        << ctype           <<
1852       "seedRow="         << fSeedingRow     <<
1853       "seedDist="        << fSeedingDist    <<
1854       "vDrift="          << vDrift          <<
1855       "zLength="         << zLength         <<
1856       "track.="          << &extSeed        <<
1857       "tOrig.="          << &extTrack       <<
1858       "seedLabel="       << seedLabel       <<
1859       "\n";
1860
1861     }
1862   }
1863 }
1864 //____________________________________________________________________________________
1865 void AliToyMCReconstruction::MakeSeeds(TObjArray * /*arr*/, Int_t sec, Int_t iRow1, Int_t iRow2)
1866 {
1867   //
1868   // Create seeds between i1 and i2 (stored in arr) with TLinearFitter
1869   //
1870   // sec: sector number
1871   // iRow1:  upper row
1872   // iRow2:  lower row
1873   //
1874
1875   // Create Fitter
1876   static TLinearFitter fitter(3,"pol2");
1877
1878   // Get 3 padrows (iRow1,iMiddle=(iRow1+iRow2)/2,iRow2)
1879   Int_t iMiddle = (iRow1+iRow2)/2;
1880   const AliTPCtrackerRow& krd = fOuterSectorArray[sec][iRow2-fInnerSectorArray->GetNRows()];   // down
1881   const AliTPCtrackerRow& krm = fOuterSectorArray[sec][iMiddle-fInnerSectorArray->GetNRows()]; // middle
1882   const AliTPCtrackerRow& kru = fOuterSectorArray[sec][iRow1-fInnerSectorArray->GetNRows()];   // up
1883
1884   // Loop over 3 cluster possibilities
1885   for (Int_t iu=0; iu < kru; iu++) {
1886     for (Int_t im=0; im < krm; im++) {
1887       for (Int_t id=0; id < krd; id++) {
1888
1889         // clear all points
1890         fitter.ClearPoints();
1891
1892         // add all three points to fitter
1893         Double_t xy[2] = {kru[iu]->GetX(),kru[iu]->GetY()};
1894         Double_t z     = kru[iu]->GetZ();
1895         fitter.AddPoint(xy,z);
1896
1897         xy[0] = krm[im]->GetX();
1898         xy[1] = krm[im]->GetY();
1899         z     = krm[im]->GetZ();
1900         fitter.AddPoint(xy,z);
1901
1902         xy[0] = krd[id]->GetX();
1903         xy[1] = krd[id]->GetY();
1904         z     = krd[id]->GetZ();
1905         fitter.AddPoint(xy,z);
1906
1907         // Evaluate and get parameters
1908         fitter.Eval();
1909
1910         // how to get the other clusters now?
1911         // ...
1912
1913       }
1914     }
1915   }
1916 }
1917
1918 //____________________________________________________________________________________
1919 void AliToyMCReconstruction::InitStreamer(const char* /*addName*/, Int_t /*level*/)
1920 {
1921   //
1922   // initilise the debug streamer and set the logging level
1923   //   use a default naming convention
1924   //
1925
1926 }
1927
1928 //____________________________________________________________________________________
1929 void AliToyMCReconstruction::ConnectInputFile(const char* file)
1930 {
1931   //
1932   // connect the tree and event pointer from the input file
1933   //
1934
1935   delete fInputFile;
1936   fInputFile=0x0;
1937   fEvent=0x0;
1938   fTree=0;
1939
1940   fInputFile=TFile::Open(file);
1941   if (!fInputFile || !fInputFile->IsOpen() || fInputFile->IsZombie()) {
1942     delete fInputFile;
1943     fInputFile=0x0;
1944     printf("ERROR: couldn't open the file '%s'\n", file);
1945     return;
1946   }
1947
1948   fTree=(TTree*)fInputFile->Get("toyMCtree");
1949   if (!fTree) {
1950     printf("ERROR: couldn't read the 'toyMCtree' from file '%s'\n", file);
1951     return;
1952   }
1953
1954   fEvent=0x0;
1955   fTree->SetBranchAddress("event",&fEvent);
1956
1957   gROOT->cd();
1958
1959   SetupTrackMaps();
1960 }
1961
1962 //____________________________________________________________________________________
1963 void AliToyMCReconstruction::Cleanup()
1964 {
1965   //
1966   // Cleanup input data
1967   //
1968
1969   delete fStreamer;
1970   fStreamer=0x0;
1971
1972   delete fEvent;
1973   fEvent = 0x0;
1974
1975   delete fTree;
1976   fTree=0x0;
1977
1978   delete fInputFile;
1979   fInputFile=0x0;
1980
1981 }
1982
1983 //____________________________________________________________________________________
1984 void AliToyMCReconstruction::SetupTrackMaps()
1985 {
1986   //
1987   //
1988   //
1989
1990   fMapTrackEvent.Delete();
1991   fMapTrackTrackInEvent.Delete();
1992
1993   if (!fTree) {
1994     AliError("Tree not connected");
1995     return;
1996   }
1997
1998   for (Int_t iev=0; iev<fTree->GetEntries(); ++iev) {
1999     fTree->GetEvent(iev);
2000     if (!fEvent) continue;
2001
2002     const Int_t ntracks=fEvent->GetNumberOfTracks();
2003     if (fMapTrackEvent.GetSize()+ntracks>=fMapTrackEvent.Capacity()) fMapTrackEvent.Expand(2*fMapTrackEvent.Capacity());
2004     if (fMapTrackTrackInEvent.GetSize()+ntracks>=fMapTrackTrackInEvent.Capacity()) fMapTrackTrackInEvent.Expand(2*fMapTrackTrackInEvent.Capacity());
2005
2006     for (Int_t itrack=0; itrack<ntracks; ++itrack){
2007       Int_t label=fEvent->GetTrack(itrack)->GetUniqueID();
2008
2009       fMapTrackEvent.Add(label,iev);
2010       fMapTrackTrackInEvent.Add(label,itrack);
2011     }
2012   }
2013
2014 }