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