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 , fNmaxEvents(-1)
  47 , fTime0(-1)
  48 , fCreateT0seed(kFALSE)
  49 , fStreamer(0x0)
  50 , fInputFile(0x0)
  51 , fTree(0x0)
  52 , fEvent(0x0)
  53 , fTPCParam(0x0)
  54 , fTPCCorrection(0x0)
  55 , fkNSectorInner(18) // hard-coded to avoid loading the parameters before
  56 , fInnerSectorArray(0x0)
  57 , fkNSectorOuter(18) // hard-coded to avoid loading the parameters before
  58 , fOuterSectorArray(0x0)
  59 , fAllClusters("AliTPCclusterMI",10000)
  60 , fMapTrackEvent(10000)
  61 , fMapTrackTrackInEvent(10000)
  62 , fIsAC(kFALSE)
  63 {
  64   //
  65   //  ctor
  66   //
  67   fTPCParam=AliTPCcalibDB::Instance()->GetParameters();
  68
  69 }
  70
  71 //____________________________________________________________________________________
  72 AliToyMCReconstruction::~AliToyMCReconstruction()
  73 {
  74   //
  75   //  dtor
  76   //
  77
  78   Cleanup();
  79 }
  80
  81 //____________________________________________________________________________________
  82 void AliToyMCReconstruction::RunReco(const char* file, Int_t nmaxEv)
  83 {
  84   //
  85   // Recostruction from associated clusters
  86   //
  87
  88   ConnectInputFile(file, nmaxEv);
  89   if (!fTree) return;
  90
  91   InitStreamer(".debug");
  92
  93   gROOT->cd();
  94
  95   static AliExternalTrackParam dummySeedT0;
  96   static AliExternalTrackParam dummySeed;
  97   static AliExternalTrackParam dummyTrack;
  98
  99   AliExternalTrackParam t0seed;
 100   AliExternalTrackParam seed;
 101   AliExternalTrackParam track;
 102   AliExternalTrackParam trackITS;
 103   AliExternalTrackParam tOrig;
 104   AliExternalTrackParam tOrigITS;
 105
 106   AliExternalTrackParam *dummy;
 107
 108   Int_t maxev=fTree->GetEntries();
 109   if (nmaxEv>0&&nmaxEv<maxev) maxev=nmaxEv;
 110
 111   for (Int_t iev=0; iev<maxev; ++iev){
 112     printf("==============  Processing Event %6d =================\n",iev);
 113     fTree->GetEvent(iev);
 114     for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
 115 //       printf(" > ======  Processing Track %6d ========  \n",itr);
 116       const AliToyMCTrack *tr=fEvent->GetTrack(itr);
 117       tOrig = *tr;
 118
 119
 120       // set dummy
 121       t0seed    = dummySeedT0;
 122       seed      = dummySeed;
 123       track     = dummyTrack;
 124
 125       Float_t z0=fEvent->GetZ();
 126       Float_t t0=fEvent->GetT0();
 127
 128       Float_t vDrift=GetVDrift();
 129       Float_t zLength=GetZLength(0);
 130
 131       // crate time0 seed, steered by fCreateT0seed
 132 //       printf("t0 seed\n");
 133       fTime0=-1.;
 134       fCreateT0seed=kTRUE;
 135       dummy = GetSeedFromTrack(tr);
 136
 137       if (dummy) {
 138         t0seed = *dummy;
 139         delete dummy;
 140
 141         // crate real seed using the time 0 from the first seed
 142         // set fCreateT0seed now to false to get the seed in z coordinates
 143         fTime0 = t0seed.GetZ()-zLength/vDrift;
 144         fCreateT0seed = kFALSE;
 145 //         printf("seed (%.2g)\n",fTime0);
 146         dummy  = GetSeedFromTrack(tr);
 147         if (dummy) {
 148           seed = *dummy;
 149           delete dummy;
 150
 151           // create fitted track
 152           if (fDoTrackFit){
 153 //             printf("track\n");
 154             dummy = GetFittedTrackFromSeed(tr, &seed);
 155             track = *dummy;
 156             delete dummy;
 157           }
 158
 159           // Copy original track and fitted track
 160           // for extrapolation to ITS last layer
 161           tOrigITS = *tr;
 162           trackITS = track;
 163
 164           // propagate seed to 0
 165           const Double_t kMaxSnp = 0.85;
 166           const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 167           AliTrackerBase::PropagateTrackTo(&seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
 168
 169           // propagate original track to ITS last layer
 170           Double_t lastLayerITS = 43.0; // same as in AliToyMCEventGenerator::MakeITSClusters (hard coded)
 171           AliTrackerBase::PropagateTrackTo(&tOrigITS,lastLayerITS,kMass,1,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
 172
 173           // rotate fitted track to the frame of the original track and propagate to same reference
 174           AliTrackerBase::PropagateTrackTo(&trackITS,lastLayerITS,kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
 175           trackITS.Rotate(tOrigITS.GetAlpha());
 176           AliTrackerBase::PropagateTrackTo(&trackITS,lastLayerITS,kMass,1,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
 177
 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         "trackITS.="   << &trackITS       <<
 202         "tOrigITS.="   << &tOrigITS       <<
 203         "\n";
 204       }
 205
 206
 207     }
 208   }
 209
 210   Cleanup();
 211 }
 212
 213
 214 //____________________________________________________________________________________
 215 void AliToyMCReconstruction::RunRecoAllClusters(const char* file, Int_t nmaxEv)
 216 {
 217   //
 218   // Reconstruction for seed from associated clusters, but array of clusters:
 219   // Step 1) Filling of cluster arrays
 220   // Step 2) Seeding from clusters associated to tracks
 221   // Step 3) Free track reconstruction using all clusters
 222   //
 223
 224   TFile f(file);
 225   if (!f.IsOpen() || f.IsZombie()) {
 226     printf("ERROR: couldn't open the file '%s'\n", file);
 227     return;
 228   }
 229
 230  fTree=(TTree*)f.Get("toyMCtree");
 231   if (!fTree) {
 232     printf("ERROR: couldn't read the 'toyMCtree' from file '%s'\n", file);
 233     return;
 234   }
 235
 236   fEvent=0x0;
 237   fTree->SetBranchAddress("event",&fEvent);
 238
 239   // read spacecharge from the Userinfo ot the tree
 240   InitSpaceCharge();
 241
 242   TString debugName=file;
 243   debugName.ReplaceAll(".root","");
 244   debugName.Append(Form(".%1d.%1d_%1d_%1d_%03d_%02d",
 245                         fUseMaterial,fIdealTracking,fClusterType,
 246                         Int_t(fCorrectionType),fSeedingRow,fSeedingDist));
 247   debugName.Append(".allClusters.debug.root");
 248
 249   gSystem->Exec(Form("test -f %s && rm %s", debugName.Data(), debugName.Data()));
 250   if (!fStreamer) fStreamer=new TTreeSRedirector(debugName.Data());
 251
 252   gROOT->cd();
 253
 254   static AliExternalTrackParam dummySeedT0;
 255   static AliExternalTrackParam dummySeed;
 256   static AliExternalTrackParam dummyTrack;
 257
 258   AliExternalTrackParam t0seed;
 259   AliExternalTrackParam seed;
 260   AliExternalTrackParam track;
 261   AliExternalTrackParam tOrig;
 262
 263   AliExternalTrackParam *dummy;
 264
 265   Int_t maxev=fTree->GetEntries();
 266   if (nmaxEv>0&&nmaxEv<maxev) maxev=nmaxEv;
 267
 268   // ===========================================================================================
 269   // Loop 1: Fill AliTPCtrackerSector structure
 270   // ===========================================================================================
 271   FillSectorStructure(maxev);
 272
 273   // settings (TODO: find the correct settings)
 274   AliTPCRecoParam *tpcRecoParam = new AliTPCRecoParam();
 275   tpcRecoParam->SetDoKinks(kFALSE);
 276   AliTPCcalibDB::Instance()->GetTransform()->SetCurrentRecoParam(tpcRecoParam);
 277   //tpcRecoParam->Print();
 278
 279   // need AliTPCReconstructor for parameter settings in AliTPCtracker
 280   AliTPCReconstructor *tpcRec   = new AliTPCReconstructor();
 281   tpcRec->SetRecoParam(tpcRecoParam);
 282
 283
 284   // ===========================================================================================
 285   // Loop 2: Seeding from clusters associated to tracks
 286   // TODO: Implement tracking from given seed!
 287   // ===========================================================================================
 288   for (Int_t iev=0; iev<maxev; ++iev){
 289     printf("==============  Processing Event %6d =================\n",iev);
 290     fTree->GetEvent(iev);
 291     for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
 292       printf(" > ======  Processing Track %6d ========  \n",itr);
 293       const AliToyMCTrack *tr=fEvent->GetTrack(itr);
 294       tOrig = *tr;
 295
 296
 297       // set dummy
 298       t0seed    = dummySeedT0;
 299       seed      = dummySeed;
 300       track     = dummyTrack;
 301
 302       Float_t z0=fEvent->GetZ();
 303       Float_t t0=fEvent->GetT0();
 304
 305       Float_t vDrift=GetVDrift();
 306       Float_t zLength=GetZLength(0);
 307
 308       Int_t nClus = 0;
 309
 310       // crate time0 seed, steered by fCreateT0seed
 311       printf("t0 seed\n");
 312       fTime0=-1.;
 313       fCreateT0seed=kTRUE;
 314       dummy = GetSeedFromTrack(tr);
 315
 316       if (dummy) {
 317         t0seed = *dummy;
 318         delete dummy;
 319
 320         // crate real seed using the time 0 from the first seed
 321         // set fCreateT0seed now to false to get the seed in z coordinates
 322         fTime0 = t0seed.GetZ()-zLength/vDrift;
 323         fCreateT0seed = kFALSE;
 324         printf("seed (%.2g)\n",fTime0);
 325         dummy  = GetSeedFromTrack(tr);
 326         if (dummy) {
 327           seed = *dummy;
 328           delete dummy;
 329
 330           // create fitted track
 331           if (fDoTrackFit){
 332             printf("track\n");
 333             dummy = GetFittedTrackFromSeedAllClusters(tr, &seed,nClus);
 334             track = *dummy;
 335             delete dummy;
 336           }
 337
 338           // propagate seed to 0
 339           const Double_t kMaxSnp = 0.85;
 340           const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 341           AliTrackerBase::PropagateTrackTo(&seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
 342
 343         }
 344       }
 345
 346       Int_t ctype(fCorrectionType);
 347
 348       if (fStreamer) {
 349         (*fStreamer) << "Tracks" <<
 350         "iev="         << iev             <<
 351         "z0="          << z0              <<
 352         "t0="          << t0              <<
 353         "fTime0="      << fTime0          <<
 354         "itr="         << itr             <<
 355         "clsType="     << fClusterType    <<
 356         "corrType="    << ctype           <<
 357         "seedRow="     << fSeedingRow     <<
 358         "seedDist="    << fSeedingDist    <<
 359         "vDrift="      << vDrift          <<
 360         "zLength="     << zLength         <<
 361         "nClus="       << nClus           <<
 362         "t0seed.="     << &t0seed         <<
 363         "seed.="       << &seed           <<
 364         "track.="      << &track          <<
 365         "tOrig.="      << &tOrig          <<
 366         "\n";
 367       }
 368
 369
 370     }
 371   }
 372
 373
 374   delete fStreamer;
 375   fStreamer=0x0;
 376
 377   delete fEvent;
 378   fEvent = 0x0;
 379
 380   delete fTree;
 381   fTree=0x0;
 382   f.Close();
 383 }
 384
 385 //____________________________________________________________________________________
 386 void AliToyMCReconstruction::RunRecoAllClustersStandardTracking(const char* file, Int_t nmaxEv)
 387 {
 388   //
 389   // Reconstruction for seed from associated clusters, but array of clusters
 390   // Step 1) Filling of cluster arrays
 391   // Step 2) Use the standard tracking: AliTPCtracker::Clusters2Tracks();
 392   //
 393
 394   TFile f(file);
 395   if (!f.IsOpen() || f.IsZombie()) {
 396     printf("ERROR: couldn't open the file '%s'\n", file);
 397     return;
 398   }
 399
 400  fTree=(TTree*)f.Get("toyMCtree");
 401   if (!fTree) {
 402     printf("ERROR: couldn't read the 'toyMCtree' from file '%s'\n", file);
 403     return;
 404   }
 405
 406   fEvent=0x0;
 407   fTree->SetBranchAddress("event",&fEvent);
 408
 409   // read spacecharge from the Userinfo ot the tree
 410   InitSpaceCharge();
 411
 412   TString debugName=file;
 413   debugName.ReplaceAll(".root","");
 414   debugName.Append(Form(".%1d.%1d_%1d_%1d_%03d_%02d",
 415                         fUseMaterial,fIdealTracking,fClusterType,
 416                         Int_t(fCorrectionType),fSeedingRow,fSeedingDist));
 417   debugName.Append(".allClusters.debug.root");
 418
 419   gSystem->Exec(Form("test -f %s && rm %s", debugName.Data(), debugName.Data()));
 420   if (!fStreamer) fStreamer=new TTreeSRedirector(debugName.Data());
 421
 422   gROOT->cd();
 423
 424   AliExternalTrackParam t0seed;
 425   AliExternalTrackParam seed;
 426   AliExternalTrackParam track;
 427   AliExternalTrackParam tOrig;
 428   AliToyMCTrack tOrigToy;
 429
 430   AliExternalTrackParam *dummy;
 431   AliTPCseed            *seedBest;
 432   AliTPCseed            *seedTmp;
 433   AliTPCclusterMI       *cluster;
 434
 435   Int_t maxev=fTree->GetEntries();
 436   if (nmaxEv>0&&nmaxEv<maxev) maxev=nmaxEv;
 437
 438
 439   // ===========================================================================================
 440   // Loop 1: Fill AliTPCtrackerSector structure
 441   // ===========================================================================================
 442   FillSectorStructure(maxev);
 443
 444   // ===========================================================================================
 445   // Loop 2: Use the TPC tracker for seeding (MakeSeeds3)
 446   // TODO: - check tracking configuration
 447   //       - add clusters and original tracks to output (how?)
 448   // ===========================================================================================
 449
 450   // settings (TODO: find the correct settings)
 451   AliTPCRecoParam *tpcRecoParam = new AliTPCRecoParam();
 452   tpcRecoParam->SetDoKinks(kFALSE);
 453   AliTPCcalibDB::Instance()->GetTransform()->SetCurrentRecoParam(tpcRecoParam);
 454   //tpcRecoParam->Print();
 455
 456   // need AliTPCReconstructor for parameter settings in AliTPCtracker
 457   AliTPCReconstructor *tpcRec   = new AliTPCReconstructor();
 458   tpcRec->SetRecoParam(tpcRecoParam);
 459
 460   // AliTPCtracker
 461   AliTPCtracker *tpcTracker = new AliTPCtracker(fTPCParam);
 462   tpcTracker->SetDebug(10);
 463
 464   // set sector arrays
 465   tpcTracker->SetTPCtrackerSectors(fInnerSectorArray,fOuterSectorArray);
 466   tpcTracker->LoadInnerSectors();
 467   tpcTracker->LoadOuterSectors();
 468
 469   // seeding
 470   static TObjArray arrTracks;
 471   TObjArray * arr    = &arrTracks;
 472   TObjArray * seeds  = new TObjArray;
 473
 474   // cuts for seeding
 475 //   Float_t cuts[4];
 476 //   cuts[0]=0.0070;  // cuts[0]   - fP4 cut
 477 //   cuts[1] = 1.5;   // cuts[1]   - tan(phi) cut
 478 //   cuts[2] = 3.;    // cuts[2]   - zvertex cut
 479 //   cuts[3] = 3.;    // cuts[3]   - fP3 cut
 480
 481   // rows for seeding
 482   Int_t lowerRow = fSeedingRow;
 483   Int_t upperRow = fSeedingRow+2*fSeedingDist;
 484   const AliTPCROC * roc      = AliTPCROC::Instance();
 485   const Int_t kNRowsInnerTPC = roc->GetNRows(0);
 486   const Int_t kNRowsTPC      = kNRowsInnerTPC + roc->GetNRows(36);
 487   if(lowerRow < kNRowsInnerTPC){
 488     Printf("Seeding row requested (%d) is lower than kNRowsInnerTPC --> use %d",lowerRow,kNRowsInnerTPC);
 489     lowerRow = kNRowsInnerTPC;
 490     upperRow = lowerRow + 20;
 491   }
 492   if(upperRow >= kNRowsTPC){
 493     Printf("Seeding row requested (%d) is larger than kNRowsTPC --> use %d",upperRow,kNRowsTPC-1);
 494     upperRow = kNRowsTPC-1;
 495     lowerRow = upperRow-20;
 496   }
 497
 498   // do the seeding
 499   for (Int_t sec=0;sec<fkNSectorOuter;sec++){
 500     //
 501     //tpcTracker->MakeSeeds3(arr, sec,upperRow,lowerRow,cuts,-1,1);
 502     MakeSeeds(arr, sec,upperRow,lowerRow); // own function (based on TLinearFitter)
 503     //tpcTracker->SumTracks(seeds,arr);
 504     //tpcTracker->SignClusters(seeds,3.0,3.0);
 505
 506   }
 507
 508   Printf("After seeding we have %d tracks",seeds->GetEntriesFast());
 509
 510   // Standard tracking
 511   tpcTracker->SetSeeds(seeds);
 512   tpcTracker->PropagateForward();
 513   Printf("After trackinging we have %d tracks",seeds->GetEntriesFast());
 514 return;
 515
 516   // Loop over all input tracks and connect to found seeds
 517   for (Int_t iev=0; iev<maxev; ++iev){
 518     printf("==============  Fill Tracks: Processing Event %6d  =================\n",iev);
 519     fTree->GetEvent(iev);
 520     for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
 521       printf(" > ======  Fill Tracks: Processing Track %6d  ========  \n",itr);
 522       const AliToyMCTrack *tr=fEvent->GetTrack(itr);
 523       tOrig = *tr;
 524       tOrigToy = *tr;
 525
 526       Float_t z0=fEvent->GetZ();
 527       Float_t t0=fEvent->GetT0();
 528       Float_t vDrift=GetVDrift();
 529       Float_t zLength=GetZLength(0);
 530
 531       // find the corresponding seed (and track)
 532       Int_t trackID            = tr->GetUniqueID();
 533       Int_t nClustersMC        = tr->GetNumberOfSpacePoints();      // number of findable clusters (ideal)
 534       if(fClusterType==1)
 535             nClustersMC        = tr->GetNumberOfDistSpacePoints();  // number of findable clusters (distorted)
 536 //       Int_t idxSeed            = 0; // index of best seed (best is with maximum number of clusters with correct ID)
 537       Int_t nSeeds             = 0; // number of seeds for MC track
 538       Int_t nSeedClusters      = 0; // number of clusters for best seed
 539       Int_t nSeedClustersTmp   = 0; // number of clusters for current seed
 540       Int_t nSeedClustersID    = 0; // number of clusters with correct ID for best seed
 541       Int_t nSeedClustersIDTmp = 0; // number of clusters with correct ID for current seed
 542       for(Int_t iSeeds = 0; iSeeds < seeds->GetEntriesFast(); ++iSeeds){
 543
 544         // set current seed and reset counters
 545         seedTmp = (AliTPCseed*)(seeds->At(iSeeds));
 546         nSeedClustersTmp   = 0;
 547         nSeedClustersIDTmp = 0;
 548
 549         if(!seedTmp) continue;
 550
 551         // loop over all rows
 552         for(Int_t iRow = seedTmp->GetRow(); iRow < seedTmp->GetRow() + seedTmp->GetNumberOfClustersIndices(); iRow++ ){
 553
 554           // get cluster and increment counters
 555           cluster = seedTmp->GetClusterFast(iRow);
 556           if(cluster){
 557             nSeedClustersTmp++;
 558             if(cluster->GetLabel(0)==trackID){
 559               nSeedClustersIDTmp++;
 560             }
 561           }
 562         }
 563
 564         // if number of corresponding clusters > 0,
 565         // increase nSeeds
 566         if(nSeedClustersTmp > 0){
 567           nSeeds++;
 568         }
 569
 570         // if number of corresponding clusters bigger than current nSeedClusters,
 571         // take this seed as the best one
 572         if(nSeedClustersIDTmp > nSeedClustersID){
 573 //        idxSeed  = iSeeds;
 574           seedBest = seedTmp;
 575           nSeedClusters   = nSeedClustersTmp;   // number of correctly assigned clusters
 576           nSeedClustersID = nSeedClustersIDTmp; // number of all clusters
 577         }
 578
 579       }
 580
 581       // cluster to track association (commented out, when used standard tracking)
 582       if (nSeeds>0&&nSeedClusters>0) {
 583         t0seed = (AliExternalTrackParam)*seedBest;
 584 //              fTime0 = t0seed.GetZ()-zLength/vDrift;
 585         // get the refitted track from the seed
 586         // this will also set the fTime0 from the seed extrapolation
 587         dummy=GetRefittedTrack(*seedBest);
 588         track=*dummy;
 589         delete dummy;
 590         //printf("seed (%.2g): %d seeds with %d clusters\n",fTime0,nSeeds,nSeedClusters);
 591
 592         //      // cluster to track association for all good seeds
 593         //      // set fCreateT0seed to true to get the seed in time coordinates
 594         //      fCreateT0seed = kTRUE;
 595         //      dummy = ClusterToTrackAssociation(seedBest,trackID,nClus);
 596
 597         //// propagate track to 0
 598         //const Double_t kMaxSnp = 0.85;
 599         //const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 600         //AliTrackerBase::PropagateTrackTo(&track,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
 601
 602       }
 603
 604       Int_t ctype(fCorrectionType);
 605
 606       if (fStreamer) {
 607         (*fStreamer) << "Tracks" <<
 608           "iev="             << iev             <<
 609           "z0="              << z0              <<
 610           "t0="              << t0              <<
 611           "fTime0="          << fTime0          <<
 612           "itr="             << itr             <<
 613           "clsType="         << fClusterType    <<
 614           "corrType="        << ctype           <<
 615           "seedRow="         << fSeedingRow     <<
 616           "seedDist="        << fSeedingDist    <<
 617           "vDrift="          << vDrift          <<
 618           "zLength="         << zLength         <<
 619           "nClustersMC="     << nClustersMC     <<
 620           "nSeeds="          << nSeeds          <<
 621           "nSeedClusters="   << nSeedClusters   <<
 622           "nSeedClustersID=" << nSeedClustersID <<
 623           "t0seed.="         << &t0seed         <<
 624           "track.="          << &track          <<
 625           "tOrig.="          << &tOrig          <<
 626           "tOrigToy.="       << &tOrigToy       <<
 627           "\n";
 628       }
 629     }
 630   }
 631
 632
 633   delete fStreamer;
 634   fStreamer=0x0;
 635
 636   delete fEvent;
 637   fEvent = 0x0;
 638
 639   delete fTree;
 640   fTree=0x0;
 641   f.Close();
 642 }
 643
 644
 645 //____________________________________________________________________________________
 646 void AliToyMCReconstruction::RunFullTracking(const char* file, Int_t nmaxEv)
 647 {
 648   //
 649   //
 650   //
 651
 652   ConnectInputFile(file,nmaxEv);
 653   if (!fTree) return;
 654
 655   InitStreamer(".fullTracking");
 656
 657   FillSectorStructureAC();
 658
 659   AliTPCReconstructor::SetStreamLevel(0);
 660
 661   TObjArray seeds;
 662   seeds.SetOwner();
 663   Int_t lowerRow=130;
 664   Int_t upperRow=150;
 665   if (lowerRow>upperRow){
 666     Int_t tmp=lowerRow;
 667     lowerRow=upperRow;
 668     upperRow=tmp;
 669   }
 670
 671   // seeding.
 672   // NOTE: the z position is set to GetTimeBin*vDrift
 673   //       therefore it is not possible to simply propagate
 674   //       the track using AliTrackerBase::Propagate, since a
 675   //       wrong B-Field will be assinged...
 676   printf("Run Seeding in %3d - %3d\n",lowerRow,upperRow);
 677   for (Int_t sec=0;sec<36;sec++){
 678     printf(" in sector: %d\n",sec);
 679     Int_t nAdded=MakeSeeds2(&seeds, sec,lowerRow,upperRow);
 680     printf("  -> Added Seeds: %d\n",nAdded);
 681     nAdded=MakeSeeds2(&seeds, sec,lowerRow-2,upperRow-2);
 682     printf("  -> Added Seeds: %d\n",nAdded);
 683     nAdded=MakeSeeds2(&seeds, sec,lowerRow-4,upperRow-4);
 684     printf("  -> Added Seeds: %d\n",nAdded);
 685   }
 686
 687   printf("Run Tracking in %3d - %3d\n",lowerRow,upperRow);
 688   Int_t firstSeed=0;
 689   for (Int_t iseed=firstSeed; iseed<seeds.GetEntriesFast();++iseed) ClusterToTrackAssociation(*(AliTPCseed*)seeds.UncheckedAt(iseed));
 690   //first seed is used to not run the tracking twice on a seed
 691   firstSeed=seeds.GetEntriesFast();
 692 //   DumpTrackInfo(&seeds);
 693
 694   lowerRow=110;
 695   upperRow=130;
 696
 697   printf("Run Seeding in %3d - %3d\n",lowerRow,upperRow);
 698   for (Int_t sec=0;sec<36;sec++){
 699     printf(" in sector: %d\n",sec);
 700     Int_t nAdded=MakeSeeds2(&seeds, sec,lowerRow,upperRow);
 701     printf("  -> Added Seeds: %d\n",nAdded);
 702   }
 703   printf("Run Tracking in %3d - %3d\n",lowerRow,upperRow);
 704   for (Int_t iseed=firstSeed; iseed<seeds.GetEntriesFast();++iseed) ClusterToTrackAssociation(*(AliTPCseed*)seeds.UncheckedAt(iseed));
 705   firstSeed=seeds.GetEntriesFast();
 706
 707   //now seeding also at more central rows with shorter seeds
 708   lowerRow=70;
 709   upperRow=90;
 710
 711   printf("Run Seeding in %3d - %3d\n",lowerRow,upperRow);
 712   for (Int_t sec=0;sec<36;sec++){
 713     printf(" in sector: %d\n",sec);
 714     Int_t nAdded=MakeSeeds2(&seeds, sec,lowerRow,upperRow);
 715     printf("  -> Added Seeds: %d\n",nAdded);
 716   }
 717   printf("Run Tracking in %3d - %3d\n",lowerRow,upperRow);
 718   for (Int_t iseed=firstSeed; iseed<seeds.GetEntriesFast();++iseed) ClusterToTrackAssociation(*(AliTPCseed*)seeds.UncheckedAt(iseed));
 719   firstSeed=seeds.GetEntriesFast();
 720
 721   //shorter seeds
 722   Int_t startUpper=upperRow-10;
 723   Int_t startLower=lowerRow-5;
 724   for (Int_t sec=0;sec<36;sec++){
 725     upperRow=startUpper;
 726     lowerRow=startLower;
 727     printf(" in sector: %d\n",sec);
 728     while (lowerRow>0){
 729       printf("Run Seeding in %3d - %3d\n",lowerRow,upperRow);
 730       Int_t nAdded=MakeSeeds2(&seeds, sec,lowerRow,upperRow);
 731       printf("  -> Added Seeds: %d\n",nAdded);
 732       for (Int_t iseed=firstSeed; iseed<seeds.GetEntriesFast();++iseed) ClusterToTrackAssociation(*(AliTPCseed*)seeds.UncheckedAt(iseed));
 733       firstSeed=seeds.GetEntriesFast();
 734       lowerRow-=5;
 735       upperRow-=5;
 736     }
 737   }
 738
 739   //track remaining
 740
 741   DumpTrackInfo(&seeds);
 742
 743 //   TObjArray seedsCentral2;
 744 //   lowerRow=45;
 745 //   upperRow=62;
 746 //
 747 //   for (Int_t sec=0;sec<36;sec++){
 748 //     Int_t nAdded=MakeSeeds2(&seedsCentral2, sec,lowerRow,upperRow);
 749 //     printf("  -> Added Seeds: %d\n",nAdded);
 750 //     nAdded=MakeSeeds2(&seedsCentral2, sec,lowerRow-2,upperRow-2);
 751 //     printf("  -> Added Seeds: %d\n",nAdded);
 752 //     nAdded=MakeSeeds2(&seedsCentral2, sec,lowerRow-4,upperRow-4);
 753 //     printf("  -> Added Seeds: %d\n",nAdded);
 754 //   }
 755 //   for (Int_t iseed=0; iseed<seedsCentral2.GetEntriesFast();++iseed) ClusterToTrackAssociation(*(AliTPCseed*)seedsCentral2.UncheckedAt(iseed));
 756 //   DumpTrackInfo(&seedsCentral2);
 757
 758   //dump clusters
 759   (*fStreamer) << "clusters" <<
 760   "cl.=" << &fAllClusters << "\n";
 761
 762   Cleanup();
 763 }
 764
 765 //____________________________________________________________________________________
 766 AliExternalTrackParam* AliToyMCReconstruction::GetSeedFromTrack(const AliToyMCTrack * const tr, Bool_t forceSeed)
 767 {
 768   //
 769   // if we don't have a valid time0 informaion (fTime0) available yet
 770   // assume we create a seed for the time0 estimate
 771   //
 772
 773   // seed point informaion
 774   AliTrackPoint    seedPoint[3];
 775   const AliTPCclusterMI *seedCluster[3]={0x0,0x0,0x0};
 776
 777   // number of clusters to loop over
 778   const Int_t ncls=(fClusterType==0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
 779   if (ncls<3){
 780     AliError(Form("Not enough points to create a seed: %d",ncls));
 781     return 0x0;
 782   }
 783   UChar_t nextSeedRow=fSeedingRow;
 784   Int_t   nseeds=0;
 785
 786   //assumes sorted clusters
 787   if (forceSeed){
 788     // force the seed creation, using the first, middle and last cluster
 789     Int_t npoints[3]={0,ncls/2,ncls-1};
 790     for (Int_t icl=0;icl<3;++icl){
 791       const AliTPCclusterMI *cl=tr->GetSpacePoint(npoints[icl]);
 792       if (fClusterType==1) cl=tr->GetDistortedSpacePoint(npoints[icl]);
 793       seedCluster[nseeds]=cl;
 794       SetTrackPointFromCluster(cl, seedPoint[nseeds]);
 795       ++nseeds;
 796     }
 797   }else{
 798     // create seeds according to the reco settings
 799     for (Int_t icl=0;icl<ncls;++icl) {
 800       const AliTPCclusterMI *cl=tr->GetSpacePoint(icl);
 801       if (fClusterType==1) cl=tr->GetDistortedSpacePoint(icl);
 802       if (!cl) continue;
 803       // use row in sector
 804       const UChar_t row=cl->GetRow() + 63*(cl->GetDetector()>35);
 805       // skip clusters without proper pad row
 806       if (row>200) continue;
 807
 808       //check seeding row
 809       // if we are in the last row and still miss a seed we use the last row
 810       //   even if the row spacing will not be equal
 811       if (row>=nextSeedRow || icl==ncls-1){
 812         seedCluster[nseeds]=cl;
 813         SetTrackPointFromCluster(cl, seedPoint[nseeds]);
 814       ++nseeds;
 815       nextSeedRow+=fSeedingDist;
 816
 817       if (nseeds==3) break;
 818       }
 819     }
 820   }
 821
 822   // check we really have 3 seeds
 823   if (nseeds!=3) {
 824     AliError(Form("Seeding failed for parameters %d, %d\n",fSeedingDist,fSeedingRow));
 825     return 0x0;
 826   }
 827
 828   // do cluster correction for fCorrectionType:
 829   //   0 - no correction
 830   //   1 - TPC center
 831   //   2 - average eta
 832   //   3 - ideal
 833   // assign the cluster abs time as z component to all seeds
 834   for (Int_t iseed=0; iseed<3; ++iseed) {
 835     Float_t xyz[3]={0,0,0};
 836     seedPoint[iseed].GetXYZ(xyz);
 837     const Float_t r=TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
 838
 839     const Int_t sector=seedCluster[iseed]->GetDetector();
 840     const Int_t sign=1-2*((sector/18)%2);
 841
 842     if ( (fClusterType == 1) && (fCorrectionType != kNoCorrection) ) {
 843 //       printf("correction type: %d\n",(Int_t)fCorrectionType);
 844
 845       // the settings below are for the T0 seed
 846       // for known T0 the z position is already calculated in SetTrackPointFromCluster
 847       if ( fCreateT0seed ){
 848         if ( fCorrectionType == kTPCCenter  ) xyz[2] = 125.*sign;
 849         //!!! TODO: is this the correct association?
 850         if ( fCorrectionType == kAverageEta ) xyz[2] = TMath::Tan(45./2.*TMath::DegToRad())*r*sign;
 851       }
 852
 853       if ( fCorrectionType == kIdeal      ) xyz[2] = seedCluster[iseed]->GetZ();
 854
 855       //!!! TODO: to be replaced with the proper correction
 856       fTPCCorrection->CorrectPoint(xyz, seedCluster[iseed]->GetDetector());
 857     }
 858
 859     // after the correction set the time bin as z-Position in case of a T0 seed
 860     if ( fCreateT0seed )
 861       xyz[2]=seedCluster[iseed]->GetTimeBin();
 862
 863     seedPoint[iseed].SetXYZ(xyz);
 864   }
 865
 866   const Double_t kMaxSnp = 0.85;
 867   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 868
 869   AliExternalTrackParam *seed = AliTrackerBase::MakeSeed(seedPoint[0], seedPoint[1], seedPoint[2]);
 870   seed->ResetCovariance(10);
 871
 872   if (fCreateT0seed){
 873     // if fTime0 < 0 we assume that we create a seed for the T0 estimate
 874     AliTrackerBase::PropagateTrackTo(seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
 875     if (TMath::Abs(seed->GetX())>3) {
 876 //       printf("Could not propagate track to 0, %.2f, %.2f, %.2f\n",seed->GetX(),seed->GetAlpha(),seed->Pt());
 877     }
 878   }
 879
 880   return seed;
 881
 882 }
 883
 884
 885 //____________________________________________________________________________________
 886 void AliToyMCReconstruction::SetTrackPointFromCluster(const AliTPCclusterMI *cl, AliTrackPoint &p )
 887 {
 888   //
 889   // make AliTrackPoint out of AliTPCclusterMI
 890   //
 891
 892   if (!cl) return;
 893     Float_t xyz[3]={0.,0.,0.};
 894   //   ClusterToSpacePoint(cl,xyz);
 895   //   cl->GetGlobalCov(cov);
 896   //TODO: what to do with the covariance matrix???
 897   //TODO: the problem is that it is used in GetAngle in AliTrackPoint
 898   //TODO: which is used by AliTrackerBase::MakeSeed to get alpha correct ...
 899   //TODO: for the moment simply assign 1 permill squared
 900   // in AliTrackPoint the cov is xx, xy, xz, yy, yz, zz
 901   //   Float_t cov[6]={xyz[0]*xyz[0]*1e-6,xyz[0]*xyz[1]*1e-6,xyz[0]*xyz[2]*1e-6,
 902   //                   xyz[1]*xyz[1]*1e-6,xyz[1]*xyz[2]*1e-6,xyz[2]*xyz[2]*1e-6};
 903   //   cl->GetGlobalXYZ(xyz);
 904   //   cl->GetGlobalCov(cov);
 905   // voluem ID to add later ....
 906   //   p.SetXYZ(xyz);
 907   //   p.SetCov(cov);
 908 //   AliTrackPoint *tp=const_cast<AliTPCclusterMI*>(cl)->MakePoint(p);
 909 //   p=*tp;
 910 //   delete tp;
 911   const_cast<AliTPCclusterMI*>(cl)->MakePoint(p);
 912   //   cl->Print();
 913   //   p.Print();
 914   p.SetVolumeID(cl->GetDetector());
 915
 916
 917   if ( !fCreateT0seed && !fIdealTracking ) {
 918     p.GetXYZ(xyz);
 919     const Int_t sector=cl->GetDetector();
 920     const Int_t sign=1-2*((sector/18)%2);
 921     const Float_t zT0=( GetZLength(sector) - (cl->GetTimeBin()-fTime0)*GetVDrift() )*sign;
 922 //     printf(" z:  %.2f  %.2f\n",xyz[2],zT0);
 923     xyz[2]=zT0;
 924     p.SetXYZ(xyz);
 925   }
 926
 927
 928   //   p.Rotate(p.GetAngle()).Print();
 929 }
 930
 931 //____________________________________________________________________________________
 932 void AliToyMCReconstruction::ClusterToSpacePoint(const AliTPCclusterMI *cl, Float_t xyz[3])
 933 {
 934   //
 935   // convert the cluster to a space point in global coordinates
 936   //
 937   if (!cl) return;
 938   xyz[0]=cl->GetRow();
 939   xyz[1]=cl->GetPad();
 940   xyz[2]=cl->GetTimeBin(); // this will not be correct at all
 941   Int_t i[3]={0,cl->GetDetector(),cl->GetRow()};
 942   //   printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
 943   fTPCParam->Transform8to4(xyz,i);
 944   //   printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
 945   fTPCParam->Transform4to3(xyz,i);
 946   //   printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
 947   fTPCParam->Transform2to1(xyz,i);
 948   //   printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
 949 }
 950
 951 //____________________________________________________________________________________
 952 AliExternalTrackParam* AliToyMCReconstruction::GetFittedTrackFromSeed(const AliToyMCTrack *tr, const AliExternalTrackParam *seed)
 953 {
 954   //
 955   //
 956   //
 957
 958   // create track
 959   AliExternalTrackParam *track = new AliExternalTrackParam(*seed);
 960
 961   Int_t ncls=(fClusterType == 0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
 962
 963   const AliTPCROC * roc = AliTPCROC::Instance();
 964
 965   const Double_t kRTPC0  = roc->GetPadRowRadii(0,0);
 966   const Double_t kRTPC1  = roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
 967   const Double_t kMaxSnp = 0.85;
 968   const Double_t kMaxR   = 500.;
 969   const Double_t kMaxZ   = 500.;
 970
 971   //   const Double_t kMaxZ0=220;
 972 //   const Double_t kZcut=3;
 973
 974   const Double_t refX = tr->GetX();
 975
 976   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 977
 978   // loop over all other points and add to the track
 979   for (Int_t ipoint=ncls-1; ipoint>=0; --ipoint){
 980     AliTrackPoint pIn;
 981     const AliTPCclusterMI *cl=tr->GetSpacePoint(ipoint);
 982     if (fClusterType == 1) cl=tr->GetDistortedSpacePoint(ipoint);
 983     SetTrackPointFromCluster(cl, pIn);
 984     if (fCorrectionType != kNoCorrection){
 985       Float_t xyz[3]={0,0,0};
 986       pIn.GetXYZ(xyz);
 987 //       if ( fCorrectionType == kIdeal ) xyz[2] = cl->GetZ();
 988       fTPCCorrection->CorrectPoint(xyz, cl->GetDetector());
 989       pIn.SetXYZ(xyz);
 990     }
 991     // rotate the cluster to the local detector frame
 992     track->Rotate(((cl->GetDetector()%18)*20+10)*TMath::DegToRad());
 993     AliTrackPoint prot = pIn.Rotate(track->GetAlpha());   // rotate to the local frame - non distoted  point
 994     if (TMath::Abs(prot.GetX())<kRTPC0) continue;
 995     if (TMath::Abs(prot.GetX())>kRTPC1) continue;
 996     //
 997     Bool_t res=AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
 998
 999     if (!res) break;
1000
1001     if (TMath::Abs(track->GetZ())>kMaxZ) break;
1002     if (TMath::Abs(track->GetX())>kMaxR) break;
1003 //     if (TMath::Abs(track->GetZ())<kZcut)continue;
1004     //
1005     Double_t pointPos[2]={0,0};
1006     Double_t pointCov[3]={0,0,0};
1007     pointPos[0]=prot.GetY();//local y
1008     pointPos[1]=prot.GetZ();//local z
1009     pointCov[0]=prot.GetCov()[3];//simay^2
1010     pointCov[1]=prot.GetCov()[4];//sigmayz
1011     pointCov[2]=prot.GetCov()[5];//sigmaz^2
1012
1013     if (!track->Update(pointPos,pointCov)) {printf("no update\n"); break;}
1014   }
1015
1016   AliTrackerBase::PropagateTrackTo2(track,refX,kMass,5.,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
1017
1018   // rotate fittet track to the frame of the original track and propagate to same reference
1019   track->Rotate(tr->GetAlpha());
1020
1021   AliTrackerBase::PropagateTrackTo2(track,refX,kMass,1.,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
1022
1023   return track;
1024 }
1025
1026
1027 //____________________________________________________________________________________
1028 AliExternalTrackParam* AliToyMCReconstruction::GetFittedTrackFromSeedAllClusters(const AliToyMCTrack *tr, const AliExternalTrackParam *seed, Int_t &nClus)
1029 {
1030   //
1031   // Tracking for given seed on an array of clusters
1032   //
1033
1034   // create track
1035   AliExternalTrackParam *track = new AliExternalTrackParam(*seed);
1036
1037   const AliTPCROC * roc = AliTPCROC::Instance();
1038
1039   const Double_t kRTPC0    = roc->GetPadRowRadii(0,0);
1040   const Double_t kRTPC1    = roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
1041   const Int_t kNRowsTPC    = roc->GetNRows(0) + roc->GetNRows(36) - 1;
1042   const Int_t kIRowsTPC    = roc->GetNRows(0) - 1;
1043   const Double_t kMaxSnp   = 0.85;
1044   const Double_t kMaxR     = 500.;
1045   const Double_t kMaxZ     = 500.;
1046   const Double_t roady     = 100.;
1047   const Double_t roadz     = 100.;
1048
1049   const Double_t refX = tr->GetX();
1050
1051   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1052
1053   Int_t  secCur   = -1;
1054   UInt_t indexCur = 0;
1055   Double_t xCur, yCur, zCur = 0.;
1056
1057   Float_t vDrift = GetVDrift();
1058
1059   // first propagate seed to outermost row
1060   AliTrackerBase::PropagateTrackTo(track,kRTPC1,kMass,5,kFALSE,kMaxSnp);
1061
1062   // Loop over rows and find the cluster candidates
1063   for( Int_t iRow = kNRowsTPC; iRow >= 0; --iRow ){
1064
1065     // inner or outer sector
1066     Bool_t bInnerSector = kTRUE;
1067     if(iRow > kIRowsTPC) bInnerSector = kFALSE;
1068
1069     // nearest track point/cluster (to be found)
1070     AliTrackPoint nearestPoint;
1071     AliTPCclusterMI *nearestCluster = NULL;
1072
1073     // Inner Sector
1074     if(bInnerSector){
1075
1076       // Propagate to center of pad row and extract parameters
1077       AliTrackerBase::PropagateTrackTo(track,roc->GetPadRowRadii(0,iRow),kMass,5,kFALSE,kMaxSnp);
1078       xCur   = track->GetX();
1079       yCur   = track->GetY();
1080       zCur   = track->GetZ();
1081       if ( !fIdealTracking ) {
1082         zCur = zCur/vDrift + fTime0; // Look at time, not at z!
1083       }
1084       secCur = GetSector(track);
1085
1086       // Find the nearest cluster (TODO: correct road settings!)
1087       Printf("inner tracking here: x = %.2f, y = %.2f, z = %.2f (Row %d Sector %d)",xCur,yCur,zCur,iRow,secCur);
1088       nearestCluster = fInnerSectorArray[secCur%fkNSectorInner][iRow].FindNearest2(yCur,zCur,roady,roadz,indexCur);
1089
1090       // Move to next row if now cluster found
1091       if(!nearestCluster) continue;
1092       //Printf("Nearest Clusters = %d (of %d) ",indexCur,fInnerSectorArray[secCur%fkNSectorInner][iRow].GetN());
1093
1094     }
1095
1096     // Outer sector
1097     else{
1098
1099       // Propagate to center of pad row and extract parameters
1100       AliTrackerBase::PropagateTrackTo(track,roc->GetPadRowRadii(36,iRow-kIRowsTPC-1),kMass,5,kFALSE,kMaxSnp);
1101       xCur   = track->GetX();
1102       yCur   = track->GetY();
1103       zCur   = track->GetZ();
1104       if ( !fIdealTracking ) {
1105         zCur = zCur/vDrift + fTime0; // Look at time, not at z!
1106       }
1107       secCur = GetSector(track);
1108
1109       // Find the nearest cluster (TODO: correct road settings!)
1110       Printf("outer tracking here: x = %.2f, y = %.2f, z = %.2f (Row %d Sector %d)",xCur,yCur,zCur,iRow,secCur);
1111       nearestCluster = fOuterSectorArray[(secCur-fkNSectorInner*2)%fkNSectorOuter][iRow-kIRowsTPC-1].FindNearest2(yCur,zCur,roady,roadz,indexCur);
1112
1113       // Move to next row if now cluster found
1114       if(!nearestCluster) continue;
1115       //Printf("Nearest Clusters = %d (of %d)",indexCur,fOuterSectorArray[(secCur-fkNSectorInner*2)%fkNSectorOuter][iRow-kIRowsTPC-1].GetN());
1116
1117     }
1118
1119     // create track point from cluster
1120     SetTrackPointFromCluster(nearestCluster,nearestPoint);
1121
1122     //Printf("Track point = %.2f %.2f %.2f",nearestPoint.GetX(),nearestPoint.GetY(),nearestPoint.GetZ());
1123
1124     // correction
1125     // TODO: also correction when looking for the next cluster?
1126     if (fCorrectionType != kNoCorrection){
1127       Float_t xyz[3]={0,0,0};
1128       nearestPoint.GetXYZ(xyz);
1129       fTPCCorrection->CorrectPoint(xyz, nearestCluster->GetDetector());
1130       nearestPoint.SetXYZ(xyz);
1131     }
1132
1133     // rotate the cluster to the local detector frame
1134     track->Rotate(((nearestCluster->GetDetector()%18)*20+10)*TMath::DegToRad());
1135     AliTrackPoint prot = nearestPoint.Rotate(track->GetAlpha());   // rotate to the local frame - non distoted  point
1136     if (TMath::Abs(prot.GetX())<kRTPC0) continue;
1137     if (TMath::Abs(prot.GetX())>kRTPC1) continue;
1138
1139
1140     //Printf("Rotated Track point = %.2f %.2f %.2f",prot.GetX(),prot.GetY(),prot.GetZ());
1141
1142     // update track with the nearest track point
1143     Bool_t res=AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
1144
1145     if (!res) break;
1146
1147     if (TMath::Abs(track->GetZ())>kMaxZ) break;
1148     if (TMath::Abs(track->GetX())>kMaxR) break;
1149     //if (TMath::Abs(track->GetZ())<kZcut)continue;
1150
1151       Double_t pointPos[2]={0,0};
1152       Double_t pointCov[3]={0,0,0};
1153       pointPos[0]=prot.GetY();//local y
1154       pointPos[1]=prot.GetZ();//local z
1155       pointCov[0]=prot.GetCov()[3];//simay^2
1156       pointCov[1]=prot.GetCov()[4];//sigmayz
1157       pointCov[2]=prot.GetCov()[5];//sigmaz^2
1158
1159       if (!track->Update(pointPos,pointCov)) {printf("no update\n");}
1160
1161       ++nClus;
1162   }
1163
1164
1165   // propagation to refX
1166   AliTrackerBase::PropagateTrackTo(track,refX,kMass,5.,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
1167
1168   // rotate fittet track to the frame of the original track and propagate to same reference
1169   track->Rotate(tr->GetAlpha());
1170
1171   AliTrackerBase::PropagateTrackTo(track,refX,kMass,1.,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
1172
1173   Printf("We have %d clusters in this track!",nClus);
1174
1175   return track;
1176 }
1177
1178 //____________________________________________________________________________________
1179 AliExternalTrackParam* AliToyMCReconstruction::ClusterToTrackAssociation(const AliTPCseed *seed, Int_t trackID, Int_t &nClus)
1180 {
1181   //
1182   // Cluster to track association for given seed on an array of clusters
1183   //
1184
1185   // create track
1186   AliExternalTrackParam *track = new AliExternalTrackParam(*seed);
1187
1188   const AliTPCROC * roc = AliTPCROC::Instance();
1189
1190   const Double_t kRTPC0    = roc->GetPadRowRadii(0,0);
1191   const Double_t kRTPC1    = roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
1192   const Int_t kNRowsTPC    = roc->GetNRows(0) + roc->GetNRows(36) - 1;
1193   const Int_t kIRowsTPC    = roc->GetNRows(0) - 1;
1194   const Double_t kMaxSnp   = 0.85;
1195   const Double_t kMaxR     = 500.;
1196   const Double_t kMaxZ     = 500.;
1197   const Double_t roady     = 0.1;
1198   const Double_t roadz     = 0.01;
1199
1200   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1201
1202   Int_t  secCur, secOld   = -1;
1203   UInt_t indexCur = 0;
1204   Double_t xCur, yCur, zCur = 0.;
1205
1206 //   Float_t vDrift = GetVDrift();
1207
1208   // first propagate seed to outermost row
1209   Bool_t res0=AliTrackerBase::PropagateTrackTo(track,kRTPC1,kMass,5,kFALSE,kMaxSnp);
1210
1211   // Loop over rows and find the cluster candidates
1212   for( Int_t iRow = kNRowsTPC; iRow >= 0; --iRow ){
1213
1214     // inner or outer sector
1215     Bool_t bInnerSector = kTRUE;
1216     if(iRow > kIRowsTPC) bInnerSector = kFALSE;
1217
1218     // nearest track point/cluster (to be found)
1219     AliTrackPoint nearestPoint;
1220     AliTPCclusterMI *nearestCluster = NULL;
1221
1222     // Inner Sector
1223     if(bInnerSector){
1224
1225       // Propagate to center of pad row and extract parameters
1226       AliTrackerBase::PropagateTrackTo(track,roc->GetPadRowRadii(0,iRow),kMass,5,kFALSE,kMaxSnp);
1227       xCur   = track->GetX();
1228       yCur   = track->GetY();
1229       zCur   = track->GetZ();
1230       secCur = GetSector(track);
1231
1232       // Find the nearest cluster (TODO: correct road settings!)
1233       //Printf("inner tracking here: x = %.2f, y = %.2f, z = %.6f (Row %d Sector %d)",xCur,yCur,zCur,iRow,secCur);
1234       nearestCluster = fInnerSectorArray[secCur%fkNSectorInner][iRow].FindNearest2(yCur,zCur,roady,roadz,indexCur);
1235
1236       // Look again at -y if nothing found here and the sector changed with respect to the last found cluster
1237       // Increase also the road in this case
1238       if(!nearestCluster && secCur != secOld && secOld > -1){
1239         //Printf("inner tracking here 2: x = %.2f, y = %.2f, z = %.6f (Row %d Sector %d)",xCur,-yCur,zCur,iRow,secCur);
1240         nearestCluster = fInnerSectorArray[secCur%fkNSectorInner][iRow].FindNearest2(-yCur,zCur,roady*100,roadz,indexCur);
1241       }
1242
1243       // Move to next row if now cluster found
1244       if(!nearestCluster) continue;
1245       //Printf("Nearest Clusters = %d (of %d) ",indexCur,fInnerSectorArray[secCur%fkNSectorInner][iRow].GetN());
1246
1247     }
1248
1249     // Outer sector
1250     else{
1251
1252       // Propagate to center of pad row and extract parameters
1253       AliTrackerBase::PropagateTrackTo(track,roc->GetPadRowRadii(36,iRow-kIRowsTPC-1),kMass,5,kFALSE,kMaxSnp);
1254       xCur   = track->GetX();
1255       yCur   = track->GetY();
1256       zCur   = track->GetZ();
1257       secCur = GetSector(track);
1258
1259       // Find the nearest cluster (TODO: correct road settings!)
1260       Printf("res0 = %d, outer tracking here: x = %.2f, y = %.2f, z = %.6f (Row %d Sector %d)",res0,xCur,yCur,zCur,iRow,secCur);
1261       nearestCluster = fOuterSectorArray[(secCur-fkNSectorInner*2)%fkNSectorOuter][iRow-kIRowsTPC-1].FindNearest2(yCur,zCur,roady,roadz,indexCur);
1262
1263       // Look again at -y if nothing found here and the sector changed with respect to the last found cluster
1264       // Increase also the road in this case
1265       if(!nearestCluster && secCur != secOld && secOld > -1){
1266         Printf("outer tracking here 2: x = %.2f, y = %.2f, z = %.6f (Row %d Sector %d)",xCur,-yCur,zCur,iRow,secCur);
1267         nearestCluster = fOuterSectorArray[(secCur-fkNSectorInner*2)%fkNSectorOuter][iRow-kIRowsTPC-1].FindNearest2(-yCur,zCur,roady*100,roadz,indexCur);
1268       }
1269
1270
1271       // Move to next row if now cluster found
1272       if(!nearestCluster) continue;
1273       Printf("Nearest Clusters = %d (of %d)",indexCur,fOuterSectorArray[(secCur-fkNSectorInner*2)%fkNSectorOuter][iRow-kIRowsTPC-1].GetN());
1274
1275     }
1276
1277     // create track point from cluster
1278     SetTrackPointFromCluster(nearestCluster,nearestPoint);
1279
1280     //Printf("Track point = %.2f %.2f %.2f",nearestPoint.GetX(),nearestPoint.GetY(),nearestPoint.GetZ());
1281
1282     // rotate the cluster to the local detector frame
1283     track->Rotate(((nearestCluster->GetDetector()%18)*20+10)*TMath::DegToRad());
1284     AliTrackPoint prot = nearestPoint.Rotate(track->GetAlpha());   // rotate to the local frame - non distoted  point
1285     if (TMath::Abs(prot.GetX())<kRTPC0) continue;
1286     if (TMath::Abs(prot.GetX())>kRTPC1) continue;
1287
1288
1289     //Printf("Rotated Track point = %.2f %.2f %.2f",prot.GetX(),prot.GetY(),prot.GetZ());
1290
1291     // update track with the nearest track point
1292     Bool_t res=AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
1293
1294     if (!res) break;
1295
1296     if (TMath::Abs(track->GetZ())>kMaxZ) break;
1297     if (TMath::Abs(track->GetX())>kMaxR) break;
1298     //if (TMath::Abs(track->GetZ())<kZcut)continue;
1299
1300       Double_t pointPos[2]={0,0};
1301       Double_t pointCov[3]={0,0,0};
1302       pointPos[0]=prot.GetY();//local y
1303       pointPos[1]=prot.GetZ();//local z
1304       pointCov[0]=prot.GetCov()[3];//simay^2
1305       pointCov[1]=prot.GetCov()[4];//sigmayz
1306       pointCov[2]=prot.GetCov()[5];//sigmaz^2
1307
1308       if (!track->Update(pointPos,pointCov)) {printf("no update\n");}
1309       secOld = secCur;
1310
1311       //Printf("Cluster belongs to track = %d",nearestCluster->GetLabel(0));
1312
1313       // only count as associate cluster if it belongs to correct track!
1314       if(nearestCluster->GetLabel(0) == trackID)
1315         ++nClus;
1316   }
1317
1318   Printf("We have %d clusters in this track!",nClus);
1319
1320   return track;
1321 }
1322
1323 //____________________________________________________________________________________
1324 void AliToyMCReconstruction::AssociateClusters(AliTPCseed &seed, Int_t firstRow, Int_t lastRow, Bool_t direction)
1325 {
1326   //
1327   // do cluster to track association from first to last row
1328   // direction 0: outwards; 1: inwards
1329   //
1330
1331   Double_t roady=10.;
1332   Double_t roadz=10.;
1333
1334   AliRieman rieman1(160);
1335   AliRieman rieman2(160);
1336   SetRieman(seed,rieman1);
1337   CopyRieman(rieman1,rieman2);
1338
1339   Int_t sec=seed.GetSector();
1340   Int_t noLastPoint=0;
1341   //TODO: change to inward and outwar search?
1342   //      -> better handling of non consecutive points
1343   if (direction){
1344     firstRow*=-1;
1345     lastRow*=-1;
1346   }
1347
1348   //always from inside out
1349   if (firstRow>lastRow){
1350     Int_t tmp=firstRow;
1351     firstRow=lastRow;
1352     lastRow=tmp;
1353   }
1354
1355   for (Int_t row=firstRow; row<=lastRow && noLastPoint<3;++row) {
1356     Int_t iRow=TMath::Abs(row);
1357     const AliTPCclusterMI *cl=seed.GetClusterPointer(iRow);
1358     if (cl) continue;
1359
1360     const Int_t secrow = iRow<63?iRow:iRow-63;
1361
1362     AliTPCtrackerSector *arrSec=(iRow<63)?fInnerSectorArray:fOuterSectorArray;
1363     const AliTPCtrackerRow& kr  = arrSec[sec%36][secrow];
1364     const Double_t maxy=arrSec[sec%36].GetMaxY(secrow);
1365
1366     Double_t y=rieman1.GetYat(kr.GetX());
1367     Double_t z=rieman1.GetZat(kr.GetX());
1368
1369     if (TMath::Abs(y)>maxy) {
1370       AliError("Tracking over sector boundaries not implemented, yet");
1371       continue;
1372     }
1373
1374     AliTPCclusterMI *n=kr.FindNearest(y,z,roady,roadz);
1375     if (!n || n->IsUsed()) {
1376       ++noLastPoint;
1377       continue;
1378     }
1379     // check for quality of the cluster
1380     // TODO: better?
1381     rieman2.AddPoint(n->GetX(), n->GetY(), n->GetZ(),
1382                      TMath::Sqrt(n->GetSigmaY2()), TMath::Sqrt(n->GetSigmaZ2()));
1383     rieman2.Update();
1384 //     printf("      Riemann results: row=%d valid=%d, Chi2=%.2f (%.2f) %d (%d)",
1385 //            iRow, rieman2.IsValid(), rieman2.GetChi2(), rieman1.GetChi2(), n->GetLabel(0),seed.GetLabel());
1386     Double_t limit=2*rieman1.GetChi2();
1387     if (fClusterType==0) limit=1000;
1388     if (rieman2.GetChi2()>limit) {
1389       CopyRieman(rieman1,rieman2);
1390       ++noLastPoint;
1391 //       printf("\n");
1392       continue;
1393     }
1394 //     printf("  +++ \n");
1395
1396     noLastPoint=0;
1397     //use point
1398     rieman1.AddPoint(n->GetX(), n->GetY(), n->GetZ(),
1399                      TMath::Sqrt(n->GetSigmaY2()), TMath::Sqrt(n->GetSigmaZ2()));
1400     rieman1.Update();
1401
1402     seed.SetClusterPointer(iRow,n);
1403     //     if (iRow<seed.GetSeed1()) seed.SetSeed1(iRow);
1404     //     if (iRow>seed.GetSeed2()) seed.SetSeed2(iRow);
1405     n->Use();
1406
1407   }
1408 }
1409
1410 //____________________________________________________________________________________
1411 void AliToyMCReconstruction::ClusterToTrackAssociation(AliTPCseed &seed)
1412 {
1413   //
1414   //
1415   //
1416
1417 //   printf("\n ============ \nnext Seed: %d\n",seed.GetLabel());
1418   //assume seed is within one sector
1419   Int_t iMiddle=(seed.GetSeed1()+seed.GetSeed2())/2;
1420   //outward
1421   AssociateClusters(seed,iMiddle+1,158,kFALSE);
1422   //inward
1423   AssociateClusters(seed,0,iMiddle,kTRUE);
1424   seed.SetIsSeeding(kFALSE);
1425
1426   CookLabel(&seed,.6);
1427 }
1428
1429
1430 //____________________________________________________________________________________
1431 void AliToyMCReconstruction::InitSpaceCharge()
1432 {
1433   //
1434   // Init the space charge map
1435   //
1436
1437   TString filename="$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps5_50kHz_precal.root";
1438   if (fTree) {
1439     TList *l=fTree->GetUserInfo();
1440     for (Int_t i=0; i<l->GetEntries(); ++i) {
1441       TString s(l->At(i)->GetName());
1442       if (s.Contains("SC_")) {
1443         filename=s;
1444         break;
1445       }
1446     }
1447   }
1448
1449   printf("Initialising the space charge map using the file: '%s'\n",filename.Data());
1450   TFile f(filename.Data());
1451   fTPCCorrection=(AliTPCCorrection*)f.Get("map");
1452
1453   //   fTPCCorrection = new AliTPCSpaceCharge3D();
1454   //   fTPCCorrection->SetSCDataFileName("$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps10_50kHz.root");
1455   //   fTPCCorrection->SetOmegaTauT1T2(0.325,1,1); // Ne CO2
1456   // //   fTPCCorrection->SetOmegaTauT1T2(0.41,1,1.05); // Ar CO2
1457   //   fTPCCorrection->InitSpaceCharge3DDistortion();
1458
1459 }
1460
1461 //____________________________________________________________________________________
1462 Double_t AliToyMCReconstruction::GetVDrift() const
1463 {
1464   //
1465   //
1466   //
1467   return fTPCParam->GetDriftV();
1468 }
1469
1470 //____________________________________________________________________________________
1471 Double_t AliToyMCReconstruction::GetZLength(Int_t roc) const
1472 {
1473   //
1474   //
1475   //
1476   if (roc<0 || roc>71) return -1;
1477   return fTPCParam->GetZLength(roc);
1478 }
1479
1480 //____________________________________________________________________________________
1481 TTree* AliToyMCReconstruction::ConnectTrees (const char* files) {
1482   TString s=gSystem->GetFromPipe(Form("ls %s",files));
1483
1484   TTree *tFirst=0x0;
1485   TObjArray *arrFiles=s.Tokenize("\n");
1486
1487   for (Int_t ifile=0; ifile<arrFiles->GetEntriesFast(); ++ifile){
1488     TString name(arrFiles->At(ifile)->GetName());
1489
1490     TPRegexp reg(".*([0-9]_[0-9]_[0-9]_[0-9]{3}_[0-9]{2}).debug.root");
1491     TObjArray *arrMatch=0x0;
1492     arrMatch=reg.MatchS(name);
1493     TString matchName;
1494     if (arrMatch && arrMatch->At(1)) matchName=arrMatch->At(1)->GetName();
1495     else matchName=Form("%02d",ifile);
1496     delete arrMatch;
1497
1498     if (!tFirst) {
1499       TFile *f=TFile::Open(name.Data());
1500       if (!f) continue;
1501       TTree *t=(TTree*)f->Get("Tracks");
1502       if (!t) {
1503         delete f;
1504         continue;
1505       }
1506
1507       t->SetName(matchName.Data());
1508       tFirst=t;
1509     } else {
1510       tFirst->AddFriend(Form("t%s=Tracks",matchName.Data()), name.Data());
1511 //       tFirst->AddFriend(Form("t%d=Tracks",ifile), name.Data());
1512     }
1513   }
1514
1515   tFirst->GetListOfFriends()->Print();
1516   return tFirst;
1517 }
1518
1519 //____________________________________________________________________________________
1520 Int_t AliToyMCReconstruction::LoadOuterSectors() {
1521   //-----------------------------------------------------------------
1522   // This function fills outer TPC sectors with clusters.
1523   // Copy and paste from AliTPCtracker
1524   //-----------------------------------------------------------------
1525   Int_t nrows = fOuterSectorArray->GetNRows();
1526   UInt_t index=0;
1527   for (Int_t sec = 0;sec<fkNSectorOuter;sec++)
1528     for (Int_t row = 0;row<nrows;row++){
1529       AliTPCtrackerRow*  tpcrow = &(fOuterSectorArray[sec%fkNSectorOuter][row]);
1530       Int_t sec2 = sec+2*fkNSectorInner;
1531       //left
1532       Int_t ncl = tpcrow->GetN1();
1533       while (ncl--) {
1534         AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1535         index=(((sec2<<8)+row)<<16)+ncl;
1536         tpcrow->InsertCluster(c,index);
1537       }
1538       //right
1539       ncl = tpcrow->GetN2();
1540       while (ncl--) {
1541         AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1542         index=((((sec2+fkNSectorOuter)<<8)+row)<<16)+ncl;
1543         tpcrow->InsertCluster(c,index);
1544       }
1545       //
1546       // write indexes for fast acces
1547       //
1548       for (Int_t i=0;i<510;i++)
1549         tpcrow->SetFastCluster(i,-1);
1550       for (Int_t i=0;i<tpcrow->GetN();i++){
1551         Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1552         tpcrow->SetFastCluster(zi,i);  // write index
1553       }
1554       Int_t last = 0;
1555       for (Int_t i=0;i<510;i++){
1556         if (tpcrow->GetFastCluster(i)<0)
1557           tpcrow->SetFastCluster(i,last);
1558         else
1559           last = tpcrow->GetFastCluster(i);
1560       }
1561     }
1562   return 0;
1563 }
1564
1565
1566 //____________________________________________________________________________________
1567 Int_t  AliToyMCReconstruction::LoadInnerSectors() {
1568   //-----------------------------------------------------------------
1569   // This function fills inner TPC sectors with clusters.
1570   // Copy and paste from AliTPCtracker
1571   //-----------------------------------------------------------------
1572   Int_t nrows = fInnerSectorArray->GetNRows();
1573   UInt_t index=0;
1574   for (Int_t sec = 0;sec<fkNSectorInner;sec++)
1575     for (Int_t row = 0;row<nrows;row++){
1576       AliTPCtrackerRow*  tpcrow = &(fInnerSectorArray[sec%fkNSectorInner][row]);
1577       //
1578       //left
1579       Int_t ncl = tpcrow->GetN1();
1580       while (ncl--) {
1581         AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1582         index=(((sec<<8)+row)<<16)+ncl;
1583         tpcrow->InsertCluster(c,index);
1584       }
1585       //right
1586       ncl = tpcrow->GetN2();
1587       while (ncl--) {
1588         AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1589         index=((((sec+fkNSectorInner)<<8)+row)<<16)+ncl;
1590         tpcrow->InsertCluster(c,index);
1591       }
1592       //
1593       // write indexes for fast acces
1594       //
1595       for (Int_t i=0;i<510;i++)
1596         tpcrow->SetFastCluster(i,-1);
1597       for (Int_t i=0;i<tpcrow->GetN();i++){
1598         Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1599         tpcrow->SetFastCluster(zi,i);  // write index
1600       }
1601       Int_t last = 0;
1602       for (Int_t i=0;i<510;i++){
1603         if (tpcrow->GetFastCluster(i)<0)
1604           tpcrow->SetFastCluster(i,last);
1605         else
1606           last = tpcrow->GetFastCluster(i);
1607       }
1608
1609     }
1610   return 0;
1611 }
1612
1613 //____________________________________________________________________________________
1614 Int_t  AliToyMCReconstruction::GetSector(AliExternalTrackParam *track) {
1615   //-----------------------------------------------------------------
1616   // This function returns the sector number for a given track
1617   //-----------------------------------------------------------------
1618
1619   Int_t sector = -1;
1620
1621   // get the sector number
1622   // rotate point to global system (track is already global!)
1623   Double_t xd[3];
1624   track->GetXYZ(xd);
1625   //track->Local2GlobalPosition(xd,track->GetAlpha());
1626
1627   // use TPCParams to get the sector number
1628   Float_t xyz[3] = {xd[0],xd[1],xd[2]};
1629   Int_t   i[3]   = {0,0,0};
1630   if(fTPCParam){
1631     sector  = fTPCParam->Transform0to1(xyz,i);
1632   }
1633
1634   return sector;
1635 }
1636
1637 //____________________________________________________________________________________
1638 void  AliToyMCReconstruction::FillSectorStructure(Int_t maxev) {
1639   //-----------------------------------------------------------------
1640   // This function fills the sector structure of AliToyMCReconstruction
1641   //-----------------------------------------------------------------
1642
1643   // cluster array of all sectors
1644   fInnerSectorArray = new AliTPCtrackerSector[fkNSectorInner];
1645   fOuterSectorArray = new AliTPCtrackerSector[fkNSectorOuter];
1646
1647   for (Int_t i=0; i<fkNSectorInner; ++i) fInnerSectorArray[i].Setup(fTPCParam,0);
1648   for (Int_t i=0; i<fkNSectorOuter; ++i) fOuterSectorArray[i].Setup(fTPCParam,1);
1649
1650   Int_t count[72][96] = { {0} , {0} };
1651
1652   for (Int_t iev=0; iev<maxev; ++iev){
1653     printf("==============  Fill Clusters: Processing Event %6d  =================\n",iev);
1654     fTree->GetEvent(iev);
1655     for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
1656       printf(" > ======  Fill Clusters: Processing Track %6d ========  \n",itr);
1657       const AliToyMCTrack *tr=fEvent->GetTrack(itr);
1658
1659       // number of clusters to loop over
1660       const Int_t ncls=(fClusterType==0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
1661
1662       for(Int_t icl=0; icl<ncls; ++icl){
1663
1664         AliTPCclusterMI *cl=const_cast<AliTPCclusterMI *>(tr->GetSpacePoint(icl));
1665         if (fClusterType==1) cl=const_cast<AliTPCclusterMI *>(tr->GetDistortedSpacePoint(icl));
1666         if (!cl) continue;
1667
1668         Int_t sec = cl->GetDetector();
1669         Int_t row = cl->GetRow();
1670
1671         // set Q of the cluster to 1, Q==0 does not work for the seeding
1672         cl->SetQ(1);
1673
1674         // set cluster time to cluster Z (if not ideal tracking)
1675         if ( !fIdealTracking ) {
1676           // a 'valid' position in z is needed for the seeding procedure
1677 //           cl->SetZ(cl->GetTimeBin()*GetVDrift());
1678           cl->SetZ(cl->GetTimeBin());
1679         }
1680         //Printf("Fill clusters (sector %d row %d): %.2f %.2f %.2f %.2f",sec,row,cl->GetX(),cl->GetY(),cl->GetZ(),cl->GetTimeBin());
1681
1682         // fill arrays for inner and outer sectors (A/C side handled internally)
1683         if (sec<fkNSectorInner*2){
1684           fInnerSectorArray[sec%fkNSectorInner].InsertCluster(cl, count[sec][row], fTPCParam);
1685         }
1686         else{
1687           fOuterSectorArray[(sec-fkNSectorInner*2)%fkNSectorOuter].InsertCluster(cl, count[sec][row], fTPCParam);
1688         }
1689
1690         ++count[sec][row];
1691       }
1692     }
1693   }
1694
1695   // fill the arrays completely
1696   // LoadOuterSectors();
1697   // LoadInnerSectors();
1698
1699   // // check the arrays
1700   // for (Int_t i=0; i<fkNSectorInner; ++i){
1701   //   for (Int_t j=0; j<fInnerSectorArray[i].GetNRows(); ++j){
1702   //     if(fInnerSectorArray[i][j].GetN()>0){
1703   //    Printf("Inner: Sector %d Row %d : %d",i,j,fInnerSectorArray[i][j].GetN());
1704   //     }
1705   //   }
1706   // }
1707   // for (Int_t i=0; i<fkNSectorInner; ++i){
1708   //   for (Int_t j=0; j<fOuterSectorArray[i].GetNRows(); ++j){
1709   //     if(fOuterSectorArray[i][j].GetN()>0){
1710   //    Printf("Outer: Sector %d Row %d : %d",i,j,fOuterSectorArray[i][j].GetN());
1711   //     }
1712   //   }
1713   // }
1714 }
1715
1716 //____________________________________________________________________________________
1717 void  AliToyMCReconstruction::FillSectorStructureAC() {
1718   //-----------------------------------------------------------------
1719   // This function fills the sector structure of AliToyMCReconstruction
1720   //-----------------------------------------------------------------
1721
1722   /*
1723    my god is the AliTPCtrackerSector stuff complicated!!!
1724    Ok, so here we will not fill the fClusters1 and fClusters2 of AliTPCtrackerRow,
1725    using InsertCluster of AliTPCtrackerSector, but only the fClusters via InsertCluster
1726    of AliTPCtrackerRow itself which then will not be owner, but we create an array in
1727    here (fAllClusters) which owns all clusters ...
1728   */
1729
1730   fIsAC=kTRUE;
1731   // cluster array of all sectors
1732   fInnerSectorArray = new AliTPCtrackerSector[2*fkNSectorInner];
1733   fOuterSectorArray = new AliTPCtrackerSector[2*fkNSectorOuter];
1734
1735   for (Int_t i=0; i<2*fkNSectorInner; ++i) {
1736     fInnerSectorArray[i].Setup(fTPCParam,0);
1737   }
1738
1739   for (Int_t i=0; i<2*fkNSectorOuter; ++i) {
1740     fOuterSectorArray[i].Setup(fTPCParam,1);
1741   }
1742
1743   Int_t count[72][96] = { {0} , {0} };
1744
1745   for (Int_t iev=0; iev<fNmaxEvents; ++iev){
1746     printf("==============  Fill Clusters: Processing Event %6d  =================\n",iev);
1747     fTree->GetEvent(iev);
1748     for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
1749 //       printf(" > ======  Fill Clusters: Processing Track %6d ========  \n",itr);
1750       const AliToyMCTrack *tr=fEvent->GetTrack(itr);
1751
1752       // number of clusters to loop over
1753       const Int_t ncls=(fClusterType==0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
1754
1755       // check if expansion of the cluster arrays is needed.
1756       if (fAllClusters.GetEntriesFast()+ncls>=fAllClusters.Capacity()) fAllClusters.Expand(2*fAllClusters.Capacity());
1757       for(Int_t icl=0; icl<ncls; ++icl){
1758
1759         AliTPCclusterMI *cl=const_cast<AliTPCclusterMI *>(tr->GetSpacePoint(icl));
1760         if (fClusterType==1) cl=const_cast<AliTPCclusterMI *>(tr->GetDistortedSpacePoint(icl));
1761         if (!cl) continue;
1762
1763         // register copy to the cluster array
1764         cl = new(fAllClusters[fAllClusters.GetEntriesFast()]) AliTPCclusterMI(*cl);
1765
1766         Int_t sec = cl->GetDetector();
1767         Int_t row = cl->GetRow();
1768
1769         // set Q of the cluster to 1, Q==0 does not work for the seeding
1770         cl->SetQ(1);
1771
1772         // set cluster time to cluster Z (if not ideal tracking)
1773         if ( !fIdealTracking ) {
1774           // a 'valid' position in z is needed for the seeding procedure
1775           Double_t sign=1;
1776           if (((sec/18)%2)==1) sign=-1;
1777           cl->SetZ(cl->GetTimeBin()*GetVDrift()*sign);
1778           //mark cluster to be time*vDrift by setting the type to 1
1779           cl->SetType(1);
1780 //           cl->SetZ(cl->GetTimeBin());
1781         }
1782         //Printf("Fill clusters (sector %d row %d): %.2f %.2f %.2f %.2f",sec,row,cl->GetX(),cl->GetY(),cl->GetZ(),cl->GetTimeBin());
1783
1784         // fill arrays for inner and outer sectors (A/C side handled internally)
1785         if (sec<fkNSectorInner*2){
1786           fInnerSectorArray[sec][row].InsertCluster(cl, count[sec][row]);
1787         }
1788         else{
1789           fOuterSectorArray[(sec-fkNSectorInner*2)][row].InsertCluster(cl, count[sec][row]);
1790         }
1791
1792         ++count[sec][row];
1793       }
1794     }
1795   }
1796
1797 }
1798
1799 //____________________________________________________________________________________
1800 AliToyMCTrack *AliToyMCReconstruction::ConvertTPCSeedToToyMCTrack(const AliTPCseed &seed)
1801 {
1802   //
1803   //
1804   //
1805
1806
1807   AliToyMCTrack *tToy = new AliToyMCTrack(seed);
1808
1809   for (Int_t icl=0; icl<159; ++icl){
1810     const AliTPCclusterMI *cl=seed.GetClusterFast(icl);
1811     if (!cl) continue;
1812     if (fClusterType==0){
1813       tToy->AddSpacePoint(*cl);
1814     } else {
1815       tToy->AddDistortedSpacePoint(*cl);
1816     }
1817   }
1818
1819   return tToy;
1820 }
1821
1822 //____________________________________________________________________________________
1823 AliExternalTrackParam* AliToyMCReconstruction::GetRefittedTrack(const AliTPCseed &seed)
1824 {
1825   //
1826   //
1827   //
1828
1829   AliExternalTrackParam *track=0x0;
1830
1831   const Float_t vDrift=GetVDrift();
1832   const Float_t zLength=GetZLength(0);
1833   const Double_t kMaxSnp = 0.85;
1834   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
1835
1836   AliToyMCTrack *toyTrack = ConvertTPCSeedToToyMCTrack(seed);
1837
1838   fTime0 = 0;
1839
1840   //get t0 estimate
1841   fCreateT0seed = kTRUE;
1842   AliExternalTrackParam *t0seed = GetSeedFromTrack(toyTrack,kTRUE);
1843   if (!t0seed) return 0x0;
1844
1845   fTime0 = t0seed->GetZ()-zLength/vDrift;
1846   delete t0seed;
1847   t0seed=0x0;
1848
1849   fCreateT0seed = kFALSE;
1850   AliExternalTrackParam *dummy = GetSeedFromTrack(toyTrack,kTRUE);
1851   track = GetFittedTrackFromSeed(toyTrack, dummy);
1852   delete dummy;
1853   // propagate seed to 0
1854   AliTrackerBase::PropagateTrackTo(track,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
1855
1856   return track;
1857 }
1858
1859 //____________________________________________________________________________________
1860 AliTPCclusterMI* AliToyMCReconstruction::FindMiddleCluster(const AliTPCclusterMI *clInner, const AliTPCclusterMI *clOuter,
1861                                                            Double_t roady, Double_t roadz,
1862                                                            AliRieman &seedFit)
1863 {
1864   //
1865   //
1866   //
1867
1868   const Int_t rocInner = clInner->GetDetector();
1869   const Int_t rocOuter = clOuter->GetDetector();
1870
1871   if ( (rocInner%18) != (rocOuter%18) ){
1872     AliError("Currently only same Sector implemented");
1873     return 0x0;
1874   }
1875
1876   const Int_t innerDet=clInner->GetDetector();
1877   const Int_t outerDet=clOuter->GetDetector();
1878   const Int_t globalRowInner  = clInner->GetRow() +(innerDet >35)*63;
1879   const Int_t globalRowOuter  = clOuter->GetRow() +(outerDet >35)*63;
1880
1881   AliTPCclusterMI *n=0x0;
1882
1883   // allow flexibility of +- nRowsGrace Rows to find a middle cluster
1884   const Int_t nRowsGrace = 0;
1885   for (Int_t iter=0; iter<2*nRowsGrace+1; ++iter){
1886     Int_t iMiddle  = (globalRowInner+globalRowOuter)/2;
1887     iMiddle+=((iter+1)/2)*(1-2*((iter+1)%2));
1888
1889     Int_t localRow=iMiddle;
1890     Int_t roc = innerDet;
1891     if (iMiddle>62){
1892       localRow-=63;
1893       roc=outerDet;
1894     }
1895
1896     AliTPCtrackerSector *arrRow=(iMiddle<63)?fInnerSectorArray:fOuterSectorArray;
1897     const AliTPCtrackerRow& krMiddle = arrRow[roc%36][localRow]; // middle
1898     // initial guess use simple linear interpolation
1899     Double_t y=(clInner->GetY()+clOuter->GetY())/2;
1900     Double_t z=(clInner->GetZ()+clOuter->GetZ())/2;
1901     if (seedFit.IsValid()) {
1902       // update values once the fit was performed
1903       y=seedFit.GetYat(krMiddle.GetX());
1904       z=seedFit.GetZat(krMiddle.GetX());
1905     }
1906
1907     n=krMiddle.FindNearest(y,z,roady,roadz);
1908     if (n) break;
1909   }
1910
1911 //   if (n)
1912 //     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,
1913 //          n->GetX(), n->GetY(),n->GetZ(),n->GetLabel(0),
1914 //            clInner->GetX(), clInner->GetY(),clInner->GetZ(),clInner->GetLabel(0),
1915 //            clOuter->GetX(), clOuter->GetY(),clOuter->GetZ(),clOuter->GetLabel(0)
1916 //         );
1917   return n;
1918 }
1919
1920 //____________________________________________________________________________________
1921 void AliToyMCReconstruction::AddMiddleClusters(AliTPCseed *seed,
1922                                                const AliTPCclusterMI *clInner, const AliTPCclusterMI *clOuter,
1923                                                Double_t roady, Double_t roadz,
1924                                                Int_t &nTotalClusters, AliRieman &seedFit)
1925 {
1926   //
1927   // Iteratively add the middle clusters
1928   //
1929
1930   // update rieman fit with every second point added
1931   AliTPCclusterMI *clMiddle=FindMiddleCluster(clInner,clOuter,roady,roadz,seedFit);
1932
1933   // break iterative process
1934   if (!clMiddle || clMiddle->IsUsed()) return;
1935
1936   const Int_t globalRowInner  = clInner->GetRow() +(clInner->GetDetector() >35)*63;
1937   const Int_t globalRowMiddle = clMiddle->GetRow()+(clMiddle->GetDetector()>35)*63;
1938   const Int_t globalRowOuter  = clOuter->GetRow() +(clOuter->GetDetector() >35)*63;
1939
1940   seed->SetClusterPointer(globalRowMiddle,clMiddle);
1941   ++nTotalClusters;
1942 //   printf("    > Total clusters: %d\n",nTotalClusters);
1943   seedFit.AddPoint(clMiddle->GetX(), clMiddle->GetY(), clMiddle->GetZ(),
1944                    TMath::Sqrt(clMiddle->GetSigmaY2()), TMath::Sqrt(clMiddle->GetSigmaZ2()));
1945
1946   if (seedFit.GetN()>3) {
1947 //     printf("      call update: %d (%d)\n",seedFit.GetN(),nTotalClusters);
1948 //     printf("      Riemann results: valid=%d, Chi2=%.2f, Chi2Y=%.2f, Chi2Z=%.2f -- %d\n",
1949 //            seedFit.IsValid(), seedFit.GetChi2(), seedFit.GetChi2Y(), seedFit.GetChi2Z(), clMiddle->GetLabel(0));
1950     seedFit.Update();
1951   }
1952   if ( seedFit.IsValid() && seedFit.GetChi2()>1000 ) return;
1953
1954   // Add middle clusters towards outer radius
1955   if (globalRowMiddle+1<globalRowOuter) AddMiddleClusters(seed,clMiddle,clOuter,roady,roadz,nTotalClusters,seedFit);
1956
1957   // Add middle clusters towards innter radius
1958   if (globalRowInner+1<globalRowMiddle) AddMiddleClusters(seed,clInner,clMiddle,roady,roadz,nTotalClusters,seedFit);
1959
1960   return;
1961 }
1962
1963 //____________________________________________________________________________________
1964 Int_t AliToyMCReconstruction::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t iRowInner, Int_t iRowOuter)
1965 {
1966   //
1967   // find seeds in a sector, requires iRowInner < iRowOuter
1968   // iRowXXX runs from 0-159, so over IROC and OROC
1969   //
1970
1971   if (!fIsAC) {
1972     AliError("This function requires the sector arrays filled for AC tracking");
1973     return 0;
1974   }
1975
1976   // swap rows in case they are in the wrong order
1977   if (iRowInner>iRowOuter) {
1978     Int_t tmp=iRowInner;
1979     iRowInner=iRowOuter;
1980     iRowOuter=tmp;
1981   }
1982
1983   if (iRowOuter>158) iRowOuter=158;
1984   if (iRowInner<0)   iRowInner=0;
1985
1986   // Define the search roads:
1987   // timeRoadCombinatorics - the maximum time difference used for the
1988   //    combinatorics. Since we don't have a 'z-Vertex' estimate this will
1989   //    reduce the combinatorics significantly when iterating on one TF
1990   //    use a little more than one full drift length of the TPC to allow for
1991   //    distortions
1992   //
1993   // timeRoad - the time difference allowed when associating the cluster
1994   //    in the middle of the other two used for the initial search
1995   //
1996   // padRoad  - the local y difference allowed when associating the middle cluster
1997
1998 //   Double_t timeRoadCombinatorics = 270./vDrift;
1999 //   Double_t timeRoad = 20./vDrift;
2000   Double_t timeRoadCombinatorics = 270.;
2001   Double_t timeRoad = 10.;
2002   Double_t  padRoad = 5.;
2003
2004
2005   // fOuterSectorArray runs from 0-95, so from iRowXXX which runs from 0-159
2006   //   the number of rows in the IROC has to be subtracted
2007   const Int_t innerRows=fInnerSectorArray->GetNRows();
2008
2009   AliTPCtrackerSector *arrInnerRow=(iRowInner<63)?fInnerSectorArray:fOuterSectorArray;
2010   AliTPCtrackerSector *arrOuterRow=(iRowOuter<63)?fInnerSectorArray:fOuterSectorArray;
2011
2012   const AliTPCtrackerRow& krInner  = arrInnerRow[sec][iRowInner - (iRowInner>62)*innerRows];   // up
2013   const AliTPCtrackerRow& krOuter  = arrOuterRow[sec][iRowOuter - (iRowOuter>62)*innerRows];   // down
2014
2015   AliTPCseed *seed = new AliTPCseed;
2016
2017   const Int_t nMaxClusters=iRowOuter-iRowInner+1;
2018 //   Int_t nScannedClusters = 0;
2019
2020   Int_t nseedAdded=0;
2021   // loop over all points in the firstand last search row
2022   for (Int_t iOuter=0; iOuter < krOuter; iOuter++) {
2023     const AliTPCclusterMI *clOuter = krOuter[iOuter];
2024     if (clOuter->IsUsed()) continue;
2025     for (Int_t iInner=0; iInner < krInner; iInner++) {
2026       const AliTPCclusterMI *clInner = krInner[iInner];
2027       if (clInner->IsUsed()) continue;
2028 // printf("\n\n Check combination %d (%d), %d (%d) -- %d (%d) -- %d\n",iOuter, iInner, clOuter->GetLabel(0), clInner->GetLabel(0),iRowOuter,iRowInner,sec);
2029       // check maximum distance for combinatorics
2030       if (TMath::Abs(clOuter->GetZ()-clInner->GetZ())>timeRoadCombinatorics) continue;
2031 // printf("  Is inside one drift\n");
2032
2033       // use rieman fit for seed description
2034       AliRieman seedFit(159);
2035       // add first two points
2036       seedFit.AddPoint(clInner->GetX(), clInner->GetY(), clInner->GetZ(),
2037                        TMath::Sqrt(clInner->GetSigmaY2()), TMath::Sqrt(clInner->GetSigmaZ2()));
2038       seedFit.AddPoint(clOuter->GetX(), clOuter->GetY(), clOuter->GetZ(),
2039                        TMath::Sqrt(clOuter->GetSigmaY2()), TMath::Sqrt(clOuter->GetSigmaZ2()));
2040
2041       // Iteratively add all clusters in the respective middle
2042       Int_t nFoundClusters=2;
2043       AddMiddleClusters(seed,clInner,clOuter,padRoad,timeRoad,nFoundClusters,seedFit);
2044 //       printf("  Clusters attached: %d\n",nFoundClusters);
2045       if (nFoundClusters>2) seedFit.Update();
2046 //       printf("  Riemann results: valid=%d, Chi2=%.2f, Chi2Y=%.2f, Chi2Z=%.2f\n",
2047 //              seedFit.IsValid(), seedFit.GetChi2(), seedFit.GetChi2Y(), seedFit.GetChi2Z());
2048
2049       // check for minimum number of assigned clusters and a decent chi2
2050       if ( nFoundClusters<0.5*nMaxClusters || seedFit.GetChi2()>1000 ){
2051         seed->Reset();
2052         continue;
2053       }
2054 //       printf("  >>> Seed accepted\n");
2055       // add original outer clusters
2056       const Int_t globalRowInner  = clInner->GetRow() +(clInner->GetDetector() >35)*63;
2057       const Int_t globalRowOuter  = clOuter->GetRow() +(clOuter->GetDetector() >35)*63;
2058
2059       seed->SetClusterPointer(globalRowInner,const_cast<AliTPCclusterMI*>(clInner));
2060       seed->SetClusterPointer(globalRowOuter,const_cast<AliTPCclusterMI*>(clOuter));
2061
2062       // set parameters retrieved from AliRieman
2063       Double_t params[5]={0};
2064       Double_t covar[15]={0};
2065       const Double_t alpha=TMath::DegToRad()*(clInner->GetDetector()%18*20.+10.);
2066       const Double_t x=clInner->GetX();
2067       seedFit.GetExternalParameters(x,params,covar);
2068
2069       seed->Set(x,alpha,params,covar);
2070
2071       // set label of the seed. At least 60% of the clusters need the correct label
2072       CookLabel(seed,.6);
2073 //       printf("  - Label: %d\n",seed->GetLabel());
2074       // mark clusters as being used
2075       MarkClustersUsed(seed);
2076
2077       seed->SetSeed1(iRowInner);
2078       seed->SetSeed2(iRowOuter);
2079       seed->SetSector(sec);
2080       ++nseedAdded;
2081
2082       seed->SetUniqueID(arr->GetEntriesFast());
2083       seed->SetIsSeeding(kTRUE);
2084
2085       arr->Add(seed);
2086       seed=new AliTPCseed;
2087
2088       break;
2089     }
2090   }
2091   //delete surplus seed
2092   delete seed;
2093   seed=0x0;
2094
2095   return nseedAdded;
2096 }
2097 //____________________________________________________________________________________
2098 void AliToyMCReconstruction::MakeSeeds(TObjArray * /*arr*/, Int_t sec, Int_t iRow1, Int_t iRow2)
2099 {
2100   //
2101   // Create seeds between i1 and i2 (stored in arr) with TLinearFitter
2102   //
2103   // sec: sector number
2104   // iRow1:  upper row
2105   // iRow2:  lower row
2106   //
2107
2108   // Create Fitter
2109   static TLinearFitter fitter(3,"pol2");
2110
2111   // Get 3 padrows (iRow1,iMiddle=(iRow1+iRow2)/2,iRow2)
2112   Int_t iMiddle = (iRow1+iRow2)/2;
2113   const AliTPCtrackerRow& krd = fOuterSectorArray[sec][iRow2-fInnerSectorArray->GetNRows()];   // down
2114   const AliTPCtrackerRow& krm = fOuterSectorArray[sec][iMiddle-fInnerSectorArray->GetNRows()]; // middle
2115   const AliTPCtrackerRow& kru = fOuterSectorArray[sec][iRow1-fInnerSectorArray->GetNRows()];   // up
2116
2117   // Loop over 3 cluster possibilities
2118   for (Int_t iu=0; iu < kru; iu++) {
2119     for (Int_t im=0; im < krm; im++) {
2120       for (Int_t id=0; id < krd; id++) {
2121
2122         // clear all points
2123         fitter.ClearPoints();
2124
2125         // add all three points to fitter
2126         Double_t xy[2] = {kru[iu]->GetX(),kru[iu]->GetY()};
2127         Double_t z     = kru[iu]->GetZ();
2128         fitter.AddPoint(xy,z);
2129
2130         xy[0] = krm[im]->GetX();
2131         xy[1] = krm[im]->GetY();
2132         z     = krm[im]->GetZ();
2133         fitter.AddPoint(xy,z);
2134
2135         xy[0] = krd[id]->GetX();
2136         xy[1] = krd[id]->GetY();
2137         z     = krd[id]->GetZ();
2138         fitter.AddPoint(xy,z);
2139
2140         // Evaluate and get parameters
2141         fitter.Eval();
2142
2143         // how to get the other clusters now?
2144         // ...
2145
2146       }
2147     }
2148   }
2149 }
2150
2151 //____________________________________________________________________________________
2152 void AliToyMCReconstruction::InitStreamer(TString addName, Int_t level)
2153 {
2154   //
2155   // initilise the debug streamer and set the logging level
2156   //   use a default naming convention
2157   //
2158
2159   delete fStreamer;
2160   fStreamer=0x0;
2161
2162   if (addName.IsNull()) addName=".dummy";
2163
2164   if (!fTree) return;
2165
2166   TString debugName=fInputFile->GetName();
2167   debugName.ReplaceAll(".root","");
2168   debugName.Append(Form(".%1d.%1d_%1d_%1d_%03d_%02d",
2169                         fUseMaterial,fIdealTracking,fClusterType,
2170                         Int_t(fCorrectionType),fSeedingRow,fSeedingDist));
2171   debugName.Append(addName);
2172   debugName.Append(".root");
2173
2174   gSystem->Exec(Form("test -f %s && rm %s", debugName.Data(), debugName.Data()));
2175   fStreamer=new TTreeSRedirector(debugName.Data());
2176   fStreamer->SetUniqueID(level);
2177
2178   gROOT->cd();
2179 }
2180
2181 //____________________________________________________________________________________
2182 void AliToyMCReconstruction::ConnectInputFile(const char* file, Int_t nmaxEv)
2183 {
2184   //
2185   // connect the tree and event pointer from the input file
2186   //
2187
2188   delete fInputFile;
2189   fInputFile=0x0;
2190   fEvent=0x0;
2191   fTree=0;
2192
2193   fInputFile=TFile::Open(file);
2194   if (!fInputFile || !fInputFile->IsOpen() || fInputFile->IsZombie()) {
2195     delete fInputFile;
2196     fInputFile=0x0;
2197     AliError(Form("ERROR: couldn't open the file '%s'\n", file));
2198     return;
2199   }
2200
2201   fTree=(TTree*)fInputFile->Get("toyMCtree");
2202   if (!fTree) {
2203     AliError(Form("ERROR: couldn't read the 'toyMCtree' from file '%s'\n", file));
2204     return;
2205   }
2206
2207   fEvent=0x0;
2208   fTree->SetBranchAddress("event",&fEvent);
2209
2210   gROOT->cd();
2211
2212   fNmaxEvents=fTree->GetEntries();
2213   if (nmaxEv>-1) fNmaxEvents=TMath::Min(nmaxEv,fNmaxEvents);
2214
2215   // setup space charge map from Userinfo of the tree
2216   InitSpaceCharge();
2217
2218   // setup the track maps
2219   SetupTrackMaps();
2220 }
2221
2222 //____________________________________________________________________________________
2223 void AliToyMCReconstruction::Cleanup()
2224 {
2225   //
2226   // Cleanup input data
2227   //
2228
2229   if (fStreamer) delete fStreamer;
2230   fStreamer=0x0;
2231
2232   delete fEvent;
2233   fEvent = 0x0;
2234
2235 //   delete fTree;
2236   fTree=0x0;
2237
2238   delete fInputFile;
2239   fInputFile=0x0;
2240
2241 }
2242
2243 //____________________________________________________________________________________
2244 void AliToyMCReconstruction::SetupTrackMaps()
2245 {
2246   //
2247   //
2248   //
2249
2250   fMapTrackEvent.Delete();
2251   fMapTrackTrackInEvent.Delete();
2252
2253   if (!fTree) {
2254     AliError("Tree not connected");
2255     return;
2256   }
2257
2258   Int_t nmaxEv=fTree->GetEntries();
2259   if (fNmaxEvents>-1) nmaxEv=fNmaxEvents;
2260
2261   for (Int_t iev=0; iev<nmaxEv; ++iev) {
2262     fTree->GetEvent(iev);
2263     if (!fEvent) continue;
2264
2265     const Int_t ntracks=fEvent->GetNumberOfTracks();
2266     if (fMapTrackEvent.GetSize()+ntracks>=fMapTrackEvent.Capacity()) fMapTrackEvent.Expand(2*fMapTrackEvent.Capacity());
2267     if (fMapTrackTrackInEvent.GetSize()+ntracks>=fMapTrackTrackInEvent.Capacity()) fMapTrackTrackInEvent.Expand(2*fMapTrackTrackInEvent.Capacity());
2268
2269     for (Int_t itrack=0; itrack<ntracks; ++itrack){
2270       Int_t label=fEvent->GetTrack(itrack)->GetUniqueID();
2271
2272       fMapTrackEvent.Add(label,iev);
2273       fMapTrackTrackInEvent.Add(label,itrack);
2274     }
2275   }
2276
2277 }
2278
2279 //____________________________________________________________________________________
2280 void AliToyMCReconstruction::CookLabel(AliTPCseed *seed, Double_t fraction, Int_t info[5])
2281 {
2282   //
2283   //
2284   //
2285
2286   Int_t labels[159]={0};
2287 //   Long64_t posMaxLabel=-1;
2288   Int_t nMaxLabel=0;  // clusters from maximum label
2289   Int_t nMaxLabel2=0; // clusters from second maximum
2290   Int_t nlabels=0;
2291   Int_t maxLabel=-1;  // label with most clusters
2292   Int_t maxLabel2=-1; // label with second most clusters
2293   Int_t nclusters=0;
2294   TExMap labelMap(159);
2295
2296   for (Int_t icl=0; icl<159; ++icl) {
2297     const AliTPCclusterMI *cl=seed->GetClusterPointer(icl);
2298     if (!cl) continue;
2299     ++nclusters;
2300
2301     const Int_t clLabel=cl->GetLabel(0);
2302     // a not assinged value returns 0, so we need to add 1 and subtract it afterwards
2303     Long64_t labelPos=labelMap.GetValue(clLabel);
2304
2305     if (!labelPos) {
2306       labelPos=nlabels+1;
2307       labelMap.Add(clLabel,labelPos);
2308       ++nlabels;
2309     }
2310     --labelPos;
2311
2312     const Int_t nCurrentLabel=++labels[labelPos];
2313     if (nCurrentLabel>nMaxLabel) {
2314       nMaxLabel2=nMaxLabel;
2315       nMaxLabel=nCurrentLabel;
2316 //       posMaxLabel=labelPos;
2317       maxLabel2=maxLabel;
2318       maxLabel=clLabel;
2319     }
2320   }
2321
2322   if (Double_t(nMaxLabel)/nclusters<fraction) maxLabel=-maxLabel;
2323
2324   seed->SetLabel(maxLabel);
2325
2326   if (info) {
2327     info[0]=nMaxLabel;
2328     info[1]=nMaxLabel2;
2329     info[2]=maxLabel2;
2330     info[3]=nclusters;
2331     info[4]=nlabels;
2332   }
2333 }
2334
2335
2336 //____________________________________________________________________________________
2337 void AliToyMCReconstruction::DumpSeedInfo(TObjArray *arr)
2338 {
2339
2340   // for debugging
2341   if (!fStreamer || !fTree) return;
2342   // swap rows in case they are in the wrong order
2343   AliInfo(Form("Number of Seeds: %d (out of %d tracks)",arr->GetEntriesFast(),fMapTrackEvent.GetSize()));
2344
2345   //loop over all events and tracks and try to associate the seed to the track
2346   for (Int_t iseed=0; iseed<arr->GetEntriesFast(); ++iseed){
2347     AliTPCseed *seed = static_cast<AliTPCseed*>(arr->UncheckedAt(iseed));
2348
2349     // get original track
2350     Int_t seedLabel=seed->GetLabel();
2351     Int_t iev=fMapTrackEvent.GetValue(TMath::Abs(seedLabel));
2352     Int_t itr=fMapTrackTrackInEvent.GetValue(TMath::Abs(seedLabel));
2353
2354     fTree->GetEvent(iev);
2355     const AliToyMCTrack *toyTrack = fEvent->GetTrack(itr);
2356
2357     DumpSeedInfo(toyTrack,seed);
2358   }
2359 }
2360
2361 //____________________________________________________________________________________
2362 void AliToyMCReconstruction::DumpTrackInfo(TObjArray *arr)
2363 {
2364
2365   // for debugging
2366   if (!fStreamer || !fTree) return;
2367   // swap rows in case they are in the wrong order
2368   AliInfo(Form("Number of Seeds: %d (out of %d tracks)",arr->GetEntriesFast(),fMapTrackEvent.GetSize()));
2369
2370   //loop over all events and tracks and try to associate the seed to the track
2371   AliTPCseed dummySeed;
2372   for (Int_t iev=0; iev<fNmaxEvents; ++iev) {
2373     fTree->GetEvent(iev);
2374     const Int_t ntracks=fEvent->GetNumberOfTracks();
2375     for (Int_t itr=0; itr<ntracks; ++itr) {
2376       const AliToyMCTrack *toyTrack = fEvent->GetTrack(itr);
2377
2378       Bool_t foundSeed=kFALSE;
2379       for (Int_t iseed=0; iseed<arr->GetEntriesFast(); ++iseed) {
2380         AliTPCseed *seed = static_cast<AliTPCseed*>(arr->UncheckedAt(iseed));
2381         const UInt_t tmpLabel=TMath::Abs(seed->GetLabel());
2382         if (toyTrack->GetUniqueID()!=tmpLabel) continue;
2383
2384         // dump all seeds with the correct label
2385         DumpSeedInfo(toyTrack,seed);
2386         foundSeed=kTRUE;
2387       }
2388
2389       if (!foundSeed) DumpSeedInfo(toyTrack,&dummySeed);
2390
2391     }
2392   }
2393 }
2394
2395 //____________________________________________________________________________________
2396 void AliToyMCReconstruction::DumpSeedInfo(const AliToyMCTrack *toyTrack, AliTPCseed *seed)
2397 {
2398   //
2399   //
2400   //
2401
2402   const Double_t kMaxSnp = 0.85;
2403   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
2404   Float_t vDrift=GetVDrift();
2405   Float_t zLength=GetZLength(0);
2406
2407   AliExternalTrackParam tOrig(*toyTrack);
2408   AliExternalTrackParam tOrigITS(*toyTrack);
2409
2410   // propagate original track to ITS last layer
2411   Double_t lastLayerITS = 43.0; // same as in AliToyMCEventGenerator::MakeITSClusters (hard coded)
2412   AliTrackerBase::PropagateTrackTo(&tOrigITS,lastLayerITS,kMass,1,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
2413
2414   AliExternalTrackParam dummyParam;
2415   Bool_t isDummy=kFALSE;
2416   //create refitted track, this will also set the fTime0
2417   AliExternalTrackParam *track=GetRefittedTrack(*seed);
2418   if (!track) {
2419     track=&dummyParam;
2420     isDummy=kTRUE;
2421   }
2422   AliTrackerBase::PropagateTrackTo(track,0,kMass,5.,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
2423   track->Rotate(tOrig.GetAlpha());
2424   AliTrackerBase::PropagateTrackTo(track,0,kMass,1.,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
2425
2426   // rotate fitted track to the frame of the original track and propagate to same reference
2427   AliExternalTrackParam trackITS(*track);
2428   AliTrackerBase::PropagateTrackTo(&trackITS,lastLayerITS,kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
2429   trackITS.Rotate(tOrigITS.GetAlpha());
2430   AliTrackerBase::PropagateTrackTo(&trackITS,lastLayerITS,kMass,1,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
2431
2432   Int_t seedSec=seed->GetSector();
2433   Int_t seedID =seed->GetUniqueID();
2434   //
2435   //count findable and found clusters in the seed
2436   //
2437   Int_t iRowInner=seed->GetSeed1();
2438   Int_t iRowOuter=seed->GetSeed2();
2439
2440   Int_t nClustersSeedMax=iRowOuter-iRowInner+1;
2441   Int_t nClustersFindable=0;
2442   Int_t nClustersSeed=0;
2443
2444   Int_t ncls=(fClusterType==0)?toyTrack->GetNumberOfSpacePoints():toyTrack->GetNumberOfDistSpacePoints();
2445
2446   Int_t rowInner=iRowInner-(iRowInner>62)*63;
2447   Int_t rowOuter=iRowOuter-(iRowOuter>62)*63;
2448
2449   //findable in the current seed sector
2450   Int_t lastROC=-1;
2451   Int_t rocMaxCl=-1;
2452   Int_t nCrossedROCs=0;
2453   Int_t nMaxClROC=0;
2454   Int_t nclROC=0;
2455   Int_t row1=-1;
2456   Int_t row2=-1;
2457   Int_t row1Maxcl=-1;
2458   Int_t row2Maxcl=-1;
2459   for (Int_t icl=0; icl<ncls; ++icl) {
2460     const AliTPCclusterMI *cl=(fClusterType==0)?toyTrack->GetSpacePoint(icl):toyTrack->GetDistortedSpacePoint(icl);
2461     const Int_t roc=cl->GetDetector();
2462     const Int_t row=cl->GetRow();
2463     const Int_t rowGlobal=row+(roc>35)*63;
2464
2465     AliTPCclusterMI *seedCl = seed->GetClusterPointer(rowGlobal);
2466     if (seedCl) {
2467       AliTPCclusterMI *clc=const_cast<AliTPCclusterMI*>(cl);
2468       if (seedCl->GetDetector()==roc&&seedCl->IsUsed()) clc->Use();
2469       clc->SetLabel(seedID,1);
2470
2471     }
2472
2473 //     if (row1<0) row1=rowGlobal;
2474
2475     if ( (roc%36) != (lastROC%36)) {
2476       ++nCrossedROCs;
2477       if (nclROC>nMaxClROC) {
2478         nMaxClROC=nclROC;
2479         rocMaxCl=lastROC;
2480         row1Maxcl=row1;
2481         row2Maxcl=row2;
2482       }
2483       lastROC=roc%36;
2484       nclROC=0;
2485       row1=rowGlobal;
2486     }
2487     ++nclROC;
2488     row2=rowGlobal;
2489
2490     if ( (roc%36)!=(seedSec%36) ) continue;
2491 //     if ( (row<rowInner) || (row>rowOuter) ) continue;
2492     ++nClustersFindable;
2493
2494   }
2495
2496   if (nclROC>nMaxClROC) {
2497     rocMaxCl=lastROC;
2498     nMaxClROC=nclROC;
2499     row1Maxcl=row1;
2500     row2Maxcl=row2;
2501   }
2502
2503   Int_t firstRow=160;
2504   Int_t lastRow=0;
2505   Int_t nClustersInTrack=0;
2506   //found in seed
2507   for (Int_t icl=0; icl<159; ++icl) {
2508     const AliTPCclusterMI *cl=seed->GetClusterPointer(icl);
2509     if (!cl) continue;
2510     ++nClustersInTrack;
2511     const Int_t row=cl->GetRow();
2512     const Int_t rowGlobal=row+(cl->GetDetector()>35)*63;
2513     if (rowGlobal>lastRow)  lastRow=rowGlobal;
2514     if (rowGlobal<firstRow) firstRow=rowGlobal;
2515     if ( (row<rowInner) || (row>rowOuter) ) continue;
2516     ++nClustersSeed;
2517   }
2518
2519   Float_t z0=fEvent->GetZ();
2520   Float_t t0=fEvent->GetT0();
2521
2522   Int_t ctype(fCorrectionType);
2523
2524   Int_t info[5]={0};
2525   CookLabel(seed,.6,info);
2526   Int_t seedLabel=seed->GetLabel();
2527
2528   Int_t labelOrig=toyTrack->GetUniqueID();
2529
2530   AliToyMCTrack *tr2 = const_cast<AliToyMCTrack*>(toyTrack);
2531
2532   (*fStreamer) << "Seeds" <<
2533   //   "iev="             << iev               <<
2534   //   "iseed="           << iseed             <<
2535   //   "itr="             << itr               <<
2536   "z0="             << z0                <<
2537   "t0="             << t0                <<
2538   "vDrift="         << vDrift            <<
2539   "zLength="        << zLength           <<
2540   "fTime0="         << fTime0            <<
2541   "clsType="        << fClusterType      <<
2542   "corrType="       << ctype             <<
2543
2544   "tOrig.="         << &tOrig            <<
2545   "tOrigITS.="      << &tOrigITS         <<
2546
2547   "to.nclFindable=" << nClustersFindable <<
2548   "to.nclTot="      << ncls              <<
2549   "to.label="       << labelOrig         <<
2550   "to.nCrossedROCs="<< nCrossedROCs      <<
2551   "to.rocMax="      << rocMaxCl          <<
2552   "to.rocMaxNcl="   << nMaxClROC         <<
2553   "to.row1Max="     << row1Maxcl         <<
2554   "to.row2Max="     << row2Maxcl         <<
2555
2556   "track.="         << track             <<
2557   "trackITS.="      << &trackITS         <<
2558
2559   "s.RowInner="     << iRowInner         <<
2560   "s.RowOuter="     << iRowOuter         <<
2561   "s.nclMax="       << nClustersSeedMax  <<
2562   "s.ncl="          << nClustersSeed     <<
2563   "s.ID="           << seedID            <<
2564
2565   "tr.firstClRow="  << firstRow          <<
2566   "tr.lastClRow="   << lastRow           <<
2567   "tr.ncl="         << nClustersInTrack  <<
2568   "tr.label="       << seedLabel         <<
2569
2570   "tr.LabelNcl="    << info[0]           <<
2571   "tr.Label2Ncl="   << info[1]           <<
2572   "tr.Label2="      << info[2]           <<
2573   "tr.nclTot="      << info[3]           <<
2574   "tr.Nlabels="     << info[4]           <<
2575
2576   "tr.Sec="         << seedSec           <<
2577
2578   "seed.="           << seed             <<
2579   "toyTrack.="       << tr2              <<
2580   "\n";
2581
2582   if (!isDummy) delete track;
2583 }
2584
2585 //____________________________________________________________________________________
2586 void AliToyMCReconstruction::MarkClustersUsed(AliTPCseed *seed)
2587 {
2588   //
2589   //
2590   //
2591
2592   for (Int_t icl=0; icl<159; ++icl) {
2593     AliTPCclusterMI *cl=seed->GetClusterPointer(icl);
2594     if (cl) cl->Use();
2595   }
2596 }
2597
2598 //____________________________________________________________________________________
2599 void AliToyMCReconstruction::ResetClustersZtoTime(AliTPCseed *seed)
2600 {
2601   //
2602   //
2603   //
2604
2605   for (Int_t icl=0; icl<159; ++icl) {
2606     AliTPCclusterMI *cl=seed->GetClusterPointer(icl);
2607     if (cl) cl->SetZ(cl->GetTimeBin());
2608   }
2609 }
2610
2611 //____________________________________________________________________________________
2612 void AliToyMCReconstruction::DumpTracksToTree(const char* file)
2613 {
2614   //
2615   //
2616   //
2617   ConnectInputFile(file);
2618   if (!fTree) return;
2619
2620   delete fStreamer;
2621   fStreamer=0x0;
2622
2623   TString debugName=fInputFile->GetName();
2624   debugName.ReplaceAll(".root",".AllTracks.root");
2625
2626   gSystem->Exec(Form("test -f %s && rm %s", debugName.Data(), debugName.Data()));
2627   fStreamer=new TTreeSRedirector(debugName.Data());
2628
2629   for (Int_t iev=0;iev<fNmaxEvents;++iev){
2630     fTree->GetEvent(iev);
2631     for (Int_t itr=0; itr<fEvent->GetNumberOfTracks();++itr){
2632       AliToyMCTrack *toyTrack=const_cast<AliToyMCTrack*>(fEvent->GetTrack(itr));
2633       Int_t trackID=toyTrack->GetUniqueID();
2634
2635       (*fStreamer) << "Tracks" <<
2636       "iev="  << iev <<
2637       "itr="  << itr <<
2638       "trackID=" << trackID <<
2639       "track.="  << toyTrack <<
2640       "\n";
2641
2642     }
2643   }
2644
2645   Cleanup();
2646 }
2647
2648 //____________________________________________________________________________________
2649 void AliToyMCReconstruction::SetRieman(const AliTPCseed &seed, AliRieman &rieman)
2650 {
2651   //
2652   //
2653   //
2654
2655   rieman.Reset();
2656   for (Int_t icl=0; icl<159; ++icl) {
2657     const AliTPCclusterMI *cl=seed.GetClusterPointer(icl);
2658     if (!cl) continue;
2659     rieman.AddPoint(cl->GetX(), cl->GetY(), cl->GetZ(),
2660                     TMath::Sqrt(cl->GetSigmaY2()), TMath::Sqrt(cl->GetSigmaZ2()));
2661   }
2662   rieman.Update();
2663 }
2664
2665 //____________________________________________________________________________________
2666 void AliToyMCReconstruction::CopyRieman(const AliRieman &from, AliRieman &to)
2667 {
2668   //
2669   //
2670   //
2671
2672   if (to.GetCapacity()<from.GetCapacity()) return;
2673   to.Reset();
2674
2675   for (Int_t i=0;i<from.GetN();++i) to.AddPoint(from.GetX()[i],from.GetY()[i],from.GetZ()[i],from.GetSy()[i],from.GetSz()[i]);
2676 }
2677