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