]>
Commit | Line | Data |
---|---|---|
10a5e8ec | 1 | #ifndef __CINT__ |
2 | #include <iostream.h> | |
3 | #include "AliTRDtracker.h" | |
4 | ||
5 | #include "AliTRDcluster.h" | |
6 | #include "AliTRDhit.h" | |
7 | #include "AliTRDv1.h" | |
8 | #include "AliTRDgeometry.h" | |
9 | #include "AliTRDparameter.h" | |
10 | #include "alles.h" | |
11 | #include "AliTRDmcTrack.h" | |
12 | ||
13 | #include "AliTRDtrack.h" | |
14 | #include "AliTrackReference.h" | |
15 | ||
16 | #include "TFile.h" | |
17 | #include "TParticle.h" | |
18 | #include "TStopwatch.h" | |
19 | #endif | |
20 | ||
21 | ||
22 | void AliTRDsaveTrackableSeeds() { | |
23 | ||
24 | TObjArray mctracks(2000); | |
25 | TObjArray *TRDmcTracks = &mctracks; | |
26 | ||
27 | TFile *geofile =TFile::Open("AliTRDclusters.root"); | |
28 | AliTRDtracker *Tracker = new AliTRDtracker(geofile); | |
29 | Int_t nEvent = 0; | |
30 | Tracker->SetEventNumber(nEvent); | |
31 | ||
32 | AliTRDgeometry *fGeo = (AliTRDgeometry*) geofile->Get("TRDgeometry"); | |
33 | //AliTRDparameter *fPar = (AliTRDparameter*) geofile->Get("TRDparameter"); | |
34 | ||
35 | Char_t *alifile = "galice.root"; | |
36 | ||
37 | // Connect the AliRoot file containing Geometry, Kine, Hits, and Digits | |
38 | TFile *gafl = (TFile*) gROOT->GetListOfFiles()->FindObject(alifile); | |
39 | if (!gafl) { | |
40 | cout << "Open the ALIROOT-file " << alifile << endl; | |
41 | gafl = new TFile(alifile); | |
42 | } | |
43 | else { | |
44 | cout << alifile << " is already open" << endl; | |
45 | } | |
46 | ||
47 | // Get AliRun object from file or create it if not on file | |
48 | gAlice = (AliRun*) gafl->Get("gAlice"); | |
49 | if (gAlice) | |
50 | cout << "AliRun object found on file" << endl; | |
51 | else | |
52 | gAlice = new AliRun("gAlice","Alice test program"); | |
53 | ||
54 | AliTRDv1 *fTRD = (AliTRDv1*) gAlice->GetDetector("TRD"); | |
55 | ||
56 | // Import the Trees for the event nEvent in the file | |
57 | const Int_t nparticles = gAlice->GetEvent(nEvent); | |
58 | ||
59 | printf("found %d particles in event %d \n", nparticles, nEvent); | |
60 | ||
61 | // Create TRDmcTracks for each tpc seed | |
62 | Int_t label = -1, charge = 0; | |
63 | Bool_t primary; | |
64 | Float_t mass; | |
65 | Int_t pdg_code; | |
66 | ||
67 | TDirectory *savedir=gDirectory; | |
68 | ||
69 | TFile *in=TFile::Open("AliTPCBackTracks.root"); | |
70 | if (!in->IsOpen()) { | |
71 | cerr<<"can't open file AliTPCBackTracks.root !\n"; return; | |
72 | } | |
73 | ||
74 | char tname[100]; | |
75 | sprintf(tname,"seedsTPCtoTRD_%d",nEvent); | |
76 | TTree *seedTree=(TTree*)in->Get(tname); | |
77 | if (!seedTree) { | |
78 | cerr<<"AliTRDtracker::PropagateBack(): "; | |
79 | cerr<<"can't get a tree with seeds from TPC !\n"; | |
80 | } | |
81 | ||
82 | AliTPCtrack *seed=new AliTPCtrack; | |
83 | seedTree->SetBranchAddress("tracks",&seed); | |
84 | ||
85 | Int_t nSeeds = (Int_t) seedTree->GetEntries(); | |
86 | for (Int_t is=0; is<nSeeds; is++) { | |
87 | seedTree->GetEvent(is); | |
88 | Int_t lbl = seed->GetLabel(); | |
89 | if(TMath::Abs(lbl) > nparticles) { | |
90 | printf("extra high seed label %d \n", lbl); | |
91 | continue; | |
92 | } | |
93 | TParticle *p = gAlice->Particle(TMath::Abs(lbl)); | |
94 | ||
95 | primary = kTRUE; | |
96 | if (p->GetFirstMother()>=0) primary = kFALSE; | |
97 | ||
98 | pdg_code = (Int_t) p->GetPdgCode(); | |
99 | ||
100 | if ((pdg_code == 10010020) || | |
101 | (pdg_code == 10010030) || | |
102 | (pdg_code == 50000050) || | |
103 | (pdg_code == 50000051) || | |
104 | (pdg_code == 10020040)) { | |
105 | mass = 0.; | |
106 | charge = 0; | |
107 | } | |
108 | else { | |
109 | TParticlePDG *pdg = p->GetPDG(); | |
110 | charge = (Int_t) pdg->Charge(); | |
111 | mass=pdg->Mass(); | |
112 | } | |
113 | ||
114 | AliTRDmcTrack *t = new AliTRDmcTrack(TMath::Abs(lbl), lbl, primary, | |
115 | mass, charge, pdg_code); | |
116 | TRDmcTracks->AddLast(t); | |
117 | } | |
118 | delete seed; | |
119 | delete seedTree; | |
120 | ||
121 | savedir->cd(); | |
122 | ||
123 | Int_t i, mctIndex[nparticles]; | |
124 | for(i=0; i<nparticles; i++) mctIndex[i]=-1; | |
125 | ||
126 | Int_t nMCtracks = TRDmcTracks->GetEntriesFast(); | |
127 | AliTRDmcTrack *mct = 0; | |
128 | for(i = 0; i < nMCtracks; i++) { | |
129 | mct = (AliTRDmcTrack*) TRDmcTracks->UncheckedAt(i); | |
130 | mctIndex[mct->GetTrackIndex()] = i; | |
131 | } | |
132 | ||
133 | // Loop through the TRD rec points and set XYZ and Pin and Pout | |
134 | ||
135 | Double_t xin[6], xout[6]; | |
136 | Double_t Px, Py, Pz; | |
137 | ||
138 | Float_t pos[3], rot[3]; | |
139 | ||
140 | for(Int_t j = 0; j < 6; j++) { | |
141 | xin[j] = Tracker->GetX(0,j,0)-3; | |
142 | xout[j] = Tracker->GetX(0,j,0); | |
143 | } | |
144 | ||
145 | TTree *trRefTree = gAlice->TreeTR(); | |
146 | if (!trRefTree) { | |
147 | cout << "<AliTRDreadTrackRef> No TR tree found" << endl; | |
148 | return; | |
149 | } | |
150 | ||
151 | Int_t nBytes = 0; | |
152 | Bool_t fEntrance, fExit; | |
153 | ||
154 | Int_t nTrack = (Int_t) trRefTree->GetEntries(); | |
155 | cout << "<AliTRDreadTrackRef> Found " << nTrack | |
156 | << " primary particles with track refs" << endl; | |
157 | ||
158 | for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) { | |
159 | ||
160 | printf("TrackRefs for track %d out of %d \n",iTrack,nTrack); | |
161 | gAlice->ResetTrackReferences(); | |
162 | nBytes += trRefTree->GetEvent(iTrack); | |
163 | ||
164 | AliTrackReference *tr = 0; | |
165 | tr = (AliTrackReference*) fTRD->FirstTrackReference(-1); | |
166 | ||
167 | while (tr) { | |
168 | ||
169 | Int_t track = tr->GetTrack(); | |
170 | if(mctIndex[track] >= 0) { | |
171 | mct = (AliTRDmcTrack*) TRDmcTracks->UncheckedAt(mctIndex[track]); | |
172 | ||
173 | pos[0] = tr->X(); | |
174 | pos[1] = tr->Y(); | |
175 | pos[2] = tr->Z(); | |
176 | ||
177 | Double_t phi = TMath::ATan2(pos[1],pos[0]); | |
178 | if(phi < 0) phi = 2 * TMath::Pi() + phi; | |
179 | ||
180 | Int_t sec = ((Int_t) (phi*180/TMath::Pi())) / 20; | |
181 | fGeo->Rotate(fGeo->GetDetector(0,0,17-sec), pos, rot); | |
182 | ||
183 | fEntrance = kFALSE; fExit = kFALSE; | |
184 | for(i = 0; i < 6; i++) { | |
185 | if(TMath::Abs(xin[i] - rot[0]) < 1.4) { fEntrance = kTRUE; break; } | |
186 | if(TMath::Abs(xout[i] - rot[0]) < 1.4) { fExit = kTRUE; break; } | |
187 | } | |
188 | ||
189 | Px = tr->Px(); | |
190 | Py = tr->Py(); | |
191 | Pz = tr->Pz(); | |
192 | ||
193 | if(fEntrance) { | |
194 | mct->SetXYZin(i, rot[0], rot[1], rot[2]); | |
195 | mct->SetPin(i, Px, Py, Pz); | |
196 | /* | |
197 | printf("entr plane %d: rp_x = %f vs %f = xin \n",i,rot[0],xin[i]); | |
198 | printf(" : Px, Py, Pz = %f, %f, %f \n",Px,Py,Pz); | |
199 | printf(" : y, z = %f, %f\n",rot[1],rot[2]); | |
200 | */ | |
201 | } | |
202 | if(fExit) { | |
203 | mct->SetXYZout(i, rot[0], rot[1], rot[2]); | |
204 | mct->SetPout(i, Px, Py, Pz); | |
205 | /* | |
206 | printf("exit plane %d: rp_x = %f vs %f = xout \n",i,rot[0],xout[i]); | |
207 | printf(" : Px, Py, Pz = %f, %f, %f \n",Px,Py,Pz); | |
208 | printf(" : y, z = %f, %f\n",rot[1],rot[2]); | |
209 | */ | |
210 | } | |
211 | } | |
212 | tr = (AliTrackReference *) fTRD->NextTrackReference(); | |
213 | } | |
214 | } | |
215 | ||
216 | // Loop through TRD clusters and assign indexes to MC tracks | |
217 | ||
218 | Char_t *clusterfile = "AliTRDclusters.root"; | |
219 | TObjArray carray(2000); | |
220 | TObjArray *ClustersArray = &carray; | |
221 | Tracker->ReadClusters(ClustersArray,clusterfile); | |
222 | ||
223 | printf("done with ReadClusters()\n"); | |
224 | ||
225 | Int_t nClusters = carray.GetEntriesFast(); | |
226 | ||
227 | Int_t track, det, det0, det1, ltb, plane, ind0, ind1, ncl; | |
228 | Double_t y, y0, y1, q, q0, q1; | |
229 | ||
230 | AliTRDcluster *c0 = NULL, *c1 = NULL; | |
231 | printf("nClusters = %d \n", nClusters); | |
232 | ||
233 | for (Int_t i = 0; i < nClusters; i++) { | |
234 | ||
235 | printf("\r assigning cluster %d out of %d", i, nClusters); | |
236 | AliTRDcluster *cl = (AliTRDcluster *) ClustersArray->UncheckedAt(i); | |
237 | ||
238 | for(Int_t j=0; j<3; j++) { | |
239 | track = cl->GetLabel(j); | |
240 | if(track < 0) continue; | |
241 | if(track >= nparticles) { | |
242 | printf("Track index %d is larger than total number of particles %d\n" | |
243 | ,track,nparticles); | |
244 | continue; | |
245 | } | |
246 | if(mctIndex[track] < 0) continue; | |
247 | mct = (AliTRDmcTrack*) TRDmcTracks->UncheckedAt(mctIndex[track]); | |
248 | ||
249 | label = mct->GetTrackIndex(); | |
250 | if(label != track) { | |
251 | printf("Got label %d while expected %d !\n", label, track); | |
252 | continue; | |
253 | } | |
254 | ||
255 | ncl = mct->GetNumberOfClusters(); | |
256 | mct->SetNumberOfClusters(ncl+1); | |
257 | ||
258 | det=cl->GetDetector(); | |
259 | plane = fGeo->GetPlane(det); | |
260 | ltb=cl->GetLocalTimeBin(); | |
261 | if((ltb < 0) || (ltb >= kMAX_TB)) continue; | |
262 | ||
263 | ind0 = mct->GetClusterIndex(ltb, plane, 0); | |
264 | ind1 = mct->GetClusterIndex(ltb, plane, 1); | |
265 | if(ind0 < 0) { | |
266 | mct->Update(ltb,plane,0,i); | |
267 | } else { | |
268 | c0 = (AliTRDcluster *) ClustersArray->UncheckedAt(ind0); | |
269 | det0 = c0->GetDetector(); | |
270 | y = cl->GetY(); | |
271 | y0 = c0->GetY(); | |
272 | q = TMath::Abs(cl->GetQ()); | |
273 | q0 = TMath::Abs(c0->GetQ()); | |
274 | if((det == det0) && (TMath::Abs(y-y0) < 1.5)) { | |
275 | if(q > q0) mct->Update(ltb,plane,0,i); | |
276 | } else { | |
277 | if (ind1 < 0) { | |
278 | mct->Update(ltb,plane,1,i); | |
279 | } | |
280 | else { | |
281 | c1 = (AliTRDcluster *) ClustersArray->UncheckedAt(ind1); | |
282 | det1 = c1->GetDetector(); | |
283 | y1 = c1->GetY(); | |
284 | q1 = TMath::Abs(c1->GetQ()); | |
285 | if((det == det0) && (TMath::Abs(y-y0) < 1.5)) { | |
286 | if(q > q1) mct->Update(ltb,plane,1,i); | |
287 | } | |
288 | } | |
289 | } | |
290 | } | |
291 | } | |
292 | } | |
293 | ||
294 | // Write TRDmcTracks in output file | |
295 | ||
296 | savedir=gDirectory; | |
297 | Char_t *filename = "AliTRDtrackableSeeds.root"; | |
298 | TFile *out = new TFile(filename,"RECREATE"); | |
299 | ||
300 | TTree *tracktree = new TTree("MCtracks","TRD MC tracks"); | |
301 | ||
302 | AliTRDmcTrack *iotrack=0; | |
303 | tracktree->Branch("MCtracks","AliTRDmcTrack",&iotrack,32000,0); | |
304 | ||
305 | Int_t ntracks = TRDmcTracks->GetEntriesFast(); | |
306 | ||
307 | for (Int_t i=0; i<ntracks; i++) { | |
308 | AliTRDmcTrack *pt=(AliTRDmcTrack*)TRDmcTracks->UncheckedAt(i); | |
309 | iotrack=pt; | |
310 | tracktree->Fill(); | |
311 | printf("Put track with label %d and %d clusters in the output tree \n", | |
312 | pt->GetTrackIndex(),pt->GetNumberOfClusters()); | |
313 | } | |
314 | ||
315 | tracktree->Write(); | |
316 | out->Close(); | |
317 | savedir->cd(); | |
318 | ||
319 | return; | |
320 | } | |
321 |