Update timestamps for new AMANDA simulation (17/02/2015)
[u/mrichter/AliRoot.git] / TPC / AliTPCCmpNG.C
CommitLineData
a7a1dd76 1/// \file AliTPCCmpNG.C
3a4edebe 2///
a7a1dd76 3/// version: 1.0
4/// description:
5/// define a class TPCGenTrack
6/// save TPC related properties of tracks into a single tree
3a4edebe 7///
a7a1dd76 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
3a4edebe 15///
a7a1dd76 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();
3a4edebe 26///
a7a1dd76 27/// Details:
3a4edebe 28///
a7a1dd76 29/// Step 1 - summurize information from simulation
3a4edebe 30///
a7a1dd76 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):
3a4edebe 59///
a7a1dd76 60/// Step 2 - compare reconstructed tracks with simulated
3a4edebe 61///
a7a1dd76 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();
3a4edebe 73///
a7a1dd76 74/// Step 3 - study the results
3a4edebe 75///
a7a1dd76 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")
3a4edebe 81///
a7a1dd76 82/// History:
3a4edebe 83///
a7a1dd76 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
3a4edebe 95///
a7a1dd76 96/// \author Jiri Chudoba
97/// \date 24.09.2002
e8558587 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
153void WaitForReturn();
154Bool_t ImportgAlice(TFile *file);
155TFile* OpenAliceFile(char *fn, Bool_t importgAlice = kFALSE, char *mode = "read");
156Int_t Chain(TString baseDir, TString subDirNameMask, TString fn, TChain& chain);
157
158#endif
159
160////////////////////////////////////////////////////////////////////////
161void WaitForReturn()
162{
3a4edebe 163/// wait until user press return;
164
e8558587 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
182Int_t Chain(TString baseDir, TString subDirNameMask, TString fn, TChain& chain) {
3a4edebe 183/// chain all files fn in the subdirectories which match subDirName
184/// return number of chained files
e8558587 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////////////////////////////////////////////////////////////////////////
212Bool_t ImportgAlice(TFile *file) {
3a4edebe 213/// read in gAlice object from the file
214
e8558587 215 gAlice = (AliRun*)file->Get("gAlice");
216 if (!gAlice) return kFALSE;
217 return kTRUE;
218}
219////////////////////////////////////////////////////////////////////////
220TFile* 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////////////////////////////////////////////////////////////////////////
239class TPCGenInfo: public TObject {
240
241public:
242 TPCGenInfo();
a7a1dd76 243 AliTrackReference *fTrackRef; ///< track reference saved in the output tree
244 TParticle *fParticle; ///< generated particle
245 Int_t fLabel; ///< track label
e8558587 246
a7a1dd76 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
e8558587 251
252 ClassDef(TPCGenInfo,1) // container for
253};
254ClassImp(TPCGenInfo)
255////////////////////////////////////////////////////////////////////////
256TPCGenInfo::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////////////////////////////////////////////////////////////////////////
276const Int_t kgRowBytes = 32;
277
278class digitRow: public TObject {
279
280public:
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};
297ClassImp(digitRow)
298////////////////////////////////////////////////////////////////////////
299digitRow::digitRow()
300{
301 Reset();
302}
303////////////////////////////////////////////////////////////////////////
304digitRow & 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////////////////////////////////////////////////////////////////////////
310void 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////////////////////////////////////////////////////////////////////////
322Bool_t digitRow::TestRow(Int_t row)
323{
a7a1dd76 324/// return kTRUE if row is on
325
e8558587 326 Int_t iC = row/8;
327 Int_t iB = row%8;
328 return TESTBIT(fDig[iC],iB);
329}
330////////////////////////////////////////////////////////////////////////
331Int_t digitRow::RowsOn(Int_t upto)
332{
a7a1dd76 333/// returns number of rows with a digit
334/// count only rows less equal row number upto
335
e8558587 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////////////////////////////////////////////////////////////////////////
346void digitRow::Reset()
347{
a7a1dd76 348/// resets all rows to zero
349
e8558587 350 for (Int_t i = 0; i<kgRowBytes; i++) {
351 fDig[i] <<= 8;
352 }
353}
354////////////////////////////////////////////////////////////////////////
355Int_t digitRow::Last()
356{
a7a1dd76 357/// returns the last row number with a digit
358/// returns -1 if now digits
359
e8558587 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////////////////////////////////////////////////////////////////////////
368Int_t digitRow::First()
369{
a7a1dd76 370/// returns the first row number with a digit
371/// returns -1 if now digits
372
e8558587 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
392class TPCFindGenTracks {
393
394public:
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
416public:
a7a1dd76 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
e8558587 441 // the fVDist[3] contains size of the 3-vector
a7a1dd76 442 TParticle *fParticle; //!< generated particle
e8558587 443
a7a1dd76 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
e8558587 448
449private:
450
451// some constants for the original non-pareller tracking (by Y.Belikov)
a7a1dd76 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
d8d3b5b8 456 static const Double_t kRaddeg = 180./TMath::Pi();
e8558587 457
a7a1dd76 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
e8558587 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};
466ClassImp(TPCFindGenTracks)
467
468////////////////////////////////////////////////////////////////////////
469TPCFindGenTracks::TPCFindGenTracks()
470{
471 Reset();
472}
473
474////////////////////////////////////////////////////////////////////////
475TPCFindGenTracks::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////////////////////////////////////////////////////////////////////////
488void 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////////////////////////////////////////////////////////////////////////
514TPCFindGenTracks::~TPCFindGenTracks()
515{
516 ;
517}
518////////////////////////////////////////////////////////////////////////
519Int_t TPCFindGenTracks::Exec(Int_t nEvents, Int_t firstEventNr)
520{
521 fNEvents = nEvents;
522 fFirstEventNr = firstEventNr;
523 return Exec();
524}
525
526////////////////////////////////////////////////////////////////////////
527Int_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////////////////////////////////////////////////////////////////////////
590void 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////////////////////////////////////////////////////////////////////////
617void 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////////////////////////////////////////////////////////////////////////
632Int_t TPCFindGenTracks::TreeKLoop()
633{
a7a1dd76 634/// open the file with treeK
635/// loop over all entries there and save information about some tracks
636
e8558587 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////////////////////////////////////////////////////////////////////////
702Int_t TPCFindGenTracks::TreeDLoop()
703{
a7a1dd76 704/// open the file with treeD
705/// loop over all entries there and save information about some tracks
e8558587 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////////////////////////////////////////////////////////////////////////
769Int_t TPCFindGenTracks::TreeTRLoop()
770{
a7a1dd76 771/// loop over TrackReferences and store the first one for each track
772
e8558587 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////////////////////////////////////////////////////////////////////////
833Float_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
852class TPCCmpTr {
853
854public:
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
878private:
a7a1dd76 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
e8558587 884
a7a1dd76 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
e8558587 888
a7a1dd76 889 char *fFnGenTracks; //!< input file name with gen tracks
e8558587 890 TFile *fFileGenTracks;
891 TTree *fTreeGenTracks;
892
a7a1dd76 893 char *fFnHits; //!< input file name with gAlice object (needed for B)
894 TFile *fFileHits; //!< input file with gAlice
e8558587 895
a7a1dd76 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
e8558587 899
a7a1dd76 900 AliTPCtrack *fTPCTrack; //!< pointer to TPC track to connect branch
901 Int_t *fIndexRecTracks; //!< index of particle label in the TreeT_TPC
e8558587 902
a7a1dd76 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
e8558587 909 // the fVDist[3] contains size of the 3-vector
a7a1dd76 910 AliTrackReference *fTrackRef; //!< track reference saved in the output tree
911 Int_t fReconstructed; //!< flag if track was reconstructed
912 AliTPCParam* fParamTPC; //!< AliTPCParam
e8558587 913
a7a1dd76 914 Int_t fNParticles; //!< number of particles in the input tree genTracks
915 Int_t fDebug; //!< debug flag
e8558587 916
a7a1dd76 917 Int_t fNextTreeGenEntryToRead; //!< last entry already read from genTracks tree
918 TPCGenInfo *fGenInfo; //!< container for all the details
e8558587 919
a7a1dd76 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
e8558587 924
925
926 ClassDef(TPCCmpTr,1) // class which creates and fills tree with TPCGenTrack objects
927};
3a4edebe 928/// \cond CLASSIMP
e8558587 929ClassImp(TPCCmpTr)
3a4edebe 930/// \endcond
e8558587 931
932////////////////////////////////////////////////////////////////////////
933TPCCmpTr::TPCCmpTr()
934{
935 Reset();
936}
937
938////////////////////////////////////////////////////////////////////////
939TPCCmpTr::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////////////////////////////////////////////////////////////////////////
955void 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////////////////////////////////////////////////////////////////////////
983TPCCmpTr::~TPCCmpTr()
984{
985 ;
986}
987////////////////////////////////////////////////////////////////////////
988Int_t TPCCmpTr::Exec(Int_t nEvents, Int_t firstEventNr)
989{
990 fNEvents = nEvents;
991 fFirstEventNr = firstEventNr;
992 return Exec();
993}
994
995////////////////////////////////////////////////////////////////////////
996Int_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////////////////////////////////////////////////////////////////////////
1063Bool_t TPCCmpTr::ConnectGenTree()
1064{
a7a1dd76 1065/// connect all branches from the genTracksTree
1066/// use the same variables as for the new cmp tree, it may work
1067
e8558587 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////////////////////////////////////////////////////////////////////////
1097void 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////////////////////////////////////////////////////////////////////////
1127void 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
1142Float_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
1153Int_t TPCCmpTr::TreeTLoop(Int_t eventNr)
1154{
a7a1dd76 1155/// loop over all TPC reconstructed tracks and store info in memory
e8558587 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////////////////////////////////////////////////////////////////////////
1191Int_t TPCCmpTr::TreeGenLoop(Int_t eventNr)
1192{
a7a1dd76 1193/// loop over all entries for a given event, find corresponding
1194/// rec. track and store in the fTreeCmp
1195
e8558587 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();
d8d3b5b8 1234 if (fRecPhi<0) fRecPhi+=2*TMath::Pi();
1235 if (fRecPhi>=2*TMath::Pi()) fRecPhi-=2*TMath::Pi();
e8558587 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////////////////////////////////////////////////////////////////////////