2 // Main authors: Matevz Tadel & Alja Mrak-Tadel: 2006, 2007
4 /**************************************************************************
5 * Copyright(c) 1998-2008, ALICE Experiment at CERN, all rights reserved. *
6 * See http://aliceinfo.cern.ch/Offline/AliRoot/License.html for *
7 * full copyright notice. *
8 **************************************************************************/
11 #include "AliEveVSDCreator.h"
13 #include "AliEveEventManager.h"
16 #include <AliITSLoader.h>
17 #include <AliTPCTrackHitsV2.h>
20 #include <AliESDEvent.h>
22 #include <AliTPCclusterMI.h>
23 #include <AliTPCClustersRow.h>
24 #include <AliITSclusterV2.h>
25 #include <AliESDkink.h>
26 #include <AliESDtrack.h>
29 #include <AliRunLoader.h>
30 #include <AliTPCParam.h>
32 #include <TDatabasePDG.h>
36 //______________________________________________________________________________
38 // Create VSD file from ALICE data.
40 ClassImp(AliEveVSDCreator)
42 AliEveVSDCreator::AliEveVSDCreator(const Text_t* name, const Text_t* title) :
54 // Particles not in ROOT's PDG database occuring in ALICE
55 AliPDG::AddParticlesToPdgDataBase();
58 /******************************************************************************/
60 void AliEveVSDCreator::CreateVSD(const Text_t* vsdFile)
62 // Create the VSD for current event in AliEveEventManager.
63 // Result is stored in directory "Event%04d".
65 static const TEveException kEH("AliEveVSDCreator::CreateVSD ");
67 AliEveEventManager::AssertGeometry();
69 fRunLoader = AliEveEventManager::AssertRunLoader();
72 printf("%s open seems ok. Now loading sim data.\n", kEH.Data());
74 fRunLoader->LoadHeader();
75 fRunLoader->LoadKinematics();
76 fRunLoader->LoadTrackRefs();
77 fRunLoader->LoadHits();
80 printf("%s opening output TEveVSD.\n", kEH.Data());
82 TFile* file = TFile::Open(vsdFile, "UPDATE", "ALICE Visualization Summary Data");
84 Int_t curEvent = AliEveEventManager::CurrentEventId();
86 TString eventDir; eventDir.Form("Event%04d", curEvent);
88 if (file->Get(eventDir))
90 Warning(kEH, "Event-dir '%s' existed -- replacing.", eventDir.Data());
91 file->Delete(eventDir + ";*");
94 fDirectory = new TDirectoryFile(eventDir, Form("Data for event %d", curEvent));
97 printf("%s creating trees now ...\n", kEH.Data());
102 printf("%s trees created, closing files.\n", kEH.Data());
109 // clean after the TEveVSD data was sucessfuly written
122 printf("%s all done.\n", kEH.Data());
125 void AliEveVSDCreator::CreateTrees()
127 // Create and fill the output trees by calling all the
128 // ConvertXyzz() functions.
129 // Exceptions from individual functions are displayed as warnings.
131 static const TEveException kEH("AliEveVSDCreator::CreateTrees ");
134 throw kEH + "output directory not set.";
139 printf("%sConvertKinematics.\n", kEH.Data());
142 catch(TEveException& exc)
144 Warning(kEH, "%s", exc.Data());
150 printf("%sConvertHits.\n", kEH.Data());
153 catch(TEveException& exc)
155 Warning(kEH, "%s", exc.Data());
161 printf("%sConvertClusters.\n", kEH.Data());
164 catch(TEveException& exc)
166 Warning(kEH, "%s", exc.Data());
172 printf("%sConvertRecTracks.\n", kEH.Data());
175 catch(TEveException& exc)
177 Warning(kEH, "Executing %s skipping AliEveV0 extraction.", exc.Data());
178 goto end_esd_processing;
181 // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
183 Warning(kEH, "Explicitly abandoning further conversion.");
186 // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
191 printf("%sConvertV0.\n", kEH.Data());
194 catch(TEveException& exc)
202 printf("%sConvertKinks.\n", kEH.Data());
205 catch(TEveException& exc)
215 printf("%sConvertGenInfo.\n", kEH.Data());
218 catch(TEveException& exc)
220 Warning(kEH, "%s", exc.Data());
226 /******************************************************************************/
228 /******************************************************************************/
230 void AliEveVSDCreator::ConvertKinematics()
232 // Convert kinematics.
233 // Track references are not stored, they should be.
235 static const TEveException kEH("AliEveVSDCreator::ConvertKinematics ");
238 throw kEH + "kinematics already converted";
240 AliStack* stack = fRunLoader->Stack();
242 throw kEH + "stack is null.";
245 fTreeK = new TTree("Kinematics", "Particles created during simulation.");
247 Int_t nentries = stack->GetNtrack();
248 std::vector<TEveMCTrack> vmc(nentries);
249 for (Int_t idx = 0; idx < nentries; ++idx)
251 TParticle* tp = stack->Particle(idx);
253 vmc[idx].fLabel = idx;
256 // read track refrences
257 // functionality now in AliEveKineTools.
259 TTree* fTreeTR = fRunLoader->TreeTR();
262 Warning(kEH, "no TrackRefs; some data will not be available.");
264 TClonesArray* RunArrayTR = 0;
265 fTreeTR->SetBranchAddress("AliRun", &RunArrayTR);
267 Int_t nPrimaries = (Int_t) fTreeTR->GetEntries();
268 for (Int_t iPrimPart = 0; iPrimPart<nPrimaries; iPrimPart++) {
269 // printf("T0 fTreeTR->GetEntry(%d) \n",iPrimPart);
270 fTreeTR->GetEntry(iPrimPart);
271 // printf("END fTreeTR->GetEntry(%d) \n",iPrimPart);
273 for (Int_t iTrackRef = 0; iTrackRef < RunArrayTR->GetEntriesFast(); iTrackRef++) {
274 AliTrackReference *trackRef = (AliTrackReference*)RunArrayTR->At(iTrackRef);
275 Int_t track = trackRef->GetTrack();
276 if(track < nentries && track > 0){
277 TEveMCTrack& mct = vmc[track];
278 if(trackRef->TestBit(kNotDeleted)) {
280 mct.t_decay = trackRef->GetTime();
281 mct.V_decay.x = trackRef->X();
282 mct.V_decay.y = trackRef->Y();
283 mct.V_decay.z = trackRef->Z();
284 mct.P_decay.x = trackRef->Px();
285 mct.P_decay.y = trackRef->Py();
286 mct.P_decay.z = trackRef->Pz();
287 if(TMath::Abs(mct.GetPdgCode()) == 11)
288 mct.decayed = false; // a bug in TreeTR
296 fTreeK->Branch("K", "TEveMCTrack", &fpK, fBuffSize);
298 printf("sizeofvmc = %d\n", (Int_t) vmc.size());
299 for(std::vector<TEveMCTrack>::iterator k = vmc.begin(); k != vmc.end(); ++k)
301 TEveMCTrack& mct = *k;
305 Int_t mi = mct.fLabel;
307 while (m->GetMother(0) != -1)
310 printf("cnt %d mi=%d, mo=%d\n", cnt, mi, m->GetMother(0));
312 mi = m->GetMother(0);
321 fTreeK->BuildIndex("fLabel");
325 //==============================================================================
327 //==============================================================================
329 void AliEveVSDCreator::ConvertHits()
332 // TPC hits are handled specially as they are compressed - only mayor
335 static const TEveException kEH("AliEveVSDCreator::ConvertHits ");
338 throw kEH + "hits already converted.";
341 fTreeH = new TTree("Hits", "Simulated energy depositions in detector.");
342 fTreeH->Branch("H", "TEveHit", &fpH, fBuffSize);
346 ConvertAnyHits("ITS", "AliITShit", 0, 0);
347 ConvertTPCHits( 1, fTPCHitRes*fTPCHitRes);
348 ConvertAnyHits("TRD", "AliTRDhit", 2, fTRDHitRes*fTRDHitRes);
349 ConvertAnyHits("TOF", "AliTOFhit", 3, 0);
352 catch(TEveException& exc)
354 Warning(kEH, "%s", exc.Data());
358 //------------------------------------------------------------------------------
360 void AliEveVSDCreator::ConvertAnyHits(const TString& detector,
361 const TString& clones_class,
362 Int_t det_id, Float_t minDistSqr)
364 // Convert hits for detector.
366 static const TEveException kEH("AliEveVSDCreator::ConvertAnyHits ");
368 Float_t x,y,z, x1,y1,z1;
370 TTree* treeh = fRunLoader->GetTreeH(detector, kFALSE);
373 Warning(kEH, "no hits for %s.", detector.Data());
377 TClonesArray *arr = new TClonesArray(clones_class);
378 treeh->SetBranchAddress(detector, &arr);
379 Int_t np = treeh->GetEntries();
381 // In TreeH files hits are grouped in clones arrays
382 // each eva particle has its own clone array.
384 std::map<Int_t, Int_t> hmap;
386 for (Int_t i = 0; i < np; ++i)
389 Int_t evaIdx = np - i - 1;
390 Int_t nh = arr->GetEntriesFast();
392 // printf("%d entry %d hits for primary %d \n", i, nh, evaIdx);
393 for (Int_t j = 0; j < nh; ++j)
395 AliHit* aliHit = (AliHit*)arr->UncheckedAt(j);
397 x1 = aliHit->X(); y1 = aliHit->Y(); z1 = aliHit->Z();
398 if (minDistSqr == 0 || (x-x1)*(x-x1) + (y-y1)*(y-y1) + (z-z1)*(z-z1) > minDistSqr)
402 fH.fLabel = aliHit->GetTrack();
403 fH.fEvaLabel = evaIdx;
404 fH.fV.Set(aliHit->X(), aliHit->Y(), aliHit->Z());
417 for (std::map<Int_t, Int_t>::iterator j = hmap.begin(); j != hmap.end(); ++j)
419 GetGeninfo(j->first)->fNHits += j->second;
423 //------------------------------------------------------------------------------
425 void AliEveVSDCreator::ConvertTPCHits(Int_t det_id, Float_t minDistSqr)
427 // Convert hits for TPC.
429 static const TEveException kEH("AliEveVSDCreator::ConvertTPCHits ");
431 Float_t x,y,z, x1,y1,z1;
433 TTree* treeh = fRunLoader->GetTreeH("TPC", false);
436 Warning(kEH, "no hits for %s.", "TPC");
440 AliTPCTrackHitsV2 hv2, *hv2p = &hv2;
441 treeh->SetBranchAddress("TPC2", &hv2p);
442 Int_t np = treeh->GetEntries();
444 std::map<Int_t, Int_t> hmap;
446 for (Int_t i = 0; i < np; ++i)
449 Int_t evaIdx = np -i -1;
450 if (hv2.First() == 0) continue;
454 AliHit* ah = hv2.GetHit();
455 x1 = ah->X(); y1 = ah->Y(); z1 = ah->Z();
456 if (minDistSqr == 0 || (x-x1)*(x-x1) + (y-y1)*(y-y1) + (z-z1)*(z-z1) > minDistSqr)
460 fH.fLabel = ah->Track();
461 fH.fEvaLabel = evaIdx;
462 fH.fV.fX = x1; fH.fV.fY = y1; fH.fV.fZ = z1;
468 x = x1; y = y1; z = z1;
470 } while (hv2.Next());
474 for (std::map<Int_t, Int_t>::iterator j = hmap.begin(); j != hmap.end(); ++j)
476 GetGeninfo(j->first)->fNHits += j->second;
481 //==============================================================================
483 //==============================================================================
485 void AliEveVSDCreator::ConvertClusters()
489 // Only supported for ITS and TPC at the moment, see dedicated
490 // functions ConvertITSClusters() and ConvertTPCClusters().
492 // It should be possible now to do this in a general manner (with
493 // the alignment framework).
495 static const TEveException kEH("AliEveVSDCreator::ConvertClusters ");
498 throw kEH + "clusters already converted.";
501 fTreeC = new TTree("Clusters", "Reconstructed points of particle passage.");
502 fTreeC->Branch("C", "TEveCluster", &fpC, fBuffSize);
506 ConvertAnyClusters("ITS", "ITSRecPoints", 0);
507 ConvertTPCClusters( 1);
508 ConvertAnyClusters("TRD", "TRDcluster", 2);
509 ConvertAnyClusters("TOF", "TOF", 3);
512 catch(TEveException& exc)
514 Warning(kEH, "%s", exc.Data());
518 //------------------------------------------------------------------------------
520 void AliEveVSDCreator::ConvertAnyClusters(const TString& detector,
521 const TString& branch_name,
524 // Convert clusters for detector and transform them to global coordinates.
526 static const TEveException kEH("AliEveVSDCreator::ConvertITSClusters ");
528 fRunLoader->LoadRecPoints(detector);
529 TTree* tree = fRunLoader->GetTreeR(detector, false);
531 throw kEH + "'TreeR' not found.";
533 TClonesArray *cl = 0;
534 TBranch *branch = tree->GetBranch(branch_name);
535 branch->SetAddress(&cl);
537 Int_t nmods = branch->GetEntries();
539 std::map<Int_t, Int_t> cmap;
541 for (Int_t mod = 0; mod < nmods; ++mod)
543 branch->GetEntry(mod);
544 Int_t nc = cl->GetEntriesFast();
545 for (Int_t j = 0; j < nc; ++j)
547 AliCluster *c = (AliCluster*) cl->UncheckedAt(j);
551 fC.fLabel[0] = c->GetLabel(0);
552 fC.fLabel[1] = c->GetLabel(1);
553 fC.fLabel[2] = c->GetLabel(2);
555 c->GetGlobalXYZ(fC.fV.Arr());
561 while (i < 3 && fC.fLabel[i])
562 cmap[fC.fLabel[i++]]++;
568 for (std::map<Int_t, Int_t>::iterator j=cmap.begin(); j!=cmap.end(); ++j)
570 GetGeninfo(j->first)->fNClus += j->second;
574 //------------------------------------------------------------------------------
576 void AliEveVSDCreator::ConvertTPCClusters(Int_t det_id)
578 // Convert TPC clusters and transform them to global coordinates.
580 static const TEveException kEH("AliEveVSDCreator::ConvertTPCClusters ");
582 const Int_t kMaxCl = 100*160;
584 fRunLoader->LoadRecPoints("TPC");
585 TTree* tree = fRunLoader->GetTreeR("TPC", false);
587 throw kEH + "'TreeR' not found.";
589 AliTPCClustersRow *clrow = new AliTPCClustersRow;
590 clrow->SetClass("AliTPCclusterMI");
591 clrow->SetArray(kMaxCl);
592 tree->SetBranchAddress("Segment", &clrow);
594 Int_t nEnt = tree->GetEntries();
596 std::map<Int_t, Int_t> cmap;
598 for (Int_t j = 0; j < nEnt; ++j)
600 if (!tree->GetEvent(j)) continue;
602 TClonesArray *cl = clrow->GetArray();
603 Int_t ncl = cl->GetEntriesFast();
607 AliCluster *c = (AliCluster*) cl->UncheckedAt(ncl);
611 fC.fLabel[0] = c->GetLabel(0);
612 fC.fLabel[1] = c->GetLabel(1);
613 fC.fLabel[2] = c->GetLabel(2);
615 c->GetGlobalXYZ(fC.fV.Arr());
621 while (i < 3 && fC.fLabel[i])
622 cmap[fC.fLabel[i++]]++;
631 for (std::map<Int_t, Int_t>::iterator j=cmap.begin(); j!=cmap.end(); ++j)
633 GetGeninfo(j->first)->fNClus += j->second;
638 /******************************************************************************/
640 /******************************************************************************/
642 void AliEveVSDCreator::ConvertRecTracks()
644 // Convert reconstructed tracks.
646 static const TEveException kEH("AliEveVSDCreator::ConvertRecTracks ");
649 throw kEH + "tracks already converted.";
651 AliESDEvent* esdEvent = AliEveEventManager::AssertESD();
654 fTreeR = new TTree("RecTracks", "Reconstructed particle trajectories.");
655 fTreeR->Branch("R", "TEveRecTrack", &fpR, 512*1024);
657 // reconstructed tracks
658 AliESDtrack* esdTrack;
660 for (Int_t n = 0; n < esdEvent->GetNumberOfTracks(); ++n)
662 esdTrack = esdEvent->GetTrack(n);
664 fR.fLabel = esdTrack->GetLabel();
665 fR.fStatus = (Int_t) esdTrack->GetStatus();
666 fR.fSign = (Int_t) esdTrack->GetSign();
667 esdTrack->GetXYZ(dbuf); fR.fV.Set(dbuf);
668 esdTrack->GetPxPyPz(dbuf); fR.fP.Set(dbuf);
669 Double_t ep = esdTrack->GetP();
670 fR.fBeta = ep/TMath::Sqrt(ep*ep + TMath::C()*TMath::C()*esdTrack->GetMass()*esdTrack->GetMass());
674 fTreeR->BuildIndex("fLabel");
677 /******************************************************************************/
679 void AliEveVSDCreator::ConvertV0()
681 // Convert reconstructed V0s.
683 static const TEveException kEH("AliEveVSDCreator::ConvertV0 ");
686 throw kEH + "AliEveV0 already converted.";
688 AliESDEvent* esdEvent = AliEveEventManager::AssertESD();
691 fTreeV0 = new TTree("AliEveV0", "AliEveV0 points");
693 fTreeV0->Branch("AliEveV0", "TEveRecV0", &fpV0, 512*1024,1);
695 for (Int_t n = 0; n < esdEvent->GetNumberOfV0s(); ++n)
697 AliESDv0 *av = esdEvent->GetV0(n);
698 AliESDtrack *trackN = esdEvent->GetTrack(av->GetNindex()); // negative daughter
699 AliESDtrack *trackP = esdEvent->GetTrack(av->GetPindex()); // positive daughter
703 fV0.fStatus = av->GetStatus();
704 // Point of closest approach
705 av->GetXYZ(pos[0],pos[1],pos[2]);
706 fV0.fVCa.fX = pos[0];
707 fV0.fVCa.fY = pos[1];
708 fV0.fVCa.fZ = pos[2];
709 // set birth vertex of neutral particle
710 av->GetXYZ(pos[0], pos[1], pos[2]);
711 fV0.fV0Birth.Set(pos);
713 // momentum and position of negative particle
714 av->GetParamN()->GetPxPyPz(pos);
716 av->GetParamN()->GetXYZ(pos);
719 // momentum and position of positive particle
720 av->GetParamP()->GetPxPyPz(pos);
722 av->GetParamP()->GetXYZ(pos);
725 fV0.fLabel = 0; // !!!! mother label unknown
726 fV0.fPdg = av->GetPdgCode();
729 fV0.fDLabel[0] = TMath::Abs(trackN->GetLabel());
730 fV0.fDLabel[1] = TMath::Abs(trackP->GetLabel());
732 // printf("AliEveV0 convert labels(%d,%d) index(%d,%d)\n",
733 // fV0.d_label[0], fV0.d_label[1],
734 // av->GetNIndex(), av->GetPIndex());
738 // if (esdEvent->GetNumberOfV0s()) fTreeV0->BuildIndex("label");
742 /******************************************************************************/
744 void AliEveVSDCreator::ConvertKinks()
746 // Convert reconstructed kinks.
748 static const TEveException kEH("AliEveVSDCreator::ConvertKinks ");
751 throw kEH + "Kinks already converted.";
753 throw kEH + "Currently non-supported - TEveRecKink being updated.";
756 AliESDEvent* esdEvent = AliEveEventManager::AssertESD();
759 fTreeKK = new TTree("Kinks", "ESD Kinks");
761 fTreeKK->Branch("KK", "TEveRecKink", &fpKK, fBuffSize);
763 // printf("CONVERT KINK Read %d entries in tree kinks \n", esdEvent->GetNumberOfKinks());
764 for (Int_t n = 0; n < esdEvent->GetNumberOfKinks(); ++n)
766 AliESDkink* kk = esdEvent->GetKink(n);
770 fKK.fLabel = kk->GetLabel(0);
771 fKK.fStatus = 0; // status is Char_t[12] ... have no idea how/what to extract.
773 // reconstructed kink position
774 fKK.fLabelSec = kk->GetLabel(1);
775 fKK.fVKink.Set(kk->GetPosition());
777 const AliExternalTrackParam& tpMother = kk->RefParamMother();
778 // momentum and position of mother
779 tpMother.GetPxPyPz(pos);
781 tpMother.GetXYZ(pos);
783 const Double_t* par = tpMother.GetParameter();
784 // printf("KINK Pt %f, %f \n",1/tpMother.Pt(),par[4] );
785 fKK.fSign = (par[4] < 0) ? -1 : 1;
787 const AliExternalTrackParam& tpDaughter = kk->RefParamDaughter();
788 // momentum and position of daughter
789 tpDaughter.GetPxPyPz(pos);
791 tpDaughter.GetXYZ(pos);
796 if (esdEvent->GetNumberOfKinks()) fTreeKK->BuildIndex("label");
800 /******************************************************************************/
802 /******************************************************************************/
804 void AliEveVSDCreator::ConvertGenInfo()
806 // Build simulation-reconstruction cross-reference table.
807 // In a rather poor state at the moment.
809 static const TEveException kEH("AliEveVSDCreator::ConvertGenInfo ");
812 throw kEH + "GI already converted.";
815 fTreeGI = new TTree("TEveMCRecCrossRef", "Objects prepared for cross querry");
817 TEveMCRecCrossRef::Class()->IgnoreTObjectStreamer(true);
818 fTreeGI->Branch("GI", "TEveMCRecCrossRef", &fpGI, fBuffSize);
819 fTreeGI->Branch("K.", "TEveMCTrack", &fpK);
820 fTreeGI->Branch("R.", "TEveRecTrack", &fpR);
822 for (std::map<Int_t, TEveMCRecCrossRef*>::iterator j=fGenInfoMap.begin(); j!=fGenInfoMap.end(); ++j)
825 fGI.fLabel = j->first;
826 fTreeK->GetEntry(j->first);
829 Int_t re = fTreeR->GetEntryNumberWithIndex(j->first);
833 // Int_t hasV0 = fTreeV0->GetEntryNumberWithIndex(j->first);
835 // fGI.has_AliEveV0 = true;
838 Int_t hasKk = fTreeKK->GetEntryNumberWithIndex(j->first);
847 /******************************************************************************/
848 /******************************************************************************/
850 /******************************************************************************/
851 /******************************************************************************/
853 TEveMCRecCrossRef* AliEveVSDCreator::GetGeninfo(Int_t label)
855 // Return the cross-reference structure for given label.
856 // If the object does not exist it is created.
858 // printf("get_geninfo %d\n", label);
859 TEveMCRecCrossRef* gi;
860 std::map<Int_t, TEveMCRecCrossRef*>::iterator i = fGenInfoMap.find(label);
861 if (i == fGenInfoMap.end())
863 gi = new TEveMCRecCrossRef();
864 fGenInfoMap[label] = gi;