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,v 1.1 2007/05/25 07:30:47 gunter Exp $ 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 "G4SystemOfUnits.hh" << 48 #include "Randomize.hh" 47 #include "Randomize.hh" 49 48 50 #include "G4QGSDiffractiveExcitation.hh" 49 #include "G4QGSDiffractiveExcitation.hh" 51 #include "G4LorentzRotation.hh" 50 #include "G4LorentzRotation.hh" 52 #include "G4ThreeVector.hh" 51 #include "G4ThreeVector.hh" 53 #include "G4ParticleDefinition.hh" 52 #include "G4ParticleDefinition.hh" 54 #include "G4VSplitableHadron.hh" 53 #include "G4VSplitableHadron.hh" 55 #include "G4ExcitedString.hh" 54 #include "G4ExcitedString.hh" 56 //#include "G4ios.hh" 55 //#include "G4ios.hh" 57 56 58 #include "G4Exp.hh" << 57 G4QGSDiffractiveExcitation::G4QGSDiffractiveExcitation() // Uzhi 59 #include "G4Log.hh" << 60 #include "G4Pow.hh" << 61 << 62 //============================================ << 63 << 64 //#define debugDoubleDiffraction << 65 << 66 //============================================ << 67 << 68 G4QGSDiffractiveExcitation::G4QGSDiffractiveEx << 69 { 58 { 70 } 59 } 71 60 72 G4QGSDiffractiveExcitation::~G4QGSDiffractiveE << 61 G4bool G4QGSDiffractiveExcitation:: >> 62 ExciteParticipants(G4VSplitableHadron *projectile, G4VSplitableHadron *target) const 73 { 63 { 74 } << 75 64 >> 65 G4LorentzVector Pprojectile=projectile->Get4Momentum(); 76 66 77 G4bool G4QGSDiffractiveExcitation:: << 67 // -------------------- Projectile parameters ----------------------------------- 78 ExciteParticipants(G4VSplitableHadron *project << 68 G4bool PutOnMassShell=0; 79 { << 69 80 #ifdef debugDoubleDiffraction << 70 // G4double M0projectile=projectile->GetDefinition()->GetPDGMass(); // With de-excitation 81 G4cout<<G4endl<<"G4QGSDiffractiveExcitatio << 71 G4double M0projectile = Pprojectile.mag(); // Without de-excitation 82 G4cout<<"Proj Targ "<<projectile->GetDefin << 72 83 G4cout<<"Proj 4 Mom "<<projectile->Get4Mom << 73 if(M0projectile < projectile->GetDefinition()->GetPDGMass()) 84 G4cout<<"Targ 4 Mom "<<target->Get4Momentu << 74 { 85 #endif << 75 PutOnMassShell=1; 86 << 76 M0projectile=projectile->GetDefinition()->GetPDGMass(); 87 G4LorentzVector Pprojectile=projectile->Get4 << 77 } 88 << 78 89 // -------------------- Projectile parameter << 79 G4double Mprojectile2 = M0projectile * M0projectile; 90 G4bool PutOnMassShell=0; << 80 91 << 81 G4int PDGcode=projectile->GetDefinition()->GetPDGEncoding(); 92 G4double M0projectile = Pprojectile.mag(); << 82 G4int absPDGcode=std::abs(PDGcode); 93 << 83 G4double ProjectileDiffCut; 94 if(M0projectile < projectile->GetDefinition( << 84 G4double AveragePt2; 95 { << 85 96 PutOnMassShell=1; << 86 if( absPDGcode > 1000 ) //------Projectile is baryon -------- 97 M0projectile=projectile->GetDefinition()-> << 87 { 98 } << 88 ProjectileDiffCut = 1.1; // GeV 99 << 89 AveragePt2 = 0.3; // GeV^2 100 // -------------------- Target parameters -- << 90 } 101 G4LorentzVector Ptarget=target->Get4Momentum << 91 else if( absPDGcode == 211 || PDGcode == 111) //------Projectile is Pion ----------- 102 << 92 { 103 G4double M0target = Ptarget.mag(); << 93 ProjectileDiffCut = 1.0; // GeV 104 << 94 AveragePt2 = 0.3; // GeV^2 105 if(M0target < target->GetDefinition()->GetPD << 95 } 106 { << 96 else if( absPDGcode == 321 || PDGcode == -311) //------Projectile is Kaon ----------- 107 PutOnMassShell=1; << 97 { 108 M0target=target->GetDefinition()->GetPDGMa << 98 ProjectileDiffCut = 1.1; // GeV 109 } << 99 AveragePt2 = 0.3; // GeV^2 110 << 100 } 111 G4LorentzVector Psum=Pprojectile+Ptarget; << 101 else //------Projectile is undefined, Nucleon assumed 112 G4double S=Psum.mag2(); << 102 { 113 G4double SqrtS=std::sqrt(S); << 103 ProjectileDiffCut = 1.1; // GeV 114 << 104 AveragePt2 = 0.3; // GeV^2 115 if (SqrtS < M0projectile + M0target) {return << 105 }; 116 << 106 117 << 107 ProjectileDiffCut = ProjectileDiffCut * GeV; 118 G4double Mprojectile2 = M0projectile * M0pro << 108 AveragePt2 = AveragePt2 * GeV*GeV; 119 G4double Mtarget2 = M0target * M0tar << 109 120 << 110 // -------------------- Target parameters ---------------------------------------------- 121 // Transform momenta to cms and then rotate << 111 G4LorentzVector Ptarget=target->Get4Momentum(); 122 << 112 123 G4LorentzRotation toCms(-1*Psum.boostVector( << 113 G4double M0target = Ptarget.mag(); 124 << 114 125 G4LorentzVector Ptmp=toCms*Pprojectile; << 115 if(M0target < target->GetDefinition()->GetPDGMass()) 126 << 116 { 127 if ( Ptmp.pz() <= 0. ) << 117 PutOnMassShell=1; 128 { << 118 M0target=target->GetDefinition()->GetPDGMass(); 129 // "String" moving backwards in CMS, abor << 119 } 130 //G4cout << " abort Collision!! " << G4end << 120 131 return false; << 121 G4double Mtarget2 = M0target * M0target; //Ptarget.mag2(); // for AA-inter. 132 } << 122 133 << 123 G4double NuclearNucleonDiffCut = 1.1*GeV; 134 toCms.rotateZ(-1*Ptmp.phi()); << 124 135 toCms.rotateY(-1*Ptmp.theta()); << 125 G4double ProjectileDiffCut2 = ProjectileDiffCut * ProjectileDiffCut; 136 << 126 G4double NuclearNucleonDiffCut2 = NuclearNucleonDiffCut * NuclearNucleonDiffCut; 137 G4LorentzRotation toLab(toCms.inverse()); << 127 138 << 128 // Transform momenta to cms and then rotate parallel to z axis; 139 Pprojectile.transform(toCms); << 129 140 Ptarget.transform(toCms); << 130 G4LorentzVector Psum; 141 << 131 Psum=Pprojectile+Ptarget; 142 G4double PZcms2=(S*S+Mprojectile2*Mprojectil << 132 143 2*S*Mprojectile2-2*S*Mtarget2-2*Mproj << 133 G4LorentzRotation toCms(-1*Psum.boostVector()); 144 << 134 145 if (PZcms2 < 0) {return false;} // It can << 135 G4LorentzVector Ptmp=toCms*Pprojectile; 146 << 136 147 G4double PZcms = std::sqrt(PZcms2); << 137 if ( Ptmp.pz() <= 0. ) 148 << 138 { 149 if (PutOnMassShell) << 139 // "String" moving backwards in CMS, abort collision !! 150 { << 140 //G4cout << " abort Collision!! " << G4endl; 151 if (Pprojectile.z() > 0.) << 141 return false; 152 { << 142 } 153 Pprojectile.setPz( PZcms); << 143 154 Ptarget.setPz( -PZcms); << 144 toCms.rotateZ(-1*Ptmp.phi()); 155 } << 145 toCms.rotateY(-1*Ptmp.theta()); 156 else << 146 157 { << 147 G4LorentzRotation toLab(toCms.inverse()); 158 Pprojectile.setPz(-PZcms); << 148 159 Ptarget.setPz( PZcms); << 149 Pprojectile.transform(toCms); 160 }; << 150 Ptarget.transform(toCms); 161 << 151 162 Pprojectile.setE(std::sqrt(Mprojectile2+sq << 152 G4double Pt2; 163 Ptarget.setE( std::sqrt( Mtarget2+sq << 153 G4double ProjMassT2, ProjMassT; 164 } << 154 G4double TargMassT2, TargMassT; 165 << 155 G4double PZcms2, PZcms; 166 G4double maxPtSquare = PZcms2; << 156 G4double PMinusNew, TPlusNew; 167 << 157 168 #ifdef debugDoubleDiffraction << 158 G4double S=Psum.mag2(); 169 G4cout << "Pprojectile after boost to CMS: << 159 G4double SqrtS=std::sqrt(S); 170 G4cout << "Ptarget after boost to CMS: << 160 171 #endif << 161 if(SqrtS < 2200*MeV) {return false;} // The model cannot work for pp-interactions 172 << 162 // at Plab < 1.3 GeV/c. Uzhi 173 G4int PrPDGcode=projectile->GetDefinition << 163 174 G4int absPrPDGcode=std::abs(PrPDGcode); << 164 PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2- 175 G4double MinPrDiffMass(0.); << 165 2*S*Mprojectile2-2*S*Mtarget2-2*Mprojectile2*Mtarget2)/4./S; 176 G4double AveragePt2(0.); << 166 if(PZcms2 < 0) 177 << 167 {return false;} // It can be in an interaction with off-shell nuclear nucleon 178 if (M0projectile <= projectile->GetDefinitio << 168 179 { // Normal projectile << 169 PZcms = std::sqrt(PZcms2); 180 if( absPrPDGcode > 1000 ) << 170 181 { << 171 if(PutOnMassShell) 182 if ( absPrPDGcode > 4000 && absPrPDGcode << 172 { 183 { << 173 if(Pprojectile.z() > 0.) 184 MinPrDiffMass = projectile->GetDefinit << 174 { 185 AveragePt2 = 0.3; << 175 Pprojectile.setPz( PZcms); 186 } << 176 Ptarget.setPz( -PZcms); 187 else << 177 } 188 { << 178 else 189 MinPrDiffMass = 1.16; // << 179 { 190 AveragePt2 = 0.3; << 180 Pprojectile.setPz(-PZcms); 191 } << 181 Ptarget.setPz( PZcms); 192 } << 182 }; 193 else if( absPrPDGcode == 211 || PrPDGcode << 183 194 { << 184 Pprojectile.setE(std::sqrt(Mprojectile2+ 195 MinPrDiffMass = 1.0; << 185 Pprojectile.x()*Pprojectile.x()+ 196 AveragePt2 = 0.3; << 186 Pprojectile.y()*Pprojectile.y()+ 197 } << 187 PZcms2)); 198 else if( absPrPDGcode == 321 || absPrPDGco << 188 Ptarget.setE(std::sqrt( Mtarget2 + 199 { << 189 Ptarget.x()*Ptarget.x()+ 200 MinPrDiffMass = 1.1; << 190 Ptarget.y()*Ptarget.y()+ 201 AveragePt2 = 0.3; << 191 PZcms2)); 202 } << 192 } 203 else if( absPrPDGcode > 400 && absPrPDGcod << 193 204 { << 194 G4double maxPtSquare = PZcms2; 205 MinPrDiffMass = projectile->GetDefinitio << 195 206 AveragePt2 = 0.3; << 196 //G4cout << "Pprojectile aft boost : " << Pprojectile << G4endl; 207 } << 197 //G4cout << "Ptarget aft boost : " << Ptarget << G4endl; 208 else << 198 // G4cout << "cms aft boost : " << (Pprojectile+ Ptarget) << G4endl; 209 { << 199 210 MinPrDiffMass = 1.16; << 200 // G4cout << " Projectile Xplus / Xminus : " << 211 AveragePt2 = 0.3; << 201 // Pprojectile.plus() << " / " << Pprojectile.minus() << G4endl; 212 } << 202 // G4cout << " Target Xplus / Xminus : " << 213 } << 203 // Ptarget.plus() << " / " << Ptarget.minus() << G4endl; 214 else << 204 215 { // Excited projectile << 205 G4LorentzVector Qmomentum; 216 MinPrDiffMass = M0projectile + 220.0*MeV; << 206 G4double Qminus, Qplus; 217 AveragePt2 = 0.3; << 207 218 } << 208 // /* Vova 219 << 209 G4int whilecount=0; 220 MinPrDiffMass = MinPrDiffMass * GeV; << 210 do { 221 AveragePt2 = AveragePt2 * GeV*GeV; << 211 // Generate pt 222 //------------------------------------------ << 212 223 G4double MinTrDiffMass = 1.16*GeV; << 213 if (whilecount++ >= 500 && (whilecount%100)==0) 224 << 214 // G4cout << "G4QGSDiffractiveExcitation::ExciteParticipants possibly looping" 225 if (SqrtS < MinPrDiffMass + MinTrDiffMass) { << 215 // << ", loop count/ maxPtSquare : " 226 << 216 // << whilecount << " / " << maxPtSquare << G4endl; 227 G4double MinPrDiffMass2 = MinPrDiffMass * Mi << 217 if (whilecount > 1000 ) 228 G4double MinTrDiffMass2 = MinTrDiffMass * Mi << 218 { 229 << 219 Qmomentum=G4LorentzVector(0.,0.,0.,0.); 230 G4double Pt2; << 220 return false; // Ignore this interaction 231 G4double ProjMassT2, ProjMassT; << 221 } 232 G4double TargMassT2, TargMassT; << 222 233 G4double PMinusNew, TPlusNew; << 223 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0); 234 << 224 235 G4LorentzVector Qmomentum; << 225 //G4cout << "generated Pt " << Qmomentum << G4endl; 236 G4double Qminus, Qplus; << 226 //G4cout << "Pprojectile with pt : " << Pprojectile+Qmomentum << G4endl; 237 << 227 //G4cout << "Ptarget with pt : " << Ptarget-Qmomentum << G4endl; 238 G4int whilecount=0; << 228 239 do { << 229 // Momentum transfer 240 if (whilecount++ >= 500 && (whilecount%100 << 230 /* // Uzhi 241 if (whilecount > 1000 ) { << 231 G4double Xmin = minmass / ( Pprojectile.e() + Ptarget.e() ); 242 Qmomentum=G4LorentzVector(0.,0.,0.,0.); << 232 G4double Xmax=1.; 243 return false; // Ignore this interactio << 233 G4double Xplus =ChooseX(Xmin,Xmax); 244 } << 234 G4double Xminus=ChooseX(Xmin,Xmax); 245 << 235 246 // Generate pt << 236 // G4cout << " X-plus " << Xplus << G4endl; 247 Qmomentum=G4LorentzVector(GaussianPt(Avera << 237 // G4cout << " X-minus " << Xminus << G4endl; 248 << 238 249 Pt2=G4ThreeVector(Qmomentum.vect()).mag2() << 239 G4double pt2=G4ThreeVector(Qmomentum.vect()).mag2(); 250 ProjMassT2=MinPrDiffMass2+Pt2; << 240 G4double Qplus =-1 * pt2 / Xminus/Ptarget.minus(); 251 ProjMassT =std::sqrt(ProjMassT2); << 241 G4double Qminus= pt2 / Xplus /Pprojectile.plus(); 252 << 242 */ // Uzhi * 253 TargMassT2=MinTrDiffMass2+Pt2; << 243 254 TargMassT =std::sqrt(TargMassT2); << 244 Pt2=G4ThreeVector(Qmomentum.vect()).mag2(); 255 << 245 ProjMassT2=Mprojectile2+Pt2; 256 if (SqrtS < ProjMassT + TargMassT) continu << 246 ProjMassT =std::sqrt(ProjMassT2); 257 << 247 258 PZcms2=(S*S+ProjMassT2*ProjMassT2+TargMass << 248 TargMassT2=Mtarget2+Pt2; 259 2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjM << 249 TargMassT =std::sqrt(TargMassT2); 260 if (PZcms2 < 0 ) {PZcms2=0;}; << 250 261 PZcms =std::sqrt(PZcms2); << 251 PZcms2=(S*S+ProjMassT2*ProjMassT2+ 262 << 252 TargMassT2*TargMassT2- 263 G4double PMinusMin=std::sqrt(ProjMassT2+PZ << 253 2.*S*ProjMassT2-2.*S*TargMassT2- 264 G4double PMinusMax=SqrtS-TargMassT; << 254 2.*ProjMassT2*TargMassT2)/4./S; 265 << 255 if(PZcms2 < 0 ) {PZcms2=0;}; 266 PMinusNew=ChooseP(PMinusMin,PMinusMax); << 256 PZcms =std::sqrt(PZcms2); 267 Qminus=PMinusNew-Pprojectile.minus(); << 257 268 << 258 G4double PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms; 269 G4double TPlusMin=std::sqrt(TargMassT2+PZc << 259 G4double PMinusMax=SqrtS-TargMassT; 270 G4double TPlusMax=SqrtS-ProjMassT; << 260 271 << 261 PMinusNew=ChooseP(PMinusMin,PMinusMax); 272 TPlusNew=ChooseP(TPlusMin, TPlusMax); << 262 Qminus=PMinusNew-Pprojectile.minus(); 273 Qplus=-(TPlusNew-Ptarget.plus()); << 263 274 << 264 G4double TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms; 275 Qmomentum.setPz( (Qplus-Qminus)/2 ); << 265 G4double TPlusMax=SqrtS-ProjMassT; 276 Qmomentum.setE( (Qplus+Qminus)/2 ); << 266 277 << 267 TPlusNew=ChooseP(TPlusMin, TPlusMax); 278 } while ( (Pprojectile+Qmomentum).mag2() < << 268 Qplus=-(TPlusNew-Ptarget.plus()); 279 (Ptarget -Qmomentum).mag2() < MinTrD << 269 280 << 270 Qmomentum.setPz( (Qplus-Qminus)/2 ); 281 Pprojectile += Qmomentum; << 271 Qmomentum.setE( (Qplus+Qminus)/2 ); 282 Ptarget -= Qmomentum; << 272 283 << 273 //G4cout << "Qplus / Qminus " << Qplus << " / " << Qminus<<G4endl; 284 // Transform back and update SplitableHadron << 274 // G4cout << "pt2" << pt2 << G4endl; 285 Pprojectile.transform(toLab); << 275 // G4cout << "Qmomentum " << Qmomentum << G4endl; 286 Ptarget.transform(toLab); << 276 // G4cout << " Masses (P/T) : " << (Pprojectile+Qmomentum).mag() << 287 << 277 // " / " << (Ptarget-Qmomentum).mag() << G4endl; 288 #ifdef debugDoubleDiffraction << 278 /* // Uzhi 289 G4cout << "Pprojectile after boost to Lab: << 279 } while ( (Pprojectile+Qmomentum).mag2() <= Mprojectile2 || 290 G4cout << "Ptarget after boost to Lab: << 280 (Ptarget-Qmomentum).mag2() <= Mtarget2 ); 291 #endif << 281 */ // Uzhi * >> 282 >> 283 >> 284 } while (( (Pprojectile+Qmomentum).mag2() < Mprojectile2 || // Uzhi No without excitation >> 285 (Ptarget -Qmomentum).mag2() < Mtarget2 ) || // Uzhi >> 286 ( (Pprojectile+Qmomentum).mag2() < ProjectileDiffCut2 && // Uzhi No double Diffraction >> 287 (Ptarget -Qmomentum).mag2() < NuclearNucleonDiffCut2) );// Uzhi >> 288 >> 289 if((Ptarget-Qmomentum).mag2() < NuclearNucleonDiffCut2) // Uzhi Projectile diffraction >> 290 { >> 291 G4double TMinusNew=SqrtS-PMinusNew; >> 292 Qminus=Ptarget.minus()-TMinusNew; >> 293 TPlusNew=TargMassT2/TMinusNew; >> 294 Qplus=Ptarget.plus()-TPlusNew; >> 295 >> 296 Qmomentum.setPz( (Qplus-Qminus)/2 ); >> 297 Qmomentum.setE( (Qplus+Qminus)/2 ); >> 298 } >> 299 else if((Pprojectile+Qmomentum).mag2() < ProjectileDiffCut2) // Uzhi Target diffraction >> 300 { >> 301 G4double PPlusNew=SqrtS-TPlusNew; >> 302 Qplus=PPlusNew-Pprojectile.plus(); >> 303 PMinusNew=ProjMassT2/PPlusNew; >> 304 Qminus=PMinusNew-Pprojectile.minus(); >> 305 >> 306 Qmomentum.setPz( (Qplus-Qminus)/2 ); >> 307 Qmomentum.setE( (Qplus+Qminus)/2 ); >> 308 }; >> 309 >> 310 Pprojectile += Qmomentum; >> 311 Ptarget -= Qmomentum; >> 312 >> 313 // Vova >> 314 >> 315 /* >> 316 Pprojectile.setPz(0.); >> 317 Pprojectile.setE(SqrtS-M0target); >> 318 >> 319 Ptarget.setPz(0.); >> 320 Ptarget.setE(M0target); >> 321 */ >> 322 >> 323 //G4cout << "Pprojectile with Q : " << Pprojectile << G4endl; >> 324 //G4cout << "Ptarget with Q : " << Ptarget << G4endl; >> 325 >> 326 // G4cout << "Projectile back: " << toLab * Pprojectile << G4endl; >> 327 // G4cout << "Target back: " << toLab * Ptarget << G4endl; >> 328 >> 329 // Transform back and update SplitableHadron Participant. >> 330 Pprojectile.transform(toLab); >> 331 Ptarget.transform(toLab); >> 332 >> 333 //G4cout << "Pprojectile with Q M: " << Pprojectile<<" "<< Pprojectile.mag() << G4endl; >> 334 //G4cout << "Ptarget with Q M: " << Ptarget <<" "<< Ptarget.mag() << G4endl; >> 335 >> 336 //G4cout << "Target mass " << Ptarget.mag() << G4endl; >> 337 >> 338 target->Set4Momentum(Ptarget); >> 339 >> 340 //G4cout << "Projectile mass " << Pprojectile.mag() << G4endl; 292 341 293 target->Set4Momentum(Ptarget); << 342 projectile->Set4Momentum(Pprojectile); 294 projectile->Set4Momentum(Pprojectile); << 295 343 296 return true; << 344 return true; 297 } 345 } 298 346 299 347 300 G4ExcitedString * G4QGSDiffractiveExcitation:: 348 G4ExcitedString * G4QGSDiffractiveExcitation:: 301 String(G4VSplitableHadron * hadron, G4bool isP << 349 String(G4VSplitableHadron * hadron, G4bool isProjectile) const 302 { 350 { 303 hadron->SplitUp(); << 351 hadron->SplitUp(); 304 G4Parton *start= hadron->GetNextParton(); << 352 G4Parton *start= hadron->GetNextParton(); 305 if ( start==NULL) << 353 if ( start==NULL) 306 { G4cout << " G4QGSDiffractiveExcitation::St << 354 { G4cout << " G4FTFModel::String() Error:No start parton found"<< G4endl; 307 return NULL; << 355 return NULL; 308 } << 356 } 309 G4Parton *end = hadron->GetNextParton(); << 357 G4Parton *end = hadron->GetNextParton(); 310 if ( end==NULL) << 358 if ( end==NULL) 311 { G4cout << " G4QGSDiffractiveExcitation::St << 359 { G4cout << " G4FTFModel::String() Error:No end parton found"<< G4endl; 312 return NULL; << 360 return NULL; 313 } << 361 } 314 << 362 315 G4ExcitedString * string; << 363 G4ExcitedString * string; 316 if ( isProjectile ) << 364 if ( isProjectile ) 317 { << 365 { 318 string= new G4ExcitedString(end,start, +1) << 366 string= new G4ExcitedString(end,start, +1); 319 } else { << 367 } else { 320 string= new G4ExcitedString(start,end, -1) << 368 string= new G4ExcitedString(start,end, -1); 321 } << 369 } 322 << 370 323 string->SetPosition(hadron->GetPosition()); << 371 string->SetPosition(hadron->GetPosition()); 324 << 372 325 // momenta of string ends << 373 // momenta of string ends 326 << 374 G4double ptSquared= hadron->Get4Momentum().perp2(); 327 G4double maxAvailMomentumSquared=sqr(hadron- << 375 G4double transverseMassSquared= hadron->Get4Momentum().plus() 328 << 376 * hadron->Get4Momentum().minus(); 329 G4double widthOfPtSquare = 0.5*sqr(GeV); << 377 330 G4ThreeVector pt=GaussianPt(widthOfPtSquare, << 378 331 << 379 G4double maxAvailMomentumSquared= 332 G4LorentzVector Pstart(G4LorentzVector(pt,0. << 380 sqr( std::sqrt(transverseMassSquared) - std::sqrt(ptSquared) ); 333 G4LorentzVector Pend; << 381 334 Pend.setPx(hadron->Get4Momentum().px() - pt. << 382 G4double widthOfPtSquare = 0.25; // Uzhi <Pt^2>=0.25 ?????????????????? 335 Pend.setPy(hadron->Get4Momentum().py() - pt. << 383 G4ThreeVector pt=GaussianPt(widthOfPtSquare,maxAvailMomentumSquared); 336 << 384 337 G4double tm1=hadron->Get4Momentum().minus() << 385 G4LorentzVector Pstart(G4LorentzVector(pt,0.)); 338 ( Pend.perp2()-Pstart.perp2() ) / had << 386 G4LorentzVector Pend; 339 << 387 Pend.setPx(hadron->Get4Momentum().px() - pt.x()); 340 G4double tm2= std::sqrt( std::max(0., sqr(tm << 388 Pend.setPy(hadron->Get4Momentum().py() - pt.y()); 341 4. * Pend.perp2() * hadron->Get4Momen << 389 342 / hadron->Get4Momentum().plus() )); << 390 G4double tm1=hadron->Get4Momentum().minus() + 343 << 391 ( Pend.perp2()-Pstart.perp2() ) / hadron->Get4Momentum().plus(); 344 G4int Sign= isProjectile ? -1 : 1; << 392 345 << 393 G4double tm2= std::sqrt( std::max(0., sqr(tm1) - 346 G4double endMinus = 0.5 * (tm1 + Sign*tm2); << 394 4. * Pend.perp2() * hadron->Get4Momentum().minus() 347 G4double startMinus= hadron->Get4Momentum(). << 395 / hadron->Get4Momentum().plus() )); 348 << 396 349 G4double startPlus= Pstart.perp2() / startM << 397 G4int Sign= isProjectile ? -1 : 1; 350 G4double endPlus = hadron->Get4Momentum().p << 398 351 << 399 G4double endMinus = 0.5 * (tm1 + Sign*tm2); 352 Pstart.setPz(0.5*(startPlus - startMinus)); << 400 G4double startMinus= hadron->Get4Momentum().minus() - endMinus; 353 Pstart.setE(0.5*(startPlus + startMinus)); << 401 354 << 402 G4double startPlus= Pstart.perp2() / startMinus; 355 Pend.setPz(0.5*(endPlus - endMinus)); << 403 G4double endPlus = hadron->Get4Momentum().plus() - startPlus; 356 Pend.setE(0.5*(endPlus + endMinus)); << 404 357 << 405 Pstart.setPz(0.5*(startPlus - startMinus)); 358 start->Set4Momentum(Pstart); << 406 Pstart.setE(0.5*(startPlus + startMinus)); 359 end->Set4Momentum(Pend); << 407 360 << 408 Pend.setPz(0.5*(endPlus - endMinus)); 361 #ifdef debugQGSdiffExictation << 409 Pend.setE(0.5*(endPlus + endMinus)); 362 G4cout << " generated string flavors << 410 363 G4cout << " generated string momenta: qu << 411 start->Set4Momentum(Pstart); 364 G4cout << " generated string momenta: Diqu << 412 end->Set4Momentum(Pend); 365 G4cout << " sum of ends << 413 366 G4cout << " Original << 414 #ifdef G4_FTFDEBUG 367 #endif << 415 G4cout << " generated string flavors " << start->GetPDGcode() << " / " << end->GetPDGcode() << G4endl; 368 << 416 G4cout << " generated string momenta: quark " << start->Get4Momentum() << "mass : " <<start->Get4Momentum().mag()<< G4endl; 369 return string; << 417 G4cout << " generated string momenta: Diquark " << end ->Get4Momentum() << "mass : " <<end->Get4Momentum().mag()<< G4endl; >> 418 G4cout << " sum of ends " << Pstart+Pend << G4endl; >> 419 G4cout << " Original " << hadron->Get4Momentum() << G4endl; >> 420 #endif >> 421 >> 422 return string; 370 } 423 } 371 424 372 425 373 // --------- private methods ----------------- 426 // --------- private methods ---------------------- 374 427 375 G4double G4QGSDiffractiveExcitation::ChooseP(G << 428 G4double G4QGSDiffractiveExcitation::ChooseP(G4double Pmin, G4double Pmax) const // Uzhi 376 { 429 { 377 // choose an x between Xmin and Xmax with P( << 430 // choose an x between Xmin and Xmax with P(x) ~ 1/x 378 // to be improved... << 379 431 380 G4double range=Pmax-Pmin; << 432 // to be improved... 381 433 382 if ( Pmin <= 0. || range <=0. ) << 434 G4double range=Pmax-Pmin; // Uzhi 383 { << 435 384 G4cout << " Pmin, range : " << Pmin << " , << 436 if ( Pmin <= 0. || range <=0. ) 385 throw G4HadronicException(__FILE__, __LINE << 437 { 386 } << 438 G4cout << " Pmin, range : " << Pmin << " , " << range << G4endl; 387 << 439 throw G4HadronicException(__FILE__, __LINE__, "G4QGSDiffractiveExcitation::ChooseP : Invalid arguments "); 388 G4double P; << 440 } 389 P=Pmin * G4Pow::GetInstance()->powA(Pmax/Pmi << 441 390 //debug-hpw cout << "DiffractiveX "<<x<<G4en << 442 G4double P; 391 return P; << 443 /* // Uzhi 392 } << 444 do { >> 445 x=Xmin + G4UniformRand() * range; >> 446 } while ( Xmin/x < G4UniformRand() ); >> 447 */ // Uzhi 393 448 >> 449 P=Pmin * std::pow(Pmax/Pmin,G4UniformRand()); // Uzhi >> 450 >> 451 //debug-hpw cout << "DiffractiveX "<<x<<G4endl; >> 452 return P; >> 453 } 394 454 395 G4ThreeVector G4QGSDiffractiveExcitation::Gaus << 455 G4ThreeVector G4QGSDiffractiveExcitation::GaussianPt(G4double AveragePt2, G4double maxPtSquare) const // Uzhi 396 { // @@ this method is used in FTF 456 { // @@ this method is used in FTFModel as well. Should go somewhere common! 397 G4double Pt2; << 457 >> 458 G4double Pt2; >> 459 /* // Uzhi >> 460 do { >> 461 pt2=widthSquare * std::log( G4UniformRand() ); >> 462 } while ( pt2 > maxPtSquare); >> 463 */ // Uzhi >> 464 >> 465 Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() * (std::exp(-maxPtSquare/AveragePt2)-1.));// Uzhi >> 466 >> 467 G4double Pt=std::sqrt(Pt2); >> 468 >> 469 G4double phi=G4UniformRand() * twopi; >> 470 >> 471 return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.); >> 472 } >> 473 >> 474 G4QGSDiffractiveExcitation::G4QGSDiffractiveExcitation(const G4QGSDiffractiveExcitation &) >> 475 { >> 476 throw G4HadronicException(__FILE__, __LINE__, "G4QGSDiffractiveExcitation copy contructor not meant to be called"); >> 477 } 398 478 399 Pt2 = -AveragePt2 * G4Log(1. + G4UniformRand << 400 479 401 G4double Pt=std::sqrt(Pt2); << 480 G4QGSDiffractiveExcitation::~G4QGSDiffractiveExcitation() >> 481 { >> 482 } 402 483 403 G4double phi=G4UniformRand() * twopi; << 404 484 405 return G4ThreeVector (Pt*std::cos(phi), Pt*s << 485 const G4QGSDiffractiveExcitation & G4QGSDiffractiveExcitation::operator=(const G4QGSDiffractiveExcitation &) >> 486 { >> 487 throw G4HadronicException(__FILE__, __LINE__, "G4QGSDiffractiveExcitation = operator meant to be called"); >> 488 return *this; >> 489 } >> 490 >> 491 >> 492 int G4QGSDiffractiveExcitation::operator==(const G4QGSDiffractiveExcitation &) const >> 493 { >> 494 throw G4HadronicException(__FILE__, __LINE__, "G4QGSDiffractiveExcitation == operator meant to be called"); >> 495 return false; >> 496 } >> 497 >> 498 int G4QGSDiffractiveExcitation::operator!=(const G4QGSDiffractiveExcitation &) const >> 499 { >> 500 throw G4HadronicException(__FILE__, __LINE__, "G4QGSDiffractiveExcitation != operator meant to be called"); >> 501 return true; 406 } 502 } 407 503