/* $Id$ */
+//
+// Realisation of the TVirtualMC interface for the FLUKA code
+// (See official web side http://www.fluka.org/).
+//
+// This implementation makes use of the TGeo geometry modeller.
+// User configuration is via automatic generation of FLUKA input cards.
+//
+// Authors:
+// A. Fasso
+// E. Futo
+// A. Gheata
+// A. Morsch
+//
+
#include <Riostream.h>
-#include "AliModule.h"
-#include "AliRun.h"
+//#include "AliModule.h"
+//#include "AliRun.h"
#include "TClonesArray.h"
#include "TFlukaGeo.h"
#include "TCallf77.h" //For the fortran calls
TFluka::TFluka()
:TVirtualMC(),
fVerbosityLevel(0),
- sInputFileName("")
+ fInputFileName("")
{
//
// Default constructor
fCurrentFlukaRegion = -1;
fGeom = 0;
fMaterials = 0;
+ fDummyBoundary = 0;
+ fFieldFlag = 1;
}
//______________________________________________________________________________
TFluka::TFluka(const char *title, Int_t verbosity, Bool_t isRootGeometrySupported)
:TVirtualMC("TFluka",title, isRootGeometrySupported),
fVerbosityLevel(verbosity),
- sInputFileName(""),
+ fInputFileName(""),
fTrackIsEntering(0),
fTrackIsExiting(0),
fTrackIsNew(0)
fNVolumes = 0;
fCurrentFlukaRegion = -1;
+ fDummyBoundary = 0;
+ fFieldFlag = 1;
fGeom = new TFlukaMCGeometry("geom", "ALICE geometry");
+ if (verbosity > 2) fGeom->SetDebugMode(kTRUE);
fMaterials = 0;
}
//______________________________________________________________________________
TFluka::~TFluka() {
- if (fVerbosityLevel >=3)
- cout << "==> TFluka::~TFluka() destructor called." << endl;
-
+// Destructor
delete fGeom;
-
if (fVerbosityLevel >=3)
cout << "<== TFluka::~TFluka() destructor called." << endl;
}
// TFluka control methods
//______________________________________________________________________________
void TFluka::Init() {
-
- if (fVerbosityLevel >=3)
- cout << "==> TFluka::Init() called." << endl;
+//
+// Geometry initialisation
+//
+ if (fVerbosityLevel >=3) cout << "==> TFluka::Init() called." << endl;
if (!gGeoManager) new TGeoManager("geom", "FLUKA geometry");
fApplication->ConstructGeometry();
gGeoManager->CloseGeometry("di");
gGeoManager->DefaultColors(); // to be removed
fNVolumes = fGeom->NofVolumes();
- printf("== Number of volumes: %i\n ==", fNVolumes);
fGeom->CreateFlukaMatFile("flukaMat.inp");
- cout << "\t* InitPhysics() - Prepare input file to be called" << endl;
+ 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
}
//
// Build-up table with region to medium correspondance
//
- if (fVerbosityLevel >=3)
+ if (fVerbosityLevel >=3) {
cout << "==> TFluka::FinishGeometry() called." << endl;
-
- printf("----FinishGeometry - nothing to do with TGeo\n");
-
- if (fVerbosityLevel >=3)
+ printf("----FinishGeometry - nothing to do with TGeo\n");
cout << "<== TFluka::FinishGeometry() called." << endl;
+ }
}
//______________________________________________________________________________
void TFluka::BuildPhysics() {
+//
+// Prepare FLUKA input files and call FLUKA physics initialisation
+//
+
if (fVerbosityLevel >=3)
cout << "==> TFluka::BuildPhysics() called." << endl;
- InitPhysics(); // prepare input file with the current physics settings
+// Prepare input file with the current physics settings
+ InitPhysics();
cout << "\t* InitPhysics() - Prepare input file was called" << endl;
if (fVerbosityLevel >=2)
GLOBAL.lfdrtr = true;
if (fVerbosityLevel >=2)
- cout << "\t* Opening file " << sInputFileName << endl;
- const char* fname = sInputFileName;
+ cout << "\t* Opening file " << fInputFileName << endl;
+ const char* fname = fInputFileName;
fluka_openinp(lunin, PASSCHARA(fname));
if (fVerbosityLevel >=2)
flukam(1);
if (fVerbosityLevel >=2)
- cout << "\t* Closing file " << sInputFileName << endl;
+ cout << "\t* Closing file " << fInputFileName << endl;
fluka_closeinp(lunin);
FinishGeometry();
//______________________________________________________________________________
void TFluka::ProcessEvent() {
+//
+// Process one event
+//
if (fVerbosityLevel >=3)
cout << "==> TFluka::ProcessEvent() called." << endl;
fApplication->GeneratePrimaries();
}
//______________________________________________________________________________
-void TFluka::ProcessRun(Int_t nevent) {
+#if ROOT_VERSION_CODE >= 262150
+Bool_t TFluka::ProcessRun(Int_t nevent) {
+#else
+void TFluka::ProcessRun(Int_t nevent) {
+#endif
+//
+// Run steering
+//
+
if (fVerbosityLevel >=3)
cout << "==> TFluka::ProcessRun(" << nevent << ") called."
<< endl;
cout << "\t* GLOBAL.fdrtr = " << (GLOBAL.lfdrtr?'T':'F') << endl;
cout << "\t* Calling flukam again..." << endl;
}
- fApplication->InitGeometry();
- fApplication->BeginEvent();
- ProcessEvent();
+
+ Int_t todo = TMath::Abs(nevent);
+
+ for (Int_t ev = 0; ev < todo; ev++) {
+ fApplication->InitGeometry();
+ fApplication->BeginEvent();
+ ProcessEvent();
+ }
+
fApplication->FinishEvent();
if (fVerbosityLevel >=3)
cout << "<== TFluka::ProcessRun(" << nevent << ") called."
<< endl;
+ #if ROOT_VERSION_CODE >= 262150
+ return kTRUE;
+ #endif
}
//_____________________________________________________________________________
// Is it needed with TGeo ??? - to clear-up
//
- printf("Gstpar called with %6d %5s %12.4e %6d\n", itmed, param, parval, fGeom->GetFlukaMaterial(itmed));
+ if (fVerbosityLevel >=3) printf("Gstpar called with %6d %5s %12.4e %6d\n", itmed, param, parval, fGeom->GetFlukaMaterial(itmed));
Bool_t process = kFALSE;
if (strncmp(param, "DCAY", 4) == 0 ||
} else {
SetCut(param, parval, fGeom->GetFlukaMaterial(itmed));
}
-
-
-
}
// functions from GGEOM
return 50000050;
}
// Error id
- if (id == 0) {
+ if (id == 0 || id < -6 || id > 250) {
if (fVerbosityLevel >= 1)
printf("PDGFromId: Error id = 0\n");
return -1;
//
// set methods
//
-
-void TFluka::SetProcess(const char* flagName, Int_t flagValue, Int_t imed)
+ void TFluka::SetProcess(const char* flagName, Int_t flagValue, Int_t imat)
{
+// Set process user flag for material imat
+//
strcpy(&fProcessFlag[fNbOfProc][0],flagName);
fProcessValue[fNbOfProc] = flagValue;
- fProcessMedium[fNbOfProc] = imed;
+ fProcessMaterial[fNbOfProc] = imat;
fNbOfProc++;
}
//______________________________________________________________________________
-void TFluka::SetProcess(const char* flagName, Int_t flagValue)
+
+#if ROOT_VERSION_CODE >= 262151
+ Bool_t TFluka::SetProcess(const char* flagName, Int_t flagValue)
+#else
+ void TFluka::SetProcess(const char* flagName, Int_t flagValue)
+#endif
{
+// Set process user flag
+//
+
Int_t i;
if (fNbOfProc < 100) {
for (i=0; i<fNbOfProc; i++) {
if (strcmp(&fProcessFlag[i][0],flagName) == 0) {
fProcessValue[fNbOfProc] = flagValue;
- return;
+ fProcessMaterial[fNbOfProc] = -1;
+#if ROOT_VERSION_CODE >= 262151
+ return kFALSE;
+#else
+ return
+#endif
}
}
strcpy(&fProcessFlag[fNbOfProc][0],flagName);
- fProcessValue[fNbOfProc++] = flagValue;
+ fProcessMaterial[fNbOfProc] = -1;
+ fProcessValue[fNbOfProc++] = flagValue;
+
}
else
cout << "Nb of SetProcess calls exceeds 100 - ignored" << endl;
+#if ROOT_VERSION_CODE >= 262151
+ return kTRUE;
+#endif
}
//______________________________________________________________________________
-void TFluka::SetCut(const char* cutName, Double_t cutValue, Int_t imed)
+ void TFluka::SetCut(const char* cutName, Double_t cutValue, Int_t imed)
{
+// Set user cut value for material imed
+//
strcpy(&fCutFlag[fNbOfCut][0],cutName);
fCutValue[fNbOfCut] = cutValue;
- fCutMedium[fNbOfCut] = imed;
+ fCutMaterial[fNbOfCut] = imed;
fNbOfCut++;
}
//______________________________________________________________________________
-void TFluka::SetCut(const char* cutName, Double_t cutValue)
+#if ROOT_VERSION_CODE >= 262151
+ Bool_t TFluka::SetCut(const char* cutName, Double_t cutValue)
+#else
+ void TFluka::SetCut(const char* cutName, Double_t cutValue)
+#endif
{
+// Set user cut value
+//
Int_t i;
if (fNbOfCut < 100) {
for (i=0; i<fNbOfCut; i++) {
if (strcmp(&fCutFlag[i][0],cutName) == 0) {
fCutValue[fNbOfCut] = cutValue;
+#if ROOT_VERSION_CODE >= 262151
+ return kFALSE;
+#else
return;
+
+#endif
}
}
strcpy(&fCutFlag[fNbOfCut][0],cutName);
- fCutMedium[fNbOfCut] = -1;
+ fCutMaterial[fNbOfCut] = -1;
fCutValue[fNbOfCut++] = cutValue;
}
else
cout << "Nb of SetCut calls exceeds 100 - ignored" << endl;
+#if ROOT_VERSION_CODE >= 262151
+ return kTRUE;
+#endif
}
//______________________________________________________________________________
//______________________________________________________________________________
void TFluka::InitPhysics()
{
+//
+// Physics initialisation with preparation of FLUKA input cards
+//
+ printf("=>InitPhysics\n");
Int_t i, j, k;
Double_t fCut;
- FILE *pAliceCoreInp, *pAliceFlukaMat, *pAliceInp, *pGaliceCuts;
+ FILE *pAliceCoreInp, *pAliceFlukaMat, *pAliceInp;
Double_t zero = 0.0;
Double_t one = 1.0;
Double_t three = 3.0;
Float_t fLastMaterial = fGeom->GetLastMaterialIndex();
- printf(" last FLUKA material is %g\n", fLastMaterial);
+ if (fVerbosityLevel >= 3) printf(" last FLUKA material is %g\n", fLastMaterial);
// Prepare Cerenkov
TList *matList = gGeoManager->GetListOfMaterials();
Int_t nmaterial = matList->GetSize();
fMaterials = new Int_t[nmaterial];
-
-
// construct file names
TString sAliceCoreInp = getenv("ALICE_ROOT");
- TString sGaliceCuts = sAliceCoreInp;
sAliceCoreInp +="/TFluka/input/";
TString sAliceTmp = "flukaMat.inp";
TString sAliceInp = GetInputFileName();
sAliceCoreInp += GetCoreInputFileName();
- sGaliceCuts +="/data/galice.cuts";
// open files
exit(1);
}
- if ((pGaliceCuts = fopen(sGaliceCuts.Data(),"r")) == NULL) {
- printf("\nCannot open file %s\n",sGaliceCuts.Data());
- exit(1);
- }
-
// copy core input file
Char_t sLine[255];
Float_t fEventsPerRun;
Float_t matMin = three;
Float_t matMax = fLastMaterial;
Bool_t global = kTRUE;
- if (fProcessMedium[i] != -1) {
- matMin = Float_t(fProcessMedium[i]);
+ if (fProcessMaterial[i] != -1) {
+ matMin = Float_t(fProcessMaterial[i]);
matMax = matMin;
global = kFALSE;
}
for (j=0; j<fNbOfProc; j++) {
if ((strncmp(&fProcessFlag[j][0],"BREM",4) == 0) &&
(fProcessValue[j] == 1 || fProcessValue[j] == 2) &&
- (fProcessMedium[j] == fProcessMedium[i])) {
+ (fProcessMaterial[j] == fProcessMaterial[i])) {
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");
fCut = 0.0;
for (k=0; k<fNbOfCut; k++) {
if (strncmp(&fCutFlag[k][0],"PPCUTM",6) == 0 &&
- (fCutMedium[k] == fProcessMedium[i])) fCut = fCutValue[k];
+ (fCutMaterial[k] == fProcessMaterial[i])) fCut = fCutValue[k];
}
fprintf(pAliceInp,"%10.4g",fCut);
// fCut; = e+, e- kinetic energy threshold (in GeV) for explicit pair production.
fCut = 0.0;
for (k=0; k<fNbOfCut; k++) {
if (strncmp(&fCutFlag[k][0],"BCUTM",5) == 0 &&
- (fCutMedium[k] == fProcessMedium[i])) fCut = fCutValue[k];
+ (fCutMaterial[k] == fProcessMaterial[i])) fCut = fCutValue[k];
}
fprintf(pAliceInp,"%10.4g%10.1f%10.1f\n",fCut,matMin,matMax);
// fCut = photon energy threshold (GeV) for explicit bremsstrahlung production
fCut = -1.0;
for (k=0; k<fNbOfCut; k++) {
if (strncmp(&fCutFlag[k][0],"BCUTE",5) == 0 &&
- (fCutMedium[k] == fProcessMedium[i])) fCut = fCutValue[k];
+ (fCutMaterial[k] == fProcessMaterial[i])) fCut = fCutValue[k];
}
//fCut = kinetic energy threshold (GeV) for e+/e- bremsstrahlung (resets to default=0)
// zero = not used
fCut = -1.0;
for (k=0; k<fNbOfCut; k++) {
if (strncmp(&fCutFlag[k][0],"CUTGAM",6) == 0 &&
- (fCutMedium[k] == fProcessMedium[i])) fCut = fCutValue[k];
+ (fCutMaterial[k] == fProcessMaterial[i])) fCut = fCutValue[k];
}
// 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
fCut = -1.0;
for (j=0; j<fNbOfCut; j++) {
if (strncmp(&fCutFlag[j][0],"CUTGAM",6) == 0 &&
- (fCutMedium[j] == fProcessMedium[i])) fCut = fCutValue[j];
+ (fCutMaterial[j] == fProcessMaterial[i])) fCut = fCutValue[j];
}
// zero = energy threshold (GeV) for Compton scattering (= 0.0 : ignored)
// zero = energy threshold (GeV) for Photoelectric (= 0.0 : ignored)
for (j = 0; j < fNbOfProc; j++) {
if ((strncmp(&fProcessFlag[j][0],"PAIR",4) == 0) &&
fProcessValue[j] == 1 &&
- (fProcessMedium[j] == fProcessMedium[i])) goto NOBREM;
+ (fProcessMaterial[j] == fProcessMaterial[i])) goto NOBREM;
}
if (fProcessValue[i] == 1 || fProcessValue[i] == 2) {
fprintf(pAliceInp,"*\n*Bremsstrahlung by muons and charged hadrons is activated\n");
fCut = 0.0;
for (j=0; j<fNbOfCut; j++) {
if (strncmp(&fCutFlag[j][0],"BCUTM",5) == 0 &&
- (fCutMedium[j] == fProcessMedium[i])) fCut = fCutValue[j];
+ (fCutMaterial[j] == fProcessMaterial[i])) fCut = fCutValue[j];
}
// fCut = photon energy threshold (GeV) for explicit bremsstrahlung production
// matMin = lower bound of the material indices in which the respective thresholds apply
TGeoMaterial* material = dynamic_cast<TGeoMaterial*> (matList->At(im));
Int_t idmat = material->GetIndex();
- if (!global && idmat != fProcessMedium[i]) continue;
+ if (!global && idmat != fProcessMaterial[i]) continue;
fMaterials[idmat] = im;
// Skip media with no Cerenkov properties
fCut = 1.0e+6;
for (j = 0; j < fNbOfCut; j++) {
if (strncmp(&fCutFlag[j][0],"DCUTM",5) == 0 &&
- fCutMedium[j] == fProcessMedium[i]) fCut = fCutValue[j];
+ fCutMaterial[j] == fProcessMaterial[i]) fCut = fCutValue[j];
}
// fCut = kinetic energy threshold (GeV) for delta ray production (discrete energy transfer)
// zero = ignored
Float_t matMin = three;
Float_t matMax = fLastMaterial;
Bool_t global = kTRUE;
- if (fCutMedium[i] != -1) {
- matMin = Float_t(fCutMedium[i]);
- matMax = matMin;
- global = kFALSE;
+ if (fCutMaterial[i] != -1) {
+ matMin = Float_t(fCutMaterial[i]);
+ matMax = matMin;
+ global = kFALSE;
}
- // cuts handled in SetProcess calls
+
+ // cuts handled in SetProcess calls
if (strncmp(&fCutFlag[i][0],"BCUTM",5) == 0) continue;
else if (strncmp(&fCutFlag[i][0],"BCUTE",5) == 0) continue;
else if (strncmp(&fCutFlag[i][0],"DCUTM",5) == 0) continue;
else if (strncmp(&fCutFlag[i][0],"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(&fCutFlag[i][0],"DCUTE",5) == 0) {
+ fprintf(pAliceInp,"*\n*Cut for delta rays by electrons\n");
+ fprintf(pAliceInp,"*Generated from call: SetCut('DCUTE',cut);\n");
+ // -fCutValue[i];
+ // 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; k<nreg; k++) {
+ ireg = reglist[k];
+ fprintf(pAliceInp,"EMFCUT %10.4g%10.1f%10.1f%10.1f%10.1f\n",-fCutValue[i],zero,zero,ireg,ireg);
+ }
+ }
+ fprintf(pAliceInp,"DELTARAY %10.4g%10.3f%10.3f%10.1f%10.1f%10.1f\n",fCutValue[i], 100., 1.03, matMin, matMax, 1.0);
+ fprintf(pAliceInp,"STEPSIZE %10.4g%10.3f%10.3f%10.1f%10.1f\n", 0.1, 1.0, 1.00,
+ Float_t(gGeoManager->GetListOfUVolumes()->GetEntriesFast()-1), 1.0);
+ } // end of if for delta-rays by electrons
+
+
// gammas
// G4 particles: "gamma"
// G3 default value: 0.001 GeV
// gMC ->SetCut("CUTGAM",cut); // cut for gammas
else if (strncmp(&fCutFlag[i][0],"CUTGAM",6) == 0 && global) {
- fprintf(pAliceInp,"*\n*Cut for gamma\n");
- fprintf(pAliceInp,"*Generated from call: SetCut('CUTGAM',cut);\n");
- // -fCutValue[i];
- // 7.0 = lower bound of the particle id-numbers to which the cut-off
- fprintf(pAliceInp,"PART-THR %10.4g%10.1f\n",-fCutValue[i],7.0);
+ fprintf(pAliceInp,"*\n*Cut for gamma\n");
+ fprintf(pAliceInp,"*Generated from call: SetCut('CUTGAM',cut);\n");
+ // -fCutValue[i];
+ // 7.0 = lower bound of the particle id-numbers to which the cut-off
+ fprintf(pAliceInp,"PART-THR %10.4g%10.1f\n",-fCutValue[i],7.0);
}
-
else if (strncmp(&fCutFlag[i][0],"CUTGAM",6) == 0 && !global) {
- fprintf(pAliceInp,"*\n*Cut specific to material for gamma\n");
- fprintf(pAliceInp,"*Generated from call: SetCut('CUTGAM',cut);\n");
- // -fCutValue[i];
- // 7.0 = lower bound of the particle id-numbers to which the cut-off
- fprintf(pAliceInp,"EMFCUT %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", 0., fCutValue[i], zero, matMin, matMax, one);
- }
-
+ fprintf(pAliceInp,"*\n*Cut specific to material for gamma\n");
+ fprintf(pAliceInp,"*Generated from call: SetCut('CUTGAM',cut);\n");
+ // fCutValue[i];
+ // loop over materials for EMFCUT FLUKA cards
+ for (j=0; j < matMax-matMin+1; j++) {
+ Int_t nreg, imat, *reglist;
+ Float_t ireg;
+ imat = (Int_t) matMin + j;
+ reglist = fGeom->GetMaterialList(imat, nreg);
+ // loop over regions of a given material
+ for (Int_t k=0; k<nreg; k++) {
+ ireg = reglist[k];
+ fprintf(pAliceInp,"EMFCUT %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", zero, fCutValue[i], 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(&fCutFlag[i][0],"CUTELE",6) == 0 && global) {
- fprintf(pAliceInp,"*\n*Cut for electrons\n");
- fprintf(pAliceInp,"*Generated from call: SetCut('CUTELE',cut);\n");
- // -fCutValue[i];
- // 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",-fCutValue[i],three,4.0,one);
+ fprintf(pAliceInp,"*\n*Cut for electrons\n");
+ fprintf(pAliceInp,"*Generated from call: SetCut('CUTELE',cut);\n");
+ // -fCutValue[i];
+ // 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",-fCutValue[i],three,4.0,one);
}
-
else if (strncmp(&fCutFlag[i][0],"CUTELE",6) == 0 && !global) {
- fprintf(pAliceInp,"*\n*Cut specific to material for electrons\n");
- fprintf(pAliceInp,"*Generated from call: SetCut('CUTELE',cut);\n");
- // -fCutValue[i];
- // 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,"EMFCUT %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", -fCutValue[i], zero, zero, matMin, matMax, one);
- }
-
+ fprintf(pAliceInp,"*\n*Cut specific to material for electrons\n");
+ fprintf(pAliceInp,"*Generated from call: SetCut('CUTELE',cut);\n");
+ // -fCutValue[i];
+ // 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; k<nreg; k++) {
+ ireg = reglist[k];
+ fprintf(pAliceInp,"EMFCUT %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", -fCutValue[i], zero, zero, ireg, ireg, one);
+ }
+ }
+ } // end of else if for electrons
+
+
// 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(&fCutFlag[i][0],"CUTNEU",6) == 0 && global) {
- fprintf(pAliceInp,"*\n*Cut for neutral hadrons\n");
- fprintf(pAliceInp,"*Generated from call: SetCut('CUTNEU',cut);\n");
+ 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",-fCutValue[i],8.0,9.0);
+ // 8.0 = Neutron
+ // 9.0 = Antineutron
+ fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],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",-fCutValue[i],12.0,12.0);
+ // 12.0 = Kaon zero long
+ // 12.0 = Kaon zero long
+ fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],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",-fCutValue[i],17.0,19.0);
+ // 17.0 = Lambda, 18.0 = Antilambda
+ // 19.0 = Kaon zero short
+ fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],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",-fCutValue[i],22.0,25.0);
+ // 22.0 = Sigma zero, Pion zero, Kaon zero
+ // 25.0 = Antikaon zero
+ fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],22.0,25.0);
- // 32.0 = Antisigma zero
- // 32.0 = Antisigma zero
- fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],32.0,32.0);
+ // 32.0 = Antisigma zero
+ // 32.0 = Antisigma zero
+ fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],32.0,32.0);
- // 34.0 = Xi zero
- // 35.0 = AntiXi zero
- fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],34.0,35.0);
+ // 34.0 = Xi zero
+ // 35.0 = AntiXi zero
+ fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],34.0,35.0);
- // 47.0 = D zero
- // 48.0 = AntiD zero
- fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],47.0,48.0);
+ // 47.0 = D zero
+ // 48.0 = AntiD zero
+ fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],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",-fCutValue[i],53.0,53.0);
+ // 53.0 = Xi_c zero
+ // 53.0 = Xi_c zero
+ fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],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",-fCutValue[i],55.0,56.0);
+ // 55.0 = Xi'_c zero
+ // 56.0 = Omega_c zero
+ fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],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",-fCutValue[i],59.0,59.0);
+ // 59.0 = AntiXi_c zero
+ // 59.0 = AntiXi_c zero
+ fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],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",-fCutValue[i],61.0,62.0);
+ // 61.0 = AntiXi'_c zero
+ // 62.0 = AntiOmega_c zero
+ fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],61.0,62.0);
}
// charged hadrons
// G3 default value: 0.01 GeV
//gMC ->SetCut("CUTHAD",cut); // cut for charged hadrons
else if (strncmp(&fCutFlag[i][0],"CUTHAD",6) == 0 && global) {
- fprintf(pAliceInp,"*\n*Cut for charged hadrons\n");
- fprintf(pAliceInp,"*Generated from call: SetCut('CUTHAD',cut);\n");
+ 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",-fCutValue[i],1.0,2.0);
+ // 1.0 = Proton
+ // 2.0 = Antiproton
+ fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],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",-fCutValue[i],13.0,16.0);
+ // 13.0 = Positive Pion, Negative Pion, Positive Kaon
+ // 16.0 = Negative Kaon
+ fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],13.0,16.0);
- // 20.0 = Negative Sigma
- // 21.0 = Positive Sigma
- fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],20.0,21.0);
+ // 20.0 = Negative Sigma
+ // 21.0 = Positive Sigma
+ fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],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",-fCutValue[i],31.0,33.0,2.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",-fCutValue[i],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",-fCutValue[i],36.0,39.0);
+ // 36.0 = Negative Xi, Positive Xi, Omega minus
+ // 39.0 = Antiomega
+ fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],36.0,39.0);
- // 45.0 = D plus
- // 46.0 = D minus
- fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],45.0,46.0);
+ // 45.0 = D plus
+ // 46.0 = D minus
+ fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],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",-fCutValue[i],49.0,52.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",-fCutValue[i],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",-fCutValue[i],54.0,60.0,6.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",-fCutValue[i],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",-fCutValue[i],57.0,58.0);
+ // 57.0 = Antilambda_c minus
+ // 58.0 = AntiXi_c minus
+ fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],57.0,58.0);
}
// muons
// G3 default value: 0.01 GeV
//gMC ->SetCut("CUTMUO",cut); // cut for mu+, mu-
else if (strncmp(&fCutFlag[i][0],"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",-fCutValue[i],10.0,11.0);
+ 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",-fCutValue[i],10.0,11.0);
}
- // 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(&fCutFlag[i][0],"DCUTE",5) == 0) {
- fprintf(pAliceInp,"*\n*Cut for delta rays by electrons ????????????\n");
- fprintf(pAliceInp,"*Generated from call: SetCut('DCUTE',cut);\n");
- // -fCutValue[i];
- // zero = ignored
- // zero = ignored
- // matMin = lower bound of the material indices in which the respective thresholds apply
- // fLastMaterial = upper bound of the material indices in which the respective thresholds apply
- fprintf(pAliceInp,"EMFCUT %10.4g%10.1f%10.1f%10.1f%10.1f\n",-fCutValue[i],zero,zero,matMin,matMax);
- fprintf(pAliceInp,"DELTARAY %10.4g%10.3f%10.3f%10.1f%10.1f%10.1f\n",fCutValue[i], 100., 1.03, matMin, matMax, 1.0);
- fprintf(pAliceInp,"STEPSIZE %10.4g%10.3f%10.3f%10.1f%10.1f\n", 0.1, 1.0, 1.00,
- Float_t(gGeoManager->GetListOfUVolumes()->GetEntriesFast()-1), 1.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(&fCutFlag[i][0],"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",fCutValue[i]*1.e9,zero,zero,-6.0,64.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",fCutValue[i]*1.e9,zero,zero,-6.0,64.0);
}
else if (global){
- cout << "SetCut for flag=" << &fCutFlag[i][0] << " value=" << fCutValue[i] << " not yet implemented!" << endl;
+ cout << "SetCut for flag=" << &fCutFlag[i][0] << " value=" << fCutValue[i] << " not yet implemented!" << endl;
}
else {
- cout << "SetCut for flag=" << &fCutFlag[i][0] << " value=" << fCutValue[i] << " (material specific) not yet implemented!" << endl;
+ cout << "SetCut for flag=" << &fCutFlag[i][0] << " value=" << fCutValue[i] << " (material specific) not yet implemented!" << endl;
}
} //end of loop over SetCut calls
fclose(pAliceCoreInp);
fclose(pAliceFlukaMat);
fclose(pAliceInp);
- fclose(pGaliceCuts);
} // end of InitPhysics
//______________________________________________________________________________
Bool_t TFluka::IsTrackExiting() const
{
+// True if track is exiting volume
+//
Int_t caller = GetCaller();
if (caller == 12) // bxdraw exiting
return 1;
//______________________________________________________________________________
Int_t TFluka::NSecondaries() const
+
+{
// Number of secondary particles generated in the current step
// FINUC.np = number of secondaries except light and heavy ions
// FHEAVY.npheav = number of secondaries for light and heavy secondary ions
-{
Int_t caller = GetCaller();
if (caller == 6) // valid only after usdraw
return FINUC.np + FHEAVY.npheav;
void TFluka::GetSecondary(Int_t isec, Int_t& particleId,
TLorentzVector& position, TLorentzVector& momentum)
{
+// Copy particles from secondary stack to vmc stack
+//
+
Int_t caller = GetCaller();
if (caller == 6) { // valid only after usdraw
if (isec >= 0 && isec < FINUC.np) {
//______________________________________________________________________________
TMCProcess TFluka::ProdProcess(Int_t) const
+
+{
// Name of the process that has produced the secondary particles
// in the current step
-{
const TMCProcess kIpNoProc = kPNoProcess;
const TMCProcess kIpPDecay = kPDecay;
const TMCProcess kIpPPair = kPPair;
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff