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