[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

//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.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;
  }

// Get the pointer to the VMC
  TFluka* fluka =  (TFluka*) gMC;
  Int_t verbosityLevel = fluka->GetVerbosityLevel();
  Bool_t debug = (verbosityLevel>=3)?kTRUE:kFALSE;
  fluka->SetTrackIsNew(kTRUE);
//  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; kp <= EMFSTK.npemf - 1; kp++) {

//* 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;

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