- Correct setting of FUDGEM parameter.
[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 <TList.h>
25 #include <TObjArray.h>
26 #include <TVirtualMC.h>
27 #include <TGeoMaterial.h>
28 #include <TGeoMedium.h>
29 #include <TGeoManager.h>
30 #include <TGeoMedium.h>
31
32 Float_t            TFlukaConfigOption::fgMatMin(-1.);
33 Float_t            TFlukaConfigOption::fgMatMax(-1.);
34 FILE*              TFlukaConfigOption::fgFile(0x0);
35 TFlukaMCGeometry*  TFlukaConfigOption::fgGeom(0x0);
36
37 Double_t TFlukaConfigOption::fgDCutValue[11];
38 Int_t    TFlukaConfigOption::fgDProcessFlag[15];
39
40
41 ClassImp(TFlukaConfigOption)
42
43
44 TFlukaConfigOption::TFlukaConfigOption()
45 {
46     // Default constructor
47     fMedium  = -1;
48     fCMatMin = -1.;
49     fCMatMax = -1.;    
50     Int_t i;
51     for (i = 0; i < 11; i++) fCutValue[i]    = -1.;
52     for (i = 0; i < 15; i++) fProcessFlag[i] = -1;
53 }
54
55
56 TFlukaConfigOption::TFlukaConfigOption(Int_t medium)
57 {
58     // Constructor
59     fMedium = medium;
60     fCMatMin = -1.;
61     fCMatMax = -1.;    
62     Int_t i;
63     for (i = 0; i < 11; i++) fCutValue[i]    = -1.;
64     for (i = 0; i < 15; i++) fProcessFlag[i] = -1;
65 }
66
67 void TFlukaConfigOption::SetCut(const char* flagname, Double_t val)
68 {
69     // Set a cut value
70     const TString cuts[11] = 
71         {"CUTGAM", "CUTELE", "CUTNEU", "CUTHAD", "CUTMUO", "BCUTE", "BCUTM", "DCUTE", "DCUTM", "PPCUTM", "TOFMAX"};
72     Int_t i;
73     for (i = 0; i < 11; i++) {
74         if (cuts[i].CompareTo(flagname) == 0) {
75             fCutValue[i] = val;
76             if (fMedium == -1) fgDCutValue[i] = val;
77             break;
78         }
79     }
80 }
81
82 void TFlukaConfigOption::SetProcess(const char* flagname, Int_t flag)
83 {
84     // Set a process flag
85     const TString process[15] = 
86         {"DCAY", "PAIR", "COMP", "PHOT", "PFIS", "DRAY", "ANNI", "BREM", "MUNU", "CKOV", 
87          "HADR", "LOSS", "MULS", "RAYL", "STRA"};
88     
89     Int_t i;
90     for (i = 0; i < 15; i++) {
91         if (process[i].CompareTo(flagname) == 0) {
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     //
105     // Check if global option or medium specific
106
107     Bool_t mediumIsSensitive = kFALSE;
108     TGeoMedium*   med    = 0x0;
109     TGeoMedium*   medium = 0x0;    
110     TGeoMaterial* mat    = 0x0;
111
112     if (fMedium != -1) {
113         TFluka* fluka = (TFluka*) gMC;
114         fMedium = fgGeom->GetFlukaMaterial(fMedium);
115         //
116         // Check if material is actually used
117         Int_t* reglist;
118         Int_t nreg;     
119         reglist = fgGeom->GetMaterialList(fMedium, nreg);
120         if (nreg == 0) {
121             printf("Material not used !\n");
122             return;
123         }
124         //
125         // Find material
126         TObjArray *matList = fluka->GetFlukaMaterials();
127         Int_t nmaterial =  matList->GetEntriesFast();
128         fCMaterial = 0;
129         for (Int_t im = 0; im < nmaterial; im++)
130         {
131             fCMaterial = dynamic_cast<TGeoMaterial*> (matList->At(im));
132             Int_t idmat = fCMaterial->GetIndex();
133             if (idmat == fMedium) break;            
134         } // materials
135         //
136         // Find medium
137         TList *medlist = gGeoManager->GetListOfMedia();
138         TIter next(medlist);
139         while((med = (TGeoMedium*)next()))
140         {
141             mat = med->GetMaterial();
142             if (mat->GetIndex() == fMedium) {
143                 medium = med;
144                 break;
145             }
146         } // media
147         //
148         // Check if sensitive
149         if (medium->GetParam(0) != 0.) mediumIsSensitive = kTRUE;
150
151
152         fprintf(fgFile,"*\n*Material specific process and cut settings for #%8d %s\n", fMedium, fCMaterial->GetName());
153         fCMatMin = fMedium;
154         fCMatMax = fMedium;
155     } else {
156         fprintf(fgFile,"*\n*Global process and cut settings \n");
157         fCMatMin = fgMatMin;
158         fCMatMax = fgMatMax;
159     }
160
161 //
162 // Handle Process Flags 
163 //
164 //
165 //  First make sure that all cuts are taken into account
166     if (DefaultProcessFlag(kPAIR) > 0 && fProcessFlag[kPAIR] == -1 && (fCutValue[kCUTELE] >= 0. || fCutValue[kPPCUTM] >= 0.)) 
167         fProcessFlag[kPAIR] = DefaultProcessFlag(kPAIR);
168     if (DefaultProcessFlag(kBREM) > 0 && fProcessFlag[kBREM] == -1 && (fCutValue[kBCUTE]  >= 0. || fCutValue[kBCUTM] >= 0.)) 
169         fProcessFlag[kBREM] = DefaultProcessFlag(kBREM);
170     if (DefaultProcessFlag(kDRAY) > 0 && fProcessFlag[kDRAY] == -1 && (fCutValue[kDCUTE]  >= 0. || fCutValue[kDCUTM] >= 0.)) 
171         fProcessFlag[kDRAY] = DefaultProcessFlag(kDRAY);
172 //    
173 //
174     if (fProcessFlag[kDCAY] != -1) ProcessDCAY();
175     if (fProcessFlag[kPAIR] != -1) ProcessPAIR();
176     if (fProcessFlag[kBREM] != -1) ProcessBREM();
177     if (fProcessFlag[kCOMP] != -1) ProcessCOMP();
178     if (fProcessFlag[kPHOT] != -1) ProcessPHOT();
179     if (fProcessFlag[kPFIS] != -1) ProcessPFIS();
180     if (fProcessFlag[kANNI] != -1) ProcessANNI();
181     if (fProcessFlag[kMUNU] != -1) ProcessMUNU();
182     if (fProcessFlag[kHADR] != -1) ProcessHADR();
183     if (fProcessFlag[kMULS] != -1) ProcessMULS();
184     if (fProcessFlag[kRAYL] != -1) ProcessRAYL();
185
186     if (fProcessFlag[kLOSS] != -1 || fProcessFlag[kDRAY] != -1)                                    ProcessLOSS();
187     if ((fMedium == -1 && fProcessFlag[kCKOV] > 0) || (fMedium > -1 && fProcessFlag[kCKOV] != -1)) ProcessCKOV();
188     
189 //
190 // Handle Cuts
191 //
192     if (fCutValue[kCUTGAM] >= 0.) ProcessCUTGAM();
193     if (fCutValue[kCUTELE] >= 0.) ProcessCUTELE();
194     if (fCutValue[kCUTNEU] >= 0.) ProcessCUTNEU();
195     if (fCutValue[kCUTHAD] >= 0.) ProcessCUTHAD();
196     if (fCutValue[kCUTMUO] >= 0.) ProcessCUTMUO();
197 //
198 //  Time of flight 
199     if (fCutValue[kTOFMAX] >= 0.) ProcessTOFMAX();
200
201 //
202 //  Tracking precission
203     if (mediumIsSensitive) ProcessSensitiveMedium();
204 }
205
206 void TFlukaConfigOption::ProcessDCAY()
207 {
208     // Process DCAY option
209     fprintf(fgFile,"*\n* --- DCAY --- Decays. Flag = %5d\n", fProcessFlag[kDCAY]);
210     if (fProcessFlag[kDCAY] == 0) {
211         printf("Decays cannot be switched off \n");
212     } else {
213         fprintf(fgFile, "* Decays are on by default\n");
214     }
215 }
216
217
218 void TFlukaConfigOption::ProcessPAIR()
219 {
220     // Process PAIR option
221     fprintf(fgFile,"*\n* --- PAIR --- Pair production by gammas, muons and hadrons. Flag = %5d, PPCUTM = %13.4g, PPCUTE = %13.4g\n", 
222             fProcessFlag[kPAIR], fCutValue[kCUTELE], fCutValue[kPPCUTM]);
223     //
224     // gamma -> e+ e-
225     //
226     if (fProcessFlag[kPAIR] > 0) {
227         fprintf(fgFile,"EMFCUT    %10.1f%10.1f%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",0., 0., 0.0, fCMatMin, fCMatMax, 1.);
228     } else {
229         fprintf(fgFile,"EMFCUT    %10.1f%10.1f%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",0., 0., 1e10,  fCMatMin, fCMatMax, 1.);
230     }
231     
232     //
233     // Direct pair production by Muons and Hadrons
234     //
235     //
236     // Attention ! This card interferes with BREM
237     //
238
239     if (fProcessFlag[kBREM] == -1 ) fProcessFlag[kBREM] = fgDProcessFlag[kBREM];
240     if (fCutValue[kBCUTM]   == -1.) fCutValue[kBCUTM]   = fgDCutValue[kBCUTM];  
241
242
243     Float_t flag = -3.;    
244     if (fProcessFlag[kPAIR] >  0 && fProcessFlag[kBREM] == 0) flag =  1.;
245     if (fProcessFlag[kPAIR] == 0 && fProcessFlag[kBREM]  > 0) flag =  2.;
246     if (fProcessFlag[kPAIR] >  0 && fProcessFlag[kBREM]  > 0) flag =  3.;    
247     if (fProcessFlag[kPAIR] == 0 && fProcessFlag[kBREM] == 0) flag = -3.;
248     // Flag BREM card as handled
249     fProcessFlag[kBREM] = - fProcessFlag[kBREM];
250     
251     //
252     // Energy cut for pair prodution
253     //
254     Float_t cutP = fCutValue[kPPCUTM];
255     if (fCutValue[kPPCUTM]   == -1.) cutP = fgDCutValue[kPPCUTM];       
256     // In G3 this is the cut on the total energy of the e+e- pair
257     // In FLUKA the cut is on the kinetic energy of the electron and poistron
258     cutP = cutP / 2. - 0.51099906e-3;
259     if (cutP < 0.) cutP = 0.;
260     // No explicite generation of e+/e-
261     if (fProcessFlag[kPAIR] == 2) cutP = -1.;
262     //
263     // Energy cut for bremsstrahlung
264     //
265     Float_t cutB = 0.;
266     if (flag > 1.) {
267         fprintf(fgFile,"*\n* +++ BREM --- Bremsstrahlung by muons/hadrons. Flag = %5d, BCUTM = %13.4g \n",
268             fProcessFlag[kBREM], fCutValue[kBCUTM]);
269
270         cutB =  fCutValue[kBCUTM];
271         // No explicite production of gammas
272         if (fProcessFlag[kBREM] == 2) cutB = -1.;
273     }
274
275     fprintf(fgFile,"PAIRBREM  %10.1f%10.4g%10.4g%10.1f%10.1f\n",flag, cutP, cutB, fCMatMin, fCMatMax);
276 }
277
278
279 void TFlukaConfigOption::ProcessBREM()
280 {
281     // Process BREM option
282     fprintf(fgFile,"*\n* --- BREM --- Bremsstrahlung by e+/- and muons/hadrons. Flag = %5d, BCUTE = %13.4g, BCUTM = %13.4g \n",
283             fProcessFlag[kBREM], fCutValue[kBCUTE], fCutValue[kBCUTM]);
284
285     //
286     // e+/- -> e+/- gamma
287     //
288     Float_t cutB = fCutValue[kBCUTE];
289     if (fCutValue[kBCUTE]   == -1.) cutB = fgDCutValue[kBCUTE]; 
290     
291     
292     if (TMath::Abs(fProcessFlag[kBREM]) > 0) {
293         fprintf(fgFile,"EMFCUT    %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fELPO-THR\n",cutB,  0., 0.,  fCMatMin, fCMatMax, 1.);
294     } else {
295         fprintf(fgFile,"EMFCUT    %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fELPO-THR\n",1.e10, 0., 0.,  fCMatMin, fCMatMax, 1.);
296     }
297
298     //
299     // Bremsstrahlung by muons and hadrons
300     //
301     cutB = fCutValue[kBCUTM];
302     if (fCutValue[kBCUTM]   == -1.) cutB = fgDCutValue[kBCUTM]; 
303     if (fProcessFlag[kBREM] == 2) cutB = -1.;
304     Float_t flag = 2.;
305     if (fProcessFlag[kBREM] == 0) flag = -2.;
306     
307     fprintf(fgFile,"PAIRBREM  %10.1f%10.4g%10.4g%10.1f%10.1f\n", flag, 0., cutB, fCMatMin, fCMatMax);
308 }
309
310 void TFlukaConfigOption::ProcessCOMP()
311 {
312     // Process COMP option
313     fprintf(fgFile,"*\n* --- COMP --- Compton scattering  Flag = %5d \n", fProcessFlag[kCOMP]);
314
315     //
316     // Compton scattering
317     //
318
319     if (fProcessFlag[kCOMP] > 0) {
320         fprintf(fgFile,"EMFCUT    %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",0.   , 0., 0.,  fCMatMin, fCMatMax, 1.);
321     } else {
322         fprintf(fgFile,"EMFCUT    %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",1.e10, 0., 0.,  fCMatMin, fCMatMax, 1.);
323     }
324 }
325
326 void TFlukaConfigOption::ProcessPHOT()
327 {
328     // Process PHOS option
329     fprintf(fgFile,"*\n* --- PHOT --- Photoelectric effect. Flag = %5d\n", fProcessFlag[kPHOT]);
330
331     //
332     // Photoelectric effect
333     //
334
335     if (fProcessFlag[kPHOT] > 0) {
336         fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",0., 0., 0.,  fCMatMin, fCMatMax, 1.);
337     } else {
338         fprintf(fgFile,"EMFCUT    %10.1f%10.4g%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",0., 1.e10, 0.,  fCMatMin, fCMatMax, 1.);
339     }
340 }
341
342 void TFlukaConfigOption::ProcessANNI()
343 {
344     // Process ANNI option
345     fprintf(fgFile,"*\n* --- ANNI --- Positron annihilation. Flag = %5d \n", fProcessFlag[kANNI]);
346
347     //
348     // Positron annihilation
349     //
350
351     if (fProcessFlag[kANNI] > 0) {
352         fprintf(fgFile,"EMFCUT    %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fANNH-THR\n",0.   , 0., 0.,  fCMatMin, fCMatMax, 1.);
353     } else {
354         fprintf(fgFile,"EMFCUT    %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fANNH-THR\n",1.e10, 0., 0.,  fCMatMin, fCMatMax, 1.);
355     }
356 }
357
358
359 void TFlukaConfigOption::ProcessPFIS()
360 {
361     // Process PFIS option
362     fprintf(fgFile,"*\n* --- PFIS --- Photonuclear interaction  Flag = %5d\n", fProcessFlag[kPFIS]);
363
364     //
365     // Photonuclear interactions
366     //
367
368     if (fProcessFlag[kPFIS] > 0) {
369         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.);
370     } else {
371         fprintf(fgFile,"PHOTONUC  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",-1.                          , 0., 0.,  fCMatMin, fCMatMax, 1.);
372     }
373 }
374
375 void TFlukaConfigOption::ProcessMUNU()
376 {
377     // Process MUNU option
378     fprintf(fgFile,"*\n* --- MUNU --- Muon nuclear interaction. Flag = %5d\n", fProcessFlag[kMUNU]);
379     
380     //
381     // Muon nuclear interactions
382     //
383     if (fProcessFlag[kMUNU] > 0) {
384         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.);
385     } else {
386         fprintf(fgFile,"MUPHOTON  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",-1.                          , 0.,   0.,    fCMatMin, fCMatMax, 1.);
387     }
388 }
389
390 void TFlukaConfigOption::ProcessRAYL()
391 {
392     // Process RAYL option
393     fprintf(fgFile,"*\n* --- RAYL --- Rayleigh Scattering. Flag = %5d\n", fProcessFlag[kRAYL]);
394
395     //
396     // Rayleigh scattering
397     //
398     Int_t nreg;
399     Int_t* reglist = fgGeom->GetMaterialList(fMedium, nreg);
400     //
401     // Loop over regions of a given material
402     for (Int_t k = 0; k < nreg; k++) {
403         Float_t ireg = reglist[k];
404         if (fProcessFlag[kRAYL] > 0) {
405             fprintf(fgFile,"EMFRAY    %10.1f%10.1f%10.1f%10.1f\n", 1., ireg, ireg, 1.);
406         } else {
407             fprintf(fgFile,"EMFRAY    %10.1f%10.1f%10.1f%10.1f\n", 3., ireg, ireg, 1.);
408         }
409     }
410 }
411
412 void TFlukaConfigOption::ProcessCKOV()
413 {
414     // Process CKOV option
415     fprintf(fgFile,"*\n* --- CKOV --- Cerenkov Photon production.  %5d\n", fProcessFlag[kCKOV]);
416
417     //
418     // Cerenkov photon production
419     //
420
421     TFluka* fluka = (TFluka*) gMC;
422     TObjArray *matList = fluka->GetFlukaMaterials();
423     Int_t nmaterial =  matList->GetEntriesFast();
424     for (Int_t im = 0; im < nmaterial; im++)
425     {
426         TGeoMaterial* material = dynamic_cast<TGeoMaterial*> (matList->At(im));
427         Int_t idmat = material->GetIndex();
428 //
429 // Check if global option
430         if (fMedium != -1 && idmat != fMedium) continue;
431         
432         TFlukaCerenkov* cerenkovProp;
433         if (!(cerenkovProp = dynamic_cast<TFlukaCerenkov*>(material->GetCerenkovProperties()))) continue;
434         //
435         // This medium has Cerenkov properties 
436         //
437         //
438         if (fMedium == -1 || (fMedium != -1 && fProcessFlag[kCKOV] > 0)) {
439             // Write OPT-PROD card for each medium 
440             Float_t  emin  = cerenkovProp->GetMinimumEnergy();
441             Float_t  emax  = cerenkovProp->GetMaximumEnergy();        
442             fprintf(fgFile, "OPT-PROD  %10.4g%10.4g%10.4g%10.4g%10.4g%10.4gCERENKOV\n", emin, emax, 0., 
443                     Float_t(idmat), Float_t(idmat), 0.); 
444             //
445             // Write OPT-PROP card for each medium 
446             // Forcing FLUKA to call user routines (queffc.cxx, rflctv.cxx, rfrndx.cxx)
447             //
448             fprintf(fgFile, "OPT-PROP  %10.4g%10.4g%10.4g%10.1f%10.1f%10.1fWV-LIMIT\n",  
449                     cerenkovProp->GetMinimumWavelength(), cerenkovProp->GetMaximumWavelength(), cerenkovProp->GetMaximumWavelength(), 
450                     Float_t(idmat), Float_t(idmat), 0.0);
451             
452
453             fprintf(fgFile, "OPT-PROP  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n", -100., -100., -100., 
454                     Float_t(idmat), Float_t(idmat), 0.0);
455             
456             for (Int_t j = 0; j < 3; j++) {
457                 fprintf(fgFile, "OPT-PROP  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f&\n", -100., -100., -100., 
458                         Float_t(idmat), Float_t(idmat), 0.0);
459             }
460
461
462             // Photon detection efficiency user defined     
463             if (cerenkovProp->IsSensitive())
464                 fprintf(fgFile, "OPT-PROP  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fSENSITIV\n", -100., -100., -100., 
465                         Float_t(idmat), Float_t(idmat), 0.0);
466             // Material has a reflective surface
467             if (cerenkovProp->IsMetal())
468                 fprintf(fgFile, "OPT-PROP  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fMETAL\n", -100., -100., -100., 
469                         Float_t(idmat), Float_t(idmat), 0.0);
470
471         } else {
472             fprintf(fgFile,"OPT-PROD  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fCERE-OFF\n",0., 0., 0., fCMatMin, fCMatMax, 1.);
473         }
474     }
475 }
476
477
478 void TFlukaConfigOption::ProcessHADR()
479 {
480     // Process HADR option
481     fprintf(fgFile,"*\n* --- HADR --- Hadronic interactions. Flag = %5d\n", fProcessFlag[kHADR]);
482     
483     if (fProcessFlag[kHADR] > 0) {
484         fprintf(fgFile,"*\n*Hadronic interaction is ON by default in FLUKA\n");
485     } else {
486         if (fMedium != -1) printf("Hadronic interactions cannot be switched off material by material !\n");
487         fprintf(fgFile,"THRESHOL  %10.1f%10.1f%10.1f%10.1e%10.1f\n",0., 0., 0., 1.e10, 0.);
488     }
489 }
490
491
492
493 void TFlukaConfigOption::ProcessMULS()
494 {
495     // Process MULS option
496     fprintf(fgFile,"*\n* --- MULS --- Muliple Scattering. Flag = %5d\n", fProcessFlag[kMULS]);
497     //
498     // Multiple scattering
499     //
500     if (fProcessFlag[kMULS] > 0) {
501         fprintf(fgFile,"*\n*Multiple scattering is  ON by default in FLUKA\n");
502     } else {
503         fprintf(fgFile,"MULSOPT   %10.1f%10.1f%10.1f%10.1f%10.1f\n",0., 3., 3., fCMatMin, fCMatMax);
504     }
505 }
506
507 void TFlukaConfigOption::ProcessLOSS()
508 {
509     // Process LOSS option
510     fprintf(fgFile,"*\n* --- LOSS --- Ionisation energy loss. Flags: LOSS = %5d, DRAY = %5d, STRA = %5d; Cuts: DCUTE = %13.4g, DCUTM = %13.4g \n",
511             fProcessFlag[kLOSS], fProcessFlag[kDRAY], fProcessFlag[kSTRA], fCutValue[kDCUTE], fCutValue[kDCUTM]);
512     //
513     // Ionisation energy loss
514     //
515     //
516     // Impose consistency
517     
518     if (fProcessFlag[kLOSS] == 1 || fProcessFlag[kLOSS] == 3) {
519         fProcessFlag[kDRAY] = 1;
520     } else if (fProcessFlag[kLOSS] == 2) {
521         fProcessFlag[kDRAY] = 0;
522         fCutValue[kDCUTE] = 1.e10;
523         fCutValue[kDCUTM] = 1.e10;      
524     } else {
525         if (fProcessFlag[kDRAY] == 1) {
526             fProcessFlag[kLOSS] = 1;
527         } else if (fProcessFlag[kDRAY] == 0) {
528             fProcessFlag[kLOSS] = 2;
529             fCutValue[kDCUTE] = 1.e10;
530             fCutValue[kDCUTM] = 1.e10;  
531         }
532     }
533     
534     if (fCutValue[kDCUTE] == -1.) fCutValue[kDCUTE] = fgDCutValue[kDCUTE];
535     if (fCutValue[kDCUTM] == -1.) fCutValue[kDCUTM] = fgDCutValue[kDCUTM];    
536     
537     Float_t cutM =  fCutValue[kDCUTM];    
538         
539
540     if (fProcessFlag[kSTRA] == -1) fProcessFlag[kSTRA] = fgDProcessFlag[kSTRA];
541     if (fProcessFlag[kSTRA] ==  0) fProcessFlag[kSTRA] = 1;
542     Float_t stra = (Float_t) fProcessFlag[kSTRA];
543     
544
545     if (fProcessFlag[kLOSS] == 1 || fProcessFlag[kLOSS] == 3) {
546 //
547 // Restricted energy loss fluctuations 
548 //
549         fprintf(fgFile,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f\n", 1., 1., stra, fCMatMin, fCMatMax);
550         fprintf(fgFile,"DELTARAY  %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n", cutM, 0., 0., fCMatMin, fCMatMax, 1.);
551     } else if (fProcessFlag[kLOSS] == 4) {
552 //
553 // No fluctuations
554 //
555         fprintf(fgFile,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f\n",-1., -1., stra, fCMatMin, fCMatMax);        
556         fprintf(fgFile,"DELTARAY  %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n", 1.e10, 0., 0., fCMatMin, fCMatMax, 1.);      
557     }  else {
558 //
559 // Full fluctuations
560 //
561         fprintf(fgFile,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f\n",1., 1., stra, fCMatMin, fCMatMax);  
562         fprintf(fgFile,"DELTARAY  %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n", 1.e10, 0., 0., fCMatMin, fCMatMax, 1.);      
563     }
564 }
565
566
567 void TFlukaConfigOption::ProcessCUTGAM()
568 {
569 // Cut on gammas
570 //
571     fprintf(fgFile,"*\n*Cut for Gammas. CUTGAM = %13.4g\n", fCutValue[kCUTGAM]);
572     if (fMedium == -1) {
573         fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", 
574                 0., fCutValue[kCUTGAM], 0., 0., Float_t(fgGeom->NofVolumes()), 1.);
575
576     } else {
577         Int_t nreg, *reglist;
578         Float_t ireg;
579         reglist = fgGeom->GetMaterialList(fMedium, nreg);
580         // Loop over regions of a given material
581         for (Int_t k = 0; k < nreg; k++) {
582             ireg = reglist[k];
583             fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", 0.,fCutValue[kCUTGAM], 0., ireg, ireg, 1.);
584         }
585     }
586     
587     // Transport production cut used for pemf
588     //
589     //  FUDGEM paramter. The parameter takes into account th contribution of atomic electrons to multiple scattering.
590     //  For production and transport cut-offs larger than 100 keV it must be set = 1.0, while in the keV region it must be
591     Float_t parFudgem = (fCutValue[kCUTGAM] > 1.e-4)? 1.0 : 0.0 ;
592     fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1fPROD-CUT\n", 
593             0., fCutValue[kCUTGAM], parFudgem, fCMatMin, fCMatMax, 1.);
594 }
595
596 void TFlukaConfigOption::ProcessCUTELE()
597 {
598 // Cut on e+/e-
599 //
600     fprintf(fgFile,"*\n*Cut for e+/e-. CUTELE = %13.4g\n", fCutValue[kCUTELE]);
601     if (fMedium == -1) {
602         fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", 
603                 -fCutValue[kCUTELE], 0., 0., 0., Float_t(fgGeom->NofVolumes()), 1.);
604     } else {
605         Int_t nreg, *reglist;
606         Float_t ireg;
607         reglist = fgGeom->GetMaterialList(fMedium, nreg);
608         // Loop over regions of a given material
609         for (Int_t k = 0; k < nreg; k++) {
610             ireg = reglist[k];
611             fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", -fCutValue[kCUTELE], 0., 0., ireg, ireg, 1.);
612         }
613     }
614     // Transport production cut used for pemf
615     //
616     //  FUDGEM paramter. The parameter takes into account th contribution of atomic electrons to multiple scattering.
617     //  For production and transport cut-offs larger than 100 keV it must be set = 1.0, while in the keV region it must be
618     Float_t parFudgem = (fCutValue[kCUTELE] > 1.e-4)? 1.0 : 0.0;
619     fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1fPROD-CUT\n", 
620             -fCutValue[kCUTELE], 0., parFudgem, fCMatMin, fCMatMax, 1.);
621 }
622
623 void TFlukaConfigOption::ProcessCUTNEU()
624 {
625   // Cut on neutral hadrons
626   fprintf(fgFile,"*\n*Cut for neutral hadrons. CUTNEU = %13.4g\n", fCutValue[kCUTNEU]);
627   
628   // Cut on neutral hadrons
629   fprintf(fgFile,"*\n*Cut for neutral hadrons. CUTNEU = %13.4g\n", fCutValue[kCUTNEU]);
630   
631   // Energy group structure of the 72-neutron ENEA data set:
632   const Float_t neutronGroupUpLimit[72] = {
633     1.9600E-02, 1.7500E-02, 1.4918E-02, 1.3499E-02, 1.2214E-02, 1.1052E-02, 1.0000E-02, 9.0484E-03,
634     8.1873E-03, 7.4082E-03, 6.7032E-03, 6.0653E-03, 5.4881E-03, 4.9659E-03, 4.4933E-03, 4.0657E-03,
635     3.6788E-03, 3.3287E-03, 3.0119E-03, 2.7253E-03, 2.4660E-03, 2.2313E-03, 2.0190E-03, 1.8268E-03,
636     1.6530E-03, 1.4957E-03, 1.3534E-03, 1.2246E-03, 1.1080E-03, 1.0026E-03, 9.0718E-04, 8.2085E-04,
637     7.4274E-04, 6.0810E-04, 4.9787E-04, 4.0762E-04, 3.3373E-04, 2.7324E-04, 2.2371E-04, 1.8316E-04,
638     1.4996E-04, 1.2277E-04, 8.6517E-05, 5.2475E-05, 3.1828E-05, 2.1852E-05, 1.5034E-05, 1.0332E-05,
639     7.1018E-06, 4.8809E-06, 3.3546E-06, 2.3054E-06, 1.5846E-06, 1.0446E-06, 6.8871E-07, 4.5400E-07,
640     2.7537E-07, 1.6702E-07, 1.0130E-07, 6.1442E-08, 3.7267E-08, 2.2603E-08, 1.5535E-08, 1.0677E-08,
641     7.3375E-09, 5.0435E-09, 3.4662E-09, 2.3824E-09, 1.6374E-09, 1.1254E-09, 6.8257E-10, 4.1400E-10
642   };
643
644   Float_t cut = fCutValue[kCUTNEU];
645   //
646   // 8.0 = Neutron
647   // 9.0 = Antineutron
648   // Find the FLUKA neutron group corresponding to the cut
649   //
650   Float_t neutronCut = cut;
651   Int_t groupCut = 1; // if cut is > 19.6 MeV not low energy neutron transport is done
652   if (neutronCut < 0.0196) {
653     neutronCut = 0.0196;
654     // Search the group cutoff for the energy cut
655     Int_t i;
656     for( i=0; i<72; i++ ) {
657       if (cut > neutronGroupUpLimit[i]) break;
658     }
659     groupCut = i+1;
660   } 
661   
662   if (fMedium == -1) {
663         Float_t cut = fCutValue[kCUTNEU];
664         // 8.0 = Neutron
665         // 9.0 = Antineutron
666         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -neutronCut,  8.0,  9.0);
667         fprintf(fgFile,"LOW-BIAS  %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n",
668                 Float_t(groupCut), 73.0, 0.95, 2., Float_t(fgGeom->NofVolumes()), 1.);
669         //
670         //
671         // 12.0 = Kaon zero long
672         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 12.0, 12.0);
673         // 17.0 = Lambda, 18.0 = Antilambda
674         // 19.0 = Kaon zero short
675         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 17.0, 19.0);
676         // 22.0 = Sigma zero, Pion zero, Kaon zero
677         // 25.0 = Antikaon zero
678         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 22.0, 25.0);
679         // 32.0 = Antisigma zero
680         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 32.0, 32.0);
681         // 34.0 = Xi zero
682         // 35.0 = AntiXi zero
683         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 34.0, 35.0);
684         // 47.0 = D zero
685         // 48.0 = AntiD zero
686         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 47.0, 48.0);
687         // 53.0 = Xi_c zero
688         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 53.0, 53.0);
689         // 55.0 = Xi'_c zero
690         // 56.0 = Omega_c zero
691         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 55.0, 56.0);
692         // 59.0 = AntiXi_c zero
693         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 59.0, 59.0);
694         // 61.0 = AntiXi'_c zero
695         // 62.0 = AntiOmega_c zero
696         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 61.0, 62.0);
697     } else {
698         Int_t nreg, *reglist;
699         Float_t ireg;
700         reglist = fgGeom->GetMaterialList(fMedium, nreg);
701         // Loop over regions of a given material
702         for (Int_t k = 0; k < nreg; k++) {
703           ireg = reglist[k];
704           fprintf(fgFile,"LOW-BIAS  %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n",
705                   Float_t(groupCut), 73.0, 0.95, ireg, ireg, 1.);
706         }
707
708         Warning("ProcessCUTNEU", 
709                 "Material #%4d %s: Cut on neutral hadrons (Ekin > %9.3e) material by material only implemented for low-energy neutrons !\n", 
710                 fMedium, fCMaterial->GetName(), cut);
711     }
712 }
713
714 void TFlukaConfigOption::ProcessCUTHAD()
715 {
716     // Cut on charged hadrons
717     fprintf(fgFile,"*\n*Cut for charge hadrons. CUTHAD = %13.4g\n", fCutValue[kCUTHAD]);
718     Float_t cut = fCutValue[kCUTHAD];
719     if (fMedium == -1) {
720         // 1.0 = Proton
721         // 2.0 = Antiproton
722         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut,  1.0,  2.0);
723         // 13.0 = Positive Pion, Negative Pion, Positive Kaon
724         // 16.0 = Negative Kaon
725         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 13.0, 16.0);
726         // 20.0 = Negative Sigma
727         // 21.0 = Positive Sigma
728         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 20.0, 21.0);
729         // 31.0 = Antisigma minus
730         // 33.0 = Antisigma plus
731         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 31.0, 31.0);
732         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 33.0, 33.0);
733         // 36.0 = Negative Xi, Positive Xi, Omega minus
734         // 39.0 = Antiomega
735         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 36.0, 39.0);
736         // 45.0 = D plus
737         // 46.0 = D minus
738         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 45.0, 46.0);
739         // 49.0 = D_s plus, D_s minus, Lambda_c plus
740         // 52.0 = Xi_c plus
741         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 49.0, 52.0);
742         // 54.0 = Xi'_c plus
743         // 60.0 = AntiXi'_c minus
744         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 54.0, 54.0);
745         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 60.0, 60.0);
746         // 57.0 = Antilambda_c minus
747         // 58.0 = AntiXi_c minus
748         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 57.0, 58.0);
749     } else {
750       Warning("ProcessCUTHAD", 
751               "Material #%4d %s: Cut on charged hadrons (Ekin > 9.3e) material by material not yet implemented !\n", 
752               fMedium, fCMaterial->GetName(), cut); 
753     }
754 }
755
756 void TFlukaConfigOption::ProcessCUTMUO()
757 {
758     // Cut on muons
759     fprintf(fgFile,"*\n*Cut for muons. CUTMUO = %13.4g\n", fCutValue[kCUTMUO]);
760     Float_t cut = fCutValue[kCUTMUO];
761     if (fMedium == -1) {
762         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n",-cut, 10.0, 11.0);
763     } else {
764         Warning("ProcessCUTMUO", "Material #%4d %s: Cut on muons (Ekin > %9.3e) material by material not yet implemented !\n", 
765                 fMedium, fCMaterial->GetName(), cut);
766     }
767     
768     
769 }
770
771 void TFlukaConfigOption::ProcessTOFMAX()
772 {
773     // Cut time of flight
774     Float_t cut = fCutValue[kTOFMAX];
775     fprintf(fgFile,"*\n*Cut on time of flight. TOFMAX = %13.4g\n", fCutValue[kTOFMAX]);
776     fprintf(fgFile,"TIME-CUT  %10.4g%10.1f%10.1f%10.1f%10.1f\n",cut*1.e9,0.,0.,-6.0,64.0);
777 }
778
779 void  TFlukaConfigOption::ProcessSensitiveMedium()
780 {
781     //
782     // Special options for sensitive media
783     //
784
785     fprintf(fgFile,"*\n*Options for sensitive medium \n");
786     //
787     // EMFFIX
788     fprintf(fgFile,"EMFFIX    %10.1f%10.3f%10.1f%10.1f%10.1f%10.1f\n", fCMatMin, 0.05, 0., 0., 0., 0.);
789     //
790     // FLUKAFIX
791     fprintf(fgFile,"FLUKAFIX  %10.3f                    %10.3f\n", 0.05, fCMatMin);
792 }