[1199] | 1 | #include "G4BinaryLightIonReaction.hh" |
---|
| 2 | #include "G4Nucleus.hh" |
---|
| 3 | #include "G4Proton.hh" |
---|
| 4 | #include "G4PionPlus.hh" |
---|
| 5 | #include "G4ThreeVector.hh" |
---|
| 6 | #include "G4LorentzVector.hh" |
---|
| 7 | #include "G4DynamicParticle.hh" |
---|
| 8 | |
---|
| 9 | #include "G4LeptonConstructor.hh" |
---|
| 10 | #include "G4IonConstructor.hh" |
---|
| 11 | #include "G4Gamma.hh" |
---|
| 12 | |
---|
| 13 | G4V3DNucleus * global3DNucleus; |
---|
| 14 | |
---|
| 15 | |
---|
| 16 | void Init(G4int code, vector<G4ParticleDefinition *> & particles); |
---|
| 17 | void GetParam(G4int & A, G4int & Z, G4double & pMin, G4double & pMax, |
---|
| 18 | G4int &nEvents, G4int & particleCode); |
---|
| 19 | G4ThreeVector GetSpherePoint(G4double r); |
---|
| 20 | G4ThreeVector GetSurfacePoint(G4double r); |
---|
| 21 | |
---|
| 22 | int main(int argc, char * argv[]) |
---|
| 23 | { |
---|
| 24 | G4String a=G4Proton::ProtonDefinition()->GetParticleName(); |
---|
| 25 | G4cout << " At Start " << a << G4endl; |
---|
| 26 | G4cout.setf( std::ios::scientific, std::ios::floatfield ); |
---|
| 27 | // totalRes.Start(); |
---|
| 28 | G4LeptonConstructor leptons; |
---|
| 29 | leptons.ConstructParticle(); |
---|
| 30 | G4IonConstructor ions; |
---|
| 31 | ions.ConstructParticle(); |
---|
| 32 | |
---|
| 33 | G4int A, Z, nEvents, particleCode; |
---|
| 34 | G4double pMin, pMax; |
---|
| 35 | vector<G4ParticleDefinition *> particles; |
---|
| 36 | if(argc == 7) |
---|
| 37 | { |
---|
| 38 | A = atoi(argv[1]); |
---|
| 39 | Z = atoi(argv[2]); |
---|
| 40 | pMin = G4double(atoi(argv[3]))*MeV; |
---|
| 41 | pMax = G4double(atoi(argv[4]))*MeV; |
---|
| 42 | nEvents = atoi(argv[5]); |
---|
| 43 | particleCode = atoi(argv[6]); |
---|
| 44 | } |
---|
| 45 | else |
---|
| 46 | GetParam(A, Z, pMin, pMax, nEvents, particleCode); |
---|
| 47 | |
---|
| 48 | if(pMin > pMax) |
---|
| 49 | { |
---|
| 50 | G4double pTmp = pMin; |
---|
| 51 | pMin = pMax; |
---|
| 52 | pMax = pTmp; |
---|
| 53 | } |
---|
| 54 | |
---|
| 55 | Init(particleCode, particles); |
---|
| 56 | |
---|
| 57 | G4BinaryCascade hkm; |
---|
| 58 | hkm.SetDeExcitation(0); |
---|
| 59 | |
---|
| 60 | for(G4int event = 0; event < nEvents; ++event) |
---|
| 61 | { |
---|
| 62 | G4cerr << " ---- event " << event << "----------------------"<< G4endl; |
---|
| 63 | G4V3DNucleus * fancyNucleus = new G4Fancy3DNucleus; // for Propagate() |
---|
| 64 | fancyNucleus->Init(16, 8); |
---|
| 65 | G4V3DNucleus * projectile = new G4Fancy3DNucleus; |
---|
| 66 | projectile->Init(12, 6); |
---|
| 67 | G4double radius = fancyNucleus->GetOuterRadius(); |
---|
| 68 | G4double impactMax = fancyNucleus->GetOuterRadius()+projectile->GetOuterRadius(); |
---|
| 69 | // position |
---|
| 70 | G4double aX=(2.*G4UniformRand()-1.)*impactMax; |
---|
| 71 | G4double aY=(2.*G4UniformRand()-1.)*impactMax; |
---|
| 72 | |
---|
| 73 | G4ThreeVector pos(aX, aY, -20.*fermi); |
---|
| 74 | // momentum |
---|
| 75 | G4double mass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(6, 12); |
---|
| 76 | G4double p = 280*MeV; |
---|
| 77 | G4double e = sqrt(mass*mass+p*p); |
---|
| 78 | G4LorentzVector mom(0, p, 0, e); |
---|
| 79 | |
---|
| 80 | // build the projectile |
---|
| 81 | G4KineticTrackVector * ktv = new G4KineticTrackVector; |
---|
| 82 | G4ThreeVector beta = mom.beta(); |
---|
| 83 | projectile->DoLorentzBoost(beta); |
---|
| 84 | projectile->StartLoop(); |
---|
| 85 | G4Nucleon * aNuc; |
---|
| 86 | while( (aNuc=projectile->GetNextNucleon()) ) |
---|
| 87 | { |
---|
| 88 | G4KineticTrack * it = |
---|
| 89 | new G4KineticTrack(aNuc->GetDefinition(), 0, |
---|
| 90 | aNuc->GetPosition()+pos, const_cast<G4LorentzVector &>(aNuc->GetMomentum()) ); |
---|
| 91 | ktv->push_back(it); |
---|
| 92 | } |
---|
| 93 | |
---|
| 94 | // do the reactions |
---|
| 95 | G4ReactionProductVector *result=hkm.Propagate(ktv, fancyNucleus); |
---|
| 96 | |
---|
| 97 | delete fancyNucleus; |
---|
| 98 | delete projectile; |
---|
| 99 | |
---|
| 100 | G4ReactionProductVector::iterator result_iter; |
---|
| 101 | G4int isec=0; |
---|
| 102 | for(result_iter = result->begin();result_iter != result->end(); |
---|
| 103 | ++result_iter) |
---|
| 104 | { |
---|
| 105 | G4ThreeVector pFinal= (*result_iter)->GetMomentum(); |
---|
| 106 | |
---|
| 107 | G4cout << " Final particle # "<< ++isec << ": " |
---|
| 108 | << pFinal.x() << " " |
---|
| 109 | << pFinal.y() << " " |
---|
| 110 | << pFinal.z() << " " |
---|
| 111 | << (*result_iter)->GetTotalEnergy() << " " |
---|
| 112 | << G4endl; |
---|
| 113 | } |
---|
| 114 | } |
---|
| 115 | cout.flush(); |
---|
| 116 | return 0; |
---|
| 117 | } |
---|
| 118 | |
---|
| 119 | /* |
---|
| 120 | |
---|
| 121 | void PrintResourceUsage(ResourceMgr & resMgr) |
---|
| 122 | { |
---|
| 123 | double userTime = resMgr.GetUserTime(); |
---|
| 124 | double sysTime = resMgr.GetSysTime(); |
---|
| 125 | double cpuTime = userTime+sysTime; |
---|
| 126 | cout << "cpu time: " << cpuTime << " (" << userTime << "+" |
---|
| 127 | << sysTime << ")" << endl; |
---|
| 128 | } |
---|
| 129 | |
---|
| 130 | |
---|
| 131 | */ |
---|
| 132 | |
---|
| 133 | void Init(G4int code, vector<G4ParticleDefinition *> & particles) |
---|
| 134 | { |
---|
| 135 | if(code == 0) |
---|
| 136 | { |
---|
| 137 | particles.push_back(G4Gamma::Gamma()); |
---|
| 138 | return; |
---|
| 139 | } |
---|
| 140 | |
---|
| 141 | G4int count = 0; |
---|
| 142 | while(code != 0) |
---|
| 143 | { |
---|
| 144 | if(code & 1) |
---|
| 145 | { |
---|
| 146 | switch (count) |
---|
| 147 | { |
---|
| 148 | case 0: |
---|
| 149 | particles.push_back(G4Proton::Proton()); |
---|
| 150 | break; |
---|
| 151 | case 1: |
---|
| 152 | particles.push_back(G4Neutron::Neutron()); |
---|
| 153 | break; |
---|
| 154 | case 2: |
---|
| 155 | particles.push_back(G4AntiProton::AntiProton()); |
---|
| 156 | break; |
---|
| 157 | case 3: |
---|
| 158 | particles.push_back(G4PionPlus::PionPlus()); |
---|
| 159 | break; |
---|
| 160 | case 4: |
---|
| 161 | particles.push_back(G4PionZero::PionZero()); |
---|
| 162 | break; |
---|
| 163 | case 5: |
---|
| 164 | particles.push_back(G4PionMinus::PionMinus()); |
---|
| 165 | break; |
---|
| 166 | case 6: |
---|
| 167 | particles.push_back(G4KaonPlus::KaonPlus()); |
---|
| 168 | break; |
---|
| 169 | case 7: |
---|
| 170 | particles.push_back(G4KaonZero::KaonZero()); |
---|
| 171 | break; |
---|
| 172 | case 8: |
---|
| 173 | particles.push_back(G4KaonMinus::KaonMinus()); |
---|
| 174 | break; |
---|
| 175 | case 9: |
---|
| 176 | particles.push_back(G4SigmaPlus::SigmaPlus()); |
---|
| 177 | break; |
---|
| 178 | case 10: |
---|
| 179 | particles.push_back(G4SigmaZero::SigmaZero()); |
---|
| 180 | break; |
---|
| 181 | case 11: |
---|
| 182 | particles.push_back(G4SigmaMinus::SigmaMinus()); |
---|
| 183 | break; |
---|
| 184 | } |
---|
| 185 | } |
---|
| 186 | ++count; |
---|
| 187 | code = code>>1; |
---|
| 188 | } |
---|
| 189 | } |
---|
| 190 | |
---|
| 191 | |
---|
| 192 | void GetParam(G4int & A, G4int & Z, G4double & pMin, G4double & pMax, |
---|
| 193 | G4int &nEvents, G4int & particleCode) |
---|
| 194 | { |
---|
| 195 | // A and Z |
---|
| 196 | cout << "Input Nucleus A and Z: "; |
---|
| 197 | cin >> A >> Z; |
---|
| 198 | // number of events |
---|
| 199 | cout << "Input number of events: "; |
---|
| 200 | cin >> nEvents; |
---|
| 201 | } |
---|
| 202 | |
---|
| 203 | |
---|
| 204 | |
---|
| 205 | |
---|
| 206 | G4ThreeVector GetSpherePoint(G4double r) |
---|
| 207 | { |
---|
| 208 | // random point INSIDE a sphere |
---|
| 209 | |
---|
| 210 | G4double x = r*(2*G4UniformRand()-1); |
---|
| 211 | G4double tmp = sqrt(r*r-x*x); |
---|
| 212 | G4double y = tmp*(2*G4UniformRand()-1); |
---|
| 213 | tmp = sqrt(r*r-x*x-y*y); |
---|
| 214 | G4double z = tmp*(2*G4UniformRand()-1); |
---|
| 215 | G4ThreeVector point; |
---|
| 216 | point.setX(x); |
---|
| 217 | point.setY(y); |
---|
| 218 | point.setZ(z); |
---|
| 219 | return point; |
---|
| 220 | } |
---|
| 221 | G4ThreeVector GetSurfacePoint(G4double r) |
---|
| 222 | { |
---|
| 223 | // random point oin surface of a sphere |
---|
| 224 | |
---|
| 225 | G4double x = r*(2*G4UniformRand()-1); |
---|
| 226 | G4double tmp = sqrt(r*r-x*x); |
---|
| 227 | G4double y = tmp*(2*G4UniformRand()-1); |
---|
| 228 | G4double z = -sqrt(r*r-x*x-y*y); |
---|
| 229 | G4ThreeVector point; |
---|
| 230 | point.setX(x); |
---|
| 231 | point.setY(y); |
---|
| 232 | point.setZ(z); |
---|
| 233 | return point; |
---|
| 234 | } |
---|