CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
G4ProcessHelper.cc
Go to the documentation of this file.
1 #include"G4ParticleTable.hh"
2 #include "Randomize.hh"
3 
4 #include<iostream>
5 #include<fstream>
6 #include <string>
7 
11 
12 #include"SimG4Core/CustomPhysics/interface/G4ProcessHelper.hh"
15 
16 G4ProcessHelper::G4ProcessHelper(const edm::ParameterSet & p){
17 
18  particleTable = G4ParticleTable::GetParticleTable();
19 
20  theProton = particleTable->FindParticle("proton");
21  theNeutron = particleTable->FindParticle("neutron");
22 
23  G4String line;
24 
25  edm::FileInPath fp = p.getParameter<edm::FileInPath>("processesDef");
26  std::string processDefFilePath = fp.fullPath();
27  std::ifstream process_stream (processDefFilePath.c_str());
28 
29  resonant = p.getParameter<bool>("resonant");
30  ek_0 = p.getParameter<double>("resonanceEnergy")*GeV;
31  gamma = p.getParameter<double>("gamma")*GeV;
32  amplitude = p.getParameter<double>("amplitude")*millibarn;
33  suppressionfactor = p.getParameter<double>("reggeSuppression");
34  hadronlifetime = p.getParameter<double>("hadronLifeTime");
35  reggemodel = p.getParameter<bool>("reggeModel");
36  mixing = p.getParameter<double>("mixing");
37 
38  edm::LogInfo("CustomPhysics")<<"Read in physics parameters:"<<G4endl;
39  edm::LogInfo("CustomPhysics")<<"Resonant = "<< resonant <<G4endl;
40  edm::LogInfo("CustomPhysics")<<"ResonanceEnergy = "<<ek_0/GeV<<" GeV"<<G4endl;
41  edm::LogInfo("CustomPhysics")<<"Gamma = "<<gamma/GeV<<" GeV"<<G4endl;
42  edm::LogInfo("CustomPhysics")<<"Amplitude = "<<amplitude/millibarn<<" millibarn"<<G4endl;
43  edm::LogInfo("CustomPhysics")<<"ReggeSuppression = "<<100*suppressionfactor<<" %"<<G4endl;
44  edm::LogInfo("CustomPhysics")<<"HadronLifeTime = "<<hadronlifetime<<" s"<<G4endl;
45  edm::LogInfo("CustomPhysics")<<"ReggeModel = "<< reggemodel <<G4endl;
46  edm::LogInfo("CustomPhysics")<<"Mixing = "<< mixing*100 <<" %"<<G4endl;
47 
48 
49  checkfraction = 0;
50  n_22 = 0;
51  n_23 = 0;
52 
53 
54 
55  while(getline(process_stream,line)){
56  std::vector<G4String> tokens;
57  //Getting a line
58  ReadAndParse(line,tokens,"#");
59  //Important info
60  G4String incident = tokens[0];
61 
62  G4ParticleDefinition* incidentDef = particleTable->FindParticle(incident);
63  //particleTable->DumpTable();
64  G4int incidentPDG = incidentDef->GetPDGEncoding();
65  known_particles[incidentDef]=true;
66 
67  G4String target = tokens[1];
68  edm::LogInfo("CustomPhysics")<<"Incident: "<<incident
69  <<" Target: "<<target<<G4endl;
70 
71  // Making a ReactionProduct
72  ReactionProduct prod;
73  for (size_t i = 2; i != tokens.size();i++){
74  G4String part = tokens[i];
75  if (particleTable->contains(part))
76  {
77  prod.push_back(particleTable->FindParticle(part)->GetPDGEncoding());
78  } else {
79  edm::LogInfo("CustomPhysics")<<"Particle: "<<part<<" is unknown."<<G4endl;
80  G4Exception("G4ProcessHelper", "UnkownParticle", FatalException,
81  "Initialization: The reaction product list contained an unknown particle");
82  }
83  }
84  if (target == "proton")
85  {
86  pReactionMap[incidentPDG].push_back(prod);
87  } else if (target == "neutron") {
88  nReactionMap[incidentPDG].push_back(prod);
89  } else {
90  G4Exception("G4ProcessHelper", "IllegalTarget", FatalException,
91  "Initialization: The reaction product list contained an illegal target particle");
92  }
93  }
94 
95  process_stream.close();
96 
97  G4ParticleTable::G4PTblDicIterator* theParticleIterator;
98  theParticleIterator = particleTable->GetIterator();
99 
100  theParticleIterator->reset();
101  while( (*theParticleIterator)() ){
102  CustomParticle* particle = dynamic_cast<CustomParticle*>(theParticleIterator->value());
103  std::string name = theParticleIterator->value()->GetParticleName();
104  G4DecayTable* table = theParticleIterator->value()->GetDecayTable();
105  if(particle!=0&&table!=0&&name.find("cloud")>name.size()&&hadronlifetime > 0)
106  {
107  particle->SetPDGLifeTime(hadronlifetime*s);
108  particle->SetPDGStable(false);
109  edm::LogInfo("CustomPhysics")<<"Lifetime of: "<<name<<" set to: "<<particle->GetPDGLifeTime()/s<<" s."<<G4endl;
110  edm::LogInfo("CustomPhysics")<<"Stable: "<<particle->GetPDGStable()<<G4endl;
111  }
112  }
113  theParticleIterator->reset();
114 
115 
116 }
117 
118 G4bool G4ProcessHelper::ApplicabilityTester(const G4ParticleDefinition& aPart){
119  const G4ParticleDefinition* aP = &aPart;
120  if (known_particles[aP]) return true;
121  return false;
122 }
123 
124 G4double G4ProcessHelper::GetInclusiveCrossSection(const G4DynamicParticle *aParticle,
125  const G4Element *anElement){
126 
127  //We really do need a dedicated class to handle the cross sections. They might not always be constant
128 
129 
130  //Disassemble the PDG-code
131 
132  G4int thePDGCode = aParticle->GetDefinition()->GetPDGEncoding();
133  double boost = (aParticle->GetKineticEnergy()+aParticle->GetMass())/aParticle->GetMass();
134  // G4cout<<"thePDGCode: "<<thePDGCode<<G4endl;
135  G4double theXsec = 0;
136  G4String name = aParticle->GetDefinition()->GetParticleName();
137  if(!reggemodel)
138  {
139  //Flat cross section
140  if(CustomPDGParser::s_isRGlueball(thePDGCode)) {
141  theXsec = 24 * millibarn;
142  } else {
143  std::vector<G4int> nq=CustomPDGParser::s_containedQuarks(thePDGCode);
144  // edm::LogInfo("CustomPhysics")<<"Number of quarks: "<<nq.size()<<G4endl;
145  for (std::vector<G4int>::iterator it = nq.begin();
146  it != nq.end();
147  it++)
148  {
149  // edm::LogInfo("CustomPhysics")<<"Quarkvector: "<<*it<<G4endl;
150  if (*it == 1 || *it == 2) theXsec += 12 * millibarn;
151  if (*it == 3) theXsec += 6 * millibarn;
152  }
153  }
154  }
155  else
156  { //reggemodel
157  double R = Regge(boost);
158  double P = Pom(boost);
159  if(thePDGCode>0)
160  {
161  if(CustomPDGParser::s_isMesonino(thePDGCode)) theXsec=(P+R)*millibarn;
162  if(CustomPDGParser::s_isSbaryon(thePDGCode)) theXsec=2*P*millibarn;
163  if(CustomPDGParser::s_isRMeson(thePDGCode)||CustomPDGParser::s_isRGlueball(thePDGCode)) theXsec=(R+2*P)*millibarn;
164  if(CustomPDGParser::s_isRBaryon(thePDGCode)) theXsec=3*P*millibarn;
165  }
166  else
167  {
168  if(CustomPDGParser::s_isMesonino(thePDGCode)) theXsec=P*millibarn;
169  if(CustomPDGParser::s_isSbaryon(thePDGCode)) theXsec=(2*(P+R)+30/sqrt(boost))*millibarn;
170  if(CustomPDGParser::s_isRMeson(thePDGCode)||CustomPDGParser::s_isRGlueball(thePDGCode)) theXsec=(R+2*P)*millibarn;
171  if(CustomPDGParser::s_isRBaryon(thePDGCode)) theXsec=3*P*millibarn;
172  }
173  }
174 
175 
176 
177  //Adding resonance
178 
179  if(resonant)
180  {
181  double e_0 = ek_0 + aParticle->GetDefinition()->GetPDGMass(); //Now total energy
182 
183  e_0 = sqrt(aParticle->GetDefinition()->GetPDGMass()*aParticle->GetDefinition()->GetPDGMass()
184  + theProton->GetPDGMass()*theProton->GetPDGMass()
185  + 2.*e_0*theProton->GetPDGMass());
186  // edm::LogInfo("CustomPhysics")<<e_0/GeV<<G4endl;
187  // edm::LogInfo("CustomPhysics")<<ek_0/GeV<<" "<<aParticle->GetDefinition()->GetPDGMass()/GeV<<" "<<theProton->GetPDGMass()/GeV<<G4endl;
188  double sqrts=sqrt(aParticle->GetDefinition()->GetPDGMass()*aParticle->GetDefinition()->GetPDGMass()
189  + theProton->GetPDGMass()*theProton->GetPDGMass() + 2*aParticle->GetTotalEnergy()*theProton->GetPDGMass());
190 
191  double res_result = amplitude*(gamma*gamma/4.)/((sqrts-e_0)*(sqrts-e_0)+(gamma*gamma/4.));//Non-relativistic Breit Wigner
192 
193  theXsec += res_result;
194  // if(fabs(aParticle->GetKineticEnergy()/GeV-200)<10) std::cout<<sqrts/GeV<<" "<<theXsec/millibarn<<std::endl;
195  }
196 
197 
198  // std::cout<<"Xsec/nucleon: "<<theXsec/millibarn<<"millibarn, Total Xsec: "<<theXsec * anElement->GetN()/millibarn<<" millibarn"<<std::endl;
199  // return theXsec * anElement->GetN();// * 0.523598775598299;
200  return theXsec * pow(anElement->GetN(),0.7)*1.25;// * 0.523598775598299;
201 
202 }
203 
204 ReactionProduct G4ProcessHelper::GetFinalState(const G4Track& aTrack, G4ParticleDefinition*& aTarget){
205 
206  const G4DynamicParticle* aDynamicParticle = aTrack.GetDynamicParticle();
207 
208  //-----------------------------------------------
209  // Choose n / p as target
210  // and get ReactionProductList pointer
211  //-----------------------------------------------
212 
213  G4Material* aMaterial = aTrack.GetMaterial();
214  const G4ElementVector* theElementVector = aMaterial->GetElementVector() ;
215  const G4double* NbOfAtomsPerVolume = aMaterial->GetVecNbOfAtomsPerVolume();
216 
217  G4double NumberOfProtons=0;
218  G4double NumberOfNucleons=0;
219 
220  for ( size_t elm=0 ; elm < aMaterial->GetNumberOfElements() ; elm++ )
221  {
222  //Summing number of protons per unit volume
223  NumberOfProtons += NbOfAtomsPerVolume[elm]*(*theElementVector)[elm]->GetZ();
224  //Summing nucleons (not neutrons)
225  NumberOfNucleons += NbOfAtomsPerVolume[elm]*(*theElementVector)[elm]->GetN();
226  }
227 
228  if(G4UniformRand()<NumberOfProtons/NumberOfNucleons)
229  {
230  theReactionMap = &pReactionMap;
231  theTarget = theProton;
232  } else {
233  theReactionMap = &nReactionMap;
234  theTarget = theNeutron;
235  }
236  aTarget = theTarget;
237 
238  G4int theIncidentPDG = aDynamicParticle->GetDefinition()->GetPDGEncoding();
239 
240  if(reggemodel
241  &&CustomPDGParser::s_isMesonino(theIncidentPDG)
242  && G4UniformRand()*mixing>0.5
243  &&aDynamicParticle->GetDefinition()->GetPDGCharge()==0.
244  )
245  {
246  // G4cout<<"Oscillating..."<<G4endl;
247  theIncidentPDG *= -1;
248  }
249 
250 
251  bool baryonise=false;
252 
253  if(reggemodel
254  && G4UniformRand()>0.9
255  &&(
256  (CustomPDGParser::s_isMesonino(theIncidentPDG)&&theIncidentPDG>0)
257  ||
258  CustomPDGParser::s_isRMeson(theIncidentPDG)
259  )
260  )
261  baryonise=true;
262 
263  //Making a pointer directly to the ReactionProductList we are looking at. Makes life easier :-)
264  ReactionProductList* aReactionProductList = &((*theReactionMap)[theIncidentPDG]);
265 
266  //-----------------------------------------------
267  // Count processes
268  // kinematic check
269  // compute number of 2 -> 2 and 2 -> 3 processes
270  //-----------------------------------------------
271 
272  G4int N22 = 0; //Number of 2 -> 2 processes
273  G4int N23 = 0; //Number of 2 -> 3 processes. Couldn't think of more informative names
274 
275  //This is the list to be populated
276  ReactionProductList theReactionProductList;
277  std::vector<bool> theChargeChangeList;
278 
279  for (ReactionProductList::iterator prod_it = aReactionProductList->begin();
280  prod_it != aReactionProductList->end();
281  prod_it++){
282  G4int secondaries = prod_it->size();
283  // If the reaction is not possible we will not consider it
284 // if(ReactionIsPossible(*prod_it,aDynamicParticle)){
285  if(ReactionIsPossible(*prod_it,aDynamicParticle)
286  && (!reggemodel ||
287  (baryonise&&ReactionGivesBaryon(*prod_it))
288  ||
289  (!baryonise&&!ReactionGivesBaryon(*prod_it))
290  ||
291  (CustomPDGParser::s_isSbaryon(theIncidentPDG))
292  ||
293  (CustomPDGParser::s_isRBaryon(theIncidentPDG))
294  )
295  )
296  {
297  // The reaction is possible. Let's store and count it
298  theReactionProductList.push_back(*prod_it);
299  if (secondaries == 2){
300  N22++;
301  } else if (secondaries ==3) {
302  N23++;
303  } else {
304  G4cerr << "ReactionProduct has unsupported number of secondaries: "<<secondaries<<G4endl;
305  }
306  } /*else {
307  edm::LogInfo("CustomPhysics")<<"There was an impossible process"<<G4endl;
308  }*/
309  }
310  // edm::LogInfo("CustomPhysics")<<"The size of the ReactionProductList is: "<<theReactionProductList.size()<<G4endl;
311 
312  if (theReactionProductList.size()==0) G4Exception("G4ProcessHelper", "NoProcessPossible", FatalException,
313  "GetFinalState: No process could be selected from the given list.");
314 
315  // For the Regge model no phase space considerations. We pick a process at random
316  if(reggemodel)
317  {
318  int n_rps = theReactionProductList.size();
319  int select = (int)( G4UniformRand()*n_rps);
320  // G4cout<<"Possible: "<<n_rps<<", chosen: "<<select<<G4endl;
321  return theReactionProductList[select];
322  }
323 
324  // Fill a probability map. Remember total probability
325  // 2->2 is 0.15*1/n_22 2->3 uses phase space
326  G4double p22 = 0.15;
327  G4double p23 = 1-p22; // :-)
328 
329  std::vector<G4double> Probabilities;
330  std::vector<G4bool> TwotoThreeFlag;
331 
332  G4double CumulatedProbability = 0;
333 
334  // To each ReactionProduct we assign a cumulated probability and a flag
335  // discerning between 2 -> 2 and 2 -> 3
336  for (unsigned int i = 0; i != theReactionProductList.size(); i++){
337  if (theReactionProductList[i].size() == 2) {
338  CumulatedProbability += p22/N22;
339  TwotoThreeFlag.push_back(false);
340  } else {
341  CumulatedProbability += p23/N23;
342  TwotoThreeFlag.push_back(true);
343  }
344  Probabilities.push_back(CumulatedProbability);
345  // edm::LogInfo("CustomPhysics")<<"Pushing back cumulated probability: "<<CumulatedProbability<<G4endl;
346  }
347 
348  //Renormalising probabilities
349  // edm::LogInfo("CustomPhysics")<<"Probs: ";
350  for (std::vector<G4double>::iterator it = Probabilities.begin();
351  it != Probabilities.end();
352  it++)
353  {
354  *it /= CumulatedProbability;
355  // edm::LogInfo("CustomPhysics")<<*it<<" ";
356  }
357  // edm::LogInfo("CustomPhysics")<<G4endl;
358 
359  // Choosing ReactionProduct
360 
361  G4bool selected = false;
362  G4int tries = 0;
363  // ReactionProductList::iterator prod_it;
364 
365  //Keep looping over the list until we have a choice, or until we have tried 100 times
366  unsigned int i;
367  while(!selected && tries < 100){
368  i=0;
369  G4double dice = G4UniformRand();
370  // edm::LogInfo("CustomPhysics")<<"What's the dice?"<<dice<<G4endl;
371  while(dice>Probabilities[i] && i<theReactionProductList.size()){
372  // edm::LogInfo("CustomPhysics")<<"i: "<<i<<G4endl;
373  i++;
374  }
375 
376  // edm::LogInfo("CustomPhysics")<<"Chosen i: "<<i<<G4endl;
377 
378  if(!TwotoThreeFlag[i]) {
379  // 2 -> 2 processes are chosen immediately
380  selected = true;
381  } else {
382  // 2 -> 3 processes require a phase space lookup
383  if (PhaseSpace(theReactionProductList[i],aDynamicParticle)>G4UniformRand()) selected = true;
384  //selected = true;
385  }
386  // double suppressionfactor=0.5;
387  if(selected&&particleTable->FindParticle(theReactionProductList[i][0])->GetPDGCharge()!=aDynamicParticle->GetDefinition()->GetPDGCharge())
388  {
389  /*
390  edm::LogInfo("CustomPhysics")<<"Incoming particle "<<aDynamicParticle->GetDefinition()->GetParticleName()
391  <<" has charge "<<aDynamicParticle->GetDefinition()->GetPDGCharge()<<G4endl;
392  edm::LogInfo("CustomPhysics")<<"Suggested particle "<<particleTable->FindParticle(theReactionProductList[i][0])->GetParticleName()
393  <<" has charge "<<particleTable->FindParticle(theReactionProductList[i][0])->GetPDGCharge()<<G4endl;
394  */
395  if(G4UniformRand()<suppressionfactor) selected = false;
396  }
397  tries++;
398  // edm::LogInfo("CustomPhysics")<<"Tries: "<<tries<<G4endl;
399  }
400  if(tries>=100) G4cerr<<"Could not select process!!!!"<<G4endl;
401 
402  // edm::LogInfo("CustomPhysics")<<"So far so good"<<G4endl;
403  // edm::LogInfo("CustomPhysics")<<"Sec's: "<<theReactionProductList[i].size()<<G4endl;
404 
405  //Updating checkfraction:
406  if (theReactionProductList[i].size()==2) {
407  n_22++;
408  } else {
409  n_23++;
410  }
411 
412  checkfraction = (1.0*n_22)/(n_22+n_23);
413  // edm::LogInfo("CustomPhysics")<<"n_22: "<<n_22<<" n_23: "<<n_23<<" Checkfraction: "<<checkfraction<<G4endl;
414  // edm::LogInfo("CustomPhysics") <<"Biig number: "<<n_22+n_23<<G4endl;
415  //Return the chosen ReactionProduct
416  return theReactionProductList[i];
417 }
418 
419 G4double G4ProcessHelper::ReactionProductMass(const ReactionProduct& aReaction,const G4DynamicParticle* aDynamicParticle){
420  // Incident energy:
421  G4double E_incident = aDynamicParticle->GetTotalEnergy();
422  //edm::LogInfo("CustomPhysics")<<"Total energy: "<<E_incident<<" Kinetic: "<<aDynamicParticle->GetKineticEnergy()<<G4endl;
423  // sqrt(s)= sqrt(m_1^2 + m_2^2 + 2 E_1 m_2)
424  G4double m_1 = aDynamicParticle->GetDefinition()->GetPDGMass();
425  G4double m_2 = theTarget->GetPDGMass();
426  //edm::LogInfo("CustomPhysics")<<"M_R: "<<m_1/GeV<<" GeV, M_np: "<<m_2/GeV<<" GeV"<<G4endl;
427  G4double sqrts = sqrt(m_1*m_1 + m_2*(m_2 + 2 * E_incident));
428  //edm::LogInfo("CustomPhysics")<<"sqrt(s) = "<<sqrts/GeV<<" GeV"<<G4endl;
429  // Sum of rest masses after reaction:
430  G4double M_after = 0;
431  for (ReactionProduct::const_iterator r_it = aReaction.begin(); r_it !=aReaction.end(); r_it++){
432  //edm::LogInfo("CustomPhysics")<<"Mass contrib: "<<(particleTable->FindParticle(*r_it)->GetPDGMass())/MeV<<" MeV"<<G4endl;
433  M_after += particleTable->FindParticle(*r_it)->GetPDGMass();
434  }
435  //edm::LogInfo("CustomPhysics")<<"Intending to return this ReactionProductMass: "<<(sqrts - M_after)/MeV<<" MeV"<<G4endl;
436  return sqrts - M_after;
437 }
438 
439 G4bool G4ProcessHelper::ReactionIsPossible(const ReactionProduct& aReaction,const G4DynamicParticle* aDynamicParticle){
440  if (ReactionProductMass(aReaction,aDynamicParticle)>0) return true;
441  return false;
442 }
443 
444 G4bool G4ProcessHelper::ReactionGivesBaryon(const ReactionProduct& aReaction){
445  for (ReactionProduct::const_iterator it = aReaction.begin();it!=aReaction.end();it++)
447  return false;
448 }
449 
450 G4double G4ProcessHelper::PhaseSpace(const ReactionProduct& aReaction,const G4DynamicParticle* aDynamicParticle){
451  G4double qValue = ReactionProductMass(aReaction,aDynamicParticle);
452  G4double phi = sqrt(1+qValue/(2*0.139*GeV))*pow(qValue/(1.1*GeV),3./2.);
453  return (phi/(1+phi));
454 }
455 
456 void G4ProcessHelper::ReadAndParse(const G4String& str,
457  std::vector<G4String>& tokens,
458  const G4String& delimiters)
459 {
460  // Skip delimiters at beginning.
461  G4String::size_type lastPos = str.find_first_not_of(delimiters, 0);
462  // Find first "non-delimiter".
463  G4String::size_type pos = str.find_first_of(delimiters, lastPos);
464 
465  while (G4String::npos != pos || G4String::npos != lastPos)
466  {
467  //Skipping leading / trailing whitespaces
468  G4String temp = str.substr(lastPos, pos - lastPos);
469  while(temp.c_str()[0] == ' ') temp.erase(0,1);
470  while(temp[temp.size()-1] == ' ') temp.erase(temp.size()-1,1);
471  // Found a token, add it to the vector.
472  tokens.push_back(temp);
473  // Skip delimiters. Note the "not_of"
474  lastPos = str.find_first_not_of(delimiters, pos);
475  // Find next "non-delimiter"
476  pos = str.find_first_of(delimiters, lastPos);
477  }
478 }
479 
480 double G4ProcessHelper::Regge(const double boost)
481 {
482  double a=2.165635078566177;
483  double b=0.1467453738547229;
484  double c=-0.9607903711871166;
485  return 1.5*exp(a+b/boost+c*log(boost));
486 }
487 
488 
489 double G4ProcessHelper::Pom(const double boost)
490 {
491  double a=4.138224000651535;
492  double b=1.50377557581421;
493  double c=-0.05449742257808247;
494  double d=0.0008221235048211401;
495  return a + b*sqrt(boost) + c*boost + d*pow(boost,1.5);
496 }
497 
list table
Definition: asciidump.py:386
T getParameter(std::string const &) const
int i
Definition: DBlmapReader.cc:9
static bool s_isRGlueball(int pdg)
static bool s_isSbaryon(int pdg)
#define P
static bool s_isRMeson(int pdg)
uint16_t size_type
static bool s_isMesonino(int pdg)
T sqrt(T t)
Definition: SSEVec.h:46
static bool s_isRBaryon(int pdg)
part
Definition: HCALResponse.h:21
static std::vector< int > s_containedQuarks(int pdg)
double b
Definition: hdecay.h:120
double a
Definition: hdecay.h:121
std::string fullPath() const
Definition: FileInPath.cc:171
tuple size
Write out results.
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
Definition: DDAxes.h:10