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