//***** Generated by Geant4 Geometry Editor at Thu Apr 04 19:40:34 GMT+03:00 //2002 ***** //------HeaderFile- #include "KM3Detector.hh" //#include "G4UnitsTable.hh" //#include "G4PhysicalConstants.hh" //#include "G4SystemOfUnits.hh" //#include "aanet.hh" #include "G4VUserDetectorConstruction.hh" #include "KM3SD.hh" #include "KM3StackingAction.hh" #include "globals.hh" #include "G4Material.hh" #include "G4MaterialTable.hh" #include "G4Element.hh" #include "G4ElementTable.hh" #include "G4Box.hh" #include "G4Sphere.hh" #include "G4Tubs.hh" #include "G4LogicalVolume.hh" #include "G4ThreeVector.hh" #include "G4PVPlacement.hh" #include "G4PVReplica.hh" #include "G4SDManager.hh" #include "G4VisAttributes.hh" #include "G4Colour.hh" #include "G4UserLimits.hh" #include "G4RegionStore.hh" // newgeant #include "G4Processor/GDMLProcessor.h" #include "G4GDMLParser.hh" //newgeant #include "G4LogicalVolumeStore.hh" #include "G4PhysicalVolumeStore.hh" #include "G4GeometryManager.hh" static const string versn = GIT_VERSION; static const string dater = BUILD_DATE; KM3Detector::KM3Detector() { allCathods = new KM3Cathods(); allStoreys = new std::vector; allOMs = new std::vector; } KM3Detector::~KM3Detector() { // newgeant sxp.Finalize(); delete allCathods; for (size_t i = 0; i < allOMs->size(); i++) { (*allOMs)[i]->CathodsIDs->clear(); delete (*allOMs)[i]->CathodsIDs; free((*allOMs)[i]); } allOMs->clear(); delete allOMs; for (size_t i = 0; i < allStoreys->size(); i++) { (*allStoreys)[i]->BenthosIDs->clear(); delete (*allStoreys)[i]->BenthosIDs; free((*allStoreys)[i]); } allStoreys->clear(); delete allStoreys; } void KM3Detector::FindDetectorRadius() { G4double absdetectorRadius = 0.0; lowestStorey = 50.0 * m; //this initial value has no particular reason! I just wanted to make sure that it always will be higher than the actual lowest DOM highestStorey = 0.0; outerStorey = 0.0; detectorMaxz = 0.0; detectorMaxRho = 0.0; //std::vector x_storeys; //std::vector y_storeys; G4double x_storeys = 0.0; G4double y_storeys = 0.0; for (size_t isto = 0; isto < allStoreys->size(); isto++) { G4ThreeVector pos = (*(allStoreys))[isto]->position; // x_storeys.push_back(pos[0]); // y_storeys.push_back(pos[1]); x_storeys+=pos[0]; y_storeys+=pos[1]; if (pos[2] < lowestStorey) lowestStorey = pos[2]; if (pos[2] > highestStorey) highestStorey = pos[2]; } G4double x_center = x_storeys/allStoreys->size(); G4double y_center = y_storeys/allStoreys->size(); G4double z_center = lowestStorey + (highestStorey-lowestStorey)/2.; detectorCenter = G4ThreeVector(x_center, y_center, z_center); for (size_t isto = 0; isto < allStoreys->size(); isto++) { G4ThreeVector pos = (*(allStoreys))[isto]->position; G4double radius = (*(allStoreys))[isto]->radius; G4ThreeVector pos_shifted = pos-detectorCenter; G4double dist = pos_shifted.mag() + radius; if (absdetectorRadius < dist) absdetectorRadius = dist; G4double distRho = sqrt(pos_shifted[0] * pos_shifted[0] + pos_shifted[1] * pos_shifted[1]) + radius; if (outerStorey < distRho) outerStorey = distRho; } //detectorCenter = G4ThreeVector(0.0, 0.0, lowestStorey + (highestStorey-lowestStorey)/2.); detectorMaxRho = outerStorey + MaxAbsDist; detectorRadius = MaxAbsDist + absdetectorRadius; bottomPosition = 0.0 * m; detectorMaxz = highestStorey + MaxAbsDist; if(verbosity>0){ G4cout << "Detector radius (m) "<PutFromDetector(detectorCenter, detectorMaxRho, detectorMaxz, bottomPosition); //cout<<"FLAG TEST DETECTOR DIMENSIONS"<AddElement(elementO, 0.481); Crust->AddElement(elementSi, 0.277); Crust->AddElement(elementAl, 0.081); Crust->AddElement(elementFe, 0.050); Crust->AddElement(elementCa, 0.036); Crust->AddElement(elementNa, 0.028); Crust->AddElement(elementK, 0.026); Crust->AddElement(elementMg, 0.021); // WATER---here is added the depedence of the salinity and density on the // depth (UNESCO)-- // At the surface the seawater density is approximately 1.029gr/cm3 // (E.G.Anassontzis, "The NESTOR Site", // Proceedings of the 3rd NESTOR International Workshop, October 19-21, 1993, // Pylos Greece, pp614-630) // based on the composition of ocean seawater for 35o/oo // www.soest.hawaii.edu/oceanography we have G4double abundanceNa = 469.0e-3 * mole / kg; G4double abundanceMg = 52.8e-3 * mole / kg; G4double abundanceCa = 10.3e-3 * mole / kg; G4double abundanceK = 10.2e-3 * mole / kg; G4double abundanceCl = 545.9e-3 * mole / kg; G4double abundanceSO4 = 28.2e-3 * mole / kg; G4double abundanceHCO3 = 2.33e-3 * mole / kg; // The Mediterranean Salinity is about 40o/oo (gr/kgr of seawater) // we scale these abundances for 40o/oo salinity G4double scale = 40.0 / 35.0; abundanceNa *= scale; abundanceMg *= scale; abundanceCa *= scale; abundanceK *= scale; abundanceCl *= scale; abundanceSO4 *= scale; abundanceHCO3 *= scale; // we convert to gr/kgr of the solution (seawater) for each element or // composite abundanceNa *= weightNa; abundanceMg *= weightMg; abundanceCa *= weightCa; abundanceK *= weightK; abundanceCl *= weightCl; abundanceSO4 *= (weightS + 4 * weightO); abundanceHCO3 *= (weightH + weightC + 3 * weightO); G4double abundanceH2O = 1.0 - (abundanceNa + abundanceMg + abundanceCa + abundanceK + abundanceCl + abundanceSO4 + abundanceHCO3); //------------------the composites of sea //water------------------------------------ // if the material is not used to fill a specific logical volume the density // is not used G4Material *H2O = new G4Material("H2O", 1.0 * g / cm3, 2); H2O->AddElement(elementH, 2); H2O->AddElement(elementO, 1); G4Material *SO4 = new G4Material("SO4", 1.0 * g / cm3, 2); SO4->AddElement(elementS, 1); SO4->AddElement(elementO, 4); G4Material *HCO3 = new G4Material("HCO3", 1.0 * g / cm3, 3); HCO3->AddElement(elementH, 1); HCO3->AddElement(elementC, 1); HCO3->AddElement(elementO, 3); // calculate the density of sea water using the compressibility (Apostolis // Thesis appendix C) // and the detector depth (in meters : not water equivalent) G4double Compressibility = 4.29e-5 / bar; G4double gravity = 9.8 * m / (s * s); G4double surfaceDensity = 1.029 * g / cm3; G4double seawaterDensity = surfaceDensity * exp(gravity * surfaceDensity * Compressibility * detectorDepth); //----------------------------------------------------------------- G4Material *Water = new G4Material("Water", seawaterDensity, 8, kStateLiquid, 287.15 * kelvin, 1.0 * atmosphere); Water->AddMaterial(H2O, abundanceH2O); Water->AddMaterial(SO4, abundanceSO4); Water->AddMaterial(HCO3, abundanceHCO3); Water->AddElement(elementCl, abundanceCl); Water->AddElement(elementMg, abundanceMg); Water->AddElement(elementNa, abundanceNa); Water->AddElement(elementCa, abundanceCa); Water->AddElement(elementK, abundanceK); // GLASS and GELL change (27/5/2005) of the composition of borosilicate glass // according to // US Standards (http://www.udel.edu/chem/GlassShop/PhysicalProperties.htm) // 80.6% SiO2, 13.0% B2O3, 4.0% Na2O, 2.4% Al2O3 G4Material *materialSiO2 = new G4Material("SiO2", 1.0 * g / cm3, 2); materialSiO2->AddElement(elementSi, 1); materialSiO2->AddElement(elementO, 2); G4Material *materialB2O3 = new G4Material("B2O3", 1.0 * g / cm3, 2); materialB2O3->AddElement(elementB, 2); materialB2O3->AddElement(elementO, 3); G4Material *materialNa2O = new G4Material("Na2O", 1.0 * g / cm3, 2); materialNa2O->AddElement(elementNa, 2); materialNa2O->AddElement(elementO, 1); G4Material *materialAl2O3 = new G4Material("Al2O3", 1.0 * g / cm3, 2); materialAl2O3->AddElement(elementAl, 2); materialAl2O3->AddElement(elementO, 3); G4Material *Glass = new G4Material("Glass", 2.23 * g / cm3, 4, kStateSolid, 287.15 * kelvin, 1.0 * atmosphere); Glass->AddMaterial(materialSiO2, 0.806); Glass->AddMaterial(materialB2O3, 0.130); Glass->AddMaterial(materialNa2O, 0.040); Glass->AddMaterial(materialAl2O3, 0.024); // Silicone Gel Material. // density is 0.97g/cm3 // (http://www.wackersilicones.com/documents/techdatasheets/silgel612.pdf) // Polydimethylsiloxane polymeres -- (C2H6OSi)n G4Material *Gell = new G4Material("Gell", 0.97 * g / cm3, 4, kStateSolid, 287.15 * kelvin, 1.0 * atmosphere); Gell->AddElement(elementC, 2); Gell->AddElement(elementH, 6); Gell->AddElement(elementO, 1); Gell->AddElement(elementSi, 1); // AIR // ----------------------------------------------------------------------------- G4Material *Air = new G4Material("Air", 1.29e-03 * g / cm3, 2); Air->AddElement(elementN, .7); Air->AddElement(elementO, .3); // CATHOD // ------------------------------------------------------------------------------- G4Material *Cathod = new G4Material("Cathod", 22, 47.867 * g / mole, 4.507 * g / cm3, kStateSolid, 287.15 * kelvin, 1.0 * atmosphere); // G4cout<<*(G4Material::GetMaterialTable())<AddProperty("RINDEX", PPCKOV, RINDEX_WATER, NUMENTRIES); Properties_Water->AddProperty("ABSLENGTH", &ABSORPTION_WATER_E[0], &ABSORPTION_WATER[0], ABSORPTION_WATER_E.size()); Properties_Water->AddProperty("MIELENGTH", PPCKOV, SCATTER_WATER, NUMENTRIES); Properties_Water->AddConstProperty("MIEETA", MieModel); Water->SetMaterialPropertiesTable(Properties_Water); // GLASS ----------------------------------------------------------- G4MaterialPropertiesTable *Properties_Glass = new G4MaterialPropertiesTable(); Properties_Glass->AddProperty("ABSLENGTH", PPCKOV, ABSORPTION_GLASS, NUMENTRIES); Properties_Glass->AddProperty("RINDEX", PPCKOV, RINDEX_GLASS, NUMENTRIES); Glass->SetMaterialPropertiesTable(Properties_Glass); // GELL ----------------------------------------------------------- G4MaterialPropertiesTable *Properties_Gell = new G4MaterialPropertiesTable(); Properties_Gell->AddProperty("ABSLENGTH", PPCKOV, ABSORPTION_GELL, NUMENTRIES); Properties_Gell->AddProperty("RINDEX", PPCKOV, RINDEX_GELL, NUMENTRIES); Gell->SetMaterialPropertiesTable(Properties_Gell); // AIR ----------------------------------------------------------- // absorption length of air set to 1 meter to prevent trapping of a photon // inside air material G4MaterialPropertiesTable *Properties_Air = new G4MaterialPropertiesTable(); Properties_Air->AddProperty("ABSLENGTH", PPCKOV, ABSORPTION_AIR, NUMENTRIES); Properties_Air->AddProperty("RINDEX", PPCKOV, RINDEX_AIR, NUMENTRIES); Air->SetMaterialPropertiesTable(Properties_Air); // CATHOD ---------------------------------------------------------- G4MaterialPropertiesTable *Properties_Cath = new G4MaterialPropertiesTable(); Properties_Cath->AddProperty("ABSLENGTH", PPCKOV, ABSORPTION_CATH, NUMENTRIES); Properties_Cath->AddProperty("RINDEX", PPCKOV, RINDEX_CATH, NUMENTRIES); Properties_Cath->AddProperty("Q_EFF", PPCKOV, Q_EFF, NUMENTRIES); Properties_Cath->AddProperty("ANGULAR_ACCEPTANCE", COSANGLES, ACCEPTANCE, NUMENTRIES_ANGLEACC); Cathod->SetMaterialPropertiesTable(Properties_Cath); } #include "G4VoxelLimits.hh" G4int KM3Detector::TotalPMTEntities(const G4VPhysicalVolume *aPVolume, const G4GDMLParser* parser) const { static G4int Cathods = 0; static G4int Storeys = 0; static G4int OMs = 0; static G4AffineTransform AffineTrans; static G4RotationMatrix RotationMatr; static G4int Depth = 0; static G4int Hist[20]; static std::vector *aBenthosIDs; // static in order to load the benthos (OMs) static std::vector *aCathodsIDs; // static in order to load the cathods (Cathods) //G4cout <GetCopyNo()<<" "<GetName()<GetCopyNo(); Depth++; RotationMatr = RotationMatr * aPVolume->GetObjectRotationValue(); // if(aPVolume->GetName() == "CathodVolume_PV"){ //for newgeant add "_PV" at // the end of physical volume name if ((aPVolume->GetName()).contains("CathodVolume")) { // for newgeant add // "_PV" at the end of // physical volume name G4ThreeVector Position = AffineTrans.TransformPoint(aPVolume->GetObjectTranslation()); G4ThreeVector Direction = RotationMatr(G4ThreeVector(0.0, 0.0, 1.0)); G4Transform3D trans(RotationMatr, Position); // estimate cathod radius/////////////////////////////// G4double CathodRadius = 0.0; G4double CathodHeight = -1.0 * mm; // set default to begative, since it is // not applicable to spherical cathods if (aPVolume->GetLogicalVolume()->GetSolid()->GetEntityType() == G4String("G4Sphere")) { CathodRadius = ((G4Sphere *)aPVolume->GetLogicalVolume()->GetSolid()) ->GetOuterRadius(); G4double InnerRadius = ((G4Sphere *)aPVolume->GetLogicalVolume()->GetSolid()) ->GetInnerRadius(); if ((CathodRadius - InnerRadius) < 1.001 * mm) CathodRadius = 0.5 * (CathodRadius + InnerRadius); // applicable mainly to shell type cathods (EM) } else if (aPVolume->GetLogicalVolume()->GetSolid()->GetEntityType() == G4String("G4Tubs")) { CathodRadius = ((G4Tubs *)aPVolume->GetLogicalVolume()->GetSolid()) ->GetOuterRadius(); // applicable to thin tube cathods // (normal run) CathodHeight = ((G4Tubs *)aPVolume->GetLogicalVolume()->GetSolid()) ->GetZHalfLength(); CathodHeight *= 2.0; // full height } /////////////////////////////////////////// string cath = "CathodID_"; auto result = cath + std::to_string(Cathods); G4int detx_ID = parser->GetConstant(result); allCathods->addCathod(trans, Position, Direction, CathodRadius, CathodHeight, Depth - 1, detx_ID); for (G4int i = 1; i < Depth; i++){ allCathods->addToTree(Hist[i]); // G4cout << "Depth: "<< Depth << " Hist["<push_back(Cathods); Cathods++; } else { if ((aPVolume->GetName()).contains("OMVolume")) { // for newgeant add "_PV" // at the end of physical // volume name OMPositions *aOM = (OMPositions *)malloc(sizeof(OMPositions)); aOM->position = AffineTrans.TransformPoint(aPVolume->GetObjectTranslation()); aCathodsIDs = new std::vector; aOM->CathodsIDs = aCathodsIDs; // if OM is sphere then set the outer radius as radius, // if it is tubs then set the proper radius // else set the geometrical sum of the extend on the three axis //(maximum extend. Exact only for Boxes) if (aPVolume->GetLogicalVolume()->GetSolid()->GetEntityType() == G4String("G4Sphere")) { aOM->radius = ((G4Sphere *)aPVolume->GetLogicalVolume()->GetSolid()) ->GetOuterRadius(); } else if (aPVolume->GetLogicalVolume()->GetSolid()->GetEntityType() == G4String("G4Tubs")) { G4double zLength = ((G4Tubs *)aPVolume->GetLogicalVolume()->GetSolid()) ->GetZHalfLength(); G4double oRadius = ((G4Tubs *)aPVolume->GetLogicalVolume()->GetSolid()) ->GetOuterRadius(); aOM->radius = sqrt(zLength * zLength + oRadius * oRadius); } else { G4VoxelLimits voxelLimits; // Defaults to "infinite" limits. G4AffineTransform affineTransform; // no transform G4double xmin, xmax, ymin, ymax, zmin, zmax; aPVolume->GetLogicalVolume()->GetSolid()->CalculateExtent( kXAxis, voxelLimits, affineTransform, xmin, xmax); aPVolume->GetLogicalVolume()->GetSolid()->CalculateExtent( kYAxis, voxelLimits, affineTransform, ymin, ymax); aPVolume->GetLogicalVolume()->GetSolid()->CalculateExtent( kZAxis, voxelLimits, affineTransform, zmin, zmax); xmax = fmax(fabs(xmax), fabs(xmin)); ymax = fmax(fabs(ymax), fabs(ymin)); zmax = fmax(fabs(zmax), fabs(zmin)); aOM->radius = sqrt(xmax * xmax + ymax * ymax + zmax * zmax); } allOMs->push_back(aOM); aBenthosIDs->push_back(OMs); OMs++; } if ((aPVolume->GetName()) .contains("StoreyVolume")) { // for newgeant add "_PV" at the end // of physical volume name StoreysPositions *aStorey = (StoreysPositions *)malloc(sizeof(StoreysPositions)); aStorey->position = AffineTrans.TransformPoint(aPVolume->GetObjectTranslation()); aBenthosIDs = new std::vector; aStorey->BenthosIDs = aBenthosIDs; allStoreys->push_back(aStorey); Storeys++; } G4AffineTransform tempoaffine(aPVolume->GetObjectRotationValue().inverse(), aPVolume->GetObjectTranslation()); AffineTrans = tempoaffine * AffineTrans; for (G4int i = 0; i < aPVolume->GetLogicalVolume()->GetNoDaughters(); i++) { // the following is to fix new GDML that does not apply a copyNumber to // physical volumes aPVolume->GetLogicalVolume()->GetDaughter(i)->SetCopyNo(i); TotalPMTEntities(aPVolume->GetLogicalVolume()->GetDaughter(i), parser); } AffineTrans = tempoaffine.Inverse() * AffineTrans; } RotationMatr = RotationMatr * aPVolume->GetObjectRotationValue().inverse(); Depth--; return Cathods; } G4VPhysicalVolume *KM3Detector::Construct() { SetUpVariables(); ConstructMaterials(); //G4GeometryManager::GetInstance()->SetWorldMaximumExtent(1000.0 * m); G4GDMLParser parser; // newgeant parser.Read(Geometry_File); // newgeant fWorld = (G4VPhysicalVolume *)parser.GetWorldVolume(); // newgeant if (fWorld == 0) G4Exception("World volume not set properly check your setup selection " "criteria or GDML input!", "", FatalException, ""); G4int Cathods_number = TotalPMTEntities(fWorld, &parser); if(verbosity>0) G4cout << "Total Cathods " << Cathods_number << G4endl; // G4int History[10]={20,20,18,0,2,0,0,0,0,0}; // G4int dep=5; // G4cout <<"123456 "<(dep,History)<GetCathodId(dep,History)< *aLogicalStore; G4String aString1("CathodVolume"); G4String aString2("DeadVolume"); size_t theSize = G4LogicalVolumeStore::GetInstance()->size(); aLogicalStore = G4LogicalVolumeStore::GetInstance(); size_t PhysicalSize = G4PhysicalVolumeStore::GetInstance()->size(); std::vector *aPhysicalStore; aPhysicalStore = G4PhysicalVolumeStore::GetInstance(); G4VPhysicalVolume* aPhysicalVolume; for(size_t k = 0; k< PhysicalSize; k++){ aPhysicalVolume = (*aPhysicalStore)[k]; if(aPhysicalVolume->GetName().contains("CrustVolume")){ const G4ThreeVector vect = aPhysicalVolume->GetTranslation(); aLogicalVolume = aPhysicalVolume->GetLogicalVolume(); G4cout<<"Crust Volume Z position is: "<GetSolid(); G4cout<<"Sea bottom position is: "<GetZHalfLength()/m <GetName() << G4endl; G4cout << "Enter Color" << G4endl; scanf("%s", Colour); G4String colour(Colour); if (colour == G4String("none")) aLogicalVolume->SetVisAttributes(G4VisAttributes::Invisible); else { G4VisAttributes *MyVis; G4bool ok = true; if (colour == G4String("white")) MyVis = new G4VisAttributes(G4Colour::White()); else if (colour == G4String("gray") || colour == G4String("grey")) MyVis = new G4VisAttributes(G4Colour::Gray()); else if (colour == G4String("black")) MyVis = new G4VisAttributes(G4Colour::Black()); else if (colour == G4String("red")) MyVis = new G4VisAttributes(G4Colour::Red()); else if (colour == G4String("green")) MyVis = new G4VisAttributes(G4Colour::Green()); else if (colour == G4String("blue")) MyVis = new G4VisAttributes(G4Colour::Blue()); else if (colour == G4String("cyan")) MyVis = new G4VisAttributes(G4Colour::Cyan()); else if (colour == G4String("magenta")) MyVis = new G4VisAttributes(G4Colour::Magenta()); else if (colour == G4String("yellow")) MyVis = new G4VisAttributes(G4Colour::Yellow()); else ok = false; if (ok) { MyVis->SetVisibility(true); MyVis->SetLineStyle(MyVis->unbroken); MyVis->SetLineWidth(1.); // MyVis->SetForceWireFrame(true); MyVis->SetForceSolid(true); aLogicalVolume->SetVisAttributes(MyVis); } else aLogicalVolume->SetVisAttributes(G4VisAttributes::Invisible); } } } // fully adjustable benthos and storey linked list G4cout << "Total World Volume Entities= " << fWorld->GetLogicalVolume()->TotalVolumeEntities() << G4endl; // Next find from OM positions and radius in each store the storey radius for (size_t istorey = 0; istorey < allStoreys->size(); istorey++) { G4ThreeVector storeyposition = (*allStoreys)[istorey]->position; G4double MAXdist = 0.0; size_t OMnumber = (*allStoreys)[istorey]->BenthosIDs->size(); for (size_t iom = 0; iom < OMnumber; iom++) { G4int iOM = (*((*allStoreys)[istorey]->BenthosIDs))[iom]; G4ThreeVector OMposition = (*allOMs)[iOM]->position; G4double OMradius = (*allOMs)[iOM]->radius; G4double dist = (storeyposition - OMposition).mag() + OMradius; if (dist > MAXdist) MAXdist = dist; } (*allStoreys)[istorey]->radius = MAXdist; } // find detector radius and detector center from the Storeys FindDetectorRadius(); //--------Write the header of the outfile and the Cathods Position, Direction //and History Tree G4int nnn0 = 0; G4int nnn1 = 1; G4int nben = allCathods->GetNumberOfCathods(); // here in case of antares evt format I should write something general about // simulation //if (useANTARESformat) evtAANETWriter->ReadRunHeader(); if (useANTARESformat || useAANETformat) { ostringstream tagstream; tagstream << "KM3Sim " << Geometry_File; evtAANETWriter->AddHeaderInfo("detector", tagstream.str()); tagstream.str(""); //seed: progname ilux iseed k1 k2 //default ilux is 3 in HEP::CLHEP::setTheSeed (long seed, int lux=3) tagstream << "KM3Sim " << CLHEP::HepRandom::getTheSeed() ; evtAANETWriter->AddHeaderInfo("seed",tagstream.str()); tagstream.str(""); tagstream << "KM3Sim "<currentDateTime(); //tagstream << "KM3Sim "<AddHeaderInfo("simul",tagstream.str()); evtAANETWriter->WriteRunHeader(); } // find the total photocathod area on a OM G4int CaPerOM = (*allOMs)[0]->CathodsIDs->size(); TotCathodArea = CaPerOM * pi * allCathods->GetCathodRadius(0) * allCathods->GetCathodRadius(0); // this is valid only if at simulation // level (not EM or HA param) all cathods // have the same radius. Easy to change to // account for a detector with varius // cathod types // return the physical World return fWorld; } void KM3Detector::ConstructSDandField(){ //------------------------------------------------ // Sensitive detectors //------------------------------------------------ G4SDManager *SDman = G4SDManager::GetSDMpointer(); G4String MySDname = "mydetector1/MySD"; KM3SD *aMySD = new KM3SD(MySDname); aMySD->SetVerboseLevel(verbosity); aMySD->myStDetector = this; SDman->AddNewDetector(aMySD); G4LogicalVolume *aLogicalVolume; std::vector *aLogicalStore; G4String aString1("CathodVolume"); G4String aString2("DeadVolume"); size_t theSize = G4LogicalVolumeStore::GetInstance()->size(); aLogicalStore = G4LogicalVolumeStore::GetInstance(); for (size_t i = 0; i < theSize; i++) { aLogicalVolume = (*aLogicalStore)[i]; if (((aLogicalVolume->GetName()).contains(aString1)) || ((aLogicalVolume->GetName()).contains(aString2))) { aLogicalVolume->SetSensitiveDetector(aMySD); } } }