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

Last change on this file since 967 was 962, checked in by garnier, 15 years ago

update processes
File size: 30.8 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.4 2007/07/06 07:38:36 mkossov Exp $
28// GEANT4 tag $Name: geant4-09-02-ref-02 $
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
41//#define debug
42
43#include "G4QString.hh"
44#include <algorithm>
45// Static parameters definition
46G4double G4QString::MassCut=350.*MeV;     // minimum mass cut for the string
47G4double G4QString::ClusterMass=150.*MeV; // minimum cluster mass for fragmentation
48G4double G4QString::SigmaQT=0.5*GeV;      // quarkTransverseMomentum distribution parameter
49G4double G4QString::DiquarkSuppress=0.1;  // is Diquark suppression parameter 
50G4double G4QString::DiquarkBreakProb=0.1; // is Diquark breaking probability
51G4double G4QString::SmoothParam=0.9;      // QGS model parameter
52G4double G4QString::StrangeSuppress=0.435;// Strangeness suppression (u:d:s=1:1:0.3 ?M.K.)
53G4double G4QString::widthOfPtSquare=-0.72*GeV*GeV; // pt -width2 forStringExcitation
54G4int G4QString::StringLoopInterrupt=999; // String fragmentation LOOP limit
55G4int G4QString::ClusterLoopInterrupt=500;// Cluster fragmentation LOOP limit
56
57G4QString::G4QString() : theDirection(0),thePosition(G4ThreeVector(0.,0.,0.)),
58 hadronizer(new G4QHadronBuilder){}
59
60G4QString::G4QString(G4QParton* Color, G4QParton* AntiColor, G4int Direction) :
61  hadronizer(new G4QHadronBuilder)
62{
63  ExciteString(Color, AntiColor, Direction);
64}
65
66G4QString::G4QString(G4QParton* QCol,G4QParton* Gluon,G4QParton* QAntiCol,G4int Direction):
67         theDirection(Direction),thePosition(QCol->GetPosition()),hadronizer(new G4QHadronBuilder)
68{
69  thePartons.push_back(QCol);
70  thePartons.push_back(Gluon);
71  thePartons.push_back(QAntiCol);
72}
73
74G4QString::G4QString(const G4QString &right) : theDirection(right.GetDirection()),
75thePosition(right.GetPosition()), hadronizer(new G4QHadronBuilder){}
76
77G4QString::~G4QString()
78 {if(thePartons.size())std::for_each(thePartons.begin(),thePartons.end(),DeleteQParton());}
79
80G4LorentzVector G4QString::Get4Momentum() const
81{
82         G4LorentzVector momentum(0.,0.,0.,0.);
83         for(unsigned i=0; i<thePartons.size(); i++) momentum += thePartons[i]->Get4Momentum();
84         return momentum;
85}
86
87void G4QString::LorentzRotate(const G4LorentzRotation & rotation)
88{
89         for(unsigned i=0; i<thePartons.size(); i++)
90            thePartons[i]->Set4Momentum(rotation*thePartons[i]->Get4Momentum());
91}
92
93void G4QString::InsertParton(G4QParton* aParton, const G4QParton* addafter)
94{
95  G4QPartonVector::iterator insert_index;   // Begin by default (? M.K.)
96       
97         if(addafter != NULL) 
98         {
99           insert_index=std::find(thePartons.begin(), thePartons.end(), addafter);
100           if (insert_index == thePartons.end())                // No object addafter in thePartons
101           {
102                                        G4cerr<<"***G4QString::InsertParton: Address Parton is not found"<<G4endl;
103      G4Exception("G4QString::InsertParton:","72",FatalException,"NoAddressParton");
104           }
105         }
106         thePartons.insert(insert_index+1, aParton);
107} 
108
109G4LorentzRotation G4QString::TransformToCenterOfMass()
110{
111         G4LorentzRotation toCms(-1*Get4Momentum().boostVector());
112         for(unsigned i=0; i<thePartons.size(); i++)
113                             thePartons[i]->Set4Momentum(toCms*thePartons[i]->Get4Momentum());
114         return toCms;
115}
116
117G4LorentzRotation G4QString::TransformToAlignedCms()
118{
119         G4LorentzVector momentum=Get4Momentum();
120         G4LorentzRotation toAlignedCms(-1*momentum.boostVector());
121         momentum= toAlignedCms* thePartons[0]->Get4Momentum();
122         toAlignedCms.rotateZ(-1*momentum.phi());
123         toAlignedCms.rotateY(-1*momentum.theta());
124         for(unsigned index=0; index<thePartons.size(); index++)
125         {
126           momentum=toAlignedCms * thePartons[index]->Get4Momentum();
127           thePartons[index]->Set4Momentum(momentum);
128         }
129         return toAlignedCms;
130}
131
132void G4QString::Boost(G4ThreeVector& Velocity)
133{
134  for(unsigned cParton=0; cParton<thePartons.size(); cParton++)
135  {
136    G4LorentzVector Mom = thePartons[cParton]->Get4Momentum();
137    Mom.boost(Velocity);
138    thePartons[cParton]->Set4Momentum(Mom);
139  }
140}
141
142G4QParton* G4QString::GetColorParton(void) const
143{
144  G4QParton* start = *(thePartons.begin());
145  G4QParton* end = *(thePartons.end()-1);
146  G4int Encoding = start->GetPDGCode();
147  if (Encoding<-1000 || (Encoding  < 1000 && Encoding > 0)) return start;
148  return end; 
149}
150
151G4QParton* G4QString::GetAntiColorParton(void) const
152{
153  G4QParton* start = *(thePartons.begin());
154  G4QParton* end = *(thePartons.end()-1);
155  G4int Encoding = start->GetPDGCode();
156  if (Encoding < -1000 || (Encoding  < 1000 && Encoding > 0)) return end; 
157  return start; 
158}
159
160// Fill parameters
161void G4QString::SetParameters(G4double mCut,G4double clustM, G4double sigQT,G4double DQSup,
162    G4double DQBU, G4double smPar, G4double SSup, G4double SigPt, G4int SLmax, G4int CLmax)
163{//  =============================================================================
164  MassCut         = mCut;           // minimum mass cut for the string
165  ClusterMass     = clustM;         // minimum cluster mass for the fragmentation
166  SigmaQT         = sigQT;          // quark transverse momentum distribution parameter
167  DiquarkSuppress = DQSup;          // is Diquark suppression parameter 
168  DiquarkBreakProb= DQBU;           // is Diquark breaking probability
169  SmoothParam     = smPar;          // QGS model parameter
170  StrangeSuppress = SSup;           // Strangeness suppression parameter
171         widthOfPtSquare = -2*SigPt*SigPt;      // width^2 of pt for string excitation
172  StringLoopInterrupt = SLmax;      // String fragmentation LOOP limit
173  ClusterLoopInterrupt= CLmax;      // Cluster fragmentation LOOP limit
174}
175
176// Pt distribution @@ one can use 1/(1+A*Pt^2)^B
177G4ThreeVector G4QString::GaussianPt(G4double widthSquare, G4double maxPtSquare) const
178{
179         G4double pt2; do{pt2=widthSquare*std::log(G4UniformRand());} while (pt2>maxPtSquare);
180                pt2=std::sqrt(pt2);
181         G4double phi=G4UniformRand()*twopi;
182         return G4ThreeVector(pt2*std::cos(phi),pt2*std::sin(phi),0.);   
183} // End of GaussianPt
184
185// Diffractively excite the string
186void G4QString::DiffString(G4QHadron* hadron, G4bool isProjectile)
187{
188         hadron->SplitUp();
189         G4QParton* start = hadron->GetNextParton();
190         if( start==NULL) 
191         {
192    G4cerr<<"***G4QString::DiffString: No start parton found, nothing is done"<<G4endl;
193           return;
194        }
195         G4QParton* end = hadron->GetNextParton();
196         if( end==NULL) 
197         {
198    G4cerr<<"***G4QString::DiffString: No end parton found, nothing is done"<<G4endl;
199           return;
200         }
201         if(isProjectile)       ExciteString(end, start, 1); //  1 = Projectile
202         else                   ExciteString(start, end,-1); // -1 = Target
203         SetPosition(hadron->GetPosition());
204  // momenta of string ends
205         G4double ptSquared= hadron->Get4Momentum().perp2();
206  G4double hmins=hadron->Get4Momentum().minus();
207  G4double hplus=hadron->Get4Momentum().plus();
208         G4double transMassSquared=hplus*hmins;
209        G4double maxMomentum = std::sqrt(transMassSquared) - std::sqrt(ptSquared);
210         G4double maxAvailMomentumSquared = maxMomentum * maxMomentum;
211        G4ThreeVector pt=GaussianPt(widthOfPtSquare,maxAvailMomentumSquared);
212        G4LorentzVector Pstart(G4LorentzVector(pt,0.));
213         G4LorentzVector Pend(hadron->Get4Momentum().px(), hadron->Get4Momentum().py(), 0.);
214         Pend-=Pstart;
215         G4double tm1=hmins+(Pend.perp2()-Pstart.perp2())/hplus;
216         G4double tm2=std::sqrt( std::max(0., tm1*tm1-4*Pend.perp2()*hmins/hplus ) );
217        G4int Sign= isProjectile ? TARGET : PROJECTILE;
218        G4double endMinus  = 0.5 * (tm1 + Sign*tm2);
219         G4double startMinus= hmins - endMinus;
220        G4double startPlus = Pstart.perp2() / startMinus;
221         G4double endPlus   = hplus - startPlus;
222         Pstart.setPz(0.5*(startPlus - startMinus));
223         Pstart.setE (0.5*(startPlus + startMinus));
224         Pend.setPz  (0.5*(endPlus - endMinus));
225         Pend.setE   (0.5*(endPlus + endMinus));
226         start->Set4Momentum(Pstart);
227         end->Set4Momentum(Pend);
228#ifdef debug
229                G4cout<<"G4QString::DiffString: StartQ="<<start->GetPDGCode()<<start->Get4Momentum()<<"("
230        <<start->Get4Momentum().mag()<<"), EndQ="<<end->GetPDGCode()<<end ->Get4Momentum()
231        <<"("<<end->Get4Momentum().mag()<<"), sumOfEnds="<<Pstart+Pend<<", H4M(original)="
232        <<hadron->Get4Momentum()<<G4endl;
233#endif
234} // End of DiffString (The string is excited diffractively)
235
236// Excite the string by two partons
237void G4QString::ExciteString(G4QParton* Color, G4QParton* AntiColor, G4int Direction)
238{
239                theDirection = Direction;
240  thePosition  = Color->GetPosition();
241  thePartons.push_back(Color);
242  thePartons.push_back(AntiColor);
243} // End of ExciteString
244
245// LUND Longitudinal fragmentation
246G4double G4QString::GetLundLightConeZ(G4double zmin, G4double zmax, G4int , // @@ ? M.K.
247                                      G4QHadron* pHadron, G4double Px, G4double Py)
248{
249    static const G4double  alund = 0.7/GeV/GeV; 
250    // If blund get restored, you MUST adapt the calculation of zOfMaxyf.
251    //static const G4double  blund = 1;
252    G4double z, yf;
253    G4double Mt2 = Px*Px + Py*Py + pHadron->GetMass2();
254    G4double zOfMaxyf=alund*Mt2/(alund*Mt2+1.);
255    G4double maxYf=(1.-zOfMaxyf)/zOfMaxyf * std::exp(-alund*Mt2/zOfMaxyf);
256    do
257    {
258       z = zmin + G4UniformRand()*(zmax-zmin);
259       // yf = std::pow(1. - z, blund)/z*std::exp(-alund*Mt2/z);
260              yf = (1-z)/z * std::exp(-alund*Mt2/z);
261    } while (G4UniformRand()*maxYf>yf); 
262    return z;
263} // End of GetLundLightConeZ
264
265
266// QGSM Longitudinal fragmentation
267G4double G4QString::GetQGSMLightConeZ(G4double zmin, G4double zmax, G4int PartonEncoding,
268                                                                                                                                                                                                                                                                                                                G4QHadron* , G4double, G4double) // @@ M.K.
269{
270  static const G4double arho=0.5;
271  static const G4double aphi=0.;
272  static const G4double an=-0.5;
273  static const G4double ala=-0.75;
274  static const G4double aksi=-1.;
275  static const G4double alft=0.5;
276  G4double z;   
277  G4double theA(0), d2, yf;
278  G4int absCode = std::abs(PartonEncoding);
279  // @@ Crazy algorithm is used for simple power low...
280  if (absCode < 10)                            // Quarks (@@ 9 can be a gluon)
281  { 
282    if(absCode == 1 || absCode == 2) theA = arho;
283    else if(absCode == 3)            theA = aphi;
284    else
285    {
286     G4cerr<<"***G4QString::GetQGSMLightConeZ: CHIPS is SU(3), quakCode="<<absCode<<G4endl;
287     G4Exception("G4QString::GetQGSMLightConeZ:","72",FatalException,"WrongQuark");
288    }
289    do 
290    {
291      z  = zmin + G4UniformRand()*(zmax - zmin);
292      d2 = alft - theA;
293      yf = std::pow( 1.-z, d2);
294    } 
295    while (G4UniformRand()>yf);
296  }
297  else                                     // Di-quarks (@@ Crazy codes are not checked)
298  {       
299    if (absCode==3101||absCode==3103||absCode==3201||absCode==3203)   d2=alft-ala-ala+arho;
300    else if(absCode==1103||absCode==2101||absCode==2203||absCode==2103) d2=alft-an-an+arho;
301    else                                                            d2=alft-aksi-aksi+arho;
302    do 
303    {
304      z = zmin + G4UniformRand()*(zmax - zmin);
305      yf = std::pow(1.-z, d2);
306    } 
307    while (G4UniformRand()>yf); 
308  }
309  return z;
310} // End of GetQGSMLightConeZ
311
312// Diffractively excite the string
313G4QHadronVector* G4QString::FragmentString(G4bool QL)
314{
315  // Can no longer modify Parameters for Fragmentation.
316         // PastInitPhase=true;                                     // Now static
317       
318  //    check if string has enough mass to fragment...
319         G4QHadronVector* LeftVector=LightFragmentationTest();
320         if(LeftVector) return LeftVector;
321       
322         LeftVector = new G4QHadronVector;
323         G4QHadronVector* RightVector = new G4QHadronVector;
324
325  // this should work but it's only a semi deep copy.
326  // %GF        G4ExcitedString theStringInCMS(theString);
327  G4QString* theStringInCMS=CPExcited();                      // must be deleted
328         G4LorentzRotation toCms=theStringInCMS->TransformToAlignedCms();
329         G4bool success=false;
330         G4bool inner_sucess=true;
331         G4int attempt=0;
332         while (!success && attempt++<StringLoopInterrupt)              //@@ It's Loop with break
333         {
334                  G4QString* currentString = new G4QString(*theStringInCMS);
335                  if(LeftVector->size())
336    {
337      std::for_each(LeftVector->begin(), LeftVector->end(), DeleteQHadron());
338                    LeftVector->clear();
339    }
340                  if(RightVector->size())
341                  {
342      std::for_each(RightVector->begin(), RightVector->end(), DeleteQHadron());
343                    RightVector->clear();
344                  }
345                  inner_sucess=true;  // set false on failure..
346                  while (!StopFragmentation())
347                  {  // Split current string into hadron + new string
348                        G4QString* newString=0;  // used as output from SplitUp...
349                           G4QHadron* Hadron=Splitup(QL);
350                           if(Hadron && IsFragmentable())
351                           {
352                             if(currentString->GetDecayDirection()>0) LeftVector->push_back(Hadron);
353        else RightVector->push_back(Hadron);
354                             delete currentString;
355                             currentString=newString;
356                           }
357      else
358      {
359                             // abandon ... start from the beginning
360                             if (newString) delete newString;
361                             if (Hadron)    delete Hadron;
362                             inner_sucess=false;
363                             break;
364                           }
365                  } 
366                // Split current string into 2 final Hadrons
367                  if( inner_sucess && SplitLast(currentString, LeftVector, RightVector) )       success=true;
368                  delete currentString;
369         }
370         delete theStringInCMS;
371         if (!success)
372         {
373                  if(RightVector->size())
374                  {
375      std::for_each(RightVector->begin(), RightVector->end(), DeleteQHadron());
376                    RightVector->clear();
377                  }
378                  delete RightVector;
379                  if(LeftVector->size())
380    {
381      std::for_each(LeftVector->begin(), LeftVector->end(), DeleteQHadron());
382                    LeftVector->clear();
383    }
384                  return LeftVector;
385         }             
386         // Join Left and Right Vectors into LeftVector in correct order.
387         while(!RightVector->empty())
388         {
389            LeftVector->push_back(RightVector->back());
390            RightVector->erase(RightVector->end()-1);
391         }
392         delete RightVector;
393         CalculateHadronTimePosition(Get4Momentum().mag(), LeftVector);
394         G4LorentzRotation toObserverFrame(toCms.inverse());
395         for(unsigned C1 = 0; C1 < LeftVector->size(); C1++)
396         {
397           G4QHadron* Hadron = (*LeftVector)[C1];
398           G4LorentzVector Momentum = Hadron->Get4Momentum();
399           Momentum = toObserverFrame*Momentum;
400           Hadron->Set4Momentum(Momentum);
401           G4LorentzVector Coordinate(Hadron->GetPosition(), Hadron->GetFormationTime());
402           Momentum = toObserverFrame*Coordinate;
403           Hadron->SetFormationTime(Momentum.e());
404           Hadron->SetPosition(GetPosition()+Momentum.vect());
405         }
406         return LeftVector;
407} // End of FragmentLundString
408
409// Creates a string, using only the end-partons of the string
410G4QString* G4QString::CPExcited()
411{
412        G4QParton* LeftParton = new G4QParton(GetLeftParton());
413        G4QParton* RightParton= new G4QParton(GetRightParton());
414        return new G4QString(LeftParton,RightParton,GetDirection());
415} // End of CPExcited
416
417// Simple decay of the string
418G4QHadronVector* G4QString::LightFragmentationTest()
419{
420  // Check string decay threshold
421               
422         G4QHadronVector* result=0;  // return 0 when string exceeds the mass cut
423       
424         G4QHadronPair hadrons((G4QHadron*)0, (G4QHadron*)0); // pair of hadrons
425         if(sqr(FragmentationMass(0,&hadrons)+MassCut)<Mass2()) return 0; //Par MassCut
426       
427         result = new G4QHadronVector;
428       
429         if(hadrons.second==0)                   // Second hadron exists
430         {
431           // Substitute string by a light hadron, Note that Energy is not conserved here! @@
432           G4ThreeVector Mom3 = Get4Momentum().vect();
433           G4LorentzVector Mom(Mom3, std::sqrt(Mom3.mag2() + hadrons.first->GetMass2()) );
434    result->push_back(new G4QHadron(hadrons.first, 0, GetPosition(), Mom));
435         }
436  else 
437         {
438           //... string was qq--qqbar type: Build two stable hadrons,
439           G4LorentzVector  Mom1, Mom2;
440           Sample4Momentum(&Mom1, hadrons.first->GetMass(), &Mom2, hadrons.second->GetMass(),
441                                                                                  Get4Momentum().mag());
442           result->push_back(new G4QHadron(hadrons.first, 0, GetPosition(), Mom1));
443           result->push_back(new G4QHadron(hadrons.second,0, GetPosition(), Mom2));
444    G4ThreeVector Velocity = Get4Momentum().boostVector();
445    G4int L=result->size(); if(L) for(G4int i=0; i<L; i++) (*result)[i]->Boost(Velocity);
446         }
447         return result;
448} // End of LightFragmentationTest
449
450// Calculate Fragmentation Mass (if pdefs # 0, returns two hadrons)
451G4double G4QString::FragmentationMass(G4QHcreate build, G4QHadronPair* pdefs)
452{
453  G4double mass;
454  // Example how to use an interface to different member functions
455         if(build==0) build=&G4QHadronBuilder::BuildLowSpin; // @@ Build S Hadrons?
456  G4QHadron* Hadron1 = 0;                            // @@ Not initialized
457  G4QHadron* Hadron2 = 0;
458  if(!FourQuarkString() )
459  {
460    // spin 0 meson or spin 1/2 barion will be built
461    Hadron1 = (hadronizer->*build)(GetLeftParton(), GetRightParton());
462    mass    = Hadron1->GetMass();
463  }
464  else
465  {
466    // string is qq--qqbar: Build two stable hadrons, with extra uubar or ddbar quark pair
467           G4int iflc = (G4UniformRand() < 0.5)? 1 : 2;
468           if (GetLeftParton()->GetPDGCode() < 0) iflc = -iflc;
469           //... theSpin = 4; spin 3/2 baryons will be built
470           Hadron1 = (hadronizer->*build)(GetLeftParton(),CreateParton(iflc));
471           Hadron2 = (hadronizer->*build)(GetRightParton(),CreateParton(-iflc));
472    mass    = Hadron1->GetMass() + Hadron2->GetMass();
473  }
474         if(pdefs) // need to return hadrons as well....
475         {
476           pdefs->first  = Hadron1;
477           pdefs->second = Hadron2;
478         } 
479  return mass;
480} // End of FragmentationMass
481
482 // Checks that the string is qq-(qq)bar string
483G4bool G4QString::FourQuarkString() const
484{
485  return    GetLeftParton()->GetParticleSubType() == "di_quark"
486         && GetRightParton()->GetParticleSubType()== "di_quark";
487} // End of FourQuarkString
488
489void G4QString::CalculateHadronTimePosition(G4double StringMass, G4QHadronVector* Hadrons)
490{
491  // `yo-yo` formation time
492  static const G4double kappa = 1.0 * GeV/fermi;
493  static const G4double dkappa = kappa+kappa;
494  for(unsigned c1 = 0; c1 < Hadrons->size(); c1++)
495  {
496    G4double SumPz = 0; 
497    G4double SumE  = 0;
498    for(unsigned c2 = 0; c2 < c1; c2++)
499    {
500      G4LorentzVector hc2M=(*Hadrons)[c2]->Get4Momentum();
501      SumPz += hc2M.pz();
502      SumE  += hc2M.e();   
503    }
504    G4QHadron* hc1=(*Hadrons)[c1];
505    G4LorentzVector hc1M=hc1->Get4Momentum();
506    G4double HadronE = hc1M.e();
507    G4double HadronPz= hc1M.pz();
508    hc1->SetFormationTime((StringMass-SumPz-SumPz+HadronE-HadronPz)/dkappa);
509    hc1->SetPosition(G4ThreeVector(0,0,(StringMass-SumE-SumE-HadronE+HadronPz)/dkappa));
510  } 
511} // End of CalculateHadronTimePosition
512
513void G4QString::SetLeftPartonStable()
514{
515     theStableParton=GetLeftParton();
516     theDecayParton=GetRightParton();
517     decaying=Right;
518}
519
520void G4QString::SetRightPartonStable()
521{
522     theStableParton=GetRightParton();
523     theDecayParton=GetLeftParton();
524     decaying=Left;
525}
526
527G4int G4QString::GetDecayDirection() const
528{
529         if      (decaying == Left ) return +1;
530         else if (decaying == Right) return -1;
531         else
532         {
533                                G4cerr<<"***G4QString::GetDecayDirection: wrong DecayDirection="<<decaying<<G4endl;
534    G4Exception("G4QString::GetDecayDirection:","72",FatalException,"WrongDecayDirection");
535         }
536         return 0;
537}
538
539G4ThreeVector G4QString::StablePt()
540{
541         if (decaying == Left ) return Ptright;
542         else if (decaying == Right ) return Ptleft;
543         else
544         {
545                                G4cerr<<"***G4QString::StablePt: wrong DecayDirection="<<decaying<<G4endl;
546    G4Exception("G4QString::StablePt:","72",FatalException,"WrongDecayDirection");
547         }
548         return G4ThreeVector();
549}
550
551G4ThreeVector G4QString::DecayPt()
552{
553         if      (decaying == Left  ) return Ptleft;
554         else if (decaying == Right ) return Ptright;
555         else
556         {
557                                G4cerr<<"***G4QString::DecayPt: wrong DecayDirection="<<decaying<<G4endl;
558    G4Exception("G4QString::DecayPt:","72",FatalException,"WrongDecayDirection");
559         }
560         return G4ThreeVector();
561}
562
563G4double G4QString::LightConeDecay()
564{
565         if      (decaying == Left  ) return Pplus;
566         else if (decaying == Right ) return Pminus;
567         else
568         {
569                                G4cerr<<"***G4QString::LightConeDecay: wrong DecayDirection="<<decaying<<G4endl;
570    G4Exception("G4QString::LightConeDecay:","72",FatalException,"WrongDecayDirection");
571         }
572         return 0;
573}
574
575G4LorentzVector G4QString::GetFragmentation4Mom() const
576{
577  G4LorentzVector momentum(Ptleft+Ptright,0.5*(Pplus+Pminus));
578         momentum.setPz(0.5*(Pplus-Pminus));
579         return momentum;
580}
581
582// Random choice of string end to use it for creating the hadron (decay)   
583G4QHadron* G4QString::Splitup(G4bool QL)
584{
585  SideOfDecay = (G4UniformRand() < 0.5)? 1: -1;
586  if(SideOfDecay<0) SetLeftPartonStable();  // Decay Right parton
587  else              SetRightPartonStable(); // Decay Left parton
588  G4QParton* newStringEnd;
589  G4QHadron* Hadron;
590  if(DecayIsQuark()) Hadron=QuarkSplitup(GetDecayParton(), newStringEnd);    // MF1
591  else               Hadron= DiQuarkSplitup(GetDecayParton(), newStringEnd); // MF2
592  // create new String from old, ie. keep Left and Right order, but replace decay
593  G4LorentzVector* HadronMomentum=SplitEandP(Hadron, QL);                    // MF3
594  if(HadronMomentum)
595  {   
596           G4ThreeVector   Pos(0.,0.,0.);
597           Hadron->Set4Momentum(*HadronMomentum);
598           UpdateString(newStringEnd, HadronMomentum);
599           delete HadronMomentum;
600  }     
601  return Hadron;
602} // End of Splitup
603
604void G4QString::UpdateString(G4QParton* decay, const G4LorentzVector* mom)
605{
606         decaying=None;
607         if(decaying == Left)
608         {
609    G4QParton* theFirst = thePartons[0];
610                                delete theFirst;
611                  theFirst = decay;
612                  Ptleft  -= mom->vect();
613                  Ptleft.setZ(0.);
614         }
615  else if (decaying == Right)
616         {
617    G4QParton* theLast = thePartons[thePartons.size()-1];
618    delete theLast;
619                  theLast  = decay;
620                  Ptright -= mom->vect();
621                  Ptright.setZ(0.);
622         }
623  else
624         {
625                                G4cerr<<"***G4QString::UpdateString: wrong oldDecay="<<decaying<<G4endl;
626    G4Exception("G4QString::UpdateString","72",FatalException,"WrongDecayDirection");
627         }
628         Pplus  -= mom->e() + mom->pz();
629         Pminus -= mom->e() - mom->pz();
630} // End of UpdateString
631
632// QL=true for QGSM and QL=false for Lund fragmentation
633G4LorentzVector* G4QString::SplitEandP(G4QHadron* pHadron, G4bool QL)
634{
635  G4double HadronMass = pHadron->GetMass();
636  // calculate and assign hadron transverse momentum component HadronPx andHadronPy
637  G4ThreeVector HadronPt = SampleQuarkPt() + DecayPt();
638  HadronPt.setZ(0.);
639  //...  sample z to define hadron longitudinal momentum and energy
640  //... but first check the available phase space
641  G4double DecayQuarkMass2  = sqr(GetDecayParton()->GetMass()); // Mass of quark? M.K.
642  G4double HadronMass2T = HadronMass*HadronMass + HadronPt.mag2();
643  if (DecayQuarkMass2 + HadronMass2T >= SmoothParam*Mass2() )  return 0;                // restart!
644  //... then compute allowed z region  z_min <= z <= z_max
645  G4double zMin = HadronMass2T/Mass2();
646  G4double zMax = 1. - DecayQuarkMass2/Mass2();
647  if (zMin >= zMax) return 0;           // have to start all over! 
648  G4double z=0;
649  if(QL) z = GetQGSMLightConeZ(zMin, zMax, GetDecayParton()->GetPDGCode(), pHadron,
650                               HadronPt.x(), HadronPt.y());     
651  else   z = GetLundLightConeZ(zMin, zMax, GetDecayParton()->GetPDGCode(), pHadron,
652                               HadronPt.x(), HadronPt.y());     
653  //... now compute hadron longitudinal momentum and energy
654  // longitudinal hadron momentum component HadronPz
655  G4double zl= z*LightConeDecay();
656  G4double HadronE = (zl + HadronMass2T/zl)/2;
657  HadronPt.setZ( GetDecayDirection() * (zl - HadronE) );
658  G4LorentzVector* a4Momentum= new G4LorentzVector(HadronPt,HadronE);
659
660  return a4Momentum;
661}
662
663G4ThreeVector G4QString::SampleQuarkPt()
664{
665   G4double Pt = SigmaQT * std::sqrt( -std::log(G4UniformRand()));
666   G4double phi = twopi*G4UniformRand();
667   return G4ThreeVector(Pt * std::cos(phi),Pt * std::sin(phi),0);
668}
669
670void G4QString::Sample4Momentum(G4LorentzVector* Mom, G4double Mass,
671                         G4LorentzVector* AntiMom, G4double AntiMass, G4double InitialMass)
672{
673  G4double m2 = Mass*Mass;
674  G4double am2= AntiMass*AntiMass;
675  G4double dub=InitialMass*InitialMass - m2 - am2;
676  G4double r_val = dub - 4*m2*am2;
677  G4double Pabs = (r_val > 0.) ? std::sqrt(r_val)/(InitialMass*InitialMass) : 0;
678  //... sample unit vector       
679  G4double r  = G4UniformRand();                    // @@ G4RandomDirection()
680  G4double pz = 1. - r - r;                         // cos(theta)
681  G4double st = std::sqrt(1. - pz * pz) * Pabs;
682  G4double phi= twopi*G4UniformRand();
683  G4double px = st*std::cos(phi);
684  G4double py = st*std::sin(phi);
685  pz *= Pabs;
686  G4double p2=Pabs*Pabs;
687  Mom->setPx(px); Mom->setPy(py); Mom->setPz(pz);
688  Mom->setE(std::sqrt(p2 + Mass*Mass));
689  AntiMom->setPx(-px); AntiMom->setPy(-py); AntiMom->setPz(-pz);
690  AntiMom->setE (std::sqrt(Pabs*Pabs + AntiMass*AntiMass));
691} // End of Sample4Momentum
692
693G4bool G4QString::SplitLast(G4QString* string, G4QHadronVector* LeftVector,
694                                               G4QHadronVector* RightVector)
695{
696  //... perform last cluster decay
697  G4ThreeVector ClusterVel =string->Get4Momentum().boostVector();
698  G4double ResidualMass    =string->Mass(); 
699  G4double ClusterMassCut = ClusterMass;
700  G4int cClusterInterrupt = 0;
701  G4QHadron* LeftHadron;
702  G4QHadron* RightHadron;
703  do
704  {
705    if(cClusterInterrupt++ >= ClusterLoopInterrupt) return false;
706           G4QParton* quark = 0;
707           string->SetLeftPartonStable(); // to query quark contents..
708           if(string->DecayIsQuark() && string->StableIsQuark()) // There're quarks on clusterEnds
709                    LeftHadron= QuarkSplitup(string->GetLeftParton(), quark);
710           else
711    {
712             //... there is a Diquark on cluster ends
713                    G4int IsParticle;
714                    if(string->StableIsQuark())IsParticle=(string->GetLeftParton()->GetPDGCode()>0)?-1:1;
715                    else                       IsParticle=(string->GetLeftParton()->GetPDGCode()>0)?1:-1;
716      G4QPartonPair QuarkPair = CreatePartonPair(IsParticle,false);  // no diquarks wanted
717      quark = QuarkPair.GetParton2();
718      LeftHadron=hadronizer->Build(QuarkPair.GetParton1(), string->GetLeftParton());
719           }
720    RightHadron = hadronizer->Build(string->GetRightParton(), quark);
721    //... repeat procedure, if mass of cluster is too low to produce hadrons
722    //... ClusterMassCut = 0.15*GeV model parameter
723           if ( quark->GetParticleSubType()== "quark" ) ClusterMassCut = 0.;
724           else                                         ClusterMassCut = ClusterMass;
725  } while(ResidualMass <= LeftHadron->GetMass() + RightHadron->GetMass() + ClusterMassCut);
726  //... compute hadron momenta and energies   
727  G4LorentzVector  LeftMom, RightMom;
728  G4ThreeVector    Pos;
729  Sample4Momentum(&LeftMom,LeftHadron->GetMass(),&RightMom,RightHadron->GetMass(),
730                                                                             ResidualMass);
731  LeftMom.boost(ClusterVel);
732  RightMom.boost(ClusterVel);
733  LeftVector->push_back(new G4QHadron(LeftHadron, 0, Pos, LeftMom));
734  RightVector->push_back(new G4QHadron(RightHadron, 0, Pos, RightMom));
735
736  return true;
737} // End of SplitLast
738
739G4QHadron*      G4QString::QuarkSplitup(G4QParton*      decay, G4QParton *&created)
740{
741  G4int IsParticle=(decay->GetPDGCode()>0) ? -1 : +1; // if we have a quark, we need antiquark (or diquark)
742  G4QPartonPair QuarkPair = CreatePartonPair(IsParticle);
743  created = QuarkPair.GetParton2();
744  return hadronizer->Build(QuarkPair.GetParton1(), decay);
745} // End of QuarkSplitup
746
747//
748G4QHadron* G4QString::DiQuarkSplitup(G4QParton* decay, G4QParton* &created)
749{
750  //... can Diquark break or not?
751  if (G4UniformRand() < DiquarkBreakProb )
752  {
753    //... Diquark break
754    G4int stableQuarkEncoding = decay->GetPDGCode()/1000;
755    G4int decayQuarkEncoding = (decay->GetPDGCode()/100)%10;
756    if (G4UniformRand() < 0.5)
757    {
758      G4int Swap = stableQuarkEncoding;
759      stableQuarkEncoding = decayQuarkEncoding;
760      decayQuarkEncoding = Swap;
761    }
762    G4int IsParticle=(decayQuarkEncoding>0) ? -1 : +1; 
763           // if we have a quark, we need antiquark)
764    G4QPartonPair QuarkPair = CreatePartonPair(IsParticle,false);  // no diquarks wanted
765    //... Build new Diquark
766    G4int QuarkEncoding=QuarkPair.GetParton2()->GetPDGCode();
767    G4int i10  = std::max(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
768    G4int i20  = std::min(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
769    G4int spin = (i10 != i20 && G4UniformRand() <= 0.5)? 1 : 3;
770    G4int NewDecayEncoding = -1*IsParticle*(i10 * 1000 + i20 * 100 + spin);
771    created = CreateParton(NewDecayEncoding);
772    G4QParton* decayQuark=CreateParton(decayQuarkEncoding);
773    return hadronizer->Build(QuarkPair.GetParton1(), decayQuark);
774  }
775  else
776  {
777    //... Diquark does not break
778    G4int IsParticle=(decay->GetPDGCode()>0) ? +1 : -1; 
779        // if we have a diquark, we need quark)
780    G4QPartonPair QuarkPair = CreatePartonPair(IsParticle,false);  // no diquarks wanted
781    created = QuarkPair.GetParton2();
782    return hadronizer->Build(QuarkPair.GetParton1(), decay);
783 }
784} // End of DiQuarkSplitup
785
786G4QPartonPair G4QString::CreatePartonPair(G4int NeedParticle, G4bool AllowDiquarks)
787{
788  //  NeedParticle = {+1 for Particle, -1 for AntiParticle}
789  if(AllowDiquarks && G4UniformRand() < DiquarkSuppress)
790  {
791    // Create a Diquark - AntiDiquark pair , first in pair is anti to IsParticle
792    G4int q1  = SampleQuarkFlavor();
793    G4int q2  = SampleQuarkFlavor();
794    G4int spin = (q1 != q2 && G4UniformRand() <= 0.5)? 1 : 3;
795    // Convention: quark with higher PDG number is first
796    G4int PDGcode = (std::max(q1,q2) * 1000 + std::min(q1,q2) * 100 + spin) * NeedParticle;
797    return G4QPartonPair(CreateParton(-PDGcode),CreateParton(PDGcode));
798  }
799  else
800  {
801    // Create a Quark - AntiQuark pair, first in pair  IsParticle
802    G4int PDGcode=SampleQuarkFlavor()*NeedParticle;
803    return G4QPartonPair(CreateParton(PDGcode),CreateParton(-PDGcode));
804  }
805} // End of CreatePartonPair
Note: See TracBrowser for help on using the repository browser.