1 /// \file AliTPCCmpNG.C
5 /// define a class TPCGenTrack
6 /// save TPC related properties of tracks into a single tree
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
19 /// TPCFindGenTracks *t = new TPCFindGenTracks("galice.root","tpc.digits.root")
24 /// TPCCmpTr *t2 = new TPCCmpTr("tpc.tracks.root","genTracks.root","cmpTracks.root");
29 /// Step 1 - summurize information from simulation
31 /// Compile macro with ACLIC:
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:
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):
60 /// Step 2 - compare reconstructed tracks with simulated
62 /// Load (and compile) the macro:
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:
74 /// Step 3 - study the results
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")
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
96 /// \author Jiri Chudoba
99 #if !defined(__CINT__) || defined(__MAKECINT__)
100 #include "iostream.h"
107 #include "TBenchmark.h"
108 #include "TStopwatch.h"
109 #include "TParticle.h"
111 #include "AliStack.h"
112 #include "AliTPCtrack.h"
113 #include "AliSimDigits.h"
114 #include "AliTPCParam.h"
115 #include "TParticle.h"
117 #include "AliDetector.h"
118 #include "AliTrackReference.h"
121 // #include "AliConst.h"
124 #include "AliTPCTracking.C"
127 ////////////////////////////////////////////////////////////////////////
129 // name: ML.C (Macro Library - collection of functions used in diff macros
131 // last update: 08.05.2002
132 // author: Jiri Chudoba
138 // 08.05.02 - first version
140 ////////////////////////////////////////////////////////////////////////
142 #if !defined(__CINT__) || defined(__MAKECINT__)
143 #include "iostream.h"
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);
160 ////////////////////////////////////////////////////////////////////////
163 /// wait until user press return;
166 Bool_t done = kFALSE;
167 TTimer *timer = new TTimer("gSystem->ProcessEvents();", 50, kFALSE);
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: ");
176 if (input) done = kTRUE;
180 ////////////////////////////////////////////////////////////////////////
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
186 // open baseDir, loop over subdirs
188 void *dirp = gSystem->OpenDirectory(baseDir);
190 cerr<<"Could not open base directory "<<baseDir.Data()<<endl;
194 Long_t id, size, flag, modTime;
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);
207 gSystem->FreeDirectory(dirp);
208 // cerr<<"rc = "<<rc<<endl;
211 ////////////////////////////////////////////////////////////////////////
212 Bool_t ImportgAlice(TFile *file) {
213 /// read in gAlice object from the file
215 gAlice = (AliRun*)file->Get("gAlice");
216 if (!gAlice) return kFALSE;
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 "
227 if (!importgAlice) return file;
228 if (ImportgAlice(file)) return file;
231 ////////////////////////////////////////////////////////////////////////
234 ////////////////////////////////////////////////////////////////////////
236 // Start of implementation of the class TPCGenInfo
238 ////////////////////////////////////////////////////////////////////////
239 class TPCGenInfo: public TObject {
243 AliTrackReference *fTrackRef; ///< track reference saved in the output tree
244 TParticle *fParticle; ///< generated particle
245 Int_t fLabel; ///< track label
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
252 ClassDef(TPCGenInfo,1) // container for
255 ////////////////////////////////////////////////////////////////////////
256 TPCGenInfo::TPCGenInfo()
261 ////////////////////////////////////////////////////////////////////////
263 ////////////////////////////////////////////////////////////////////////
265 // End of implementation of the class TPCGenInfo
267 ////////////////////////////////////////////////////////////////////////
271 ////////////////////////////////////////////////////////////////////////
273 // Start of implementation of the class digitRow
275 ////////////////////////////////////////////////////////////////////////
276 const Int_t kgRowBytes = 32;
278 class digitRow: public TObject {
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);
293 UChar_t fDig[kgRowBytes];
295 ClassDef(digitRow,1) // container for digit pattern
298 ////////////////////////////////////////////////////////////////////////
303 ////////////////////////////////////////////////////////////////////////
304 digitRow & digitRow::operator=(const digitRow &digOld)
306 for (Int_t i = 0; i<kgRowBytes; i++) fDig[i] = digOld.fDig[i];
309 ////////////////////////////////////////////////////////////////////////
310 void digitRow::SetRow(Int_t row)
312 if (row >= 8*kgRowBytes) {
313 cerr<<"digitRow::SetRow: index "<<row<<" out of bounds."<<endl;
321 ////////////////////////////////////////////////////////////////////////
322 Bool_t digitRow::TestRow(Int_t row)
324 /// return kTRUE if row is on
328 return TESTBIT(fDig[iC],iB);
330 ////////////////////////////////////////////////////////////////////////
331 Int_t digitRow::RowsOn(Int_t upto)
333 /// returns number of rows with a digit
334 /// count only rows less equal row number upto
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++;
345 ////////////////////////////////////////////////////////////////////////
346 void digitRow::Reset()
348 /// resets all rows to zero
350 for (Int_t i = 0; i<kgRowBytes; i++) {
354 ////////////////////////////////////////////////////////////////////////
355 Int_t digitRow::Last()
357 /// returns the last row number with a digit
358 /// returns -1 if now digits
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;
367 ////////////////////////////////////////////////////////////////////////
368 Int_t digitRow::First()
370 /// returns the first row number with a digit
371 /// returns -1 if now digits
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;
381 ////////////////////////////////////////////////////////////////////////
383 // end of implementation of a class digitRow
385 ////////////////////////////////////////////////////////////////////////
386 ////////////////////////////////////////////////////////////////////////
388 // Start of implementation of the class TPCFindGenTracks
390 ////////////////////////////////////////////////////////////////////////
392 class 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();
403 Int_t Exec(Int_t nEvents, Int_t firstEventNr);
404 void CreateTreeGenTracks();
405 void CloseOutputFile();
409 void SetFirstEventNr(Int_t i) {fFirstEventNr = i;}
410 void SetNEvents(Int_t i) {fNEvents = i;}
411 void SetDebug(Int_t level) {fDebug = level;}
413 Float_t TR2LocalX(AliTrackReference *trackRef,
414 AliTPCParam *paramTPC);
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
437 AliTPCParam* fParamTPC; //!< AliTPCParam
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
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
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();
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.;
464 ClassDef(TPCFindGenTracks,1) // class which creates and fills tree with TPCGenTrack objects
466 ClassImp(TPCFindGenTracks)
468 ////////////////////////////////////////////////////////////////////////
469 TPCFindGenTracks::TPCFindGenTracks()
474 ////////////////////////////////////////////////////////////////////////
475 TPCFindGenTracks::TPCFindGenTracks(char* fnHits, char* fnDigits,
477 Int_t nEvents, Int_t firstEvent)
480 fFirstEventNr = firstEvent;
481 fEventNr = firstEvent;
485 fFnDigits = fnDigits;
487 ////////////////////////////////////////////////////////////////////////
488 void TPCFindGenTracks::Reset()
494 fFnRes = "genTracks.root";
495 fFnHits = "rfio:galice.root";
496 fFnDigits = "rfio:its.tpc.trd.digits.root";
500 fContainerDigitRow = 0;
513 ////////////////////////////////////////////////////////////////////////
514 TPCFindGenTracks::~TPCFindGenTracks()
518 ////////////////////////////////////////////////////////////////////////
519 Int_t TPCFindGenTracks::Exec(Int_t nEvents, Int_t firstEventNr)
522 fFirstEventNr = firstEventNr;
526 ////////////////////////////////////////////////////////////////////////
527 Int_t TPCFindGenTracks::Exec()
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;
540 fFileTreeD = TFile::Open(fFnDigits,"read");
541 if (!fFileTreeD->IsOpen()) {
542 cerr<<"Cannot open file "<<fFnDigits<<endl;
546 fParamTPC = LoadTPCParam(fFileTreeD);
548 cerr<<"TPC parameters not found and could not be created"<<endl;
552 for (fEventNr = fFirstEventNr; fEventNr < fFirstEventNr+fNEvents;
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++) {
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;
572 delete [] fContainerDigitRow;
573 delete [] fContainerTR;
579 cerr<<"Exec finished"<<endl;
589 ////////////////////////////////////////////////////////////////////////
590 void TPCFindGenTracks::CreateTreeGenTracks()
592 fFileGenTracks = TFile::Open(fFnRes,"RECREATE");
593 if (!fFileGenTracks) {
594 cerr<<"Error in CreateTreeGenTracks: cannot open file "<<fFnRes<<endl;
597 fTreeGenTracks = new TTree("genTracksTree","genTracksTree");
598 TBranch *branchBits = fTreeGenTracks->Branch("bitsBranch", "digitRow", &fDigitRow, 4000, 0);
600 cerr<<"Error in CreateTreeGenTracks: cannot create branch."<<endl;
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");
610 fTreeGenTracks->Branch("Particle","TParticle",&fParticle);
611 fTreeGenTracks->Branch("fVDist",&fVDist,"fVDist[4]/D");
612 fTreeGenTracks->Branch("TR","AliTrackReference",&fTrackRef);
614 fTreeGenTracks->AutoSave();
616 ////////////////////////////////////////////////////////////////////////
617 void TPCFindGenTracks::CloseOutputFile()
619 if (!fFileGenTracks) {
620 cerr<<"File "<<fFnRes<<" not found as an open file."<<endl;
623 fFileGenTracks->cd();
624 fTreeGenTracks->Write();
625 delete fTreeGenTracks;
626 fFileGenTracks->Close();
627 delete fFileGenTracks;
631 ////////////////////////////////////////////////////////////////////////
632 Int_t TPCFindGenTracks::TreeKLoop()
634 /// open the file with treeK
635 /// loop over all entries there and save information about some tracks
639 cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
642 AliStack * stack = gAlice->Stack();
643 if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
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();
652 for (Int_t iParticle = 0; iParticle < fNParticles; iParticle++) {
653 // for (Int_t iParticle = 0; iParticle < fDebug; iParticle++) {
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<<" ";
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;
678 if (fDigitRow->TestRow(seedRow11) && fDigitRow->TestRow(seedRow12))
680 if (fDigitRow->TestRow(seedRow21) && fDigitRow->TestRow(seedRow22))
683 if (fIndexTR[iParticle] >= 0) {
684 fTrackRef = &(fContainerTR[fIndexTR[iParticle]]);
685 // cerr<<"Debug: fTrackRef->X() = "<<fTrackRef->X()<<endl;
687 fTrackRef->SetTrack(-1);
690 fTreeGenTracks->Fill();
693 fTreeGenTracks->AutoSave();
694 // delete gAlice; gAlice = 0;
695 // fFileHits->Close();
697 if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
701 ////////////////////////////////////////////////////////////////////////
702 Int_t TPCFindGenTracks::TreeDLoop()
704 /// open the file with treeD
705 /// loop over all entries there and save information about some tracks
708 // Int_t nrow_up=fParamTPC->GetNRowUp();
709 // Int_t nrows=fParamTPC->GetNRowLow()+nrow_up;
710 Int_t nInnerSector = fParamTPC->GetNInnerSector();
712 Int_t zero=fParamTPC->GetZeroSup();
713 // Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
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);
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;
727 fParamTPC->AdjustSectorRow(digits->GetID(),sec,row);
728 if (fDebug > 1) cout<<sec<<' '<<row<<" \r";
729 // cerr<<sec<<' '<<row<<endl;
731 // here I expect that upper sectors follow lower sectors
732 if (sec > nInnerSector) rowShift = fParamTPC->GetNRowLow();
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) {
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;
748 if (label >= 0 && label <= fNParticles) {
749 // if (label >= 0 && label <= fDebug) {
751 cout<<"Inner loop: sector, iRow, iColumn, label, value, row "
753 <<iRow<<" "<<iColumn<<" "<<label<<" "<<digitValue
756 fContainerDigitRow[label].SetRow(row+rowShift);
760 } while (digits->Next());
763 if (fDebug > 2) cerr<<"end of TreeDLoop"<<endl;
768 ////////////////////////////////////////////////////////////////////////
769 Int_t TPCFindGenTracks::TreeTRLoop()
771 /// loop over TrackReferences and store the first one for each track
773 TTree *treeTR=gAlice->TreeTR();
775 cerr<<"TreeTR not found"<<endl;
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");
782 cerr<<"TPC branch in TR not found"<<endl;
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();
793 // save info in the fContainerTR
794 if (label<0 || label > fNParticles) {
795 cerr<<"Wrong label: "<<label<<endl;
800 // store only if we do not have any track reference yet for this label
801 if (fIndexTR[label] >= 0) continue;
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;
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;
816 fIndexTR[label] = ++fLastIndexTR;
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());
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;
832 ////////////////////////////////////////////////////////////////////////
833 Float_t TPCFindGenTracks::TR2LocalX(AliTrackReference *trackRef,
834 AliTPCParam *paramTPC) {
836 Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
838 paramTPC->Transform0to1(x,index);
839 paramTPC->Transform1to2(x,index);
842 ////////////////////////////////////////////////////////////////////////
846 ////////////////////////////////////////////////////////////////////////
848 // Start of implementation of the class TPCCmpTr
850 ////////////////////////////////////////////////////////////////////////
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);
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;}
874 // tmp method, should go to TrackReferenceTPC
875 Float_t TR2LocalX(AliTrackReference *trackRef,
876 AliTPCParam *paramTPC);
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
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
889 char *fFnGenTracks; //!< input file name with gen tracks
890 TFile *fFileGenTracks;
891 TTree *fTreeGenTracks;
893 char *fFnHits; //!< input file name with gAlice object (needed for B)
894 TFile *fFileHits; //!< input file with gAlice
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
900 AliTPCtrack *fTPCTrack; //!< pointer to TPC track to connect branch
901 Int_t *fIndexRecTracks; //!< index of particle label in the TreeT_TPC
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
914 Int_t fNParticles; //!< number of particles in the input tree genTracks
915 Int_t fDebug; //!< debug flag
917 Int_t fNextTreeGenEntryToRead; //!< last entry already read from genTracks tree
918 TPCGenInfo *fGenInfo; //!< container for all the details
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
926 ClassDef(TPCCmpTr,1) // class which creates and fills tree with TPCGenTrack objects
932 ////////////////////////////////////////////////////////////////////////
938 ////////////////////////////////////////////////////////////////////////
939 TPCCmpTr::TPCCmpTr(char* fnRecTracks,
943 Int_t nEvents, Int_t firstEvent)
946 fFnRecTracks = fnRecTracks;
947 fFnGenTracks = fnGenTracks;
950 fFirstEventNr = firstEvent;
951 fEventNr = firstEvent;
954 ////////////////////////////////////////////////////////////////////////
955 void TPCCmpTr::Reset()
961 fFnCmp = "cmpTracks.root";
962 fFnHits = "galice.root";
968 fRowsWithDigitsInn = 0;
970 fRowsTrackLength = 0;
976 fFnRecTracks = "tpc.tracks.root";
982 ////////////////////////////////////////////////////////////////////////
983 TPCCmpTr::~TPCCmpTr()
987 ////////////////////////////////////////////////////////////////////////
988 Int_t TPCCmpTr::Exec(Int_t nEvents, Int_t firstEventNr)
991 fFirstEventNr = firstEventNr;
995 ////////////////////////////////////////////////////////////////////////
996 Int_t TPCCmpTr::Exec()
1001 fDigitRow = new digitRow();
1004 cerr<<"output tree not created"<<endl;
1007 fFileHits = OpenAliceFile(fFnHits,kTRUE,"read"); //gAlice is read here
1009 cerr<<"Cannot open file with gAlice object "<<fFnHits<<endl;
1015 fParamTPC = LoadTPCParam(fFileHits);
1017 cerr<<"TPC parameters not found and could not be created"<<endl;
1021 if (!ConnectGenTree()) {
1022 cerr<<"Cannot connect tree with generated tracks"<<endl;
1026 fNextTreeGenEntryToRead = 0;
1027 cerr<<"fFirstEventNr, fNEvents: "<<fFirstEventNr<<" "<<fNEvents<<endl;
1028 for (Int_t eventNr = fFirstEventNr; eventNr < fFirstEventNr+fNEvents;
1030 fNParticles = gAlice->GetEvent(fEventNr);
1031 fIndexRecTracks = new Int_t[fNParticles];
1032 for (Int_t i = 0; i<fNParticles; i++) {
1033 fIndexRecTracks[i] = -1;
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;
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;
1045 delete fTreeRecTracks;
1047 delete [] fIndexRecTracks;
1052 cerr<<"Exec finished"<<endl;
1062 ////////////////////////////////////////////////////////////////////////
1063 Bool_t TPCCmpTr::ConnectGenTree()
1065 /// connect all branches from the genTracksTree
1066 /// use the same variables as for the new cmp tree, it may work
1068 fFileGenTracks = TFile::Open(fFnGenTracks,"READ");
1069 if (!fFileGenTracks) {
1070 cerr<<"Error in ConnectGenTree: cannot open file "<<fFnGenTracks<<endl;
1073 fTreeGenTracks = (TTree*)fFileGenTracks->Get("genTracksTree");
1074 if (!fTreeGenTracks) {
1075 cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
1076 <<fFnGenTracks<<endl;
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);
1090 cout<<"Number of gen. tracks with TR: "<<fTreeGenTracks->GetEntries()<<endl;
1096 ////////////////////////////////////////////////////////////////////////
1097 void TPCCmpTr::CreateTreeCmp()
1099 fFileCmp = TFile::Open(fFnCmp,"RECREATE");
1101 cerr<<"Error in CreateTreeCmp: cannot open file "<<fFnCmp<<endl;
1104 fTreeCmp = new TTree("TPCcmpTracks","TPCcmpTracks");
1105 TBranch *branchBits = fTreeCmp->Branch("bitsBranch", "digitRow", &fDigitRow, 4000, 0);
1107 cerr<<"Error in CreateTreeCmp: cannot create branch."<<endl;
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");
1117 fTreeCmp->Branch("fReconstructed",&fReconstructed,"fReconstructed/I");
1118 fTreeCmp->Branch("fTPCTrack","AliTPCtrack",&fTPCTrack);
1120 fTreeCmp->Branch("Particle","TParticle",&fParticle);
1121 fTreeCmp->Branch("fVDist",&fVDist,"fVDist[4]/D");
1122 fTreeCmp->Branch("TR","AliTrackReference",&fTrackRef);
1124 fTreeCmp->AutoSave();
1126 ////////////////////////////////////////////////////////////////////////
1127 void TPCCmpTr::CloseOutputFile()
1130 cerr<<"File "<<fFnCmp<<" not found as an open file."<<endl;
1140 ////////////////////////////////////////////////////////////////////////
1142 Float_t TPCCmpTr::TR2LocalX(AliTrackReference *trackRef,
1143 AliTPCParam *paramTPC) {
1145 Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
1147 paramTPC->Transform0to1(x,index);
1148 paramTPC->Transform1to2(x,index);
1151 ////////////////////////////////////////////////////////////////////////
1153 Int_t TPCCmpTr::TreeTLoop(Int_t eventNr)
1155 /// loop over all TPC reconstructed tracks and store info in memory
1158 if (!fFileRecTracks) fFileRecTracks = TFile::Open(fFnRecTracks,"read");
1159 if (!fFileRecTracks->IsOpen()) {
1160 cerr<<"Cannot open file "<<fFnRecTracks<<endl;
1164 char treeNameBase[11] = "TreeT_TPC_";
1166 sprintf(treeName,"%s%d",treeNameBase,eventNr);
1168 fTreeRecTracks=(TTree*)fFileRecTracks->Get(treeName);
1169 if (!fTreeRecTracks) {
1170 cerr<<"Can't get a tree with TPC rec. tracks named "<<treeName<<endl;
1174 Int_t nEntries = (Int_t) fTreeRecTracks->GetEntries();
1175 if (fDebug > 2) cout<<"Event, rec. tracks: "<<eventNr<<" "
1177 TBranch * br= fTreeRecTracks->GetBranch("tracks");
1178 br->SetAddress(&fTPCTrack);
1180 for (Int_t iEntry=0; iEntry<nEntries;iEntry++){
1181 br->GetEntry(iEntry);
1182 Int_t label = fTPCTrack->GetLabel();
1183 fIndexRecTracks[label] = iEntry;
1186 if (fDebug > 2) cerr<<"end of TreeTLoop"<<endl;
1190 ////////////////////////////////////////////////////////////////////////
1191 Int_t TPCCmpTr::TreeGenLoop(Int_t eventNr)
1193 /// loop over all entries for a given event, find corresponding
1194 /// rec. track and store in the fTreeCmp
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);
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;
1213 fRecPhi = fLambda = fRecPt_1 = 0.;
1216 cerr<<"fLabel, fIndexRecTracks[fLabel] "<<fLabel<<" "<<fIndexRecTracks[fLabel]<<endl;
1218 if (fIndexRecTracks[fLabel] >= 0) {
1219 Int_t nBytes = fTreeRecTracks->GetEvent(fIndexRecTracks[fLabel]);
1222 fdEdx = fTPCTrack->GetdEdx();
1223 Double_t Localx = TR2LocalX(fTrackRef,fParamTPC);
1225 cerr<<"Track local X before prolongation: "<<fTPCTrack->GetX()<<endl;
1227 fTPCTrack->PropagateTo(Localx);
1230 cerr<<"Track local X after prolongation: "<<fTPCTrack->GetX()<<endl;
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]);
1246 fTreeCmp->AutoSave();
1247 // delete gAlice; gAlice = 0;
1248 // fFileHits->Close();
1250 if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
1254 ////////////////////////////////////////////////////////////////////////