Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4QChipolino.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id$
28 //
29 // ---------------- G4QChipolino ----------------
30 // by Mikhail Kossov, Sept 1999.
31 // class for Quasmon initiated Chipolinos generated by CHIPS Model
32 // --------------------------------------------------------------------
33 // Short description: In the CHIPS model not only hadrons are considered,
34 // but the di-hadrons, which can not be convereged to the quark content
35 // of only one hadron (e.g. pi+pi+, K+p, Delta++p etc). This kind of
36 // hadronic states, which can be easily decayed in two hadrons, is called
37 // Chipolino-particle in the model.
38 // ----------------------------------------------------------------------
39 
40 //#define debug
41 //#define pdebug
42 
43 #include "G4QChipolino.hh"
44 #include <cstdlib>
45 using namespace std;
46 
48 {
49  // @@ Does not work as static const ??
50  G4QPDGCode Pi0(111);
51  G4double mPi0 = Pi0.GetMass();
52  G4QContent Pi0QC = Pi0.GetQuarkContent();
53  G4int ban =QCont.GetBaryonNumber();
54  G4int tban=abs(3*ban);
55  G4int tot=QCont.GetTot(); // Initial total number of quarks in QC
56  G4int tod=tot%2; // tot is even - meson or dibaryon-nucleus
57  if ( (!tod && (tot < 4 || (ban && tot < tban) ) ) || (tod && tot < tban+2) )
58  QCont.IncQAQ(1,0.); // Add quark-pair
59 
60  G4QContent rQC=QCont; // Copy for possible reduction ("annihilation" of q-aq pairs)
61  tot=rQC.GetTot(); // New total number of quarks in QC (temporary)
62  if (tot%2)rQC.DecQAQ(-tban-2); // Reduce pairs, keep only 5 quarks (baryon case)
63  else if(ban)rQC.DecQAQ(-tban); // Reduce pairs, keep only 6 quarks (dibaryon case)
64  else rQC.DecQAQ(-4); // Reduce pairs, keep only 4 quarks (meson case)
65  tot=rQC.GetTot(); // Final total number of quarks (updated)
66 #ifdef debug
67  cout<<"G4QChipolino is called with QC="<<QCont<<",rQC="<<rQC<<",tot="<<tot<<G4endl;
68 #endif
69  minM=1000000.; // Prototype of minimal mass (@@ just a big number)
70  theQPDG1 = Pi0;
71  theQPDG2 = Pi0;
72  theQCont1 = Pi0QC;
73  if (!tot) // Should not be here, just in case (strange input)
74  {
75  G4cerr<<"***G4QChipolino: shouldn't be here 1 QC="<<rQC<<G4endl;
76  }
77  else if (tot==2 || tot==3) // Should not be here (basic octet/singlet states)
78  {
79  G4cerr<<"***G4QChipolino: shouldn't be here 2 QC="<<rQC<<G4endl;
80  theQCont1= rQC;
81  theQPDG1.InitByQCont(rQC);
82  theQCont = rQC+Pi0QC;
83  }
84  else if (tot==4) // Two possible combinations for the meson
85  {
86  G4QContent bQC=rQC.IndQ();
87 #ifdef debug
88  G4cout<<"G4QChipolino: tot=4,rQC="<<rQC<<",bQC="<<bQC<<G4endl;
89 #endif
90  for(int j=0; j<2; j++)
91  {
92  G4QContent aQC=rQC.IndAQ(j);
93  G4QContent cQC=bQC+aQC;
94  G4QPDGCode cQPDG(cQC);
95  G4double M1=cQPDG.GetMass();
96  if(cQPDG.GetPDGCode()==221) M1=mPi0;
97  G4QContent oQC=rQC-cQC;
98 #ifdef debug
99  cout<<"G4QChipolino: aQC="<<aQC<<", cQC="<<cQC<<", oQC="<<oQC<<G4endl;
100 #endif
101  G4QPDGCode oQPDG(oQC);
102  G4double M2=oQPDG.GetMass();
103  if(oQPDG.GetPDGCode()==221) M2=mPi0;
104  G4double m_value=M1+M2;
105 #ifdef debug
106  cout<<"G4QChipolino: c="<<cQPDG<<",cM="<<M1<<",o="<<oQPDG<<",oM="<<M2
107  <<",cM+0M="<<m_value<<", curMinM="<<minM<<G4endl;
108 #endif
109  if(m_value<minM)
110  {
111  minM=m_value;
112  theQPDG1 = cQPDG;
113  theQCont1 = cQC;
114  theQPDG2 = oQPDG;
115  }
116  }
117  }
118  else if (tot==5) // Four possible combinations for the baryon
119  {
120  G4int nQ=rQC.GetQ();
121  G4int nA=rQC.GetAQ();
122  G4bool fl=nA>nQ; // Flag of antibaryon case
123 #ifdef pdebug
124  cout<<"G4QChipolino: Baryon case nQ="<<nQ<<",nA="<<nA<<",QC="<<rQC
125  <<",fl="<<fl<<G4endl;
126 #endif
127  G4QContent bQC;
128  if (fl) bQC=rQC.IndQ(); // Antibaryon case
129  else bQC=rQC.IndAQ(); // Baryon case - QC of antiquark
130  for (int i=0; i<4; i++)
131  {
132  G4QContent cQC;
133  if (fl) cQC=bQC+rQC.IndAQ(i);
134  else cQC=bQC+rQC.IndQ(i);// Make mesonout of anti-quark
135  G4QPDGCode cQPDG(cQC); // Make QPDG particle
136  G4double M1=cQPDG.GetMass(); // Get meson mass
137  if(cQPDG.GetPDGCode()==221) M1=mPi0; // Make pi0 out of eta
138  G4QContent oQC=rQC-cQC; // Make residual baryon
139  G4QPDGCode oQPDG(oQC); // Make QPDG of residual baryon
140  G4double M2=oQPDG.GetMass(); // Get baryon mass
141  if(oQPDG.GetPDGCode()==221) M2=mPi0; // @@ Never !!
142  G4double m_value=M1+M2;
143  if(m_value<minM)
144  {
145  minM=m_value;
146  theQPDG1 = cQPDG;
147  theQCont1 = cQC;
148  theQPDG2 = oQPDG;
149  }
150  }
151 #ifdef pdebug
152  cout<<"G4QChipolino: Baryon case minM="<<minM<<", M="<<theQCont1<<theQPDG1
153  <<", B="<<theQPDG2<<G4endl;
154 #endif
155  }
156  else if (tot==6) // Four possible combinations for the di-baryon
157  {
158  if(ban)
159  {
160  G4int nQ=rQC.GetQ();
161  G4int nA=rQC.GetAQ();
162  G4bool fl=nA>nQ; // Flag of anti-dibaryon case
163 #ifdef debug
164  cout<<"G4QChipolino: Di-Bar. case nQ="<<nQ<<",nA="<<nA<<",QC="<<rQC<<",fl="<<fl<<G4endl;
165 #endif
166  for (int i=0; i<4; i++)
167  {
168  G4QContent aQC;
169  if (fl) aQC=rQC.IndAQ(i);
170  else aQC=rQC.IndQ(i);
171  for (int j=i+1; j<5; j++)
172  {
173  G4QContent bQC;
174  if (fl) bQC=aQC+rQC.IndAQ(j);
175  else bQC=aQC+rQC.IndQ(j);
176  for (int k=j+1; k<6; k++)
177  {
178  G4QContent cQC;
179  if (fl) cQC=bQC+rQC.IndAQ(k);
180  else cQC=bQC+rQC.IndQ(k);
181  G4QPDGCode cQPDG(cQC);
182  G4double M1=cQPDG.GetMass();
183  if(cQPDG.GetPDGCode()==221) M1=mPi0;
184  G4QContent oQC=rQC-cQC;
185  G4QPDGCode oQPDG=(oQC);
186  G4double M2=oQPDG.GetMass();
187  if(oQPDG.GetPDGCode()==221) M2=mPi0;
188  G4double m_value=M1+M2;
189  if(m_value<minM)
190  {
191  minM=m_value;
192  theQPDG1 = cQPDG;
193  theQCont1 = cQC;
194  theQPDG2 = oQPDG;
195  }
196  }
197  }
198  }
199  }
200  else // Baryon-AntiBaryon
201  {
202  theQCont1 = rQC.IndQ(0)+rQC.IndQ(1)+rQC.IndQ(2);
203  theQPDG1.InitByQCont(theQCont1);
204  theQPDG2.InitByQCont(rQC.IndAQ(0)+rQC.IndAQ(1)+rQC.IndAQ(2));
205  }
206  }
207  else if(((rQC.GetU() )>(rQC.GetS() -4) && (rQC.GetD() )>(rQC.GetS() -4)) ||
208  ((rQC.GetAU())>(rQC.GetAS()-4) && (rQC.GetAD())>(rQC.GetAS()-4)) )
209  {
210  G4int kD=rQC.GetD();
211  G4int kU=rQC.GetU();
212  G4int kS=rQC.GetS();
213  G4int mD=rQC.GetAD();
214  G4int mU=rQC.GetAU();
215  G4int mS=rQC.GetAS();
216  G4int nQ=rQC.GetQ();
217  G4int nA=rQC.GetAQ();
218  G4bool fl=nA>nQ; // Flag of anti-fragment case
219 #ifdef debug
220  G4cout<<"G4QChipolino: NucFragment case nQ="<<nQ<<",nAQ="<<nA<<", QC="<<rQC<<",fl="<<fl
221  <<G4endl;
222 #endif
223  if( (fl && kS>1) || (!fl && mS>1))
224  {
225 #ifdef debug
226  G4cerr<<"***G4QChipolino: ***Overfowed by strange quarks*** rQC="<<rQC<<G4endl;
227  //throw G4QException("G4QChipolino: NuclearFragment is overflowed by strangeQuarks");
228 #endif
229  }
230  else if(fl) // ===> Anti-fragment
231  {
232  //G4cerr<<"***G4QChipolino: ***Anti-nuclear fragments*** rQC="<<rQC<<G4endl;
233  //throw G4QException("G4QChipolino: Antinuclear fragments are not yet supported");
234  if(!mS) // No strange quarks
235  {
236  G4int nI=mU-mD; // Isotopic shift
237  G4int nN=(mU+mD-nI*3)/6;
238  if(!kS) // No kaons
239  {
240  if((nI>=0&&nN>=0)||(nI<0&&nN>=-nI)) // Delta isn't necessary
241  {
242  if(nI>0) // Excess of antiprotons
243  {
244  theQPDG1 = G4QPDGCode(-(90000000+1000*(nN+nI-1)+nN));// A Fragment-AProton
245  theQPDG2 = G4QPDGCode(-2212); // An Anti-Proton
246  }
247  else // Excess of a-neutrons
248  {
249  theQPDG1 = G4QPDGCode(-(90000000+1000*(nN+nI)+nN-1));// A Fragment-ANeutron
250  theQPDG2 = G4QPDGCode(-2112); // An Anti-Neutron
251  }
252  }
253  else if((nI>=0&&nN>-2)||(nI<0&&nN>-nI-2)) // Delta can be a part
254  {
255  if(nI>0) // Excess of au-quarks
256  {
257  theQPDG1=G4QPDGCode(-(90000000+1000*(nN+nI-2)+nN+1));// A Fragment-AProton
258  theQPDG2=G4QPDGCode(-2224); // An Anti-Delta++
259  }
260  else // Excess of ad-quarks
261  {
262  theQPDG1=G4QPDGCode(-(90000000+1000*(nN+nI+1)+nN-2));// A Fragment-ANeutron
263  theQPDG2=G4QPDGCode(-1114); // An Anti-Delta-
264  }
265  }
266  else
267  {
268  G4cerr<<"***G4QChipolino:**A**IsotopicAsymmetry (without S),rQC="<<rQC<<G4endl;
269  //throw G4QException("G4QChipolino: IsotopicAsymmety of AntiMultyBar Quasmon");
270  }
271  }
272  else if(kS<2) // NucFrag+K is possible
273  {
274  nN =(mU+mD-4-nI*3)/6;
275  if(nI>0) // Excess of au-quarks
276  {
277  nN+=1;
278  theQPDG1 = G4QPDGCode(-(90000000+1000*(nN+nI-1)+nN)); // An Anti-Fragment
279  theQPDG2 = G4QPDGCode(-321); // A K- meson
280  }
281  else
282  {
283  theQPDG1 = G4QPDGCode(-(90000000+1000*(nN+nI+1)+nN)); // An AntiFragment
284  theQPDG2 = G4QPDGCode(-311); // An Anti-K0 meson
285  }
286  }
287  else
288  {
289  G4cerr<<"***G4QChipolino: ***Too many kaons are needed*** rQC="<<rQC<<G4endl;
290  //throw G4QException("G4QChipolino: Too Many Kaons are needed for AntiNucFragm");
291  }
292  }
293  else // Fragment with strangeness
294  {
295  if(mS<=mU&&mS<=mD) // Fragment consisting of Neutrons, Protons & Lambrdas only
296  {
297  G4int nI=mU-mD; // Isotopic shift
298  G4int nN=(mU+mD-mS-mS-nI*3)/6;
299  if((nI>=0&&nN>=0)||(nI<0&&nN>=-nI)) // Delta isn't necessary
300  {
301  if(nI>0) // Excess of protons
302  {
303  theQPDG1 = G4QPDGCode(-(90000000+1000*(kS*1000+nN+nI-1)+nN));// Fragm-AProton
304  theQPDG2 = G4QPDGCode(-2212); // An Anti-Proton
305  }
306  else // Excess of neutrons
307  {
308  theQPDG1 = G4QPDGCode(-(90000000+1000*(kS*1000+nN+nI)+nN-1));//Fragm-ANeutron
309  theQPDG2 = G4QPDGCode(-2112); // An Anti-Neutron
310  }
311  }
312  else if((nI>=0&&nN>-2)||(nI<0&&nN>-nI-2)) // Delta can be a part
313  {
314  if(nI>0) // Excess of au-quarks
315  {
316  theQPDG1=G4QPDGCode(-(90000000+1000*(kS*1000+nN+nI-2)+nN+1));//Fragm-AProton
317  theQPDG2=G4QPDGCode(-2224); // An Anti-Delta++
318  }
319  else // Excess of ad-quarks
320  {
321  theQPDG1=G4QPDGCode(-(90000000+1000*(kS*1000+nN+nI+1)+nN-2));//Fragm-ANeutron
322  theQPDG2=G4QPDGCode(-1114); // An Anti-Delta-
323  }
324  }
325  else
326  {
327  G4cerr<<"***G4QChipolino:**A**IsotopicAssimetry (with S)*** rQC="<<rQC<<G4endl;
328  //throw G4QException("G4QChipolino: Isotopics of Strange AntiMultyBarQuasmon");
329  }
330  }
331  else // Excess of s-quarks
332  {
333  G4int lam=mU; // A#of Anti-Lambdas
334  if (lam>mD) lam=mD;
335  G4int lD=mD-lam; // Residual ad-quarks
336  G4int lU=mU-lam; // Residual au-quarks
337  G4int lS=mS-lam; // Residual as-quarks
338  if(lD+lU+lS!=3||lD<0||lU<0||lS<0)
339  {
340  G4cerr<<"***G4QChipolino:*AntiFragment* rQC="<<rQC<<",s="<<lS<<",u="<<lU<<",d"
341  <<lD<<G4endl;
342  //throw G4QException("G4QChipolino: Exotic superstrange AntiMultyBaryon");
343  }
344  if ( !lD && lU==2) theQPDG2=G4QPDGCode(-3222); // Anti-Sigma+
345  else if( !lU && lD==2) theQPDG2=G4QPDGCode(-3112); // Anti-Sigma-
346  else if( !lD && lU==1) theQPDG2=G4QPDGCode(-3322); // Anti-Ksi0
347  else if( !lU && lD==1) theQPDG2=G4QPDGCode(-3312); // Anti-Ksi-
348  else theQPDG2=G4QPDGCode(-3334); // Anti-Omega-
349  theQPDG1=G4QPDGCode(-(90+lam)*1000000); // Anti Strange Matter
350  }
351  theQCont1 = rQC-theQPDG2.GetQuarkContent(); // QCont of Fragment-H
352  theQCont = rQC; // QCont of Chipolino
353  }
354  }
355  else // ===> Nuclear Fragment
356  {
357  if(!kS) // No strange quarks
358  {
359  G4int nI=kU-kD; // Isotopic shift
360  G4int nN=(kU+kD-nI*3)/6;
361  if(!mS) // No kaons
362  {
363  if((nI>=0&&nN>=0)||(nI<0&&nN>=-nI)) // Delta isn't necessary
364  {
365  if(nI>0) // Excess of protons
366  {
367  theQPDG1 = G4QPDGCode(90000000+1000*(nN+nI-1)+nN); // A Fragment-Proton
368  theQPDG2 = G4QPDGCode(2212); // A Proton
369  }
370  else // Excess of neutrons
371  {
372  theQPDG1 = G4QPDGCode(90000000+1000*(nN+nI)+nN-1); // A Fragment-Neutron
373  theQPDG2 = G4QPDGCode(2112); // A Neutron
374  }
375  }
376  else if((nI>=0&&nN>-2)||(nI<0&&nN>-nI-2)) // Delta can be a part
377  {
378  if(nI>0) // Excess of u-quarks
379  {
380  theQPDG1=G4QPDGCode(90000000+1000*(nN+nI-2)+nN+1); // A Fragment-Proton
381  theQPDG2=G4QPDGCode(2224); // A Delta++
382  }
383  else // Excess of d-quarks
384  {
385  theQPDG1=G4QPDGCode(90000000+1000*(nN+nI+1)+nN-2); // A Fragment-Neutron
386  theQPDG2=G4QPDGCode(1114); // A Delta-
387  }
388  }
389  else
390  {
391  G4cerr<<"***G4QChipolino:***Isotopic assimetry (without S), rQC="<<rQC<<G4endl;
392  //throw G4QException("G4QChipolino:ExoticIsotopicAssimety of MultyBarQuasmon");
393  }
394  }
395  else if(mS<2) // NucFrag+K is possible
396  {
397  nN =(kU+kD-4-nI*3)/6;
398  if(nI>0) // Excess of u-quarks
399  {
400  nN+=1;
401  theQPDG1 = G4QPDGCode(90000000+1000*(nN+nI-1)+nN); // A Fragment
402  theQPDG2 = G4QPDGCode(321); // A K+ meson
403  }
404  else
405  {
406  theQPDG1 = G4QPDGCode(90000000+1000*(nN+nI+1)+nN); // A Fragment
407  theQPDG2 = G4QPDGCode(311); // A K0 meson
408  }
409  }
410  else
411  {
412  G4cerr<<"***G4QChipolino: ***Too many kaons are needed*** rQC="<<rQC<<G4endl;
413  //throw G4QException("G4QChipolino: More than one Kaon is needed for NuclFragm");
414  }
415  }
416  else // Fragment with strangeness
417  {
418  if(kS<=kU&&kS<=kD) // Fragment consisting of Neutrons, Protons & Lambrdas only
419  {
420  G4int nI=kU-kD; // Isotopic shift
421  G4int nN=(kU+kD-kS-kS-nI*3)/6;
422  if((nI>=0&&nN>=0)||(nI<0&&nN>=-nI)) // Delta isn't necessary
423  {
424  if(nI>0) // Excess of protons
425  {
426  theQPDG1 = G4QPDGCode(90000000+1000*(kS*1000+nN+nI-1)+nN);// Fragment-Proton
427  theQPDG2 = G4QPDGCode(2212); // A Proton
428  }
429  else // Excess of neutrons
430  {
431  theQPDG1 = G4QPDGCode(90000000+1000*(kS*1000+nN+nI)+nN-1);// Fragment-Neutron
432  theQPDG2 = G4QPDGCode(2112); // A Neutron
433  }
434  }
435  else if((nI>=0&&nN>-2)||(nI<0&&nN>-nI-2)) // Delta can be a part
436  {
437  if(nI>0) // Excess of u-quarks
438  {
439  theQPDG1=G4QPDGCode(90000000+1000*(kS*1000+nN+nI-2)+nN+1);// Fragment-Proton
440  theQPDG2=G4QPDGCode(2224); // A Delta++
441  }
442  else // Excess of d-quarks
443  {
444  theQPDG1=G4QPDGCode(90000000+1000*(kS*1000+nN+nI+1)+nN-2);// Fragment-Neutron
445  theQPDG2=G4QPDGCode(1114); // A Delta-
446  }
447  }
448  else
449  {
450  G4cerr<<"***G4QChipolino: ***Isotopic assimetry (with S)*** rQC="<<rQC<<G4endl;
451  //throw G4QException("G4QChipolino:IsotopicAssimety of StrangeMultyBar Quasm");
452  }
453  }
454  else // Excess of s-quarks
455  {
456  G4int lam=kU; // A#of Lambda
457  if (lam>kD) lam=kD;
458  G4int lD=kD-lam; // Residual d-quarks
459  G4int lU=kU-lam; // Residual u-quarks
460  G4int lS=kS-lam; // Residual s-quarks
461  if(lD+lU+lS!=3||lD<0||lU<0||lS<0)
462  {
463  G4cerr<<"***G4QChipolino:*Fragment*rQC="<<rQC<<",s="<<lS<<",u="<<lU<<",d"
464  <<lD<<G4endl;
465  //throw G4QException("G4QChipolino: Exotic superstrange Multy Baryon");
466  }
467  if ( !lD && lU==2) theQPDG2=G4QPDGCode(3222); // Sigma+
468  else if( !lU && lD==2) theQPDG2=G4QPDGCode(3112); // Sigma-
469  else if( !lD && lU==1) theQPDG2=G4QPDGCode(3322); // Ksi0
470  else if( !lU && lD==1) theQPDG2=G4QPDGCode(3312); // Ksi-
471  else theQPDG2=G4QPDGCode(3334); // Omega-
472  theQPDG1=G4QPDGCode((90+lam)*1000000); // Strange Matter
473  }
474  theQCont1 = rQC-theQPDG2.GetQuarkContent(); // QCont of Fragment-H
475  theQCont = rQC; // QCont of Chipolino
476  }
477  }
478  }
479  else
480  {
481  G4cerr<<"***G4QChipolino: ***Exotics*** rQC="<<rQC<<G4endl;
482  //throw G4QException("G4QChipolino: can't be constructed for exotic baryon or meson");
483  }
484 }
485 
487 {
488  theQPDG1 = right.theQPDG1;
489  theQPDG2 = right.theQPDG2;
490  theQCont = right.theQCont;
491  theQCont1 = right.theQCont1;
492  minM = right.minM;
493 }
494 
496 {
497  theQPDG1 = right->theQPDG1;
498  theQPDG2 = right->theQPDG2;
499  theQCont = right->theQCont;
500  theQCont1 = right->theQCont1;
501  minM = right->minM;
502 }
503 
505 {
506  if(this != &right) // Beware of self assignment
507  {
508  theQPDG1 = right.theQPDG1;
509  theQPDG2 = right.theQPDG2;
510  theQCont = right.theQCont;
511  theQCont1 = right.theQCont1;
512  minM = right.minM;
513  }
514  return *this;
515 }
516 
518 
519 // Standard output for G4QChipolino
520 ostream& operator<<(ostream& lhs, G4QChipolino& rhs)
521 {
522  lhs<<"{1="<<rhs.GetQPDG1()<<",2="<<rhs.GetQPDG2()<< "}";
523  return lhs;
524 }
525 
526 
527 
528 
529 
530