G4AdjointPrimaryGeneratorAction.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$
00027 //
00029 //      Class Name:     G4AdjointPrimaryGeneratorAction
00030 //      Author:         L. Desorgher
00031 //      Organisation:   SpaceIT GmbH
00032 //      Contract:       ESA contract 21435/08/NL/AT
00033 //      Customer:       ESA/ESTEC
00035 
00036 #include "G4AdjointPrimaryGeneratorAction.hh"
00037 
00038 #include "G4PhysicalConstants.hh"
00039 #include "G4Event.hh"
00040 #include "G4ParticleTable.hh"
00041 #include "G4ParticleDefinition.hh"
00042 #include "G4ParticleTable.hh" 
00043 #include "G4AdjointSimManager.hh" 
00044 #include "G4AdjointPrimaryGenerator.hh"
00046 //
00047 G4AdjointPrimaryGeneratorAction::G4AdjointPrimaryGeneratorAction()
00048   : Emin(0.), Emax(0.), EminIon(0.), EmaxIon(0.), NbOfAdjointPrimaryTypes(0),
00049     index_particle(100000), last_generated_part_was_adjoint(false),
00050     radius_spherical_source(0.), fwd_ion(0), adj_ion(0), 
00051     ion_name("not_defined")
00052 {
00053   theAdjointPrimaryGenerator= new G4AdjointPrimaryGenerator();
00054 
00055   PrimariesConsideredInAdjointSim[G4String("e-")]=false;
00056   PrimariesConsideredInAdjointSim[G4String("gamma")]=false;
00057   PrimariesConsideredInAdjointSim[G4String("proton")]=false;
00058   PrimariesConsideredInAdjointSim[G4String("ion")]=false;
00059 
00060   ListOfPrimaryFwdParticles.clear();
00061   ListOfPrimaryAdjParticles.clear();
00062 }
00064 //
00065 G4AdjointPrimaryGeneratorAction::~G4AdjointPrimaryGeneratorAction()
00066 {
00067   delete theAdjointPrimaryGenerator;
00068 }
00070 //
00071 void G4AdjointPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
00072 {
00073    if ( !last_generated_part_was_adjoint ) {
00074          
00075          index_particle++;
00076          if (index_particle >= ListOfPrimaryAdjParticles.size()) index_particle =0;
00077          
00078         
00079          G4double E1=Emin;
00080          G4double E2=Emax;
00081          if (!ListOfPrimaryAdjParticles[index_particle]) UpdateListOfPrimaryParticles();//ion has not been created yet
00082         
00083          if (ListOfPrimaryAdjParticles[index_particle]->GetParticleName() == "adj_proton") {
00084                 E1=EminIon;
00085                 E2=EmaxIon;
00086          }
00087          if (ListOfPrimaryAdjParticles[index_particle]->GetParticleType() == "adjoint_nucleus") {
00088                 G4int A= ListOfPrimaryAdjParticles[index_particle]->GetAtomicMass();
00089                 E1=EminIon*A;
00090                 E2=EmaxIon*A;
00091          }
00092          theAdjointPrimaryGenerator->GenerateAdjointPrimaryVertex(anEvent,
00093                                                                   ListOfPrimaryAdjParticles[index_particle],
00094                                                                   E1,E2);
00095          G4PrimaryVertex* aPrimVertex = anEvent->GetPrimaryVertex();
00096   
00097   
00098           p=aPrimVertex->GetPrimary()->GetMomentum();
00099           pos=aPrimVertex->GetPosition();
00100           G4double pmag=p.mag();
00101           
00102           G4double m0=ListOfPrimaryAdjParticles[index_particle]->GetPDGMass();
00103           G4double ekin=std::sqrt( m0*m0 + pmag*pmag) -m0;
00104         
00105   
00106           //The factor pi is to normalise the weight to the directional flux
00107           G4double adjoint_source_area = G4AdjointSimManager::GetInstance()->GetAdjointSourceArea();
00108           G4double adjoint_weight = ComputeEnergyDistWeight(ekin,E1,E2)*adjoint_source_area*pi;
00109 
00110           aPrimVertex->SetWeight(adjoint_weight);
00111 
00112           last_generated_part_was_adjoint =true;
00113           G4AdjointSimManager::GetInstance()->SetAdjointTrackingMode(true);
00114           G4AdjointSimManager::GetInstance()->RegisterAdjointPrimaryWeight(adjoint_weight);
00115    }
00116    else {
00117           //fwd particle equivalent to the last generated adjoint particle ios generated        
00118           G4PrimaryVertex* aPrimVertex = new G4PrimaryVertex();
00119           aPrimVertex->SetPosition(pos.x(),pos.y(),pos.z());
00120           aPrimVertex->SetT0(0.);
00121           G4PrimaryParticle* aPrimParticle = new G4PrimaryParticle(ListOfPrimaryFwdParticles[index_particle],
00122                                                            -p.x(),-p.y(),-p.z()); 
00123         
00124           aPrimVertex->SetPrimary(aPrimParticle);
00125           anEvent->AddPrimaryVertex(aPrimVertex);                                                  
00126           last_generated_part_was_adjoint =false;
00127           G4AdjointSimManager::GetInstance()->SetAdjointTrackingMode(false);
00128    }               
00129 }
00131 //
00132 void G4AdjointPrimaryGeneratorAction::SetEmin(G4double val)
00133 {
00134   Emin=val;
00135   EminIon=val;
00136 }
00138 //
00139 void G4AdjointPrimaryGeneratorAction::SetEmax(G4double val)
00140 {
00141   Emax=val;
00142   EmaxIon=val;
00143 }
00145 //
00146 void G4AdjointPrimaryGeneratorAction::SetEminIon(G4double val)
00147 {
00148   EminIon=val;
00149 }
00151 //
00152 void G4AdjointPrimaryGeneratorAction::SetEmaxIon(G4double val)
00153 {
00154   EmaxIon=val;
00155 }
00157 //
00158 G4double G4AdjointPrimaryGeneratorAction::ComputeEnergyDistWeight(G4double E ,G4double E1, G4double E2)
00159 {
00160  //  We generate N numbers of primaries with  a 1/E energy law distribution.
00161  //  We have therefore  an energy distribution function    
00162  //             f(E)=C/E  (1)
00163  //  with   C a constant that is such that
00164  //             N=Integral(f(E),E1,E2)=C.std::log(E2/E1)  (2)
00165  //  Therefore from (2) we get          
00166  //             C=N/ std::log(E2/E1) (3) 
00167  //  and        
00168  //             f(E)=N/ std::log(E2/E1)/E (4)
00169  //For the adjoint simulation we need a energy distribution  f'(E)=1..
00170  //To get that we need therefore to apply a weight to the primary
00171  //             W=1/f(E)=E*std::log(E2/E1)/N  
00172  //
00173   return std::log(E2/E1)*E/G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun();
00174  
00175 }
00177 //
00178 void G4AdjointPrimaryGeneratorAction::SetSphericalAdjointPrimarySource(G4double radius, G4ThreeVector center_pos)
00179 { 
00180   radius_spherical_source = radius;
00181   center_spherical_source = center_pos;
00182   type_of_adjoint_source ="Spherical";
00183   theAdjointPrimaryGenerator->SetSphericalAdjointPrimarySource(radius,center_pos);
00184 }
00186 //
00187 void G4AdjointPrimaryGeneratorAction::SetAdjointPrimarySourceOnAnExtSurfaceOfAVolume(const G4String& volume_name)
00188 {
00189   type_of_adjoint_source ="ExternalSurfaceOfAVolume";
00190   theAdjointPrimaryGenerator->SetAdjointPrimarySourceOnAnExtSurfaceOfAVolume(volume_name);
00191 }
00193 //
00194 void G4AdjointPrimaryGeneratorAction::ConsiderParticleAsPrimary(const G4String& particle_name)
00195 {
00196   if (PrimariesConsideredInAdjointSim.find(particle_name) != PrimariesConsideredInAdjointSim.end()){
00197         PrimariesConsideredInAdjointSim[particle_name]=true;
00198   }
00199   UpdateListOfPrimaryParticles();
00200 }
00202 //
00203 void G4AdjointPrimaryGeneratorAction::NeglectParticleAsPrimary(const G4String& particle_name)
00204 {
00205   if (PrimariesConsideredInAdjointSim.find(particle_name) != PrimariesConsideredInAdjointSim.end()){
00206         PrimariesConsideredInAdjointSim[particle_name]= false;
00207   }
00208   UpdateListOfPrimaryParticles();
00209 }
00211 //
00212 void G4AdjointPrimaryGeneratorAction::UpdateListOfPrimaryParticles()
00213 {
00214     G4ParticleTable* theParticleTable = G4ParticleTable::GetParticleTable();
00215     ListOfPrimaryFwdParticles.clear();
00216     ListOfPrimaryAdjParticles.clear();
00217     std::map<G4String, G4bool>::iterator iter;
00218     for( iter = PrimariesConsideredInAdjointSim.begin(); iter != PrimariesConsideredInAdjointSim.end(); ++iter ) {
00219         if(iter->second) {
00220                 G4String fwd_particle_name = iter->first;
00221                 if ( fwd_particle_name != "ion") {
00222                         G4String adj_particle_name = G4String("adj_") + fwd_particle_name;
00223                         ListOfPrimaryFwdParticles.push_back(theParticleTable->FindParticle(fwd_particle_name));
00224                         ListOfPrimaryAdjParticles.push_back(theParticleTable->FindParticle(adj_particle_name));
00225                 }
00226                 else {
00227                         if (fwd_ion ){
00228                                 ion_name=fwd_ion->GetParticleName();
00229                                 G4String adj_ion_name=G4String("adj_") +ion_name;
00230                                 ListOfPrimaryFwdParticles.push_back(fwd_ion);
00231                                 ListOfPrimaryAdjParticles.push_back(adj_ion);
00232                         }
00233                         else {
00234                                 ListOfPrimaryFwdParticles.push_back(0);
00235                                 ListOfPrimaryAdjParticles.push_back(0);
00236                                 
00237                         }       
00238                 }
00239         }       
00240    }
00241 }
00243 //
00244 void G4AdjointPrimaryGeneratorAction::SetPrimaryIon(G4ParticleDefinition* adjointIon, G4ParticleDefinition* fwdIon)
00245 {
00246   fwd_ion = fwdIon;
00247   adj_ion = adjointIon;
00248   UpdateListOfPrimaryParticles();
00249 }
00250 

Generated on Mon May 27 17:47:37 2013 for Geant4 by  doxygen 1.4.7