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