X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TFluka%2FTFluka.cxx;h=d55b2165cea127b0e77505291213f5450d5141f1;hb=a9d747801228ed288829aab75a516f3bf6552b00;hp=e7828a205e86c13b64c29e048e66dd3e698c11cb;hpb=83ec523afd647569fd0272a2e80c4587d7ad0427;p=u%2Fmrichter%2FAliRoot.git diff --git a/TFluka/TFluka.cxx b/TFluka/TFluka.cxx index e7828a205e8..d55b2165cea 100644 --- a/TFluka/TFluka.cxx +++ b/TFluka/TFluka.cxx @@ -35,7 +35,7 @@ #include "TCallf77.h" //For the fortran calls #include "Fdblprc.h" //(DBLPRC) fluka common #include "Fepisor.h" //(EPISOR) fluka common -#include "Ffinuc.h" //(FINUC) fluka common +#include "Ffinuc.h" //(FINUC) fluka common #include "Fiounit.h" //(IOUNIT) fluka common #include "Fpaprop.h" //(PAPROP) fluka common #include "Fpart.h" //(PART) fluka common @@ -43,7 +43,9 @@ #include "Fpaprop.h" //(PAPROP) fluka common #include "Ffheavy.h" //(FHEAVY) fluka common #include "Fopphst.h" //(OPPHST) fluka common -#include "Fstack.h" //(STACK) fluka common +#include "Fstack.h" //(STACK) fluka common +#include "Fstepsz.h" //(STEPSZ) fluka common +#include "Fopphst.h" //(OPPHST) fluka common #include "TVirtualMC.h" #include "TMCProcess.h" @@ -54,21 +56,27 @@ #include "TGeoMCGeometry.h" #include "TFlukaCerenkov.h" #include "TFlukaConfigOption.h" +#include "TFlukaScoringOption.h" #include "TLorentzVector.h" +#include "TArrayI.h" // Fluka methods that may be needed. #ifndef WIN32 # define flukam flukam_ # define fluka_openinp fluka_openinp_ +# define fluka_openout fluka_openout_ # define fluka_closeinp fluka_closeinp_ # define mcihad mcihad_ # define mpdgha mpdgha_ +# define newplo newplo_ #else # define flukam FLUKAM # define fluka_openinp FLUKA_OPENINP +# define fluka_openout FLUKA_OPENOUT # define fluka_closeinp FLUKA_CLOSEINP # define mcihad MCIHAD # define mpdgha MPDGHA +# define newplo NEWPLO #endif extern "C" @@ -77,7 +85,9 @@ extern "C" // Prototypes for FLUKA functions // void type_of_call flukam(const int&); + void type_of_call newplo(); void type_of_call fluka_openinp(const int&, DEFCHARA); + void type_of_call fluka_openout(const int&, DEFCHARA); void type_of_call fluka_closeinp(const int&); int type_of_call mcihad(const int&); int type_of_call mpdgha(const int&); @@ -96,8 +106,7 @@ TFluka::TFluka() :TVirtualMC(), fVerbosityLevel(0), fInputFileName(""), - fProcesses(0), - fCuts(0), + fUserConfig(0), fUserScore(0) { // @@ -112,6 +121,9 @@ TFluka::TFluka() fDummyBoundary = 0; fFieldFlag = 1; fStopped = 0; + fStopEvent = 0; + fStopRun = 0; + fNEvent = 0; } //______________________________________________________________________________ @@ -122,24 +134,29 @@ 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 - if (fVerbosityLevel >=3) - cout << "<== TFluka::TFluka(" << title << ") constructor called." << endl; - + if (fVerbosityLevel >=3) + cout << "<== TFluka::TFluka(" << title << ") constructor called." << endl; + SetCoreInputFileName(); + SetInputFileName(); + SetGeneratePemf(kFALSE); fNVolumes = 0; fCurrentFlukaRegion = -1; fDummyBoundary = 0; fFieldFlag = 1; fGeneratePemf = kFALSE; fMCGeo = new TGeoMCGeometry("MCGeo", "TGeo Implementation of VirtualMCGeometry", kTRUE); - fGeom = new TFlukaMCGeometry("geom", "ALICE geometry"); + fGeom = new TFlukaMCGeometry("geom", "FLUKA VMC Geometry"); if (verbosity > 2) fGeom->SetDebugMode(kTRUE); fMaterials = 0; fStopped = 0; + fStopEvent = 0; + fStopRun = 0; + fNEvent = 0; + PrintHeader(); } //______________________________________________________________________________ @@ -151,17 +168,15 @@ TFluka::~TFluka() { delete fGeom; delete fMCGeo; - if (fCuts) { - fCuts->Delete(); - delete fCuts; + if (fUserConfig) { + fUserConfig->Delete(); + delete fUserConfig; } - - if (fProcesses) { - fProcesses->Delete(); - delete fProcesses; + + if (fUserScore) { + fUserScore->Delete(); + delete fUserScore; } - - } // @@ -180,14 +195,19 @@ void TFluka::Init() { gGeoManager->SetTopVolume(top); gGeoManager->CloseGeometry("di"); gGeoManager->DefaultColors(); // to be removed + + // Now we have TGeo geometry created and we have to patch FlukaVmc.inp + // with the material mapping file FlukaMat.inp + fNVolumes = fGeom->NofVolumes(); fGeom->CreateFlukaMatFile("flukaMat.inp"); if (fVerbosityLevel >=3) { printf("== Number of volumes: %i\n ==", fNVolumes); cout << "\t* InitPhysics() - Prepare input file to be called" << endl; - } - // now we have TGeo geometry created and we have to patch alice.inp - // with the material mapping file FlukaMat.inp + } + + fApplication->InitGeometry(); + } @@ -211,8 +231,34 @@ void TFluka::BuildPhysics() { if (fVerbosityLevel >=3) cout << "==> TFluka::BuildPhysics() called." << endl; -// Prepare input file with the current physics settings + + + if (fVerbosityLevel >=3) { + TList *medlist = gGeoManager->GetListOfMedia(); + TIter next(medlist); + TGeoMedium* med = 0x0; + TGeoMaterial* mat = 0x0; + Int_t ic = 0; + + while((med = (TGeoMedium*)next())) + { + mat = med->GetMaterial(); + printf("Medium %5d %12s %5d %5d\n", ic, (med->GetName()), med->GetId(), mat->GetIndex()); + ic++; + } + } + + // + // At this stage we have the information on materials and cuts available. + // Now create the pemf file + + if (fGeneratePemf) fGeom->CreatePemfFile(); + + // + // Prepare input file with the current physics settings + InitPhysics(); + cout << "\t* InitPhysics() - Prepare input file was called" << endl; if (fVerbosityLevel >=2) @@ -223,7 +269,9 @@ void TFluka::BuildPhysics() { if (fVerbosityLevel >=2) cout << "\t* Opening file " << fInputFileName << endl; const char* fname = fInputFileName; + fluka_openinp(lunin, PASSCHARA(fname)); + fluka_openout(11, PASSCHARA("fluka.out")); if (fVerbosityLevel >=2) cout << "\t* Calling flukam..." << endl; @@ -238,7 +286,6 @@ void TFluka::BuildPhysics() { if (fVerbosityLevel >=3) cout << "<== TFluka::Init() called." << endl; - if (fVerbosityLevel >=3) cout << "<== TFluka::BuildPhysics() called." << endl; } @@ -248,13 +295,23 @@ void TFluka::ProcessEvent() { // // Process one event // - if (fVerbosityLevel >=3) - cout << "==> TFluka::ProcessEvent() called." << endl; - fApplication->GeneratePrimaries(); - EPISOR.lsouit = true; - flukam(1); - if (fVerbosityLevel >=3) - cout << "<== TFluka::ProcessEvent() called." << endl; + if (fStopRun) { + printf("User Run Abortion: No more events handled !\n"); + fNEvent += 1; + return; + } + + if (fVerbosityLevel >=3) + cout << "==> TFluka::ProcessEvent() called." << endl; + fApplication->GeneratePrimaries(); + EPISOR.lsouit = true; + flukam(1); + if (fVerbosityLevel >=3) + cout << "<== TFluka::ProcessEvent() called." << endl; + // + // Increase event number + // + fNEvent += 1; } //______________________________________________________________________________ @@ -272,7 +329,6 @@ Bool_t TFluka::ProcessRun(Int_t nevent) { cout << "\t* Calling flukam again..." << endl; } - fApplication->InitGeometry(); Int_t todo = TMath::Abs(nevent); for (Int_t ev = 0; ev < todo; ev++) { fApplication->BeginEvent(); @@ -283,6 +339,9 @@ Bool_t TFluka::ProcessRun(Int_t nevent) { if (fVerbosityLevel >=3) cout << "<== TFluka::ProcessRun(" << nevent << ") called." << endl; + // Write fluka specific scoring output + newplo(); + return kTRUE; } @@ -350,6 +409,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) { @@ -366,6 +426,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)); @@ -501,7 +563,8 @@ void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat, Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd, Double_t stemax, Double_t deemax, Double_t epsil, Double_t stmin, Float_t* ubuf, Int_t nbuf) { - // + // Define a medium + // kmed = gGeoManager->GetListOfMedia()->GetSize()+1; fMCGeo->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax, epsil, stmin, ubuf, nbuf); @@ -512,7 +575,8 @@ void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat, Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd, Double_t stemax, Double_t deemax, Double_t epsil, Double_t stmin, Double_t* ubuf, Int_t nbuf) { - // + // Define a medium + // kmed = gGeoManager->GetListOfMedia()->GetSize()+1; fMCGeo->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax, epsil, stmin, ubuf, nbuf); @@ -531,9 +595,17 @@ void TFluka::Matrix(Int_t& krot, Double_t thetaX, Double_t phiX, void TFluka::Gstpar(Int_t itmed, const char* param, Double_t parval) { // // - - if (fVerbosityLevel >=3) printf("Gstpar called with %6d %5s %12.4e %6d\n", itmed, param, parval, fGeom->GetFlukaMaterial(itmed)); - +// Check if material is used + if (fVerbosityLevel >= 3) + printf("Gstpar called with %6d %5s %12.4e %6d\n", itmed, param, parval, fGeom->GetFlukaMaterial(itmed)); + Int_t* reglist; + Int_t nreg; + reglist = fGeom->GetMaterialList(fGeom->GetFlukaMaterial(itmed), nreg); + if (nreg == 0) { + return; + } + +// Bool_t process = kFALSE; if (strncmp(param, "DCAY", 4) == 0 || strncmp(param, "PAIR", 4) == 0 || @@ -672,14 +744,42 @@ void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t* ppckov, medium->SetCerenkovProperties(cerenkovProperties); } +void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t* ppckov, + Float_t* absco, Float_t* effic, Float_t* rindex, Float_t* rfl) { +// +// Set Cerenkov properties for medium itmed +// +// npckov: number of sampling points +// ppckov: energy values +// absco: absorption length +// effic: quantum efficiency +// rindex: refraction index +// rfl: reflectivity for boundary to medium itmed +// +// +// Create object holding Cerenkov properties +// + TFlukaCerenkov* cerenkovProperties = new TFlukaCerenkov(npckov, ppckov, absco, effic, rindex, rfl); +// +// Pass object to medium + TGeoMedium* medium = gGeoManager->GetMedium(itmed); + medium->SetCerenkovProperties(cerenkovProperties); +} + + //______________________________________________________________________________ void TFluka::SetCerenkov(Int_t /*itmed*/, Int_t /*npckov*/, Double_t * /*ppckov*/, Double_t * /*absco*/, Double_t * /*effic*/, Double_t * /*rindex*/) { // -// Not implemented with TGeo - what G4 did ? Any FLUKA card generated? - Warning("SetCerenkov", "Not implemented with TGeo"); +// Double_t version not implemented } - + +void TFluka::SetCerenkov(Int_t /*itmed*/, Int_t /*npckov*/, Double_t* /*ppckov*/, + Double_t* /*absco*/, Double_t* /*effic*/, Double_t* /*rindex*/, Double_t* /*rfl*/) { +// +// // Double_t version not implemented +} + // Euclid //______________________________________________________________________________ void TFluka::WriteEuclid(const char* /*fileName*/, const char* /*topVol*/, @@ -781,37 +881,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; } @@ -820,8 +914,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); } //______________________________________________________________________________ @@ -830,26 +935,27 @@ Bool_t TFluka::SetCut(const char* cutName, Double_t cutValue) // Set user cut value // // -// Update if already in the list -// + SetCut(cutName, cutValue, -1); + return kTRUE; +} - 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; +void TFluka::SetUserScoring(const char* option, Int_t npr, char* outfile, Float_t* what) +{ +// +// Adds a user scoring option to the list +// + TFlukaScoringOption* opt = new TFlukaScoringOption(option, "User Scoring", npr,outfile,what); + fUserScore->Add(opt); +} +//______________________________________________________________________________ +void TFluka::SetUserScoring(const char* option, Int_t npr, char* outfile, Float_t* what, const char* det1, const char* det2, const char* det3) +{ +// +// Adds a user scoring option to the list +// + TFlukaScoringOption* opt = new TFlukaScoringOption(option, "User Scoring", npr, outfile, what, det1, det2, det3); + fUserScore->Add(opt); } //______________________________________________________________________________ @@ -865,1136 +971,154 @@ void TFluka::InitPhysics() // // Physics initialisation with preparation of FLUKA input cards // - printf("=>InitPhysics\n"); - Int_t j, k; - Double_t fCut; - - FILE *pAliceCoreInp, *pAliceFlukaMat, *pAliceInp; - - 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 sAliceCoreInp = getenv("ALICE_ROOT"); - sAliceCoreInp +="/TFluka/input/"; - TString sAliceTmp = "flukaMat.inp"; - TString sAliceInp = GetInputFileName(); - sAliceCoreInp += GetCoreInputFileName(); - -// open files - - if ((pAliceCoreInp = fopen(sAliceCoreInp.Data(),"r")) == NULL) { - printf("\nCannot open file %s\n",sAliceCoreInp.Data()); - exit(1); - } - if ((pAliceFlukaMat = fopen(sAliceTmp.Data(),"r")) == NULL) { - printf("\nCannot open file %s\n",sAliceTmp.Data()); - exit(1); - } - if ((pAliceInp = fopen(sAliceInp.Data(),"w")) == NULL) { - printf("\nCannot open file %s\n",sAliceInp.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,pAliceCoreInp)) != NULL) { - if (strncmp(sLine,"GEOEND",6) != 0) - fprintf(pAliceInp,"%s",sLine); // copy until GEOEND card - else { - fprintf(pAliceInp,"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,pAliceFlukaMat)) != NULL) { // copy flukaMat.inp file - fprintf(pAliceInp,"%s\n",sLine); - } - - while ((fgets(sLine,255,pAliceCoreInp)) != NULL) { - if (strncmp(sLine,"START",5) != 0) - fprintf(pAliceInp,"%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(pAliceInp,"*----------------------------------------------------------------------------- \n"); - fprintf(pAliceInp,"*----- The following data are generated from SetProcess and SetCut calls ----- \n"); - fprintf(pAliceInp,"*----------------------------------------------------------------------------- \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) { - matMin = Float_t(proc->Medium()); - matMax = matMin; - global = kFALSE; - } - - // 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(pAliceInp,"*\n*Kinetic energy threshold (GeV) for e+ annihilation - resets to default=0.\n"); - fprintf(pAliceInp,"*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(pAliceInp,"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(pAliceInp,"*\n*No annihilation - no FLUKA card generated\n"); - fprintf(pAliceInp,"*Generated from call: SetProcess('ANNI',0)\n"); - } - else { - fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('ANNI',?) call.\n"); - fprintf(pAliceInp,"*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 - // gMC ->SetProcess("PAIR",1); // PAIRBREM 1. 0. 0. 3. lastmat - // EMFCUT 0. 0. -1. 3. lastmat 0. PHOT-THR - 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->Flag() == 2) && - (procp->Medium() == proc->Medium())) { - fprintf(pAliceInp,"*\n*Bremsstrahlung and pair production by muons and charged hadrons both activated\n"); - fprintf(pAliceInp,"*Generated from call: SetProcess('BREM',1) and SetProcess('PAIR',1)\n"); - fprintf(pAliceInp,"*Energy threshold set by call SetCut('BCUTM',cut) or set to 0.\n"); - fprintf(pAliceInp,"*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(pAliceInp,"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 - fCut = 0.0; - nextc.Reset(); - while ((cut = (TFlukaConfigOption*)nextc())) { - if (strncmp(cut->GetName(), "PPCUTM", 6) == 0 && - (cut->Medium() == proc->Medium())) fCut = cut->Cut(); - } - fprintf(pAliceInp,"%10.4g",fCut); - // fCut; = 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 - fCut = 0.0; - nextc.Reset(); - while ((cut = (TFlukaConfigOption*)nextc())) { - if (strncmp(cut->GetName(), "BCUTM", 5) == 0 && - (cut->Medium() == proc->Medium())) fCut = cut->Cut(); - } - fprintf(pAliceInp,"%10.4g%10.1f%10.1f\n",fCut,matMin,matMax); - // fCut = 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(pAliceInp,"*\n*Kinetic energy threshold (GeV) for e+/e- bremsstrahlung - resets to default=0.\n"); - fprintf(pAliceInp,"*Generated from call: SetProcess('BREM',1);\n"); - fCut = -1.0; - nextc.Reset(); - while ((cut = (TFlukaConfigOption*)nextc())) { - if (strncmp(cut->GetName(), "BCUTE", 5) == 0 && - (cut->Medium() == proc->Medium())) fCut = cut->Cut(); - } - //fCut = 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(pAliceInp,"EMFCUT %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fELPO-THR\n",fCut,zero,zero,matMin,matMax,one); - - // for e+ and e- - fprintf(pAliceInp,"*\n*Pair production by electrons is activated\n"); - fprintf(pAliceInp,"*Generated from call: SetProcess('PAIR',1);\n"); - fCut = -1.0; - nextc.Reset(); - while ((cut = (TFlukaConfigOption*)nextc())) { - if (strncmp(cut->GetName(), "CUTGAM", 6) == 0 && - (cut->Medium() == proc->Medium())) fCut = cut->Cut(); - } - // fCut = 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(pAliceInp,"EMFCUT %10.1f%10.1f%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",zero,zero,fCut,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(pAliceInp,"*\n*Pair production by muons and charged hadrons is activated\n"); - fprintf(pAliceInp,"*Generated from call: SetProcess('PAIR',1) or SetProcess('PAIR',2)\n"); - fprintf(pAliceInp,"*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(pAliceInp,"PAIRBREM %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,zero,zero,matMin,matMax); - - // for e+ and e- - fprintf(pAliceInp,"*\n*Pair production by electrons is activated\n"); - fprintf(pAliceInp,"*Generated from call: SetProcess('PAIR',1) or SetProcess('PAIR',2)\n"); - fCut = -1.0; - nextc.Reset(); - while ((cut = (TFlukaConfigOption*)nextc())) { - if (strncmp(cut->GetName(), "CUTGAM", 6) == 0 && - (cut->Medium() == proc->Medium())) fCut = 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) - // fCut = 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(pAliceInp,"EMFCUT %10.1f%10.1f%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",zero,zero,fCut,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 || proc->Flag() == 2) { - fprintf(pAliceInp,"*\n*Bremsstrahlung by muons and charged hadrons is activated\n"); - fprintf(pAliceInp,"*Generated from call: SetProcess('BREM',1) or SetProcess('BREM',2)\n"); - fprintf(pAliceInp,"*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 - fCut = 0.0; - nextc.Reset(); - while ((cut = (TFlukaConfigOption*)nextc())) { - if (strncmp(cut->GetName(), "BCUTM", 5) == 0 && - (cut->Medium() == proc->Medium())) fCut = cut->Cut(); - } - // fCut = 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(pAliceInp,"PAIRBREM %10.1f%10.1f%10.4g%10.1f%10.1f\n",two,zero,fCut,matMin,matMax); - - // for e+ and e- - fprintf(pAliceInp,"*\n*Kinetic energy threshold (GeV) for e+/e- bremsstrahlung - resets to default=0.\n"); - fprintf(pAliceInp,"*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"; - fprintf(pAliceInp,"EMFCUT %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fELPO-THR\n",-one,zero,zero,matMin,matMax,one); - } - else if (proc->Flag() == 0) { - fprintf(pAliceInp,"*\n*No bremsstrahlung - no FLUKA card generated\n"); - fprintf(pAliceInp,"*Generated from call: SetProcess('BREM',0)\n"); - } - else { - fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('BREM',?) call.\n"); - fprintf(pAliceInp,"*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(pAliceInp, "* \n"); - fprintf(pAliceInp, "*Cerenkov photon generation\n"); - fprintf(pAliceInp, "*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(pAliceInp, "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(pAliceInp, "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(pAliceInp, "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(pAliceInp, "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(pAliceInp, "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(pAliceInp, "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(pAliceInp,"*\n*No Cerenkov photon generation\n"); - fprintf(pAliceInp,"*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(pAliceInp,"OPT-PROD %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fCERE-OFF\n",zero,zero,zero,matMin,matMax,one); - } - else { - fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('CKOV',?) call.\n"); - fprintf(pAliceInp,"*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(pAliceInp,"*\n*Energy threshold (GeV) for Compton scattering - resets to default=0.\n"); - fprintf(pAliceInp,"*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"; - fprintf(pAliceInp,"EMFCUT %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",-one,zero,zero,matMin,matMax,one); - } - else if (proc->Flag() == 0) { - fprintf(pAliceInp,"*\n*No Compton scattering - no FLUKA card generated\n"); - fprintf(pAliceInp,"*Generated from call: SetProcess('COMP',0)\n"); - } - else { - fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('COMP',?) call.\n"); - fprintf(pAliceInp,"*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",1); // 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; - - // 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(pAliceInp,"*\n*Kinetic energy threshold (GeV) for delta ray production\n"); - fprintf(pAliceInp,"*Generated from call: SetProcess('DRAY',0) or SetProcess('DRAY',4)\n"); - fprintf(pAliceInp,"*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(pAliceInp,"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(pAliceInp,"*\n*Kinetic energy threshold (GeV) for delta ray production\n"); - fprintf(pAliceInp,"*Generated from call: SetProcess('DRAY',flag), flag=1,2,3\n"); - fprintf(pAliceInp,"*Delta ray production by muons switched on\n"); - fprintf(pAliceInp,"*Energy threshold set by call SetCut('DCUTM',cut) or set to 1.0e+6.\n"); - fCut = 1.0e+6; - nextc.Reset(); - while ((cut = (TFlukaConfigOption*)nextc())) { - if (strncmp(cut->GetName(), "DCUTM", 5) == 0 && - cut->Medium() == proc->Medium()) fCut = cut->Cut(); - } - // fCut = 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(pAliceInp,"DELTARAY %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n",fCut,zero,zero,matMin,matMax,one); - } - else { - fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('DRAY',?) call.\n"); - fprintf(pAliceInp,"*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(pAliceInp,"*\n*Hadronic interaction is ON by default in FLUKA\n"); - fprintf(pAliceInp,"*No FLUKA card generated\n"); - } - else if (proc->Flag() == 0) { - fprintf(pAliceInp,"*\n*Hadronic interaction is set OFF\n"); - fprintf(pAliceInp,"*Generated from call: SetProcess('HADR',0);\n"); - fprintf(pAliceInp,"*Switching off hadronic interactions not foreseen in FLUKA\n"); - fprintf(pAliceInp,"THRESHOL %10.1f%10.1f%10.1f%10.1e%10.1f\n",zero, zero, zero, 1.e10, zero); - } - else { - fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('HADR',?) call.\n"); - fprintf(pAliceInp,"*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() == 2) { // complete energy loss fluctuations - fprintf(pAliceInp,"*\n*Complete energy loss fluctuations do not exist in FLUKA\n"); - fprintf(pAliceInp,"*Generated from call: SetProcess('LOSS',2);\n"); - fprintf(pAliceInp,"*flag=2=complete energy loss fluctuations\n"); - fprintf(pAliceInp,"*No FLUKA card generated\n"); - } - else if (proc->Flag() == 1 || proc->Flag() == 3) { // restricted energy loss fluctuations - fprintf(pAliceInp,"*\n*Restricted energy loss fluctuations\n"); - fprintf(pAliceInp,"*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(pAliceInp,"IONFLUCT %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,one,one,matMin,matMax); - } - else if (proc->Flag() == 4) { // no energy loss fluctuations - fprintf(pAliceInp,"*\n*No energy loss fluctuations\n"); - fprintf(pAliceInp,"*\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(pAliceInp,"IONFLUCT %10.1f%10.1f%10.1f%10.1f%10.1f\n",-one,-one,one,matMin,matMax); - } - else { - fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('LOSS',?) call.\n"); - fprintf(pAliceInp,"*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(pAliceInp,"*\n*Multiple scattering is ON by default for e+e- and for hadrons/muons\n"); - fprintf(pAliceInp,"*No FLUKA card generated\n"); - } - else if (proc->Flag() == 0) { - fprintf(pAliceInp,"*\n*Multiple scattering is set OFF\n"); - fprintf(pAliceInp,"*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(pAliceInp,"MULSOPT %10.1f%10.1f%10.1f%10.1f%10.1f\n",zero,three,three,matMin,matMax); - } - else { - fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('MULS',?) call.\n"); - fprintf(pAliceInp,"*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(pAliceInp,"*\n*Muon nuclear interactions with production of secondary hadrons\n"); - fprintf(pAliceInp,"*\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(pAliceInp,"MUPHOTON %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,zero,zero,matMin,matMax); - } - else if (proc->Flag() == 2) { - fprintf(pAliceInp,"*\n*Muon nuclear interactions without production of secondary hadrons\n"); - fprintf(pAliceInp,"*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(pAliceInp,"MUPHOTON %10.1f%10.1f%10.1f%10.1f%10.1f\n",two,zero,zero,matMin,matMax); - } - else if (proc->Flag() == 0) { - fprintf(pAliceInp,"*\n*No muon nuclear interaction - no FLUKA card generated\n"); - fprintf(pAliceInp,"*Generated from call: SetProcess('MUNU',0)\n"); - } - else { - fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('MUNU',?) call.\n"); - fprintf(pAliceInp,"*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(pAliceInp,"*\n*No photonuclear interactions\n"); - fprintf(pAliceInp,"*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(pAliceInp,"PHOTONUC %10.1f%10.1f%10.1f%10.1f%10.1f\n",-one,zero,zero,matMin,matMax); - } - else if (proc->Flag() == 1) { - fprintf(pAliceInp,"*\n*Photon nuclear interactions are activated at all energies\n"); - fprintf(pAliceInp,"*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(pAliceInp,"PHOTONUC %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,zero,zero,matMin,matMax); - } - else if (proc->Flag() == 0) { - fprintf(pAliceInp,"*\n*No photofission - no FLUKA card generated\n"); - fprintf(pAliceInp,"*Generated from call: SetProcess('PFIS',0)\n"); - } - else { - fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('PFIS',?) call.\n"); - fprintf(pAliceInp,"*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 || proc->Flag() == 2) { - fprintf(pAliceInp,"*\n*Photo electric effect is activated\n"); - fprintf(pAliceInp,"*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"; - fprintf(pAliceInp,"EMFCUT %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",zero,-one,zero,matMin,matMax,one); - } - else if (proc->Flag() == 0) { - fprintf(pAliceInp,"*\n*No photo electric effect - no FLUKA card generated\n"); - fprintf(pAliceInp,"*Generated from call: SetProcess('PHOT',0)\n"); - } - else { - fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('PHOT',?) call.\n"); - fprintf(pAliceInp,"*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(pAliceInp,"*\n*Rayleigh scattering is ON by default in FLUKA\n"); - fprintf(pAliceInp,"*No FLUKA card generated\n"); - } - else if (proc->Flag() == 0) { - fprintf(pAliceInp,"*\n*Rayleigh scattering is set OFF\n"); - fprintf(pAliceInp,"*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(pAliceInp,"EMFRAY %10.1f%10.1f%10.1f%10.1f\n",-one,three,matMin,matMax); - } - else { - fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('RAYL',?) call.\n"); - fprintf(pAliceInp,"*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(pAliceInp,"*\n*Synchrotron radiation generation is NOT implemented in FLUKA\n"); - fprintf(pAliceInp,"*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(pAliceInp,"*\n*Automatic calculation of tracking medium parameters is always ON in FLUKA\n"); - fprintf(pAliceInp,"*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(pAliceInp,"*\n*Ionization energy losses calculation is activated\n"); - fprintf(pAliceInp,"*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(pAliceInp,"IONFLUCT %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,one,one,matMin,matMax); - } - else { - fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('STRA',?) call.\n"); - fprintf(pAliceInp,"*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 + } //end of while until START card + + fin: - -// 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) { - matMin = Float_t(cut->Medium()); - matMax = matMin; - global = kFALSE; - } - - // 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; - - // delta-rays by electrons - // G4 particles: "e-" - // G3 default value: 10**4 GeV - // gMC ->SetCut("DCUTE",cut); // cut for deltarays by electrons - else if (strncmp(cut->GetName(),"DCUTE",5) == 0) { - fprintf(pAliceInp,"*\n*Cut for delta rays by electrons\n"); - fprintf(pAliceInp,"*Generated from call: SetCut('DCUTE',cut);\n"); - // -cut->Cut(); - // 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 - // 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); - } - } - fprintf(pAliceInp,"DELTARAY %10.4g%10.3f%10.3f%10.1f%10.1f%10.1f\n",cut->Cut(), 100., 1.03, matMin, matMax, 1.0); - } // end of if for delta-rays by electrons +// Pass information to configuration objects + + Float_t fLastMaterial = fGeom->GetLastMaterialIndex(); + TFlukaConfigOption::SetStaticInfo(pFlukaVmcInp, 3, fLastMaterial, fGeom); + + TIter next(fUserConfig); + TFlukaConfigOption* proc; + while((proc = dynamic_cast (next()))) proc->WriteFlukaInputCards(); +// +// Process Fluka specific scoring options +// + TFlukaScoringOption::SetStaticInfo(pFlukaVmcInp, fGeom); + Float_t loginp = 49.0; + Int_t inp = 0; + Int_t nscore = fUserScore->GetEntries(); + + TFlukaScoringOption *mopo = 0x0; + TFlukaScoringOption *mopi = 0x0; - // 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(pAliceInp,"*\n*Cut for gamma\n"); - fprintf(pAliceInp,"*Generated from call: SetCut('CUTGAM',cut);\n"); - // -cut->Cut(); - // 7.0 = lower bound of the particle id-numbers to which the cut-off - fprintf(pAliceInp,"PART-THR %10.4g%10.1f\n",-cut->Cut(),7.0); - } - else if (strncmp(cut->GetName(),"CUTGAM",6) == 0 && !global) { - fprintf(pAliceInp,"*\n*Cut specific to material for gamma\n"); - fprintf(pAliceInp,"*Generated from call: SetCut('CUTGAM',cut);\n"); - // cut->Cut(); - // 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; kCut(), 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(pAliceInp,"*\n*Cut for electrons\n"); - fprintf(pAliceInp,"*Generated from call: SetCut('CUTELE',cut);\n"); - // -cut->Cut(); - // three = lower bound of the particle id-numbers to which the cut-off - // 4.0 = upper bound of the particle id-numbers to which the cut-off - // one = step length in assigning numbers - fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f%10.1f\n",-cut->Cut(),three,4.0,one); - } - else if (strncmp(cut->GetName(),"CUTELE",6) == 0 && !global) { - fprintf(pAliceInp,"*\n*Cut specific to material for electrons\n"); - fprintf(pAliceInp,"*Generated from call: SetCut('CUTELE',cut);\n"); - // -cut->Cut(); - // 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 + for (Int_t isc = 0; isc < nscore; isc++) + { + mopo = dynamic_cast (fUserScore->At(isc)); + char* fileName = mopo->GetFileName(); + Int_t size = strlen(fileName); + Float_t lun = -1.; +// +// Check if new output file has to be opened + for (Int_t isci = 0; isci < isc; isci++) { + + mopi = dynamic_cast (fUserScore->At(isci)); + if(strncmp(mopi->GetFileName(), fileName, size)==0) { + // + // No, the file already exists + lun = mopi->GetLun(); + mopo->SetLun(lun); + break; + } + } // inner loop + + if (lun == -1.) { + // Open new output file + inp++; + mopo->SetLun(loginp + inp); + mopo->WriteOpenFlukaFile(); + } + mopo->WriteFlukaInputCards(); + } - // 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(pAliceInp,"*\n*Cut for neutral hadrons\n"); - fprintf(pAliceInp,"*Generated from call: SetCut('CUTNEU',cut);\n"); - - // 8.0 = Neutron - // 9.0 = Antineutron - fprintf(pAliceInp,"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(pAliceInp,"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(pAliceInp,"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(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),22.0,25.0); - - // 32.0 = Antisigma zero - // 32.0 = Antisigma zero - fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),32.0,32.0); - - // 34.0 = Xi zero - // 35.0 = AntiXi zero - fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),34.0,35.0); - - // 47.0 = D zero - // 48.0 = AntiD zero - fprintf(pAliceInp,"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(pAliceInp,"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(pAliceInp,"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(pAliceInp,"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(pAliceInp,"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(pAliceInp,"*\n*Cut for charged hadrons\n"); - fprintf(pAliceInp,"*Generated from call: SetCut('CUTHAD',cut);\n"); - - // 1.0 = Proton - // 2.0 = Antiproton - fprintf(pAliceInp,"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(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),13.0,16.0); - - // 20.0 = Negative Sigma - // 21.0 = Positive Sigma - fprintf(pAliceInp,"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(pAliceInp,"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(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),36.0,39.0); - - // 45.0 = D plus - // 46.0 = D minus - fprintf(pAliceInp,"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(pAliceInp,"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(pAliceInp,"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(pAliceInp,"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(pAliceInp,"*\n*Cut for muons\n"); - fprintf(pAliceInp,"*Generated from call: SetCut('CUTMUO',cut);\n"); - // 10.0 = Muon+ - // 11.0 = Muon- - fprintf(pAliceInp,"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(pAliceInp,"*\n*Time of flight cuts in seconds\n"); - fprintf(pAliceInp,"*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(pAliceInp,"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 - // Add START and STOP card - fprintf(pAliceInp,"START %10.1f\n",fEventsPerRun); - fprintf(pAliceInp,"STOP \n"); + fprintf(pFlukaVmcInp,"START %10.1f\n",fEventsPerRun); + fprintf(pFlukaVmcInp,"STOP \n"); // Close files - - fclose(pAliceCoreInp); - fclose(pAliceFlukaMat); - fclose(pAliceInp); - + 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 //______________________________________________________________________________ -void TFluka::SetMaxStep(Double_t) +void TFluka::SetMaxStep(Double_t step) { -// SetMaxStep is dummy procedure in TFluka ! - if (fVerbosityLevel >=3) - cout << "SetMaxStep is dummy procedure in TFluka !" << endl; +// Set the maximum step size + if (step > 1.e4) return; + + Int_t mreg, latt; + fGeom->GetCurrentRegion(mreg, latt); + STEPSZ.stepmx[mreg - 1] = step; +} + + +Double_t TFluka::MaxStep() const +{ +// Return the maximum for current medium + Int_t mreg, latt; + fGeom->GetCurrentRegion(mreg, latt); + return (STEPSZ.stepmx[mreg - 1]); } //______________________________________________________________________________ @@ -2026,13 +1150,13 @@ void TFluka::TrackPosition(TLorentzVector& position) const // TRACKR.ytrack = y-position of the last point // TRACKR.ztrack = z-position of the last point Int_t caller = GetCaller(); - if (caller == 3 || caller == 6 || caller == 11 || caller == 12) { //bxdraw,endraw,usdraw + if (caller == 3 || caller == 6 || caller == 11 || caller == 12 || caller == 50) { //bxdraw,endraw,usdraw,ckov position.SetX(GetXsco()); position.SetY(GetYsco()); position.SetZ(GetZsco()); position.SetT(TRACKR.atrack); } - else if (caller == 4) { // mgdraw + else if (caller == 4) { // mgdraw,mgdraw resuming position.SetX(TRACKR.xtrack[TRACKR.ntrack]); position.SetY(TRACKR.ytrack[TRACKR.ntrack]); position.SetZ(TRACKR.ztrack[TRACKR.ntrack]); @@ -2043,6 +1167,11 @@ void TFluka::TrackPosition(TLorentzVector& position) const position.SetY(TRACKR.ytrack[TRACKR.ntrack]); position.SetZ(TRACKR.ztrack[TRACKR.ntrack]); position.SetT(0); + } else if (caller == 40) { // mgdraw resuming transport + position.SetX(TRACKR.spausr[0]); + position.SetY(TRACKR.spausr[1]); + position.SetZ(TRACKR.spausr[2]); + position.SetT(TRACKR.spausr[3]); } else Warning("TrackPosition","position not available"); @@ -2058,16 +1187,21 @@ void TFluka::TrackPosition(Double_t& x, Double_t& y, Double_t& z) const // TRACKR.ytrack = y-position of the last point // TRACKR.ztrack = z-position of the last point Int_t caller = GetCaller(); - if (caller == 3 || caller == 6 || caller == 11 || caller == 12) { //bxdraw,endraw,usdraw + if (caller == 3 || caller == 6 || caller == 11 || caller == 12 || caller == 50) { //bxdraw,endraw,usdraw,ckov x = GetXsco(); y = GetYsco(); z = GetZsco(); } - else if (caller == 4 || caller == 5) { // mgdraw, sodraw + else if (caller == 4 || caller == 5) { // mgdraw, sodraw, mgdraw resuming x = TRACKR.xtrack[TRACKR.ntrack]; y = TRACKR.ytrack[TRACKR.ntrack]; z = TRACKR.ztrack[TRACKR.ntrack]; } + else if (caller == 40) { // mgdraw resuming transport + x = TRACKR.spausr[0]; + y = TRACKR.spausr[1]; + z = TRACKR.spausr[2]; + } else Warning("TrackPosition","position not available"); } @@ -2084,7 +1218,7 @@ void TFluka::TrackMomentum(TLorentzVector& momentum) const // TRACKR.jtrack = identity number of the particle // PAPROP.am[TRACKR.jtrack] = particle mass in gev Int_t caller = GetCaller(); - if (caller != 2) { // not eedraw + if (caller != 2 && caller != 40) { // not eedraw or mgdraw resuming if (TRACKR.ptrack >= 0) { momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck); momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck); @@ -2100,6 +1234,12 @@ void TFluka::TrackMomentum(TLorentzVector& momentum) const momentum.SetE(TRACKR.etrack); return; } + } else if (caller == 40) { // mgdraw resuming transport + momentum.SetPx(TRACKR.spausr[4]); + momentum.SetPy(TRACKR.spausr[5]); + momentum.SetPz(TRACKR.spausr[6]); + momentum.SetE (TRACKR.spausr[7]); + return; } else Warning("TrackMomentum","momentum not available"); @@ -2117,7 +1257,7 @@ void TFluka::TrackMomentum(Double_t& px, Double_t& py, Double_t& pz, Double_t& e // TRACKR.jtrack = identity number of the particle // PAPROP.am[TRACKR.jtrack] = particle mass in gev Int_t caller = GetCaller(); - if (caller != 2) { // not eedraw + if (caller != 2 && caller != 40) { // not eedraw and not mgdraw resuming if (TRACKR.ptrack >= 0) { px = TRACKR.ptrack*TRACKR.cxtrck; py = TRACKR.ptrack*TRACKR.cytrck; @@ -2133,6 +1273,12 @@ void TFluka::TrackMomentum(Double_t& px, Double_t& py, Double_t& pz, Double_t& e e = TRACKR.etrack; return; } + } else if (caller == 40) { // mgdraw resuming transport + px = TRACKR.spausr[4]; + py = TRACKR.spausr[5]; + pz = TRACKR.spausr[6]; + e = TRACKR.spausr[7]; + return; } else Warning("TrackMomentum","momentum not available"); @@ -2144,12 +1290,14 @@ Double_t TFluka::TrackStep() const // Return the length in centimeters of the current step // TRACKR.ctrack = total curved path Int_t caller = GetCaller(); - if (caller == 11 || caller==12 || caller == 3 || caller == 6) //bxdraw,endraw,usdraw + if (caller == 11 || caller==12 || caller == 3 || caller == 6 || caller == 50 || caller == 40) //bxdraw,endraw,usdraw, ckov return 0.0; else if (caller == 4) //mgdraw return TRACKR.ctrack; - else - return -1.0; + else { + Warning("TrackStep", "track step not available"); + return 0.0; + } } //______________________________________________________________________________ @@ -2157,10 +1305,14 @@ Double_t TFluka::TrackLength() const { // TRACKR.cmtrck = cumulative curved path since particle birth Int_t caller = GetCaller(); - if (caller == 11 || caller==12 || caller == 3 || caller == 4 || caller == 6) //bxdraw,endraw,mgdraw,usdraw + if (caller == 11 || caller==12 || caller == 3 || caller == 4 || caller == 6 || caller == 50) //bxdraw,endraw,mgdraw,usdraw,ckov return TRACKR.cmtrck; - else - return -1.0; + else if (caller == 40) // mgdraw resuming transport + return TRACKR.spausr[8]; + else { + Warning("TrackLength", "track length not available"); + return 0.0; + } } //______________________________________________________________________________ @@ -2169,10 +1321,14 @@ Double_t TFluka::TrackTime() const // Return the current time of flight of the track being transported // TRACKR.atrack = age of the particle Int_t caller = GetCaller(); - if (caller == 11 || caller==12 || caller == 3 || caller == 4 || caller == 6) //bxdraw,endraw,mgdraw,usdraw + if (caller == 11 || caller==12 || caller == 3 || caller == 4 || caller == 6 || caller == 50) //bxdraw,endraw,mgdraw,usdraw,ckov return TRACKR.atrack; - else - return -1; + else if (caller == 40) + return TRACKR.spausr[3]; + else { + Warning("TrackTime", "track time not available"); + return 0.0; + } } //______________________________________________________________________________ @@ -2186,19 +1342,21 @@ Double_t TFluka::Edep() const // -->no energy loss along the track // if TRACKR.ntrack > 0, TRACKR.mtrack > 0: // -->energy loss distributed along the track -// TRACKR.dtrack = energy deposition of the jth deposition even +// TRACKR.dtrack = energy deposition of the jth deposition event // If coming from bxdraw we have 2 steps of 0 length and 0 edep + // If coming from usdraw we just signal particle production - no edep + // If just first time after resuming, no edep for the primary Int_t caller = GetCaller(); - if (caller == 11 || caller==12) return 0.0; + if (caller == 11 || caller==12 || caller==6 || caller == 40) return 0.0; Double_t sum = 0; for ( Int_t j=0;j= 0 && isec < FINUC.np) { - particleId = PDGFromId(FINUC.kpart[isec]); - position.SetX(fXsco); - position.SetY(fYsco); - position.SetZ(fZsco); - position.SetT(TRACKR.atrack); - momentum.SetPx(FINUC.plr[isec]*FINUC.cxr[isec]); - momentum.SetPy(FINUC.plr[isec]*FINUC.cyr[isec]); - momentum.SetPz(FINUC.plr[isec]*FINUC.czr[isec]); - momentum.SetE(FINUC.tki[isec] + PAPROP.am[FINUC.kpart[isec]+6]); - } - else if (isec >= FINUC.np && isec < FINUC.np + FHEAVY.npheav) { - Int_t jsec = isec - FINUC.np; - particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!! - position.SetX(fXsco); - position.SetY(fYsco); - position.SetZ(fZsco); - position.SetT(TRACKR.atrack); - momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]); - momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]); - momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]); - if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6) - momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]); - else if (FHEAVY.tkheav[jsec] > 6) - momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!! + Int_t caller = GetCaller(); + if (caller == 6) { // valid only after usdraw + if (FINUC.np > 0) { + // Hadronic interaction + if (isec >= 0 && isec < FINUC.np) { + particleId = PDGFromId(FINUC.kpart[isec]); + position.SetX(fXsco); + position.SetY(fYsco); + position.SetZ(fZsco); + position.SetT(TRACKR.atrack); + momentum.SetPx(FINUC.plr[isec]*FINUC.cxr[isec]); + momentum.SetPy(FINUC.plr[isec]*FINUC.cyr[isec]); + momentum.SetPz(FINUC.plr[isec]*FINUC.czr[isec]); + momentum.SetE(FINUC.tki[isec] + PAPROP.am[FINUC.kpart[isec]+6]); + } + else if (isec >= FINUC.np && isec < FINUC.np + FHEAVY.npheav) { + Int_t jsec = isec - FINUC.np; + particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!! + position.SetX(fXsco); + position.SetY(fYsco); + position.SetZ(fZsco); + position.SetT(TRACKR.atrack); + momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]); + momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]); + momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]); + if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6) + momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]); + else if (FHEAVY.tkheav[jsec] > 6) + momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!! + } + else + Warning("GetSecondary","isec out of range"); + } + } else if (caller == 50) { + Int_t index = OPPHST.lstopp - isec; + position.SetX(OPPHST.xoptph[index]); + position.SetY(OPPHST.yoptph[index]); + position.SetZ(OPPHST.zoptph[index]); + position.SetT(OPPHST.agopph[index]); + Double_t p = OPPHST.poptph[index]; + + momentum.SetPx(p * OPPHST.txopph[index]); + momentum.SetPy(p * OPPHST.tyopph[index]); + momentum.SetPz(p * OPPHST.tzopph[index]); + momentum.SetE(p); } else - Warning("GetSecondary","isec out of range"); - } - else - Warning("GetSecondary","no secondaries available"); + Warning("GetSecondary","no secondaries available"); + } // end of GetSecondary + //______________________________________________________________________________ TMCProcess TFluka::ProdProcess(Int_t) const @@ -2449,26 +1627,66 @@ TMCProcess TFluka::ProdProcess(Int_t) const Int_t mugamma = (TRACKR.jtrack == 7 || TRACKR.jtrack == 10 || TRACKR.jtrack == 11); - if (fIcode == 102) return kPDecay; + if (fIcode == 102) return kPDecay; else if (fIcode == 104 || fIcode == 217) return kPPair; - else if (fIcode == 219) return kPCompton; - else if (fIcode == 221) return kPPhotoelectric; + else if (fIcode == 219) return kPCompton; + else if (fIcode == 221) return kPPhotoelectric; else if (fIcode == 105 || fIcode == 208) return kPBrem; else if (fIcode == 103 || fIcode == 400) return kPDeltaRay; else if (fIcode == 210 || fIcode == 212) return kPDeltaRay; else if (fIcode == 214 || fIcode == 215) return kPAnnihilation; - else if (fIcode == 101) return kPHadronic; + else if (fIcode == 101) return kPHadronic; else if (fIcode == 101) { - if (!mugamma) return kPHadronic; - else if (TRACKR.jtrack == 7) return kPPhotoFission; - else return kPMuonNuclear; + if (!mugamma) return kPHadronic; + else if (TRACKR.jtrack == 7) return kPPhotoFission; + else return kPMuonNuclear; } - else if (fIcode == 225) return kPRayleigh; + else if (fIcode == 225) return kPRayleigh; // Fluka codes 100, 300 and 400 still to be investigasted - else return kPNoProcess; + else return kPNoProcess; } +Int_t TFluka::StepProcesses(TArrayI &proc) const +{ + // + // Return processes active in the current step + // + proc.Set(1); + TMCProcess iproc; + switch (fIcode) { + case 15: + case 24: + case 33: + case 41: + case 52: + iproc = kPTOFlimit; + break; + case 12: + case 14: + case 21: + case 22: + case 23: + case 31: + case 32: + case 40: + case 51: + iproc = kPStop; + break; + case 50: + iproc = kPLightAbsorption; + break; + case 59: + iproc = kPLightRefraction; + case 20: + iproc = kPPhotoelectric; + break; + default: + iproc = ProdProcess(0); + } + proc[0] = iproc; + return 1; +} //______________________________________________________________________________ Int_t TFluka::VolId2Mate(Int_t id) const { @@ -2495,7 +1713,12 @@ Int_t TFluka::VolId(const Text_t* volName) const // Time consuming. (Only used during set-up) // Could be replaced by hash-table // - return fMCGeo->VolId(volName); + char sname[20]; + Int_t len; + strncpy(sname, volName, len = strlen(volName)); + sname[len] = 0; + while (sname[len - 1] == ' ') sname[--len] = 0; + return fMCGeo->VolId(sname); } //______________________________________________________________________________ @@ -2550,15 +1773,23 @@ const char* TFluka::CurrentVolOffName(Int_t off) const } //______________________________________________________________________________ -Int_t TFluka::CurrentMaterial(Float_t & /*a*/, Float_t & /*z*/, - Float_t & /*dens*/, Float_t & /*radl*/, Float_t & /*absl*/) const +Int_t TFluka::CurrentMaterial(Float_t & a, Float_t & z, + Float_t & dens, Float_t & radl, Float_t & absl) const { // -// Return the current medium number ??? what about material properties +// Return the current medium number and material properties // Int_t copy; Int_t id = TFluka::CurrentVolID(copy); Int_t med = TFluka::VolId2Mate(id); + TGeoVolume* vol = gGeoManager->GetCurrentVolume(); + TGeoMaterial* mat = vol->GetMaterial(); + a = mat->GetA(); + z = mat->GetZ(); + dens = mat->GetDensity(); + radl = mat->GetRadLen(); + absl = mat->GetIntLen(); + return med; } @@ -2646,7 +1877,66 @@ void TFluka::SetMreg(Int_t l) } + + +TString TFluka::ParticleName(Int_t pdg) const +{ + // Return particle name for particle with pdg code pdg. + Int_t ifluka = IdFromPDG(pdg); + return TString((CHPPRP.btype[ifluka+6]), 8); +} + + +Double_t TFluka::ParticleMass(Int_t pdg) const +{ + // Return particle mass for particle with pdg code pdg. + Int_t ifluka = IdFromPDG(pdg); + return (PAPROP.am[ifluka+6]); +} + +Double_t TFluka::ParticleCharge(Int_t pdg) const +{ + // Return particle charge for particle with pdg code pdg. + Int_t ifluka = IdFromPDG(pdg); + return Double_t(PAPROP.ichrge[ifluka+6]); +} + +Double_t TFluka::ParticleLifeTime(Int_t pdg) const +{ + // Return particle lifetime for particle with pdg code pdg. + Int_t ifluka = IdFromPDG(pdg); + return (PAPROP.thalf[ifluka+6]); +} + +void TFluka::Gfpart(Int_t pdg, char* name, Int_t& type, Float_t& mass, Float_t& charge, Float_t& tlife) +{ + // Retrieve particle properties for particle with pdg code pdg. + + strcpy(name, ParticleName(pdg).Data()); + type = ParticleMCType(pdg); + mass = ParticleMass(pdg); + charge = ParticleCharge(pdg); + tlife = ParticleLifeTime(pdg); +} + +void TFluka::PrintHeader() +{ + // + // Print a header + printf("\n"); + printf("\n"); + printf("------------------------------------------------------------------------------\n"); + printf("- You are using the TFluka Virtual Monte Carlo Interface to FLUKA. -\n"); + printf("- Please see the file fluka.out for FLUKA output and licensing information. -\n"); + printf("------------------------------------------------------------------------------\n"); + printf("\n"); + printf("\n"); +} + + + #define pushcerenkovphoton pushcerenkovphoton_ +#define usersteppingckv usersteppingckv_ extern "C" { @@ -2667,6 +1957,21 @@ extern "C" { polx, poly, polz, kPCerenkov, ntr, wgt, 0); } -} + void usersteppingckv(Int_t & nphot, Int_t & mreg, Double_t & x, Double_t & y, Double_t & z) + { + // + // Calls stepping in order to signal cerenkov production + // + TFluka *fluka = (TFluka*)gMC; + fluka->SetMreg(mreg); + fluka->SetXsco(x); + fluka->SetYsco(y); + fluka->SetZsco(z); + fluka->SetNCerenkov(nphot); + fluka->SetCaller(50); + printf("userstepping ckv: %10d %10d %13.3f %13.3f %13.2f %s\n", nphot, mreg, x, y, z, fluka->CurrentVolName()); + (TVirtualMCApplication::Instance())->Stepping(); + } +}