[u/mrichter/AliRoot.git] / TFluka / stupre.cxx

#include <Riostream.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
#include "Femfrgn.h"  //(EFMRGN) fluka common

//Virtual MC
#include "TFluka.h"
#include "TFlukaCodes.h"
#include "TVirtualMCStack.h"
#include "TVirtualMCApplication.h"
#include "TParticle.h"
#include "TVector3.h"

extern "C" {
void stupre()
{
//*----------------------------------------------------------------------*
//*                                                                      *
//*  SeT User PRoperties for Emf particles                               *
//*                                                                      *
//*----------------------------------------------------------------------*
// Get the pointer to the VMC
  TFluka* fluka =  (TFluka*) gMC;
  static Double_t emassmev = PAPROP.am[kFLUKAelectron + 6] * 1000.;
    
  Int_t lbhabh = 0;
  if (EVTFLG.ldltry == 1) {
    if (EMFSTK.ichemf[EMFSTK.npemf-1] * EMFSTK.ichemf[EMFSTK.npemf-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.npemf-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;
  }


  Int_t verbosityLevel = fluka->GetVerbosityLevel();
  Bool_t debug = (verbosityLevel>=3)?kTRUE:kFALSE;
  fluka->SetTrackIsNew(kTRUE);

// 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; kp <= EMFSTK.npemf - 1; kp++) {
    
// Ckeck transport cut first
    Int_t    ireg   = EMFSTK.iremf[kp];
    Double_t cut    = (TMath::Abs(EMFSTK.ichemf[kp]) == 1) ? EMFRGN.elethr[ireg-1] :  EMFRGN.phothr[ireg-1];
    Double_t e      = EMFSTK.etemf[kp];

    if ((e < cut) 
	&& ( 
	    (EMFSTK.ichemf[kp] ==  0) ||
	    (EMFSTK.ichemf[kp] == -1) ||
	    (EMFSTK.ichemf[kp] ==  1 &&  EMFRGN.phothr[ireg-1] > emassmev)
	    )
	)
    {
	EMFSTK.iespak[kp][mkbmx2-1] = TRACKR.ispusr[mkbmx2-1];
	EMFSTK.iespak[kp][mkbmx2-2] =  0;
	EMFSTK.iespak[kp][mkbmx2-3] = TRACKR.jtrack;
	continue;
    }

// Save the parent track number and reset it at each loop
    Int_t done = 0;
    Int_t parent =  TRACKR.ispusr[mkbmx2-1];
    Int_t flukaid = 0;
    
// Identify particle type
    if      (EMFSTK.ichemf[kp] == -1)  flukaid = kFLUKAelectron;
    else if (EMFSTK.ichemf[kp] ==  0)  flukaid = kFLUKAphoton;
    else if (EMFSTK.ichemf[kp] ==  1)  flukaid = kFLUKApositron;
    

    e *= emvgev;
    Int_t    pdg    = fluka->PDGFromId(flukaid);
    Double_t p      = sqrt(e * e - PAPROP.am[flukaid+6] * PAPROP.am[flukaid+6]);
    Double_t px     = p * EMFSTK.uemf[kp];
    Double_t py     = p * EMFSTK.vemf[kp];
    Double_t pz     = p * EMFSTK.wemf[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.xemf[kp];
    Double_t vy     = EMFSTK.yemf[kp];
    Double_t vz     = EMFSTK.zemf[kp];
    Double_t weight = EMFSTK.wtemf[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->PushTrack(done, parent, pdg,
			   px, py, pz, e, vx, vy, vz, tof,
			   polx, poly, polz, mech, ntr, weight, is);
	if (debug) cout << endl << " !!! stupre (PAIR, ..) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << "energy " << e-PAPROP.am[flukaid+6] << endl;

	EMFSTK.iespak[kp][mkbmx2-1] = ntr;
	EMFSTK.iespak[kp][mkbmx2-2] = 0;
    } // end of lpairp, lphoel, lannfl, lannrs
    
//* Compton: secondary is true only if charged (e+, e-)
    else if ((EVTFLG.lcmptn == 1)) {

	if (EMFSTK.ichemf[kp] != 0) {
	    mech = kPCompton;
	    cppstack->PushTrack(done, parent, pdg,
			       px, py, pz, e, vx, vy, vz, tof,
			       polx, poly, polz, mech, ntr, weight, is);
	    if (debug) cout << endl << " !!! stupre (COMPTON) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
	    EMFSTK.iespak[kp][mkbmx2-1] = ntr;
	    EMFSTK.iespak[kp][mkbmx2-2] = 0;
	}
    } // end of lcmptn
    
//* Bremsstrahlung: true secondary only if charge = 0 (photon)
    else if ((EVTFLG.lbrmsp == 1)) {
	if (EMFSTK.ichemf[kp] == 0) {
	    mech = kPBrem;
	    cppstack->PushTrack(done, parent, pdg,
			       px, py, pz, e, vx, vy, vz, tof,
			       polx, poly, polz, mech, ntr, weight, is);
	    if (debug) cout << endl << " !!! stupre (BREMS) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
	    EMFSTK.iespak[kp][mkbmx2-1] = ntr;
	    EMFSTK.iespak[kp][mkbmx2-2] = 0;
	}
    } // end of lbrmsp
    
//* Delta ray: If Bhabha, true secondary only if negative (electron)
    else if ((EVTFLG.ldltry == 1)) {
	if (lbhabh == 1) {
	    if (EMFSTK.ichemf[kp] == -1) {
		mech = kPDeltaRay;
		cppstack->PushTrack(done, parent, pdg,
				   px, py, pz, e, vx, vy, vz, tof,
				   polx, poly, polz, mech, ntr, weight, is);
		EMFSTK.iespak[kp][mkbmx2-1] = ntr;
		EMFSTK.iespak[kp][mkbmx2-2] = 0;
	   if (debug) cout << endl << " !!! stupre (BHABA) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
	    } // end of Bhabha
	} // lbhabh == 1
	
//* Delta ray: Otherwise Moller: true secondary is the electron with
//*            lower energy, which has been put higher in the stack
	else if (kp == EMFSTK.npemf-1) {
	    mech = kPDeltaRay;
	    cppstack->PushTrack(done, parent, pdg,
			       px, py, pz, e, vx, vy, vz, tof,
			       polx, poly, polz, mech, ntr, weight, is);
	    if (debug) cout << endl << " !!! stupre (Moller) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
	    EMFSTK.iespak[kp][mkbmx2-1] = ntr;
	    EMFSTK.iespak[kp][mkbmx2-2] = 0;
	} // end of Delta ray
    } // end of ldltry
    
  } // end of loop
} // end of stupre
} // end of extern "C"
Commit	Line	Data
f827183d	1	#include <Riostream.h>
f827183d	2	#ifndef WIN32
	3	# define stupre stupre_
	4	#else
	5	# define stupre STUPRE
	6	#endif
	7	//
	8	// Fluka include
	9	#include "Fdimpar.h" //(DIMPAR) fluka include
e3926667	10
f827183d	11	// Fluka commons
	12	#include "Fdblprc.h" //(DBLPRC) fluka common
	13	#include "Femfstk.h" //(EMFSTK) fluka common
	14	#include "Fevtflg.h" //(EVTFLG) fluka common
	15	#include "Fpaprop.h" //(PAPROP) fluka common
	16	#include "Ftrackr.h" //(TRACKR) fluka common
e3926667	17	#include "Femfrgn.h" //(EFMRGN) fluka common
f827183d	18
	19	//Virtual MC
	20	#include "TFluka.h"
d566901f	21	#include "TFlukaCodes.h"
f827183d	22	#include "TVirtualMCStack.h"
	23	#include "TVirtualMCApplication.h"
	24	#include "TParticle.h"
	25	#include "TVector3.h"
	26
	27	extern "C" {
	28	void stupre()
	29	{
	30	//----------------------------------------------------------------------
	31	//* *
	32	//* SeT User PRoperties for Emf particles *
	33	//* *
	34	//----------------------------------------------------------------------
d566901f	35	// Get the pointer to the VMC
	36	TFluka* fluka = (TFluka*) gMC;
	37	static Double_t emassmev = PAPROP.am[kFLUKAelectron + 6] * 1000.;
e3926667	38
f827183d	39	Int_t lbhabh = 0;
f827183d	40	if (EVTFLG.ldltry == 1) {
7dac99f1	41	if (EMFSTK.ichemf[EMFSTK.npemf-1] * EMFSTK.ichemf[EMFSTK.npemf-2] < 0) lbhabh = 1;
f827183d	42	}
d566901f	43
e3926667	44	// mkbmx1 = dimension for kwb real spare array in fluka stack in DIMPAR
e3926667	45	// mkbmx2 = dimension for kwb int. spare array in fluka stack in DIMPAR
f827183d	46	// EMFSTK.espark = spare real variables available for
f827183d	47	// EMFSTK.iespak = spare integer variables available for
e3926667	48	// TRACKR.spausr = user defined spare variables for the current particle
	49	// TRACKR.ispusr = user defined spare flags for the current particle
	50	// EMFSTK.louemf = user flag
	51	// TRACKR.llouse = user defined flag for the current particle
f827183d	52
f827183d	53	Int_t npnw, ispr;
e3926667	54	for (npnw = EMFSTK.npstrt-1; npnw <= EMFSTK.npemf-1; npnw++) {
d566901f	55
	56	for (ispr = 0; ispr <= mkbmx1-1; ispr++)
	57	EMFSTK.espark[npnw][ispr] = TRACKR.spausr[ispr];
	58
	59	for (ispr = 0; ispr <= mkbmx2-1; ispr++)
	60	EMFSTK.iespak[npnw][ispr] = TRACKR.ispusr[ispr];
e3926667	61
d566901f	62	EMFSTK.louemf[npnw] = TRACKR.llouse;
f827183d	63	}
f827183d	64
d566901f	65
2bc4c610	66	Int_t verbosityLevel = fluka->GetVerbosityLevel();
2bc4c610	67	Bool_t debug = (verbosityLevel>=3)?kTRUE:kFALSE;
fbf08100	68	fluka->SetTrackIsNew(kTRUE);
e3926667	69
f827183d	70	// Get the stack produced from the generator
	71	TVirtualMCStack* cppstack = fluka->GetStack();
	72
	73	// EVTFLG.ntrcks = track number
	74	// Increment the track number and put it into the last flag
	75
	76	Int_t kp;
e3926667	77
7dac99f1	78	for (kp = EMFSTK.npstrt - 1; kp <= EMFSTK.npemf - 1; kp++) {
e3926667	79
e3926667	80	// Ckeck transport cut first
d566901f	81	Int_t ireg = EMFSTK.iremf[kp];
d566901f	82	Double_t cut = (TMath::Abs(EMFSTK.ichemf[kp]) == 1) ? EMFRGN.elethr[ireg-1] : EMFRGN.phothr[ireg-1];
e3926667	83	Double_t e = EMFSTK.etemf[kp];
b7940702	84
e3926667	85	if ((e < cut)
	86	&& (
	87	(EMFSTK.ichemf[kp] == 0) \|\|
	88	(EMFSTK.ichemf[kp] == -1) \|\|
11e28b34	89	(EMFSTK.ichemf[kp] == 1 && EMFRGN.phothr[ireg-1] > emassmev)
e3926667	90	)
	91	)
	92	{
b7940702	93	EMFSTK.iespak[kp][mkbmx2-1] = TRACKR.ispusr[mkbmx2-1];
e3926667	94	EMFSTK.iespak[kp][mkbmx2-2] = 0;
18e0cabb	95	EMFSTK.iespak[kp][mkbmx2-3] = TRACKR.jtrack;
e3926667	96	continue;
e3926667	97	}
f827183d	98
d566901f	99	// Save the parent track number and reset it at each loop
610cd189	100	Int_t done = 0;
81ee7392	101	Int_t parent = TRACKR.ispusr[mkbmx2-1];
f827183d	102	Int_t flukaid = 0;
3b8c325d	103
e3926667	104	// Identify particle type
d566901f	105	if (EMFSTK.ichemf[kp] == -1) flukaid = kFLUKAelectron;
	106	else if (EMFSTK.ichemf[kp] == 0) flukaid = kFLUKAphoton;
	107	else if (EMFSTK.ichemf[kp] == 1) flukaid = kFLUKApositron;
e3926667	108
e3926667	109
d566901f	110	e *= emvgev;
e3926667	111	Int_t pdg = fluka->PDGFromId(flukaid);
d15554e0	112	Double_t p = sqrt(e * e - PAPROP.am[flukaid+6] * PAPROP.am[flukaid+6]);
81f1d030	113	Double_t px = p * EMFSTK.uemf[kp];
b7940702	114	Double_t py = p * EMFSTK.vemf[kp];
b7940702	115	Double_t pz = p * EMFSTK.wemf[kp];
d15554e0	116	Double_t tof = EMFSTK.agemf[kp];
	117	Double_t polx = EMFSTK.upol[kp];
	118	Double_t poly = EMFSTK.vpol[kp];
	119	Double_t polz = EMFSTK.wpol[kp];
81f1d030	120	Double_t vx = EMFSTK.xemf[kp];
	121	Double_t vy = EMFSTK.yemf[kp];
	122	Double_t vz = EMFSTK.zemf[kp];
7dac99f1	123	Double_t weight = EMFSTK.wtemf[kp];
7dac99f1	124
f827183d	125	Int_t ntr;
	126	TMCProcess mech;
	127	Int_t is = 0;
	128
	129	//* case of no parent left (pair, photoelectric, annihilation):
	130	//* all secondaries are true
	131	if ((EVTFLG.lpairp == 1) \|\| (EVTFLG.lphoel == 1) \|\|
	132	(EVTFLG.lannfl == 1) \|\| (EVTFLG.lannrs == 1)) {
81ee7392	133
	134	if (EVTFLG.lpairp == 1) mech = kPPair;
	135	else if (EVTFLG.lphoel == 1) mech = kPPhotoelectric;
	136	else mech = kPAnnihilation;
642f15cf	137	cppstack->PushTrack(done, parent, pdg,
81ee7392	138	px, py, pz, e, vx, vy, vz, tof,
81ee7392	139	polx, poly, polz, mech, ntr, weight, is);
e3926667	140	if (debug) cout << endl << " !!! stupre (PAIR, ..) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << "energy " << e-PAPROP.am[flukaid+6] << endl;
d15554e0	141
81ee7392	142	EMFSTK.iespak[kp][mkbmx2-1] = ntr;
48aa0331	143	EMFSTK.iespak[kp][mkbmx2-2] = 0;
f827183d	144	} // end of lpairp, lphoel, lannfl, lannrs
	145
	146	//* Compton: secondary is true only if charged (e+, e-)
	147	else if ((EVTFLG.lcmptn == 1)) {
7dac99f1	148
7dac99f1	149	if (EMFSTK.ichemf[kp] != 0) {
81ee7392	150	mech = kPCompton;
642f15cf	151	cppstack->PushTrack(done, parent, pdg,
81ee7392	152	px, py, pz, e, vx, vy, vz, tof,
81ee7392	153	polx, poly, polz, mech, ntr, weight, is);
2bc4c610	154	if (debug) cout << endl << " !!! stupre (COMPTON) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
81ee7392	155	EMFSTK.iespak[kp][mkbmx2-1] = ntr;
48aa0331	156	EMFSTK.iespak[kp][mkbmx2-2] = 0;
81ee7392	157	}
f827183d	158	} // end of lcmptn
81ee7392	159
f827183d	160	//* Bremsstrahlung: true secondary only if charge = 0 (photon)
f827183d	161	else if ((EVTFLG.lbrmsp == 1)) {
7dac99f1	162	if (EMFSTK.ichemf[kp] == 0) {
81ee7392	163	mech = kPBrem;
642f15cf	164	cppstack->PushTrack(done, parent, pdg,
81ee7392	165	px, py, pz, e, vx, vy, vz, tof,
81ee7392	166	polx, poly, polz, mech, ntr, weight, is);
2bc4c610	167	if (debug) cout << endl << " !!! stupre (BREMS) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
81ee7392	168	EMFSTK.iespak[kp][mkbmx2-1] = ntr;
48aa0331	169	EMFSTK.iespak[kp][mkbmx2-2] = 0;
81ee7392	170	}
f827183d	171	} // end of lbrmsp
81ee7392	172
f827183d	173	//* Delta ray: If Bhabha, true secondary only if negative (electron)
f827183d	174	else if ((EVTFLG.ldltry == 1)) {
81ee7392	175	if (lbhabh == 1) {
7dac99f1	176	if (EMFSTK.ichemf[kp] == -1) {
81ee7392	177	mech = kPDeltaRay;
642f15cf	178	cppstack->PushTrack(done, parent, pdg,
81ee7392	179	px, py, pz, e, vx, vy, vz, tof,
	180	polx, poly, polz, mech, ntr, weight, is);
	181	EMFSTK.iespak[kp][mkbmx2-1] = ntr;
48aa0331	182	EMFSTK.iespak[kp][mkbmx2-2] = 0;
2bc4c610	183	if (debug) cout << endl << " !!! stupre (BHABA) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
81ee7392	184	} // end of Bhabha
d15554e0	185	} // lbhabh == 1
81ee7392	186
f827183d	187	//* Delta ray: Otherwise Moller: true secondary is the electron with
f827183d	188	//* lower energy, which has been put higher in the stack
7dac99f1	189	else if (kp == EMFSTK.npemf-1) {
81ee7392	190	mech = kPDeltaRay;
642f15cf	191	cppstack->PushTrack(done, parent, pdg,
81ee7392	192	px, py, pz, e, vx, vy, vz, tof,
81ee7392	193	polx, poly, polz, mech, ntr, weight, is);
2bc4c610	194	if (debug) cout << endl << " !!! stupre (Moller) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
81ee7392	195	EMFSTK.iespak[kp][mkbmx2-1] = ntr;
48aa0331	196	EMFSTK.iespak[kp][mkbmx2-2] = 0;
81ee7392	197	} // end of Delta ray
f827183d	198	} // end of ldltry
81ee7392	199
f827183d	200	} // end of loop
f827183d	201	} // end of stupre
	202	} // end of extern "C"
	203