Update timestamps for new AMANDA simulation (17/02/2015)
[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 /// wait until user press return;
164
165   char    *input;
166   Bool_t done = kFALSE;
167   TTimer  *timer = new TTimer("gSystem->ProcessEvents();", 50, kFALSE);
168
169   do {
170     timer->TurnOn();
171     timer->Reset();
172     // Now let's read the input, we can use here any
173     // stdio or iostream reading methods. like std::cin >> myinputl;
174     input = Getline("Type <return> to continue: "); 
175     timer->TurnOff();
176     if (input) done = kTRUE;
177   } while (!done);
178 }
179
180 ////////////////////////////////////////////////////////////////////////
181
182 Int_t Chain(TString baseDir, TString subDirNameMask, TString fn, TChain& chain) {
183 /// chain all files fn in the subdirectories which match subDirName
184 /// return number of chained files
185
186 // open baseDir, loop over subdirs
187
188   void *dirp = gSystem->OpenDirectory(baseDir); 
189   if (!dirp) {
190     cerr<<"Could not open base directory "<<baseDir.Data()<<endl;
191     return 0;
192   }
193   const char *afile;
194   Long_t id, size, flag, modTime;
195   Int_t rc = 0;
196   char relName[100];
197   while((afile = gSystem->GetDirEntry(dirp))) {
198 //    printf("file: %s\n",afile);
199     if (strstr(afile,subDirNameMask.Data())) {
200       sprintf(relName,"%s/%s/%s",baseDir.Data(),afile,fn.Data());
201 //      cerr<<"relName = "<<relName<<endl;
202       if(!gSystem->GetPathInfo(relName, &id, &size, &flag, &modTime)) { 
203         rc += chain.Add(relName);      
204       }
205     }
206   }
207   gSystem->FreeDirectory(dirp);
208 //  cerr<<"rc = "<<rc<<endl;
209   return rc;
210 }
211 ////////////////////////////////////////////////////////////////////////
212 Bool_t ImportgAlice(TFile *file) {
213 /// read in gAlice object from the file
214
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 /// \cond CLASSIMP
929 ClassImp(TPCCmpTr)
930 /// \endcond
931   
932 ////////////////////////////////////////////////////////////////////////
933 TPCCmpTr::TPCCmpTr()
934 {
935   Reset();
936 }
937
938 ////////////////////////////////////////////////////////////////////////
939 TPCCmpTr::TPCCmpTr(char* fnRecTracks,
940                    char* fnGenTracks,
941                    char* fnCmp,
942                    char* fnGalice,
943                    Int_t nEvents, Int_t firstEvent)
944 {
945   Reset();
946   fFnRecTracks = fnRecTracks;
947   fFnGenTracks = fnGenTracks;
948   fFnCmp = fnCmp;
949   fFnHits = fnGalice;
950   fFirstEventNr = firstEvent;
951   fEventNr = firstEvent;
952   fNEvents = nEvents;
953 }
954 ////////////////////////////////////////////////////////////////////////
955 void TPCCmpTr::Reset()
956 {
957   fDigitRow = 0;
958   fEventNr = 0;
959   fNEvents = 0;
960   fTreeCmp = 0;
961   fFnCmp = "cmpTracks.root";
962   fFnHits = "galice.root";
963   fFileGenTracks = 0;
964   fFileHits =0;
965   fParticle = 0;
966   fTrackRef = 0;
967
968   fRowsWithDigitsInn = 0;
969   fRowsWithDigits = 0;
970   fRowsTrackLength = 0;
971   fDigitsInSeed = 0;
972
973   fDebug = 0;
974
975   fParamTPC = 0;
976   fFnRecTracks = "tpc.tracks.root";
977   fTreeRecTracks = 0;
978   fFileRecTracks = 0;
979   fTPCTrack = 0; 
980   fGenInfo = 0; 
981 }
982 ////////////////////////////////////////////////////////////////////////
983 TPCCmpTr::~TPCCmpTr()
984 {
985   ;
986 }
987 ////////////////////////////////////////////////////////////////////////
988 Int_t TPCCmpTr::Exec(Int_t nEvents, Int_t firstEventNr)
989 {
990   fNEvents = nEvents;
991   fFirstEventNr = firstEventNr;
992   return Exec();
993 }
994
995 ////////////////////////////////////////////////////////////////////////
996 Int_t TPCCmpTr::Exec()
997 {
998   TStopwatch timer;
999   timer.Start();
1000
1001   fDigitRow = new digitRow();
1002   CreateTreeCmp();
1003   if (!fTreeCmp) {
1004     cerr<<"output tree not created"<<endl;
1005     return 1;
1006   }
1007   fFileHits = OpenAliceFile(fFnHits,kTRUE,"read"); //gAlice is read here
1008   if (!fFileHits) {
1009     cerr<<"Cannot open file with gAlice object "<<fFnHits<<endl;
1010     return 1;  
1011   }
1012   fFileHits->cd();
1013   SetFieldFactor(); 
1014
1015   fParamTPC = LoadTPCParam(fFileHits);
1016   if (!fParamTPC) {
1017     cerr<<"TPC parameters not found and could not be created"<<endl;
1018     return 1;
1019   }
1020
1021   if (!ConnectGenTree()) {
1022     cerr<<"Cannot connect tree with generated tracks"<<endl;
1023     return 1;
1024   }
1025   fFileHits->cd();
1026   fNextTreeGenEntryToRead = 0;
1027   cerr<<"fFirstEventNr, fNEvents: "<<fFirstEventNr<<" "<<fNEvents<<endl;
1028   for (Int_t eventNr = fFirstEventNr; eventNr < fFirstEventNr+fNEvents;
1029        eventNr++) {
1030     fNParticles = gAlice->GetEvent(fEventNr);    
1031     fIndexRecTracks = new Int_t[fNParticles];
1032     for (Int_t i = 0; i<fNParticles; i++) {
1033       fIndexRecTracks[i] = -1;
1034     }
1035   
1036     cout<<"Start to process event "<<fEventNr<<endl;
1037     cout<<"\tfNParticles = "<<fNParticles<<endl;
1038     if (fDebug>2) cout<<"\tStart loop over TreeT"<<endl;
1039     if (TreeTLoop(eventNr)>0) return 1;
1040
1041     if (fDebug>2) cout<<"\tStart loop over tree genTracks"<<endl;
1042     if (TreeGenLoop(eventNr)>0) return 1;
1043     if (fDebug>2) cout<<"\tEnd loop over tree genTracks"<<endl;
1044
1045     delete fTreeRecTracks;
1046
1047     delete [] fIndexRecTracks;
1048   }
1049
1050   CloseOutputFile();
1051
1052   cerr<<"Exec finished"<<endl;
1053   fFileHits->cd();
1054   delete gAlice;
1055   fFileHits->Close();
1056   delete fFileHits;
1057
1058   timer.Stop();
1059   timer.Print();
1060   return 0;
1061 }
1062 ////////////////////////////////////////////////////////////////////////
1063 Bool_t TPCCmpTr::ConnectGenTree()
1064 {
1065 /// connect all branches from the genTracksTree
1066 /// use the same variables as for the new cmp tree, it may work
1067
1068   fFileGenTracks = TFile::Open(fFnGenTracks,"READ");
1069   if (!fFileGenTracks) {
1070     cerr<<"Error in ConnectGenTree: cannot open file "<<fFnGenTracks<<endl;
1071     return kFALSE;
1072   }
1073   fTreeGenTracks = (TTree*)fFileGenTracks->Get("genTracksTree");
1074   if (!fTreeGenTracks) {
1075     cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
1076         <<fFnGenTracks<<endl;
1077     return kFALSE;
1078   }
1079   fTreeGenTracks->SetBranchAddress("fEventNr",&fEventNr);
1080   fTreeGenTracks->SetBranchAddress("fLabel",&fLabel);
1081   fTreeGenTracks->SetBranchAddress("fRowsWithDigitsInn",&fRowsWithDigitsInn);
1082   fTreeGenTracks->SetBranchAddress("fRowsWithDigits",&fRowsWithDigits);
1083   fTreeGenTracks->SetBranchAddress("fRowsTrackLength",&fRowsTrackLength);
1084   fTreeGenTracks->SetBranchAddress("fDigitsInSeed",&fDigitsInSeed);
1085   fTreeGenTracks->SetBranchAddress("Particle",&fParticle);
1086   fTreeGenTracks->SetBranchAddress("fVDist",fVDist);
1087   fTreeGenTracks->SetBranchAddress("TR",&fTrackRef);
1088
1089   if (fDebug > 1) {
1090     cout<<"Number of gen. tracks with TR: "<<fTreeGenTracks->GetEntries()<<endl;
1091   }
1092   return kTRUE;
1093 }
1094
1095
1096 ////////////////////////////////////////////////////////////////////////
1097 void TPCCmpTr::CreateTreeCmp() 
1098 {
1099   fFileCmp = TFile::Open(fFnCmp,"RECREATE");
1100   if (!fFileCmp) {
1101     cerr<<"Error in CreateTreeCmp: cannot open file "<<fFnCmp<<endl;
1102     return;
1103   }
1104   fTreeCmp = new TTree("TPCcmpTracks","TPCcmpTracks");
1105   TBranch *branchBits = fTreeCmp->Branch("bitsBranch", "digitRow", &fDigitRow, 4000, 0);
1106   if (!branchBits) {
1107     cerr<<"Error in CreateTreeCmp: cannot create branch."<<endl;
1108     return;
1109   }
1110   fTreeCmp->Branch("fEventNr",&fEventNr,"fEventNr/I");
1111   fTreeCmp->Branch("fLabel",&fLabel,"fLabel/I");
1112   fTreeCmp->Branch("fRowsWithDigitsInn",&fRowsWithDigitsInn,"fRowsWithDigitsInn/I");
1113   fTreeCmp->Branch("fRowsWithDigits",&fRowsWithDigits,"fRowsWithDigits/I");
1114   fTreeCmp->Branch("fRowsTrackLength",&fRowsTrackLength,"fRowsTrackLength/I");
1115   fTreeCmp->Branch("fDigitsInSeed",&fDigitsInSeed,"fDigitsInSeed/I");
1116
1117   fTreeCmp->Branch("fReconstructed",&fReconstructed,"fReconstructed/I");
1118   fTreeCmp->Branch("fTPCTrack","AliTPCtrack",&fTPCTrack);
1119
1120   fTreeCmp->Branch("Particle","TParticle",&fParticle);
1121   fTreeCmp->Branch("fVDist",&fVDist,"fVDist[4]/D");
1122   fTreeCmp->Branch("TR","AliTrackReference",&fTrackRef);
1123
1124   fTreeCmp->AutoSave();
1125 }
1126 ////////////////////////////////////////////////////////////////////////
1127 void TPCCmpTr::CloseOutputFile() 
1128 {
1129   if (!fFileCmp) {
1130     cerr<<"File "<<fFnCmp<<" not found as an open file."<<endl;
1131     return;
1132   }
1133   fFileCmp->cd();
1134   fTreeCmp->Write();  
1135   delete fTreeCmp;
1136   fFileCmp->Close();
1137   delete fFileCmp;
1138   return;
1139 }
1140 ////////////////////////////////////////////////////////////////////////
1141
1142 Float_t TPCCmpTr::TR2LocalX(AliTrackReference *trackRef,
1143                             AliTPCParam *paramTPC) {
1144
1145   Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
1146   Int_t index[4];
1147   paramTPC->Transform0to1(x,index);
1148   paramTPC->Transform1to2(x,index);
1149   return x[0];
1150 }
1151 ////////////////////////////////////////////////////////////////////////
1152
1153 Int_t TPCCmpTr::TreeTLoop(Int_t eventNr)
1154 {
1155 /// loop over all TPC reconstructed tracks and store info in memory
1156
1157   
1158   if (!fFileRecTracks) fFileRecTracks = TFile::Open(fFnRecTracks,"read");
1159   if (!fFileRecTracks->IsOpen()) {
1160     cerr<<"Cannot open file "<<fFnRecTracks<<endl;
1161     return 1;
1162   }
1163
1164   char treeNameBase[11] = "TreeT_TPC_";
1165   char treeName[20];
1166   sprintf(treeName,"%s%d",treeNameBase,eventNr);
1167
1168   fTreeRecTracks=(TTree*)fFileRecTracks->Get(treeName);
1169   if (!fTreeRecTracks) {
1170     cerr<<"Can't get a tree with TPC rec. tracks named "<<treeName<<endl;
1171     return 1;
1172   }
1173   
1174   Int_t nEntries = (Int_t) fTreeRecTracks->GetEntries();
1175   if (fDebug > 2) cout<<"Event, rec. tracks: "<<eventNr<<" "
1176                       <<nEntries<<endl;
1177   TBranch * br= fTreeRecTracks->GetBranch("tracks");
1178   br->SetAddress(&fTPCTrack);
1179
1180   for (Int_t iEntry=0; iEntry<nEntries;iEntry++){
1181     br->GetEntry(iEntry);
1182     Int_t label = fTPCTrack->GetLabel();
1183     fIndexRecTracks[label] =  iEntry; 
1184   }  
1185
1186   if (fDebug > 2) cerr<<"end of TreeTLoop"<<endl;
1187
1188   return 0;
1189 }
1190 ////////////////////////////////////////////////////////////////////////
1191 Int_t TPCCmpTr::TreeGenLoop(Int_t eventNr)
1192 {
1193 /// loop over all entries for a given event, find corresponding
1194 /// rec. track and store in the fTreeCmp
1195
1196   fFileGenTracks->cd();
1197   Int_t entry = fNextTreeGenEntryToRead;
1198   Double_t nParticlesTR = fTreeGenTracks->GetEntriesFast();
1199   cerr<<"fNParticles, nParticlesTR, fNextTreeGenEntryToRead: "<<fNParticles<<" "
1200       <<nParticlesTR<<" "<<fNextTreeGenEntryToRead<<endl;
1201   while (entry < nParticlesTR) {
1202     fTreeGenTracks->GetEntry(entry);
1203     entry++;
1204     if (fEventNr < eventNr) continue;
1205     if (fEventNr > eventNr) break;
1206     fNextTreeGenEntryToRead = entry-1;
1207     if (fDebug > 2 && fLabel < 10) {
1208       cerr<<"Fill track with a label "<<fLabel<<endl;
1209     }
1210
1211     fReconstructed = 0;
1212     fdEdx = 0.;
1213     fRecPhi = fLambda = fRecPt_1 = 0.;
1214
1215     if (fDebug > 2) {
1216       cerr<<"fLabel, fIndexRecTracks[fLabel] "<<fLabel<<" "<<fIndexRecTracks[fLabel]<<endl;
1217     }
1218     if (fIndexRecTracks[fLabel] >= 0) {
1219       Int_t nBytes = fTreeRecTracks->GetEvent(fIndexRecTracks[fLabel]);
1220       if (nBytes > 0) {
1221         fReconstructed = 1;
1222         fdEdx = fTPCTrack->GetdEdx();
1223         Double_t Localx = TR2LocalX(fTrackRef,fParamTPC);
1224         if (fDebug > 3) {
1225           cerr<<"Track local X before prolongation: "<<fTPCTrack->GetX()<<endl;
1226         }
1227         fTPCTrack->PropagateTo(Localx);
1228         Double_t par[5];
1229         if (fDebug > 3) {
1230           cerr<<"Track local X after prolongation: "<<fTPCTrack->GetX()<<endl;
1231         }
1232         fTPCTrack->GetExternalParameters(Localx,par);
1233         fRecPhi=TMath::ASin(par[2]) + fTPCTrack->GetAlpha();
1234         if (fRecPhi<0) fRecPhi+=2*TMath::Pi();
1235         if (fRecPhi>=2*TMath::Pi()) fRecPhi-=2*TMath::Pi();
1236 //        fRecPhi = (fRecPhi)*kRaddeg;
1237         fLambda = TMath::ATan(par[3]);
1238         fRecPt_1 = TMath::Abs(par[4]);
1239       }
1240     }
1241
1242
1243     fTreeCmp->Fill();
1244
1245   }
1246   fTreeCmp->AutoSave();
1247 //  delete gAlice; gAlice = 0;
1248 //  fFileHits->Close();
1249
1250   if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
1251
1252   return 0;
1253 }
1254 ////////////////////////////////////////////////////////////////////////