Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4DNARuddIonisationExtendedModel.cc
Go to the documentation of this file.
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// $Id$
27// GEANT4 tag $Name: $
28//
29
30
31// Modified by Z. Francis to handle HZE && inverse rudd function sampling 26-10-2010
32
35#include "G4SystemOfUnits.hh"
37#include "G4LossTableManager.hh"
40
41//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
42
43using namespace std;
44
45//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
46
48 const G4String& nam)
49 :G4VEmModel(nam),isInitialised(false)
50{
51 // nistwater = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
52 fpWaterDensity = 0;
53
54 slaterEffectiveCharge[0]=0.;
55 slaterEffectiveCharge[1]=0.;
56 slaterEffectiveCharge[2]=0.;
57 sCoefficient[0]=0.;
58 sCoefficient[1]=0.;
59 sCoefficient[2]=0.;
60
61 lowEnergyLimitForA[1] = 0 * eV;
62 lowEnergyLimitForA[2] = 0 * eV;
63 lowEnergyLimitForA[3] = 0 * eV;
64 lowEnergyLimitOfModelForA[1] = 100 * eV;
65 lowEnergyLimitOfModelForA[4] = 1 * keV;
66 lowEnergyLimitOfModelForA[5] = 0.5 * MeV; // For A = 3 or above, limit is MeV/uma
67 killBelowEnergyForA[1] = lowEnergyLimitOfModelForA[1];
68 killBelowEnergyForA[4] = lowEnergyLimitOfModelForA[4];
69 killBelowEnergyForA[5] = lowEnergyLimitOfModelForA[5];
70
71
72 verboseLevel= 0;
73 // Verbosity scale:
74 // 0 = nothing
75 // 1 = warning for energy non-conservation
76 // 2 = details of energy budget
77 // 3 = calculation of cross sections, file openings, sampling of atoms
78 // 4 = entering in methods
79
80 if( verboseLevel>0 )
81 {
82 G4cout << "Rudd ionisation model is constructed " << G4endl;
83 }
84
85 //Mark this model as "applicable" for atomic deexcitation
87 fAtomDeexcitation = 0;
89}
90
91//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
92
94{
95 // Cross section
96
97 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
98 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
99 {
100 G4DNACrossSectionDataSet* table = pos->second;
101 delete table;
102 }
103
104 // The following removal is forbidden G4VEnergyLossModel takes care of deletion
105 // however coverity will signal this as an error
106 //if (fAtomDeexcitation) {delete fAtomDeexcitation;}
107
108}
109
110//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
111
113 const G4DataVector& /*cuts*/)
114{
115
116 if (verboseLevel > 3)
117 G4cout << "Calling G4DNARuddIonisationExtendedModel::Initialise()" << G4endl;
118
119 // Energy limits
120
121 G4String fileProton("dna/sigma_ionisation_p_rudd");
122 G4String fileHydrogen("dna/sigma_ionisation_h_rudd");
123 G4String fileAlphaPlusPlus("dna/sigma_ionisation_alphaplusplus_rudd");
124 G4String fileAlphaPlus("dna/sigma_ionisation_alphaplus_rudd");
125 G4String fileHelium("dna/sigma_ionisation_he_rudd");
126 G4String fileCarbon("dna/sigma_ionisation_c_rudd");
127 G4String fileNitrogen("dna/sigma_ionisation_n_rudd");
128 G4String fileOxygen("dna/sigma_ionisation_o_rudd");
129 G4String fileIron("dna/sigma_ionisation_fe_rudd");
130
131 G4DNAGenericIonsManager *instance;
134 G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen");
135 G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
136 G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
137 G4ParticleDefinition* heliumDef = instance->GetIon("helium");
138 G4ParticleDefinition* carbonDef = instance->GetIon("carbon");
139 G4ParticleDefinition* nitrogenDef = instance->GetIon("nitrogen");
140 G4ParticleDefinition* oxygenDef = instance->GetIon("oxygen");
141 G4ParticleDefinition* ironDef = instance->GetIon("iron");
142
143 G4String proton;
144 G4String hydrogen;
145 G4String alphaPlusPlus;
146 G4String alphaPlus;
147 G4String helium;
148 G4String carbon;
149 G4String nitrogen;
150 G4String oxygen;
151 G4String iron;
152
153 G4double scaleFactor = 1 * m*m;
154
155 // LIMITS AND DATA
156
157 proton = protonDef->GetParticleName();
158 tableFile[proton] = fileProton;
159 lowEnergyLimit[proton] = lowEnergyLimitForA[1];
160 highEnergyLimit[proton] = 500. * keV;
161
162 // Cross section
163
165 eV,
166 scaleFactor );
167 tableProton->LoadData(fileProton);
168 tableData[proton] = tableProton;
169
170 // **********************************************************************************************
171
172 hydrogen = hydrogenDef->GetParticleName();
173 tableFile[hydrogen] = fileHydrogen;
174
175 lowEnergyLimit[hydrogen] = lowEnergyLimitForA[1];
176 highEnergyLimit[hydrogen] = 100. * MeV;
177
178 // Cross section
179
181 eV,
182 scaleFactor );
183 tableHydrogen->LoadData(fileHydrogen);
184
185 tableData[hydrogen] = tableHydrogen;
186
187 // **********************************************************************************************
188
189 alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
190 tableFile[alphaPlusPlus] = fileAlphaPlusPlus;
191
192 lowEnergyLimit[alphaPlusPlus] = lowEnergyLimitForA[4];
193 highEnergyLimit[alphaPlusPlus] = 400. * MeV;
194
195 // Cross section
196
198 eV,
199 scaleFactor );
200 tableAlphaPlusPlus->LoadData(fileAlphaPlusPlus);
201
202 tableData[alphaPlusPlus] = tableAlphaPlusPlus;
203
204 // **********************************************************************************************
205
206 alphaPlus = alphaPlusDef->GetParticleName();
207 tableFile[alphaPlus] = fileAlphaPlus;
208
209 lowEnergyLimit[alphaPlus] = lowEnergyLimitForA[4];
210 highEnergyLimit[alphaPlus] = 400. * MeV;
211
212 // Cross section
213
215 eV,
216 scaleFactor );
217 tableAlphaPlus->LoadData(fileAlphaPlus);
218 tableData[alphaPlus] = tableAlphaPlus;
219
220 // **********************************************************************************************
221
222 helium = heliumDef->GetParticleName();
223 tableFile[helium] = fileHelium;
224
225 lowEnergyLimit[helium] = lowEnergyLimitForA[4];
226 highEnergyLimit[helium] = 400. * MeV;
227
228 // Cross section
229
231 eV,
232 scaleFactor );
233 tableHelium->LoadData(fileHelium);
234 tableData[helium] = tableHelium;
235
236 // **********************************************************************************************
237
238 carbon = carbonDef->GetParticleName();
239 tableFile[carbon] = fileCarbon;
240
241 lowEnergyLimit[carbon] = lowEnergyLimitForA[5] * particle->GetAtomicMass();
242 highEnergyLimit[carbon] = 1e6* particle->GetAtomicMass() * MeV;
243
244 // Cross section
245
247 eV,
248 scaleFactor );
249 tableCarbon->LoadData(fileCarbon);
250 tableData[carbon] = tableCarbon;
251
252 // **********************************************************************************************
253
254 oxygen = oxygenDef->GetParticleName();
255 tableFile[oxygen] = fileOxygen;
256
257 lowEnergyLimit[oxygen] = lowEnergyLimitForA[5]* particle->GetAtomicMass();
258 highEnergyLimit[oxygen] = 1e6* particle->GetAtomicMass()* MeV;
259
260 // Cross section
261
263 eV,
264 scaleFactor );
265 tableOxygen->LoadData(fileOxygen);
266 tableData[oxygen] = tableOxygen;
267
268 // **********************************************************************************************
269
270 nitrogen = nitrogenDef->GetParticleName();
271 tableFile[nitrogen] = fileNitrogen;
272
273 lowEnergyLimit[nitrogen] = lowEnergyLimitForA[5]* particle->GetAtomicMass();
274 highEnergyLimit[nitrogen] = 1e6* particle->GetAtomicMass()* MeV;
275
276 // Cross section
277
279 eV,
280 scaleFactor );
281 tableNitrogen->LoadData(fileNitrogen);
282 tableData[nitrogen] = tableNitrogen;
283
284 // **********************************************************************************************
285
286 iron = ironDef->GetParticleName();
287 tableFile[iron] = fileIron;
288
289 lowEnergyLimit[iron] = lowEnergyLimitForA[5]* particle->GetAtomicMass();
290 highEnergyLimit[iron] = 1e6* particle->GetAtomicMass()* MeV;
291
292 // Cross section
293
295 eV,
296 scaleFactor );
297 tableIron->LoadData(fileIron);
298 tableData[iron] = tableIron;
299
300 // ZF Following lines can be replaced by:
301 SetLowEnergyLimit(lowEnergyLimit[particle->GetParticleName()]);
302 SetHighEnergyLimit(highEnergyLimit[particle->GetParticleName()]);
303 // at least for HZE
304
305 /*
306 if (particle==protonDef)
307 {
308 SetLowEnergyLimit(lowEnergyLimit[proton]);
309 SetHighEnergyLimit(highEnergyLimit[proton]);
310 }
311
312 if (particle==hydrogenDef)
313 {
314 SetLowEnergyLimit(lowEnergyLimit[hydrogen]);
315 SetHighEnergyLimit(highEnergyLimit[hydrogen]);
316 }
317
318 if (particle==heliumDef)
319 {
320 SetLowEnergyLimit(lowEnergyLimit[helium]);
321 SetHighEnergyLimit(highEnergyLimit[helium]);
322 }
323
324 if (particle==alphaPlusDef)
325 {
326 SetLowEnergyLimit(lowEnergyLimit[alphaPlus]);
327 SetHighEnergyLimit(highEnergyLimit[alphaPlus]);
328 }
329
330 if (particle==alphaPlusPlusDef)
331 {
332 SetLowEnergyLimit(lowEnergyLimit[alphaPlusPlus]);
333 SetHighEnergyLimit(highEnergyLimit[alphaPlusPlus]);
334 }
335
336 if (particle==carbonDef)
337 {
338 SetLowEnergyLimit(lowEnergyLimit[carbon]);
339 SetHighEnergyLimit(highEnergyLimit[carbon]);
340 }
341
342 if (particle==oxygenDef)
343 {
344 SetLowEnergyLimit(lowEnergyLimit[oxygen]);
345 SetHighEnergyLimit(highEnergyLimit[oxygen]);
346 }*/
347
348 //----------------------------------------------------------------------
349
350 if( verboseLevel>0 )
351 {
352 G4cout << "Rudd ionisation model is initialized " << G4endl
353 << "Energy range: "
354 << LowEnergyLimit() / eV << " eV - "
355 << HighEnergyLimit() / keV << " keV for "
356 << particle->GetParticleName()
357 << G4endl;
358 }
359
360 // Initialize water density pointer
362
363 //
364
365 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
366
367 if (isInitialised) { return; }
369 isInitialised = true;
370}
371
372//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
373
375 const G4ParticleDefinition* particleDefinition,
376 G4double k,
377 G4double,
378 G4double)
379{
380 if (verboseLevel > 3)
381 G4cout << "Calling CrossSectionPerVolume() of G4DNARuddIonisationExtendedModel" << G4endl;
382
383 // Calculate total cross section for model
384
385 G4DNAGenericIonsManager *instance;
387
388 if (
389 particleDefinition != G4Proton::ProtonDefinition()
390 &&
391 particleDefinition != instance->GetIon("hydrogen")
392 &&
393 particleDefinition != instance->GetIon("alpha++")
394 &&
395 particleDefinition != instance->GetIon("alpha+")
396 &&
397 particleDefinition != instance->GetIon("helium")
398 &&
399 particleDefinition != instance->GetIon("carbon")
400 &&
401 particleDefinition != instance->GetIon("nitrogen")
402 &&
403 particleDefinition != instance->GetIon("oxygen")
404 &&
405 particleDefinition != instance->GetIon("iron")
406 )
407
408 return 0;
409
410 G4double lowLim = 0;
411
412 if ( particleDefinition == G4Proton::ProtonDefinition()
413 || particleDefinition == instance->GetIon("hydrogen")
414 )
415
416 lowLim = lowEnergyLimitOfModelForA[1];
417
418 else if ( particleDefinition == instance->GetIon("alpha++")
419 || particleDefinition == instance->GetIon("alpha+")
420 || particleDefinition == instance->GetIon("helium")
421 )
422
423 lowLim = lowEnergyLimitOfModelForA[4];
424
425 else lowLim = lowEnergyLimitOfModelForA[5];
426
427 G4double highLim = 0;
428 G4double sigma=0;
429
430
431 G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
432
433 if(waterDensity!= 0.0)
434// if (material == nistwater || material->GetBaseMaterial() == nistwater)
435 {
436 const G4String& particleName = particleDefinition->GetParticleName();
437
438 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
439 pos2 = highEnergyLimit.find(particleName);
440
441 if (pos2 != highEnergyLimit.end())
442 {
443 highLim = pos2->second;
444 }
445
446 if (k <= highLim)
447 {
448
449 //SI : XS must not be zero otherwise sampling of secondaries method ignored
450
451 if (k < lowLim) k = lowLim;
452
453 //
454
455 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
456 pos = tableData.find(particleName);
457
458 if (pos != tableData.end())
459 {
460 G4DNACrossSectionDataSet* table = pos->second;
461 if (table != 0)
462 {
463 sigma = table->FindValue(k);
464 }
465 }
466 else
467 {
468 G4Exception("G4DNARuddIonisationExtendedModel::CrossSectionPerVolume","em0002",
469 FatalException,"Model not applicable to particle type.");
470 }
471
472 } // if (k >= lowLim && k < highLim)
473
474 if (verboseLevel > 2)
475 {
476 G4cout << "__________________________________" << G4endl;
477 G4cout << "°°° G4DNARuddIonisationExtendedModel - XS INFO START" << G4endl;
478 G4cout << "°°° Kinetic energy(eV)=" << k/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
479 G4cout << "°°° Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
480 G4cout << "°°° Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
481 // G4cout << " - Cross section per water molecule (cm^-1)=" << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
482 G4cout << "°°° G4DNARuddIonisationExtendedModel - XS INFO END" << G4endl;
483
484 }
485
486 } // if (waterMaterial)
487
488 return sigma*waterDensity;
489// return sigma*material->GetAtomicNumDensityVector()[1];
490
491}
492
493//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
494
495void G4DNARuddIonisationExtendedModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
496 const G4MaterialCutsCouple* /*couple*/,
497 const G4DynamicParticle* particle,
498 G4double,
499 G4double)
500{
501 if (verboseLevel > 3)
502 G4cout << "Calling SampleSecondaries() of G4DNARuddIonisationExtendedModel" << G4endl;
503
504 G4double lowLim = 0;
505 G4double highLim = 0;
506
507 // ZF: the following line summarizes the commented part
508 if(particle->GetDefinition()->GetAtomicMass() <= 4) lowLim = killBelowEnergyForA[particle->GetDefinition()->GetAtomicMass()];
509 else lowLim = killBelowEnergyForA[5]*particle->GetDefinition()->GetAtomicMass();
510
511 /*
512 if(particle->GetDefinition()->GetAtomicMass() >= 5) lowLim = killBelowEnergyForA[5]*particle->GetDefinition()->GetAtomicMass();
513
514
515 if ( particle->GetDefinition() == G4Proton::ProtonDefinition()
516 || particle->GetDefinition() == instance->GetIon("hydrogen")
517 )
518
519 lowLim = killBelowEnergyForA[1];
520
521 if ( particle->GetDefinition() == instance->GetIon("alpha++")
522 || particle->GetDefinition() == instance->GetIon("alpha+")
523 || particle->GetDefinition() == instance->GetIon("helium")
524 )
525
526 lowLim = killBelowEnergyForA[4];*/
527
528 G4double k = particle->GetKineticEnergy();
529
530 const G4String& particleName = particle->GetDefinition()->GetParticleName();
531
532 // SI - the following is useless since lowLim is already defined
533 /*
534 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
535 pos1 = lowEnergyLimit.find(particleName);
536
537 if (pos1 != lowEnergyLimit.end())
538 {
539 lowLim = pos1->second;
540 }
541 */
542
543 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
544 pos2 = highEnergyLimit.find(particleName);
545
546 if (pos2 != highEnergyLimit.end())highLim = pos2->second;
547
548 if (k >= lowLim && k < highLim)
549 {
550 G4ParticleDefinition* definition = particle->GetDefinition();
551 G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
552 /*
553 G4double particleMass = definition->GetPDGMass();
554 G4double totalEnergy = k + particleMass;
555 G4double pSquare = k*(totalEnergy+particleMass);
556 G4double totalMomentum = std::sqrt(pSquare);
557 */
558
559 G4int ionizationShell = RandomSelect(k,particleName);
560
561 // sample deexcitation
562 // here we assume that H_{2}O electronic levels are the same of Oxigen.
563 // this can be considered true with a rough 10% error in energy on K-shell,
564
565 G4int secNumberInit = 0; // need to know at a certain point the enrgy of secondaries
566 G4int secNumberFinal = 0; // So I'll make the diference and then sum the energies
567 G4double bindingEnergy = 0;
568 bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
569
570 if(fAtomDeexcitation) {
571 G4int Z = 8;
573
574 if (ionizationShell <5 && ionizationShell >1)
575 {
576 as = G4AtomicShellEnumerator(4-ionizationShell);
577 }
578 else if (ionizationShell <2)
579 {
581 }
582
583
584 // DEBUG
585 // if (ionizationShell == 4) {
586 //
587 // G4cout << "Z: " << Z << " as: " << as
588 // << " ionizationShell: " << ionizationShell << " bindingEnergy: "<< bindingEnergy/eV << G4endl;
589 // G4cout << "Press <Enter> key to continue..." << G4endl;
590 // G4cin.ignore();
591 // }
592
593
594 const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
595 secNumberInit = fvect->size();
596 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
597 secNumberFinal = fvect->size();
598 }
599
600
601 G4double secondaryKinetic = RandomizeEjectedElectronEnergy(definition,k,ionizationShell);
602
603 G4double cosTheta = 0.;
604 G4double phi = 0.;
605 RandomizeEjectedElectronDirection(definition, k,secondaryKinetic, cosTheta, phi, ionizationShell);
606
607 G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
608 G4double dirX = sinTheta*std::cos(phi);
609 G4double dirY = sinTheta*std::sin(phi);
610 G4double dirZ = cosTheta;
611 G4ThreeVector deltaDirection(dirX,dirY,dirZ);
612 deltaDirection.rotateUz(primaryDirection);
613
614 // Ignored for ions on electrons
615 /*
616 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
617
618 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
619 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
620 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
621 G4double finalMomentum = std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz);
622 finalPx /= finalMomentum;
623 finalPy /= finalMomentum;
624 finalPz /= finalMomentum;
625
626 G4ThreeVector direction;
627 direction.set(finalPx,finalPy,finalPz);
628
629 fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()) ;
630 */
631
633 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
634 G4double deexSecEnergy = 0;
635 for (G4int j=secNumberInit; j < secNumberFinal; j++) {
636
637 deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();
638
639 }
640
642 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy-secondaryKinetic-deexSecEnergy);
643
644 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic) ;
645 fvect->push_back(dp);
646
647 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
649 ionizationShell,
650 theIncomingTrack);
651 }
652
653 // SI - not useful since low energy of model is 0 eV
654
655 if (k < lowLim)
656 {
660 }
661
662}
663
664//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
665
666G4double G4DNARuddIonisationExtendedModel::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
667 G4double k,
668 G4int shell)
669{
670 //-- Fast sampling method -----
671 G4double proposed_energy;
672 G4double random1;
673 G4double value_sampling;
674 G4double max1;
675
676 do
677 {
678 proposed_energy = ProposedSampledEnergy(particleDefinition, k, shell); // Proposed energy by inverse function sampling
679
680 max1=0.;
681
682 for(G4double en=0.; en<20.; en+=1.) if(RejectionFunction(particleDefinition, k, en, shell) > max1)
683 max1=RejectionFunction(particleDefinition, k, en, shell);
684
685 random1 = G4UniformRand()*max1;
686
687 value_sampling = RejectionFunction(particleDefinition, k, proposed_energy, shell);
688
689 } while(random1 > value_sampling);
690
691 return(proposed_energy);
692}
693
694//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
695
696
697void G4DNARuddIonisationExtendedModel::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
698 G4double k,
699 G4double secKinetic,
700 G4double & cosTheta,
701 G4double & phi,
702 G4int shell )
703{
704 G4double maxSecKinetic = 0.;
705 G4double maximumEnergyTransfer = 0.;
706
707 /* if (particleDefinition == G4Proton::ProtonDefinition()
708 || particleDefinition == instance->GetIon("hydrogen"))
709 {
710 if(k/MeV < 100.)maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
711 else {
712 G4double beta2 = 1.-(1.+k*);
713 maxSecKinetic =
714 }
715 }
716
717 if (particleDefinition == instance->GetIon("helium")
718 || particleDefinition == instance->GetIon("alpha+")
719 || particleDefinition == instance->GetIon("alpha++"))
720 {
721 maxSecKinetic = 4.* (0.511 / 3728) * k;
722 }
723
724 if (particleDefinition == G4Carbon::Carbon())
725 {
726 maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k / 12.;
727 }*/
728
729 // ZF. generalized & relativistic version
730
731 if( (k/MeV)/(particleDefinition->GetPDGMass()/MeV) <= 0.1 )
732 {
733 maximumEnergyTransfer= 4.* (electron_mass_c2 / particleDefinition->GetPDGMass()) * k;
734 maximumEnergyTransfer+=waterStructure.IonisationEnergy(shell);
735 }
736 else
737 {
738 G4double approx_nuc_number = particleDefinition->GetPDGMass() / proton_mass_c2;
739 G4double en_per_nucleon = k/approx_nuc_number;
740 G4double beta2 = 1. - 1./pow( (1.+(en_per_nucleon/electron_mass_c2)*(electron_mass_c2/proton_mass_c2)), 2.);
741 G4double gamma = 1./sqrt(1.-beta2);
742 maximumEnergyTransfer = 2.*electron_mass_c2*(gamma*gamma-1.)/(1.+2.*gamma*(electron_mass_c2/particleDefinition->GetPDGMass())+pow(electron_mass_c2/particleDefinition->GetPDGMass(), 2.) );
743 maximumEnergyTransfer+=waterStructure.IonisationEnergy(shell);
744 }
745
746 maxSecKinetic = maximumEnergyTransfer-waterStructure.IonisationEnergy(shell);
747
748 phi = twopi * G4UniformRand();
749
750 if (secKinetic>100*eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic);
751 else cosTheta = (2.*G4UniformRand())-1.;
752
753}
754
755//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
756
757G4double G4DNARuddIonisationExtendedModel::RejectionFunction(G4ParticleDefinition* particleDefinition,
758 G4double k,
759 G4double proposed_ws,
760 G4int ionizationLevelIndex)
761{
762 const G4int j=ionizationLevelIndex;
763 G4double Bj_energy, alphaConst;
764 G4double Ry = 13.6*eV;
765 const G4double Gj[5] = {0.99, 1.11, 1.11, 0.52, 1.};
766
767 // const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV}; //Ding Paper
768
769 // Following values provided by M. Dingfelder (priv. comm)
770 const G4double Bj[5] = {12.60*eV, 14.70*eV, 18.40*eV, 32.20*eV, 540*eV};
771
772 if (j == 4)
773 {
774 alphaConst = 0.66;
775 //---Note that the following (j==4) cases are provided by M. Dingfelder (priv. comm)
776 Bj_energy = waterStructure.IonisationEnergy(ionizationLevelIndex);
777 //---
778 }
779 else
780 {
781 alphaConst = 0.64;
782 Bj_energy = Bj[ionizationLevelIndex];
783 }
784
785 G4double energyTransfer = proposed_ws + Bj_energy;
786 proposed_ws/=Bj_energy;
787 G4DNAGenericIonsManager *instance;
789 G4double tau = 0.;
790 G4double A_ion = 0.;
791 tau = (electron_mass_c2 / particleDefinition->GetPDGMass()) * k;
792 A_ion = particleDefinition->GetAtomicMass();
793
794 G4double v2;
795 G4double beta2;
796
797 if((tau/MeV)<5.447761194e-2)
798 {
799 v2 = tau / Bj_energy;
800 beta2 = 2.*tau / electron_mass_c2;
801 }
802 // Relativistic
803 else
804 {
805 v2 = (electron_mass_c2 / 2. / Bj_energy) * (1. - (1./ pow( (1.+ (tau/electron_mass_c2)),2) ));
806 beta2 =1. - 1./(1.+ (tau/electron_mass_c2/A_ion))/(1.+ (tau/electron_mass_c2/A_ion));
807 }
808
809 G4double v = std::sqrt(v2);
810 G4double wc = 4.*v2 - 2.*v - (Ry/(4.*Bj_energy));
811 G4double rejection_term = 1.+std::exp(alphaConst*(proposed_ws - wc) / v);
812 rejection_term = (1./rejection_term)*CorrectionFactor(particleDefinition,k,ionizationLevelIndex) * Gj[j];
813 //* (S/Bj_energy) ; Not needed anymore
814
815 G4bool isHelium = false;
816
817 if ( particleDefinition == G4Proton::ProtonDefinition()
818 || particleDefinition == instance->GetIon("hydrogen")
819 )
820 {
821 return(rejection_term);
822 }
823
824 else if(particleDefinition->GetAtomicMass() > 4) // anything above Helium
825 {
826 G4double Z = particleDefinition->GetAtomicNumber();
827 G4double x = 100.*std::sqrt(beta2)/std::pow(Z,(2./3.));
828 G4double Zeffion = Z*(1.-std::exp(-1.316*x+0.112*x*x-0.0650*x*x*x));
829 rejection_term*=Zeffion*Zeffion;
830 }
831
832 else if (particleDefinition == instance->GetIon("alpha++") )
833 {
834 isHelium = true;
835 slaterEffectiveCharge[0]=0.;
836 slaterEffectiveCharge[1]=0.;
837 slaterEffectiveCharge[2]=0.;
838 sCoefficient[0]=0.;
839 sCoefficient[1]=0.;
840 sCoefficient[2]=0.;
841 }
842
843 else if (particleDefinition == instance->GetIon("alpha+") )
844 {
845 isHelium = true;
846 slaterEffectiveCharge[0]=2.0;
847 // The following values are provided by M. Dingfelder (priv. comm)
848 slaterEffectiveCharge[1]=2.0;
849 slaterEffectiveCharge[2]=2.0;
850 //
851 sCoefficient[0]=0.7;
852 sCoefficient[1]=0.15;
853 sCoefficient[2]=0.15;
854 }
855
856 else if (particleDefinition == instance->GetIon("helium") )
857 {
858 isHelium = true;
859 slaterEffectiveCharge[0]=1.7;
860 slaterEffectiveCharge[1]=1.15;
861 slaterEffectiveCharge[2]=1.15;
862 sCoefficient[0]=0.5;
863 sCoefficient[1]=0.25;
864 sCoefficient[2]=0.25;
865 }
866
867// if ( particleDefinition == instance->GetIon("helium")
868// || particleDefinition == instance->GetIon("alpha+")
869// || particleDefinition == instance->GetIon("alpha++")
870// )
871 if (isHelium)
872 {
873
874 G4double zEff = particleDefinition->GetPDGCharge() / eplus + particleDefinition->GetLeptonNumber();
875
876 zEff -= ( sCoefficient[0] * S_1s(k, energyTransfer, slaterEffectiveCharge[0], 1.) +
877 sCoefficient[1] * S_2s(k, energyTransfer, slaterEffectiveCharge[1], 2.) +
878 sCoefficient[2] * S_2p(k, energyTransfer, slaterEffectiveCharge[2], 2.) );
879
880 rejection_term*= zEff * zEff;
881 }
882
883 return (rejection_term);
884}
885
886
887//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
888
889
890G4double G4DNARuddIonisationExtendedModel::ProposedSampledEnergy(G4ParticleDefinition* particle,
891 G4double k,
892 G4int ionizationLevelIndex)
893{
894
895 const G4int j=ionizationLevelIndex;
896
897 G4double A1, B1, C1, D1, E1, A2, B2, C2, D2;
898 //G4double alphaConst ;
899 G4double Bj_energy;
900
901 // const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV}; //Ding Paper
902 // Following values provided by M. Dingfelder (priv. comm)
903
904 const G4double Bj[5] = {12.60*eV, 14.70*eV, 18.40*eV, 32.20*eV, 540*eV};
905
906 if (j == 4)
907 {
908 //Data For Liquid Water K SHELL from Dingfelder (Protons in Water)
909 A1 = 1.25;
910 B1 = 0.5;
911 C1 = 1.00;
912 D1 = 1.00;
913 E1 = 3.00;
914 A2 = 1.10;
915 B2 = 1.30;
916 C2 = 1.00;
917 D2 = 0.00;
918 //alphaConst = 0.66;
919 //---Note that the following (j==4) cases are provided by M. Dingfelder (priv. comm)
920 Bj_energy = waterStructure.IonisationEnergy(ionizationLevelIndex);
921 //---
922 }
923 else
924 {
925 //Data For Liquid Water from Dingfelder (Protons in Water)
926 A1 = 1.02;
927 B1 = 82.0;
928 C1 = 0.45;
929 D1 = -0.80;
930 E1 = 0.38;
931 A2 = 1.07;
932 //B2 = 14.6; From Ding Paper
933 // Value provided by M. Dingfelder (priv. comm)
934 B2 = 11.6;
935 //
936 C2 = 0.60;
937 D2 = 0.04;
938 //alphaConst = 0.64;
939
940 Bj_energy = Bj[ionizationLevelIndex];
941 }
942
943 G4double tau = 0.;
944 G4double A_ion = 0.;
945 tau = (electron_mass_c2 / particle->GetPDGMass()) * k;
946
947 A_ion = particle->GetAtomicMass();
948
949 G4double v2;
950 G4double beta2;
951 if((tau/MeV)<5.447761194e-2)
952 {
953 v2 = tau / Bj_energy;
954 beta2 = 2.*tau / electron_mass_c2;
955 }
956 // Relativistic
957 else
958 {
959 v2 = (electron_mass_c2 / 2. / Bj_energy) * (1. - (1./ pow( (1.+ (tau/electron_mass_c2)),2) ));
960 beta2 =1. - 1./(1.+ (tau/electron_mass_c2/A_ion))/(1.+ (tau/electron_mass_c2/A_ion));
961 }
962
963 G4double v = std::sqrt(v2);
964 //G4double wc = 4.*v2 - 2.*v - (Ry/(4.*Bj_energy));
965 G4double L1 = (C1* std::pow(v,(D1))) / (1.+ E1*std::pow(v, (D1+4.)));
966 G4double L2 = C2*std::pow(v,(D2));
967 G4double H1 = (A1*std::log(1.+v2)) / (v2+(B1/v2));
968 G4double H2 = (A2/v2) + (B2/(v2*v2));
969 G4double F1 = L1+H1;
970 G4double F2 = (L2*H2)/(L2+H2);
971
972 // ZF. generalized & relativistic version
973 G4double maximumEnergy;
974
975 //---- maximum kinetic energy , non relativistic ------
976 if( (k/MeV)/(particle->GetPDGMass()/MeV) <= 0.1 )
977 {
978 maximumEnergy = 4.* (electron_mass_c2 / particle->GetPDGMass()) * k;
979 }
980 //---- relativistic -----------------------------------
981 else
982 {
983 G4double gamma = 1./sqrt(1.-beta2);
984 maximumEnergy = 2.*electron_mass_c2*(gamma*gamma-1.)/
985 (1.+2.*gamma*(electron_mass_c2/particle->GetPDGMass())+pow(electron_mass_c2/particle->GetPDGMass(), 2.) );
986 }
987
988 //either it is transfered energy or secondary electron energy ...
989 //maximumEnergy-=Bj_energy;
990
991 //-----------------------------------------------------
992 G4double wmax = maximumEnergy/Bj_energy;
993 G4double c = wmax*(F2*wmax+F1*(2.+wmax))/(2.*(1.+wmax)*(1.+wmax));
994 c=1./c; //!!!!!!!!!!! manual calculus leads to c=1/c
995 G4double randVal = G4UniformRand();
996 G4double proposed_ws = F1*F1*c*c + 2.*F2*c*randVal - 2.*F1*c*randVal;
997 proposed_ws = -F1*c+2.*randVal+std::sqrt(proposed_ws);
998 // proposed_ws = -F1*c+2.*randVal-std::sqrt(proposed_ws);
999 proposed_ws/= ( F1*c + F2*c - 2.*randVal );
1000 proposed_ws*=Bj_energy;
1001
1002 return(proposed_ws);
1003
1004}
1005
1006//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1007
1008G4double G4DNARuddIonisationExtendedModel::S_1s(G4double t,
1009 G4double energyTransferred,
1010 G4double slaterEffectiveChg,
1011 G4double shellNumber)
1012{
1013 // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2)
1014 // Dingfelder, in Chattanooga 2005 proceedings, formula (7)
1015
1016 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
1017 G4double value = 1. - std::exp(-2 * r) * ( ( 2. * r + 2. ) * r + 1. );
1018
1019 return value;
1020}
1021
1022//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1023
1024G4double G4DNARuddIonisationExtendedModel::S_2s(G4double t,
1025 G4double energyTransferred,
1026 G4double slaterEffectiveChg,
1027 G4double shellNumber)
1028{
1029 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4)
1030 // Dingfelder, in Chattanooga 2005 proceedings, formula (8)
1031
1032 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
1033 G4double value = 1. - std::exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
1034
1035 return value;
1036
1037}
1038
1039//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1040
1041G4double G4DNARuddIonisationExtendedModel::S_2p(G4double t,
1042 G4double energyTransferred,
1043 G4double slaterEffectiveChg,
1044 G4double shellNumber)
1045{
1046 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^3 + 2/3 r^4)
1047 // Dingfelder, in Chattanooga 2005 proceedings, formula (9)
1048
1049 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
1050 G4double value = 1. - std::exp(-2 * r) * (((( 2./3. * r + 4./3.) * r + 2.) * r + 2.) * r + 1.);
1051
1052 return value;
1053}
1054
1055//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1056
1057G4double G4DNARuddIonisationExtendedModel::R(G4double t,
1058 G4double energyTransferred,
1059 G4double slaterEffectiveChg,
1060 G4double shellNumber)
1061{
1062 // tElectron = m_electron / m_alpha * t
1063 // Dingfelder, in Chattanooga 2005 proceedings, p 4
1064
1065 G4double tElectron = 0.511/3728. * t;
1066 // The following values are provided by M. Dingfelder (priv. comm)
1067 G4double H = 2.*13.60569172 * eV;
1068 G4double value = std::sqrt ( 2. * tElectron / H ) / ( energyTransferred / H ) * (slaterEffectiveChg/shellNumber);
1069
1070 return value;
1071}
1072
1073//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1074
1075G4double G4DNARuddIonisationExtendedModel::CorrectionFactor(G4ParticleDefinition* particleDefinition, G4double k, G4int shell)
1076{
1077 // ZF Shortened
1078 G4DNAGenericIonsManager *instance;
1080
1081 if (particleDefinition == instance->GetIon("hydrogen") && shell < 4)
1082 {
1083 G4double value = (std::log10(k/eV)-4.2)/0.5;
1084 // The following values are provided by M. Dingfelder (priv. comm)
1085 return((0.6/(1+std::exp(value))) + 0.9);
1086 }
1087 else
1088 {
1089 return(1.);
1090 }
1091}
1092
1093//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1094
1095G4int G4DNARuddIonisationExtendedModel::RandomSelect(G4double k, const G4String& particle )
1096{
1097
1098 G4int level = 0;
1099
1100 // Retrieve data table corresponding to the current particle type
1101
1102 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
1103 pos = tableData.find(particle);
1104
1105 if (pos != tableData.end())
1106 {
1107 G4DNACrossSectionDataSet* table = pos->second;
1108
1109 if (table != 0)
1110 {
1111 G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
1112
1113 const size_t n(table->NumberOfComponents());
1114 size_t i(n);
1115 G4double value = 0.;
1116
1117 while (i>0)
1118 {
1119 i--;
1120 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
1121
1122 value += valuesBuffer[i];
1123 }
1124
1125 value *= G4UniformRand();
1126
1127 i = n;
1128
1129 while (i > 0)
1130 {
1131 i--;
1132
1133 if (valuesBuffer[i] > value)
1134 {
1135 delete[] valuesBuffer;
1136 return i;
1137 }
1138 value -= valuesBuffer[i];
1139 }
1140
1141 if (valuesBuffer) delete[] valuesBuffer;
1142
1143 }
1144 }
1145 else
1146 {
1147 G4Exception("G4DNARuddIonisationExtendedModel::RandomSelect","em0002",
1148 FatalException,"Model not applicable to particle type.");
1149 }
1150
1151 return level;
1152}
1153
1154//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1155
1156G4double G4DNARuddIonisationExtendedModel::PartialCrossSection(const G4Track& track )
1157{
1158 G4double sigma = 0.;
1159
1160 const G4DynamicParticle* particle = track.GetDynamicParticle();
1161 G4double k = particle->GetKineticEnergy();
1162
1163 G4double lowLim = 0;
1164 G4double highLim = 0;
1165
1166 const G4String& particleName = particle->GetDefinition()->GetParticleName();
1167
1168 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
1169 pos1 = lowEnergyLimit.find(particleName);
1170
1171 if (pos1 != lowEnergyLimit.end())
1172 {
1173 lowLim = pos1->second;
1174 }
1175
1176 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
1177 pos2 = highEnergyLimit.find(particleName);
1178
1179 if (pos2 != highEnergyLimit.end())
1180 {
1181 highLim = pos2->second;
1182 }
1183
1184 if (k >= lowLim && k <= highLim)
1185 {
1186 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
1187 pos = tableData.find(particleName);
1188
1189 if (pos != tableData.end())
1190 {
1191 G4DNACrossSectionDataSet* table = pos->second;
1192 if (table != 0)
1193 {
1194 sigma = table->FindValue(k);
1195 }
1196 }
1197 else
1198 {
1199 G4Exception("G4DNARuddIonisationExtendedModel::PartialCrossSection","em0002",
1200 FatalException,"Model not applicable to particle type.");
1201 }
1202 }
1203
1204 return sigma;
1205}
1206
1207//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1208
1209G4double G4DNARuddIonisationExtendedModel::Sum(G4double /* energy */, const G4String& /* particle */)
1210{
1211 return 0;
1212}
1213
G4AtomicShellEnumerator
@ eIonizedMolecule
@ FatalException
@ fStopAndKill
double G4double
Definition: G4Types.hh:64
int G4int
Definition: G4Types.hh:66
bool G4bool
Definition: G4Types.hh:67
#define G4endl
Definition: G4ios.hh:52
G4DLLIMPORT std::ostream G4cout
#define C1
#define G4UniformRand()
Definition: Randomize.hh:53
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:72
static G4DNAChemistryManager * Instance()
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
virtual G4double FindValue(G4double e, G4int componentId=0) const
virtual size_t NumberOfComponents(void) const
virtual const G4VEMDataSet * GetComponent(G4int componentId) const
virtual G4bool LoadData(const G4String &argFileName)
static G4DNAGenericIonsManager * Instance(void)
G4ParticleDefinition * GetIon(const G4String &name)
static G4DNAMolecularMaterial * Instance()
const std::vector< double > * GetNumMolPerVolTableFor(const G4Material *) const
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
G4DNARuddIonisationExtendedModel(const G4ParticleDefinition *p=0, const G4String &nam="DNARuddIonisationExtendedModel")
const G4ThreeVector & GetMomentumDirection() const
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
static G4Electron * Electron()
Definition: G4Electron.cc:94
static G4LossTableManager * Instance()
G4VAtomDeexcitation * AtomDeexcitation()
size_t GetIndex() const
Definition: G4Material.hh:261
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:576
const G4Track * GetCurrentTrack() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4int GetAtomicNumber() const
G4int GetAtomicMass() const
G4double GetPDGCharge() const
const G4String & GetParticleName() const
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
const G4DynamicParticle * GetDynamicParticle() const
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:585
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:109
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:529
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:522
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:592
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:641
void ProposeTrackStatus(G4TrackStatus status)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41