[bf6ed57] | 1 | #include <set>
|
---|
| 2 | #include <map>
|
---|
| 3 | #include <utility>
|
---|
| 4 | #include <vector>
|
---|
[df634866] | 5 | #include <algorithm>
|
---|
[d2ac201] | 6 | #include <sstream>
|
---|
[bf6ed57] | 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 "TF2.h"
|
---|
| 18 | #include "TH1F.h"
|
---|
| 19 | #include "TMath.h"
|
---|
| 20 | #include "TSystem.h"
|
---|
| 21 |
|
---|
| 22 | using namespace std;
|
---|
| 23 |
|
---|
| 24 | // TODO: asymmetric detector
|
---|
| 25 | class Delphes3DGeometry {
|
---|
| 26 | public:
|
---|
| 27 | Delphes3DGeometry(TGeoManager *geom = NULL);
|
---|
| 28 | ~Delphes3DGeometry() {}
|
---|
| 29 |
|
---|
[a0b6d15] | 30 | void readFile(const char* filename, const char* ParticlePropagator="ParticlePropagator",
|
---|
| 31 | const char* TrackingEfficiency="ChargedHadronTrackingEfficiency",
|
---|
| 32 | const char* MuonEfficiency="MuonEfficiency",
|
---|
| 33 | const char* Calorimeters="Calorimeter");
|
---|
[bf6ed57] | 34 |
|
---|
| 35 | void setContingency(Double_t contingency) { contingency_ = contingency; }
|
---|
| 36 | void setCaloBarrelThickness(Double_t thickness) { calo_barrel_thickness_ = thickness; }
|
---|
| 37 | void setCaloEndcapThickness(Double_t thickness) { calo_endcap_thickness_ = thickness; }
|
---|
| 38 | void setMuonSystemThickness(Double_t thickness) { muonSystem_thickn_ = thickness; }
|
---|
| 39 |
|
---|
| 40 | TGeoVolume* getDetector(bool withTowers = true);
|
---|
| 41 |
|
---|
| 42 | private:
|
---|
[df634866] | 43 | std::pair<Double_t, Double_t> addTracker(TGeoVolume *top);
|
---|
| 44 | 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);
|
---|
| 45 | std::pair<Double_t, Double_t> addMuonDets(TGeoVolume *top, const char* name, Double_t innerBarrelRadius, Double_t innerBarrelLength);
|
---|
| 46 | void addCaloTowers(TGeoVolume *top, const char* name, Double_t innerBarrelRadius, Double_t innerBarrelLength, set< pair<Double_t, Int_t> >& caloBinning);
|
---|
[bf6ed57] | 47 |
|
---|
| 48 | private:
|
---|
| 49 |
|
---|
| 50 | TGeoManager *geom_;
|
---|
| 51 |
|
---|
| 52 | TGeoMedium *vacuum_;
|
---|
| 53 | TGeoMedium *tkmed_;
|
---|
| 54 | TGeoMedium *calomed_;
|
---|
| 55 | TGeoMedium *mudetmed_;
|
---|
| 56 |
|
---|
| 57 | Double_t contingency_;
|
---|
| 58 | Double_t calo_barrel_thickness_;
|
---|
| 59 | Double_t calo_endcap_thickness_;
|
---|
| 60 | Double_t muonSystem_thickn_;
|
---|
| 61 | Double_t tk_radius_;
|
---|
| 62 | Double_t tk_length_;
|
---|
[a0b6d15] | 63 | Double_t tk_etamax_;
|
---|
[df634866] | 64 |
|
---|
| 65 | std::vector<std::string> calorimeters_;
|
---|
| 66 | std::vector<std::string> muondets_;
|
---|
| 67 |
|
---|
| 68 | std::map<std::string, Double_t> muonSystem_etamax_;
|
---|
| 69 | std::map<std::string, set< pair<Double_t, Int_t> > > caloBinning_;
|
---|
[bf6ed57] | 70 |
|
---|
| 71 | };
|
---|
| 72 |
|
---|
| 73 | Delphes3DGeometry::Delphes3DGeometry(TGeoManager *geom) {
|
---|
| 74 |
|
---|
| 75 | //--- the geometry manager
|
---|
| 76 | geom_ = geom==NULL? gGeoManager : geom;
|
---|
| 77 |
|
---|
| 78 | //--- define some materials
|
---|
| 79 | TGeoMaterial *matVacuum = new TGeoMaterial("Vacuum", 0,0,0);
|
---|
| 80 | TGeoMaterial *matAl = new TGeoMaterial("Al", 26.98,13,2.7); // placeholder
|
---|
| 81 |
|
---|
| 82 | //--- define some media
|
---|
| 83 | TGeoMedium *Vacuum = new TGeoMedium("Vacuum",1, matVacuum);
|
---|
| 84 | TGeoMedium *Al = new TGeoMedium("Root Material",2, matAl);
|
---|
| 85 | vacuum_ = Vacuum;
|
---|
| 86 | tkmed_ = Vacuum; // placeholder
|
---|
| 87 | calomed_ = Al; // placeholder
|
---|
| 88 | mudetmed_ = Al; // placeholder
|
---|
| 89 |
|
---|
| 90 | // custom parameters
|
---|
| 91 | contingency_ = 10.;
|
---|
| 92 | calo_barrel_thickness_ = 50.;
|
---|
| 93 | calo_endcap_thickness_ = 75.;
|
---|
[df634866] | 94 | muonSystem_thickn_ = 10.;
|
---|
[bf6ed57] | 95 |
|
---|
| 96 | // read these parameters from the Delphes Card (with default values)
|
---|
| 97 | tk_radius_ = 120.;
|
---|
| 98 | tk_length_ = 150.;
|
---|
[a0b6d15] | 99 | tk_etamax_ = 3.0;
|
---|
[bf6ed57] | 100 | }
|
---|
| 101 |
|
---|
[a0b6d15] | 102 | void Delphes3DGeometry::readFile(const char *configFile,
|
---|
| 103 | const char* ParticlePropagator, const char* TrackingEfficiency,
|
---|
| 104 | const char* MuonEfficiency, const char* Calorimeters) {
|
---|
[bf6ed57] | 105 |
|
---|
| 106 | ExRootConfReader *confReader = new ExRootConfReader;
|
---|
| 107 | confReader->ReadFile(configFile);
|
---|
| 108 |
|
---|
[a0b6d15] | 109 | tk_radius_ = confReader->GetDouble(Form("%s::Radius",ParticlePropagator), 1.0)*100; // tk_radius
|
---|
| 110 | tk_length_ = confReader->GetDouble(Form("%s::HalfLength",ParticlePropagator), 3.0)*100; // tk_length
|
---|
[bf6ed57] | 111 |
|
---|
[a0b6d15] | 112 | {
|
---|
| 113 | TString tkEffFormula = confReader->GetString(Form("%s::EfficiencyFormula",TrackingEfficiency),"abs(eta)<3.0");
|
---|
| 114 | tkEffFormula.ReplaceAll("pt","x");
|
---|
| 115 | tkEffFormula.ReplaceAll("eta","y");
|
---|
| 116 | tkEffFormula.ReplaceAll("phi","0.");
|
---|
| 117 | TF2* tkEffFunction = new TF2("tkEff",tkEffFormula,0,1000,-10,10);
|
---|
| 118 | TH1F etaHisto("eta","eta",100,5.,-5.);
|
---|
| 119 | Double_t pt,eta;
|
---|
| 120 | for(int i=0;i<1000;++i) {
|
---|
| 121 | tkEffFunction->GetRandom2(pt,eta);
|
---|
| 122 | etaHisto.Fill(eta);
|
---|
| 123 | }
|
---|
| 124 | Int_t bin = -1;
|
---|
| 125 | bin = etaHisto.FindFirstBinAbove(0.5);
|
---|
| 126 | Double_t etamin = (bin>-1) ? etaHisto.GetBinLowEdge(bin) : -10.;
|
---|
| 127 | bin = etaHisto.FindLastBinAbove(0.5);
|
---|
| 128 | Double_t etamax = (bin>-1) ? etaHisto.GetBinLowEdge(bin+1) : -10.;
|
---|
| 129 | tk_etamax_ = TMath::Max(fabs(etamin),fabs(etamax)); // tk_etamax
|
---|
| 130 | delete tkEffFunction;
|
---|
| 131 | }
|
---|
| 132 |
|
---|
| 133 | {
|
---|
[df634866] | 134 | muondets_.push_back("muons");
|
---|
[a0b6d15] | 135 | TString muonEffFormula = confReader->GetString(Form("%s::EfficiencyFormula",MuonEfficiency),"abs(eta)<2.0");
|
---|
[bf6ed57] | 136 | muonEffFormula.ReplaceAll("pt","x");
|
---|
| 137 | muonEffFormula.ReplaceAll("eta","y");
|
---|
| 138 | muonEffFormula.ReplaceAll("phi","0.");
|
---|
| 139 | TF2* muEffFunction = new TF2("muEff",muonEffFormula,0,1000,-10,10);
|
---|
[a0b6d15] | 140 | TH1F etaHisto("eta2","eta2",100,5.,-5.);
|
---|
[bf6ed57] | 141 | Double_t pt,eta;
|
---|
| 142 | for(int i=0;i<1000;++i) {
|
---|
| 143 | muEffFunction->GetRandom2(pt,eta);
|
---|
| 144 | etaHisto.Fill(eta);
|
---|
| 145 | }
|
---|
| 146 | Int_t bin = -1;
|
---|
| 147 | bin = etaHisto.FindFirstBinAbove(0.5);
|
---|
| 148 | Double_t etamin = (bin>-1) ? etaHisto.GetBinLowEdge(bin) : -10.;
|
---|
| 149 | bin = etaHisto.FindLastBinAbove(0.5);
|
---|
| 150 | Double_t etamax = (bin>-1) ? etaHisto.GetBinLowEdge(bin+1) : -10.;
|
---|
[df634866] | 151 | muonSystem_etamax_["muons"] = TMath::Max(fabs(etamin),fabs(etamax)); // muonSystem_etamax
|
---|
[bf6ed57] | 152 | delete muEffFunction;
|
---|
[a0b6d15] | 153 | }
|
---|
[df634866] | 154 |
|
---|
| 155 | std::string s(Calorimeters);
|
---|
| 156 | std::replace( s.begin(), s.end(), ',', ' ' );
|
---|
| 157 | std::istringstream stream( s );
|
---|
| 158 | std::string word;
|
---|
| 159 | while (stream >> word) calorimeters_.push_back(word);
|
---|
| 160 |
|
---|
[bf6ed57] | 161 | caloBinning_.clear(); // calo binning
|
---|
[df634866] | 162 | for(std::vector<std::string>::const_iterator calo=calorimeters_.begin();calo!=calorimeters_.end(); ++calo) {
|
---|
| 163 | set< pair<Double_t, Int_t> > caloBinning;
|
---|
| 164 | ExRootConfParam paramEtaBins, paramPhiBins;
|
---|
| 165 | ExRootConfParam param = confReader->GetParam(Form("%s::EtaPhiBins",calo->c_str()));
|
---|
| 166 | Int_t size = param.GetSize();
|
---|
| 167 | for(int i = 0; i < size/2; ++i) {
|
---|
| 168 | paramEtaBins = param[i*2];
|
---|
| 169 | paramPhiBins = param[i*2+1];
|
---|
| 170 | assert(paramEtaBins.GetSize()==1);
|
---|
| 171 | caloBinning.insert(std::make_pair(paramEtaBins[0].GetDouble(),paramPhiBins.GetSize()-1));
|
---|
| 172 | }
|
---|
| 173 | caloBinning_[*calo] = caloBinning;
|
---|
[bf6ed57] | 174 | }
|
---|
| 175 |
|
---|
| 176 | delete confReader;
|
---|
| 177 |
|
---|
| 178 | }
|
---|
| 179 |
|
---|
| 180 | TGeoVolume* Delphes3DGeometry::getDetector(bool withTowers) {
|
---|
[1d42548] | 181 | // compute the envelope
|
---|
| 182 | Double_t system_radius = tk_radius_+calo_barrel_thickness_+3*contingency_;
|
---|
[df634866] | 183 | Double_t system_length = tk_length_+contingency_+(contingency_+calo_endcap_thickness_)*calorimeters_.size()+contingency_;
|
---|
[1d42548] | 184 | // the detector volume
|
---|
| 185 | TGeoVolume *top = geom_->MakeBox("Delphes3DGeometry", vacuum_, system_radius, system_radius, system_length);
|
---|
| 186 | // build the detector
|
---|
[df634866] | 187 | std::pair<Double_t, Double_t> limits = addTracker(top);
|
---|
| 188 | Double_t radius = limits.first;
|
---|
| 189 | Double_t length = limits.second;
|
---|
| 190 | for(std::vector<std::string>::const_iterator calo = calorimeters_.begin(); calo != calorimeters_.end(); ++calo) {
|
---|
| 191 | limits = addCalorimeter(top,calo->c_str(),radius,length,caloBinning_[*calo]);
|
---|
| 192 | if (withTowers) {
|
---|
| 193 | addCaloTowers(top,calo->c_str(),radius,length,caloBinning_[*calo]);
|
---|
| 194 | }
|
---|
| 195 | radius = limits.first;
|
---|
| 196 | length = limits.second;
|
---|
| 197 | }
|
---|
| 198 | for(std::vector<std::string>::const_iterator muon = muondets_.begin(); muon != muondets_.end(); ++muon) {
|
---|
| 199 | limits = addMuonDets(top, muon->c_str(), radius, length);
|
---|
| 200 | radius = limits.first;
|
---|
| 201 | length = limits.second;
|
---|
[bf6ed57] | 202 | }
|
---|
[1d42548] | 203 | // return the result
|
---|
[bf6ed57] | 204 | return top;
|
---|
| 205 | }
|
---|
| 206 |
|
---|
[df634866] | 207 | std::pair<Double_t, Double_t> Delphes3DGeometry::addTracker(TGeoVolume *top) {
|
---|
[1d42548] | 208 | // tracker: a cylinder with two cones substracted
|
---|
[d2ac201] | 209 | new TGeoCone("forwardTkAcceptance",(tk_length_/2.+0.05),0.,tk_radius_,(tk_length_)*2.*exp(-tk_etamax_)/(1-exp(-2.*tk_etamax_)),tk_radius_);
|
---|
[1d42548] | 210 | TGeoTranslation *tr1 = new TGeoTranslation("tkacc1",0., 0., tk_length_/2.);
|
---|
| 211 | tr1->RegisterYourself();
|
---|
| 212 | TGeoRotation *negz = new TGeoRotation("tknegz",0,180,0);
|
---|
| 213 | negz->RegisterYourself();
|
---|
| 214 | TGeoCombiTrans *tr2 = new TGeoCombiTrans("tkacc2",0.,0.,-tk_length_/2.,negz);
|
---|
| 215 | tr2->RegisterYourself();
|
---|
| 216 | TGeoCompositeShape* tracker_cs = new TGeoCompositeShape("tracker_cs","forwardTkAcceptance:tkacc1+forwardTkAcceptance:tkacc2");
|
---|
| 217 | TGeoVolume *tracker = new TGeoVolume("tracker",tracker_cs,tkmed_);
|
---|
[bf6ed57] | 218 | tracker->SetLineColor(kYellow);
|
---|
| 219 | top->AddNode(tracker,1);
|
---|
[df634866] | 220 | return std::make_pair(tk_radius_,tk_length_);
|
---|
[bf6ed57] | 221 | }
|
---|
| 222 |
|
---|
[df634866] | 223 | std::pair<Double_t, Double_t> Delphes3DGeometry::addCalorimeter(TGeoVolume *top, const char* name,
|
---|
| 224 | Double_t innerBarrelRadius, Double_t innerBarrelLength, set< pair<Double_t, Int_t> >& caloBinning) {
|
---|
| 225 | // parameters derived from the inputs
|
---|
| 226 | Double_t calo_endcap_etamax = TMath::Max(fabs(caloBinning.begin()->first),fabs(caloBinning.rbegin()->first));
|
---|
| 227 | Double_t calo_barrel_innerRadius = innerBarrelRadius+contingency_;
|
---|
| 228 | Double_t calo_barrel_length = innerBarrelLength + calo_barrel_thickness_;
|
---|
| 229 | Double_t calo_endcap_etamin = -log(innerBarrelRadius/(2*innerBarrelLength));
|
---|
| 230 | Double_t calo_endcap_innerRadius1 = innerBarrelLength*2.*exp(-calo_endcap_etamax)/(1-exp(-2.*calo_endcap_etamax));
|
---|
| 231 | Double_t calo_endcap_innerRadius2 = (innerBarrelLength+calo_endcap_thickness_)*2.*exp(-calo_endcap_etamax)/(1-exp(-2.*calo_endcap_etamax));
|
---|
| 232 | Double_t calo_endcap_outerRadius1 = innerBarrelRadius;
|
---|
| 233 | Double_t calo_endcap_outerRadius2 = innerBarrelRadius+calo_barrel_thickness_;
|
---|
| 234 | Double_t calo_endcap_coneThickness = TMath::Min(calo_barrel_thickness_ * (1-exp(-2.*calo_endcap_etamin)) / (2.*exp(-calo_endcap_etamin)), calo_endcap_thickness_);
|
---|
| 235 | Double_t calo_endcap_diskThickness = TMath::Max(0.,calo_endcap_thickness_-calo_endcap_coneThickness);
|
---|
| 236 |
|
---|
[bf6ed57] | 237 | // calorimeters: tube truncated in eta + cones
|
---|
[df634866] | 238 | new TGeoTube(Form("%s_barrel_cylinder",name),calo_barrel_innerRadius,calo_barrel_innerRadius+calo_barrel_thickness_,calo_barrel_length);
|
---|
| 239 | new TGeoCone(Form("%s_endcap_cone",name),calo_endcap_coneThickness/2.,calo_endcap_innerRadius1,calo_endcap_outerRadius1,calo_endcap_innerRadius2,calo_endcap_outerRadius2);
|
---|
| 240 | new TGeoTube(Form("%s_endcap_disk",name),calo_endcap_innerRadius2,tk_radius_+calo_barrel_thickness_,calo_endcap_diskThickness/2.);
|
---|
| 241 | TGeoTranslation *tr1 = new TGeoTranslation(Form("%s_tr1",name),0., 0., (calo_endcap_coneThickness+calo_endcap_diskThickness)/2.);
|
---|
[bf6ed57] | 242 | tr1->RegisterYourself();
|
---|
[df634866] | 243 | TGeoCompositeShape *calo_endcap_cs = new TGeoCompositeShape(Form("%s_endcap_cs",name),Form("%s_endcap_cone+%s_endcap_disk:%s_tr1",name,name,name));
|
---|
| 244 | TGeoTranslation *trc1 = new TGeoTranslation(Form("%s_endcap1_position",name),0.,0., innerBarrelLength+calo_endcap_coneThickness/2.);
|
---|
[bf6ed57] | 245 | trc1->RegisterYourself();
|
---|
[df634866] | 246 | TGeoRotation *negz = new TGeoRotation(Form("%s_negz",name),0,180,0);
|
---|
| 247 | TGeoCombiTrans *trc2 = new TGeoCombiTrans(Form("%s_endcap2_position",name),0.,0.,-(innerBarrelLength+calo_endcap_coneThickness/2.),negz);
|
---|
[bf6ed57] | 248 | trc2->RegisterYourself();
|
---|
[df634866] | 249 | TGeoTranslation *trc1c = new TGeoTranslation(Form("%s_endcap1_position_cont",name),0.,0., innerBarrelLength+calo_endcap_coneThickness/2.+contingency_);
|
---|
[bf6ed57] | 250 | trc1c->RegisterYourself();
|
---|
[df634866] | 251 | TGeoCombiTrans *trc2c = new TGeoCombiTrans(Form("%s_endcap2_position_cont",name),0.,0.,-(innerBarrelLength+calo_endcap_coneThickness/2.)-contingency_,negz);
|
---|
[bf6ed57] | 252 | trc2c->RegisterYourself();
|
---|
[df634866] | 253 | TGeoVolume *calo_endcap = new TGeoVolume(Form("%s_endcap",name),calo_endcap_cs,calomed_);
|
---|
| 254 | TGeoCompositeShape *calo_barrel_cs = new TGeoCompositeShape(Form("%s_barrel_cs",name),
|
---|
| 255 | Form("%s_barrel_cylinder-%s_endcap_cs:%s_endcap1_position-%s_endcap_cs:%s_endcap2_position",name,name,name,name,name));
|
---|
| 256 | TGeoVolume *calo_barrel = new TGeoVolume(Form("%s_barrel",name),calo_barrel_cs,calomed_);
|
---|
[bf6ed57] | 257 | calo_endcap->SetLineColor(kViolet);
|
---|
| 258 | calo_endcap->SetFillColor(kViolet);
|
---|
| 259 | calo_barrel->SetLineColor(kRed);
|
---|
| 260 | top->AddNode(calo_endcap,1,trc1c);
|
---|
| 261 | top->AddNode(calo_endcap,2,trc2c);
|
---|
| 262 | top->AddNode(calo_barrel,1);
|
---|
[df634866] | 263 | return std::make_pair(calo_barrel_innerRadius+calo_barrel_thickness_,innerBarrelLength+calo_endcap_thickness_+contingency_);
|
---|
[bf6ed57] | 264 | }
|
---|
| 265 |
|
---|
[df634866] | 266 | std::pair<Double_t, Double_t> Delphes3DGeometry::addMuonDets(TGeoVolume *top, const char* name, Double_t innerBarrelRadius, Double_t innerBarrelLength) {
|
---|
[bf6ed57] | 267 | // muon system: tube + disks
|
---|
[df634866] | 268 | Double_t muonSystem_radius = innerBarrelRadius + contingency_;
|
---|
| 269 | Double_t muonSystem_length = innerBarrelLength + contingency_;
|
---|
| 270 | Double_t muonSystem_rmin = muonSystem_length*2.*exp(-muonSystem_etamax_[name])/(1-exp(-2.*muonSystem_etamax_[name]));
|
---|
| 271 | TGeoVolume *muon_barrel = geom_->MakeTube(Form("%s_barrel",name),mudetmed_,muonSystem_radius,muonSystem_radius+muonSystem_thickn_,muonSystem_length);
|
---|
[bf6ed57] | 272 | muon_barrel->SetLineColor(kBlue);
|
---|
| 273 | top->AddNode(muon_barrel,1);
|
---|
[df634866] | 274 | TGeoVolume *muon_endcap = geom_->MakeTube(Form("%s_endcap",name),mudetmed_,muonSystem_rmin,muonSystem_radius+muonSystem_thickn_,muonSystem_thickn_/2.);
|
---|
[bf6ed57] | 275 | muon_endcap->SetLineColor(kBlue);
|
---|
[df634866] | 276 | TGeoTranslation *trm1 = new TGeoTranslation(Form("%sEndcap1_position",name),0.,0.,muonSystem_length);
|
---|
[bf6ed57] | 277 | trm1->RegisterYourself();
|
---|
[df634866] | 278 | TGeoTranslation *trm2 = new TGeoTranslation(Form("%sEndcap2_position",name),0.,0.,-muonSystem_length);
|
---|
[bf6ed57] | 279 | trm1->RegisterYourself();
|
---|
| 280 | top->AddNode(muon_endcap,1,trm1);
|
---|
[164f032] | 281 | top->AddNode(muon_endcap,2,trm2);
|
---|
[df634866] | 282 | return std::make_pair(muonSystem_radius,muonSystem_length);
|
---|
[bf6ed57] | 283 | }
|
---|
| 284 |
|
---|
[df634866] | 285 | void Delphes3DGeometry::addCaloTowers(TGeoVolume *top, const char* name,
|
---|
| 286 | Double_t innerBarrelRadius, Double_t innerBarrelLength, set< pair<Double_t, Int_t> >& caloBinning) {
|
---|
[bf6ed57] | 287 |
|
---|
[df634866] | 288 | TGeoVolume* calo_endcap = top->GetNode(Form("%s_endcap_1",name))->GetVolume();
|
---|
| 289 | TGeoVolume* calo_barrel = top->GetNode(Form("%s_barrel_1",name))->GetVolume();
|
---|
| 290 | Double_t calo_endcap_etamin = -log(innerBarrelRadius/(2*innerBarrelLength));
|
---|
| 291 | Double_t calo_endcap_coneThickness = TMath::Min(calo_barrel_thickness_ * (1-exp(-2.*calo_endcap_etamin)) / (2.*exp(-calo_endcap_etamin)), calo_endcap_thickness_);
|
---|
[bf6ed57] | 292 |
|
---|
| 293 | // calo towers in the barrel
|
---|
| 294 | Double_t vertices[16] = {0.,0.,0.,0.,0.,0.,0.,0.}; // summit of the pyramid
|
---|
[df634866] | 295 | Double_t R = tk_radius_ + contingency_+(contingency_+calo_barrel_thickness_)*calorimeters_.size(); // radius of the muons system = height of the pyramid
|
---|
| 296 | Int_t nEtaBins = caloBinning.size();
|
---|
[bf6ed57] | 297 | // this rotation is to make the tower point "up"
|
---|
[df634866] | 298 | TGeoRotation* initTowerRot = new TGeoRotation(Form("%s_initTowerRot",name),0.,90.,0.);
|
---|
| 299 | TGeoCombiTrans* initTower = new TGeoCombiTrans(Form("%s_initTower",name),0.,-R/2.,0.,initTowerRot);
|
---|
[bf6ed57] | 300 | initTower->RegisterYourself();
|
---|
| 301 | // eta bins... we build one pyramid per eta slice and then translate it nphi times.
|
---|
| 302 | // phi bins represented by rotations around z
|
---|
| 303 | Double_t *y = new Double_t[nEtaBins];
|
---|
| 304 | Double_t *dx = new Double_t[nEtaBins];
|
---|
| 305 | Int_t *nphi = new Int_t[nEtaBins];
|
---|
| 306 | Int_t etaslice = 0;
|
---|
| 307 | std::map<std::pair<int,int>, TGeoRotation*> phirotations;
|
---|
[df634866] | 308 | for(set< pair<Double_t, Int_t> >::const_iterator bin=caloBinning.begin(); bin!=caloBinning.end();++bin) {
|
---|
| 309 | if(abs(bin->first)>calo_endcap_etamin) continue; // only in the barrel
|
---|
[bf6ed57] | 310 | nphi[etaslice] = bin->second;
|
---|
| 311 | y[etaslice] = 0.5*R*(1-exp(-2*bin->first))/exp(-bin->first);
|
---|
| 312 | Double_t phiRotationAngle = 360./nphi[etaslice];
|
---|
| 313 | dx[etaslice] = R*tan(TMath::Pi()*phiRotationAngle/360.);
|
---|
| 314 | for(int phislice=0;phislice<nphi[etaslice];++phislice) {
|
---|
[df634866] | 315 | phirotations[make_pair(etaslice,phislice)] = new TGeoRotation(Form("%s_phi%d_%d",name,etaslice,phislice),phiRotationAngle*phislice,0.,0.);
|
---|
[bf6ed57] | 316 | phirotations[make_pair(etaslice,phislice)]->RegisterYourself();
|
---|
| 317 | }
|
---|
| 318 | ++etaslice;
|
---|
| 319 | }
|
---|
[d2ac201] | 320 | nEtaBins = etaslice;
|
---|
[bf6ed57] | 321 | for(int i=0;i<nEtaBins-1;++i) { // loop on the eta slices
|
---|
| 322 | vertices[8] = -dx[i]; vertices[9] = y[i];
|
---|
| 323 | vertices[10] = -dx[i]; vertices[11] = y[i+1];
|
---|
| 324 | vertices[12] = dx[i]; vertices[13] = y[i+1];
|
---|
| 325 | vertices[14] = dx[i]; vertices[15] = y[i];
|
---|
[df634866] | 326 | new TGeoArb8(Form("%s_tower%d",name,i),R/2., vertices); // tower in the proper eta slice, at phi=0
|
---|
[bf6ed57] | 327 | // intersection between the tower and the calo_barrel
|
---|
[df634866] | 328 | 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));
|
---|
| 329 | TGeoVolume *finaltower = new TGeoVolume(Form("%s_ftower%d",name,i),finaltower_cs,calomed_);
|
---|
[bf6ed57] | 330 | finaltower->SetLineColor(kRed);
|
---|
| 331 | for(int j=0;j<nphi[i];++j) { // loop on the phi slices
|
---|
| 332 | calo_barrel->AddNode(finaltower,j,phirotations[make_pair(i,j)]);
|
---|
| 333 | }
|
---|
| 334 | }
|
---|
| 335 | delete[] y;
|
---|
| 336 | delete[] dx;
|
---|
| 337 | delete[] nphi;
|
---|
| 338 | //the towers in the forward region
|
---|
[df634866] | 339 | R = tk_length_+contingency_+(contingency_+calo_endcap_thickness_)*calorimeters_.size(); // Z of the muons system = height of the pyramid
|
---|
| 340 | nEtaBins = caloBinning.size();
|
---|
| 341 | // translation to bring the origin of the tower to (0,0,0) (well, not really as the endcap is not yet in place)
|
---|
| 342 | TGeoTranslation* towerdz = new TGeoTranslation(Form("%s_towerdz",name),0.,0.,R/2.-(innerBarrelLength+calo_endcap_coneThickness/2.));
|
---|
[bf6ed57] | 343 | towerdz->RegisterYourself();
|
---|
| 344 | // eta bins... we build one pyramid per eta slice and then translate it nphi times.
|
---|
| 345 | Double_t *r = new Double_t[nEtaBins];
|
---|
| 346 | nphi = new Int_t[nEtaBins];
|
---|
| 347 | etaslice = 0;
|
---|
| 348 | phirotations.clear();
|
---|
[df634866] | 349 | for(set< pair<Double_t, Int_t> >::const_iterator bin=caloBinning.begin(); bin!=caloBinning.end();++bin) {
|
---|
| 350 | if(bin->first<calo_endcap_etamin) continue; // only in the + endcap
|
---|
[bf6ed57] | 351 | r[etaslice] = R*2*exp(-bin->first)/(1-exp(-2*bin->first));
|
---|
| 352 | nphi[etaslice] = bin->second;
|
---|
| 353 | Double_t phiRotationAngle = 360./nphi[etaslice];
|
---|
| 354 | for(int phislice=0;phislice<nphi[etaslice];++phislice) {
|
---|
[df634866] | 355 | phirotations[make_pair(etaslice,phislice)] = new TGeoRotation(Form("%s_forward_phi%d_%d",name,etaslice,phislice),phiRotationAngle*phislice,0.,0.);
|
---|
[bf6ed57] | 356 | phirotations[make_pair(etaslice,phislice)]->RegisterYourself();
|
---|
| 357 | }
|
---|
| 358 | ++etaslice;
|
---|
| 359 | }
|
---|
| 360 | nEtaBins = etaslice;
|
---|
[d2ac201] | 361 | for(int i=0;i<nEtaBins-1;++i) { // loop on the eta slices
|
---|
| 362 | vertices[8] = -r[i+1]*sin(TMath::Pi()/nphi[i]); vertices[9] = r[i+1]*cos(TMath::Pi()/nphi[i]);
|
---|
| 363 | vertices[10] = -r[i]*sin(TMath::Pi()/nphi[i]); vertices[11] = r[i]*cos(TMath::Pi()/nphi[i]);
|
---|
| 364 | vertices[12] = r[i]*sin(TMath::Pi()/nphi[i]); vertices[13] = r[i]*cos(TMath::Pi()/nphi[i]);
|
---|
| 365 | vertices[14] = r[i+1]*sin(TMath::Pi()/nphi[i]); vertices[15] = r[i+1]*cos(TMath::Pi()/nphi[i]);
|
---|
[df634866] | 366 | new TGeoArb8(Form("%sfwdtower%d",name,i),R/2., vertices); // tower in the proper eta slice, at phi=0
|
---|
[bf6ed57] | 367 | // intersection between the tower and the calo_endcap
|
---|
[df634866] | 368 | TGeoCompositeShape *finalfwdtower_cs = new TGeoCompositeShape(Form("%sffwdtower%d_cs",name,i),Form("%sfwdtower%d:%s_towerdz*%s_endcap_cs",name,i,name,name));
|
---|
| 369 | TGeoVolume *finalfwdtower = new TGeoVolume(Form("%sffwdtower%d",name,i),finalfwdtower_cs,calomed_);
|
---|
[bf6ed57] | 370 | finalfwdtower->SetLineColor(kViolet);
|
---|
| 371 | for(int j=0;j<nphi[i];++j) { // loop on the phi slices
|
---|
| 372 | calo_endcap->AddNode(finalfwdtower,j,phirotations[make_pair(i,j)]);
|
---|
| 373 | }
|
---|
| 374 | }
|
---|
| 375 | delete[] r;
|
---|
| 376 | delete[] nphi;
|
---|
| 377 | }
|
---|
| 378 |
|
---|
[a0b6d15] | 379 | void geometry(const char* filename = "delphes_card_CMS.tcl", const char* ParticlePropagator="ParticlePropagator",
|
---|
| 380 | const char* TrackingEfficiency="ChargedHadronTrackingEfficiency",
|
---|
| 381 | const char* MuonEfficiency="MuonEfficiency",
|
---|
| 382 | const char* Calorimeters="Calorimeter")
|
---|
[bf6ed57] | 383 | {
|
---|
| 384 | gSystem->Load("libGeom");
|
---|
| 385 | gSystem->Load("../libDelphes");
|
---|
| 386 | TGeoManager *geom = new TGeoManager("delphes", "Delphes geometry");
|
---|
| 387 |
|
---|
| 388 | // make the top container volume -> designed to contain a "big" detector (ATLAS)
|
---|
| 389 | TGeoVolume *top = geom->MakeBox("TOP", 0, 1500, 1500, 2300);
|
---|
| 390 | geom->SetTopVolume(top);
|
---|
| 391 |
|
---|
| 392 | // build the detector
|
---|
| 393 | Delphes3DGeometry det3D;
|
---|
[a0b6d15] | 394 | det3D.readFile(filename,ParticlePropagator, TrackingEfficiency, MuonEfficiency, Calorimeters);
|
---|
[bf6ed57] | 395 | top->AddNode(det3D.getDetector(true),1);
|
---|
| 396 |
|
---|
| 397 | // draw it
|
---|
| 398 | geom->CloseGeometry();
|
---|
[164f032] | 399 | //top->Draw();
|
---|
| 400 | //TFile* file = new TFile("DelpheGeom.root","RECREATE");
|
---|
| 401 | //top->Write("DelphesGeometry");
|
---|
| 402 | //file->Close();
|
---|
| 403 |
|
---|
| 404 | TEveManager::Create();
|
---|
| 405 |
|
---|
| 406 | //TFile::SetCacheFileDir(".");
|
---|
| 407 | //gGeoManager = gEve->GetGeometry("DelpheGeom.root");
|
---|
| 408 | gGeoManager->DefaultColors();
|
---|
| 409 |
|
---|
| 410 | TGeoVolume* top = gGeoManager->GetTopVolume()->FindNode("Delphes3DGeometry_1")->GetVolume();
|
---|
| 411 |
|
---|
| 412 | TEveGeoTopNode* trk = new TEveGeoTopNode(gGeoManager, top->FindNode("tracker_1"));
|
---|
| 413 | trk->SetVisLevel(6);
|
---|
| 414 | gEve->AddGlobalElement(trk);
|
---|
| 415 |
|
---|
| 416 | TEveGeoTopNode* calo = new TEveGeoTopNode(gGeoManager, top->FindNode("Calorimeter_barrel_1"));
|
---|
| 417 | calo->SetVisLevel(3);
|
---|
| 418 | gEve->AddGlobalElement(calo);
|
---|
| 419 | calo = new TEveGeoTopNode(gGeoManager, top->FindNode("Calorimeter_endcap_1"));
|
---|
| 420 | calo->SetVisLevel(3);
|
---|
| 421 | calo->UseNodeTrans();
|
---|
| 422 | gEve->AddGlobalElement(calo);
|
---|
| 423 | calo = new TEveGeoTopNode(gGeoManager, top->FindNode("Calorimeter_endcap_2"));
|
---|
| 424 | calo->SetVisLevel(3);
|
---|
| 425 | calo->UseNodeTrans();
|
---|
| 426 | gEve->AddGlobalElement(calo);
|
---|
| 427 |
|
---|
| 428 | TEveGeoTopNode* muon = new TEveGeoTopNode(gGeoManager, top->FindNode("muons_barrel_1"));
|
---|
| 429 | muon->SetVisLevel(4);
|
---|
| 430 | gEve->AddGlobalElement(muon);
|
---|
| 431 | muon = new TEveGeoTopNode(gGeoManager, top->FindNode("muons_endcap_1"));
|
---|
| 432 | muon->SetVisLevel(4);
|
---|
| 433 | muon->UseNodeTrans();
|
---|
| 434 | gEve->AddGlobalElement(muon);
|
---|
| 435 | muon = new TEveGeoTopNode(gGeoManager, top->FindNode("muons_endcap_2"));
|
---|
| 436 | muon->SetVisLevel(4);
|
---|
| 437 | muon->UseNodeTrans();
|
---|
| 438 | gEve->AddGlobalElement(muon);
|
---|
| 439 |
|
---|
| 440 | gEve->FullRedraw3D(kTRUE);
|
---|
| 441 |
|
---|
| 442 | // EClipType not exported to CINT (see TGLUtil.h):
|
---|
| 443 | // 0 - no clip, 1 - clip plane, 2 - clip box
|
---|
| 444 | TGLViewer *v = gEve->GetDefaultGLViewer();
|
---|
| 445 | v->GetClipSet()->SetClipType(1);
|
---|
| 446 | v->ColorSet().Background().SetColor(kMagenta+4);
|
---|
| 447 | v->SetGuideState(TGLUtil::kAxesEdge, kTRUE, kFALSE, 0);
|
---|
| 448 | v->RefreshPadEditor(v);
|
---|
| 449 |
|
---|
| 450 | v->CurrentCamera().RotateRad(-1.2, 0.5);
|
---|
| 451 | v->DoDraw();
|
---|
| 452 |
|
---|
[bf6ed57] | 453 | }
|
---|
| 454 |
|
---|