source: trunk/source/processes/hadronic/models/parton_string/hadronization/src/G4HadronBuilder.cc

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

update ti head
File size: 9.4 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: G4HadronBuilder.cc,v 1.11 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:
34//             Gunter Folger, August/September 2001
35//               Create class; algorithm previously in G4VLongitudinalStringDecay.
36// -----------------------------------------------------------------------------
37
38#include "G4HadronBuilder.hh"
39#include "G4HadronicException.hh"
40#include "Randomize.hh"
41#include "G4ParticleTable.hh"
42
43G4HadronBuilder::G4HadronBuilder(G4double mesonMix, G4double barionMix,
44                     std::vector<double> scalarMesonMix,
45                     std::vector<double> vectorMesonMix)
46{
47        mesonSpinMix=mesonMix;       
48        barionSpinMix=barionMix;
49        scalarMesonMixings=scalarMesonMix;
50        vectorMesonMixings=vectorMesonMix;
51}
52
53G4ParticleDefinition * G4HadronBuilder::Build(G4ParticleDefinition * black, G4ParticleDefinition * white)
54{
55
56        if (black->GetParticleSubType()== "di_quark" || white->GetParticleSubType()== "di_quark" ) {
57
58//    Barion
59           Spin spin = (G4UniformRand() < barionSpinMix) ? SpinHalf : SpinThreeHalf;
60           return Barion(black,white,spin);
61
62        } else {       
63
64//    Meson
65           Spin spin = (G4UniformRand() < mesonSpinMix) ? SpinZero : SpinOne;
66           return Meson(black,white,spin);
67
68        }
69}
70
71//-------------------------------------------------------------------------
72
73G4ParticleDefinition * G4HadronBuilder::BuildLowSpin(G4ParticleDefinition * black, G4ParticleDefinition * white)
74{
75        if ( black->GetParticleSubType()== "quark" && white->GetParticleSubType()== "quark" ) {
76                return Meson(black,white, SpinZero);
77        } else {
78//                    will return a SpinThreeHalf Barion if all quarks the same
79                return Barion(black,white, SpinHalf); 
80        }       
81}
82
83//-------------------------------------------------------------------------
84
85G4ParticleDefinition * G4HadronBuilder::BuildHighSpin(G4ParticleDefinition * black, G4ParticleDefinition * white)
86{
87        if ( black->GetParticleSubType()== "quark" && white->GetParticleSubType()== "quark" ) {
88                return Meson(black,white, SpinOne);
89        } else {
90                return Barion(black,white,SpinThreeHalf);
91        }
92}
93
94//-------------------------------------------------------------------------
95
96G4ParticleDefinition * G4HadronBuilder::Meson(G4ParticleDefinition * black, 
97                                              G4ParticleDefinition * white, Spin theSpin)
98{
99#ifdef G4VERBOSE
100//  Verify Input Charge
101   
102   G4double charge =  black->GetPDGCharge() 
103                    + white->GetPDGCharge();     
104   if (std::abs(charge) > 2 || std::abs(3.*charge - 3*G4int(charge*1.001)) > perCent )   // 1.001 to avoid int(.9999) -> 0
105        {
106            G4cerr << " G4HadronBuilder::Build()" << G4endl;
107            G4cerr << "    Invalid total charge found for on input: " 
108                        << charge<< G4endl;
109            G4cerr << "    PGDcode input quark1/quark2 : " <<
110                        black->GetPDGEncoding() << " / "<< 
111                        white->GetPDGEncoding() << G4endl;
112            G4cerr << G4endl;
113        } 
114#endif 
115       
116        G4int id1= black->GetPDGEncoding();
117        G4int id2= white->GetPDGEncoding();
118//      G4int ifl1= std::max(std::abs(id1), std::abs(id2));
119        if ( std::abs(id1) < std::abs(id2) )
120           {
121           G4int xchg = id1; 
122           id1 = id2; 
123           id2 = xchg;
124           }
125       
126        if (std::abs(id1) > 3 ) 
127           throw G4HadronicException(__FILE__, __LINE__, "G4HadronBuilder::Meson : Illegal Quark content as input");
128       
129        G4int PDGEncoding=0;
130
131        if (id1 + id2 == 0) {     
132           G4double rmix = G4UniformRand();
133           G4int    imix = 2*std::abs(id1) - 1;
134           if(theSpin == SpinZero) {
135              PDGEncoding = 110*(1 + (G4int)(rmix + scalarMesonMixings[imix - 1])
136                                   + (G4int)(rmix + scalarMesonMixings[imix])
137                                ) +  theSpin;
138           } else {
139              PDGEncoding = 110*(1 + (G4int)(rmix + vectorMesonMixings[imix - 1])
140                                   + (G4int)(rmix + vectorMesonMixings[imix])
141                                ) +  theSpin;
142           }
143        } else {
144           PDGEncoding = 100 * std::abs(id1) + 10 * std::abs(id2) +  theSpin; 
145           G4bool IsUp = (std::abs(id1)&1) == 0;        // quark 1 up type quark (u or c)
146           G4bool IsAnti = id1 < 0;             // quark 1 is antiquark?
147           if( (IsUp && IsAnti ) || (!IsUp && !IsAnti ) ) 
148              PDGEncoding = - PDGEncoding;
149        }
150           
151           
152        G4ParticleDefinition * MesonDef=
153                G4ParticleTable::GetParticleTable()->FindParticle(PDGEncoding);
154#ifdef G4VERBOSE
155        if (MesonDef == 0 ) {
156                G4cerr << " G4HadronBuilder - Warning: No particle for PDGcode= "
157                       << PDGEncoding << G4endl;
158        }
159        if  ( (  black->GetPDGCharge() 
160               + white->GetPDGCharge()
161               - MesonDef->GetPDGCharge() ) > perCent   ) {
162                G4cerr << " G4HadronBuilder - Warning: Incorrect Charge : "
163                        << " Quark1/2 = " 
164                        << black->GetParticleName() << " / "
165                        << white->GetParticleName() 
166                        << " resulting Hadron " << MesonDef->GetParticleName() 
167                        << G4endl;
168        }
169#endif
170
171        return MesonDef;
172}
173
174
175G4ParticleDefinition * G4HadronBuilder::Barion(G4ParticleDefinition * black, 
176                                              G4ParticleDefinition * white,Spin theSpin)
177{
178
179#ifdef G4VERBOSE
180//  Verify Input Charge
181   G4double charge =  black->GetPDGCharge() 
182                    + white->GetPDGCharge();     
183   if (std::abs(charge) > 2 || std::abs(3.*charge - 3*G4int(charge*1.001)) > perCent )
184        {
185            G4cerr << " G4HadronBuilder::Build()" << G4endl;
186            G4cerr << "    Invalid total charge found for on input: " 
187                        << charge<< G4endl;
188            G4cerr << "    PGDcode input quark1/quark2 : " <<
189                        black->GetPDGEncoding() << " / "<< 
190                        white->GetPDGEncoding() << G4endl;
191            G4cerr << G4endl;
192        } 
193#endif 
194        G4int id1= black->GetPDGEncoding();
195        G4int id2= white->GetPDGEncoding();
196        if ( std::abs(id1) < std::abs(id2) )
197           {
198           G4int xchg = id1; 
199           id1 = id2; 
200           id2 = xchg;
201           }
202
203        if (std::abs(id1) < 1000 || std::abs(id2) > 3 ) 
204           throw G4HadronicException(__FILE__, __LINE__, "G4HadronBuilder::Barion: Illegal quark content as input");   
205
206        G4int ifl1= std::abs(id1)/1000;
207        G4int ifl2 = (std::abs(id1) - ifl1 * 1000)/100;
208        G4int diquarkSpin = std::abs(id1)%10; 
209        G4int ifl3 = id2;
210        if (id1 < 0)
211           {
212           ifl1 = - ifl1;
213           ifl2 = - ifl2;
214           }
215        //... Construct barion, distinguish Lambda and Sigma barions.
216        G4int kfla = std::abs(ifl1);
217        G4int kflb = std::abs(ifl2);
218        G4int kflc = std::abs(ifl3);
219
220        G4int kfld = std::max(kfla,kflb);
221              kfld = std::max(kfld,kflc);
222        G4int kflf = std::min(kfla,kflb);
223              kflf = std::min(kflf,kflc);
224
225        G4int kfle = kfla + kflb + kflc - kfld - kflf;
226
227        //... barion with content uuu or ddd or sss has always spin = 3/2
228        theSpin = (kfla == kflb && kflb == kflc)? SpinThreeHalf : theSpin;   
229
230        G4int kfll = 0;
231        if(theSpin == SpinHalf && kfld > kfle && kfle > kflf) { 
232// Spin J=1/2 and all three quarks different
233// Two states exist: (uds -> lambda or sigma0)
234//   -  lambda: s(ud)0 s : 3122; ie. reverse the two lighter quarks
235//   -  sigma0: s(ud)1 s : 3212
236           if(diquarkSpin == 1 ) {
237              if ( kfla == kfld) {   // heaviest quark in diquark
238                 kfll = 1;
239              } else {
240                 kfll = (G4int)(0.25 + G4UniformRand());
241              }
242           }   
243           if(diquarkSpin == 3 && kfla != kfld)
244               kfll = (G4int)(0.75 + G4UniformRand());
245        }
246       
247        G4int PDGEncoding;
248        if (kfll == 1)
249           PDGEncoding = 1000 * kfld + 100 * kflf + 10 * kfle + theSpin;
250        else   
251           PDGEncoding = 1000 * kfld + 100 * kfle + 10 * kflf + theSpin;
252
253        if (id1 < 0)
254           PDGEncoding = -PDGEncoding;
255
256
257        G4ParticleDefinition * BarionDef=
258                G4ParticleTable::GetParticleTable()->FindParticle(PDGEncoding);
259#ifdef G4VERBOSE
260        if (BarionDef == 0 ) {
261                G4cerr << " G4HadronBuilder - Warning: No particle for PDGcode= "
262                       << PDGEncoding << G4endl;
263        }
264        if  ( (  black->GetPDGCharge() 
265               + white->GetPDGCharge()
266               - BarionDef->GetPDGCharge() ) > perCent   ) {
267                G4cerr << " G4HadronBuilder - Warning: Incorrect Charge : "
268                        << " DiQuark/Quark = " 
269                        << black->GetParticleName() << " / "
270                        << white->GetParticleName() 
271                        << " resulting Hadron " << BarionDef->GetParticleName() 
272                        << G4endl;
273        }
274#endif
275
276        return BarionDef;
277}
Note: See TracBrowser for help on using the repository browser.