- Timestamp:
- Apr 6, 2009, 12:30:29 PM (16 years ago)
- Location:
- trunk/source/processes/hadronic/models/parton_string/diffraction
- Files:
-
- 11 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/models/parton_string/diffraction/History
r819 r962 1 $Id: History,v 1. 1 2007/04/24 10:32:59 gunterExp $1 $Id: History,v 1.5 2008/12/09 10:43:47 vuzhinsk Exp $ 2 2 ------------------------------------------------------------------- 3 3 … … 15 15 * Please list in reverse chronological order (last date on top) 16 16 --------------------------------------------------------------- 17 9 December 08, V. Uzhinsky (hadr-string-diff-V09-01-04) 18 - Improvement of delete operators in FTF 19 20 5 December 08, V. Uzhinsky (hadr-string-diff-V09-01-03) 21 - Some objects did not erase in FTFModel desructor. These lead to memory 22 leak. 17 23 18 24 25 2 June 08, G.Folger (hadr-string-diff-V09-01-02) 26 - on branch geant4-09-01-patches_branch, fix compilation warning for unused 27 variables in G4FTFModel.cc 28 29 30 31 ?????????????????????????????? (hadr-string-diff-V09-01-01) 32 33 31 March 2008 V. Uzhinsky Tag : hadr-string-diff-V09-01-00 34 - G4FTFParameters.cc and G4FTFParameters.hh were copied from G4FTFCrossSection 35 corresponding files. 36 - New files - G4ElasticHNScattering have been added. They implement elastic 37 scattering of hadron in intra-nuclear collisions in FTF model. 38 - The corresponding changes have been done in G4FTFModel.cc and 39 G4FTFParticipants.cc 40 41 19 42 24 Apr 2007 Gunter Folger (hadr-string-diff-V08-02-00) 20 43 - merge in change done by ftf development. 21 44 - Created History file. 45 -
trunk/source/processes/hadronic/models/parton_string/diffraction/include/G4DiffractiveExcitation.hh
r819 r962 25 25 // 26 26 // 27 // $Id: G4DiffractiveExcitation.hh,v 1. 1 2007/05/25 06:56:53 vuzhinsk Exp $27 // $Id: G4DiffractiveExcitation.hh,v 1.2 2008/04/25 14:20:13 vuzhinsk Exp $ 28 28 29 29 #ifndef G4DiffractiveExcitation_h … … 42 42 class G4VSplitableHadron; 43 43 class G4ExcitedString; 44 #include "G4FTFParameters.hh" // Uzhi 19.04.08 44 45 #include "G4ThreeVector.hh" 45 46 … … 52 53 virtual ~G4DiffractiveExcitation(); 53 54 54 virtual G4bool ExciteParticipants (G4VSplitableHadron *aPartner, G4VSplitableHadron * bPartner) const; 55 virtual G4bool ExciteParticipants (G4VSplitableHadron *aPartner, 56 G4VSplitableHadron * bPartner, 57 G4FTFParameters *theParameters) const; 58 55 59 virtual G4ExcitedString * String(G4VSplitableHadron * aHadron, G4bool isProjectile) const; 56 57 // void SetPtWidth(G4double aValue) { widthOfPtSquare = aValue*aValue; }58 // void SetExtraMass(G4double aValue) { minExtraMass = aValue; }59 // void SetMinimumMass(G4double aValue) { minmass = aValue; }60 61 60 62 61 private: … … 64 63 G4DiffractiveExcitation(const G4DiffractiveExcitation &right); 65 64 66 // G4double ChooseX(G4double Xmin, G4double Xmax) const;// Uzhi65 G4ThreeVector GaussianPt(G4double AveragePt2, G4double maxPtSquare) const; // Uzhi 67 66 G4double ChooseP(G4double Pmin, G4double Pmax) const; // Uzhi 68 69 // G4ThreeVector GaussianPt(G4double widthSquare, G4double maxPtSquare) const;70 G4ThreeVector GaussianPt(G4double AveragePt2, G4double maxPtSquare) const; // Uzhi71 67 72 68 const G4DiffractiveExcitation & operator=(const G4DiffractiveExcitation &right); … … 74 70 int operator!=(const G4DiffractiveExcitation &right) const; 75 71 76 private:77 // Model Parameters:78 /* // Uzhi79 const G4double widthOfPtSquare; // width^2 of pt for string excitation80 const G4double minExtraMass; // minimum excitation mass81 const G4double minmass; // mean pion transverse mass; used for Xmin82 */ // Uzhi83 72 }; 84 73 -
trunk/source/processes/hadronic/models/parton_string/diffraction/include/G4DiffractiveHHScatterer.hh
r819 r962 25 25 // 26 26 // 27 // $Id: G4DiffractiveHHScatterer.hh,v 1. 3 2006/06/29 20:54:25 gunterExp $27 // $Id: G4DiffractiveHHScatterer.hh,v 1.4 2008/04/25 14:20:13 vuzhinsk Exp $ 28 28 29 29 #ifndef G4DiffractiveHHScatterer_h … … 44 44 class G4KineticTrack; 45 45 #include "G4KineticTrackVector.hh" 46 #include "G4FTFParameters.hh" // Uzhi 21.04.08 46 47 47 48 class G4DiffractiveHHScatterer … … 57 58 const G4DiffractiveExcitation * theExcitation; 58 59 G4LundStringFragmentation * theStringFragmentation; 60 G4FTFParameters *theParameters; // Uzhi 21.04.08 59 61 }; 60 62 -
trunk/source/processes/hadronic/models/parton_string/diffraction/include/G4DiffractiveSplitableHadron.hh
r819 r962 26 26 // 27 27 // $Id: G4DiffractiveSplitableHadron.hh,v 1.4 2006/06/29 20:54:27 gunter Exp $ 28 // GEANT4 tag $Name: $28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 -
trunk/source/processes/hadronic/models/parton_string/diffraction/include/G4FTFModel.hh
r819 r962 25 25 // 26 26 // 27 // $Id: G4FTFModel.hh,v 1. 5 2007/04/24 10:32:59 gunterExp $28 // GEANT4 tag $Name: geant4-09-0 1-patch-02 $27 // $Id: G4FTFModel.hh,v 1.7 2008/04/25 14:20:13 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // Class Description … … 52 52 class G4VSplitableHadron; 53 53 class G4ExcitedString; 54 55 #include "G4FTFParameters.hh" // Uzhi 29.03.08 54 56 #include "G4FTFParticipants.hh" 55 57 56 58 #include "G4ExcitedStringVector.hh" 57 59 #include "G4DiffractiveExcitation.hh" 58 60 #include "G4ElasticHNScattering.hh" 59 61 60 62 class G4FTFModel : public G4VPartonStringModel … … 84 86 85 87 private: 86 88 89 G4ReactionProduct theProjectile; 87 90 G4FTFParticipants theParticipants; 88 G4ReactionProduct theProjectile; 89 91 92 G4FTFParameters *theParameters; // Uzhi 29.03.08 90 93 G4DiffractiveExcitation * theExcitation; 91 94 G4ElasticHNScattering * theElastic; // Uzhi 29.03.08 92 95 93 96 94 97 }; 95 98 99 // ------------------------------------------------------------ 96 100 inline 97 101 G4V3DNucleus * G4FTFModel::GetWoundedNucleus() const -
trunk/source/processes/hadronic/models/parton_string/diffraction/include/G4FTFParticipants.hh
r819 r962 25 25 // 26 26 // 27 // $Id: G4FTFParticipants.hh,v 1. 4 2006/06/29 20:54:32 gunterExp $28 // GEANT4 tag $Name: geant4-09-0 1-patch-02 $27 // $Id: G4FTFParticipants.hh,v 1.5 2008/03/31 15:34:01 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 … … 41 41 42 42 #include "G4VParticipants.hh" 43 #include "G4FTFParameters.hh" // Uzhi 29.03.08 43 44 #include <vector> 44 45 #include "G4Nucleon.hh" … … 60 61 int operator!=(const G4FTFParticipants &right) const; 61 62 62 void BuildInteractions(const G4ReactionProduct &thePrimary); 63 void GetList(const G4ReactionProduct &thePrimary, 64 G4FTFParameters *theParameters); // Uzhi 29.03.08 65 66 void StartLoop(); 63 67 G4bool Next(); 64 68 const G4InteractionContent & GetInteraction() const; 65 66 void StartLoop();67 69 70 std::vector<G4InteractionContent *> theInteractions; 68 71 private: 69 72 70 std::vector<G4InteractionContent *> theInteractions;73 // std::vector<G4InteractionContent *> theInteractions; 71 74 72 75 G4int currentInteraction; … … 95 98 96 99 #endif 97 98 -
trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4DiffractiveExcitation.cc
r819 r962 25 25 // 26 26 // 27 // $Id: G4DiffractiveExcitation.cc,v 1. 1 2007/05/25 06:56:53 vuzhinskExp $27 // $Id: G4DiffractiveExcitation.cc,v 1.7 2008/12/18 13:01:58 gunter Exp $ 28 28 // ------------------------------------------------------------ 29 29 // GEANT 4 class implemetation file … … 31 31 // ---------------- G4DiffractiveExcitation -------------- 32 32 // by Gunter Folger, October 1998. 33 // 34 // 35 // 33 // diffractive Excitation used by strings models 34 // Take a projectile and a target 35 // excite the projectile and target 36 36 // Essential changed by V. Uzhinsky in November - December 2006 37 37 // in order to put it in a correspondence with original FRITIOF 38 38 // model. Variant of FRITIOF with nucleon de-excitation is implemented. 39 39 // Other changes by V.Uzhinsky in May 2007 were introduced to fit 40 // meson-nucleon interactions 40 // meson-nucleon interactions. Additional changes by V. Uzhinsky 41 // were introduced in December 2006. They treat diffraction dissociation 42 // processes more exactly. 41 43 // --------------------------------------------------------------------- 42 44 … … 51 53 #include "G4VSplitableHadron.hh" 52 54 #include "G4ExcitedString.hh" 55 #include "G4FTFParameters.hh" // Uzhi 19.04.08 53 56 //#include "G4ios.hh" 54 55 G4DiffractiveExcitation::G4DiffractiveExcitation() // Uzhi 56 { 57 } 58 57 //#include "UZHI_diffraction.hh" 58 59 G4DiffractiveExcitation::G4DiffractiveExcitation() 60 { 61 } 62 63 // --------------------------------------------------------------------- 59 64 G4bool G4DiffractiveExcitation:: 60 ExciteParticipants(G4VSplitableHadron *projectile, G4VSplitableHadron *target) const 61 { 62 63 G4LorentzVector Pprojectile=projectile->Get4Momentum(); 64 65 // -------------------- Projectile parameters ----------------------------------- 66 G4bool PutOnMassShell=0; 67 68 // With de-excitation 69 // G4double M0projectile=projectile->GetDefinition()->GetPDGMass(); 70 // Without de-excitation 71 G4double M0projectile = Pprojectile.mag(); 72 73 if(M0projectile < projectile->GetDefinition()->GetPDGMass()) 74 { 75 PutOnMassShell=1; 76 M0projectile=projectile->GetDefinition()->GetPDGMass(); 77 } 78 79 G4double Mprojectile2 = M0projectile * M0projectile; 80 81 G4int PDGcode=projectile->GetDefinition()->GetPDGEncoding(); 82 G4int absPDGcode=std::abs(PDGcode); 83 G4double ProjectileDiffCut; 84 G4double ProjectileBackgroundWeight; // Uzhi May 2007 85 86 G4double TargetBackgroundWeight; // Uzhi May 2007 87 88 G4double AveragePt2; 89 90 if( absPDGcode > 1000 ) //------Projectile is baryon -------- 91 { 92 ProjectileDiffCut = 1.1; // GeV 93 ProjectileBackgroundWeight=0.; // Uzhi May 2007 94 AveragePt2 = 0.3; // GeV^2 95 } 96 else if( absPDGcode == 211 || PDGcode == 111) //------Projectile is Pion ----------- 97 { 98 ProjectileDiffCut = 0.6; // GeV Uzhi May 2007 1.0->0.6 99 ProjectileBackgroundWeight=0.9; // Uzhi May 2007 100 AveragePt2 = 0.3; // GeV^2 101 } 102 else if( absPDGcode == 321 || PDGcode == -311) //------Projectile is Kaon ----------- 103 { 104 ProjectileDiffCut = 0.7; // GeV 1.1 105 ProjectileBackgroundWeight=0.5; // Uzhi May 2007 ??? 106 AveragePt2 = 0.3; // GeV^2 107 } 108 else //------Projectile is undefined, 109 //------Nucleon assumed 110 { 111 ProjectileDiffCut = 1.1; // GeV 112 ProjectileBackgroundWeight=0.; // Uzhi May 2007 113 AveragePt2 = 0.3; // GeV^2 114 }; 115 116 ProjectileDiffCut = ProjectileDiffCut * GeV; 117 AveragePt2 = AveragePt2 * GeV*GeV; 118 119 // -------------------- Target parameters ---------------------------------------------- 120 G4LorentzVector Ptarget=target->Get4Momentum(); 121 122 G4double M0target = Ptarget.mag(); 123 124 if(M0target < target->GetDefinition()->GetPDGMass()) 125 { 126 PutOnMassShell=1; 127 M0target=target->GetDefinition()->GetPDGMass(); 128 } 129 130 G4double Mtarget2 = M0target * M0target; //Ptarget.mag2(); 131 // for AA-inter. 132 133 G4double NuclearNucleonDiffCut = 1.1*GeV; 134 TargetBackgroundWeight=0.; // Uzhi May 2007 135 136 G4double ProjectileDiffCut2 = ProjectileDiffCut * ProjectileDiffCut; 137 G4double NuclearNucleonDiffCut2 = NuclearNucleonDiffCut * NuclearNucleonDiffCut; 65 ExciteParticipants(G4VSplitableHadron *projectile, 66 G4VSplitableHadron *target, 67 G4FTFParameters *theParameters) const 68 { 69 G4bool PutOnMassShell=0; 70 71 // -------------------- Projectile parameters ----------------------- 72 73 G4LorentzVector Pprojectile=projectile->Get4Momentum(); 74 // G4double M0projectile=projectile->GetDefinition()->GetPDGMass(); // With de-excitation 75 G4double M0projectile = Pprojectile.mag(); // Without de-excitation 76 /* 77 G4cout<<"ExciteParticipants-------------------"<<G4endl; 78 G4cout<<"Mom "<<Pprojectile<<" mass "<<M0projectile<<G4endl; 79 */ 80 if(M0projectile < projectile->GetDefinition()->GetPDGMass()) 81 { 82 PutOnMassShell=1; 83 M0projectile=projectile->GetDefinition()->GetPDGMass(); 84 } 85 86 G4double M0projectile2 = M0projectile * M0projectile; 87 88 G4int PDGcode=projectile->GetDefinition()->GetPDGEncoding(); 89 G4int absPDGcode=std::abs(PDGcode); 90 91 G4double ProjectileDiffStateMinMass=theParameters->GetProjMinDiffMass(); 92 G4double ProjectileNonDiffStateMinMass=theParameters->GetProjMinNonDiffMass(); 93 G4double ProbProjectileDiffraction=theParameters->GetProbabilityOfProjDiff(); 94 /* 95 G4cout<<ProjectileDiffStateMinMass<<" "<<ProjectileNonDiffStateMinMass<<" "<<ProbProjectileDiffraction<<G4endl; 96 */ 97 // -------------------- Target paraExciteParticipantsmeters ------------------------- 98 G4LorentzVector Ptarget=target->Get4Momentum(); 99 G4double M0target = Ptarget.mag(); 100 101 //G4cout<<"Mom "<<Ptarget<<" mass "<<M0target<<G4endl; 102 103 if(M0target < target->GetDefinition()->GetPDGMass()) 104 { 105 PutOnMassShell=1; 106 M0target=target->GetDefinition()->GetPDGMass(); 107 } 108 109 G4double M0target2 = M0target * M0target; //Ptarget.mag2(); 110 // for AA-inter. 111 G4double TargetDiffStateMinMass=theParameters->GetTarMinDiffMass(); 112 G4double TargetNonDiffStateMinMass=theParameters->GetTarMinNonDiffMass(); 113 G4double ProbTargetDiffraction=theParameters->GetProbabilityOfTarDiff(); 114 /* 115 G4cout<<TargetDiffStateMinMass<<" "<<TargetNonDiffStateMinMass<<" "<<ProbTargetDiffraction<<G4endl; 116 */ 117 G4double AveragePt2=theParameters->GetAveragePt2(); 118 119 // Kinematical properties of the interactions -------------- 120 G4LorentzVector Psum; // 4-momentum in CMS 121 Psum=Pprojectile+Ptarget; 122 G4double S=Psum.mag2(); 123 124 //G4cout<<" sqrt(s) "<<std::sqrt(S)<<G4endl; 125 126 // ------------------------------------------------------------------ 127 128 //ProbProjectileDiffraction=1.; 129 //ProbTargetDiffraction =1.; 130 G4double ProbOfDiffraction=ProbProjectileDiffraction + 131 ProbTargetDiffraction; 132 133 if(ProbOfDiffraction!=0.) 134 { 135 ProbProjectileDiffraction/=ProbOfDiffraction; 136 } 137 else 138 { 139 ProbProjectileDiffraction=0.; 140 } 141 // ProbTargetDiffraction /=ProbOfDiffraction; 142 143 //G4cout<<"ProbOfDiffraction "<<ProbOfDiffraction<<"ProbProjectileDiffraction "<<ProbProjectileDiffraction<<G4endl; // Vova 144 145 G4double ProjectileDiffStateMinMass2 = ProjectileDiffStateMinMass * 146 ProjectileDiffStateMinMass; 147 G4double ProjectileNonDiffStateMinMass2 = ProjectileNonDiffStateMinMass * 148 ProjectileNonDiffStateMinMass; 149 150 G4double TargetDiffStateMinMass2 = TargetDiffStateMinMass * 151 TargetDiffStateMinMass; 152 G4double TargetNonDiffStateMinMass2 = TargetNonDiffStateMinMass * 153 TargetNonDiffStateMinMass; 138 154 139 155 // Transform momenta to cms and then rotate parallel to z axis; 140 156 141 G4LorentzVector Psum; 142 Psum=Pprojectile+Ptarget; 143 144 G4LorentzRotation toCms(-1*Psum.boostVector()); 145 146 G4LorentzVector Ptmp=toCms*Pprojectile; 147 148 if ( Ptmp.pz() <= 0. ) 149 { 157 // G4LorentzVector Psum; 158 // Psum=Pprojectile+Ptarget; 159 160 G4LorentzRotation toCms(-1*Psum.boostVector()); 161 162 G4LorentzVector Ptmp=toCms*Pprojectile; 163 if ( Ptmp.pz() <= 0. ) 164 { 150 165 // "String" moving backwards in CMS, abort collision !! 151 166 //G4cout << " abort Collision!! " << G4endl; 152 return false; 153 } 154 155 toCms.rotateZ(-1*Ptmp.phi()); 156 toCms.rotateY(-1*Ptmp.theta()); 157 158 G4LorentzRotation toLab(toCms.inverse()); 159 160 Pprojectile.transform(toCms); 161 Ptarget.transform(toCms); 162 163 G4double Pt2; 164 G4double ProjMassT2, ProjMassT; 165 G4double TargMassT2, TargMassT; 166 G4double PZcms2, PZcms; 167 G4double PMinusNew, TPlusNew; 168 169 G4double S=Psum.mag2(); 170 G4double SqrtS=std::sqrt(S); 171 172 if(absPDGcode > 1000 && SqrtS < 2200*MeV) 173 {return false;} // The model cannot work for 174 // p+p-interactions 175 // at Plab < 1.3 GeV/c. 176 177 if(( absPDGcode == 211 || PDGcode == 111) && SqrtS < 1600*MeV) 178 {return false;} // The model cannot work for 179 // Pi+p-interactions 180 // at Plab < 1. GeV/c. 181 182 if(( absPDGcode == 321 || PDGcode == -311) && SqrtS < 1600*MeV) 183 {return false;} // The model cannot work for 184 // K+p-interactions 185 // at Plab < ??? GeV/c. ??? 186 187 PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2- 188 2*S*Mprojectile2-2*S*Mtarget2-2*Mprojectile2*Mtarget2)/4./S; 189 if(PZcms2 < 0) 190 {return false;} // It can be in an interaction with off-shell nuclear nucleon 191 192 PZcms = std::sqrt(PZcms2); 193 194 if(PutOnMassShell) 167 return false; 168 } 169 170 toCms.rotateZ(-1*Ptmp.phi()); 171 toCms.rotateY(-1*Ptmp.theta()); 172 173 G4LorentzRotation toLab(toCms.inverse()); 174 175 Pprojectile.transform(toCms); 176 Ptarget.transform(toCms); 177 178 G4double Pt2; 179 G4double ProjMassT2, ProjMassT; 180 G4double TargMassT2, TargMassT; 181 G4double PZcms2, PZcms; 182 G4double PMinusMin, PMinusMax; 183 // G4double PPlusMin , PPlusMax; 184 G4double TPlusMin , TPlusMax; 185 G4double PMinusNew, PPlusNew, TPlusNew, TMinusNew; 186 187 // G4double S=Psum.mag2(); 188 G4double SqrtS=std::sqrt(S); 189 190 if(absPDGcode > 1000 && SqrtS < 2200*MeV) 191 {return false;} // The model cannot work for 192 // p+p-interactions 193 // at Plab < 1.3 GeV/c. 194 195 if(( absPDGcode == 211 || PDGcode == 111) && SqrtS < 1600*MeV) 196 {return false;} // The model cannot work for 197 // Pi+p-interactions 198 // at Plab < 1. GeV/c. 199 200 if(( absPDGcode == 321 || PDGcode == -311) && SqrtS < 1600*MeV) 201 {return false;} // The model cannot work for 202 // K+p-interactions 203 // at Plab < ??? GeV/c. ??? 204 205 PZcms2=(S*S+M0projectile2*M0projectile2+M0target2*M0target2- 206 2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2) 207 /4./S; 208 209 if(PZcms2 < 0) 210 {return false;} // It can be in an interaction with off-shell nuclear nucleon 211 212 PZcms = std::sqrt(PZcms2); 213 214 if(PutOnMassShell) 215 { 216 if(Pprojectile.z() > 0.) 217 { 218 Pprojectile.setPz( PZcms); 219 Ptarget.setPz( -PZcms); 220 } 221 else 222 { 223 Pprojectile.setPz(-PZcms); 224 Ptarget.setPz( PZcms); 225 }; 226 227 Pprojectile.setE(std::sqrt(M0projectile2 + 228 Pprojectile.x()*Pprojectile.x()+ 229 Pprojectile.y()*Pprojectile.y()+ 230 PZcms2)); 231 Ptarget.setE(std::sqrt(M0target2 + 232 Ptarget.x()*Ptarget.x()+ 233 Ptarget.y()*Ptarget.y()+ 234 PZcms2)); 235 } 236 237 G4double maxPtSquare; // = PZcms2; 238 /* 239 G4cout << "Pprojectile aft boost : " << Pprojectile <<" "<<Pprojectile.mag()<< G4endl; 240 G4cout << "Ptarget aft boost : " << Ptarget <<" "<<Ptarget.mag()<< G4endl; 241 G4cout << "cms aft boost : " << (Pprojectile+ Ptarget) << G4endl; 242 G4cout << " Projectile Xplus / Xminus : " << 243 Pprojectile.plus() << " / " << Pprojectile.minus() << G4endl; 244 G4cout << " Target Xplus / Xminus : " << Ptarget.plus() << " / " << Ptarget.minus() << G4endl; 245 G4cout<<"maxPtSquare "<<maxPtSquare<<G4endl; 246 */ 247 G4LorentzVector Qmomentum; 248 G4double Qminus, Qplus; 249 250 G4int whilecount=0; 251 // Choose a process 252 253 if(G4UniformRand() < ProbOfDiffraction) 254 { 255 if(G4UniformRand() < ProbProjectileDiffraction) 256 { //-------- projectile diffraction --------------- 257 //G4cout<<" Projectile diffraction"<<G4endl; 258 //Uzhi_projectilediffraction++; 259 do { 260 // Generate pt 261 // if (whilecount++ >= 500 && (whilecount%100)==0) 262 // G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping" 263 // << ", loop count/ maxPtSquare : " 264 // << whilecount << " / " << maxPtSquare << G4endl; 265 if (whilecount > 1000 ) 195 266 { 196 if(Pprojectile.z() > 0.) 197 { 198 Pprojectile.setPz( PZcms); 199 Ptarget.setPz( -PZcms); 200 } 201 else 202 { 203 Pprojectile.setPz(-PZcms); 204 Ptarget.setPz( PZcms); 205 }; 206 207 Pprojectile.setE(std::sqrt(Mprojectile2+ 208 Pprojectile.x()*Pprojectile.x()+ 209 Pprojectile.y()*Pprojectile.y()+ 210 PZcms2)); 211 Ptarget.setE(std::sqrt( Mtarget2 + 212 Ptarget.x()*Ptarget.x()+ 213 Ptarget.y()*Ptarget.y()+ 214 PZcms2)); 215 } 216 217 G4double maxPtSquare = PZcms2; 218 219 //G4cout << "Pprojectile aft boost : " << Pprojectile << G4endl; 220 //G4cout << "Ptarget aft boost : " << Ptarget << G4endl; 221 // G4cout << "cms aft boost : " << (Pprojectile+ Ptarget) << G4endl; 222 223 // G4cout << " Projectile Xplus / Xminus : " << 224 // Pprojectile.plus() << " / " << Pprojectile.minus() << G4endl; 225 // G4cout << " Target Xplus / Xminus : " << 226 // Ptarget.plus() << " / " << Ptarget.minus() << G4endl; 227 228 G4LorentzVector Qmomentum; 229 G4double Qminus, Qplus; 230 231 G4int whilecount=0; 232 do { 233 // Generate pt 234 235 if (whilecount++ >= 500 && (whilecount%100)==0) 267 Qmomentum=G4LorentzVector(0.,0.,0.,0.); 268 return false; // Ignore this interaction 269 }; 270 // --------------- Check that the interaction is possible ----------- 271 ProjMassT2=ProjectileDiffStateMinMass2; 272 ProjMassT =ProjectileDiffStateMinMass; 273 274 TargMassT2=M0target2; 275 TargMassT =M0target; 276 277 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2- 278 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) 279 /4./S; 280 //G4cout<<" Pt2 Mpt Mtt Pz2 "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl; 281 282 if(PZcms2 < 0 ) 283 { 284 /* 285 G4cout<<"whilecount "<<whilecount<<" "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl; 286 G4int Uzhi; G4cin>>Uzhi; 287 */ 288 return false; 289 }; 290 maxPtSquare=PZcms2; 291 292 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0); 293 Pt2=G4ThreeVector(Qmomentum.vect()).mag2(); 294 295 ProjMassT2=ProjectileDiffStateMinMass2+Pt2; 296 ProjMassT =std::sqrt(ProjMassT2); 297 298 TargMassT2=M0target2+Pt2; 299 TargMassT =std::sqrt(TargMassT2); 300 301 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2- 302 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) 303 /4./S; 304 //G4cout<<" Pt2 Mpt Mtt Pz2 "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl; 305 306 // if(PZcms2 < 0 ) {PZcms2=0;}; 307 if(PZcms2 < 0 ) continue; 308 PZcms =std::sqrt(PZcms2); 309 310 PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms; 311 PMinusMax=SqrtS-TargMassT; 312 //G4cout<<" SqrtS P+mim max "<<SqrtS<<" "<<PMinusMin<<" "<<PMinusMax<<G4endl; 313 314 PMinusNew=ChooseP(PMinusMin, PMinusMax); 315 // PMinusNew=1./sqrt(1./PMinusMin-G4UniformRand()*(1./PMinusMin-1./PMinusMax)); 316 317 TMinusNew=SqrtS-PMinusNew; 318 Qminus=Ptarget.minus()-TMinusNew; 319 TPlusNew=TargMassT2/TMinusNew; 320 Qplus=Ptarget.plus()-TPlusNew; 321 322 Qmomentum.setPz( (Qplus-Qminus)/2 ); 323 Qmomentum.setE( (Qplus+Qminus)/2 ); 324 } while ( 325 ((Pprojectile+Qmomentum).mag2() < ProjectileDiffStateMinMass2) || //No without excitation 326 ((Ptarget -Qmomentum).mag2() < M0target2 )); 327 } 328 else 329 { // -------------- Target diffraction ---------------- 330 //G4cout<<" Target difraction"<<G4endl; 331 //Uzhi_targetdiffraction++; 332 do { 333 // Generate pt 334 // if (whilecount++ >= 500 && (whilecount%100)==0) 236 335 // G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping" 237 // << ", loop count/ maxPtSquare : " 336 // << ", loop count/ maxPtSquare : " 238 337 // << whilecount << " / " << maxPtSquare << G4endl; 239 if (whilecount > 1000 ) 240 { 241 Qmomentum=G4LorentzVector(0.,0.,0.,0.); 242 return false; // Ignore this interaction 243 } 244 245 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0); 246 247 //G4cout << "generated Pt " << Qmomentum << G4endl; 248 //G4cout << "Pprojectile with pt : " << Pprojectile+Qmomentum << G4endl; 249 //G4cout << "Ptarget with pt : " << Ptarget-Qmomentum << G4endl; 250 251 // Momentum transfer 252 /* // Uzhi 253 G4double Xmin = minmass / ( Pprojectile.e() + Ptarget.e() ); 254 G4double Xmax=1.; 255 G4double Xplus =ChooseX(Xmin,Xmax); 256 G4double Xminus=ChooseX(Xmin,Xmax); 257 258 // G4cout << " X-plus " << Xplus << G4endl; 259 // G4cout << " X-minus " << Xminus << G4endl; 260 261 G4double pt2=G4ThreeVector(Qmomentum.vect()).mag2(); 262 G4double Qplus =-1 * pt2 / Xminus/Ptarget.minus(); 263 G4double Qminus= pt2 / Xplus /Pprojectile.plus(); 264 */ // Uzhi * 265 266 Pt2=G4ThreeVector(Qmomentum.vect()).mag2(); 267 ProjMassT2=Mprojectile2+Pt2; 268 ProjMassT =std::sqrt(ProjMassT2); 269 270 TargMassT2=Mtarget2+Pt2; 271 TargMassT =std::sqrt(TargMassT2); 272 273 PZcms2=(S*S+ProjMassT2*ProjMassT2+ 274 TargMassT2*TargMassT2- 275 2.*S*ProjMassT2-2.*S*TargMassT2- 276 2.*ProjMassT2*TargMassT2)/4./S; 277 if(PZcms2 < 0 ) {PZcms2=0;}; 278 PZcms =std::sqrt(PZcms2); 279 280 G4double PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms; 281 G4double PMinusMax=SqrtS-TargMassT; 282 283 if(G4UniformRand() < ProjectileBackgroundWeight) // Uzhi May 2007 284 { 285 PMinusNew=PMinusMax-(PMinusMax-PMinusMin)*G4UniformRand(); // Uzhi May 2007 286 } // Uzhi May 2007 287 else // Uzhi May 2007 288 { // Uzhi May 2007 289 PMinusNew=ChooseP(PMinusMin, PMinusMax); // Uzhi May 2007 290 }; // Uzhi May 2007 291 // PMinusNew=ChooseP(PMinusMin, PMinusMax); // Uzhi May 2007 292 Qminus=PMinusNew-Pprojectile.minus(); 293 294 G4double TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms; 295 G4double TPlusMax=SqrtS-PMinusNew; // Uzhi May 2007 296 // G4double TPlusMax=SqrtS-ProjMassT; // Uzhi May 2007 297 298 if(G4UniformRand() < TargetBackgroundWeight) // Uzhi May 2007 299 { 300 TPlusNew=TPlusMax-(TPlusMax-TPlusMin)*G4UniformRand(); // Uzhi May 2007 301 } // Uzhi May 2007 302 else // Uzhi May 2007 303 { // Uzhi May 2007 304 TPlusNew=ChooseP(TPlusMin, TPlusMax); // Uzhi May 2007 305 }; // Uzhi May 2007 306 307 // TPlusNew=ChooseP(TPlusMin, TPlusMax); // Uzhi May 2007 308 Qplus=-(TPlusNew-Ptarget.plus()); 309 310 Qmomentum.setPz( (Qplus-Qminus)/2 ); 311 Qmomentum.setE( (Qplus+Qminus)/2 ); 312 313 //G4cout << "Qplus / Qminus " << Qplus << " / " << Qminus<<G4endl; 314 // G4cout << "pt2" << pt2 << G4endl; 315 // G4cout << "Qmomentum " << Qmomentum << G4endl; 316 // G4cout << " Masses (P/T) : " << (Pprojectile+Qmomentum).mag() << 317 // " / " << (Ptarget-Qmomentum).mag() << G4endl; 318 319 } while ( 320 ( (Pprojectile+Qmomentum).mag2() < Mprojectile2 || // Uzhi No without excitation 321 (Ptarget -Qmomentum).mag2() < Mtarget2 ) || // Uzhi 322 ( (Pprojectile+Qmomentum).mag2() < ProjectileDiffCut2 && // Uzhi No double Diffraction 323 (Ptarget -Qmomentum).mag2() < NuclearNucleonDiffCut2) );// Uzhi 324 325 if((Ptarget-Qmomentum).mag2() < NuclearNucleonDiffCut2) // Uzhi 326 //Projectile diffraction 338 if (whilecount > 1000 ) 327 339 { 328 G4double TMinusNew=SqrtS-PMinusNew; 329 Qminus=Ptarget.minus()-TMinusNew; 330 TPlusNew=TargMassT2/TMinusNew; 331 Qplus=Ptarget.plus()-TPlusNew; 332 333 Qmomentum.setPz( (Qplus-Qminus)/2 ); 334 Qmomentum.setE( (Qplus+Qminus)/2 ); 335 } 336 else if((Pprojectile+Qmomentum).mag2() < ProjectileDiffCut2) // Uzhi 337 //Target diffraction 340 Qmomentum=G4LorentzVector(0.,0.,0.,0.); 341 return false; // Ignore this interaction 342 }; 343 // --------------- Check that the interaction is possible ----------- 344 ProjMassT2=M0projectile2; 345 ProjMassT =M0projectile; 346 347 TargMassT2=TargetDiffStateMinMass2; 348 TargMassT =TargetDiffStateMinMass; 349 350 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2- 351 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) 352 /4./S; 353 //G4cout<<" Pt2 Mpt Mtt Pz2 "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl; 354 355 if(PZcms2 < 0 ) 356 { 357 /* 358 G4cout<<"whilecount "<<whilecount<<" "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl; 359 G4int Uzhi; G4cin>>Uzhi; 360 */ 361 return false; 362 }; 363 maxPtSquare=PZcms2; 364 365 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0); 366 Pt2=G4ThreeVector(Qmomentum.vect()).mag2(); 367 368 ProjMassT2=M0projectile2+Pt2; 369 ProjMassT =std::sqrt(ProjMassT2); 370 371 TargMassT2=TargetDiffStateMinMass2+Pt2; 372 TargMassT =std::sqrt(TargMassT2); 373 374 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2- 375 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) 376 /4./S; 377 /* 378 if(PZcms2 < 0 ) 379 { 380 G4cout<<"whilecount "<<whilecount<<" "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl; 381 G4int Uzhi; G4cin>>Uzhi; 382 return false; 383 }; 384 */ 385 if(PZcms2 < 0 ) continue; 386 PZcms =std::sqrt(PZcms2); 387 388 TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms; 389 TPlusMax=SqrtS-ProjMassT; 390 391 //G4cout<<" Tmin max "<<TPlusMin<<" "<<TPlusMax<<G4endl; 392 393 TPlusNew=ChooseP(TPlusMin, TPlusMax); 394 395 //TPlusNew=TPlusMax; 396 //G4cout<<"T+new "<<TPlusNew<<G4endl; 397 398 PPlusNew=SqrtS-TPlusNew; 399 Qplus=PPlusNew-Pprojectile.plus(); 400 PMinusNew=ProjMassT2/PPlusNew; 401 Qminus=PMinusNew-Pprojectile.minus(); 402 403 Qmomentum.setPz( (Qplus-Qminus)/2 ); 404 Qmomentum.setE( (Qplus+Qminus)/2 ); 405 406 } while ( 407 ((Pprojectile+Qmomentum).mag2() < M0projectile2 ) || //No without excitation 408 ((Ptarget -Qmomentum).mag2() < TargetDiffStateMinMass2)); 409 } 410 } 411 else //----------- Non-diffraction process ------------ 412 { 413 //G4cout<<" Non-difraction"<<G4endl; 414 do { 415 // Generate pt 416 // if (whilecount++ >= 500 && (whilecount%100)==0) 417 // G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping" 418 // << ", loop count/ maxPtSquare : " 419 // << whilecount << " / " << maxPtSquare << G4endl; 420 if (whilecount > 1000 ) 338 421 { 339 G4double PPlusNew=SqrtS-TPlusNew; 340 Qplus=PPlusNew-Pprojectile.plus(); 341 PMinusNew=ProjMassT2/PPlusNew; 342 Qminus=PMinusNew-Pprojectile.minus(); 343 344 Qmomentum.setPz( (Qplus-Qminus)/2 ); 345 Qmomentum.setE( (Qplus+Qminus)/2 ); 422 Qmomentum=G4LorentzVector(0.,0.,0.,0.); 423 return false; // Ignore this interaction 346 424 }; 425 // --------------- Check that the interaction is possible ----------- 426 ProjMassT2=ProjectileNonDiffStateMinMass2; 427 ProjMassT =ProjectileNonDiffStateMinMass; 428 429 TargMassT2=TargetNonDiffStateMinMass2; 430 TargMassT =TargetNonDiffStateMinMass; 431 432 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2- 433 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) 434 /4./S; 435 //G4cout<<" Pt2 Mpt Mtt Pz2 "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl; 436 437 if(PZcms2 < 0 ) 438 { 439 /* 440 G4cout<<"whilecount "<<whilecount<<" "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl; 441 G4int Uzhi; G4cin>>Uzhi; 442 */ 443 return false; 444 }; 445 maxPtSquare=PZcms2; 446 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0); 447 Pt2=G4ThreeVector(Qmomentum.vect()).mag2(); 448 449 ProjMassT2=ProjectileNonDiffStateMinMass2+Pt2; 450 ProjMassT =std::sqrt(ProjMassT2); 451 452 TargMassT2=TargetNonDiffStateMinMass2+Pt2; 453 TargMassT =std::sqrt(TargMassT2); 454 455 PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2- 456 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2) 457 /4./S; 458 /* 459 G4cout<<"ProjectileNonDiffStateMinMass2 "<<ProjectileNonDiffStateMinMass2<<G4endl; 460 G4cout<<"TargetNonDiffStateMinMass2 "<<TargetNonDiffStateMinMass2<<G4endl; 461 G4cout<<"Mt "<<ProjMassT<<" "<<TargMassT<<" "<<Pt2<<" "<<PZcms2<<G4endl<<G4endl; 462 */ 463 // if(PZcms2 < 0 ) {PZcms2=0;}; 464 if(PZcms2 < 0 ) continue; 465 PZcms =std::sqrt(PZcms2); 466 467 PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms; 468 PMinusMax=SqrtS-TargMassT; 469 470 PMinusNew=ChooseP(PMinusMin, PMinusMax); 471 // PMinusNew=1./sqrt(1./PMinusMin-G4UniformRand()*(1./PMinusMin-1./PMinusMax)); 472 473 //G4cout<<"Proj "<<PMinusMin<<" "<<PMinusMax<<" "<<PMinusNew<<G4endl; 474 475 //PMinusNew=PMinusMax; //+++++++++++++++++++++++++++++++++++ Vova 476 477 Qminus=PMinusNew-Pprojectile.minus(); 478 479 TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms; 480 // TPlusMax=SqrtS-PMinusNew; // Vova 481 TPlusMax=SqrtS-ProjMassT; // Vova 482 483 TPlusNew=ChooseP(TPlusMin, TPlusMax); 484 485 //G4cout<<"Targ "<<TPlusMin<<" "<<TPlusMax<<" "<<TPlusNew<<G4endl; 486 //G4cout<<PMinusNew<<" "<<TPlusNew<<G4endl; 487 488 Qplus=-(TPlusNew-Ptarget.plus()); 489 490 Qmomentum.setPz( (Qplus-Qminus)/2 ); 491 Qmomentum.setE( (Qplus+Qminus)/2 ); 492 /* 493 G4cout << "Qplus / Qminus " << Qplus << " / " << Qminus<<G4endl; 494 G4cout << "pt2" << pt2 << G4endl; 495 G4cout << "Qmomentum " << Qmomentum << G4endl; 496 G4cout << " Masses (P/T) : " << (Pprojectile+Qmomentum).mag() << 497 " / " << (Ptarget-Qmomentum).mag() << G4endl; // mag() 498 G4cout<<"Mprojectile "<<std::sqrt(M0projectile2)<<G4endl; 499 G4cout<<"Mtarget "<<std::sqrt(M0target2 )<<G4endl; 500 G4cout<<"ProjectileDiffStateMinMass "<<std::sqrt(ProjectileDiffStateMinMass2)<<G4endl; 501 G4cout<<"TargetDiffStateMinMass "<<std::sqrt(TargetDiffStateMinMass2)<<G4endl; 502 */ 503 } while ( 504 ((Pprojectile+Qmomentum).mag2() < ProjectileNonDiffStateMinMass2) || //No double Diffraction 505 ((Ptarget -Qmomentum).mag2() < TargetNonDiffStateMinMass2 )); 506 } 507 508 //G4int Uzhiinp; G4cin>>Uzhiinp; // Vova 347 509 348 510 Pprojectile += Qmomentum; 349 511 Ptarget -= Qmomentum; 350 351 /* 352 Pprojectile.setPz(0.); 353 Pprojectile.setE(SqrtS-M0target); 354 355 Ptarget.setPz(0.); 356 Ptarget.setE(M0target); 357 */ 358 359 //G4cout << "Pprojectile with Q : " << Pprojectile << G4endl; 360 //G4cout << "Ptarget with Q : " << Ptarget << G4endl; 361 362 // G4cout << "Projectile back: " << toLab * Pprojectile << G4endl; 363 // G4cout << "Target back: " << toLab * Ptarget << G4endl; 512 /* 513 G4cout << "Pprojectile with Q : " << Pprojectile << G4endl; 514 G4cout << "Ptarget with Q : " << Ptarget << G4endl; 515 G4cout << "Target mass " << Ptarget.mag() << G4endl; 516 G4cout << "Projectile mass " << Pprojectile.mag() << G4endl; 517 // 518 //G4cout << "Projectile back: " << toLab * Pprojectile << G4endl; 519 //G4cout << "Target back: " << toLab * Ptarget << G4endl; 520 */ 521 //-------------- Flip if projectale moves in backward direction ------------ 522 //G4bool Flip=Pprojectile.pz()< 0.; 523 364 524 365 525 // Transform back and update SplitableHadron Participant. … … 369 529 //G4cout << "Pprojectile with Q M: " << Pprojectile<<" "<< Pprojectile.mag() << G4endl; 370 530 //G4cout << "Ptarget with Q M: " << Ptarget <<" "<< Ptarget.mag() << G4endl; 371 372 //G4cout << "Target mass " << Ptarget.mag() << G4endl; 373 531 //G4cout << "Target mass " << Ptarget.mag() << G4endl; 532 //G4cout << "Projectile mass " << Pprojectile.mag() << G4endl; 533 534 /* 535 if(!Flip){ 536 projectile->Set4Momentum(Pprojectile); 374 537 target->Set4Momentum(Ptarget); 375 376 //G4cout << "Projectile mass " << Pprojectile.mag() << G4endl; 538 } 539 else { 540 G4ParticleDefinition * t_Definition=projectile->GetDefinition(); 541 projectile->SetDefinition(target->GetDefinition()); 542 projectile->Set4Momentum(Ptarget); 543 target->SetDefinition(t_Definition); 544 target->Set4Momentum(Pprojectile); 545 } 546 */ 547 // 548 /* 549 if(G4UniformRand() < 1.) { 550 G4ParticleDefinition * t_Definition=projectile->GetDefinition(); 551 projectile->SetDefinition(target->GetDefinition()); 552 target->SetDefinition(t_Definition); 553 } 554 */ // For flip, for HARP 555 556 G4double ZcoordinateOfCurrentInteraction = target->GetPosition().z(); 557 // It is assumed that nucleon z-coordinates are ordered on increasing ----------- 558 559 G4double betta_z=projectile->Get4Momentum().pz()/projectile->Get4Momentum().e(); 560 561 G4double ZcoordinateOfPreviousCollision=projectile->GetPosition().z(); 562 if(projectile->GetSoftCollisionCount()==0) { 563 projectile->SetTimeOfCreation(0.); 564 target->SetTimeOfCreation(0.); 565 ZcoordinateOfPreviousCollision=ZcoordinateOfCurrentInteraction; 566 } 567 568 G4ThreeVector thePosition(projectile->GetPosition().x(), 569 projectile->GetPosition().y(), 570 ZcoordinateOfCurrentInteraction); 571 projectile->SetPosition(thePosition); 572 573 G4double TimeOfPreviousCollision=projectile->GetTimeOfCreation(); 574 G4double TimeOfCurrentCollision=TimeOfPreviousCollision+ 575 (ZcoordinateOfCurrentInteraction-ZcoordinateOfPreviousCollision)/betta_z; 576 577 projectile->SetTimeOfCreation(TimeOfCurrentCollision); 578 target->SetTimeOfCreation(TimeOfCurrentCollision); 377 579 378 580 projectile->Set4Momentum(Pprojectile); 581 target->Set4Momentum(Ptarget); 582 583 projectile->IncrementCollisionCount(1); 584 target->IncrementCollisionCount(1); 585 586 // 587 //G4cout<<"Out of Excitation --------------------"<<G4endl; 588 //G4int Uzhiinp; G4cin>>Uzhiinp; // Vova 379 589 380 590 return true; 381 591 } 382 592 383 593 // --------------------------------------------------------------------- 384 594 G4ExcitedString * G4DiffractiveExcitation:: 385 595 String(G4VSplitableHadron * hadron, G4bool isProjectile) const 386 596 { 597 598 //G4cout<<"G4DiffractiveExcitation::String isProj"<<isProjectile<<G4endl; 599 387 600 hadron->SplitUp(); 388 601 G4Parton *start= hadron->GetNextParton(); 389 if ( start==NULL) 602 if ( start==NULL) 390 603 { G4cout << " G4FTFModel::String() Error:No start parton found"<< G4endl; 391 604 return NULL; 392 605 } 393 606 G4Parton *end = hadron->GetNextParton(); 394 if ( end==NULL) 607 if ( end==NULL) 395 608 { G4cout << " G4FTFModel::String() Error:No end parton found"<< G4endl; 396 609 return NULL; 397 610 } 398 611 399 612 G4ExcitedString * string; 400 if ( isProjectile ) 613 if ( isProjectile ) 401 614 { 402 615 string= new G4ExcitedString(end,start, +1); … … 404 617 string= new G4ExcitedString(start,end, -1); 405 618 } 406 619 // Uzhi 620 //G4cout<<"G4ExcitedString * G4DiffractiveExcitation::String"<<G4endl; 621 //G4cout<<hadron->GetTimeOfCreation()<<" "<<hadron->GetPosition()/fermi<<G4endl; 622 623 string->SetTimeOfCreation(hadron->GetTimeOfCreation()); 407 624 string->SetPosition(hadron->GetPosition()); 408 625 409 626 // momenta of string ends 627 // 628 G4double Momentum=hadron->Get4Momentum().vect().mag(); 629 G4double Plus=hadron->Get4Momentum().e() + Momentum; 630 G4double Minus=hadron->Get4Momentum().e() - Momentum; 631 632 G4ThreeVector tmp; 633 if(Momentum > 0.) 634 { 635 tmp.set(hadron->Get4Momentum().px(), 636 hadron->Get4Momentum().py(), 637 hadron->Get4Momentum().pz()); 638 tmp/=Momentum; 639 } 640 else 641 { 642 tmp.set(0.,0.,1.); 643 }; 644 645 G4LorentzVector Pstart(tmp,0.); 646 G4LorentzVector Pend(tmp,0.); 647 648 if(isProjectile) 649 { 650 Pstart*=(-1.)*Minus/2.; 651 Pend *=(+1.)*Plus /2.; 652 } 653 else 654 { 655 Pstart*=(+1.)*Plus/2.; 656 Pend *=(-1.)*Minus/2.; 657 }; 658 659 Momentum=-Pstart.mag(); 660 Pstart.setT(Momentum); // It is assumed that quark has m=0. 661 662 Momentum=-Pend.mag(); 663 Pend.setT(Momentum); // It is assumed that di-quark has m=0. 664 // 665 /* Uzhi 410 666 G4double ptSquared= hadron->Get4Momentum().perp2(); 411 667 G4double transverseMassSquared= hadron->Get4Momentum().plus() … … 416 672 sqr( std::sqrt(transverseMassSquared) - std::sqrt(ptSquared) ); 417 673 418 G4double widthOfPtSquare = 0.25 ; // Uzhi<Pt^2>=0.25 ??????????????????674 G4double widthOfPtSquare = 0.25*GeV*GeV; // Uzhi 11.07 <Pt^2>=0.25 ?????????????????? 419 675 G4ThreeVector pt=GaussianPt(widthOfPtSquare,maxAvailMomentumSquared); 420 676 … … 432 688 433 689 G4int Sign= isProjectile ? -1 : 1; 434 690 435 691 G4double endMinus = 0.5 * (tm1 + Sign*tm2); 436 692 G4double startMinus= hadron->Get4Momentum().minus() - endMinus; … … 444 700 Pend.setPz(0.5*(endPlus - endMinus)); 445 701 Pend.setE(0.5*(endPlus + endMinus)); 446 702 */ // Uzhi 447 703 start->Set4Momentum(Pstart); 448 704 end->Set4Momentum(Pend); 449 705 /* 706 G4cout<<"G4DiffractiveExcitation::String hadro"<<hadron->Get4Momentum()<<" "<<hadron->Get4Momentum().mag2()<<G4endl; 707 708 G4cout<<"G4DiffractiveExcitation::String start"<<start->Get4Momentum()<<" "<<start->GetPDGcode()<<G4endl; 709 710 G4cout<<"G4DiffractiveExcitation::String end "<< end->Get4Momentum()<<" "<< end->GetPDGcode()<<G4endl; 711 G4int Uzhi; G4cin>>Uzhi; 712 */ 450 713 #ifdef G4_FTFDEBUG 451 G4cout << " generated string flavors " << start->GetPDGcode() << " / " << end->GetPDGcode() << G4endl; 452 G4cout << " generated string momenta: quark " << start->Get4Momentum() << "mass : " <<start->Get4Momentum().mag()<< G4endl; 453 G4cout << " generated string momenta: Diquark " << end ->Get4Momentum() << "mass : " <<end->Get4Momentum().mag()<< G4endl; 714 G4cout << " generated string flavors " 715 << start->GetPDGcode() << " / " 716 << end->GetPDGcode() << G4endl; 717 G4cout << " generated string momenta: quark " 718 << start->Get4Momentum() << "mass : " 719 <<start->Get4Momentum().mag() << G4endl; 720 G4cout << " generated string momenta: Diquark " 721 << end ->Get4Momentum() 722 << "mass : " <<end->Get4Momentum().mag()<< G4endl; 454 723 G4cout << " sum of ends " << Pstart+Pend << G4endl; 455 724 G4cout << " Original " << hadron->Get4Momentum() << G4endl; 456 725 #endif 457 458 return string; 726 727 return string; 459 728 } 460 729 … … 462 731 // --------- private methods ---------------------- 463 732 733 // --------------------------------------------------------------------- 464 734 G4double G4DiffractiveExcitation::ChooseP(G4double Pmin, G4double Pmax) const // Uzhi 465 735 { 466 736 // choose an x between Xmin and Xmax with P(x) ~ 1/x 467 468 737 // to be improved... 469 738 470 739 G4double range=Pmax-Pmin; // Uzhi 471 472 if ( Pmin <= 0. || range <=0. ) 740 741 if ( Pmin <= 0. || range <=0. ) 473 742 { 474 743 G4cout << " Pmin, range : " << Pmin << " , " << range << G4endl; … … 489 758 } 490 759 491 G4ThreeVector G4DiffractiveExcitation::GaussianPt(G4double AveragePt2, G4double maxPtSquare) const // Uzhi 760 // --------------------------------------------------------------------- 761 G4ThreeVector G4DiffractiveExcitation::GaussianPt(G4double AveragePt2, 762 G4double maxPtSquare) const // Uzhi 492 763 { // @@ this method is used in FTFModel as well. Should go somewhere common! 493 764 494 765 G4double Pt2; 495 766 /* // Uzhi … … 499 770 */ // Uzhi 500 771 501 Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() * (std::exp(-maxPtSquare/AveragePt2)-1.));// Uzhi 502 772 Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() * 773 (std::exp(-maxPtSquare/AveragePt2)-1.));// Uzhi 774 503 775 G4double Pt=std::sqrt(Pt2); 504 776 505 777 G4double phi=G4UniformRand() * twopi; 506 507 return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.); 508 } 509 778 779 return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.); 780 } 781 782 // --------------------------------------------------------------------- 510 783 G4DiffractiveExcitation::G4DiffractiveExcitation(const G4DiffractiveExcitation &) 511 784 { -
trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4DiffractiveHHScatterer.cc
r819 r962 30 30 #include "G4KineticTrack.hh" 31 31 #include "G4DiffractiveSplitableHadron.hh" 32 #include "G4FTFParameters.hh" // Uzhi 21.04.08 32 33 33 34 G4DiffractiveHHScatterer::G4DiffractiveHHScatterer() … … 37 38 {} 38 39 40 // ------------------------------------------------------------------- 39 41 G4KineticTrackVector * G4DiffractiveHHScatterer:: 40 42 Scatter(const G4KineticTrack & aTrack, const G4KineticTrack & bTrack) … … 44 46 G4DiffractiveSplitableHadron aHadron(& aTrack); 45 47 G4DiffractiveSplitableHadron bHadron(& bTrack); 46 if ( ! theExcitation->ExciteParticipants(& aHadron, & bHadron)) 48 theParameters = new G4FTFParameters(aHadron.GetDefinition(), // -------- Uzhi 21.04.08 49 1.,1., 100.); 50 //s);// ------------------------- Uzhi 21.04.08 51 if ( ! theExcitation->ExciteParticipants(& aHadron, 52 & bHadron, 53 theParameters)) // -------- Uzhi 21.04.08 47 54 { 48 55 return NULL; -
trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4DiffractiveSplitableHadron.cc
r819 r962 25 25 // 26 26 // 27 // $Id: G4DiffractiveSplitableHadron.cc,v 1. 6 2006/06/29 20:54:36 gunterExp $28 // GEANT4 tag $Name: geant4-09-0 1-patch-02 $27 // $Id: G4DiffractiveSplitableHadron.cc,v 1.7 2008/03/31 15:34:01 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 -
trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4FTFModel.cc
r819 r962 25 25 // 26 26 // 27 // $Id: G4FTFModel.cc,v 1. 7 2007/04/24 10:32:59 gunterExp $28 // GEANT4 tag $Name: geant4-09-0 1-patch-02 $27 // $Id: G4FTFModel.cc,v 1.13 2008/12/09 10:40:52 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 … … 38 38 39 39 #include "G4FTFModel.hh" 40 #include "G4FTFParameters.hh" // Uzhi 29.03.08 40 41 #include "G4FTFParticipants.hh" 41 42 #include "G4InteractionContent.hh" … … 43 44 #include "G4ParticleDefinition.hh" 44 45 #include "G4ios.hh" 46 #include <utility> // Uzhi 29.03.08 45 47 46 48 // Class G4FTFModel 47 49 48 G4FTFModel::G4FTFModel():theExcitation(new G4DiffractiveExcitation()) // Uzhi 50 G4FTFModel::G4FTFModel():theExcitation(new G4DiffractiveExcitation()), 51 theElastic(new G4ElasticHNScattering()) // Uzhi 29.03.08 49 52 { 50 53 G4VPartonStringModel::SetThisPointer(this); 51 } 52 53 G4FTFModel::G4FTFModel(G4double a, G4double b, G4double c):theExcitation(new G4DiffractiveExcitation()) 54 theParameters=0; // Uzhi 9.12.08 55 } 56 57 /* 58 G4FTFModel::G4FTFModel(G4double , G4double , G4double ):theExcitation(new // Uzhi 9.12.08 G4DiffractiveExcitation()) 54 59 { 55 60 G4VPartonStringModel::SetThisPointer(this); … … 62 67 G4VPartonStringModel::SetThisPointer(this); 63 68 } 64 69 */ 65 70 66 71 67 72 G4FTFModel::~G4FTFModel() 68 {} 73 { 74 if( theParameters != 0 ) delete theParameters; // Uzhi 5.12.08 75 // Because FTF model can be called for various particles 76 // theParameters must be erased at the end of each call. 77 // Thus the delete is olso in G4FTFModel::GetStrings() method 78 if( theExcitation != 0 ) delete theExcitation; // Uzhi 5.12.08 79 if( theElastic != 0 ) delete theElastic; // Uzhi 5.12.08 80 } 69 81 70 82 … … 86 98 } 87 99 100 // ------------------------------------------------------------ 88 101 void G4FTFModel::Init(const G4Nucleus & aNucleus, const G4DynamicParticle & aProjectile) 89 102 { 90 theParticipants.Init(aNucleus.GetN(),aNucleus.GetZ()); // Uzhi N-mass number Z-charge91 92 103 theProjectile = aProjectile; 93 } 94 104 //G4cout<<"G4FTFModel::Init "<<aNucleus.GetN()<<" "<<aNucleus.GetZ()<<G4endl; 105 theParticipants.Init(aNucleus.GetN(),aNucleus.GetZ()); 106 // Uzhi N-mass number Z-charge ------------------------- Uzhi 29.03.08 107 108 // --- cms energy 109 110 G4double s = sqr( theProjectile.GetMass() ) + 111 sqr( G4Proton::Proton()->GetPDGMass() ) + 112 2*theProjectile.GetTotalEnergy()*G4Proton::Proton()->GetPDGMass(); 113 /* 114 G4cout << " primary Total E (GeV): " << theProjectile.GetTotalEnergy()/GeV << G4endl; 115 G4cout << " primary Mass (GeV): " << theProjectile.GetMass() /GeV << G4endl; 116 G4cout << "cms std::sqrt(s) (GeV) = " << std::sqrt(s) / GeV << G4endl; 117 */ 118 119 if( theParameters != 0 ) delete theParameters; // Uzhi 9.12.08 120 theParameters = new G4FTFParameters(theProjectile.GetDefinition(), 121 aNucleus.GetN(),aNucleus.GetZ(), 122 s);// ------------------------- Uzhi 19.04.08 123 //theParameters->SetProbabilityOfElasticScatt(0.); // To turn on/off (1/0) elastic scattering 124 } 125 126 // ------------------------------------------------------------ 95 127 G4ExcitedStringVector * G4FTFModel::GetStrings() 96 128 { 97 theParticipants.BuildInteractions(theProjectile); 98 129 //G4cout<<"theParticipants.GetList"<<G4endl; 130 theParticipants.GetList(theProjectile,theParameters); 131 //G4cout<<"ExciteParticipants()"<<G4endl; 99 132 if (! ExciteParticipants()) return NULL;; 100 133 //G4cout<<"theStrings = BuildStrings()"<<G4endl; 101 134 G4ExcitedStringVector * theStrings = BuildStrings(); 102 135 //G4cout<<"Return to theStrings "<<G4endl; 136 if( theParameters != 0 ) // Uzhi 9.12.08 137 { // Uzhi 9.12.08 138 delete theParameters; // Uzhi 9.12.08 139 theParameters=0; // Uzhi 9.12.08 140 } // Uzhi 9.12.08 103 141 return theStrings; 104 142 } 105 143 144 // ------------------------------------------------------------ 106 145 struct DeleteVSplitableHadron { void operator()(G4VSplitableHadron * aH){delete aH;} }; 107 146 108 G4ExcitedStringVector * G4FTFModel::BuildStrings() 109 { 110 111 // Loop over all collisions; find all primaries, and all target ( targets may 112 // be duplicate in the List ( to unique G4VSplitableHadrons) 113 114 G4ExcitedStringVector * strings; 115 strings = new G4ExcitedStringVector(); 116 117 std::vector<G4VSplitableHadron *> primaries; 118 std::vector<G4VSplitableHadron *> targets; 119 120 theParticipants.StartLoop(); // restart a loop 121 while ( theParticipants.Next() ) 122 { 123 const G4InteractionContent & interaction=theParticipants.GetInteraction(); 124 // do not allow for duplicates ... 125 if ( primaries.end() == std::find(primaries.begin(), primaries.end(), interaction.GetProjectile()) ) 126 primaries.push_back(interaction.GetProjectile()); 127 128 if ( targets.end() == std::find(targets.begin(), targets.end(),interaction.GetTarget()) ) 129 targets.push_back(interaction.GetTarget()); 130 } 131 132 133 // G4cout << "BuildStrings prim/targ " << primaries.entries() << " , " << 134 // targets.entries() << G4endl; 135 136 137 unsigned int ahadron; 138 for ( ahadron=0; ahadron < primaries.size() ; ahadron++) 139 { 140 G4bool isProjectile=true; 141 strings->push_back(theExcitation->String(primaries[ahadron], isProjectile)); 142 } 143 for ( ahadron=0; ahadron < targets.size() ; ahadron++) 144 { 145 G4bool isProjectile=false; 146 strings->push_back(theExcitation->String(targets[ahadron], isProjectile)); 147 } 148 149 std::for_each(primaries.begin(), primaries.end(), DeleteVSplitableHadron()); 150 primaries.clear(); 151 std::for_each(targets.begin(), targets.end(), DeleteVSplitableHadron()); 152 targets.clear(); 153 154 return strings; 155 } 156 147 // ------------------------------------------------------------ 157 148 G4bool G4FTFModel::ExciteParticipants() 158 149 { 159 150 /* // Uzhi 29.03.08 For elastic Scatt. 151 G4cout<<" In ExciteParticipants() "<<theParticipants.theInteractions.size()<<G4endl; 152 G4cout<<" test Params Tot "<<theParameters->GetTotalCrossSection()<<G4endl; 153 G4cout<<" test Params Ela "<<theParameters->GetElasticCrossSection()<<G4endl; 154 155 G4int counter=0; 156 */ // Uzhi 29.03.08 157 158 159 160 160 while (theParticipants.Next()) 161 161 { 162 162 const G4InteractionContent & collision=theParticipants.GetInteraction(); 163 164 //G4cout << " soft colls : " << collision.GetNumberOfSoftCollisions() << G4endl; // Uzhi no match 163 /* 164 counter++; 165 G4cout<<" Inter # "<<counter<<G4endl; 166 */ 165 167 G4VSplitableHadron * projectile=collision.GetProjectile(); 166 168 G4VSplitableHadron * target=collision.GetTarget(); 167 169 168 if ( ! theExcitation->ExciteParticipants(projectile, target) ) 170 // // Uzhi 29.03.08 171 G4bool Successfull; 172 if(G4UniformRand()< theParameters->GetProbabilityOfElasticScatt()) 173 { 174 //G4cout<<"Elastic"<<G4endl; 175 Successfull=theElastic->ElasticScattering(projectile, target, theParameters); 176 } 177 else 178 { 179 //G4cout<<"Inelastic"<<G4endl; 180 Successfull=theExcitation->ExciteParticipants(projectile, target, theParameters); 181 } 182 // if(!Successfull) 183 // // Uzhi 29.03.08 184 185 // if ( ! theExcitation->ExciteParticipants(projectile, target) ) 186 if(!Successfull) 169 187 { 170 188 // give up, clean up … … 176 194 const G4InteractionContent & interaction=theParticipants.GetInteraction(); 177 195 // do not allow for duplicates ... 178 if ( primaries.end() == std::find(primaries.begin(), primaries.end(), interaction.GetProjectile()) ) 196 if ( primaries.end() == std::find(primaries.begin(), primaries.end(), 197 interaction.GetProjectile()) ) 179 198 primaries.push_back(interaction.GetProjectile()); 180 199 181 if ( targets.end() == std::find(targets.begin(), targets.end(),interaction.GetTarget()) ) 200 if ( targets.end() == std::find(targets.begin(), targets.end(), 201 interaction.GetTarget()) ) 182 202 targets.push_back(interaction.GetTarget()); 183 203 } … … 188 208 189 209 return false; 190 } 191 210 } // End of the loop Uzhi 211 } 212 return true; 213 } 214 // ------------------------------------------------------------ 215 G4ExcitedStringVector * G4FTFModel::BuildStrings() 216 { 217 // Loop over all collisions; find all primaries, and all target ( targets may 218 // be duplicate in the List ( to unique G4VSplitableHadrons) 219 220 G4ExcitedStringVector * strings; 221 strings = new G4ExcitedStringVector(); 222 223 std::vector<G4VSplitableHadron *> primaries; 224 std::vector<G4VSplitableHadron *> targets; 225 226 theParticipants.StartLoop(); // restart a loop 227 while ( theParticipants.Next() ) 228 { 229 const G4InteractionContent & interaction=theParticipants.GetInteraction(); 230 // do not allow for duplicates ... 231 if ( primaries.end() == std::find(primaries.begin(), primaries.end(), 232 interaction.GetProjectile()) ) 233 primaries.push_back(interaction.GetProjectile()); 234 235 if ( targets.end() == std::find(targets.begin(), targets.end(), 236 interaction.GetTarget()) ) 237 targets.push_back(interaction.GetTarget()); 192 238 } 193 return true; 194 } 239 240 241 // G4cout << "BuildStrings prim/targ " << primaries.size() << " , " << 242 // targets.size() << G4endl; 243 244 unsigned int ahadron; 245 // Only for hA-interactions Uzhi ------------------------------------- 246 for ( ahadron=0; ahadron < primaries.size() ; ahadron++) 247 { 248 //G4ThreeVector aPosition=primaries[ahadron]->GetPosition(); 249 //G4cout<<"Proj Build "<<aPosition<<" "<<primaries[ahadron]->GetTimeOfCreation()<<G4endl; 250 G4bool isProjectile=true; 251 strings->push_back(theExcitation->String(primaries[ahadron], isProjectile)); 252 } 253 254 for ( ahadron=0; ahadron < targets.size() ; ahadron++) 255 { 256 //G4ThreeVector aPosition=targets[ahadron]->GetPosition(); 257 //G4cout<<"Targ Build "<<aPosition<<" "<<targets[ahadron]->GetTimeOfCreation()<<G4endl; 258 G4bool isProjectile=false; 259 strings->push_back(theExcitation->String(targets[ahadron], isProjectile)); 260 } 261 262 std::for_each(primaries.begin(), primaries.end(), DeleteVSplitableHadron()); 263 primaries.clear(); 264 std::for_each(targets.begin(), targets.end(), DeleteVSplitableHadron()); 265 targets.clear(); 266 267 return strings; 268 } 269 // ------------------------------------------------------------ -
trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4FTFParticipants.cc
r819 r962 25 25 // 26 26 // 27 // $Id: G4FTFParticipants.cc,v 1. 7 2007/04/24 10:33:00 gunterExp $28 // GEANT4 tag $Name: geant4-09-0 1-patch-02 $27 // $Id: G4FTFParticipants.cc,v 1.9 2008/06/13 12:49:23 vuzhinsk Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // ------------------------------------------------------------ … … 38 38 // ------------------------------------------------------------ 39 39 40 #include "G4FTFParameters.hh" // Uzhi 29.03.08 40 41 #include "G4FTFParticipants.hh" 41 42 #include "G4DiffractiveSplitableHadron.hh" 42 43 #include "G4VSplitableHadron.hh" 43 //#include "G4PomeronCrossSection.hh" // Uzhi44 #include "G4FTFCrossSection.hh" // Uzhi45 44 #include "Randomize.hh" 46 #include <utility> 47 45 #include <utility> // Uzhi 29.03.08 48 46 49 47 // Class G4FTFParticipants 50 51 52 48 53 49 G4FTFParticipants::G4FTFParticipants() … … 75 71 //{} 76 72 77 void G4FTFParticipants::BuildInteractions(const G4ReactionProduct &thePrimary) 73 void G4FTFParticipants::GetList(const G4ReactionProduct &thePrimary, 74 G4FTFParameters *theParameters) // Uzhi 29.03.08 78 75 { 79 76 … … 83 80 theInteractions.clear(); 84 81 85 // --- cms energy 86 87 G4double s = sqr( thePrimary.GetMass() ) + 88 sqr( G4Proton::Proton()->GetPDGMass() ) + 89 2*thePrimary.GetTotalEnergy()*G4Proton::Proton()->GetPDGMass(); 90 91 // G4cout << " primary Total E (GeV): " << thePrimary.GetTotalEnergy()/GeV << G4endl; 92 // G4cout << " primary Mass (GeV): " << thePrimary.GetMass() /GeV << G4endl; 93 // G4cout << "cms std::sqrt(s) (GeV) = " << std::sqrt(s) / GeV << G4endl; 94 95 // G4PomeronCrossSection theCrossSection(thePrimary.GetDefinition()); // Uzhi 96 G4FTFCrossSection theCrossSection(thePrimary.GetDefinition(),s); // Uzhi 97 98 99 G4double deltaxy=2 * fermi; 82 G4double deltaxy=2 * fermi; // Extra nuclear radius 100 83 101 84 G4VSplitableHadron * primarySplitable=new G4DiffractiveSplitableHadron(thePrimary); 102 85 103 G4double xyradius; 104 xyradius =theNucleus->GetOuterRadius() + deltaxy; 105 86 G4double xyradius; 87 xyradius =theNucleus->GetOuterRadius() + deltaxy; // Impact parameter sampling 88 // radius 106 89 G4bool nucleusNeedsShift = true; 107 90 … … 115 98 theNucleus->StartLoop(); 116 99 G4Nucleon * nucleon; 117 while ( (nucleon=theNucleus->GetNextNucleon()) ) 100 //G4int InterNumber=0; // Uzhi 101 //while ( (nucleon=theNucleus->GetNextNucleon())&& (InterNumber < 1) ) // Uzhi 102 while ( (nucleon=theNucleus->GetNextNucleon()) ) // Uzhi 118 103 { 119 104 G4double impact2= sqr(impactX - nucleon->GetPosition().x()) + 120 sqr(impactY - nucleon->GetPosition().y()); 121 // if ( theCrossSection.GetInelasticProbability(s,impact2) // Uzhi 122 if ( theCrossSection.GetInelasticProbability( impact2/fermi/fermi) // Uzhi 105 sqr(impactY - nucleon->GetPosition().y()); 106 107 // if ( theParameters->GetInelasticProbability(impact2/fermi/fermi) // Uzhi 29.03.08 108 if ( theParameters->GetProbabilityOfInteraction(impact2/fermi/fermi) // Uzhi 29.03.08 123 109 > G4UniformRand() ) 124 110 { 111 //InterNumber++; 125 112 if ( nucleusNeedsShift ) 126 113 { // on the first hit, shift nucleus … … 136 123 nucleon->Hit(targetSplitable); 137 124 } 138 G4InteractionContent * aInteraction = new G4InteractionContent(primarySplitable); 125 G4InteractionContent * aInteraction = 126 new G4InteractionContent(primarySplitable); 139 127 aInteraction->SetTarget(targetSplitable); 140 128 theInteractions.push_back(aInteraction); … … 152 140 153 141 // Implementation (private) methods 154 155 156 157 158
Note: See TracChangeset
for help on using the changeset viewer.