TFluka/TFlukaConfigOption.cxx

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