00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #include "G4MoleculeCounter.hh"
00029 #include "G4UIcommand.hh"
00030 #include "G4UnitsTable.hh"
00031 #include <iomanip>
00032
00033 double compDoubleWithPrecision::fPrecision = 5e-13;
00034
00035 G4MoleculeCounter::G4MoleculeCounter()
00036 {
00037 fUse = FALSE;
00038 fVerbose = 0 ;
00039 }
00040
00041 G4MoleculeCounter* G4MoleculeCounter::fpInstance = 0;
00042
00043 G4MoleculeCounter* G4MoleculeCounter::GetMoleculeCounter()
00044 {
00045 if(!fpInstance)
00046 fpInstance = new G4MoleculeCounter();
00047
00048 return fpInstance;
00049 }
00050
00051 void G4MoleculeCounter::DeleteInstance()
00052 {
00053 if(fpInstance)
00054 {
00055 delete fpInstance;
00056 fpInstance = 0;
00057 }
00058 }
00059
00060 void G4MoleculeCounter::AddAMoleculeAtTime(const G4Molecule& molecule, G4double time)
00061 {
00062 if(fVerbose > 1)
00063 {
00064 G4cout<<"G4MoleculeCounter::AddAMoleculeAtTime : "<< molecule.GetName()
00065 << " at time : " << G4BestUnit(time, "Time") <<G4endl;
00066 }
00067
00068 if(fDontRegister[molecule.GetDefinition()]) return ;
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088 if(fVerbose)
00089 {
00090 G4cout<<"-------------------------"<<G4endl ;
00091 G4cout << "G4MoleculeCounter::AddAMoleculeAtTime " << G4endl;
00092 G4cout<<"!! Molecule = " << molecule.GetName() << G4endl;
00093 G4cout<<"!! At Time = "<< G4BestUnit(time, "Time") <<G4endl;
00094 }
00095
00096 CounterMapType::iterator counterMap_i = fCounterMap.find(molecule) ;
00097
00098 if(counterMap_i == fCounterMap.end())
00099 {
00100 if(fVerbose) G4cout << " !! ***** Map is empty " << G4endl;
00101 fCounterMap[molecule][time] = 1;
00102 }
00103 else if(counterMap_i->second.empty())
00104 {
00105 if(fVerbose) G4cout << " !! ***** Map is empty " << G4endl;
00106 counterMap_i->second[time] = 1;
00107 }
00108 else
00109 {
00110 NbMoleculeAgainstTime::iterator end = counterMap_i->second.end();
00111 end--;
00112
00113 if(fVerbose)
00114 G4cout<<"!! End Time = "<< G4BestUnit(end->first, "Time") <<G4endl;
00115
00116 if(end->first <= time)
00117 {
00118 counterMap_i->second[time]=end->second + 1;
00119 }
00120 else
00121 {
00122 NbMoleculeAgainstTime::iterator it = counterMap_i->second.lower_bound(time);
00123
00124 while(it->first > time && it!=counterMap_i->second.begin())
00125 {
00126 if(fVerbose)
00127 G4cout<<"!! ********** Is going back!!!!"<<G4endl;
00128 it--;
00129 }
00130
00131 if(it==counterMap_i->second.begin() && it->first > time)
00132 {
00133 if(fVerbose)
00134 G4cout<<"!! ********** Illegal !!!!"<<G4endl;
00135 return ;
00136 }
00137
00138 if(fVerbose)
00139 {
00140 G4cout<<"!! PREVIOUS NB = "<< it->second <<G4endl;
00141 G4cout<<"!! PREVIOUS TIME = "<< G4BestUnit(it->first,"Time") <<G4endl;
00142 }
00143 counterMap_i->second[time]=it->second + 1;
00144 }
00145 }
00146
00147 if(fVerbose)
00148 G4cout<<"!! NB = "<< fCounterMap[molecule][time]<<G4endl;
00149 }
00150
00151 void G4MoleculeCounter::RemoveAMoleculeAtTime(const G4Molecule& molecule, G4double time)
00152 {
00153 if(fVerbose > 1)
00154 {
00155 G4cout<<"-------------------------"<<G4endl ;
00156 G4cout<<"G4MoleculeCounter::RemoveAMoleculeAtTime : "<< molecule.GetName()
00157 << " at time : " << G4BestUnit(time,"Time") <<G4endl;
00158 G4cout<<"-------------------------"<<G4endl ;
00159 }
00160
00161 if(fDontRegister[molecule.GetDefinition()]) return ;
00162
00163 NbMoleculeAgainstTime& nbMolPerTime = fCounterMap[molecule];
00164
00165 if(nbMolPerTime.empty())
00166 {
00167 molecule.PrintState();
00168 G4String errMsg = "You are trying to remove molecule "
00169 + molecule.GetName()
00170 +" from the counter while this kind of molecules has not been registered yet";
00171 G4Exception("G4MoleculeCounter::RemoveAMoleculeAtTime","",FatalErrorInArgument, errMsg);
00172
00173 return;
00174 }
00175 else
00176 {
00177 NbMoleculeAgainstTime::iterator it ;
00178
00179 if(nbMolPerTime.size() == 1)
00180 {
00181 it = nbMolPerTime.begin() ;
00182 if(fVerbose)
00183 G4cout << "!! fCounterMap[molecule].size() == 1" << G4endl;
00184 }
00185 else
00186 {
00187 it = nbMolPerTime.lower_bound(time);
00188 }
00189
00190 if(it==nbMolPerTime.end())
00191 {
00192 if(fVerbose)
00193 G4cout << " ********** NO ITERATOR !!!!!!!!! " << G4endl;
00194 it--;
00195
00196 if(time<it->first)
00197 {
00198 G4String errMsg = "There was no "+ molecule.GetName()
00199 + " record at the time or even before the time asked";
00200 G4Exception("G4MoleculeCounter::RemoveAMoleculeAtTime","",FatalErrorInArgument, errMsg);
00201 }
00202 }
00203
00204 if(fVerbose)
00205 {
00206
00207 G4cout<<"!! Molecule = " << molecule.GetName() << G4endl;
00208 G4cout<<"!! At Time = "<< G4BestUnit(time,"Time") <<G4endl;
00209 G4cout<<"!! PREVIOUS TIME = "<< G4BestUnit(it->first,"Time")<<G4endl;
00210 G4cout<<"!! PREVIOUS Nb = "<< it->second <<G4endl;
00211 }
00212
00213
00214
00215 if(nbMolPerTime.value_comp()(*it, *nbMolPerTime.begin()))
00216 {
00217 if(fVerbose)
00218 G4cout<<"!! ***** In value_comp ... " << G4endl;
00219 it++;
00220 if(time<it->first)
00221 {
00222 G4String errMsg = "There was no "+ molecule.GetName()
00223 + " record at the time or even before the time asked";
00224 G4Exception("G4MoleculeCounter::RemoveAMoleculeAtTime","",FatalErrorInArgument, errMsg);
00225 }
00226 }
00227
00228 while(it->first - time > compDoubleWithPrecision::fPrecision && it!=nbMolPerTime.begin())
00229 {
00230 if(fVerbose)
00231 {
00232 G4cout<<"!! ***** Is going back!!!!"<<G4endl;
00233 G4cout<<"!! PREVIOUS TIME = "<< G4BestUnit(it-> first,"Time") <<G4endl;
00234 }
00235 it--;
00236 }
00237
00238 if(it==nbMolPerTime.begin() && it->first > time)
00239 {
00240 if(fVerbose)
00241 G4cout<<"!! ********** Illegal !!!!"<<G4endl;
00242 return ;
00243 }
00244
00245 if(fVerbose)
00246 {
00247 G4cout<<"!! PREVIOUS NB = "<< (*it).second <<G4endl;
00248 G4cout<<"!! PREVIOUS TIME = "<< G4BestUnit(it->first,"Time")<<G4endl;
00249 }
00250 nbMolPerTime[time]=it->second - 1;
00251 }
00252 if(fVerbose)
00253 {
00254 G4cout<<"!! NB = "<< nbMolPerTime[time]<<G4endl;
00255 }
00256 }
00257
00258 std::auto_ptr<vector<G4Molecule> > G4MoleculeCounter::GetRecordedMolecules()
00259 {
00260 if(fVerbose > 1)
00261 {
00262 G4cout<<"Entering in G4MoleculeCounter::RecordMolecules"<<G4endl;
00263 }
00264
00265 CounterMapType::iterator it;
00266 std::auto_ptr< vector<G4Molecule> > output (new vector<G4Molecule>) ;
00267
00268 for(it = fCounterMap.begin() ; it != fCounterMap.end() ; it++)
00269 {
00270 output->push_back((*it).first);
00271 }
00272 return output;
00273 }
00274