Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4ParticleHPEnAngCorrelation.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 // particle_hp -- source file
27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron transport model.
29 //
30 // 100413 Fix bug in incidence energy by T. Koi
31 //
32 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
33 //
35 #include "G4LorentzRotation.hh"
36 #include "G4LorentzVector.hh"
37 #include "G4RotationMatrix.hh"
38 #include "G4IonTable.hh"
39 
41 {
43 
44  // do we have an appropriate distribution
45  if(nProducts!=1) throw G4HadronicException(__FILE__, __LINE__, "More than one product in SampleOne");
46 
47  // get the result
48  G4ReactionProductVector * temp=0;
49  G4int i=0;
50 
51  G4int icounter=0;
52  G4int icounter_max=1024;
53  while(temp == 0) {
54  icounter++;
55  if ( icounter > icounter_max ) {
56  G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
57  break;
58  }
59  temp = theProducts[i++].Sample(anEnergy,1);
60  }
61 
62  // is the multiplicity correct
63  if(temp->size()!=1) throw G4HadronicException(__FILE__, __LINE__, "SampleOne: Yield not correct");
64 
65  // fill result
66  result = temp->operator[](0);
67 
68  // some garbage collection
69  delete temp;
70 
71  // return result
72  return result;
73 }
74 
76 {
78  G4int i;
80  G4ReactionProduct theCMS;
82 
83  if(frameFlag==2
84  || frameFlag==3) // Added for particle HP
85  {
86  // simplify and double check @
87  G4ThreeVector the3IncidentPart = fCache.Get().theProjectileRP->GetMomentum(); //theProjectileRP has value in LAB
88  G4double nEnergy = fCache.Get().theProjectileRP->GetTotalEnergy();
89  G4ThreeVector the3Target = fCache.Get().theTarget->GetMomentum(); //theTarget has value in LAB
90  G4double tEnergy = fCache.Get().theTarget->GetTotalEnergy();
91  G4double totE = nEnergy+tEnergy;
92  G4ThreeVector the3CMS = the3Target+the3IncidentPart;
93  theCMS.SetMomentum(the3CMS);
94  G4double cmsMom = std::sqrt(the3CMS*the3CMS);
95  G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
96  theCMS.SetMass(sqrts);
97  theCMS.SetTotalEnergy(totE);
98  G4ReactionProduct aIncidentPart;
99  aIncidentPart.Lorentz(*fCache.Get().theProjectileRP, theCMS);
100  //TKDB 100413
101  //ENDF-6 Formats Manual ENDF-102
102  //CHAPTER 6. FILE 6: PRODUCT ENERGY-ANGLE DISTRIBUTIONS
103  //LCT Reference system for secondary energy and angle (incident energy is always given in the LAB system)
104  //anEnergy = aIncidentPart.GetKineticEnergy();
105  anEnergy = fCache.Get().theProjectileRP->GetKineticEnergy(); //should be same argumment of "anEnergy"
106 
107  G4LorentzVector Ptmp (aIncidentPart.GetMomentum(), aIncidentPart.GetTotalEnergy());
108 
109  toZ.rotateZ(-1*Ptmp.phi());
110  toZ.rotateY(-1*Ptmp.theta());
111  }
112  fCache.Get().theTotalMeanEnergy=0;
113  G4LorentzRotation toLab(toZ.inverse()); //toLab only change axis NOT to LAB system
114  //- get first number of particles, to check if sum of Z and N is not bigger than target values
115  std::vector<int> nParticles;
116  bool bNPOK = true;
117 //TKDB_PHP_150507
118 #ifdef PHP_AS_HP
119 #endif
120 //TKDB_PHP_161107
121  G4int iTry(0);
122 //TKDB_PHP_161107
123 //TKDB_PHP_150507
124  do {
125  G4int sumZ = 0;
126  G4int sumA = 0;
127  nParticles.clear();
128  for(i=0; i<nProducts; i++)
129  {
130  G4int massCode = G4int(theProducts[i].GetMassCode());
131  G4int nPart;
132  nPart = theProducts[i].GetMultiplicity(anEnergy);
133  sumZ += massCode/1000 * nPart;
134  sumA += massCode % 1000 * nPart;
135 #ifdef G4PHPDEBUG
136  if( getenv("G4ParticleHPDebug") ) G4cout << i << " G4ParticleHPEnAngCorrelation::MULTIPLICITY " << massCode << " sumZ " << sumZ << " sumA " << sumA << " NPART " << nPart << G4endl;
137 #endif
138  nParticles.push_back( nPart );
139  }
140  bNPOK = true;
141  double targetZ = fCache.Get().theTarget->GetDefinition()->GetAtomicNumber();
142  double targetA = fCache.Get().theTarget->GetDefinition()->GetAtomicMass();
143  targetZ += fCache.Get().theProjectileRP->GetDefinition()->GetAtomicNumber();
144  targetA += fCache.Get().theProjectileRP->GetDefinition()->GetAtomicMass();
145  if ( bAdjustFinalState ) {
146 /*
147 G4cout << "TKDB G4ParticleHPEnAngCorrelation::Sample 1" << G4endl;
148 G4cout << "TKDB "
149 << "targetZ = " << targetZ
150 << ", targetA = " << targetA
151 << ", sumZ = " << sumZ
152 << ", sumA = " << sumA
153 << ", int( targetZ-sumZ ) = " << int( targetZ-sumZ )
154 << ", int( targetA-sumA ) = " << int( targetA-sumA )
155 //<< ", G4IonTable::GetIonTable()->GetIon ( int(targetZ - sumZ), (int)(targetA - sumA), 0.0 ) = " << G4IonTable::GetIonTable()->GetIon ( int(targetZ - sumZ), (int)(targetA - sumA), 0.0 )
156 << G4endl;
157 */
158  //if ( (sumZ != targetZ || sumA != targetA ) &&
159  // (sumZ > targetZ || sumA > targetA
160  // || ! G4IonTable::GetIonTable()->GetIon ( int(targetZ - sumZ), (int)(targetA - sumA), 0.0 ) ) ){ // e.g. Z=3, A=2
161  if ( ( sumZ != targetZ || sumA != targetA )
162  && ( sumZ > targetZ || sumA > targetA || (targetZ-sumZ) >= (targetA-sumA) ) ) {
163  // e.g. Z=3, A=2
164  bNPOK = false;
165  //nParticles.clear();
166 #ifdef G4PHPDEBUG
167  if ( getenv("G4ParticleHPDebug") )
168  G4cerr << " WRONG MULTIPLICITY Z= " << sumZ
169  << " > " << targetZ
170  << " A= " << sumA
171  << " > " << targetA << G4endl;
172 #endif
173  }
174  }
175 //TKDB_PHP_150507
176 #ifdef PHP_AS_HP
177 #endif
178 //TKDB_PHP_161107
179  iTry++;
180  if ( iTry > 1024 ) {
181  G4Exception("G4ParticleHPEnAngCorrelation::Sample",
182  "Warning",
183  JustWarning,
184  "Too many trials were done. Exiting current loop by force. You may have Probably, the result violating (baryon number) conservation law will be obtained.");
185  bNPOK=true;
186  }
187 //TKDB_PHP_161107
188 //TKDB_PHP_150507
189 
190  }while(!bNPOK); // Loop checking, 11.05.2015, T. Koi
191 
192  for(i=0; i<nProducts; i++)
193  {
194  //- if( nParticles[i] == 0 ) continue;
195  it = theProducts[i].Sample(anEnergy,nParticles[i]);
197  // if( getenv("G4PHPTEST") ) G4cout << " EnAnG energy sampled " << it->operator[](0)->GetKineticEnergy() << " aMeanEnergy " << aMeanEnergy << G4endl; // GDEB
198  //if(aMeanEnergy>0)
199  //151120 TK Modified for solving reproducibility problem
200  //This change may have side effect.
201  if(aMeanEnergy>=0)
202  {
203  fCache.Get().theTotalMeanEnergy += aMeanEnergy;
204  }
205  else
206  {
207  fCache.Get().theTotalMeanEnergy = anEnergy/nProducts+theProducts[i].GetQValue();
208  }
209  if(it!=0)
210  {
211  for(unsigned int ii=0; ii<it->size(); ii++)
212  {
213  //if(!getenv("G4PHP_NO_LORENTZ_BOOST")) {
214  G4LorentzVector pTmp1 (it->operator[](ii)->GetMomentum(),
215  it->operator[](ii)->GetTotalEnergy());
216  pTmp1 = toLab*pTmp1;
217  if( getenv("G4PHPTEST") ) G4cout << " G4particleHPEnAngCorrelation COS THETA " << std::cos(it->operator[](ii)->GetMomentum().theta()) << G4endl;
218  it->operator[](ii)->SetMomentum(pTmp1.vect());
219  it->operator[](ii)->SetTotalEnergy(pTmp1.e());
220  if( getenv("G4PHPTEST") ) G4cout << " G4particleHPEnAngCorrelation COS THETA after toLab " << std::cos(it->operator[](ii)->GetMomentum().theta()) << G4endl;
221 
222  if(frameFlag==1) // target rest //TK 100413 should be LAB?
223  {
224  it->operator[](ii)->Lorentz(*(it->operator[](ii)), -1.*(*fCache.Get().theTarget)); //TK 100413 Is this really need?
225  }
226  else if(frameFlag==2 ) // CMS
227  {
228 #ifdef G4PHPDEBUG
229  if( getenv("G4ParticleHPDebug") )
230  G4cout <<"G4ParticleHPEnAngCorrelation: before Lorentz boost "<<
231  it->at(ii)->GetKineticEnergy()<<" "<<
232  it->at(ii)->GetMomentum()<<G4endl;
233 #endif
234  it->operator[](ii)->Lorentz(*(it->operator[](ii)), -1.*theCMS);
235 #ifdef G4PHPDEBUG
236  if( getenv("G4ParticleHPDebug") )
237  G4cout <<"G4ParticleHPEnAngCorrelation: after Lorentz boost "<<
238  it->at(ii)->GetKineticEnergy()<<" "<<
239  it->at(ii)->GetMomentum()<<G4endl;
240 #endif
241  }
242  //TK120515 migrate frameFlag (MF6 LCT) = 3
243  else if(frameFlag==3) // CMS A<=4 other LAB
244  {
245  if ( theProducts[i].GetMassCode() > 4 ) //Alpha AWP 3.96713
246  {
247  //LAB
248  it->operator[](ii)->Lorentz(*(it->operator[](ii)), -1.*(*fCache.Get().theTarget)); //TK 100413 Is this really need?
249 #ifdef G4PHPDEBUG
250  if( getenv("G4ParticleHPDebug") )
251  G4cout <<"G4ParticleHPEnAngCorrelation: after Lorentz boost "<<
252  it->at(ii)->GetKineticEnergy()<<" "<<
253  it->at(ii)->GetMomentum()<<G4endl;
254 #endif
255  }
256  else
257  {
258  //CMS
259  it->operator[](ii)->Lorentz(*(it->operator[](ii)), -1.*theCMS);
260 #ifdef G4PHPDEBUG
261  if( getenv("G4ParticleHPDebug") )
262  G4cout <<"G4ParticleHPEnAngCorrelation: after Lorentz boost "<<
263  it->at(ii)->GetKineticEnergy()<<" "<<
264  it->at(ii)->GetMomentum()<<G4endl;
265 #endif
266  }
267  }
268  else
269  {
270  throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPEnAngCorrelation::Sample: The frame of the finalstate is not specified");
271  }
272  if( getenv("G4PHPTEST") ) G4cout << frameFlag << " G4particleHPEnAngCorrelation COS THETA after Lorentz " << std::cos(it->operator[](ii)->GetMomentum().theta()) << G4endl;
273 
274  // }//getenv("G4PHP_NO_LORENTZ_BOOST"))
275  // G4cout << ii << " EnAnG energy after boost " << it->operator[](ii)->GetKineticEnergy() << G4endl; //GDEB
276  result->push_back(it->operator[](ii));
277  }
278  delete it;
279  }
280  }
281 
282  return result;
283 }
284 
G4int GetMultiplicity(G4double anEnergy)
G4ReactionProductVector * Sample(G4double anEnergy, G4int nParticles)
void SetMomentum(const G4double x, const G4double y, const G4double z)
#define G4endl
Definition: G4ios.hh:61
G4ReactionProductVector * Sample(G4double anEnergy)
G4ReactionProduct * SampleOne(G4double anEnergy)
void SetTotalEnergy(const G4double en)
void SetMass(const G4double mas)
G4double GetTotalEnergy() const
double G4double
Definition: G4Types.hh:76
HepLorentzRotation & rotateZ(double delta)
HepLorentzRotation inverse() const
G4ThreeVector GetMomentum() const
std::vector< G4ReactionProduct * > G4ReactionProductVector
G4GLOB_DLL std::ostream G4cerr
G4double G4ParticleHPJENDLHEData::G4double result
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
G4double MeanEnergyOfThisInteraction()
HepLorentzRotation & rotateY(double delta)
void Lorentz(const G4ReactionProduct &p1, const G4ReactionProduct &p2)