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