Fork me on GitHub

source: git/display/DelphesEventDisplay.cc@ 30bb83a

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 30bb83a was 30bb83a, checked in by Christophe Delaere <christophe.delaere@…>, 10 years ago

Some small cleaning
  • Property mode set to 100644
File size: 17.8 KB
Line 
1#include <cassert>
2#include <iostream>
3#include <utility>
4#include <algorithm>
5#include "TGeoManager.h"
6#include "TGeoVolume.h"
7#include "external/ExRootAnalysis/ExRootConfReader.h"
8#include "external/ExRootAnalysis/ExRootTreeReader.h"
9#include "display/DelphesCaloData.h"
10#include "display/DelphesBranchElement.h"
11#include "display/Delphes3DGeometry.h"
12#include "display/DelphesEventDisplay.h"
13#include "classes/DelphesClasses.h"
14#include "TEveElement.h"
15#include "TEveJetCone.h"
16#include "TEveTrack.h"
17#include "TEveTrackPropagator.h"
18#include "TEveCalo.h"
19#include "TEveManager.h"
20#include "TEveGeoNode.h"
21#include "TEveTrans.h"
22#include "TEveViewer.h"
23#include "TEveBrowser.h"
24#include "TEveArrow.h"
25#include "TMath.h"
26#include "TSystem.h"
27#include "TRootBrowser.h"
28#include "TGButton.h"
29#include "TClonesArray.h"
30
31DelphesEventDisplay::DelphesEventDisplay()
32{
33 event_id_ = 0;
34 tkRadius_ = 1.29;
35 totRadius_ = 2.0;
36 tkHalfLength_ = 3.0;
37 bz_ = 3.8;
38 chain_ = new TChain("Delphes");
39 treeReader_ = 0;
40 delphesDisplay_ = 0;
41}
42
43DelphesEventDisplay::~DelphesEventDisplay()
44{
45 delete chain_;
46}
47
48DelphesEventDisplay::DelphesEventDisplay(const char *configFile, const char *inputFile, Delphes3DGeometry& det3D)
49{
50 event_id_ = 0;
51 tkRadius_ = 1.29;
52 totRadius_ = 2.0;
53 tkHalfLength_ = 3.0;
54 bz_ = 3.8;
55 chain_ = new TChain("Delphes");
56 treeReader_ = 0;
57 delphesDisplay_ = 0;
58
59 // initialize the application
60 TEveManager::Create(kTRUE, "IV");
61 TGeoManager* geom = gGeoManager;
62
63 // build the detector
64 tkRadius_ = det3D.getTrackerRadius();
65 totRadius_ = det3D.getDetectorRadius();
66 tkHalfLength_ = det3D.getTrackerHalfLength();
67 bz_ = det3D.getBField();
68 TGeoVolume* top = det3D.getDetector(false);
69 geom->SetTopVolume(top);
70 TEveElementList *geometry = new TEveElementList("Geometry");
71 TObjArray* nodes = top->GetNodes();
72 TIter itNodes(nodes);
73 TGeoNode* nodeobj;
74 TEveGeoTopNode* node;
75 while((nodeobj = (TGeoNode*)itNodes.Next())) {
76 node = new TEveGeoTopNode(gGeoManager,nodeobj);
77 node->UseNodeTrans();
78 geometry->AddElement(node);
79 }
80
81 // Create chain of root trees
82 chain_->Add(inputFile);
83
84 // Create object of class ExRootTreeReader
85 printf("*** Opening Delphes data file ***\n");
86 treeReader_ = new ExRootTreeReader(chain_);
87
88 // prepare data collections
89 readConfig(configFile, det3D, elements_, arrays_);
90 for(std::vector<DelphesBranchBase*>::iterator element = elements_.begin(); element<elements_.end(); ++element) {
91 DelphesBranchElement<TEveTrackList>* item_v1 = dynamic_cast<DelphesBranchElement<TEveTrackList>*>(*element);
92 DelphesBranchElement<TEveElementList>* item_v2 = dynamic_cast<DelphesBranchElement<TEveElementList>*>(*element);
93 if(item_v1) gEve->AddElement(item_v1->GetContainer());
94 if(item_v2) gEve->AddElement(item_v2->GetContainer());
95 }
96
97 // viewers and scenes
98 delphesDisplay_ = new DelphesDisplay;
99 gEve->AddGlobalElement(geometry);
100 delphesDisplay_->ImportGeomRPhi(geometry);
101 delphesDisplay_->ImportGeomRhoZ(geometry);
102 // find the first calo data and use that to initialize the calo display
103 for(std::vector<DelphesBranchBase*>::iterator data=elements_.begin();data<elements_.end();++data) {
104 if(TString((*data)->GetType())=="tower") { // we could also use GetClassName()=="DelphesCaloData"
105 DelphesCaloData* container = dynamic_cast<DelphesBranchElement<DelphesCaloData>*>((*data))->GetContainer();
106 assert(container);
107 TEveCalo3D *calo3d = new TEveCalo3D(container);
108 calo3d->SetBarrelRadius(tkRadius_);
109 calo3d->SetEndCapPos(tkHalfLength_);
110 gEve->AddGlobalElement(calo3d);
111 delphesDisplay_->ImportCaloRPhi(calo3d);
112 delphesDisplay_->ImportCaloRhoZ(calo3d);
113 TEveCaloLego *lego = new TEveCaloLego(container);
114 lego->InitMainTrans();
115 lego->RefMainTrans().SetScale(TMath::TwoPi(), TMath::TwoPi(), TMath::Pi());
116 lego->SetAutoRebin(kFALSE);
117 lego->Set2DMode(TEveCaloLego::kValSizeOutline);
118 delphesDisplay_->ImportCaloLego(lego);
119 break;
120 }
121 }
122
123 make_gui();
124 load_event();
125 gEve->Redraw3D(kTRUE);
126
127}
128
129// function that parses the config to extract the branches of interest and prepare containers
130void DelphesEventDisplay::readConfig(const char *configFile, Delphes3DGeometry& det3D, std::vector<DelphesBranchBase*>& elements, std::vector<TClonesArray*>& arrays) {
131 ExRootConfReader *confReader = new ExRootConfReader;
132 confReader->ReadFile(configFile);
133 Double_t tk_radius = det3D.getTrackerRadius();
134 Double_t tk_length = det3D.getTrackerHalfLength();
135 Double_t tk_Bz = det3D.getBField();
136 Double_t mu_radius = det3D.getDetectorRadius();
137 Double_t mu_length = det3D.getDetectorHalfLength();
138 TAxis* etaAxis = det3D.getCaloAxes().first;
139 TAxis* phiAxis = det3D.getCaloAxes().second;
140 ExRootConfParam branches = confReader->GetParam("TreeWriter::Branch");
141 Int_t nBranches = branches.GetSize()/3;
142 DelphesBranchElement<TEveTrackList>* tlist;
143 DelphesBranchElement<DelphesCaloData>* clist;
144 DelphesBranchElement<TEveElementList>* elist;
145 // first loop with all but tracks
146 for(Int_t b = 0; b<nBranches; ++b) {
147 TString input = branches[b*3].GetString();
148 TString name = branches[b*3+1].GetString();
149 TString className = branches[b*3+2].GetString();
150 if(className=="Tower") {
151 if(input.Contains("eflow",TString::kIgnoreCase) || name.Contains("eflow",TString::kIgnoreCase)) continue; //no eflow
152 clist = new DelphesBranchElement<DelphesCaloData>(name,"tower",kBlack);
153 clist->GetContainer()->SetEtaBins(etaAxis);
154 clist->GetContainer()->SetPhiBins(phiAxis);
155 elements.push_back(clist);
156 } else if(className=="Jet") {
157 if(input.Contains("GenJetFinder")) {
158 elist = new DelphesBranchElement<TEveElementList>(name,"jet",kCyan);
159 elist->GetContainer()->SetRnrSelf(false);
160 elist->GetContainer()->SetRnrChildren(false);
161 elements.push_back(elist);
162 } else {
163 elements.push_back(new DelphesBranchElement<TEveElementList>(name,"jet",kYellow));
164 }
165 } else if(className=="Electron") {
166 tlist = new DelphesBranchElement<TEveTrackList>(name,"electron",kRed);
167 elements.push_back(tlist);
168 TEveTrackPropagator *trkProp = tlist->GetContainer()->GetPropagator();
169 trkProp->SetMagField(0., 0., -tk_Bz);
170 trkProp->SetMaxR(tk_radius);
171 trkProp->SetMaxZ(tk_length);
172 } else if(className=="Photon") {
173 tlist = new DelphesBranchElement<TEveTrackList>(name,"photon",kYellow);
174 elements.push_back(tlist);
175 TEveTrackPropagator *trkProp = tlist->GetContainer()->GetPropagator();
176 trkProp->SetMagField(0., 0., 0.);
177 trkProp->SetMaxR(tk_radius);
178 trkProp->SetMaxZ(tk_length);
179 } else if(className=="Muon") {
180 tlist = new DelphesBranchElement<TEveTrackList>(name,"muon",kGreen);
181 elements.push_back(tlist);
182 TEveTrackPropagator *trkProp = tlist->GetContainer()->GetPropagator();
183 trkProp->SetMagField(0., 0., -tk_Bz);
184 trkProp->SetMaxR(mu_radius);
185 trkProp->SetMaxZ(mu_length);
186 } else if(className=="MissingET") {
187 elements.push_back(new DelphesBranchElement<TEveElementList>(name,"vector",kViolet));
188 } else if(className=="GenParticle") {
189 tlist = new DelphesBranchElement<TEveTrackList>(name,"genparticle",kCyan);
190 elements.push_back(tlist);
191 tlist->GetContainer()->SetRnrSelf(false);
192 tlist->GetContainer()->SetRnrChildren(false);
193 TEveTrackPropagator *trkProp = tlist->GetContainer()->GetPropagator();
194 trkProp->SetMagField(0., 0., -tk_Bz);
195 trkProp->SetMaxR(tk_radius);
196 trkProp->SetMaxZ(tk_length);
197 } else {
198 continue;
199 }
200 arrays.push_back(treeReader_->UseBranch(name));
201 }
202 // second loop for tracks
203 for(Int_t b = 0; b<nBranches; ++b) {
204 TString input = branches[b*3].GetString();
205 TString name = branches[b*3+1].GetString();
206 TString className = branches[b*3+2].GetString();
207 if(className=="Track") {
208 if(input.Contains("eflow",TString::kIgnoreCase) || name.Contains("eflow",TString::kIgnoreCase)) continue; //no eflow
209 tlist = new DelphesBranchElement<TEveTrackList>(name,"track",kBlue);
210 elements.push_back(tlist);
211 TEveTrackPropagator *trkProp = tlist->GetContainer()->GetPropagator();
212 trkProp->SetMagField(0., 0., -tk_Bz);
213 trkProp->SetMaxR(tk_radius);
214 trkProp->SetMaxZ(tk_length);
215 arrays.push_back(treeReader_->UseBranch(name));
216 }
217 }
218}
219
220void DelphesEventDisplay::load_event()
221{
222 // Load event specified in global event_id_.
223 // The contents of previous event are removed.
224
225 // safety
226 if(event_id_ >= treeReader_->GetEntries() || event_id_<0 ) return;
227
228 //TODO move this to the status bar ???
229 printf("Loading event %d.\n", event_id_);
230
231 // clear the previous event
232 gEve->GetViewers()->DeleteAnnotations();
233 for(std::vector<DelphesBranchBase*>::iterator data=elements_.begin();data<elements_.end();++data) {
234 (*data)->Reset();
235 }
236
237 // Load selected branches with data from specified event
238 treeReader_->ReadEntry(event_id_);
239
240 // loop over selected branches, and apply the proper recipe to fill the collections.
241 // this is basically to loop on arrays_ to fill elements_.
242 std::vector<TClonesArray*>::iterator data = arrays_.begin();
243 std::vector<DelphesBranchBase*>::iterator element = elements_.begin();
244 for(; data<arrays_.end() && element<elements_.end(); ++data, ++element) {
245 TString type = (*element)->GetType();
246 // branch on the element type
247 if(type=="tower") delphes_read_towers(*data,*element);
248 else if(type=="track" || type=="photon" || type=="electron" || type=="muon" || type=="genparticle") delphes_read_tracks(*data,*element);
249 else if(type=="jet") delphes_read_jets(*data,*element);
250 else if(type=="vector") delphes_read_vectors(*data,*element);
251 }
252
253 // update display
254 TEveElement* top = (TEveElement*)gEve->GetCurrentEvent();
255 delphesDisplay_->DestroyEventRPhi();
256 delphesDisplay_->ImportEventRPhi(top);
257 delphesDisplay_->DestroyEventRhoZ();
258 delphesDisplay_->ImportEventRhoZ(top);
259 //update_html_summary();
260
261 gEve->Redraw3D(kFALSE, kTRUE);
262}
263
264void DelphesEventDisplay::delphes_read_towers(TClonesArray* data, DelphesBranchBase* element) {
265 DelphesCaloData* container = dynamic_cast<DelphesBranchElement<DelphesCaloData>*>(element)->GetContainer();
266 assert(container);
267 // Loop over all towers
268 TIter itTower(data);
269 Tower *tower;
270 while((tower = (Tower *) itTower.Next()))
271 {
272 container->AddTower(tower->Edges[0], tower->Edges[1], tower->Edges[2], tower->Edges[3]);
273 container->FillSlice(0, tower->Eem);
274 container->FillSlice(1, tower->Ehad);
275 }
276 container->DataChanged();
277}
278
279void DelphesEventDisplay::delphes_read_tracks(TClonesArray* data, DelphesBranchBase* element) {
280 TEveTrackList* container = dynamic_cast<DelphesBranchElement<TEveTrackList>*>(element)->GetContainer();
281 assert(container);
282 TString type = element->GetType();
283 TIter itTrack(data);
284 Int_t counter = 0;
285 TEveTrack *eveTrack;
286 TEveTrackPropagator *trkProp = container->GetPropagator();
287 if(type=="track") {
288 // Loop over all tracks
289 Track *track;
290 while((track = (Track *) itTrack.Next())) {
291 TParticle pb(track->PID, 1, 0, 0, 0, 0,
292 track->P4().Px(), track->P4().Py(),
293 track->P4().Pz(), track->P4().E(),
294 track->X, track->Y, track->Z, 0.0);
295 eveTrack = new TEveTrack(&pb, counter, trkProp);
296 eveTrack->SetName(Form("%s [%d]", pb.GetName(), counter++));
297 eveTrack->SetStdTitle();
298 eveTrack->SetAttLineAttMarker(container);
299 container->AddElement(eveTrack);
300 eveTrack->SetLineColor(element->GetColor());
301 eveTrack->MakeTrack();
302 }
303 } else if(type=="electron") {
304 // Loop over all electrons
305 Electron *electron;
306 while((electron = (Electron *) itTrack.Next())) {
307 TParticle pb(electron->Charge<0?11:-11, 1, 0, 0, 0, 0,
308 electron->P4().Px(), electron->P4().Py(),
309 electron->P4().Pz(), electron->P4().E(),
310 0., 0., 0., 0.);
311
312 eveTrack = new TEveTrack(&pb, counter, trkProp);
313 eveTrack->SetName(Form("%s [%d]", pb.GetName(), counter++));
314 eveTrack->SetStdTitle();
315 eveTrack->SetAttLineAttMarker(container);
316 container->AddElement(eveTrack);
317 eveTrack->SetLineColor(element->GetColor());
318 eveTrack->MakeTrack();
319 }
320 } else if(type=="muon") {
321 // Loop over all muons
322 Muon *muon;
323 while((muon = (Muon *) itTrack.Next())) {
324 TParticle pb(muon->Charge<0?13:-13, 1, 0, 0, 0, 0,
325 muon->P4().Px(), muon->P4().Py(),
326 muon->P4().Pz(), muon->P4().E(),
327 0., 0., 0., 0.);
328
329 eveTrack = new TEveTrack(&pb, counter, trkProp);
330 eveTrack->SetName(Form("%s [%d]", pb.GetName(), counter++));
331 eveTrack->SetStdTitle();
332 eveTrack->SetAttLineAttMarker(container);
333 container->AddElement(eveTrack);
334 eveTrack->SetLineColor(element->GetColor());
335 eveTrack->MakeTrack();
336 }
337 } else if(type=="photon") {
338 // Loop over all photons
339 Photon *photon;
340 while((photon = (Photon *) itTrack.Next())) {
341 TParticle pb(22, 1, 0, 0, 0, 0,
342 photon->P4().Px(), photon->P4().Py(),
343 photon->P4().Pz(), photon->P4().E(),
344 0., 0., 0., 0.);
345 eveTrack = new TEveTrack(&pb, counter, trkProp);
346 eveTrack->SetName(Form("%s [%d]", pb.GetName(), counter++));
347 eveTrack->SetStdTitle();
348 eveTrack->SetAttLineAttMarker(container);
349 eveTrack->SetLineStyle(7);
350 container->AddElement(eveTrack);
351 eveTrack->SetLineColor(element->GetColor());
352 eveTrack->MakeTrack();
353 }
354 } else if(type=="genparticle") {
355 // Loop over all particles
356 GenParticle *particle;
357 while((particle = (GenParticle *) itTrack.Next())) {
358 TParticle pb(particle->PID, particle->Status, particle->M1, particle->M2, particle->D1, particle->D2,
359 particle->P4().Px(), particle->P4().Py(),
360 particle->P4().Pz(), particle->P4().E(),
361 particle->X, particle->Y, particle->Z, particle->T);
362 eveTrack = new TEveTrack(&pb, counter, trkProp);
363 eveTrack->SetName(Form("%s [%d]", pb.GetName(), counter++));
364 eveTrack->SetStdTitle();
365 eveTrack->SetAttLineAttMarker(container);
366 container->AddElement(eveTrack);
367 eveTrack->SetLineColor(element->GetColor());
368 if(particle->Charge==0) eveTrack->SetLineStyle(7);
369 eveTrack->MakeTrack();
370 }
371 }
372}
373
374void DelphesEventDisplay::delphes_read_jets(TClonesArray* data, DelphesBranchBase* element) {
375 TEveElementList* container = dynamic_cast<DelphesBranchElement<TEveElementList>*>(element)->GetContainer();
376 assert(container);
377 TIter itJet(data);
378 Jet *jet;
379 TEveJetCone *eveJetCone;
380 // Loop over all jets
381 Int_t counter = 0;
382 while((jet = (Jet *) itJet.Next()))
383 {
384 eveJetCone = new TEveJetCone();
385 eveJetCone->SetTitle(Form("jet [%d]: Pt=%f, Eta=%f, \nPhi=%f, M=%f",counter,jet->PT, jet->Eta, jet->Phi, jet->Mass));
386 eveJetCone->SetName(Form("jet [%d]", counter++));
387 eveJetCone->SetMainTransparency(60);
388 eveJetCone->SetLineColor(element->GetColor());
389 eveJetCone->SetFillColor(element->GetColor());
390 eveJetCone->SetCylinder(tkRadius_ - 10, tkHalfLength_ - 10);
391 eveJetCone->SetPickable(kTRUE);
392 eveJetCone->AddEllipticCone(jet->Eta, jet->Phi, jet->DeltaEta, jet->DeltaPhi);
393 container->AddElement(eveJetCone);
394 }
395}
396
397void DelphesEventDisplay::delphes_read_vectors(TClonesArray* data, DelphesBranchBase* element) {
398 TEveElementList* container = dynamic_cast<DelphesBranchElement<TEveElementList>*>(element)->GetContainer();
399 assert(container);
400 TIter itMet(data);
401 MissingET *MET;
402 TEveArrow *eveMet;
403 // Missing Et
404 Double_t maxPt = 50.;
405 // TODO to be changed as we don't have access to maxPt anymore. MET scale could be a general parameter set in GUI
406 while((MET = (MissingET*) itMet.Next())) {
407 eveMet = new TEveArrow((tkRadius_ * MET->MET/maxPt)*cos(MET->Phi), (tkRadius_ * MET->MET/maxPt)*sin(MET->Phi), 0., 0., 0., 0.);
408 eveMet->SetMainColor(element->GetColor());
409 eveMet->SetTubeR(0.04);
410 eveMet->SetConeR(0.08);
411 eveMet->SetConeL(0.10);
412 eveMet->SetPickable(kTRUE);
413 eveMet->SetName("Missing Et");
414 eveMet->SetTitle(Form("Missing Et (%.1f GeV)",MET->MET));
415 container->AddElement(eveMet);
416 }
417}
418
419/******************************************************************************/
420// GUI
421/******************************************************************************/
422
423void DelphesEventDisplay::make_gui()
424{
425 // Create minimal GUI for event navigation.
426 // TODO: better GUI could be made based on the ch15 of the manual (Writing a GUI)
427
428 // add a tab on the left
429 TEveBrowser* browser = gEve->GetBrowser();
430 browser->StartEmbedding(TRootBrowser::kLeft);
431
432 // set the main title
433 TGMainFrame* frmMain = new TGMainFrame(gClient->GetRoot(), 1000, 600);
434 frmMain->SetWindowName("Delphes Event Display");
435 frmMain->SetCleanup(kDeepCleanup);
436
437 // build the navigation menu
438 TGHorizontalFrame* hf = new TGHorizontalFrame(frmMain);
439 {
440 TString icondir;
441 if(gSystem->Getenv("ROOTSYS"))
442 icondir = Form("%s/icons/", gSystem->Getenv("ROOTSYS"));
443 if(!gSystem->OpenDirectory(icondir))
444 icondir = Form("%s/icons/", (const char*)gSystem->GetFromPipe("root-config --etcdir") );
445 TGPictureButton* b = 0;
446
447 b = new TGPictureButton(hf, gClient->GetPicture(icondir+"GoBack.gif"));
448 hf->AddFrame(b);
449 b->Connect("Clicked()", "DelphesEventDisplay", this, "Bck()");
450
451 b = new TGPictureButton(hf, gClient->GetPicture(icondir+"GoForward.gif"));
452 hf->AddFrame(b);
453 b->Connect("Clicked()", "DelphesEventDisplay", this, "Fwd()");
454 }
455 frmMain->AddFrame(hf);
456 frmMain->MapSubwindows();
457 frmMain->Resize();
458 frmMain->MapWindow();
459 browser->StopEmbedding();
460 browser->SetTabTitle("Event Control", 0);
461}
462
Note: See TracBrowser for help on using the repository browser.