- Correct card for optical properties diielectroc-METAL boundary
[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
386             fprintf(fgFile, "OPT-PROP  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n", -100., -100., -100., 
387                     Float_t(idmat), Float_t(idmat), 0.0);
388             
389             for (Int_t j = 0; j < 3; j++) {
390                 fprintf(fgFile, "OPT-PROP  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f&\n", -100., -100., -100., 
391                         Float_t(idmat), Float_t(idmat), 0.0);
392             }
393
394
395             // Photon detection efficiency user defined     
396             if (cerenkovProp->IsSensitive())
397                 fprintf(fgFile, "OPT-PROP  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fSENSITIV\n", -100., -100., -100., 
398                         Float_t(idmat), Float_t(idmat), 0.0);
399             // Material has a reflective surface
400             if (cerenkovProp->IsMetal())
401                 fprintf(fgFile, "OPT-PROP  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fMETAL\n", -100., -100., -100., 
402                         Float_t(idmat), Float_t(idmat), 0.0);
403
404         } else {
405             fprintf(fgFile,"OPT-PROD  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fCERE-OFF\n",0., 0., 0., fCMatMin, fCMatMax, 1.);
406         }
407     }
408 }
409
410
411 void TFlukaConfigOption::ProcessHADR()
412 {
413     // Process HADR option
414     fprintf(fgFile,"*\n* --- HADR --- Hadronic interactions. Flag = %5d\n", fProcessFlag[kHADR]);
415     
416     if (fProcessFlag[kHADR] > 0) {
417         fprintf(fgFile,"*\n*Hadronic interaction is ON by default in FLUKA\n");
418     } else {
419         if (fMedium != -1) printf("Hadronic interactions cannot be switched off material by material !\n");
420         fprintf(fgFile,"THRESHOL  %10.1f%10.1f%10.1f%10.1e%10.1f\n",0., 0., 0., 1.e10, 0.);
421     }
422 }
423
424
425
426 void TFlukaConfigOption::ProcessMULS()
427 {
428     // Process MULS option
429     fprintf(fgFile,"*\n* --- MULS --- Muliple Scattering. Flag = %5d\n", fProcessFlag[kMULS]);
430     //
431     // Multiple scattering
432     //
433     if (fProcessFlag[kMULS] > 0) {
434         fprintf(fgFile,"*\n*Multiple scattering is  ON by default in FLUKA\n");
435     } else {
436         fprintf(fgFile,"MULSOPT   %10.1f%10.1f%10.1f%10.1f%10.1f\n",0., 3., 3., fCMatMin, fCMatMax);
437     }
438 }
439
440 void TFlukaConfigOption::ProcessLOSS()
441 {
442     // Process LOSS option
443     fprintf(fgFile,"*\n* --- LOSS --- Ionisation energy loss. Flags: LOSS = %5d, DRAY = %5d, STRA = %5d; Cuts: DCUTE = %13.4g, DCUTM = %13.4g \n",
444             fProcessFlag[kLOSS], fProcessFlag[kDRAY], fProcessFlag[kSTRA], fCutValue[kDCUTE], fCutValue[kDCUTM]);
445     //
446     // Ionisation energy loss
447     //
448     //
449     // Impose consistency
450     
451     if (fProcessFlag[kLOSS] == 1 || fProcessFlag[kLOSS] == 3) {
452         fProcessFlag[kDRAY] = 1;
453     } else if (fProcessFlag[kLOSS] == 2) {
454         fProcessFlag[kDRAY] = 0;
455         fCutValue[kDCUTE] = 1.e10;
456         fCutValue[kDCUTM] = 1.e10;      
457     } else {
458         if (fProcessFlag[kDRAY] == 1) {
459             fProcessFlag[kLOSS] = 1;
460         } else if (fProcessFlag[kDRAY] == 0) {
461             fProcessFlag[kLOSS] = 2;
462             fCutValue[kDCUTE] = 1.e10;
463             fCutValue[kDCUTM] = 1.e10;  
464         }
465     }
466     
467     if (fCutValue[kDCUTE] == -1.) fCutValue[kDCUTE] = fgDCutValue[kDCUTE];
468     if (fCutValue[kDCUTM] == -1.) fCutValue[kDCUTM] = fgDCutValue[kDCUTM];    
469     
470     Float_t cutM =  fCutValue[kDCUTM];    
471         
472
473     if (fProcessFlag[kSTRA] == -1) fProcessFlag[kSTRA] = fgDProcessFlag[kSTRA];
474     if (fProcessFlag[kSTRA] ==  0) fProcessFlag[kSTRA] = 1;
475     Float_t stra = (Float_t) fProcessFlag[kSTRA];
476     
477
478     if (fProcessFlag[kLOSS] == 1 || fProcessFlag[kLOSS] == 3) {
479 //
480 // Restricted energy loss fluctuations 
481 //
482         fprintf(fgFile,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f\n", 1., 1., stra, fCMatMin, fCMatMax);
483         fprintf(fgFile,"DELTARAY  %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n", cutM, 0., 0., fCMatMin, fCMatMax, 1.);
484     } else if (fProcessFlag[kLOSS] == 4) {
485 //
486 // No fluctuations
487 //
488         fprintf(fgFile,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f\n",-1., -1., stra, fCMatMin, fCMatMax);        
489         fprintf(fgFile,"DELTARAY  %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n", 1.e10, 0., 0., fCMatMin, fCMatMax, 1.);      
490     }  else {
491 //
492 // Full fluctuations
493 //
494         fprintf(fgFile,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f\n",1., 1., stra, fCMatMin, fCMatMax);  
495         fprintf(fgFile,"DELTARAY  %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n", 1.e10, 0., 0., fCMatMin, fCMatMax, 1.);      
496     }
497 }
498
499
500 void TFlukaConfigOption::ProcessCUTGAM()
501 {
502 // Cut on gammas
503 //
504     fprintf(fgFile,"*\n*Cut for Gammas. CUTGAM = %13.4g\n", fCutValue[kCUTGAM]);
505     if (fMedium == -1) {
506         fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", 
507                 0., fCutValue[kCUTGAM], 0., 0., Float_t(fgGeom->NofVolumes()), 1.);
508     } else {
509         Int_t nreg, *reglist;
510         Float_t ireg;
511         reglist = fgGeom->GetMaterialList(fMedium, nreg);
512         // Loop over regions of a given material
513         for (Int_t k = 0; k < nreg; k++) {
514             ireg = reglist[k];
515             fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", 0.,fCutValue[kCUTGAM], 0., ireg, ireg, 1.);
516         }
517     }
518 }
519
520 void TFlukaConfigOption::ProcessCUTELE()
521 {
522 // Cut on e+/e-
523 //
524     fprintf(fgFile,"*\n*Cut for e+/e-. CUTELE = %13.4g\n", fCutValue[kCUTELE]);
525     if (fMedium == -1) {
526         fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", 
527                 -fCutValue[kCUTELE], 0., 0., 0., Float_t(fgGeom->NofVolumes()), 1.);
528     } else {
529         Int_t nreg, *reglist;
530         Float_t ireg;
531         reglist = fgGeom->GetMaterialList(fMedium, nreg);
532         // Loop over regions of a given material
533         for (Int_t k = 0; k < nreg; k++) {
534             ireg = reglist[k];
535             fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", -fCutValue[kCUTELE], 0., 0., ireg, ireg, 1.);
536         }
537     }
538 }
539
540 void TFlukaConfigOption::ProcessCUTNEU()
541 {
542     // Cut on neutral hadrons
543     fprintf(fgFile,"*\n*Cut for neutral hadrons. CUTNEU = %13.4g\n", fCutValue[kCUTNEU]);
544     if (fMedium == -1) {
545         Float_t cut = fCutValue[kCUTNEU];
546         //
547         // 8.0 = Neutron
548         // 9.0 = Antineutron
549         //
550         // If the cut is > 19.6 MeV it is assumed the low energy neutron transport is requested.
551         // In this case the cut has to coincide with the upper  limit of the first energy group.
552         //
553         Float_t neutronCut = cut;
554         if (neutronCut < 0.0196) {
555             neutronCut = 0.0196;
556             printf("Cut on neutron lower than upper limit if first energy group.\n");
557             printf("Cut reset to 19.6 MeV !\n");
558         }
559         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -neutronCut,  8.0,  9.0);
560         //
561         // 12.0 = Kaon zero long
562         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 12.0, 12.0);
563         // 17.0 = Lambda, 18.0 = Antilambda
564         // 19.0 = Kaon zero short
565         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 17.0, 19.0);
566         // 22.0 = Sigma zero, Pion zero, Kaon zero
567         // 25.0 = Antikaon zero
568         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 22.0, 25.0);
569         // 32.0 = Antisigma zero
570         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 32.0, 32.0);
571         // 34.0 = Xi zero
572         // 35.0 = AntiXi zero
573         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 34.0, 35.0);
574         // 47.0 = D zero
575         // 48.0 = AntiD zero
576         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 47.0, 48.0);
577         // 53.0 = Xi_c zero
578         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 53.0, 53.0);
579         // 55.0 = Xi'_c zero
580         // 56.0 = Omega_c zero
581         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 55.0, 56.0);
582         // 59.0 = AntiXi_c zero
583         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 59.0, 59.0);
584         // 61.0 = AntiXi'_c zero
585         // 62.0 = AntiOmega_c zero
586         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 61.0, 62.0);
587     } else {
588         printf("Cuts on neutral hadrons material by material not yet implemented !\n");
589     }
590 }
591
592 void TFlukaConfigOption::ProcessCUTHAD()
593 {
594     // Cut on charged hadrons
595     fprintf(fgFile,"*\n*Cut for charge hadrons. CUTHAD = %13.4g\n", fCutValue[kCUTHAD]);
596     if (fMedium == -1) {
597         Float_t cut = fCutValue[kCUTHAD];
598         // 1.0 = Proton
599         // 2.0 = Antiproton
600         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut,  1.0,  2.0);
601         // 13.0 = Positive Pion, Negative Pion, Positive Kaon
602         // 16.0 = Negative Kaon
603         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 13.0, 16.0);
604         // 20.0 = Negative Sigma
605         // 21.0 = Positive Sigma
606         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 20.0, 21.0);
607         // 31.0 = Antisigma minus
608         // 33.0 = Antisigma plus
609         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 31.0, 31.0);
610         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 33.0, 33.0);
611         // 36.0 = Negative Xi, Positive Xi, Omega minus
612         // 39.0 = Antiomega
613         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 36.0, 39.0);
614         // 45.0 = D plus
615         // 46.0 = D minus
616         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 45.0, 46.0);
617         // 49.0 = D_s plus, D_s minus, Lambda_c plus
618         // 52.0 = Xi_c plus
619         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 49.0, 52.0);
620         // 54.0 = Xi'_c plus
621         // 60.0 = AntiXi'_c minus
622         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 54.0, 54.0);
623         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 60.0, 60.0);
624         // 57.0 = Antilambda_c minus
625         // 58.0 = AntiXi_c minus
626         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 57.0, 58.0);
627     } else {
628         printf("Cuts on charged hadrons material by material not yet implemented !\n");
629     }
630 }
631
632 void TFlukaConfigOption::ProcessCUTMUO()
633 {
634     // Cut on muons
635     fprintf(fgFile,"*\n*Cut for muons. CUTMUO = %13.4g\n", fCutValue[kCUTMUO]);
636     Float_t cut = fCutValue[kCUTMUO];
637     if (fMedium == -1) {
638         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n",-cut, 10.0, 11.0);
639     } else {
640         printf("Cuts on muons material by material not yet implemented !\n");
641     }
642     
643     
644 }
645
646 void TFlukaConfigOption::ProcessTOFMAX()
647 {
648     // Cut time of flight
649     Float_t cut = fCutValue[kTOFMAX];
650     fprintf(fgFile,"*\n*Cut on time of flight. TOFMAX = %13.4g\n", fCutValue[kTOFMAX]);
651     fprintf(fgFile,"TIME-CUT  %10.4g%10.1f%10.1f%10.1f%10.1f\n",cut*1.e9,0.,0.,-6.0,64.0);
652 }