Changeset 962 for trunk/source/processes/hadronic/models/parton_string

trunk/source/processes/hadronic/models/parton_string/History

-                      r819
+                      r962
      * Please list in reverse chronological order (last date on top)
      ---------------------------------------------------------------
+June 2008 V. Uzhinsky (hadr-prtn-V09-01-02)
+   G4QGSMFragmentation.cc and G4LundStringFragmentation.cc were
+   decoupled at calculation of formation time due to adding methods
+   and member in G4VLongitudinalStringDecay class for manipulation
+   with string tension.
+   Memory leak was erased in G4QGSMFragmentation.cc and
+   G4VLongitudinalStringDecay.cc thank to Gunter.
+June 2008 V. Uzhinsky (hadr-prtn-V09-01-01)
+   Changes in G4ExcitedString class were tagged. They are needed for
+   operation of FTF.
+June 2008 V. Uzhinsky    (hadr-prtn-V09-01-00)
+-----------------------------------------------
+. String fragmentation was revised, and parameters were tuned.
+. FTF parameters were tuned for proton-proton interaction
+. FTF parameters for pion-nucleon interactions were determined very rouhgly
+. Quiasi-elastic hadron-nucleus scattering was implemented in FTF
+. Formation time was implemented in FTF, and string tension was tuned
 May 2007 G.Folger       (hadr-prtn-V08-02-02)
 …
 - D. Wright created History file for parton_string directory
+March 2008 V. Uzhinsky (hadr-string-diff-V09-01-00)
+--------------------------------------------------------
+- Elastic hadron intra-nuclear nucleon scattering was inserted in
+  FTF model. This allows to simulate quasi-elastic and multi-particles
+  production together.
+- Small re-orangement of FTF model was done. G4FTFCrossSection modules
+  were re-named into G4FTFParameters and moved to /diffraction

trunk/source/processes/hadronic/models/parton_string/diffraction/History

-                      r819
+                      r962
 $Id: History,v 1.1 2007/04/24 10:32:59 gunter Exp $
+$Id: History,v 1.5 2008/12/09 10:43:47 vuzhinsk Exp $
 -------------------------------------------------------------------
 …
      * Please list in reverse chronological order (last date on top)
      ---------------------------------------------------------------
+December 08, V. Uzhinsky (hadr-string-diff-V09-01-04)
+- Improvement of delete operators in FTF
+December 08, V. Uzhinsky (hadr-string-diff-V09-01-03)
+- Some objects did not erase in FTFModel desructor. These lead to memory
+  leak.
+June 08, G.Folger             (hadr-string-diff-V09-01-02)
+-  on branch geant4-09-01-patches_branch, fix compilation warning for unused
+     variables in G4FTFModel.cc
+??????????????????????????????      (hadr-string-diff-V09-01-01)
+March 2008 V. Uzhinsky           Tag : hadr-string-diff-V09-01-00
+- G4FTFParameters.cc and G4FTFParameters.hh were copied from G4FTFCrossSection
+  corresponding files.
+- New files - G4ElasticHNScattering have been added. They implement elastic
+  scattering of hadron in intra-nuclear collisions in FTF model.
+- The corresponding changes have been done in G4FTFModel.cc and
+  G4FTFParticipants.cc
 Apr 2007 Gunter Folger  (hadr-string-diff-V08-02-00)
 -  merge in change done by ftf development.
 - Created History file.

trunk/source/processes/hadronic/models/parton_string/diffraction/include/G4DiffractiveExcitation.hh

-                      r819
+                      r962
 //
 //
 // $Id: G4DiffractiveExcitation.hh,v 1.1 2007/05/25 06:56:53 vuzhinsk Exp $
+// $Id: G4DiffractiveExcitation.hh,v 1.2 2008/04/25 14:20:13 vuzhinsk Exp $
 #ifndef G4DiffractiveExcitation_h
 …
 class G4VSplitableHadron;
 class G4ExcitedString;
+#include "G4FTFParameters.hh"                            // Uzhi 19.04.08
 #include "G4ThreeVector.hh"
 …
       virtual ~G4DiffractiveExcitation();
+      virtual G4bool ExciteParticipants (G4VSplitableHadron *aPartner, G4VSplitableHadron * bPartner) const;
+      virtual G4bool ExciteParticipants (G4VSplitableHadron *aPartner,
+                                         G4VSplitableHadron * bPartner,
+                                         G4FTFParameters *theParameters) const;
       virtual G4ExcitedString * String(G4VSplitableHadron * aHadron, G4bool isProjectile) const;
-//      void SetPtWidth(G4double aValue) { widthOfPtSquare = aValue*aValue; }
-//      void SetExtraMass(G4double aValue) { minExtraMass = aValue; }
-//      void SetMinimumMass(G4double aValue) { minmass = aValue; }
   private:
 …
       G4DiffractiveExcitation(const G4DiffractiveExcitation &right);
 //    G4double ChooseX(G4double Xmin, G4double Xmax) const;                       // Uzhi
+      G4ThreeVector GaussianPt(G4double  AveragePt2, G4double maxPtSquare) const; // Uzhi
       G4double ChooseP(G4double Pmin, G4double Pmax) const;                       // Uzhi
-//    G4ThreeVector GaussianPt(G4double widthSquare, G4double maxPtSquare) const;
-      G4ThreeVector GaussianPt(G4double  AveragePt2, G4double maxPtSquare) const;  // Uzhi
       const G4DiffractiveExcitation & operator=(const G4DiffractiveExcitation &right);
 …
       int operator!=(const G4DiffractiveExcitation &right) const;
-  private:
-// Model Parameters:
-/*                                                                                // Uzhi
-        const G4double widthOfPtSquare; // width^2 of pt for string excitation
-        const G4double minExtraMass;    // minimum excitation mass
-        const G4double minmass; // mean pion transverse mass; used for Xmin
-*/                                                                                // Uzhi
 };

trunk/source/processes/hadronic/models/parton_string/diffraction/include/G4DiffractiveHHScatterer.hh

-                      r819
+                      r962
 //
 //
 // $Id: G4DiffractiveHHScatterer.hh,v 1.3 2006/06/29 20:54:25 gunter Exp $
+// $Id: G4DiffractiveHHScatterer.hh,v 1.4 2008/04/25 14:20:13 vuzhinsk Exp $
 #ifndef G4DiffractiveHHScatterer_h
 …
 class G4KineticTrack;
 #include "G4KineticTrackVector.hh"
+#include "G4FTFParameters.hh"                            // Uzhi 21.04.08
 class G4DiffractiveHHScatterer
 …
 const G4DiffractiveExcitation * theExcitation;
 G4LundStringFragmentation * theStringFragmentation;
+G4FTFParameters  *theParameters;                       // Uzhi  21.04.08
 };

trunk/source/processes/hadronic/models/parton_string/diffraction/include/G4DiffractiveSplitableHadron.hh

r819	r962
26	26	//
27	27	// $Id: G4DiffractiveSplitableHadron.hh,v 1.4 2006/06/29 20:54:27 gunter Exp $
28		// GEANT4 tag $Name: $
	28	// GEANT4 tag $Name: geant4-09-02-ref-02 $
29	29	//
30	30

trunk/source/processes/hadronic/models/parton_string/diffraction/include/G4FTFModel.hh

-                      r819
+                      r962
 //
 //
 // $Id: G4FTFModel.hh,v 1.5 2007/04/24 10:32:59 gunter Exp $
 // GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4FTFModel.hh,v 1.7 2008/04/25 14:20:13 vuzhinsk Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // Class Description
 …
 class G4VSplitableHadron;
 class G4ExcitedString;
+#include "G4FTFParameters.hh"                            // Uzhi 29.03.08
 #include "G4FTFParticipants.hh"
 #include "G4ExcitedStringVector.hh"
 #include "G4DiffractiveExcitation.hh"
+#include "G4ElasticHNScattering.hh"
 class G4FTFModel : public G4VPartonStringModel
 …
   private:
+       G4ReactionProduct theProjectile;
        G4FTFParticipants theParticipants;
+       G4ReactionProduct theProjectile;
+       G4FTFParameters  *theParameters;        // Uzhi  29.03.08
        G4DiffractiveExcitation * theExcitation;
+       G4ElasticHNScattering   * theElastic;   // Uzhi 29.03.08
 };
+// ------------------------------------------------------------
 inline
 G4V3DNucleus * G4FTFModel::GetWoundedNucleus() const

trunk/source/processes/hadronic/models/parton_string/diffraction/include/G4FTFParticipants.hh

-                      r819
+                      r962
 //
 //
 // $Id: G4FTFParticipants.hh,v 1.4 2006/06/29 20:54:32 gunter Exp $
 // GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4FTFParticipants.hh,v 1.5 2008/03/31 15:34:01 vuzhinsk Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 …
 #include "G4VParticipants.hh"
+#include "G4FTFParameters.hh"   // Uzhi 29.03.08
 #include <vector>
 #include "G4Nucleon.hh"
 …
       int operator!=(const G4FTFParticipants &right) const;
+      void BuildInteractions(const G4ReactionProduct  &thePrimary);
+      void GetList(const G4ReactionProduct  &thePrimary,
+                         G4FTFParameters    *theParameters); // Uzhi 29.03.08
+      void StartLoop();
       G4bool Next();
       const G4InteractionContent & GetInteraction() const;
-      void StartLoop();
+      std::vector<G4InteractionContent *> theInteractions;
   private:
       std::vector<G4InteractionContent *> theInteractions;
+//      std::vector<G4InteractionContent *> theInteractions;
       G4int currentInteraction;
 …
 #endif

trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4DiffractiveExcitation.cc

-                      r819
+                      r962
 //
 //
 // $Id: G4DiffractiveExcitation.cc,v 1.1 2007/05/25 06:56:53 vuzhinsk Exp $
+// $Id: G4DiffractiveExcitation.cc,v 1.7 2008/12/18 13:01:58 gunter Exp $
 // ------------------------------------------------------------
 //      GEANT 4 class implemetation file
 …
 //      ---------------- G4DiffractiveExcitation --------------
 //             by Gunter Folger, October 1998.
 //         diffractive Excitation used by strings models
 //         Take a projectile and a target
 //         excite the projectile and target
+//        diffractive Excitation used by strings models
+//             Take a projectile and a target
+//             excite the projectile and target
 //  Essential changed by V. Uzhinsky in November - December 2006
 //  in order to put it in a correspondence with original FRITIOF
 //  model. Variant of FRITIOF with nucleon de-excitation is implemented.
 //  Other changes by V.Uzhinsky in May 2007 were introduced to fit
+//  meson-nucleon interactions
+//  meson-nucleon interactions. Additional changes by V. Uzhinsky
+//  were introduced in December 2006. They treat diffraction dissociation
+//  processes more exactly.
 // ---------------------------------------------------------------------
 …
 #include "G4VSplitableHadron.hh"
 #include "G4ExcitedString.hh"
+#include "G4FTFParameters.hh"                            // Uzhi 19.04.08
 //#include "G4ios.hh"
+G4DiffractiveExcitation::G4DiffractiveExcitation()                         // Uzhi
+{
+}
+//#include "UZHI_diffraction.hh"
+G4DiffractiveExcitation::G4DiffractiveExcitation()
+{
+}
+// ---------------------------------------------------------------------
 G4bool G4DiffractiveExcitation::
+  ExciteParticipants(G4VSplitableHadron *projectile, G4VSplitableHadron *target) const
+{
+           G4LorentzVector Pprojectile=projectile->Get4Momentum();
+// -------------------- Projectile parameters -----------------------------------
+           G4bool PutOnMassShell=0;
+// With de-excitation
+//         G4double M0projectile=projectile->GetDefinition()->GetPDGMass();
+// Without de-excitation
+           G4double M0projectile = Pprojectile.mag();
+           if(M0projectile < projectile->GetDefinition()->GetPDGMass())
+             {
+              PutOnMassShell=1;
+              M0projectile=projectile->GetDefinition()->GetPDGMass();
+             }
+           G4double Mprojectile2 = M0projectile * M0projectile;
+           G4int    PDGcode=projectile->GetDefinition()->GetPDGEncoding();
+           G4int    absPDGcode=std::abs(PDGcode);
+           G4double ProjectileDiffCut;
+           G4double ProjectileBackgroundWeight;           // Uzhi May 2007
+           G4double TargetBackgroundWeight;               // Uzhi May 2007
+           G4double AveragePt2;
+           if( absPDGcode > 1000 )                        //------Projectile is baryon --------
+             {
+              ProjectileDiffCut = 1.1;                    // GeV
+              ProjectileBackgroundWeight=0.;              // Uzhi May 2007
+              AveragePt2 = 0.3;                           // GeV^2
+             }
+           else if( absPDGcode == 211 || PDGcode ==  111) //------Projectile is Pion -----------
+             {
+              ProjectileDiffCut = 0.6;                   // GeV Uzhi May 2007 1.0->0.6
+              ProjectileBackgroundWeight=0.9;            // Uzhi May 2007
+              AveragePt2 = 0.3;                          // GeV^2
+             }
+           else if( absPDGcode == 321 || PDGcode == -311) //------Projectile is Kaon -----------
+             {
+              ProjectileDiffCut = 0.7;                    // GeV 1.1
+              ProjectileBackgroundWeight=0.5;             // Uzhi May 2007   ???
+              AveragePt2 = 0.3;                           // GeV^2
+             }
+           else                                           //------Projectile is undefined,
+                                                          //------Nucleon assumed
+             {
+              ProjectileDiffCut = 1.1;                    // GeV
+              ProjectileBackgroundWeight=0.;              // Uzhi May 2007
+              AveragePt2 = 0.3;                           // GeV^2
+             };
+           ProjectileDiffCut = ProjectileDiffCut * GeV;
+           AveragePt2 = AveragePt2 * GeV*GeV;
+// -------------------- Target parameters ----------------------------------------------
+           G4LorentzVector Ptarget=target->Get4Momentum();
+           G4double M0target = Ptarget.mag();
+           if(M0target < target->GetDefinition()->GetPDGMass())
+             {
+              PutOnMassShell=1;
+              M0target=target->GetDefinition()->GetPDGMass();
+             }
+           G4double Mtarget2 = M0target * M0target;                      //Ptarget.mag2();
+                                                                         // for AA-inter.
+           G4double NuclearNucleonDiffCut = 1.1*GeV;
+           TargetBackgroundWeight=0.;                           // Uzhi May 2007
+           G4double ProjectileDiffCut2 = ProjectileDiffCut * ProjectileDiffCut;
+           G4double NuclearNucleonDiffCut2 = NuclearNucleonDiffCut * NuclearNucleonDiffCut;
+  ExciteParticipants(G4VSplitableHadron *projectile,
+                     G4VSplitableHadron *target,
+                     G4FTFParameters    *theParameters) const
+{
+     G4bool PutOnMassShell=0;
+// -------------------- Projectile parameters -----------------------
+     G4LorentzVector Pprojectile=projectile->Get4Momentum();
+//   G4double M0projectile=projectile->GetDefinition()->GetPDGMass(); // With de-excitation
+     G4double M0projectile = Pprojectile.mag();                       // Without de-excitation
+/*
+G4cout<<"ExciteParticipants-------------------"<<G4endl;
+G4cout<<"Mom "<<Pprojectile<<" mass "<<M0projectile<<G4endl;
+*/
+     if(M0projectile < projectile->GetDefinition()->GetPDGMass())
+     {
+      PutOnMassShell=1;
+      M0projectile=projectile->GetDefinition()->GetPDGMass();
+     }
+     G4double M0projectile2 = M0projectile * M0projectile;
+     G4int    PDGcode=projectile->GetDefinition()->GetPDGEncoding();
+     G4int    absPDGcode=std::abs(PDGcode);
+     G4double ProjectileDiffStateMinMass=theParameters->GetProjMinDiffMass();
+     G4double ProjectileNonDiffStateMinMass=theParameters->GetProjMinNonDiffMass();
+     G4double ProbProjectileDiffraction=theParameters->GetProbabilityOfProjDiff();
+/*
+G4cout<<ProjectileDiffStateMinMass<<" "<<ProjectileNonDiffStateMinMass<<" "<<ProbProjectileDiffraction<<G4endl;
+*/
+// -------------------- Target paraExciteParticipantsmeters -------------------------
+     G4LorentzVector Ptarget=target->Get4Momentum();
+     G4double M0target = Ptarget.mag();
+//G4cout<<"Mom "<<Ptarget<<" mass "<<M0target<<G4endl;
+     if(M0target < target->GetDefinition()->GetPDGMass())
+     {
+      PutOnMassShell=1;
+      M0target=target->GetDefinition()->GetPDGMass();
+     }
+     G4double M0target2 = M0target * M0target;             //Ptarget.mag2();
+                                                           // for AA-inter.
+     G4double TargetDiffStateMinMass=theParameters->GetTarMinDiffMass();
+     G4double TargetNonDiffStateMinMass=theParameters->GetTarMinNonDiffMass();
+     G4double ProbTargetDiffraction=theParameters->GetProbabilityOfTarDiff();
+/*
+G4cout<<TargetDiffStateMinMass<<" "<<TargetNonDiffStateMinMass<<" "<<ProbTargetDiffraction<<G4endl;
+*/
+     G4double AveragePt2=theParameters->GetAveragePt2();
+// Kinematical properties of the interactions --------------
+     G4LorentzVector Psum;      // 4-momentum in CMS
+     Psum=Pprojectile+Ptarget;
+     G4double S=Psum.mag2();
+//G4cout<<" sqrt(s) "<<std::sqrt(S)<<G4endl;
+// ------------------------------------------------------------------
+//ProbProjectileDiffraction=1.;
+//ProbTargetDiffraction    =1.;
+     G4double ProbOfDiffraction=ProbProjectileDiffraction +
+                                ProbTargetDiffraction;
+     if(ProbOfDiffraction!=0.)
+     {
+      ProbProjectileDiffraction/=ProbOfDiffraction;
+     }
+     else
+     {
+      ProbProjectileDiffraction=0.;
+     }
+//   ProbTargetDiffraction    /=ProbOfDiffraction;
+//G4cout<<"ProbOfDiffraction "<<ProbOfDiffraction<<"ProbProjectileDiffraction "<<ProbProjectileDiffraction<<G4endl;   // Vova
+     G4double ProjectileDiffStateMinMass2    = ProjectileDiffStateMinMass    *
+                                               ProjectileDiffStateMinMass;
+     G4double ProjectileNonDiffStateMinMass2 = ProjectileNonDiffStateMinMass *
+                                               ProjectileNonDiffStateMinMass;
+     G4double TargetDiffStateMinMass2        = TargetDiffStateMinMass        *
+                                               TargetDiffStateMinMass;
+     G4double TargetNonDiffStateMinMass2     = TargetNonDiffStateMinMass     *
+                                               TargetNonDiffStateMinMass;
 // Transform momenta to cms and then rotate parallel to z axis;
+           G4LorentzVector Psum;
+           Psum=Pprojectile+Ptarget;
+           G4LorentzRotation toCms(-1*Psum.boostVector());
+           G4LorentzVector Ptmp=toCms*Pprojectile;
+           if ( Ptmp.pz() <= 0. )
+           {
+//         G4LorentzVector Psum;
+//         Psum=Pprojectile+Ptarget;
+     G4LorentzRotation toCms(-1*Psum.boostVector());
+     G4LorentzVector Ptmp=toCms*Pprojectile;
+     if ( Ptmp.pz() <= 0. )
+        {
            // "String" moving backwards in  CMS, abort collision !!
            //G4cout << " abort Collision!! " << G4endl;
+                   return false;
+           }
+           toCms.rotateZ(-1*Ptmp.phi());
+           toCms.rotateY(-1*Ptmp.theta());
+           G4LorentzRotation toLab(toCms.inverse());
+           Pprojectile.transform(toCms);
+           Ptarget.transform(toCms);
+           G4double Pt2;
+           G4double ProjMassT2, ProjMassT;
+           G4double TargMassT2, TargMassT;
+           G4double PZcms2, PZcms;
+           G4double PMinusNew, TPlusNew;
+           G4double S=Psum.mag2();
+           G4double SqrtS=std::sqrt(S);
+           if(absPDGcode > 1000 && SqrtS < 2200*MeV)
+             {return false;}                         // The model cannot work for
+                                                     // p+p-interactions
+                                                     // at Plab < 1.3 GeV/c.
+           if(( absPDGcode == 211 || PDGcode ==  111) && SqrtS < 1600*MeV)
+             {return false;}                         // The model cannot work for
+                                                     // Pi+p-interactions
+                                                     // at Plab < 1. GeV/c.
+           if(( absPDGcode == 321 || PDGcode == -311) && SqrtS < 1600*MeV)
+             {return false;}                         // The model cannot work for
+                                                     // K+p-interactions
+                                                     // at Plab < ??? GeV/c.  ???
+           PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2-
+*S*Mprojectile2-2*S*Mtarget2-2*Mprojectile2*Mtarget2)/4./S;
+           if(PZcms2 < 0)
+             {return false;}   // It can be in an interaction with off-shell nuclear nucleon
+           PZcms = std::sqrt(PZcms2);
+           if(PutOnMassShell)
+          return false;
+         }
+     toCms.rotateZ(-1*Ptmp.phi());
+     toCms.rotateY(-1*Ptmp.theta());
+     G4LorentzRotation toLab(toCms.inverse());
+     Pprojectile.transform(toCms);
+     Ptarget.transform(toCms);
+     G4double Pt2;
+     G4double ProjMassT2, ProjMassT;
+     G4double TargMassT2, TargMassT;
+     G4double PZcms2, PZcms;
+     G4double PMinusMin, PMinusMax;
+//     G4double PPlusMin , PPlusMax;
+     G4double TPlusMin , TPlusMax;
+     G4double PMinusNew, PPlusNew, TPlusNew, TMinusNew;
+//   G4double S=Psum.mag2();
+     G4double SqrtS=std::sqrt(S);
+     if(absPDGcode > 1000 && SqrtS < 2200*MeV)
+     {return false;}                         // The model cannot work for
+                                             // p+p-interactions
+                                             // at Plab < 1.3 GeV/c.
+     if(( absPDGcode == 211 || PDGcode ==  111) && SqrtS < 1600*MeV)
+     {return false;}                         // The model cannot work for
+                                             // Pi+p-interactions
+                                             // at Plab < 1. GeV/c.
+     if(( absPDGcode == 321 || PDGcode == -311) && SqrtS < 1600*MeV)
+     {return false;}                         // The model cannot work for
+                                             // K+p-interactions
+                                             // at Plab < ??? GeV/c.  ???
+     PZcms2=(S*S+M0projectile2*M0projectile2+M0target2*M0target2-
+*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2)
+             /4./S;
+     if(PZcms2 < 0)
+     {return false;}   // It can be in an interaction with off-shell nuclear nucleon
+     PZcms = std::sqrt(PZcms2);
+     if(PutOnMassShell)
+     {
+      if(Pprojectile.z() > 0.)
+      {
+       Pprojectile.setPz( PZcms);
+       Ptarget.setPz(    -PZcms);
+      }
+      else
+      {
+       Pprojectile.setPz(-PZcms);
+       Ptarget.setPz(     PZcms);
+      };
+      Pprojectile.setE(std::sqrt(M0projectile2                  +
+                                 Pprojectile.x()*Pprojectile.x()+
+                                 Pprojectile.y()*Pprojectile.y()+
+                                 PZcms2));
+      Ptarget.setE(std::sqrt(M0target2              +
+                             Ptarget.x()*Ptarget.x()+
+                             Ptarget.y()*Ptarget.y()+
+                             PZcms2));
+     }
+     G4double maxPtSquare; // = PZcms2;
+/*
+G4cout << "Pprojectile aft boost : " << Pprojectile <<" "<<Pprojectile.mag()<< G4endl;
+G4cout << "Ptarget aft boost : " << Ptarget <<" "<<Ptarget.mag()<< G4endl;
+G4cout << "cms aft boost : " << (Pprojectile+ Ptarget) << G4endl;
+G4cout << " Projectile Xplus / Xminus : " <<
+            Pprojectile.plus() << " / " << Pprojectile.minus() << G4endl;
+G4cout << " Target Xplus / Xminus : " <<           Ptarget.plus() << " / " << Ptarget.minus() << G4endl;
+G4cout<<"maxPtSquare "<<maxPtSquare<<G4endl;
+*/
+     G4LorentzVector Qmomentum;
+     G4double Qminus, Qplus;
+     G4int whilecount=0;
+//  Choose a process
+     if(G4UniformRand() < ProbOfDiffraction)
+       {
+        if(G4UniformRand() < ProbProjectileDiffraction)
+        { //-------- projectile diffraction ---------------
+//G4cout<<"   Projectile diffraction"<<G4endl;
+//Uzhi_projectilediffraction++;
+         do {
+//             Generate pt
+//             if (whilecount++ >= 500 && (whilecount%100)==0)
+//               G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping"
+//               << ", loop count/ maxPtSquare : "
+//               << whilecount << " / " << maxPtSquare << G4endl;
+             if (whilecount > 1000 )
+             {
+              if(Pprojectile.z() > 0.)
+                {
+                 Pprojectile.setPz( PZcms);
+                 Ptarget.setPz(    -PZcms);
+                }
+              else
+                 {
+                 Pprojectile.setPz(-PZcms);
+                 Ptarget.setPz(     PZcms);
+                };
+               Pprojectile.setE(std::sqrt(Mprojectile2+
+                                                       Pprojectile.x()*Pprojectile.x()+
+                                                       Pprojectile.y()*Pprojectile.y()+
+                                                       PZcms2));
+               Ptarget.setE(std::sqrt(    Mtarget2    +
+                                                       Ptarget.x()*Ptarget.x()+
+                                                       Ptarget.y()*Ptarget.y()+
+                                                       PZcms2));
+             }
+           G4double maxPtSquare = PZcms2;
+//G4cout << "Pprojectile aft boost : " << Pprojectile << G4endl;
+//G4cout << "Ptarget aft boost : " << Ptarget << G4endl;
+//         G4cout << "cms aft boost : " << (Pprojectile+ Ptarget) << G4endl;
+//         G4cout << " Projectile Xplus / Xminus : " <<
+//              Pprojectile.plus() << " / " << Pprojectile.minus() << G4endl;
+//         G4cout << " Target Xplus / Xminus : " <<
+//              Ptarget.plus() << " / " << Ptarget.minus() << G4endl;
+           G4LorentzVector Qmomentum;
+           G4double Qminus, Qplus;
+           G4int whilecount=0;
+           do {
+//  Generate pt
+               if (whilecount++ >= 500 && (whilecount%100)==0)
+              Qmomentum=G4LorentzVector(0.,0.,0.,0.);
+              return false;    //  Ignore this interaction
+             };
+// --------------- Check that the interaction is possible -----------
+             ProjMassT2=ProjectileDiffStateMinMass2;
+             ProjMassT =ProjectileDiffStateMinMass;
+             TargMassT2=M0target2;
+             TargMassT =M0target;
+             PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
+.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
+                    /4./S;
+//G4cout<<" Pt2 Mpt Mtt Pz2 "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl;
+if(PZcms2 < 0 )
+{
+/*
+G4cout<<"whilecount "<<whilecount<<" "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl;
+G4int Uzhi; G4cin>>Uzhi;
+*/
+return false;
+};
+             maxPtSquare=PZcms2;
+             Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
+             Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
+             ProjMassT2=ProjectileDiffStateMinMass2+Pt2;
+             ProjMassT =std::sqrt(ProjMassT2);
+             TargMassT2=M0target2+Pt2;
+             TargMassT =std::sqrt(TargMassT2);
+             PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
+.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
+                    /4./S;
+//G4cout<<" Pt2 Mpt Mtt Pz2 "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl;
+//             if(PZcms2 < 0 ) {PZcms2=0;};
+             if(PZcms2 < 0 ) continue;
+             PZcms =std::sqrt(PZcms2);
+             PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
+             PMinusMax=SqrtS-TargMassT;
+//G4cout<<" SqrtS P+mim max "<<SqrtS<<" "<<PMinusMin<<" "<<PMinusMax<<G4endl;
+             PMinusNew=ChooseP(PMinusMin, PMinusMax);
+// PMinusNew=1./sqrt(1./PMinusMin-G4UniformRand()*(1./PMinusMin-1./PMinusMax));
+             TMinusNew=SqrtS-PMinusNew;
+             Qminus=Ptarget.minus()-TMinusNew;
+             TPlusNew=TargMassT2/TMinusNew;
+             Qplus=Ptarget.plus()-TPlusNew;
+             Qmomentum.setPz( (Qplus-Qminus)/2 );
+             Qmomentum.setE(  (Qplus+Qminus)/2 );
+          } while (
+((Pprojectile+Qmomentum).mag2() <  ProjectileDiffStateMinMass2) ||  //No without excitation
+((Ptarget    -Qmomentum).mag2() <  M0target2                  ));
+        }
+        else
+        { // -------------- Target diffraction ----------------
+//G4cout<<"   Target difraction"<<G4endl;
+//Uzhi_targetdiffraction++;
+         do {
+//             Generate pt
+//             if (whilecount++ >= 500 && (whilecount%100)==0)
 //               G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping"
 //               << ", loop count/ maxPtSquare : "
+//               << ", loop count/ maxPtSquare : "
 //               << whilecount << " / " << maxPtSquare << G4endl;
+               if (whilecount > 1000 )
+               {
+                 Qmomentum=G4LorentzVector(0.,0.,0.,0.);
+                 return false;    //  Ignore this interaction
+               }
+               Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
+//G4cout << "generated Pt " << Qmomentum << G4endl;
+//G4cout << "Pprojectile with pt : " << Pprojectile+Qmomentum << G4endl;
+//G4cout << "Ptarget with pt : " << Ptarget-Qmomentum << G4endl;
+//  Momentum transfer
+/*                                                                          // Uzhi
+               G4double Xmin = minmass / ( Pprojectile.e() + Ptarget.e() );
+               G4double Xmax=1.;
+               G4double Xplus =ChooseX(Xmin,Xmax);
+               G4double Xminus=ChooseX(Xmin,Xmax);
+//             G4cout << " X-plus  " << Xplus << G4endl;
+//             G4cout << " X-minus " << Xminus << G4endl;
+               G4double pt2=G4ThreeVector(Qmomentum.vect()).mag2();
+               G4double Qplus =-1 * pt2 / Xminus/Ptarget.minus();
+               G4double Qminus=     pt2 / Xplus /Pprojectile.plus();
+*/                                                                          // Uzhi    *
+               Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
+               ProjMassT2=Mprojectile2+Pt2;
+               ProjMassT =std::sqrt(ProjMassT2);
+               TargMassT2=Mtarget2+Pt2;
+               TargMassT =std::sqrt(TargMassT2);
+               PZcms2=(S*S+ProjMassT2*ProjMassT2+
+                           TargMassT2*TargMassT2-
+.*S*ProjMassT2-2.*S*TargMassT2-
+.*ProjMassT2*TargMassT2)/4./S;
+               if(PZcms2 < 0 ) {PZcms2=0;};
+               PZcms =std::sqrt(PZcms2);
+               G4double PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
+               G4double PMinusMax=SqrtS-TargMassT;
+               if(G4UniformRand() < ProjectileBackgroundWeight)              // Uzhi May 2007
+                 {
+                  PMinusNew=PMinusMax-(PMinusMax-PMinusMin)*G4UniformRand(); // Uzhi May 2007
+                 }                                                           // Uzhi May 2007
+               else                                                          // Uzhi May 2007
+                 {                                                           // Uzhi May 2007
+                  PMinusNew=ChooseP(PMinusMin, PMinusMax);                   // Uzhi May 2007
+                 };                                                          // Uzhi May 2007
+//               PMinusNew=ChooseP(PMinusMin, PMinusMax);                    // Uzhi May 2007
+               Qminus=PMinusNew-Pprojectile.minus();
+               G4double TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
+               G4double TPlusMax=SqrtS-PMinusNew;                            // Uzhi May 2007
+//               G4double TPlusMax=SqrtS-ProjMassT;                          // Uzhi May 2007
+               if(G4UniformRand() < TargetBackgroundWeight)                  // Uzhi May 2007
+                 {
+                  TPlusNew=TPlusMax-(TPlusMax-TPlusMin)*G4UniformRand();     // Uzhi May 2007
+                 }                                                           // Uzhi May 2007
+               else                                                          // Uzhi May 2007
+                 {                                                           // Uzhi May 2007
+                  TPlusNew=ChooseP(TPlusMin, TPlusMax);                      // Uzhi May 2007
+                 };                                                          // Uzhi May 2007
+//               TPlusNew=ChooseP(TPlusMin, TPlusMax);                       // Uzhi May 2007
+               Qplus=-(TPlusNew-Ptarget.plus());
+               Qmomentum.setPz( (Qplus-Qminus)/2 );
+               Qmomentum.setE(  (Qplus+Qminus)/2 );
+//G4cout << "Qplus / Qminus " << Qplus << " / " << Qminus<<G4endl;
+//           G4cout << "pt2" << pt2 << G4endl;
+//           G4cout << "Qmomentum " << Qmomentum << G4endl;
+//           G4cout << " Masses (P/T) : " << (Pprojectile+Qmomentum).mag() <<
+//                             " / " << (Ptarget-Qmomentum).mag() << G4endl;
+           } while (
+( (Pprojectile+Qmomentum).mag2() <  Mprojectile2       ||      // Uzhi No without excitation
+  (Ptarget    -Qmomentum).mag2() <  Mtarget2             ) ||  // Uzhi
+( (Pprojectile+Qmomentum).mag2() <  ProjectileDiffCut2 &&      // Uzhi No double Diffraction
+  (Ptarget    -Qmomentum).mag2()  <  NuclearNucleonDiffCut2) );// Uzhi
+           if((Ptarget-Qmomentum).mag2() < NuclearNucleonDiffCut2)           // Uzhi
+//Projectile diffraction
+             if (whilecount > 1000 )
+             {
+              G4double TMinusNew=SqrtS-PMinusNew;
+              Qminus=Ptarget.minus()-TMinusNew;
+              TPlusNew=TargMassT2/TMinusNew;
+              Qplus=Ptarget.plus()-TPlusNew;
+              Qmomentum.setPz( (Qplus-Qminus)/2 );
+              Qmomentum.setE(  (Qplus+Qminus)/2 );
+             }
+           else if((Pprojectile+Qmomentum).mag2() <  ProjectileDiffCut2)    // Uzhi
+//Target diffraction
+              Qmomentum=G4LorentzVector(0.,0.,0.,0.);
+              return false;    //  Ignore this interaction
+             };
+// --------------- Check that the interaction is possible -----------
+             ProjMassT2=M0projectile2;
+             ProjMassT =M0projectile;
+             TargMassT2=TargetDiffStateMinMass2;
+             TargMassT =TargetDiffStateMinMass;
+             PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
+.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
+                    /4./S;
+//G4cout<<" Pt2 Mpt Mtt Pz2 "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl;
+if(PZcms2 < 0 )
+{
+/*
+G4cout<<"whilecount "<<whilecount<<" "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl;
+G4int Uzhi; G4cin>>Uzhi;
+*/
+return false;
+};
+             maxPtSquare=PZcms2;
+             Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
+             Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
+             ProjMassT2=M0projectile2+Pt2;
+             ProjMassT =std::sqrt(ProjMassT2);
+             TargMassT2=TargetDiffStateMinMass2+Pt2;
+             TargMassT =std::sqrt(TargMassT2);
+             PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
+.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
+                    /4./S;
+/*
+if(PZcms2 < 0 )
+{
+G4cout<<"whilecount "<<whilecount<<" "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl;
+G4int Uzhi; G4cin>>Uzhi;
+return false;
+};
+*/
+             if(PZcms2 < 0 ) continue;
+             PZcms =std::sqrt(PZcms2);
+             TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
+             TPlusMax=SqrtS-ProjMassT;
+//G4cout<<" Tmin max "<<TPlusMin<<" "<<TPlusMax<<G4endl;
+             TPlusNew=ChooseP(TPlusMin, TPlusMax);
+//TPlusNew=TPlusMax;
+//G4cout<<"T+new "<<TPlusNew<<G4endl;
+             PPlusNew=SqrtS-TPlusNew;
+             Qplus=PPlusNew-Pprojectile.plus();
+             PMinusNew=ProjMassT2/PPlusNew;
+             Qminus=PMinusNew-Pprojectile.minus();
+             Qmomentum.setPz( (Qplus-Qminus)/2 );
+             Qmomentum.setE(  (Qplus+Qminus)/2 );
+          } while (
+ ((Pprojectile+Qmomentum).mag2() <  M0projectile2          ) ||  //No without excitation
+ ((Ptarget    -Qmomentum).mag2() <  TargetDiffStateMinMass2));
+         }
+        }
+        else  //----------- Non-diffraction process ------------
+        {
+//G4cout<<"   Non-difraction"<<G4endl;
+         do {
+//             Generate pt
+//             if (whilecount++ >= 500 && (whilecount%100)==0)
+//               G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping"
+//               << ", loop count/ maxPtSquare : "
+//               << whilecount << " / " << maxPtSquare << G4endl;
+             if (whilecount > 1000 )
+             {
+              G4double PPlusNew=SqrtS-TPlusNew;
+              Qplus=PPlusNew-Pprojectile.plus();
+              PMinusNew=ProjMassT2/PPlusNew;
+              Qminus=PMinusNew-Pprojectile.minus();
+              Qmomentum.setPz( (Qplus-Qminus)/2 );
+              Qmomentum.setE(  (Qplus+Qminus)/2 );
+              Qmomentum=G4LorentzVector(0.,0.,0.,0.);
+              return false;    //  Ignore this interaction
              };
+// --------------- Check that the interaction is possible -----------
+             ProjMassT2=ProjectileNonDiffStateMinMass2;
+             ProjMassT =ProjectileNonDiffStateMinMass;
+             TargMassT2=TargetNonDiffStateMinMass2;
+             TargMassT =TargetNonDiffStateMinMass;
+             PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
+.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
+                   /4./S;
+//G4cout<<" Pt2 Mpt Mtt Pz2 "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl;
+if(PZcms2 < 0 )
+{
+/*
+G4cout<<"whilecount "<<whilecount<<" "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl;
+G4int Uzhi; G4cin>>Uzhi;
+*/
+return false;
+};
+             maxPtSquare=PZcms2;
+             Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
+             Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
+             ProjMassT2=ProjectileNonDiffStateMinMass2+Pt2;
+             ProjMassT =std::sqrt(ProjMassT2);
+             TargMassT2=TargetNonDiffStateMinMass2+Pt2;
+             TargMassT =std::sqrt(TargMassT2);
+             PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
+.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
+                   /4./S;
+/*
+G4cout<<"ProjectileNonDiffStateMinMass2 "<<ProjectileNonDiffStateMinMass2<<G4endl;
+G4cout<<"TargetNonDiffStateMinMass2     "<<TargetNonDiffStateMinMass2<<G4endl;
+G4cout<<"Mt "<<ProjMassT<<" "<<TargMassT<<" "<<Pt2<<" "<<PZcms2<<G4endl<<G4endl;
+*/
+//             if(PZcms2 < 0 ) {PZcms2=0;};
+             if(PZcms2 < 0 ) continue;
+             PZcms =std::sqrt(PZcms2);
+             PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
+             PMinusMax=SqrtS-TargMassT;
+             PMinusNew=ChooseP(PMinusMin, PMinusMax);
+// PMinusNew=1./sqrt(1./PMinusMin-G4UniformRand()*(1./PMinusMin-1./PMinusMax));
+//G4cout<<"Proj "<<PMinusMin<<" "<<PMinusMax<<" "<<PMinusNew<<G4endl;
+//PMinusNew=PMinusMax; //+++++++++++++++++++++++++++++++++++ Vova
+             Qminus=PMinusNew-Pprojectile.minus();
+             TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
+//             TPlusMax=SqrtS-PMinusNew;                      // Vova
+             TPlusMax=SqrtS-ProjMassT;                        // Vova
+             TPlusNew=ChooseP(TPlusMin, TPlusMax);
+//G4cout<<"Targ "<<TPlusMin<<" "<<TPlusMax<<" "<<TPlusNew<<G4endl;
+//G4cout<<PMinusNew<<" "<<TPlusNew<<G4endl;
+             Qplus=-(TPlusNew-Ptarget.plus());
+             Qmomentum.setPz( (Qplus-Qminus)/2 );
+             Qmomentum.setE(  (Qplus+Qminus)/2 );
+/*
+G4cout << "Qplus / Qminus " << Qplus << " / " << Qminus<<G4endl;
+G4cout << "pt2" << pt2 << G4endl;
+G4cout << "Qmomentum " << Qmomentum << G4endl;
+G4cout << " Masses (P/T) : " << (Pprojectile+Qmomentum).mag() <<
+           " / " << (Ptarget-Qmomentum).mag() << G4endl;   // mag()
+G4cout<<"Mprojectile "<<std::sqrt(M0projectile2)<<G4endl;
+G4cout<<"Mtarget     "<<std::sqrt(M0target2    )<<G4endl;
+G4cout<<"ProjectileDiffStateMinMass "<<std::sqrt(ProjectileDiffStateMinMass2)<<G4endl;
+G4cout<<"TargetDiffStateMinMass "<<std::sqrt(TargetDiffStateMinMass2)<<G4endl;
+*/
+       } while (
+ ((Pprojectile+Qmomentum).mag2() <  ProjectileNonDiffStateMinMass2) || //No double Diffraction
+ ((Ptarget    -Qmomentum).mag2() <  TargetNonDiffStateMinMass2    ));
+         }
+//G4int Uzhiinp; G4cin>>Uzhiinp;    // Vova
            Pprojectile += Qmomentum;
            Ptarget     -= Qmomentum;
+/*
+Pprojectile.setPz(0.);
+Pprojectile.setE(SqrtS-M0target);
+Ptarget.setPz(0.);
+Ptarget.setE(M0target);
+*/
+//G4cout << "Pprojectile with Q : " << Pprojectile << G4endl;
+//G4cout << "Ptarget with Q : " << Ptarget << G4endl;
+//         G4cout << "Projectile back: " << toLab * Pprojectile << G4endl;
+//         G4cout << "Target back: " << toLab * Ptarget << G4endl;
+/*
+G4cout << "Pprojectile with Q : " << Pprojectile << G4endl;
+G4cout << "Ptarget with Q : " << Ptarget << G4endl;
+G4cout << "Target        mass  " <<  Ptarget.mag() << G4endl;
+G4cout << "Projectile mass  " <<  Pprojectile.mag() << G4endl;
+//
+//G4cout << "Projectile back: " << toLab * Pprojectile << G4endl;
+//G4cout << "Target back: " << toLab * Ptarget << G4endl;
+*/
+//-------------- Flip if projectale moves in backward direction ------------
+//G4bool Flip=Pprojectile.pz()< 0.;
 // Transform back and update SplitableHadron Participant.
 …
 //G4cout << "Pprojectile with Q M: " << Pprojectile<<" "<<  Pprojectile.mag() << G4endl;
 //G4cout << "Ptarget     with Q M: " << Ptarget    <<" "<<  Ptarget.mag()     << G4endl;
+//G4cout << "Target      mass  " <<  Ptarget.mag() << G4endl;
+//G4cout << "Target      mass  " <<  Ptarget.mag() << G4endl;
+//G4cout << "Projectile mass  " <<  Pprojectile.mag() << G4endl;
+/*
+if(!Flip){
+           projectile->Set4Momentum(Pprojectile);
            target->Set4Momentum(Ptarget);
+//G4cout << "Projectile mass  " <<  Pprojectile.mag() << G4endl;
+         }
+else     {
+           G4ParticleDefinition * t_Definition=projectile->GetDefinition();
+           projectile->SetDefinition(target->GetDefinition());
+           projectile->Set4Momentum(Ptarget);
+           target->SetDefinition(t_Definition);
+           target->Set4Momentum(Pprojectile);
+          }
+*/
+//
+/*
+if(G4UniformRand() < 1.) {
+           G4ParticleDefinition * t_Definition=projectile->GetDefinition();
+           projectile->SetDefinition(target->GetDefinition());
+           target->SetDefinition(t_Definition);
+}
+*/ // For flip, for HARP
+           G4double ZcoordinateOfCurrentInteraction = target->GetPosition().z();
+// It is assumed that nucleon z-coordinates are ordered on increasing -----------
+           G4double betta_z=projectile->Get4Momentum().pz()/projectile->Get4Momentum().e();
+           G4double ZcoordinateOfPreviousCollision=projectile->GetPosition().z();
+           if(projectile->GetSoftCollisionCount()==0) {
+              projectile->SetTimeOfCreation(0.);
+              target->SetTimeOfCreation(0.);
+              ZcoordinateOfPreviousCollision=ZcoordinateOfCurrentInteraction;
+           }
+           G4ThreeVector thePosition(projectile->GetPosition().x(),
+                                     projectile->GetPosition().y(),
+                                     ZcoordinateOfCurrentInteraction);
+           projectile->SetPosition(thePosition);
+           G4double TimeOfPreviousCollision=projectile->GetTimeOfCreation();
+           G4double TimeOfCurrentCollision=TimeOfPreviousCollision+
+                    (ZcoordinateOfCurrentInteraction-ZcoordinateOfPreviousCollision)/betta_z;
+           projectile->SetTimeOfCreation(TimeOfCurrentCollision);
+           target->SetTimeOfCreation(TimeOfCurrentCollision);
            projectile->Set4Momentum(Pprojectile);
+           target->Set4Momentum(Ptarget);
+           projectile->IncrementCollisionCount(1);
+           target->IncrementCollisionCount(1);
+//
+//G4cout<<"Out of Excitation --------------------"<<G4endl;
+//G4int Uzhiinp; G4cin>>Uzhiinp;    // Vova
            return true;
+}
+// ---------------------------------------------------------------------
 G4ExcitedString * G4DiffractiveExcitation::
            String(G4VSplitableHadron * hadron, G4bool isProjectile) const
+{
+//G4cout<<"G4DiffractiveExcitation::String isProj"<<isProjectile<<G4endl;
         hadron->SplitUp();
         G4Parton *start= hadron->GetNextParton();
         if ( start==NULL)
+        if ( start==NULL)
         { G4cout << " G4FTFModel::String() Error:No start parton found"<< G4endl;
           return NULL;
+        }
         G4Parton *end  = hadron->GetNextParton();
         if ( end==NULL)
+        if ( end==NULL)
         { G4cout << " G4FTFModel::String() Error:No end parton found"<< G4endl;
           return NULL;
+        }
         G4ExcitedString * string;
         if ( isProjectile )
+        if ( isProjectile )
+        {
                 string= new G4ExcitedString(end,start, +1);
 …
                 string= new G4ExcitedString(start,end, -1);
+        }
+// Uzhi
+//G4cout<<"G4ExcitedString * G4DiffractiveExcitation::String"<<G4endl;
+//G4cout<<hadron->GetTimeOfCreation()<<" "<<hadron->GetPosition()/fermi<<G4endl;
+        string->SetTimeOfCreation(hadron->GetTimeOfCreation());
         string->SetPosition(hadron->GetPosition());
 // momenta of string ends
+//
+        G4double Momentum=hadron->Get4Momentum().vect().mag();
+        G4double  Plus=hadron->Get4Momentum().e() + Momentum;
+        G4double Minus=hadron->Get4Momentum().e() - Momentum;
+        G4ThreeVector tmp;
+        if(Momentum > 0.)
+        {
+         tmp.set(hadron->Get4Momentum().px(),
+                 hadron->Get4Momentum().py(),
+                 hadron->Get4Momentum().pz());
+         tmp/=Momentum;
+         }
+         else
+         {
+          tmp.set(0.,0.,1.);
+         };
+        G4LorentzVector Pstart(tmp,0.);
+        G4LorentzVector   Pend(tmp,0.);
+        if(isProjectile)
+        {
+         Pstart*=(-1.)*Minus/2.;
+         Pend  *=(+1.)*Plus /2.;
+         }
+        else
+         {
+          Pstart*=(+1.)*Plus/2.;
+          Pend  *=(-1.)*Minus/2.;
+         };
+         Momentum=-Pstart.mag();
+         Pstart.setT(Momentum);  // It is assumed that quark has m=0.
+         Momentum=-Pend.mag();
+         Pend.setT(Momentum);    // It is assumed that di-quark has m=0.
+//
+/* Uzhi
         G4double ptSquared= hadron->Get4Momentum().perp2();
         G4double transverseMassSquared= hadron->Get4Momentum().plus()
 …
                  sqr( std::sqrt(transverseMassSquared) - std::sqrt(ptSquared) );
         G4double widthOfPtSquare = 0.25;          // Uzhi <Pt^2>=0.25 ??????????????????
+        G4double widthOfPtSquare = 0.25*GeV*GeV;       // Uzhi 11.07 <Pt^2>=0.25 ??????????????????
         G4ThreeVector pt=GaussianPt(widthOfPtSquare,maxAvailMomentumSquared);
 …
         G4int Sign= isProjectile ? -1 : 1;
         G4double endMinus  = 0.5 * (tm1 + Sign*tm2);
         G4double startMinus= hadron->Get4Momentum().minus() - endMinus;
 …
         Pend.setPz(0.5*(endPlus - endMinus));
         Pend.setE(0.5*(endPlus + endMinus));
+*/ // Uzhi
         start->Set4Momentum(Pstart);
         end->Set4Momentum(Pend);
+/*
+G4cout<<"G4DiffractiveExcitation::String hadro"<<hadron->Get4Momentum()<<" "<<hadron->Get4Momentum().mag2()<<G4endl;
+G4cout<<"G4DiffractiveExcitation::String start"<<start->Get4Momentum()<<" "<<start->GetPDGcode()<<G4endl;
+G4cout<<"G4DiffractiveExcitation::String end  "<<  end->Get4Momentum()<<" "<<  end->GetPDGcode()<<G4endl;
+G4int Uzhi; G4cin>>Uzhi;
+*/
 #ifdef G4_FTFDEBUG
+        G4cout << " generated string flavors          " << start->GetPDGcode() << " / " << end->GetPDGcode() << G4endl;
+        G4cout << " generated string momenta:   quark " << start->Get4Momentum() << "mass : " <<start->Get4Momentum().mag()<< G4endl;
+        G4cout << " generated string momenta: Diquark " << end ->Get4Momentum() << "mass : " <<end->Get4Momentum().mag()<< G4endl;
+        G4cout << " generated string flavors          "
+               << start->GetPDGcode() << " / "
+               << end->GetPDGcode()                    << G4endl;
+        G4cout << " generated string momenta:   quark "
+               << start->Get4Momentum() << "mass : "
+               <<start->Get4Momentum().mag()           << G4endl;
+        G4cout << " generated string momenta: Diquark "
+               << end ->Get4Momentum()
+               << "mass : " <<end->Get4Momentum().mag()<< G4endl;
         G4cout << " sum of ends                       " << Pstart+Pend << G4endl;
         G4cout << " Original                          " << hadron->Get4Momentum() << G4endl;
 #endif
         return string;
+        return string;
+}
 …
 // --------- private methods ----------------------
+// ---------------------------------------------------------------------
 G4double G4DiffractiveExcitation::ChooseP(G4double Pmin, G4double Pmax) const // Uzhi
+{
 // choose an x between Xmin and Xmax with P(x) ~ 1/x
 //  to be improved...
         G4double range=Pmax-Pmin;                                         // Uzhi
         if ( Pmin <= 0. || range <=0. )
+        if ( Pmin <= 0. || range <=0. )
+        {
                 G4cout << " Pmin, range : " << Pmin << " , " << range << G4endl;
 …
+}
+G4ThreeVector G4DiffractiveExcitation::GaussianPt(G4double AveragePt2, G4double maxPtSquare) const // Uzhi
+// ---------------------------------------------------------------------
+G4ThreeVector G4DiffractiveExcitation::GaussianPt(G4double AveragePt2,
+                                                  G4double maxPtSquare) const // Uzhi
 {            //  @@ this method is used in FTFModel as well. Should go somewhere common!
         G4double Pt2;
 /*                                                                          // Uzhi
 …
 */                                                                          // Uzhi
+        Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() * (std::exp(-maxPtSquare/AveragePt2)-1.));// Uzhi
+        Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() *
+                           (std::exp(-maxPtSquare/AveragePt2)-1.));// Uzhi
         G4double Pt=std::sqrt(Pt2);
         G4double phi=G4UniformRand() * twopi;
+        return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);
+}
+        return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);
+}
+// ---------------------------------------------------------------------
 G4DiffractiveExcitation::G4DiffractiveExcitation(const G4DiffractiveExcitation &)
+{

trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4DiffractiveHHScatterer.cc

-                      r819
+                      r962
 #include "G4KineticTrack.hh"
 #include "G4DiffractiveSplitableHadron.hh"
+#include "G4FTFParameters.hh"                            // Uzhi 21.04.08
 G4DiffractiveHHScatterer::G4DiffractiveHHScatterer()
 …
 {}
+// -------------------------------------------------------------------
 G4KineticTrackVector * G4DiffractiveHHScatterer::
 Scatter(const G4KineticTrack & aTrack, const G4KineticTrack & bTrack)
 …
   G4DiffractiveSplitableHadron aHadron(& aTrack);
   G4DiffractiveSplitableHadron bHadron(& bTrack);
+  if ( ! theExcitation->ExciteParticipants(& aHadron, & bHadron))
+  theParameters = new G4FTFParameters(aHadron.GetDefinition(), // -------- Uzhi 21.04.08
+.,1., 100.);
+                                          //s);// ------------------------- Uzhi 21.04.08
+  if ( ! theExcitation->ExciteParticipants(& aHadron,
+                                           & bHadron,
+                                           theParameters))     // -------- Uzhi 21.04.08
+  {
         return NULL;

trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4DiffractiveSplitableHadron.cc

-                      r819
+                      r962
 //
 //
 // $Id: G4DiffractiveSplitableHadron.cc,v 1.6 2006/06/29 20:54:36 gunter Exp $
 // GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4DiffractiveSplitableHadron.cc,v 1.7 2008/03/31 15:34:01 vuzhinsk Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //

trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4FTFModel.cc

-                      r819
+                      r962
 //
 //
 // $Id: G4FTFModel.cc,v 1.7 2007/04/24 10:32:59 gunter Exp $
 // GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4FTFModel.cc,v 1.13 2008/12/09 10:40:52 vuzhinsk Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 …
 #include "G4FTFModel.hh"
+#include "G4FTFParameters.hh"                            // Uzhi 29.03.08
 #include "G4FTFParticipants.hh"
 #include "G4InteractionContent.hh"
 …
 #include "G4ParticleDefinition.hh"
 #include "G4ios.hh"
+#include <utility>                                        // Uzhi 29.03.08
 // Class G4FTFModel
+G4FTFModel::G4FTFModel():theExcitation(new G4DiffractiveExcitation()) // Uzhi
+G4FTFModel::G4FTFModel():theExcitation(new G4DiffractiveExcitation()),
+                         theElastic(new G4ElasticHNScattering()) // Uzhi 29.03.08
+{
         G4VPartonStringModel::SetThisPointer(this);
+}
+G4FTFModel::G4FTFModel(G4double a, G4double b, G4double c):theExcitation(new G4DiffractiveExcitation())
+        theParameters=0;                                         // Uzhi 9.12.08
+}
+/*
+G4FTFModel::G4FTFModel(G4double , G4double , G4double ):theExcitation(new // Uzhi 9.12.08 G4DiffractiveExcitation())
+{
         G4VPartonStringModel::SetThisPointer(this);
 …
         G4VPartonStringModel::SetThisPointer(this);
+}
+*/
 G4FTFModel::~G4FTFModel()
+{}
+{
+   if( theParameters != 0 ) delete theParameters;             // Uzhi 5.12.08
+// Because FTF model can be called for various particles
+// theParameters must be erased at the end of each call.
+// Thus the delete is olso in G4FTFModel::GetStrings() method
+   if( theExcitation != 0 ) delete theExcitation;             // Uzhi 5.12.08
+   if( theElastic    != 0 ) delete theElastic;                // Uzhi 5.12.08
+}
 …
+}
+// ------------------------------------------------------------
 void G4FTFModel::Init(const G4Nucleus & aNucleus, const G4DynamicParticle & aProjectile)
+{
-        theParticipants.Init(aNucleus.GetN(),aNucleus.GetZ()); // Uzhi N-mass number Z-charge
         theProjectile = aProjectile;
+}
+//G4cout<<"G4FTFModel::Init "<<aNucleus.GetN()<<" "<<aNucleus.GetZ()<<G4endl;
+        theParticipants.Init(aNucleus.GetN(),aNucleus.GetZ());
+// Uzhi N-mass number Z-charge ------------------------- Uzhi 29.03.08
+// --- cms energy
+        G4double s = sqr( theProjectile.GetMass() ) +
+                     sqr( G4Proton::Proton()->GetPDGMass() ) +
+*theProjectile.GetTotalEnergy()*G4Proton::Proton()->GetPDGMass();
+/*
+G4cout << " primary Total E (GeV): " << theProjectile.GetTotalEnergy()/GeV << G4endl;
+G4cout << " primary Mass    (GeV): " << theProjectile.GetMass() /GeV << G4endl;
+G4cout << "cms std::sqrt(s) (GeV) = " << std::sqrt(s) / GeV << G4endl;
+*/
+      if( theParameters != 0 ) delete theParameters;                    // Uzhi 9.12.08
+      theParameters = new G4FTFParameters(theProjectile.GetDefinition(),
+                                          aNucleus.GetN(),aNucleus.GetZ(),
+                                          s);// ------------------------- Uzhi 19.04.08
+//theParameters->SetProbabilityOfElasticScatt(0.); // To turn on/off (1/0) elastic scattering
+}
+// ------------------------------------------------------------
 G4ExcitedStringVector * G4FTFModel::GetStrings()
+{
+        theParticipants.BuildInteractions(theProjectile);
+//G4cout<<"theParticipants.GetList"<<G4endl;
+        theParticipants.GetList(theProjectile,theParameters);
+//G4cout<<"ExciteParticipants()"<<G4endl;
         if (! ExciteParticipants()) return NULL;;
+//G4cout<<"theStrings = BuildStrings()"<<G4endl;
         G4ExcitedStringVector * theStrings = BuildStrings();
+//G4cout<<"Return to theStrings "<<G4endl;
+        if( theParameters != 0 )                              // Uzhi 9.12.08
+        {                                                     // Uzhi 9.12.08
+          delete theParameters;                               // Uzhi 9.12.08
+          theParameters=0;                                    // Uzhi 9.12.08
+        }                                                     // Uzhi 9.12.08
         return theStrings;
+}
+// ------------------------------------------------------------
 struct DeleteVSplitableHadron { void operator()(G4VSplitableHadron * aH){delete aH;} };
+G4ExcitedStringVector * G4FTFModel::BuildStrings()
+{
+// Loop over all collisions; find all primaries, and all target ( targets may
+//  be duplicate in the List ( to unique G4VSplitableHadrons)
+        G4ExcitedStringVector * strings;
+        strings = new G4ExcitedStringVector();
+        std::vector<G4VSplitableHadron *> primaries;
+        std::vector<G4VSplitableHadron *> targets;
+        theParticipants.StartLoop();    // restart a loop
+        while ( theParticipants.Next() )
+        {
+            const G4InteractionContent & interaction=theParticipants.GetInteraction();
+                 //  do not allow for duplicates ...
+            if ( primaries.end() == std::find(primaries.begin(), primaries.end(), interaction.GetProjectile()) )
+                primaries.push_back(interaction.GetProjectile());
+            if ( targets.end() == std::find(targets.begin(), targets.end(),interaction.GetTarget()) )
+                targets.push_back(interaction.GetTarget());
+        }
+//      G4cout << "BuildStrings prim/targ " << primaries.entries() << " , " <<
+//                                           targets.entries() << G4endl;
+        unsigned int ahadron;
+        for ( ahadron=0; ahadron < primaries.size() ; ahadron++)
+        {
+            G4bool isProjectile=true;
+            strings->push_back(theExcitation->String(primaries[ahadron], isProjectile));
+        }
+        for ( ahadron=0; ahadron < targets.size() ; ahadron++)
+        {
+            G4bool isProjectile=false;
+            strings->push_back(theExcitation->String(targets[ahadron], isProjectile));
+        }
+        std::for_each(primaries.begin(), primaries.end(), DeleteVSplitableHadron());
+        primaries.clear();
+        std::for_each(targets.begin(), targets.end(), DeleteVSplitableHadron());
+        targets.clear();
+        return strings;
+}
+// ------------------------------------------------------------
 G4bool G4FTFModel::ExciteParticipants()
+{
+/*    // Uzhi 29.03.08                     For elastic Scatt.
+G4cout<<"  In ExciteParticipants() "<<theParticipants.theInteractions.size()<<G4endl;
+G4cout<<" test Params Tot "<<theParameters->GetTotalCrossSection()<<G4endl;
+G4cout<<" test Params Ela "<<theParameters->GetElasticCrossSection()<<G4endl;
+G4int counter=0;
+*/   // Uzhi 29.03.08
         while (theParticipants.Next())
+        {
            const G4InteractionContent & collision=theParticipants.GetInteraction();
+//G4cout << " soft colls : " << collision.GetNumberOfSoftCollisions() << G4endl; // Uzhi no match
+/*
+counter++;
+G4cout<<" Inter # "<<counter<<G4endl;
+*/
            G4VSplitableHadron * projectile=collision.GetProjectile();
            G4VSplitableHadron * target=collision.GetTarget();
+           if ( ! theExcitation->ExciteParticipants(projectile, target) )
+//   // Uzhi 29.03.08
+           G4bool Successfull;
+           if(G4UniformRand()< theParameters->GetProbabilityOfElasticScatt())
+           {
+//G4cout<<"Elastic"<<G4endl;
+            Successfull=theElastic->ElasticScattering(projectile, target, theParameters);
+           }
+           else
+           {
+//G4cout<<"Inelastic"<<G4endl;
+            Successfull=theExcitation->ExciteParticipants(projectile, target, theParameters);
+           }
+//           if(!Successfull)
+// // Uzhi 29.03.08
+//         if ( ! theExcitation->ExciteParticipants(projectile, target) )
+           if(!Successfull)
+           {
 //           give up, clean up
 …
                     const G4InteractionContent & interaction=theParticipants.GetInteraction();
                          //  do not allow for duplicates ...
+                    if ( primaries.end() == std::find(primaries.begin(), primaries.end(), interaction.GetProjectile()) )
+                    if ( primaries.end() == std::find(primaries.begin(), primaries.end(),
+                                                      interaction.GetProjectile()) )
                         primaries.push_back(interaction.GetProjectile());
+                    if ( targets.end() == std::find(targets.begin(), targets.end(),interaction.GetTarget()) )
+                    if ( targets.end()   == std::find(targets.begin(), targets.end(),
+                                                      interaction.GetTarget()) )
                         targets.push_back(interaction.GetTarget());
+                }
 …
                 return false;
+           }
+           }  // End of the loop Uzhi
+        }
+        return true;
+}
+// ------------------------------------------------------------
+G4ExcitedStringVector * G4FTFModel::BuildStrings()
+{
+// Loop over all collisions; find all primaries, and all target ( targets may
+//  be duplicate in the List ( to unique G4VSplitableHadrons)
+        G4ExcitedStringVector * strings;
+        strings = new G4ExcitedStringVector();
+        std::vector<G4VSplitableHadron *> primaries;
+        std::vector<G4VSplitableHadron *> targets;
+        theParticipants.StartLoop();    // restart a loop
+        while ( theParticipants.Next() )
+        {
+            const G4InteractionContent & interaction=theParticipants.GetInteraction();
+                 //  do not allow for duplicates ...
+            if ( primaries.end() == std::find(primaries.begin(), primaries.end(),
+                                                interaction.GetProjectile()) )
+                primaries.push_back(interaction.GetProjectile());
+            if ( targets.end()   == std::find(targets.begin(), targets.end(),
+                                                interaction.GetTarget()) )
+                targets.push_back(interaction.GetTarget());
+        }
+        return true;
+}
+//      G4cout << "BuildStrings prim/targ " << primaries.size() << " , " <<
+//                                             targets.size() << G4endl;
+        unsigned int ahadron;
+// Only for hA-interactions Uzhi -------------------------------------
+        for ( ahadron=0; ahadron < primaries.size() ; ahadron++)
+        {
+//G4ThreeVector aPosition=primaries[ahadron]->GetPosition();
+//G4cout<<"Proj Build "<<aPosition<<" "<<primaries[ahadron]->GetTimeOfCreation()<<G4endl;
+            G4bool isProjectile=true;
+            strings->push_back(theExcitation->String(primaries[ahadron], isProjectile));
+        }
+        for ( ahadron=0; ahadron < targets.size() ; ahadron++)
+        {
+//G4ThreeVector aPosition=targets[ahadron]->GetPosition();
+//G4cout<<"Targ Build "<<aPosition<<" "<<targets[ahadron]->GetTimeOfCreation()<<G4endl;
+            G4bool isProjectile=false;
+            strings->push_back(theExcitation->String(targets[ahadron], isProjectile));
+        }
+        std::for_each(primaries.begin(), primaries.end(), DeleteVSplitableHadron());
+        primaries.clear();
+        std::for_each(targets.begin(), targets.end(), DeleteVSplitableHadron());
+        targets.clear();
+        return strings;
+}
+// ------------------------------------------------------------

trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4FTFParticipants.cc

-                      r819
+                      r962
 //
 //
 // $Id: G4FTFParticipants.cc,v 1.7 2007/04/24 10:33:00 gunter Exp $
 // GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4FTFParticipants.cc,v 1.9 2008/06/13 12:49:23 vuzhinsk Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // ------------------------------------------------------------
 …
 // ------------------------------------------------------------
+#include "G4FTFParameters.hh"                            // Uzhi 29.03.08
 #include "G4FTFParticipants.hh"
 #include "G4DiffractiveSplitableHadron.hh"
 #include "G4VSplitableHadron.hh"
-//#include "G4PomeronCrossSection.hh"                      // Uzhi
-#include "G4FTFCrossSection.hh"                            // Uzhi
 #include "Randomize.hh"
+#include <utility>
+#include <utility>                                        // Uzhi 29.03.08
 // Class G4FTFParticipants
 G4FTFParticipants::G4FTFParticipants()
 …
 //{}
+void G4FTFParticipants::BuildInteractions(const G4ReactionProduct  &thePrimary)
+void G4FTFParticipants::GetList(const G4ReactionProduct  &thePrimary,
+                                      G4FTFParameters    *theParameters) // Uzhi 29.03.08
+{
 …
     theInteractions.clear();
+// --- cms energy
+    G4double s = sqr( thePrimary.GetMass() ) +
+                 sqr( G4Proton::Proton()->GetPDGMass() ) +
+*thePrimary.GetTotalEnergy()*G4Proton::Proton()->GetPDGMass();
+//    G4cout << " primary Total E (GeV): " << thePrimary.GetTotalEnergy()/GeV << G4endl;
+//    G4cout << " primary Mass    (GeV): " << thePrimary.GetMass() /GeV << G4endl;
+//    G4cout << "cms std::sqrt(s) (GeV) = " << std::sqrt(s) / GeV << G4endl;
+//    G4PomeronCrossSection theCrossSection(thePrimary.GetDefinition());      // Uzhi
+      G4FTFCrossSection theCrossSection(thePrimary.GetDefinition(),s);        // Uzhi
+    G4double deltaxy=2 * fermi;
+    G4double deltaxy=2 * fermi;                       // Extra nuclear radius
     G4VSplitableHadron * primarySplitable=new G4DiffractiveSplitableHadron(thePrimary);
     G4double xyradius;
     xyradius =theNucleus->GetOuterRadius() + deltaxy;
+    G4double xyradius;
+    xyradius =theNucleus->GetOuterRadius() + deltaxy; // Impact parameter sampling
+                                                      // radius
     G4bool nucleusNeedsShift = true;
 …
         theNucleus->StartLoop();
         G4Nucleon * nucleon;
+        while ( (nucleon=theNucleus->GetNextNucleon()) )
+//G4int InterNumber=0;           // Uzhi
+//while ( (nucleon=theNucleus->GetNextNucleon())&& (InterNumber < 1) ) // Uzhi
+        while ( (nucleon=theNucleus->GetNextNucleon()) ) // Uzhi
+        {
            G4double impact2= sqr(impactX - nucleon->GetPosition().x()) +
+                    sqr(impactY - nucleon->GetPosition().y());
+//         if ( theCrossSection.GetInelasticProbability(s,impact2)                       // Uzhi
+           if ( theCrossSection.GetInelasticProbability(  impact2/fermi/fermi)           // Uzhi
+                             sqr(impactY - nucleon->GetPosition().y());
+//         if ( theParameters->GetInelasticProbability(impact2/fermi/fermi) // Uzhi 29.03.08
+           if ( theParameters->GetProbabilityOfInteraction(impact2/fermi/fermi) // Uzhi 29.03.08
                 > G4UniformRand() )
+           {
+//InterNumber++;
                 if ( nucleusNeedsShift )
                 {                       // on the first hit, shift nucleus
 …
                     nucleon->Hit(targetSplitable);
+                }
+                G4InteractionContent * aInteraction = new G4InteractionContent(primarySplitable);
+                G4InteractionContent * aInteraction =
+                                       new G4InteractionContent(primarySplitable);
                 aInteraction->SetTarget(targetSplitable);
                 theInteractions.push_back(aInteraction);
 …
 // Implementation (private) methods

trunk/source/processes/hadronic/models/parton_string/hadronization/include/G4ExcitedStringDecay.hh

r819	r962
26	26	//
27	27	// $Id: G4ExcitedStringDecay.hh,v 1.7 2007/05/03 22:06:17 gunter Exp $
28		// GEANT4 tag $Name: $
	28	// GEANT4 tag $Name: geant4-09-02-ref-02 $
29	29	//
30	30	#ifndef G4ExcitedStringDecay_h

trunk/source/processes/hadronic/models/parton_string/hadronization/include/G4FragmentingString.hh

-                      r819
+                      r962
 //
 //
 // $Id: G4FragmentingString.hh,v 1.3 2006/06/29 20:54:44 gunter Exp $
 // GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4FragmentingString.hh,v 1.4 2007/12/20 15:38:06 vuzhinsk Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 …
                           G4ParticleDefinition * newdecay,
                           const G4LorentzVector *momentum);
+      G4FragmentingString(const G4FragmentingString &old,      // Uzhi
+                          G4ParticleDefinition * newdecay);    // Uzhi
       ~G4FragmentingString();
 …
       G4double Mass() const;
       G4double Mass2() const;
+      G4double MassT2() const;
       G4ParticleDefinition* GetLeftParton(void) const;

trunk/source/processes/hadronic/models/parton_string/hadronization/include/G4HadronBuilder.hh

r819	r962
26	26	//
27	27	// $Id: G4HadronBuilder.hh,v 1.3 2006/06/29 20:54:46 gunter Exp $
28		// GEANT4 tag $Name: $
	28	// GEANT4 tag $Name: geant4-09-02-ref-02 $
29	29	//
30	30	// -----------------------------------------------------------------------------

trunk/source/processes/hadronic/models/parton_string/hadronization/include/G4LundStringFragmentation.hh

-                      r819
+                      r962
 //
 //
 // $Id: G4LundStringFragmentation.hh,v 1.4 2007/04/24 14:55:23 gunter Exp $
 // GEANT4 tag $Name: geant4-09-01-patch-02 $ Maxim Komogorov
+// $Id: G4LundStringFragmentation.hh,v 1.6 2008/04/25 14:20:14 vuzhinsk Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $ Maxim Komogorov
 //
 // -----------------------------------------------------------------------------
 …
+    {
 public:
     G4LundStringFragmentation();
+//    G4LundStringFragmentation(G4double sigmaPt);
+     G4LundStringFragmentation(const G4LundStringFragmentation &right);
+     virtual ~G4LundStringFragmentation();
+     virtual G4KineticTrackVector* FragmentString(const G4ExcitedString& theString);
+    G4LundStringFragmentation(const G4LundStringFragmentation &right);
+    virtual ~G4LundStringFragmentation();
-public:
     const G4LundStringFragmentation & operator=(const G4LundStringFragmentation &right);
     int operator==(const G4LundStringFragmentation &right) const;
     int operator!=(const G4LundStringFragmentation &right) const;
+    virtual G4KineticTrackVector* FragmentString(const G4ExcitedString& theString);
 private:
+   virtual G4double GetLightConeZ(G4double zmin, G4double zmax, G4int PartonEncoding,  G4ParticleDefinition* pHadron, G4double Px, G4double Py);
+   void SetMinimalStringMass(const G4FragmentingString  * const string);
+   void SetMinimalStringMass2(const G4double aValue);
-   virtual void Sample4Momentum(G4LorentzVector* Mom, G4double Mass, G4LorentzVector* AntiMom, G4double AntiMass, G4double InitialMass);
    virtual G4bool StopFragmenting(const G4FragmentingString  * const string);
    virtual G4bool IsFragmentable(const G4FragmentingString * const string);
+   virtual G4LorentzVector * SplitEandP(G4ParticleDefinition * pHadron, G4FragmentingString * string);
    virtual G4bool SplitLast(G4FragmentingString * string,
+                    G4KineticTrackVector * LeftVector,
+                    G4KineticTrackVector * RightVector);
+   void SetMinimalStringMass(const G4FragmentingString  * const string);
+   void SetMinimalStringMass2(const G4double aValue);
+                            G4KineticTrackVector * LeftVector,
+                            G4KineticTrackVector * RightVector);
+   virtual void Sample4Momentum(G4LorentzVector* Mom,     G4double Mass,
+                                G4LorentzVector* AntiMom, G4double AntiMass,
+                                G4double InitialMass);
+   virtual G4LorentzVector * SplitEandP(G4ParticleDefinition * pHadron,
+                                        G4FragmentingString * string,
+                                        G4FragmentingString * newString); // Uzhi
+   virtual G4double GetLightConeZ(G4double zmin, G4double zmax,
+                                  G4int PartonEncoding,
+                                  G4ParticleDefinition* pHadron,
+                                  G4double Px, G4double Py);
 private:
+// ------ For estimation of a minimal string mass ---------------
+   G4double Mass_of_light_quark;
+   G4double Mass_of_heavy_quark;
+   G4double Mass_of_string_junction;
+// ------ An estimated minimal string mass ----------------------
    G4double MinimalStringMass;
    G4double MinimalStringMass2;
+// ------ Minimal invariant mass used at a string fragmentation -
    G4double WminLUND;
 };

trunk/source/processes/hadronic/models/parton_string/hadronization/include/G4QGSMFragmentation.hh

-                      r819
+                      r962
 //
 //
 // $Id: G4QGSMFragmentation.hh,v 1.4 2007/04/24 14:55:23 gunter Exp $
 // GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4QGSMFragmentation.hh,v 1.5 2007/12/20 15:38:07 vuzhinsk Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // -----------------------------------------------------------------------------
 …
    virtual G4bool StopFragmenting(const G4FragmentingString  * const string);
    virtual G4bool IsFragmentable(const G4FragmentingString * const string);
+   virtual G4LorentzVector * SplitEandP(G4ParticleDefinition * pHadron, G4FragmentingString * string);
+   virtual G4LorentzVector * SplitEandP(G4ParticleDefinition * pHadron,
+                                        G4FragmentingString * string,     // Uzhi
+                                        G4FragmentingString * newString); // Uzhi
    virtual G4bool SplitLast(G4FragmentingString * string,
                     G4KineticTrackVector * LeftVector,

trunk/source/processes/hadronic/models/parton_string/hadronization/include/G4VKinkyStringDecay.hh

r819	r962
26	26	//
27	27	// $Id: G4VKinkyStringDecay.hh,v 1.3 2006/06/29 20:54:53 gunter Exp $
28		// GEANT4 tag $Name: $
	28	// GEANT4 tag $Name: geant4-09-02-ref-02 $
29	29	// Maxim Komogorov
30	30	//

trunk/source/processes/hadronic/models/parton_string/hadronization/include/G4VLongitudinalStringDecay.hh

-                      r819
+                      r962
 //
 //
 // $Id: G4VLongitudinalStringDecay.hh,v 1.4 2007/04/24 14:55:23 gunter Exp $
 // GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4VLongitudinalStringDecay.hh,v 1.6 2008/06/23 08:35:54 vuzhinsk Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 // Maxim Komogorov
 //
 …
 class G4FragmentingString;
 //**********************************************************************************************
+//**************************************************************************************
 class G4VLongitudinalStringDecay
+   {
 public:
    G4VLongitudinalStringDecay();
+            G4VLongitudinalStringDecay();
    virtual ~G4VLongitudinalStringDecay();
 private:
+//  G4VLongitudinalStringDecay(const G4VLongitudinalStringDecay &right);
+//  const G4VLongitudinalStringDecay & operator=(const G4VLongitudinalStringDecay &right);
    int operator==(const G4VLongitudinalStringDecay &right) const;
    int operator!=(const G4VLongitudinalStringDecay &right) const;
 …
 public:
    virtual G4KineticTrackVector* FragmentString(const G4ExcitedString& theString)=0;
+protected:
+// For changing Mass Cut used for selection of very small mass strings
+   virtual void SetMassCut(G4double aValue);
+// For handling a string with very low mass
+   G4KineticTrackVector * LightFragmentationTest(const G4ExcitedString * const theString);
+// To store created quarks or 2 last hadrons
+   typedef std::pair<G4ParticleDefinition*, G4ParticleDefinition*> pDefPair;
+// For creation of hadrons from given quark pair
+   typedef G4ParticleDefinition * (G4HadronBuilder::*Pcreate)
+                                        (G4ParticleDefinition*, G4ParticleDefinition*);
+//-----------------------------------------------------------------------------
+// Used by LightFragmentationTest for estimation of lowest possible mass of
+// given quark system
+   G4double FragmentationMass(const G4FragmentingString * const string,
+                              Pcreate build=0,
+                              pDefPair * pdefs=0);
+   G4ParticleDefinition* FindParticle(G4int Encoding);
+   virtual void Sample4Momentum(G4LorentzVector* Mom,     G4double Mass,
+                                G4LorentzVector* AntiMom, G4double AntiMass,
+                                G4double InitialMass)=0;
+//-----------------------------------------------------------------------------
+// For decision on continue or stop string fragmentation
+   virtual G4bool StopFragmenting(const G4FragmentingString  * const string)=0;
+   virtual G4bool IsFragmentable(const G4FragmentingString * const string)=0;
+// If a string can not fragment, make last break into 2 hadrons
+   virtual G4bool SplitLast(G4FragmentingString * string,
+                    G4KineticTrackVector * LeftVector,
+                    G4KineticTrackVector * RightVector)=0;
+//-----------------------------------------------------------------------------
+// If a string fragments, do the following
+// For transver of a string to its CMS frame
+   G4ExcitedString *CPExcited(const G4ExcitedString& string);
+   G4KineticTrack * Splitup(G4FragmentingString *string,
+                            G4FragmentingString *&newString);
+   G4ParticleDefinition * QuarkSplitup(G4ParticleDefinition* decay,
+                                       G4ParticleDefinition *&created);
+   G4ParticleDefinition * DiQuarkSplitup(G4ParticleDefinition* decay,
+                                         G4ParticleDefinition *&created);
+   pDefPair CreatePartonPair(G4int NeedParticle, G4bool AllowDiquarks=true);
+public:
+//   used by G4VKinkyStringDecy..
+   G4int SampleQuarkFlavor(void);
+   G4ThreeVector SampleQuarkPt();
+protected:
+//-----------------------------------------------------------------------------
+// For determination of kinematical properties of created hadron
+//   virtual G4LorentzVector * SplitEandP(G4ParticleDefinition * pHadron,    // Uzhi
+//                                        G4FragmentingString * string  )=0; // Uzhi
+   virtual G4LorentzVector * SplitEandP(G4ParticleDefinition * pHadron,      // Uzhi
+                                        G4FragmentingString * string,        // Uzhi
+                                        G4FragmentingString * newString  )=0;// Uzhi
+   virtual G4double GetLightConeZ(G4double zmin, G4double zmax,
+                                  G4int PartonEncoding,
+                                  G4ParticleDefinition* pHadron,
+                                  G4double Px, G4double Py       ) = 0;
+   void CalculateHadronTimePosition(G4double theInitialStringMass,
+                                    G4KineticTrackVector *);
+// Used for some test purposes ------------------------------------------------
+   void ConstructParticle();
+   G4ParticleDefinition* CreateHadron(G4int id1, G4int id2,
+                                      G4bool theGivenSpin, G4int theSpin);
+//-----------------------------------------------------------------------------
+public:
    G4KineticTrackVector* DecayResonans (G4KineticTrackVector* aHadrons);
    void SetSigmaTransverseMomentum(G4double aQT);
    void SetStrangenessSuppression(G4double aValue);
 …
    void SetVectorMesonMixings( std::vector<G4double> aVector);
+//   used by G4VKinkyStringDecy..
+   G4int SampleQuarkFlavor(void);
+   G4ThreeVector SampleQuarkPt();
+   void SetStringTensionParameter(G4double aValue);            // Uzhi 20 June 08
 //private:
 protected:
 …
    G4double GetClusterMass()            {return ClusterMass;};
    G4int    GetClusterLoopInterrupt()   {return ClusterLoopInterrupt;};
+   G4ParticleDefinition* CreateHadron(G4int id1, G4int id2, G4bool theGivenSpin, G4int theSpin);
+   virtual void Sample4Momentum(G4LorentzVector* Mom, G4double Mass, G4LorentzVector* AntiMom, G4double AntiMass, G4double InitialMass)=0;
+protected:
+   // Additional protected declarations
+   virtual G4double GetLightConeZ(G4double zmin, G4double zmax, G4int PartonEncoding,  G4ParticleDefinition* pHadron, G4double Px, G4double Py) = 0;
+   G4double GetStringTensionParameter() {return Kappa;};       // Uzhi 20 June 08
 //private:
 protected:
 …
    G4HadronBuilder *hadronizer;
-   void ConstructParticle();
    G4double pspin_meson;
    G4double pspin_barion;
 …
    G4bool    PastInitPhase;
+   G4KineticTrackVector * LightFragmentationTest(const G4ExcitedString * const theString);
+   virtual G4bool StopFragmenting(const G4FragmentingString  * const string)=0;
+   virtual G4bool IsFragmentable(const G4FragmentingString * const string)=0;
+   G4double Kappa; // String tension parameter                 // Uzhi 20 June 08
 //   G4double MinFragmentationMass(G4ExcitedString * theString,
 //                              G4ParticleDefinition*& Hadron1,
 //                              G4ParticleDefinition*& Hadron2);
-   typedef std::pair<G4ParticleDefinition*, G4ParticleDefinition*> pDefPair;
-   typedef G4ParticleDefinition * (G4HadronBuilder::*Pcreate)
-                                        (G4ParticleDefinition*, G4ParticleDefinition*);
-   G4double FragmentationMass(
-                const G4FragmentingString * const string,
-                Pcreate build=0,
-                pDefPair * pdefs=0);
-   G4KineticTrack * Splitup(G4FragmentingString *string, G4FragmentingString *&newString);
-   virtual G4LorentzVector * SplitEandP(G4ParticleDefinition * pHadron, G4FragmentingString * string)=0;
-   virtual G4bool SplitLast(G4FragmentingString * string,
-                    G4KineticTrackVector * LeftVector,
-                    G4KineticTrackVector * RightVector)=0;
-   void CalculateHadronTimePosition(G4double theInitialStringMass, G4KineticTrackVector *);
-   G4ExcitedString *CPExcited(const G4ExcitedString& string);
-   G4ParticleDefinition* FindParticle(G4int Encoding);
-   // Additional Implementation Declarations
-   G4ParticleDefinition * QuarkSplitup(G4ParticleDefinition* decay,
-                                G4ParticleDefinition *&created);
-   G4ParticleDefinition * DiQuarkSplitup(G4ParticleDefinition* decay,
-                                        G4ParticleDefinition *&created);
-   pDefPair CreatePartonPair(G4int NeedParticle, G4bool AllowDiquarks=true);
 };
+//**********************************************************************************************
+//*************************************************************************************
 // Class G4VLongitudinalStringDecay
 #endif

trunk/source/processes/hadronic/models/parton_string/hadronization/src/G4FragmentingString.cc

-                      r819
+                      r962
         if ( old.decaying == Left )
+        {
+//G4cout<<" Left "<<G4endl;
+//G4cout<<"Pt right "<<Ptright<<G4endl;
+//G4cout<<"Pt left  "<<Ptleft <<G4endl;
                 RightParton= old.RightParton;
                 Ptright    = old.Ptright;
 …
                 Ptleft     = old.Ptleft - momentum->vect();
                 Ptleft.setZ(0.);
+//G4cout<<"Pt right "<<Ptright<<G4endl;
+//G4cout<<"Pt left  "<<Ptleft <<G4endl;
         } else if ( old.decaying == Right )
+        {
+//G4cout<<" Right "<<G4endl;
+//G4cout<<"Pt right "<<Ptright<<G4endl;
+//G4cout<<"Pt left  "<<Ptleft <<G4endl;
                 RightParton = newdecay;
                 Ptright     = old.Ptright - momentum->vect();
 …
                 LeftParton  = old.LeftParton;
                 Ptleft      = old.Ptleft;
+//G4cout<<"Pt right "<<Ptright<<G4endl;
+//G4cout<<"Pt left  "<<Ptleft <<G4endl;
         } else
+        {
 …
 //---------------------------------------------------------------------------------
+G4FragmentingString::G4FragmentingString(const G4FragmentingString &old,  // Uzhi
+                                         G4ParticleDefinition * newdecay) // Uzhi
+{                                                                         // Uzhi
+        decaying=None;                                                    // Uzhi
+        if ( old.decaying == Left )                                       // Uzhi
+        {                                                                 // Uzhi
+                RightParton= old.RightParton;                             // Uzhi
+                LeftParton = newdecay;                                    // Uzhi
+        } else if ( old.decaying == Right )                               // Uzhi
+        {                                                                 // Uzhi
+                RightParton = newdecay;                                   // Uzhi
+                LeftParton  = old.LeftParton;                             // Uzhi
+        } else                                                            // Uzhi
+        {
+                throw G4HadronicException(__FILE__, __LINE__, "G4FragmentingString::G4FragmentingString: no decay Direction defined");
+        }
+}
+//---------------------------------------------------------------------------------
 G4FragmentingString::~G4FragmentingString()
 {}
 …
         return std::sqrt(this->Mass2());
+}
+G4double G4FragmentingString::MassT2() const
+{
+        return Pplus*Pminus;
+}

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

-                      r819
+                      r962
 //
 //
 // $Id: G4HadronBuilder.cc,v 1.6 2006/06/29 20:55:01 gunter Exp $
 // GEANT4 tag $Name:  $
+// $Id: G4HadronBuilder.cc,v 1.7 2008/04/25 14:20:14 vuzhinsk Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // -----------------------------------------------------------------------------

trunk/source/processes/hadronic/models/parton_string/hadronization/src/G4LundStringFragmentation.cc

-                      r819
+                      r962
 //
 //
 // $Id: G4LundStringFragmentation.cc,v 1.7 2007/04/24 14:55:23 gunter Exp $
 // GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4LundStringFragmentation.cc,v 1.13 2008/06/23 09:17:10 gcosmo Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $ 1.8
 //
 // -----------------------------------------------------------------------------
 …
 // Class G4LundStringFragmentation
 //****************************************************************************************
+//*************************************************************************************
 G4LundStringFragmentation::G4LundStringFragmentation()
+   {
+   MinimalStringMass  = 0.;             // Uzhi
+   MinimalStringMass2 = 0.;             // Uzhi
+   WminLUND = 1.*GeV;                   // Uzhi
+   SmoothParam  = 0.2;                  // Uzhi
+// ------ For estimation of a minimal string mass ---------------
+    Mass_of_light_quark    =140.*MeV;
+    Mass_of_heavy_quark    =500.*MeV;
+    Mass_of_string_junction=720.*MeV;
+// ------ An estimated minimal string mass ----------------------
+    MinimalStringMass  = 0.;
+    MinimalStringMass2 = 0.;
+// ------ Minimal invariant mass used at a string fragmentation -
+    WminLUND = 0.7*GeV;                   // Uzhi 0.8 1.5
+// ------ Smooth parameter used at a string fragmentation for ---
+// ------ smearinr sharp mass cut-off ---------------------------
+    SmoothParam  = 0.2;
+    SetStringTensionParameter(0.25);                           // Uzhi 20 June 08
+   }
+// G4LundStringFragmentation::G4LundStringFragmentation(G4double sigmaPt)
+// : G4VLongitudinalStringDecay(sigmaPt)
+//    {
+//    }
+// --------------------------------------------------------------
 G4LundStringFragmentation::G4LundStringFragmentation(const G4LundStringFragmentation &) : G4VLongitudinalStringDecay()
+   {
+   }
 G4LundStringFragmentation::~G4LundStringFragmentation()
+   {
+   }
 //****************************************************************************************
+//*************************************************************************************
 const G4LundStringFragmentation & G4LundStringFragmentation::operator=(const G4LundStringFragmentation &)
 …
+   }
+//****************************************************************************************
+//----------------------------------------------------------------------------------------------------------
+G4KineticTrackVector* G4LundStringFragmentation::FragmentString(const G4ExcitedString& theString)
+//--------------------------------------------------------------------------------------
+void G4LundStringFragmentation::SetMinimalStringMass(const G4FragmentingString  * const string)  // Uzhi
+{
+/*
+G4cout<<"In SetMinMass -------------------"<<std::sqrt(string->Mass2())<<G4endl;
+G4cout<<string->GetLeftParton()->GetPDGEncoding()<<" "<<
+        string->GetRightParton()->GetPDGEncoding()<<G4endl;
+*/
+  G4double EstimatedMass=0.;
+  G4int Number_of_quarks=0;
+  G4int Qleft =std::abs(string->GetLeftParton()->GetPDGEncoding());
+  if( Qleft > 1000)
+    {
+     Number_of_quarks+=2;
+     G4int q1=Qleft/1000;
+     if( q1 < 3) {EstimatedMass +=Mass_of_light_quark;}
+     if( q1 > 2) {EstimatedMass +=Mass_of_heavy_quark;}
+     G4int q2=(Qleft/100)%10;
+     if( q2 < 3) {EstimatedMass +=Mass_of_light_quark;}
+     if( q2 > 2) {EstimatedMass +=Mass_of_heavy_quark;}
+     EstimatedMass +=Mass_of_string_junction;
+    }
+  else
+    {
+     Number_of_quarks++;
+     if( Qleft < 3) {EstimatedMass +=Mass_of_light_quark;}
+     if( Qleft > 2) {EstimatedMass +=Mass_of_heavy_quark;}
+    }
+  G4int Qright=std::abs(string->GetRightParton()->GetPDGEncoding());
+  if( Qright > 1000)
+    {
+     Number_of_quarks+=2;
+     G4int q1=Qright/1000;
+     if( q1 < 3) {EstimatedMass +=Mass_of_light_quark;}
+     if( q1 > 2) {EstimatedMass +=Mass_of_heavy_quark;}
+     G4int q2=(Qright/100)%10;
+     if( q2 < 3) {EstimatedMass +=Mass_of_light_quark;}
+     if( q2 > 2) {EstimatedMass +=Mass_of_heavy_quark;}
+     EstimatedMass +=Mass_of_string_junction;
+    }
+  else
+    {
+     Number_of_quarks++;
+     if( Qright < 3) {EstimatedMass +=Mass_of_light_quark;}
+     if( Qright > 2) {EstimatedMass +=Mass_of_heavy_quark;}
+    }
+  if(Number_of_quarks==2){EstimatedMass +=100.*MeV;}
+  if(Number_of_quarks==3){EstimatedMass += 20.*MeV;}
+  if(Number_of_quarks==4){EstimatedMass -=2.*Mass_of_string_junction;
+                          if(EstimatedMass <= 1600.*MeV){EstimatedMass-=200.*MeV;}
+                          else                          {EstimatedMass+=100.*MeV;}
+                         }
+  MinimalStringMass=EstimatedMass;
+  SetMinimalStringMass2(EstimatedMass);
+//G4cout<<"Out SetMinimalStringMass "<<MinimalStringMass<<G4endl;
+}
+//--------------------------------------------------------------------------------------
+void G4LundStringFragmentation::SetMinimalStringMass2(
+                                 const G4double aValue)
+{
+  MinimalStringMass2=aValue * aValue;
+}
+//--------------------------------------------------------------------------------------
+G4KineticTrackVector* G4LundStringFragmentation::FragmentString(
+                                const G4ExcitedString& theString)
+{
 …
 //      check if string has enough mass to fragment...
+        SetMassCut(160.*MeV); // For LightFragmentationTest it is required
+                              // that no one pi-meson can be produced
+/*
+G4cout<<G4endl<<"G4LundStringFragmentation::"<<G4endl;
+G4cout<<"FragmentString Position"<<theString.GetPosition()/fermi<<" "<<
+theString.GetTimeOfCreation()/fermi<<G4endl;
+G4cout<<"FragmentString Momentum"<<theString.Get4Momentum()<<theString.Get4Momentum().mag()<<G4endl;
+*/
         G4KineticTrackVector * LeftVector=LightFragmentationTest(&theString);
         if ( LeftVector != 0 ) {
+//G4cout<<"Return single hadron from string"<<G4endl;
+                                return LeftVector;}
+        LeftVector = new G4KineticTrackVector;
+//G4cout<<"Return single hadron insted of string"<<G4endl;
+// Uzhi insert 6.05.08 start
+          if(LeftVector->size() == 1){
+ // One hadron is saved in the interaction
+            LeftVector->operator[](0)->SetFormationTime(theString.GetTimeOfCreation());
+            LeftVector->operator[](0)->SetPosition(theString.GetPosition());
+/* // To set large formation time open *
+            LeftVector->operator[](0)->SetFormationTime(theString.GetTimeOfCreation()+100.*fermi);
+            LeftVector->operator[](0)->SetPosition(theString.GetPosition());
+            G4ThreeVector aPosition(theString.GetPosition().x(),
+                                    theString.GetPosition().y(),
+                                    theString.GetPosition().z()+100.*fermi);
+            LeftVector->operator[](0)->SetPosition(aPosition);
+*/
+//G4cout<<"Single hadron "<<LeftVector->operator[](0)->GetPosition()<<" "<<LeftVector->operator[](0)->GetFormationTime()<<G4endl;
+          } else {    // 2 hadrons created from qq-qqbar are stored
+            LeftVector->operator[](0)->SetFormationTime(theString.GetTimeOfCreation());
+            LeftVector->operator[](0)->SetPosition(theString.GetPosition());
+            LeftVector->operator[](1)->SetFormationTime(theString.GetTimeOfCreation());
+            LeftVector->operator[](1)->SetPosition(theString.GetPosition());
+          }
+// Uzhi insert 6.05.08 end
+          return LeftVector;
+        }
+//--------------------- The string can fragment -------------------------------
+//--------------- At least two particles can be produced ----------------------
+                               LeftVector =new G4KineticTrackVector;
         G4KineticTrackVector * RightVector=new G4KineticTrackVector;
-// this should work but its only a semi deep copy. %GF  G4ExcitedString theStringInCMS(theString);
         G4ExcitedString *theStringInCMS=CPExcited(theString);
         G4LorentzRotation toCms=theStringInCMS->TransformToAlignedCms();
 …
         G4int attempt=0;
         while ( !success && attempt++ < StringLoopInterrupt )
+        {
+        { // If the string fragmentation do not be happend, repeat the fragmentation---
                 G4FragmentingString *currentString=new G4FragmentingString(*theStringInCMS);
 //G4cout<<"Main FragmentString cur M2  "<<std::sqrt(currentString->Mass2())<<G4endl;
                 std::for_each(LeftVector->begin(), LeftVector->end(), DeleteKineticTrack());
+//G4cout<<"FragmentString cur M2  "<<std::sqrt(currentString->Mass2())<<G4endl;
+          // Cleaning up the previously produced hadrons ------------------------------
+                std::for_each(LeftVector->begin() , LeftVector->end() , DeleteKineticTrack());
                 LeftVector->clear();
                 std::for_each(RightVector->begin(), RightVector->end(), DeleteKineticTrack());
                 RightVector->clear();
+          // Main fragmentation loop until the string will not be able to fragment ----
                 inner_sucess=true;  // set false on failure..
                 while (! StopFragmenting(currentString) )
                 {  // Split current string into hadron + new string
-//                      G4FragmentingString *PreviousString=currentString;   // Uzhi
                         G4FragmentingString *newString=0;  // used as output from SplitUp...
 //G4cout<<"FragmentString to Splitup ===================================="<<G4endl;
+//G4cout<<"++++++++++++++++++++++++++ Enter num--------------------------"<<G4endl;
 //G4int Uzhi; G4cin>>Uzhi;                         // Uzhi
                         G4KineticTrack * Hadron=Splitup(currentString,newString);
 //G4cout<<" Hadron "<<Hadron<<G4endl;
+//                      if ( Hadron != 0 && IsFragmentable(newString))       // Uzhi
+                        if ( Hadron != 0 )                                   // Uzhi
+                        if ( Hadron != 0 )  // Store the hadron
+                        {
                            if ( currentString->GetDecayDirection() > 0 )
 …
                            delete currentString;
                            currentString=newString;
+                        } /* else {                                // Uzhi
+                         // abandon ... start from the beginning
+                           if (newString) delete newString;                  // ??? Uzhi local?
+                           if (Hadron)    delete Hadron;
+//                           currentString = PreviousString;                 // Uzhi
+                           inner_sucess=false;
+                           break;
+                        } */                                      // Uzhi
+//                        delete PreviousString;                             // ??? Uzhi local?
+                        }
                 };
                 // Split current string into 2 final Hadrons
+        // Split remaining string into 2 final Hadrons ------------------------
 //G4cout<<"FragmentString to SplitLast if inner_sucess#0"<<inner_sucess<<G4endl;
                 if ( inner_sucess &&                                         // Uzhi
+                if ( inner_sucess &&
                      SplitLast(currentString,LeftVector, RightVector) )
+                {
 …
+                }
                 delete currentString;
+        }
+        }  // End of the loop in attemps to fragment the string
         delete theStringInCMS;
 …
         G4LorentzRotation toObserverFrame(toCms.inverse());
+//        LeftVector->operator[](0)->SetFormationTime(theString.GetTimeOfCreation());
+//        LeftVector->operator[](0)->SetPosition(theString.GetPosition());
+        G4double      TimeOftheStringCreation=theString.GetTimeOfCreation();
+        G4ThreeVector PositionOftheStringCreation(theString.GetPosition());
+/*  // For large formation time open *
+        G4double      TimeOftheStringCreation=theString.GetTimeOfCreation()+100*fermi;
+        G4ThreeVector PositionOftheStringCreation(theString.GetPosition().x(),
+                                                  theString.GetPosition().y(),
+                                                  theString.GetPosition().z()+100*fermi);
+*/
+/*
+        if(theString.GetPosition().y() > 100.*fermi){
+// It is a projectile-like string -------------------------------------
+          G4double Zmin=theString.GetPosition().y()-1000.*fermi;
+          G4double Zmax=theString.GetPosition().z();
+          TimeOftheStringCreation=
+          (Zmax-Zmin)*theString.Get4Momentum().e()/theString.Get4Momentum().z();
+          G4ThreeVector aPosition(0.,0.,Zmax);
+          PositionOftheStringCreation=aPosition;
+        }
+*/
         for(size_t C1 = 0; C1 < LeftVector->size(); C1++)
+        {
 …
            Momentum = toObserverFrame*Momentum;
            Hadron->Set4Momentum(Momentum);
            G4LorentzVector Coordinate(Hadron->GetPosition(), Hadron->GetFormationTime());
            Momentum = toObserverFrame*Coordinate;
            Hadron->SetFormationTime(Momentum.e());
+           Hadron->SetFormationTime(TimeOftheStringCreation+Momentum.e());
            G4ThreeVector aPosition(Momentum.vect());
+           Hadron->SetPosition(theString.GetPosition()+aPosition);
+        }
+//         Hadron->SetPosition(theString.GetPosition()+aPosition);
+           Hadron->SetPosition(PositionOftheStringCreation+aPosition);
+//G4cout<<"Hadron "<<C1<<" "<<Hadron->GetPosition()/fermi<<" "<<Hadron->GetFormationTime()/fermi<<G4endl;
+        };
 //G4cout<<"Out FragmentString"<<G4endl;
         return LeftVector;
+}
+//----------------------------------------------------------------------------------
+G4bool G4LundStringFragmentation::IsFragmentable(const G4FragmentingString * const string)
+{
+//G4cout<<"In IsFragmentable"<<G4endl;
+  SetMinimalStringMass(string);                                                     // Uzhi
+//G4cout<<"Out IsFragmentable MinMass"<<MinimalStringMass<<" String Mass"<<std::sqrt(string->Get4Momentum().mag2())<<G4endl;
+  return sqr(MinimalStringMass + WminLUND) < string->Get4Momentum().mag2();         // Uzhi
+//      return sqr(FragmentationMass(string)+MassCut) <                             // Uzhi
+//                      string->Mass2();                                            // Uzhi
+}
+//----------------------------------------------------------------------------------------
+G4bool G4LundStringFragmentation::StopFragmenting(const G4FragmentingString * const string)
+{
+//G4cout<<"StopFragmenting"<<G4endl;
+  SetMinimalStringMass(string);
+//G4cout<<"StopFragm MinMass "<<MinimalStringMass<<" String Mass "<<std::sqrt(string->Get4Momentum().mag2())<<G4endl;
+  return (MinimalStringMass + WminLUND)*
+             (1 + SmoothParam * (1.-2*G4UniformRand())) >
+                   string->Mass();
+}
+//----------------------------------------------------------------------------------------
+G4bool G4LundStringFragmentation::SplitLast(G4FragmentingString * string,
+                                             G4KineticTrackVector * LeftVector,
+                                             G4KineticTrackVector * RightVector)
+{
+    //... perform last cluster decay
+/*
+G4cout<<"SplitLast String mass "<<string->Mass()<<G4endl;
+G4cout<<string->GetLeftParton()->GetPDGEncoding()<<" "<<G4endl;
+G4cout<<string->GetRightParton()->GetPDGEncoding()<<" "<<G4endl;
+*/
+    G4LorentzVector Str4Mom=string->Get4Momentum();
+//G4cout<<"String 4 momentum "<<Str4Mom<<G4endl;
+    G4ThreeVector ClusterVel =string->Get4Momentum().boostVector();
+    G4double ResidualMass   = string->Mass();
+    G4ParticleDefinition * LeftHadron, * RightHadron;
+    G4int cClusterInterrupt = 0;
+    do
+    {
+//G4cout<<" Cicle "<<cClusterInterrupt<<" "<< ClusterLoopInterrupt<<G4endl;
+        if (cClusterInterrupt++ >= ClusterLoopInterrupt)
+        {
+          return false;
+        }
+        G4ParticleDefinition * quark = NULL;
+        string->SetLeftPartonStable(); // to query quark contents..
+        if (!string->FourQuarkString() )
+        {
+            // The string is q-qbar, or q-qq, or qbar-qqbar type
+           if (string->DecayIsQuark() && string->StableIsQuark() )
+           {
+            //... there are quarks on cluster ends
+              LeftHadron= QuarkSplitup(string->GetLeftParton(), quark);
+           } else
+           {
+           //... there is a Diquark on one of the cluster ends
+              G4int IsParticle;
+              if ( string->StableIsQuark() )
+              {
+                 IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? -1 : +1;
+              } else
+              {
+                 IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? +1 : -1;
+              }
+              pDefPair QuarkPair = CreatePartonPair(IsParticle,false);  // no diquarks wanted
+              quark = QuarkPair.second;
+              LeftHadron=hadronizer->Build(QuarkPair.first, string->GetLeftParton());
+           }
+           RightHadron = hadronizer->Build(string->GetRightParton(), quark);
+        } else
+        {
+            // The string is qq-qqbar type. Diquarks are on the string ends
+            G4int LiftQuark1= string->GetLeftParton()->GetPDGEncoding()/1000;
+            G4int LiftQuark2=(string->GetLeftParton()->GetPDGEncoding()/100)%10;
+            G4int RightQuark1= string->GetRightParton()->GetPDGEncoding()/1000;
+            G4int RightQuark2=(string->GetRightParton()->GetPDGEncoding()/100)%10;
+           if(G4UniformRand()<0.5)
+           {
+              LeftHadron =hadronizer->Build(FindParticle( LiftQuark1),
+                                            FindParticle(RightQuark1));
+              RightHadron=hadronizer->Build(FindParticle( LiftQuark2),
+                                            FindParticle(RightQuark2));
+           } else
+           {
+              LeftHadron =hadronizer->Build(FindParticle( LiftQuark1),
+                                            FindParticle(RightQuark2));
+              RightHadron=hadronizer->Build(FindParticle( LiftQuark2),
+                                            FindParticle(RightQuark1));
+           }
+        }
+/*
+G4cout<<"SplitLast Left Right hadrons "<<LeftHadron->GetPDGEncoding()<<" "<<RightHadron->GetPDGEncoding()<<G4endl;
+G4cout<<"SplitLast Left Right hadrons "<<LeftHadron->GetPDGMass()<<" "<<RightHadron->GetPDGMass()<<G4endl;
+G4cout<<"Sum H mass Str Mass "<<LeftHadron->GetPDGMass() + RightHadron->GetPDGMass()<<" "<<ResidualMass<<G4endl;
+*/
+       //... repeat procedure, if mass of cluster is too low to produce hadrons
+       //... ClusterMassCut = 0.15*GeV model parameter
+    }
+    while (ResidualMass <= LeftHadron->GetPDGMass() + RightHadron->GetPDGMass());// UzhiVOVA
+    //... compute hadron momenta and energies
+    G4LorentzVector  LeftMom, RightMom;
+    G4ThreeVector    Pos;
+    Sample4Momentum(&LeftMom,  LeftHadron->GetPDGMass(),
+                    &RightMom, RightHadron->GetPDGMass(),
+                    ResidualMass);
+    LeftMom.boost(ClusterVel);
+    RightMom.boost(ClusterVel);
+    LeftVector->push_back(new G4KineticTrack(LeftHadron, 0, Pos, LeftMom));
+    RightVector->push_back(new G4KineticTrack(RightHadron, 0, Pos, RightMom));
+//G4cout<<"Out SplitLast "<<G4endl;
+    return true;
+}
 //----------------------------------------------------------------------------------------------------------
+//G4double G4LundStringFragmentation::GetLightConeZ(G4double zmin, G4double zmax,           // Uzhi
+//                                                  G4int ,  G4ParticleDefinition* pHadron, // Uzhi
+G4double G4LundStringFragmentation::GetLightConeZ(G4double zmin, G4double zmax,
+                                                  G4int,  G4ParticleDefinition* pHadron, // Uzhi
+G4double Px, G4double Py)
+void G4LundStringFragmentation::Sample4Momentum(G4LorentzVector* Mom, G4double Mass, G4LorentzVector* AntiMom, G4double AntiMass, G4double InitialMass)
+  {
+// ------ Sampling of momenta of 2 last produced hadrons --------------------
+     G4ThreeVector Pt;
+     G4double MassMt2, AntiMassMt2;
+     G4double AvailablePz, AvailablePz2;
+//G4cout<<"Sample4Momentum "<<G4endl;
+//G4cout<<"Sample4Momentum Mass"<<Mass<<" "<<AntiMass<<" "<<InitialMass<<G4endl;
+    if(Mass > 930. || AntiMass > 930.)              // If there is a baryon
+    {
+    const G4double  alund = 0.7/GeV/GeV;
+//    If blund get restored, you MUST adapt the calculation of zOfMaxyf.
+//    const G4double  blund = 1;
+    G4double z, yf;
+    G4double Mass = pHadron->GetPDGMass();
+     // ----------------- Isotripic decay ------------------------------------
+       G4double r_val = sqr(InitialMass*InitialMass - Mass*Mass - AntiMass*AntiMass) -
+                          sqr(2.*Mass*AntiMass);
+       G4double Pabs = (r_val > 0.)? std::sqrt(r_val)/(2.*InitialMass) : 0;
+     //... sample unit vector
+       G4double pz = 1. - 2.*G4UniformRand();
+       G4double st     = std::sqrt(1. - pz * pz)*Pabs;
+       G4double phi    = 2.*pi*G4UniformRand();
+       G4double px = st*std::cos(phi);
+       G4double py = st*std::sin(phi);
+       pz *= Pabs;
+    G4double Mt2 = Px*Px + Py*Py + Mass*Mass;
+    G4double zOfMaxyf=alund*Mt2/(alund*Mt2 + 1.);
+    G4double maxYf=(1-zOfMaxyf)/zOfMaxyf * std::exp(-alund*Mt2/zOfMaxyf);
+//    G4double N=1.;                                                 // Uzhi
+//    G4double OverN=1./N;                                           // Uzhi
+//    G4double ZminN=std::pow(zmin,N);                               // Uzhi
+//    G4double ZmaxN=std::pow(zmax,N);                               // Uzhi
+//    G4double Brac=ZmaxN-ZminN;                                     // Uzhi
+//G4cout<<" ZminN ZmaxN Brac Code "<<ZminN<<" "<< ZmaxN<<" "<<Brac<<" "<<PartonEncoding<<G4endl;
+//    if(std::abs(PartonEncoding) < 1000)                            // Uzhi
+      {                                                            // Uzhi q or q-bar
+//G4cout<<" quark "<<G4endl; // Vova
+       do                                                          // Uzhi
+         {
+          z = zmin + G4UniformRand()*(zmax-zmin);
+//        yf = std::pow(1. - z, blund)/z*std::exp(-alund*Mt2/z);
+          yf = (1-z)/z * std::exp(-alund*Mt2/z);
+         }
+       while (G4UniformRand()*maxYf > yf);
+      }                                                            // Uzhi
+//    else                                                           // Uzhi
+//      {                                                            // Uzhi qq or qq-bar
+// //G4cout<<"Di-quark"<<G4endl; // Vova
+//       z = std::pow(Brac * G4UniformRand() + ZminN, OverN);        // Uzhi
+//      };                                                           // Uzhi
+//
+//G4cout<<" test z "<<std::pow(2.,3.)<<" "<<z<<G4endl; // Vova
+    return z;
+       Mom->setPx(px); Mom->setPy(py); Mom->setPz(pz);
+       Mom->setE(std::sqrt(Pabs*Pabs + Mass*Mass));
+       AntiMom->setPx(-px); AntiMom->setPy(-py); AntiMom->setPz(-pz);
+       AntiMom->setE (std::sqrt(Pabs*Pabs + AntiMass*AntiMass));
+    }
+    else
+   {
+      do
+      {
+         Pt=SampleQuarkPt(); Pt.setZ(0); G4double Pt2=Pt.mag2();
+//G4cout<<"Sample4Momentum Pt x y "<<Pt.getX()<<" "<<Pt.getY()<<G4endl;
+         MassMt2    =     Mass *     Mass + Pt2;
+         AntiMassMt2= AntiMass * AntiMass + Pt2;
+//G4cout<<"Mts "<<MassMt2<<" "<<AntiMassMt2<<" "<<InitialMass*InitialMass<<G4endl;
+         AvailablePz2= sqr(InitialMass*InitialMass - MassMt2 - AntiMassMt2) -
+.*MassMt2*AntiMassMt2;
+      }
+      while(AvailablePz2 < 0.);
+      AvailablePz2 /=(4.*InitialMass*InitialMass);
+      AvailablePz = std::sqrt(AvailablePz2);
+//G4cout<<"AvailablePz "<<AvailablePz<<G4endl;
+      G4double Px=Pt.getX();
+      G4double Py=Pt.getY();
+//if(Mass > AntiMass){AvailablePz=-AvailablePz;}  // May30                   // Uzhi
+      Mom->setPx(Px); Mom->setPy(Py); Mom->setPz(AvailablePz);
+      Mom->setE(std::sqrt(MassMt2+AvailablePz2));
+//G4cout<<" 1 part "<<Px<<"  "<<Py<<" "<<AvailablePz<<" "<<std::sqrt(MassMt2+AvailablePz2)<<G4endl;
+     AntiMom->setPx(-Px); AntiMom->setPy(-Py); AntiMom->setPz(-AvailablePz);
+     AntiMom->setE (std::sqrt(AntiMassMt2+AvailablePz2));
+//G4cout<<" 2 part "<<-Px<<"  "<<-Py<<" "<<-AvailablePz<<" "<<std::sqrt(AntiMassMt2+AvailablePz2)<<G4endl;
+    }
+//-----------------------------------------------------------------------------------------
+//G4cout<<"Out Sample4Momentum "<<G4endl;
+  }
+//-----------------------------------------------------------------------------
 G4LorentzVector * G4LundStringFragmentation::SplitEandP(G4ParticleDefinition * pHadron,
+        G4FragmentingString * string)
+{
+        G4FragmentingString * string, G4FragmentingString * newString)
+{
+/*
+G4cout<<"SplitEandP "<<G4endl;
+G4cout<<"SplitEandP string mass "<<string->Mass()<<G4endl;
+G4cout<<string->GetLeftParton()->GetPDGEncoding()<<" "
+      <<string->GetRightParton()->GetPDGEncoding()<<" "<<G4endl;
+G4cout<<G4endl;
+G4cout<<newString->GetLeftParton()->GetPDGEncoding()<<" "<<G4endl;
+G4cout<<newString->GetRightParton()->GetPDGEncoding()<<" "<<G4endl;
+*/
+       G4LorentzVector String4Momentum=string->Get4Momentum();
+       G4double StringMT2=string->Get4Momentum().mt2();
+//G4cout<<"StringMt2 "<<StringMT2<<G4endl;
        G4double HadronMass = pHadron->GetPDGMass();
+       SetMinimalStringMass(string);                                             // Uzhi
+       G4double  StringMass2 = string->Mass2();                                  // Uzhi
+//G4cout<<"SplitEandP string mass "<<string->Mass()<<" Hadron mass "<<HadronMass<<pHadron->GetParticleName()<<G4endl;   // Uzhi
 //G4cout<<string->GetLeftParton()->GetPDGEncoding()<<" "<<G4endl;
 //G4cout<<string->GetRightParton()->GetPDGEncoding()<<" "<<G4endl;
 //G4cout<<" Min string mass "<<MinimalStringMass<<G4endl;
        // calculate and assign hadron transverse momentum component HadronPx andHadronPy
+//G4cout<<"Hadron mass "<<HadronMass<<" "<<pHadron->GetParticleName()<<G4endl;
+       SetMinimalStringMass(newString);
+       String4Momentum.setPz(0.);
+       G4ThreeVector StringPt=String4Momentum.vect();
+//G4cout<<"StringPt "<<StringPt<<G4endl<<G4endl;
+//G4cout<<"Min string mass "<<MinimalStringMass<<G4endl;
+// calculate and assign hadron transverse momentum component HadronPx and HadronPy
        G4ThreeVector thePt;
        thePt=SampleQuarkPt();
 …
        G4ThreeVector HadronPt = thePt +string->DecayPt();
        HadronPt.setZ(0);
+//G4cout<<"Hadron Pt"<<HadronPt<<G4endl;
+       G4ThreeVector RemSysPt = StringPt - HadronPt;
+//G4cout<<"RemSys Pt"<<RemSysPt<<G4endl;
        //...  sample z to define hadron longitudinal momentum and energy
        //... but first check the available phase space
+//       G4double DecayQuarkMass2  = sqr(string->GetDecayParton()->GetPDGMass());
+//G4cout<<" QuarkMass "<<string->GetDecayParton()->GetPDGMass()<<G4endl;              // Uzhi
+       G4double HadronMass2T = sqr(HadronMass) + HadronPt.mag2();
+//       G4double ResidualMass2T=sqr(MinimalStringMass + WminLUND) + HadronPt.mag2(); // Uzhi
+       G4double ResidualMass2T=sqr(MinimalStringMass + WminLUND) + HadronPt.mag2(); // Uzhi
+//G4cout<<" Mt h res str "<<std::sqrt(HadronMass2T)<<" "<<std::sqrt(ResidualMass2T)<<" srt mass"<<string->Mass()<<G4endl;
+//       if (DecayQuarkMass2 + HadronMass2T >= SmoothParam*(string->Mass2()) )      // Uzhi
+        G4double Pz2 = (sqr(StringMass2 - HadronMass2T - ResidualMass2T) -                  // Uzhi
+*HadronMass2T * ResidualMass2T)/4./StringMass2; // Uzhi
+//G4cout<<" Pz**2 "<<Pz2<<G4endl;
+        if(Pz2 < 0 ) {return 0;}          // have to start all over!                // Uzhi
+       G4double HadronMassT2 = sqr(HadronMass) + HadronPt.mag2();
+       G4double ResidualMassT2=sqr(MinimalStringMass) + RemSysPt.mag2();
+//G4cout<<"Mt h res str "<<std::sqrt(HadronMassT2)<<" "<<std::sqrt(ResidualMassT2)<<" srt mass"<<StringMT2<<G4endl;
+        G4double Pz2 = (sqr(StringMT2 - HadronMassT2 - ResidualMassT2) -
+*HadronMassT2 * ResidualMassT2)/4./StringMT2;
+//G4cout<<"Pz**2 "<<Pz2<<G4endl;
+        if(Pz2 < 0 ) {return 0;}          // have to start all over!
        //... then compute allowed z region  z_min <= z <= z_max
+       G4double Pz = std::sqrt(Pz2);                                               // Uzhi
+       G4double zMin = (std::sqrt(HadronMass2T+Pz2) - Pz)/std::sqrt(StringMass2);  // Uzhi
+       G4double zMax = (std::sqrt(HadronMass2T+Pz2) + Pz)/std::sqrt(StringMass2);  // Uzhi
+//G4cout<<" Zmin max "<<zMin<<" "<<zMax<<G4endl;               // Uzhi
+//       G4double zMax = 1. - DecayQuarkMass2/(string->Mass2());                   // Uzhi
+       G4double Pz = std::sqrt(Pz2);
+       G4double zMin = (std::sqrt(HadronMassT2+Pz2) - Pz)/std::sqrt(StringMT2);
+       G4double zMax = (std::sqrt(HadronMassT2+Pz2) + Pz)/std::sqrt(StringMT2);
+//G4cout<<"Zmin max "<<zMin<<" "<<zMax<<G4endl;               // Uzhi
        if (zMin >= zMax) return 0;              // have to start all over!
 …
         HadronPt.setZ(0.5* string->GetDecayDirection() *
                         (z * string->LightConeDecay() -
                          HadronMass2T/(z * string->LightConeDecay())));
+                         HadronMassT2/(z * string->LightConeDecay())));
         G4double HadronE  = 0.5* (z * string->LightConeDecay() +
                                   HadronMass2T/(z * string->LightConeDecay()));
+                                  HadronMassT2/(z * string->LightConeDecay()));
        G4LorentzVector * a4Momentum= new G4LorentzVector(HadronPt,HadronE);
+//G4cout<<"Out of SplitEandP Pz E "<<HadronPt.getZ()<<" "<<0.5* (z * string->LightConeDecay() + HadronMass2T/(z * string->LightConeDecay()))<<G4endl;
+//G4cout<<"Hadron Pt"<<HadronPt<<G4endl;
+//G4cout<<"Out of SplitEandP Pz E "<<HadronPt.getZ()<<" "<<HadronE<<G4endl;
        return a4Momentum;
+}
 //-----------------------------------------------------------------------------------------
 G4bool G4LundStringFragmentation::SplitLast(G4FragmentingString * string,
                                              G4KineticTrackVector * LeftVector,
                                              G4KineticTrackVector * RightVector)
+G4double G4LundStringFragmentation::GetLightConeZ(G4double zmin, G4double zmax,
+                                                  G4int PDGEncodingOfDecayParton,
+                                                  G4ParticleDefinition* pHadron,
+                                                  G4double Px, G4double Py)
+{
+    //... perform last cluster decay
+//G4cout<<"SplitLast String mass "<<string->Mass()<<G4endl;
+//G4cout<<string->GetLeftParton()->GetPDGEncoding()<<" "<<G4endl;
+//G4cout<<string->GetRightParton()->GetPDGEncoding()<<" "<<G4endl;
+    G4ThreeVector ClusterVel =string->Get4Momentum().boostVector();
+    G4double ResidualMass   = string->Mass();
+//    G4double ClusterMassCut = ClusterMass;
+    G4int cClusterInterrupt = 0;
+    G4ParticleDefinition * LeftHadron, * RightHadron;
+    do
+    {
+//G4cout<<" Cicle "<<cClusterInterrupt<<" "<< ClusterLoopInterrupt<<G4endl;
+        if (cClusterInterrupt++ >= ClusterLoopInterrupt)
+        {
+          return false;
+        }
+        G4ParticleDefinition * quark = NULL;
+        string->SetLeftPartonStable(); // to query quark contents..
+        if (string->DecayIsQuark() && string->StableIsQuark() )
+        {
+           //... there are quarks on cluster ends
+                LeftHadron= QuarkSplitup(string->GetLeftParton(), quark);
+        } else {
+           //... there is a Diquark on cluster ends
+                G4int IsParticle;
+                if ( string->StableIsQuark() ) {
+                  IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? -1 : +1;
+                } else {
+                  IsParticle=(string->GetLeftParton()->GetPDGEncoding()>0) ? +1 : -1;
+                }
+                pDefPair QuarkPair = CreatePartonPair(IsParticle,false);  // no diquarks wanted
+                quark = QuarkPair.second;
+                LeftHadron=hadronizer->Build(QuarkPair.first, string->GetLeftParton());
+        }
+        RightHadron = hadronizer->Build(string->GetRightParton(), quark);
+//G4cout<<"SplitLast Left Right hadrons "<<LeftHadron->GetPDGEncoding()<<" "<<RightHadron->GetPDGEncoding()<<G4endl;
+//G4cout<<"SplitLast Left Right hadrons "<<LeftHadron->GetPDGMass()<<" "<<RightHadron->GetPDGMass()<<G4endl;
+//G4cout<<" Sum H mass Str Mass "<<LeftHadron->GetPDGMass() + RightHadron->GetPDGMass()<<" "<<ResidualMass<<G4endl;
+       //... repeat procedure, if mass of cluster is too low to produce hadrons
+       //... ClusterMassCut = 0.15*GeV model parameter
+//      if ( quark->GetParticleSubType()== "quark" ) {ClusterMassCut = 0.;}        // Uzhi
+//      else {ClusterMassCut = ClusterMass;}                                       // Uzhi
+    }
+    while (ResidualMass <= LeftHadron->GetPDGMass() + RightHadron->GetPDGMass());         // Uzhi   VOVA
+//    while (ResidualMass <= LeftHadron->GetPDGMass() + RightHadron->GetPDGMass()  + ClusterMassCut); // Uzhi
+    //... compute hadron momenta and energies
+    G4LorentzVector  LeftMom, RightMom;
+    G4ThreeVector    Pos;
+//G4cout<<"Sample4Momentum"<<G4endl;
+    Sample4Momentum(&LeftMom, LeftHadron->GetPDGMass(), &RightMom, RightHadron->GetPDGMass(), ResidualMass);
+    LeftMom.boost(ClusterVel);
+    RightMom.boost(ClusterVel);
+    LeftVector->push_back(new G4KineticTrack(LeftHadron, 0, Pos, LeftMom));
+    RightVector->push_back(new G4KineticTrack(RightHadron, 0, Pos, RightMom));
+    return true;
+}
+//----------------------------------------------------------------------------------------------------------
+G4bool G4LundStringFragmentation::IsFragmentable(const G4FragmentingString * const string)
+{
+//G4cout<<"In IsFragmentable"<<G4endl;
+  SetMinimalStringMass(string);                                                            // Uzhi
+//G4cout<<"Out IsFragmentable MinMass"<<MinimalStringMass<<" String Mass"<<std::sqrt(string->Get4Momentum().mag2())<<G4endl;
+  return sqr(MinimalStringMass + WminLUND) < string->Get4Momentum().mag2();                // Uzhi
+//      return sqr(FragmentationMass(string)+MassCut) <                                    // Uzhi
+//                      string->Mass2();                                                   // Uzhi
+}
+//----------------------------------------------------------------------------------------------------------
+G4bool G4LundStringFragmentation::StopFragmenting(const G4FragmentingString * const string)
+{
+//G4cout<<"StopFragmenting"<<G4endl;
+  SetMinimalStringMass(string);                                                            // Uzhi
+//G4cout<<"StopFragm MinMass "<<MinimalStringMass<<" String Mass "<<std::sqrt(string->Get4Momentum().mag2())<<G4endl;
+  return sqr((MinimalStringMass + WminLUND)*(1 + SmoothParam * (1.-2*G4UniformRand()))) >        // Uzhi
+                   string->Get4Momentum().mag2();                                          // Uzhi
+//       sqr(FragmentationMass(string,&G4HadronBuilder::BuildHighSpin)+MassCut) >          // Uzhi
+//       string->Get4Momentum().mag2();                                                    // Uzhi
+}
+//----------------------------------------------------------------------------------------------------------
+void G4LundStringFragmentation::Sample4Momentum(G4LorentzVector* Mom, G4double Mass, G4LorentzVector* AntiMom, G4double AntiMass, G4double InitialMass)
+    {
+     G4ThreeVector Pt;                                                      // Uzhi
+     G4double MassMt2, AntiMassMt2;                                         // Uzhi
+     G4double AvailablePz, AvailablePz2;                                    // Uzhi
+//G4cout<<" Smpl4Mom "<<Mass<<" "<<AntiMass<<" "<<InitialMass<<G4endl;
+                                                                            // Uzhi
+    do                                                                      // Uzhi
+      {                                                                     // Uzhi
+       Pt=SampleQuarkPt(); Pt.setZ(0); G4double Pt2=Pt.mag2();              // Uzhi
+//G4cout<<"Sample4Momentum Pt x y "<<Pt.getX()<<" "<<Pt.getY()<<G4endl;
+       MassMt2    =     Mass *     Mass + Pt2;                              // Uzhi
+       AntiMassMt2= AntiMass * AntiMass + Pt2;                              // Uzhi
+//G4cout<<"Mts "<<MassMt2<<" "<<AntiMassMt2<<" "<<InitialMass*InitialMass<<G4endl;
+       AvailablePz2= sqr(InitialMass*InitialMass - MassMt2 - AntiMassMt2) -
+.*MassMt2*AntiMassMt2;                                // Uzhi
+      }                                                                     // Uzhi
+    while(AvailablePz2 < 0.);                                               // Uzhi
+                                                                            // Uzhi
+    AvailablePz2 /=(4.*InitialMass*InitialMass);                            // Uzhi
+                                                                            // Uzhi
+    AvailablePz = std::sqrt(AvailablePz2);                                  // Uzhi
+//G4cout<<"AvailablePz "<<AvailablePz<<G4endl;
+    G4double Px=Pt.getX();                                                  // Uzhi
+    G4double Py=Pt.getY();                                                  // Uzhi
+                                                                            // Uzhi
+    Mom->setPx(Px); Mom->setPy(Py); Mom->setPz(AvailablePz);                // Uzhi
+    Mom->setE(std::sqrt(MassMt2+AvailablePz2));                             // Uzhi
+//G4cout<<" 1 part "<<Px<<"  "<<Py<<" "<<AvailablePz<<" "<<std::sqrt(MassMt2+AvailablePz2)<<G4endl;
+                                                                            // Uzhi
+    AntiMom->setPx(-Px); AntiMom->setPy(-Py); AntiMom->setPz(-AvailablePz); // Uzhi
+    AntiMom->setE (std::sqrt(AntiMassMt2+AvailablePz2));                    // Uzhi
+//G4cout<<" 2 part "<<-Px<<"  "<<-Py<<" "<<-AvailablePz<<" "<<std::sqrt(AntiMassMt2+AvailablePz2)<<G4endl;
+// Maybe it must be inversed!                                               // Uzhi
+/*                                                                          // Uzhi
+    G4double r_val = sqr(InitialMass*InitialMass - Mass*Mass - AntiMass*AntiMass) - sqr(2.*Mass*AntiMass);
+    G4double Pabs = (r_val > 0.)? std::sqrt(r_val)/(2.*InitialMass) : 0;
+    //... sample unit vector
+    G4double pz = 1. - 2.*G4UniformRand();
+    G4double st     = std::sqrt(1. - pz * pz)*Pabs;
+    G4double phi    = 2.*pi*G4UniformRand();
+    G4double px = st*std::cos(phi);
+    G4double py = st*std::sin(phi);
+    pz *= Pabs;
+    G4double  alund;
+//    If blund get restored, you MUST adapt the calculation of zOfMaxyf.
+//    const G4double  blund = 1;
+    G4double z, yf;
+    G4double Mass = pHadron->GetPDGMass();
+//  G4int HadronEncoding=pHadron->GetPDGEncoding();
+    Mom->setPx(px); Mom->setPy(py); Mom->setPz(pz);
+    Mom->setE(std::sqrt(Pabs*Pabs + Mass*Mass));
+    AntiMom->setPx(-px); AntiMom->setPy(-py); AntiMom->setPz(-pz);
+    AntiMom->setE (std::sqrt(Pabs*Pabs + AntiMass*AntiMass));
+*/                                                                           // Uzhi
+    G4double Mt2 = Px*Px + Py*Py + Mass*Mass;
+    if(std::abs(PDGEncodingOfDecayParton) < 1000)
+    {
+    // ---------------- Quark fragmentation ----------------------
+       alund=0.35/GeV/GeV; // Instead of 0.7 because kinks are not considered
+       G4double zOfMaxyf=alund*Mt2/(alund*Mt2 + 1.);
+       G4double maxYf=(1-zOfMaxyf)/zOfMaxyf * std::exp(-alund*Mt2/zOfMaxyf);
+       do
+         {
+          z = zmin + G4UniformRand()*(zmax-zmin);
+//        yf = std::pow(1. - z, blund)/z*std::exp(-alund*Mt2/z);
+          yf = (1-z)/z * std::exp(-alund*Mt2/z);
+         }
+       while (G4UniformRand()*maxYf > yf);
+    }
+    else
+    {
+     // ---------------- Di-quark fragmentation ----------------------
+//G4cout<<"Di-quark"<<G4endl; // Vova
+       alund=0.7/GeV/GeV;    // 0.7 2.0
+       G4double zOfMaxyf=alund*Mt2/(alund*Mt2 + 1.);
+       G4double maxYf=(1-zOfMaxyf)/zOfMaxyf * std::exp(-alund*Mt2/zOfMaxyf);
+       do
+         {
+          z = zmin + G4UniformRand()*(zmax-zmin);
+//        yf = std::pow(1. - z, blund)/z*std::exp(-alund*Mt2/z);
+          yf = (1-z)/z * std::exp(-alund*Mt2/z);
+         }
+       while (G4UniformRand()*maxYf > yf);
+    };
+//G4cout<<" test z "<<std::pow(2.,3.)<<" "<<z<<G4endl;
+    return z;
+    }
-void G4LundStringFragmentation::SetMinimalStringMass(const G4FragmentingString  * const string)  // Uzhi
+{
-//G4cout<<"In SetMinMass -------------------"<<std::sqrt(string->Mass2())<<G4endl;
-//G4cout<<string->GetLeftParton()->GetPDGEncoding()<<G4endl;
-//G4cout<<string->GetRightParton()->GetPDGEncoding()<<G4endl;
-  G4double EstimatedMass=0.750* GeV;  // 2*m_q
-  G4int Qleft =std::abs(string->GetLeftParton()->GetPDGEncoding());
-  if( Qleft > 1000)
+    {
-     G4int q1=Qleft/1000;
-     if( q1 < 3) {EstimatedMass += 0.325* GeV;}
-     if( q1 > 2) {EstimatedMass += 0.500* GeV;}
-     G4int q2=(Qleft/100)%10;
-     if( q2 < 3) {EstimatedMass += 0.325* GeV;}
-     if( q2 > 2) {EstimatedMass += 0.500* GeV;}
+    }
-  else
+    {
-     if( Qleft < 3) {EstimatedMass += 0.325* GeV;}
-     if( Qleft > 2) {EstimatedMass += 0.500* GeV;}
+    }
-  G4int Qright=std::abs(string->GetRightParton()->GetPDGEncoding());
-  if( Qright > 1000)
+    {
-     G4int q1=Qright/1000;
-     if( q1 < 3) {EstimatedMass += 0.325* GeV;}
-     if( q1 > 2) {EstimatedMass += 0.500* GeV;}
-     G4int q2=(Qright/100)%10;
-     if( q2 < 3) {EstimatedMass += 0.325* GeV;}
-     if( q2 > 2) {EstimatedMass += 0.500* GeV;}
+    }
-  else
+    {
-     if( Qright < 3) {EstimatedMass += 0.325* GeV;}
-     if( Qright > 2) {EstimatedMass += 0.500* GeV;}
+    }
-  MinimalStringMass=EstimatedMass;
-  SetMinimalStringMass2(EstimatedMass);
-/*
-  Pcreate build=&G4HadronBuilder::BuildLowSpin;
-  G4ParticleDefinition *Hadron1, *Hadron2=0;
-  G4int iflc = (G4UniformRand() < 0.5)? 1 : 2;
-  if (string->GetLeftParton()->GetParticleSubType() == "quark") iflc = -iflc;
-  if (string->GetLeftParton()->GetPDGEncoding() < 0) iflc = -iflc;
-  // 1/2 baryon (anti-baryon) and scalar meson     (QQ-q or QbarQbar-Qbar),
-  // or 2 scalar mesons                            (Q-Qbar),
-  // or 2 1/2 baryons (anti-baryons) will be built (QQ-QbarQbar)
-//G4cout<<"In SetMinMass -------------------"<<std::sqrt(string->Mass2())<<G4endl;
-//G4cout<<string->GetLeftParton()->GetPDGEncoding()<<" "<<FindParticle(iflc)->GetPDGEncoding()<<G4endl;
-//G4cout<<string->GetRightParton()->GetPDGEncoding()<<" "<<FindParticle(-iflc)->GetPDGEncoding()<<G4endl;
-  Hadron1 = (hadronizer->*build)(string->GetLeftParton(),FindParticle(iflc));
-  Hadron2 =(hadronizer->*build)(string->GetRightParton(),FindParticle(-iflc));
-  MinimalStringMass = (Hadron1)->GetPDGMass() + (Hadron2)->GetPDGMass();
-//G4cout<<(Hadron1)->GetPDGEncoding()<<" "<<(Hadron2)->GetPDGEncoding()<<G4endl;
-//G4cout<<"Out SetMinMass "<<MinimalStringMass<<G4endl;
-*/
-//  SetMinimalStringMass2(MinimalStringMass);
+}
-//*******************************************************************************************************
-void G4LundStringFragmentation::SetMinimalStringMass2(const G4double aValue)                     // Uzhi
+{
-  MinimalStringMass2=aValue * aValue;
+}
-//*******************************************************************************************************
-//****************************************************************************************

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

-                      r819
+                      r962
 //
 //
 // $Id: G4QGSMFragmentation.cc,v 1.6 2007/04/24 14:55:23 gunter Exp $
 // GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4QGSMFragmentation.cc,v 1.9 2008/06/23 08:35:55 vuzhinsk Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // -----------------------------------------------------------------------------
 …
                         } else {
                          // abandon ... start from the beginning
                            if (newString) delete newString;
+                           if (newString) delete newString;          // Uzhi restore 20.06.08
                            if (Hadron)    delete Hadron;
                            inner_sucess=false;
 …
 G4LorentzVector * G4QGSMFragmentation::SplitEandP(G4ParticleDefinition * pHadron,
+        G4FragmentingString * string)
+                                                  G4FragmentingString * string,    // Uzhi
+                                                  G4FragmentingString * ) // Uzhi
+{
        G4double HadronMass = pHadron->GetPDGMass();

trunk/source/processes/hadronic/models/parton_string/hadronization/src/G4VKinkyStringDecay.cc

-                      r819
+                      r962
 //
 //
 // $Id: G4VKinkyStringDecay.cc,v 1.3 2006/06/29 20:55:07 gunter Exp $
 // GEANT4 tag $Name:  $
+// $Id: G4VKinkyStringDecay.cc,v 1.4 2008/04/25 14:20:14 vuzhinsk Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //  Maxim Komogorov
 //

trunk/source/processes/hadronic/models/parton_string/hadronization/src/G4VLongitudinalStringDecay.cc

-                      r819
+                      r962
 //
 //
 // $Id: G4VLongitudinalStringDecay.cc,v 1.8 2007/04/24 14:55:23 gunter Exp $
 // GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4VLongitudinalStringDecay.cc,v 1.13 2008/06/23 08:35:55 vuzhinsk Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // -----------------------------------------------------------------------------
 …
    StrangeSuppress  = 0.44;    //  27 % strange quarks produced, ie. u:d:s=1:1:0.27
    DiquarkSuppress  = 0.1;
+   DiquarkSuppress  = 0.07;
    DiquarkBreakProb = 0.1;
 …
    hadronizer = new G4HadronBuilder(pspin_meson,pspin_barion,
                                 scalarMesonMix,vectorMesonMix);
+   Kappa = 1.0 * GeV/fermi;
+}
 …
+   }
 //=============================================================================================-------------
+//=============================================================================
 // Operators
 …
 //    }
 //----------------------------------------------------------------------------------------------------------
+//-----------------------------------------------------------------------------
 int G4VLongitudinalStringDecay::operator==(const G4VLongitudinalStringDecay &) const
 …
+    }
 //----------------------------------------------------------------------------------------------------------
+//-------------------------------------------------------------------------------------
 int G4VLongitudinalStringDecay::operator!=(const G4VLongitudinalStringDecay &) const
 …
+    }
+//==========================================================================================================
+//***********************************************************************************
+// For changing Mass Cut used for selection of very small mass strings
+void G4VLongitudinalStringDecay::SetMassCut(G4double aValue){MassCut=aValue;}
+//-----------------------------------------------------------------------------
+// For handling a string with very low mass
+G4KineticTrackVector* G4VLongitudinalStringDecay::LightFragmentationTest(const
+                G4ExcitedString * const string)
+{
+   // Check string decay threshold
+        G4KineticTrackVector * result=0;  // return 0 when string exceeds the mass cut
+        pDefPair hadrons((G4ParticleDefinition *)0,(G4ParticleDefinition *)0);
+        G4FragmentingString aString(*string);
+        if ( sqr(FragmentationMass(&aString,0,&hadrons)+MassCut) < aString.Mass2()) {
+                return 0;
+        }
+// The string mass is very low ---------------------------
+        result=new G4KineticTrackVector;
+        if ( hadrons.second ==0 )
+        {
+// Substitute string by light hadron, Note that Energy is not conserved here!
+/*
+#ifdef DEBUG_LightFragmentationTest
+               G4cout << "VlongSF Warning replacing string by single hadron " <<G4endl;
+               G4cout << hadrons.first->GetParticleName()
+                      << "string .. " << string->Get4Momentum() << " "
+                      << string->Get4Momentum().m() << G4endl;
+#endif
+G4cout << "VlongSF Warning replacing string by single hadron " <<G4endl;
+G4cout << hadrons.first->GetParticleName() << " string .. " <<G4endl;
+G4cout << string->Get4Momentum() << " " << string->Get4Momentum().m() << G4endl;
+*/
+               G4ThreeVector   Mom3 = string->Get4Momentum().vect();
+               G4LorentzVector Mom(Mom3,
+                                   std::sqrt(Mom3.mag2() +
+                                             sqr(hadrons.first->GetPDGMass())));
+               result->push_back(new G4KineticTrack(hadrons.first, 0,
+                                                  string->GetPosition(),
+                                                          Mom));
+        } else
+        {
+//... string was qq--qqbar type: Build two stable hadrons,
+#ifdef DEBUG_LightFragmentationTest
+               G4cout << "VlongSF Warning replacing qq-qqbar string by TWO hadrons "
+                      << hadrons.first->GetParticleName() << " / "
+                      << hadrons.second->GetParticleName()
+                      << "string .. " << string->Get4Momentum() << " "
+                      << string->Get4Momentum().m() << G4endl;
+#endif
+// Uzhi Formation time in the case???
+               G4LorentzVector  Mom1, Mom2;
+               Sample4Momentum(&Mom1, hadrons.first->GetPDGMass(),
+                               &Mom2,hadrons.second->GetPDGMass(),
+                               string->Get4Momentum().mag());
+               result->push_back(new G4KineticTrack(hadrons.first, 0,
+                                                    string->GetPosition(),
+                                                            Mom1));
+               result->push_back(new G4KineticTrack(hadrons.second, 0,
+                                                    string->GetPosition(),
+                                                    Mom2));
+               G4ThreeVector Velocity = string->Get4Momentum().boostVector();
+               result->Boost(Velocity);
+        }
+        return result;
+}
+//----------------------------------------------------------------------------------------
+G4double G4VLongitudinalStringDecay::FragmentationMass(
+            const G4FragmentingString * const string,
+                Pcreate build, pDefPair * pdefs       )
+{
+        G4double mass;
+        static G4bool NeedInit(true);
+        static std::vector<double> nomix;
+        static G4HadronBuilder * minMassHadronizer;
+        if ( NeedInit )
+        {
+           NeedInit = false;
+           nomix.resize(6);
+           for ( G4int i=0; i<6 ; i++ ) nomix[i]=0;
+//         minMassHadronizer=new G4HadronBuilder(pspin_meson,pspin_barion,nomix,nomix);
+           minMassHadronizer=hadronizer;
+        }
+        if ( build==0 ) build=&G4HadronBuilder::BuildLowSpin;
+        G4ParticleDefinition *Hadron1, *Hadron2=0;
+        if (!string->FourQuarkString() )
+        {
+           // spin 0 meson or spin 1/2 barion will be built
+           Hadron1 = (minMassHadronizer->*build)(string->GetLeftParton(),
+                                                 string->GetRightParton());
+           mass= (Hadron1)->GetPDGMass();
+        } else
+        {
+           //... string is qq--qqbar: Build two stable hadrons,
+           //...    with extra uubar or ddbar quark pair
+           G4int iflc = (G4UniformRand() < 0.5)? 1 : 2;
+           if (string->GetLeftParton()->GetPDGEncoding() < 0) iflc = -iflc;
+           //... theSpin = 4; spin 3/2 baryons will be built
+           Hadron1 = (minMassHadronizer->*build)(string->GetLeftParton(),
+                                                 FindParticle(iflc)       );
+           Hadron2 = (minMassHadronizer->*build)(string->GetRightParton(),
+                                                 FindParticle(-iflc)      );
+           mass = (Hadron1)->GetPDGMass() + (Hadron2)->GetPDGMass();
+        }
+        if ( pdefs != 0 )
+        { // need to return hadrons as well....
+           pdefs->first  = Hadron1;
+           pdefs->second = Hadron2;
+        }
+        return mass;
+}
+//----------------------------------------------------------------------------
+G4ParticleDefinition* G4VLongitudinalStringDecay::FindParticle(G4int Encoding)
+   {
+   G4ParticleDefinition* ptr = G4ParticleTable::GetParticleTable()->FindParticle(Encoding);
+      if (ptr == NULL)
+       {
+       G4cout << "Particle with encoding "<<Encoding<<" does not exist!!!"<<G4endl;
+       throw G4HadronicException(__FILE__, __LINE__, "Check your particle table");
+       }
+   return ptr;
+   }
+//-----------------------------------------------------------------------------
+//   virtual void Sample4Momentum(G4LorentzVector* Mom,     G4double Mass,
+//                                G4LorentzVector* AntiMom, G4double AntiMass,
+//                                G4double InitialMass)=0;
+//-----------------------------------------------------------------------------
+//*********************************************************************************
+//   For decision on continue or stop string fragmentation
+//   virtual G4bool StopFragmenting(const G4FragmentingString  * const string)=0;
+//   virtual G4bool IsFragmentable(const G4FragmentingString * const string)=0;
+//   If a string can not fragment, make last break into 2 hadrons
+//   virtual G4bool SplitLast(G4FragmentingString * string,
+//                            G4KineticTrackVector * LeftVector,
+//                            G4KineticTrackVector * RightVector)=0;
+//-----------------------------------------------------------------------------
+//
+//   If a string fragments, do the following
+//
+//   For transver of a string to its CMS frame
+//-----------------------------------------------------------------------------
+G4ExcitedString *G4VLongitudinalStringDecay::CPExcited(const G4ExcitedString & in)
+{
+        G4Parton *Left=new G4Parton(*in.GetLeftParton());
+        G4Parton *Right=new G4Parton(*in.GetRightParton());
+        return new G4ExcitedString(Left,Right,in.GetDirection());
+}
+//-----------------------------------------------------------------------------
+G4KineticTrack * G4VLongitudinalStringDecay::Splitup(
+                        G4FragmentingString *string,
+                        G4FragmentingString *&newString)
+{
+//G4cout<<"In G4VLong String Dec######################"<<G4endl;
+       //... random choice of string end to use for creating the hadron (decay)
+       SideOfDecay = (G4UniformRand() < 0.5)? 1: -1;
+       if (SideOfDecay < 0)
+       {
+          string->SetLeftPartonStable();
+       } else
+       {
+          string->SetRightPartonStable();
+       }
+       G4ParticleDefinition *newStringEnd;
+       G4ParticleDefinition * HadronDefinition;
+       if (string->DecayIsQuark())
+       {
+           HadronDefinition= QuarkSplitup(string->GetDecayParton(), newStringEnd);
+       } else {
+           HadronDefinition= DiQuarkSplitup(string->GetDecayParton(), newStringEnd);
+       }
+// create new String from old, ie. keep Left and Right order, but replace decay
+       newString=new G4FragmentingString(*string,newStringEnd); // To store possible
+                                                                // quark containt of new string
+       G4LorentzVector* HadronMomentum=SplitEandP(HadronDefinition, string, newString);
+       delete newString; newString=0;                          // Uzhi 20.06.08
+       G4KineticTrack * Hadron =0;
+       if ( HadronMomentum != 0 ) {
+           G4ThreeVector   Pos;
+           Hadron = new G4KineticTrack(HadronDefinition, 0,Pos, *HadronMomentum);
+           newString=new G4FragmentingString(*string,newStringEnd,
+                                        HadronMomentum);
+//G4cout<<"Out G4VLong String Dec######################"<<G4endl;
+//G4cout<<"newString 4Mom"<<newString->Get4Momentum()<<G4endl;
+//G4cout<<"newString Pl  "<<newString->LightConePlus()<<G4endl;
+//G4cout<<"newString Mi  "<<newString->LightConeMinus()<<G4endl;
+//G4cout<<"newString Pts "<<newString->StablePt()<<G4endl;
+//G4cout<<"newString Ptd "<<newString->DecayPt()<<G4endl;
+//G4cout<<"newString M2  "<<newString->Mass2()<<G4endl;
+//G4cout<<"newString M2  "<<newString->Mass()<<G4endl;
+           delete HadronMomentum;
+       }
+       return Hadron;
+}
+//--------------------------------------------------------------------------------------
+G4ParticleDefinition *
+                G4VLongitudinalStringDecay::QuarkSplitup(G4ParticleDefinition*
+                decay, G4ParticleDefinition *&created)
+{
+    G4int IsParticle=(decay->GetPDGEncoding()>0) ? -1 : +1; // if we have a quark,
+                                                            // we need antiquark
+                                                            // (or diquark)
+    pDefPair QuarkPair = CreatePartonPair(IsParticle);
+    created = QuarkPair.second;
+    return hadronizer->Build(QuarkPair.first, decay);
+}
+//-----------------------------------------------------------------------------
+G4ParticleDefinition *G4VLongitudinalStringDecay::DiQuarkSplitup(
+                                                        G4ParticleDefinition* decay,
+                                                        G4ParticleDefinition *&created)
+{
+   //... can Diquark break or not?
+   if (G4UniformRand() < DiquarkBreakProb ){
+   //... Diquark break
+      G4int stableQuarkEncoding = decay->GetPDGEncoding()/1000;
+      G4int decayQuarkEncoding = (decay->GetPDGEncoding()/100)%10;
+      if (G4UniformRand() < 0.5)
+         {
+         G4int Swap = stableQuarkEncoding;
+         stableQuarkEncoding = decayQuarkEncoding;
+         decayQuarkEncoding = Swap;
+         }
+      G4int IsParticle=(decayQuarkEncoding>0) ? -1 : +1;
+                        // if we have a quark, we need antiquark)
+      pDefPair QuarkPair = CreatePartonPair(IsParticle,false);  // no diquarks wanted
+      //... Build new Diquark
+      G4int QuarkEncoding=QuarkPair.second->GetPDGEncoding();
+      G4int i10  = std::max(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
+      G4int i20  = std::min(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
+      G4int spin = (i10 != i20 && G4UniformRand() <= 0.5)? 1 : 3;
+      G4int NewDecayEncoding = -1*IsParticle*(i10 * 1000 + i20 * 100 + spin);
+      created = FindParticle(NewDecayEncoding);
+      G4ParticleDefinition * decayQuark=FindParticle(decayQuarkEncoding);
+      G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decayQuark);
+      return had;
+//      return hadronizer->Build(QuarkPair.first, decayQuark);
+   } else {
+   //... Diquark does not break
+      G4int IsParticle=(decay->GetPDGEncoding()>0) ? +1 : -1;
+                        // if we have a diquark, we need quark)
+      pDefPair QuarkPair = CreatePartonPair(IsParticle,false);  // no diquarks wanted
+      created = QuarkPair.second;
+      G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decay);
+      return had;
+//      return G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decay);
+   }
+}
+//-----------------------------------------------------------------------------
 G4int G4VLongitudinalStringDecay::SampleQuarkFlavor(void)
 …
+   }
 //----------------------------------------------------------------------------------------------------------
+//-----------------------------------------------------------------------------
 G4VLongitudinalStringDecay::pDefPair G4VLongitudinalStringDecay::CreatePartonPair(G4int NeedParticle,G4bool AllowDiquarks)
 …
+}
 //----------------------------------------------------------------------------------------------------------
+//-----------------------------------------------------------------------------
 // G4ThreeVector G4VLongitudinalStringDecay::SampleQuarkPt()
 …
 //    return G4ThreeVector(Pt * std::cos(phi),Pt * std::sin(phi),0);
 //    }
 G4ThreeVector G4VLongitudinalStringDecay::SampleQuarkPt()
+   {
 …
+   }
 //----------------------------------------------------------------------------------------------------------
+//******************************************************************************
 void G4VLongitudinalStringDecay::CalculateHadronTimePosition(G4double theInitialStringMass, G4KineticTrackVector* Hadrons)
+   {
    // `yo-yo` formation time
+   const G4double kappa = 1.0 * GeV/fermi;
+//   const G4double kappa = 1.0 * GeV/fermi/4.;      // Uzhi String tension 1.06.08
+     G4double kappa = GetStringTensionParameter();
+//G4cout<<"Kappa "<<kappa<<G4endl;                   // Uzhi 20.06.08
+//G4int Uzhi; G4cin>>Uzhi;                           // Uzhi 20.06.08
    for(size_t c1 = 0; c1 < Hadrons->size(); c1++)
+      {
 …
+   }
 //----------------------------------------------------------------------------------------------------------
+//-----------------------------------------------------------------------------
 /*
 …
 */
+//*****************************************************************************
+void G4VLongitudinalStringDecay::SetSigmaTransverseMomentum(G4double aValue)
+{
+        if ( PastInitPhase ) {
+                throw G4HadronicException(__FILE__, __LINE__, "4VLongitudinalStringDecay::SetSigmaTransverseMomentum after FragmentString() not allowed");
+        } else {
+                SigmaQT = aValue;
+        }
+}
 //----------------------------------------------------------------------------------------------------------
+G4ParticleDefinition *
+                G4VLongitudinalStringDecay::QuarkSplitup(G4ParticleDefinition*
+                decay, G4ParticleDefinition *&created)
+{
+    G4int IsParticle=(decay->GetPDGEncoding()>0) ? -1 : +1; // if we have a quark, we need antiquark (or diquark)
+    pDefPair QuarkPair = CreatePartonPair(IsParticle);
+    created = QuarkPair.second;
+    return hadronizer->Build(QuarkPair.first, decay);
+void G4VLongitudinalStringDecay::SetStrangenessSuppression(G4double aValue)
+{
+        if ( PastInitPhase ) {
+                throw G4HadronicException(__FILE__, __LINE__, "4VLongitudinalStringDecay::SetStrangenessSuppression after FragmentString() not allowed");
+        } else {
+                StrangeSuppress = aValue;
+        }
+}
 //----------------------------------------------------------------------------------------------------------
-G4ParticleDefinition *G4VLongitudinalStringDecay::DiQuarkSplitup(
-                                                        G4ParticleDefinition* decay,
-                                                        G4ParticleDefinition *&created)
+{
-   //... can Diquark break or not?
-   if (G4UniformRand() < DiquarkBreakProb ){
-   //... Diquark break
-      G4int stableQuarkEncoding = decay->GetPDGEncoding()/1000;
-      G4int decayQuarkEncoding = (decay->GetPDGEncoding()/100)%10;
-      if (G4UniformRand() < 0.5)
+         {
-         G4int Swap = stableQuarkEncoding;
-         stableQuarkEncoding = decayQuarkEncoding;
-         decayQuarkEncoding = Swap;
+         }
-      G4int IsParticle=(decayQuarkEncoding>0) ? -1 : +1;
-                        // if we have a quark, we need antiquark)
-      pDefPair QuarkPair = CreatePartonPair(IsParticle,false);  // no diquarks wanted
-      //... Build new Diquark
-      G4int QuarkEncoding=QuarkPair.second->GetPDGEncoding();
-      G4int i10  = std::max(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
-      G4int i20  = std::min(std::abs(QuarkEncoding), std::abs(stableQuarkEncoding));
-      G4int spin = (i10 != i20 && G4UniformRand() <= 0.5)? 1 : 3;
-      G4int NewDecayEncoding = -1*IsParticle*(i10 * 1000 + i20 * 100 + spin);
-      created = FindParticle(NewDecayEncoding);
-      G4ParticleDefinition * decayQuark=FindParticle(decayQuarkEncoding);
-      G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decayQuark);
-      return had;
-//      return hadronizer->Build(QuarkPair.first, decayQuark);
-   } else {
-   //... Diquark does not break
-      G4int IsParticle=(decay->GetPDGEncoding()>0) ? +1 : -1;
-                        // if we have a diquark, we need quark)
-      pDefPair QuarkPair = CreatePartonPair(IsParticle,false);  // no diquarks wanted
-      created = QuarkPair.second;
-      G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decay);
-      return had;
-//      return G4ParticleDefinition * had=hadronizer->Build(QuarkPair.first, decay);
+   }
+}
-//-----------------------------------------------------------------------------------------
-G4KineticTrack * G4VLongitudinalStringDecay::Splitup(
-                        G4FragmentingString *string,
-                        G4FragmentingString *&newString)
+{
-       //... random choice of string end to use for creating the hadron (decay)
-       SideOfDecay = (G4UniformRand() < 0.5)? 1: -1;
-       if (SideOfDecay < 0)
+       {
-          string->SetLeftPartonStable();
-       } else
+       {
-          string->SetRightPartonStable();
+       }
-       G4ParticleDefinition *newStringEnd;
-       G4ParticleDefinition * HadronDefinition;
-       if (string->DecayIsQuark())
+       {
-           HadronDefinition= QuarkSplitup(string->GetDecayParton(), newStringEnd);
-       } else {
-           HadronDefinition= DiQuarkSplitup(string->GetDecayParton(), newStringEnd);
+       }
-// create new String from old, ie. keep Left and Right order, but replace decay
-       G4LorentzVector* HadronMomentum=SplitEandP(HadronDefinition, string);
-       G4KineticTrack * Hadron =0;
-       if ( HadronMomentum != 0 ) {
-           G4ThreeVector   Pos;
-           Hadron = new G4KineticTrack(HadronDefinition, 0,Pos, *HadronMomentum);
-           newString=new G4FragmentingString(*string,newStringEnd,
-                                        HadronMomentum);
-           delete HadronMomentum;
+       }
-       return Hadron;
+}
-//----------------------------------------------------------------------------------------------------------
-G4ExcitedString *G4VLongitudinalStringDecay::CPExcited(const G4ExcitedString & in)
+{
-        G4Parton *Left=new G4Parton(*in.GetLeftParton());
-        G4Parton *Right=new G4Parton(*in.GetRightParton());
-        return new G4ExcitedString(Left,Right,in.GetDirection());
+}
-G4double G4VLongitudinalStringDecay::FragmentationMass(
-                const G4FragmentingString *
-                const string,
-                Pcreate build,
-                pDefPair * pdefs)
+{
-        G4double mass;
-        static G4bool NeedInit(true);
-        static std::vector<double> nomix;
-        static G4HadronBuilder * minMassHadronizer;
-        if ( NeedInit )
+        {
-           NeedInit = false;
-           nomix.resize(6);
-           for ( G4int i=0; i<6 ; i++ ) nomix[i]=0;
-//         minMassHadronizer=new G4HadronBuilder(pspin_meson,pspin_barion,nomix,nomix);
-           minMassHadronizer=hadronizer;
+        }
-        if ( build==0 ) build=&G4HadronBuilder::BuildLowSpin;
-        G4ParticleDefinition *Hadron1, *Hadron2=0;
-        if (!string->FourQuarkString() )
+        {
-           // spin 0 meson or spin 1/2 barion will be built
-           Hadron1 = (minMassHadronizer->*build)(string->GetLeftParton(),
-                                      string->GetRightParton());
-           mass= (Hadron1)->GetPDGMass();
-        } else
+        {
-           //... string is qq--qqbar: Build two stable hadrons,
-           //...    with extra uubar or ddbar quark pair
-           G4int iflc = (G4UniformRand() < 0.5)? 1 : 2;
-           if (string->GetLeftParton()->GetPDGEncoding() < 0) iflc = -iflc;
-           //... theSpin = 4; spin 3/2 baryons will be built
-           Hadron1 = (minMassHadronizer->*build)(string->GetLeftParton(),FindParticle(iflc));
-           Hadron2 =(minMassHadronizer->*build)(string->GetRightParton(),FindParticle(-iflc));
-           mass = (Hadron1)->GetPDGMass() + (Hadron2)->GetPDGMass();
+        }
-        if ( pdefs != 0 )
-        { // need to return hadrons as well....
-           pdefs->first  = Hadron1;
-           pdefs->second = Hadron2;
+        }
-        return mass;
+}
-G4KineticTrackVector* G4VLongitudinalStringDecay::LightFragmentationTest(const
-                G4ExcitedString * const string)
+{
-   // Check string decay threshold
-        G4KineticTrackVector * result=0;  // return 0 when string exceeds the mass cut
-        pDefPair hadrons((G4ParticleDefinition *)0,(G4ParticleDefinition *)0);
-        G4FragmentingString aString(*string);
-        if ( sqr(FragmentationMass(&aString,0,&hadrons)+MassCut) < aString.Mass2()) {
-                return 0;
+        }
-        result=new G4KineticTrackVector;
-        if ( hadrons.second ==0 )
+        {
-                 // Substitute string by light hadron, Note that Energy is not conserved here!
-#ifdef DEBUG_LightFragmentationTest
-               G4cout << "VlongSF Warning replacing string by single hadron "
-                      << hadrons.first->GetParticleName()
-                      << "string .. " << string->Get4Momentum() << " "
-                      << string->Get4Momentum().m() << G4endl;
-#endif
-               G4ThreeVector Mom3 = string->Get4Momentum().vect();
-               G4LorentzVector Mom(Mom3,
-                                   std::sqrt(Mom3.mag2() + sqr(hadrons.first->GetPDGMass())));
-               result->push_back(new G4KineticTrack(hadrons.first, 0, string->GetPosition(), Mom));
-        } else
+        {
-           //... string was qq--qqbar type: Build two stable hadrons,
-#ifdef DEBUG_LightFragmentationTest
-               G4cout << "VlongSF Warning replacing qq-qqbar string by TWO hadrons "
-                      << hadrons.first->GetParticleName() << " / "
-                      << hadrons.second->GetParticleName()
-                      << "string .. " << string->Get4Momentum() << " "
-                      << string->Get4Momentum().m() << G4endl;
-#endif
-               G4LorentzVector  Mom1, Mom2;
-               Sample4Momentum(&Mom1, hadrons.first->GetPDGMass(),
-                               &Mom2,hadrons.second->GetPDGMass(),
-                                string->Get4Momentum().mag());
-               result->push_back(new G4KineticTrack(hadrons.first, 0, string->GetPosition(), Mom1));
-               result->push_back(new G4KineticTrack(hadrons.second, 0, string->GetPosition(), Mom2));
-               G4ThreeVector Velocity = string->Get4Momentum().boostVector();
-               result->Boost(Velocity);
+        }
-        return result;
+}
-//----------------------------------------------------------------------------------------------------------
-G4ParticleDefinition* G4VLongitudinalStringDecay::FindParticle(G4int Encoding)
+   {
-   G4ParticleDefinition* ptr = G4ParticleTable::GetParticleTable()->FindParticle(Encoding);
-      if (ptr == NULL)
+       {
-       G4cout << "Particle with encoding "<<Encoding<<" does not exist!!!"<<G4endl;
-       throw G4HadronicException(__FILE__, __LINE__, "Check your particle table");
+       }
-   return ptr;
+   }
-//----------------------------------------------------------------------------------------------------------
-void G4VLongitudinalStringDecay::SetSigmaTransverseMomentum(G4double aValue)
+{
-        if ( PastInitPhase ) {
-                throw G4HadronicException(__FILE__, __LINE__, "4VLongitudinalStringDecay::SetSigmaTransverseMomentum after FragmentString() not allowed");
-        } else {
-                SigmaQT = aValue;
+        }
+}
-//----------------------------------------------------------------------------------------------------------
-void G4VLongitudinalStringDecay::SetStrangenessSuppression(G4double aValue)
+{
-        if ( PastInitPhase ) {
-                throw G4HadronicException(__FILE__, __LINE__, "4VLongitudinalStringDecay::SetStrangenessSuppression after FragmentString() not allowed");
-        } else {
-                StrangeSuppress = aValue;
+        }
+}
-//----------------------------------------------------------------------------------------------------------
 void G4VLongitudinalStringDecay::SetDiquarkSuppression(G4double aValue)
+{
 …
+        }
+}
+//----------------------------------------------------------------------------------------
 void G4VLongitudinalStringDecay::SetDiquarkBreakProbability(G4double aValue)
 …
+        }
+}
+//-------------------------------------------------------------------------------------------
+void G4VLongitudinalStringDecay::SetStringTensionParameter(G4double aValue)// Uzhi 20 June 08
+{
+          Kappa = aValue * GeV/fermi;
+}
+//*******************************************************************************************************
+//**************************************************************************************

trunk/source/processes/hadronic/models/parton_string/management/History

-                      r819
+                      r962
 $Id: History,v 1.3 2007/04/24 10:37:10 gunter Exp $
+$Id: History,v 1.4 2008/04/01 08:20:25 vuzhinsk Exp $
 -------------------------------------------------------------------
 …
      model
+-March-2008 V. Uzhinsky               Tag : had-partonstring-mgt-V09-01-00
+  - G4FTFCrossSection.cc and G4FTFCrossSection.hh were re-named into
+    G4FTFParameters.cc and .hh, and moved to diffraction directory.
+    The corresponding class was re-named too. All of these characterize
+    the content of the files more exactly.

trunk/source/processes/hadronic/models/parton_string/management/include/G4EventGenerator.hh

r819	r962
26	26	//
27	27	// $Id: G4EventGenerator.hh,v 1.3 2006/06/29 20:55:13 gunter Exp $
28		// GEANT4 tag $Name: $
	28	// GEANT4 tag $Name: geant4-09-02-ref-02 $
29	29	//
30	30	#ifndef G4EventGenerator_h

trunk/source/processes/hadronic/models/parton_string/management/include/G4InteractionCode.hh

r819	r962
26	26	//
27	27	// $Id: G4InteractionCode.hh,v 1.3 2006/06/29 20:55:15 gunter Exp $
28		// GEANT4 tag $Name: $
	28	// GEANT4 tag $Name: geant4-09-02-ref-02 $
29	29	//
30	30	#ifndef G4InteractionCode_h

trunk/source/processes/hadronic/models/parton_string/management/include/G4InteractionContent.hh

r819	r962
26	26	//
27	27	// $Id: G4InteractionContent.hh,v 1.4 2007/01/24 10:28:54 gunter Exp $
28		// GEANT4 tag $Name: $
	28	// GEANT4 tag $Name: geant4-09-02-ref-02 $
29	29	//
30	30

trunk/source/processes/hadronic/models/parton_string/management/include/G4PomeronCrossSection.hh

r819	r962
28	28	//
29	29	// $Id: G4PomeronCrossSection.hh,v 1.3 2006/06/29 20:55:19 gunter Exp $
30		// GEANT4 tag $Name: $
	30	// GEANT4 tag $Name: geant4-09-02-ref-02 $
31	31	//
32	32	#include "G4Proton.hh"

trunk/source/processes/hadronic/models/parton_string/management/include/G4StringModel.hh

r819	r962
26	26	//
27	27	// $Id: G4StringModel.hh,v 1.3 2006/06/29 20:55:23 gunter Exp $
28		// GEANT4 tag $Name: $
	28	// GEANT4 tag $Name: geant4-09-02-ref-02 $
29	29	//
30	30	#ifndef G4StringModel_h

trunk/source/processes/hadronic/models/parton_string/management/include/G4VParticipants.hh

-                      r819
+                      r962
 //
 //
 // $Id: G4VParticipants.hh,v 1.3 2006/06/29 20:55:25 gunter Exp $
 // GEANT4 tag $Name:  $
+// $Id: G4VParticipants.hh,v 1.4 2008/05/19 13:03:20 vuzhinsk Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 …
         if ( theNucleus == NULL ) theNucleus = new G4Fancy3DNucleus();
         theNucleus->Init(theA, theZ);
+        theNucleus->SortNucleonsInZ();    // Uzhi 16.05.08 Sorting of nucleon-Z
+}

trunk/source/processes/hadronic/models/parton_string/management/include/G4VPartonStringModel.hh

r819	r962
26	26	//
27	27	// $Id: G4VPartonStringModel.hh,v 1.3 2006/06/29 20:55:27 gunter Exp $
28		// GEANT4 tag $Name: $
	28	// GEANT4 tag $Name: geant4-09-02-ref-02 $
29	29	//
30	30	#ifndef G4VPartonStringModel_h

trunk/source/processes/hadronic/models/parton_string/management/include/G4VSplitableHadron.hh

-                      r819
+                      r962
 //
 //
 // $Id: G4VSplitableHadron.hh,v 1.3 2006/06/29 20:55:29 gunter Exp $
 // GEANT4 tag $Name:  $
+// $Id: G4VSplitableHadron.hh,v 1.4 2008/05/19 13:03:20 vuzhinsk Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 …
       void SetCollisionCount(G4int aCount);
+      void SetTimeOfCreation(G4double aTime);           // Uzhi 7.05.08
+      G4double GetTimeOfCreation();                     // Uzhi 7.05.08
       void SetPosition(const G4ThreeVector &aPosition);
       const G4ThreeVector & GetPosition() const;
 …
       G4bool IsSplit() { return isSplit;}
+      G4int GetSoftCollisionCount();
+  protected:
-  protected:
-      G4int GetSoftCollisionCount();
       void Splitting() {isSplit = true;}
 …
       G4LorentzVector the4Momentum;
+      G4double TimeOfCreation;    // Uzhi 7.05.08
       G4ThreeVector thePosition;
       G4int theCollisionCount;
 …
+}
+inline void G4VSplitableHadron::SetTimeOfCreation(G4double aTime)  // Uzhi 7.05.08
+{
+        TimeOfCreation=aTime;
+}
+inline G4double G4VSplitableHadron::GetTimeOfCreation()           // Uzhi 7.05.08
+{
+        return TimeOfCreation;
+}
 inline void G4VSplitableHadron::SetPosition(const G4ThreeVector &aPosition)

trunk/source/processes/hadronic/models/parton_string/management/include/G4VStringFragmentation.hh

r819	r962
26	26	//
27	27	// $Id: G4VStringFragmentation.hh,v 1.3 2006/06/29 20:55:31 gunter Exp $
28		// GEANT4 tag $Name: $
	28	// GEANT4 tag $Name: geant4-09-02-ref-02 $
29	29	//
30	30	#ifndef G4VStringFragmentation_h

trunk/source/processes/hadronic/models/parton_string/management/include/G4VertexCode.hh

r819	r962
26	26	//
27	27	// $Id: G4VertexCode.hh,v 1.3 2006/06/29 20:55:33 gunter Exp $
28		// GEANT4 tag $Name: $
	28	// GEANT4 tag $Name: geant4-09-02-ref-02 $
29	29	//
30	30	#ifndef G4VertexCode_h

trunk/source/processes/hadronic/models/parton_string/management/src/G4EventGenerator.cc

r819	r962
26	26	//
27	27	// $Id: G4EventGenerator.cc,v 1.4 2006/06/29 20:55:37 gunter Exp $
28		// GEANT4 tag $Name: $
	28	// GEANT4 tag $Name: geant4-09-02-ref-02 $
29	29	//
30	30	// G4EventGenerator

trunk/source/processes/hadronic/models/parton_string/management/src/G4InteractionContent.cc

r819	r962
26	26	//
27	27	// $Id: G4InteractionContent.cc,v 1.4 2006/06/29 20:55:39 gunter Exp $
28		// GEANT4 tag $Name: $
	28	// GEANT4 tag $Name: geant4-09-02-ref-02 $
29	29	//
30	30	// ------------------------------------------------------------

trunk/source/processes/hadronic/models/parton_string/management/src/G4PomeronCrossSection.cc

r819	r962
26	26	//
27	27	// $Id: G4PomeronCrossSection.cc,v 1.6 2006/11/07 12:51:39 gunter Exp $
28		// GEANT4 tag $Name: $
	28	// GEANT4 tag $Name: geant4-09-02-ref-02 $
29	29	//
30	30

trunk/source/processes/hadronic/models/parton_string/management/src/G4StringModel.cc

r819	r962
26	26	//
27	27	// $Id: G4StringModel.cc,v 1.4 2006/06/29 20:55:45 gunter Exp $
28		// GEANT4 tag $Name: $
	28	// GEANT4 tag $Name: geant4-09-02-ref-02 $
29	29	//
30	30	// G4StringModel

trunk/source/processes/hadronic/models/parton_string/management/src/G4VParticipants.cc

r819	r962
26	26	//
27	27	// $Id: G4VParticipants.cc,v 1.3 2006/06/29 20:55:47 gunter Exp $
28		// GEANT4 tag $Name: $
	28	// GEANT4 tag $Name: geant4-09-02-ref-02 $
29	29	//
30	30	// ------------------------------------------------------------

trunk/source/processes/hadronic/models/parton_string/management/src/G4VPartonStringModel.cc

r819	r962
26	26	//
27	27	// $Id: G4VPartonStringModel.cc,v 1.5 2007/01/24 10:29:30 gunter Exp $
28		// GEANT4 tag $Name: $
	28	// GEANT4 tag $Name: geant4-09-02-ref-02 $
29	29	//
30	30	//// ------------------------------------------------------------

trunk/source/processes/hadronic/models/parton_string/management/src/G4VSplitableHadron.cc

-                      r819
+                      r962
 //
 //
 // $Id: G4VSplitableHadron.cc,v 1.4 2006/06/29 20:55:51 gunter Exp $
 // GEANT4 tag $Name:  $
+// $Id: G4VSplitableHadron.cc,v 1.5 2008/05/19 13:03:20 vuzhinsk Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 …
 G4VSplitableHadron::G4VSplitableHadron()
       :  theDefinition(NULL), theCollisionCount(0), isSplit(false)
+      :  theDefinition(NULL), TimeOfCreation(0.), theCollisionCount(0), isSplit(false) // Uzhi 8.05.08
+{
+}
 G4VSplitableHadron::G4VSplitableHadron(const G4ReactionProduct & aPrimary)
       :  theCollisionCount(0), isSplit(false)
+      :  TimeOfCreation(0.), theCollisionCount(0), isSplit(false)                     // Uzhi 8.05.08
+{
         theDefinition=aPrimary.GetDefinition();
 …
 G4VSplitableHadron::G4VSplitableHadron(const G4Nucleon & aNucleon)
+{
+        theCollisionCount=0;
+  isSplit = false;
+        theDefinition=aNucleon.GetParticleType();
+        the4Momentum=aNucleon.GetMomentum();
+        thePosition=aNucleon.GetPosition();
+        TimeOfCreation   = 0.;   // Uzhi 8.05.08
+        theCollisionCount= 0;
+        isSplit          = false;
+        theDefinition    =aNucleon.GetParticleType();
+        the4Momentum     =aNucleon.GetMomentum();
+        thePosition      =aNucleon.GetPosition();
+}
 G4VSplitableHadron::G4VSplitableHadron(const G4VKineticNucleon * aNucleon)
+{
+        theCollisionCount=0;
+  isSplit = false;
+        theDefinition=aNucleon->GetDefinition();
+        the4Momentum=aNucleon->Get4Momentum();
+        thePosition=aNucleon->GetPosition();
+        TimeOfCreation   = 0.;   // Uzhi 8.05.08
+        theCollisionCount= 0;
+        isSplit          = false;
+        theDefinition    =aNucleon->GetDefinition();
+        the4Momentum     =aNucleon->Get4Momentum();
+        thePosition      =aNucleon->GetPosition();
+}
 G4VSplitableHadron::G4VSplitableHadron(const G4VSplitableHadron &right)
+{
+        theCollisionCount=0;
+  isSplit = false;
+        theDefinition= right.GetDefinition();
+        the4Momentum= right.Get4Momentum();
+        thePosition=  right.GetPosition();
+        TimeOfCreation   = 0.;   // Uzhi 8.05.08
+        theCollisionCount= 0;
+        isSplit          = false;
+        theDefinition    = right.GetDefinition();
+        the4Momentum     = right.Get4Momentum();
+        thePosition      =  right.GetPosition();
+}

trunk/source/processes/hadronic/models/parton_string/management/src/G4VStringFragmentation.cc

r819	r962
26	26	//
27	27	// $Id: G4VStringFragmentation.cc,v 1.4 2006/06/29 20:55:53 gunter Exp $
28		// GEANT4 tag $Name: $
	28	// GEANT4 tag $Name: geant4-09-02-ref-02 $
29	29	//
30	30	// G4VStringFragmentation

trunk/source/processes/hadronic/models/parton_string/qgsm/History

-                      r819
+                      r962
 $Id: History,v 1.4.2.1 2008/04/23 15:25:01 gcosmo Exp $
+$Id: History,v 1.6 2008/09/19 09:54:23 gunter Exp $
 -------------------------------------------------------------------
 …
      ---------------------------------------------------------------
+Apr 2008 Gabriele Cosmo (hadr-qgsm-V09-00-00)
+Sep 2008 G.Folger       (hadr-qgsm-V09-01-01)
 ------------------------------------------------
+- Tag for release 9.1.p02.
+- Fix for bug found on windows in G4QGSParticipants.cc, bug 1018:
+   decrement of iterator fails, improve logic to not decrement.
 Mar 2008 Dennis Wright (hadr-qgsm-V09-01-00)
 -----------------------------------------------
 - Fix gcc-4.3 compiler warnings at lines 293, 395 of G4QGSMSplittableHadron.cc
+-  fix gcc-4.3 compiler warnings at lines 293, 395 of G4QGSMSplittableHadron.cc
 Apr 2007 Gunter Folger  (hadr-qgsm-V08-02-02)
 ------------------------------------------------
 - Merge in change done by ftf dev; ie. in G4QGSParticipants, theDiffExcitaton
   is constructed with default arguments.
+-  merge in change done by ftf dev; ie. in G4QGSParticipants, theDiffExcitaton
+   is constructed with default arguments.
 Jan 2007 Gunter Folger  (hadr-qgsm-V08-02-01)
 …
 Jan 2007 Gunter Folger  (hadr-qgsm-V08-02-00)
 ------------------------------------------------
 - Correct E-p non-conservation in QGS. In 4QGSMSplitableHadron.cc the smaller
   of the lightcone momenta Q+/Q- was ignored.
 - G4QGSMSplitableHadron correct divide by 0 in SampleX()
 - Add debugging output to several classes
+-  Correct E-p non-conservation in QGS. In 4QGSMSplitableHadron.cc the smaller
+   of the lightcone momenta Q+/Q- was ignored.
+-  G4QGSMSplitableHadron correct divide by 0 in SampleX()
+-  Add debugging output to several classes
 Nov 2005 Gabriele Cosmo (hadr-qgsm-V07-01-00)

trunk/source/processes/hadronic/models/parton_string/qgsm/src/G4QGSParticipants.cc

-                      r819
+                      r962
   std::vector<G4InteractionContent*>::iterator i;
   G4LorentzVector str4Mom;
+  for(i = theInteractions.begin(); i != theInteractions.end(); i++)
+  i = theInteractions.begin();
+  while ( i != theInteractions.end() )
+  {
     G4InteractionContent* anIniteraction = *i;
 …
+      }
       delete *i;
+      i=theInteractions.erase(i);
+      i--;
+    }
+      i=theInteractions.erase(i);    // i now points to the next interaction
+    } else {
+      i++;
+    }
+  }
 #ifdef debug_G4QGSPart_PSoftColl

Context Navigation

Legend:

Download in other formats: