source: trunk/source/processes/hadronic/stopping/test/MuonMinusCaptureAtRestTest.cc@ 1331

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

nvx fichiers dans CVS
File size: 10.7 KB
RevLine 
[1199]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//
27// $Id: MuonMinusCaptureAtRestTest.cc,v 1.2 2008/12/18 13:02:42 gunter Exp $
28// GEANT4 tag $Name: geant4-09-03-cand-01 $
29//
30// Johannes Peter Wellisch, 22.Apr 1997: full test-suite coded.
31#include "G4ios.hh"
32#include <fstream>
33#include <iomanip>
34
35#include "G4Material.hh"
36
37#include "G4GRSVolume.hh"
38#include "G4ProcessManager.hh"
39#include "G4HadronInelasticProcess.hh"
40
41#include "G4MuonMinusCaptureAtRest.hh"
42
43#include "G4DynamicParticle.hh"
44#include "G4LeptonConstructor.hh"
45#include "G4BaryonConstructor.hh"
46#include "G4MesonConstructor.hh"
47#include "G4IonConstructor.hh"
48
49#include "G4Box.hh"
50#include "G4PVPlacement.hh"
51
52#include "G4Step.hh"
53#include "G4StepPoint.hh"
54
55#include "G4ExcitationHandler.hh"
56
57int main()
58{
59 G4String name, symbol;
60 G4double a, iz, z, density;
61 G4int nEl, nIso;
62
63 // constructing the particles
64
65 G4LeptonConstructor aC1;
66 G4BaryonConstructor aC2;
67 G4MesonConstructor aC3;
68 G4IonConstructor aC4;
69
70 aC1.ConstructParticle();
71 aC2.ConstructParticle();
72 aC3.ConstructParticle();
73 aC4.ConstructParticle();
74
75 G4int numberOfMaterials=1;
76 G4Material* theMaterials[2000];
77
78 G4cout << "Please specify A, Z"<<G4endl;
79 G4cin >> a >> iz;
80 G4Material *theMaterial = new G4Material(name="stuff", density=18.95*g/cm3, nEl=1);
81 G4Element *theElement = new G4Element (name="stuff", symbol="Z", iz, a*g/mole);
82 theMaterial->AddElement( theElement, 1 );
83 theMaterials[25] = theMaterial;
84
85 G4cout << "Please enter material number"<<G4endl;
86 G4int inputNumber = 25;
87 theMaterials[0]=theMaterials[inputNumber];
88 G4cout << "Active Material = " << theMaterials[0]->GetName()<<G4endl;
89 G4cin >> inputNumber;
90 // ----------- here all material have been defined -----------
91
92// G4Element::DumpInfo();
93// G4Material::DumpInfo();
94
95 // ----------- the following is needed for building a track...... ------------
96
97 static const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
98 G4int imat = 0;
99 G4Box* theFrame = new G4Box ("Frame",10*m, 10*m, 10*m);
100
101 G4LogicalVolume* LogicalFrame = new G4LogicalVolume(theFrame,
102 (*theMaterialTable)[imat],
103 "LFrame", 0, 0, 0);
104
105 G4PVPlacement* PhysicalFrame = new G4PVPlacement(0,G4ThreeVector(),
106 "PFrame",LogicalFrame,0,false,0);
107 G4RotationMatrix theNull;
108 G4ThreeVector theCenter(0,0,0);
109 G4GRSVolume * theTouchable = new G4GRSVolume(PhysicalFrame, &theNull, theCenter);
110 G4TouchableHandle theTouchableHandle(theTouchable);
111 // ----------- now get all particles of interest ---------
112 G4int numberOfParticles = 1;
113 G4ParticleDefinition* theParticles[1];
114 G4ParticleDefinition* theMuon = G4MuonMinus::MuonMinusDefinition();
115 theParticles[0]=theMuon;
116
117 //------ here all the particles are Done ----------
118 G4cout << "Done with all the particles" << G4endl;
119 G4cout << "Starting process definitions" << G4endl;
120 //--------- Processes definitions ---------
121 G4MuonMinusCaptureAtRest* theProcesses[1];
122
123 G4ProcessManager* theMuonProcessManager = new G4ProcessManager(theMuon);
124 theMuon->SetProcessManager(theMuonProcessManager);
125
126 G4MuonMinusCaptureAtRest theInelasticProcess;
127 theProcesses[0] = &theInelasticProcess;
128
129 G4ForceCondition* condition = new G4ForceCondition;
130 *condition = NotForced;
131
132 G4cout << "Done with all the process definitions"<<G4endl;
133 // G4cout << "Building the CrossSection Tables. This will take a while"<<G4endl;
134
135 // ----------- define energies of interest ----------------
136
137 // ----------- needed to Build a Step, and a track -------------------
138
139 G4ParticleMomentum theDirection(1.,0.,0.);
140 G4ThreeVector aPosition(0.,0.,0.);
141 G4double aTime = 0.0;
142
143 // --------- Test the GetMeanFreePath
144
145 G4StepPoint aStepPoint;
146 G4Step aStep;
147 aStep.SetPreStepPoint(&aStepPoint);
148 G4double meanFreePath;
149 G4double incomingEnergy;
150 G4int k, i, l, hpw = 0;
151 // --------- Test the PostStepDoIt now, 10 events each --------------
152 G4cout << "Entering the DoIt loops!!!!!"<< G4endl;
153 G4ParticleChange* aFinalState;
154 G4cout << "Test the DoIt: please enter the number of events"<<G4endl;
155 G4int ll0;
156 G4cin >> ll0;
157// G4cout <<"Now debug the DoIt: enter the problem event number"<< G4endl;
158 G4int debugThisOne=1;
159// G4cin >> debugThisOne;
160 G4cout << "Please enter the Gamma energy [MeV]"<<G4endl;
161 G4cin >> incomingEnergy;
162 G4int errorOne;
163 G4cout << "Please enter the problematic event number"<<G4endl;
164 G4cin >> errorOne;
165 for (i=0; i<numberOfParticles; i++)
166 {
167 for ( G4int k=0; k<numberOfMaterials; k++)
168 {
169 LogicalFrame->SetMaterial(theMaterials[0]);
170// for( G4int j=0; j<numberOfEnergies; ++j )
171int j = 0;
172 {
173 for( G4int l=0; l<ll0; ++l )
174 {
175// incomingEnergy=theEnergies[j];
176 G4DynamicParticle* aParticle =
177 new G4DynamicParticle( theParticles[i], theDirection, incomingEnergy );
178 G4Track* aTrack = new G4Track( aParticle, aTime, aPosition );
179 aTrack->SetTouchableHandle(theTouchableHandle);
180 aStep.SetTrack( aTrack );
181 aStepPoint.SetTouchableHandle(theTouchableHandle);
182 aStepPoint.SetMaterial(theMaterials[0]);
183 aStep.SetPreStepPoint(&aStepPoint);
184 aStep.SetPostStepPoint(&aStepPoint);
185 aTrack->SetStep(&aStep);
186 ++hpw;
187 if(hpw == 1000*(hpw/1000))
188 G4cerr << "FINAL EVENTCOUNTER=" << hpw
189 << " current energy: " << incomingEnergy
190 << " of particle " << aParticle->GetDefinition()->GetParticleName()
191 << " in material " << theMaterials[k]->GetName() << G4endl;
192 if(hpw == errorOne)
193 {
194 hpw --;
195 hpw ++;
196 }
197 G4cout << "Last chance before DoIt call: "
198// << theMuonModel.GetNiso()
199 <<G4endl;
200 aFinalState = (G4ParticleChange*) (theProcesses[i]->AtRestDoIt( *aTrack, aStep ));
201 G4cout << "NUMBER OF SECONDARIES="<<aFinalState->GetNumberOfSecondaries();
202 G4double theFSEnergy = aFinalState->GetEnergyChange();
203 G4ThreeVector * theFSMomentum= const_cast<G4ThreeVector *>
204 (aFinalState->GetMomentumChange());
205 G4cout << "FINAL STATE = "<<theFSEnergy<<" ";
206 G4cout <<*theFSMomentum<<G4endl;
207 G4Track * second;
208 G4DynamicParticle * aSec;
209 G4int isec;
210 G4double QValue = 0;
211 G4double QValueM1 = 0;
212 G4double QValueM2 = 0;
213 if(aFinalState->GetStatusChange() == fAlive)
214 {
215 QValue += aFinalState->GetEnergyChange();
216 if(theParticles[i]->GetBaryonNumber()<0) QValue+= 2.*theParticles[i]->GetPDGMass();
217 }
218 for(isec=0;isec<aFinalState->GetNumberOfSecondaries();isec++)
219 {
220 second = aFinalState->GetSecondary(isec);
221 aSec = const_cast<G4DynamicParticle *> (second->GetDynamicParticle());
222 G4cout << "SECONDARIES info";
223 G4cout << aSec->GetTotalEnergy();
224 G4cout << aSec->GetMomentum();
225 G4cout << aSec->GetDefinition()->GetPDGEncoding()<<" ";
226 G4cout << (1-isec)*aFinalState->GetNumberOfSecondaries()<<" ";
227 G4cout << aSec->GetDefinition()->GetParticleName()<<" ";
228 G4cout << G4endl;
229 if(aSec->GetDefinition()->GetParticleType() == "baryon")
230 {
231 if(aSec->GetDefinition()->GetBaryonNumber() < 0)
232 {
233 QValue += aSec->GetTotalEnergy();
234 QValue += G4Neutron::Neutron()->GetPDGMass();
235 if(isec!=0) QValueM1 += aSec->GetTotalEnergy();
236 if(isec!=0) QValueM1 += aSec->GetDefinition()->GetPDGMass();
237 if(isec>1) QValueM2 += aSec->GetTotalEnergy();
238 if(isec>1) QValueM2 += aSec->GetDefinition()->GetPDGMass();
239 G4cout << "found an anti !!!" <<G4endl;
240 }
241 else
242 {
243 G4double ss = 0;
244 ss +=aSec->GetDefinition()->GetPDGMass();
245 if(aSec->GetDefinition() == G4Proton::Proton())
246 {
247 ss -=G4Proton::Proton()->GetPDGMass();
248 }
249 else
250 {
251 ss -=G4Neutron::Neutron()->GetPDGMass();
252 }
253 ss += aSec->GetKineticEnergy();
254 QValue += ss;
255 if(isec!=0) QValueM1 += ss;
256 if(isec>1) QValueM2 += ss;
257 G4cout << "found a Baryon !!!" <<G4endl;
258 }
259 }
260 else if(aSec->GetDefinition()->GetPDGEncoding() == 0)
261 {
262 QValue += aSec->GetKineticEnergy();
263 G4cout << "found a ion !!!" <<G4endl;
264 }
265 else
266 {
267 QValue += aSec->GetTotalEnergy();
268 if(isec!=0) QValueM1 += aSec->GetTotalEnergy();
269 if(isec>1) QValueM2 += aSec->GetKineticEnergy();
270 G4cout << "found a Meson !!!" <<G4endl;
271 }
272 delete second;
273 }
274 G4cout << "QVALUE = "<<QValue<<" "<<hpw<<" ";
275// G4cout <<QValueM1<<" "<<QValueM2<<" "
276 G4cout <<G4endl;
277 delete aParticle;
278 delete aTrack;
279 aFinalState->Clear();
280 } // event loop
281 } // energy loop
282 } // material loop
283 } // particle loop
284 return EXIT_SUCCESS;
285}
286
287
Note: See TracBrowser for help on using the repository browser.