Changeset 962 for trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4DiffractiveExcitation.cc
- Timestamp:
- Apr 6, 2009, 12:30:29 PM (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
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 {
Note: See TracChangeset
for help on using the changeset viewer.
