| 1 | /*
|
|---|
| 2 | * Delphes: a framework for fast simulation of a generic collider experiment
|
|---|
| 3 | * Copyright (C) 2012-2014 Universite catholique de Louvain (UCL), Belgium
|
|---|
| 4 | *
|
|---|
| 5 | * This program is free software: you can redistribute it and/or modify
|
|---|
| 6 | * it under the terms of the GNU General Public License as published by
|
|---|
| 7 | * the Free Software Foundation, either version 3 of the License, or
|
|---|
| 8 | * (at your option) any later version.
|
|---|
| 9 | *
|
|---|
| 10 | * This program is distributed in the hope that it will be useful,
|
|---|
| 11 | * but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|---|
| 12 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|---|
| 13 | * GNU General Public License for more details.
|
|---|
| 14 | *
|
|---|
| 15 | * You should have received a copy of the GNU General Public License
|
|---|
| 16 | * along with this program. If not, see <http://www.gnu.org/licenses/>.
|
|---|
| 17 | */
|
|---|
| 18 |
|
|---|
| 19 | #include <set>
|
|---|
| 20 | #include <map>
|
|---|
| 21 | #include <utility>
|
|---|
| 22 | #include <vector>
|
|---|
| 23 | #include <algorithm>
|
|---|
| 24 | #include <sstream>
|
|---|
| 25 | #include <cassert>
|
|---|
| 26 |
|
|---|
| 27 | #include "TAxis.h"
|
|---|
| 28 | #include "TGeoManager.h"
|
|---|
| 29 | #include "TGeoVolume.h"
|
|---|
| 30 | #include "TGeoMedium.h"
|
|---|
| 31 | #include "TGeoNode.h"
|
|---|
| 32 | #include "TGeoCompositeShape.h"
|
|---|
| 33 | #include "TGeoMatrix.h"
|
|---|
| 34 | #include "TGeoTube.h"
|
|---|
| 35 | #include "TGeoCone.h"
|
|---|
| 36 | #include "TGeoArb8.h"
|
|---|
| 37 | #include "TF2.h"
|
|---|
| 38 | #include "TFormula.h"
|
|---|
| 39 | #include "TH1F.h"
|
|---|
| 40 | #include "TMath.h"
|
|---|
| 41 | #include "TString.h"
|
|---|
| 42 |
|
|---|
| 43 | #include "display/Delphes3DGeometry.h"
|
|---|
| 44 |
|
|---|
| 45 | #include "classes/DelphesClasses.h"
|
|---|
| 46 | #include "external/ExRootAnalysis/ExRootConfReader.h"
|
|---|
| 47 |
|
|---|
| 48 | using namespace std;
|
|---|
| 49 |
|
|---|
| 50 | Delphes3DGeometry::Delphes3DGeometry(TGeoManager *geom, bool transp) {
|
|---|
| 51 |
|
|---|
| 52 | //--- the geometry manager
|
|---|
| 53 | geom_ = geom==NULL? gGeoManager : geom;
|
|---|
| 54 | //gGeoManager->DefaultColors();
|
|---|
| 55 |
|
|---|
| 56 | //--- define some materials
|
|---|
| 57 | TGeoMaterial *matVacuum = new TGeoMaterial("Vacuum", 0,0,0);
|
|---|
| 58 | TGeoMaterial *matAl = new TGeoMaterial("Al", 26.98,13,2.7); // placeholder
|
|---|
| 59 | if(transp) {
|
|---|
| 60 | matVacuum->SetTransparency(85);
|
|---|
| 61 | matAl->SetTransparency(85);
|
|---|
| 62 | }
|
|---|
| 63 |
|
|---|
| 64 | //--- define some media
|
|---|
| 65 | TGeoMedium *Vacuum = new TGeoMedium("Vacuum",1, matVacuum);
|
|---|
| 66 | TGeoMedium *Al = new TGeoMedium("Root Material",2, matAl);
|
|---|
| 67 | vacuum_ = Vacuum;
|
|---|
| 68 | tkmed_ = Vacuum; // placeholder
|
|---|
| 69 | calomed_ = Al; // placeholder
|
|---|
| 70 | mudetmed_ = Al; // placeholder
|
|---|
| 71 |
|
|---|
| 72 | // custom parameters
|
|---|
| 73 | contingency_ = 10.;
|
|---|
| 74 | calo_barrel_thickness_ = 50.;
|
|---|
| 75 | calo_endcap_thickness_ = 75.;
|
|---|
| 76 | muonSystem_thickn_ = 10.;
|
|---|
| 77 |
|
|---|
| 78 | // read these parameters from the Delphes Card (with default values)
|
|---|
| 79 | etaAxis_ = NULL;
|
|---|
| 80 | phiAxis_ = NULL;
|
|---|
| 81 | tk_radius_ = 120.;
|
|---|
| 82 | tk_length_ = 150.;
|
|---|
| 83 | tk_etamax_ = 3.0;
|
|---|
| 84 | tk_Bz_ = 1.;
|
|---|
| 85 | muonSystem_radius_ = 200.;
|
|---|
| 86 | }
|
|---|
| 87 |
|
|---|
| 88 | void Delphes3DGeometry::readFile(const char *configFile,
|
|---|
| 89 | const char* ParticlePropagator, const char* TrackingEfficiency,
|
|---|
| 90 | const char* MuonEfficiency, const char* Calorimeters) {
|
|---|
| 91 |
|
|---|
| 92 | ExRootConfReader *confReader = new ExRootConfReader;
|
|---|
| 93 | confReader->ReadFile(configFile);
|
|---|
| 94 |
|
|---|
| 95 | tk_radius_ = confReader->GetDouble(Form("%s::Radius",ParticlePropagator), 1.0)*100.; // tk_radius
|
|---|
| 96 | tk_length_ = confReader->GetDouble(Form("%s::HalfLength",ParticlePropagator), 3.0)*100.; // tk_length
|
|---|
| 97 | tk_Bz_ = confReader->GetDouble("ParticlePropagator::Bz", 0.0); // tk_Bz
|
|---|
| 98 |
|
|---|
| 99 | TString buffer;
|
|---|
| 100 | const char *it;
|
|---|
| 101 |
|
|---|
| 102 |
|
|---|
| 103 | {
|
|---|
| 104 | TString tkEffFormula = confReader->GetString(Form("%s::EfficiencyFormula",TrackingEfficiency),"abs(eta)<3.0");
|
|---|
| 105 | tkEffFormula.ReplaceAll("pt","x");
|
|---|
| 106 | tkEffFormula.ReplaceAll("eta","y");
|
|---|
| 107 | tkEffFormula.ReplaceAll("phi","0.");
|
|---|
| 108 |
|
|---|
| 109 | buffer.Clear();
|
|---|
| 110 | for(it = tkEffFormula.Data(); *it; ++it)
|
|---|
| 111 | {
|
|---|
| 112 | if(*it == ' ' || *it == '\t' || *it == '\r' || *it == '\n' || *it == '\\' ) continue;
|
|---|
| 113 | buffer.Append(*it);
|
|---|
| 114 | }
|
|---|
| 115 |
|
|---|
| 116 | TF2* tkEffFunction = new TF2("tkEff",buffer,0,1000,-10,10);
|
|---|
| 117 | TH1F etaHisto("eta","eta",100,5.,-5.);
|
|---|
| 118 | Double_t pt,eta;
|
|---|
| 119 | for(int i=0;i<1000;++i) {
|
|---|
| 120 | tkEffFunction->GetRandom2(pt,eta);
|
|---|
| 121 | etaHisto.Fill(eta);
|
|---|
| 122 | }
|
|---|
| 123 | Int_t bin = -1;
|
|---|
| 124 | bin = etaHisto.FindFirstBinAbove(0.5);
|
|---|
| 125 | Double_t etamin = (bin>-1) ? etaHisto.GetBinLowEdge(bin) : -10.;
|
|---|
| 126 | bin = etaHisto.FindLastBinAbove(0.5);
|
|---|
| 127 | Double_t etamax = (bin>-1) ? etaHisto.GetBinLowEdge(bin+1) : -10.;
|
|---|
| 128 | tk_etamax_ = TMath::Max(fabs(etamin),fabs(etamax)); // tk_etamax
|
|---|
| 129 | delete tkEffFunction;
|
|---|
| 130 | }
|
|---|
| 131 |
|
|---|
| 132 | {
|
|---|
| 133 | muondets_.push_back("muons");
|
|---|
| 134 | TString muonEffFormula = confReader->GetString(Form("%s::EfficiencyFormula",MuonEfficiency),"abs(eta)<2.0");
|
|---|
| 135 | muonEffFormula.ReplaceAll("pt","x");
|
|---|
| 136 | muonEffFormula.ReplaceAll("eta","y");
|
|---|
| 137 | muonEffFormula.ReplaceAll("phi","0.");
|
|---|
| 138 |
|
|---|
| 139 | buffer.Clear();
|
|---|
| 140 | for(it = muonEffFormula.Data(); *it; ++it)
|
|---|
| 141 | {
|
|---|
| 142 | if(*it == ' ' || *it == '\t' || *it == '\r' || *it == '\n' || *it == '\\' ) continue;
|
|---|
| 143 | buffer.Append(*it);
|
|---|
| 144 | }
|
|---|
| 145 |
|
|---|
| 146 | TF2* muEffFunction = new TF2("muEff",buffer,0,1000,-10,10);
|
|---|
| 147 | TH1F etaHisto("eta2","eta2",100,5.,-5.);
|
|---|
| 148 | Double_t pt,eta;
|
|---|
| 149 | for(int i=0;i<1000;++i) {
|
|---|
| 150 | muEffFunction->GetRandom2(pt,eta);
|
|---|
| 151 | etaHisto.Fill(eta);
|
|---|
| 152 | }
|
|---|
| 153 | Int_t bin = -1;
|
|---|
| 154 | bin = etaHisto.FindFirstBinAbove(0.5);
|
|---|
| 155 | Double_t etamin = (bin>-1) ? etaHisto.GetBinLowEdge(bin) : -10.;
|
|---|
| 156 | bin = etaHisto.FindLastBinAbove(0.5);
|
|---|
| 157 | Double_t etamax = (bin>-1) ? etaHisto.GetBinLowEdge(bin+1) : -10.;
|
|---|
| 158 | muonSystem_etamax_["muons"] = TMath::Max(fabs(etamin),fabs(etamax)); // muonSystem_etamax
|
|---|
| 159 | delete muEffFunction;
|
|---|
| 160 | }
|
|---|
| 161 |
|
|---|
| 162 | std::string s(Calorimeters);
|
|---|
| 163 | std::replace( s.begin(), s.end(), ',', ' ' );
|
|---|
| 164 | std::istringstream stream( s );
|
|---|
| 165 | std::string word;
|
|---|
| 166 | while (stream >> word) calorimeters_.push_back(word);
|
|---|
| 167 |
|
|---|
| 168 | caloBinning_.clear(); // calo binning
|
|---|
| 169 | for(std::vector<std::string>::const_iterator calo=calorimeters_.begin();calo!=calorimeters_.end(); ++calo) {
|
|---|
| 170 | set< pair<Double_t, Int_t> > caloBinning;
|
|---|
| 171 | ExRootConfParam paramEtaBins, paramPhiBins;
|
|---|
| 172 | ExRootConfParam param = confReader->GetParam(Form("%s::EtaPhiBins",calo->c_str()));
|
|---|
| 173 | Int_t size = param.GetSize();
|
|---|
| 174 | for(int i = 0; i < size/2; ++i) {
|
|---|
| 175 | paramEtaBins = param[i*2];
|
|---|
| 176 | paramPhiBins = param[i*2+1];
|
|---|
| 177 | assert(paramEtaBins.GetSize()==1);
|
|---|
| 178 | caloBinning.insert(std::make_pair(paramEtaBins[0].GetDouble(),paramPhiBins.GetSize()-1));
|
|---|
| 179 | }
|
|---|
| 180 | caloBinning_[*calo] = caloBinning;
|
|---|
| 181 | }
|
|---|
| 182 |
|
|---|
| 183 | set< pair<Double_t, Int_t> > caloBinning = caloBinning_[*calorimeters_.begin()];
|
|---|
| 184 | Double_t *etaBins = new Double_t[caloBinning.size()]; // note that this is the eta binning of the first calo
|
|---|
| 185 | unsigned int ii = 0;
|
|---|
| 186 | for(set< pair<Double_t, Int_t> >::const_iterator itEtaSet = caloBinning.begin(); itEtaSet != caloBinning.end(); ++itEtaSet) {
|
|---|
| 187 | etaBins[ii++] = itEtaSet->first;
|
|---|
| 188 | }
|
|---|
| 189 | etaAxis_ = new TAxis(caloBinning.size() - 1, etaBins);
|
|---|
| 190 | 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
|
|---|
| 191 |
|
|---|
| 192 | muonSystem_radius_ = tk_radius_ + contingency_ + (contingency_+calo_barrel_thickness_)*calorimeters_.size() + muonSystem_thickn_;
|
|---|
| 193 | muonSystem_length_ = tk_length_ + contingency_ + (contingency_+calo_endcap_thickness_)*calorimeters_.size() + muonSystem_thickn_;
|
|---|
| 194 |
|
|---|
| 195 | delete confReader;
|
|---|
| 196 |
|
|---|
| 197 | }
|
|---|
| 198 |
|
|---|
| 199 | TGeoVolume* Delphes3DGeometry::getDetector(bool withTowers) {
|
|---|
| 200 | // compute the envelope
|
|---|
| 201 | Double_t system_radius = tk_radius_+calo_barrel_thickness_+3*contingency_;
|
|---|
| 202 | Double_t system_length = tk_length_+contingency_+(contingency_+calo_endcap_thickness_)*calorimeters_.size()+contingency_;
|
|---|
| 203 | // the detector volume
|
|---|
| 204 | TGeoVolume *top = geom_->MakeBox("Delphes3DGeometry", vacuum_, system_radius, system_radius, system_length);
|
|---|
| 205 | // build the detector
|
|---|
| 206 | std::pair<Double_t, Double_t> limits = addTracker(top);
|
|---|
| 207 | Double_t radius = limits.first;
|
|---|
| 208 | Double_t length = limits.second;
|
|---|
| 209 | for(std::vector<std::string>::const_iterator calo = calorimeters_.begin(); calo != calorimeters_.end(); ++calo) {
|
|---|
| 210 | limits = addCalorimeter(top,calo->c_str(),radius,length,caloBinning_[*calo]);
|
|---|
| 211 | if (withTowers) {
|
|---|
| 212 | addCaloTowers(top,calo->c_str(),radius,length,caloBinning_[*calo]);
|
|---|
| 213 | }
|
|---|
| 214 | radius = limits.first;
|
|---|
| 215 | length = limits.second;
|
|---|
| 216 | }
|
|---|
| 217 | for(std::vector<std::string>::const_iterator muon = muondets_.begin(); muon != muondets_.end(); ++muon) {
|
|---|
| 218 | limits = addMuonDets(top, muon->c_str(), radius, length);
|
|---|
| 219 | radius = limits.first;
|
|---|
| 220 | length = limits.second;
|
|---|
| 221 | }
|
|---|
| 222 | // return the result
|
|---|
| 223 | return top;
|
|---|
| 224 | }
|
|---|
| 225 |
|
|---|
| 226 | std::pair<Double_t, Double_t> Delphes3DGeometry::addTracker(TGeoVolume *top) {
|
|---|
| 227 | // tracker: a cylinder with two cones substracted
|
|---|
| 228 | new TGeoCone("forwardTkAcceptance",(tk_length_/2.+0.05),0.,tk_radius_,(tk_length_)*2.*exp(-tk_etamax_)/(1-exp(-2.*tk_etamax_)),tk_radius_);
|
|---|
| 229 | TGeoTranslation *tr1 = new TGeoTranslation("tkacc1",0., 0., tk_length_/2.);
|
|---|
| 230 | tr1->RegisterYourself();
|
|---|
| 231 | TGeoRotation *negz = new TGeoRotation("tknegz",0,180,0);
|
|---|
| 232 | negz->RegisterYourself();
|
|---|
| 233 | TGeoCombiTrans *tr2 = new TGeoCombiTrans("tkacc2",0.,0.,-tk_length_/2.,negz);
|
|---|
| 234 | tr2->RegisterYourself();
|
|---|
| 235 | TGeoCompositeShape* tracker_cs = new TGeoCompositeShape("tracker_cs","forwardTkAcceptance:tkacc1+forwardTkAcceptance:tkacc2");
|
|---|
| 236 | TGeoVolume *tracker = new TGeoVolume("tracker",tracker_cs,tkmed_);
|
|---|
| 237 | tracker->SetLineColor(kYellow);
|
|---|
| 238 | top->AddNode(tracker,1);
|
|---|
| 239 | return std::make_pair(tk_radius_,tk_length_);
|
|---|
| 240 | }
|
|---|
| 241 |
|
|---|
| 242 | std::pair<Double_t, Double_t> Delphes3DGeometry::addCalorimeter(TGeoVolume *top, const char* name,
|
|---|
| 243 | Double_t innerBarrelRadius, Double_t innerBarrelLength, set< pair<Double_t, Int_t> >& caloBinning) {
|
|---|
| 244 | // parameters derived from the inputs
|
|---|
| 245 | Double_t calo_endcap_etamax = TMath::Max(fabs(caloBinning.begin()->first),fabs(caloBinning.rbegin()->first));
|
|---|
| 246 | Double_t calo_barrel_innerRadius = innerBarrelRadius+contingency_;
|
|---|
| 247 | Double_t calo_barrel_length = innerBarrelLength + calo_barrel_thickness_;
|
|---|
| 248 | Double_t calo_endcap_etamin = -log(innerBarrelRadius/(2*innerBarrelLength));
|
|---|
| 249 | Double_t calo_endcap_innerRadius1 = innerBarrelLength*2.*exp(-calo_endcap_etamax)/(1-exp(-2.*calo_endcap_etamax));
|
|---|
| 250 | Double_t calo_endcap_innerRadius2 = (innerBarrelLength+calo_endcap_thickness_)*2.*exp(-calo_endcap_etamax)/(1-exp(-2.*calo_endcap_etamax));
|
|---|
| 251 | Double_t calo_endcap_outerRadius1 = innerBarrelRadius;
|
|---|
| 252 | Double_t calo_endcap_outerRadius2 = innerBarrelRadius+calo_barrel_thickness_;
|
|---|
| 253 | Double_t calo_endcap_coneThickness = TMath::Min(calo_barrel_thickness_ * (1-exp(-2.*calo_endcap_etamin)) / (2.*exp(-calo_endcap_etamin)), calo_endcap_thickness_);
|
|---|
| 254 | Double_t calo_endcap_diskThickness = TMath::Max(0.,calo_endcap_thickness_-calo_endcap_coneThickness);
|
|---|
| 255 |
|
|---|
| 256 | // calorimeters: tube truncated in eta + cones
|
|---|
| 257 | new TGeoTube(Form("%s_barrel_cylinder",name),calo_barrel_innerRadius,calo_barrel_innerRadius+calo_barrel_thickness_,calo_barrel_length);
|
|---|
| 258 | new TGeoCone(Form("%s_endcap_cone",name),calo_endcap_coneThickness/2.,calo_endcap_innerRadius1,calo_endcap_outerRadius1,calo_endcap_innerRadius2,calo_endcap_outerRadius2);
|
|---|
| 259 | new TGeoTube(Form("%s_endcap_disk",name),calo_endcap_innerRadius2,tk_radius_+calo_barrel_thickness_,calo_endcap_diskThickness/2.);
|
|---|
| 260 | TGeoTranslation *tr1 = new TGeoTranslation(Form("%s_tr1",name),0., 0., (calo_endcap_coneThickness+calo_endcap_diskThickness)/2.);
|
|---|
| 261 | tr1->RegisterYourself();
|
|---|
| 262 | TGeoCompositeShape *calo_endcap_cs = new TGeoCompositeShape(Form("%s_endcap_cs",name),Form("%s_endcap_cone+%s_endcap_disk:%s_tr1",name,name,name));
|
|---|
| 263 | TGeoTranslation *trc1 = new TGeoTranslation(Form("%s_endcap1_position",name),0.,0., innerBarrelLength+calo_endcap_coneThickness/2.);
|
|---|
| 264 | trc1->RegisterYourself();
|
|---|
| 265 | TGeoRotation *negz = new TGeoRotation(Form("%s_negz",name),0,180,0);
|
|---|
| 266 | TGeoCombiTrans *trc2 = new TGeoCombiTrans(Form("%s_endcap2_position",name),0.,0.,-(innerBarrelLength+calo_endcap_coneThickness/2.),negz);
|
|---|
| 267 | trc2->RegisterYourself();
|
|---|
| 268 | TGeoTranslation *trc1c = new TGeoTranslation(Form("%s_endcap1_position_cont",name),0.,0., innerBarrelLength+calo_endcap_coneThickness/2.+contingency_);
|
|---|
| 269 | trc1c->RegisterYourself();
|
|---|
| 270 | TGeoCombiTrans *trc2c = new TGeoCombiTrans(Form("%s_endcap2_position_cont",name),0.,0.,-(innerBarrelLength+calo_endcap_coneThickness/2.)-contingency_,negz);
|
|---|
| 271 | trc2c->RegisterYourself();
|
|---|
| 272 | TGeoVolume *calo_endcap = new TGeoVolume(Form("%s_endcap",name),calo_endcap_cs,calomed_);
|
|---|
| 273 | TGeoCompositeShape *calo_barrel_cs = new TGeoCompositeShape(Form("%s_barrel_cs",name),
|
|---|
| 274 | Form("%s_barrel_cylinder-%s_endcap_cs:%s_endcap1_position-%s_endcap_cs:%s_endcap2_position",name,name,name,name,name));
|
|---|
| 275 | TGeoVolume *calo_barrel = new TGeoVolume(Form("%s_barrel",name),calo_barrel_cs,calomed_);
|
|---|
| 276 | calo_endcap->SetLineColor(kViolet);
|
|---|
| 277 | calo_endcap->SetFillColor(kViolet);
|
|---|
| 278 | calo_barrel->SetLineColor(kRed);
|
|---|
| 279 | top->AddNode(calo_endcap,1,trc1c);
|
|---|
| 280 | top->AddNode(calo_endcap,2,trc2c);
|
|---|
| 281 | top->AddNode(calo_barrel,1);
|
|---|
| 282 | return std::make_pair(calo_barrel_innerRadius+calo_barrel_thickness_,innerBarrelLength+calo_endcap_thickness_+contingency_);
|
|---|
| 283 | }
|
|---|
| 284 |
|
|---|
| 285 | std::pair<Double_t, Double_t> Delphes3DGeometry::addMuonDets(TGeoVolume *top, const char* name, Double_t innerBarrelRadius, Double_t innerBarrelLength) {
|
|---|
| 286 | // muon system: tube + disks
|
|---|
| 287 | Double_t muonSystem_radius = innerBarrelRadius + contingency_;
|
|---|
| 288 | Double_t muonSystem_length = innerBarrelLength + contingency_;
|
|---|
| 289 | Double_t muonSystem_rmin = muonSystem_length*2.*exp(-muonSystem_etamax_[name])/(1-exp(-2.*muonSystem_etamax_[name]));
|
|---|
| 290 | TGeoVolume *muon_barrel = geom_->MakeTube(Form("%s_barrel",name),mudetmed_,muonSystem_radius,muonSystem_radius+muonSystem_thickn_,muonSystem_length);
|
|---|
| 291 | muon_barrel->SetLineColor(kBlue);
|
|---|
| 292 | top->AddNode(muon_barrel,1);
|
|---|
| 293 | TGeoVolume *muon_endcap = geom_->MakeTube(Form("%s_endcap",name),mudetmed_,muonSystem_rmin,muonSystem_radius+muonSystem_thickn_,muonSystem_thickn_/2.);
|
|---|
| 294 | muon_endcap->SetLineColor(kBlue);
|
|---|
| 295 | TGeoTranslation *trm1 = new TGeoTranslation(Form("%sEndcap1_position",name),0.,0.,muonSystem_length);
|
|---|
| 296 | trm1->RegisterYourself();
|
|---|
| 297 | TGeoTranslation *trm2 = new TGeoTranslation(Form("%sEndcap2_position",name),0.,0.,-muonSystem_length);
|
|---|
| 298 | trm1->RegisterYourself();
|
|---|
| 299 | top->AddNode(muon_endcap,1,trm1);
|
|---|
| 300 | top->AddNode(muon_endcap,2,trm2);
|
|---|
| 301 | return std::make_pair(muonSystem_radius,muonSystem_length);
|
|---|
| 302 | }
|
|---|
| 303 |
|
|---|
| 304 | void Delphes3DGeometry::addCaloTowers(TGeoVolume *top, const char* name,
|
|---|
| 305 | Double_t innerBarrelRadius, Double_t innerBarrelLength, set< pair<Double_t, Int_t> >& caloBinning) {
|
|---|
| 306 |
|
|---|
| 307 | TGeoVolume* calo_endcap = top->GetNode(Form("%s_endcap_1",name))->GetVolume();
|
|---|
| 308 | TGeoVolume* calo_barrel = top->GetNode(Form("%s_barrel_1",name))->GetVolume();
|
|---|
| 309 | Double_t calo_endcap_etamin = -log(innerBarrelRadius/(2*innerBarrelLength));
|
|---|
| 310 | Double_t calo_endcap_coneThickness = TMath::Min(calo_barrel_thickness_ * (1-exp(-2.*calo_endcap_etamin)) / (2.*exp(-calo_endcap_etamin)), calo_endcap_thickness_);
|
|---|
| 311 |
|
|---|
| 312 | // calo towers in the barrel
|
|---|
| 313 | Double_t vertices[16] = {0.,0.,0.,0.,0.,0.,0.,0.}; // summit of the pyramid
|
|---|
| 314 | Double_t R = tk_radius_ + contingency_+(contingency_+calo_barrel_thickness_)*calorimeters_.size(); // radius of the muons system = height of the pyramid
|
|---|
| 315 | Int_t nEtaBins = caloBinning.size();
|
|---|
| 316 | // this rotation is to make the tower point "up"
|
|---|
| 317 | TGeoRotation* initTowerRot = new TGeoRotation(Form("%s_initTowerRot",name),0.,90.,0.);
|
|---|
| 318 | TGeoCombiTrans* initTower = new TGeoCombiTrans(Form("%s_initTower",name),0.,-R/2.,0.,initTowerRot);
|
|---|
| 319 | initTower->RegisterYourself();
|
|---|
| 320 | // eta bins... we build one pyramid per eta slice and then translate it nphi times.
|
|---|
| 321 | // phi bins represented by rotations around z
|
|---|
| 322 | Double_t *y = new Double_t[nEtaBins];
|
|---|
| 323 | Double_t *dx = new Double_t[nEtaBins];
|
|---|
| 324 | Int_t *nphi = new Int_t[nEtaBins];
|
|---|
| 325 | Int_t etaslice = 0;
|
|---|
| 326 | std::map<std::pair<int,int>, TGeoRotation*> phirotations;
|
|---|
| 327 | for(set< pair<Double_t, Int_t> >::const_iterator bin=caloBinning.begin(); bin!=caloBinning.end();++bin) {
|
|---|
| 328 | if(abs(bin->first)>calo_endcap_etamin) continue; // only in the barrel
|
|---|
| 329 | nphi[etaslice] = bin->second;
|
|---|
| 330 | y[etaslice] = 0.5*R*(1-exp(-2*bin->first))/exp(-bin->first);
|
|---|
| 331 | Double_t phiRotationAngle = 360./nphi[etaslice];
|
|---|
| 332 | dx[etaslice] = R*tan(TMath::Pi()*phiRotationAngle/360.);
|
|---|
| 333 | for(int phislice=0;phislice<nphi[etaslice];++phislice) {
|
|---|
| 334 | phirotations[make_pair(etaslice,phislice)] = new TGeoRotation(Form("%s_phi%d_%d",name,etaslice,phislice),phiRotationAngle*phislice,0.,0.);
|
|---|
| 335 | phirotations[make_pair(etaslice,phislice)]->RegisterYourself();
|
|---|
| 336 | }
|
|---|
| 337 | ++etaslice;
|
|---|
| 338 | }
|
|---|
| 339 | nEtaBins = etaslice;
|
|---|
| 340 | for(int i=0;i<nEtaBins-1;++i) { // loop on the eta slices
|
|---|
| 341 | vertices[8] = -dx[i]; vertices[9] = y[i];
|
|---|
| 342 | vertices[10] = -dx[i]; vertices[11] = y[i+1];
|
|---|
| 343 | vertices[12] = dx[i]; vertices[13] = y[i+1];
|
|---|
| 344 | vertices[14] = dx[i]; vertices[15] = y[i];
|
|---|
| 345 | new TGeoArb8(Form("%s_tower%d",name,i),R/2., vertices); // tower in the proper eta slice, at phi=0
|
|---|
| 346 | // intersection between the tower and the calo_barrel
|
|---|
| 347 | 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));
|
|---|
| 348 | TGeoVolume *finaltower = new TGeoVolume(Form("%s_ftower%d",name,i),finaltower_cs,calomed_);
|
|---|
| 349 | finaltower->SetLineColor(kRed);
|
|---|
| 350 | for(int j=0;j<nphi[i];++j) { // loop on the phi slices
|
|---|
| 351 | calo_barrel->AddNode(finaltower,j,phirotations[make_pair(i,j)]);
|
|---|
| 352 | }
|
|---|
| 353 | }
|
|---|
| 354 | delete[] y;
|
|---|
| 355 | delete[] dx;
|
|---|
| 356 | delete[] nphi;
|
|---|
| 357 | //the towers in the forward region
|
|---|
| 358 | R = tk_length_+contingency_+(contingency_+calo_endcap_thickness_)*calorimeters_.size(); // Z of the muons system = height of the pyramid
|
|---|
| 359 | nEtaBins = caloBinning.size();
|
|---|
| 360 | // translation to bring the origin of the tower to (0,0,0) (well, not really as the endcap is not yet in place)
|
|---|
| 361 | TGeoTranslation* towerdz = new TGeoTranslation(Form("%s_towerdz",name),0.,0.,R/2.-(innerBarrelLength+calo_endcap_coneThickness/2.));
|
|---|
| 362 | towerdz->RegisterYourself();
|
|---|
| 363 | // eta bins... we build one pyramid per eta slice and then translate it nphi times.
|
|---|
| 364 | Double_t *r = new Double_t[nEtaBins];
|
|---|
| 365 | nphi = new Int_t[nEtaBins];
|
|---|
| 366 | etaslice = 0;
|
|---|
| 367 | phirotations.clear();
|
|---|
| 368 | for(set< pair<Double_t, Int_t> >::const_iterator bin=caloBinning.begin(); bin!=caloBinning.end();++bin) {
|
|---|
| 369 | if(bin->first<calo_endcap_etamin) continue; // only in the + endcap
|
|---|
| 370 | r[etaslice] = R*2*exp(-bin->first)/(1-exp(-2*bin->first));
|
|---|
| 371 | nphi[etaslice] = bin->second;
|
|---|
| 372 | Double_t phiRotationAngle = 360./nphi[etaslice];
|
|---|
| 373 | for(int phislice=0;phislice<nphi[etaslice];++phislice) {
|
|---|
| 374 | phirotations[make_pair(etaslice,phislice)] = new TGeoRotation(Form("%s_forward_phi%d_%d",name,etaslice,phislice),phiRotationAngle*phislice,0.,0.);
|
|---|
| 375 | phirotations[make_pair(etaslice,phislice)]->RegisterYourself();
|
|---|
| 376 | }
|
|---|
| 377 | ++etaslice;
|
|---|
| 378 | }
|
|---|
| 379 | nEtaBins = etaslice;
|
|---|
| 380 | for(int i=0;i<nEtaBins-1;++i) { // loop on the eta slices
|
|---|
| 381 | vertices[8] = -r[i+1]*sin(TMath::Pi()/nphi[i]); vertices[9] = r[i+1]*cos(TMath::Pi()/nphi[i]);
|
|---|
| 382 | vertices[10] = -r[i]*sin(TMath::Pi()/nphi[i]); vertices[11] = r[i]*cos(TMath::Pi()/nphi[i]);
|
|---|
| 383 | vertices[12] = r[i]*sin(TMath::Pi()/nphi[i]); vertices[13] = r[i]*cos(TMath::Pi()/nphi[i]);
|
|---|
| 384 | vertices[14] = r[i+1]*sin(TMath::Pi()/nphi[i]); vertices[15] = r[i+1]*cos(TMath::Pi()/nphi[i]);
|
|---|
| 385 | new TGeoArb8(Form("%sfwdtower%d",name,i),R/2., vertices); // tower in the proper eta slice, at phi=0
|
|---|
| 386 | // intersection between the tower and the calo_endcap
|
|---|
| 387 | TGeoCompositeShape *finalfwdtower_cs = new TGeoCompositeShape(Form("%sffwdtower%d_cs",name,i),Form("%sfwdtower%d:%s_towerdz*%s_endcap_cs",name,i,name,name));
|
|---|
| 388 | TGeoVolume *finalfwdtower = new TGeoVolume(Form("%sffwdtower%d",name,i),finalfwdtower_cs,calomed_);
|
|---|
| 389 | finalfwdtower->SetLineColor(kViolet);
|
|---|
| 390 | for(int j=0;j<nphi[i];++j) { // loop on the phi slices
|
|---|
| 391 | calo_endcap->AddNode(finalfwdtower,j,phirotations[make_pair(i,j)]);
|
|---|
| 392 | }
|
|---|
| 393 | }
|
|---|
| 394 | delete[] r;
|
|---|
| 395 | delete[] nphi;
|
|---|
| 396 | }
|
|---|
| 397 |
|
|---|
