Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RPGFragmentation.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 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id$
27 //
28 
29 #include <iostream>
30 #include <signal.h>
31 
32 #include "G4RPGFragmentation.hh"
33 #include "G4PhysicalConstants.hh"
34 #include "G4SystemOfUnits.hh"
35 #include "G4AntiProton.hh"
36 #include "G4AntiNeutron.hh"
37 #include "Randomize.hh"
38 #include "G4Poisson.hh"
40 
42  : G4RPGReaction()
43 {
44  for (G4int i = 0; i < 20; i++) dndl[i] = 0.0;
45 }
46 
47 
50 {
51  pt = std::max( 0.001, pt );
52  G4double dx = 1./(19.*pt);
53  G4double x;
54  G4double term1;
55  G4double term2;
56 
57  for (G4int i = 1; i < 20; i++) {
58  x = (G4double(i) - 0.5)*dx;
59  term1 = 1. + parMass*parMass*x*x;
60  term2 = pt*x*et*pt*x*et + pt*pt + secMass*secMass;
61  dndl[i] = dx / std::sqrt( term1*term1*term1*term2 )
62  + dndl[i-1];
63  }
64 }
65 
66 
68 ReactionStage(const G4HadProjectile* originalIncident,
69  G4ReactionProduct& modifiedOriginal,
70  G4bool& incidentHasChanged,
71  const G4DynamicParticle* originalTarget,
72  G4ReactionProduct& targetParticle,
73  G4bool& targetHasChanged,
74  const G4Nucleus& targetNucleus,
75  G4ReactionProduct& currentParticle,
77  G4int& vecLen,
78  G4bool leadFlag,
79  G4ReactionProduct& leadingStrangeParticle)
80 {
81  //
82  // Based on H. Fesefeldt's original FORTRAN code GENXPT
83  //
84  // Generation of x- and pT- values for incident, target, and all secondary
85  // particles using a simple single variable description E D3S/DP3= F(Q)
86  // with Q^2 = (M*X)^2 + PT^2. Final state kinematics are produced by an
87  // FF-type iterative cascade method
88  //
89  // Internal units are GeV
90  //
91 
92  // Protection in case no secondary has been created. In that case use
93  // two-body scattering
94  //
95  if (vecLen == 0) return false;
96 
102 
103  G4int i, l;
104  G4bool veryForward = false;
105 
106  const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/GeV;
107  const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
108  const G4double mOriginal = modifiedOriginal.GetMass()/GeV;
109  const G4double pOriginal = modifiedOriginal.GetMomentum().mag()/GeV;
110  G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/GeV;
111  G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
112  targetMass*targetMass +
113  2.0*targetMass*etOriginal ); // GeV
114  G4double currentMass = currentParticle.GetMass()/GeV;
115  targetMass = targetParticle.GetMass()/GeV;
116 
117  // Randomize the order of the secondary particles.
118  // Note that the current and target particles are not affected.
119 
120  for (i=0; i<vecLen; ++i) {
121  G4int itemp = G4int( G4UniformRand()*vecLen );
122  G4ReactionProduct pTemp = *vec[itemp];
123  *vec[itemp] = *vec[i];
124  *vec[i] = pTemp;
125  }
126 
127  if (currentMass == 0.0 && targetMass == 0.0) {
128  // Target and projectile have annihilated. Replace them with the first
129  // two secondaries in the list. Current particle KE is maintained.
130 
131  G4double ek = currentParticle.GetKineticEnergy();
132  G4ThreeVector mom = currentParticle.GetMomentum();
133  currentParticle = *vec[0];
134  currentParticle.SetSide(1);
135  targetParticle = *vec[1];
136  targetParticle.SetSide(-1);
137  for( i=0; i<(vecLen-2); ++i )*vec[i] = *vec[i+2];
138  G4ReactionProduct *temp = vec[vecLen-1];
139  delete temp;
140  temp = vec[vecLen-2];
141  delete temp;
142  vecLen -= 2;
143  currentMass = currentParticle.GetMass()/GeV;
144  targetMass = targetParticle.GetMass()/GeV;
145  incidentHasChanged = true;
146  targetHasChanged = true;
147  currentParticle.SetKineticEnergy( ek );
148  currentParticle.SetMomentum(mom);
149  veryForward = true;
150  }
151  const G4double atomicWeight = targetNucleus.GetA_asInt();
152  const G4double atomicNumber = targetNucleus.GetZ_asInt();
153  const G4double protonMass = aProton->GetPDGMass();
154 
155  if (originalIncident->GetDefinition()->GetParticleSubType() == "kaon"
156  && G4UniformRand() >= 0.7) {
157  G4ReactionProduct temp = currentParticle;
158  currentParticle = targetParticle;
159  targetParticle = temp;
160  incidentHasChanged = true;
161  targetHasChanged = true;
162  currentMass = currentParticle.GetMass()/GeV;
163  targetMass = targetParticle.GetMass()/GeV;
164  }
165  const G4double afc = std::min( 0.75,
166  0.312+0.200*std::log(std::log(centerofmassEnergy*centerofmassEnergy))+
167  std::pow(centerofmassEnergy*centerofmassEnergy,1.5)/6000.0 );
168 
169  G4double freeEnergy = centerofmassEnergy-currentMass-targetMass;
170  G4double forwardEnergy = freeEnergy/2.;
171  G4int forwardCount = 1; // number of particles in forward hemisphere
172 
173  G4double backwardEnergy = freeEnergy/2.;
174  G4int backwardCount = 1; // number of particles in backward hemisphere
175 
176  if(veryForward) {
177  if(currentParticle.GetSide()==-1)
178  {
179  forwardEnergy += currentMass;
180  forwardCount --;
181  backwardEnergy -= currentMass;
182  backwardCount ++;
183  }
184  if(targetParticle.GetSide()!=-1)
185  {
186  backwardEnergy += targetMass;
187  backwardCount --;
188  forwardEnergy -= targetMass;
189  forwardCount ++;
190  }
191  }
192 
193  for (i=0; i<vecLen; ++i) {
194  if( vec[i]->GetSide() == -1 )
195  {
196  ++backwardCount;
197  backwardEnergy -= vec[i]->GetMass()/GeV;
198  } else {
199  ++forwardCount;
200  forwardEnergy -= vec[i]->GetMass()/GeV;
201  }
202  }
203 
204  // Check that sum of forward particle masses does not exceed forwardEnergy,
205  // and similarly for backward particles. If so, move particles from one
206  // hemisphere to another.
207 
208  if (backwardEnergy < 0.0) {
209  for (i = 0; i < vecLen; ++i) {
210  if (vec[i]->GetSide() == -1) {
211  backwardEnergy += vec[i]->GetMass()/GeV;
212  --backwardCount;
213  vec[i]->SetSide(1);
214  forwardEnergy -= vec[i]->GetMass()/GeV;
215  ++forwardCount;
216  if (backwardEnergy > 0.0) break;
217  }
218  }
219  }
220 
221  if (forwardEnergy < 0.0) {
222  for (i = 0; i < vecLen; ++i) {
223  if (vec[i]->GetSide() == 1) {
224  forwardEnergy += vec[i]->GetMass()/GeV;
225  --forwardCount;
226  vec[i]->SetSide(-1);
227  backwardEnergy -= vec[i]->GetMass()/GeV;
228  ++backwardCount;
229  if (forwardEnergy > 0.0) break;
230  }
231  }
232  }
233 
234  // Special cases for reactions near threshold
235 
236  // 1. There is only one secondary
237  if (forwardEnergy > 0.0 && backwardEnergy < 0.0) {
238  forwardEnergy += backwardEnergy;
239  backwardEnergy = 0;
240  }
241 
242  // 2. Nuclear binding energy is large
243  if (forwardEnergy + backwardEnergy < 0.0) return false;
244 
245 
246  // forwardEnergy now represents the total energy in the forward reaction
247  // hemisphere which is available for kinetic energy and particle creation.
248  // Similarly for backwardEnergy.
249 
250  // Add particles from the intranuclear cascade.
251  // nuclearExcitationCount = number of new secondaries produced by nuclear
252  // excitation
253  // extraCount = number of nucleons within these new secondaries
254  //
255  // Note: eventually have to make sure that enough nucleons are available
256  // in the case of small target nuclei
257 
258  G4double xtarg;
259  if( centerofmassEnergy < (2.0+G4UniformRand()) )
260  xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2.0*backwardCount+vecLen+2)/2.0;
261  else
262  xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2.0*backwardCount);
263  if( xtarg <= 0.0 )xtarg = 0.01;
264  G4int nuclearExcitationCount = G4Poisson( xtarg );
265  // To do: try reducing nuclearExcitationCount with increasing energy
266  // to simulate cut-off of cascade
267  if(atomicWeight<1.0001) nuclearExcitationCount = 0;
268  G4int extraNucleonCount = 0;
269  G4double extraNucleonMass = 0.0;
270 
271  if (nuclearExcitationCount > 0) {
272  const G4double nucsup[] = { 1.00, 0.7, 0.5, 0.4, 0.35, 0.3 };
273  const G4double psup[] = { 3., 6., 20., 50., 100., 1000. };
274  G4int momentumBin = 0;
275  while( (momentumBin < 6) &&
276  (modifiedOriginal.GetTotalMomentum()/GeV > psup[momentumBin]) )
277  ++momentumBin;
278  momentumBin = std::min( 5, momentumBin );
279 
280  // NOTE: in GENXPT, these new particles were given negative codes
281  // here I use NewlyAdded = true instead
282  //
283 
284  for (i = 0; i < nuclearExcitationCount; ++i) {
285  G4ReactionProduct * pVec = new G4ReactionProduct();
286  if (G4UniformRand() < nucsup[momentumBin]) {
287 
288  if (G4UniformRand() > 1.0-atomicNumber/atomicWeight)
289  pVec->SetDefinition( aProton );
290  else
291  pVec->SetDefinition( aNeutron );
292 
293  // nucleon comes from nucleus -
294  // do not subtract its mass from backward energy
295  pVec->SetSide( -2 ); // -2 means backside nucleon
296  ++extraNucleonCount;
297  extraNucleonMass += pVec->GetMass()/GeV;
298  // To do: Remove chosen nucleon from target nucleus
299  pVec->SetNewlyAdded( true );
300  vec.SetElement( vecLen++, pVec );
301  ++backwardCount;
302 
303  } else {
304 
305  G4double ran = G4UniformRand();
306  if( ran < 0.3181 )
307  pVec->SetDefinition( aPiPlus );
308  else if( ran < 0.6819 )
309  pVec->SetDefinition( aPiZero );
310  else
311  pVec->SetDefinition( aPiMinus );
312 
313  if (backwardEnergy > pVec->GetMass()/GeV) {
314  backwardEnergy -= pVec->GetMass()/GeV; // pion mass comes from free energy
315  ++backwardCount;
316  pVec->SetSide( -1 ); // backside particle, but not a nucleon
317  pVec->SetNewlyAdded( true );
318  vec.SetElement( vecLen++, pVec );
319  }
320 
321  // To do: Change proton to neutron (or vice versa) in target nucleus depending
322  // on pion charge
323  }
324  }
325  }
326 
327  // Define initial state vectors for Lorentz transformations
328  // The pseudoParticles have non-standard masses, hence the "pseudo"
329 
330  G4ReactionProduct pseudoParticle[8];
331  for (i = 0; i < 8; ++i) pseudoParticle[i].SetZero();
332 
333  pseudoParticle[0].SetMass( mOriginal*GeV );
334  pseudoParticle[0].SetMomentum( 0.0, 0.0, pOriginal*GeV );
335  pseudoParticle[0].SetTotalEnergy(
336  std::sqrt( pOriginal*pOriginal + mOriginal*mOriginal )*GeV );
337 
338  pseudoParticle[1].SetMass(protonMass); // this could be targetMass
339  pseudoParticle[1].SetTotalEnergy(protonMass);
340 
341  pseudoParticle[3].SetMass(protonMass*(1+extraNucleonCount) );
342  pseudoParticle[3].SetTotalEnergy(protonMass*(1+extraNucleonCount) );
343 
344  pseudoParticle[2] = pseudoParticle[0] + pseudoParticle[1];
345  pseudoParticle[3] = pseudoParticle[3] + pseudoParticle[0];
346 
347  pseudoParticle[0].Lorentz( pseudoParticle[0], pseudoParticle[2] );
348  pseudoParticle[1].Lorentz( pseudoParticle[1], pseudoParticle[2] );
349 
350  // Main loop for 4-momentum generation
351  // See Pitha-report (Aachen) for a detailed description of the method
352 
353  G4double aspar, pt, et, x, pp, pp1, wgt;
354  G4int innerCounter, outerCounter;
355  G4bool eliminateThisParticle, resetEnergies, constantCrossSection;
356 
357  G4double forwardKinetic = 0.0;
358  G4double backwardKinetic = 0.0;
359 
360  // Process the secondary particles in reverse order
361  // The incident particle is done after the secondaries
362  // Nucleons, including the target, in the backward hemisphere are also
363  // done later
364 
365  G4int backwardNucleonCount = 0; // number of nucleons in backward hemisphere
366  G4double totalEnergy, kineticEnergy, vecMass;
367  G4double phi;
368 
369  for (i = vecLen-1; i >= 0; --i) {
370 
371  if (vec[i]->GetNewlyAdded()) { // added from intranuclear cascade
372  if (vec[i]->GetSide() == -2) { // its a nucleon
373  if (backwardNucleonCount < 18) {
374  if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
375  for (G4int j = 0; j < vecLen; j++) delete vec[j];
376  vecLen = 0;
377  throw G4HadReentrentException(__FILE__, __LINE__,
378  "G4RPGFragmentation::ReactionStage : a pion has been counted as a backward nucleon");
379  }
380  vec[i]->SetSide(-3);
381  ++backwardNucleonCount;
382  continue; // Don't generate momenta for the first 17 backward
383  // cascade nucleons. This gets done by the cluster
384  // model later on.
385  }
386  }
387  }
388 
389  // Set pt and phi values, they are changed somewhat in the iteration loop
390  // Set mass parameter for lambda fragmentation model
391 
392  vecMass = vec[i]->GetMass()/GeV;
393  G4double ran = -std::log(1.0-G4UniformRand())/3.5;
394 
395  if (vec[i]->GetSide() == -2) { // backward nucleon
396  aspar = 0.20;
397  pt = std::sqrt( std::pow( ran, 1.2 ) );
398 
399  } else { // not a backward nucleon
400  if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
401  aspar = 0.75;
402  pt = std::sqrt( std::pow( ran, 1.7 ) );
403  } else if (vec[i]->GetDefinition()->GetParticleSubType() == "kaon") {
404  aspar = 0.70;
405  pt = std::sqrt( std::pow( ran, 1.7 ) );
406  } else { // vec[i] must be a baryon or ion
407  aspar = 0.65;
408  pt = std::sqrt( std::pow( ran, 1.5 ) );
409  }
410  }
411 
412  pt = std::max( 0.001, pt );
413  phi = G4UniformRand()*twopi;
414  vec[i]->SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
415  if (vec[i]->GetSide() > 0)
416  et = pseudoParticle[0].GetTotalEnergy()/GeV;
417  else
418  et = pseudoParticle[1].GetTotalEnergy()/GeV;
419 
420  //
421  // Start of outer iteration loop
422  //
423  outerCounter = 0;
424  eliminateThisParticle = true;
425  resetEnergies = true;
426  dndl[0] = 0.0;
427 
428  while (++outerCounter < 3) {
429 
430  FragmentationIntegral(pt, et, aspar, vecMass);
431  innerCounter = 0;
432  vec[i]->SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
433 
434  // Start of inner iteration loop
435 
436  while (++innerCounter < 7) {
437 
438  ran = G4UniformRand()*dndl[19];
439  l = 1;
440  while( ( ran > dndl[l] ) && ( l < 19 ) ) l++;
441  x = (G4double(l-1) + G4UniformRand())/19.;
442  if (vec[i]->GetSide() < 0) x *= -1.;
443  vec[i]->SetMomentum( x*et*GeV ); // set the z-momentum
444  totalEnergy = std::sqrt( x*et*x*et + pt*pt + vecMass*vecMass );
445  vec[i]->SetTotalEnergy( totalEnergy*GeV );
446  kineticEnergy = vec[i]->GetKineticEnergy()/GeV;
447 
448  if (vec[i]->GetSide() > 0) { // forward side
449  if( (forwardKinetic+kineticEnergy) < 0.95*forwardEnergy ) {
450  // Leave at least 5% of the forward free energy for the projectile
451 
452  pseudoParticle[4] = pseudoParticle[4] + (*vec[i]);
453  forwardKinetic += kineticEnergy;
454  outerCounter = 2; // leave outer loop
455  eliminateThisParticle = false; // don't eliminate this particle
456  resetEnergies = false;
457  break; // leave inner loop
458  }
459  if( innerCounter > 5 )break; // leave inner loop
460  if( backwardEnergy >= vecMass ) // switch sides
461  {
462  vec[i]->SetSide(-1);
463  forwardEnergy += vecMass;
464  backwardEnergy -= vecMass;
465  ++backwardCount;
466  }
467  } else { // backward side
468  // if (extraNucleonCount > 19) x = 0.999;
469  // G4double xxx = 0.95+0.05*extraNucleonCount/20.0;
470  // DHW: I think above lines were meant to be as follows:
471  G4double xxx = 0.999;
472  if (extraNucleonCount < 20) xxx = 0.95+0.05*extraNucleonCount/20.0;
473 
474  if ((backwardKinetic+kineticEnergy) < xxx*backwardEnergy) {
475  pseudoParticle[5] = pseudoParticle[5] + (*vec[i]);
476  backwardKinetic += kineticEnergy;
477  outerCounter = 2; // leave outer loop
478  eliminateThisParticle = false; // don't eliminate this particle
479  resetEnergies = false;
480  break; // leave inner loop
481  }
482  if (innerCounter > 5) break; // leave inner loop
483  if (forwardEnergy >= vecMass) { // switch sides
484  vec[i]->SetSide(1);
485  forwardEnergy -= vecMass;
486  backwardEnergy += vecMass;
487  backwardCount--;
488  }
489  }
490  G4ThreeVector momentum = vec[i]->GetMomentum();
491  vec[i]->SetMomentum( momentum.x() * 0.9, momentum.y() * 0.9 );
492  pt *= 0.9;
493  dndl[19] *= 0.9;
494  } // closes inner loop
495 
496  // If we get here, the inner loop has been done 6 times.
497  // If necessary, reduce energies of the previously done particles if
498  // they are lighter than protons or are in the forward hemisphere.
499  // Then continue with outer loop.
500 
501  if (resetEnergies)
502  ReduceEnergiesOfSecondaries(i+1, forwardKinetic, backwardKinetic,
503  vec, vecLen,
504  pseudoParticle[4], pseudoParticle[5],
505  pt);
506 
507  } // closes outer loop
508 
509  if (eliminateThisParticle && vec[i]->GetMayBeKilled()) {
510  // not enough energy, eliminate this particle
511 
512  if (vec[i]->GetSide() > 0) {
513  --forwardCount;
514  forwardEnergy += vecMass;
515  } else {
516  --backwardCount;
517  if (vec[i]->GetSide() == -2) {
518  --extraNucleonCount;
519  extraNucleonMass -= vecMass;
520  } else {
521  backwardEnergy += vecMass;
522  }
523  }
524 
525  for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up
526  G4ReactionProduct* temp = vec[vecLen-1];
527  delete temp;
528  // To do: modify target nucleus according to particle eliminated
529 
530  if( --vecLen == 0 ){
531  G4cout << " FALSE RETURN DUE TO ENERGY BALANCE " << G4endl;
532  return false;
533  } // all the secondaries have been eliminated
534  }
535  } // closes main loop over secondaries
536 
537  // Re-balance forward and backward energy if possible and if necessary
538 
539  G4double forwardKEDiff = forwardEnergy - forwardKinetic;
540  G4double backwardKEDiff = backwardEnergy - backwardKinetic;
541 
542  if (forwardKEDiff < 0.0 || backwardKEDiff < 0.0) {
543  ReduceEnergiesOfSecondaries(0, forwardKinetic, backwardKinetic,
544  vec, vecLen,
545  pseudoParticle[4], pseudoParticle[5],
546  pt);
547 
548  forwardKEDiff = forwardEnergy - forwardKinetic;
549  backwardKEDiff = backwardEnergy - backwardKinetic;
550  if (backwardKEDiff < 0.0) {
551  if (forwardKEDiff + backwardKEDiff > 0.0) {
552  backwardEnergy = backwardKinetic;
553  forwardEnergy += backwardKEDiff;
554  forwardKEDiff = forwardEnergy - forwardKinetic;
555  backwardKEDiff = 0.0;
556  } else {
557  G4cout << " False return due to insufficient backward energy " << G4endl;
558  return false;
559  }
560  }
561 
562  if (forwardKEDiff < 0.0) {
563  if (forwardKEDiff + backwardKEDiff > 0.0) {
564  forwardEnergy = forwardKinetic;
565  backwardEnergy += forwardKEDiff;
566  backwardKEDiff = backwardEnergy - backwardKinetic;
567  forwardKEDiff = 0.0;
568  } else {
569  G4cout << " False return due to insufficient forward energy " << G4endl;
570  return false;
571  }
572  }
573  }
574 
575  // Generate momentum for incident (current) particle, which was placed
576  // in the forward hemisphere.
577  // Set mass parameter for lambda fragmentation model.
578  // Set pt and phi values, which are changed somewhat in the iteration loop
579 
580  G4double ran = -std::log(1.0-G4UniformRand());
581  if (currentParticle.GetDefinition()->GetParticleSubType() == "pi") {
582  aspar = 0.60;
583  pt = std::sqrt( std::pow( ran/6.0, 1.7 ) );
584  } else if (currentParticle.GetDefinition()->GetParticleSubType() == "kaon") {
585  aspar = 0.50;
586  pt = std::sqrt( std::pow( ran/5.0, 1.4 ) );
587  } else {
588  aspar = 0.40;
589  pt = std::sqrt( std::pow( ran/4.0, 1.2 ) );
590  }
591 
592  phi = G4UniformRand()*twopi;
593  currentParticle.SetMomentum(pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV);
594  et = pseudoParticle[0].GetTotalEnergy()/GeV;
595  dndl[0] = 0.0;
596  vecMass = currentParticle.GetMass()/GeV;
597 
598  FragmentationIntegral(pt, et, aspar, vecMass);
599 
600  ran = G4UniformRand()*dndl[19];
601  l = 1;
602  while( ( ran > dndl[l] ) && ( l < 19 ) ) l++;
603  x = (G4double(l-1) + G4UniformRand())/19.;
604  currentParticle.SetMomentum( x*et*GeV ); // set the z-momentum
605 
606  if (forwardEnergy < forwardKinetic) {
607  totalEnergy = vecMass + 0.04*std::fabs(normal());
608  G4cout << " Not enough forward energy: forwardEnergy = "
609  << forwardEnergy << " forwardKinetic = "
610  << forwardKinetic << " total energy left = "
611  << backwardKEDiff + forwardKEDiff << G4endl;
612  } else {
613  totalEnergy = vecMass + forwardEnergy - forwardKinetic;
614  forwardKinetic = forwardEnergy;
615  }
616  currentParticle.SetTotalEnergy( totalEnergy*GeV );
617  pp = std::sqrt(std::abs( totalEnergy*totalEnergy - vecMass*vecMass) )*GeV;
618  pp1 = currentParticle.GetMomentum().mag();
619 
620  if (pp1 < 1.0e-6*GeV) {
621  G4ThreeVector iso = Isotropic(pp);
622  currentParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
623  } else {
624  currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
625  }
626  pseudoParticle[4] = pseudoParticle[4] + currentParticle;
627 
628  // Current particle now finished
629 
630  // Begin target particle
631 
632  if (backwardNucleonCount < 18) {
633  targetParticle.SetSide(-3);
634  ++backwardNucleonCount;
635 
636  } else {
637  // Set pt and phi values, they are changed somewhat in the iteration loop
638  // Set mass parameter for lambda fragmentation model
639 
640  vecMass = targetParticle.GetMass()/GeV;
641  ran = -std::log(1.0-G4UniformRand());
642  aspar = 0.40;
643  pt = std::max( 0.001, std::sqrt( std::pow( ran/4.0, 1.2 ) ) );
644  phi = G4UniformRand()*twopi;
645  targetParticle.SetMomentum(pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV);
646  et = pseudoParticle[1].GetTotalEnergy()/GeV;
647  outerCounter = 0;
648  innerCounter = 0;
649  G4bool marginalEnergy = true;
650  dndl[0] = 0.0;
651  G4double xxx = 0.999;
652  if( extraNucleonCount < 20 ) xxx = 0.95+0.05*extraNucleonCount/20.0;
653  G4ThreeVector momentum;
654 
655  while (++outerCounter < 4) {
656  FragmentationIntegral(pt, et, aspar, vecMass);
657 
658  for (innerCounter = 0; innerCounter < 6; innerCounter++) {
659  ran = G4UniformRand()*dndl[19];
660  l = 1;
661  while( ( ran > dndl[l] ) && ( l < 19 ) ) l++;
662  x = -(G4double(l-1) + G4UniformRand())/19.;
663  targetParticle.SetMomentum( x*et*GeV ); // set the z-momentum
664  totalEnergy = std::sqrt(x*et*x*et + pt*pt + vecMass*vecMass);
665  targetParticle.SetTotalEnergy( totalEnergy*GeV );
666 
667  if ((backwardKinetic+totalEnergy-vecMass) < xxx*backwardEnergy) {
668  pseudoParticle[5] = pseudoParticle[5] + targetParticle;
669  backwardKinetic += totalEnergy - vecMass;
670  outerCounter = 3; // leave outer loop
671  marginalEnergy = false;
672  break; // leave inner loop
673  }
674  momentum = targetParticle.GetMomentum();
675  targetParticle.SetMomentum(momentum.x() * 0.9, momentum.y() * 0.9);
676  pt *= 0.9;
677  dndl[19] *= 0.9;
678  }
679  }
680 
681  if (marginalEnergy) {
682  G4cout << " Extra backward kinetic energy = "
683  << 0.999*backwardEnergy - backwardKinetic << G4endl;
684  totalEnergy = vecMass + 0.999*backwardEnergy - backwardKinetic;
685  targetParticle.SetTotalEnergy(totalEnergy*GeV);
686  pp = std::sqrt(std::abs(totalEnergy*totalEnergy - vecMass*vecMass) )*GeV;
687  targetParticle.SetMomentum(momentum.x()/0.9, momentum.y()/0.9);
688  pp1 = targetParticle.GetMomentum().mag();
689  targetParticle.SetMomentum(targetParticle.GetMomentum() * pp/pp1 );
690  pseudoParticle[5] = pseudoParticle[5] + targetParticle;
691  backwardKinetic = 0.999*backwardEnergy;
692  }
693 
694  }
695 
696  if (backwardEnergy < backwardKinetic)
697  G4cout << " Backward Edif = " << backwardEnergy - backwardKinetic << G4endl;
698  if (forwardEnergy != forwardKinetic)
699  G4cout << " Forward Edif = " << forwardEnergy - forwardKinetic << G4endl;
700 
701  // Target particle finished.
702  // Now produce backward nucleons with a cluster model
703  // ps[2] = CM frame of projectile + target
704  // ps[3] = sum of projectile + nucleon cluster in lab frame
705  // ps[6] = proj + cluster 4-vector boosted into CM frame of proj + targ
706  // with secondaries, current and target particles subtracted
707  // = total 4-momentum of backward nucleon cluster
708 
709  pseudoParticle[6].Lorentz( pseudoParticle[3], pseudoParticle[2] );
710  pseudoParticle[6] = pseudoParticle[6] - pseudoParticle[4];
711  pseudoParticle[6] = pseudoParticle[6] - pseudoParticle[5];
712 
713  if (backwardNucleonCount == 1) {
714  // Target particle is the only backward nucleon. Give it the remainder
715  // of the backward kinetic energy.
716 
717  G4double ekin =
718  std::min(backwardEnergy-backwardKinetic, centerofmassEnergy/2.0-protonMass/GeV);
719 
720  if( ekin < 0.04 )ekin = 0.04 * std::fabs( normal() );
721  vecMass = targetParticle.GetMass()/GeV;
722  totalEnergy = ekin + vecMass;
723  targetParticle.SetTotalEnergy( totalEnergy*GeV );
724  pp = std::sqrt(std::abs(totalEnergy*totalEnergy - vecMass*vecMass) )*GeV;
725  pp1 = pseudoParticle[6].GetMomentum().mag();
726  if (pp1 < 1.0e-6*GeV) {
727  G4ThreeVector iso = Isotropic(pp);
728  targetParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
729  } else {
730  targetParticle.SetMomentum( pseudoParticle[6].GetMomentum() * (pp/pp1));
731  }
732  pseudoParticle[5] = pseudoParticle[5] + targetParticle;
733 
734  } else if (backwardNucleonCount > 1) {
735  // Share remaining energy with up to 17 backward nucleons
736 
737  G4int tempCount = 5;
738  if (backwardNucleonCount < 5) tempCount = backwardNucleonCount;
739  tempCount -= 2;
740 
741  G4double clusterMass = 0.;
742  if (targetParticle.GetSide() == -3)
743  clusterMass = targetParticle.GetMass()/GeV;
744  for (i = 0; i < vecLen; ++i)
745  if (vec[i]->GetSide() == -3) clusterMass += vec[i]->GetMass()/GeV;
746  clusterMass += backwardEnergy - backwardKinetic;
747 
748  totalEnergy = pseudoParticle[6].GetTotalEnergy()/GeV;
749  pseudoParticle[6].SetMass(clusterMass*GeV);
750 
751  pp = std::sqrt(std::abs(totalEnergy*totalEnergy -
752  clusterMass*clusterMass) )*GeV;
753  pp1 = pseudoParticle[6].GetMomentum().mag();
754  if (pp1 < 1.0e-6*GeV) {
755  G4ThreeVector iso = Isotropic(pp);
756  pseudoParticle[6].SetMomentum(iso.x(), iso.y(), iso.z());
757  } else {
758  pseudoParticle[6].SetMomentum(pseudoParticle[6].GetMomentum() * (-pp/pp1));
759  }
760 
761  std::vector<G4ReactionProduct*> tempList; // ptrs to backward nucleons
762  if (targetParticle.GetSide() == -3) tempList.push_back(&targetParticle);
763  for (i = 0; i < vecLen; ++i)
764  if (vec[i]->GetSide() == -3) tempList.push_back(vec[i]);
765 
766  constantCrossSection = true;
767 
768  if (tempList.size() > 1) {
769  G4int n_entry = 0;
770  wgt = GenerateNBodyEventT(pseudoParticle[6].GetMass(),
771  constantCrossSection, tempList);
772 
773  if (targetParticle.GetSide() == -3) {
774  targetParticle = *tempList[0];
775  targetParticle.Lorentz(targetParticle, pseudoParticle[6]);
776  n_entry++;
777  }
778 
779  for (i = 0; i < vecLen; ++i) {
780  if (vec[i]->GetSide() == -3) {
781  *vec[i] = *tempList[n_entry];
782  vec[i]->Lorentz(*vec[i], pseudoParticle[6]);
783  n_entry++;
784  }
785  }
786  }
787  } else return false;
788 
789  if (vecLen == 0) return false; // all the secondaries have been eliminated
790 
791  // Lorentz transformation to lab frame
792 
793  currentParticle.Lorentz( currentParticle, pseudoParticle[1] );
794  targetParticle.Lorentz( targetParticle, pseudoParticle[1] );
795  for (i = 0; i < vecLen; ++i) vec[i]->Lorentz(*vec[i], pseudoParticle[1]);
796 
797  // Set leading strange particle flag.
798  // leadFlag will be true if original particle and incident particle are
799  // both strange, in which case the incident particle becomes the leading
800  // particle.
801  // leadFlag will also be true if the target particle is strange, but the
802  // incident particle is not, in which case the target particle becomes the
803  // leading particle.
804 
805  G4bool leadingStrangeParticleHasChanged = true;
806  if (leadFlag)
807  {
808  if (currentParticle.GetDefinition() == leadingStrangeParticle.GetDefinition())
809  leadingStrangeParticleHasChanged = false;
810  if (leadingStrangeParticleHasChanged &&
811  (targetParticle.GetDefinition() == leadingStrangeParticle.GetDefinition()) )
812  leadingStrangeParticleHasChanged = false;
813  if( leadingStrangeParticleHasChanged )
814  {
815  for( i=0; i<vecLen; i++ )
816  {
817  if( vec[i]->GetDefinition() == leadingStrangeParticle.GetDefinition() )
818  {
819  leadingStrangeParticleHasChanged = false;
820  break;
821  }
822  }
823  }
824  if( leadingStrangeParticleHasChanged )
825  {
826  G4bool leadTest =
827  (leadingStrangeParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
828  leadingStrangeParticle.GetDefinition()->GetParticleSubType() == "pi");
829  G4bool targetTest =
830  (targetParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
831  targetParticle.GetDefinition()->GetParticleSubType() == "pi");
832 
833  // following modified by JLC 22-Oct-97
834 
835  if( (leadTest&&targetTest) || !(leadTest||targetTest) ) // both true or both false
836  {
837  targetParticle.SetDefinitionAndUpdateE( leadingStrangeParticle.GetDefinition() );
838  targetHasChanged = true;
839  }
840  else
841  {
842  currentParticle.SetDefinitionAndUpdateE( leadingStrangeParticle.GetDefinition() );
843  incidentHasChanged = false;
844  }
845  }
846  } // end of if( leadFlag )
847 
848  // Get number of final state nucleons and nucleons remaining in
849  // target nucleus
850 
851  std::pair<G4int, G4int> finalStateNucleons =
852  GetFinalStateNucleons(originalTarget, vec, vecLen);
853 
854  G4int protonsInFinalState = finalStateNucleons.first;
855  G4int neutronsInFinalState = finalStateNucleons.second;
856 
857  G4int numberofFinalStateNucleons =
858  protonsInFinalState + neutronsInFinalState;
859 
860  if (currentParticle.GetDefinition()->GetBaryonNumber() == 1 &&
861  targetParticle.GetDefinition()->GetBaryonNumber() == 1 &&
862  originalIncident->GetDefinition()->GetPDGMass() <
864  numberofFinalStateNucleons++;
865 
866  numberofFinalStateNucleons = std::max(1, numberofFinalStateNucleons);
867 
868  G4int PinNucleus = std::max(0,
869  G4int(targetNucleus.GetZ_asInt()) - protonsInFinalState);
870  G4int NinNucleus = std::max(0,
871  G4int(targetNucleus.GetA_asInt()-targetNucleus.GetZ_asInt()) - neutronsInFinalState);
872 
873  pseudoParticle[3].SetMomentum( 0.0, 0.0, pOriginal*GeV );
874  pseudoParticle[3].SetMass( mOriginal*GeV );
875  pseudoParticle[3].SetTotalEnergy(
876  std::sqrt( pOriginal*pOriginal + mOriginal*mOriginal )*GeV );
877 
878  G4ParticleDefinition * aOrgDef = modifiedOriginal.GetDefinition();
879  G4int diff = 0;
880  if(aOrgDef == G4Proton::Proton() || aOrgDef == G4Neutron::Neutron() ) diff = 1;
881  if(numberofFinalStateNucleons == 1) diff = 0;
882  pseudoParticle[4].SetMomentum( 0.0, 0.0, 0.0 );
883  pseudoParticle[4].SetMass( protonMass*(numberofFinalStateNucleons-diff) );
884  pseudoParticle[4].SetTotalEnergy( protonMass*(numberofFinalStateNucleons-diff) );
885 
886  G4double theoreticalKinetic =
887  pseudoParticle[3].GetTotalEnergy() + pseudoParticle[4].GetTotalEnergy() -
888  currentParticle.GetMass() - targetParticle.GetMass();
889  for (i = 0; i < vecLen; ++i) theoreticalKinetic -= vec[i]->GetMass();
890 
891  G4double simulatedKinetic =
892  currentParticle.GetKineticEnergy() + targetParticle.GetKineticEnergy();
893  for (i = 0; i < vecLen; ++i)
894  simulatedKinetic += vec[i]->GetKineticEnergy();
895 
896  pseudoParticle[5] = pseudoParticle[3] + pseudoParticle[4];
897  pseudoParticle[3].Lorentz( pseudoParticle[3], pseudoParticle[5] );
898  pseudoParticle[4].Lorentz( pseudoParticle[4], pseudoParticle[5] );
899 
900  pseudoParticle[7].SetZero();
901  pseudoParticle[7] = pseudoParticle[7] + currentParticle;
902  pseudoParticle[7] = pseudoParticle[7] + targetParticle;
903  for (i = 0; i < vecLen; ++i)
904  pseudoParticle[7] = pseudoParticle[7] + *vec[i];
905 
906  /*
907  // This code does not appear to do anything to the energy or angle spectra
908  if( vecLen <= 16 && vecLen > 0 )
909  {
910  // must create a new set of ReactionProducts here because GenerateNBody will
911  // modify the momenta for the particles, and we don't want to do this
912  //
913  G4ReactionProduct tempR[130];
914  tempR[0] = currentParticle;
915  tempR[1] = targetParticle;
916  for( i=0; i<vecLen; ++i )tempR[i+2] = *vec[i];
917  G4FastVector<G4ReactionProduct,256> tempV1;
918  tempV1.Initialize( vecLen+2 );
919  G4int tempLen1 = 0;
920  for( i=0; i<vecLen+2; ++i )tempV1.SetElement( tempLen1++, &tempR[i] );
921  constantCrossSection = true;
922 
923  wgt = GenerateNBodyEvent(pseudoParticle[3].GetTotalEnergy() +
924  pseudoParticle[4].GetTotalEnergy(),
925  constantCrossSection, tempV1, tempLen1);
926  if (wgt == -1) {
927  G4double Qvalue = 0;
928  for (i = 0; i < tempLen1; i++) Qvalue += tempV1[i]->GetMass();
929  wgt = GenerateNBodyEvent(Qvalue,
930  constantCrossSection, tempV1, tempLen1);
931  }
932  if(wgt>-.5)
933  {
934  theoreticalKinetic = 0.0;
935  for( i=0; i<tempLen1; ++i )
936  {
937  pseudoParticle[6].Lorentz( *tempV1[i], pseudoParticle[4] );
938  theoreticalKinetic += pseudoParticle[6].GetKineticEnergy();
939  }
940  }
941  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
942  }
943  */
944 
945  //
946  // Make sure that the kinetic energies are correct
947  //
948 
949  if (simulatedKinetic != 0.0) {
950  wgt = theoreticalKinetic/simulatedKinetic;
951  theoreticalKinetic = currentParticle.GetKineticEnergy() * wgt;
952  simulatedKinetic = theoreticalKinetic;
953  currentParticle.SetKineticEnergy(theoreticalKinetic);
954  pp = currentParticle.GetTotalMomentum();
955  pp1 = currentParticle.GetMomentum().mag();
956  if (pp1 < 1.0e-6*GeV) {
957  G4ThreeVector iso = Isotropic(pp);
958  currentParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
959  } else {
960  currentParticle.SetMomentum(currentParticle.GetMomentum() * (pp/pp1));
961  }
962 
963  theoreticalKinetic = targetParticle.GetKineticEnergy() * wgt;
964  targetParticle.SetKineticEnergy(theoreticalKinetic);
965  simulatedKinetic += theoreticalKinetic;
966  pp = targetParticle.GetTotalMomentum();
967  pp1 = targetParticle.GetMomentum().mag();
968 
969  if (pp1 < 1.0e-6*GeV) {
970  G4ThreeVector iso = Isotropic(pp);
971  targetParticle.SetMomentum(iso.x(), iso.y(), iso.z() );
972  } else {
973  targetParticle.SetMomentum(targetParticle.GetMomentum() * (pp/pp1) );
974  }
975 
976  for (i = 0; i < vecLen; ++i ) {
977  theoreticalKinetic = vec[i]->GetKineticEnergy() * wgt;
978  simulatedKinetic += theoreticalKinetic;
979  vec[i]->SetKineticEnergy(theoreticalKinetic);
980  pp = vec[i]->GetTotalMomentum();
981  pp1 = vec[i]->GetMomentum().mag();
982  if( pp1 < 1.0e-6*GeV ) {
983  G4ThreeVector iso = Isotropic(pp);
984  vec[i]->SetMomentum(iso.x(), iso.y(), iso.z() );
985  } else {
986  vec[i]->SetMomentum(vec[i]->GetMomentum() * (pp/pp1) );
987  }
988  }
989  }
990 
991  // Rotate(numberofFinalStateNucleons, pseudoParticle[3].GetMomentum(),
992  // modifiedOriginal, originalIncident, targetNucleus,
993  // currentParticle, targetParticle, vec, vecLen );
994 
995  // Add black track particles
996  // the total number of particles produced is restricted to 198
997  // this may have influence on very high energies
998 
999  if( atomicWeight >= 1.5 )
1000  {
1001  // npnb is number of proton/neutron black track particles
1002  // ndta is the number of deuterons, tritons, and alphas produced
1003  // epnb is the kinetic energy available for proton/neutron black track
1004  // particles
1005  // edta is the kinetic energy available for deuteron/triton/alpha particles
1006 
1007  G4int npnb = 0;
1008  G4int ndta = 0;
1009 
1010  G4double epnb, edta;
1011  if (veryForward) {
1012  epnb = targetNucleus.GetAnnihilationPNBlackTrackEnergy();
1013  edta = targetNucleus.GetAnnihilationDTABlackTrackEnergy();
1014  } else {
1015  epnb = targetNucleus.GetPNBlackTrackEnergy();
1016  edta = targetNucleus.GetDTABlackTrackEnergy();
1017  }
1018  /*
1019  G4ReactionProduct* fudge = new G4ReactionProduct();
1020  fudge->SetDefinition( aProton );
1021  G4double TT = epnb + edta;
1022  G4double MM = fudge->GetMass()/GeV;
1023  fudge->SetTotalEnergy(MM*GeV + TT*GeV);
1024  G4double pzz = std::sqrt(TT*(TT + 2.*MM));
1025  fudge->SetMomentum(0.0, 0.0, pzz*GeV);
1026  vec.SetElement(vecLen++, fudge);
1027  // G4cout << " Fudge = " << vec[vecLen-1]->GetKineticEnergy()/GeV
1028  << G4endl;
1029  */
1030 
1031  const G4double pnCutOff = 0.001;
1032  const G4double dtaCutOff = 0.001;
1033  // const G4double kineticMinimum = 1.e-6;
1034  // const G4double kineticFactor = -0.010;
1035  // G4double sprob = 0.0; // sprob = probability of self-absorption in
1036  // heavy molecules
1037  // Not currently used (DHW 9 June 2008) const G4double ekIncident = originalIncident->GetKineticEnergy()/GeV;
1038  // if (ekIncident >= 5.0) sprob = std::min(1.0, 0.6*std::log(ekIncident-4.0));
1039  if (epnb > pnCutOff)
1040  {
1041  npnb = G4Poisson((1.5+1.25*numberofFinalStateNucleons)*epnb/(epnb+edta));
1042  if (numberofFinalStateNucleons + npnb > atomicWeight)
1043  npnb = G4int(atomicWeight+0.00001 - numberofFinalStateNucleons);
1044  npnb = std::min( npnb, 127-vecLen );
1045  }
1046  if( edta >= dtaCutOff )
1047  {
1048  ndta = G4Poisson((1.5+1.25*numberofFinalStateNucleons)*edta/(epnb+edta));
1049  ndta = std::min( ndta, 127-vecLen );
1050  }
1051  if (npnb == 0 && ndta == 0) npnb = 1;
1052 
1053  AddBlackTrackParticles(epnb, npnb, edta, ndta, modifiedOriginal,
1054  PinNucleus, NinNucleus, targetNucleus,
1055  vec, vecLen);
1056  }
1057 
1058  // if( centerofmassEnergy <= (4.0+G4UniformRand()) )
1059  // MomentumCheck( modifiedOriginal, currentParticle, targetParticle,
1060  // vec, vecLen );
1061  //
1062  // calculate time delay for nuclear reactions
1063  //
1064 
1065  if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) )
1066  currentParticle.SetTOF(
1067  1.0-500.0*std::exp(-ekOriginal/0.04)*std::log(G4UniformRand()) );
1068  else
1069  currentParticle.SetTOF( 1.0 );
1070  return true;
1071 
1072 }
1073 
1074 
1075 void G4RPGFragmentation::
1076 ReduceEnergiesOfSecondaries(G4int startingIndex,
1077  G4double& forwardKinetic,
1078  G4double& backwardKinetic,
1080  G4int& vecLen,
1081  G4ReactionProduct& forwardPseudoParticle,
1082  G4ReactionProduct& backwardPseudoParticle,
1083  G4double& pt)
1084 {
1085  // Reduce energies and pt of secondaries in order to maintain
1086  // energy conservation
1087 
1088  G4double totalEnergy;
1089  G4double pp;
1090  G4double pp1;
1091  G4double px;
1092  G4double py;
1093  G4double mass;
1094  G4ReactionProduct* pVec;
1095  G4int i;
1096 
1097  forwardKinetic = 0.0;
1098  backwardKinetic = 0.0;
1099  forwardPseudoParticle.SetZero();
1100  backwardPseudoParticle.SetZero();
1101 
1102  for (i = startingIndex; i < vecLen; i++) {
1103  pVec = vec[i];
1104  if (pVec->GetSide() != -3) {
1105  mass = pVec->GetMass();
1106  totalEnergy = 0.95*pVec->GetTotalEnergy() + 0.05*mass;
1107  pVec->SetTotalEnergy(totalEnergy);
1108  pp = std::sqrt( std::abs( totalEnergy*totalEnergy - mass*mass ) );
1109  pp1 = pVec->GetMomentum().mag();
1110  if (pp1 < 1.0e-6*GeV) {
1111  G4ThreeVector iso = Isotropic(pp);
1112  pVec->SetMomentum( iso.x(), iso.y(), iso.z() );
1113  } else {
1114  pVec->SetMomentum(pVec->GetMomentum() * (pp/pp1) );
1115  }
1116 
1117  px = pVec->GetMomentum().x();
1118  py = pVec->GetMomentum().y();
1119  pt = std::max(1.0, std::sqrt( px*px + py*py ) )/GeV;
1120  if (pVec->GetSide() > 0) {
1121  forwardKinetic += pVec->GetKineticEnergy()/GeV;
1122  forwardPseudoParticle = forwardPseudoParticle + (*pVec);
1123  } else {
1124  backwardKinetic += pVec->GetKineticEnergy()/GeV;
1125  backwardPseudoParticle = backwardPseudoParticle + (*pVec);
1126  }
1127  }
1128  }
1129 }
1130 
1131  /* end of file */