Fork me on GitHub

Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • display/Delphes3DGeometry.cc

    rf53a4d2 r77e9ae1  
    1717 */
    1818
     19#include <algorithm>
     20#include <cassert>
     21#include <map>
    1922#include <set>
    20 #include <map>
     23#include <sstream>
    2124#include <utility>
    2225#include <vector>
    23 #include <algorithm>
    24 #include <sstream>
    25 #include <cassert>
    2626
    2727#include "TAxis.h"
     28#include "TF2.h"
     29#include "TFormula.h"
     30#include "TGeoArb8.h"
     31#include "TGeoCompositeShape.h"
     32#include "TGeoCone.h"
    2833#include "TGeoManager.h"
    29 #include "TGeoVolume.h"
     34#include "TGeoMatrix.h"
    3035#include "TGeoMedium.h"
    3136#include "TGeoNode.h"
    32 #include "TGeoCompositeShape.h"
    33 #include "TGeoMatrix.h"
    3437#include "TGeoTube.h"
    35 #include "TGeoCone.h"
    36 #include "TGeoArb8.h"
    37 #include "TF2.h"
    38 #include "TFormula.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    //--- 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.;
     50Delphes3DGeometry::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.;
    8688}
    8789
    8890void 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 
     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
     205TGeoVolume *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
     236std::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
     253std::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
     297std::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
     317void 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}
Note: See TracChangeset for help on using the changeset viewer.