Ignore:
Timestamp:
Apr 6, 2009, 12:30:29 PM (17 years ago)
Author:
garnier
Message:

update processes

Location:
trunk/source/processes/hadronic/models/incl/include
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/models/incl/include/G4Abla.hh

    r819 r962  
    2424// ********************************************************************
    2525//
    26 // $Id: G4Abla.hh,v 1.7 2007/12/03 19:36:05 miheikki Exp $
     26// $Id: G4Abla.hh,v 1.11 2008/06/25 17:20:03 kaitanie Exp $
    2727// Translation of INCL4.2/ABLA V3
    2828// Pekka Kaitaniemi, HIP (translation)
    … …  
    3333#include "globals.hh"
    3434
     35#include "G4InclRandomNumbers.hh"
    3536#include "G4AblaDataDefs.hh"
    3637#include "G4InclDataDefs.hh"
     38#include "G4AblaFissionBase.hh"
     39
     40#ifndef G4Abla_hh
     41#define G4Abla_hh 1
    3742
    3843/**
    … …  
    4449public:
    4550  /**
    46    *   
    47    * We support Doxygen with JavaDoc style.
    48    *
    49    * \author{pekka.kaitaniemi@helsinki.fi}
    50    */
    51 
    52   /**
    53    * A constructor.
     51   * Basic constructor.
    5452   */
    5553  G4Abla();
    … …  
    6563
    6664  /**
    67    *
     65   * Constructor that is to be used only for testing purposes.
     66   * @param aHazard random seeds
     67   * @param aVolant data structure for ABLA output   
    6868   */
    6969  G4Abla(G4Hazard *hazard, G4Volant *volant);
    7070
    7171  /**
    72    * A destructor.
    73    * A more elaborate description of the destructor.
     72   * Basic destructor.
    7473   */
    7574  ~G4Abla();
    7675
    7776  /**
    78    *
    79    */
    80   void setVerboseLevel(G4int level);
     77   * Set verbosity level.
     78   */
     79  void setVerboseLevel(G4int level) {
     80    verboseLevel = level;
     81  }
     82
     83  /**
     84   * Get the internal output data structure pointer.
     85   */
     86  G4Volant* getVolant() {
     87    return volant;
     88  }
    8189
    8290  /**
    … …  
    102110  /**
    103111   * Initialize ABLA evaporation code.
     112   *
    104113   */
    105114  void initEvapora();
    106115
    107116  /**
    108    * qrot including damping                                                
     117   * Coefficient of collective enhancement including damping                         
    109118   * Input: z,a,bet,sig,u                                                 
    110119   * Output: qr - collective enhancement factor                           
    111120   * See  junghans et al., nucl. phys. a 629 (1998) 635                   
     121   * @param z charge number
     122   * @param a mass number
     123   * @param bet beta deformation
     124   * @param sig perpendicular spin cut-off factor
     125   * @param u Energy
     126   * @return Coefficient of collective enhancement   
    112127   */
    113128  void qrot(G4double z, G4double a, G4double bet, G4double sig, G4double u, G4double *qr);
    … …  
    127142   *
    128143   */
    129   //  G4double spdef(G4int a, G4int z, G4int optxfis);
     144  G4double spdef(G4int a, G4int z, G4int optxfis);
    130145
    131146  /**
    132147   * Calculation of fissility parameter
    133148   */
    134   //  G4double fissility(int a,int z, int optxfis);
     149  G4double fissility(int a,int z, int optxfis);
    135150
    136151  /**
    … …  
    224239
    225240  /**
     241   * Random numbers.
     242   */
     243  G4double haz(G4int k);
     244  void standardRandom(G4double *rndm, G4long *seed);
     245
     246  /**
    226247   * TIRAGE ALEATOIRE DANS UNE EXPONENTIELLLE : Y=EXP(-X/T)
    227248   */
    … …  
    252273   */
    253274  void guet(G4double *x_par, G4double *z_par, G4double *find_par);
    254 
    255   // Fission
    256 public:
    257   /**
    258    *
    259    */
    260   G4double spdef(G4int a, G4int z, G4int optxfis);
    261 
    262   /**
    263    *
    264    */
    265   G4double fissility(G4int a, G4int z, G4int optxfis);
    266 
    267 //   void evapora(G4double zprf, G4double aprf, G4double ee, G4double jprf,
    268 //             G4double *zf_par, G4double *af_par, G4double *mtota_par,
    269 //             G4double *pleva_par, G4double *pxeva_par);
    270 //  G4double bfms67(G4double zms, G4double ams);
    271   //  void lpoly(G4double x, G4int n, G4double pl[]);
    272   //  G4double expohaz(G4int k, G4double T);
    273   //  G4double fd(G4double E);
    274   //  G4double f(G4double E);
    275   //  G4double fmaxhaz(G4double k, G4double T);
    276   void even_odd(G4double r_origin,G4double r_even_odd,G4int &i_out);
    277   G4double umass(G4double z,G4double n,G4double beta);
    278   G4double ecoul(G4double z1,G4double n1,G4double beta1,G4double z2,G4double n2,G4double beta2,G4double d);
    279   void fissionDistri(G4double &a,G4double &z,G4double &e,
    280                      G4double &a1,G4double &z1,G4double &e1,G4double &v1,
    281                      G4double &a2,G4double &z2,G4double &e2,G4double &v2);
    282   void standardRandom(G4double *rndm, G4long *seed);
    283   G4double haz(G4int k);
    284   G4double gausshaz(int k, double xmoy, double sig);
    285 
    286275   
    287276public:
    … …  
    312301  G4double utilabs(G4double a);
    313302  G4double dmin1(G4double a, G4double b, G4double c);
     303  G4Ec2sub* getFrldmTable() {
     304    return ec2sub;
     305  }
    314306
    315307private:
    316308  G4int verboseLevel;
    317 
     309  G4int ilast;
     310
     311  G4AblaFissionBase *fissionModel;
     312  G4InclRandomInterface *randomGenerator;
    318313  G4Pace *pace;
    319314  G4Hazard *hazard;
    320315  G4Ald *ald;
    321   G4Ablamain *ablamain;
    322   G4Emdpar *emdpar;
    323316  G4Eenuc *eenuc;
    324317  G4Ec2sub *ec2sub;
    … …  
    330323  G4VarNtp *varntp; 
    331324};
     325
     326#endif

  • trunk/source/processes/hadronic/models/incl/include/G4AblaDataDefs.hh

    r819 r962  
    2424// ********************************************************************
    2525//
    26 // $Id: G4AblaDataDefs.hh,v 1.6 2007/12/03 19:36:06 miheikki Exp $
     26// $Id: G4AblaDataDefs.hh,v 1.9 2008/06/25 17:20:04 kaitanie Exp $
    2727// Translation of INCL4.2/ABLA V3
    2828// Pekka Kaitaniemi, HIP (translation)
    … …  
    7676
    7777  G4double ecnz[EC2SUBROWS][EC2SUBCOLS];
     78
     79  /**
     80   * Dump the contents of the ecnz data table.
     81   */
     82  void dump() {
     83    for(G4int i = 0; i < EC2SUBROWS; i++) {
     84      for(G4int j = 0; j < EC2SUBCOLS; j++) {
     85        G4cout << ecnz[i][j] << " ";
     86      }
     87      G4cout << G4endl;
     88    }
     89  }
    7890};
    7991
    … …  
    8799 
    88100  G4double av,as,ak,optafan;
    89 };
    90 
    91 class G4Ablamain {
    92 public:
    93   G4Ablamain() {};
    94   ~G4Ablamain() {};
    95  
    96   G4double ap,zp,at,zt,eap,beta,bmaxnuc,crtot,crnuc,r_0, r_p,r_t,pi,bfpro,snpro,sppro,shell;
    97   G4int imax, inum;
    98101};
    99102
    … …  
    185188};
    186189
    187 #define EMDPARSIZE 1000
    188 /**
    189  * Energies widths and cross sections for em excitation.
    190  */
    191 
    192 class G4Emdpar {
    193 
    194 public:
    195   G4Emdpar() {};
    196   ~G4Emdpar() {};
    197  
    198   G4double egdr,egqr,fwhmgdr,fwhmgqr,cremde1,cremde2;                 
    199   G4double ae1[EMDPARSIZE],be1[EMDPARSIZE],ce1[EMDPARSIZE],ae2[EMDPARSIZE];     
    200   G4double be2[EMDPARSIZE],ce2[EMDPARSIZE],sre1[EMDPARSIZE],sre2[EMDPARSIZE];
    201   G4double xre1[EMDPARSIZE],xre2[EMDPARSIZE],ds1,ds2;                             
    202 };
    203 
    204190//#define VOLANTSIZE 200
    205191#define VOLANTSIZE 2000
    … …  
    211197 
    212198public:
    213   G4Volant() {};
     199  G4Volant()
     200  {
     201    clear();
     202  }
     203
    214204  ~G4Volant() {};
    215205
     206  void clear()
     207  {
     208    for(G4int i = 0; i < VOLANTSIZE; i++) {
     209      copied[i] = false;
     210      acv[i] = 0;
     211      zpcv[i] = 0;
     212      pcv[i] = 0;
     213      xcv[i] = 0;
     214      ycv[i] = 0;
     215      zcv[i] = 0;
     216      iv = 0;
     217    }
     218  }
     219
     220  G4double getTotalMass()
     221  {
     222    G4double total = 0.0;
     223    for(G4int i = 0; i <= iv; i++) {
     224      total += acv[i];
     225    }
     226    return total;
     227  }
     228
    216229  void dump()
    217230  {
     231    G4double totA = 0.0, totZ = 0.0, totP = 0.0;
    218232    G4cout <<"i \t ACV \t ZPCV \t PCV" << G4endl;
    219233    for(G4int i = 0; i <= iv; i++) {
     234      if(i == 0 && acv[i] != 0) {
     235        G4cout <<"G4Volant: Particle stored at index " << i << G4endl;
     236      }
     237      totA += acv[i];
     238      totZ += zpcv[i];
     239      totP += pcv[i];
    220240      G4cout << "volant" << i << "\t" << acv[i] << " \t " << zpcv[i] << " \t " << pcv[i] << G4endl;
    221241    }
     242    G4cout <<"Particle count index (iv) = " << iv << G4endl;
     243    G4cout <<"ABLA Total: A = " << totA << " Z = " << totZ <<  " momentum = " << totP << G4endl;
    222244  }
    223245
    224246  G4double acv[VOLANTSIZE],zpcv[VOLANTSIZE],pcv[VOLANTSIZE],xcv[VOLANTSIZE];
    225247  G4double ycv[VOLANTSIZE],zcv[VOLANTSIZE];
     248  G4bool copied[VOLANTSIZE];
    226249  G4int iv;
    227250};

  • trunk/source/processes/hadronic/models/incl/include/G4Incl.hh

    r819 r962  
    2424// ********************************************************************
    2525//
    26 // $Id: G4Incl.hh,v 1.10 2007/12/03 19:36:06 miheikki Exp $
     26// $Id: G4Incl.hh,v 1.13 2008/06/25 17:20:04 kaitanie Exp $
    2727// Translation of INCL4.2/ABLA V3
    2828// Pekka Kaitaniemi, HIP (translation)
    … …  
    3535
    3636#include "globals.hh"
     37#include "G4InclRandomNumbers.hh"
    3738#include "G4InclDataDefs.hh"
    3839#include "G4Abla.hh"
    … …  
    816817  G4Volant *volant;
    817818  G4Abla *abla;
     819  G4InclRandomInterface *randomGenerator;
    818820};
    819821

  • trunk/source/processes/hadronic/models/incl/include/G4InclAblaVirtualData.hh

    r819 r962  
    2424// ********************************************************************
    2525//
    26 // $Id: G4InclAblaVirtualData.hh,v 1.3 2007/12/03 19:36:06 miheikki Exp $
     26// $Id: G4InclAblaVirtualData.hh,v 1.5 2008/06/25 17:20:04 kaitanie Exp $
    2727// Translation of INCL4.2/ABLA V3
    2828// Pekka Kaitaniemi, HIP (translation)
    … …  
    9999private:
    100100
    101   static const G4int alphaRows = 155;
    102   static const G4int alphaCols = 100;
     101  static const G4int alphaRows = 154;
     102  static const G4int alphaCols = 99;
    103103
    104104  static const G4int paceRows = 500;

  • trunk/source/processes/hadronic/models/incl/include/G4InclDataDefs.hh

    r819 r962  
    2424// ********************************************************************
    2525//
    26 // $Id: G4InclDataDefs.hh,v 1.2 2007/09/11 13:18:43 miheikki Exp $
     26// $Id: G4InclDataDefs.hh,v 1.5 2008/06/25 17:20:04 kaitanie Exp $
    2727// Translation of INCL4.2/ABLA V3
    2828// Pekka Kaitaniemi, HIP (translation)
    … …  
    564564
    565565  /**
     566   * Clear and initialize all variables and arrays.
     567   */
     568  void clear() {
     569    particleIndex = 0;
     570    projType = 0;
     571    projEnergy = 0.0;
     572    targetA = 0;
     573    targetZ = 0;
     574    massini = 0;
     575    mzini = 0;
     576    exini = 0;
     577    pcorem = 0;
     578    mcorem = 0;
     579    pxrem = 0;
     580    pyrem = 0;
     581    pzrem = 0;
     582    mulncasc = 0;
     583    mulnevap = 0;
     584    mulntot = 0;
     585    bimpact = 0.0;
     586    jremn = 0;
     587    kfis = 0;
     588    estfis = 0;
     589    izfis = 0;
     590    iafis = 0;
     591    ntrack = 0;
     592    for(G4int i = 0; i < VARNTPSIZE; i++) {
     593      itypcasc[i] = 0;
     594      avv[i] = 0;
     595      zvv[i] = 0;
     596      enerj[i] = 0.0;
     597      plab[i] = 0.0;
     598      tetlab[i] = 0.0;
     599      philab[i] = 0.0;
     600      full[i] = false;
     601    }
     602  }
     603
     604  void addParticle(G4double A, G4double Z, G4double E, G4double P, G4double theta, G4double phi) {
     605    if(full[particleIndex]) {
     606      G4cout <<"G4VarNtp: Error. Index i = " << particleIndex << " is already occupied by particle:" << G4endl;
     607      G4cout <<"A = " << avv[particleIndex] << " Z = " << zvv[particleIndex] << G4endl;
     608      G4cout <<"Tried to replace it with:" << G4endl;
     609      G4cout <<"A = " << Z << " Z = " << Z << G4endl;
     610    } else {
     611      avv[particleIndex] = (int) A;
     612      zvv[particleIndex] = (int) Z;
     613      enerj[particleIndex] = E;
     614      plab[particleIndex] = P;
     615      tetlab[particleIndex] = theta;
     616      philab[particleIndex] = phi;
     617      full[particleIndex] = true;
     618      ntrack = particleIndex + 1;
     619      particleIndex++;
     620    }
     621  }
     622
     623  /**
     624   * Baryon number conservation check.
     625   */
     626  G4int getTotalBaryonNumber() {
     627    G4int baryonNumber = 0;
     628    for(G4int i = 0; i < ntrack; i++) {
     629      if(avv[i] > 0) {
     630        baryonNumber += avv[i];
     631      }
     632    }
     633    return baryonNumber;
     634  }
     635
     636  /**
     637   * Return total energy.
     638   */
     639  G4double getTotalEnergy() {
     640    G4double energy = 0.0;
     641    for(G4int i = 0; i < ntrack; i++) {
     642      energy += std::sqrt(std::pow(plab[i], 2) + std::pow(getMass(i), 2)); // E^2 = p^2 + m^2
     643    }
     644
     645    return energy;
     646  }
     647
     648  /**
     649   * Return total three momentum.
     650   */
     651  G4double getTotalThreeMomentum() {
     652    G4double momentum = 0;
     653    for(G4int i = 0; i < ntrack; i++) {
     654      momentum += plab[i];
     655    }
     656    return momentum;
     657  }
     658
     659  G4double getMomentumSum() {
     660    G4double momentum = 0;
     661    for(G4int i = 0; i < ntrack; i++) {
     662      momentum += plab[i];
     663    }
     664    return momentum;
     665  }
     666
     667  G4double getMass(G4int particle) {
     668    const G4double protonMass = 938.272;
     669    const G4double neutronMass = 939.565;
     670    const G4double pionMass = 139.57;
     671
     672    G4double mass = 0.0;
     673    if(avv[particle] ==  1 && zvv[particle] ==  1) mass = protonMass;
     674    if(avv[particle] ==  1 && zvv[particle] ==  0) mass = neutronMass;
     675    if(avv[particle] == -1)                        mass = pionMass;
     676    if(avv[particle] > 1)
     677      mass = avv[particle] * protonMass + zvv[particle] * neutronMass;
     678    return mass;
     679  }
     680
     681  /**
     682   * Dump debugging output.
     683   */
     684  void dump() {
     685    G4int nProton = 0, nNeutron = 0;
     686    G4int nPiPlus = 0, nPiZero = 0, nPiMinus = 0;
     687    G4int nH2 = 0, nHe3 = 0, nAlpha = 0;
     688    G4int nFragments = 0;
     689    G4int nParticles = 0;
     690    G4cout <<"Particles produced in the event (" << ntrack << "):" << G4endl;
     691    G4cout <<"A \t Z \t Ekin \t Ptot \t Theta \t Phi" << G4endl;
     692    for(G4int i = 0; i < ntrack; i++) {
     693      nParticles++;
     694      if(avv[i] ==  1 && zvv[i] ==  1) nProton++;  // Count multiplicities
     695      if(avv[i] ==  1 && zvv[i] ==  0) nNeutron++;
     696      if(avv[i] == -1 && zvv[i] ==  1) nPiPlus++;
     697      if(avv[i] == -1 && zvv[i] ==  0) nPiZero++;
     698      if(avv[i] == -1 && zvv[i] == -1) nPiMinus++;
     699      if(avv[i] ==  2 && zvv[i] ==  1) nH2++;
     700      if(avv[i] ==  3 && zvv[i] ==  2) nHe3++;
     701      if(avv[i] ==  4 && zvv[i] ==  2) nAlpha++;
     702      if(                zvv[i] >   2) nFragments++;
     703
     704      G4cout << i << " \t " << avv[i] << " \t " << zvv[i] << " \t " << enerj[i] << " \t "
     705             << plab[i] << " \t " << tetlab[i] << " \t " << philab[i] << G4endl;
     706    }
     707
     708    G4cout <<"Summary of event: " << G4endl;
     709    G4cout <<"Projectile type: " << projType <<" Energy: " << projEnergy << G4endl;
     710    G4cout <<"Target A = " << targetA << " Z = " << targetZ << G4endl;
     711    G4cout <<"Remnant from cascade: " << G4endl;
     712    G4cout <<"A = " << massini << " Z = " << mzini << " excitation E = " << exini << G4endl;
     713    G4cout <<"Particle multiplicities:" << G4endl;
     714    G4cout <<"Protons: " << nProton << " Neutrons:  " << nNeutron << G4endl;
     715    G4cout <<"pi+: " << nPiPlus << " pi0: " << nPiZero << " pi-: " << nPiMinus << G4endl;
     716    G4cout <<"H2: " << nH2 << " He3: " << nHe3 << " Alpha: " << nAlpha << G4endl;
     717    G4cout <<"Nucleus fragments = " << nFragments << G4endl;
     718    G4cout <<"Conservation laws:" << G4endl;
     719    G4cout <<"Baryon number = " <<  getTotalBaryonNumber() << G4endl;
     720    G4cout <<"Number of particles = " << nParticles << G4endl;
     721  }
     722
     723  /**
     724   * Projectile type.
     725   */
     726  G4int projType;
     727
     728  /**
     729   * Projectile energy.
     730   */
     731  G4double projEnergy;
     732
     733  /**
     734   * Target mass number.
     735   */
     736  G4int targetA;
     737
     738  /**
     739   * Target charge number.
     740   */
     741  G4int targetZ;
     742
     743  /**
    566744   * A of the remnant.
    567745   */
    … …  
    578756  G4double exini;
    579757
     758  G4double pcorem, mcorem, pxrem, pyrem, pzrem;
     759
    580760  /**
    581761   * Cascade n multip.
    … …  
    629809
    630810  /**
     811   * The state of the index:
     812   * true = reserved
     813   * false = free
     814   */
     815  G4bool full[VARNTPSIZE];
     816
     817  /**
    631818   * emitted in cascade (0) or evaporation (1).
    632819   */
    … …  
    663850   */
    664851  G4double philab[VARNTPSIZE];
     852
     853private:
     854  G4int particleIndex;
    665855};
    666856
Note: See TracChangeset for help on using the changeset viewer.