1 ////////////////////////////////////////////////////////////////////////
6 // author: Jiri Chudoba
9 // define a class TPCGenTrack
10 // save TPC related properties of tracks into a single tree
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
23 // TPCFindGenTracks *t = new TPCFindGenTracks("galice.root","tpc.digits.root")
28 // TPCCmpTr *t2 = new TPCCmpTr("tpc.tracks.root","genTracks.root","cmpTracks.root");
34 // Step 1 - summurize information from simulation
35 // ===============================================
36 // Compile macro with ACLIC:
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:
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):
65 // Step 2 - compare reconstructed tracks with simulated
66 // ====================================================
68 // Load (and compile) the macro:
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:
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")
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
99 // 03.02.03 - rename to AliTPCCmpNG.C, remove the part with rec. tracks
100 // (will be addded in a macro AliTPCCmpTr.C
103 ////////////////////////////////////////////////////////////////////////
105 #if !defined(__CINT__) || defined(__MAKECINT__)
106 #include "iostream.h"
113 #include "TBenchmark.h"
114 #include "TStopwatch.h"
115 #include "TParticle.h"
117 #include "AliStack.h"
118 #include "AliTPCtrack.h"
119 #include "AliSimDigits.h"
120 #include "AliTPCParam.h"
121 #include "TParticle.h"
123 #include "AliDetector.h"
124 #include "AliTrackReference.h"
127 // #include "AliConst.h"
130 #include "AliTPCTracking.C"
133 ////////////////////////////////////////////////////////////////////////
135 // name: ML.C (Macro Library - collection of functions used in diff macros
137 // last update: 08.05.2002
138 // author: Jiri Chudoba
144 // 08.05.02 - first version
146 ////////////////////////////////////////////////////////////////////////
148 #if !defined(__CINT__) || defined(__MAKECINT__)
149 #include "iostream.h"
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);
166 ////////////////////////////////////////////////////////////////////////
170 // wait until user press return;
173 Bool_t done = kFALSE;
174 TTimer *timer = new TTimer("gSystem->ProcessEvents();", 50, kFALSE);
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: ");
183 if (input) done = kTRUE;
187 ////////////////////////////////////////////////////////////////////////
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
193 // open baseDir, loop over subdirs
195 void *dirp = gSystem->OpenDirectory(baseDir);
197 cerr<<"Could not open base directory "<<baseDir.Data()<<endl;
201 Long_t id, size, flag, modTime;
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);
214 gSystem->FreeDirectory(dirp);
215 // cerr<<"rc = "<<rc<<endl;
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;
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 "
233 if (!importgAlice) return file;
234 if (ImportgAlice(file)) return file;
237 ////////////////////////////////////////////////////////////////////////
240 ////////////////////////////////////////////////////////////////////////
242 // Start of implementation of the class TPCGenInfo
244 ////////////////////////////////////////////////////////////////////////
245 class TPCGenInfo: public TObject {
249 AliTrackReference *fTrackRef; // track reference saved in the output tree
250 TParticle *fParticle; // generated particle
251 Int_t fLabel; // track label
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
258 ClassDef(TPCGenInfo,1) // container for
261 ////////////////////////////////////////////////////////////////////////
262 TPCGenInfo::TPCGenInfo()
267 ////////////////////////////////////////////////////////////////////////
269 ////////////////////////////////////////////////////////////////////////
271 // End of implementation of the class TPCGenInfo
273 ////////////////////////////////////////////////////////////////////////
277 ////////////////////////////////////////////////////////////////////////
279 // Start of implementation of the class digitRow
281 ////////////////////////////////////////////////////////////////////////
282 const Int_t kgRowBytes = 32;
284 class digitRow: public TObject {
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);
299 UChar_t fDig[kgRowBytes];
301 ClassDef(digitRow,1) // container for digit pattern
304 ////////////////////////////////////////////////////////////////////////
309 ////////////////////////////////////////////////////////////////////////
310 digitRow & digitRow::operator=(const digitRow &digOld)
312 for (Int_t i = 0; i<kgRowBytes; i++) fDig[i] = digOld.fDig[i];
315 ////////////////////////////////////////////////////////////////////////
316 void digitRow::SetRow(Int_t row)
318 if (row >= 8*kgRowBytes) {
319 cerr<<"digitRow::SetRow: index "<<row<<" out of bounds."<<endl;
327 ////////////////////////////////////////////////////////////////////////
328 Bool_t digitRow::TestRow(Int_t row)
331 // return kTRUE if row is on
335 return TESTBIT(fDig[iC],iB);
337 ////////////////////////////////////////////////////////////////////////
338 Int_t digitRow::RowsOn(Int_t upto)
341 // returns number of rows with a digit
342 // count only rows less equal row number upto
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++;
353 ////////////////////////////////////////////////////////////////////////
354 void digitRow::Reset()
357 // resets all rows to zero
359 for (Int_t i = 0; i<kgRowBytes; i++) {
363 ////////////////////////////////////////////////////////////////////////
364 Int_t digitRow::Last()
367 // returns the last row number with a digit
368 // returns -1 if now digits
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;
377 ////////////////////////////////////////////////////////////////////////
378 Int_t digitRow::First()
381 // returns the first row number with a digit
382 // returns -1 if now digits
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;
392 ////////////////////////////////////////////////////////////////////////
394 // end of implementation of a class digitRow
396 ////////////////////////////////////////////////////////////////////////
397 ////////////////////////////////////////////////////////////////////////
399 // Start of implementation of the class TPCFindGenTracks
401 ////////////////////////////////////////////////////////////////////////
403 class 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();
414 Int_t Exec(Int_t nEvents, Int_t firstEventNr);
415 void CreateTreeGenTracks();
416 void CloseOutputFile();
420 void SetFirstEventNr(Int_t i) {fFirstEventNr = i;}
421 void SetNEvents(Int_t i) {fNEvents = i;}
422 void SetDebug(Int_t level) {fDebug = level;}
424 Float_t TR2LocalX(AliTrackReference *trackRef,
425 AliTPCParam *paramTPC);
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
448 AliTPCParam* fParamTPC; //! AliTPCParam
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
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
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();
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.;
475 ClassDef(TPCFindGenTracks,1) // class which creates and fills tree with TPCGenTrack objects
477 ClassImp(TPCFindGenTracks)
479 ////////////////////////////////////////////////////////////////////////
480 TPCFindGenTracks::TPCFindGenTracks()
485 ////////////////////////////////////////////////////////////////////////
486 TPCFindGenTracks::TPCFindGenTracks(char* fnHits, char* fnDigits,
488 Int_t nEvents, Int_t firstEvent)
491 fFirstEventNr = firstEvent;
492 fEventNr = firstEvent;
496 fFnDigits = fnDigits;
498 ////////////////////////////////////////////////////////////////////////
499 void TPCFindGenTracks::Reset()
505 fFnRes = "genTracks.root";
506 fFnHits = "rfio:galice.root";
507 fFnDigits = "rfio:its.tpc.trd.digits.root";
511 fContainerDigitRow = 0;
524 ////////////////////////////////////////////////////////////////////////
525 TPCFindGenTracks::~TPCFindGenTracks()
529 ////////////////////////////////////////////////////////////////////////
530 Int_t TPCFindGenTracks::Exec(Int_t nEvents, Int_t firstEventNr)
533 fFirstEventNr = firstEventNr;
537 ////////////////////////////////////////////////////////////////////////
538 Int_t TPCFindGenTracks::Exec()
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;
551 fFileTreeD = TFile::Open(fFnDigits,"read");
552 if (!fFileTreeD->IsOpen()) {
553 cerr<<"Cannot open file "<<fFnDigits<<endl;
557 fParamTPC = LoadTPCParam(fFileTreeD);
559 cerr<<"TPC parameters not found and could not be created"<<endl;
563 for (fEventNr = fFirstEventNr; fEventNr < fFirstEventNr+fNEvents;
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++) {
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;
583 delete [] fContainerDigitRow;
584 delete [] fContainerTR;
590 cerr<<"Exec finished"<<endl;
600 ////////////////////////////////////////////////////////////////////////
601 void TPCFindGenTracks::CreateTreeGenTracks()
603 fFileGenTracks = TFile::Open(fFnRes,"RECREATE");
604 if (!fFileGenTracks) {
605 cerr<<"Error in CreateTreeGenTracks: cannot open file "<<fFnRes<<endl;
608 fTreeGenTracks = new TTree("genTracksTree","genTracksTree");
609 TBranch *branchBits = fTreeGenTracks->Branch("bitsBranch", "digitRow", &fDigitRow, 4000, 0);
611 cerr<<"Error in CreateTreeGenTracks: cannot create branch."<<endl;
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");
621 fTreeGenTracks->Branch("Particle","TParticle",&fParticle);
622 fTreeGenTracks->Branch("fVDist",&fVDist,"fVDist[4]/D");
623 fTreeGenTracks->Branch("TR","AliTrackReference",&fTrackRef);
625 fTreeGenTracks->AutoSave();
627 ////////////////////////////////////////////////////////////////////////
628 void TPCFindGenTracks::CloseOutputFile()
630 if (!fFileGenTracks) {
631 cerr<<"File "<<fFnRes<<" not found as an open file."<<endl;
634 fFileGenTracks->cd();
635 fTreeGenTracks->Write();
636 delete fTreeGenTracks;
637 fFileGenTracks->Close();
638 delete fFileGenTracks;
642 ////////////////////////////////////////////////////////////////////////
643 Int_t TPCFindGenTracks::TreeKLoop()
646 // open the file with treeK
647 // loop over all entries there and save information about some tracks
652 cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
655 AliStack * stack = gAlice->Stack();
656 if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
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();
665 for (Int_t iParticle = 0; iParticle < fNParticles; iParticle++) {
666 // for (Int_t iParticle = 0; iParticle < fDebug; iParticle++) {
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<<" ";
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;
691 if (fDigitRow->TestRow(seedRow11) && fDigitRow->TestRow(seedRow12))
693 if (fDigitRow->TestRow(seedRow21) && fDigitRow->TestRow(seedRow22))
696 if (fIndexTR[iParticle] >= 0) {
697 fTrackRef = &(fContainerTR[fIndexTR[iParticle]]);
698 // cerr<<"Debug: fTrackRef->X() = "<<fTrackRef->X()<<endl;
700 fTrackRef->SetTrack(-1);
703 fTreeGenTracks->Fill();
706 fTreeGenTracks->AutoSave();
707 // delete gAlice; gAlice = 0;
708 // fFileHits->Close();
710 if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
714 ////////////////////////////////////////////////////////////////////////
715 Int_t TPCFindGenTracks::TreeDLoop()
718 // open the file with treeD
719 // loop over all entries there and save information about some tracks
723 // Int_t nrow_up=fParamTPC->GetNRowUp();
724 // Int_t nrows=fParamTPC->GetNRowLow()+nrow_up;
725 Int_t nInnerSector = fParamTPC->GetNInnerSector();
727 Int_t zero=fParamTPC->GetZeroSup();
728 // Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
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);
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;
742 fParamTPC->AdjustSectorRow(digits->GetID(),sec,row);
743 if (fDebug > 1) cout<<sec<<' '<<row<<" \r";
744 // cerr<<sec<<' '<<row<<endl;
746 // here I expect that upper sectors follow lower sectors
747 if (sec > nInnerSector) rowShift = fParamTPC->GetNRowLow();
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) {
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;
763 if (label >= 0 && label <= fNParticles) {
764 // if (label >= 0 && label <= fDebug) {
766 cout<<"Inner loop: sector, iRow, iColumn, label, value, row "
768 <<iRow<<" "<<iColumn<<" "<<label<<" "<<digitValue
771 fContainerDigitRow[label].SetRow(row+rowShift);
775 } while (digits->Next());
778 if (fDebug > 2) cerr<<"end of TreeDLoop"<<endl;
783 ////////////////////////////////////////////////////////////////////////
784 Int_t TPCFindGenTracks::TreeTRLoop()
787 // loop over TrackReferences and store the first one for each track
790 TTree *treeTR=gAlice->TreeTR();
792 cerr<<"TreeTR not found"<<endl;
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");
799 cerr<<"TPC branch in TR not found"<<endl;
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();
810 // save info in the fContainerTR
811 if (label<0 || label > fNParticles) {
812 cerr<<"Wrong label: "<<label<<endl;
817 // store only if we do not have any track reference yet for this label
818 if (fIndexTR[label] >= 0) continue;
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;
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;
833 fIndexTR[label] = ++fLastIndexTR;
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());
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;
849 ////////////////////////////////////////////////////////////////////////
850 Float_t TPCFindGenTracks::TR2LocalX(AliTrackReference *trackRef,
851 AliTPCParam *paramTPC) {
853 Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
855 paramTPC->Transform0to1(x,index);
856 paramTPC->Transform1to2(x,index);
859 ////////////////////////////////////////////////////////////////////////
863 ////////////////////////////////////////////////////////////////////////
865 // Start of implementation of the class TPCCmpTr
867 ////////////////////////////////////////////////////////////////////////
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);
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;}
891 // tmp method, should go to TrackReferenceTPC
892 Float_t TR2LocalX(AliTrackReference *trackRef,
893 AliTPCParam *paramTPC);
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
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
906 char *fFnGenTracks; //! input file name with gen tracks
907 TFile *fFileGenTracks;
908 TTree *fTreeGenTracks;
910 char *fFnHits; //! input file name with gAlice object (needed for B)
911 TFile *fFileHits; //! input file with gAlice
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
917 AliTPCtrack *fTPCTrack; //! pointer to TPC track to connect branch
918 Int_t *fIndexRecTracks; //! index of particle label in the TreeT_TPC
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
931 Int_t fNParticles; //! number of particles in the input tree genTracks
932 Int_t fDebug; //! debug flag
934 Int_t fNextTreeGenEntryToRead; //! last entry already read from genTracks tree
935 TPCGenInfo *fGenInfo; //! container for all the details
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
943 ClassDef(TPCCmpTr,1) // class which creates and fills tree with TPCGenTrack objects
947 ////////////////////////////////////////////////////////////////////////
953 ////////////////////////////////////////////////////////////////////////
954 TPCCmpTr::TPCCmpTr(char* fnRecTracks,
958 Int_t nEvents, Int_t firstEvent)
961 fFnRecTracks = fnRecTracks;
962 fFnGenTracks = fnGenTracks;
965 fFirstEventNr = firstEvent;
966 fEventNr = firstEvent;
969 ////////////////////////////////////////////////////////////////////////
970 void TPCCmpTr::Reset()
976 fFnCmp = "cmpTracks.root";
977 fFnHits = "galice.root";
983 fRowsWithDigitsInn = 0;
985 fRowsTrackLength = 0;
991 fFnRecTracks = "tpc.tracks.root";
997 ////////////////////////////////////////////////////////////////////////
998 TPCCmpTr::~TPCCmpTr()
1002 ////////////////////////////////////////////////////////////////////////
1003 Int_t TPCCmpTr::Exec(Int_t nEvents, Int_t firstEventNr)
1006 fFirstEventNr = firstEventNr;
1010 ////////////////////////////////////////////////////////////////////////
1011 Int_t TPCCmpTr::Exec()
1016 fDigitRow = new digitRow();
1019 cerr<<"output tree not created"<<endl;
1022 fFileHits = OpenAliceFile(fFnHits,kTRUE,"read"); //gAlice is read here
1024 cerr<<"Cannot open file with gAlice object "<<fFnHits<<endl;
1030 fParamTPC = LoadTPCParam(fFileHits);
1032 cerr<<"TPC parameters not found and could not be created"<<endl;
1036 if (!ConnectGenTree()) {
1037 cerr<<"Cannot connect tree with generated tracks"<<endl;
1041 fNextTreeGenEntryToRead = 0;
1042 cerr<<"fFirstEventNr, fNEvents: "<<fFirstEventNr<<" "<<fNEvents<<endl;
1043 for (Int_t eventNr = fFirstEventNr; eventNr < fFirstEventNr+fNEvents;
1045 fNParticles = gAlice->GetEvent(fEventNr);
1046 fIndexRecTracks = new Int_t[fNParticles];
1047 for (Int_t i = 0; i<fNParticles; i++) {
1048 fIndexRecTracks[i] = -1;
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;
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;
1060 delete fTreeRecTracks;
1062 delete [] fIndexRecTracks;
1067 cerr<<"Exec finished"<<endl;
1077 ////////////////////////////////////////////////////////////////////////
1078 Bool_t TPCCmpTr::ConnectGenTree()
1081 // connect all branches from the genTracksTree
1082 // use the same variables as for the new cmp tree, it may work
1084 fFileGenTracks = TFile::Open(fFnGenTracks,"READ");
1085 if (!fFileGenTracks) {
1086 cerr<<"Error in ConnectGenTree: cannot open file "<<fFnGenTracks<<endl;
1089 fTreeGenTracks = (TTree*)fFileGenTracks->Get("genTracksTree");
1090 if (!fTreeGenTracks) {
1091 cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
1092 <<fFnGenTracks<<endl;
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);
1106 cout<<"Number of gen. tracks with TR: "<<fTreeGenTracks->GetEntries()<<endl;
1112 ////////////////////////////////////////////////////////////////////////
1113 void TPCCmpTr::CreateTreeCmp()
1115 fFileCmp = TFile::Open(fFnCmp,"RECREATE");
1117 cerr<<"Error in CreateTreeCmp: cannot open file "<<fFnCmp<<endl;
1120 fTreeCmp = new TTree("TPCcmpTracks","TPCcmpTracks");
1121 TBranch *branchBits = fTreeCmp->Branch("bitsBranch", "digitRow", &fDigitRow, 4000, 0);
1123 cerr<<"Error in CreateTreeCmp: cannot create branch."<<endl;
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");
1133 fTreeCmp->Branch("fReconstructed",&fReconstructed,"fReconstructed/I");
1134 fTreeCmp->Branch("fTPCTrack","AliTPCtrack",&fTPCTrack);
1136 fTreeCmp->Branch("Particle","TParticle",&fParticle);
1137 fTreeCmp->Branch("fVDist",&fVDist,"fVDist[4]/D");
1138 fTreeCmp->Branch("TR","AliTrackReference",&fTrackRef);
1140 fTreeCmp->AutoSave();
1142 ////////////////////////////////////////////////////////////////////////
1143 void TPCCmpTr::CloseOutputFile()
1146 cerr<<"File "<<fFnCmp<<" not found as an open file."<<endl;
1156 ////////////////////////////////////////////////////////////////////////
1158 Float_t TPCCmpTr::TR2LocalX(AliTrackReference *trackRef,
1159 AliTPCParam *paramTPC) {
1161 Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
1163 paramTPC->Transform0to1(x,index);
1164 paramTPC->Transform1to2(x,index);
1167 ////////////////////////////////////////////////////////////////////////
1169 Int_t TPCCmpTr::TreeTLoop(Int_t eventNr)
1172 // loop over all TPC reconstructed tracks and store info in memory
1176 if (!fFileRecTracks) fFileRecTracks = TFile::Open(fFnRecTracks,"read");
1177 if (!fFileRecTracks->IsOpen()) {
1178 cerr<<"Cannot open file "<<fFnRecTracks<<endl;
1182 char treeNameBase[11] = "TreeT_TPC_";
1184 sprintf(treeName,"%s%d",treeNameBase,eventNr);
1186 fTreeRecTracks=(TTree*)fFileRecTracks->Get(treeName);
1187 if (!fTreeRecTracks) {
1188 cerr<<"Can't get a tree with TPC rec. tracks named "<<treeName<<endl;
1192 Int_t nEntries = (Int_t) fTreeRecTracks->GetEntries();
1193 if (fDebug > 2) cout<<"Event, rec. tracks: "<<eventNr<<" "
1195 TBranch * br= fTreeRecTracks->GetBranch("tracks");
1196 br->SetAddress(&fTPCTrack);
1198 for (Int_t iEntry=0; iEntry<nEntries;iEntry++){
1199 br->GetEntry(iEntry);
1200 Int_t label = fTPCTrack->GetLabel();
1201 fIndexRecTracks[label] = iEntry;
1204 if (fDebug > 2) cerr<<"end of TreeTLoop"<<endl;
1208 ////////////////////////////////////////////////////////////////////////
1209 Int_t TPCCmpTr::TreeGenLoop(Int_t eventNr)
1212 // loop over all entries for a given event, find corresponding
1213 // rec. track and store in the fTreeCmp
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);
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;
1233 fRecPhi = fLambda = fRecPt_1 = 0.;
1236 cerr<<"fLabel, fIndexRecTracks[fLabel] "<<fLabel<<" "<<fIndexRecTracks[fLabel]<<endl;
1238 if (fIndexRecTracks[fLabel] >= 0) {
1239 Int_t nBytes = fTreeRecTracks->GetEvent(fIndexRecTracks[fLabel]);
1242 fdEdx = fTPCTrack->GetdEdx();
1243 Double_t Localx = TR2LocalX(fTrackRef,fParamTPC);
1245 cerr<<"Track local X before prolongation: "<<fTPCTrack->GetX()<<endl;
1247 fTPCTrack->PropagateTo(Localx);
1250 cerr<<"Track local X after prolongation: "<<fTPCTrack->GetX()<<endl;
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]);
1266 fTreeCmp->AutoSave();
1267 // delete gAlice; gAlice = 0;
1268 // fFileHits->Close();
1270 if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
1274 ////////////////////////////////////////////////////////////////////////