]>
Commit | Line | Data |
---|---|---|
5a5a1232 | 1 | // $Header$ |
2 | ||
3 | #include "VSDCreator.h" | |
4 | ||
5 | #include <Reve/TTreeTools.h> | |
6 | ||
7 | #include <AliStack.h> | |
8 | #include <AliITSLoader.h> | |
9 | #include <AliTPCTrackHitsV2.h> | |
10 | #include <AliPDG.h> | |
11 | #include <AliHit.h> | |
12 | #include <AliSimDigits.h> | |
13 | #include <AliKalmanTrack.h> | |
14 | #include <AliESD.h> | |
15 | #include <AliESDV0MI.h> | |
16 | #include <AliTPCclusterMI.h> | |
17 | #include <AliTPCClustersRow.h> | |
18 | #include <AliITS.h> | |
19 | #include <AliITSclusterV2.h> | |
20 | #include <AliTrackReference.h> | |
21 | #include <AliESDkink.h> | |
22 | ||
23 | #include <AliRun.h> | |
24 | #include <AliTPCParam.h> | |
25 | ||
26 | #include <TSystem.h> | |
27 | #include <TFile.h> | |
28 | ||
29 | using namespace Reve; | |
30 | using namespace Alieve; | |
31 | ||
32 | using namespace std; | |
33 | ||
34 | //______________________________________________________________________ | |
35 | // VSDCreator | |
36 | // | |
37 | ||
38 | ClassImp(VSDCreator) | |
39 | ||
40 | void VSDCreator::init() | |
41 | { | |
42 | mKineType = KT_Standard; | |
43 | mDataDir = "."; | |
44 | ||
45 | mTPCHitRes = 2; | |
46 | mTRDHitRes = 2; | |
47 | ||
48 | pRunLoader = 0; | |
49 | ||
50 | // Particles not in ROOT's PDG database occuring in ALICE | |
51 | AliPDG::AddParticlesToPdgDataBase(); | |
52 | { | |
53 | TDatabasePDG *pdgDB = TDatabasePDG::Instance(); | |
54 | // const Int_t kspe=50000000; | |
55 | const Int_t kion=10000000; | |
56 | ||
57 | const Double_t kAu2Gev=0.9314943228; | |
58 | const Double_t khSlash = 1.0545726663e-27; | |
59 | const Double_t kErg2Gev = 1/1.6021773349e-3; | |
60 | const Double_t khShGev = khSlash*kErg2Gev; | |
61 | const Double_t kYear2Sec = 3600*24*365.25; | |
62 | ||
63 | pdgDB->AddParticle("Deuteron","Deuteron",2*kAu2Gev+8.071e-3,kTRUE, | |
64 | 0,1,"Ion",kion+10020); | |
65 | pdgDB->AddParticle("Triton","Triton",3*kAu2Gev+14.931e-3,kFALSE, | |
66 | khShGev/(12.33*kYear2Sec),1,"Ion",kion+10030); | |
67 | pdgDB->AddParticle("Alpha","Alpha",4*kAu2Gev+2.424e-3,kTRUE, | |
68 | khShGev/(12.33*kYear2Sec),2,"Ion",kion+20040); | |
69 | pdgDB->AddParticle("HE3","HE3",3*kAu2Gev+14.931e-3,kFALSE, | |
70 | 0,2,"Ion",kion+20030); | |
71 | } | |
72 | ||
73 | // AliKalmanTrack::SetConvConst(1); | |
74 | } | |
75 | ||
76 | VSDCreator::VSDCreator() | |
77 | { | |
78 | init(); | |
79 | } | |
80 | ||
81 | VSDCreator::VSDCreator(const Text_t* name, const Text_t* title) : | |
82 | VSD(name, title) | |
83 | { | |
84 | init(); | |
85 | } | |
86 | ||
87 | /**************************************************************************/ | |
88 | ||
89 | void VSDCreator::CreateVSD(const Text_t* data_dir, Int_t event, | |
90 | const Text_t* vsd_file) | |
91 | { | |
92 | static const Exc_t eH("VSDCreator::CreateVSD "); | |
93 | ||
94 | mDataDir = data_dir; | |
95 | mEvent = event; | |
96 | ||
97 | string galice_file (Form("%s/galice.root", mDataDir.Data())); | |
98 | ||
99 | if(mDebugLevel > 0) | |
100 | printf("%s opening %s \n", eH.Data(), galice_file.c_str()); | |
101 | ||
102 | if(gSystem->AccessPathName(galice_file.c_str(), kReadPermission)) { | |
103 | throw(eH + "Can not read file '" + galice_file + "'."); | |
104 | } | |
105 | pRunLoader = AliRunLoader::Open(galice_file.c_str()); | |
106 | if(pRunLoader == 0) | |
107 | throw(eH + "AliRunLoader::Open failed."); | |
108 | ||
109 | pRunLoader->LoadgAlice(); | |
110 | Int_t status = pRunLoader->GetEvent(mEvent); | |
111 | if(status) | |
112 | throw(eH + Form("GetEvent(%d) failed, exit code %s.", mEvent, status)); | |
113 | ||
114 | if(mDebugLevel > 0) | |
115 | printf("%s open seems ok. Now loading sim data.\n", eH.Data()); | |
116 | ||
117 | pRunLoader->LoadHeader(); | |
118 | pRunLoader->LoadKinematics(); | |
119 | pRunLoader->LoadTrackRefs(); | |
120 | pRunLoader->LoadHits(); | |
121 | ||
122 | // GledNS::PushFD(); | |
123 | ||
124 | if(mDebugLevel > 0) | |
125 | printf("%s opening output VSD.\n", eH.Data()); | |
126 | ||
127 | TFile* file = TFile::Open(vsd_file, "RECREATE", "ALICE VisualizationDataSummary"); | |
128 | mDirectory = new TDirectory("Event0", ""); | |
129 | ||
130 | if(mDebugLevel > 0) | |
131 | printf("%s creating trees now ...\n", eH.Data()); | |
132 | ||
133 | CreateTrees(); | |
134 | ||
135 | if(mDebugLevel > 0) | |
136 | printf("%s trees created, closing files.\n", eH.Data()); | |
137 | ||
138 | // file->Write(); | |
139 | file->Close(); | |
140 | delete file; | |
141 | mDirectory =0; | |
142 | ||
143 | //GledNS::PopFD(); | |
144 | ||
145 | // clean after the VSD data was sucessfuly written | |
146 | mTreeK = 0; | |
147 | mTreeH = 0; | |
148 | //mTreeTR = 0; | |
149 | mTreeC = 0; | |
150 | mTreeV0 = 0; | |
151 | mTreeKK = 0; | |
152 | mTreeR = 0; | |
153 | mTreeGI = 0; | |
154 | ||
155 | pRunLoader->UnloadAll(); | |
156 | delete pRunLoader; | |
157 | if(gAlice) { | |
158 | delete gAlice; gAlice = 0; | |
159 | } | |
160 | pRunLoader = 0; | |
161 | ||
162 | if(mDebugLevel > 0) | |
163 | printf("%s all done.\n", eH.Data()); | |
164 | } | |
165 | ||
166 | void VSDCreator::CreateTrees() | |
167 | { | |
168 | static const Exc_t eH("VSDCreator::CreateTrees "); | |
169 | ||
170 | if(mDirectory == 0) | |
171 | throw(eH + "output directory not set."); | |
172 | ||
173 | try { | |
174 | if(mDebugLevel > 1) | |
175 | printf("%s ConvertKinematics.\n", eH.Data()); | |
176 | ConvertKinematics(); | |
177 | } catch(Exc_t& exc) { WarnCaller(exc); } | |
178 | ||
179 | try { | |
180 | if(mDebugLevel > 1) | |
181 | printf("%s ConvertHits.\n", eH.Data()); | |
182 | ConvertHits(); | |
183 | } catch(Exc_t& exc) { WarnCaller(exc); } | |
184 | ||
185 | try { | |
186 | if(mDebugLevel > 1) | |
187 | printf("%s ConvertClusters.\n", eH.Data()); | |
188 | ConvertClusters(); | |
189 | } catch(Exc_t& exc) { WarnCaller(exc); } | |
190 | ||
191 | ||
192 | printf("############# EXITING, had incompatible ESDTrack class #####\n"); | |
193 | goto end_geninfo_processing; | |
194 | ||
195 | ||
196 | try { | |
197 | if(mDebugLevel > 1) | |
198 | printf("%s ConvertRecTracks.\n", eH.Data()); | |
199 | ConvertRecTracks(); | |
200 | } catch(Exc_t& exc) { | |
201 | WarnCaller(exc + " Skipping V0 extraction."); | |
202 | goto end_esd_processing; | |
203 | } | |
204 | ||
205 | try { | |
206 | if(mDebugLevel > 1) | |
207 | printf("%s ConvertV0.\n", eH.Data()); | |
208 | ConvertV0(); | |
209 | } catch(Exc_t& exc) { WarnCaller(exc); } | |
210 | ||
211 | try { | |
212 | if(mDebugLevel > 1) | |
213 | printf("%s ConvertKinks.\n", eH.Data()); | |
214 | ConvertKinks(); | |
215 | } catch(Exc_t& exc) { WarnCaller(exc); } | |
216 | ||
217 | end_esd_processing: | |
218 | ||
219 | try { | |
220 | if(mDebugLevel > 1) | |
221 | printf("%s ConvertGenInfo.\n", eH.Data()); | |
222 | ConvertGenInfo(); | |
223 | } catch(Exc_t& exc) { WarnCaller(exc); } | |
224 | ||
225 | end_geninfo_processing: | |
226 | return; | |
227 | } | |
228 | ||
229 | /**************************************************************************/ | |
230 | // Kinematics | |
231 | /**************************************************************************/ | |
232 | ||
233 | void VSDCreator::ConvertKinematics() | |
234 | { | |
235 | static const Exc_t eH("VSDCreator::ConvertKinematics "); | |
236 | ||
237 | if(mTreeK != 0) | |
238 | throw (eH + "kinematics already converted"); | |
239 | ||
240 | AliStack* stack = pRunLoader->Stack(); | |
241 | if(stack == 0) | |
242 | throw(eH + "stack is null."); | |
243 | ||
244 | mDirectory->cd(); | |
245 | mTreeK = new TTree("Kinematics", "TParticles sorted by Label"); | |
246 | ||
247 | Int_t nentries = stack->GetNtrack(); | |
248 | vector<MCTrack> vmc(nentries); | |
249 | for (Int_t idx=0; idx<nentries; idx++) { | |
250 | TParticle* tp = stack->Particle(idx); | |
251 | vmc[idx] = *tp; | |
252 | vmc[idx].label = idx; | |
253 | } | |
254 | ||
255 | // read track refrences | |
256 | TTree* mTreeTR = pRunLoader->TreeTR(); | |
257 | ||
258 | if(mTreeTR == 0) { | |
259 | WarnCaller(eH + "no TrackRefs; some data will not be available."); | |
260 | } else { | |
261 | TClonesArray* RunArrayTR = 0; | |
262 | mTreeTR->SetBranchAddress("AliRun", &RunArrayTR); | |
263 | ||
264 | Int_t nPrimaries = (Int_t) mTreeTR->GetEntries(); | |
265 | for (Int_t iPrimPart = 0; iPrimPart<nPrimaries; iPrimPart++) { | |
266 | // printf("START mTreeTR->GetEntry(%d) \n",iPrimPart); | |
267 | mTreeTR->GetEntry(iPrimPart); | |
268 | // printf("END mTreeTR->GetEntry(%d) \n",iPrimPart); | |
269 | ||
270 | for (Int_t iTrackRef = 0; iTrackRef < RunArrayTR->GetEntriesFast(); iTrackRef++) { | |
271 | AliTrackReference *trackRef = (AliTrackReference*)RunArrayTR->At(iTrackRef); | |
272 | Int_t track = trackRef->GetTrack(); | |
273 | if(track < nentries && track > 0){ | |
274 | MCTrack& mct = vmc[track]; | |
275 | if(trackRef->TestBit(kNotDeleted)) { | |
276 | mct.decayed = true; | |
277 | mct.t_decay = trackRef->GetTime(); | |
278 | mct.V_decay.x = trackRef->X(); | |
279 | mct.V_decay.y = trackRef->Y(); | |
280 | mct.V_decay.z = trackRef->Z(); | |
281 | mct.P_decay.x = trackRef->Px(); | |
282 | mct.P_decay.y = trackRef->Py(); | |
283 | mct.P_decay.z = trackRef->Pz(); | |
284 | if(TMath::Abs(mct.GetPdgCode()) == 11) | |
285 | mct.decayed = false; // a bug in TreeTR | |
286 | } | |
287 | } | |
288 | } | |
289 | } | |
290 | } | |
291 | ||
292 | mTreeK->Branch("K", "Reve::MCTrack", &mpK, fBuffSize); | |
293 | ||
294 | printf("sizeofvmc = %d\n", vmc.size()); | |
295 | for(vector<MCTrack>::iterator k=vmc.begin(); k!=vmc.end(); ++k) { | |
296 | MCTrack& mct = *k; | |
297 | mK = mct; | |
298 | ||
299 | TParticle* m = &mct; | |
300 | Int_t mi = mct.label; | |
301 | int cnt = 0; | |
302 | while(m->GetMother(0) != -1) { | |
303 | if(cnt > 100) { | |
304 | printf("cnt %d mi=%d, mo=%d\n", cnt, mi, m->GetMother(0)); | |
305 | } | |
306 | mi = m->GetMother(0); | |
307 | m = &vmc[mi]; | |
308 | ++cnt; | |
309 | } | |
310 | mK.eva_label = mi; | |
311 | ||
312 | mTreeK->Fill(); | |
313 | } | |
314 | ||
315 | mTreeK->BuildIndex("label"); | |
316 | } | |
317 | ||
318 | /**************************************************************************/ | |
319 | // Hits | |
320 | /**************************************************************************/ | |
321 | ||
322 | namespace { | |
323 | ||
324 | struct Detector | |
325 | { | |
326 | const char* name; | |
327 | const char* hitbranch; | |
328 | unsigned char detidx; | |
329 | }; | |
330 | ||
331 | Detector detects[] = { | |
332 | { "ITS", "AliITShit", 0 }, | |
333 | { "TPC", "AliTPCTrackHitsV2", 1 }, | |
334 | { "TRD", "AliTRDhit", 2 }, | |
335 | { "TOF", "AliTOFhit", 3 } | |
336 | // { "RICH", "AliRICHhit", 4 }, | |
337 | }; | |
338 | ||
339 | } | |
340 | ||
341 | /**************************************************************************/ | |
342 | ||
343 | void VSDCreator::ConvertHits() | |
344 | { | |
345 | static const Exc_t eH("VSDCreator::ConvertHits "); | |
346 | ||
347 | if(mTreeH != 0) | |
348 | throw(eH + "hits already converted."); | |
349 | ||
350 | mDirectory->cd(); | |
351 | mTreeH = new TTree("Hits", "Combined detector hits."); | |
352 | mTreeH->Branch("H", "Reve::Hit", &mpH, fBuffSize); | |
353 | ||
354 | map<Int_t, Int_t> hmap; | |
355 | // parameters for ITS, TPC hits filtering | |
356 | Float_t x,y,z, x1,y1,z1; | |
357 | Float_t tpc_sqr_res = mTPCHitRes*mTPCHitRes; | |
358 | Float_t trd_sqr_res = mTRDHitRes*mTRDHitRes; | |
359 | ||
360 | int l=0; | |
361 | // load hits from the rest of detectors | |
362 | while(detects[l].name != 0) { | |
363 | Detector& det = detects[l++]; | |
364 | ||
365 | switch(det.detidx) { | |
366 | case 1: { | |
367 | Int_t count = 0; | |
368 | TTree* treeh = pRunLoader->GetTreeH(det.name, false); | |
369 | if(treeh == 0) { | |
370 | WarnCaller(eH + "no hits for "+ det.name +"."); | |
371 | continue; | |
372 | } | |
373 | AliTPCTrackHitsV2 hv2, *_hv2=&hv2; | |
374 | treeh->SetBranchAddress("TPC2", &_hv2); | |
375 | Int_t np = treeh->GetEntries(); | |
376 | for(Int_t i=0; i<np; i++){ | |
377 | treeh->GetEntry(i); | |
378 | Int_t eva_idx = np -i -1; | |
379 | if (hv2.First() == 0) continue; | |
380 | x = y = z = 0; | |
381 | do { | |
382 | AliHit* ah = hv2.GetHit(); | |
383 | x1=ah->X();y1=ah->Y();z1=ah->Z(); | |
384 | if((x-x1)*(x-x1)+(y-y1)*(y-y1)+(z-z1)*(z-z1) > tpc_sqr_res) { | |
385 | mH.det_id = det.detidx; | |
386 | mH.subdet_id = 0; | |
387 | mH.label = ah->Track(); | |
388 | mH.eva_label = eva_idx; | |
389 | mH.V.x = x1; mH.V.y = y1; mH.V.z = z1; | |
390 | mTreeH->Fill(); | |
391 | hmap[mH.label]++; | |
392 | x = x1; y = y1; z = z1; | |
393 | count++; | |
394 | } | |
395 | } while (hv2.Next()); | |
396 | } | |
397 | // printf("%d entries in TPChits \n",count); | |
398 | break; | |
399 | } | |
400 | default: { | |
401 | TTree* treeh = pRunLoader->GetTreeH(det.name, false); | |
402 | if(treeh == 0) { | |
403 | WarnCaller(eH + "no hits for "+ det.name +"."); | |
404 | continue; | |
405 | } | |
406 | TClonesArray *arr = new TClonesArray(det.hitbranch); | |
407 | treeh->SetBranchAddress(det.name, &arr); | |
408 | Int_t np = treeh->GetEntries(); | |
409 | // in TreeH files hits are grouped in clones arrays | |
410 | // each eva particle has its own clone array | |
411 | for (Int_t i=0; i<np; i++) { | |
412 | treeh->GetEntry(i); | |
413 | Int_t eva_idx = np -i -1; | |
414 | Int_t nh=arr->GetEntriesFast(); | |
415 | x = y = z = 0; | |
416 | // printf("%d entry %d hits for primary %d \n", i, nh, eva_idx); | |
417 | for (Int_t j=0; j<nh; j++) { | |
418 | AliHit* ali_hit = (AliHit*)arr->UncheckedAt(j); | |
419 | mH.det_id = det.detidx; | |
420 | mH.subdet_id = 0; | |
421 | mH.label = ali_hit->GetTrack(); | |
422 | mH.eva_label = eva_idx; | |
423 | mH.V.Set(ali_hit->X(), ali_hit->Y(), ali_hit->Z()); | |
424 | if(det.detidx == 2) { | |
425 | x1=ali_hit->X();y1=ali_hit->Y();z1=ali_hit->Z(); | |
426 | if((x-x1)*(x-x1)+(y-y1)*(y-y1)+(z-z1)*(z-z1) < trd_sqr_res) continue; | |
427 | x=x1; y=y1; z=z1; | |
428 | } | |
429 | hmap[mH.label]++; | |
430 | mTreeH->Fill(); | |
431 | } | |
432 | } | |
433 | delete arr; | |
434 | break; | |
435 | } // end default | |
436 | } // end switch | |
437 | } // end while | |
438 | ||
439 | ||
440 | //set geninfo | |
441 | for(map<Int_t, Int_t>::iterator j=hmap.begin(); j!=hmap.end(); ++j) { | |
442 | GetGeninfo(j->first)->n_hits += j->second; | |
443 | } | |
444 | ||
445 | } | |
446 | ||
447 | /**************************************************************************/ | |
448 | // Clusters | |
449 | /**************************************************************************/ | |
450 | ||
451 | void VSDCreator::ConvertClusters() | |
452 | { | |
453 | static const Exc_t eH("VSDCreator::ConvertClusters "); | |
454 | ||
455 | if(mTreeC != 0) | |
456 | throw(eH + "clusters already converted."); | |
457 | ||
458 | mDirectory->cd(); | |
459 | mTreeC = new TTree("Clusters", "rec clusters"); | |
460 | mTreeC->Branch("C", "Reve::Cluster", &mpC, fBuffSize); | |
461 | ||
462 | try { | |
463 | ConvertITSClusters(); | |
464 | } catch(Exc_t& exc) { WarnCaller(exc); } | |
465 | ||
466 | try { | |
467 | ConvertTPCClusters(); | |
468 | } catch(Exc_t& exc) { WarnCaller(exc); } | |
469 | } | |
470 | ||
471 | /**************************************************************************/ | |
472 | ||
473 | void VSDCreator::ConvertTPCClusters() | |
474 | { | |
475 | static const Exc_t eH("VSDCreator::ConvertTPCClusters "); | |
476 | ||
477 | auto_ptr<TFile> f | |
478 | ( TFile::Open(Form("%s/TPC.RecPoints.root", mDataDir.Data())) ); | |
479 | if(!f.get()) | |
480 | throw(eH + "can not open 'TPC.RecPoints.root' file."); | |
481 | ||
482 | auto_ptr<TDirectory> d | |
483 | ( (TDirectory*) f->Get(Form("Event%d", mEvent)) ); | |
484 | if(!d.get()) | |
485 | throw(eH + Form("event directory '%d' not found.", 0)); | |
486 | ||
487 | auto_ptr<TTree> tree( (TTree*) d->Get("TreeR") ); | |
488 | if(!tree.get()) | |
489 | throw(eH + "'TreeR' not found."); | |
490 | ||
491 | auto_ptr<AliTPCParam> par( GetTpcParam(eH) ); | |
492 | ||
493 | AliTPCClustersRow clrow, *_clrow=&clrow; | |
494 | AliTPCclusterMI *cl; | |
495 | _clrow->SetClass("AliTPCclusterMI"); | |
496 | tree->SetBranchAddress("Segment", &_clrow); | |
497 | ||
498 | // count clusters | |
499 | Int_t nClusters = 0; | |
500 | Int_t n_ent = tree->GetEntries(); | |
501 | for (Int_t n=0; n<n_ent; n++) { | |
502 | tree->GetEntry(n); | |
503 | nClusters += _clrow->GetArray()->GetEntriesFast(); | |
504 | } | |
505 | ||
506 | // calculate xyz for a cluster and add it to container | |
507 | Double_t x,y,z; | |
508 | Float_t cs, sn, tmp; | |
509 | map<Int_t, Int_t> cmap; | |
510 | ||
511 | for (Int_t n=0; n<tree->GetEntries(); n++) { | |
512 | tree->GetEntry(n); | |
513 | Int_t ncl = _clrow->GetArray()->GetEntriesFast(); | |
514 | if(ncl > 0) { | |
515 | Int_t sec,row; | |
516 | par->AdjustSectorRow(_clrow->GetID(),sec,row); | |
517 | while (ncl--) { | |
518 | if(_clrow->GetArray()) { | |
519 | // cl = new AliTPCclusterMI(*(AliTPCclusterMI*)_clrow->GetArray()->UncheckedAt(ncl)); | |
520 | cl = (AliTPCclusterMI*)_clrow->GetArray()->UncheckedAt(ncl); | |
521 | if(cl->GetLabel(0) >= 0){ | |
522 | x = par->GetPadRowRadii(sec,row); y = cl->GetY(); z = cl->GetZ(); | |
523 | par->AdjustCosSin(sec,cs,sn); | |
524 | tmp = x*cs-y*sn; y= x*sn+y*cs; x=tmp; | |
525 | ||
526 | mC.det_id = 1; | |
527 | mC.subdet_id = 0; | |
528 | mC.label[0] = cl->GetLabel(0); | |
529 | mC.label[1] = cl->GetLabel(1); | |
530 | mC.label[2] = cl->GetLabel(2); | |
531 | mC.V.Set(x, y, z); | |
532 | ||
533 | mTreeC->Fill(); | |
534 | { int i = 0; | |
535 | while(i < 3 && mC.label[i]) | |
536 | cmap[mC.label[i++]]++; | |
537 | } | |
538 | } | |
539 | } | |
540 | } | |
541 | } | |
542 | } | |
543 | //set geninfo | |
544 | for(map<Int_t, Int_t>::iterator j=cmap.begin(); j!=cmap.end(); ++j) { | |
545 | GetGeninfo(j->first)->n_clus += j->second; | |
546 | } | |
547 | } | |
548 | ||
549 | /**************************************************************************/ | |
550 | ||
551 | void VSDCreator::ConvertITSClusters() | |
552 | { | |
553 | static const Exc_t eH("VSDCreator::ConvertITSClusters "); | |
554 | ||
555 | auto_ptr<TFile> f | |
556 | ( TFile::Open(Form("%s/ITS.RecPoints.root", mDataDir.Data())) ); | |
557 | if(!f.get()) | |
558 | throw(eH + "can not open 'ITS.RecPoints.root' file."); | |
559 | ||
560 | auto_ptr<TDirectory> d | |
561 | ( (TDirectory*) f->Get(Form("Event%d", mEvent)) ); | |
562 | if(!d.get()) | |
563 | throw(eH + Form("event directory '%d' not found.", 0)); | |
564 | ||
565 | auto_ptr<TTree> tree( (TTree*) d->Get("TreeR") ); | |
566 | if(!tree.get()) | |
567 | throw(eH + "'TreeR' not found."); | |
568 | ||
569 | AliITSLoader* ITSld = (AliITSLoader*) pRunLoader->GetLoader("ITSLoader"); | |
570 | //AliITS* pITS = ITSld->GetITS(); | |
571 | AliITSgeom* geom = ITSld->GetITSgeom(); | |
572 | //AliITSgeom* geom = new AliITSgeom(); | |
573 | //geom->ReadNewFile("/home/aljam/ITSgeometry.det"); | |
574 | ||
575 | //printf("alice ITS geom %p \n",geom ); | |
576 | ||
577 | if(!geom) | |
578 | throw(eH + "can not find ITS geometry"); | |
579 | ||
580 | TClonesArray *arr = new TClonesArray("AliITSclusterV2"); | |
581 | tree->SetBranchAddress("Clusters", &arr); | |
582 | Int_t nmods = tree->GetEntries(); | |
583 | Float_t gc[3]; | |
584 | map<Int_t, Int_t> cmap; | |
585 | ||
586 | for (Int_t mod=0; mod<nmods; mod++) { | |
587 | tree->GetEntry(mod); | |
588 | Int_t nc=arr->GetEntriesFast(); | |
589 | for (Int_t j=0; j<nc; j++) { | |
590 | AliITSclusterV2* recp = (AliITSclusterV2*)arr->UncheckedAt(j); | |
591 | ||
592 | Double_t rot[9]; | |
593 | geom->GetRotMatrix(mod,rot); | |
594 | Int_t lay,lad,det; | |
595 | geom->GetModuleId(mod,lay,lad,det); | |
596 | Float_t tx,ty,tz; | |
597 | geom->GetTrans(lay,lad,det,tx,ty,tz); | |
598 | ||
599 | Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi(); | |
600 | Double_t phi1=TMath::Pi()/2+alpha; | |
601 | if(lay==1) phi1+=TMath::Pi(); | |
602 | ||
603 | Float_t cp=TMath::Cos(phi1), sp=TMath::Sin(phi1); | |
604 | Float_t r=tx*cp+ty*sp; | |
605 | gc[0]= r*cp - recp->GetY()*sp; | |
606 | gc[1]= r*sp + recp->GetY()*cp; | |
607 | gc[2]= recp->GetZ(); | |
608 | ||
609 | mC.det_id = 0; | |
610 | mC.subdet_id = 0; | |
611 | mC.label[0] = recp->GetLabel(0); | |
612 | mC.label[1] = recp->GetLabel(1); | |
613 | mC.label[2] = recp->GetLabel(2); | |
614 | mC.V.x = r*cp - recp->GetY()*sp; | |
615 | mC.V.y = r*sp + recp->GetY()*cp; | |
616 | mC.V.z = recp->GetZ(); | |
617 | mTreeC->Fill(); | |
618 | { int i = 0; | |
619 | while(i < 3 && mC.label[i]) | |
620 | cmap[mC.label[i++]]++; | |
621 | } | |
622 | } | |
623 | ||
624 | for(map<Int_t, Int_t>::iterator j=cmap.begin(); j!=cmap.end(); ++j) { | |
625 | GetGeninfo(j->first)->n_clus += j->second; | |
626 | } | |
627 | } | |
628 | delete arr; | |
629 | } | |
630 | ||
631 | /**************************************************************************/ | |
632 | // ESD | |
633 | /**************************************************************************/ | |
634 | ||
635 | void VSDCreator::ConvertRecTracks() | |
636 | { | |
637 | static const Exc_t eH("VSDCreator::ConvertRecTracks "); | |
638 | ||
639 | if(mTreeR != 0) | |
640 | throw(eH + "tracks already converted."); | |
641 | ||
642 | mDirectory->cd(); | |
643 | mTreeR = new TTree("RecTracks", "rec tracks"); | |
644 | ||
645 | mTreeR->Branch("R", "Reve::RecTrack", &mpR, 512*1024,1); | |
646 | ||
647 | TFile f(Form("%s/AliESDs.root", mDataDir.Data())); | |
648 | if(!f.IsOpen()) | |
649 | throw(eH + "no AliESDs.root file."); | |
650 | ||
651 | TTree* tree = (TTree*) f.Get("esdTree"); | |
652 | if (tree == 0) | |
653 | throw(eH + "no esdTree."); | |
654 | ||
655 | ||
656 | AliESD *fEvent=0; | |
657 | tree->SetBranchAddress("ESD", &fEvent); | |
658 | tree->GetEntry(mEvent); | |
659 | ||
660 | ||
661 | // reconstructed tracks | |
662 | AliESDtrack* esd_t; | |
663 | Double_t dbuf[3]; | |
664 | for (Int_t n=0; n<fEvent->GetNumberOfTracks(); n++) { | |
665 | esd_t = fEvent->GetTrack(n); | |
666 | ||
667 | mR.label = esd_t->GetLabel(); | |
668 | mR.status = (Int_t) esd_t->GetStatus(); | |
669 | mR.sign = (Int_t) esd_t->GetSign(); | |
670 | esd_t->GetXYZ(dbuf); mR.V.Set(dbuf); | |
671 | esd_t->GetPxPyPz(dbuf); mR.P.Set(dbuf); | |
672 | Double_t ep = esd_t->GetP(); | |
673 | mR.beta = ep/TMath::Sqrt(ep*ep + TMath::C()*TMath::C()*esd_t->GetMass()*esd_t->GetMass()); | |
674 | mTreeR->Fill(); | |
675 | } | |
676 | mTreeR->BuildIndex("label"); | |
677 | } | |
678 | ||
679 | /**************************************************************************/ | |
680 | ||
681 | void VSDCreator::ConvertV0() | |
682 | { | |
683 | static const Exc_t eH("VSDCreator::ConvertV0 "); | |
684 | ||
685 | if(mTreeV0 != 0) | |
686 | throw(eH + "V0 already converted."); | |
687 | ||
688 | mDirectory->cd(); | |
689 | mTreeV0 = new TTree("V0", "V0 points"); | |
690 | ||
691 | mTreeV0->Branch("V0", "Reve::RecV0", &mpV0, 512*1024,1); | |
692 | ||
693 | TFile f(Form("%s/AliESDs.root", mDataDir.Data())); | |
694 | if(!f.IsOpen()){ | |
695 | throw(eH + "no AliESDs.root file."); | |
696 | } | |
697 | ||
698 | TTree* tree = (TTree*) f.Get("esdTree"); | |
699 | if (tree == 0) | |
700 | throw(eH + "no esdTree."); | |
701 | ||
702 | AliESD *fEvent=0; | |
703 | tree->SetBranchAddress("ESD", &fEvent); | |
704 | tree->GetEntry(mEvent); | |
705 | ||
706 | for (Int_t n =0; n< fEvent->GetNumberOfV0MIs(); n++) { | |
707 | AliESDV0MI* av = fEvent->GetV0MI(n); | |
708 | Double_t pos[3]; | |
709 | ||
710 | mV0.status = av->GetStatus(); | |
711 | // Point of closest approach | |
712 | mV0.V_ca.x = av->GetXr(0); | |
713 | mV0.V_ca.y = av->GetXr(1); | |
714 | mV0.V_ca.z = av->GetXr(2); | |
715 | // set birth vertex of neutral particle | |
716 | { Double_t p[3]; | |
717 | av->GetXYZ(p[0], p[1], p[2]); | |
718 | mV0.V0_birth.Set(p); | |
719 | } | |
720 | ||
721 | // momentum and position of negative particle | |
722 | mV0.P_neg.Set(av->GetPMp()); | |
723 | av->GetParamM()->GetXYZ(pos); | |
724 | mV0.V_neg.Set(pos); | |
725 | ||
726 | // momentum and position of positive particle | |
727 | mV0.P_pos.Set(av->GetPPp()); | |
728 | av->GetParamP()->GetXYZ(pos); | |
729 | mV0.V_pos.Set(pos); | |
730 | ||
731 | mV0.label = 0; // !!!! mother label unknown | |
732 | mV0.pdg = av->GetPdgCode(); | |
733 | // daughter indices | |
734 | mV0.d_label[0] = av->GetLab(0); | |
735 | mV0.d_label[1] = av->GetLab(1); | |
736 | ||
737 | // printf("V0 convert labels(%d,%d) index(%d,%d)\n", | |
738 | // av->GetLab(0), av->GetLab(1), | |
739 | // av->GetIndex(0), av->GetIndex(1)); | |
740 | ||
741 | mTreeV0->Fill(); | |
742 | } | |
743 | // if(fEvent->GetNumberOfV0MIs()) mTreeV0->BuildIndex("label"); | |
744 | } | |
745 | ||
746 | /**************************************************************************/ | |
747 | ||
748 | void VSDCreator::ConvertKinks() | |
749 | { | |
750 | static const Exc_t eH("VSDCreator::ConvertKinks "); | |
751 | ||
752 | if(mTreeKK != 0) | |
753 | throw(eH + "Kinks already converted."); | |
754 | ||
755 | mDirectory->cd(); | |
756 | mTreeKK = new TTree("Kinks", "ESD Kinks"); | |
757 | ||
758 | mTreeKK->Branch("KK", "Reve::RecKink", &mpKK, fBuffSize); | |
759 | ||
760 | TFile f(Form("%s/AliESDs.root", mDataDir.Data())); | |
761 | if(!f.IsOpen()){ | |
762 | throw(eH + "no AliESDs.root file."); | |
763 | } | |
764 | ||
765 | TTree* tree = (TTree*) f.Get("esdTree"); | |
766 | if (tree == 0) | |
767 | throw(eH + "no esdTree."); | |
768 | ||
769 | AliESD *fEvent=0; | |
770 | tree->SetBranchAddress("ESD", &fEvent); | |
771 | tree->GetEntry(mEvent); | |
772 | ||
773 | // printf("CONVERT KINK Read %d entries in tree kinks \n", fEvent->GetNumberOfKinks()); | |
774 | for (Int_t n =0; n< fEvent->GetNumberOfKinks(); n++) { | |
775 | AliESDkink* kk = fEvent->GetKink(n); | |
776 | ||
777 | Double_t pos[3]; | |
778 | ||
779 | ||
780 | mKK.label = kk->GetLabel(0); | |
781 | mKK.status = Int_t(kk->GetStatus(1) << 8 + kk->GetStatus(2)); | |
782 | ||
783 | // reconstructed kink position | |
784 | mKK.label_sec = kk->GetLabel(1); | |
785 | mKK.V_kink.Set(kk->GetPosition()); | |
786 | ||
787 | const AliExternalTrackParam& tp_mother = kk->RefParamMother(); | |
788 | // momentum and position of mother | |
789 | tp_mother.GetPxPyPz(pos); | |
790 | mKK.P.Set(pos); | |
791 | tp_mother.GetXYZ(pos); | |
792 | mKK.V.Set(pos); | |
793 | const Double_t* par = tp_mother.GetParameter(); | |
794 | // printf("KINK Pt %f, %f \n",1/tp_mother.Pt(),par[4] ); | |
795 | mKK.sign = (par[4] < 0) ? -1 : 1; | |
796 | ||
797 | const AliExternalTrackParam& tp_daughter = kk->RefParamDaughter(); | |
798 | // momentum and position of daughter | |
799 | tp_daughter.GetPxPyPz(pos); | |
800 | mKK.P_sec.Set(pos); | |
801 | tp_daughter.GetXYZ(pos); | |
802 | mKK.V_end.Set(pos); | |
803 | ||
804 | mTreeKK->Fill(); | |
805 | } | |
806 | if(fEvent->GetNumberOfKinks()) mTreeKK->BuildIndex("label"); | |
807 | } | |
808 | /**************************************************************************/ | |
809 | // GenInfo | |
810 | /**************************************************************************/ | |
811 | ||
812 | void VSDCreator::ConvertGenInfo() | |
813 | { | |
814 | static const Exc_t eH("VSDCreator::ConvertGenInfo "); | |
815 | ||
816 | if(mTreeGI != 0) | |
817 | throw(eH + "GI already converted."); | |
818 | ||
819 | mDirectory->cd(); | |
820 | mTreeGI = new TTree("GenInfo", "Objects prepared for cross querry"); | |
821 | ||
822 | GenInfo::Class()->IgnoreTObjectStreamer(true); | |
823 | mTreeGI->Branch("GI", "Reve::GenInfo", &mpGI, fBuffSize); | |
824 | mTreeGI->Branch("K.", "Reve::MCTrack", &mpK); | |
825 | mTreeGI->Branch("R.", "Reve::RecTrack", &mpR); | |
826 | ||
827 | for(map<Int_t, GenInfo*>::iterator j=mGenInfoMap.begin(); j!=mGenInfoMap.end(); ++j) { | |
828 | mGI = *(j->second); | |
829 | mGI.label = j->first; | |
830 | mTreeK->GetEntry(j->first); | |
831 | ||
832 | if(mTreeR) { | |
833 | Int_t re = mTreeR->GetEntryNumberWithIndex(j->first); | |
834 | if(re != -1) | |
835 | mGI.is_rec = true; | |
836 | } | |
837 | // Int_t has_v0 = mTreeV0->GetEntryNumberWithIndex(j->first); | |
838 | //if (has_v0 != -1) | |
839 | // mGI.has_V0 = true; | |
840 | if (mTreeKK) { | |
841 | Int_t has_kk = mTreeKK->GetEntryNumberWithIndex(j->first); | |
842 | if (has_kk != -1) | |
843 | mGI.has_kink = true; | |
844 | } | |
845 | mTreeGI->Fill(); | |
846 | } | |
847 | mGenInfoMap.clear(); | |
848 | } | |
849 | ||
850 | /**************************************************************************/ | |
851 | /**************************************************************************/ | |
852 | // Protected methods | |
853 | /**************************************************************************/ | |
854 | /**************************************************************************/ | |
855 | ||
856 | AliTPCParam* VSDCreator::GetTpcParam(const Exc_t& eh) | |
857 | { | |
858 | auto_ptr<TFile> fp( TFile::Open(Form("%s/galice.root", mDataDir.Data())) ); | |
859 | if(!fp.get()) | |
860 | throw(eh + "can not open 'galice.root' file."); | |
861 | AliTPCParam* par = (AliTPCParam *) fp->Get("75x40_100x60_150x60"); | |
862 | if(!par) | |
863 | throw(eh + "TPC data not found."); | |
864 | return par; | |
865 | } | |
866 | ||
867 | ||
868 | ||
869 | GenInfo* VSDCreator::GetGeninfo(Int_t label) | |
870 | { | |
871 | // printf("get_geninfo %d\n", label); | |
872 | GenInfo* gi; | |
873 | map<Int_t, GenInfo*>::iterator i = mGenInfoMap.find(label); | |
874 | if(i == mGenInfoMap.end()) { | |
875 | gi = new GenInfo(); | |
876 | mGenInfoMap[label] = gi; | |
877 | } else { | |
878 | gi = i->second; | |
879 | } | |
880 | return gi; | |
881 | } |