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