Fork me on GitHub

Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • display/Delphes3DGeometry.cc

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