From f827183d078cbc9f74ac64894f7160aced2905a1 Mon Sep 17 00:00:00 2001 From: morsch Date: Mon, 3 Mar 2003 09:36:27 +0000 Subject: [PATCH 1/1] User routines stupre and stuprf added. --- TFluka/Femfstk.h | 57 +++++++++++++++ TFluka/Fevtflg.h | 39 ++++++++++ TFluka/libTFluka.pkg | 2 +- TFluka/stupre.cxx | 170 +++++++++++++++++++++++++++++++++++++++++++ TFluka/stuprf.cxx | 112 ++++++++++++++++++++++++++++ 5 files changed, 379 insertions(+), 1 deletion(-) create mode 100644 TFluka/Femfstk.h create mode 100644 TFluka/Fevtflg.h create mode 100644 TFluka/stupre.cxx create mode 100644 TFluka/stuprf.cxx diff --git a/TFluka/Femfstk.h b/TFluka/Femfstk.h new file mode 100644 index 00000000000..44a51af3e64 --- /dev/null +++ b/TFluka/Femfstk.h @@ -0,0 +1,57 @@ +#ifndef FEMFSTK +#define FEMFSTK_H 1 + +#include "Rtypes.h" +#include "cfortran.h" + +#include "Fdimpar.h" +extern "C" { +//*$ create emfstk.add +//*copy emfstk +//* +//*=== emfstk ===========================================================* +//* +//*----------------------------------------------------------------------* +//* * +//* common emfstk (emf stack) for emf * +//* * +//* last change on 08-oct-97 by alfredo ferrari * +//* * +//*----------------------------------------------------------------------* +//* +const Int_t idmemf = mestck; + +typedef struct { + Double_t e[idmemf]; // total energy in MeV + Double_t x[idmemf]; // particle x-coordinate + Double_t y[idmemf]; // particle y-coordinate + Double_t z[idmemf]; // particle z-coordinate + Double_t u[idmemf]; // x direction cosine + Double_t v[idmemf]; // y direction cosine + Double_t w[idmemf]; // z direction cosine + Double_t dnear[idmemf]; // equivalent to GEANT "safety" + Double_t upol[idmemf]; // polarisation in x direction + Double_t vpol[idmemf]; // polarisation in y direction + Double_t wpol[idmemf]; // polarisation in z direction + Double_t usnrml[idmemf]; + Double_t vsnrml[idmemf]; + Double_t wsnrml[idmemf]; + Double_t wt[idmemf]; // weight + Double_t agemf[idmemf]; // age + Double_t espark[idmemf][mkbmx1]; + Int_t iespak[idmemf][mkbmx2]; + Int_t iq[idmemf]; // charge + Int_t ir[idmemf]; // region + Int_t irlatt[idmemf]; // lattice cell + Int_t nhpemf[idmemf]; + Int_t lloemf[idmemf]; // generation number + Int_t louemf[idmemf]; + Int_t np; // number of particles in stack + Int_t npstrt; // EMF stack index before the interaction (since + // the projectile disappears it is also the starting + // index of secondaries) +} emfstkCommon; +#define EMFSTK COMMON_BLOCK(EMFSTK,emfstk) +COMMON_BLOCK_DEF(emfstkCommon,EMFSTK); +} +#endif diff --git a/TFluka/Fevtflg.h b/TFluka/Fevtflg.h new file mode 100644 index 00000000000..4b2f84ad99b --- /dev/null +++ b/TFluka/Fevtflg.h @@ -0,0 +1,39 @@ +extern "C" { +//*$ create evtflg.add +//*copy evtflg +//* +//*=== evtflg ===========================================================* +//* +//*----------------------------------------------------------------------* +//* * +//* event flags: * +//* * +//* created on 19 may 1998 by alfredo ferrari & paola sala * +//* infn - milan * +//* * +//* last change on 13-aug-99 by alfredo ferrari * +//* * +//* * +//*----------------------------------------------------------------------* +//* +//* + +typedef struct { + Int_t lelevt; + Int_t linevt; + Int_t ldecay; + Int_t ldltry; + Int_t lpairp; + Int_t lbrmsp; + Int_t lannrs; + Int_t lannfl; + Int_t lphoel; + Int_t lcmptn; + Int_t lcohsc; + Int_t llensc; + Int_t loppsc; + Int_t ntrcks; +} evtflgCommon; +#define EVTFLG COMMON_BLOCK(EVTFLG,evtflg) +COMMON_BLOCK_DEF(evtflgCommon,EVTFLG); +} diff --git a/TFluka/libTFluka.pkg b/TFluka/libTFluka.pkg index a5a3ef8800a..58d7f31127b 100644 --- a/TFluka/libTFluka.pkg +++ b/TFluka/libTFluka.pkg @@ -1,5 +1,5 @@ # Sources -SRCS:= TFluka.cxx magfld.cxx source.cxx mgdraw.cxx bxdraw.cxx eedraw.cxx endraw.cxx sodraw.cxx usdraw.cxx FlukaVolume.cxx +SRCS:= TFluka.cxx magfld.cxx source.cxx mgdraw.cxx bxdraw.cxx eedraw.cxx endraw.cxx sodraw.cxx usdraw.cxx FlukaVolume.cxx stupre.cxx stuprf.cxx FSRCS:= FLUKA_input.f # Headers diff --git a/TFluka/stupre.cxx b/TFluka/stupre.cxx new file mode 100644 index 00000000000..d4a582f60a8 --- /dev/null +++ b/TFluka/stupre.cxx @@ -0,0 +1,170 @@ +#include +#include "AliRun.h" +#include "TFluka.h" +#ifndef WIN32 +# define stupre stupre_ +#else +# define stupre STUPRE +#endif +// +// Fluka include +#include "Fdimpar.h" //(DIMPAR) fluka include +// Fluka commons +#include "Fdblprc.h" //(DBLPRC) fluka common +#include "Femfstk.h" //(EMFSTK) fluka common +#include "Fevtflg.h" //(EVTFLG) fluka common +#include "Fpaprop.h" //(PAPROP) fluka common +#include "Ftrackr.h" //(TRACKR) fluka common + +//Virtual MC +#include "TFluka.h" +#include "TVirtualMCStack.h" +#include "TVirtualMCApplication.h" +#include "TParticle.h" +#include "TVector3.h" + +extern "C" { +void stupre() +{ +//*----------------------------------------------------------------------* +//* * +//* SeT User PRoperties for Emf particles * +//* * +//*----------------------------------------------------------------------* + + Int_t lbhabh = 0; + if (EVTFLG.ldltry == 1) { + if (EMFSTK.iq[EMFSTK.np-1] * EMFSTK.iq[EMFSTK.np-2] < 0) lbhabh = 1; + } + +// mkbmx1 = dimension for kwb real spare array in fluka stack in DIMPAR +// mkbmx2 = dimension for kwb int. spare array in fluka stack in DIMPAR +// EMFSTK.espark = spare real variables available for +// EMFSTK.iespak = spare integer variables available for +// TRACKR.spausr = user defined spare variables for the current particle +// TRACKR.ispusr = user defined spare flags for the current particle +// EMFSTK.louemf = user flag +// TRACKR.llouse = user defined flag for the current particle + + Int_t npnw, ispr; + for (npnw=EMFSTK.npstrt-1; npnw<=EMFSTK.np-1; npnw++) { + + for (ispr=0; ispr<=mkbmx1-1; ispr++) + EMFSTK.espark[npnw][ispr] = TRACKR.spausr[ispr]; + + for (ispr=0; ispr<=mkbmx2-1; ispr++) + EMFSTK.iespak[npnw][ispr] = TRACKR.ispusr[ispr]; + + EMFSTK.louemf[npnw] = TRACKR.llouse; + } + +// Get the pointer to the VMC + TVirtualMC* fluka = TFluka::GetMC(); +// Get the stack produced from the generator + TVirtualMCStack* cppstack = fluka->GetStack(); + +// EVTFLG.ntrcks = track number +// Increment the track number and put it into the last flag + + Int_t kp; + for (kp=EMFSTK.npstrt-1; npnw<=EMFSTK.np-1; kp++) { + +//* save the parent track number and reset it at each loop + Int_t ncrtrk = EVTFLG.ntrcks; + + Int_t done = 0; + Int_t parent = ncrtrk; + Int_t flukaid = 0; + if (EMFSTK.iq[kp] == -1) flukaid = 3; + else if (EMFSTK.iq[kp] == 0) flukaid = 7; + else if (EMFSTK.iq[kp] == 0) flukaid = 4; + Int_t pdg = fluka->PDGFromId(flukaid); + Double_t e = EMFSTK.e[kp]*emvgev; + Double_t p = sqrt(e*e - PAPROP.am[flukaid+6] * PAPROP.am[flukaid+6]); + Double_t px = p * EMFSTK.u[kp]; + Double_t pz = p * EMFSTK.v[kp]; + Double_t py = p * EMFSTK.w[kp]; + Double_t tof = EMFSTK.agemf[kp]; + Double_t polx = EMFSTK.upol[kp]; + Double_t poly = EMFSTK.vpol[kp]; + Double_t polz = EMFSTK.wpol[kp]; + Double_t vx = EMFSTK.x[kp]; + Double_t vy = EMFSTK.y[kp]; + Double_t vz = EMFSTK.z[kp]; + Double_t weight = EMFSTK.wt[kp]; + Int_t ntr; + TMCProcess mech; + Int_t is = 0; + +//* case of no parent left (pair, photoelectric, annihilation): +//* all secondaries are true + if ((EVTFLG.lpairp == 1) || (EVTFLG.lphoel == 1) || + (EVTFLG.lannfl == 1) || (EVTFLG.lannrs == 1)) { + + if (EVTFLG.lpairp == 1) mech = kPPair; + else if (EVTFLG.lphoel == 1) mech = kPPhotoelectric; + else mech = kPAnnihilation; + cppstack->SetTrack(done, parent, pdg, + px, py, pz, e, vx, vy, vz, tof, + polx, poly, polz, mech, ntr, weight, is); + +cout << endl << " !!! stupre: ntr=" << ntr << endl; + EMFSTK.iespak[kp][mkbmx2-1] = ntr; + } // end of lpairp, lphoel, lannfl, lannrs + +//* Compton: secondary is true only if charged (e+, e-) + else if ((EVTFLG.lcmptn == 1)) { + if (EMFSTK.iq[kp] != 0) { + mech = kPCompton; + cppstack->SetTrack(done, parent, pdg, + px, py, pz, e, vx, vy, vz, tof, + polx, poly, polz, mech, ntr, weight, is); +cout << endl << " !!! stupre: ntr=" << ntr << endl; + EMFSTK.iespak[kp][mkbmx2-1] = ntr; + } + } // end of lcmptn + +//* Bremsstrahlung: true secondary only if charge = 0 (photon) + else if ((EVTFLG.lbrmsp == 1)) { + if (EMFSTK.iq[kp] == 0) { + mech = kPBrem; + cppstack->SetTrack(done, parent, pdg, + px, py, pz, e, vx, vy, vz, tof, + polx, poly, polz, mech, ntr, weight, is); +cout << endl << " !!! stupre: ntr=" << ntr << endl; + EMFSTK.iespak[kp][mkbmx2-1] = ntr; + } + } // end of lbrmsp + +//* Delta ray: If Bhabha, true secondary only if negative (electron) + else if ((EVTFLG.ldltry == 1)) { + if (lbhabh == 1) { + if (EMFSTK.iq[kp] == 0) { + mech = kPDeltaRay; + cppstack->SetTrack(done, parent, pdg, + px, py, pz, e, vx, vy, vz, tof, + polx, poly, polz, mech, ntr, weight, is); + EMFSTK.iespak[kp][mkbmx2-1] = ntr; + } // end of Bhabha + } + +//* Delta ray: Otherwise Moller: true secondary is the electron with +//* lower energy, which has been put higher in the stack + else if (kp == EMFSTK.np) { + mech = kPDeltaRay; + cppstack->SetTrack(done, parent, pdg, + px, py, pz, e, vx, vy, vz, tof, + polx, poly, polz, mech, ntr, weight, is); +cout << endl << " !!! stupre: ntr=" << ntr << endl; + EMFSTK.iespak[kp][mkbmx2-1] = ntr; + } // end of Delta ray + } // end of ldltry + + } // end of loop + +// !!! TO BE CONFIRMED !!! + EVTFLG.ntrcks = EMFSTK.iespak[EMFSTK.np-1][mkbmx2-1]; + +} // end of stupre +} // end of extern "C" + diff --git a/TFluka/stuprf.cxx b/TFluka/stuprf.cxx new file mode 100644 index 00000000000..49ab75b1c7d --- /dev/null +++ b/TFluka/stuprf.cxx @@ -0,0 +1,112 @@ +#include +#include "AliRun.h" +#include "TFluka.h" +#ifndef WIN32 +# define stuprf stuprf_ +#else +# define stuprf STUPRF +#endif +// +// Fluka include +#include "Fdimpar.h" //(DIMPAR) fluka include +// Fluka commons +#include "Fdblprc.h" //(DBLPRC) fluka common +#include "Fevtflg.h" //(EVTFLG) fluka common +#include "Fpaprop.h" //(PAPROP) fluka common +#include "Fstack.h" //(STACK) fluka common +#include "Ftrackr.h" //(TRACKR) fluka common +#include "Ffinuc.h" //(FINUC) fluka common + +//Virtual MC +#include "TFluka.h" +#include "TVirtualMCStack.h" +#include "TVirtualMCApplication.h" +#include "TParticle.h" +#include "TVector3.h" + +extern "C" { +void stuprf(Int_t& ij, Int_t& mreg, + Double_t& xx, Double_t& yy, Double_t& zz, + Int_t& numsec, Int_t& npprmr) +{ +//*----------------------------------------------------------------------* +//* * +//* SeT User PRoperties for Fluka particles * +//* * +//*----------------------------------------------------------------------* + +// STACK.lstack = stack pointer +// STACK.louse = user flag +// TRACKR.llouse = user defined flag for the current particle + STACK.louse[STACK.lstack] = TRACKR.llouse; + +// mkbmx1 = dimension for kwb real spare array in fluka stack in DIMPAR +// mkbmx2 = dimension for kwb int. spare array in fluka stack in DIMPAR +// STACK.sparek = spare real variables available for k.w.burn +// STACK.ispark = spare integer variables available for k.w.burn +// TRACKR.spausr = user defined spare variables for the current particle +// TRACKR.ispusr = user defined spare flags for the current particle + Int_t ispr; + for (ispr=0; ispr<=mkbmx1-1; ispr++) { + STACK.sparek[STACK.lstack][ispr] = TRACKR.spausr[ispr]; + } + for (ispr=0; ispr<=mkbmx2-1; ispr++) { + STACK.ispark[STACK.lstack][ispr] = TRACKR.ispusr[ispr]; + } + +// Get the pointer to the VMC + TVirtualMC* fluka = TFluka::GetMC(); +// Get the stack produced from the generator + TVirtualMCStack* cppstack = fluka->GetStack(); + +// EVTFLG.ntrcks = track number +// Increment the track number and put it into the last flag + if (numsec-1 > npprmr) { +// Now call the SetTrack(...) + Int_t done = 0; + Int_t parent = TRACKR.ispusr[mkbmx2-1]; + Int_t pdg = fluka->PDGFromId(ij); + + Double_t px = FINUC.plr[numsec-1]*FINUC.cxr[numsec-1]; + Double_t pz = FINUC.plr[numsec-1]*FINUC.cyr[numsec-1]; + Double_t py = FINUC.plr[numsec-1]*FINUC.czr[numsec-1]; + Double_t e = FINUC.tki[numsec-1] + PAPROP.am[FINUC.kpart[numsec-1]+6]; + Double_t vx = xx; + Double_t vy = yy; + Double_t vz = zz; + Double_t tof = TRACKR.atrack; + Double_t polx = FINUC.cxrpol[numsec-1]; + Double_t poly = FINUC.cyrpol[numsec-1]; + Double_t polz = FINUC.czrpol[numsec-1]; + + TMCProcess mech = kPHadronic; + if (EVTFLG.ldecay == 1) mech = kPDecay; + else if (EVTFLG.ldltry == 1) mech = kPDeltaRay; + else if (EVTFLG.lpairp == 1) mech = kPPair; + else if (EVTFLG.lbrmsp == 1) mech = kPBrem; + Double_t weight = FINUC.wei[numsec-1]; + Int_t is = 0; + Int_t ntr; + +//virtual void SetTrack(Int_t done, Int_t parent, Int_t pdg, +//Double_t px, Double_t py, Double_t pz, Double_t e, +//Double_t vx, Double_t vy, Double_t vz, Double_t tof, +//Double_t polx, Double_t poly, Double_t polz, +//TMCProcess mech, Int_t& ntr, Double_t weight, +//Int_t is) = 0; + + + cppstack->SetTrack(done, parent, pdg, + px, py, pz, e, + vx, vy, vz, tof, + polx, poly, polz, + mech, ntr, weight, is); + +cout << endl << " !!! stuprf: ntr=" << ntr << endl; + EVTFLG.ntrcks = ntr; + STACK.ispark[STACK.lstack][mkbmx2-1] = EVTFLG.ntrcks; + } // end of if (numsec-1 > npprmr) + +} // end of stuprf +} // end of extern "C" + -- 2.39.3