]>
Commit | Line | Data |
---|---|---|
88ce87da | 1 | /**************************************************************************** |
2 | * This macro is supposed to do reconstruction in the ITS via Kalman * | |
3 | * tracker V2. The ITStracker is feeded with parametrized TPC tracks * | |
4 | * * | |
5 | * It does the following steps: * | |
6 | * 1) TPC tracking parameterization * | |
7 | * 2) ITS cluster finding V2 (via fast points !) * | |
8 | * 3) ITS track finding V2 * | |
9 | * 4) Create a reference file with simulation info (p,PDG...) * | |
10 | * * | |
11 | * (Origin: A.Dainese, Padova, andrea.dainese@pd,infn.it * | |
12 | * from AliBarrelReconstruction.C I.Belikov, CERN, Jouri.Belikov@cern.ch) * | |
13 | ****************************************************************************/ | |
14 | ||
15 | #ifndef __CINT__ | |
16 | #include "alles.h" | |
17 | #include "AliMagF.h" | |
18 | #include "AliITS.h" | |
19 | #include "AliITSgeom.h" | |
20 | #include "AliITSRecPoint.h" | |
21 | #include "AliITSclusterV2.h" | |
22 | #include "AliITSsimulationFastPoints.h" | |
23 | #include "AliITStrackerV2.h" | |
24 | #include "AliTPCtrackerParam.h" | |
25 | #endif | |
26 | ||
27 | typedef struct { | |
28 | Int_t lab; | |
29 | Int_t pdg; | |
30 | Int_t mumpdg; | |
31 | Float_t Vx,Vy,Vz; | |
32 | Float_t Px,Py,Pz; | |
33 | } RECTRACK; | |
34 | ||
35 | ||
36 | Int_t TPCParamTracks(const Char_t *galice,const Char_t *outname,const Int_t coll,const Double_t Bfield,Int_t n); | |
37 | Int_t ITSFindClusters(const Char_t *inname,const Char_t *outname,Int_t n); | |
38 | Int_t ITSFindTracks(const Char_t *galice,const Char_t *inname,const Char_t *inname2,const Char_t *outname,Int_t n); | |
39 | Int_t ITSMakeRefFile(const Char_t *galice, const Char_t *inname, const Char_t *outname, Int_t n); | |
40 | ||
41 | Int_t AliBarrelRec_TPCparam(Int_t n=1) { | |
42 | const Char_t *TPCtrkNameS="AliTPCtracksParam.root"; | |
43 | const Char_t *galiceName="galice.root"; | |
44 | const Char_t *ITSclsName="AliITSclustersV2.root"; | |
45 | const Char_t *ITStrkName="AliITStracksV2.root"; | |
46 | const Char_t *ITSrefName="ITStracksRefFile.root"; | |
47 | ||
48 | // set here the code for the type of collision (needed for TPC tracking | |
49 | // parameterization). available collisions: | |
50 | // | |
51 | // coll = 0 -> PbPb6000 (HIJING with b<2fm) | |
52 | const Int_t collcode = 0; | |
53 | // set here the value of the magnetic field | |
54 | const Double_t BfieldValue = 0.4; | |
55 | ||
56 | ||
57 | ||
58 | AliKalmanTrack::SetConvConst(100/0.299792458/BfieldValue); | |
59 | ||
60 | ||
61 | // ********** Build TPC tracks with parameterization *********** // | |
62 | if (TPCParamTracks(galiceName,TPCtrkNameS,collcode,BfieldValue,n)) { | |
63 | cerr<<"Failed to get TPC hits !\n"; | |
64 | return 1; | |
65 | } | |
66 | ||
67 | ||
68 | // ********** Find ITS clusters *********** // | |
69 | if (ITSFindClusters(galiceName,ITSclsName,n)) { | |
70 | cerr<<"Failed to get ITS clusters !\n"; | |
71 | return 1; | |
72 | } | |
73 | ||
74 | ||
75 | // ********* Find ITS tracks *********** // | |
76 | if (ITSFindTracks(galiceName,TPCtrkNameS,ITSclsName,ITStrkName,n)) { | |
77 | cerr<<"Failed to get ITS tracks !\n"; | |
78 | return 1; | |
79 | } | |
80 | ||
81 | ||
82 | // ********* Make ITS tracks reference file *********** // | |
83 | if (ITSMakeRefFile(galiceName,ITStrkName,ITSrefName,n)) { | |
84 | cerr<<"Failed to get ITS tracks ref file!\n"; | |
85 | return 1; | |
86 | } | |
87 | ||
88 | ||
89 | return 0; | |
90 | } | |
91 | ||
92 | ||
93 | Int_t TPCParamTracks(const Char_t *galice, const Char_t *outname, | |
94 | const Int_t coll,const Double_t Bfield,Int_t n) { | |
95 | ||
96 | cerr<<"\n*******************************************************************\n"; | |
97 | Int_t rc; | |
98 | ||
99 | const Char_t *name="TPCParamTracks"; | |
100 | cerr<<'\n'<<name<<"...\n"; | |
101 | gBenchmark->Start(name); | |
102 | ||
103 | TFile *outfile=TFile::Open(outname,"recreate"); | |
104 | TFile *infile =TFile::Open(galice); | |
105 | ||
106 | AliTPCtrackerParam tracker(coll,Bfield); | |
107 | rc = tracker.BuildTPCtracks(infile,outfile,n); | |
108 | ||
109 | delete gAlice; gAlice=0; | |
110 | ||
111 | infile->Close(); | |
112 | outfile->Close(); | |
113 | ||
114 | gBenchmark->Stop(name); | |
115 | gBenchmark->Show(name); | |
116 | ||
117 | return rc; | |
118 | } | |
119 | ||
120 | Int_t ITSFindClusters(const Char_t *inname, const Char_t *outname, Int_t n) { | |
121 | ||
122 | ||
123 | cerr<<"\n*******************************************************************\n"; | |
124 | ||
125 | Int_t rc=0; | |
126 | const Char_t *name="ITSFindClusters"; | |
127 | cerr<<'\n'<<name<<"...\n"; | |
128 | gBenchmark->Start(name); | |
129 | TFile *out=TFile::Open(outname,"recreate"); | |
130 | TFile *in =TFile::Open(inname,"update"); | |
131 | ||
132 | ||
133 | if (!(gAlice=(AliRun*)in->Get("gAlice"))) { | |
134 | cerr<<"Can't get gAlice !\n"; | |
135 | return 1; | |
136 | } | |
137 | ||
138 | ||
139 | AliITS *ITS = (AliITS*)gAlice->GetModule("ITS"); | |
140 | if (!ITS) { cerr<<"Can't get the ITS !\n"; return 1;} | |
141 | AliITSgeom *geom=ITS->GetITSgeom(); | |
142 | out->cd(); | |
143 | geom->Write(); | |
144 | ||
145 | Int_t ev=0; | |
146 | for (ev = 0; ev<n; ev++){ | |
147 | in->cd(); // !!!! MI directory must point to galice. - othervise problem with Tree -connection | |
148 | gAlice->GetEvent(ev); | |
149 | //gAlice->TreeR()->Reset(); //reset reconstructed tree | |
150 | ||
151 | ||
152 | TTree *pTree=gAlice->TreeR(); | |
153 | if (!pTree){ | |
154 | gAlice->MakeTree("R"); | |
155 | pTree = gAlice->TreeR(); | |
156 | } | |
157 | TBranch *branch=pTree->GetBranch("ITSRecPoints"); | |
158 | if (!branch) { | |
159 | //if not reconstructed ITS branch do reconstruction | |
160 | ITS->MakeBranch("R",0); | |
161 | //////////////// Taken from ITSHitsToFastPoints.C /////////////////////// | |
162 | for (Int_t i=0;i<3;i++) { | |
163 | ITS->SetSimulationModel(i,new AliITSsimulationFastPoints()); | |
164 | } | |
165 | Int_t nsignal=25; | |
166 | Int_t size=-1; | |
167 | Int_t bgr_ev=Int_t(ev/nsignal); | |
168 | ITS->HitsToFastRecPoints(ev,bgr_ev,size," ","All"," "); | |
169 | /////////////////////////////////////////////////////////////////////// | |
170 | gAlice->GetEvent(ev); //MI comment - in HitsToFast... they reset treeR to 0 | |
171 | //they overwrite full reconstructed event ???? ... so lets connect TreeR one more | |
172 | //time | |
173 | } | |
174 | ||
175 | ||
176 | ||
177 | out->cd(); | |
178 | TClonesArray *clusters=new TClonesArray("AliITSclusterV2",10000); | |
179 | char cname[100]; | |
180 | sprintf(cname,"TreeC_ITS_%d",ev); | |
181 | ||
182 | TTree *cTree=new TTree(cname,"ITS clusters"); | |
183 | cTree->Branch("Clusters",&clusters); | |
184 | ||
185 | pTree=gAlice->TreeR(); | |
186 | if (!pTree) { cerr<<"Can't get TreeR !\n"; return 1; } | |
187 | branch=pTree->GetBranch("ITSRecPoints"); | |
188 | if (!branch) { cerr<<"Can't get ITSRecPoints branch !\n"; return 1;} | |
189 | TClonesArray *points=new TClonesArray("AliITSRecPoint",10000); | |
190 | branch->SetAddress(&points); | |
191 | ||
192 | TClonesArray &cl=*clusters; | |
193 | Int_t nclusters=0; | |
194 | Int_t nentr=(Int_t)pTree->GetEntries(); | |
195 | AliITSgeom *geom=ITS->GetITSgeom(); | |
196 | ||
197 | for (Int_t i=0; i<nentr; i++) { | |
198 | if (!pTree->GetEvent(i)) {cTree->Fill(); continue;} | |
199 | Int_t lay,lad,det; geom->GetModuleId(i,lay,lad,det); | |
200 | Float_t x,y,zshift; geom->GetTrans(lay,lad,det,x,y,zshift); | |
201 | Double_t rot[9]; geom->GetRotMatrix(lay,lad,det,rot); | |
202 | Double_t yshift = x*rot[0] + y*rot[1]; | |
203 | Int_t ndet=(lad-1)*geom->GetNdetectors(lay) + (det-1); | |
204 | Int_t ncl=points->GetEntriesFast(); | |
205 | nclusters+=ncl; | |
206 | Float_t lp[5]; | |
207 | Int_t lab[6]; | |
208 | for (Int_t j=0; j<ncl; j++) { | |
209 | AliITSRecPoint *p=(AliITSRecPoint*)points->UncheckedAt(j); | |
210 | lp[0]=-p->GetX()-yshift; if (lay==1) lp[0]=-lp[0]; | |
211 | lp[1]=p->GetZ()+zshift; | |
212 | lp[2]=p->GetSigmaX2(); | |
213 | lp[3]=p->GetSigmaZ2(); | |
214 | lp[4]=p->GetQ(); | |
215 | lab[0]=p->GetLabel(0);lab[1]=p->GetLabel(1);lab[2]=p->GetLabel(2); | |
216 | lab[3]=ndet; | |
217 | ||
218 | Int_t label=lab[0]; | |
219 | TParticle *part=(TParticle*)gAlice->Particle(label); | |
220 | label=-3; | |
221 | while (part->P() < 0.005) { | |
222 | Int_t m=part->GetFirstMother(); | |
223 | if (m<0) {cerr<<"Primary momentum: "<<part->P()<<endl; break;} | |
224 | label=m; | |
225 | part=(TParticle*)gAlice->Particle(label); | |
226 | } | |
227 | if (lab[1]<0) lab[1]=label; | |
228 | else if (lab[2]<0) lab[2]=label; | |
229 | else cerr<<"No empty labels !\n"; | |
230 | ||
231 | new(cl[j]) AliITSclusterV2(lab,lp); | |
232 | } | |
233 | cTree->Fill(); clusters->Delete(); | |
234 | points->Delete(); | |
235 | } | |
236 | cTree->Write(); | |
237 | cerr<<"Number of clusters: "<<nclusters<<endl; | |
238 | delete cTree; delete clusters; delete points; | |
239 | ||
240 | } | |
241 | ||
242 | ||
243 | delete gAlice; gAlice=0; | |
244 | in->Close(); | |
245 | out->Close(); | |
246 | gBenchmark->Stop(name); | |
247 | gBenchmark->Show(name); | |
248 | ||
249 | return rc; | |
250 | } | |
251 | ||
252 | Int_t ITSFindTracks(const Char_t *galice, const Char_t * inname, | |
253 | const Char_t *inname2, const Char_t *outname, | |
254 | Int_t n) { | |
255 | ||
256 | ||
257 | cerr<<"\n*******************************************************************\n"; | |
258 | ||
259 | Int_t rc=0; | |
260 | const Char_t *name="ITSFindTracks"; | |
261 | cerr<<'\n'<<name<<"...\n"; | |
262 | gBenchmark->Start(name); | |
263 | ||
264 | ||
265 | TFile *out=TFile::Open(outname,"recreate"); | |
266 | TFile *in =TFile::Open(inname); | |
267 | TFile *in2 =TFile::Open(inname2); | |
268 | ||
269 | AliITSgeom *geom=(AliITSgeom*)gFile->Get("AliITSgeom"); | |
270 | if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;} | |
271 | ||
272 | ||
273 | for (Int_t ev=0; ev<n; ev++){ | |
274 | AliITStrackerV2 tracker(geom,ev); | |
275 | ||
276 | TArrayF o; | |
277 | o.Set(3); | |
278 | TFile* vtxFile = new TFile(galice); | |
279 | vtxFile->cd(); | |
280 | AliHeader* header = 0; | |
281 | TTree* treeE = (TTree*)gDirectory->Get("TE"); | |
282 | treeE->SetBranchAddress("Header",&header); | |
283 | treeE->GetEntry(ev); | |
284 | AliGenEventHeader* genHeader = header->GenEventHeader(); | |
285 | if(genHeader) { | |
286 | // get primary vertex position | |
287 | genHeader->PrimaryVertex(o); | |
288 | vtxFile->Close(); | |
289 | delete header; | |
290 | // set position of primary vertex | |
291 | Double_t vtx[3]; | |
292 | vtx[0]=o[0]; vtx[1]=o[1]; vtx[2]=o[2]; | |
293 | cerr<<"+++\n+++ Reading primary vertex position from galice.root\n+++\n+++ Setting primary vertex z = "<<vtx[2]<<" cm for ITS tracking\n+++\n"; | |
294 | tracker.SetVertex(vtx); | |
295 | }else { | |
296 | cerr<<"+++\n+++ Event header not found in galice.root:\n+++ Primary vertex in (0,0,0) [default]\n+++\n"; | |
297 | } | |
298 | ||
299 | // setup vertex constraint in the two tracking passes | |
300 | Int_t flags[2]; | |
301 | flags[0]=0; | |
302 | tracker.SetupFirstPass(flags); | |
303 | flags[0]=-1; | |
304 | tracker.SetupSecondPass(flags); | |
305 | ||
306 | rc=tracker.Clusters2Tracks(in,out); | |
307 | ||
308 | } | |
309 | ||
310 | delete gAlice; gAlice=0; | |
311 | in->Close(); | |
312 | in2->Close(); | |
313 | out->Close(); | |
314 | ||
315 | gBenchmark->Stop(name); | |
316 | gBenchmark->Show(name); | |
317 | ||
318 | return rc; | |
319 | } | |
320 | ||
321 | ||
322 | Int_t ITSMakeRefFile(const Char_t *galice, const Char_t *inname, | |
323 | const Char_t *outname, Int_t n) { | |
324 | ||
325 | ||
326 | cerr<<"\n*******************************************************************\n"; | |
327 | ||
328 | Int_t rc=0; | |
329 | const Char_t *name="ITSMakeRefFile"; | |
330 | cerr<<'\n'<<name<<"...\n"; | |
331 | gBenchmark->Start(name); | |
332 | ||
333 | ||
334 | TFile *out = TFile::Open(outname,"recreate"); | |
335 | TFile *trk = TFile::Open(inname); | |
336 | TFile *kin = TFile::Open(galice); | |
337 | ||
338 | ||
339 | // Get gAlice object from file | |
340 | if(!(gAlice=(AliRun*)kin->Get("gAlice"))) { | |
341 | cerr<<"gAlice has not been found on galice.root !\n"; | |
342 | return 1; | |
343 | } | |
344 | ||
345 | ||
346 | ||
347 | Int_t label; | |
348 | TParticle *Part; | |
349 | TParticle *Mum; | |
350 | static RECTRACK rectrk; | |
351 | ||
352 | ||
353 | for(Int_t event=0; event<n; event++){ | |
354 | ||
355 | AliITStrackV2 *itstrack=0; | |
356 | ||
357 | Int_t nparticles=gAlice->GetEvent(event); | |
358 | ||
359 | trk->cd(); | |
360 | ||
361 | // Tree with ITS tracks | |
362 | char tname[100]; | |
363 | sprintf(tname,"TreeT_ITS_%d",event); | |
364 | ||
365 | TTree *tracktree=(TTree*)trk->Get(tname); | |
366 | TBranch *tbranch=tracktree->GetBranch("tracks"); | |
367 | Int_t nentr=(Int_t)tracktree->GetEntries(); | |
368 | ||
369 | // Tree for true track parameters | |
370 | char ttname[100]; | |
371 | sprintf(ttname,"Tree_Ref_%d",event); | |
372 | TTree *reftree = new TTree(ttname,"Tree with true track params"); | |
373 | reftree->Branch("rectracks",&rectrk,"lab/I:pdg:Vx/F:Vy:Vz:Px:Py:Pz"); | |
374 | ||
375 | for (Int_t i=0; i<nentr; i++) { | |
376 | itstrack=new AliITStrackV2; | |
377 | tbranch->SetAddress(&itstrack); | |
378 | tracktree->GetEvent(i); | |
379 | label = TMath::Abs(itstrack->GetLabel()); | |
380 | ||
381 | Part = (TParticle*)gAlice->Particle(label); | |
382 | rectrk.lab=label; | |
383 | rectrk.pdg=Part->GetPdgCode(); | |
384 | if(Part->GetFirstMother()>=0) { | |
385 | Mum = (TParticle*)gAlice->Particle(Part->GetFirstMother()); | |
386 | rectrk.mumpdg=Mum->GetPdgCode(); | |
387 | } else { | |
388 | rectrk.mumpdg=-1; | |
389 | } | |
390 | rectrk.Vx=Part->Vx(); | |
391 | rectrk.Vy=Part->Vy(); | |
392 | rectrk.Vz=Part->Vz(); | |
393 | rectrk.Px=Part->Px(); | |
394 | rectrk.Py=Part->Py(); | |
395 | rectrk.Pz=Part->Pz(); | |
396 | ||
397 | reftree->Fill(); | |
398 | } | |
399 | ||
400 | out->cd(); | |
401 | reftree->Write(); | |
402 | ||
403 | } | |
404 | ||
405 | trk->Close(); | |
406 | kin->Close(); | |
407 | out->Close(); | |
408 | ||
409 | gBenchmark->Stop(name); | |
410 | gBenchmark->Show(name); | |
411 | ||
412 | ||
413 | return rc; | |
414 | ||
415 | } | |
416 | ||
417 | ||
418 | ||
419 | ||
420 | ||
421 | ||
422 | ||
423 | ||
424 | ||
425 | ||
426 | ||
427 |