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

Last change on this file since 1199 was 1199, checked in by garnier, 15 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.