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: G4FTFAnnihilation.cc 97625 2016-06-06 13:35:58Z gcosmo $ 27 // 28 // 28 29 29 // ------------------------------------------- 30 // ------------------------------------------------------------ 30 // GEANT 4 class implemetation file 31 // GEANT 4 class implemetation file 31 // 32 // 32 // ---------------- G4FTFAnnihilation --- 33 // ---------------- G4FTFAnnihilation -------------- 33 // by V. Uzhinsky, Spring 2011. 34 // by V. Uzhinsky, Spring 2011. 34 // Take a projectile and a targ 35 // Take a projectile and a target 35 // make annihilation or re-orangement o 36 // make annihilation or re-orangement of quarks and anti-quarks. 36 // Ideas of Quark-Gluon-String model my A. 37 // Ideas of Quark-Gluon-String model my A. Capella and A.B. Kaidalov 37 // are implemented. 38 // are implemented. 38 // ------------------------------------------- 39 // --------------------------------------------------------------------- 39 40 40 #include "globals.hh" 41 #include "globals.hh" 41 #include "Randomize.hh" 42 #include "Randomize.hh" 42 #include "G4PhysicalConstants.hh" 43 #include "G4PhysicalConstants.hh" 43 #include "G4SystemOfUnits.hh" 44 #include "G4SystemOfUnits.hh" 44 45 45 #include "G4DiffractiveSplitableHadron.hh" 46 #include "G4DiffractiveSplitableHadron.hh" 46 #include "G4DiffractiveExcitation.hh" 47 #include "G4DiffractiveExcitation.hh" 47 #include "G4FTFParameters.hh" 48 #include "G4FTFParameters.hh" 48 #include "G4ElasticHNScattering.hh" 49 #include "G4ElasticHNScattering.hh" 49 #include "G4FTFAnnihilation.hh" 50 #include "G4FTFAnnihilation.hh" 50 51 51 #include "G4LorentzRotation.hh" 52 #include "G4LorentzRotation.hh" 52 #include "G4RotationMatrix.hh" 53 #include "G4RotationMatrix.hh" 53 #include "G4ThreeVector.hh" 54 #include "G4ThreeVector.hh" 54 #include "G4ParticleDefinition.hh" 55 #include "G4ParticleDefinition.hh" 55 #include "G4VSplitableHadron.hh" 56 #include "G4VSplitableHadron.hh" 56 #include "G4ExcitedString.hh" 57 #include "G4ExcitedString.hh" 57 #include "G4ParticleTable.hh" 58 #include "G4ParticleTable.hh" 58 #include "G4Neutron.hh" 59 #include "G4Neutron.hh" 59 #include "G4ParticleDefinition.hh" 60 #include "G4ParticleDefinition.hh" 60 61 61 #include "G4Exp.hh" 62 #include "G4Exp.hh" 62 #include "G4Log.hh" 63 #include "G4Log.hh" 63 #include "G4Pow.hh" 64 #include "G4Pow.hh" 64 65 >> 66 //#include "G4ios.hh" 65 //#include "UZHI_diffraction.hh" 67 //#include "UZHI_diffraction.hh" 66 68 67 #include "G4ParticleTable.hh" << 69 #include "G4ParticleTable.hh" // Uzhi March 2016 68 << 69 //============================================ 70 //============================================================================ 70 71 71 //#define debugFTFannih 72 //#define debugFTFannih 72 73 73 74 74 //============================================ 75 //============================================================================ 75 76 76 G4FTFAnnihilation::G4FTFAnnihilation() {} 77 G4FTFAnnihilation::G4FTFAnnihilation() {} 77 78 78 79 79 //============================================ 80 //============================================================================ 80 81 81 G4FTFAnnihilation::~G4FTFAnnihilation() {} 82 G4FTFAnnihilation::~G4FTFAnnihilation() {} 82 83 83 84 84 //============================================ 85 //============================================================================ 85 86 86 G4bool G4FTFAnnihilation::Annihilate( G4VSplit 87 G4bool G4FTFAnnihilation::Annihilate( G4VSplitableHadron* projectile, 87 G4VSplit 88 G4VSplitableHadron* target, 88 G4VSplit 89 G4VSplitableHadron*& AdditionalString, 89 G4FTFPar 90 G4FTFParameters* theParameters ) const { 90 91 >> 92 //theParameters->SetProbabilityOfAnnihilation( 0.0 ); // Uzhi March 2016 ??? for other Anti_bar annih. >> 93 91 #ifdef debugFTFannih 94 #ifdef debugFTFannih 92 G4cout << "---------------------------- Anni 95 G4cout << "---------------------------- Annihilation----------------" << G4endl; 93 #endif 96 #endif 94 97 95 CommonVariables common; << 96 << 97 // Projectile parameters 98 // Projectile parameters 98 common.Pprojectile = projectile->Get4Momentu << 99 G4LorentzVector Pprojectile = projectile->Get4Momentum(); 99 G4int ProjectilePDGcode = projectile->GetDef 100 G4int ProjectilePDGcode = projectile->GetDefinition()->GetPDGEncoding(); 100 if ( ProjectilePDGcode > 0 ) { 101 if ( ProjectilePDGcode > 0 ) { 101 target->SetStatus( 3 ); << 102 target->SetStatus( 3 ); // 2->3 Uzhi Oct 2014 102 return false; 103 return false; 103 } 104 } 104 G4double M0projectile2 = common.Pprojectile. << 105 //G4double M0projectile = Pprojectile.mag(); >> 106 //G4double M0projectile2 = projectile->GetDefinition()->GetPDGMass() * >> 107 // projectile->GetDefinition()->GetPDGMass(); >> 108 G4double M0projectile2 = Pprojectile.mag2(); 105 109 106 // Target parameters 110 // Target parameters 107 G4int TargetPDGcode = target->GetDefinition( 111 G4int TargetPDGcode = target->GetDefinition()->GetPDGEncoding(); 108 common.Ptarget = target->Get4Momentum(); << 112 G4LorentzVector Ptarget = target->Get4Momentum(); 109 G4double M0target2 = common.Ptarget.mag2(); << 113 //G4double M0target = Ptarget.mag(); >> 114 //G4double M0target2 = target->GetDefinition()->GetPDGMass() * >> 115 // target->GetDefinition()->GetPDGMass(); >> 116 G4double M0target2 = Ptarget.mag2(); 110 117 111 #ifdef debugFTFannih 118 #ifdef debugFTFannih 112 G4cout << "PDG codes " << ProjectilePDGcode 119 G4cout << "PDG codes " << ProjectilePDGcode << " " << TargetPDGcode << G4endl 113 << "Pprojec " << common.Pprojectile < << 120 << "Pprojec " << Pprojectile << " " << Pprojectile.mag() << G4endl 114 << "Ptarget " << common.Ptarget << << 121 << "Ptarget " << Ptarget << " " << Ptarget.mag() << G4endl 115 << "M0 proj target " << std::sqrt( M0 122 << "M0 proj target " << std::sqrt( M0projectile2 ) 116 << " " << std::sqrt( M0target2 ) << G 123 << " " << std::sqrt( M0target2 ) << G4endl; 117 #endif 124 #endif 118 125 >> 126 G4double AveragePt2 = theParameters->GetAveragePt2(); >> 127 119 // Kinematical properties of the interaction 128 // Kinematical properties of the interactions 120 G4LorentzVector Psum = common.Pprojectile + << 129 G4LorentzVector Psum; // 4-momentum in CMS 121 common.S = Psum.mag2(); << 130 Psum = Pprojectile + Ptarget; 122 common.SqrtS = std::sqrt( common.S ); << 131 G4double S = Psum.mag2(); >> 132 123 #ifdef debugFTFannih 133 #ifdef debugFTFannih 124 G4cout << "Psum SqrtS S " << Psum << " " << << 134 G4cout << "Psum SqrtS S " << Psum << " " << std::sqrt( S ) << " " << S << G4endl; 125 #endif 135 #endif 126 136 127 // Transform momenta to cms and then rotate 137 // Transform momenta to cms and then rotate parallel to z axis 128 G4LorentzRotation toCms( -1*Psum.boostVector 138 G4LorentzRotation toCms( -1*Psum.boostVector() ); 129 G4LorentzVector Ptmp( toCms*common.Pprojecti << 139 G4LorentzVector Ptmp = toCms*Pprojectile; 130 toCms.rotateZ( -1*Ptmp.phi() ); 140 toCms.rotateZ( -1*Ptmp.phi() ); 131 toCms.rotateY( -1*Ptmp.theta() ); 141 toCms.rotateY( -1*Ptmp.theta() ); 132 common.toLab = toCms.inverse(); << 142 G4LorentzRotation toLab( toCms.inverse() ); 133 143 134 if ( G4UniformRand() <= G4Pow::GetInstance() << 144 G4double SqrtS = std::sqrt( S ); 135 common.RotateStrings = true; << 145 G4double maxPtSquare; 136 common.RandomRotation.rotateZ( 2.0*pi*G4Un << 146 G4double X_a( 0.0 ), X_b( 0.0 ), X_c( 0.0 ), X_d( 0.0 ); 137 common.RandomRotation.rotateY( std::acos( << 138 common.RandomRotation.rotateZ( 2.0*pi*G4Un << 139 } << 140 147 141 G4double MesonProdThreshold = projectile->Ge 148 G4double MesonProdThreshold = projectile->GetDefinition()->GetPDGMass() + 142 target->GetDef 149 target->GetDefinition()->GetPDGMass() + 143 ( 2.0*140.0 + 150 ( 2.0*140.0 + 16.0 )*MeV; // 2 Mpi + DeltaE 144 G4double Prel2 = sqr(common.S) + sqr(M0proje << 151 G4double Prel2 = S*S + M0projectile2*M0projectile2 + M0target2*M0target2 - 145 - 2.0*( common.S*(M0project << 152 2.0*S*M0projectile2 - 2.0*S*M0target2 - 2.0*M0projectile2*M0target2; 146 Prel2 /= common.S; << 153 Prel2 /= S; 147 G4double X_a = 0.0, X_b = 0.0, X_c = 0.0, X_ << 154 //G4cout << "Prel2 " << Prel2 << G4endl; 148 if ( Prel2 <= 0.0 ) { << 155 if ( Prel2 <= 0.0 ) { // *MeV*MeV 1600. 149 // Annihilation at rest! Values are copied 156 // Annihilation at rest! Values are copied from Parameters 150 X_a = 625.1; // mb // 3-shirt diagram << 157 X_a = 625.1; // mb // 3-shirt diagram 151 X_b = 0.0; // mb // anti-quark-quark << 158 X_b = 0.0; // 9.780 12 Dec. 2012; // mb // anti-quark-quark annihilation 152 X_c = 49.989; // mb // 2 Q-Qbar string << 159 X_c = 49.989; // mb 153 X_d = 6.614; // mb // One Q-Qbar strin << 160 X_d = 6.614; // mb >> 161 154 #ifdef debugFTFannih 162 #ifdef debugFTFannih 155 G4cout << "Annih at Rest X a b c d " << X_ 163 G4cout << "Annih at Rest X a b c d " << X_a << " " << X_b << " " << X_c << " " << X_d 156 << G4endl; 164 << G4endl; 157 #endif 165 #endif >> 166 158 } else { // Annihilation in flight! 167 } else { // Annihilation in flight! 159 G4double FlowF = 1.0 / std::sqrt( Prel2 )* 168 G4double FlowF = 1.0 / std::sqrt( Prel2 )*GeV; >> 169 160 // Process cross sections 170 // Process cross sections 161 X_a = 25.0*FlowF; // mb 3-shirt diagram 171 X_a = 25.0*FlowF; // mb 3-shirt diagram 162 if ( common.SqrtS < MesonProdThreshold ) { << 172 if ( SqrtS < MesonProdThreshold ) { 163 X_b = 3.13 + 140.0*G4Pow::GetInstance()- << 173 X_b = 3.13 + 140.0*G4Pow::GetInstance()->powA( ( MesonProdThreshold - SqrtS )/GeV, 2.5 ); 164 } else { 174 } else { 165 X_b = 6.8*GeV / common.SqrtS; // mb ant << 175 X_b = 6.8*GeV / SqrtS; // mb anti-quark-quark annihilation 166 } 176 } 167 if ( projectile->GetDefinition()->GetPDGMa 177 if ( projectile->GetDefinition()->GetPDGMass() + target->GetDefinition()->GetPDGMass() 168 > common.SqrtS ) { << 178 > SqrtS ) { 169 X_b = 0.0; 179 X_b = 0.0; 170 } 180 } 171 // This can be in an interaction of low en 181 // This can be in an interaction of low energy anti-baryon with off-shell nuclear nucleon 172 X_c = 2.0 * FlowF * sqr( projectile->GetDe 182 X_c = 2.0 * FlowF * sqr( projectile->GetDefinition()->GetPDGMass() + 173 target->GetDefini << 183 target->GetDefinition()->GetPDGMass() ) / S; // mb re-arrangement of 174 << 184 // 2 quarks and 2 anti-quarks 175 X_d = 23.3*GeV*GeV / common.S; // mb anti- << 185 X_d = 23.3*GeV*GeV / S; // mb anti-quark-quark string creation >> 186 176 #ifdef debugFTFannih 187 #ifdef debugFTFannih 177 G4cout << "Annih in Flight X a b c d " << 188 G4cout << "Annih in Flight X a b c d " << X_a << " " << X_b << " " << X_c << " " << X_d 178 << G4endl << "SqrtS MesonProdThresh << 189 << G4endl << "SqrtS MesonProdThreshold " << SqrtS << " " << MesonProdThreshold 179 << G4endl; 190 << G4endl; 180 #endif 191 #endif >> 192 181 } 193 } 182 194 183 G4bool isUnknown = false; << 195 if ((ProjectilePDGcode == -2212 || ProjectilePDGcode == -2214)&& ( TargetPDGcode == 2212 || TargetPDGcode == 2214 ) ) { 184 if ( TargetPDGcode == 2212 || TargetPDGcod << 196 X_b *= 5.0; X_c *= 5.0; X_d *= 6.0; // Pbar P 185 if ( ProjectilePDGcode == -2212 || << 197 } else if ((ProjectilePDGcode == -2212 || ProjectilePDGcode == -2214)&& ( TargetPDGcode == 2112 || TargetPDGcode == 2114 ) ) { 186 X_b *= 5.0; X_c *= 5.0; X_d *= 6.0; // << 198 X_b *= 4.0; X_c *= 4.0; X_d *= 4.0; // Pbar N 187 } else if ( ProjectilePDGcode == -2112 || << 199 } else if ((ProjectilePDGcode == -2112 || ProjectilePDGcode == -2114)&& ( TargetPDGcode == 2212 || TargetPDGcode == 2214 ) ) { 188 X_b *= 4.0; X_c *= 4.0; X_d *= 4.0; // << 200 X_b *= 4.0; X_c *= 4.0; X_d *= 4.0; // NeutrBar P 189 } else if ( ProjectilePDGcode == -3122 ) { << 201 } else if ((ProjectilePDGcode == -2112 || ProjectilePDGcode == -2114)&& ( TargetPDGcode == 2112 || TargetPDGcode == 2114 ) ) { 190 X_b *= 3.0; X_c *= 3.0; X_d *= 2.0; // << 202 X_b *= 5.0; X_c *= 5.0; X_d *= 6.0; // NeutrBar N 191 } else if ( ProjectilePDGcode == -3112 ) { << 203 } else if ((ProjectilePDGcode == -3122 || ProjectilePDGcode == -3124)&& ( TargetPDGcode == 2212 || TargetPDGcode == 2214 ) ) { 192 X_b *= 2.0; X_c *= 2.0; X_d *= 0.0; // << 204 X_b *= 3.0; X_c *= 3.0; X_d *= 2.0; // LambdaBar P 193 } else if ( ProjectilePDGcode == -3212 ) { << 205 } else if ((ProjectilePDGcode == -3122 || ProjectilePDGcode == -3124)&& ( TargetPDGcode == 2112 || TargetPDGcode == 2114 ) ) { 194 X_b *= 3.0; X_c *= 3.0; X_d *= 2.0; // << 206 X_b *= 3.0; X_c *= 3.0; X_d *= 2.0; // LambdaBar N 195 } else if ( ProjectilePDGcode == -3222 ) { << 207 } else if ((ProjectilePDGcode == -3112 || ProjectilePDGcode == -3114)&& ( TargetPDGcode == 2212 || TargetPDGcode == 2214 ) ) { 196 X_b *= 4.0; X_c *= 4.0; X_d *= 2.0; // << 208 X_b *= 2.0; X_c *= 2.0; X_d *= 0.0; // Sigma-Bar P 197 } else if ( ProjectilePDGcode == -3312 ) { << 209 } else if ((ProjectilePDGcode == -3112 || ProjectilePDGcode == -3114)&& ( TargetPDGcode == 2112 || TargetPDGcode == 2114 ) ) { 198 X_b *= 1.0; X_c *= 1.0; X_d *= 0.0; // << 210 X_b *= 4.0; X_c *= 4.0; X_d *= 2.0; // Sigma-Bar N 199 } else if ( ProjectilePDGcode == -3322 ) { << 211 } else if ((ProjectilePDGcode == -3212 || ProjectilePDGcode == -3214)&& ( TargetPDGcode == 2212 || TargetPDGcode == 2214 ) ) { 200 X_b *= 2.0; X_c *= 2.0; X_d *= 0.0; // << 212 X_b *= 3.0; X_c *= 3.0; X_d *= 2.0; // Sigma0Bar P 201 } else if ( ProjectilePDGcode == -3334 ) { << 213 } else if ((ProjectilePDGcode == -3212 || ProjectilePDGcode == -3214)&& ( TargetPDGcode == 2112 || TargetPDGcode == 2114 ) ) { 202 X_b *= 0.0; X_c *= 0.0; X_d *= 0.0; // << 214 X_b *= 3.0; X_c *= 3.0; X_d *= 2.0; // Sigma0Bar N 203 } else { << 215 } else if ((ProjectilePDGcode == -3222 || ProjectilePDGcode == -3224)&& ( TargetPDGcode == 2212 || TargetPDGcode == 2214 ) ) { 204 isUnknown = true; << 216 X_b *= 4.0; X_c *= 4.0; X_d *= 2.0; // Sigma+Bar P 205 } << 217 } else if ((ProjectilePDGcode == -3222 || ProjectilePDGcode == -3224)&& ( TargetPDGcode == 2112 || TargetPDGcode == 2114 ) ) { 206 } else if ( TargetPDGcode == 2112 || Targe << 218 X_b *= 2.0; X_c *= 2.0; X_d *= 0.0; // Sigma+Bar N 207 if ( ProjectilePDGcode == -2212 || << 219 } else if ((ProjectilePDGcode == -3312 || ProjectilePDGcode == -3314)&& ( TargetPDGcode == 2212 || TargetPDGcode == 2214 ) ) { 208 X_b *= 4.0; X_c *= 4.0; X_d *= 4.0; // << 220 X_b *= 1.0; X_c *= 1.0; X_d *= 0.0; // Xi-Bar P 209 } else if ( ProjectilePDGcode == -2112 || << 221 } else if ((ProjectilePDGcode == -3312 || ProjectilePDGcode == -3314)&& ( TargetPDGcode == 2112 || TargetPDGcode == 2114 ) ) { 210 X_b *= 5.0; X_c *= 5.0; X_d *= 6.0; // << 222 X_b *= 2.0; X_c *= 2.0; X_d *= 0.0; // Xi-Bar N 211 } else if ( ProjectilePDGcode == -3122 ) { << 223 } else if ((ProjectilePDGcode == -3322 || ProjectilePDGcode == -3324)&& ( TargetPDGcode == 2212 || TargetPDGcode == 2214 ) ) { 212 X_b *= 3.0; X_c *= 3.0; X_d *= 2.0; // << 224 X_b *= 2.0; X_c *= 2.0; X_d *= 0.0; // Xi0Bar P 213 } else if ( ProjectilePDGcode == -3112 ) { << 225 } else if ((ProjectilePDGcode == -3322 || ProjectilePDGcode == -3324)&& ( TargetPDGcode == 2112 || TargetPDGcode == 2114 ) ) { 214 X_b *= 4.0; X_c *= 4.0; X_d *= 2.0; // << 226 X_b *= 1.0; X_c *= 1.0; X_d *= 0.0; // Xi0Bar N 215 } else if ( ProjectilePDGcode == -3212 ) { << 227 } else if ( ProjectilePDGcode == -3334 && ( TargetPDGcode == 2212 || TargetPDGcode == 2214 ) ) { 216 X_b *= 3.0; X_c *= 3.0; X_d *= 2.0; // << 228 X_b *= 0.0; X_c *= 0.0; X_d *= 0.0; // Omega-Bar P 217 } else if ( ProjectilePDGcode == -3222 ) { << 229 } else if ( ProjectilePDGcode == -3334 && ( TargetPDGcode == 2112 || TargetPDGcode == 2114 ) ) { 218 X_b *= 2.0; X_c *= 2.0; X_d *= 0.0; // << 230 X_b *= 0.0; X_c *= 0.0; X_d *= 0.0; // Omega-Bar N 219 } else if ( ProjectilePDGcode == -3312 ) { << 220 X_b *= 2.0; X_c *= 2.0; X_d *= 0.0; // << 221 } else if ( ProjectilePDGcode == -3322 ) { << 222 X_b *= 1.0; X_c *= 1.0; X_d *= 0.0; // << 223 } else if ( ProjectilePDGcode == -3334 ) { << 224 X_b *= 0.0; X_c *= 0.0; X_d *= 0.0; // << 225 } else { << 226 isUnknown = true; << 227 } << 228 } else { 231 } else { 229 isUnknown = true; << 230 } << 231 if ( isUnknown ) { << 232 G4cout << "Unknown anti-baryon for FTF ann 232 G4cout << "Unknown anti-baryon for FTF annihilation: PDGcodes - " 233 << ProjectilePDGcode << " " << Targ 233 << ProjectilePDGcode << " " << TargetPDGcode << G4endl; 234 } 234 } 235 235 236 #ifdef debugFTFannih 236 #ifdef debugFTFannih 237 G4cout << "Annih Actual X a b c d " << X_a < 237 G4cout << "Annih Actual X a b c d " << X_a << " " << X_b << " " << X_c << " " << X_d << G4endl; 238 #endif 238 #endif 239 239 240 G4double Xannihilation = X_a + X_b + X_c + X 240 G4double Xannihilation = X_a + X_b + X_c + X_d; >> 241 //X_a=0.; // Uzhi >> 242 //X_b=0.; >> 243 //X_c=0.; >> 244 //X_d=0.; >> 245 //Xannihilation = X_a + X_b + X_c + X_d; 241 246 242 // Projectile unpacking 247 // Projectile unpacking 243 UnpackBaryon( ProjectilePDGcode, common.AQ[0 << 248 G4int AQ[3]; >> 249 UnpackBaryon( ProjectilePDGcode, AQ[0], AQ[1], AQ[2] ); 244 250 245 // Target unpacking 251 // Target unpacking 246 UnpackBaryon( TargetPDGcode, common.Q[0], co << 252 G4int Q[3]; >> 253 UnpackBaryon( TargetPDGcode, Q[0], Q[1], Q[2] ); 247 254 248 G4double Ksi = G4UniformRand(); 255 G4double Ksi = G4UniformRand(); 249 256 250 if ( Ksi < X_a / Xannihilation ) { 257 if ( Ksi < X_a / Xannihilation ) { 251 return Create3QuarkAntiQuarkStrings( proje << 252 } << 253 258 254 G4int resultCode = 99; << 259 // Simulation of 3 anti-quark-quark strings creation 255 if ( Ksi < (X_a + X_b) / Xannihilation ) { << 260 // Sampling of anti-quark order in projectile 256 resultCode = Create1DiquarkAntiDiquarkStri << 257 if ( resultCode == 0 ) { << 258 return true; << 259 } else if ( resultCode == 99 ) { << 260 return false; << 261 } << 262 } << 263 261 264 if ( Ksi < ( X_a + X_b + X_c ) / Xannihilati << 262 #ifdef debugFTFannih 265 resultCode = Create2QuarkAntiQuarkStrings( << 263 G4cout << "Process a, 3 shirt diagram" << G4endl; 266 if ( resultCode == 0 ) { << 264 #endif 267 return true; << 268 } else if ( resultCode == 99 ) { << 269 return false; << 270 } << 271 } << 272 265 273 if ( Ksi < ( X_a + X_b + X_c + X_d ) / Xanni << 266 G4int SampledCase = G4RandFlat::shootInt( G4long( 6 ) ); 274 return Create1QuarkAntiQuarkString( projec << 275 } << 276 267 277 return true; << 268 G4int Tmp1( 0 ), Tmp2( 0 ); 278 } << 269 if ( SampledCase == 0 ) { >> 270 } else if ( SampledCase == 1 ) { >> 271 Tmp1 = AQ[1]; AQ[1] = AQ[2]; AQ[2] = Tmp1; >> 272 } else if ( SampledCase == 2 ) { >> 273 Tmp1 = AQ[0]; AQ[0] = AQ[1]; AQ[1] = Tmp1; >> 274 } else if ( SampledCase == 3 ) { >> 275 Tmp1 = AQ[0]; Tmp2 = AQ[1]; AQ[0] = AQ[2]; AQ[1] = Tmp1; AQ[2] = Tmp2; >> 276 } else if ( SampledCase == 4 ) { >> 277 Tmp1 = AQ[0]; Tmp2 = AQ[1]; AQ[0] = Tmp2; AQ[1] = AQ[2]; AQ[2] = Tmp1; >> 278 } else if ( SampledCase == 5 ) { >> 279 Tmp1 = AQ[0]; Tmp2 = AQ[1]; AQ[0] = AQ[2]; AQ[1] = Tmp2; AQ[2] = Tmp1; >> 280 } 279 281 >> 282 // Set the string properties 280 283 281 //-------------------------------------------- << 284 //G4cout << "String 1 " << AQ[0] << " " << Q[0] << G4endl; >> 285 projectile->SplitUp(); >> 286 projectile->SetFirstParton( AQ[0] ); >> 287 projectile->SetSecondParton( Q[0] ); >> 288 projectile->SetStatus( 0 ); 282 289 283 G4bool G4FTFAnnihilation:: << 290 // Uzhi March 2016 start 284 Create3QuarkAntiQuarkStrings( G4VSplitableHadr << 291 G4int aAQ, aQ; 285 G4VSplitableHadr << 292 aAQ=std::abs( AQ[0] ); aQ=std::abs( Q[0] ); 286 G4VSplitableHadr << 293 G4int NewCode; 287 G4FTFParameters* << 294 G4double aKsi = G4UniformRand(); 288 G4FTFAnnihilatio << 295 289 // Simulation of 3 anti-quark - quark string << 296 if ( aAQ == aQ ) >> 297 { >> 298 if ( aAQ != 3 ) >> 299 { >> 300 NewCode = 111; // Pi0-meson >> 301 if ( aKsi < 0.5 ) >> 302 { >> 303 NewCode = 221; // Eta -meson >> 304 if ( aKsi < 0.25 ) {NewCode = 331;} // Eta'-meson >> 305 } >> 306 } else >> 307 { >> 308 NewCode = 221; // Eta -meson >> 309 if( aKsi < 0.5 ) {NewCode = 331;} // Eta'-meson >> 310 } >> 311 } else >> 312 { >> 313 if ( aAQ > aQ ){ NewCode = aAQ*100 + aQ*10 + 1; NewCode *= aAQ/AQ[0]; } >> 314 else { NewCode = aQ*100 + aAQ*10 + 1; NewCode *= aQ/Q[0]; } >> 315 } 290 316 291 #ifdef debugFTFannih << 317 G4ParticleDefinition* TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewCode ); 292 G4cout << "Process a, 3 shirt diagram" << G4 << 318 if(!TestParticle) return false; 293 #endif << 319 projectile->SetDefinition( TestParticle ); >> 320 >> 321 theParameters->SetProjMinDiffMass( 0.5 ); // Uzhi 2016 M+140 ?? >> 322 theParameters->SetProjMinNonDiffMass( 0.5 ); // Uzhi 2016 M+140 ?? >> 323 // Uzhi March 2016 end >> 324 >> 325 //G4cout << "String 2 " << Q[1] << " " << AQ[1] << G4endl; >> 326 target->SplitUp(); >> 327 target->SetFirstParton( Q[1] ); >> 328 target->SetSecondParton( AQ[1] ); >> 329 target->SetStatus( 0 ); >> 330 >> 331 // Uzhi March 2016 Start >> 332 aAQ=std::abs( AQ[1] ); aQ=std::abs( Q[1] ); aKsi = G4UniformRand(); >> 333 if ( aAQ == aQ ) >> 334 { >> 335 if ( aAQ != 3 ) >> 336 { >> 337 NewCode = 111; // Pi0-meson >> 338 if ( aKsi < 0.5 ) >> 339 { >> 340 NewCode = 221; // Eta -meson >> 341 if ( aKsi < 0.25 ) {NewCode = 331;} // Eta'-meson >> 342 } >> 343 } else >> 344 { >> 345 NewCode = 221; // Eta -meson >> 346 if( aKsi < 0.5 ) {NewCode = 331;} // Eta'-meson >> 347 } >> 348 } else >> 349 { >> 350 if ( aAQ > aQ ){ NewCode = aAQ*100 + aQ*10 + 1; NewCode *= aAQ/AQ[1]; } >> 351 else { NewCode = aQ*100 + aAQ*10 + 1; NewCode *= aQ/Q[1]; } >> 352 } >> 353 TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewCode ); >> 354 if(!TestParticle) return false; >> 355 target->SetDefinition( TestParticle ); >> 356 >> 357 theParameters->SetTarMinDiffMass( 0.5 ); // Uzhi 2016 M+140 ?? >> 358 theParameters->SetTarMinNonDiffMass( 0.5 ); // Uzhi 2016 M+140 ?? >> 359 // Uzhi March 2016 end >> 360 >> 361 >> 362 //G4cout << "String 3 " << AQ[2] << " " << Q[2] << G4endl; >> 363 AdditionalString = new G4DiffractiveSplitableHadron(); >> 364 >> 365 // Uzhi March 2016 start >> 366 aAQ=std::abs( AQ[2] ); aQ=std::abs( Q[2] ); aKsi = G4UniformRand(); >> 367 >> 368 if ( aAQ == aQ ) >> 369 { >> 370 if ( aAQ != 3 ) >> 371 { >> 372 NewCode = 111; // Pi0-meson >> 373 if ( aKsi < 0.5 ) >> 374 { >> 375 NewCode = 221; // Eta -meson >> 376 if ( aKsi < 0.25 ) {NewCode = 331;} // Eta'-meson >> 377 } >> 378 } else >> 379 { >> 380 NewCode = 221; // Eta -meson >> 381 if( aKsi < 0.5 ) {NewCode = 331;} // Eta'-meson >> 382 } >> 383 } else >> 384 { >> 385 if ( aAQ > aQ ){ NewCode = aAQ*100 + aQ*10 + 1; NewCode *= aAQ/AQ[2]; } >> 386 else { NewCode = aQ*100 + aAQ*10 + 1; NewCode *= aQ/Q[2]; } >> 387 } >> 388 >> 389 TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewCode ); >> 390 if(!TestParticle) return false; >> 391 AdditionalString->SetDefinition( TestParticle ); >> 392 // Uzhi March 2016 end >> 393 >> 394 AdditionalString->SplitUp(); >> 395 AdditionalString->SetFirstParton( AQ[2] ); >> 396 AdditionalString->SetSecondParton( Q[2] ); >> 397 AdditionalString->SetStatus( 0 ); >> 398 //G4cout << G4endl << "*AdditionalString in Annih" << AdditionalString << G4endl; >> 399 >> 400 // Sampling kinematical properties >> 401 // 1 string AQ[0]-Q[0]// 2 string AQ[1]-Q[1]// 3 string AQ[2]-Q[2] >> 402 >> 403 G4ThreeVector Quark_Mom[6]; >> 404 G4double ModMom2[6]; //ModMom[6] >> 405 >> 406 AveragePt2 = 200.0*200.0; maxPtSquare = S; >> 407 >> 408 G4double SumMt( 0.0 ); >> 409 G4double MassQ2 = 0.0; // 100.0*100.0*MeV*MeV; >> 410 G4int NumberOfTries( 0 ); >> 411 G4double ScaleFactor( 1.0 ); >> 412 >> 413 const G4int maxNumberOfLoops = 1000; >> 414 G4int loopCounter = 0; >> 415 do { >> 416 NumberOfTries++; >> 417 if ( NumberOfTries == 100*(NumberOfTries/100) ) { >> 418 // At large number of tries it would be better to reduce the values of <Pt^2> >> 419 ScaleFactor /= 2.0; >> 420 AveragePt2 *= ScaleFactor; >> 421 } >> 422 G4ThreeVector PtSum( 0.0, 0.0, 0.0 ); >> 423 for ( G4int i = 0; i < 6; i++ ) { >> 424 Quark_Mom [i] = GaussianPt( AveragePt2, maxPtSquare ); >> 425 PtSum += Quark_Mom[i]; >> 426 } >> 427 PtSum /= 6.0; >> 428 SumMt = 0.0; >> 429 for( G4int i = 0; i < 6; i++ ) { >> 430 Quark_Mom[i] -= PtSum; >> 431 //ModMom[i] = Quark_Mom[i].mag(); >> 432 ModMom2[i] = Quark_Mom[i].mag2(); >> 433 SumMt += std::sqrt( ModMom2[i] + MassQ2 ); >> 434 } >> 435 } while ( ( SumMt > SqrtS ) && >> 436 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */ >> 437 if ( loopCounter >= maxNumberOfLoops ) { >> 438 return false; >> 439 } 294 440 295 // Sampling kinematical properties of quark. << 441 G4double WminusTarget( 0.0 ), WplusProjectile( 0.0 ); 296 442 297 const G4int maxNumberOfLoops = 1000; << 443 // Closed is variant with sampling of Xs at minimum 298 G4double MassQ2 = 0.0; // Simp << 444 //G4double SumMod_anti = ModMom[0] + ModMom[1] + ModMom[2]; 299 // In p << 445 //Quark_Mom[0].setZ( ModMom[0]/SumMod_anti ); 300 G4double Quark_Xs[6]; << 446 //Quark_Mom[1].setZ( ModMom[1]/SumMod_anti ); 301 G4ThreeVector Quark_Mom[6]; << 447 //Quark_Mom[2].setZ( ModMom[2]/SumMod_anti ); 302 << 448 //G4double SumMod_bary = ModMom[3] + ModMom[4] + ModMom[5]; 303 G4double Alfa_R = 0.5; << 449 //Quark_Mom[3].setZ( ModMom[3]/SumMod_bary ); 304 G4double AveragePt2 = 200.0*200.0, maxPtSqua << 450 //Quark_Mom[4].setZ( ModMom[4]/SumMod_bary ); 305 G4double ScaleFactor = 1.0; << 451 //Quark_Mom[5].setZ( ModMom[5]/SumMod_bary ); 306 G4double Alfa = 0.0, Beta = 0.0; << 452 //G4double Alfa = SumMod_anti*SumMod_anti; 307 << 453 //G4double Beta = SumMod_bary*SumMod_bary; 308 G4int NumberOfTries = 0, loopCounter = 0; << 454 //G4double DecayMomentum2 = S*S + Alfa*Alfa + Beta*Beta 309 << 455 // - 2.0*S*Alfa - 2.0*S*Beta - 2.0*Alfa*Beta; 310 do { << 456 //WminusTarget = ( S - Alfa + Beta + std::sqrt( DecayMomentum2 ) )/2.0/SqrtS; 311 // Sampling X's of anti-baryon and baryon << 457 //WplusProjectile = SqrtS - Beta/WminusTarget; 312 G4double x1 = 0.0, x2 = 0.0, x3 = 0.0; << 458 // Closed is variant with sampling of Xs at minimum 313 G4double Product = 1.0; << 459 314 for ( G4int iCase = 0; iCase < 2; ++iCase << 460 // Sampling X's of anti-baryon >> 461 G4double Alfa_R = 0.5; >> 462 NumberOfTries = 0; >> 463 ScaleFactor = 1.0; >> 464 G4bool Succes( true ); >> 465 >> 466 loopCounter = 0; >> 467 do { >> 468 >> 469 Succes = true; >> 470 NumberOfTries++; >> 471 if ( NumberOfTries == 100*(NumberOfTries/100) ) { >> 472 // At large number of tries it would be better to reduce the values of Pt's >> 473 ScaleFactor /= 2.0; >> 474 } 315 475 316 G4double r1 = G4UniformRand(), r2 = G4Un << 317 if ( Alfa_R == 1.0 ) { 476 if ( Alfa_R == 1.0 ) { 318 x1 = 1.0 - std::sqrt( r1 ); << 477 G4double Xaq1 = 1.0 - std::sqrt( G4UniformRand() ); 319 x2 = (1.0 - x1) * r2; << 478 G4double Xaq2 = (1.0 - Xaq1) * G4UniformRand(); >> 479 G4double Xaq3 = 1.0 - Xaq1 - Xaq2; >> 480 Quark_Mom[0].setZ( Xaq1 ); Quark_Mom[1].setZ( Xaq2 ); Quark_Mom[2].setZ( Xaq3 ); 320 } else { 481 } else { 321 x1 = sqr( r1 ); << 482 G4double Xaq1 = sqr( G4UniformRand() ); 322 x2 = (1.0 - x1) * sqr( std::sin( pi/2. << 483 G4double Xaq2 = (1.0 - Xaq1)*sqr( std::sin( pi/2.0*G4UniformRand() ) ); >> 484 G4double Xaq3 = 1.0 - Xaq1 - Xaq2; >> 485 Quark_Mom[0].setZ( Xaq1 ); Quark_Mom[1].setZ( Xaq2 ); Quark_Mom[2].setZ( Xaq3 ); 323 } 486 } 324 x3 = 1.0 - x1 - x2; << 325 << 326 G4int index = iCase*3; // 0 for anti-ba << 327 Quark_Xs[index] = x1; Quark_Xs[index+1] << 328 Product *= (x1*x2*x3); << 329 } << 330 487 331 if ( Product == 0.0 ) continue; << 488 // Sampling X's of baryon >> 489 if ( Alfa_R == 1.0 ) { >> 490 G4double Xq1 = 1.0 - std::sqrt( G4UniformRand() ); >> 491 G4double Xq2 = (1.0 - Xq1) * G4UniformRand(); >> 492 G4double Xq3 = 1.0 - Xq1 - Xq2; >> 493 Quark_Mom[3].setZ( Xq1 ); Quark_Mom[4].setZ( Xq2 ); Quark_Mom[5].setZ( Xq3 ); >> 494 } else { >> 495 G4double Xq1 = sqr( G4UniformRand() ); >> 496 G4double Xq2 = (1.0 - Xq1) * sqr( std::sin( pi/2.0*G4UniformRand() ) ); >> 497 G4double Xq3 = 1.0 - Xq1 - Xq2; >> 498 Quark_Mom[3].setZ( Xq1 ); Quark_Mom[4].setZ( Xq2 ); Quark_Mom[5].setZ( Xq3 ); >> 499 } 332 500 333 ++NumberOfTries; << 501 G4double Alfa( 0.0 ), Beta( 0.0 ); 334 if ( NumberOfTries == 100*(NumberOfTries/1 << 502 for ( G4int i = 0; i < 3; i++ ) { // For Anti-baryon 335 // After a large number of tries, it is << 503 if ( Quark_Mom[i].getZ() != 0.0 ) { 336 ScaleFactor /= 2.0; << 504 Alfa += ( ScaleFactor * ModMom2[i] + MassQ2 ) / Quark_Mom[i].getZ(); 337 AveragePt2 *= ScaleFactor; << 505 } else { 338 } << 506 Succes = false; >> 507 } >> 508 } >> 509 for ( G4int i = 3; i < 6; i++ ) { // For baryon >> 510 if ( Quark_Mom[i].getZ() != 0.0 ) { >> 511 Beta += ( ScaleFactor * ModMom2[i] + MassQ2 ) / Quark_Mom[i].getZ(); >> 512 } else { >> 513 Succes = false; >> 514 } >> 515 } 339 516 340 G4ThreeVector PtSum( 0.0, 0.0, 0.0 ); << 517 if ( ! Succes ) continue; 341 for ( G4int i = 0; i < 6; ++i ) { << 342 Quark_Mom [i] = GaussianPt( AveragePt2, << 343 PtSum += Quark_Mom[i]; << 344 } << 345 518 346 PtSum /= 6.0; << 519 if ( std::sqrt( Alfa ) + std::sqrt( Beta ) > SqrtS ) { 347 Alfa = 0.0; Beta = 0.0; << 520 Succes = false; >> 521 continue; >> 522 } 348 523 349 for ( G4int i = 0; i < 6; ++i ) { // Loop << 524 G4double DecayMomentum2 = S*S + Alfa*Alfa + Beta*Beta 350 Quark_Mom[i] -= PtSum; << 525 - 2.0*S*Alfa - 2.0*S*Beta - 2.0*Alfa*Beta; >> 526 >> 527 WminusTarget = ( S - Alfa + Beta + std::sqrt( DecayMomentum2 ) ) / 2.0 / SqrtS; >> 528 WplusProjectile = SqrtS - Beta/WminusTarget; >> 529 >> 530 } while ( ( ! Succes ) && >> 531 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */ >> 532 if ( loopCounter >= maxNumberOfLoops ) { >> 533 return false; >> 534 } 351 535 352 G4double val = ( Quark_Mom[i].mag2() + M << 536 G4double SqrtScaleF = std::sqrt( ScaleFactor ); 353 if ( i < 3 ) { // anti-baryon << 537 for ( G4int i = 0; i < 3; i++ ) { 354 Alfa += val; << 538 G4double Pz = WplusProjectile * Quark_Mom[i].getZ() / 2.0 - 355 } else { // baryon (iCase == 1) << 539 ( ScaleFactor * ModMom2[i] + MassQ2 ) / 356 Beta += val; << 540 ( 2.0 * WplusProjectile * Quark_Mom[i].getZ() ); >> 541 Quark_Mom[i].setZ( Pz ); >> 542 if ( ScaleFactor != 1.0 ) { >> 543 Quark_Mom[i].setX( SqrtScaleF * Quark_Mom[i].getX() ); >> 544 Quark_Mom[i].setY( SqrtScaleF * Quark_Mom[i].getY() ); >> 545 } >> 546 } >> 547 for ( G4int i = 3; i < 6; i++ ) { >> 548 G4double Pz = -WminusTarget * Quark_Mom[i].getZ() / 2.0 + >> 549 ( ScaleFactor * ModMom2[i] + MassQ2 ) / >> 550 ( 2.0 * WminusTarget * Quark_Mom[i].getZ() ); >> 551 Quark_Mom[i].setZ( Pz ); >> 552 if ( ScaleFactor != 1.0 ) { >> 553 Quark_Mom[i].setX( SqrtScaleF * Quark_Mom[i].getX() ); >> 554 Quark_Mom[i].setY( SqrtScaleF * Quark_Mom[i].getY() ); 357 } 555 } 358 } 556 } >> 557 //G4cout << "Sum AQ " << Quark_Mom[0] + Quark_Mom[1] + Quark_Mom[2] << G4endl >> 558 // << "Sum Q " << Quark_Mom[3] + Quark_Mom[4] + Quark_Mom[5] << G4endl; 359 559 360 } while ( ( std::sqrt( Alfa ) + std::sqrt( B << 560 G4ThreeVector tmp = Quark_Mom[0] + Quark_Mom[3]; 361 ++loopCounter < maxNumberOfLoops ) << 561 G4LorentzVector Pstring1( tmp, std::sqrt( Quark_Mom[0].mag2() + MassQ2 ) + >> 562 std::sqrt( Quark_Mom[3].mag2() + MassQ2 ) ); >> 563 G4double Ystring1 = Pstring1.rapidity(); >> 564 >> 565 //G4cout << "Mom 1 string " << G4endl << Quark_Mom[0] << G4endl << Quark_Mom[3] << G4endl >> 566 // << tmp << " " << tmp.mag() << G4endl; >> 567 //G4cout << "1 str " << Pstring1 << " " << Pstring1.mag() << " " << Ystring1 << G4endl; >> 568 >> 569 tmp = Quark_Mom[1] + Quark_Mom[4]; >> 570 G4LorentzVector Pstring2( tmp, std::sqrt( Quark_Mom[1].mag2() + MassQ2 ) + >> 571 std::sqrt( Quark_Mom[4].mag2() + MassQ2 ) ); >> 572 G4double Ystring2 = Pstring2.rapidity(); >> 573 >> 574 //G4cout << "Mom 2 string " << G4endl << Quark_Mom[1] << G4endl << Quark_Mom[4] << G4endl >> 575 // << tmp << " " << tmp.mag() << G4endl; >> 576 //G4cout << "2 str " << Pstring2 << " " << Pstring2.mag() << " " << Ystring2 << G4endl; >> 577 >> 578 tmp = Quark_Mom[2] + Quark_Mom[5]; >> 579 G4LorentzVector Pstring3( tmp, std::sqrt( Quark_Mom[2].mag2() + MassQ2 ) + >> 580 std::sqrt( Quark_Mom[5].mag2() + MassQ2 ) ); >> 581 G4double Ystring3 = Pstring3.rapidity(); >> 582 >> 583 //G4cout << "Mom 3 string " << G4endl << Quark_Mom[2] << G4endl << Quark_Mom[5] << G4endl >> 584 // << tmp << " " << tmp.mag() << G4endl; >> 585 //G4cout << "3 str " << Pstring3 << " " << Pstring3.mag() << " " << Ystring3 << G4endl >> 586 // << "SumE " << Pstring1.e() + Pstring2.e() + Pstring3.e() << G4endl >> 587 // << Pstring1.mag() << " " <<Pstring2.mag() << " " << Pstring3.mag() << G4endl; >> 588 //G4int Uzhi; G4cin >> Uzhi; >> 589 >> 590 G4LorentzVector LeftString( 0.0, 0.0, 0.0, 0.0 ); >> 591 if ( Ystring1 > Ystring2 && Ystring2 > Ystring3 ) { >> 592 Pprojectile = Pstring1; >> 593 LeftString = Pstring2; >> 594 Ptarget = Pstring3; >> 595 } >> 596 if ( Ystring1 > Ystring3 && Ystring3 > Ystring2 ) { >> 597 Pprojectile = Pstring1; >> 598 LeftString = Pstring3; >> 599 Ptarget = Pstring2; >> 600 } >> 601 >> 602 if ( Ystring2 > Ystring1 && Ystring1 > Ystring3 ) { >> 603 Pprojectile = Pstring2; >> 604 LeftString = Pstring1; >> 605 Ptarget = Pstring3; >> 606 } >> 607 if ( Ystring2 > Ystring3 && Ystring3 > Ystring1 ) { >> 608 Pprojectile = Pstring2; >> 609 LeftString = Pstring3; >> 610 Ptarget = Pstring1; >> 611 } >> 612 >> 613 if ( Ystring3 > Ystring1 && Ystring1 > Ystring2 ) { >> 614 Pprojectile = Pstring3; >> 615 LeftString = Pstring1; >> 616 Ptarget = Pstring2; >> 617 } >> 618 if ( Ystring3 > Ystring2 && Ystring2 > Ystring1 ) { >> 619 Pprojectile = Pstring3; >> 620 LeftString = Pstring2; >> 621 Ptarget = Pstring1; >> 622 } >> 623 //G4cout << "SumP " << Pprojectile + LeftString + Ptarget << " " << SqrtS << G4endl; >> 624 >> 625 Pprojectile.transform( toLab ); >> 626 LeftString.transform( toLab ); >> 627 Ptarget.transform( toLab ); >> 628 //G4cout << "SumP " << Pprojectile + LeftString + Ptarget << " " << SqrtS << G4endl; 362 629 363 if ( loopCounter >= maxNumberOfLoops ) { << 630 // Calculation of the creation time 364 return false; << 631 projectile->SetTimeOfCreation( target->GetTimeOfCreation() ); 365 } << 632 projectile->SetPosition( target->GetPosition() ); >> 633 AdditionalString->SetTimeOfCreation( target->GetTimeOfCreation() ); >> 634 AdditionalString->SetPosition( target->GetPosition() ); >> 635 // Creation time and position of target nucleon were determined in >> 636 // ReggeonCascade() of G4FTFModel >> 637 >> 638 //G4cout << "Mproj " << Pprojectile.mag() << G4endl << "Mtarg " << Ptarget.mag() << G4endl; >> 639 projectile->Set4Momentum( Pprojectile ); >> 640 AdditionalString->Set4Momentum( LeftString ); >> 641 target->Set4Momentum( Ptarget ); >> 642 projectile->IncrementCollisionCount( 1 ); >> 643 AdditionalString->IncrementCollisionCount( 1 ); >> 644 target->IncrementCollisionCount( 1 ); 366 645 367 G4double DecayMomentum2 = sqr(common.S) + sq << 646 theParameters->SetProbabilityOfAnnihilation( 0.0 ); // Uzhi March 2016 368 - 2.0*( common.S*( << 369 647 370 G4double WminusTarget = 0.0, WplusProjectile << 648 return true; 371 WminusTarget = ( common.S - Alfa + Beta + st << 372 WplusProjectile = common.SqrtS - Beta/Wminus << 373 << 374 for ( G4int iCase = 0; iCase < 2; ++iCase ) << 375 G4int index = iCase*3; // << 376 G4double w = WplusProjectile; // << 377 if ( iCase == 1 ) w = - WminusTarget; // << 378 for ( G4int i = 0; i < 3; ++i ) { << 379 G4double Pz = w * Quark_Xs[index+i] / 2. << 380 ( Quark_Mom[index+i].mag2( << 381 ( 2.0 * w * Quark_Xs[index << 382 Quark_Mom[index+i].setZ( Pz ); << 383 } << 384 } << 385 649 386 // Sampling of anti-quark order in projectil << 650 } // End of if ( Ksi < X_a / Xannihilation ) 387 G4int SampledCase = (G4int)G4RandFlat::shoot << 388 G4int Tmp1 = 0, Tmp2 = 0; << 389 switch ( SampledCase ) { << 390 case 1 : Tmp1 = common.AQ[1]; common.AQ[1] << 391 case 2 : Tmp1 = common.AQ[0]; common.AQ[0] << 392 case 3 : Tmp1 = common.AQ[0]; Tmp2 << 393 common.AQ[1] = Tmp1; comm << 394 case 4 : Tmp1 = common.AQ[0]; Tmp2 << 395 common.AQ[1] = common.AQ[2]; comm << 396 case 5 : Tmp1 = common.AQ[0]; Tmp2 << 397 common.AQ[1] = Tmp2; comm << 398 } << 399 651 400 // Set the string properties << 652 // Simulation of anti-diquark-diquark string creation 401 // An anti quark - quark pair can have the q << 402 // or a vector meson: the last digit of the << 403 // For simplicity only scalar is considered << 404 G4int NewCode = 0, antiQuark = 0, quark = 0; << 405 G4ParticleDefinition* TestParticle = nullptr << 406 for ( G4int iString = 0; iString < 3; ++iStr << 407 if ( iString == 0 ) { << 408 antiQuark = common.AQ[0]; quark = commo << 409 projectile->SetFirstParton( antiQuark ); << 410 projectile->SetSecondParton( quark ); << 411 projectile->SetStatus( 0 ); << 412 } else if ( iString == 1 ) { << 413 quark = common.Q[1]; antiQuark = common << 414 target->SetFirstParton( quark ); << 415 target->SetSecondParton( antiQuark ); << 416 target->SetStatus( 0 ); << 417 } else { // iString == 2 << 418 antiQuark = common.AQ[2]; quark = commo << 419 } << 420 G4int absAntiQuark = std::abs( antiQuark ) << 421 G4double aKsi = G4UniformRand(); << 422 if ( absAntiQuark == absQuark ) { << 423 if ( absAntiQuark != 3 ) { // Not yet c << 424 NewCode = 111; // Pi0-meson << 425 if ( aKsi < 0.5 ) { << 426 NewCode = 221; // Eta -meso << 427 if ( aKsi < 0.25 ) { << 428 NewCode = 331; // Eta'-meso << 429 } << 430 } << 431 } else { << 432 NewCode = 221; // Eta -meso << 433 if ( aKsi < 0.5 ) { << 434 NewCode = 331; // Eta'-meso << 435 } << 436 } << 437 } else { // Vector mesons - rho, omega, p << 438 if ( absAntiQuark > absQuark ) { << 439 NewCode = absAntiQuark*100 + absQuark* << 440 } else { << 441 NewCode = absQuark*100 + absAntiQuark* << 442 } << 443 } << 444 if ( iString == 2 ) AdditionalString = new << 445 TestParticle = G4ParticleTable::GetParticl << 446 if ( ! TestParticle ) return false; << 447 if ( iString == 0 ) { << 448 projectile->SetDefinition( TestParticle << 449 theParameters->SetProjMinDiffMass( 0.5 ) << 450 theParameters->SetProjMinNonDiffMass( 0. << 451 } else if ( iString == 1 ) { << 452 target->SetDefinition( TestParticle ); << 453 theParameters->SetTarMinDiffMass( 0.5 ); << 454 theParameters->SetTarMinNonDiffMass( 0.5 << 455 } else { // iString == 2 << 456 AdditionalString->SetDefinition( TestPar << 457 AdditionalString->SetFirstParton( common << 458 AdditionalString->SetSecondParton( commo << 459 AdditionalString->SetStatus( 0 ); << 460 } << 461 } // End of the for loop over the 3 string << 462 653 463 // 1st string AQ[0]-Q[0], 2nd string AQ[1]-Q << 654 if ( Ksi < (X_a + X_b) / Xannihilation ) { >> 655 >> 656 #ifdef debugFTFannih >> 657 G4cout << "Process b, quark - anti-quark annihilation, di-q - anti-di-q string" << G4endl; >> 658 #endif 464 659 465 G4LorentzVector Pstring1, Pstring2, Pstring3 << 660 G4int CandidatsN( 0 ), CandAQ[9][2], CandQ[9][2]; 466 G4int QuarkOrder[3] = { 0 }; << 661 G4int LeftAQ1( 0 ), LeftAQ2( 0 ), LeftQ1( 0 ), LeftQ2( 0 ); 467 G4double YstringMax = 0.0, YstringMin = 0.0; << 662 468 for ( G4int i = 0; i < 3; ++i ) { << 663 for ( G4int iAQ = 0; iAQ < 3; iAQ++ ) { 469 G4ThreeVector tmp = Quark_Mom[i] + Quark_M << 664 for ( G4int iQ = 0; iQ < 3; iQ++ ) { 470 G4LorentzVector Pstring( tmp, std::sqrt( Q << 665 if ( -AQ[iAQ] == Q[iQ] ) { 471 std::sqrt( Q << 666 if ( iAQ == 0 ) { CandAQ[CandidatsN][0] = 1; CandAQ[CandidatsN][1] = 2; } 472 // Add protection for rapidity = 0.5*ln( << 667 if ( iAQ == 1 ) { CandAQ[CandidatsN][0] = 0; CandAQ[CandidatsN][1] = 2; } 473 G4double Ystring = 0.0; << 668 if ( iAQ == 2 ) { CandAQ[CandidatsN][0] = 0; CandAQ[CandidatsN][1] = 1; } 474 if ( Pstring.e() > 1.0e-30 ) { << 669 if ( iQ == 0 ) { CandQ[CandidatsN][0] = 1; CandQ[CandidatsN][1] = 2; } 475 if ( Pstring.e() + Pstring.pz() < 1.0e-3 << 670 if ( iQ == 1 ) { CandQ[CandidatsN][0] = 0; CandQ[CandidatsN][1] = 2; } 476 Ystring = -1.0e30; // A very large n << 671 if ( iQ == 2 ) { CandQ[CandidatsN][0] = 0; CandQ[CandidatsN][1] = 1; } 477 if ( Pstring.e() - Pstring.pz() < 1.0e << 672 CandidatsN++; 478 Ystring = 1.0e30; // A very large p << 479 } else { // Normal case << 480 Ystring = Pstring.rapidity(); << 481 } 673 } 482 } 674 } 483 } 675 } 484 // Keep ordering in rapidity: "1" highest, << 676 //G4cout << "CandidatsN " << CandidatsN << G4endl; 485 if ( i == 0 ) { << 677 486 Pstring1 = Pstring; YstringMax = Yst << 678 if ( CandidatsN != 0 ) { 487 QuarkOrder[0] = 0; << 679 G4int SampledCase = G4RandFlat::shootInt( G4long( CandidatsN ) ); 488 } else if ( i == 1 ) { << 680 LeftAQ1 = AQ[ CandAQ[SampledCase][0] ]; 489 if ( Ystring > YstringMax ) { << 681 LeftAQ2 = AQ[ CandAQ[SampledCase][1] ]; 490 Pstring2 = Pstring1; YstringMin = Yst << 682 LeftQ1 = Q[ CandQ[SampledCase][0] ]; 491 Pstring1 = Pstring; YstringMax = Yst << 683 LeftQ2 = Q[ CandQ[SampledCase][1] ]; 492 QuarkOrder[0] = 1; QuarkOrder[1] = 0; << 684 >> 685 // Build anti-diquark and diquark >> 686 G4int Anti_DQ( 0 ), DQ( 0 ); >> 687 if ( std::abs( LeftAQ1 ) > std::abs( LeftAQ2 ) ) { >> 688 Anti_DQ = 1000*LeftAQ1 + 100*LeftAQ2 - 3; // 1 493 } else { 689 } else { 494 Pstring2 = Pstring; YstringMin = Yst << 690 Anti_DQ = 1000*LeftAQ2 + 100*LeftAQ1 - 3; // 1 495 QuarkOrder[1] = 1; << 496 } 691 } 497 } else { // i == 2 << 692 //if ( G4UniformRand() > 0.5 ) Anti_DQ -= 2; 498 if ( Ystring > YstringMax ) { << 693 if ( std::abs( LeftQ1 ) > std::abs( LeftQ2 ) ) { 499 Pstring3 = Pstring2; << 694 DQ = 1000*LeftQ1 + 100*LeftQ2 + 3; // 1 500 Pstring2 = Pstring1; << 501 Pstring1 = Pstring; << 502 QuarkOrder[1] = QuarkOrder[0]; << 503 QuarkOrder[2] = QuarkOrder[1]; << 504 QuarkOrder[0] = 2; << 505 } else if ( Ystring > YstringMin ) { << 506 Pstring3 = Pstring2; << 507 Pstring2 = Pstring; << 508 } else { 695 } else { 509 Pstring3 = Pstring; << 696 DQ = 1000*LeftQ2 + 100*LeftQ1 + 3; // 1 510 QuarkOrder[2] = 2; << 511 } 697 } 512 } << 698 // if ( G4UniformRand() > 0.5 ) DQ += 2; 513 } << 514 << 515 G4LorentzVector Quark_4Mom[6]; << 516 for ( G4int i = 0; i < 6; ++i ) { << 517 Quark_4Mom[i] = G4LorentzVector( Quark_Mom << 518 if ( common.RotateStrings ) Quark_4Mom[i] << 519 Quark_4Mom[i].transform( common.toLab ); << 520 } << 521 699 522 projectile->Splitting(); << 700 // Set the string properties 523 projectile->GetNextAntiParton()->Set4Momentu << 701 //G4cout << "Left ADiQ DiQ " << Anti_DQ << " " << DQ << G4endl; 524 projectile->GetNextParton()->Set4Momentum( Q << 702 projectile->SplitUp(); 525 << 703 //projectile->SetFirstParton( Anti_DQ ); 526 target->Splitting(); << 704 //projectile->SetSecondParton( DQ ); 527 target->GetNextParton()->Set4Momentum( Quark << 705 projectile->SetFirstParton( DQ ); 528 target->GetNextAntiParton()->Set4Momentum( Q << 706 projectile->SetSecondParton( Anti_DQ ); 529 << 707 projectile->SetStatus( 0 ); 530 AdditionalString->Splitting(); << 708 target->SetStatus( 4 ); // The target nucleon has annihilated 3->4 Uzhi Oct 2014 531 AdditionalString->GetNextAntiParton()->Set4M << 709 Pprojectile.setPx( 0.0 ); 532 AdditionalString->GetNextParton()->Set4Momen << 710 Pprojectile.setPy( 0.0 ); 533 << 711 Pprojectile.setPz( 0.0 ); 534 common.Pprojectile = Pstring1; // << 712 Pprojectile.setE( SqrtS ); 535 common.Ptarget = Pstring3; // << 713 Pprojectile.transform( toLab ); 536 G4LorentzVector LeftString( Pstring2 ); // << 714 // Uzhi March 2016 if QQ_QQbar will interact Set Mmin, MdifMin 537 << 715 538 if ( common.RotateStrings ) { << 716 // Calculation of the creation time 539 common.Pprojectile *= common.RandomRotatio << 717 projectile->SetTimeOfCreation( target->GetTimeOfCreation() ); 540 common.Ptarget *= common.RandomRotatio << 718 projectile->SetPosition( target->GetPosition() ); 541 LeftString *= common.RandomRotatio << 719 // Creation time and position of target nucleon were determined in 542 } << 720 // ReggeonCascade() of G4FTFModel 543 << 721 544 common.Pprojectile.transform( common.toLab ) << 722 //G4cout << "Mproj " << Pprojectile.mag() << G4endl 545 common.Ptarget.transform( common.toLab ); << 723 // << "Mtarg " << Ptarget.mag() << G4endl; 546 LeftString.transform( common.toLab ); << 724 projectile->Set4Momentum( Pprojectile ); 547 << 548 // Calculation of the creation time << 549 // Creation time and position of target nucl << 550 projectile->SetTimeOfCreation( target->GetTi << 551 projectile->SetPosition( target->GetPosition << 552 AdditionalString->SetTimeOfCreation( target- << 553 AdditionalString->SetPosition( target->GetPo << 554 << 555 projectile->Set4Momentum( common.Pprojectile << 556 AdditionalString->Set4Momentum( LeftString ) << 557 target->Set4Momentum( common.Ptarget ); << 558 << 559 projectile->IncrementCollisionCount( 1 ); << 560 AdditionalString->IncrementCollisionCount( 1 << 561 target->IncrementCollisionCount( 1 ); << 562 << 563 return true; << 564 } << 565 << 566 << 567 //-------------------------------------------- << 568 << 569 G4int G4FTFAnnihilation:: << 570 Create1DiquarkAntiDiquarkString( G4VSplitableH << 571 G4VSplitableH << 572 G4FTFAnnihila << 573 // Simulation of anti-diquark-diquark string << 574 // This method returns an integer code - ins << 575 // "0" : successfully ended and nothing el << 576 // "1" : successfully completed, but the w << 577 // "99" : unsuccessfully ended, nothing els << 578 725 579 #ifdef debugFTFannih << 726 projectile->IncrementCollisionCount( 1 ); 580 G4cout << "Process b, quark - anti-quark ann << 727 target->IncrementCollisionCount( 1 ); 581 #endif << 582 728 583 G4int CandidatsN = 0, CandAQ[9][2] = {}, Can << 729 //theParameters->SetProbabilityOfAnnihilation( 0.0 ); // Uzhi March 2016 584 for ( G4int iAQ = 0; iAQ < 3; ++iAQ ) { // << 730 // In the case baryon and anti-baryon are created. Thus the antibaryon can annihilate latter. 585 for ( G4int iQ = 0; iQ < 3; ++iQ ) { // << 586 if ( -common.AQ[iAQ] == common.Q[iQ] ) { << 587 // Here "0", "1", "2" means, respectiv << 588 // of the (anti-baryon) projectile or << 589 if ( iAQ == 0 ) { CandAQ[CandidatsN][0 << 590 if ( iAQ == 1 ) { CandAQ[CandidatsN][0 << 591 if ( iAQ == 2 ) { CandAQ[CandidatsN][0 << 592 if ( iQ == 0 ) { CandQ[CandidatsN][0] << 593 if ( iQ == 1 ) { CandQ[CandidatsN][0] << 594 if ( iQ == 2 ) { CandQ[CandidatsN][0] << 595 ++CandidatsN; << 596 } << 597 } << 598 } << 599 731 600 // Remaining two (anti-)quarks that form the << 732 return true; 601 G4int LeftAQ1 = 0, LeftAQ2 = 0, LeftQ1 = 0, << 602 if ( CandidatsN != 0 ) { << 603 G4int SampledCase = (G4int)G4RandFlat::sho << 604 LeftAQ1 = common.AQ[ CandAQ[SampledCase][0 << 605 LeftAQ2 = common.AQ[ CandAQ[SampledCase][1 << 606 LeftQ1 = common.Q[ CandQ[SampledCase][0] << 607 LeftQ2 = common.Q[ CandQ[SampledCase][1] << 608 << 609 // Build anti-diquark and diquark : the la << 610 // of anti-quark - anti-quark and quark - << 611 // or quarks are different. For simplicity << 612 G4int Anti_DQ = 0, DQ = 0; << 613 if ( std::abs( LeftAQ1 ) > std::abs( LeftA << 614 Anti_DQ = 1000*LeftAQ1 + 100*LeftAQ2 - 3 << 615 } else { << 616 Anti_DQ = 1000*LeftAQ2 + 100*LeftAQ1 - 3 << 617 } << 618 if ( std::abs( LeftQ1 ) > std::abs( LeftQ2 << 619 DQ = 1000*LeftQ1 + 100*LeftQ2 + 3; << 620 } else { << 621 DQ = 1000*LeftQ2 + 100*LeftQ1 + 3; << 622 } 733 } 623 734 624 // Set the string properties << 735 } // End of if ( Ksi < (X_a + X_b) / Xannihilation ) 625 projectile->SetFirstParton( DQ ); << 626 projectile->SetSecondParton( Anti_DQ ); << 627 << 628 // It is assumed that quark and di-quark m << 629 G4LorentzVector Pquark = G4LorentzVector( << 630 G4LorentzVector Paquark = G4LorentzVector( << 631 << 632 if ( common.RotateStrings ) { << 633 Pquark *= common.RandomRotation; << 634 Paquark *= common.RandomRotation; << 635 } << 636 736 637 Pquark.transform( common.toLab ); << 737 if ( Ksi < ( X_a + X_b + X_c ) / Xannihilation ) { 638 Paquark.transform( common.toLab ); << 639 738 640 projectile->GetNextParton()->Set4Momentum( << 739 // Simulation of 2 anti-quark-quark strings creation 641 projectile->GetNextAntiParton()->Set4Momen << 642 740 643 projectile->Splitting(); << 741 #ifdef debugFTFannih >> 742 G4cout << "Process c, quark - anti-quark and string junctions annihilation, 2 strings left." >> 743 << G4endl; >> 744 #endif 644 745 645 projectile->SetStatus( 0 ); << 746 G4int CandidatsN( 0 ), CandAQ[9][2], CandQ[9][2]; 646 target->SetStatus( 4 ); // The target nuc << 747 G4int LeftAQ1( 0 ), LeftAQ2( 0 ), LeftQ1( 0 ), LeftQ2( 0 ); 647 common.Pprojectile.setPx( 0.0 ); << 648 common.Pprojectile.setPy( 0.0 ); << 649 common.Pprojectile.setPz( 0.0 ); << 650 common.Pprojectile.setE( common.SqrtS ); << 651 common.Pprojectile.transform( common.toLab << 652 748 653 // Calculation of the creation time << 749 for ( G4int iAQ = 0; iAQ < 3; iAQ++ ) { 654 // Creation time and position of target nu << 750 for ( G4int iQ = 0; iQ < 3; iQ++ ) { 655 projectile->SetTimeOfCreation( target->Get << 751 if ( -AQ[iAQ] == Q[iQ] ) { 656 projectile->SetPosition( target->GetPositi << 752 if ( iAQ == 0 ) { CandAQ[CandidatsN][0] = 1; CandAQ[CandidatsN][1] = 2; } 657 projectile->Set4Momentum( common.Pprojecti << 753 if ( iAQ == 1 ) { CandAQ[CandidatsN][0] = 0; CandAQ[CandidatsN][1] = 2; } >> 754 if ( iAQ == 2 ) { CandAQ[CandidatsN][0] = 0; CandAQ[CandidatsN][1] = 1; } >> 755 if ( iQ == 0 ) { CandQ[CandidatsN][0] = 1; CandQ[CandidatsN][1] = 2; } >> 756 if ( iQ == 1 ) { CandQ[CandidatsN][0] = 0; CandQ[CandidatsN][1] = 2; } >> 757 if ( iQ == 2 ) { CandQ[CandidatsN][0] = 0; CandQ[CandidatsN][1] = 1; } >> 758 CandidatsN++; >> 759 } >> 760 } >> 761 } >> 762 //G4cout << "CandidatsN " << CandidatsN << G4endl; 658 763 659 projectile->IncrementCollisionCount( 1 ); << 764 if ( CandidatsN != 0 ) { 660 target->IncrementCollisionCount( 1 ); << 765 G4int SampledCase = G4RandFlat::shootInt( G4long( CandidatsN ) ); >> 766 LeftAQ1 = AQ[ CandAQ[SampledCase][0] ]; >> 767 LeftAQ2 = AQ[ CandAQ[SampledCase][1] ]; >> 768 if ( G4UniformRand() < 0.5 ) { >> 769 LeftQ1 = Q[ CandQ[SampledCase][0] ]; >> 770 LeftQ2 = Q[ CandQ[SampledCase][1] ]; >> 771 } else { >> 772 LeftQ2 = Q[ CandQ[SampledCase][0] ]; >> 773 LeftQ1 = Q[ CandQ[SampledCase][1] ]; >> 774 } 661 775 662 return 0; // Completed successfully: noth << 776 // Set the string properties 663 } // End of if ( CandidatsN != 0 ) << 777 //G4cout << "String 1 " << LeftAQ1 << " " << LeftQ1 << G4endl; >> 778 projectile->SplitUp(); >> 779 projectile->SetFirstParton( LeftAQ1 ); >> 780 projectile->SetSecondParton( LeftQ1 ); >> 781 projectile->SetStatus( 0 ); 664 782 665 // If we allow the string to interact with o << 783 // Uzhi March 2016 start 666 // set up MinDiffrMass in Parameters, and as << 784 G4int aAQ, aQ; 667 << 785 aAQ=std::abs( LeftAQ1 ); aQ=std::abs( LeftQ1 ); 668 return 1; // Successfully ended, but the wo << 786 >> 787 G4int NewCode; >> 788 G4double aKsi = G4UniformRand(); >> 789 >> 790 if ( aAQ == aQ ) >> 791 { >> 792 if ( aAQ != 3 ) >> 793 { >> 794 NewCode = 111; // Pi0-meson >> 795 if ( aKsi < 0.5 ) >> 796 { >> 797 NewCode = 221; // Eta -meson >> 798 if ( aKsi < 0.25 ) {NewCode = 331;} // Eta'-meson >> 799 } >> 800 } else >> 801 { >> 802 NewCode = 221; // Eta -meson >> 803 if( aKsi < 0.5 ) {NewCode = 331;} // Eta'-meson >> 804 } >> 805 } else >> 806 { >> 807 if ( aAQ > aQ ){ NewCode = aAQ*100 + aQ*10 + 1; NewCode *= aAQ/LeftAQ1; } >> 808 else { NewCode = aQ*100 + aAQ*10 + 1; NewCode *= aQ/LeftQ1; } 669 } 809 } 670 810 >> 811 G4ParticleDefinition* TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewCode ); >> 812 if(!TestParticle) return false; >> 813 projectile->SetDefinition( TestParticle ); >> 814 theParameters->SetProjMinDiffMass( 0.5 ); // (0.5) // GeV Uzhi March 2016 ??? >> 815 theParameters->SetProjMinNonDiffMass( 0.5 ); >> 816 // Uzhi March 2016 end >> 817 >> 818 //G4cout << "String 2 " << LeftAQ2 << " " << LeftQ2 << G4endl; >> 819 target->SplitUp(); >> 820 target->SetFirstParton( LeftQ2 ); >> 821 target->SetSecondParton( LeftAQ2 ); >> 822 target->SetStatus( 0 ); 671 823 672 //-------------------------------------------- << 824 // Uzhi March 2016 start 673 << 825 aAQ=std::abs( LeftAQ2 ); aQ=std::abs( LeftQ2 ); aKsi = G4UniformRand(); 674 G4int G4FTFAnnihilation:: << 675 Create2QuarkAntiQuarkStrings( G4VSplitableHadr << 676 G4VSplitableHadr << 677 G4FTFParameters* << 678 G4FTFAnnihilatio << 679 // Simulation of 2 anti-quark-quark strings << 680 // This method returns an integer code - ins << 681 // "0" : successfully ended and nothing el << 682 // "1" : successfully completed, but the w << 683 // "99" : unsuccessfully ended, nothing els << 684 826 685 #ifdef debugFTFannih << 827 if ( aAQ == aQ ) 686 G4cout << "Process c, quark - anti-quark and << 828 { 687 << G4endl; << 829 if ( aAQ != 3 ) 688 #endif << 830 { >> 831 NewCode = 111; // Pi0-meson >> 832 if ( aKsi < 0.5 ) >> 833 { >> 834 NewCode = 221; // Eta -meson >> 835 if ( aKsi < 0.25 ) {NewCode = 331;} // Eta'-meson >> 836 } >> 837 } else >> 838 { >> 839 NewCode = 221; // Eta -meson >> 840 if( aKsi < 0.5 ) {NewCode = 331;} // Eta'-meson >> 841 } >> 842 } else >> 843 { >> 844 if ( aAQ > aQ ){ NewCode = aAQ*100 + aQ*10 + 1; NewCode *= aAQ/LeftAQ2; } >> 845 else { NewCode = aQ*100 + aAQ*10 + 1; NewCode *= aQ/LeftQ2; } >> 846 } 689 847 690 // Sampling kinematical properties: 1st stri << 848 TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewCode ); 691 G4ThreeVector Quark_Mom[4]; << 849 if(!TestParticle) return false; 692 G4double Quark_Xs[4]; << 850 target->SetDefinition( TestParticle ); 693 G4double AveragePt2 = 200.0*200.0, maxPtSqua << 851 theParameters->SetTarMinDiffMass( 0.5 ); // Uzhi March 2016 ??? 694 G4int NumberOfTries = 0, loopCounter = 0; << 852 theParameters->SetTarMinNonDiffMass( 0.5 ); 695 const G4int maxNumberOfLoops = 1000; << 853 // Uzhi March 2016 696 G4double Alfa = 0.0, Beta = 0.0; << 854 697 G4double WminusTarget = 0.0, WplusProjectile << 855 // Sampling kinematical properties 698 do { << 856 // 1 string LeftAQ1-LeftQ1// 2 string LeftAQ2-LeftQ2 699 // Sampling X's of the 2 quarks and 2 anti << 857 G4ThreeVector Quark_Mom[4]; 700 << 858 G4double ModMom2[4]; //ModMom[4], 701 G4double Product = 1.0; << 859 702 for ( G4int iCase = 0; iCase < 2; ++iCase << 860 AveragePt2 = 200.0*200.0; maxPtSquare = S; 703 G4double x = 0.0, r = G4UniformRand(); << 861 704 if ( Alfa_R == 1.0 ) { << 862 G4double SumMt( 0.0 ); 705 if ( iCase == 0 ) { // first string << 863 G4double MassQ2 = 0.0; //100.0*100.0*MeV*MeV; 706 x = std::sqrt( r ); << 864 G4int NumberOfTries( 0 ); 707 } else { // second string << 865 G4double ScaleFactor( 1.0 ); 708 x = 1.0 - std::sqrt( r ); << 866 >> 867 const G4int maxNumberOfLoops = 1000; >> 868 G4int loopCounter = 0; >> 869 do { >> 870 NumberOfTries++; >> 871 if ( NumberOfTries == 100*(NumberOfTries/100) ) { >> 872 // At large number of tries it would be better to reduce the values of <Pt^2> >> 873 ScaleFactor /= 2.0; >> 874 AveragePt2 *= ScaleFactor; 709 } 875 } 710 } else { << 876 G4ThreeVector PtSum( 0.0, 0.0, 0.0 ); 711 x = sqr( std::sin( pi/2.0*r ) ); << 877 for( G4int i = 0; i < 4; i++ ) { >> 878 Quark_Mom[i] = GaussianPt( AveragePt2, maxPtSquare ); >> 879 PtSum += Quark_Mom[i]; >> 880 } >> 881 PtSum /= 4.0; >> 882 SumMt = 0.0; >> 883 for ( G4int i = 0; i < 4; i++ ) { >> 884 Quark_Mom[i] -= PtSum; >> 885 //ModMom[i] = Quark_Mom[i].mag(); >> 886 ModMom2[i] = Quark_Mom[i].mag2(); >> 887 SumMt += std::sqrt( ModMom2[i] + MassQ2 ); >> 888 } >> 889 } while ( ( SumMt > SqrtS ) && >> 890 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */ >> 891 if ( loopCounter >= maxNumberOfLoops ) { >> 892 return false; 712 } 893 } 713 G4int index = iCase*2; // 0 for the fir << 714 Quark_Xs[index] = x ; Quark_Xs[index+1] << 715 Product *= x*(1.0-x); << 716 } << 717 894 718 if ( Product == 0.0 ) continue; << 895 G4double WminusTarget( 0.0 ), WplusProjectile( 0.0 ); 719 896 720 ++NumberOfTries; << 897 // Sampling X's of anti-baryon 721 if ( NumberOfTries == 100*(NumberOfTries/1 << 898 G4double Alfa_R = 0.5; 722 // After a large number of tries, it is << 899 NumberOfTries = 0; 723 ScaleFactor /= 2.0; << 900 ScaleFactor = 1.0; 724 AveragePt2 *= ScaleFactor; << 901 G4bool Succes( true ); 725 } << 902 >> 903 loopCounter = 0; >> 904 do { >> 905 >> 906 Succes = true; >> 907 NumberOfTries++; >> 908 if ( NumberOfTries == 100*(NumberOfTries/100) ) { >> 909 // At large number of tries it would be better to reduce the values of Pt's >> 910 ScaleFactor /= 2.0; >> 911 } 726 912 727 G4ThreeVector PtSum( 0.0, 0.0, 0.0 ); << 913 if ( Alfa_R == 1.0 ) { 728 for( G4int i = 0; i < 4; ++i ) { << 914 G4double Xaq1 = std::sqrt( G4UniformRand() ); 729 Quark_Mom[i] = GaussianPt( AveragePt2, m << 915 G4double Xaq2 = 1.0 - Xaq1; 730 PtSum += Quark_Mom[i]; << 916 Quark_Mom[0].setZ( Xaq1 ); Quark_Mom[1].setZ( Xaq2 ); 731 } << 917 } else { >> 918 G4double Xaq1 = sqr( std::sin( pi/2.0*G4UniformRand() ) ); >> 919 G4double Xaq2 = 1.0 - Xaq1; >> 920 Quark_Mom[0].setZ( Xaq1 ); Quark_Mom[1].setZ( Xaq2 ); >> 921 } 732 922 733 PtSum /= 4.0; << 923 // Sampling X's of baryon ------------ 734 for ( G4int i = 0; i < 4; ++i ) { << 924 if ( Alfa_R == 1.0 ) { 735 Quark_Mom[i] -= PtSum; << 925 G4double Xq1 = 1.0 - std::sqrt( G4UniformRand() ); 736 } << 926 G4double Xq2 = 1.0 - Xq1; >> 927 Quark_Mom[2].setZ( Xq1 ); Quark_Mom[3].setZ( Xq2 ); >> 928 } else { >> 929 G4double Xq1 = sqr( std::sin( pi/2.0*G4UniformRand() ) ); >> 930 G4double Xq2 = 1.0 - Xq1; >> 931 Quark_Mom[2].setZ( Xq1 ); Quark_Mom[3].setZ( Xq2 ); >> 932 } 737 933 738 Alfa = 0.0; Beta = 0.0; << 934 G4double Alfa( 0.0 ), Beta( 0.0 ); 739 for ( G4int iCase = 0; iCase < 2; ++iCase << 935 for ( G4int i = 0; i < 2; i++ ) { // For Anti-baryon 740 G4int index = iCase * 2; << 936 if ( Quark_Mom[i].getZ() != 0.0 ) { 741 for ( G4int i = 0; i < 2; ++i ) { << 937 Alfa += ( ScaleFactor * ModMom2[i] + MassQ2 ) / Quark_Mom[i].getZ(); 742 G4double val = ( Quark_Mom[index+i]. << 938 } else { 743 if ( iCase == 0 ) { // first string << 939 Succes = false; 744 Alfa += val; << 745 } else { // second strin << 746 Beta += val; << 747 } 940 } 748 } << 941 } 749 } << 942 for ( G4int i = 2; i < 4; i++ ) { // For baryon 750 << 943 if ( Quark_Mom[i].getZ() != 0.0 ) { 751 } while ( ( std::sqrt( Alfa ) + std::sqrt( B << 944 Beta += ( ScaleFactor * ModMom2[i] + MassQ2 ) / Quark_Mom[i].getZ(); 752 ++loopCounter < maxNumberOfLoops ) << 945 } else { >> 946 Succes = false; >> 947 } >> 948 } 753 949 754 if ( loopCounter >= maxNumberOfLoops ) { << 950 if ( ! Succes ) continue; 755 return 99; // unsuccessfully ended, nothi << 756 } << 757 951 758 G4double DecayMomentum2 = sqr(common.S) + sq << 952 if ( std::sqrt( Alfa ) + std::sqrt( Beta ) > SqrtS ) { 759 - 2.0*( common.S*( << 953 Succes = false; 760 WminusTarget = ( common.S - Alfa + Beta + st << 954 continue; 761 WplusProjectile = common.SqrtS - Beta/Wminus << 955 } 762 << 763 for ( G4int iCase = 0; iCase < 2; ++iCase ) << 764 G4int index = iCase*2; // 0 for the first << 765 for ( G4int i = 0; i < 2; ++i ) { << 766 G4double w = WplusProjectile; / << 767 if ( iCase == 1 ) w = - WminusTarget; / << 768 G4double Pz = w * Quark_Xs[index+i] / 2. << 769 - ( Quark_Mom[index+i].mag << 770 ( 2.0 * w * Quark_Xs[ind << 771 Quark_Mom[index+i].setZ( Pz ); << 772 } << 773 } << 774 956 775 G4int CandidatsN = 0, CandAQ[9][2] = {}, Can << 957 G4double DecayMomentum2 = S*S + Alfa*Alfa + Beta*Beta 776 G4int LeftAQ1 = 0, LeftAQ2 = 0, LeftQ1 = 0, << 958 - 2.0*S*Alfa - 2.0*S*Beta - 2.0*Alfa*Beta; 777 for ( G4int iAQ = 0; iAQ < 3; ++iAQ ) { // << 959 WminusTarget = ( S - Alfa + Beta + std::sqrt( DecayMomentum2 ) ) / 2.0 / SqrtS; 778 for ( G4int iQ = 0; iQ < 3; ++iQ ) { // << 960 WplusProjectile = SqrtS - Beta/WminusTarget; 779 if ( -common.AQ[iAQ] == common.Q[iQ] ) { << 961 780 // Here "0", "1", "2" means, respectiv << 962 } while ( ( ! Succes ) && 781 // of the (anti-baryon) projectile or << 963 ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */ 782 if ( iAQ == 0 ) { CandAQ[CandidatsN][0 << 964 if ( loopCounter >= maxNumberOfLoops ) { 783 if ( iAQ == 1 ) { CandAQ[CandidatsN][0 << 965 return false; 784 if ( iAQ == 2 ) { CandAQ[CandidatsN][0 << 785 if ( iQ == 0 ) { CandQ[CandidatsN][0] << 786 if ( iQ == 1 ) { CandQ[CandidatsN][0] << 787 if ( iQ == 2 ) { CandQ[CandidatsN][0] << 788 ++CandidatsN; << 789 } 966 } 790 } << 791 } << 792 967 793 if ( CandidatsN != 0 ) { << 968 G4double SqrtScaleF = std::sqrt( ScaleFactor ); 794 G4int SampledCase = (G4int)G4RandFlat::sho << 795 LeftAQ1 = common.AQ[ CandAQ[SampledCase][0 << 796 LeftAQ2 = common.AQ[ CandAQ[SampledCase][1 << 797 if ( G4UniformRand() < 0.5 ) { << 798 LeftQ1 = common.Q[ CandQ[SampledCase][0] << 799 LeftQ2 = common.Q[ CandQ[SampledCase][1] << 800 } else { << 801 LeftQ2 = common.Q[ CandQ[SampledCase][0] << 802 LeftQ1 = common.Q[ CandQ[SampledCase][1] << 803 } << 804 969 805 // Set the string properties << 970 for ( G4int i = 0; i < 2; i++ ) { 806 // An anti quark - quark pair can have the << 971 G4double Pz = WplusProjectile * Quark_Mom[i].getZ() / 2.0 - 807 // or a vector meson: the last digit of th << 972 ( ScaleFactor * ModMom2[i] + MassQ2 ) / 808 // For simplicity only scalar is considere << 973 ( 2.0 * WplusProjectile * Quark_Mom[i].getZ() ); 809 G4int NewCode = 0, antiQuark = 0, quark = << 974 Quark_Mom[i].setZ( Pz ); 810 G4ParticleDefinition* TestParticle = nullp << 975 if ( ScaleFactor != 1.0 ) { 811 for ( G4int iString = 0; iString < 2; ++iS << 976 Quark_Mom[i].setX( SqrtScaleF * Quark_Mom[i].getX() ); 812 if ( iString == 0 ) { << 977 Quark_Mom[i].setY( SqrtScaleF * Quark_Mom[i].getY() ); 813 antiQuark = LeftAQ1; quark = LeftQ1; << 814 projectile->SetFirstParton( antiQuark << 815 projectile->SetSecondParton( quark ); << 816 projectile->SetStatus( 0 ); << 817 } else { // iString == 1 << 818 quark = LeftQ2; antiQuark = LeftAQ2; << 819 target->SetFirstParton( quark ); << 820 target->SetSecondParton( antiQuark ); << 821 target->SetStatus( 0 ); << 822 } << 823 G4int absAntiQuark = std::abs( antiQuark << 824 G4double aKsi = G4UniformRand(); << 825 if ( absAntiQuark == absQuark ) { << 826 if ( absAntiQuark != 3 ) { << 827 NewCode = 111; // Pi0-meson << 828 if ( aKsi < 0.5 ) { << 829 NewCode = 221; // Eta -meso << 830 if ( aKsi < 0.25 ) { << 831 NewCode = 331; // Eta'-meso << 832 } << 833 } << 834 } else { << 835 NewCode = 221; // Eta -meso << 836 if ( aKsi < 0.5 ) { << 837 NewCode = 331; // Eta'-meso << 838 } << 839 } << 840 } else { << 841 if ( absAntiQuark > absQuark ) { << 842 NewCode = absAntiQuark*100 + absQuar << 843 } else { << 844 NewCode = absQuark*100 + absAntiQuar << 845 } 978 } >> 979 //G4cout << "Anti Q " << i << " " << Quark_Mom[i] << G4endl; 846 } 980 } 847 TestParticle = G4ParticleTable::GetParti << 981 for ( G4int i = 2; i < 4; i++ ) { 848 if ( ! TestParticle ) return 99; // uns << 982 G4double Pz = -WminusTarget * Quark_Mom[i].getZ() / 2.0 + 849 if ( iString == 0 ) { << 983 ( ScaleFactor * ModMom2[i] + MassQ2 ) / 850 projectile->SetDefinition( TestParticl << 984 ( 2.0 * WminusTarget * Quark_Mom[i].getZ() ); 851 theParameters->SetProjMinDiffMass( 0.5 << 985 Quark_Mom[i].setZ( Pz ); 852 theParameters->SetProjMinNonDiffMass( << 986 if ( ScaleFactor != 1.0 ) { 853 } else { // iString == 1 << 987 Quark_Mom[i].setX( SqrtScaleF * Quark_Mom[i].getX() ); 854 target->SetDefinition( TestParticle ); << 988 Quark_Mom[i].setY( SqrtScaleF * Quark_Mom[i].getY() ); 855 theParameters->SetTarMinDiffMass( 0.5 << 856 theParameters->SetTarMinNonDiffMass( 0 << 857 } << 858 } // End of loop over the 2 string cases << 859 << 860 G4int QuarkOrder[2]; << 861 G4LorentzVector Pstring1, Pstring2; << 862 G4double Ystring1 = 0.0, Ystring2 = 0.0; << 863 << 864 for ( G4int iCase = 0; iCase < 2; ++iCase << 865 G4ThreeVector tmp = Quark_Mom[iCase] + Q << 866 G4LorentzVector Pstring( tmp, std::sqrt( << 867 std::sqrt( << 868 // Add protection for rapidity = 0.5*ln << 869 G4double Ystring = 0.0; << 870 if ( Pstring.e() > 1.0e-30 ) { << 871 if ( Pstring.e() + Pstring.pz() < 1.0e << 872 Ystring = -1.0e30; // A very large << 873 if ( Pstring.e() - Pstring.pz() < 1. << 874 Ystring = 1.0e30; // A very large << 875 } else { // Normal case << 876 Ystring = Pstring.rapidity(); << 877 } << 878 } 989 } >> 990 //G4cout << "Bary Q " << i << " " << Quark_Mom[i] << G4endl; 879 } 991 } 880 if ( iCase == 0 ) { // For the first st << 992 //G4cout << "Sum AQ " << Quark_Mom[0] + Quark_Mom[1] << G4endl 881 Pstring1 = Pstring; Ystring1 = Ystring << 993 // << "Sum Q " << Quark_Mom[2] + Quark_Mom[3] << G4endl; 882 } else { // For the second s << 883 Pstring2 = Pstring; Ystring2 = Ystring << 884 } << 885 } << 886 if ( Ystring1 > Ystring2 ) { << 887 common.Pprojectile = Pstring1; common.P << 888 QuarkOrder[0] = 0; QuarkOrder[1] = 1; << 889 } else { << 890 common.Pprojectile = Pstring2; common.P << 891 QuarkOrder[0] = 1; QuarkOrder[1] = 0; << 892 } << 893 994 894 if ( common.RotateStrings ) { << 995 G4ThreeVector tmp = Quark_Mom[0] + Quark_Mom[2]; 895 common.Pprojectile *= common.RandomRotat << 996 G4LorentzVector Pstring1( tmp, std::sqrt( Quark_Mom[0].mag2() + MassQ2 ) + 896 common.Ptarget *= common.RandomRotat << 997 std::sqrt( Quark_Mom[2].mag2() + MassQ2 ) ); 897 } << 998 G4double Ystring1 = Pstring1.rapidity(); >> 999 >> 1000 //G4cout << "Mom 1 string " << G4endl << Quark_Mom[0] << G4endl << Quark_Mom[2] << G4endl >> 1001 // << tmp << " " << tmp.mag() << G4endl; >> 1002 //G4cout << "1 str " << Pstring1 << " " << Pstring1.mag() << " " << Ystring1 << G4endl; >> 1003 >> 1004 tmp = Quark_Mom[1] + Quark_Mom[3]; >> 1005 G4LorentzVector Pstring2( tmp, std::sqrt( Quark_Mom[1].mag2() + MassQ2 ) + >> 1006 std::sqrt( Quark_Mom[3].mag2() + MassQ2 ) ); >> 1007 G4double Ystring2 = Pstring2.rapidity(); >> 1008 >> 1009 //G4cout << "Mom 2 string " << G4endl <<Quark_Mom[1] << G4endl << Quark_Mom[3] << G4endl >> 1010 // << tmp << " " << tmp.mag() << G4endl; >> 1011 //G4cout << "2 str " << Pstring2 << " " << Pstring2.mag() << " " << Ystring2 << G4endl; >> 1012 >> 1013 if ( Ystring1 > Ystring2 ) { >> 1014 Pprojectile = Pstring1; >> 1015 Ptarget = Pstring2; >> 1016 } else { >> 1017 Pprojectile = Pstring2; >> 1018 Ptarget = Pstring1; >> 1019 } 898 1020 899 common.Pprojectile.transform( common.toLab << 1021 //G4cout << "SumP CMS " << Pprojectile + Ptarget << " " << SqrtS << G4endl; 900 common.Ptarget.transform( common.toLab ); << 1022 Pprojectile.transform( toLab ); 901 << 1023 Ptarget.transform( toLab ); 902 G4LorentzVector Quark_4Mom[4]; << 1024 //G4cout << " SumP Lab " << Pprojectile + Ptarget << " " << SqrtS << G4endl; 903 for ( G4int i = 0; i < 4; ++i ) { << 1025 904 Quark_4Mom[i] = G4LorentzVector( Quark_M << 1026 // Calculation of the creation time 905 if ( common.RotateStrings ) Quark_4Mom[i << 1027 projectile->SetTimeOfCreation( target->GetTimeOfCreation() ); 906 Quark_4Mom[i].transform( common.toLab ); << 1028 projectile->SetPosition( target->GetPosition() ); 907 } << 1029 // Creation time and position of target nucleon were determined in >> 1030 // ReggeonCascade() of G4FTFModel >> 1031 //G4cout << "Mproj " << Pprojectile.mag() << G4endl << "Mtarg " << Ptarget.mag() << G4endl; >> 1032 projectile->Set4Momentum( Pprojectile ); >> 1033 target->Set4Momentum( Ptarget ); >> 1034 projectile->IncrementCollisionCount( 1 ); >> 1035 target->IncrementCollisionCount( 1 ); 908 1036 909 projectile->Splitting(); << 1037 theParameters->SetProbabilityOfAnnihilation( 0.0 ); // Uzhi March 2016 910 projectile->GetNextAntiParton()->Set4Momen << 911 projectile->GetNextParton()->Set4Momentum( << 912 << 913 target->Splitting(); << 914 target->GetNextParton()->Set4Momentum( Qua << 915 target->GetNextAntiParton()->Set4Momentum( << 916 1038 917 // Calculation of the creation time << 1039 return true; 918 // Creation time and position of target nu << 919 projectile->SetTimeOfCreation( target->Get << 920 projectile->SetPosition( target->GetPositi << 921 projectile->Set4Momentum( common.Pprojecti << 922 target->Set4Momentum( common.Ptarget ); << 923 1040 924 projectile->IncrementCollisionCount( 1 ); << 1041 } // End of if ( CandidatsN != 0 ) 925 target->IncrementCollisionCount( 1 ); << 926 1042 927 return 0; // Completed successfully: noth << 1043 } // End of if ( Ksi < ( X_a + X_b + X_c ) / Xannihilation ) 928 } // End of if ( CandidatsN != 0 ) << 929 1044 930 return 1; // Successfully ended, but the wo << 1045 // Simulation of anti-quark-quark string creation 931 } << 932 1046 >> 1047 if ( Ksi < ( X_a + X_b + X_c + X_d ) / Xannihilation ) { 933 1048 934 //-------------------------------------------- << 1049 #ifdef debugFTFannih >> 1050 G4cout << "Process d, only 1 quark - anti-quark string" << G4endl; >> 1051 #endif 935 1052 936 G4bool G4FTFAnnihilation:: << 1053 G4int CandidatsN( 0 ), CandAQ[36], CandQ[36]; 937 Create1QuarkAntiQuarkString( G4VSplitableHadro << 1054 G4int LeftAQ( 0 ), LeftQ( 0 ); 938 G4VSplitableHadro << 939 G4FTFParameters* << 940 G4FTFAnnihilation << 941 // Simulation of anti-quark - quark string c << 942 1055 943 #ifdef debugFTFannih << 1056 for ( G4int iAQ1 = 0; iAQ1 < 3; iAQ1++ ) { 944 G4cout << "Process d, only 1 quark - anti-qu << 1057 for ( G4int iAQ2 = 0; iAQ2 < 3; iAQ2++ ) { 945 #endif << 1058 if ( iAQ1 != iAQ2 ) { 946 << 1059 for ( G4int iQ1 = 0; iQ1 < 3; iQ1++ ) { 947 // Determine the set of candidates anti-quar << 1060 for ( G4int iQ2 = 0; iQ2 < 3; iQ2++ ) { 948 // Here "0", "1", "2" means, respectively, " << 1061 if ( iQ1 != iQ2 ) { 949 // of the (anti-baryon) projectile or (nucle << 1062 if ( -AQ[iAQ1] == Q[iQ1] && -AQ[iAQ2] == Q[iQ2] ) { 950 G4int CandidatsN = 0, CandAQ[36], CandQ[36]; << 1063 if ( iAQ1 == 0 && iAQ2 == 1 ) { CandAQ[CandidatsN] = 2; } 951 G4int LeftAQ = 0, LeftQ = 0; << 1064 if ( iAQ1 == 1 && iAQ2 == 0 ) { CandAQ[CandidatsN] = 2; } 952 for ( G4int iAQ1 = 0; iAQ1 < 3; ++iAQ1 ) { << 1065 953 for ( G4int iAQ2 = 0; iAQ2 < 3; ++iAQ2 ) { << 1066 if ( iAQ1 == 0 && iAQ2 == 2 ) { CandAQ[CandidatsN] = 1; } 954 if ( iAQ1 != iAQ2 ) { << 1067 if ( iAQ1 == 2 && iAQ2 == 0 ) { CandAQ[CandidatsN] = 1; } 955 for ( G4int iQ1 = 0; iQ1 < 3; ++iQ1 ) << 1068 956 for ( G4int iQ2 = 0; iQ2 < 3; ++iQ2 << 1069 if ( iAQ1 == 1 && iAQ2 == 2 ) { CandAQ[CandidatsN] = 0; } 957 if ( iQ1 != iQ2 ) { << 1070 if ( iAQ1 == 2 && iAQ2 == 1 ) { CandAQ[CandidatsN] = 0; } 958 if ( -common.AQ[iAQ1] == common. << 1071 959 if ( ( iAQ1 == 0 && i << 1072 if ( iQ1 == 0 && iQ2 == 1 ) { CandQ[CandidatsN] = 2; } 960 CandAQ[CandidatsN] = 2; << 1073 if ( iQ1 == 1 && iQ2 == 0 ) { CandQ[CandidatsN] = 2; } 961 } else if ( ( iAQ1 == 0 && i << 1074 962 CandAQ[CandidatsN] = 1; << 1075 if ( iQ1 == 0 && iQ2 == 2 ) { CandQ[CandidatsN] = 1; } 963 } else if ( ( iAQ1 == 1 && i << 1076 if ( iQ1 == 2 && iQ2 == 0 ) { CandQ[CandidatsN] = 1; } 964 CandAQ[CandidatsN] = 0; << 1077 >> 1078 if ( iQ1 == 1 && iQ2 == 2 ) { CandQ[CandidatsN] = 0; } >> 1079 if ( iQ1 == 2 && iQ2 == 1 ) { CandQ[CandidatsN] = 0; } >> 1080 CandidatsN++; 965 } 1081 } 966 if ( ( iQ1 == 0 && << 967 CandQ[CandidatsN] = 2; << 968 } else if ( ( iQ1 == 0 && << 969 CandQ[CandidatsN] = 1; << 970 } else if ( ( iQ1 == 1 && << 971 CandQ[CandidatsN] = 0; << 972 } << 973 ++CandidatsN; << 974 } 1082 } 975 } 1083 } 976 } 1084 } 977 } 1085 } 978 } 1086 } 979 } 1087 } 980 } << 981 << 982 if ( CandidatsN != 0 ) { << 983 G4int SampledCase = (G4int)G4RandFlat::sho << 984 LeftAQ = common.AQ[ CandAQ[SampledCase] ]; << 985 LeftQ = common.Q[ CandQ[SampledCase] ]; << 986 << 987 // Set the string properties << 988 projectile->SetFirstParton( LeftQ ); << 989 projectile->SetSecondParton( LeftAQ ); << 990 projectile->SetStatus( 0 ); << 991 G4int aAQ = std::abs( LeftAQ ), aQ = std:: << 992 G4int NewCode = 0; << 993 G4double aKsi = G4UniformRand(); << 994 // The string can have the quantum number << 995 // of the PDG code is, respectively, 1 and << 996 if ( aAQ == aQ ) { << 997 if ( aAQ != 3 ) { << 998 NewCode = 111; // Pi0-meson << 999 if ( aKsi < 0.5 ) { << 1000 NewCode = 221; // Eta -meson << 1001 if ( aKsi < 0.25 ) { << 1002 NewCode = 331; // Eta'-meson << 1003 } << 1004 } << 1005 } else { << 1006 NewCode = 221; // Eta -meson << 1007 if ( aKsi < 0.5 ) { << 1008 NewCode = 331; // Eta'-meson << 1009 } << 1010 } << 1011 } else { << 1012 if ( aAQ > aQ ) { << 1013 NewCode = aAQ*100 + aQ*10 + 1; NewCod << 1014 } else { << 1015 NewCode = aQ*100 + aAQ*10 + 1; NewCod << 1016 } << 1017 } << 1018 1088 1019 G4ParticleDefinition* TestParticle = G4Pa << 1089 if ( CandidatsN != 0 ) { 1020 if ( ! TestParticle ) return false; << 1090 G4int SampledCase = G4RandFlat::shootInt( G4long( CandidatsN ) ); 1021 projectile->SetDefinition( TestParticle ) << 1091 LeftAQ = AQ[ CandAQ[SampledCase] ]; 1022 theParameters->SetProjMinDiffMass( 0.5 ); << 1092 LeftQ = Q[ CandQ[SampledCase] ]; 1023 theParameters->SetProjMinNonDiffMass( 0.5 << 1093 //G4cout << "Left Aq Q " << LeftAQ << " " << LeftQ << G4endl; 1024 << 1094 1025 target->SetStatus( 4 ); // The target nu << 1095 // Set the string properties 1026 common.Pprojectile.setPx( 0.0 ); << 1096 projectile->SplitUp(); 1027 common.Pprojectile.setPy( 0.0 ); << 1097 //projectile->SetFirstParton( LeftAQ ); 1028 common.Pprojectile.setPz( 0.0 ); << 1098 //projectile->SetSecondParton( LeftQ ); 1029 common.Pprojectile.setE( common.SqrtS ); << 1099 projectile->SetFirstParton( LeftQ ); >> 1100 projectile->SetSecondParton( LeftAQ ); >> 1101 projectile->SetStatus( 0 ); 1030 1102 1031 common.Pprojectile.transform( common.toLa << 1103 // Uzhi March 2016 start >> 1104 G4int aAQ, aQ; >> 1105 aAQ=std::abs( LeftAQ ); aQ=std::abs( LeftQ ); >> 1106 >> 1107 G4int NewCode; >> 1108 G4double aKsi = G4UniformRand(); >> 1109 >> 1110 if ( aAQ == aQ ) >> 1111 { >> 1112 if ( aAQ != 3 ) >> 1113 { >> 1114 NewCode = 111; // Pi0-meson >> 1115 if ( aKsi < 0.5 ) >> 1116 { >> 1117 NewCode = 221; // Eta -meson >> 1118 if ( aKsi < 0.25 ) {NewCode = 331;} // Eta'-meson >> 1119 } >> 1120 } else >> 1121 { >> 1122 NewCode = 221; // Eta -meson >> 1123 if( aKsi < 0.5 ) {NewCode = 331;} // Eta'-meson >> 1124 } >> 1125 } else >> 1126 { >> 1127 if ( aAQ > aQ ){ NewCode = aAQ*100 + aQ*10 + 1; NewCode *= aAQ/LeftAQ; } >> 1128 else { NewCode = aQ*100 + aAQ*10 + 1; NewCode *= aQ/LeftQ; } >> 1129 } 1032 1130 1033 G4LorentzVector Pquark = G4LorentzVector << 1131 G4ParticleDefinition* TestParticle = G4ParticleTable::GetParticleTable()->FindParticle( NewCode ); 1034 G4LorentzVector Paquark = G4LorentzVector << 1132 if(!TestParticle) return false; >> 1133 projectile->SetDefinition( TestParticle ); >> 1134 theParameters->SetProjMinDiffMass( 0.5 ); // (0.5) // GeV Uzhi March 2016 >> 1135 theParameters->SetProjMinNonDiffMass( 0.5 ); >> 1136 // Uzhi March 2016 end >> 1137 >> 1138 target->SetStatus( 4 ); // The target nucleon has annihilated 3->4 Uzhi Oct 2014 >> 1139 Pprojectile.setPx( 0.0 ); >> 1140 Pprojectile.setPy( 0.0 ); >> 1141 Pprojectile.setPz( 0.0 ); >> 1142 Pprojectile.setE( SqrtS ); >> 1143 Pprojectile.transform( toLab ); >> 1144 >> 1145 // Calculation of the creation time >> 1146 projectile->SetTimeOfCreation( target->GetTimeOfCreation() ); >> 1147 projectile->SetPosition( target->GetPosition() ); >> 1148 // Creation time and position of target nucleon were determined in >> 1149 // ReggeonCascade() of G4FTFModel 1035 1150 1036 if ( common.RotateStrings ) { << 1151 //G4cout << "Mproj " << Pprojectile.mag() << G4endl << "Mtarg " << Ptarget.mag() << G4endl; 1037 Pquark *= common.RandomRotation; Paquar << 1152 projectile->Set4Momentum( Pprojectile ); 1038 } << 1039 Pquark.transform(common.toLab); projecti << 1040 Paquark.transform(common.toLab); projecti << 1041 1153 1042 projectile->Splitting(); << 1154 projectile->IncrementCollisionCount( 1 ); >> 1155 target->IncrementCollisionCount( 1 ); 1043 1156 1044 // Calculation of the creation time << 1157 theParameters->SetProbabilityOfAnnihilation( 0.0 ); // Uzhi March 2016 1045 // Creation time and position of target n << 1046 projectile->SetTimeOfCreation( target->Ge << 1047 projectile->SetPosition( target->GetPosit << 1048 projectile->Set4Momentum( common.Pproject << 1049 1158 1050 projectile->IncrementCollisionCount( 1 ); << 1159 return true; 1051 target->IncrementCollisionCount( 1 ); << 1160 } 1052 1161 1053 return true; << 1162 } // End of if ( Ksi < ( X_a + X_b + X_c + X_d ) / Xannihilation ) 1054 } // End of if ( CandidatsN != 0 ) << 1055 1163 >> 1164 //G4cout << "Pr Y " << Pprojectile.rapidity() << " Tr Y " << Ptarget.rapidity() << G4endl; 1056 return true; 1165 return true; 1057 } 1166 } 1058 1167 1059 1168 1060 //=========================================== 1169 //============================================================================ 1061 1170 1062 G4double G4FTFAnnihilation::ChooseX( G4double 1171 G4double G4FTFAnnihilation::ChooseX( G4double /* Alpha */, G4double /* Beta */ ) const { 1063 // If for sampling Xs other values of Alfa 1172 // If for sampling Xs other values of Alfa and Beta instead of 0.5 will be 1064 // chosen the method will be implemented 1173 // chosen the method will be implemented 1065 //G4double tmp = Alpha*Beta; 1174 //G4double tmp = Alpha*Beta; 1066 //tmp *= 1.0; 1175 //tmp *= 1.0; 1067 return 0.5; 1176 return 0.5; 1068 } 1177 } 1069 1178 1070 1179 >> 1180 1071 //=========================================== 1181 //============================================================================ 1072 1182 1073 G4ThreeVector G4FTFAnnihilation::GaussianPt( 1183 G4ThreeVector G4FTFAnnihilation::GaussianPt( G4double AveragePt2, G4double maxPtSquare ) const { 1074 // @@ this method is used in FTFModel as w 1184 // @@ this method is used in FTFModel as well. Should go somewhere common! 1075 G4double Pt2 = 0.0; << 1185 G4double Pt2( 0.0 ); 1076 if ( AveragePt2 <= 0.0 ) { 1186 if ( AveragePt2 <= 0.0 ) { 1077 Pt2 = 0.0; 1187 Pt2 = 0.0; 1078 } else { 1188 } else { 1079 Pt2 = -AveragePt2 * G4Log( 1.0 + G4Unifor 1189 Pt2 = -AveragePt2 * G4Log( 1.0 + G4UniformRand() * 1080 ( G4E 1190 ( G4Exp( -maxPtSquare/AveragePt2 ) -1.0 ) ); 1081 } 1191 } 1082 G4double Pt = std::sqrt( Pt2 ); 1192 G4double Pt = std::sqrt( Pt2 ); 1083 G4double phi = G4UniformRand() * twopi; 1193 G4double phi = G4UniformRand() * twopi; 1084 return G4ThreeVector ( Pt*std::cos( phi ), 1194 return G4ThreeVector ( Pt*std::cos( phi ), Pt*std::sin( phi ), 0.0 ); 1085 } 1195 } 1086 1196 1087 1197 1088 //=========================================== 1198 //============================================================================ 1089 1199 1090 void G4FTFAnnihilation::UnpackBaryon( G4int I 1200 void G4FTFAnnihilation::UnpackBaryon( G4int IdPDG, G4int& Q1, G4int& Q2, G4int& Q3 ) const { 1091 G4int AbsId = std::abs( IdPDG ); 1201 G4int AbsId = std::abs( IdPDG ); 1092 Q1 = AbsId / 1000; 1202 Q1 = AbsId / 1000; 1093 Q2 = ( AbsId % 1000 ) / 100; 1203 Q2 = ( AbsId % 1000 ) / 100; 1094 Q3 = ( AbsId % 100 ) / 10; 1204 Q3 = ( AbsId % 100 ) / 10; 1095 if ( IdPDG < 0 ) { Q1 = -Q1; Q2 = -Q2; Q3 = 1205 if ( IdPDG < 0 ) { Q1 = -Q1; Q2 = -Q2; Q3 = -Q3; } // Anti-baryon 1096 return; 1206 return; 1097 } 1207 } 1098 1208 1099 1209 1100 //=========================================== 1210 //============================================================================ 1101 1211 1102 G4FTFAnnihilation::G4FTFAnnihilation( const G 1212 G4FTFAnnihilation::G4FTFAnnihilation( const G4FTFAnnihilation& ) { 1103 throw G4HadronicException( __FILE__, __LINE 1213 throw G4HadronicException( __FILE__, __LINE__, 1104 "G4FTFAnnihilati << 1214 "G4FTFAnnihilation copy contructor not meant to be called" ); 1105 } 1215 } 1106 1216 1107 1217 1108 //=========================================== 1218 //============================================================================ 1109 1219 1110 const G4FTFAnnihilation & G4FTFAnnihilation:: 1220 const G4FTFAnnihilation & G4FTFAnnihilation::operator=( const G4FTFAnnihilation& ) { 1111 throw G4HadronicException( __FILE__, __LINE 1221 throw G4HadronicException( __FILE__, __LINE__, 1112 "G4FTFAnnihilati 1222 "G4FTFAnnihilation = operator not meant to be called" ); 1113 } 1223 } 1114 1224 1115 1225 1116 //=========================================== 1226 //============================================================================ 1117 1227 1118 G4bool G4FTFAnnihilation::operator==( const G << 1228 int G4FTFAnnihilation::operator==( const G4FTFAnnihilation& ) const { 1119 throw G4HadronicException( __FILE__, __LINE 1229 throw G4HadronicException( __FILE__, __LINE__, 1120 "G4FTFAnnihilati 1230 "G4FTFAnnihilation == operator not meant to be called" ); 1121 } 1231 } 1122 1232 1123 1233 1124 //=========================================== 1234 //============================================================================ 1125 1235 1126 G4bool G4FTFAnnihilation::operator!=( const G << 1236 int G4FTFAnnihilation::operator!=( const G4FTFAnnihilation& ) const { 1127 throw G4HadronicException( __FILE__, __LINE 1237 throw G4HadronicException( __FILE__, __LINE__, 1128 "G4DiffractiveEx 1238 "G4DiffractiveExcitation != operator not meant to be called" ); 1129 } 1239 } 1130 1240