CLASS  1.1
src/CLASS.cxx
Aller à la documentation de ce fichier.
00001 #include "CLASS.hxx"
00002 #include "Reactor.hxx"
00003 #include "TreatmentFactory.hxx"
00004 #include "Defines.hxx"
00005 #include "IsotopicVector.hxx"
00006 #include "LogFile.hxx"
00007 
00008 #include <ctime>
00009 #include "time.h"
00010 #include <cmath>
00011 #include <iomanip>
00012 #include <fstream>
00013 #include <iostream>
00014 #include <algorithm>
00015 #include <omp.h>
00016 #include <TFile.h>
00017 #include <TTree.h>
00018 //#include <gsl/gsl_vector.h>
00019 //#include <gsl/gsl_multimin.h>
00020 
00021 
00022 //________________________________________________________________________
00023 //
00024 //              CLASS
00025 //
00026 //
00027 //
00028 //
00029 //________________________________________________________________________
00030 
00031 
00032 float random(float a, float b) //peak random numebr between a and b
00033 {
00034     float range = pow(2., 31);
00035     srand(time(NULL)); //initialize the srand
00036     return (float)a + (float)(b-a)*rand()/range;
00037 }
00038 
00039 //________________________________________________________________________
00040 CLASS::CLASS()
00041 {
00042 DBGL;
00043         fPrintStep = (long int)(3600*24*365.4);  // One Step per Month
00044         fAbsoluteTime = 0;
00045         fTreatmentFactoryLastIndex = 0;
00046         fReactorLastIndex = 0;
00047         fStockManagement = true;
00048         fBuildingMethod = 0;
00049         
00050         fOutputName = "CLASS_Default.root";
00051         string logname = "CLASS.log";
00052         fLog = new LogFile("CLASS.log");
00053 
00054 DBGL;
00055 }
00056 
00057 //________________________________________________________________________
00058 CLASS::CLASS(long int abstime)
00059 {
00060 DBGL;
00061         fPrintStep = (long int)(3600*24*365.4);  // One Step per Month
00062         fAbsoluteTime = abstime;
00063         fTreatmentFactoryLastIndex = 0;
00064         fReactorLastIndex = 0;
00065         fStockManagement = true;
00066         fBuildingMethod = 0;
00067         
00068         fOutputName = "CLASS_Default.root";
00069         string logname = "CLASS.log";
00070         fLog = new LogFile("CLASS.log");
00071 
00072 DBGL;
00073 }
00074 
00075 //________________________________________________________________________
00076 CLASS::CLASS(string name, long int abstime)
00077 {
00078 DBGL;
00079         fPrintStep = (long int)(3600*24*365.4);  // One Step per Month
00080         fAbsoluteTime = abstime;
00081         fTreatmentFactoryLastIndex = 0;
00082         fReactorLastIndex = 0;
00083         fStockManagement = true;
00084         fBuildingMethod = 0;
00085         fOutputName = name;
00086         string logname = "CLASS.log";
00087         fLog = new LogFile("CLASS.log");
00088 
00089 DBGL;
00090 }
00091 //________________________________________________________________________
00092 CLASS::~CLASS()
00093 {
00094 DBGL;
00095 DBGL;
00096 }
00097 
00098 //________________________________________________________________________
00099 void CLASS::AddTreatmentFactory(TreatmentFactory* treatmentfactory)
00100 {
00101 DBGL;
00102         
00103         fTreatmentFactory.push_back(treatmentfactory);
00104         fTreatmentFactory.back()->SetParc(this);
00105         fTreatmentFactory.back()->SetLog(fLog);
00106         fTreatmentFactory.back()->SetStockManagement(fStockManagement);
00107         fTreatmentFactoryLastIndex++;
00108         fTreatmentFactoryIndex.push_back(fTreatmentFactoryLastIndex);
00109 DBGL;
00110 }
00111 
00112 //________________________________________________________________________
00113 void CLASS::AddReactor(Reactor* reactor)
00114 {
00115 DBGL;
00116         fReactor.push_back(reactor);
00117         fReactor.back()->SetParc(this);
00118         fReactor.back()->SetLog(fLog);
00119         fReactorLastIndex++;
00120         fReactorIndex.push_back(fReactorLastIndex);
00121 DBGL;
00122 }
00123 
00124 //________________________________________________________________________
00125 void CLASS::RemoveReactor()
00126 {
00127 DBGL;
00128         for(int i = (int)fReactor.size()-1; i >= 0; i--)
00129         {
00130                 if(fAbsoluteTime == fReactor[i]->GetCreationTime() + fReactor[i]->GetLifeTime())
00131                 {
00132                         fReactor[i]->GetAssociedTreatmentFactory()->AddIVCooling(fReactor[i]->GetIVReactor());
00133                         delete (fReactor[i]);
00134                         fReactor.erase(fReactor.begin()+i);
00135                         fReactorIndex.erase(fReactorIndex.begin()+i);
00136                 }
00137         }
00138 DBGL;
00139 }
00140 
00141 //________________________________________________________________________
00142 //____________________________ Building Method ___________________________
00143 //________________________________________________________________________
00144 
00145 void CLASS::UpdateParcStock()
00146 {
00147 DBGL;
00148         
00149         fParcStock.clear();
00150         
00151         for(int i=0; i< (int)fTreatmentFactory.size(); i++ )
00152         {
00153                 for(int j = 0; j < (int)fTreatmentFactory[i]->GetIVStock().size(); j++ )
00154                 {
00155                         fParcStock.insert( pair< pair< int ,int > , IsotopicVector > 
00156                                 ( pair< int ,int >(i, j), fTreatmentFactory[i]->GetIVStock()[j]) ) ;
00157                 }
00158         }
00159         
00160 DBGL;
00161 }
00162 
00163 //________________________________________________________________________
00164 IsotopicVector CLASS::BuildIsotopicVector(IsotopicVector isotopicvector)
00165 {
00166 DBGL;
00167         UpdateParcStock();
00168         
00169         IsotopicVector BuildIV;
00170         if (fStockManagement == false )
00171                 return BuildIV;
00172 
00173         if(fBuildingMethod == 0)
00174                 return MonteCarloBuild(isotopicvector);
00175         
00176 //      else if(fBuildingMethod == 1)
00177 //              return MinimizationBuild(isotopicvector);
00178         else 
00179         {
00180                 cout            << "!!Warning!! !!!CLASS!!! Wrong BUilding Method !\n" << endl;
00181                 fLog->fLog      << "!!Warning!! !!!CLASS!!! Wrong BUilding Method !\n" << endl;
00182                 exit(1);
00183         }
00184         
00185 DBGL;
00186 }
00187 
00188 
00189 //________________________________________________________________________
00190 void CLASS::DumpParcStock()
00191 {
00192 DBGL;
00193         if(fStockManagement == false) return;
00194         for(int i = 0; i < (int)fTreatmentFactory.size(); i++)
00195                 fTreatmentFactory[i]->ClearIVStock();
00196         
00197         
00198         map< pair<int,int>, IsotopicVector >::iterator it;
00199         for(it = fParcStock.begin(); it != fParcStock.end(); it++)
00200                 fTreatmentFactory[(*it).first.first]->AddIVStock((*it).second);
00201         
00202 DBGL;   
00203 }
00204 
00205 
00206 
00207 
00208 
00209 //________________________________________________________________________
00210 
00211 //________________________________________________________________________
00212 
00213 /*struct MinimizationParameter {
00214         int z;
00215         vector<int> FactoryIndex;
00216         vector<int> StockIndex;
00217         map< pair<int,int>, IsotopicVector > ParcStock;
00218         IsotopicVector isotopicvector;
00219         IsotopicVector TotalStock;
00220 
00221 };
00222 
00223 
00224 double  my_f (const gsl_vector *v, void *params)
00225 {
00226         
00227         MinimizationParameter *par = (MinimizationParameter *)params;
00228         int z = (*par).z;
00229         vector<int> FactoryIndex = (*par).FactoryIndex;
00230         vector<int> StockIndex   = (*par).StockIndex;
00231         map< pair<int,int>, IsotopicVector > ParcStock = (*par).ParcStock;
00232         IsotopicVector isotopicvector = (*par).isotopicvector;
00233         IsotopicVector TotalStock = (*par).TotalStock;
00234 
00235         double x[FactoryIndex.size()];
00236         
00237                 IsotopicVector IVSum;
00238         for (int i =0; i < (int)FactoryIndex.size(); i++)
00239         {
00240                 x[i] = gsl_vector_get(v, i);
00241                 if(x[i]<0) x[i] = -x[i];
00242         }
00243 
00244         
00245         for (int i = 0; i < (int)FactoryIndex.size(); i++)
00246         {
00247                 map< pair<int,int>, IsotopicVector >::iterator it;
00248                 it = ParcStock.find(pair<int,int>(FactoryIndex[i], StockIndex[i] ));
00249                 IVSum += (*it).second.GetAtomicComposition(z)*x[i]/Norme((*it).second.GetAtomicComposition(z));
00250         }
00251         return Distance(IVSum,isotopicvector.GetAtomicComposition(z)); 
00252 }
00253 
00254 // The gradient of f, df = (df/dx, df/dy). 
00255 void my_df (const gsl_vector *v, void *params, gsl_vector *df)
00256 {
00257         MinimizationParameter *par = (MinimizationParameter *)params;
00258         int z = (*par).z;
00259         vector<int> FactoryIndex = (*par).FactoryIndex;
00260         vector<int> StockIndex   = (*par).StockIndex;
00261         map< pair<int,int>, IsotopicVector > ParcStock = (*par).ParcStock;
00262         IsotopicVector isotopicvector = (*par).isotopicvector;
00263         IsotopicVector TotalStock = (*par).TotalStock;
00264 
00265 
00266         double x[FactoryIndex.size()];  
00267         
00268         for (int i =0; i < (int)FactoryIndex.size(); i++)
00269         {
00270                 x[i] = gsl_vector_get(v, i);
00271                 if(x[i]<0) x[i] = -x[i];
00272         }
00273         for (int i = 0; i < (int)FactoryIndex.size(); i++)
00274         {
00275                 double Sum = 0;
00276                 map< pair<int,int>, IsotopicVector >::iterator it_i;
00277                 it_i = ParcStock.find(pair<int,int>(FactoryIndex[i], StockIndex[i] ));
00278                 map<ZAI ,double > isotopicquantity = TotalStock.GetAtomicComposition(z).GetIsotopicQuantity();
00279                 for(map<ZAI ,double >::iterator it = isotopicquantity.begin(); it != isotopicquantity.end(); it++)
00280                 {       
00281                         Sum += 2*(*it_i).second.GetAtomicComposition(z).GetZAIIsotopicQuantity((*it).first)
00282                         *isotopicvector.GetAtomicComposition(z).GetZAIIsotopicQuantity((*it).first);
00283                         
00284                         for (int j = 0; j < (int)FactoryIndex.size(); j++) 
00285                         {
00286                                 //if(i == j) j++;
00287 
00288                                 map< pair<int,int>, IsotopicVector >::iterator it_j;
00289                                 it_j = ParcStock.find(pair<int,int>(FactoryIndex[j], StockIndex[j] ));
00290                                 
00291                                 Sum -= 2*x[j]*(*it_i).second.GetZAIIsotopicQuantity((*it).first)
00292                                 *(*it_j).second.GetZAIIsotopicQuantity((*it).first);
00293                         }
00294                 }
00295                 gsl_vector_set(df, i, Sum);
00296                 
00297         }
00298 }
00299 
00300 // Compute both f and df together. 
00301 void my_fdf (const gsl_vector *x, void *params, double *f, gsl_vector *df) 
00302 {
00303         *f = my_f(x, params); 
00304         my_df(x, params, df);
00305 }
00306 
00307 
00308 
00309 IsotopicVector CLASS::MinimizationBuild(IsotopicVector isotopicvector )
00310 {
00311 DBGL;
00312         IsotopicVector BuildIV;
00313         vector<int>             FactoryIndex;
00314         vector<int>             StockIndex;
00315         
00316         map< pair<int,int>, IsotopicVector >::iterator  it;
00317 
00318         for(it = fParcStock.begin(); it != fParcStock.end(); it++ )
00319         {
00320                         FactoryIndex.push_back((*it).first.first);
00321                         StockIndex.push_back((*it).first.second);
00322         }
00323         
00324         
00325         if( (int)FactoryIndex.size() <= 3) return BuildIV;
00326         
00327         for(int k = 0; k < (int)isotopicvector.GetAtomicSpecies().size(); k++ ) // Loop on the Atomic Species
00328         {
00329                 int z = isotopicvector.GetAtomicSpecies()[k];           // Get the Atomic Species Number
00330                 IsotopicVector BuildIVz;
00331                 IsotopicVector IVAtomicCompo = isotopicvector.GetAtomicComposition(z);
00332                 //double delta = 1e-3*Norme(IVAtomicCompo);
00333 
00334                 
00335                 
00336                 
00337                                 
00338                 size_t iter = 0;
00339                 int status;
00340                 
00341                 const gsl_multimin_fdfminimizer_type *T;
00342                 gsl_multimin_fdfminimizer *s;
00343                 
00344                 // Index of IV 
00345                 MinimizationParameter par ;
00346                 par.z = z;
00347                 par.FactoryIndex = FactoryIndex;
00348                 par.StockIndex   = StockIndex;
00349                 par.ParcStock    = fParcStock;
00350                 par.isotopicvector = isotopicvector;
00351                 par.TotalStock = fTotalStock;
00352                 //double par[5] = { 1.0, 2.0, 10.0, 20.0, 30.0 };
00353 
00354                 gsl_vector *x;
00355                 gsl_multimin_function_fdf my_func;
00356 
00357                 my_func.n = (int)FactoryIndex.size();  // number of function components 
00358                 my_func.f = my_f;
00359                 my_func.df = my_df;
00360                 my_func.fdf = my_fdf;
00361                 my_func.params = &par;
00362                 
00363                 // Starting point, x = (5,7) 
00364                 x = gsl_vector_alloc (FactoryIndex.size());
00365                 for (int i = 0; i < (int)FactoryIndex.size(); i++) 
00366                 {
00367                         map< pair<int,int>, IsotopicVector >::iterator  it;
00368                         it = fParcStock.find(pair<int,int>(FactoryIndex[i], StockIndex[i] ));
00369                         gsl_vector_set (x, i, Norme((*it).second.GetAtomicComposition(z))/2. );
00370                 }
00371                 T = gsl_multimin_fdfminimizer_conjugate_fr;
00372                 s = gsl_multimin_fdfminimizer_alloc (T, FactoryIndex.size());
00373                 
00374                 gsl_multimin_fdfminimizer_set (s, &my_func, x, 0.01, 1e-1);
00375                 
00376                 do
00377                 {
00378                         iter++;
00379                         status = gsl_multimin_fdfminimizer_iterate (s);
00380                         
00381                         if (status)
00382                                 break;
00383                         
00384                         status = gsl_multimin_test_gradient (s->gradient, 1e-3);
00385                         //status = gsl_multimin_test_size (0, 1e-3);
00386                         if (status == GSL_SUCCESS)
00387                                 printf ("Minimum found !\n");
00388                         
00389                 }
00390                 while (status == GSL_CONTINUE && iter < 100);
00391                 
00392                 
00393                 for (int i = 0; i < (int)FactoryIndex.size(); i++)
00394                 {
00395                         map< pair<int,int>, IsotopicVector >::iterator it;
00396                         it = fParcStock.find(pair<int,int>(FactoryIndex[i], StockIndex[i] ));
00397                         double Contribution_i = gsl_vector_get (s->x, i);
00398                         BuildIVz += (*it).second.GetAtomicComposition(z)*Contribution_i/Norme((*it).second.GetAtomicComposition(z));
00399                         (*it).second -= (*it).second.GetAtomicComposition(z)*Contribution_i/Norme((*it).second.GetAtomicComposition(z));
00400                 }
00401                 
00402                 
00403                 gsl_multimin_fdfminimizer_free (s);
00404                 gsl_vector_free (x);
00405                 
00406                 BuildIV += BuildIVz;
00407         }
00408         return BuildIV;
00409 DBGL;   
00410 }
00411 */
00412 
00413 //________________________________________________________________________
00414 IsotopicVector CLASS::MonteCarloBuild(IsotopicVector isotopicvector )
00415 {
00416 
00417         IsotopicVector BuildIV;
00418         
00419         
00420         int z = 0;
00421         
00422         vector<int>             FactoryIndex;
00423         vector<int>             StockIndex;
00424         map< pair<int,int>, IsotopicVector >::iterator  it;
00425 
00426         for(it = fParcStock.begin(); it != fParcStock.end(); it++ )
00427         {
00428                         FactoryIndex.push_back((*it).first.first);
00429                         StockIndex.push_back((*it).first.second);
00430         }
00431         
00432         if( (int)FactoryIndex.size() == 0) return BuildIV;
00433 
00434         for(int k = 0; k < (int)isotopicvector.GetAtomicSpecies().size(); k++ ) // Loop on the Atomic Species
00435         {
00436                 z = isotopicvector.GetAtomicSpecies()[k];               // Get the Atomic Species Number
00437                 
00438                 IsotopicVector IVAtomicCompo = isotopicvector.GetAtomicComposition(z);
00439                 IsotopicVector BuildIVz;
00440                 double delta = 1e-3*Norme(IVAtomicCompo);
00441                 
00442                 
00443                 
00444                 vector<int>             FactoryIndextmp = FactoryIndex;
00445                 vector<int>             StockIndextmp   = StockIndex;
00446                 
00447                 double DistanceN = Norme(IVAtomicCompo);        // Distance at the Nst iteration
00448                 double DistanceN_1 = Norme(IVAtomicCompo);      // Distance at the N-1 st iteration
00449                 
00450                 while(FactoryIndextmp.size() !=0)
00451                 {
00452                         int IVIndex = FactoryIndextmp.size();
00453                         
00454                         // Find a random IV contribution
00455                         while(IVIndex == (int)FactoryIndextmp.size()) IVIndex = (int)random(0,FactoryIndextmp.size());
00456                         
00457                         map< pair<int,int>, IsotopicVector >::iterator it;
00458                         it = fParcStock.find(pair<int,int>(FactoryIndextmp[IVIndex], StockIndextmp[IVIndex] )); 
00459                         IsotopicVector IVToAdd = (*it).second.GetAtomicComposition(z) ; 
00460                         
00461                         
00462                         IsotopicVector BuildIVzTmp = BuildIVz + IVToAdd/Norme(IVToAdd) * delta; // Build the new IVBuild 
00463                         
00464                         DistanceN = Distance( BuildIVzTmp, isotopicvector.GetAtomicComposition(z) ); //Progress determination
00465                         
00466                         if(DistanceN > DistanceN_1 || delta > Norme(IVToAdd) ) // Not Good Peak remove it from Index
00467                         {
00468                                 FactoryIndextmp.erase(FactoryIndextmp.begin()+IVIndex);
00469                                 StockIndextmp.erase(StockIndextmp.begin()+IVIndex);
00470                         }
00471                         else    // Good Peak, Add the contribution to IVBuild and remove from stock...
00472                         {
00473                                 // Validation of Values....
00474                                 DistanceN_1 = DistanceN;
00475                                 BuildIVz = BuildIVzTmp;
00476                                 
00477                                 // remove contribution from stock
00478                                 map< pair<int,int>, IsotopicVector >::iterator it;
00479                                 it = fParcStock.find(pair<int,int>(FactoryIndextmp[IVIndex], StockIndextmp[IVIndex] ));
00480                                 
00481                                 (*it).second -= IVToAdd/Norme(IVToAdd) * delta;
00482                         }
00483                 }
00484                 BuildIV += BuildIVz;
00485         }
00486         
00487 DBGL;
00488         return BuildIV;
00489 }
00490 
00491 
00492 //________________________________________________________________________
00493 void CLASS::BuildTimeVector(long int t)
00494 {
00495 DBGL;
00496         fTimeStep.clear();
00497         fTimeStep.insert( pair<long int ,int>(t,1) );
00498 
00499 //********* Printing Step *********//
00500         {
00501                 long int step = 0;
00502                 if(step >= fAbsoluteTime        )
00503                         fTimeStep.insert( pair<long int ,int>(step,1) );
00504                 step += fPrintStep;
00505                 do 
00506                 {
00507                         
00508                         if(step >= fAbsoluteTime        )
00509                                 fTimeStep.insert( pair<long int ,int>(step,1) );
00510                         step += fPrintStep;
00511                 }
00512                 while( step < t );
00513         }
00514 
00515         for(int i = 0; i < (int)fReactor.size();i++)
00516         {
00517                 long int step = fReactor[i]->GetCreationTime();
00518                 long int coolingstep = fReactor[i]->GetAssociedTreatmentFactory()->GetCoolingTime();
00519                 long int seprationstep = fReactor[i]->GetAssociedTreatmentFactory()->GetSeparationTime();
00520 
00521 //********* Reactor Evolution Step *********//
00522                 // set destruction of a reactor
00523                 if( (fReactor[i]->GetCreationTime() + fReactor[i]->GetLifeTime() > fAbsoluteTime) &&
00524                     (fReactor[i]->GetCreationTime() + fReactor[i]->GetLifeTime()<t) )
00525                 {
00526                         pair< map<long int, int>::iterator, bool > IResult ;
00527                         IResult = fTimeStep.insert( pair<long int ,int>(fReactor[i]->GetCreationTime() + fReactor[i]->GetLifeTime(),2) );
00528                         if( IResult.second == false ) IResult.first->second |= 2;
00529                 }
00530                 
00531                 // Set End of reactor cycle
00532                 if(step >= fAbsoluteTime &&  step <= t && step < fReactor[i]->GetCreationTime() + fReactor[i]->GetLifeTime())           
00533                 {
00534                         pair< map<long int, int>::iterator, bool > IResult = fTimeStep.insert( pair<long int ,int>(step,4) );
00535                         if( IResult.second == false ) IResult.first->second |= 4;
00536                 }
00537                 
00538 
00539 //********* TF Evolution Step *********//
00540                 step += fReactor[i]->GetCycleTime();
00541                 do 
00542                 {
00543                         
00544                         if(step > fAbsoluteTime && step <= t && step < fReactor[i]->GetCreationTime() + fReactor[i]->GetLifeTime())
00545                         {                                               // Set End of reactor cycle
00546                                 pair< map<long int, int>::iterator, bool > IResult = fTimeStep.insert( pair<long int ,int>(step,4) );
00547                                 if( IResult.second == false ) IResult.first->second  |= 4;
00548                         }
00549 
00550                         if(step >= fAbsoluteTime && step + coolingstep <= t)                    // Set End of Cooling
00551                         {
00552                                 pair< map<long int, int>::iterator, bool > IResult = fTimeStep.insert( pair<long int ,int>(step+coolingstep,8) );
00553                                 if( IResult.second == false ) IResult.first->second |= 8;
00554                         }
00555 
00556                         if(step >= fAbsoluteTime && step + coolingstep + seprationstep <= t)    // Set end of Separation
00557                         {
00558                                 pair< map<long int, int>::iterator, bool > IResult;
00559                                 IResult = fTimeStep.insert( pair<long int ,int>(step+coolingstep+seprationstep,16) );
00560                                 if( IResult.second == false ) IResult.first->second |= 16;
00561                         }
00562                         step += fReactor[i]->GetCycleTime();
00563                 }
00564                 while(step <= t && step <= fReactor[i]->GetCreationTime() + fReactor[i]->GetLifeTime() );
00565         }
00566         
00567 //********* TF Evolution Step *********//
00568 //*** In Case of Evolution Restart ****//
00569         for(int i =0; i < (int)fTreatmentFactory.size(); i++) 
00570         {
00571                 
00572                 for(int j = 0; j<(int)fTreatmentFactory[i]->GetIVCooling().size(); j++ )// Set End of Cooling
00573                 {
00574                         if(fTreatmentFactory[i]->GetCoolingStartingTime()[j] +  fTreatmentFactory[i]->GetCoolingTime() > fAbsoluteTime )
00575                         {
00576                                 pair< map<long int, int>::iterator, bool > IResult;
00577                                 IResult = fTimeStep.insert( pair<long int ,int>(fTreatmentFactory[i]->GetCoolingStartingTime()[j] +  fTreatmentFactory[i]->GetCoolingTime(),8) );
00578                                 if( IResult.second == false ) IResult.first->second |= 8;
00579                         }
00580                 }
00581                 for(int j = 0; j< (int)fTreatmentFactory[i]->GetIVSeparating().size(); j++ )// Set end of Separation
00582                 {
00583                         if(fTreatmentFactory[i]->GetSeparatingStartingTime()[j] +  fTreatmentFactory[i]->GetSeparationTime() > fAbsoluteTime )
00584                         {
00585                                 pair< map<long int, int>::iterator, bool > IResult = fTimeStep.insert( pair<long int ,int>(fTreatmentFactory[i]->GetSeparatingStartingTime()[j] +  fTreatmentFactory[i]->GetSeparationTime(),16) );
00586                                 if( IResult.second == false ) IResult.first->second |= 16;
00587                         }
00588                 }
00589         }
00590 
00591         
00592         
00593 //****** Print the Time Index ******//
00594         ofstream TimeStepfile("CLASS_TimeStep", ios_base::app);         // Open the File
00595         
00596         if(!TimeStepfile)
00597         {
00598                 cout            << "!!Warning!! !!!CLASS!!! Can't open \" CLASS_TimeStep \"\n" << endl;
00599                 fLog->fLog      << "!!Warning!! !!!CLASS!!! Can't open \" CLASS_TimeStep \"\n" << endl;
00600         }
00601         for(map<long int ,int >::iterator it = fTimeStep.begin(); it != fTimeStep.end(); it++)
00602                 TimeStepfile << (*it).first/365.4/3600/24 << " " << (*it).second << endl;
00603         
00604         DBGL;   
00605 }
00606 
00607 
00608 
00609 //________________________________________________________________________
00610 //___________________________ Evolution Method ___________________________
00611 //________________________________________________________________________
00612 
00613 void CLASS::TreatmentEvolution()
00614 {
00615 DBGL;
00616         for(int i = 0; i < (int) fTreatmentFactory.size();i++)
00617                 fTreatmentFactory[i]->Evolution(fAbsoluteTime);
00618 DBGL;
00619 }
00620 
00621 //________________________________________________________________________
00622 void CLASS::ReactorEvolution()
00623 {
00624 DBGL;
00625         
00626 
00627  #pragma omp sections
00628         {
00629         #pragma omp section
00630                 {TreatmentEvolution();}
00631 
00632         #pragma omp section
00633         {
00634                 #pragma omp parallel for
00635                         for(int i = 0; i < (int)fReactor.size(); i++)
00636                                 fReactor[i]->Evolution(fAbsoluteTime);
00637         }
00638         }
00639         
00640 
00641         for(int i = 0; i < (int) fTreatmentFactory.size();i++)
00642                 fTreatmentFactory[i]->Dump();
00643 
00644         UpdateParcStock();
00645         for(int i = 0; i < (int)fReactor.size(); i++)
00646                 fReactor[i]->Dump();
00647         
00648         DumpParcStock();
00649         
00650 DBGL;
00651 }
00652 
00653 //________________________________________________________________________
00654 void CLASS::Evolution(long int t)
00655 {
00656 DBGL;
00657         OpenOutputTree();
00658         OutAttach();
00659         BuildTimeVector(t);
00660 
00661         int Start = time(NULL);
00662         for(map<long int ,int >::iterator it = fTimeStep.begin(); it != fTimeStep.end(); it++)
00663         {
00664 
00665                 ResetQuantity();
00666                 fAbsoluteTime = (*it).first;
00667                 //if(fAbsoluteTime > t) return;
00668                 
00669                 if( (*it).second & 2 || (*it).second & 1 )
00670                 {
00671                         ReactorEvolution();
00672                         if((*it).second & 2 )
00673                                 RemoveReactor();
00674                         
00675                         if((*it).second & 2 )
00676                                 (*it).second ^= 2;
00677                         if((*it).second & 4 )
00678                                 (*it).second ^= 4;
00679                         if((*it).second & 8 )
00680                                 (*it).second ^= 8;
00681                         if((*it).second & 16 )
00682                                 (*it).second ^= 16;
00683                 }
00684                 
00685                 if( (*it).second & 4 )
00686                 {
00687                         ReactorEvolution();
00688                         if((*it).second & 4 )
00689                                 (*it).second ^= 4;
00690                         if((*it).second & 8 )
00691                                 (*it).second ^= 8;
00692                         if((*it).second & 16 )
00693                                 (*it).second ^= 16;
00694                 }
00695                 
00696                 if( (*it).second & 8 || (*it).second & 16 )
00697                 {
00698                         TreatmentEvolution();
00699                         for(int i = 0; i < (int) fTreatmentFactory.size();i++)
00700                                 fTreatmentFactory[i]->Dump();
00701 
00702                         if((*it).second & 8 )
00703                                 (*it).second ^= 8;
00704                         if((*it).second & 16 )
00705                                 (*it).second ^= 16;
00706                 }
00707                 
00708 
00709                 
00710                 if( (*it).second & 1 )
00711                         #pragma omp single
00712                         {
00713 //                              Print();
00714 //                              Write();
00715                                 fOutT->Fill();
00716                                 //fTotalStock.Print();
00717                                 ProgressPrintout(Start, t);
00718                         }
00719         }
00720 #pragma omp single
00721         {CloseOutputTree();}
00722 DBGL;
00723 }
00724 void CLASS::ProgressPrintout(int starttime, long int t)
00725 {
00726 DBGL;
00727         int TimeNow = time(NULL);
00728         int Spent = difftime(TimeNow, starttime);
00729         double Time = fAbsoluteTime/3600/24/365.4;
00730         double Total = t/3600/24/365.4;
00731         double Remain =  (Spent/Time * Total - Spent);
00732         cout << "                                                                                             " << flush ;
00733         cout << "\rProcessed " << setprecision(4) << Time << " / " << setprecision(4) << Total << " STEP "
00734              << setprecision(4) << (Time/Total*100.0) << "%) ---  I still need : "
00735         << (long int)Remain/60 << " min "<<  (long int)(Remain)%60 << " sec to finish ! \r" << flush;
00736 DBGL;
00737 }
00738 
00739 //________________________________________________________________________
00740 //______________________________ Out Method ______________________________
00741 //________________________________________________________________________
00742 
00743 void CLASS::ResetQuantity()
00744 {
00745 DBGL;
00746         IsotopicVector EmptyIV;
00747         fTotalInReactor         = EmptyIV;
00748         fTotalWaste             = EmptyIV;
00749         fTotalStock             = EmptyIV;
00750 //      fTotalGodIncome         = EmptyIV;
00751         fTotalCooling           = EmptyIV;
00752         fTotalSeparating        = EmptyIV;
00753         fIVInCycleTotal         = EmptyIV;
00754         fIVTotal                = EmptyIV;
00755 DBGL;
00756 }
00757 
00758 //________________________________________________________________________
00759 void CLASS::OpenOutputTree()
00760 {
00761 DBGL;
00762 
00763         cout << "Opening : " << fOutputName << " ...";
00764         fLog->fLog << "Opening : " << fOutputName << " ...";
00765         fOutTree = new TFile(fOutputName.c_str(),"UPDATE");
00766 
00767         if(!fOutTree)
00768         {
00769                 cout << "\nCould not open " << fOutputName <<endl;
00770                 fLog->fLog << "\nCould not open " << fOutputName <<endl;
00771                 exit(-1);
00772         }
00773         cout << "\t ...O.K." << endl;
00774 
00775         fOutT = new TTree("Data","Data Tree");
00776         cout << "Creating Data Tree...";
00777         fLog->fLog << "Creating Data Tree...";
00778         if(!fOutTree)
00779         {
00780                 cout << "\nCould not create Data Tree in " << fOutputName << endl;
00781                 fLog->fLog << "\nCould not create Data Tree in " << fOutputName << endl;
00782                 exit(-1);
00783         }
00784         cout << "\t ...O.K!" << endl;
00785         fLog->fLog <<  "\t ...O.K!" << endl;
00786         // fOutT->SetAutoFlush(1000000);
00787         // fOutT->SetAutoSave(100000);
00788 }
00789 void CLASS::CloseOutputTree()
00790 {
00791 DBGL;
00792 
00793         fOutTree->ls();
00794         cout << "Writing outTree " << fOutputName << endl;
00795         fLog->fLog << "Writing outTree " << fOutputName << endl;
00796         fOutTree->Write();
00797 
00798         if(fOutTree->IsOpen()) {
00799                 cout << "Deleting outTree : " << endl;
00800                 fLog->fLog << "Deleting outTree : " << endl;
00801                 delete fOutT;
00802                 cout << "Closing file : " << fOutputName <<endl;
00803                 fLog->fLog << "Closing file : " << fOutputName <<endl;
00804                 fOutTree-> Close();
00805                 delete fOutTree;
00806         } else {
00807                 cout << "File was not opened " << fOutputName << endl;
00808                 fLog->fLog << "File was not opened " << fOutputName << endl;
00809                 exit(-1);
00810         }
00811 }
00812 //________________________________________________________________________
00813 void CLASS::OutAttach()
00814 {
00815 DBGL;
00816         ResetQuantity();
00817         //Branch Absolut Time
00818         fOutT->Branch("AbsTime",&fAbsoluteTime,"AbsoluteTime/L");
00819         
00820         // Branch the Sum IV
00821 
00822         
00823         fOutT->Branch("TF_WASTE.", "IsotopicVector", &fTotalWaste);
00824         fOutT->Branch("TF_STOCK.", "IsotopicVector", &fTotalStock);
00825         fOutT->Branch("TF_SEPARATING.", "IsotopicVector", &fTotalSeparating);
00826         fOutT->Branch("TF_COOLING.", "IsotopicVector", &fTotalCooling);
00827         fOutT->Branch("GODINCOME.", "IsotopicVector", &fTotalGodIncome);
00828         fOutT->Branch("INTOTALINREACTOR.", "IsotopicVector", &fTotalInReactor);
00829         fOutT->Branch("INCYCLE.", "IsotopicVector", &fIVInCycleTotal);
00830         fOutT->Branch("TOTAL.", "IsotopicVector", &fIVTotal);
00831         
00832         fOutT->Branch("TreatmentFactory", "TreatmentFactory", &fTreatmentFactory.at(0));
00833 //      fOutT->Branch("Reactor", "Reactor", &fReactor.at(0));
00834         
00835 DBGL;
00836 }
00837 
00838 //________________________________________________________________________
00839 void CLASS::Write()
00840 {       
00841 DBGL;
00842         string basename = "CLASS";
00843 
00844         //****** Total *****//
00845         IsotopicVector TotalInReactor;
00846         for(int i = 0; i < (int)fReactor.size(); i++)
00847         {       
00848                 TotalInReactor += fReactor[i]->GetIVReactor();
00849         }
00850         string RTotal = basename + "_R_TOTAL";
00851         TotalInReactor.Write(RTotal, fAbsoluteTime);
00852 
00853         //****** Total *****//
00854         IsotopicVector TotalWaste;
00855         IsotopicVector TotalGodIncome;
00856         IsotopicVector TotalStock;
00857         IsotopicVector TotalCooling;
00858         IsotopicVector TotalSeparating;
00859 
00860 
00861         for(int i = 0; i < (int) fTreatmentFactory.size();i++)
00862         {
00863                 TotalWaste += fTreatmentFactory[i]->GetIVWaste();
00864                 TotalStock += fTreatmentFactory[i]->GetIVFullStock();
00865                 TotalGodIncome += fTreatmentFactory[i]->GetIVGodIncome();
00866                 
00867                 for(int j=0; j < (int)fTreatmentFactory[i]->GetIVCooling().size(); j++)
00868                         TotalCooling += fTreatmentFactory[i]->GetIVCooling()[j];
00869                 for(int j=0; j < (int)fTreatmentFactory[i]->GetIVSeparating().size(); j++)
00870                         TotalSeparating += fTreatmentFactory[i]->GetIVSeparating()[j];
00871         }
00872 
00873         string TFWaste = basename + "_TF_WASTE";
00874         string TFStock = basename + "_TF_STOCK";
00875         string TFGodIncome  = basename + "_TF_GODINCOME";
00876         string TFSeparating = basename + "_TF_SPERATING";
00877 
00878         TotalWaste.Write(TFWaste, fAbsoluteTime);
00879         TotalStock.Write(TFStock, fAbsoluteTime);
00880         TotalGodIncome.Write(TFGodIncome, fAbsoluteTime);
00881         TotalSeparating.Write(TFSeparating, fAbsoluteTime);
00882 
00883 //---- Total Writing ----/
00884 
00885         string InCycleTotal     = basename + "_IN_CYCLE_TOTAL";
00886         string Total            = basename + "_TOTAL";
00887 
00888         IsotopicVector IVInCycleTotal;
00889         IsotopicVector IVTotal;
00890 
00891         IVInCycleTotal  = TotalCooling + TotalInReactor + TotalSeparating + TotalStock;
00892         IVTotal         = IVInCycleTotal + TotalWaste;
00893 
00894         IVInCycleTotal.Write(InCycleTotal, fAbsoluteTime);
00895         IVTotal.Write(Total, fAbsoluteTime);
00896 
00897 DBGL;
00898 }
00899 
00900 //________________________________________________________________________
00901 void CLASS::Print()
00902 {
00903 DBGL;
00904         for(int i = 0; i < (int) fTreatmentFactory.size();i++)
00905         {
00906                 cout << "!!!!!!!!!STEP : " << fAbsoluteTime/(int)(3600*24*365.4) << endl;
00907                 cout << "TreatmentFactory : " << endl;
00908                 cout << "Cooling ";
00909                 cout << fTreatmentFactory[i]->GetIVCooling().size()<< endl;
00910                 cout << "Waste " << endl;
00911                 fTreatmentFactory[i]->GetIVWaste().Print();
00912                 cout << "Stock :" << endl;
00913                 for(int j=0; j < (int)fTreatmentFactory[i]->GetIVStock().size(); j++)
00914                 {
00915                 cout << j << endl;
00916                 fTreatmentFactory[i]->GetIVStock()[j].Print();
00917                 }
00918         }
00919 
00920         for(int i = 0; i < (int)fReactor.size(); i++)
00921         {       
00922                 cout << "Reactor" << endl;
00923                 fReactor[i]->GetIVReactor().Print();
00924         }
00925 DBGL;
00926 }
 Tout Classes Fichiers Fonctions Variables Macros