diff --git a/source/trunk/include/DataBank.hxx b/source/trunk/include/DataBank.hxx
index 551dea050b9c1e906df1e0b00111271188d7a1d8..8bd594eb06131324d869b0d975abe8bd467adb1b 100755
--- a/source/trunk/include/DataBank.hxx
+++ b/source/trunk/include/DataBank.hxx
@@ -83,6 +83,7 @@ public :
void LoadFPYield(string SponfaneusYield, string ReactionYield); // Build Fision Yields maps;
void SetWeightedDistanceCalculation(bool val=true) { fWeightedDistance = val;} // Set weighted Distance calculation
+ void SetEvolutionDataInterpolation(bool val=true) { fEvolutionDataInterpolation = val;} // Set weighted Distance calculation
//********* Modification Method *********//
@@ -139,6 +140,7 @@ protected :
bool fUseRK4EvolutionMethod; ///< if true use RK4 calculation, mtriciel unstead
bool fOldReadMethod; ///< use old DB format
bool fWeightedDistance; ///< USe XS weighted distance calculation
+ bool fEvolutionDataInterpolation; ///< USe XS weighted distance calculation
string fFuelType;
diff --git a/source/trunk/include/EvolutionData.hxx b/source/trunk/include/EvolutionData.hxx
index 35c885eb4cb5ffbfd3e01cb677c85b5e228cf8e9..4fb36cb9b673812f6d5e038e1c975b526cd89175 100755
--- a/source/trunk/include/EvolutionData.hxx
+++ b/source/trunk/include/EvolutionData.hxx
@@ -31,8 +31,10 @@ typedef long long int cSecond;
EvolutionData operator*(EvolutionData const& evol, double F);
EvolutionData operator*(double F, EvolutionData const& evol);
EvolutionData operator/(EvolutionData const& evol, double F);
+EvolutionData operator+(EvolutionData const& evol1, EvolutionData const& evol2);
-
+double Distance(IsotopicVector IV1, EvolutionData Evd1 );
+double Distance(EvolutionData Evd1, IsotopicVector IV1 );
class EvolutionData : public CLSSObject
{
@@ -57,6 +59,9 @@ public :
void SetFuelType(string fueltype) { fFuelType = fueltype; }
void SetPower(double power) { fPower = power; }
void SetFlux(TGraph* flux ) { fFlux = flux; }
+
+
+ void SetEvolutionData(map<ZAI, TGraph*> maptoinsert) { fEvolutionData = maptoinsert;}
void SetFissionXS(map<ZAI, TGraph*> maptoinsert) { fFissionXS = maptoinsert;}
void SetCaptureXS(map<ZAI, TGraph*> maptoinsert) { fCaptureXS = maptoinsert;}
void Setn2nXS(map<ZAI, TGraph*> maptoinsert) { fn2nXS = maptoinsert;}
@@ -146,8 +151,7 @@ protected :
ClassDef(EvolutionData,0);
};
-double Distance(IsotopicVector IV1, EvolutionData Evd1 );
-double Distance(EvolutionData Evd1, IsotopicVector IV1 );
+
#endif
diff --git a/source/trunk/src/DataBank.cxx b/source/trunk/src/DataBank.cxx
index 8c416add4a2f17934e12ef013c7803df8e6e4b30..61009214f6b8e905d4ee2761464eb4914979c4dc 100755
--- a/source/trunk/src/DataBank.cxx
+++ b/source/trunk/src/DataBank.cxx
@@ -206,6 +206,9 @@ DataBank<IsotopicVector>::DataBank():DynamicalSystem()
fTheMatrix = 0;
fWeightedDistance = false;
+ fEvolutionDataInterpolation = false;
+
+
fOldReadMethod = true;
fUseRK4EvolutionMethod = true;
@@ -230,6 +233,9 @@ DataBank<IsotopicVector>::DataBank(LogFile* Log, string DB_index_file, bool setl
IsLog(setlog);
fWeightedDistance = false;
+ fEvolutionDataInterpolation = false;
+
+
fTheNucleiVector = 0;
fTheMatrix = 0;
@@ -1297,12 +1303,59 @@ EvolutionData DataBank<IsotopicVector>::GetClosest(IsotopicVector isotopicvector
map<IsotopicVector, EvolutionData > evolutiondb = fDataBank;
double distance = 0;
+ map<IsotopicVector, EvolutionData >::iterator it_close = evolutiondb.begin();
+
+
+ map<IsotopicVector, EvolutionData >::iterator it;
+
+
if(fWeightedDistance)
{
Distance(isotopicvector.GetActinidesComposition()
* evolutiondb.begin()->second.GetIsotopicVectorAt(t).GetActinidesComposition().GetSumOfAll()
/ isotopicvector.GetActinidesComposition().GetSumOfAll(),
evolutiondb.begin()->second);
+
+
+ for( it = evolutiondb.begin(); it != evolutiondb.end(); it++ )
+ {
+ double D = 0;
+ D = Distance(isotopicvector.GetActinidesComposition()
+ * evolutiondb.begin()->second.GetIsotopicVectorAt(t).GetActinidesComposition().GetSumOfAll()
+ / isotopicvector.GetActinidesComposition().GetSumOfAll(),
+ evolutiondb.begin()->second);
+
+
+ if (D< distance)
+ {
+ distance = D;
+ it_close = it;
+ }
+ }
+
+ return (*it_close).second;
+ }
+ else if (fEvolutionDataInterpolation)
+ {
+ map<double, EvolutionData> distance_map;
+ map<double, EvolutionData>::iterator it_distance;
+ int NClose = 64;
+ int Nstep = 0;
+ EvolutionData EvolInterpolate;
+ double SumOfDistance = 0;
+ for( it_distance = distance_map.begin(); Nstep < NClose; it_distance++)
+ {
+ if(Nstep == 0)
+ EvolInterpolate = 1./(*it_distance).first * (*it_distance).second;
+ else
+ EvolInterpolate = EvolInterpolate + 1./(*it_distance).first * (*it_distance).second;
+
+ SumOfDistance += 1./(*it_distance).first;
+ Nstep++;
+
+ }
+ return 1/SumOfDistance * EvolInterpolate;
+
}
else
{
@@ -1311,39 +1364,26 @@ EvolutionData DataBank<IsotopicVector>::GetClosest(IsotopicVector isotopicvector
/ evolutiondb.begin()->second.GetIsotopicVectorAt(t).GetActinidesComposition().GetSumOfAll()
* isotopicvector.GetActinidesComposition().GetSumOfAll(),
fDistanceType, fDistanceParameter);
- }
-
- map<IsotopicVector, EvolutionData >::iterator it_close = evolutiondb.begin();
- map<IsotopicVector, EvolutionData >::iterator it;
- for( it = evolutiondb.begin(); it != evolutiondb.end(); it++ )
- {
double D = 0;
- if(fWeightedDistance)
- {
- D = Distance(isotopicvector.GetActinidesComposition()
- * evolutiondb.begin()->second.GetIsotopicVectorAt(t).GetActinidesComposition().GetSumOfAll()
- / isotopicvector.GetActinidesComposition().GetSumOfAll(),
- evolutiondb.begin()->second);
- }
- else
- {
- D = Distance(isotopicvector.GetActinidesComposition(),
- evolutiondb.begin()->second.GetIsotopicVectorAt(t).GetActinidesComposition()
- / evolutiondb.begin()->second.GetIsotopicVectorAt(t).GetActinidesComposition().GetSumOfAll()
- * isotopicvector.GetActinidesComposition().GetSumOfAll(),
- fDistanceType, fDistanceParameter);
- }
+ D = Distance(isotopicvector.GetActinidesComposition(),
+ evolutiondb.begin()->second.GetIsotopicVectorAt(t).GetActinidesComposition()
+ / evolutiondb.begin()->second.GetIsotopicVectorAt(t).GetActinidesComposition().GetSumOfAll()
+ * isotopicvector.GetActinidesComposition().GetSumOfAll(),
+ fDistanceType, fDistanceParameter);
+
if (D< distance)
{
distance = D;
it_close = it;
}
+ return (*it_close).second;
+
}
- return (*it_close).second;
+
}
@@ -1366,24 +1406,26 @@ void DataBank<IsotopicVector>::CalculateDistanceParameter()
//We calculate the weight for the distance calculation.
map<IsotopicVector ,EvolutionData >::iterator it;
map<IsotopicVector ,EvolutionData > databank = (*this).GetDataBank();
- int NevolutionDatainDataBank=0;
+ int NevolutionDatainDataBank = 0;
- for( it = databank.begin(); it != databank.end(); it++ ){
+ for( it = databank.begin(); it != databank.end(); it++ )
+ {
NevolutionDatainDataBank++;
map<ZAI ,double>::iterator itit;
map<ZAI ,double> isovector=(*it).first.GetIsotopicQuantity();
- for(itit=isovector.begin(); itit != isovector.end(); itit++){//Boucle sur ZAI
- ZAI TmpZAI=(*itit).first;
+ for(itit=isovector.begin(); itit != isovector.end(); itit++) //Boucle sur ZAI
+ {
double TmpXS=0;
- for(int i=1;i<4;i++){ //Loop on Reactions 1==fission, 2==capture, 3==n2n
- TmpXS+= (*it).second.GetXSForAt(0,TmpZAI,i);
- }
- fDistanceParameter.Add(TmpZAI,TmpXS);
+
+ for( int i=1; i<4; i++ ) //Loop on Reactions 1==fission, 2==capture, 3==n2n
+ TmpXS+= (*it).second.GetXSForAt(0, (*itit).first, i);
+
+ fDistanceParameter.Add((*itit).first,TmpXS);
}
}
- fDistanceParameter.Multiply((double)1.0/NevolutionDatainDataBank);
+ fDistanceParameter.Multiply( (double)1.0/NevolutionDatainDataBank );
GetLog()->fLog <<"!!INFO!! Distance Parameters "<<endl;
diff --git a/source/trunk/src/EvolutionData.cxx b/source/trunk/src/EvolutionData.cxx
index 4e1a82878f34bc84771803388972d9e58a717588..52f893a3e015e26f607da29793bfb53bf1634cae 100755
--- a/source/trunk/src/EvolutionData.cxx
+++ b/source/trunk/src/EvolutionData.cxx
@@ -110,6 +110,130 @@ EvolutionData operator*(EvolutionData const& evol, double F)
return evoltmp;
}
+EvolutionData operator+(EvolutionData const& evol1, EvolutionData const& evol2)
+{
+
+
+ EvolutionData EvolSum = evol1;
+ map<ZAI ,TGraph* > EvolutionData1 = evol1.GetEvolutionData();
+ map<ZAI ,TGraph* > EvolutionData2 = evol2.GetEvolutionData();
+ map<ZAI ,TGraph* >::iterator it;
+
+ for(it = EvolutionData2.begin(); it != EvolutionData2.end(); it++)
+ {
+ pair<map<ZAI, TGraph*>::iterator, bool> IResult;
+
+ IResult = EvolutionData1.insert( pair<ZAI, TGraph*> ( *it ) );
+
+ if(!(IResult.second) )
+ {
+ double X[(*it).second->GetN()];
+ double Y[(*it).second->GetN()];
+ map<ZAI ,TGraph* >::iterator it2 = EvolutionData1.find( (*it).first );
+
+
+ for(int i = 0; i < (*it).second->GetN(); i++)
+ {
+ double y;
+ (*it).second->GetPoint( i, X[i], y );
+ Y[i] = y + (*it2).second->Eval(X[i]);
+ }
+ (*it).second = new TGraph((*it).second->GetN(), X, Y);
+ }
+ }
+
+ EvolSum.SetEvolutionData(EvolutionData1);
+
+
+ EvolutionData1 = evol1.GetFissionXS();
+ EvolutionData2 = evol2.GetFissionXS();
+
+ for(it = EvolutionData2.begin(); it != EvolutionData2.end(); it++)
+ {
+ pair<map<ZAI, TGraph*>::iterator, bool> IResult;
+
+ IResult = EvolutionData1.insert( pair<ZAI, TGraph*> ( *it ) );
+
+ if(!(IResult.second) )
+ {
+ double X[(*it).second->GetN()];
+ double Y[(*it).second->GetN()];
+ map<ZAI ,TGraph* >::iterator it2 = EvolutionData1.find( (*it).first );
+
+
+ for(int i = 0; i < (*it).second->GetN(); i++)
+ {
+ double y;
+ (*it).second->GetPoint( i, X[i], y );
+ Y[i] = y + (*it2).second->Eval(X[i]);
+ }
+ (*it).second = new TGraph((*it).second->GetN(), X, Y);
+ }
+ }
+ EvolSum.SetFissionXS(EvolutionData1);
+
+
+
+ EvolutionData1 = evol1.GetCaptureXS();
+ EvolutionData2 = evol2.GetCaptureXS();
+
+ for(it = EvolutionData2.begin(); it != EvolutionData2.end(); it++)
+ {
+ pair<map<ZAI, TGraph*>::iterator, bool> IResult;
+
+ IResult = EvolutionData1.insert( pair<ZAI, TGraph*> ( *it ) );
+
+ if(!(IResult.second) )
+ {
+ double X[(*it).second->GetN()];
+ double Y[(*it).second->GetN()];
+ map<ZAI ,TGraph* >::iterator it2 = EvolutionData1.find( (*it).first );
+
+
+ for(int i = 0; i < (*it).second->GetN(); i++)
+ {
+ double y;
+ (*it).second->GetPoint( i, X[i], y );
+ Y[i] = y + (*it2).second->Eval(X[i]);
+ }
+ (*it).second = new TGraph((*it).second->GetN(), X, Y);
+ }
+ }
+ EvolSum.SetCaptureXS(EvolutionData1);
+
+
+
+ EvolutionData1 = evol1.Getn2nXS();
+ EvolutionData2 = evol2.Getn2nXS();
+
+ for(it = EvolutionData2.begin(); it != EvolutionData2.end(); it++)
+ {
+ pair<map<ZAI, TGraph*>::iterator, bool> IResult;
+
+ IResult = EvolutionData1.insert( pair<ZAI, TGraph*> ( *it ) );
+
+ if(!(IResult.second) )
+ {
+ double X[(*it).second->GetN()];
+ double Y[(*it).second->GetN()];
+ map<ZAI ,TGraph* >::iterator it2 = EvolutionData1.find( (*it).first );
+
+
+ for(int i = 0; i < (*it).second->GetN(); i++)
+ {
+ double y;
+ (*it).second->GetPoint( i, X[i], y );
+ Y[i] = y + (*it2).second->Eval(X[i]);
+ }
+ (*it).second = new TGraph((*it).second->GetN(), X, Y);
+ }
+ }
+ EvolSum.Setn2nXS(EvolutionData1);
+
+
+ return EvolSum;
+}
+
//________________________________________________________________________
EvolutionData operator*(double F, EvolutionData const& evol)