TPC/AliTPCCmpNG.C

   1 /// \file AliTPCCmpNG.C
   2 ///
   3 /// version: 1.0
   4 /// description:
   5 ///        define a class TPCGenTrack
   6 ///        save TPC related properties of tracks into a single tree
   7 ///
   8 /// input:
   9 ///        Int_t nEvents      ... nr of events to process
  10 ///        Int_t firstEventNr ... first event number (starts from 0)
  11 ///        char* fnRecTracks .. name of file with reconstructed tracks
  12 ///        char* fnHits ... name of file with hits and Kine Tree
  13 ///        char* fnDigits  ... name of file with digits
  14 ///        char* fnTracks .. output file name, default genTracks.root
  15 ///
  16 /// How to use:
  17 ///  Typical usage:
  18 ///    .L AliTPCCmpNG.C+
  19 ///    TPCFindGenTracks *t = new TPCFindGenTracks("galice.root","tpc.digits.root")
  20 ///    t->Exec();
  21 ///    .q
  22 ///    aliroot
  23 ///    .L  AliTPCCmpNG.C+
  24 ///    TPCCmpTr *t2 = new TPCCmpTr("tpc.tracks.root","genTracks.root","cmpTracks.root");
  25 ///    t2->Exec();
  26 ///
  27 ///  Details:
  28 ///
  29 ///  Step 1 - summurize information from simulation
  30 ///
  31 ///  Compile macro with ACLIC:
  32 ///     .L AliTPCCmpNG.C+
  33 ///  create an object TPCFindGenTracks, which processes information
  34 ///  from simulations. As input it needs:
  35 ///     object gAlice: to get magnetic field
  36 ///     TreeK: to get parameters of generated particles
  37 ///     TreeTR: to get track parameters at the TPC entry
  38 ///     TreeD: to get number of digits and digits pattern
  39 ///                for a given track in TPC
  40 ///  These input objects can be in different files, gAlice, TreeK and
  41 ///  TreeTR are in the file fnHits, TreeD in the file fnDigits (can be
  42 ///  the same as fnHits. Output is written to the file fnRes
  43 ///  ("genTracks.root" by default). Use can specify number of
  44 ///  events to process and the first event number:
  45 ///    TPCFindGenTracks *t = new TPCFindGenTracks("galice.root","tpc.digits.root","genTracks.root",1,0)
  46 ///  The filenames in the example on previous line are defaults, user can
  47 ///  specify just the file name with gAlice object (and TreeTR and TreeK),
  48 ///  so equivalent shorter initialization is:
  49 ///    TPCFindGenTracks *t = new TPCFindGenTracks("galice.root")
  50 ///  The task is done by calling Exec() method:
  51 ///    t->Exec();
  52 ///  User can set different debug levels by invoking:
  53 ///    t->SetDebug(debugLevel)
  54 ///  Number of events to process and the first event number can be
  55 ///  specified as parameters to Exec:
  56 ///    t->Exec(nEvents, firstEvent)
  57 ///  Then you have to quit root to get rid of problems with deleting gAlice
  58 ///  object (it is not deleted, but read again in the following step):
  59 ///
  60 ///  Step 2 - compare reconstructed tracks with simulated
  61 ///
  62 ///  Load (and compile) the macro:
  63 ///   .L AliTPCCmpNG.C+
  64 ///  Create object TPCCmpTr, which does the comparison. As input it requires
  65 ///  name of the file with reconstructed TPC tracks. You can specify
  66 ///  name of the file with genTrack tree (summarized info about simulation),
  67 ///  file with gAlice object, output file name, number of events to process
  68 ///  and first event number:
  69 ///  TPCCmpTr *t2 = new TPCCmpTr("tpc.tracks.root","genTracks.root","cmpTracks.root","galice.root",1,0);
  70 ///  The interface is quite similar to the TPCFindGenTracks class.
  71 ///  Then just invoke Exec() method:
  72 ///  t2->Exec();
  73 ///
  74 ///  Step 3 - study the results
  75 ///
  76 ///  Load the outoput TTree and you can do Draw(), Scan() or other
  77 ///  usual things to do with TTree:
  78 ///  TFile *f = new TFile("cmpTracks.root")
  79 ///  TTree *t = (TTree*)f->Get("TPCcmpTracks")
  80 ///  t->Draw("fVDist[3]","fReconstructed")
  81 ///
  82 /// History:
  83 ///
  84 /// 24.09.02 - first version
  85 /// 24.01.03 - v7, before change from TPC Special Hits to TrackReferences
  86 /// 26.01.03 - change from TPC Special Hits to TrackReferences
  87 ///            (loop over TreeTR instead of TreeH)
  88 /// 28.01.03 - v8 last version before removing TPC special point
  89 /// 28.01.03 - remove TPC special point, loop over TreeH
  90 ///            store TParticle and AliTrack
  91 /// 29.01.03 - v9 last version before moving the loop over rec. tracks
  92 ///            into separate step
  93 /// 03.02.03 - rename to AliTPCCmpNG.C, remove the part with rec. tracks
  94 ///            (will be addded in a macro AliTPCCmpTr.C
  95 ///
  96 /// \author Jiri Chudoba
  97 /// \date 24.09.2002
  98
  99 #if !defined(__CINT__) || defined(__MAKECINT__)
 100 #include "iostream.h"
 101 #include <stdio.h>
 102 #include <string.h>
 103 #include "Rtypes.h"
 104 #include "TFile.h"
 105 #include "TTree.h"
 106 #include "TString.h"
 107 #include "TBenchmark.h"
 108 #include "TStopwatch.h"
 109 #include "TParticle.h"
 110 #include "AliRun.h"
 111 #include "AliStack.h"
 112 #include "AliTPCtrack.h"
 113 #include "AliSimDigits.h"
 114 #include "AliTPCParam.h"
 115 #include "TParticle.h"
 116 #include "AliTPC.h"
 117 #include "AliDetector.h"
 118 #include "AliTrackReference.h"
 119 #include "TSystem.h"
 120 #include "TTimer.h"
 121 // #include "AliConst.h"
 122 #endif
 123
 124 #include "AliTPCTracking.C"
 125
 126 // include ML.C:
 127 ////////////////////////////////////////////////////////////////////////
 128 //
 129 // name: ML.C  (Macro Library - collection of functions used in diff macros
 130 // date: 08.05.2002
 131 // last update: 08.05.2002
 132 // author: Jiri Chudoba
 133 // version: 1.0
 134 // description:
 135 //
 136 // History:
 137 //
 138 // 08.05.02 - first version
 139 //
 140 ////////////////////////////////////////////////////////////////////////
 141
 142 #if !defined(__CINT__) || defined(__MAKECINT__)
 143 #include "iostream.h"
 144 #include "Rtypes.h"
 145 #include "TSystem.h"
 146 #include "TTimer.h"
 147 #include "Getline.h"
 148 #include "TChain.h"
 149 #include "TString.h"
 150 #include "TFile.h"
 151 #include "AliRun.h"
 152
 153 void WaitForReturn();
 154 Bool_t ImportgAlice(TFile *file);
 155 TFile* OpenAliceFile(char *fn, Bool_t importgAlice = kFALSE, char *mode = "read");
 156 Int_t Chain(TString baseDir, TString subDirNameMask, TString fn, TChain& chain);
 157
 158 #endif
 159
 160 ////////////////////////////////////////////////////////////////////////
 161 void WaitForReturn()
 162 {
 163 //
 164 // wait until user press return;
 165 //
 166   char    *input;
 167   Bool_t done = kFALSE;
 168   TTimer  *timer = new TTimer("gSystem->ProcessEvents();", 50, kFALSE);
 169
 170   do {
 171     timer->TurnOn();
 172     timer->Reset();
 173     // Now let's read the input, we can use here any
 174     // stdio or iostream reading methods. like std::cin >> myinputl;
 175     input = Getline("Type <return> to continue: ");
 176     timer->TurnOff();
 177     if (input) done = kTRUE;
 178   } while (!done);
 179 }
 180
 181 ////////////////////////////////////////////////////////////////////////
 182
 183 Int_t Chain(TString baseDir, TString subDirNameMask, TString fn, TChain& chain) {
 184 // chain all files fn in the subdirectories which match subDirName
 185 // return number of chained files
 186
 187 // open baseDir, loop over subdirs
 188
 189   void *dirp = gSystem->OpenDirectory(baseDir);
 190   if (!dirp) {
 191     cerr<<"Could not open base directory "<<baseDir.Data()<<endl;
 192     return 0;
 193   }
 194   const char *afile;
 195   Long_t id, size, flag, modTime;
 196   Int_t rc = 0;
 197   char relName[100];
 198   while((afile = gSystem->GetDirEntry(dirp))) {
 199 //    printf("file: %s\n",afile);
 200     if (strstr(afile,subDirNameMask.Data())) {
 201       sprintf(relName,"%s/%s/%s",baseDir.Data(),afile,fn.Data());
 202 //      cerr<<"relName = "<<relName<<endl;
 203       if(!gSystem->GetPathInfo(relName, &id, &size, &flag, &modTime)) {
 204         rc += chain.Add(relName);
 205       }
 206     }
 207   }
 208   gSystem->FreeDirectory(dirp);
 209 //  cerr<<"rc = "<<rc<<endl;
 210   return rc;
 211 }
 212 ////////////////////////////////////////////////////////////////////////
 213 Bool_t ImportgAlice(TFile *file) {
 214 // read in gAlice object from the file
 215   gAlice = (AliRun*)file->Get("gAlice");
 216   if (!gAlice)  return kFALSE;
 217   return kTRUE;
 218 }
 219 ////////////////////////////////////////////////////////////////////////
 220 TFile* OpenAliceFile(char *fn, Bool_t importgAlice, char *mode) {
 221   TFile *file = TFile::Open(fn,mode);
 222   if (!file->IsOpen()) {
 223     cerr<<"OpenAliceFile: cannot open file "<<fn<<" in mode "
 224         <<mode<<endl;
 225     return 0;
 226   }
 227   if (!importgAlice) return file;
 228   if (ImportgAlice(file)) return file;
 229   return 0;
 230 }
 231 ////////////////////////////////////////////////////////////////////////
 232
 233
 234 ////////////////////////////////////////////////////////////////////////
 235 //
 236 // Start of implementation of the class TPCGenInfo
 237 //
 238 ////////////////////////////////////////////////////////////////////////
 239 class TPCGenInfo: public TObject {
 240
 241 public:
 242   TPCGenInfo();
 243   AliTrackReference *fTrackRef;   ///< track reference saved in the output tree
 244   TParticle *fParticle;           ///< generated particle
 245   Int_t fLabel;                   ///< track label
 246
 247   Int_t fRowsWithDigitsInn;    ///< number of rows with digits in the inner sectors
 248   Int_t fRowsWithDigits;       ///< number of rows with digits in the outer sectors
 249   Int_t fRowsTrackLength;      ///< last - first row with digit
 250   Int_t fDigitsInSeed;         ///< digits in the default seed rows
 251
 252   ClassDef(TPCGenInfo,1)  // container for
 253 };
 254 ClassImp(TPCGenInfo)
 255 ////////////////////////////////////////////////////////////////////////
 256 TPCGenInfo::TPCGenInfo()
 257 {
 258   fTrackRef = 0;
 259   fParticle = 0;
 260 }
 261 ////////////////////////////////////////////////////////////////////////
 262
 263 ////////////////////////////////////////////////////////////////////////
 264 //
 265 // End of implementation of the class TPCGenInfo
 266 //
 267 ////////////////////////////////////////////////////////////////////////
 268
 269
 270
 271 ////////////////////////////////////////////////////////////////////////
 272 //
 273 // Start of implementation of the class digitRow
 274 //
 275 ////////////////////////////////////////////////////////////////////////
 276 const Int_t kgRowBytes = 32;
 277
 278 class digitRow: public TObject {
 279
 280 public:
 281   digitRow();
 282 //  digitRow(){;}
 283   virtual ~digitRow(){;}
 284   void SetRow(Int_t row);
 285   Bool_t TestRow(Int_t row);
 286   digitRow & operator=(const digitRow &digOld);
 287   Int_t RowsOn(Int_t upto=8*kgRowBytes);
 288   Int_t Last();
 289   Int_t First();
 290   void Reset();
 291
 292 //private:
 293   UChar_t fDig[kgRowBytes];
 294
 295   ClassDef(digitRow,1)  // container for digit pattern
 296 };
 297 ClassImp(digitRow)
 298 ////////////////////////////////////////////////////////////////////////
 299 digitRow::digitRow()
 300 {
 301   Reset();
 302 }
 303 ////////////////////////////////////////////////////////////////////////
 304 digitRow & digitRow::operator=(const digitRow &digOld)
 305 {
 306   for (Int_t i = 0; i<kgRowBytes; i++) fDig[i] = digOld.fDig[i];
 307   return (*this);
 308 }
 309 ////////////////////////////////////////////////////////////////////////
 310 void digitRow::SetRow(Int_t row)
 311 {
 312   if (row >= 8*kgRowBytes) {
 313     cerr<<"digitRow::SetRow: index "<<row<<" out of bounds."<<endl;
 314     return;
 315   }
 316   Int_t iC = row/8;
 317   Int_t iB = row%8;
 318   SETBIT(fDig[iC],iB);
 319 }
 320
 321 ////////////////////////////////////////////////////////////////////////
 322 Bool_t digitRow::TestRow(Int_t row)
 323 {
 324 /// return kTRUE if row is on
 325
 326   Int_t iC = row/8;
 327   Int_t iB = row%8;
 328   return TESTBIT(fDig[iC],iB);
 329 }
 330 ////////////////////////////////////////////////////////////////////////
 331 Int_t digitRow::RowsOn(Int_t upto)
 332 {
 333 /// returns number of rows with a digit
 334 /// count only rows less equal row number upto
 335
 336   Int_t total = 0;
 337   for (Int_t i = 0; i<kgRowBytes; i++) {
 338     for (Int_t j = 0; j < 8; j++) {
 339       if (i*8+j > upto) return total;
 340       if (TESTBIT(fDig[i],j))  total++;
 341     }
 342   }
 343   return total;
 344 }
 345 ////////////////////////////////////////////////////////////////////////
 346 void digitRow::Reset()
 347 {
 348 /// resets all rows to zero
 349
 350   for (Int_t i = 0; i<kgRowBytes; i++) {
 351     fDig[i] <<= 8;
 352   }
 353 }
 354 ////////////////////////////////////////////////////////////////////////
 355 Int_t digitRow::Last()
 356 {
 357 /// returns the last row number with a digit
 358 /// returns -1 if now digits
 359
 360   for (Int_t i = kgRowBytes-1; i>=0; i--) {
 361     for (Int_t j = 7; j >= 0; j--) {
 362       if TESTBIT(fDig[i],j) return i*8+j;
 363     }
 364   }
 365   return -1;
 366 }
 367 ////////////////////////////////////////////////////////////////////////
 368 Int_t digitRow::First()
 369 {
 370 /// returns the first row number with a digit
 371 /// returns -1 if now digits
 372
 373   for (Int_t i = 0; i<kgRowBytes; i++) {
 374     for (Int_t j = 0; j < 8; j++) {
 375       if (TESTBIT(fDig[i],j)) return i*8+j;
 376     }
 377   }
 378   return -1;
 379 }
 380
 381 ////////////////////////////////////////////////////////////////////////
 382 //
 383 // end of implementation of a class digitRow
 384 //
 385 ////////////////////////////////////////////////////////////////////////
 386 ////////////////////////////////////////////////////////////////////////
 387 //
 388 // Start of implementation of the class TPCFindGenTracks
 389 //
 390 ////////////////////////////////////////////////////////////////////////
 391
 392 class TPCFindGenTracks {
 393
 394 public:
 395   TPCFindGenTracks();
 396   TPCFindGenTracks(char* fnHits,
 397                    char* fnDigits ="tpc.digits.root",
 398                    char* fnRes    ="genTracks.root",
 399                    Int_t nEvents=1, Int_t firstEvent=0);
 400   virtual ~TPCFindGenTracks();
 401   void Reset();
 402   Int_t Exec();
 403   Int_t Exec(Int_t nEvents, Int_t firstEventNr);
 404   void CreateTreeGenTracks();
 405   void CloseOutputFile();
 406   Int_t TreeKLoop();
 407   Int_t TreeTRLoop();
 408   Int_t TreeDLoop();
 409   void SetFirstEventNr(Int_t i) {fFirstEventNr = i;}
 410   void SetNEvents(Int_t i) {fNEvents = i;}
 411   void SetDebug(Int_t level) {fDebug = level;}
 412
 413   Float_t TR2LocalX(AliTrackReference *trackRef,
 414                     AliTPCParam *paramTPC);
 415
 416 public:
 417   Int_t fDebug;                   //!< debug flag
 418   Int_t fEventNr;                 //!< current event number
 419   Int_t fLabel;                   //!< track label
 420   Int_t fNEvents;                 //!< number of events to process
 421   Int_t fFirstEventNr;            //!< first event to process
 422   Int_t fNParticles;              //!< number of particles in TreeK
 423   TTree *fTreeGenTracks;          //!< output tree with generated tracks
 424   char *fFnRes;                   //!< output file name with stored tracks
 425   char *fFnHits;                  //!< input file name with hits
 426   char *fFnDigits;                //!< input file name with digits
 427   TFile *fFileGenTracks;             //!< output file with stored fTreeGenTracks
 428   TFile *fFileHits;               //!< input file with hits
 429   TFile *fFileTreeD;              //!< input file with digits
 430   digitRow *fDigitRow;            //!< pointer to the object saved in Branch
 431   digitRow *fContainerDigitRow;   //!< big container for partial information
 432   AliTrackReference *fTrackRef;   //!< track reference saved in the output tree
 433   AliTrackReference *fContainerTR;//!< big container for partial information
 434   Int_t *fIndexTR;                //!< index of particle label in the fContainerTR
 435   Int_t fLastIndexTR;             //!< last used index in fIndexTR
 436
 437   AliTPCParam* fParamTPC;         //!< AliTPCParam
 438
 439   Double_t fVPrim[3];             //!< primary vertex position
 440   Double_t fVDist[4];             //!< distance of the particle vertex from primary vertex
 441                                   // the fVDist[3] contains size of the 3-vector
 442   TParticle *fParticle;           //!< generated particle
 443
 444   Int_t fRowsWithDigitsInn;       //!< number of rows with digits in the inner sectors
 445   Int_t fRowsWithDigits;          //!< number of rows with digits in the outer sectors
 446   Int_t fRowsTrackLength;         //!< last - first row with digit
 447   Int_t fDigitsInSeed;            //!< digits in the default seed rows
 448
 449 private:
 450
 451 // some constants for the original non-pareller tracking (by Y.Belikov)
 452   static const Int_t seedRow11 = 158;  ///< nRowUp - 1
 453   static const Int_t seedRow12 = 139;  ///< nRowUp - 1 - (Int_t) 0.125*nRowUp
 454   static const Int_t seedRow21 = 149;  ///< seedRow11 - shift
 455   static const Int_t seedRow22 = 130;  ///< seedRow12 - shift
 456   static const Double_t kRaddeg = 180./TMath::Pi();
 457
 458   static const Int_t fgMaxIndexTR = 50000; ///< maximum number of tracks with a track ref
 459   static const Int_t fgMaxParticles = 2000000; ///< maximum number of generated particles
 460   static const Double_t fgPtCut = .001; ///< do not store particles with generated pT less than this
 461   static const Float_t fgTrackRefLocalXMax = 82.95;
 462   static const Float_t fgTrackRefLocalXMaxDelta = 500.;
 463
 464   ClassDef(TPCFindGenTracks,1)    // class which creates and fills tree with TPCGenTrack objects
 465 };
 466 ClassImp(TPCFindGenTracks)
 467
 468 ////////////////////////////////////////////////////////////////////////
 469 TPCFindGenTracks::TPCFindGenTracks()
 470 {
 471   Reset();
 472 }
 473
 474 ////////////////////////////////////////////////////////////////////////
 475 TPCFindGenTracks::TPCFindGenTracks(char* fnHits, char* fnDigits,
 476                                    char* fnRes,
 477                                    Int_t nEvents, Int_t firstEvent)
 478 {
 479   Reset();
 480   fFirstEventNr = firstEvent;
 481   fEventNr = firstEvent;
 482   fNEvents = nEvents;
 483   fFnRes = fnRes;
 484   fFnHits = fnHits;
 485   fFnDigits = fnDigits;
 486 }
 487 ////////////////////////////////////////////////////////////////////////
 488 void TPCFindGenTracks::Reset()
 489 {
 490   fDigitRow = 0;
 491   fEventNr = 0;
 492   fNEvents = 0;
 493   fTreeGenTracks = 0;
 494   fFnRes = "genTracks.root";
 495   fFnHits = "rfio:galice.root";
 496   fFnDigits = "rfio:its.tpc.trd.digits.root";
 497   fFileGenTracks = 0;
 498   fFileHits =0;
 499   fFileTreeD =0;
 500   fContainerDigitRow = 0;
 501   fContainerTR = 0;
 502   fIndexTR = 0;
 503   fLastIndexTR = -1;
 504   fParticle = 0;
 505   fTrackRef = 0;
 506   fDebug = 0;
 507   fVPrim[0] = -1000.;
 508   fVPrim[1] = -1000.;
 509   fVPrim[2] = -1000.;
 510   fParamTPC = 0;
 511
 512 }
 513 ////////////////////////////////////////////////////////////////////////
 514 TPCFindGenTracks::~TPCFindGenTracks()
 515 {
 516   ;
 517 }
 518 ////////////////////////////////////////////////////////////////////////
 519 Int_t TPCFindGenTracks::Exec(Int_t nEvents, Int_t firstEventNr)
 520 {
 521   fNEvents = nEvents;
 522   fFirstEventNr = firstEventNr;
 523   return Exec();
 524 }
 525
 526 ////////////////////////////////////////////////////////////////////////
 527 Int_t TPCFindGenTracks::Exec()
 528 {
 529   TStopwatch timer;
 530   timer.Start();
 531
 532   fDigitRow = new digitRow();
 533   CreateTreeGenTracks();
 534   if (!fTreeGenTracks) return 1;
 535   fFileHits = OpenAliceFile(fFnHits,kTRUE,"read"); //gAlice is read here
 536   if (!fFileHits) return 1;
 537   fFileHits->cd();
 538   SetFieldFactor();
 539
 540   fFileTreeD = TFile::Open(fFnDigits,"read");
 541   if (!fFileTreeD->IsOpen()) {
 542     cerr<<"Cannot open file "<<fFnDigits<<endl;
 543     return 1;
 544   }
 545
 546   fParamTPC = LoadTPCParam(fFileTreeD);
 547   if (!fParamTPC) {
 548     cerr<<"TPC parameters not found and could not be created"<<endl;
 549     return 1;
 550   }
 551
 552   for (fEventNr = fFirstEventNr; fEventNr < fFirstEventNr+fNEvents;
 553        fEventNr++) {
 554     fNParticles = gAlice->GetEvent(fEventNr);
 555     fContainerDigitRow = new digitRow[fNParticles];
 556     fContainerTR = new AliTrackReference[fgMaxIndexTR];
 557     fIndexTR = new Int_t[fNParticles];
 558     for (Int_t i = 0; i<fNParticles; i++) {
 559       fIndexTR[i] = -1;
 560     }
 561
 562     cout<<"Start to process event "<<fEventNr<<endl;
 563     cout<<"\tfNParticles = "<<fNParticles<<endl;
 564     if (fDebug>2) cout<<"\tStart loop over TreeD"<<endl;
 565     if (TreeDLoop()>0) return 1;
 566     if (fDebug>2) cout<<"\tStart loop over TreeTR"<<endl;
 567     if (TreeTRLoop()>0) return 1;
 568     if (fDebug>2) cout<<"\tStart loop over TreeK"<<endl;
 569     if (TreeKLoop()>0) return 1;
 570     if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
 571
 572     delete [] fContainerDigitRow;
 573     delete [] fContainerTR;
 574     delete [] fIndexTR;
 575   }
 576
 577   CloseOutputFile();
 578
 579   cerr<<"Exec finished"<<endl;
 580   fFileHits->cd();
 581   delete gAlice;
 582   fFileHits->Close();
 583   delete fFileHits;
 584
 585   timer.Stop();
 586   timer.Print();
 587   return 0;
 588 }
 589 ////////////////////////////////////////////////////////////////////////
 590 void TPCFindGenTracks::CreateTreeGenTracks()
 591 {
 592   fFileGenTracks = TFile::Open(fFnRes,"RECREATE");
 593   if (!fFileGenTracks) {
 594     cerr<<"Error in CreateTreeGenTracks: cannot open file "<<fFnRes<<endl;
 595     return;
 596   }
 597   fTreeGenTracks = new TTree("genTracksTree","genTracksTree");
 598   TBranch *branchBits = fTreeGenTracks->Branch("bitsBranch", "digitRow", &fDigitRow, 4000, 0);
 599   if (!branchBits) {
 600     cerr<<"Error in CreateTreeGenTracks: cannot create branch."<<endl;
 601     return;
 602   }
 603   fTreeGenTracks->Branch("fEventNr",&fEventNr,"fEventNr/I");
 604   fTreeGenTracks->Branch("fLabel",&fLabel,"fLabel/I");
 605   fTreeGenTracks->Branch("fRowsWithDigitsInn",&fRowsWithDigitsInn,"fRowsWithDigitsInn/I");
 606   fTreeGenTracks->Branch("fRowsWithDigits",&fRowsWithDigits,"fRowsWithDigits/I");
 607   fTreeGenTracks->Branch("fRowsTrackLength",&fRowsTrackLength,"fRowsTrackLength/I");
 608   fTreeGenTracks->Branch("fDigitsInSeed",&fDigitsInSeed,"fDigitsInSeed/I");
 609
 610   fTreeGenTracks->Branch("Particle","TParticle",&fParticle);
 611   fTreeGenTracks->Branch("fVDist",&fVDist,"fVDist[4]/D");
 612   fTreeGenTracks->Branch("TR","AliTrackReference",&fTrackRef);
 613
 614   fTreeGenTracks->AutoSave();
 615 }
 616 ////////////////////////////////////////////////////////////////////////
 617 void TPCFindGenTracks::CloseOutputFile()
 618 {
 619   if (!fFileGenTracks) {
 620     cerr<<"File "<<fFnRes<<" not found as an open file."<<endl;
 621     return;
 622   }
 623   fFileGenTracks->cd();
 624   fTreeGenTracks->Write();
 625   delete fTreeGenTracks;
 626   fFileGenTracks->Close();
 627   delete fFileGenTracks;
 628   return;
 629 }
 630
 631 ////////////////////////////////////////////////////////////////////////
 632 Int_t TPCFindGenTracks::TreeKLoop()
 633 {
 634 /// open the file with treeK
 635 /// loop over all entries there and save information about some tracks
 636
 637   fFileHits->cd();
 638   if (fDebug > 0) {
 639     cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
 640         <<fEventNr<<endl;
 641   }
 642   AliStack * stack = gAlice->Stack();
 643   if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
 644
 645 // not all generators give primary vertex position. Take the vertex
 646 // of the particle 0 as primary vertex.
 647   fParticle = stack->ParticleFromTreeK(0);
 648   fVPrim[0] = fParticle->Vx();
 649   fVPrim[1] = fParticle->Vy();
 650   fVPrim[2] = fParticle->Vz();
 651
 652   for (Int_t iParticle = 0; iParticle < fNParticles; iParticle++) {
 653 //  for (Int_t iParticle = 0; iParticle < fDebug; iParticle++) {
 654
 655 // load only particles with TR
 656     if (fIndexTR[iParticle] < 0) continue;
 657     fParticle = stack->ParticleFromTreeK(iParticle);
 658     if (fDebug > 3 && iParticle < 10) {
 659       cout<<"processing particle "<<iParticle<<" ";
 660       fParticle->Print();
 661     }
 662
 663 // fill the tree
 664
 665     fLabel = iParticle;
 666     fVDist[0] = fParticle->Vx()-fVPrim[0];
 667     fVDist[1] = fParticle->Vy()-fVPrim[1];
 668     fVDist[2] = fParticle->Vz()-fVPrim[2];
 669     fVDist[3] = TMath::Sqrt(fVDist[0]*fVDist[0]+fVDist[1]*fVDist[1]+fVDist[2]*fVDist[2]);
 670     fDigitRow = &(fContainerDigitRow[iParticle]);
 671     fRowsWithDigitsInn = fDigitRow->RowsOn(63); // 63 = number of inner rows
 672     fRowsWithDigits = fDigitRow->RowsOn();
 673     fRowsTrackLength = fDigitRow->Last() - fDigitRow->First();
 674     if (fDebug > 2 && iParticle < 10) {
 675       cerr<<"Fill track with a label "<<iParticle<<endl;
 676     }
 677     fDigitsInSeed = 0;
 678     if (fDigitRow->TestRow(seedRow11) && fDigitRow->TestRow(seedRow12))
 679       fDigitsInSeed = 1;
 680     if (fDigitRow->TestRow(seedRow21) && fDigitRow->TestRow(seedRow22))
 681       fDigitsInSeed += 10;
 682
 683     if (fIndexTR[iParticle] >= 0) {
 684       fTrackRef = &(fContainerTR[fIndexTR[iParticle]]);
 685 //      cerr<<"Debug: fTrackRef->X() = "<<fTrackRef->X()<<endl;
 686     } else {
 687       fTrackRef->SetTrack(-1);
 688     }
 689
 690     fTreeGenTracks->Fill();
 691
 692   }
 693   fTreeGenTracks->AutoSave();
 694 //  delete gAlice; gAlice = 0;
 695 //  fFileHits->Close();
 696
 697   if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
 698
 699   return 0;
 700 }
 701 ////////////////////////////////////////////////////////////////////////
 702 Int_t TPCFindGenTracks::TreeDLoop()
 703 {
 704 /// open the file with treeD
 705 /// loop over all entries there and save information about some tracks
 706
 707
 708 //  Int_t nrow_up=fParamTPC->GetNRowUp();
 709 //  Int_t nrows=fParamTPC->GetNRowLow()+nrow_up;
 710   Int_t nInnerSector = fParamTPC->GetNInnerSector();
 711   Int_t rowShift = 0;
 712   Int_t zero=fParamTPC->GetZeroSup();
 713 //  Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
 714
 715   char treeDName[100];
 716   sprintf(treeDName,"TreeD_75x40_100x60_150x60_%d",fEventNr);
 717   TTree *treeD=(TTree*)fFileTreeD->Get(treeDName);
 718   AliSimDigits digitsAddress, *digits=&digitsAddress;
 719   treeD->GetBranch("Segment")->SetAddress(&digits);
 720
 721   Int_t sectorsByRows=(Int_t)treeD->GetEntries();
 722   if (fDebug > 1) cout<<"\tsectorsByRows = "<<sectorsByRows<<endl;
 723   for (Int_t i=0; i<sectorsByRows; i++) {
 724 //  for (Int_t i=5720; i<sectorsByRows; i++) {
 725     if (!treeD->GetEvent(i)) continue;
 726     Int_t sec,row;
 727     fParamTPC->AdjustSectorRow(digits->GetID(),sec,row);
 728     if (fDebug > 1) cout<<sec<<' '<<row<<"                          \r";
 729 //    cerr<<sec<<' '<<row<<endl;
 730
 731 // here I expect that upper sectors follow lower sectors
 732     if (sec > nInnerSector) rowShift = fParamTPC->GetNRowLow();
 733     digits->First();
 734     do {
 735       Int_t iRow=digits->CurrentRow();
 736       Int_t iColumn=digits->CurrentColumn();
 737       Short_t digitValue = digits->CurrentDigit();
 738 //      cout<<"Inner loop: sector, iRow, iColumn "
 739 //        <<sec<<" "<<iRow<<" "<<iColumn<<endl;
 740       if (digitValue >= zero) {
 741         Int_t label;
 742         for (Int_t j = 0; j<3; j++) {
 743           label = digits->GetTrackID(iRow,iColumn,j);
 744           if (label >= fNParticles) {
 745             cerr<<"particle label too big: fNParticles, label "
 746                 <<fNParticles<<" "<<label<<endl;
 747           }
 748           if (label >= 0 && label <= fNParticles) {
 749 //        if (label >= 0 && label <= fDebug) {
 750             if (fDebug > 6 ) {
 751               cout<<"Inner loop: sector, iRow, iColumn, label, value, row "
 752                   <<sec<<" "
 753                   <<iRow<<" "<<iColumn<<" "<<label<<" "<<digitValue
 754                   <<" "<<row<<endl;
 755             }
 756             fContainerDigitRow[label].SetRow(row+rowShift);
 757           }
 758         }
 759       }
 760     } while (digits->Next());
 761   }
 762
 763   if (fDebug > 2) cerr<<"end of TreeDLoop"<<endl;
 764
 765   return 0;
 766 }
 767
 768 ////////////////////////////////////////////////////////////////////////
 769 Int_t TPCFindGenTracks::TreeTRLoop()
 770 {
 771 /// loop over TrackReferences and store the first one for each track
 772
 773   TTree *treeTR=gAlice->TreeTR();
 774   if (!treeTR) {
 775     cerr<<"TreeTR not found"<<endl;
 776     return 1;
 777   }
 778   Int_t nPrimaries = (Int_t) treeTR->GetEntries();
 779   if (fDebug > 1) cout<<"There are "<<nPrimaries<<" entries in TreeTR"<<endl;
 780   TBranch *TPCBranchTR  = treeTR->GetBranch("TPC");
 781   if (!TPCBranchTR) {
 782     cerr<<"TPC branch in TR not found"<<endl;
 783     return 1;
 784   }
 785   TClonesArray* TPCArrayTR = new TClonesArray("AliTrackReference");
 786   TPCBranchTR->SetAddress(&TPCArrayTR);
 787   for (Int_t iPrimPart = 0; iPrimPart<nPrimaries; iPrimPart++) {
 788     TPCBranchTR->GetEntry(iPrimPart);
 789     for (Int_t iTrackRef = 0; iTrackRef < TPCArrayTR->GetEntriesFast(); iTrackRef++) {
 790       AliTrackReference *trackRef = (AliTrackReference*)TPCArrayTR->At(iTrackRef);
 791       Int_t label = trackRef->GetTrack();
 792
 793 // save info in the fContainerTR
 794       if (label<0  || label > fNParticles) {
 795         cerr<<"Wrong label: "<<label<<endl;
 796         continue;
 797 //      return 1;
 798       }
 799
 800 // store only if we do not have any track reference yet for this label
 801       if (fIndexTR[label] >= 0) continue;
 802
 803 // store only references with localX < fgTrackRefLocalXMax +- fgTrackRefLocalXMaxDelta
 804 // and the pT > fgPtCut
 805       Float_t localX = TR2LocalX(trackRef,fParamTPC);
 806       if (TMath::Abs(localX-fgTrackRefLocalXMax)>fgTrackRefLocalXMaxDelta) continue;
 807       if (trackRef->Pt() < fgPtCut) continue;
 808
 809
 810 //      cout<<"label, xg "<<label<<" "<<xg<<endl;
 811       if (fLastIndexTR >= fgMaxIndexTR-1) {
 812         cerr<<"Too many tracks with track reference. Increase the constant"
 813             <<" fgMaxIndexTR"<<endl;
 814         return 1;
 815       }
 816       fIndexTR[label] =  ++fLastIndexTR;
 817
 818 // someone was lazy to implement copy ctor in AliTrackReference, so we have to do
 819 //  it here "manually"
 820       fContainerTR[fIndexTR[label]].SetPosition(trackRef->X(),trackRef->Y(),trackRef->Z());
 821
 822       fContainerTR[fIndexTR[label]].SetMomentum(trackRef->Px(),trackRef->Py(),trackRef->Pz());
 823       fContainerTR[fIndexTR[label]].SetTrack(trackRef->GetTrack());
 824       fContainerTR[fIndexTR[label]].SetLength(trackRef->GetLength());
 825 //      cerr<<"Debug: trackRef->X(), stored: "<<trackRef->X()<<" "
 826 //        << fContainerTR[fIndexTR[label]].X()<<endl;
 827
 828     }
 829   }
 830   return 0;
 831 }
 832 ////////////////////////////////////////////////////////////////////////
 833 Float_t TPCFindGenTracks::TR2LocalX(AliTrackReference *trackRef,
 834                                     AliTPCParam *paramTPC) {
 835
 836   Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
 837   Int_t index[4];
 838   paramTPC->Transform0to1(x,index);
 839   paramTPC->Transform1to2(x,index);
 840   return x[0];
 841 }
 842 ////////////////////////////////////////////////////////////////////////
 843
 844
 845
 846 ////////////////////////////////////////////////////////////////////////
 847 //
 848 // Start of implementation of the class TPCCmpTr
 849 //
 850 ////////////////////////////////////////////////////////////////////////
 851
 852 class TPCCmpTr {
 853
 854 public:
 855   TPCCmpTr();
 856   TPCCmpTr(char* fnRecTracks,
 857            char* fnGenTracks   ="genTracks.root",
 858            char* fnCmpRes      ="cmpTracks.root",
 859            char* fnGalice      ="galice.root",
 860            Int_t nEvents=1, Int_t firstEvent=0);
 861   virtual ~TPCCmpTr();
 862   void Reset();
 863   Int_t Exec();
 864   Int_t Exec(Int_t nEvents, Int_t firstEventNr);
 865   void CreateTreeCmp();
 866   void CloseOutputFile();
 867   Bool_t ConnectGenTree();
 868   Int_t TreeGenLoop(Int_t eventNr);
 869   Int_t TreeTLoop(Int_t eventNr);
 870   void SetFirstEventNr(Int_t i) {fFirstEventNr = i;}
 871   void SetNEvents(Int_t i) {fNEvents = i;}
 872   void SetDebug(Int_t level) {fDebug = level;}
 873
 874 // tmp method, should go to TrackReferenceTPC
 875   Float_t TR2LocalX(AliTrackReference *trackRef,
 876                     AliTPCParam *paramTPC);
 877
 878 private:
 879   digitRow *fDigitRow;            //!< pointer to the object saved in Branch
 880   Int_t fEventNr;                 //!< current event number
 881   Int_t fLabel;                   //!< track label
 882   Int_t fNEvents;                 //!< number of events to process
 883   Int_t fFirstEventNr;            //!< first event to process
 884
 885   char *fFnCmp;                   //!< output file name with cmp tracks
 886   TFile *fFileCmp;                //!< output file with cmp tracks
 887   TTree *fTreeCmp;                //!< output tree with cmp tracks
 888
 889   char *fFnGenTracks;             //!< input file name with gen tracks
 890   TFile *fFileGenTracks;
 891   TTree *fTreeGenTracks;
 892
 893   char *fFnHits;                  //!< input file name with gAlice object (needed for B)
 894   TFile *fFileHits;               //!< input file with gAlice
 895
 896   char *fFnRecTracks;             //!< input file name with tpc rec. tracks
 897   TFile *fFileRecTracks;          //!< input file with reconstructed tracks
 898   TTree *fTreeRecTracks;          //!< tree with reconstructed tracks
 899
 900   AliTPCtrack *fTPCTrack;         //!< pointer to TPC track to connect branch
 901   Int_t *fIndexRecTracks;         //!< index of particle label in the TreeT_TPC
 902
 903   Int_t fRowsWithDigitsInn;       //!< number of rows with digits in the inner sectors
 904   Int_t fRowsWithDigits;          //!< number of rows with digits in the outer sectors
 905   Int_t fRowsTrackLength;         //!< last - first row with digit
 906   Int_t fDigitsInSeed;            //!< digits in the default seed rows
 907   TParticle *fParticle;           //!< generated particle
 908   Double_t fVDist[4];             //!< distance of the particle vertex from primary vertex
 909                                   // the fVDist[3] contains size of the 3-vector
 910   AliTrackReference *fTrackRef;   //!< track reference saved in the output tree
 911   Int_t   fReconstructed;         //!< flag if track was reconstructed
 912   AliTPCParam* fParamTPC;         //!< AliTPCParam
 913
 914   Int_t fNParticles;              //!< number of particles in the input tree genTracks
 915   Int_t fDebug;                   //!< debug flag
 916
 917   Int_t fNextTreeGenEntryToRead;    //!< last entry already read from genTracks tree
 918   TPCGenInfo *fGenInfo;           //!< container for all the details
 919
 920   Double_t fRecPhi;         ///< reconstructed phi angle (0;2*kPI)
 921   Double_t fLambda;         ///< reconstructed
 922   Double_t fRecPt_1;        ///< reconstructed
 923   Float_t fdEdx;           ///< reconstructed  dEdx
 924
 925
 926   ClassDef(TPCCmpTr,1)    // class which creates and fills tree with TPCGenTrack objects
 927 };
 928 ClassImp(TPCCmpTr)
 929
 930 ////////////////////////////////////////////////////////////////////////
 931 TPCCmpTr::TPCCmpTr()
 932 {
 933   Reset();
 934 }
 935
 936 ////////////////////////////////////////////////////////////////////////
 937 TPCCmpTr::TPCCmpTr(char* fnRecTracks,
 938                    char* fnGenTracks,
 939                    char* fnCmp,
 940                    char* fnGalice,
 941                    Int_t nEvents, Int_t firstEvent)
 942 {
 943   Reset();
 944   fFnRecTracks = fnRecTracks;
 945   fFnGenTracks = fnGenTracks;
 946   fFnCmp = fnCmp;
 947   fFnHits = fnGalice;
 948   fFirstEventNr = firstEvent;
 949   fEventNr = firstEvent;
 950   fNEvents = nEvents;
 951 }
 952 ////////////////////////////////////////////////////////////////////////
 953 void TPCCmpTr::Reset()
 954 {
 955   fDigitRow = 0;
 956   fEventNr = 0;
 957   fNEvents = 0;
 958   fTreeCmp = 0;
 959   fFnCmp = "cmpTracks.root";
 960   fFnHits = "galice.root";
 961   fFileGenTracks = 0;
 962   fFileHits =0;
 963   fParticle = 0;
 964   fTrackRef = 0;
 965
 966   fRowsWithDigitsInn = 0;
 967   fRowsWithDigits = 0;
 968   fRowsTrackLength = 0;
 969   fDigitsInSeed = 0;
 970
 971   fDebug = 0;
 972
 973   fParamTPC = 0;
 974   fFnRecTracks = "tpc.tracks.root";
 975   fTreeRecTracks = 0;
 976   fFileRecTracks = 0;
 977   fTPCTrack = 0;
 978   fGenInfo = 0;
 979 }
 980 ////////////////////////////////////////////////////////////////////////
 981 TPCCmpTr::~TPCCmpTr()
 982 {
 983   ;
 984 }
 985 ////////////////////////////////////////////////////////////////////////
 986 Int_t TPCCmpTr::Exec(Int_t nEvents, Int_t firstEventNr)
 987 {
 988   fNEvents = nEvents;
 989   fFirstEventNr = firstEventNr;
 990   return Exec();
 991 }
 992
 993 ////////////////////////////////////////////////////////////////////////
 994 Int_t TPCCmpTr::Exec()
 995 {
 996   TStopwatch timer;
 997   timer.Start();
 998
 999   fDigitRow = new digitRow();
1000   CreateTreeCmp();
1001   if (!fTreeCmp) {
1002     cerr<<"output tree not created"<<endl;
1003     return 1;
1004   }
1005   fFileHits = OpenAliceFile(fFnHits,kTRUE,"read"); //gAlice is read here
1006   if (!fFileHits) {
1007     cerr<<"Cannot open file with gAlice object "<<fFnHits<<endl;
1008     return 1;
1009   }
1010   fFileHits->cd();
1011   SetFieldFactor();
1012
1013   fParamTPC = LoadTPCParam(fFileHits);
1014   if (!fParamTPC) {
1015     cerr<<"TPC parameters not found and could not be created"<<endl;
1016     return 1;
1017   }
1018
1019   if (!ConnectGenTree()) {
1020     cerr<<"Cannot connect tree with generated tracks"<<endl;
1021     return 1;
1022   }
1023   fFileHits->cd();
1024   fNextTreeGenEntryToRead = 0;
1025   cerr<<"fFirstEventNr, fNEvents: "<<fFirstEventNr<<" "<<fNEvents<<endl;
1026   for (Int_t eventNr = fFirstEventNr; eventNr < fFirstEventNr+fNEvents;
1027        eventNr++) {
1028     fNParticles = gAlice->GetEvent(fEventNr);
1029     fIndexRecTracks = new Int_t[fNParticles];
1030     for (Int_t i = 0; i<fNParticles; i++) {
1031       fIndexRecTracks[i] = -1;
1032     }
1033
1034     cout<<"Start to process event "<<fEventNr<<endl;
1035     cout<<"\tfNParticles = "<<fNParticles<<endl;
1036     if (fDebug>2) cout<<"\tStart loop over TreeT"<<endl;
1037     if (TreeTLoop(eventNr)>0) return 1;
1038
1039     if (fDebug>2) cout<<"\tStart loop over tree genTracks"<<endl;
1040     if (TreeGenLoop(eventNr)>0) return 1;
1041     if (fDebug>2) cout<<"\tEnd loop over tree genTracks"<<endl;
1042
1043     delete fTreeRecTracks;
1044
1045     delete [] fIndexRecTracks;
1046   }
1047
1048   CloseOutputFile();
1049
1050   cerr<<"Exec finished"<<endl;
1051   fFileHits->cd();
1052   delete gAlice;
1053   fFileHits->Close();
1054   delete fFileHits;
1055
1056   timer.Stop();
1057   timer.Print();
1058   return 0;
1059 }
1060 ////////////////////////////////////////////////////////////////////////
1061 Bool_t TPCCmpTr::ConnectGenTree()
1062 {
1063 /// connect all branches from the genTracksTree
1064 /// use the same variables as for the new cmp tree, it may work
1065
1066   fFileGenTracks = TFile::Open(fFnGenTracks,"READ");
1067   if (!fFileGenTracks) {
1068     cerr<<"Error in ConnectGenTree: cannot open file "<<fFnGenTracks<<endl;
1069     return kFALSE;
1070   }
1071   fTreeGenTracks = (TTree*)fFileGenTracks->Get("genTracksTree");
1072   if (!fTreeGenTracks) {
1073     cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
1074         <<fFnGenTracks<<endl;
1075     return kFALSE;
1076   }
1077   fTreeGenTracks->SetBranchAddress("fEventNr",&fEventNr);
1078   fTreeGenTracks->SetBranchAddress("fLabel",&fLabel);
1079   fTreeGenTracks->SetBranchAddress("fRowsWithDigitsInn",&fRowsWithDigitsInn);
1080   fTreeGenTracks->SetBranchAddress("fRowsWithDigits",&fRowsWithDigits);
1081   fTreeGenTracks->SetBranchAddress("fRowsTrackLength",&fRowsTrackLength);
1082   fTreeGenTracks->SetBranchAddress("fDigitsInSeed",&fDigitsInSeed);
1083   fTreeGenTracks->SetBranchAddress("Particle",&fParticle);
1084   fTreeGenTracks->SetBranchAddress("fVDist",fVDist);
1085   fTreeGenTracks->SetBranchAddress("TR",&fTrackRef);
1086
1087   if (fDebug > 1) {
1088     cout<<"Number of gen. tracks with TR: "<<fTreeGenTracks->GetEntries()<<endl;
1089   }
1090   return kTRUE;
1091 }
1092
1093
1094 ////////////////////////////////////////////////////////////////////////
1095 void TPCCmpTr::CreateTreeCmp()
1096 {
1097   fFileCmp = TFile::Open(fFnCmp,"RECREATE");
1098   if (!fFileCmp) {
1099     cerr<<"Error in CreateTreeCmp: cannot open file "<<fFnCmp<<endl;
1100     return;
1101   }
1102   fTreeCmp = new TTree("TPCcmpTracks","TPCcmpTracks");
1103   TBranch *branchBits = fTreeCmp->Branch("bitsBranch", "digitRow", &fDigitRow, 4000, 0);
1104   if (!branchBits) {
1105     cerr<<"Error in CreateTreeCmp: cannot create branch."<<endl;
1106     return;
1107   }
1108   fTreeCmp->Branch("fEventNr",&fEventNr,"fEventNr/I");
1109   fTreeCmp->Branch("fLabel",&fLabel,"fLabel/I");
1110   fTreeCmp->Branch("fRowsWithDigitsInn",&fRowsWithDigitsInn,"fRowsWithDigitsInn/I");
1111   fTreeCmp->Branch("fRowsWithDigits",&fRowsWithDigits,"fRowsWithDigits/I");
1112   fTreeCmp->Branch("fRowsTrackLength",&fRowsTrackLength,"fRowsTrackLength/I");
1113   fTreeCmp->Branch("fDigitsInSeed",&fDigitsInSeed,"fDigitsInSeed/I");
1114
1115   fTreeCmp->Branch("fReconstructed",&fReconstructed,"fReconstructed/I");
1116   fTreeCmp->Branch("fTPCTrack","AliTPCtrack",&fTPCTrack);
1117
1118   fTreeCmp->Branch("Particle","TParticle",&fParticle);
1119   fTreeCmp->Branch("fVDist",&fVDist,"fVDist[4]/D");
1120   fTreeCmp->Branch("TR","AliTrackReference",&fTrackRef);
1121
1122   fTreeCmp->AutoSave();
1123 }
1124 ////////////////////////////////////////////////////////////////////////
1125 void TPCCmpTr::CloseOutputFile()
1126 {
1127   if (!fFileCmp) {
1128     cerr<<"File "<<fFnCmp<<" not found as an open file."<<endl;
1129     return;
1130   }
1131   fFileCmp->cd();
1132   fTreeCmp->Write();
1133   delete fTreeCmp;
1134   fFileCmp->Close();
1135   delete fFileCmp;
1136   return;
1137 }
1138 ////////////////////////////////////////////////////////////////////////
1139
1140 Float_t TPCCmpTr::TR2LocalX(AliTrackReference *trackRef,
1141                             AliTPCParam *paramTPC) {
1142
1143   Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
1144   Int_t index[4];
1145   paramTPC->Transform0to1(x,index);
1146   paramTPC->Transform1to2(x,index);
1147   return x[0];
1148 }
1149 ////////////////////////////////////////////////////////////////////////
1150
1151 Int_t TPCCmpTr::TreeTLoop(Int_t eventNr)
1152 {
1153 /// loop over all TPC reconstructed tracks and store info in memory
1154
1155
1156   if (!fFileRecTracks) fFileRecTracks = TFile::Open(fFnRecTracks,"read");
1157   if (!fFileRecTracks->IsOpen()) {
1158     cerr<<"Cannot open file "<<fFnRecTracks<<endl;
1159     return 1;
1160   }
1161
1162   char treeNameBase[11] = "TreeT_TPC_";
1163   char treeName[20];
1164   sprintf(treeName,"%s%d",treeNameBase,eventNr);
1165
1166   fTreeRecTracks=(TTree*)fFileRecTracks->Get(treeName);
1167   if (!fTreeRecTracks) {
1168     cerr<<"Can't get a tree with TPC rec. tracks named "<<treeName<<endl;
1169     return 1;
1170   }
1171
1172   Int_t nEntries = (Int_t) fTreeRecTracks->GetEntries();
1173   if (fDebug > 2) cout<<"Event, rec. tracks: "<<eventNr<<" "
1174                       <<nEntries<<endl;
1175   TBranch * br= fTreeRecTracks->GetBranch("tracks");
1176   br->SetAddress(&fTPCTrack);
1177
1178   for (Int_t iEntry=0; iEntry<nEntries;iEntry++){
1179     br->GetEntry(iEntry);
1180     Int_t label = fTPCTrack->GetLabel();
1181     fIndexRecTracks[label] =  iEntry;
1182   }
1183
1184   if (fDebug > 2) cerr<<"end of TreeTLoop"<<endl;
1185
1186   return 0;
1187 }
1188 ////////////////////////////////////////////////////////////////////////
1189 Int_t TPCCmpTr::TreeGenLoop(Int_t eventNr)
1190 {
1191 /// loop over all entries for a given event, find corresponding
1192 /// rec. track and store in the fTreeCmp
1193
1194   fFileGenTracks->cd();
1195   Int_t entry = fNextTreeGenEntryToRead;
1196   Double_t nParticlesTR = fTreeGenTracks->GetEntriesFast();
1197   cerr<<"fNParticles, nParticlesTR, fNextTreeGenEntryToRead: "<<fNParticles<<" "
1198       <<nParticlesTR<<" "<<fNextTreeGenEntryToRead<<endl;
1199   while (entry < nParticlesTR) {
1200     fTreeGenTracks->GetEntry(entry);
1201     entry++;
1202     if (fEventNr < eventNr) continue;
1203     if (fEventNr > eventNr) break;
1204     fNextTreeGenEntryToRead = entry-1;
1205     if (fDebug > 2 && fLabel < 10) {
1206       cerr<<"Fill track with a label "<<fLabel<<endl;
1207     }
1208
1209     fReconstructed = 0;
1210     fdEdx = 0.;
1211     fRecPhi = fLambda = fRecPt_1 = 0.;
1212
1213     if (fDebug > 2) {
1214       cerr<<"fLabel, fIndexRecTracks[fLabel] "<<fLabel<<" "<<fIndexRecTracks[fLabel]<<endl;
1215     }
1216     if (fIndexRecTracks[fLabel] >= 0) {
1217       Int_t nBytes = fTreeRecTracks->GetEvent(fIndexRecTracks[fLabel]);
1218       if (nBytes > 0) {
1219         fReconstructed = 1;
1220         fdEdx = fTPCTrack->GetdEdx();
1221         Double_t Localx = TR2LocalX(fTrackRef,fParamTPC);
1222         if (fDebug > 3) {
1223           cerr<<"Track local X before prolongation: "<<fTPCTrack->GetX()<<endl;
1224         }
1225         fTPCTrack->PropagateTo(Localx);
1226         Double_t par[5];
1227         if (fDebug > 3) {
1228           cerr<<"Track local X after prolongation: "<<fTPCTrack->GetX()<<endl;
1229         }
1230         fTPCTrack->GetExternalParameters(Localx,par);
1231         fRecPhi=TMath::ASin(par[2]) + fTPCTrack->GetAlpha();
1232         if (fRecPhi<0) fRecPhi+=2*TMath::Pi();
1233         if (fRecPhi>=2*TMath::Pi()) fRecPhi-=2*TMath::Pi();
1234 //        fRecPhi = (fRecPhi)*kRaddeg;
1235         fLambda = TMath::ATan(par[3]);
1236         fRecPt_1 = TMath::Abs(par[4]);
1237       }
1238     }
1239
1240
1241     fTreeCmp->Fill();
1242
1243   }
1244   fTreeCmp->AutoSave();
1245 //  delete gAlice; gAlice = 0;
1246 //  fFileHits->Close();
1247
1248   if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
1249
1250   return 0;
1251 }
1252 ////////////////////////////////////////////////////////////////////////