source: trunk/source/processes/hadronic/models/parton_string/hadronization/src/G4QGSMFragmentation.cc @ 1340

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

update ti head
File size: 13.7 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: G4QGSMFragmentation.cc,v 1.10 2010/09/20 12:46:23 vuzhinsk Exp $
28// GEANT4 tag $Name: geant4-09-03-ref-09 $
29//
30// -----------------------------------------------------------------------------
31//      GEANT 4 class implementation file
32//
33//      History: first implementation, Maxim Komogorov, 10-Jul-1998
34// -----------------------------------------------------------------------------
35#include "G4QGSMFragmentation.hh"
36#include "G4FragmentingString.hh"
37#include "G4DiQuarks.hh"
38#include "G4Quarks.hh"
39
40#include "Randomize.hh"
41#include "G4ios.hh"
42
43// Class G4QGSMFragmentation
44//****************************************************************************************
45 
46G4QGSMFragmentation::G4QGSMFragmentation() :
47arho(0.5), aphi(0.), an(-0.5), ala(-0.75), aksi(-1.), alft(0.5)
48   {
49   }
50
51G4QGSMFragmentation::G4QGSMFragmentation(const G4QGSMFragmentation &) : G4VLongitudinalStringDecay(),
52arho(0.5), aphi(0.), an(-0.5), ala(-0.75), aksi(-1.), alft(0.5)
53   {
54   }
55
56G4QGSMFragmentation::~G4QGSMFragmentation()
57   {
58   }
59
60//****************************************************************************************
61
62const G4QGSMFragmentation & G4QGSMFragmentation::operator=(const G4QGSMFragmentation &)
63   {
64    throw G4HadronicException(__FILE__, __LINE__, "G4QGSMFragmentation::operator= meant to not be accessable");
65    return *this;
66   }
67
68int G4QGSMFragmentation::operator==(const G4QGSMFragmentation &right) const
69   {
70   return !memcmp(this, &right, sizeof(G4QGSMFragmentation));
71   }
72
73int G4QGSMFragmentation::operator!=(const G4QGSMFragmentation &right) const
74   {
75   return memcmp(this, &right, sizeof(G4QGSMFragmentation));
76   }
77 
78//****************************************************************************************
79//----------------------------------------------------------------------------------------------------------
80
81G4KineticTrackVector* G4QGSMFragmentation::FragmentString(const G4ExcitedString& theString)
82{
83//    Can no longer modify Parameters for Fragmentation.
84        PastInitPhase=true;
85       
86//      check if string has enough mass to fragment...
87        G4KineticTrackVector * LeftVector=LightFragmentationTest(&theString);
88        if ( LeftVector != 0 ) return LeftVector;
89       
90        LeftVector = new G4KineticTrackVector;
91        G4KineticTrackVector * RightVector=new G4KineticTrackVector;
92
93// this should work but its only a semi deep copy. %GF  G4ExcitedString theStringInCMS(theString);
94        G4ExcitedString *theStringInCMS=CPExcited(theString);
95        G4LorentzRotation toCms=theStringInCMS->TransformToAlignedCms();
96
97        G4bool success=false, inner_sucess=true;
98        G4int attempt=0;
99        while ( !success && attempt++ < StringLoopInterrupt )
100        {
101                G4FragmentingString *currentString=new G4FragmentingString(*theStringInCMS);
102
103                std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
104                LeftVector->clear();
105                std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
106                RightVector->clear();
107               
108                inner_sucess=true;  // set false on failure..
109                while (! StopFragmenting(currentString) )
110                {  // Split current string into hadron + new string
111                        G4FragmentingString *newString=0;  // used as output from SplitUp...
112                        G4KineticTrack * Hadron=Splitup(currentString,newString);
113                        if ( Hadron != 0 && IsFragmentable(newString)) 
114                        {
115                           if ( currentString->GetDecayDirection() > 0 )
116                                   LeftVector->push_back(Hadron);
117                           else
118                                   RightVector->push_back(Hadron);
119                           delete currentString;
120                           currentString=newString;
121                        } else {
122                         // abandon ... start from the beginning
123                           if (newString) delete newString;          // Uzhi restore 20.06.08
124                           if (Hadron)    delete Hadron;
125                           inner_sucess=false;
126                           break;
127                        }
128                } 
129                // Split current string into 2 final Hadrons
130                if ( inner_sucess && 
131                     SplitLast(currentString,LeftVector, RightVector) ) 
132                {
133                        success=true;
134                }
135                delete currentString;
136        }
137       
138        delete theStringInCMS;
139       
140        if ( ! success )
141        {
142                std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
143                LeftVector->clear();
144                std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
145                delete RightVector;
146                return LeftVector;
147        }
148               
149        // Join Left- and RightVector into LeftVector in correct order.
150        while(!RightVector->empty())
151        {
152            LeftVector->push_back(RightVector->back());
153            RightVector->erase(RightVector->end()-1);
154        }
155        delete RightVector;
156
157        CalculateHadronTimePosition(theString.Get4Momentum().mag(), LeftVector);
158
159        G4LorentzRotation toObserverFrame(toCms.inverse());
160
161        for(size_t C1 = 0; C1 < LeftVector->size(); C1++)
162        {
163           G4KineticTrack* Hadron = LeftVector->operator[](C1);
164           G4LorentzVector Momentum = Hadron->Get4Momentum();
165           Momentum = toObserverFrame*Momentum;
166           Hadron->Set4Momentum(Momentum);
167           G4LorentzVector Coordinate(Hadron->GetPosition(), Hadron->GetFormationTime());
168           Momentum = toObserverFrame*Coordinate;
169           Hadron->SetFormationTime(Momentum.e());
170           G4ThreeVector aPosition(Momentum.vect());
171           Hadron->SetPosition(theString.GetPosition()+aPosition);
172        }
173        return LeftVector;
174               
175
176
177}
178
179//----------------------------------------------------------------------------------------------------------
180
181G4double G4QGSMFragmentation::GetLightConeZ(G4double zmin, G4double zmax, G4int PartonEncoding,  G4ParticleDefinition* , G4double , G4double )
182{   
183  G4double z;   
184  G4double theA(0), d1, d2, yf;
185  G4int absCode = std::abs( PartonEncoding );
186  if (absCode < 10)
187  { 
188    if(absCode == 1 || absCode == 2) theA = arho;
189    else if(absCode == 3) theA = aphi;
190    else throw G4HadronicException(__FILE__, __LINE__, "Unknown PDGencoding in G4QGSMFragmentation::G4LightConeZ");
191
192    do 
193    {
194      z  = zmin + G4UniformRand() * (zmax - zmin);
195      d1 =  (1. - z);
196      d2 =  (alft - theA);
197      yf = std::pow(d1, d2);
198    } 
199    while (G4UniformRand() > yf);
200  }
201  else
202  {       
203    if(absCode == 1103 || absCode == 2101 || 
204       absCode == 2203 || absCode == 2103)
205    {
206      d2 =  (alft - (2.*an - arho));
207    }
208    else if(absCode == 3101 || absCode == 3103 ||
209            absCode == 3201 || absCode == 3203)
210    {
211      d2 =  (alft - (2.*ala - arho));
212    }
213    else
214    {
215      d2 =  (alft - (2.*aksi - arho));
216    }
217
218    do 
219    {
220      z = zmin + G4UniformRand() * (zmax - zmin);
221      d1 =  (1. - z);
222      yf = std::pow(d1, d2);
223    } 
224    while (G4UniformRand() > yf); 
225  }
226  return z;
227}
228//-----------------------------------------------------------------------------------------
229
230G4LorentzVector * G4QGSMFragmentation::SplitEandP(G4ParticleDefinition * pHadron,
231                                                  G4FragmentingString * string,    // Uzhi
232                                                  G4FragmentingString * ) // Uzhi
233{
234       G4double HadronMass = pHadron->GetPDGMass();
235
236       // calculate and assign hadron transverse momentum component HadronPx andHadronPy
237       G4ThreeVector thePt;
238       thePt=SampleQuarkPt();
239       G4ThreeVector HadronPt = thePt +string->DecayPt();
240       HadronPt.setZ(0);
241       //...  sample z to define hadron longitudinal momentum and energy
242       //... but first check the available phase space
243       G4double DecayQuarkMass2  = sqr(string->GetDecayParton()->GetPDGMass());
244       G4double HadronMass2T = sqr(HadronMass) + HadronPt.mag2();
245       if (DecayQuarkMass2 + HadronMass2T >= SmoothParam*(string->Mass2()) ) 
246          return 0;             // have to start all over!
247
248       //... then compute allowed z region  z_min <= z <= z_max
249 
250       G4double zMin = HadronMass2T/(string->Mass2());
251       G4double zMax = 1. - DecayQuarkMass2/(string->Mass2());
252       if (zMin >= zMax) return 0;              // have to start all over!
253       
254       G4double z = GetLightConeZ(zMin, zMax,
255                       string->GetDecayParton()->GetPDGEncoding(), pHadron,
256                       HadronPt.x(), HadronPt.y());     
257       
258       //... now compute hadron longitudinal momentum and energy
259       // longitudinal hadron momentum component HadronPz
260
261        HadronPt.setZ(0.5* string->GetDecayDirection() *
262                        (z * string->LightConeDecay() - 
263                         HadronMass2T/(z * string->LightConeDecay())));
264        G4double HadronE  = 0.5* (z * string->LightConeDecay() + 
265                                  HadronMass2T/(z * string->LightConeDecay()));
266
267       G4LorentzVector * a4Momentum= new G4LorentzVector(HadronPt,HadronE);
268
269       return a4Momentum;
270}
271
272
273//-----------------------------------------------------------------------------------------
274
275G4bool G4QGSMFragmentation::SplitLast(G4FragmentingString * string,
276                                             G4KineticTrackVector * LeftVector,
277                                             G4KineticTrackVector * RightVector)
278{
279    //... perform last cluster decay
280    G4ThreeVector ClusterVel =string->Get4Momentum().boostVector();
281    G4double ResidualMass    =string->Mass(); 
282    G4double ClusterMassCut = ClusterMass;
283    G4int cClusterInterrupt = 0;
284    G4ParticleDefinition * LeftHadron, * RightHadron;
285    do
286    {
287        if (cClusterInterrupt++ >= ClusterLoopInterrupt)
288        {
289          return false;
290        }
291        G4ParticleDefinition * quark = NULL;
292        string->SetLeftPartonStable(); // to query quark contents..
293        if (string->DecayIsQuark() && string->StableIsQuark() ) 
294        {
295           //... there are quarks on cluster ends
296                LeftHadron= QuarkSplitup(string->GetLeftParton(), quark);
297        } else {
298           //... there is a Diquark on cluster ends
299                G4int IsParticle;
300                if ( string->StableIsQuark() ) {
301                  IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? -1 : +1; 
302                } else {
303                  IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? +1 : -1;
304                }
305                pDefPair QuarkPair = CreatePartonPair(IsParticle,false);  // no diquarks wanted
306                quark = QuarkPair.second;
307                LeftHadron=hadronizer->Build(QuarkPair.first, string->GetLeftParton());
308        }
309        RightHadron = hadronizer->Build(string->GetRightParton(), quark);
310
311       //... repeat procedure, if mass of cluster is too low to produce hadrons
312       //... ClusterMassCut = 0.15*GeV model parameter
313        if ( quark->GetParticleSubType()== "quark" ) {ClusterMassCut = 0.;}
314        else {ClusterMassCut = ClusterMass;}
315    } 
316    while (ResidualMass <= LeftHadron->GetPDGMass() + RightHadron->GetPDGMass()  + ClusterMassCut);
317
318    //... compute hadron momenta and energies   
319    G4LorentzVector  LeftMom, RightMom;
320    G4ThreeVector    Pos;
321    Sample4Momentum(&LeftMom, LeftHadron->GetPDGMass(), &RightMom, RightHadron->GetPDGMass(), ResidualMass);
322    LeftMom.boost(ClusterVel);
323    RightMom.boost(ClusterVel);
324    LeftVector->push_back(new G4KineticTrack(LeftHadron, 0, Pos, LeftMom));
325    RightVector->push_back(new G4KineticTrack(RightHadron, 0, Pos, RightMom));
326
327    return true;
328
329}
330
331//----------------------------------------------------------------------------------------------------------
332
333G4bool G4QGSMFragmentation::IsFragmentable(const G4FragmentingString * const string)
334{
335        return sqr(FragmentationMass(string)+MassCut) <
336                        string->Mass2();
337}
338
339//----------------------------------------------------------------------------------------------------------
340
341G4bool G4QGSMFragmentation::StopFragmenting(const G4FragmentingString * const string)
342{
343        return
344         sqr(FragmentationMass(string,&G4HadronBuilder::BuildHighSpin)+MassCut) >
345         string->Get4Momentum().mag2();
346}
347
348//----------------------------------------------------------------------------------------------------------
349
350//----------------------------------------------------------------------------------------------------------
351
352void G4QGSMFragmentation::Sample4Momentum(G4LorentzVector* Mom, G4double Mass, G4LorentzVector* AntiMom, G4double AntiMass, G4double InitialMass) 
353    {
354    G4double r_val = sqr(InitialMass*InitialMass - Mass*Mass - AntiMass*AntiMass) - sqr(2.*Mass*AntiMass);
355    G4double Pabs = (r_val > 0.)? std::sqrt(r_val)/(2.*InitialMass) : 0;
356
357    //... sample unit vector       
358    G4double pz = 1. - 2.*G4UniformRand(); 
359    G4double st     = std::sqrt(1. - pz * pz)*Pabs;
360    G4double phi    = 2.*pi*G4UniformRand();
361    G4double px = st*std::cos(phi);
362    G4double py = st*std::sin(phi);
363    pz *= Pabs;
364   
365    Mom->setPx(px); Mom->setPy(py); Mom->setPz(pz);
366    Mom->setE(std::sqrt(Pabs*Pabs + Mass*Mass));
367
368    AntiMom->setPx(-px); AntiMom->setPy(-py); AntiMom->setPz(-pz);
369    AntiMom->setE (std::sqrt(Pabs*Pabs + AntiMass*AntiMass));
370    }
371
372   
373//*********************************************************************************************
Note: See TracBrowser for help on using the repository browser.