Context Navigation

source: trunk/source/processes/electromagnetic/lowenergy/test/G4PenelopeComptonTest.cc @ 1199

Last change on this file since 1199 was 1199, checked in by garnier, 15 years ago
nvx fichiers dans CVS
File size: 23.9 KB

Line
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: G4PenelopeComptonTest.cc,v 1.10 2008/03/27 13:44:04 pandola Exp $
28	// GEANT4 tag $Name: geant4-09-03-cand-01 $
29	//
30	// -------------------------------------------------------------------
31	// GEANT 4 class file --- Copyright CERN 1998
32	// CERN Geneva Switzerland
33	//
34	//
35	// File name: G4PenelopeComptonTest.cc
36	//
37	// Author: Francesco Longo
38	//
39	// Creation date: 04 january 2001
40	//
41	// Modifications: Luciano Pandola (27 november 2002)
42	// Adapted in order to test G4PenelopeCompton
43	// Minor modification in n-tuple filling
44	// Updated analysis to AIDA 3.0
45	//
46	// -------------------------------------------------------------------
47
48	#include "globals.hh"
49	#include "G4ios.hh"
50	#include <fstream>
51	#include <iomanip>
52
53	#include "G4ParticleDefinition.hh"
54	#include "G4ParticleTypes.hh"
55	#include "G4ParticleTable.hh"
56	#include "G4Material.hh"
57	#include "G4MaterialTable.hh"
58	#include "G4VDiscreteProcess.hh"
59	#include "G4VProcess.hh"
60	#include "G4ProcessManager.hh"
61	#include "G4LowEnergyCompton.hh"
62	#include "G4PenelopeCompton.hh"
63	#include "G4EnergyLossTables.hh"
64	#include "G4VParticleChange.hh"
65	#include "G4ParticleChange.hh"
66	#include "G4DynamicParticle.hh"
67	#include "G4ForceCondition.hh"
68	#include "G4RunManager.hh"
69	#include "G4RegionStore.hh"
70	#include "G4LowEnergyBremsstrahlung.hh"
71	#include "G4LowEnergyIonisation.hh"
72	#include "G4eIonisation.hh"
73	#include "G4MultipleScattering.hh"
74	#include "G4eIonisation.hh"
75	#include "G4eBremsstrahlung.hh"
76	#include "G4eplusAnnihilation.hh"
77	#include "G4StateManager.hh"
78	#include "G4ApplicationState.hh"
79	#include "G4ComptonScattering52.hh"
80	#include "G4PhotoElectricEffect.hh"
81
82	#include "G4Electron.hh"
83	#include "G4Positron.hh"
84	#include "G4Gamma.hh"
85
86	#include "G4GRSVolume.hh"
87	#include "G4Box.hh"
88	#include "G4PVPlacement.hh"
89	#include "G4Step.hh"
90	#include "G4ProductionCutsTable.hh"
91	#include "G4MaterialCutsCouple.hh"
92
93	#include "G4UnitsTable.hh"
94	#include "AIDA/IManagedObject.h"
95
96	#include <memory>
97	#include "AIDA/IAnalysisFactory.h"
98	#include "AIDA/ITreeFactory.h"
99	#include "AIDA/ITree.h"
100	#include "AIDA/IHistogramFactory.h"
101	#include "AIDA/IHistogram1D.h"
102	#include "AIDA/IHistogram2D.h"
103	#include "AIDA/IHistogram3D.h"
104	#include "AIDA/ITupleFactory.h"
105	#include "AIDA/ITuple.h"
106
107
108	G4int main()
109	{
110
111	// Setup
112
113	G4int nIterations = 100000;
114	G4int materialId = 3;
115
116	//G4cout.setf(std::ios::scientific,std::ios::floatfield );
117
118	// -------------------------------------------------------------------
119
120	// ---- HBOOK initialization
121
122	G4String fileName;
123	G4cout << "Filename?" << G4endl;
124	G4cin >> fileName;
125
126	//build up the factories
127	AIDA::IAnalysisFactory *af = AIDA_createAnalysisFactory();
128	G4cout << "Analysis factory created" << G4endl;
129	//parameters for the TreeFactory
130	G4bool fileExists = true;
131	G4bool readOnly = false;
132	AIDA::ITreeFactory *tf = af->createTreeFactory();
133	AIDA::ITree *tree = tf->create(fileName,"hbook",readOnly,fileExists);
134	G4cout << "Tree store : " << tree->storeName() << G4endl;
135	G4cout << " Booked Hbook File " << G4endl;
136
137	//AIDA::IHistogramFactory hf = af->createHistogramFactory(tree);
138	AIDA::ITupleFactory tpf = af->createTupleFactory(tree );
139
140	// ---- primary ntuple ------
141	AIDA::ITuple* ntuple1 = tpf->create("1",
142	"Primary Ntuple","double eprimary,energyf,de,eloc,dedx,pxch,pych,pzch,pch,thetach");
143
144	// ---- secondary ntuple ------
145	AIDA::ITuple* ntuple2 = tpf->create("2",
146	"Secondary Ntuple","double eprimary,p_el,e_el,theta_el,costheta_el,ekin_el,ekin_rel");
147
148	// ---- table ntuple ------
149	AIDA::ITuple* ntuple3 = tpf->create("3","Mean Free Path Ntuple","double kinen,mfp");
150
151	// ---- fluorescence ntuple -------
152	AIDA::ITuple* ntuple4 = tpf->create("4","Fluorescence Ntuple","double eprimary,px,py,pz,e,partType");
153
154	//--------- Materials definition ---------
155	G4Element* SiEl = new G4Element ("SiEl","Si",14.,28.055*g/mole);
156	G4Element* FeEl = new G4Element ("FeEl","Fe",26.,55.58*g/mole);
157	G4Element* CuEl = new G4Element ("CuEl","Cu",29.,63.55*g/mole);
158	G4Element* WEl = new G4Element ("WEl","W",74.,183.85*g/mole);
159	G4Element* PbEl = new G4Element ("PbEl","Pb",82.,207.19*g/mole);
160	G4Element* UEl = new G4Element ("UEl","U",92.,238.03*g/mole);
161	G4Element* H = new G4Element ("Hydrogen", "H", 1. , 1.01*g/mole);
162	G4Element* O = new G4Element ("Oxygen" , "O", 8. , 16.00*g/mole);
163	G4Element* C = new G4Element ("Carbon" , "C", 6. , 12.00*g/mole);
164	G4Element* Cs = new G4Element ("Cesium" , "Cs", 55. , 132.905*g/mole);
165	G4Element* I = new G4Element ("Iodine" , "I", 53. , 126.9044*g/mole);
166	G4Element* N = new G4Element ("Nitrogen","N",7.,14.0*g/mole);
167	G4Element* AuEl = new G4Element ("Gold","Au",79.,196.97*g/mole);
168	G4Element* AlEl = new G4Element ("Aluminum","Al",13,26.981*g/mole);
169
170	G4Material* Si = new G4Material("Silicon", 2.33*g/cm3,1);
171	Si->AddElement(SiEl,1);
172	G4Material* Fe = new G4Material("Iron",7.87*g/cm3,1);
173	Fe->AddElement(FeEl,1);
174	G4Material* Cu = new G4Material("Copper", 8.96*g/cm3,1);
175	Cu->AddElement(CuEl,1);
176	G4Material* W = new G4Material("Tungsten",19.30*g/cm3,1);
177	W->AddElement(WEl,1);
178	G4Material* Pb = new G4Material("Lead",11.35*g/cm3,1);
179	Pb->AddElement(PbEl,1);
180	G4Material* U = new G4Material("Uranium",18.95*g/cm3,1);
181	U->AddElement(UEl,1);
182	G4Material* maO = new G4Material("Oxygen", 1.1*g/cm3,1);
183	maO->AddElement(O,2);
184	G4Material* water = new G4Material ("Water" , 1.*g/cm3, 2);
185	water->AddElement(H,2);
186	water->AddElement(O,1);
187	G4Material* ethane = new G4Material ("Ethane" , 0.4241*g/cm3, 2);
188	ethane->AddElement(H,6);
189	ethane->AddElement(C,2);
190	G4Material* csi = new G4Material ("CsI" , 4.53*g/cm3, 2);
191	csi->AddElement(Cs,1);
192	csi->AddElement(I,1);
193	G4Material* Air = new G4Material
194	("Air", 1.2929kg/m3, 2, kStateGas, 300.00kelvin, 1.0*atmosphere);
195	Air->AddElement(N,0.8);
196	Air->AddElement(O,0.2);
197	G4Material* Gold = new G4Material("Gold",19.3*g/cm3,1);
198	Gold->AddElement(AuEl,1);
199	G4Material* Al = new G4Material("Aluminum",2.70*g/cm3,1);
200	Al->AddElement(AlEl,1);
201
202
203	// Interactive set-up
204
205	G4cout << "How many interactions? " << G4endl;
206	G4cin >> nIterations;
207
208	if (nIterations <= 0) G4Exception("Wrong input");
209
210	G4double initEnergy = 10*keV;
211	G4double initX = 0.;
212	G4double initY = 0.;
213	G4double initZ = 1.;
214
215	G4cout << "Enter the initial particle energy E (MeV)" << G4endl;
216	G4cin >> initEnergy ;
217
218	initEnergy = initEnergy*MeV;
219
220	if (initEnergy <= 0.) G4Exception("Wrong input");
221
222	static const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
223
224	G4int nMaterials = G4Material::GetNumberOfMaterials();
225
226	G4cout << "Available materials are: " << G4endl;
227	for (G4int mat = 0; mat < nMaterials; mat++)
228	{
229	G4cout << mat << ") "
230	<< (*theMaterialTable)[mat]->GetName()
231	<< G4endl;
232	}
233
234	G4cout << "Which material? " << G4endl;
235	G4cin >> materialId;
236
237	G4Material* material = (*theMaterialTable)[materialId] ;
238
239	G4cout << "The selected material is: "
240	<< material->GetName()
241	<< G4endl;
242
243	G4double dimX = 1*mm;
244	G4double dimY = 1*mm;
245	G4double dimZ = 1*mm;
246
247	// Geometry
248
249	G4Box* theFrame = new G4Box ("Frame",dimX, dimY, dimZ);
250
251	G4LogicalVolume* logicalFrame = new G4LogicalVolume(theFrame,
252	(*theMaterialTable)[materialId],
253	"LFrame", 0, 0, 0);
254	logicalFrame->SetMaterial(material);
255
256	G4PVPlacement* physicalFrame = new G4PVPlacement(0,G4ThreeVector(),
257	"PFrame",logicalFrame,0,false,0);
258
259	G4RunManager* rm = new G4RunManager();
260
261	rm->GeometryHasBeenModified();
262	rm->DefineWorldVolume(physicalFrame);
263
264	G4cout << "[OK] World is defined " << G4endl;
265	G4StateManager::GetStateManager()->SetNewState(G4State_Idle);
266	rm->DumpRegion("DefaultRegionForTheWorld"); //this forces the region update!
267
268	// Particle definitions
269
270	G4ParticleDefinition* gamma = G4Gamma::GammaDefinition();
271	G4ParticleDefinition* electron = G4Electron::ElectronDefinition();
272	G4ParticleDefinition* positron = G4Positron::PositronDefinition();
273
274	G4ProductionCutsTable* cutsTable = G4ProductionCutsTable::GetProductionCutsTable();
275	G4ProductionCuts* cuts = cutsTable->GetDefaultProductionCuts();
276	G4double cutG=1*nanometer;
277	G4double cutE=1*nanometer;
278	cuts->SetProductionCut(cutG, 0); //gammas
279	cuts->SetProductionCut(cutE, 1); //electrons
280	cuts->SetProductionCut(cutE, 2); //positrons
281	G4cout << "Cuts are defined " << G4endl;
282
283	//G4Gamma::SetEnergyRange(2.5e-4MeV,1e5MeV);
284	//G4Electron::SetEnergyRange(2.5e-4MeV,1e5MeV);
285	//G4Positron::SetEnergyRange(2.5e-4MeV,1e5MeV);
286
287	G4MaterialCutsCouple* theCouple = new G4MaterialCutsCouple(material,cuts);
288
289	logicalFrame->SetMaterialCutsCouple(theCouple);
290
291	//G4cout << "La coppia: " << theCouple << G4endl;
292	G4cout << "Recalculation needed: " << theCouple->IsRecalcNeeded() << G4endl;
293
294	cutsTable->UpdateCoupleTable(physicalFrame);
295	cutsTable->PhysicsTableUpdated();
296	cutsTable->DumpCouples();
297
298
299	// Processes
300
301	G4int processType;
302	G4cout << "LowEnergy[1], Penelope [2] or Standard-5.2 [3] Compton?" << G4endl;
303	G4cin >> processType;
304	if ( !(processType == 1 \|\| processType == 2 \|\| processType == 3))
305	{
306	G4Exception("Wrong input");
307	}
308
309	G4VDiscreteProcess* gammaProcess = NULL;
310
311	if (processType == 1)
312	{
313	gammaProcess = new G4LowEnergyCompton();
314	G4cout << "The selected model is Low Energy" << G4endl;
315	}
316	else if (processType == 2)
317	{
318	G4PenelopeCompton* thePenelopeCompton = new G4PenelopeCompton();
319	//thePenelopeCompton->SetUseAtomicDeexcitation(false);
320	gammaProcess = thePenelopeCompton;
321	G4cout << "The selected model is Penelope" << G4endl;
322	//G4cout << "Fluorescence is disabled" << G4endl;
323	}
324	else if (processType == 3)
325	{
326	gammaProcess = new G4ComptonScattering52();
327	G4cout << "The selected model is Standard" << G4endl;
328	}
329
330	G4VProcess* theeminusMultipleScattering = new G4MultipleScattering();
331	G4VProcess* theeminusIonisation = new G4eIonisation();
332	G4VProcess* theeminusBremsstrahlung = new G4eBremsstrahlung();
333	G4VProcess* theeplusMultipleScattering = new G4MultipleScattering();
334	G4VProcess* theeplusIonisation = new G4eIonisation();
335	G4VProcess* theeplusBremsstrahlung = new G4eBremsstrahlung();
336	G4VProcess* theeplusAnnihilation = new G4eplusAnnihilation();
337
338	//----------------
339	// process manager
340	//----------------
341
342	// gamma
343
344	G4ProcessManager* gProcessManager = new G4ProcessManager(gamma);
345	gamma->SetProcessManager(gProcessManager);
346	gProcessManager->AddProcess(gammaProcess);
347	G4ForceCondition* condition; //l'ho fissata a zero! E' onesto??
348
349	//electron
350
351	G4ProcessManager* eProcessManager = new G4ProcessManager(electron);
352	electron->SetProcessManager(eProcessManager);
353	eProcessManager->AddProcess(theeminusMultipleScattering);
354	eProcessManager->AddProcess(theeminusIonisation);
355	eProcessManager->AddProcess(theeminusBremsstrahlung);
356
357	//positron
358
359	G4ProcessManager* pProcessManager = new G4ProcessManager(positron);
360	positron->SetProcessManager(pProcessManager);
361	pProcessManager->AddProcess(theeplusMultipleScattering);
362	pProcessManager->AddProcess(theeplusIonisation);
363	pProcessManager->AddProcess(theeplusBremsstrahlung);
364	pProcessManager->AddProcess(theeplusAnnihilation);
365
366	//--------------
367	// set ordering
368	//--------------
369
370
371	eProcessManager->
372	SetProcessOrdering(theeminusMultipleScattering, idxAlongStep,1);
373	eProcessManager->
374	SetProcessOrdering(theeminusIonisation, idxAlongStep,2);
375
376	eProcessManager->
377	SetProcessOrdering(theeminusMultipleScattering, idxPostStep,1);
378	eProcessManager->
379	SetProcessOrdering(theeminusIonisation, idxPostStep,2);
380	eProcessManager->
381	SetProcessOrdering(theeminusBremsstrahlung, idxPostStep,3);
382
383
384
385	pProcessManager->SetProcessOrderingToFirst(theeplusAnnihilation, idxAtRest);
386	pProcessManager->
387	SetProcessOrdering(theeplusMultipleScattering, idxAlongStep,1);
388	pProcessManager->
389	SetProcessOrdering(theeplusIonisation, idxAlongStep,2);
390
391	pProcessManager->
392	SetProcessOrdering(theeplusMultipleScattering, idxPostStep,1);
393	pProcessManager->
394	SetProcessOrdering(theeplusIonisation, idxPostStep,2);
395	pProcessManager->
396	SetProcessOrdering(theeplusBremsstrahlung, idxPostStep,3);
397	pProcessManager->
398	SetProcessOrdering(theeplusAnnihilation, idxPostStep,4);
399
400	// G4LowEnergyIonisation IonisationProcess;
401	// eProcessManager->AddProcess(&IonisationProcess);
402	// eProcessManager->SetProcessOrdering(&IonisationProcess,idxAlongStep,1);
403	// eProcessManager->SetProcessOrdering(&IonisationProcess,idxPostStep, 1);
404
405	// G4LowEnergyBremsstrahlung BremstrahlungProcess;
406	// eProcessManager->AddProcess(&BremstrahlungProcess);
407	// eProcessManager->SetProcessOrdering(&BremstrahlungProcess,idxAlongStep,1);
408	// eProcessManager->SetProcessOrdering(&BremstrahlungProcess,idxPostStep, 1);
409
410	// G4eIonisation IonisationPlusProcess;
411	// pPositronProcessManager->AddProcess(&IonisationPlusProcess);
412	// pProcessManager->
413	// SetProcessOrdering(&IonisationPlusProcess,idxAlongStep,1);
414	// pProcessManager->SetProcessOrdering(&IonisationPlusProcess,idxPostStep,1);
415
416
417
418	// Create a DynamicParticle
419
420	G4double eEnergy = initEnergy*MeV;
421	G4ParticleMomentum eDirection(initX,initY,initZ);
422	G4DynamicParticle dynamicGamma(G4Gamma::Gamma(),eDirection,eEnergy);
423
424	//dynamicGamma.DumpInfo(0);
425
426	// Track
427
428	G4ThreeVector aPosition(0.,0.,0.);
429	G4double aTime = 0. ;
430
431	G4Track* gTrack = new G4Track(&dynamicGamma,aTime,aPosition);
432
433	G4GRSVolume* touche = new G4GRSVolume(physicalFrame, NULL, aPosition);
434	//gTrack->SetTouchableHandle(touche); //verificare!!!!!!!!!!!!
435
436
437	// Step
438
439	G4Step* step = new G4Step();
440	step->SetTrack(gTrack);
441
442	G4StepPoint* aPoint = new G4StepPoint();
443	aPoint->SetPosition(aPosition);
444	aPoint->SetMaterial(material);
445	aPoint->SetMaterialCutsCouple(theCouple);
446	G4double safety = 10000.*cm;
447	aPoint->SetSafety(safety);
448	step->SetPreStepPoint(aPoint);
449
450	// PostStep point
451	G4StepPoint* newPoint = new G4StepPoint();
452	G4ThreeVector newPosition(0.,0.,1.*mm);
453	newPoint->SetPosition(newPosition);
454	newPoint->SetMaterial(material);
455	newPoint->SetMaterialCutsCouple(theCouple);
456	newPoint->SetSafety(safety);
457	step->SetPostStepPoint(newPoint);
458
459	// Check applicability
460
461	if (! (gammaProcess->IsApplicable(*gamma)))
462	{
463	G4Exception("Not Applicable");
464	}
465	else
466	{
467	G4cout<< "applicability OK" << G4endl;
468	}
469
470	// Initialize the physics tables (in which material?)
471	gammaProcess->BuildPhysicsTable(*gamma);
472
473	/*
474	theeminusMultipleScattering->BuildPhysicsTable(*electron);
475	theeminusIonisation->BuildPhysicsTable(*electron);
476	theeminusBremsstrahlung->BuildPhysicsTable(*electron);
477	theeplusMultipleScattering->BuildPhysicsTable(*positron);
478	theeplusIonisation->BuildPhysicsTable(*positron);
479	theeplusBremsstrahlung->BuildPhysicsTable(*positron);
480	theeplusAnnihilation->BuildPhysicsTable(*positron) ;
481	*/
482
483	G4cout<< "table OK" << G4endl;
484
485	G4Material* apttoMaterial ;
486	G4String MaterialName ;
487
488	G4double minArg = 100eV,maxArg = 100GeV, argStp;
489	const G4int pntNum = 300;
490	G4double Tkin[pntNum+1];
491	G4double meanFreePath=0. ;
492
493	argStp = (std::log10(maxArg)-std::log10(minArg))/pntNum;
494
495	for(G4int d = 0; d < pntNum+1; d++)
496	{
497	Tkin[d] = std::pow(10,(std::log10(minArg) + d*argStp));
498	}
499
500	G4double sti = 1.*mm;
501	step->SetStepLength(sti);
502
503	apttoMaterial = (*theMaterialTable)[materialId] ;
504	MaterialName = apttoMaterial->GetName() ;
505	logicalFrame->SetMaterial(apttoMaterial);
506
507	gTrack->SetStep(step);
508
509	G4LowEnergyCompton* gammaLowEProcess2 = (G4LowEnergyCompton*) gammaProcess;
510	G4PenelopeCompton* gammaLowEProcess = (G4PenelopeCompton*) gammaProcess;
511
512	for (G4int i=0 ; i<pntNum; i++)
513	{
514	dynamicGamma.SetKineticEnergy(Tkin[i]);
515	if (processType == 1)
516	{
517	meanFreePath=gammaLowEProcess2
518	->DumpMeanFreePath(*gTrack, sti, condition);
519	}
520	else if (processType == 2)
521	{
522	meanFreePath=gammaLowEProcess
523	->DumpMeanFreePath(*gTrack, sti, condition);
524	}
525
526	ntuple3->fill(ntuple3->findColumn("kinen"),std::log10(Tkin[i]));
527	ntuple3->fill(ntuple3->findColumn("mfp"),std::log10(meanFreePath/cm));
528	ntuple3->addRow();
529	}
530	G4cout << "Mean Free Path OK" << G4endl;
531
532	// --------- Test the DoIt
533
534	G4cout << "DoIt in " << material->GetName() << G4endl;
535
536
537	dynamicGamma.SetKineticEnergy(eEnergy);
538	G4int iter;
539	for (iter=0; iter<nIterations; iter++)
540	{
541
542	step->SetStepLength(1*nanometer);
543
544	G4cout << "Iteration = " << iter
545	<< " - Step Length = "
546	<< step->GetStepLength()/mm << " mm "
547	<< G4endl;
548
549
550	gTrack->SetStep(step);
551
552	G4ParticleChange* particleChange = (G4ParticleChange) gammaProcess->PostStepDoIt(gTrack,*step);
553
554	// Primary physical quantities
555
556	G4double energyChange = particleChange->GetEnergy();
557
558	G4double dedx = initEnergy - energyChange ;
559	G4double dedxNow = dedx / (step->GetStepLength());
560
561	G4ThreeVector eChange =
562	particleChange->CalcMomentum(energyChange,
563	(*particleChange->GetMomentumDirection()),
564	particleChange->GetMass());
565
566	G4double pxChange = eChange.x();
567	G4double pyChange = eChange.y();
568	G4double pzChange = eChange.z();
569	G4double pChange =
570	std::sqrt(pxChangepxChange + pyChangepyChange + pzChange*pzChange);
571	G4double thetaChange = eChange.theta();
572	thetaChange = thetaChange/deg; //conversion in degrees
573	G4double xChange = particleChange->GetPosition()->x();
574	G4double yChange = particleChange->GetPosition()->y();
575	G4double zChange = particleChange->GetPosition()->z();
576	G4double localEnergyDeposit = particleChange->GetLocalEnergyDeposit();
577
578	//G4double thetaChange = particleChange->GetPositionChange()->theta();
579
580	G4cout << "---- Primary after the step ---- " << G4endl;
581
582	G4cout << "Position (x,y,z) = "
583	<< xChange << " "
584	<< yChange << " "
585	<< zChange << " "
586	<< G4endl;
587
588	G4cout << "---- Energy: " << energyChange/MeV << " MeV, "
589	<< "(px,py,pz): ("
590	<< pxChange/MeV << ","
591	<< pyChange/MeV << ","
592	<< pzChange/MeV << ") MeV"
593	<< G4endl;
594
595	G4cout << "---- Energy loss (dE) = " << dedx/keV << " keV" << G4endl;
596	// G4cout << "Stopping power (dE/dx)=" << dedxNow << G4endl;
597
598	ntuple1->fill(ntuple1->findColumn("eprimary"),initEnergy/MeV);
599	ntuple1->fill(ntuple1->findColumn("energyf"),energyChange/initEnergy);
600	ntuple1->fill(ntuple1->findColumn("de"),dedx/MeV);
601	ntuple1->fill(ntuple1->findColumn("eloc"),localEnergyDeposit/MeV);
602	ntuple1->fill(ntuple1->findColumn("dedx"),dedxNow/(MeV/cm));
603	ntuple1->fill(ntuple1->findColumn("pxch"),pxChange/initEnergy);
604	ntuple1->fill(ntuple1->findColumn("pych"),pyChange/initEnergy);
605	ntuple1->fill(ntuple1->findColumn("pzch"),pzChange/initEnergy);
606	ntuple1->fill(ntuple1->findColumn("pch"),pChange/initEnergy);
607	ntuple1->fill(ntuple1->findColumn("thetach"),thetaChange);
608	ntuple1->addRow();
609
610	// Secondaries physical quantities
611
612	// Secondaries
613	G4cout << " Secondaries " <<
614	particleChange->GetNumberOfSecondaries() << G4endl;
615
616	G4double px_el=0,py_el=0,pz_el=0,p_el=0,e_el=0,theta_el=0,eKin_el=0;
617	G4double costheta_el = 0;
618
619	for (G4int i = 0; i < (particleChange->GetNumberOfSecondaries()); i++)
620	{
621	// The following two items should be filled per event, not
622	// per secondary; filled here just for convenience, to avoid
623	// complicated logic to dump ntuple when there are no secondaries
624
625	G4Track* finalParticle = particleChange->GetSecondary(i) ;
626
627	G4double e = finalParticle->GetTotalEnergy();
628	G4double eKin = finalParticle->GetKineticEnergy();
629	G4double px = (finalParticle->GetMomentum()).x();
630	G4double py = (finalParticle->GetMomentum()).y();
631	G4double pz = (finalParticle->GetMomentum()).z();
632	G4double theta = (finalParticle->GetMomentum()).theta();
633	G4double cosTheta = (finalParticle->GetMomentum()).cosTheta();
634	G4double p = std::sqrt(pxpx+pypy+pz*pz);
635	theta = theta/deg; //conversion in degrees
636	if (eKin > initEnergy)
637	{
638	G4cout << "WARNING: eFinal > eInit " << G4endl;
639	}
640
641	G4String particleName =
642	finalParticle->GetDefinition()->GetParticleName();
643	G4cout << "==== Final "
644	<< particleName << " "
645	<< "energy: " << e/MeV << " MeV, "
646	<< "eKin: " << eKin/MeV << " MeV, "
647	<< "(px,py,pz): ("
648	<< px/MeV << ","
649	<< py/MeV << ","
650	<< pz/MeV << ") MeV "
651	<< G4endl;
652
653	G4int partType=-1;
654	if (particleName == "e-") {
655	partType = 1;
656	px_el=px;
657	py_el=py;
658	pz_el=pz;
659	p_el=p;
660	e_el=e;
661	theta_el=theta;
662	costheta_el = cosTheta;
663	eKin_el=eKin;
664	}
665	else if (particleName == "gamma") partType = 2;
666
667
668	delete particleChange->GetSecondary(i);
669
670
671	// Fill the secondaries ntuple
672
673	// Normalize all to the energy of primary
674	// for gammas initEnergy=initP
675	if (i==0) {
676	ntuple2->fill(ntuple2->findColumn("eprimary"),initEnergy);
677	ntuple2->fill(ntuple2->findColumn("p_el"),p_el/initEnergy);
678	ntuple2->fill(ntuple2->findColumn("e_el"),e_el/(initEnergy+electron_mass_c2));
679	ntuple2->fill(ntuple2->findColumn("theta_el"),theta_el);
680	ntuple2->fill(ntuple2->findColumn("costheta_el"),costheta_el);
681	ntuple2->fill(ntuple2->findColumn("ekin_el"),eKin_el);
682	ntuple2->fill(ntuple2->findColumn("ekin_rel"),eKin_el/initEnergy);
683	ntuple2->addRow();
684	}
685	else //fluorescence
686	{
687	ntuple4->fill(ntuple4->findColumn("eprimary"),initEnergy);
688	ntuple4->fill(ntuple4->findColumn("px"),px/keV);
689	ntuple4->fill(ntuple4->findColumn("py"),py/keV);
690	ntuple4->fill(ntuple4->findColumn("pz"),pz/keV);
691	ntuple4->fill(ntuple4->findColumn("e"),eKin/keV);
692	ntuple4->fill(ntuple4->findColumn("partType"),(G4double) partType);
693	ntuple4->addRow();
694	}
695
696	}
697	particleChange->Clear();
698
699	}
700
701
702
703	G4cout << "Iteration number: " << iter << G4endl;
704
705	G4cout << "Committing.............." << G4endl;
706	tree->commit();
707	G4cout << "Closing the tree........" << G4endl;
708	tree->close();
709
710	delete step;
711	delete tpf;
712	delete tree;
713	//delete af;
714
715
716	G4cout << "END OF THE MAIN PROGRAM" << G4endl;
717	return 0;
718	}
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735

Note: See TracBrowser for help on using the repository browser.

Download in other formats: