sInputFileName(""),
fTrackIsEntering(0),
fTrackIsExiting(0),
+ fTrackIsNew(0),
fDetector(0),
fCurrentFlukaRegion(-1)
{
//____________________________________________________________________________
void TFluka::Init() {
+ FGeometryInit* geominit = FGeometryInit::GetInstance();
if (fVerbosityLevel >=3)
cout << "==> TFluka::Init() called." << endl;
- cout << "\t* InitPhysics() - Prepare input file to be called" << endl;
- InitPhysics(); // prepare input file
- cout << "\t* InitPhysics() - Prepare input file called" << endl;
+ cout << "\t* InitPhysics() - Prepare input file to be called" << endl;
+ geominit->Init();
+ // now we have G4 geometry created and we have to patch alice.inp
+ // with the material mapping file FlukaMat.inp
+ InitPhysics(); // prepare input file with the current physics settings
+ cout << "\t* InitPhysics() - Prepare input file was called" << endl;
if (fVerbosityLevel >=2)
cout << "\t* Changing lfdrtr = (" << (GLOBAL.lfdrtr?'T':'F')
Int_t TFluka::PDGFromId(Int_t id) const
{
-
//
// Return PDG code and pseudo ENDF code from Fluka code
- //IPTOKP array goes from official to internal
+ // IPTOKP array goes from official to internal
if (id == -1) {
// Cerenkov photon
printf("\n PDGFromId: Cerenkov Photon \n");
return 50000050;
}
-
+// Error id
if (id == 0) {
if (fVerbosityLevel >= 1)
printf("PDGFromId: Error id = 0\n");
return -1;
}
-
+// Good id
Int_t intfluka = GetFlukaIPTOKP(id);
if (intfluka == 0) {
if (fVerbosityLevel >= 1)
}
if (fVerbosityLevel >= 3)
printf("mpdgha called with %d %d \n", id, intfluka);
+ // MPDGHA() goes from fluka internal to pdg.
return mpdgha(intfluka);
}
// !!! it should be available from Flugg !!!
Int_t i, j, k;
Double_t fCut;
- Float_t fLastMaterial = 31.0;
+ FGeometryInit* geominit = FGeometryInit::GetInstance();
+ Float_t fLastMaterial = geominit->GetLastMaterialIndex();
+ printf(" last FLUKA material is %g\n", fLastMaterial);
// construct file names
- TString sAliceInp = getenv("ALICE_ROOT");
- sAliceInp +="/TFluka/input/";
- TString sAliceCoreInp = sAliceInp;
- sAliceInp += GetInputFileName();
+ TString sAliceCoreInp = getenv("ALICE_ROOT");
+ sAliceCoreInp +="/TFluka/input/";
+ TString sAliceTmp = "flukaMat.inp";
+ TString sAliceInp = GetInputFileName();
sAliceCoreInp += GetCoreInputFileName();
ifstream AliceCoreInp(sAliceCoreInp.Data());
+ ifstream AliceFlukaMat(sAliceTmp.Data());
ofstream AliceInp(sAliceInp.Data());
-// copy core input file until (not included) START card
+// copy core input file
Char_t sLine[255];
Float_t fEventsPerRun;
+
+ while (AliceCoreInp.getline(sLine,255)) {
+ if (strncmp(sLine,"GEOEND",6) != 0)
+ AliceInp << sLine << endl; // copy until GEOEND card
+ else {
+ AliceInp << "GEOEND" << endl; // add GEOEND card
+ goto flukamat;
+ }
+ } // end of while until GEOEND card
+
+flukamat:
+ while (AliceFlukaMat.getline(sLine,255)) { // copy flukaMat.inp file
+ AliceInp << sLine << endl;
+ }
+
while (AliceCoreInp.getline(sLine,255)) {
if (strncmp(sLine,"START",5) != 0)
AliceInp << sLine << endl;
sscanf(sLine+10,"%10f",&fEventsPerRun);
goto fin;
}
- } //end of while
+ } //end of while until START card
fin:
// in G3 the process control values meaning can be different for
}
} // else if (strncmp(&sProcessFlag[i][0],"PHOT",4) == 0)
+
// Rayleigh scattering
// G3 default value: 0
// G4 process: G4OpRayleigh
AliceInp << endl;
}
} // end of else if (strncmp(&sProcessFlag[i][0],"RAYL",4) == 0)
-
- else { // processes not yet treated
+
+ // 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(&sProcessFlag[i][0],"SYNC",4) == 0) {
+ AliceInp << "*";
+ AliceInp << endl;
+ AliceInp << "*Synchrotron radiation generation is NOT implemented in FLUKA";
+ AliceInp << endl;
+ AliceInp << "*No FLUKA card generated";
+ AliceInp << endl;
+ }
+
// 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(&sProcessFlag[i][0],"AUTO",4) == 0) {
+ AliceInp << "*";
+ AliceInp << endl;
+ AliceInp << "*Automatic calculation of tracking medium parameters is always ON in FLUKA";
+ AliceInp << endl;
+ AliceInp << "*No FLUKA card generated";
+ AliceInp << endl;
+ }
+
+
+ // 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(&sProcessFlag[i][0],"STRA",4) == 0) {
+ if (iProcessValue[i] == 0 || iProcessValue[i] == 2 || iProcessValue[i] == 3) {
+ AliceInp << "*";
+ AliceInp << endl;
+ AliceInp << "*Ionization energy losses calculation is activated";
+ AliceInp << endl;
+ AliceInp << "*Generated from call: SetProcess('STRA',n);, n=0,1,2";
+ AliceInp << endl;
+ AliceInp << setw(10) << "IONFLUCT ";
+ AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
+ AliceInp << setw(10) << 1.0; // restricted energy loss fluctuations
+ // (for hadrons and muons) switched on
+ AliceInp << setw(10) << 1.0; // restricted energy loss fluctuations
+ // (for e+ and e-) switched on
+ AliceInp << setw(10) << 1.0; // minimal accuracy
+ AliceInp << setw(10) << 3.0; // upper bound of the material indices in
+ // which the respective thresholds apply
+ AliceInp << setprecision(2);
+ AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply
+ AliceInp << setprecision(1);
+ AliceInp << setw(10) << 1.0; // step length in assigning indices
+ AliceInp << endl;
+ }
+ else {
+ AliceInp << "*";
+ AliceInp << endl;
+ AliceInp << "*Illegal flag value in SetProcess('STRA',?) call.";
+ AliceInp << endl;
+ AliceInp << "*No FLUKA card generated";
+ AliceInp << endl;
+ }
+ } // else if (strncmp(&sProcessFlag[i][0],"STRA",4) == 0)
+
+
+
+
+ else { // processes not yet treated
// light photon absorption (Cerenkov photons)
// it is turned on when Cerenkov process is turned on
// gMC ->SetProcess("LABS",2); // ??? Cerenkov light absorption
- // 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
-
- // 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
cout << "SetProcess for flag=" << &sProcessFlag[i][0] << " value=" << iProcessValue[i] << " not yet implemented!" << endl;
}
AliceInp << setw(10) << "STOP ";
AliceInp << endl;
-}
+} // end of InitPhysics
void TFluka::SetMaxStep(Double_t)
// TRACKR.ytrack = y-position of the last point
// TRACKR.ztrack = z-position of the last point
Int_t caller = GetCaller();
- if (caller == 1 || caller == 3 || caller == 6 || caller == 11 || caller == 12) { //bxdraw,endraw,usdraw
+ if (caller == 3 || caller == 6 || caller == 11 || caller == 12) { //bxdraw,endraw,usdraw
position.SetX(GetXsco());
position.SetY(GetYsco());
position.SetZ(GetZsco());
// TRACKR.ytrack = y-position of the last point
// TRACKR.ztrack = z-position of the last point
Int_t caller = GetCaller();
- if (caller == 1 || caller == 3 || caller == 6 || caller == 11 || caller == 12) { //bxdraw,endraw,usdraw
+ if (caller == 3 || caller == 6 || caller == 11 || caller == 12) { //bxdraw,endraw,usdraw
x = GetXsco();
y = GetYsco();
z = GetZsco();
// Return the length in centimeters of the current step
// TRACKR.ctrack = total curved path
Int_t caller = GetCaller();
- if (caller == 1 || caller == 3 || caller == 6) //bxdraw,endraw,usdraw
+ if (caller == 11 || caller==12 || caller == 3 || caller == 6) //bxdraw,endraw,usdraw
return 0.0;
else if (caller == 4) //mgdraw
return TRACKR.ctrack;
Double_t TFluka::TrackLength() const
{
-// Still wrong !!!
-// This is the sum of substeps !!!
-// TRACKR.ctrack = total curved path of the current step
-// Sum of the substeps is identical to TRACKR.ctrack if the is no mag. field
-// The sum of all step length starting from the beginning of the track
-// for the time being returns only the length in centimeters of the current step
- Double_t sum = 0;
+// TRACKR.cmtrck = cumulative curved path since particle birth
Int_t caller = GetCaller();
- if (caller == 1 || caller == 3 || caller == 4 || caller == 6) { //bxdraw,endraw,mgdraw,usdraw
- for ( Int_t j=0;j<TRACKR.ntrack;j++) {
- sum +=TRACKR.ttrack[j];
- }
- return sum;
- }
+ if (caller == 111 || caller==12 || caller == 3 || caller == 4 || caller == 6) //bxdraw,endraw,mgdraw,usdraw
+ return TRACKR.cmtrck;
else
return -1.0;
}
// Return the current time of flight of the track being transported
// TRACKR.atrack = age of the particle
Int_t caller = GetCaller();
- if (caller == 1 || caller == 3 || caller == 4 || caller == 6) //bxdraw,endraw,mgdraw,usdraw
+ if (caller == 11 || caller==12 || caller == 3 || caller == 4 || caller == 6) //bxdraw,endraw,mgdraw,usdraw
return TRACKR.atrack;
else
return -1;
// if TRACKR.ntrack > 0, TRACKR.mtrack > 0:
// -->energy loss distributed along the track
// TRACKR.dtrack = energy deposition of the jth deposition even
+
+ // If coming from bxdraw we have 2 steps of 0 length and 0 edep
+ Int_t caller = GetCaller();
+ if (caller == 11 || caller==12) return 0.0;
Double_t sum = 0;
for ( Int_t j=0;j<TRACKR.mtrack;j++) {
sum +=TRACKR.dtrack[j];
// PAPROP.am = particle mass in GeV
// TRACKR.jtrack = identity number of the particle
Int_t caller = GetCaller();
- if (caller != 2) // not eedraw
+ if (caller != 2) { // not eedraw
+// cout << "JTRACK=" << TRACKR.jtrack << " mass=" << PAPROP.am[TRACKR.jtrack+6] << endl;
return PAPROP.am[TRACKR.jtrack+6];
+ }
else
return -1000.0;
}
//
Bool_t TFluka::IsNewTrack() const
{
-// ???????????????,
-// True if the track is not at the boundary of the current volume
-// Not true in some cases in bxdraw - to be solved
+// Return true for the first call of Stepping()
+/*
+// True if the track has positive cummulative length
Int_t caller = GetCaller();
- if (caller == 1)
- return 1; // how to handle double step ?????????????
+ if (caller != 2) { // not eedraw
+ if (TRACKR.cmtrck > 0.0)
+ return 1;
+ else
+ return 0;
+ }
else
- return 0; // ??????????????
+ return 0;
+*/
+ return fTrackIsNew;
}
Bool_t TFluka::IsTrackInside() const
// it will be shortened to reach only the boundary.
// Therefore IsTrackInside() is always true.
Int_t caller = GetCaller();
- if (caller == 1) // bxdraw
+ if (caller == 11 || caller==12) // bxdraw
return 0;
else
return 1;
Warning("GetSecondary","no secondaries available");
} // end of GetSecondary
-TMCProcess TFluka::ProdProcess() const
+TMCProcess TFluka::ProdProcess(Int_t) const
// Name of the process that has produced the secondary particles
// in the current step
{
return name;
}
-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
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff