Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4QString.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 // ------------------------------------------------------------
30 // GEANT 4 class implementation file
31 //
32 // ---------------- G4QString ----------------
33 // by Mikhail Kossov Oct, 2006
34 // class for excited string used by Parton String Models
35 // For comparison mirror member functions are taken from G4 classes:
36 // G4FragmentingString
37 // G4ExcitedStringDecay
38 // ---------------------------------------------------------------------
39 // Short description: If partons from the G4QPartonPair are close in
40 // rapidity, they create Quasmons, but if they are far in the rapidity
41 // space, they can not interact directly. Say the bottom parton (quark)
42 // has rapidity 0, and the top parton (antiquark) has rapidity 8, then
43 // the top quark splits in two by radiating gluon, and each part has
44 // rapidity 4, then the gluon splits in quark-antiquark pair (rapidity
45 // 2 each), and then the quark gadiates anothe gluon and reachs rapidity
46 // 1. Now it can interact with the bottom antiquark, creating a Quasmon
47 // or a hadron. The intermediate partons is the string ladder.
48 // ---------------------------------------------------------------------
49 
50 //#define debug
51 //#define pdebug
52 //#define edebug
53 
54 #include <algorithm>
55 
56 #include "G4QString.hh"
57 #include "G4PhysicalConstants.hh"
58 #include "G4SystemOfUnits.hh"
59 
60 // Static parameters definition
61 G4double G4QString::MassCut=350.*MeV; // minimum mass cut for the string
62 G4double G4QString::SigmaQT=0.5*GeV; // quarkTransverseMomentum distribution parameter
63 G4double G4QString::DiquarkSuppress=0.1; // is Diquark suppression parameter
64 G4double G4QString::DiquarkBreakProb=0.1; // is Diquark breaking probability
65 G4double G4QString::SmoothParam=0.9; // QGS model parameter
66 G4double G4QString::StrangeSuppress=0.435;// Strangeness suppression (u:d:s=1:1:0.3 ?M.K.)
67 G4double G4QString::widthOfPtSquare=-0.72*GeV*GeV; // pt -width2 forStringExcitation
68 
69 G4QString::G4QString() : theDirection(0), thePosition(G4ThreeVector(0.,0.,0.)),
70  theStableParton(0), theDecayParton(0){}
71 
72 G4QString::G4QString(G4QParton* Color, G4QParton* AntiColor, G4int Direction)
73  : SideOfDecay(0)
74 {
75 #ifdef debug
76  G4cout<<"G4QString::PPD-Constructor: Direction="<<Direction<<G4endl;
77 #endif
78  ExciteString(Color, AntiColor, Direction);
79 #ifdef debug
80  G4cout<<"G4QString::PPD-Constructor: ------>> String is excited"<<G4endl;
81 #endif
82 }
83 
84 G4QString::G4QString(G4QPartonPair* CAC): SideOfDecay(0)
85 {
86 #ifdef debug
87  G4cout<<"G4QString::PartonPair-Constructor: Is CALLED"<<G4endl;
88 #endif
89  ExciteString(CAC->GetParton1(), CAC->GetParton2(), CAC->GetDirection());
90 #ifdef debug
91  G4cout<<"G4QString::PartonPair-Constructor: ------>> String is excited"<<G4endl;
92 #endif
93 }
94 
95 G4QString::G4QString(G4QParton* QCol,G4QParton* Gluon,G4QParton* QAntiCol,G4int Direction)
96  : theDirection(Direction), thePosition(QCol->GetPosition()), SideOfDecay(0)
97 {
98  thePartons.push_back(QCol);
99  G4LorentzVector sum=QCol->Get4Momentum();
100  thePartons.push_back(Gluon);
101  sum+=Gluon->Get4Momentum();
102  thePartons.push_back(QAntiCol);
103  sum+=QAntiCol->Get4Momentum();
104  Pplus =sum.e() + sum.pz();
105  Pminus=sum.e() - sum.pz();
106  decaying=None;
107 }
108 
109 G4QString::G4QString(const G4QString &right) : theDirection(right.GetDirection()),
110  thePosition(right.GetPosition()), SideOfDecay(0)
111 {
112  //LeftParton=right.LeftParton;
113  //RightParton=right.RightParton;
114  Ptleft=right.Ptleft;
115  Ptright=right.Ptright;
116  Pplus=right.Pplus;
117  Pminus=right.Pminus;
118  decaying=right.decaying;
119 }
120 
121 G4QString::~G4QString()
122 {if(thePartons.size()) std::for_each(thePartons.begin(),thePartons.end(),DeleteQParton());}
123 
124 G4LorentzVector G4QString::Get4Momentum() const
125 {
126  G4LorentzVector momentum(0.,0.,0.,0.);
127  for(unsigned i=0; i<thePartons.size(); i++) momentum += thePartons[i]->Get4Momentum();
128  return momentum;
129 }
130 
131 void G4QString::LorentzRotate(const G4LorentzRotation & rotation)
132 {
133  for(unsigned i=0; i<thePartons.size(); i++)
134  thePartons[i]->Set4Momentum(rotation*thePartons[i]->Get4Momentum());
135 }
136 
137 //void G4QString::InsertParton(G4QParton* aParton, const G4QParton* addafter)
138 //{
139 // G4QPartonVector::iterator insert_index; // Begin by default (? M.K.)
140 // if(addafter != NULL)
141 // {
142 // insert_index=std::find(thePartons.begin(), thePartons.end(), addafter);
143 // if (insert_index == thePartons.end()) // No object addafter in thePartons
144 // {
145 // G4cerr<<"***G4QString::InsertParton: Addressed Parton is not found"<<G4endl;
146 // G4Exception("G4QString::InsertParton:","72",FatalException,"NoAddressParton");
147 // }
148 // }
149 // thePartons.insert(insert_index+1, aParton);
150 //}
151 
152 void G4QString::Boost(G4ThreeVector& Velocity)
153 {
154  for(unsigned cParton=0; cParton<thePartons.size(); cParton++)
155  {
156  G4LorentzVector Mom = thePartons[cParton]->Get4Momentum();
157  Mom.boost(Velocity);
158  thePartons[cParton]->Set4Momentum(Mom);
159  }
160 }
161 
162 //G4QParton* G4QString::GetColorParton(void) const
163 //{
164 // G4QParton* start = *(thePartons.begin());
165 // G4QParton* end = *(thePartons.end()-1);
166 // G4int Encoding = start->GetPDGCode();
167 // if (Encoding<-1000 || (Encoding < 1000 && Encoding > 0)) return start;
168 // return end;
169 //}
170 
171 //G4QParton* G4QString::GetAntiColorParton(void) const
172 //{
173 // G4QParton* start = *(thePartons.begin());
174 // G4QParton* end = *(thePartons.end()-1);
175 // G4int Encoding = start->GetPDGCode();
176 // if (Encoding < -1000 || (Encoding < 1000 && Encoding > 0)) return end;
177 // return start;
178 //}
179 
180 // Fill parameters
181 void G4QString::SetParameters(G4double mCut, G4double sigQT, G4double DQSup, G4double DQBU,
182  G4double smPar, G4double SSup, G4double SigPt)
183 {
184  MassCut = mCut; // minimum mass cut for the string
185  SigmaQT = sigQT; // quark transverse momentum distribution parameter
186  DiquarkSuppress = DQSup; // is Diquark suppression parameter
187  DiquarkBreakProb= DQBU; // is Diquark breaking probability
188  SmoothParam = smPar; // QGS model parameter
189  StrangeSuppress = SSup; // Strangeness suppression parameter
190  widthOfPtSquare = -2*SigPt*SigPt; // width^2 of pt for string excitation
191 }
192 
193 // Pt distribution @@ one can use 1/(1+A*Pt^2)^B
194 G4ThreeVector G4QString::GaussianPt(G4double widthSquare, G4double maxPtSquare) const
195 {
196  G4double pt2; do{pt2=widthSquare*std::log(G4UniformRand());} while (pt2>maxPtSquare);
197  pt2=std::sqrt(pt2);
198  G4double phi=G4UniformRand()*twopi;
199  return G4ThreeVector(pt2*std::cos(phi),pt2*std::sin(phi),0.);
200 } // End of GaussianPt
201 
202 // Diffractively excite the string
203 //void G4QString::DiffString(G4QHadron* hadron, G4bool isProjectile)
204 //{
205 // hadron->SplitUp();
206 // G4QParton* start = hadron->GetNextParton();
207 // if( start==NULL)
208 // {
209 // G4cerr<<"***G4QString::DiffString: No start parton found, nothing is done"<<G4endl;
210 // return;
211 // }
212 // G4QParton* end = hadron->GetNextParton();
213 // if( end==NULL)
214 // {
215 // G4cerr<<"***G4QString::DiffString: No end parton found, nothing is done"<<G4endl;
216 // return;
217 // }
218 // if(isProjectile) ExciteString(end, start, 1); // 1 = Projectile
219 // else ExciteString(start, end,-1); // -1 = Target
220 // SetPosition(hadron->GetPosition());
221 // // momenta of string ends
222 // G4double ptSquared= hadron->Get4Momentum().perp2();
223 // G4double hmins=hadron->Get4Momentum().minus();
224 // G4double hplus=hadron->Get4Momentum().plus();
225 // G4double transMassSquared=hplus*hmins;
226 // G4double maxMomentum = std::sqrt(transMassSquared) - std::sqrt(ptSquared);
227 // G4double maxAvailMomentumSquared = maxMomentum * maxMomentum;
228 // G4ThreeVector pt=GaussianPt(widthOfPtSquare,maxAvailMomentumSquared);
229 // G4LorentzVector Pstart(G4LorentzVector(pt,0.));
230 // G4LorentzVector Pend(hadron->Get4Momentum().px(), hadron->Get4Momentum().py(), 0.);
231 // Pend-=Pstart;
232 // G4double tm1=hmins+(Pend.perp2()-Pstart.perp2())/hplus;
233 // G4double tm2=std::sqrt( std::max(0., tm1*tm1-4*Pend.perp2()*hmins/hplus ) );
234 // G4int Sign= isProjectile ? TARGET : PROJECTILE;
235 // G4double endMinus = 0.5 * (tm1 + Sign*tm2);
236 // G4double startMinus= hmins - endMinus;
237 // G4double startPlus = Pstart.perp2() / startMinus;
238 // G4double endPlus = hplus - startPlus;
239 // Pstart.setPz(0.5*(startPlus - startMinus));
240 // Pstart.setE (0.5*(startPlus + startMinus));
241 // Pend.setPz (0.5*(endPlus - endMinus));
242 // Pend.setE (0.5*(endPlus + endMinus));
243 // start->Set4Momentum(Pstart);
244 // end->Set4Momentum(Pend);
245 //#ifdef debug
246 // G4cout<<"G4QString::DiffString: StartQ="<<start->GetPDGCode()<<start->Get4Momentum()<<"("
247 // <<start->Get4Momentum().mag()<<"), EndQ="<<end->GetPDGCode()<<end ->Get4Momentum()
248 // <<"("<<end->Get4Momentum().mag()<<"), sumOfEnds="<<Pstart+Pend<<", H4M(original)="
249 // <<hadron->Get4Momentum()<<G4endl;
250 //#endif
251 //} // End of DiffString (The string is excited diffractively)
252 
253 // Excite the string by two partons
254 void G4QString::ExciteString(G4QParton* Color, G4QParton* AntiColor, G4int Direction)
255 {
256 #ifdef debug
257  G4cout<<"G4QString::ExciteString: **Called**, direction="<<Direction<<G4endl;
258 #endif
259  if(thePartons.size()) std::for_each(thePartons.begin(),thePartons.end(),DeleteQParton());
260  thePartons.clear();
261  theDirection = Direction;
262  thePosition = Color->GetPosition();
263 #ifdef debug
264  G4cout<<"G4QString::ExciteString: ColourPosition = "<<thePosition<<", beg="<<Color->GetPDGCode()
265  <<",end="<<AntiColor->GetPDGCode()<<G4endl;
266 #endif
267  thePartons.push_back(Color);
268  G4LorentzVector sum=Color->Get4Momentum();
269  thePartons.push_back(AntiColor); // @@ Complain of Valgrind
270  sum+=AntiColor->Get4Momentum();
271  Pplus =sum.e() + sum.pz();
272  Pminus=sum.e() - sum.pz();
273  decaying=None;
274 #ifdef debug
275  G4cout<<"G4QString::ExciteString: ***Done***, beg="<<(*thePartons.begin())->GetPDGCode()
276  <<",end="<<(*(thePartons.end()-1))->GetPDGCode()<<G4endl;
277 #endif
278 } // End of ExciteString
279 
280 // LUND Longitudinal fragmentation
281 G4double G4QString::GetLundLightConeZ(G4double zmin, G4double zmax, G4int , // @@ ? M.K.
282  G4QHadron* pHadron, G4double Px, G4double Py)
283 {
284  static const G4double alund = 0.7/GeV/GeV;
285  // If blund get restored, you MUST adapt the calculation of zOfMaxyf.
286  //static const G4double blund = 1;
287  G4double z, yf;
288  G4double Mt2 = Px*Px + Py*Py + pHadron->GetMass2();
289  G4double zOfMaxyf=alund*Mt2/(alund*Mt2+1.);
290  G4double maxYf=(1.-zOfMaxyf)/zOfMaxyf * std::exp(-alund*Mt2/zOfMaxyf);
291  do
292  {
293  z = zmin + G4UniformRand()*(zmax-zmin);
294  // yf = std::pow(1. - z, blund)/z*std::exp(-alund*Mt2/z);
295  yf = (1-z)/z * std::exp(-alund*Mt2/z);
296  } while (G4UniformRand()*maxYf>yf);
297  return z;
298 } // End of GetLundLightConeZ
299 
300 
301 // QGSM Longitudinal fragmentation
302 G4double G4QString::GetQGSMLightConeZ(G4double zmin, G4double zmax, G4int PartonEncoding,
303  G4QHadron* , G4double, G4double) // @@ M.K.
304 {
305  static const G4double arho=0.5;
306  static const G4double aphi=0.;
307  static const G4double an=-0.5;
308  static const G4double ala=-0.75;
309  static const G4double aksi=-1.;
310  static const G4double alft=0.5;
311  G4double z;
312  G4double theA(0), d2, yf;
313  G4int absCode = std::abs(PartonEncoding);
314  // @@ Crazy algorithm is used for simple power low...
315  if (absCode < 10) // Quarks (@@ 9 can be a gluon)
316  {
317  if(absCode == 1 || absCode == 2) theA = arho;
318  else if(absCode == 3) theA = aphi;
319  else
320  {
321  G4cerr<<"***G4QString::GetQGSMLightConeZ: CHIPS is SU(3), quakCode="<<absCode<<G4endl;
322  G4Exception("G4QString::GetQGSMLightConeZ:","72",FatalException,"WrongQuark");
323  }
324  do
325  {
326  z = zmin + G4UniformRand()*(zmax - zmin);
327  d2 = alft - theA;
328  yf = std::pow( 1.-z, d2);
329  }
330  while (G4UniformRand()>yf);
331  }
332  else // Di-quarks (@@ Crazy codes are not checked)
333  {
334  if (absCode==3101||absCode==3103||absCode==3201||absCode==3203) d2=alft-ala-ala+arho;
335  else if(absCode==1103||absCode==2101||absCode==2203||absCode==2103) d2=alft-an-an+arho;
336  else d2=alft-aksi-aksi+arho;
337  do
338  {
339  z = zmin + G4UniformRand()*(zmax - zmin);
340  yf = std::pow(1.-z, d2);
341  }
342  while (G4UniformRand()>yf);
343  }
344  return z;
345 } // End of GetQGSMLightConeZ
346 
347 // Diffractively excite the string (QL=true - QGS Light Cone, =false - Lund Light Cone)
348 G4QHadronVector* G4QString::FragmentString(G4bool QL)
349 {
350  // Can no longer modify Parameters for Fragmentation.
351 #ifdef edebug
352  G4LorentzVector string4M=Get4Momentum(); // Just for Energy-Momentum ConservCheck
353 #endif
354 #ifdef debug
355  G4cout<<"G4QString::FragmentString:-->Called,QL="<<QL<<", M="<<Get4Momentum().m()<<", L="
356  <<GetLeftParton()->Get4Momentum()<<",R="<<GetRightParton()->Get4Momentum()<<G4endl;
357 #endif
358  // check if string has enough mass to fragment. If not, convert to one or two hadrons
359  G4QHadronVector* LeftVector = LightFragmentationTest();
360  if(LeftVector)
361  {
362 #ifdef edebug
363  G4LorentzVector sL=string4M;
364  for(unsigned L = 0; L < LeftVector->size(); L++)
365  {
366  G4QHadron* LH = (*LeftVector)[L];
367  G4LorentzVector L4M=LH->Get4Momentum();
368  sL-=L4M;
369  G4cout<<"-EMC-G4QStr::FragStr:L#"<<L<<",PDG="<<LH->GetPDGCode()<<",4M="<<L4M<<G4endl;
370  }
371  G4cout<<"-EMC-G4QString::FragmentString:---LightFragmentation---> Res4M="<<sL<<G4endl;
372 #endif
373  return LeftVector; //@@ Just decay in 2 or 1 (?) hadron, if below theCut
374  }
375 #ifdef debug
376  G4cout<<"G4QString::FragmentString:OUTPUT is not yet defined, define Left/Right"<<G4endl;
377 #endif
378  LeftVector = new G4QHadronVector; // Valgrind complain to LeftVector
379  G4QHadronVector* RightVector = new G4QHadronVector;
380  // Remember 4-momenta of the string ends (@@ only for the two-parton string, no gluons)
381  G4LorentzVector left4M=GetLeftParton()->Get4Momentum(); // For recovery when failed
382  G4LorentzVector right4M=GetRightParton()->Get4Momentum();
383 #ifdef debug
384  G4cout<<"G4QString::FragmString: ***Remember*** L4M="<<left4M<<", R4M="<<right4M<<G4endl;
385 #endif
386  G4int leftPDG=GetLeftParton()->GetPDGCode();
387  G4int rightPDG=GetRightParton()->GetPDGCode();
388  // Transform string to CMS
389  G4LorentzVector momentum=Get4Momentum();
390  G4LorentzRotation toCms(-(momentum.boostVector()));
391  momentum= toCms * thePartons[0]->Get4Momentum();
392  toCms.rotateZ(-1*momentum.phi());
393  toCms.rotateY(-1*momentum.theta());
394  for(unsigned index=0; index<thePartons.size(); index++)
395  {
396  momentum = toCms * thePartons[index]->Get4Momentum(); // @@ reuse of the momentum
397  thePartons[index]->Set4Momentum(momentum);
398  }
399  // Copy the string for independent attempts
400  G4QParton* LeftParton = new G4QParton(GetLeftParton());
401  G4QParton* RightParton= new G4QParton(GetRightParton());
402  G4QString* theStringInCMS = new G4QString(LeftParton,RightParton,GetDirection());
403 #ifdef debug
404  G4cout<<"G4QString::FragmentString: Copy with nP="<<theStringInCMS->thePartons.size()
405  <<", beg="<<(*(theStringInCMS->thePartons.begin()))->GetPDGCode()
406  <<", end="<<(*(theStringInCMS->thePartons.end()-1))->GetPDGCode()<<G4endl;
407 #endif
408  G4bool success=false;
409  G4bool inner_sucess=true;
410  G4int attempt=0;
411  G4int StringLoopInterrupt=27; // String fragmentation LOOP limit
412 #ifdef debug
413  G4cout<<"G4QString::FragmentString: BeforeWhileLOOP, max = "<<StringLoopInterrupt
414  <<", nP="<<thePartons.size()<<", beg="<<(*thePartons.begin())->GetPDGCode()
415  <<",end="<<(*(thePartons.end()-1))->GetPDGCode()<<G4endl;
416 #endif
417 #ifdef edebug
418  G4LorentzVector cmS4M=theStringInCMS->Get4Momentum();
419  G4cout<<"-EMC-G4QString::FragmString: c4M="<<cmS4M<<",Max="<<StringLoopInterrupt<<G4endl;
420 #endif
421  while (!success && attempt++ < StringLoopInterrupt) // Try fragment String till success
422  {
423  // Recover the CMS String
424  LeftParton = new G4QParton(theStringInCMS->GetLeftParton());
425  RightParton= new G4QParton(theStringInCMS->GetRightParton());
426  ExciteString(LeftParton, RightParton, theStringInCMS->GetDirection());
427 #ifdef edebug
428  G4LorentzVector cm4M=cmS4M; // Copy the full momentum for the reduction and check
429  G4cout<<"-EMC-.G4QString::FragmentString: CHEK "<<cm4M<<" ?= "<<Get4Momentum()<<G4endl;
430 #endif
431 #ifdef debug
432  G4cout<<"G4QString::FragmentString:===>LOOP, attempt = "<<attempt<<", nP="
433  <<thePartons.size()<<", beg="<<(*thePartons.begin())->GetPDGCode()
434  <<",end="<<(*(thePartons.end()-1))->GetPDGCode()<<G4endl;
435 #endif
436  // Now clean up all hadrons in the Left and Right vectors for the new attempt
437  if(LeftVector->size())
438  {
439  std::for_each(LeftVector->begin(), LeftVector->end(), DeleteQHadron());
440  LeftVector->clear();
441  }
442  //delete LeftVector; // @@ Valgrind ?
443  if(RightVector->size())
444  {
445  std::for_each(RightVector->begin(), RightVector->end(), DeleteQHadron());
446  RightVector->clear();
447  }
448  //delete RightVector; // @@ Valgrind ?
449  inner_sucess=true; // set false on failure
450  while (!StopFragmentation()) // LOOP with break
451  { // Split current string into hadron + new string state
452 #ifdef debug
453  G4cout<<"G4QString::FragmentString:-->Begin LOOP/LOOP, decaying="<<decaying<<G4endl;
454 #endif
455  G4QHadron* Hadron=Splitup(QL); // MAIN: where the hadron is split from the string
456 #ifdef edebug
457  cm4M-=Hadron->Get4Momentum();
458  G4cout<<"-EMC-G4QString::FragmentString:LOOP/LOOP,d4M="<<cm4M-Get4Momentum()<<G4endl;
459 #endif
460  G4bool canBeFrag=IsFragmentable();
461 #ifdef debug
462  G4cout<<"G4QString::FragmentString: LOOP/LOOP, canBeFrag="<<canBeFrag<<", decay="
463  <<decaying<<", H="<<Hadron<<", newStringMass="<<Get4Momentum().m()<<G4endl;
464 #endif
465  if(Hadron && canBeFrag)
466  {
467 #ifdef debug
468  G4cout<<">>G4QString::FragmentString: LOOP/LOOP-NO FRAGM-, dec="<<decaying<<G4endl;
469 #endif
470  if(GetDecayDirection()>0) LeftVector->push_back(Hadron);
471  else RightVector->push_back(Hadron);
472  }
473  else
474  {
475  // @@ Try to convert to the resonance and decay, abandon if M<mGS+mPI0
476  // abandon ... start from the beginning
477 #ifdef debug
478  G4cout<<"G4QString::FragmentString: LOOP/LOOP, Start from scratch"<<G4endl;
479 #endif
480  if (Hadron) delete Hadron;
481  inner_sucess=false;
482  break;
483  }
484 #ifdef debug
485  G4cout<<"G4QString::FragmentString: LOOP/LOOP End, nP="
486  <<thePartons.size()<<", beg="<<(*thePartons.begin())->GetPDGCode()
487  <<",end="<<(*(thePartons.end()-1))->GetPDGCode()<<G4endl;
488 #endif
489  }
490 #ifdef edebug
491  G4LorentzVector fLR=cmS4M-Get4Momentum();
492  for(unsigned L = 0; L < LeftVector->size(); L++)
493  {
494  G4QHadron* LH = (*LeftVector)[L];
495  G4LorentzVector L4M=LH->Get4Momentum();
496  fLR-=L4M;
497  G4cout<<"-EMC-G4QStr::FrStr:L#"<<L<<",PDG="<<LH->GetPDGCode()<<",4M="<<L4M<<G4endl;
498  }
499  for(unsigned R = 0; R < RightVector->size(); R++)
500  {
501  G4QHadron* RH = (*RightVector)[R];
502  G4LorentzVector R4M=RH->Get4Momentum();
503  fLR-=R4M;
504  G4cout<<"-EMC-G4QStr::FrStr:R#"<<R<<",PDG="<<RH->GetPDGCode()<<",4M="<<R4M<<G4endl;
505  }
506  G4cout<<"-EMC-G4QString::FragmentString:L/R_BeforLast->r4M/M2="<<fLR<<fLR.m2()<<G4endl;
507 #endif
508  // Split current string into 2 final Hadrons
509 #ifdef debug
510  G4cout<<"G4QString::FragmentString: inner_success = "<<inner_sucess<<G4endl;
511 #endif
512  if(inner_sucess)
513  {
514  success=true; // Default prototype
515  //... perform last cluster decay
516  G4LorentzVector tot4M = Get4Momentum();
517  G4double totM = tot4M.m();
518 #ifdef debug
519  G4cout<<"G4QString::FragmString: string4M="<<tot4M<<totM<<G4endl;
520 #endif
521  G4QHadron* LeftHadron;
522  G4QHadron* RightHadron;
523  G4QParton* RQuark = 0;
524  SetLeftPartonStable(); // to query quark contents
525  if(DecayIsQuark() && StableIsQuark()) // There're quarks on clusterEnds
526  {
527 #ifdef debug
528  G4cout<<"G4QString::FragmentString: LOOP Quark Algorithm"<<G4endl;
529 #endif
530  LeftHadron= QuarkSplitup(GetLeftParton(), RQuark);
531  }
532  else
533  {
534 #ifdef debug
535  G4cout<<"G4QString::FragmentString: LOOP Di-Quark Algorithm"<<G4endl;
536 #endif
537  //... there is a Diquark on cluster ends
538  G4int IsParticle;
539  if(StableIsQuark()) IsParticle=(GetLeftParton()->GetPDGCode()>0)?-1:1;
540  else IsParticle=(GetLeftParton()->GetPDGCode()>0)?1:-1;
541  G4QPartonPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks
542  RQuark = QuarkPair.GetParton2();
543  G4QParton* LQuark = QuarkPair.GetParton1();
544  LeftHadron = CreateHadron(LQuark, GetLeftParton()); // Create Left Hadron
545  delete LQuark; // Delete the temporaryParton
546  }
547  RightHadron = CreateHadron(GetRightParton(), RQuark); // Create Right Hadron
548  delete RQuark; // Delete the temporaryParton
549  //... repeat procedure, if mass of cluster is too low to produce hadrons
550  G4double LhM=LeftHadron->GetMass();
551  G4double RhM=RightHadron->GetMass();
552 #ifdef debug
553  G4cout<<"G4QStr::FrSt:L="<<LeftHadron->GetPDGCode()<<",R="<<RightHadron->GetPDGCode()
554  <<",ML="<<LhM<<",MR="<<RhM<<",SumM="<<LhM+RhM<<",tM="<<totM<<G4endl;
555 #endif
556  if(totM < LhM + RhM) success=false;
557  //... compute hadron momenta and energies
558  if(success)
559  {
560  G4ThreeVector Pos=GetPosition();
561  G4LorentzVector Lh4M(0.,0.,0.,LhM);
562  G4LorentzVector Rh4M(0.,0.,0.,RhM);
563  if(G4QHadron(tot4M).DecayIn2(Lh4M,Rh4M))
564  {
565  LeftVector->push_back(new G4QHadron(LeftHadron, 0, Pos, Lh4M));
566  delete LeftHadron;
567  RightVector->push_back(new G4QHadron(RightHadron, 0, Pos, Rh4M));
568  delete RightHadron;
569  }
570 #ifdef debug
571  G4cout<<"->>G4QStr::FragString:HFilled (L) PDG="<<LeftHadron->GetPDGCode()<<", 4M="
572  <<Lh4M<<", (R) PDG="<<RightHadron->GetPDGCode()<<", 4M="<<Rh4M<<G4endl;
573 #endif
574 #ifdef edebug
575  G4cout<<"-EMC-G4QString::FragmentString: Residual4M="<<tot4M-Lh4M-Rh4M<<G4endl;
576 #endif
577  }
578  else
579  {
580  if(LeftHadron) delete LeftHadron;
581  if(RightHadron) delete RightHadron;
582  }
583  } // End of inner success
584  } // End of while
585  delete theStringInCMS;
586 #ifdef debug
587  G4cout<<"G4QString::FragmentString: LOOP/LOOP, success="<<success<<G4endl;
588 #endif
589  if (!success)
590  {
591  if(RightVector->size())
592  {
593  std::for_each(RightVector->begin(), RightVector->end(), DeleteQHadron());
594  RightVector->clear();
595  }
596  delete RightVector;
597  if(LeftVector->size())
598  {
599  std::for_each(LeftVector->begin(), LeftVector->end(), DeleteQHadron());
600  LeftVector->clear();
601  }
602  delete LeftVector;
603 #ifdef debug
604  G4cout<<"G4QString::FragmString:StringNotFragm,L4M="<<left4M<<",R4M="<<right4M<<G4endl;
605 #endif
606  // Recover the Left/Right partons 4-moms of the String in ZLS
607  GetLeftParton()->SetPDGCode(leftPDG);
608  GetRightParton()->SetPDGCode(rightPDG);
609  GetLeftParton()->Set4Momentum(left4M);
610  GetRightParton()->Set4Momentum(right4M);
611  return 0; // The string can not be fragmented
612  }
613  // @@@@@ Print collected Left and Right Hadrons (decay resonances!)
614 #ifdef edebug
615  G4LorentzVector sLR=cmS4M;
616  for(unsigned L = 0; L < LeftVector->size(); L++)
617  {
618  G4QHadron* LH = (*LeftVector)[L];
619  G4LorentzVector L4M=LH->Get4Momentum();
620  sLR-=L4M;
621  G4cout<<"-EMC-G4QStr::FragmStri:L#"<<L<<",PDG="<<LH->GetPDGCode()<<",4M="<<L4M<<G4endl;
622  }
623  for(unsigned R = 0; R < RightVector->size(); R++)
624  {
625  G4QHadron* RH = (*RightVector)[R];
626  G4LorentzVector R4M=RH->Get4Momentum();
627  sLR-=R4M;
628  G4cout<<"-EMC-G4QStr::FragmStri:R#"<<R<<",PDG="<<RH->GetPDGCode()<<",4M="<<R4M<<G4endl;
629  }
630  G4cout<<"-EMC-G4QString::FragmentString:---L/R_BeforeMerge---> Res4M="<<sLR<<G4endl;
631 #endif
632  // Join Left and Right Vectors into LeftVector in correct order.
633  while(!RightVector->empty())
634  {
635  LeftVector->push_back(RightVector->back());
636  RightVector->erase(RightVector->end()-1);
637  }
638  delete RightVector;
639  // @@ A trick, the real bug should be found !!
640  G4QHadronVector::iterator ilv; // @@
641  for(ilv = LeftVector->begin(); ilv < LeftVector->end(); ilv++) // @@
642  {
643  G4ThreeVector CV=(*ilv)->Get4Momentum().vect(); // @@
644  if(CV.x()==0. && CV.y()==0. && CV.z()==0.) LeftVector->erase(ilv); // @@
645  }
646  // Calculate time and position of hadrons with @@ very rough formation time
647  G4double StringMass=Get4Momentum().mag();
648  static const G4double dkappa = 2.0 * GeV/fermi; // @@ 2*kappa kappa=1 GeV/fermi (?)
649  for(unsigned c1 = 0; c1 < LeftVector->size(); c1++)
650  {
651  G4double SumPz = 0;
652  G4double SumE = 0;
653  for(unsigned c2 = 0; c2 < c1; c2++)
654  {
655  G4LorentzVector hc2M=(*LeftVector)[c2]->Get4Momentum();
656  SumPz += hc2M.pz();
657  SumE += hc2M.e();
658  }
659  G4QHadron* hc1=(*LeftVector)[c1];
660  G4LorentzVector hc1M=hc1->Get4Momentum();
661  G4double HadronE = hc1M.e();
662  G4double HadronPz= hc1M.pz();
663  hc1->SetFormationTime((StringMass-SumPz-SumPz+HadronE-HadronPz)/dkappa);
664  hc1->SetPosition(G4ThreeVector(0,0,(StringMass-SumE-SumE-HadronE+HadronPz)/dkappa));
665  }
666  G4LorentzRotation toObserverFrame(toCms.inverse());
667 #ifdef debug
668  G4cout<<"G4QString::FragmentString: beforeLoop LVsize = "<<LeftVector->size()<<G4endl;
669 #endif
670  for(unsigned C1 = 0; C1 < LeftVector->size(); C1++)
671  {
672  G4QHadron* Hadron = (*LeftVector)[C1];
673  G4LorentzVector Momentum = Hadron->Get4Momentum();
674  Momentum = toObserverFrame*Momentum;
675  Hadron->Set4Momentum(Momentum);
676  G4LorentzVector Coordinate(Hadron->GetPosition(), Hadron->GetFormationTime());
677  Momentum = toObserverFrame*Coordinate;
678  Hadron->SetFormationTime(Momentum.e());
679  Hadron->SetPosition(GetPosition()+Momentum.vect());
680  }
681 #ifdef edebug
682  G4LorentzVector sLA=string4M;
683  for(unsigned L = 0; L < LeftVector->size(); L++)
684  {
685  G4QHadron* LH = (*LeftVector)[L];
686  G4LorentzVector L4M=LH->Get4Momentum();
687  sLA-=L4M;
688  G4cout<<"-EMC-G4QStr::FragmStri:L#"<<L<<",PDG="<<LH->GetPDGCode()<<",4M="<<L4M<<G4endl;
689  }
690  G4cout<<"-EMC-G4QString::FragmentString:---LSAfterMerge&Conv---> Res4M="<<sLA<<G4endl;
691 #endif
692 #ifdef debug
693  G4cout<<"G4QString::FragmentString: *** Done *** "<<G4endl;
694 #endif
695  return LeftVector; // Should be deleted by user (@@ Valgrind complain ?)
696 } // End of FragmentString
697 
698 // Simple decay of the string if the excitation mass is too small for HE fragmentation
699 // !! If the mass is below the single hadron threshold, make warning (improve) and convert
700 // the string to the single S-hadron breaking energy conservation (temporary) and improve,
701 // taking the threshold into account on the level of the String creation (merge strings) !!
702 G4QHadronVector* G4QString::LightFragmentationTest()
703 {
704  // Check string decay threshold
705  G4LorentzVector tot4M=Get4Momentum();
706 #ifdef debug
707  G4cout<<"G4QString::LightFragmentationTest: ***Called***, string4M="<<tot4M<<G4endl;
708 #endif
709  G4QHadronVector* result=0; // return 0 when string exceeds the mass cut or below mh1+mh2
710 
711  G4QHadronPair hadrons((G4QHadron*)0, (G4QHadron*)0); // pair of hadrons for output of FrM
712  G4double fragMass = FragmentationMass(0,&hadrons); // Minimum mass to decay the string
713 #ifdef debug
714  G4cout<<"G4QString::LightFragTest: before check nP="<<thePartons.size()<<", MS2="
715  <<Mass2()<<", MCut="<<MassCut<<", beg="<<(*thePartons.begin())->GetPDGCode()
716  <<",end="<<(*(thePartons.end()-1))->GetPDGCode()<<", fM="<<fragMass<<G4endl;
717 #endif
718  if(Mass2() > sqr(fragMass+MassCut))// Big enough to fragment in a lader (avoid the decay)
719  {
720  if(hadrons.first) delete hadrons.first;
721  if(hadrons.second) delete hadrons.second;
722 #ifdef debug
723  G4cout<<"G4QString::LightFragTest:NO,M2="<<Mass2()<<">"<<sqr(fragMass+MassCut)<<G4endl;
724 #endif
725  return result; // =0. Depends on the parameter of the Mass Cut
726  }
727  G4double totM= tot4M.m();
728  G4QHadron* h1=hadrons.first;
729  G4QHadron* h2=hadrons.second;
730  if(h1 && h2)
731  {
732  G4double h1M = h1->GetMass();
733  G4double h2M = h2->GetMass();
734 #ifdef debug
735  G4cout<<"G4QString::LightFragTest:tM="<<totM<<","<<h1M<<"+"<<h2M<<"+"<<h1M+h2M<<G4endl;
736 #endif
737  if(h1M + h2M <= totM) // The string can decay in these two hadrons
738  {
739  // Create two stable hadrons
740  G4LorentzVector h4M1(0.,0.,0.,h1M);
741  G4LorentzVector h4M2(0.,0.,0.,h2M);
742  if(G4QHadron(tot4M).DecayIn2(h4M1,h4M2))
743  {
744  result = new G4QHadronVector;
745  result->push_back(new G4QHadron(hadrons.first, 0, GetPosition(), h4M1));
746  result->push_back(new G4QHadron(hadrons.second,0, GetPosition(), h4M2));
747  }
748 #ifdef edebug
749  G4int L=result->size(); if(L) for(G4int i=0; i<L; i++)
750  {
751  tot4M-=(*result)[i]->Get4Momentum();
752  G4cout<<"-EMC-G4QString::LightFragTest: i="<<i<<", residual4M="<<tot4M<<G4endl;
753  }
754 #endif
755  }
756 #ifdef debug
757  else G4cout<<"-Warning-G4QString::LightFragTest: TooBigHadronMasses to decay"<<G4endl;
758 #endif
759  }
760 #ifdef debug
761  else G4cout<<"-Warning-G4QString::LightFragTest: No Hadrons have been proposed"<<G4endl;
762 #endif
763  delete hadrons.first;
764  delete hadrons.second;
765  return result;
766 } // End of LightFragmentationTest
767 
768 // Calculate Fragmentation Mass (if pdefs # 0, returns two hadrons)
769 G4double G4QString::FragmentationMass(G4int HighSpin, G4QHadronPair* pdefs)
770 {
771  G4double mass=0.;
772 #ifdef debug
773  G4cout<<"G4QString::FragmMass: ***Called***, s4M="<<Get4Momentum()<<G4endl;
774 #endif
775  // Example how to use an interface to different member functions
776  G4QHadron* Hadron1 = 0;
777  G4QHadron* Hadron2 = 0;
778 #ifdef debug
779  G4cout<<"G4QString::FragmentationMass: Create spin-0 or spin-1/2 hadron: nP="
780  <<thePartons.size()<<", beg="<<(*thePartons.begin())->GetPDGCode()
781  <<",end="<<(*(thePartons.end()-1))->GetPDGCode()<<G4endl;
782 #endif
783  G4int iflc = (G4UniformRand() < 0.5) ? 1 : 2; // Create additional Q-antiQ pair @@ No S
784  G4int LPDG= GetLeftParton()->GetPDGCode();
785  G4int LT = GetLeftParton()->GetType();
786  if ( (LPDG > 0 && LT == 1) || (LPDG < 0 && LT == 2) ) iflc = -iflc; // anti-quark
787  G4QParton* piflc = new G4QParton( iflc);
788  G4QParton* miflc = new G4QParton(-iflc);
789  if(HighSpin)
790  {
791  Hadron1 = CreateHighSpinHadron(GetLeftParton(),piflc);
792  Hadron2 = CreateHighSpinHadron(GetRightParton(),miflc);
793 #ifdef debug
794  G4cout<<"G4QString::FragmentationMass: High, PDG1="<<Hadron1->GetPDGCode()
795  <<", PDG2="<<Hadron2->GetPDGCode()<<G4endl;
796 #endif
797  }
798  else
799  {
800  Hadron1 = CreateLowSpinHadron(GetLeftParton(),piflc);
801  Hadron2 = CreateLowSpinHadron(GetRightParton(),miflc);
802 #ifdef debug
803  G4cout<<"G4QString::FragmentationMass: Low, PDG1="<<Hadron1->GetPDGCode()
804  <<", PDG2="<<Hadron2->GetPDGCode()<<G4endl;
805 #endif
806  }
807  mass = Hadron1->GetMass() + Hadron2->GetMass();
808  if(pdefs) // need to return hadrons as well as the mass estimate
809  {
810  pdefs->first = Hadron1; // To be deleted by the calling program if not zero
811  pdefs->second = Hadron2; // To be deleted by the calling program if not zero
812  }
813  else // Forget about the hadrons
814  {
815  if(Hadron1) delete Hadron1;
816  if(Hadron2) delete Hadron2;
817  }
818  delete piflc;
819  delete miflc;
820 #ifdef debug
821  G4cout<<"G4QString::FragmentationMass: ***Done*** mass="<<mass<<G4endl;
822 #endif
823  return mass;
824 } // End of FragmentationMass
825 
826 void G4QString::SetLeftPartonStable()
827 {
828  theStableParton=GetLeftParton();
829  theDecayParton=GetRightParton();
830  decaying=Right;
831 }
832 
833 void G4QString::SetRightPartonStable()
834 {
835  theStableParton=GetRightParton();
836  theDecayParton=GetLeftParton();
837  decaying=Left;
838 }
839 
840 G4int G4QString::GetDecayDirection() const
841 {
842  if (decaying == Left ) return +1;
843  else if (decaying == Right) return -1;
844  else
845  {
846  G4cerr<<"***G4QString::GetDecayDirection: wrong DecayDirection="<<decaying<<G4endl;
847  G4Exception("G4QString::GetDecayDirection:","72",FatalException,"WrongDecayDirection");
848  }
849  return 0;
850 }
851 
852 //G4ThreeVector G4QString::StablePt()
853 //{
854 // if (decaying == Left ) return Ptright;
855 // else if (decaying == Right ) return Ptleft;
856 // else
857 // {
858 // G4cerr<<"***G4QString::StablePt: wrong DecayDirection="<<decaying<<G4endl;
859 // G4Exception("G4QString::StablePt:","72",FatalException,"WrongDecayDirection");
860 // }
861 // return G4ThreeVector();
862 //}
863 
864 G4ThreeVector G4QString::DecayPt()
865 {
866  if (decaying == Left ) return Ptleft;
867  else if (decaying == Right ) return Ptright;
868  else
869  {
870  G4cerr<<"***G4QString::DecayPt: wrong DecayDirection="<<decaying<<G4endl;
871  G4Exception("G4QString::DecayPt:","72",FatalException,"WrongDecayDirection");
872  }
873  return G4ThreeVector();
874 }
875 
876 // Random choice of string end to use it for creating the hadron (decay)
877 G4QHadron* G4QString::Splitup(G4bool QL)
878 {
879  SideOfDecay = (G4UniformRand() < 0.5) ? 1: -1;
880 #ifdef debug
881  G4cout<<"G4QString::Splitup:**Called**,s="<<SideOfDecay<<",s4M="<<Get4Momentum()<<G4endl;
882 #endif
883  if(SideOfDecay<0) SetLeftPartonStable(); // Decay Right parton
884  else SetRightPartonStable(); // Decay Left parton
885  G4QParton* newStringEnd;
886  G4QHadron* Hadron;
887  if(DecayIsQuark()) Hadron= QuarkSplitup(theDecayParton, newStringEnd); // Split Quark
888  else Hadron= DiQuarkSplitup(theDecayParton, newStringEnd); // Split DiQuark
889 #ifdef debug
890  G4cout<<"G4QString::Splitup: newStringEndPDG="<<newStringEnd->GetPDGCode()<<", nP="
891  <<thePartons.size()<<", beg="<<(*thePartons.begin())->GetPDGCode()
892  <<",end="<<(*(thePartons.end()-1))->GetPDGCode()<<G4endl;
893 #endif
894  // create new String from the old one: keep Left and Right order, but replace decay
895  G4LorentzVector* HadronMomentum=SplitEandP(Hadron, QL);//The decayed parton isn't changed
896 #ifdef debug
897  G4cout<<"G4QString::Splitup: HM="<<HadronMomentum<<", nP="
898  <<thePartons.size()<<", beg="<<(*thePartons.begin())->GetPDGCode()
899  <<",end="<<(*(thePartons.end()-1))->GetPDGCode()<<G4endl;
900 #endif
901  if(HadronMomentum) // The decay succeeded, now the new 4-mon can be set to NewStringEnd
902  {
903 #ifdef pdebug
904  G4cout<<"---->>G4QString::Splitup: HFilled 4M="<<*HadronMomentum<<",PDG="
905  <<Hadron->GetPDGCode()<<",s4M-h4M="<<Get4Momentum()-*HadronMomentum<<G4endl;
906 #endif
907  newStringEnd->Set4Momentum(theDecayParton->Get4Momentum()-*HadronMomentum);
908  Hadron->Set4Momentum(*HadronMomentum);
909  Hadron->SetPosition(GetPosition());
910  if(decaying == Left)
911  {
912  G4QParton* theFirst = thePartons[0]; // Substitute for the First Parton
913  delete theFirst; // The OldParton instance is deleted
914  thePartons[0] = newStringEnd; // Delete equivalent for newStringEnd
915 #ifdef debug
916  G4cout<<"G4QString::Splitup: theFirstPDG="<<theFirst->GetPDGCode()<<G4endl;
917 #endif
918  Ptleft -= HadronMomentum->vect();
919  Ptleft.setZ(0.); // @@ (Z is anyway ignored) M.K. (?)
920  }
921  else if (decaying == Right)
922  {
923  G4QParton* theLast = thePartons[thePartons.size()-1]; // Substitute for theLastParton
924  delete theLast; // The OldParton instance is deleted
925  thePartons[thePartons.size()-1] = newStringEnd; // Delete equivalent for newStringEnd
926 #ifdef debug
927  G4cout<<"G4QString::Splitup: theLastPDG="<<theLast->GetPDGCode()<<", nP="
928  <<thePartons.size()<<", beg="<<thePartons[0]->GetPDGCode()<<",end="
929  <<thePartons[thePartons.size()-1]->GetPDGCode()<<",P="<<theLast
930  <<"="<<thePartons[thePartons.size()-1]<<G4endl;
931 #endif
932  Ptright -= HadronMomentum->vect();
933  Ptright.setZ(0.); // @@ (Z is anyway ignored) M.K. (?)
934  }
935  else
936  {
937  G4cerr<<"***G4QString::Splitup: wrong oldDecay="<<decaying<<G4endl;
938  G4Exception("G4QString::Splitup","72",FatalException,"WrongDecayDirection");
939  }
940  Pplus -= HadronMomentum->e() + HadronMomentum->pz();// Reduce Pplus ofTheString (Left)
941  Pminus -= HadronMomentum->e() - HadronMomentum->pz();// Reduce Pminus ofTheString(Rite)
942 #ifdef debug
943  G4cout<<"G4QString::Splitup: P+="<<Pplus<<",P-="<<Pminus<<", nP="
944  <<thePartons.size()<<", beg="<<(*thePartons.begin())->GetPDGCode()
945  <<",end="<<(*(thePartons.end()-1))->GetPDGCode()<<G4endl;
946  G4cout<<">...>G4QString::Splitup: NewString4M="<<Get4Momentum()<<G4endl;
947 #endif
948  delete HadronMomentum;
949  }
950 #ifdef debug
951  G4cout<<"G4QString::Splitup: ***Done*** H4M="<<Hadron->Get4Momentum()<<", nP="
952  <<thePartons.size()<<", beg="<<(*thePartons.begin())->GetPDGCode()
953  <<",end="<<(*(thePartons.end()-1))->GetPDGCode()<<G4endl;
954 #endif
955  return Hadron;
956 } // End of Splitup
957 
958 // QL=true for QGSM and QL=false for Lund fragmentation
959 G4LorentzVector* G4QString::SplitEandP(G4QHadron* pHadron, G4bool QL)
960 {
961  G4double HadronMass = pHadron->GetMass();
962 #ifdef debug
963  G4cout<<"G4QString::SplitEandP: ***Called*** HMass="<<HadronMass<<G4endl;
964 #endif
965  // calculate and assign hadron transverse momentum component HadronPx andHadronPy
966  G4ThreeVector HadronPt = SampleQuarkPt() + DecayPt(); // @@ SampleQuarkPt & DecayPt once
967  HadronPt.setZ(0.);
968  //... sample z to define hadron longitudinal momentum and energy
969  //... but first check the available phase space
970  G4double HadronMass2T = HadronMass*HadronMass + HadronPt.mag2();
971  if (HadronMass2T >= SmoothParam*Mass2() ) return 0; // restart!
972  //... then compute allowed z region z_min <= z <= z_max
973  G4double zMin = HadronMass2T/Mass2();
974  G4double zMax = 1.;
975 #ifdef debug
976  G4cout<<"G4QString::SplitEandP: zMin="<<zMin<<", zMax"<<zMax<<G4endl;
977 #endif
978  if (zMin >= zMax) return 0; // have to start all over!
979  G4double z=0;
980  if(QL) z = GetQGSMLightConeZ(zMin, zMax, theDecayParton->GetPDGCode(), pHadron,
981  HadronPt.x(), HadronPt.y());
982  else z = GetLundLightConeZ(zMin, zMax, theDecayParton->GetPDGCode(), pHadron,
983  HadronPt.x(), HadronPt.y());
984  //... now compute hadron longitudinal momentum and energy
985  // longitudinal hadron momentum component HadronPz
986  G4double zl= z;
987  if (decaying == Left ) zl*=Pplus;
988  else if (decaying == Right ) zl*=Pminus;
989  else // @@ Is that possible?
990  {
991  G4cerr<<"***G4QString::SplitEandP: wrong DecayDirection="<<decaying<<G4endl;
992  G4Exception("G4QString::SplitEandP:","72",FatalException,"WrongDecayDirection");
993  }
994  G4double HadronE = (zl + HadronMass2T/zl)/2;
995  HadronPt.setZ( GetDecayDirection() * (zl - HadronE) );
996  G4LorentzVector* a4Momentum= new G4LorentzVector(HadronPt,HadronE);
997  return a4Momentum;
998 }
999 
1000 G4ThreeVector G4QString::SampleQuarkPt()
1001 {
1002  G4double Pt = SigmaQT * std::sqrt( -std::log(G4UniformRand()));
1003  G4double phi = twopi*G4UniformRand();
1004  return G4ThreeVector(Pt * std::cos(phi),Pt * std::sin(phi),0);
1005 }
1006 
1007 G4QHadron* G4QString::QuarkSplitup(G4QParton* decay, G4QParton* &created)// VGComplTo decay
1008 {
1009  G4int IsParticle=(decay->GetPDGCode()>0) ? -1 : +1; // a quark needs antiquark or diquark
1010  G4QPartonPair QuarkPair = CreatePartonPair(IsParticle);
1011  created = QuarkPair.GetParton2(); // New Parton after splitting
1012 #ifdef debug
1013  G4cout<<"G4QString::QuarkSplitup: ***Called*** crP="<<created->GetPDGCode()<<G4endl;
1014 #endif
1015  G4QParton* P1=QuarkPair.GetParton1();
1016  G4QHadron* result=CreateHadron(P1, decay); // New Hadron after splitting
1017  delete P1; // Clean up the temporary parton
1018  return result;
1019 } // End of QuarkSplitup
1020 
1021 //
1022 G4QHadron* G4QString::DiQuarkSplitup(G4QParton* decay, G4QParton* &created)
1023 {
1024  //... can Diquark break or not?
1025  if (G4UniformRand() < DiquarkBreakProb )
1026  {
1027  //... Diquark break
1028  G4int stableQuarkEncoding = decay->GetPDGCode()/1000;
1029  G4int decayQuarkEncoding = (decay->GetPDGCode()/100)%10;
1030  if (G4UniformRand() < 0.5)
1031  {
1032  G4int Swap = stableQuarkEncoding;
1033  stableQuarkEncoding = decayQuarkEncoding;
1034  decayQuarkEncoding = Swap;
1035  }
1036  G4int IsParticle=(decayQuarkEncoding>0) ? -1 : +1;// Diquark is equivalent to antiquark
1037  G4QPartonPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
1038  G4QParton* P2=QuarkPair.GetParton2();
1039  G4int QuarkEncoding=P2->GetPDGCode();
1040  delete P2;
1041  G4int i10 = std::max(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
1042  G4int i20 = std::min(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
1043  G4int spin = (i10 != i20 && G4UniformRand() <= 0.5) ? 1 : 3;
1044  G4int NewDecayEncoding = -1*IsParticle*(i10 * 1000 + i20 * 100 + spin);
1045  created = new G4QParton(NewDecayEncoding);
1046 #ifdef debug
1047  G4cout<<"G4QString::DiQuarkSplitup: inside, crP="<<created->GetPDGCode()<<G4endl;
1048 #endif
1049  G4QParton* decayQuark= new G4QParton(decayQuarkEncoding);
1050  G4QParton* P1=QuarkPair.GetParton1();
1051  G4QHadron* newH=CreateHadron(P1, decayQuark);
1052  delete P1;
1053  delete decayQuark;
1054  return newH;
1055  }
1056  else
1057  {
1058  //... Diquark does not break
1059  G4int IsParticle=(decay->GetPDGCode()>0) ? +1 : -1;
1060  G4QPartonPair QuarkPair = CreatePartonPair(IsParticle,false); // no diquarks wanted
1061  created = QuarkPair.GetParton2();
1062 #ifdef debug
1063  G4cout<<"G4QString::DiQuarkSplitup: diQ not break, crP="<<created->GetPDGCode()<<G4endl;
1064 #endif
1065  G4QParton* P1=QuarkPair.GetParton1();
1066  G4QHadron* newH=CreateHadron(P1, decay);
1067  delete P1;
1068  return newH;
1069  }
1070 } // End of DiQuarkSplitup
1071 
1072 G4QPartonPair G4QString::CreatePartonPair(G4int NeedParticle, G4bool AllowDiquarks)
1073 {
1074 #ifdef debug
1075  G4cout<<"G4QString::CreatePartonPair: ***Called***, P="<<NeedParticle<<", ALLOWdQ="
1076  <<AllowDiquarks<<G4endl;
1077 #endif
1078  // NeedParticle = {+1 for Particle, -1 for AntiParticle}
1079  if(AllowDiquarks && G4UniformRand() < DiquarkSuppress)
1080  {
1081  // Create a Diquark - AntiDiquark pair , first in pair is anti to IsParticle
1082  G4int q1 = SampleQuarkFlavor();
1083  G4int q2 = SampleQuarkFlavor();
1084  G4int spin = (q1 != q2 && G4UniformRand() <= 0.5) ? 1 : 3; // @@ 0.5 M.K.?
1085  // Convention: quark with higher PDG number is first
1086  G4int PDGcode = (std::max(q1,q2) * 1000 + std::min(q1,q2) * 100 + spin) * NeedParticle;
1087 #ifdef debug
1088  G4cout<<"G4QString::CreatePartonPair: Created dQ-AdQ, PDG="<<PDGcode<<G4endl;
1089 #endif
1090  return G4QPartonPair(new G4QParton(-PDGcode), new G4QParton(PDGcode));
1091  }
1092  else
1093  {
1094  // Create a Quark - AntiQuark pair, first in pair is a Particle
1095  G4int PDGcode=SampleQuarkFlavor()*NeedParticle;
1096 #ifdef debug
1097  G4cout<<"G4QString::CreatePartonPair: Created Q-aQ, PDG="<<PDGcode<<G4endl;
1098 #endif
1099  return G4QPartonPair(new G4QParton(PDGcode), new G4QParton(-PDGcode));
1100  }
1101 } // End of CreatePartonPair
1102 
1103 // Creation of the Meson out of two partons (q, anti-q)
1104 G4QHadron* G4QString::CreateMeson(G4QParton* black, G4QParton* white, Spin theSpin)
1105 {
1106  static G4double scalarMesonMixings[6]={0.5, 0.25, 0.5, 0.25, 1.0, 0.5};
1107  static G4double vectorMesonMixings[6]={0.5, 0.0, 0.5, 0.0, 1.0, 1.0};
1108  G4int id1= black->GetPDGCode();
1109  G4int id2= white->GetPDGCode();
1110 #ifdef debug
1111  G4cout<<"G4QString::CreateMeson: bT="<<black->GetType()<<"("<<id1<<"), wT="
1112  <<white->GetType()<<"("<<id2<<")"<<G4endl;
1113 #endif
1114  if (std::abs(id1) < std::abs(id2)) // exchange black and white
1115  {
1116  G4int xchg = id1;
1117  id1 = id2;
1118  id2 = xchg;
1119  }
1120  if(std::abs(id1)>3)
1121  {
1122  G4cerr<<"***G4QString::CreateMeson: q1="<<id1<<", q2="<<id2
1123  <<" while CHIPS is only SU(3)"<<G4endl;
1124  G4Exception("G4QString::CreateMeson:","72",FatalException,"HeavyQuarkFound");
1125  }
1126  G4int PDGEncoding=0;
1127  if(!(id1+id2)) // annihilation case (neutral)
1128  {
1129  G4double rmix = G4UniformRand();
1130  G4int imix = 2*std::abs(id1) - 1;
1131  if(theSpin == SpinZero)
1132  PDGEncoding = 110*(1 + G4int(rmix + scalarMesonMixings[imix - 1])
1133  + G4int(rmix + scalarMesonMixings[imix] ) ) + theSpin;
1134  else
1135  PDGEncoding = 110*(1 + G4int(rmix + vectorMesonMixings[imix - 1])
1136  + G4int(rmix + vectorMesonMixings[imix] ) ) + theSpin;
1137  }
1138  else
1139  {
1140  PDGEncoding = 100 * std::abs(id1) + 10 * std::abs(id2) + theSpin;
1141  G4bool IsUp = (std::abs(id1)&1) == 0; // quark 1 is up type quark (u or c?)
1142  G4bool IsAnti = id1 < 0; // quark 1 is an antiquark?
1143  if( (IsUp && IsAnti) || (!IsUp && !IsAnti) ) PDGEncoding = - PDGEncoding;
1144  // Correction for the true neutral mesons
1145  if( PDGEncoding == -111 || PDGEncoding == -113 || PDGEncoding == -223 ||
1146  PDGEncoding == -221 || PDGEncoding == -331 || PDGEncoding == -333 )
1147  PDGEncoding = - PDGEncoding;
1148  }
1149  G4QHadron* Meson= new G4QHadron(PDGEncoding);
1150 #ifdef debug
1151  G4cout<<"G4QString::CreateBaryon: Meson is created with PDG="<<PDGEncoding<<G4endl;
1152 #endif
1153  //delete black; // It is better to delete here and consider
1154  //delete white; // the hadron creation as a delete equivalent
1155  return Meson;
1156 }
1157 
1158 // Creation of the Baryon out of two partons (q, di-q), (anti-q, anti-di-q)
1159 G4QHadron* G4QString::CreateBaryon(G4QParton* black, G4QParton* white, Spin theSpin)
1160 {
1161  G4int id1= black->GetPDGCode();
1162  G4int id2= white->GetPDGCode();
1163 #ifdef debug
1164  G4cout<<"G4QString::CreateBaryon: bT="<<black->GetType()<<"("<<id1<<"), wT="
1165  <<white->GetType()<<"("<<id2<<")"<<G4endl;
1166 #endif
1167  if(std::abs(id1) < std::abs(id2))
1168  {
1169  G4int xchg = id1;
1170  id1 = id2;
1171  id2 = xchg;
1172  }
1173  if(std::abs(id1)<1000 || std::abs(id2)> 3)
1174  {
1175  G4cerr<<"***G4QString::CreateBaryon: q1="<<id1<<", q2="<<id2
1176  <<" can't create a Baryon"<<G4endl;
1177  G4Exception("G4QString::CreateBaryon:","72",FatalException,"WrongQdQSequence");
1178  }
1179  G4int ifl1= std::abs(id1)/1000;
1180  G4int ifl2 = (std::abs(id1) - ifl1 * 1000)/100;
1181  G4int diquarkSpin = std::abs(id1)%10;
1182  G4int ifl3 = id2;
1183  if (id1 < 0) {ifl1 = - ifl1; ifl2 = - ifl2;}
1184  //... Construct baryon, distinguish Lambda and Sigma baryons.
1185  G4int kfla = std::abs(ifl1);
1186  G4int kflb = std::abs(ifl2);
1187  G4int kflc = std::abs(ifl3);
1188  G4int kfld = std::max(kfla,kflb);
1189  kfld = std::max(kfld,kflc);
1190  G4int kflf = std::min(kfla,kflb);
1191  kflf = std::min(kflf,kflc);
1192  G4int kfle = kfla + kflb + kflc - kfld - kflf;
1193  //... baryon with content uuu or ddd or sss has always spin = 3/2
1194  if(kfla==kflb && kflb==kflc) theSpin=SpinThreeHalf;
1195 
1196  G4int kfll = 0;
1197  if(theSpin == SpinHalf && kfld > kfle && kfle > kflf)
1198  {
1199  // Spin J=1/2 and all three quarks different
1200  // Two states exist: (uds -> lambda or sigma0)
1201  // - lambda: s(ud)0 s : 3122; ie. reverse the two lighter quarks
1202  // - sigma0: s(ud)1 s : 3212
1203  if(diquarkSpin == 1 )
1204  {
1205  if ( kfla == kfld) kfll = 1; // heaviest quark in diquark
1206  else kfll = G4int(0.25 + G4UniformRand());
1207  }
1208  if(diquarkSpin==3 && kfla!=kfld) kfll = G4int(0.75+G4UniformRand());
1209  }
1210  G4int PDGEncoding=0;
1211  if (kfll == 1) PDGEncoding = 1000 * kfld + 100 * kflf + 10 * kfle + theSpin;
1212  else PDGEncoding = 1000 * kfld + 100 * kfle + 10 * kflf + theSpin;
1213  if (id1 < 0) PDGEncoding = -PDGEncoding;
1214  G4QHadron* Baryon= new G4QHadron(PDGEncoding);
1215 #ifdef debug
1216  G4cout<<"G4QString::CreateBaryon: Baryon is created with PDG="<<PDGEncoding<<G4endl;
1217 #endif
1218  //delete black; // It is better to delete here and consider
1219  //delete white; // the hadron creation as a delete equivalent
1220  return Baryon;
1221 } // End of CreateBaryon
1222 
1223 G4QHadron* G4QString::CreateHadron(G4QParton* black, G4QParton* white)
1224 {
1225  //static G4double mesonLowSpin = 0.25; // probability to create scalar meson (2s+1)
1226  //static G4double baryonLowSpin= 1./3.; // probability to create 1/2 baryon (2s+1)
1227  static G4double mesonLowSpin = 0.5; // probability to create scalar meson (spFlip)
1228  static G4double baryonLowSpin= 0.5; // probability to create 1/2 baryon (spinFlip)
1229  G4int bT=black->GetType();
1230  G4int wT=white->GetType();
1231 #ifdef debug
1232  G4cout<<"G4QString::CreateHadron: bT="<<bT<<"("<<black->GetPDGCode()<<"), wT="<<wT<<"("
1233  <<white->GetPDGCode()<<")"<<G4endl;
1234 #endif
1235  if(bT==2 || wT==2)
1236  {
1237  // Baryon consists of quark and at least one di-quark
1238  Spin spin = (G4UniformRand() < baryonLowSpin) ? SpinHalf : SpinThreeHalf;
1239 #ifdef debug
1240  G4cout<<"G4QString::CreateHadron: ----> Baryon is under creation"<<G4endl;
1241 #endif
1242  return CreateBaryon(black, white, spin);
1243  }
1244  else
1245  {
1246  // Meson consists of quark and abnti-quark
1247  Spin spin = (G4UniformRand() < mesonLowSpin) ? SpinZero : SpinOne;
1248 #ifdef debug
1249  G4cout<<"G4QString::CreateHadron: ----> Meson is under creation"<<G4endl;
1250 #endif
1251  return CreateMeson(black, white, spin);
1252  }
1253 } // End of Create Hadron
1254 
1255 // Creation of only High Spin (2,3/2) hadrons
1256 G4QHadron* G4QString::CreateLowSpinHadron(G4QParton* black, G4QParton* white)
1257 {
1258  G4int bT=black->GetType();
1259  G4int wT=white->GetType();
1260 #ifdef debug
1261  G4cout<<"G4QString::CreateLowSpinHadron: ***Called***, bT="<<bT<<"("<<black->GetPDGCode()
1262  <<"), wT="<<wT<<"("<<white->GetPDGCode()<<")"<<G4endl;
1263 #endif
1264  if(bT == 1 && wT == 1)
1265  {
1266 #ifdef debug
1267  G4cout<<"G4QString::CreateLowSpinHadron: ----> Meson is under creation"<<G4endl;
1268 #endif
1269  return CreateMeson(black, white, SpinZero);
1270  }
1271  else // returns a SpinThreeHalf Baryon if all quarks are the same
1272  {
1273 #ifdef debug
1274  G4cout<<"G4QString::CreateLowSpinHadron: ----> Baryon is under creation"<<G4endl;
1275 #endif
1276  return CreateBaryon(black, white, SpinHalf);
1277  }
1278 } // End of CreateLowSpinHadron
1279 
1280 // Creation of only High Spin (2,3/2) hadrons
1281 G4QHadron* G4QString::CreateHighSpinHadron(G4QParton* black, G4QParton* white)
1282 {
1283  G4int bT=black->GetType();
1284  G4int wT=white->GetType();
1285 #ifdef debug
1286  G4cout<<"G4QString::CreateHighSpinHadron:***Called***, bT="<<bT<<"("<<black->GetPDGCode()
1287  <<"), wT="<<wT<<"("<<white->GetPDGCode()<<")"<<G4endl;
1288 #endif
1289  if(bT == 1 && wT == 1)
1290  {
1291 #ifdef debug
1292  G4cout<<"G4QString::CreateHighSpinHadron: ----> Meson is created"<<G4endl;
1293 #endif
1294  return CreateMeson(black,white, SpinOne);
1295  }
1296  else
1297  {
1298 #ifdef debug
1299  G4cout<<"G4QString::CreateHighSpinHadron: ----> Baryon is created"<<G4endl;
1300 #endif
1301  return CreateBaryon(black,white,SpinThreeHalf);
1302  }
1303 } // End of CreateHighSpinHadron
Generated on Sat May 25 2013 14:33:47 for Geant4 by   doxygen 1.8.4