Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // 27 // G4MicroElecInelasticModel_new.cc, 2011/08/2 28 // 2020/05/20 P. Caro 29 // Q. Gibaru is with CEA [a] 30 // M. Raine and D. Lambert a 31 // 32 // A part of this work has been funded by the 33 // [a] CEA, DAM, DIF - 91297 ARPAJON, France 34 // [b] ONERA - DPHY, 2 avenue E.Belin, 31055 T 35 // [c] CNES, 18 av.E.Belin, 31401 Toulouse CED 36 // 37 // Based on the following publications 38 // - A.Valentin, M. Raine, 39 // Inelastic cross-sections of low energy e 40 // for the simulation of heavy ion trac 41 // NSS Conf. Record 2010, pp. 80-85 42 // https://doi.org/10.1109/NSSMIC. 43 // 44 // - A.Valentin, M. Raine, M.Gaillardin, 45 // Geant4 physics processes for microdo 46 // very low energy electromagnetic mode 47 // https://doi.org/10.1016/j.nimb. 48 // NIM B, vol. 288, pp. 66-73, 2012, pa 49 // heavy ions in Si, NIM B, vol. 287, p 50 // https://doi.org/10.1016/j.nimb. 51 // 52 // - M. Raine, M. Gaillardin, P. Paillet 53 // Geant4 physics processes for silicon 54 // Improvements and extension of the en 55 // NIM B, vol. 325, pp. 97-100, 2014 56 // https://doi.org/10.1016/j.nimb. 57 // 58 // - J. Pierron, C. Inguimbert, M. Belhaj 59 // Electron emission yield for low ener 60 // Monte Carlo simulation and experimen 61 // Journal of Applied Physics 121 (2017 62 // https://doi.org/10.1063/1.498 63 // 64 // - P. Caron, 65 // Study of Electron-Induced Single-Eve 66 // PHD, 16th October 2019 67 // 68 // - Q.Gibaru, C.Inguimbert, P.Caron, M.Raine 69 // Geant4 physics processes for microdo 70 // Extension of MicroElec to very low e 71 // NIM B, 2020, in review. 72 // 73 // 74 //....oooOO0OOooo........oooOO0OOooo........oo 75 76 77 #include "globals.hh" 78 #include "G4MicroElecInelasticModel_new.hh" 79 #include "G4PhysicalConstants.hh" 80 #include "G4SystemOfUnits.hh" 81 #include "G4ios.hh" 82 #include "G4UnitsTable.hh" 83 #include "G4UAtomicDeexcitation.hh" 84 #include "G4LossTableManager.hh" 85 #include "G4ionEffectiveCharge.hh" 86 #include "G4MicroElecMaterialStructure.hh" 87 #include "G4DeltaAngle.hh" 88 89 #include <sstream> 90 91 //....oooOO0OOooo........oooOO0OOooo........oo 92 93 using namespace std; 94 95 //....oooOO0OOooo........oooOO0OOooo........oo 96 97 G4MicroElecInelasticModel_new::G4MicroElecInel 98 const G4ParticleDefinition*, const G4St 99 :G4VEmModel(nam),isInitialised(false) 100 { 101 102 verboseLevel= 0; 103 // Verbosity scale: 104 // 0 = nothing 105 // 1 = warning for energy non-conservation 106 // 2 = details of energy budget 107 // 3 = calculation of cross sections, file o 108 // 4 = entering in methods 109 110 if( verboseLevel>0 ) 111 { 112 G4cout << "MicroElec inelastic model is co 113 } 114 115 //Mark this model as "applicable" for atomic 116 SetDeexcitationFlag(true); 117 fAtomDeexcitation = nullptr; 118 fParticleChangeForGamma = nullptr; 119 120 // default generator 121 SetAngularDistribution(new G4DeltaAngle()); 122 123 // Selection of computation method 124 fasterCode = true; 125 SEFromFermiLevel = false; 126 } 127 128 //....oooOO0OOooo........oooOO0OOooo........oo 129 130 G4MicroElecInelasticModel_new::~G4MicroElecIne 131 { 132 // Cross section 133 // (0) 134 for (auto const& p : tableTCS) { 135 MapData* tableData = p.second; 136 for (auto const& pos : *tableData) { delet 137 delete tableData; 138 } 139 tableTCS.clear(); 140 141 // (1) 142 for (auto const & obj : eDiffDatatable) { 143 auto ptr = obj.second; 144 ptr->clear(); 145 delete ptr; 146 } 147 148 for (auto const & obj : pDiffDatatable) { 149 auto ptr = obj.second; 150 ptr->clear(); 151 delete ptr; 152 } 153 154 // (2) 155 for (auto const & obj : eNrjTransStorage) { 156 delete obj.second; 157 } 158 for (auto const & obj : pNrjTransStorage) { 159 delete obj.second; 160 } 161 162 // (3) 163 for (auto const& p : eProbaShellStorage) { 164 delete p.second; 165 } 166 167 for (auto const& p : pProbaShellStorage) { 168 delete p.second; 169 } 170 171 // (4) 172 for (auto const& p : eIncidentEnergyStorage) 173 delete p.second; 174 } 175 176 for (auto const& p : pIncidentEnergyStorage) 177 delete p.second; 178 } 179 180 // (5) 181 for (auto const& p : eVecmStorage) { 182 delete p.second; 183 } 184 185 for (auto const& p : pVecmStorage) { 186 delete p.second; 187 } 188 189 // (6) 190 for (auto const& p : tableMaterialsStructure 191 delete p.second; 192 } 193 } 194 195 //....oooOO0OOooo........oooOO0OOooo........oo 196 197 void G4MicroElecInelasticModel_new::Initialise 198 const G4DataVector& /*cuts*/) 199 { 200 if (isInitialised) { return; } 201 202 if (verboseLevel > 3) 203 G4cout << "Calling G4MicroElecInelasticMod 204 205 const char* path = G4FindDataDir("G4LEDATA") 206 if (!path) 207 { 208 G4Exception("G4MicroElecElasticModel_new 209 return; 210 } 211 212 G4String modelName = "mermin"; 213 G4cout << "****************************" << 214 G4cout << modelName << " model loaded !" << 215 216 // Energy limits 217 G4ParticleDefinition* electronDef = G4Electr 218 G4ParticleDefinition* protonDef = G4Proton:: 219 G4String electron = electronDef->GetParticle 220 G4String proton = protonDef->GetParticleName 221 222 G4double scaleFactor = 1.0; 223 224 // *** ELECTRON 225 lowEnergyLimit[electron] = 2 * eV; 226 highEnergyLimit[electron] = 10.0 * MeV; 227 228 // Cross section 229 G4ProductionCutsTable* theCoupleTable = G4Pr 230 G4int numOfCouples = (G4int)theCoupleTable-> 231 232 for (G4int i = 0; i < numOfCouples; ++i) { 233 const G4Material* material = theCoupleTabl 234 G4cout << "Material " << i + 1 << " / " << 235 if (material->GetName() == "Vacuum") conti 236 G4String mat = material->GetName().substr( 237 MapData* tableData = new MapData; 238 currentMaterialStructure = new G4MicroElec 239 240 tableMaterialsStructures[mat] = currentMat 241 if (particle == electronDef) { 242 //TCS 243 G4String fileElectron("Inelastic/" + mod 244 G4cout << fileElectron << G4endl; 245 G4MicroElecCrossSectionDataSet_new* tabl 246 tableE->LoadData(fileElectron); 247 tableData->insert(make_pair(electron, ta 248 249 // DCS 250 std::ostringstream eFullFileName; 251 if (fasterCode) { 252 eFullFileName << path << "/microelec/Inelast 253 G4cout << "Faster code = true" << G4endl; 254 G4cout << "Inelastic/cumulated_" + modelName 255 } 256 else { 257 eFullFileName << path << "/microelec/Inelast 258 G4cout << "Faster code = false" << G4endl; 259 G4cout << "Inelastic/" + modelName + "_sigma 260 } 261 262 std::ifstream eDiffCrossSection(eFullFil 263 if (!eDiffCrossSection) 264 { 265 std::stringstream ss; 266 ss << "Missing data " << eFullFileName.str 267 std::string sortieString = ss.str(); 268 269 if (fasterCode) G4Exception("G4MicroElecIn 270 FatalException, sortieString.c_s 271 else { 272 G4Exception("G4MicroElecInelasticModel_n 273 FatalException, "Missing data file:/micr 274 } 275 } 276 277 // Clear the arrays for re-initializatio 278 // Octobre 22nd, 2014 - Melanie Raine 279 //Creating vectors of maps for DCS and C 280 //Each vector is storing one map for eac 281 G4bool isUsed1 = false; 282 vector<TriDimensionMap>* eDiffCrossSecti 283 new vector<TriDimensionMap>; //Storage of [I 284 vector<TriDimensionMap>* eNrjTransfData 285 new vector<TriDimensionMap>; //Storage of po 286 vector<VecMap>* eProbaShellMap = new vec 287 vector<G4double>* eTdummyVec = new vecto 288 VecMap* eVecm = new VecMap; //Transfered 289 290 for (G4int j = 0; j < currentMaterialStr 291 { 292 eDiffCrossSectionData->push_back(TriDimens 293 eNrjTransfData->push_back(TriDimensionMap( 294 eProbaShellMap->push_back(VecMap()); 295 } 296 297 eTdummyVec->push_back(0.); 298 while (!eDiffCrossSection.eof()) 299 { 300 G4double tDummy; //incident energy 301 G4double eDummy; //transfered energy 302 eDiffCrossSection >> tDummy >> eDummy; 303 if (tDummy != eTdummyVec->back()) eTdummyV 304 305 G4double tmp; //probability 306 for (G4int j = 0; j < currentMaterialStruc 307 { 308 eDiffCrossSection >> tmp; 309 (*eDiffCrossSectionData)[j][tDummy][eD 310 311 if (fasterCode) 312 { 313 (*eNrjTransfData)[j][tDummy][(*eDiffCros 314 (*eProbaShellMap)[j][tDummy].push_back(( 315 } 316 else { // SI - only if eof is not rea 317 (*eDiffCrossSectionData)[j][tDummy][eDummy 318 (*eVecm)[tDummy].push_back(eDummy); 319 } 320 } 321 } 322 // 323 //G4cout << "add to material vector" << 324 325 //Filing maps for the current material i 326 if (fasterCode) { 327 isUsed1 = true; 328 eNrjTransStorage[mat] = eNrjTransfData; 329 eProbaShellStorage[mat] = eProbaShellMap; 330 } 331 else { 332 eDiffDatatable[mat] = eDiffCrossSectionData; 333 eVecmStorage[mat] = eVecm; 334 } 335 eIncidentEnergyStorage[mat] = eTdummyVec 336 337 //Cleanup support vectors 338 if (!isUsed1) { 339 delete eProbaShellMap; 340 delete eNrjTransfData; 341 } else { 342 delete eDiffCrossSectionData; 343 delete eVecm; 344 } 345 } 346 347 // *** PROTON 348 if (particle == protonDef) 349 { 350 // Cross section 351 G4String fileProton("Inelastic/" + modelName 352 G4MicroElecCrossSectionDataSet_new* tableP = 353 tableP->LoadData(fileProton); 354 tableData->insert(make_pair(proton, tableP)) 355 356 // DCS 357 std::ostringstream pFullFileName; 358 if (fasterCode) { 359 pFullFileName << path << "/microelec/Inela 360 G4cout << "Faster code = true" << G4endl; 361 G4cout << "Inelastic/cumulated_" + modelNa 362 } 363 else { 364 pFullFileName << path << "/microelec/Inela 365 G4cout << "Faster code = false" << G4endl; 366 G4cout << "Inelastic/" + modelName + "_sig 367 } 368 369 std::ifstream pDiffCrossSection(pFullFileNam 370 if (!pDiffCrossSection) 371 { 372 if (fasterCode) G4Exception("G4MicroElec 373 FatalException, "Missing data file:/ 374 else { 375 G4Exception("G4MicroElecInelasticModel 376 FatalException, "Missing data file:/mi 377 } 378 } 379 380 // 381 // Clear the arrays for re-initialization ca 382 // Octobre 22nd, 2014 - Melanie Raine 383 //Creating vectors of maps for DCS and Cumul 384 //Each vector is storing one map for each sh 385 386 G4bool isUsed1 = false; 387 vector<TriDimensionMap>* pDiffCrossSectionDa 388 new vector<TriDimensionMap>; //Storage of 389 vector<TriDimensionMap>* pNrjTransfData = 390 new vector<TriDimensionMap>; //Storage of 391 vector<VecMap>* pProbaShellMap = 392 new vector<VecMap>; //Storage of the vecto 393 vector<G4double>* pTdummyVec = 394 new vector<G4double>; //Storage of inciden 395 VecMap* pVecm = new VecMap; //Transfered ene 396 397 for (G4int j = 0; j < currentMaterialStructu 398 //Filling the map vectors with an empty ma 399 { 400 pDiffCrossSectionData->push_back(TriDime 401 pNrjTransfData->push_back(TriDimensionMa 402 pProbaShellMap->push_back(VecMap()); 403 } 404 405 pTdummyVec->push_back(0.); 406 while (!pDiffCrossSection.eof()) 407 { 408 G4double tDummy; //incident energy 409 G4double eDummy; //transfered energy 410 pDiffCrossSection >> tDummy >> eDummy; 411 if (tDummy != pTdummyVec->back()) pTdumm 412 413 G4double tmp; //probability 414 for (G4int j = 0; j < currentMaterialStr 415 { 416 pDiffCrossSection >> tmp; 417 (*pDiffCrossSectionData)[j][tDummy][eDummy 418 // ArrayofMaps[j] -> fill with 3DMap(incid 419 // 2Dmap (transferedEnergy,proba=tmp) ) -> 420 // with proba for transfered energy eDummy 421 422 if (fasterCode) 423 { 424 (*pNrjTransfData)[j][tDummy][(*pDiffCr 425 (*pProbaShellMap)[j][tDummy].push_back 426 } 427 else { // SI - only if eof is not reached 428 (*pDiffCrossSectionData)[j][tDummy][eDum 429 (*pVecm)[tDummy].push_back(eDummy); 430 } 431 } 432 } 433 434 //Filing maps for the current material into 435 if (fasterCode) { 436 isUsed1 = true; 437 pNrjTransStorage[mat] = pNrjTransfData; 438 pProbaShellStorage[mat] = pProbaShellMap; 439 } 440 else { 441 pDiffDatatable[mat] = pDiffCrossSectionDat 442 pVecmStorage[mat] = pVecm; 443 } 444 pIncidentEnergyStorage[mat] = pTdummyVec; 445 446 //Cleanup support vectors 447 if (!isUsed1) { 448 delete pProbaShellMap; 449 delete pNrjTransfData; 450 } else { 451 delete pDiffCrossSectionData; 452 delete pVecm; 453 } 454 } 455 tableTCS[mat] = tableData; 456 } 457 if (particle==electronDef) 458 { 459 SetLowEnergyLimit(lowEnergyLimit[electro 460 SetHighEnergyLimit(highEnergyLimit[elect 461 } 462 463 if (particle==protonDef) 464 { 465 SetLowEnergyLimit(100*eV); 466 SetHighEnergyLimit(10*MeV); 467 } 468 469 if( verboseLevel>1 ) 470 { 471 G4cout << "MicroElec Inelastic model is 472 << "Energy range: " 473 << LowEnergyLimit() / keV << " keV - " 474 << HighEnergyLimit() / MeV << " MeV for 475 << particle->GetParticleName() 476 << " with mass (amu) " << particle->Get 477 << " and charge " << particle->GetPDGCh 478 << G4endl << G4endl ; 479 } 480 481 fAtomDeexcitation = G4LossTableManager::Ins 482 483 fParticleChangeForGamma = GetParticleChangeF 484 isInitialised = true; 485 } 486 487 //....oooOO0OOooo........oooOO0OOooo........oo 488 489 G4double G4MicroElecInelasticModel_new::CrossS 490 const G4ParticleDefinition* par 491 G4double ekin, 492 G4double, 493 G4double) 494 { 495 if (verboseLevel > 3) G4cout << "Calling Cro 496 497 G4double density = material->GetTotNbOfAtoms 498 currentMaterial = material->GetName().substr 499 500 MapStructure::iterator structPos; 501 structPos = tableMaterialsStructures.find(cu 502 503 // Calculate total cross section for model 504 TCSMap::iterator tablepos; 505 tablepos = tableTCS.find(currentMaterial); 506 507 if (tablepos == tableTCS.end() ) 508 { 509 G4String str = "Material "; 510 str += currentMaterial + " TCS Table not 511 G4Exception("G4MicroElecInelasticModel_n 512 return 0; 513 } 514 else if(structPos == tableMaterialsStructure 515 { 516 G4String str = "Material "; 517 str += currentMaterial + " Structure not 518 G4Exception("G4MicroElecInelasticModel_n 519 return 0; 520 } 521 else { 522 MapData* tableData = tablepos->second; 523 currentMaterialStructure = structPos->seco 524 525 G4double sigma = 0; 526 527 const G4String& particleName = particleDef 528 G4String nameLocal = particleName; 529 G4int pdg = particleDefinition->GetPDGEnco 530 G4int Z = particleDefinition->GetAtomicNum 531 532 G4double Zeff = 1.0, Zeff2 = Zeff*Zeff; 533 G4double Mion_c2 = particleDefinition->Get 534 535 if (Mion_c2 > proton_mass_c2) 536 { 537 ekin *= proton_mass_c2 / Mion_c2; 538 nameLocal = "proton"; 539 } 540 541 G4double lowLim = currentMaterialStructure 542 G4double highLim = currentMaterialStructur 543 544 if (ekin >= lowLim && ekin < highLim) 545 { 546 std::map< G4String, G4MicroElecCrossSectionD 547 pos = tableData->find(nameLocal); //find par 548 549 if (pos != tableData->end()) 550 { 551 G4MicroElecCrossSectionDataSet_new* tabl 552 if (table != 0) 553 { 554 sigma = table->FindValue(ekin); 555 556 if (Mion_c2 > proton_mass_c2) { 557 sigma = 0.; 558 for (G4int i = 0; i < currentMaterialStr 559 Zeff = BKZ(ekin / (proton_mass_c2 / Mi 560 Zeff2 = Zeff*Zeff; 561 sigma += Zeff2*table->FindShellValue(e 562 // il faut utiliser le ekin mis à l'echelle p 563 // valeur dans les tables proton 564 565 } 566 } 567 else { 568 sigma = table->FindValue(ekin); 569 } 570 } 571 } 572 else 573 { 574 G4Exception("G4MicroElecInelasticModel_n 575 "em0002", FatalExcepti 576 "Model not applicable 577 } 578 } 579 else 580 { 581 return 1 / DBL_MAX; 582 } 583 584 if (verboseLevel > 3) 585 { 586 G4cout << "---> Kinetic energy (eV)=" << eki 587 G4cout << " - Cross section per Si atom (cm^ 588 G4cout << " - Cross section per Si atom (cm^ 589 } 590 591 return (sigma)*density; 592 } 593 } 594 595 //....oooOO0OOooo........oooOO0OOooo........oo 596 597 void G4MicroElecInelasticModel_new::SampleSeco 598 599 600 601 602 { 603 604 if (verboseLevel > 3) 605 G4cout << "Calling SampleSecondaries() of 606 607 G4int pdg = particle->GetParticleDefinition( 608 G4double lowLim = currentMaterialStructure-> 609 G4double highLim = currentMaterialStructure- 610 611 G4double ekin = particle->GetKineticEnergy() 612 G4double k = ekin ; 613 614 G4ParticleDefinition* PartDef = particle->Ge 615 const G4String& particleName = PartDef->GetP 616 G4String nameLocal2 = particleName ; 617 G4double particleMass = particle->GetDefinit 618 G4double originalMass = particleMass; // a p 619 G4int originalZ = particle->GetDefinition()- 620 621 if (particleMass > proton_mass_c2) 622 { 623 k *= proton_mass_c2/particleMass ; 624 PartDef = G4Proton::ProtonDefinition(); 625 nameLocal2 = "proton" ; 626 } 627 628 if (k >= lowLim && k < highLim) 629 { 630 G4ParticleMomentum primaryDirection = pa 631 G4double totalEnergy = ekin + particleMa 632 G4double pSquare = ekin * (totalEnergy + 633 G4double totalMomentum = std::sqrt(pSqua 634 635 G4int Shell = 1; 636 637 Shell = RandomSelect(k,nameLocal2,origin 638 639 G4double bindingEnergy = currentMaterial 640 G4double limitEnergy = currentMaterialSt 641 642 G4bool weaklyBound = currentMaterialStruct 643 644 if (verboseLevel > 3) 645 { 646 G4cout << "---> Kinetic energy (eV)=" << k 647 G4cout << "Shell: " << Shell << ", energy: 648 } 649 650 651 if (k<limitEnergy) { 652 if (weaklyBound && k > currentMaterialStru 653 limitEnergy = currentMaterialStructure-> 654 } 655 else return; } 656 657 // sample deexcitation 658 659 std::size_t secNumberInit = 0; // need 660 std::size_t secNumberFinal = 0; // So I' 661 662 // G4cout << currentMaterial << G4endl; 663 G4int Z = currentMaterialStructure->GetZ 664 G4int shellEnum = currentMaterialStructu 665 if (currentMaterialStructure->IsShellWea 666 667 if(fAtomDeexcitation && shellEnum >=0) 668 { 669 // G4cout << "enter if deex and shell 0 670 G4AtomicShellEnumerator as = G4AtomicShell 671 const G4AtomicShell* shell = fAtomDeexcita 672 secNumberInit = fvect->size(); 673 fAtomDeexcitation->GenerateParticles(fvect 674 secNumberFinal = fvect->size(); 675 } 676 677 G4double secondaryKinetic=-1000*eV; 678 SEFromFermiLevel = false; 679 if (!fasterCode) 680 { 681 secondaryKinetic = RandomizeEjectedElectro 682 } 683 else 684 { 685 secondaryKinetic = RandomizeEjectedElectro 686 } 687 688 if (verboseLevel > 3) 689 { 690 G4cout << "Ionisation process" << G4endl; 691 G4cout << "Shell: " << Shell << " Kin. ene 692 << " Sec. energy (eV)=" << secondaryKinet 693 } 694 G4ThreeVector deltaDirection = 695 GetAngularDistribution()->SampleDirectionFor 696 Z, Shell, 697 couple->GetMaterial()); 698 699 if (particle->GetDefinition() == G4Elect 700 { 701 G4double deltaTotalMomentum = std::sqrt(se 702 703 G4double finalPx = totalMomentum*primaryDi 704 G4double finalPy = totalMomentum*primaryDi 705 G4double finalPz = totalMomentum*primaryDi 706 G4double finalMomentum = std::sqrt(finalPx 707 finalPx /= finalMomentum; 708 finalPy /= finalMomentum; 709 finalPz /= finalMomentum; 710 711 G4ThreeVector direction; 712 direction.set(finalPx,finalPy,finalPz); 713 714 fParticleChangeForGamma->ProposeMomentumDi 715 } 716 else fParticleChangeForGamma->ProposeMom 717 718 // note that secondaryKinetic is the ene 719 G4double deexSecEnergy = 0; 720 for (std::size_t j=secNumberInit; j < se 721 deexSecEnergy = deexSecEnergy + (*fvec 722 } 723 // correction CI 12/01/2023 limit energy 724 //if (SEFromFermiLevel) limitEnergy = curr 725 //fParticleChangeForGamma->SetProposedKi 726 //fParticleChangeForGamma->ProposeLocalE 727 728 // correction CI 09/03/2022 limit energy 729 //if (!SEFromFermiLevel && weaklyBound) li 730 //if (!SEFromFermiLevel && weaklyBound) li 731 fParticleChangeForGamma->SetProposedKineti 732 fParticleChangeForGamma->ProposeLocalEnerg 733 734 if (secondaryKinetic>0) 735 { 736 G4DynamicParticle* dp = new G4DynamicParti 737 fvect->push_back(dp); 738 } 739 } 740 } 741 742 //....oooOO0OOooo........oooOO0OOooo........oo 743 744 G4double G4MicroElecInelasticModel_new::Random 745 const G4ParticleDefinition* particleD 746 G4double k, G4int shell, G4double originalM 747 { 748 G4double secondaryElectronKineticEnergy=0.; 749 if (particleDefinition == G4Electron::Electr 750 { 751 G4double maximumEnergyTransfer=k; 752 G4double crossSectionMaximum = 0.; 753 G4double minEnergy = currentMaterialStru 754 G4double maxEnergy = maximumEnergyTransf 755 G4int nEnergySteps = 100; 756 757 G4double value(minEnergy); 758 G4double stpEnergy(std::pow(maxEnergy/va 759 G4int step(nEnergySteps); 760 while (step>0) 761 { 762 --step; 763 G4double differentialCrossSection = 764 DifferentialCrossSection(particleDefinit 765 crossSectionMaximum = std::max(crossSectio 766 value*=stpEnergy; 767 } 768 769 do 770 { 771 secondaryElectronKineticEnergy = G4Uniform 772 (maximumEnergyTransfer-currentMaterialSt 773 } while(G4UniformRand()*crossSectionMaximum 774 DifferentialCrossSection(particleDefinitio 775 (secondaryElectronKineticEnergy+currentM 776 // added 12/01/2023 777 return secondaryElectronKineticEnergy; 778 } 779 else if (particleDefinition == G4Proton::Pro 780 { 781 G4double maximumEnergyTransfer = 782 ComputeElasticQmax(k/(proton_mass_c2/origina 783 currentMaterialStructure->Energy(shel 784 originalMass/c_squared, electron_mass 785 786 G4double crossSectionMaximum = 0.; 787 788 G4double minEnergy = currentMaterialStru 789 G4double maxEnergy = maximumEnergyTransf 790 G4int nEnergySteps = 100; 791 792 G4double value(minEnergy); 793 G4double stpEnergy(std::pow(maxEnergy/va 794 G4int step(nEnergySteps); 795 796 while (step>0) 797 { 798 --step; 799 G4double differentialCrossSection = 800 DifferentialCrossSection(particleDefinit 801 crossSectionMaximum = std::max(crossSectio 802 value*=stpEnergy; 803 } 804 805 G4double energyTransfer = 0.; 806 do 807 { 808 energyTransfer = G4UniformRand() * maximum 809 } while(G4UniformRand()*crossSectionMaximum 810 DifferentialCrossSection(particleDefinitio 811 812 secondaryElectronKineticEnergy = 813 energyTransfer-currentMaterialStructure->Get 814 815 } 816 return std::max(secondaryElectronKineticEner 817 } 818 819 //....oooOO0OOooo........oooOO0OOooo........oo 820 821 G4double G4MicroElecInelasticModel_new::Random 822 const G4ParticleDefinition* particleD 823 { 824 G4double secondaryElectronKineticEnergy = 0. 825 G4double random = G4UniformRand(); 826 G4bool weaklyBound = currentMaterialStructur 827 G4double transf = TransferedEnergy(particleD 828 if (!weaklyBound) { 829 secondaryElectronKineticEnergy = transf - 830 if(secondaryElectronKineticEnergy <= 0.) { 831 secondaryElectronKineticEnergy = 0.0; 832 } 833 } 834 else { 835 secondaryElectronKineticEnergy = transf - 836 // for weaklybound electrons = gap + aver 837 if (secondaryElectronKineticEnergy <= 0.) 838 secondaryElectronKineticEnergy = 0.0; 839 SEFromFermiLevel = true; 840 } 841 } 842 //corrections CI 07/02/2022 - added 843 return secondaryElectronKineticEnergy; 844 } 845 846 //....oooOO0OOooo........oooOO0OOooo........oo 847 848 G4double G4MicroElecInelasticModel_new::Transf 849 const G4ParticleDefinition* particleDefinit 850 G4double k, 851 G4int ionizationLevelIndex, 852 G4double random) 853 { 854 G4double nrj = 0.; 855 G4double valueK1 = 0; 856 G4double valueK2 = 0; 857 G4double valuePROB21 = 0; 858 G4double valuePROB22 = 0; 859 G4double valuePROB12 = 0; 860 G4double valuePROB11 = 0; 861 G4double nrjTransf11 = 0; 862 G4double nrjTransf12 = 0; 863 G4double nrjTransf21 = 0; 864 G4double nrjTransf22 = 0; 865 866 G4double maximumEnergyTransfer1 = 0; 867 G4double maximumEnergyTransfer2 = 0; 868 G4double maximumEnergyTransferP = 4.* (elect 869 G4double bindingEnergy = currentMaterialStru 870 871 if (particleDefinition == G4Electron::Electr 872 { 873 dataDiffCSMap::iterator iterator_Nrj; 874 iterator_Nrj = eNrjTransStorage.find(cur 875 876 dataProbaShellMap::iterator iterator_Pro 877 iterator_Proba = eProbaShellStorage.find 878 879 incidentEnergyMap::iterator iterator_Tdu 880 iterator_Tdummy = eIncidentEnergyStorage 881 882 if(iterator_Nrj == eNrjTransStorage.end( 883 iterator_Tdummy == eIncidentEnergyStorage.e 884 { 885 G4String str = "Material "; 886 str += currentMaterial + " not found!"; 887 G4Exception("G4MicroElecInelasticModel_new 888 FatalException, str); 889 } 890 else { 891 vector<TriDimensionMap>* eNrjTransfData = it 892 vector<VecMap>* eProbaShellMap = iterator_Pr 893 vector<G4double>* eTdummyVec = iterator_Tdum 894 895 // k should be in eV auto, :std::vector<doub 896 auto k2 = std::upper_bound(eTdummyVec->begin 897 eTdummyVec->end(), 898 k); 899 auto k1 = k2 - 1; 900 901 // SI : the following condition avoids situa 902 if (random <= (*eProbaShellMap)[ionizationLe 903 && random <= (*eProbaShellMap)[ionizatio 904 { 905 auto prob12 = 906 std::upper_bound((*eProbaShellMap)[ion 907 (*eProbaShellMap)[ionizationLevel 908 random); 909 910 auto prob11 = prob12 - 1; 911 912 auto prob22 = 913 std::upper_bound((*eProbaShellMap)[ion 914 (*eProbaShellMap)[ionizationLevel 915 random); 916 917 auto prob21 = prob22 - 1; 918 919 valueK1 = *k1; 920 valueK2 = *k2; 921 valuePROB21 = *prob21; 922 valuePROB22 = *prob22; 923 valuePROB12 = *prob12; 924 valuePROB11 = *prob11; 925 926 // The following condition avoid getting 927 if (valuePROB11 == 0) nrjTransf11 = bind 928 else nrjTransf11 = (*eNrjTransfData)[ion 929 if (valuePROB12 == 1) 930 { 931 if ((valueK1 + bindingEnergy) / 2. > value 932 maximumEnergyTransfer1 = valueK1; 933 else 934 maximumEnergyTransfer1 = (valueK1 + bind 935 936 nrjTransf12 = maximumEnergyTransfer1; 937 } 938 else 939 nrjTransf12 = (*eNrjTransfData)[ioniza 940 941 if (valuePROB21 == 0) nrjTransf21 = bind 942 else nrjTransf21 = (*eNrjTransfData)[ion 943 if (valuePROB22 == 1) 944 { 945 if ((valueK2 + bindingEnergy) / 2. > value 946 else maximumEnergyTransfer2 = (valueK2 + b 947 948 nrjTransf22 = maximumEnergyTransfer2; 949 } 950 else nrjTransf22 = (*eNrjTransfData)[ion 951 952 } 953 // Avoids cases where cum xs is zero for k1 954 if (random > (*eProbaShellMap)[ionizationLev 955 { 956 auto prob22 = 957 std::upper_bound((*eProbaShellMap)[ion 958 (*eProbaShellMap)[ionizationLevel 959 random); 960 auto prob21 = prob22 - 1; 961 962 valueK1 = *k1; 963 valueK2 = *k2; 964 valuePROB21 = *prob21; 965 valuePROB22 = *prob22; 966 967 nrjTransf21 = (*eNrjTransfData)[ionizati 968 nrjTransf22 = (*eNrjTransfData)[ionizati 969 970 G4double interpolatedvalue2 = Interpolat 971 valuePROB22, 972 random, 973 nrjTransf21, 974 nrjTransf22); 975 976 // zeros are explicitly set 977 G4double value = Interpolate(valueK1, va 978 979 return value; 980 } 981 } 982 } 983 else if (particleDefinition == G4Proton::Pro 984 { 985 // k should be in eV 986 dataDiffCSMap::iterator iterator_Nrj; 987 iterator_Nrj = pNrjTransStorage.find(cur 988 989 dataProbaShellMap::iterator iterator_Pro 990 iterator_Proba = pProbaShellStorage.find 991 992 incidentEnergyMap::iterator iterator_Tdu 993 iterator_Tdummy = pIncidentEnergyStorage 994 995 if (iterator_Nrj == pNrjTransStorage.end() 996 iterator_Tdummy == pIncidentEnergyStorag 997 { 998 G4String str = "Material "; 999 str += currentMaterial + " not found!"; 1000 G4Exception("G4MicroElecInelasticModel_ 1001 FatalException, str); 1002 } 1003 else 1004 { 1005 vector<TriDimensionMap>* pNrjTransfData = 1006 vector<VecMap>* pProbaShellMap = iterator 1007 vector<G4double>* pTdummyVec = iterator_T 1008 1009 auto k2 = std::upper_bound(pTdummyVec->be 1010 pTdummyVec->end(), 1011 k); 1012 1013 auto k1 = k2 - 1; 1014 1015 // SI : the following condition avoids si 1016 // for eg. when the last element is 1017 if (random <= (*pProbaShellMap)[ionizatio 1018 && random <= (*pProbaShellMap)[ioniza 1019 { 1020 auto prob12 = 1021 std::upper_bound((*pProbaShellMap)[ioniza 1022 (*pProbaShellMap)[ionizationLevelInd 1023 random); 1024 auto prob11 = prob12 - 1; 1025 auto prob22 = 1026 std::upper_bound((*pProbaShellMap)[ioniza 1027 (*pProbaShellMap)[ionizationLevelInd 1028 random); 1029 auto prob21 = prob22 - 1; 1030 1031 valueK1 = *k1; 1032 valueK2 = *k2; 1033 valuePROB21 = *prob21; 1034 valuePROB22 = *prob22; 1035 valuePROB12 = *prob12; 1036 valuePROB11 = *prob11; 1037 1038 // The following condition avoid gett 1039 // and forces cumxs = 1 for maximum e 1040 if (valuePROB11 == 0) nrjTransf11 = b 1041 else nrjTransf11 = (*pNrjTransfData)[ 1042 1043 if (valuePROB12 == 1) nrjTransf12 = m 1044 else nrjTransf12 = (*pNrjTransfData)[ 1045 1046 if (valuePROB21 == 0) nrjTransf21 = b 1047 else nrjTransf21 = (*pNrjTransfData)[ 1048 1049 if (valuePROB22 == 1) nrjTransf22 = m 1050 else nrjTransf22 = (*pNrjTransfData)[ 1051 1052 } 1053 1054 // Avoids cases where cum xs is zero for 1055 if (random > (*pProbaShellMap)[ionization 1056 { 1057 auto prob22 = 1058 std::upper_bound((*pProbaShellMap)[ioniza 1059 (*pProbaShellMap)[ionizationLevelInd 1060 random); 1061 1062 auto prob21 = prob22 - 1; 1063 1064 valueK1 = *k1; 1065 valueK2 = *k2; 1066 valuePROB21 = *prob21; 1067 valuePROB22 = *prob22; 1068 1069 nrjTransf21 = (*pNrjTransfData)[ioniz 1070 nrjTransf22 = (*pNrjTransfData)[ioniz 1071 1072 G4double interpolatedvalue2 = Interpo 1073 valuePROB22, 1074 random, 1075 nrjTransf21, 1076 nrjTransf22); 1077 1078 // zeros are explicitly set 1079 G4double value = Interpolate(valueK1, 1080 return value; 1081 } 1082 } 1083 } 1084 // End electron and proton cases 1085 1086 G4double nrjTransfProduct = nrjTransf11 * n 1087 1088 if (nrjTransfProduct != 0.) 1089 { 1090 nrj = QuadInterpolator(valuePROB11, 1091 valuePROB12, 1092 valuePROB21, 1093 valuePROB22, 1094 nrjTransf11, 1095 nrjTransf12, 1096 nrjTransf21, 1097 nrjTransf22, 1098 valueK1, 1099 valueK2, 1100 k, 1101 random); 1102 } 1103 1104 return nrj; 1105 } 1106 1107 //....oooOO0OOooo........oooOO0OOooo........o 1108 1109 G4double G4MicroElecInelasticModel_new::Diffe 1110 const G4ParticleDefinition * particl 1111 G4double k, 1112 G4double energyTransfer, 1113 G4int LevelIndex) 1114 { 1115 G4double sigma = 0.; 1116 1117 if (energyTransfer >= currentMaterialStruct 1118 { 1119 G4double valueT1 = 0; 1120 G4double valueT2 = 0; 1121 G4double valueE21 = 0; 1122 G4double valueE22 = 0; 1123 G4double valueE12 = 0; 1124 G4double valueE11 = 0; 1125 1126 G4double xs11 = 0; 1127 G4double xs12 = 0; 1128 G4double xs21 = 0; 1129 G4double xs22 = 0; 1130 1131 if (particleDefinition == G4Electron::E 1132 { 1133 1134 dataDiffCSMap::iterator iterator_Proba; 1135 iterator_Proba = eDiffDatatable.find(curr 1136 1137 incidentEnergyMap::iterator iterator_Nrj; 1138 iterator_Nrj = eIncidentEnergyStorage.fin 1139 1140 TranfEnergyMap::iterator iterator_TransfN 1141 iterator_TransfNrj = eVecmStorage.find(cu 1142 1143 if (iterator_Proba != eDiffDatatable.end( 1144 && iterator_TransfNrj!= eVecmStorage. 1145 { 1146 vector<TriDimensionMap>* eDiffCrossSe 1147 vector<G4double>* eTdummyVec = iterat 1148 VecMap* eVecm = iterator_TransfNrj->s 1149 1150 // k should be in eV and energy trans 1151 auto t2 = std::upper_bound(eTdummyVec 1152 auto t1 = t2 - 1; 1153 // SI : the following condition avoid 1154 if (energyTransfer <= (*eVecm)[(*t1)] 1155 { 1156 auto e12 = std::upper_bound((*eVecm)[(* 1157 auto e11 = e12 - 1; 1158 auto e22 = std::upper_bound((*eVecm)[(* 1159 auto e21 = e22 - 1; 1160 1161 valueT1 = *t1; 1162 valueT2 = *t2; 1163 valueE21 = *e21; 1164 valueE22 = *e22; 1165 valueE12 = *e12; 1166 valueE11 = *e11; 1167 1168 xs11 = (*eDiffCrossSectionData)[LevelIn 1169 xs12 = (*eDiffCrossSectionData)[LevelIn 1170 xs21 = (*eDiffCrossSectionData)[LevelIn 1171 xs22 = (*eDiffCrossSectionData)[LevelIn 1172 } 1173 } 1174 else { 1175 G4String str = "Material "; 1176 str += currentMaterial + " not found!"; 1177 G4Exception("G4MicroElecDielectricModel 1178 } 1179 } 1180 1181 if (particleDefinition == G4Proton::Pro 1182 { 1183 dataDiffCSMap::iterator iterator_Proba; 1184 iterator_Proba = pDiffDatatable.find(curr 1185 1186 incidentEnergyMap::iterator iterator_Nrj; 1187 iterator_Nrj = pIncidentEnergyStorage.fin 1188 1189 TranfEnergyMap::iterator iterator_TransfN 1190 iterator_TransfNrj = pVecmStorage.find(cu 1191 1192 if (iterator_Proba != pDiffDatatable.end( 1193 && iterator_TransfNrj != pVecmStorage 1194 { 1195 vector<TriDimensionMap>* pDiffCrossSe 1196 vector<G4double>* pTdummyVec = iterat 1197 VecMap* pVecm = iterator_TransfNrj->s 1198 1199 // k should be in eV and energy trans 1200 auto t2 = 1201 std::upper_bound(pTdummyVec->begin(), pTd 1202 auto t1 = t2 - 1; 1203 if (energyTransfer <= (*pVecm)[(*t1)] 1204 { 1205 auto e12 = std::upper_bound((*pVecm)[(* 1206 auto e11 = e12 - 1; 1207 auto e22 = std::upper_bound((*pVecm)[(* 1208 auto e21 = e22 - 1; 1209 1210 valueT1 = *t1; 1211 valueT2 = *t2; 1212 valueE21 = *e21; 1213 valueE22 = *e22; 1214 valueE12 = *e12; 1215 valueE11 = *e11; 1216 1217 xs11 = (*pDiffCrossSectionData)[LevelIn 1218 xs12 = (*pDiffCrossSectionData)[LevelIn 1219 xs21 = (*pDiffCrossSectionData)[LevelIn 1220 xs22 = (*pDiffCrossSectionData)[LevelIn 1221 } 1222 } 1223 else { 1224 G4String str = "Material "; 1225 str += currentMaterial + " not found!"; 1226 G4Exception("G4MicroElecDielectricModel 1227 } 1228 } 1229 1230 G4double xsProduct = xs11 * xs12 * xs21 1231 if (xsProduct != 0.) 1232 { 1233 sigma = QuadInterpolator( valueE11, v 1234 valueE21, valueE22, 1235 xs11, xs12, 1236 xs21, xs22, 1237 valueT1, valueT2, 1238 k, energyTransfer); 1239 } 1240 1241 } 1242 1243 return sigma; 1244 } 1245 1246 //....oooOO0OOooo........oooOO0OOooo........o 1247 1248 1249 G4double G4MicroElecInelasticModel_new::Inter 1250 G4double e2, 1251 G4double e, 1252 G4double xs1, 1253 G4double xs2) 1254 { 1255 G4double value = 0.; 1256 if (e == e1 || e1 == e2 ) { return xs1; } 1257 if (e == e2) { return xs2; } 1258 1259 // Log-log interpolation by default 1260 if (e1 > 0. && e2 > 0. && xs1 > 0. && xs2 > 1261 { 1262 G4double a = std::log(xs2/xs1)/ std::lo 1263 G4double b = std::log(xs2) - a * std::l 1264 G4double sigma = a * std::log(e) + b; 1265 value = (std::exp(sigma)); 1266 } 1267 1268 // Switch to log-lin interpolation for fast 1269 else if (xs1 > 0. && xs2 > 0. && fasterCode 1270 { 1271 G4double d1 = std::log(xs1); 1272 G4double d2 = std::log(xs2); 1273 value = std::exp((d1 + (d2 - d1) * (e - 1274 } 1275 1276 // Switch to lin-lin interpolation for fast 1277 // in case one of xs1 or xs2 (=cum proba) v 1278 else 1279 { 1280 G4double d1 = xs1; 1281 G4double d2 = xs2; 1282 value = (d1 + (d2 - d1) * (e - e1) / (e 1283 } 1284 1285 return value; 1286 } 1287 1288 //....oooOO0OOooo........oooOO0OOooo........o 1289 1290 G4double G4MicroElecInelasticModel_new::QuadI 1291 G4double e21, G4double e22, 1292 G4double xs11, G4double xs12, 1293 G4double xs21, G4double xs22, 1294 G4double t1, G4double t2, 1295 G4double t, G4double e) 1296 { 1297 G4double interpolatedvalue1 = Interpolate(e 1298 G4double interpolatedvalue2 = Interpolate(e 1299 G4double value = Interpolate(t1, t2, t, int 1300 return value; 1301 } 1302 1303 //....oooOO0OOooo........oooOO0OOooo........o 1304 1305 G4int G4MicroElecInelasticModel_new::RandomSe 1306 { 1307 G4int level = 0; 1308 1309 TCSMap::iterator tablepos; 1310 tablepos = tableTCS.find(currentMaterial); 1311 if (tablepos == tableTCS.end()) { 1312 G4Exception("G4MicroElecInelasticModel_ne 1313 return level; 1314 } 1315 1316 MapData* tableData = tablepos->second; 1317 1318 std::map< G4String,G4MicroElecCrossSectionD 1319 pos = tableData->find(particle); 1320 1321 std::vector<G4double> Zeff(currentMaterialS 1322 if(originalMass>proton_mass_c2) { 1323 for(G4int nl=0;nl<currentMaterialStructur 1324 Zeff[nl] = BKZ(k/(proton_mass_c2/origin 1325 } 1326 } 1327 1328 if (pos != tableData->end()) 1329 { 1330 G4MicroElecCrossSectionDataSet_new* tab 1331 1332 if (table != 0) 1333 { 1334 G4double* valuesBuffer = new G4double[tab 1335 const G4int n = (G4int)table->NumberOfCom 1336 G4int i = (G4int)n; 1337 G4double value = 0.; 1338 1339 while (i>0) 1340 { 1341 --i; 1342 valuesBuffer[i] = table->GetComponent 1343 value += valuesBuffer[i]; 1344 } 1345 value *= G4UniformRand(); 1346 1347 i = n; 1348 1349 while (i > 0) 1350 { 1351 --i; 1352 1353 if (valuesBuffer[i] > value) 1354 { 1355 delete[] valuesBuffer; 1356 return i; 1357 } 1358 value -= valuesBuffer[i]; 1359 } 1360 1361 if (valuesBuffer) delete[] valuesBuffer; 1362 1363 } 1364 } 1365 else 1366 { 1367 G4Exception("G4MicroElecInelasticModel_ 1368 } 1369 1370 return level; 1371 } 1372 1373 //....oooOO0OOooo........oooOO0OOooo........o 1374 1375 G4double G4MicroElecInelasticModel_new::Compu 1376 G4double x = E/mass; 1377 return c_light*std::sqrt(x*(x + 2.0))/(x + 1378 } 1379 1380 //....oooOO0OOooo........oooOO0OOooo........o 1381 1382 G4double G4MicroElecInelasticModel_new::Compu 1383 G4double v1i = ComputeRelativistVelocity(T1 1384 G4double v2i = ComputeRelativistVelocity(T2 1385 1386 G4double v2f = 2*M1/(M1+M2)*v1i + (M2-M1)/( 1387 G4double vtransfer2a = v2f*v2f-v2i*v2i; 1388 1389 v2f = 2*M1/(M1+M2)*v1i + (M2-M1)/(M1+M2)*v2 1390 G4double vtransfer2b = v2f*v2f-v2i*v2i; 1391 1392 G4double vtransfer2 = std::max(vtransfer2a, 1393 return 0.5*M2*vtransfer2; 1394 } 1395 1396 //....oooOO0OOooo........oooOO0OOooo........o 1397 1398 G4double G4MicroElecInelasticModel_new::stepF 1399 return (x < 0.) ? 1.0 : 0.0; 1400 } 1401 1402 //....oooOO0OOooo........oooOO0OOooo........o 1403 1404 G4double G4MicroElecInelasticModel_new::vrkre 1405 { 1406 G4double r = vF*( std::pow(v/vF+1., 3.) - f 1407 + 4.*(v/vF)*(v/vF) ) + stepFunc(v/vF- 1408 4.*(v/vF)*(v/vF) + 3.*std::po 1409 - 0.5*std::pow(v/vF, 5.)); 1410 return r/(10.*v/vF); 1411 } 1412 1413 //....oooOO0OOooo........oooOO0OOooo........o 1414 1415 G4double G4MicroElecInelasticModel_new::BKZ(G 1416 { 1417 // need atomic unit conversion 1418 G4double hbar = hbar_Planck, hbar2 = hbar*h 1419 G4double hartree = 2*Ry, a0 = Bohr_radius, 1420 G4double vp = ComputeRelativistVelocity(Ep, 1421 1422 vp /= velocity; 1423 1424 G4double wp = Eplasmon/hartree; 1425 G4double a = std::pow(4./9./CLHEP::pi, 1./3 1426 G4double vF = std::pow(wp*wp/(3.*a*a*a), 1. 1427 G4double c = 0.9; 1428 G4double vr = vrkreussler(vp /*in u.a*/, vF 1429 G4double yr = vr/std::pow(Zp, 2./3.); 1430 G4double q = 0.; 1431 if(Zp==2) q = 1-exp(-c*vr/(Zp-5./16.)); 1432 else q = 1.-exp(-c*(yr-0.07)); 1433 G4double Neq = Zp*(1.-q); 1434 G4double l0 = 0.; 1435 if(Neq<=2) l0 = 3./(Zp-0.3*(Neq-1))/2.; 1436 else l0 = 0.48*std::pow(Neq, 2./3.)/(Zp-Neq 1437 if(Zp==2) c = 1.0; 1438 else c = 3./2.; 1439 return Zp*(q + c*(1.-q)/vF/vF/2.0 * log(1.+ 1440 } 1441 1442 //....oooOO0OOooo........oooOO0OOooo........o 1443