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

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

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