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