Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // 26 // >> 27 // $Id: G4QGSDiffractiveExcitation.cc 100828 2016-11-02 15:25:59Z gcosmo $ 27 // ------------------------------------------- 28 // ------------------------------------------------------------ 28 // GEANT 4 class implemetation file 29 // GEANT 4 class implemetation file 29 // 30 // 30 // ---------------- G4QGSDiffractiveExcit 31 // ---------------- G4QGSDiffractiveExcitation -------------- 31 // by Gunter Folger, October 1998. 32 // by Gunter Folger, October 1998. 32 // diffractive Excitation used by stri 33 // diffractive Excitation used by strings models 33 // Take a projectile and a target 34 // Take a projectile and a target 34 // excite the projectile and target 35 // excite the projectile and target 35 // Essential changed by V. Uzhinsky in Novemb 36 // Essential changed by V. Uzhinsky in November - December 2006 36 // in order to put it in a correspondence wit 37 // in order to put it in a correspondence with original FRITIOF 37 // model. Variant of FRITIOF with nucleon de- 38 // model. Variant of FRITIOF with nucleon de-excitation is implemented. 38 // ------------------------------------------- 39 // --------------------------------------------------------------------- 39 40 40 // Modified: 41 // Modified: 41 // 25-05-07 : G.Folger 42 // 25-05-07 : G.Folger 42 // move from management/G4DiffractiveExc 43 // move from management/G4DiffractiveExcitation to to qgsm/G4QGSDiffractiveExcitation 43 // 44 // 44 45 45 #include "globals.hh" 46 #include "globals.hh" 46 #include "G4PhysicalConstants.hh" 47 #include "G4PhysicalConstants.hh" 47 #include "G4SystemOfUnits.hh" 48 #include "G4SystemOfUnits.hh" 48 #include "Randomize.hh" 49 #include "Randomize.hh" 49 50 50 #include "G4QGSDiffractiveExcitation.hh" 51 #include "G4QGSDiffractiveExcitation.hh" 51 #include "G4LorentzRotation.hh" 52 #include "G4LorentzRotation.hh" 52 #include "G4ThreeVector.hh" 53 #include "G4ThreeVector.hh" 53 #include "G4ParticleDefinition.hh" 54 #include "G4ParticleDefinition.hh" 54 #include "G4VSplitableHadron.hh" 55 #include "G4VSplitableHadron.hh" 55 #include "G4ExcitedString.hh" 56 #include "G4ExcitedString.hh" 56 //#include "G4ios.hh" 57 //#include "G4ios.hh" 57 58 58 #include "G4Exp.hh" 59 #include "G4Exp.hh" 59 #include "G4Log.hh" 60 #include "G4Log.hh" 60 #include "G4Pow.hh" 61 #include "G4Pow.hh" 61 62 62 //============================================ << 63 << 64 //#define debugDoubleDiffraction << 65 << 66 //============================================ << 67 63 68 G4QGSDiffractiveExcitation::G4QGSDiffractiveEx 64 G4QGSDiffractiveExcitation::G4QGSDiffractiveExcitation() 69 { 65 { 70 } 66 } 71 67 72 G4QGSDiffractiveExcitation::~G4QGSDiffractiveE 68 G4QGSDiffractiveExcitation::~G4QGSDiffractiveExcitation() 73 { 69 { 74 } 70 } 75 71 76 72 77 G4bool G4QGSDiffractiveExcitation:: 73 G4bool G4QGSDiffractiveExcitation:: 78 ExciteParticipants(G4VSplitableHadron *project << 74 ExciteParticipants(G4VSplitableHadron *projectile, G4VSplitableHadron *target) const 79 { 75 { 80 #ifdef debugDoubleDiffraction << 81 G4cout<<G4endl<<"G4QGSDiffractiveExcitatio << 82 G4cout<<"Proj Targ "<<projectile->GetDefin << 83 G4cout<<"Proj 4 Mom "<<projectile->Get4Mom << 84 G4cout<<"Targ 4 Mom "<<target->Get4Momentu << 85 #endif << 86 << 87 G4LorentzVector Pprojectile=projectile->Get4 76 G4LorentzVector Pprojectile=projectile->Get4Momentum(); 88 77 89 // -------------------- Projectile parameter 78 // -------------------- Projectile parameters ----------------------------------- 90 G4bool PutOnMassShell=0; 79 G4bool PutOnMassShell=0; 91 80 >> 81 //G4double M0projectile=projectile->GetDefinition()->GetPDGMass(); // With de-excitation 92 G4double M0projectile = Pprojectile.mag(); 82 G4double M0projectile = Pprojectile.mag(); // Without de-excitation 93 83 94 if(M0projectile < projectile->GetDefinition( 84 if(M0projectile < projectile->GetDefinition()->GetPDGMass()) 95 { 85 { 96 PutOnMassShell=1; 86 PutOnMassShell=1; 97 M0projectile=projectile->GetDefinition()-> 87 M0projectile=projectile->GetDefinition()->GetPDGMass(); 98 } 88 } 99 89 >> 90 G4double Mprojectile2 = M0projectile * M0projectile; >> 91 >> 92 G4int PDGcode=projectile->GetDefinition()->GetPDGEncoding(); >> 93 G4int absPDGcode=std::abs(PDGcode); >> 94 G4double ProjectileDiffCut; >> 95 G4double AveragePt2; >> 96 >> 97 if( absPDGcode > 1000 ) //------Projectile is baryon -------- >> 98 { >> 99 ProjectileDiffCut = 1.1; // GeV >> 100 AveragePt2 = 0.3; // GeV^2 >> 101 } >> 102 else if( absPDGcode == 211 || PDGcode == 111) //------Projectile is Pion ----------- >> 103 { >> 104 ProjectileDiffCut = 1.0; // GeV >> 105 AveragePt2 = 0.3; // GeV^2 >> 106 } >> 107 else if( absPDGcode == 321 || PDGcode == -311) //------Projectile is Kaon ----------- >> 108 { >> 109 ProjectileDiffCut = 1.1; // GeV >> 110 AveragePt2 = 0.3; // GeV^2 >> 111 } >> 112 else //------Projectile is undefined, Nucleon assumed >> 113 { >> 114 ProjectileDiffCut = 1.1; // GeV >> 115 AveragePt2 = 0.3; // GeV^2 >> 116 }; >> 117 >> 118 ProjectileDiffCut = ProjectileDiffCut * GeV; >> 119 AveragePt2 = AveragePt2 * GeV*GeV; >> 120 100 // -------------------- Target parameters -- 121 // -------------------- Target parameters ---------------------------------------------- 101 G4LorentzVector Ptarget=target->Get4Momentum 122 G4LorentzVector Ptarget=target->Get4Momentum(); 102 123 103 G4double M0target = Ptarget.mag(); 124 G4double M0target = Ptarget.mag(); 104 125 105 if(M0target < target->GetDefinition()->GetPD 126 if(M0target < target->GetDefinition()->GetPDGMass()) 106 { 127 { 107 PutOnMassShell=1; 128 PutOnMassShell=1; 108 M0target=target->GetDefinition()->GetPDGMa 129 M0target=target->GetDefinition()->GetPDGMass(); 109 } 130 } 110 131 111 G4LorentzVector Psum=Pprojectile+Ptarget; << 132 G4double Mtarget2 = M0target * M0target; //Ptarget.mag2(); // for AA-inter. 112 G4double S=Psum.mag2(); << 113 G4double SqrtS=std::sqrt(S); << 114 133 115 if (SqrtS < M0projectile + M0target) {return << 134 G4double NuclearNucleonDiffCut = 1.1*GeV; 116 << 117 135 118 G4double Mprojectile2 = M0projectile * M0pro << 136 G4double ProjectileDiffCut2 = ProjectileDiffCut * ProjectileDiffCut; 119 G4double Mtarget2 = M0target * M0tar << 137 G4double NuclearNucleonDiffCut2 = NuclearNucleonDiffCut * NuclearNucleonDiffCut; 120 138 121 // Transform momenta to cms and then rotate 139 // Transform momenta to cms and then rotate parallel to z axis; 122 140 >> 141 G4LorentzVector Psum; >> 142 Psum=Pprojectile+Ptarget; >> 143 123 G4LorentzRotation toCms(-1*Psum.boostVector( 144 G4LorentzRotation toCms(-1*Psum.boostVector()); 124 145 125 G4LorentzVector Ptmp=toCms*Pprojectile; 146 G4LorentzVector Ptmp=toCms*Pprojectile; 126 147 127 if ( Ptmp.pz() <= 0. ) 148 if ( Ptmp.pz() <= 0. ) 128 { 149 { 129 // "String" moving backwards in CMS, abor 150 // "String" moving backwards in CMS, abort collision !! 130 //G4cout << " abort Collision!! " << G4end 151 //G4cout << " abort Collision!! " << G4endl; 131 return false; 152 return false; 132 } 153 } 133 154 134 toCms.rotateZ(-1*Ptmp.phi()); 155 toCms.rotateZ(-1*Ptmp.phi()); 135 toCms.rotateY(-1*Ptmp.theta()); 156 toCms.rotateY(-1*Ptmp.theta()); 136 157 137 G4LorentzRotation toLab(toCms.inverse()); 158 G4LorentzRotation toLab(toCms.inverse()); 138 159 139 Pprojectile.transform(toCms); 160 Pprojectile.transform(toCms); 140 Ptarget.transform(toCms); 161 Ptarget.transform(toCms); 141 162 142 G4double PZcms2=(S*S+Mprojectile2*Mprojectil << 163 G4double Pt2; 143 2*S*Mprojectile2-2*S*Mtarget2-2*Mproj << 164 G4double ProjMassT2, ProjMassT; >> 165 G4double TargMassT2, TargMassT; >> 166 G4double PZcms2, PZcms; >> 167 G4double PMinusNew, TPlusNew; 144 168 145 if (PZcms2 < 0) {return false;} // It can << 169 G4double S=Psum.mag2(); >> 170 G4double SqrtS=std::sqrt(S); 146 171 147 G4double PZcms = std::sqrt(PZcms2); << 172 if(SqrtS < 2200*MeV) {return false;} // The model cannot work for pp-interactions >> 173 // at Plab < 1.3 GeV/c. Uzhi 148 174 149 if (PutOnMassShell) << 175 PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2- >> 176 2*S*Mprojectile2-2*S*Mtarget2-2*Mprojectile2*Mtarget2)/4./S; >> 177 if(PZcms2 < 0) >> 178 {return false;} // It can be in an interaction with off-shell nuclear nucleon >> 179 >> 180 PZcms = std::sqrt(PZcms2); >> 181 >> 182 if(PutOnMassShell) 150 { 183 { 151 if (Pprojectile.z() > 0.) << 184 if(Pprojectile.z() > 0.) 152 { 185 { 153 Pprojectile.setPz( PZcms); 186 Pprojectile.setPz( PZcms); 154 Ptarget.setPz( -PZcms); 187 Ptarget.setPz( -PZcms); 155 } << 188 } else 156 else << 157 { 189 { 158 Pprojectile.setPz(-PZcms); 190 Pprojectile.setPz(-PZcms); 159 Ptarget.setPz( PZcms); 191 Ptarget.setPz( PZcms); 160 }; 192 }; 161 193 162 Pprojectile.setE(std::sqrt(Mprojectile2+sq << 194 Pprojectile.setE(std::sqrt(Mprojectile2+ 163 Ptarget.setE( std::sqrt( Mtarget2+sq << 195 Pprojectile.x()*Pprojectile.x()+ >> 196 Pprojectile.y()*Pprojectile.y()+ >> 197 PZcms2)); >> 198 Ptarget.setE(std::sqrt( Mtarget2 + >> 199 Ptarget.x()*Ptarget.x()+ >> 200 Ptarget.y()*Ptarget.y()+ >> 201 PZcms2)); 164 } 202 } 165 203 166 G4double maxPtSquare = PZcms2; 204 G4double maxPtSquare = PZcms2; 167 205 168 #ifdef debugDoubleDiffraction << 206 //G4cout << "Pprojectile aft boost : " << Pprojectile << G4endl; 169 G4cout << "Pprojectile after boost to CMS: << 207 //G4cout << "Ptarget aft boost : " << Ptarget << G4endl; 170 G4cout << "Ptarget after boost to CMS: << 208 // G4cout << "cms aft boost : " << (Pprojectile+ Ptarget) << G4endl; 171 #endif << 209 // G4cout << " Projectile Xplus / Xminus : " << 172 << 210 // Pprojectile.plus() << " / " << Pprojectile.minus() << G4endl; 173 G4int PrPDGcode=projectile->GetDefinition << 211 // G4cout << " Target Xplus / Xminus : " << 174 G4int absPrPDGcode=std::abs(PrPDGcode); << 212 // Ptarget.plus() << " / " << Ptarget.minus() << G4endl; 175 G4double MinPrDiffMass(0.); << 176 G4double AveragePt2(0.); << 177 << 178 if (M0projectile <= projectile->GetDefinitio << 179 { // Normal projectile << 180 if( absPrPDGcode > 1000 ) << 181 { << 182 if ( absPrPDGcode > 4000 && absPrPDGcode << 183 { << 184 MinPrDiffMass = projectile->GetDefinit << 185 AveragePt2 = 0.3; << 186 } << 187 else << 188 { << 189 MinPrDiffMass = 1.16; // << 190 AveragePt2 = 0.3; << 191 } << 192 } << 193 else if( absPrPDGcode == 211 || PrPDGcode << 194 { << 195 MinPrDiffMass = 1.0; << 196 AveragePt2 = 0.3; << 197 } << 198 else if( absPrPDGcode == 321 || absPrPDGco << 199 { << 200 MinPrDiffMass = 1.1; << 201 AveragePt2 = 0.3; << 202 } << 203 else if( absPrPDGcode > 400 && absPrPDGcod << 204 { << 205 MinPrDiffMass = projectile->GetDefinitio << 206 AveragePt2 = 0.3; << 207 } << 208 else << 209 { << 210 MinPrDiffMass = 1.16; << 211 AveragePt2 = 0.3; << 212 } << 213 } << 214 else << 215 { // Excited projectile << 216 MinPrDiffMass = M0projectile + 220.0*MeV; << 217 AveragePt2 = 0.3; << 218 } << 219 << 220 MinPrDiffMass = MinPrDiffMass * GeV; << 221 AveragePt2 = AveragePt2 * GeV*GeV; << 222 //------------------------------------------ << 223 G4double MinTrDiffMass = 1.16*GeV; << 224 << 225 if (SqrtS < MinPrDiffMass + MinTrDiffMass) { << 226 << 227 G4double MinPrDiffMass2 = MinPrDiffMass * Mi << 228 G4double MinTrDiffMass2 = MinTrDiffMass * Mi << 229 << 230 G4double Pt2; << 231 G4double ProjMassT2, ProjMassT; << 232 G4double TargMassT2, TargMassT; << 233 G4double PMinusNew, TPlusNew; << 234 213 235 G4LorentzVector Qmomentum; 214 G4LorentzVector Qmomentum; 236 G4double Qminus, Qplus; 215 G4double Qminus, Qplus; 237 216 238 G4int whilecount=0; 217 G4int whilecount=0; 239 do { 218 do { >> 219 // Generate pt >> 220 240 if (whilecount++ >= 500 && (whilecount%100 221 if (whilecount++ >= 500 && (whilecount%100)==0) 241 if (whilecount > 1000 ) { << 222 //G4cout << "G4QGSDiffractiveExcitation::ExciteParticipants possibly looping" 242 Qmomentum=G4LorentzVector(0.,0.,0.,0.); << 223 // << ", loop count/ maxPtSquare : " 243 return false; // Ignore this interactio << 224 // << whilecount << " / " << maxPtSquare << G4endl; 244 } << 225 >> 226 if (whilecount > 1000 ) >> 227 { >> 228 Qmomentum=G4LorentzVector(0.,0.,0.,0.); >> 229 return false; // Ignore this interaction >> 230 } 245 231 246 // Generate pt << 247 Qmomentum=G4LorentzVector(GaussianPt(Avera 232 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0); 248 233 >> 234 //G4cout << "generated Pt " << Qmomentum << G4endl; >> 235 //G4cout << "Pprojectile with pt : " << Pprojectile+Qmomentum << G4endl; >> 236 //G4cout << "Ptarget with pt : " << Ptarget-Qmomentum << G4endl; >> 237 >> 238 // Momentum transfer >> 239 /* >> 240 G4double Xmin = minmass / ( Pprojectile.e() + Ptarget.e() ); >> 241 G4double Xmax=1.; >> 242 G4double Xplus =ChooseX(Xmin,Xmax); >> 243 G4double Xminus=ChooseX(Xmin,Xmax); >> 244 >> 245 //G4cout << " X-plus " << Xplus << G4endl; >> 246 //G4cout << " X-minus " << Xminus << G4endl; >> 247 >> 248 G4double pt2=G4ThreeVector(Qmomentum.vect()).mag2(); >> 249 G4double Qplus =-1 * pt2 / Xminus/Ptarget.minus(); >> 250 G4double Qminus= pt2 / Xplus /Pprojectile.plus(); >> 251 */ >> 252 249 Pt2=G4ThreeVector(Qmomentum.vect()).mag2() 253 Pt2=G4ThreeVector(Qmomentum.vect()).mag2(); 250 ProjMassT2=MinPrDiffMass2+Pt2; << 254 ProjMassT2=Mprojectile2+Pt2; 251 ProjMassT =std::sqrt(ProjMassT2); 255 ProjMassT =std::sqrt(ProjMassT2); 252 256 253 TargMassT2=MinTrDiffMass2+Pt2; << 257 TargMassT2=Mtarget2+Pt2; 254 TargMassT =std::sqrt(TargMassT2); 258 TargMassT =std::sqrt(TargMassT2); 255 259 256 if (SqrtS < ProjMassT + TargMassT) continu << 260 PZcms2=(S*S+ProjMassT2*ProjMassT2+ 257 << 261 TargMassT2*TargMassT2- 258 PZcms2=(S*S+ProjMassT2*ProjMassT2+TargMass << 262 2.*S*ProjMassT2-2.*S*TargMassT2- 259 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjM << 263 2.*ProjMassT2*TargMassT2)/4./S; 260 if (PZcms2 < 0 ) {PZcms2=0;}; << 264 if(PZcms2 < 0 ) {PZcms2=0;}; 261 PZcms =std::sqrt(PZcms2); 265 PZcms =std::sqrt(PZcms2); 262 266 263 G4double PMinusMin=std::sqrt(ProjMassT2+PZ 267 G4double PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms; 264 G4double PMinusMax=SqrtS-TargMassT; 268 G4double PMinusMax=SqrtS-TargMassT; 265 269 266 PMinusNew=ChooseP(PMinusMin,PMinusMax); 270 PMinusNew=ChooseP(PMinusMin,PMinusMax); 267 Qminus=PMinusNew-Pprojectile.minus(); 271 Qminus=PMinusNew-Pprojectile.minus(); 268 272 269 G4double TPlusMin=std::sqrt(TargMassT2+PZc 273 G4double TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms; 270 G4double TPlusMax=SqrtS-ProjMassT; 274 G4double TPlusMax=SqrtS-ProjMassT; 271 275 272 TPlusNew=ChooseP(TPlusMin, TPlusMax); 276 TPlusNew=ChooseP(TPlusMin, TPlusMax); 273 Qplus=-(TPlusNew-Ptarget.plus()); 277 Qplus=-(TPlusNew-Ptarget.plus()); 274 278 275 Qmomentum.setPz( (Qplus-Qminus)/2 ); 279 Qmomentum.setPz( (Qplus-Qminus)/2 ); 276 Qmomentum.setE( (Qplus+Qminus)/2 ); 280 Qmomentum.setE( (Qplus+Qminus)/2 ); 277 281 278 } while ( (Pprojectile+Qmomentum).mag2() < << 282 //G4cout << "Qplus / Qminus " << Qplus << " / " << Qminus<<G4endl; 279 (Ptarget -Qmomentum).mag2() < MinTrD << 283 //G4cout << "pt2" << pt2 << G4endl; >> 284 //G4cout << "Qmomentum " << Qmomentum << G4endl; >> 285 //G4cout << " Masses (P/T) : " << (Pprojectile+Qmomentum).mag() << >> 286 // " / " << (Ptarget-Qmomentum).mag() << G4endl; >> 287 /* >> 288 } while ( (Pprojectile+Qmomentum).mag2() <= Mprojectile2 || >> 289 (Ptarget-Qmomentum).mag2() <= Mtarget2 ); >> 290 */ >> 291 } while ( /* Loop checking, 26.10.2015, A.Ribon */ >> 292 ( (Pprojectile+Qmomentum).mag2() < Mprojectile2 || // No without excitation >> 293 (Ptarget -Qmomentum).mag2() < Mtarget2 ) || >> 294 ( (Pprojectile+Qmomentum).mag2() < ProjectileDiffCut2 && // No double Diffraction >> 295 (Ptarget -Qmomentum).mag2() < NuclearNucleonDiffCut2) ); >> 296 >> 297 if((Ptarget-Qmomentum).mag2() < NuclearNucleonDiffCut2) // Uzhi Projectile diffraction >> 298 { >> 299 G4double TMinusNew=SqrtS-PMinusNew; >> 300 Qminus=Ptarget.minus()-TMinusNew; >> 301 TPlusNew=TargMassT2/TMinusNew; >> 302 Qplus=Ptarget.plus()-TPlusNew; >> 303 >> 304 Qmomentum.setPz( (Qplus-Qminus)/2 ); >> 305 Qmomentum.setE( (Qplus+Qminus)/2 ); >> 306 } >> 307 else if((Pprojectile+Qmomentum).mag2() < ProjectileDiffCut2) // Uzhi Target diffraction >> 308 { >> 309 G4double PPlusNew=SqrtS-TPlusNew; >> 310 Qplus=PPlusNew-Pprojectile.plus(); >> 311 PMinusNew=ProjMassT2/PPlusNew; >> 312 Qminus=PMinusNew-Pprojectile.minus(); >> 313 >> 314 Qmomentum.setPz( (Qplus-Qminus)/2 ); >> 315 Qmomentum.setE( (Qplus+Qminus)/2 ); >> 316 }; 280 317 281 Pprojectile += Qmomentum; 318 Pprojectile += Qmomentum; 282 Ptarget -= Qmomentum; 319 Ptarget -= Qmomentum; 283 320 >> 321 // Vova >> 322 >> 323 /* >> 324 Pprojectile.setPz(0.); >> 325 Pprojectile.setE(SqrtS-M0target); >> 326 Ptarget.setPz(0.); >> 327 Ptarget.setE(M0target); >> 328 */ >> 329 >> 330 //G4cout << "Pprojectile with Q : " << Pprojectile << G4endl; >> 331 //G4cout << "Ptarget with Q : " << Ptarget << G4endl; >> 332 //G4cout << "Projectile back: " << toLab * Pprojectile << G4endl; >> 333 //G4cout << "Target back: " << toLab * Ptarget << G4endl; >> 334 284 // Transform back and update SplitableHadron 335 // Transform back and update SplitableHadron Participant. 285 Pprojectile.transform(toLab); 336 Pprojectile.transform(toLab); 286 Ptarget.transform(toLab); 337 Ptarget.transform(toLab); 287 338 288 #ifdef debugDoubleDiffraction << 339 //G4cout << "Pprojectile with Q M: " << Pprojectile<<" "<< Pprojectile.mag() << G4endl; 289 G4cout << "Pprojectile after boost to Lab: << 340 //G4cout << "Ptarget with Q M: " << Ptarget <<" "<< Ptarget.mag() << G4endl; 290 G4cout << "Ptarget after boost to Lab: << 341 //G4cout << "Target mass " << Ptarget.mag() << G4endl; 291 #endif << 292 342 293 target->Set4Momentum(Ptarget); 343 target->Set4Momentum(Ptarget); >> 344 >> 345 //G4cout << "Projectile mass " << Pprojectile.mag() << G4endl; >> 346 294 projectile->Set4Momentum(Pprojectile); 347 projectile->Set4Momentum(Pprojectile); 295 348 296 return true; 349 return true; 297 } 350 } 298 351 299 352 300 G4ExcitedString * G4QGSDiffractiveExcitation:: 353 G4ExcitedString * G4QGSDiffractiveExcitation:: 301 String(G4VSplitableHadron * hadron, G4bool isP 354 String(G4VSplitableHadron * hadron, G4bool isProjectile) const 302 { 355 { 303 hadron->SplitUp(); 356 hadron->SplitUp(); 304 G4Parton *start= hadron->GetNextParton(); 357 G4Parton *start= hadron->GetNextParton(); 305 if ( start==NULL) << 358 if ( start==NULL) { 306 { G4cout << " G4QGSDiffractiveExcitation::St << 359 G4cout << " G4FTFModel::String() Error:No start parton found"<< G4endl; 307 return NULL; 360 return NULL; 308 } 361 } 309 G4Parton *end = hadron->GetNextParton(); 362 G4Parton *end = hadron->GetNextParton(); 310 if ( end==NULL) << 363 if ( end==NULL) { 311 { G4cout << " G4QGSDiffractiveExcitation::St << 364 G4cout << " G4FTFModel::String() Error:No end parton found"<< G4endl; 312 return NULL; 365 return NULL; 313 } 366 } 314 367 315 G4ExcitedString * string; 368 G4ExcitedString * string; 316 if ( isProjectile ) << 369 if ( isProjectile ) { 317 { << 318 string= new G4ExcitedString(end,start, +1) 370 string= new G4ExcitedString(end,start, +1); 319 } else { 371 } else { 320 string= new G4ExcitedString(start,end, -1) 372 string= new G4ExcitedString(start,end, -1); 321 } 373 } 322 374 323 string->SetPosition(hadron->GetPosition()); 375 string->SetPosition(hadron->GetPosition()); 324 376 325 // momenta of string ends 377 // momenta of string ends >> 378 G4double ptSquared= hadron->Get4Momentum().perp2(); >> 379 G4double transverseMassSquared= hadron->Get4Momentum().plus() >> 380 * hadron->Get4Momentum().minus(); 326 381 327 G4double maxAvailMomentumSquared=sqr(hadron- << 382 G4double maxAvailMomentumSquared= >> 383 sqr( std::sqrt(transverseMassSquared) - std::sqrt(ptSquared) ); 328 384 329 G4double widthOfPtSquare = 0.5*sqr(GeV); << 385 G4double widthOfPtSquare = 0.25; // Uzhi <Pt^2>=0.25 ??? 330 G4ThreeVector pt=GaussianPt(widthOfPtSquare, 386 G4ThreeVector pt=GaussianPt(widthOfPtSquare,maxAvailMomentumSquared); 331 387 332 G4LorentzVector Pstart(G4LorentzVector(pt,0. 388 G4LorentzVector Pstart(G4LorentzVector(pt,0.)); 333 G4LorentzVector Pend; 389 G4LorentzVector Pend; 334 Pend.setPx(hadron->Get4Momentum().px() - pt. 390 Pend.setPx(hadron->Get4Momentum().px() - pt.x()); 335 Pend.setPy(hadron->Get4Momentum().py() - pt. 391 Pend.setPy(hadron->Get4Momentum().py() - pt.y()); 336 392 337 G4double tm1=hadron->Get4Momentum().minus() 393 G4double tm1=hadron->Get4Momentum().minus() + 338 ( Pend.perp2()-Pstart.perp2() ) / had 394 ( Pend.perp2()-Pstart.perp2() ) / hadron->Get4Momentum().plus(); 339 395 340 G4double tm2= std::sqrt( std::max(0., sqr(tm 396 G4double tm2= std::sqrt( std::max(0., sqr(tm1) - 341 4. * Pend.perp2() * hadron->Get4Momen 397 4. * Pend.perp2() * hadron->Get4Momentum().minus() 342 / hadron->Get4Momentum().plus() )); 398 / hadron->Get4Momentum().plus() )); 343 399 344 G4int Sign= isProjectile ? -1 : 1; 400 G4int Sign= isProjectile ? -1 : 1; 345 401 346 G4double endMinus = 0.5 * (tm1 + Sign*tm2); 402 G4double endMinus = 0.5 * (tm1 + Sign*tm2); 347 G4double startMinus= hadron->Get4Momentum(). 403 G4double startMinus= hadron->Get4Momentum().minus() - endMinus; 348 404 349 G4double startPlus= Pstart.perp2() / startM 405 G4double startPlus= Pstart.perp2() / startMinus; 350 G4double endPlus = hadron->Get4Momentum().p 406 G4double endPlus = hadron->Get4Momentum().plus() - startPlus; 351 407 352 Pstart.setPz(0.5*(startPlus - startMinus)); 408 Pstart.setPz(0.5*(startPlus - startMinus)); 353 Pstart.setE(0.5*(startPlus + startMinus)); 409 Pstart.setE(0.5*(startPlus + startMinus)); 354 410 355 Pend.setPz(0.5*(endPlus - endMinus)); 411 Pend.setPz(0.5*(endPlus - endMinus)); 356 Pend.setE(0.5*(endPlus + endMinus)); 412 Pend.setE(0.5*(endPlus + endMinus)); 357 413 358 start->Set4Momentum(Pstart); 414 start->Set4Momentum(Pstart); 359 end->Set4Momentum(Pend); 415 end->Set4Momentum(Pend); 360 416 361 #ifdef debugQGSdiffExictation << 417 #ifdef G4_FTFDEBUG 362 G4cout << " generated string flavors << 418 G4cout << " generated string flavors " << start->GetPDGcode() << " / " << end->GetPDGcode() << G4endl; 363 G4cout << " generated string momenta: qu << 419 G4cout << " generated string momenta: quark " << start->Get4Momentum() << "mass : " <<start->Get4Momentum().mag()<< G4endl; 364 G4cout << " generated string momenta: Diqu << 420 G4cout << " generated string momenta: Diquark " << end ->Get4Momentum() << "mass : " <<end->Get4Momentum().mag()<< G4endl; 365 G4cout << " sum of ends << 421 G4cout << " sum of ends " << Pstart+Pend << G4endl; 366 G4cout << " Original << 422 G4cout << " Original " << hadron->Get4Momentum() << G4endl; 367 #endif 423 #endif 368 424 369 return string; 425 return string; 370 } 426 } 371 427 372 428 373 // --------- private methods ----------------- 429 // --------- private methods ---------------------- 374 430 375 G4double G4QGSDiffractiveExcitation::ChooseP(G 431 G4double G4QGSDiffractiveExcitation::ChooseP(G4double Pmin, G4double Pmax) const 376 { 432 { 377 // choose an x between Xmin and Xmax with P( 433 // choose an x between Xmin and Xmax with P(x) ~ 1/x 378 // to be improved... 434 // to be improved... 379 435 380 G4double range=Pmax-Pmin; 436 G4double range=Pmax-Pmin; 381 437 382 if ( Pmin <= 0. || range <=0. ) 438 if ( Pmin <= 0. || range <=0. ) 383 { 439 { 384 G4cout << " Pmin, range : " << Pmin << " , 440 G4cout << " Pmin, range : " << Pmin << " , " << range << G4endl; 385 throw G4HadronicException(__FILE__, __LINE 441 throw G4HadronicException(__FILE__, __LINE__, "G4QGSDiffractiveExcitation::ChooseP : Invalid arguments "); 386 } 442 } 387 443 388 G4double P; 444 G4double P; >> 445 /* >> 446 do { >> 447 x=Xmin + G4UniformRand() * range; >> 448 } while ( Xmin/x < G4UniformRand() ); >> 449 */ >> 450 389 P=Pmin * G4Pow::GetInstance()->powA(Pmax/Pmi 451 P=Pmin * G4Pow::GetInstance()->powA(Pmax/Pmin,G4UniformRand()); >> 452 390 //debug-hpw cout << "DiffractiveX "<<x<<G4en 453 //debug-hpw cout << "DiffractiveX "<<x<<G4endl; 391 return P; 454 return P; 392 } 455 } 393 456 394 << 395 G4ThreeVector G4QGSDiffractiveExcitation::Gaus 457 G4ThreeVector G4QGSDiffractiveExcitation::GaussianPt(G4double AveragePt2, G4double maxPtSquare) const 396 { // @@ this method is used in FTF 458 { // @@ this method is used in FTFModel as well. Should go somewhere common! >> 459 397 G4double Pt2; 460 G4double Pt2; >> 461 /* >> 462 do { >> 463 pt2=widthSquare * G4Log( G4UniformRand() ); >> 464 } while ( pt2 > maxPtSquare); >> 465 */ 398 466 399 Pt2 = -AveragePt2 * G4Log(1. + G4UniformRand 467 Pt2 = -AveragePt2 * G4Log(1. + G4UniformRand() * (G4Exp(-maxPtSquare/AveragePt2)-1.)); 400 468 401 G4double Pt=std::sqrt(Pt2); 469 G4double Pt=std::sqrt(Pt2); 402 470 403 G4double phi=G4UniformRand() * twopi; 471 G4double phi=G4UniformRand() * twopi; 404 472 405 return G4ThreeVector (Pt*std::cos(phi), Pt*s 473 return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.); 406 } 474 } >> 475 407 476