File indexing completed on 2025-02-23 09:22:04
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042 #include "Run.hh"
0043
0044 #include "DetectorConstruction.hh"
0045 #include "PrimaryGeneratorAction.hh"
0046
0047 #include "G4AnalysisManager.hh"
0048 #include "G4EmCalculator.hh"
0049 #include "G4H2O.hh"
0050 #include "G4Molecule.hh"
0051 #include "G4MoleculeCounter.hh"
0052 #include "G4MoleculeGun.hh"
0053 #include "G4MoleculeTable.hh"
0054 #include "G4SystemOfUnits.hh"
0055 #include "G4UnitsTable.hh"
0056 #include "G4ios.hh"
0057
0058 #include <G4Scheduler.hh>
0059 #include <iomanip>
0060
0061
0062
0063 Run::Run(DetectorConstruction* det) : G4Run(), fDetector(det), fParticle(0)
0064 {
0065 fSoma3DEdep = new G4double[fDetector->GetnbSomacomp()];
0066 for (G4int i = 0; i < fDetector->GetnbSomacomp(); i++) {
0067 fSoma3DEdep[i] = 0.;
0068 }
0069 fDend3DEdep = new G4double[fDetector->GetnbDendritecomp()];
0070 for (G4int i = 0; i < fDetector->GetnbDendritecomp(); i++) {
0071 fDend3DEdep[i] = 0.;
0072 }
0073 fAxon3DEdep = new G4double[fDetector->GetnbAxoncomp()];
0074 for (G4int i = 0; i < fDetector->GetnbAxoncomp(); i++) {
0075 fAxon3DEdep[i] = 0.;
0076 }
0077
0078 fEdepAll = fEdepAll_err = fEdepMedium = fEdepMedium_err = fEdepSlice = fEdepSlice_err =
0079 fEdepSoma = fEdepSoma_err = 0.;
0080 fEdepDend = fEdepDend_err = fEdepAxon = fEdepAxon_err = fEdepNeuron = fEdepNeuron_err = 0.;
0081 }
0082
0083
0084
0085 Run::~Run()
0086 {
0087 delete[] fSoma3DEdep;
0088 delete[] fDend3DEdep;
0089 delete[] fAxon3DEdep;
0090 }
0091
0092
0093
0094 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
0095 {
0096 fParticle = particle;
0097 fEkin = energy;
0098 }
0099
0100
0101
0102 void Run::AddPrimaryLET(G4double let)
0103 {
0104 fLET += let;
0105 fLET2 += let * let;
0106 }
0107
0108
0109
0110 void Run::SetTrackLength(G4double t)
0111 {
0112 ftrackLength = t;
0113 fTrackLen += t;
0114 fTrackLen2 += t * t;
0115 }
0116
0117
0118
0119 void Run::CountProcesses(const G4VProcess* process)
0120 {
0121 G4String procName = process->GetProcessName();
0122 std::map<G4String, G4int>::iterator it = fProcCounter.find(procName);
0123 if (it == fProcCounter.end()) {
0124 fProcCounter[procName] = 1;
0125 }
0126 else {
0127 fProcCounter[procName]++;
0128 }
0129 }
0130
0131
0132
0133 void Run::ParticleCount(G4String name, G4double Ekin)
0134 {
0135 std::map<G4String, ParticleData>::iterator it = fParticleDataMap1.find(name);
0136 if (it == fParticleDataMap1.end()) {
0137 fParticleDataMap1[name] = ParticleData(1, Ekin, Ekin, Ekin);
0138 }
0139 else {
0140 ParticleData& data = it->second;
0141 data.fCount++;
0142 data.fEmean += Ekin;
0143
0144 G4double emin = data.fEmin;
0145 if (Ekin < emin) data.fEmin = Ekin;
0146 G4double emax = data.fEmax;
0147 if (Ekin > emax) data.fEmax = Ekin;
0148 }
0149 }
0150
0151
0152
0153 void Run::ParticleCountNeuron(G4String name, G4double Ekin)
0154 {
0155 std::map<G4String, ParticleData>::iterator it = fParticleDataMap2.find(name);
0156 if (it == fParticleDataMap2.end()) {
0157 fParticleDataMap2[name] = ParticleData(1, Ekin, Ekin, Ekin);
0158 }
0159 else {
0160 ParticleData& data = it->second;
0161 data.fCount++;
0162 data.fEmean += Ekin;
0163
0164 G4double emin = data.fEmin;
0165 if (Ekin < emin) data.fEmin = Ekin;
0166 G4double emax = data.fEmax;
0167 if (Ekin > emax) data.fEmax = Ekin;
0168 }
0169 }
0170
0171
0172
0173
0174
0175
0176
0177
0178
0179
0180
0181 void Run::MoleculeCountNeuron(G4Molecule* molecule)
0182 {
0183 G4String moleculename = molecule->GetName();
0184 std::map<G4String, G4int>::iterator it = fMoleculeNumber.find(moleculename);
0185 if (it == fMoleculeNumber.end()) {
0186 fMoleculeNumber[moleculename] = 1;
0187 }
0188 else {
0189 fMoleculeNumber[moleculename]++;
0190 }
0191 }
0192
0193
0194
0195 void Run::AddEflow(G4double eflow)
0196 {
0197 fEnergyFlow += eflow;
0198 fEnergyFlow2 += eflow * eflow;
0199 }
0200
0201
0202
0203 void Run::Merge(const G4Run* run)
0204 {
0205 const Run* localRun = static_cast<const Run*>(run);
0206
0207
0208
0209 fParticle = localRun->fParticle;
0210 fEkin = localRun->fEkin;
0211 ftrackLength = localRun->ftrackLength;
0212 fTrackLen += localRun->fTrackLen;
0213 fTrackLen2 += localRun->fTrackLen2;
0214 fLET += localRun->fLET;
0215 fLET2 += localRun->fLET2;
0216
0217
0218 std::map<G4String, G4int>::const_iterator itp;
0219 for (itp = localRun->fProcCounter.begin(); itp != localRun->fProcCounter.end(); ++itp) {
0220 G4String procName = itp->first;
0221 G4int localCount = itp->second;
0222 if (fProcCounter.find(procName) == fProcCounter.end()) {
0223 fProcCounter[procName] = localCount;
0224 }
0225 else {
0226 fProcCounter[procName] += localCount;
0227 }
0228 }
0229
0230
0231 std::map<G4String, ParticleData>::const_iterator itc;
0232 for (itc = localRun->fParticleDataMap1.begin(); itc != localRun->fParticleDataMap1.end(); ++itc) {
0233 G4String name = itc->first;
0234 const ParticleData& localData = itc->second;
0235 if (fParticleDataMap1.find(name) == fParticleDataMap1.end()) {
0236 fParticleDataMap1[name] =
0237 ParticleData(localData.fCount, localData.fEmean, localData.fEmin, localData.fEmax);
0238 }
0239 else {
0240 ParticleData& data = fParticleDataMap1[name];
0241 data.fCount += localData.fCount;
0242 data.fEmean += localData.fEmean;
0243 G4double emin = localData.fEmin;
0244 if (emin < data.fEmin) data.fEmin = emin;
0245 G4double emax = localData.fEmax;
0246 if (emax > data.fEmax) data.fEmax = emax;
0247 }
0248 }
0249
0250
0251 std::map<G4String, ParticleData>::const_iterator itn;
0252 for (itn = localRun->fParticleDataMap2.begin(); itn != localRun->fParticleDataMap2.end(); ++itn) {
0253 G4String name = itn->first;
0254 const ParticleData& localData = itn->second;
0255 if (fParticleDataMap2.find(name) == fParticleDataMap2.end()) {
0256 fParticleDataMap2[name] =
0257 ParticleData(localData.fCount, localData.fEmean, localData.fEmin, localData.fEmax);
0258 }
0259 else {
0260 ParticleData& data = fParticleDataMap2[name];
0261 data.fCount += localData.fCount;
0262 data.fEmean += localData.fEmean;
0263 G4double emin = localData.fEmin;
0264 if (emin < data.fEmin) data.fEmin = emin;
0265 G4double emax = localData.fEmax;
0266 if (emax > data.fEmax) data.fEmax = emax;
0267 }
0268 }
0269
0270 std::map<G4String, G4int>::const_iterator itm;
0271 for (itm = localRun->fMoleculeNumber.begin(); itm != localRun->fMoleculeNumber.end(); ++itm) {
0272 G4String moleculeName = itm->first;
0273 G4int localCount = itm->second;
0274 if (fMoleculeNumber.find(moleculeName) == fMoleculeNumber.end()) {
0275 fMoleculeNumber[moleculeName] = localCount;
0276 }
0277 else {
0278 fMoleculeNumber[moleculeName] += localCount;
0279 }
0280 }
0281
0282
0283
0284 for (G4int i = 0; i < fDetector->GetnbSomacomp(); i++) {
0285 fSoma3DEdep[i] += localRun->fSoma3DEdep[i];
0286 }
0287 for (G4int i = 0; i < fDetector->GetnbDendritecomp(); i++) {
0288 fDend3DEdep[i] += localRun->fDend3DEdep[i];
0289 }
0290 for (G4int i = 0; i < fDetector->GetnbAxoncomp(); i++) {
0291 fAxon3DEdep[i] += localRun->fAxon3DEdep[i];
0292 }
0293
0294
0295
0296 fEdepAll += localRun->fEdepAll;
0297 fEdepAll_err += localRun->fEdepAll_err;
0298 fEdepMedium += localRun->fEdepMedium;
0299 fEdepMedium_err += localRun->fEdepMedium_err;
0300 fEdepSlice += localRun->fEdepSlice;
0301 fEdepSlice_err += localRun->fEdepSlice_err;
0302 fEdepSoma += localRun->fEdepSoma;
0303 fEdepSoma_err += localRun->fEdepSoma_err;
0304 fEdepDend += localRun->fEdepDend;
0305 fEdepDend_err += localRun->fEdepDend_err;
0306 fEdepAxon += localRun->fEdepAxon;
0307 fEdepAxon_err += localRun->fEdepAxon_err;
0308 fEdepNeuron += localRun->fEdepNeuron;
0309 fEdepNeuron_err += localRun->fEdepNeuron_err;
0310
0311 fEnergyFlow += localRun->fEnergyFlow;
0312 fEnergyFlow2 += localRun->fEnergyFlow2;
0313
0314 G4Run::Merge(run);
0315 }
0316
0317
0318
0319 void Run::EndOfRun()
0320 {
0321 G4int prec = 5, wid = prec + 2;
0322 G4int dfprec = G4cout.precision(prec);
0323
0324
0325 G4String Particle = fParticle->GetParticleName();
0326 G4double GunArea;
0327
0328 G4Material* material = fDetector->GetTargetMaterial();
0329
0330
0331
0332 fTrackLen /= numberOfEvent;
0333 fTrackLen2 /= numberOfEvent;
0334 G4double rms = fTrackLen2 - fTrackLen * fTrackLen;
0335 if (rms > 0.)
0336 rms = std::sqrt(rms);
0337 else
0338 rms = 0.;
0339
0340 G4int TotNbofEvents = numberOfEvent;
0341 G4double EdepTotal = fEdepAll;
0342 G4double EdepTotal2 = fEdepAll_err;
0343 EdepTotal /= TotNbofEvents;
0344 EdepTotal2 /= TotNbofEvents;
0345 G4double rmst = EdepTotal2 - EdepTotal * EdepTotal;
0346 if (rmst > 0.)
0347 rmst = std::sqrt(rmst);
0348 else
0349 rmst = 0.;
0350
0351
0352 G4EmCalculator emCalculator;
0353 G4double dEdxFull = 0.;
0354 if (fParticle->GetPDGCharge() != 0.) {
0355 dEdxFull = emCalculator.ComputeTotalDEDX(fEkin, fParticle, material);
0356 }
0357
0358
0359 G4double meandEdx = (EdepTotal / fTrackLen) / (keV / um);
0360 G4double meandEdxerr = (rmst / fTrackLen) / (keV / um);
0361 G4int ncols = 0;
0362 G4float tmp, En;
0363 G4int Ntrav = 0;
0364 FILE* fp = fopen("OutputPerEvent.out", "r");
0365 while (1) {
0366 ncols = fscanf(fp, " %f %f %f %f %f %f %f %f", &tmp, &tmp, &tmp, &tmp, &En, &tmp, &tmp, &tmp);
0367 if (ncols < 0) break;
0368 if (En > 0) Ntrav++;
0369 }
0370 fclose(fp);
0371
0372 GunArea = fDetector->GetTotSurfMedium();
0373
0374 G4double FluenceDoseperBeam = 0.160218 * (dEdxFull) / (GunArea * std::pow(10, 18));
0375
0376 G4cout << "\n ======= The summary of simulation results 'neuron' ========\n";
0377 G4cout << "\n Primary particle = " << Particle
0378 << "\n Kinetic energy of beam = " << fEkin / MeV << " A*MeV"
0379 << "\n Particle traversals the neuron = " << Ntrav << " of " << numberOfEvent
0380 << "\n Full LET of beam as formulas = " << dEdxFull / (keV / um) << " keV/um"
0381 << "\n Mean LET of beam as simulation = " << meandEdx << " +- " << meandEdxerr
0382 << " keV/um"
0383 << "\n Mean track length of beam = " << fTrackLen / um << " +- " << rms << " um"
0384 << "\n Particle fluence = " << numberOfEvent / (GunArea / (cm * cm))
0385 << " particles/cm^2"
0386 << "\n Fluence dose (full) = "
0387 << numberOfEvent * FluenceDoseperBeam / (joule / kg) << " Gy"
0388 << "\n Fluence dose ber beam = " << FluenceDoseperBeam / (joule / kg) << " Gy"
0389 << G4endl;
0390
0391 if (numberOfEvent == 0) {
0392 G4cout.precision(dfprec);
0393 return;
0394 }
0395
0396
0397
0398 G4cout << "\n List of generated physical process:" << G4endl;
0399
0400 G4int index = 0;
0401 std::map<G4String, G4int>::iterator it;
0402 for (it = fProcCounter.begin(); it != fProcCounter.end(); it++) {
0403 G4String procName = it->first;
0404 G4int count = it->second;
0405 G4String space = " ";
0406 if (++index % 1 == 0) space = "\n";
0407 G4cout << " " << std::setw(20) << procName << "=" << std::setw(7) << count << space;
0408 }
0409 G4cout << G4endl;
0410
0411
0412
0413 G4cout << "\n List of generated particles outside neuron structure:" << G4endl;
0414
0415 std::map<G4String, ParticleData>::iterator itc;
0416 for (itc = fParticleDataMap1.begin(); itc != fParticleDataMap1.end(); itc++) {
0417 G4String name = itc->first;
0418 ParticleData data = itc->second;
0419 G4int count = data.fCount;
0420 G4double eMean = data.fEmean / count;
0421 G4double eMin = data.fEmin;
0422 G4double eMax = data.fEmax;
0423
0424 G4double Eflow = data.fEmean / TotNbofEvents;
0425
0426 G4cout << " " << std::setw(13) << name << " : " << std::setw(7) << count
0427 << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy") << "\t( "
0428 << G4BestUnit(eMin, "Energy") << " --> " << G4BestUnit(eMax, "Energy")
0429 << ") \tEflow/event = " << G4BestUnit(Eflow, "Energy") << G4endl;
0430 }
0431
0432
0433
0434 G4cout << "\n Number of secondary particles inside neuron structure:" << G4endl;
0435
0436 std::map<G4String, ParticleData>::iterator itn;
0437 for (itn = fParticleDataMap2.begin(); itn != fParticleDataMap2.end(); itn++) {
0438 G4String name = itn->first;
0439 ParticleData data = itn->second;
0440 G4int count = data.fCount;
0441
0442 G4cout << " " << std::setw(13) << name << " : " << std::setw(7) << count << G4endl;
0443 }
0444
0445
0446
0447 G4cout << "\n Number of molecular products inside neuron structure:"
0448 << "\n time: 1 ps - 10 ps " << G4endl;
0449
0450 G4int ind = 0;
0451 std::map<G4String, G4int>::iterator itm;
0452 for (itm = fMoleculeNumber.begin(); itm != fMoleculeNumber.end(); itm++) {
0453 G4String moleculeName = itm->first;
0454 G4int count = itm->second;
0455 G4String space = " ";
0456 if (++ind % 3 == 0) space = "\n";
0457
0458 G4cout << " " << std::setw(13) << moleculeName << " : " << std::setw(7) << count << G4endl;
0459 }
0460
0461
0462 G4cout << "\n Total energy (MeV) deposition :" << G4endl;
0463
0464 G4cout << " " << std::setw(13) << "All volume: " << std::setw(7) << fEdepAll / MeV << "\n "
0465 << " " << std::setw(13) << "Bounding slice: " << std::setw(7)
0466 << (fEdepSlice + fEdepNeuron) / MeV << "\n "
0467 << " " << std::setw(13) << "Neuron: " << std::setw(7) << fEdepNeuron / MeV << "\n "
0468 << " " << std::setw(13) << "Soma: " << std::setw(7) << fEdepSoma / MeV << "\n "
0469 << " " << std::setw(13) << "Dendrites: " << std::setw(7) << fEdepDend / MeV << "\n "
0470 << " " << std::setw(13) << "Axon: " << std::setw(7) << fEdepAxon / MeV << G4endl;
0471
0472
0473
0474
0475 G4int somaCompHit = 0;
0476 G4double somaCompEdep = 0.;
0477 G4double somaCompDose = 0.;
0478 G4double somaCompEdep2 = 0.;
0479 G4double somaCompDose2 = 0.;
0480
0481 remove("Soma3DEdep.out");
0482 for (G4int i = 0; i < fDetector->GetnbSomacomp(); i++) {
0483 if (fSoma3DEdep[i] > 0.0) {
0484 somaCompHit++;
0485 somaCompEdep += fSoma3DEdep[i];
0486 somaCompEdep2 += fSoma3DEdep[i] * fSoma3DEdep[i];
0487
0488 std::ofstream WriteDataInSoma("Soma3DEdep.out", std::ios::app);
0489
0490 WriteDataInSoma << fDetector->GetPosSomacomp(i).x() << '\t' << " "
0491 << fDetector->GetPosSomacomp(i).y() << '\t' << " "
0492 << fDetector->GetPosSomacomp(i).z() << '\t'
0493 << " "
0494
0495 << fSoma3DEdep[i] / keV << '\t' << " " << G4endl;
0496 }
0497 }
0498
0499
0500 G4double rmsEdepS = 0.;
0501 G4double rmsDoseS = 0.;
0502 if (somaCompHit > 0) {
0503 somaCompEdep /= somaCompHit;
0504 somaCompEdep2 /= somaCompHit;
0505 rmsEdepS = somaCompEdep2 - somaCompEdep * somaCompEdep;
0506 if (rmsEdepS > 0.)
0507 rmsEdepS = std::sqrt(rmsEdepS);
0508 else
0509 rmsEdepS = 0.;
0510 somaCompDose /= somaCompHit;
0511 somaCompDose2 /= somaCompHit;
0512 rmsDoseS = somaCompDose2 - somaCompDose * somaCompDose;
0513 if (rmsDoseS > 0.)
0514 rmsDoseS = std::sqrt(rmsDoseS);
0515 else
0516 rmsDoseS = 0.;
0517 }
0518
0519 G4int DendCompHit = 0;
0520 G4double DendCompEdep = 0.;
0521 G4double DendCompDose = 0.;
0522 G4double DendCompEdep2 = 0.;
0523 G4double DendCompDose2 = 0.;
0524 remove("Dend3DEdep.out");
0525 std::ofstream WriteDataInDend("Dend3DEdep.out", std::ios::app);
0526 for (G4int i = 0; i < fDetector->GetnbDendritecomp(); i++) {
0527 if (fDend3DEdep[i] > 0.0) {
0528 DendCompHit++;
0529 DendCompEdep += fDend3DEdep[i];
0530 DendCompEdep2 += fDend3DEdep[i] * fDend3DEdep[i];
0531
0532 WriteDataInDend
0533
0534 << fDetector->GetPosDendcomp(i).x() << '\t' << " " << fDetector->GetPosDendcomp(i).y()
0535 << '\t' << " " << fDetector->GetPosDendcomp(i).z() << '\t' << " "
0536 << fDetector->GetDistADendSoma(i) << '\t' << " " << fDetector->GetDistBDendSoma(i) << '\t'
0537 << " " << fDend3DEdep[i] / keV << '\t' << " " << G4endl;
0538 }
0539 }
0540
0541 G4double rmsEdepD = 0.;
0542 G4double rmsDoseD = 0.;
0543 if (DendCompHit > 0) {
0544 DendCompEdep /= DendCompHit;
0545 DendCompEdep2 /= DendCompHit;
0546 rmsEdepD = DendCompEdep2 - DendCompEdep * DendCompEdep;
0547 if (rmsEdepD > 0.)
0548 rmsEdepD = std::sqrt(rmsEdepD);
0549 else
0550 rmsEdepD = 0.;
0551 DendCompDose /= DendCompHit;
0552 DendCompDose2 /= DendCompHit;
0553 rmsDoseD = DendCompDose2 - DendCompDose * DendCompDose;
0554 if (rmsDoseD > 0.)
0555 rmsDoseD = std::sqrt(rmsDoseD);
0556 else
0557 rmsDoseD = 0.;
0558 }
0559
0560 G4int AxonCompHit = 0;
0561 G4double AxonCompEdep = 0.;
0562 G4double AxonCompDose = 0.;
0563 G4double AxonCompEdep2 = 0.;
0564 G4double AxonCompDose2 = 0.;
0565 remove("Axon3DEdep.out");
0566 std::ofstream WriteDataInAxon("Axon3DEdep.out", std::ios::app);
0567 for (G4int i = 0; i < fDetector->GetnbAxoncomp(); i++) {
0568 if (fAxon3DEdep[i] > 0.0) {
0569 AxonCompHit++;
0570 AxonCompEdep += fAxon3DEdep[i];
0571 AxonCompEdep2 += fAxon3DEdep[i] * fAxon3DEdep[i];
0572
0573 WriteDataInAxon
0574
0575 << fDetector->GetPosAxoncomp(i).x() << '\t' << " " << fDetector->GetPosAxoncomp(i).y()
0576 << '\t' << " " << fDetector->GetPosAxoncomp(i).z() << '\t' << " "
0577 << fDetector->GetDistAxonsoma(i) << '\t' << " " << fAxon3DEdep[i] / keV << '\t' << " "
0578 << G4endl;
0579 }
0580 }
0581
0582 G4double rmsEdepA = 0.;
0583 G4double rmsDoseA = 0.;
0584 if (AxonCompHit > 0) {
0585 AxonCompEdep /= AxonCompHit;
0586 AxonCompEdep2 /= AxonCompHit;
0587 rmsEdepA = AxonCompEdep2 - AxonCompEdep * AxonCompEdep;
0588 if (rmsEdepA > 0.)
0589 rmsEdepA = std::sqrt(rmsEdepA);
0590 else
0591 rmsEdepA = 0.;
0592 AxonCompDose /= AxonCompHit;
0593 AxonCompDose2 /= AxonCompHit;
0594 rmsDoseA = AxonCompDose2 - AxonCompDose * AxonCompDose;
0595 if (rmsDoseA > 0.)
0596 rmsDoseA = std::sqrt(rmsDoseA);
0597 else
0598 rmsDoseA = 0.;
0599 }
0600
0601 G4cout << "\n Number of compartments traversed by particle tracks :" << G4endl;
0602 G4cout << " " << std::setw(13) << "Soma: " << std::setw(7) << somaCompHit
0603 << " of total: " << fDetector->GetnbSomacomp() << "\n "
0604 << " " << std::setw(13) << "Dendrites: " << std::setw(7) << DendCompHit
0605 << " of total: " << fDetector->GetnbDendritecomp() << "\n "
0606 << " " << std::setw(13) << "Axon: " << std::setw(7) << AxonCompHit
0607 << " of total: " << fDetector->GetnbAxoncomp() << "\n " << G4endl;
0608 G4cout << "\n Dendritic (or Axon) compartmental energy deposits \n"
0609 << " at the distance from Soma have been written into *.out files:"
0610 << "\n Dend3DEdep.out, Axon3DEdep.out, Soma3DEdep.out"
0611 << "\n Outputs of energy deposition per event written in data file:"
0612 << "\n OutputPerEvent.out" << G4endl;
0613
0614
0615 fProcCounter.clear();
0616 fParticleDataMap2.clear();
0617
0618
0619 G4cout.precision(dfprec);
0620 }
0621
0622
