#include <G4Absorber.hh>
Public Member Functions | |
G4Absorber (G4double cutOnP) | |
~G4Absorber () | |
G4bool | WillBeAbsorbed (const G4KineticTrack &kt) |
G4bool | Absorb (G4KineticTrack &kt, G4KineticTrackVector &tgt) |
G4KineticTrackVector * | GetAbsorbers () |
G4KineticTrackVector * | GetProducts () |
G4bool | FindAbsorbers (G4KineticTrack &kt, G4KineticTrackVector &tgt) |
G4bool | FindProducts (G4KineticTrack &kt) |
Definition at line 34 of file G4Absorber.hh.
G4Absorber::G4Absorber | ( | G4double | cutOnP | ) |
Definition at line 38 of file G4Absorber.cc.
00039 { 00040 theCutOnP = cutOnP; 00041 theAbsorbers = new G4KineticTrackVector; 00042 theProducts = new G4KineticTrackVector; 00043 }
G4Absorber::~G4Absorber | ( | ) |
G4bool G4Absorber::Absorb | ( | G4KineticTrack & | kt, | |
G4KineticTrackVector & | tgt | |||
) |
Definition at line 72 of file G4Absorber.cc.
References FindAbsorbers(), and FindProducts().
00073 { 00074 if(!FindAbsorbers(kt, tgt)) 00075 return false; 00076 return FindProducts(kt); 00077 }
G4bool G4Absorber::FindAbsorbers | ( | G4KineticTrack & | kt, | |
G4KineticTrackVector & | tgt | |||
) |
Definition at line 80 of file G4Absorber.cc.
References DBL_MAX, G4KineticTrack::GetDefinition(), G4ParticleDefinition::GetPDGCharge(), and G4KineticTrack::GetPosition().
Referenced by Absorb().
00082 { 00083 // Find a closest ( in space) pair of Nucleons capable to absorb pi+/pi- 00084 // pi+ can be absorbed on np or nn resulting in pp or np 00085 // pi- can be absorbed on np or pp resulting in nn or np 00086 00087 // @GF: FindAbsorbers is unused, logic is seriously wrong 00088 00089 G4KineticTrack * kt1 = NULL; 00090 G4KineticTrack * kt2 = NULL; 00091 G4double dist1 = DBL_MAX; // dist to closest nucleon 00092 G4double dist2 = DBL_MAX; // dist to next close 00093 G4double charge1 = 0; 00094 // G4double charge2 = 0; // charge2 is only assigned to, never used 00095 G4double charge0 = kt.GetDefinition()->GetPDGCharge(); 00096 G4ThreeVector pos = kt.GetPosition(); 00097 00098 std::vector<G4KineticTrack *>::iterator iter; 00099 for(iter = tgt.begin(); iter != tgt.end(); ++iter) 00100 { 00101 G4KineticTrack * curr = *iter; 00102 G4double dist = (pos-curr->GetPosition()).mag(); 00103 if(dist >= dist2) 00104 continue; 00105 if(dist < dist1) 00106 { 00107 if(dist1 == DBL_MAX) // accept 1st as a candidate, 00108 { 00109 kt1 = curr; 00110 charge1 = kt1->GetDefinition()->GetPDGCharge(); 00111 dist1 = dist; 00112 continue; 00113 } 00114 if(dist2 == DBL_MAX) // accept the candidate and shift kt1 to kt2 00115 { // @GF: should'nt we check if compatible? 00116 kt2 = kt1; 00117 // charge2 = charge1; 00118 dist2 = dist1; 00119 kt1 = curr; 00120 charge1 = kt1->GetDefinition()->GetPDGCharge(); 00121 dist1 = dist; 00122 continue; 00123 } 00124 // test the compatibility with charge conservation for new config 00125 G4double charge = curr->GetDefinition()->GetPDGCharge(); 00126 if((charge0+charge1+charge < 0.) || //test config (curr,kt1) 00127 (charge0+charge1+charge) > 2*eplus) 00128 { // incompatible: change kt1 with curr. 00129 kt1 = curr; 00130 charge1 = charge; 00131 dist1 = dist; 00132 } 00133 else 00134 { // compatible: change kt1 with curr and kt2 with kt1 00135 kt2 = kt1; 00136 // charge2 = charge1; 00137 dist2 = dist1; 00138 kt1 = curr; 00139 charge1 = charge; 00140 dist1 = dist; 00141 } 00142 continue; 00143 } 00144 // here if dist1 < dist < dist2 00145 if(dist2 == DBL_MAX) // accept the candidate 00146 { 00147 kt2 = curr; 00148 // charge2 = kt2->GetDefinition()->GetPDGCharge(); 00149 dist2 = dist; 00150 continue; 00151 } 00152 // test the compatibility with charge conservation 00153 G4double charge = curr->GetDefinition()->GetPDGCharge(); 00154 if((charge0+charge1+charge < 0.) || 00155 (charge0+charge1+charge) > 2*eplus) 00156 continue; // incomatible: do nothing 00157 // compatible: change kt2 with curr 00158 kt2 = curr; 00159 // charge2 = charge; 00160 dist2 = dist; 00161 } 00162 00163 theAbsorbers->clear(); // do not delete tracks in theAbsorbers vector! 00164 if((kt1 == NULL) || (kt2 == NULL)) 00165 return false; 00166 00167 theAbsorbers->push_back(kt1); 00168 theAbsorbers->push_back(kt2); 00169 return true; 00170 }
G4bool G4Absorber::FindProducts | ( | G4KineticTrack & | kt | ) |
Definition at line 174 of file G4Absorber.cc.
References G4KineticTrack::Get4Momentum(), G4KineticTrack::GetDefinition(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), G4KineticTrack::GetPosition(), G4Neutron::Neutron(), and G4Proton::Proton().
Referenced by Absorb().
00175 { 00176 // Choose the products type 00177 G4ParticleDefinition * prod1; 00178 G4ParticleDefinition * prod2; 00179 G4KineticTrack * abs1 = (*theAbsorbers)[0]; 00180 G4KineticTrack * abs2 = (*theAbsorbers)[1]; 00181 00182 G4double charge = kt.GetDefinition()->GetPDGCharge(); 00183 if(charge == eplus) 00184 { // a neutron become proton 00185 prod1 = G4Proton::Proton(); 00186 if(abs1->GetDefinition() == G4Neutron::Neutron()) 00187 prod2 = abs2->GetDefinition(); 00188 else 00189 prod2 = G4Proton::Proton(); 00190 } 00191 else if(charge == -eplus) 00192 { // a proton become neutron 00193 prod1 = G4Neutron::Neutron(); 00194 if(abs1->GetDefinition() == G4Proton::Proton()) 00195 prod2 = abs2->GetDefinition(); 00196 else 00197 prod2 = G4Neutron::Neutron(); 00198 } 00199 else // charge = 0: leave particle types unchenged 00200 { 00201 prod1 = abs1->GetDefinition(); 00202 prod2 = abs2->GetDefinition(); 00203 } 00204 00205 // Translate to the CMS frame 00206 G4LorentzVector momLab = kt.Get4Momentum()+abs1->Get4Momentum()+ 00207 abs2->Get4Momentum(); 00208 G4LorentzRotation toCMSFrame((-1)*momLab.boostVector()); 00209 G4LorentzRotation toLabFrame(momLab.boostVector()); 00210 G4LorentzVector momCMS = toCMSFrame*momLab; 00211 00212 // Evaluate the final momentum of products 00213 G4double ms1 = prod1->GetPDGMass(); 00214 G4double ms2 = prod2->GetPDGMass(); 00215 G4double e0 = momCMS.e(); 00216 G4double squareP = (e0*e0*e0*e0-2*e0*e0*(ms1*ms1+ms2*ms2)+ 00217 (ms2*ms2-ms1*ms1)*(ms2*ms2-ms1*ms1))/(4*e0*e0); 00218 // if(squareP < 0) // should never happen 00219 // squareP = 0; 00220 G4ThreeVector mom1CMS = GetRandomDirection(); 00221 mom1CMS = std::sqrt(squareP)*mom1CMS; 00222 G4LorentzVector final4Mom1CMS(mom1CMS, std::sqrt(squareP+ms1*ms1)); 00223 G4LorentzVector final4Mom2CMS((-1)*mom1CMS, std::sqrt(squareP+ms2*ms2)); 00224 00225 // Go back to the lab frame 00226 G4LorentzVector mom1 = toLabFrame*final4Mom1CMS; 00227 G4LorentzVector mom2 = toLabFrame*final4Mom2CMS; 00228 00229 // ------ debug 00230 /* 00231 G4LorentzVector temp = mom1+mom2; 00232 00233 cout << (1/MeV)*momLab.x() << " " << (1/MeV)*momLab.y() << " " 00234 << (1/MeV)*momLab.z() << " " << (1/MeV)*momLab.t() << " " 00235 << (1/MeV)*momLab.vect().mag() << " " << (1/MeV)*momLab.mag() << " " 00236 << (1/MeV)*temp.x() << " " << (1/MeV)*temp.y() << " " 00237 << (1/MeV)*temp.z() << " " << (1/MeV)*temp.t() << " " 00238 << (1/MeV)*temp.vect().mag() << " " << (1/MeV)*temp.mag() << " " 00239 << (1/MeV)*std::sqrt(squareP) << endl; 00240 00241 */ 00242 // ------ end debug 00243 00244 // Build two new kinetic tracks and add to products 00245 G4KineticTrack * kt1 = new G4KineticTrack(prod1, 0., abs1->GetPosition(), 00246 mom1); 00247 G4KineticTrack * kt2 = new G4KineticTrack(prod2, 0., abs2->GetPosition(), 00248 mom2); 00249 // ------ debug 00250 /* 00251 G4LorentzVector initialMom1 = abs1->Get4Momentum(); 00252 G4LorentzVector initialMom2 = abs2->Get4Momentum(); 00253 G4LorentzVector pion4MomCMS = toCMSFrame*kt.Get4Momentum(); 00254 cout << (1/MeV)*initialMom1.x() << " " << (1/MeV)*initialMom1.y() << " " 00255 << (1/MeV)*initialMom1.z() << " " << (1/MeV)*initialMom1.e() << " " 00256 << (1/MeV)*initialMom1.vect().mag() << " " 00257 << (1/MeV)*initialMom2.x() << " " << (1/MeV)*initialMom2.y() << " " 00258 << (1/MeV)*initialMom2.z() << " " << (1/MeV)*initialMom2.e() << " " 00259 << (1/MeV)*initialMom2.vect().mag() << " " 00260 << (1/MeV)*mom1.x() << " " << (1/MeV)*mom1.y() << " " 00261 << (1/MeV)*mom1.z() << " " << (1/MeV)*mom1.e() << " " 00262 << (1/MeV)*mom1.vect().mag() << " " 00263 << (1/MeV)*mom2.x() << " " << (1/MeV)*mom2.y() << " " 00264 << (1/MeV)*mom2.z() << " " << (1/MeV)*mom2.e() << " " 00265 << (1/MeV)*mom2.vect().mag() << " " 00266 << (1/MeV)*pion4MomCMS.x() << " " << (1/MeV)*pion4MomCMS.y() << " " 00267 << (1/MeV)*pion4MomCMS.z() << " " << (1/MeV)*pion4MomCMS.e() << " " 00268 << (1/MeV)*pion4MomCMS.vect().mag() << " " 00269 << (1/MeV)*final4Mom1CMS.x() << " " << (1/MeV)*final4Mom1CMS.y() << " " 00270 << (1/MeV)*final4Mom1CMS.z() << " " << (1/MeV)*final4Mom1CMS.e() << " " 00271 << (1/MeV)*final4Mom1CMS.vect().mag() << " " 00272 << (1/MeV)*final4Mom2CMS.x() << " " << (1/MeV)*final4Mom2CMS.y() << " " 00273 << (1/MeV)*final4Mom2CMS.z() << " " << (1/MeV)*final4Mom2CMS.e() << " " 00274 << (1/MeV)*final4Mom2CMS.vect().mag() << endl; 00275 */ 00276 // ------ end debug 00277 00278 theProducts->clear(); 00279 theProducts->push_back(kt1); 00280 theProducts->push_back(kt2); 00281 return true; 00282 }
G4KineticTrackVector * G4Absorber::GetAbsorbers | ( | ) | [inline] |
G4KineticTrackVector * G4Absorber::GetProducts | ( | ) | [inline] |
bool G4Absorber::WillBeAbsorbed | ( | const G4KineticTrack & | kt | ) |
Definition at line 53 of file G4Absorber.cc.
References G4KineticTrack::Get4Momentum(), G4KineticTrack::GetActualMass(), G4KineticTrack::GetDefinition(), G4PionMinus::PionMinus(), G4PionPlus::PionPlus(), and G4PionZero::PionZero().
00054 { 00055 // FixMe: actually only for pions 00056 // if(kt.Get4Momentum().vect().mag() < theCutOnP) 00057 // Cut on kinetic Energy... 00058 if (kt.Get4Momentum().e() - kt.GetActualMass() < theCutOnP) 00059 { 00060 if(kt.GetDefinition() == G4PionPlus::PionPlus() || 00061 kt.GetDefinition() == G4PionZero::PionZero() || 00062 kt.GetDefinition() == G4PionMinus::PionMinus()) 00063 { 00064 return true; 00065 } 00066 } 00067 return false; 00068 }