Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-11-22 02:45:26 

0001 #include "G4ProcessManager.hh"
0002 #include "G4ParticleTable.hh"
0003 #include "G4HadronicException.hh"
0004 
0005 #include "SimG4Core/CustomPhysics/interface/FullModelHadronicProcess.h"
0006 #include "SimG4Core/CustomPhysics/interface/G4ProcessHelper.h"
0007 #include "SimG4Core/CustomPhysics/interface/Decay3Body.h"
0008 #include "SimG4Core/CustomPhysics/interface/CustomPDGParser.h"
0009 #include "SimG4Core/CustomPhysics/interface/CustomParticle.h"
0010 
0011 using namespace CLHEP;
0012 
0013 FullModelHadronicProcess::FullModelHadronicProcess(G4ProcessHelper* aHelper, const G4String& processName)
0014     : G4VDiscreteProcess(processName), theHelper(aHelper) {}
0015 
0016 FullModelHadronicProcess::~FullModelHadronicProcess() {}
0017 
0018 G4bool FullModelHadronicProcess::IsApplicable(const G4ParticleDefinition& aP) {
0019   return theHelper->ApplicabilityTester(aP);
0020 }
0021 
0022 G4double FullModelHadronicProcess::GetMicroscopicCrossSection(const G4DynamicParticle* aParticle,
0023                                                               const G4Element* anElement,
0024                                                               G4double aTemp) {
0025   //Get the cross section for this particle/element combination from the ProcessHelper
0026   G4double InclXsec = theHelper->GetInclusiveCrossSection(aParticle, anElement);
0027   //  G4cout<<"Returned cross section from helper was: "<<InclXsec/millibarn<<" millibarn"<<G4endl;
0028   return InclXsec;
0029 }
0030 
0031 G4double FullModelHadronicProcess::GetMeanFreePath(const G4Track& aTrack, G4double, G4ForceCondition*) {
0032   G4Material* aMaterial = aTrack.GetMaterial();
0033   const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
0034   G4double sigma = 0.0;
0035 
0036   G4int nElements = aMaterial->GetNumberOfElements();
0037 
0038   const G4double* theAtomicNumDensityVector = aMaterial->GetAtomicNumDensityVector();
0039   G4double aTemp = aMaterial->GetTemperature();
0040 
0041   for (G4int i = 0; i < nElements; ++i) {
0042     G4double xSection = GetMicroscopicCrossSection(aParticle, (*aMaterial->GetElementVector())[i], aTemp);
0043     sigma += theAtomicNumDensityVector[i] * xSection;
0044   }
0045   G4double res = DBL_MAX;
0046   if (sigma > 0.0) {
0047     res = 1. / sigma;
0048   }
0049   return res;
0050 }
0051 
0052 G4VParticleChange* FullModelHadronicProcess::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) {
0053   //  G4cout<<"**** Entering FullModelHadronicProcess::PostStepDoIt       ******"<<G4endl;
0054   const G4TouchableHandle& thisTouchable(aTrack.GetTouchableHandle());
0055 
0056   // A little setting up
0057   aParticleChange.Initialize(aTrack);
0058   const G4DynamicParticle* IncidentRhadron = aTrack.GetDynamicParticle();
0059   CustomParticle* CustomIncident = static_cast<CustomParticle*>(IncidentRhadron->GetDefinition());
0060   const G4ThreeVector& aPosition = aTrack.GetPosition();
0061   //  G4cout<<"G: "<<aStep.GetStepLength()/cm<<G4endl;
0062   const G4int theIncidentPDG = IncidentRhadron->GetDefinition()->GetPDGEncoding();
0063   G4ParticleTable* theParticleTable = G4ParticleTable::GetParticleTable();
0064   std::vector<G4ParticleDefinition*> theParticleDefinitions;
0065   G4bool IncidentSurvives = false;
0066   G4bool TargetSurvives = false;
0067   G4Nucleus targetNucleus(aTrack.GetMaterial());
0068   G4ParticleDefinition* outgoingRhadron = nullptr;
0069   G4ParticleDefinition* outgoingCloud = nullptr;
0070   G4ParticleDefinition* outgoingTarget = nullptr;
0071 
0072   G4ThreeVector p_0 = IncidentRhadron->GetMomentum();
0073   G4double e_kin_0 = IncidentRhadron->GetKineticEnergy();
0074   //  G4cout<<e_kin_0/GeV<<G4endl;
0075 
0076   G4DynamicParticle* cloudParticle = new G4DynamicParticle();
0077   /*
0078   if(CustomPDGParser::s_isRMeson(theIncidentPDG))
0079     G4cout<<"Rmeson"<<G4endl;
0080   if(CustomPDGParser::s_isRBaryon(theIncidentPDG)) 
0081     G4cout<<"Rbaryon"<<G4endl;
0082   */
0083   cloudParticle->SetDefinition(CustomIncident->GetCloud());
0084 
0085   if (cloudParticle->GetDefinition() == nullptr) {
0086     G4cout << "FullModelHadronicProcess::PostStepDoIt  Definition of particle cloud not available!!" << G4endl;
0087   }
0088   /*
0089   G4cout<<"Incoming particle was "<<IncidentRhadron->GetDefinition()->GetParticleName()
0090   <<". Corresponding cloud is "<<cloudParticle->GetDefinition()->GetParticleName()<<G4endl;
0091   G4cout<<"Kinetic energy was: "<<IncidentRhadron->GetKineticEnergy()/GeV<<" GeV"<<G4endl;
0092   */
0093   double scale = cloudParticle->GetDefinition()->GetPDGMass() / IncidentRhadron->GetDefinition()->GetPDGMass();
0094   //  G4cout<<"Mass ratio: "<<scale<<G4endl;
0095   G4LorentzVector cloudMomentum(IncidentRhadron->GetMomentum() * scale, cloudParticle->GetDefinition()->GetPDGMass());
0096   G4LorentzVector gluinoMomentum(IncidentRhadron->GetMomentum() * (1. - scale),
0097                                  CustomIncident->GetSpectator()->GetPDGMass());
0098 
0099   //These two for getting CMS transforms later (histogramming purposes...)
0100   G4LorentzVector FullRhadron4Momentum = IncidentRhadron->Get4Momentum();
0101   const G4LorentzVector& Cloud4Momentum = cloudMomentum;
0102 
0103   cloudParticle->Set4Momentum(cloudMomentum);
0104 
0105   G4DynamicParticle* OrgPart = cloudParticle;
0106 
0107   /*
0108   G4cout<<"Original momentum: "<<IncidentRhadron->Get4Momentum().v().mag()/GeV
0109   <<" GeV, corresponding to gamma: "
0110   <<IncidentRhadron->GetTotalEnergy()/IncidentRhadron->GetDefinition()->GetPDGMass()<<G4endl;
0111   
0112   G4cout<<"Cloud momentum: "<<cloudParticle->Get4Momentum().v().mag()/GeV
0113   <<" GeV, corresponding to gamma: "
0114   <<cloudParticle->GetTotalEnergy()/cloudParticle->GetDefinition()->GetPDGMass()<<G4endl;
0115   */
0116 
0117   double E_0 = IncidentRhadron->GetKineticEnergy();
0118   G4double ek = OrgPart->GetKineticEnergy();
0119   G4double amas = OrgPart->GetDefinition()->GetPDGMass();
0120   G4ThreeVector dir = (OrgPart->GetMomentum()).unit();
0121   G4double tkin = targetNucleus.Cinema(ek);
0122   ek += tkin;
0123 
0124   // calculate black track energies
0125   tkin = targetNucleus.EvaporationEffects(ek);
0126   ek -= tkin;
0127 
0128   if (ek + gluinoMomentum.e() - gluinoMomentum.m() <= 0.1 * MeV || ek <= 0.) {
0129     //Very rare event...
0130     G4cout << "Kinetic energy is sick" << G4endl;
0131     G4cout << "Full R-hadron: " << (ek + gluinoMomentum.e() - gluinoMomentum.m()) / MeV << " MeV" << G4endl;
0132     G4cout << "Quark system: " << ek / MeV << " MeV" << G4endl;
0133     aParticleChange.ProposeTrackStatus(fStopButAlive);  // AR_NEWCODE_IMPORT
0134     return &aParticleChange;
0135   }
0136   OrgPart->SetKineticEnergy(ek);
0137   G4double p = std::sqrt(ek * (ek + 2 * amas));
0138   OrgPart->SetMomentum(dir * p);
0139 
0140   //Get the final state particles
0141   G4ParticleDefinition* aTarget;
0142   ReactionProduct rp = theHelper->GetFinalState(aTrack, aTarget);
0143   G4bool force2to2 = false;
0144   //  G4cout<<"Trying to get final state..."<<G4endl;
0145   while (rp.size() != 2 && force2to2) {
0146     rp = theHelper->GetFinalState(aTrack, aTarget);
0147   }
0148   G4int NumberOfSecondaries = rp.size();
0149   //  G4cout<<"Multiplicity of selected final state: "<<rp.size()<<G4endl;
0150 
0151   //Getting CMS transforms. Boosting is done at histogram filling
0152   G4LorentzVector Target4Momentum(0., 0., 0., aTarget->GetPDGMass());
0153 
0154   G4LorentzVector psum_full = FullRhadron4Momentum + Target4Momentum;
0155   G4LorentzVector psum_cloud = Cloud4Momentum + Target4Momentum;
0156   G4ThreeVector trafo_full_cms = (-1) * psum_full.boostVector();
0157   G4ThreeVector trafo_cloud_cms = (-1) * psum_cloud.boostVector();
0158 
0159   // OK Let's make some particles :-)
0160   // We're not using them for anything yet, I know, but let's make sure the machinery is there
0161   for (ReactionProduct::iterator it = rp.begin(); it != rp.end(); ++it) {
0162     G4ParticleDefinition* tempDef = theParticleTable->FindParticle(*it);
0163     CustomParticle* tempCust = dynamic_cast<CustomParticle*>(tempDef);
0164     if (tempDef == aTarget)
0165       TargetSurvives = true;
0166 
0167     //      if (tempDef->GetParticleType()=="rhadron")
0168     if (tempCust != nullptr) {
0169       outgoingRhadron = tempDef;
0170       //Setting outgoing cloud definition
0171       outgoingCloud = tempCust->GetCloud();
0172       if (outgoingCloud == nullptr) {
0173         G4cout << "FullModelHadronicProcess::PostStepDoIt  Definition of outgoing particle cloud not available!"
0174                << G4endl;
0175       }
0176       /*
0177     G4cout<<"Outgoing Rhadron is: "<<outgoingRhadron->GetParticleName()<<G4endl;
0178     G4cout<<"Outgoing cloud is: "<<outgoingCloud->GetParticleName()<<G4endl;
0179       */
0180     }
0181 
0182     if (tempDef == G4Proton::Proton() || tempDef == G4Neutron::Neutron())
0183       outgoingTarget = tempDef;
0184     if (tempCust == nullptr && rp.size() == 2)
0185       outgoingTarget = tempDef;
0186     if (tempDef->GetPDGEncoding() == theIncidentPDG) {
0187       IncidentSurvives = true;
0188     } else {
0189       theParticleDefinitions.push_back(tempDef);
0190     }
0191   }
0192 
0193   if (outgoingTarget == nullptr)
0194     outgoingTarget = theParticleTable->FindParticle(rp[1]);
0195 
0196   // A little debug information
0197   /*
0198   G4cout<<"The particles coming out of this reaction will be: ";
0199   for (std::vector<G4DynamicParticle*>::iterator it = theDynamicParticles.begin();
0200        it != theDynamicParticles.end();
0201        it++){
0202     G4cout<< (*it)->GetDefinition()->GetParticleName()<<" ";
0203   }
0204   G4cout<<G4endl;
0205   */
0206   // If the incident particle survives it is not a "secondary", although
0207   // it would probably be easier to fStopAndKill the track every time.
0208   if (IncidentSurvives)
0209     NumberOfSecondaries--;
0210   aParticleChange.SetNumberOfSecondaries(NumberOfSecondaries);
0211 
0212   // ADAPTED FROM G4LEPionMinusInelastic::ApplyYourself
0213   // It is bloody ugly, but such is the way of cut 'n' paste
0214 
0215   // Set up the incident
0216   // This is where rotation to z-axis is done
0217   const G4HadProjectile* originalIncident = new G4HadProjectile(*OrgPart);
0218 
0219   //Maybe I need to calculate trafo from R-hadron... Bollocks! Labframe does not move. CMS does.
0220   G4HadProjectile* org2 = new G4HadProjectile(*OrgPart);
0221   G4LorentzRotation trans = org2->GetTrafoToLab();
0222   delete org2;
0223 
0224   // create the target particle
0225 
0226   G4DynamicParticle* originalTarget = new G4DynamicParticle;
0227   originalTarget->SetDefinition(aTarget);
0228 
0229   G4ReactionProduct targetParticle(aTarget);
0230 
0231   G4ReactionProduct currentParticle(const_cast<G4ParticleDefinition*>(originalIncident->GetDefinition()));
0232   currentParticle.SetMomentum(originalIncident->Get4Momentum().vect());
0233   currentParticle.SetKineticEnergy(originalIncident->GetKineticEnergy());
0234 
0235   /*
0236   G4cout<<"After creation:"<<G4endl;
0237   G4cout<<"currentParticle: "<<currentParticle.GetMomentum()/GeV<<" GeV vs. "<<OrgPart->Get4Momentum()/GeV<<" GeV"<<G4endl;
0238   G4cout<<"targetParticle: "<<targetParticle.GetMomentum()/GeV<<" GeV"<<G4endl;
0239   G4cout<<"Fourmomentum from originalIncident: "<<originalIncident->Get4Momentum()<<" vs "<<OrgPart->Get4Momentum()<<G4endl;
0240   */
0241 
0242   G4ReactionProduct modifiedOriginal = currentParticle;
0243 
0244   currentParticle.SetSide(1);  // incident always goes in forward hemisphere
0245   targetParticle.SetSide(-1);  // target always goes in backward hemisphere
0246   G4bool incidentHasChanged = false;
0247   if (!IncidentSurvives)
0248     incidentHasChanged = true;  //I wonder if I am supposed to do this...
0249   G4bool targetHasChanged = false;
0250   if (!TargetSurvives)
0251     targetHasChanged = true;  //Ditto here
0252   G4bool quasiElastic = false;
0253   if (rp.size() == 2)
0254     quasiElastic = true;                                //Oh well...
0255   G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> vec;  // vec will contain the secondary particles
0256   G4int vecLen = 0;
0257   vec.Initialize(0);
0258 
0259   // I hope my understanding of "secondary" is correct here
0260   // I think that it entails only what is not a surviving incident of target
0261 
0262   for (G4int i = 0; i != NumberOfSecondaries; i++) {
0263     if (theParticleDefinitions[i] != aTarget && theParticleDefinitions[i] != originalIncident->GetDefinition() &&
0264         theParticleDefinitions[i] != outgoingRhadron && theParticleDefinitions[i] != outgoingTarget) {
0265       G4ReactionProduct* pa = new G4ReactionProduct;
0266       pa->SetDefinition(theParticleDefinitions[i]);
0267       (G4UniformRand() < 0.5) ? pa->SetSide(-1) : pa->SetSide(1);
0268       vec.SetElement(vecLen++, pa);
0269     }
0270   }
0271 
0272   if (incidentHasChanged)
0273     currentParticle.SetDefinitionAndUpdateE(outgoingCloud);
0274   if (incidentHasChanged)
0275     modifiedOriginal.SetDefinition(
0276         outgoingCloud);  //Is this correct? It solves the "free energy" checking in ReactionDynamics
0277   if (targetHasChanged)
0278     targetParticle.SetDefinitionAndUpdateE(outgoingTarget);
0279 
0280   //  G4cout<<"Calling CalculateMomenta... "<<G4endl;
0281   /*
0282   G4cout<<"Current particle starts as: "<<currentParticle.GetDefinition()->GetParticleName()<<G4endl;
0283   G4cout<<"with momentum: "<<currentParticle.GetMomentum()/GeV<<" GeV"<<G4endl;
0284   G4cout<<"May be killed?: "<<currentParticle.GetMayBeKilled()<<G4endl;
0285   */
0286 
0287   CalculateMomenta(vec,
0288                    vecLen,
0289                    originalIncident,
0290                    originalTarget,
0291                    modifiedOriginal,
0292                    targetNucleus,
0293                    currentParticle,
0294                    targetParticle,
0295                    incidentHasChanged,
0296                    targetHasChanged,
0297                    quasiElastic);
0298 
0299   //  G4cout <<"Cloud loss: "<<(e_temp-currentParticle.GetKineticEnergy())/GeV<<" GeV"<<G4endl;
0300 
0301   G4String cPname = currentParticle.GetDefinition()->GetParticleName();
0302 
0303   //  if(cPname!="rhadronmesoncloud"&&cPname!="rhadronbaryoncloud")
0304   //    {
0305   /*
0306   G4cout<<"Current particle is now: "<<cPname <<G4endl;
0307   G4cout<<"with momentum: "<<currentParticle.GetMomentum()/GeV<<" GeV"<<G4endl;
0308   G4cout<<"and kinetic energy: "<<currentParticle.GetKineticEnergy()/GeV<<" GeV"<<G4endl;
0309   G4cout<<"May be killed?: "<<currentParticle.GetMayBeKilled()<<G4endl;
0310   G4cout<<"Modified original is: "<<modifiedOriginal.GetDefinition()->GetParticleName()<<G4endl;
0311   G4cout<<"with momentum: "<<modifiedOriginal.GetMomentum()/GeV<<" GeV"<<G4endl;
0312   G4cout<<"and kinetic energy: "<<modifiedOriginal.GetKineticEnergy()/GeV<<" GeV"<<G4endl;
0313   G4cout<<"May be killed?: "<<modifiedOriginal.GetMayBeKilled()<<G4endl;
0314   G4cout<<"Target particle is: "<<targetParticle.GetDefinition()->GetParticleName()<<G4endl;
0315   G4cout<<"with momentum: "<<targetParticle.GetMomentum()/GeV<<" GeV"<<G4endl;
0316   G4cout<<"and kinetic energy: "<<targetParticle.GetKineticEnergy()/GeV<<" GeV"<<G4endl;
0317   G4cout<<"May be killed?: "<<targetParticle.GetMayBeKilled()<<G4endl;
0318   G4cout<<"incidentHasChanged: "<<incidentHasChanged<<G4endl;
0319   G4cout<<"targetHasChanged: "<<targetHasChanged<<G4endl;
0320   G4cout<<"Particles in vec:"<<G4endl;
0321   for(int i=0; i<vecLen; ++i )
0322     {
0323       G4cout<< vec[i]->GetDefinition()->GetParticleName()<<G4endl;
0324     }
0325   */
0326   // G4cout<<"Done!"<<G4endl;
0327 
0328   aParticleChange.SetNumberOfSecondaries(vecLen + NumberOfSecondaries);
0329   G4double e_kin = 0;
0330   G4LorentzVector cloud_p4_new;  //Cloud 4-momentum in lab after collision
0331   //  n++;
0332   //  G4cout << n << G4endl;
0333   /*
0334   if(cPname!="rhadronmesoncloud"&&cPname!="rhadronbaryoncloud") {
0335     G4cout<<"Cloud deleted!!! AAARRRRGGGHHH!!!"<<G4endl;
0336     G4cout<<"Cloud name: "<<cPname<<G4endl;
0337     G4cout<<"E_kin_0: "<<e_kin_0/GeV<<" GeV"<<G4endl;
0338     //    G4cout<<"n: "<<n<<G4endl;
0339     //    n=0;
0340   }
0341   */
0342   cloud_p4_new.setVectM(currentParticle.GetMomentum(), currentParticle.GetMass());
0343   cloud_p4_new *= trans;
0344 
0345   G4LorentzVector cloud_p4_old_full = Cloud4Momentum;  //quark system in CMS BEFORE collision
0346   cloud_p4_old_full.boost(trafo_full_cms);
0347   G4LorentzVector cloud_p4_old_cloud = Cloud4Momentum;  //quark system in cloud CMS BEFORE collision
0348   cloud_p4_old_cloud.boost(trafo_cloud_cms);
0349   G4LorentzVector cloud_p4_full = cloud_p4_new;  //quark system in CMS AFTER collision
0350   cloud_p4_full.boost(trafo_full_cms);
0351   G4LorentzVector cloud_p4_cloud = cloud_p4_new;  //quark system in cloud CMS AFTER collision
0352   cloud_p4_cloud.boost(trafo_cloud_cms);
0353 
0354   G4LorentzVector p_g_cms = gluinoMomentum;  //gluino in CMS BEFORE collision
0355   p_g_cms.boost(trafo_full_cms);
0356 
0357   G4LorentzVector p4_new(cloud_p4_new.v() + gluinoMomentum.v(), outgoingRhadron->GetPDGMass());
0358   //  G4cout<<"P4-diff: "<<(p4_new-cloud_p4_new-gluinoMomentum)/GeV<<", magnitude: "
0359   // <<(p4_new-cloud_p4_new-gluinoMomentum).m()/MeV<<" MeV" <<G4endl;
0360 
0361   G4ThreeVector p_new = p4_new.vect();
0362 
0363   aParticleChange.ProposeLocalEnergyDeposit((p4_new - cloud_p4_new - gluinoMomentum).m());
0364 
0365   if (incidentHasChanged) {
0366     G4DynamicParticle* p0 = new G4DynamicParticle;
0367     p0->SetDefinition(outgoingRhadron);
0368     p0->SetMomentum(p_new);
0369 
0370     // May need to run SetDefinitionAndUpdateE here...
0371     G4Track* Track0 = new G4Track(p0, aTrack.GetGlobalTime(), aPosition);
0372     Track0->SetTouchableHandle(thisTouchable);
0373     aParticleChange.AddSecondary(Track0);
0374     /*
0375       G4cout<<"Adding a particle "<<p0->GetDefinition()->GetParticleName()<<G4endl;
0376       G4cout<<"with momentum: "<<p0->GetMomentum()/GeV<<" GeV"<<G4endl;
0377       G4cout<<"and kinetic energy: "<<p0->GetKineticEnergy()/GeV<<" GeV"<<G4endl;
0378       */
0379     if (p0->GetKineticEnergy() > e_kin_0) {
0380       G4cout << "ALAAAAARM!!! (incident changed from " << IncidentRhadron->GetDefinition()->GetParticleName() << " to "
0381              << p0->GetDefinition()->GetParticleName() << ")" << G4endl;
0382       G4cout << "Energy loss: " << (e_kin_0 - p0->GetKineticEnergy()) / GeV << " GeV (should be positive)" << G4endl;
0383       //Turns out problem is only in 2 -> 3 (Won't fix 2 -> 2 energy deposition)
0384       if (rp.size() != 3)
0385         G4cout << "DOUBLE ALAAAAARM!!!" << G4endl;
0386     } /*else {
0387     G4cout<<"NO ALAAAAARM!!!"<<G4endl;
0388     }*/
0389     if (std::abs(p0->GetKineticEnergy() - e_kin_0) > 100 * GeV) {
0390       G4cout << "Diff. too big" << G4endl;
0391     }
0392     aParticleChange.ProposeTrackStatus(fStopAndKill);
0393   } else {
0394     G4double p = p_new.mag();
0395     if (p > DBL_MIN)
0396       aParticleChange.ProposeMomentumDirection(p_new.x() / p, p_new.y() / p, p_new.z() / p);
0397     else
0398       aParticleChange.ProposeMomentumDirection(1.0, 0.0, 0.0);
0399   }
0400 
0401   //    return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep);
0402   if (targetParticle.GetMass() > 0.0)  // targetParticle can be eliminated in TwoBody
0403   {
0404     G4DynamicParticle* p1 = new G4DynamicParticle;
0405     p1->SetDefinition(targetParticle.GetDefinition());
0406     //G4cout<<"Target secondary: "<<targetParticle.GetDefinition()->GetParticleName()<<G4endl;
0407     G4ThreeVector momentum = targetParticle.GetMomentum();
0408     momentum = momentum.rotate(cache, what);
0409     p1->SetMomentum(momentum);
0410     p1->SetMomentum((trans * p1->Get4Momentum()).vect());
0411     G4Track* Track1 = new G4Track(p1, aTrack.GetGlobalTime(), aPosition);
0412     Track1->SetTouchableHandle(thisTouchable);
0413     aParticleChange.AddSecondary(Track1);
0414   }
0415   G4DynamicParticle* pa;
0416   /*
0417     G4cout<<"vecLen: "<<vecLen<<G4endl;
0418     G4cout<<"#sec's: "<<aParticleChange.GetNumberOfSecondaries()<<G4endl;
0419   */
0420 
0421   for (int i = 0; i < vecLen; ++i) {
0422     pa = new G4DynamicParticle();
0423     pa->SetDefinition(vec[i]->GetDefinition());
0424     pa->SetMomentum(vec[i]->GetMomentum());
0425     pa->Set4Momentum(trans * (pa->Get4Momentum()));
0426     G4ThreeVector pvec = pa->GetMomentum();
0427     G4Track* Trackn = new G4Track(pa, aTrack.GetGlobalTime(), aPosition);
0428     Trackn->SetTouchableHandle(thisTouchable);
0429     aParticleChange.AddSecondary(Trackn);
0430 
0431     delete vec[i];
0432   }
0433 
0434   // Histogram filling
0435   const G4DynamicParticle* theRhadron = FindRhadron(&aParticleChange);
0436 
0437   if (theRhadron != nullptr || IncidentSurvives) {
0438     double E_new;
0439     if (IncidentSurvives) {
0440       E_new = e_kin;
0441     } else {
0442       E_new = theRhadron->GetKineticEnergy();
0443       if (CustomPDGParser::s_isRMeson(theRhadron->GetDefinition()->GetPDGEncoding()) !=
0444               CustomPDGParser::s_isRMeson(theIncidentPDG) ||
0445           CustomPDGParser::s_isMesonino(theRhadron->GetDefinition()->GetPDGEncoding()) !=
0446               CustomPDGParser::s_isMesonino(theIncidentPDG)) {
0447         G4cout << "Rm: " << CustomPDGParser::s_isRMeson(theRhadron->GetDefinition()->GetPDGEncoding())
0448                << " vs: " << CustomPDGParser::s_isRMeson(theIncidentPDG) << G4endl;
0449         G4cout << "Sm: " << CustomPDGParser::s_isMesonino(theRhadron->GetDefinition()->GetPDGEncoding())
0450                << " vs: " << CustomPDGParser::s_isMesonino(theIncidentPDG) << G4endl;
0451       }
0452     }
0453 
0454     //Calculating relevant scattering angles.
0455     G4LorentzVector p4_old_full = FullRhadron4Momentum;  //R-hadron in CMS BEFORE collision
0456     p4_old_full.boost(trafo_full_cms);
0457     G4LorentzVector p4_old_cloud = FullRhadron4Momentum;  //R-hadron in cloud CMS BEFORE collision
0458     p4_old_cloud.boost(trafo_cloud_cms);
0459     G4LorentzVector p4_full = p4_new;  //R-hadron in CMS AFTER collision
0460     //      G4cout<<p4_full.v()/GeV<<G4endl;
0461     p4_full = p4_full.boost(trafo_full_cms);
0462     // G4cout<<p4_full.m()<<" / "<<(cloud_p4_new+gluinoMomentum).boost(trafo_full_cms).m()<<G4endl;
0463     G4LorentzVector p4_cloud = p4_new;  //R-hadron in cloud CMS AFTER collision
0464     p4_cloud.boost(trafo_cloud_cms);
0465 
0466     G4double AbsDeltaE = E_0 - E_new;
0467     //      G4cout <<"Energy loss: "<<AbsDeltaE/GeV<<G4endl;
0468     if (AbsDeltaE > 10 * GeV) {
0469       G4cout << "Energy loss larger than 10 GeV..." << G4endl;
0470       G4cout << "E_0: " << E_0 / GeV << " GeV" << G4endl;
0471       G4cout << "E_new: " << E_new / GeV << " GeV" << G4endl;
0472       G4cout << "Gamma: " << IncidentRhadron->GetTotalEnergy() / IncidentRhadron->GetDefinition()->GetPDGMass()
0473              << G4endl;
0474       G4cout << "x: " << aPosition.x() / cm << " cm" << G4endl;
0475     }
0476   }
0477   delete originalIncident;
0478   delete originalTarget;
0479   //  aParticleChange.DumpInfo();
0480   //  G4cout << "Exiting FullModelHadronicProcess::PostStepDoIt"<<G4endl;
0481 
0482   //clear interaction length
0483   ClearNumberOfInteractionLengthLeft();
0484 
0485   return &aParticleChange;
0486 }
0487 
0488 void FullModelHadronicProcess::CalculateMomenta(
0489     G4FastVector<G4ReactionProduct, MYGHADLISTSIZE>& vec,
0490     G4int& vecLen,
0491     const G4HadProjectile* originalIncident,  // the original incident particle
0492     const G4DynamicParticle* originalTarget,
0493     G4ReactionProduct& modifiedOriginal,  // Fermi motion and evap. effects included
0494     G4Nucleus& targetNucleus,
0495     G4ReactionProduct& currentParticle,
0496     G4ReactionProduct& targetParticle,
0497     G4bool& incidentHasChanged,
0498     G4bool& targetHasChanged,
0499     G4bool quasiElastic) {
0500   FullModelReactionDynamics theReactionDynamics;
0501 
0502   cache = 0;
0503   what = originalIncident->Get4Momentum().vect();
0504 
0505   if (quasiElastic) {
0506     //      G4cout<<"We are calling TwoBody..."<<G4endl;
0507     theReactionDynamics.TwoBody(
0508         vec, vecLen, modifiedOriginal, originalTarget, currentParticle, targetParticle, targetNucleus, targetHasChanged);
0509 
0510     return;
0511   }
0512 
0513   //If ProduceStrangeParticlePairs is commented out, let's cut this one as well
0514   G4ReactionProduct leadingStrangeParticle;
0515   G4bool leadFlag = MarkLeadingStrangeParticle(currentParticle, targetParticle, leadingStrangeParticle);
0516 
0517   //
0518   // Note: the number of secondaries can be reduced in GenerateXandPt and TwoCluster
0519   //
0520   G4bool finishedGenXPt = false;
0521   G4bool annihilation = false;
0522   if (originalIncident->GetDefinition()->GetPDGEncoding() < 0 && currentParticle.GetMass() == 0.0 &&
0523       targetParticle.GetMass() == 0.0) {
0524     // original was an anti-particle and annihilation has taken place
0525     annihilation = true;
0526     G4double ekcor = 1.0;
0527     G4double ek = originalIncident->GetKineticEnergy();
0528     G4double ekOrg = ek;
0529 
0530     const G4double tarmas = originalTarget->GetDefinition()->GetPDGMass();
0531     if (ek > 1.0 * GeV)
0532       ekcor = 1. / (ek / GeV);
0533     const G4double atomicWeight = G4double(targetNucleus.GetN_asInt());
0534     ek = 2 * tarmas + ek * (1. + ekcor / atomicWeight);
0535 
0536     G4double tkin = targetNucleus.Cinema(ek);
0537     //ek += tkin;
0538     ekOrg += tkin;
0539     modifiedOriginal.SetKineticEnergy(ekOrg);
0540   }
0541 
0542   const G4double twsup[] = {1.0, 0.7, 0.5, 0.3, 0.2, 0.1};
0543   G4double rand1 = G4UniformRand();
0544   G4double rand2 = G4UniformRand();
0545   if ((annihilation || (vecLen >= 6) || (modifiedOriginal.GetKineticEnergy() / GeV >= 1.0)) &&
0546       (((originalIncident->GetDefinition() == G4KaonPlus::KaonPlus()) ||
0547         (originalIncident->GetDefinition() == G4KaonMinus::KaonMinus()) ||
0548         (originalIncident->GetDefinition() == G4KaonZeroLong::KaonZeroLong()) ||
0549         (originalIncident->GetDefinition() == G4KaonZeroShort::KaonZeroShort())) &&
0550        ((rand1 < 0.5) || (rand2 > twsup[vecLen]))))
0551     finishedGenXPt = theReactionDynamics.GenerateXandPt(vec,
0552                                                         vecLen,
0553                                                         modifiedOriginal,
0554                                                         originalIncident,
0555                                                         currentParticle,
0556                                                         targetParticle,
0557                                                         targetNucleus,
0558                                                         incidentHasChanged,
0559                                                         targetHasChanged,
0560                                                         leadFlag,
0561                                                         leadingStrangeParticle);
0562   if (finishedGenXPt) {
0563     Rotate(vec, vecLen);
0564     return;
0565   }
0566 
0567   G4bool finishedTwoClu = false;
0568   if (modifiedOriginal.GetTotalMomentum() / MeV < 1.0) {
0569     for (G4int i = 0; i < vecLen; i++)
0570       delete vec[i];
0571     vecLen = 0;
0572   } else {
0573     theReactionDynamics.SuppressChargedPions(vec,
0574                                              vecLen,
0575                                              modifiedOriginal,
0576                                              currentParticle,
0577                                              targetParticle,
0578                                              targetNucleus,
0579                                              incidentHasChanged,
0580                                              targetHasChanged);
0581 
0582     try {
0583       finishedTwoClu = theReactionDynamics.TwoCluster(vec,
0584                                                       vecLen,
0585                                                       modifiedOriginal,
0586                                                       originalIncident,
0587                                                       currentParticle,
0588                                                       targetParticle,
0589                                                       targetNucleus,
0590                                                       incidentHasChanged,
0591                                                       targetHasChanged,
0592                                                       leadFlag,
0593                                                       leadingStrangeParticle);
0594     } catch (G4HadronicException& aR) {
0595       G4ExceptionDescription ed;
0596       aR.Report(ed);
0597       G4Exception("FullModelHadronicProcess::CalculateMomenta", "had066", FatalException, ed);
0598     }
0599   }
0600   if (finishedTwoClu) {
0601     Rotate(vec, vecLen);
0602     return;
0603   }
0604 
0605   //
0606   // PNBlackTrackEnergy is the kinetic energy available for
0607   //   proton/neutron black track particles [was enp(1) in fortran code]
0608   // DTABlackTrackEnergy is the kinetic energy available for
0609   //   deuteron/triton/alpha particles      [was enp(3) in fortran code]
0610   // the atomic weight of the target nucleus is >= 1.5            AND
0611   //   neither the incident nor the target particles have changed  AND
0612   //     there is no kinetic energy available for either proton/neutron
0613   //     or for deuteron/triton/alpha black track particles
0614   // For diffraction scattering on heavy nuclei use elastic routines instead
0615 
0616   theReactionDynamics.TwoBody(
0617       vec, vecLen, modifiedOriginal, originalTarget, currentParticle, targetParticle, targetNucleus, targetHasChanged);
0618 }
0619 
0620 G4bool FullModelHadronicProcess::MarkLeadingStrangeParticle(const G4ReactionProduct& currentParticle,
0621                                                             const G4ReactionProduct& targetParticle,
0622                                                             G4ReactionProduct& leadParticle) {
0623   // the following was in GenerateXandPt and TwoCluster
0624   // add a parameter to the GenerateXandPt function telling it about the strange particle
0625   //
0626   // assumes that the original particle was a strange particle
0627   //
0628   G4bool lead = false;
0629   if ((currentParticle.GetMass() >= G4KaonPlus::KaonPlus()->GetPDGMass()) &&
0630       (currentParticle.GetDefinition() != G4Proton::Proton()) &&
0631       (currentParticle.GetDefinition() != G4Neutron::Neutron())) {
0632     lead = true;
0633     leadParticle = currentParticle;  //  set lead to the incident particle
0634   } else if ((targetParticle.GetMass() >= G4KaonPlus::KaonPlus()->GetPDGMass()) &&
0635              (targetParticle.GetDefinition() != G4Proton::Proton()) &&
0636              (targetParticle.GetDefinition() != G4Neutron::Neutron())) {
0637     lead = true;
0638     leadParticle = targetParticle;  //   set lead to the target particle
0639   }
0640   return lead;
0641 }
0642 
0643 void FullModelHadronicProcess::Rotate(G4FastVector<G4ReactionProduct, MYGHADLISTSIZE>& vec, G4int& vecLen) {
0644   G4double rotation = 2. * pi * G4UniformRand();
0645   cache = rotation;
0646   G4int i;
0647   for (i = 0; i < vecLen; ++i) {
0648     G4ThreeVector momentum = vec[i]->GetMomentum();
0649     momentum = momentum.rotate(rotation, what);
0650     vec[i]->SetMomentum(momentum);
0651   }
0652 }
0653 
0654 const G4DynamicParticle* FullModelHadronicProcess::FindRhadron(G4ParticleChange* aParticleChange) {
0655   G4int nsec = aParticleChange->GetNumberOfSecondaries();
0656   if (nsec == 0)
0657     return nullptr;
0658   int i = 0;
0659   G4bool found = false;
0660   while (i != nsec && !found) {
0661     //    G4cout<<"Checking "<<aParticleChange->GetSecondary(i)->GetDynamicParticle()->GetDefinition()->GetParticleName()<<G4endl;
0662     //    if (aParticleChange->GetSecondary(i)->GetDynamicParticle()->GetDefinition()->GetParticleType()=="rhadron") found = true;
0663     if (dynamic_cast<CustomParticle*>(aParticleChange->GetSecondary(i)->GetDynamicParticle()->GetDefinition()) !=
0664         nullptr)
0665       found = true;
0666     i++;
0667   }
0668   i--;
0669   if (found)
0670     return aParticleChange->GetSecondary(i)->GetDynamicParticle();
0671   return nullptr;
0672 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!