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