MURE
Loading...
Searching...
No Matches
MCNPConnector_detectors_legacy.cxx
Go to the documentation of this file.
1/*
2 This file is part of MURE,
3 Copyright (C) 2007-2021 MURE developers.
4
5 MURE is free software: you can redistribute it and/or modify
6 it under the terms of the GNU Lesser General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
9
10 MURE is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU Lesser General Public License for more details.
14
15 You should have received a copy of the GNU Lesser General Public License
16 along with MURE. If not, see <http://www.gnu.org/licenses/>.
17*/
18
19#include <cstring>
20#include <typeinfo>
21#include <map>
22using namespace std;
23
24#include "MUREGlobal.hxx"
25#include "MiscFunction.hxx"
26#include <libmctal/TMTally.hxx>
27
28#include "MCNPConnector.hxx"
29#include "ReadXSFile.hxx"
30#include "MCNPTally.hxx"
31
32#include "MCNPNode.hxx"
33#include "MCNPCylinder.hxx"
34
35namespace MCNP
36{
37
38//________________________________________________________________________
39void Connector::BuildDetectors()
40{
41 if(fRebuildDetector)
42 {
43 LOG_ERROR("USE OF Connector::BuildDetectors() is deprecated in case of Evolution. PLEASE REMOVE IT FROM YOUR CODE.");
44 }
45 LOG_INFO("---------------------------------------------------------------------------");
46 LOG_INFO("Connector::BuildDetectors : Building tallies ...");
47
48 fRebuildDetector = true;
49
50 gMURE->GetEvolutionSolver()->BuildPerturbativeMaterials();
51
52 map<int, int> TheTemperatureMap;
53 int num = 1;
54 for(int i = 0 ; i < int(gMURE->GetEvolutionSolver()->GetEvolutiveSystemVector().size()); i++)
55 {
56 if(!(gMURE->GetEvolutionSolver()->GetEvolutiveSystemVector()[i]->GetCell()->GetNumber() < 0))
57 {
58 double T = gMURE->GetEvolutionSolver()->GetEvolutiveSystemVector()[i]->GetCell()->GetMaterial()->GetTemperature();
59 int TempIndex = gMURE->GetTemperatureMap()->AddUsedTemp(T);
60 if(TempIndex < 0)
61 {
62 LOG_ERROR("Temperature " << T << " K is not in the fTemperatureMap with the current precision " << gMURE->GetTemperatureMap()->GetDeltaTPrecision() << "K."
63 << "\t This temperature is at " << -TempIndex << " K from the closest fTemperatureMap index.");
64 }
65 if(TheTemperatureMap[TempIndex] == 0)
66 {
67 TheTemperatureMap[TempIndex] = num;
68 num++;
69 }
70 }
71 }
72
73 vector<Tally *> Flux(gMURE->GetTemperatureMap()->GetNumberOfUsedfTemperatures());
74 vector<Tally *> FM(gMURE->GetTemperatureMap()->GetNumberOfUsedfTemperatures()); // sigma*phi tallies
75 for(int i = 0; i < int(Flux.size()); i++)
76 {
77 Flux[i] = new Tally;
78 FM[i] = new Tally;
79 }
80 for(int i = 0 ; i < int(gMURE->GetEvolutionSolver()->GetEvolutiveSystemVector().size()); i++)
81 {
82 Cell *TheCell = gMURE->GetEvolutionSolver()->GetEvolutiveSystemVector()[i]->GetCell();
83 if(!(TheCell->GetNumber() < 0))
84 {
85 double T = TheCell->GetMaterial()->GetTemperature();
86 int TempIndex = gMURE->GetTemperatureMap()->AddUsedTemp(T);
87 int index = TheTemperatureMap[TempIndex] - 1;
88 Flux[index]->Add(TheCell);
89 FM[index]->Add(TheCell);
90 TheCell->SetTallyNumber(Flux[index]->GetNumber());
91 TheCell->SetTallyBinNumber(int(Flux[index]->GetBin().size()) - 1);
92 Material *TheMaterial = TheCell->GetMaterial();
93 LOG_DEBUG("for tally " << Flux[index]->GetNumber() << " bin is " << int(Flux[index]->GetBin().size()) - 1 << " for cell" << TheCell->GetNumber());
94 for(int j = 0; j < int(TheMaterial->GetComposition().size()); j++)
95 {
96 Nucleus_ptr nucleus = gMURE->GetGlobalNucleiVector()[TheMaterial->GetComposition()[j]->GetGlobalNucleiIndex()];
97 int TallyBinNum = int(Flux[index]->GetBin().size()) - 1;
98 //nucleus->AddCellnTallyBinNumber(TheCell->GetNumber(),TallyBinNum);
99 nucleus->AddNucleusMCRecord(FM[index]->GetNumber(), TheCell->GetNumber(), TallyBinNum);
100 }
101 }
102
103 }
104
105 int mbin = 0;
106 int OldTallyNum = 0;
107 for(int i = 0; i < int(gMURE->GetEvolutionSolver()->GetPerturbativeMaterial().size()) ; i++)
108 {
109 if(!(gMURE->GetEvolutionSolver()->GetPerturbativeMaterial()[i]->IsOutCore()))
110 {
111 Nucleus_ptr nucleus = gMURE->GetEvolutionSolver()->GetPerturbativeMaterial()[i]->GetComposition()[0];
112 int TempIndex = gMURE->GetTemperatureMap()->AddUsedTemp(nucleus->GetTemperature());
113 int tallypos = TheTemperatureMap[TempIndex] - 1;
114 int tallynum = FM[tallypos]->GetNumber();
115 for(int j = 0; j < int(nucleus->GetNucleusMCRecord().size()); j++) //it should be only 1!
116 {
117 NucleusMCRecord *TheRecord = nucleus->GetNucleusMCRecord()[j];
118 //nucleus->AddTallyNumber(tallynum);
119 if(tallynum != OldTallyNum)
120 {
121 OldTallyNum = tallynum;
122 mbin = 0;
123 }
124 if (nucleus->GetZAI()->GetReactionDaughterSize() > 0)
125 {
126 vector<int> XScode = TheRecord->AddMultiplicatorBinIndex(mbin, nucleus.Get());
127 for(int k = 0; k < int(XScode.size()); k++)
128 {
129 FM[tallypos]->AddMultiplicator(gMURE->GetEvolutionSolver()->GetPerturbativeMaterial()[i], XScode[k]);
130 }
131 }
132 FM[tallypos]->AddMultiplicator(gMURE->GetEvolutionSolver()->GetPerturbativeMaterial()[i], 1);
133 TheRecord->SetTotalXSMultiplicatorBinIndex(mbin);
134 mbin++;
135 }
136 }
137 }
138 //
139 // Build Global tallies for control material, rods,...
140 //
141 BuildGlobalTallies();
142
143 LOG_INFO("Detectors are built.");
144 LOG_INFO("---------------------------------------------------------------------------");
145
146}
147
148//________________________________________________________________________
149void Connector::UpdateSigmaPhiDetectors()
150{
151 LOG_INFO("---------------------------------------------------------------------------");
152 LOG_INFO("Updating Detectors from " << GetMCDetectorOutputFileName(gMURE->GetRealMCInputFileName()) << " ...");
153
154 ReadDetectorFile();
155
156 if(gMURE->GetSource()->GetKcode())
157 {
158 UpdateKeff();
159 }
160
161 int NumberOfTallies = fMCOutput->GetNbTal();
162
163 //
164 //flux for cells
165 //
166 for(int i = 0 ; i < int(gMURE->GetEvolutionSolver()->GetEvolutiveSystemVector().size()); i++)
167 {
168 Cell *TheCell = gMURE->GetEvolutionSolver()->GetEvolutiveSystemVector()[i]->GetCell();
169 if(!(TheCell->GetNumber() < 0))
170 {
171 int TallyNum = TheCell->GetTallyNumber();
172 int TallyBinNum = TheCell->GetTallyBinNumber();
173 for(int j = 0; j < NumberOfTallies; j++)
174 {
175 TMTally *TheTally = fMCOutput->GetTally(j);
176 //find the the Tally associated to the Cell
177 if(int(TheTally->fNo) == TallyNum)
178 {
179 ValErr_t UnNormalizedPhi = TheTally->fVal[TallyBinNum][0][0][0][0][0][0][0];
180 TheCell->SetFlux(UnNormalizedPhi.fVal);
181 TheCell->SetMCFlux(UnNormalizedPhi);
182 //cout<<"PHHHHI="<<UnNormalizedPhi<<endl;
183 }
184
185 }
186 }
187 }
188 //
189 //sigma*flux for nucleus
190 //
191
192 for(int i = 0 ; i < int(gMURE->GetEvolutionSolver()->GetPerturbativeMaterial().size()); i++)
193 {
194 if(!(gMURE->GetEvolutionSolver()->GetPerturbativeMaterial()[i]->IsOutCore()))
195 {
196 Nucleus_ptr nucleus = gMURE->GetEvolutionSolver()->GetPerturbativeMaterial()[i]->GetComposition()[0];
197 LOG_DEBUG("fPerturbativeMaterial[i]->GetComposition()[0] :" << gMURE->GetEvolutionSolver()->GetPerturbativeMaterial()[i]->GetComposition()[0]->GetZAI()->GetFullName());
198 int CellIdx = 0;
199 for(int RecIdx = 0; RecIdx < int(nucleus->GetNucleusMCRecord().size()) ; RecIdx++)
200 {
201 NucleusMCRecord *Record = nucleus->GetNucleusMCRecord()[RecIdx];
202
203 TMTally *TheTally = 0;
204 for(int j = 0; j < NumberOfTallies; j++)
205 {
206 if(int(fMCOutput->GetTally(j)->fNo) == Record->GetTallyNumber())
207 {
208 TheTally = fMCOutput->GetTally(j);
209 break;
210 }
211 }
212 if(TheTally == 0)
213 {
214 LOG_ERROR("Can't find the tally " << Record->GetTallyNumber());
215 }
216 if (nucleus->GetZAI()->GetReactionDaughterSize() > 0)
217 {
218 for(int j = 0; j < Record->GetNumberOfCellBin(); j++)
219 {
220 int n = Record->GetTallyBinIndex(j);
221 for(int k = 0; k < int(Record->GetNumberOfMultiplicatorBin()); k++)
222 {
223 int m = Record->GetMultiplicatorBinIndex()[k];
224
225 ValErr_t UnNormalizedSigmaPhi = TheTally->fVal[n][0][0][0][m][0][0][0];
226 //if(nucleus->Z()==92 && nucleus->A()==235)
227 //{
228 // cout<<"FM bin"<<m<<" bin="<<n<<" sigfi="<<UnNormalizedSigmaPhi.fVal<<endl;
229
230 //}
231 if(UnNormalizedSigmaPhi.fErr == 0) UnNormalizedSigmaPhi.fErr = 1.e-10;
232
233 if(gMURE->GetEvolutionSolver()->IsPartialMCRun() && gMURE->GetEvolutionSolver()->IsPartialParticleMCRun())
234 {
235 nucleus->SetMCSigmaPhi2(CellIdx, k, UnNormalizedSigmaPhi);
236 }
237 else
238 {
239 nucleus->SetMCSigmaPhi(CellIdx, k, UnNormalizedSigmaPhi);
240 }
241 }
242
243 int m = Record->GetTotalXSMultiplicatorBinIndex();
244 double UnNormalizedTotalSigmaPhi = TheTally->fVal[n][0][0][0][m][0][0][0].fVal;
245 /* if(nucleus->Z()==92 && nucleus->A()==235)
246 {
247 cout<<"for tall ="<<TheTally->fNo<<" cell idx="<<CellIdx<<endl;
248 cout<<" for nuc "<<nucleus->Print()<<" setting TOTAL XS("<<UnNormalizedTotalSigmaPhi<<") from cell bin "<<CellIdx<<"("<<n<<") for fm bin ("<<m<<")"<<endl;
249 } */
250 nucleus->SetTotalSigmaPhi(CellIdx, UnNormalizedTotalSigmaPhi);
251 CellIdx++;
252 }
253 }
254 }
255 }
256 }
257 //
258 // Reading in the Global tally results for every 'real cell'
259 //
260 UpdateGlobalTallies();
261
262 //
263 // The file is read and every thing is updated : delete MCNP reader
264 //
265 delete fMCOutput;
266
267 LOG_INFO("All Detectors read and updated.");
268 LOG_INFO("---------------------------------------------------------------------------");
269}
270
271//________________________________________________________________________
272void Connector::BuildMultiGroupDetectors()
273{
274 LOG_INFO("---------------------------------------------------------------------------");
275 LOG_INFO("Connector::BuildMultiGroupDetectors: Building tallies ...");
276
277 gMURE->GetEvolutionSolver()->BuildPerturbativeMaterials();
278 vector<Material *> MaterialVector = gMURE->GetEvolutionSolver()->GetPerturbativeMaterial();
279 BuildMultiGroupDetectors(MaterialVector);
280
281 LOG_INFO("Multigroup Detectors are built.");
282 LOG_INFO("---------------------------------------------------------------------------");
283}
284
285//________________________________________________________________________
286void Connector::BuildMultiGroupDetectors(vector<Material *> &MaterialVector)
287{
288
289 if(fRebuildDetector)
290 {
291 LOG_ERROR("USE OF Connector::BuildDetectors is deprecated in case of Evolution. PLEASE REMOVE IT FROM YOUR CODE.");
292 }
293
294 fRebuildDetector = true;
295
296
297 map<int, int> TheTemperatureMap;
298 int num = 1;
299 for(int i = 0 ; i < int(gMURE->GetEvolutionSolver()->GetEvolutiveSystemVector().size()); i++)
300 {
301 if(!(gMURE->GetEvolutionSolver()->GetEvolutiveSystemVector()[i]->GetCell()->GetNumber() < 0))
302 {
303 double T = gMURE->GetEvolutionSolver()->GetEvolutiveSystemVector()[i]->GetCell()->GetMaterial()->GetTemperature();
304 int TempIndex = gMURE->GetTemperatureMap()->AddUsedTemp(T);
305 if(TempIndex < 0)
306 {
307 LOG_ERROR("Temperature " << T << " K is not in the fTemperatureMap with the current precision " << gMURE->GetTemperatureMap()->GetDeltaTPrecision() << "K."
308 << "\t This temperature is at " << -TempIndex << " K from the closest fTemperatureMap index.");
309 }
310 if(TheTemperatureMap[TempIndex] == 0)
311 {
312 TheTemperatureMap[TempIndex] = num;
313 num++;
314 }
315 }
316 }
317
318 vector<Tally *> Flux(1);
319 vector<Tally *> FM(1);
320 vector<Tally *> FM8(1);
321
322 Flux[0] = new Tally;
323 FM[0] = new Tally;
324 FM[0]->SetPrintable(false);
325
326 if(gMURE->GetEvolutionSolver()->IsU8StdTally()) //std tallies are used for U-238 only
327 {
328 FM8.resize(gMURE->GetTemperatureMap()->GetNumberOfUsedfTemperatures());
329 for(int i = 0; i < int(FM8.size()); i++)
330 FM8[i] = new Tally;
331 }
332 //
333 // Here we define the groups...
334 // It can be wrong for cross-sections >20MeV or very thermal spectra
335 //
336 Flux[0]->AddMultiGroupEnergy(gMURE->GetEvolutionSolver()->GetMultiGroupLowerEnergy(), gMURE->GetEvolutionSolver()->GetMultiGroupUpperEnergy(), gMURE->GetEvolutionSolver()->GetMultiGroupDecadeMultiplicator());
337
338 for(int i = 0 ; i < int(gMURE->GetEvolutionSolver()->GetEvolutiveSystemVector().size()); i++)
339 {
340 Cell *TheCell = gMURE->GetEvolutionSolver()->GetEvolutiveSystemVector()[i]->GetCell();
341 if(!(TheCell->GetNumber() < 0))
342 {
343 double T = TheCell->GetMaterial()->GetTemperature();
344 int TempIndex = gMURE->GetTemperatureMap()->AddUsedTemp(T);
345 Flux[0]->Add(TheCell);
346 FM[0]->Add(TheCell, false);
347 int index = TheTemperatureMap[TempIndex] - 1;
348 if(gMURE->GetEvolutionSolver()->IsU8StdTally())
349 {
350 FM8[index]->Add(TheCell);
351 }
352 TheCell->SetTallyNumber(Flux[0]->GetNumber());
353 TheCell->SetTallyBinNumber(int(Flux[0]->GetBin().size()) - 1);
354 Material *TheMaterial = TheCell->GetMaterial();
355 LOG_DEBUG("for tally " << Flux[0]->GetNumber() << " bin is " << int(Flux[0]->GetBin().size()) - 1 << " for cell" << TheCell->GetNumber());
356 for(int j = 0; j < int(TheMaterial->GetComposition().size()); j++)
357 {
358 Nucleus_ptr nucleus = gMURE->GetGlobalNucleiVector()[TheMaterial->GetComposition()[j]->GetGlobalNucleiIndex()];
359 int TallyBinNum = int(Flux[0]->GetBin().size()) - 1;
360 int TallyNum = FM[0]->GetNumber();
361 if(gMURE->GetEvolutionSolver()->IsU8StdTally() && nucleus->Z() == 92 && nucleus->A() == 238)
362 {
363 TallyBinNum = int(FM8[index]->GetBin().size()) - 1;
364 TallyNum = FM8[index]->GetNumber();
365 }
366 //if(gMURE->GetEvolutionSolver()->IsU8StdTally() && nucleus->Z()==92 && nucleus->A()==238)
367 // nucleus->SetBinNumber(int(FM8[TheTemperatureMap[TempIndex]-1]->GetBin().size())-1);//here the bin number is the FM bin num
368 //else
369 // nucleus->SetBinNumber(int(Flux[0]->GetBin().size())-1);
370 //nucleus->AddCellNumber(TheCell->GetNumber());
371 //nucleus->AddCellnTallyBinNumber(TheCell->GetNumber(),TallyBinNum);
372
373 nucleus->AddNucleusMCRecord(TallyNum, TheCell->GetNumber(), TallyBinNum);
374 /* if(nucleus->Z()==92 && (nucleus->A()==235 || nucleus->A()==238) )
375 {
376 cout<<"adding cell "<<TheCell->GetNumber()<<" and bin "<<TallyBinNum<<" to "<<nucleus->Print()<<endl;
377 } */
378 }
379 }
380
381 }
382 int mbin = 0;
383 int OldTallyNum = 0;
384 int mbin8 = 0;
385 int OldTallyNum8 = 0;
386 for(int i = 0; i < int(MaterialVector.size()) ; i++)
387 {
388 if(!(MaterialVector[i]->IsOutCore()))
389 {
390 for(int kkk = 0; kkk < int(MaterialVector[i]->GetComposition().size()); kkk++)
391 {
392 Nucleus_ptr nucleus = MaterialVector[i]->GetComposition()[kkk];
393 int TempIndex = gMURE->GetTemperatureMap()->AddUsedTemp(nucleus->GetTemperature());
394 int tallypos = 0;
395 int tallynum = 0;
396 for(int j = 0; j < int(nucleus->GetNucleusMCRecord().size()); j++)
397 {
398 NucleusMCRecord *TheRecord = nucleus->GetNucleusMCRecord()[j];
399 if(gMURE->GetEvolutionSolver()->IsU8StdTally() && nucleus->Z() == 92 && nucleus->A() == 238) //std tallies are used for U-238 only
400 {
401 tallypos = TheTemperatureMap[TempIndex] - 1;
402 tallynum = FM8[tallypos]->GetNumber();
403 //nucleus->AddTallyNumber(tallynum);
404 /* if(nucleus->Z()==92 && nucleus->A()==238)
405 cout<<"adding tally "<<tallynum<<" to "<<nucleus->Print()<<" in bin "<< mbin8<<endl; */
406 if(tallynum != OldTallyNum8)
407 {
408 OldTallyNum8 = tallynum;
409 mbin8 = 0;
410 }
411 }
412 else //not U-238 or full multi-group
413 {
414 tallynum = FM[tallypos]->GetNumber();
415 //nucleus->AddTallyNumber(tallynum);
416 /* if(nucleus->Z()==92 && nucleus->A()==235)
417 cout<<"adding tally "<<tallynum<<" to "<<nucleus->Print()<<" in bin "<< mbin<<endl;
418 */
419 if(tallynum != OldTallyNum)
420 {
421 OldTallyNum = tallynum;
422 mbin = 0;
423 }
424 }
425 if (nucleus->GetZAI()->GetReactionDaughterSize() > 0)
426 {
427 int FMbin = mbin;
428 if(gMURE->GetEvolutionSolver()->IsU8StdTally() && nucleus->Z() == 92 && nucleus->A() == 238) //std tallies are used for U-238 only
429 FMbin = mbin8;
430
431 //here AddFMBinNumber modify the FMbin var ; thus need to affect either mbin or mbin8
432 vector<int> XScode = TheRecord->AddMultiplicatorBinIndex(FMbin, nucleus.Get());
433
434 if(gMURE->GetEvolutionSolver()->IsU8StdTally() && nucleus->Z() == 92 && nucleus->A() == 238) //std tallies are used for U-238 only
435 mbin8 = FMbin;
436 else//not U-238 or full multi-group
437 mbin = FMbin;
438
439 for(int k = 0; k < int(XScode.size()); k++)
440 {
441 if(gMURE->GetEvolutionSolver()->IsU8StdTally() && nucleus->Z() == 92 && nucleus->A() == 238) //std tallies are used for U-238 only
442 {
443 FM8[tallypos]->AddMultiplicator(MaterialVector[i], XScode[k]);
444 /* if(XScode[k]==-6)
445 cout<< nucleus->Print()<<" adding -6 to bin "<<FMbin<<" in tal="<<FM8[tallypos]->GetNumber()<<endl;
446 */
447 }
448 else //not U-238 or full multi-group
449 {
450 if(XScode[k] == -6 && nucleus->Z() == 92 && nucleus->A() == 235)
451 /*cout<< nucleus->Print()<<" adding -6 to bin "<<FMbin<<" in tal="<<FM[tallypos]->GetNumber()<<endl;
452 */ FM[tallypos]->AddMultiplicator(MaterialVector[i], XScode[k]);
453 }
454 }
455 nucleus->SetXSCode(XScode);
456 }
457 if(gMURE->GetEvolutionSolver()->IsU8StdTally() && nucleus->Z() == 92 && nucleus->A() == 238) //std tallies are used for U-238 only
458 {
459 FM8[tallypos]->AddMultiplicator(MaterialVector[i], 1);
460
461 TheRecord->SetTotalXSMultiplicatorBinIndex(mbin8);
462 //cout<< nucleus->Print()<<" adding totalFM bin "<<FM8[tallypos]->GetNumber()<<" in bin "<<mbin8<<endl;
463 mbin8++;
464 }
465 else //not U-238 or full multi-group
466 {
467 FM[tallypos]->AddMultiplicator(MaterialVector[i], 1);
468 TheRecord->SetTotalXSMultiplicatorBinIndex(mbin);
469 //if(nucleus->Z()==92 && nucleus->A()==235)
470 // cout<<"adding total FM bin "<<FM[tallypos]->GetNumber()<<" to "<<nucleus->Print()<<" in bin "<< mbin<<endl;
471 mbin++;
472 }
473 }
474 }
475 }
476 }
477 //
478 // Build Global tallies for control material, rods,...
479 //
480 BuildGlobalTallies();
481
482
483}
484//________________________________________________________________________
485void Connector::ReadDetectorFile()
486{
487
488 string Mfile = gMURE->GetMCRunDirectory() + "/" + GetMCDetectorOutputFileName(gMURE->GetRealMCInputFileName());
489 ifstream MyFile(Mfile.c_str());
490
491 fMCOutput = new TMctal;
492
493
494 if(!fMCOutput->Read(MyFile, false))
495 {
496 LOG_ERROR("Error while reading tallies from m file : could not open file " << Mfile);
497 }
498
499}
500
501//________________________________________________________________________
502void Connector::UpdateMultiGroupSigmaPhiDetectors()
503{
504 LOG_INFO("---------------------------------------------------------------------------");
505 LOG_INFO("Updating MultiGroup detectors from " << GetMCDetectorOutputFileName(gMURE->GetRealMCInputFileName()) << " ...");
506
507 vector<Material *> MaterialVector = gMURE->GetEvolutionSolver()->GetPerturbativeMaterial();
508 UpdateMultiGroupSigmaPhiDetectors(MaterialVector);
509
510 LOG_INFO("All Multigroup Detector read and updated.");
511 LOG_INFO("---------------------------------------------------------------------------");
512}
513
514//________________________________________________________________________
515void Connector::UpdateMultiGroupSigmaPhiDetectors(vector<Material *> &MaterialVector)
516{
517
518 ReadDetectorFile();
519
520 if(gMURE->GetSource()->GetKcode())
521 {
522 UpdateKeff();
523 }
524
525 int NumberOfTallies = fMCOutput->GetNbTal();
526
527 //
528 //flux for cells
529 //
530 for(int i = 0 ; i < int(gMURE->GetEvolutionSolver()->GetEvolutiveSystemVector().size()); i++)
531 {
532 Cell *TheCell = gMURE->GetEvolutionSolver()->GetEvolutiveSystemVector()[i]->GetCell();
533 if(!(TheCell->GetNumber() < 0))
534 {
535 int TallyNum = TheCell->GetTallyNumber();
536 int TallyBinNum = TheCell->GetTallyBinNumber();
537 //if(TheCell->GetNumber()==11)
538 //cout<<"Reading Flux for cell "<<TheCell->GetNumber()<<" in tall="<<TallyNum<<" from bin "<<TallyBinNum<<endl;
539 for(int j = 0; j < NumberOfTallies; j++)
540 {
541 TMTally *TheTally = fMCOutput->GetTally(j);
542 //find the the Tally associated to the Cell
543 if(int(TheTally->fNo) == TallyNum)
544 {
545 float *EnergyBins = TheTally->fEV;
546 vector<float> NRJ;
547 NRJ.assign(EnergyBins, EnergyBins + (TheTally->fEN - 1));
548 TheCell->SetEnergyGroup(NRJ);
549 TheCell->BuildMultiGroupFlux();
550 //double Correlation=0;
551 //double Sum=0;
552 for(int e = 0; e < int(TheTally->fEN) - 1; e++)
553 {
554 ValErr_t UnNormalizedPhi = TheTally->fVal[TallyBinNum][0][0][0][0][0][e][0];
555 TheCell->SetMultiGroupFlux(e, UnNormalizedPhi);
556 //Correlation+=UnNormalizedPhi.fErr*UnNormalizedPhi.fErr;
557 //Sum+=UnNormalizedPhi.fVal;
558 }
559
560 ValErr_t Phi = TheTally->fVal[TallyBinNum][0][0][0][0][0][int(TheTally->fEN) - 1][0];
561 //Correlation=Phi.fErr/sqrt(Correlation);
562 //cout<<"cell="<<TheCell->GetNumber()<<" flux="<<Phi<<" sum="<<Sum<<endl;
563 TheCell->SetFlux(Phi.fVal);
564 TheCell->SetMCFlux(Phi);
565 }
566 }
567 }
568 }
569 //
570 // sigma*flux for nucleus
571 //
572 int kkk = 0;
573 char Wheele[8] = {'|', '/', '-', '\\', '|', '/', '-', '\\'};
574 cout << "Reading Sigma from ACE file. Please Wait.... " << flush;
575 #pragma omp parallel for reduction(+:kkk) schedule(static)
576 for(int i = 0 ; i < int(MaterialVector.size()); i++)
577 {
578 bool IsEvolvingOrPerturbative = (MaterialVector[i]->IsEvolving() || MaterialVector[i]->IsPerturbative());
579 if(!(MaterialVector[i]->IsOutCore()) && IsEvolvingOrPerturbative )
580 {
581 //int jjj=0;
582 //cout<<" maat num="<<MaterialVector[i]->GetNumber()<<" size="<<MaterialVector[i]->GetComposition().size()<<" ctr="<<MaterialVector[i]->IsControlMaterial()<<endl;
583 for(int jjj = 0; jjj < int(MaterialVector[i]->GetComposition().size()); jjj++)
584 {
585 if(MaterialVector[i]->GetComposition(jjj)->GetZAI()->IsInDataBase())
586 {
587 //
588 // just to see a rotating char to wait without panic
589 //
590 printf("%c%c", char(8), Wheele[kkk]);
591 fflush(stdout);
592 kkk = (kkk + 1) % 8;
593
594 //
595 // Now the job
596 //
597 Nucleus_ptr nucleus = MaterialVector[i]->GetComposition()[jjj];
598
599 //LOG_DEBUG("MaterialVector[i]->GetComposition()[0]"<<MaterialVector[i]->GetComposition()[jjj]->GetZAI()->GetFullName());
600 //cout<<"MaterialVector[i]->GetComposition()[0]="<<MaterialVector[i]->GetComposition()[jjj]->GetZAI()->GetFullName()<<endl;
601
602 // OOOOOOOOOOOO Starting of std tallies treatment for U-238 oooooooooooooooo
603 if(gMURE->GetEvolutionSolver()->IsU8StdTally() && nucleus->Z() == 92 && nucleus->A() == 238) //std MCNP tallies for U-238 only
604 {
605 int CellIdx = 0;
606 for(int RecIdx = 0; RecIdx < int(nucleus->GetNucleusMCRecord().size()); RecIdx++)
607 {
608 NucleusMCRecord *Record = nucleus->GetNucleusMCRecord()[RecIdx];
609
610 TMTally *TheTally = 0;
611 for(int j = 0; j < NumberOfTallies; j++)
612 {
613 if(int(fMCOutput->GetTally(j)->fNo) == Record->GetTallyNumber())
614 {
615 TheTally = fMCOutput->GetTally(j);
616 break;
617 }
618 }
619 //cout<<endl<<fPerturbativeMaterial[i]->PrintAll()<<endl;
620
621 if(TheTally == 0)
622 {
623 LOG_ERROR("Can't find the tally " << Record->GetTallyNumber());
624 }
625 if (nucleus->GetZAI()->GetReactionDaughterSize() > 0)
626 {
627 for(int j = 0 ; j < Record->GetNumberOfCellBin() ; j++)
628 {
629 int n = Record->GetTallyBinIndex(j);
630 for(int k = 0; k < Record->GetNumberOfMultiplicatorBin(); k++)
631 {
632 int m = Record->GetMultiplicatorBinIndex()[k];
633
634 ValErr_t UnNormalizedSigmaPhi = TheTally->fVal[n][0][0][0][m][0][0][0];
635 if(UnNormalizedSigmaPhi.fErr == 0)
636 UnNormalizedSigmaPhi.fErr = 1.e-10;
637 if(gMURE->GetEvolutionSolver()->IsPartialMCRun() && gMURE->GetEvolutionSolver()->IsPartialParticleMCRun())
638 {
639 nucleus->SetMCSigmaPhi2(CellIdx, k, UnNormalizedSigmaPhi);
640 }
641 else
642 nucleus->SetMCSigmaPhi(CellIdx, k, UnNormalizedSigmaPhi);
643 }
644 int m = Record->GetTotalXSMultiplicatorBinIndex();
645 double UnNormalizedTotalSigmaPhi = TheTally->fVal[n][0][0][0][m][0][0][0].fVal;
646 nucleus->SetTotalSigmaPhi(CellIdx, UnNormalizedTotalSigmaPhi);
647 CellIdx++;
648 }
649 }
650 }
651 // OOOOOOOOOOOO end of std tallies treatment for U-238 oooooooooooooooo
652 }
653 else //not U-238 or full multi-group
654 {
655 XSDIRLine *good_line = new XSDIRLine(nucleus->GetBestLine());
656 // Test on directory name of temporary XS files (if this directory == _xs_tmp pearhaps the adress directory is not good)
657 ifstream thebase;
658 ifstream temp((good_line->GetXSFileName()).c_str());
659 if(!temp.good())
660 {
661 stringstream KeyWord;
662 string XSTMPDirectory = good_line->GetXSFileName();
663 KeyWord << XSTMPDirectory[0] << XSTMPDirectory[1] << XSTMPDirectory[2] << XSTMPDirectory[3] << XSTMPDirectory[4] << XSTMPDirectory[5] << XSTMPDirectory[6];
664 string TestWord = "_xs_tmp";
665 if(!strcmp(TestWord.c_str(), KeyWord.str().c_str())) // CAREFUL if string are identical : return 0
666 {
667 KeyWord.str("");
668 KeyWord << gMURE->GetMCRunDirectory() << "/" << good_line->GetXSFileName();
669 thebase.open(KeyWord.str().c_str());
670 }
671 else
672 LOG_ERROR("Can't find " << good_line->GetXSFileName() << " file");
673 }
674 else
675 thebase.open((good_line->GetXSFileName()).c_str());
676 temp.close();
677
678 int RCL = good_line->GetRecordLength();
679#ifdef MCNP5
680 if(RCL == 1024)
681 RCL *= 4; // HERE because with fortran default (intel, ...) rcl are in word and not in bytes
682#endif
683 ReadXSFile *xsf;
684 if(good_line->GetDataType() == 2) // we are reading a Binary library
685 xsf = new ReadXSFile(&thebase, good_line->GetStartRecord(), good_line->GetTableLength(),
686 RCL, good_line->GetEntriesPerRecord(), good_line->GetXSFileName());
687 else //it is a ASCII library
688 xsf = new ReadXSFile(&thebase, good_line->GetStartRecord(), good_line->GetTableLength(), good_line->GetXSFileName());
689
690 if (nucleus->GetZAI()->GetReactionDaughterSize() > 0)
691 {
692 //cout<<nucleus->GetBestLine()<<endl;
693 int CellIdx = 0;
694
695 for(int RecIdx = 0; RecIdx < int(nucleus->GetNucleusMCRecord().size()); RecIdx++)
696 {
697 NucleusMCRecord *Record = nucleus->GetNucleusMCRecord()[RecIdx];
698 for(int j = 0; j < Record->GetNumberOfCellBin(); j++)
699 {
700 Cell *TheCell = gMURE->FindCell(Record->GetCellNumber(j));
701
702 for(int ii = 0; ii < int(nucleus->GetXSCode().size()); ii++)
703 {
704 int ReactionCode = nucleus->GetXSCode()[ii];
705 int k = nucleus-> FindFMBin(ReactionCode);
706 //int n = nucleus->GetCellnTallyBinNumber()[j].second; // ?????? //nucleus->GetBinNumber()[j];
707 //cout<<"For cell num="<<TheCell->GetNumber()<<" NG="<<TheCell->GetNEnergyGroup()<<" MultiGroupFlux="<<TheCell->GetMultiGroupFlux().size()<<" "<<TheCell->GetEnergyGroups()<<endl;
708 //cout<< "For cell, tally num="<<TheCell->GetTallyNumber()<<endl;
709 ValErr_t UnNormalizedSigmaPhi = xsf->SigmaPhi(TheCell->GetNEnergyGroup(), TheCell->GetMultiGroupFlux(),
710 TheCell->GetEnergyGroups().data(), ReactionCode);
711
712 /*if(ReactionCode==-6)
713 {
714
715 if(nucleus->Z()==92 && nucleus->A()==235)
716 {
717 cout<<nucleus->Print()<<" fission is set to j="<<j<<" k="<<k<<" of cell="<<Record->GetCellNumber(j)<<" bin="<<Record->GetTallyBinIndex(j)<<endl;
718 cout<<"XS="<<UnNormalizedSigmaPhi<<endl;
719 }
720 }*/
721
722 if(UnNormalizedSigmaPhi.fErr == 0)
723 UnNormalizedSigmaPhi.fErr = 1.e-10;
724 if(gMURE->GetEvolutionSolver()->IsPartialMCRun() && gMURE->GetEvolutionSolver()->IsPartialParticleMCRun())
725 {
726 nucleus->SetMCSigmaPhi2(CellIdx, k, UnNormalizedSigmaPhi);
727 }
728 else
729 nucleus->SetMCSigmaPhi(CellIdx, k, UnNormalizedSigmaPhi);
730 }
731 if(gMURE->IsBetaCalculation()) // Test if Nucleus is used in Beta calculation
732 if(gMURE->IsNucleusUsedForBeta(nucleus->GetZAI()->Z(), nucleus->GetZAI()->A()))
733 {
734 ValErr_t UnNormalizedNuSigmaFisPhi = xsf->NuSigmaFisPhi(TheCell->GetNEnergyGroup(), TheCell->GetMultiGroupFlux(),
735 TheCell->GetEnergyGroups().data());
736 nucleus->SetNuSigmaFisPhi(CellIdx, UnNormalizedNuSigmaFisPhi.fVal); // sets nu tot sigma phi for an energy bin
737 }
738
739 double UnNormalizedTotalSigmaPhi = xsf->SigmaPhi(TheCell->GetNEnergyGroup(), TheCell->GetMultiGroupFlux(),
740 TheCell->GetEnergyGroups().data(), 1).fVal;
741 nucleus->SetTotalSigmaPhi(CellIdx, UnNormalizedTotalSigmaPhi);
742 /* if(nucleus->Z()==92 && nucleus->A()==235)
743 {
744 cout<<nucleus->Print()<<" Total bin is set to CellIdx="<<CellIdx<<" of cell="<<Record->GetCellNumber(j)<<endl;
745 cout<<"XS="<<UnNormalizedTotalSigmaPhi<<endl;
746 } */
747 CellIdx++;
748 }
749 }
750 }
751 delete xsf;
752 }
753 } //end of "if is in data base"
754 }//end of "for on material compo"
755 }//end of "if !IsOutCore"
756 }//end of "for Material"
757
758 //
759 // Reading in the Global tally results for every 'true cell'
760 //
761 //
762 UpdateGlobalTallies();
763
764 //
765 // The file is read and every thing is updated : delete MCNP reader
766 //
767 delete fMCOutput;
768
769}
770
771
772//________________________________________________________________________
773void Connector::UpdateKeff()
774{
775 bool IsMCOutfileRead = true;
776 if(fMCOutput == 0)
777 {
778 IsMCOutfileRead = false;
779 ReadDetectorFile();
780 }
781
782 TMKcode *TheKcode = fMCOutput->GetKcode();
783 unsigned long Taille = TheKcode->GetTaille();
784 unsigned long LastCycle = TheKcode->GetTotalCycle() - 1;
785
786 if(Taille == 19) gMURE->SetKeff(TheKcode->GetCycle19(LastCycle)->fKEffAll.fVal);
787 if(Taille == 19) gMURE->SetKeff_Err(TheKcode->GetCycle19(LastCycle)->fKEffAll.fErr);
788
789 if(!IsMCOutfileRead) //this is just an keff update; one as to delete TMctal object
790 {
791 delete fMCOutput;
792 }
793
794}
795
796//________________________________________________________________________
797void Connector::BuildGlobalTallies()
798{
799 //
800 //Tallies for all real cells(real=non void cell really present in the geometry)
801 //
802 for(int i = 0; i < int(gMURE->GetCellVector().size()); i++)
803 {
804 if(gMURE->GetCellVector()[i]->GetMaterial() && !(gMURE->GetCellVector()[i]->GetNumber() < 0) )
805 //No Tallies build on virtual cells linked to OutCore Materials!
806 {
807 gMURE->GetTrueCellVector().push_back(gMURE->GetCellVector()[i]);
808 gMURE->GetCellVector()[i]->SetTrueCell();
809 }
810 }
811 //out<<"in global"<<endl;
812 //
813 //build global tallies(fission,absorption,...)if needed (i.e.for evolution control)
814 //
815 // MaterialControl (Heavy Nucleus, Boron, ...)
816 if(gMURE->GetEvolutionSolver()->IsMaterialControl())
817 {
818 if(gMURE->GetTrueCellVector().size())
819 {
820 Reaction *Nufiss = new Reaction(-6);
821 Nufiss->Multiply(-7);
822
823 vector<Tally *> GlobalTally(gMURE->GetTrueCellVector().size());
824 for(int i = 0; i < int(GlobalTally.size()); i++)
825 {
826 GlobalTally[i] = new Tally;
827 }
828 for(int i = 0; i < int(gMURE->GetTrueCellVector().size()); i++)
829 {
830 GlobalTally[i]->Add(gMURE->GetTrueCellVector()[i], false);
831 //cout<<i<<" cell "<<gMURE->GetTrueCellVector()[i]->GetNumber()<<endl;
832 gMURE->GetTrueCellVector()[i]->SetGlobalTallyNumber(GlobalTally[i]->GetNumber());
833 Material *Mat = gMURE->GetTrueCellVector()[i]->GetMaterial();
834 GlobalTally[i]->SetBinVolume(0, 1e-6); //SD card must be=1!
835 GlobalTally[i]->AddMultiplicator(Mat, -2, -1); //Order of this multiplicator is not random!
836 GlobalTally[i]->AddMultiplicator(Mat, 16, -1); //0=abs,1=n2n,2=n3n,4=fiss,5=Nu*fiss
837 GlobalTally[i]->AddMultiplicator(Mat, 17, -1);
838 GlobalTally[i]->AddMultiplicator(Mat, -6, -1);
839 GlobalTally[i]->AddMultiplicator(Mat, Nufiss, -1);
840
841 if(dynamic_cast<ControlMaterial *>(Mat))
842 {
843 //First find the perturbative material poison
844 Material *TheControlMat = gMURE->GetEvolutionSolver()->FindPerturbativeMaterial(Mat);
845 if(!TheControlMat)
846 LOG_ERROR("No perturbative ControlMaterial found.");
847 //Here are the global(over volume)reaction rates of ControlMaterial (Poison,Fissile,...)
848 for(int j = 0; j < int(((ControlMaterial *)Mat)->GetControlReactions().size()); j++)
849 {
850 GlobalTally[i]->AddMultiplicator(TheControlMat, &((ControlMaterial *)Mat)->GetControlReactions()[j], -1);
851 }
852 }
853 }
854 }
855 //neutron current through the most outershape
856 Tally *GlobalLosses = new Tally(1);
857 if(gMURE->GetOutermostShape().Get() == 0)
858 {
859 LOG_ERROR("Most outer shape has not been specified for control material loss evaluation.Give it with gMURE->SetOutermostShape(Ext) where Ext is the exterior Shape.");
860 }
861 GlobalLosses->Add(gMURE->GetOutermostShape());
862 GlobalLosses->SetBinSurface(0, 1e-4); //SD card must be=1cm2(=1e-4m2)!
863 gMURE->GetEvolutionSolver()->SetGlobalLossesTallyNumber(GlobalLosses->GetNumber());
864 }
865
866 //Control Rod
867 if(gMURE->GetEvolutionSolver()->GetControlRodCell().size())
868 {
869 //Tally*GlobalTally=newTally[gMURE->GetEvolutionSolver()->GetControlRodCell().size()];
870 vector<Tally *> GlobalTally(gMURE->GetEvolutionSolver()->GetControlRodCell().size());
871 for(int i = 0; i < int(GlobalTally.size()); i++)
872 {
873 GlobalTally[i] = new Tally;
874 }
875
876 for(int i = 0; i < int(gMURE->GetEvolutionSolver()->GetControlRodCell().size()); i++)
877 {
878 PinCell *ControlRodPinCell = dynamic_cast<PinCell *>(gMURE->GetEvolutionSolver()->GetControlRodCell()[i]);
879 if (ControlRodPinCell != NULL)
880 {
881 int ControlLayerNumber = ControlRodPinCell->GetControlLayerNumber();
882 if(ControlLayerNumber == -1)
883 {
884 LOG_ERROR("One of the ControlRod is a PinCell and the layer containing control material in this PinCell has not been defined.Give it with gMURE->SetControlLayerNumber(i) where i is the number of the layer, 0 being the central one.");
885 }
886 GlobalTally[i]->Add(ControlRodPinCell, ControlLayerNumber);
887 }
888 else
889 {
890 GlobalTally[i]->Add(gMURE->GetEvolutionSolver()->GetControlRodCell()[i]);
891 }
892 gMURE->GetEvolutionSolver()->GetControlRodCell()[i]->SetGlobalTallyNumber(GlobalTally[i]->GetNumber());
893 GlobalTally[i]->SetBinVolume(0, 1e-6);
894 Material *Mat ;
895
896 if (ControlRodPinCell != NULL)
897 {
898 int ControlLayerNumber = ControlRodPinCell->GetControlLayerNumber();
899 Mat = ControlRodPinCell->GetLayerMaterial(ControlLayerNumber);
900 }
901 else
902 {
903 Mat = gMURE->GetEvolutionSolver()->GetControlRodCell()[i]->GetMaterial();
904 }
905
906 GlobalTally[i]->AddMultiplicator(Mat, -2, -1);
907 }
908
909 }
910
911 /*commented because not needed in a fully mirrored geometry and error prone...
912 Tally* GlobalLosses=new Tally(2);
913 if(fOutermostShape.Get()==0)
914 {
915 cerr<<"==========================================================="<<endl;
916 cerr<<"ERROR in MURE::BuildTallies()"<<endl;
917 cerr<<"\t Most outer shape has not been specified for fissile"<<endl;
918 cerr<<"\t control losse valuation.Give it with"<<endl;
919 cerr<<"\t gMURE->SetOutermostShape(Ext)where Ext is the exterior Shape"<<endl;
920 exit(1);
921 }
922 GlobalLosses->Add(GetOutermostShape());
923 GlobalLosses->SetBinSurface(0,1e-4);//SD card must be=1cm2!
924 fGlobalLossesTallyNumber=GlobalLosses->GetNumber();
925 */
926 //
927 // Genaral Energy Spectrum (use for a better understanding of the physics)
928 //
929 if(gMURE->IsESpectrum())
930 {
931 Tally *ESpectrum = new Tally(4);
932 for(int i = 0; i < (int)gMURE->GetTrueCellVector().size(); i++)
933 ESpectrum->Add(gMURE->GetTrueCellVector()[i]);
934 if(gMURE->IsUseEnergyBinsFile())
935 ESpectrum->AddEnergy(gMURE->GetEnergyBins());
936 else
937 ESpectrum->AddEnergy(10000, 0.00011, 19640330.00);
938
939 gMURE->SetESpectrumTallyNumber(ESpectrum->GetNumber());
940 }
941
942}
943
944//________________________________________________________________________
945void Connector::UpdateGlobalTallies()
946{
947
948 int NumberOfTallies = fMCOutput->GetNbTal();
949
950 if(gMURE->GetEvolutionSolver()->IsMaterialControl())
951 {
952
953 LOG_DEBUG("Material control: finding corresponding tallies...");
954 for(int i = 0 ; i < int(gMURE->GetTrueCellVector().size()); i++)
955 {
956 Cell *TheCell = gMURE->GetTrueCellVector()[i];
957 int GlobalTallyNum = TheCell->GetGlobalTallyNumber();
958
959 for(int j = 0; j < NumberOfTallies; j++)
960 {
961 TMTally *TheTally = fMCOutput->GetTally(j);
962 if(int(TheTally->fNo) == GlobalTallyNum)
963 {
964 //if(TheCell->IsFissile())
965 for(int m = 0; m < 5; m++) //Order of this multiplicator is not random!
966 {
967 // 0 = abs, 1= n2n, 2=n3n, 3=fiss, 4=Nu*fiss
968 ValErr_t UnNormalizedRate = TheTally->fVal[0][0][0][0][m][0][0][0];
969
970 if(m == 0) TheCell->SetCellAbsorptions(UnNormalizedRate);
971 if(m == 1) TheCell->SetCellN2N(UnNormalizedRate);
972 if(m == 2) TheCell->SetCellN3N(UnNormalizedRate);
973 if(m == 3) TheCell->SetCellFissions(UnNormalizedRate);
974 if(m == 4) TheCell->SetCellNuFissions(UnNormalizedRate);
975 }
976
977 if(dynamic_cast<ControlMaterial *>(TheCell->GetMaterial()))
978 {
979 LOG_DEBUG("Control material CELL FOUND");
980 vector<Reaction> React = ((ControlMaterial *)TheCell->GetMaterial())->GetControlReactions();
981 for(int m = 5 ; m < 5 + int(React.size()); m++)
982 {
983 ValErr_t UnNormalizedRate = TheTally->fVal[0][0][0][0][m][0][0][0];
984 TheCell->SetControlRate(React[m - 5], UnNormalizedRate);
985 }
986 }
987 }
988 if(int(TheTally->fNo) == gMURE->GetEvolutionSolver()->GetGlobalLossesTallyNumber())
989 gMURE->GetEvolutionSolver()->SetNeutronLosses(TheTally->fVal[0][0][0][0][0][0][0][0]);
990 }
991 }
992 gMURE->GetEvolutionSolver()->SetGlobalRates();
993 }
994
995
996 if(gMURE->GetEvolutionSolver()->GetControlRodCell().size())
997 {
998 for(int i = 0 ; i < int(gMURE->GetEvolutionSolver()->GetControlRodCell().size()); i++)
999 {
1000 Cell *TheCell = gMURE->GetEvolutionSolver()->GetControlRodCell()[i];
1001 int GlobalTallyNum = TheCell->GetGlobalTallyNumber();
1002
1003 for(int j = 0; j < NumberOfTallies; j++)
1004 {
1005 TMTally *TheTally = fMCOutput->GetTally(j);
1006 if(int(TheTally->fNo) == GlobalTallyNum)
1007 {
1008 ValErr_t UnNormalizedRate = TheTally->fVal[0][0][0][0][0][0][0][0];
1009 TheCell->SetCellAbsorptions(UnNormalizedRate);
1010 }
1011 }
1012 }
1013 }
1014
1015 if(gMURE->IsESpectrum())
1016 {
1017
1018 int NumberOfCells = gMURE->GetTrueCellVector().size();
1019 double Zero = 1.e-18 ;
1020 //WriteEnergySpectrum()
1021 stringstream name;
1022 name << "ESPECTRUM";
1023 string filename = gMURE->GetMCRunDirectory() + "/" + name.str();
1024 ofstream eoutfile(filename.c_str(), ios::app);
1025
1026 for(int j = 0; j < NumberOfTallies; j++)
1027 {
1028 TMTally *TheTally = fMCOutput->GetTally(j);
1029 if(int(TheTally->fNo) == gMURE->GetESpectrumTallyNumber())
1030 {
1031 LOG_DEBUG("ESpectrumTallyNum = " << gMURE->GetESpectrumTallyNumber());
1032 float *EnergyBins = TheTally->fEV;
1033
1034
1035 for(int f = 0; f < NumberOfCells; f++)
1036 {
1037 for(int e = 0; e < int(TheTally->fEN) - 1; e++)
1038 {
1039 double Flux_E = TheTally->fVal[f][0][0][0][0][0][e][0].fVal;
1040 double Flux_E_err = TheTally->fVal[f][0][0][0][0][0][e][0].fErr;
1041 if ( Flux_E == 0. ) Flux_E = Zero;
1042 eoutfile << EnergyBins[e] << " " << Flux_E << " " << Flux_E_err << endl;
1043 }
1044 eoutfile << "&" << endl;
1045 }
1046 }
1047 }
1048 }
1049
1050}
1051
1052//________________________________________________________________________
1053void Connector::BuildTHMultiGroupDetectors()
1054{
1055 LOG_INFO("---------------------------------------------------------------------------");
1056 LOG_INFO("Thermal Hydraulics Tallies building.");
1057 LOG_IMP_WARN("No perturbative materials built : this method must be used carefully");
1058
1059 vector<Material *> MaterialVector = gMURE->GetMaterialVector();
1060 fTHMultigroupTallies = true;
1061 BuildMultiGroupDetectors(MaterialVector);
1062
1063 LOG_INFO("Tallies built.");
1064 LOG_INFO("---------------------------------------------------------------------------");
1065}
1066
1067//________________________________________________________________________
1068void Connector::UpdateTHMultiGroupSigmaPhiDetectors()
1069{
1070 LOG_INFO("---------------------------------------------------------------------------");
1071 LOG_INFO("Updating Thermal Hydraulics SigmaPhiMultiGroup from " << GetMCDetectorOutputFileName(gMURE->GetRealMCInputFileName()) << " ...");
1072
1073 vector<Material *> MaterialVector = gMURE->GetMaterialVector();
1074 UpdateMultiGroupSigmaPhiDetectors(MaterialVector);
1075
1076 LOG_INFO("All Thermal Hydraulics Tallies read and updated.");
1077 LOG_INFO("---------------------------------------------------------------------------");
1078}
1079
1080
1081
1082}
MCNPConnector.hxx
Header file for MCNP::Connector class.
MCNPCylinder.hxx
Header file for MCNP::Cylinder class.
MCNPNode.hxx
Header file for MCNP::Node class.
MCNPTally.hxx
Header file for MCNP SimpleBin, GroupBin, LatticeBin and Tally classes.
LOG_DEBUG
#define LOG_DEBUG(msg)
Definition MUREGlobal.hxx:45
LOG_ERROR
#define LOG_ERROR(msg)
Print an error and terminate the program by calling exit(999).
Definition MUREGlobal.hxx:55
gMURE
std::unique_ptr< MURE > gMURE
LOG_IMP_WARN
#define LOG_IMP_WARN(msg)
Print a important warning (something happened, but MURE took care)
Definition MUREGlobal.hxx:51
LOG_INFO
#define LOG_INFO(msg)
Print info message (Now MURE is doing this ...)
Definition MUREGlobal.hxx:47
MiscFunction.hxx
Header file for vector operation, Regression and FindFitParameters functions.
ReadXSFile.hxx
Header file for ReadXSFile class.
Cell
A Cell is composed from a Shape and a Material.
Definition Cell.hxx:84
Cell::GetNEnergyGroup
int GetNEnergyGroup()
Definition Cell.hxx:350
Cell::SetCellFissions
void SetCellFissions(ValErr_t rate)
Definition Cell.hxx:400
Cell::SetMultiGroupFlux
void SetMultiGroupFlux(int i, ValErr_t Flux)
Definition Cell.hxx:338
Cell::SetFlux
void SetFlux(double phi)
Definition Cell.hxx:320
Cell::GetNumber
int GetNumber()
Definition Cell.hxx:158
Cell::GetTallyNumber
int GetTallyNumber()
Definition Cell.hxx:276
Cell::SetCellN3N
void SetCellN3N(ValErr_t rate)
Definition Cell.hxx:396
Cell::SetEnergyGroup
void SetEnergyGroup(vector< float > &nrj)
Definition Cell.hxx:354
Cell::SetMCFlux
void SetMCFlux(ValErr_t phi)
Definition Cell.hxx:324
Cell::SetCellN2N
void SetCellN2N(ValErr_t rate)
Definition Cell.hxx:392
Cell::SetTallyNumber
void SetTallyNumber(int num)
Definition Cell.hxx:272
Cell::GetMultiGroupFlux
ValErr_t GetMultiGroupFlux(int i)
Definition Cell.hxx:342
Cell::GetMaterial
Material * GetMaterial()
Definition Cell.hxx:135
Cell::SetCellNuFissions
void SetCellNuFissions(ValErr_t rate)
Definition Cell.hxx:404
Cell::GetEnergyGroups
vector< float > & GetEnergyGroups()
Definition Cell.hxx:359
Cell::SetControlRate
void SetControlRate(Reaction r, ValErr_t rate)
Definition Cell.hxx:417
Cell::GetTallyBinNumber
int GetTallyBinNumber()
Definition Cell.hxx:284
Cell::BuildMultiGroupFlux
void BuildMultiGroupFlux()
Build the group flux for multigroup run.
Definition Cell.cxx:720
Cell::SetCellAbsorptions
void SetCellAbsorptions(ValErr_t rate)
Definition Cell.hxx:388
Cell::SetTallyBinNumber
void SetTallyBinNumber(int num)
Definition Cell.hxx:280
Cell::GetGlobalTallyNumber
int GetGlobalTallyNumber()
Definition Cell.hxx:383
ConnectorPlugin::fRebuildDetector
bool fRebuildDetector
whether or not rebuild detectors (for thermo-hydraulics)
Definition ConnectorPlugin.hxx:186
ConnectorPlugin::fTHMultigroupTallies
bool fTHMultigroupTallies
True in Thermal hydraulics multigroup calculation.
Definition ConnectorPlugin.hxx:192
ControlMaterial
A ControlMaterial.
Definition ControlMaterial.hxx:51
MCNP::Connector::BuildTHMultiGroupDetectors
void BuildTHMultiGroupDetectors() override
Definition MCNPConnector_detectors_legacy.cxx:1053
MCNP::Connector::fMCOutput
TMctal * fMCOutput
A result tallies from a MCNP m's files.
Definition MCNPConnector.hxx:234
MCNP::Connector::UpdateKeff
void UpdateKeff() override
read from MC output the keff and give it to MURE
Definition MCNPConnector_detectors_legacy.cxx:773
MCNP::Connector::UpdateTHMultiGroupSigmaPhiDetectors
void UpdateTHMultiGroupSigmaPhiDetectors() override
Definition MCNPConnector_detectors_legacy.cxx:1068
MCNP::Connector::BuildGlobalTallies
void BuildGlobalTallies()
build global tallies for global reaction rates (Control Rod, ...) for special evolution controls
Definition MCNPConnector_detectors_legacy.cxx:797
MCNP::Connector::UpdateSigmaPhiDetectors
void UpdateSigmaPhiDetectors() override
update (read from MC output) automatic tallies for standard evolution
Definition MCNPConnector_detectors_legacy.cxx:149
MCNP::Connector::BuildDetectors
void BuildDetectors() override
build automatic tallies for standard evolution
Definition MCNPConnector_detectors_legacy.cxx:39
MCNP::Connector::BuildMultiGroupDetectors
void BuildMultiGroupDetectors() override
build automatic tallies for in a multigroup evolution
Definition MCNPConnector_detectors_legacy.cxx:272
MCNP::Connector::ReadDetectorFile
void ReadDetectorFile()
Definition MCNPConnector_detectors_legacy.cxx:485
MCNP::Connector::UpdateGlobalTallies
void UpdateGlobalTallies()
update the global tallies
Definition MCNPConnector_detectors_legacy.cxx:945
MCNP::Connector::UpdateMultiGroupSigmaPhiDetectors
void UpdateMultiGroupSigmaPhiDetectors() override
update (read from MC output) automatic tallies for multigroup evolution
Definition MCNPConnector_detectors_legacy.cxx:502
MCNP::Connector::GetMCDetectorOutputFileName
string GetMCDetectorOutputFileName(string InputFileName, bool OnlyKeff=false) override
Definition MCNPConnector.hxx:90
MCNP::Tally
Define a MCNP Tally.
Definition MCNPTally.hxx:161
MCNP::Tally::SetBinSurface
void SetBinSurface(int i, double S)
set the ith bin surface of Tally bin
Definition MCNPTally.cxx:399
Material
A Material constituing a Cell.
Definition Material.hxx:83
Material::GetTemperature
double GetTemperature()
Definition Material.hxx:271
Material::GetComposition
vector< Nucleus_ptr > & GetComposition()
Definition Material.hxx:182
MureTally::SetPrintable
void SetPrintable(bool flag=true)
Definition MureTally.hxx:282
MureTally::AddMultiGroupEnergy
void AddMultiGroupEnergy(double Emin, double Emax, double decad_mult=1)
Definition MureTally.cxx:320
MureTally::GetNumber
int GetNumber()
Definition MureTally.hxx:154
MureTally::AddEnergy
void AddEnergy(int NE, double *E)
Definition MureTally.cxx:236
MureTally::Add
void Add(Cell *C, bool SetCellTally=true)
Add a cell bin.
Definition MureTally.cxx:438
NucleusMCRecord
A NucleusMCRecord allows to record data for writing &reading reaction rates tallies.
Definition NucleusMCRecord.hxx:54
NucleusMCRecord::GetMultiplicatorBinIndex
vector< int > & GetMultiplicatorBinIndex()
Definition NucleusMCRecord.hxx:87
NucleusMCRecord::GetTallyNumber
int GetTallyNumber()
Definition NucleusMCRecord.hxx:73
NucleusMCRecord::AddMultiplicatorBinIndex
vector< int > AddMultiplicatorBinIndex(int &FMBinIdx, Nucleus *nuc)
Set Multiplicator tally bin index for reaction rates.
Definition NucleusMCRecord.cxx:67
NucleusMCRecord::GetNumberOfMultiplicatorBin
int GetNumberOfMultiplicatorBin()
Definition NucleusMCRecord.hxx:110
NucleusMCRecord::GetCellNumber
int GetCellNumber(int i)
Definition NucleusMCRecord.hxx:115
NucleusMCRecord::SetTotalXSMultiplicatorBinIndex
void SetTotalXSMultiplicatorBinIndex(int FMbin)
Definition NucleusMCRecord.hxx:91
NucleusMCRecord::GetTallyBinIndex
int GetTallyBinIndex(int i)
Definition NucleusMCRecord.hxx:119
NucleusMCRecord::GetTotalXSMultiplicatorBinIndex
int GetTotalXSMultiplicatorBinIndex()
Definition NucleusMCRecord.hxx:95
NucleusMCRecord::GetNumberOfCellBin
int GetNumberOfCellBin()
Definition NucleusMCRecord.hxx:106
PinCell
PinCell class allows to create cylindrical cell set included as Matrioshka.
Definition PinCell.hxx:53
PinCell::GetLayerMaterial
Material * GetLayerMaterial(int i)
return the material of the layer i (i=0, innerest layer)
Definition PinCell.cxx:162
PinCell::GetControlLayerNumber
int GetControlLayerNumber()
Definition PinCell.hxx:226
Reaction
Define a Reaction list for Tally multiplicator inputs.
Definition Reaction.hxx:41
Reaction::Multiply
void Multiply(int code)
multiply reaction (e.g. nu*fission)
Definition Reaction.cxx:77
ReadXSFile
Allows to read XS file in the MCNP ACE format.
Definition ReadXSFile.hxx:83
ReadXSFile::NuSigmaFisPhi
ValErr_t NuSigmaFisPhi(int N, vector< ValErr_t > &Phi, float *Ephi)
return
Definition ReadXSFile.cxx:770
ReadXSFile::SigmaPhi
ValErr_t SigmaPhi(int N, vector< ValErr_t > &Phi, float *Ephi, int MT, double emin=0, double emax=1e40)
return
Definition ReadXSFile.cxx:794
Reference_ptr
Handle dynamical object creation and pointer affectation.
Definition TReference.hxx:74
Reference_ptr::Get
T * Get() const
Definition TReference.hxx:102
XSDIRLine
Extract all parameters from an XSDIR line.
Definition XSDIRLine.hxx:43
XSDIRLine::GetTableLength
int GetTableLength()
Definition XSDIRLine.hxx:79
XSDIRLine::GetDataType
int GetDataType()
Definition XSDIRLine.hxx:71
XSDIRLine::GetStartRecord
int GetStartRecord()
Definition XSDIRLine.hxx:75
XSDIRLine::GetRecordLength
int GetRecordLength()
Definition XSDIRLine.hxx:83
XSDIRLine::GetXSFileName
string GetXSFileName()
Definition XSDIRLine.hxx:131
XSDIRLine::GetEntriesPerRecord
int GetEntriesPerRecord()
Definition XSDIRLine.hxx:87
MCNP
This MCNP (concrete) class is used to built a reactor assembly.
Definition MCNPBrick.hxx:41
std
the namespace of the Standard C++
MURE Project, documentation generated by Doxygen 1.9.7 - Fri Jan 19 2024