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