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