From fb2cbbecf3f0d2e56028a85236a6d54e670dfd8e Mon Sep 17 00:00:00 2001 From: morsch Date: Thu, 16 Dec 2004 09:42:39 +0000 Subject: [PATCH] Writing of FLUKA input cards for physics configuration is now the responsibility of TFlukaConfigOptions. One TFlukaConfigOption per material with special process flags and cuts is created. One object with material id = -1 is created to hold the global settings. --- TFluka/TFluka.cxx | 1258 +++------------------------------ TFluka/TFluka.h | 6 +- TFluka/TFlukaConfigOption.cxx | 618 +++++++++++++++- TFluka/TFlukaConfigOption.h | 73 +- 4 files changed, 764 insertions(+), 1191 deletions(-) diff --git a/TFluka/TFluka.cxx b/TFluka/TFluka.cxx index 0f48827c017..261289ca38d 100644 --- a/TFluka/TFluka.cxx +++ b/TFluka/TFluka.cxx @@ -103,8 +103,7 @@ TFluka::TFluka() :TVirtualMC(), fVerbosityLevel(0), fInputFileName(""), - fProcesses(0), - fCuts(0), + fUserConfig(0), fUserScore(0) { // @@ -132,8 +131,7 @@ TFluka::TFluka(const char *title, Int_t verbosity, Bool_t isRootGeometrySupporte fTrackIsEntering(0), fTrackIsExiting(0), fTrackIsNew(0), - fProcesses(new TObjArray(100)), - fCuts(new TObjArray(100)), + fUserConfig(new TObjArray(100)), fUserScore(new TObjArray(100)) { // create geometry interface @@ -148,7 +146,7 @@ TFluka::TFluka(const char *title, Int_t verbosity, Bool_t isRootGeometrySupporte fFieldFlag = 1; fGeneratePemf = kFALSE; fMCGeo = new TGeoMCGeometry("MCGeo", "TGeo Implementation of VirtualMCGeometry", kTRUE); - fGeom = new TFlukaMCGeometry("geom", "FLUKA VMC Geometry"); + fGeom = new TFlukaMCGeometry("geom", "FLUKA VMC Geometry"); if (verbosity > 2) fGeom->SetDebugMode(kTRUE); fMaterials = 0; fStopped = 0; @@ -166,14 +164,9 @@ TFluka::~TFluka() { delete fGeom; delete fMCGeo; - if (fCuts) { - fCuts->Delete(); - delete fCuts; - } - - if (fProcesses) { - fProcesses->Delete(); - delete fProcesses; + if (fUserConfig) { + fUserConfig->Delete(); + delete fUserConfig; } @@ -378,6 +371,7 @@ void TFluka::Material(Int_t& kmat, const char* name, Double_t a, Double_t z, Double_t dens, Double_t radl, Double_t absl, Double_t* /*buf*/, Int_t /*nwbuf*/) { // +// Define a material TGeoMaterial *mat; kmat = gGeoManager->GetListOfMaterials()->GetSize(); if ((z-Int_t(z)) > 1E-3) { @@ -394,6 +388,8 @@ void TFluka::Material(Int_t& kmat, const char* name, Double_t a, //______________________________________________________________________________ void TFluka::Mixture(Int_t& kmat, const char *name, Float_t *a, Float_t *z, Double_t dens, Int_t nlmat, Float_t *wmat) { +// +// Define a material mixture // Double_t* da = fGeom->CreateDoubleArray(a, TMath::Abs(nlmat)); Double_t* dz = fGeom->CreateDoubleArray(z, TMath::Abs(nlmat)); @@ -819,37 +815,31 @@ void TFluka::StopTrack() void TFluka::SetProcess(const char* flagName, Int_t flagValue, Int_t imed) { // Set process user flag for material imat -// - TFlukaConfigOption* proc = new TFlukaConfigOption(flagName, flagValue, imed); - fProcesses->Add(proc); -} - -//______________________________________________________________________________ -Bool_t TFluka::SetProcess(const char* flagName, Int_t flagValue) -{ -// Set process user flag // // // Update if already in the list // - - TIter next(fProcesses); + TIter next(fUserConfig); TFlukaConfigOption* proc; while((proc = (TFlukaConfigOption*)next())) { - if (strcmp(proc->GetName(), flagName) == 0) { - proc->SetFlag(flagValue); - proc->SetMedium(-1); - return kTRUE; - } + if (proc->Medium() == imed) { + proc->SetProcess(flagName, flagValue); + return; + } } + proc = new TFlukaConfigOption(imed); + proc->SetProcess(flagName, flagValue); + fUserConfig->Add(proc); +} + +//______________________________________________________________________________ +Bool_t TFluka::SetProcess(const char* flagName, Int_t flagValue) +{ +// Set process user flag // -// If not create a new process // - - proc = new TFlukaConfigOption(flagName, flagValue); - fProcesses->Add(proc); - + SetProcess(flagName, flagValue, -1); return kTRUE; } @@ -858,8 +848,19 @@ void TFluka::SetCut(const char* cutName, Double_t cutValue, Int_t imed) { // Set user cut value for material imed // - TFlukaConfigOption* cut = new TFlukaConfigOption(cutName, cutValue, imed); - fCuts->Add(cut); + TIter next(fUserConfig); + TFlukaConfigOption* proc; + while((proc = (TFlukaConfigOption*)next())) + { + if (proc->Medium() == imed) { + proc->SetCut(cutName, cutValue); + return; + } + } + + proc = new TFlukaConfigOption(imed); + proc->SetCut(cutName, cutValue); + fUserConfig->Add(proc); } //______________________________________________________________________________ @@ -868,26 +869,8 @@ Bool_t TFluka::SetCut(const char* cutName, Double_t cutValue) // Set user cut value // // -// Update if already in the list -// - - TIter next(fCuts); - TFlukaConfigOption* cut; - while((cut = (TFlukaConfigOption*)next())) - { - if (strcmp(cut->GetName(), cutName) == 0) { - cut->SetCut(cutValue); - return kTRUE; - } - } -// -// If not create a new process -// - - cut = new TFlukaConfigOption(cutName, cutValue); - fCuts->Add(cut); - - return kTRUE; + SetCut(cutName, cutValue, -1); + return kTRUE; } void TFluka::SetUserScoring(const char* option, Int_t npar, Float_t what[12]) @@ -913,1117 +896,92 @@ void TFluka::InitPhysics() // // Physics initialisation with preparation of FLUKA input cards // - printf("=>InitPhysics\n"); - Int_t j, k; - Double_t theCut; - - FILE *pFlukaVmcCoreInp, *pFlukaVmcFlukaMat, *pFlukaVmcInp; - - Double_t zero = 0.0; - Double_t one = 1.0; - Double_t two = 2.0; - Double_t three = 3.0; - - Float_t fLastMaterial = fGeom->GetLastMaterialIndex(); - if (fVerbosityLevel >= 3) printf(" last FLUKA material is %g\n", fLastMaterial); - - // Prepare Cerenkov - TObjArray *matList = GetFlukaMaterials(); - Int_t nmaterial = matList->GetEntriesFast(); - fMaterials = new Int_t[nmaterial+3]; - -// construct file names + printf("=>InitPhysics\n"); - TString sFlukaVmcCoreInp = getenv("ALICE_ROOT"); - sFlukaVmcCoreInp +="/TFluka/input/"; - TString sFlukaVmcTmp = "flukaMat.inp"; - TString sFlukaVmcInp = GetInputFileName(); - sFlukaVmcCoreInp += GetCoreInputFileName(); - -// open files - - if ((pFlukaVmcCoreInp = fopen(sFlukaVmcCoreInp.Data(),"r")) == NULL) { - printf("\nCannot open file %s\n",sFlukaVmcCoreInp.Data()); - exit(1); - } - if ((pFlukaVmcFlukaMat = fopen(sFlukaVmcTmp.Data(),"r")) == NULL) { - printf("\nCannot open file %s\n",sFlukaVmcTmp.Data()); - exit(1); - } - if ((pFlukaVmcInp = fopen(sFlukaVmcInp.Data(),"w")) == NULL) { - printf("\nCannot open file %s\n",sFlukaVmcInp.Data()); - exit(1); - } +// Construct file names + FILE *pFlukaVmcCoreInp, *pFlukaVmcFlukaMat, *pFlukaVmcInp; + TString sFlukaVmcCoreInp = getenv("ALICE_ROOT"); + sFlukaVmcCoreInp +="/TFluka/input/"; + TString sFlukaVmcTmp = "flukaMat.inp"; + TString sFlukaVmcInp = GetInputFileName(); + sFlukaVmcCoreInp += GetCoreInputFileName(); + +// Open files + if ((pFlukaVmcCoreInp = fopen(sFlukaVmcCoreInp.Data(),"r")) == NULL) { + printf("\nCannot open file %s\n",sFlukaVmcCoreInp.Data()); + exit(1); + } + if ((pFlukaVmcFlukaMat = fopen(sFlukaVmcTmp.Data(),"r")) == NULL) { + printf("\nCannot open file %s\n",sFlukaVmcTmp.Data()); + exit(1); + } + if ((pFlukaVmcInp = fopen(sFlukaVmcInp.Data(),"w")) == NULL) { + printf("\nCannot open file %s\n",sFlukaVmcInp.Data()); + exit(1); + } -// copy core input file - Char_t sLine[255]; - Float_t fEventsPerRun; - - while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) { - if (strncmp(sLine,"GEOEND",6) != 0) - fprintf(pFlukaVmcInp,"%s",sLine); // copy until GEOEND card - else { - fprintf(pFlukaVmcInp,"GEOEND\n"); // add GEOEND card - goto flukamat; - } - } // end of while until GEOEND card - +// Copy core input file + Char_t sLine[255]; + Float_t fEventsPerRun; + + while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) { + if (strncmp(sLine,"GEOEND",6) != 0) + fprintf(pFlukaVmcInp,"%s",sLine); // copy until GEOEND card + else { + fprintf(pFlukaVmcInp,"GEOEND\n"); // add GEOEND card + goto flukamat; + } + } // end of while until GEOEND card + flukamat: - while ((fgets(sLine,255,pFlukaVmcFlukaMat)) != NULL) { // copy flukaMat.inp file - fprintf(pFlukaVmcInp,"%s\n",sLine); - } - - while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) { - if (strncmp(sLine,"START",5) != 0) - fprintf(pFlukaVmcInp,"%s\n",sLine); - else { - sscanf(sLine+10,"%10f",&fEventsPerRun); - goto fin; - } - } //end of while until START card - -fin: -// in G3 the process control values meaning can be different for -// different processes, but for most of them is: -// 0 process is not activated -// 1 process is activated WITH generation of secondaries -// 2 process is activated WITHOUT generation of secondaries -// if process does not generate secondaries => 1 same as 2 -// -// Exceptions: -// MULS: also 3 -// LOSS: also 3, 4 -// RAYL: only 0,1 -// HADR: may be > 2 -// - -// Loop over number of SetProcess calls - fprintf(pFlukaVmcInp,"*----------------------------------------------------------------------------- \n"); - fprintf(pFlukaVmcInp,"*----- The following data are generated from SetProcess and SetCut calls ----- \n"); - fprintf(pFlukaVmcInp,"*----------------------------------------------------------------------------- \n"); - -// Outer loop over processes - TIter next(fProcesses); - TFlukaConfigOption *proc; -// Inner loop over processes - TIter nextp(fProcesses); - TFlukaConfigOption *procp; -// Loop over cuts - TIter nextc(fCuts); - TFlukaConfigOption *cut = 0x0; - - while((proc = (TFlukaConfigOption*)next())) { - Float_t matMin = three; - Float_t matMax = fLastMaterial; - Bool_t global = kTRUE; - if (proc->Medium() != -1) { - Int_t mat; - if ((mat = proc->Medium()) >= GetFlukaMaterials()->GetEntries()) continue; - matMin = Float_t(mat); - matMax = matMin; - global = kFALSE; - - fprintf(pFlukaVmcInp,"*\n*Material specific process setting for #%8d \n", mat); - } - - // annihilation - // G3 default value: 1 - // G4 processes: G4eplusAnnihilation/G4IeplusAnnihilation - // Particles: e+ - // Physics: EM - // flag = 0 no annihilation - // flag = 1 annihilation, decays processed - // flag = 2 annihilation, no decay product stored - // gMC ->SetProcess("ANNI",1); // EMFCUT -1. 0. 0. 3. lastmat 0. ANNH-THR - if (strncmp(proc->GetName(),"ANNI",4) == 0) { - if (proc->Flag() == 1 || proc->Flag() == 2) { - fprintf(pFlukaVmcInp,"*\n*Kinetic energy threshold (GeV) for e+ annihilation - resets to default=0.\n"); - fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('ANNI',1) or SetProcess('ANNI',2)\n"); - // -one = kinetic energy threshold (GeV) for e+ annihilation (resets to default=0) - // zero = not used - // zero = not used - // matMin = lower bound of the material indices in which the respective thresholds apply - // matMax = upper bound of the material indices in which the respective thresholds apply - // one = step length in assigning indices - // "ANNH-THR"; - fprintf(pFlukaVmcInp,"EMFCUT %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fANNH-THR\n",-one,zero,zero,matMin,matMax,one); - } - else if (proc->Flag() == 0) { - fprintf(pFlukaVmcInp,"*\n*No annihilation - no FLUKA card generated\n"); - fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('ANNI',0)\n"); - } - else { - fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('ANNI',?) call.\n"); - fprintf(pFlukaVmcInp,"*No FLUKA card generated\n"); - } - } + while ((fgets(sLine,255,pFlukaVmcFlukaMat)) != NULL) { // copy flukaMat.inp file + fprintf(pFlukaVmcInp,"%s\n",sLine); + } - // bremsstrahlung and pair production are both activated - // G3 default value: 1 - // G4 processes: G4eBremsstrahlung/G4IeBremsstrahlung, - // G4MuBremsstrahlung/G4IMuBremsstrahlung, - // G4LowEnergyBremstrahlung - // Particles: e-/e+; mu+/mu- - // Physics: EM - // flag = 0 no bremsstrahlung - // flag = 1 bremsstrahlung, photon processed - // flag = 2 bremsstrahlung, no photon stored - // gMC ->SetProcess("BREM",1); // PAIRBREM 2. 0. 0. 3. lastmat - // EMFCUT -1. 0. 0. 3. lastmat 0. ELPO-THR - // G3 default value: 1 - // G4 processes: G4GammaConversion, - // G4MuPairProduction/G4IMuPairProduction - // G4LowEnergyGammaConversion - // Particles: gamma, mu - // Physics: EM - // flag = 0 no delta rays - // flag = 1 delta rays, secondaries processed - // flag = 2 delta rays, no secondaries stored - - else if ((strncmp(proc->GetName(),"PAIR",4) == 0) && (proc->Flag() == 1 || proc->Flag() == 2)) { - - nextp.Reset(); - - while ((procp = (TFlukaConfigOption*)nextp())) { - if ((strncmp(procp->GetName(),"BREM",4) == 0) && - (proc->Flag() == 1) && - (procp->Medium() == proc->Medium())) { - fprintf(pFlukaVmcInp,"*\n*Bremsstrahlung and pair production by muons and charged hadrons both activated\n"); - fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('BREM',1) and SetProcess('PAIR',1)\n"); - fprintf(pFlukaVmcInp,"*Energy threshold set by call SetCut('BCUTM',cut) or set to 0.\n"); - fprintf(pFlukaVmcInp,"*Energy threshold set by call SetCut('PPCUTM',cut) or set to 0.\n"); - // three = bremsstrahlung and pair production by muons and charged hadrons both are activated - fprintf(pFlukaVmcInp,"PAIRBREM %10.1f",three); - // direct pair production by muons - // G4 particles: "e-", "e+" - // G3 default value: 0.01 GeV - //gMC ->SetCut("PPCUTM",cut); // total energy cut for direct pair prod. by muons - theCut = 0.0; - nextc.Reset(); - while ((cut = (TFlukaConfigOption*)nextc())) { - if (strncmp(cut->GetName(), "PPCUTM", 6) == 0 && - (cut->Medium() == proc->Medium())) theCut = cut->Cut(); - } - fprintf(pFlukaVmcInp,"%10.4g",theCut); - // theCut; = e+, e- kinetic energy threshold (in GeV) for explicit pair production. - // muon and hadron bremsstrahlung - // G4 particles: "gamma" - // G3 default value: CUTGAM=0.001 GeV - //gMC ->SetCut("BCUTM",cut); // cut for muon and hadron bremsstrahlung - theCut = 0.0; - nextc.Reset(); - while ((cut = (TFlukaConfigOption*)nextc())) { - if (strncmp(cut->GetName(), "BCUTM", 5) == 0 && - (cut->Medium() == proc->Medium())) theCut = cut->Cut(); - } - fprintf(pFlukaVmcInp,"%10.4g%10.1f%10.1f\n",theCut,matMin,matMax); - // theCut = photon energy threshold (GeV) for explicit bremsstrahlung production - // matMin = lower bound of the material indices in which the respective thresholds apply - // matMax = upper bound of the material indices in which the respective thresholds apply - - // for e+ and e- - fprintf(pFlukaVmcInp,"*\n*Kinetic energy threshold (GeV) for e+/e- bremsstrahlung - resets to default=0.\n"); - fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('BREM',1);\n"); - theCut = -1.0; - nextc.Reset(); - while ((cut = (TFlukaConfigOption*)nextc())) { - if (strncmp(cut->GetName(), "BCUTE", 5) == 0 && - (cut->Medium() == proc->Medium())) theCut = cut->Cut(); - } - //theCut = kinetic energy threshold (GeV) for e+/e- bremsstrahlung (resets to default=0) - // zero = not used - // zero = not used - // matMin = lower bound of the material indices in which the respective thresholds apply - // matMax = upper bound of the material indices in which the respective thresholds apply - // one = step length in assigning indices - // "ELPO-THR"; - fprintf(pFlukaVmcInp,"EMFCUT %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fELPO-THR\n",theCut,zero,zero,matMin,matMax,one); - - // for gamma -> e+ and e- - fprintf(pFlukaVmcInp,"*\n*Pair production by photons is activated\n"); - fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PAIR',1);\n"); - theCut = -1.0; - nextc.Reset(); - while ((cut = (TFlukaConfigOption*)nextc())) { - if (strncmp(cut->GetName(), "CUTELE", 6) == 0 && - (cut->Medium() == proc->Medium())) theCut = cut->Cut(); - } - // theCut = energy threshold (GeV) for gamma pair production (< 0.0 : resets to default, = 0.0 : ignored) - // matMin = lower bound of the material indices in which the respective thresholds apply - // matMax = upper bound of the material indices in which the respective thresholds apply - // one = step length in assigning indices - fprintf(pFlukaVmcInp,"EMFCUT %10.1f%10.1f%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",zero,zero,theCut,matMin,matMax,one); - goto BOTH; - } // end of if for BREM - } // end of loop for BREM - - // only pair production by muons and charged hadrons is activated - fprintf(pFlukaVmcInp,"*\n*Pair production by muons and charged hadrons is activated\n"); - fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PAIR',1) or SetProcess('PAIR',2)\n"); - fprintf(pFlukaVmcInp,"*Energy threshold set by call SetCut('PPCUTM',cut) or set to 0.\n"); - // direct pair production by muons - // G4 particles: "e-", "e+" - // G3 default value: 0.01 GeV - //gMC ->SetCut("PPCUTM",cut); // total energy cut for direct pair prod. by muons - // one = pair production by muons and charged hadrons is activated - // zero = e+, e- kinetic energy threshold (in GeV) for explicit pair production. - // zero = no explicit bremsstrahlung production is simulated - // matMin = lower bound of the material indices in which the respective thresholds apply - // matMax = upper bound of the material indices in which the respective thresholds apply - fprintf(pFlukaVmcInp,"PAIRBREM %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,zero,zero,matMin,matMax); - - // for gamma -> e+ and e- - fprintf(pFlukaVmcInp,"*\n*Pair production by electrons is activated\n"); - fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PAIR',1) or SetProcess('PAIR',2)\n"); - theCut = -1.0; - nextc.Reset(); - while ((cut = (TFlukaConfigOption*)nextc())) { - if (strncmp(cut->GetName(), "CUTELE", 6) == 0 && - (cut->Medium() == proc->Medium())) theCut = cut->Cut(); + while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) { + if (strncmp(sLine,"START",5) != 0) + fprintf(pFlukaVmcInp,"%s\n",sLine); + else { + sscanf(sLine+10,"%10f",&fEventsPerRun); + goto fin; } - // zero = energy threshold (GeV) for Compton scattering (= 0.0 : ignored) - // zero = energy threshold (GeV) for Photoelectric (= 0.0 : ignored) - // theCut = energy threshold (GeV) for gamma pair production (< 0.0 : resets to default, = 0.0 : ignored) - // matMin = lower bound of the material indices in which the respective thresholds apply - // matMax = upper bound of the material indices in which the respective thresholds apply - // one = step length in assigning indices - fprintf(pFlukaVmcInp,"EMFCUT %10.1f%10.1f%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",zero,zero,theCut,matMin,matMax,one); - - BOTH: - k = 0; - } // end of if for PAIR - - - - // bremsstrahlung - // G3 default value: 1 - // G4 processes: G4eBremsstrahlung/G4IeBremsstrahlung, - // G4MuBremsstrahlung/G4IMuBremsstrahlung, - // G4LowEnergyBremstrahlung - // Particles: e-/e+; mu+/mu- - // Physics: EM - // flag = 0 no bremsstrahlung - // flag = 1 bremsstrahlung, photon processed - // flag = 2 bremsstrahlung, no photon stored - // gMC ->SetProcess("BREM",1); // PAIRBREM 2. 0. 0. 3. lastmat - // EMFCUT -1. 0. 0. 3. lastmat 0. ELPO-THR - else if (strncmp(proc->GetName(),"BREM",4) == 0) { - nextp.Reset(); - while((procp = (TFlukaConfigOption*)nextp())) { - if ((strncmp(procp->GetName(),"PAIR",4) == 0) && - procp->Flag() == 1 && - (procp->Medium() == proc->Medium())) goto NOBREM; - } - if (proc->Flag() == 1) { - fprintf(pFlukaVmcInp,"*\n*Bremsstrahlung by muons and charged hadrons is activated\n"); - fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('BREM',1) or SetProcess('BREM',2)\n"); - fprintf(pFlukaVmcInp,"*Energy threshold set by call SetCut('BCUTM',cut) or set to 0.\n"); - // two = bremsstrahlung by muons and charged hadrons is activated - // zero = no meaning - // muon and hadron bremsstrahlung - // G4 particles: "gamma" - // G3 default value: CUTGAM=0.001 GeV - //gMC ->SetCut("BCUTM",cut); // cut for muon and hadron bremsstrahlung - theCut = 0.0; - nextc.Reset(); - while ((cut = (TFlukaConfigOption*)nextc())) { - if (strncmp(cut->GetName(), "BCUTM", 5) == 0 && - (cut->Medium() == proc->Medium())) theCut = cut->Cut(); - } - // theCut = photon energy threshold (GeV) for explicit bremsstrahlung production - // matMin = lower bound of the material indices in which the respective thresholds apply - // matMax = upper bound of the material indices in which the respective thresholds apply - fprintf(pFlukaVmcInp,"PAIRBREM %10.1f%10.1f%10.4g%10.1f%10.1f\n",two,zero,theCut,matMin,matMax); - - // for e+ and e- - fprintf(pFlukaVmcInp,"*\n*Kinetic energy threshold (GeV) for e+/e- bremsstrahlung - resets to default=0.\n"); - fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('BREM',1);"); - // - one = kinetic energy threshold (GeV) for e+/e- bremsstrahlung (resets to default=0) - // zero = not used - // zero = not used - // matMin = lower bound of the material indices in which the respective thresholds apply - // matMax = upper bound of the material indices in which the respective thresholds apply - // one = step length in assigning indices - //"ELPO-THR"; - theCut = -1.0; - nextc.Reset(); - while ((cut = (TFlukaConfigOption*)nextc())) { - if (strncmp(cut->GetName(), "CUTGAM", 6) == 0 && - (cut->Medium() == proc->Medium())) theCut = cut->Cut(); - } - fprintf(pFlukaVmcInp,"EMFCUT %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fELPO-THR\n", theCut,zero,zero,matMin,matMax,one); - } - else if (proc->Flag() == 0) { - fprintf(pFlukaVmcInp,"*\n*No bremsstrahlung - no FLUKA card generated\n"); - fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('BREM',0)\n"); - } - else { - fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('BREM',?) call.\n"); - fprintf(pFlukaVmcInp,"*No FLUKA card generated\n"); - } - NOBREM: - j = 0; - } // end of else if (strncmp(proc->GetName(),"BREM",4) == 0) - - // Cerenkov photon generation - // G3 default value: 0 - // G4 process: G4Cerenkov - // - // Particles: charged - // Physics: Optical - // flag = 0 no Cerenkov photon generation - // flag = 1 Cerenkov photon generation - // flag = 2 Cerenkov photon generation with primary stopped at each step - //xx gMC ->SetProcess("CKOV",1); // ??? Cerenkov photon generation - - else if (strncmp(proc->GetName(),"CKOV",4) == 0) { - if ((proc->Flag() == 1 || proc->Flag() == 2) && global) { - // Write comments - fprintf(pFlukaVmcInp, "* \n"); - fprintf(pFlukaVmcInp, "*Cerenkov photon generation\n"); - fprintf(pFlukaVmcInp, "*Generated from call: SetProcess('CKOV',1) or SetProcess('CKOV',2)\n"); - // Loop over media - for (Int_t im = 0; im < nmaterial; im++) - { - TGeoMaterial* material = dynamic_cast (matList->At(im)); - Int_t idmat = material->GetIndex(); - - if (!global && idmat != proc->Medium()) continue; - - fMaterials[idmat] = im; - // Skip media with no Cerenkov properties - TFlukaCerenkov* cerenkovProp; - if (!(cerenkovProp = dynamic_cast(material->GetCerenkovProperties()))) continue; - // - // This medium has Cerenkov properties - // - // - // Write OPT-PROD card for each medium - Float_t emin = cerenkovProp->GetMinimumEnergy(); - Float_t emax = cerenkovProp->GetMaximumEnergy(); - fprintf(pFlukaVmcInp, "OPT-PROD %10.4g%10.4g%10.4g%10.4g%10.4g%10.4gCERENKOV\n", emin, emax, 0., - Float_t(idmat), Float_t(idmat), 0.); - // - // Write OPT-PROP card for each medium - // Forcing FLUKA to call user routines (queffc.cxx, rflctv.cxx, rfrndx.cxx) - // - fprintf(pFlukaVmcInp, "OPT-PROP %10.4g%10.4g%10.4g%10.1f%10.1f%10.1fWV-LIMIT\n", - cerenkovProp->GetMinimumWavelength(), - cerenkovProp->GetMaximumWavelength(), - cerenkovProp->GetMaximumWavelength(), - Float_t(idmat), Float_t(idmat), 0.0); - - if (cerenkovProp->IsMetal()) { - fprintf(pFlukaVmcInp, "OPT-PROP %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fMETAL\n", - -100., -100., -100., - Float_t(idmat), Float_t(idmat), 0.0); - } else { - fprintf(pFlukaVmcInp, "OPT-PROP %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n", - -100., -100., -100., - Float_t(idmat), Float_t(idmat), 0.0); - } - - - for (Int_t j = 0; j < 3; j++) { - fprintf(pFlukaVmcInp, "OPT-PROP %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f&\n", - -100., -100., -100., - Float_t(idmat), Float_t(idmat), 0.0); - } - // Photon detection efficiency user defined - - if (cerenkovProp->IsSensitive()) - fprintf(pFlukaVmcInp, "OPT-PROP %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fSENSITIV\n", - -100., -100., -100., - Float_t(idmat), Float_t(idmat), 0.0); - - } // materials - } else if (proc->Flag() == 0) { - fprintf(pFlukaVmcInp,"*\n*No Cerenkov photon generation\n"); - fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('CKOV',0)\n"); - // zero = not used - // zero = not used - // zero = not used - // matMin = lower bound of the material indices in which the respective thresholds apply - // matMax = upper bound of the material indices in which the respective thresholds apply - // one = step length in assigning indices - //"CERE-OFF"; - fprintf(pFlukaVmcInp,"OPT-PROD %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fCERE-OFF\n",zero,zero,zero,matMin,matMax,one); - } - else { - fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('CKOV',?) call.\n"); - fprintf(pFlukaVmcInp,"*No FLUKA card generated\n"); - } - } // end of else if (strncmp(proc->GetName(),"CKOV",4) == 0) - - // Compton scattering - // G3 default value: 1 - // G4 processes: G4ComptonScattering, - // G4LowEnergyCompton, - // G4PolarizedComptonScattering - // Particles: gamma - // Physics: EM - // flag = 0 no Compton scattering - // flag = 1 Compton scattering, electron processed - // flag = 2 Compton scattering, no electron stored - // gMC ->SetProcess("COMP",1); // EMFCUT -1. 0. 0. 3. lastmat 0. PHOT-THR - else if (strncmp(proc->GetName(),"COMP",4) == 0) { - if (proc->Flag() == 1 || proc->Flag() == 2) { - fprintf(pFlukaVmcInp,"*\n*Energy threshold (GeV) for Compton scattering - resets to default=0.\n"); - fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('COMP',1);\n"); - // - one = energy threshold (GeV) for Compton scattering - resets to default=0. - // zero = not used - // zero = not used - // matMin = lower bound of the material indices in which the respective thresholds apply - // matMax = upper bound of the material indices in which the respective thresholds apply - // one = step length in assigning indices - //"PHOT-THR"; - theCut = -1.0; - nextc.Reset(); - while ((cut = (TFlukaConfigOption*)nextc())) { - if (strncmp(cut->GetName(), "CUTELE", 6) == 0 && - (cut->Medium() == proc->Medium())) theCut = cut->Cut(); - } - fprintf(pFlukaVmcInp,"EMFCUT %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",theCut,zero,zero,matMin,matMax,one); - } - else if (proc->Flag() == 0) { - fprintf(pFlukaVmcInp,"*\n*No Compton scattering - no FLUKA card generated\n"); - fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('COMP',0)\n"); - } - else { - fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('COMP',?) call.\n"); - fprintf(pFlukaVmcInp,"*No FLUKA card generated\n"); - } - } // end of else if (strncmp(proc->GetName(),"COMP",4) == 0) - - // decay - // G3 default value: 1 - // G4 process: G4Decay - // - // Particles: all which decay is applicable for - // Physics: General - // flag = 0 no decays - // flag = 1 decays, secondaries processed - // flag = 2 decays, no secondaries stored - //gMC ->SetProcess("DCAY",0); // not available - else if ((strncmp(proc->GetName(),"DCAY",4) == 0) && proc->Flag() == 0) - cout << "SetProcess for flag =" << proc->GetName() << " value=" << proc->Flag() << " not avaliable!" << endl; - else if ((strncmp(proc->GetName(),"DCAY",4) == 0) && proc->Flag() == 1) { - // Nothing to do decays are switched on by default - } - - - // delta-ray - // G3 default value: 2 - // !! G4 treats delta rays in different way - // G4 processes: G4eIonisation/G4IeIonization, - // G4MuIonisation/G4IMuIonization, - // G4hIonisation/G4IhIonisation - // Particles: charged - // Physics: EM - // flag = 0 no energy loss - // flag = 1 restricted energy loss fluctuations - // flag = 2 complete energy loss fluctuations - // flag = 3 same as 1 - // flag = 4 no energy loss fluctuations - // gMC ->SetProcess("DRAY",0); // DELTARAY 1.E+6 0. 0. 3. lastmat 0. - else if (strncmp(proc->GetName(),"DRAY",4) == 0) { - if (proc->Flag() == 0 || proc->Flag() == 4) { - fprintf(pFlukaVmcInp,"*\n*Kinetic energy threshold (GeV) for delta ray production\n"); - fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('DRAY',0) or SetProcess('DRAY',4)\n"); - fprintf(pFlukaVmcInp,"*No delta ray production by muons - threshold set artificially high\n"); - Double_t emin = 1.0e+6; // kinetic energy threshold (GeV) for delta ray production (discrete energy transfer) - // zero = ignored - // zero = ignored - // matMin = lower bound of the material indices in which the respective thresholds apply - // matMax = upper bound of the material indices in which the respective thresholds apply - // one = step length in assigning indices - fprintf(pFlukaVmcInp,"DELTARAY %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n",emin,zero,zero,matMin,matMax,one); - } - else if (proc->Flag() == 1 || proc->Flag() == 2 || proc->Flag() == 3) { - fprintf(pFlukaVmcInp,"*\n*Kinetic energy threshold (GeV) for delta ray production\n"); - fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('DRAY',flag), flag=1,2,3\n"); - fprintf(pFlukaVmcInp,"*Delta ray production by muons switched on\n"); - fprintf(pFlukaVmcInp,"*Energy threshold set by call SetCut('DCUTM',cut) or set to 1.0e+6.\n"); - theCut = 1.0e+6; - nextc.Reset(); - // - // Check cut one delta-rays from electrons - // - while ((cut = (TFlukaConfigOption*)nextc())) { - if (strncmp(cut->GetName(), "DCUTM", 5) == 0 && - cut->Medium() == proc->Medium()) theCut = cut->Cut(); - } - // theCut = kinetic energy threshold (GeV) for delta ray production (discrete energy transfer) - // zero = ignored - // zero = ignored - // matMin = lower bound of the material indices in which the respective thresholds apply - // matMax = upper bound of the material indices in which the respective thresholds apply - // one = step length in assigning indices - fprintf(pFlukaVmcInp,"DELTARAY %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n",theCut,zero,zero,matMin,matMax,one); - } - else { - fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('DRAY',?) call.\n"); - fprintf(pFlukaVmcInp,"*No FLUKA card generated\n"); - } - } // end of else if (strncmp(proc->GetName(),"DRAY",4) == 0) - - // hadronic process - // G3 default value: 1 - // G4 processes: all defined by TG4PhysicsConstructorHadron - // - // Particles: hadrons - // Physics: Hadron - // flag = 0 no multiple scattering - // flag = 1 hadronic interactions, secondaries processed - // flag = 2 hadronic interactions, no secondaries stored - // gMC ->SetProcess("HADR",1); // ??? hadronic process - //Select pure GEANH (HADR 1) or GEANH/NUCRIN (HADR 3) ????? - else if (strncmp(proc->GetName(),"HADR",4) == 0) { - if (proc->Flag() == 1 || proc->Flag() == 2) { - fprintf(pFlukaVmcInp,"*\n*Hadronic interaction is ON by default in FLUKA\n"); - fprintf(pFlukaVmcInp,"*No FLUKA card generated\n"); - } - else if (proc->Flag() == 0) { - fprintf(pFlukaVmcInp,"*\n*Hadronic interaction is set OFF\n"); - fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('HADR',0);\n"); - fprintf(pFlukaVmcInp,"*Switching off hadronic interactions not foreseen in FLUKA\n"); - fprintf(pFlukaVmcInp,"THRESHOL %10.1f%10.1f%10.1f%10.1e%10.1f\n",zero, zero, zero, 1.e10, zero); - } - else { - fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('HADR',?) call.\n"); - fprintf(pFlukaVmcInp,"*No FLUKA card generated\n"); - } - } // end of else if (strncmp(proc->GetName(),"HADR",4) == 0) - - - // energy loss - // G3 default value: 2 - // G4 processes: G4eIonisation/G4IeIonization, - // G4MuIonisation/G4IMuIonization, - // G4hIonisation/G4IhIonisation - // - // Particles: charged - // Physics: EM - // flag=0 no energy loss - // flag=1 restricted energy loss fluctuations - // flag=2 complete energy loss fluctuations - // flag=3 same as 1 - // flag=4 no energy loss fluctuations - // If the value ILOSS is changed, then (in G3) cross-sections and energy - // loss tables must be recomputed via the command 'PHYSI' - // gMC ->SetProcess("LOSS",2); // ??? IONFLUCT ? energy loss - else if (strncmp(proc->GetName(),"LOSS",4) == 0) { - if (proc->Flag() > 0 || proc->Flag() < 4) { // restricted energy loss fluctuations - fprintf(pFlukaVmcInp,"*\n*Restricted energy loss fluctuations\n"); - fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('LOSS',1) or SetProcess('LOSS',3)\n"); - // one = restricted energy loss fluctuations (for hadrons and muons) switched on - // one = restricted energy loss fluctuations (for e+ and e-) switched on - // one = minimal accuracy - // matMin = lower bound of the material indices in which the respective thresholds apply - // upper bound of the material indices in which the respective thresholds apply - fprintf(pFlukaVmcInp,"IONFLUCT %10.1f%10.1f%10.1f%10.1f%10.1f\n",one, one, 4., matMin, matMax); - } - else if (proc->Flag() == 4) { // no energy loss fluctuations - fprintf(pFlukaVmcInp,"*\n*No energy loss fluctuations\n"); - fprintf(pFlukaVmcInp,"*\n*Generated from call: SetProcess('LOSS',4)\n"); - // - one = restricted energy loss fluctuations (for hadrons and muons) switched off - // - one = restricted energy loss fluctuations (for e+ and e-) switched off - // one = minimal accuracy - // matMin = lower bound of the material indices in which the respective thresholds apply - // matMax = upper bound of the material indices in which the respective thresholds apply - fprintf(pFlukaVmcInp,"IONFLUCT %10.1f%10.1f%10.1f%10.1f%10.1f\n",-one,-one,one,matMin,matMax); - } - else { - fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('LOSS',?) call.\n"); - fprintf(pFlukaVmcInp,"*No FLUKA card generated\n"); - } - } // end of else if (strncmp(proc->GetName(),"LOSS",4) == 0) - - - // multiple scattering - // G3 default value: 1 - // G4 process: G4MultipleScattering/G4IMultipleScattering - // - // Particles: charged - // Physics: EM - // flag = 0 no multiple scattering - // flag = 1 Moliere or Coulomb scattering - // flag = 2 Moliere or Coulomb scattering - // flag = 3 Gaussian scattering - // gMC ->SetProcess("MULS",1); // MULSOPT multiple scattering - else if (strncmp(proc->GetName(),"MULS",4) == 0) { - if (proc->Flag() == 1 || proc->Flag() == 2 || proc->Flag() == 3) { - fprintf(pFlukaVmcInp,"*\n*Multiple scattering is ON by default for e+e- and for hadrons/muons\n"); - fprintf(pFlukaVmcInp,"*No FLUKA card generated\n"); - } - else if (proc->Flag() == 0) { - fprintf(pFlukaVmcInp,"*\n*Multiple scattering is set OFF\n"); - fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('MULS',0);\n"); - // zero = ignored - // three = multiple scattering for hadrons and muons is completely suppressed - // three = multiple scattering for e+ and e- is completely suppressed - // matMin = lower bound of the material indices in which the respective thresholds apply - // matMax = upper bound of the material indices in which the respective thresholds apply - fprintf(pFlukaVmcInp,"MULSOPT %10.1f%10.1f%10.1f%10.1f%10.1f\n",zero,three,three,matMin,matMax); - } - else { - fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('MULS',?) call.\n"); - fprintf(pFlukaVmcInp,"*No FLUKA card generated\n"); - } - } // end of else if (strncmp(proc->GetName(),"MULS",4) == 0) - - - // muon nuclear interaction - // G3 default value: 0 - // G4 processes: G4MuNuclearInteraction, - // G4MuonMinusCaptureAtRest - // - // Particles: mu - // Physics: Not set - // flag = 0 no muon-nuclear interaction - // flag = 1 nuclear interaction, secondaries processed - // flag = 2 nuclear interaction, secondaries not processed - // gMC ->SetProcess("MUNU",1); // MUPHOTON 1. 0. 0. 3. lastmat - else if (strncmp(proc->GetName(),"MUNU",4) == 0) { - if (proc->Flag() == 1) { - fprintf(pFlukaVmcInp,"*\n*Muon nuclear interactions with production of secondary hadrons\n"); - fprintf(pFlukaVmcInp,"*\n*Generated from call: SetProcess('MUNU',1);\n"); - // one = full simulation of muon nuclear interactions and production of secondary hadrons - // zero = ratio of longitudinal to transverse virtual photon cross-section - Default = 0.25. - // zero = fraction of rho-like interactions ( must be < 1) - Default = 0.75. - // matMin = lower bound of the material indices in which the respective thresholds apply - // matMax = upper bound of the material indices in which the respective thresholds apply - fprintf(pFlukaVmcInp,"MUPHOTON %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,zero,zero,matMin,matMax); - } - else if (proc->Flag() == 2) { - fprintf(pFlukaVmcInp,"*\n*Muon nuclear interactions without production of secondary hadrons\n"); - fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('MUNU',2);\n"); - // two = full simulation of muon nuclear interactions and production of secondary hadrons - // zero = ratio of longitudinal to transverse virtual photon cross-section - Default = 0.25. - // zero = fraction of rho-like interactions ( must be < 1) - Default = 0.75. - // matMin = lower bound of the material indices in which the respective thresholds apply - // matMax = upper bound of the material indices in which the respective thresholds apply - fprintf(pFlukaVmcInp,"MUPHOTON %10.1f%10.1f%10.1f%10.1f%10.1f\n",two,zero,zero,matMin,matMax); - } - else if (proc->Flag() == 0) { - fprintf(pFlukaVmcInp,"*\n*No muon nuclear interaction - no FLUKA card generated\n"); - fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('MUNU',0)\n"); - } - else { - fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('MUNU',?) call.\n"); - fprintf(pFlukaVmcInp,"*No FLUKA card generated\n"); - } - } // end of else if (strncmp(proc->GetName(),"MUNU",4) == 0) - - - // photofission - // G3 default value: 0 - // G4 process: ?? - // - // Particles: gamma - // Physics: ?? - // gMC ->SetProcess("PFIS",0); // PHOTONUC -1. 0. 0. 3. lastmat 0. - // flag = 0 no photon fission - // flag = 1 photon fission, secondaries processed - // flag = 2 photon fission, no secondaries stored - else if (strncmp(proc->GetName(),"PFIS",4) == 0) { - if (proc->Flag() == 0) { - fprintf(pFlukaVmcInp,"*\n*No photonuclear interactions\n"); - fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PFIS',0);\n"); - // - one = no photonuclear interactions - // zero = not used - // zero = not used - // matMin = lower bound of the material indices in which the respective thresholds apply - // matMax = upper bound of the material indices in which the respective thresholds apply - fprintf(pFlukaVmcInp,"PHOTONUC %10.1f%10.1f%10.1f%10.1f%10.1f\n",-one,zero,zero,matMin,matMax); - } - else if (proc->Flag() == 1) { - fprintf(pFlukaVmcInp,"*\n*Photon nuclear interactions are activated at all energies\n"); - fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PFIS',1);\n"); - // one = photonuclear interactions are activated at all energies - // zero = not used - // zero = not used - // matMin = lower bound of the material indices in which the respective thresholds apply - // matMax = upper bound of the material indices in which the respective thresholds apply - fprintf(pFlukaVmcInp,"PHOTONUC %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,zero,zero,matMin,matMax); - } - else if (proc->Flag() == 0) { - fprintf(pFlukaVmcInp,"*\n*No photofission - no FLUKA card generated\n"); - fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PFIS',0)\n"); - } - else { - fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('PFIS',?) call.\n"); - fprintf(pFlukaVmcInp,"*No FLUKA card generated\n"); - } - } - - - // photo electric effect - // G3 default value: 1 - // G4 processes: G4PhotoElectricEffect - // G4LowEnergyPhotoElectric - // Particles: gamma - // Physics: EM - // flag = 0 no photo electric effect - // flag = 1 photo electric effect, electron processed - // flag = 2 photo electric effect, no electron stored - // gMC ->SetProcess("PHOT",1); // EMFCUT 0. -1. 0. 3. lastmat 0. PHOT-THR - else if (strncmp(proc->GetName(),"PHOT",4) == 0) { - if (proc->Flag() == 1) { - fprintf(pFlukaVmcInp,"*\n*Photo electric effect is activated\n"); - fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PHOT',1);\n"); - // zero = ignored - // - one = resets to default=0. - // zero = ignored - // matMin = lower bound of the material indices in which the respective thresholds apply - // matMax = upper bound of the material indices in which the respective thresholds apply - // one = step length in assigning indices - //"PHOT-THR"; - theCut = -1.0; - nextc.Reset(); - while ((cut = (TFlukaConfigOption*)nextc())) { - if (strncmp(cut->GetName(), "CUTELE", 6) == 0 && - (cut->Medium() == proc->Medium())) theCut = cut->Cut(); - } - fprintf(pFlukaVmcInp,"EMFCUT %10.1f%10.4g%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",zero,theCut,zero,matMin,matMax,one); - } - else if (proc->Flag() == 0) { - fprintf(pFlukaVmcInp,"*\n*No photo electric effect - no FLUKA card generated\n"); - fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PHOT',0)\n"); - } - else { - fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('PHOT',?) call.\n"); - fprintf(pFlukaVmcInp,"*No FLUKA card generated\n"); - } - } // else if (strncmp(proc->GetName(),"PHOT",4) == 0) - - - // Rayleigh scattering - // G3 default value: 0 - // G4 process: G4OpRayleigh - // - // Particles: optical photon - // Physics: Optical - // flag = 0 Rayleigh scattering off - // flag = 1 Rayleigh scattering on - //xx gMC ->SetProcess("RAYL",1); - else if (strncmp(proc->GetName(),"RAYL",4) == 0) { - if (proc->Flag() == 1) { - fprintf(pFlukaVmcInp,"*\n*Rayleigh scattering is ON by default in FLUKA\n"); - fprintf(pFlukaVmcInp,"*No FLUKA card generated\n"); - } - else if (proc->Flag() == 0) { - fprintf(pFlukaVmcInp,"*\n*Rayleigh scattering is set OFF\n"); - fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('RAYL',0);\n"); - // - one = no Rayleigh scattering and no binding corrections for Compton - // matMin = lower bound of the material indices in which the respective thresholds apply - // matMax = upper bound of the material indices in which the respective thresholds apply - fprintf(pFlukaVmcInp,"EMFRAY %10.1f%10.1f%10.1f%10.1f\n",-one,three,matMin,matMax); - } - else { - fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('RAYL',?) call.\n"); - fprintf(pFlukaVmcInp,"*No FLUKA card generated\n"); - } - } // end of else if (strncmp(proc->GetName(),"RAYL",4) == 0) - - - // synchrotron radiation in magnetic field - // G3 default value: 0 - // G4 process: G4SynchrotronRadiation - // - // Particles: ?? - // Physics: Not set - // flag = 0 no synchrotron radiation - // flag = 1 synchrotron radiation - //xx gMC ->SetProcess("SYNC",1); // synchrotron radiation generation - else if (strncmp(proc->GetName(),"SYNC",4) == 0) { - fprintf(pFlukaVmcInp,"*\n*Synchrotron radiation generation is NOT implemented in FLUKA\n"); - fprintf(pFlukaVmcInp,"*No FLUKA card generated\n"); - } - - - // Automatic calculation of tracking medium parameters - // flag = 0 no automatic calculation - // flag = 1 automatic calculation - //xx gMC ->SetProcess("AUTO",1); // ??? automatic computation of the tracking medium parameters - else if (strncmp(proc->GetName(),"AUTO",4) == 0) { - fprintf(pFlukaVmcInp,"*\n*Automatic calculation of tracking medium parameters is always ON in FLUKA\n"); - fprintf(pFlukaVmcInp,"*No FLUKA card generated\n"); - } - - - // To control energy loss fluctuation model - // flag = 0 Urban model - // flag = 1 PAI model - // flag = 2 PAI+ASHO model (not active at the moment) - //xx gMC ->SetProcess("STRA",1); // ??? energy fluctuation model - else if (strncmp(proc->GetName(),"STRA",4) == 0) { - if (proc->Flag() == 0 || proc->Flag() == 2 || proc->Flag() == 3) { - fprintf(pFlukaVmcInp,"*\n*Ionization energy losses calculation is activated\n"); - fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('STRA',n);, n=0,1,2\n"); - // one = restricted energy loss fluctuations (for hadrons and muons) switched on - // one = restricted energy loss fluctuations (for e+ and e-) switched on - // one = minimal accuracy - // matMin = lower bound of the material indices in which the respective thresholds apply - // matMax = upper bound of the material indices in which the respective thresholds apply - fprintf(pFlukaVmcInp,"IONFLUCT %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,one,one,matMin,matMax); - } - else { - fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('STRA',?) call.\n"); - fprintf(pFlukaVmcInp,"*No FLUKA card generated\n"); - } - } // else if (strncmp(proc->GetName(),"STRA",4) == 0) - - - - - else { // processes not yet treated - - // light photon absorption (Cerenkov photons) - // it is turned on when Cerenkov process is turned on - // G3 default value: 0 - // G4 process: G4OpAbsorption, G4OpBoundaryProcess - // - // Particles: optical photon - // Physics: Optical - // flag = 0 no absorption of Cerenkov photons - // flag = 1 absorption of Cerenkov photons - // gMC ->SetProcess("LABS",2); // ??? Cerenkov light absorption - - - - cout << "SetProcess for flag=" << proc->GetName() << " value=" << proc->Flag() << " not yet implemented!" << endl; - } - } //end of loop number of SetProcess calls - - -// Loop over number of SetCut calls - - nextc.Reset(); - while ((cut = (TFlukaConfigOption*)nextc())) { - Float_t matMin = three; - Float_t matMax = fLastMaterial; - Bool_t global = kTRUE; - if (cut->Medium() != -1) { - Int_t mat; - if ((mat = cut->Medium()) >= GetFlukaMaterials()->GetEntries()) continue; - matMin = Float_t(mat); - matMax = matMin; - global = kFALSE; - TGeoMaterial* material = (TGeoMaterial*) (GetFlukaMaterials())->At(GetMaterialIndex(mat)); - fprintf(pFlukaVmcInp,"*\n*Material specific cut setting for #%8d %s %s %13.3e\n", - mat, material->GetName(), cut->GetName(), cut->Cut()); - - } - - // cuts handled in SetProcess calls - if (strncmp(cut->GetName(),"BCUTM",5) == 0) continue; - else if (strncmp(cut->GetName(),"BCUTE",5) == 0) continue; - else if (strncmp(cut->GetName(),"DCUTM",5) == 0) continue; - else if (strncmp(cut->GetName(),"PPCUTM",6) == 0) continue; - - // gammas - // G4 particles: "gamma" - // G3 default value: 0.001 GeV - // gMC ->SetCut("CUTGAM",cut); // cut for gammas - - else if (strncmp(cut->GetName(),"CUTGAM",6) == 0 && global) { - fprintf(pFlukaVmcInp,"*\n*Cut for gamma\n"); - fprintf(pFlukaVmcInp,"*Generated from call: SetCut('CUTGAM',cut);\n"); - fprintf(pFlukaVmcInp,"EMFCUT %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", - zero, cut->Cut(), zero, zero, Float_t(fGeom->NofVolumes()), one); - } - else if (strncmp(cut->GetName(),"CUTGAM",6) == 0 && !global) { - // loop over materials for EMFCUT FLUKA cards - for (j=0; j < matMax-matMin+1; j++) { - Int_t nreg, imat, *reglist; - Float_t ireg; - imat = (Int_t) matMin + j; - reglist = fGeom->GetMaterialList(imat, nreg); - // loop over regions of a given material - for (Int_t k = 0; k < nreg; k++) { - ireg = reglist[k]; - fprintf(pFlukaVmcInp,"EMFCUT %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", zero, cut->Cut(), zero, ireg, ireg, one); - } - } - } // end of else if for gamma - - - // electrons - // G4 particles: "e-" - // ?? positrons - // G3 default value: 0.001 GeV - //gMC ->SetCut("CUTELE",cut); // cut for e+,e- - else if (strncmp(cut->GetName(),"CUTELE",6) == 0 && global) { - fprintf(pFlukaVmcInp,"*\n*Cut for electrons\n"); - fprintf(pFlukaVmcInp,"*Generated from call: SetCut('CUTELE',cut);\n"); - fprintf(pFlukaVmcInp,"EMFCUT %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", - -cut->Cut(), zero, zero, zero, Float_t(fGeom->NofVolumes()), one); - } - else if (strncmp(cut->GetName(),"CUTELE",6) == 0 && !global) { - // loop over materials for EMFCUT FLUKA cards - for (j=0; j < matMax-matMin+1; j++) { - Int_t nreg, imat, *reglist; - Float_t ireg; - imat = (Int_t) matMin + j; - reglist = fGeom->GetMaterialList(imat, nreg); - // loop over regions of a given material - for (k=0; kCut(), zero, zero, ireg, ireg, one); - } - } - } // end of else if for electrons + } //end of while until START card + + fin: + // Pass information to configuration objects - // neutral hadrons - // G4 particles: of type "baryon", "meson", "nucleus" with zero charge - // G3 default value: 0.01 GeV - //gMC ->SetCut("CUTNEU",cut); // cut for neutral hadrons - else if (strncmp(cut->GetName(),"CUTNEU",6) == 0 && global) { - fprintf(pFlukaVmcInp,"*\n*Cut for neutral hadrons\n"); - fprintf(pFlukaVmcInp,"*Generated from call: SetCut('CUTNEU',cut);\n"); - - // 8.0 = Neutron - // 9.0 = Antineutron - fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),8.0,9.0); - - // 12.0 = Kaon zero long - // 12.0 = Kaon zero long - fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),12.0,12.0); - - // 17.0 = Lambda, 18.0 = Antilambda - // 19.0 = Kaon zero short - fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),17.0,19.0); - - // 22.0 = Sigma zero, Pion zero, Kaon zero - // 25.0 = Antikaon zero - fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),22.0,25.0); - - // 32.0 = Antisigma zero - // 32.0 = Antisigma zero - fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),32.0,32.0); - - // 34.0 = Xi zero - // 35.0 = AntiXi zero - fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),34.0,35.0); - - // 47.0 = D zero - // 48.0 = AntiD zero - fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),47.0,48.0); - - // 53.0 = Xi_c zero - // 53.0 = Xi_c zero - fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),53.0,53.0); - - // 55.0 = Xi'_c zero - // 56.0 = Omega_c zero - fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),55.0,56.0); - - // 59.0 = AntiXi_c zero - // 59.0 = AntiXi_c zero - fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),59.0,59.0); - - // 61.0 = AntiXi'_c zero - // 62.0 = AntiOmega_c zero - fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),61.0,62.0); - } - - // charged hadrons - // G4 particles: of type "baryon", "meson", "nucleus" with non-zero charge - // G3 default value: 0.01 GeV - //gMC ->SetCut("CUTHAD",cut); // cut for charged hadrons - else if (strncmp(cut->GetName(),"CUTHAD",6) == 0 && global) { - fprintf(pFlukaVmcInp,"*\n*Cut for charged hadrons\n"); - fprintf(pFlukaVmcInp,"*Generated from call: SetCut('CUTHAD',cut);\n"); - - // 1.0 = Proton - // 2.0 = Antiproton - fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),1.0,2.0); - - // 13.0 = Positive Pion, Negative Pion, Positive Kaon - // 16.0 = Negative Kaon - fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),13.0,16.0); - - // 20.0 = Negative Sigma - // 21.0 = Positive Sigma - fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),20.0,21.0); - - // 31.0 = Antisigma minus - // 33.0 = Antisigma plus - // 2.0 = step length - fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f%10.1f\n",-cut->Cut(),31.0,33.0,2.0); - - // 36.0 = Negative Xi, Positive Xi, Omega minus - // 39.0 = Antiomega - fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),36.0,39.0); - - // 45.0 = D plus - // 46.0 = D minus - fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),45.0,46.0); - - // 49.0 = D_s plus, D_s minus, Lambda_c plus - // 52.0 = Xi_c plus - fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),49.0,52.0); - - // 54.0 = Xi'_c plus - // 60.0 = AntiXi'_c minus - // 6.0 = step length - fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f%10.1f\n",-cut->Cut(),54.0,60.0,6.0); - - // 57.0 = Antilambda_c minus - // 58.0 = AntiXi_c minus - fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),57.0,58.0); - } - - // muons - // G4 particles: "mu+", "mu-" - // G3 default value: 0.01 GeV - //gMC ->SetCut("CUTMUO",cut); // cut for mu+, mu- - else if (strncmp(cut->GetName(),"CUTMUO",6)== 0 && global) { - fprintf(pFlukaVmcInp,"*\n*Cut for muons\n"); - fprintf(pFlukaVmcInp,"*Generated from call: SetCut('CUTMUO',cut);\n"); - // 10.0 = Muon+ - // 11.0 = Muon- - fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),10.0,11.0); - } - - // - // time of flight cut in seconds - // G4 particles: all - // G3 default value: 0.01 GeV - //gMC ->SetCut("TOFMAX",tofmax); // time of flight cuts in seconds - else if (strncmp(cut->GetName(),"TOFMAX",6) == 0) { - fprintf(pFlukaVmcInp,"*\n*Time of flight cuts in seconds\n"); - fprintf(pFlukaVmcInp,"*Generated from call: SetCut('TOFMAX',tofmax);\n"); - // zero = ignored - // zero = ignored - // -6.0 = lower bound of the particle numbers for which the transport time cut-off and/or the start signal is to be applied - // 64.0 = upper bound of the particle numbers for which the transport time cut-off and/or the start signal is to be applied - fprintf(pFlukaVmcInp,"TIME-CUT %10.4g%10.1f%10.1f%10.1f%10.1f\n",cut->Cut()*1.e9,zero,zero,-6.0,64.0); - } - - else if (global){ - cout << "SetCut for flag=" << cut->GetName() << " value=" << cut->Cut() << " not yet implemented!" << endl; - } - else { - cout << "SetCut for flag=" << cut->GetName() << " value=" << cut->Cut() << " (material specific) not yet implemented!" << endl; - } - - } //end of loop over SetCut calls - + Float_t fLastMaterial = fGeom->GetLastMaterialIndex(); + TFlukaConfigOption::SetStaticInfo(pFlukaVmcInp, 3, fLastMaterial, fGeom); + + TIter next(fUserConfig); + TFlukaConfigOption* proc; + while((proc = (TFlukaConfigOption*)next())) proc->WriteFlukaInputCards(); + // Add START and STOP card fprintf(pFlukaVmcInp,"START %10.1f\n",fEventsPerRun); fprintf(pFlukaVmcInp,"STOP \n"); // Close files - fclose(pFlukaVmcCoreInp); fclose(pFlukaVmcFlukaMat); fclose(pFlukaVmcInp); - + + +// +// Initialisation needed for Cerenkov photon production and transport + TObjArray *matList = GetFlukaMaterials(); + Int_t nmaterial = matList->GetEntriesFast(); + fMaterials = new Int_t[nmaterial+3]; + + for (Int_t im = 0; im < nmaterial; im++) + { + TGeoMaterial* material = dynamic_cast (matList->At(im)); + Int_t idmat = material->GetIndex(); + fMaterials[idmat] = im; + } } // end of InitPhysics diff --git a/TFluka/TFluka.h b/TFluka/TFluka.h index 2a117102635..f40cee50542 100644 --- a/TFluka/TFluka.h +++ b/TFluka/TFluka.h @@ -123,8 +123,7 @@ class TFluka : public TVirtualMC { virtual void SetProcess(const char* flagName, Int_t flagValue, Int_t imed); virtual Bool_t SetCut(const char* cutName, Double_t cutValue); virtual void SetCut(const char* cutName, Double_t cutValue, Int_t imed); - virtual TObjArray* GetListOfProceses() {return fProcesses;} - virtual TObjArray* GetListOfCuts() {return fCuts;} + virtual TObjArray* GetListOfUserConfigs() {return fUserConfig;} virtual Double_t Xsec(char*, Double_t, Int_t, Int_t); // Particle table usage @@ -377,8 +376,7 @@ class TFluka : public TVirtualMC { TGeoMCGeometry *fMCGeo; // Interface to TGeo builder // SetProcess, SetCut and user Scoring dynamic storage - TObjArray* fProcesses; // List of processes - TObjArray* fCuts; // List of cuts + TObjArray* fUserConfig; // List of user physics configuration TObjArray* fUserScore; // List of user scoring options ClassDef(TFluka,1) //C++ interface to Fluka montecarlo diff --git a/TFluka/TFlukaConfigOption.cxx b/TFluka/TFlukaConfigOption.cxx index d2eecff94f6..7c39dff54f6 100644 --- a/TFluka/TFlukaConfigOption.cxx +++ b/TFluka/TFlukaConfigOption.cxx @@ -16,44 +16,624 @@ /* $Id$*/ #include "TFlukaConfigOption.h" +#include "TFlukaMCGeometry.h" +#include "TFluka.h" +#include "TFlukaCerenkov.h" + +#include +#include +#include +#include + +Float_t TFlukaConfigOption::fgMatMin(-1.); +Float_t TFlukaConfigOption::fgMatMax(-1.); +FILE* TFlukaConfigOption::fgFile(0x0); +TFlukaMCGeometry* TFlukaConfigOption::fgGeom(0x0); + +Double_t TFlukaConfigOption::fgDCutValue[11]; +Int_t TFlukaConfigOption::fgDProcessFlag[15]; + + ClassImp(TFlukaConfigOption) TFlukaConfigOption::TFlukaConfigOption() { // Default constructor + fMedium = -1; + fCMatMin = -1.; + fCMatMax = -1.; + Int_t i; + for (i = 0; i < 11; i++) fCutValue[i] = -1.; + for (i = 0; i < 15; i++) fProcessFlag[i] = -1; } -TFlukaConfigOption::TFlukaConfigOption(const char* cutName, Double_t cut) - : TNamed(cutName, "Cut") +TFlukaConfigOption::TFlukaConfigOption(Int_t medium) { // Constructor - fCutValue = cut; - fMedium = -1; + fMedium = medium; + fCMatMin = -1.; + fCMatMax = -1.; + Int_t i; + for (i = 0; i < 11; i++) fCutValue[i] = -1.; + for (i = 0; i < 15; i++) fProcessFlag[i] = -1; } -TFlukaConfigOption::TFlukaConfigOption(const char* cutName, Double_t cut, Int_t imed) - : TNamed(cutName, "Cut") +void TFlukaConfigOption::SetCut(const char* flagname, Double_t val) { - // Constructor - fCutValue = cut; - fMedium = imed; + // Set a cut value + const TString cuts[11] = + {"CUTGAM", "CUTELE", "CUTNEU", "CUTHAD", "CUTMUO", "BCUTE", "BCUTM", "DCUTE", "DCUTM", "PPCUTM", "TOFMAX"}; + Int_t i; + for (i = 0; i < 11; i++) { + if (cuts[i].CompareTo(flagname) == 0) { + fCutValue[i] = val; + if (fMedium == -1) fgDCutValue[i] = val; + break; + } + } } +void TFlukaConfigOption::SetProcess(const char* flagname, Int_t flag) +{ + // Set a process flag + const TString process[15] = + {"DCAY", "PAIR", "COMP", "PHOT", "PFIS", "DRAY", "ANNI", "BREM", "MUNU", "CKOV", + "HADR", "LOSS", "MULS", "RAYL", "STRA"}; + Int_t i; + for (i = 0; i < 15; i++) { + if (process[i].CompareTo(flagname) == 0) { + fProcessFlag[i] = flag; + if (fMedium == -1) fgDProcessFlag[i] = flag; + break; + } + } +} -TFlukaConfigOption::TFlukaConfigOption(const char* procName, Int_t flag) - : TNamed(procName, "Process") +void TFlukaConfigOption::WriteFlukaInputCards() { - // Constructor - fProcessFlag = flag; - fMedium = -1; + // Write the FLUKA input cards for the set of process flags and cuts + // + // + if (fMedium > -1) { + fprintf(fgFile,"*\n*Material specific process and cut settings for #%8d \n", fMedium); + fCMatMin = fMedium; + fCMatMax = fMedium; + } else { + fprintf(fgFile,"*\n*Global process and cut settings \n"); + fCMatMin = fgMatMin; + fCMatMax = fgMatMax; + } + +// +// Handle Process Flags +// + if (fProcessFlag[kDCAY] != -1) ProcessDCAY(); + if (fProcessFlag[kPAIR] != -1) ProcessPAIR(); + if (fProcessFlag[kBREM] != -1) ProcessBREM(); + if (fProcessFlag[kCOMP] != -1) ProcessCOMP(); + if (fProcessFlag[kPHOT] != -1) ProcessPHOT(); + if (fProcessFlag[kPFIS] != -1) ProcessPFIS(); + if (fProcessFlag[kANNI] != -1) ProcessANNI(); + if (fProcessFlag[kMUNU] != -1) ProcessMUNU(); + if (fProcessFlag[kHADR] != -1) ProcessHADR(); + if (fProcessFlag[kMULS] != -1) ProcessMULS(); + if (fProcessFlag[kRAYL] != -1) ProcessRAYL(); + + if (fProcessFlag[kLOSS] != -1 || fProcessFlag[kDRAY] != -1) ProcessLOSS(); + if ((fMedium == -1 && fProcessFlag[kCKOV] > 0) || (fMedium > -1 && fProcessFlag[kCKOV] != -1)) ProcessCKOV(); + +// +// Handle Cuts +// + if (fCutValue[kCUTGAM] >= 0.) ProcessCUTGAM(); + if (fCutValue[kCUTELE] >= 0.) ProcessCUTELE(); + if (fCutValue[kCUTNEU] >= 0.) ProcessCUTNEU(); + if (fCutValue[kCUTHAD] >= 0.) ProcessCUTHAD(); + if (fCutValue[kCUTMUO] >= 0.) ProcessCUTMUO(); + + if (fCutValue[kTOFMAX] >= 0.) ProcessTOFMAX(); } -TFlukaConfigOption::TFlukaConfigOption(const char* procName, Int_t flag, Int_t imed) - : TNamed(procName, "Process") +void TFlukaConfigOption::ProcessDCAY() { - // Constructor - fProcessFlag = flag; - fMedium = imed; + // Process DCAY option + fprintf(fgFile,"*\n* --- DCAY --- Decays. Flag = %5d\n", fProcessFlag[kDCAY]); + if (fProcessFlag[kDCAY] == 0) { + printf("Decays cannot be switched off \n"); + } else { + fprintf(fgFile, "* Decays are on by default\n"); + } +} + + +void TFlukaConfigOption::ProcessPAIR() +{ + // Process PAIR option + fprintf(fgFile,"*\n* --- PAIR --- Pair production by gammas, muons and hadrons. Flag = %5d, PPCUTM = %13.4g \n", + fProcessFlag[kPAIR], fCutValue[kPPCUTM]); + // + // gamma -> e+ e- + // + if (fProcessFlag[kPAIR] > 0) { + fprintf(fgFile,"EMFCUT %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",0., 0., 0., fCMatMin, fCMatMax, 1.); + } else { + fprintf(fgFile,"EMFCUT %10.1f%10.1f%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",0., 0., 1e10, fCMatMin, fCMatMax, 1.); + } + + // + // Direct pair production by Muons and Hadrons + // + // + // Attention ! This card interferes with BREM + // + + if (fProcessFlag[kBREM] == -1 ) fProcessFlag[kBREM] = fgDProcessFlag[kBREM]; + if (fCutValue[kBCUTM] == -1.) fCutValue[kBCUTM] = fgDCutValue[kBCUTM]; + + + Float_t flag = -3.; + if (fProcessFlag[kPAIR] > 0 && fProcessFlag[kBREM] == 0) flag = 1.; + if (fProcessFlag[kPAIR] == 0 && fProcessFlag[kBREM] > 0) flag = 2.; + if (fProcessFlag[kPAIR] > 0 && fProcessFlag[kBREM] > 0) flag = 3.; + if (fProcessFlag[kPAIR] == 0 && fProcessFlag[kBREM] == 0) flag = -3.; + // Flag BREM card as handled + fProcessFlag[kBREM] = -1; + + // + // Energy cut for pair prodution + // + Float_t cutP = fCutValue[kPPCUTM]; + if (fCutValue[kPPCUTM] == -1.) cutP = fgDCutValue[kPPCUTM]; + // In G3 this is the cut on the total energy of the e+e- pair + // In FLUKA the cut is on the kinetic energy of the electron and poistron + cutP = cutP / 2. - 0.51099906e-3; + if (cutP < 0.) cutP = 0.; + // No explicite generation of e+/e- + if (fProcessFlag[kPAIR] == 2) cutP = -1.; + // + // Energy cut for bremsstrahlung + // + Float_t cutB = 0.; + if (flag > 1.) { + fprintf(fgFile,"*\n* +++ BREM --- Bremsstrahlung by muons/hadrons. Flag = %5d, BCUTM = %13.4g \n", + fProcessFlag[kBREM], fCutValue[kBCUTM]); + + cutB = fCutValue[kBCUTM]; + // No explicite production of gammas + if (fProcessFlag[kBREM] == 2) cutB = -1.; + } + + fprintf(fgFile,"PAIRBREM %10.1f%10.4g%10.4g%10.1f%10.1f\n",flag, cutP, cutB, fCMatMin, fCMatMax); +} + + +void TFlukaConfigOption::ProcessBREM() +{ + // Process BREM option + fprintf(fgFile,"*\n* --- BREM --- Bremsstrahlung by e+/- and muons/hadrons. Flag = %5d, BCUTE = %13.4g, BCUTM = %13.4g \n", + fProcessFlag[kBREM], fCutValue[kBCUTE], fCutValue[kBCUTM]); + + // + // e+/- -> e+/- gamma + // + Float_t cutB = fCutValue[kBCUTE]; + if (fCutValue[kBCUTE] == -1.) cutB = fgDCutValue[kBCUTE]; + + + if (fProcessFlag[kBREM] > 0) { + fprintf(fgFile,"EMFCUT %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fELPO-THR\n",cutB, 0., 0., fCMatMin, fCMatMax, 1.); + } else { + fprintf(fgFile,"EMFCUT %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fELPO-THR\n",1.e10, 0., 0., fCMatMin, fCMatMax, 1.); + } + + // + // Bremsstrahlung by muons and hadrons + // + cutB = fCutValue[kBCUTM]; + if (fCutValue[kBCUTM] == -1.) cutB = fgDCutValue[kBCUTM]; + if (fProcessFlag[kBREM] == 2) cutB = -1.; + Float_t flag = 2.; + if (fProcessFlag[kBREM] == 0) flag = -2.; + + fprintf(fgFile,"PAIRBREM %10.1f%10.4g%10.4g%10.1f%10.1f\n", flag, 0., cutB, fCMatMin, fCMatMax); +} + +void TFlukaConfigOption::ProcessCOMP() +{ + // Process COMP option + fprintf(fgFile,"*\n* --- COMP --- Compton scattering Flag = %5d \n", fProcessFlag[kCOMP]); + + // + // Compton scattering + // + + if (fProcessFlag[kCOMP] > 0) { + fprintf(fgFile,"EMFCUT %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",0. , 0., 0., fCMatMin, fCMatMax, 1.); + } else { + fprintf(fgFile,"EMFCUT %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",1.e10, 0., 0., fCMatMin, fCMatMax, 1.); + } +} + +void TFlukaConfigOption::ProcessPHOT() +{ + // Process PHOS option + fprintf(fgFile,"*\n* --- PHOT --- Photoelectric effect. Flag = %5d\n", fProcessFlag[kPHOT]); + + // + // Photoelectric effect + // + + if (fProcessFlag[kPHOT] > 0) { + fprintf(fgFile,"EMFCUT %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",0. , 0., 0., fCMatMin, fCMatMax, 1.); + } else { + fprintf(fgFile,"EMFCUT %10.1f%10.4g%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",0., 1.e10, 0., fCMatMin, fCMatMax, 1.); + } +} + +void TFlukaConfigOption::ProcessANNI() +{ + // Process ANNI option + fprintf(fgFile,"*\n* --- ANNI --- Positron annihilation. Flag = %5d \n", fProcessFlag[kANNI]); + + // + // Positron annihilation + // + + if (fProcessFlag[kANNI] > 0) { + fprintf(fgFile,"EMFCUT %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fANNH-THR\n",0. , 0., 0., fCMatMin, fCMatMax, 1.); + } else { + fprintf(fgFile,"EMFCUT %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fANNH-THR\n",1.e10, 0., 0., fCMatMin, fCMatMax, 1.); + } +} + + +void TFlukaConfigOption::ProcessPFIS() +{ + // Process PFIS option + fprintf(fgFile,"*\n* --- PFIS --- Photonuclear interaction Flag = %5d\n", fProcessFlag[kPFIS]); + + // + // Photonuclear interactions + // + + if (fProcessFlag[kPFIS] > 0) { + 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.); + } else { + fprintf(fgFile,"PHOTONUC %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",-1. , 0., 0., fCMatMin, fCMatMax, 1.); + } +} + +void TFlukaConfigOption::ProcessMUNU() +{ + // Process MUNU option + fprintf(fgFile,"*\n* --- MUNU --- Muon nuclear interaction. Flag = %5d\n", fProcessFlag[kMUNU]); + + // + // Muon nuclear interactions + // + if (fProcessFlag[kMUNU] > 0) { + 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.); + } else { + fprintf(fgFile,"MUPHOTON %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",-1. , 0., 0., fCMatMin, fCMatMax, 1.); + } +} + +void TFlukaConfigOption::ProcessRAYL() +{ + // Process RAYL option + fprintf(fgFile,"*\n* --- RAYL --- Rayleigh Scattering. Flag = %5d\n", fProcessFlag[kRAYL]); + + // + // Rayleigh scattering + // + Int_t nreg; + Int_t* reglist = fgGeom->GetMaterialList(fMedium, nreg); + // + // Loop over regions of a given material + for (Int_t k = 0; k < nreg; k++) { + Float_t ireg = reglist[k]; + if (fProcessFlag[kRAYL] > 0) { + fprintf(fgFile,"EMFRAY %10.1f%10.1f%10.1f%10.1f\n", 1., ireg, ireg, 1.); + } else { + fprintf(fgFile,"EMFRAY %10.1f%10.1f%10.1f%10.1f\n", 3., ireg, ireg, 1.); + } + } +} + +void TFlukaConfigOption::ProcessCKOV() +{ + // Process CKOV option + fprintf(fgFile,"*\n* --- CKOV --- Cerenkov Photon production. %5d\n", fProcessFlag[kCKOV]); + + // + // Cerenkov photon production + // + + TFluka* fluka = (TFluka*) gMC; + TObjArray *matList = fluka->GetFlukaMaterials(); + Int_t nmaterial = matList->GetEntriesFast(); + for (Int_t im = 0; im < nmaterial; im++) + { + TGeoMaterial* material = dynamic_cast (matList->At(im)); + Int_t idmat = material->GetIndex(); +// +// Check if global option + if (fMedium != -1 && idmat != fMedium) continue; + + TFlukaCerenkov* cerenkovProp; + if (!(cerenkovProp = dynamic_cast(material->GetCerenkovProperties()))) continue; + // + // This medium has Cerenkov properties + // + // + if (fMedium == -1 || (fMedium != -1 && fProcessFlag[kCKOV] > 0)) { + // Write OPT-PROD card for each medium + Float_t emin = cerenkovProp->GetMinimumEnergy(); + Float_t emax = cerenkovProp->GetMaximumEnergy(); + fprintf(fgFile, "OPT-PROD %10.4g%10.4g%10.4g%10.4g%10.4g%10.4gCERENKOV\n", emin, emax, 0., + Float_t(idmat), Float_t(idmat), 0.); + // + // Write OPT-PROP card for each medium + // Forcing FLUKA to call user routines (queffc.cxx, rflctv.cxx, rfrndx.cxx) + // + fprintf(fgFile, "OPT-PROP %10.4g%10.4g%10.4g%10.1f%10.1f%10.1fWV-LIMIT\n", + cerenkovProp->GetMinimumWavelength(), cerenkovProp->GetMaximumWavelength(), cerenkovProp->GetMaximumWavelength(), + Float_t(idmat), Float_t(idmat), 0.0); + + if (cerenkovProp->IsMetal()) { + fprintf(fgFile, "OPT-PROP %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fMETAL\n", -100., -100., -100., + Float_t(idmat), Float_t(idmat), 0.0); + } else { + fprintf(fgFile, "OPT-PROP %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n", -100., -100., -100., + Float_t(idmat), Float_t(idmat), 0.0); + } + + + for (Int_t j = 0; j < 3; j++) { + fprintf(fgFile, "OPT-PROP %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f&\n", -100., -100., -100., + Float_t(idmat), Float_t(idmat), 0.0); + } + // Photon detection efficiency user defined + + if (cerenkovProp->IsSensitive()) + fprintf(fgFile, "OPT-PROP %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fSENSITIV\n", -100., -100., -100., + Float_t(idmat), Float_t(idmat), 0.0); + } else { + fprintf(fgFile,"OPT-PROD %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fCERE-OFF\n",0., 0., 0., fCMatMin, fCMatMax, 1.); + } + } +} + + +void TFlukaConfigOption::ProcessHADR() +{ + // Process HADR option + fprintf(fgFile,"*\n* --- HADR --- Hadronic interactions. Flag = %5d\n", fProcessFlag[kHADR]); + + if (fProcessFlag[kHADR] > 0) { + fprintf(fgFile,"*\n*Hadronic interaction is ON by default in FLUKA\n"); + } else { + if (fMedium != -1) printf("Hadronic interactions cannot be switched off material by material !\n"); + fprintf(fgFile,"THRESHOL %10.1f%10.1f%10.1f%10.1e%10.1f\n",0., 0., 0., 1.e10, 0.); + } +} + + + +void TFlukaConfigOption::ProcessMULS() +{ + // Process MULS option + fprintf(fgFile,"*\n* --- MULS --- Muliple Scattering. Flag = %5d\n", fProcessFlag[kMULS]); + // + // Multiple scattering + // + if (fProcessFlag[kMULS] > 0) { + fprintf(fgFile,"*\n*Multiple scattering is ON by default in FLUKA\n"); + } else { + fprintf(fgFile,"MULSOPT %10.1f%10.1f%10.1f%10.1f%10.1f\n",0., 3., 3., fCMatMin, fCMatMax); + } +} + +void TFlukaConfigOption::ProcessLOSS() +{ + // Process LOSS option + fprintf(fgFile,"*\n* --- LOSS --- Ionisation energy loss. Flags: LOSS = %5d, DRAY = %5d, STRA = %5d; Cuts: DCUTE = %13.4g, DCUTM = %13.4g \n", + fProcessFlag[kLOSS], fProcessFlag[kDRAY], fProcessFlag[kSTRA], fCutValue[kDCUTE], fCutValue[kDCUTM]); + // + // Ionisation energy loss + // + // + // Impose consistency + + if (fProcessFlag[kLOSS] == 1 || fProcessFlag[kLOSS] == 3) { + fProcessFlag[kDRAY] = 1; + } else if (fProcessFlag[kLOSS] == 2) { + fProcessFlag[kDRAY] = 0; + fCutValue[kDCUTE] = 1.e10; + fCutValue[kDCUTM] = 1.e10; + } else { + if (fProcessFlag[kDRAY] == 1) { + fProcessFlag[kLOSS] = 1; + } else if (fProcessFlag[kDRAY] == 0) { + fProcessFlag[kLOSS] = 2; + fCutValue[kDCUTE] = 1.e10; + fCutValue[kDCUTM] = 1.e10; + } + } + + if (fCutValue[kDCUTE] == -1.) fCutValue[kDCUTE] = fgDCutValue[kDCUTE]; + if (fCutValue[kDCUTM] == -1.) fCutValue[kDCUTM] = fgDCutValue[kDCUTM]; + + Float_t cutM = fCutValue[kDCUTM]; + + + if (fProcessFlag[kSTRA] == -1) fProcessFlag[kSTRA] = fgDProcessFlag[kSTRA]; + if (fProcessFlag[kSTRA] == 0) fProcessFlag[kSTRA] = 1; + Float_t stra = (Float_t) fProcessFlag[kSTRA]; + + + if (fProcessFlag[kLOSS] == 1 || fProcessFlag[kLOSS] == 3) { +// +// Restricted energy loss fluctuations +// + fprintf(fgFile,"IONFLUCT %10.1f%10.1f%10.1f%10.1f%10.1f\n", 1., 1., stra, fCMatMin, fCMatMax); + fprintf(fgFile,"DELTARAY %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n", cutM, 0., 0., fCMatMin, fCMatMax, 1.); + } else if (fProcessFlag[kLOSS] == 4) { +// +// No fluctuations +// + fprintf(fgFile,"IONFLUCT %10.1f%10.1f%10.1f%10.1f%10.1f\n",-1., -1., stra, fCMatMin, fCMatMax); + fprintf(fgFile,"DELTARAY %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n", 1.e10, 0., 0., fCMatMin, fCMatMax, 1.); + } else { +// +// Full fluctuations +// + fprintf(fgFile,"IONFLUCT %10.1f%10.1f%10.1f%10.1f%10.1f\n",1., 1., stra, fCMatMin, fCMatMax); + fprintf(fgFile,"DELTARAY %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n", 1.e10, 0., 0., fCMatMin, fCMatMax, 1.); + } +} + + +void TFlukaConfigOption::ProcessCUTGAM() +{ +// Cut on gammas +// + fprintf(fgFile,"*\n*Cut for Gammas. CUTGAM = %13.4g\n", fCutValue[kCUTGAM]); + if (fMedium == -1) { + fprintf(fgFile,"EMFCUT %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", + 0., fCutValue[kCUTGAM], 0., 0., Float_t(fgGeom->NofVolumes()), 1.); + } else { + Int_t nreg, *reglist; + Float_t ireg; + reglist = fgGeom->GetMaterialList(fMedium, nreg); + // Loop over regions of a given material + for (Int_t k = 0; k < nreg; k++) { + ireg = reglist[k]; + fprintf(fgFile,"EMFCUT %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", 0.,fCutValue[kCUTGAM], 0., ireg, ireg, 1.); + } + } +} + +void TFlukaConfigOption::ProcessCUTELE() +{ +// Cut on e+/e- +// + fprintf(fgFile,"*\n*Cut for e+/e-. CUTELE = %13.4g\n", fCutValue[kCUTELE]); + if (fMedium == -1) { + fprintf(fgFile,"EMFCUT %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", + -fCutValue[kCUTELE], 0., 0., 0., Float_t(fgGeom->NofVolumes()), 1.); + } else { + Int_t nreg, *reglist; + Float_t ireg; + reglist = fgGeom->GetMaterialList(fMedium, nreg); + // Loop over regions of a given material + for (Int_t k = 0; k < nreg; k++) { + ireg = reglist[k]; + fprintf(fgFile,"EMFCUT %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", -fCutValue[kCUTELE], 0., 0., ireg, ireg, 1.); + } + } +} + +void TFlukaConfigOption::ProcessCUTNEU() +{ + // Cut on neutral hadrons + fprintf(fgFile,"*\n*Cut for neutal hadrons. CUTNEU = %13.4g\n", fCutValue[kCUTNEU]); + if (fMedium == -1) { + Float_t cut = fCutValue[kCUTNEU]; + // 8.0 = Neutron + // 9.0 = Antineutron + fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 8.0, 9.0); + // 12.0 = Kaon zero long + fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 12.0, 12.0); + // 17.0 = Lambda, 18.0 = Antilambda + // 19.0 = Kaon zero short + fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 17.0, 19.0); + // 22.0 = Sigma zero, Pion zero, Kaon zero + // 25.0 = Antikaon zero + fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 22.0, 25.0); + // 32.0 = Antisigma zero + fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 32.0, 32.0); + // 34.0 = Xi zero + // 35.0 = AntiXi zero + fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 34.0, 35.0); + // 47.0 = D zero + // 48.0 = AntiD zero + fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 47.0, 48.0); + // 53.0 = Xi_c zero + fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 53.0, 53.0); + // 55.0 = Xi'_c zero + // 56.0 = Omega_c zero + fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 55.0, 56.0); + // 59.0 = AntiXi_c zero + fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 59.0, 59.0); + // 61.0 = AntiXi'_c zero + // 62.0 = AntiOmega_c zero + fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 61.0, 62.0); + } else { + printf("Cuts on neutral hadrons material by material not yet implemented !\n"); + } +} + +void TFlukaConfigOption::ProcessCUTHAD() +{ + // Cut on charged hadrons + fprintf(fgFile,"*\n*Cut for charge hadrons. CUTHAD = %13.4g\n", fCutValue[kCUTHAD]); + if (fMedium == -1) { + Float_t cut = fCutValue[kCUTHAD]; + // 1.0 = Proton + // 2.0 = Antiproton + fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 1.0, 2.0); + // 13.0 = Positive Pion, Negative Pion, Positive Kaon + // 16.0 = Negative Kaon + fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 13.0, 16.0); + // 20.0 = Negative Sigma + // 21.0 = Positive Sigma + fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 20.0, 21.0); + // 31.0 = Antisigma minus + // 33.0 = Antisigma plus + fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 31.0, 31.0); + fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 33.0, 33.0); + // 36.0 = Negative Xi, Positive Xi, Omega minus + // 39.0 = Antiomega + fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 36.0, 39.0); + // 45.0 = D plus + // 46.0 = D minus + fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 45.0, 46.0); + // 49.0 = D_s plus, D_s minus, Lambda_c plus + // 52.0 = Xi_c plus + fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 49.0, 52.0); + // 54.0 = Xi'_c plus + // 60.0 = AntiXi'_c minus + fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 54.0, 54.0); + fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 60.0, 60.0); + // 57.0 = Antilambda_c minus + // 58.0 = AntiXi_c minus + fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n", -cut, 57.0, 58.0); + } else { + printf("Cuts on charged hadrons material by material not yet implemented !\n"); + } +} + +void TFlukaConfigOption::ProcessCUTMUO() +{ + // Cut on muons + fprintf(fgFile,"*\n*Cut for muons. CUTMUO = %13.4g\n", fCutValue[kCUTMUO]); + Float_t cut = fCutValue[kCUTMUO]; + if (fMedium == -1) { + fprintf(fgFile,"PART-THR %10.4g%10.1f%10.1f\n",-cut, 10.0, 11.0); + } else { + printf("Cuts on muons material by material not yet implemented !\n"); + } + + +} + +void TFlukaConfigOption::ProcessTOFMAX() +{ + // Cut time of flight + Float_t cut = fCutValue[kTOFMAX]; + fprintf(fgFile,"*\n*Cut on time of flight. TOFMAX = %13.4g\n", fCutValue[kTOFMAX]); + fprintf(fgFile,"TIME-CUT %10.4g%10.1f%10.1f%10.1f%10.1f\n",cut*1.e9,0.,0.,-6.0,64.0); } diff --git a/TFluka/TFlukaConfigOption.h b/TFluka/TFlukaConfigOption.h index 3153d58aa44..1b0fe05e351 100644 --- a/TFluka/TFlukaConfigOption.h +++ b/TFluka/TFlukaConfigOption.h @@ -18,36 +18,73 @@ // // /////////////////////////////////////////////////////////////////////////////// // -// -// -#include +// +#include -class TFlukaConfigOption : public TNamed +typedef enum {kDCAY, kPAIR, kCOMP, kPHOT, kPFIS, kDRAY, kANNI, kBREM, + kMUNU, kCKOV, kHADR, kLOSS, kMULS, kRAYL, kSTRA} FlukaProcessOption_t; +typedef enum {kCUTGAM, kCUTELE, kCUTNEU, kCUTHAD, kCUTMUO, kBCUTE, kBCUTM, kDCUTE, kDCUTM, kPPCUTM, kTOFMAX} FlukaCutOption_t; +class TFlukaMCGeometry; + +class TFlukaConfigOption : public TObject { public: // Constructors TFlukaConfigOption(); - TFlukaConfigOption(const char* cutName, Double_t cut); - TFlukaConfigOption(const char* cutName, Double_t cut, Int_t imed); - TFlukaConfigOption(const char* procName, Int_t flag); - TFlukaConfigOption(const char* procName, Int_t flag, Int_t imed); + TFlukaConfigOption(Int_t imed); // Getters - Double_t Cut() const {return fCutValue;} - Int_t Flag() const {return fProcessFlag;} - Int_t Medium() {return fMedium;} - Bool_t HasMediumAssigned() {return (fMedium > -1);} + Double_t Cut(FlukaCutOption_t i) const {return fCutValue[i];} + Int_t Flag(FlukaProcessOption_t i) const {return fProcessFlag[i];} + Int_t Medium() const {return fMedium;} // Setters - void SetCut(Double_t val) {fCutValue = val;} - void SetFlag(Int_t val) {fProcessFlag = val;} - void SetMedium(Int_t imed) {fMedium = imed;} - + void SetCut(const char* flagname, Double_t val); + void SetProcess(const char* flagname, Int_t flagValue); + void SetMedium(Int_t imed) {fMedium = imed;} + // + void WriteFlukaInputCards(); + void ProcessDCAY(); + void ProcessPAIR(); + void ProcessBREM(); + void ProcessCOMP(); + void ProcessPHOT(); + void ProcessANNI(); + void ProcessPFIS(); + void ProcessMUNU(); + void ProcessRAYL(); + void ProcessCKOV(); + void ProcessHADR(); + void ProcessMULS(); + void ProcessLOSS(); + void ProcessSTRA(); + + + void ProcessCUTGAM(); + void ProcessCUTELE(); + void ProcessCUTNEU(); + void ProcessCUTHAD(); + void ProcessCUTMUO(); + void ProcessTOFMAX(); + + // + static void SetStaticInfo(FILE* file, Float_t matMin, Float_t matMax, TFlukaMCGeometry* geom) + {fgFile = file; fgMatMin = matMin; fgMatMax = matMax; fgGeom = geom;} protected: - Double_t fCutValue; // User cut - Int_t fProcessFlag; // User flag assigned to processes + Double_t fCutValue[11]; // User cut + Int_t fProcessFlag[15]; // User flag assigned to processes Int_t fMedium; // Materials assigned to user settings + Float_t fCMatMin; // Minimum material number used for current card + Float_t fCMatMax; // Maximum material number used for current card + + // static + static Double_t fgDCutValue[11]; // User default cut + static Int_t fgDProcessFlag[15]; // User default flag assigned to processes + static Float_t fgMatMin; // Minimum material number + static Float_t fgMatMax; // Maximum meterial number + static FILE* fgFile; // Output file + static TFlukaMCGeometry* fgGeom; // Pointer to geometry ClassDef(TFlukaConfigOption, 1) // Fluka Configuration Option }; -- 2.31.1