]> 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
Small bug fixed for some class of tracks (Sacha)
[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{
501e3a8d 574 // Cut on neutral hadrons
575 fprintf(fgFile,"*\n*Cut for neutral hadrons. CUTNEU = %13.4g\n", fCutValue[kCUTNEU]);
576
577 // Cut on neutral hadrons
578 fprintf(fgFile,"*\n*Cut for neutral hadrons. CUTNEU = %13.4g\n", fCutValue[kCUTNEU]);
579
580 // Energy group structure of the 72-neutron ENEA data set:
581 const Float_t neutronGroupUpLimit[72] = {
582 1.9600E-02, 1.7500E-02, 1.4918E-02, 1.3499E-02, 1.2214E-02, 1.1052E-02, 1.0000E-02, 9.0484E-03,
583 8.1873E-03, 7.4082E-03, 6.7032E-03, 6.0653E-03, 5.4881E-03, 4.9659E-03, 4.4933E-03, 4.0657E-03,
584 3.6788E-03, 3.3287E-03, 3.0119E-03, 2.7253E-03, 2.4660E-03, 2.2313E-03, 2.0190E-03, 1.8268E-03,
585 1.6530E-03, 1.4957E-03, 1.3534E-03, 1.2246E-03, 1.1080E-03, 1.0026E-03, 9.0718E-04, 8.2085E-04,
586 7.4274E-04, 6.0810E-04, 4.9787E-04, 4.0762E-04, 3.3373E-04, 2.7324E-04, 2.2371E-04, 1.8316E-04,
587 1.4996E-04, 1.2277E-04, 8.6517E-05, 5.2475E-05, 3.1828E-05, 2.1852E-05, 1.5034E-05, 1.0332E-05,
588 7.1018E-06, 4.8809E-06, 3.3546E-06, 2.3054E-06, 1.5846E-06, 1.0446E-06, 6.8871E-07, 4.5400E-07,
589 2.7537E-07, 1.6702E-07, 1.0130E-07, 6.1442E-08, 3.7267E-08, 2.2603E-08, 1.5535E-08, 1.0677E-08,
590 7.3375E-09, 5.0435E-09, 3.4662E-09, 2.3824E-09, 1.6374E-09, 1.1254E-09, 6.8257E-10, 4.1400E-10
591 };
592
593 Float_t cut = fCutValue[kCUTNEU];
594 //
595 // 8.0 = Neutron
596 // 9.0 = Antineutron
597 // Find the FLUKA neutron group corresponding to the cut
598 //
599 Float_t neutronCut = cut;
600 Int_t groupCut = 1; // if cut is > 19.6 MeV not low energy neutron transport is done
601 if (neutronCut < 0.0196) {
602 neutronCut = 0.0196;
603 // Search the group cutoff for the energy cut
604 Int_t i;
605 for( i=0; i<72; i++ ) {
606 if (cut > neutronGroupUpLimit[i]) break;
607 }
608 groupCut = i+1;
609 }
610
611 if (fMedium == -1) {
612 Float_t cut = fCutValue[kCUTNEU];
613 // 8.0 = Neutron
614 // 9.0 = Antineutron
615 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -neutronCut, 8.0, 9.0);
616 fprintf(fgFile,"LOW-BIAS %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n",
617 Float_t(groupCut), 73.0, 0.95, 2., Float_t(fgGeom->NofVolumes()), 1.);
618 //
2906553d 619 //
fb2cbbec 620 // 12.0 = Kaon zero long
621 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 12.0, 12.0);
622 // 17.0 = Lambda, 18.0 = Antilambda
623 // 19.0 = Kaon zero short
624 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 17.0, 19.0);
625 // 22.0 = Sigma zero, Pion zero, Kaon zero
626 // 25.0 = Antikaon zero
627 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 22.0, 25.0);
628 // 32.0 = Antisigma zero
629 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 32.0, 32.0);
630 // 34.0 = Xi zero
631 // 35.0 = AntiXi zero
632 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 34.0, 35.0);
633 // 47.0 = D zero
634 // 48.0 = AntiD zero
635 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 47.0, 48.0);
636 // 53.0 = Xi_c zero
637 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 53.0, 53.0);
638 // 55.0 = Xi'_c zero
639 // 56.0 = Omega_c zero
640 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 55.0, 56.0);
641 // 59.0 = AntiXi_c zero
642 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 59.0, 59.0);
643 // 61.0 = AntiXi'_c zero
644 // 62.0 = AntiOmega_c zero
645 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 61.0, 62.0);
646 } else {
501e3a8d 647 Int_t nreg, *reglist;
648 Float_t ireg;
649 reglist = fgGeom->GetMaterialList(fMedium, nreg);
650 // Loop over regions of a given material
651 for (Int_t k = 0; k < nreg; k++) {
652 ireg = reglist[k];
653 fprintf(fgFile,"LOW-BIAS %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n",
654 Float_t(groupCut), 73.0, 0.95, ireg, ireg, 1.);
655 }
656
657 printf("Cuts on neutral hadrons material by material only implemented for low-energy neutrons !\n");
fb2cbbec 658 }
659}
660
661void TFlukaConfigOption::ProcessCUTHAD()
662{
663 // Cut on charged hadrons
664 fprintf(fgFile,"*\n*Cut for charge hadrons. CUTHAD = %13.4g\n", fCutValue[kCUTHAD]);
665 if (fMedium == -1) {
666 Float_t cut = fCutValue[kCUTHAD];
667 // 1.0 = Proton
668 // 2.0 = Antiproton
669 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 1.0, 2.0);
670 // 13.0 = Positive Pion, Negative Pion, Positive Kaon
671 // 16.0 = Negative Kaon
672 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 13.0, 16.0);
673 // 20.0 = Negative Sigma
674 // 21.0 = Positive Sigma
675 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 20.0, 21.0);
676 // 31.0 = Antisigma minus
677 // 33.0 = Antisigma plus
678 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 31.0, 31.0);
679 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 33.0, 33.0);
680 // 36.0 = Negative Xi, Positive Xi, Omega minus
681 // 39.0 = Antiomega
682 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 36.0, 39.0);
683 // 45.0 = D plus
684 // 46.0 = D minus
685 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 45.0, 46.0);
686 // 49.0 = D_s plus, D_s minus, Lambda_c plus
687 // 52.0 = Xi_c plus
688 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 49.0, 52.0);
689 // 54.0 = Xi'_c plus
690 // 60.0 = AntiXi'_c minus
691 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 54.0, 54.0);
692 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 60.0, 60.0);
693 // 57.0 = Antilambda_c minus
694 // 58.0 = AntiXi_c minus
695 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 57.0, 58.0);
696 } else {
501e3a8d 697 printf("Cuts on charged hadrons material by material not yet implemented !\n");
fb2cbbec 698 }
699}
700
701void TFlukaConfigOption::ProcessCUTMUO()
702{
703 // Cut on muons
704 fprintf(fgFile,"*\n*Cut for muons. CUTMUO = %13.4g\n", fCutValue[kCUTMUO]);
705 Float_t cut = fCutValue[kCUTMUO];
706 if (fMedium == -1) {
707 fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n",-cut, 10.0, 11.0);
708 } else {
709 printf("Cuts on muons material by material not yet implemented !\n");
710 }
711
712
713}
714
715void TFlukaConfigOption::ProcessTOFMAX()
716{
717 // Cut time of flight
718 Float_t cut = fCutValue[kTOFMAX];
719 fprintf(fgFile,"*\n*Cut on time of flight. TOFMAX = %13.4g\n", fCutValue[kTOFMAX]);
720 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 721}