Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/parton_string/hadronization/src/G4ExcitedStringDecay.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/hadronic/models/parton_string/hadronization/src/G4ExcitedStringDecay.cc (Version 11.3.0) and /processes/hadronic/models/parton_string/hadronization/src/G4ExcitedStringDecay.cc (Version 10.6.p2)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 // Historic fragment from M.Komogorov; clean-u     26 // Historic fragment from M.Komogorov; clean-up still necessary @@@
 27                                                    27 
 28 #include "G4ExcitedStringDecay.hh"                 28 #include "G4ExcitedStringDecay.hh"
 29 #include "G4SystemOfUnits.hh"                      29 #include "G4SystemOfUnits.hh"
 30 #include "G4KineticTrack.hh"                       30 #include "G4KineticTrack.hh"
 31 #include "G4LundStringFragmentation.hh"        << 
 32 #include "G4HadronicInteractionRegistry.hh"    << 
 33                                                    31 
 34 #include "G4SampleResonance.hh"                    32 #include "G4SampleResonance.hh"
 35                                                    33 
 36 //#define debug_G4ExcitedStringDecay               34 //#define debug_G4ExcitedStringDecay
 37 //#define debug_G4ExcitedStringCorr                35 //#define debug_G4ExcitedStringCorr
 38                                                    36 
 39 G4ExcitedStringDecay::G4ExcitedStringDecay(G4V <<  37 
 40   : G4VStringFragmentation(), theStringDecay(p <<  38 G4ExcitedStringDecay::G4ExcitedStringDecay() : G4VStringFragmentation(),theStringDecay(0)
                                                   >>  39 {}
                                                   >>  40 
                                                   >>  41 
                                                   >>  42 G4ExcitedStringDecay::G4ExcitedStringDecay(G4VLongitudinalStringDecay * aStringDecay)
                                                   >>  43 : G4VStringFragmentation(),
                                                   >>  44   theStringDecay(aStringDecay)
                                                   >>  45 {}
                                                   >>  46 
                                                   >>  47 
                                                   >>  48 G4ExcitedStringDecay::G4ExcitedStringDecay(const G4ExcitedStringDecay &)
                                                   >>  49 : G4VStringFragmentation(),
                                                   >>  50   theStringDecay(0)
 41 {                                                  51 {
 42   if(!ptr) {                                   <<  52   throw G4HadronicException(__FILE__, __LINE__, "G4ExcitedStringDecay::copy ctor not accessible");
 43     G4HadronicInteraction* p =                 <<  53 } 
 44       G4HadronicInteractionRegistry::Instance( <<  54 
 45     theStringDecay = static_cast<G4VLongitudin << 
 46     if(!theStringDecay) { theStringDecay = new << 
 47   }                                            << 
 48   SetModelName(theStringDecay->GetModelName()) << 
 49 }                                              << 
 50                                                    55 
 51 G4ExcitedStringDecay::~G4ExcitedStringDecay()      56 G4ExcitedStringDecay::~G4ExcitedStringDecay()
 52 {}                                             <<  57 {
                                                   >>  58 }
                                                   >>  59 
                                                   >>  60 
                                                   >>  61 const G4ExcitedStringDecay & G4ExcitedStringDecay::operator=(const G4ExcitedStringDecay &)
                                                   >>  62 {
                                                   >>  63   throw G4HadronicException(__FILE__, __LINE__, "G4ExcitedStringDecay::operator= meant to not be accessable");
                                                   >>  64   return *this;
                                                   >>  65 }
                                                   >>  66 
                                                   >>  67 
                                                   >>  68 G4bool G4ExcitedStringDecay::operator==(const G4ExcitedStringDecay &) const
                                                   >>  69 {
                                                   >>  70   return false;
                                                   >>  71 }
                                                   >>  72 
                                                   >>  73 
                                                   >>  74 G4bool G4ExcitedStringDecay::operator!=(const G4ExcitedStringDecay &) const
                                                   >>  75 {
                                                   >>  76   return true;
                                                   >>  77 }
                                                   >>  78 
 53                                                    79 
 54 G4KineticTrackVector *G4ExcitedStringDecay::Fr     80 G4KineticTrackVector *G4ExcitedStringDecay::FragmentString(const G4ExcitedString &theString)
 55 {                                                  81 {
 56   return theStringDecay->FragmentString(theStr <<  82  if ( theStringDecay == NULL ) theStringDecay=new G4LundStringFragmentation();
                                                   >>  83  return theStringDecay->FragmentString(theString);
 57 }                                                  84 }
 58                                                    85 
                                                   >>  86   
 59 G4KineticTrackVector *G4ExcitedStringDecay::Fr     87 G4KineticTrackVector *G4ExcitedStringDecay::FragmentStrings(const G4ExcitedStringVector * theStrings)
 60 {                                                  88 {
 61   G4LorentzVector KTsum(0.,0.,0.,0.);              89   G4LorentzVector KTsum(0.,0.,0.,0.);
 62                                                    90 
 63   #ifdef debug_G4ExcitedStringDecay                91   #ifdef debug_G4ExcitedStringDecay
 64   G4cout<<G4endl;                                  92   G4cout<<G4endl;
 65   G4cout<<"--------------------------- G4Excit     93   G4cout<<"--------------------------- G4ExcitedStringDecay ----------------------"<<G4endl;
 66   G4cout<<"Hadronization of Excited Strings: t     94   G4cout<<"Hadronization of Excited Strings: theStrings->size() "<<theStrings->size()<<G4endl;
 67   #endif                                           95   #endif
 68                                                    96 
 69   for ( unsigned int astring=0; astring < theS     97   for ( unsigned int astring=0; astring < theStrings->size(); astring++)
 70   // for ( unsigned int astring=0; astring < 1     98   // for ( unsigned int astring=0; astring < 1; astring++)
 71   {                                                99   {
 72     // G4cout<<"theStrings->operator[](astring    100     // G4cout<<"theStrings->operator[](astring)->IsExcited() "<<" "<<astring<<" "<<theStrings->operator[](astring)->IsExcited()<<G4endl;
 73     if ( theStrings->operator[](astring)->IsEx    101     if ( theStrings->operator[](astring)->IsExcited() )
 74          {KTsum+= theStrings->operator[](astri    102          {KTsum+= theStrings->operator[](astring)->Get4Momentum();}
 75     else {KTsum+=theStrings->operator[](astrin    103     else {KTsum+=theStrings->operator[](astring)->GetKineticTrack()->Get4Momentum();}
 76   }                                               104   }
 77                                                   105 
 78   G4LorentzRotation toCms( -1 * KTsum.boostVec    106   G4LorentzRotation toCms( -1 * KTsum.boostVector() );
 79   G4LorentzRotation toLab(toCms.inverse());       107   G4LorentzRotation toLab(toCms.inverse());
 80   G4LorentzVector Ptmp;                           108   G4LorentzVector Ptmp;
 81   KTsum=G4LorentzVector(0.,0.,0.,0.);             109   KTsum=G4LorentzVector(0.,0.,0.,0.);
 82                                                   110 
 83   for ( unsigned int astring=0; astring < theS    111   for ( unsigned int astring=0; astring < theStrings->size(); astring++)
 84   // for ( unsigned int astring=0; astring < 1    112   // for ( unsigned int astring=0; astring < 1; astring++)
 85   {                                               113   {
 86     if ( theStrings->operator[](astring)->IsEx    114     if ( theStrings->operator[](astring)->IsExcited() )
 87     {                                             115     {
 88       Ptmp=toCms * theStrings->operator[](astr    116       Ptmp=toCms * theStrings->operator[](astring)->GetLeftParton()->Get4Momentum();
 89       theStrings->operator[](astring)->GetLeft    117       theStrings->operator[](astring)->GetLeftParton()->Set4Momentum(Ptmp);
 90                                                   118 
 91       Ptmp=toCms * theStrings->operator[](astr    119       Ptmp=toCms * theStrings->operator[](astring)->GetRightParton()->Get4Momentum();
 92       theStrings->operator[](astring)->GetRigh    120       theStrings->operator[](astring)->GetRightParton()->Set4Momentum(Ptmp);
 93                                                   121 
 94       KTsum+= theStrings->operator[](astring)-    122       KTsum+= theStrings->operator[](astring)->Get4Momentum();
 95     }                                             123     }
 96     else                                          124     else
 97     {                                             125     {
 98       Ptmp=toCms * theStrings->operator[](astr    126       Ptmp=toCms * theStrings->operator[](astring)->GetKineticTrack()->Get4Momentum();
 99       theStrings->operator[](astring)->GetKine    127       theStrings->operator[](astring)->GetKineticTrack()->Set4Momentum(Ptmp);
100       KTsum+= theStrings->operator[](astring)-    128       KTsum+= theStrings->operator[](astring)->GetKineticTrack()->Get4Momentum();
101     }                                             129     }
102   }                                               130   }
103                                                   131 
104   G4SampleResonance BrW;                          132   G4SampleResonance BrW;
105   const G4ParticleDefinition* TrackDefinition=    133   const G4ParticleDefinition* TrackDefinition=0;
106                                                   134 
107   G4KineticTrackVector * theResult = new G4Kin    135   G4KineticTrackVector * theResult = new G4KineticTrackVector;
108   G4int attempts(0);                              136   G4int attempts(0);
109   G4bool success=false;                           137   G4bool success=false;
110   G4bool NeedEnergyCorrector=false;               138   G4bool NeedEnergyCorrector=false;
111   do {                                            139   do {
112         #ifdef debug_G4ExcitedStringDecay         140         #ifdef debug_G4ExcitedStringDecay  
113         G4cout<<"New try No "<<attempts<<" to     141         G4cout<<"New try No "<<attempts<<" to hadronize strings"<<G4endl;
114         #endif                                    142         #endif
115                                                   143 
116   std::for_each(theResult->begin() , theResult    144   std::for_each(theResult->begin() , theResult->end() , DeleteKineticTrack());
117   theResult->clear();                             145   theResult->clear();
118                                                   146 
119   attempts++;                                     147   attempts++;
120                                                   148 
121   G4LorentzVector KTsecondaries(0.,0.,0.,0.);     149   G4LorentzVector KTsecondaries(0.,0.,0.,0.);
122   NeedEnergyCorrector=false;                      150   NeedEnergyCorrector=false;
123                                                   151 
124   for ( unsigned int astring=0; astring < theS    152   for ( unsigned int astring=0; astring < theStrings->size(); astring++)
125         // for ( unsigned int astring=0; astri    153         // for ( unsigned int astring=0; astring < 1; astring++)  // Proj. Last Str. Decay for FTF
126         // for ( unsigned int astring=1; astri    154         // for ( unsigned int astring=1; astring < 2; astring++)  // Proj. Last Str. Decay for QGS
127         // for ( unsigned int astring=0; astri    155         // for ( unsigned int astring=0; astring < 1; astring++)  // For testing purposes
128   {                                               156   {
129           #ifdef debug_G4ExcitedStringDecay       157           #ifdef debug_G4ExcitedStringDecay  
130           G4cout<<"String No "<<astring+1<<" E    158           G4cout<<"String No "<<astring+1<<" Excited?  "<<theStrings->operator[](astring)->IsExcited()<<G4endl;
131                                                   159 
132           G4cout<<"String No "<<astring+1<<" 4    160           G4cout<<"String No "<<astring+1<<" 4Momentum "<<theStrings->operator[](astring)->Get4Momentum()
133           <<" "<<theStrings->operator[](astrin    161           <<" "<<theStrings->operator[](astring)->Get4Momentum().mag()<<G4endl;
134           #endif                                  162           #endif
135                                                   163 
136           G4KineticTrackVector * generatedKine    164           G4KineticTrackVector * generatedKineticTracks = NULL;
137     if ( theStrings->operator[](astring)->IsEx    165     if ( theStrings->operator[](astring)->IsExcited() )
138     {                                             166     {
139              #ifdef debug_G4ExcitedStringDecay    167              #ifdef debug_G4ExcitedStringDecay  
140              G4cout<<"Fragment String with par    168              G4cout<<"Fragment String with partons: "
141                    <<theStrings->operator[](as    169                    <<theStrings->operator[](astring)->GetLeftParton()->GetPDGcode() <<" "
142                    <<theStrings->operator[](as    170                    <<theStrings->operator[](astring)->GetRightParton()->GetPDGcode()<<" "
143                    <<"Direction "<<theStrings-    171                    <<"Direction "<<theStrings->operator[](astring)->GetDirection()<<G4endl;
144              #endif                               172              #endif
145          generatedKineticTracks=FragmentString    173          generatedKineticTracks=FragmentString(*theStrings->operator[](astring));
146              #ifdef debug_G4ExcitedStringDecay    174              #ifdef debug_G4ExcitedStringDecay  
147              G4cout<<"(G4ExcitedStringDecay) N    175              G4cout<<"(G4ExcitedStringDecay) Number of produced hadrons = "
148                    <<generatedKineticTracks->s    176                    <<generatedKineticTracks->size()<<G4endl;
149              #endif                               177              #endif
150     } else {                                      178     } else {
151              #ifdef debug_G4ExcitedStringDecay    179              #ifdef debug_G4ExcitedStringDecay  
152              G4cout<<"   GetTrack from the Str    180              G4cout<<"   GetTrack from the String"<<G4endl;
153              #endif                               181              #endif
154              G4LorentzVector Mom=theStrings->o    182              G4LorentzVector Mom=theStrings->operator[](astring)->GetKineticTrack()->Get4Momentum();
155              G4KineticTrack * aTrack= new G4Ki    183              G4KineticTrack * aTrack= new G4KineticTrack(
156                             theStrings->operat    184                             theStrings->operator[](astring)->GetKineticTrack()->GetDefinition(),
157                             theStrings->operat    185                             theStrings->operator[](astring)->GetKineticTrack()->GetFormationTime(),
158                             G4ThreeVector(0),     186                             G4ThreeVector(0), Mom);
159                                                   187 
160              aTrack->SetPosition(theStrings->o    188              aTrack->SetPosition(theStrings->operator[](astring)->GetKineticTrack()->GetPosition());
161                                                   189 
162              #ifdef debug_G4ExcitedStringDecay    190              #ifdef debug_G4ExcitedStringDecay  
163              G4cout<<"   A particle stored in     191              G4cout<<"   A particle stored in the track is "<<aTrack->GetDefinition()->GetParticleName()<<G4endl;
164              #endif                               192              #endif
165                                                   193 
166        generatedKineticTracks = new G4KineticT    194        generatedKineticTracks = new G4KineticTrackVector;
167        generatedKineticTracks->push_back(aTrac    195        generatedKineticTracks->push_back(aTrack);
168     }                                             196     }    
169                                                   197 
170     if (generatedKineticTracks == nullptr || g << 198     if (generatedKineticTracks->size() == 0)
171     {                                             199     {
172              // G4cerr << "G4VPartonStringMode    200              // G4cerr << "G4VPartonStringModel:No KineticTracks produced" << G4endl; 
173        delete generatedKineticTracks;          << 201              // continue;                                                             
174        generatedKineticTracks = nullptr;       << 202              success=false; NeedEnergyCorrector=false; break;                      
175              continue;                         << 
176     }                                             203     }
177                                                   204 
178           G4LorentzVector KTsum1(0.,0.,0.,0.);    205           G4LorentzVector KTsum1(0.,0.,0.,0.);
179           for ( unsigned int aTrack=0; aTrack<    206           for ( unsigned int aTrack=0; aTrack<generatedKineticTracks->size();aTrack++)
180     {                                             207     {
181              #ifdef debug_G4ExcitedStringDecay    208              #ifdef debug_G4ExcitedStringDecay  
182              G4cout<<"Prod part No. "<<aTrack+    209              G4cout<<"Prod part No. "<<aTrack+1<<" "
183                    <<(*generatedKineticTracks)    210                    <<(*generatedKineticTracks)[aTrack]->GetDefinition()->GetParticleName()<<" "
184                    <<(*generatedKineticTracks)    211                    <<(*generatedKineticTracks)[aTrack]->Get4Momentum()
185                    <<(*generatedKineticTracks)    212                    <<(*generatedKineticTracks)[aTrack]->Get4Momentum().mag()<<G4endl;
186              #endif                               213              #endif
187              // --------------- Sampling mass     214              // --------------- Sampling mass of unstable hadronic resonances ----------------
188              TrackDefinition = (*generatedKine    215              TrackDefinition = (*generatedKineticTracks)[aTrack]->GetDefinition();
189                                                   216 
190        if (TrackDefinition->IsShortLived())       217        if (TrackDefinition->IsShortLived())
191              {                                    218              {
192                G4double NewTrackMass =            219                G4double NewTrackMass = 
193                  BrW.SampleMass( TrackDefiniti    220                  BrW.SampleMass( TrackDefinition->GetPDGMass(), TrackDefinition->GetPDGWidth(),
194                                  BrW.GetMinimu    221                                  BrW.GetMinimumMass( TrackDefinition ) + 10.0*MeV, 
195                                  TrackDefiniti    222                                  TrackDefinition->GetPDGMass() + 5.0*TrackDefinition->GetPDGWidth() );
196                G4LorentzVector Tmp=G4LorentzVe    223                G4LorentzVector Tmp=G4LorentzVector((*generatedKineticTracks)[aTrack]->Get4Momentum());
197                Tmp.setE(std::sqrt(sqr(NewTrack    224                Tmp.setE(std::sqrt(sqr(NewTrackMass) + Tmp.vect().mag2()));
198                                                   225 
199                (*generatedKineticTracks)[aTrac    226                (*generatedKineticTracks)[aTrack]->Set4Momentum(Tmp);
200                                                   227 
201                #ifdef debug_G4ExcitedStringDec    228                #ifdef debug_G4ExcitedStringDecay  
202                G4cout<<"Resonance *** "<<aTrac    229                G4cout<<"Resonance *** "<<aTrack+1<<" "
203                      <<(*generatedKineticTrack    230                      <<(*generatedKineticTracks)[aTrack]->GetDefinition()->GetParticleName()<<" "
204                      <<(*generatedKineticTrack    231                      <<(*generatedKineticTracks)[aTrack]->Get4Momentum()
205                      <<(*generatedKineticTrack    232                      <<(*generatedKineticTracks)[aTrack]->Get4Momentum().mag()<<G4endl;
206                #endif                             233                #endif
207            }                                      234            }
208              //-------------------------------    235              //-------------------------------------------------------------------------------
209                                                   236 
210              theResult->push_back(generatedKin    237              theResult->push_back(generatedKineticTracks->operator[](aTrack));
211              KTsum1+= (*generatedKineticTracks    238              KTsum1+= (*generatedKineticTracks)[aTrack]->Get4Momentum();
212     }                                             239     }
213     KTsecondaries+=KTsum1;                        240     KTsecondaries+=KTsum1;
214                                                   241   
215           #ifdef debug_G4ExcitedStringDecay       242           #ifdef debug_G4ExcitedStringDecay  
216           G4cout << "String secondaries(" <<ge    243           G4cout << "String secondaries(" <<generatedKineticTracks->size()<< ")"<<G4endl
217                  <<"Init  string  momentum: "<    244                  <<"Init  string  momentum: "<< theStrings->operator[](astring)->Get4Momentum()<<G4endl
218                  <<"Final hadrons momentum: "<    245                  <<"Final hadrons momentum: "<< KTsum1 << G4endl;
219           #endif                                  246           #endif
220                                                   247 
221     if  ( KTsum1.e() > 0 &&                       248     if  ( KTsum1.e() > 0 && 
222                 std::abs((KTsum1.e()-theString    249                 std::abs((KTsum1.e()-theStrings->operator[](astring)->Get4Momentum().e()) / KTsum1.e()) > perMillion )
223     {                                             250     {
224       NeedEnergyCorrector=true;                   251       NeedEnergyCorrector=true;
225     }                                             252     }
226                                                   253 
227           #ifdef debug_G4ExcitedStringDecay       254           #ifdef debug_G4ExcitedStringDecay  
228           G4cout<<"NeedEnergyCorrection yes/no    255           G4cout<<"NeedEnergyCorrection yes/no "<<NeedEnergyCorrector<<G4endl;
229           #endif                                  256           #endif
230                                                   257 
231           // clean up                             258           // clean up
232     delete generatedKineticTracks;                259     delete generatedKineticTracks;
233     success=true;                                 260     success=true;                          
234   }                                               261   }
235                                                   262 
236   if ( NeedEnergyCorrector ) success=EnergyAnd    263   if ( NeedEnergyCorrector ) success=EnergyAndMomentumCorrector(theResult, KTsum);
237   } while (!success && (attempts < 100));  /*     264   } while (!success && (attempts < 100));  /* Loop checking, 07.08.2015, A.Ribon */
238                                                   265 
239   for ( unsigned int aTrack=0; aTrack<theResul    266   for ( unsigned int aTrack=0; aTrack<theResult->size();aTrack++)
240   {                                               267   {
241     Ptmp=(*theResult)[aTrack]->Get4Momentum();    268     Ptmp=(*theResult)[aTrack]->Get4Momentum();
242     Ptmp.transform( toLab);                       269     Ptmp.transform( toLab);
243     (*theResult)[aTrack]->Set4Momentum(Ptmp);     270     (*theResult)[aTrack]->Set4Momentum(Ptmp);
244   }                                               271   }
245                                                   272 
246   #ifdef debug_G4ExcitedStringDecay               273   #ifdef debug_G4ExcitedStringDecay  
247   G4cout<<"End of the strings fragmentation (G    274   G4cout<<"End of the strings fragmentation (G4ExcitedStringDecay)"<<G4endl;
248                                                   275 
249   G4LorentzVector  KTsum1(0.,0.,0.,0.);           276   G4LorentzVector  KTsum1(0.,0.,0.,0.); 
250                                                   277 
251   for ( unsigned int aTrack=0; aTrack<theResul    278   for ( unsigned int aTrack=0; aTrack<theResult->size();aTrack++)
252   {                                               279   {
253       G4cout << " corrected tracks .. " << (*t    280       G4cout << " corrected tracks .. " << (*theResult)[aTrack]->GetDefinition()->GetParticleName()
254       <<"  " << (*theResult)[aTrack]->Get4Mome    281       <<"  " << (*theResult)[aTrack]->Get4Momentum() 
255       <<"  " << (*theResult)[aTrack]->Get4Mome    282       <<"  " << (*theResult)[aTrack]->Get4Momentum().mag()<< G4endl;
256       KTsum1+= (*theResult)[aTrack]->Get4Momen    283       KTsum1+= (*theResult)[aTrack]->Get4Momentum();
257   }                                               284   }
258                                                   285 
259   G4cout << "Needcorrector/success " << NeedEn    286   G4cout << "Needcorrector/success " << NeedEnergyCorrector << "/" << success 
260          << ", Corrected total  4 momentum " <    287          << ", Corrected total  4 momentum " << KTsum1  << G4endl;
261   if ( ! success ) G4cout << "failed to correc    288   if ( ! success ) G4cout << "failed to correct E/p" << G4endl;  
262                                                   289 
263   G4cout<<"End of the Hadronization (G4Excited    290   G4cout<<"End of the Hadronization (G4ExcitedStringDecay)"<<G4endl;
264   #endif                                          291   #endif
265                                                   292 
266   if (!success)                                   293   if (!success)
267   {                                               294   {
268     if (theResult->size() != 0)                   295     if (theResult->size() != 0)
269     {                                             296     {
270       std::for_each(theResult->begin() , theRe    297       std::for_each(theResult->begin() , theResult->end() , DeleteKineticTrack());
271       theResult->clear();                         298       theResult->clear();
272       delete theResult; theResult=0;              299       delete theResult; theResult=0;
273     }                                             300     }
274     for ( unsigned int astring=0; astring < th    301     for ( unsigned int astring=0; astring < theStrings->size(); astring++)
275     // for ( unsigned int astring=0; astring <    302     // for ( unsigned int astring=0; astring < 1; astring++) // Need more correct. For testing purposes.
276     {                                             303     {
277       if ( theStrings->operator[](astring)->Is    304       if ( theStrings->operator[](astring)->IsExcited() )
278       {                                           305       {
279         Ptmp=theStrings->operator[](astring)->    306         Ptmp=theStrings->operator[](astring)->GetLeftParton()->Get4Momentum();
280         Ptmp.transform( toLab);                   307         Ptmp.transform( toLab);
281         theStrings->operator[](astring)->GetLe    308         theStrings->operator[](astring)->GetLeftParton()->Set4Momentum(Ptmp);
282                                                   309 
283         Ptmp=theStrings->operator[](astring)->    310         Ptmp=theStrings->operator[](astring)->GetRightParton()->Get4Momentum();
284         Ptmp.transform( toLab);                   311         Ptmp.transform( toLab);
285         theStrings->operator[](astring)->GetRi    312         theStrings->operator[](astring)->GetRightParton()->Set4Momentum(Ptmp);
286       }                                           313       }
287       else                                        314       else
288       {                                           315       {
289         Ptmp=theStrings->operator[](astring)->    316         Ptmp=theStrings->operator[](astring)->GetKineticTrack()->Get4Momentum();
290         Ptmp.transform( toLab);                   317         Ptmp.transform( toLab);
291         theStrings->operator[](astring)->GetKi    318         theStrings->operator[](astring)->GetKineticTrack()->Set4Momentum(Ptmp);
292       }                                           319       }
293     }                                             320     }
294   }                                               321   }
295   return theResult;                               322   return theResult;
296 }                                                 323 }
297                                                   324 
298                                                   325 
299 G4bool G4ExcitedStringDecay::                     326 G4bool G4ExcitedStringDecay::
300 EnergyAndMomentumCorrector(G4KineticTrackVecto    327 EnergyAndMomentumCorrector(G4KineticTrackVector* Output, G4LorentzVector& TotalCollisionMom)   
301 {                                                 328 {
302     const int    nAttemptScale = 500;             329     const int    nAttemptScale = 500;
303     const double ErrLimit = 1.E-5;                330     const double ErrLimit = 1.E-5;
304     if (Output->empty()) return TRUE;             331     if (Output->empty()) return TRUE;
305     G4LorentzVector SumMom;                       332     G4LorentzVector SumMom;
306     G4double        SumMass = 0;                  333     G4double        SumMass = 0;
307     G4double        TotalCollisionMass = Total    334     G4double        TotalCollisionMass = TotalCollisionMom.m();
308                                                   335 
309     std::vector<G4double> HadronMass; G4double    336     std::vector<G4double> HadronMass; G4double HadronM(0.);
310                                                   337 
311     #ifdef debug_G4ExcitedStringCorr              338     #ifdef debug_G4ExcitedStringCorr
312     G4cout<<G4endl<<"EnergyAndMomentumCorrecto    339     G4cout<<G4endl<<"EnergyAndMomentumCorrector. Number of particles: "<<Output->size()<<G4endl;
313     #endif                                        340     #endif
314     // Calculate sum hadron 4-momenta and summ    341     // Calculate sum hadron 4-momenta and summing hadron mass
315     unsigned int cHadron;                         342     unsigned int cHadron;
316     for (cHadron = 0; cHadron < Output->size()    343     for (cHadron = 0; cHadron < Output->size(); cHadron++)
317     {                                             344     {
318         SumMom  += Output->operator[](cHadron)    345         SumMom  += Output->operator[](cHadron)->Get4Momentum();
319         HadronM=Output->operator[](cHadron)->G    346         HadronM=Output->operator[](cHadron)->Get4Momentum().mag(); HadronMass.push_back(HadronM);
320         SumMass += Output->operator[](cHadron)    347         SumMass += Output->operator[](cHadron)->Get4Momentum().mag();  //GetDefinition()->GetPDGMass();
321     }                                             348     }
322                                                   349 
323     #ifdef debug_G4ExcitedStringCorr              350     #ifdef debug_G4ExcitedStringCorr
324     G4cout<<"Sum part mom "<<SumMom<<" "<<SumM    351     G4cout<<"Sum part mom "<<SumMom<<" "<<SumMom.mag()<<G4endl
325           <<"Sum str  mom "<<TotalCollisionMom    352           <<"Sum str  mom "<<TotalCollisionMom<<" "<<TotalCollisionMom.mag()<<G4endl;
326     G4cout<<"SumMass TotalCollisionMass "<<Sum    353     G4cout<<"SumMass TotalCollisionMass "<<SumMass<<" "<<TotalCollisionMass<<G4endl;
327     #endif                                        354     #endif
328                                                   355 
329     // Cannot correct a single particle           356     // Cannot correct a single particle
330     if (Output->size() < 2) return FALSE;         357     if (Output->size() < 2) return FALSE;
331                                                   358 
332     if (SumMass > TotalCollisionMass) return F    359     if (SumMass > TotalCollisionMass) return FALSE;
333     SumMass = SumMom.m2();                        360     SumMass = SumMom.m2();
334     if (SumMass < 0) return FALSE;                361     if (SumMass < 0) return FALSE;
335     SumMass = std::sqrt(SumMass);                 362     SumMass = std::sqrt(SumMass);
336                                                   363 
337     // Compute c.m.s. hadron velocity and boos    364     // Compute c.m.s. hadron velocity and boost KTV to hadron c.m.s.
338     // G4ThreeVector Beta = -SumMom.boostVecto    365     // G4ThreeVector Beta = -SumMom.boostVector();          
339     G4ThreeVector Beta = -TotalCollisionMom.bo    366     G4ThreeVector Beta = -TotalCollisionMom.boostVector(); 
340     Output->Boost(Beta);                          367     Output->Boost(Beta);
341                                                   368 
342     // Scale total c.m.s. hadron energy (hadro    369     // Scale total c.m.s. hadron energy (hadron system mass).
343     // It should be equal interaction mass        370     // It should be equal interaction mass
344     G4double Scale = 1;                           371     G4double Scale = 1;
345     G4int cAttempt = 0;                           372     G4int cAttempt = 0;
346     G4double Sum = 0;                             373     G4double Sum = 0;
347     G4bool success = false;                       374     G4bool success = false;
348     for (cAttempt = 0; cAttempt < nAttemptScal    375     for (cAttempt = 0; cAttempt < nAttemptScale; cAttempt++)
349     {                                             376     {
350       Sum = 0;                                    377       Sum = 0;
351       for (cHadron = 0; cHadron < Output->size    378       for (cHadron = 0; cHadron < Output->size(); cHadron++)
352       {                                           379       {
353         HadronM = HadronMass.at(cHadron);         380         HadronM = HadronMass.at(cHadron);
354         G4LorentzVector HadronMom = Output->op    381         G4LorentzVector HadronMom = Output->operator[](cHadron)->Get4Momentum();
355         HadronMom.setVect(Scale*HadronMom.vect    382         HadronMom.setVect(Scale*HadronMom.vect());
356         G4double E = std::sqrt(HadronMom.vect(    383         G4double E = std::sqrt(HadronMom.vect().mag2() + sqr(HadronM));
357                                                   384                                                        //sqr(Output->operator[](cHadron)->GetDefinition()->GetPDGMass()));
358         HadronMom.setE(E);                        385         HadronMom.setE(E);
359         Output->operator[](cHadron)->Set4Momen    386         Output->operator[](cHadron)->Set4Momentum(HadronMom);
360         Sum += E;                                 387         Sum += E;
361       }                                           388       } 
362       Scale = TotalCollisionMass/Sum;             389       Scale = TotalCollisionMass/Sum;
363       #ifdef debug_G4ExcitedStringCorr            390       #ifdef debug_G4ExcitedStringCorr
364       G4cout << "Scale-1=" << Scale -1            391       G4cout << "Scale-1=" << Scale -1 
365              << ",  TotalCollisionMass=" << To    392              << ",  TotalCollisionMass=" << TotalCollisionMass
366              << ",  Sum=" << Sum                  393              << ",  Sum=" << Sum
367        << G4endl;                                 394        << G4endl;
368       #endif                                      395       #endif
369       if (std::fabs(Scale - 1) <= ErrLimit)       396       if (std::fabs(Scale - 1) <= ErrLimit) 
370       {                                           397       {
371         success = true;                           398         success = true;
372   break;                                          399   break;
373       }                                           400       }
374     }                                             401     }
375                                                   402 
376     #ifdef debug_G4ExcitedStringCorr              403     #ifdef debug_G4ExcitedStringCorr
377     if (!success)                                 404     if (!success)
378     {                                             405     {
379       G4cout << "G4ExcitedStringDecay::EnergyA    406       G4cout << "G4ExcitedStringDecay::EnergyAndMomentumCorrector - Warning"<<G4endl;
380       G4cout << "   Scale not unity at end of     407       G4cout << "   Scale not unity at end of iteration loop: "<<TotalCollisionMass<<" "<<Sum<<" "<<Scale<<G4endl;
381       G4cout << "   Number of secondaries: " <    408       G4cout << "   Number of secondaries: " << Output->size() << G4endl;
382       G4cout << "   Wanted total energy: " <<     409       G4cout << "   Wanted total energy: " <<  TotalCollisionMom.e() << G4endl; 
383       G4cout << "   Increase number of attempt    410       G4cout << "   Increase number of attempts or increase ERRLIMIT"<<G4endl;
384       // throw G4HadronicException(__FILE__, _    411       // throw G4HadronicException(__FILE__, __LINE__, "G4ExcitedStringDecay failed to correct...");
385     }                                             412     }
386     #endif                                        413     #endif     
387                                                   414 
388     // Compute c.m.s. interaction velocity and    415     // Compute c.m.s. interaction velocity and KTV back boost   
389     Beta = TotalCollisionMom.boostVector();       416     Beta = TotalCollisionMom.boostVector();
390     Output->Boost(Beta);                          417     Output->Boost(Beta);
391                                                   418 
392     return success;                               419     return success;
393 }                                                 420 }
394                                                   421 
395                                                   422