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.4)


  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 102435 2017-01-27 08:28:15Z 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();
 59 }                                              <<  64   nb_fwd_gammas_per_event = 1;
                                                   >>  65   nb_adj_primary_gammas_per_event = 1;
                                                   >>  66   nb_adj_primary_electrons_per_event = 1;
 60                                                    67 
 61 // ------------------------------------------- <<  68 }
                                                   >>  69 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 62 //                                                 70 //
 63 G4AdjointPrimaryGeneratorAction::~G4AdjointPri     71 G4AdjointPrimaryGeneratorAction::~G4AdjointPrimaryGeneratorAction()
 64 {                                                  72 {
 65   delete theAdjointPrimaryGenerator;               73   delete theAdjointPrimaryGenerator;
 66 }                                                  74 }
 67                                                <<  75 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 68 // ------------------------------------------- << 
 69 //                                                 76 //
 70 void G4AdjointPrimaryGeneratorAction::Generate     77 void G4AdjointPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
 71 {                                                  78 {
 72   G4int evt_id = anEvent->GetEventID();        << 
 73   std::size_t n = ListOfPrimaryAdjParticles.si << 
 74   index_particle = std::size_t(evt_id) - n * ( << 
 75                                                << 
 76   G4double E1 = Emin;                          << 
 77   G4double E2 = Emax;                          << 
 78   if (ListOfPrimaryAdjParticles[index_particle << 
 79     UpdateListOfPrimaryParticles();  // ion ha << 
 80                                                << 
 81   if (ListOfPrimaryAdjParticles[index_particle << 
 82     E1 = EminIon;                              << 
 83     E2 = EmaxIon;                              << 
 84   }                                            << 
 85   if (ListOfPrimaryAdjParticles[index_particle << 
 86     G4int A = ListOfPrimaryAdjParticles[index_ << 
 87     E1 = EminIon * A;                          << 
 88     E2 = EmaxIon * A;                          << 
 89   }                                            << 
 90   // Generate first the forwrad primaries      << 
 91   theAdjointPrimaryGenerator->GenerateFwdPrima << 
 92     anEvent, ListOfPrimaryFwdParticles[index_p << 
 93   G4PrimaryVertex* fwdPrimVertex = anEvent->Ge << 
 94                                                << 
 95   p = fwdPrimVertex->GetPrimary()->GetMomentum << 
 96   pos = fwdPrimVertex->GetPosition();          << 
 97   G4double pmag = p.mag();                     << 
 98   G4double m0 = ListOfPrimaryFwdParticles[inde << 
 99   G4double ekin = std::sqrt(m0 * m0 + pmag * p << 
100                                                << 
101   G4double weight_correction = 1.;             << 
102   // For gamma generate the particle along the << 
103   G4ThreeVector dir = -p / p.mag();            << 
104                                                << 
105   weight_correction = 1.;                      << 
106                                                << 
107   if (ListOfPrimaryFwdParticles[index_particle << 
108   {                                            << 
109     G4double weight = (1. / nb_fwd_gammas_per_ << 
110     fwdPrimVertex->SetWeight(weight);          << 
111     for (G4int i = 0; i < nb_fwd_gammas_per_ev << 
112       auto newFwdPrimVertex = new G4PrimaryVer << 
113       newFwdPrimVertex->SetPosition(pos.x(), p << 
114       newFwdPrimVertex->SetT0(0.);             << 
115       auto aPrimParticle =                     << 
116         new G4PrimaryParticle(ListOfPrimaryFwd << 
117       newFwdPrimVertex->SetPrimary(aPrimPartic << 
118       newFwdPrimVertex->SetWeight(weight);     << 
119       anEvent->AddPrimaryVertex(newFwdPrimVert << 
120     }                                          << 
121   }                                            << 
122                                                    79 
123   // Now generate the adjoint primaries        <<  80    G4int evt_id=anEvent->GetEventID();
124   auto adjPrimVertex = new G4PrimaryVertex();  <<  81    size_t n=ListOfPrimaryAdjParticles.size();
125   adjPrimVertex->SetPosition(pos.x(), pos.y(), <<  82    index_particle=size_t(evt_id)-n*(size_t(evt_id)/n);
126   adjPrimVertex->SetT0(0.);                    <<  83 
127   auto aPrimParticle =                         <<  84 
128     new G4PrimaryParticle(ListOfPrimaryAdjPart <<  85    G4double E1=Emin;
129                                                <<  86    G4double E2=Emax;
130   adjPrimVertex->SetPrimary(aPrimParticle);    <<  87    if (!ListOfPrimaryAdjParticles[index_particle]) UpdateListOfPrimaryParticles();//ion has not been created yet
131   anEvent->AddPrimaryVertex(adjPrimVertex);    <<  88 
132                                                <<  89    if (ListOfPrimaryAdjParticles[index_particle]->GetParticleName() == "adj_proton") {
133   // The factor pi is to normalise the weight  <<  90      E1=EminIon;
134   G4double adjoint_source_area = G4AdjointSimM <<  91      E2=EmaxIon;
135   G4double adjoint_weight =                    <<  92    }
136     weight_correction * ComputeEnergyDistWeigh <<  93    if (ListOfPrimaryAdjParticles[index_particle]->GetParticleType() == "adjoint_nucleus") {
137   if (ListOfPrimaryAdjParticles[index_particle <<  94      G4int A= ListOfPrimaryAdjParticles[index_particle]->GetAtomicMass();
138     // The weight will be corrected at the end <<  95      E1=EminIon*A;
139     // are used                                <<  96      E2=EmaxIon*A;
140     adjoint_weight = adjoint_weight / nb_adj_p <<  97    }
141     for (G4int i = 0; i < nb_adj_primary_gamma <<  98    //Generate first the forwrad primaries
142       auto newAdjPrimVertex = new G4PrimaryVer <<  99    theAdjointPrimaryGenerator->GenerateFwdPrimaryVertex(anEvent,ListOfPrimaryFwdParticles[index_particle],E1,E2);
143       newAdjPrimVertex->SetPosition(pos.x(), p << 100    G4PrimaryVertex* fwdPrimVertex = anEvent->GetPrimaryVertex();
144       newAdjPrimVertex->SetT0(0.);             << 101 
145       aPrimParticle =                          << 102    p=fwdPrimVertex->GetPrimary()->GetMomentum();
146         new G4PrimaryParticle(ListOfPrimaryAdj << 103    pos=fwdPrimVertex->GetPosition();
147       newAdjPrimVertex->SetPrimary(aPrimPartic << 104    G4double pmag=p.mag();
148       newAdjPrimVertex->SetWeight(adjoint_weig << 105    G4double m0=ListOfPrimaryFwdParticles[index_particle]->GetPDGMass();
149       anEvent->AddPrimaryVertex(newAdjPrimVert << 106    G4double ekin=std::sqrt( m0*m0 + pmag*pmag) -m0;
150     }                                          << 107 
151   }                                            << 108    G4double weight_correction=1.;
152   else if (ListOfPrimaryAdjParticles[index_par << 109    //For gamma generate the particle along the backward ray
153     // The weight will be corrected at the end << 110    G4ThreeVector dir=-p/p.mag();
154     // are used                                << 111 
155     adjoint_weight = adjoint_weight / nb_adj_p << 112    /*if (ListOfPrimaryAdjParticles[index_particle]->GetParticleName() == "adj_gamma"){
156     for (G4int i = 0; i < nb_adj_primary_elect << 113 
157       auto newAdjPrimVertex = new G4PrimaryVer << 114       theAdjointPrimaryGenerator
158       newAdjPrimVertex->SetPosition(pos.x(), p << 115               ->ComputeAccumulatedDepthVectorAlongBackRay(pos,dir,ekin,ListOfPrimaryAdjParticles[index_particle]);
159       newAdjPrimVertex->SetT0(0.);             << 116 
160       aPrimParticle =                          << 117 
161         new G4PrimaryParticle(ListOfPrimaryAdj << 118       G4double distance = theAdjointPrimaryGenerator
162       newAdjPrimVertex->SetPrimary(aPrimPartic << 119        ->SampleDistanceAlongBackRayAndComputeWeightCorrection(weight_correction);
163       newAdjPrimVertex->SetWeight(adjoint_weig << 120 
164       anEvent->AddPrimaryVertex(newAdjPrimVert << 121       //pos=pos+dir*distance;
165     }                                          << 122       //fwdPrimVertex->SetPosition(pos[0],pos[1],pos[2]);
166   }                                            << 123    }
167   adjPrimVertex->SetWeight(adjoint_weight);    << 124    */
                                                   >> 125    weight_correction=1.;
                                                   >> 126 
                                                   >> 127    if (ListOfPrimaryFwdParticles[index_particle]  ==G4Gamma::Gamma() && nb_fwd_gammas_per_event>1 ){
                                                   >> 128      G4double weight = (1./nb_fwd_gammas_per_event);
                                                   >> 129      fwdPrimVertex->SetWeight(weight);
                                                   >> 130      for (int i=0;i<nb_fwd_gammas_per_event-1;i++){
                                                   >> 131          G4PrimaryVertex* newFwdPrimVertex = new G4PrimaryVertex();
                                                   >> 132          newFwdPrimVertex->SetPosition(pos.x(),pos.y(),pos.z());
                                                   >> 133          newFwdPrimVertex->SetT0(0.);
                                                   >> 134          G4PrimaryParticle* aPrimParticle = new G4PrimaryParticle(ListOfPrimaryFwdParticles[index_particle],
                                                   >> 135                p.x(),p.y(),p.z());
                                                   >> 136          newFwdPrimVertex->SetPrimary(aPrimParticle);
                                                   >> 137          newFwdPrimVertex->SetWeight(weight);
                                                   >> 138          anEvent->AddPrimaryVertex(newFwdPrimVertex);
                                                   >> 139        }
                                                   >> 140    }
                                                   >> 141 
                                                   >> 142    //Now generate the adjoint primaries
                                                   >> 143    G4PrimaryVertex* adjPrimVertex = new G4PrimaryVertex();
                                                   >> 144    adjPrimVertex->SetPosition(pos.x(),pos.y(),pos.z());
                                                   >> 145    adjPrimVertex->SetT0(0.);
                                                   >> 146    G4PrimaryParticle* aPrimParticle = new G4PrimaryParticle(ListOfPrimaryAdjParticles[index_particle],
                                                   >> 147                          -p.x(),-p.y(),-p.z());
                                                   >> 148 
                                                   >> 149    adjPrimVertex->SetPrimary(aPrimParticle);
                                                   >> 150    anEvent->AddPrimaryVertex(adjPrimVertex);
                                                   >> 151 
                                                   >> 152    //The factor pi is to normalise the weight to the directional flux
                                                   >> 153    G4double adjoint_source_area = G4AdjointSimManager::GetInstance()->GetAdjointSourceArea();
                                                   >> 154    G4double adjoint_weight = weight_correction*ComputeEnergyDistWeight(ekin,E1,E2)*adjoint_source_area*pi;
                                                   >> 155    //if (ListOfPrimaryFwdParticles[index_particle]  ==G4Gamma::Gamma()) adjoint_weight = adjoint_weight/3.;
                                                   >> 156    if (ListOfPrimaryAdjParticles[index_particle]->GetParticleName() == "adj_gamma") {
                                                   >> 157      //The weight will be corrected at the end of the track if splitted tracks are used
                                                   >> 158      adjoint_weight = adjoint_weight/nb_adj_primary_gammas_per_event;
                                                   >> 159      for (int i=0;i<nb_adj_primary_gammas_per_event-1;i++){
                                                   >> 160           G4PrimaryVertex* newAdjPrimVertex = new G4PrimaryVertex();
                                                   >> 161         newAdjPrimVertex->SetPosition(pos.x(),pos.y(),pos.z());
                                                   >> 162         newAdjPrimVertex->SetT0(0.);
                                                   >> 163         aPrimParticle = new G4PrimaryParticle(ListOfPrimaryAdjParticles[index_particle],
                                                   >> 164                                 -p.x(),-p.y(),-p.z());
                                                   >> 165         newAdjPrimVertex->SetPrimary(aPrimParticle);
                                                   >> 166         newAdjPrimVertex->SetWeight(adjoint_weight);
                                                   >> 167         anEvent->AddPrimaryVertex(newAdjPrimVertex);
                                                   >> 168      }
                                                   >> 169    }
                                                   >> 170    else if  (ListOfPrimaryAdjParticles[index_particle]->GetParticleName() == "adj_electron") {
                                                   >> 171      //The weight will be corrected at the end of the track if splitted tracks are used
                                                   >> 172      adjoint_weight = adjoint_weight/nb_adj_primary_electrons_per_event;
                                                   >> 173      for (int i=0;i<nb_adj_primary_electrons_per_event-1;i++){
                                                   >> 174           G4PrimaryVertex* newAdjPrimVertex = new G4PrimaryVertex();
                                                   >> 175         newAdjPrimVertex->SetPosition(pos.x(),pos.y(),pos.z());
                                                   >> 176         newAdjPrimVertex->SetT0(0.);
                                                   >> 177         aPrimParticle = new G4PrimaryParticle(ListOfPrimaryAdjParticles[index_particle],
                                                   >> 178                                 -p.x(),-p.y(),-p.z());
                                                   >> 179         newAdjPrimVertex->SetPrimary(aPrimParticle);
                                                   >> 180         newAdjPrimVertex->SetWeight(adjoint_weight);
                                                   >> 181         anEvent->AddPrimaryVertex(newAdjPrimVertex);
                                                   >> 182      }
                                                   >> 183    }
                                                   >> 184    adjPrimVertex->SetWeight(adjoint_weight);
                                                   >> 185 
                                                   >> 186    //Call some methods of G4AdjointSimManager
                                                   >> 187    G4AdjointSimManager::GetInstance()->SetAdjointTrackingMode(true);
                                                   >> 188    G4AdjointSimManager::GetInstance()->ClearEndOfAdjointTrackInfoVectors();
                                                   >> 189    G4AdjointSimManager::GetInstance()->ResetDidOneAdjPartReachExtSourceDuringEvent();
                                                   >> 190 
                                                   >> 191   /* if ( !last_generated_part_was_adjoint ) {
                                                   >> 192    
                                                   >> 193    index_particle++;
                                                   >> 194    if (index_particle >= ListOfPrimaryAdjParticles.size()) index_particle =0;
                                                   >> 195    
                                                   >> 196   
                                                   >> 197    G4double E1=Emin;
                                                   >> 198    G4double E2=Emax;
                                                   >> 199    if (!ListOfPrimaryAdjParticles[index_particle]) UpdateListOfPrimaryParticles();//ion has not been created yet
                                                   >> 200   
                                                   >> 201    if (ListOfPrimaryAdjParticles[index_particle]->GetParticleName() == "adj_proton") {
                                                   >> 202     E1=EminIon;
                                                   >> 203     E2=EmaxIon;
                                                   >> 204    }
                                                   >> 205    if (ListOfPrimaryAdjParticles[index_particle]->GetParticleType() == "adjoint_nucleus") {
                                                   >> 206     G4int A= ListOfPrimaryAdjParticles[index_particle]->GetAtomicMass();
                                                   >> 207     E1=EminIon*A;
                                                   >> 208     E2=EmaxIon*A;
                                                   >> 209    }
                                                   >> 210    theAdjointPrimaryGenerator->GenerateAdjointPrimaryVertex(anEvent,
                                                   >> 211                   ListOfPrimaryAdjParticles[index_particle],
                                                   >> 212                   E1,E2);
                                                   >> 213    G4PrimaryVertex* aPrimVertex = anEvent->GetPrimaryVertex();
                                                   >> 214   
                                                   >> 215   
                                                   >> 216       p=aPrimVertex->GetPrimary()->GetMomentum();
                                                   >> 217     pos=aPrimVertex->GetPosition();
                                                   >> 218     G4double pmag=p.mag();
                                                   >> 219     
                                                   >> 220     G4double m0=ListOfPrimaryAdjParticles[index_particle]->GetPDGMass();
                                                   >> 221     G4double ekin=std::sqrt( m0*m0 + pmag*pmag) -m0;
                                                   >> 222     
                                                   >> 223   
                                                   >> 224       //The factor pi is to normalise the weight to the directional flux
                                                   >> 225     G4double adjoint_source_area = G4AdjointSimManager::GetInstance()->GetAdjointSourceArea();
                                                   >> 226       G4double adjoint_weight = ComputeEnergyDistWeight(ekin,E1,E2)*adjoint_source_area*pi;
                                                   >> 227 
                                                   >> 228       aPrimVertex->SetWeight(adjoint_weight);
                                                   >> 229 
                                                   >> 230     last_generated_part_was_adjoint =true;
                                                   >> 231     G4AdjointSimManager::GetInstance()->SetAdjointTrackingMode(true);
                                                   >> 232     G4AdjointSimManager::GetInstance()->RegisterAdjointPrimaryWeight(adjoint_weight);
                                                   >> 233    }
                                                   >> 234    else {
                                                   >> 235           //fwd particle equivalent to the last generated adjoint particle ios generated  
                                                   >> 236       G4PrimaryVertex* aPrimVertex = new G4PrimaryVertex();
                                                   >> 237     aPrimVertex->SetPosition(pos.x(),pos.y(),pos.z());
                                                   >> 238       aPrimVertex->SetT0(0.);
                                                   >> 239     G4PrimaryParticle* aPrimParticle = new G4PrimaryParticle(ListOfPrimaryFwdParticles[index_particle],
                                                   >> 240                    -p.x(),-p.y(),-p.z()); 
                                                   >> 241   
                                                   >> 242     aPrimVertex->SetPrimary(aPrimParticle);
                                                   >> 243     anEvent->AddPrimaryVertex(aPrimVertex);                
                                                   >> 244     last_generated_part_was_adjoint =false;
                                                   >> 245     G4AdjointSimManager::GetInstance()->SetAdjointTrackingMode(false);
                                                   >> 246   */
168                                                   247 
169   // Call some methods of G4AdjointSimManager  << 
170   G4AdjointSimManager::GetInstance()->SetAdjoi << 
171   G4AdjointSimManager::GetInstance()->ClearEnd << 
172   G4AdjointSimManager::GetInstance()->ResetDid << 
173 }                                                 248 }
174                                                << 249 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
175 // ------------------------------------------- << 
176 //                                                250 //
177 void G4AdjointPrimaryGeneratorAction::SetEmin(    251 void G4AdjointPrimaryGeneratorAction::SetEmin(G4double val)
178 {                                                 252 {
179   Emin = val;                                  << 253   Emin=val;
180   EminIon = val;                               << 254   EminIon=val;
181 }                                                 255 }
182                                                << 256 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
183 // ------------------------------------------- << 
184 //                                                257 //
185 void G4AdjointPrimaryGeneratorAction::SetEmax(    258 void G4AdjointPrimaryGeneratorAction::SetEmax(G4double val)
186 {                                                 259 {
187   Emax = val;                                  << 260   Emax=val;
188   EmaxIon = val;                               << 261   EmaxIon=val;
189 }                                                 262 }
190                                                << 263 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
191 // ------------------------------------------- << 
192 //                                                264 //
193 void G4AdjointPrimaryGeneratorAction::SetEminI    265 void G4AdjointPrimaryGeneratorAction::SetEminIon(G4double val)
194 {                                                 266 {
195   EminIon = val;                               << 267   EminIon=val;
196 }                                                 268 }
197                                                << 269 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
198 // ------------------------------------------- << 
199 //                                                270 //
200 void G4AdjointPrimaryGeneratorAction::SetEmaxI    271 void G4AdjointPrimaryGeneratorAction::SetEmaxIon(G4double val)
201 {                                                 272 {
202   EmaxIon = val;                               << 273   EmaxIon=val;
203 }                                                 274 }
204                                                << 275 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
205 // ------------------------------------------- << 
206 //                                                276 //
207 G4double G4AdjointPrimaryGeneratorAction::Comp << 277 G4double G4AdjointPrimaryGeneratorAction::ComputeEnergyDistWeight(G4double E ,G4double E1, G4double E2)
208                                                << 
209 {                                                 278 {
210   // We generate N numbers of primaries with   << 279  //  We generate N numbers of primaries with  a 1/E energy law distribution.
211   // We have therefore  an energy distribution << 280  //  We have therefore  an energy distribution function    
212   //   f(E)=C/E  (1)                           << 281  //     f(E)=C/E  (1)
213   // with   C a constant that is such that     << 282  //  with   C a constant that is such that
214   //   N=Integral(f(E),E1,E2)=C.std::log(E2/E1 << 283  //   N=Integral(f(E),E1,E2)=C.std::log(E2/E1)  (2)
215   // Therefore from (2) we get                 << 284  //  Therefore from (2) we get    
216   //   C=N/ std::log(E2/E1) (3)                << 285  //     C=N/ std::log(E2/E1) (3) 
217   // and                                       << 286  //  and        
218   //   f(E)=N/ std::log(E2/E1)/E (4)           << 287  //   f(E)=N/ std::log(E2/E1)/E (4)
219   // For the adjoint simulation we need a ener << 288  //For the adjoint simulation we need a energy distribution  f'(E)=1..
220   // To get that we need therefore to apply a  << 289  //To get that we need therefore to apply a weight to the primary
221   //   W=1/f(E)=E*std::log(E2/E1)/N            << 290  //   W=1/f(E)=E*std::log(E2/E1)/N  
222   //                                           << 291  //
223   return std::log(E2 / E1) * E / G4AdjointSimM << 292   return std::log(E2/E1)*E/G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun();
                                                   >> 293  
224 }                                                 294 }
225                                                << 295 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
226 // ------------------------------------------- << 
227 //                                                296 //
228 void G4AdjointPrimaryGeneratorAction::SetSpher << 297 void G4AdjointPrimaryGeneratorAction::SetSphericalAdjointPrimarySource(G4double radius, G4ThreeVector center_pos)
229                                                << 298 { 
230 {                                              << 
231   radius_spherical_source = radius;               299   radius_spherical_source = radius;
232   center_spherical_source = center_pos;           300   center_spherical_source = center_pos;
233   type_of_adjoint_source = "Spherical";        << 301   type_of_adjoint_source ="Spherical";
234   theAdjointPrimaryGenerator->SetSphericalAdjo << 302   theAdjointPrimaryGenerator->SetSphericalAdjointPrimarySource(radius,center_pos);
235 }                                                 303 }
236                                                << 304 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
237 // ------------------------------------------- << 
238 //                                                305 //
239 void G4AdjointPrimaryGeneratorAction::SetAdjoi << 306 void G4AdjointPrimaryGeneratorAction::SetAdjointPrimarySourceOnAnExtSurfaceOfAVolume(const G4String& volume_name)
240   const G4String& volume_name)                 << 
241 {                                                 307 {
242   type_of_adjoint_source = "ExternalSurfaceOfA << 308   type_of_adjoint_source ="ExternalSurfaceOfAVolume";
243   theAdjointPrimaryGenerator->SetAdjointPrimar    309   theAdjointPrimaryGenerator->SetAdjointPrimarySourceOnAnExtSurfaceOfAVolume(volume_name);
244 }                                                 310 }
245                                                << 311 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
246 // ------------------------------------------- << 
247 //                                                312 //
248 void G4AdjointPrimaryGeneratorAction::Consider    313 void G4AdjointPrimaryGeneratorAction::ConsiderParticleAsPrimary(const G4String& particle_name)
249 {                                                 314 {
250   if (PrimariesConsideredInAdjointSim.find(par << 315   if (PrimariesConsideredInAdjointSim.find(particle_name) != PrimariesConsideredInAdjointSim.end()){
251   {                                            << 316     PrimariesConsideredInAdjointSim[particle_name]=true;
252     PrimariesConsideredInAdjointSim[particle_n << 
253   }                                               317   }
254   UpdateListOfPrimaryParticles();                 318   UpdateListOfPrimaryParticles();
255 }                                                 319 }
256                                                << 320 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
257 // ------------------------------------------- << 
258 //                                                321 //
259 void G4AdjointPrimaryGeneratorAction::NeglectP    322 void G4AdjointPrimaryGeneratorAction::NeglectParticleAsPrimary(const G4String& particle_name)
260 {                                                 323 {
261   if (PrimariesConsideredInAdjointSim.find(par << 324   if (PrimariesConsideredInAdjointSim.find(particle_name) != PrimariesConsideredInAdjointSim.end()){
262   {                                            << 325     PrimariesConsideredInAdjointSim[particle_name]= false;
263     PrimariesConsideredInAdjointSim[particle_n << 
264   }                                               326   }
265   UpdateListOfPrimaryParticles();                 327   UpdateListOfPrimaryParticles();
266 }                                                 328 }
267                                                << 329 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
268 // ------------------------------------------- << 
269 //                                                330 //
270 void G4AdjointPrimaryGeneratorAction::UpdateLi    331 void G4AdjointPrimaryGeneratorAction::UpdateListOfPrimaryParticles()
271 {                                                 332 {
272   G4ParticleTable* theParticleTable = G4Partic << 333     G4ParticleTable* theParticleTable = G4ParticleTable::GetParticleTable();
273   ListOfPrimaryFwdParticles.clear();           << 334     ListOfPrimaryFwdParticles.clear();
274   ListOfPrimaryAdjParticles.clear();           << 335     ListOfPrimaryAdjParticles.clear();
275   for (const auto& iter : PrimariesConsideredI << 336     std::map<G4String, G4bool>::iterator iter;
276     if (iter.second) {                         << 337     for( iter = PrimariesConsideredInAdjointSim.begin(); iter != PrimariesConsideredInAdjointSim.end(); ++iter ) {
277       G4String fwd_particle_name = iter.first; << 338   if(iter->second) {
278       if (fwd_particle_name != "ion") {        << 339     G4String fwd_particle_name = iter->first;
279         G4String adj_particle_name = G4String( << 340     if ( fwd_particle_name != "ion") {
280         ListOfPrimaryFwdParticles.push_back(th << 341       G4String adj_particle_name = G4String("adj_") + fwd_particle_name;
281         ListOfPrimaryAdjParticles.push_back(th << 342       ListOfPrimaryFwdParticles.push_back(theParticleTable->FindParticle(fwd_particle_name));
282       }                                        << 343       ListOfPrimaryAdjParticles.push_back(theParticleTable->FindParticle(adj_particle_name));
283       else {                                   << 344       /*
284         if (fwd_ion != nullptr) {              << 345       if ( fwd_particle_name == "gamma") {
285           ion_name = fwd_ion->GetParticleName( << 346         for (G4int i=0;i<2;i++){
286           G4String adj_ion_name = G4String("ad << 347           ListOfPrimaryFwdParticles.push_back(theParticleTable->FindParticle(fwd_particle_name));
287           ListOfPrimaryFwdParticles.push_back( << 348           ListOfPrimaryAdjParticles.push_back(theParticleTable->FindParticle(adj_particle_name));
288           ListOfPrimaryAdjParticles.push_back( << 349         }
289         }                                      << 350       }
290         else {                                 << 351       */
291           ListOfPrimaryFwdParticles.push_back( << 352     }
292           ListOfPrimaryAdjParticles.push_back( << 353     else {
293         }                                      << 354       if (fwd_ion ){
294       }                                        << 355         ion_name=fwd_ion->GetParticleName();
295     }                                          << 356         G4String adj_ion_name=G4String("adj_") +ion_name;
296   }                                            << 357         ListOfPrimaryFwdParticles.push_back(fwd_ion);
                                                   >> 358         ListOfPrimaryAdjParticles.push_back(adj_ion);
                                                   >> 359       }
                                                   >> 360       else {
                                                   >> 361         ListOfPrimaryFwdParticles.push_back(0);
                                                   >> 362         ListOfPrimaryAdjParticles.push_back(0);
                                                   >> 363         
                                                   >> 364       } 
                                                   >> 365     }
                                                   >> 366   } 
                                                   >> 367    }
297 }                                                 368 }
298                                                << 369 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
299 // ------------------------------------------- << 
300 //                                                370 //
301 void G4AdjointPrimaryGeneratorAction::SetPrima << 371 void G4AdjointPrimaryGeneratorAction::SetPrimaryIon(G4ParticleDefinition* adjointIon, G4ParticleDefinition* fwdIon)
302                                                << 
303 {                                                 372 {
304   fwd_ion = fwdIon;                               373   fwd_ion = fwdIon;
305   adj_ion = adjointIon;                           374   adj_ion = adjointIon;
306   UpdateListOfPrimaryParticles();                 375   UpdateListOfPrimaryParticles();
307 }                                                 376 }
                                                   >> 377 
308                                                   378