G4NeutronHPFinalState.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 //
00027 //080721 Create adjust_final_state method by T. Koi
00028 //080801 Residual reconstruction with theNDLDataA,Z (A, Z, and momentum are adjusted) by T. Koi
00029 //101110 Set lower limit for gamma energy(1keV) by T. Koi
00030 
00031 #include "G4NeutronHPFinalState.hh"
00032 #include "G4PhysicalConstants.hh"
00033 #include "G4SystemOfUnits.hh"
00034 #include "G4ParticleTable.hh"
00035 #include "G4Gamma.hh"
00036 #include "G4Neutron.hh"
00037 
00038 void G4NeutronHPFinalState::adjust_final_state ( G4LorentzVector init_4p_lab )
00039 {
00040 
00041    G4double minimum_energy = 1*keV;
00042 
00043    if ( adjustResult != true ) return;
00044 
00045    G4int nSecondaries = theResult.GetNumberOfSecondaries();
00046 
00047    G4int sum_Z = 0;
00048    G4int sum_A = 0;
00049    G4int max_SecZ = 0;
00050    G4int max_SecA = 0;
00051    G4int imaxA = -1;
00052    for ( int i = 0 ; i < nSecondaries ; i++ )
00053    {
00054       sum_Z += theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicNumber();
00055       max_SecZ = std::max ( max_SecZ , theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicNumber() );
00056       sum_A += theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicMass();
00057       max_SecA = std::max ( max_SecA , theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicMass() );
00058       if ( theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicMass() == max_SecA ) imaxA = i;
00059    }
00060 
00061    G4ParticleDefinition* resi_pd = NULL;
00062    G4bool needOneMoreSec = false;
00063    G4ParticleDefinition* oneMoreSec_pd = NULL;
00064    if ( (int)(theBaseZ - sum_Z) == 0 && (int)(theBaseA + 1 - sum_A) == 0 )
00065    {
00066       //All secondaries are already created;
00067       resi_pd = theResult.GetSecondary( imaxA )->GetParticle()->GetDefinition(); 
00068    }
00069    else
00070    {
00071       if ( max_SecA > int(theBaseA + 1 - sum_A) ) 
00072       {
00073          //Most heavy secondary is interpreted as residual
00074          resi_pd = theResult.GetSecondary( imaxA )->GetParticle()->GetDefinition(); 
00075          needOneMoreSec = true;
00076       }
00077       else 
00078       {
00079          //creation of residual is requierd 
00080          resi_pd = G4ParticleTable::GetParticleTable()->GetIon ( int(theBaseZ - sum_Z) , (int)(theBaseA + 1 - sum_A) , 0.0 );
00081       }
00082 
00083       if ( needOneMoreSec )
00084       {
00085          if ( int(theBaseZ - sum_Z) == 0 && (int)(theBaseA + 1 - sum_A) > 0 )
00086          {
00087             if ( int(theBaseA + 1 - sum_A) > 1 ) G4cout << "More than one neutron is required for the balance of baryon number!" << G4endl;
00088             oneMoreSec_pd = G4Neutron::Neutron();
00089          }
00090          else 
00091          {
00092             oneMoreSec_pd = G4ParticleTable::GetParticleTable()->GetIon ( int(theBaseZ - sum_Z) , (int)(theBaseA + 1 - sum_A) , 0.0 );
00093          }
00094       }
00095   
00096       if ( resi_pd == NULL )
00097       {
00098          // theNDLDataZ,A has the Z and A of used NDL file
00099          if ( (int)(theNDLDataZ - sum_Z) == 0 && (int)(theNDLDataA + 1 - sum_A) == 0 )
00100          {
00101             G4int dif_Z = ( int ) ( theNDLDataZ - theBaseZ );
00102             G4int dif_A = ( int ) ( theNDLDataA - theBaseA );
00103             resi_pd = G4ParticleTable::GetParticleTable()->GetIon ( max_SecZ - dif_Z , max_SecA - dif_A , 0.0 );
00104             for ( int i = 0 ; i < nSecondaries ; i++ )
00105             {
00106                if ( theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicNumber() == max_SecZ  
00107                  && theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetAtomicMass() == max_SecA )
00108                {
00109                   G4ThreeVector p = theResult.GetSecondary( i )->GetParticle()->GetMomentum();
00110                   p = p * resi_pd->GetPDGMass()/ G4ParticleTable::GetParticleTable()->GetIon ( max_SecZ , max_SecA , 0.0 )->GetPDGMass(); 
00111                   theResult.GetSecondary( i )->GetParticle()->SetDefinition( resi_pd );
00112                   theResult.GetSecondary( i )->GetParticle()->SetMomentum( p );
00113                } 
00114             }
00115          }
00116       }
00117    }
00118 
00119 
00120    G4LorentzVector secs_4p_lab( 0.0 );
00121 
00122    G4int n_sec = theResult.GetNumberOfSecondaries();
00123    G4double fast = 0;
00124    G4double slow = 1;
00125    G4int ifast = 0;
00126    G4int islow = 0;
00127    G4int ires = -1;
00128 
00129    for ( G4int i = 0 ; i < n_sec ; i++ )
00130    {
00131 
00132       //G4cout << "HP_DB " << i 
00133       //       << " " << theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetParticleName() 
00134       //       << " 4p " << theResult.GetSecondary( i )->GetParticle()->Get4Momentum() 
00135       //       << " ke " << theResult.GetSecondary( i )->GetParticle()->Get4Momentum().e() - theResult.GetSecondary( i )->GetParticle()->GetDefinition()->GetPDGMass() 
00136       //       << G4endl;
00137 
00138       secs_4p_lab += theResult.GetSecondary( i )->GetParticle()->Get4Momentum();
00139 
00140       G4double beta = 0;
00141       if ( theResult.GetSecondary( i )->GetParticle()->GetDefinition() != G4Gamma::Gamma() )
00142       {
00143          beta = theResult.GetSecondary( i )->GetParticle()->Get4Momentum().beta(); 
00144       }
00145       else
00146       {
00147          beta = 1;
00148       }
00149 
00150       if ( theResult.GetSecondary( i )->GetParticle()->GetDefinition() == resi_pd ) ires = i; 
00151 
00152       if ( slow > beta && beta != 0 ) 
00153       {
00154          slow = beta; 
00155          islow = i;
00156       }
00157 
00158       if ( fast <= beta ) 
00159       {
00160          if ( fast != 1 )
00161          {
00162             fast = beta; 
00163             ifast = i;
00164          } 
00165          else
00166          {
00167 //          fast is already photon then check E
00168             G4double e = theResult.GetSecondary( i )->GetParticle()->Get4Momentum().e();
00169             if ( e > theResult.GetSecondary( ifast )->GetParticle()->Get4Momentum().e() )   
00170             {
00171 //             among photons, the highest E becomes the fastest 
00172                ifast = i; 
00173             }
00174          } 
00175       }
00176    }
00177 
00178 
00179    G4LorentzVector dif_4p = init_4p_lab - secs_4p_lab;
00180 
00181    //G4cout << "HP_DB dif_4p " << init_4p_lab - secs_4p_lab << G4endl;
00182    //G4cout << "HP_DB dif_3p mag " << ( dif_4p.v() ).mag() << G4endl;
00183    //G4cout << "HP_DB dif_e " << dif_4p.e() - ( dif_4p.v() ).mag()<< G4endl;
00184 
00185    G4LorentzVector p4(0);
00186    if ( ires == -1 ) 
00187    {
00188 //    Create and Add Residual Nucleus 
00189       ires = nSecondaries;
00190       nSecondaries += 1;
00191 
00192       G4DynamicParticle* res = new G4DynamicParticle ( resi_pd , dif_4p.v() );    
00193       theResult.AddSecondary ( res );    
00194 
00195       p4 = res->Get4Momentum(); 
00196       if ( slow > p4.beta() ) 
00197       {
00198          slow = p4.beta(); 
00199          islow = ires;
00200       }
00201       dif_4p = init_4p_lab - ( secs_4p_lab + p4 );  
00202    }
00203 
00204    if ( needOneMoreSec )
00205    {
00206       nSecondaries += 1;
00207       G4DynamicParticle* one = new G4DynamicParticle ( oneMoreSec_pd , dif_4p.v() );    
00208       theResult.AddSecondary ( one );    
00209       p4 = one->Get4Momentum(); 
00210       if ( slow > p4.beta() ) 
00211       {
00212          slow = p4.beta(); 
00213          islow = nSecondaries-1; //Because the first is 0th, so the last becomes "nSecondaries-1"
00214       }
00215       dif_4p = init_4p_lab - ( secs_4p_lab + p4 );  
00216    }
00217 
00218    //Which is bigger dif_p or dif_e
00219 
00220    if ( dif_4p.v().mag() < std::abs( dif_4p.e() ) )
00221    {
00222 
00223       // Adjust p
00224       //if ( dif_4p.v().mag() < 1*MeV )
00225       if ( minimum_energy < dif_4p.v().mag() && dif_4p.v().mag() < 1*MeV )
00226       {
00227 
00228          nSecondaries += 1;
00229          theResult.AddSecondary ( new G4DynamicParticle ( G4Gamma::Gamma() , dif_4p.v() ) );    
00230 
00231       }
00232       else
00233       {
00234          //G4cout << "HP_DB Difference in dif_p is too large (>1MeV) or too small(<1keV) to adjust, so that give up tuning" << G4endl;
00235       }
00236 
00237    }
00238    else
00239    {
00240 
00241       // dif_p > dif_e 
00242       // at first momentum 
00243       // Move residual momentum
00244 
00245       p4 = theResult.GetSecondary( ires )->GetParticle()->Get4Momentum();
00246       theResult.GetSecondary( ires )->GetParticle()->SetMomentum( p4.v() + dif_4p.v() ); 
00247       dif_4p = init_4p_lab - ( secs_4p_lab - p4 + theResult.GetSecondary( ires )->GetParticle()->Get4Momentum()  );  
00248 
00249       //G4cout << "HP_DB new residual kinetic energy " << theResult.GetSecondary( ires )->GetParticle()->GetKineticEnergy() << G4endl;
00250 
00251    }
00252 
00253    G4double dif_e = dif_4p.e() - ( dif_4p.v() ).mag();
00254    //G4cout << "HP_DB dif_e " << dif_e << G4endl;
00255 
00256    if ( dif_e > 0 )
00257    {
00258 
00259 //    create 2 gamma 
00260 
00261       nSecondaries += 2;
00262       G4double e1 = ( dif_4p.e() -dif_4p.v().mag() ) / 2;
00263       
00264       if ( minimum_energy < e1 )  
00265       {
00266          G4double costh = 2.*G4UniformRand()-1.;
00267          G4double phi = twopi*G4UniformRand();
00268          G4ThreeVector dir( std::sin(std::acos(costh))*std::cos(phi), 
00269                             std::sin(std::acos(costh))*std::sin(phi),
00270                             costh);
00271          theResult.AddSecondary ( new G4DynamicParticle ( G4Gamma::Gamma() , e1*dir ) );    
00272          theResult.AddSecondary ( new G4DynamicParticle ( G4Gamma::Gamma() , -e1*dir ) );    
00273       }
00274       else
00275       {
00276          //G4cout << "HP_DB Difference is too small(<1keV) to adjust, so that neglect it" << G4endl;
00277       }
00278 
00279    }
00280    else    //dif_e < 0
00281    {
00282 
00283 //    At first reduce KE of the fastest secondary;
00284       G4double ke0 = theResult.GetSecondary( ifast )->GetParticle()->GetKineticEnergy();
00285       G4ThreeVector p0 = theResult.GetSecondary( ifast )->GetParticle()->GetMomentum();
00286       G4ThreeVector dir = ( theResult.GetSecondary( ifast )->GetParticle()->GetMomentum() ).unit();
00287 
00288       //G4cout << "HP_DB ifast " << ifast << " ke0 " << ke0 << G4endl;
00289 
00290       if ( ke0 + dif_e > 0 )
00291       {
00292          theResult.GetSecondary( ifast )->GetParticle()->SetKineticEnergy( ke0 + dif_e ); 
00293          G4ThreeVector dp = p0 - theResult.GetSecondary( ifast )->GetParticle()->GetMomentum();  
00294 
00295          G4ThreeVector p = theResult.GetSecondary( islow )->GetParticle()->GetMomentum();
00296          //theResult.GetSecondary( islow )->GetParticle()->SetMomentum( p - dif_e*dir ); 
00297          theResult.GetSecondary( islow )->GetParticle()->SetMomentum( p + dp ); 
00298       }
00299       else
00300       {
00301          //G4cout << "HP_DB Difference in dif_e too large ( <0MeV ) to adjust, so that give up tuning" << G4endl;
00302       }
00303 
00304    }
00305 
00306 }

Generated on Mon May 27 17:49:01 2013 for Geant4 by  doxygen 1.4.7