//* tinpz = direction cosine of the beam polariz. with respect to*
//* z-axis *
//* polfra = polarization fraction *
-//* nforce = number of the region of forced interaction *
-//* xfor = x-coord. of the starting point of the region nforce *
-//* yfor = y-coord. of the starting point of the region nforce *
-//* zfor = z-coord. of the starting point of the region nforce *
-//* disfor = thickness of the region nforce in cm *
-//* wfor = relative weight of the particle due to forcing *
+//* beawei = weight of the beam particles *
+//* bmaxis(j,i) = j_th component of the i_th axis used to define the*
+//* conventional x,y,z beam reference frame *
+//*!!!!! ATTENTION: in C++ it is the component bmaxis(i,j) !!!!! *
//* ijbeam = beam particle type (see btype in /paprop/) *
//* ijhion = heavy ion type if ijbeam = -2 *
-//* ipbite = flag describing the shape of the momentum *
-//* distribution of the beam *
-//* 0=rectangular, 1=gaussian *
-//* idiv = flag describing the shape of the angular *
-//* divergence distribution of the beam *
-//* 0=rectangular, 1=gaussian *
-//* ixspot = flag describing the shape of the spatial *
-//* distribution of the beam spot in x-direction *
-//* 0=rectangular, 1=gaussian *
-//* iyspot = flag describing the shape of the spatial *
-//* distribution of the beam spot in y-direction *
-//* 0=rectangular, 1=gaussian *
-//* beawei = weight of the beam particles *
+//* ldpgss = true for a gaussian momentum distribution of the *
+//* beam particles, false for a rectangular one *
+//* ldvgss = true for a gaussian angular divergence distribution *
+//* of the beam particles, false for a rectangular one *
+//* ldxgss = true for a gaussian spatial distribution of the beam *
+//* spot in the x-direction, false for a rectangular one *
+//* ldygss = true for a gaussian spatial distribution of the beam *
+//* spot in the y-direction, false for a rectangular one *
//* lbeamc = flag for an annular beam *
//* lpperp = flag for polar. perp. to the beam direction *
-//* lpfrac = flag for interpreting the polar. fraction *
+//* lpfrac = flag for interpreting the polar. fraction
+//* lbaxis = logical flag for using a beam axis frame different *
+//* from the standard one *
//* *
//*----------------------------------------------------------------------*
-// LOGICAL LBEAMC, LPPERP, LPFRAC
+// LBEAMC, LPPERP, LPFRAC, LDPGSS, LDVGSS, LDXGSS, LDYGSS,
+// & LBAXIS
typedef struct {
Double_t pbeam;
Double_t tinpz;
Double_t polfra;
Double_t beawei;
- Double_t xfor;
- Double_t yfor;
- Double_t zfor;
- Double_t disfor;
- Double_t wfor;
+ Double_t bmaxis[3][3];
Int_t ijbeam;
Int_t ijhion;
- Int_t ipbite;
- Int_t idiv;
- Int_t ixspot;
- Int_t iyspot;
- Int_t nforce;
+ Int_t ldpgss;
+ Int_t ldvgss;
+ Int_t ldxgss;
+ Int_t ldygss;
Int_t lbeamc;
Int_t lpperp;
Int_t lpfrac;
+ Int_t lbaxis;
} beamCommon;
#define BEAM COMMON_BLOCK(BEAM,beam)
COMMON_BLOCK_DEF(beamCommon,BEAM);
#ifndef FDIMPAR_H
-#define FDIMPAR_H 1
+#define FDIMPAR_H
#include "cfortran.h"
#include "Rtypes.h"
//* mostck = stack dimension for optical photons *
//* mxprsn = secondary stack dimension for resonance generator *
//* mxpdpm = secondary stack dimension for dpm generators *
+//* Mxpscs = secondary stack dimension overall *
//* mxoutu = maximum number of output units *
//* nallwp = number of allowed particles *
//* nelemx = number of maximum allowed elements of a compound *
//*----------------------------------------------------------------------*
//* *
const Int_t mxxrgn = 10000;
-const Int_t mxxmdf = 210;
-const Int_t mxxmde = 202;
+const Int_t mxxmdf = 510;
+const Int_t mxxmde = 502;
const Int_t mfstck = 2500;
const Int_t mestck = 100;
const Int_t mostck = 2000;
const Int_t mxprsn = 100;
const Int_t mxpdpm = 800;
+const Int_t mxpscs = 1999;
const Int_t mxoutu = 50;
const Int_t nallwp = 64;
const Int_t nelemx = 80;
Double_t wsnrml[mestck];
Double_t wtemf[mestck]; // weight
Double_t agemf[mestck]; // age
+ Double_t cmpemf[mestck];
Double_t espark[mestck][mkbmx1];
Int_t iespak[mestck][mkbmx2];
Int_t ichemf[mestck]; // charge
//*----------------------------------------------------------------------*
//*
const Int_t mxheav = 100;
-const Int_t kxheav = 12;
+const Int_t kxheav = 30;
typedef struct {
Double_t cxheav[mxheav];
+#ifndef FFINUC
+#define FFINUC_H 1
+
+#include "Rtypes.h"
+#include "cfortran.h"
+
+#include "Fdimpar.h"
extern "C" {
//*$ create finuc.add
//*copy finuc
//* *
//*----------------------------------------------------------------------*
//*
-const Int_t mxp = 999;
+const Int_t mxp = mxpscs;
//*
typedef struct {
#define FINUC COMMON_BLOCK(FINUC,finuc)
COMMON_BLOCK_DEF(finucCommon,FINUC);
}
+#endif
//*----------------------------------------------------------------------*
//*
+const Int_t mxgnpr = 33;
typedef struct {
Double_t am[nallwp+7]; //(-6:NALLWP)
Double_t amdisc[nallwp+7]; //(-6:NALLWP)
typedef struct {
Char_t btype[nallwp+7][8]; //(-6:NALLWP)
- Char_t genpar[30][8]; //(30)
+ Char_t genpar[mxgnpr][8]; //(30)
} chpprpCommon;
#define CHPPRP COMMON_BLOCK(CHPPRP,chpprp)
COMMON_BLOCK_DEF(chpprpCommon,CHPPRP);
//* raddly = delay (s) in production wrt the nominal primary "0" *
//* time for particle produced in radioactive decays *
//* (i.e. those coming from decays of daughter isotopes) *
+//* cmpath = cumulative path travelled by the particle since it *
+//* was produced (cm)
//* sparek = spare real variables available for k.w.burn *
//* ispark = spare integer variables available for k.w.burn *
//* ilo = type of the particle (see btype in /paprop/) *
Double_t agestk[mfstck+1]; //(0:MFSTCK)
Double_t aknshr[mfstck+1]; //(0:MFSTCK)
Double_t raddly[mfstck+1]; //(0:MFSTCK)
+ Double_t cmpath[mfstck+1]; //(0:MFSTCK)
Double_t sparek[mfstck+1][mkbmx1]; //(MKBMX1,0:MFSTCK)
Int_t ispark[mfstck+1][mkbmx2]; //(MKBMX2,0:MFSTCK)
Int_t ilo[mfstck+1]; //(0:MFSTCK)
//* ttrack = length of the ith track segment *
//* 1 < j < mtrack *
//* dtrack = energy deposition of the jth deposition event *
+//* dptrck = momentum loss of the jth deposition event *
//* *
//* jtrack = identity number of the particle *
//* etrack = total energy of the particle *
//* biasing techniques are used (for example inelastic *
//* interaction length biasing) *
//* ctrack = total curved path *
+//* cmtrck = cumulative curved path since particle birth *
//* zfftrk = <z_eff> of the particle *
//* zfrttk = actual z_eff of the particle *
//* atrack = age of the particle *
Double_t ztrack[mxtrck+1];
Double_t ttrack[mxtrck];
Double_t dtrack[mxtrck];
+ Double_t dptrck[mxtrck][3];
Double_t etrack;
Double_t ptrack;
Double_t cxtrck;
Double_t zfrttk;
Double_t atrack;
Double_t ctrack;
+ Double_t cmtrck;
Double_t akshrt;
Double_t aklong;
Double_t wscrng;
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 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);
}
}
} // 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;
}
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 == 1 || caller == 3 || caller == 4 || caller == 6) //bxdraw,endraw,mgdraw,usdraw
+ return TRACKR.cmtrck;
else
return -1.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
+// 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;
}
Bool_t TFluka::IsTrackInside() const
(TVirtualMCApplication::Instance())->Stepping();
fluka->SetCaller(11);
fluka->SetTrackIsEntering();
- printf("bxdraw (en) \n");
+ printf("bxdraw (en) mreg=%d newreg=%d \n",mreg,newreg);
fluka->SetMreg(newreg);
(TVirtualMCApplication::Instance())->Stepping();
// fluka->SetCaller(1);
// projectile
gener->SetProjectile("P", 1, 1);
gener->SetTarget ("A", 16, 8);
+ gener->SetBoostLHC(0);
+
// tell hijing to keep the full parent child chain
gener->KeepFullEvent();
// enable jet quenching
gener->Generate();
Int_t npart = stack->GetNprimary();
// Vertex
- Float_t za = 4000. * gRandom->Rndm() - 2000.;
+ Float_t za = 4000. * gRandom->Rndm() -2000.;
// Direction
- Float_t dir = (gRandom->Rndm() < 0.5) ? 1. : -1.;
-
+// Float_t dir = (za < 0.) ? 1. : -1.;
+ Float_t dir = (gRandom->Rndm() < 0.5)? 1. : -1;
for (Int_t part=0; part<npart; part++) {
particle = stack->Particle(part);
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff