]> 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
Issue AliError instead of AliFatal in case of inconsistent number of tracks
[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>
24#include <TObjArray.h>
25#include <TVirtualMC.h>
26#include <TGeoMaterial.h>
27
28Float_t TFlukaConfigOption::fgMatMin(-1.);
29Float_t TFlukaConfigOption::fgMatMax(-1.);
30FILE* TFlukaConfigOption::fgFile(0x0);
31TFlukaMCGeometry* TFlukaConfigOption::fgGeom(0x0);
32
33Double_t TFlukaConfigOption::fgDCutValue[11];
34Int_t TFlukaConfigOption::fgDProcessFlag[15];
35
36
1df5fa54 37ClassImp(TFlukaConfigOption)
38
39
40TFlukaConfigOption::TFlukaConfigOption()
41{
42 // Default constructor
fb2cbbec 43 fMedium = -1;
44 fCMatMin = -1.;
45 fCMatMax = -1.;
46 Int_t i;
47 for (i = 0; i < 11; i++) fCutValue[i] = -1.;
48 for (i = 0; i < 15; i++) fProcessFlag[i] = -1;
1df5fa54 49}
50
51
fb2cbbec 52TFlukaConfigOption::TFlukaConfigOption(Int_t medium)
1df5fa54 53{
54 // Constructor
fb2cbbec 55 fMedium = medium;
56 fCMatMin = -1.;
57 fCMatMax = -1.;
58 Int_t i;
59 for (i = 0; i < 11; i++) fCutValue[i] = -1.;
60 for (i = 0; i < 15; i++) fProcessFlag[i] = -1;
1df5fa54 61}
62
fb2cbbec 63void TFlukaConfigOption::SetCut(const char* flagname, Double_t val)
1df5fa54 64{
fb2cbbec 65 // Set a cut value
66 const TString cuts[11] =
67 {"CUTGAM", "CUTELE", "CUTNEU", "CUTHAD", "CUTMUO", "BCUTE", "BCUTM", "DCUTE", "DCUTM", "PPCUTM", "TOFMAX"};
68 Int_t i;
69 for (i = 0; i < 11; i++) {
70 if (cuts[i].CompareTo(flagname) == 0) {
71 fCutValue[i] = val;
72 if (fMedium == -1) fgDCutValue[i] = val;
73 break;
74 }
75 }
1df5fa54 76}
77
fb2cbbec 78void TFlukaConfigOption::SetProcess(const char* flagname, Int_t flag)
79{
80 // Set a process flag
81 const TString process[15] =
82 {"DCAY", "PAIR", "COMP", "PHOT", "PFIS", "DRAY", "ANNI", "BREM", "MUNU", "CKOV",
83 "HADR", "LOSS", "MULS", "RAYL", "STRA"};
37656ca5 84
fb2cbbec 85 Int_t i;
86 for (i = 0; i < 15; i++) {
87 if (process[i].CompareTo(flagname) == 0) {
88 fProcessFlag[i] = flag;
89 if (fMedium == -1) fgDProcessFlag[i] = flag;
90 break;
91 }
92 }
93}
1df5fa54 94
fb2cbbec 95void TFlukaConfigOption::WriteFlukaInputCards()
1df5fa54 96{
fb2cbbec 97 // Write the FLUKA input cards for the set of process flags and cuts
98 //
99 //
81f1d030 100 if (fMedium != -1) {
89aa6d28 101 TFluka* fluka = (TFluka*) gMC;
81f1d030 102 fMedium = fgGeom->GetFlukaMaterial(fMedium);
103 //
104 // Check if material is actually used
105 Int_t* reglist;
106 Int_t nreg;
107 reglist = fgGeom->GetMaterialList(fMedium, nreg);
108 if (nreg == 0) {
109 printf("Material not used !\n");
110 return;
111 }
112
89aa6d28 113 TObjArray *matList = fluka->GetFlukaMaterials();
114 Int_t nmaterial = matList->GetEntriesFast();
115 TGeoMaterial* material = 0;
116 for (Int_t im = 0; im < nmaterial; im++)
117 {
118 material = dynamic_cast<TGeoMaterial*> (matList->At(im));
119 Int_t idmat = material->GetIndex();
120 if (idmat == fMedium) break;
121 }
122
123
124//
125// Check if global option
126
127 fprintf(fgFile,"*\n*Material specific process and cut settings for #%8d %s\n", fMedium, material->GetName());
fb2cbbec 128 fCMatMin = fMedium;
129 fCMatMax = fMedium;
130 } else {
131 fprintf(fgFile,"*\n*Global process and cut settings \n");
132 fCMatMin = fgMatMin;
133 fCMatMax = fgMatMax;
134 }
135
136//
137// Handle Process Flags
138//
139 if (fProcessFlag[kDCAY] != -1) ProcessDCAY();
140 if (fProcessFlag[kPAIR] != -1) ProcessPAIR();
141 if (fProcessFlag[kBREM] != -1) ProcessBREM();
142 if (fProcessFlag[kCOMP] != -1) ProcessCOMP();
143 if (fProcessFlag[kPHOT] != -1) ProcessPHOT();
144 if (fProcessFlag[kPFIS] != -1) ProcessPFIS();
145 if (fProcessFlag[kANNI] != -1) ProcessANNI();
146 if (fProcessFlag[kMUNU] != -1) ProcessMUNU();
147 if (fProcessFlag[kHADR] != -1) ProcessHADR();
148 if (fProcessFlag[kMULS] != -1) ProcessMULS();
149 if (fProcessFlag[kRAYL] != -1) ProcessRAYL();
150
151 if (fProcessFlag[kLOSS] != -1 || fProcessFlag[kDRAY] != -1) ProcessLOSS();
152 if ((fMedium == -1 && fProcessFlag[kCKOV] > 0) || (fMedium > -1 && fProcessFlag[kCKOV] != -1)) ProcessCKOV();
153
154//
155// Handle Cuts
156//
157 if (fCutValue[kCUTGAM] >= 0.) ProcessCUTGAM();
158 if (fCutValue[kCUTELE] >= 0.) ProcessCUTELE();
159 if (fCutValue[kCUTNEU] >= 0.) ProcessCUTNEU();
160 if (fCutValue[kCUTHAD] >= 0.) ProcessCUTHAD();
161 if (fCutValue[kCUTMUO] >= 0.) ProcessCUTMUO();
162
163 if (fCutValue[kTOFMAX] >= 0.) ProcessTOFMAX();
1df5fa54 164}
165
fb2cbbec 166void TFlukaConfigOption::ProcessDCAY()
1df5fa54 167{
fb2cbbec 168 // Process DCAY option
169 fprintf(fgFile,"*\n* --- DCAY --- Decays. Flag = %5d\n", fProcessFlag[kDCAY]);
170 if (fProcessFlag[kDCAY] == 0) {
171 printf("Decays cannot be switched off \n");
172 } else {
173 fprintf(fgFile, "* Decays are on by default\n");
174 }
175}
176
177
178void TFlukaConfigOption::ProcessPAIR()
179{
180 // Process PAIR option
81f1d030 181 fprintf(fgFile,"*\n* --- PAIR --- Pair production by gammas, muons and hadrons. Flag = %5d, PPCUTM = %13.4g, PPCUTE = %13.4g\n",
182 fProcessFlag[kPAIR], fCutValue[kCUTELE], fCutValue[kPPCUTM]);
fb2cbbec 183 //
184 // gamma -> e+ e-
185 //
186 if (fProcessFlag[kPAIR] > 0) {
81f1d030 187 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 188 } else {
189 fprintf(fgFile,"EMFCUT %10.1f%10.1f%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",0., 0., 1e10, fCMatMin, fCMatMax, 1.);
190 }
191
192 //
193 // Direct pair production by Muons and Hadrons
194 //
195 //
196 // Attention ! This card interferes with BREM
197 //
198
199 if (fProcessFlag[kBREM] == -1 ) fProcessFlag[kBREM] = fgDProcessFlag[kBREM];
200 if (fCutValue[kBCUTM] == -1.) fCutValue[kBCUTM] = fgDCutValue[kBCUTM];
201
202
203 Float_t flag = -3.;
204 if (fProcessFlag[kPAIR] > 0 && fProcessFlag[kBREM] == 0) flag = 1.;
205 if (fProcessFlag[kPAIR] == 0 && fProcessFlag[kBREM] > 0) flag = 2.;
206 if (fProcessFlag[kPAIR] > 0 && fProcessFlag[kBREM] > 0) flag = 3.;
207 if (fProcessFlag[kPAIR] == 0 && fProcessFlag[kBREM] == 0) flag = -3.;
208 // Flag BREM card as handled
37656ca5 209 fProcessFlag[kBREM] = - fProcessFlag[kBREM];
fb2cbbec 210
211 //
212 // Energy cut for pair prodution
213 //
214 Float_t cutP = fCutValue[kPPCUTM];
215 if (fCutValue[kPPCUTM] == -1.) cutP = fgDCutValue[kPPCUTM];
216 // In G3 this is the cut on the total energy of the e+e- pair
217 // In FLUKA the cut is on the kinetic energy of the electron and poistron
218 cutP = cutP / 2. - 0.51099906e-3;
219 if (cutP < 0.) cutP = 0.;
220 // No explicite generation of e+/e-
221 if (fProcessFlag[kPAIR] == 2) cutP = -1.;
222 //
223 // Energy cut for bremsstrahlung
224 //
225 Float_t cutB = 0.;
226 if (flag > 1.) {
227 fprintf(fgFile,"*\n* +++ BREM --- Bremsstrahlung by muons/hadrons. Flag = %5d, BCUTM = %13.4g \n",
228 fProcessFlag[kBREM], fCutValue[kBCUTM]);
229
230 cutB = fCutValue[kBCUTM];
231 // No explicite production of gammas
232 if (fProcessFlag[kBREM] == 2) cutB = -1.;
233 }
234
235 fprintf(fgFile,"PAIRBREM %10.1f%10.4g%10.4g%10.1f%10.1f\n",flag, cutP, cutB, fCMatMin, fCMatMax);
236}
237
238
239void TFlukaConfigOption::ProcessBREM()
240{
241 // Process BREM option
242 fprintf(fgFile,"*\n* --- BREM --- Bremsstrahlung by e+/- and muons/hadrons. Flag = %5d, BCUTE = %13.4g, BCUTM = %13.4g \n",
243 fProcessFlag[kBREM], fCutValue[kBCUTE], fCutValue[kBCUTM]);
244
245 //
246 // e+/- -> e+/- gamma
247 //
248 Float_t cutB = fCutValue[kBCUTE];
249 if (fCutValue[kBCUTE] == -1.) cutB = fgDCutValue[kBCUTE];
250
251
37656ca5 252 if (TMath::Abs(fProcessFlag[kBREM]) > 0) {
fb2cbbec 253 fprintf(fgFile,"EMFCUT %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fELPO-THR\n",cutB, 0., 0., fCMatMin, fCMatMax, 1.);
254 } else {
255 fprintf(fgFile,"EMFCUT %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fELPO-THR\n",1.e10, 0., 0., fCMatMin, fCMatMax, 1.);
256 }
257
258 //
259 // Bremsstrahlung by muons and hadrons
260 //
261 cutB = fCutValue[kBCUTM];
262 if (fCutValue[kBCUTM] == -1.) cutB = fgDCutValue[kBCUTM];
263 if (fProcessFlag[kBREM] == 2) cutB = -1.;
264 Float_t flag = 2.;
265 if (fProcessFlag[kBREM] == 0) flag = -2.;
266
267 fprintf(fgFile,"PAIRBREM %10.1f%10.4g%10.4g%10.1f%10.1f\n", flag, 0., cutB, fCMatMin, fCMatMax);
268}
269
270void TFlukaConfigOption::ProcessCOMP()
271{
272 // Process COMP option
273 fprintf(fgFile,"*\n* --- COMP --- Compton scattering Flag = %5d \n", fProcessFlag[kCOMP]);
274
275 //
276 // Compton scattering
277 //
278
279 if (fProcessFlag[kCOMP] > 0) {
280 fprintf(fgFile,"EMFCUT %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",0. , 0., 0., fCMatMin, fCMatMax, 1.);
281 } else {
282 fprintf(fgFile,"EMFCUT %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",1.e10, 0., 0., fCMatMin, fCMatMax, 1.);
283 }
284}
285
286void TFlukaConfigOption::ProcessPHOT()
287{
288 // Process PHOS option
289 fprintf(fgFile,"*\n* --- PHOT --- Photoelectric effect. Flag = %5d\n", fProcessFlag[kPHOT]);
290
291 //
292 // Photoelectric effect
293 //
294
295 if (fProcessFlag[kPHOT] > 0) {
81f1d030 296 fprintf(fgFile,"EMFCUT %10.4g%10.4g%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",0., 0., 0., fCMatMin, fCMatMax, 1.);
fb2cbbec 297 } else {
298 fprintf(fgFile,"EMFCUT %10.1f%10.4g%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",0., 1.e10, 0., fCMatMin, fCMatMax, 1.);
299 }
300}
301
302void TFlukaConfigOption::ProcessANNI()
303{
304 // Process ANNI option
305 fprintf(fgFile,"*\n* --- ANNI --- Positron annihilation. Flag = %5d \n", fProcessFlag[kANNI]);
306
307 //
308 // Positron annihilation
309 //
310
311 if (fProcessFlag[kANNI] > 0) {
312 fprintf(fgFile,"EMFCUT %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fANNH-THR\n",0. , 0., 0., fCMatMin, fCMatMax, 1.);
313 } else {
314 fprintf(fgFile,"EMFCUT %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fANNH-THR\n",1.e10, 0., 0., fCMatMin, fCMatMax, 1.);
315 }
316}
317
318
319void TFlukaConfigOption::ProcessPFIS()
320{
321 // Process PFIS option
322 fprintf(fgFile,"*\n* --- PFIS --- Photonuclear interaction Flag = %5d\n", fProcessFlag[kPFIS]);
323
324 //
325 // Photonuclear interactions
326 //
327
328 if (fProcessFlag[kPFIS] > 0) {
329 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.);
330 } else {
331 fprintf(fgFile,"PHOTONUC %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",-1. , 0., 0., fCMatMin, fCMatMax, 1.);
332 }
333}
334
335void TFlukaConfigOption::ProcessMUNU()
336{
337 // Process MUNU option
338 fprintf(fgFile,"*\n* --- MUNU --- Muon nuclear interaction. Flag = %5d\n", fProcessFlag[kMUNU]);
339
340 //
341 // Muon nuclear interactions
342 //
343 if (fProcessFlag[kMUNU] > 0) {
344 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.);
345 } else {
346 fprintf(fgFile,"MUPHOTON %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",-1. , 0., 0., fCMatMin, fCMatMax, 1.);
347 }
348}
349
350void TFlukaConfigOption::ProcessRAYL()
351{
352 // Process RAYL option
353 fprintf(fgFile,"*\n* --- RAYL --- Rayleigh Scattering. Flag = %5d\n", fProcessFlag[kRAYL]);
354
355 //
356 // Rayleigh scattering
357 //
358 Int_t nreg;
359 Int_t* reglist = fgGeom->GetMaterialList(fMedium, nreg);
360 //
361 // Loop over regions of a given material
362 for (Int_t k = 0; k < nreg; k++) {
363 Float_t ireg = reglist[k];
364 if (fProcessFlag[kRAYL] > 0) {
365 fprintf(fgFile,"EMFRAY %10.1f%10.1f%10.1f%10.1f\n", 1., ireg, ireg, 1.);
366 } else {
367 fprintf(fgFile,"EMFRAY %10.1f%10.1f%10.1f%10.1f\n", 3., ireg, ireg, 1.);
368 }
369 }
370}
371
372void TFlukaConfigOption::ProcessCKOV()
373{
374 // Process CKOV option
375 fprintf(fgFile,"*\n* --- CKOV --- Cerenkov Photon production. %5d\n", fProcessFlag[kCKOV]);
376
377 //
378 // Cerenkov photon production
379 //
380
381 TFluka* fluka = (TFluka*) gMC;
382 TObjArray *matList = fluka->GetFlukaMaterials();
383 Int_t nmaterial = matList->GetEntriesFast();
384 for (Int_t im = 0; im < nmaterial; im++)
385 {
386 TGeoMaterial* material = dynamic_cast<TGeoMaterial*> (matList->At(im));
387 Int_t idmat = material->GetIndex();
388//
389// Check if global option
390 if (fMedium != -1 && idmat != fMedium) continue;
391
392 TFlukaCerenkov* cerenkovProp;
393 if (!(cerenkovProp = dynamic_cast<TFlukaCerenkov*>(material->GetCerenkovProperties()))) continue;
394 //
395 // This medium has Cerenkov properties
396 //
397 //
398 if (fMedium == -1 || (fMedium != -1 && fProcessFlag[kCKOV] > 0)) {
399 // Write OPT-PROD card for each medium
400 Float_t emin = cerenkovProp->GetMinimumEnergy();
401 Float_t emax = cerenkovProp->GetMaximumEnergy();
402 fprintf(fgFile, "OPT-PROD %10.4g%10.4g%10.4g%10.4g%10.4g%10.4gCERENKOV\n", emin, emax, 0.,
403 Float_t(idmat), Float_t(idmat), 0.);
404 //
405 // Write OPT-PROP card for each medium
406 // Forcing FLUKA to call user routines (queffc.cxx, rflctv.cxx, rfrndx.cxx)
407 //
408 fprintf(fgFile, "OPT-PROP %10.4g%10.4g%10.4g%10.1f%10.1f%10.1fWV-LIMIT\n",
409 cerenkovProp->GetMinimumWavelength(), cerenkovProp->GetMaximumWavelength(), cerenkovProp->GetMaximumWavelength(),
410 Float_t(idmat), Float_t(idmat), 0.0);
411
2906553d 412
413 fprintf(fgFile, "OPT-PROP %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n", -100., -100., -100.,
414 Float_t(idmat), Float_t(idmat), 0.0);
fb2cbbec 415
416 for (Int_t j = 0; j < 3; j++) {
417 fprintf(fgFile, "OPT-PROP %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f&\n", -100., -100., -100.,
418 Float_t(idmat), Float_t(idmat), 0.0);
419 }
2906553d 420
421
422 // Photon detection efficiency user defined
fb2cbbec 423 if (cerenkovProp->IsSensitive())
424 fprintf(fgFile, "OPT-PROP %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fSENSITIV\n", -100., -100., -100.,
425 Float_t(idmat), Float_t(idmat), 0.0);
2906553d 426 // Material has a reflective surface
427 if (cerenkovProp->IsMetal())
428 fprintf(fgFile, "OPT-PROP %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fMETAL\n", -100., -100., -100.,
429 Float_t(idmat), Float_t(idmat), 0.0);
430
fb2cbbec 431 } else {
432 fprintf(fgFile,"OPT-PROD %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fCERE-OFF\n",0., 0., 0., fCMatMin, fCMatMax, 1.);
433 }
434 }
435}
436
437
438void TFlukaConfigOption::ProcessHADR()
439{
440 // Process HADR option
441 fprintf(fgFile,"*\n* --- HADR --- Hadronic interactions. Flag = %5d\n", fProcessFlag[kHADR]);
442
443 if (fProcessFlag[kHADR] > 0) {
444 fprintf(fgFile,"*\n*Hadronic interaction is ON by default in FLUKA\n");
445 } else {
446 if (fMedium != -1) printf("Hadronic interactions cannot be switched off material by material !\n");
447 fprintf(fgFile,"THRESHOL %10.1f%10.1f%10.1f%10.1e%10.1f\n",0., 0., 0., 1.e10, 0.);
448 }
449}
450
451
452
453void TFlukaConfigOption::ProcessMULS()
454{
455 // Process MULS option
456 fprintf(fgFile,"*\n* --- MULS --- Muliple Scattering. Flag = %5d\n", fProcessFlag[kMULS]);
457 //
458 // Multiple scattering
459 //
460 if (fProcessFlag[kMULS] > 0) {
461 fprintf(fgFile,"*\n*Multiple scattering is ON by default in FLUKA\n");
462 } else {
463 fprintf(fgFile,"MULSOPT %10.1f%10.1f%10.1f%10.1f%10.1f\n",0., 3., 3., fCMatMin, fCMatMax);
464 }
465}
466
467void TFlukaConfigOption::ProcessLOSS()
468{
469 // Process LOSS option
470 fprintf(fgFile,"*\n* --- LOSS --- Ionisation energy loss. Flags: LOSS = %5d, DRAY = %5d, STRA = %5d; Cuts: DCUTE = %13.4g, DCUTM = %13.4g \n",
471 fProcessFlag[kLOSS], fProcessFlag[kDRAY], fProcessFlag[kSTRA], fCutValue[kDCUTE], fCutValue[kDCUTM]);
472 //
473 // Ionisation energy loss
474 //
475 //
476 // Impose consistency
477
478 if (fProcessFlag[kLOSS] == 1 || fProcessFlag[kLOSS] == 3) {
479 fProcessFlag[kDRAY] = 1;
480 } else if (fProcessFlag[kLOSS] == 2) {
481 fProcessFlag[kDRAY] = 0;
482 fCutValue[kDCUTE] = 1.e10;
483 fCutValue[kDCUTM] = 1.e10;
484 } else {
485 if (fProcessFlag[kDRAY] == 1) {
486 fProcessFlag[kLOSS] = 1;
487 } else if (fProcessFlag[kDRAY] == 0) {
488 fProcessFlag[kLOSS] = 2;
489 fCutValue[kDCUTE] = 1.e10;
490 fCutValue[kDCUTM] = 1.e10;
491 }
492 }
493
494 if (fCutValue[kDCUTE] == -1.) fCutValue[kDCUTE] = fgDCutValue[kDCUTE];
495 if (fCutValue[kDCUTM] == -1.) fCutValue[kDCUTM] = fgDCutValue[kDCUTM];
496
497 Float_t cutM = fCutValue[kDCUTM];
498
499
500 if (fProcessFlag[kSTRA] == -1) fProcessFlag[kSTRA] = fgDProcessFlag[kSTRA];
501 if (fProcessFlag[kSTRA] == 0) fProcessFlag[kSTRA] = 1;
502 Float_t stra = (Float_t) fProcessFlag[kSTRA];
503
504
505 if (fProcessFlag[kLOSS] == 1 || fProcessFlag[kLOSS] == 3) {
506//
507// Restricted energy loss fluctuations
508//
509 fprintf(fgFile,"IONFLUCT %10.1f%10.1f%10.1f%10.1f%10.1f\n", 1., 1., stra, fCMatMin, fCMatMax);
2906553d 510 fprintf(fgFile,"DELTARAY %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n", cutM, 0., 0., fCMatMin, fCMatMax, 1.);
fb2cbbec 511 } else if (fProcessFlag[kLOSS] == 4) {
512//
513// No fluctuations
514//
515 fprintf(fgFile,"IONFLUCT %10.1f%10.1f%10.1f%10.1f%10.1f\n",-1., -1., stra, fCMatMin, fCMatMax);
516 fprintf(fgFile,"DELTARAY %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n", 1.e10, 0., 0., fCMatMin, fCMatMax, 1.);
517 } else {
518//
519// Full fluctuations
520//
521 fprintf(fgFile,"IONFLUCT %10.1f%10.1f%10.1f%10.1f%10.1f\n",1., 1., stra, fCMatMin, fCMatMax);
522 fprintf(fgFile,"DELTARAY %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n", 1.e10, 0., 0., fCMatMin, fCMatMax, 1.);
523 }
524}
525
526
527void TFlukaConfigOption::ProcessCUTGAM()
528{
529// Cut on gammas
530//
531 fprintf(fgFile,"*\n*Cut for Gammas. CUTGAM = %13.4g\n", fCutValue[kCUTGAM]);
532 if (fMedium == -1) {
533 fprintf(fgFile,"EMFCUT %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n",
534 0., fCutValue[kCUTGAM], 0., 0., Float_t(fgGeom->NofVolumes()), 1.);
46bee1bf 535
fb2cbbec 536 } else {
537 Int_t nreg, *reglist;
538 Float_t ireg;
539 reglist = fgGeom->GetMaterialList(fMedium, nreg);
540 // Loop over regions of a given material
541 for (Int_t k = 0; k < nreg; k++) {
542 ireg = reglist[k];
543 fprintf(fgFile,"EMFCUT %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", 0.,fCutValue[kCUTGAM], 0., ireg, ireg, 1.);
544 }
545 }
46bee1bf 546 fprintf(fgFile,"EMFCUT %10.4g%10.4g%10.1f%10.1f%10.1f%10.1fPROD-CUT\n",
37656ca5 547 0., fCutValue[kCUTGAM], 1., fCMatMin, fCMatMax, 1.);
fb2cbbec 548}
549
550void TFlukaConfigOption::ProcessCUTELE()
551{
552// Cut on e+/e-
553//
554 fprintf(fgFile,"*\n*Cut for e+/e-. CUTELE = %13.4g\n", fCutValue[kCUTELE]);
555 if (fMedium == -1) {
556 fprintf(fgFile,"EMFCUT %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n",
557 -fCutValue[kCUTELE], 0., 0., 0., Float_t(fgGeom->NofVolumes()), 1.);
558 } else {
559 Int_t nreg, *reglist;
560 Float_t ireg;
561 reglist = fgGeom->GetMaterialList(fMedium, nreg);
562 // Loop over regions of a given material
563 for (Int_t k = 0; k < nreg; k++) {
564 ireg = reglist[k];
565 fprintf(fgFile,"EMFCUT %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", -fCutValue[kCUTELE], 0., 0., ireg, ireg, 1.);
566 }
567 }
46bee1bf 568 fprintf(fgFile,"EMFCUT %10.4g%10.4g%10.1f%10.1f%10.1f%10.1fPROD-CUT\n",
37656ca5 569 -fCutValue[kCUTELE], 0., 1., fCMatMin, fCMatMax, 1.);
fb2cbbec 570}
571
572void TFlukaConfigOption::ProcessCUTNEU()
573{
574 // Cut on neutral hadrons
2906553d 575 fprintf(fgFile,"*\n*Cut for neutral hadrons. CUTNEU = %13.4g\n", fCutValue[kCUTNEU]);
fb2cbbec 576 if (fMedium == -1) {
577 Float_t cut = fCutValue[kCUTNEU];
2906553d 578 //
fb2cbbec 579 // 8.0 = Neutron
580 // 9.0 = Antineutron
2906553d 581 //
582 // If the cut is > 19.6 MeV it is assumed the low energy neutron transport is requested.
583 // In this case the cut has to coincide with the upper limit of the first energy group.
584 //
585 Float_t neutronCut = cut;
586 if (neutronCut < 0.0196) {
587 neutronCut = 0.0196;
89aa6d28 588 printf("Cut on neutron lower than upper limit of first energy group.\n");
2906553d 589 printf("Cut reset to 19.6 MeV !\n");
590 }
591 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -neutronCut, 8.0, 9.0);
592 //
fb2cbbec 593 // 12.0 = Kaon zero long
594 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 12.0, 12.0);
595 // 17.0 = Lambda, 18.0 = Antilambda
596 // 19.0 = Kaon zero short
597 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 17.0, 19.0);
598 // 22.0 = Sigma zero, Pion zero, Kaon zero
599 // 25.0 = Antikaon zero
600 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 22.0, 25.0);
601 // 32.0 = Antisigma zero
602 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 32.0, 32.0);
603 // 34.0 = Xi zero
604 // 35.0 = AntiXi zero
605 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 34.0, 35.0);
606 // 47.0 = D zero
607 // 48.0 = AntiD zero
608 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 47.0, 48.0);
609 // 53.0 = Xi_c zero
610 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 53.0, 53.0);
611 // 55.0 = Xi'_c zero
612 // 56.0 = Omega_c zero
613 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 55.0, 56.0);
614 // 59.0 = AntiXi_c zero
615 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 59.0, 59.0);
616 // 61.0 = AntiXi'_c zero
617 // 62.0 = AntiOmega_c zero
618 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 61.0, 62.0);
619 } else {
620 printf("Cuts on neutral hadrons material by material not yet implemented !\n");
621 }
622}
623
624void TFlukaConfigOption::ProcessCUTHAD()
625{
626 // Cut on charged hadrons
627 fprintf(fgFile,"*\n*Cut for charge hadrons. CUTHAD = %13.4g\n", fCutValue[kCUTHAD]);
628 if (fMedium == -1) {
629 Float_t cut = fCutValue[kCUTHAD];
630 // 1.0 = Proton
631 // 2.0 = Antiproton
632 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 1.0, 2.0);
633 // 13.0 = Positive Pion, Negative Pion, Positive Kaon
634 // 16.0 = Negative Kaon
635 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 13.0, 16.0);
636 // 20.0 = Negative Sigma
637 // 21.0 = Positive Sigma
638 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 20.0, 21.0);
639 // 31.0 = Antisigma minus
640 // 33.0 = Antisigma plus
641 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 31.0, 31.0);
642 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 33.0, 33.0);
643 // 36.0 = Negative Xi, Positive Xi, Omega minus
644 // 39.0 = Antiomega
645 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 36.0, 39.0);
646 // 45.0 = D plus
647 // 46.0 = D minus
648 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 45.0, 46.0);
649 // 49.0 = D_s plus, D_s minus, Lambda_c plus
650 // 52.0 = Xi_c plus
651 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 49.0, 52.0);
652 // 54.0 = Xi'_c plus
653 // 60.0 = AntiXi'_c minus
654 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 54.0, 54.0);
655 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 60.0, 60.0);
656 // 57.0 = Antilambda_c minus
657 // 58.0 = AntiXi_c minus
658 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 57.0, 58.0);
659 } else {
660 printf("Cuts on charged hadrons material by material not yet implemented !\n");
661 }
662}
663
664void TFlukaConfigOption::ProcessCUTMUO()
665{
666 // Cut on muons
667 fprintf(fgFile,"*\n*Cut for muons. CUTMUO = %13.4g\n", fCutValue[kCUTMUO]);
668 Float_t cut = fCutValue[kCUTMUO];
669 if (fMedium == -1) {
670 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n",-cut, 10.0, 11.0);
671 } else {
672 printf("Cuts on muons material by material not yet implemented !\n");
673 }
674
675
676}
677
678void TFlukaConfigOption::ProcessTOFMAX()
679{
680 // Cut time of flight
681 Float_t cut = fCutValue[kTOFMAX];
682 fprintf(fgFile,"*\n*Cut on time of flight. TOFMAX = %13.4g\n", fCutValue[kTOFMAX]);
683 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 684}