GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
BaryonResonanceDecayer.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2024, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5 
6  Costas Andreopoulos <c.andreopoulos \at cern.ch>
7  University of Liverpool
8 */
9 //____________________________________________________________________________
10 #include <cmath>
11 
12 #include <cmath>
13 
14 #include <TClonesArray.h>
15 #include <TDecayChannel.h>
16 #include <TMath.h>
17 
34 #include "Math/GSLMinimizer.h"
35 #include <Math/WrappedParamFunction.h>
36 
37 using namespace genie;
38 using namespace genie::controls;
39 using namespace genie::constants;
40 //____________________________________________________________________________
42 Decayer("genie::BaryonResonanceDecayer")
43 {
44  this->Initialize();
45 }
46 //____________________________________________________________________________
48 Decayer("genie::BaryonResonanceDecayer", config)
49 {
50  this->Initialize();
51 }
52 //____________________________________________________________________________
54 {
55 
56  for ( unsigned int i = 0; i < fRParams.size() ; ++i ) {
57  delete fRParams[i] ;
58  }
59 
60 }
61 //____________________________________________________________________________
63 {
64  LOG("ResonanceDecay", pINFO)
65  << "Running resonance decayer "
66  << ((fRunBefHadroTransp) ? "*before*" : "*after*") << " FSI";
67 
68  // Loop over particles, find unstable ones and decay them
69  TObjArrayIter piter(event);
70  GHepParticle * p = 0;
71  int ipos = -1;
72 
73  while( (p = (GHepParticle *) piter.Next()) ) {
74 
75  ipos++;
76  LOG("ResonanceDecay", pDEBUG) << "Checking: " << p->Name();
77 
78  int pdg_code = p->Pdg();
79  GHepStatus_t status_code = p->Status();
80 
81  // std::cout << "Decaing particle " << ipos << " with PDG " << pdg_code << std::endl ;
82 
83  if(!this->IsHandled (pdg_code)) continue;
84  if(!this->ToBeDecayed(pdg_code, status_code)) continue;
85 
86  LOG("ResonanceDecay", pNOTICE)
87  << "Decaying unstable particle: " << p->Name();
88 
89  bool decayed = this->Decay(ipos, event);
90  if ( ! decayed ) {
91  LOG("ResonanceDecay", pWARN) << "Resonance not decayed!" ;
92  LOG("ResonanceDecay", pWARN) << "Quitting the current event generation thread" ;
93 
94  event -> EventFlags() -> SetBitNumber(kHadroSysGenErr, true);
95 
97  exception.SetReason("Resonance not decayed");
98  exception.SwitchOnFastForward();
99  throw exception;
100 
101  return ;
102  }
103 
104  } // loop over particles
105 
106  LOG("ResonanceDecay", pNOTICE)
107  << "Done finding & decaying baryon resonances";
108 }
109 //____________________________________________________________________________
111  int decay_particle_id, GHepRecord * event) const
112 {
113  // Reset previous decay weight
114  fWeight = 1.;
115 
116  // Get particle to be decayed
117  GHepParticle * decay_particle = event->Particle(decay_particle_id);
118  if( ! decay_particle) {
119  LOG("ResonanceDecay", pERROR)
120  << "Particle to be decayed not in the event record. Particle ud: " << decay_particle_id ;
121  return false;
122  }
123 
124  bool to_be_deleted ;
125 
126  // Select a decay channel
127  TDecayChannel * selected_decay_channel =
128  this->SelectDecayChannel(decay_particle_id, event, to_be_deleted ) ;
129 
130  if(!selected_decay_channel) {
131  LOG("ResonanceDecay", pERROR)
132  << "No decay channel for particle " << decay_particle_id ;
133  LOG("ResonanceDecay", pERROR)
134  << *event ;
135 
136  return false;
137  }
138 
139  // Decay the exclusive state and copy daughters in the event record
140  bool decayed = this->DecayExclusive(decay_particle_id, event, selected_decay_channel);
141 
142  if ( to_be_deleted )
143  delete selected_decay_channel ;
144 
145  if ( ! decayed ) return false ;
146 
147  // Update the event weight for each weighted particle decay
148  double weight = event->Weight() * fWeight;
149  event->SetWeight(weight);
150 
151  // Mark input particle as a 'decayed state' & add its daughter links
152  decay_particle->SetStatus(kIStDecayedState);
153 
154  return true;
155 }
156 //____________________________________________________________________________
157 TDecayChannel * BaryonResonanceDecayer::SelectDecayChannel( int decay_particle_id,
158  GHepRecord * event,
159  bool & to_be_deleted ) const
160 {
161  // Get particle to be decayed
162  GHepParticle * decay_particle = event->Particle(decay_particle_id);
163  if(!decay_particle) return 0;
164 
165  // Get the particle 4-momentum and PDG code
166  TLorentzVector decay_particle_p4 = *(decay_particle->P4());
167  int decay_particle_pdg_code = decay_particle->Pdg();
168 
169  // Find the particle in the PDG library & quit if it does not exist
170  TParticlePDG * mother =
171  PDGLibrary::Instance()->Find(decay_particle_pdg_code);
172  if(!mother) {
173  LOG("ResonanceDecay", pERROR)
174  << "\n *** The particle with PDG code = " << decay_particle_pdg_code
175  << " was not found in PDGLibrary";
176  return 0;
177  }
178  LOG("ResonanceDecay", pINFO)
179  << "Decaying a " << mother->GetName()
180  << " with P4 = " << utils::print::P4AsString(&decay_particle_p4);
181 
182  // Get the resonance mass W (generally different from the mass associated
183  // with the input PDG code, since the it is produced off the mass shell)
184  double W = decay_particle_p4.M();
185  LOG("ResonanceDecay", pINFO) << "Available mass W = " << W;
186 
187  // Get all decay channels
188  TObjArray * original_decay_list = mother->DecayList();
189 
190  unsigned int nch = original_decay_list -> GetEntries();
191  LOG("ResonanceDecay", pINFO)
192  << mother->GetName() << " has: " << nch << " decay channels";
193 
194  // Loop over the decay channels (dc) and write down the branching
195  // ratios to be used for selecting a decay channel.
196  // Since a baryon resonance can be created at W < Mres, explicitly
197  // check and inhibit decay channels for which W > final-state-mass
198 
199  bool has_evolved_brs = BaryonResonanceDecayer::HasEvolvedBRs( decay_particle_pdg_code ) ;
200 
201  TObjArray * actual_decay_list = nullptr ;
202 
203  if ( has_evolved_brs ) {
204  actual_decay_list = EvolveDeltaBR( decay_particle_pdg_code, original_decay_list, W ) ;
205  if ( ! actual_decay_list ) return nullptr ;
206  nch = actual_decay_list -> GetEntries() ;
207  to_be_deleted = true ;
208  }
209  else {
210  actual_decay_list = original_decay_list ;
211  to_be_deleted = false ;
212  }
213 
214  double BR[nch], tot_BR = 0;
215 
216  for(unsigned int ich = 0; ich < nch; ich++) {
217 
218  TDecayChannel * ch = (TDecayChannel *) actual_decay_list -> At(ich);
219 
220  double fsmass = this->FinalStateMass(ch) ;
221  if ( fsmass < W ) {
222 
223  SLOG("ResonanceDecay", pDEBUG)
224  << "Using channel: " << ich
225  << " with final state mass = " << fsmass << " GeV";
226 
227  tot_BR += ch->BranchingRatio();
228 
229  } else {
230  SLOG("ResonanceDecay", pINFO)
231  << "Suppresing channel: " << ich
232  << " with final state mass = " << fsmass << " GeV";
233  } // final state mass
234 
235  BR[ich] = tot_BR;
236  }//channel loop
237 
238  if( tot_BR <= 0. ) {
239  SLOG("ResonanceDecay", pWARN)
240  << "None of the " << nch << " decay channels is available @ W = " << W;
241  return 0;
242  }
243 
244  // Select a resonance based on the branching ratios
245  unsigned int ich = 0, sel_ich; // id of selected decay channel
246  RandomGen * rnd = RandomGen::Instance();
247  double x = tot_BR * rnd->RndDec().Rndm();
248  do {
249  sel_ich = ich;
250  } while (x > BR[ich++]);
251 
252  TDecayChannel * sel_ch = (TDecayChannel *) actual_decay_list -> At(sel_ich);
253 
254  LOG("ResonanceDecay", pINFO)
255  << "Selected " << sel_ch->NDaughters() << "-particle decay channel ("
256  << sel_ich << ") has BR = " << sel_ch->BranchingRatio();
257 
258  if ( has_evolved_brs ) {
259 
260  for ( unsigned int i = 0; i < nch; ++i) {
261  if ( sel_ich != i ) delete actual_decay_list -> At(i);
262  }
263 
264  delete actual_decay_list ;
265  }
266 
267  return sel_ch;
268 }
269 //____________________________________________________________________________
271  int decay_particle_id, GHepRecord * event, TDecayChannel * ch) const
272 {
273  // Find the particle to be decayed in the event record
274  GHepParticle * decay_particle = event->Particle(decay_particle_id);
275  if(!decay_particle) return false ;
276 
277  // Get the decayed particle 4-momentum, 4-position and PDG code
278  TLorentzVector decay_particle_p4 = *(decay_particle->P4());
279  TLorentzVector decay_particle_x4 = *(decay_particle->X4());
280  int decay_particle_pdg_code = decay_particle->Pdg();
281 
282  // Get the final state mass spectrum and the particle codes
283  // for the selected decay channel
284  unsigned int nd = ch->NDaughters();
285 
286  int pdgc[nd];
287  double mass[nd];
288 
289  for(unsigned int iparticle = 0; iparticle < nd; iparticle++) {
290 
291  int daughter_code = ch->DaughterPdgCode(iparticle);
292  TParticlePDG * daughter = PDGLibrary::Instance()->Find(daughter_code);
293  assert(daughter);
294 
295  pdgc[iparticle] = daughter_code;
296  mass[iparticle] = daughter->Mass();
297 
298  SLOG("ResonanceDecay", pINFO)
299  << "+ daughter[" << iparticle << "]: "
300  << daughter->GetName() << " (pdg-code = "
301  << pdgc[iparticle] << ", mass = " << mass[iparticle] << ")";
302  }
303 
304  // Check whether the expected channel is Delta->pion+nucleon
305  bool is_delta = (decay_particle_pdg_code == kPdgP33m1232_DeltaPP ||
306  decay_particle_pdg_code == kPdgP33m1232_DeltaP ||
307  decay_particle_pdg_code == kPdgP33m1232_Delta0);
308 
309  bool is_delta_N_Pi_decay = is_delta && this->IsPiNDecayChannel(ch);
310 
311  // Decay the resonance using an N-body phase space generator
312  // The particle will be decayed in its rest frame and then the daughters
313  // will be boosted back to the original frame.
314 
315  bool is_permitted = fPhaseSpaceGenerator.SetDecay(decay_particle_p4, nd, mass);
316  if ( ! is_permitted ) return false ;
317 
318  // Find the maximum phase space decay weight
319  // double wmax = fPhaseSpaceGenerator.GetWtMax();
320  double wmax = -1;
321  for(int i=0; i<50; i++) {
322  double w = fPhaseSpaceGenerator.Generate();
323  wmax = TMath::Max(wmax,w);
324  }
325  assert(wmax>0);
326  LOG("ResonanceDecay", pINFO)
327  << "Max phase space gen. weight for current decay: " << wmax;
328 
330  {
331  // Generating weighted decays
332  // Do a single draw of momentum 4-vectors and then stop,
333  // taking into account the weight for this particular draw
334  double w = fPhaseSpaceGenerator.Generate();
335  fWeight *= TMath::Max(w/wmax, 1.);
336  }
337  else
338  {
339  // Generating un-weighted decays
340  RandomGen * rnd = RandomGen::Instance();
341  wmax *= 2;
342  bool accept_decay=false;
343  unsigned int itry=0;
344 
345  while(!accept_decay)
346  {
347  itry++;
348  assert(itry<kMaxUnweightDecayIterations);
349 
350  double w = fPhaseSpaceGenerator.Generate();
351  double gw = wmax * rnd->RndDec().Rndm();
352 
353  if(w>wmax) {
354  LOG("ResonanceDecay", pWARN)
355  << "Current decay weight = " << w << " > wmax = " << wmax;
356  }
357  LOG("ResonanceDecay", pINFO)
358  << "Current decay weight = " << w << " / R = " << gw;
359 
360  accept_decay = (gw<=w);
361 
362  // Extra logic that applies only for Delta -> N + pi
363  if( accept_decay && is_delta_N_Pi_decay ) {
364 
365  // We don't want the decay Delta -> Pi + N to be isotropic in the Delta referemce frame
366  // as generated by the simple phase space generator.
367  // In order to sample pion distribution according to W(Theta, phi) in the Delta resonance decay,
368  // we make use of the following.
369  // Note that Theta and Phi are defined in a reference frame which depends on the whole event
370  // For each event generated from a Delta -> N + Pi event with Pi emitted at
371  // at angles Theta and Phi (in the Delta rest frame), attach a random number to it.
372  // then we calculate W(Theta, Phi).
373  // Each possible final state is used to evaluate (Theta, Phi),
374  // then a random number is thrown, if the the random number is higher than W(Theta, Phi) drop this event and go
375  // back to re-generate an event and repeat the procedure.
376  // Otherwise keep this event to the record.
377  // For efficiency reasons the maxium of the function is Q2 dependent
378 
379  // Locate the pion in the decay products
380  // at this point we already know that the pion is unique so the first pion we find is our pion
381  unsigned int pi_id = 0 ;
382 
383  for(unsigned int iparticle = 0; iparticle < nd; iparticle++) {
384 
385  if ( genie::pdg::IsPion( ch->DaughterPdgCode(iparticle) ) ) {
386  pi_id = iparticle ;
387  break ;
388  }
389  }//iparticle
390 
391  TLorentzVector * lab_pion = fPhaseSpaceGenerator.GetDecay(pi_id);
392 
393  accept_decay = AcceptPionDecay( *lab_pion, decay_particle_id, event) ;
394 
395  } //if it is a Delta -> N + Pi
396 
397  }//accept_decay
398 
399  }//fGenerateWeighted
400 
401  // A decay was generated - Copy to the event record
402 
403  // Check whether the interaction is off a nuclear target or free nucleon
404  // Depending on whether this module is run before or after the hadron
405  // transport module it would affect the daughters status code
406  GHepParticle * target_nucleus = event->TargetNucleus();
407  bool in_nucleus = (target_nucleus!=0);
408 
409  // Loop over daughter list and add corresponding GHepParticles
410  for(unsigned int id = 0; id < nd; id++) {
411 
412  int daughter_pdg_code = pdgc[id];
413 
414  TLorentzVector * daughter_p4 = fPhaseSpaceGenerator.GetDecay(id);
415  LOG("ResonanceDecay", pDEBUG)
416  << "Adding daughter particle with PDG code = " << pdgc[id]
417  << " and mass = " << mass[id] << " GeV";
418 
419  bool is_hadron = pdg::IsHadron(daughter_pdg_code);
420  bool hadron_in_nuc = (in_nucleus && is_hadron && fRunBefHadroTransp);
421 
422  GHepStatus_t daughter_status_code = (hadron_in_nuc) ?
424 
425  event->AddParticle(
426  daughter_pdg_code, daughter_status_code, decay_particle_id,-1,-1,-1,
427  *daughter_p4, decay_particle_x4);
428  }
429 
430  return true ;
431 }
432 //__________________________________________________________________________________
433 TObjArray * BaryonResonanceDecayer::EvolveDeltaBR(int dec_part_pdgc, TObjArray * decay_list, double W) const {
434 
435  unsigned int nch = decay_list -> GetEntries();
436 
437  std::vector<double> widths( nch, 0. ) ;
438  double tot = 0. ;
439 
440  TDecayChannel * temp = nullptr ;
441 
442  for ( unsigned int i = 0 ; i < nch ; ++i ) {
443 
444  temp = (TDecayChannel*) decay_list -> At(i) ;
445  tot += widths[i] = EvolveDeltaDecayWidth(dec_part_pdgc, temp, W ) ;
446 
447  }
448 
449  if ( tot <= 0. ) return nullptr ;
450 
451  TObjArray * new_list = new TObjArray() ;
452 
453  TDecayChannel * update = nullptr ;
454 
455  for ( unsigned int i = 0 ; i < nch ; ++i ) {
456 
457  if ( widths[i] <= 0. ) continue ;
458 
459  temp = (TDecayChannel*) decay_list -> At(i) ;
460 
461  unsigned int nd = temp -> NDaughters() ;
462  std::vector<Int_t> ds( nd, 0 ) ;
463  for ( unsigned int d = 0 ; d < nd; ++d ) {
464  ds[d] = temp -> DaughterPdgCode(d) ;
465  }
466 
467  update = new TDecayChannel(
468  temp -> Number(),
469  temp -> MatrixElementCode(),
470  widths[i] / tot,
471  nd,
472  & ds[0]
473  ) ;
474 
475  new_list -> Add( update ) ;
476  }
477 
478  new_list -> SetOwner(kFALSE);
479 
480  return new_list ;
481 }
482 
483 //____________________________________________________________________________
484 double BaryonResonanceDecayer::EvolveDeltaDecayWidth(int dec_part_pdgc, TDecayChannel * ch, double W) const {
485 
486  /*
487  * The decay widths of the Delta in Pions or in N gammas are not constant.
488  * They depend on the actual mass of the decaying delta (W) they need to be evolved accordingly.
489  * This method tweaks the Delta branching ratios as a function of the W and
490  * returns the proper one depending on the specific decay channel.
491  */
492 
493  // identify the decay channel
494  // The delta decays only in 3 ways
495  // Delta -> Charged Pi + N
496  // Delta -> Pi0 + N
497  // Delta -> Gamma + N
498 
499  // They have evolution as a function of W that are different if the final state has pions or not
500  // so having tagged the pion is enough for the purpose of this method.
501 
502  bool has_pion = false ;
503  int pion_id = -1 ;
504  int nucleon_id = -1 ;
505  unsigned int nd = ch -> NDaughters() ;
506  for(unsigned int i = 0 ; i < nd; ++i ) {
507  if ( genie::pdg::IsPion( ch -> DaughterPdgCode(i) ) ) {
508  has_pion = true ;
509  pion_id = i ;
510  }
511 
512  if ( genie::pdg::IsNucleon( ch -> DaughterPdgCode(i) ) ) {
513  nucleon_id = i ;
514  }
515  }
516 
517 
518  // The first and most trivial evolution of the Width as a function of W
519  // is that if W is lower then the final state mass the width collapses to 0.
520 
521  if ( W < this -> FinalStateMass( ch ) ) {
522 
523  return 0. ;
524 
525  }
526 
527  // At this point, W is high enough to assume the decay of the delta in this channel
528  //
529  // The amplitude dependencies on W scales with the momentum of the pion or the photon respectivelly
530  // following these relationships
531  //
532  // (p_pi(W))^3
533  // Ampl_pi(W) = Ampl_pi(def)x---------------
534  // (p_pi(def))^3
535  //
536  //
537  // (p_ga(W))^3 (F_ga(W))^2
538  // Ampl_ga(W) = Ampl_ga(def)x--------------- x ---------------
539  // (p_ga(def))^3 (F_ga(def))^2
540  //
541  // where the "def" stand for the nominal value of the Delta mass.
542  // - pi_* are the momentum of the gamma and of the pion coming from the decay
543  // - F_ga is the form factor
544  //
545  // So the new amplitudes are evaluated and the proper value is returned
546 
547  // Getting the delta resonance from GENIE database
548  Resonance_t res = genie::utils::res::FromPdgCode( dec_part_pdgc ) ;
549 
550  double m = genie::utils::res::Mass( res ) ;
551  double m_2 = TMath::Power(m, 2);
552 
553  double mN = genie::pdg::IsProton( ch -> DaughterPdgCode( nucleon_id ) ) ? genie::constants::kProtonMass : genie::constants::kNucleonMass ;
554  double mN_2 = TMath::Power( mN, 2);
555 
556  double W_2 = TMath::Power(W, 2);
557 
558  double scaling = 0. ;
559 
560  if ( has_pion ) {
561 
562  double mPion = TMath::Abs( ch -> DaughterPdgCode( pion_id ) ) == kPdgPiP ? genie::constants::kPionMass : genie::constants::kPi0Mass ;
563  double m_aux1= TMath::Power( mN + mPion, 2) ;
564  double m_aux2= TMath::Power( mN - mPion, 2) ;
565 
566  // momentum of the pion in the Delta reference frame
567  double pPi_W = TMath::Sqrt((W_2-m_aux1)*(W_2-m_aux2))/(2*W); // at W
568  double pPi_m = TMath::Sqrt((m_2-m_aux1)*(m_2-m_aux2))/(2*m); // at the default Delta mass
569 
570  scaling = TMath::Power( pPi_W / pPi_m , 3 ) ;
571 
572  }
573  else {
574 
575  // momentum of the photon in the Delta Reference frame = Energy of the photon
576  double Egamma_W = (W_2-mN_2)/(2*W); // at W
577  double Egamma_m = (m_2-mN_2)/(2*m); // at the default Delta mass
578 
579  // form factor of the photon production
580  double fgamma_W = 1./(TMath::Power(1+Egamma_W*Egamma_W/fFFScaling, 2));
581  double fgamma_m = 1./(TMath::Power(1+Egamma_m*Egamma_m/fFFScaling, 2));
582 
583  scaling = TMath::Power( Egamma_W / Egamma_m, 3 ) * TMath::Power( fgamma_W / fgamma_m , 2 ) ;
584  }
585 
586  // get the width of the delta and obtain the width of the decay in the channel we are evolving
587  // evaluated at the nominal mass of the delta
588  double defChWidth = ch -> BranchingRatio() * genie::utils::res::Width( res ) ;
589 
590  return defChWidth * scaling ;
591 
592 }
593 //____________________________________________________________________________
594 bool BaryonResonanceDecayer::AcceptPionDecay( TLorentzVector pion,
595  int dec_part_id,
596  const GHepRecord * event ) const {
597 
598  // This evaluate the function W(theta, phi) as a function of the emitted pion and of the status of
599  // the Delta to be decayed and the whole event
600 
601  // in its simplest form W(theta) is
602  // W(Theta) = 1 − P[ 3/2 ] x L_2(cos Theta) + P[ 1/2 ] x L_2(cos Theta)
603  // where
604  // L_2 is the second Legendre polynomial L_2(x) = (3x^2 -1)/2
605  // and P[3/2] and P[1/2] have to some up to 1.
606  // But the code has been extended to include a phi dependence
607 
608  // Get the delta 4-momentum
609  GHepParticle * decay_particle = event->Particle( dec_part_id );
610  TLorentzVector delta_p4 = *(decay_particle->P4() );
611 
612  // find incoming lepton
613  TLorentzVector in_lep_p4( * (event -> Probe()-> GetP4()) ) ;
614 
615  // find outgoing lepton
616  Interaction * interaction = event->Summary();
617  TLorentzVector out_lep_p4 = interaction->KinePtr()->FSLeptonP4() ;
618 
619  TLorentzVector q = in_lep_p4 - out_lep_p4 ;
620 
621  pion.Boost(-delta_p4.BoostVector() ); // this gives us the pion in the Delta reference frame
622  q.Boost(-delta_p4.BoostVector() ); // this gives us the transferred momentm in the Delta reference frame
623 
624  TVector3 pion_dir = pion.Vect().Unit() ;
625  TVector3 z_axis = q.Vect().Unit() ;
626 
627 
628 
629  unsigned int q2_index = 0 ;
630 
631  // find out Q2 region for values
632  // note that Q2 is a lorentz invariant so it does not matter it is evaluated in the lab frame
633  // like in this case or in the Delta reference frame
634  double Q2 = - q.Mag2() ;
635  while( q2_index < fQ2Thresholds.size() ) {
636  if ( Q2 < fQ2Thresholds[q2_index] ) ++q2_index ;
637  else break ;
638  }
639 
640  double w_function ;
641 
642  if ( fDeltaThetaOnly ) {
643 
644  double c_t = pion_dir*z_axis; // cos theta
645 
646  w_function = 1. - (fR33[q2_index] - 0.5)*(3.*c_t*c_t - 1.) ;
647 
648  }
649 
650  else {
651 
652  double x[2] ; // 0 : theta, 1 : phi
653 
654  x[0] = pion_dir.Angle( z_axis ) ; // theta
655 
656  in_lep_p4.Boost(-delta_p4.BoostVector() ) ;
657  out_lep_p4.Boost( -delta_p4.BoostVector() ) ;
658 
659  // evaluate reference frame -> define x axis
660  TVector3 y_axis = in_lep_p4.Vect().Cross( out_lep_p4.Vect() ).Unit() ;
661  TVector3 x_axis = y_axis.Cross(z_axis);
662 
663  TVector3 pion_perp( z_axis.Cross( pion_dir.Cross( z_axis ).Unit() ) ) ;
664 
665  x[1] = pion_perp.Angle( x_axis ) ; // phi
666 
667  w_function = BaryonResonanceDecayer::PionAngularDist( x, fRParams[q2_index] ) ;
668 
669  }
670 
671  if ( fW_max[q2_index] <= 0. ) {
672  const_cast<genie::BaryonResonanceDecayer*>( this ) -> fW_max[q2_index] = FindDistributionExtrema( q2_index, true ) ;
673  }
674 
675  double aidrnd = fW_max[q2_index] * RandomGen::Instance()-> RndDec().Rndm();
676 
677  return ( aidrnd <= w_function ) ;
678 
679 }
680 //____________________________________________________________________________
682 {
683  return fWeight;
684 }
685 //____________________________________________________________________________
686 double BaryonResonanceDecayer::FinalStateMass( TDecayChannel * ch ) const
687 {
688 // Computes the total mass of the final state system
689 
690  double mass = 0;
691  unsigned int nd = ch->NDaughters();
692 
693  for(unsigned int iparticle = 0; iparticle < nd; iparticle++) {
694 
695  int daughter_code = ch->DaughterPdgCode(iparticle);
696  TParticlePDG * daughter = PDGLibrary::Instance()->Find(daughter_code);
697  assert(daughter);
698 
699  double md = daughter->Mass();
700 
701  // hack to switch off channels giving rare occurences of |1114| that has
702  // no decay channels in the pdg table (08/2007)
703  if(TMath::Abs(daughter_code) == 1114) {
704  LOG("ResonanceDecay", pNOTICE)
705  << "Disabling decay channel containing resonance 1114";;
706  md = 999999999;
707  }
708  mass += md;
709  }
710  return mass;
711 }
712 //____________________________________________________________________________
713 bool BaryonResonanceDecayer::IsPiNDecayChannel( TDecayChannel * ch ) const
714 {
715  if(!ch) return false;
716 
717  unsigned int nd = ch->NDaughters();
718  if(nd != 2) return false;
719 
720  int npi = 0;
721  int nnucl = 0;
722 
723  for(unsigned int iparticle = 0; iparticle < nd; iparticle++) {
724 
725  int daughter_code = ch->DaughterPdgCode(iparticle);
726 
727  if( genie::pdg::IsPion( daughter_code ) )
728  npi++;
729 
730  if ( genie::pdg::IsNucleon(daughter_code ) )
731  nnucl++;
732 
733  }//iparticle
734 
735  if(npi == 1 && nnucl == 1) return true;
736 
737  return false;
738 }
739 //____________________________________________________________________________
741 {
742 
743 }
744 //____________________________________________________________________________
745 bool BaryonResonanceDecayer::IsHandled(int pdg_code) const
746 {
747  bool is_handled = utils::res::IsBaryonResonance(pdg_code);
748 
749  LOG("ResonanceDecay", pDEBUG)
750  << "Can decay particle with PDG code = " << pdg_code
751  << "? " << ((is_handled)? "Yes" : "No");
752 
753  return is_handled ;
754 }
755 //____________________________________________________________________________
756 void BaryonResonanceDecayer::InhibitDecay(int pdgc, TDecayChannel * dc) const
757 {
758  if(! this->IsHandled(pdgc)) return;
759  if(!dc) return;
760 
761  //
762  // Not implemented
763  //
764 }
765 //____________________________________________________________________________
766 void BaryonResonanceDecayer::UnInhibitDecay(int pdgc, TDecayChannel * dc) const
767 {
768  if(! this->IsHandled(pdgc)) return;
769  if(!dc) return;
770 
771  //
772  // Not implemented
773  //
774 }
775 //____________________________________________________________________________
776 bool BaryonResonanceDecayer::IsDelta( int dec_part_pdgc ) {
777 
778  dec_part_pdgc = abs( dec_part_pdgc ) ;
779 
780  return ( dec_part_pdgc == kPdgP33m1232_DeltaM ||
781  dec_part_pdgc == kPdgP33m1232_Delta0 ||
782  dec_part_pdgc == kPdgP33m1232_DeltaP ||
783  dec_part_pdgc == kPdgP33m1232_DeltaPP ) ;
784 }
785 //____________________________________________________________________________
786 bool BaryonResonanceDecayer::HasEvolvedBRs( int dec_part_pdgc ) {
787 
788  dec_part_pdgc = abs( dec_part_pdgc ) ;
789 
790  // the evolution of the Delta BR as a function of W is meaningful only when there are
791  // more than one decay channels.
792  // Delta++ and Delta- have only one decay channel bacause of baryon number and charge conservation
793 
794  return ( dec_part_pdgc == kPdgP33m1232_Delta0 ||
795  dec_part_pdgc == kPdgP33m1232_DeltaP ) ;
796 }
797 //____________________________________________________________________________
798 double BaryonResonanceDecayer::PionAngularDist( const double* x, const double * par ) {
799 
800  double c_t = TMath::Cos( x[0] ) ;
801  double s_t = TMath::Sin( x[0] ) ;
802 
803  double c_phi = TMath::Cos( x[1 ] );
804 
805  double theta_dep_only = 1. - ( par[0] - 0.5 )*( 3.*c_t*c_t - 1. ) ;
806  double phi_dependency = kSqrt3 *( 2.*par[1]*s_t*c_t*c_phi + par[2]*s_t*(2.*c_phi*c_phi-1.) ) ;
807 
808  return theta_dep_only - phi_dependency ;
809 
810 }
811 //____________________________________________________________________________
812 double BaryonResonanceDecayer::FindDistributionExtrema( unsigned int q2_bin, bool find_maximum ) const {
813 
814  // Choose method upon creation between:
815  // kConjugateFR, kConjugatePR, kVectorBFGS,
816  // kVectorBFGS2, kSteepestDescent
817  ROOT::Math::GSLMinimizer min( ROOT::Math::kVectorBFGS );
818 
819  min.SetMaxFunctionCalls(1000);
820  min.SetMaxIterations(1000);
821  min.SetTolerance( fMaxTolerance );
822 
823  ROOT::Math::WrappedParamFunction<ROOT::Math::FreeParamMultiFunctionPtr> f( ( find_maximum ?
826  2, 3, fRParams[q2_bin] ) ;
827 
828  // Steps are defined as the same fraction of the range of each variable
829  double step[2] = { 0.00005 * kPi, 0.00005 * 2 * kPi } ;
830 
831  min.SetFunction(f);
832 
833  do {
834 
835  // Set the free variables to be minimized!
836  // it has been observed that some initial values are making the minimization to fail.
837  // e.g. ( pi/2, pi/2 )
838  // at the same time, the absolute extrema seems very stable with respect to changes of the initial variable
839  // hence the minimization is done inside a loop where the variables a initialized randomly
840 
841  min.SetLimitedVariable( 0, "theta", kPi * RandomGen::Instance()-> RndDec().Rndm(), step[0], 0., kPi );
842  min.SetLimitedVariable( 1, "phi", 2*kPi * RandomGen::Instance()-> RndDec().Rndm(), step[1], 0., 2*kPi );
843 
844  } while( ! min.Minimize() ) ;
845 
846  const double *xs = min.X();
847 
848  double result = find_maximum ? -1 * f( xs ) : f( xs ) ;
849 
850  LOG("BaryonResonanceDecayer", pINFO) << (find_maximum ? "Maximum " : "Minimum ")
851  << "of angular distribution found in ( "
852  << xs[0] << ", " << xs[1] << " ): "
853  << result ;
854 
855  LOG("BaryonResonanceDecayer", pDEBUG) << "Minimum found in " << min.NCalls() << " calls" ;
856 
857  return result ;
858 
859 }
860 
861 
862 //____________________________________________________________________________
864 
866 
867  this -> GetParam( "FFScaling", fFFScaling ) ;
868 
869  this -> GetParamDef( "Delta-ThetaOnly", fDeltaThetaOnly, true ) ;
870 
871  this -> GetParamDef( "DeltaDecayMaximumTolerance", fMaxTolerance, 0.0005 ) ;
872 
873  bool invalid_configuration = false ;
874 
875  // load R33 parameters
876  this -> GetParamVect( "Delta-R33", fR33 ) ;
877 
878  // load Q2 thresholds if necessary
879  if ( fR33.size() > 1 ) {
880  this -> GetParamVect("Delta-Q2", fQ2Thresholds ) ;
881  }
882  else {
883  fQ2Thresholds.clear() ;
884  }
885 
886  // check if the number of Q2 matches the number of R33
887  if ( fQ2Thresholds.size() != fR33.size() -1 ) {
888  invalid_configuration = true ;
889  LOG("BaryonResonanceDecayer", pFATAL) << "Delta-Q2 and Delta-R33 have wrong sizes" ;
890  LOG("BaryonResonanceDecayer", pFATAL) << "Delta-Q2 -> " << fQ2Thresholds.size() ;
891  LOG("BaryonResonanceDecayer", pFATAL) << "Delta-R33 -> " << fR33.size() ;
892  }
893 
894  if ( fDeltaThetaOnly ) {
895 
896  // check the parameters validity
897  for ( unsigned int i = 0 ; i < fR33.size(); ++i ) {
898  if ( (fR33[i] < -0.5) || (fR33[i] > 1.) ) {
899  invalid_configuration = true ;
900  LOG("BaryonResonanceDecayer", pFATAL) << "Delta-R33[" << i << "] out of valid range [-0.5 ; 1 ]" ;
901  LOG("BaryonResonanceDecayer", pFATAL) << "Delta-R33[" << i << "] = " << fR33[i] ;
902  break ;
903  }
904  }
905 
906  // set appropriate maxima
907  fW_max.resize( fR33.size(), 0. ) ;
908  for ( unsigned int i = 0 ; i < fR33.size(); ++i ) {
909  fW_max[i] = ( fR33[i] < 0.5 ? 2. * ( 1. - fR33[i] ) : fR33[i] + 0.5 ) + fMaxTolerance ;
910  }
911 
912  } // Delta Theta Only
913 
914  else {
915 
916  // load R31 and R3m1 parameters
917  this -> GetParamVect( "Delta-R31", fR31 ) ;
918 
919  this -> GetParamVect( "Delta-R3m1", fR3m1 ) ;
920 
921  // check if they match the numbers of R33
922  if ( (fR31.size() != fR33.size()) || (fR3m1.size() != fR33.size()) ) {
923  LOG("BaryonResonanceDecayer", pFATAL) << "Delta-R31 or Delta-R3m1 sizes don't match Delta-R33" ;
924  LOG("BaryonResonanceDecayer", pFATAL) << "R31: " << fR31.size()
925  << ", R3m1: " << fR31.size()
926  << " while R33: " << fR33.size() ;
927  invalid_configuration = true ;
928  }
929 
930  for ( unsigned int i = 0; i < fRParams.size() ; ++i ) {
931  delete [] fRParams[i] ;
932  }
933  fRParams.clear() ;
934  // fill the container by Q2 bin instead of the parmaeter bin
935  for ( unsigned int i = 0 ; i < fR33.size(); ++i ) {
936  fRParams.push_back( new double[3]{ fR33[i], fR31[i], fR3m1[i] } ) ;
937  }
938 
939 
940  // check if they are physical
941  fW_max.resize( fR33.size(), 0. ) ;
942  for ( unsigned int i = 0 ; i < fR33.size(); ++i ) {
943 
944  double temp_min = FindDistributionExtrema( i, false ) ;
945  if ( temp_min < 0. ) {
946  LOG("BaryonResonanceDecayer", pFATAL) << "pion angular distribution minimum is negative for Q2 bin " << i ;
947  invalid_configuration = true ;
948  break ;
949  }
950 
951  double temp_max = FindDistributionExtrema( i, true ) ;
952  if ( temp_max <= 0. ) {
953  LOG("BaryonResonanceDecayer", pFATAL) << "pion angular distribution maximum is non positive for Q2 bin " << i ;
954  invalid_configuration = true ;
955  break ;
956  }
957 
958  fW_max[i] = temp_max + fMaxTolerance ;
959 
960  }
961 
962  }
963 
964  if ( invalid_configuration ) {
965 
966  LOG("BaryonResonanceDecayer", pFATAL)
967  << "Invalid configuration: Exiting" ;
968 
969  // From the FreeBSD Library Functions Manual
970  //
971  // EX_CONFIG (78) Something was found in an unconfigured or miscon-
972  // figured state.
973 
974  exit( 78 ) ;
975 
976  }
977 
978 }
979 
Baryon resonance decayer module.
virtual bool ToBeDecayed(int pdgc, GHepStatus_t ist) const
Definition: Decayer.cxx:51
const int kPdgP33m1232_DeltaPP
Definition: PDGCodes.h:107
bool IsPion(int pdgc)
Definition: PDGUtils.cxx:326
std::vector< double * > fRParams
bool Decay(int dec_part_id, GHepRecord *event) const
#define pERROR
Definition: Messenger.h:59
static const double kNucleonMass
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
void UnInhibitDecay(int pdgc, TDecayChannel *ch=0) const
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:78
static bool IsDelta(int dec_part_pdgc)
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
bool IsNucleon(int pdgc)
Definition: PDGUtils.cxx:346
static bool HasEvolvedBRs(int dec_part_pdgc)
#define pFATAL
Definition: Messenger.h:56
int GetParamVect(const std::string &comm_name, std::vector< T > &v, bool is_top_call=true) const
Handle to load vectors of parameters.
string P4AsString(const TLorentzVector *p)
Definition: PrintUtils.cxx:27
static const unsigned int kMaxUnweightDecayIterations
Definition: Controls.h:61
bool DecayExclusive(int dec_part_id, GHepRecord *event, TDecayChannel *ch) const
double Mass(Resonance_t res)
resonance mass (GeV)
double Width(Resonance_t res)
resonance width (GeV)
double EvolveDeltaDecayWidth(int dec_part_pdgc, TDecayChannel *ch, double W) const
enum genie::EResonance Resonance_t
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
GHepStatus_t Status(void) const
Definition: GHepParticle.h:64
const int kPdgP33m1232_DeltaP
Definition: PDGCodes.h:106
const int kPdgP33m1232_DeltaM
Definition: PDGCodes.h:104
TDecayChannel * SelectDecayChannel(int dec_part_id, GHepRecord *event, bool &to_be_deleted) const
double W(const Interaction *const i)
Definition: KineUtils.cxx:1101
int Pdg(void) const
Definition: GHepParticle.h:63
string Name(void) const
Name that corresponds to the PDG code.
bool IsHadron(int pdgc)
Definition: PDGUtils.cxx:392
Summary information for an interaction.
Definition: Interaction.h:56
double FindDistributionExtrema(unsigned int i, bool find_maximum=false) const
An exception thrown by EventRecordVisitorI when the normal processing sequence has to be disrupted (f...
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
const TLorentzVector & FSLeptonP4(void) const
Definition: Kinematics.h:65
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
bool IsPiNDecayChannel(TDecayChannel *ch) const
void ProcessEventRecord(GHepRecord *event) const
bool fRunBefHadroTransp
is invoked before or after FSI?
Definition: Decayer.h:57
const int kPdgPiP
Definition: PDGCodes.h:158
TObjArray * EvolveDeltaBR(int dec_part_pdgc, TObjArray *decay_list, double W) const
double FinalStateMass(TDecayChannel *ch) const
#define pINFO
Definition: Messenger.h:62
bool fGenerateWeighted
generate weighted or unweighted decays?
Definition: Decayer.h:56
Resonance_t FromPdgCode(int pdgc)
PDG code -&gt; resonance id.
#define pWARN
Definition: Messenger.h:60
const int kPdgP33m1232_Delta0
Definition: PDGCodes.h:105
static double MinusPionAngularDist(const double *x, const double *par)
void InhibitDecay(int pdgc, TDecayChannel *ch=0) const
virtual void LoadConfig(void)
Definition: Decayer.cxx:135
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:79
void SetStatus(GHepStatus_t s)
Definition: GHepParticle.h:126
static double PionAngularDist(const double *x, const double *par)
Base class for decayer classes. Implements common configuration, allowing users to toggle on/off flag...
Definition: Decayer.h:34
bool AcceptPionDecay(TLorentzVector lab_pion, int dec_part_id, const GHepRecord *event) const
bool IsBaryonResonance(int pdgc)
is input a baryon resonance?
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
#define pNOTICE
Definition: Messenger.h:61
bool GetParamDef(const RgKey &name, T &p, const T &def) const
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:84
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:45
static constexpr double m
Definition: Units.h:71
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
TRandom3 & RndDec(void) const
rnd number generator used by decay models
Definition: RandomGen.h:56
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
enum genie::EGHepStatus GHepStatus_t
#define pDEBUG
Definition: Messenger.h:63