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

Last change on this file since 1071 was 1055, checked in by garnier, 15 years ago

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