TPC/Upgrade/AliToyMCReconstruction.cxx

   1
   2 #include <TDatabasePDG.h>
   3 #include <TString.h>
   4 #include <TSystem.h>
   5 #include <TROOT.h>
   6 #include <TFile.h>
   7 #include <TPRegexp.h>
   8
   9 #include <AliExternalTrackParam.h>
  10 #include <AliTPCcalibDB.h>
  11 #include <AliTPCclusterMI.h>
  12 #include <AliTPCSpaceCharge3D.h>
  13 #include <AliTrackerBase.h>
  14 #include <AliTrackPointArray.h>
  15 #include <AliLog.h>
  16 #include <AliTPCParam.h>
  17 #include <AliTPCROC.h>
  18 #include <TTreeStream.h>
  19 #include <AliTPCReconstructor.h>
  20 #include <AliTPCTransform.h>
  21 #include <AliTPCtracker.h>
  22 #include <AliTPCtrackerSector.h>
  23
  24 #include "AliToyMCTrack.h"
  25 #include "AliToyMCEvent.h"
  26
  27 #include "AliToyMCReconstruction.h"
  28
  29 /*
  30
  31
  32
  33 */
  34
  35 //____________________________________________________________________________________
  36 AliToyMCReconstruction::AliToyMCReconstruction() : TObject()
  37 , fSeedingRow(140)
  38 , fSeedingDist(10)
  39 , fClusterType(0)
  40 , fCorrectionType(kNoCorrection)
  41 , fDoTrackFit(kTRUE)
  42 , fUseMaterial(kFALSE)
  43 , fIdealTracking(kFALSE)
  44 , fTime0(-1)
  45 , fCreateT0seed(kFALSE)
  46 , fStreamer(0x0)
  47 , fTree(0x0)
  48 , fEvent(0x0)
  49 , fTPCParam(0x0)
  50 , fSpaceCharge(0x0)
  51 , fkNSectorInner(18) // hard-coded to avoid loading the parameters before
  52 , fInnerSectorArray(0x0)
  53 , fkNSectorOuter(18) // hard-coded to avoid loading the parameters before
  54 , fOuterSectorArray(0x0)
  55 {
  56   //
  57   //  ctor
  58   //
  59   fTPCParam=AliTPCcalibDB::Instance()->GetParameters();
  60
  61 }
  62
  63 //____________________________________________________________________________________
  64 AliToyMCReconstruction::~AliToyMCReconstruction()
  65 {
  66   //
  67   //  dtor
  68   //
  69
  70   delete fStreamer;
  71 //   delete fTree;
  72 }
  73
  74 //____________________________________________________________________________________
  75 void AliToyMCReconstruction::RunReco(const char* file, Int_t nmaxEv)
  76 {
  77   //
  78   // Recostruction from associated clusters
  79   //
  80
  81   TFile f(file);
  82   if (!f.IsOpen() || f.IsZombie()) {
  83     printf("ERROR: couldn't open the file '%s'\n", file);
  84     return;
  85   }
  86
  87  fTree=(TTree*)f.Get("toyMCtree");
  88   if (!fTree) {
  89     printf("ERROR: couldn't read the 'toyMCtree' from file '%s'\n", file);
  90     return;
  91   }
  92
  93   fEvent=0x0;
  94   fTree->SetBranchAddress("event",&fEvent);
  95
  96   // read spacecharge from the Userinfo ot the tree
  97   InitSpaceCharge();
  98
  99   TString debugName=file;
 100   debugName.ReplaceAll(".root","");
 101   debugName.Append(Form(".%1d.%1d_%1d_%1d_%03d_%02d",
 102                         fUseMaterial,fIdealTracking,fClusterType,
 103                         Int_t(fCorrectionType),fSeedingRow,fSeedingDist));
 104   debugName.Append(".debug.root");
 105
 106   gSystem->Exec(Form("test -f %s && rm %s", debugName.Data(), debugName.Data()));
 107   if (!fStreamer) fStreamer=new TTreeSRedirector(debugName.Data());
 108
 109   gROOT->cd();
 110
 111   static AliExternalTrackParam dummySeedT0;
 112   static AliExternalTrackParam dummySeed;
 113   static AliExternalTrackParam dummyTrack;
 114
 115   AliExternalTrackParam t0seed;
 116   AliExternalTrackParam seed;
 117   AliExternalTrackParam track;
 118   AliExternalTrackParam tOrig;
 119
 120   AliExternalTrackParam *dummy;
 121
 122   Int_t maxev=fTree->GetEntries();
 123   if (nmaxEv>0&&nmaxEv<maxev) maxev=nmaxEv;
 124
 125   for (Int_t iev=0; iev<maxev; ++iev){
 126     printf("==============  Processing Event %6d =================\n",iev);
 127     fTree->GetEvent(iev);
 128     for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
 129       printf(" > ======  Processing Track %6d ========  \n",itr);
 130       const AliToyMCTrack *tr=fEvent->GetTrack(itr);
 131       tOrig = *tr;
 132
 133
 134       // set dummy
 135       t0seed    = dummySeedT0;
 136       seed      = dummySeed;
 137       track     = dummyTrack;
 138
 139       Float_t z0=fEvent->GetZ();
 140       Float_t t0=fEvent->GetT0();
 141
 142       Float_t vdrift=GetVDrift();
 143       Float_t zLength=GetZLength(0);
 144
 145       // crate time0 seed, steered by fCreateT0seed
 146 //       printf("t0 seed\n");
 147       fTime0=-1.;
 148       fCreateT0seed=kTRUE;
 149       dummy = GetSeedFromTrack(tr);
 150
 151       if (dummy) {
 152         t0seed = *dummy;
 153         delete dummy;
 154
 155         // crate real seed using the time 0 from the first seed
 156         // set fCreateT0seed now to false to get the seed in z coordinates
 157         fTime0 = t0seed.GetZ()-zLength/vdrift;
 158         fCreateT0seed = kFALSE;
 159 //         printf("seed (%.2g)\n",fTime0);
 160         dummy  = GetSeedFromTrack(tr);
 161         if (dummy) {
 162           seed = *dummy;
 163           delete dummy;
 164
 165           // create fitted track
 166           if (fDoTrackFit){
 167 //             printf("track\n");
 168             dummy = GetFittedTrackFromSeed(tr, &seed);
 169             track = *dummy;
 170             delete dummy;
 171           }
 172
 173           // propagate seed to 0
 174           const Double_t kMaxSnp = 0.85;
 175           const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 176           AliTrackerBase::PropagateTrackTo(&seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
 177
 178         }
 179       }
 180
 181       Int_t ctype(fCorrectionType);
 182
 183       if (fStreamer) {
 184         (*fStreamer) << "Tracks" <<
 185         "iev="         << iev             <<
 186         "z0="          << z0              <<
 187         "t0="          << t0              <<
 188         "fTime0="      << fTime0          <<
 189         "itr="         << itr             <<
 190         "clsType="     << fClusterType    <<
 191         "corrType="    << ctype           <<
 192         "seedRow="     << fSeedingRow     <<
 193         "seedDist="    << fSeedingDist    <<
 194         "vdrift="      << vdrift          <<
 195         "zLength="     << zLength         <<
 196         "t0seed.="     << &t0seed         <<
 197         "seed.="       << &seed           <<
 198         "track.="      << &track          <<
 199         "tOrig.="      << &tOrig          <<
 200         "\n";
 201       }
 202
 203
 204     }
 205   }
 206
 207   delete fStreamer;
 208   fStreamer=0x0;
 209
 210   delete fEvent;
 211   fEvent = 0x0;
 212
 213   delete fTree;
 214   fTree=0x0;
 215   f.Close();
 216 }
 217
 218
 219 //____________________________________________________________________________________
 220 void AliToyMCReconstruction::RunRecoAllClusters(const char* file, Int_t nmaxEv)
 221 {
 222   //
 223   // Reconstruction for seed from associated clusters, but array of clusters:
 224   // Step 1) Filling of cluster arrays
 225   // Step 2) Seeding from clusters associated to tracks
 226   // Step 3) Free track reconstruction using all clusters
 227   //
 228
 229   TFile f(file);
 230   if (!f.IsOpen() || f.IsZombie()) {
 231     printf("ERROR: couldn't open the file '%s'\n", file);
 232     return;
 233   }
 234
 235  fTree=(TTree*)f.Get("toyMCtree");
 236   if (!fTree) {
 237     printf("ERROR: couldn't read the 'toyMCtree' from file '%s'\n", file);
 238     return;
 239   }
 240
 241   fEvent=0x0;
 242   fTree->SetBranchAddress("event",&fEvent);
 243
 244   // read spacecharge from the Userinfo ot the tree
 245   InitSpaceCharge();
 246
 247   TString debugName=file;
 248   debugName.ReplaceAll(".root","");
 249   debugName.Append(Form(".%1d.%1d_%1d_%1d_%03d_%02d",
 250                         fUseMaterial,fIdealTracking,fClusterType,
 251                         Int_t(fCorrectionType),fSeedingRow,fSeedingDist));
 252   debugName.Append(".allClusters.debug.root");
 253
 254   gSystem->Exec(Form("test -f %s && rm %s", debugName.Data(), debugName.Data()));
 255   if (!fStreamer) fStreamer=new TTreeSRedirector(debugName.Data());
 256
 257   gROOT->cd();
 258
 259   static AliExternalTrackParam dummySeedT0;
 260   static AliExternalTrackParam dummySeed;
 261   static AliExternalTrackParam dummyTrack;
 262
 263   AliExternalTrackParam t0seed;
 264   AliExternalTrackParam seed;
 265   AliExternalTrackParam track;
 266   AliExternalTrackParam tOrig;
 267
 268   AliExternalTrackParam *dummy;
 269
 270   Int_t maxev=fTree->GetEntries();
 271   if (nmaxEv>0&&nmaxEv<maxev) maxev=nmaxEv;
 272
 273   // ===========================================================================================
 274   // Loop 1: Fill AliTPCtrackerSector structure
 275   // ===========================================================================================
 276   FillSectorStructure(maxev);
 277
 278   // settings (TODO: find the correct settings)
 279   AliTPCRecoParam *tpcRecoParam = new AliTPCRecoParam();
 280   tpcRecoParam->SetDoKinks(kFALSE);
 281   AliTPCcalibDB::Instance()->GetTransform()->SetCurrentRecoParam(tpcRecoParam);
 282   //tpcRecoParam->Print();
 283
 284   // need AliTPCReconstructor for parameter settings in AliTPCtracker
 285   AliTPCReconstructor *tpcRec   = new AliTPCReconstructor();
 286   tpcRec->SetRecoParam(tpcRecoParam);
 287
 288
 289   // ===========================================================================================
 290   // Loop 2: Seeding from clusters associated to tracks
 291   // TODO: Implement tracking from given seed!
 292   // ===========================================================================================
 293   for (Int_t iev=0; iev<maxev; ++iev){
 294     printf("==============  Processing Event %6d =================\n",iev);
 295     fTree->GetEvent(iev);
 296     for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
 297       printf(" > ======  Processing Track %6d ========  \n",itr);
 298       const AliToyMCTrack *tr=fEvent->GetTrack(itr);
 299       tOrig = *tr;
 300
 301
 302       // set dummy
 303       t0seed    = dummySeedT0;
 304       seed      = dummySeed;
 305       track     = dummyTrack;
 306
 307       Float_t z0=fEvent->GetZ();
 308       Float_t t0=fEvent->GetT0();
 309
 310       Float_t vdrift=GetVDrift();
 311       Float_t zLength=GetZLength(0);
 312
 313       Int_t nClus = 0;
 314
 315       // crate time0 seed, steered by fCreateT0seed
 316       printf("t0 seed\n");
 317       fTime0=-1.;
 318       fCreateT0seed=kTRUE;
 319       dummy = GetSeedFromTrack(tr);
 320
 321       if (dummy) {
 322         t0seed = *dummy;
 323         delete dummy;
 324
 325         // crate real seed using the time 0 from the first seed
 326         // set fCreateT0seed now to false to get the seed in z coordinates
 327         fTime0 = t0seed.GetZ()-zLength/vdrift;
 328         fCreateT0seed = kFALSE;
 329         printf("seed (%.2g)\n",fTime0);
 330         dummy  = GetSeedFromTrack(tr);
 331         if (dummy) {
 332           seed = *dummy;
 333           delete dummy;
 334
 335           // create fitted track
 336           if (fDoTrackFit){
 337             printf("track\n");
 338             dummy = GetFittedTrackFromSeedAllClusters(tr, &seed,nClus);
 339             track = *dummy;
 340             delete dummy;
 341           }
 342
 343           // propagate seed to 0
 344           const Double_t kMaxSnp = 0.85;
 345           const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 346           AliTrackerBase::PropagateTrackTo(&seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
 347
 348         }
 349       }
 350
 351       Int_t ctype(fCorrectionType);
 352
 353       if (fStreamer) {
 354         (*fStreamer) << "Tracks" <<
 355         "iev="         << iev             <<
 356         "z0="          << z0              <<
 357         "t0="          << t0              <<
 358         "fTime0="      << fTime0          <<
 359         "itr="         << itr             <<
 360         "clsType="     << fClusterType    <<
 361         "corrType="    << ctype           <<
 362         "seedRow="     << fSeedingRow     <<
 363         "seedDist="    << fSeedingDist    <<
 364         "vdrift="      << vdrift          <<
 365         "zLength="     << zLength         <<
 366         "nClus="       << nClus           <<
 367         "t0seed.="     << &t0seed         <<
 368         "seed.="       << &seed           <<
 369         "track.="      << &track          <<
 370         "tOrig.="      << &tOrig          <<
 371         "\n";
 372       }
 373
 374
 375     }
 376   }
 377
 378
 379   delete fStreamer;
 380   fStreamer=0x0;
 381
 382   delete fEvent;
 383   fEvent = 0x0;
 384
 385   delete fTree;
 386   fTree=0x0;
 387   f.Close();
 388 }
 389
 390 //____________________________________________________________________________________
 391 void AliToyMCReconstruction::RunRecoAllClustersStandardTracking(const char* file, Int_t nmaxEv)
 392 {
 393   //
 394   // Reconstruction for seed from associated clusters, but array of clusters
 395   // Step 1) Filling of cluster arrays
 396   // Step 2) Use the standard tracking: AliTPCtracker::Clusters2Tracks();
 397   //
 398
 399   TFile f(file);
 400   if (!f.IsOpen() || f.IsZombie()) {
 401     printf("ERROR: couldn't open the file '%s'\n", file);
 402     return;
 403   }
 404
 405  fTree=(TTree*)f.Get("toyMCtree");
 406   if (!fTree) {
 407     printf("ERROR: couldn't read the 'toyMCtree' from file '%s'\n", file);
 408     return;
 409   }
 410
 411   fEvent=0x0;
 412   fTree->SetBranchAddress("event",&fEvent);
 413
 414   // read spacecharge from the Userinfo ot the tree
 415   InitSpaceCharge();
 416
 417   TString debugName=file;
 418   debugName.ReplaceAll(".root","");
 419   debugName.Append(Form(".%1d.%1d_%1d_%1d_%03d_%02d",
 420                         fUseMaterial,fIdealTracking,fClusterType,
 421                         Int_t(fCorrectionType),fSeedingRow,fSeedingDist));
 422   debugName.Append(".allClusters.debug.root");
 423
 424   gSystem->Exec(Form("test -f %s && rm %s", debugName.Data(), debugName.Data()));
 425   if (!fStreamer) fStreamer=new TTreeSRedirector(debugName.Data());
 426
 427   gROOT->cd();
 428
 429   AliExternalTrackParam *dummy;
 430   AliExternalTrackParam *track;
 431
 432   Int_t maxev=fTree->GetEntries();
 433   if (nmaxEv>0&&nmaxEv<maxev) maxev=nmaxEv;
 434
 435
 436   // ===========================================================================================
 437   // Loop 1: Fill AliTPCtrackerSector structure
 438   // ===========================================================================================
 439   FillSectorStructure(maxev);
 440
 441   // ===========================================================================================
 442   // Loop 2: Use the full TPC tracker
 443   // TODO: - check tracking configuration
 444   //       - add clusters and original tracks to output (how?)
 445   // ===========================================================================================
 446
 447   // settings (TODO: find the correct settings)
 448   AliTPCRecoParam *tpcRecoParam = new AliTPCRecoParam();
 449   tpcRecoParam->SetDoKinks(kFALSE);
 450   AliTPCcalibDB::Instance()->GetTransform()->SetCurrentRecoParam(tpcRecoParam);
 451   //tpcRecoParam->Print();
 452
 453   // need AliTPCReconstructor for parameter settings in AliTPCtracker
 454   AliTPCReconstructor *tpcRec   = new AliTPCReconstructor();
 455   tpcRec->SetRecoParam(tpcRecoParam);
 456
 457   // AliTPCtracker
 458   AliTPCtracker *tpcTracker = new AliTPCtracker(fTPCParam);
 459   tpcTracker->SetDebug(10);
 460
 461   // set sector arrays
 462   tpcTracker->SetTPCtrackerSectors(fInnerSectorArray,fOuterSectorArray);
 463   tpcTracker->LoadInnerSectors();
 464   tpcTracker->LoadOuterSectors();
 465
 466   // tracking
 467   tpcTracker->Clusters2Tracks();
 468   //tpcTracker->PropagateForward();
 469   TObjArray *trackArray = tpcTracker->GetSeeds();
 470
 471    for(Int_t iTracks = 0; iTracks < trackArray->GetEntriesFast(); ++iTracks){
 472      printf(" > ======  Fill Track %6d ========  \n",iTracks);
 473
 474      track = (AliExternalTrackParam*)(trackArray->At(iTracks));
 475
 476      if (fStreamer) {
 477        (*fStreamer) << "Tracks" <<
 478          "track.="      << track          <<
 479          "\n";
 480      }
 481    }
 482
 483   delete fStreamer;
 484   fStreamer=0x0;
 485
 486   delete fEvent;
 487   fEvent = 0x0;
 488
 489   delete fTree;
 490   fTree=0x0;
 491   f.Close();
 492 }
 493
 494
 495 //____________________________________________________________________________________
 496 AliExternalTrackParam* AliToyMCReconstruction::GetSeedFromTrack(const AliToyMCTrack * const tr)
 497 {
 498   //
 499   // if we don't have a valid time0 informaion (fTime0) available yet
 500   // assume we create a seed for the time0 estimate
 501   //
 502
 503   // seed point informaion
 504   AliTrackPoint    seedPoint[3];
 505   const AliTPCclusterMI *seedCluster[3]={0x0,0x0,0x0};
 506
 507   // number of clusters to loop over
 508   const Int_t ncls=(fClusterType==0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
 509
 510   UChar_t nextSeedRow=fSeedingRow;
 511   Int_t   nseeds=0;
 512
 513   //assumes sorted clusters
 514   for (Int_t icl=0;icl<ncls;++icl) {
 515     const AliTPCclusterMI *cl=tr->GetSpacePoint(icl);
 516     if (fClusterType==1) cl=tr->GetDistortedSpacePoint(icl);
 517     if (!cl) continue;
 518     // use row in sector
 519     const UChar_t row=cl->GetRow() + 63*(cl->GetDetector()>35);
 520     // skip clusters without proper pad row
 521     if (row>200) continue;
 522
 523     //check seeding row
 524     // if we are in the last row and still miss a seed we use the last row
 525     //   even if the row spacing will not be equal
 526     if (row>=nextSeedRow || icl==ncls-1){
 527       seedCluster[nseeds]=cl;
 528       SetTrackPointFromCluster(cl, seedPoint[nseeds]);
 529       ++nseeds;
 530       nextSeedRow+=fSeedingDist;
 531
 532       if (nseeds==3) break;
 533     }
 534   }
 535
 536   // check we really have 3 seeds
 537   if (nseeds!=3) {
 538     AliError(Form("Seeding failed for parameters %d, %d\n",fSeedingDist,fSeedingRow));
 539     return 0x0;
 540   }
 541
 542   // do cluster correction for fCorrectionType:
 543   //   0 - no correction
 544   //   1 - TPC center
 545   //   2 - average eta
 546   //   3 - ideal
 547   // assign the cluster abs time as z component to all seeds
 548   for (Int_t iseed=0; iseed<3; ++iseed) {
 549     Float_t xyz[3]={0,0,0};
 550     seedPoint[iseed].GetXYZ(xyz);
 551     const Float_t r=TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
 552
 553     const Int_t sector=seedCluster[iseed]->GetDetector();
 554     const Int_t sign=1-2*((sector/18)%2);
 555
 556     if ( (fClusterType == 1) && (fCorrectionType != kNoCorrection) ) {
 557 //       printf("correction type: %d\n",(Int_t)fCorrectionType);
 558
 559       // the settings below are for the T0 seed
 560       // for known T0 the z position is already calculated in SetTrackPointFromCluster
 561       if ( fCreateT0seed ){
 562         if ( fCorrectionType == kTPCCenter  ) xyz[2] = 125.*sign;
 563         //!!! TODO: is this the correct association?
 564         if ( fCorrectionType == kAverageEta ) xyz[2] = TMath::Tan(45./2.*TMath::DegToRad())*r*sign;
 565       }
 566
 567       if ( fCorrectionType == kIdeal      ) xyz[2] = seedCluster[iseed]->GetZ();
 568
 569       //!!! TODO: to be replaced with the proper correction
 570       fSpaceCharge->CorrectPoint(xyz, seedCluster[iseed]->GetDetector());
 571     }
 572
 573     // after the correction set the time bin as z-Position in case of a T0 seed
 574     if ( fCreateT0seed )
 575       xyz[2]=seedCluster[iseed]->GetTimeBin();
 576
 577     seedPoint[iseed].SetXYZ(xyz);
 578   }
 579
 580   const Double_t kMaxSnp = 0.85;
 581   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 582
 583   AliExternalTrackParam *seed = AliTrackerBase::MakeSeed(seedPoint[0], seedPoint[1], seedPoint[2]);
 584   seed->ResetCovariance(10);
 585
 586   if (fCreateT0seed){
 587     // if fTime0 < 0 we assume that we create a seed for the T0 estimate
 588     AliTrackerBase::PropagateTrackTo(seed,0,kMass,5,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
 589     if (TMath::Abs(seed->GetX())>3) {
 590 //       printf("Could not propagate track to 0, %.2f, %.2f, %.2f\n",seed->GetX(),seed->GetAlpha(),seed->Pt());
 591     }
 592   }
 593
 594   return seed;
 595
 596 }
 597
 598 //____________________________________________________________________________________
 599 void AliToyMCReconstruction::SetTrackPointFromCluster(const AliTPCclusterMI *cl, AliTrackPoint &p )
 600 {
 601   //
 602   // make AliTrackPoint out of AliTPCclusterMI
 603   //
 604
 605   if (!cl) return;
 606     Float_t xyz[3]={0.,0.,0.};
 607   //   ClusterToSpacePoint(cl,xyz);
 608   //   cl->GetGlobalCov(cov);
 609   //TODO: what to do with the covariance matrix???
 610   //TODO: the problem is that it is used in GetAngle in AliTrackPoint
 611   //TODO: which is used by AliTrackerBase::MakeSeed to get alpha correct ...
 612   //TODO: for the moment simply assign 1 permill squared
 613   // in AliTrackPoint the cov is xx, xy, xz, yy, yz, zz
 614   //   Float_t cov[6]={xyz[0]*xyz[0]*1e-6,xyz[0]*xyz[1]*1e-6,xyz[0]*xyz[2]*1e-6,
 615   //                   xyz[1]*xyz[1]*1e-6,xyz[1]*xyz[2]*1e-6,xyz[2]*xyz[2]*1e-6};
 616   //   cl->GetGlobalXYZ(xyz);
 617   //   cl->GetGlobalCov(cov);
 618   // voluem ID to add later ....
 619   //   p.SetXYZ(xyz);
 620   //   p.SetCov(cov);
 621   AliTrackPoint *tp=const_cast<AliTPCclusterMI*>(cl)->MakePoint();
 622   p=*tp;
 623   delete tp;
 624   //   cl->Print();
 625   //   p.Print();
 626   p.SetVolumeID(cl->GetDetector());
 627
 628
 629   if ( !fCreateT0seed && !fIdealTracking ) {
 630     p.GetXYZ(xyz);
 631     const Int_t sector=cl->GetDetector();
 632     const Int_t sign=1-2*((sector/18)%2);
 633     const Float_t zT0=( GetZLength(sector) - (cl->GetTimeBin()-fTime0)*GetVDrift() )*sign;
 634 //     printf(" z:  %.2f  %.2f\n",xyz[2],zT0);
 635     xyz[2]=zT0;
 636     p.SetXYZ(xyz);
 637   }
 638
 639
 640   //   p.Rotate(p.GetAngle()).Print();
 641 }
 642
 643 //____________________________________________________________________________________
 644 void AliToyMCReconstruction::ClusterToSpacePoint(const AliTPCclusterMI *cl, Float_t xyz[3])
 645 {
 646   //
 647   // convert the cluster to a space point in global coordinates
 648   //
 649   if (!cl) return;
 650   xyz[0]=cl->GetRow();
 651   xyz[1]=cl->GetPad();
 652   xyz[2]=cl->GetTimeBin(); // this will not be correct at all
 653   Int_t i[3]={0,cl->GetDetector(),cl->GetRow()};
 654   //   printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
 655   fTPCParam->Transform8to4(xyz,i);
 656   //   printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
 657   fTPCParam->Transform4to3(xyz,i);
 658   //   printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
 659   fTPCParam->Transform2to1(xyz,i);
 660   //   printf("%.2f, %.2f, %.2f - %d, %d, %d\n",xyz[0],xyz[1],xyz[2],i[0],i[1],i[2]);
 661 }
 662
 663 //____________________________________________________________________________________
 664 AliExternalTrackParam* AliToyMCReconstruction::GetFittedTrackFromSeed(const AliToyMCTrack *tr, const AliExternalTrackParam *seed)
 665 {
 666   //
 667   //
 668   //
 669
 670   // create track
 671   AliExternalTrackParam *track = new AliExternalTrackParam(*seed);
 672
 673   Int_t ncls=(fClusterType == 0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
 674
 675   const AliTPCROC * roc = AliTPCROC::Instance();
 676
 677   const Double_t kRTPC0  = roc->GetPadRowRadii(0,0);
 678   const Double_t kRTPC1  = roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
 679   const Double_t kMaxSnp = 0.85;
 680   const Double_t kMaxR   = 500.;
 681   const Double_t kMaxZ   = 500.;
 682
 683   //   const Double_t kMaxZ0=220;
 684 //   const Double_t kZcut=3;
 685
 686   const Double_t refX = tr->GetX();
 687
 688   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 689
 690   // loop over all other points and add to the track
 691   for (Int_t ipoint=ncls-1; ipoint>=0; --ipoint){
 692     AliTrackPoint pIn;
 693     const AliTPCclusterMI *cl=tr->GetSpacePoint(ipoint);
 694     if (fClusterType == 1) cl=tr->GetDistortedSpacePoint(ipoint);
 695     SetTrackPointFromCluster(cl, pIn);
 696     if (fCorrectionType != kNoCorrection){
 697       Float_t xyz[3]={0,0,0};
 698       pIn.GetXYZ(xyz);
 699 //       if ( fCorrectionType == kIdeal ) xyz[2] = cl->GetZ();
 700       fSpaceCharge->CorrectPoint(xyz, cl->GetDetector());
 701       pIn.SetXYZ(xyz);
 702     }
 703     // rotate the cluster to the local detector frame
 704     track->Rotate(((cl->GetDetector()%18)*20+10)*TMath::DegToRad());
 705     AliTrackPoint prot = pIn.Rotate(track->GetAlpha());   // rotate to the local frame - non distoted  point
 706     if (TMath::Abs(prot.GetX())<kRTPC0) continue;
 707     if (TMath::Abs(prot.GetX())>kRTPC1) continue;
 708     //
 709     Bool_t res=AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
 710
 711     if (!res) break;
 712
 713     if (TMath::Abs(track->GetZ())>kMaxZ) break;
 714     if (TMath::Abs(track->GetX())>kMaxR) break;
 715 //     if (TMath::Abs(track->GetZ())<kZcut)continue;
 716     //
 717     Double_t pointPos[2]={0,0};
 718     Double_t pointCov[3]={0,0,0};
 719     pointPos[0]=prot.GetY();//local y
 720     pointPos[1]=prot.GetZ();//local z
 721     pointCov[0]=prot.GetCov()[3];//simay^2
 722     pointCov[1]=prot.GetCov()[4];//sigmayz
 723     pointCov[2]=prot.GetCov()[5];//sigmaz^2
 724
 725     if (!track->Update(pointPos,pointCov)) {printf("no update\n"); break;}
 726   }
 727
 728   AliTrackerBase::PropagateTrackTo2(track,refX,kMass,5.,kTRUE,kMaxSnp,0,kFALSE,fUseMaterial);
 729
 730   // rotate fittet track to the frame of the original track and propagate to same reference
 731   track->Rotate(tr->GetAlpha());
 732
 733   AliTrackerBase::PropagateTrackTo2(track,refX,kMass,1.,kFALSE,kMaxSnp,0,kFALSE,fUseMaterial);
 734
 735   return track;
 736 }
 737
 738
 739 //____________________________________________________________________________________
 740 AliExternalTrackParam* AliToyMCReconstruction::GetFittedTrackFromSeedAllClusters(const AliToyMCTrack *tr, const AliExternalTrackParam *seed, Int_t &nClus)
 741 {
 742   //
 743   // Tracking for given seed on an array of clusters
 744   //
 745
 746   // create track
 747   AliExternalTrackParam *track = new AliExternalTrackParam(*seed);
 748
 749   const AliTPCROC * roc = AliTPCROC::Instance();
 750
 751   const Double_t kRTPC0    = roc->GetPadRowRadii(0,0);
 752   const Double_t kRTPC1    = roc->GetPadRowRadii(36,roc->GetNRows(36)-1);
 753   const Int_t kNRowsTPC    = roc->GetNRows(0) + roc->GetNRows(36) - 1;
 754   const Int_t kIRowsTPC    = roc->GetNRows(0) - 1;
 755   const Double_t kMaxSnp   = 0.85;
 756   const Double_t kMaxR     = 500.;
 757   const Double_t kMaxZ     = 500.;
 758   const Double_t roady     = 100.;
 759   const Double_t roadz     = 100.;
 760
 761   const Double_t refX = tr->GetX();
 762
 763   const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 764
 765   Int_t  secCur   = -1;
 766   UInt_t indexCur = 0;
 767   Double_t xCur, yCur, zCur = 0.;
 768
 769   // first propagate seed to outermost row
 770   AliTrackerBase::PropagateTrackTo(track,kRTPC1,kMass,5,kFALSE,kMaxSnp);
 771
 772   // Loop over rows and find the cluster candidates
 773   for( Int_t iRow = kNRowsTPC; iRow >= 0; --iRow ){
 774
 775     // inner or outer sector
 776     Bool_t bInnerSector = kTRUE;
 777     if(iRow > kIRowsTPC) bInnerSector = kFALSE;
 778
 779     // nearest track point/cluster (to be found)
 780     AliTrackPoint nearestPoint;
 781     AliTPCclusterMI *nearestCluster = NULL;
 782
 783     // Inner Sector
 784     if(bInnerSector){
 785
 786       // Propagate to center of pad row and extract parameters
 787       AliTrackerBase::PropagateTrackTo(track,roc->GetPadRowRadii(0,iRow),kMass,5,kFALSE,kMaxSnp);
 788       xCur   = track->GetX();
 789       yCur   = track->GetY();
 790       zCur   = track->GetZ();
 791
 792       secCur = GetSector(track);
 793
 794       // Find the nearest cluster (TODO: correct road settings!)
 795       //Printf("inner tracking here: x = %.2f, y = %.2f, z = %.2f (Row %d Sector %d)",xCur,yCur,zCur,iRow,secCur);
 796       nearestCluster = fInnerSectorArray[secCur%fkNSectorInner][iRow].FindNearest2(yCur,zCur,roady,roadz,indexCur);
 797
 798       // Move to next row if now cluster found
 799       if(!nearestCluster) continue;
 800       //Printf("Nearest Clusters = %d (of %d) ",indexCur,fInnerSectorArray[secCur%fkNSectorInner][iRow].GetN());
 801
 802     }
 803
 804     // Outer sector
 805     else{
 806
 807       // Propagate to center of pad row and extract parameters
 808       AliTrackerBase::PropagateTrackTo(track,roc->GetPadRowRadii(36,iRow-kIRowsTPC-1),kMass,5,kFALSE,kMaxSnp);
 809       xCur   = track->GetX();
 810       yCur   = track->GetY();
 811       zCur   = track->GetZ();
 812
 813       secCur = GetSector(track);
 814
 815       // Find the nearest cluster (TODO: correct road settings!)
 816       //Printf("outer tracking here: x = %.2f, y = %.2f, z = %.2f (Row %d Sector %d)",xCur,yCur,zCur,iRow,secCur);
 817       nearestCluster = fOuterSectorArray[(secCur-fkNSectorInner*2)%fkNSectorOuter][iRow-kIRowsTPC-1].FindNearest2(yCur,zCur,roady,roadz,indexCur);
 818
 819       // Move to next row if now cluster found
 820       if(!nearestCluster) continue;
 821       //Printf("Nearest Clusters = %d (of %d)",indexCur,fOuterSectorArray[(secCur-fkNSectorInner*2)%fkNSectorOuter][iRow-kIRowsTPC-1].GetN());
 822
 823     }
 824
 825     // create track point from cluster
 826     SetTrackPointFromCluster(nearestCluster,nearestPoint);
 827
 828     //Printf("Track point = %.2f %.2f %.2f",nearestPoint.GetX(),nearestPoint.GetY(),nearestPoint.GetZ());
 829
 830     // rotate the cluster to the local detector frame
 831     track->Rotate(((nearestCluster->GetDetector()%18)*20+10)*TMath::DegToRad());
 832     AliTrackPoint prot = nearestPoint.Rotate(track->GetAlpha());   // rotate to the local frame - non distoted  point
 833     if (TMath::Abs(prot.GetX())<kRTPC0) continue;
 834     if (TMath::Abs(prot.GetX())>kRTPC1) continue;
 835
 836     //Printf("Rotated Track point = %.2f %.2f %.2f",prot.GetX(),prot.GetY(),prot.GetZ());
 837
 838     // update track with the nearest track point
 839     Bool_t res=kTRUE;
 840     if (fUseMaterial) res=AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp);
 841     else res=AliTrackerBase::PropagateTrackTo(track,prot.GetX(),kMass,5,kFALSE,kMaxSnp,0,kFALSE,kFALSE);
 842
 843     if (!res) break;
 844
 845     if (TMath::Abs(track->GetZ())>kMaxZ) break;
 846     if (TMath::Abs(track->GetX())>kMaxR) break;
 847     //if (TMath::Abs(track->GetZ())<kZcut)continue;
 848
 849       Double_t pointPos[2]={0,0};
 850       Double_t pointCov[3]={0,0,0};
 851       pointPos[0]=prot.GetY();//local y
 852       pointPos[1]=prot.GetZ();//local z
 853       pointCov[0]=prot.GetCov()[3];//simay^2
 854       pointCov[1]=prot.GetCov()[4];//sigmayz
 855       pointCov[2]=prot.GetCov()[5];//sigmaz^2
 856
 857       if (!track->Update(pointPos,pointCov)) {printf("no update\n");}
 858
 859       ++nClus;
 860   }
 861
 862
 863   // propagation to refX
 864   if (fUseMaterial) AliTrackerBase::PropagateTrackTo2(track,refX,kMass,5.,kTRUE,kMaxSnp);
 865   else AliTrackerBase::PropagateTrackTo2(track,refX,kMass,5.,kTRUE,kMaxSnp,0,kFALSE,kFALSE);
 866
 867   // rotate fittet track to the frame of the original track and propagate to same reference
 868   track->Rotate(tr->GetAlpha());
 869
 870   if (fUseMaterial) AliTrackerBase::PropagateTrackTo2(track,refX,kMass,1.,kFALSE,kMaxSnp);
 871   else AliTrackerBase::PropagateTrackTo2(track,refX,kMass,1.,kFALSE,kMaxSnp,0,kFALSE,kFALSE);
 872
 873   Printf("We have %d clusters in this track!",nClus);
 874
 875   return track;
 876 }
 877
 878
 879 //____________________________________________________________________________________
 880 void AliToyMCReconstruction::InitSpaceCharge()
 881 {
 882   //
 883   // Init the space charge map
 884   //
 885
 886   TString filename="$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps5_50kHz_precal.root";
 887   if (fTree) {
 888     TList *l=fTree->GetUserInfo();
 889     for (Int_t i=0; i<l->GetEntries(); ++i) {
 890       TString s(l->At(i)->GetName());
 891       if (s.Contains("SC_")) {
 892         filename=s;
 893         break;
 894       }
 895     }
 896   }
 897
 898   printf("Initialising the space charge map using the file: '%s'\n",filename.Data());
 899   TFile f(filename.Data());
 900   fSpaceCharge=(AliTPCSpaceCharge3D*)f.Get("map");
 901
 902   //   fSpaceCharge = new AliTPCSpaceCharge3D();
 903   //   fSpaceCharge->SetSCDataFileName("$ALICE_ROOT/TPC/Calib/maps/SC_NeCO2_eps10_50kHz.root");
 904   //   fSpaceCharge->SetOmegaTauT1T2(0.325,1,1); // Ne CO2
 905   // //   fSpaceCharge->SetOmegaTauT1T2(0.41,1,1.05); // Ar CO2
 906   //   fSpaceCharge->InitSpaceCharge3DDistortion();
 907
 908 }
 909
 910 //____________________________________________________________________________________
 911 Double_t AliToyMCReconstruction::GetVDrift() const
 912 {
 913   //
 914   //
 915   //
 916   return fTPCParam->GetDriftV();
 917 }
 918
 919 //____________________________________________________________________________________
 920 Double_t AliToyMCReconstruction::GetZLength(Int_t roc) const
 921 {
 922   //
 923   //
 924   //
 925   if (roc<0 || roc>71) return -1;
 926   return fTPCParam->GetZLength(roc);
 927 }
 928
 929 //____________________________________________________________________________________
 930 TTree* AliToyMCReconstruction::ConnectTrees (const char* files) {
 931   TString s=gSystem->GetFromPipe(Form("ls %s",files));
 932
 933   TTree *tFirst=0x0;
 934   TObjArray *arrFiles=s.Tokenize("\n");
 935
 936   for (Int_t ifile=0; ifile<arrFiles->GetEntriesFast(); ++ifile){
 937     TString name(arrFiles->At(ifile)->GetName());
 938
 939     TPRegexp reg(".*([0-9]_[0-9]_[0-9]_[0-9]{3}_[0-9]{2}).debug.root");
 940     TObjArray *arrMatch=0x0;
 941     arrMatch=reg.MatchS(name);
 942     TString matchName;
 943     if (arrMatch && arrMatch->At(1)) matchName=arrMatch->At(1)->GetName();
 944     else matchName=Form("%02d",ifile);
 945     delete arrMatch;
 946
 947     if (!tFirst) {
 948       TFile *f=TFile::Open(name.Data());
 949       if (!f) continue;
 950       TTree *t=(TTree*)f->Get("Tracks");
 951       if (!t) {
 952         delete f;
 953         continue;
 954       }
 955
 956       t->SetName(matchName.Data());
 957       tFirst=t;
 958     } else {
 959       tFirst->AddFriend(Form("t%s=Tracks",matchName.Data()), name.Data());
 960 //       tFirst->AddFriend(Form("t%d=Tracks",ifile), name.Data());
 961     }
 962   }
 963
 964   tFirst->GetListOfFriends()->Print();
 965   return tFirst;
 966 }
 967
 968 //_____________________________________________________________________________
 969 Int_t AliToyMCReconstruction::LoadOuterSectors() {
 970   //-----------------------------------------------------------------
 971   // This function fills outer TPC sectors with clusters.
 972   // Copy and paste from AliTPCtracker
 973   //-----------------------------------------------------------------
 974   Int_t nrows = fOuterSectorArray->GetNRows();
 975   UInt_t index=0;
 976   for (Int_t sec = 0;sec<fkNSectorOuter;sec++)
 977     for (Int_t row = 0;row<nrows;row++){
 978       AliTPCtrackerRow*  tpcrow = &(fOuterSectorArray[sec%fkNSectorOuter][row]);
 979       Int_t sec2 = sec+2*fkNSectorInner;
 980       //left
 981       Int_t ncl = tpcrow->GetN1();
 982       while (ncl--) {
 983         AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
 984         index=(((sec2<<8)+row)<<16)+ncl;
 985         tpcrow->InsertCluster(c,index);
 986       }
 987       //right
 988       ncl = tpcrow->GetN2();
 989       while (ncl--) {
 990         AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
 991         index=((((sec2+fkNSectorOuter)<<8)+row)<<16)+ncl;
 992         tpcrow->InsertCluster(c,index);
 993       }
 994       //
 995       // write indexes for fast acces
 996       //
 997       for (Int_t i=0;i<510;i++)
 998         tpcrow->SetFastCluster(i,-1);
 999       for (Int_t i=0;i<tpcrow->GetN();i++){
1000         Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1001         tpcrow->SetFastCluster(zi,i);  // write index
1002       }
1003       Int_t last = 0;
1004       for (Int_t i=0;i<510;i++){
1005         if (tpcrow->GetFastCluster(i)<0)
1006           tpcrow->SetFastCluster(i,last);
1007         else
1008           last = tpcrow->GetFastCluster(i);
1009       }
1010     }
1011   return 0;
1012 }
1013
1014
1015 //_____________________________________________________________________________
1016 Int_t  AliToyMCReconstruction::LoadInnerSectors() {
1017   //-----------------------------------------------------------------
1018   // This function fills inner TPC sectors with clusters.
1019   // Copy and paste from AliTPCtracker
1020   //-----------------------------------------------------------------
1021   Int_t nrows = fInnerSectorArray->GetNRows();
1022   UInt_t index=0;
1023   for (Int_t sec = 0;sec<fkNSectorInner;sec++)
1024     for (Int_t row = 0;row<nrows;row++){
1025       AliTPCtrackerRow*  tpcrow = &(fInnerSectorArray[sec%fkNSectorInner][row]);
1026       //
1027       //left
1028       Int_t ncl = tpcrow->GetN1();
1029       while (ncl--) {
1030         AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1031         index=(((sec<<8)+row)<<16)+ncl;
1032         tpcrow->InsertCluster(c,index);
1033       }
1034       //right
1035       ncl = tpcrow->GetN2();
1036       while (ncl--) {
1037         AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1038         index=((((sec+fkNSectorInner)<<8)+row)<<16)+ncl;
1039         tpcrow->InsertCluster(c,index);
1040       }
1041       //
1042       // write indexes for fast acces
1043       //
1044       for (Int_t i=0;i<510;i++)
1045         tpcrow->SetFastCluster(i,-1);
1046       for (Int_t i=0;i<tpcrow->GetN();i++){
1047         Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1048         tpcrow->SetFastCluster(zi,i);  // write index
1049       }
1050       Int_t last = 0;
1051       for (Int_t i=0;i<510;i++){
1052         if (tpcrow->GetFastCluster(i)<0)
1053           tpcrow->SetFastCluster(i,last);
1054         else
1055           last = tpcrow->GetFastCluster(i);
1056       }
1057
1058     }
1059   return 0;
1060 }
1061
1062 //_____________________________________________________________________________
1063 Int_t  AliToyMCReconstruction::GetSector(AliExternalTrackParam *track) {
1064   //-----------------------------------------------------------------
1065   // This function returns the sector number for a given track
1066   //-----------------------------------------------------------------
1067
1068   Int_t sector = -1;
1069
1070   // get the sector number
1071   // rotate point to global system (track is already global!)
1072   Double_t xd[3];
1073   track->GetXYZ(xd);
1074   //track->Local2GlobalPosition(xd,track->GetAlpha());
1075
1076   // use TPCParams to get the sector number
1077   Float_t xyz[3] = {xd[0],xd[1],xd[2]};
1078   Int_t   i[3]   = {0,0,0};
1079   if(fTPCParam){
1080     sector  = fTPCParam->Transform0to1(xyz,i);
1081   }
1082
1083   return sector;
1084 }
1085
1086 //_____________________________________________________________________________
1087 void  AliToyMCReconstruction::FillSectorStructure(Int_t maxev) {
1088   //-----------------------------------------------------------------
1089   // This function fills the sector structure of AliToyMCReconstruction
1090   //-----------------------------------------------------------------
1091
1092   // cluster array of all sectors
1093   fInnerSectorArray = new AliTPCtrackerSector[fkNSectorInner];
1094   fOuterSectorArray = new AliTPCtrackerSector[fkNSectorOuter];
1095
1096   for (Int_t i=0; i<fkNSectorInner; ++i) fInnerSectorArray[i].Setup(fTPCParam,0);
1097   for (Int_t i=0; i<fkNSectorOuter; ++i) fOuterSectorArray[i].Setup(fTPCParam,1);
1098
1099   Int_t count[72][96] = { {0} , {0} };
1100
1101   for (Int_t iev=0; iev<maxev; ++iev){
1102     printf("==============  Fill Clusters: Processing Event %6d  =================\n",iev);
1103     fTree->GetEvent(iev);
1104     for (Int_t itr=0; itr<fEvent->GetNumberOfTracks(); ++itr){
1105       printf(" > ======  Fill Clusters: Processing Track %6d ========  \n",itr);
1106       const AliToyMCTrack *tr=fEvent->GetTrack(itr);
1107
1108       // number of clusters to loop over
1109       const Int_t ncls=(fClusterType==0)?tr->GetNumberOfSpacePoints():tr->GetNumberOfDistSpacePoints();
1110
1111       for(Int_t icl=0; icl<ncls; ++icl){
1112
1113         AliTPCclusterMI *cl=const_cast<AliTPCclusterMI *>(tr->GetSpacePoint(icl));
1114         if (fClusterType==1) cl=const_cast<AliTPCclusterMI *>(tr->GetDistortedSpacePoint(icl));
1115         if (!cl) continue;
1116
1117         Int_t sec = cl->GetDetector();
1118         Int_t row = cl->GetRow();
1119
1120         // set cluster time to cluster Z
1121         cl->SetZ(cl->GetTimeBin());
1122
1123         //Printf("Fill clusters (sector %d row %d): %.2f %.2f %.2f %.2f",sec,row,cl->GetX(),cl->GetY(),cl->GetZ(),cl->GetTimeBin());
1124
1125         // fill arrays for inner and outer sectors (A/C side handled internally)
1126         if (sec<fkNSectorInner*2){
1127           fInnerSectorArray[sec%fkNSectorInner].InsertCluster(const_cast<AliTPCclusterMI*>(cl), count[sec][row], fTPCParam);
1128         }
1129         else{
1130           fOuterSectorArray[(sec-fkNSectorInner*2)%fkNSectorOuter].InsertCluster(const_cast<AliTPCclusterMI*>(cl), count[sec][row], fTPCParam);
1131         }
1132
1133         ++count[sec][row];
1134       }
1135     }
1136   }
1137
1138   // fill the arrays completely
1139   LoadOuterSectors();
1140   LoadInnerSectors();
1141
1142   // // check the arrays
1143   // for (Int_t i=0; i<fkNSectorInner; ++i){
1144   //   for (Int_t j=0; j<fInnerSectorArray[i].GetNRows(); ++j){
1145   //     if(fInnerSectorArray[i][j].GetN()>0){
1146   //    Printf("Inner: Sector %d Row %d : %d",i,j,fInnerSectorArray[i][j].GetN());
1147   //     }
1148   //   }
1149   // }
1150   // for (Int_t i=0; i<fkNSectorInner; ++i){
1151   //   for (Int_t j=0; j<fOuterSectorArray[i].GetNRows(); ++j){
1152   //     if(fOuterSectorArray[i][j].GetN()>0){
1153   //    Printf("Outer: Sector %d Row %d : %d",i,j,fOuterSectorArray[i][j].GetN());
1154   //     }
1155   //   }
1156   // }
1157 }