G4MoleculeCounter.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id: G4MoleculeCounter.cc 65022 2012-11-12 16:43:12Z gcosmo $
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     // G4double pstime = time/picosecond;
00071     //
00072     // G4double fractpart=-1, intpart=-1;
00073     // fractpart = modf (pstime , &intpart);
00074     //
00075     // if(pstime<10.1)
00076     // {
00077     // fractpart *= 10 ;
00078     // fractpart = floor(fractpart)/10;
00079     // pstime = intpart+fractpart;
00080     // }
00081     //
00082     // else
00083     // {
00084     // pstime = intpart;
00085     // }
00086     // time = pstime*picosecond ;
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 //            G4cout << "G4MoleculeCounter::RemoveAMoleculeAtTime " << G4endl;
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         // If valgrind problem on the line below, it means that the pointer "it"
00214         // points nowhere
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 

Generated on Mon May 27 17:48:53 2013 for Geant4 by  doxygen 1.4.7