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