Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

G4LowEnergyPolarizedComptonTest.cc @ 1350

Last change on this file since 1350 was 1350, checked in by garnier, 13 years ago
update to last version 4.9.4
File size: 22.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: G4LowEnergyPolarizedComptonTest.cc,v 1.8 2006/06/29 19:44:05 gunter Exp $
28	// GEANT4 tag $Name: geant4-09-04-ref-00 $
29	//
30	// -------------------------------------------------------------------
31	// GEANT 4 class file --- Copyright CERN 1998
32	// CERN Geneva Switzerland
33	//
34	//
35	// File name: G4ComptonScatteringTest.cc
36	//
37	// Author: Francesco Longo & Gerardo Depaola
38	//
39	// Creation date: 23 january 2001
40	//
41	// Modifications:
42	//
43	// -------------------------------------------------------------------
44
45	#include "globals.hh"
46	#include "G4ios.hh"
47	#include <fstream>
48	#include <iomanip>
49
50	#include "G4ParticleDefinition.hh"
51	#include "G4ParticleTypes.hh"
52	#include "G4ParticleTable.hh"
53	#include "G4Material.hh"
54	#include "G4MaterialTable.hh"
55	#include "G4VDiscreteProcess.hh"
56	#include "G4VProcess.hh"
57	#include "G4ProcessManager.hh"
58
59	#include "G4ComptonScattering.hh"
60	#include "G4PolarizedComptonScattering.hh"
61	#include "G4LowEnergyCompton.hh"
62	#include "G4LowEnergyPolarizedCompton.hh"
63
64	#include "G4EnergyLossTables.hh"
65	#include "G4VParticleChange.hh"
66	#include "G4ParticleChange.hh"
67	#include "G4DynamicParticle.hh"
68
69	#include "G4LowEnergyBremsstrahlung.hh"
70	#include "G4LowEnergyIonisation.hh"
71	#include "G4eIonisation.hh"
72	#include "G4MultipleScattering.hh"
73	#include "G4eIonisation.hh"
74	#include "G4eBremsstrahlung.hh"
75	#include "G4eplusAnnihilation.hh"
76
77	#include "G4Electron.hh"
78	#include "G4Positron.hh"
79	#include "G4Gamma.hh"
80
81	#include "G4GRSVolume.hh"
82	#include "G4Box.hh"
83	#include "G4PVPlacement.hh"
84	#include "G4Step.hh"
85
86	#include "G4UnitsTable.hh"
87
88	#include "CLHEP/Hist/TupleManager.h"
89	#include "CLHEP/Hist/HBookFile.h"
90	#include "CLHEP/Hist/Histogram.h"
91	#include "CLHEP/Hist/Tuple.h"
92
93
94	HepTupleManager* hbookManager;
95
96	G4int main()
97	{
98
99	// Setup
100
101	G4int nIterations = 100000;
102	G4int materialId = 3;
103
104	// G4cout.setf( ios::scientific, ios::floatfield );
105
106	// -------------------------------------------------------------------
107
108	// ---- HBOOK initialization
109
110
111	hbookManager = new HBookFile("comptontest.hbook", 58);
112	assert (hbookManager != 0);
113
114	// ---- Book a histogram and ntuples
115
116	G4cout<<"Hbook file name: "<<((HBookFile*) hbookManager)->filename()<<G4endl;
117
118
119	G4double initEnergy = 1*MeV;
120	G4double initX = 0.;
121	G4double initY = 0.;
122	G4double initZ = 1.;
123
124	G4cout << "Enter the initial particle energy E (keV)" << G4endl;
125	G4cin >> initEnergy ;
126	initEnergy = initEnergy * keV;
127	G4double limit = initEnergy/keV;
128
129	G4cout << limit << G4endl;
130
131	if (initEnergy <= 0.) G4Exception("Wrong input");
132
133	// ---- primary ntuple ------
134	HepTuple* ntuple1 = hbookManager->ntuple("Primary Ntuple");
135	assert (ntuple1 != 0);
136
137	// ---- secondary ntuple ------
138	HepTuple* ntuple2 = hbookManager->ntuple("Secondary Ntuple");
139	assert (ntuple2 != 0);
140
141	/*
142	// ---- table ntuple ------
143	HepTuple* ntuple3 = hbookManager->ntuple("Mean Free Path Ntuple");
144	assert (ntuple3 != 0);
145	*/
146
147	// ---- secondaries histos ----
148
149	HepHistogram* heETot;
150	heETot = hbookManager->histogram("Electron Total Energy", 100,0.,limit);
151	assert (heETot != 0);
152
153	HepHistogram* heP;
154	heP = hbookManager->histogram("Electron Momentum", 100,0.,limit);
155	assert (heP != 0);
156
157	HepHistogram* hgETot;
158	hgETot = hbookManager->histogram("Gamma Total Energy", 100,0.,limit);
159	assert (hgETot != 0);
160
161	HepHistogram* hgP;
162	hgP = hbookManager->histogram("Gamma Momentum", 100,0.,limit);
163	assert (hgP != 0);
164
165	HepHistogram* hgTheta;
166	hgTheta = hbookManager->histogram("Theta Scattered Gamma ", 100,0.,4.);
167	assert (hgTheta != 0);
168
169	HepHistogram* hgPhi;
170	hgPhi = hbookManager->histogram("Phi Scattered Gamma ", 100,-4.,4.);
171	assert (hgPhi != 0);
172
173	HepHistogram* hSumE;
174	hSumE = hbookManager->histogram("Total Energy", 100,0.,2*limit);
175	assert (hSumE != 0);
176
177	HepHistogram* hgRapp;
178	hgRapp = hbookManager->histogram("Energy Theta Relation", 100,0.,2.);
179	assert (hgRapp != 0);
180
181	HepHistogram* hNSec;
182	hNSec = hbookManager->histogram("Number of secondaries", 100,0.,10.);
183	assert (hNSec != 0);
184
185	HepHistogram* hDebug;
186	hDebug = hbookManager->histogram("Debug", 100,0.,limit);
187	assert (hDebug != 0);
188
189
190	//--------- Materials definition ---------
191
192	G4Material* Si = new G4Material("Silicon", 14., 28.055g/mole, 2.33g/cm3);
193	G4Material* Fe = new G4Material("Iron", 26., 55.85g/mole, 7.87g/cm3);
194	G4Material* Cu = new G4Material("Copper", 29., 63.55g/mole, 8.96g/cm3);
195	G4Material* W = new G4Material("Tungsten", 74., 183.85g/mole, 19.30g/cm3);
196	G4Material* Pb = new G4Material("Lead", 82., 207.19g/mole, 11.35g/cm3);
197	G4Element* H = new G4Element ("Hydrogen", "H", 1. , 1.01*g/mole);
198	G4Element* O = new G4Element ("Oxygen" , "O", 8. , 16.00*g/mole);
199	G4Element* C = new G4Element ("Carbon" , "C", 6. , 12.00*g/mole);
200	G4Element* Cs = new G4Element ("Cesium" , "Cs", 55. , 132.905*g/mole);
201	G4Element* I = new G4Element ("Iodide" , "I", 53. , 126.9044*g/mole);
202	G4Material* maO = new G4Material("Oxygen", 8., 16.00g/mole, 1.1g/cm3);
203
204	G4Material* csi = new G4Material ("CsI" , 4.53*g/cm3, 2);
205	csi->AddElement(Cs,1);
206	csi->AddElement(I,1);
207
208
209	// Interactive set-up
210
211	G4cout << "How many interactions? " << G4endl;
212	G4cin >> nIterations;
213
214	if (nIterations <= 0) G4Exception("Wrong input");
215
216
217
218	static const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
219
220	G4int nMaterials = G4Material::GetNumberOfMaterials();
221
222	G4cout << "Available materials are: " << G4endl;
223	for (G4int mat = 0; mat < nMaterials; mat++)
224	{
225	G4cout << mat << ") "
226	<< (*theMaterialTable)[mat]->GetName()
227	<< G4endl;
228	}
229
230	G4cout << "Which material? " << G4endl;
231	G4cin >> materialId;
232
233	G4Material* material = (*theMaterialTable)[materialId] ;
234
235	G4cout << "The selected material is: "
236	<< material->GetName()
237	<< G4endl;
238
239	G4double dimX = 1*mm;
240	G4double dimY = 1*mm;
241	G4double dimZ = 1*mm;
242
243	// Geometry
244
245	G4Box* theFrame = new G4Box ("Frame",dimX, dimY, dimZ);
246
247	G4LogicalVolume* logicalFrame = new G4LogicalVolume(theFrame,
248	(*theMaterialTable)[materialId],
249	"LFrame", 0, 0, 0);
250	logicalFrame->SetMaterial(material);
251
252	G4PVPlacement* physicalFrame = new G4PVPlacement(0,G4ThreeVector(),
253	"PFrame",logicalFrame,0,false,0);
254
255	// Particle definitions
256
257	G4ParticleDefinition* gamma = G4Gamma::GammaDefinition();
258	G4ParticleDefinition* electron = G4Electron::ElectronDefinition();
259	G4ParticleDefinition* positron = G4Positron::PositronDefinition();
260
261	gamma->SetCuts(0.1*micrometer);
262	electron->SetCuts(0.1*micrometer);
263	positron->SetCuts(0.1*micrometer);
264
265	G4Gamma::SetEnergyRange(2.5e-4MeV,1e5MeV);
266	G4Electron::SetEnergyRange(2.5e-4MeV,1e5MeV);
267	G4Positron::SetEnergyRange(2.5e-4MeV,1e5MeV);
268
269	G4cout<<"the cut in energy for gamma in: "<<
270	(*theMaterialTable)[materialId]->GetName()
271	<<" is: "<<(G4Gamma::GetCutsInEnergy()[materialId])/keV
272	<<" keV" <<G4endl;
273	G4cout<<"the cut in energy for e- in: "<<
274	(*theMaterialTable)[materialId]->GetName()
275	<<" is: "<<(G4Electron::GetCutsInEnergy()[materialId])/keV
276	<<" keV" <<G4endl;
277
278	// Processes
279
280
281	G4int processType;
282	G4cout << "LowEnergy [1] or Standard [2] Compton or Standard PolarizedCompton[3] or LowEnergyPolarizedCompton [4]" << G4endl;
283	G4cin >> processType;
284	if ( !(processType == 1 \|\| processType == 2
285	\|\| processType == 3 \|\| processType == 4))
286	{
287	G4Exception("Wrong input");
288	}
289
290	G4VDiscreteProcess* gammaProcess;
291
292
293	if (processType == 1)
294	{
295	gammaProcess = new G4LowEnergyCompton();
296	}
297	else if (processType == 2)
298	{
299	gammaProcess = new G4ComptonScattering();
300	}
301	else if (processType == 3)
302	{
303	gammaProcess = new G4PolarizedComptonScattering();
304	}
305	else if (processType == 4)
306	{
307	gammaProcess = new G4LowEnergyPolarizedCompton();
308	}
309
310	G4VProcess* theeminusMultipleScattering = new G4MultipleScattering();
311	G4VProcess* theeminusIonisation = new G4eIonisation();
312	G4VProcess* theeminusBremsstrahlung = new G4eBremsstrahlung();
313	G4VProcess* theeplusMultipleScattering = new G4MultipleScattering();
314	G4VProcess* theeplusIonisation = new G4eIonisation();
315	G4VProcess* theeplusBremsstrahlung = new G4eBremsstrahlung();
316	G4VProcess* theeplusAnnihilation = new G4eplusAnnihilation();
317
318	//----------------
319	// process manager
320	//----------------
321
322	// gamma
323
324	G4ProcessManager* gProcessManager = new G4ProcessManager(gamma);
325	gamma->SetProcessManager(gProcessManager);
326	gProcessManager->AddDiscreteProcess(gammaProcess);
327	G4ForceCondition* condition;
328
329	//electron
330
331	G4ProcessManager* eProcessManager = new G4ProcessManager(electron);
332	electron->SetProcessManager(eProcessManager);
333	eProcessManager->AddProcess(theeminusMultipleScattering);
334	eProcessManager->AddProcess(theeminusIonisation);
335	eProcessManager->AddProcess(theeminusBremsstrahlung);
336
337	//positron
338
339	G4ProcessManager* pProcessManager = new G4ProcessManager(positron);
340	positron->SetProcessManager(pProcessManager);
341	pProcessManager->AddProcess(theeplusMultipleScattering);
342	pProcessManager->AddProcess(theeplusIonisation);
343	pProcessManager->AddProcess(theeplusBremsstrahlung);
344	pProcessManager->AddProcess(theeplusAnnihilation);
345
346	//--------------
347	// set ordering
348	//--------------
349
350
351	eProcessManager->
352	SetProcessOrdering(theeminusMultipleScattering, idxAlongStep,1);
353	eProcessManager->
354	SetProcessOrdering(theeminusIonisation, idxAlongStep,2);
355
356	eProcessManager->
357	SetProcessOrdering(theeminusMultipleScattering, idxPostStep,1);
358	eProcessManager->
359	SetProcessOrdering(theeminusIonisation, idxPostStep,2);
360	eProcessManager->
361	SetProcessOrdering(theeminusBremsstrahlung, idxPostStep,3);
362
363	pProcessManager->SetProcessOrderingToFirst(theeplusAnnihilation, idxAtRest);
364	pProcessManager->
365	SetProcessOrdering(theeplusMultipleScattering, idxAlongStep,1);
366	pProcessManager->
367	SetProcessOrdering(theeplusIonisation, idxAlongStep,2);
368
369	pProcessManager->
370	SetProcessOrdering(theeplusMultipleScattering, idxPostStep,1);
371	pProcessManager->
372	SetProcessOrdering(theeplusIonisation, idxPostStep,2);
373	pProcessManager->
374	SetProcessOrdering(theeplusBremsstrahlung, idxPostStep,3);
375	pProcessManager->
376	SetProcessOrdering(theeplusAnnihilation, idxPostStep,4);
377
378	// Create a DynamicParticle
379
380	// G4double eEnergy = initEnergy*keV;
381	G4double eEnergy = initEnergy;
382
383	// G4cout << eEnergy/keV << " INIT ENERGY (keV)" << G4endl;
384
385	G4ParticleMomentum eDirection(initX,initY,initZ);
386	G4DynamicParticle dynamicGamma(G4Gamma::Gamma(),eDirection,eEnergy);
387
388	G4cout << eDirection << " Direction" << G4endl;
389
390
391	// if (processType == 3 \|\| processType == 4)
392	// {
393	G4double PolX, PolY, PolZ;
394	G4cout << "Polarization Vector" << G4endl;
395	G4cin >> PolX >> PolY >> PolZ;
396	dynamicGamma.SetPolarization(PolX, PolY, PolZ);
397	//G4cout << "polarization" << dynamicGamma.GetPolarization() << G4endl;
398	// }
399
400
401	dynamicGamma.DumpInfo();
402
403	// Track
404
405	G4ThreeVector aPosition(0.,0.,0.);
406	G4double aTime = 0. ;
407
408	G4Track* gTrack = new G4Track(&dynamicGamma,aTime,aPosition);
409
410	G4GRSVolume* touche = new G4GRSVolume(physicalFrame, NULL, aPosition);
411	gTrack->SetTouchable(touche);
412
413	// Step
414
415	G4Step* step = new G4Step();
416	step->SetTrack(gTrack);
417
418	G4StepPoint* aPoint = new G4StepPoint();
419	aPoint->SetPosition(aPosition);
420	aPoint->SetMaterial(material);
421	G4double safety = 10000.*cm;
422	aPoint->SetSafety(safety);
423	step->SetPreStepPoint(aPoint);
424
425	// Check applicability
426
427	if (! (gammaProcess->IsApplicable(*gamma)))
428	{
429	G4Exception("Not Applicable");
430	}
431	else
432	{
433	G4cout<< "applicability OK" << G4endl;
434	}
435
436	// Initialize the physics tables (in which material?)
437
438	gammaProcess->BuildPhysicsTable(*gamma);
439
440	theeminusMultipleScattering->BuildPhysicsTable(*electron);
441	theeminusIonisation->BuildPhysicsTable(*electron);
442	theeminusBremsstrahlung->BuildPhysicsTable(*electron);
443	theeplusMultipleScattering->BuildPhysicsTable(*positron);
444	theeplusIonisation->BuildPhysicsTable(*positron);
445	theeplusBremsstrahlung->BuildPhysicsTable(*positron);
446	theeplusAnnihilation->BuildPhysicsTable(*positron) ;
447
448	// G4cout<< "table OK" << endl;
449
450	/*
451
452	// Test GetMeanFreePath()
453
454	G4Material* apttoMaterial ;
455	G4String MaterialName ;
456
457	G4double minArg = 100eV,maxArg = 100GeV, argStp;
458	const G4int pntNum = 300;
459	G4double Tkin[pntNum+1];
460	G4double meanFreePath=0. ;
461
462	argStp = (std::log10(maxArg)-std::log10(minArg))/pntNum;
463
464	for(G4int d = 0; d < pntNum+1; d++)
465	{
466	Tkin[d] = std::pow(10,(std::log10(minArg) + d*argStp));
467	}
468
469	G4double sti = 1.*mm;
470	step->SetStepLength(sti);
471
472	// for ( G4int J = 0 ; J < nMaterials ; J++ )
473	// {
474	apttoMaterial = (*theMaterialTable)[materialId] ;
475	MaterialName = apttoMaterial->GetName() ;
476	logicalFrame->SetMaterial(apttoMaterial);
477
478	gTrack->SetStep(step);
479
480	G4LowEnergyCompton* gammaLowEProcess =
481	(G4LowEnergyCompton*) gammaProcess;
482	G4ComptonScattering* gammaStdProcess =
483	(G4ComptonScattering*) gammaProcess;
484
485
486	for (G4int i=0 ; i<pntNum; i++)
487	{
488	dynamicGamma.SetKineticEnergy(Tkin[i]);
489	if (processType == 1)
490	{
491	meanFreePath=gammaLowEProcess
492	->GetMeanFreePath(*gTrack, sti, condition);
493	}
494	else
495	{
496	meanFreePath=gammaStdProcess
497	->GetMeanFreePath(*gTrack, sti, condition);
498	}
499
500	ntuple3->column("kinen",Tkin[i]);
501	ntuple3->column("mfp",meanFreePath/cm);
502	ntuple3->dumpData();
503
504	// G4cout << meanFreePath/cm << G4endl;
505
506	}
507	G4cout << "Mean Free Path OK" << G4endl;
508	*/
509
510	// --------- Test the DoIt
511
512	G4cout << "DoIt in " << material->GetName() << G4endl;
513
514	dynamicGamma.SetKineticEnergy(eEnergy);
515	dynamicGamma.SetMomentumDirection(initX,initY,initZ);
516
517	for (G4int iter=0; iter<nIterations; iter++)
518	{
519
520	step->SetStepLength(1*micrometer);
521
522
523
524	G4cout << "Iteration = " << iter
525	<< " - Step Length = "
526	<< step->GetStepLength()/mm << " mm "
527	<< G4endl;
528
529	gTrack->SetStep(step);
530
531	G4StepPoint* preStep = step->GetPreStepPoint();
532	G4StepPoint* postStep = step->GetPostStepPoint();
533	G4ThreeVector prePosition = preStep->GetPosition();
534	G4ThreeVector postPosition = postStep->GetPosition();
535
536	//G4cout << prePosition << "pre step point "<< G4endl;
537	//G4cout << postPosition << "post step point "<< G4endl;
538
539	G4ThreeVector polInitial=dynamicGamma.GetPolarization();
540
541	G4cout << polInitial << " Initial Polarization" << G4endl;
542
543	G4VParticleChange* dummy;
544	dummy = gammaProcess->PostStepDoIt(gTrack, step);
545	G4ParticleChange* particleChange = (G4ParticleChange*) dummy;
546
547
548	// Primary physical quantities
549
550
551	// particleChange->DumpInfo();
552
553	G4double energyChange = particleChange->GetEnergyChange();
554
555	G4double dedx = initEnergy - energyChange ;
556
557	G4double dedxNow = dedx / (step->GetStepLength());
558
559	G4ThreeVector eChange =
560	particleChange->CalcMomentum(energyChange,
561	(*particleChange->GetMomentumChange()),
562	particleChange->GetMassChange());
563
564	G4double pxChange = eChange.x();
565	G4double pyChange = eChange.y();
566	G4double pzChange = eChange.z();
567	G4double pChange =
568	std::sqrt(pxChangepxChange + pyChangepyChange + pzChange*pzChange);
569	G4double thetaChange = eChange.theta();
570
571
572	const G4ThreeVector* momChange =particleChange->GetMomentumDirectionChange();
573
574	G4cout << (momChange->x()) << " " << (momChange->y()) << " " << (momChange->z()) << " " << G4endl;
575
576	G4cout << eChange << "newdir" << G4endl;
577	G4double phiChange = eChange.phi();
578
579	G4double xChange = particleChange->GetPositionChange()->x();
580	G4double yChange = particleChange->GetPositionChange()->y();
581	G4double zChange = particleChange->GetPositionChange()->z();
582
583	//G4cout << "Theta " << thetaChange << G4endl;
584	//G4cout << "Phi " << phiChange << G4endl;
585
586	G4cout << "---- Primary after the step ---- " << G4endl;
587
588	G4cout << "Position (x,y,z) = "
589	<< xChange << " "
590	<< yChange << " "
591	<< zChange << " "
592	<< G4endl;
593
594	G4cout << " Initial Energy " << initEnergy/keV << " keV" << G4endl;
595	G4cout << "---- Energy: " << energyChange/MeV << " MeV, "
596	<< "(px,py,pz): ("
597	<< pxChange/keV << ","
598	<< pyChange/keV << ","
599	<< pzChange/keV << ") keV"
600	<< G4endl;
601	/* G4cout << "---- Energy loss (dE) = " << dedx/keV << " keV" << G4endl;
602	G4cout << "Stopping power (dE/dx)=" << dedxNow << G4endl;
603	*/
604
605
606	G4double electronMass = 511.22*keV; // da inserire la definizione
607
608	G4double Ratio = energyChange/
609	(initEnergy/(1 + (initEnergy*(1-std::cos(thetaChange))/electronMass)));
610	// testenergy
611
612	//G4cout << Ratio << "RATIO" << G4endl;
613	//G4cout << energyChange/keV << "ENERGY (keV)" << G4endl;
614
615
616	const G4ThreeVector* polChange=particleChange->GetPolarizationChange();
617
618
619	//G4cout << pxChange/pChange << "X" << G4endl;
620	//G4cout << pyChange/pChange << "Y" << G4endl;
621	//G4cout << pzChange/pChange << "Z" << G4endl;
622
623	//G4cout << polChange->x() << "pol X" << G4endl;
624	//G4cout << polChange->y() << "pol Y" << G4endl;
625	//G4cout << polChange->z() << "pol Z" << G4endl;
626	//G4cout << polChange->mag() << "pol mag" << G4endl;
627
628
629	G4double ScalarProduct = (polChange->x())*(pxChange/pChange)+
630	(pyChange/pChange)*(polChange->y())+
631	(pzChange/pChange)*(polChange->z());
632
633	//G4cout << ScalarProduct << "scalar product" << G4endl;
634
635	hgETot->accumulate(energyChange/keV);
636	hgP->accumulate(pChange/keV);
637	hgTheta->accumulate(thetaChange);
638	hgPhi->accumulate(phiChange);
639	hgRapp->accumulate(Ratio);
640
641	// Secondaries
642
643	ntuple1->column("eprimary", initEnergy/keV);
644	ntuple1->column("energyf", energyChange/keV);
645	ntuple1->column("de", dedx/keV);
646	ntuple1->column("dedx", dedxNow/keV);
647	ntuple1->column("pxch", pxChange);
648	ntuple1->column("pych", pyChange);
649	ntuple1->column("pzch", pzChange);
650	ntuple1->column("pch", pChange);
651	ntuple1->column("polx",(polInitial.x()));
652	ntuple1->column("poly",(polInitial.y()));
653	ntuple1->column("polz",(polInitial.z()));
654	ntuple1->column("polchx",(polChange->x()));
655	ntuple1->column("polchy",(polChange->y()));
656	ntuple1->column("polchz",(polChange->z()));
657	ntuple1->column("thetach", thetaChange);
658	ntuple1->column("phich", phiChange);
659	ntuple1->dumpData();
660
661	// Secondaries physical quantities
662
663	hNSec->accumulate(particleChange->GetNumberOfSecondaries());
664	hDebug->accumulate(particleChange->GetLocalEnergyDeposit());
665
666	G4cout << " secondaries " <<
667	particleChange->GetNumberOfSecondaries() << G4endl;
668
669	G4double Etotal = 0.;
670	Etotal += energyChange;
671
672	//G4cout << " Total energy" << Etotal << G4endl;
673
674	for (G4int i = 0; i < (particleChange->GetNumberOfSecondaries()); i++)
675	{
676	// The following two items should be filled per event, not
677	// per secondary; filled here just for convenience, to avoid
678	// complicated logic to dump ntuple when there are no secondaries
679
680	G4Track* finalParticle = particleChange->GetSecondary(i) ;
681
682	G4double e = finalParticle->GetTotalEnergy();
683	G4double eKin = finalParticle->GetKineticEnergy();
684	G4double px = (finalParticle->GetMomentum()).x();
685	G4double py = (finalParticle->GetMomentum()).y();
686	G4double pz = (finalParticle->GetMomentum()).z();
687	G4double theta = (finalParticle->GetMomentum()).theta();
688	G4double p = std::sqrt(pxpx+pypy+pz*pz);
689
690	if (eKin > initEnergy)
691	{
692	G4cout << "WARNING: eKinFinal > eKinInit " << G4endl;
693	// << e
694	// << " > " initEnergy
695
696	}
697
698	G4String particleName =
699	finalParticle->GetDefinition()->GetParticleName();
700	G4cout << "==== Final "
701	<< particleName << " "
702	<< "energy: " << e/keV << " keV, "
703	<< "eKin: " << eKin/keV << " keV, "
704	<< "(px,py,pz): ("
705	<< px/keV << ","
706	<< py/keV << ","
707	<< pz/keV << ") keV "
708	<< G4endl;
709
710	// G4cout << " energia secondaria" << e << G4endl;
711
712	heETot->accumulate(eKin/keV);
713	heP->accumulate(p/keV);
714
715	Etotal += eKin;
716	//G4cout << " energia totale" << Etotal << G4endl;
717
718	G4int partType;
719	if (particleName == "e-") partType = 1;
720	else if (particleName == "e+") partType = 2;
721	else if (particleName == "gamma") partType = 3;
722
723	// Fill the secondaries ntuple
724
725	ntuple2->column("event",iter);
726	ntuple2->column("eprimary",initEnergy/keV);
727	ntuple2->column("px", px);
728	ntuple2->column("py", py);
729	ntuple2->column("pz", pz);
730	ntuple2->column("p", p);
731	ntuple2->column("e", e/keV);
732	ntuple2->column("theta", theta);
733	ntuple2->column("ekin", eKin/keV);
734	ntuple2->column("type", partType);
735
736	ntuple2->dumpData();
737
738	delete particleChange->GetSecondary(i);
739	}
740
741	// G4cout << Etotal/keV << " E total (keV) " << G4endl;
742	hSumE->accumulate(Etotal/keV);
743	particleChange->Clear();
744
745	}
746
747
748	// G4cout << "Iteration number: " << iter << G4endl;
749	hbookManager->write();
750	delete hbookManager;
751
752	delete step;
753
754	G4cout << "END OF THE MAIN PROGRAM" << G4endl;
755	}
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772

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

Context Navigation

source: trunk/source/processes/electromagnetic/lowenergy/test/G4LowEnergyPolarizedComptonTest.cc @ 1350

Download in other formats: