ch2/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 /// \file DetectorConstruction.cc
0027 /// \brief Implementation of the DetectorConstruction class
0028 //
0029 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0030
0031 #include "DetectorConstruction.hh"
0032
0033 #include "G4RunManager.hh"
0034 #include "G4NistManager.hh"
0035 #include "G4Box.hh"
0036 #include "G4LogicalVolume.hh"
0037 #include "G4RegionStore.hh"
0038 #include "G4VisAttributes.hh"
0039 #include "PrimaryGeneratorAction.hh"
0040
0041 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0042
0043 DetectorConstruction::DetectorConstruction()
0044 {
0045     //instantiate the messenger
0046     fMessenger = new DetectorConstructionMessenger(this);
0047
0048     //Crystal size
0049     fCrystalSize.setX(20*CLHEP::mm);
0050     fCrystalSize.setY(20*CLHEP::mm);
0051     fCrystalSize.setZ(0.0305*CLHEP::mm);
0052
0053     //Crystal planes or axes considered
0054     fLattice = "(111)";
0055
0056     //Detector size
0057     fDetectorSize.setX(10*CLHEP::cm);
0058     fDetectorSize.setY(10*CLHEP::cm);
0059     fDetectorSize.setZ(0.1*CLHEP::mm);
0060 }
0061
0062 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0063
0064 G4VPhysicalVolume* DetectorConstruction::Construct()
0065 {
0066     //Check overlap option
0067     G4bool checkOverlaps = true;
0068
0069     //Materials
0070     G4NistManager* nist = G4NistManager::Instance();
0071     G4Material* world_mat = nist->FindOrBuildMaterial("G4_Galactic");
0072     G4Material* silicon = nist->FindOrBuildMaterial("G4_Si");
0073
0074     //World
0075     G4Box* solidWorld = new G4Box("World", 0.2*CLHEP::m, 0.2*CLHEP::m, 10.*CLHEP::m);
0076     G4LogicalVolume* logicWorld = new G4LogicalVolume(solidWorld, world_mat, "World");
0077     G4VPhysicalVolume* physWorld = new G4PVPlacement
0078                                     (0,                     // no rotation
0079                                      G4ThreeVector(),       // centre position
0080                                      logicWorld,            // its logical volume
0081                                      "World",               // its name
0082                                      0,                     // its mother volume
0083                                      false,                 // no boolean operation
0084                                      0,                     // copy number
0085                                      checkOverlaps);        // overlaps checking
0086     logicWorld->SetVisAttributes(G4VisAttributes::GetInvisible());
0087
0088
0089     // --------------- Crystal ------------------------------------
0090
0091     //Select crystal material
0092     fCrystalMaterial = nist->FindOrBuildMaterial(fCrystalMaterialStr);
0093
0094     //Crystal rotation angle (also the angle of crystal planes vs the beam)
0095     G4RotationMatrix* crystalRotationMatrix = new G4RotationMatrix;
0096     crystalRotationMatrix->rotateY(-fAngleX);
0097     crystalRotationMatrix->rotateX(-fAngleY);
0098
0099     //Setting crystal position
0100     G4ThreeVector posCrystal = G4ThreeVector(0., 0., fCrystalSize.z()/2.);
0101
0102     //crystal volume
0103     G4Box* solidCrystal = new G4Box("Crystal",
0104                                     fCrystalSize.x()/2,
0105                                     fCrystalSize.y()/2,
0106                                     fCrystalSize.z()/2.);
0107
0108     fLogicCrystal = new G4LogicalVolume(solidCrystal,
0109                                         fCrystalMaterial,
0110                                         "Crystal");
0111     new G4PVPlacement(crystalRotationMatrix,
0112                       posCrystal,
0113                       fLogicCrystal,
0114                       "Crystal",
0115                       logicWorld,
0116                       false,
0117                       0,
0118                       checkOverlaps);
0119
0120     if (fActivateChannelingModel)
0121     {
0122       //crystal region (necessary for the FastSim model)
0123       G4Region* regionCh = new G4Region("Crystal");
0124       regionCh->AddRootLogicalVolume(fLogicCrystal);
0125     }
0126
0127     //visualization attributes
0128     G4VisAttributes* crystalVisAttribute =
0129         new G4VisAttributes(G4Colour(1., 0., 0.));
0130     crystalVisAttribute->SetForceSolid(true);
0131     fLogicCrystal->SetVisAttributes(crystalVisAttribute);
0132
0133     //print Crystal info
0134     G4cout << "Crystal material: " << fCrystalMaterial->GetName() << G4endl;
0135     G4cout << "Crystal size: " << fCrystalSize.x()/CLHEP::mm
0136            << "x" << fCrystalSize.y()/CLHEP::mm
0137            << "x" << fCrystalSize.z()/CLHEP::mm << " mm3" << G4endl;
0138     if (fActivateChannelingModel)
0139     {
0140         G4cout << "G4ChannelingFastSimModel activated" << G4endl;
0141         G4cout << "Crystal bending angle: " << fBendingAngle << " rad" << G4endl;
0142         G4cout << "Crystal Lattice: " << fLattice << G4endl;
0143         G4cout << "Crystal angleX: " << fAngleX << " rad" << G4endl;
0144         G4cout << "Crystal angleY: " << fAngleY << " rad" << G4endl;
0145         G4cout << "ActivateRadiationModel: " << fActivateRadiationModel << G4endl;
0146
0147         if (fCrystallineUndulatorAmplitude > DBL_EPSILON &&
0148             fCrystallineUndulatorPeriod > DBL_EPSILON) {
0149             G4cout << "Crystalline undulator activated: " << G4endl;
0150             G4cout << "undulator amplitude: "
0151                    << fCrystallineUndulatorAmplitude/CLHEP::nm
0152                    << " nm" << G4endl;
0153             G4cout << "undulator period: "
0154                    << fCrystallineUndulatorPeriod/CLHEP::mm
0155                    << " mm" << G4endl;
0156             G4cout << "undulator phase: "
0157                    << fCrystallineUndulatorPhase
0158                    << " rad" << G4endl;
0159         }
0160
0161         G4cout << G4endl;
0162     }
0163     else
0164     {
0165         G4cout << "G4ChannelingFastSimModel is not activated" << G4endl << G4endl;
0166     }
0167
0168     // --------------- Detector -----------------------------------
0169     //Setting detector position
0170     G4ThreeVector posDetector =
0171             G4ThreeVector(0, 0, fDetectorFrontPosZ+fDetectorSize.z()/2.);
0172
0173     //particle detector volume
0174     G4Box* detector = new G4Box("Detector",
0175                                 fDetectorSize.x()/2,
0176                                 fDetectorSize.y()/2,
0177                                 fDetectorSize.z()/2.);
0178
0179     G4LogicalVolume* logicDetector = new G4LogicalVolume(detector,
0180                                                          silicon,
0181                                                          "Detector");
0182     new G4PVPlacement(0,
0183                       posDetector,
0184                       logicDetector,
0185                       "Detector",
0186                       logicWorld,
0187                       false,
0188                       0,
0189                       checkOverlaps);
0190
0191     //visualization attributes
0192     G4VisAttributes* detectorVisAttribute =
0193         new G4VisAttributes(G4Colour(1., 0., 1., 0.3));
0194     detectorVisAttribute->SetForceSolid(true);
0195     logicDetector->SetVisAttributes(detectorVisAttribute);
0196
0197     //always return the physical World
0198     return physWorld;
0199 }
0200
0201 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0202
0203 void DetectorConstruction::ConstructSDandField()
0204 {
0205   if (fActivateChannelingModel)
0206   {
0207     // --------------- fast simulation ----------------------------
0208     //extract the region of the crystal from the store
0209     G4RegionStore* regionStore = G4RegionStore::GetInstance();
0210     G4Region* regionCh = regionStore->GetRegion("Crystal");
0211
0212     //create the channeling model for this region
0213     G4ChannelingFastSimModel* channelingModel =
0214         new G4ChannelingFastSimModel("ChannelingModel", regionCh);
0215
0216     //activate the channeling model
0217     channelingModel->Input(fCrystalMaterial, fLattice, fPotentialPath);
0218     //setting bending angle of the crystal planes (default is 0)
0219     channelingModel->GetCrystalData()->SetBendingAngle(fBendingAngle, fLogicCrystal);
0220
0221     //setting crystalline undulator parameters
0222     //NOTE: they are incompatible with a bent crystal
0223     if (fCrystallineUndulatorAmplitude > DBL_EPSILON &&
0224         fCrystallineUndulatorPeriod > DBL_EPSILON)
0225     {
0226         channelingModel->GetCrystalData()->SetCrystallineUndulatorParameters(
0227             fCrystallineUndulatorAmplitude,
0228             fCrystallineUndulatorPeriod,
0229             fCrystallineUndulatorPhase,
0230             fLogicCrystal);
0231     }
0232
0233     /*
0234     Set the multiple of critical channeling angles (Lindhard angles) which defines
0235     the angular cut of the model (otherwise standard Geant4 is active):
0236     a number too low reduces the accuracy of the model;
0237     a number too high sometimes drastically reduces the simulation speed.
0238     The default value is 100, while for many problems 10-20 is ok.
0239     CAUTION: If you set this value to 1 meaning the angular cut = 1*Lindhard angle,
0240     this will cut off the physics of overbarrier motion, still important even if
0241     the particle is not in channeling. This will provide incorrect results.
0242     */
0243     channelingModel->SetDefaultLindhardAngleNumberHighLimit(fLindhardAngles);
0244     //you may set a particular limit for a certain particle type
0245     //(has a priority vs default):
0246     channelingModel->SetLindhardAngleNumberHighLimit(fLindhardAnglesProton,"proton");
0247     channelingModel->SetLindhardAngleNumberHighLimit(fLindhardAnglesAntiproton,
0248                                                      "anti_proton");
0249     channelingModel->SetLindhardAngleNumberHighLimit(fLindhardAnglesPiPlus,"pi+");
0250     channelingModel->SetLindhardAngleNumberHighLimit(fLindhardAnglesPiMinus,"pi-");
0251     channelingModel->SetLindhardAngleNumberHighLimit(fLindhardAnglesPositron,"e+");
0252     channelingModel->SetLindhardAngleNumberHighLimit(fLindhardAnglesElectron,"e-");
0253     channelingModel->SetLindhardAngleNumberHighLimit(fLindhardAnglesMuPlus,"mu+");
0254     channelingModel->SetLindhardAngleNumberHighLimit(fLindhardAnglesMuMinus,"mu-");
0255     //channelingModel->SetLindhardAngleNumberHighLimit(your_value,"your_particle");
0256
0257     /*
0258     Set the low kinetic energy cut for the model:
0259     too low energy may reduce the simulation speed (sometimes drastically),
0260     too high energy can cut off a useful physics for radiation losses.
0261     A recommended value depends a lot on the case. Usually it should be above 100 MeV,
0262     since below quantum channeling effects may become important, though lower energies
0263     are not forbidden. For energies considerably above 1 GeV and a crystal thick enough
0264     for multiphoton radiation emission,a lower energy limit of 300-500 MeV is recommended.
0265     */
0266     channelingModel->SetDefaultLowKineticEnergyLimit(fParticleMinKinEnergy);
0267     //you may set a particular limit for a certain particle type
0268     //(has a priority vs default):
0269     channelingModel->SetLowKineticEnergyLimit(fProtonMinKinEnergy,"proton");
0270     channelingModel->SetLowKineticEnergyLimit(fAntiprotonMinKinEnergy,"anti_proton");
0271     channelingModel->SetLowKineticEnergyLimit(fPiPlusMinKinEnergy,"pi+");
0272     channelingModel->SetLowKineticEnergyLimit(fPiMinusMinKinEnergy,"pi-");
0273     channelingModel->SetLowKineticEnergyLimit(fPositronMinKinEnergy,"e+");
0274     channelingModel->SetLowKineticEnergyLimit(fElectronMinKinEnergy,"e-");
0275     channelingModel->SetLowKineticEnergyLimit(fMuPlusMinKinEnergy,"mu+");
0276     channelingModel->SetLowKineticEnergyLimit(fMuMinusMinKinEnergy,"mu-");
0277     //channelingModel->SetLowKineticEnergyLimit(your_value,"your_particle");
0278
0279     /*
0280     activate the radiation model (do it only when you want to take into account the
0281     radiation production in an oriented crystal; it reduces simulation speed.)
0282     */
0283     if (fActivateRadiationModel)
0284     {
0285         channelingModel->RadiationModelActivate();
0286         G4cout << "Radiation model activated" << G4endl;
0287
0288         /*
0289         Set the number of the photons used in Monte Carlo integral in Baier-Katkov:
0290         too low number reduces the accuracy of the model;
0291         too high number reduces the calculation speed.
0292         In most of the cases 150 is a minimal secure number.
0293         */
0294         channelingModel->GetRadiationModel()->
0295                 SetSamplingPhotonsNumber(fSamplingPhotonsNumber);
0296
0297         /*
0298         Increase the statistics of sampling photons in a certain energy range.
0299         By default it is not active! In many cases it is not necessary.
0300         It is very useful for a soft spectrum part, when it is considerably below
0301         the charged particle energy, in order to increase the accuracy. It is possible
0302         to apply as many ranges as you want.
0303         NOTE: usually important for crystalline undulator
0304         CAUTION: insert only an integer number as a multiple of a photon statistics
0305         and ONLY > 1 .
0306         CAUTION: this energy range must not be beyond the range of radiation energies
0307         (i.e. below minimum photon energy (see below)) and must not intersect another
0308         range with an increased statistics if any.
0309         CAUTION: this is a multiple of the statistics of sampling photons randomly get in
0310         this energy range => make sure the total number of sampling photons is high enough
0311         to regularly get in this energy range, otherwise the statistics will not increase.
0312         */
0313         if(fTimesPhotonStatistics>1)
0314         {channelingModel->GetRadiationModel()->
0315                     AddStatisticsInPhotonEnergyRegion(fMinPhotonEnergyAddStat,
0316                                                       fMaxPhotonEnergyAddStat,
0317                                                       fTimesPhotonStatistics);}
0318
0319         /*
0320         Adjust the angular distribution of the sampling photons by
0321         changing the multiple of the opening radiation angle 1/gamma:
0322         the model should work correctly in the range 2-5;
0323         too small multiple reduces the accuracy;
0324         too high value requires more sampling photons.
0325         */
0326         channelingModel->GetRadiationModel()->
0327                 SetRadiationAngleFactor(fRadiationAngleFactor);
0328
0329         /*
0330         Set the minimal energy of radiated photon: generally it depends on which part
0331         of the spectrum you are interested in:
0332         too low number vs the charged particle energy may require more sampling photons
0333         in a total or in a particular energy range;
0334         too high number may reduce the accuracy in radiation energy loss.
0335         Generally 1 MeV is a recommended value.
0336         */
0337         channelingModel->GetRadiationModel()->SetMinPhotonEnergy(fMinPhotonEnergy);
0338
0339         /*
0340         Set the number of trajectory steps after which the radiation probability
0341         check (whether the probability is below or above of the threshold) is performed;
0342         at the first iteration also the sampling photons are generated by using
0343         the angles of this first part of the trajectory as an argument:
0344         too higher number of steps reduces the accuracy due to considerable excess
0345         of the single radiation probability threshold;
0346         too low number may reduce the accuracy of the angular distribution of
0347         sampling photons.
0348         Generally the range between 1000-10000 steps is recommended.
0349         */
0350         channelingModel->GetRadiationModel()->
0351                 SetNSmallTrajectorySteps(fNSmallTrajectorySteps);
0352     }
0353   }
0354 }
0355
0356 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0357