X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=TFluka%2FTFlukaConfigOption.cxx;h=d6dde3fc6dfc07e0b0389c6fd3fdf79ee0373fab;hp=7c39dff54f668fd31818b112afd7b7cfe7ef3a7f;hb=8134b16b635b1138c79cde7f0924a127bf487e21;hpb=fb2cbbecf3f0d2e56028a85236a6d54e670dfd8e diff --git a/TFluka/TFlukaConfigOption.cxx b/TFluka/TFlukaConfigOption.cxx index 7c39dff54f6..d6dde3fc6df 100644 --- a/TFluka/TFlukaConfigOption.cxx +++ b/TFluka/TFlukaConfigOption.cxx @@ -21,9 +21,13 @@ #include "TFlukaCerenkov.h" #include +#include #include #include #include +#include +#include +#include Float_t TFlukaConfigOption::fgMatMin(-1.); Float_t TFlukaConfigOption::fgMatMax(-1.); @@ -81,6 +85,7 @@ void TFlukaConfigOption::SetProcess(const char* flagname, Int_t 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) { @@ -96,8 +101,55 @@ void TFlukaConfigOption::WriteFlukaInputCards() // 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); + // + // Check if global option or medium specific + + Bool_t mediumIsSensitive = kFALSE; + TGeoMedium* med = 0x0; + TGeoMedium* medium = 0x0; + TGeoMaterial* mat = 0x0; + + if (fMedium != -1) { + TFluka* fluka = (TFluka*) gMC; + fMedium = fgGeom->GetFlukaMaterial(fMedium); + // + // Check if material is actually used + Int_t* reglist; + Int_t nreg; + reglist = fgGeom->GetMaterialList(fMedium, nreg); + if (nreg == 0) { + printf("Material not used !\n"); + return; + } + // + // Find material + TObjArray *matList = fluka->GetFlukaMaterials(); + Int_t nmaterial = matList->GetEntriesFast(); + fCMaterial = 0; + for (Int_t im = 0; im < nmaterial; im++) + { + fCMaterial = dynamic_cast (matList->At(im)); + Int_t idmat = fCMaterial->GetIndex(); + if (idmat == fMedium) break; + } // materials + // + // Find medium + TList *medlist = gGeoManager->GetListOfMedia(); + TIter next(medlist); + while((med = (TGeoMedium*)next())) + { + mat = med->GetMaterial(); + if (mat->GetIndex() == fMedium) { + medium = med; + break; + } + } // media + // + // Check if sensitive + if (medium->GetParam(0) != 0.) mediumIsSensitive = kTRUE; + + + fprintf(fgFile,"*\n*Material specific process and cut settings for #%8d %s\n", fMedium, fCMaterial->GetName()); fCMatMin = fMedium; fCMatMax = fMedium; } else { @@ -108,7 +160,17 @@ void TFlukaConfigOption::WriteFlukaInputCards() // // Handle Process Flags +// +// +// First make sure that all cuts are taken into account + if (DefaultProcessFlag(kPAIR) > 0 && fProcessFlag[kPAIR] == -1 && (fCutValue[kCUTELE] >= 0. || fCutValue[kPPCUTM] >= 0.)) + fProcessFlag[kPAIR] = DefaultProcessFlag(kPAIR); + if (DefaultProcessFlag(kBREM) > 0 && fProcessFlag[kBREM] == -1 && (fCutValue[kBCUTE] >= 0. || fCutValue[kBCUTM] >= 0.)) + fProcessFlag[kBREM] = DefaultProcessFlag(kBREM); + if (DefaultProcessFlag(kDRAY) > 0 && fProcessFlag[kDRAY] == -1 && (fCutValue[kDCUTE] >= 0. || fCutValue[kDCUTM] >= 0.)) + fProcessFlag[kDRAY] = DefaultProcessFlag(kDRAY); // +// if (fProcessFlag[kDCAY] != -1) ProcessDCAY(); if (fProcessFlag[kPAIR] != -1) ProcessPAIR(); if (fProcessFlag[kBREM] != -1) ProcessBREM(); @@ -132,8 +194,13 @@ void TFlukaConfigOption::WriteFlukaInputCards() if (fCutValue[kCUTNEU] >= 0.) ProcessCUTNEU(); if (fCutValue[kCUTHAD] >= 0.) ProcessCUTHAD(); if (fCutValue[kCUTMUO] >= 0.) ProcessCUTMUO(); - +// +// Time of flight if (fCutValue[kTOFMAX] >= 0.) ProcessTOFMAX(); + +// +// Tracking precission + if (mediumIsSensitive) ProcessSensitiveMedium(); } void TFlukaConfigOption::ProcessDCAY() @@ -151,13 +218,13 @@ void TFlukaConfigOption::ProcessDCAY() 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]); + fprintf(fgFile,"*\n* --- PAIR --- Pair production by gammas, muons and hadrons. Flag = %5d, PPCUTM = %13.4g, PPCUTE = %13.4g\n", + fProcessFlag[kPAIR], fCutValue[kCUTELE], 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.); + fprintf(fgFile,"EMFCUT %10.1f%10.1f%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",0., 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.); } @@ -179,7 +246,7 @@ void TFlukaConfigOption::ProcessPAIR() 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; + fProcessFlag[kBREM] = - fProcessFlag[kBREM]; // // Energy cut for pair prodution @@ -222,7 +289,7 @@ void TFlukaConfigOption::ProcessBREM() if (fCutValue[kBCUTE] == -1.) cutB = fgDCutValue[kBCUTE]; - if (fProcessFlag[kBREM] > 0) { + if (TMath::Abs(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.); @@ -266,7 +333,7 @@ void TFlukaConfigOption::ProcessPHOT() // 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.); + fprintf(fgFile,"EMFCUT %10.4g%10.4g%10.4g%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.); } @@ -382,24 +449,25 @@ void TFlukaConfigOption::ProcessCKOV() 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); - } - + + 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 - + + + // 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); + // Material has a reflective surface + 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-PROD %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fCERE-OFF\n",0., 0., 0., fCMatMin, fCMatMax, 1.); } @@ -479,7 +547,7 @@ void TFlukaConfigOption::ProcessLOSS() // 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.); + 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 @@ -504,6 +572,7 @@ void TFlukaConfigOption::ProcessCUTGAM() 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; @@ -514,6 +583,14 @@ void TFlukaConfigOption::ProcessCUTGAM() fprintf(fgFile,"EMFCUT %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", 0.,fCutValue[kCUTGAM], 0., ireg, ireg, 1.); } } + + // Transport production cut used for pemf + // + // FUDGEM paramter. The parameter takes into account th contribution of atomic electrons to multiple scattering. + // For production and transport cut-offs larger than 100 keV it must be set = 1.0, while in the keV region it must be + Float_t parFudgem = (fCutValue[kCUTGAM] > 1.e-4)? 1.0 : 0.0 ; + fprintf(fgFile,"EMFCUT %10.4g%10.4g%10.1f%10.1f%10.1f%10.1fPROD-CUT\n", + 0., fCutValue[kCUTGAM], parFudgem, fCMatMin, fCMatMax, 1.); } void TFlukaConfigOption::ProcessCUTELE() @@ -534,17 +611,63 @@ void TFlukaConfigOption::ProcessCUTELE() fprintf(fgFile,"EMFCUT %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", -fCutValue[kCUTELE], 0., 0., ireg, ireg, 1.); } } + // Transport production cut used for pemf + // + // FUDGEM paramter. The parameter takes into account th contribution of atomic electrons to multiple scattering. + // For production and transport cut-offs larger than 100 keV it must be set = 1.0, while in the keV region it must be + Float_t parFudgem = (fCutValue[kCUTELE] > 1.e-4)? 1.0 : 0.0; + fprintf(fgFile,"EMFCUT %10.4g%10.4g%10.1f%10.1f%10.1f%10.1fPROD-CUT\n", + -fCutValue[kCUTELE], 0., parFudgem, fCMatMin, fCMatMax, 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); + // Cut on neutral hadrons + fprintf(fgFile,"*\n*Cut for neutral hadrons. CUTNEU = %13.4g\n", fCutValue[kCUTNEU]); + + // Cut on neutral hadrons + fprintf(fgFile,"*\n*Cut for neutral hadrons. CUTNEU = %13.4g\n", fCutValue[kCUTNEU]); + + // Energy group structure of the 72-neutron ENEA data set: + const Float_t neutronGroupUpLimit[72] = { + 1.9600E-02, 1.7500E-02, 1.4918E-02, 1.3499E-02, 1.2214E-02, 1.1052E-02, 1.0000E-02, 9.0484E-03, + 8.1873E-03, 7.4082E-03, 6.7032E-03, 6.0653E-03, 5.4881E-03, 4.9659E-03, 4.4933E-03, 4.0657E-03, + 3.6788E-03, 3.3287E-03, 3.0119E-03, 2.7253E-03, 2.4660E-03, 2.2313E-03, 2.0190E-03, 1.8268E-03, + 1.6530E-03, 1.4957E-03, 1.3534E-03, 1.2246E-03, 1.1080E-03, 1.0026E-03, 9.0718E-04, 8.2085E-04, + 7.4274E-04, 6.0810E-04, 4.9787E-04, 4.0762E-04, 3.3373E-04, 2.7324E-04, 2.2371E-04, 1.8316E-04, + 1.4996E-04, 1.2277E-04, 8.6517E-05, 5.2475E-05, 3.1828E-05, 2.1852E-05, 1.5034E-05, 1.0332E-05, + 7.1018E-06, 4.8809E-06, 3.3546E-06, 2.3054E-06, 1.5846E-06, 1.0446E-06, 6.8871E-07, 4.5400E-07, + 2.7537E-07, 1.6702E-07, 1.0130E-07, 6.1442E-08, 3.7267E-08, 2.2603E-08, 1.5535E-08, 1.0677E-08, + 7.3375E-09, 5.0435E-09, 3.4662E-09, 2.3824E-09, 1.6374E-09, 1.1254E-09, 6.8257E-10, 4.1400E-10 + }; + + Float_t cut = fCutValue[kCUTNEU]; + // + // 8.0 = Neutron + // 9.0 = Antineutron + // Find the FLUKA neutron group corresponding to the cut + // + Float_t neutronCut = cut; + Int_t groupCut = 1; // if cut is > 19.6 MeV not low energy neutron transport is done + if (neutronCut < 0.0196) { + neutronCut = 0.0196; + // Search the group cutoff for the energy cut + Int_t i; + for( i=0; i<72; i++ ) { + if (cut > neutronGroupUpLimit[i]) break; + } + groupCut = i+1; + } + + 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", -neutronCut, 8.0, 9.0); + fprintf(fgFile,"LOW-BIAS %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", + Float_t(groupCut), 73.0, 0.95, 2., Float_t(fgGeom->NofVolumes()), 1.); + // + // // 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 @@ -572,7 +695,19 @@ void TFlukaConfigOption::ProcessCUTNEU() // 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"); + 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,"LOW-BIAS %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", + Float_t(groupCut), 73.0, 0.95, ireg, ireg, 1.); + } + + Warning("ProcessCUTNEU", + "Material #%4d %s: Cut on neutral hadrons (Ekin > %9.3e) material by material only implemented for low-energy neutrons !\n", + fMedium, fCMaterial->GetName(), cut); } } @@ -580,8 +715,8 @@ void TFlukaConfigOption::ProcessCUTHAD() { // Cut on charged hadrons fprintf(fgFile,"*\n*Cut for charge hadrons. CUTHAD = %13.4g\n", fCutValue[kCUTHAD]); + Float_t cut = 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); @@ -612,7 +747,9 @@ void TFlukaConfigOption::ProcessCUTHAD() // 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"); + Warning("ProcessCUTHAD", + "Material #%4d %s: Cut on charged hadrons (Ekin > 9.3e) material by material not yet implemented !\n", + fMedium, fCMaterial->GetName(), cut); } } @@ -624,7 +761,8 @@ void TFlukaConfigOption::ProcessCUTMUO() 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"); + Warning("ProcessCUTMUO", "Material #%4d %s: Cut on muons (Ekin > %9.3e) material by material not yet implemented !\n", + fMedium, fCMaterial->GetName(), cut); } @@ -637,3 +775,18 @@ void TFlukaConfigOption::ProcessTOFMAX() 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); } + +void TFlukaConfigOption::ProcessSensitiveMedium() +{ + // + // Special options for sensitive media + // + + fprintf(fgFile,"*\n*Options for sensitive medium \n"); + // + // EMFFIX + fprintf(fgFile,"EMFFIX %10.1f%10.3f%10.1f%10.1f%10.1f%10.1f\n", fCMatMin, 0.05, 0., 0., 0., 0.); + // + // FLUKAFIX + fprintf(fgFile,"FLUKAFIX %10.3f %10.3f\n", 0.05, fCMatMin); +}