CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
FullModelHadronicProcess.cc
Go to the documentation of this file.
1 #include "G4HadReentrentException.hh"
2 #include "G4ProcessManager.hh"
3 #include "G4ParticleTable.hh"
4 
5 #include "SimG4Core/CustomPhysics/interface/FullModelHadronicProcess.hh"
6 #include "SimG4Core/CustomPhysics/interface/G4ProcessHelper.hh"
7 #include "SimG4Core/CustomPhysics/interface/Decay3Body.h"
8 #include "SimG4Core/CustomPhysics/interface/FullModelReactionDynamics.hh"
9 #include "SimG4Core/CustomPhysics/interface/CustomPDGParser.h"
10 #include "SimG4Core/CustomPhysics/interface/CustomParticle.h"
11 
12 FullModelHadronicProcess::FullModelHadronicProcess(G4ProcessHelper * aHelper, const G4String& processName) :
13  G4VDiscreteProcess(processName), theHelper(aHelper)
14 {}
15 
16 
17 FullModelHadronicProcess::~FullModelHadronicProcess(){}
18 
19 G4bool FullModelHadronicProcess::IsApplicable(const G4ParticleDefinition& aP)
20 {
21  return theHelper->ApplicabilityTester(aP);
22 }
23 
24 G4double FullModelHadronicProcess::GetMicroscopicCrossSection(const G4DynamicParticle *aParticle,
25  const G4Element *anElement,
26  G4double aTemp)
27 {
28  //Get the cross section for this particle/element combination from the ProcessHelper
29  G4double InclXsec = theHelper->GetInclusiveCrossSection(aParticle,anElement);
30  // G4cout<<"Returned cross section from helper was: "<<InclXsec/millibarn<<" millibarn"<<G4endl;
31 
32  // Need to provide Set-methods for these in time
33  G4double HighestEnergyLimit = 10 * TeV ;
34  G4double LowestEnergyLimit = 1 * eV;
35 
36  G4double ParticleEnergy = aParticle->GetKineticEnergy();
37 
38 
39  if (ParticleEnergy > HighestEnergyLimit || ParticleEnergy < LowestEnergyLimit){
40  return 0;
41  } else {
42  // G4cout << "Microscopic Cross Section: "<<InclXsec / millibarn<<" millibarn"<<G4endl;
43  return InclXsec;
44  }
45 
46 }
47 
48 G4double FullModelHadronicProcess::GetMeanFreePath(const G4Track& aTrack, G4double, G4ForceCondition*)
49 {
50  G4Material *aMaterial = aTrack.GetMaterial();
51  const G4DynamicParticle *aParticle = aTrack.GetDynamicParticle();
52  G4double sigma = 0.0;
53 
54  G4int nElements = aMaterial->GetNumberOfElements();
55 
56  const G4double *theAtomicNumDensityVector =
57  aMaterial->GetAtomicNumDensityVector();
58  G4double aTemp = aMaterial->GetTemperature();
59 
60  for( G4int i=0; i<nElements; ++i )
61  {
62  G4double xSection =
63  GetMicroscopicCrossSection( aParticle, (*aMaterial->GetElementVector())[i], aTemp);
64  sigma += theAtomicNumDensityVector[i] * xSection;
65  }
66 
67  return 1.0/sigma;
68 
69 }
70 
71 G4VParticleChange* FullModelHadronicProcess::PostStepDoIt(const G4Track& aTrack,
72  const G4Step& aStep)
73 {
74  // G4cout<<"***************** Entering FullModelHadronicProcess::PostStepDoIt **********************"<<G4endl;
75 
76  const G4TouchableHandle thisTouchable(aTrack.GetTouchableHandle());
77 
78  // A little setting up
79  aParticleChange.Initialize(aTrack);
80  // G4DynamicParticle* OrgPart = const_cast<G4DynamicParticle*>(aTrack.GetDynamicParticle());
81  G4DynamicParticle* IncidentRhadron = const_cast<G4DynamicParticle*>(aTrack.GetDynamicParticle());
82  CustomParticle* CustomIncident = static_cast<CustomParticle*>(IncidentRhadron->GetDefinition());
83  const G4ThreeVector aPosition = aTrack.GetPosition();
84  // std::cout<<"G: "<<aStep.GetStepLength()/cm<<std::endl;
85  const G4int theIncidentPDG = IncidentRhadron->GetDefinition()->GetPDGEncoding();
86  G4ParticleTable* theParticleTable = G4ParticleTable::GetParticleTable();
87  std::vector<G4ParticleDefinition*> theParticleDefinitions;
88  // std::vector<G4DynamicParticle*> *theDynamicParticles;//These are probably redundant, but they can easily be removed :-)
89  G4bool IncidentSurvives = false;
90  G4bool TargetSurvives = false;
91  G4Nucleus targetNucleus(aTrack.GetMaterial());
92  G4ParticleDefinition* outgoingRhadron=0;
93  G4ParticleDefinition* outgoingCloud=0;
94  G4ParticleDefinition* outgoingTarget=0;
95 // double gamma = IncidentRhadron->GetTotalEnergy()/IncidentRhadron->GetDefinition()->GetPDGMass();
96 
97 
98  G4ThreeVector p_0 = IncidentRhadron->GetMomentum();
99 
100  // static int n=0;
101 
102  G4double e_kin_0 = IncidentRhadron->GetKineticEnergy();
103  // G4cout<<e_kin_0/GeV<<G4endl;
104 
105  G4DynamicParticle* cloudParticle = new G4DynamicParticle();
106  /*
107  if(CustomPDGParser::s_isRMeson(theIncidentPDG))
108  std::cout<<"Rmeson"<<std::endl;
109  if(CustomPDGParser::s_isRBaryon(theIncidentPDG))
110  std::cout<<"Rbaryon"<<std::endl;
111  */
112  /*
113  if(CustomPDGParser::s_isRBaryon(theIncidentPDG))
114  cloudParticle->SetDefinition(theParticleTable->FindParticle("rhadronbaryoncloud"));
115  if(CustomPDGParser::s_isRMeson(theIncidentPDG) || CustomPDGParser::s_isRGlueball(theIncidentPDG) )
116  cloudParticle->SetDefinition(theParticleTable->FindParticle("rhadronmesoncloud"));
117  */
118  cloudParticle->SetDefinition(CustomIncident->GetCloud());
119 
120  if(cloudParticle->GetDefinition() == 0)
121  {
122  std::cout << "ToyModelHadronicProcess::PostStepDoIt Definition of particle cloud not available!!" << std::endl;
123  }
124  /*
125  std::cout<<"Incoming particle was "<<IncidentRhadron->GetDefinition()->GetParticleName()<<". Corresponding cloud is "<<cloudParticle->GetDefinition()->GetParticleName()<<std::endl;
126  G4cout<<"Kinetic energy was: "<<IncidentRhadron->GetKineticEnergy()/GeV<<" GeV"<<G4endl;
127  */
128  double scale=cloudParticle->GetDefinition()->GetPDGMass()/IncidentRhadron->GetDefinition()->GetPDGMass();
129  // std::cout<<"Mass ratio: "<<scale<<std::endl;
130  G4LorentzVector cloudMomentum;
131  cloudMomentum.setVectM(IncidentRhadron->GetMomentum()*scale,cloudParticle->GetDefinition()->GetPDGMass());
132  G4LorentzVector gluinoMomentum;
133  // gluinoMomentum.setVectM(IncidentRhadron->GetMomentum()*(1.-scale),theParticleTable->FindParticle("~g")->GetPDGMass());
134  gluinoMomentum.setVectM(IncidentRhadron->GetMomentum()*(1.-scale),CustomIncident->GetSpectator()->GetPDGMass());
135  /*
136  G4cout <<"Are these the same?"<<G4endl;
137  G4cout<<gluinoMomentum<<G4endl;
138  G4cout<<(1.-scale) * IncidentRhadron->Get4Momentum()<<G4endl;
139  */
140  //These two for getting CMS transforms later (histogramming purposes...)
141  G4LorentzVector FullRhadron4Momentum = IncidentRhadron->Get4Momentum();
142  G4LorentzVector Cloud4Momentum = cloudMomentum;
143 
144  cloudParticle->Set4Momentum(cloudMomentum);
145 
146  G4DynamicParticle* OrgPart = cloudParticle;
147 
148  /*
149  std::cout<<"Original momentum: "<<IncidentRhadron->Get4Momentum().v().mag()/GeV<<" GeV, corresponding to gamma: "
150  <<IncidentRhadron->GetTotalEnergy()/IncidentRhadron->GetDefinition()->GetPDGMass()<<std::endl;
151 
152  std::cout<<"Cloud momentum: "<<cloudParticle->Get4Momentum().v().mag()/GeV<<" GeV, corresponding to gamma: "
153  <<cloudParticle->GetTotalEnergy()/cloudParticle->GetDefinition()->GetPDGMass()<<std::endl;
154  */
155 
156  double E_0 = IncidentRhadron->GetKineticEnergy();
157 
158 
159  G4double ek = OrgPart->GetKineticEnergy();
160  G4double amas = OrgPart->GetDefinition()->GetPDGMass();
161 
162  G4double tkin = targetNucleus.Cinema( ek );
163  ek += tkin;
164  OrgPart->SetKineticEnergy( ek );
165  G4double et = ek + amas;
166  G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
167  G4double pp = OrgPart->GetMomentum().mag();
168  if( pp > 0.0 )
169  {
170  G4ThreeVector momentum = OrgPart->GetMomentum();
171  OrgPart->SetMomentum( momentum * (p/pp) );
172  }
173 
174  // calculate black track energies
175 
176  tkin = targetNucleus.EvaporationEffects( ek );
177  ek -= tkin;
178  if(ek+gluinoMomentum.e()-gluinoMomentum.m()<=0.1*MeV||ek<=0.) {
179  //Very rare event...
180  G4cout<<"Kinetic energy is sick"<<G4endl;
181  G4cout<<"Full R-hadron: "<<(ek+gluinoMomentum.e()-gluinoMomentum.m())/MeV<<" MeV" <<G4endl;
182  G4cout<<"Quark system: "<<ek/MeV<<" MeV"<<G4endl;
183  aParticleChange.ProposeTrackStatus( fStopAndKill );
184  return &aParticleChange;
185  }
186  OrgPart->SetKineticEnergy( ek );
187  et = ek + amas;
188  p = std::sqrt( std::abs((et-amas)*(et+amas)) );
189  pp = OrgPart->GetMomentum().mag();
190 
191  if( pp > 0.0 )
192  {
193  G4ThreeVector momentum = OrgPart->GetMomentum();
194  OrgPart->SetMomentum( momentum * (p/pp) );
195  }
196 
197 
198 
199  //Get the final state particles
200 
201  G4ParticleDefinition* aTarget;
202  ReactionProduct rp = theHelper->GetFinalState(aTrack,aTarget);
203  G4bool force2to2 = false;
204  // G4cout<<"Trying to get final state..."<<G4endl;
205  while(rp.size()!=2 && force2to2){
206  rp = theHelper->GetFinalState(aTrack,aTarget);
207  }
208  G4int NumberOfSecondaries = rp.size();
209  // G4cout<<"Multiplicity of selected final state: "<<rp.size()<<G4endl;
210 
211  //Getting CMS transforms. Boosting is done at histogram filling
212  G4LorentzVector Target4Momentum;
213  Target4Momentum.setVectM(0.,aTarget->GetPDGMass());
214  // Target4Momentum.setVectM(0.,targetNucleus.GetN()*GeV);
215  G4LorentzVector psum_full,psum_cloud;
216  psum_full = FullRhadron4Momentum + Target4Momentum;
217  psum_cloud = Cloud4Momentum + Target4Momentum;
218  G4ThreeVector trafo_full_cms = (-1)*psum_full.boostVector();
219  G4ThreeVector trafo_cloud_cms = (-1)*psum_cloud.boostVector();
220  /*
221  psum_full.boost(trafo_full_cms);
222  psum_cloud.boost(trafo_cloud_cms);
223  std::cout<<"Checking that the momenta are in deed zero:"<<psum_full.vect()<<std::endl;
224  */
225 
226  // OK Let's make some particles :-)
227  // We're not using them for anything yet, I know, but let's make sure the machinery is there
228 
229  for(ReactionProduct::iterator it = rp.begin();
230  it != rp.end();
231  it++)
232  {
233  G4ParticleDefinition* tempDef = theParticleTable->FindParticle(*it);
234  CustomParticle* tempCust = dynamic_cast<CustomParticle*>(tempDef);
235  if (tempDef==aTarget) TargetSurvives = true;
236 
237  // if (tempDef->GetParticleType()=="rhadron")
238  if (tempCust!=0)
239  {
240  outgoingRhadron = tempDef;
241  //Setting outgoing cloud definition
242  /*
243  if(CustomPDGParser::s_isRBaryon(*it))
244  outgoingCloud=theParticleTable->FindParticle("rhadronbaryoncloud");
245  if(CustomPDGParser::s_isRMeson(*it) || CustomPDGParser::s_isRGlueball(*it) )
246  outgoingCloud=theParticleTable->FindParticle("rhadronmesoncloud");
247  */
248  outgoingCloud=tempCust->GetCloud();
249  if(outgoingCloud == 0)
250  {
251  std::cout << "ToyModelHadronicProcess::PostStepDoIt Definition of outgoing particle cloud not available!!" << std::endl;
252  }
253  /*
254  std::cout<<"Outgoing Rhadron is: "<<outgoingRhadron->GetParticleName()<<std::endl;
255  std::cout<<"Outgoing cloud is: "<<outgoingCloud->GetParticleName()<<std::endl;
256  */
257  }
258 
259  if (tempDef==G4Proton::Proton()||tempDef==G4Neutron::Neutron()) outgoingTarget = tempDef;
260  // if (tempDef->GetParticleType()!="rhadron"&&rp.size()==2) outgoingTarget = tempDef;
261  if (tempCust==0&&rp.size()==2) outgoingTarget = tempDef;
262  if (tempDef->GetPDGEncoding()==theIncidentPDG) {
263  IncidentSurvives = true;
264  } else {
265  theParticleDefinitions.push_back(tempDef);
266  /*
267  G4DynamicParticle* tempDyn = new G4DynamicParticle();
268  tempDyn->SetDefinition(tempDef);
269  theDynamicParticles->push_back(tempDyn);
270  */
271  }
272  }
273 
274  //Not using this, so...
275  // delete theDynamicParticles;
276 
277  if (outgoingTarget==0) outgoingTarget = theParticleTable->FindParticle(rp[1]);
278 
279  // A little debug information
280  /*
281  G4cout<<"The particles coming out of this reaction will be: ";
282  for (std::vector<G4DynamicParticle*>::iterator it = theDynamicParticles.begin();
283  it != theDynamicParticles.end();
284  it++){
285  G4cout<< (*it)->GetDefinition()->GetParticleName()<<" ";
286  }
287  G4cout<<G4endl;
288  */
289  // If the incident particle survives it is not a "secondary", although
290  // it would probably be easier to fStopAndKill the track every time.
291  if(IncidentSurvives) NumberOfSecondaries--;
292  aParticleChange.SetNumberOfSecondaries(NumberOfSecondaries);
293 
294 
295  // ADAPTED FROM G4LEPionMinusInelastic::ApplyYourself
296  // It is bloody ugly, but such is the way of cut 'n' paste
297 
298 
299  // Set up the incident
300  const G4HadProjectile* originalIncident = new G4HadProjectile(*OrgPart);//This is where rotation to z-axis is done (Aarrggh!)
301 
302 
303  //Maybe I need to calculate trafo from R-hadron... Bollocks! Labframe does not move. CMS does.
304  G4HadProjectile* org2 = new G4HadProjectile(*OrgPart);
305  G4LorentzRotation trans = org2->GetTrafoToLab();
306  delete org2;
307 
308  // if (originalIncident->GetKineticEnergy()<= 0.1*MeV) { //Needs rescaling. The kinetic energy of the quarksystem is the relevant quantity
309 
310  /*
311  G4cout<<"Kinetic energies: "<<G4endl;
312  G4cout<<"True kinetic energy: "<<originalIncident->GetKineticEnergy()/GeV<<" GeV"<<G4endl;
313  G4cout<<"Mass: "<<originalIncident->GetDefinition()->GetPDGMass()/GeV<<" GeV"<<G4endl;
314 
315  G4double e_kin_rescaled = targetNucleus.EvaporationEffects(originalIncident->GetTotalEnergy()-originalIncident->GetDefinition()->GetPDGMass());
316 
317  G4cout<<"Rescaled kinetic energy: "<<e_kin_rescaled<<G4endl;
318 
319  const G4double cutOff = 0.1*MeV;
320 
321  if ( e_kin_rescaled < cutOff )
322  {
323  aParticleChange.ProposeTrackStatus( fStopAndKill );//If the dice decides not to cascade I stop the particle
324  return &aParticleChange;
325  }
326  */
327  // create the target particle
328 
329  G4DynamicParticle *originalTarget = new G4DynamicParticle;
330  originalTarget->SetDefinition(aTarget);
331 
332  G4ReactionProduct targetParticle(aTarget);
333 
334 
335  G4ReactionProduct currentParticle(const_cast<G4ParticleDefinition *>(originalIncident->GetDefinition() ) );
336  currentParticle.SetMomentum( originalIncident->Get4Momentum().vect() );
337  currentParticle.SetKineticEnergy( originalIncident->GetKineticEnergy() );
338 
339  /*
340  G4cout<<"After creation:"<<G4endl;
341  G4cout<<"currentParticle: "<<currentParticle.GetMomentum()/GeV<<" GeV vs. "<<OrgPart->Get4Momentum()/GeV<<" GeV"<<G4endl;
342  G4cout<<"targetParticle: "<<targetParticle.GetMomentum()/GeV<<" GeV"<<G4endl;
343  G4cout<<"Fourmomentum from originalIncident: "<<originalIncident->Get4Momentum()<<" vs "<<OrgPart->Get4Momentum()<<G4endl;
344  */
345 
346 
347  G4ReactionProduct modifiedOriginal = currentParticle;
348 
349  currentParticle.SetSide( 1 ); // incident always goes in forward hemisphere
350  targetParticle.SetSide( -1 ); // target always goes in backward hemisphere
351  G4bool incidentHasChanged = false;
352  if (!IncidentSurvives) incidentHasChanged = true; //I wonder if I am supposed to do this...
353  G4bool targetHasChanged = false;
354  if (!TargetSurvives) targetHasChanged = true; //Ditto here
355  G4bool quasiElastic = false;
356  if (rp.size()==2) quasiElastic = true; //Oh well...
357  G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec; // vec will contain the secondary particles
358  G4int vecLen = 0;
359  vec.Initialize( 0 );
360 
361 
362 
363  // I hope my understanding of "secondary" is correct here
364  // I think that it entails only what is not a surviving incident of target
365 
366  for (G4int i = 0; i!=NumberOfSecondaries;i++){
367  if(theParticleDefinitions[i]!=aTarget
368  && theParticleDefinitions[i]!=originalIncident->GetDefinition()
369  && theParticleDefinitions[i]!=outgoingRhadron
370  && theParticleDefinitions[i]!=outgoingTarget)
371  {
372  G4ReactionProduct* pa = new G4ReactionProduct;
373  pa->SetDefinition( theParticleDefinitions[i] );
374  (G4UniformRand() < 0.5) ? pa->SetSide( -1 ) : pa->SetSide( 1 );
375  vec.SetElement( vecLen++, pa );
376  }
377  }
378 
379  // if (incidentHasChanged) currentParticle.SetDefinitionAndUpdateE( outgoingRhadron );
380 
381  if (incidentHasChanged) currentParticle.SetDefinitionAndUpdateE( outgoingCloud );
382  if (incidentHasChanged) modifiedOriginal.SetDefinition( outgoingCloud );//Is this correct? It solves the "free energy" checking in ReactionDynamics
383  if (targetHasChanged) targetParticle.SetDefinitionAndUpdateE( outgoingTarget );
384 
385  // G4cout<<"Calling CalculateMomenta... "<<G4endl;
386  /*
387  G4cout<<"Current particle starts as: "<<currentParticle.GetDefinition()->GetParticleName()<<G4endl;
388  G4cout<<"with momentum: "<<currentParticle.GetMomentum()/GeV<<" GeV"<<G4endl;
389  G4cout<<"May be killed?: "<<currentParticle.GetMayBeKilled()<<G4endl;
390  */
391  // G4double e_temp = currentParticle.GetKineticEnergy();
392 
393  CalculateMomenta( vec, vecLen,
394  originalIncident, originalTarget, modifiedOriginal,
395  targetNucleus, currentParticle, targetParticle,
396  incidentHasChanged, targetHasChanged, quasiElastic );
397 
398  // G4cout <<"Cloud loss: "<<(e_temp-currentParticle.GetKineticEnergy())/GeV<<" GeV"<<G4endl;
399 
400  G4String cPname = currentParticle.GetDefinition()->GetParticleName();
401 
402  // if(cPname!="rhadronmesoncloud"&&cPname!="rhadronbaryoncloud")
403  // {
404  /*
405  G4cout<<"Current particle is now: "<<cPname <<G4endl;
406  G4cout<<"with momentum: "<<currentParticle.GetMomentum()/GeV<<" GeV"<<G4endl;
407  G4cout<<"and kinetic energy: "<<currentParticle.GetKineticEnergy()/GeV<<" GeV"<<G4endl;
408  G4cout<<"May be killed?: "<<currentParticle.GetMayBeKilled()<<G4endl;
409  G4cout<<"Modified original is: "<<modifiedOriginal.GetDefinition()->GetParticleName()<<G4endl;
410  G4cout<<"with momentum: "<<modifiedOriginal.GetMomentum()/GeV<<" GeV"<<G4endl;
411  G4cout<<"and kinetic energy: "<<modifiedOriginal.GetKineticEnergy()/GeV<<" GeV"<<G4endl;
412  G4cout<<"May be killed?: "<<modifiedOriginal.GetMayBeKilled()<<G4endl;
413  G4cout<<"Target particle is: "<<targetParticle.GetDefinition()->GetParticleName()<<G4endl;
414  G4cout<<"with momentum: "<<targetParticle.GetMomentum()/GeV<<" GeV"<<G4endl;
415  G4cout<<"and kinetic energy: "<<targetParticle.GetKineticEnergy()/GeV<<" GeV"<<G4endl;
416  G4cout<<"May be killed?: "<<targetParticle.GetMayBeKilled()<<G4endl;
417  G4cout<<"incidentHasChanged: "<<incidentHasChanged<<G4endl;
418  G4cout<<"targetHasChanged: "<<targetHasChanged<<G4endl;
419  G4cout<<"Particles in vec:"<<G4endl;
420  for(int i=0; i<vecLen; ++i )
421  {
422  G4cout<< vec[i]->GetDefinition()->GetParticleName()<<G4endl;
423  }
424  */
425 
426  // }
427 
428  // G4cout<<"Done!"<<G4endl;
429  // const G4LorentzRotation& trans(originalIncident->GetTrafoToLab());
430  // G4cout<<"Check aParticleChange.GetNumberOfSecondaries(): "<<aParticleChange.GetNumberOfSecondaries()<<G4endl;
431 
432  aParticleChange.SetNumberOfSecondaries(vecLen+NumberOfSecondaries);
433  G4double e_kin=0;
434  G4LorentzVector cloud_p4_new; //Cloud 4-momentum in lab after collision
435  // n++;
436  // G4cout << n << G4endl;
437  /*
438  if(cPname!="rhadronmesoncloud"&&cPname!="rhadronbaryoncloud") {
439  G4cout<<"Cloud deleted!!! AAARRRRGGGHHH!!!"<<G4endl;
440  G4cout<<"Cloud name: "<<cPname<<G4endl;
441  G4cout<<"E_kin_0: "<<e_kin_0/GeV<<" GeV"<<G4endl;
442  // G4cout<<"n: "<<n<<G4endl;
443  // n=0;
444  }
445  */
446  cloud_p4_new.setVectM(currentParticle.GetMomentum(),currentParticle.GetMass());
447  cloud_p4_new *= trans;
448 
449  G4LorentzVector cloud_p4_old_full = Cloud4Momentum; //quark system in CMS BEFORE collision
450  cloud_p4_old_full.boost(trafo_full_cms);
451  G4LorentzVector cloud_p4_old_cloud = Cloud4Momentum; //quark system in cloud CMS BEFORE collision
452  cloud_p4_old_cloud.boost(trafo_cloud_cms);
453  G4LorentzVector cloud_p4_full = cloud_p4_new; //quark system in CMS AFTER collision
454  cloud_p4_full.boost(trafo_full_cms);
455  G4LorentzVector cloud_p4_cloud = cloud_p4_new; //quark system in cloud CMS AFTER collision
456  cloud_p4_cloud.boost(trafo_cloud_cms);
457 
458  G4LorentzVector p_g_cms = gluinoMomentum; //gluino in CMS BEFORE collision
459  p_g_cms.boost(trafo_full_cms);
460 
461  G4LorentzVector p4_new;
462  // p4_new.setVectM(cloud_p4_full.v()+p_g_cms.v(),outgoingRhadron->GetPDGMass());
463  // p4_new.boost(-trafo_full_cms);
464  // p4_new = cloud_p4_new + gluinoMomentum;
465  p4_new.setVectM( cloud_p4_new.v() + gluinoMomentum.v(), outgoingRhadron->GetPDGMass() );
466  // G4cout<<"P4-diff: "<<(p4_new-cloud_p4_new-gluinoMomentum)/GeV<<", magnitude: "<<(p4_new-cloud_p4_new-gluinoMomentum).m()/MeV<<" MeV" <<G4endl;
467 
468  G4ThreeVector p_new;
469  p_new = p4_new.vect();
470 
471  aParticleChange.ProposeLocalEnergyDeposit((p4_new-cloud_p4_new-gluinoMomentum).m());
472 
473  if( incidentHasChanged )
474  {
475  G4DynamicParticle* p0 = new G4DynamicParticle;
476  p0->SetDefinition( outgoingRhadron );
477  p0->SetMomentum( p_new );
478 
479  // May need to run SetDefinitionAndUpdateE here...
480  G4Track* Track0 = new G4Track(p0,
481  aTrack.GetGlobalTime(),
482  aPosition);
483  Track0->SetTouchableHandle(thisTouchable);
484  aParticleChange.AddSecondary(Track0);
485  /*
486  G4cout<<"Adding a particle "<<p0->GetDefinition()->GetParticleName()<<G4endl;
487  G4cout<<"with momentum: "<<p0->GetMomentum()/GeV<<" GeV"<<G4endl;
488  G4cout<<"and kinetic energy: "<<p0->GetKineticEnergy()/GeV<<" GeV"<<G4endl;
489  */
490  if(p0->GetKineticEnergy()>e_kin_0) {
491  G4cout<<"ALAAAAARM!!! (incident changed from "
492  <<IncidentRhadron->GetDefinition()->GetParticleName()
493  <<" to "<<p0->GetDefinition()->GetParticleName()<<")"<<G4endl;
494  G4cout<<"Energy loss: "<<(e_kin_0-p0->GetKineticEnergy())/GeV<<" GeV (should be positive)"<<G4endl;
495  //Turns out problem is only in 2 -> 3 (Won't fix 2 -> 2 energy deposition)
496  if(rp.size()!=3) G4cout<<"DOUBLE ALAAAAARM!!!"<<G4endl;
497  } /*else {
498  G4cout<<"NO ALAAAAARM!!!"<<G4endl;
499  }*/
500  if(std::abs(p0->GetKineticEnergy()-e_kin_0)>100*GeV) {
501  G4cout<<"Diff. too big"<<G4endl;
502  }
503 
504  aParticleChange.ProposeTrackStatus( fStopAndKill );
505  }
506  else
507  {
508 
509  G4double p = p_new.mag();
510  if( p > DBL_MIN )
511  aParticleChange.ProposeMomentumDirection( p_new.x()/p, p_new.y()/p, p_new.z()/p );
512  else
513  aParticleChange.ProposeMomentumDirection( 1.0, 0.0, 0.0 );
514 
515  G4double aE = sqrt(p*p+(outgoingRhadron->GetPDGMass()*outgoingRhadron->GetPDGMass()) );
516  e_kin = aE - outgoingRhadron->GetPDGMass();
517  /*
518  G4cout<<"New momentum: "<<m/GeV<<" GeV"<<G4endl;
519  G4cout<<"Kinetic energy: "<<e_kin/GeV<<" GeV"<<G4endl;
520  */
521  if(e_kin>e_kin_0) {
522  G4cout<<"ALAAAAARM!!!"<<G4endl;
523  G4cout<<"Energy loss: "<<(e_kin_0-e_kin)/GeV<<" GeV (should be positive)"<<G4endl;
524  if(rp.size()!=3) G4cout<<"DOUBLE ALAAAAARM!!!"<<G4endl;
525  }
526  if(std::abs(e_kin-e_kin_0)>100*GeV) {
527  G4cout<<"Diff. too big"<<G4endl;
528  }
529 
530  if (std::fabs(aE)<.1*eV) aE=.1*eV;
531  aParticleChange.ProposeEnergy( aE- outgoingRhadron->GetPDGMass() ); //I do not know if this is correct...
532  if(std::abs(e_kin-e_kin_0)>100*GeV) {
533  G4cout<<"Diff. too big"<<G4endl;
534  }
535 
536  }
537 
538  // return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep);
539  if( targetParticle.GetMass() > 0.0 ) // targetParticle can be eliminated in TwoBody
540  {
541  G4DynamicParticle *p1 = new G4DynamicParticle;
542  p1->SetDefinition( targetParticle.GetDefinition() );
543  // G4cout<<"Target secondary: "<<targetParticle.GetDefinition()->GetParticleName()<<G4endl;
544  G4ThreeVector momentum = targetParticle.GetMomentum();
545  momentum = momentum.rotate(cache,what);
546  p1->SetMomentum( momentum );
547  p1->SetMomentum( (trans*p1->Get4Momentum()).vect());
548  G4Track* Track1 = new G4Track(p1,
549  aTrack.GetGlobalTime(),
550  aPosition);
551  Track1->SetTouchableHandle(thisTouchable);
552  aParticleChange.AddSecondary(Track1);
553  }
554  G4DynamicParticle *pa;
555  /*
556  G4cout<<"vecLen: "<<vecLen<<G4endl;
557  G4cout<<"#sec's: "<<aParticleChange.GetNumberOfSecondaries()<<G4endl;
558  */
559 
560 
561 
562  for(int i=0; i<vecLen; ++i )
563  {
564  pa = new G4DynamicParticle();
565  pa->SetDefinition( vec[i]->GetDefinition() );
566  pa->SetMomentum( vec[i]->GetMomentum() );
567  pa->Set4Momentum(trans*(pa->Get4Momentum()));
568  G4ThreeVector pvec = pa->GetMomentum();
569  G4Track* Trackn = new G4Track(pa,
570  aTrack.GetGlobalTime(),
571  aPosition);
572  Trackn->SetTouchableHandle(thisTouchable);
573  aParticleChange.AddSecondary(Trackn);
574 
575  // debug
576 
577 // G4cerr << "FullModelHadronicProcess: New secondary " << i
578 // << " ID " << Trackn->GetTrackID()
579 // << " PDG " << Trackn->GetDefinition()->GetParticleName()
580 // << " position " << Trackn->GetPosition()
581 // << " volume " << Trackn->GetTouchable()
582 // << " handle " << Trackn->GetTouchableHandle() << G4endl;
583 
584  delete vec[i];
585  }
586 
587  // Histogram filling
588  const G4DynamicParticle* theRhadron = FindRhadron(&aParticleChange);
589 
590 
591  if (theRhadron!=NULL||IncidentSurvives)
592  {
593 
594  double E_new;
595  if(IncidentSurvives)
596  {
597  // E_new = currentParticle.GetKineticEnergy();
598  E_new = e_kin;
599  } else {
600  E_new = theRhadron->GetKineticEnergy();
601  if(CustomPDGParser::s_isRMeson(theRhadron->GetDefinition()->GetPDGEncoding())
602  !=CustomPDGParser::s_isRMeson(theIncidentPDG)
603  ||
604  CustomPDGParser::s_isMesonino(theRhadron->GetDefinition()->GetPDGEncoding())
605  !=CustomPDGParser::s_isMesonino(theIncidentPDG)
606  ) {
607 
608  G4cout<<"Rm: "<<CustomPDGParser::s_isRMeson(theRhadron->GetDefinition()->GetPDGEncoding())
609  <<" vs: "<<CustomPDGParser::s_isRMeson(theIncidentPDG)<<G4endl;
610  G4cout<<"Sm: "<<CustomPDGParser::s_isMesonino(theRhadron->GetDefinition()->GetPDGEncoding())
611  <<" vs: "<<CustomPDGParser::s_isMesonino(theIncidentPDG)<<G4endl;
612 
613  }
614  }
615 
616  //Calculating relevant scattering angles.
617  G4LorentzVector p4_old_full = FullRhadron4Momentum; //R-hadron in CMS BEFORE collision
618  p4_old_full.boost(trafo_full_cms);
619  G4LorentzVector p4_old_cloud = FullRhadron4Momentum; //R-hadron in cloud CMS BEFORE collision
620  p4_old_cloud.boost(trafo_cloud_cms);
621  G4LorentzVector p4_full = p4_new; //R-hadron in CMS AFTER collision
622  // G4cout<<p4_full.v()/GeV<<G4endl;
623  p4_full=p4_full.boost(trafo_full_cms);
624  // G4cout<<p4_full.m()<<" / "<<(cloud_p4_new+gluinoMomentum).boost(trafo_full_cms).m()<<G4endl;
625  G4LorentzVector p4_cloud = p4_new; //R-hadron in cloud CMS AFTER collision
626  p4_cloud.boost(trafo_cloud_cms);
627 
628 
629  /*
630  G4double dtheta_fullcms = p4_full.vect().angle(p4_old_full.vect());
631  G4double dtheta_cloudcms = p4_cloud.vect().angle(p4_old_cloud.vect());
632  G4double dtheta_lab = p_new.angle(p_0);//acos(p_0*p_new/(p_0.mag()*p_new.mag()));
633 
634  G4double cloud_dtheta_fullcms = cloud_p4_full.vect().angle(cloud_p4_old_full.vect());
635  G4double cloud_dtheta_cloudcms = cloud_p4_cloud.vect().angle(p4_old_cloud.vect());
636  G4double cloud_dtheta_lab = cloud_p4_new.vect().angle(p_0);
637  //Writing out momenta for manual check of boosts:
638  G4cout<<"******************************************"<<G4endl;
639 
640  G4cout<<"R-hadron, before: "<<G4endl;
641  G4cout<<"Lab: "<<FullRhadron4Momentum.v().mag()/GeV<<" GeV"<<G4endl;
642  G4cout<<"CMS: "<<p4_old_full.v().mag()/GeV<<" GeV"<<G4endl;
643  G4cout<<"after: "<<G4endl;
644  G4cout<<"Lab: "<<p4_new.v().mag()/GeV<<" GeV"<<G4endl;
645  G4cout<<"CMS: "<<p4_full.v().mag()/GeV<<" GeV"<<G4endl;
646  G4cout<<"cos(theta*): "<<cos(dtheta_fullcms)<<G4endl;
647  G4cout<<"Gluino: "<<G4endl;
648  G4cout<<"Lab: "<<gluinoMomentum.v().mag()/GeV<<" GeV"<<G4endl;
649  G4cout<<"CMS: "<<p_g_cms.v().mag()/GeV<<" GeV"<<G4endl;
650 
651  G4cout<<"Cloud, before: "<<G4endl;
652  G4cout<<"Lab: "<<Cloud4Momentum.v().mag()/GeV<<" GeV"<<G4endl;
653  G4cout<<"CMS: "<<cloud_p4_old_full.v().mag()/GeV<<" GeV"<<G4endl;
654  G4cout<<"after: "<<G4endl;
655  G4cout<<"CMS: "<<cloud_p4_full.v().mag()/GeV<<" GeV"<<G4endl;
656  G4cout<<"cloud cos(theta*): "<<cos(cloud_dtheta_fullcms)<<G4endl;
657  G4cout<<"Longitudinal: "<<cos(cloud_dtheta_fullcms)*cloud_p4_full.v().mag()/GeV<<" GeV"<<G4endl;
658  G4cout<<"~Combined longitudinal: "<<(cos(cloud_dtheta_fullcms)*cloud_p4_full.v().mag()+p4_old_full.v().mag())/GeV<<" GeV"<<G4endl;
659  if ((cos(cloud_dtheta_fullcms)*cloud_p4_full.v().mag()+p4_old_full.v().mag())<0.) G4cout<<"ALAAARM"<<G4endl;
660  G4cout<<"Psum cos(theta*): "<<cos((p_g_cms.v()+cloud_p4_full.v()).angle(p4_old_full.v()))<<G4endl;
661  G4cout<<"True R-hadron (CMS): "<<(p_g_cms.v()+cloud_p4_full.v()).mag()/GeV<<" GeV"<<G4endl;
662  G4cout<<"******************************************"<<G4endl;
663  G4cout<<"Fucking manual fucking calculation:"<<G4endl;
664  G4cout<<"Momenta:"<<G4endl;
665  G4cout<<"Cloud, lab: "<<cloud_p4_new.v()/GeV<<" + gluino: "<<gluinoMomentum.v()/GeV
666  <<" = "<<(cloud_p4_new.v()+gluinoMomentum.v())/GeV
667  <<". Boost to CMS: "<<(cloud_p4_new+gluinoMomentum).boost(trafo_full_cms).v()/GeV<<G4endl;
668  G4cout<<"Cloud, CMS: "<<cloud_p4_full.v()/GeV<<" + gluino: "<<p_g_cms.v()/GeV
669  <<" = "<<(cloud_p4_full.v()+p_g_cms.v())/GeV
670  <<". Boost to Lab: "<<(cloud_p4_full+p_g_cms).boost(-trafo_full_cms).v()/GeV<<G4endl;
671  G4cout<<"Ref: "<<p4_new.v()/GeV<<" / "<<p4_full.v()/GeV<<G4endl;
672  G4cout<<"******************************************"<<G4endl;
673  */
674 
675 
676  // std::cout<<"Lab, fullcms, cloudcms: "<<dtheta_lab<<", "<<dtheta_fullcms<<", "<<dtheta_cloudcms<<std::endl;
677  G4double AbsDeltaE = E_0-E_new;
678  // G4cout <<"Energy loss: "<<AbsDeltaE/GeV<<G4endl;
679  if (AbsDeltaE > 10*GeV) {
680  G4cout<<"Energy loss larger than 10 GeV..."<<G4endl;
681  G4cout<<"E_0: "<<E_0/GeV<<" GeV"<<G4endl;
682  G4cout<<"E_new: "<<E_new/GeV<<" GeV"<<G4endl;
683  G4cout<<"Gamma: "<<IncidentRhadron->GetTotalEnergy()/IncidentRhadron->GetDefinition()->GetPDGMass()<<G4endl;
684  G4cout<<"x: "<<aPosition.x()/cm<<" cm"<<G4endl;
685  }
686  }
687  delete originalIncident;
688  delete originalTarget;
689  // aParticleChange.DumpInfo();
690  // G4cout << "Exiting FullModelHadronicProcess::PostStepDoIt"<<G4endl;
691 
692  //clear interaction length
693  ClearNumberOfInteractionLengthLeft();
694 
695 
696  return &aParticleChange;
697 
698 }
699 
700 
701 void FullModelHadronicProcess::CalculateMomenta(
702  G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
703  G4int &vecLen,
704  const G4HadProjectile *originalIncident, // the original incident particle
705  const G4DynamicParticle *originalTarget,
706  G4ReactionProduct &modifiedOriginal, // Fermi motion and evap. effects included
707  G4Nucleus &targetNucleus,
708  G4ReactionProduct &currentParticle,
709  G4ReactionProduct &targetParticle,
710  G4bool &incidentHasChanged,
711  G4bool &targetHasChanged,
712  G4bool quasiElastic )
713 {
714  FullModelReactionDynamics theReactionDynamics;
715 
716  cache = 0;
717  what = originalIncident->Get4Momentum().vect();
718 
719  //Commented out like in casqmesmin.F where CALL STPAIR is commented out
720  /*
721  theReactionDynamics.ProduceStrangeParticlePairs( vec, vecLen,
722  modifiedOriginal, originalTarget,
723  currentParticle, targetParticle,
724  incidentHasChanged, targetHasChanged );
725  */
726 
727  if( quasiElastic )
728  {
729  // G4cout<<"We are calling TwoBody..."<<G4endl;
730  theReactionDynamics.TwoBody( vec, vecLen,
731  modifiedOriginal, originalTarget,
732  currentParticle, targetParticle,
733  targetNucleus, targetHasChanged );
734 
735  return;
736  }
737 
738  //If ProduceStrangeParticlePairs is commented out, let's cut this one as well
739  G4ReactionProduct leadingStrangeParticle;
740  G4bool leadFlag = MarkLeadingStrangeParticle( currentParticle,
741  targetParticle,
742  leadingStrangeParticle );
743 
744 
745 
746  //
747  // Note: the number of secondaries can be reduced in GenerateXandPt and TwoCluster
748  //
749  G4bool finishedGenXPt = false;
750  G4bool annihilation = false;
751  if( originalIncident->GetDefinition()->GetPDGEncoding() < 0 &&
752  currentParticle.GetMass() == 0.0 && targetParticle.GetMass() == 0.0 )
753  {
754  // original was an anti-particle and annihilation has taken place
755  annihilation = true;
756  G4double ekcor = 1.0;
757  G4double ek = originalIncident->GetKineticEnergy();
758  G4double ekOrg = ek;
759 
760  const G4double tarmas = originalTarget->GetDefinition()->GetPDGMass();
761  if( ek > 1.0*GeV )ekcor = 1./(ek/GeV);
762  const G4double atomicWeight = targetNucleus.GetN();
763  ek = 2*tarmas + ek*(1.+ekcor/atomicWeight);
764 
765  G4double tkin = targetNucleus.Cinema( ek );
766  ek += tkin;
767  ekOrg += tkin;
768  modifiedOriginal.SetKineticEnergy( ekOrg );
769  }
770 
771  const G4double twsup[] = { 1.0, 0.7, 0.5, 0.3, 0.2, 0.1 };
772  G4double rand1 = G4UniformRand();
773  G4double rand2 = G4UniformRand();
774  if( (annihilation || (vecLen >= 6) || (modifiedOriginal.GetKineticEnergy()/GeV >= 1.0)) &&
775  (((originalIncident->GetDefinition() == G4KaonPlus::KaonPlus()) ||
776  (originalIncident->GetDefinition() == G4KaonMinus::KaonMinus()) ||
777  (originalIncident->GetDefinition() == G4KaonZeroLong::KaonZeroLong()) ||
778  (originalIncident->GetDefinition() == G4KaonZeroShort::KaonZeroShort())) &&
779  ((rand1 < 0.5) || (rand2 > twsup[vecLen]))))
780  finishedGenXPt =
781  theReactionDynamics.GenerateXandPt( vec, vecLen,
782  modifiedOriginal, originalIncident,
783  currentParticle, targetParticle,
784  targetNucleus, incidentHasChanged,
785  targetHasChanged, leadFlag,
786  leadingStrangeParticle );
787  if( finishedGenXPt )
788  {
789  Rotate(vec, vecLen);
790  return;
791  }
792 
793  G4bool finishedTwoClu = false;
794  if( modifiedOriginal.GetTotalMomentum()/MeV < 1.0 )
795  {
796 
797  for(G4int i=0; i<vecLen; i++) delete vec[i];
798  vecLen = 0;
799  }
800  else
801  {
802 
803  theReactionDynamics.SuppressChargedPions( vec, vecLen,
804  modifiedOriginal, currentParticle,
805  targetParticle, targetNucleus,
806  incidentHasChanged, targetHasChanged );
807 
808  try
809  {
810  finishedTwoClu = theReactionDynamics.TwoCluster( vec, vecLen,
811  modifiedOriginal, originalIncident,
812  currentParticle, targetParticle,
813  targetNucleus, incidentHasChanged,
814  targetHasChanged, leadFlag,
815  leadingStrangeParticle );
816  }
817  catch(G4HadReentrentException aC)
818  {
819  aC.Report(G4cout);
820  throw G4HadReentrentException(__FILE__, __LINE__, "Failing to calculate momenta");
821  }
822  }
823  if( finishedTwoClu )
824  {
825  Rotate(vec, vecLen);
826  return;
827  }
828 
829  //
830  // PNBlackTrackEnergy is the kinetic energy available for
831  // proton/neutron black track particles [was enp(1) in fortran code]
832  // DTABlackTrackEnergy is the kinetic energy available for
833  // deuteron/triton/alpha particles [was enp(3) in fortran code]
834  //const G4double pnCutOff = 0.1;
835  //const G4double dtaCutOff = 0.1;
836  //if( (targetNucleus.GetN() >= 1.5)
837  // && !(incidentHasChanged || targetHasChanged)
838  // && (targetNucleus.GetPNBlackTrackEnergy()/MeV <= pnCutOff)
839  // && (targetNucleus.GetDTABlackTrackEnergy()/MeV <= dtaCutOff) )
840  //{
841  // the atomic weight of the target nucleus is >= 1.5 AND
842  // neither the incident nor the target particles have changed AND
843  // there is no kinetic energy available for either proton/neutron
844  // or for deuteron/triton/alpha black track particles
845  // For diffraction scattering on heavy nuclei use elastic routines instead
846  //G4cerr << "*** Error in G4InelasticInteraction::CalculateMomenta" << G4endl;
847  //G4cerr << "*** the elastic scattering would be better here ***" <<G4endl;
848  //}
849  theReactionDynamics.TwoBody( vec, vecLen,
850  modifiedOriginal, originalTarget,
851  currentParticle, targetParticle,
852  targetNucleus, targetHasChanged );
853 }
854 
855 
856 G4bool FullModelHadronicProcess::MarkLeadingStrangeParticle(
857  const G4ReactionProduct &currentParticle,
858  const G4ReactionProduct &targetParticle,
859  G4ReactionProduct &leadParticle )
860 {
861  // the following was in GenerateXandPt and TwoCluster
862  // add a parameter to the GenerateXandPt function telling it about the strange particle
863  //
864  // assumes that the original particle was a strange particle
865  //
866  G4bool lead = false;
867  if( (currentParticle.GetMass() >= G4KaonPlus::KaonPlus()->GetPDGMass()) &&
868  (currentParticle.GetDefinition() != G4Proton::Proton()) &&
869  (currentParticle.GetDefinition() != G4Neutron::Neutron()) )
870  {
871  lead = true;
872  leadParticle = currentParticle; // set lead to the incident particle
873  }
874  else if( (targetParticle.GetMass() >= G4KaonPlus::KaonPlus()->GetPDGMass()) &&
875  (targetParticle.GetDefinition() != G4Proton::Proton()) &&
876  (targetParticle.GetDefinition() != G4Neutron::Neutron()) )
877  {
878  lead = true;
879  leadParticle = targetParticle; // set lead to the target particle
880  }
881  return lead;
882 }
883 
884 void FullModelHadronicProcess::Rotate(G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec, G4int &vecLen)
885 {
886  G4double rotation = 2.*pi*G4UniformRand();
887  cache = rotation;
888  G4int i;
889  for( i=0; i<vecLen; ++i )
890  {
891  G4ThreeVector momentum = vec[i]->GetMomentum();
892  momentum = momentum.rotate(rotation, what);
893  vec[i]->SetMomentum(momentum);
894  }
895 }
896 
897 const G4DynamicParticle* FullModelHadronicProcess::FindRhadron(G4ParticleChange* aParticleChange)
898 {
899  G4int nsec = aParticleChange->GetNumberOfSecondaries();
900  if (nsec==0) return 0;
901  int i = 0;
902  G4bool found = false;
903  while (i!=nsec && !found){
904  // G4cout<<"Checking "<<aParticleChange->GetSecondary(i)->GetDynamicParticle()->GetDefinition()->GetParticleName()<<G4endl;
905  // if (aParticleChange->GetSecondary(i)->GetDynamicParticle()->GetDefinition()->GetParticleType()=="rhadron") found = true;
906  if (dynamic_cast<CustomParticle*>(aParticleChange->GetSecondary(i)->GetDynamicParticle()->GetDefinition())!=0) found = true;
907  i++;
908  }
909  i--;
910  if(found) return aParticleChange->GetSecondary(i)->GetDynamicParticle();
911  return 0;
912 }
913 
i
int i
Definition: DBlmapReader.cc:9
createTree.pp
tuple pp
Definition: createTree.py:15
abs
#define abs(x)
Definition: mlp_lapack.h:159
CustomParticle::GetCloud
G4ParticleDefinition * GetCloud()
Definition: CustomParticle.h:38
NULL
#define NULL
Definition: scimark2.h:8
CustomPDGParser::s_isRMeson
static bool s_isRMeson(int pdg)
Definition: CustomPDGParser.cc:43
CustomPDGParser::s_isMesonino
static bool s_isMesonino(int pdg)
Definition: CustomPDGParser.cc:50
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:28
CustomParticle::GetSpectator
G4ParticleDefinition * GetSpectator()
Definition: CustomParticle.h:40
L1TEmulatorMonitor_cff.p
tuple p
Definition: L1TEmulatorMonitor_cff.py:48
p1
double p1[4]
Definition: TauolaWrapper.h:89
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:41
pi
double pi
Definition: Vx3DHLTAnalyzer.h:63
Generated for CMSSW Reference Manual by doxygen 1.8.5