Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4QContent.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id$
28 //
29 // ---------------- G4QContent ----------------
30 // by Mikhail Kossov, Sept 1999.
31 // class for Quasmon initiated Contents used by CHIPS Model
32 // --------------------------------------------------------------------
33 // Short description: This is the basic class of the CHIPS model. It
34 // describes the quark content of the Quasmon, which is a generalized
35 // hadronic state. All Quasmons are bags, characterized by the quark
36 // Content (QContent), but the spin is not fixed and only light (u,d,s)
37 // quarks are considered (SU(3)). The hadrons are the ground states for
38 // the corresponding quasmons. The Chipolino (G4QChipolino) or nuclear
39 // cluster are examples for another Quark Content.
40 // --------------------------------------------------------------------
41 // @@ In future total spin & c,b,t of the Hadron can be added @@ M.K.@@
42 // --------------------------------------------------------------------
43 
44 //#define debug
45 //#define pdebug
46 //#define ppdebug
47 //#define erdebug
48 
49 #include "G4QContent.hh"
50 #include <cmath>
51 using namespace std;
52 
53 // Constructors
54 
55 // Initialize by the full quark content
56 G4QContent::G4QContent(G4int d, G4int u, G4int s_value, G4int ad, G4int au, G4int as):
57  nD(d),nU(u),nS(s_value),nAD(ad),nAU(au),nAS(as)
58 {
59  // Bug report 1131 identifies conditional to have no effect.
60  // Remove it.
61  // D.H. Wright 11 November 2010
62  /*
63  if(d<0||u<0||s_value<0||ad<0||au<0||as<0)
64  {
65 #ifdef erdebug
66  G4cerr<<"***G4QContent:"<<d<<","<<u<<","<<s_value<<","<<ad<<","<<au<<","<<as<<G4endl;
67 #endif
68  if(d<0) ad-=d;
69  if(u<0) au-=u;
70  if(s_value<0) as-=s_value;
71  if(ad<0) d-=ad;
72  if(au<0) u-=au;
73  if(as<0) s_value-=as;
74  }
75  */
76 }
77 
78 // Initialize by a pair of partons
79 G4QContent::G4QContent(std::pair<G4int,G4int> PP): nD(0),nU(0),nS(0),nAD(0),nAU(0),nAS(0)
80 {
81  G4int P1=PP.first;
82  G4int P2=PP.second;
83  if(!P1 || !P2)
84  {
85  // G4cerr<<"***G4QContent::Constr(pair): Zero parton P1="<<P1<<", P2="<<P2<<G4endl;
86  // G4Exception("G4QContent::Constructor(pair):","72",FatalException,"WrongPartonPair");
87  G4ExceptionDescription ed;
88  ed << "Wrong parton pair: zero parton P1=" << P1 << ", P2=" << P2 << G4endl;
89  G4Exception("G4QContent::G4QContent()", "HAD_CHPS_0072", FatalException, ed);
90  }
91 #ifdef debug
92  G4cout<<"G4QContent::PairConstr: P1="<<P1<<", P2="<<P2<<G4endl;
93 #endif
94  G4bool suc=true;
95  G4int A1=P1;
96  if (P1 > 7) A1= P1/100;
97  else if(P1 <-7) A1=(-P1)/100;
98  else if(P1 < 0) A1= -P1;
99  G4int P11=0;
100  G4int P12=0;
101  if(A1>7)
102  {
103  P11=A1/10;
104  P12=A1%10;
105  }
106  if(P1>0)
107  {
108  if(!P11)
109  {
110  if (A1==1) ++nD;
111  else if(A1==2) ++nU;
112  else if(A1==3) ++nS;
113  else suc=false;
114  }
115  else
116  {
117  if (P11==1) ++nD;
118  else if(P11==2) ++nU;
119  else if(P11==3) ++nS;
120  else suc=false;
121  if (P12==1) ++nD;
122  else if(P12==2) ++nU;
123  else if(P12==3) ++nS;
124  else suc=false;
125  }
126  }
127  else // negative parton
128  {
129  if(!P11)
130  {
131  if (A1==1) ++nAD;
132  else if(A1==2) ++nAU;
133  else if(A1==3) ++nAS;
134  else suc=false;
135  }
136  else
137  {
138  if (P11==1) ++nAD;
139  else if(P11==2) ++nAU;
140  else if(P11==3) ++nAS;
141  else suc=false;
142  if (P12==1) ++nAD;
143  else if(P12==2) ++nAU;
144  else if(P12==3) ++nAS;
145  else suc=false;
146  }
147  }
148 #ifdef debug
149  G4cout<<"G4QContent::PCo:1:"<<nD<<","<<nU<<","<<nS<<","<<nAD<<","<<nAU<<","<<nAS<<G4endl;
150 #endif
151  G4int A2=P2;
152  if (P2 > 7) A2= P2/100;
153  else if(P2 <-7) A2=(-P2)/100;
154  else if(P2 < 0) A2= -P2;
155  G4int P21=0;
156  G4int P22=0;
157  if(A2>7)
158  {
159  P21=A2/10;
160  P22=A2%10;
161  }
162  if(P2>0)
163  {
164  if(!P21)
165  {
166  if (A2==1) ++nD;
167  else if(A2==2) ++nU;
168  else if(A2==3) ++nS;
169  else suc=false;
170  }
171  else
172  {
173  if (P21==1) ++nD;
174  else if(P21==2) ++nU;
175  else if(P21==3) ++nS;
176  else suc=false;
177  if (P22==1) ++nD;
178  else if(P22==2) ++nU;
179  else if(P22==3) ++nS;
180  else suc=false;
181  }
182  }
183  else // negative parton
184  {
185  if(!P21)
186  {
187  if (A2==1) ++nAD;
188  else if(A2==2) ++nAU;
189  else if(A2==3) ++nAS;
190  else suc=false;
191  }
192  else
193  {
194  if (P21==1) ++nAD;
195  else if(P21==2) ++nAU;
196  else if(P21==3) ++nAS;
197  else suc=false;
198  if (P22==1) ++nAD;
199  else if(P22==2) ++nAU;
200  else if(P22==3) ++nAS;
201  else suc=false;
202  }
203  }
204  if(!suc)
205  {
206  // G4cerr<<"***G4QContent::Constr(pair): Impossible partons, P1="<<P1<<",P2="<<P2<<G4endl;
207  // G4Exception("G4QContent::Constructor(pair):","72",FatalException,"ImpossibPartonPair");
208  G4ExceptionDescription ed;
209  ed << "Impossible parton pair, P1=" << P1 << ",P2=" << P2 << G4endl;
210  G4Exception("G4QContent::G4QContent()", "HAD_CHPS_0073", FatalException, ed);
211  }
212 #ifdef debug
213  G4cout<<"G4QContent::PCo:2:"<<nD<<","<<nU<<","<<nS<<","<<nAD<<","<<nAU<<","<<nAS<<G4endl;
214 #endif
215 }
216 
217 G4QContent::G4QContent(const G4QContent& right)
218 {
219  nU = right.nU;
220  nD = right.nD;
221  nS = right.nS;
222  nAU = right.nAU;
223  nAD = right.nAD;
224  nAS = right.nAS;
225 }
226 
227 G4QContent::G4QContent(G4QContent* right)
228 {
229  nU = right->nU;
230  nD = right->nD;
231  nS = right->nS;
232  nAU = right->nAU;
233  nAD = right->nAD;
234  nAS = right->nAS;
235 }
236 
237 // Assignment operator (copy stile for possible Vector extention)
238 const G4QContent& G4QContent::operator=(const G4QContent &right)
239 {
240  if(this != &right) // Beware of self assignment
241  {
242  nU = right.nU;
243  nD = right.nD;
244  nS = right.nS;
245  nAU = right.nAU;
246  nAD = right.nAD;
247  nAS = right.nAS;
248  }
249  return *this;
250 }
251 
252 // Standard output for QC {d,u,s,ad,au,as}
253 ostream& operator<<(ostream& lhs, G4QContent& rhs)
254 {
255  lhs << "{" << rhs.GetD() << "," << rhs.GetU() << "," << rhs.GetS() << ","
256  << rhs.GetAD() << "," << rhs.GetAU() << "," << rhs.GetAS() << "}";
257  return lhs;
258 }
259 
260 // Standard output for const QC {d,u,s,ad,au,as}
261 ostream& operator<<(ostream& lhs, const G4QContent& rhs)
262 {
263  lhs << "{" << rhs.GetD() << "," << rhs.GetU() << "," << rhs.GetS() << ","
264  << rhs.GetAD() << "," << rhs.GetAU() << "," << rhs.GetAS() << "}";
265  return lhs;
266 }
267 
268 // Subtract Quark Content
269 G4QContent G4QContent::operator-=(const G4QContent& rhs)
270 {
271 #ifdef debug
272  G4cout<<"G4QC::-=(const): is called:"<<G4endl;
273 #endif
274  G4int rD=rhs.nD;
275  G4int rU=rhs.nU;
276  G4int rS=rhs.nS;
277  G4int rAD=rhs.nAD;
278  G4int rAU=rhs.nAU;
279  G4int rAS=rhs.nAS;
284  if(nU<rU||nAU<rAU||nD<rD||nAD<rAD)
285  {
286  G4int dU=rU-nU;
287  G4int dAU=rAU-nAU;
288  if(dU>0||dAU>0)
289  {
290  G4int kU=dU;
291  if(kU<dAU) kU=dAU; // Get biggest difference
292  G4int mU=rU;
293  if(rAU<mU) mU=rAU; // Get a#of possible SS pairs
294  if(kU<=mU) // Total compensation
295  {
296  rU-=kU;
297  rAU-=kU;
298  rD+=kU;
299  rAD+=kU;
300  }
301  else // Partial compensation
302  {
303  rU-=mU;
304  rAU-=mU;
305  rD+=mU;
306  rAD+=mU;
307  }
308  }
309  G4int dD=rD-nD;
310  G4int dAD=rAD-nAD;
311  if(dD>0||dAD>0)
312  {
313  G4int kD=dD;
314  if(kD<dAD) kD=dAD; // Get biggest difference
315  G4int mD=rD;
316  if(rAD<mD) mD=rAD; // Get a#of possible SS pairs
317  if(kD<=mD) // Total compensation
318  {
319  rD-=kD;
320  rAD-=kD;
321  rU+=kD;
322  rAU+=kD;
323  }
324  else // Partial compensation
325  {
326  rD-=mD;
327  rAD-=mD;
328  rU+=mD;
329  rAU+=mD;
330  }
331  }
332  }
333 #ifdef debug
334  G4cout<<"G4QC::-=:comp: "<<rD<<","<<rU<<","<<rS<<","<<rAD<<","<<rAU<<","<<rAS<<G4endl;
335 #endif
336  if(rS==1 && rAS==1 && (nS<1 || nAS<1)) // Eta case, switch quark pairs (?)
337  {
338  rS =0;
339  rAS=0;
340  if(nU>rU&&nAU>rAU)
341  {
342  rU +=1;
343  rAU+=1;
344  }
345  else
346  {
347  rD +=1;
348  rAD+=1;
349  }
350  }
351  nD -= rD;
352  if (nD<0)
353  {
354  nAD -= nD;
355  nD = 0;
356  }
357  nU -= rU;
358  if (nU<0)
359  {
360  nAU -= nU;
361  nU = 0;
362  }
363  nS -= rS;
364  if (nS<0)
365  {
366  nAS -= nS;
367  nS = 0;
368  }
369  nAD -= rAD;
370  if (nAD<0)
371  {
372  nD -= nAD;
373  nAD = 0;
374  }
375  nAU -= rAU;
376  if (nAU<0)
377  {
378  nU -= nAU;
379  nAU = 0;
380  }
381  nAS -= rAS;
382  if (nAS<0)
383  {
384  nS -= nAS;
385  nAS = 0;
386  }
387  return *this;
388 }
389 
390 // Subtract Quark Content
391 G4QContent G4QContent::operator-=(G4QContent& rhs)
392 {
393 #ifdef debug
394  G4cout<<"G4QC::-=: is called:"<<G4endl;
395 #endif
396  G4int rD=rhs.nD;
397  G4int rU=rhs.nU;
398  G4int rS=rhs.nS;
399  G4int rAD=rhs.nAD;
400  G4int rAU=rhs.nAU;
401  G4int rAS=rhs.nAS;
402  G4int rQ =rD+rU+rS;
403  G4int rAQ=rAD+rAU+rAS;
404  G4int nQ =nD+nU+nS;
405  G4int nAQ=nAD+nAU+nAS;
406  if(nQ<rQ||nAQ<rAQ)
407  {
408  G4int dU=rU-nU;
409  G4int dAU=rAU-nAU;
410  if(dU>0||dAU>0)
411  {
412  G4int kU=dU;
413  if(kU<dAU) kU=dAU; // Get biggest difference
414  G4int mU=rU;
415  if(rAU<mU) mU=rAU; // Get a#of possible SS pairs
416  if(kU<=mU) // Total compensation
417  {
418  rU-=kU;
419  rAU-=kU;
420  rD+=kU;
421  rAD+=kU;
422  }
423  else // Partial compensation
424  {
425  rU-=mU;
426  rAU-=mU;
427  rD+=mU;
428  rAD+=mU;
429  }
430  }
431  G4int dD=rD-nD;
432  G4int dAD=rAD-nAD;
433  if(dD>0||dAD>0)
434  {
435  G4int kD=dD;
436  if(kD<dAD) kD=dAD; // Get biggest difference
437  G4int mD=rD;
438  if(rAD<mD) mD=rAD; // Get a#of possible SS pairs
439  if(kD<=mD) // Total compensation
440  {
441  rD-=kD;
442  rAD-=kD;
443  rU+=kD;
444  rAU+=kD;
445  }
446  else // Partial compensation
447  {
448  rD-=mD;
449  rAD-=mD;
450  rU+=mD;
451  rAU+=mD;
452  }
453  }
454  }
455  if(rS==1 && rAS==1 && (nS<1 || nAS<1)) // Eta case, switch quark pairs (?)
456  {
457  rS =0;
458  rAS=0;
459  if(nU>rU&&nAU>rAU)
460  {
461  rU +=1;
462  rAU+=1;
463  }
464  else
465  {
466  rD +=1;
467  rAD+=1;
468  }
469  }
470  nD -= rD;
471  if (nD<0)
472  {
473  nAD -= nD;
474  nD = 0;
475  }
476  nU -= rU;
477  if (nU<0)
478  {
479  nAU -= nU;
480  nU = 0;
481  }
482  nS -= rS;
483  if (nS<0)
484  {
485  nAS -= nS;
486  nS = 0;
487  }
488  nAD -= rAD;
489  if (nAD<0)
490  {
491  nD -= nAD;
492  nAD = 0;
493  }
494  nAU -= rAU;
495  if (nAU<0)
496  {
497  nU -= nAU;
498  nAU = 0;
499  }
500  nAS -= rAS;
501  if (nAS<0)
502  {
503  nS -= nAS;
504  nAS = 0;
505  }
506  return *this;
507 }
508 
509 // Overloading of QC addition
510 G4QContent operator+(const G4QContent& lhs, const G4QContent& rhs)
511 {
512  G4QContent s_value = lhs;
513  return s_value += rhs;
514 }
515 
516 // Overloading of QC subtraction
517 G4QContent operator-(const G4QContent& lhs, const G4QContent& rhs)
518 {
519  G4QContent s_value = lhs;
520  return s_value -= rhs;
521 }
522 
523 // Overloading of QC multiplication by Int
524 G4QContent operator*(const G4QContent& lhs, const G4int& rhs)
525 {
526  G4QContent s_value = lhs;
527  return s_value *= rhs;
528 }
529 
530 // Overloading of Int times QC multiplication
531 G4QContent operator*(const G4int& lhs, const G4QContent& rhs)
532 {
533  G4QContent s_value = rhs;
534  return s_value *= lhs;
535 }
536 
537 // Destructor - - - - - - -
538 G4QContent::~G4QContent() {}
539 
540 // Subtract neutral pion from Quark Content (with possible hidden strangeness)
541 G4bool G4QContent::SubtractPi0()
542 {
543 #ifdef debug
544  G4cout<<"G4QC::SubtractPi0: U="<<nU<<", AU="<<nAU<<", D="<<nD<<", AD="<<nAD<<G4endl;
545 #endif
546  G4int tot=GetTot();
547  G4int ab =GetBaryonNumber();
548  if(ab){if(tot<3*ab+2) return false;}
549  else if(tot<4) return false;
550 
551  if(nU>0 && nAU>0)
552  {
553  nU--;
554  nAU--;
555  return true;
556  }
557  else if(nD>0 && nAD>0)
558  {
559  nD--;
560  nAD--;
561  return true;
562  }
563  return false;
564 }
565 
566 // Subtract charged pion from Quark Content
567 G4bool G4QContent::SubtractPion()
568 {
569 #ifdef debug
570  G4cout<<"G4QC::SubtractPion: U="<<nU<<", AU="<<nAU<<", D="<<nD<<", AD="<<nAD<<G4endl;
571 #endif
572  G4int tot=GetTot();
573  G4int ab =GetBaryonNumber();
574  if(ab){if(tot<3*ab+2) return false;}
575  else if(tot<4) return false;
576 
577  if((nU>nAU || (ab && nU>0))&& nAD>0)
578  {
579  nU--;
580  nAD--;
581  return true;
582  }
583  else if((nAU>nU || (ab && nAU>0)) && nD>0)
584  {
585  nAU--;
586  nD--;
587  return true;
588  }
589  return false;
590 }
591 
592 // Subtract Hadron from Quark Content
593 G4bool G4QContent::SubtractHadron(G4QContent h)
594 {
595 #ifdef debug
596  G4cout<<"G4QC::SubtractHadron "<<h<<" is called for QC="<<GetThis()<<G4endl;
597 #endif
598  if(h.GetU()<=nU && h.GetD()<=nD && h.GetS()<=nS&&
599  h.GetAU()<=nAU&& h.GetAD()<=nAD&& h.GetAS()<=nAS) return true;
600  return false;
601 }
602 
603 // Subtract Kaon from Quark Content
604 G4bool G4QContent::SubtractKaon(G4double mQ)
605 {
606 #ifdef debug
607  G4cout<<"G4QC::SubtractKaon is called: QC="<<GetThis()<<G4endl;
608 #endif
609  if(mQ<640.) return false;
610  G4int tot=GetTot();
611  G4int ab =GetBaryonNumber();
612  if(ab){if(tot<3*ab+2) return false;}
613  else if(tot<4) return false;
614 
615  if((nS>nAS || (ab && nS>0)) && (nAD>0 || nAU>0))
616  {
617  nS--;
618  if (nAU>0) nAU--;
619  else nAD--;
620  return true;
621  }
622  else if((nAS>nS || (ab && nAS>0)) && (nD>0 || nU>0))
623  {
624  nAS--;
625  if (nU>0) nU--;
626  else nD--;
627  return true;
628  }
629 #ifdef debug
630  G4cout<<"G4QCont::SubtractKaon Can't SubtractKaon: QC="<<GetThis()<<G4endl;
631 #endif
632  return false;
633 }
634 
635 // Split any hadronic system in two hadrons
636 G4QContent G4QContent::SplitChipo (G4double mQ)
637 {
638  G4QContent Pi(0,1,0,1,0,0);
639  if (nU>0&&nAU>0) Pi=G4QContent(0,1,0,0,1,0);
640  else if (nD>0&&nAD>0) Pi=G4QContent(1,0,0,1,0,0);
641  else if (nD>=nU&&nAU>=nAD) Pi=G4QContent(1,0,0,0,1,0);
642  G4int bn=GetBaryonNumber();
643  G4int b =abs(bn);
644  if(!b && mQ<545.&&nS>0&&nAS>0) // Cancel strange sea
645  {
646  G4int ss= nS;
647  if(nAS<nS) ss=nAS;
648  nS -=ss;
649  nAS-=ss;
650  }
651  if (!b)DecQAQ(-4);
652  else if(b==1)DecQAQ(-5);
653  else DecQAQ(0);
654  G4int tot=GetTot();
655  G4int q=GetQ();
656  G4int aq=GetAQ();
657  G4QContent r=Pi; // Pion prototype of returned value
658  if((tot!=4||q!=2) && (tot!=5||(q!=1&&aq!=1)) && (tot!=6||abs(b)!=2))
659  {
660  //#ifdef erdebug
661  G4cerr<<"***G4QCont::SplitChipo: QC="<<GetThis()<<G4endl;
662  //#endif
663  }
664  else if(tot==4) // Mesonic (eight possibilities)
665  {
666  r=GetThis();
667  if (r.SubtractPi0()) ; // Try any trivial algorithm of splitting
668  else if(r.SubtractPion()) ;
669  else if(r.SubtractKaon(mQ)) ;
670  else
671  {
672  //#ifdef debug
673  G4cerr<<"***G4QCont::SplitChipo:MesTot="<<tot<<",b="<<b<<",q="<<q<<",a="<<aq<<G4endl;
674  //#endif
675  }
676  }
677  else if(b==1&&tot==5) // Baryonic (four possibilities)
678  {
679  if(nU==3)
680  {
681  r.SetU(1);
682  r+=IndAQ();
683  }
684  else if(nD==3)
685  {
686  r.SetD(1);
687  r+=IndAQ();
688  }
689  else if(nS==3)
690  {
691  r.SetS(1);
692  r+=IndAQ();
693  }
694  else if(nAU==3)
695  {
696  r.SetAU(1);
697  r+=IndQ();
698  }
699  else if(nAD==3)
700  {
701  r.SetAD(1);
702  r+=IndQ();
703  }
704  else if(nAS==3)
705  {
706  r.SetAS(1);
707  r+=IndQ();
708  }
709  else if(q==1&&nU)
710  {
711  r.SetU(1);
712  if(nAU) r.SetAU(1);
713  else r.SetAD(1);
714  }
715  else if(q==1&&nD)
716  {
717  r.SetD(1);
718  if(nAD) r.SetAD(1);
719  else r.SetAU(1);
720  }
721  else if(q==1&&nS)
722  {
723  r.SetS(1);
724  if(nAS) r.SetAS(1);
725  else r.SetAU(1);
726  }
727  else if(aq==1&&nAU)
728  {
729  r.SetAU(1);
730  if(nU) r.SetU(1);
731  else r.SetD(1);
732  }
733  else if(aq==1&&nAD)
734  {
735  r.SetAD(1);
736  if(nD) r.SetD(1);
737  else r.SetU(1);
738  }
739  else if(aq==1&&nAS)
740  {
741  r.SetAS(1);
742  if(nS) r.SetS(1);
743  else r.SetU(1);
744  }
745  else
746  {
747  //#ifdef erdebug
748  G4cerr<<"***G4QCont::SplitChipo: Baryonic tot=5,b=1,qCont="<<GetThis()<<G4endl;
749  //#endif
750  }
751  }
752  else if(tot==b*3) // MultyBaryon cace
753  {
754  r=GetThis();
755  if (bn>0) // baryonium
756  {
757  G4QContent la(1,1,1,0,0,0);
758  G4QContent nt(2,1,0,0,0,0);
759  G4QContent pr(1,2,0,0,0,0);
760  G4QContent ks(0,1,2,0,0,0);
761  if (nD>nU) ks=G4QContent(1,0,2,0,0,0);
762  G4QContent dm(3,0,0,0,0,0);
763  G4QContent dp(0,3,0,0,0,0);
764  G4QContent om(0,0,3,0,0,0);
765  if (nU>=nD&&nU>=nS)
766  {
767  if (r.SubtractHadron(pr)) r-=pr; // These functions only check
768  else if(r.SubtractHadron(dp)) r-=dp;
769  else if(r.SubtractHadron(nt)) r-=nt;
770  else if(r.SubtractHadron(la)) r-=la;
771  else if(r.SubtractHadron(dm)) r-=dm;
772  else
773  {
774  //#ifdef erdebug
775  G4cerr<<"***G4QCont::SplitChipo:Dibar (1) tot=6, b=2, qCont="<<GetThis()<<G4endl;
776  //#endif
777  }
778  }
779  else if(nD>=nU&&nD>=nS)
780  {
781  if (r.SubtractHadron(nt)) r-=nt; // These functions only check
782  else if(r.SubtractHadron(dm)) r-=dm;
783  else if(r.SubtractHadron(pr)) r-=pr;
784  else if(r.SubtractHadron(dp)) r-=dp;
785  else if(r.SubtractHadron(la)) r-=la;
786  else
787  {
788  //#ifdef erdebug
789  G4cerr<<"***G4QContent::SplitChipo:Dib(2) tot=6, b=2, qCont="<<GetThis()<<G4endl;
790  //#endif
791  }
792  }
793  else
794  {
795  if (r.SubtractHadron(la)) r-=la; // These functions only check
796  else if(r.SubtractHadron(ks)) r-=ks;
797  else if(r.SubtractHadron(om)) r-=om;
798  else if(r.SubtractHadron(pr)) r-=pr;
799  else if(r.SubtractHadron(nt)) r-=nt;
800  else
801  {
802  //#ifdef erdebug
803  G4cerr<<"***G4QContent::SplitChipo:Dib(3) tot=6, b=2, qCont="<<GetThis()<<G4endl;
804  //#endif
805  }
806  }
807  }
808  else // Anti-baryonium
809  {
810  G4QContent la(0,0,0,1,1,1);
811  G4QContent pr(0,0,0,1,2,0);
812  G4QContent nt(0,0,0,2,1,0);
813  G4QContent ks(0,1,2,0,0,0);
814  if(nAD>nAU)ks=G4QContent(0,0,0,1,0,2);
815  G4QContent dm(0,0,0,3,0,0);
816  G4QContent dp(0,0,0,0,3,0);
817  G4QContent om(0,0,0,0,0,3);
818  if (nAU>=nAD&&nAU>=nAS)
819  {
820  if (r.SubtractHadron(pr)) r-=pr; // These functions only check
821  else if(r.SubtractHadron(dp)) r-=dp;
822  else if(r.SubtractHadron(nt)) r-=nt;
823  else if(r.SubtractHadron(la)) r-=la;
824  else if(r.SubtractHadron(dm)) r-=dm;
825  else
826  {
827  //#ifdef erdebug
828  G4cerr<<"***G4QContent::SplitChipo:ADib(1) tot=6,b=2, qCont="<<GetThis()<<G4endl;
829  //#endif
830  }
831  }
832  else if(nAD>=nAU&&nAD>=nAS)
833  {
834  if (r.SubtractHadron(nt)) r-=nt; // These functions only check
835  else if(r.SubtractHadron(dm)) r-=dm;
836  else if(r.SubtractHadron(pr)) r-=pr;
837  else if(r.SubtractHadron(dp)) r-=dp;
838  else if(r.SubtractHadron(la)) r-=la;
839  else
840  {
841  //#ifdef erdebug
842  G4cerr<<"***G4QContent::SplitChipo:ADib(2) tot=6,b=2, qCont="<<GetThis()<<G4endl;
843  //#endif
844  }
845  }
846  else
847  {
848  if (r.SubtractHadron(la)) r-=la; // These functions only check
849  else if(r.SubtractHadron(ks)) r-=ks;
850  else if(r.SubtractHadron(om)) r-=om;
851  else if(r.SubtractHadron(pr)) r-=pr;
852  else if(r.SubtractHadron(nt)) r-=nt;
853  else
854  {
855  //#ifdef erdebug
856  G4cerr<<"***G4QContent::SplitChipo:ADib(3) tot=6,b=2, qCont="<<GetThis()<<G4endl;
857  //#endif
858  }
859  }
860  }
861  }
862  else // More than Dibaryon (@@ can use the same algorithm as for dibaryon)
863  {
864  //#ifdef erdebug
865  G4cerr<<"*G4QContent::SplitChipolino:UnknownHadron with QuarkCont="<<GetThis()<<G4endl;
866  //#endif
867  }
868  return r;
869 }// End of G4QContent::SplitChipolino
870 
871 // Return one-quark QC using index (a kind of iterator)
872 G4QContent G4QContent::IndQ (G4int index)
873 {
874 #ifdef debug
875  G4cout << "G4QC::IndQ is called"<<G4endl;
876 #endif
877  if(index<nD) return G4QContent(1,0,0,0,0,0);
878  else if(index<nD+nU) return G4QContent(0,1,0,0,0,0);
879  else if(index<nD+nU+nS) return G4QContent(0,0,1,0,0,0);
880  //#ifdef erdebug
881  else G4cerr<<"***G4QC::IndQ:index="<<index<<" for the QuarkContent="<<GetThis()<<G4endl;
882  //throw G4QException("***G4QC::IndQ: Index exceeds the total number of quarks");
883  //#endif
884  return G4QContent(0,0,0,0,0,0);
885 }
886 
887 // Return one-antiquark QC using index (a kind of iterator)
888 G4QContent G4QContent::IndAQ (G4int index)
889 {
890 #ifdef debug
891  G4cout << "G4QC::IndAQ is called"<<G4endl;
892 #endif
893  if(index<nAD) return G4QContent(0,0,0,1,0,0);
894  else if(index<nAD+nAU) return G4QContent(0,0,0,0,1,0);
895  else if(index<nAD+nAU+nAS) return G4QContent(0,0,0,0,0,1);
896  //#ifdef erdebug
897  else G4cerr<<"***G4QC::IndAQ:index="<<index<<" for the QuarkContent="<<GetThis()<<G4endl;
898  //throw G4QException("***G4QC::IndAQ: Index exceeds the total number of antiquarks");
899  //#endif
900  return G4QContent(0,0,0,0,0,0);
901 }
902 
903 // Reduces number (if nQAQ<0:all) of valence Q-Qbar pairs, returns #of pairs to reduce more
904 G4int G4QContent::DecQAQ(const G4int& nQAQ)
905 {
906 #ifdef debug
907  G4cout<<"G4QCont::DecQC: n="<<nQAQ<<","<<GetThis()<<G4endl;
908 #endif
909  G4int ban = GetBaryonNumber();
910  G4int tot = GetTot(); // Total number of quarks in QC
911  if (tot==ban*3) return 0;// Nothing to reduce & nothing to reduce more
912  G4int nUP=0; // U/AU min factor (a#of u which can be subtracted)
913  if (nU>=nAU) nUP=nAU;
914  else nUP= nU;
915 
916  G4int nDP=0; // D/AD min factor (a#of d which can be subtracted)
917  if (nD>=nAD) nDP=nAD;
918  else nDP= nD;
919 
920  G4int nSP=0; // S/AS min factor (a#of s which can be subtracted)
921  if (nS>=nAS) nSP=nAS;
922  else nSP= nS;
923 
924  G4int nLP =nUP+nDP; // a#of light quark pairs which can be subtracted
925  G4int nTotP=nLP+nSP; // total#of existing pairs which can be subtracted
926  G4int nReal=nQAQ; // demanded #of pairs for reduction (by demand)
927  G4int nRet =nTotP-nQAQ; // a#of additional pairs for further reduction
928  if (nQAQ<0) // === Limited reduction case @@ not tuned for baryons !!
929  {
930  G4int res=tot+nQAQ;
931 #ifdef debug
932  G4cout<<"G4QC::DecQC: tot="<<tot<<", nTP="<<nTotP<<", res="<<res<<G4endl;
933 #endif
934  if(res<0)
935  {
936  IncQAQ(1,0.); // Increment by one not strange pair to get the minimum
937  return 1;
938  }
939  res -=nTotP+nTotP;
940  if(res<0) nReal=nTotP+res/2;
941  else nReal=nTotP;
942  nRet = tot/2-nReal;
943  }
944  else if(!nQAQ)
945  {
946  nReal=nTotP;
947  nRet =0;
948  }
949  else if(nRet<0) nReal=nTotP;
950 
951  if (!nReal) return nRet; // Now nothing to be done
952  // ---------- Decrimenting by nReal pairs
953 #ifdef debug
954  G4cout<<"G4QC::DecQC: demanded "<<nQAQ<<" pairs, executed "<<nReal<<" pairs"<<G4endl;
955 #endif
956  for (G4int i=0; i<nReal; i++)
958  {
959  G4double base = nTotP;
960  //if (nRet && nSP==1 && !nQAQ) base = nLP; // Keep S-Sbar pair if possible
961  G4int j = static_cast<int>(base*G4UniformRand()); // Random integer "SortOfQuark"
962  if (nUP && j<nUP && (nRet>2 || nUP>1 || (nD<2 && nS<2)))// --- U-Ubar pair
963  {
964 #ifdef debug
965  G4cout<<"G4QC::DecQC: decrementing UAU pair UP="<<nUP<<", QC="<<GetThis()<<G4endl;
966 #endif
967  nU--;
968  nAU--;
969  nUP--;
970  nLP--;
971  nTotP--;
972  }
973  else if (nDP && j<nLP && (nRet>2 || nDP>1 || (nU<2 && nS<2)))// --- D-Ubar pair
974  {
975 #ifdef debug
976  G4cout<<"G4QC::DecQC: decrementing DAD pair DP="<<nDP<<", QC="<<GetThis()<<G4endl;
977 #endif
978  nD--;
979  nAD--;
980  nDP--;
981  nLP--;
982  nTotP--;
983  }
984  else if (nSP&& (nRet>2 || nSP>1 || (nU<2 && nD<2))) // --- S-Sbar pair
985  {
986 #ifdef debug
987  G4cout<<"G4QC::DecQC: decrementing SAS pair SP="<<nSP<<", QC="<<GetThis()<<G4endl;
988 #endif
989  nS--;
990  nAS--;
991  nSP--;
992  nTotP--;
993  }
994  else if (nUP) // --- U-Ubar pair cancelation (final)
995  {
996 #ifdef debug
997  G4cout<<"G4QC::DecQC:Decrement UAU pair (final) UP="<<nUP<<",QC="<<GetThis()<<G4endl;
998 #endif
999  nU--;
1000  nAU--;
1001  nUP--;
1002  nLP--;
1003  nTotP--;
1004  }
1005  else if (nDP) // --- D-Ubar pair cancelation (final)
1006  {
1007 #ifdef debug
1008  G4cout<<"G4QC::DecQC:Decrement DAD pair (final) DP="<<nDP<<",QC="<<GetThis()<<G4endl;
1009 #endif
1010  nD--;
1011  nAD--;
1012  nDP--;
1013  nLP--;
1014  nTotP--;
1015  }
1016  else if (nSP) // --- S-Sbar pair cancelation (final)
1017  {
1018 #ifdef debug
1019  G4cout<<"G4QC::DecQC: decrementing SAS pair SP="<<nSP<<", QC="<<GetThis()<<G4endl;
1020 #endif
1021  nS--;
1022  nAS--;
1023  nSP--;
1024  nTotP--;
1025  }
1026  else G4cout<<"***G4QC::DecQC:i="<<i<<",j="<<j<<",D="<<nDP<<",U="<<nUP<<",S="<<nSP
1027  <<",T="<<nTotP<<",nRet="<<nRet<<", QC="<<GetThis()<<G4endl;
1028  }
1029 #ifdef debug
1030  G4cout<<"G4QC::DecQC: >->-> OUT <-<-< nRet="<<nRet<<", QC="<<GetThis()<<G4endl;
1031 #endif
1032  return nRet;
1033 }
1034 
1035 // Increment quark pairs
1036 void G4QContent::IncQAQ(const G4int& nQAQ, const G4double& sProb)
1037 {
1038  G4int tot = GetTot();
1039  G4QContent mQC = GetThis();
1040  for (int i=0; i<nQAQ; i++)
1041  {
1042  G4int j = static_cast<int>((2.+sProb)*G4UniformRand()); // 0-U, 1-D, 2-S
1043 #ifdef debug
1044  G4cout<<"IncQC:out QC="<<GetThis()<<",j="<<j<<" for i="<<i<<G4endl;
1045 #endif
1046  //if (!j)
1047  if ( !j && (nU<=nD || nU<=nS))
1048  {
1049  nU++;
1050  nAU++;
1051  tot+=2;
1052  }
1053  //else if (j==1)
1054  else if (j==1 && (nD<=nU || nD<=nS))
1055  {
1056  nD++;
1057  nAD++;
1058  tot+=2;
1059  }
1060  //else
1061  else if (j>1&& (nS<=nU || nS<=nD))
1062  {
1063  nS++;
1064  nAS++;
1065  tot+=2;
1066  }
1067  else if (!j)
1068  {
1069  nD++;
1070  nAD++;
1071  tot+=2;
1072  }
1073  else if (j==1)
1074  {
1075  nU++;
1076  nAU++;
1077  tot+=2;
1078  }
1079  else
1080  {
1081  nS++;
1082  nAS++;
1083  tot+=2;
1084  }
1085  //else if (nD<=nU)
1086  //{
1087  // nD++;
1088  // nAD++;
1089  // tot+=2;
1090  //}
1091  //else
1092  //{
1093  // nU++;
1094  // nAU++;
1095  // tot+=2;
1096  //}
1097  }
1098 }
1099 
1100 // Calculate a#of protons in the QC (for nuclei)
1101 G4int G4QContent::GetP() const
1102 {
1103  G4int rD=nD-nAD; // Constituent d-quarks
1104  G4int rU=nU-nAU; // Constituent u-quarks
1105  G4int rS=nS-nAS; // Constituent s-quarks
1106  G4int dQ=rD-rU; // Isotopic assimetry
1107  G4int b3=rD+rU+rS; // (Baryon number) * 3
1108  return (b3-3*(rS+dQ))/6;
1109 }
1110 
1111 // Calculate a#of neutrons in the QC (for nuclei)
1112 G4int G4QContent::GetN() const
1113 {
1114  G4int rD=nD-nAD; // Constituent d-quarks
1115  G4int rU=nU-nAU; // Constituent u-quarks
1116  G4int rS=nS-nAS; // Constituent s-quarks
1117  G4int dQ=rD-rU; // Isotopic assimetry
1118  G4int b3=rD+rU+rS; // (Baryon number) * 3
1119  return (b3-3*(rS-dQ))/6;
1120 }
1121 
1122 // Calculate a#of lambdas in the QC (for nuclei)
1123 G4int G4QContent::GetL() const
1124 {
1125  G4int rS=nS-nAS; // Constituent s-quarks
1126  return rS;
1127 }
1128 
1129 // Calculate a#of anti-protons in the QC (for anti-nuclei)
1130 G4int G4QContent::GetAP() const
1131 {
1132  G4int rD=nAD-nD; // Constituent anti-d-quarks
1133  G4int rU=nAU-nU; // Constituent anti-u-quarks
1134  G4int rS=nAS-nS; // Constituent anti-s-quarks
1135  G4int dQ=rD-rU; // Isotopic assimetry
1136  G4int b3=rD+rU+rS; // - (Baryon number) * 3
1137  return (b3-3*(rS+dQ))/6;
1138 }
1139 
1140 // Calculate a#of anti-neutrons in the QC (for anti-nuclei)
1141 G4int G4QContent::GetAN() const
1142 {
1143  G4int rD=nAD-nD; // Constituent anti-d-quarks
1144  G4int rU=nAU-nU; // Constituent anti-u-quarks
1145  G4int rS=nAS-nS; // Constituent anti-s-quarks
1146  G4int dQ=rD-rU; // Isotopic assimetry
1147  G4int b3=rD+rU+rS; // - (Baryon number) * 3
1148  return (b3-3*(rS-dQ))/6;
1149 }
1150 
1151 // Calculate a#of anti-lambdas in the QC (for anti-nuclei)
1152 G4int G4QContent::GetAL() const
1153 {
1154  G4int rS=nAS-nS; // Constituent anti-s-quarks
1155  return rS;
1156 }
1157 
1158 // Calculate charge for the QC
1159 G4int G4QContent::GetCharge() const
1160 {
1161  static const G4int cU = 2;
1162  static const G4int cD =-1;
1163  static const G4int cS =-1;
1164  static const G4int cAU =-2;
1165  static const G4int cAD = 1;
1166  static const G4int cAS = 1;
1167 
1168  G4int c=0;
1169  if(nU) c+=nU*cU;
1170  if(nD) c+=nD*cD;
1171  if(nS) c+=nS*cS;
1172  if(nAU)c+=nAU*cAU;
1173  if(nAD)c+=nAD*cAD;
1174  if(nAS)c+=nAS*cAS;
1175  //#ifdef erdebug
1176  if(c%3) G4cerr<<"***G4QCont:GetCharge:c="<<c<<"/3 isn't integer, QC="<<GetThis()<<G4endl;
1177  //#endif
1178  return c/3;
1179 }
1180 
1181 // Calculate a Baryon Number for the QC
1182 G4int G4QContent::GetBaryonNumber() const
1183 {
1184 #ifdef pdebug
1185  G4cout<<"G4QContent::GetBarNum: U="<<nU<<", D="<<nD<<", S="<<nS<<", AU="<<nAU<<", AD="
1186  <<nAD<<", AS="<<nAS<<G4endl;
1187 #endif
1188  G4int b=nU+nD+nS-nAU-nAD-nAS;
1189  //#ifdef erdebug
1190  if(b%3)
1191  {
1192  // G4cerr<<"-Warning-G4QContent::GetBaryonNumber="<<b<<"/3 isn't anIntegerValue"<<G4endl;
1193  // G4Exception("G4QContent::GetBaryonNumber:","72",FatalException,"Wrong Baryon Number");
1194  G4ExceptionDescription ed;
1195  ed << "Wrong Baryon Number: warning " << b << "/3 isn't an integer"
1196  << G4endl;
1197  G4Exception("G4QContent::GetBaryonNumber()", "HAD_CHPS_0072", FatalException, ed);
1198  }
1199  //#endif
1200  return b/3;
1201 }
1202 
1203 // Gives the PDG of the lowest (in mass) S-wave hadron or Chipolino (=10) for double hadron
1204 G4int G4QContent::GetSPDGCode() const
1205 {
1206  G4int p = 0; // Prototype of output SPDG
1207  G4int n = GetTot(); // Total number of quarks
1208  if(!n) return 22; // Photon does not have any Quark Content
1209  G4int mD=nD; // A # of D quarks or anti U quarks
1210  if (nD<=0) mD=nAD;
1211  G4int mU=nU; // A # of U quarks or anti U quarks
1212  if (nU<=0) mU=nAU;
1213  G4int mS=nS; // A # of S quarks or anti U quarks
1214  if (nS<=0) mS= nAS;
1215  // ---------------------- Cancelation of q-qbar pairs in case of an excess
1216  if ( nU>nAU && nAU>0)
1217  {
1218  mU=nU-nAU;
1219  n-=nAU+nAU;
1220  }
1221  if (nAU>nU && nU>0)
1222  {
1223  mU=nAU-nU;
1224  n-=nU+nU;
1225  }
1226  if ( nD>nAD && nAD>0)
1227  {
1228  mD=nD-nAD;
1229  n-=nAD+nAD;
1230  }
1231  if (nAD>nD && nD>0)
1232  {
1233  mD=nAD-nD;
1234  n-=nD+nD;
1235  }
1236  if ( nS>nAS && nAS>0)
1237  {
1238  mS=nS-nAS;
1239  n-=nAS+nAS;
1240  }
1241  if (nAS>nS && nS>0)
1242  {
1243  mS= nAS-nS;
1244  n-=nS+nS;
1245  }
1246  // ---------------------- Cancelation of q-qbar pairs in case of an equality
1247  if (nAD==nD && nD>0)
1248  {
1249  G4int dD=nD+nD;
1250  if(n>dD)
1251  {
1252  mD=0;
1253  n-=dD;
1254  }
1255  else if (n==dD)
1256  {
1257  mD=2;
1258  n=2;
1259  }
1260  else
1261  {
1262 #ifdef debug
1263  G4cout<<"***G4QC::SPDG:CanD U="<<mU<<",D="<<mD<<",S="<<mS<<",QC="<<GetThis()<<G4endl;
1264 #endif
1265  return 0;
1266  }
1267  }
1268  if (nAU==nU && nU>0)
1269  {
1270  G4int dU=nU+nU;
1271  if(n>dU)
1272  {
1273  mU=0;
1274  n-=dU;
1275  }
1276  else if (n==dU)
1277  {
1278  mU=2;
1279  n=2;
1280  }
1281  else
1282  {
1283 #ifdef debug
1284  G4cout<<"***G4QC::SPDG:CanU U="<<mU<<",D="<<mD<<",S="<<mS<<",QC="<<GetThis()<<G4endl;
1285 #endif
1286  return 0;
1287  }
1288  }
1289  if (nAS==nS && nS>0) //@@ Starts with S-quarks - should be randomized and mass limited
1290  {
1291  G4int dS=nS+nS;
1292  if(n>dS)
1293  {
1294  mS=0;
1295  n-=dS;
1296  }
1297  else if (n==dS)
1298  {
1299  mS=2;
1300  n=2;
1301  }
1302  else
1303  {
1304 #ifdef debug
1305  G4cout<<"***G4QC::SPDG:CanS U="<<mU<<",D="<<mD<<",S="<<mS<<",QC="<<GetThis()<<G4endl;
1306 #endif
1307  return 0;
1308  }
1309  }
1310 
1311  G4int b=GetBaryonNumber();
1312  G4int c=GetCharge();
1313  G4int s_value=GetStrangeness();
1314 #ifdef debug
1315  G4cout<<"G4QC::SPDGC:bef. b="<<b<<",n="<<n<<",c="<<c<<",s="<<s_value<<",Q="<<GetThis()<<G4endl;
1316 #endif
1317  if (b) // =---------------------= Baryon case
1318  {
1319 
1320  G4int ab=abs(b);
1321  if(ab>=2 && n>=6) // Multi-Baryonium (NuclearFragment)
1322  {
1323  G4int mI=nU-nAU-nD+nAD;
1324  //if (abs(mI)>3||mS>3||(b>0&&s_value<-1)||(b<0&&s_value>1)) return 0;
1325  //else if(abs(mI)>2||mS>2||(b>0&&s_value< 0)||(b<0&&s_value>0)) return 10;
1326  if ( (b > 0 && s_value < -1) || (b < 0 && s_value > 1) ) return 10;
1327  else if (abs(mI) > 2 || mS > 2
1328  || (b > 0 && s_value < 0)
1329  || (b < 0 && s_value > 0)) return GetZNSPDGCode();
1330  else if(mU>=mS&&mD>=mS&&mU+mD+mS==3*b) // Possible Unary Nuclear Cluster
1331  {
1332  G4int mZ=(mU+mD-mS-mS+3*mI)/6;
1333  p = 90000000+1000*(1000*mS+mZ)+mZ-mI;
1334  if(b>0) return p;
1335  else return -p;
1336  }
1337  else return 10;
1338  }
1339  // Normal One Baryon States: Heavy quark should come first
1340  if(n>5) return GetZNSPDGCode(); //B+M+M Tripolino etc
1341  if(n==5) return 10; //B+M Chipolino
1342  if(mS>0) // Strange Baryons
1343  {
1344  p=3002;
1345  if (mS==3) p+=332; // Decuplet
1346  else if (mS==2)
1347  {
1348  if (mU==1 && mD==0) p+=320;
1349  else if (mU==0 && mD==1) p+=310;
1350  else
1351  {
1352 #ifdef debug
1353  G4cout<<"**G4QC::SPDG:ExoticBSS,U="<<mU<<",D="<<mD<<",S="<<mS<<GetThis()<<G4endl;
1354 #endif
1355  return GetZNSPDGCode();
1356  }
1357  }
1358  else if (mS==1)
1359  {
1360  if (mU==2 && mD==0) p+=220;
1361  else if (mU==1 && mD==1) p+=120; // Lambda (M_Lambda<M_Sigma0) PDG_Sigma=3212
1362  else if (mU==0 && mD==2) p+=110;
1363  else
1364  {
1365 #ifdef debug
1366  G4cout<<"***G4QC::SPDG:ExoticBS,U="<<mU<<",D="<<mD<<",S="<<mS<<GetThis()<<G4endl;
1367 #endif
1368  return GetZNSPDGCode();
1369  }
1370  }
1371  else // Superstrange case
1372  {
1373 #ifdef debug
1374  G4cout<<"***G4QC::GetSPDG:ExoBarS,U="<<mU<<",D="<<mD<<",S="<<mS<<GetThis()<<G4endl;
1375 #endif
1376  return GetZNSPDGCode();
1377  }
1378  }
1379  else if (mU>0) // Not Strange Baryons
1380  {
1381  p=2002;
1382  if (mU==3 && mD==0) p+=222; // Decuplet
1383  else if (mU==2 && mD==1) p+=210;
1384  else if (mU==1 && mD==2) p+=110; // There is a higher Delta S31 (PDG=1212)
1385  else
1386  {
1387 #ifdef debug
1388  G4cout<<"***G4QC::SPDG:ExoBaryonU,U="<<mU<<",D="<<mD<<",S="<<mS<<GetThis()<<G4endl;
1389 #endif
1390  return GetZNSPDGCode();
1391  }
1392  }
1393  else if (mD==3) p=1114; // Decuplet
1394  else
1395  {
1396 #ifdef debug
1397  G4cout<<"**G4QC::SPDG:ExoticBaryonD,U="<<mU<<",D="<<mD<<",S="<<mS<<GetThis()<<G4endl;
1398 #endif
1399  return GetZNSPDGCode();
1400  }
1401  if (b<0) p=-p;
1402  }
1403  else // ------------------------->>------------------------>> Meson case
1404  {
1405 #ifdef debug
1406  G4cout<<"G4QC::SPDG:mDUS="<<mD<<","<<mU<<","<<mS<<",b,c,s="<<b<<","<<c<<","<<s_value<<G4endl;
1407 #endif
1408  if(n>4) // Super Exotics
1409  {
1410 #ifdef debug
1411  G4cout<<"G4QC::SPDG:n>4 SEx:U="<<mU<<",D="<<mD<<",S="<<mS<<",QC="<<GetThis()<<G4endl;
1412 #endif
1413  return 0;
1414  }
1415  if(n==4) return 10; // M+M Chipolino
1416  if(abs(s_value)>1)
1417  {
1418 #ifdef debug
1419  G4cout<<"**G4QC::SPDG:Stran="<<s_value<<",QC="<<GetThis()<<" - Superstrange Meson"<<G4endl;
1420 #endif
1421  return 0;
1422  }
1423  // Heavy quark should come first
1424  if(mS>0) // Strange Mesons
1425  {
1426  p=301;
1427  if (mS==2)
1428  {
1429  //if (G4UniformRand()<0.333) p=221; // eta
1430  //else p+=30; // eta'
1431  p=221;
1432  }
1433  else if (mU==1 && mD==0) p+=20;
1434  else if (mU==0 && mD==1) p+=10;
1435  else
1436  {
1437 #ifdef debug
1438  G4cout<<"*G4QC::SPDG:ExMS U="<<mU<<",D="<<mD<<",S="<<mS<<",QC="<<GetThis()<<G4endl;
1439 #endif
1440  return 0;
1441  }
1442  }
1443  else if (mU>0) // Isotopic Mesons
1444  {
1445  p=201;
1446  //if (mU==2 && mD==0) p=221; // Performance problems
1447  if (mU==2 && mD==0) p=111;
1448  else if (mU==1 && mD==1) p+=10;
1449  else
1450  {
1451 #ifdef debug
1452  G4cout<<"*G4QC::SPDG:ExMU U="<<mU<<",D="<<mD<<",S="<<mS<<",QC="<<GetThis()<<G4endl;
1453 #endif
1454  return 0;
1455  }
1456  }
1457  else if (mD==2) p=111;
1458  else
1459  {
1460 #ifdef debug
1461  G4cout<<"***G4QC::SPDG:ExMD U="<<mU<<",D="<<mD<<",S="<<mS<<",QC="<<GetThis()<<G4endl;
1462 #endif
1463  return 0;
1464  }
1465  if (c<0 || (c==0 && mS>0 && s_value>0)) p=-p;
1466  }
1467 #ifdef debug
1468  G4cout<<"G4QC::GetSPDG:output SPDGcode="<<p<<" for the QuarkContent="<<GetThis()<<G4endl;
1469 #endif
1470  return p;
1471 }
1472 
1473 // === Calculate a number of combinations of rhc out of lhc ==
1474 G4int G4QContent::NOfCombinations(const G4QContent& rhs) const
1475 {
1476  G4int c=1; // Default number of combinations?
1477 #ifdef ppdebug
1478  G4cout<<"G4QContent::NOfComb: This="<<GetThis()<<", selectQC="<<rhs<<G4endl;
1479 #endif
1480  G4int mD=rhs.GetD();
1481  G4int mU=rhs.GetU();
1482  G4int mS=rhs.GetS();
1483  G4int mAD=rhs.GetAD();
1484  G4int mAU=rhs.GetAU();
1485  G4int mAS=rhs.GetAS();
1486  G4int mN=mD+mU+mS-mAD-mAU-mAS;
1488  if (( ((nD < mD || nAD < mAD) && !(mD-mAD)) ||
1489  ((nU < mU || nAU < mAU) && !(mU-mAU)) ||
1490  ((nS < mS || nAS < mAS) && !(mS-mAS)) ) && !mN) return 1;
1491  if(mD>0)
1492  {
1493  int j=nD;
1494  if (j<=0) return 0;
1495  if(mD>1||j>1) for (int i=1; i<=mD; i++)
1496  {
1497  if(!j) return 0;
1498  c*=j/i;
1499  j--;
1500  }
1501  }
1502  if(mU>0)
1503  {
1504  int j=nU;
1505  if (j<=0) return 0;
1506  if(mU>1||j>1) for (int i=1; i<=mU; i++)
1507  {
1508  if(!j) return 0;
1509  c*=j/i;
1510  j--;
1511  }
1512  }
1513  if(mS>0)
1514  {
1515  int j=nS;
1516  if (j<=0) return 0;
1517  if(mS>1||j>1) for (int i=1; i<=mS; i++)
1518  {
1519  if(!j) return 0;
1520  c*=j/i;
1521  j--;
1522  }
1523  }
1524  if(mAD>0)
1525  {
1526  int j=nAD;
1527  if (j<=0) return 0;
1528  if(mAD>1||j>1) for (int i=1; i<=mAD; i++)
1529  {
1530  if(!j) return 0;
1531  c*=j/i;
1532  j--;
1533  }
1534  }
1535  if(mAU>0)
1536  {
1537  int j=nAU;
1538  if (j<=0) return 0;
1539  if(mAU>1||j>1) for (int i=1; i<=mAU; i++)
1540  {
1541  if(!j) return 0;
1542  c*=j/i;
1543  j--;
1544  }
1545  }
1546  if(mAS>0)
1547  {
1548  int j=nAS;
1549  if (j<=0) return 0;
1550  if(mAS>1||j>1) for (int i=1; i<=mAS; i++)
1551  {
1552  if(!j) return 0;
1553  c*=j/i;
1554  j--;
1555  }
1556  }
1557  return c;
1558 }
1559 
1560 // Make PDG's of PartonPairs for Mesons & Baryons (only)
1561 std::pair<G4int,G4int> G4QContent::MakePartonPair() const
1562 {
1563  G4double S=0.;
1564  S+=nD;
1565  G4double dP=S;
1566  S+=nU;
1567  G4double uP=S;
1568  S+=nS;
1569  G4double sP=S;
1570  S+=nAD;
1571  G4double dA=S;
1572  S+=nAU;
1573  G4double uA=S;
1574  S+=nAS;
1575  if(!S)
1576  {
1577  G4int f= static_cast<int>(1.+2.3*G4UniformRand()); // Random flavor @@ a Parameter
1578  return std::make_pair(f,-f);
1579  }
1580  G4int f=0;
1581  G4double R=S*G4UniformRand();
1582  if (R<dP) f=1;
1583  else if(R<uP) f=2;
1584  else if(R<sP) f=3;
1585  else if(R<dA) f=-1;
1586  else if(R<uA) f=-2;
1587  else f=-3;
1588 
1589  if (f < 0) { // anti-quark
1590  if(nD || nU || nS) // a Meson
1591  {
1592  if (nD) return std::make_pair(1,f);
1593  else if(nU) return std::make_pair(2,f);
1594  else return std::make_pair(3,f);
1595  }
1596  else // Anti-Baryon
1597  {
1598  // @@ Can be improved, taking into acount weights (i,i): w=3, (i,j#i): w=4(s=0 + s=1)
1599  G4int AD=nAD;
1600  if(f==-1) AD--;
1601  G4int AU=nAU;
1602  if(f==-1) AU--;
1603  G4int AS=nAS;
1604  if(f==-1) AS--;
1605  if (AS)
1606  {
1607  if (AS==2) return std::make_pair(-3303,f); // 3301 does not exist
1608  else if(AU) return std::make_pair(-3201,f); // @@ only lightest
1609  else return std::make_pair(-3101,f); // @@ only lightest
1610  }
1611  else if(AU)
1612  {
1613  if (AU==2) return std::make_pair(-2203,f); // 2201 does not exist
1614  else return std::make_pair(-2101,f); // @@ only lightest
1615  }
1616  else return std::make_pair(-1103,f); // 1101 does not exist
1617  }
1618 
1619  } else { // quark (f is a PDG code of the quark)
1620 
1621  if(nAD || nAU || nAS) // a Meson
1622  {
1623  if (nAD) return std::make_pair(f,-1);
1624  else if(nAU) return std::make_pair(f,-2);
1625  else return std::make_pair(f,-3);
1626  }
1627  else // Anti-Baryon
1628  {
1629  // @@ Can be improved, taking into acount weights (i,i): w=3, (i,j#i): w=4(s=0 + s=1)
1630  // G4int AD=nD; AD is unused
1631  // if(f==-1) AD--; dead code: f >= 0 in this branch (DHW 11 Nov 2010)
1632  G4int AU=nU;
1633  // if(f==-1) AU--; dead code: f >= 0 in this branch (DHW 11 Nov 2010)
1634  G4int AS=nS;
1635  // if(f==-1) AS--; dead code: f >= 0 in this branch (DHW 11 Nov 2010)
1636 
1637  if (AS) {
1638  if (AS==2) return std::make_pair(f,3303); // 3301 does not exist
1639  else if(AU) return std::make_pair(f,3201); // @@ only lightest
1640  else return std::make_pair(f,3101); // @@ only lightest
1641 
1642  } else if(AU) {
1643  if (AU==2) return std::make_pair(f,2203); // 2201 does not exist
1644  else return std::make_pair(f,2101); // @@ only lightest
1645  }
1646  else return std::make_pair(f,1103); // 1101 does not exist
1647  }
1648  } // test on f
1649 }
1650 
1651 
1652 // Add parton (pPDG) to the hadron (this QC) & get parton PDG
1653 // (Baryons,Mesons,Anti-Baryons)
1654 
1655 G4int G4QContent::AddParton(G4int pPDG) const
1656 {
1657 #ifdef debug
1658  G4cout<<"G4QContent::AddParton: This="<<GetThis()<<", pPDG="<<pPDG<<G4endl;
1659 #endif
1660  if(!pPDG || pPDG==9 || pPDG==21)
1661  {
1662 #ifdef debug
1663  G4cout<<"-Warning-G4QContent::AddParton: ImpossibleToAdd PartonWithPDG="<<pPDG<<G4endl;
1664 #endif
1665  return 0;
1666  }
1667  G4int aPDG = std::abs(pPDG);
1668  if( (aPDG>3 && aPDG<1101) || pPDG>3303) // @@ 1101 does not exist
1669  {
1670 #ifdef debug
1671  G4cout<<"-Warning-G4QContent::AddParton: Impossible Parton with PDG="<<pPDG<<G4endl;
1672 #endif
1673  return 0;
1674  }
1675  G4int HBN = GetBaryonNumber();
1676  if( HBN > 1 || HBN <-1)
1677  {
1678 #ifdef debug
1679  G4cout<<"-Warning-G4QContent::AddParton: Impossible Hadron with BaryonN="<<HBN<<G4endl;
1680 #endif
1681  return 0;
1682  }
1683  G4int AD=nAD;
1684  G4int AU=nAU;
1685  G4int AS=nAS;
1686  G4int QD=nD;
1687  G4int QU=nU;
1688  G4int QS=nS;
1689  if(aPDG>99) // Parton is DiQuark/antiDiQuark
1690  {
1691  G4int rPDG=aPDG/100;
1692  G4int P1=rPDG/10; // First quark
1693  G4int P2=rPDG%10; // Second quark
1694 #ifdef debug
1695  G4cout<<"G4QContent::AddParton: DiQuark/AntiDiQuark, P1="<<P1<<", P2="<<P2<<G4endl;
1696 #endif
1697  if(pPDG>0) // -- DiQuark
1698  {
1699 #ifdef debug
1700  G4cout<<"G4QContent::AddParton: DiQuark, P1="<<P1<<", P2="<<P2<<",HBN="<<HBN<<G4endl;
1701 #endif
1702  if (P1==3 && P2==3) // ----> ss DiQuark
1703  {
1704  if(HBN<0 && AS>1) AS-=2; // >> Annihilation of ss-DiQuark with anti-Baryon
1705  else if(!HBN && AS==1)
1706  {
1707  AS=0;
1708  ++QS;
1709  }
1710  else if(HBN || (!HBN && !AS)) return 0;
1711  }
1712  else if(P1==3 && P2==2) // ----> su DiQuark
1713  {
1714  if(HBN<0 && AS && AU) // >> Annihilation of su-DiQuark with anti-Baryon
1715  {
1716  --AS;
1717  --AU;
1718  }
1719  else if(!HBN && (AS || AU))
1720  {
1721  if(AS)
1722  {
1723  --AS;
1724  ++QU;
1725  }
1726  else
1727  {
1728  --AU;
1729  ++QS;
1730  }
1731  }
1732  else if(HBN || (!HBN && !AS && !AU)) return 0;
1733  }
1734  else if(P1==3 && P2==1) // ----> sd DiQuark
1735  {
1736  if(HBN<0 && AS && AD) // >> Annihilation of sd-DiQuark with anti-Baryon
1737  {
1738  --AS;
1739  --AD;
1740  }
1741  else if(!HBN && (AS || AD))
1742  {
1743  if(AS)
1744  {
1745  --AS;
1746  ++QD;
1747  }
1748  else
1749  {
1750  --AD;
1751  ++QS;
1752  }
1753  }
1754  else if(HBN || (!HBN && !AS && !AD)) return 0;
1755  }
1756  else if(P1==2 && P2==2) // ----> uu DiQuark
1757  {
1758  if(HBN<0 && AU>1) AU-=2; // >> Annihilation of uu-DiQuark with anti-Baryon
1759  else if(!HBN && AU==1)
1760  {
1761  AU=0;
1762  ++QU;
1763  }
1764  else if(HBN || (!HBN && !AU)) return 0;
1765  }
1766  else if(P1==2 && P2==1) // ----> ud DiQuark
1767  {
1768  if(HBN<0 && AD && AU) // >> Annihilation of ud-DiQuark with anti-Baryon
1769  {
1770  --AD;
1771  --AU;
1772  }
1773  else if(!HBN && (AD || AU))
1774  {
1775  if(AD)
1776  {
1777  --AD;
1778  ++QU;
1779  }
1780  else
1781  {
1782  --AU;
1783  ++QD;
1784  }
1785  }
1786  else if(HBN || (!HBN && !AU && !AD)) return 0;
1787  }
1788  else // ----> dd DiQuark
1789  {
1790  if(HBN<0 && AD>1) AD-=2; // >> Annihilation of dd-DiQuark with anti-Baryon
1791  else if(!HBN && AD==1)
1792  {
1793  AD=0;
1794  ++QD;
1795  }
1796  else if(HBN || (!HBN && !AD)) return 0;
1797  }
1798 #ifdef debug
1799  G4cout<<"G4QContent::AddParton: DQ, QC="<<QD<<","<<QU<<","<<QS<<","<<AD<<","<<AU<<","
1800  <<AS<<G4endl;
1801 #endif
1802  if (HBN<0) // ....... Hadron is an Anti-Baryon
1803  {
1804  if (AD) return -1; // ----->> Answer is anti-d
1805  else if(AU) return -2; // ----->> Answer is anti-u
1806  else return -3; // ----->> Answer is anti-s
1807  }
1808  else // ... Hadron is aMeson with annihilatedAntiQuark
1809  {
1810  if (QS) // --------- There is an s-quark
1811  {
1812  if (QS==2) return 3303; // ----->> Answer is ss (3301 does not exist)
1813  else if(QU) return 3201; // ----->> Answer is su (@@ only lightest)
1814  else return 3101; // ----->> Answer is sd (@@ only lightest)
1815  }
1816  else if(QU) // --------- There is an u quark
1817  {
1818  if (QU==2) return 2203; // ----->> Answer is uu (2201 does not exist)
1819  else return 2101; // ----->> Answer is ud (@@ only lightest)
1820  }
1821  else return 1103; // ----->> Answer is dd (1101 does not exist)
1822  }
1823  }
1824  else // -- antiDiQuark
1825  {
1826 #ifdef debug
1827  G4cout<<"G4QContent::AddParton: AntiDiQuark,P1="<<P1<<",P2="<<P2<<",B="<<HBN<<G4endl;
1828 #endif
1829  if (P1==3 && P2==3) // ----> anti-s-anti-s DiQuark
1830  {
1831  if(HBN>0 && QS>1) QS-=2; // >> Annihilation of anti-ss-DiQuark with Baryon
1832  else if(!HBN && QS==1)
1833  {
1834  QS=0;
1835  ++AS;
1836  }
1837  else if(HBN || (!HBN && !QS)) return 0;
1838  }
1839  else if(P1==3 && P2==2) // ----> anti-s-anti-u DiQuark
1840  {
1841  if(HBN>0 && QS && QU) // >> Annihilation of anti-su-DiQuark with Baryon
1842  {
1843  --QS;
1844  --QU;
1845  }
1846  else if(!HBN && (QS || QU))
1847  {
1848  if(QS)
1849  {
1850  --QS;
1851  ++AU;
1852  }
1853  else
1854  {
1855  --QU;
1856  ++AS;
1857  }
1858  }
1859  else if(HBN || (!HBN && !QS && !QU)) return 0;
1860  }
1861  else if(P1==3 && P2==1) // ----> anti-s-anti-d DiQuark
1862  {
1863  if(HBN>0 && QS && QD) // >> Annihilation of anti-sd-DiQuark with Baryon
1864  {
1865  --QS;
1866  --QD;
1867  }
1868  else if(!HBN && (QS || QD))
1869  {
1870  if(QS)
1871  {
1872  --QS;
1873  ++AD;
1874  }
1875  else
1876  {
1877  --QD;
1878  ++AS;
1879  }
1880  }
1881  else if(HBN || (!HBN && !QS && !QD)) return 0;
1882  }
1883  else if(P1==2 && P2==2) // ----> anti-u-anti-u DiQuark
1884  {
1885  if(HBN>0 && QU>1) QU-=2; // >> Annihilation of anti-uu-DiQuark with Baryon
1886  else if(!HBN && QU==1)
1887  {
1888  QU=0;
1889  ++AU;
1890  }
1891  else if(HBN || (!HBN && !QU)) return 0;
1892  }
1893  else if(P1==2 && P2==1) // ----> anti-u-anti-d DiQuark
1894  {
1895  if(HBN>0 && QU && QD) // >> Annihilation of anti-ud-DiQuark with Baryon
1896  {
1897  --QU;
1898  --QD;
1899  }
1900  else if(!HBN && (QU || QD))
1901  {
1902  if(QU)
1903  {
1904  --QU;
1905  ++AD;
1906  }
1907  else
1908  {
1909  --QD;
1910  ++AU;
1911  }
1912  }
1913  else if(HBN || (!HBN && !QU && !QD)) return 0;
1914  }
1915  else // ----> anti-d=anti-d DiQuark
1916  {
1917  if(HBN>0 && QD>1) QD-=2; // >> Annihilation of anti-dd-DiQuark with Baryon
1918  else if(!HBN && QD==1)
1919  {
1920  QD=0;
1921  ++AD;
1922  }
1923  else if(HBN || (!HBN && !QD)) return 0;
1924  }
1925 #ifdef debug
1926  G4cout<<"G4QContent::AddParton:ADQ, QC="<<QD<<","<<QU<<","<<QS<<","<<AD<<","<<AU<<","
1927  <<AS<<G4endl;
1928 #endif
1929  if (HBN>0) // ....... Hadron is an Baryon
1930  {
1931  if (QD) return 1; // ----->> Answer is d
1932  else if(QU) return 2; // ----->> Answer is u
1933  else return 3; // ----->> Answer is s
1934  }
1935  else // ....... Meson with annihilated Anti-Quark
1936  {
1937  if (AS) // --------- There is an anti-s quark
1938  {
1939  if (AS==2) return -3303; // ----->> Answer is anti-ss (3301 does not exist)
1940  else if(AU) return -3201; // ----->> Answer is anti-su (@@ only lightest)
1941  else return -3101; // ----->> Answer is anti-sd (@@ only lightest)
1942  }
1943  else if(AU) // --------- There is an anti-u quark
1944  {
1945  if (AU==2) return -2203; // ----->> Answer is anti-uu (2201 does not exist)
1946  else return -2101; // ----->> Answer is anti-ud (@@ only lightest)
1947  }
1948  else return -1103; // ----->> Answer is anti-dd (1101 does not exist)
1949  }
1950  }
1951  }
1952  else // Parton is Quark/antiQuark
1953  {
1954  if(pPDG>0) // -- Quark
1955  {
1956 #ifdef debug
1957  G4cout<<"G4QContent::AddParton: Quark, A="<<AD<<","<<AU<<","<<AS<<",B="<<HBN<<G4endl;
1958 #endif
1959  if (aPDG==1) // ----> d quark
1960  {
1961  if(HBN<0 && AD) AD--; // ****> Annihilation of d-quark with anti-Baryon
1962  else if(HBN || (!HBN && !AD)) return 0;
1963  }
1964  else if(aPDG==2) // ----> u quark
1965  {
1966  if(HBN<0 && AU) AU--; // ****> Annihilation of u-quark with anti-Baryon
1967  else if(HBN || (!HBN && !AU)) return 0;
1968  }
1969  else // ----> s quark
1970  {
1971  if(HBN<0 && AS) AS--; // ****> Annihilation of s-quark with anti-Baryon
1972  else if(HBN || (!HBN && !AS)) return 0;
1973  }
1974 #ifdef debug
1975  G4cout<<"G4QContent::AddParton: Q, QC="<<QD<<","<<QU<<","<<QS<<","<<AD<<","<<AU<<","
1976  <<AS<<G4endl;
1977 #endif
1978  if (!HBN) // ....... Hadron is a Meson (passingThrougAbove)
1979  {
1980  if (QD) return 1; // ----->> Answer is d
1981  else if(QU) return 2; // ----->> Answer is u
1982  else return 3; // ----->> Answer is s
1983  }
1984  else // ....... AntiBaryon with annihilated AntiQuark
1985  {
1986  if (AS) // --------- There is an anti-s quark
1987  {
1988  if (AS==2) return -3303; // ----->> Answer is ss (3301 does not exist)
1989  else if(AU) return -3201; // ----->> Answer is su (@@ only lightest)
1990  else return -3101; // ----->> Answer is sd (@@ only lightest)
1991  }
1992  else if(AU)
1993  {
1994  if (AU==2) return -2203; // ----->> Answer is uu (2201 does not exist)
1995  else return -2101; // ----->> Answer is ud (@@ only lightest)
1996  }
1997  else return -1103; // ----->> Answer is dd (1101 does not exist)
1998  }
1999  }
2000  else // -- antiQuark
2001  {
2002 #ifdef debug
2003  G4cout<<"G4QContent::AddParton: antiQ, Q="<<QD<<","<<QU<<","<<QS<<",B="<<HBN<<G4endl;
2004 #endif
2005  if (aPDG==1) // ---->> anti-d quark
2006  {
2007  if(HBN>0 && QD) QD--; // ****> Annihilation of anti-d-quark with Baryon
2008  else if(HBN || (!HBN && !QD)) return 0;
2009  }
2010  else if(aPDG==2) // ----> anti-u quark
2011  {
2012  if(HBN>0 && QU) QU--; // ****> Annihilation of anti-u-quark with Baryon
2013  else if(HBN || (!HBN && !QU)) return 0;
2014  }
2015  else // ----> anti-s quark
2016  {
2017  if(HBN>0 && QS) QS--; // ****> Annihilation of anti-s-quark with Baryon
2018  else if(HBN || (!HBN && !QS)) return 0;
2019  }
2020 #ifdef debug
2021  G4cout<<"G4QContent::AddParton: AQ, QC="<<QD<<","<<QU<<","<<QS<<","<<AD<<","<<AU<<","
2022  <<AS<<G4endl;
2023 #endif
2024  if (!HBN) // ....... Hadron is a Meson (passingThrougAbove)
2025  {
2026  if (AD) return -1; // ----->> Answer is anti-d
2027  else if(AU) return -2; // ----->> Answer is anti-u
2028  else return -3; // ----->> Answer is anti-s
2029  }
2030  else // ....... Baryon with annihilated Quark
2031  {
2032  if (QS) // --------- There is an anti-s quark
2033  {
2034  if (QS==2) return 3303; // ----->> Answer is ss (3301 does not exist)
2035  else if(QU) return 3201; // ----->> Answer is su (@@ only lightest)
2036  else return 3101; // ----->> Answer is sd (@@ only lightest)
2037  }
2038  else if(QU)
2039  {
2040  if (QU==2) return 2203; // ----->> Answer is uu (2201 does not exist)
2041  else return 2101; // ----->> Answer is ud (@@ only lightest)
2042  }
2043  else return 1103; // ----->> Answer is dd (1101 does not exist)
2044  }
2045  }
2046  }
2047 }
Generated on Sat May 25 2013 14:33:45 for Geant4 by   doxygen 1.8.4