ch5/src/DetectorConstruction.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // gpaterno, October 2025
0027 //
0028 /// \file DetectorConstruction.cc
0029 /// \brief Implementation of the DetectorConstruction class
0030 //
0031 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0032
0033 #include "DetectorConstruction.hh"
0034 #include "SensitiveDetector.hh"
0035
0036 #include "G4RunManager.hh"
0037 #include "G4NistManager.hh"
0038 #include "G4Transform3D.hh"
0039 #include "G4Box.hh"
0040 #include "G4Tubs.hh"
0041 #include "G4Cons.hh"
0042 #include "G4Orb.hh"
0043 #include "G4Sphere.hh"
0044 #include "G4Trd.hh"
0045 #include "G4VSolid.hh"
0046 #include "G4LogicalVolume.hh"
0047 #include "G4SystemOfUnits.hh"
0048 #include "G4AnalysisManager.hh"
0049 #include "G4RegionStore.hh"
0050 #include "G4VisAttributes.hh"
0051 #include "G4SystemOfUnits.hh"
0052 #include "globals.hh"
0053
0054 #include "G4SubtractionSolid.hh"
0055
0056 #include "G4SDManager.hh"
0057 #include "G4VSensitiveDetector.hh"
0058 #include "G4MultiFunctionalDetector.hh"
0059 #include "G4PSDoseDeposit.hh"
0060 #include "G4PSEnergyDeposit.hh"
0061 #include "G4SDParticleFilter.hh"
0062
0063 #include "G4UniformMagField.hh"
0064 #include "G4FieldManager.hh"
0065 #include "G4TransportationManager.hh"
0066 #include <G4ChordFinder.hh>
0067
0068 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0069
0070 DetectorConstruction::DetectorConstruction()
0071 {
0072     //instantiate the messenger
0073     fMessenger = new DetectorConstructionMessenger(this);
0074 }
0075
0076 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0077
0078 G4VPhysicalVolume* DetectorConstruction::Construct()
0079 {
0080     G4cout << G4endl << "### DetectorConstruction::Construct() ###" << G4endl
0081            << G4endl;
0082
0083
0084     //sanity check
0085     if (fHybridSource) {
0086         G4cout << "This is a hybrid positron source!" << G4endl << G4endl;
0087     } else {
0088         fConverter = false;
0089         fRadiatorConverterSepDistance = 0.;
0090         G4cout << "This is a single-crystal positron source!" << G4endl
0091                << G4endl;
0092     }
0093
0094
0095     //check overlap option
0096     G4bool checkOverlaps = true;
0097
0098
0099     //Materials
0100     G4NistManager* nist = G4NistManager::Instance();
0101     G4Material* Silicon = nist->FindOrBuildMaterial("G4_Si");
0102     G4Material* PWO = nist->FindOrBuildMaterial("G4_PbWO4");
0103     G4Material* Diamond = nist->FindOrBuildMaterial("G4_C");
0104     G4Material* Tungsten = nist->FindOrBuildMaterial("G4_W");
0105     G4Material* Iridium = nist->FindOrBuildMaterial("G4_W");
0106     G4Material* Germanium = nist->FindOrBuildMaterial("G4_Ge");
0107     G4Element* elBi = nist->FindOrBuildElement("Bi");
0108     G4Element* elGe = nist->FindOrBuildElement("Ge");
0109     G4Element* elO = nist->FindOrBuildElement("O");
0110     G4Material* BGO = new G4Material("G4_BGO", 7.13*g/cm3, 3);
0111     BGO->AddElement(elBi, 0.671054);
0112     BGO->AddElement(elGe, 0.17482);
0113     BGO->AddElement(elO, 0.154126);
0114
0115     //Vacuum
0116     G4double z = 7.;
0117     G4double a = 14.007 * CLHEP::g/CLHEP::mole;
0118     G4double density = CLHEP::universe_mean_density;
0119     G4double pressure = 1.E-6 * 1.E-3 * CLHEP::bar; //10-6 mbar
0120     G4double temperature = 300. * CLHEP::kelvin; //300 K
0121     G4Material* Vacuum = new G4Material("Vacuum",
0122                                         z,
0123                                         a,
0124                                         density,
0125                                         kStateGas,
0126                                         temperature,
0127                                         pressure);
0128
0129
0130     // ------------------- World -------------------------
0131     G4Box* solidWorld = new G4Box("World", 3.*m, 3.*m, 20.*m);
0132
0133     G4LogicalVolume* logicWorld = new G4LogicalVolume(solidWorld,
0134                                                       Vacuum,
0135                                                       "World");
0136
0137     logicWorld->SetVisAttributes(G4VisAttributes::GetInvisible());
0138
0139     G4VPhysicalVolume* physWorld = new G4PVPlacement
0140                                    (0,              // no rotation
0141                                    G4ThreeVector(), // centre position
0142                                    logicWorld,      // its logical volume
0143                                    "World",         // its name
0144                                    0,               // its mother volume
0145                                    false,           // no boolean operation
0146                                    0,               // copy number
0147                                    checkOverlaps);  // overlaps checking
0148
0149
0150     // --------------- Crystal (Radiator) -------------------------
0151     //set visualization attributes
0152     G4VisAttributes* CrystalVisAttribute =
0153         new G4VisAttributes(G4Colour(0., 0., 1., 1.));
0154     CrystalVisAttribute->SetForceSolid(true);
0155
0156     //select crystal material
0157     if (fCrystalMaterialStr == "PWO") {
0158         fCrystalMaterial = PWO;
0159     } else if (fCrystalMaterialStr == "BGO") {
0160         fCrystalMaterial = BGO;
0161     } else if (fCrystalMaterialStr == "C") {
0162         fCrystalMaterial = Diamond;
0163     } else if (fCrystalMaterialStr == "W") {
0164         fCrystalMaterial = Tungsten;
0165     } else if (fCrystalMaterialStr == "Ir") {
0166         fCrystalMaterial = Iridium;
0167     } else if (fCrystalMaterialStr == "Ge") {
0168         fCrystalMaterial = Germanium;
0169     } else {
0170         fCrystalMaterial = Silicon;
0171     }
0172
0173     //Crystal solid and logic volumes
0174     G4Box* crystalSolid = new G4Box("Crystal",
0175                                     fCrystalSize.x()*0.5,
0176                                     fCrystalSize.y()*0.5,
0177                                     fCrystalSize.z()*0.5);
0178
0179     fCrystalLogic = new G4LogicalVolume(crystalSolid,
0180                                         fCrystalMaterial,
0181                                         "Crystal");
0182
0183     //Set actually or not the radiatior crystal
0184     fCrystalLogic->SetVisAttributes(CrystalVisAttribute);
0185
0186     //Crystal position
0187     G4double tollfCrystalZ = 0.*mm;
0188     if (fAngleX > 0 || fAngleY > 0) {
0189         tollfCrystalZ = 0.15*mm;
0190         //this allows us to tolerate misalignements up to 30 mrad on 1 cm thick crystals
0191     }
0192     fCrystalZ = -fCrystalSize.z()*0.5 - tollfCrystalZ;
0193     G4ThreeVector posCrystal = G4ThreeVector(0.*mm, 0.*mm, fCrystalZ);
0194
0195     //Crystal rotation angle (also the angle of crystal planes vs the beam)
0196     G4RotationMatrix* crystalRotationMatrix = new G4RotationMatrix;
0197     crystalRotationMatrix->rotateY(-fAngleX);
0198     crystalRotationMatrix->rotateX(-fAngleY);
0199
0200     //Crystal placement
0201     new G4PVPlacement(crystalRotationMatrix,
0202                       posCrystal,
0203                       fCrystalLogic,
0204                       "Crystal",
0205                       logicWorld,
0206                       false,
0207                       0,
0208                       checkOverlaps);
0209
0210     //Crystal region (necessary for the FastSim model)
0211     fCrystalRegion = new G4Region("Crystal");
0212     fCrystalRegion->AddRootLogicalVolume(fCrystalLogic);
0213
0214     //Print Crystal info
0215     G4cout << "Radiator Crystal set!" << G4endl;
0216     G4cout << "Crystal material: " << fCrystalMaterial->GetName() << G4endl;
0217     G4cout << "Crystal size: " << fCrystalSize.x()/mm
0218                                << "x" << fCrystalSize.y()/mm
0219            << "x" << fCrystalSize.z()/mm << " mm3" << G4endl;
0220     G4cout << "RadiatorZ: " << fCrystalZ/mm << " mm" << G4endl;
0221     G4cout << G4endl;
0222
0223
0224     // -------- volume for the magnetic field ------------
0225     if (fSetMagneticField &&
0226         fRadiatorConverterSepDistance > fFieldRegionLength + 0.1*mm) {
0227         G4double magFieldRegionZ = fRadiatorConverterSepDistance*0.5;
0228
0229         G4Tubs* tub1 = new G4Tubs("MFvolume",
0230                                   0.*cm,
0231                                   fFieldRegionLength*0.5,
0232                                   fFieldRegionLength*0.5,
0233                                   0.*deg,
0234                                   360.*deg);
0235
0236         G4RotationMatrix *xRot = new G4RotationMatrix;
0237         xRot->rotateX(90.*deg);
0238
0239         fMFlogic = new G4LogicalVolume(tub1,
0240                                        Vacuum,
0241                                        "MFvolume");
0242
0243         G4VisAttributes* MFregionVisAttribute =
0244             new G4VisAttributes(G4Colour(1., 0., 0., 0.35));
0245         MFregionVisAttribute->SetForceSolid(true);
0246         fMFlogic->SetVisAttributes(MFregionVisAttribute);
0247         new G4PVPlacement(xRot,
0248                           G4ThreeVector(0., 0., magFieldRegionZ),
0249                           fMFlogic,
0250                           "MFvolume",
0251                           logicWorld,
0252                           false,
0253                           0,
0254                           checkOverlaps);
0255
0256         G4cout << "Magnetic Field set!" << G4endl;
0257         G4cout << "Field value: " << fFieldValue/tesla << " T" << G4endl;
0258         G4cout << "Field Region Diameter: " << fFieldRegionLength/mm
0259                << " mm" << G4endl;
0260         G4cout << "Field Region Length: " << fFieldRegionLength/mm
0261                << " mm" << G4endl;
0262         G4cout << G4endl;
0263
0264         fSetCollimator = false;
0265     }
0266
0267
0268     // ------------------ Collimator ---------------------
0269     else if (fSetCollimator &&
0270              fRadiatorConverterSepDistance > fCollimatorThickness + 0.1*mm) {
0271         G4double CollimatorZ = fRadiatorCollimatorSepDistance +
0272                                fCollimatorThickness*0.5;
0273
0274         G4Box* outerCollimSolid = new G4Box("outerCollimSolid",
0275                                             fCollimatorSide*0.5,
0276                                             fCollimatorSide*0.5,
0277                                             fCollimatorThickness*0.5);
0278
0279         G4VSolid* innerCollimSolid;
0280         if (fCollimatorHole == "circular") {
0281             innerCollimSolid = new G4Tubs("innerCollimSolid",
0282                                           0.*cm,
0283                                           fCollimatorAperture*0.5,
0284                                           fCollimatorThickness*0.51,
0285                                           0.*deg,
0286                                           360.*deg);
0287         } else {
0288             innerCollimSolid = new G4Box("innerCollimSolid",
0289                                            fCollimatorAperture*0.5,
0290                                          fCollimatorAperture*0.5,
0291                                          fCollimatorThickness*0.6);
0292         }
0293
0294         G4SubtractionSolid* CollimSolid = new G4SubtractionSolid("Collimator",
0295                                                                  outerCollimSolid,
0296                                                                  innerCollimSolid);
0297
0298         fCollimatorLogic = new G4LogicalVolume(CollimSolid,
0299                                                Tungsten,
0300                                                "Collimator");
0301
0302         G4VisAttributes* CollimatorVisAttribute =
0303             new G4VisAttributes(G4Colour(0.5, 0.5, 0.5));
0304         CollimatorVisAttribute->SetForceSolid(true);
0305         fCollimatorLogic->SetVisAttributes(CollimatorVisAttribute);
0306         new G4PVPlacement(0,
0307                           G4ThreeVector(0., 0., CollimatorZ),
0308                           fCollimatorLogic,
0309                           "Collimator",
0310                           logicWorld,
0311                           false,
0312                           0,
0313                           checkOverlaps);
0314
0315         G4cout << "Collimator set!" << G4endl;
0316         G4cout << "Collimator aperture: " << fCollimatorAperture/mm
0317                << " mm" << G4endl;
0318         G4cout << "Collimator hole shape: " << fCollimatorHole << G4endl;
0319         G4cout << "Collimator thickness: " << fCollimatorThickness/mm
0320                << " mm" << G4endl;
0321         G4cout << "CollimatorZ: " << CollimatorZ/mm << " mm" << G4endl;
0322         G4cout << G4endl;
0323     }
0324
0325
0326     // -------- Converter (Target) randomly oriented -----
0327     //visualization attributes
0328     G4VisAttributes* ConverterVisAttribute =
0329         new G4VisAttributes(G4Colour(0., 0.2, 0.8, 0.8));
0330
0331     //select Converter material
0332     G4Material* sphereMaterial;
0333     if (fConverterMaterialStr == "PWO") {
0334         fConverterMaterial = PWO;
0335         sphereMaterial = PWO;
0336     } else if (fConverterMaterialStr == "BGO") {
0337         fConverterMaterial = BGO;
0338         sphereMaterial = BGO;
0339     } else if (fConverterMaterialStr == "Ir") {
0340         fConverterMaterial = Iridium;
0341         sphereMaterial = Iridium;
0342     } else {
0343         fConverterMaterial = Tungsten; //defualt is W
0344         sphereMaterial = Tungsten;
0345     }
0346
0347     if (fGranularConverter) {
0348         fConverterMaterial = Vacuum;
0349         //Converter will be a vacuum box filled with spheres (of fConverterMaterial)
0350         ConverterVisAttribute->SetForceSolid(false);
0351     } else {
0352         ConverterVisAttribute->SetForceSolid(true);
0353     }
0354
0355     //set Converter size
0356     G4double ConverterWidth = fConverterSize.x();
0357     G4double ConverterHeight = fConverterSize.y();
0358     G4double ConverterThickness = fConverterSize.z();
0359
0360     //Converter position
0361     G4double toll = fVirtualDetectorSize.z();
0362     if (fRadiatorConverterSepDistance > 0) {
0363         toll = 2.*fVirtualDetectorSize.z();
0364     }
0365     fConverterZ = fRadiatorConverterSepDistance + ConverterThickness*0.5 +
0366                   fVirtualDetectorSize.z() + toll;
0367     G4ThreeVector posConverter = G4ThreeVector(0.*mm, 0.*mm, fConverterZ);
0368
0369     //Converter volume
0370     G4Box* ConverterSolid = new G4Box("Converter",
0371                                       ConverterWidth*0.5,
0372                                       ConverterHeight*0.5,
0373                                       ConverterThickness*0.5);
0374
0375     fConverterLogic = new G4LogicalVolume(ConverterSolid,
0376                                           fConverterMaterial,
0377                                           "Converter");
0378
0379     fConverterLogic->SetVisAttributes(ConverterVisAttribute);
0380
0381     if (fConverter) {
0382         new G4PVPlacement(0,
0383                           posConverter,
0384                           fConverterLogic,
0385                           "Converter",
0386                           logicWorld,
0387                           false,
0388                           0,
0389                           checkOverlaps);
0390
0391         //print Converter info
0392         G4cout << "RadiatorConverterSepDistance: " << fRadiatorConverterSepDistance
0393                << " mm " << G4endl;
0394         G4cout << "ConverterZ: " << fConverterZ/mm << " mm" << G4endl;
0395         G4cout << "Converter material: " << fConverterMaterial->GetName() << G4endl;
0396         G4cout << "Converter size: " << ConverterWidth/mm
0397                << "x" << ConverterHeight/mm
0398                << "x" << ConverterThickness/mm << " mm3" << G4endl;
0399         G4cout << G4endl;
0400
0401         //Is it a glanular converter?
0402         if (fGranularConverter) {
0403             G4cout << "Granular Converter set!" << G4endl;
0404             G4cout << "sphere material: " << sphereMaterial->GetName() << G4endl;
0405             G4cout << "sphere radius: " << fSphereRadius/mm << " mm" << G4endl;
0406
0407             G4int numLayers = G4int(ConverterThickness/(std::sqrt(2.)*fSphereRadius));
0408             G4int NintX = G4int(ConverterWidth/(2.*fSphereRadius)) - 0;
0409             G4int NintY = G4int(ConverterHeight/(2.*fSphereRadius)) - 0;
0410             G4double spacing = 2*fSphereRadius;
0411
0412             //numLayers = 10;
0413             //NintX = NintY = 10;
0414
0415             G4cout << "numLayers: " << numLayers
0416                    << ", NintX: " << NintX
0417                    << ", NintY: " << NintY << G4endl;
0418
0419             G4Sphere* sphereSolid = new G4Sphere("Sphere",
0420                                                  0,
0421                                                  fSphereRadius,
0422                                                  0,
0423                                                  360*deg,
0424                                                  0,
0425                                                  180*deg);
0426
0427             G4VisAttributes* sphereVisAtt =
0428                 new G4VisAttributes(G4Colour(0., 0., 1., 0.7));
0429             sphereVisAtt->SetForceSolid(true);
0430
0431             G4int k = 0;
0432             for (int layer = 0; layer < numLayers; layer++) {
0433                 //alternating number of spheres in each layer so as
0434                 //they are even for odd layers and odd for even layers.
0435                 G4int numSpheresX = (layer % 2 == 0) ?
0436                     (NintX % 2 == 0) ? NintX-1 : NintX :
0437                     (NintX % 2 == 0) ? NintX : NintX-1;
0438                 G4int numSpheresY = (layer % 2 == 0) ?
0439                     (NintY % 2 == 0) ? NintY-1 : NintY :
0440                     (NintY % 2 == 0) ? NintY : NintY-1;
0441
0442                 //calculate the total width and height of the sphere arrangement in each layer
0443                 G4double totalWidth = numSpheresX*(2.*fSphereRadius);
0444                 G4double totalHeight = numSpheresY*(2.*fSphereRadius);
0445
0446                 //calculate the starting position for the first sphere in each layer
0447                 G4double startX = -totalWidth/2. + fSphereRadius;
0448                 G4double startY = -totalHeight/2. + fSphereRadius;
0449
0450                 //sphere positioning
0451                 for (int i = 0; i < numSpheresX; ++i) {
0452                     for (int j = 0; j < numSpheresY; ++j) {
0453                         std::string sphereName = "Sphere_" + std::to_string(layer) +
0454                                                        "_" + std::to_string(i) +
0455                                                        "_" + std::to_string(j);
0456
0457                         fSphereLogic[k] = new G4LogicalVolume(sphereSolid,
0458                                                               sphereMaterial,
0459                                                               sphereName);
0460                         fSphereLogic[k]->SetVisAttributes(sphereVisAtt);
0461                         fScoringVolume.push_back(fSphereLogic[k]);
0462
0463                         G4double xPosition = startX + i*spacing;
0464                         G4double yPosition = startY + j*spacing;
0465                         G4double zPosition = -ConverterThickness/2. + fSphereRadius
0466                                              + layer*std::sqrt(2)*fSphereRadius;
0467                         G4ThreeVector position(xPosition, yPosition, zPosition);
0468
0469                         new G4PVPlacement(0,
0470                                           position,
0471                                           fSphereLogic[k],
0472                                           sphereName,
0473                                           fConverterLogic,
0474                                           false,
0475                                           k,
0476                                           checkOverlaps);
0477
0478                         k++;
0479                     }
0480                 }
0481             }
0482
0483             fNSpheres = k;
0484             G4cout << "Nspheres positioned: " << fNSpheres << G4endl << G4endl;
0485
0486         }
0487
0488     }
0489
0490
0491     // --------------- virtual Detectors -----------------
0492     //position
0493     G4double VirtualDetector0Z = fConverterZ - ConverterThickness*0.5
0494                                              - fVirtualDetectorSize.z()*0.5;
0495     G4ThreeVector posVirtualDetector0 = G4ThreeVector(0, 0, VirtualDetector0Z);
0496     G4ThreeVector frontVirtualDetector0 = G4ThreeVector(0, 0, VirtualDetector0Z -
0497                                                               fVirtualDetectorSize.z()*0.5);
0498     G4cout << "VirtualDetector0Z: " << VirtualDetector0Z/mm << " mm" << G4endl;
0499     fVirtualDetectorPositionVector.push_back(frontVirtualDetector0);
0500
0501     G4double VirtualDetector1Z = fConverterZ + ConverterThickness*0.5
0502                                              + fVirtualDetectorSize.z()*0.5;
0503     G4ThreeVector posVirtualDetector1 = G4ThreeVector(0, 0, VirtualDetector1Z);
0504     G4ThreeVector frontVirtualDetector1 = G4ThreeVector(0, 0, VirtualDetector1Z -
0505                                                               fVirtualDetectorSize.z()*0.5);
0506     if (fHybridSource) {
0507         G4cout << "VirtualDetector1Z: " << VirtualDetector1Z/mm << " mm" << G4endl;
0508         fVirtualDetectorPositionVector.push_back(frontVirtualDetector1);
0509     }
0510
0511     G4double VirtualDetector2Z = fVirtualDetectorSize.z()*0.5;
0512     G4ThreeVector posVirtualDetector2 = G4ThreeVector(0, 0, VirtualDetector2Z);
0513     G4ThreeVector frontVirtualDetector2 = G4ThreeVector(0, 0, VirtualDetector2Z -
0514                                                               fVirtualDetectorSize.z()*0.5);
0515     if (fRadiatorConverterSepDistance > 0) {
0516         G4cout << "VirtualDetector2Z: " << VirtualDetector2Z/mm << " mm" << G4endl;
0517         fVirtualDetectorPositionVector.push_back(frontVirtualDetector2);
0518     }
0519
0520     //virtual Detector volume
0521     G4Box* VirtualDetectorSolid = new G4Box("VirtualDetector",
0522                                             fVirtualDetectorSize.x()*0.5,
0523                                             fVirtualDetectorSize.y()*0.5,
0524                                             fVirtualDetectorSize.z()*0.5);
0525
0526     fVirtualDetectorLogic0 = new G4LogicalVolume(VirtualDetectorSolid,
0527                                                  Vacuum,
0528                                                  "VirtualDetector0");
0529
0530     fVirtualDetectorLogic1 = new G4LogicalVolume(VirtualDetectorSolid,
0531                                                  Vacuum,
0532                                                  "VirtualDetector1");
0533
0534     fVirtualDetectorLogic2 = new G4LogicalVolume(VirtualDetectorSolid,
0535                                                  Vacuum,
0536                                                  "VirtualDetector2");
0537
0538     G4VisAttributes* VirtualDetectorVisAttribute =
0539         new G4VisAttributes(G4Colour(1., 1., 1.));
0540     VirtualDetectorVisAttribute->SetForceSolid(false);
0541     fVirtualDetectorLogic0->SetVisAttributes(VirtualDetectorVisAttribute);
0542     fVirtualDetectorLogic1->SetVisAttributes(VirtualDetectorVisAttribute);
0543     fVirtualDetectorLogic2->SetVisAttributes(VirtualDetectorVisAttribute);
0544
0545     new G4PVPlacement(0,
0546                       posVirtualDetector0,
0547                       fVirtualDetectorLogic0,
0548                       "VirtualDetector0",
0549                       logicWorld,
0550                       false,
0551                       0,
0552                       checkOverlaps);
0553
0554     if (fHybridSource) {
0555         new G4PVPlacement(0,
0556                           posVirtualDetector1,
0557                           fVirtualDetectorLogic1,
0558                           "VirtualDetector1",
0559                           logicWorld,
0560                           false,
0561                           1,
0562                           checkOverlaps);
0563     }
0564
0565     if (fRadiatorConverterSepDistance > 0) {
0566         new G4PVPlacement(0,
0567                           posVirtualDetector2,
0568                           fVirtualDetectorLogic2,
0569                           "VirtualDetector2",
0570                           logicWorld,
0571                           false,
0572                           2,
0573                           checkOverlaps);
0574     }
0575
0576
0577     //always return the physical World
0578     return physWorld;
0579 }
0580
0581 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0582
0583 void DetectorConstruction::ConstructSDandField()
0584 {
0585     //activate OC effects through FastSim model
0586     if (fActivateOCeffects) {
0587         G4RegionStore* regionStore = G4RegionStore::GetInstance();
0588         G4Region* RegionCh = regionStore->GetRegion("Crystal");
0589
0590         G4ChannelingFastSimModel* ChannelingModel =
0591             new G4ChannelingFastSimModel("ChannelingModel", RegionCh);
0592
0593         ChannelingModel->Input(fCrystalLogic->GetMaterial(), fLattice, fPotentialPath);
0594         ChannelingModel->GetCrystalData()->SetBendingAngle(fBendingAngle, fCrystalLogic);
0595
0596         G4double fParticleLEth = 1.*GeV; //deafult 200.*MeV (5.*GeV -> much faster)
0597         G4double fLindhardAngles = 10; //default 100
0598         ChannelingModel->SetLowKineticEnergyLimit(fParticleLEth, "e-");
0599         ChannelingModel->SetLowKineticEnergyLimit(fParticleLEth, "e+");
0600         ChannelingModel->SetLindhardAngleNumberHighLimit(fLindhardAngles, "e-");
0601         ChannelingModel->SetLindhardAngleNumberHighLimit(fLindhardAngles, "e+");
0602
0603         G4cout << G4endl;
0604         G4cout << "Oriented Crystal effects set through FastSim model" << G4endl;
0605         G4cout << "Crystal bending angle: " << fBendingAngle << " rad" << G4endl;
0606         G4cout << "Crystal Lattice: " << fLattice << G4endl;
0607         G4cout << "Crystal AngleX: " << fAngleX << " rad" << G4endl;
0608         G4cout << "Crystal AngleY: " << fAngleY << " rad" << G4endl;
0609         G4cout << "fParticleLEth: " << fParticleLEth/MeV << " MeV" << G4endl;
0610         G4cout << "fLindhardAngles: " << fLindhardAngles << G4endl;
0611         G4cout << "ActivateRadiationModel: " << fActivateRadiationModel << G4endl;
0612
0613         if (fActivateRadiationModel) {
0614             ChannelingModel->RadiationModelActivate();
0615             G4int fSamplingPhotonsNumber = 150; //default 150
0616             G4int fNSmallTrajectorySteps = 10000; //default 10000
0617             G4double fRadiactionAngleFactor = 4.; //deafult 4
0618             G4double fSinglePhotonRadProbLimit = 0.25; //default 0.25
0619             G4double fLEthreshold = 1.*MeV;
0620             ChannelingModel->GetRadiationModel()->
0621                 SetSamplingPhotonsNumber(fSamplingPhotonsNumber);
0622             ChannelingModel->GetRadiationModel()->
0623                 SetNSmallTrajectorySteps(fNSmallTrajectorySteps);
0624             ChannelingModel->GetRadiationModel()->
0625                 SetRadiationAngleFactor(fRadiactionAngleFactor);
0626             ChannelingModel->GetRadiationModel()
0627                 ->SetSinglePhotonRadiationProbabilityLimit(fSinglePhotonRadProbLimit);
0628             ChannelingModel->GetRadiationModel()
0629                 ->SetSpectrumEnergyRange(fLEthreshold, 20.*GeV, 100);
0630
0631             G4cout << "SamplingPhotonsNumber: "
0632                    << fSamplingPhotonsNumber << G4endl;
0633             G4cout << "NSmallTrajectorySteps: "
0634                    << fNSmallTrajectorySteps << G4endl;
0635             G4cout << "fRadiactionAngleFactor: "
0636                    << fRadiactionAngleFactor << G4endl;
0637             G4cout << "fSinglePhotonRadProbLimit: "
0638                    << fSinglePhotonRadProbLimit << G4endl;
0639             G4cout << "Low Eenergy threshold to emit photons and record their energy: "
0640                    << fLEthreshold/MeV << " MeV" << G4endl << G4endl;
0641         } else {
0642             G4cout << G4endl;
0643         }
0644     }
0645
0646
0647     //built-in Edep Scorer in the Radiator Crystal
0648     G4MultiFunctionalDetector* multisd = new G4MultiFunctionalDetector("multisd");
0649     G4VPrimitiveScorer* edepscorer = new G4PSEnergyDeposit("edep");
0650     multisd->RegisterPrimitive(edepscorer);
0651     SetSensitiveDetector(fCrystalLogic->GetName(), multisd);
0652     G4SDManager::GetSDMpointer()->AddNewDetector(multisd);
0653
0654
0655     //Sensitive Volumes (Virtual Detectors)
0656     G4VSensitiveDetector* vDetector = new SensitiveDetector("det");
0657     G4SDManager::GetSDMpointer()->AddNewDetector(vDetector);
0658     fVirtualDetectorLogic0->SetSensitiveDetector(vDetector);
0659     if (fHybridSource) {
0660         fVirtualDetectorLogic1->SetSensitiveDetector(vDetector);
0661     }
0662     if (fRadiatorConverterSepDistance > 0) {
0663         fVirtualDetectorLogic2->SetSensitiveDetector(vDetector);
0664     }
0665
0666
0667     //Magnetic field
0668     if (fSetMagneticField & fHybridSource) {
0669         G4UniformMagField* myField =
0670             new G4UniformMagField(G4ThreeVector(0., fFieldValue, 0.));
0671         G4FieldManager* localfieldMgr = new G4FieldManager(myField);
0672         localfieldMgr->CreateChordFinder(myField);
0673         fMFlogic->SetFieldManager(localfieldMgr, true);
0674     }
0675
0676
0677     G4cout << "### End of DetectorConstruction ###" << G4endl << G4endl << G4endl;
0678 }
0679
0680 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0681