Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

testG4VXTRenergyLoss.cc @ 1199

Last change on this file since 1199 was 1199, checked in by garnier, 15 years ago
nvx fichiers dans CVS
File size: 25.6 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	//
28	//
29	//
30	//
31	// Test routine for G4VXTRenergyLoss class and inherited from it code
32	//
33	// History:
34	//
35	// 20.01.05, V. Grichine
36
37	#include "G4ios.hh"
38	#include <fstream>
39	#include <cmath>
40	#include "globals.hh"
41	#include "Randomize.hh"
42	#include "G4UnitsTable.hh"
43
44
45	#include <iomanip>
46
47	#include "G4Isotope.hh"
48	#include "G4Element.hh"
49	#include "G4Material.hh"
50	#include "G4MaterialCutsCouple.hh"
51	#include "G4Region.hh"
52	#include "G4ProductionCuts.hh"
53	#include "G4RegionStore.hh"
54	#include "G4MaterialTable.hh"
55
56	#include "G4Box.hh"
57	#include "G4LogicalVolume.hh"
58
59	#include "G4VXTRenergyLoss.hh"
60	#include "G4RegularXTRadiator.hh"
61	#include "G4TransparentRegXTRadiator.hh"
62	#include "G4GammaXTRadiator.hh"
63	#include "G4StrawTubeXTRadiator.hh"
64
65	#include "G4XTRGammaRadModel.hh"
66	#include "G4XTRRegularRadModel.hh"
67	#include "G4XTRTransparentRegRadModel.hh"
68
69	#include "G4SynchrotronRadiation.hh"
70
71	#include "G4ParticleDefinition.hh"
72	#include "G4Proton.hh"
73
74
75
76	int main()
77	{
78
79	/*
80	std::ofstream outdEdx("XTRdEdx.out", std::ios::out ) ;
81	outdEdx.setf( std::ios::scientific, std::ios::floatfield );
82
83	std::ofstream outdNdx("XTRdNdx.out", std::ios::out ) ;
84	outdNdx.setf( std::ios::scientific, std::ios::floatfield );
85
86	std::ofstream outXsc("InitXsc.out", std::ios::out ) ;
87	outXsc.setf( std::ios::scientific, std::ios::floatfield );
88
89	// std::ifstream fileRead("exp.dat", std::ios::out ) ;
90	// fileRead.setf( std::ios::scientific, std::ios::floatfield );
91
92	std::ofstream fileWrite1("mpXTR.dat", std::ios::out ) ;
93	fileWrite1.setf( std::ios::scientific, std::ios::floatfield );
94	*/
95
96
97	/////////////////////////////////////////////////////////////////
98	//
99	// Create materials
100
101
102	G4String name, symbol ; //a =mass of a mole;
103	G4double a, z ; //z =mean number of protons;
104	G4double density, foilDensity, gasDensity, totDensity ;
105	G4double fractionFoil, fractionGas ;
106	G4int nel ;
107
108	//G4int ncomponents, natoms;
109	G4int ncomponents;
110	//G4double abundance, fractionmass;
111	G4double fractionmass;
112	//G4double temperature, pressure;
113
114	/////////////////////////////////////
115	//
116	// define Elements
117
118
119	a = 1.01*g/mole;
120	G4Element* elH = new G4Element(name="Hydrogen",symbol="H" , z= 1., a);
121
122	a = 6.94*g/mole;
123	G4Element* elLi = new G4Element(name="Lithium",symbol="Li" , z= 3., a);
124
125	a = 9.01*g/mole;
126	G4Element* elBe = new G4Element(name="Berillium",symbol="Be" , z= 4., a);
127
128	a = 12.01*g/mole;
129	G4Element* elC = new G4Element(name="Carbon", symbol="C", z=6., a);
130
131	a = 14.01*g/mole;
132	G4Element* elN = new G4Element(name="Nitrogen",symbol="N" , z= 7., a);
133
134	a = 16.00*g/mole;
135	G4Element* elO = new G4Element(name="Oxygen" ,symbol="O" , z= 8., a);
136
137	a = 39.948*g/mole;
138	G4Element* elAr = new G4Element(name="Argon", symbol="Ar", z=18., a);
139
140	a = 131.29*g/mole;
141	G4Element* elXe = new G4Element(name="Xenon", symbol="Xe", z=54., a);
142
143	a = 19.00*g/mole;
144	G4Element* elF = new G4Element(name="Fluorine", symbol="F", z=9., a);
145
146	/////////////////////////////////////////////////////////////////
147	//
148	// Detector windows, electrodes
149	// Al for electrodes
150
151	density = 2.700*g/cm3;
152	a = 26.98*g/mole;
153	G4Material* Al = new G4Material(name="Aluminium", z=13., a, density);
154
155
156	////////////////////////////////////////////////////////////////////////////
157	//
158	// Materials for popular X-ray TR radiators
159	//
160
161	// TRT_CH2
162
163	density = 0.935*g/cm3;
164	G4Material* TRT_CH2 = new G4Material(name="TRT_CH2",density, nel=2);
165	TRT_CH2->AddElement(elC,1);
166	TRT_CH2->AddElement(elH,2);
167
168	// Radiator
169
170	density = 0.059*g/cm3;
171	G4Material* Radiator = new G4Material(name="Radiator",density, nel=2);
172	Radiator->AddElement(elC,1);
173	Radiator->AddElement(elH,2);
174	// Carbon Fiber
175
176	density = 0.145*g/cm3;
177	G4Material* CarbonFiber = new G4Material(name="CarbonFiber",density, nel=1);
178	CarbonFiber->AddElement(elC,1);
179
180	// Lithium
181
182	density = 0.534*g/cm3;
183	G4Material* Li = new G4Material(name="Li",density, nel=1);
184	Li->AddElement(elLi,1);
185
186	// Beryllium
187
188	density = 1.848*g/cm3;
189	G4Material* Be = new G4Material(name="Be",density, nel=1);
190	Be->AddElement(elBe,1);
191
192	// Mylar
193
194	density = 1.39*g/cm3;
195	G4Material* Mylar = new G4Material(name="Mylar", density, nel=3);
196	Mylar->AddElement(elO,2);
197	Mylar->AddElement(elC,5);
198	Mylar->AddElement(elH,4);
199
200	// Kapton (polyimide) ??? since = Mylar C5H4O2
201
202	density = 1.39*g/cm3;
203	G4Material* Kapton = new G4Material(name="Kapton", density, nel=3);
204	Kapton->AddElement(elO,2);
205	Kapton->AddElement(elC,5);
206	Kapton->AddElement(elH,4);
207
208	// Polypropelene
209
210	G4Material* CH2 = new G4Material ("CH2" , 0.91*g/cm3, 2);
211	CH2->AddElement(elH,2);
212	CH2->AddElement(elC,1);
213
214	//////////////////////////////////////////////////////////////////////////
215	//
216	// Noble gases , STP conditions
217
218	// Helium as detector gas, STP
219
220	density = 0.178*mg/cm3 ;
221	a = 4.0026*g/mole ;
222	G4Material* He = new G4Material(name="He",z=2., a, density );
223
224	// Neon as detector gas, STP
225
226	density = 0.900*mg/cm3 ;
227	a = 20.179*g/mole ;
228	G4Material* Ne = new G4Material(name="Ne",z=10., a, density );
229
230	// Argon as detector gas, STP
231
232	density = 1.7836*mg/cm3 ; // STP
233	G4Material* Argon = new G4Material(name="Argon" , density, ncomponents=1);
234	Argon->AddElement(elAr, 1);
235
236	// Krypton as detector gas, STP
237
238	density = 3.700*mg/cm3 ;
239	a = 83.80*g/mole ;
240	G4Material* Kr = new G4Material(name="Kr",z=36., a, density );
241
242	// Xenon as detector gas, STP
243
244	density = 5.858*mg/cm3 ;
245	a = 131.29*g/mole ;
246	G4Material* Xe = new G4Material(name="Xenon",z=54., a, density );
247
248	/////////////////////////////////////////////////////////////////////////////
249	//
250	// Hydrocarbones, metane and others
251
252	// Metane, STP
253
254	density = 0.7174*mg/cm3 ;
255	G4Material* metane = new G4Material(name="CH4",density,nel=2) ;
256	metane->AddElement(elC,1) ;
257	metane->AddElement(elH,4) ;
258
259	// Propane, STP
260
261	density = 2.005*mg/cm3 ;
262	G4Material* propane = new G4Material(name="C3H8",density,nel=2) ;
263	propane->AddElement(elC,3) ;
264	propane->AddElement(elH,8) ;
265
266	// iso-Butane (methylpropane), STP
267
268	density = 2.67*mg/cm3 ;
269	G4Material* isobutane = new G4Material(name="isoC4H10",density,nel=2) ;
270	isobutane->AddElement(elC,4) ;
271	isobutane->AddElement(elH,10) ;
272
273	///////////////////////////////////////////////////////////////////////////
274	//
275	// Molecular gases
276
277	// Carbon dioxide, STP
278
279	density = 1.977*mg/cm3;
280	G4Material* CO2 = new G4Material(name="CO2", density, nel=2,
281	kStateGas,273.15kelvin,1.atmosphere);
282	CO2->AddElement(elC,1);
283	CO2->AddElement(elO,2);
284
285	// Carbon dioxide, STP
286
287	density = 1.977*mg/cm3;
288	G4Material* CarbonDioxide = new G4Material(name="CO2", density, nel=2);
289	CarbonDioxide->AddElement(elC,1);
290	CarbonDioxide->AddElement(elO,2);
291
292
293	// Nitrogen, STP
294
295	density = 1.25053*mg/cm3 ; // STP
296	G4Material* Nitrogen = new G4Material(name="N2" , density, ncomponents=1);
297	Nitrogen->AddElement(elN, 2);
298
299	// Oxygen, STP
300
301	density = 1.4289*mg/cm3 ; // STP
302	G4Material* Oxygen = new G4Material(name="O2" , density, ncomponents=1);
303	Oxygen->AddElement(elO, 2);
304
305	/* *****************************
306
307	density = 1.25053*mg/cm3 ; // STP
308	a = 14.01*g/mole ; // get atomic weight !!!
309	// a = 28.016*g/mole;
310	G4Material* N2 = new G4Material(name="Nitrogen", z= 7.,a,density) ;
311
312	density = 1.25053*mg/cm3 ; // STP
313	G4Material* anotherN2 = new G4Material(name="anotherN2", density,ncomponents=2);
314	anotherN2->AddElement(elN, 1);
315	anotherN2->AddElement(elN, 1);
316
317	// air made from oxigen and nitrogen only
318
319	density = 1.290*mg/cm3; // old air from elements
320	G4Material* air = new G4Material(name="air" , density, ncomponents=2);
321	air->AddElement(elN, fractionmass=0.7);
322	air->AddElement(elO, fractionmass=0.3);
323
324	******************************************** */
325
326	// Dry Air (average composition), STP
327
328	density = 1.2928*mg/cm3 ; // STP
329	G4Material* Air = new G4Material(name="Air" , density, ncomponents=3);
330	Air->AddMaterial( Nitrogen, fractionmass = 0.7557 ) ;
331	Air->AddMaterial( Oxygen, fractionmass = 0.2315 ) ;
332	Air->AddMaterial( Argon, fractionmass = 0.0128 ) ;
333
334	////////////////////////////////////////////////////////////////////////////
335	//
336	// MWPC mixtures
337
338	// 80% Xe + 20% CO2, STP
339
340	density = 5.0818*mg/cm3 ;
341	G4Material* Xe20CO2 = new G4Material(name="Xe20CO2" , density, ncomponents=2);
342	Xe20CO2->AddMaterial( Xe, fractionmass = 0.922 ) ;
343	Xe20CO2->AddMaterial( CarbonDioxide, fractionmass = 0.078 ) ;
344
345	// 80% Kr + 20% CO2, STP
346
347	density = 3.601*mg/cm3 ;
348	G4Material* Kr20CO2 = new G4Material(name="Kr20CO2", density,
349	ncomponents=2);
350	Kr20CO2->AddMaterial( Kr, fractionmass = 0.89 ) ;
351	Kr20CO2->AddMaterial( CarbonDioxide, fractionmass = 0.11 ) ;
352
353	// Xe + 55% He + 15% CH4 ; NIM A294 (1990) 465-472; STP
354
355	density = 1.963*mg/cm3;
356	G4Material* Xe55He15CH4 = new G4Material(name="Xe55He15CH4",density,
357	ncomponents=3);
358	Xe55He15CH4->AddMaterial(Xe, 0.895);
359	Xe55He15CH4->AddMaterial(He, 0.050);
360	Xe55He15CH4->AddMaterial(metane,0.055);
361
362	// 90% Xe + 10% CH4, STP ; NIM A248 (1986) 379-388
363
364	density = 5.344*mg/cm3 ;
365	G4Material* Xe10CH4 = new G4Material(name="Xe10CH4" , density,
366	ncomponents=2);
367	Xe10CH4->AddMaterial( Xe, fractionmass = 0.987 ) ;
368	Xe10CH4->AddMaterial( metane, fractionmass = 0.013 ) ;
369
370	// 95% Xe + 5% CH4, STP ; NIM A214 (1983) 261-268
371
372	density = 5.601*mg/cm3 ;
373	G4Material* Xe5CH4 = new G4Material(name="Xe5CH4" , density,
374	ncomponents=2);
375	Xe5CH4->AddMaterial( Xe, fractionmass = 0.994 ) ;
376	Xe5CH4->AddMaterial( metane, fractionmass = 0.006 ) ;
377
378	// 80% Xe + 20% CH4, STP ; NIM A253 (1987) 235-244
379
380	density = 4.83*mg/cm3 ;
381	G4Material* Xe20CH4 = new G4Material(name="Xe20CH4" , density,
382	ncomponents=2);
383	Xe20CH4->AddMaterial( Xe, fractionmass = 0.97 ) ;
384	Xe20CH4->AddMaterial( metane, fractionmass = 0.03 ) ;
385
386	// 93% Ar + 7% CH4, STP ; NIM 107 (1973) 413-422
387
388	density = 1.709*mg/cm3 ;
389	G4Material* Ar7CH4 = new G4Material(name="Ar7CH4" , density,
390	ncomponents=2);
391	Ar7CH4->AddMaterial( Argon, fractionmass = 0.971 ) ;
392	Ar7CH4->AddMaterial( metane, fractionmass = 0.029 ) ;
393
394	// 93% Kr + 7% CH4, STP ; NIM 107 (1973) 413-422
395
396	density = 3.491*mg/cm3 ;
397	G4Material* Kr7CH4 = new G4Material(name="Kr7CH4" , density,
398	ncomponents=2);
399	Kr7CH4->AddMaterial( Kr, fractionmass = 0.986 ) ;
400	Kr7CH4->AddMaterial( metane, fractionmass = 0.014 ) ;
401
402	// 0.5(95% Xe + 5% CH4)+0.5(93% Ar + 7% CH4), STP ; NIM A214 (1983) 261-268
403
404	density = 3.655*mg/cm3 ;
405	G4Material* XeArCH4 = new G4Material(name="XeArCH4" , density,
406	ncomponents=2);
407	XeArCH4->AddMaterial( Xe5CH4, fractionmass = 0.766 ) ;
408	XeArCH4->AddMaterial( Ar7CH4, fractionmass = 0.234 ) ;
409
410
411	////////////////////////////////////////////////////////////
412	//
413	// Geometry
414
415
416	G4double foilThick = 0.02mm ; // 25micrometer ;
417	G4double gasGap = 0.50mm ; // 1500micrometer ;
418	G4double foilGasRatio = foilThick/(foilThick+gasGap) ;
419	G4int foilNumber = 120 ; // 188 ;
420	G4double detGap = 0.01*mm ;
421
422
423	G4double alphaPlate = 2.0 ;
424	G4double alphaGas = 10.0 ;
425	G4int modelNumber = 0 ;
426
427	// TR radiator envelope
428
429	G4double radThick = foilNumber*(foilThick + gasGap) - gasGap + detGap;
430
431	G4double absorberRadius = 10.*cm;
432
433	G4double absorberThickness = 15.0mm ; // 40.0mm ;
434
435	// Preparation of mixed radiator material
436
437	foilDensity = 0.91g/cm3 ;// CH2 1.39g/cm3; // Mylar 0.534*g/cm3; //Li
438	gasDensity = 1.2928mg/cm3 ; // Air 0.178mg/cm3 ; // He 1.977*mg/cm3; // CO2
439
440	totDensity = foilDensityfoilGasRatio + gasDensity(1.0-foilGasRatio) ;
441
442	fractionFoil = foilDensity*foilGasRatio/totDensity ;
443	fractionGas = gasDensity*(1.0-foilGasRatio)/totDensity ;
444
445	G4Material* radiatorMat = new G4Material(name="radiatorMat" , totDensity,
446	ncomponents=2);
447	radiatorMat->AddMaterial( CH2, fractionFoil ) ;
448	radiatorMat->AddMaterial( Air, fractionGas ) ;
449
450	// G4cout << *(G4Material::GetMaterialTable()) << G4endl;
451	// default materials of the calorimeter and TR radiator
452
453	G4Material* fRadiatorMat = radiatorMat ; // CH2 Mylar ;
454	G4Material* foilMat = CH2 ; // Li ; // CH2 ;
455	G4Material* gasMat = Air ; // He; // CO2 ;
456	G4Material* absMat = Xe20CO2 ; // He ;// CO2 ;
457
458	// fWindowMat = Mylar ;
459	// fElectrodeMat = Al ;
460
461	// AbsorberMaterial = Ar7CH4; // Xe10CH4; // Xe55He15CH4;
462
463
464	// fGapMat = Xe10CH4 ;
465
466	// WorldMaterial = Air ; // CO2 ;
467
468
469
470	///////////////////////
471
472	G4int i, j, k, nBin, numOfMaterials, iSan, nbOfElements, sanIndex, row ;
473
474	const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable() ;
475
476	numOfMaterials = theMaterialTable->size();
477
478	G4String testName;
479
480	for( k = 0; k < numOfMaterials; k++ )
481	{
482	// if((*theMaterialTable)[k]->GetName() != testName) continue ;
483
484	// outFile << "Material : " <<(*theMaterialTable)[k]->GetName() << G4endl ;
485	// G4cout <<k<<"\t"<< "Material : " <<(*theMaterialTable)[k]->GetName() << G4endl ;
486	}
487
488	// G4cout<<"Enter material name for test : "<<std::flush ;
489	// G4cin>>testName ;
490
491
492	// G4Region* regGasDet = new G4Region("VertexDetector");
493	// regGasDet->AddRootLogicalVolume(logicAbsorber);
494
495	G4ProductionCuts* cuts = new G4ProductionCuts();
496	cuts->SetProductionCut(10.*mm,"gamma");
497	cuts->SetProductionCut(1.*mm,"e-");
498	cuts->SetProductionCut(1.*mm,"e+");
499
500	// regGasDet->SetProductionCuts(cuts);
501
502	G4cout.precision(4);
503
504	// G4MaterialCutsCouple* matCC = new G4MaterialCutsCouple(
505	// (*theMaterialTable)[k], cuts);
506
507	G4cout<<"radThick = " <<radThick/mm<<" mm"<<G4endl ;
508	G4cout<<"foilNumber = " <<foilNumber<<G4endl ;
509	G4cout<<"radiatorMat = " <<radiatorMat->GetName()<<G4endl ;
510
511
512	G4Box* solidRadiator = new G4Box("Radiator",1.1*absorberRadius ,
513	1.1*absorberRadius,
514	0.5*radThick ) ;
515
516	G4LogicalVolume* logicRadiator = new G4LogicalVolume(solidRadiator,
517	radiatorMat,
518	"Radiator");
519
520
521	// const G4RegionStore* theRegionStore = G4RegionStore::GetInstance();
522	// G4Region* gas = theRegionStore->GetRegion("XTRdEdxDetector");
523
524	G4VXTRenergyLoss* processXTR;
525
526	const G4ParticleDefinition proton(
527	name, 0.9382723GeV, 0.0MeV, eplus,
528	1, +1, 0,
529	1, +1, 0,
530	"baryon", 0, +1, 2212,
531	true, -1.0, NULL,
532	false, "neucleon"
533	);
534
535	G4ParticleDefinition* theProton = G4Proton::ProtonDefinition();
536	// *proton = theProton;
537
538	// G4String fXTRModel = "transpR";
539	G4String fXTRModel = "transpM";
540
541	// G4String fXTRModel = "regR";
542	// G4String fXTRModel = "regM";
543
544	// G4String fXTRModel = "gammaR";
545	// G4String fXTRModel = "gammaM";
546
547	// G4String fXTRModel = "strawR";
548
549	if(fXTRModel == "gammaR" )
550	{
551	// G4GammaXTRadiator*
552	processXTR = new G4GammaXTRadiator(logicRadiator,
553	1000.,
554	100.,
555	foilMat,
556	gasMat,
557	foilThick,
558	gasGap,
559	foilNumber,
560	"GammaXTRadiator");
561	}
562	else if(fXTRModel == "gammaM" )
563	{
564	// G4XTRGammaRadModel*
565	processXTR = new G4XTRGammaRadModel(logicRadiator,
566	1000.,
567	100.,
568	foilMat,
569	gasMat,
570	foilThick,
571	gasGap,
572	foilNumber,
573	"GammaXTRmodel");
574	}
575	else if(fXTRModel == "strawR" )
576	{
577
578	// G4StrawTubeXTRadiator*
579	processXTR = new G4StrawTubeXTRadiator(logicRadiator,
580	foilMat,
581	gasMat,
582	0.53, // foilThick,
583	3.14159, // gasGap,
584	absMat,
585	true,
586	"strawXTRadiator");
587	}
588	else if(fXTRModel == "regR" )
589	{
590	// G4RegularXTRadiator*
591	processXTR = new G4RegularXTRadiator(logicRadiator,
592	foilMat,
593	gasMat,
594	foilThick,
595	gasGap,
596	foilNumber,
597	"RegularXTRadiator");
598	}
599	else if(fXTRModel == "transpR" )
600	{
601	// G4TransparentRegXTRadiator*
602	processXTR = new G4TransparentRegXTRadiator(logicRadiator,
603	foilMat,
604	gasMat,
605	foilThick,
606	gasGap,
607	foilNumber,
608	"TranspRegXTRadiator");
609	}
610	else if(fXTRModel == "regM" )
611	{
612	// G4XTRRegularRadModel*
613	processXTR = new G4XTRRegularRadModel(logicRadiator,
614	foilMat,
615	gasMat,
616	foilThick,
617	gasGap,
618	foilNumber,
619	"RegularXTRmodel");
620
621	}
622	else if(fXTRModel == "transpM" )
623	{
624	// G4XTRTransparentRegRadModel*
625	processXTR = new G4XTRTransparentRegRadModel(logicRadiator,
626	foilMat,
627	gasMat,
628	foilThick,
629	gasGap,
630	foilNumber,
631	"TranspRegXTRmodel");
632	}
633	else
634	{
635	G4Exception("Invalid XTR model name", "InvalidSetup",
636	FatalException, "XTR model name is out of the name list");
637	}
638	processXTR->SetVerboseLevel(1);
639	// processXTR->SetAngleRadDistr(true);
640
641	// processXTR->BuildPhysicsTable(proton);
642
643	// processXTR->SetVerboseLevel(1);
644
645	static G4int totBin = processXTR->GetTotBin();
646	nBin = totBin;
647	G4cout<<"totBin = "<<totBin<<G4endl;
648
649	// test of XTR table step do-it
650
651
652	G4double energyTR = 10*keV, cofAngle = 5.1, angle2, dNdA, xCompton, lambdaC;
653	G4double charge = 1.0;
654	G4double chargeSq = charge*charge ;
655	G4double gamma = 4.e4;
656	G4cout<<"gamma = "<<gamma<<G4endl;
657	G4cout<<"energyTR = "<<energyTR/keV<<" keV"<<G4endl;
658
659	processXTR->SetGamma(gamma);
660
661	// processXTR->GetAngleVector(energyTR,nBin);
662
663	// xCompton = processXTR->GetGasCompton(energyTR);
664
665	// lambdaC = 1./xCompton;
666
667	// G4cout<<"lambdaC = "<<lambdaC/m <<" m; for energy = "<<energyTR/keV<<" keV"<<G4endl;
668
669	// G4double dNdA = processXTR->SpectralXTRdEdx(energyTR);
670
671	// G4double angle2 = cofAngle*cofAngle/gamma/gamma;
672
673	// G4double dNdAngle = processXTR-> AngleXTRdEdx(angle2);
674	/*
675	for(i = 0; i < 40; i++ )
676	{
677	angle2 = processXTR->GetRandomAngle(energyTR,40);
678	G4cout<<"random thetagamma = "<<std::sqrt(angle2)gamma<<G4endl;
679	}
680	*/
681
682	/*
683	for(i = 0; i < 40; i++ )
684	{
685	cofAngle = 0.5*i;
686	G4double angle2 = cofAngle*cofAngle/gamma/gamma;
687	G4double dNdAngle = processXTR-> AngleXTRdEdx(angle2);
688	dNdAngle *=fine_structure_const/pi;
689	G4cout<<"cofAngle = "<<cofAngle<<"; angle = "<<cofAngle/gamma
690	<<"; dNdAngle = "<<dNdAngle<<G4endl;
691	}
692	*/
693
694
695	G4int iTkin;
696	G4cout<<"gamma = "<<gamma<<G4endl;
697
698	G4double TkinScaled = (gamma - 1.)*proton_mass_c2;
699
700	/*
701	for( iTkin = 0; iTkin < totBin; iTkin++ )
702	{
703	if(TkinScaled < processXTR->GetProtonVector()->
704	GetLowEdgeEnergy(iTkin)) break;
705	}
706
707	G4double xtrEnergy[100];
708	G4int spectrum[100];
709
710
711	for( k = 0; k < 100; k++ )
712	{
713	xtrEnergy[k] = (1.0+ 1.0k)keV;
714	spectrum[k] = 0;
715	}
716
717
718	for( i = 0; i < 10000; i++ )
719	{
720	energyTR = processXTR->GetXTRrandomEnergy(TkinScaled,iTkin);
721
722	for( k = 0; k < 100; k++ )
723	{
724	if( energyTR <= xtrEnergy[k] ) break;
725	}
726	spectrum[k] += 1;
727	}
728
729	// output to file
730
731
732	if(fXTRModel == "gammaR" )
733	{
734	std::ofstream fileWrite("gammaR.dat", std::ios::out ) ;
735	fileWrite.setf( std::ios::scientific, std::ios::floatfield );
736
737	for( k = 0; k < 41; k++ )
738	{
739	G4cout<<k<<"\t"<<xtrEnergy[k]/keV<<"\t"<<spectrum[k]<<G4endl;
740	fileWrite<<xtrEnergy[k]/keV<<"\t"<<spectrum[k]<<G4endl;
741	}
742	}
743	else if(fXTRModel == "gammaM" )
744	{
745	std::ofstream fileWrite("gammaM.dat", std::ios::out ) ;
746	fileWrite.setf( std::ios::scientific, std::ios::floatfield );
747
748	for( k = 0; k < 41; k++ )
749	{
750	G4cout<<k<<"\t"<<xtrEnergy[k]/keV<<"\t"<<spectrum[k]<<G4endl;
751	fileWrite<<xtrEnergy[k]/keV<<"\t"<<spectrum[k]<<G4endl;
752	}
753	}
754	else if(fXTRModel == "strawR" )
755	{
756	std::ofstream fileWrite("strawR.dat", std::ios::out ) ;
757	fileWrite.setf( std::ios::scientific, std::ios::floatfield );
758
759	for( k = 0; k < 41; k++ )
760	{
761	G4cout<<k<<"\t"<<xtrEnergy[k]/keV<<"\t"<<spectrum[k]<<G4endl;
762	fileWrite<<xtrEnergy[k]/keV<<"\t"<<spectrum[k]<<G4endl;
763	}
764	}
765	else if(fXTRModel == "regR" )
766	{
767	std::ofstream fileWrite("regR.dat", std::ios::out ) ;
768	fileWrite.setf( std::ios::scientific, std::ios::floatfield );
769
770	for( k = 0; k < 41; k++ )
771	{
772	G4cout<<k<<"\t"<<xtrEnergy[k]/keV<<"\t"<<spectrum[k]<<G4endl;
773	fileWrite<<xtrEnergy[k]/keV<<"\t"<<spectrum[k]<<G4endl;
774	}
775	}
776	else if(fXTRModel == "transpR" )
777	{
778	std::ofstream fileWrite("transpR.dat", std::ios::out ) ;
779	fileWrite.setf( std::ios::scientific, std::ios::floatfield );
780
781	for( k = 0; k < 41; k++ )
782	{
783	G4cout<<k<<"\t"<<xtrEnergy[k]/keV<<"\t"<<spectrum[k]<<G4endl;
784	fileWrite<<xtrEnergy[k]/keV<<"\t"<<spectrum[k]<<G4endl;
785	}
786	}
787	else if(fXTRModel == "regM" )
788	{
789	std::ofstream fileWrite("regM.dat", std::ios::out ) ;
790	fileWrite.setf( std::ios::scientific, std::ios::floatfield );
791
792	for( k = 0; k < 41; k++ )
793	{
794	G4cout<<k<<"\t"<<xtrEnergy[k]/keV<<"\t"<<spectrum[k]<<G4endl;
795	fileWrite<<xtrEnergy[k]/keV<<"\t"<<spectrum[k]<<G4endl;
796	}
797	}
798	else if(fXTRModel == "transpM" )
799	{
800	std::ofstream fileWrite("transpM.dat", std::ios::out ) ;
801	fileWrite.setf( std::ios::scientific, std::ios::floatfield );
802
803	for( k = 0; k < 41; k++ )
804	{
805	G4cout<<k<<"\t"<<xtrEnergy[k]/keV<<"\t"<<spectrum[k]<<G4endl;
806	fileWrite<<xtrEnergy[k]/keV<<"\t"<<spectrum[k]<<G4endl;
807	}
808	}
809	else
810	{
811	G4Exception("Invalid XTR model name, no output file", "InvalidSetup",
812	FatalException, "XTR model name is out of the name list");
813	}
814	*/
815
816
817
818	G4cout.precision(12);
819	G4double ksi, prob;
820	G4SynchrotronRadiation* sr = new G4SynchrotronRadiation();
821	// sr->SetRootNumber(100);
822	// ksi = 1.e-8;
823	// ksi = 0.;
824	prob = sr->GetIntProbSR( ksi);
825	G4cout<<"ksi = "<<ksi<<"; SR probability = "<<prob<<G4endl;
826	return 1 ;
827	}
828

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

Context Navigation

source: trunk/source/processes/electromagnetic/xrays/test/testG4VXTRenergyLoss.cc @ 1199

Download in other formats: