source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/body/src/G4QContent.cc @ 1340

Last change on this file since 1340 was 1340, checked in by garnier, 14 years ago

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