Ignore:
Timestamp:
Dec 22, 2010, 3:52:27 PM (13 years ago)
Author:
garnier
Message:

geant4 tag 9.4

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/models/incl/src/G4Incl.cc

    r1340 r1347  
    2424// ********************************************************************
    2525//
    26 // $Id: G4Incl.cc,v 1.33 2010/10/29 06:48:43 gunter Exp $
     26// $Id: G4Incl.cc,v 1.37 2010/12/15 07:41:31 gunter Exp $
    2727// Translation of INCL4.2/ABLA V3
    2828// Pekka Kaitaniemi, HIP (translation)
    … …  
    5757  randomGenerator = new G4InclGeant4Random();
    5858  //  randomGenerator = new G4Ranecu();
     59  clearState();
    5960}
    6061
    … …  
    8586  //  randomGenerator = new G4Ranecu();
    8687  randomGenerator = new G4InclGeant4Random();
     88  clearState();
    8789}
    8890
    … …  
    149151
    150152  theLogger = 0;
     153  clearState();
    151154}
    152155
    153156G4Incl::~G4Incl()
    154157{
     158  delete be;
     159  delete ps;
     160  delete fermi;
     161  delete qvp;
     162
    155163  delete randomGenerator;
    156164  delete light_gaus_nuc;
    … …  
    385393void G4Incl::processEventIncl(G4InclInput *input)
    386394{
     395  //G4cout <<"Starting event " << eventnumber << G4endl;
    387396  if(input == 0) {
    388397    G4cerr <<"G4Incl fatal error: NULL pointer passed as input!" << G4endl;
    … …  
    393402  const G4double uma = 931.4942;
    394403  const G4double melec = 0.511;
    395 
    396   G4double   pcorem = 0.0;
    397   G4double   pxrem = 0.0;
    398   G4double   pyrem = 0.0;
    399   G4double   pzrem = 0.0;
     404  const G4double fmp = 938.2796;
     405
     406  G4double pcorem = 0.0;
     407  G4double pxrem = 0.0;
     408  G4double pyrem = 0.0;
     409  G4double pzrem = 0.0;
    400410
    401411  G4double ap = 0.0, zp = 0.0, mprojo = 0.0, pbeam = 0.0;
    402412
    403413  varntp->clear();
    404   if(calincl->bulletType() == 3) { // pi+
     414
     415  if(calincl->bulletType() == -12) {
     416    be->ia_be = std::abs(calincl->bulletType());
     417    be->iz_be = 6;
     418  } else if(calincl->bulletType() == -666) {
     419    be->iz_be = calincl->extendedProjectileZ();
     420    be->ia_be = calincl->extendedProjectileA();
     421  }
     422
     423  if(calincl->isExtendedProjectile() == false && calincl->bulletType() < -max_a_proj) {
     424  //  if(calincl->bulletType() < -max_a_proj) {
     425    G4cout <<"max a of composite projectile is: " << max_a_proj << G4endl;
     426    exit(0);
     427  }
     428  if(calincl->bulletType() < 0) {
     429    //    calincl->bulletType() = std::floor(calincl->bulletType() + 0.1); WTF???
     430    be->pms_be=100.;
     431    G4int i_tabled=0;
     432    if(be->iz_be == 3 && be->ia_be == 6) {
     433      be->rms_be=2.56;
     434      be->bind_be=32.0;
     435      i_tabled=1;
     436    } else if(be->iz_be == 3 && be->ia_be == 7) { // TODO: Check the values!
     437      be->rms_be=2.56;
     438      be->bind_be=32.0;
     439      i_tabled=1;
     440    } else if(be->iz_be == 3 && be->ia_be == 8) {
     441      be->rms_be=2.40;
     442      be->bind_be=39.25;
     443      i_tabled=1;
     444    } else if(be->iz_be == 4 && be->ia_be == 7) {
     445      be->rms_be=2.51;
     446      be->bind_be=58.17;
     447      i_tabled=1;
     448    } else if(be->iz_be == 4 && be->ia_be == 9) {
     449      be->rms_be=2.51;
     450      be->bind_be=58.17;
     451      i_tabled=1;
     452    } else if(be->iz_be == 4 && be->ia_be == 10) {
     453      be->rms_be=2.45;
     454      be->bind_be=64.75;
     455      i_tabled=1;
     456    } else if(be->iz_be == 5 && be->ia_be == 10) {
     457      be->rms_be=2.45;
     458      be->bind_be=64.75;
     459      i_tabled=1;
     460    } else if(be->iz_be == 5 && be->ia_be == 11) {
     461      be->rms_be=2.40;
     462      be->bind_be=76.21;
     463      i_tabled=1;
     464    } else if(be->iz_be == 6 && be->ia_be == 9) { // TODO: Check the values!
     465      be->rms_be=2.44;
     466      be->bind_be=92.17;
     467      i_tabled=1;
     468    } else if(be->iz_be == 6 && be->ia_be == 10) { // TODO: Check the values!
     469      be->rms_be=2.44;
     470      be->bind_be=92.17;
     471      i_tabled=1;
     472    } else if(be->iz_be == 6 && be->ia_be == 11) { // TODO: Check the values!
     473      be->rms_be=2.44;
     474      be->bind_be=92.17;
     475      i_tabled=1;
     476    } else if(be->iz_be == 6 && calincl->bulletType() == -12) { // Special Carbon case
     477      G4cout <<"Carbon 12 (special) selected." << G4endl;
     478      be->rms_be=2.44;
     479      be->bind_be=92.17;
     480      i_tabled=1;
     481    } else if(be->iz_be == 6 && be->ia_be == 12) {
     482      be->rms_be=2.44;
     483      be->bind_be=92.17;
     484      i_tabled=1;
     485    } else if(be->iz_be == 7 && be->ia_be == 16) {
     486      be->rms_be=2.73;
     487      be->bind_be=127.62;
     488      i_tabled=1;
     489    } else {
     490      G4cout <<"Warning: No rms and binding for projectile ion A = " << be->ia_be << " Z = " << be->iz_be << G4endl;
     491      be->rms_be=2.44;
     492      be->bind_be=92.17;
     493      G4cout <<"Warning: Using probably bad values rms = " << be->rms_be << " binding = " << be->bind_be << G4endl;
     494      i_tabled=1;     
     495    }
     496     
     497    if(i_tabled == 0) {
     498      G4cout <<"This heavy ion (a,z)= " << be->ia_be << " " << be->iz_be << " is not defined as beam in INCL" << G4endl;
     499      exit(0);
     500    }
     501     
     502    //    G4cout <<"z projectile, rms_r, rms_p (gaussian model)" << be->iz_be << " " << be->rms_be << " " << be->pms_be << G4endl;
     503    //    G4cout <<"binding energy (mev):" << be->bind_be << G4endl;
     504    //    G4cout <<"fermi-breakup dresner below a=" << calincl->f[11] << G4endl;
     505  }     
     506  //  G4cout <<"Target Mass and Charge: " << calincl->targetA() << " " << calincl->targetZ() << G4endl;
     507  //  calincl->f[10] = 0; // No clusters
     508
     509  if(calincl->bulletType() == -12) {  // C12 special case
     510    mprojo=fmp*std::abs(calincl->bulletType()) - be->bind_be;
     511    pbeam=std::sqrt(calincl->bulletE()*(calincl->bulletE()+2.*mprojo));
     512    ap=std::abs(calincl->bulletType());
     513    zp=be->iz_be;
     514  } else if(calincl->bulletType() == -666) { // Generic extended projectile
     515    mprojo=fmp*be->ia_be - be->bind_be;
     516    pbeam=std::sqrt(calincl->bulletE()*(calincl->bulletE()+2.*mprojo));
     517    ap=be->ia_be;
     518    zp=be->iz_be;
     519  }
     520  // pi+
     521  if(calincl->bulletType() == 3) {
    405522    mprojo = 139.56995;
    406523    ap = 0.0;
    … …  
    408525  }
    409526
    410   if(calincl->bulletType() == 4) { // pi0
     527  // pi0
     528  if(calincl->bulletType() == 4) {
    411529    mprojo = 134.9764;
    412530    ap = 0.0;
    … …  
    414532  }
    415533
    416   if(calincl->bulletType() == 5) { // pi-
     534  // pi-
     535  if(calincl->bulletType() == 5) {
    417536    mprojo = 139.56995;
    418537    ap = 0.0;
    … …  
    420539  }
    421540
    422   // Coulomb en entree seulement pour les particules ci-dessous.
    423   if(calincl->bulletType() == 1) { // proton
     541  // coulomb en entree seulement pour les particules ci-dessous
     542
     543  // proton
     544  if(calincl->bulletType() == 1) {
    424545    mprojo = 938.27231;
    425546    ap = 1.0;
    426547    zp = 1.0;
    427548  }
    428  
    429   if(calincl->bulletType() == 2) { // neutron 
     549
     550  // neutron 
     551  if(calincl->bulletType() == 2) {
    430552    mprojo = 939.56563;
    431553    ap = 1.0;
    432554    zp = 0.0;
    433555  }
    434  
    435   if(calincl->bulletType() == 6) { // deuteron
     556
     557  // deuteron
     558  if(calincl->bulletType() == 6) {
    436559    mprojo = 1875.61276;
    437560    ap = 2.0;
    438561    zp = 1.0;
    439562  }
    440  
    441   if(calincl->bulletType() == 7) { // triton
     563
     564  // triton
     565  if(calincl->bulletType() == 7) {
    442566    mprojo = 2808.95;
    443567    ap = 3.0;
    444568    zp = 1.0;
    445569  }
    446  
    447   if(calincl->bulletType() == 8) { // He3
     570
     571  // He3
     572  if(calincl->bulletType() == 8) {
    448573    mprojo = 2808.42;
    449574    ap = 3.0;
    … …  
    451576  }
    452577
    453   if(calincl->bulletType() == 9) { // Alpha
     578  // Alpha
     579  if(calincl->bulletType() == 9) {
    454580    mprojo = 3727.42;
    455581    ap = 4.0;
    … …  
    457583  }
    458584
    459   if(calincl->bulletType() == -12) { // Carbon
    460     mprojo = 12.0 * 938.27231;
     585  // Carbon
     586  if(calincl->bulletType() == -12) {
     587    mprojo = 6.0*938.27231 + 6.0*939.56563;
    461588    ap = 12.0;
    462589    zp = 6.0;
    463   } else if(calincl->bulletType() == -666) { // Generic extended ion projectile
    464     mprojo = calincl->extendedProjectileA() * 938.27231;
    465     ap = calincl->extendedProjectileA();
    466     zp = calincl->extendedProjectileZ();
    467590  }
    468591
    … …  
    482605  G4double bimpac = 0.0;
    483606  G4int jrem = 0;
     607  G4double xjrem = 0.0, yjrem = 0.0, zjrem = 0.0;
    484608  G4double alrem = 0.0;
    485609
    486   /**
    487    * Coulomb barrier treatment.
    488    */
     610  // Coulomb barrier
     611 
    489612  G4double probaTrans = 0.0;
    490613  G4double rndm = 0.0;
    491614
    492615  if((calincl->bulletType() == 1) || (calincl->bulletType() >= 6)) {
     616    //    probaTrans = coulombTransm(calincl->bulletE(),apro,zpro,calincl->targetA(),calincl->targetZ());
    493617    probaTrans = coulombTransm(calincl->bulletE(),ap,zp,calincl->targetA(),calincl->targetZ());
    494618    standardRandom(&rndm, &(hazard->ial));
    … …  
    499623  }
    500624
    501   /**
    502    * Call the actual INCL routine.
    503    */
    504 
     625  //  G4cout <<"Before PNU:" << G4endl;
    505626  //  randomGenerator->printSeeds();
    506   G4double jremx, jremy, jremz;
    507   pnu(&ibert, &nopart,&izrem,&iarem,&esrem,&erecrem,&alrem,&berem,&garem,&bimpac,&jrem,
    508       &jremx, &jremy, &jremz);
     627  // Call the actual INCL routine:
     628  pnu(&ibert, &nopart,&izrem,&iarem,&esrem,&erecrem,&alrem,&berem,&garem,&bimpac,
     629      &jrem, &xjrem, &yjrem, &zjrem);
     630  //  G4cout <<"After PNU:" << G4endl;
     631  //  randomGenerator->printSeeds();
     632  G4double mrem = int(zjrem/197.328); // CHECK
     633  if (mrem > jrem) mrem=jrem;
     634  if (mrem < -jrem) mrem=-jrem;
     635
     636//   nopart=1;
     637//   kind[0]=1;
     638//   ep[0]=799.835;
     639//   alpha[0]=0.08716;
     640//   beta[0]=0.;
     641//   gam[0]=0.99619;
     642//   izrem=82;
     643//   iarem=208;
     644//   esrem=200.;
     645//   erecrem=0.18870;
     646//   alrem=-0.47101;
     647//   berem=0.;
     648//   garem=0.88213;
     649//   bimpac=2.;
    509650  forceAbsor(&nopart, &iarem, &izrem, &esrem, &erecrem, &alrem, &berem, &garem, &jrem);
    510   G4double aprf = double(iarem); // mass number of the prefragment
    511   G4double jprf = 0.0;           // angular momentum of the prefragment
    512 
    513   // Mean angular momentum of prefragment.                                 
    514   jprf = 0.165 * std::pow(at,(2.0/3.0)) * aprf * (at - aprf) / (at - 1.0);                               
     651  G4double aprf = double(iarem);    // mass number of the prefragment
     652  G4double jprf = 0.0;                // angular momentum of the prefragment
     653
     654  // Mean angular momentum of prefragment                                 
     655  jprf = 0.165 * std::pow(at,(2.0/3.0)) * aprf*(at - aprf)/(at - 1.0);                               
    515656  if (jprf < 0) {
    516657    jprf = 0.0;
    517658  }
    518659
    519   // Reference M. de Jong, Ignatyuk, Schmidt Nuc. Phys A 613, p442, 7th line
     660  // check m.de jong, ignatyuk, schmidt nuc.phys a 613, pg442, 7th line
    520661  jprf = std::sqrt(2*jprf);
     662
    521663  jprf = jrem;
    522   varntp->jremn = jrem; // Copy jrem to output tuple.
    523 
    524   G4double numpi = 0;  // Number  of pions.
    525   G4double multn = 0;  // Number (multiplicity) of neutrons.
    526   G4double multp = 0;  // Number (multiplicity) of protons.
    527 
    528   // Ecriture dans le ntuple des particules de cascade (sauf remnant).     
    529   varntp->ntrack = nopart; // Nombre de particules pour ce tir.
     664  varntp->jremn = jrem;      // jrem copie dans le ntuple
     665
     666  G4double numpi = 0;  // compteurs de pions, neutrons protons
     667  G4double multn = 0;
     668  G4double multp = 0;
     669
     670  // ecriture dans le ntuple des particules de cascade (sauf remnant)     
     671  varntp->ntrack = nopart;          // nombre de particules pour ce tir
     672  if(varntp->ntrack >= VARNTPSIZE) {
     673    if(verboseLevel > 2) {
     674      G4cout <<"G4Incl error: Output data structure not big enough." << G4endl;
     675    }
     676  }
    530677  varntp->massini = iarem;
    531678  varntp->mzini = izrem;
    … …  
    533680  varntp->bimpact = bimpac;
    534681 
    535   /**
    536    * Three ways to compute the mass of the remnant:
    537    * -from the output of the cascade and the canonic mass
    538    * -from energy balance (input - all emitted energies)
    539    * -following the approximations of the cugnon code (esrem...)
    540    */
    541   G4double f0 = calincl->targetA();
    542   G4double f1 = calincl->targetZ();
    543   G4double f2 = calincl->bulletE();
    544   G4double mcorem = mprojo + f2 + abla->pace2(f0, f1) + f0 * uma - f1 * melec;
     682  //  three ways to compute the mass of the remnant:
     683  //                -from the output of the cascade and the canonic mass
     684  //                -from energy balance (input - all emitted energies)
     685  //                -following the approximations of the cugnon code (esrem...)
     686  G4double mcorem = mprojo + calincl->bulletE() + abla->pace2(double(calincl->targetA()),double(calincl->targetZ()))
     687    + calincl->targetA()*uma - calincl->targetZ()*melec;
    545688
    546689  G4double pxbil = 0.0;
    … …  
    549692
    550693  if(nopart > -1) {
    551     for(G4int j = 0; j < nopart; j++) {
    552       varntp->itypcasc[j] = 1;
     694    // Fill the projectile spectator variables
     695    varntp->masp = ps->a_projspec;
     696    varntp->mzsp = ps->z_projspec;
     697    varntp->exsp = ps->ex_projspec;
     698    varntp->spectatorP1 = ps->p1_projspec;
     699    varntp->spectatorP2 = ps->p2_projspec;
     700    varntp->spectatorP3 = ps->p3_projspec;
     701    varntp->spectatorT = ps->t_projspec;
     702    for(G4int j = 0; j <= nopart; j++) {
     703      if(ep[j] < 0.0) { // Workaround to avoid negative energies (and taking std::sqrt of a negative number).
     704        G4cout <<"G4Incl: Not registering particle with energy: " << ep[j] << G4endl;
     705        continue;
     706      }
     707      if(kind[j] == 0) continue; // Empty particle rows are sometimes produced by lurking indexing problems. We can simply skip those "bad" entries...
     708      if(gam[j] > CLHEP::pi) {
     709        if(verboseLevel > 2) {
     710          G4cout <<"G4Incl: Just avoided floating point exception by using an ugly hack..." << G4endl;
     711        }
     712        continue; // Avoid floating point exception
     713      }
     714
     715      varntp->itypcasc[j] = 1; // Particle was produced by the cascade
     716      // Spectators of composite projectiles (7/2006, AB)
     717      // (kind is negative in that case)
     718      if(kind[j] <= 0) { // Particle is a projectile spectator that comes directly from the cascade
     719        kind[j] *= -1;
     720        varntp->itypcasc[j]=-1;
     721        //      G4cout <<"Spectator registered!" << G4endl;
     722        //      continue;
     723      }
     724       
    553725      // kind(): 1=proton, 2=neutron, 3=pi+, 4=pi0, 5=pi -     
    554726      if(kind[j] == 1) {
    555727        varntp->avv[j] = 1;
    556728        varntp->zvv[j] = 1;
    557         varntp->plab[j] = std::sqrt(ep[j] * (ep[j] + 1876.5592)); // cugnon
     729        varntp->plab[j] = std::sqrt(ep[j]*(ep[j]+1876.5592)); // cugnon
    558730        multp = multp + 1;
    559731        mcorem = mcorem - ep[j] - 938.27231;
    … …  
    568740        varntp->zvv[j] = 0;
    569741        varntp->plab[j] = std::sqrt(ep[j]*(ep[j]+1876.5592)); // cugnon
    570         //varntp->plab[j] = std::sqrt(ep[j] * (ep[j] + 1879.13126)); // PK mass check
    571742        multn = multn + 1;
    572743        mcorem = mcorem - ep[j] - 939.56563;
    … …  
    576747        }
    577748      }
    578 
     749     
    579750      if(kind[j] == 3) {
    580751        varntp->avv[j] = -1;
    … …  
    656827        mcorem = mcorem - ep[j] - 3724.818;
    657828        if(verboseLevel > 3) {
    658           G4cout <<"G4Incl: He4 produced! " << G4endl;
     829          G4cout <<"G4Incl: He3 produced! " << G4endl;
    659830          G4cout <<"G4Incl: Momentum: "<< varntp->plab[j] << G4endl;
    660831        }
    … …  
    669840
    670841      if(verboseLevel > 3) {
    671         G4cout <<"Momentum: " << varntp->plab[j] << G4endl;
    672         G4cout <<"Theta: " << varntp->tetlab[j] << G4endl;
    673         G4cout <<"Phi: " << varntp->philab[j] << G4endl;
     842        G4cout <<"Momentum: " << varntp->plab[j]   << G4endl;
     843        G4cout <<"Theta: "    << varntp->tetlab[j] << G4endl;
     844        G4cout <<"Phi: "      << varntp->philab[j] << G4endl;
    674845      }
    675846    }
    676847
    677848    // calcul de la masse (impulsion) du remnant coherente avec la conservation d'energie:
    678     pcorem=std::sqrt(erecrem*(erecrem +2.*938.2796*iarem));   // cugnon
    679     mcorem = 938.2796*iarem;                               // cugnon
    680                
     849    pcorem = std::sqrt(erecrem*(erecrem + 2.0 * 938.2796 * iarem));   // cugnon
     850    mcorem = 938.2796 * iarem;                               // cugnon
     851    varntp->pcorem = pcorem;
     852    varntp->mcorem = mcorem;
    681853    // Note: Il faut negliger l'energie d'excitation (ESREM) pour que le bilan
    682854    // d'impulsion soit correct a la sortie de la cascade.....et prendre la
    683855    // masse MCOREM comme ci-dessus (fausse de ~1GeV par rapport aux tables...)       
    684     pxrem=pcorem*alrem;
    685     pyrem=pcorem*berem;
    686     pzrem=pcorem*garem;
    687        
    688     pxbil=pxbil+pxrem;
    689     pybil=pybil+pyrem;
    690     pzbil=pzbil+pzrem;
    691 
    692     if((std::fabs(pzbil-pbeam) > 5.0) || (std::sqrt(std::pow(pxbil,2)+std::pow(pybil,2)) >= 3.0)) {
     856    pxrem = pcorem * alrem;
     857    pyrem = pcorem * berem;
     858    pzrem = pcorem * garem;
     859    varntp->pxrem = pxrem;
     860    varntp->pyrem = pyrem;
     861    varntp->pzrem = pzrem;
     862    pxbil = pxbil + pxrem;
     863    pybil = pybil + pyrem;
     864    pzbil = pzbil + pzrem;
     865
     866    // If on purpose, add here the spectator nuclei:   
     867    if(calincl->bulletType() < 0 && ps->a_projspec != 0) {
     868      pxbil=pxbil+ps->p1_projspec;
     869      pybil=pybil+ps->p2_projspec;
     870      pzbil=pzbil+ps->p3_projspec;
     871    }
     872
     873    if((std::fabs(pzbil - pbeam) > 5.0) || (std::sqrt(std::pow(pxbil,2) + std::pow(pybil,2)) >= 3.0)) {
    693874      if(verboseLevel > 3) {
    694875        G4cout <<"Bad momentum conservation after INCL:" << G4endl;
    … …  
    704885    varntp->iafis = 0;
    705886
    706     //  varntp->ntrack = varntp->ntrack + 1;  // on recopie le remnant dans le ntuple
     887    // on recopie le remnant dans le ntuple
     888    varntp->ntrack = varntp->ntrack + 1;
    707889    varntp->massini = iarem;
    708890    varntp->mzini = izrem;
    709891    varntp->exini = esrem;
    710     varntp->itypcasc[varntp->ntrack] = 1;
    711     varntp->avv[varntp->ntrack] = iarem;
    712     varntp->zvv[varntp->ntrack]= izrem;
    713     varntp->plab[varntp->ntrack] = pcorem;
    714     varntp->enerj[varntp->ntrack] = std::sqrt(std::pow(pcorem,2) + std::pow(mcorem,2)) - mcorem;
    715     varntp->tetlab[varntp->ntrack] = 180.0*std::acos(garem)/3.141592654;
    716     varntp->philab[varntp->ntrack] = 180.0*std::atan2(berem,alrem)/3.141592654;
    717     varntp->ntrack++;
    718     varntp->mulncasc = varntp->ntrack;
    719     varntp->mulnevap = 0;
    720     varntp->mulntot = varntp->mulncasc + varntp->mulnevap;
    721     if(verboseLevel > 3) {
    722       G4cout <<"G4Incl: Returning nucleus fragment. " << G4endl;
    723       G4cout <<"G4Incl: Fragment A = " << varntp->avv[varntp->ntrack] << " Z = " << varntp->zvv[varntp->ntrack] << G4endl;
    724       G4cout <<"Energy: " << varntp->enerj[varntp->ntrack] << G4endl;
    725       G4cout <<"Momentum: " << varntp->plab[varntp->ntrack] << G4endl;
    726       G4cout <<"Theta: " << varntp->tetlab[varntp->ntrack] << G4endl;
    727       G4cout <<"Phi: " << varntp->philab[varntp->ntrack] << G4endl;
    728     }
    729   }
    730   else {
    731     if(nopart == -2) {
    732       varntp->ntrack = -2; //FIX: Error flag to remove events containing unphysical events (Ekin > Ebullet).
    733     }
    734     else {
    735       varntp->ntrack = -1;
    736     }
     892    varntp->pxrem = pxrem;
     893    varntp->pyrem = pyrem;
     894    varntp->pzrem = pzrem;
     895    varntp->mcorem = mcorem;
     896    varntp->erecrem = pcorem;
     897    varntp->erecrem = erecrem;
     898
     899#ifdef G4INCLDEBUG
     900    theLogger->fillHistogram1D("bimpact", varntp->bimpact);
     901#endif
     902
     903#ifdef G4INCLDEBUG
     904    theLogger->fillHistogram1D("mzini", varntp->mzini);
     905#endif
     906  }
     907  if(nopart == -2) {
     908    varntp->ntrack = -2; //FIX: Error flag to remove events containing unphysical events (Ekin > Ebullet).
     909    evaporationResult->ntrack = -2; //FIX: Error flag to remove events containing unphysical events (Ekin > Ebullet).
     910  }
     911  else if(nopart == -1) {
     912    varntp->ntrack = -1;
     913    evaporationResult->ntrack = -1;
     914  }
     915  if(verboseLevel > 2) {
     916    G4cout << __FILE__ << ":" << __LINE__ << "Dump varntp after combining: " << G4endl;
     917    varntp->dump();
    737918  }
    738919}
    … …  
    81508331}
    81518332
    8152 void G4Incl::standardRandom(G4double *rndm, G4long *seed)
    8153 {
    8154   (*seed) = (*seed); // Avoid warning during compilation.
     8333void G4Incl::standardRandom(G4double *rndm, G4long*)
     8334{
    81558335  // Use Geant4 G4UniformRand
    81568336  //  (*rndm) = G4UniformRand();
    … …  
    85508730    } // enddo
    85518731    G4double p_spec2=std::pow(p1_spec,2)+std::pow(p2_spec,2)+std::pow(p3_spec,2);
    8552     G4double s_spec = sqrt(std::pow(e_spec,2)-p_spec2);
     8732    G4double s_spec = std::sqrt(std::pow(e_spec,2)-p_spec2);
    85538733
    85548734    // no projectile spectator if a>=4 and a=z or a=n (no dresner breakup)
    … …  
    88449024}
    88459025
     9026// Initialization helper
     9027void G4Incl::clearState() {
     9028  G4int epSize = 300;
     9029  for(G4int i = 0; i < epSize; ++i) {
     9030    kind[i] = 0;
     9031    ep[i] = 0.0;
     9032    alpha[i] = 0.0;
     9033    beta[i] = 0.0;
     9034    gam[i] = 0.0;
     9035  }
     9036  inside_step = 0;
     9037}
     9038
    88469039// C------------------------------------------------------------------------
    88479040// C Logging
    … …  
    89239116//       COMMON/BL2/CROIS(19900),K,IND(20000),JND(20000)
    89249117//       COMMON/BL3/R1,R2,X1(300),X2(300),X3(300),IA1,IA2,RAB2             P-N22270
     9118
     9119  if(index < 0) {
     9120    G4cout <<";; Error! index < 0! index = " << index << G4endl;
     9121    return;
     9122  }
    89259123
    89269124  G4int i_ind1 = bl1->ind1[bl2->ind[index]];
Note: See TracChangeset for help on using the changeset viewer.