ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPDeExGammas.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ParticleHPDeExGammas.cc
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 // neutron_hp -- source file
27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron transport model.
29 //
30 // 080801 Prohibit level transition to oneself by T. Koi
31 //
32 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
33 //
35 #include "G4SystemOfUnits.hh"
36 
37 void G4ParticleHPDeExGammas::Init(std::istream & aDataFile)
38 {
39  // G4cout << this << "ExGammas Init LEVEL " << G4endl; //GDEB
40  G4ParticleHPGamma ** theGammas = new G4ParticleHPGamma * [50];
41  G4int nGammas = 0;
42  G4int nBuff = 50;
43  for(;;)
44  {
45  G4ParticleHPGamma * theNew = new G4ParticleHPGamma;
46  if(!theNew->Init(aDataFile))
47  {
48  delete theNew;
49  break;
50  }
51  else
52  {
53  if(nGammas==nBuff)
54  {
55  nBuff+=50;
56  G4ParticleHPGamma ** buffer = new G4ParticleHPGamma * [nBuff];
57  for(G4int i=0;i<nGammas;i++) buffer[i] = theGammas[i];
58  delete [] theGammas;
59  theGammas = buffer;
60  }
61  theGammas[nGammas] = theNew;
62  nGammas++;
63  }
64  }
65  // all gammas are in. Now sort them into levels.
66 
67  // count the levels
68 
69  G4double currentE = 0;
70  G4double nextE = 0;
71  G4int i;
72  G4double epsilon = 0.01*keV;
73  for(i=0; i<nGammas; i++)
74  {
75  nextE = theGammas[i]->GetLevelEnergy();
76  if(std::abs(currentE-nextE)>epsilon) nLevels++;
77  currentE = nextE;
78  }
79 
80  // G4cout << this << "LEVEL " << nLevels << G4endl; //GDEB
81  // Build the levels
82 
84  levelStart = new G4int [nLevels];
85  levelSize = new G4int [nLevels];
86 
87  // fill the levels
88 
89  currentE = 0;
90  nextE = 0;
91  G4int levelCounter=-1;
92  for(i=0; i<nGammas; i++)
93  {
94  nextE = theGammas[i]->GetLevelEnergy();
95  if(std::abs(currentE-nextE)>epsilon)
96  {
97  levelCounter++;
98  levelStart[levelCounter] = i;
99  levelSize[levelCounter] = 0;
100  }
101  levelSize[levelCounter]++;
102  currentE = nextE;
103  }
104 
105  for(i=0; i<nLevels; i++)
106  {
108  for(G4int ii=levelStart[i]; ii<levelStart[i]+levelSize[i]; ii++)
109  {
110  theLevels[i].SetGamma(ii-levelStart[i], theGammas[ii]);
111  }
112  }
113 
114 // set the next relation in the gammas.
115  G4double levelE, gammaE, currentLevelE;
116  G4double min;
117  for(i=0; i<nGammas; i++)
118  {
119  G4int it=-1;
120  gammaE = theGammas[i]->GetGammaEnergy();
121  currentLevelE = theGammas[i]->GetLevelEnergy();
122  min = currentLevelE-gammaE-epsilon;
123  for(G4int ii=0; ii<nLevels; ii++)
124  {
125  levelE = theLevels[ii].GetLevelEnergy();
126  if(std::abs(currentLevelE-(levelE+gammaE))<min)
127  {
128  min = std::abs(currentLevelE-(levelE+gammaE));
129  it = ii;
130  }
131  }
132 //080728
133  if ( it != -1 && currentLevelE == theLevels[it].GetLevelEnergy() )
134  {
135  //TK Comment; Some data file in /Inelastic/Gammas has inconsistent level data (no level to transit)
136  //G4cout << "DeExGammas Transition level error: it " << it << " " << currentLevelE << " " << gammaE << " " << theLevels[it-1].GetLevelEnergy() << " " << currentLevelE - theLevels[it-1].GetLevelEnergy() << G4endl;
137  // Forced to connect the next(previous) level
138  it +=-1;
139  }
140 //080728
141  if(it!=-1) theGammas[i]->SetNext(&theLevels[it]);
142  }
143  // some garbage collection
144 
145  delete [] theGammas;
146 
147  // and we are Done.
148 }