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