]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TFluka/TFlukaConfigOption.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Warning corrected.
[u/mrichter/AliRoot.git] / TFluka / TFlukaConfigOption.cxx
CommitLineData
1df5fa54 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id$*/
17
18#include "TFlukaConfigOption.h"
fb2cbbec 19#include "TFlukaMCGeometry.h"
20#include "TFluka.h"
21#include "TFlukaCerenkov.h"
22
23#include <TString.h>
f79639f2 24#include <TList.h>
fb2cbbec 25#include <TObjArray.h>
26#include <TVirtualMC.h>
27#include <TGeoMaterial.h>
f79639f2 28#include <TGeoMedium.h>
29#include <TGeoManager.h>
30#include <TGeoMedium.h>
fb2cbbec 31
32Float_t TFlukaConfigOption::fgMatMin(-1.);
33Float_t TFlukaConfigOption::fgMatMax(-1.);
34FILE* TFlukaConfigOption::fgFile(0x0);
35TFlukaMCGeometry* TFlukaConfigOption::fgGeom(0x0);
36
37Double_t TFlukaConfigOption::fgDCutValue[11];
38Int_t TFlukaConfigOption::fgDProcessFlag[15];
39
40
1df5fa54 41ClassImp(TFlukaConfigOption)
42
43
44TFlukaConfigOption::TFlukaConfigOption()
45{
46 // Default constructor
fb2cbbec 47 fMedium = -1;
48 fCMatMin = -1.;
49 fCMatMax = -1.;
50 Int_t i;
51 for (i = 0; i < 11; i++) fCutValue[i] = -1.;
52 for (i = 0; i < 15; i++) fProcessFlag[i] = -1;
1df5fa54 53}
54
55
fb2cbbec 56TFlukaConfigOption::TFlukaConfigOption(Int_t medium)
1df5fa54 57{
58 // Constructor
fb2cbbec 59 fMedium = medium;
60 fCMatMin = -1.;
61 fCMatMax = -1.;
62 Int_t i;
63 for (i = 0; i < 11; i++) fCutValue[i] = -1.;
64 for (i = 0; i < 15; i++) fProcessFlag[i] = -1;
1df5fa54 65}
66
fb2cbbec 67void TFlukaConfigOption::SetCut(const char* flagname, Double_t val)
1df5fa54 68{
fb2cbbec 69 // Set a cut value
70 const TString cuts[11] =
71 {"CUTGAM", "CUTELE", "CUTNEU", "CUTHAD", "CUTMUO", "BCUTE", "BCUTM", "DCUTE", "DCUTM", "PPCUTM", "TOFMAX"};
72 Int_t i;
73 for (i = 0; i < 11; i++) {
74 if (cuts[i].CompareTo(flagname) == 0) {
75 fCutValue[i] = val;
76 if (fMedium == -1) fgDCutValue[i] = val;
77 break;
78 }
79 }
1df5fa54 80}
81
fb2cbbec 82void TFlukaConfigOption::SetProcess(const char* flagname, Int_t flag)
83{
84 // Set a process flag
85 const TString process[15] =
86 {"DCAY", "PAIR", "COMP", "PHOT", "PFIS", "DRAY", "ANNI", "BREM", "MUNU", "CKOV",
87 "HADR", "LOSS", "MULS", "RAYL", "STRA"};
37656ca5 88
fb2cbbec 89 Int_t i;
90 for (i = 0; i < 15; i++) {
91 if (process[i].CompareTo(flagname) == 0) {
92 fProcessFlag[i] = flag;
93 if (fMedium == -1) fgDProcessFlag[i] = flag;
94 break;
95 }
96 }
97}
1df5fa54 98
fb2cbbec 99void TFlukaConfigOption::WriteFlukaInputCards()
1df5fa54 100{
fb2cbbec 101 // Write the FLUKA input cards for the set of process flags and cuts
102 //
103 //
f79639f2 104 //
105 // Check if global option or medium specific
106
107 Bool_t mediumIsSensitive = kFALSE;
108 TGeoMedium* med = 0x0;
109 TGeoMedium* medium = 0x0;
110 TGeoMaterial* mat = 0x0;
111
81f1d030 112 if (fMedium != -1) {
89aa6d28 113 TFluka* fluka = (TFluka*) gMC;
81f1d030 114 fMedium = fgGeom->GetFlukaMaterial(fMedium);
115 //
116 // Check if material is actually used
117 Int_t* reglist;
118 Int_t nreg;
119 reglist = fgGeom->GetMaterialList(fMedium, nreg);
120 if (nreg == 0) {
cbd7b44e 121 // Material not used -- return
81f1d030 122 return;
123 }
f79639f2 124 //
125 // Find material
89aa6d28 126 TObjArray *matList = fluka->GetFlukaMaterials();
127 Int_t nmaterial = matList->GetEntriesFast();
8134b16b 128 fCMaterial = 0;
89aa6d28 129 for (Int_t im = 0; im < nmaterial; im++)
130 {
8134b16b 131 fCMaterial = dynamic_cast<TGeoMaterial*> (matList->At(im));
132 Int_t idmat = fCMaterial->GetIndex();
89aa6d28 133 if (idmat == fMedium) break;
f79639f2 134 } // materials
135 //
136 // Find medium
137 TList *medlist = gGeoManager->GetListOfMedia();
138 TIter next(medlist);
139 while((med = (TGeoMedium*)next()))
140 {
141 mat = med->GetMaterial();
142 if (mat->GetIndex() == fMedium) {
143 medium = med;
144 break;
145 }
146 } // media
147 //
148 // Check if sensitive
149 if (medium->GetParam(0) != 0.) mediumIsSensitive = kTRUE;
150
89aa6d28 151
8134b16b 152 fprintf(fgFile,"*\n*Material specific process and cut settings for #%8d %s\n", fMedium, fCMaterial->GetName());
fb2cbbec 153 fCMatMin = fMedium;
154 fCMatMax = fMedium;
155 } else {
156 fprintf(fgFile,"*\n*Global process and cut settings \n");
157 fCMatMin = fgMatMin;
158 fCMatMax = fgMatMax;
159 }
160
161//
162// Handle Process Flags
f79639f2 163//
164//
165// First make sure that all cuts are taken into account
166 if (DefaultProcessFlag(kPAIR) > 0 && fProcessFlag[kPAIR] == -1 && (fCutValue[kCUTELE] >= 0. || fCutValue[kPPCUTM] >= 0.))
167 fProcessFlag[kPAIR] = DefaultProcessFlag(kPAIR);
168 if (DefaultProcessFlag(kBREM) > 0 && fProcessFlag[kBREM] == -1 && (fCutValue[kBCUTE] >= 0. || fCutValue[kBCUTM] >= 0.))
169 fProcessFlag[kBREM] = DefaultProcessFlag(kBREM);
170 if (DefaultProcessFlag(kDRAY) > 0 && fProcessFlag[kDRAY] == -1 && (fCutValue[kDCUTE] >= 0. || fCutValue[kDCUTM] >= 0.))
171 fProcessFlag[kDRAY] = DefaultProcessFlag(kDRAY);
fb2cbbec 172//
f79639f2 173//
fb2cbbec 174 if (fProcessFlag[kDCAY] != -1) ProcessDCAY();
175 if (fProcessFlag[kPAIR] != -1) ProcessPAIR();
176 if (fProcessFlag[kBREM] != -1) ProcessBREM();
177 if (fProcessFlag[kCOMP] != -1) ProcessCOMP();
178 if (fProcessFlag[kPHOT] != -1) ProcessPHOT();
179 if (fProcessFlag[kPFIS] != -1) ProcessPFIS();
180 if (fProcessFlag[kANNI] != -1) ProcessANNI();
181 if (fProcessFlag[kMUNU] != -1) ProcessMUNU();
182 if (fProcessFlag[kHADR] != -1) ProcessHADR();
183 if (fProcessFlag[kMULS] != -1) ProcessMULS();
184 if (fProcessFlag[kRAYL] != -1) ProcessRAYL();
185
186 if (fProcessFlag[kLOSS] != -1 || fProcessFlag[kDRAY] != -1) ProcessLOSS();
187 if ((fMedium == -1 && fProcessFlag[kCKOV] > 0) || (fMedium > -1 && fProcessFlag[kCKOV] != -1)) ProcessCKOV();
188
189//
190// Handle Cuts
191//
192 if (fCutValue[kCUTGAM] >= 0.) ProcessCUTGAM();
193 if (fCutValue[kCUTELE] >= 0.) ProcessCUTELE();
194 if (fCutValue[kCUTNEU] >= 0.) ProcessCUTNEU();
195 if (fCutValue[kCUTHAD] >= 0.) ProcessCUTHAD();
196 if (fCutValue[kCUTMUO] >= 0.) ProcessCUTMUO();
f79639f2 197//
198// Time of flight
fb2cbbec 199 if (fCutValue[kTOFMAX] >= 0.) ProcessTOFMAX();
f79639f2 200
201//
202// Tracking precission
203 if (mediumIsSensitive) ProcessSensitiveMedium();
1df5fa54 204}
205
fb2cbbec 206void TFlukaConfigOption::ProcessDCAY()
1df5fa54 207{
fb2cbbec 208 // Process DCAY option
209 fprintf(fgFile,"*\n* --- DCAY --- Decays. Flag = %5d\n", fProcessFlag[kDCAY]);
210 if (fProcessFlag[kDCAY] == 0) {
211 printf("Decays cannot be switched off \n");
212 } else {
213 fprintf(fgFile, "* Decays are on by default\n");
214 }
215}
216
217
218void TFlukaConfigOption::ProcessPAIR()
219{
220 // Process PAIR option
81f1d030 221 fprintf(fgFile,"*\n* --- PAIR --- Pair production by gammas, muons and hadrons. Flag = %5d, PPCUTM = %13.4g, PPCUTE = %13.4g\n",
222 fProcessFlag[kPAIR], fCutValue[kCUTELE], fCutValue[kPPCUTM]);
fb2cbbec 223 //
224 // gamma -> e+ e-
225 //
226 if (fProcessFlag[kPAIR] > 0) {
81f1d030 227 fprintf(fgFile,"EMFCUT %10.1f%10.1f%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",0., 0., 0.0, fCMatMin, fCMatMax, 1.);
fb2cbbec 228 } else {
229 fprintf(fgFile,"EMFCUT %10.1f%10.1f%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",0., 0., 1e10, fCMatMin, fCMatMax, 1.);
230 }
231
232 //
233 // Direct pair production by Muons and Hadrons
234 //
235 //
236 // Attention ! This card interferes with BREM
237 //
238
239 if (fProcessFlag[kBREM] == -1 ) fProcessFlag[kBREM] = fgDProcessFlag[kBREM];
240 if (fCutValue[kBCUTM] == -1.) fCutValue[kBCUTM] = fgDCutValue[kBCUTM];
241
242
243 Float_t flag = -3.;
244 if (fProcessFlag[kPAIR] > 0 && fProcessFlag[kBREM] == 0) flag = 1.;
245 if (fProcessFlag[kPAIR] == 0 && fProcessFlag[kBREM] > 0) flag = 2.;
246 if (fProcessFlag[kPAIR] > 0 && fProcessFlag[kBREM] > 0) flag = 3.;
247 if (fProcessFlag[kPAIR] == 0 && fProcessFlag[kBREM] == 0) flag = -3.;
248 // Flag BREM card as handled
37656ca5 249 fProcessFlag[kBREM] = - fProcessFlag[kBREM];
fb2cbbec 250
251 //
252 // Energy cut for pair prodution
253 //
254 Float_t cutP = fCutValue[kPPCUTM];
255 if (fCutValue[kPPCUTM] == -1.) cutP = fgDCutValue[kPPCUTM];
256 // In G3 this is the cut on the total energy of the e+e- pair
257 // In FLUKA the cut is on the kinetic energy of the electron and poistron
258 cutP = cutP / 2. - 0.51099906e-3;
259 if (cutP < 0.) cutP = 0.;
260 // No explicite generation of e+/e-
261 if (fProcessFlag[kPAIR] == 2) cutP = -1.;
262 //
263 // Energy cut for bremsstrahlung
264 //
265 Float_t cutB = 0.;
266 if (flag > 1.) {
267 fprintf(fgFile,"*\n* +++ BREM --- Bremsstrahlung by muons/hadrons. Flag = %5d, BCUTM = %13.4g \n",
268 fProcessFlag[kBREM], fCutValue[kBCUTM]);
269
270 cutB = fCutValue[kBCUTM];
271 // No explicite production of gammas
272 if (fProcessFlag[kBREM] == 2) cutB = -1.;
273 }
274
275 fprintf(fgFile,"PAIRBREM %10.1f%10.4g%10.4g%10.1f%10.1f\n",flag, cutP, cutB, fCMatMin, fCMatMax);
276}
277
278
279void TFlukaConfigOption::ProcessBREM()
280{
281 // Process BREM option
282 fprintf(fgFile,"*\n* --- BREM --- Bremsstrahlung by e+/- and muons/hadrons. Flag = %5d, BCUTE = %13.4g, BCUTM = %13.4g \n",
283 fProcessFlag[kBREM], fCutValue[kBCUTE], fCutValue[kBCUTM]);
284
285 //
286 // e+/- -> e+/- gamma
287 //
288 Float_t cutB = fCutValue[kBCUTE];
289 if (fCutValue[kBCUTE] == -1.) cutB = fgDCutValue[kBCUTE];
290
291
37656ca5 292 if (TMath::Abs(fProcessFlag[kBREM]) > 0) {
fb2cbbec 293 fprintf(fgFile,"EMFCUT %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fELPO-THR\n",cutB, 0., 0., fCMatMin, fCMatMax, 1.);
294 } else {
295 fprintf(fgFile,"EMFCUT %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fELPO-THR\n",1.e10, 0., 0., fCMatMin, fCMatMax, 1.);
296 }
297
298 //
299 // Bremsstrahlung by muons and hadrons
300 //
301 cutB = fCutValue[kBCUTM];
302 if (fCutValue[kBCUTM] == -1.) cutB = fgDCutValue[kBCUTM];
303 if (fProcessFlag[kBREM] == 2) cutB = -1.;
304 Float_t flag = 2.;
305 if (fProcessFlag[kBREM] == 0) flag = -2.;
306
307 fprintf(fgFile,"PAIRBREM %10.1f%10.4g%10.4g%10.1f%10.1f\n", flag, 0., cutB, fCMatMin, fCMatMax);
308}
309
310void TFlukaConfigOption::ProcessCOMP()
311{
312 // Process COMP option
313 fprintf(fgFile,"*\n* --- COMP --- Compton scattering Flag = %5d \n", fProcessFlag[kCOMP]);
314
315 //
316 // Compton scattering
317 //
318
319 if (fProcessFlag[kCOMP] > 0) {
320 fprintf(fgFile,"EMFCUT %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",0. , 0., 0., fCMatMin, fCMatMax, 1.);
321 } else {
322 fprintf(fgFile,"EMFCUT %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",1.e10, 0., 0., fCMatMin, fCMatMax, 1.);
323 }
324}
325
326void TFlukaConfigOption::ProcessPHOT()
327{
328 // Process PHOS option
329 fprintf(fgFile,"*\n* --- PHOT --- Photoelectric effect. Flag = %5d\n", fProcessFlag[kPHOT]);
330
331 //
332 // Photoelectric effect
333 //
334
335 if (fProcessFlag[kPHOT] > 0) {
81f1d030 336 fprintf(fgFile,"EMFCUT %10.4g%10.4g%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",0., 0., 0., fCMatMin, fCMatMax, 1.);
fb2cbbec 337 } else {
338 fprintf(fgFile,"EMFCUT %10.1f%10.4g%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",0., 1.e10, 0., fCMatMin, fCMatMax, 1.);
339 }
340}
341
342void TFlukaConfigOption::ProcessANNI()
343{
344 // Process ANNI option
345 fprintf(fgFile,"*\n* --- ANNI --- Positron annihilation. Flag = %5d \n", fProcessFlag[kANNI]);
346
347 //
348 // Positron annihilation
349 //
350
351 if (fProcessFlag[kANNI] > 0) {
352 fprintf(fgFile,"EMFCUT %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fANNH-THR\n",0. , 0., 0., fCMatMin, fCMatMax, 1.);
353 } else {
354 fprintf(fgFile,"EMFCUT %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fANNH-THR\n",1.e10, 0., 0., fCMatMin, fCMatMax, 1.);
355 }
356}
357
358
359void TFlukaConfigOption::ProcessPFIS()
360{
361 // Process PFIS option
362 fprintf(fgFile,"*\n* --- PFIS --- Photonuclear interaction Flag = %5d\n", fProcessFlag[kPFIS]);
363
364 //
365 // Photonuclear interactions
366 //
367
368 if (fProcessFlag[kPFIS] > 0) {
369 fprintf(fgFile,"PHOTONUC %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",(Float_t) fProcessFlag[kPFIS], 0., 0., fCMatMin, fCMatMax, 1.);
370 } else {
371 fprintf(fgFile,"PHOTONUC %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",-1. , 0., 0., fCMatMin, fCMatMax, 1.);
372 }
373}
374
375void TFlukaConfigOption::ProcessMUNU()
376{
377 // Process MUNU option
378 fprintf(fgFile,"*\n* --- MUNU --- Muon nuclear interaction. Flag = %5d\n", fProcessFlag[kMUNU]);
379
380 //
381 // Muon nuclear interactions
382 //
383 if (fProcessFlag[kMUNU] > 0) {
384 fprintf(fgFile,"MUPHOTON %10.1f%10.3f%10.3f%10.1f%10.1f%10.1f\n",(Float_t )fProcessFlag[kMUNU], 0.25, 0.75, fCMatMin, fCMatMax, 1.);
385 } else {
386 fprintf(fgFile,"MUPHOTON %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",-1. , 0., 0., fCMatMin, fCMatMax, 1.);
387 }
388}
389
390void TFlukaConfigOption::ProcessRAYL()
391{
392 // Process RAYL option
393 fprintf(fgFile,"*\n* --- RAYL --- Rayleigh Scattering. Flag = %5d\n", fProcessFlag[kRAYL]);
394
395 //
396 // Rayleigh scattering
397 //
398 Int_t nreg;
399 Int_t* reglist = fgGeom->GetMaterialList(fMedium, nreg);
400 //
401 // Loop over regions of a given material
402 for (Int_t k = 0; k < nreg; k++) {
403 Float_t ireg = reglist[k];
404 if (fProcessFlag[kRAYL] > 0) {
405 fprintf(fgFile,"EMFRAY %10.1f%10.1f%10.1f%10.1f\n", 1., ireg, ireg, 1.);
406 } else {
407 fprintf(fgFile,"EMFRAY %10.1f%10.1f%10.1f%10.1f\n", 3., ireg, ireg, 1.);
408 }
409 }
410}
411
412void TFlukaConfigOption::ProcessCKOV()
413{
414 // Process CKOV option
415 fprintf(fgFile,"*\n* --- CKOV --- Cerenkov Photon production. %5d\n", fProcessFlag[kCKOV]);
416
417 //
418 // Cerenkov photon production
419 //
420
421 TFluka* fluka = (TFluka*) gMC;
422 TObjArray *matList = fluka->GetFlukaMaterials();
423 Int_t nmaterial = matList->GetEntriesFast();
424 for (Int_t im = 0; im < nmaterial; im++)
425 {
426 TGeoMaterial* material = dynamic_cast<TGeoMaterial*> (matList->At(im));
427 Int_t idmat = material->GetIndex();
428//
429// Check if global option
430 if (fMedium != -1 && idmat != fMedium) continue;
431
432 TFlukaCerenkov* cerenkovProp;
433 if (!(cerenkovProp = dynamic_cast<TFlukaCerenkov*>(material->GetCerenkovProperties()))) continue;
434 //
435 // This medium has Cerenkov properties
436 //
437 //
438 if (fMedium == -1 || (fMedium != -1 && fProcessFlag[kCKOV] > 0)) {
439 // Write OPT-PROD card for each medium
440 Float_t emin = cerenkovProp->GetMinimumEnergy();
441 Float_t emax = cerenkovProp->GetMaximumEnergy();
442 fprintf(fgFile, "OPT-PROD %10.4g%10.4g%10.4g%10.4g%10.4g%10.4gCERENKOV\n", emin, emax, 0.,
443 Float_t(idmat), Float_t(idmat), 0.);
444 //
445 // Write OPT-PROP card for each medium
446 // Forcing FLUKA to call user routines (queffc.cxx, rflctv.cxx, rfrndx.cxx)
447 //
448 fprintf(fgFile, "OPT-PROP %10.4g%10.4g%10.4g%10.1f%10.1f%10.1fWV-LIMIT\n",
449 cerenkovProp->GetMinimumWavelength(), cerenkovProp->GetMaximumWavelength(), cerenkovProp->GetMaximumWavelength(),
450 Float_t(idmat), Float_t(idmat), 0.0);
451
2906553d 452
453 fprintf(fgFile, "OPT-PROP %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n", -100., -100., -100.,
454 Float_t(idmat), Float_t(idmat), 0.0);
fb2cbbec 455
456 for (Int_t j = 0; j < 3; j++) {
457 fprintf(fgFile, "OPT-PROP %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f&\n", -100., -100., -100.,
458 Float_t(idmat), Float_t(idmat), 0.0);
459 }
2906553d 460
461
462 // Photon detection efficiency user defined
fb2cbbec 463 if (cerenkovProp->IsSensitive())
464 fprintf(fgFile, "OPT-PROP %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fSENSITIV\n", -100., -100., -100.,
465 Float_t(idmat), Float_t(idmat), 0.0);
2906553d 466 // Material has a reflective surface
467 if (cerenkovProp->IsMetal())
468 fprintf(fgFile, "OPT-PROP %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fMETAL\n", -100., -100., -100.,
469 Float_t(idmat), Float_t(idmat), 0.0);
470
fb2cbbec 471 } else {
472 fprintf(fgFile,"OPT-PROD %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fCERE-OFF\n",0., 0., 0., fCMatMin, fCMatMax, 1.);
473 }
474 }
475}
476
477
478void TFlukaConfigOption::ProcessHADR()
479{
480 // Process HADR option
481 fprintf(fgFile,"*\n* --- HADR --- Hadronic interactions. Flag = %5d\n", fProcessFlag[kHADR]);
482
483 if (fProcessFlag[kHADR] > 0) {
484 fprintf(fgFile,"*\n*Hadronic interaction is ON by default in FLUKA\n");
485 } else {
486 if (fMedium != -1) printf("Hadronic interactions cannot be switched off material by material !\n");
487 fprintf(fgFile,"THRESHOL %10.1f%10.1f%10.1f%10.1e%10.1f\n",0., 0., 0., 1.e10, 0.);
488 }
489}
490
491
492
493void TFlukaConfigOption::ProcessMULS()
494{
495 // Process MULS option
496 fprintf(fgFile,"*\n* --- MULS --- Muliple Scattering. Flag = %5d\n", fProcessFlag[kMULS]);
497 //
498 // Multiple scattering
499 //
500 if (fProcessFlag[kMULS] > 0) {
501 fprintf(fgFile,"*\n*Multiple scattering is ON by default in FLUKA\n");
502 } else {
503 fprintf(fgFile,"MULSOPT %10.1f%10.1f%10.1f%10.1f%10.1f\n",0., 3., 3., fCMatMin, fCMatMax);
504 }
505}
506
507void TFlukaConfigOption::ProcessLOSS()
508{
509 // Process LOSS option
510 fprintf(fgFile,"*\n* --- LOSS --- Ionisation energy loss. Flags: LOSS = %5d, DRAY = %5d, STRA = %5d; Cuts: DCUTE = %13.4g, DCUTM = %13.4g \n",
511 fProcessFlag[kLOSS], fProcessFlag[kDRAY], fProcessFlag[kSTRA], fCutValue[kDCUTE], fCutValue[kDCUTM]);
512 //
513 // Ionisation energy loss
514 //
515 //
516 // Impose consistency
517
518 if (fProcessFlag[kLOSS] == 1 || fProcessFlag[kLOSS] == 3) {
519 fProcessFlag[kDRAY] = 1;
520 } else if (fProcessFlag[kLOSS] == 2) {
521 fProcessFlag[kDRAY] = 0;
522 fCutValue[kDCUTE] = 1.e10;
523 fCutValue[kDCUTM] = 1.e10;
524 } else {
525 if (fProcessFlag[kDRAY] == 1) {
526 fProcessFlag[kLOSS] = 1;
527 } else if (fProcessFlag[kDRAY] == 0) {
528 fProcessFlag[kLOSS] = 2;
529 fCutValue[kDCUTE] = 1.e10;
530 fCutValue[kDCUTM] = 1.e10;
531 }
532 }
533
534 if (fCutValue[kDCUTE] == -1.) fCutValue[kDCUTE] = fgDCutValue[kDCUTE];
535 if (fCutValue[kDCUTM] == -1.) fCutValue[kDCUTM] = fgDCutValue[kDCUTM];
536
537 Float_t cutM = fCutValue[kDCUTM];
538
539
540 if (fProcessFlag[kSTRA] == -1) fProcessFlag[kSTRA] = fgDProcessFlag[kSTRA];
541 if (fProcessFlag[kSTRA] == 0) fProcessFlag[kSTRA] = 1;
542 Float_t stra = (Float_t) fProcessFlag[kSTRA];
543
544
545 if (fProcessFlag[kLOSS] == 1 || fProcessFlag[kLOSS] == 3) {
546//
547// Restricted energy loss fluctuations
548//
549 fprintf(fgFile,"IONFLUCT %10.1f%10.1f%10.1f%10.1f%10.1f\n", 1., 1., stra, fCMatMin, fCMatMax);
2906553d 550 fprintf(fgFile,"DELTARAY %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n", cutM, 0., 0., fCMatMin, fCMatMax, 1.);
fb2cbbec 551 } else if (fProcessFlag[kLOSS] == 4) {
552//
553// No fluctuations
554//
555 fprintf(fgFile,"IONFLUCT %10.1f%10.1f%10.1f%10.1f%10.1f\n",-1., -1., stra, fCMatMin, fCMatMax);
556 fprintf(fgFile,"DELTARAY %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n", 1.e10, 0., 0., fCMatMin, fCMatMax, 1.);
557 } else {
558//
559// Full fluctuations
560//
561 fprintf(fgFile,"IONFLUCT %10.1f%10.1f%10.1f%10.1f%10.1f\n",1., 1., stra, fCMatMin, fCMatMax);
562 fprintf(fgFile,"DELTARAY %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n", 1.e10, 0., 0., fCMatMin, fCMatMax, 1.);
563 }
564}
565
566
567void TFlukaConfigOption::ProcessCUTGAM()
568{
569// Cut on gammas
570//
571 fprintf(fgFile,"*\n*Cut for Gammas. CUTGAM = %13.4g\n", fCutValue[kCUTGAM]);
572 if (fMedium == -1) {
573 fprintf(fgFile,"EMFCUT %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n",
574 0., fCutValue[kCUTGAM], 0., 0., Float_t(fgGeom->NofVolumes()), 1.);
46bee1bf 575
fb2cbbec 576 } else {
577 Int_t nreg, *reglist;
578 Float_t ireg;
579 reglist = fgGeom->GetMaterialList(fMedium, nreg);
580 // Loop over regions of a given material
581 for (Int_t k = 0; k < nreg; k++) {
582 ireg = reglist[k];
583 fprintf(fgFile,"EMFCUT %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", 0.,fCutValue[kCUTGAM], 0., ireg, ireg, 1.);
584 }
585 }
8134b16b 586
587 // Transport production cut used for pemf
588 //
589 // FUDGEM paramter. The parameter takes into account th contribution of atomic electrons to multiple scattering.
590 // For production and transport cut-offs larger than 100 keV it must be set = 1.0, while in the keV region it must be
591 Float_t parFudgem = (fCutValue[kCUTGAM] > 1.e-4)? 1.0 : 0.0 ;
46bee1bf 592 fprintf(fgFile,"EMFCUT %10.4g%10.4g%10.1f%10.1f%10.1f%10.1fPROD-CUT\n",
8134b16b 593 0., fCutValue[kCUTGAM], parFudgem, fCMatMin, fCMatMax, 1.);
fb2cbbec 594}
595
596void TFlukaConfigOption::ProcessCUTELE()
597{
598// Cut on e+/e-
599//
600 fprintf(fgFile,"*\n*Cut for e+/e-. CUTELE = %13.4g\n", fCutValue[kCUTELE]);
601 if (fMedium == -1) {
602 fprintf(fgFile,"EMFCUT %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n",
603 -fCutValue[kCUTELE], 0., 0., 0., Float_t(fgGeom->NofVolumes()), 1.);
604 } else {
605 Int_t nreg, *reglist;
606 Float_t ireg;
607 reglist = fgGeom->GetMaterialList(fMedium, nreg);
608 // Loop over regions of a given material
609 for (Int_t k = 0; k < nreg; k++) {
610 ireg = reglist[k];
611 fprintf(fgFile,"EMFCUT %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", -fCutValue[kCUTELE], 0., 0., ireg, ireg, 1.);
612 }
613 }
8134b16b 614 // Transport production cut used for pemf
615 //
616 // FUDGEM paramter. The parameter takes into account th contribution of atomic electrons to multiple scattering.
617 // For production and transport cut-offs larger than 100 keV it must be set = 1.0, while in the keV region it must be
618 Float_t parFudgem = (fCutValue[kCUTELE] > 1.e-4)? 1.0 : 0.0;
46bee1bf 619 fprintf(fgFile,"EMFCUT %10.4g%10.4g%10.1f%10.1f%10.1f%10.1fPROD-CUT\n",
8134b16b 620 -fCutValue[kCUTELE], 0., parFudgem, fCMatMin, fCMatMax, 1.);
fb2cbbec 621}
622
623void TFlukaConfigOption::ProcessCUTNEU()
624{
501e3a8d 625 // Cut on neutral hadrons
626 fprintf(fgFile,"*\n*Cut for neutral hadrons. CUTNEU = %13.4g\n", fCutValue[kCUTNEU]);
627
628 // Cut on neutral hadrons
629 fprintf(fgFile,"*\n*Cut for neutral hadrons. CUTNEU = %13.4g\n", fCutValue[kCUTNEU]);
630
631 // Energy group structure of the 72-neutron ENEA data set:
632 const Float_t neutronGroupUpLimit[72] = {
633 1.9600E-02, 1.7500E-02, 1.4918E-02, 1.3499E-02, 1.2214E-02, 1.1052E-02, 1.0000E-02, 9.0484E-03,
634 8.1873E-03, 7.4082E-03, 6.7032E-03, 6.0653E-03, 5.4881E-03, 4.9659E-03, 4.4933E-03, 4.0657E-03,
635 3.6788E-03, 3.3287E-03, 3.0119E-03, 2.7253E-03, 2.4660E-03, 2.2313E-03, 2.0190E-03, 1.8268E-03,
636 1.6530E-03, 1.4957E-03, 1.3534E-03, 1.2246E-03, 1.1080E-03, 1.0026E-03, 9.0718E-04, 8.2085E-04,
637 7.4274E-04, 6.0810E-04, 4.9787E-04, 4.0762E-04, 3.3373E-04, 2.7324E-04, 2.2371E-04, 1.8316E-04,
638 1.4996E-04, 1.2277E-04, 8.6517E-05, 5.2475E-05, 3.1828E-05, 2.1852E-05, 1.5034E-05, 1.0332E-05,
639 7.1018E-06, 4.8809E-06, 3.3546E-06, 2.3054E-06, 1.5846E-06, 1.0446E-06, 6.8871E-07, 4.5400E-07,
640 2.7537E-07, 1.6702E-07, 1.0130E-07, 6.1442E-08, 3.7267E-08, 2.2603E-08, 1.5535E-08, 1.0677E-08,
641 7.3375E-09, 5.0435E-09, 3.4662E-09, 2.3824E-09, 1.6374E-09, 1.1254E-09, 6.8257E-10, 4.1400E-10
642 };
643
644 Float_t cut = fCutValue[kCUTNEU];
645 //
646 // 8.0 = Neutron
647 // 9.0 = Antineutron
648 // Find the FLUKA neutron group corresponding to the cut
649 //
650 Float_t neutronCut = cut;
651 Int_t groupCut = 1; // if cut is > 19.6 MeV not low energy neutron transport is done
652 if (neutronCut < 0.0196) {
653 neutronCut = 0.0196;
654 // Search the group cutoff for the energy cut
655 Int_t i;
656 for( i=0; i<72; i++ ) {
657 if (cut > neutronGroupUpLimit[i]) break;
658 }
659 groupCut = i+1;
660 }
661
662 if (fMedium == -1) {
663 Float_t cut = fCutValue[kCUTNEU];
664 // 8.0 = Neutron
665 // 9.0 = Antineutron
666 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -neutronCut, 8.0, 9.0);
667 fprintf(fgFile,"LOW-BIAS %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n",
668 Float_t(groupCut), 73.0, 0.95, 2., Float_t(fgGeom->NofVolumes()), 1.);
669 //
2906553d 670 //
fb2cbbec 671 // 12.0 = Kaon zero long
672 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 12.0, 12.0);
673 // 17.0 = Lambda, 18.0 = Antilambda
674 // 19.0 = Kaon zero short
675 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 17.0, 19.0);
676 // 22.0 = Sigma zero, Pion zero, Kaon zero
677 // 25.0 = Antikaon zero
678 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 22.0, 25.0);
679 // 32.0 = Antisigma zero
680 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 32.0, 32.0);
681 // 34.0 = Xi zero
682 // 35.0 = AntiXi zero
683 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 34.0, 35.0);
684 // 47.0 = D zero
685 // 48.0 = AntiD zero
686 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 47.0, 48.0);
687 // 53.0 = Xi_c zero
688 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 53.0, 53.0);
689 // 55.0 = Xi'_c zero
690 // 56.0 = Omega_c zero
691 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 55.0, 56.0);
692 // 59.0 = AntiXi_c zero
693 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 59.0, 59.0);
694 // 61.0 = AntiXi'_c zero
695 // 62.0 = AntiOmega_c zero
696 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 61.0, 62.0);
697 } else {
501e3a8d 698 Int_t nreg, *reglist;
699 Float_t ireg;
700 reglist = fgGeom->GetMaterialList(fMedium, nreg);
701 // Loop over regions of a given material
702 for (Int_t k = 0; k < nreg; k++) {
703 ireg = reglist[k];
704 fprintf(fgFile,"LOW-BIAS %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n",
705 Float_t(groupCut), 73.0, 0.95, ireg, ireg, 1.);
706 }
707
8134b16b 708 Warning("ProcessCUTNEU",
709 "Material #%4d %s: Cut on neutral hadrons (Ekin > %9.3e) material by material only implemented for low-energy neutrons !\n",
710 fMedium, fCMaterial->GetName(), cut);
fb2cbbec 711 }
712}
713
714void TFlukaConfigOption::ProcessCUTHAD()
715{
716 // Cut on charged hadrons
717 fprintf(fgFile,"*\n*Cut for charge hadrons. CUTHAD = %13.4g\n", fCutValue[kCUTHAD]);
8134b16b 718 Float_t cut = fCutValue[kCUTHAD];
fb2cbbec 719 if (fMedium == -1) {
fb2cbbec 720 // 1.0 = Proton
721 // 2.0 = Antiproton
722 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 1.0, 2.0);
723 // 13.0 = Positive Pion, Negative Pion, Positive Kaon
724 // 16.0 = Negative Kaon
725 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 13.0, 16.0);
726 // 20.0 = Negative Sigma
727 // 21.0 = Positive Sigma
728 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 20.0, 21.0);
729 // 31.0 = Antisigma minus
730 // 33.0 = Antisigma plus
731 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 31.0, 31.0);
732 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 33.0, 33.0);
733 // 36.0 = Negative Xi, Positive Xi, Omega minus
734 // 39.0 = Antiomega
735 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 36.0, 39.0);
736 // 45.0 = D plus
737 // 46.0 = D minus
738 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 45.0, 46.0);
739 // 49.0 = D_s plus, D_s minus, Lambda_c plus
740 // 52.0 = Xi_c plus
741 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 49.0, 52.0);
742 // 54.0 = Xi'_c plus
743 // 60.0 = AntiXi'_c minus
744 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 54.0, 54.0);
745 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 60.0, 60.0);
746 // 57.0 = Antilambda_c minus
747 // 58.0 = AntiXi_c minus
748 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 57.0, 58.0);
749 } else {
8134b16b 750 Warning("ProcessCUTHAD",
751 "Material #%4d %s: Cut on charged hadrons (Ekin > 9.3e) material by material not yet implemented !\n",
752 fMedium, fCMaterial->GetName(), cut);
fb2cbbec 753 }
754}
755
756void TFlukaConfigOption::ProcessCUTMUO()
757{
758 // Cut on muons
759 fprintf(fgFile,"*\n*Cut for muons. CUTMUO = %13.4g\n", fCutValue[kCUTMUO]);
760 Float_t cut = fCutValue[kCUTMUO];
761 if (fMedium == -1) {
762 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n",-cut, 10.0, 11.0);
763 } else {
8134b16b 764 Warning("ProcessCUTMUO", "Material #%4d %s: Cut on muons (Ekin > %9.3e) material by material not yet implemented !\n",
765 fMedium, fCMaterial->GetName(), cut);
fb2cbbec 766 }
767
768
769}
770
771void TFlukaConfigOption::ProcessTOFMAX()
772{
773 // Cut time of flight
774 Float_t cut = fCutValue[kTOFMAX];
775 fprintf(fgFile,"*\n*Cut on time of flight. TOFMAX = %13.4g\n", fCutValue[kTOFMAX]);
776 fprintf(fgFile,"TIME-CUT %10.4g%10.1f%10.1f%10.1f%10.1f\n",cut*1.e9,0.,0.,-6.0,64.0);
1df5fa54 777}
f79639f2 778
779void TFlukaConfigOption::ProcessSensitiveMedium()
780{
781 //
782 // Special options for sensitive media
783 //
784
785 fprintf(fgFile,"*\n*Options for sensitive medium \n");
786 //
787 // EMFFIX
788 fprintf(fgFile,"EMFFIX %10.1f%10.3f%10.1f%10.1f%10.1f%10.1f\n", fCMatMin, 0.05, 0., 0., 0., 0.);
789 //
790 // FLUKAFIX
791 fprintf(fgFile,"FLUKAFIX %10.3f %10.3f\n", 0.05, fCMatMin);
792}