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