| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 /// \file DetectorConstruction.cc
27 /// \brief Implementation of the DetectorConstruction class
28 //
29 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
30
31 #include "DetectorConstruction.hh"
32
33 #include "G4RunManager.hh"
34 #include "G4NistManager.hh"
35 #include "G4Box.hh"
36 #include "G4LogicalVolume.hh"
37 #include "G4RegionStore.hh"
38 #include "G4VisAttributes.hh"
39 #include "PrimaryGeneratorAction.hh"
40
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
42
43 DetectorConstruction::DetectorConstruction()
44 {
45 //instantiate the messenger
46 fMessenger = new DetectorConstructionMessenger(this);
47
48 //Crystal size
49 fCrystalSize.setX(20*CLHEP::mm);
50 fCrystalSize.setY(20*CLHEP::mm);
51 fCrystalSize.setZ(0.0305*CLHEP::mm);
52
53 //Crystal planes or axes considered
54 fLattice = "(111)";
55
56 //Detector size
57 fDetectorSize.setX(10*CLHEP::cm);
58 fDetectorSize.setY(10*CLHEP::cm);
59 fDetectorSize.setZ(0.1*CLHEP::mm);
60 }
61
62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
63
64 G4VPhysicalVolume* DetectorConstruction::Construct()
65 {
66 //Check overlap option
67 G4bool checkOverlaps = true;
68
69 //Materials
70 G4NistManager* nist = G4NistManager::Instance();
71 G4Material* world_mat = nist->FindOrBuildMaterial("G4_Galactic");
72 G4Material* silicon = nist->FindOrBuildMaterial("G4_Si");
73
74 //World
75 G4Box* solidWorld = new G4Box("World", 0.2*CLHEP::m, 0.2*CLHEP::m, 10.*CLHEP::m);
76 G4LogicalVolume* logicWorld = new G4LogicalVolume(solidWorld, world_mat, "World");
77 G4VPhysicalVolume* physWorld = new G4PVPlacement
78 (0, // no rotation
79 G4ThreeVector(), // centre position
80 logicWorld, // its logical volume
81 "World", // its name
82 0, // its mother volume
83 false, // no boolean operation
84 0, // copy number
85 checkOverlaps); // overlaps checking
86 logicWorld->SetVisAttributes(G4VisAttributes::GetInvisible());
87
88
89 // --------------- Crystal ------------------------------------
90
91 //Select crystal material
92 fCrystalMaterial = nist->FindOrBuildMaterial(fCrystalMaterialStr);
93
94 //Crystal rotation angle (also the angle of crystal planes vs the beam)
95 G4RotationMatrix* crystalRotationMatrix = new G4RotationMatrix;
96 crystalRotationMatrix->rotateY(-fAngleX);
97 crystalRotationMatrix->rotateX(-fAngleY);
98
99 //Setting crystal position
100 G4ThreeVector posCrystal = G4ThreeVector(0., 0., fCrystalSize.z()/2.);
101
102 //crystal volume
103 G4Box* solidCrystal = new G4Box("Crystal",
104 fCrystalSize.x()/2,
105 fCrystalSize.y()/2,
106 fCrystalSize.z()/2.);
107
108 fLogicCrystal = new G4LogicalVolume(solidCrystal,
109 fCrystalMaterial,
110 "Crystal");
111 new G4PVPlacement(crystalRotationMatrix,
112 posCrystal,
113 fLogicCrystal,
114 "Crystal",
115 logicWorld,
116 false,
117 0,
118 checkOverlaps);
119
120 if (fActivateChannelingModel)
121 {
122 //crystal region (necessary for the FastSim model)
123 G4Region* regionCh = new G4Region("Crystal");
124 regionCh->AddRootLogicalVolume(fLogicCrystal);
125 }
126
127 //visualization attributes
128 G4VisAttributes* crystalVisAttribute =
129 new G4VisAttributes(G4Colour(1., 0., 0.));
130 crystalVisAttribute->SetForceSolid(true);
131 fLogicCrystal->SetVisAttributes(crystalVisAttribute);
132
133 //print Crystal info
134 G4cout << "Crystal material: " << fCrystalMaterial->GetName() << G4endl;
135 G4cout << "Crystal size: " << fCrystalSize.x()/CLHEP::mm
136 << "x" << fCrystalSize.y()/CLHEP::mm
137 << "x" << fCrystalSize.z()/CLHEP::mm << " mm3" << G4endl;
138 if (fActivateChannelingModel)
139 {
140 G4cout << "G4ChannelingFastSimModel activated" << G4endl;
141 G4cout << "Crystal bending angle: " << fBendingAngle << " rad" << G4endl;
142 G4cout << "Crystal Lattice: " << fLattice << G4endl;
143 G4cout << "Crystal angleX: " << fAngleX << " rad" << G4endl;
144 G4cout << "Crystal angleY: " << fAngleY << " rad" << G4endl;
145 G4cout << "ActivateRadiationModel: " << fActivateRadiationModel << G4endl;
146
147 if (fCrystallineUndulatorAmplitude > DBL_EPSILON &&
148 fCrystallineUndulatorPeriod > DBL_EPSILON) {
149 G4cout << "Crystalline undulator activated: " << G4endl;
150 G4cout << "undulator amplitude: "
151 << fCrystallineUndulatorAmplitude/CLHEP::nm
152 << " nm" << G4endl;
153 G4cout << "undulator period: "
154 << fCrystallineUndulatorPeriod/CLHEP::mm
155 << " mm" << G4endl;
156 G4cout << "undulator phase: "
157 << fCrystallineUndulatorPhase
158 << " rad" << G4endl;
159 }
160
161 G4cout << G4endl;
162 }
163 else
164 {
165 G4cout << "G4ChannelingFastSimModel is not activated" << G4endl << G4endl;
166 }
167
168 // --------------- Detector -----------------------------------
169 //Setting detector position
170 G4ThreeVector posDetector =
171 G4ThreeVector(0, 0, fDetectorFrontPosZ+fDetectorSize.z()/2.);
172
173 //particle detector volume
174 G4Box* detector = new G4Box("Detector",
175 fDetectorSize.x()/2,
176 fDetectorSize.y()/2,
177 fDetectorSize.z()/2.);
178
179 G4LogicalVolume* logicDetector = new G4LogicalVolume(detector,
180 silicon,
181 "Detector");
182 new G4PVPlacement(0,
183 posDetector,
184 logicDetector,
185 "Detector",
186 logicWorld,
187 false,
188 0,
189 checkOverlaps);
190
191 //visualization attributes
192 G4VisAttributes* detectorVisAttribute =
193 new G4VisAttributes(G4Colour(1., 0., 1., 0.3));
194 detectorVisAttribute->SetForceSolid(true);
195 logicDetector->SetVisAttributes(detectorVisAttribute);
196
197 //always return the physical World
198 return physWorld;
199 }
200
201 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
202
203 void DetectorConstruction::ConstructSDandField()
204 {
205 if (fActivateChannelingModel)
206 {
207 // --------------- fast simulation ----------------------------
208 //extract the region of the crystal from the store
209 G4RegionStore* regionStore = G4RegionStore::GetInstance();
210 G4Region* regionCh = regionStore->GetRegion("Crystal");
211
212 //create the channeling model for this region
213 G4ChannelingFastSimModel* channelingModel =
214 new G4ChannelingFastSimModel("ChannelingModel", regionCh);
215
216 //activate the channeling model
217 channelingModel->Input(fCrystalMaterial, fLattice, fPotentialPath);
218 //setting bending angle of the crystal planes (default is 0)
219 channelingModel->GetCrystalData()->SetBendingAngle(fBendingAngle, fLogicCrystal);
220
221 //setting crystalline undulator parameters
222 //NOTE: they are incompatible with a bent crystal
223 if (fCrystallineUndulatorAmplitude > DBL_EPSILON &&
224 fCrystallineUndulatorPeriod > DBL_EPSILON)
225 {
226 channelingModel->GetCrystalData()->SetCrystallineUndulatorParameters(
227 fCrystallineUndulatorAmplitude,
228 fCrystallineUndulatorPeriod,
229 fCrystallineUndulatorPhase,
230 fLogicCrystal);
231 }
232
233 /*
234 Set the multiple of critical channeling angles (Lindhard angles) which defines
235 the angular cut of the model (otherwise standard Geant4 is active):
236 a number too low reduces the accuracy of the model;
237 a number too high sometimes drastically reduces the simulation speed.
238 The default value is 100, while for many problems 10-20 is ok.
239 CAUTION: If you set this value to 1 meaning the angular cut = 1*Lindhard angle,
240 this will cut off the physics of overbarrier motion, still important even if
241 the particle is not in channeling. This will provide incorrect results.
242 */
243 channelingModel->SetDefaultLindhardAngleNumberHighLimit(fLindhardAngles);
244 //you may set a particular limit for a certain particle type
245 //(has a priority vs default):
246 channelingModel->SetLindhardAngleNumberHighLimit(fLindhardAnglesProton,"proton");
247 channelingModel->SetLindhardAngleNumberHighLimit(fLindhardAnglesAntiproton,
248 "anti_proton");
249 channelingModel->SetLindhardAngleNumberHighLimit(fLindhardAnglesPiPlus,"pi+");
250 channelingModel->SetLindhardAngleNumberHighLimit(fLindhardAnglesPiMinus,"pi-");
251 channelingModel->SetLindhardAngleNumberHighLimit(fLindhardAnglesPositron,"e+");
252 channelingModel->SetLindhardAngleNumberHighLimit(fLindhardAnglesElectron,"e-");
253 channelingModel->SetLindhardAngleNumberHighLimit(fLindhardAnglesMuPlus,"mu+");
254 channelingModel->SetLindhardAngleNumberHighLimit(fLindhardAnglesMuMinus,"mu-");
255 //channelingModel->SetLindhardAngleNumberHighLimit(your_value,"your_particle");
256
257 /*
258 Set the low kinetic energy cut for the model:
259 too low energy may reduce the simulation speed (sometimes drastically),
260 too high energy can cut off a useful physics for radiation losses.
261 A recommended value depends a lot on the case. Usually it should be above 100 MeV,
262 since below quantum channeling effects may become important, though lower energies
263 are not forbidden. For energies considerably above 1 GeV and a crystal thick enough
264 for multiphoton radiation emission,a lower energy limit of 300-500 MeV is recommended.
265 */
266 channelingModel->SetDefaultLowKineticEnergyLimit(fParticleMinKinEnergy);
267 //you may set a particular limit for a certain particle type
268 //(has a priority vs default):
269 channelingModel->SetLowKineticEnergyLimit(fProtonMinKinEnergy,"proton");
270 channelingModel->SetLowKineticEnergyLimit(fAntiprotonMinKinEnergy,"anti_proton");
271 channelingModel->SetLowKineticEnergyLimit(fPiPlusMinKinEnergy,"pi+");
272 channelingModel->SetLowKineticEnergyLimit(fPiMinusMinKinEnergy,"pi-");
273 channelingModel->SetLowKineticEnergyLimit(fPositronMinKinEnergy,"e+");
274 channelingModel->SetLowKineticEnergyLimit(fElectronMinKinEnergy,"e-");
275 channelingModel->SetLowKineticEnergyLimit(fMuPlusMinKinEnergy,"mu+");
276 channelingModel->SetLowKineticEnergyLimit(fMuMinusMinKinEnergy,"mu-");
277 //channelingModel->SetLowKineticEnergyLimit(your_value,"your_particle");
278
279 /*
280 activate the radiation model (do it only when you want to take into account the
281 radiation production in an oriented crystal; it reduces simulation speed.)
282 */
283 if (fActivateRadiationModel)
284 {
285 channelingModel->RadiationModelActivate();
286 G4cout << "Radiation model activated" << G4endl;
287
288 /*
289 Set the number of the photons used in Monte Carlo integral in Baier-Katkov:
290 too low number reduces the accuracy of the model;
291 too high number reduces the calculation speed.
292 In most of the cases 150 is a minimal secure number.
293 */
294 channelingModel->GetRadiationModel()->
295 SetSamplingPhotonsNumber(fSamplingPhotonsNumber);
296
297 /*
298 Increase the statistics of sampling photons in a certain energy range.
299 By default it is not active! In many cases it is not necessary.
300 It is very useful for a soft spectrum part, when it is considerably below
301 the charged particle energy, in order to increase the accuracy. It is possible
302 to apply as many ranges as you want.
303 NOTE: usually important for crystalline undulator
304 CAUTION: insert only an integer number as a multiple of a photon statistics
305 and ONLY > 1 .
306 CAUTION: this energy range must not be beyond the range of radiation energies
307 (i.e. below minimum photon energy (see below)) and must not intersect another
308 range with an increased statistics if any.
309 CAUTION: this is a multiple of the statistics of sampling photons randomly get in
310 this energy range => make sure the total number of sampling photons is high enough
311 to regularly get in this energy range, otherwise the statistics will not increase.
312 */
313 if(fTimesPhotonStatistics>1)
314 {channelingModel->GetRadiationModel()->
315 AddStatisticsInPhotonEnergyRegion(fMinPhotonEnergyAddStat,
316 fMaxPhotonEnergyAddStat,
317 fTimesPhotonStatistics);}
318
319 /*
320 Adjust the angular distribution of the sampling photons by
321 changing the multiple of the opening radiation angle 1/gamma:
322 the model should work correctly in the range 2-5;
323 too small multiple reduces the accuracy;
324 too high value requires more sampling photons.
325 */
326 channelingModel->GetRadiationModel()->
327 SetRadiationAngleFactor(fRadiationAngleFactor);
328
329 /*
330 Set the minimal energy of radiated photon: generally it depends on which part
331 of the spectrum you are interested in:
332 too low number vs the charged particle energy may require more sampling photons
333 in a total or in a particular energy range;
334 too high number may reduce the accuracy in radiation energy loss.
335 Generally 1 MeV is a recommended value.
336 */
337 channelingModel->GetRadiationModel()->SetMinPhotonEnergy(fMinPhotonEnergy);
338
339 /*
340 Set the number of trajectory steps after which the radiation probability
341 check (whether the probability is below or above of the threshold) is performed;
342 at the first iteration also the sampling photons are generated by using
343 the angles of this first part of the trajectory as an argument:
344 too higher number of steps reduces the accuracy due to considerable excess
345 of the single radiation probability threshold;
346 too low number may reduce the accuracy of the angular distribution of
347 sampling photons.
348 Generally the range between 1000-10000 steps is recommended.
349 */
350 channelingModel->GetRadiationModel()->
351 SetNSmallTrajectorySteps(fNSmallTrajectorySteps);
352 }
353 }
354 }
355
356 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
357
358