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

Last change on this file since 836 was 819, checked in by garnier, 16 years ago

import all except CVS
File size: 35.7 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.43.4.1 2008/04/23 14:47:22 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-01-patch-02 $
29//
30//      ---------------- G4QContent ----------------
31//             by Mikhail Kossov, Sept 1999.
32//      class for Quasmon initiated Contents used by CHIPS Model
33// --------------------------------------------------------------------
34// @@ In future total spin & c,b,t of the Hadron can be added @@ M.K.@@
35
36//#define debug
37//#define pdebug
38//#define ppdebug
39//#define erdebug
40
41#include "G4QContent.hh"
42#include <cmath>
43using namespace std;
44
45// Constructors
46G4QContent::G4QContent(G4int d, G4int u, G4int s, G4int ad, G4int au, G4int as):
47  nD(d),nU(u),nS(s),nAD(ad),nAU(au),nAS(as)
48{
49  if(d<0||u<0||s<0||ad<0||au<0||as<0)
50  {
51#ifdef erdebug
52    G4cerr<<"***G4QContent:"<<d<<","<<u<<","<<s<<","<<ad<<","<<au<<","<<as<<G4endl;
53    //throw G4QException("***G4QContent::Constructor: Negative Quark Content");
54#endif
55    if(d<0) ad-=d;
56    if(u<0) au-=u;
57    if(s<0) as-=s;
58    if(ad<0) d-=ad;
59    if(au<0) u-=au;
60    if(as<0) s-=as;
61  }
62}
63
64G4QContent::G4QContent(const G4QContent& right)
65{
66  nU  = right.nU;
67  nD  = right.nD;
68  nS  = right.nS;
69  nAU = right.nAU;
70  nAD = right.nAD;
71  nAS = right.nAS;
72}
73
74G4QContent::G4QContent(G4QContent* right)
75{
76  nU  = right->nU;
77  nD  = right->nD;
78  nS  = right->nS;
79  nAU = right->nAU;
80  nAD = right->nAD;
81  nAS = right->nAS;
82}
83
84// Assignment operator (copy stile for possible Vector extention)
85const G4QContent& G4QContent::operator=(const G4QContent &right)
86{//               ==============================================
87  if(this != &right)                          // Beware of self assignment
88  {
89    nU  = right.nU;
90    nD  = right.nD;
91    nS  = right.nS;
92    nAU = right.nAU;
93    nAD = right.nAD;
94    nAS = right.nAS;
95                }
96  return *this;
97}
98
99// Standard output for QC {d,u,s,ad,au,as}
100ostream& operator<<(ostream& lhs, G4QContent& rhs)
101{//      =========================================
102  lhs << "{" << rhs.GetD() << "," << rhs.GetU() << "," << rhs.GetS() << ","
103      << rhs.GetAD() << "," << rhs.GetAU() << "," << rhs.GetAS() << "}";
104  return lhs;
105}
106
107// Standard output for const QC {d,u,s,ad,au,as}
108ostream& operator<<(ostream& lhs, const G4QContent& rhs)
109{//      ===============================================
110  lhs << "{" << rhs.GetD() << "," << rhs.GetU() << "," << rhs.GetS() << ","
111      << rhs.GetAD() << "," << rhs.GetAU() << "," << rhs.GetAS() << "}";
112  return lhs;
113}
114
115// Subtract Quark Content
116G4QContent G4QContent::operator-=(const G4QContent& rhs)
117//         =============================================
118{
119#ifdef debug
120  G4cout<<"G4QC::-=(const): is called:"<<G4endl;
121#endif
122  G4int rD=rhs.nD;
123  G4int rU=rhs.nU;
124  G4int rS=rhs.nS;
125  G4int rAD=rhs.nAD;
126  G4int rAU=rhs.nAU;
127  G4int rAS=rhs.nAS;
128  ///////////G4int rQ =rD+rU+rS;
129  ///////////G4int rAQ=rAD+rAU+rAS;
130  ///////////G4int nQ =nD+nU+nS;
131  ///////////G4int nAQ=nAD+nAU+nAS;
132  if(nU<rU||nAU<rAU||nD<rD||nAD<rAD)
133  {
134    G4int dU=rU-nU;
135    G4int dAU=rAU-nAU;
136    if(dU>0||dAU>0)
137        {
138      G4int kU=dU;
139      if(kU<dAU) kU=dAU;                  // Get biggest difference
140      G4int mU=rU;
141      if(rAU<mU) mU=rAU;                  // Get a#of possible SS pairs
142      if(kU<=mU)                          // Total compensation
143             {
144        rU-=kU;
145        rAU-=kU;
146        rD+=kU;
147        rAD+=kU;
148      }
149      else                                // Partial compensation
150      {
151        rU-=mU;
152        rAU-=mU;
153        rD+=mU;
154        rAD+=mU;
155      }
156    }
157    G4int dD=rD-nD;
158    G4int dAD=rAD-nAD;
159    if(dD>0||dAD>0)
160           {
161      G4int kD=dD;
162      if(kD<dAD) kD=dAD;                  // Get biggest difference
163      G4int mD=rD;
164      if(rAD<mD) mD=rAD;                  // Get a#of possible SS pairs
165      if(kD<=mD)                          // Total compensation
166             {
167        rD-=kD;
168        rAD-=kD;
169        rU+=kD;
170        rAU+=kD;
171      }
172      else                                // Partial compensation
173      {
174        rD-=mD;
175        rAD-=mD;
176        rU+=mD;
177        rAU+=mD;
178      }
179    }
180  }
181#ifdef debug
182  G4cout<<"G4QC::-=:comp: "<<rD<<","<<rU<<","<<rS<<","<<rAD<<","<<rAU<<","<<rAS<<G4endl;
183#endif
184  if(rS==1 && rAS==1 && (nS<1 || nAS<1))  // Eta case, switch quark pairs (?)
185  {
186    rS =0;
187    rAS=0;
188    if(nU>rU&&nAU>rAU)
189    {
190      rU +=1;
191      rAU+=1;
192           }
193    else
194    {
195      rD +=1;
196      rAD+=1;
197           }
198  }
199  nD -= rD;
200  if (nD<0)
201  {
202    nAD -= nD;
203    nD   = 0;
204  }
205  nU -= rU;
206  if (nU<0)
207  {
208    nAU -= nU;
209    nU   = 0;
210  }
211  nS -= rS;
212  if (nS<0)
213  {
214    nAS -= nS;
215    nS   = 0;
216  }
217  nAD -= rAD;
218  if (nAD<0)
219  {
220    nD -= nAD;
221    nAD = 0;
222  }
223  nAU -= rAU;
224  if (nAU<0)
225  {
226    nU -= nAU;
227    nAU = 0;
228  }
229  nAS -= rAS;
230  if (nAS<0)
231  {
232    nS -= nAS;
233    nAS = 0;
234  }
235  return *this;
236} 
237
238// Subtract Quark Content
239G4QContent G4QContent::operator-=(G4QContent& rhs)
240//         =======================================
241{
242#ifdef debug
243  G4cout<<"G4QC::-=: is called:"<<G4endl;
244#endif
245  G4int rD=rhs.nD;
246  G4int rU=rhs.nU;
247  G4int rS=rhs.nS;
248  G4int rAD=rhs.nAD;
249  G4int rAU=rhs.nAU;
250  G4int rAS=rhs.nAS;
251  G4int rQ =rD+rU+rS;
252  G4int rAQ=rAD+rAU+rAS;
253  G4int nQ =nD+nU+nS;
254  G4int nAQ=nAD+nAU+nAS;
255  if(nQ<rQ||nAQ<rAQ)
256  {
257    G4int dU=rU-nU;
258    G4int dAU=rAU-nAU;
259    if(dU>0||dAU>0)
260           {
261      G4int kU=dU;
262      if(kU<dAU) kU=dAU;                  // Get biggest difference
263      G4int mU=rU;
264      if(rAU<mU) mU=rAU;                  // Get a#of possible SS pairs
265      if(kU<=mU)                          // Total compensation
266             {
267        rU-=kU;
268        rAU-=kU;
269        rD+=kU;
270        rAD+=kU;
271      }
272      else                                // Partial compensation
273      {
274        rU-=mU;
275        rAU-=mU;
276        rD+=mU;
277        rAD+=mU;
278      }
279    }
280    G4int dD=rD-nD;
281    G4int dAD=rAD-nAD;
282    if(dD>0||dAD>0)
283           {
284      G4int kD=dD;
285      if(kD<dAD) kD=dAD;                  // Get biggest difference
286      G4int mD=rD;
287      if(rAD<mD) mD=rAD;                  // Get a#of possible SS pairs
288      if(kD<=mD)                          // Total compensation
289             {
290        rD-=kD;
291        rAD-=kD;
292        rU+=kD;
293        rAU+=kD;
294      }
295      else                                // Partial compensation
296      {
297        rD-=mD;
298        rAD-=mD;
299        rU+=mD;
300        rAU+=mD;
301      }
302    }
303  }
304  if(rS==1 && rAS==1 && (nS<1 || nAS<1))  // Eta case, switch quark pairs (?)
305  {
306    rS =0;
307    rAS=0;
308    if(nU>rU&&nAU>rAU)
309    {
310      rU +=1;
311      rAU+=1;
312           }
313    else
314    {
315      rD +=1;
316      rAD+=1;
317           }
318  }
319  nD -= rD;
320  if (nD<0)
321  {
322    nAD -= nD;
323    nD   = 0;
324  }
325  nU -= rU;
326  if (nU<0)
327  {
328    nAU -= nU;
329    nU   = 0;
330  }
331  nS -= rS;
332  if (nS<0)
333  {
334    nAS -= nS;
335    nS   = 0;
336  }
337  nAD -= rAD;
338  if (nAD<0)
339  {
340    nD -= nAD;
341    nAD = 0;
342  }
343  nAU -= rAU;
344  if (nAU<0)
345  {
346    nU -= nAU;
347    nAU = 0;
348  }
349  nAS -= rAS;
350  if (nAS<0)
351  {
352    nS -= nAS;
353    nAS = 0;
354  }
355  return *this;
356} 
357
358// Overloading of QC addition
359G4QContent operator+(const G4QContent& lhs, const G4QContent& rhs)
360{//        =======================================================
361  G4QContent s  = lhs;
362  return     s += rhs;
363}
364
365// Overloading of QC subtraction
366G4QContent operator-(const G4QContent& lhs, const G4QContent& rhs)
367{//        =======================================================
368  G4QContent s  = lhs;
369  return     s -= rhs;
370}
371
372// Overloading of QC multiplication by Int
373G4QContent operator*(const G4QContent& lhs, const G4int&      rhs)
374{//        =======================================================
375  G4QContent s  = lhs;
376  return     s *= rhs;
377}
378
379// Overloading of Int times QC multiplication
380G4QContent operator*(const G4int&      lhs, const G4QContent& rhs)
381{//        =======================================================
382  G4QContent s  = rhs;
383  return     s *= lhs;
384}
385
386// Destructor - - - - - - -
387G4QContent::~G4QContent() {}
388
389// Subtract neutral pion from Quark Content (with possible hidden strangeness)
390G4bool G4QContent::SubtractPi0()
391{//    =========================
392#ifdef debug
393  G4cout<<"G4QC::SubtractPi0: U="<<nU<<", AU="<<nAU<<", D="<<nD<<", AD="<<nAD<<G4endl;
394#endif
395  G4int tot=GetTot();
396  G4int ab =GetBaryonNumber();
397  if(ab){if(tot<3*ab+2) return false;}
398  else if(tot<4) return false;
399
400  if(nU>0 && nAU>0)
401  {
402    nU--;
403    nAU--;
404    return true;
405  }
406  else if(nD>0 && nAD>0)
407  {
408    nD--;
409    nAD--;
410    return true;
411  }
412  return false;
413}
414
415// Subtract charged pion from Quark Content
416G4bool G4QContent::SubtractPion()
417{//    ==========================
418#ifdef debug
419  G4cout<<"G4QC::SubtractPion: U="<<nU<<", AU="<<nAU<<", D="<<nD<<", AD="<<nAD<<G4endl;
420#endif
421  G4int tot=GetTot();
422  G4int ab =GetBaryonNumber();
423  if(ab){if(tot<3*ab+2) return false;}
424  else if(tot<4) return false;
425
426  if((nU>nAU || (ab && nU>0))&& nAD>0)
427  {
428    nU--;
429    nAD--;
430    return true;
431  }
432  else if((nAU>nU || (ab && nAU>0)) && nD>0)
433  {
434    nAU--;
435    nD--;
436    return true;
437  }
438  return false;
439}
440
441// Subtract Hadron from Quark Content
442G4bool G4QContent::SubtractHadron(G4QContent h)
443{//    ========================================
444#ifdef debug
445  G4cout<<"G4QC::SubtractHadron "<<h<<" is called for QC="<<GetThis()<<G4endl;
446#endif
447  if(h.GetU()<=nU  && h.GetD()<=nD  && h.GetS()<=nS&&
448     h.GetAU()<=nAU&& h.GetAD()<=nAD&& h.GetAS()<=nAS) return true;
449  return false;
450}
451
452// Subtract Kaon from Quark Content
453G4bool G4QContent::SubtractKaon(G4double mQ)
454{//    =====================================
455#ifdef debug
456  G4cout<<"G4QC::SubtractKaon is called: QC="<<GetThis()<<G4endl;
457#endif
458  if(mQ<640.) return false;
459  G4int tot=GetTot();
460  G4int ab =GetBaryonNumber();
461  if(ab){if(tot<3*ab+2) return false;}
462  else if(tot<4) return false;
463
464  if((nS>nAS || (ab && nS>0)) && (nAD>0 || nAU>0))
465  {
466    nS--;
467    if (nAU>0) nAU--;
468    else       nAD--;
469    return true;
470  }
471  else if((nAS>nS || (ab && nAS>0)) && (nD>0 || nU>0))
472  {
473    nAS--;
474    if (nU>0) nU--;
475    else      nD--;
476    return true;
477  }
478#ifdef debug
479  G4cout<<"G4QCont::SubtractKaon Can't SubtractKaon: QC="<<GetThis()<<G4endl;
480#endif
481  return false;
482}
483
484// Split any hadronic system in two hadrons
485G4QContent G4QContent::SplitChipo (G4double mQ)
486{//        ====================================
487  G4QContent Pi(0,1,0,1,0,0);
488  if      (nU>0&&nAU>0) Pi=G4QContent(0,1,0,0,1,0);
489  else if (nD>0&&nAD>0) Pi=G4QContent(1,0,0,1,0,0);
490  else if (nD>=nU&&nAU>=nAD) Pi=G4QContent(1,0,0,0,1,0);
491  G4int bn=GetBaryonNumber();
492  G4int b =abs(bn);
493  if(!b && mQ<545.&&nS>0&&nAS>0) // Cancel strange sea
494  {
495    G4int      ss= nS;
496    if(nAS<nS) ss=nAS;
497    nS -=ss;
498    nAS-=ss;
499  }
500  if       (!b)DecQAQ(-4);
501  else if(b==1)DecQAQ(-5);
502  else         DecQAQ(0);
503  G4int tot=GetTot();
504  G4int q=GetQ();
505  G4int aq=GetAQ();
506  G4QContent r=Pi;     // Pion prototype of returned value
507  if((tot!=4||q!=2) && (tot!=5||(q!=1&&aq!=1)) && (tot!=6||abs(b)!=2))
508  {
509           //#ifdef erdebug
510    G4cerr<<"***G4QCont::SplitChipo: QC="<<GetThis()<<G4endl;
511           //#endif
512  }
513  else if(tot==4)      // Mesonic (eight possibilities)
514  {
515    r=GetThis();
516    if     (r.SubtractPi0())    ; // Try any trivial algorithm of splitting
517    else if(r.SubtractPion())   ;
518    else if(r.SubtractKaon(mQ)) ;
519    else
520    {
521             //#ifdef debug
522      G4cerr<<"***G4QCont::SplitChipo:MesTot="<<tot<<",b="<<b<<",q="<<q<<",a="<<aq<<G4endl;
523             //#endif
524    }
525  }
526  else if(b==1&&tot==5)      // Baryonic (four possibilities)
527  {
528    if(nU==3)
529           {
530      r.SetU(1);
531      r+=IndAQ();
532    }
533    else if(nD==3)
534           {
535      r.SetD(1);
536      r+=IndAQ();
537    }
538    else if(nS==3)
539           {
540      r.SetS(1);
541      r+=IndAQ();
542    }
543    else if(nAU==3)
544           {
545      r.SetAU(1);
546      r+=IndQ();
547    }
548    else if(nAD==3)
549           {
550      r.SetAD(1);
551      r+=IndQ();
552    }
553    else if(nAS==3)
554           {
555      r.SetAS(1);
556      r+=IndQ();
557    }
558    else if(q==1&&nU)
559           {
560      r.SetU(1);
561      if(nAU) r.SetAU(1);
562      else    r.SetAD(1);
563    }
564    else if(q==1&&nD)
565           {
566      r.SetD(1);
567      if(nAD) r.SetAD(1);
568      else    r.SetAU(1);
569    }
570    else if(q==1&&nS)
571           {
572      r.SetS(1);
573      if(nAS) r.SetAS(1);
574      else    r.SetAU(1);
575    }
576    else if(aq==1&&nAU)
577           {
578      r.SetAU(1);
579      if(nU) r.SetU(1);
580      else   r.SetD(1);
581    }
582    else if(aq==1&&nAD)
583           {
584      r.SetAD(1);
585      if(nD) r.SetD(1);
586      else   r.SetU(1);
587    }
588    else if(aq==1&&nAS)
589           {
590      r.SetAS(1);
591      if(nS) r.SetS(1);
592      else   r.SetU(1);
593    }
594    else
595    {
596             //#ifdef erdebug
597      G4cerr<<"***G4QCont::SplitChipo: Baryonic tot=5,b=1,qCont="<<GetThis()<<G4endl;
598      //#endif
599    }
600  }
601  else if(tot==b*3)          // MultyBaryon cace
602  {
603    r=GetThis();
604    if (bn>0)                // baryonium
605           {
606      G4QContent la(1,1,1,0,0,0);
607      G4QContent nt(2,1,0,0,0,0);
608      G4QContent pr(1,2,0,0,0,0);
609      G4QContent ks(0,1,2,0,0,0);
610      if (nD>nU) ks=G4QContent(1,0,2,0,0,0);
611      G4QContent dm(3,0,0,0,0,0);
612      G4QContent dp(0,3,0,0,0,0);
613      G4QContent om(0,0,3,0,0,0);
614      if     (nU>=nD&&nU>=nS)
615      {
616        if     (r.SubtractHadron(pr)) r-=pr; // These functions only check
617        else if(r.SubtractHadron(dp)) r-=dp;
618        else if(r.SubtractHadron(nt)) r-=nt;
619        else if(r.SubtractHadron(la)) r-=la;
620        else if(r.SubtractHadron(dm)) r-=dm;
621        else
622        {
623          //#ifdef erdebug
624          G4cerr<<"***G4QCont::SplitChipo:Dibar (1) tot=6, b=2, qCont="<<GetThis()<<G4endl;
625          //#endif
626        }
627      }
628      else if(nD>=nU&&nD>=nS)
629      {
630        if     (r.SubtractHadron(nt)) r-=nt; // These functions only check
631        else if(r.SubtractHadron(dm)) r-=dm;
632        else if(r.SubtractHadron(pr)) r-=pr;
633        else if(r.SubtractHadron(dp)) r-=dp;
634        else if(r.SubtractHadron(la)) r-=la;
635        else
636        {
637            //#ifdef erdebug
638          G4cerr<<"***G4QContent::SplitChipo:Dib(2) tot=6, b=2, qCont="<<GetThis()<<G4endl;
639          //#endif
640        }
641      }
642      else
643      {
644        if     (r.SubtractHadron(la)) r-=la; // These functions only check
645        else if(r.SubtractHadron(ks)) r-=ks;
646        else if(r.SubtractHadron(om)) r-=om;
647        else if(r.SubtractHadron(pr)) r-=pr;
648        else if(r.SubtractHadron(nt)) r-=nt;
649        else
650        {
651            //#ifdef erdebug
652          G4cerr<<"***G4QContent::SplitChipo:Dib(3) tot=6, b=2, qCont="<<GetThis()<<G4endl;
653          //#endif
654        }
655      }
656           }
657    else                     // Anti-baryonium
658           {
659      G4QContent la(0,0,0,1,1,1);
660      G4QContent pr(0,0,0,1,2,0);
661      G4QContent nt(0,0,0,2,1,0);
662      G4QContent ks(0,1,2,0,0,0);
663      if(nAD>nAU)ks=G4QContent(0,0,0,1,0,2);
664      G4QContent dm(0,0,0,3,0,0);
665      G4QContent dp(0,0,0,0,3,0);
666      G4QContent om(0,0,0,0,0,3);
667      if     (nAU>=nAD&&nAU>=nAS)
668      {
669        if     (r.SubtractHadron(pr)) r-=pr; // These functions only check
670        else if(r.SubtractHadron(dp)) r-=dp;
671        else if(r.SubtractHadron(nt)) r-=nt;
672        else if(r.SubtractHadron(la)) r-=la;
673        else if(r.SubtractHadron(dm)) r-=dm;
674        else
675        {
676            //#ifdef erdebug
677          G4cerr<<"***G4QContent::SplitChipo:ADib(1) tot=6,b=2, qCont="<<GetThis()<<G4endl;
678          //#endif
679        }
680      }
681      else if(nAD>=nAU&&nAD>=nAS)
682      {
683        if     (r.SubtractHadron(nt)) r-=nt; // These functions only check
684        else if(r.SubtractHadron(dm)) r-=dm;
685        else if(r.SubtractHadron(pr)) r-=pr;
686        else if(r.SubtractHadron(dp)) r-=dp;
687        else if(r.SubtractHadron(la)) r-=la;
688        else
689        {
690            //#ifdef erdebug
691          G4cerr<<"***G4QContent::SplitChipo:ADib(2) tot=6,b=2, qCont="<<GetThis()<<G4endl;
692          //#endif
693        }
694      }
695      else
696      {
697        if     (r.SubtractHadron(la)) r-=la; // These functions only check
698        else if(r.SubtractHadron(ks)) r-=ks;
699        else if(r.SubtractHadron(om)) r-=om;
700        else if(r.SubtractHadron(pr)) r-=pr;
701        else if(r.SubtractHadron(nt)) r-=nt;
702        else
703        {
704            //#ifdef erdebug
705          G4cerr<<"***G4QContent::SplitChipo:ADib(3) tot=6,b=2, qCont="<<GetThis()<<G4endl;
706          //#endif
707        }
708      }
709           }
710  }
711  else // More than Dibaryon (@@ can use the same algorithm as for dibaryon)
712  {
713    //#ifdef erdebug
714    G4cerr<<"*G4QContent::SplitChipolino:UnknownHadron with QuarkCont="<<GetThis()<<G4endl;
715    //#endif
716  }
717  return r;
718}// End of G4QContent::SplitChipolino
719
720// Return one-quark QC using index (a kind of iterator)
721G4QContent G4QContent::IndQ (G4int index)
722{//        ==============================
723#ifdef debug
724  G4cout << "G4QC::IndQ is called"<<G4endl;
725#endif
726  if(index<nD) return G4QContent(1,0,0,0,0,0);
727  else if(index<nD+nU) return G4QContent(0,1,0,0,0,0);
728  else if(index<nD+nU+nS) return G4QContent(0,0,1,0,0,0);
729  //#ifdef erdebug
730  else G4cerr<<"***G4QC::IndQ:index="<<index<<" for the QuarkContent="<<GetThis()<<G4endl;
731  //throw G4QException("***G4QC::IndQ: Index exceeds the total number of quarks");
732  //#endif
733  return G4QContent(0,0,0,0,0,0);
734}
735
736// Return one-antiquark QC using index (a kind of iterator)
737G4QContent G4QContent::IndAQ (G4int index)
738{//        ==============================
739#ifdef debug
740  G4cout << "G4QC::IndAQ is called"<<G4endl;
741#endif
742  if(index<nAD) return G4QContent(0,0,0,1,0,0);
743  else if(index<nAD+nAU) return G4QContent(0,0,0,0,1,0);
744  else if(index<nAD+nAU+nAS) return G4QContent(0,0,0,0,0,1);
745  //#ifdef erdebug
746  else G4cerr<<"***G4QC::IndAQ:index="<<index<<" for the QuarkContent="<<GetThis()<<G4endl;
747  //throw G4QException("***G4QC::IndAQ: Index exceeds the total number of antiquarks");
748  //#endif
749  return G4QContent(0,0,0,0,0,0);
750}
751
752// Reduces number (if nQAQ<0:all) of valence Q-Qbar pairs, returns #of pairs to reduce more
753G4int G4QContent::DecQAQ(const G4int& nQAQ)
754{//   =====================================
755#ifdef pdebug
756  G4cout<<"G4QCont::DecQC: n="<<nQAQ<<","<<GetThis()<<G4endl;
757#endif
758  G4int ban = GetBaryonNumber();
759  G4int tot = GetTot();    // Total number of quarks in QC
760  if (tot==ban*3) return 0;// Nothing to reduce & nothing to reduce more
761  G4int nUP=0;             // U/AU min factor (a#of u which can be subtracted)
762  if (nU>=nAU) nUP=nAU;
763  else         nUP= nU;
764
765  G4int nDP=0;             // D/AD min factor (a#of d which can be subtracted)
766  if (nD>=nAD) nDP=nAD;
767  else         nDP= nD;
768
769  G4int nSP=0;             // S/AS min factor (a#of s which can be subtracted)
770  if (nS>=nAS) nSP=nAS;
771  else         nSP= nS;
772
773  G4int nLP  =nUP+nDP;     // a#of light quark pairs which can be subtracted
774  G4int nTotP=nLP+nSP;     // total#of existing pairs which can be subtracted
775  G4int nReal=nQAQ;        // demanded #of pairs for reduction (by demand)
776  G4int nRet =nTotP-nQAQ;  // a#of additional pairs for further reduction
777  if (nQAQ<0)              // === Limited reduction case @@ not tuned for baryons !!
778  {
779    G4int res=tot+nQAQ;
780#ifdef pdebug
781        G4cout<<"G4QC::DecQC: tot="<<tot<<", nTP="<<nTotP<<", res="<<res<<G4endl;
782#endif
783    if(res<0)
784    {
785      IncQAQ(1,0.);        // Increment by one not strange pair to get the minimum
786      return 1;
787    }
788    res -=nTotP+nTotP;
789    if(res<0) nReal=nTotP+res/2;
790    else nReal=nTotP;
791    nRet = tot/2-nReal;
792  }
793  else if(!nQAQ)
794  {
795    nReal=nTotP;
796    nRet =0;
797  }
798  else if(nRet<0) nReal=nTotP;
799
800  if (!nReal) return nRet; // Now nothing to be done
801  // ---------- Decrimenting by nReal pairs
802#ifdef pdebug
803  G4cout<<"G4QC::DecQC: demanded "<<nQAQ<<" pairs, executed "<<nReal<<" pairs"<<G4endl;
804#endif
805  ///////////G4int nt = tot - nTotP - nTotP;
806  for (G4int i=0; i<nReal; i++)
807  {
808    G4double base = nTotP;
809    //if (nRet && nSP==1 && !nQAQ) base = nLP;              // Keep S-Sbar pair if possible
810    G4int j = static_cast<int>(base*G4UniformRand());       // Random integer "SortOfQuark"
811    if (nUP && j<nUP && (nRet>2 || nUP>1 || (nD<2 && nS<2)))// --- U-Ubar pair
812           {
813#ifdef pdebug
814      G4cout<<"G4QC::DecQC: decrementing UAU pair UP="<<nUP<<", QC="<<GetThis()<<G4endl;
815#endif
816      nU--;
817      nAU--;
818      nUP--;
819      nLP--;
820      nTotP--;
821           } 
822    else if (nDP && j<nLP && (nRet>2 || nDP>1 || (nU<2 && nS<2)))// --- D-Ubar pair
823           {
824#ifdef pdebug
825      G4cout<<"G4QC::DecQC: decrementing DAD pair DP="<<nDP<<", QC="<<GetThis()<<G4endl;
826#endif
827      nD--;
828      nAD--;
829      nDP--;
830      nLP--;
831      nTotP--;
832           } 
833    else if (nSP&& (nRet>2 || nSP>1 || (nU<2 && nD<2)))          // --- S-Sbar pair
834           {
835#ifdef pdebug
836      G4cout<<"G4QC::DecQC: decrementing SAS pair SP="<<nSP<<", QC="<<GetThis()<<G4endl;
837#endif
838      nS--;
839      nAS--;
840      nSP--;
841      nTotP--;
842           }
843    else if (nUP)                                  // --- U-Ubar pair cancelation (final)
844           {
845#ifdef pdebug
846      G4cout<<"G4QC::DecQC:Decrement UAU pair (final) UP="<<nUP<<",QC="<<GetThis()<<G4endl;
847#endif
848      nU--;
849      nAU--;
850      nUP--;
851      nLP--;
852      nTotP--;
853           } 
854    else if (nDP)                                 // --- D-Ubar pair cancelation (final)
855           {
856#ifdef pdebug
857      G4cout<<"G4QC::DecQC:Decrement DAD pair (final) DP="<<nDP<<",QC="<<GetThis()<<G4endl;
858#endif
859      nD--;
860      nAD--;
861      nDP--;
862      nLP--;
863      nTotP--;
864           } 
865    else if (nSP)                                 // --- S-Sbar pair cancelation (final)
866           {
867#ifdef pdebug
868      G4cout<<"G4QC::DecQC: decrementing SAS pair SP="<<nSP<<", QC="<<GetThis()<<G4endl;
869#endif
870      nS--;
871      nAS--;
872      nSP--;
873      nTotP--;
874           }
875    else G4cout<<"***G4QC::DecQC:i="<<i<<",j="<<j<<",D="<<nDP<<",U="<<nUP<<",S="<<nSP
876               <<",T="<<nTotP<<",nRet="<<nRet<<", QC="<<GetThis()<<G4endl;
877  }
878#ifdef pdebug
879  G4cout<<"G4QC::DecQC: >>>OUT<<< nRet="<<nRet<<", QC="<<GetThis()<<G4endl;
880#endif
881  return nRet;
882}
883
884// Increment quark pairs
885void G4QContent::IncQAQ(const G4int& nQAQ, const G4double& sProb)
886{//  ============================================================
887  G4int tot = GetTot();
888  G4QContent mQC = GetThis();
889  for (int i=0; i<nQAQ; i++)
890  {
891    G4int j = static_cast<int>((2.+sProb)*G4UniformRand()); // 0-U, 1-D, 2-S
892#ifdef debug
893    G4cout<<"IncQC:out QC="<<GetThis()<<",j="<<j<<" for i="<<i<<G4endl;
894#endif
895    //if      (!j)
896           if      ( !j && (nU<=nD || nU<=nS))
897    {
898      nU++;
899      nAU++;
900      tot+=2;
901           }
902    //else if (j==1)
903           else if (j==1 && (nD<=nU || nD<=nS))
904           {
905      nD++;
906      nAD++;
907      tot+=2;
908           }
909    //else
910           else if (j>1&& (nS<=nU || nS<=nD))
911    {
912      nS++;
913      nAS++;
914      tot+=2;
915           }
916    else if (!j)
917           {
918      nD++;
919      nAD++;
920      tot+=2;
921           }
922    else if (j==1)
923           {
924      nU++;
925      nAU++;
926      tot+=2;
927    }     
928    else
929    {
930      nS++;
931      nAS++;
932      tot+=2;
933           }
934    //else if (nD<=nU)
935           //{
936    //  nD++;
937    //  nAD++;
938    //  tot+=2;
939           //}
940    //else
941           //{
942    //  nU++;
943    //  nAU++;
944    //  tot+=2;
945    //}     
946  }
947}
948
949// Calculate a#of protons in the QC (for nuclei)
950G4int G4QContent::GetP() const
951{//   ========================
952  G4int rD=nD-nAD;                                   // Constituent d-quarks
953  G4int rU=nU-nAU;                                   // Constituent u-quarks
954  G4int rS=nS-nAS;                                   // Constituent s-quarks
955  G4int dQ=rD-rU;                                    // Isotopic assimetry
956  G4int b3=rD+rU+rS;                                 // (Baryon number) * 3
957  return (b3-3*(rS+dQ))/6;
958}
959
960// Calculate a#of neutrons in the QC (for nuclei)
961G4int G4QContent::GetN() const
962{//   ========================
963  G4int rD=nD-nAD;                                   // Constituent d-quarks
964  G4int rU=nU-nAU;                                   // Constituent u-quarks
965  G4int rS=nS-nAS;                                   // Constituent s-quarks
966  G4int dQ=rD-rU;                                    // Isotopic assimetry
967  G4int b3=rD+rU+rS;                                 // (Baryon number) * 3
968  return (b3-3*(rS-dQ))/6;
969}
970
971// Calculate a#of lambdas in the QC (for nuclei)
972G4int G4QContent::GetL() const
973{//   ========================
974  G4int rS=nS-nAS;                                   // Constituent s-quarks
975  return rS;
976}
977
978// Calculate a#of anti-protons in the QC (for anti-nuclei)
979G4int G4QContent::GetAP() const
980{//   =========================
981  G4int rD=nAD-nD;                                   // Constituent anti-d-quarks
982  G4int rU=nAU-nU;                                   // Constituent anti-u-quarks
983  G4int rS=nAS-nS;                                   // Constituent anti-s-quarks
984  G4int dQ=rD-rU;                                    // Isotopic assimetry
985  G4int b3=rD+rU+rS;                                 // - (Baryon number) * 3
986  return (b3-3*(rS+dQ))/6;
987}
988
989// Calculate a#of anti-neutrons in the QC (for anti-nuclei)
990G4int G4QContent::GetAN() const
991{//   =========================
992  G4int rD=nAD-nD;                                   // Constituent anti-d-quarks
993  G4int rU=nAU-nU;                                   // Constituent anti-u-quarks
994  G4int rS=nAS-nS;                                   // Constituent anti-s-quarks
995  G4int dQ=rD-rU;                                    // Isotopic assimetry
996  G4int b3=rD+rU+rS;                                 // - (Baryon number) * 3
997  return (b3-3*(rS-dQ))/6;
998}
999
1000// Calculate a#of anti-lambdas in the QC (for anti-nuclei)
1001G4int G4QContent::GetAL() const
1002{//   =========================
1003  G4int rS=nAS-nS;                                   // Constituent anti-s-quarks
1004  return rS;
1005}
1006
1007// Calculate charge for the QC
1008G4int G4QContent::GetCharge() const
1009{//   =============================
1010  static const G4int cU  = 2;
1011  static const G4int cD  =-1;
1012  static const G4int cS  =-1;
1013  static const G4int cAU =-2;
1014  static const G4int cAD = 1;
1015  static const G4int cAS = 1;
1016
1017  G4int c=0;
1018  if(nU) c+=nU*cU;
1019  if(nD) c+=nD*cD;
1020  if(nS) c+=nS*cS;
1021  if(nAU)c+=nAU*cAU;
1022  if(nAD)c+=nAD*cAD;
1023  if(nAS)c+=nAS*cAS;
1024  //#ifdef erdebug
1025  if(c%3) G4cerr<<"***G4QCont:GetCharge:c="<<c<<"/3 isn't integer, QC="<<GetThis()<<G4endl;
1026  //#endif
1027  return c/3;
1028}
1029
1030// Calculate a Baryon Number for the QC
1031G4int G4QContent::GetBaryonNumber() const
1032{//   ===================================
1033  G4int b=nU+nD+nS-nAU-nAD-nAS;
1034  //#ifdef erdebug
1035  if(b%3) G4cerr<<"-Warn-G4QContent:BaryonNumber="<<b<<"/3 isn't an integer value"<<G4endl;
1036  //#endif
1037  return b/3;
1038}
1039
1040// Gives the PDG of the lowest (in mass) S-wave hadron or Chipolino (=10) for double hadron
1041G4int G4QContent::GetSPDGCode() const
1042{//   ===============================
1043  G4int p = 0;           // Prototype of output SPDG
1044  G4int n = GetTot();    // Total number of quarks
1045  if(!n) return 22;      // Photon does not have any Quark Content
1046  G4int mD=nD;           // A # of D quarks or anti U quarks
1047  if (nD<=0) mD=nAD;
1048  G4int mU=nU;           // A # of U quarks or anti U quarks
1049  if (nU<=0) mU=nAU;
1050  G4int mS=nS;           // A # of S quarks or anti U quarks
1051  if (nS<=0) mS= nAS;
1052  // ---------------------- Cancelation of q-qbar pairs in case of an excess
1053  if ( nU>nAU && nAU>0)
1054  {
1055    mU=nU-nAU;
1056    n-=nAU+nAU;
1057  }
1058  if (nAU>nU  &&  nU>0)
1059  {
1060    mU=nAU-nU;
1061    n-=nU+nU;
1062  }
1063  if ( nD>nAD && nAD>0)
1064  {
1065    mD=nD-nAD;
1066    n-=nAD+nAD;
1067  }
1068  if (nAD>nD  &&  nD>0)
1069  {
1070    mD=nAD-nD;
1071    n-=nD+nD;
1072  }
1073  if ( nS>nAS && nAS>0)
1074  {
1075    mS=nS-nAS;
1076    n-=nAS+nAS;
1077  }
1078  if (nAS>nS  &&  nS>0)
1079  {
1080    mS= nAS-nS;
1081    n-=nS+nS;
1082  }
1083  // ---------------------- Cancelation of q-qbar pairs in case of an equality
1084  if (nAD==nD  &&  nD>0)
1085  {
1086    G4int dD=nD+nD;
1087    if(n>dD)
1088          {
1089      mD=0;
1090      n-=dD;
1091           }
1092    else if (n==dD)
1093          {
1094      mD=2;
1095      n=2;
1096           }
1097    else
1098    {
1099#ifdef pdebug
1100      G4cout<<"***G4QC::SPDG:CanD U="<<mU<<",D="<<mD<<",S="<<mS<<",QC="<<GetThis()<<G4endl;
1101#endif
1102      return 0;
1103    }
1104  }
1105  if (nAU==nU  &&  nU>0)
1106  {
1107    G4int dU=nU+nU;
1108    if(n>dU)
1109           {
1110      mU=0;
1111      n-=dU;
1112           }
1113    else if (n==dU)
1114          {
1115      mU=2;
1116      n=2;
1117           }
1118    else
1119    {
1120#ifdef pdebug
1121      G4cout<<"***G4QC::SPDG:CanU U="<<mU<<",D="<<mD<<",S="<<mS<<",QC="<<GetThis()<<G4endl;
1122#endif
1123      return 0;
1124    }
1125  }
1126  if (nAS==nS  &&  nS>0) //@@ Starts with S-quarks - should be randomized and mass limited
1127  {
1128    G4int dS=nS+nS;
1129    if(n>dS)
1130           {
1131      mS=0;
1132      n-=dS;
1133           }
1134    else if (n==dS)
1135          {
1136      mS=2;
1137      n=2;
1138           }
1139    else
1140    {
1141#ifdef pdebug
1142      G4cout<<"***G4QC::SPDG:CanS U="<<mU<<",D="<<mD<<",S="<<mS<<",QC="<<GetThis()<<G4endl;
1143#endif
1144      return 0;
1145    }
1146  }
1147
1148  G4int b=GetBaryonNumber();
1149  G4int c=GetCharge();
1150  G4int s=GetStrangeness();
1151#ifdef pdebug
1152  G4cout<<"G4QC::SPDGC:bef. b="<<b<<",n="<<n<<",c="<<c<<",s="<<s<<",Q="<<GetThis()<<G4endl;
1153#endif
1154  if (b)                                         // ==================== Baryon case
1155  {
1156    G4int ab=abs(b);
1157    if(ab>=2 && n>=6)                            // Multi-Baryonium (NuclearFragment)
1158           {
1159      G4int mI=nU-nAU-nD+nAD;
1160      //if     (abs(mI)>3||mS>3||(b>0&&s<-1)||(b<0&&s>1)) return  0;
1161      //else if(abs(mI)>2||mS>2||(b>0&&s< 0)||(b<0&&s>0)) return 10;
1162      if ( (b > 0 && s == -1) || (b < 0 && s == 1) ) return 10;
1163      else if (abs(mI) > 2 || mS > 2 
1164                           || (b > 0 && s < 0) 
1165                           || (b < 0 && s > 0)) return GetZNSPDGCode();
1166      else if(mU>=mS&&mD>=mS&&mU+mD+mS==3*b)     // Possible Unary Nuclear Cluster
1167             {
1168        G4int mZ=(mU+mD-mS-mS+3*mI)/6;
1169        p = 90000000+1000*(1000*mS+mZ)+mZ-mI;
1170        if(b>0) return  p;
1171        else    return -p;
1172             }
1173      else return 10;
1174           }
1175    // Normal One Baryon States: Heavy quark should come first
1176    if(n>5) return GetZNSPDGCode();            //B+M+M Tripolino etc
1177    if(n==5) return 10;                        //B+M Chipolino
1178    if(mS>0)                                   // Strange Baryons
1179           {
1180      p=3002;
1181      if      (mS==3)            p+=332;       // Decuplet
1182      else if (mS==2)
1183      {
1184        if      (mU==1 && mD==0) p+=320;
1185        else if (mU==0 && mD==1) p+=310;
1186        else
1187                      {
1188#ifdef debug
1189          G4cout<<"**G4QC::SPDG:ExoticBSS,U="<<mU<<",D="<<mD<<",S="<<mS<<GetThis()<<G4endl;
1190#endif
1191          return GetZNSPDGCode();
1192        }
1193             }
1194      else if (mS==1)
1195      {
1196        if      (mU==2 && mD==0) p+=220;
1197        else if (mU==1 && mD==1) p+=120;       // Lambda (M_Lambda<M_Sigma0) PDG_Sigma=3212
1198        else if (mU==0 && mD==2) p+=110;
1199        else
1200        {
1201#ifdef debug
1202          G4cout<<"***G4QC::SPDG:ExoticBS,U="<<mU<<",D="<<mD<<",S="<<mS<<GetThis()<<G4endl;
1203#endif
1204          return GetZNSPDGCode();
1205        }
1206             }
1207      else                                     // Superstrange case
1208      {
1209#ifdef debug
1210        G4cout<<"***G4QC::GetSPDG:ExoBarS,U="<<mU<<",D="<<mD<<",S="<<mS<<GetThis()<<G4endl;
1211#endif
1212        return GetZNSPDGCode();
1213      }
1214           }
1215    else if (mU>0)                               // Not Strange Baryons
1216           {
1217      p=2002;
1218      if      (mU==3 && mD==0) p+=222;           // Decuplet
1219      else if (mU==2 && mD==1) p+=210;
1220      else if (mU==1 && mD==2) p+=110;           // There is a higher Delta S31 (PDG=1212)
1221      else
1222      {
1223#ifdef debug
1224        G4cout<<"***G4QC::SPDG:ExoBaryonU,U="<<mU<<",D="<<mD<<",S="<<mS<<GetThis()<<G4endl;
1225#endif
1226        return GetZNSPDGCode();
1227      }
1228           }
1229    else if (mD==3) p=1114;                      // Decuplet
1230    else
1231    {
1232#ifdef debug
1233      G4cout<<"**G4QC::SPDG:ExoticBaryonD,U="<<mU<<",D="<<mD<<",S="<<mS<<GetThis()<<G4endl;
1234#endif
1235      return GetZNSPDGCode();
1236           }
1237    if (b<0) p=-p;
1238  }
1239  else             // ====================>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> Meson case
1240  {
1241#ifdef debug
1242    G4cout<<"G4QC::SPDG:mDUS="<<mD<<","<<mU<<","<<mS<<",b,c,s="<<b<<","<<c<<","<<s<<G4endl;
1243#endif
1244    if(n>4)                                      // Super Exotics
1245    {
1246#ifdef debug
1247      G4cout<<"G4QC::SPDG:n>4 SEx:U="<<mU<<",D="<<mD<<",S="<<mS<<",QC="<<GetThis()<<G4endl;
1248#endif
1249      return 0;
1250           }
1251    if(n==4) return 10;                          // M+M Chipolino
1252    if(abs(s)>1)
1253           {
1254#ifdef debug
1255      G4cout<<"**G4QC::SPDG:Stran="<<s<<",QC="<<GetThis()<<" - Superstrange Meson"<<G4endl;
1256#endif
1257      return 0;
1258           }
1259    // Heavy quark should come first
1260    if(mS>0)                                     // Strange Mesons
1261           {
1262      p=301;
1263      if      (mS==2)
1264      {
1265        //if (G4UniformRand()<0.333) p=221;        // eta
1266        //else                       p+=30;        // eta'
1267        p=221;
1268      }
1269      else if (mU==1 && mD==0) p+=20;
1270      else if (mU==0 && mD==1) p+=10;
1271      else
1272      {
1273#ifdef debug
1274        G4cout<<"*G4QC::SPDG:ExMS U="<<mU<<",D="<<mD<<",S="<<mS<<",QC="<<GetThis()<<G4endl;
1275#endif
1276        return 0;
1277      }
1278           }
1279    else if (mU>0)                               // Isotopic Mesons
1280           {
1281      p=201;
1282      //if      (mU==2 && mD==0) p=221; // Performance problems
1283      if      (mU==2 && mD==0) p=111;
1284      else if (mU==1 && mD==1) p+=10;
1285      else
1286             {
1287#ifdef debug
1288        G4cout<<"*G4QC::SPDG:ExMU U="<<mU<<",D="<<mD<<",S="<<mS<<",QC="<<GetThis()<<G4endl;
1289#endif
1290        return 0;
1291      }
1292           }
1293    else if (mD==2) p=111;
1294    else
1295    {
1296#ifdef debug
1297      G4cout<<"***G4QC::SPDG:ExMD U="<<mU<<",D="<<mD<<",S="<<mS<<",QC="<<GetThis()<<G4endl;
1298#endif
1299      return 0;
1300    }
1301    if (c<0 || (c==0 && mS>0 && s>0)) p=-p;
1302  }
1303#ifdef debug
1304  G4cout<<"G4QC::GetSPDG:output SPDGcode="<<p<<" for the QuarkContent="<<GetThis()<<G4endl;
1305#endif
1306  return p;
1307}
1308
1309// === Calculate a number of combinations of rhc out of lhc ==
1310G4int G4QContent::NOfCombinations(const G4QContent& rhs) const
1311{//   ========================================================
1312  G4int c=1; // Default number of combinations?
1313#ifdef ppdebug
1314  G4cout<<"G4QContent::NOfComb: This="<<GetThis()<<", selectQC="<<rhs<<G4endl;
1315#endif
1316  G4int mD=rhs.GetD();
1317  G4int mU=rhs.GetU();
1318  G4int mS=rhs.GetS();
1319  G4int mAD=rhs.GetAD();
1320  G4int mAU=rhs.GetAU();
1321  G4int mAS=rhs.GetAS();
1322  G4int mN=mD+mU+mS-mAD-mAU-mAS;
1323  ////////////G4int PDG=abs(GetSPDGCode());
1324  if (( ((nD < mD || nAD < mAD) && !(mD-mAD)) || 
1325        ((nU < mU || nAU < mAU) && !(mU-mAU)) ||
1326        ((nS < mS || nAS < mAS) && !(mS-mAS)) ) && !mN) return 1;
1327  if(mD>0)
1328  {
1329    int j=nD;
1330    if (j<=0) return 0;
1331    if(mD>1||j>1) for (int i=1; i<=mD; i++)
1332           {
1333      if(!j) return 0;
1334      c*=j/i;
1335      j--;
1336    }
1337  };
1338  if(mU>0)
1339  {
1340    int j=nU;
1341    if (j<=0) return 0;
1342    if(mU>1||j>1) for (int i=1; i<=mU; i++)
1343           {
1344      if(!j) return 0;
1345      c*=j/i;
1346      j--;
1347    }
1348  };
1349  if(mS>0)
1350  {
1351    int j=nS;
1352    if (j<=0) return 0;
1353    if(mS>1||j>1) for (int i=1; i<=mS; i++)
1354           {
1355      if(!j) return 0;
1356      c*=j/i;
1357      j--;
1358    }
1359  };
1360  if(mAD>0)
1361  {
1362    int j=nAD;
1363    if (j<=0) return 0;
1364    if(mAD>1||j>1) for (int i=1; i<=mAD; i++)
1365           {
1366      if(!j) return 0;
1367      c*=j/i;
1368      j--;
1369    }
1370  };
1371  if(mAU>0)
1372  {
1373    int j=nAU;
1374    if (j<=0) return 0;
1375    if(mAU>1||j>1) for (int i=1; i<=mAU; i++)
1376           {
1377      if(!j) return 0;
1378      c*=j/i;
1379      j--;
1380    }
1381  };
1382  if(mAS>0)
1383  {
1384    int j=nAS;
1385    if (j<=0) return 0;
1386    if(mAS>1||j>1) for (int i=1; i<=mAS; i++)
1387           {
1388      if(!j) return 0;
1389      c*=j/i;
1390      j--;
1391    }
1392  };
1393 
1394  return c;
1395} 
Note: See TracBrowser for help on using the repository browser.