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