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

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

nvx fichiers dans CVS
File size: 25.6 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//
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
76int 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.15*kelvin,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.02*mm ; // 25*micrometer ;
417 G4double gasGap = 0.50*mm ; // 1500*micrometer ;
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.0*mm ; // 40.0*mm ;
434
435 // Preparation of mixed radiator material
436
437 foilDensity = 0.91*g/cm3 ;// CH2 1.39*g/cm3; // Mylar 0.534*g/cm3; //Li
438 gasDensity = 1.2928*mg/cm3 ; // Air 0.178*mg/cm3 ; // He 1.977*mg/cm3; // CO2
439
440 totDensity = foilDensity*foilGasRatio + 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.9382723*GeV, 0.0*MeV, 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 theta*gamma = "<<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.0*k)*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.