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>
35 //______________________________________________________________________________
37 // Create VSD file from ALICE data.
39 ClassImp(AliEveVSDCreator)
41 AliEveVSDCreator::AliEveVSDCreator(const Text_t* name, const Text_t* title) :
53 // Particles not in ROOT's PDG database occuring in ALICE
54 AliPDG::AddParticlesToPdgDataBase();
56 TDatabasePDG *pdgDB = TDatabasePDG::Instance();
57 // const Int_t kspe=50000000;
58 const Int_t kion=10000000;
60 const Double_t kAu2Gev=0.9314943228;
61 const Double_t khSlash = 1.0545726663e-27;
62 const Double_t kErg2Gev = 1/1.6021773349e-3;
63 const Double_t khShGev = khSlash*kErg2Gev;
64 const Double_t kYear2Sec = 3600*24*365.25;
66 pdgDB->AddParticle("Deuteron","Deuteron",2*kAu2Gev+8.071e-3,kTRUE,
67 0,1,"Ion",kion+10020);
68 pdgDB->AddParticle("Triton","Triton",3*kAu2Gev+14.931e-3,kFALSE,
69 khShGev/(12.33*kYear2Sec),1,"Ion",kion+10030);
70 pdgDB->AddParticle("Alpha","Alpha",4*kAu2Gev+2.424e-3,kTRUE,
71 khShGev/(12.33*kYear2Sec),2,"Ion",kion+20040);
72 pdgDB->AddParticle("HE3","HE3",3*kAu2Gev+14.931e-3,kFALSE,
73 0,2,"Ion",kion+20030);
76 // AliKalmanTrack::SetConvConst(1);
79 /******************************************************************************/
81 void AliEveVSDCreator::CreateVSD(const Text_t* vsdFile)
83 // Create the VSD for current event in AliEveEventManager.
84 // Result is stored in vsdFile.
86 // Needs to be extended to support conversion of multiple events.
88 static const TEveException kEH("AliEveVSDCreator::CreateVSD ");
90 fRunLoader = AliEveEventManager::AssertRunLoader();
93 printf("%s open seems ok. Now loading sim data.\n", kEH.Data());
95 fRunLoader->LoadHeader();
96 fRunLoader->LoadKinematics();
97 fRunLoader->LoadTrackRefs();
98 fRunLoader->LoadHits();
103 printf("%s opening output TEveVSD.\n", kEH.Data());
105 TFile* file = TFile::Open(vsdFile, "RECREATE", "ALICE Visualization Summary Data");
106 fDirectory = new TDirectoryFile("Event0", "");
109 printf("%s creating trees now ...\n", kEH.Data());
114 printf("%s trees created, closing files.\n", kEH.Data());
123 // clean after the TEveVSD data was sucessfuly written
136 printf("%s all done.\n", kEH.Data());
139 void AliEveVSDCreator::CreateTrees()
141 // Create and fill the output trees by calling all the
142 // ConvertXyzz() functions.
143 // Exceptions from individual functions are displayed as warnings.
145 static const TEveException kEH("AliEveVSDCreator::CreateTrees ");
148 throw(kEH + "output directory not set.");
152 printf("%sConvertKinematics.\n", kEH.Data());
154 } catch(TEveException& exc) { Warning(kEH, exc); }
156 Warning(kEH, "Explicitly abandoning further conversion.");
161 printf("%sConvertHits.\n", kEH.Data());
163 } catch(TEveException& exc) { Warning(kEH, exc); }
167 printf("%sConvertClusters.\n", kEH.Data());
169 } catch(TEveException& exc) { Warning(kEH, exc); }
173 printf("%sConvertRecTracks.\n", kEH.Data());
175 } catch(TEveException& exc) {
176 Warning(exc, "skipping AliEveV0 extraction.");
177 goto end_esd_processing;
182 printf("%sConvertV0.\n", kEH.Data());
184 } catch(TEveException& exc) { Warning(kEH, exc); }
188 printf("%sConvertKinks.\n", kEH.Data());
190 } catch(TEveException& exc) { Warning(kEH, exc); }
196 printf("%sConvertGenInfo.\n", kEH.Data());
198 } catch(TEveException& exc) { Warning(kEH, exc); }
203 /******************************************************************************/
205 /******************************************************************************/
207 void AliEveVSDCreator::ConvertKinematics()
209 // Convert kinematics.
210 // Track references are not stored, they should be.
212 static const TEveException kEH("AliEveVSDCreator::ConvertKinematics ");
215 throw (kEH + "kinematics already converted");
217 AliStack* stack = fRunLoader->Stack();
219 throw(kEH + "stack is null.");
222 fTreeK = new TTree("Kinematics", "TParticles sorted by Label");
224 Int_t nentries = stack->GetNtrack();
225 std::vector<TEveMCTrack> vmc(nentries);
226 for (Int_t idx=0; idx<nentries; idx++) {
227 TParticle* tp = stack->Particle(idx);
229 vmc[idx].fLabel = idx;
232 // read track refrences
233 // functionality now in AliEveKineTools.
235 TTree* fTreeTR = fRunLoader->TreeTR();
238 Warning(kEH, "no TrackRefs; some data will not be available.");
240 TClonesArray* RunArrayTR = 0;
241 fTreeTR->SetBranchAddress("AliRun", &RunArrayTR);
243 Int_t nPrimaries = (Int_t) fTreeTR->GetEntries();
244 for (Int_t iPrimPart = 0; iPrimPart<nPrimaries; iPrimPart++) {
245 // printf("T0 fTreeTR->GetEntry(%d) \n",iPrimPart);
246 fTreeTR->GetEntry(iPrimPart);
247 // printf("END fTreeTR->GetEntry(%d) \n",iPrimPart);
249 for (Int_t iTrackRef = 0; iTrackRef < RunArrayTR->GetEntriesFast(); iTrackRef++) {
250 AliTrackReference *trackRef = (AliTrackReference*)RunArrayTR->At(iTrackRef);
251 Int_t track = trackRef->GetTrack();
252 if(track < nentries && track > 0){
253 TEveMCTrack& mct = vmc[track];
254 if(trackRef->TestBit(kNotDeleted)) {
256 mct.t_decay = trackRef->GetTime();
257 mct.V_decay.x = trackRef->X();
258 mct.V_decay.y = trackRef->Y();
259 mct.V_decay.z = trackRef->Z();
260 mct.P_decay.x = trackRef->Px();
261 mct.P_decay.y = trackRef->Py();
262 mct.P_decay.z = trackRef->Pz();
263 if(TMath::Abs(mct.GetPdgCode()) == 11)
264 mct.decayed = false; // a bug in TreeTR
272 fTreeK->Branch("K", "TEveMCTrack", &fpK, fBuffSize);
274 printf("sizeofvmc = %d\n", (Int_t) vmc.size());
275 for(std::vector<TEveMCTrack>::iterator k=vmc.begin(); k!=vmc.end(); ++k) {
276 TEveMCTrack& mct = *k;
280 Int_t mi = mct.fLabel;
282 while(m->GetMother(0) != -1) {
284 printf("cnt %d mi=%d, mo=%d\n", cnt, mi, m->GetMother(0));
286 mi = m->GetMother(0);
295 fTreeK->BuildIndex("fLabel");
298 /******************************************************************************/
300 /******************************************************************************/
306 const char* fName; // Detector name.
307 const char* fHitbranch; // Name of branch containing hits.
308 unsigned char fDetidx; // Index identifying the detector internally.
311 Detector_t fgDetectors[] = {
312 { "ITS", "AliITShit", 0 },
313 { "TPC", "AliTPCTrackHitsV2", 1 },
314 { "TRD", "AliTRDhit", 2 },
315 { "TOF", "AliTOFhit", 3 }
316 // { "HMPID", "AliHMPIDhit", 4 },
321 /******************************************************************************/
323 void AliEveVSDCreator::ConvertHits()
326 // TPC hits are handled specially as they are compressed - only mayor
329 static const TEveException kEH("AliEveVSDCreator::ConvertHits ");
332 throw(kEH + "hits already converted.");
335 fTreeH = new TTree("Hits", "Combined detector hits.");
336 fTreeH->Branch("H", "TEveHit", &fpH, fBuffSize);
338 std::map<Int_t, Int_t> hmap;
339 // parameters for ITS, TPC hits filtering
340 Float_t x,y,z, x1,y1,z1;
341 Float_t tpcSqrRes = fTPCHitRes*fTPCHitRes;
342 Float_t trdSqrRes = fTRDHitRes*fTRDHitRes;
345 // load hits from the rest of detectors
346 while (fgDetectors[l].fName != 0)
348 Detector_t& det = fgDetectors[l++];
355 TTree* treeh = fRunLoader->GetTreeH(det.fName, false);
357 Warning(kEH, Form("no hits for %s.", det.fName));
360 AliTPCTrackHitsV2 hv2, *hv2p = &hv2;
361 treeh->SetBranchAddress("TPC2", &hv2p);
362 Int_t np = treeh->GetEntries();
363 for (Int_t i = 0; i < np; ++i)
366 Int_t evaIdx = np -i -1;
367 if (hv2.First() == 0) continue;
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) > tpcSqrRes)
374 fH.fDetId = det.fDetidx;
376 fH.fLabel = ah->Track();
377 fH.fEvaLabel = evaIdx;
378 fH.fV.fX = x1; fH.fV.fY = y1; fH.fV.fZ = z1;
381 x = x1; y = y1; z = z1;
384 } while (hv2.Next());
386 // printf("%d entries in TPChits \n",count);
391 TTree* treeh = fRunLoader->GetTreeH(det.fName, false);
393 Warning(kEH, Form("no hits for %s.", det.fName));
396 TClonesArray *arr = new TClonesArray(det.fHitbranch);
397 treeh->SetBranchAddress(det.fName, &arr);
398 Int_t np = treeh->GetEntries();
399 // in TreeH files hits are grouped in clones arrays
400 // each eva particle has its own clone array
401 for (Int_t i = 0; i < np; ++i)
404 Int_t evaIdx = np -i -1;
405 Int_t nh=arr->GetEntriesFast();
407 // printf("%d entry %d hits for primary %d \n", i, nh, evaIdx);
408 for (Int_t j = 0; j < nh; ++j)
410 AliHit* aliHit = (AliHit*)arr->UncheckedAt(j);
411 fH.fDetId = det.fDetidx;
413 fH.fLabel = aliHit->GetTrack();
414 fH.fEvaLabel = evaIdx;
415 fH.fV.Set(aliHit->X(), aliHit->Y(), aliHit->Z());
416 if (det.fDetidx == 2)
418 x1=aliHit->X();y1=aliHit->Y();z1=aliHit->Z();
419 if ((x-x1)*(x-x1)+(y-y1)*(y-y1)+(z-z1)*(z-z1) < trdSqrRes) continue;
434 for(std::map<Int_t, Int_t>::iterator j=hmap.begin(); j!=hmap.end(); ++j)
436 GetGeninfo(j->first)->fNHits += j->second;
440 /******************************************************************************/
442 /******************************************************************************/
444 void AliEveVSDCreator::ConvertClusters()
448 // Only supported for ITS and TPC at the moment, see dedicated
449 // functions ConvertITSClusters() and ConvertTPCClusters().
451 // It should be possible now to do this in a general manner (with
452 // the alignment framework).
454 static const TEveException kEH("AliEveVSDCreator::ConvertClusters ");
457 throw(kEH + "clusters already converted.");
460 fTreeC = new TTree("Clusters", "rec clusters");
461 fTreeC->Branch("C", "TEveCluster", &fpC, fBuffSize);
464 ConvertITSClusters();
465 } catch(TEveException& exc) { Warning(kEH, exc); }
468 ConvertTPCClusters();
469 } catch(TEveException& exc) { Warning(kEH, exc); }
472 /******************************************************************************/
474 void AliEveVSDCreator::ConvertTPCClusters()
476 // Convert TPC clusters and transform them to global coordinates.
478 static const TEveException kEH("AliEveVSDCreator::ConvertTPCClusters ");
480 fRunLoader->LoadRecPoints("TPC");
481 TTree* tree = fRunLoader->GetTreeR("TPC", false);
483 throw(kEH + "'TreeR' not found.");
485 AliTPCClustersRow clrow, *clrowp = &clrow;
487 clrow.SetClass("AliTPCclusterMI");
488 tree->SetBranchAddress("Segment", &clrowp);
492 Int_t nEnt = tree->GetEntries();
493 for (Int_t n = 0; n < nEnt; ++n)
496 nClusters += clrow.GetArray()->GetEntriesFast();
499 std::map<Int_t, Int_t> cmap;
501 for (Int_t n = 0; n < tree->GetEntries(); ++n)
504 Int_t ncl = clrow.GetArray()->GetEntriesFast();
509 if (clrow.GetArray())
511 // cl = new AliTPCclusterMI(*(AliTPCclusterMI*)clrow.GetArray()->UncheckedAt(ncl));
512 cl = (AliTPCclusterMI*)clrow.GetArray()->UncheckedAt(ncl);
513 if (cl->GetLabel(0) >= 0)
515 Float_t g[3]; //global coordinates
519 fC.fSubdetId = clrow.GetID();
520 fC.fLabel[0] = cl->GetLabel(0);
521 fC.fLabel[1] = cl->GetLabel(1);
522 fC.fLabel[2] = cl->GetLabel(2);
528 while (i < 3 && fC.fLabel[i])
529 cmap[fC.fLabel[i++]]++;
537 for (std::map<Int_t, Int_t>::iterator j=cmap.begin(); j!=cmap.end(); ++j)
539 GetGeninfo(j->first)->fNClus += j->second;
543 /******************************************************************************/
545 void AliEveVSDCreator::ConvertITSClusters()
547 // Convert ITS clusters and transform them to global coordinates.
549 static const TEveException kEH("AliEveVSDCreator::ConvertITSClusters ");
551 fRunLoader->LoadRecPoints("ITS");
552 TTree* tree = fRunLoader->GetTreeR("ITS", false);
554 throw(kEH + "'TreeR' not found.");
557 AliITSLoader *itsLd = (AliITSLoader*) fRunLoader->GetLoader("ITSLoader");
558 AliITSgeom *geom = itsLd->GetITSgeom();
560 //printf("alice ITS geom %p \n",geom );
563 throw(kEH + "can not find ITS geometry");
565 TClonesArray *arr = new TClonesArray("AliITSclusterV2");
566 tree->SetBranchAddress("Clusters", &arr);
567 Int_t nmods = tree->GetEntries();
569 std::map<Int_t, Int_t> cmap;
571 for (Int_t mod = 0; mod < nmods; ++mod)
574 Int_t nc=arr->GetEntriesFast();
575 for (Int_t j = 0; j < nc; ++j)
577 AliITSclusterV2* recp = (AliITSclusterV2*) arr->UncheckedAt(j);
580 geom->GetRotMatrix(mod,rot);
582 geom->GetModuleId(mod,lay,lad,det);
584 geom->GetTrans(lay,lad,det,tx,ty,tz);
586 Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
587 Double_t phi1=TMath::Pi()/2+alpha;
588 if (lay == 1) phi1+=TMath::Pi();
590 Float_t cp=TMath::Cos(phi1), sp=TMath::Sin(phi1);
591 Float_t r=tx*cp+ty*sp;
592 gc[0] = r*cp - recp->GetY()*sp;
593 gc[1] = r*sp + recp->GetY()*cp;
594 gc[2] = recp->GetZ();
598 fC.fLabel[0] = recp->GetLabel(0);
599 fC.fLabel[1] = recp->GetLabel(1);
600 fC.fLabel[2] = recp->GetLabel(2);
601 fC.fV.fX = r*cp - recp->GetY()*sp;
602 fC.fV.fY = r*sp + recp->GetY()*cp;
603 fC.fV.fZ = recp->GetZ();
606 while(i < 3 && fC.fLabel[i])
607 cmap[fC.fLabel[i++]]++;
611 for (std::map<Int_t, Int_t>::iterator j=cmap.begin(); j!=cmap.end(); ++j)
613 GetGeninfo(j->first)->fNClus += j->second;
619 /******************************************************************************/
621 /******************************************************************************/
623 void AliEveVSDCreator::ConvertRecTracks()
625 // Convert reconstructed tracks.
627 static const TEveException kEH("AliEveVSDCreator::ConvertRecTracks ");
630 throw(kEH + "tracks already converted.");
632 AliESDEvent* esdEvent = AliEveEventManager::AssertESD();
635 fTreeR = new TTree("RecTracks", "rec tracks");
637 fTreeR->Branch("R", "TEveRecTrack", &fpR, 512*1024,1);
639 // reconstructed tracks
640 AliESDtrack* esdTrack;
642 for (Int_t n = 0; n < esdEvent->GetNumberOfTracks(); ++n)
644 esdTrack = esdEvent->GetTrack(n);
646 fR.fLabel = esdTrack->GetLabel();
647 fR.fStatus = (Int_t) esdTrack->GetStatus();
648 fR.fSign = (Int_t) esdTrack->GetSign();
649 esdTrack->GetXYZ(dbuf); fR.fV.Set(dbuf);
650 esdTrack->GetPxPyPz(dbuf); fR.fP.Set(dbuf);
651 Double_t ep = esdTrack->GetP();
652 fR.fBeta = ep/TMath::Sqrt(ep*ep + TMath::C()*TMath::C()*esdTrack->GetMass()*esdTrack->GetMass());
655 fTreeR->BuildIndex("label");
659 /******************************************************************************/
661 void AliEveVSDCreator::ConvertV0()
663 // Convert reconstructed V0s.
665 static const TEveException kEH("AliEveVSDCreator::ConvertV0 ");
668 throw(kEH + "AliEveV0 already converted.");
670 AliESDEvent* esdEvent = AliEveEventManager::AssertESD();
673 fTreeV0 = new TTree("AliEveV0", "AliEveV0 points");
675 fTreeV0->Branch("AliEveV0", "TEveRecV0", &fpV0, 512*1024,1);
677 for (Int_t n = 0; n < esdEvent->GetNumberOfV0s(); ++n)
679 AliESDv0 *av = esdEvent->GetV0(n);
680 AliESDtrack *trackN = esdEvent->GetTrack(av->GetNindex()); // negative daughter
681 AliESDtrack *trackP = esdEvent->GetTrack(av->GetPindex()); // positive daughter
685 fV0.fStatus = av->GetStatus();
686 // Point of closest approach
687 av->GetXYZ(pos[0],pos[1],pos[2]);
688 fV0.fVCa.fX = pos[0];
689 fV0.fVCa.fY = pos[1];
690 fV0.fVCa.fZ = pos[2];
691 // set birth vertex of neutral particle
692 av->GetXYZ(pos[0], pos[1], pos[2]);
693 fV0.fV0Birth.Set(pos);
695 // momentum and position of negative particle
696 av->GetParamN()->GetPxPyPz(pos);
698 av->GetParamN()->GetXYZ(pos);
701 // momentum and position of positive particle
702 av->GetParamP()->GetPxPyPz(pos);
704 av->GetParamP()->GetXYZ(pos);
707 fV0.fLabel = 0; // !!!! mother label unknown
708 fV0.fPdg = av->GetPdgCode();
711 fV0.fDLabel[0] = TMath::Abs(trackN->GetLabel());
712 fV0.fDLabel[1] = TMath::Abs(trackP->GetLabel());
714 // printf("AliEveV0 convert labels(%d,%d) index(%d,%d)\n",
715 // fV0.d_label[0], fV0.d_label[1],
716 // av->GetNIndex(), av->GetPIndex());
720 // if (esdEvent->GetNumberOfV0s()) fTreeV0->BuildIndex("label");
724 /******************************************************************************/
726 void AliEveVSDCreator::ConvertKinks()
728 // Convert reconstructed kinks.
730 static const TEveException kEH("AliEveVSDCreator::ConvertKinks ");
733 throw(kEH + "Kinks already converted.");
735 throw kEH + "Currently non-supported - TEveRecKink being updated.";
738 AliESDEvent* esdEvent = AliEveEventManager::AssertESD();
741 fTreeKK = new TTree("Kinks", "ESD Kinks");
743 fTreeKK->Branch("KK", "TEveRecKink", &fpKK, fBuffSize);
745 // printf("CONVERT KINK Read %d entries in tree kinks \n", esdEvent->GetNumberOfKinks());
746 for (Int_t n = 0; n < esdEvent->GetNumberOfKinks(); ++n)
748 AliESDkink* kk = esdEvent->GetKink(n);
752 fKK.fLabel = kk->GetLabel(0);
753 fKK.fStatus = 0; // status is Char_t[12] ... have no idea how/what to extract.
755 // reconstructed kink position
756 fKK.fLabelSec = kk->GetLabel(1);
757 fKK.fVKink.Set(kk->GetPosition());
759 const AliExternalTrackParam& tpMother = kk->RefParamMother();
760 // momentum and position of mother
761 tpMother.GetPxPyPz(pos);
763 tpMother.GetXYZ(pos);
765 const Double_t* par = tpMother.GetParameter();
766 // printf("KINK Pt %f, %f \n",1/tpMother.Pt(),par[4] );
767 fKK.fSign = (par[4] < 0) ? -1 : 1;
769 const AliExternalTrackParam& tpDaughter = kk->RefParamDaughter();
770 // momentum and position of daughter
771 tpDaughter.GetPxPyPz(pos);
773 tpDaughter.GetXYZ(pos);
778 if (esdEvent->GetNumberOfKinks()) fTreeKK->BuildIndex("label");
782 /******************************************************************************/
784 /******************************************************************************/
786 void AliEveVSDCreator::ConvertGenInfo()
788 // Build simulation-reconstruction cross-reference table.
789 // In a rather poor state at the moment.
791 static const TEveException kEH("AliEveVSDCreator::ConvertGenInfo ");
794 throw(kEH + "GI already converted.");
797 fTreeGI = new TTree("TEveMCRecCrossRef", "Objects prepared for cross querry");
799 TEveMCRecCrossRef::Class()->IgnoreTObjectStreamer(true);
800 fTreeGI->Branch("GI", "TEveMCRecCrossRef", &fpGI, fBuffSize);
801 fTreeGI->Branch("K.", "TEveMCTrack", &fpK);
802 fTreeGI->Branch("R.", "TEveRecTrack", &fpR);
804 for (std::map<Int_t, TEveMCRecCrossRef*>::iterator j=fGenInfoMap.begin(); j!=fGenInfoMap.end(); ++j) {
806 fGI.fLabel = j->first;
807 fTreeK->GetEntry(j->first);
810 Int_t re = fTreeR->GetEntryNumberWithIndex(j->first);
814 // Int_t hasV0 = fTreeV0->GetEntryNumberWithIndex(j->first);
816 // fGI.has_AliEveV0 = true;
818 Int_t hasKk = fTreeKK->GetEntryNumberWithIndex(j->first);
827 /******************************************************************************/
828 /******************************************************************************/
830 /******************************************************************************/
831 /******************************************************************************/
833 TEveMCRecCrossRef* AliEveVSDCreator::GetGeninfo(Int_t label)
835 // Return the cross-reference structure for given label.
836 // If the object does not exist it is created.
838 // printf("get_geninfo %d\n", label);
839 TEveMCRecCrossRef* gi;
840 std::map<Int_t, TEveMCRecCrossRef*>::iterator i = fGenInfoMap.find(label);
841 if (i == fGenInfoMap.end()) {
842 gi = new TEveMCRecCrossRef();
843 fGenInfoMap[label] = gi;