Fork me on GitHub

Changeset cfc3160 in git for examples/geometry.C


Ignore:
Timestamp:
Oct 15, 2014, 4:01:52 PM (10 years ago)
Author:
Christophe Delaere <christophe.delaere@…>
Branches:
ImprovedOutputFile, Timing, dual_readout, llp, master
Children:
84dd1c8
Parents:
4d999a57
Message:

Migration of the code from script to library

The code still crashes, but debugging will be much easier. Also much
faster!

File:
1 edited

Legend:

Unmodified
Added
Removed

  • examples/geometry.C

    r4d999a57 rcfc3160  
    1515#include "TGeoCone.h"
    1616#include "TGeoArb8.h"
    17 //#include "external/ExRootAnalysis/ExRootConfReader.h"
    18 //#include "external/ExRootAnalysis/ExRootTreeReader.h"
    19 //#include "display/DelphesCaloData.h"
    20 //#include "display/DelphesDisplay.h"
    21 #include "../display/DelphesBranchElement.h"
    22 //#include "classes/DelphesClasses.h"
     17#include "external/ExRootAnalysis/ExRootConfReader.h"
     18#include "external/ExRootAnalysis/ExRootTreeReader.h"
     19#include "display/DelphesCaloData.h"
     20#include "display/DelphesDisplay.h"
     21#include "display/DelphesBranchElement.h"
     22#include "display/Delphes3DGeometry.h"
     23#include "classes/DelphesClasses.h"
    2324#include "TF2.h"
    2425#include "TH1F.h"
    … …  
    8384
    8485DelphesDisplay *gDelphesDisplay = 0;
    85 
    86 /******************************************************************************/
    87 // Construction of the geometry
    88 /******************************************************************************/
    89 
    90 // TODO: asymmetric detector
    91 
    92 class Delphes3DGeometry {
    93    public:
    94      Delphes3DGeometry(TGeoManager *geom = NULL);
    95      ~Delphes3DGeometry() {}
    96 
    97      void readFile(const char* filename, const char* ParticlePropagator="ParticlePropagator",
    98                                          const char* TrackingEfficiency="ChargedHadronTrackingEfficiency",
    99                                          const char* MuonEfficiency="MuonEfficiency",
    100                                          const char* Calorimeters="Calorimeter");
    101 
    102      void loadFromFile(const char* filename, const char* name="DelphesGeometry");
    103      void save(const char* filename, const char* name="DelphesGeometry");
    104 
    105      void setContingency(Double_t contingency) { contingency_ = contingency; }
    106      void setCaloBarrelThickness(Double_t thickness) { calo_barrel_thickness_ = thickness; }
    107      void setCaloEndcapThickness(Double_t thickness) { calo_endcap_thickness_ = thickness; }
    108      void setMuonSystemThickness(Double_t thickness) { muonSystem_thickn_ = thickness; }
    109 
    110      TGeoVolume* getDetector(bool withTowers = true);
    111 
    112      Double_t getTrackerRadius() const { return tk_radius_; }
    113      Double_t getDetectorRadius() const { return muonSystem_radius_; }
    114      Double_t getTrackerHalfLength() const { return tk_length_; }
    115      Double_t getDetectorHalfLength() const { return muonSystem_length_; }
    116      Double_t getBField() const { return tk_Bz_; }
    117      std::pair<TAxis*, TAxis*> getCaloAxes() { return std::make_pair(etaAxis_,phiAxis_); }
    118 
    119    private:
    120      std::pair<Double_t, Double_t> addTracker(TGeoVolume *top);
    121      std::pair<Double_t, Double_t> addCalorimeter(TGeoVolume *top, const char* name, Double_t innerBarrelRadius, Double_t innerBarrelLength, set< pair<Double_t, Int_t> >& caloBinning);
    122      std::pair<Double_t, Double_t> addMuonDets(TGeoVolume *top, const char* name, Double_t innerBarrelRadius, Double_t innerBarrelLength);
    123      void addCaloTowers(TGeoVolume *top, const char* name, Double_t innerBarrelRadius, Double_t innerBarrelLength, set< pair<Double_t, Int_t> >& caloBinning);
    124 
    125    private:
    126 
    127      TGeoManager *geom_;
    128 
    129      TGeoMedium *vacuum_;
    130      TGeoMedium *tkmed_;
    131      TGeoMedium *calomed_;
    132      TGeoMedium *mudetmed_;
    133 
    134      TAxis* etaAxis_;
    135      TAxis* phiAxis_;
    136 
    137      Double_t contingency_;
    138      Double_t calo_barrel_thickness_;
    139      Double_t calo_endcap_thickness_;
    140      Double_t muonSystem_thickn_;
    141      Double_t muonSystem_radius_;
    142      Double_t muonSystem_length_;
    143      Double_t tk_radius_;
    144      Double_t tk_length_;
    145      Double_t tk_etamax_;
    146      Double_t tk_Bz_;
    147 
    148      std::vector<std::string> calorimeters_;
    149      std::vector<std::string> muondets_;
    150 
    151      std::map<std::string, Double_t> muonSystem_etamax_;
    152      std::map<std::string, set< pair<Double_t, Int_t> > > caloBinning_;
    153      
    154 };
    155 
    156 Delphes3DGeometry::Delphes3DGeometry(TGeoManager *geom) {
    157 
    158    //--- the geometry manager
    159    geom_ = geom==NULL? gGeoManager : geom;
    160    //gGeoManager->DefaultColors();
    161 
    162    //--- define some materials
    163    TGeoMaterial *matVacuum = new TGeoMaterial("Vacuum", 0,0,0);
    164    TGeoMaterial *matAl = new TGeoMaterial("Al", 26.98,13,2.7); // placeholder
    165    //TODO: create different materials for different subdetectors???
    166    matVacuum->SetTransparency(85); //TODO: tune
    167    matAl->SetTransparency(85); //TODO: tune
    168 
    169    //--- define some media
    170    TGeoMedium *Vacuum = new TGeoMedium("Vacuum",1, matVacuum);
    171    TGeoMedium *Al = new TGeoMedium("Root Material",2, matAl);
    172    vacuum_ = Vacuum;
    173    tkmed_ = Vacuum; // placeholder
    174    calomed_ = Al;   // placeholder
    175    mudetmed_ = Al;  // placeholder
    176 
    177    // custom parameters
    178    contingency_ = 10.;
    179    calo_barrel_thickness_ = 50.;
    180    calo_endcap_thickness_ = 75.;
    181    muonSystem_thickn_ = 10.;
    182 
    183    // read these parameters from the Delphes Card (with default values)
    184    etaAxis_   = NULL;
    185    phiAxis_   = NULL;
    186    tk_radius_ = 120.;
    187    tk_length_ = 150.;
    188    tk_etamax_ = 3.0;
    189    tk_Bz_     = 1.;
    190    muonSystem_radius_ = 200.;
    191 }
    192 
    193 void Delphes3DGeometry::readFile(const char *configFile,
    194                                  const char* ParticlePropagator, const char* TrackingEfficiency,
    195                                  const char* MuonEfficiency, const char* Calorimeters) {
    196 
    197    ExRootConfReader *confReader = new ExRootConfReader;
    198    confReader->ReadFile(configFile);
    199 
    200    tk_radius_ = confReader->GetDouble(Form("%s::Radius",ParticlePropagator), 1.0)*100.;         // tk_radius
    201    tk_length_ = confReader->GetDouble(Form("%s::HalfLength",ParticlePropagator), 3.0)*100.;     // tk_length
    202    tk_Bz_     = confReader->GetDouble("ParticlePropagator::Bz", 0.0);                           // tk_Bz
    203 
    204    {
    205    TString tkEffFormula = confReader->GetString(Form("%s::EfficiencyFormula",TrackingEfficiency),"abs(eta)<3.0");
    206    tkEffFormula.ReplaceAll("pt","x");
    207    tkEffFormula.ReplaceAll("eta","y");
    208    tkEffFormula.ReplaceAll("phi","0.");
    209    TF2* tkEffFunction = new TF2("tkEff",tkEffFormula,0,1000,-10,10);
    210    TH1F etaHisto("eta","eta",100,5.,-5.);
    211    Double_t pt,eta;
    212    for(int i=0;i<1000;++i) {
    213      tkEffFunction->GetRandom2(pt,eta);
    214      etaHisto.Fill(eta);
    215    }
    216    Int_t bin = -1;
    217    bin = etaHisto.FindFirstBinAbove(0.5);
    218    Double_t etamin = (bin>-1) ? etaHisto.GetBinLowEdge(bin) : -10.;
    219    bin = etaHisto.FindLastBinAbove(0.5);
    220    Double_t etamax = (bin>-1) ? etaHisto.GetBinLowEdge(bin+1) : -10.;
    221    tk_etamax_ = TMath::Max(fabs(etamin),fabs(etamax));                                          // tk_etamax
    222    delete tkEffFunction;
    223    }
    224 
    225    {
    226    muondets_.push_back("muons");
    227    TString muonEffFormula = confReader->GetString(Form("%s::EfficiencyFormula",MuonEfficiency),"abs(eta)<2.0");
    228    muonEffFormula.ReplaceAll("pt","x");
    229    muonEffFormula.ReplaceAll("eta","y");
    230    muonEffFormula.ReplaceAll("phi","0.");
    231    TF2* muEffFunction = new TF2("muEff",muonEffFormula,0,1000,-10,10);
    232    TH1F etaHisto("eta2","eta2",100,5.,-5.);
    233    Double_t pt,eta;
    234    for(int i=0;i<1000;++i) {
    235      muEffFunction->GetRandom2(pt,eta);
    236      etaHisto.Fill(eta);
    237    }
    238    Int_t bin = -1;
    239    bin = etaHisto.FindFirstBinAbove(0.5);
    240    Double_t etamin = (bin>-1) ? etaHisto.GetBinLowEdge(bin) : -10.;
    241    bin = etaHisto.FindLastBinAbove(0.5);
    242    Double_t etamax = (bin>-1) ? etaHisto.GetBinLowEdge(bin+1) : -10.;
    243    muonSystem_etamax_["muons"] = TMath::Max(fabs(etamin),fabs(etamax));                 // muonSystem_etamax
    244    delete muEffFunction;
    245    }
    246 
    247    std::string s(Calorimeters);
    248    std::replace( s.begin(), s.end(), ',', ' ' );
    249    std::istringstream stream( s );
    250    std::string word;
    251    while (stream >> word) calorimeters_.push_back(word);
    252 
    253    caloBinning_.clear();                                                                // calo binning
    254    for(std::vector<std::string>::const_iterator calo=calorimeters_.begin();calo!=calorimeters_.end(); ++calo) {
    255      set< pair<Double_t, Int_t> > caloBinning;
    256      ExRootConfParam paramEtaBins, paramPhiBins;
    257      ExRootConfParam param = confReader->GetParam(Form("%s::EtaPhiBins",calo->c_str()));
    258      Int_t size = param.GetSize();
    259      for(int i = 0; i < size/2; ++i) {
    260        paramEtaBins = param[i*2];
    261        paramPhiBins = param[i*2+1];
    262        assert(paramEtaBins.GetSize()==1);
    263        caloBinning.insert(std::make_pair(paramEtaBins[0].GetDouble(),paramPhiBins.GetSize()-1));
    264      }
    265      caloBinning_[*calo] = caloBinning;
    266    }
    267 
    268    set< pair<Double_t, Int_t> > caloBinning = caloBinning_[*calorimeters_.begin()];
    269    Double_t *etaBins = new Double_t[caloBinning.size()]; // note that this is the eta binning of the first calo
    270    unsigned int ii = 0;
    271    for(set< pair<Double_t, Int_t> >::const_iterator itEtaSet = caloBinning.begin(); itEtaSet != caloBinning.end(); ++itEtaSet) {
    272      etaBins[ii++] = itEtaSet->first;
    273    }
    274    etaAxis_ = new TAxis(caloBinning.size() - 1, etaBins);
    275    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
    276 
    277    muonSystem_radius_ = tk_radius_ + contingency_ + (contingency_+calo_barrel_thickness_)*calorimeters_.size() + muonSystem_thickn_;
    278    muonSystem_length_ = tk_length_ + contingency_ + (contingency_+calo_endcap_thickness_)*calorimeters_.size() + muonSystem_thickn_;
    279 
    280    delete confReader;
    281 
    282 }
    283 
    284 TGeoVolume* Delphes3DGeometry::getDetector(bool withTowers) {
    285    // compute the envelope
    286    Double_t system_radius = tk_radius_+calo_barrel_thickness_+3*contingency_;
    287    Double_t system_length = tk_length_+contingency_+(contingency_+calo_endcap_thickness_)*calorimeters_.size()+contingency_;
    288    // the detector volume
    289    TGeoVolume *top = geom_->MakeBox("Delphes3DGeometry", vacuum_, system_radius, system_radius, system_length);
    290    // build the detector
    291    std::pair<Double_t, Double_t> limits = addTracker(top);
    292    Double_t radius = limits.first;
    293    Double_t length = limits.second;
    294    for(std::vector<std::string>::const_iterator calo = calorimeters_.begin(); calo != calorimeters_.end(); ++calo) {
    295      limits = addCalorimeter(top,calo->c_str(),radius,length,caloBinning_[*calo]);
    296      if (withTowers) {
    297        addCaloTowers(top,calo->c_str(),radius,length,caloBinning_[*calo]);
    298      }
    299      radius = limits.first;
    300      length = limits.second;
    301    }
    302    for(std::vector<std::string>::const_iterator muon = muondets_.begin(); muon != muondets_.end(); ++muon) {
    303      limits = addMuonDets(top, muon->c_str(), radius, length);
    304      radius = limits.first;
    305      length = limits.second;
    306    }
    307    // return the result
    308    return top;
    309 }
    310 
    311 std::pair<Double_t, Double_t> Delphes3DGeometry::addTracker(TGeoVolume *top) {
    312    // tracker: a cylinder with two cones substracted
    313    new TGeoCone("forwardTkAcceptance",(tk_length_/2.+0.05),0.,tk_radius_,(tk_length_)*2.*exp(-tk_etamax_)/(1-exp(-2.*tk_etamax_)),tk_radius_);
    314    TGeoTranslation *tr1  = new TGeoTranslation("tkacc1",0., 0., tk_length_/2.);
    315    tr1->RegisterYourself();
    316    TGeoRotation *negz    = new TGeoRotation("tknegz",0,180,0);
    317    negz->RegisterYourself();
    318    TGeoCombiTrans  *tr2  = new TGeoCombiTrans("tkacc2",0.,0.,-tk_length_/2.,negz);
    319    tr2->RegisterYourself();
    320    TGeoCompositeShape* tracker_cs = new TGeoCompositeShape("tracker_cs","forwardTkAcceptance:tkacc1+forwardTkAcceptance:tkacc2");
    321    TGeoVolume *tracker = new TGeoVolume("tracker",tracker_cs,tkmed_);   
    322    tracker->SetLineColor(kYellow);
    323    top->AddNode(tracker,1);
    324    return std::make_pair(tk_radius_,tk_length_);
    325 }
    326 
    327 std::pair<Double_t, Double_t> Delphes3DGeometry::addCalorimeter(TGeoVolume *top, const char* name,
    328                                                                 Double_t innerBarrelRadius, Double_t innerBarrelLength, set< pair<Double_t, Int_t> >& caloBinning) {
    329    // parameters derived from the inputs
    330    Double_t calo_endcap_etamax        = TMath::Max(fabs(caloBinning.begin()->first),fabs(caloBinning.rbegin()->first));
    331    Double_t calo_barrel_innerRadius   = innerBarrelRadius+contingency_;
    332    Double_t calo_barrel_length        = innerBarrelLength + calo_barrel_thickness_;
    333    Double_t calo_endcap_etamin        = -log(innerBarrelRadius/(2*innerBarrelLength));
    334    Double_t calo_endcap_innerRadius1  = innerBarrelLength*2.*exp(-calo_endcap_etamax)/(1-exp(-2.*calo_endcap_etamax));
    335    Double_t calo_endcap_innerRadius2  = (innerBarrelLength+calo_endcap_thickness_)*2.*exp(-calo_endcap_etamax)/(1-exp(-2.*calo_endcap_etamax));
    336    Double_t calo_endcap_outerRadius1  = innerBarrelRadius;
    337    Double_t calo_endcap_outerRadius2  = innerBarrelRadius+calo_barrel_thickness_;
    338    Double_t calo_endcap_coneThickness = TMath::Min(calo_barrel_thickness_ * (1-exp(-2.*calo_endcap_etamin)) / (2.*exp(-calo_endcap_etamin)), calo_endcap_thickness_);
    339    Double_t calo_endcap_diskThickness = TMath::Max(0.,calo_endcap_thickness_-calo_endcap_coneThickness);
    340 
    341    // calorimeters: tube truncated in eta + cones
    342    new TGeoTube(Form("%s_barrel_cylinder",name),calo_barrel_innerRadius,calo_barrel_innerRadius+calo_barrel_thickness_,calo_barrel_length);
    343    new TGeoCone(Form("%s_endcap_cone",name),calo_endcap_coneThickness/2.,calo_endcap_innerRadius1,calo_endcap_outerRadius1,calo_endcap_innerRadius2,calo_endcap_outerRadius2);
    344    new TGeoTube(Form("%s_endcap_disk",name),calo_endcap_innerRadius2,tk_radius_+calo_barrel_thickness_,calo_endcap_diskThickness/2.);
    345    TGeoTranslation *tr1 = new TGeoTranslation(Form("%s_tr1",name),0., 0., (calo_endcap_coneThickness+calo_endcap_diskThickness)/2.);
    346    tr1->RegisterYourself();
    347    TGeoCompositeShape *calo_endcap_cs = new TGeoCompositeShape(Form("%s_endcap_cs",name),Form("%s_endcap_cone+%s_endcap_disk:%s_tr1",name,name,name));
    348    TGeoTranslation *trc1 = new TGeoTranslation(Form("%s_endcap1_position",name),0.,0., innerBarrelLength+calo_endcap_coneThickness/2.);
    349    trc1->RegisterYourself();
    350    TGeoRotation *negz = new TGeoRotation(Form("%s_negz",name),0,180,0);
    351    TGeoCombiTrans  *trc2 = new TGeoCombiTrans(Form("%s_endcap2_position",name),0.,0.,-(innerBarrelLength+calo_endcap_coneThickness/2.),negz);
    352    trc2->RegisterYourself();
    353    TGeoTranslation *trc1c = new TGeoTranslation(Form("%s_endcap1_position_cont",name),0.,0., innerBarrelLength+calo_endcap_coneThickness/2.+contingency_);
    354    trc1c->RegisterYourself();
    355    TGeoCombiTrans  *trc2c = new TGeoCombiTrans(Form("%s_endcap2_position_cont",name),0.,0.,-(innerBarrelLength+calo_endcap_coneThickness/2.)-contingency_,negz);
    356    trc2c->RegisterYourself();
    357    TGeoVolume *calo_endcap = new TGeoVolume(Form("%s_endcap",name),calo_endcap_cs,calomed_);
    358    TGeoCompositeShape *calo_barrel_cs = new TGeoCompositeShape(Form("%s_barrel_cs",name),
    359                                                                Form("%s_barrel_cylinder-%s_endcap_cs:%s_endcap1_position-%s_endcap_cs:%s_endcap2_position",name,name,name,name,name));
    360    TGeoVolume *calo_barrel = new TGeoVolume(Form("%s_barrel",name),calo_barrel_cs,calomed_);
    361    calo_endcap->SetLineColor(kViolet);
    362    calo_endcap->SetFillColor(kViolet);
    363    calo_barrel->SetLineColor(kRed);
    364    top->AddNode(calo_endcap,1,trc1c);
    365    top->AddNode(calo_endcap,2,trc2c);
    366    top->AddNode(calo_barrel,1);
    367    return std::make_pair(calo_barrel_innerRadius+calo_barrel_thickness_,innerBarrelLength+calo_endcap_thickness_+contingency_);
    368 }
    369 
    370 std::pair<Double_t, Double_t> Delphes3DGeometry::addMuonDets(TGeoVolume *top, const char* name, Double_t innerBarrelRadius, Double_t innerBarrelLength) {
    371    // muon system: tube + disks
    372    Double_t muonSystem_radius = innerBarrelRadius + contingency_;
    373    Double_t muonSystem_length = innerBarrelLength + contingency_;
    374    Double_t muonSystem_rmin   = muonSystem_length*2.*exp(-muonSystem_etamax_[name])/(1-exp(-2.*muonSystem_etamax_[name]));
    375    TGeoVolume *muon_barrel = geom_->MakeTube(Form("%s_barrel",name),mudetmed_,muonSystem_radius,muonSystem_radius+muonSystem_thickn_,muonSystem_length);
    376    muon_barrel->SetLineColor(kBlue);
    377    top->AddNode(muon_barrel,1);
    378    TGeoVolume *muon_endcap = geom_->MakeTube(Form("%s_endcap",name),mudetmed_,muonSystem_rmin,muonSystem_radius+muonSystem_thickn_,muonSystem_thickn_/2.);
    379    muon_endcap->SetLineColor(kBlue);
    380    TGeoTranslation *trm1 = new TGeoTranslation(Form("%sEndcap1_position",name),0.,0.,muonSystem_length);
    381    trm1->RegisterYourself();
    382    TGeoTranslation *trm2 = new TGeoTranslation(Form("%sEndcap2_position",name),0.,0.,-muonSystem_length);
    383    trm1->RegisterYourself();
    384    top->AddNode(muon_endcap,1,trm1);
    385    top->AddNode(muon_endcap,2,trm2);
    386    return std::make_pair(muonSystem_radius,muonSystem_length);
    387 }
    388 
    389 void Delphes3DGeometry::addCaloTowers(TGeoVolume *top, const char* name,
    390                                                        Double_t innerBarrelRadius, Double_t innerBarrelLength, set< pair<Double_t, Int_t> >& caloBinning) {
    391 
    392    TGeoVolume* calo_endcap = top->GetNode(Form("%s_endcap_1",name))->GetVolume();
    393    TGeoVolume* calo_barrel = top->GetNode(Form("%s_barrel_1",name))->GetVolume();
    394    Double_t calo_endcap_etamin = -log(innerBarrelRadius/(2*innerBarrelLength));
    395    Double_t calo_endcap_coneThickness = TMath::Min(calo_barrel_thickness_ * (1-exp(-2.*calo_endcap_etamin)) / (2.*exp(-calo_endcap_etamin)), calo_endcap_thickness_);
    396 
    397    // calo towers in the barrel
    398    Double_t vertices[16] = {0.,0.,0.,0.,0.,0.,0.,0.}; // summit of the pyramid
    399    Double_t R  = tk_radius_ + contingency_+(contingency_+calo_barrel_thickness_)*calorimeters_.size(); // radius of the muons system = height of the pyramid
    400    Int_t nEtaBins = caloBinning.size();
    401    // this rotation is to make the tower point "up"
    402    TGeoRotation* initTowerRot = new TGeoRotation(Form("%s_initTowerRot",name),0.,90.,0.);
    403    TGeoCombiTrans* initTower  = new TGeoCombiTrans(Form("%s_initTower",name),0.,-R/2.,0.,initTowerRot);
    404    initTower->RegisterYourself();
    405    // eta bins... we build one pyramid per eta slice and then translate it nphi times.
    406    // phi bins represented by rotations around z
    407    Double_t *y = new Double_t[nEtaBins];
    408    Double_t *dx = new Double_t[nEtaBins];
    409    Int_t *nphi = new Int_t[nEtaBins];
    410    Int_t etaslice = 0;
    411    std::map<std::pair<int,int>, TGeoRotation*> phirotations;
    412    for(set< pair<Double_t, Int_t> >::const_iterator bin=caloBinning.begin(); bin!=caloBinning.end();++bin) {
    413      if(abs(bin->first)>calo_endcap_etamin) continue; // only in the barrel
    414      nphi[etaslice] = bin->second;
    415      y[etaslice] = 0.5*R*(1-exp(-2*bin->first))/exp(-bin->first);
    416      Double_t phiRotationAngle = 360./nphi[etaslice];
    417      dx[etaslice] = R*tan(TMath::Pi()*phiRotationAngle/360.);
    418      for(int phislice=0;phislice<nphi[etaslice];++phislice) {
    419        phirotations[make_pair(etaslice,phislice)] = new TGeoRotation(Form("%s_phi%d_%d",name,etaslice,phislice),phiRotationAngle*phislice,0.,0.);
    420        phirotations[make_pair(etaslice,phislice)]->RegisterYourself();
    421      }
    422      ++etaslice;
    423    }
    424    nEtaBins = etaslice;
    425    for(int i=0;i<nEtaBins-1;++i) { // loop on the eta slices
    426      vertices[8]  = -dx[i]; vertices[9]  = y[i];
    427      vertices[10] = -dx[i]; vertices[11] = y[i+1];
    428      vertices[12] =  dx[i]; vertices[13] = y[i+1];
    429      vertices[14] =  dx[i]; vertices[15] = y[i];
    430      new TGeoArb8(Form("%s_tower%d",name,i),R/2., vertices); // tower in the proper eta slice, at phi=0
    431      // intersection between the tower and the calo_barrel
    432      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));
    433      TGeoVolume *finaltower = new TGeoVolume(Form("%s_ftower%d",name,i),finaltower_cs,calomed_);
    434      finaltower->SetLineColor(kRed);
    435      for(int j=0;j<nphi[i];++j) { // loop on the phi slices
    436        calo_barrel->AddNode(finaltower,j,phirotations[make_pair(i,j)]);
    437      }
    438    }
    439    delete[] y;
    440    delete[] dx;
    441    delete[] nphi;
    442    //the towers in the forward region
    443    R  = tk_length_+contingency_+(contingency_+calo_endcap_thickness_)*calorimeters_.size(); // Z of the muons system = height of the pyramid
    444    nEtaBins = caloBinning.size();
    445    // translation to bring the origin of the tower to (0,0,0) (well, not really as the endcap is not yet in place)
    446    TGeoTranslation* towerdz = new TGeoTranslation(Form("%s_towerdz",name),0.,0.,R/2.-(innerBarrelLength+calo_endcap_coneThickness/2.));
    447    towerdz->RegisterYourself();
    448    // eta bins... we build one pyramid per eta slice and then translate it nphi times.
    449    Double_t *r = new Double_t[nEtaBins];
    450    nphi = new Int_t[nEtaBins];
    451    etaslice = 0;
    452    phirotations.clear();
    453    for(set< pair<Double_t, Int_t> >::const_iterator bin=caloBinning.begin(); bin!=caloBinning.end();++bin) {
    454      if(bin->first<calo_endcap_etamin) continue; // only in the + endcap
    455      r[etaslice] = R*2*exp(-bin->first)/(1-exp(-2*bin->first));
    456      nphi[etaslice] = bin->second;
    457      Double_t phiRotationAngle = 360./nphi[etaslice];
    458      for(int phislice=0;phislice<nphi[etaslice];++phislice) {
    459        phirotations[make_pair(etaslice,phislice)] = new TGeoRotation(Form("%s_forward_phi%d_%d",name,etaslice,phislice),phiRotationAngle*phislice,0.,0.);
    460        phirotations[make_pair(etaslice,phislice)]->RegisterYourself();
    461      }
    462      ++etaslice;
    463    }
    464    nEtaBins = etaslice;
    465    for(int i=0;i<nEtaBins-1;++i) { // loop on the eta slices
    466      vertices[8]  = -r[i+1]*sin(TMath::Pi()/nphi[i]); vertices[9]  = r[i+1]*cos(TMath::Pi()/nphi[i]);
    467      vertices[10] = -r[i]*sin(TMath::Pi()/nphi[i]);   vertices[11] = r[i]*cos(TMath::Pi()/nphi[i]);
    468      vertices[12] =  r[i]*sin(TMath::Pi()/nphi[i]);   vertices[13] = r[i]*cos(TMath::Pi()/nphi[i]);
    469      vertices[14] =  r[i+1]*sin(TMath::Pi()/nphi[i]); vertices[15] = r[i+1]*cos(TMath::Pi()/nphi[i]);
    470      new TGeoArb8(Form("%sfwdtower%d",name,i),R/2., vertices); // tower in the proper eta slice, at phi=0
    471      // intersection between the tower and the calo_endcap
    472      TGeoCompositeShape *finalfwdtower_cs = new TGeoCompositeShape(Form("%sffwdtower%d_cs",name,i),Form("%sfwdtower%d:%s_towerdz*%s_endcap_cs",name,i,name,name));
    473      TGeoVolume *finalfwdtower = new TGeoVolume(Form("%sffwdtower%d",name,i),finalfwdtower_cs,calomed_);
    474      finalfwdtower->SetLineColor(kViolet);
    475      for(int j=0;j<nphi[i];++j) { // loop on the phi slices
    476        calo_endcap->AddNode(finalfwdtower,j,phirotations[make_pair(i,j)]);
    477      }
    478    }
    479    delete[] r;
    480    delete[] nphi;
    481 }
    48286
    48387/******************************************************************************/
    … …  
    940544   // load the libraries
    941545   gSystem->Load("libGeom");
    942    gSystem->Load("../libDelphes");
     546   //gSystem->Load("../libDelphes");
    943547   gSystem->Load("../libDelphesDisplay");
    944548
Note: See TracChangeset for help on using the changeset viewer.