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