]>
Commit | Line | Data |
---|---|---|
b3d9d240 | 1 | #include <fstream.h> |
2 | ||
3 | #include <TClassTable.h> | |
4 | #include <TFile.h> | |
5 | #include <TTree.h> | |
6 | #include <TString.h> | |
7 | #include <TClonesArray.h> | |
8 | #include <TParticle.h> | |
9 | ||
10 | #include "AliRun.h" | |
11 | #include "AliITS.h" | |
12 | #include "AliITSgeom.h" | |
13 | #include "AliITSRecPoint.h" | |
14 | ||
15 | Int_t AliITSStoreFindableTracksCompiled | |
16 | (Int_t nMinClusters = 5, const Text_t *evname = "galice", Int_t evnum = 0) | |
17 | { | |
18 | // Make sure that ALICE objects are loaded | |
19 | if (gAlice) { | |
20 | delete gAlice; | |
21 | gAlice = 0; | |
22 | } | |
23 | ||
24 | // Define the names of all involved files | |
25 | TString strEventFile(evname); | |
26 | TString strOutputFile(evname); | |
27 | strEventFile.Append(".root"); | |
28 | strOutputFile.Append("_tracks_"); | |
29 | strOutputFile += evnum; | |
30 | strOutputFile.Append(".root"); | |
31 | ||
32 | // Connect the Root Galice file containing Geometry, Kine and Hits | |
33 | TFile *fileEvent = (TFile*)gROOT->GetListOfFiles()->FindObject(strEventFile); | |
34 | if (!fileEvent) fileEvent = new TFile(strEventFile,"UPDATE"); | |
35 | ||
36 | // Get AliRun object from file | |
37 | gAlice = (AliRun*)fileEvent->Get("gAlice"); | |
38 | if (gAlice) cout << "OK, found an AliRun object in file" << endl; | |
39 | ||
40 | // Get ITS related objects and data | |
41 | AliITS* ITS =(AliITS *)gAlice->GetDetector("ITS"); | |
42 | if (!ITS) { | |
43 | cerr << "ITS object not found!" << endl; | |
44 | return 1; | |
45 | } | |
46 | AliITSgeom *geometry = ITS->GetITSgeom(); | |
47 | if (!geometry) { | |
48 | cerr << "ITS geometry object not found!" << endl; | |
49 | return 2; | |
50 | } | |
51 | ||
52 | // Count the number of modules per layer | |
53 | Int_t nLadders, nDetectors, mod_min[6], mod_max[6]; | |
54 | for(Int_t i = 0; i < 6; i++) { | |
55 | nLadders = geometry->GetNladders(i + 1); | |
56 | nDetectors = geometry->GetNdetectors(i + 1); | |
57 | mod_min[i] = geometry->GetModuleIndex(i + 1, 1, 1); | |
58 | mod_max[i] = geometry->GetModuleIndex(i + 1, nLadders, nDetectors); | |
59 | } | |
60 | ||
61 | // Load event and ITS recpoints | |
62 | Int_t nParticles = gAlice->GetEvent(evnum); | |
63 | cout << "Event number: " << evnum << endl; | |
64 | cout << "# particles : " << nParticles <<endl; | |
65 | if (nParticles <= 0) { | |
66 | cerr << "Can't have <= 0 particles!" << endl; | |
67 | return 3; | |
68 | } | |
69 | AliITSRecPoint *recp = 0; | |
70 | TClonesArray *recPoints = ITS->RecPoints(); | |
71 | TObjArray *particles = gAlice->Particles(); | |
72 | Int_t nTracks = gAlice->GetNtrack(); //FCA correction | |
73 | Bool_t *hitITSLayer[6]; | |
74 | for (Int_t i = 0; i < 6; i++) { | |
75 | hitITSLayer[i] = new Bool_t[nTracks]; | |
76 | for (Int_t j = 0; j < nTracks; j++) hitITSLayer[i][j] = kFALSE; | |
77 | } | |
78 | ||
79 | // Load recpoints in event | |
80 | TTree *TR = gAlice->TreeR(); | |
81 | if (!TR) { | |
82 | cerr << "TreeR object not found!" << endl; | |
83 | return 4; | |
84 | } | |
85 | ||
86 | // Scan recpoints and define findable tracks | |
87 | Int_t nModules = (Int_t)TR->GetEntries(), nPoints = 0, nEmpty = 0; | |
88 | cout << "Found " << nModules; | |
89 | cout << " entries in the TreeR (must be one per module!)" << endl; | |
90 | for (Int_t layer = 1; layer <= 6; layer++) { | |
91 | for (Int_t mod = mod_min[layer - 1]; mod <= mod_max[layer - 1]; mod++) { | |
92 | ITS->ResetRecPoints(); | |
93 | TR->GetEntry(mod); | |
94 | nPoints = recPoints->GetEntries(); | |
95 | if(!nPoints) { | |
96 | nEmpty++; | |
97 | continue; | |
98 | } | |
99 | for (Int_t point = 0; point < nPoints; point++) { | |
100 | recp = (AliITSRecPoint*)recPoints->UncheckedAt(point); | |
101 | for (Int_t it = 0; it < 3; it++) { | |
102 | Int_t track = recp->GetLabel(it); | |
103 | if(track < 0) continue; | |
104 | if(track > nTracks) { | |
105 | cout << "Found track index " << track; | |
106 | cout << " whilw gAlice->GetNtrack() = " << nTracks << endl; | |
107 | continue; | |
108 | } | |
109 | hitITSLayer[layer - 1][track] = kTRUE; | |
110 | } // loop over recpoint labels | |
111 | } //loop over points | |
112 | } //loop over modules | |
113 | } //loop over layers | |
114 | cout << "Found " << nEmpty << " empty modules" << endl; | |
115 | ||
116 | // Scan the file of tracks in TPC to retrieve the findable TPC tracks | |
117 | TString strLabelsTPC; | |
118 | Int_t label, pdg_code, nFindablesTPC = 0; | |
119 | Double_t dummy; | |
120 | ifstream tpc("good_tracks_tpc"); | |
121 | while (tpc >> label >> pdg_code) { | |
122 | for (Int_t i = 0; i < 6; i++) tpc >> dummy; | |
123 | nFindablesTPC++; | |
124 | strLabelsTPC.Append(Form("[%d]", label)); | |
125 | } | |
126 | ||
127 | // Define the TTree with tracks data by means of a set of variables | |
128 | Int_t nFindablesITS = 0, nFindablesITSTPC = 0; | |
129 | Int_t nhits, tpc_ok, mother, entry = 0; | |
130 | Double_t vx, vy, vz; | |
131 | Double_t px, py, pz, pt; | |
132 | ||
133 | TTree *tree = new TTree("Tracks", "Findable tracks in ITS"); | |
134 | ||
135 | tree->Branch("vx", &vx, "vx/D"); | |
136 | tree->Branch("vy", &vy, "vy/D"); | |
137 | tree->Branch("vz", &vz, "vz/D"); | |
138 | tree->Branch("px", &px, "px/D"); | |
139 | tree->Branch("py", &py, "py/D"); | |
140 | tree->Branch("pz", &pz, "pz/D"); | |
141 | tree->Branch("pt", &pt, "pt/D"); | |
142 | tree->Branch("label", &label, "label/I"); | |
143 | tree->Branch("entry", &entry, "entry/I"); | |
144 | tree->Branch("mother", &mother, "mother/I"); | |
145 | tree->Branch("pdg_code", &pdg_code, "pdg_code/I"); | |
146 | tree->Branch("nhits", &nhits, "nhits/I"); | |
147 | tree->Branch("tpc_ok", &tpc_ok, "tpc_ok/I"); | |
148 | ||
149 | // Fill the tree | |
150 | cout << endl; | |
151 | TParticle *p = 0; | |
152 | for (Int_t i = 0; i < nTracks; i++) { | |
153 | nhits = 0; | |
154 | for (Int_t j = 0; j < 6; j++) if (hitITSLayer[j][i]) nhits++; | |
155 | if (nhits < nMinClusters) continue; | |
156 | p = gAlice->Particle(i); | |
157 | px = p->Px(); | |
158 | py = p->Py(); | |
159 | pz = p->Pz(); | |
160 | pt = p->Pt(); | |
161 | vx = p->Vx(); | |
162 | vy = p->Vy(); | |
163 | vz = p->Vz(); | |
164 | mother = p->GetFirstMother(); | |
165 | cout << "Track " << i << " stored\r" << flush; | |
166 | tpc_ok = (strLabelsTPC.Contains(Form("[%d]", i))); | |
167 | pdg_code = p->GetPdgCode(); | |
168 | label = i; | |
169 | nFindablesITS++; | |
170 | if (tpc_ok) nFindablesITSTPC++; | |
171 | tree->Fill(); | |
172 | entry++; | |
173 | } | |
174 | ||
175 | // Save into a file | |
176 | TFile *fileOutput = new TFile(strOutputFile, "recreate"); | |
177 | tree->Write(); | |
178 | fileOutput->Close(); | |
179 | ||
180 | cout << "# findable tracks in TPC : " << nFindablesTPC << endl; | |
181 | cout << "# findable tracks in ITS : " << nFindablesITS << endl; | |
182 | cout << "# findable tracks in ITS+TPC : " << nFindablesITSTPC << endl; | |
183 | ||
184 | return 0; | |
185 | } |