Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/run/src/G4AdjointPrimaryGeneratorAction.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /run/src/G4AdjointPrimaryGeneratorAction.cc (Version 11.3.0) and /run/src/G4AdjointPrimaryGeneratorAction.cc (Version 10.2.p2)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 // G4AdjointPrimaryGeneratorAction implementat <<  26 // $Id: G4AdjointPrimaryGeneratorAction.cc 81773 2014-06-05 08:35:38Z gcosmo $
 27 //                                                 27 //
 28 // ------------------------------------------- <<  28 /////////////////////////////////////////////////////////////////////////////
 29 //   Class Name:   G4AdjointPrimaryGeneratorAc <<  29 //      Class Name: G4AdjointPrimaryGeneratorAction
 30 //   Author:       L. Desorgher, 2007-2009     <<  30 //  Author:         L. Desorgher
 31 //   Organisation: SpaceIT GmbH                <<  31 //  Organisation:   SpaceIT GmbH
 32 //   Contract:     ESA contract 21435/08/NL/AT <<  32 //  Contract: ESA contract 21435/08/NL/AT
 33 //   Customer:     ESA/ESTEC                   <<  33 //  Customer:       ESA/ESTEC
 34 // ------------------------------------------- <<  34 /////////////////////////////////////////////////////////////////////////////
 35                                                    35 
 36 #include "G4AdjointPrimaryGeneratorAction.hh"      36 #include "G4AdjointPrimaryGeneratorAction.hh"
 37                                                    37 
 38 #include "G4AdjointPrimaryGenerator.hh"        <<  38 #include "G4PhysicalConstants.hh"
 39 #include "G4AdjointSimManager.hh"              << 
 40 #include "G4Event.hh"                              39 #include "G4Event.hh"
 41 #include "G4Gamma.hh"                          << 
 42 #include "G4ParticleDefinition.hh"             << 
 43 #include "G4ParticleTable.hh"                      40 #include "G4ParticleTable.hh"
 44 #include "G4PhysicalConstants.hh"              <<  41 #include "G4ParticleDefinition.hh"
                                                   >>  42 #include "G4ParticleTable.hh" 
                                                   >>  43 #include "G4AdjointSimManager.hh" 
                                                   >>  44 #include "G4AdjointPrimaryGenerator.hh"
                                                   >>  45 #include "G4Gamma.hh"
 45                                                    46 
 46 // ------------------------------------------- <<  47 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 47 //                                                 48 //
 48 G4AdjointPrimaryGeneratorAction::G4AdjointPrim     49 G4AdjointPrimaryGeneratorAction::G4AdjointPrimaryGeneratorAction()
 49 {                                              <<  50   : Emin(0.), Emax(0.), EminIon(0.), EmaxIon(0.),
 50   theAdjointPrimaryGenerator = new G4AdjointPr <<  51     index_particle(100000),
 51                                                <<  52     radius_spherical_source(0.), fwd_ion(0), adj_ion(0), 
 52   PrimariesConsideredInAdjointSim[G4String("e- <<  53     ion_name("not_defined")
 53   PrimariesConsideredInAdjointSim[G4String("ga <<  54 {
 54   PrimariesConsideredInAdjointSim[G4String("pr <<  55   theAdjointPrimaryGenerator= new G4AdjointPrimaryGenerator();
 55   PrimariesConsideredInAdjointSim[G4String("io <<  56 
                                                   >>  57   PrimariesConsideredInAdjointSim[G4String("e-")]=false;
                                                   >>  58   PrimariesConsideredInAdjointSim[G4String("gamma")]=false;
                                                   >>  59   PrimariesConsideredInAdjointSim[G4String("proton")]=false;
                                                   >>  60   PrimariesConsideredInAdjointSim[G4String("ion")]=false;
 56                                                    61 
 57   ListOfPrimaryFwdParticles.clear();               62   ListOfPrimaryFwdParticles.clear();
 58   ListOfPrimaryAdjParticles.clear();               63   ListOfPrimaryAdjParticles.clear();
                                                   >>  64   nb_fwd_gammas_per_event = 10;
 59 }                                                  65 }
 60                                                <<  66 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 61 // ------------------------------------------- << 
 62 //                                                 67 //
 63 G4AdjointPrimaryGeneratorAction::~G4AdjointPri     68 G4AdjointPrimaryGeneratorAction::~G4AdjointPrimaryGeneratorAction()
 64 {                                                  69 {
 65   delete theAdjointPrimaryGenerator;               70   delete theAdjointPrimaryGenerator;
 66 }                                                  71 }
 67                                                <<  72 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 68 // ------------------------------------------- << 
 69 //                                                 73 //
 70 void G4AdjointPrimaryGeneratorAction::Generate     74 void G4AdjointPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
 71 {                                              <<  75 {  G4int evt_id=anEvent->GetEventID();
 72   G4int evt_id = anEvent->GetEventID();        <<  76    size_t n=ListOfPrimaryAdjParticles.size();
 73   std::size_t n = ListOfPrimaryAdjParticles.si <<  77    index_particle=size_t(evt_id)-n*(size_t(evt_id)/n);
 74   index_particle = std::size_t(evt_id) - n * ( <<  78 
 75                                                <<  79 
 76   G4double E1 = Emin;                          <<  80    G4double E1=Emin;
 77   G4double E2 = Emax;                          <<  81    G4double E2=Emax;
 78   if (ListOfPrimaryAdjParticles[index_particle <<  82    if (!ListOfPrimaryAdjParticles[index_particle]) UpdateListOfPrimaryParticles();//ion has not been created yet
 79     UpdateListOfPrimaryParticles();  // ion ha <<  83 
 80                                                <<  84    if (ListOfPrimaryAdjParticles[index_particle]->GetParticleName() == "adj_proton") {
 81   if (ListOfPrimaryAdjParticles[index_particle <<  85      E1=EminIon;
 82     E1 = EminIon;                              <<  86      E2=EmaxIon;
 83     E2 = EmaxIon;                              <<  87    }
 84   }                                            <<  88    if (ListOfPrimaryAdjParticles[index_particle]->GetParticleType() == "adjoint_nucleus") {
 85   if (ListOfPrimaryAdjParticles[index_particle <<  89      G4int A= ListOfPrimaryAdjParticles[index_particle]->GetAtomicMass();
 86     G4int A = ListOfPrimaryAdjParticles[index_ <<  90      E1=EminIon*A;
 87     E1 = EminIon * A;                          <<  91      E2=EmaxIon*A;
 88     E2 = EmaxIon * A;                          <<  92    }
 89   }                                            <<  93    theAdjointPrimaryGenerator->GenerateFwdPrimaryVertex(anEvent,ListOfPrimaryFwdParticles[index_particle],E1,E2);
 90   // Generate first the forwrad primaries      <<  94    G4PrimaryVertex* fwdPrimVertex = anEvent->GetPrimaryVertex();
 91   theAdjointPrimaryGenerator->GenerateFwdPrima <<  95 
 92     anEvent, ListOfPrimaryFwdParticles[index_p <<  96    p=fwdPrimVertex->GetPrimary()->GetMomentum();
 93   G4PrimaryVertex* fwdPrimVertex = anEvent->Ge <<  97    pos=fwdPrimVertex->GetPosition();
 94                                                <<  98    G4double pmag=p.mag();
 95   p = fwdPrimVertex->GetPrimary()->GetMomentum <<  99    G4double m0=ListOfPrimaryFwdParticles[index_particle]->GetPDGMass();
 96   pos = fwdPrimVertex->GetPosition();          << 100    G4double ekin=std::sqrt( m0*m0 + pmag*pmag) -m0;
 97   G4double pmag = p.mag();                     << 101 
 98   G4double m0 = ListOfPrimaryFwdParticles[inde << 102    G4double weight_correction=1.;
 99   G4double ekin = std::sqrt(m0 * m0 + pmag * p << 103    //For gamma generate the particle along the backward ray
100                                                << 104    G4ThreeVector dir=-p/p.mag();
101   G4double weight_correction = 1.;             << 105 
102   // For gamma generate the particle along the << 106    /*if (ListOfPrimaryAdjParticles[index_particle]->GetParticleName() == "adj_gamma"){
103   G4ThreeVector dir = -p / p.mag();            << 107 
104                                                << 108       theAdjointPrimaryGenerator
105   weight_correction = 1.;                      << 109               ->ComputeAccumulatedDepthVectorAlongBackRay(pos,dir,ekin,ListOfPrimaryAdjParticles[index_particle]);
106                                                << 110 
107   if (ListOfPrimaryFwdParticles[index_particle << 111 
108   {                                            << 112       G4double distance = theAdjointPrimaryGenerator
109     G4double weight = (1. / nb_fwd_gammas_per_ << 113        ->SampleDistanceAlongBackRayAndComputeWeightCorrection(weight_correction);
110     fwdPrimVertex->SetWeight(weight);          << 114 
111     for (G4int i = 0; i < nb_fwd_gammas_per_ev << 115       //pos=pos+dir*distance;
112       auto newFwdPrimVertex = new G4PrimaryVer << 116       //fwdPrimVertex->SetPosition(pos[0],pos[1],pos[2]);
113       newFwdPrimVertex->SetPosition(pos.x(), p << 117    }
114       newFwdPrimVertex->SetT0(0.);             << 118    */
115       auto aPrimParticle =                     << 119    weight_correction=1.;
116         new G4PrimaryParticle(ListOfPrimaryFwd << 120 
117       newFwdPrimVertex->SetPrimary(aPrimPartic << 121    if (ListOfPrimaryFwdParticles[index_particle]  ==G4Gamma::Gamma() && nb_fwd_gammas_per_event>1 ){
118       newFwdPrimVertex->SetWeight(weight);     << 122      G4double weight = (1./nb_fwd_gammas_per_event);
119       anEvent->AddPrimaryVertex(newFwdPrimVert << 123      fwdPrimVertex->SetWeight(weight);
120     }                                          << 124      for (int i=0;i<nb_fwd_gammas_per_event-1;i++){
121   }                                            << 125          G4PrimaryVertex* newFwdPrimVertex = new G4PrimaryVertex();
                                                   >> 126          newFwdPrimVertex->SetPosition(pos.x(),pos.y(),pos.z());
                                                   >> 127          newFwdPrimVertex->SetT0(0.);
                                                   >> 128          G4PrimaryParticle* aPrimParticle = new G4PrimaryParticle(ListOfPrimaryFwdParticles[index_particle],
                                                   >> 129                p.x(),p.y(),p.z());
                                                   >> 130          newFwdPrimVertex->SetPrimary(aPrimParticle);
                                                   >> 131          newFwdPrimVertex->SetWeight(weight);
                                                   >> 132          anEvent->AddPrimaryVertex(newFwdPrimVertex);
                                                   >> 133        }
                                                   >> 134    }
                                                   >> 135 
                                                   >> 136 
                                                   >> 137    G4PrimaryVertex* adjPrimVertex = new G4PrimaryVertex();
                                                   >> 138    adjPrimVertex->SetPosition(pos.x(),pos.y(),pos.z());
                                                   >> 139    adjPrimVertex->SetT0(0.);
                                                   >> 140    G4PrimaryParticle* aPrimParticle = new G4PrimaryParticle(ListOfPrimaryAdjParticles[index_particle],
                                                   >> 141                          -p.x(),-p.y(),-p.z());
                                                   >> 142 
                                                   >> 143    adjPrimVertex->SetPrimary(aPrimParticle);
                                                   >> 144    anEvent->AddPrimaryVertex(adjPrimVertex);
                                                   >> 145 
                                                   >> 146    //The factor pi is to normalise the weight to the directional flux
                                                   >> 147    G4double adjoint_source_area = G4AdjointSimManager::GetInstance()->GetAdjointSourceArea();
                                                   >> 148    G4double adjoint_weight = weight_correction*ComputeEnergyDistWeight(ekin,E1,E2)*adjoint_source_area*pi;
                                                   >> 149    if (ListOfPrimaryFwdParticles[index_particle]  ==G4Gamma::Gamma()) adjoint_weight = adjoint_weight/3.;
                                                   >> 150    adjPrimVertex->SetWeight(adjoint_weight);
                                                   >> 151 
                                                   >> 152 
                                                   >> 153 
                                                   >> 154 
                                                   >> 155    G4AdjointSimManager::GetInstance()->SetAdjointTrackingMode(true);
                                                   >> 156 
                                                   >> 157 
                                                   >> 158 
                                                   >> 159   /* if ( !last_generated_part_was_adjoint ) {
                                                   >> 160    
                                                   >> 161    index_particle++;
                                                   >> 162    if (index_particle >= ListOfPrimaryAdjParticles.size()) index_particle =0;
                                                   >> 163    
                                                   >> 164   
                                                   >> 165    G4double E1=Emin;
                                                   >> 166    G4double E2=Emax;
                                                   >> 167    if (!ListOfPrimaryAdjParticles[index_particle]) UpdateListOfPrimaryParticles();//ion has not been created yet
                                                   >> 168   
                                                   >> 169    if (ListOfPrimaryAdjParticles[index_particle]->GetParticleName() == "adj_proton") {
                                                   >> 170     E1=EminIon;
                                                   >> 171     E2=EmaxIon;
                                                   >> 172    }
                                                   >> 173    if (ListOfPrimaryAdjParticles[index_particle]->GetParticleType() == "adjoint_nucleus") {
                                                   >> 174     G4int A= ListOfPrimaryAdjParticles[index_particle]->GetAtomicMass();
                                                   >> 175     E1=EminIon*A;
                                                   >> 176     E2=EmaxIon*A;
                                                   >> 177    }
                                                   >> 178    theAdjointPrimaryGenerator->GenerateAdjointPrimaryVertex(anEvent,
                                                   >> 179                   ListOfPrimaryAdjParticles[index_particle],
                                                   >> 180                   E1,E2);
                                                   >> 181    G4PrimaryVertex* aPrimVertex = anEvent->GetPrimaryVertex();
                                                   >> 182   
                                                   >> 183   
                                                   >> 184       p=aPrimVertex->GetPrimary()->GetMomentum();
                                                   >> 185     pos=aPrimVertex->GetPosition();
                                                   >> 186     G4double pmag=p.mag();
                                                   >> 187     
                                                   >> 188     G4double m0=ListOfPrimaryAdjParticles[index_particle]->GetPDGMass();
                                                   >> 189     G4double ekin=std::sqrt( m0*m0 + pmag*pmag) -m0;
                                                   >> 190     
                                                   >> 191   
                                                   >> 192       //The factor pi is to normalise the weight to the directional flux
                                                   >> 193     G4double adjoint_source_area = G4AdjointSimManager::GetInstance()->GetAdjointSourceArea();
                                                   >> 194       G4double adjoint_weight = ComputeEnergyDistWeight(ekin,E1,E2)*adjoint_source_area*pi;
                                                   >> 195 
                                                   >> 196       aPrimVertex->SetWeight(adjoint_weight);
                                                   >> 197 
                                                   >> 198     last_generated_part_was_adjoint =true;
                                                   >> 199     G4AdjointSimManager::GetInstance()->SetAdjointTrackingMode(true);
                                                   >> 200     G4AdjointSimManager::GetInstance()->RegisterAdjointPrimaryWeight(adjoint_weight);
                                                   >> 201    }
                                                   >> 202    else {
                                                   >> 203           //fwd particle equivalent to the last generated adjoint particle ios generated  
                                                   >> 204       G4PrimaryVertex* aPrimVertex = new G4PrimaryVertex();
                                                   >> 205     aPrimVertex->SetPosition(pos.x(),pos.y(),pos.z());
                                                   >> 206       aPrimVertex->SetT0(0.);
                                                   >> 207     G4PrimaryParticle* aPrimParticle = new G4PrimaryParticle(ListOfPrimaryFwdParticles[index_particle],
                                                   >> 208                    -p.x(),-p.y(),-p.z()); 
                                                   >> 209   
                                                   >> 210     aPrimVertex->SetPrimary(aPrimParticle);
                                                   >> 211     anEvent->AddPrimaryVertex(aPrimVertex);                
                                                   >> 212     last_generated_part_was_adjoint =false;
                                                   >> 213     G4AdjointSimManager::GetInstance()->SetAdjointTrackingMode(false);
                                                   >> 214   */
122                                                   215 
123   // Now generate the adjoint primaries        << 
124   auto adjPrimVertex = new G4PrimaryVertex();  << 
125   adjPrimVertex->SetPosition(pos.x(), pos.y(), << 
126   adjPrimVertex->SetT0(0.);                    << 
127   auto aPrimParticle =                         << 
128     new G4PrimaryParticle(ListOfPrimaryAdjPart << 
129                                                << 
130   adjPrimVertex->SetPrimary(aPrimParticle);    << 
131   anEvent->AddPrimaryVertex(adjPrimVertex);    << 
132                                                << 
133   // The factor pi is to normalise the weight  << 
134   G4double adjoint_source_area = G4AdjointSimM << 
135   G4double adjoint_weight =                    << 
136     weight_correction * ComputeEnergyDistWeigh << 
137   if (ListOfPrimaryAdjParticles[index_particle << 
138     // The weight will be corrected at the end << 
139     // are used                                << 
140     adjoint_weight = adjoint_weight / nb_adj_p << 
141     for (G4int i = 0; i < nb_adj_primary_gamma << 
142       auto newAdjPrimVertex = new G4PrimaryVer << 
143       newAdjPrimVertex->SetPosition(pos.x(), p << 
144       newAdjPrimVertex->SetT0(0.);             << 
145       aPrimParticle =                          << 
146         new G4PrimaryParticle(ListOfPrimaryAdj << 
147       newAdjPrimVertex->SetPrimary(aPrimPartic << 
148       newAdjPrimVertex->SetWeight(adjoint_weig << 
149       anEvent->AddPrimaryVertex(newAdjPrimVert << 
150     }                                          << 
151   }                                            << 
152   else if (ListOfPrimaryAdjParticles[index_par << 
153     // The weight will be corrected at the end << 
154     // are used                                << 
155     adjoint_weight = adjoint_weight / nb_adj_p << 
156     for (G4int i = 0; i < nb_adj_primary_elect << 
157       auto newAdjPrimVertex = new G4PrimaryVer << 
158       newAdjPrimVertex->SetPosition(pos.x(), p << 
159       newAdjPrimVertex->SetT0(0.);             << 
160       aPrimParticle =                          << 
161         new G4PrimaryParticle(ListOfPrimaryAdj << 
162       newAdjPrimVertex->SetPrimary(aPrimPartic << 
163       newAdjPrimVertex->SetWeight(adjoint_weig << 
164       anEvent->AddPrimaryVertex(newAdjPrimVert << 
165     }                                          << 
166   }                                            << 
167   adjPrimVertex->SetWeight(adjoint_weight);    << 
168                                                << 
169   // Call some methods of G4AdjointSimManager  << 
170   G4AdjointSimManager::GetInstance()->SetAdjoi << 
171   G4AdjointSimManager::GetInstance()->ClearEnd << 
172   G4AdjointSimManager::GetInstance()->ResetDid << 
173 }                                                 216 }
174                                                << 217 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
175 // ------------------------------------------- << 
176 //                                                218 //
177 void G4AdjointPrimaryGeneratorAction::SetEmin(    219 void G4AdjointPrimaryGeneratorAction::SetEmin(G4double val)
178 {                                                 220 {
179   Emin = val;                                  << 221   Emin=val;
180   EminIon = val;                               << 222   EminIon=val;
181 }                                                 223 }
182                                                << 224 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
183 // ------------------------------------------- << 
184 //                                                225 //
185 void G4AdjointPrimaryGeneratorAction::SetEmax(    226 void G4AdjointPrimaryGeneratorAction::SetEmax(G4double val)
186 {                                                 227 {
187   Emax = val;                                  << 228   Emax=val;
188   EmaxIon = val;                               << 229   EmaxIon=val;
189 }                                                 230 }
190                                                << 231 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
191 // ------------------------------------------- << 
192 //                                                232 //
193 void G4AdjointPrimaryGeneratorAction::SetEminI    233 void G4AdjointPrimaryGeneratorAction::SetEminIon(G4double val)
194 {                                                 234 {
195   EminIon = val;                               << 235   EminIon=val;
196 }                                                 236 }
197                                                << 237 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
198 // ------------------------------------------- << 
199 //                                                238 //
200 void G4AdjointPrimaryGeneratorAction::SetEmaxI    239 void G4AdjointPrimaryGeneratorAction::SetEmaxIon(G4double val)
201 {                                                 240 {
202   EmaxIon = val;                               << 241   EmaxIon=val;
203 }                                                 242 }
204                                                << 243 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
205 // ------------------------------------------- << 
206 //                                                244 //
207 G4double G4AdjointPrimaryGeneratorAction::Comp << 245 G4double G4AdjointPrimaryGeneratorAction::ComputeEnergyDistWeight(G4double E ,G4double E1, G4double E2)
208                                                << 
209 {                                                 246 {
210   // We generate N numbers of primaries with   << 247  //  We generate N numbers of primaries with  a 1/E energy law distribution.
211   // We have therefore  an energy distribution << 248  //  We have therefore  an energy distribution function    
212   //   f(E)=C/E  (1)                           << 249  //     f(E)=C/E  (1)
213   // with   C a constant that is such that     << 250  //  with   C a constant that is such that
214   //   N=Integral(f(E),E1,E2)=C.std::log(E2/E1 << 251  //   N=Integral(f(E),E1,E2)=C.std::log(E2/E1)  (2)
215   // Therefore from (2) we get                 << 252  //  Therefore from (2) we get    
216   //   C=N/ std::log(E2/E1) (3)                << 253  //     C=N/ std::log(E2/E1) (3) 
217   // and                                       << 254  //  and        
218   //   f(E)=N/ std::log(E2/E1)/E (4)           << 255  //   f(E)=N/ std::log(E2/E1)/E (4)
219   // For the adjoint simulation we need a ener << 256  //For the adjoint simulation we need a energy distribution  f'(E)=1..
220   // To get that we need therefore to apply a  << 257  //To get that we need therefore to apply a weight to the primary
221   //   W=1/f(E)=E*std::log(E2/E1)/N            << 258  //   W=1/f(E)=E*std::log(E2/E1)/N  
222   //                                           << 259  //
223   return std::log(E2 / E1) * E / G4AdjointSimM << 260   return std::log(E2/E1)*E/G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun();
                                                   >> 261  
224 }                                                 262 }
225                                                << 263 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
226 // ------------------------------------------- << 
227 //                                                264 //
228 void G4AdjointPrimaryGeneratorAction::SetSpher << 265 void G4AdjointPrimaryGeneratorAction::SetSphericalAdjointPrimarySource(G4double radius, G4ThreeVector center_pos)
229                                                << 266 { 
230 {                                              << 
231   radius_spherical_source = radius;               267   radius_spherical_source = radius;
232   center_spherical_source = center_pos;           268   center_spherical_source = center_pos;
233   type_of_adjoint_source = "Spherical";        << 269   type_of_adjoint_source ="Spherical";
234   theAdjointPrimaryGenerator->SetSphericalAdjo << 270   theAdjointPrimaryGenerator->SetSphericalAdjointPrimarySource(radius,center_pos);
235 }                                                 271 }
236                                                << 272 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
237 // ------------------------------------------- << 
238 //                                                273 //
239 void G4AdjointPrimaryGeneratorAction::SetAdjoi << 274 void G4AdjointPrimaryGeneratorAction::SetAdjointPrimarySourceOnAnExtSurfaceOfAVolume(const G4String& volume_name)
240   const G4String& volume_name)                 << 
241 {                                                 275 {
242   type_of_adjoint_source = "ExternalSurfaceOfA << 276   type_of_adjoint_source ="ExternalSurfaceOfAVolume";
243   theAdjointPrimaryGenerator->SetAdjointPrimar    277   theAdjointPrimaryGenerator->SetAdjointPrimarySourceOnAnExtSurfaceOfAVolume(volume_name);
244 }                                                 278 }
245                                                << 279 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
246 // ------------------------------------------- << 
247 //                                                280 //
248 void G4AdjointPrimaryGeneratorAction::Consider    281 void G4AdjointPrimaryGeneratorAction::ConsiderParticleAsPrimary(const G4String& particle_name)
249 {                                                 282 {
250   if (PrimariesConsideredInAdjointSim.find(par << 283   if (PrimariesConsideredInAdjointSim.find(particle_name) != PrimariesConsideredInAdjointSim.end()){
251   {                                            << 284     PrimariesConsideredInAdjointSim[particle_name]=true;
252     PrimariesConsideredInAdjointSim[particle_n << 
253   }                                               285   }
254   UpdateListOfPrimaryParticles();                 286   UpdateListOfPrimaryParticles();
255 }                                                 287 }
256                                                << 288 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
257 // ------------------------------------------- << 
258 //                                                289 //
259 void G4AdjointPrimaryGeneratorAction::NeglectP    290 void G4AdjointPrimaryGeneratorAction::NeglectParticleAsPrimary(const G4String& particle_name)
260 {                                                 291 {
261   if (PrimariesConsideredInAdjointSim.find(par << 292   if (PrimariesConsideredInAdjointSim.find(particle_name) != PrimariesConsideredInAdjointSim.end()){
262   {                                            << 293     PrimariesConsideredInAdjointSim[particle_name]= false;
263     PrimariesConsideredInAdjointSim[particle_n << 
264   }                                               294   }
265   UpdateListOfPrimaryParticles();                 295   UpdateListOfPrimaryParticles();
266 }                                                 296 }
267                                                << 297 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
268 // ------------------------------------------- << 
269 //                                                298 //
270 void G4AdjointPrimaryGeneratorAction::UpdateLi    299 void G4AdjointPrimaryGeneratorAction::UpdateListOfPrimaryParticles()
271 {                                                 300 {
272   G4ParticleTable* theParticleTable = G4Partic << 301     G4ParticleTable* theParticleTable = G4ParticleTable::GetParticleTable();
273   ListOfPrimaryFwdParticles.clear();           << 302     ListOfPrimaryFwdParticles.clear();
274   ListOfPrimaryAdjParticles.clear();           << 303     ListOfPrimaryAdjParticles.clear();
275   for (const auto& iter : PrimariesConsideredI << 304     std::map<G4String, G4bool>::iterator iter;
276     if (iter.second) {                         << 305     for( iter = PrimariesConsideredInAdjointSim.begin(); iter != PrimariesConsideredInAdjointSim.end(); ++iter ) {
277       G4String fwd_particle_name = iter.first; << 306   if(iter->second) {
278       if (fwd_particle_name != "ion") {        << 307     G4String fwd_particle_name = iter->first;
279         G4String adj_particle_name = G4String( << 308     if ( fwd_particle_name != "ion") {
280         ListOfPrimaryFwdParticles.push_back(th << 309       G4String adj_particle_name = G4String("adj_") + fwd_particle_name;
281         ListOfPrimaryAdjParticles.push_back(th << 310       ListOfPrimaryFwdParticles.push_back(theParticleTable->FindParticle(fwd_particle_name));
282       }                                        << 311       ListOfPrimaryAdjParticles.push_back(theParticleTable->FindParticle(adj_particle_name));
283       else {                                   << 312       if ( fwd_particle_name == "gamma") {
284         if (fwd_ion != nullptr) {              << 313         for (G4int i=0;i<2;i++){
285           ion_name = fwd_ion->GetParticleName( << 314           ListOfPrimaryFwdParticles.push_back(theParticleTable->FindParticle(fwd_particle_name));
286           G4String adj_ion_name = G4String("ad << 315           ListOfPrimaryAdjParticles.push_back(theParticleTable->FindParticle(adj_particle_name));
287           ListOfPrimaryFwdParticles.push_back( << 316         }
288           ListOfPrimaryAdjParticles.push_back( << 317       }
289         }                                      << 318     }
290         else {                                 << 319     else {
291           ListOfPrimaryFwdParticles.push_back( << 320       if (fwd_ion ){
292           ListOfPrimaryAdjParticles.push_back( << 321         ion_name=fwd_ion->GetParticleName();
293         }                                      << 322         G4String adj_ion_name=G4String("adj_") +ion_name;
294       }                                        << 323         ListOfPrimaryFwdParticles.push_back(fwd_ion);
295     }                                          << 324         ListOfPrimaryAdjParticles.push_back(adj_ion);
296   }                                            << 325       }
                                                   >> 326       else {
                                                   >> 327         ListOfPrimaryFwdParticles.push_back(0);
                                                   >> 328         ListOfPrimaryAdjParticles.push_back(0);
                                                   >> 329         
                                                   >> 330       } 
                                                   >> 331     }
                                                   >> 332   } 
                                                   >> 333    }
297 }                                                 334 }
298                                                << 335 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
299 // ------------------------------------------- << 
300 //                                                336 //
301 void G4AdjointPrimaryGeneratorAction::SetPrima << 337 void G4AdjointPrimaryGeneratorAction::SetPrimaryIon(G4ParticleDefinition* adjointIon, G4ParticleDefinition* fwdIon)
302                                                << 
303 {                                                 338 {
304   fwd_ion = fwdIon;                               339   fwd_ion = fwdIon;
305   adj_ion = adjointIon;                           340   adj_ion = adjointIon;
306   UpdateListOfPrimaryParticles();                 341   UpdateListOfPrimaryParticles();
307 }                                                 342 }
                                                   >> 343 
308                                                   344