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

#include <Riostream.h>
#include "TCallf77.h"      //For the fortran calls
#ifndef WIN32
# define stuprf stuprf_
# define usrdci usrdci_
#else
# define stuprf STUPRF
# define usrdci USRDCI
#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 "Fflkstk.h"  //(FLKSTK)  fluka common
#include "Ftrackr.h"  //(TRACKR) fluka common
#include "Fgenstk.h"  //(GENSTK)  fluka common
#include "Ffheavy.h"  //(FHEAVY)  fluka common


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

extern "C" {

    void type_of_call usrdci(int&, int&, int&, int&);

    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                             *
//*                                                                      *
//*----------------------------------------------------------------------*
// Get the pointer to the VMC
  TFluka* fluka =  (TFluka*) gMC;

// FLKSTK.npflka  = stack pointer
// FLKSTK.louse   = user flag
// TRACKR.llouse = user defined flag for the current particle
    
   FLKSTK.louse[FLKSTK.npflka] = 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
// FLKSTK.sparek  = spare real variables available for k.w.burn
// FLKSTK.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;
  if (numsec <= npprmr) {
      for (ispr = 0; ispr <= mkbmx1 - 1; ispr++) {
	  FLKSTK.sparek[FLKSTK.npflka][ispr] = TRACKR.spausr[ispr];
      }  
      for (ispr = 0; ispr <= mkbmx2 - 1; ispr++) {
	  FLKSTK.ispark[FLKSTK.npflka][ispr] = TRACKR.ispusr[ispr];
      }  
  }
  
  // save parent info
  FLKSTK.ispark[FLKSTK.npflka][mkbmx2 - 3] = TRACKR.jtrack;              // fluka particle id
  FLKSTK.ispark[FLKSTK.npflka][mkbmx2 - 4] = TRACKR.ispusr[mkbmx2 - 1];  // current track number
  FLKSTK.ispark[FLKSTK.npflka][mkbmx2 - 5] = npprmr;                     // flag npprmr>0


  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
// was numsec -1
// clarify with Alberto
  
//  npprmr > 0, the secondary being loaded is actually still the interacting
//  particle (it can happen in some biasing situations)

  if (numsec > npprmr) {
// Now call the PushTrack(...)
    Int_t done = 0;

    Int_t parent  = TRACKR.ispusr[mkbmx2-1];
    Int_t kpart   = GENSTK.kpart[numsec-1];
    Int_t pdg;
    Int_t a = 1;
    Int_t z = 1;
    Int_t ism = 0;
    if (kpart < -6) {
//	kpart = -kpart;
//	z = kpart / 100000;
//	a = kpart - z * 100000;
//	a /= 100;
	usrdci(kpart, a, z, ism);
	pdg = fluka->GetIonPdg(z, a);
	fluka->AddIon(a, z);
    } else {
	pdg  = fluka->PDGFromId(kpart);
    }
    
    Double_t px = GENSTK.plr[numsec-1] * GENSTK.cxr[numsec-1];
    Double_t py = GENSTK.plr[numsec-1] * GENSTK.cyr[numsec-1];
    Double_t pz = GENSTK.plr[numsec-1] * GENSTK.czr[numsec-1];
    Double_t e;
    
    if (a == 1) {
	e  = GENSTK.tki[numsec-1] + PAPROP.am[GENSTK.kpart[numsec-1]+6];
    } else {
	e  = GENSTK.tki[numsec-1] + Double_t(a) * 0.93149 + 8.071e-3;
    }
    
    
    Double_t vx = xx;
    Double_t vy = yy;
    Double_t vz = zz;

    Double_t tof  = TRACKR.atrack;
    Double_t polx = GENSTK.cxrpol[numsec-1];
    Double_t poly = GENSTK.cyrpol[numsec-1];
    Double_t polz = GENSTK.czrpol[numsec-1];
    
    TMCProcess mech = kPHadronic;
    if (EVTFLG.ldecay == 1) {
        mech = kPDecay;
        if (debug) cout << endl << "Decay" << endl;
        FLKSTK.nlattc[FLKSTK.npflka] = TRACKR.lt1trk;
    } else if (EVTFLG.ldltry == 1) {
        mech = kPDeltaRay;
        if( fluka->GetIcode() == kKASHEA ) {
           //  For all interactions secondaries are put on GENSTK common (kp=1,np)
           //  but for KASHEA delta ray generation where only the secondary elec-
           //  tron is present and stacked on FLKSTK common for kp=lstack
           pdg  = fluka->PDGFromId( FLKSTK.iloflk[FLKSTK.npflka] );
           px   = FLKSTK.pmoflk[FLKSTK.npflka] * FLKSTK.txflk[FLKSTK.npflka];
           py   = FLKSTK.pmoflk[FLKSTK.npflka] * FLKSTK.tyflk[FLKSTK.npflka];
           pz   = FLKSTK.pmoflk[FLKSTK.npflka] * FLKSTK.tzflk[FLKSTK.npflka];
           e    = FLKSTK.tkeflk[FLKSTK.npflka] + PAPROP.am[FLKSTK.iloflk[FLKSTK.npflka]+6];
           polx = FLKSTK.txpol[FLKSTK.npflka];
           poly = FLKSTK.typol[FLKSTK.npflka];
           polz = FLKSTK.tzpol[FLKSTK.npflka];
           if (debug) cout << endl << "Delta Ray from KASHEA...." << " pdg from FLKSTK=" << pdg << endl;
        } else {
	    if (debug) cout << endl << "Delta Ray" << endl;
	}
    } else if (EVTFLG.lpairp == 1) {
        mech = kPPair;
        if (debug) cout << endl << "Pair Production" << endl;
    } else if (EVTFLG.lbrmsp == 1) {
        mech = kPBrem;
        if (debug) cout << endl << "Bremsstrahlung" << endl;
    }

    Double_t weight = GENSTK.wei[numsec-1];
    Int_t is = 0;
    Int_t ntr;
    // 
    // Save particle in VMC stack
    cppstack->PushTrack(done, parent, pdg,
			px, py, pz, e,
			vx, vy, vz, tof,
			polx, poly, polz,
			mech, ntr, weight, is);
    if (debug) {
	
	cout << endl << " !!! stuprf: ntr=" << ntr << " pdg " << pdg << " parent=" << parent
	     << " parent_pdg="<< fluka->PDGFromId(TRACKR.jtrack) << " numsec "
	     << numsec << " npprmr " << npprmr << " icode=" << fluka->GetIcode() 
	     << "kin. energy [keV] = " <<  GENSTK.tki[numsec-1] * 1.e6 << endl << endl;
	
    }
    

//
//  Save current track number
    FLKSTK.ispark[FLKSTK.npflka][mkbmx2-1] = ntr;
    FLKSTK.ispark[FLKSTK.npflka][mkbmx2-2] = 0;
  } // end of if (numsec > npprmr)
//  else {
//     if(debug) {
//        cout << endl << " !!! stuprf: skipping pushtrack   track=" << TRACKR.ispusr[mkbmx2-1]
//              << " pdg " << fluka->PDGFromId(TRACKR.jtrack) << " numsec=" << numsec<< " npprmr=" << npprmr
//              << " GENSTK pdg=" << fluka->PDGFromId(GENSTK.kpart[numsec-1]) << endl;
//     }
//  }
} // end of stuprf

} // end of extern "C"
Commit	Line	Data
	1	#include <Riostream.h>
	2	#include "TCallf77.h" //For the fortran calls
	3	#ifndef WIN32
	4	# define stuprf stuprf_
	5	# define usrdci usrdci_
	6	#else
	7	# define stuprf STUPRF
	8	# define usrdci USRDCI
	9	#endif
	10	//
	11	// Fluka include
	12	#include "Fdimpar.h" //(DIMPAR) fluka include
	13	// Fluka commons
	14	#include "Fdblprc.h" //(DBLPRC) fluka common
	15	#include "Fevtflg.h" //(EVTFLG) fluka common
	16	#include "Fpaprop.h" //(PAPROP) fluka common
	17	#include "Fflkstk.h" //(FLKSTK) fluka common
	18	#include "Ftrackr.h" //(TRACKR) fluka common
	19	#include "Fgenstk.h" //(GENSTK) fluka common
	20	#include "Ffheavy.h" //(FHEAVY) fluka common
	21
	22
	23	//Virtual MC
	24	#include "TFluka.h"
	25	#include "TVirtualMCStack.h"
	26	#include "TVirtualMCApplication.h"
	27	#include "TParticle.h"
	28	#include "TVector3.h"
	29
	30	extern "C" {
	31
	32	void type_of_call usrdci(int&, int&, int&, int&);
	33
	34	void stuprf(Int_t& /ij/, Int_t& /mreg/,
	35	Double_t& xx, Double_t& yy, Double_t& zz,
	36	Int_t& numsec, Int_t& npprmr)
	37	{
	38	//----------------------------------------------------------------------
	39	//* *
	40	//* SeT User PRoperties for Fluka particles *
	41	//* *
	42	//----------------------------------------------------------------------
	43	// Get the pointer to the VMC
	44	TFluka* fluka = (TFluka*) gMC;
	45
	46	// FLKSTK.npflka = stack pointer
	47	// FLKSTK.louse = user flag
	48	// TRACKR.llouse = user defined flag for the current particle
	49
	50	FLKSTK.louse[FLKSTK.npflka] = TRACKR.llouse;
	51
	52	// mkbmx1 = dimension for kwb real spare array in fluka stack in DIMPAR
	53	// mkbmx2 = dimension for kwb int. spare array in fluka stack in DIMPAR
	54	// FLKSTK.sparek = spare real variables available for k.w.burn
	55	// FLKSTK.ispark = spare integer variables available for k.w.burn
	56	// TRACKR.spausr = user defined spare variables for the current particle
	57	// TRACKR.ispusr = user defined spare flags for the current particle
	58	Int_t ispr;
	59	if (numsec <= npprmr) {
	60	for (ispr = 0; ispr <= mkbmx1 - 1; ispr++) {
	61	FLKSTK.sparek[FLKSTK.npflka][ispr] = TRACKR.spausr[ispr];
	62	}
	63	for (ispr = 0; ispr <= mkbmx2 - 1; ispr++) {
	64	FLKSTK.ispark[FLKSTK.npflka][ispr] = TRACKR.ispusr[ispr];
	65	}
	66	}
	67
	68	// save parent info
	69	FLKSTK.ispark[FLKSTK.npflka][mkbmx2 - 3] = TRACKR.jtrack; // fluka particle id
	70	FLKSTK.ispark[FLKSTK.npflka][mkbmx2 - 4] = TRACKR.ispusr[mkbmx2 - 1]; // current track number
	71	FLKSTK.ispark[FLKSTK.npflka][mkbmx2 - 5] = npprmr; // flag npprmr>0
	72
	73
	74	Int_t verbosityLevel = fluka->GetVerbosityLevel();
	75	Bool_t debug = (verbosityLevel>=3)? kTRUE:kFALSE;
	76
	77	fluka->SetTrackIsNew(kTRUE);
	78	// TVirtualMC* fluka = TFluka::GetMC();
	79	// Get the stack produced from the generator
	80	TVirtualMCStack* cppstack = fluka->GetStack();
	81
	82	// EVTFLG.ntrcks = track number
	83	// Increment the track number and put it into the last flag
	84	// was numsec -1
	85	// clarify with Alberto
	86
	87	// npprmr > 0, the secondary being loaded is actually still the interacting
	88	// particle (it can happen in some biasing situations)
	89
	90	if (numsec > npprmr) {
	91	// Now call the PushTrack(...)
	92	Int_t done = 0;
	93
	94	Int_t parent = TRACKR.ispusr[mkbmx2-1];
	95	Int_t kpart = GENSTK.kpart[numsec-1];
	96	Int_t pdg;
	97	Int_t a = 1;
	98	Int_t z = 1;
	99	Int_t ism = 0;
	100	if (kpart < -6) {
	101	// kpart = -kpart;
	102	// z = kpart / 100000;
	103	// a = kpart - z * 100000;
	104	// a /= 100;
	105	usrdci(kpart, a, z, ism);
	106	pdg = fluka->GetIonPdg(z, a);
	107	fluka->AddIon(a, z);
	108	} else {
	109	pdg = fluka->PDGFromId(kpart);
	110	}
	111
	112	Double_t px = GENSTK.plr[numsec-1] * GENSTK.cxr[numsec-1];
	113	Double_t py = GENSTK.plr[numsec-1] * GENSTK.cyr[numsec-1];
	114	Double_t pz = GENSTK.plr[numsec-1] * GENSTK.czr[numsec-1];
	115	Double_t e;
	116
	117	if (a == 1) {
	118	e = GENSTK.tki[numsec-1] + PAPROP.am[GENSTK.kpart[numsec-1]+6];
	119	} else {
	120	e = GENSTK.tki[numsec-1] + Double_t(a) * 0.93149 + 8.071e-3;
	121	}
	122
	123
	124	Double_t vx = xx;
	125	Double_t vy = yy;
	126	Double_t vz = zz;
	127
	128	Double_t tof = TRACKR.atrack;
	129	Double_t polx = GENSTK.cxrpol[numsec-1];
	130	Double_t poly = GENSTK.cyrpol[numsec-1];
	131	Double_t polz = GENSTK.czrpol[numsec-1];
	132
	133	TMCProcess mech = kPHadronic;
	134	if (EVTFLG.ldecay == 1) {
	135	mech = kPDecay;
	136	if (debug) cout << endl << "Decay" << endl;
	137	FLKSTK.nlattc[FLKSTK.npflka] = TRACKR.lt1trk;
	138	} else if (EVTFLG.ldltry == 1) {
	139	mech = kPDeltaRay;
	140	if( fluka->GetIcode() == kKASHEA ) {
	141	// For all interactions secondaries are put on GENSTK common (kp=1,np)
	142	// but for KASHEA delta ray generation where only the secondary elec-
	143	// tron is present and stacked on FLKSTK common for kp=lstack
	144	pdg = fluka->PDGFromId( FLKSTK.iloflk[FLKSTK.npflka] );
	145	px = FLKSTK.pmoflk[FLKSTK.npflka] * FLKSTK.txflk[FLKSTK.npflka];
	146	py = FLKSTK.pmoflk[FLKSTK.npflka] * FLKSTK.tyflk[FLKSTK.npflka];
	147	pz = FLKSTK.pmoflk[FLKSTK.npflka] * FLKSTK.tzflk[FLKSTK.npflka];
	148	e = FLKSTK.tkeflk[FLKSTK.npflka] + PAPROP.am[FLKSTK.iloflk[FLKSTK.npflka]+6];
	149	polx = FLKSTK.txpol[FLKSTK.npflka];
	150	poly = FLKSTK.typol[FLKSTK.npflka];
	151	polz = FLKSTK.tzpol[FLKSTK.npflka];
	152	if (debug) cout << endl << "Delta Ray from KASHEA...." << " pdg from FLKSTK=" << pdg << endl;
	153	} else {
	154	if (debug) cout << endl << "Delta Ray" << endl;
	155	}
	156	} else if (EVTFLG.lpairp == 1) {
	157	mech = kPPair;
	158	if (debug) cout << endl << "Pair Production" << endl;
	159	} else if (EVTFLG.lbrmsp == 1) {
	160	mech = kPBrem;
	161	if (debug) cout << endl << "Bremsstrahlung" << endl;
	162	}
	163
	164	Double_t weight = GENSTK.wei[numsec-1];
	165	Int_t is = 0;
	166	Int_t ntr;
	167	//
	168	// Save particle in VMC stack
	169	cppstack->PushTrack(done, parent, pdg,
	170	px, py, pz, e,
	171	vx, vy, vz, tof,
	172	polx, poly, polz,
	173	mech, ntr, weight, is);
	174	if (debug) {
	175
	176	cout << endl << " !!! stuprf: ntr=" << ntr << " pdg " << pdg << " parent=" << parent
	177	<< " parent_pdg="<< fluka->PDGFromId(TRACKR.jtrack) << " numsec "
	178	<< numsec << " npprmr " << npprmr << " icode=" << fluka->GetIcode()
	179	<< "kin. energy [keV] = " << GENSTK.tki[numsec-1] * 1.e6 << endl << endl;
	180
	181	}
	182
	183
	184	//
	185	// Save current track number
	186	FLKSTK.ispark[FLKSTK.npflka][mkbmx2-1] = ntr;
	187	FLKSTK.ispark[FLKSTK.npflka][mkbmx2-2] = 0;
	188	} // end of if (numsec > npprmr)
	189	// else {
	190	// if(debug) {
	191	// cout << endl << " !!! stuprf: skipping pushtrack track=" << TRACKR.ispusr[mkbmx2-1]
	192	// << " pdg " << fluka->PDGFromId(TRACKR.jtrack) << " numsec=" << numsec<< " npprmr=" << npprmr
	193	// << " GENSTK pdg=" << fluka->PDGFromId(GENSTK.kpart[numsec-1]) << endl;
	194	// }
	195	// }
	196	} // end of stuprf
	197
	198	} // end of extern "C"
	199