TPC/AliTPCCmpNG.C

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