| 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 // Author: Alexei Sytov
27 // Co-author: Gianfranco Paternò (modifications & testing)
28 // On the base of the CRYSTALRAD realization of channeling model:
29 // A. I. Sytov, V. V. Tikhomirov, and L. Bandiera PRAB 22, 064601 (2019)
30
31 /// \file G4ChannelingFastSimModel.cc
32 /// \brief Implementation of the G4ChannelingFastSimModel class
33 //
34 //
35 //
36 #include "G4ChannelingFastSimModel.hh"
37
38 #include "Randomize.hh"
39
40 #include "G4TransportationManager.hh"
41 #include "G4SystemOfUnits.hh"
42
43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
44
45 G4ChannelingFastSimModel::G4ChannelingFastSimModel(const G4String& modelName,
46 G4Region* envelope)
47 : G4VFastSimulationModel(modelName, envelope)
48 {
49 }
50
51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
52
53 G4ChannelingFastSimModel::G4ChannelingFastSimModel(const G4String& modelName)
54 : G4VFastSimulationModel(modelName)
55 {
56
57 }
58
59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
60
61 G4ChannelingFastSimModel::~G4ChannelingFastSimModel()
62 {
63 }
64
65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
66
67 G4bool G4ChannelingFastSimModel::IsApplicable(const G4ParticleDefinition& particleType)
68 {
69 return std::abs(particleType.GetPDGCharge())>DBL_EPSILON;
70 }
71
72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
73
74 G4bool G4ChannelingFastSimModel::ModelTrigger(const G4FastTrack& fastTrack)
75 {
76 //default output
77 G4bool modelTrigger = false;
78
79 G4int particleDefinitionID =
80 fastTrack.GetPrimaryTrack()->GetParticleDefinition()->GetParticleDefinitionID();
81 //kinetic energy
82 G4double ekinetic = fastTrack.GetPrimaryTrack()->GetKineticEnergy();
83
84 //energy cut, at the beginning, to not check everything else
85 if(ekinetic > GetLowKineticEnergyLimit(particleDefinitionID))
86 {
87 //current logical volume
88 G4LogicalVolume* crystallogic = fastTrack.GetEnvelopeLogicalVolume();
89 fCrystalData->SetGeometryParameters(crystallogic);
90
91 G4ThreeVector momentumDirection = fastTrack.GetPrimaryTrackLocalDirection();
92 // the particle angle vs crystal plane or axis
93 G4double angle = std::atan(momentumDirection.x()/momentumDirection.z());
94 //recalculate angle into the lattice reference system
95 angle = fCrystalData->
96 AngleXFromBoxToLattice(angle,
97 (fCrystalData->CoordinatesFromBoxToLattice(
98 fastTrack.GetPrimaryTrackLocalPosition())).z());
99 if (fCrystalData->GetModel()==2)
100 {
101 angle = std::sqrt(angle*angle+
102 std::pow(std::atan(momentumDirection.y()/
103 momentumDirection.z()),2));
104 }
105
106 //particle mass
107 G4double mass = fastTrack.GetPrimaryTrack()->GetParticleDefinition()->GetPDGMass();
108 //particle total energy
109 G4double etotal = mass + ekinetic;
110 //particle charge
111 G4double charge = fastTrack.GetPrimaryTrack()->
112 GetParticleDefinition()->GetPDGCharge();
113
114 //Particle position
115 G4ThreeVector xyz0 = fastTrack.GetPrimaryTrackLocalPosition();
116 //Step estimate
117 G4double dz0 = fCrystalData->GetMaxSimulationStep(etotal,mass,charge);
118 xyz0 += 2*dz0*momentumDirection;//overestimated particle shift on the next step
119 //in channeling
120
121 //Applies the parameterisation not at the last step, only forward local direction
122 //above low energy limit and below angular limit
123
124 modelTrigger = (crystallogic->GetSolid()->
125 Inside(xyz0)==kInside) &&
126 momentumDirection.z()>0. &&
127 std::abs(angle) <
128 std::max(
129 GetLindhardAngleNumberHighLimit(particleDefinitionID) *
130 fCrystalData->GetLindhardAngle(etotal,
131 mass,
132 charge),
133 GetHighAngleLimit(particleDefinitionID));
134 }
135
136 return modelTrigger;
137 }
138
139 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
140
141 void G4ChannelingFastSimModel::DoIt(const G4FastTrack& fastTrack,
142 G4FastStep& fastStep)
143 {
144 G4double etotal;//particle total energy
145 G4double etotalPreStep;//etotal at the previous step
146 G4double etotalToSetParticleProperties;//etotal value at which
147 //SetParticleProperties is calculated
148 G4double ekinetic = 0;//kinetic energy
149 G4double eDeposited = 0.;//deposited energy along the trajectory
150 G4double elossAccum = 0;// accumulate local energy loss (not radiation)
151 G4double mass; //particle mass
152 G4double charge;//particle charge
153 G4double tGlobal; //global time
154 G4double tGlobalPreStep; //global time at the previous step
155 G4ThreeVector xyz0;// the coordinates in the local reference system of the volume
156 G4ThreeVector xyz0PreStep;// xyz at the previous step
157 G4ThreeVector xyz;// the coordinates in the co-rotating reference system within
158 //a channel (elementary periodic cell)
159 G4double x,y,z; // the coordinates in the co-rotating reference system within
160 //a channel (elementary periodic cell)
161 G4double tx0,ty0; // the angles in the local reference system of the volume
162 G4double tx,ty; // the angles in the co-rotating reference system within
163 //a channel (elementary periodic cell)
164 G4double txPreStep,tyPreStep;// tx,ty at the previous step
165 G4ThreeVector momentumDirection;
166 G4ThreeVector scatteringAnglesAndEnergyLoss;//output of scattering functions
167 G4double lindhardAngleNumberHighLimit0; //current high limit of the angle expressed in
168 //[Lindhard angle] units
169 G4double highAngleLimit0; //current absolute high limit of the angle expressed
170
171 //coordinates in Runge-Kutta calculations
172 G4double x1=0.,x2=0.,x3=0.,x4=0.,y1=0.,y2=0.,y3=0.,y4=0.;
173 //angles in Runge-Kutta calculations
174 G4double tx1=0.,tx2=0.,tx3=0.,tx4=0.,ty1=0.,ty2=0.,ty3=0.,ty4=0.;
175 //variables in Runge-Kutta calculations
176 G4double kvx1=0.,kvx2=0.,kvx3=0.,kvx4=0.,kvy1=0.,kvy2=0.,kvy3=0.,kvy4=0.;
177 //simulation step along z (internal step of the model) and its parts
178 G4double dz,dzd3,dzd8;//dzd3 = dz/3; dzd8 = dz/8;
179 //simulation step along the momentum direction
180 G4double momentumDirectionStep;
181 //effective simulation step (taking into account nuclear density along the trajectory)
182 G4double effectiveStep;
183
184 // flag, if Inside(xyz0) switches to kInside
185 G4bool inside = false;
186
187 G4LogicalVolume* crystallogic = fastTrack.GetEnvelopeLogicalVolume();
188 fCrystalData->SetGeometryParameters(crystallogic);
189
190 //set the max number of secondaries (photons) that can be added at this fastStep
191 if (fRad)
192 {
193 fastStep.SetNumberOfSecondaryTracks(fMaxPhotonsProducedPerStep);
194 //reseting the BaierKatkov integral to start it with the new trajectory
195 fBaierKatkov->ResetRadIntegral();//to avoid any memory from the previous trajectory
196 }
197
198 mass = fastTrack.GetPrimaryTrack()->GetParticleDefinition()->GetPDGMass();
199 etotal = mass + fastTrack.GetPrimaryTrack()->GetKineticEnergy();
200 charge = fastTrack.GetPrimaryTrack()->GetParticleDefinition()->GetPDGCharge();
201
202 G4String particleName =
203 fastTrack.GetPrimaryTrack()->GetParticleDefinition()->GetParticleName();
204
205 lindhardAngleNumberHighLimit0 =
206 GetLindhardAngleNumberHighLimit(fastTrack.GetPrimaryTrack()->
207 GetParticleDefinition()->GetParticleDefinitionID());
208 highAngleLimit0 = GetHighAngleLimit(fastTrack.GetPrimaryTrack()->
209 GetParticleDefinition()->GetParticleDefinitionID());
210
211 //set fCrystalData parameters depending on the particle parameters
212 fCrystalData->SetParticleProperties(etotal, mass, charge, particleName);
213
214 //global time
215 tGlobal = fastTrack.GetPrimaryTrack()->GetGlobalTime();
216
217 //coordinates in the co-rotating reference system within a channel
218 xyz0= fastTrack.GetPrimaryTrackLocalPosition();
219 xyz = fCrystalData->CoordinatesFromBoxToLattice(xyz0);
220 x=xyz.x();
221 y=xyz.y();
222 z=xyz.z();
223
224 momentumDirection=fastTrack.GetPrimaryTrackLocalDirection();
225 //angle in the co-rotating reference system within a channel
226 //(!!! ONLY FORWARD DIRECTION, momentumDirection.getZ()>0,
227 //valid for high energies defined by the standard energy cuts)
228 tx0 = std::atan(momentumDirection.x()/momentumDirection.z());
229 ty0 = std::atan(momentumDirection.y()/momentumDirection.z());
230
231 //angles in the co-rotating reference system within a channel
232 tx = fCrystalData->AngleXFromBoxToLattice(tx0,z);
233 ty = ty0;
234
235 etotalToSetParticleProperties = etotal*0.999;
236 G4bool inCrystal=true;//flag necessary to escape the cycle (at inCrystal=0;)
237 //do calculations until the particle is inside the volume
238 do
239 {
240 //remember the global time before the next step dz
241 tGlobalPreStep=tGlobal;
242 //remember the coordinates before the next step dz
243 xyz0PreStep = xyz0;
244 //remember the angles and the total energy before the step dz
245 txPreStep = tx;
246 tyPreStep = ty;
247 etotalPreStep = etotal;
248
249 dz = fCrystalData->GetSimulationStep(tx,ty);
250 dzd3=dz/3;
251 dzd8=dz/8;
252
253 //trajectory calculation:
254 //Runge-Cutt "3/8"
255 //fCrystalData->GetCurv(z)*fCrystalData->GetCorrectionZ() is due to dependence of
256 //the radius on x; GetCurv gets 1/R for the central ("central plane/axis")
257
258 //first step
259 kvx1=fCrystalData->Ex(x,y);
260 x1=x+tx*dzd3;
261 tx1=tx+(kvx1-fCrystalData->GetCurv(z)*fCrystalData->GetCorrectionZ())*dzd3;
262 if (fCrystalData->GetModel()==2)
263 {
264 kvy1=fCrystalData->Ey(x,y);
265 y1=y+ty*dzd3;
266 ty1=ty+kvy1*dzd3;
267 }
268
269 //second step
270 kvx2=fCrystalData->Ex(x1,y1);
271 x2=x-tx*dzd3+tx1*dz;
272 tx2=tx-(kvx1-fCrystalData->GetCurv(z)*fCrystalData->GetCorrectionZ())*dzd3+
273 (kvx2-fCrystalData->GetCurv(z)*fCrystalData->GetCorrectionZ())*dz;
274 if (fCrystalData->GetModel()==2)
275 {
276 kvy2=fCrystalData->Ey(x1,y1);
277 y2=y-ty*dzd3+ty1*dz;
278 ty2=ty-kvy1*dzd3+kvy2*dz;
279 }
280
281 //third step
282 kvx3=fCrystalData->Ex(x2,y2);
283 x3=x+(tx-tx1+tx2)*dz;
284 tx3=tx+(kvx1-kvx2+kvx3-fCrystalData->GetCurv(z)*fCrystalData->GetCorrectionZ())*dz;
285 if (fCrystalData->GetModel()==2)
286 {
287 kvy3=fCrystalData->Ey(x2,y2);
288 y3=y+(ty-ty1+ty2)*dz;
289 ty3=ty+(kvy1-kvy2+kvy3)*dz;
290 }
291
292 //fourth step
293 kvx4=fCrystalData->Ex(x3,y3);
294 x4=x+(tx+3.*tx1+3.*tx2+tx3)*dzd8;
295 tx4=tx+(kvx1+3.*kvx2+3.*kvx3+kvx4)*dzd8-
296 fCrystalData->GetCurv(z)*fCrystalData->GetCorrectionZ()*dz;
297 if (fCrystalData->GetModel()==2)
298 {
299 kvy4=fCrystalData->Ey(x3,y3);
300 y4=y+(ty+3.*ty1+3.*ty2+ty3)*dzd8;
301 ty4=ty+(kvy1+3.*kvy2+3.*kvy3+kvy4)*dzd8;
302 }
303 else
304 {
305 y4 =y+ty*dz;
306 ty4=ty;
307 }
308
309 x=x4;
310 tx=tx4;
311 y=y4;
312 ty=ty4;
313
314 z+=dz*fCrystalData->GetCorrectionZ();//motion along the z coordinate
315 //("central plane/axis", no current plane/axis)
316
317 xyz = fCrystalData->ChannelChange(x,y,z);
318 x=xyz.x();
319 y=xyz.y();
320 z=xyz.z();
321
322 //the coordinates in the local reference system of the volume
323 //this vector will be used in the cycle escape condition and
324 //in the radiation model (if activated)
325 xyz0=fCrystalData->CoordinatesFromLatticeToBox(xyz);
326
327 momentumDirectionStep=
328 dz*std::sqrt(1+std::pow(std::tan(tx),2)+std::pow(std::tan(ty),2));
329 tGlobal+=momentumDirectionStep/(fCrystalData->GetBeta())/CLHEP::c_light;
330
331 //default scattering and energy loss 0
332 scatteringAnglesAndEnergyLoss = G4ThreeVector(0.,0.,0.);
333
334 //calculate separately for each element of the crystal
335 for (G4int i = 0; i < fCrystalData->GetNelements(); i++)
336 {
337 //effective step taking into account nuclear density along the trajectory
338 effectiveStep = momentumDirectionStep*fCrystalData->NuclearDensity(x,y,i);
339 //Coulomb scattering on screened atomic potential (both multiple and single)
340 scatteringAnglesAndEnergyLoss += fCrystalData->
341 CoulombAtomicScattering(effectiveStep,momentumDirectionStep,i);
342
343 //Amorphous part of ionization energy losses
344 elossAccum += fCrystalData->IonizationLosses(momentumDirectionStep, i);
345 }
346 //electron scattering and coherent part of ionization energy losses
347 scatteringAnglesAndEnergyLoss += fCrystalData->CoulombElectronScattering(
348 fCrystalData->MinIonizationEnergy(x,y),
349 fCrystalData->ElectronDensity(x,y),
350 momentumDirectionStep);
351 tx += scatteringAnglesAndEnergyLoss.x();
352 ty += scatteringAnglesAndEnergyLoss.y();
353 elossAccum += scatteringAnglesAndEnergyLoss.z();
354
355 // recalculate the energy depended parameters
356 //(only if the energy decreased enough, not at each step)
357 if (etotalToSetParticleProperties>etotal)
358 {
359 fCrystalData->SetParticleProperties(etotal, mass, charge, particleName);
360 etotalToSetParticleProperties = etotal*0.999;
361 }
362
363 //chain of conditions to escape the cycle
364 // if Inside(xyz0)==kInside has been already true
365 //(a particle has been inside the crystal)
366 if (inside)
367 {
368 // if low energy
369 if (etotal-mass<=GetLowKineticEnergyLimit(fastTrack.GetPrimaryTrack()->
370 GetParticleDefinition()->
371 GetParticleDefinitionID()))
372 {inCrystal = false;}//escape the cycle
373 //check if the angle w.r.t. the axes or planes is too high =>
374 //return to standard Geant4:
375 else if (fCrystalData->GetModel()==1) //1D model, field of planes
376 {
377 //if the angle w.r.t. the planes is too high
378 if (std::abs(tx) >=
379 std::max(lindhardAngleNumberHighLimit0*
380 fCrystalData->GetLindhardAngle(),
381 highAngleLimit0))
382 {inCrystal = false;}//escape the cycle
383 }
384 else if (fCrystalData->GetModel()==2) //2D model, field of axes
385 {
386 //if the angle w.r.t. the axes is too high
387 if (std::sqrt(tx*tx+ty*ty) >=
388 std::max(lindhardAngleNumberHighLimit0*
389 fCrystalData->GetLindhardAngle(),
390 highAngleLimit0))
391 {inCrystal = false;}//escape the cycle
392 }
393
394 //radiation production & radiation energy losses
395 //works only if the radiation model is activated
396 if (fRad)
397 {
398 //back to the local reference system of the volume
399 tx0 = fCrystalData->AngleXFromLatticeToBox(tx,z);
400 ty0 = ty;
401 //xyz0 was calculated above
402
403 //running the radiation model and checking if a photon has been emitted
404 if(fBaierKatkov->DoRadiation(etotal,mass,
405 tx0,ty0,
406 scatteringAnglesAndEnergyLoss.x(),
407 scatteringAnglesAndEnergyLoss.y(),
408 momentumDirectionStep,tGlobal,xyz0,
409 crystallogic->
410 GetSolid()->
411 Inside(xyz0)!=kInside&&inCrystal))
412 // also it was checked if the particle is escaping the volume
413 // calculate the radiation integral immidiately in this case
414 {
415 //a photon has been emitted!
416 //shift the particle back into the radiation point
417 etotal = fBaierKatkov->GetParticleNewTotalEnergy();
418 tx0 = fBaierKatkov->GetParticleNewAngleX();
419 ty0 = fBaierKatkov->GetParticleNewAngleY();
420 tGlobal = fBaierKatkov->GetNewGlobalTime();
421 xyz0 = fBaierKatkov->GetParticleNewCoordinateXYZ();
422
423 //add secondary photon
424 fBaierKatkov->GeneratePhoton(fastStep);
425
426 //particle energy was changed
427 fCrystalData->SetParticleProperties(etotal, mass, charge, particleName);
428
429 //coordinates in the co-rotating reference system within a channel
430 xyz = fCrystalData->CoordinatesFromBoxToLattice(xyz0);
431 x=xyz.x();
432 y=xyz.y();
433 z=xyz.z();
434
435 //angles in the co-rotating reference system within a channel
436 tx = fCrystalData->AngleXFromBoxToLattice(tx0,z);
437 ty = ty0;
438 }
439 }
440 else
441 {
442 //we calculate deposited energy and energy losses ONLY in absence
443 //of radiation otherwise we do it only at the end of model
444 etotal -= elossAccum;
445 eDeposited += elossAccum;
446 elossAccum=0;
447 ekinetic = etotal-mass;
448 if(ekinetic<1*keV)
449 {
450 G4cout << "Warning in G4ChannelingFastSimModel: " <<
451 ekinetic << "<" << 1*keV << " !" << G4endl;
452 eDeposited-=(1*keV-ekinetic);
453 ekinetic = 1*keV;
454 G4cout << "Setting deposited energy=" <<
455 eDeposited << " & ekinetic=" << ekinetic << G4endl;
456 etotal = mass+ekinetic;
457 }
458 }
459
460 //precise check if the particle is escaping the volume
461 if (crystallogic->GetSolid()->
462 Inside(xyz0)!=kInside)
463 {
464 //one step back to remain inside the volume
465 //after the escape of the volume
466 tGlobal = tGlobalPreStep;
467 xyz0 = xyz0PreStep;
468 tx = txPreStep;
469 ty = tyPreStep;
470 etotal = etotalPreStep;
471 z-=dz*fCrystalData->GetCorrectionZ();
472 // change the flag => this particle will not enter
473 // the model before escape this volume
474
475 inCrystal = false; //escape the cycle
476 }
477 }
478 else
479 {
480 // if Inside(xyz0)==kInside we can enable checking of particle escape
481 if (crystallogic->GetSolid()->
482 Inside(xyz0)==kInside)
483 {inside = true;}
484 // a very rare case, if a particle remains
485 // on the boundary and escapes the crystal
486 else if (crystallogic->GetSolid()->
487 Inside(xyz0)==kOutside)
488 {inCrystal = false;}//escape the cycle
489 }
490 }
491 while (inCrystal);
492
493 //the angles in the local reference system of the volume
494 tx0 = fCrystalData->AngleXFromLatticeToBox(tx,z);
495 ty0 = ty;
496
497 //set global time
498 fastStep.ProposePrimaryTrackFinalTime(tGlobal);
499 //set final position
500 fastStep.ProposePrimaryTrackFinalPosition(xyz0);
501
502 //set deposited energy (due to ionization)
503 etotal -= elossAccum;
504 eDeposited += elossAccum;
505 ekinetic = etotal-mass;
506 if(ekinetic<1*keV)
507 {
508 G4cout << "Warning in G4ChannelingFastSimModel: " <<
509 ekinetic << "<" << 1*keV << " !" << G4endl;
510 eDeposited-=(1*keV-ekinetic);
511 ekinetic = 1*keV;
512 G4cout << "Setting deposited energy=" <<
513 eDeposited << " & ekinetic=" << ekinetic << G4endl;
514 }
515 fastStep.ProposeTotalEnergyDeposited(eDeposited);
516 //set final kinetic energy
517 fastStep.ProposePrimaryTrackFinalKineticEnergy(ekinetic);
518
519
520 //set final momentum direction
521 G4double momentumDirectionZ =
522 1./std::sqrt(1.+std::pow(std::tan(tx0),2)+std::pow(std::tan(ty0),2));
523 fastStep.ProposePrimaryTrackFinalMomentumDirection(
524 G4ThreeVector(momentumDirectionZ*std::tan(tx0),
525 momentumDirectionZ*std::tan(ty0),
526 momentumDirectionZ));
527 }
528
529 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
530
531 void G4ChannelingFastSimModel::Input(const G4Material *crystal,
532 const G4String &lattice,
533 const G4String &filePath)
534 {
535 //initializing the class with containing all
536 //the crystal material and crystal lattice data and
537 //Channeling scattering and ionization processes
538 fCrystalData = new G4ChannelingFastSimCrystalData();
539 //setting all the crystal material and lattice data
540 fCrystalData->SetMaterialProperties(crystal,lattice,filePath);
541
542 //setting default low energy cuts for kinetic energy
543 SetLowKineticEnergyLimit(1*GeV,"proton");
544 SetLowKineticEnergyLimit(1*GeV,"anti_proton");
545 SetLowKineticEnergyLimit(200*MeV,"e-");
546 SetLowKineticEnergyLimit(200*MeV,"e+");
547
548 //set the model high limit of the angle expressed in [Lindhard angle] units
549 SetLindhardAngleNumberHighLimit(100.,"proton");
550 SetLindhardAngleNumberHighLimit(100.,"anti_proton");
551 SetLindhardAngleNumberHighLimit(100.,"e-");
552 SetLindhardAngleNumberHighLimit(100.,"e+");
553 }
554
555 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
556
557 void G4ChannelingFastSimModel::RadiationModelActivate()
558 {
559 fRad = true;
560 //activate the Baier-Katkov radiation model
561 fBaierKatkov = new G4BaierKatkov();
562 }
563
564 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
565