Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DPMJET2_5Model.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * *
21 // * Parts of this code which have been developed by QinetiQ Ltd *
22 // * under contract to the European Space Agency (ESA) are the *
23 // * intellectual property of ESA. Rights to use, copy, modify and *
24 // * redistribute this software for general public use are granted *
25 // * in compliance with any licensing, distribution and development *
26 // * policy adopted by the Geant4 Collaboration. This code has been *
27 // * written by QinetiQ Ltd for the European Space Agency, under ESA *
28 // * contract 19770/06/NL/JD (Technology Research Programme). *
29 // * *
30 // * By using, copying, modifying or distributing the software (or *
31 // * any work based on the software) you agree to acknowledge its *
32 // * use in resulting scientific publications, and indicate your *
33 // * acceptance of all terms of the Geant4 Software license. *
34 // ********************************************************************
35 //
38 //
39 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
40 //
41 // MODULE: G4DPMJET2_5Model.cc
42 //
43 // Version: 0.B
44 // Date: 02/04/08
45 // Author: P R Truscott
46 // Organisation: QinetiQ Ltd, UK
47 // Customer: ESA/ESTEC, NOORDWIJK
48 // Contract: 19770/06/NL/JD
49 //
50 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
52 //
53 #ifdef G4_USE_DPMJET
54 
55 
56 #include "G4DPMJET2_5Model.hh"
57 #include "G4GlaubAADataSetHandler.hh"
58 
59 #include "G4ExcitationHandler.hh"
60 #include "G4Evaporation.hh"
61 #include "G4FermiBreakUp.hh"
62 #include "G4PhotonEvaporation.hh"
63 #include "G4PreCompoundModel.hh"
64 #include "G4ParticleDefinition.hh"
65 #include "G4ParticleTable.hh"
66 #include "G4DynamicParticle.hh"
67 #include "Randomize.hh"
68 #include "G4Fragment.hh"
69 #include "G4VNuclearDensity.hh"
70 #include "G4NuclearShellModelDensity.hh"
71 #include "G4NuclearFermiDensity.hh"
72 #include "G4FermiMomentum.hh"
73 #include "G4ReactionProductVector.hh"
74 #include "G4LorentzVector.hh"
75 #include "G4ParticleMomentum.hh"
76 #include "G4Poisson.hh"
77 #include "G4ParticleTable.hh"
78 #include "G4IonTable.hh"
79 #include "G4LorentzVector.hh"
80 #include "G4HadTmpUtil.hh"
81 #include "globals.hh"
82 
83 #include <fstream>
84 
85 #include "G4DPMJET2_5Interface.hh"
86 
88 //
89 // Constructor without arguments
90 //
91 // This constructor uses a default pre-compound. It initialises the
92 // variables (including the de-excitation), but note that much of the work is
93 // done in the member function Initialise(), which is dedicated to
94 // initialising variables in DPMJET-II.5 to class-default values.
95 //
96 G4DPMJET2_5Model::G4DPMJET2_5Model () : G4HadronicInteraction("DPMJET2_5")
97 {
98 //
99 // Set the default verbose level to 0 - no output.
100 //
101  SetVerboseLevel(1);
102 //
103 #ifdef G4VERBOSE
104  if (GetVerboseLevel()>0) {
105  G4cout <<"G4DPMJET2_5Model default constructor" <<G4endl;
106  }
107 #endif
108 //
109 //
110 // Send message to stdout to advise that the G4DPMJET2_5Model model is being used.
111 //
112  theInitType = DEFAULT;
113  PrintWelcomeMessage();
114 //
115 // Use the precompound model for nuclear de-excitation.
116 //
117  theExcitationHandler = 0;
118  SetDefaultPreCompoundModel();
119 //
120 //
121 // Set the minimum and maximum range for the model (despite nomanclature, this
122 // is in energy per nucleon number).
123 //
124  SetMinEnergy(5.0*GeV);
125  SetMaxEnergy(1000.0*TeV);
126 //
127 //
128 // Initialise the DPMJET model - this effectively does what the DPMJET
129 // subroutine DMINIT does without reading an input file.
130 //
131  debug = false;
132  debug_level = 0;
133  lunber = 14;
134  dpmver = 2.5;
135 
136  LFALSE = 0;
137  LTRUE = 1;
138 
139  Initialise ();
140 //
141 //
142 // Next bit directs how and how many Glauber data sets are loaded
143 // or created.
144 //
145  theGlauberDataSetHandler = G4GlaubAADataSetHandler::getInstance();
146  theGlauberDataSetHandler->SetMaxGlauberDataSets(-1);
147 
148  theParticleTable = G4ParticleTable::GetParticleTable();
149  theIonTable = const_cast <G4IonTable *> (theParticleTable->GetIonTable());
150 //
151 //
152 }
154 //
155 // Constructor with DPMJET-II.5 initialisation type
156 //
157 // This constructor uses a default pre-compound. It initialises the
158 // variables (including the de-excitation), but note that much of the work is
159 // done in the member function Initialise(), which is dedicated to
160 // initialising variables in DPMJET-II.5. The user is able to define whether to
161 // use the default values, DPMJET-II.5 settings or DPMJET-III settings.
162 //
163 G4DPMJET2_5Model::G4DPMJET2_5Model (const G4DPMJET2_5InitialisationType initType)
164 {
165 //
166 // Set the default verbose level to 0 - no output.
167 //
168  SetVerboseLevel(1);
169 //
170 #ifdef G4VERBOSE
171  if (GetVerboseLevel()>0) {
172  G4cout <<"G4DPMJET2_5Model constructor #1 " <<G4endl;
173  }
174 #endif
175 //
176 //
177 // Send message to stdout to advise that the G4DPMJET2_5Model model is being used.
178 //
179  theInitType = initType;
180  PrintWelcomeMessage();
181 //
182 // Use the precompound model for nuclear de-excitation.
183 //
184  theExcitationHandler = 0;
185  SetDefaultPreCompoundModel();
186 //
187 //
188 // Set the minimum and maximum range for the model (despite nomanclature, this
189 // is in energy per nucleon number).
190 //
191  SetMinEnergy(5.0*GeV);
192  SetMaxEnergy(1000.0*TeV);
193 //
194 //
195 // Initialise the DPMJET model - this effectively does what the DPMJET
196 // subroutine DMINIT does without reading an input file.
197 //
198  debug = false;
199  debug_level = 0;
200  lunber = 14;
201  dpmver = 2.5;
202 
203  LFALSE = 0;
204  LTRUE = 1;
205 
206  Initialise ();
207 //
208 //
209 // Next bit directs how and how many Glauber data sets are loaded
210 // or created.
211 //
212  theGlauberDataSetHandler = G4GlaubAADataSetHandler::getInstance();
213  theGlauberDataSetHandler->SetMaxGlauberDataSets (-1);
214 
215  theParticleTable = G4ParticleTable::GetParticleTable();
216  theIonTable = theParticleTable->GetIonTable();
217 //
218 //
219 }
221 //
222 // Constructor with de-excitation handler.
223 //
224 // This constructor uses the user-provided de-excitation handler. However, it
225 // is possible for the use to provide a NULL pointer, in which case, it is
226 // assumed that the user doesn't want to simulate de-excitation - USER, BEWARE!
227 //
228 // The member function initialises the variables (including the de-excitation),
229 // but note that much of the work is done in the member function Initialise(),
230 // which is dedicated to initialising variables in DPMJET-II.5.
231 //
232 G4DPMJET2_5Model::G4DPMJET2_5Model (G4ExcitationHandler *aExcitationHandler,
233  const G4DPMJET2_5InitialisationType initType)
234 {
235 //
236 // Set the default verbose level to 0 - no output.
237 //
238  SetVerboseLevel(1);
239 //
240 // Send message to stdout to advise that the G4DPMJET2_5Model model is being used.
241 //
242 #ifdef G4VERBOSE
243  if (GetVerboseLevel()>0) {
244  G4cout <<"G4DPMJET2_5Model constructor #2 " <<G4endl;
245  }
246 #endif
247  theInitType = initType;
248  PrintWelcomeMessage();
249 //
250 // The user is able to provide the excitation handler.
251 //
252  theExcitationHandler = aExcitationHandler;
253  thePreComp = 0;
254 //
255 //
256 // Set the minimum and maximum range for the model (despite nomanclature, this
257 // is in energy per nucleon number).
258 //
259  SetMinEnergy(5.0*GeV);
260  SetMaxEnergy(1000.0*TeV);
261 //
262 //
263 // Initialise the DPMJET model - this effectively does what the DPMJET
264 // subroutine DMINIT does without reading an input file.
265 //
266  debug = false;
267  debug_level = 0;
268  lunber = 14;
269  dpmver = 2.5;
270 
271  LFALSE = 0;
272  LTRUE = 1;
273 
274  Initialise ();
275 //
276 //
277 // Next bit directs how and how many Glauber data sets are loaded
278 // or created.
279 //
280  theGlauberDataSetHandler = G4GlaubAADataSetHandler::getInstance();
281  theGlauberDataSetHandler->SetMaxGlauberDataSets (-1);
282 
283  theParticleTable = G4ParticleTable::GetParticleTable();
284  theIonTable = theParticleTable->GetIonTable();
285 //
286 //
287 }
289 //
290 // Constructor with pre-compound model.
291 //
292 // This constructor uses the user-provided pre-equilibrium model. However, it
293 // is possible for the use to provide a NULL pointer, in which case, it is
294 // assumed that the user doesn't want to simulate pre-equilibrium. - USER, BEWARE!
295 //
296 // The member function initialises the variables (including the de-excitation),
297 // but note that much of the work is done in the member function Initialise(),
298 // which is dedicated to initialising variables in DPMJET-II.5.
299 //
300 G4DPMJET2_5Model::G4DPMJET2_5Model (G4VPreCompoundModel *aPreComp,
301  const G4DPMJET2_5InitialisationType initType)
302 {
303 //
304 // Set the default verbose level to 0 - no output.
305 //
306  SetVerboseLevel(1);
307 //
308 // Send message to stdout to advise that the G4DPMJET2_5Model model is being used.
309 //
310 #ifdef G4VERBOSE
311  if (GetVerboseLevel()>0) {
312  G4cout <<"G4DPMJET2_5Model constructor #3 " <<G4endl;
313  }
314 #endif
315  theInitType = initType;
316  PrintWelcomeMessage();
317 //
318 // The user is able to provide the pre-compound model.
319 //
320  theExcitationHandler = 0;
321  thePreComp = aPreComp;
322 //
323 //
324 // Set the minimum and maximum range for the model (despite nomanclature, this
325 // is in energy per nucleon number).
326 //
327  SetMinEnergy(5.0*GeV);
328  SetMaxEnergy(1000.0*TeV);
329 //
330 //
331 // Initialise the DPMJET model - this effectively does what the DPMJET
332 // subroutine DMINIT does without reading an input file.
333 //
334  debug = false;
335  debug_level = 0;
336  lunber = 14;
337  dpmver = 2.5;
338 
339  LFALSE = 0;
340  LTRUE = 1;
341 
342  Initialise ();
343 //
344 //
345 // Next bit directs how and how many Glauber data sets are loaded
346 // or created.
347 //
348  theGlauberDataSetHandler = G4GlaubAADataSetHandler::getInstance();
349  theGlauberDataSetHandler->SetMaxGlauberDataSets (-1);
350 
351  theParticleTable = G4ParticleTable::GetParticleTable();
352  theIonTable = theParticleTable->GetIonTable();
353 //
354 //
355 }
357 //
358 // Destructor
359 //
360 G4DPMJET2_5Model::~G4DPMJET2_5Model ()
361 {
362 //
363 //
364 // The destructor doesn't have to do a great deal!
365 //
366  if (theExcitationHandler) delete theExcitationHandler;
367  if (thePreComp) delete thePreComp;
368 // delete theGlauberDataSetHandler;
369  theGlauberDataSetHandler->UnloadAllGlauberData();
370 }
372 //
373 // IsApplicable
374 //
375 // This member function simply determines whether there is relevant information
376 // in Glauber data for this projectile and target, and if the nucleus is
377 // sensible.
378 //
379 G4bool G4DPMJET2_5Model::IsApplicable (
380  const G4HadProjectile &theTrack, G4Nucleus &theTarget)
381 {
382 //
383 //
384 // Get relevant information about the projectile and target (A, Z)
385 //
386  const G4ParticleDefinition *definitionP = theTrack.GetDefinition();
387  G4int AP = definitionP->GetBaryonNumber();
388  G4int ZP = G4int(definitionP->GetPDGCharge()/eplus + 0.5);
389  G4int AT = theTarget.GetA_asInt();
390  G4int ZT = theTarget.GetZ_asInt();
391 
392  if (AP >= 2 && ZP >= 1 && AT >= 2 && ZT >=1) {
393  return theGlauberDataSetHandler->IsGlauberDataSetAvailable(AP,AT);
394  }
395  else {
396  return false;
397  }
398 }
400 //
401 // ApplyYourself
402 //
403 // Member function to process an event, and get information about the products.
404 // The phases are:
405 //
406 // (1) Determine the information about the projectile and target.
407 //
408 // (2) Identify to the Glauber data set handler which data need to be used. If
409 // the GDSH finds there are no Glauber profile data for the collision,the
410 // product is the unchanged projectile.
411 //
412 // (3) Initialise further common-block variables in DPMJET-II.5, and perform
413 // FORTRAN calls (note this is taken largely from an interpretation of CORSIKA
414 // and DPMJET-II.5 FORTRAN ... I think there's duplication in some of the
415 // initialisation on top of what Initialise() does, but CORSIKA has similar
416 // duplication. I'm not confident at removing this out at the moment.
417 //
418 // (4) Call the DPMJET-II.5 FORTRAN subroutine DPMEVT.
419 //
420 // (5) Pick out the final state particles and nuclei. In the case of nuclei
421 // use the de-excitation handler if one has been defined. Transfer all these
422 // particles to the final-state vector.
423 //
424 // (6) if very-verbose output is demanded by the user, there is a print-out
425 // of the total energy and total momentum before and after collision.
426 //
427 G4HadFinalState *G4DPMJET2_5Model::ApplyYourself (
428  const G4HadProjectile &theTrack, G4Nucleus &theTarget)
429 {
430 //
431 //
432 // The secondaries will be returned in G4HadFinalState &theParticleChange -
433 // initialise this. The original track will always be discontinued and
434 // secondaries followed.
435 //
436  theParticleChange.Clear();
437  theParticleChange.SetStatusChange(stopAndKill);
438 //
439 //
440 // Get relevant information about the projectile and target (A, Z, energy/nuc,
441 // momentum, etc).
442 //
443  const G4ParticleDefinition *definitionP = theTrack.GetDefinition();
444  G4int AP = definitionP->GetBaryonNumber();
445  G4int ZP = G4int(definitionP->GetPDGCharge()/eplus+0.5);
446  G4double M = definitionP->GetPDGMass();
447  G4ThreeVector pP = theTrack.Get4Momentum().vect();
448  G4double T = theTrack.GetKineticEnergy()/G4double(AP); // Units are MeV/nuc
449  G4double E = theTrack.GetTotalEnergy()/G4double(AP); // Units are MeV/nuc
450  G4int AT = theTarget.GetA_asInt();
451  G4int ZT = theTarget.GetZ_asInt();
452  G4double mpnt = theTarget.AtomicMass(AT, ZT);
453  G4double TotalEPre = theTrack.GetTotalEnergy() + mpnt;
454  // theTarget.AtomicMass(AT, ZT) + theTarget.GetEnergyDeposit();
455 // G4LorentzRotation transformToLab =
456 // (const_cast <G4HadProjectile*> (&theTrack))->GetTrafoToLab();
457 //
458 //
459 // Output relevant information on initial conditions if verbose. Note that
460 // most of the verbsse output is dealt with through private member function
461 // calls.
462 //
463  if (verboseLevel >= 2)
464  {
465  G4cout <<"########################################"
466  <<"########################################"
467  <<G4endl;
468  G4cout.precision(6);
469  G4cout <<"IN G4DPMJET2_5Model::ApplyYourself" <<G4endl;
470  G4cout <<"START OF EVENT" <<G4endl;
471  G4cout <<"Initial projectile (A,Z) = (" <<AP <<", " <<ZP <<")" <<G4endl;
472  G4cout <<"Initial target (A,Z) = (" <<AT <<", " <<ZT <<")" <<G4endl;
473  G4cout <<"Projectile momentum = " <<pP/MeV <<" MeV/c" <<G4endl;
474  G4cout <<"Total energy = " <<TotalEPre/MeV <<" MeV" <<G4endl;
475  G4cout <<"Kinetic energy/nuc = " <<T/MeV <<" MeV" <<G4endl;
476  }
477 //
478 //
479 // Setup variables and call the DPMJET model. There is a significant amount of
480 // initialisation which still needs to be done, based on DPMJET-II.5 main
481 // program and card reader subroutine, and the CORSIKA implementation.
482 //
483  G4int AP1 = AP;
484  G4int ZP1 = ZP;
485  G4int AT1 = AT;
486  G4int ZT1 = ZT;
487 //
488 //
489 // vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
490 // ***** WARNING *****
491 // The following is a provisional "catch" for ions with A==Z. The de-excitation
492 // and precompound model can produce such nuclei, although they should decay
493 // into free protons. For the moment, DPMJET-II.5 doesn't treat them ... i.e.
494 // the FORTRAN code would crash. Therefore, return such ions without
495 // nuclear interactions.
496 //
497  if (AP1 > 1 && AP1 == ZP1) {
498  theParticleChange.SetStatusChange(isAlive);
499  theParticleChange.SetEnergyChange(theTrack.GetKineticEnergy());
500  theParticleChange.SetMomentumChange(theTrack.Get4Momentum().vect().unit());
501  if (verboseLevel >= 2) {
502  G4cout <<"PROJECTILE WITH AP = " <<AP1 <<" == ZP = " <<ZP1
503  <<" REJECTED" <<G4endl;
504  G4cout <<"########################################"
505  <<"########################################"
506  <<G4endl;
507  }
508  return &theParticleChange;
509  }
510 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
511 
512  nucc_.ibproj = 1; // IBPROJ = 1
513  nucc_.ijproj = 1; // IJPROJ = 1
514  collis_.ijprox = 1; // IJPROX = 1
515  nucc_.ip = AP1; // IP = IP_P
516  nucc_.ipz = ZP1; // IPZ = IPZ_P
517  nucc_.it = AT1; // IT = IT_P
518  nucc_.itz = ZT1; // ITZ = ITZ_P
519  collis_.ijtar = 1; // IJTAR=1
520 //
521 //
522 // Note that epn is deliberately given units of GeV/nuc, because of the units
523 // used in DPMJET-II.5.
524 //
525  G4double epn = E / GeV;
526  G4double mpn = M / (GeV*AP); // Projectile mass per nucleon in GeV/(c2*nuc)
527 // G4double gamma = epn / mpn;
528 // G4double elab = epn * M / GeV; // ELAB = EPN * AMPRO_P
529 
530  diffra_.isingd = ISINGD;
531  user2_.isingx = ISINGX;
532  user2_.idubld = IDUBLD;
533  user2_.sdfrac = SDFRAC;
534 
535 // G4double amu_c2GeV = amu_c2 / GeV;
536  G4double ppn = std::sqrt((epn-mpn)*(epn+mpn));
537  // Units of GeV/(c*nuc)
538  // PPN = SQRT( (EPN-AMPROJ)*(EPN+AMPROJ) )
539 //
540 //
541 // Before setting the remainder of the variables for DPMJET-II.5, check for
542 // appropriate Glauber data. If the value returned is false, then set the change
543 // object to the source particle.
544 //
545  if (!(theGlauberDataSetHandler->SetCurrentGlauberDataSet(AP1,AT1,ppn))) {
546  theParticleChange.SetStatusChange(isAlive);
547  theParticleChange.SetEnergyChange(theTrack.GetKineticEnergy());
548  theParticleChange.SetMomentumChange(theTrack.Get4Momentum().vect().unit());
549  return &theParticleChange;
550  }
551  dtumat_.ntaxx[0] = AT1;
552  dtumat_.nztaxx[0] = ZT1;
553  dtumat_.nprxx[0] = AP1;
554  dtumat_.nzprxx[0] = ZP1;
555 //
556 //
557 // Set the remainder of the variables for DPMJET-II.5 FORTRAN.
558 //
559  nncms_.pproj = ppn; // PPROJ = PPN
560  nncms_.eproj = epn; // EPROJ = EPN
561  mpnt /= (AT * GeV); // Mass per nuclon of target in GeV/(c2*nuc)
562  nncms_.umo = std::sqrt(mpn*mpn + mpnt*mpnt + 2.0*mpnt*epn);
563  // UMO = SQRT( AMPROJ**2 + AMTAR**2 +2.D0*AMTAR*EPROJ )
564  // Note I believe this equation is only correct
565  // if the subsequent equations (for pTthr)
566  // needs the Ecm for the NUCLEON-NUCLEON system
567  user2_.cmener = nncms_.umo; // CMENER = UMO
568  collis_.s = nncms_.umo * nncms_.umo;
569  // SS = UMO**2
570  collis_.ptthr = 3.0;
571  if (strufu_.istrut == 1)
572  {
573  collis_.ptthr = 2.1 + 0.15*std::pow(std::log10(user2_.cmener/50.),3.0);
574  // PTTHR = 2.1D0+0.15D0*(LOG10(CMENER/50.))**3
575  }
576  else if (strufu_.istrut == 2)
577  {
578  collis_.ptthr = 2.5 + 0.12*std::pow(std::log10(user2_.cmener/50.),3.0);
579  // PTTHR = 2.5D0+0.12D0*(LOG10(CMENER/50.))**3
580  }
581  collis_.ptthr2 = collis_.ptthr; // PTTHR2 = PTTHR
582  nncms_.gamcm = (epn + mpnt) / nncms_.umo;
583  // GAMCM = (EPROJ+AMTAR)/UMO
584  // Note I believe this equation is only correct
585  // if the subsequent equations (for pTthr)
586  // need the Ecm for the NUCLEON-NUCLEON system
587  nncms_.bgcm = ppn / nncms_.umo; // PPROJ/UMO
588  nncms_.pcm = nncms_.gamcm*ppn - nncms_.bgcm*epn;
589  // PCM = GAMCM*PPROJ - BGCM*EPROJ
590  sigma_.sigsof = 37.8 * std::pow(collis_.s,0.076);
591  // ALFA = 1.076D0
592  // A = 37.8D0
593  // SIGSOF = A * SS**(ALFA-1.D0)
594  seasu3_.seasq = SEASQ; // SEASQ = 0.50D0
595  xseadi_.ssmima = SSMIMA; // SSMIMA = 1.201D0
596  xseadi_.ssmimq = xseadi_.ssmima * xseadi_.ssmima;
597  // SSMIMQ = SSMIMA**2
598  taufo_.taufor = TAUFOR;
599  taufo_.ktauge = KTAUGE;
600 
601  if ( theInitType == DEFAULT ) {
602  final_.ifinal = 0; // IFINAL = 1
603  evappp_.ievap = 0; // IEVAP = 0
604  parevt_.levprt = LTRUE; // LEVPRT = .FALSE.
605  parevt_.ilvmod = 1; // ILVMOD = 1
606  parevt_.ldeexg = LFALSE; // LDEEXG = .FALSE.
607  parevt_.lheavy = LFALSE; // LHEAVY = .FALSE.
608  frbkcm_.lfrmbk = LFALSE; // LFRMBK = .FALSE.
609  inpflg_.ifiss = 0; // IFISS = 0
610  } else if ( theInitType == CORSIKA ) {
611  final_.ifinal = 0; // IFINAL = 1
612  evappp_.ievap = 0; // IEVAP = 0
613  parevt_.levprt = LTRUE; // LEVPRT = .FALSE.
614  parevt_.ilvmod = 1; // ILVMOD = 1
615  parevt_.ldeexg = LFALSE; // LDEEXG = .FALSE.
616  parevt_.lheavy = LTRUE; // LHEAVY = .FALSE.
617  frbkcm_.lfrmbk = LFALSE; // LFRMBK = .FALSE.
618  inpflg_.ifiss = 0; // IFISS = 0
619  }
620  else if ( theInitType == DPM2_5 ) {
621  final_.ifinal = 0; // IFINAL = 0
622  evappp_.ievap = 0; // IEVAP = 0
623  parevt_.levprt = LTRUE; // LEVPRT = .TRUE. NOTE: THIS IS AT ODDS WITH WHAT'S IN DPMJET-II.5, BUT IF NOT SET, ALL EVENTS GET REJECTED.
624  parevt_.ilvmod = 1; // ILVMOD = 1
625  parevt_.ldeexg = LFALSE; // LDEEXG = .FALSE.
626  parevt_.lheavy = LFALSE; // LHEAVY = .FALSE.
627  frbkcm_.lfrmbk = LFALSE; // LFRMBK = .FALSE.
628  inpflg_.ifiss = 0; // IFISS = 0
629  }
630  else if ( theInitType == DPM3 ) {
631  final_.ifinal = 0; // IFINAL = 0
632  evappp_.ievap = 0; // IEVAP = 0
633  parevt_.levprt = LTRUE; // LEVPRT = .TRUE. NOTE: THIS IS AT ODDS WITH WHAT'S IN DPMJET-II.5, BUT IF NOT SET, ALL EVENTS GET REJECTED.
634  parevt_.ilvmod = 1; // ILVMOD = 1
635  parevt_.ldeexg = LFALSE; // LDEEXG = .FALSE.
636  parevt_.lheavy = LFALSE; // LHEAVY = .FALSE.
637  frbkcm_.lfrmbk = LFALSE; // LFRMBK = .FALSE.
638  inpflg_.ifiss = 0; // IFISS = 0
639  }
640  xsecpt_.ptcut = collis_.ptthr; // PTCUT = PTTHR
641 
642  G4double dsig1[maxpro+1];
643  csj1mi_ (&xsecpt_.ptcut, &dsig1[0]); // CALL CSJ1MI(PTCUT,DSIG1)
644  xsecpt_.dsigh = dsig1[0]; // SIG1 = DSIG1(0)
645  // DSIGH = SIG1
646  G4int i = 0;
647  G4double pt = 0.0;
648  samppt_ (&i,&pt); // SAMPPT(0,PT)
649  collap_.s3 = collis_.s; // S3 = SS
650  collap_.ijproj1 = collis_.ijprox; // IJPROJ1 = IJPROX
651  collap_.ijtar1 = collis_.ijtar; // IJTAR1 = IJTAR
652  collap_.ptthr1 = collis_.ptthr; // PTTHR1 = PTTHR
653  collap_.iophrd1 = collis_.iophrd; // IOPHRD1 = IOPHRD
654  collap_.ijprlu1 = collis_.ijprlu; // IJPRLU1 = IJPRLU
655  collap_.ijtalu1 = collis_.ijtalu; // IJTALU1 = IJTALU
656  collap_.ptthr3 = collis_.ptthr2; // PTTHR3 = PTTHR2
657 
658  G4int iiipro = nucc_.ijproj; // IIPROJ = IJPROJ & IIIPRO = IIPROJ
659  G4int iiitar = nucc_.ijtarg; // IITARG = IJTARG
660  G4int kkmat = 1;
661  G4int nhkkh1 = 1; // NHKKH1 = 1
662 
663  for (i=0; i<8; i++) user1_.projty[i] = paname_.btype[iiipro][i];
664  // PROJTY=BTYPE(IPROJ)
665  G4int irej = 1;
666  G4int evtcnt = 0;
667 
668  do {
669 /* bufueh_.annvv = 0.001; // ANNVV = 0.001
670  bufueh_.annss = 0.001; // ANNSS = 0.001
671  bufueh_.annsv = 0.001; // ANNSV = 0.001
672  bufueh_.annvs = 0.001; // ANNVS = 0.001
673  bufueh_.anncc = 0.001; // ANNCC = 0.001
674  bufueh_.anndv = 0.001; // ANNDV = 0.001
675  bufueh_.annvd = 0.001; // ANNVD = 0.001
676  bufueh_.annds = 0.001; // ANNDS = 0.001
677  bufueh_.annsd = 0.001; // ANNSD = 0.001
678  bufueh_.annhh = 0.001; // ANNHH = 0.001
679  bufueh_.annzz = 0.001; // ANNZZ = 0.001
680  bufueh_.anndi = 0.001; // ANNDI = 0.001
681  bufueh_.annzd = 0.001; // ANNZD = 0.001
682  bufueh_.anndz = 0.001; // ANNDZ = 0.001
683  bufueh_.ptvv = 0.0; // PTVV = 0.
684  bufueh_.ptss = 0.0; // PTSS = 0.
685  bufueh_.ptsv = 0.0; // PTSV = 0.
686  bufueh_.ptvs = 0.0; // PTVS = 0.
687  bufueh_.ptcc = 0.0; // PTCC = 0.
688  bufueh_.ptdv = 0.0; // PTDV = 0.
689  bufueh_.ptvd = 0.0; // PTVD = 0.
690  bufueh_.ptds = 0.0; // PTDS = 0.
691  bufueh_.ptsd = 0.0; // PTSD = 0.
692  bufueh_.pthh = 0.0; // PTHH = 0.
693  bufueh_.ptzz = 0.0; // PTZZ = 0.
694  bufueh_.ptdi = 0.0; // PTDI = 0.
695  bufueh_.ptzd = 0.0; // PTZD = 0.
696  bufueh_.ptdz = 0.0; // PTDZ = 0.
697  bufueh_.eevv = 0.0; // EEVV = 0.
698  bufueh_.eess = 0.0; // EESS = 0.
699  bufueh_.eesv = 0.0; // EESV = 0.
700  bufueh_.eevs = 0.0; // EEVS = 0.
701  bufueh_.eecc = 0.0; // EECC = 0.
702  bufueh_.eedv = 0.0; // EEDV = 0.
703  bufueh_.eevd = 0.0; // EEVD = 0.
704  bufueh_.eeds = 0.0; // EEDS = 0.
705  bufueh_.eesd = 0.0; // EESD = 0.
706  bufueh_.eehh = 0.0; // EEHH = 0.
707  bufueh_.eezz = 0.0; // EEZZ = 0.
708  bufueh_.eedi = 0.0; // EEDI = 0.
709  bufueh_.eezd = 0.0; // EEZD = 0.
710  bufueh_.eedz = 0.0; // EEDZ = 0.
711  ncouch_.acouvv = 0.0; // ACOUVV = 0.
712  ncouch_.acouss = 0.0; // ACOUSS = 0.
713  ncouch_.acousv = 0.0; // ACOUSV = 0.
714  ncouch_.acouvs = 0.0; // ACOUVS = 0.
715  ncouch_.acouzz = 0.0; // ACOUZZ = 0.
716  ncouch_.acouhh = 0.0; // ACOUHH = 0.
717  ncouch_.acouds = 0.0; // ACOUDS = 0.
718  ncouch_.acousd = 0.0; // ACOUSD = 0.
719  ncouch_.acoudz = 0.0; // ACOUDZ = 0.
720  ncouch_.acouzd = 0.0; // ACOUZD = 0.
721  ncouch_.acoudi = 0.0; // ACOUDI = 0.
722  ncouch_.acoudv = 0.0; // ACOUDV = 0.
723  ncouch_.acouvd = 0.0; // ACOUVD = 0.
724  ncouch_.acoucc = 0.0; // ACOUCC = 0.*/
725  if (evtcnt > 0)
726  G4cout <<"REJECTED KKINC EVENT. RETRY # = " <<evtcnt <<G4endl;
727 //
728 //
729 // Generate an event using DPMJET II-5. NOTE this is a call to dpmevt, and
730 // could possibly be replaced by call to kkinc_. After the call,
731 // ResetCurrentGlauberDataSet is used to delete any temporary Glauber data
732 // set if one had to be set up.
733 //
734  //G4cout << "Call to kkinc_" << G4endl;
735  kkinc_ (&epn, &AT1, &ZT1, &AP1, &ZP1, &iiipro, &kkmat, &iiitar, &nhkkh1,
736  &irej);
737  // CALL KKINC(EPN,IIT,IITZ,IIP,IIPZ,IIPROJ,KKMAT,
738  // * IITARG,NHKKH1,IREJ)
739 // dpmevt_ (&elabt, &iiipro, &AP1, &ZP1, &AT1, &ZT1, &kkmat, &nhkkh1);
740  } while (irej == 1 && ++evtcnt <100);
741  //G4cout << "Call to reset G-data" << G4endl;
742 
743  theGlauberDataSetHandler->ResetCurrentGlauberDataSet();
744 //
745 //
746 // If the event has been rejected more than 100 times, then set the track
747 // as still active and return to the calling routine.
748 //
749  if (irej == 1) {
750  theParticleChange.SetStatusChange(isAlive);
751  theParticleChange.SetEnergyChange(theTrack.GetKineticEnergy());
752  theParticleChange.SetMomentumChange(theTrack.Get4Momentum().vect().unit());
753  if (verboseLevel >= 2) {
754  G4cout <<"Event rejected and original track maintained" <<G4endl;
755  G4cout <<"########################################"
756  <<"########################################"
757  <<G4endl;
758  }
759  return &theParticleChange;
760  }
761 //
762 //
763 // Determine number of final state particles (including nuclear fragments,
764 // and load into the particle change if you can identify the particles.
765 //
766  G4int n = hkkevt_.nhkk;
767  G4int m = 0;
768  G4Fragment *fragment = 0;
769  if (verboseLevel >= 2) DumpVerboseInformation1 (n);
770 //
771 //
772 // Now go through each of the secondaries and add to theParticleChange.
773 //
774  for (G4int i=0; i<n; i++)
775  {
776  if (hkkevt_.isthkk[i]==1 || hkkevt_.isthkk[i]==-1)
777  {
778 //
779 // Particle is a final state secondary and not a nucleus.
780 // Determine what this secondary particle is, and if valid, load dynamic
781 // parameters.
782 //
783  G4ParticleDefinition* theParticle =
784  theParticleTable->FindParticle(hkkevt_.idhkk[i]);
785  if (theParticle)
786  {
787  G4double px = hkkevt_.phkk[i][0] * GeV;
788  G4double py = hkkevt_.phkk[i][1] * GeV;
789  G4double pz = hkkevt_.phkk[i][2] * GeV;
790  G4double et = hkkevt_.phkk[i][3] * GeV;
791 // G4LorentzVector lv = transformToLab * G4LorentzVector(px,py,pz,et);
792  G4LorentzVector lv = G4LorentzVector(px,py,pz,et);
793 
794  G4DynamicParticle *theDynamicParticle =
795  new G4DynamicParticle(theParticle,lv);
796  theParticleChange.AddSecondary (theDynamicParticle);
797 
798  if (verboseLevel >= 2)
799  DumpVerboseInformation2 (theParticle->GetParticleName(),
800  lv.vect(), et, theDynamicParticle->GetKineticEnergy(),pP);
801  }
802  }
803  else if (hkkevt_.idhkk[i]==80000 && hkkevt_.isthkk[i]==1001)
804  {
805 //
806 //
807 // Particle is a secondary nucleus. Determine the details of the nuclear
808 // fragment prior to de-excitation. (Note that the 1 eV in the total energy
809 // is a safety factor to avoid any possibility of negative rest mass energy.)
810 // Note also that we don't full trust the energy provided by the DPMJET-II.5,
811 // and there it's based on the Geant4-determined ion rest-mass.
812 //
813  G4int nucA = extevt_.idres[i];
814  G4int nucZ = extevt_.idxres[i];
815  if (nucA>0 && nucZ>0) {
816  m++;
817  fragment = 0;
818  G4double px = hkkevt_.phkk[i][0] * GeV;
819  G4double py = hkkevt_.phkk[i][1] * GeV;
820  G4double pz = hkkevt_.phkk[i][2] * GeV;
821  G4double ionMass = theIonTable->GetIonMass(nucZ,nucA);
822 // GetIonMass(nucZ,nucA) + nucex[i]; // check how to get this energy
823  G4double dpmMass = hkkevt_.phkk[i][4] * GeV;
824  if (dpmMass > ionMass) ionMass = dpmMass;
825  G4double et = std::sqrt(px*px + py*py + pz*pz + ionMass*ionMass);
826 // G4LorentzVector lv = transformToLab * G4LorentzVector(px,py,pz,et+1.0*eV);
827  G4LorentzVector lv = G4LorentzVector(px,py,pz,et+1.0*eV);
828  fragment = new G4Fragment(nucA, nucZ, lv);
829  if (verboseLevel >= 2)
830  DumpVerboseInformation3 (m, nucA, nucZ, lv.vect(), et, et-ionMass, pP);
831 //
832 //
833 // Now we can decay the nuclear fragment if present. The secondaries are
834 // collected and boosted as well. The priority is to use a pre-compound
835 // de-excitation, otherwise the standard excitation-handler is used.
836 //
837  if (fragment && (thePreComp || theExcitationHandler))
838  {
839  G4ReactionProductVector *products = 0;
840  if (thePreComp && fragment->GetA() > 1.5)
841  products = thePreComp->DeExcite(*fragment);
842  else
843  products = theExcitationHandler->BreakItUp(*fragment);
844  delete fragment;
845  fragment = 0;
846  G4ReactionProductVector::iterator iter;
847  for (iter = products->begin(); iter != products->end(); ++iter)
848  {
849  G4DynamicParticle *secondary =
850  new G4DynamicParticle((*iter)->GetDefinition(),
851  (*iter)->GetTotalEnergy(), (*iter)->GetMomentum());
852  theParticleChange.AddSecondary (secondary);
853  G4String particleName = (*iter)->GetDefinition()->GetParticleName();
854  delete (*iter);
855  if (verboseLevel >= 2) {
856  if (particleName.find("[",0) < particleName.size())
857  DumpVerboseInformation4 (m, particleName, secondary->GetMomentum(),
858  secondary->GetTotalEnergy(), secondary->GetKineticEnergy(), pP);
859  else DumpVerboseInformation2 (particleName, secondary->GetMomentum(),
860  secondary->GetTotalEnergy(), secondary->GetKineticEnergy(), pP);
861  }
862  }
863  delete products;
864  }
865  if (fragment != 0)
866  {
867 //
868 //
869 // Add the excited fragment to the product vector. Note that this is temporary
870 // since we should at the atomic excitation to strip all electrons ... i.e. it's
871 // actually a bare nucleus not an atom in the ground state.
872 //
873  G4ParticleDefinition *theParticleDefinition = theIonTable->
874  GetIon(nucZ,nucA);
875  G4DynamicParticle *theDynamicParticle =
876  new G4DynamicParticle(theParticleDefinition,lv);
877  theParticleChange.AddSecondary (theDynamicParticle);
878  delete fragment;
879  fragment = 0;
880  }
881  }
882  }
883  }
884 
885  if (verboseLevel >= 3) {
886 //
887 //
888 // Calculate and display the energy and momenta before and after the collision.
889 // Everything is calculated for the lab frame.
890 //
891  G4double TotalEPost = 0.0;
892  G4ThreeVector TotalPPost;
893  G4double charge = 0.0;
894  G4int baryon = 0;
895  G4int lepton = 0;
896 // G4int parity = 0;
897  G4int nSecondaries = theParticleChange.GetNumberOfSecondaries();
898  for (G4int j=0; j<nSecondaries; j++) {
899  TotalEPost += theParticleChange.GetSecondary(j)->
900  GetParticle()->GetTotalEnergy();
901  TotalPPost += theParticleChange.GetSecondary(j)->
902  GetParticle()->GetMomentum();
903  G4ParticleDefinition *theParticle = theParticleChange.GetSecondary(j)->
904  GetParticle()->GetDefinition();
905  charge += theParticle->GetPDGCharge();
906  baryon += theParticle->GetBaryonNumber();
907  lepton += theParticle->GetLeptonNumber();
908 // parity += theParticle->GetPDGiParity();
909  }
910  G4cout <<"----------------------------------------"
911  <<"----------------------------------------"
912  <<G4endl;
913  G4cout <<"Total energy before collision = " <<TotalEPre/MeV
914  <<" MeV" <<G4endl;
915  G4cout <<"Total energy after collision = " <<TotalEPost/MeV
916  <<" MeV" <<G4endl;
917  G4cout <<"Total momentum before collision = " <<pP/MeV
918  <<" MeV/c" <<G4endl;
919  G4cout <<"Total momentum after collision = " <<TotalPPost/MeV
920  <<" MeV/c" <<G4endl;
921  if (verboseLevel >= 4) {
922  G4cout <<"Total charge before collision = " <<(ZP+ZT)*eplus
923  <<G4endl;
924  G4cout <<"Total charge after collision = " <<charge
925  <<G4endl;
926  G4cout <<"Total baryon number before collision = "<<AP+AT
927  <<G4endl;
928  G4cout <<"Total baryon number after collision = "<<baryon
929  <<G4endl;
930  G4cout <<"Total lepton number before collision = 0"
931  <<G4endl;
932  G4cout <<"Total lepton number after collision = "<<lepton
933  <<G4endl;
934  }
935  G4cout <<"----------------------------------------"
936  <<"----------------------------------------"
937  <<G4endl;
938  }
939 
940  if (verboseLevel >= 2)
941  G4cout <<"########################################"
942  <<"########################################"
943  <<G4endl;
944 
945  return &theParticleChange;
946 }
948 //
949 // SetNoDeexcitation
950 //
951 // Deletes an exiting de-excitation handlers and zeros the pointer. This
952 // allows the simulation to run without any de-excitation, or just pre-compound
953 // model only.
954 //
955 // Note that you need to separately SetNoPreCompoundModel, if you REALLY don't
956 // want any nuclear de-excitation. But please only do this to understand the
957 // contribution of de-excitation/pre-equilibrium to the full simulation.
958 // Running without de-excitation or pre-compound is physically unrealistic!
959 //
960 //
961 void G4DPMJET2_5Model::SetNoDeexcitation ()
962 {
963  if (theExcitationHandler)
964  {
965  delete theExcitationHandler;
966  theExcitationHandler = 0;
967  }
968 }
970 //
971 // SetNoPreCompoundModel
972 //
973 // Deletes an exiting pre-equilibrium model and zeros the pointer. This
974 // allows the simulation to run without any pre-compound, or just -de-excitation
975 // model only.
976 //
977 // Note that you need to separately SetNoDeexcitation, if you REALLY don't want
978 // any nuclear de-excitation. But please only do this to understand the
979 // contribution of de-excitation/pre-equilibrium to the full simulation.
980 // Running without de-excitation or pre-compound is physically unrealistic!
981 //
982 //
983 void G4DPMJET2_5Model::SetNoPreCompoundModel ()
984 {
985  if (thePreComp)
986  {
987  delete thePreComp;
988  thePreComp = 0;
989  }
990 }
992 //
993 // SetDefaultDeexcitation
994 //
995 // Note that this is used by the default constructor, as well as can be called
996 // directly by the user.
997 //
998 void G4DPMJET2_5Model::SetDefaultDeexcitation ()
999 {
1000 // SetNoDeexcitation();
1001 
1002  theExcitationHandler = new G4ExcitationHandler;
1003  G4Evaporation * theEvaporation = new G4Evaporation;
1004  G4FermiBreakUp * theFermiBreakUp = new G4FermiBreakUp;
1005  G4PhotonEvaporation* thePhotonEvap = new G4PhotonEvaporation;
1006  theExcitationHandler->SetEvaporation(theEvaporation);
1007  theExcitationHandler->SetFermiModel(theFermiBreakUp);
1008  theExcitationHandler->SetMaxAandZForFermiBreakUp(17, 9);
1009  theExcitationHandler->SetFermiModel(theFermiBreakUp);
1010  theExcitationHandler->SetPhotonEvaporation(thePhotonEvap);
1011 }
1013 //
1014 // SetDefaultPreCompoundModel
1015 //
1016 // Note that this is used by the default constructor, as well as can be called
1017 // directly by the user.
1018 //
1019 void G4DPMJET2_5Model::SetDefaultPreCompoundModel ()
1020 {
1021 // SetNoPreCompoundModel();
1022 
1023  G4ExcitationHandler *anExcitationHandler = new G4ExcitationHandler;
1024  G4Evaporation * theEvaporation = new G4Evaporation;
1025  G4FermiBreakUp * theFermiBreakUp = new G4FermiBreakUp;
1026  G4PhotonEvaporation* thePhotonEvap = new G4PhotonEvaporation;
1027  anExcitationHandler->SetEvaporation(theEvaporation);
1028  anExcitationHandler->SetFermiModel(theFermiBreakUp);
1029  anExcitationHandler->SetMaxAandZForFermiBreakUp(17, 9);
1030  anExcitationHandler->SetPhotonEvaporation(thePhotonEvap);
1031 
1032  thePreComp = new G4PreCompoundModel(anExcitationHandler);
1033 }
1035 //
1036 // SetVerboseFortranOutput
1037 //
1038 G4bool G4DPMJET2_5Model::SetVerboseFortranOutput (const G4String filename)
1039 {
1040  g4dpmjet_close_fort6_ ();
1041  if (filename == "" || filename == "stdo" ||
1042  filename == "stdout" || filename == "std::out" )
1043  {
1044  verboseFortranFile = "std::out";
1045  return true;
1046  }
1047  else
1048  {
1049  G4int namelen = filename.length();
1050  char *ptr = new char[namelen+1];
1051  filename.copy(ptr,namelen,0);
1052  ptr[namelen] = '\0';
1053 // char *ptr = 0;
1054 // ptr = const_cast<char*> (filename.c_str());
1055  ftnlogical opened = LFALSE;
1056  g4dpmjet_open_fort6_ (&namelen, &opened, ptr);
1057  delete [] ptr;
1058  if (opened == LTRUE) {
1059  verboseFortranFile = filename;
1060  return true;
1061  } else {
1062  verboseFortranFile = "std::out";
1063  return false;
1064  }
1065  }
1066 }
1068 //
1069 // PrintWelcomeMessage
1070 //
1071 void G4DPMJET2_5Model::PrintWelcomeMessage () const
1072 {
1073  G4cout <<G4endl;
1074  G4cout <<" *****************************************************************"
1075  <<G4endl;
1076  G4cout <<" Interface to DPMJET2.5 for nuclear-nuclear interactions activated"
1077  <<G4endl;
1078  G4cout <<" Version number : 00.00.0B File date : 23/05/08" <<G4endl;
1079  G4cout <<" (Interface written by QinetiQ Ltd for the European Space Agency)"
1080  <<G4endl;
1081  G4cout <<G4endl;
1082  G4cout <<" Initialisation of DPMJET-II.5 variables will be according to "
1083  <<theInitType <<G4endl;
1084  G4cout <<" *****************************************************************"
1085  <<G4endl;
1086  G4cout << G4endl;
1087 
1088  return;
1089 }
1091 //
1092 // DumpVerboseInformation1
1093 //
1094 // Dumps raw information about the DPMJET-II.5 simulation if verbosity set
1095 // to 4 or more.
1096 //
1097 void G4DPMJET2_5Model::DumpVerboseInformation1 (const G4int n) const
1098 {
1099  G4cout <<"----------------------------------------"
1100  <<"----------------------------------------" <<G4endl;
1101  G4cout <<n <<" INTERMEDIATE AND FINAL-STATE SECONDARIES PRODUCED" <<G4endl;
1102  if (verboseLevel >= 4)
1103  {
1104  G4cout <<"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
1105  <<"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" <<G4endl;
1106  G4cout <<"ORIGINAL DPMJET-II.5 OUTPUT FOR EVENT:" <<G4endl;
1107  G4cout <<"Note that (1) the particles are yet to be transformed according"
1108  <<G4endl;
1109  G4cout <<" to incident particle direction" <<G4endl;
1110  G4cout <<" (2) the units of energy, momentum and mass are GeV,"
1111  <<G4endl;
1112  G4cout <<" GeV/c and GeV/c^2 respectively" <<G4endl;
1113  G4cout <<" I"
1114  <<" ISTHKK"
1115  <<" IDHKK"
1116  <<" IDRES"
1117  <<" IDXRES"
1118  <<" PX"
1119  <<" PY"
1120  <<" PZ"
1121  <<" TOTAL ENERGY"
1122  <<" MASS"
1123  <<G4endl;
1124  for (G4int i=0; i<n; i++)
1125  {
1126  G4cout.unsetf(std::ios::scientific);
1127  G4cout.setf(std::ios::fixed|std::ios::right|std::ios::adjustfield);
1128  G4cout.precision(0);
1129  G4cout <<std::setw(5) <<i
1130  <<std::setw(10) <<hkkevt_.isthkk[i]
1131  <<std::setw(10) <<hkkevt_.idhkk[i]
1132  <<std::setw(10) <<extevt_.idres[i]
1133  <<std::setw(10) <<extevt_.idxres[i];
1134  G4cout.unsetf(std::ios::fixed);
1135  G4cout.setf(std::ios::scientific|std::ios::right|std::ios::adjustfield);
1136  G4cout.precision(7);
1137  G4cout <<std::setw(15) <<hkkevt_.phkk[i][0]
1138  <<std::setw(15) <<hkkevt_.phkk[i][1]
1139  <<std::setw(15) <<hkkevt_.phkk[i][2]
1140  <<std::setw(15) <<hkkevt_.phkk[i][3]
1141  <<std::setw(15) <<hkkevt_.phkk[i][4]
1142  <<G4endl;
1143  }
1144  G4cout <<"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
1145  <<"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" <<G4endl;
1146  }
1147  G4cout.setf(std::ios::fixed);
1148  G4cout <<" THE FOLLOWING LISTS ONLY THE FINAL-STATE SECONDARIES" <<G4endl;
1149 }
1151 //
1152 // Initialise
1153 // This is intended to do exactly what the main program and subroutine DMINIT
1154 // attempt to do with variable initialisation. i've included most of the
1155 // FORTRAN source code for reference. Note that I'm certain that many of these
1156 // lines relate to tallying of events for output by the standalone version
1157 // of DPMJET-II.5, but for the moment, I'd rather initialise those variables,
1158 // just in case not doing so has an adverse effect.
1159 //
1160 void G4DPMJET2_5Model::Initialise ()
1161 {
1162 //
1163 //
1164 // This first line is intended to make sure the block data statements are
1165 // executed, since we're not running from a FORTRAN main program.
1166 //
1167  g4dpmjet_initialise_block_data_ ();
1168 //
1169 //
1170  dpar_.aam[4] = 0.001; // AAM(5)=0.001D0
1171  dpar_.aam[5] = 0.001; // AAM(6)=0.001D0
1172  dpar_.aam[132] = 0.001; // AAM(133)=0.001D0
1173  dpar_.aam[133] = 0.001; // AAM(134)=0.001D0
1174  dpar_.aam[134] = 0.001; // AAM(135)=0.001D0
1175  dpar_.aam[135] = 0.001; // AAM(136)=0.001D0
1176 
1177  vxsvd_.nxsp = 0; // NXSP=0
1178  vxsvd_.nxst = 0; // NXST=0
1179  vxsvd_.nxsap = 0; // NXSAP=0
1180  vxsvd_.nxsat = 0; // NXSAT=0
1181  vxsvd_.nxvp = 0; // NXVP=0
1182  vxsvd_.nxvt = 0; // NXVT=0
1183  vxsvd_.nxdp = 0; // NXDP=0
1184  vxsvd_.nxdt = 0; // NXDT=0
1185 
1186  for (G4int i=0; i<50; i++)
1187  {
1188  vxsvd_.vxsp[i] = 1.0E-08; // VXSP(II)=1.D-8
1189  vxsvd_.vxst[i] = 1.0E-08; // VXST(II)=1.D-8
1190  vxsvd_.vxsap[i] = 1.0E-08; // VXSAP(II)=1.D-8
1191  vxsvd_.vxsat[i] = 1.0E-08; // VXSAT(II)=1.D-8
1192  vxsvd_.vxvp[i] = 1.0E-08; // VXVP(II)=1.D-8
1193  vxsvd_.vxvt[i] = 1.0E-08; // VXVT(II)=1.D-8
1194  vxsvd_.vxdp[i] = 1.0E-08; // VXDP(II)=1.D-8
1195  vxsvd_.vxst[i] = 1.0E-08; // VXST(II)=1.D-8
1196  }
1197 
1198  if (debug)
1199  {
1200 #ifdef G4VERBOSE
1201  if (GetVerboseLevel()>0) {
1202  G4cout <<"AT G4DPMJET2_5Model::Initialise:" <<G4endl;
1203  }
1204 #endif
1205  dprin_.ipri = debug_level; // IPRI = LEVLDB
1206  dprin_.ipev = debug_level; // IPEV = LEVLDB
1207  dprin_.ippa = debug_level; // IPPA = LEVLDB
1208  dprin_.ipco = debug_level; // IPCO = LEVLDB
1209  dprin_.init = debug_level; // INIT = LEVLDB
1210  dprin_.iphkk = debug_level; // IPHKK = LEVLDB
1211  pydat1_.mstu[25] = 10; // MSTU(26) = 10
1212  }
1213  else
1214  {
1215  dprin_.ipri = 0; // IPRI = 0
1216  dprin_.ipev = 0; // IPEV = 0
1217  dprin_.ippa = 0; // IPPA = 0
1218  dprin_.ipco =-2; // IPCO = -2
1219  dprin_.init = 0; // INIT = 0
1220  dprin_.iphkk = 0; // IPHKK = 0
1221  pydat1_.mstu[25] = 0; // MSTU(26) = 10
1222  }
1223 
1224  for (G4int i=0; i<7; i++)
1225  {
1226  diqrej_.idiqre[i] = 0;
1227  diqrej_.idiqrz[i] = 0;
1228  }
1229 
1230  for (G4int i=0; i<3; i++)
1231  {
1232  diqrej_.idvre[i] = 0;
1233  diqrej_.ivdre[i] = 0;
1234  diqrej_.idsre[i] = 0;
1235  diqrej_.isdre[i] = 0;
1236  diqrej_.idzre[i] = 0;
1237  diqrej_.izdre[i] = 0;
1238  }
1239 
1240  diqsum_.ndvuu = 0; // NDVUU = 0
1241  diqsum_.ndvus = 0; // NDVUS = 0
1242  diqsum_.ndvss = 0; // NDVSS = 0
1243  diqsum_.nvduu = 0; // NVDUU = 0
1244  diqsum_.nvdus = 0; // NVDUS = 0
1245  diqsum_.nvdss = 0; // NVDSS = 0
1246  diqsum_.ndsuu = 0; // NDSUU = 0
1247  diqsum_.ndsus = 0; // NDSUS = 0
1248  diqsum_.ndsss = 0; // NDSSS = 0
1249  diqsum_.nsduu = 0; // NSDUU = 0
1250  diqsum_.nsdus = 0; // NSDUS = 0
1251  diqsum_.nsdss = 0; // NSDSS = 0
1252  diqsum_.ndzuu = 0; // NDZUU = 0
1253  diqsum_.ndzus = 0; // NDZUS = 0
1254  diqsum_.ndzss = 0; // NDZSS = 0
1255  diqsum_.nzduu = 0; // NZDUU = 0
1256  diqsum_.nzdus = 0; // NZDUS = 0
1257  diqsum_.nzdss = 0; // NZDSS = 0
1258  diqsum_.nadvuu = 0; // NADVUU = 0
1259  diqsum_.nadvus = 0; // NADVUS = 0
1260  diqsum_.nadvss = 0; // NADVSS = 0
1261  diqsum_.navduu = 0; // NAVDUU = 0
1262  diqsum_.navdus = 0; // NAVDUS = 0
1263  diqsum_.navdss = 0; // NAVDSS = 0
1264  diqsum_.nadsuu = 0; // NADSUU = 0
1265  diqsum_.nadsus = 0; // NADSUS = 0
1266  diqsum_.nadsss = 0; // NADSSS = 0
1267  diqsum_.nasduu = 0; // NASDUU = 0
1268  diqsum_.nasdus = 0; // NASDUS = 0
1269  diqsum_.nasdss = 0; // NASDSS = 0
1270  diqsum_.nadzuu = 0; // NADZUU = 0
1271  diqsum_.nadzus = 0; // NADZUS = 0
1272  diqsum_.nadzss = 0; // NADZSS = 0
1273  diqsum_.nazduu = 0; // NAZDUU = 0
1274  diqsum_.nazdus = 0; // NAZDUS = 0
1275  diqsum_.nazdss = 0; // NAZDSS = 0
1276  hdjase_.nhse1 = 0; // NHSE1 = 0
1277  hdjase_.nhse2 = 0; // NHSE2 = 0
1278  hdjase_.nhse3 = 0; // NHSE3 = 0
1279  hdjase_.nhase1 = 0; // NHASE1 = 0
1280  hdjase_.nhase2 = 0; // NHASE2 = 0
1281  hdjase_.nhase3 = 0; // NHASE3 = 0
1282 //
1283 //
1284 // Parton pt distribution.
1285 //
1286  G4int i = 1;
1287  G4double pt1 = 0.0;
1288  G4double pt2 = 0.0;
1289  G4int ipt = 0;
1290  G4int nevt = 0;
1291  parpt_ (&i,&pt1,&pt2,&ipt,&nevt); // CALL PARPT(1,PT1,PT2,IPT,NEVT)
1292 //
1293 //
1294 // Initialise BAMJET, DECAY and HADRIN.
1295 //
1296  ddatar_ (); // CALL DDATAR
1297  dhadde_ (); // CALL DHADDE
1298  dchant_ (); // CALL DCHANT
1299  dchanh_ (); // CALL DCHANH
1300 
1301  G4double epn = 0.0;
1302  G4double ppn = 0.0;
1303  defaul_ (&epn,&ppn); // CALL DEFAUL(EPN,PPN)
1304  defaux_ (&epn,&ppn); // CALL DEFAUX(EPN,PPN)
1305 
1306  coulo_.icoul = 1; // ICOUL = 1
1307  nuclea_.icoull = 1; // ICOULL = 1
1308  edens_.ieden = 0; // IEDEN = 0
1309  dprin_.itopd = 0; // ITOPD = 0
1310 
1311  if ( theInitType == DEFAULT ) {
1312  TAUFOR = 5.0E+00;
1313  KTAUGE = 25;
1314  } else if ( theInitType == CORSIKA ) {
1315  TAUFOR = 5.0E+00; // TAUFOR = 5.D0
1316  KTAUGE = 25; // KTAUGE = 25
1317  } else if ( theInitType == DPM2_5 ) {
1318  TAUFOR = 105.0E+00; // TAUFOR = 105.D0
1319  KTAUGE = 10; // KTAUGE = 10
1320  } else if ( theInitType == DPM3 ) {
1321  TAUFOR = 3.5E+00;
1322  KTAUGE = 10;
1323  }
1324  taufo_.taufor = TAUFOR;
1325  taufo_.ktauge = KTAUGE;
1326 
1327  ITAUVE = 1; // ITAUVE = 1
1328  taufo_.itauve = ITAUVE;
1329  taufo_.incmod = 1; // INCMOD = 1
1330 //
1331 //
1332 // Definition of soft quark distributions, Fermi, Pauli
1333 //
1334  G4double xseaco = 1.0; // XSEACO = 1.00D0
1335  xseadi_.xseacu = 1.05 - xseaco; // XSEACU = 1.05D0-XSEACO
1336 
1337  if ( theInitType == DEFAULT ) {
1338  UNON = 3.50;
1339  UNOM = 1.11;
1340  UNOSEA = 5.0;
1341  droppt_.fermp = LTRUE;
1342  nucimp_.fermod = 0.6;
1343  } else if ( theInitType == CORSIKA ) {
1344  UNON = 3.50; // UNON = 3.50D0
1345  UNOM = 1.11; // UNOM = 1.11D0
1346  UNOSEA = 5.0; // UNOSEA = 5.0D0
1347  droppt_.fermp = LTRUE; // FERMP = .TRUE.
1348  nucimp_.fermod = 0.6; // FERMOD = 0.6D0
1349  } else if ( theInitType == DPM2_5 ) {
1350  UNON = 3.50; // UNON = 3.50D0
1351  UNOM = 1.11; // UNOM = 1.11D0
1352  UNOSEA = 5.0; // UNOSEA = 5.0D0
1353  droppt_.fermp = LTRUE; // FERMP = .TRUE.
1354  nucimp_.fermod = 0.6; // FERMOD = 0.6D0
1355  } else if ( theInitType == DPM3 ) {
1356  UNON = 2.00;
1357  UNOM = 1.5;
1358  UNOSEA = 5.0;
1359  droppt_.fermp = LTRUE;
1360  nucimp_.fermod = 0.55;
1361  }
1362  xseadi_.unon = UNON;
1363  xseadi_.unom = UNOM;
1364  xseadi_.unosea = UNOSEA;
1365 
1366  nuclea_.fermdd = 0.6; // FERMDD = 0.6D0
1367  ferfor_.iferfo = 1; // IFERFO = 1
1368  dprin_.ipaupr = 0; // IPAUPR = 0
1369  droppt_.lpauli = LTRUE; // LPAULI = .TRUE.
1370 //
1371 //
1372 // Definition of cuts for x-sampling.
1373 //
1374  if ( theInitType == DEFAULT ) {
1375  CVQ = 1.8;
1376  CDQ = 2.0;
1377  CSEA = 0.5;
1378  SSMIMA = 0.9;
1379  } else if ( theInitType == CORSIKA ) {
1380  CVQ = 1.8; // CVQ = 1.8D0
1381  CDQ = 2.0; // CDQ = 2.0D0
1382  CSEA = 0.5; // CSEA = 0.5D0
1383  SSMIMA = 0.901; // SSMIMA = 0.901D0
1384  } else if ( theInitType == DPM2_5 ) {
1385  CVQ = 1.8; // CVQ = 1.8D0
1386  CDQ = 2.0; // CDQ = 2.0D0
1387  CSEA = 0.5; // CSEA = 0.5D0
1388  SSMIMA = 1.201; // SSMIMA = 1.201D0
1389  } else if ( theInitType == DPM3 ) {
1390  CVQ = 1.0;
1391  CDQ = 2.0;
1392  CSEA = 0.1;
1393  SSMIMA = 0.14;
1394  }
1395  xseadi_.cvq = CVQ;
1396  xseadi_.cdq = CDQ;
1397  xseadi_.csea = CSEA;
1398  xseadi_.ssmima = SSMIMA;
1399 
1400  xseadi_.ssmimq = xseadi_.ssmima * xseadi_.ssmima;
1401  // SSMIMQ = SSMIMA**2
1402  if ( theInitType == DEFAULT ) {
1403  VVMTHR = 0.0;
1404  } else if ( theInitType == CORSIKA ) {
1405  VVMTHR = 0.0; // VVMTHR = 0.D0
1406  } else if ( theInitType == DPM2_5 ) {
1407  VVMTHR = 0.0; // VVMTHR = 0.D0
1408  } else if ( theInitType == DPM3 ) {
1409  VVMTHR = 2.0;
1410  }
1411  xseadi_.vvmthr = VVMTHR;
1412 //
1413 //
1414 // There is a final call. Set recombin, seasu3, coninpt, allpart and interdpm
1415 //
1416  final_.ifinal = 0; // IFINAL = 0
1417  recom_.irecom = 0; // IRECOM = 0
1418  seadiq_.lseadi = LTRUE; // LSEADI = .TRUE.
1419 
1420  if ( theInitType == DEFAULT ) {
1421  SEASQ = 0.5;
1422  MKCRON = 1;
1423  CRONCO = 0.64;
1424  } else if ( theInitType == CORSIKA ) {
1425  SEASQ = 0.5; // SEASQ = 0.50D0
1426  MKCRON = 0; // MKCRON = 0
1427  CRONCO = 0.0; // CRONCO = 0.00D0
1428  } else if ( theInitType == DPM2_5 ) {
1429  SEASQ = 0.5; // SEASQ = 0.50D0
1430  MKCRON = 1; // MKCRON = 1
1431  CRONCO = 0.64; // CRONCO = 0.64D0
1432  } else if ( theInitType == DPM3 ) {
1433  SEASQ = 1.0;
1434  MKCRON = 1;
1435  CRONCO = 0.64;
1436  }
1437  seasu3_.seasq = SEASQ;
1438  cronin_.mkcron = MKCRON;
1439  cronin_.cronco = CRONCO;
1440 
1441  droppt_.ihada = LTRUE; // IHADA = .TRUE.
1442 
1443 // inxdpm_.intdpm = 0; // INTDPM = 0
1444 // Note FORTRAN initialises the variable IROEH = 0, but this isn't used
1445 // nor is it contained within a common block.
1446 //
1447 //
1448 // Definition for popcork, casadiqu, popcorse
1449 //
1450  popcck_.pdbck = 0.0; // PDBCK = 0.D0
1451  popcck_.ijpock = 0; // IJPOCK = 0
1452  casadi_.icasad = 1; // ICASAD = 1
1453  casadi_.casaxx = 0.05; // CASAXX = 0.05D0 ! corrected Nov. 2001
1454  popcck_.pdbse = 0.45; // PDBSE = 0.45D0 ! with baryon stopping
1455  popcck_.pdbseu = 0.45; // PDBSEU = 0.45D0 ! with baryon stopping
1456 // C PDBSE = 0.D0 ! without baryon stopping
1457 // C PDBSEU = 0.D0 ! without baryon stopping
1458  popcck_.irejck = 0; // IREJCK = 0
1459  popcck_.irejse = 0; // IREJSE = 0
1460  popcck_.irejs3 = 0; // IREJS3 = 0
1461  popcck_.irejs0 = 0; // IREJS0 = 0
1462  popcck_.ick4 = 0; // ICK4 = 0
1463  popcck_.ise4 = 0; // ISE4 = 0
1464  popcck_.ise43 = 0; // ISE43 = 0
1465  popcck_.ihad4 = 0; // IHAD4 = 0
1466  popcck_.ick6 = 0; // ICK6 = 0
1467  popcck_.ise6 = 0; // ISE6 = 0
1468  popcck_.ise63 = 0; // ISE63 = 0
1469  popcck_.ihad6 = 0; // IHAD6 = 0
1470  popcck_.irejsa = 0; // IREJSA = 0
1471  popcck_.ireja3 = 0; // IREJA3 = 0
1472  popcck_.ireja0 = 0; // IREJA0 = 0
1473  popcck_.isea4 = 0; // ISEA4 = 0
1474  popcck_.isea43 = 0; // ISEA43 = 0
1475  popcck_.ihada4 = 0; // IHADA4 = 0
1476  popcck_.isea6 = 0; // ISEA6 = 0
1477  popcck_.isea63 = 0; // ISEA63 = 0
1478  popcck_.ihada6 = 0; // IHADA6 = 0
1479 
1480  if ( theInitType == DEFAULT || theInitType == CORSIKA ) {
1481  popcor_.pdb = 0.1; // PDB = 0.10D0
1482  } else if ( theInitType == DPM2_5 ) {
1483  popcor_.pdb = 0.1; // PDB = 0.10D0
1484  } else if ( theInitType == DPM3 ) {
1485  popcor_.pdb = 0.15;
1486  }
1487  popcor_.ajsdef = 0.0; // AJSDEF = 0.D0
1488 //
1489 //
1490 // Definition of fluctuat, intpt, hadroniz, diquarks, singlech, evapor
1491 // (Charmed particles set to decay : IHADRINZ>=2
1492 //
1493  fluctu_.ifluct = 0; // IFLUCT = 0
1494  droppt_.intpt = LTRUE; // INTPT = .TRUE.
1495  colle_.ihadrz = 2; // IHADRZ = 2
1496  ifragm_.ifrag = 1; // IFRAG = 1
1497  promu_.ipromu = 1; // IPROMU = 1
1498  if (colle_.ihadrz >= 2)
1499  {
1500  ifragm_.ifrag = colle_.ihadrz - 1;
1501  // IFRAG = IHADRZ-1
1502  lundin_ (); // CALL LUNDIN
1503  }
1504  diquax_.idiqua = 1; // IDIQUA = 1
1505  diquax_.idiquu = 1; // IDIQUU = 1
1506  diquax_.amedd = 0.9; // AMEDD = 0.9D0
1507  sincha_.isicha = 0; // ISICHA = 0
1508  evappp_.ievap = 0; // IEVAP = 0
1509  seaqxx_.seaqx = 0.5; // SEAQX = 0.5D0
1510  seaqxx_.seaqxn = 0.5; // SEAQXN = 0.5D0
1511  kglaub_.jglaub = 2; // JGLAUB = 2
1512  hadthr_.ehadth = 5.0; // EHADTH = 5.D0
1513 //
1514 //
1515 // Definitions for hbook, pomtable, cmhisto, central, strucfun
1516 //
1517  hboo_.ihbook = 1; // IHBOOK = 1
1518  pomtab_.ipomta = 1; // IPOMTA = 1
1519  cmhico_. cmhis = 0.0; // CMHIS = 0.0D+00 ! Lab System
1520  zentra_.icentr = 0; // ICENTR = 0
1521  user2_.istruf = 222; // ISTRUF = 222
1522  strufu_.istrum = 0; // ISTRUM = 0
1523  strufu_.istrut = user2_.istruf / 100;
1524  // ISTRUT = ISTRUF/100
1525  user2_.istruf = user2_.istruf - strufu_.istrut*100;
1526  // ISTRUF = ISTRUF-ISTRUT*100
1527  strufu_.istrum = user2_.istruf; // ISTRUM = ISTRUF
1528 
1529  if ( theInitType == DEFAULT ) {
1530  ISINGD = 1;
1531  ISINGX = 1;
1532  IDUBLD = 0;
1533  SDFRAC = 1.0;
1534  } else if ( theInitType == CORSIKA ) {
1535  ISINGD = 0; // ISINGD = 0
1536  ISINGX = 0; // ISINGX = 0
1537  IDUBLD = 0; // IDUBLD = 0
1538  SDFRAC = 0.0; // SDFRAC = 0.
1539  } else if ( theInitType == DPM2_5 ) {
1540  ISINGD = 1; // ISINGD = 1
1541  ISINGX = 1; // ISINGX = 1
1542  IDUBLD = 0; // IDUBLD = 0
1543  SDFRAC = 1.0; // SDFRAC = 1.
1544  } else if ( theInitType == DPM3 ) {
1545  ISINGD = 0;
1546  ISINGX = 1;
1547  IDUBLD = 0;
1548  SDFRAC = 1.0;
1549  }
1550  diffra_.isingd = ISINGD;
1551  user2_.isingx = ISINGX;
1552  user2_.idubld = IDUBLD;
1553  user2_.sdfrac = SDFRAC;
1554 //
1555 //
1556 // Definitions for start, inforeje, gluxplit, partev, sampt
1557 // NEVNTS passed to DTMAI as argument, NEVHAD in COMMON
1558 //
1559  colle_.nfile = 0; // NFILE = 0
1560  nstari_.nstart = 1; // NSTART = 1
1561  stars_.istar2 = 0; // ISTAR2 = 0
1562  stars_.istar3 = 0; // ISTAR3 = 0
1563  user2_.ptlar = 2.0; // PTLAR = 2.D0
1564 
1565 // G4int iglaub = 0; // IGLAUB = 0 !! Note that this is just a local variable
1566 
1567 // infore_.ifrej = 0; // IFREJ = 0 Note rejection diagnostics not required
1568  gluspl_.nugluu = 1; // NUGLUU = 1
1569  gluspl_.nsgluu = 0; // NSGLUU = 0
1570  colle_.nvers = 1; // NVERS = 1
1571  ptsamp_.isampt = 4; // ISAMPT = 4
1572 //
1573 //
1574 // Definitions for selhard, sigmapom, pshow, secinter.
1575 //
1576  dropjj_.dropjt = 0.0; // DROPJT = 0.D0
1577  collis_.iophrd = 2; // IOPHRD = 2
1578  collis_.ptthr = 3.0; // PTTHR = 3.D0
1579 // collis_.ptthr2 = collis_.ptthr; // PTTHR2 = PTTHR
1580  user2_.cmener = 100.0; // CMENER = 100.D0
1581  if (strufu_.istrut == 1)
1582  {
1583  collis_.ptthr = 2.1 + 0.15*std::pow(std::log10(user2_.cmener/50.),3.0);
1584  // PTTHR = 2.1D0+0.15D0*(LOG10(CMENER/50.))**3
1585  }
1586  else if (strufu_.istrut == 2)
1587  {
1588  collis_.ptthr = 2.5 + 0.12*std::pow(std::log10(user2_.cmener/50.),3.0);
1589  // PTTHR = 2.5D0+0.12D0*(LOG10(CMENER/50.))**3
1590  }
1591  collis_.ptthr2 = collis_.ptthr; // PTTHR2 = PTTHR
1592  pomtyp_.ipim = 2; // IPIM = 2
1593  pomtyp_.icon = 48; // ICON = 48
1594  pomtyp_.isig = 10; // ISIG = 10
1595  pomtyp_.lmax = 30; // LMAX = 30
1596  pomtyp_.mmax = 100; // MMAX = 100
1597  pomtyp_.nmax = 2; // NMAX = 2
1598  pomtyp_.difel = 0.0; // DIFEL = 0.D0
1599  pomtyp_.difnu = 1.0; // DIFNU = 1.D0
1600  pshow_.ipshow = 1; // IPSHOW = 1
1601  secint_.isecin = 0; // ISECIN = 0
1602 //
1603 //
1604 // This next bit is associated with evaporation. I'm not sure if it's needed
1605 // as IEVAP = 0, but will initialise in any case.
1606 //
1607  if ( theInitType == DEFAULT ) {
1608  parevt_.levprt = LTRUE;
1609  parevt_.ilvmod = 1;
1610  parevt_.ldeexg = LFALSE;
1611  parevt_.lheavy = LFALSE;
1612  frbkcm_.lfrmbk = LFALSE;
1613  inpflg_.ifiss = 0;
1614  } else if ( theInitType == CORSIKA ) {
1615  parevt_.levprt = LTRUE; // LEVPRT = .TRUE.
1616  parevt_.ilvmod = 1; // ILVMOD = 1
1617  parevt_.ldeexg = LTRUE; // LDEEXG = .TRUE.
1618  parevt_.lheavy = LTRUE; // LHEAVY = .TRUE.
1619  frbkcm_.lfrmbk = LTRUE; // LFRMBK = .TRUE.
1620  inpflg_.ifiss = 0; // IFISS = 0
1621  } else if ( theInitType == DPM2_5 ) {
1622  parevt_.levprt = LFALSE; // LEVPRT = .FALSE.
1623  parevt_.ilvmod = 1; // ILVMOD = 1
1624  parevt_.ldeexg = LFALSE; // LDEEXG = .FALSE.
1625  parevt_.lheavy = LFALSE; // LHEAVY = .FALSE.
1626  frbkcm_.lfrmbk = LFALSE; // LFRMBK = .FALSE.
1627  inpflg_.ifiss = 0; // IFISS = 0
1628  } else if ( theInitType == DPM3 ) {
1629  parevt_.levprt = LFALSE;
1630  parevt_.ilvmod = 1;
1631  parevt_.ldeexg = LFALSE;
1632  parevt_.lheavy = LFALSE;
1633  frbkcm_.lfrmbk = LFALSE;
1634  inpflg_.ifiss = 0;
1635  }
1636 
1637  hettp_.nbertp = lunber; // NBERTP = LUNBER
1638 
1639  verboseFortranFile = "fort.6";
1640  G4int namelen = verboseFortranFile.length();
1641  char *ptr1 = new char[namelen+1];
1642  verboseFortranFile.copy(ptr1,namelen,0);
1643  ptr1[namelen] = '\0';
1644 // ptr1 = const_cast<char*> (verboseFortranFile.c_str());
1645  ftnlogical opened = LFALSE;
1646  g4dpmjet_open_fort6_ (&namelen, &opened, ptr1);
1647  delete [] ptr1;
1648  if (opened == LFALSE)
1649  {
1650  G4cout <<"ATTEMPTED TO OPEN fort.6 TO OUTPUT VERBOSE FORTRAN TEXT" <<G4endl;
1651  G4cout <<"HOWEVER THIS WAS NOT POSSIBLE" <<G4endl;
1652  }
1653 
1654 #ifdef G4VERBOSE
1655  if (GetVerboseLevel()>0) {
1656  G4cout <<"AT G4DPMJET2_5Model::Initialise: before NUCLEAR.BIN" <<G4endl;
1657  G4cout <<"OPENING NUCLEAR.BIN ON FILE UNIT " <<lunber <<G4endl;
1658  }
1659 #endif
1660  if ( !getenv("G4DPMJET2_5DATA") )
1661  {
1662  G4cout <<"ENVIRONMENT VARIABLE G4DPMJET2_5DATA NOT SET " <<G4endl;
1663  throw G4HadronicException(__FILE__, __LINE__,
1664  "Please setenv G4DPMJET2_5DATA to point to the dpmjet2.5 data files.");
1665  }
1666  defaultDirName = G4String(getenv("G4DPMJET2_5DATA")) + "/NUCLEAR.BIN";
1667  namelen = defaultDirName.length();
1668  ptr1 = new char[namelen+1];
1669  defaultDirName.copy(ptr1,namelen,0);
1670  ptr1[namelen] = '\0';
1671 // ptr1 = const_cast<char*> (defaultDirName.c_str());
1672  opened = LFALSE;
1673  g4dpmjet_open_nuclear_bin_ (&namelen, &lunber, &opened, ptr1);
1674  delete [] ptr1;
1675  if (opened == LFALSE)
1676  {
1677 //
1678 //
1679 // Problems with locating NUCLEAR.BIN file.
1680 //
1681  G4cout <<"NUCLEAR.BIN FILE NOT FOUND IN DIRECTORY " <<defaultDirName
1682  <<G4endl;
1683  throw G4HadronicException(__FILE__, __LINE__,
1684  "NUCLEAR.BIN file not present.");
1685  }
1686 #ifdef G4VERBOSE
1687  if (GetVerboseLevel()>0)
1688  G4cout <<"AT G4DPMJET2_5Model::Initialise: after NUCLEAR.BIN" <<G4endl;
1689 #endif
1690  G4cout << "CALL BERTTP" << G4endl;
1691  berttp_ (); // CALL BERTTP
1692  G4cout << "CALL BERTTP done" << G4endl;
1693  if (evappp_.ievap == 1)
1694  {
1695 #ifdef G4VERBOSE
1696  if (GetVerboseLevel()>0)
1697  G4cout <<"AT G4DPMJET2_5Model::Initialise: before INCINI" <<G4endl;
1698 #endif
1699  incini_ (); // CALL INCINI
1700 #ifdef G4VERBOSE
1701  if (GetVerboseLevel()>0)
1702  G4cout <<"AT G4DPMJET2_5Model::Initialise: after INCINI" <<G4endl;
1703 #endif
1704  }
1705  g4dpmjet_close_nuclear_bin_ (&lunber);
1706 #ifdef G4VERBOSE
1707  if (GetVerboseLevel()>0)
1708  G4cout <<"AT G4DPMJET2_5Model::Initialise: NUCLEAR.BIN closed" <<G4endl;
1709 #endif
1710 
1711  ptlarg_.xsmax = 0.8; // XSMAX = 0.8D0
1712  dprin_.itopd = 0; // ITOPD = 0
1713 
1714  G4int iseed1 = 0;
1715  G4int iseed2 = 0;
1716  rd2out_ (&iseed1,&iseed2);
1717 //
1718 //
1719 // This next bit outputs important variables if debug is switched on.
1720 //
1721 #ifdef G4VERBOSE
1722  if (GetVerboseLevel()>0) {
1723  G4cout <<"AT G4DPMJET2_5Model::Initialise:" <<G4endl;
1724  G4cout <<"Printout of important Parameters before DPMJET2.5" <<G4endl;
1725  G4cout <<"Please note for DPMJET input all numbers are floating point!" <<G4endl;
1726  G4cout <<"PROJPAR " <<nucc_.ip <<" "
1727  <<nucc_.ipz <<G4endl;
1728  G4cout <<"TARPAR " <<nucc_.it <<" "
1729  <<nucc_.itz <<G4endl;
1730  G4cout <<"MOMENTUM " <<ppn <<G4endl;
1731  G4cout <<"ENERGY " <<epn <<G4endl;
1732  G4cout <<"CMENERGY " <<nncms_.umo <<G4endl;
1733  G4cout <<"NOFINALE " <<final_.ifinal <<G4endl;
1734  G4cout <<"EVAPORAT " <<evappp_.ievap <<G4endl;
1735  G4cout <<"OUTLEVEL " <<dprin_.ipri <<" "
1736  <<dprin_.ipev <<" "
1737  <<dprin_.ippa <<" "
1738  <<dprin_.ipco <<" "
1739  <<dprin_.init <<" "
1740  <<dprin_.iphkk <<G4endl;
1741  G4cout <<"RANDOMIZ " <<iseed1 <<" "
1742  <<iseed2 <<G4endl;
1743  G4cout <<"STRUCFUN " <<user2_.istruf+100*strufu_.istrut <<G4endl;
1744  G4cout <<"SAMPT " <<ptsamp_.isampt <<G4endl;
1745  G4cout <<"SELHARD " <<0 <<" "
1746  <<collis_.iophrd <<" "
1747  <<0 <<" "
1748  <<dropjj_.dropjt <<" "
1749  <<collis_.ptthr <<" "
1750  <<collis_.ptthr2 << G4endl;
1751  G4cout <<"SIGMAPOM " <<0 <<" "
1752  <<pomtyp_.isig <<" "
1753  <<pomtyp_.ipim + 10*pomtyp_.icon <<" "
1754  <<pomtyp_.lmax <<" "
1755  <<pomtyp_.mmax <<" "
1756  <<pomtyp_.nmax <<G4endl;
1757  G4cout <<"PSHOWER " <<pshow_.ipshow <<G4endl;
1758  G4cout <<"CENTRAL " <<zentra_.icentr <<G4endl;
1759  G4cout <<"CMHISTO " <<cmhico_.cmhis <<G4endl;
1760  G4cout <<"SEASU3 " <<seasu3_.seasq <<G4endl;
1761  G4cout <<"RECOMBIN " <<recom_.irecom <<G4endl;
1762  G4cout <<"SINGDIFF " <<diffra_.isingd <<G4endl;
1763  G4cout <<"TAUFOR " <<taufo_.taufor <<" "
1764  <<taufo_.ktauge <<" "
1765  <<taufo_.itauve <<G4endl;
1766  G4cout <<"POPCORN " <<popcor_.pdb <<G4endl;
1767  G4cout <<"POPCORCK " <<popcck_.ijpock <<" "
1768  <<popcck_.pdbck <<G4endl;
1769  G4cout <<"POPCORSE " <<popcck_.pdbse <<" "
1770  <<popcck_.pdbseu <<G4endl;
1771  G4cout <<"CASADIQU " <<casadi_.icasad <<" "
1772  <<casadi_.casaxx <<G4endl;
1773  G4cout <<"DIQUARKS " <<diquax_.idiqua <<" "
1774  <<diquax_.idiquu <<" "
1775  <<diquax_.amedd <<G4endl;
1776  G4cout <<"HADRONIZ " <<colle_.ihadrz <<G4endl;
1777  G4cout <<"INTPT " <<droppt_.intpt <<G4endl;
1778  G4cout <<"PAULI " <<droppt_.lpauli <<G4endl;
1779  G4cout <<"FERMI " <<droppt_.fermp <<" "
1780  <<nucimp_.fermod <<G4endl;
1781  G4cout <<"CRONINPT " <<cronin_.mkcron <<" "
1782  <<cronin_.cronco <<G4endl;
1783  G4cout <<"SEADISTR " <<xseadi_.xseacu+0.95 <<" "
1784  <<xseadi_.unon <<" "
1785  <<xseadi_.unom <<" "
1786  <<xseadi_.unosea <<G4endl;
1787  G4cout <<"SEAQUARK " <<seaqxx_.seaqx <<" "
1788  <<seaqxx_.seaqxn <<G4endl;
1789  G4cout <<"SECINTER " <<secint_.isecin <<G4endl;
1790  G4cout <<"XCUTS " <<xseadi_.cvq <<" "
1791  <<xseadi_.cdq <<" "
1792  <<xseadi_.csea <<" "
1793  <<xseadi_.ssmima <<G4endl;
1794  }
1795 #endif
1796 
1797  bufues_.bnnvv = 0.001; // BNNVV=0.001
1798  bufues_.bnnss = 0.001; // BNNSS=0.001
1799  bufues_.bnnsv = 0.001; // BNNSV=0.001
1800  bufues_.bnnvs = 0.001; // BNNVS=0.001
1801  bufues_.bnncc = 0.001; // BNNCC=0.001
1802  bufues_.bnndv = 0.001; // BNNDV=0.001
1803  bufues_.bnnvd = 0.001; // BNNVD=0.001
1804  bufues_.bnnds = 0.001; // BNNDS=0.001
1805  bufues_.bnnsd = 0.001; // BNNSD=0.001
1806  bufues_.bnnhh = 0.001; // BNNHH=0.001
1807  bufues_.bnnzz = 0.001; // BNNZZ=0.001
1808  bufues_.bnndi = 0.001; // BNNDI=0.001
1809  bufues_.bnnzd = 0.001; // BNNZD=0.001
1810  bufues_.bnndz = 0.001; // BNNDZ=0.001
1811  bufues_.bptvv = 0.0; // BPTVV=0.
1812  bufues_.bptss = 0.0; // BPTSS=0.
1813  bufues_.bptsv = 0.0; // BPTSV=0.
1814  bufues_.bptvs = 0.0; // BPTVS=0.
1815  bufues_.bptcc = 0.0; // BPTCC=0.
1816  bufues_.bptdv = 0.0; // BPTDV=0.
1817  bufues_.bptvd = 0.0; // BPTVD=0.
1818  bufues_.bptds = 0.0; // BPTDS=0.
1819  bufues_.bptsd = 0.0; // BPTSD=0.
1820  bufues_.bpthh = 0.0; // BPTHH=0.
1821  bufues_.bptzz = 0.0; // BPTZZ=0.
1822  bufues_.bptdi = 0.0; // BPTDI=0.
1823  bufues_.bptzd = 0.0; // BPTZD=0.
1824  bufues_.bptdz = 0.0; // BPTDZ=0.
1825  bufues_.beevv = 0.0; // BEEVV=0.
1826  bufues_.beess = 0.0; // BEESS=0.
1827  bufues_.beesv = 0.0; // BEESV=0.
1828  bufues_.beevs = 0.0; // BEEVS=0.
1829  bufues_.beecc = 0.0; // BEECC=0.
1830  bufues_.beedv = 0.0; // BEEDV=0.
1831  bufues_.beevd = 0.0; // BEEVD=0.
1832  bufues_.beeds = 0.0; // BEEDS=0.
1833  bufues_.beesd = 0.0; // BEESD=0.
1834  bufues_.beehh = 0.0; // BEEHH=0.
1835  bufues_.beezz = 0.0; // BEEZZ=0.
1836  bufues_.beedi = 0.0; // BEEDI=0.
1837  bufues_.beezd = 0.0; // BEEZD=0.
1838  bufues_.beedz = 0.0; // BEEDZ=0.
1839  ncoucs_.bcouvv = 0.0; // BCOUVV=0.
1840  ncoucs_.bcouss = 0.0; // BCOUSS=0.
1841  ncoucs_.bcousv = 0.0; // BCOUSV=0.
1842  ncoucs_.bcouvs = 0.0; // BCOUVS=0.
1843  ncoucs_.bcouzz = 0.0; // BCOUZZ=0.
1844  ncoucs_.bcouhh = 0.0; // BCOUHH=0.
1845  ncoucs_.bcouds = 0.0; // BCOUDS=0.
1846  ncoucs_.bcousd = 0.0; // BCOUSD=0.
1847  ncoucs_.bcoudz = 0.0; // BCOUDZ=0.
1848  ncoucs_.bcouzd = 0.0; // BCOUZD=0.
1849  ncoucs_.bcoudi = 0.0; // BCOUDI=0.
1850  ncoucs_.bcoudv = 0.0; // BCOUDV=0.
1851  ncoucs_.bcouvd = 0.0; // BCOUVD=0.
1852  ncoucs_.bcoucc = 0.0; // BCOUCC=0.
1853 //
1854 //
1855 // Initialisation of
1856 // ANNVV, ANNSS ... ANNDZ
1857 // PTVV, PTSS ... PTDZ
1858 // EEVV, EESS ... EEDZ
1859 // ACOUVV, ACOUSS, ... ACOUCC
1860 // now all moved to ApplyYourself member function.
1861 //
1862 // droppt_.ipadis = LFALSE; // IPADIS = .FALSE.
1863 // droppt_.ihadvv = LFALSE; // IHADVV = .FALSE.
1864 // droppt_.ihadsv = LFALSE; // IHADSV = .FALSE.
1865 // droppt_.ihadvs = LFALSE; // IHADVS = .FALSE.
1866 
1867  nucc_.ijtarg = 1; // IJTARG=1
1868 
1869 //
1870 //
1871 // The following commented out since it seems to have more to do with histogram
1872 // generation.
1873 //
1874 // i = 1;
1875 // G4int idummy;
1876 // distr_ (&i,&nucc_.ijproj,&ppn,&idummy); // CALL DISTR( 1,IJPROJ,PPN,IDUMMY )
1877  if (pomtyp_.ipim == 2) {prblm2_ (&user2_.cmener);}
1878 //
1879 //
1880 // Initialise hard scattering & transverse momentum for soft scattering
1881 //
1882  i = 0;
1883  jtdtu_ (&i); // CALL JTDTU( 0 )
1884  i = 0;
1885  G4double pt;
1886  samppt_ (&i,&pt); // CALL SAMPPT(0,PT)
1887 }
1889 //
1890 #endif
Generated on Sat May 25 2013 14:32:57 for Geant4 by   doxygen 1.8.4