00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032 #include "G4NeutronHPEnAngCorrelation.hh"
00033 #include "G4LorentzRotation.hh"
00034 #include "G4LorentzVector.hh"
00035 #include "G4RotationMatrix.hh"
00036
00037 G4ReactionProduct * G4NeutronHPEnAngCorrelation::SampleOne(G4double anEnergy)
00038 {
00039 G4ReactionProduct * result = new G4ReactionProduct;
00040
00041
00042 if(nProducts!=1) throw G4HadronicException(__FILE__, __LINE__, "More than one product in SampleOne");
00043
00044
00045 G4ReactionProductVector * temp=0;
00046 G4int i=0;
00047 while(temp == 0) temp = theProducts[i++].Sample(anEnergy);
00048
00049
00050 if(temp->size()!=1) throw G4HadronicException(__FILE__, __LINE__, "SampleOne: Yield not correct");
00051
00052
00053 result = temp->operator[](0);
00054
00055
00056 delete temp;
00057
00058
00059 return result;
00060 }
00061
00062 G4ReactionProductVector * G4NeutronHPEnAngCorrelation::Sample(G4double anEnergy)
00063 {
00064 G4ReactionProductVector * result = new G4ReactionProductVector;
00065 G4int i;
00066 G4ReactionProductVector * it;
00067 G4ReactionProduct theCMS;
00068 G4LorentzRotation toZ;
00069
00070
00071 if(frameFlag==2||frameFlag==3)
00072 {
00073
00074 G4ThreeVector the3Neutron = theNeutron.GetMomentum();
00075 G4double nEnergy = theNeutron.GetTotalEnergy();
00076 G4ThreeVector the3Target = theTarget.GetMomentum();
00077 G4double tEnergy = theTarget.GetTotalEnergy();
00078 G4double totE = nEnergy+tEnergy;
00079 G4ThreeVector the3CMS = the3Target+the3Neutron;
00080 theCMS.SetMomentum(the3CMS);
00081 G4double cmsMom = std::sqrt(the3CMS*the3CMS);
00082 G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
00083 theCMS.SetMass(sqrts);
00084 theCMS.SetTotalEnergy(totE);
00085 G4ReactionProduct aNeutron;
00086 aNeutron.Lorentz(theNeutron, theCMS);
00087
00088
00089
00090
00091
00092 anEnergy = theNeutron.GetKineticEnergy();
00093
00094 G4LorentzVector Ptmp (aNeutron.GetMomentum(), aNeutron.GetTotalEnergy());
00095
00096 toZ.rotateZ(-1*Ptmp.phi());
00097 toZ.rotateY(-1*Ptmp.theta());
00098 }
00099 theTotalMeanEnergy=0;
00100 G4LorentzRotation toLab(toZ.inverse());
00101 for(i=0; i<nProducts; i++)
00102 {
00103 it = theProducts[i].Sample(anEnergy);
00104 G4double aMeanEnergy = theProducts[i].MeanEnergyOfThisInteraction();
00105 if(aMeanEnergy>0)
00106 {
00107 theTotalMeanEnergy += aMeanEnergy;
00108 }
00109 else
00110 {
00111 theTotalMeanEnergy = anEnergy/nProducts+theProducts[i].GetQValue();
00112 }
00113 if(it!=0)
00114 {
00115 for(unsigned int ii=0; ii<it->size(); ii++)
00116 {
00117 G4LorentzVector pTmp1 (it->operator[](ii)->GetMomentum(),
00118 it->operator[](ii)->GetTotalEnergy());
00119 pTmp1 = toLab*pTmp1;
00120 it->operator[](ii)->SetMomentum(pTmp1.vect());
00121 it->operator[](ii)->SetTotalEnergy(pTmp1.e());
00122 if(frameFlag==1)
00123 {
00124 it->operator[](ii)->Lorentz(*(it->operator[](ii)), -1.*theTarget);
00125 }
00126 else if(frameFlag==2)
00127 {
00128 #ifdef G4_NHP_DEBUG
00129 cout <<"G4NeutronHPEnAngCorrelation: "<<
00130 it->at(ii)->GetTotalEnergy()<<" "<<
00131 it->at(ii)->GetMomentum()<<G4endl;
00132 #endif
00133 it->operator[](ii)->Lorentz(*(it->operator[](ii)), -1.*theCMS);
00134 }
00135
00136 else if(frameFlag==3)
00137 {
00138 if ( theProducts[i].GetMassCode() > 4 )
00139 {
00140
00141 it->operator[](ii)->Lorentz(*(it->operator[](ii)), -1.*theTarget);
00142 }
00143 else
00144 {
00145
00146 it->operator[](ii)->Lorentz(*(it->operator[](ii)), -1.*theCMS);
00147 }
00148 }
00149 else
00150 {
00151 throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPEnAngCorrelation::Sample: The frame of the finalstate is not specified");
00152 }
00153 result->push_back(it->operator[](ii));
00154 }
00155 delete it;
00156 }
00157 }
00158 return result;
00159 }
00160