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