Fork me on GitHub

source: git/examples/geometry.C@ c4d879d

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

Small cleaning
  • Property mode set to 100644
File size: 33.2 KB
Line 
1#include <set>
2#include <map>
3#include <utility>
4#include <vector>
5#include <algorithm>
6#include <sstream>
7#include "TGeoManager.h"
8#include "TGeoVolume.h"
9#include "TGeoMedium.h"
10#include "TGeoNode.h"
11#include "TGeoCompositeShape.h"
12#include "TGeoMatrix.h"
13#include "TGeoTube.h"
14#include "TGeoCone.h"
15#include "TGeoArb8.h"
16//#include "external/ExRootAnalysis/ExRootConfReader.h"
17//#include "external/ExRootAnalysis/ExRootTreeReader.h"
18//#include "display/DelphesCaloData.h"
19//#include "display/DelphesDisplay.h"
20//#include "classes/DelphesClasses.h"
21#include "TF2.h"
22#include "TH1F.h"
23#include "TChain.h"
24#include "TEveElement.h"
25#include "TEveJetCone.h"
26#include "TEveTrack.h"
27#include "TEveTrackPropagator.h"
28#include "TEveCalo.h"
29#include "TMath.h"
30#include "TSystem.h"
31#include "TEveManager.h"
32#include "TEveGeoNode.h"
33#include "TEveTrans.h"
34#include "TEveViewer.h"
35#include "TEveBrowser.h"
36#include "TRootBrowser.h"
37#include "TGLViewer.h"
38#include "TGButton.h"
39#include "TCollection.h"
40#include "TClonesArray.h"
41#include "TGLClip.h"
42
43/*
44 * assembly.C: sauvegarde as shape-extract -> implement in the geometry class (read/write)
45 * histobrowser.C: intégration d'histogrammes dans le display (on pourrait avoir Pt, eta, phi pour les principales collections)
46 * also from alice_esd: summary html table
47 *
48 */
49using namespace std;
50
51// Forward declarations.
52class Delphes3DGeometry;
53class ExRootTreeReader;
54class DelphesCaloData;
55class DelphesDisplay;
56void make_gui();
57void load_event();
58void delphes_read();
59
60// Configuration and global variables.
61Int_t event_id = 0; // Current event id.
62Double_t gRadius = 1.29;
63Double_t gHalfLength = 3.0;
64Double_t gBz = 3.8;
65
66TAxis *gEtaAxis = 0;
67TAxis *gPhiAxis = 0;
68
69TChain gChain("Delphes");
70
71ExRootTreeReader *gTreeReader = 0;
72
73TClonesArray *gBranchTower = 0;
74TClonesArray *gBranchTrack = 0;
75TClonesArray *gBranchJet = 0;
76
77DelphesCaloData *gCaloData = 0;
78TEveElementList *gJetList = 0;
79TEveTrackList *gTrackList = 0;
80
81DelphesDisplay *gDelphesDisplay = 0;
82
83/******************************************************************************/
84// Construction of the geometry
85/******************************************************************************/
86
87// TODO: asymmetric detector
88
89class Delphes3DGeometry {
90 public:
91 Delphes3DGeometry(TGeoManager *geom = NULL);
92 ~Delphes3DGeometry() {}
93
94 void readFile(const char* filename, const char* ParticlePropagator="ParticlePropagator",
95 const char* TrackingEfficiency="ChargedHadronTrackingEfficiency",
96 const char* MuonEfficiency="MuonEfficiency",
97 const char* Calorimeters="Calorimeter");
98
99 void loadFromFile(const char* filename, const char* name="DelphesGeometry");
100 void save(const char* filename, const char* name="DelphesGeometry");
101
102 void setContingency(Double_t contingency) { contingency_ = contingency; }
103 void setCaloBarrelThickness(Double_t thickness) { calo_barrel_thickness_ = thickness; }
104 void setCaloEndcapThickness(Double_t thickness) { calo_endcap_thickness_ = thickness; }
105 void setMuonSystemThickness(Double_t thickness) { muonSystem_thickn_ = thickness; }
106
107 TGeoVolume* getDetector(bool withTowers = true);
108
109 Double_t getTrackerRadius() const { return tk_radius_; }
110 Double_t getTrackerHalfLength() const { return tk_length_; }
111 Double_t getBField() const { return tk_Bz_; }
112 std::pair<TAxis*, TAxis*> getCaloAxes() { return std::make_pair(etaAxis_,phiAxis_); }
113
114 private:
115 std::pair<Double_t, Double_t> addTracker(TGeoVolume *top);
116 std::pair<Double_t, Double_t> addCalorimeter(TGeoVolume *top, const char* name, Double_t innerBarrelRadius, Double_t innerBarrelLength, set< pair<Double_t, Int_t> >& caloBinning);
117 std::pair<Double_t, Double_t> addMuonDets(TGeoVolume *top, const char* name, Double_t innerBarrelRadius, Double_t innerBarrelLength);
118 void addCaloTowers(TGeoVolume *top, const char* name, Double_t innerBarrelRadius, Double_t innerBarrelLength, set< pair<Double_t, Int_t> >& caloBinning);
119
120 private:
121
122 TGeoManager *geom_;
123
124 TGeoMedium *vacuum_;
125 TGeoMedium *tkmed_;
126 TGeoMedium *calomed_;
127 TGeoMedium *mudetmed_;
128
129 TAxis* etaAxis_;
130 TAxis* phiAxis_;
131
132 Double_t contingency_;
133 Double_t calo_barrel_thickness_;
134 Double_t calo_endcap_thickness_;
135 Double_t muonSystem_thickn_;
136 Double_t tk_radius_;
137 Double_t tk_length_;
138 Double_t tk_etamax_;
139 Double_t tk_Bz_;
140
141 std::vector<std::string> calorimeters_;
142 std::vector<std::string> muondets_;
143
144 std::map<std::string, Double_t> muonSystem_etamax_;
145 std::map<std::string, set< pair<Double_t, Int_t> > > caloBinning_;
146
147};
148
149Delphes3DGeometry::Delphes3DGeometry(TGeoManager *geom) {
150
151 //--- the geometry manager
152 geom_ = geom==NULL? gGeoManager : geom;
153 //gGeoManager->DefaultColors();
154
155 //--- define some materials
156 TGeoMaterial *matVacuum = new TGeoMaterial("Vacuum", 0,0,0);
157 TGeoMaterial *matAl = new TGeoMaterial("Al", 26.98,13,2.7); // placeholder
158 //TODO: create different materials for different subdetectors???
159 matVacuum->SetTransparency(85); //TODO: tune
160 matAl->SetTransparency(85); //TODO: tune
161
162 //--- define some media
163 TGeoMedium *Vacuum = new TGeoMedium("Vacuum",1, matVacuum);
164 TGeoMedium *Al = new TGeoMedium("Root Material",2, matAl);
165 vacuum_ = Vacuum;
166 tkmed_ = Vacuum; // placeholder
167 calomed_ = Al; // placeholder
168 mudetmed_ = Al; // placeholder
169
170 // custom parameters
171 contingency_ = 10.;
172 calo_barrel_thickness_ = 50.;
173 calo_endcap_thickness_ = 75.;
174 muonSystem_thickn_ = 10.;
175
176 // read these parameters from the Delphes Card (with default values)
177 etaAxis_ = NULL;
178 phiAxis_ = NULL;
179 tk_radius_ = 120.;
180 tk_length_ = 150.;
181 tk_etamax_ = 3.0;
182 tk_Bz_ = 1.;
183}
184
185void Delphes3DGeometry::readFile(const char *configFile,
186 const char* ParticlePropagator, const char* TrackingEfficiency,
187 const char* MuonEfficiency, const char* Calorimeters) {
188
189 ExRootConfReader *confReader = new ExRootConfReader;
190 confReader->ReadFile(configFile);
191
192 tk_radius_ = confReader->GetDouble(Form("%s::Radius",ParticlePropagator), 1.0)*100.; // tk_radius
193 tk_length_ = confReader->GetDouble(Form("%s::HalfLength",ParticlePropagator), 3.0)*100.; // tk_length
194 tk_Bz_ = confReader->GetDouble("ParticlePropagator::Bz", 0.0); // tk_Bz
195
196 {
197 TString tkEffFormula = confReader->GetString(Form("%s::EfficiencyFormula",TrackingEfficiency),"abs(eta)<3.0");
198 tkEffFormula.ReplaceAll("pt","x");
199 tkEffFormula.ReplaceAll("eta","y");
200 tkEffFormula.ReplaceAll("phi","0.");
201 TF2* tkEffFunction = new TF2("tkEff",tkEffFormula,0,1000,-10,10);
202 TH1F etaHisto("eta","eta",100,5.,-5.);
203 Double_t pt,eta;
204 for(int i=0;i<1000;++i) {
205 tkEffFunction->GetRandom2(pt,eta);
206 etaHisto.Fill(eta);
207 }
208 Int_t bin = -1;
209 bin = etaHisto.FindFirstBinAbove(0.5);
210 Double_t etamin = (bin>-1) ? etaHisto.GetBinLowEdge(bin) : -10.;
211 bin = etaHisto.FindLastBinAbove(0.5);
212 Double_t etamax = (bin>-1) ? etaHisto.GetBinLowEdge(bin+1) : -10.;
213 tk_etamax_ = TMath::Max(fabs(etamin),fabs(etamax)); // tk_etamax
214 delete tkEffFunction;
215 }
216
217 {
218 muondets_.push_back("muons");
219 TString muonEffFormula = confReader->GetString(Form("%s::EfficiencyFormula",MuonEfficiency),"abs(eta)<2.0");
220 muonEffFormula.ReplaceAll("pt","x");
221 muonEffFormula.ReplaceAll("eta","y");
222 muonEffFormula.ReplaceAll("phi","0.");
223 TF2* muEffFunction = new TF2("muEff",muonEffFormula,0,1000,-10,10);
224 TH1F etaHisto("eta2","eta2",100,5.,-5.);
225 Double_t pt,eta;
226 for(int i=0;i<1000;++i) {
227 muEffFunction->GetRandom2(pt,eta);
228 etaHisto.Fill(eta);
229 }
230 Int_t bin = -1;
231 bin = etaHisto.FindFirstBinAbove(0.5);
232 Double_t etamin = (bin>-1) ? etaHisto.GetBinLowEdge(bin) : -10.;
233 bin = etaHisto.FindLastBinAbove(0.5);
234 Double_t etamax = (bin>-1) ? etaHisto.GetBinLowEdge(bin+1) : -10.;
235 muonSystem_etamax_["muons"] = TMath::Max(fabs(etamin),fabs(etamax)); // muonSystem_etamax
236 delete muEffFunction;
237 }
238
239 std::string s(Calorimeters);
240 std::replace( s.begin(), s.end(), ',', ' ' );
241 std::istringstream stream( s );
242 std::string word;
243 while (stream >> word) calorimeters_.push_back(word);
244
245 caloBinning_.clear(); // calo binning
246 for(std::vector<std::string>::const_iterator calo=calorimeters_.begin();calo!=calorimeters_.end(); ++calo) {
247 set< pair<Double_t, Int_t> > caloBinning;
248 ExRootConfParam paramEtaBins, paramPhiBins;
249 ExRootConfParam param = confReader->GetParam(Form("%s::EtaPhiBins",calo->c_str()));
250 Int_t size = param.GetSize();
251 for(int i = 0; i < size/2; ++i) {
252 paramEtaBins = param[i*2];
253 paramPhiBins = param[i*2+1];
254 assert(paramEtaBins.GetSize()==1);
255 caloBinning.insert(std::make_pair(paramEtaBins[0].GetDouble(),paramPhiBins.GetSize()-1));
256 }
257 caloBinning_[*calo] = caloBinning;
258 }
259
260 set< pair<Double_t, Int_t> > caloBinning = caloBinning_[*calorimeters_.begin()];
261 Double_t *etaBins = new Double_t[caloBinning.size()]; // note that this is the eta binning of the first calo
262 unsigned int ii = 0;
263 for(set< pair<Double_t, Int_t> >::const_iterator itEtaSet = caloBinning.begin(); itEtaSet != caloBinning.end(); ++itEtaSet) {
264 etaBins[ii++] = itEtaSet->first;
265 }
266 etaAxis_ = new TAxis(caloBinning.size() - 1, etaBins);
267 phiAxis_ = new TAxis(72, -TMath::Pi(), TMath::Pi()); // note that this is fixed while #phibins could vary, also with eta, which doesn't seem possible in ROOT
268
269 delete confReader;
270
271}
272
273TGeoVolume* Delphes3DGeometry::getDetector(bool withTowers) {
274 // compute the envelope
275 Double_t system_radius = tk_radius_+calo_barrel_thickness_+3*contingency_;
276 Double_t system_length = tk_length_+contingency_+(contingency_+calo_endcap_thickness_)*calorimeters_.size()+contingency_;
277 // the detector volume
278 TGeoVolume *top = geom_->MakeBox("Delphes3DGeometry", vacuum_, system_radius, system_radius, system_length);
279 // build the detector
280 std::pair<Double_t, Double_t> limits = addTracker(top);
281 Double_t radius = limits.first;
282 Double_t length = limits.second;
283 for(std::vector<std::string>::const_iterator calo = calorimeters_.begin(); calo != calorimeters_.end(); ++calo) {
284 limits = addCalorimeter(top,calo->c_str(),radius,length,caloBinning_[*calo]);
285 if (withTowers) {
286 addCaloTowers(top,calo->c_str(),radius,length,caloBinning_[*calo]);
287 }
288 radius = limits.first;
289 length = limits.second;
290 }
291 for(std::vector<std::string>::const_iterator muon = muondets_.begin(); muon != muondets_.end(); ++muon) {
292 limits = addMuonDets(top, muon->c_str(), radius, length);
293 radius = limits.first;
294 length = limits.second;
295 }
296 // return the result
297 return top;
298}
299
300std::pair<Double_t, Double_t> Delphes3DGeometry::addTracker(TGeoVolume *top) {
301 // tracker: a cylinder with two cones substracted
302 new TGeoCone("forwardTkAcceptance",(tk_length_/2.+0.05),0.,tk_radius_,(tk_length_)*2.*exp(-tk_etamax_)/(1-exp(-2.*tk_etamax_)),tk_radius_);
303 TGeoTranslation *tr1 = new TGeoTranslation("tkacc1",0., 0., tk_length_/2.);
304 tr1->RegisterYourself();
305 TGeoRotation *negz = new TGeoRotation("tknegz",0,180,0);
306 negz->RegisterYourself();
307 TGeoCombiTrans *tr2 = new TGeoCombiTrans("tkacc2",0.,0.,-tk_length_/2.,negz);
308 tr2->RegisterYourself();
309 TGeoCompositeShape* tracker_cs = new TGeoCompositeShape("tracker_cs","forwardTkAcceptance:tkacc1+forwardTkAcceptance:tkacc2");
310 TGeoVolume *tracker = new TGeoVolume("tracker",tracker_cs,tkmed_);
311 tracker->SetLineColor(kYellow);
312 top->AddNode(tracker,1);
313 return std::make_pair(tk_radius_,tk_length_);
314}
315
316std::pair<Double_t, Double_t> Delphes3DGeometry::addCalorimeter(TGeoVolume *top, const char* name,
317 Double_t innerBarrelRadius, Double_t innerBarrelLength, set< pair<Double_t, Int_t> >& caloBinning) {
318 // parameters derived from the inputs
319 Double_t calo_endcap_etamax = TMath::Max(fabs(caloBinning.begin()->first),fabs(caloBinning.rbegin()->first));
320 Double_t calo_barrel_innerRadius = innerBarrelRadius+contingency_;
321 Double_t calo_barrel_length = innerBarrelLength + calo_barrel_thickness_;
322 Double_t calo_endcap_etamin = -log(innerBarrelRadius/(2*innerBarrelLength));
323 Double_t calo_endcap_innerRadius1 = innerBarrelLength*2.*exp(-calo_endcap_etamax)/(1-exp(-2.*calo_endcap_etamax));
324 Double_t calo_endcap_innerRadius2 = (innerBarrelLength+calo_endcap_thickness_)*2.*exp(-calo_endcap_etamax)/(1-exp(-2.*calo_endcap_etamax));
325 Double_t calo_endcap_outerRadius1 = innerBarrelRadius;
326 Double_t calo_endcap_outerRadius2 = innerBarrelRadius+calo_barrel_thickness_;
327 Double_t calo_endcap_coneThickness = TMath::Min(calo_barrel_thickness_ * (1-exp(-2.*calo_endcap_etamin)) / (2.*exp(-calo_endcap_etamin)), calo_endcap_thickness_);
328 Double_t calo_endcap_diskThickness = TMath::Max(0.,calo_endcap_thickness_-calo_endcap_coneThickness);
329
330 // calorimeters: tube truncated in eta + cones
331 new TGeoTube(Form("%s_barrel_cylinder",name),calo_barrel_innerRadius,calo_barrel_innerRadius+calo_barrel_thickness_,calo_barrel_length);
332 new TGeoCone(Form("%s_endcap_cone",name),calo_endcap_coneThickness/2.,calo_endcap_innerRadius1,calo_endcap_outerRadius1,calo_endcap_innerRadius2,calo_endcap_outerRadius2);
333 new TGeoTube(Form("%s_endcap_disk",name),calo_endcap_innerRadius2,tk_radius_+calo_barrel_thickness_,calo_endcap_diskThickness/2.);
334 TGeoTranslation *tr1 = new TGeoTranslation(Form("%s_tr1",name),0., 0., (calo_endcap_coneThickness+calo_endcap_diskThickness)/2.);
335 tr1->RegisterYourself();
336 TGeoCompositeShape *calo_endcap_cs = new TGeoCompositeShape(Form("%s_endcap_cs",name),Form("%s_endcap_cone+%s_endcap_disk:%s_tr1",name,name,name));
337 TGeoTranslation *trc1 = new TGeoTranslation(Form("%s_endcap1_position",name),0.,0., innerBarrelLength+calo_endcap_coneThickness/2.);
338 trc1->RegisterYourself();
339 TGeoRotation *negz = new TGeoRotation(Form("%s_negz",name),0,180,0);
340 TGeoCombiTrans *trc2 = new TGeoCombiTrans(Form("%s_endcap2_position",name),0.,0.,-(innerBarrelLength+calo_endcap_coneThickness/2.),negz);
341 trc2->RegisterYourself();
342 TGeoTranslation *trc1c = new TGeoTranslation(Form("%s_endcap1_position_cont",name),0.,0., innerBarrelLength+calo_endcap_coneThickness/2.+contingency_);
343 trc1c->RegisterYourself();
344 TGeoCombiTrans *trc2c = new TGeoCombiTrans(Form("%s_endcap2_position_cont",name),0.,0.,-(innerBarrelLength+calo_endcap_coneThickness/2.)-contingency_,negz);
345 trc2c->RegisterYourself();
346 TGeoVolume *calo_endcap = new TGeoVolume(Form("%s_endcap",name),calo_endcap_cs,calomed_);
347 TGeoCompositeShape *calo_barrel_cs = new TGeoCompositeShape(Form("%s_barrel_cs",name),
348 Form("%s_barrel_cylinder-%s_endcap_cs:%s_endcap1_position-%s_endcap_cs:%s_endcap2_position",name,name,name,name,name));
349 TGeoVolume *calo_barrel = new TGeoVolume(Form("%s_barrel",name),calo_barrel_cs,calomed_);
350 calo_endcap->SetLineColor(kViolet);
351 calo_endcap->SetFillColor(kViolet);
352 calo_barrel->SetLineColor(kRed);
353 top->AddNode(calo_endcap,1,trc1c);
354 top->AddNode(calo_endcap,2,trc2c);
355 top->AddNode(calo_barrel,1);
356 return std::make_pair(calo_barrel_innerRadius+calo_barrel_thickness_,innerBarrelLength+calo_endcap_thickness_+contingency_);
357}
358
359std::pair<Double_t, Double_t> Delphes3DGeometry::addMuonDets(TGeoVolume *top, const char* name, Double_t innerBarrelRadius, Double_t innerBarrelLength) {
360 // muon system: tube + disks
361 Double_t muonSystem_radius = innerBarrelRadius + contingency_;
362 Double_t muonSystem_length = innerBarrelLength + contingency_;
363 Double_t muonSystem_rmin = muonSystem_length*2.*exp(-muonSystem_etamax_[name])/(1-exp(-2.*muonSystem_etamax_[name]));
364 TGeoVolume *muon_barrel = geom_->MakeTube(Form("%s_barrel",name),mudetmed_,muonSystem_radius,muonSystem_radius+muonSystem_thickn_,muonSystem_length);
365 muon_barrel->SetLineColor(kBlue);
366 top->AddNode(muon_barrel,1);
367 TGeoVolume *muon_endcap = geom_->MakeTube(Form("%s_endcap",name),mudetmed_,muonSystem_rmin,muonSystem_radius+muonSystem_thickn_,muonSystem_thickn_/2.);
368 muon_endcap->SetLineColor(kBlue);
369 TGeoTranslation *trm1 = new TGeoTranslation(Form("%sEndcap1_position",name),0.,0.,muonSystem_length);
370 trm1->RegisterYourself();
371 TGeoTranslation *trm2 = new TGeoTranslation(Form("%sEndcap2_position",name),0.,0.,-muonSystem_length);
372 trm1->RegisterYourself();
373 top->AddNode(muon_endcap,1,trm1);
374 top->AddNode(muon_endcap,2,trm2);
375 return std::make_pair(muonSystem_radius,muonSystem_length);
376}
377
378void Delphes3DGeometry::addCaloTowers(TGeoVolume *top, const char* name,
379 Double_t innerBarrelRadius, Double_t innerBarrelLength, set< pair<Double_t, Int_t> >& caloBinning) {
380
381 TGeoVolume* calo_endcap = top->GetNode(Form("%s_endcap_1",name))->GetVolume();
382 TGeoVolume* calo_barrel = top->GetNode(Form("%s_barrel_1",name))->GetVolume();
383 Double_t calo_endcap_etamin = -log(innerBarrelRadius/(2*innerBarrelLength));
384 Double_t calo_endcap_coneThickness = TMath::Min(calo_barrel_thickness_ * (1-exp(-2.*calo_endcap_etamin)) / (2.*exp(-calo_endcap_etamin)), calo_endcap_thickness_);
385
386 // calo towers in the barrel
387 Double_t vertices[16] = {0.,0.,0.,0.,0.,0.,0.,0.}; // summit of the pyramid
388 Double_t R = tk_radius_ + contingency_+(contingency_+calo_barrel_thickness_)*calorimeters_.size(); // radius of the muons system = height of the pyramid
389 Int_t nEtaBins = caloBinning.size();
390 // this rotation is to make the tower point "up"
391 TGeoRotation* initTowerRot = new TGeoRotation(Form("%s_initTowerRot",name),0.,90.,0.);
392 TGeoCombiTrans* initTower = new TGeoCombiTrans(Form("%s_initTower",name),0.,-R/2.,0.,initTowerRot);
393 initTower->RegisterYourself();
394 // eta bins... we build one pyramid per eta slice and then translate it nphi times.
395 // phi bins represented by rotations around z
396 Double_t *y = new Double_t[nEtaBins];
397 Double_t *dx = new Double_t[nEtaBins];
398 Int_t *nphi = new Int_t[nEtaBins];
399 Int_t etaslice = 0;
400 std::map<std::pair<int,int>, TGeoRotation*> phirotations;
401 for(set< pair<Double_t, Int_t> >::const_iterator bin=caloBinning.begin(); bin!=caloBinning.end();++bin) {
402 if(abs(bin->first)>calo_endcap_etamin) continue; // only in the barrel
403 nphi[etaslice] = bin->second;
404 y[etaslice] = 0.5*R*(1-exp(-2*bin->first))/exp(-bin->first);
405 Double_t phiRotationAngle = 360./nphi[etaslice];
406 dx[etaslice] = R*tan(TMath::Pi()*phiRotationAngle/360.);
407 for(int phislice=0;phislice<nphi[etaslice];++phislice) {
408 phirotations[make_pair(etaslice,phislice)] = new TGeoRotation(Form("%s_phi%d_%d",name,etaslice,phislice),phiRotationAngle*phislice,0.,0.);
409 phirotations[make_pair(etaslice,phislice)]->RegisterYourself();
410 }
411 ++etaslice;
412 }
413 nEtaBins = etaslice;
414 for(int i=0;i<nEtaBins-1;++i) { // loop on the eta slices
415 vertices[8] = -dx[i]; vertices[9] = y[i];
416 vertices[10] = -dx[i]; vertices[11] = y[i+1];
417 vertices[12] = dx[i]; vertices[13] = y[i+1];
418 vertices[14] = dx[i]; vertices[15] = y[i];
419 new TGeoArb8(Form("%s_tower%d",name,i),R/2., vertices); // tower in the proper eta slice, at phi=0
420 // intersection between the tower and the calo_barrel
421 TGeoCompositeShape *finaltower_cs = new TGeoCompositeShape(Form("%s_ftower%d_cs",name,i),Form("%s_tower%d:%s_initTower*%s_barrel_cs",name,i,name,name));
422 TGeoVolume *finaltower = new TGeoVolume(Form("%s_ftower%d",name,i),finaltower_cs,calomed_);
423 finaltower->SetLineColor(kRed);
424 for(int j=0;j<nphi[i];++j) { // loop on the phi slices
425 calo_barrel->AddNode(finaltower,j,phirotations[make_pair(i,j)]);
426 }
427 }
428 delete[] y;
429 delete[] dx;
430 delete[] nphi;
431 //the towers in the forward region
432 R = tk_length_+contingency_+(contingency_+calo_endcap_thickness_)*calorimeters_.size(); // Z of the muons system = height of the pyramid
433 nEtaBins = caloBinning.size();
434 // translation to bring the origin of the tower to (0,0,0) (well, not really as the endcap is not yet in place)
435 TGeoTranslation* towerdz = new TGeoTranslation(Form("%s_towerdz",name),0.,0.,R/2.-(innerBarrelLength+calo_endcap_coneThickness/2.));
436 towerdz->RegisterYourself();
437 // eta bins... we build one pyramid per eta slice and then translate it nphi times.
438 Double_t *r = new Double_t[nEtaBins];
439 nphi = new Int_t[nEtaBins];
440 etaslice = 0;
441 phirotations.clear();
442 for(set< pair<Double_t, Int_t> >::const_iterator bin=caloBinning.begin(); bin!=caloBinning.end();++bin) {
443 if(bin->first<calo_endcap_etamin) continue; // only in the + endcap
444 r[etaslice] = R*2*exp(-bin->first)/(1-exp(-2*bin->first));
445 nphi[etaslice] = bin->second;
446 Double_t phiRotationAngle = 360./nphi[etaslice];
447 for(int phislice=0;phislice<nphi[etaslice];++phislice) {
448 phirotations[make_pair(etaslice,phislice)] = new TGeoRotation(Form("%s_forward_phi%d_%d",name,etaslice,phislice),phiRotationAngle*phislice,0.,0.);
449 phirotations[make_pair(etaslice,phislice)]->RegisterYourself();
450 }
451 ++etaslice;
452 }
453 nEtaBins = etaslice;
454 for(int i=0;i<nEtaBins-1;++i) { // loop on the eta slices
455 vertices[8] = -r[i+1]*sin(TMath::Pi()/nphi[i]); vertices[9] = r[i+1]*cos(TMath::Pi()/nphi[i]);
456 vertices[10] = -r[i]*sin(TMath::Pi()/nphi[i]); vertices[11] = r[i]*cos(TMath::Pi()/nphi[i]);
457 vertices[12] = r[i]*sin(TMath::Pi()/nphi[i]); vertices[13] = r[i]*cos(TMath::Pi()/nphi[i]);
458 vertices[14] = r[i+1]*sin(TMath::Pi()/nphi[i]); vertices[15] = r[i+1]*cos(TMath::Pi()/nphi[i]);
459 new TGeoArb8(Form("%sfwdtower%d",name,i),R/2., vertices); // tower in the proper eta slice, at phi=0
460 // intersection between the tower and the calo_endcap
461 TGeoCompositeShape *finalfwdtower_cs = new TGeoCompositeShape(Form("%sffwdtower%d_cs",name,i),Form("%sfwdtower%d:%s_towerdz*%s_endcap_cs",name,i,name,name));
462 TGeoVolume *finalfwdtower = new TGeoVolume(Form("%sffwdtower%d",name,i),finalfwdtower_cs,calomed_);
463 finalfwdtower->SetLineColor(kViolet);
464 for(int j=0;j<nphi[i];++j) { // loop on the phi slices
465 calo_endcap->AddNode(finalfwdtower,j,phirotations[make_pair(i,j)]);
466 }
467 }
468 delete[] r;
469 delete[] nphi;
470}
471
472/******************************************************************************/
473// Initialization and steering functions
474/******************************************************************************/
475
476void delphes_event_display(const char *configFile, const char *inputFile, const Delphes3DGeometry& det3D)
477{
478
479 // initialize the application
480 TEveManager::Create(kTRUE, "IV");
481 TGeoManager* geom = gGeoManager;
482
483 // build the detector
484 gRadius = det3D.getTrackerRadius();
485 gHalfLength = det3D.getTrackerHalfLength();
486 gBz = det3D.getBField();
487 gEtaAxis = det3D.getCaloAxes().first;
488 gPhiAxis = det3D.getCaloAxes().second;
489
490 TGeoVolume* top = det3D.getDetector(false); //TODO: can this be set in the GUI?
491 geom->SetTopVolume(top);
492 TEveElementList *geometry = new TEveElementList("Geometry");
493 TEveGeoTopNode* trk = new TEveGeoTopNode(gGeoManager, top->FindNode("tracker_1"));
494 trk->SetVisLevel(6);
495 geometry->AddElement(trk);
496 TEveGeoTopNode* calo = new TEveGeoTopNode(gGeoManager, top->FindNode("Calorimeter_barrel_1"));
497 calo->SetVisLevel(3);
498 geometry->AddElement(calo);
499 calo = new TEveGeoTopNode(gGeoManager, top->FindNode("Calorimeter_endcap_1"));
500 calo->SetVisLevel(3);
501 calo->UseNodeTrans();
502 geometry->AddElement(calo);
503 calo = new TEveGeoTopNode(gGeoManager, top->FindNode("Calorimeter_endcap_2"));
504 calo->SetVisLevel(3);
505 calo->UseNodeTrans();
506 geometry->AddElement(calo);
507 TEveGeoTopNode* muon = new TEveGeoTopNode(gGeoManager, top->FindNode("muons_barrel_1"));
508 muon->SetVisLevel(4);
509 geometry->AddElement(muon);
510 muon = new TEveGeoTopNode(gGeoManager, top->FindNode("muons_endcap_1"));
511 muon->SetVisLevel(4);
512 muon->UseNodeTrans();
513 geometry->AddElement(muon);
514 muon = new TEveGeoTopNode(gGeoManager, top->FindNode("muons_endcap_2"));
515 muon->SetVisLevel(4);
516 muon->UseNodeTrans();
517 geometry->AddElement(muon);
518 //gGeoManager->DefaultColors();
519
520 // Create chain of root trees
521 gChain.Add(inputFile);
522
523 // Create object of class ExRootTreeReader
524 printf("*** Opening Delphes data file ***\n");
525 gTreeReader = new ExRootTreeReader(&gChain);
526
527 // Get pointers to branches
528 gBranchTower = gTreeReader->UseBranch("Tower");
529 gBranchTrack = gTreeReader->UseBranch("Track");
530 gBranchJet = gTreeReader->UseBranch("Jet");
531
532//TODO make it configurable, for more objects (or can we guess from the config?)
533// idea: for pf objects, we could use the TEveCompound to show track + cluster ??? (nice display but little meaning)
534// for MET and SHT, show an arrow (tooltip = title)
535// for electrons and muons, create additional track collections
536// for photons, use TEveStraightLineSet
537
538/*
539 add Branch Calorimeter/eflowTracks EFlowTrack Track
540 add Branch Calorimeter/eflowPhotons EFlowPhoton Tower
541 add Branch Calorimeter/eflowNeutralHadrons EFlowNeutralHadron Tower
542
543 add Branch GenJetFinder/jets GenJet Jet
544
545 add Branch UniqueObjectFinder/jets Jet Jet
546
547 add Branch UniqueObjectFinder/electrons Electron Electron
548 add Branch UniqueObjectFinder/photons Photon Photon
549 add Branch UniqueObjectFinder/muons Muon Muon
550
551 add Branch MissingET/momentum MissingET MissingET
552 add Branch ScalarHT/energy ScalarHT ScalarHT
553*/
554
555 // data
556 gCaloData = new DelphesCaloData(2);
557 gCaloData->RefSliceInfo(0).Setup("ECAL", 0.1, kRed);
558 gCaloData->RefSliceInfo(1).Setup("HCAL", 0.1, kBlue);
559 gCaloData->SetEtaBins(gEtaAxis);
560 gCaloData->SetPhiBins(gPhiAxis);
561 gCaloData->IncDenyDestroy();
562
563 gJetList = new TEveElementList("Jets");
564 gEve->AddElement(gJetList);
565
566 gTrackList = new TEveTrackList("Tracks");
567 gTrackList->SetMainColor(kBlue);
568 gTrackList->SetMarkerColor(kRed);
569 gTrackList->SetMarkerStyle(kCircle);
570 gTrackList->SetMarkerSize(0.5);
571 gEve->AddElement(gTrackList);
572
573 TEveTrackPropagator *trkProp = gTrackList->GetPropagator();
574 trkProp->SetMagField(0.0, 0.0, -gBz);
575 trkProp->SetMaxR(gRadius);
576 trkProp->SetMaxZ(gHalfLength);
577
578 // viewers and scenes
579
580 TEveCalo3D *calo3d = new TEveCalo3D(gCaloData);
581 calo3d->SetBarrelRadius(gRadius);
582 calo3d->SetEndCapPos(gHalfLength);
583
584 //gStyle->SetPalette(1, 0);
585 TEveCaloLego *lego = new TEveCaloLego(gCaloData);
586 lego->InitMainTrans();
587 lego->RefMainTrans().SetScale(TMath::TwoPi(), TMath::TwoPi(), TMath::Pi());
588 lego->SetAutoRebin(kFALSE);
589 lego->Set2DMode(TEveCaloLego::kValSizeOutline);
590
591 gDelphesDisplay = new DelphesDisplay;
592 gEve->AddGlobalElement(geometry);
593 gEve->AddGlobalElement(calo3d);
594 gDelphesDisplay->ImportGeomRPhi(geometry);
595 gDelphesDisplay->ImportCaloRPhi(calo3d);
596 gDelphesDisplay->ImportGeomRhoZ(geometry);
597 gDelphesDisplay->ImportCaloRhoZ(calo3d);
598 gDelphesDisplay->ImportCaloLego(lego);
599 gEve->Redraw3D(kTRUE);
600
601}
602
603//______________________________________________________________________________
604void load_event()
605{
606 // Load event specified in global event_id.
607 // The contents of previous event are removed.
608
609 printf("Loading event %d.\n", event_id);
610
611 gEve->GetViewers()->DeleteAnnotations();
612
613 if(gCaloData) gCaloData->ClearTowers();
614 if(gJetList) gJetList->DestroyElements();
615 if(gTrackList) gTrackList->DestroyElements();
616
617 delphes_read();
618
619 TEveElement* top = (TEveElement*)gEve->GetCurrentEvent();
620 gDelphesDisplay->DestroyEventRPhi();
621 gDelphesDisplay->ImportEventRPhi(top);
622 gDelphesDisplay->DestroyEventRhoZ();
623 gDelphesDisplay->ImportEventRhoZ(top);
624
625 //update_html_summary();
626
627 gEve->Redraw3D(kFALSE, kTRUE);
628}
629
630void delphes_read()
631{
632
633 TIter itTower(gBranchTower);
634 TIter itTrack(gBranchTrack);
635 TIter itJet(gBranchJet);
636
637 Tower *tower;
638 Track *track;
639 Jet *jet;
640
641 TEveJetCone *eveJetCone;
642 TEveTrack *eveTrack;
643
644 Int_t counter;
645
646 TEveTrackPropagator *trkProp = gTrackList->GetPropagator();
647 if(event_id >= gTreeReader->GetEntries()) return;
648
649 // Load selected branches with data from specified event
650 gTreeReader->ReadEntry(event_id);
651 // Loop over all towers
652 itTower.Reset();
653 while((tower = (Tower *) itTower.Next()))
654 {
655 gCaloData->AddTower(tower->Edges[0], tower->Edges[1], tower->Edges[2], tower->Edges[3]);
656 gCaloData->FillSlice(0, tower->Eem);
657 gCaloData->FillSlice(1, tower->Ehad);
658 }
659 gCaloData->DataChanged();
660 // Loop over all tracks
661 itTrack.Reset();
662 counter = 0;
663 while((track = (Track *) itTrack.Next()))
664 {
665 TParticle pb(track->PID, 1, 0, 0, 0, 0,
666 track->P4().Px(), track->P4().Py(),
667 track->P4().Pz(), track->P4().E(),
668 track->X, track->Y, track->Z, 0.0);
669
670 eveTrack = new TEveTrack(&pb, counter, trkProp);
671 eveTrack->SetName(Form("%s [%d]", pb.GetName(), counter++));
672 eveTrack->SetStdTitle();
673 eveTrack->SetAttLineAttMarker(gTrackList);
674
675 switch(TMath::Abs(track->PID))
676 {
677 case 11:
678 eveTrack->SetLineColor(kRed);
679 break;
680 case 13:
681 eveTrack->SetLineColor(kGreen);
682 break;
683 default:
684 eveTrack->SetLineColor(kBlue);
685 }
686 gTrackList->AddElement(eveTrack);
687 eveTrack->MakeTrack();
688 }
689 // Loop over all jets
690 itJet.Reset();
691 counter = 0;
692 while((jet = (Jet *) itJet.Next()))
693 {
694 eveJetCone = new TEveJetCone();
695 eveJetCone->SetName(Form("jet [%d]", counter++));
696 eveJetCone->SetMainTransparency(60);
697 eveJetCone->SetLineColor(kYellow);
698 eveJetCone->SetCylinder(gRadius - 10, gHalfLength - 10);
699 eveJetCone->SetPickable(kTRUE);
700 eveJetCone->AddEllipticCone(jet->Eta, jet->Phi, jet->DeltaEta, jet->DeltaPhi);
701 gJetList->AddElement(eveJetCone);
702 }
703}
704
705
706/******************************************************************************/
707// GUI
708/******************************************************************************/
709
710//______________________________________________________________________________
711//
712// EvNavHandler class is needed to connect GUI signals.
713
714class EvNavHandler
715{
716public:
717 void Fwd()
718 {
719 if (event_id < gTreeReader->GetEntries() - 1) {
720 ++event_id;
721 load_event();
722 } else {
723 printf("Already at last event.\n");
724 }
725 }
726 void Bck()
727 {
728 if (event_id > 0) {
729 --event_id;
730 load_event();
731 } else {
732 printf("Already at first event.\n");
733 }
734 }
735};
736
737//______________________________________________________________________________
738void make_gui()
739{
740 // Create minimal GUI for event navigation.
741 // TODO: better GUI could be made based on the ch15 of the manual (Writing a GUI)
742
743 // add a tab on the left
744 TEveBrowser* browser = gEve->GetBrowser();
745 browser->StartEmbedding(TRootBrowser::kLeft);
746
747 // set the main title
748 TGMainFrame* frmMain = new TGMainFrame(gClient->GetRoot(), 1000, 600);
749 frmMain->SetWindowName("Delphes Event Display");
750 frmMain->SetCleanup(kDeepCleanup);
751
752 // build the navigation menu
753 TGHorizontalFrame* hf = new TGHorizontalFrame(frmMain);
754 {
755 TString icondir;
756 if(gSystem->Getenv("ROOTSYS"))
757 icondir = Form("%s/icons/", gSystem->Getenv("ROOTSYS"));
758 if(!gSystem->OpenDirectory(icondir))
759 icondir = Form("%s/icons/", (const char*)gSystem->GetFromPipe("root-config --etcdir") );
760 TGPictureButton* b = 0;
761 EvNavHandler *fh = new EvNavHandler;
762
763 b = new TGPictureButton(hf, gClient->GetPicture(icondir+"GoBack.gif"));
764 hf->AddFrame(b);
765 b->Connect("Clicked()", "EvNavHandler", fh, "Bck()");
766
767 b = new TGPictureButton(hf, gClient->GetPicture(icondir+"GoForward.gif"));
768 hf->AddFrame(b);
769 b->Connect("Clicked()", "EvNavHandler", fh, "Fwd()");
770 }
771 frmMain->AddFrame(hf);
772 frmMain->MapSubwindows();
773 frmMain->Resize();
774 frmMain->MapWindow();
775 browser->StopEmbedding();
776 browser->SetTabTitle("Event Control", 0);
777}
778
779/******************************************************************************/
780// MAIN
781/******************************************************************************/
782
783void geometry(const char* filename = "delphes_card_CMS.tcl", const char* ParticlePropagator="ParticlePropagator",
784 const char* TrackingEfficiency="ChargedHadronTrackingEfficiency",
785 const char* MuonEfficiency="MuonEfficiency",
786 const char* Calorimeters="Calorimeter")
787{
788
789 // load the libraries
790 gSystem->Load("libGeom");
791 gSystem->Load("../libDelphes");
792 gSystem->Load("../libDelphesDisplay");
793
794 // create the detector representation
795 Delphes3DGeometry det3D(new TGeoManager("delphes", "Delphes geometry"));
796 det3D.readFile(filename, ParticlePropagator, TrackingEfficiency, MuonEfficiency, Calorimeters);
797
798 // create the application items
799 delphes_event_display("delphes_card_CMS.tcl", "../delphes_output.root", det3D);
800 make_gui();
801 load_event();
802 gEve->Redraw3D(kTRUE); // Reset camera after the first event has been shown.
803
804 // EClipType not exported to CINT (see TGLUtil.h):
805 // 0 - no clip, 1 - clip plane, 2 - clip box
806 TGLViewer *v = gEve->GetDefaultGLViewer();
807 Double_t plane[4] = { 0., 1., 0., 0. };
808 //v->GetClipSet()->SetClipState(1,plane);
809 //v->GetClipSet()->SetClipType(1);
810 //v->ColorSet().Background().SetColor(kMagenta+4);
811 //v->SetGuideState(TGLUtil::kAxesEdge, kTRUE, kFALSE, 0);
812 v->RefreshPadEditor(v);
813 v->CurrentCamera().RotateRad(-1.2, 0.5);
814 v->DoDraw();
815
816}
817
Note: See TracBrowser for help on using the repository browser.