X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TFluka%2Fstupre.cxx;h=937c2368ba718d1b5203bec9ac90b9a73de5d455;hb=56982dd1cb7231735a60927784e42c51c115c4e5;hp=b1ccfc4078612765dbeea323d03e29c4d5664577;hpb=bcf1cd9c85dc92d768c09ac1a7e3acaf83284332;p=u%2Fmrichter%2FAliRoot.git diff --git a/TFluka/stupre.cxx b/TFluka/stupre.cxx index b1ccfc40786..937c2368ba7 100644 --- a/TFluka/stupre.cxx +++ b/TFluka/stupre.cxx @@ -1,5 +1,4 @@ #include -#include "AliRun.h" #ifndef WIN32 # define stupre stupre_ #else @@ -8,17 +7,18 @@ // // 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" @@ -32,39 +32,41 @@ 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 +// 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 +// 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; + 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(); @@ -72,31 +74,51 @@ void stupre() // 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 + +// 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; + 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; + - if (EMFSTK.ichemf[kp] == -1) flukaid = 3; - else if (EMFSTK.ichemf[kp] == 0) flukaid = 7; - else if (EMFSTK.ichemf[kp] == 0) flukaid = 4; - Int_t pdg = fluka->PDGFromId(flukaid); - Double_t e = EMFSTK.etemf[kp] * emvgev; + 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.u[kp]; - Double_t pz = p * EMFSTK.v[kp]; - Double_t py = p * EMFSTK.w[kp]; + 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.x[kp]; - Double_t vy = EMFSTK.y[kp]; - Double_t vz = EMFSTK.z[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; @@ -114,9 +136,10 @@ void stupre() 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; + 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-) @@ -129,6 +152,7 @@ void stupre() 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 @@ -141,6 +165,7 @@ void stupre() 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 @@ -153,6 +178,7 @@ void stupre() 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 @@ -166,12 +192,11 @@ void stupre() 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"