Neutron cuts according to FLUKA groups. (Ernesto Lopez Torres)
[u/mrichter/AliRoot.git] / TFluka / TFlukaConfigOption.cxx
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"
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
28 Float_t            TFlukaConfigOption::fgMatMin(-1.);
29 Float_t            TFlukaConfigOption::fgMatMax(-1.);
30 FILE*              TFlukaConfigOption::fgFile(0x0);
31 TFlukaMCGeometry*  TFlukaConfigOption::fgGeom(0x0);
32
33 Double_t TFlukaConfigOption::fgDCutValue[11];
34 Int_t    TFlukaConfigOption::fgDProcessFlag[15];
35
36
37 ClassImp(TFlukaConfigOption)
38
39
40 TFlukaConfigOption::TFlukaConfigOption()
41 {
42     // Default constructor
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;
49 }
50
51
52 TFlukaConfigOption::TFlukaConfigOption(Int_t medium)
53 {
54     // Constructor
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;
61 }
62
63 void TFlukaConfigOption::SetCut(const char* flagname, Double_t val)
64 {
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     }
76 }
77
78 void 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"};
84     
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 }
94
95 void TFlukaConfigOption::WriteFlukaInputCards()
96 {
97     // Write the FLUKA input cards for the set of process flags and cuts
98     //
99     //
100     if (fMedium != -1) {
101         TFluka* fluka = (TFluka*) gMC;
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
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());
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();
164 }
165
166 void TFlukaConfigOption::ProcessDCAY()
167 {
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
178 void TFlukaConfigOption::ProcessPAIR()
179 {
180     // Process PAIR option
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]);
183     //
184     // gamma -> e+ e-
185     //
186     if (fProcessFlag[kPAIR] > 0) {
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.);
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
209     fProcessFlag[kBREM] = - fProcessFlag[kBREM];
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
239 void 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     
252     if (TMath::Abs(fProcessFlag[kBREM]) > 0) {
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
270 void 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
286 void 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) {
296         fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",0., 0., 0.,  fCMatMin, fCMatMax, 1.);
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
302 void 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
319 void 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
335 void 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
350 void 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
372 void 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             
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);
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             }
420
421
422             // Photon detection efficiency user defined     
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);
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
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
438 void 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
453 void 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
467 void 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);
510         fprintf(fgFile,"DELTARAY  %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n", cutM, 0., 0., fCMatMin, fCMatMax, 1.);
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
527 void 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.);
535
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     }
546     fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1fPROD-CUT\n", 
547             0., fCutValue[kCUTGAM], 1., fCMatMin, fCMatMax, 1.);
548 }
549
550 void 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     }
568     fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1fPROD-CUT\n", 
569             -fCutValue[kCUTELE], 0., 1., fCMatMin, fCMatMax, 1.);
570 }
571
572 void TFlukaConfigOption::ProcessCUTNEU()
573 {
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         //
619         //
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 {
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");
658     }
659 }
660
661 void 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 {
697       printf("Cuts on charged hadrons material by material not yet implemented !\n"); 
698     }
699 }
700
701 void 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
715 void 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);
721 }