08ce290ca8423ceb48f4d3080e64c3a6bdb0f441
[u/mrichter/AliRoot.git] / 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 ////////////////////////////////////////////////////////////////////////