Fork me on GitHub

source: git/examples/geometry.C@ a0b6d15

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

Made the configuration more flexible

The names of modules used in the tcl files are now optionnally set when
loading the cfg. It allows to read the FCC file despites it being quite
different.
Still, for now only one calorimeter can be displayed.
  • Property mode set to 100644
File size: 16.6 KB
Line 
1#include <set>
2#include <map>
3#include <utility>
4#include <vector>
5#include "TGeoManager.h"
6#include "TGeoVolume.h"
7#include "TGeoMedium.h"
8#include "TGeoNode.h"
9#include "TGeoCompositeShape.h"
10#include "TGeoMatrix.h"
11#include "TGeoTube.h"
12#include "TGeoCone.h"
13#include "TGeoArb8.h"
14//#include "../external/ExRootAnalysis/ExRootConfReader.h"
15#include "TF2.h"
16#include "TH1F.h"
17#include "TMath.h"
18#include "TSystem.h"
19
20using namespace std;
21
22// TODO: asymmetric detector
23// TODO: generalize for FCC-like config: >1 calorimeter & flexibility in module names
24class Delphes3DGeometry {
25 public:
26 Delphes3DGeometry(TGeoManager *geom = NULL);
27 ~Delphes3DGeometry() {}
28
29 void readFile(const char* filename, const char* ParticlePropagator="ParticlePropagator",
30 const char* TrackingEfficiency="ChargedHadronTrackingEfficiency",
31 const char* MuonEfficiency="MuonEfficiency",
32 const char* Calorimeters="Calorimeter");
33
34 void setContingency(Double_t contingency) { contingency_ = contingency; }
35 void setCaloBarrelThickness(Double_t thickness) { calo_barrel_thickness_ = thickness; }
36 void setCaloEndcapThickness(Double_t thickness) { calo_endcap_thickness_ = thickness; }
37 void setMuonSystemThickness(Double_t thickness) { muonSystem_thickn_ = thickness; }
38
39 TGeoVolume* getDetector(bool withTowers = true);
40
41 private:
42 void addTracker(TGeoVolume *top);
43 void addCalorimeters(TGeoVolume *top);
44 void addMuonDets(TGeoVolume *top);
45 void addCaloTowers(TGeoVolume *top);
46
47 private:
48
49 TGeoManager *geom_;
50
51 TGeoMedium *vacuum_;
52 TGeoMedium *tkmed_;
53 TGeoMedium *calomed_;
54 TGeoMedium *mudetmed_;
55
56 Double_t contingency_;
57 Double_t calo_barrel_thickness_;
58 Double_t calo_endcap_thickness_;
59 Double_t muonSystem_thickn_;
60 Double_t tk_radius_;
61 Double_t tk_length_;
62 Double_t tk_etamax_;
63 Double_t calo_endcap_etamax_;
64 Double_t muonSystem_etamax_;
65 Double_t calo_barrel_innerRadius_;
66 Double_t calo_endcap_etamin_;
67 Double_t calo_endcap_innerRadius1_;
68 Double_t calo_endcap_innerRadius2_;
69 Double_t calo_endcap_outerRadius1_;
70 Double_t calo_endcap_outerRadius2_;
71 Double_t calo_endcap_coneThickness_;
72 Double_t calo_endcap_diskThickness_;
73
74 set< pair<Double_t, Int_t> > caloBinning_;
75
76};
77
78Delphes3DGeometry::Delphes3DGeometry(TGeoManager *geom) {
79
80 //--- the geometry manager
81 geom_ = geom==NULL? gGeoManager : geom;
82
83 //--- define some materials
84 TGeoMaterial *matVacuum = new TGeoMaterial("Vacuum", 0,0,0);
85 TGeoMaterial *matAl = new TGeoMaterial("Al", 26.98,13,2.7); // placeholder
86
87 //--- define some media
88 TGeoMedium *Vacuum = new TGeoMedium("Vacuum",1, matVacuum);
89 TGeoMedium *Al = new TGeoMedium("Root Material",2, matAl);
90 vacuum_ = Vacuum;
91 tkmed_ = Vacuum; // placeholder
92 calomed_ = Al; // placeholder
93 mudetmed_ = Al; // placeholder
94
95 // custom parameters
96 contingency_ = 10.;
97 calo_barrel_thickness_ = 50.;
98 calo_endcap_thickness_ = 75.;
99 muonSystem_thickn_ = contingency_;
100
101 // read these parameters from the Delphes Card (with default values)
102 tk_radius_ = 120.;
103 tk_length_ = 150.;
104 tk_etamax_ = 3.0;
105 calo_endcap_etamax_ = 2.6;
106 muonSystem_etamax_ = 2.4;
107}
108
109void Delphes3DGeometry::readFile(const char *configFile,
110 const char* ParticlePropagator, const char* TrackingEfficiency,
111 const char* MuonEfficiency, const char* Calorimeters) {
112
113 ExRootConfReader *confReader = new ExRootConfReader;
114 confReader->ReadFile(configFile);
115
116 tk_radius_ = confReader->GetDouble(Form("%s::Radius",ParticlePropagator), 1.0)*100; // tk_radius
117 tk_length_ = confReader->GetDouble(Form("%s::HalfLength",ParticlePropagator), 3.0)*100; // tk_length
118
119 {
120 TString tkEffFormula = confReader->GetString(Form("%s::EfficiencyFormula",TrackingEfficiency),"abs(eta)<3.0");
121 tkEffFormula.ReplaceAll("pt","x");
122 tkEffFormula.ReplaceAll("eta","y");
123 tkEffFormula.ReplaceAll("phi","0.");
124 TF2* tkEffFunction = new TF2("tkEff",tkEffFormula,0,1000,-10,10);
125 TH1F etaHisto("eta","eta",100,5.,-5.);
126 Double_t pt,eta;
127 for(int i=0;i<1000;++i) {
128 tkEffFunction->GetRandom2(pt,eta);
129 etaHisto.Fill(eta);
130 }
131 Int_t bin = -1;
132 bin = etaHisto.FindFirstBinAbove(0.5);
133 Double_t etamin = (bin>-1) ? etaHisto.GetBinLowEdge(bin) : -10.;
134 bin = etaHisto.FindLastBinAbove(0.5);
135 Double_t etamax = (bin>-1) ? etaHisto.GetBinLowEdge(bin+1) : -10.;
136 tk_etamax_ = TMath::Max(fabs(etamin),fabs(etamax)); // tk_etamax
137 delete tkEffFunction;
138 }
139
140 {
141 TString muonEffFormula = confReader->GetString(Form("%s::EfficiencyFormula",MuonEfficiency),"abs(eta)<2.0");
142 muonEffFormula.ReplaceAll("pt","x");
143 muonEffFormula.ReplaceAll("eta","y");
144 muonEffFormula.ReplaceAll("phi","0.");
145 TF2* muEffFunction = new TF2("muEff",muonEffFormula,0,1000,-10,10);
146 TH1F etaHisto("eta2","eta2",100,5.,-5.);
147 Double_t pt,eta;
148 for(int i=0;i<1000;++i) {
149 muEffFunction->GetRandom2(pt,eta);
150 etaHisto.Fill(eta);
151 }
152 Int_t bin = -1;
153 bin = etaHisto.FindFirstBinAbove(0.5);
154 Double_t etamin = (bin>-1) ? etaHisto.GetBinLowEdge(bin) : -10.;
155 bin = etaHisto.FindLastBinAbove(0.5);
156 Double_t etamax = (bin>-1) ? etaHisto.GetBinLowEdge(bin+1) : -10.;
157 muonSystem_etamax_ = TMath::Max(fabs(etamin),fabs(etamax)); // muonSystem_etamax
158 delete muEffFunction;
159 }
160
161 caloBinning_.clear(); // calo binning
162 ExRootConfParam paramEtaBins, paramPhiBins;
163 ExRootConfParam param = confReader->GetParam(Form("%s::EtaPhiBins",Calorimeters));
164 Int_t size = param.GetSize();
165 for(int i = 0; i < size/2; ++i) {
166 paramEtaBins = param[i*2];
167 paramPhiBins = param[i*2+1];
168 assert(paramEtaBins.GetSize()==1);
169 caloBinning_.insert(std::make_pair(paramEtaBins[0].GetDouble(),paramPhiBins.GetSize()-1));
170 }
171
172 if (size>0) calo_endcap_etamax_ = TMath::Max(fabs(caloBinning_.begin()->first),fabs(caloBinning_.rbegin()->first)); // calo_endcap_etamax_
173
174 delete confReader;
175
176 calo_barrel_innerRadius_ = tk_radius_+contingency_;
177 calo_endcap_etamin_ = -log(tk_radius_/(2*tk_length_));
178 calo_endcap_innerRadius1_ = tk_length_*2.*exp(-calo_endcap_etamax_)/(1-exp(-2.*calo_endcap_etamax_));
179 calo_endcap_innerRadius2_ = (tk_length_+calo_endcap_thickness_)*2.*exp(-calo_endcap_etamax_)/(1-exp(-2.*calo_endcap_etamax_));
180 calo_endcap_outerRadius1_ = tk_radius_;
181 calo_endcap_outerRadius2_ = tk_radius_+calo_barrel_thickness_;
182 calo_endcap_coneThickness_ = calo_barrel_thickness_ * (1-exp(-2.*calo_endcap_etamin_)) / (2.*exp(-calo_endcap_etamin_));
183 calo_endcap_diskThickness_ = TMath::Max(0.,calo_endcap_thickness_-calo_endcap_coneThickness_);
184}
185
186TGeoVolume* Delphes3DGeometry::getDetector(bool withTowers) {
187 TGeoVolume *top = geom_->MakeBox("Delphes3DGeometry", vacuum_, 1500, 1500, 2300); // determine the size from what we know about the detector TODO
188 addTracker(top);
189 addCalorimeters(top); // TODO: allow for more than one calo
190 addMuonDets(top);
191 if (withTowers) {
192 addCaloTowers(top);
193 }
194 return top;
195}
196
197//TODO: there should be a cut by two cones to limit the acceptance in eta
198//typically from ChargedHadronTrackingEfficiency
199void Delphes3DGeometry::addTracker(TGeoVolume *top) {
200 // tracker: a cylinder
201 TGeoVolume *tracker = geom_->MakeTube("tracker", tkmed_, 0., tk_radius_, tk_length_);
202 tracker->SetLineColor(kYellow);
203 top->AddNode(tracker,1);
204}
205
206void Delphes3DGeometry::addCalorimeters(TGeoVolume *top) {
207 // calorimeters: tube truncated in eta + cones
208
209 /*TGeoTube *calo_barrel_cylinder =*/ new TGeoTube("calo_barrel_cylinder",calo_barrel_innerRadius_,tk_radius_+calo_barrel_thickness_+contingency_,tk_length_+calo_barrel_thickness_);
210 /*TGeoCone *calo_endcap_cone =*/ new TGeoCone("calo_endcap_cone",calo_endcap_coneThickness_/2.,calo_endcap_innerRadius1_,calo_endcap_outerRadius1_,calo_endcap_innerRadius2_,calo_endcap_outerRadius2_);
211 /*TGeoTube *calo_endcap_disk =*/ new TGeoTube("calo_endcap_disk",calo_endcap_innerRadius2_,tk_radius_+calo_barrel_thickness_,calo_endcap_diskThickness_/2.);
212 TGeoTranslation *tr1 = new TGeoTranslation("tr1",0., 0., (calo_endcap_coneThickness_+calo_endcap_diskThickness_)/2.);
213 tr1->RegisterYourself();
214 TGeoCompositeShape *calo_endcap_cs = new TGeoCompositeShape("calo_endcap_cs","calo_endcap_cone+calo_endcap_disk:tr1");
215 TGeoTranslation *trc1 = new TGeoTranslation("calo_endcap1_position",0.,0., tk_length_+calo_endcap_coneThickness_/2.);
216 trc1->RegisterYourself();
217 TGeoRotation *negz = new TGeoRotation("negz",0,180,0);
218 TGeoCombiTrans *trc2 = new TGeoCombiTrans("calo_endcap2_position",0.,0.,-(tk_length_+calo_endcap_coneThickness_/2.),negz);
219 trc2->RegisterYourself();
220 TGeoTranslation *trc1c = new TGeoTranslation("calo_endcap1_position_cont",0.,0., tk_length_+calo_endcap_coneThickness_/2.+contingency_);
221 trc1c->RegisterYourself();
222 TGeoCombiTrans *trc2c = new TGeoCombiTrans("calo_endcap2_position_cont",0.,0.,-(tk_length_+calo_endcap_coneThickness_/2.)-contingency_,negz);
223 trc2c->RegisterYourself();
224 TGeoVolume *calo_endcap = new TGeoVolume("calo_endcap",calo_endcap_cs,calomed_);
225 TGeoCompositeShape *calo_barrel_cs = new TGeoCompositeShape("calo_barrel_cs","calo_barrel_cylinder-calo_endcap_cs:calo_endcap1_position-calo_endcap_cs:calo_endcap2_position");
226 TGeoVolume *calo_barrel = new TGeoVolume("calo_barrel",calo_barrel_cs,calomed_);
227 calo_endcap->SetLineColor(kViolet);
228 calo_endcap->SetFillColor(kViolet);
229 calo_barrel->SetLineColor(kRed);
230 top->AddNode(calo_endcap,1,trc1c);
231 top->AddNode(calo_endcap,2,trc2c);
232 top->AddNode(calo_barrel,1);
233}
234
235void Delphes3DGeometry::addMuonDets(TGeoVolume *top) {
236
237 // muon system: tube + disks
238 Double_t muonSystem_radius = tk_radius_+calo_barrel_thickness_+2*contingency_;
239 Double_t muonSystem_length = tk_length_+TMath::Max(calo_endcap_coneThickness_,calo_endcap_thickness_)+2*contingency_;
240 Double_t muonSystem_rmin = muonSystem_length*2.*exp(-muonSystem_etamax_)/(1-exp(-2.*muonSystem_etamax_));
241 TGeoVolume *muon_barrel = geom_->MakeTube("muon_barrel",mudetmed_,muonSystem_radius,muonSystem_radius+muonSystem_thickn_,muonSystem_length);
242 muon_barrel->SetLineColor(kBlue);
243 top->AddNode(muon_barrel,1);
244 TGeoVolume *muon_endcap = geom_->MakeTube("muon_endcap",mudetmed_,muonSystem_rmin,muonSystem_radius+muonSystem_thickn_,muonSystem_thickn_/2.);
245 muon_endcap->SetLineColor(kBlue);
246 TGeoTranslation *trm1 = new TGeoTranslation("muonEndcap1_position",0.,0.,muonSystem_length);
247 trm1->RegisterYourself();
248 TGeoTranslation *trm2 = new TGeoTranslation("muonEndcap2_position",0.,0.,-muonSystem_length);
249 trm1->RegisterYourself();
250 top->AddNode(muon_endcap,1,trm1);
251 top->AddNode(muon_endcap,1,trm2);
252}
253
254void Delphes3DGeometry::addCaloTowers(TGeoVolume *top) {
255
256 TGeoVolume* calo_endcap = top->GetNode("calo_endcap_1")->GetVolume();
257 TGeoVolume* calo_barrel = top->GetNode("calo_barrel_1")->GetVolume();
258
259 // calo towers in the barrel
260 Double_t vertices[16] = {0.,0.,0.,0.,0.,0.,0.,0.}; // summit of the pyramid
261 Double_t R = tk_radius_+calo_barrel_thickness_+2*contingency_; // radius of the muons system = height of the pyramid
262 Int_t nEtaBins = caloBinning_.size();
263 // this rotation is to make the tower point "up"
264 TGeoRotation* initTowerRot = new TGeoRotation("initTowerRot",0.,90.,0.);
265 TGeoCombiTrans* initTower = new TGeoCombiTrans("initTower",0.,-R/2.,0.,initTowerRot);
266 initTower->RegisterYourself();
267 // eta bins... we build one pyramid per eta slice and then translate it nphi times.
268 // phi bins represented by rotations around z
269 Double_t *y = new Double_t[nEtaBins];
270 Double_t *dx = new Double_t[nEtaBins];
271 Int_t *nphi = new Int_t[nEtaBins];
272 Int_t etaslice = 0;
273 std::map<std::pair<int,int>, TGeoRotation*> phirotations;
274 for(set< pair<Double_t, Int_t> >::const_iterator bin=caloBinning_.begin(); bin!=caloBinning_.end();++bin) {
275 if(abs(bin->first)>calo_endcap_etamin_) continue; // only in the barrel
276 nphi[etaslice] = bin->second;
277 y[etaslice] = 0.5*R*(1-exp(-2*bin->first))/exp(-bin->first);
278 Double_t phiRotationAngle = 360./nphi[etaslice];
279 dx[etaslice] = R*tan(TMath::Pi()*phiRotationAngle/360.);
280 for(int phislice=0;phislice<nphi[etaslice];++phislice) {
281 phirotations[make_pair(etaslice,phislice)] = new TGeoRotation(Form("phi%d_%d",etaslice,phislice),phiRotationAngle*phislice,0.,0.);
282 phirotations[make_pair(etaslice,phislice)]->RegisterYourself();
283 }
284 ++etaslice;
285 }
286 nEtaBins = ++etaslice;
287 for(int i=0;i<nEtaBins-1;++i) { // loop on the eta slices
288 vertices[8] = -dx[i]; vertices[9] = y[i];
289 vertices[10] = -dx[i]; vertices[11] = y[i+1];
290 vertices[12] = dx[i]; vertices[13] = y[i+1];
291 vertices[14] = dx[i]; vertices[15] = y[i];
292 /*TGeoArb8 *tower =*/ new TGeoArb8(Form("tower%d",i),R/2., vertices); // tower in the proper eta slice, at phi=0
293 // intersection between the tower and the calo_barrel
294 TGeoCompositeShape *finaltower_cs = new TGeoCompositeShape(Form("ftower%d_cs",i),Form("tower%d:initTower*calo_barrel_cs",i));
295 TGeoVolume *finaltower = new TGeoVolume(Form("ftower%d",i),finaltower_cs,calomed_);
296 finaltower->SetLineColor(kRed);
297 for(int j=0;j<nphi[i];++j) { // loop on the phi slices
298 calo_barrel->AddNode(finaltower,j,phirotations[make_pair(i,j)]);
299 }
300 }
301 delete[] y;
302 delete[] dx;
303 delete[] nphi;
304 //the towers in the forward region
305 R = tk_length_+calo_endcap_thickness_+3*contingency_; // Z of the muons system = height of the pyramid
306 nEtaBins = caloBinning_.size();
307 // translation to bring the origin of the tower to (0,0,0)
308 TGeoTranslation* towerdz = new TGeoTranslation("towerdz",0.,0.,R/2.-(tk_length_+calo_endcap_coneThickness_/2.));
309 towerdz->RegisterYourself();
310 // eta bins... we build one pyramid per eta slice and then translate it nphi times.
311 Double_t *r = new Double_t[nEtaBins];
312 nphi = new Int_t[nEtaBins];
313 etaslice = 0;
314 phirotations.clear();
315 for(set< pair<Double_t, Int_t> >::const_iterator bin=caloBinning_.begin(); bin!=caloBinning_.end();++bin) {
316 if(bin->first<calo_endcap_etamin_) continue; // only in the + endcap
317 r[etaslice] = R*2*exp(-bin->first)/(1-exp(-2*bin->first));
318 nphi[etaslice] = bin->second;
319 Double_t phiRotationAngle = 360./nphi[etaslice];
320 for(int phislice=0;phislice<nphi[etaslice];++phislice) {
321 phirotations[make_pair(etaslice,phislice)] = new TGeoRotation(Form("forward_phi%d_%d",etaslice,phislice),phiRotationAngle*phislice,0.,0.);
322 phirotations[make_pair(etaslice,phislice)]->RegisterYourself();
323 }
324 ++etaslice;
325 }
326 nEtaBins = etaslice;
327 for(int i=0;i<nEtaBins;++i) { // loop on the eta slices
328 vertices[8] = -r[i+1]*sin(TMath::Pi()/20.); vertices[9] = r[i+1]*cos(TMath::Pi()/20.);
329 vertices[10] = -r[i]*sin(TMath::Pi()/20.); vertices[11] = r[i]*cos(TMath::Pi()/20.);
330 vertices[12] = r[i]*sin(TMath::Pi()/20.); vertices[13] = r[i]*cos(TMath::Pi()/20.);
331 vertices[14] = r[i+1]*sin(TMath::Pi()/20.); vertices[15] = r[i+1]*cos(TMath::Pi()/20.);
332 /*TGeoArb8 *fwdtower =*/ new TGeoArb8(Form("fwdtower%d",i),R/2., vertices); // tower in the proper eta slice, at phi=0
333 // intersection between the tower and the calo_endcap
334 TGeoCompositeShape *finalfwdtower_cs = new TGeoCompositeShape(Form("ffwdtower%d_cs",i),Form("fwdtower%d:towerdz*calo_endcap_cs",i));
335 TGeoVolume *finalfwdtower = new TGeoVolume(Form("ffwdtower%d",i),finalfwdtower_cs,calomed_);
336 finalfwdtower->SetLineColor(kViolet);
337 for(int j=0;j<nphi[i];++j) { // loop on the phi slices
338 calo_endcap->AddNode(finalfwdtower,j,phirotations[make_pair(i,j)]);
339 }
340 }
341 delete[] r;
342 delete[] nphi;
343}
344
345void geometry(const char* filename = "delphes_card_CMS.tcl", const char* ParticlePropagator="ParticlePropagator",
346 const char* TrackingEfficiency="ChargedHadronTrackingEfficiency",
347 const char* MuonEfficiency="MuonEfficiency",
348 const char* Calorimeters="Calorimeter")
349{
350 gSystem->Load("libGeom");
351 gSystem->Load("../libDelphes");
352 TGeoManager *geom = new TGeoManager("delphes", "Delphes geometry");
353
354 // make the top container volume -> designed to contain a "big" detector (ATLAS)
355 TGeoVolume *top = geom->MakeBox("TOP", 0, 1500, 1500, 2300);
356 geom->SetTopVolume(top);
357
358 // build the detector
359 Delphes3DGeometry det3D;
360 det3D.readFile(filename,ParticlePropagator, TrackingEfficiency, MuonEfficiency, Calorimeters);
361 top->AddNode(det3D.getDetector(true),1);
362
363 // draw it
364 geom->CloseGeometry();
365 top->Draw();
366}
367
Note: See TracBrowser for help on using the repository browser.