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 }