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 }