source: trunk/source/processes/hadronic/models/binary_cascade/test/BinaryLightIonReaction.cc@ 1228

Last change on this file since 1228 was 1199, checked in by garnier, 16 years ago

nvx fichiers dans CVS
File size: 12.9 KB
Line 
1//
2// ********************************************************************
3// * DISCLAIMER *
4// * *
5// * The following disclaimer summarizes all the specific disclaimers *
6// * of contributors to this software. The specific disclaimers,which *
7// * govern, are listed with their locations in: *
8// * http://cern.ch/geant4/license *
9// * *
10// * Neither the authors of this software system, nor their employing *
11// * institutes,nor the agencies providing financial support for this *
12// * work make any representation or warranty, express or implied, *
13// * regarding this software system or assume any liability for its *
14// * use. *
15// * *
16// * This code implementation is the intellectual property of the *
17// * GEANT4 collaboration. *
18// * By copying, distributing or modifying the Program (or any work *
19// * based on the Program) you indicate your acceptance of this *
20// * statement, and all its terms. *
21// ********************************************************************
22//
23//
24// $Id: BinaryLightIonReaction.cc,v 1.11 2004/06/03 15:39:54 hpw Exp $
25// GEANT4 tag $Name: geant4-09-03-cand-01 $
26//
27// Johannes Peter Wellisch, 22.Apr 1997: full test-suite coded.
28#include "G4ios.hh"
29#include <fstream>
30#include <iomanip>
31
32#include "G4Material.hh"
33
34#include "G4GRSVolume.hh"
35#include "G4ProcessManager.hh"
36#include "G4HadronInelasticProcess.hh"
37
38#include "G4IonInelasticProcess.hh"
39#include "G4BinaryLightIonReaction.hh"
40
41#include "G4DynamicParticle.hh"
42#include "G4LeptonConstructor.hh"
43#include "G4BaryonConstructor.hh"
44#include "G4MesonConstructor.hh"
45#include "G4IonConstructor.hh"
46
47#include "G4Box.hh"
48#include "G4PVPlacement.hh"
49
50#include "G4Step.hh"
51#include "G4StepPoint.hh"
52
53#include "G4TripathiCrossSection.hh"
54
55 // forward declarations
56
57 G4int sortEnergies( const double Px, const double Py, const double Pz,
58 const double Ekin, double* sortedPx, double* sortedPy,
59 double* sortedPz, double* sortedE );
60
61 // here comes the code
62
63 G4int sortEnergies( const double Px, const double Py, const double Pz,
64 const double Ekin, double* sortedPx, double* sortedPy,
65 double* sortedPz, double* sortedE)
66 {
67 for( int i=0; i<10; ++i ) {
68 if( abs(Ekin) > sortedE[i] ) {
69 sortedE[i] = Ekin;
70 sortedPx[i] = Px;
71 sortedPy[i] = Py;
72 sortedPz[i] = Pz;
73 return 1;
74 }
75 }
76 return 0;
77 }
78
79// extern "C" int ecpt();
80
81 int main()
82 {
83
84// ecpt();
85
86 G4cout.setf( std::ios::scientific, std::ios::floatfield );
87 std::ofstream outFile( "InInelasticAlpha.listing.GetMeanFreePath", std::ios::out);
88 outFile.setf( std::ios::scientific, std::ios::floatfield );
89 std::ofstream outFile1( "InInelasticAlpha.listing.DoIt", std::ios::out);
90 outFile1.setf( std::ios::scientific, std::ios::floatfield );
91
92 G4String name, symbol;
93 G4double a, iz, z, density;
94 G4int nEl, nIso;
95
96 // constructing the particles
97
98 G4LeptonConstructor aC1;
99 G4BaryonConstructor aC2;
100 G4MesonConstructor aC3;
101 G4IonConstructor aC4;
102
103 aC1.ConstructParticle();
104 aC2.ConstructParticle();
105 aC3.ConstructParticle();
106 aC4.ConstructParticle();
107
108 G4int numberOfMaterials=1;
109 G4Material* theMaterials[2000];
110
111 G4cout << "Please specify A, Z"<<G4endl;
112 G4cin >> a >> iz;
113 G4Material *theMaterial = new G4Material(name="stuff", density=18.95*g/cm3, nEl=1);
114 G4Element *theElement = new G4Element (name="stuff", symbol="Z", iz, a*g/mole);
115 theMaterial->AddElement( theElement, 1 );
116 theMaterials[25] = theMaterial;
117
118 G4cout << "Please enter material number"<<G4endl;
119 G4int inputNumber = 25;
120 theMaterials[0]=theMaterials[inputNumber];
121 G4cout << "Active Material = " << theMaterials[0]->GetName()<<G4endl;
122 G4cin >> inputNumber;
123 // ----------- here all material have been defined -----------
124
125// G4Element::DumpInfo();
126// G4Material::DumpInfo();
127
128 // ----------- the following is needed for building a track...... ------------
129
130 static const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
131 G4int imat = 0;
132 G4Box* theFrame = new G4Box ("Frame",10*m, 10*m, 10*m);
133
134 G4LogicalVolume* LogicalFrame = new G4LogicalVolume(theFrame, (*theMaterialTable)[imat],
135 "LFrame", 0, 0, 0);
136
137 G4PVPlacement* PhysicalFrame = new G4PVPlacement(0,G4ThreeVector(),
138 "PFrame",LogicalFrame,0,false,0);
139 G4RotationMatrix theNull;
140 G4ThreeVector theCenter(0,0,0);
141 G4GRSVolume * theTouchable = new G4GRSVolume(PhysicalFrame, &theNull, theCenter);
142 G4TouchableHandle theTouchableHandle(theTouchable);
143 // ----------- now get all particles of interest ---------
144 G4int numberOfParticles = 1;
145 G4ParticleDefinition* theParticles[1];
146 G4cout << "Please enter the projectile A, Z"<<G4endl;
147 G4cin >> a >> z;
148 G4ParticleDefinition* theIon = G4ParticleTable::GetParticleTable()->FindIon(z, a, 0, z);
149 theParticles[0]=theIon;
150
151 //------ here all the particles are Done ----------
152 G4cout << "Done with all the particles" << G4endl;
153 G4cout << "Starting process definitions" << G4endl;
154 //--------- Processes definitions ---------
155 G4IonInelasticProcess* theProcesses[1];
156
157 G4ProcessManager* theIonProcessManager = new G4ProcessManager(theIon);
158 theIon->SetProcessManager(theIonProcessManager);
159 G4IonInelasticProcess theInelasticProcess;
160 G4BinaryLightIonReaction theIonModel;
161// G4TripathiCrossSection theXSec;
162// theInelasticProcess.GetCrossSectionDataStore()->AddDataSet(&theXSec);
163
164// theInelasticProcess.RegisterMe(theTheoModel);
165 theInelasticProcess.RegisterMe(&theIonModel);
166 theIonProcessManager->AddDiscreteProcess(&theInelasticProcess);
167 theProcesses[0] = &theInelasticProcess;
168
169 G4ForceCondition* condition = new G4ForceCondition;
170 *condition = NotForced;
171
172 G4cout << "Done with all the process definitions"<<G4endl;
173 // G4cout << "Building the CrossSection Tables. This will take a while"<<G4endl;
174
175 // ----------- define energies of interest ----------------
176
177 // ----------- needed to Build a Step, and a track -------------------
178
179 G4ParticleMomentum theDirection(1.,0.,0.);
180 G4ThreeVector aPosition(0.,0.,0.);
181 G4double aTime = 0.0;
182
183 // --------- Test the GetMeanFreePath
184
185 G4StepPoint * aPreStepPoint, * aPostStepPoint;
186 G4Step aStep;
187 aPreStepPoint= aStep.GetPreStepPoint();
188 G4double meanFreePath;
189 G4double incomingEnergy;
190 G4int k, i, l, hpw = 0;
191 // --------- Test the PostStepDoIt now, 10 events each --------------
192 G4cout << "Entering the DoIt loops!!!!!"<< G4endl;
193 G4ParticleChange* aFinalState;
194 G4cout << "Test the DoIt: please enter the number of events"<<G4endl;
195 G4int ll0;
196 G4cin >> ll0;
197// G4cout <<"Now debug the DoIt: enter the problem event number"<< G4endl;
198 G4int debugThisOne=1;
199// G4cin >> debugThisOne;
200 G4cout << "Please enter the projectile energy [AMeV]"<<G4endl;
201 G4cin >> incomingEnergy;
202 incomingEnergy*=a;
203 G4int errorOne;
204 G4cout << "Please enter the problematic event number"<<G4endl;
205 G4cin >> errorOne;
206 for (i=0; i<numberOfParticles; i++)
207 {
208 outFile << G4endl
209 << "New particle type: " << theParticles[i]->GetParticleName()
210 << " " << i << G4endl;
211 for ( G4int k=0; k<numberOfMaterials; k++)
212 {
213 outFile << G4endl << "Entering Material " << theMaterials[0]->GetName()
214 << " for particle " << theParticles[i]->GetParticleName() << G4endl;
215 LogicalFrame->SetMaterial(theMaterials[0]);
216// for( G4int j=0; j<numberOfEnergies; ++j )
217int j = 0;
218 {
219 for( G4int l=0; l<ll0; ++l )
220 {
221// incomingEnergy=theEnergies[j];
222 G4DynamicParticle* aParticle =
223 new G4DynamicParticle( theParticles[i], theDirection, incomingEnergy );
224 G4Track* aTrack = new G4Track( aParticle, aTime, aPosition );
225 aTrack->SetTouchableHandle(theTouchableHandle);
226 aStep.SetTrack( aTrack );
227 aPreStepPoint = aStep.GetPreStepPoint();
228 aPostStepPoint = aStep.GetPostStepPoint();
229 aPreStepPoint->SetTouchableHandle(theTouchableHandle);
230 aPreStepPoint->SetMaterial(theMaterials[0]);
231 aPostStepPoint->SetTouchableHandle(theTouchableHandle);
232 aPostStepPoint->SetMaterial(theMaterials[0]);
233// aStep.SetPreStepPoint(aPreStepPoint);
234// aStep.SetPostStepPoint(aPostStepPoint);
235 aTrack->SetStep(&aStep);
236 ++hpw;
237 if(hpw == 1*(hpw/1))
238 G4cerr << "FINAL EVENTCOUNTER=" << hpw
239 << " current energy: " << incomingEnergy
240 << " of particle " << aParticle->GetDefinition()->GetParticleName()
241 << " in material " << theMaterials[k]->GetName() << G4endl;
242 if(hpw == errorOne)
243 {
244 hpw --;
245 hpw ++;
246 }
247 G4cout << "Last chance before DoIt call: "
248// << theIonModel.GetNiso()
249 <<G4endl;
250 aFinalState = (G4ParticleChange*) (theProcesses[i]->PostStepDoIt( *aTrack, aStep ));
251 G4cout << "NUMBER OF SECONDARIES="<<aFinalState->GetNumberOfSecondaries();
252 G4double theFSEnergy = aFinalState->GetEnergyChange();
253 G4ThreeVector * theFSMomentum= const_cast<G4ThreeVector *>
254 (aFinalState->GetMomentumChange());
255 G4cout << "FINAL STATE = "<<theFSEnergy<<" ";
256 G4cout <<*theFSMomentum<<G4endl;
257 G4Track * second;
258 G4DynamicParticle * aSec;
259 G4int isec;
260 G4double QValue = 0;
261 G4double QValueM1 = 0;
262 G4double QValueM2 = 0;
263 if(aFinalState->GetStatusChange() == fAlive)
264 {
265 QValue += aFinalState->GetEnergyChange();
266 if(theParticles[i]->GetBaryonNumber()<0) QValue+= 2.*theParticles[i]->GetPDGMass();
267 }
268 G4double esec(0);
269 G4ThreeVector psec(0);
270 for(isec=0;isec<aFinalState->GetNumberOfSecondaries();isec++)
271 {
272 second = aFinalState->GetSecondary(isec);
273 aSec = const_cast<G4DynamicParticle *> (second->GetDynamicParticle());
274 esec+= aSec->GetTotalEnergy();
275 psec+= aSec->GetMomentum();
276 G4cout << aSec->GetTotalEnergy()<<" ";
277 G4cout << aSec->GetMomentum().x()<<" ";
278 G4cout << aSec->GetMomentum().y()<<" ";
279 G4cout << aSec->GetMomentum().z()<<" ";
280 G4cout << aSec->GetDefinition()->GetPDGEncoding()<<" ";
281 G4cout << (1-isec)*aFinalState->GetNumberOfSecondaries()<<" ";
282 G4cout << aSec->GetDefinition()->GetParticleName()<<" ";
283 G4cout << "SECONDARIES info";
284 G4cout << G4endl;
285 if(aSec->GetDefinition()->GetParticleType() == "baryon")
286 {
287 if(aSec->GetDefinition()->GetBaryonNumber() < 0)
288 {
289 QValue += aSec->GetTotalEnergy();
290 QValue += G4Neutron::Neutron()->GetPDGMass();
291 if(isec!=0) QValueM1 += aSec->GetTotalEnergy();
292 if(isec!=0) QValueM1 += aSec->GetDefinition()->GetPDGMass();
293 if(isec>1) QValueM2 += aSec->GetTotalEnergy();
294 if(isec>1) QValueM2 += aSec->GetDefinition()->GetPDGMass();
295// G4cout << "found an anti !!!" <<G4endl;
296 }
297 else
298 {
299 G4double ss = 0;
300 ss +=aSec->GetDefinition()->GetPDGMass();
301 if(aSec->GetDefinition() == G4Proton::Proton())
302 {
303 ss -=G4Proton::Proton()->GetPDGMass();
304 }
305 else
306 {
307 ss -=G4Neutron::Neutron()->GetPDGMass();
308 }
309 ss += aSec->GetKineticEnergy();
310 QValue += ss;
311 if(isec!=0) QValueM1 += ss;
312 if(isec>1) QValueM2 += ss;
313// G4cout << "found a Baryon !!!" <<G4endl;
314 }
315 }
316 else if(aSec->GetDefinition()->GetPDGEncoding() == 0)
317 {
318 QValue += aSec->GetKineticEnergy();
319// G4cout << "found a ion !!!" <<G4endl;
320 }
321 else
322 {
323 QValue += aSec->GetTotalEnergy();
324 if(isec!=0) QValueM1 += aSec->GetTotalEnergy();
325 if(isec>1) QValueM2 += aSec->GetKineticEnergy();
326// G4cout << "found a Meson !!!" <<G4endl;
327 }
328 delete second;
329 }
330 G4cout << "QVALUE = "<<QValue<<" "<<hpw<<" ";
331// G4cout <<QValueM1<<" "<<QValueM2<<" "
332 G4cout <<G4endl;
333 G4cout <<"Event Total e/p " << esec/MeV << " " << psec/MeV<< G4endl;
334 delete aParticle;
335 delete aTrack;
336 aFinalState->Clear();
337 } // event loop
338 } // energy loop
339 } // material loop
340 } // particle loop
341 G4cout << G4endl <<"Light-ION-SUCESS" <<G4endl;
342 return EXIT_SUCCESS;
343}
344
345
Note: See TracBrowser for help on using the repository browser.