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