CLASS
1.1
|
00001 #include "IsotopicVector.hxx" 00002 #include "ZAI.hxx" 00003 #include "LogFile.hxx" 00004 #include "Defines.hxx" 00005 00006 00007 #include <cmath> 00008 #include <iostream> 00009 #include <fstream> 00010 #include <string> 00011 #include <algorithm> 00012 00013 //________________________________________________________________________ 00014 //________________________________________________________________________ 00015 // 00016 // 00017 // 00018 // IsotopicVector 00019 // 00020 // 00021 //________________________________________________________________________ 00022 //________________________________________________________________________ 00023 00024 00025 00026 00027 //________________________________________________________________________ 00028 //__________________________Operator Overlaoding__________________________ 00029 //________________________________________________________________________ 00030 00031 //____________________________General Operator____________________________ 00032 //________________________________________________________________________ 00033 00034 ClassImp(IsotopicVector) 00035 00036 double Norme(IsotopicVector IV1) 00037 { 00038 DBGL; 00039 IsotopicVector IV; 00040 return Distance(IV1, IV); 00041 DBGL; 00042 } 00043 00044 double Distance(IsotopicVector IV1, IsotopicVector IV2 ) 00045 { 00046 DBGL; 00047 double d2 = 0; 00048 00049 IsotopicVector IVtmp = IV1 + IV2; 00050 map<ZAI ,double> IVtmpIsotopicQuantity = IVtmp.GetIsotopicQuantity(); 00051 00052 for(map<ZAI ,double >::iterator it = IVtmpIsotopicQuantity.begin(); it != IVtmpIsotopicQuantity.end(); it++) 00053 { 00054 double Z1 = IV1.GetZAIIsotopicQuantity((*it).first); 00055 double Z2 = IV2.GetZAIIsotopicQuantity((*it).first); 00056 d2 += pow(Z1-Z2 , 2 ); 00057 } 00058 00059 DBGL; 00060 return sqrt(d2); 00061 } 00062 00063 00064 IsotopicVector operator+(IsotopicVector const& IVa, IsotopicVector const& IVb) 00065 { 00066 IsotopicVector IVtmp; 00067 IVtmp = IVa; 00068 return IVtmp += IVb; 00069 } 00070 00071 //________________________________________________________________________ 00072 IsotopicVector operator-(IsotopicVector const& IVa, IsotopicVector const& IVb) 00073 { 00074 IsotopicVector IVtmp; 00075 IVtmp = IVa; 00076 return IVtmp -= IVb; 00077 } 00078 00079 00080 //________________________________________________________________________ 00081 IsotopicVector operator*(ZAI const& zai, double F) 00082 { 00083 IsotopicVector IVtmp; 00084 00085 IVtmp.Add( zai, F); 00086 00087 return IVtmp; 00088 } 00089 //________________________________________________________________________ 00090 IsotopicVector operator/(ZAI const& zai, double F) 00091 { 00092 IsotopicVector IVtmp; 00093 00094 IVtmp.Add( zai, 1./F); 00095 00096 return IVtmp; 00097 } 00098 00099 00100 //________________________________________________________________________ 00101 IsotopicVector operator*(IsotopicVector const& IVA, double F) 00102 { 00103 DBGL; 00104 IsotopicVector IV = IVA; 00105 IV.Multiply(F); 00106 return IV; 00107 } 00108 00109 //________________________________________________________________________ 00110 IsotopicVector operator/(IsotopicVector const& IVA, double F) 00111 { 00112 DBGL; 00113 IsotopicVector IV = IVA; 00114 IV.Multiply(1./F); 00115 return IV; 00116 } 00117 00118 00119 //____________________________InClass Operator____________________________ 00120 00121 //________________________________________________________________________ 00122 IsotopicVector& IsotopicVector::operator+=(const IsotopicVector& IVa) 00123 { 00124 Add(IVa); 00125 return *this; 00126 } 00127 00128 //________________________________________________________________________ 00129 IsotopicVector& IsotopicVector::operator-=(const IsotopicVector& IVa) 00130 { 00131 Remove(IVa); 00132 return *this; 00133 } 00134 00135 00136 00137 //________________________________________________________________________ 00138 //________________________Constructor & Destructor________________________ 00139 //________________________________________________________________________ 00140 IsotopicVector::IsotopicVector() 00141 { 00142 DBGL; 00143 DBGL; 00144 } 00145 00146 00147 //________________________________________________________________________ 00148 IsotopicVector::~IsotopicVector() 00149 { 00150 DBGL; 00151 } 00152 00153 00154 00155 //________________________________________________________________________ 00156 //_____________________________General Method_____________________________ 00157 //________________________________________________________________________ 00158 void IsotopicVector::Clear() 00159 { 00160 DBGL; 00161 fIsotopicQuantityNeeded.clear(); 00162 fIsotopicQuantity.clear(); 00163 DBGL; 00164 } 00165 //________________________________________________________________________ 00166 void IsotopicVector::ClearNeed() 00167 { 00168 DBGL; 00169 fIsotopicQuantityNeeded.clear(); 00170 DBGL; 00171 } 00172 00173 //________________________________________________________________________ 00174 void IsotopicVector::Multiply(double factor) 00175 { 00176 DBGL; 00177 for(map<ZAI ,double >::iterator it = fIsotopicQuantity.begin(); it != fIsotopicQuantity.end(); it++) 00178 (*it).second = (*it).second * factor; 00179 DBGL; 00180 00181 } 00182 00183 //________________________________________________________________________ 00184 void IsotopicVector::Add(const ZAI& zai, double quantity) 00185 { 00186 DBGL; 00187 if( ceil(quantity*1e6) - quantity*1e6 > quantity*1e6 - floor(quantity*1e6) ) 00188 quantity = floor(quantity*1e6)*1/1e6; 00189 else quantity = ceil(quantity*1e6)*1/1e6; 00190 00191 pair<map<ZAI, double>::iterator, bool> IResult; 00192 if(quantity > 0) 00193 { 00194 IResult = fIsotopicQuantity.insert( pair<ZAI ,double>(zai, quantity)); 00195 if(IResult.second == false) 00196 IResult.first->second += quantity; 00197 } 00198 00199 DBGL; 00200 } 00201 //________________________________________________________________________ 00202 00203 void IsotopicVector::Add(const IsotopicVector& isotopicvector) 00204 { 00205 DBGL; 00206 map<ZAI ,double> isotopicquantity = isotopicvector.GetIsotopicQuantity(); 00207 00208 for(map<ZAI ,double >::iterator it = isotopicquantity.begin(); it != isotopicquantity.end(); it++) 00209 Add( (*it).first, (*it).second); 00210 DBGL; 00211 00212 } 00213 //________________________________________________________________________ 00214 00215 void IsotopicVector::Add(const map<ZAI ,double>& quantity) 00216 { 00217 DBGL; 00218 map<ZAI ,double> isotopicquantity = quantity; 00219 00220 for(map<ZAI ,double >::iterator it = isotopicquantity.begin(); it != isotopicquantity.end(); it++) 00221 Add( (*it).first, (*it).second); 00222 DBGL; 00223 00224 } 00225 00226 00227 //________________________________________________________________________ 00228 void IsotopicVector::Remove(const ZAI& zai, double quantity) 00229 { 00230 DBGL; 00231 00232 map<ZAI ,double>::iterator it; 00233 it = fIsotopicQuantity.find(zai); 00234 00235 if( ceil(quantity*1e6) - quantity*1e6 > quantity*1e6 - floor(quantity*1e6) ) 00236 quantity = floor(quantity*1e6)*1/1e6; 00237 else quantity = ceil(quantity*1e6)*1/1e6; 00238 00239 if(quantity > 0) 00240 { 00241 if ( it != fIsotopicQuantity.end() ) 00242 { 00243 if (it->second > quantity) 00244 it->second = it->second - quantity; 00245 else 00246 { 00247 Need(zai, quantity - it->second ); 00248 it->second = 0; 00249 } 00250 } 00251 else 00252 { 00253 Need(zai, quantity); 00254 } 00255 } 00256 DBGL; 00257 00258 } 00259 00260 //________________________________________________________________________ 00261 void IsotopicVector::Remove(const IsotopicVector& isotopicvector) 00262 { 00263 DBGL; 00264 map<ZAI ,double> isotopicquantity = isotopicvector.GetIsotopicQuantity(); 00265 for(map<ZAI ,double >::iterator it = isotopicquantity.begin(); it != isotopicquantity.end(); it++) 00266 Remove( (*it).first, (*it).second); 00267 DBGL; 00268 } 00269 00270 //________________________________________________________________________ 00271 void IsotopicVector::Need(const ZAI& zai, double quantity) 00272 { 00273 DBGL; 00274 if( ceil(quantity*1e6) - quantity*1e6 > quantity*1e6 - floor(quantity*1e6) ) 00275 quantity = floor(quantity*1e6)*1/1e6; 00276 else quantity = ceil(quantity*1e6)*1/1e6; 00277 00278 00279 pair<map<ZAI, double>::iterator, bool> IResult; 00280 if(quantity > 0) 00281 { 00282 IResult = fIsotopicQuantityNeeded.insert( pair<ZAI ,double>(zai, quantity)); 00283 if(IResult.second == false) 00284 IResult.first->second += quantity; 00285 } 00286 DBGL; 00287 00288 00289 } 00290 00291 //________________________________________________________________________ 00292 void IsotopicVector::Need(const IsotopicVector& isotopicvector) 00293 { 00294 DBGL; 00295 map<ZAI ,double> isotopicquantity = isotopicvector.GetIsotopicQuantity(); 00296 00297 for(map<ZAI ,double >::iterator it = isotopicquantity.begin(); it != isotopicquantity.end(); it++) 00298 Need( (*it).first, (*it).second); 00299 DBGL; 00300 } 00301 00302 00303 //________________________________________________________________________ 00304 double IsotopicVector::GetZAIIsotopicQuantity(const ZAI& zai) const 00305 { 00306 DBGL; 00307 map<ZAI ,double> IsotopicQuantity = GetIsotopicQuantity(); 00308 00309 map<ZAI ,double>::iterator it; 00310 it = IsotopicQuantity.find(zai); 00311 DBGL; 00312 00313 if ( it != IsotopicQuantity.end() ) 00314 { 00315 return it->second; 00316 } 00317 else 00318 { 00319 return 0; 00320 } 00321 } 00322 00323 //________________________________________________________________________ 00324 double IsotopicVector::GetZAIIsotopicQuantity(const int z, const int a, const int i) const 00325 { 00326 DBGL; 00327 ZAI zai(z, a, i); 00328 return GetZAIIsotopicQuantity(zai); 00329 } 00330 00331 IsotopicVector IsotopicVector::GetAtomicComposition(int z) const 00332 { 00333 DBGL; 00334 IsotopicVector IV; 00335 map<ZAI ,double > IsotopicQuantity = GetIsotopicQuantity(); 00336 map<ZAI ,double >::iterator it; 00337 for( it = IsotopicQuantity.begin(); it != IsotopicQuantity.end(); it++) 00338 00339 if( (*it).first.Z() == z ) 00340 IV += (*it).first * (*it).second; 00341 00342 return IV; 00343 DBGL; 00344 } 00345 00346 vector<int> IsotopicVector::GetAtomicSpecies() const 00347 { 00348 DBGL; 00349 vector<int> AtomicSpecies; 00350 00351 map<ZAI ,double > IsotopicQuantity = GetIsotopicQuantity(); 00352 map<ZAI ,double >::iterator it; 00353 for( it = IsotopicQuantity.begin(); it != IsotopicQuantity.end(); it++) 00354 if( (int)AtomicSpecies.size() ==0 || (*it).first.Z() != AtomicSpecies.back() ) 00355 AtomicSpecies.push_back((*it).first.Z()); 00356 00357 DBGL; 00358 return AtomicSpecies; 00359 } 00360 00361 00362 //________________________________________________________________________ 00363 void IsotopicVector::Write(string filename, long int time) const 00364 { 00365 ofstream IVfile(filename.c_str(), ios_base::app); // Open the File 00366 if(!IVfile) 00367 cout << "!!Warning!! !!!IsotopicVector!!! \n Can't open \"" << filename << "\"\n" << endl; 00368 00369 IVfile << time << " "; 00370 00371 map<ZAI ,double> IsotopicQuantity = GetIsotopicQuantity(); 00372 for(map<ZAI ,double >::iterator it = IsotopicQuantity.begin(); 00373 it != IsotopicQuantity.end(); it++) 00374 { 00375 IVfile << (*it).first.Z() << " "; 00376 IVfile << (*it).first.A() << " "; 00377 IVfile << (*it).first.I() << " "; 00378 IVfile << (*it).second << " "; 00379 00380 if((*it).first.Z()>70) 00381 { 00382 char Z[33]; 00383 sprintf(Z,"%d",(*it).first.Z()); 00384 char A[33]; 00385 sprintf(A,"%d",(*it).first.A()); 00386 char I[33]; 00387 sprintf(I,"%d",(*it).first.I()); 00388 00389 00390 string filenameZAI = filename+ "_DIR" +"/" + Z + "_" + A + "_"+ I; 00391 string cmd = " mkdir " + filename+ "_DIR" +" 2>/dev/null"; 00392 int test = system(cmd.c_str()); 00393 if (test!=0) cout << "!!Warning!! !!!IsotopicVector!!! \n Can't create \"" << filename<< "_DIR directory " << "\"\n" << endl; 00394 ofstream IVfileZAI(filenameZAI.c_str(), ios_base::app); // Open the File 00395 00396 if(!IVfileZAI) 00397 cout << "!!Warning!! !!!IsotopicVector!!! \n Can't open \"" << filenameZAI << "\"\n" << endl; 00398 00399 IVfileZAI << time << " " << (*it).second << endl;; 00400 IVfileZAI.close(); 00401 } 00402 00403 } 00404 IVfile << endl; 00405 } 00406 //________________________________________________________________________ 00407 void IsotopicVector::Print(string option) const 00408 { 00409 DBGL; 00410 cout << "**************************" << endl; 00411 cout << "*Isotopic Vector Property*" << endl; 00412 cout << "**************************" << endl << endl; 00413 00414 bool QuantityPrint = false; 00415 bool DBPrint = false; 00416 00417 QuantityPrint = true; 00418 00419 if(QuantityPrint) 00420 { 00421 cout << "*Isotopic Vector Quantity*" << endl; 00422 map<ZAI ,double> IsotopicQuantity = GetIsotopicQuantity(); 00423 for(map<ZAI ,double >::iterator it = IsotopicQuantity.begin(); 00424 it != IsotopicQuantity.end(); it++) 00425 { 00426 cout << (*it).first.Z() << " "; 00427 cout << (*it).first.A() << " "; 00428 cout << (*it).first.I() << " "; 00429 cout << ": " << (*it).second; 00430 cout << endl; 00431 } 00432 cout << endl; 00433 cout << "*Isotopic Vector Quantity Needed*" << endl; 00434 map<ZAI ,double> IsotopicQuantityNeeded = GetIsotopicQuantityNeeded(); 00435 for(map<ZAI ,double >::iterator it = IsotopicQuantityNeeded.begin(); 00436 it != IsotopicQuantityNeeded.end(); it++) 00437 { 00438 cout << (*it).first.Z() << " "; 00439 cout << (*it).first.A() << " "; 00440 cout << (*it).first.I() << " "; 00441 cout << ": " << (*it).second; 00442 cout << endl; 00443 } 00444 cout << endl; 00445 } 00446 if(DBPrint) 00447 { 00448 cout << "****Isotopic Vector DB****" << endl; 00449 } 00450 DBGL; 00451 } 00452 00453 00454 00455 00456 00457