g4simulation/g4mvtx/PHG4EICMvtxDetector.cc

0001 #include "PHG4EICMvtxDetector.h"
0002
0003 #include "PHG4MvtxDefs.h"
0004 #include "PHG4MvtxDisplayAction.h"
0005
0006 #include <mvtx/CylinderGeom_Mvtx.h>
0007
0008 #include <g4detectors/PHG4CylinderGeomContainer.h>
0009
0010 #include <phparameter/PHParameters.h>
0011 #include <phparameter/PHParametersContainer.h>
0012
0013 #include <g4main/PHG4Detector.h>       // for PHG4Detector
0014 #include <g4main/PHG4DisplayAction.h>  // for PHG4DisplayAction
0015 #include <g4main/PHG4Subsystem.h>      // for PHG4Subsystem
0016
0017 #include <phool/PHCompositeNode.h>
0018 #include <phool/PHIODataNode.h>
0019 #include <phool/PHNode.h>          // for PHNode
0020 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0021 #include <phool/PHObject.h>        // for PHObject
0022 #include <phool/getClass.h>
0023
0024 #include <Geant4/G4AssemblyVolume.hh>
0025 #include <Geant4/G4LogicalVolume.hh>
0026 #include <Geant4/G4Material.hh>
0027 #include <Geant4/G4RotationMatrix.hh>  // for G4RotationMatrix
0028 #include <Geant4/G4String.hh>          // for G4String
0029 #include <Geant4/G4SystemOfUnits.hh>
0030 #include <Geant4/G4ThreeVector.hh>      // for G4ThreeVector
0031 #include <Geant4/G4Transform3D.hh>      // for G4Transform3D
0032 #include <Geant4/G4Types.hh>            // for G4double
0033 #include <Geant4/G4VPhysicalVolume.hh>  // for G4VPhysicalVolume
0034
0035 // xerces has some shadowed variables
0036 #pragma GCC diagnostic push
0037 #pragma GCC diagnostic ignored "-Wshadow"
0038 #include <Geant4/G4GDMLParser.hh>
0039 #include <Geant4/G4GDMLReadStructure.hh>  // for G4GDMLReadStructure
0040 #pragma GCC diagnostic pop
0041
0042 #include <cmath>
0043 #include <cstdio>    // for sprintf
0044 #include <iostream>  // for operator<<, basic...
0045 #include <memory>
0046 #include <sstream>
0047 #include <utility>  // for pair, make_pair
0048 #include <vector>   // for vector, vector<>:...
0049
0050 using namespace std;
0051
0052 PHG4EICMvtxDetector::PHG4EICMvtxDetector(PHG4Subsystem* subsys, PHCompositeNode* Node, const PHParametersContainer* _paramsContainer, const std::string& dnam)
0053   : PHG4Detector(subsys, Node, dnam)
0054   , m_DisplayAction(dynamic_cast<PHG4MvtxDisplayAction*>(subsys->GetDisplayAction()))
0055   , m_ParamsContainer(_paramsContainer)
0056   , m_StaveGeometryFile(_paramsContainer->GetParameters(PHG4MvtxDefs::GLOBAL)->get_string_param("stave_geometry_file"))
0057 {
0058   //  Verbosity(2);
0059
0060   if (Verbosity() > 0)
0061   {
0062     cout << "PHG4EICMvtxDetector constructor called" << endl;
0063   }
0064
0065   if (Verbosity() > 10)
0066   {
0067     cout << " cm " << cm << " mm " << mm << endl;
0068   }
0069   for (int ilayer = 0; ilayer < n_Layers; ++ilayer)
0070   {
0071     const PHParameters* params = m_ParamsContainer->GetParameters(ilayer);
0072     m_IsLayerActive[ilayer] = params->get_int_param("active");
0073     m_IsLayerAbsorberActive[ilayer] = params->get_int_param("absorberactive");
0074     m_IsBlackHole[ilayer] = params->get_int_param("blackhole");
0075     m_N_staves[ilayer] = params->get_int_param("N_staves");
0076     m_nominal_radius[ilayer] = params->get_double_param("layer_nominal_radius");
0077     m_nominal_phitilt[ilayer] = params->get_double_param("phitilt");
0078     m_nominal_phi0[ilayer] = params->get_double_param("phi0");
0079   }
0080   /*
0081   const PHParameters* alpide_params = m_ParamsContainer->GetParameters(PHG4MvtxDefs::ALPIDE_SEGMENTATION);
0082   m_PixelX = alpide_params->get_double_param("pixel_x");
0083   m_PixelZ = alpide_params->get_double_param("pixel_z");
0084   m_PixelThickness = alpide_params->get_double_param("pixel_thickness");
0085   */
0086   if (Verbosity() > 0)
0087   {
0088     cout << "PHG4EICMvtxDetector constructor: making Mvtx detector. " << endl;
0089   }
0090 }
0091
0092 //_______________________________________________________________
0093 //_______________________________________________________________
0094 int PHG4EICMvtxDetector::IsSensor(G4VPhysicalVolume* volume) const
0095 {
0096   // Is this volume one of the sensors?
0097   // Checks if pointer matches one of our stored sensors for this layer
0098   if (m_SensorPV.find(volume) != m_SensorPV.end())
0099   {
0100     if (Verbosity() > 0)
0101     {
0102       cout << " -- PHG4MvtxTDetector::IsSensor --" << endl;
0103       cout << " volume Name : " << volume->GetName() << endl;
0104       cout << " -----------------------------------------" << endl;
0105     }
0106     return 1;
0107   }
0108
0109   return 0;
0110 }
0111
0112 int PHG4EICMvtxDetector::IsInMvtx(G4VPhysicalVolume* volume, int& layer, int& stave) const
0113 {
0114   // Does this stave belong to this layer?
0115   // Since the Assembly volume read from GDML does not give unique pointers
0116   // to sensors, we need to check the stave, which is unique
0117   auto iter = m_StavePV.find(volume);
0118   if (iter != m_StavePV.end())
0119   {
0120     tie(layer, stave) = iter->second;
0121     if (Verbosity() > 0)
0122     {
0123       cout << " -- PHG4EICMvtxDetector::IsInMvtx --" << endl;
0124       cout << " layer: " << layer << endl;
0125       cout << " stave: " << stave << endl;
0126       cout << " volume Name : " << volume->GetName() << endl;
0127       cout << " stave Name  : " << iter->first->GetName() << endl;
0128       cout << " -----------------------------------------" << endl;
0129     }
0130     return 1;
0131   }
0132
0133   return 0;
0134 }
0135
0136 int PHG4EICMvtxDetector::get_layer(int index) const
0137 {
0138   // Get Mvtx layer from stave index in the Mvtx
0139   // Mvtx stave index start from 0, i.e index = 0 for stave 0 in layer 0
0140   int lay = 0;
0141   while (!(index < m_N_staves[lay]))
0142   {
0143     index -= m_N_staves[lay];
0144     lay++;
0145   }
0146   return lay;
0147 }
0148
0149 int PHG4EICMvtxDetector::get_stave(int index) const
0150 {
0151   // Get stave index in the layer from stave index in the Mvtx
0152   int lay = 0;
0153   while (!(index < m_N_staves[lay]))
0154   {
0155     index -= m_N_staves[lay];
0156     lay++;
0157   }
0158   return index;
0159 }
0160
0161 void PHG4EICMvtxDetector::ConstructMe(G4LogicalVolume* logicWorld)
0162 {
0163   // This is called from PHG4PhenixDetector::Construct()
0164
0165   if (Verbosity() > 0)
0166   {
0167     cout << endl
0168          << "PHG4EICMvtxDetector::Construct called for Mvtx " << endl;
0169   }
0170   // the tracking layers are placed directly in the world volume, since some layers are (touching) double layers
0171   // this reads in the ITS stave geometry from a file and constructs the layer from it
0172   ConstructMvtx(logicWorld);
0173
0174   // This object provides the strip center locations when given the ladder segment and strip internal cordinates in the sensor
0175   AddGeometryNode();
0176   return;
0177 }
0178
0179 int PHG4EICMvtxDetector::ConstructMvtx(G4LogicalVolume* trackerenvelope)
0180 {
0181   if (Verbosity() > 0)
0182   {
0183     cout << " PHG4EICMvtxDetector::ConstructMvtx:" << endl;
0184     cout << endl;
0185   }
0186   //===================================
0187   // Import the stave physical volume here
0188   //===================================
0189
0190   // import the staves from the gemetry file
0191   std::unique_ptr<G4GDMLReadStructure> reader(new G4GDMLReadStructure());
0192   G4GDMLParser gdmlParser(reader.get());
0193   gdmlParser.Read(m_StaveGeometryFile, false);
0194
0195   // figure out which assembly we want
0196   char assemblyname[500];
0197   sprintf(assemblyname, "MVTXStave");
0198
0199   if (Verbosity() > 0)
0200   {
0201     cout << "Geting the stave assembly named " << assemblyname << endl;
0202   }
0203   G4AssemblyVolume* av_ITSUStave = reader->GetAssembly(assemblyname);
0204
0205   for (unsigned short ilayer = 0; ilayer < n_Layers; ++ilayer)
0206   {
0207     if (m_IsLayerActive[ilayer])
0208     {
0209       if (Verbosity() > 0)
0210       {
0211         cout << endl;
0212         cout << " Constructing Layer " << ilayer << endl;
0213       }
0214       ConstructMvtx_Layer(ilayer, av_ITSUStave, trackerenvelope);
0215     }
0216   }
0217   FillPVArray(av_ITSUStave);
0218   SetDisplayProperty(av_ITSUStave);
0219
0220   return 0;
0221 }
0222
0223 int PHG4EICMvtxDetector::ConstructMvtx_Layer(int layer, G4AssemblyVolume* av_ITSUStave, G4LogicalVolume*& trackerenvelope)
0224 {
0225   //=========================================
0226   // Now we populate the whole layer with the staves
0227   //==========================================
0228
0229   int N_staves = m_N_staves[layer];
0230   G4double layer_nominal_radius = m_nominal_radius[layer];
0231   G4double phitilt = m_nominal_phitilt[layer];
0232   G4double phi0 = m_nominal_phi0[layer];  // YCM: azimuthal offset for the first stave
0233
0234   if (N_staves < 0)
0235   {
0236     // The user did not specify how many staves to use for this layer, so we calculate the best value
0237     // We choose a phistep that fits N_staves into this radius, but with an arclength separation of AT LEAST arcstep
0238     // ideally, the radius would be chosen so that numstaves = N_staves exactly, to get the closest spacing of staves possible without overlaps
0239     double arcstep = 12.25;
0240     double numstaves = 2.0 * M_PI * layer_nominal_radius / arcstep;  // this is just to print out
0241     N_staves = int(2.0 * M_PI * layer_nominal_radius / arcstep);     // this is the number of staves used
0242
0243     // Also suggest the ideal radius for this layer
0244     if (Verbosity() > 0)
0245     {
0246       cout << " Calculated N_staves for layer " /*<< layer*/
0247            << " layer_nominal_radius " << layer_nominal_radius
0248            << " ITS arcstep " << arcstep
0249            << " circumference divided by arcstep  " << numstaves
0250            << " N_staves " << N_staves
0251            << endl;
0252       cout << "A radius for this layer of " << (double) N_staves * arcstep / (2.0 * M_PI) + 0.01 << " or "
0253            << (double) (N_staves + 1) * arcstep / (2.0 * M_PI) + 0.01 << " would produce  perfect stave spacing" << endl;
0254     }
0255   }
0256
0257   G4double phistep = get_phistep(layer);  // this produces even stave spacing
0258   double z_location = 0.0;
0259
0260   if (Verbosity() > 0)
0261   {
0262     cout << " layer " /*<< layer*/
0263          << " layer_nominal_radius " << layer_nominal_radius
0264          << " N_staves " << N_staves
0265          << " phistep " << phistep
0266          << " phitilt " << phitilt
0267          << " phi0    " << phi0
0268          << endl;
0269   }
0270
0271   // The stave starts out at (0,0,0) and oriented so that the sensors face upward in y
0272   // So we need to rotate the sensor 90 degrees before placing it using phi_offset
0273   // note that the gdml file uses a negative phi_offset - different coord system, apparently - the following works
0274   double phi_offset = M_PI / 2.0;
0275
0276   for (int iphi = 0; iphi < N_staves; iphi++)
0277   {
0278     // Place the ladder segment envelopes at the correct z and phi
0279     // This is the azimuthal angle at which we place the stave
0280     // phi0 is the azimuthal offset for the first stave
0281     G4double phi_rotation = phi0 + (double) iphi * phistep;
0282
0283     G4RotationMatrix Ra;
0284     G4ThreeVector Ta;
0285
0286     if (Verbosity() > 0)
0287     {
0288       cout << "phi_offset = " << phi_offset << " iphi " << iphi << " phi_rotation = " << phi_rotation << " phitilt " << phitilt << endl;
0289     }
0290     // It  is first rotated in phi by the azimuthal angle phi_rotation, plus the 90 degrees needed to point the face of the sensor  at the origin,  plus the tilt (if a tilt is appropriate)
0291
0292     // note - if this is layer 0-2, phitilt is the additional tilt for clearance. Otherwise it is zero
0293     Ra.rotateZ(phi_rotation + phi_offset + phitilt);
0294     // Then translated as follows
0295
0296     Ta.setX(layer_nominal_radius * cos(phi_rotation));
0297     Ta.setY(layer_nominal_radius * sin(phi_rotation));
0298     Ta.setZ(z_location);
0299
0300     if (Verbosity() > 0)
0301     {
0302       cout << " iphi " << iphi << " phi_rotation " << phi_rotation
0303            << " x " << layer_nominal_radius * cos(phi_rotation)
0304            << " y " << layer_nominal_radius * sin(phi_rotation)
0305            << " z " << z_location
0306            << endl;
0307     }
0308     G4Transform3D Tr(Ra, Ta);
0309
0310     av_ITSUStave->MakeImprint(trackerenvelope, Tr, 0, OverlapCheck());
0311   }
0312
0313   if (Verbosity() > 0)
0314   {
0315     cout << "This layer has a total of " << N_staves << " staves" << endl;
0316   }
0317   return 0;
0318 }
0319
0320 void PHG4EICMvtxDetector::SetDisplayProperty(G4AssemblyVolume* av)
0321 {
0322   //  cout <<"SetDisplayProperty - G4AssemblyVolume w/ TotalImprintedVolumes "<<av->TotalImprintedVolumes()
0323   //   <<"/"<<av->GetImprintsCount()<<endl;
0324
0325   std::vector<G4VPhysicalVolume*>::iterator it = av->GetVolumesIterator();
0326
0327   int nDaughters = av->TotalImprintedVolumes();
0328   for (int i = 0; i < nDaughters; ++i, ++it)
0329   {
0330     //  cout <<"SetDisplayProperty - AV["<<i<<"] = "<<(*it)->GetName()<<endl;
0331     G4VPhysicalVolume* pv = (*it);
0332
0333     G4LogicalVolume* worldLogical = pv->GetLogicalVolume();
0334     SetDisplayProperty(worldLogical);
0335   }
0336 }
0337
0338 void PHG4EICMvtxDetector::SetDisplayProperty(G4LogicalVolume* lv)
0339 {
0340   string material_name(lv->GetMaterial()->GetName());
0341
0342   if (Verbosity() >= 50)
0343   {
0344     cout << "SetDisplayProperty - LV " << lv->GetName() << " built with "
0345          << material_name << endl;
0346   }
0347   vector<string> matname = {"SI", "KAPTON", "ALUMINUM", "Carbon", "M60J3K", "WATER"};
0348   bool found = false;
0349   for (const string& nam : matname)
0350   {
0351     if (material_name.find(nam) != std::string::npos)
0352     {
0353       m_DisplayAction->AddVolume(lv, nam);
0354       if (Verbosity() >= 50)
0355       {
0356         cout << "SetDisplayProperty - LV " << lv->GetName() << " display with " << nam << endl;
0357       }
0358       found = true;
0359       break;
0360     }
0361   }
0362   if (!found)
0363   {
0364     m_DisplayAction->AddVolume(lv, "ANYTHING_ELSE");
0365   }
0366   int nDaughters = lv->GetNoDaughters();
0367   for (int i = 0; i < nDaughters; ++i)
0368   {
0369     G4VPhysicalVolume* pv = lv->GetDaughter(i);
0370
0371     // cout <<"SetDisplayProperty - PV["<<i<<"] = "<<pv->GetName()<<endl;
0372
0373     G4LogicalVolume* worldLogical = pv->GetLogicalVolume();
0374     SetDisplayProperty(worldLogical);
0375   }
0376 }
0377
0378 void PHG4EICMvtxDetector::AddGeometryNode()
0379 {
0380   int active = 0;
0381   for (auto& isAct : m_IsLayerActive)
0382   {
0383     active |= isAct;
0384   }
0385   if (active)  // At least one layer is active
0386   {
0387     ostringstream geonode;
0388     geonode << "CYLINDERGEOM_" << ((m_SuperDetector != "NONE") ? m_SuperDetector : m_Detector);
0389     PHG4CylinderGeomContainer* geo = findNode::getClass<PHG4CylinderGeomContainer>(topNode(), geonode.str().c_str());
0390     if (!geo)
0391     {
0392       geo = new PHG4CylinderGeomContainer();
0393       PHNodeIterator iter(topNode());
0394       PHCompositeNode* runNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "RUN"));
0395       PHIODataNode<PHObject>* newNode = new PHIODataNode<PHObject>(geo, geonode.str().c_str(), "PHObject");
0396       runNode->addNode(newNode);
0397     }
0398     // here in the detector class we have internal units(mm), convert to cm
0399     // before putting into the geom object
0400     for (unsigned short ilayer = 0; ilayer < n_Layers; ++ilayer)
0401     {
0402       CylinderGeom_Mvtx* mygeom = new CylinderGeom_Mvtx(ilayer,
0403                                                         m_N_staves[ilayer],
0404                                                         m_nominal_radius[ilayer] / cm,
0405                                                         get_phistep(ilayer) / rad,
0406                                                         m_nominal_phitilt[ilayer] / rad,
0407                                                         m_nominal_phi0[ilayer] / rad);
0408       geo->AddLayerGeom(ilayer, mygeom);
0409     }  // loop per layers
0410     if (Verbosity())
0411     {
0412       geo->identify();
0413     }
0414   }  // is active
0415 }  // AddGeometryNode
0416
0417 void PHG4EICMvtxDetector::FillPVArray(G4AssemblyVolume* av)
0418 {
0419   if (Verbosity() > 0)
0420   {
0421     cout << "-- FillPVArray --" << endl;
0422   }
0423   std::vector<G4VPhysicalVolume*>::iterator it = av->GetVolumesIterator();
0424
0425   int nDaughters = av->TotalImprintedVolumes();
0426   for (int i = 0; i < nDaughters; ++i, ++it)
0427   {
0428     G4VPhysicalVolume* pv = (*it);
0429
0430     G4LogicalVolume* worldLogical = pv->GetLogicalVolume();
0431     // we only care about the staves, which contain the sensors, not the structures
0432     if (pv->GetName().find("MVTXHalfStave_pv") != string::npos)
0433     {
0434       int layer = get_layer(m_StavePV.size());
0435       int stave = get_stave(m_StavePV.size());
0436
0437       m_StavePV.insert(make_pair(pv, make_tuple(layer, stave)));
0438
0439       if (Verbosity() > 0)
0440       {
0441         cout << "Mvtx layer id " << layer << endl;
0442         cout << "Stave in layer id " << stave << endl;
0443         cout << "Mvtx stave count " << m_StavePV.size() << endl;
0444         cout << "FillPVArray - AV[" << i << "] = " << (*it)->GetName() << endl;
0445         cout << "              LV[" << i << "] = " << worldLogical->GetName() << endl;
0446       }
0447
0448       FindSensor(worldLogical);
0449     }
0450     else  // in case of stave structure
0451     {
0452       if (Verbosity() > 0)
0453       {
0454         cout << "FillPVArray - AV[" << i << "] = " << (*it)->GetName() << endl;
0455         cout << "              LV[" << i << "] = " << worldLogical->GetName() << endl;
0456       }
0457     }
0458   }
0459 }
0460
0461 void PHG4EICMvtxDetector::FindSensor(G4LogicalVolume* lv)
0462 {
0463   int nDaughters = lv->GetNoDaughters();
0464   for (int i = 0; i < nDaughters; ++i)
0465   {
0466     G4VPhysicalVolume* pv = lv->GetDaughter(i);
0467     if (Verbosity() > 0)
0468     {
0469       cout << "                 PV[" << i << "]: " << pv->GetName() << endl;
0470     }
0471     if (pv->GetName().find("MVTXSensor_") != string::npos)
0472     {
0473       m_SensorPV.insert(pv);
0474
0475       if (Verbosity() > 0)
0476       {
0477         cout << "                      Adding Sensor Vol <" << pv->GetName() << " (" << m_SensorPV.size() << ")>" << endl;
0478       }
0479     }
0480
0481     G4LogicalVolume* worldLogical = pv->GetLogicalVolume();
0482
0483     if (Verbosity() > 0)
0484     {
0485       cout << "                 LV[" << i << "]: " << worldLogical->GetName() << endl;
0486     }
0487     FindSensor(worldLogical);
0488   }
0489 }