From 82a3f706450e81c50b31645c33f923272a300324 Mon Sep 17 00:00:00 2001 From: morsch Date: Tue, 16 Oct 2007 16:20:59 +0000 Subject: [PATCH] Correct energy of interacting particle after BREMS, DELTA, COMPTON --- TFluka/TFluka.cxx | 51 ++++++++++++++++++++++++++++++++++++++++------- TFluka/stupre.cxx | 25 ++++++++++++++++------- 2 files changed, 62 insertions(+), 14 deletions(-) diff --git a/TFluka/TFluka.cxx b/TFluka/TFluka.cxx index 5f434c6d145..ddf3e43851b 100644 --- a/TFluka/TFluka.cxx +++ b/TFluka/TFluka.cxx @@ -150,6 +150,7 @@ TFluka::TFluka() // // Default constructor // + for (Int_t i = 0; i < 4; i++) fPint[i] = 0.; } //______________________________________________________________________________ @@ -185,6 +186,8 @@ TFluka::TFluka(const char *title, Int_t verbosity, Bool_t isRootGeometrySupporte fUserScore(new TObjArray(100)) { // create geometry interface + for (Int_t i = 0; i < 4; i++) fPint[i] = 0.; + if (fVerbosityLevel >=3) cout << "<== TFluka::TFluka(" << title << ") constructor called." << endl; SetCoreInputFileName(); @@ -1425,9 +1428,10 @@ void TFluka::TrackMomentum(TLorentzVector& momentum) const FlukaCallerCode_t caller = GetCaller(); FlukaProcessCode_t icode = GetIcode(); - if (caller != kEEDRAW && - caller != kMGResumedTrack && - caller != kSODRAW && + if (caller != kEEDRAW && + caller != kMGResumedTrack && + caller != kSODRAW && + caller != kUSDRAW && (caller != kENDRAW || (icode != kEMFSCOstopping1 && icode != kEMFSCOstopping2))) { if (TRACKR.ptrack >= 0) { momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck); @@ -1466,6 +1470,19 @@ void TFluka::TrackMomentum(TLorentzVector& momentum) const momentum.SetPy(p * FLKSTK.tyflk[ist]); momentum.SetPz(p * FLKSTK.tzflk[ist]); momentum.SetE(e); + } else if (caller == kUSDRAW) { + if (icode == 208 || icode == 210 || icode == 212 || icode == 219) { + momentum.SetPx(fPint[0]); + momentum.SetPy(fPint[1]); + momentum.SetPz(fPint[2]); + momentum.SetE(fPint[3]); + } else { + Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack - ParticleMassFPC(TRACKR.jtrack) * ParticleMassFPC(TRACKR.jtrack)); + momentum.SetPx(p*TRACKR.cxtrck); + momentum.SetPy(p*TRACKR.cytrck); + momentum.SetPz(p*TRACKR.cztrck); + momentum.SetE(TRACKR.etrack); + } } else Warning("TrackMomentum","momentum not available"); @@ -1487,6 +1504,7 @@ void TFluka::TrackMomentum(Double_t& px, Double_t& py, Double_t& pz, Double_t& e if (caller != kEEDRAW && caller != kMGResumedTrack && caller != kSODRAW && + caller != kUSDRAW && (caller != kENDRAW || (icode != kEMFSCOstopping1 && icode != kEMFSCOstopping2))) { if (TRACKR.ptrack >= 0) { px = TRACKR.ptrack*TRACKR.cxtrck; @@ -1523,6 +1541,19 @@ void TFluka::TrackMomentum(Double_t& px, Double_t& py, Double_t& pz, Double_t& e px = p * FLKSTK.txflk[ist]; py = p * FLKSTK.tyflk[ist]; pz = p * FLKSTK.tzflk[ist]; + } else if (caller == kUSDRAW) { + if (icode == 208 || icode == 210 || icode == 212 || icode == 219) { + px = fPint[0]; + py = fPint[1]; + pz = fPint[2]; + e = fPint[3]; + } else { + Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack - ParticleMassFPC(TRACKR.jtrack) * ParticleMassFPC(TRACKR.jtrack)); + px = p*TRACKR.cxtrck; + py = p*TRACKR.cytrck; + pz = p*TRACKR.cztrck; + e = TRACKR.etrack; + } } else Warning("TrackMomentum","momentum not available"); @@ -1558,8 +1589,10 @@ Double_t TFluka::TrackLength() const return TRACKR.cmtrck; else if (caller == kMGResumedTrack) return TRACKR.spausr[8]; + else if (caller == kSODRAW) + return 0.0; else { - Warning("TrackLength", "track length not available"); + Warning("TrackLength", "track length not available for caller %5d \n", caller); return 0.0; } } @@ -1704,7 +1737,7 @@ Double_t TFluka::TrackMass() const if (caller != kEEDRAW && caller != kSODRAW) return PAPROP.am[CorrectFlukaId()+6]; else if (caller == kSODRAW) { - Int_t ifl = PDGFromId(FLKSTK.iloflk[FLKSTK.npflka]); + Int_t ifl = FLKSTK.iloflk[FLKSTK.npflka]; return PAPROP.am[ifl + 6]; } else @@ -1973,12 +2006,16 @@ Int_t TFluka::StepProcesses(TArrayI &proc) const proc.Set(1); TMCProcess iproc; - if (caller == kBXEntering || caller == kBXExiting || caller == kENDRAW) { + if (caller == kBXEntering || caller == kBXExiting || caller == kEEDRAW) { iproc = kPTransportation; } else { switch (icode) { case kEMFSCO: - iproc = kPEnergyLoss; + if (Edep() > 0.) { + iproc = kPEnergyLoss; + } else { + iproc = kPTransportation; + } break; case kKASKADtimekill: case kEMFSCOtimekill: diff --git a/TFluka/stupre.cxx b/TFluka/stupre.cxx index b1f54559fec..3e62829d7d1 100644 --- a/TFluka/stupre.cxx +++ b/TFluka/stupre.cxx @@ -65,6 +65,8 @@ void stupre() Int_t verbosityLevel = fluka->GetVerbosityLevel(); Bool_t debug = (verbosityLevel>=3)?kTRUE:kFALSE; + debug = 100; + fluka->SetTrackIsNew(kTRUE); // Get the stack produced from the generator @@ -155,7 +157,9 @@ void stupre() 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; - } + } else { + fluka->SetPint(px, py, pz, e); + } } // end of lcmptn //* Bremsstrahlung: true secondary only if charge = 0 (photon) @@ -168,7 +172,10 @@ void stupre() 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; - } + } else { + fluka->SetPint(px, py, pz, e); + } + } // end of lbrmsp //* Delta ray: If Bhabha, true secondary only if negative (electron) @@ -177,12 +184,14 @@ void stupre() 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); + 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 + if (debug) cout << endl << " !!! stupre (BHABA) : ntr=" << ntr << " pdg " << pdg << " parent=" << parent << endl; + } else { + fluka->SetPint(px, py, pz, e); + } } // lbhabh == 1 //* Delta ray: Otherwise Moller: true secondary is the electron with @@ -195,7 +204,9 @@ void stupre() 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 + } else { + fluka->SetPint(px, py, pz, e); + } } // end of ldltry } // end of loop -- 2.39.3