QELEventGeneratorSM.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2017, GENIE Neutrino MC Generator Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Igor Kakorin <kakorin@jinr.ru>, Joint Institute for Nuclear Research
8  adopted from fortran code provided by
9  Konstantin Kuzmin <kkuzmin@theor.jinr.ru>,
10  Joint Institute for Nuclear Research, Institute for Theoretical and Experimental Physics
11  Vladimir Lyubushkin,
12  Joint Institute for Nuclear Research
13  Vadim Naumov <vnaumov@theor.jinr.ru>,
14  Joint Institute for Nuclear Research
15  based on code of Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
16  University of Liverpool & STFC Rutherford Appleton Lab
17 
18  For the class documentation see the corresponding header file.
19 */
20 //____________________________________________________________________________
21 
22 #include <TMath.h>
23 
41 
42 
43 #include "Framework/Utils/Range1.h"
46 #include "Framework/Utils/Cache.h"
48 
49 using namespace genie;
50 using namespace genie::controls;
51 using namespace genie::utils;
52 
53 namespace { // anonymous namespace (file only visibility)
54  const double eps = std::numeric_limits<double>::epsilon();
55 }
56 //___________________________________________________________________________
58 KineGeneratorWithCache("genie::QELEventGeneratorSM")
59 {
60 
61 }
62 //___________________________________________________________________________
64 KineGeneratorWithCache("genie::QELEventGeneratorSM", config)
65 {
66 
67 }
68 //___________________________________________________________________________
70 {
71 
72 }
73 //___________________________________________________________________________
75 {
76  LOG("QELEvent", pINFO) << "Generating QE event kinematics...";
77 
78  if(fGenerateUniformly) {
79  LOG("QELEvent", pNOTICE)
80  << "Generating kinematics uniformly over the allowed phase space";
81  }
82 
83  // Get the interaction and set the 'trust' bits
84  Interaction * interaction = evrec->Summary();
85  const InitialState & init_state = interaction -> InitState();
86  interaction->SetBit(kISkipProcessChk);
87  interaction->SetBit(kISkipKinematicChk);
88 
89  // Skip if no hit nucleon is set
90  if(! evrec->HitNucleon())
91  {
92  LOG("QELEvent", pFATAL) << "No hit nucleon was set";
93  gAbortingInErr = true;
94  exit(1);
95  }
96 
97  // Access the target from the interaction summary
98  Target * tgt = init_state.TgtPtr();
99 
100  // Get the random number generators
102 
103  // Access cross section algorithm for running thread
105  const EventGeneratorI * evg = rtinfo->RunningThread();
106  fXSecModel = evg->CrossSectionAlg();
107 
108  // heavy nucleus is nucleus that heavier than hydrogen and deuterium
109  bool isHeavyNucleus = tgt->A()>=3;
110 
111  sm_utils->SetInteraction(interaction);
112  // phase space for heavy nucleus is different from light one
113  fkps = isHeavyNucleus?kPSQ2vfE:kPSQ2fE;
115  // Try to calculate the maximum cross-section in kinematical limits
116  // if not pre-computed already
117  double xsec_max1 = (fGenerateUniformly) ? -1 : this->MaxXSec(evrec);
118  double xsec_max2 = (fGenerateUniformly) ? -1 : (rQ2.max<fQ2Min)? 0: this->MaxXSec2(evrec);// this make correct calculation of probability
119  double vmax= isHeavyNucleus?this->MaxDiffv(evrec) : 0.;
120 
121 
122  // generate Q2, v
123  double gQ2, v, xsec;
124  unsigned int iter = 0;
125  bool accept = false;
126  TLorentzVector q;
127  while(1)
128  {
129  LOG("QELEvent", pINFO) << "Attempt #: " << iter;
130  if(iter > kRjMaxIterations)
131  {
132  LOG("QELEvent", pWARN)
133  << "Couldn't select a valid Q2 after " << iter << " iterations";
134  evrec->EventFlags()->SetBitNumber(kKineGenErr, true);
136  exception.SetReason("Couldn't select kinematics");
137  exception.SwitchOnFastForward();
138  throw exception;
139  }
140 
141  // Pick Q2 and v
142  double xsec_max = 0.;
143  double pth = rnd->RndKine().Rndm();
144  //pth < prob1/(prob1+prob2), where prob1,prob2 - probabilities to generate event in area1 (Q2<fQ2Min) and area2 (Q2>fQ2Min) which are not normalized
145  if (pth <= xsec_max1*(TMath::Min(rQ2.max, fQ2Min)-rQ2.min)/(xsec_max1*(TMath::Min(rQ2.max, fQ2Min)-rQ2.min)+xsec_max2*(rQ2.max-fQ2Min)))
146  {
147  xsec_max = xsec_max1;
148  gQ2 = (rnd->RndKine().Rndm() * (TMath::Min(rQ2.max, fQ2Min)-rQ2.min)) + rQ2.min;
149  }
150  else
151  {
152  xsec_max = xsec_max2;
153  gQ2 = (rnd->RndKine().Rndm() * (rQ2.max-fQ2Min)) + fQ2Min;
154  }
155  Range1D_t rv = sm_utils->vQES_SM_lim(gQ2);
156  // for nuclei heavier than deuterium generate energy transfer in defined energy interval
157  v = 0.;
158  if (isHeavyNucleus)
159  {
160  v = vmax * rnd->RndKine().Rndm();
161  if (v > (rv.max-rv.min)){
162  continue;
163  }
164  }
165  v += rv.min;
166 
167  Kinematics * kinematics = interaction->KinePtr();
168  kinematics->SetKV(kKVQ2, gQ2);
169  kinematics->SetKV(kKVv, v);
170  xsec = fXSecModel->XSec(interaction, fkps);
171 
172  //-- Decide whether to accept the current kinematics
173  if(!fGenerateUniformly) {
174  this->AssertXSecLimits(interaction, xsec, xsec_max);
175 
176  double t = xsec_max * rnd->RndKine().Rndm();
177 
178 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
179  LOG("QELEvent", pDEBUG)
180  << "xsec= " << xsec << ", J= " << J << ", Rnd= " << t;
181 #endif
182  accept = (t < xsec);
183  }
184  else {
185  accept = (xsec>0);
186  }
187 
188  // If the generated kinematics are accepted, finish-up module's job
189  if(accept)
190  {
191 
192  break;
193  }
194  iter++;
195  }
196 
197 
198  // Z-frame is frame where momentum transfer is (v,0,0,qv)
199  double qv = TMath::Sqrt(v*v+gQ2);
200  TLorentzVector transferMom(0, 0, qv, v);
201 
202  Range1D_t rkF = sm_utils->kFQES_SM_lim(gQ2, v);
203  double kF = (rnd->RndKine().Rndm() * (rkF.max-rkF.min)) + rkF.min;
204 
205  // Momentum of initial neutrino in LAB frame
206  TLorentzVector * tempTLV = evrec->Probe()->GetP4();
207  TLorentzVector neutrinoMom = *tempTLV;
208  delete tempTLV;
209 
210  // define all angles in Z frame
211  double theta = neutrinoMom.Vect().Theta();
212  double phi = neutrinoMom.Vect().Phi();
213  double theta_k = sm_utils-> GetTheta_k(v, qv);
214  double costheta_k = TMath::Cos(theta_k);
215  double sintheta_k = TMath::Sin(theta_k);
216  double E_p; //energy of initial nucleon
217  double theta_p = sm_utils-> GetTheta_p(kF, v, qv, E_p);
218  double costheta_p = TMath::Cos(theta_p);
219  double sintheta_p = TMath::Sin(theta_p);
220  double fi_p = 2 * TMath::Pi() * rnd->RndGen().Rndm();
221  double cosfi_p = TMath::Cos(fi_p);
222  double sinfi_p = TMath::Sin(fi_p);
223  double psi = 2 * TMath::Pi() * rnd->RndGen().Rndm();
224 
225  // 4-momentum of initial neutrino in Z-frame
226  TLorentzVector neutrinoMomZ(neutrinoMom.P()*sintheta_k, 0, neutrinoMom.P()*costheta_k, neutrinoMom.E());
227 
228  // 4-momentum of final lepton in Z-frame
229  TLorentzVector outLeptonMom = neutrinoMomZ-transferMom;
230 
231  // 4-momentum of initial nucleon in Z-frame
232  TLorentzVector inNucleonMom(kF*sintheta_p*cosfi_p, kF*sintheta_p*sinfi_p, kF*costheta_p, E_p);
233 
234  // 4-momentum of final nucleon in Z-frame
235  TLorentzVector outNucleonMom = inNucleonMom+transferMom;
236 
237 
238  // Rotate all vectors from Z frame to LAB frame
239  TVector3 yvec (0.0, 1.0, 0.0);
240  TVector3 zvec (0.0, 0.0, 1.0);
241  TVector3 unit_nudir = neutrinoMom.Vect().Unit();
242 
243  outLeptonMom.Rotate(theta-theta_k, yvec);
244  outLeptonMom.Rotate(phi, zvec);
245 
246  inNucleonMom.Rotate(theta-theta_k, yvec);
247  inNucleonMom.Rotate(phi, zvec);
248 
249  outNucleonMom.Rotate(theta-theta_k, yvec);
250  outNucleonMom.Rotate(phi, zvec);
251 
252  outLeptonMom.Rotate(psi, unit_nudir);
253  inNucleonMom.Rotate(psi, unit_nudir);
254  outNucleonMom.Rotate(psi, unit_nudir);
255 
256  // set 4-momentum of struck nucleon
257  TLorentzVector * p4 = tgt->HitNucP4Ptr();
258  p4->SetPx( inNucleonMom.Px() );
259  p4->SetPy( inNucleonMom.Py() );
260  p4->SetPz( inNucleonMom.Pz() );
261  p4->SetE ( inNucleonMom.E() );
262 
263  int rpdgc = interaction->RecoilNucleonPdg();
264  assert(rpdgc);
265  double gW = PDGLibrary::Instance()->Find(rpdgc)->Mass();
266  LOG("QELEvent", pNOTICE) << "Selected: W = "<< gW;
267  double M = init_state.Tgt().HitNucP4().M();
268  double E = init_state.ProbeE(kRfHitNucRest);
269 
270  // (W,Q2) -> (x,y)
271  double gx=0, gy=0;
272  kinematics::WQ2toXY(E,M,gW,gQ2,gx,gy);
273 
274  // lock selected kinematics & clear running values
275  interaction->KinePtr()->SetQ2(gQ2, true);
276  interaction->KinePtr()->SetW (gW, true);
277  interaction->KinePtr()->Setx (gx, true);
278  interaction->KinePtr()->Sety (gy, true);
279  interaction->KinePtr()->ClearRunningValues();
280 
281  // set the cross section for the selected kinematics
282  evrec->SetDiffXSec(xsec,fkps);
283  if(fGenerateUniformly) {
284  double vol = sm_utils->PhaseSpaceVolume(fkps);
285  double totxsec = evrec->XSec();
286  double wght = (vol/totxsec)*xsec;
287  LOG("QELEvent", pNOTICE) << "Kinematics wght = "<< wght;
288 
289  // apply computed weight to the current event weight
290  wght *= evrec->Weight();
291  LOG("QELEvent", pNOTICE) << "Current event wght = " << wght;
292  evrec->SetWeight(wght);
293  }
294  TLorentzVector x4l(*(evrec->Probe())->X4());
295 
296  evrec->AddParticle(interaction->FSPrimLeptonPdg(), kIStStableFinalState, evrec->ProbePosition(),-1,-1,-1, outLeptonMom, x4l);
297 
298  GHepStatus_t ist;
300  ist = kIStStableFinalState;
301  else
303 
304  GHepParticle outNucleon(interaction->RecoilNucleonPdg(), ist, evrec->HitNucleonPosition(),-1,-1,-1, outNucleonMom , x4l);
305  evrec->AddParticle(outNucleon);
306 
307  // Store struck nucleon momentum
308  GHepParticle * nucleon = evrec->HitNucleon();
309  LOG("QELEvent",pNOTICE) << "pn: " << inNucleonMom.X() << ", " <<inNucleonMom.Y() << ", " <<inNucleonMom.Z() << ", " <<inNucleonMom.E();
310  nucleon->SetMomentum(inNucleonMom);
312 
313  // skip if not a nuclear target
314  if(evrec->Summary()->InitState().Tgt().IsNucleus())
315  // add a recoiled nucleus remnant
316  this->AddTargetNucleusRemnant(evrec);
317 
318 
319  LOG("QELEvent", pINFO) << "Done generating QE event kinematics!";
320 }
321 //___________________________________________________________________________
323 {
324 // add the remnant nuclear target at the GHEP record
325 
326  LOG("QELEvent", pINFO) << "Adding final state nucleus";
327 
328  double Px = 0;
329  double Py = 0;
330  double Pz = 0;
331  double E = 0;
332 
333  GHepParticle * nucleus = evrec->TargetNucleus();
334  int A = nucleus->A();
335  int Z = nucleus->Z();
336 
337  int fd = nucleus->FirstDaughter();
338  int ld = nucleus->LastDaughter();
339 
340  for(int id = fd; id <= ld; id++) {
341 
342  // compute A,Z for final state nucleus & get its PDG code and its mass
343  GHepParticle * particle = evrec->Particle(id);
344  assert(particle);
345  int pdgc = particle->Pdg();
346  bool is_p = pdg::IsProton (pdgc);
347  bool is_n = pdg::IsNeutron(pdgc);
348 
349  if (is_p) Z--;
350  if (is_p || is_n) A--;
351 
352  Px += particle->Px();
353  Py += particle->Py();
354  Pz += particle->Pz();
355  E += particle->E();
356 
357  }//daughters
358 
359  TParticlePDG * remn = 0;
360  int ipdgc = pdg::IonPdgCode(A, Z);
361  remn = PDGLibrary::Instance()->Find(ipdgc);
362  if(!remn) {
363  LOG("HadronicVtx", pFATAL)
364  << "No particle with [A = " << A << ", Z = " << Z
365  << ", pdgc = " << ipdgc << "] in PDGLibrary!";
366  assert(remn);
367  }
368 
369  double Mi = nucleus->Mass();
370  Px *= -1;
371  Py *= -1;
372  Pz *= -1;
373  E = Mi-E;
374 
375  // Add the nucleus to the event record
376  LOG("QELEvent", pINFO)
377  << "Adding nucleus [A = " << A << ", Z = " << Z
378  << ", pdgc = " << ipdgc << "]";
379 
380  int imom = evrec->TargetNucleusPosition();
381  evrec->AddParticle(
382  ipdgc,kIStStableFinalState, imom,-1,-1,-1, Px,Py,Pz,E, 0,0,0,0);
383 
384  LOG("QELEvent", pINFO) << "Done";
385  LOG("QELEvent", pINFO) << *evrec;
386 }
387 //___________________________________________________________________________
389 {
390  Algorithm::Configure(config);
391  this->LoadConfig();
392 }
393 //____________________________________________________________________________
395 {
396  Algorithm::Configure(config);
397 
398  Registry r( "QELEventGeneratorSM_specific", false ) ;
399  r.Set( "sm_utils_algo", RgAlg("genie::SmithMonizUtils","Default") ) ;
400 
402 
403  this->LoadConfig();
404 }
405 //____________________________________________________________________________
407 {
408 
409  // Safety factor for the maximum differential cross section
410  GetParamDef( "MaxXSec-SafetyFactor", fSafetyFactor, 1.2) ;
411  GetParamDef( "MaxDiffv-SafetyFactor",fSafetyFacor_nu, 1.1);
412 
413  // Minimum energy for which max xsec would be cached, forcing explicit
414  // calculation for lower eneries
415  GetParamDef( "Cache-MinEnergy", fEMin, 1.00) ;
416  GetParamDef( "Threshold-Q2", fQ2Min, 2.00);
417 
418  // Maximum allowed fractional cross section deviation from maxim cross
419  // section used in rejection method
420  GetParamDef( "MaxXSec-DiffTolerance", fMaxXSecDiffTolerance, 999999.);
422 
423  //-- Generate kinematics uniformly over allowed phase space and compute
424  // an event weight?
425  GetParamDef( "UniformOverPhaseSpace", fGenerateUniformly, false);
426 
427  //-- Generate kinematics uniformly over allowed phase space and compute
428  // an event weight?
429  GetParamDef( "IsNucleonInNucleus", fGenerateNucleonInNucleus, true);
430 
431  sm_utils = const_cast<genie::SmithMonizUtils *>(dynamic_cast<const genie::SmithMonizUtils *>( this -> SubAlg("sm_utils_algo") ) ) ;
432 }
433 //____________________________________________________________________________
434 double QELEventGeneratorSM::ComputeMaxXSec(const Interaction * interaction) const
435 {
436  double xsec_max = -1;
437  const int N_Q2 = 8;
438  const int N_v = 8;
439 
441  const double logQ2min = TMath::Log(TMath::Max(rQ2.min, eps));
442  const double logQ2max = TMath::Log(TMath::Min(rQ2.max, fQ2Min));
443 
444  double tmp_xsec_max = -1;
445  // Now scan through kinematical variables Q2,v
446  for (int Q2_n=0; Q2_n < N_Q2; Q2_n++)
447  { // Scan around Q2
448  double Q2 = TMath::Exp(Q2_n * (logQ2max-logQ2min)/N_Q2 + logQ2min);
449  Range1D_t rv = sm_utils->vQES_SM_lim(Q2);
450  const double logvmin = TMath::Log(TMath::Max(rv.min, eps));
451  const double logvmax = TMath::Log(TMath::Max(rv.max, TMath::Max(rv.min, eps)));
452  for (int v_n=0; v_n < N_v; v_n++)
453  { // Scan around v
454  double v = TMath::Exp(v_n * (logvmax-logvmin)/N_v + logvmin);
455  Kinematics * kinematics = interaction->KinePtr();
456  kinematics->SetKV(kKVQ2, Q2);
457  kinematics->SetKV(kKVv, v);
458  // Compute the QE cross section for the current kinematics
459  double xs = fXSecModel->XSec(interaction, fkps);
460  if (xs > tmp_xsec_max)
461  tmp_xsec_max = xs;
462  } // Done with v scan
463  }// Done with Q2 scan
464 
465 
466  xsec_max = tmp_xsec_max;
467  // Apply safety factor, since value retrieved from the cache might
468  // correspond to a slightly different value
469  xsec_max *= fSafetyFactor;
470  return xsec_max;
471 
472 }
473 //___________________________________________________________________________
474 double QELEventGeneratorSM::ComputeMaxXSec2(const Interaction * interaction) const
475 {
476  double xsec_max = -1;
477  const int N_Q2 = 8;
478  const int N_v = 8;
479 
481  if (rQ2.max<fQ2Min) return xsec_max;
482  const double logQ2min = TMath::Log(fQ2Min);
483  const double logQ2max = TMath::Log(rQ2.max);
484 
485  double tmp_xsec_max = -1;
486  // Now scan through kinematical variables Q2,v
487  for (int Q2_n=0; Q2_n < N_Q2; Q2_n++)
488  { // Scan around Q2
489  double Q2 = TMath::Exp(Q2_n * (logQ2max-logQ2min)/N_Q2 + logQ2min);
490  Range1D_t rv = sm_utils->vQES_SM_lim(Q2);
491  const double logvmin = TMath::Log(TMath::Max(rv.min, eps));
492  const double logvmax = TMath::Log(TMath::Max(rv.max, TMath::Max(rv.min, eps)));
493  for (int v_n=0; v_n < N_v; v_n++)
494  { // Scan around v
495  double v = TMath::Exp(v_n * (logvmax-logvmin)/N_v + logvmin);
496  Kinematics * kinematics = interaction->KinePtr();
497  kinematics->SetKV(kKVQ2, Q2);
498  kinematics->SetKV(kKVv, v);
499  // Compute the QE cross section for the current kinematics
500  double xs = fXSecModel->XSec(interaction, fkps);
501  if (xs > tmp_xsec_max)
502  tmp_xsec_max = xs;
503  } // Done with v scan
504  }// Done with Q2 scan
505 
506 
507  xsec_max = tmp_xsec_max;
508  // Apply safety factor, since value retrieved from the cache might
509  // correspond to a slightly different value
510  xsec_max *= fSafetyFactor;
511  return xsec_max;
512 
513 }
514 //___________________________________________________________________________
515 double QELEventGeneratorSM::MaxXSec2(GHepRecord * event_rec) const
516 {
517  LOG("Kinematics", pINFO)
518  << "Getting max. differential xsec for the rejection method";
519 
520  double xsec_max = -1;
521  Interaction * interaction = event_rec->Summary();
522 
523  LOG("Kinematics", pINFO)
524  << "Attempting to find a cached max{dxsec/dK} value";
525  xsec_max = this->FindMaxXSec2(interaction);
526  if(xsec_max>0) return xsec_max;
527 
528  LOG("Kinematics", pINFO)
529  << "Attempting to compute the max{dxsec/dK} value";
530  xsec_max = this->ComputeMaxXSec2(interaction);
531  if(xsec_max>0) {
532  LOG("Kinematics", pINFO) << "max{dxsec/dK} = " << xsec_max;
533  this->CacheMaxXSec2(interaction, xsec_max);
534  return xsec_max;
535  }
536 
537  LOG("Kinematics", pNOTICE)
538  << "Can not generate event kinematics {K} (max_xsec({K};E)<=0)";
539  // xsec for selected kinematics = 0
540  event_rec->SetDiffXSec(0,kPSNull);
541  // switch on error flag
542  event_rec->EventFlags()->SetBitNumber(kKineGenErr, true);
543  // reset 'trust' bits
544  interaction->ResetBit(kISkipProcessChk);
545  interaction->ResetBit(kISkipKinematicChk);
546  // throw exception
548  exception.SetReason("kinematics generation: max_xsec({K};E)<=0");
549  exception.SwitchOnFastForward();
550  throw exception;
551 
552  return 0;
553 }
554 //___________________________________________________________________________
556  const Interaction * interaction) const
557 {
558 // Find a cached max xsec for the specified xsec algorithm & interaction and
559 // close to the specified energy
560 
561  // get neutrino energy
562  double E = this->Energy(interaction);
563  LOG("Kinematics", pINFO) << "E = " << E;
564 
565  if(E < fEMin) {
566  LOG("Kinematics", pINFO)
567  << "Below minimum energy - Forcing explicit calculation";
568  return -1.;
569  }
570 
571  // access the the cache branch
572  CacheBranchFx * cb = this->AccessCacheBranch2(interaction);
573 
574  // if there are enough points stored in the cache buffer to build a
575  // spline, then intepolate
576  if( cb->Spl() ) {
577  if( E >= cb->Spl()->XMin() && E <= cb->Spl()->XMax()) {
578  double spl_max_xsec = cb->Spl()->Evaluate(E);
579  LOG("Kinematics", pINFO)
580  << "\nInterpolated: max xsec (E=" << E << ") = " << spl_max_xsec;
581  return spl_max_xsec;
582  }
583  LOG("Kinematics", pINFO)
584  << "Outside spline boundaries - Forcing explicit calculation";
585  return -1.;
586  }
587 
588  // if there are not enough points at the cache buffer to have a spline,
589  // look whether there is another point that is sufficiently close
590  double dE = TMath::Min(0.25, 0.05*E);
591  const map<double,double> & fmap = cb->Map();
592  map<double,double>::const_iterator iter = fmap.lower_bound(E);
593  if(iter != fmap.end()) {
594  if(TMath::Abs(E - iter->first) < dE) return iter->second;
595  }
596 
597  return -1;
598 
599 }
600 //___________________________________________________________________________
602  const Interaction * interaction, double max_xsec) const
603 {
604  LOG("Kinematics", pINFO)
605  << "Adding the computed max{dxsec/dK} value to cache";
606  CacheBranchFx * cb = this->AccessCacheBranch2(interaction);
607 
608  double E = this->Energy(interaction);
609  if(max_xsec>0) cb->AddValues(E,max_xsec);
610 
611  if(! cb->Spl() ) {
612  if( cb->Map().size() > 40 ) cb->CreateSpline();
613  }
614 
615  if( cb->Spl() ) {
616  if( E < cb->Spl()->XMin() || E > cb->Spl()->XMax() ) {
617  cb->CreateSpline();
618  }
619  }
620 }
621 //___________________________________________________________________________
623  const Interaction * interaction) const
624 {
625 // Returns the cache branch for this algorithm and this interaction. If no
626 // branch is found then one is created.
627 
628  Cache * cache = Cache::Instance();
629 
630  // build the cache branch key as: namespace::algorithm/config/interaction
631  string algkey = this->Id().Key();
632  string intkey = interaction->AsString();
633  string key = cache->CacheBranchKey(algkey, intkey, "2nd");
634 
635  CacheBranchFx * cache_branch =
636  dynamic_cast<CacheBranchFx *> (cache->FindCacheBranch(key));
637  if(!cache_branch) {
638  //-- create the cache branch at the first pass
639  LOG("Kinematics", pINFO) << "No Max d^nXSec/d{K}^n cache branch found";
640  LOG("Kinematics", pINFO) << "Creating cache branch - key = " << key;
641 
642  cache_branch = new CacheBranchFx("max[d^nXSec/d^n{K}] over phase space");
643  cache->AddCacheBranch(key, cache_branch);
644  }
645  assert(cache_branch);
646 
647  return cache_branch;
648 }
649 //___________________________________________________________________________
650 double QELEventGeneratorSM::ComputeMaxDiffv(const Interaction * /* interaction */) const
651 {
652  double max_diffv = -1;
653  const int N_Q2 = 10;
654 
656  for (int Q2_n = 0; Q2_n<N_Q2; Q2_n++) // Scan around Q2
657  {
658  double Q2 = rQ2.min + 1.*Q2_n * (rQ2.max-rQ2.min)/N_Q2;
659  Range1D_t rv = sm_utils->vQES_SM_lim(Q2);
660  if (rv.max-rv.min > max_diffv)
661  max_diffv = rv.max-rv.min;
662  } // Done with Q2 scan
663  max_diffv *= fSafetyFactor;
664  return max_diffv;
665 
666 }
667 //___________________________________________________________________________
668 double QELEventGeneratorSM::MaxDiffv(GHepRecord * event_rec) const
669 {
670  LOG("Kinematics", pINFO)
671  << "Getting max. vmax(Q2)-vmin(Q2) for the rejection method";
672 
673  double max_diffv = -1;
674  Interaction * interaction = event_rec->Summary();
675 
676  LOG("Kinematics", pINFO)
677  << "Attempting to find a cached max{vmax(Q2)-vmin(Q2)} value";
678  max_diffv = this->FindMaxDiffv(interaction);
679  if(max_diffv>0) return max_diffv;
680 
681  LOG("Kinematics", pINFO)
682  << "Attempting to compute the max{vmax(Q2)-vmin(Q2)} value";
683  max_diffv = this->ComputeMaxDiffv(interaction);
684  if(max_diffv>0) {
685  LOG("Kinematics", pINFO) << "max{vmax(Q2)-vmin(Q2)} = " << max_diffv;
686  this->CacheMaxDiffv(interaction, max_diffv);
687  return max_diffv;
688  }
689 
690  LOG("Kinematics", pNOTICE)
691  << "Can not generate event kinematics (max{vmax(Q2)-vmin(Q2);E}<=0)";
692  // xsec for selected kinematics = 0
693  event_rec->SetDiffXSec(0,kPSNull);
694  // switch on error flag
695  event_rec->EventFlags()->SetBitNumber(kKineGenErr, true);
696  // reset 'trust' bits
697  interaction->ResetBit(kISkipProcessChk);
698  interaction->ResetBit(kISkipKinematicChk);
699  // throw exception
701  exception.SetReason("kinematics generation: max{vmax(Q2)-vmin(Q2);E}<=0");
702  exception.SwitchOnFastForward();
703  throw exception;
704 
705  return 0;
706 }
707 //___________________________________________________________________________
709  const Interaction * interaction) const
710 {
711 // Find a cached maximum of vmax(Q2)-vmin(Q2) for xsec algorithm & interaction and
712 // close to the specified energy
713 
714  // get neutrino energy
715  double E = this->Energy(interaction);
716  LOG("Kinematics", pINFO) << "E = " << E;
717 
718  if(E < fEMin) {
719  LOG("Kinematics", pINFO)
720  << "Below minimum energy - Forcing explicit calculation";
721  return -1.;
722  }
723 
724  // access the the cache branch
725  CacheBranchFx * cb = this->AccessCacheBranchDiffv(interaction);
726 
727  // if there are enough points stored in the cache buffer to build a
728  // spline, then intepolate
729  if( cb->Spl() ) {
730  if( E >= cb->Spl()->XMin() && E <= cb->Spl()->XMax()) {
731  double spl_maxdiffv = cb->Spl()->Evaluate(E);
732  LOG("Kinematics", pINFO)
733  << "\nInterpolated: max vmax(Q2)-vmin(Q2) (E=" << E << ") = " << spl_maxdiffv;
734  return spl_maxdiffv;
735  }
736  LOG("Kinematics", pINFO)
737  << "Outside spline boundaries - Forcing explicit calculation";
738  return -1.;
739  }
740 
741  // if there are not enough points at the cache buffer to have a spline,
742  // look whether there is another point that is sufficiently close
743  double dE = TMath::Min(0.25, 0.05*E);
744  const map<double,double> & fmap = cb->Map();
745  map<double,double>::const_iterator iter = fmap.lower_bound(E);
746  if(iter != fmap.end()) {
747  if(TMath::Abs(E - iter->first) < dE) return iter->second;
748  }
749 
750  return -1;
751 
752 }
753 //___________________________________________________________________________
755  const Interaction * interaction, double max_diffv) const
756 {
757  LOG("Kinematics", pINFO)
758  << "Adding the computed max{vmax(Q2)-vmin(Q2)} value to cache";
759  CacheBranchFx * cb = this->AccessCacheBranchDiffv(interaction);
760 
761  double E = this->Energy(interaction);
762  if(max_diffv>0) cb->AddValues(E,max_diffv);
763 
764  if(! cb->Spl() ) {
765  if( cb->Map().size() > 40 ) cb->CreateSpline();
766  }
767 
768  if( cb->Spl() ) {
769  if( E < cb->Spl()->XMin() || E > cb->Spl()->XMax() ) {
770  cb->CreateSpline();
771  }
772  }
773 }
774 //___________________________________________________________________________
776  const Interaction * interaction) const
777 {
778 // Returns the cache branch for this algorithm and this interaction. If no
779 // branch is found then one is created.
780 
781  Cache * cache = Cache::Instance();
782 
783  // build the cache branch key as: namespace::algorithm/config/interaction
784  string algkey = this->Id().Key();
785  string intkey = interaction->AsString();
786  string key = cache->CacheBranchKey(algkey, intkey, "diffv");
787 
788  CacheBranchFx * cache_branch =
789  dynamic_cast<CacheBranchFx *> (cache->FindCacheBranch(key));
790  if(!cache_branch) {
791  //-- create the cache branch at the first pass
792  LOG("Kinematics", pINFO) << "No Max vmax(Q2)-vmin(Q2) cache branch found";
793  LOG("Kinematics", pINFO) << "Creating cache branch - key = " << key;
794 
795  cache_branch = new CacheBranchFx("max[vmax(Q2)-vmin(Q2)] over phase space");
796  cache->AddCacheBranch(key, cache_branch);
797  }
798  assert(cache_branch);
799 
800  return cache_branch;
801 }
802 //___________________________________________________________________________
void SetInteraction(const Interaction *i)
Range1D_t Q2QES_SM_lim(void) const
virtual double MaxXSec(GHepRecord *evrec) const
int Z(void) const
virtual GHepParticle * Particle(int position) const
Definition: GHepRecord.cxx:148
virtual void SetWeight(double wght)
Definition: GHepRecord.h:131
double ComputeMaxDiffv(const Interaction *in) const
double fQ2Min
Q2-threshold for seeking the second maximum.
bool fGenerateUniformly
uniform over allowed phase space + event weight?
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:141
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
double E(void) const
Get energy.
Definition: GHepParticle.h:92
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
virtual Interaction * Summary(void) const
Definition: GHepRecord.cxx:135
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
void SetQ2(double Q2, bool selected=false)
Definition: Kinematics.cxx:265
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
double fSafetyFactor
maxxsec -> maxxsec * safety_factor
int FirstDaughter(void) const
Definition: GHepParticle.h:69
virtual double XSec(const Interaction *i, KinePhaseSpace_t k=kPSfE) const =0
Compute the cross section for the input interaction.
int RecoilNucleonPdg(void) const
recoil nucleon pdg
int A(void) const
Definition: Target.h:71
A simple [min,max] interval for doubles.
Definition: Range1.h:43
virtual void AssertXSecLimits(const Interaction *in, double xsec, double xsec_max) const
enum genie::EGHepStatus GHepStatus_t
#define pFATAL
Definition: Messenger.h:57
bool IsNucleus(void) const
Definition: Target.cxx:289
double fMaxXSecDiffTolerance
max{100*(xsec-maxxsec)/.5*(xsec+maxxsec)} if xsec>maxxsec
double MaxDiffv(GHepRecord *evrec) const
Defines the EventGeneratorI interface.
Generated/set kinematical variables for an event.
Definition: Kinematics.h:40
virtual int HitNucleonPosition(void) const
Definition: GHepRecord.cxx:454
virtual double Weight(void) const
Definition: GHepRecord.h:125
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
double ComputeMaxXSec(const Interaction *in) const
void SetMomentum(const TLorentzVector &p4)
Definition: config.py:1
double Mass(void) const
Mass that corresponds to the PDG code.
double Evaluate(double x) const
Definition: Spline.cxx:362
Range1D_t vQES_SM_lim(double Q2) const
Abstract class. Provides a data caching mechanism for for concrete implementations of the EventRecord...
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
double Pz(void) const
Get Pz.
Definition: GHepParticle.h:91
double dE
string AsString(void) const
int FSPrimLeptonPdg(void) const
final state primary lepton pdg
double GetBindingEnergy(void) const
double Px(void) const
Get Px.
Definition: GHepParticle.h:89
void AddCacheBranch(string key, CacheBranchI *branch)
Definition: Cache.cxx:97
virtual GHepParticle * Probe(void) const
Definition: GHepRecord.cxx:321
CacheBranchFx * AccessCacheBranchDiffv(const Interaction *in) const
Float_t Z
Definition: plot.C:38
int Pdg(void) const
Definition: GHepParticle.h:64
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:304
Summary information for an interaction.
Definition: Interaction.h:56
void CacheMaxXSec2(const Interaction *in, double xsec) const
void AddValues(double x, double y)
const TLorentzVector & HitNucP4(void) const
Definition: Target.h:92
Range1D_t kFQES_SM_lim(double nu, double Q2) const
int LastDaughter(void) const
Definition: GHepParticle.h:70
An exception thrown by EventRecordVisitorI when the normal processing sequence has to be disrupted (f...
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:299
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
Spline * Spl(void) const
Definition: CacheBranchFx.h:48
Float_t E
Definition: plot.C:20
double ComputeMaxXSec2(const Interaction *in) const
void Configure(const Registry &config)
string CacheBranchKey(string k0, string k1="", string k2="") const
Definition: Cache.cxx:102
TLorentzVector * GetP4(void) const
void WQ2toXY(double Ev, double M, double W, double Q2, double &x, double &y)
Definition: KineUtils.cxx:1046
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:41
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
void ProcessEventRecord(GHepRecord *event_rec) const
bool fGenerateNucleonInNucleus
generate struck nucleon in nucleus
double MaxXSec2(GHepRecord *evrec) const
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:330
#define pINFO
Definition: Messenger.h:63
void CacheMaxDiffv(const Interaction *in, double xsec) const
TRandom3 & RndKine(void) const
rnd number generator used by kinematics generators
Definition: RandomGen.h:51
Misc GENIE control constants.
Double_t xsec[nknots]
Definition: testXsec.C:47
#define pWARN
Definition: Messenger.h:61
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition: Interaction.h:48
virtual const XSecAlgorithmI * CrossSectionAlg(void) const =0
CacheBranchI * FindCacheBranch(string key)
finding/adding cache branches
Definition: Cache.cxx:89
GENIE Cache Memory.
Definition: Cache.h:39
void Setx(double x, bool selected=false)
Definition: Kinematics.cxx:241
void SetKV(KineVar_t kv, double value)
Definition: Kinematics.cxx:345
double XMax(void) const
Definition: Spline.h:78
void SetRemovalEnergy(double Erm)
static RunningThreadInfo * Instance(void)
static const double A
Definition: Units.h:82
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:98
double max
Definition: Range1.h:54
void SetW(double W, bool selected=false)
Definition: Kinematics.cxx:289
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:264
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
const map< double, double > & Map(void) const
Definition: CacheBranchFx.h:47
Contains auxiliary functions for Smith-Moniz model. Is a concrete implementation of the XSecAlgorit...
virtual TBits * EventFlags(void) const
Definition: GHepRecord.h:118
void AddTargetNucleusRemnant(GHepRecord *evrec) const
add a recoiled nucleus remnant
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
double PhaseSpaceVolume(KinePhaseSpace_t ps) const
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
void Sety(double y, bool selected=false)
Definition: Kinematics.cxx:253
double epsilon
virtual GHepParticle * HitNucleon(void) const
Definition: GHepRecord.cxx:350
Target * TgtPtr(void) const
Definition: InitialState.h:68
TRandom3 & RndGen(void) const
rnd number generator for generic usage
Definition: RandomGen.h:81
double FindMaxXSec2(const Interaction *in) const
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
double fEMin
min E for which maxxsec is cached - forcing explicit calc.
virtual double Energy(const Interaction *in) const
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:69
exit(0)
double FindMaxDiffv(const Interaction *in) const
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
virtual double XSec(void) const
Definition: GHepRecord.h:127
double gQ2
Definition: gtestDISSF.cxx:55
const InitialState & InitState(void) const
Definition: Interaction.h:69
virtual void AddParticle(const GHepParticle &p)
Definition: GHepRecord.cxx:535
double min
Definition: Range1.h:53
int A(void) const
#define pNOTICE
Definition: Messenger.h:62
void ClearRunningValues(void)
Definition: Kinematics.cxx:357
bool GetParamDef(const RgKey &name, T &p, const T &def) const
const Target & Tgt(void) const
Definition: InitialState.h:67
virtual int ProbePosition(void) const
Definition: GHepRecord.cxx:389
static const unsigned int kRjMaxIterations
Definition: Controls.h:27
static Cache * Instance(void)
Definition: Cache.cxx:76
A simple cache branch storing the cached data in a TNtuple.
Definition: CacheBranchFx.h:38
double XMin(void) const
Definition: Spline.h:77
const EventGeneratorI * RunningThread(void)
void kinematics()
Definition: kinematics.C:10
bool gAbortingInErr
Definition: Messenger.cxx:56
void Set(RgIMapPair entry)
Definition: Registry.cxx:282
double ProbeE(RefFrame_t rf) const
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:46
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:40
Keep info on the event generation thread currently on charge. This is used so that event generation m...
string Key(void) const
Definition: AlgId.h:47
virtual int TargetNucleusPosition(void) const
Definition: GHepRecord.cxx:406
Root of GENIE utility namespaces.
CacheBranchFx * AccessCacheBranch2(const Interaction *in) const
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
virtual void SetDiffXSec(double xsec, KinePhaseSpace_t ps)
Definition: GHepRecord.h:134
Initial State information.
Definition: InitialState.h:49
#define pDEBUG
Definition: Messenger.h:64
double Py(void) const
Get Py.
Definition: GHepParticle.h:90
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:353