Changeset 1228 for trunk/source/processes/hadronic/models/parton_string/diffraction/src

Timestamp:

Jan 8, 2010, 11:56:51 AM (14 years ago)

Author:

garnier

Message:

update geant4.9.3 tag

Location:

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

Files:

• G4DiffractiveExcitation.cc (modified) (56 diffs)
• G4DiffractiveSplitableHadron.cc (modified) (1 diff)
• G4ElasticHNScattering.cc (modified) (7 diffs)
• G4FTFModel.cc (modified) (34 diffs)
• G4FTFParameters.cc (modified) (4 diffs)
• G4FTFParticipants.cc (modified) (5 diffs)

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

@@ -25 +25 @@
 //
 //
-// $Id: G4DiffractiveExcitation.cc,v 1.16 2009/10/06 10:10:36 vuzhinsk Exp $
+// $Id: G4DiffractiveExcitation.cc,v 1.21 2009/12/15 19:14:31 vuzhinsk Exp $
 // ------------------------------------------------------------
 //      GEANT 4 class implemetation file
@@ -73 +73 @@
                      G4VSplitableHadron *target,
                      G4FTFParameters    *theParameters,
-                     G4ElasticHNScattering *theElastic) const  // Uzhi 03.09.09
+                     G4ElasticHNScattering *theElastic) const  
 {
 // -------------------- Projectile parameters -----------------------
      G4LorentzVector Pprojectile=projectile->Get4Momentum();
-//G4cout<<"Excite P "<<Pprojectile<<G4endl;
 
      if(Pprojectile.z() < 0.)
@@ -88 +87 @@
 
      G4int    ProjectilePDGcode=projectile->GetDefinition()->GetPDGEncoding();
-//     G4ParticleDefinition * ProjectileDefinition = projectile->GetDefinition();
      G4int    absProjectilePDGcode=std::abs(ProjectilePDGcode);
 
@@ -115 +113 @@
      G4double M0target = Ptarget.mag();
 
-//G4cout<<"Excite T "<<Ptarget<<G4endl;
-//G4cout<<"PDGs "<<ProjectilePDGcode<<" "<<TargetPDGcode<<G4endl;
-//G4cout<<"Masses "<<M0projectile<<" "<<M0target<<G4endl;
-
      G4double TargetRapidity = Ptarget.rapidity();
 
@@ -149 +143 @@
 
      G4LorentzVector Ptmp=toCms*Pprojectile;
+
      if ( Ptmp.pz() <= 0. )
      {
@@ -167 +162 @@
 
      G4double SqrtS=std::sqrt(S);
-//G4cout<<" SqrtS "<<SqrtS<<" 2200 "<<G4endl;
-     if(absProjectilePDGcode > 1000 && SqrtS < 2300*MeV && SqrtS < SumMasses)
+               
+     if(absProjectilePDGcode > 1000 && (SqrtS < 2300*MeV || SqrtS < SumMasses))
      {target->SetStatus(2);  return false;}  // The model cannot work for
                                              // p+p-interactions
@@ -174 +169 @@
 
      if(( absProjectilePDGcode == 211 || ProjectilePDGcode ==  111) && 
-        (SqrtS < 1600*MeV) && (SqrtS < SumMasses))
+        ((SqrtS < 1600*MeV) || (SqrtS < SumMasses)))
      {target->SetStatus(2);  return false;}  // The model cannot work for
                                              // Pi+p-interactions
@@ -181 +176 @@
      if(( (absProjectilePDGcode == 321) || (ProjectilePDGcode == -311)   ||
           (absProjectilePDGcode == 311) || (absProjectilePDGcode == 130) ||
-          (absProjectilePDGcode == 310)) && (SqrtS < 1600*MeV) && (SqrtS < SumMasses)) 
+          (absProjectilePDGcode == 310)) && ((SqrtS < 1600*MeV) || (SqrtS < SumMasses))) 
      {target->SetStatus(2);  return false;}  // The model cannot work for
                                              // K+p-interactions
@@ -223 +218 @@
 // Charge exchange can be possible for baryons -----------------
 
-//G4cout<<1./std::exp(-2.0*(ProjectileRapidity - TargetRapidity))<<G4endl;
-//G4int Uzhi; G4cin>>Uzhi;
 // Getting the values needed for exchange ----------------------
      G4double MagQuarkExchange        =theParameters->GetMagQuarkExchange();
@@ -231 +224 @@
 
 //     G4double NucleonMass=
-              (G4ParticleTable::GetParticleTable()->FindParticle(2112))->GetPDGMass();     
+//              (G4ParticleTable::GetParticleTable()->FindParticle(2112))->GetPDGMass();     
      G4double DeltaMass=
               (G4ParticleTable::GetParticleTable()->FindParticle(2224))->GetPDGMass();
@@ -239 +232 @@
         std::exp(-SlopeQuarkExchange*(ProjectileRapidity - TargetRapidity)))
      {    
-
-//G4cout<<"Exchange "<<G4endl;
-
       G4int NewProjCode(0), NewTargCode(0);
 
-//G4bool StandardExcitation(false); // ================================= 
       G4int ProjQ1(0), ProjQ2(0), ProjQ3(0);
 
@@ -264 +253 @@
       G4int TargExchangeQ(0);
 
-//G4cout<<"Targ "<<TargQ1<<" "<<TargQ2<<" "<<TargQ3<<G4endl;
       if(absProjectilePDGcode < 1000 )
       {    // projectile is meson ----------------- 
-//G4cout<<"Proj Q1 Q2 "<<ProjQ1<<" "<<ProjQ2<<G4endl;
-
        if(ProjQ1 > 0 ) // ProjQ1 is quark
        {  
@@ -283 +269 @@
          TargExchangeQ = TargQ3;  TargQ3=ProjExchangeQ; ProjQ1=TargExchangeQ;
         }
-//G4cout<<" Pr Tr "<<ProjQ1<<" "<<TargQ1<<" "<<TargQ2<<" "<<TargQ3<<G4endl;
        } else          // ProjQ2 is quark
        {  
@@ -310 +295 @@
        }
 
-//       projectile->SetDefinition(
-//       G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode)); 
-
        NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3);
- 
+
        if( (TargQ1 == TargQ2) && (TargQ1 == TargQ3) &&
            (SqrtS > M0projectile+DeltaMass))           {NewTargCode +=2;} //Create Delta isobar
@@ -326 +308 @@
                (SqrtS > M0projectile+DeltaMass))      {NewTargCode +=2;}  //Create Delta isobar
        else                                           {}                 //Save initial nucleon
-
-       target->SetDefinition(
-       G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode)); 
-//G4cout<<"New PDGs "<<NewProjCode<<" "<<NewTargCode<<G4endl;
-//G4int Uzhi; G4cin>>Uzhi;
-       //}
+//       target->SetDefinition(                                          // Fix 15.12.09
+//       G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode));// Fix 15.12.09 
       } else 
       {    // projectile is baryon ----------------
-G4double Same=0.; //0.3; //0.5;
+       G4double Same=0.; //0.3; //0.5;
        G4bool ProjDeltaHasCreated(false);
        G4bool TargDeltaHasCreated(false);
@@ -341 +319 @@
        if(G4UniformRand() < 0.5)     // Sampling exchange quark from proj. or targ.
        {   // Sampling exchanged quark from the projectile ---
-//G4cout<<"Proj Exc"<<G4endl;
         if( Ksi < 0.333333 ) 
         {ProjExchangeQ = ProjQ1;}
@@ -348 +325 @@
         else
         {ProjExchangeQ = ProjQ3;}
-//G4cout<<"ProjExchangeQ "<<ProjExchangeQ<<G4endl;
-
-        if((ProjExchangeQ != TargQ1)||(G4UniformRand()<Same)) // Vova
+
+        if((ProjExchangeQ != TargQ1)||(G4UniformRand()<Same)) 
         {
          TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjExchangeQ=TargExchangeQ;
         } else
-        if((ProjExchangeQ != TargQ2)||(G4UniformRand()<Same)) // Vova
+        if((ProjExchangeQ != TargQ2)||(G4UniformRand()<Same)) 
         {
          TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjExchangeQ=TargExchangeQ;
@@ -361 +337 @@
          TargExchangeQ = TargQ3;  TargQ3=ProjExchangeQ; ProjExchangeQ=TargExchangeQ;
         }
-//G4cout<<"ProjExchangeQ "<<ProjExchangeQ<<G4endl;
+
         if( Ksi < 0.333333 ) 
         {ProjQ1=ProjExchangeQ;}
@@ -369 +345 @@
         {ProjQ3=ProjExchangeQ;}
 
-/*
-        NewProjCode = NewNucleonId(ProjQ1, ProjQ2, ProjQ3); // *****************************
-
-        if((ProjQ1==ProjQ2) && (ProjQ1==ProjQ3)) {NewProjCode +=2; ProjDeltaHasCreated=true;}
-        else if(projectile->GetDefinition()->GetPDGiIsospin() == 3)// Projectile was Delta
-        { if(G4UniformRand() > DeltaProbAtQuarkExchange)
-                                                 {NewProjCode +=2; ProjDeltaHasCreated=true;} 
-          else                                   {NewProjCode +=0; ProjDeltaHasCreated=false;}
-        }
-        else                                                       // Projectile was Nucleon
-        {
-         if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS > DeltaMass+M0target)) 
-                                                 {NewProjCode +=2; ProjDeltaHasCreated=true;}
-         else                                    {NewProjCode +=0; ProjDeltaHasCreated=false;}
-        } 
-
-        NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3); // *****************************
-
-        if((TargQ1==TargQ2) && (TargQ1==TargQ3)) {NewTargCode +=2; TargDeltaHasCreated=true;}  
-        else if(target->GetDefinition()->GetPDGiIsospin() == 3)    // Target was Delta
-        { if(G4UniformRand() > DeltaProbAtQuarkExchange)
-                                                 {NewTargCode +=2; TargDeltaHasCreated=true;} 
-          else                                   {NewTargCode +=0; TargDeltaHasCreated=false;}
-        }
-        else                                                       // Target was Nucleon
-        {
-         if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS > M0projectile+DeltaMass)) 
-                                                 {NewTargCode +=2; TargDeltaHasCreated=true;}
-         else                                    {NewTargCode +=0; TargDeltaHasCreated=false;}
-        }         
-
-        if((absProjectilePDGcode == NewProjCode) && (absTargetPDGcode == NewTargCode))
-        { // Nothing was changed! It is not right!?
-        }
-*/
-/*-----------------------------------------------------------------------------------------
-        if(!ProjDeltaHasCreated && !TargDeltaHasCreated) // No Deltas were created then
-        { if( G4UniformRand() < 0.5) {NewProjCode +=2; ProjDeltaHasCreated = true;}
-          else                       {NewTargCode +=2; TargDeltaHasCreated = true;} 
-        }
-
-        if(!(ProjDeltaHasCreated && TargDeltaHasCreated))
-        {
-         if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS >2600.)) 
-         { 
-G4cout<<"2 Deltas "<<G4endl;
-          if(!ProjDeltaHasCreated)
-               {NewProjCode +=2; ProjDeltaHasCreated = true;}
-          else {NewTargCode +=2; TargDeltaHasCreated = true;} // For Delta isobar
-         }
-        }
-*/
        } else
        {   // Sampling exchanged quark from the target -------
-//G4cout<<"Targ Exc"<<G4endl;
         if( Ksi < 0.333333 ) 
         {TargExchangeQ = TargQ1;}
@@ -430 +353 @@
         else
         {TargExchangeQ = TargQ3;}
-//G4cout<<"TargExchangeQ "<<TargExchangeQ<<G4endl;
-        if((TargExchangeQ != ProjQ1)||(G4UniformRand()<Same)) // Vova
+
+        if((TargExchangeQ != ProjQ1)||(G4UniformRand()<Same)) 
         {
          ProjExchangeQ = ProjQ1; ProjQ1=TargExchangeQ; TargExchangeQ=ProjExchangeQ;
         } else
-        if((TargExchangeQ != ProjQ2)||(G4UniformRand()<Same)) // Vova
+        if((TargExchangeQ != ProjQ2)||(G4UniformRand()<Same)) 
         {
          ProjExchangeQ = ProjQ2; ProjQ2=TargExchangeQ; TargExchangeQ=ProjExchangeQ;
@@ -442 +365 @@
          ProjExchangeQ = ProjQ3;  ProjQ3=TargExchangeQ; TargExchangeQ=ProjExchangeQ;
         }
-//G4cout<<"TargExchangeQ "<<TargExchangeQ<<G4endl;
+
         if( Ksi < 0.333333 ) 
         {TargQ1=TargExchangeQ;}
@@ -449 +372 @@
         else
         {TargQ3=TargExchangeQ;}
-/*
-        NewProjCode = NewNucleonId(ProjQ1, ProjQ2, ProjQ3);
-        if( (ProjQ1 == ProjQ2) && (ProjQ1 == ProjQ3) )  
-        {NewProjCode +=2; ProjDeltaHasCreated = true;}
-
-        NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3);
-        if( (TargQ1 == TargQ2) && (TargQ1 == TargQ3) )  
-        {NewTargCode +=2; TargDeltaHasCreated = true;}
-
-        if(!ProjDeltaHasCreated && !TargDeltaHasCreated)
-        {NewTargCode +=2; TargDeltaHasCreated = true;}
-
-        if(!(ProjDeltaHasCreated && TargDeltaHasCreated))
-        {
-         if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS >2600.)) 
-         { 
-G4cout<<"2 Deltas "<<G4endl;
-          if(!ProjDeltaHasCreated)
-               {NewProjCode +=2; ProjDeltaHasCreated = true;}
-          else {NewTargCode +=2; TargDeltaHasCreated = true;}
-         }
-        }
-*/
+
        } // End of sampling baryon
 
@@ -509 +410 @@
 // Forming baryons --------------------------------------------------
 
-//       if( ProjDeltaHasCreated && TargDeltaHasCreated ) {StandardExcitation = true;}
-
-//G4cout<<"New PDG "<<NewProjCode<<" "<<NewTargCode<<G4endl;
-//G4int Uzhi; G4cin>>Uzhi;
       } // End of if projectile is baryon ---------------------------
 
@@ -518 +415 @@
 // If we assume that final state hadrons after the charge exchange will be
 // in the ground states, we have to put ----------------------------------
+
       if(M0projectile < 
          (G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode))->GetPDGMass())
@@ -537 +435 @@
              2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2)
              /4./S;
-//G4cout<<"New Masses "<<M0projectile<<" "<<M0target<<G4endl;
-//G4cout<<"PZcms2 "<<PZcms2<<G4endl;
 
       if(PZcms2 < 0) {return false;}  // It can be if energy is not sufficient for Delta
@@ -556 +452 @@
       Ptarget.setE(std::sqrt(M0target2+PZcms2));
 
-//G4cout<<"Proj "<<Pprojectile<<" "<<Pprojectile.mag()<<G4endl;
-//G4cout<<"Targ "<<Ptarget<<" "<<Ptarget.mag()<<G4endl;
-
-//      if(!StandardExcitation)
       {
        Pprojectile.transform(toLab);
@@ -571 +463 @@
 
        G4bool Result= theElastic->ElasticScattering (projectile,target,theParameters);
-//G4cout<<"1 Delta result "<<Result<<"**********"<<G4endl;
+
        return Result;
       } 
-/*
-else
-      {
-       if(M0projectile > ProjectileDiffStateMinMass) 
-       { ProjectileDiffStateMinMass +=200.*MeV; ProjectileNonDiffStateMinMass +=200.*MeV;}
-
-       if(M0target > TargetDiffStateMinMass) 
-       { TargetDiffStateMinMass +=200.*MeV; TargetNonDiffStateMinMass +=200.*MeV;}
-
-       ProbOfDiffraction         = 1.;
-       ProbProjectileDiffraction =0.5;
-       ProbTargetDiffraction     =0.5;
-G4cout<<"Exc DD "<<M0projectile<<" "<<M0target<<" --------------------"<<G4endl;
-//       After that standard FTF excitation
-      }
-*/
      }  // End of charge exchange part ------------------------------
 
 // ------------------------------------------------------------------
-//     G4double ProbOfDiffraction=ProbProjectileDiffraction +
-//                                ProbTargetDiffraction;
-//G4cout<<ProbOfDiffraction<<" "<<ProbProjectileDiffraction<<" "<<ProbTargetDiffraction<<G4endl;
-//G4int Uzhi; G4cin>>Uzhi;
      if(ProbOfDiffraction!=0.)
      {
@@ -606 +478 @@
      }
 
-//ProbOfDiffraction=1.; //0.5; // Vova 
-
      G4double ProjectileDiffStateMinMass2    = ProjectileDiffStateMinMass    *
                                                ProjectileDiffStateMinMass;
@@ -617 +487 @@
      G4double TargetNonDiffStateMinMass2     = TargetNonDiffStateMinMass     *
                                                TargetNonDiffStateMinMass;
-
 
      G4double Pt2;
@@ -637 +506 @@
         if(G4UniformRand() < ProbProjectileDiffraction)
         { //-------- projectile diffraction ---------------
-//G4cout<<"Proj difr"<<G4endl;
          do {
 //             Generate pt
@@ -659 +527 @@
                      2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
                     /4./S;
-//G4cout<<"PZcms2 < 0 false "<<PZcms2<<G4endl;
+
              if(PZcms2 < 0 ) 
              {
@@ -703 +571 @@
         { // -------------- Target diffraction ----------------
          do {
-//G4cout<<"Targ difr"<<G4endl;
 //             Generate pt
 //             if (whilecount++ >= 500 && (whilecount%100)==0)
@@ -770 +637 @@
         else  //----------- Non-diffraction process ------------
         {
-//G4cout<<"Non difr"<<G4endl;
          do {
 //             Generate pt
@@ -841 +707 @@
            Ptarget     -= Qmomentum;
 
-//G4cout<<"Masses "<<Pprojectile.mag()<<" "<<Ptarget.mag()<<G4endl;
-//G4int Uzhi; G4cin>>Uzhi;
-
 // Transform back and update SplitableHadron Participant.
            Pprojectile.transform(toLab);
@@ -865 +728 @@
 
 // ---------------------------------------------------------------------
-//G4ExcitedString * G4DiffractiveExcitation::
-//           String(G4VSplitableHadron * hadron, G4bool isProjectile) const
 void G4DiffractiveExcitation::CreateStrings(G4VSplitableHadron * hadron, 
                                             G4bool isProjectile,
@@ -886 +747 @@
           return;
         }
+
         G4LorentzVector Phadron=hadron->Get4Momentum();
-/*
-G4cout<<"Create strings had "<<hadron->GetDefinition()->GetParticleName()<<" "<<Phadron<<" "<<Phadron.mag()<<G4endl;
-G4cout<<"isProjectile "<<isProjectile<<G4endl;
-G4cout<<"start Q "<<start->GetDefinition()->GetPDGEncoding()<<G4endl;
-G4cout<<"end   Q "<<end->GetDefinition()->GetPDGEncoding()<<G4endl;
-*/
+
         G4LorentzVector Pstart(0.,0.,0.,0.);
         G4LorentzVector Pend(0.,0.,0.,0.);
@@ -914 +771 @@
         G4double W = hadron->Get4Momentum().mag();
         G4double W2=W*W;
-//G4cout<<"W Wmin "<<W<<" "<<Wmin<<G4endl;
+
         G4double Pt(0.), x1(0.), x2(0.), x3(0.);
         G4bool Kink=false;
@@ -922 +779 @@
           G4double Pt2kink=theParameters->GetPt2Kink();
           Pt = std::sqrt(Pt2kink*(std::pow(W2/16./Pt2kink+1.,G4UniformRand()) - 1.));
-// Pt=0.;
-//G4cout<<"Pt2kink Pt Pt2 "<<Pt2kink<<" "<<Pt<<" "<<Pt*Pt<<G4endl;
+
           if(Pt > 500.*MeV)
           {
              G4double Ymax = std::log(W/2./Pt + std::sqrt(W2/4./Pt/Pt - 1.));
              G4double Y=Ymax*(1.- 2.*G4UniformRand());
-//G4cout<<"Ymax Y "<<Ymax<<" "<<Y<<G4endl;
+
              x1=1.-Pt/W*std::exp( Y);
              x3=1.-Pt/W*std::exp(-Y);
              x2=2.-x1-x3;
-//G4cout<<"X1 X2 X3 "<<x1<<" "<<x2<<" "<<x3<<G4endl;
+
              G4double Mass_startQ = 650.*MeV;
              if(PDGcode_startQ <  3) Mass_startQ =  325.*MeV;
@@ -947 +803 @@
      
              G4double P2_2=sqr((2.-x1-x3)*W/2.);
-//G4cout<<"P2_1 P2_2 P2_3 "<<P2_1<<" "<<P2_2<<" "<<P2_3<<G4endl;
-             if((P2_1 < 0.) || (P2_3 <0.))
+
+             if((P2_1 <= 0.) || (P2_3 <= 0.))
              { Kink=false;}
              else
@@ -955 +811 @@
                G4double P_2=std::sqrt(P2_2);
                G4double P_3=std::sqrt(P2_3);
-//G4cout<<"P_1 P_2 P_3 "<<P_1<<" "<<P_2<<" "<<P_3<<G4endl;
+
                G4double CosT12=(P2_3-P2_1-P2_2)/(2.*P_1*P_2);
                G4double CosT13=(P2_2-P2_1-P2_3)/(2.*P_1*P_3);
-//G4cout<<"CosT12 CosT13 "<<CosT12<<" "<<CosT13<<" "<<std::acos(CosT12)*180./3.14159<<" "<<std::acos(CosT13)*180./3.14159<<G4endl;
+//             Pt=P_2*std::sqrt(1.-CosT12*CosT12);  // because system was rotated 11.12.09
+
                if((std::abs(CosT12) >1.) || (std::abs(CosT13) > 1.)) 
                { Kink=false;}
@@ -964 +821 @@
                { 
                  Kink=true; 
+                 Pt=P_2*std::sqrt(1.-CosT12*CosT12);  // because system was rotated 11.12.09
                  Pstart.setPx(-Pt); Pstart.setPy(0.); Pstart.setPz(P_3*CosT13); 
                  Pend.setPx(   0.); Pend.setPy(  0.); Pend.setPz(          P_1); 
@@ -981 +839 @@
                  PkinkQ2=Pkink - PkinkQ1;
 //------------------------- Minimizing Pt1^2+Pt3^2 ------------------------------
-/*
-G4cout<<"Pstart "<<Pstart<<G4endl;
-G4cout<<"Pkink  "<<Pkink <<G4endl;
-G4cout<<"Pkink1 "<<PkinkQ1<<G4endl;
-G4cout<<"Pkink2 "<<PkinkQ2<<G4endl;
-G4cout<<"Pend   "<<Pend  <<G4endl;
-*/
+
                  G4double Cos2Psi=(sqr(x1) -sqr(x3)+2.*sqr(x3*CosT13))/
                           std::sqrt(sqr(sqr(x1)-sqr(x3)) + sqr(2.*x1*x3*CosT13));
                  G4double Psi=std::acos(Cos2Psi);
-//G4cout<<"Psi "<<Psi<<" "<<Psi*180./twopi<<G4endl;
 
 G4LorentzRotation Rotate;
@@ -997 +848 @@
 else             {Rotate.rotateY(pi-Psi);}                   
 Rotate.rotateZ(twopi*G4UniformRand());
-
-//G4cout<<"Rotate"<<G4endl;
 
 Pstart*=Rotate;
@@ -1005 +854 @@
 PkinkQ2*=Rotate;
 Pend*=Rotate;
-/*
-G4cout<<"Pstart "<<Pstart<<G4endl;
-G4cout<<"Pkink1 "<<PkinkQ1 <<G4endl;
-G4cout<<"Pkink2 "<<PkinkQ2 <<G4endl;
-G4cout<<"Pend   "<<Pend  <<G4endl;
-*/
+
                }
              }      // end of if((P2_1 < 0.) || (P2_3 <0.))
@@ -1018 +862 @@
 //--------------------------------------------------------------------------------
 
-//G4cout<<"Kink "<<Kink<<G4endl;
-
         if(Kink)
         {                                        // Kink is possible
@@ -1028 +870 @@
           for(unsigned int Iq=0; Iq <3; Iq++)
           {
-//G4cout<<Iq<<" "<<QuarkProbabilitiesAtGluonSplitUp[Iq]<<G4endl;
+
 if(Ksi > QuarkProbabilitiesAtGluonSplitUp[Iq]) QuarkInGluon++;}
 
-//G4cout<<"Gquark "<<QuarkInGluon<<G4endl;
           G4Parton * Gquark = new G4Parton(QuarkInGluon);
           G4Parton * Ganti_quark = new G4Parton(-QuarkInGluon);
-/*
-G4cout<<Gquark->GetDefinition()->GetParticleName()<<" "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->GetDefinition()->GetPDGMass()<<G4endl;
-G4cout<<Ganti_quark->GetDefinition()->GetParticleName()<<" "<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->GetDefinition()->GetPDGMass()<<G4endl;
-*/
-
-//G4int Uzhi; G4cin>>Uzhi;
 
 //-------------------------------------------------------------------------------
-/*
-DefineMomentumInZ(G4double aLightConeMomentum, G4bool aDirection); 
-      void Set4Momentum(const G4LorentzVector & aMomentum);
-      G4int PDGencoding;
-      G4ParticleDefinition * theDefinition;
-      G4LorentzVector theMomentum;
-      G4ThreeVector   thePosition;
-      
-      G4int theColour;
-      G4double theIsoSpinZ;
-      G4double theSpinZ;
-      
-      G4double theX;
-*/
-//-------------------------------------------------------------------------------
-
-//G4cout<<"Phadron "<<Phadron<<" mass "<<Phadron.mag()<<G4endl;
           G4LorentzRotation toCMS(-1*Phadron.boostVector());
-//G4cout<<"To lab"<<G4endl;
+
           G4LorentzRotation toLab(toCMS.inverse());
 
@@ -1067 +885 @@
           PkinkQ2.transform(toLab);
           Pend.transform(toLab);    end->Set4Momentum(Pend);
-/*
-G4cout<<"Pstart "<<start->GetDefinition()->GetPDGEncoding()<<Pstart<<G4endl;
-G4cout<<"Pkink1 "<<PkinkQ1 <<G4endl;
-G4cout<<"Pkink2 "<<PkinkQ2 <<G4endl;
-G4cout<<"Pend   "<<end->GetDefinition()->GetPDGEncoding()<<Pend  <<G4endl;
-*/
-// !!!    G4ExcitedString * FirstString(0); G4ExcitedString * SecondString(0);
+
           G4int absPDGcode=std::abs(hadron->GetDefinition()->GetPDGEncoding());
-
-//G4cout<<"absPDGcode "<<absPDGcode<<G4endl;
 
           if(absPDGcode < 1000)
@@ -1088 +898 @@
                 Ganti_quark->Set4Momentum(PkinkQ1);
                 Gquark->Set4Momentum(PkinkQ2);
-/*
-G4cout<<" Proj Meson end Q"<<G4endl;
-G4cout<<"First string  ============ "<<G4endl;
-G4cout<<"end  Q "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl;
-G4cout<<"G antiQ"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl;
-G4cout<<"Sum P  "<<(Ganti_quark->Get4Momentum()+end->Get4Momentum())<<G4endl;
-G4cout<<"Secondstring   ============ "<<G4endl;
-G4cout<<"G    Q "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl;
-G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl;
-
-G4cout<<"Sum P  "<<(Gquark->Get4Momentum()+start->Get4Momentum())<<G4endl;
-*/
+
               } else
               {                            // Anti_Quark on the end
@@ -1106 +905 @@
                 Gquark->Set4Momentum(PkinkQ1);
                 Ganti_quark->Set4Momentum(PkinkQ2);
-/*
-G4cout<<" Proj Meson end Qbar"<<G4endl;
-G4cout<<"First string  ============ "<<G4endl;
-G4cout<<"end  Q "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl;
-G4cout<<"G     Q"<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl;
-G4cout<<"Sum P  "<<(Gquark->Get4Momentum()+end->Get4Momentum())<<G4endl;
-G4cout<<"Secondstring   ============ "<<G4endl;
-G4cout<<"G antQ "<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl;
-G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl;
-G4cout<<"Sum P  "<<(Ganti_quark->Get4Momentum()+start->Get4Momentum())<<G4endl;
-*/
+
               }   // end of if(end->GetPDGcode() > 0)
             } else {                      // Target
@@ -1125 +914 @@
                 Ganti_quark->Set4Momentum(PkinkQ2);
                 Gquark->Set4Momentum(PkinkQ1);
-/*
-G4cout<<" Targ Meson end Q"<<G4endl;
-G4cout<<"First string   ============ "<<G4endl;
-G4cout<<"G antiQ"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl;
-G4cout<<"end  Q "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl;
-G4cout<<"Sum P  "<<(Ganti_quark->Get4Momentum()+end->Get4Momentum())<<G4endl;
-G4cout<<"Secondstring   ============ "<<G4endl;
-G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl;
-G4cout<<"G    Q "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl;
-G4cout<<"Sum P  "<<(Gquark->Get4Momentum()+start->Get4Momentum())<<G4endl;
-*/
+
               } else
               {                            // Anti_Quark on the end
@@ -1142 +921 @@
                 Gquark->Set4Momentum(PkinkQ2);
                 Ganti_quark->Set4Momentum(PkinkQ1);
-/*
-G4cout<<" Targ Meson end Qbar"<<G4endl;
-G4cout<<"First string   ============ "<<G4endl;
-G4cout<<"G     Q"<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl;
-G4cout<<"end  Q "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl;
-G4cout<<"Sum P  "<<(Gquark->Get4Momentum()+end->Get4Momentum())<<G4endl;
-G4cout<<"Secondstring   ============ "<<G4endl;
-G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl;
-G4cout<<"G antQ "<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl;
-G4cout<<"Sum P  "<<(Ganti_quark->Get4Momentum()+start->Get4Momentum())<<G4endl;
-*/
+
               }   // end of if(end->GetPDGcode() > 0)
             }     // end of if ( isProjectile )
@@ -1166 +935 @@
                 Gquark->Set4Momentum(PkinkQ1);
                 Ganti_quark->Set4Momentum(PkinkQ2);
-/*
-G4cout<<" Proj baryon end QQ"<<G4endl;
-G4cout<<"First string   ============ "<<G4endl;
-G4cout<<"end OQ "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl;
-G4cout<<"G    Q"<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl;
-G4cout<<"Sum P  "<<(Gquark->Get4Momentum()+end->Get4Momentum())<<G4endl;
-G4cout<<"Secondstring   ============ "<<G4endl;
-G4cout<<"G  Qbar"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl;
-G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl;
-G4cout<<"Sum P  "<<(Ganti_quark->Get4Momentum()+start->Get4Momentum())<<G4endl;
-*/
+
               } else
               {                            // Anti_DiQuark on the end or quark
@@ -1183 +942 @@
                 Ganti_quark->Set4Momentum(PkinkQ1);
                 Gquark->Set4Momentum(PkinkQ2);
-/*
-G4cout<<" Proj baryon end Q"<<G4endl;
-G4cout<<"First string   ============ "<<G4endl;
-G4cout<<"end OQ "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl;
-G4cout<<"G antQ"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl;
-G4cout<<"Sum P  "<<(Ganti_quark->Get4Momentum()+end->Get4Momentum())<<G4endl;
-G4cout<<"Secondstring   ============ "<<G4endl;
-G4cout<<"G  Q   "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl;
-G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl;
-G4cout<<"Sum P  "<<(Gquark->Get4Momentum()+start->Get4Momentum())<<G4endl;
-*/
+
               }   // end of if(end->GetPDGcode() > 0)
             } else {                      // Target
@@ -1205 +954 @@
                 Gquark->Set4Momentum(PkinkQ1);
                 Ganti_quark->Set4Momentum(PkinkQ2);
-/*
-G4cout<<" Targ baryon end QQ"<<G4endl;
-G4cout<<"First string   ============ "<<G4endl;
-G4cout<<"G  Q   "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl;
-G4cout<<"end OQ "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl;
-G4cout<<"Sum P  "<<(Gquark->Get4Momentum()+end->Get4Momentum())<<G4endl;
-G4cout<<"Secondstring   ============ "<<G4endl;
-G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl;
-G4cout<<"G  Qbar"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl;
-G4cout<<"Sum P  "<<(Ganti_quark->Get4Momentum()+start->Get4Momentum())<<G4endl;
-*/
+
               } else
               {                            // Anti_DiQuark on the end or Q
@@ -1222 +961 @@
                 Gquark->Set4Momentum(PkinkQ2);
                 Ganti_quark->Set4Momentum(PkinkQ1);
-/*
-G4cout<<" Targ baryon end Q"<<G4endl;
-G4cout<<"First string   ============ "<<G4endl;
-G4cout<<"G  Qbar"<<Ganti_quark->GetDefinition()->GetPDGEncoding()<<" "<<Ganti_quark->Get4Momentum()<<G4endl;
-G4cout<<"end O "<<end->GetDefinition()->GetPDGEncoding()<<" "<<end->Get4Momentum()<<G4endl;
-G4cout<<"Sum P  "<<(Ganti_quark->Get4Momentum()+end->Get4Momentum())<<G4endl;
-G4cout<<"Secondstring   ============ "<<G4endl;
-G4cout<<"startQ "<<start->GetDefinition()->GetPDGEncoding()<<" "<<start->Get4Momentum()<<G4endl;
-G4cout<<"G  Q   "<<Gquark->GetDefinition()->GetPDGEncoding()<<" "<<Gquark->Get4Momentum()<<G4endl;
-G4cout<<"Sum P  "<<(Gquark->Get4Momentum()+start->Get4Momentum())<<G4endl;
-*/
+
               }   // end of if(end->GetPDGcode() > 0)
             }     // end of if ( isProjectile )
@@ -1346 +1075 @@
 {            //  @@ this method is used in FTFModel as well. Should go somewhere common!
 
-        G4double Pt2;
-        Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() * 
-                           (std::exp(-maxPtSquare/AveragePt2)-1.));
-
+        G4double Pt2(0.);
+        if(AveragePt2 <= 0.) {Pt2=0.;}
+        else
+        {
+         Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() * 
+                            (std::exp(-maxPtSquare/AveragePt2)-1.));
+        }
         G4double Pt=std::sqrt(Pt2);
         G4double phi=G4UniformRand() * twopi;

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

@@ -26 +26 @@
 //
 // $Id: G4DiffractiveSplitableHadron.cc,v 1.8 2009/07/31 11:03:00 vuzhinsk Exp $
-// GEANT4 tag $Name: geant4-09-03-cand-01 $
+// GEANT4 tag $Name: geant4-09-03 $
 //
 

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

@@ -25 +25 @@
 //
 //
-// $Id: G4ElasticHNScattering.cc,v 1.11 2009/10/06 10:10:36 vuzhinsk Exp $
+// $Id: G4ElasticHNScattering.cc,v 1.14 2009/12/16 17:51:13 gunter Exp $
 // ------------------------------------------------------------
 //      GEANT 4 class implemetation file
@@ -58 +58 @@
                               G4FTFParameters    *theParameters) const
 {
-//G4cout<<"G4ElasticHNScattering::ElasticScattering"<<G4endl;
 // -------------------- Projectile parameters -----------------------------------
            G4LorentzVector Pprojectile=projectile->Get4Momentum();
-//G4cout<<"Elastic P "<<Pprojectile<<G4endl;
+
            if(Pprojectile.z() < 0.)
            {
@@ -68 +67 @@
            } 
 
-//           G4double ProjectileRapidity = Pprojectile.rapidity();
-//           G4int    ProjectilePDGcode=projectile->GetDefinition()->GetPDGEncoding();
-//           G4ParticleDefinition * ProjectileDefinition = projectile->GetDefinition();
-
            G4bool PutOnMassShell(false);
 
@@ -86 +81 @@
 
 // -------------------- Target parameters ----------------------------------------------
-//           G4int    TargetPDGcode=target->GetDefinition()->GetPDGEncoding();
 
            G4LorentzVector Ptarget=target->Get4Momentum();
-//G4cout<<"Elastic T "<<Ptarget<<G4endl;
-//           G4double TargetRapidity = Ptarget.rapidity();
+
            G4double M0target = Ptarget.mag();
 
@@ -149 +142 @@
             else // if(M0projectile > projectile->GetDefinition()->GetPDGMass())
             {
-             target->SetStatus(2);                                   // Uzhi 17.07.09
+             target->SetStatus(2);                                   
              return false;                   // The projectile was not excited,
                                              // but the energy was too low to put
@@ -183 +176 @@
            G4double maxPtSquare = PZcms2;
 
-//----- Charge exchange between the projectile and the target, if possible 
-//        (G4UniformRand() < 0.5*std::exp(-0.5*(ProjectileRapidity - TargetRapidity))))
-/*
-           if((ProjectilePDGcode != TargetPDGcode) && 
-             ((ProjectilePDGcode > 1000) && (TargetPDGcode > 1000)) &&
-             (G4UniformRand() < 1.0*std::exp(-0.5*(ProjectileRapidity - TargetRapidity)))) 
-           {
-              projectile->SetDefinition(target->GetDefinition());
-              target->SetDefinition(ProjectileDefinition); 
-           }
-*/
 // ------ Now we can calculate the transfered Pt --------------------------
            G4double Pt2;                                                    
@@ -255 +237 @@
 {            //  @@ this method is used in FTFModel as well. Should go somewhere common!
         
-        G4double Pt2;
-        Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() * 
+        G4double Pt2(0.);
+        if(AveragePt2 <= 0.) {Pt2=0.;}
+        else
+        {
+         Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() * 
                 (std::exp(-maxPtSquare/AveragePt2)-1.)); 
-        
+        }
         G4double Pt=std::sqrt(Pt2);
         G4double phi=G4UniformRand() * twopi;

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

@@ -25 +25 @@
 //
 //
-// $Id: G4FTFModel.cc,v 1.28 2009/10/29 14:55:33 vuzhinsk Exp $
-// GEANT4 tag $Name: geant4-09-03-cand-01 $
+// $Id: G4FTFModel.cc,v 1.34 2009/12/15 19:14:31 vuzhinsk Exp $
+// GEANT4 tag $Name: geant4-09-03 $
 //
 
@@ -56 +56 @@
         G4VPartonStringModel::SetThisPointer(this);
         theParameters=0;
+        NumberOfInvolvedNucleon=0;
 }
 
@@ -61 +62 @@
 G4FTFModel::~G4FTFModel()
 {
-   if( theParameters != 0 ) delete theParameters; 
 // Because FTF model can be called for various particles
 // theParameters must be erased at the end of each call.
 // Thus the delete is also in G4FTFModel::GetStrings() method
+   if( theParameters != 0 ) delete theParameters; 
    if( theExcitation != 0 ) delete theExcitation;
    if( theElastic    != 0 ) delete theElastic; 
+
+   if( NumberOfInvolvedNucleon != 0)
+   {
+    for(G4int i=0; i < NumberOfInvolvedNucleon; i++)
+    {
+     G4VSplitableHadron * aNucleon = TheInvolvedNucleon[i]->GetSplitableHadron();
+     if(aNucleon) delete aNucleon;
+    }
+   }
 }
 
@@ -93 +103 @@
 
 // --- cms energy
-
         G4double s = sqr( theProjectile.GetMass() ) +
                      sqr( G4Proton::Proton()->GetPDGMass() ) +
                      2*theProjectile.GetTotalEnergy()*G4Proton::Proton()->GetPDGMass();
-/*
-G4cout<<"----------------------------------------"<<G4endl;
-G4cout << " primary Total E (GeV): " << theProjectile.GetTotalEnergy()/GeV << G4endl;
-G4cout << " primary Mass    (GeV): " << theProjectile.GetMass() /GeV << G4endl;
-G4cout << " primary 3Mom           " << theProjectile.GetMomentum() << G4endl;
-G4cout << " primary space position " << theProjectile.GetPositionInNucleus() << G4endl;
-G4cout << "cms std::sqrt(s) (GeV) = " << std::sqrt(s) / GeV << G4endl;
-*/
 
       if( theParameters != 0 ) delete theParameters;
@@ -110 +111 @@
                                           aNucleus.GetN(),aNucleus.GetZ(),
                                           s);
+
 //theParameters->SetProbabilityOfElasticScatt(0.); 
 // To turn on/off (1/0) elastic scattering
@@ -116 +118 @@
 
 // ------------------------------------------------------------
+struct DeleteVSplitableHadron { void operator()(G4VSplitableHadron * aH){ delete aH;} };
+
+
+// ------------------------------------------------------------
 G4ExcitedStringVector * G4FTFModel::GetStrings()
-{
-//G4int Uzhi; G4cin>>Uzhi;
- 
-//G4cout<<"GetList"<<G4endl;
+{ 
+        G4ExcitedStringVector * theStrings(0);
+
         theParticipants.GetList(theProjectile,theParameters);
-//G4cout<<"Regge"<<G4endl;
-        ReggeonCascade();                                     // Uzhi 26 July 09 
-//G4cout<<"On mass shell"<<G4endl;
-        if (! PutOnMassShell()    ) return NULL;              // Uzhi 26 July 09
-//G4cout<<"Excite"<<G4endl;
-        if (! ExciteParticipants()) return NULL;
-//G4cout<<"Strings"<<G4endl;
-        G4ExcitedStringVector * theStrings = BuildStrings();
-//G4cout<<"Out FTF N strings "<<theStrings->size()<<G4endl;
-//G4cout<<"GetResidualNucleus"<<G4endl;
-        GetResidualNucleus();
-//G4cout<<"Out of FTF"<<G4endl;
-        if( theParameters != 0 )
-        {
-          delete theParameters;
-          theParameters=0;
+
+        ReggeonCascade(); 
+
+        G4bool Success(true);
+        if( PutOnMassShell() )
+        {
+         if( ExciteParticipants() )
+         {
+          theStrings = BuildStrings();
+
+          GetResidualNucleus();
+
+          if( theParameters != 0 )
+          {
+           delete theParameters;
+           theParameters=0;
+          }
+         } else                      // if( ExciteParticipants() )
+         {     Success=false;}
+        } else                       // if( PutOnMassShell() )
+        {      Success=false;}
+
+        if(!Success)   
+        {
+           // -------------- Erase the projectile ----------------
+         std::vector<G4VSplitableHadron *> primaries;
+
+         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());
+         }
+         std::for_each(primaries.begin(), primaries.end(), DeleteVSplitableHadron());
+         primaries.clear();
         }
-        return theStrings;
-}
-
-// ------------------------------------------------------------
-struct DeleteVSplitableHadron { void operator()(G4VSplitableHadron * aH){ delete aH;} };
-
-// ------------------------------------------------------------
-void G4FTFModel::ReggeonCascade()                             // Uzhi 26 July 2009
+
+// -------------- Cleaning of the memory --------------
+// -------------- Erase the target nucleons -----------
+        G4VSplitableHadron * aNucleon = 0;
+        for(G4int i=0; i < NumberOfInvolvedNucleon; i++)
+        {
+           aNucleon = TheInvolvedNucleon[i]->GetSplitableHadron();
+           if(aNucleon) delete aNucleon;
+        } 
+
+        NumberOfInvolvedNucleon=0;
+
+        return theStrings;
+
+}
+//-------------------------------------------------------------------
+void G4FTFModel::ReggeonCascade()                             
 { //--- Implementation of reggeon theory inspired model-------
         NumberOfInvolvedNucleon=0;
 
-//G4int PrimInt(0);
         theParticipants.StartLoop();
         while (theParticipants.Next())
-        {
-//PrimInt++;       
+        {   
            const G4InteractionContent & collision=theParticipants.GetInteraction();
            G4Nucleon * TargetNucleon=collision.GetTargetNucleon();
-//G4cout<<"Prim Nnucl "<<TargetNucleon->Get4Momentum()<<G4endl;
+
            TheInvolvedNucleon[NumberOfInvolvedNucleon]=TargetNucleon;
            NumberOfInvolvedNucleon++;
@@ -165 +199 @@
 
            theParticipants.theNucleus->StartLoop();
-           G4Nucleon * Neighbour;
-//G4int NucleonNum(0);
+           G4Nucleon * Neighbour(0);
+
            while ( (Neighbour = theParticipants.theNucleus->GetNextNucleon()) )
            {
@@ -177 +211 @@
                 std::exp(-impact2/theParameters->GetR2ofNuclearDestruction()))  
              { // The neighbour nucleon is involved in the reggeon cascade
-//G4cout<<" involved "<<NucleonNum<<" "<<Neighbour->Get4Momentum()<<G4endl;
+
               TheInvolvedNucleon[NumberOfInvolvedNucleon]=Neighbour;
               NumberOfInvolvedNucleon++;
-//G4cout<<" PrimInt"<<" "<<NumberOfInvolvedNucleon<<G4endl;
 
               G4VSplitableHadron *targetSplitable; 
@@ -186 +219 @@
 
               Neighbour->Hit(targetSplitable);
-//              Neighbour->SetBindingEnergy(3.*Neighbour->GetBindingEnergy()); // Uzhi 5.10.09
               targetSplitable->SetStatus(2);     
              }
             }  // end of if(!Neighbour->AreYouHit())
-//NucleonNum++;
            }   // end of while (theParticipant.theNucleus->GetNextNucleon())
-//G4cout<<"Prim Int N Ninvolv "<<PrimInt<<" "<<NumberOfInvolvedNucleon<<G4endl;
         }      // end of while (theParticipants.Next())
-//G4cout<<"At end "<<PrimInt<<" "<<NumberOfInvolvedNucleon<<G4endl;
 
 // ---------------- Calculation of creation time for each target nucleon -----------
@@ -240 +269 @@
 {
 // -------------- Properties of the projectile ----------------
-//G4cout<<"Put on Mass-shell"<<G4endl;
         theParticipants.StartLoop();    // restart a loop
         theParticipants.Next();
         G4VSplitableHadron * primary = theParticipants.GetInteraction().GetProjectile();
         G4LorentzVector Pprojectile=primary->Get4Momentum();
-//G4cout<<"Proj "<<Pprojectile<<G4endl;
+
+// To get original projectile particle
+
         if(Pprojectile.z() < 0.){return false;}
 
@@ -252 +282 @@
 //-------------------------------------------------------------
         G4LorentzVector Psum      = Pprojectile;
-        G4double        SumMasses = Mprojectile;
+        G4double        SumMasses = Mprojectile + 20.*MeV; // 13.12.09
+                                               // Separation energy for projectile
 
 //--------------- Target nucleus ------------------------------
@@ -259 +290 @@
         G4int ResidualMassNumber=theNucleus->GetMassNumber();
         G4int ResidualCharge    =theNucleus->GetCharge();
-//G4cout<<"Nucleus "<<ResidualMassNumber<<" "<<ResidualCharge<<G4endl;
         ResidualExcitationEnergy=0.;
         G4LorentzVector PnuclearResidual(0.,0.,0.,0.);
@@ -267 +297 @@
 
         theNucleus->StartLoop();
-//G4int NucleonNum(0);
+
         while ((aNucleon = theNucleus->GetNextNucleon()))
         {
          if(aNucleon->AreYouHit())
          {   // Involved nucleons
-//G4cout<<"          "<<NucleonNum<<" "<<aNucleon->Get4Momentum()<<G4endl;
           Psum += aNucleon->Get4Momentum();
           SumMasses += aNucleon->GetDefinition()->GetPDGMass();
+          SumMasses += 20.*MeV;   // 13.12.09 Separation energy for a nucleon
           ResidualMassNumber--;
           ResidualCharge-=(G4int) aNucleon->GetDefinition()->GetPDGCharge();
@@ -283 +313 @@
           PnuclearResidual += aNucleon->Get4Momentum();
          }  // end of if(!aNucleon->AreYouHit())
-//NucleonNum++;
         }   // end of while (theNucleus->GetNextNucleon())
 
-//G4cout<<"Nucleus "<<ResidualMassNumber<<" "<<ResidualCharge<<G4endl;
-//G4cout<<"PResid "<<PnuclearResidual<<G4endl;
-
         Psum += PnuclearResidual;
+
         G4double ResidualMass(0.);
         if(ResidualMassNumber == 0)
@@ -303 +330 @@
         }
  
-//        ResidualMass +=ResidualExcitationEnergy; // Will be given after checks
+//      ResidualMass +=ResidualExcitationEnergy; // Will be given after checks
         SumMasses += ResidualMass;
 
@@ -310 +337 @@
         G4double SqrtS=Psum.mag();
         G4double     S=Psum.mag2();
-//G4cout<<" SqrtS SumMasses "<<SqrtS <<" "<< SumMasses<<G4endl;
-//G4cout<<"Res M E* "<<ResidualMass<<" "<<ResidualExcitationEnergy<<G4endl;
+
         if(SqrtS < SumMasses)      {return false;} // It is impossible to simulate
                                                    // after putting nuclear nucleons
                                                    // on mass-shell
+
         if(SqrtS < SumMasses+ResidualExcitationEnergy) {ResidualExcitationEnergy=0.;}
 
         ResidualMass +=ResidualExcitationEnergy;
         SumMasses    +=ResidualExcitationEnergy;
-//G4cout<<"Res M E* "<<ResidualMass<<" "<<ResidualExcitationEnergy<<G4endl;
 
 //-------------------------------------------------------------
@@ -325 +351 @@
         G4int MaxNumberOfDeltas = (int)((SqrtS - SumMasses)/(400.*MeV));
         G4int NumberOfDeltas(0);
-//SumMasses=Mprojectile;
+
         if(theNucleus->GetMassNumber() != 1)
         {
-          G4double ProbDeltaIsobar(0.);  // 1. *******************************
+          G4double ProbDeltaIsobar(0.);  // 1. *** Can be set if it is needed
           for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
           {
@@ -343 +369 @@
               targetSplitable->SetDefinition(ptr);
               SumMasses+=targetSplitable->GetDefinition()->GetPDGMass();
-//G4cout<<" i "<<i<<" Delta "<<targetSplitable->GetDefinition()->GetPDGMass()<<G4endl;
             } 
-            else 
-            {
-//              SumMasses+=TheInvolvedNucleon[i]->GetSplitableHadron()->
-//                         GetDefinition()->GetPDGMass();
-//G4cout<<" i "<<i<<" Nuclo "<<TheInvolvedNucleon[i]->GetSplitableHadron()->GetDefinition()->GetPDGMass()<<G4endl;
-            }
           }   // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
         }   // end of if(theNucleus.GetMassNumber() != 1)
-//G4cout<<"MaxNumberOfDeltas NumberOfDeltas "<<MaxNumberOfDeltas<<" "<<NumberOfDeltas<<G4endl;
-//G4cout<<" SqrtS SumMasses "<<SqrtS <<" "<< SumMasses<<G4endl;
-//-------------------------------------------------------------
-
+//-------------------------------------------------------------
         G4LorentzRotation toCms(-1*Psum.boostVector());
         G4LorentzVector Ptmp=toCms*Pprojectile;
-
         if ( Ptmp.pz() <= 0. )                                
         {  // "String" moving backwards in  CMS, abort collision !!
@@ -366 +381 @@
         }
 
-        toCms.rotateZ(-1*Ptmp.phi());
-        toCms.rotateY(-1*Ptmp.theta());
+//        toCms.rotateZ(-1*Ptmp.phi());              // Uzhi 5.12.09
+//        toCms.rotateY(-1*Ptmp.theta());            // Uzhi 5.12.09
         
         G4LorentzRotation toLab(toCms.inverse());
 
-//        Pprojectile.transform(toCms);
-//G4cout<<"Proj in CMS "<<Pprojectile<<G4endl;
-
-//G4cout<<"Main work"<<G4endl;
 //-------------------------------------------------------------
 //------- Ascribing of the involved nucleons Pt and Xminus ----
         G4double Dcor        = theParameters->GetDofNuclearDestruction()/
                                                theNucleus->GetMassNumber();
-//                                                  NumberOfInvolvedNucleon;
+
         G4double AveragePt2  = theParameters->GetPt2ofNuclearDestruction();
         G4double maxPtSquare = theParameters->GetMaxPt2ofNuclearDestruction();
-//G4cout<<" D Pt2 "<<Dcor<<" "<<AveragePt2<<G4endl;
 
         G4double M2target(0.);
+        G4double WminusTarget(0.);
+        G4double WplusProjectile(0.);
+
         G4int    NumberOfTries(0);
         G4double ScaleFactor(1.);
-        do    // while (SqrtS < Mprojectile + std::sqrt(M2target))
-        {     // while (DecayMomentum < 0.)
-
-          NumberOfTries++;
-//G4cout<<"NumberOfTries "<<NumberOfTries<<G4endl;
-          if(NumberOfTries == 100*(NumberOfTries/100))
-          { // At large number of tries it would be better to reduce the values
-            ScaleFactor/=2.;
-            Dcor       *=ScaleFactor;
-            AveragePt2 *=ScaleFactor;
-//G4cout<<"NumberOfTries "<<NumberOfTries<<" "<<Dcor<<" "<<AveragePt2<<G4endl;
-//G4int Uzhi; G4cin>>Uzhi;
-          }
-//G4cout<<"Start Decay "<<G4endl; G4int Uzhi; G4cin>>Uzhi;
-          G4ThreeVector PtSum(0.,0.,0.);
-          G4double XminusSum(0.);
-          G4double Xminus(0.);
-          G4bool Success=true;
-
-          do      // while(Success == false);
-          {
-//G4cout<<"Sample Pt and X"<<" Ninv "<<NumberOfInvolvedNucleon<<G4endl; // G4int Uzhi1; G4cin>>Uzhi1;
-             Success=true;
+        G4bool OuterSuccess(true);
+        do    // while (!OuterSuccess)
+        {
+          OuterSuccess=true;
+
+          do    // while (SqrtS < Mprojectile + std::sqrt(M2target))
+          {     // while (DecayMomentum < 0.)
+
+            NumberOfTries++;
+            if(NumberOfTries == 100*(NumberOfTries/100))   // 100
+            { // At large number of tries it would be better to reduce the values
+              ScaleFactor/=2.;
+              Dcor       *=ScaleFactor;
+              AveragePt2 *=ScaleFactor;
+            }
+
+            G4ThreeVector PtSum(0.,0.,0.);
+            G4double XminusSum(0.);
+            G4double Xminus(0.);
+            G4bool InerSuccess=true;
+
+            do      // while(!InerSuccess);
+            {
+             InerSuccess=true;
 
              PtSum    =G4ThreeVector(0.,0.,0.);
@@ -420 +434 @@
                G4ThreeVector tmpPt = GaussianPt(AveragePt2, maxPtSquare);
                PtSum += tmpPt;
-
                G4ThreeVector tmpX=GaussianPt(Dcor*Dcor, 1.);
                Xminus=tmpX.x();
@@ -427 +440 @@
                G4LorentzVector tmp(tmpPt.x(),tmpPt.y(),Xminus,0.);
                aNucleon->SetMomentum(tmp);
-//G4cout<<"Nucl mom "<<aNucleon->GetMomentum()<<G4endl;
              }   // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
 
@@ -445 +457 @@
               DeltaXminus = -1./theNucleus->GetMassNumber();
              }
-//G4cout<<"Deltas "<<DeltaX<<" "<<DeltaY<<" "<<DeltaXminus<<G4endl;
 
              XminusSum=1.;
              M2target =0.;
-
 
              for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
@@ -456 +466 @@
 
                Xminus = aNucleon->Get4Momentum().pz() - DeltaXminus;
-//G4cout<<i<<" "<<Xminus<<" "<<XminusSum<<G4endl;
                XminusSum-=Xminus;               
+
                if((Xminus <= 0.)   || (Xminus > 1.) || 
-                  (XminusSum <=0.) || (XminusSum > 1.)) {Success=false; break;}
+                  (XminusSum <=0.) || (XminusSum > 1.)) {InerSuccess=false; break;}
  
                G4double Px=aNucleon->Get4Momentum().px() - DeltaX;
@@ -472 +482 @@
              }   // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
 
-             if(Success && (ResidualMassNumber != 0))
+             if(InerSuccess && (ResidualMassNumber != 0))
              {
               M2target +=(ResidualMass*ResidualMass + PtSum.mag2())/XminusSum;
              }
-//G4cout<<"Success "<<Success<<G4endl;
-//G4int Uzhi; G4cin>>Uzhi;
-          } while(Success == 0);
-
-//-------------------------------------------------------------
-//G4cout<<"SqrtS > Mprojectile + std::sqrt(M2target) "<<SqrtS<<" "<<Mprojectile<<" "<< std::sqrt(M2target)<<" "<<Mprojectile + std::sqrt(M2target)<<G4endl;
-//G4int Uzhi3; G4cin>>Uzhi3;
-        } while (SqrtS < Mprojectile + std::sqrt(M2target));
-
-//-------------------------------------------------------------
-        G4double DecayMomentum2= S*S+M2projectile*M2projectile+M2target*M2target
-                                -2.*S*M2projectile - 2.*S*M2target 
-                                       -2.*M2projectile*M2target;
-
-        G4double WminusTarget=(S-M2projectile+M2target+std::sqrt(DecayMomentum2))/2./SqrtS;
-        G4double WplusProjectile=SqrtS - M2target/WminusTarget;
-//G4cout<<"DecM W+ W- "<<DecayMomentum2<<" "<<WplusProjectile<<" "<<WminusTarget<<G4endl;
-//-------------------------------------------------------------
-        G4double Pzprojectile=WplusProjectile/2. - M2projectile/2./WplusProjectile;
-        G4double Eprojectile =WplusProjectile/2. + M2projectile/2./WplusProjectile;
-        Pprojectile.setPz(Pzprojectile);  Pprojectile.setE(Eprojectile);
-
-        Pprojectile.transform(toLab);       // The work with the projectile
-        primary->Set4Momentum(Pprojectile); // is finished at the moment.
-//G4cout<<"Proj After "<<Pprojectile<<G4endl;
-//-------------------------------------------------------------
-//Ninvolv=0;
-
-        G4ThreeVector Residual3Momentum(0.,0.,1.);
-
-        for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
-        {
+            } while(!InerSuccess);
+          } while (SqrtS < Mprojectile + std::sqrt(M2target));
+//-------------------------------------------------------------
+          G4double DecayMomentum2= S*S+M2projectile*M2projectile+M2target*M2target
+                                    -2.*S*M2projectile - 2.*S*M2target 
+                                         -2.*M2projectile*M2target;
+
+          WminusTarget=(S-M2projectile+M2target+std::sqrt(DecayMomentum2))/2./SqrtS;
+          WplusProjectile=SqrtS - M2target/WminusTarget;
+//-------------------------------------------------------------
+          for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
+          {
            G4Nucleon * aNucleon = TheInvolvedNucleon[i];
            G4LorentzVector tmp=aNucleon->Get4Momentum();
-           Residual3Momentum-=tmp.vect();
 
            G4double Mt2 = sqr(tmp.x())+sqr(tmp.y())+
@@ -520 +509 @@
            G4double E = WminusTarget*Xminus/2. + Mt2/(2.*WminusTarget*Xminus);
 
+           if( E+Pz > WplusProjectile ){OuterSuccess=false; break;}
+          }   // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
+        } while(!OuterSuccess);
+
+//-------------------------------------------------------------
+        G4double Pzprojectile=WplusProjectile/2. - M2projectile/2./WplusProjectile;
+        G4double Eprojectile =WplusProjectile/2. + M2projectile/2./WplusProjectile;
+        Pprojectile.setPz(Pzprojectile);  Pprojectile.setE(Eprojectile);
+
+        Pprojectile.transform(toLab);       // The work with the projectile
+        primary->Set4Momentum(Pprojectile); // is finished at the moment.
+
+//-------------------------------------------------------------
+        G4ThreeVector Residual3Momentum(0.,0.,1.);
+
+        for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
+        {
+           G4Nucleon * aNucleon = TheInvolvedNucleon[i];
+           G4LorentzVector tmp=aNucleon->Get4Momentum();
+           Residual3Momentum-=tmp.vect();
+
+           G4double Mt2 = sqr(tmp.x())+sqr(tmp.y())+
+                          aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass()*
+                          aNucleon->GetSplitableHadron()->GetDefinition()->GetPDGMass();
+           G4double Xminus=tmp.z();
+
+           G4double Pz=-WminusTarget*Xminus/2. + Mt2/(2.*WminusTarget*Xminus);
+           G4double E = WminusTarget*Xminus/2. + Mt2/(2.*WminusTarget*Xminus);
+
            tmp.setPz(Pz); 
            tmp.setE(E);
+
            tmp.transform(toLab);
-//G4cout<<"invol  "<<Ninvolv<<" "<<tmp<<G4endl;
-//Ninvolv++;
+
            aNucleon->SetMomentum(tmp);
+
            G4VSplitableHadron * targetSplitable=aNucleon->GetSplitableHadron();
            targetSplitable->Set4Momentum(tmp);
-//G4cout<<"nucleon M "<<aNucleon->Get4Momentum()<<G4endl;
            
         }   // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
 
         G4double Mt2Residual=sqr(ResidualMass) +
-                                 sqr(Residual3Momentum.x())+sqr(Residual3Momentum.x());
+                                 sqr(Residual3Momentum.x())+sqr(Residual3Momentum.y());
 
         G4double PzResidual=-WminusTarget*Residual3Momentum.z()/2. + 
@@ -540 +558 @@
                              Mt2Residual/(2.*WminusTarget*Residual3Momentum.z());
 
-//        G4LorentzVector Residual4Momentum(0.,0.,0.,0.);
         Residual4Momentum.setPx(Residual3Momentum.x());
         Residual4Momentum.setPy(Residual3Momentum.y());
         Residual4Momentum.setPz(PzResidual); 
         Residual4Momentum.setE(EResidual);
+
         Residual4Momentum.transform(toLab);
-
-//G4cout<<"Return Nucleus"<<G4endl;
-//-------------------------------------------------------------
-//-------------------------------------------------------------
-//-------------------------------------------------------------
-//G4int Uzhi2; G4cin>>Uzhi2;
-
+//-------------------------------------------------------------
  return true;
 }
@@ -559 +571 @@
 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
-
-
-//G4int InterNumber=0; // Vova
-
-        G4bool Successfull=false;
-        theParticipants.StartLoop();                          //Uzhi 26 July 09
-
-//      while (theParticipants.Next()&& (InterNumber < 3)) // Vova
+    G4bool Successfull(false);
+//    do {                           // } while (Successfull == false) // Closed 15.12.09
+        Successfull=false;
+        theParticipants.StartLoop();
         while (theParticipants.Next())
         {          
            const G4InteractionContent & collision=theParticipants.GetInteraction();
-//
-//counter++;
-//G4cout<<" Int num "<<counter<<G4endl;
-//
+
            G4VSplitableHadron * projectile=collision.GetProjectile();
            G4VSplitableHadron * target=collision.GetTarget();
-//         G4Nucleon * TargetNucleon=collision.GetTargetNucleon(); // Uzhi 16.07.09
-// Uzhi 16.07.09 ----------------------------
+
            if(G4UniformRand()< theParameters->GetProbabilityOfElasticScatt())
            { //   Elastic scattering -------------------------
-//G4cout<<"Elastic"<<G4endl;
             if(theElastic->ElasticScattering(projectile, target, theParameters))
             {
-             Successfull = Successfull || true;
+            Successfull = Successfull || true;
             } else
             {
-//G4cout<<"Elastic Not succes"<<G4endl;
              Successfull = Successfull || false;
              target->SetStatus(2);
-/*
-             if(target->GetStatus() == 0)                         // Uzhi 17.07.09
-             {
-              G4VSplitableHadron * aHit=0;                        // Uzhi 16.07.09
-              TargetNucleon->Hit(aHit);                           // Uzhi 16.07.09
-             };
-*/
-            };
+            }
            }
            else
            { //   Inelastic scattering ---------------------- 
-//G4cout<<"InElastic"<<G4endl;
             if(theExcitation->ExciteParticipants(projectile, target, 
                                                  theParameters, theElastic))
@@ -614 +601 @@
             } else
             {
-//G4cout<<"InElastic Non succes"<<G4endl;
              Successfull = Successfull || false;
              target->SetStatus(2);
-/*
-             if(target->GetStatus() == 0)                         // Uzhi 16.06.09
-             {
-              G4VSplitableHadron * aHit=0;                        // Uzhi 16.07.09
-              TargetNucleon->Hit(aHit);                           // Uzhi 16.07.09
-             };
-*/
-            };
+            }
            }
-        }       // end of the loop Uzhi 9.07.09
-// Uzhi 16.07.09 ----------------------------
-
-        if(!Successfull)
-        {
-//G4cout<<"Process not successfull"<<G4endl;
-//           give up, clean up
-          std::vector<G4VSplitableHadron *> primaries;
-//        std::vector<G4VSplitableHadron *> targets;        // Uzhi 31.07.09
-          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());
-/*  // Uzhi 31.07.09
-            if ( targets.end()   == std::find(targets.begin(), targets.end(),
-                                                      interaction.GetTarget()) ) 
-                targets.push_back(interaction.GetTarget());
-*/  // Uzhi 31.07.09
-          }
-          std::for_each(primaries.begin(), primaries.end(), DeleteVSplitableHadron());
-          primaries.clear();
-/*  // Uzhi 31.07.09    
-          std::for_each(targets.begin(), targets.end(), DeleteVSplitableHadron());
-          targets.clear();
-*/  // Uzhi 31.07.09
-//          theParticipants.theNucleus->StartLoop();
-
-//G4cout<<"NumberOfInvolvedNucleon "<<NumberOfInvolvedNucleon<<G4endl;
-          G4VSplitableHadron * aNucleon = 0;
-          for(G4int i=0; i < NumberOfInvolvedNucleon; i++)
-          {
-           aNucleon = TheInvolvedNucleon[i]->GetSplitableHadron();
-           if(aNucleon)
-           { 
-             if(aNucleon->GetStatus() != 0) delete aNucleon;
-//           if(aNucleon->GetStatus() == 2)  DeleteVSplitableHadron()(aNucleon);
-           }
-          } 
-
-          NumberOfInvolvedNucleon=0;
-//G4cout<<"Out of Excit"<<G4endl; G4int Uzhi; G4cin>>Uzhi;
-          return false;
-        }  // End of if(!Successfull)
-
-        return true;
+        }       // end of while (theParticipants.Next())
+//       } while (Successfull == false);                        // Closed 15.12.09
+        return Successfull;
 }
 // ------------------------------------------------------------
@@ -682 +615 @@
 //  be duplicate in the List ( to unique G4VSplitableHadrons)
 
-//G4cout<<"In build string"<<G4endl;
-
         G4ExcitedStringVector * strings;
         strings = new G4ExcitedStringVector();
         
         std::vector<G4VSplitableHadron *> primaries;
-        std::vector<G4VSplitableHadron *> targets;
-        std::vector<G4Nucleon          *> TargetNucleons;     // Uzhi 16.07.09
         
         G4ExcitedString * FirstString(0);     // If there will be a kink,
-        G4ExcitedString * SecondString(0);    // two strings will be prodused.
+        G4ExcitedString * SecondString(0);    // two strings will be produced.
 
         theParticipants.StartLoop();    // restart a loop 
-//G4int InterCount(0); // Uzhi
+
         while ( theParticipants.Next() ) 
         {
@@ -703 +632 @@
             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());
-
-            if ( TargetNucleons.end()   == std::find(TargetNucleons.begin(),     
-                                                     TargetNucleons.end(),       
-                                                interaction.GetTargetNucleon()) )
-                TargetNucleons.push_back(interaction.GetTargetNucleon());        
-//InterCount++;
+                primaries.push_back(interaction.GetProjectile());     
         }
             
-        
-//G4cout << "BuildStrings prim/targ/woundN " << primaries.size() << " , " <<targets.size() <<", "<<TargetNucleons.size()<< G4endl;
-
         unsigned int ahadron;
-// Only for hA-interactions Uzhi -------------------------------------
-//G4int StringN(0);
-//G4cout<<"Proj strings -----------------------"<<G4endl;
         for ( ahadron=0; ahadron < primaries.size() ; ahadron++)
         {
-//G4cout<<" string# "<<ahadron<<" "<<primaries[ahadron]->Get4Momentum()<<G4endl;
             G4bool isProjectile(0);
             if(primaries[ahadron]->GetStatus() == 1) {isProjectile=true; }
@@ -734 +646 @@
                                          FirstString, SecondString,
                                          theParameters);
-//G4cout<<"1str 2str "<<FirstString<<" "<<SecondString<<G4endl;
+
             if(FirstString  != 0) strings->push_back(FirstString);
             if(SecondString != 0) strings->push_back(SecondString);
-
-//StringN++; G4cout<<"Proj string "<<strings->size()<<G4endl;
         }
-//G4cout<<"Targ strings ------------------------------"<<G4endl;
-        for ( ahadron=0; ahadron < targets.size() ; ahadron++)
-        {
-//G4cout<<"targets[ahadron]->GetStatus() "<<targets[ahadron]->GetStatus()<<G4endl;
-            if(targets[ahadron]->GetStatus() == 1)   // Uzhi 17.07.09
-            {
-             G4bool isProjectile=false;
-             FirstString=0; SecondString=0;
-             theExcitation->CreateStrings(targets[ahadron], isProjectile,
-                                          FirstString, SecondString,
-                                          theParameters);
-             if(FirstString  != 0) strings->push_back(FirstString);
-             if(SecondString != 0) strings->push_back(SecondString);
-
-//StringN++; G4cout<<"Targ string"<<StringN<<G4endl;
-            } else
-            {
-             if(targets[ahadron]->GetStatus() == 0)// Uzhi 17.07.09 Nucleon was rejected
-             {
-              G4VSplitableHadron * aHit=0;          // Uzhi 16.07.09
-              TargetNucleons[ahadron]->Hit(aHit);   // Uzhi 16.07.09
-             }
-            }
-        }
-
-//G4cout<<"Proj + Targ string "<<strings->size()<<G4endl;
-//G4cout<<"Involv strings NumberOfInvolvedNucleon "<<NumberOfInvolvedNucleon<<G4endl;
+
         for (G4int ahadron=0; ahadron < NumberOfInvolvedNucleon ; ahadron++)
         {
-/*
-G4cout<<"ahadron "<<ahadron<<" Status "<<
-TheInvolvedNucleon[ahadron]->GetSplitableHadron()->GetStatus()<<
-TheInvolvedNucleon[ahadron]->GetSplitableHadron()->Get4Momentum()<<G4endl;
-*/
-            if(TheInvolvedNucleon[ahadron]->GetSplitableHadron()->GetStatus() == 2)
-            {
-//StringN++; G4cout<<"Invol string "<<StringN<<G4endl;
+            if(TheInvolvedNucleon[ahadron]->GetSplitableHadron()->GetStatus() !=0) //== 2)
+            {
              G4bool isProjectile=false;
              FirstString=0; SecondString=0;
@@ -786 +664 @@
              if(FirstString  != 0) strings->push_back(FirstString);
              if(SecondString != 0) strings->push_back(SecondString);
-
-//           strings->push_back(theExcitation->String(
-//                      TheInvolvedNucleon[ahadron]->GetSplitableHadron(), 
-//                                                           isProjectile));
-            } 
-//G4cout<<"Proj + Targ+Involved string "<<strings->size()<<G4endl;
-/*
-else
-            {
-G4cout<<"Else ????????????"<<G4endl; 
-             if(targets[ahadron]->GetStatus() == 0)// Uzhi 17.07.09 Nucleon was rejected
-             {
-              G4VSplitableHadron * aHit=0;          // Uzhi 16.07.09
-              TargetNucleons[ahadron]->Hit(aHit);   // Uzhi 16.07.09
-             }
             }
-*/
         }
-
-//G4int Uzhi; G4cin>>Uzhi;
 
         std::for_each(primaries.begin(), primaries.end(), DeleteVSplitableHadron());
         primaries.clear();
-
-        std::for_each(targets.begin(), targets.end(), DeleteVSplitableHadron());
-        targets.clear();
 
         return strings;
@@ -826 +683 @@
         {
          G4Nucleon * aNucleon = TheInvolvedNucleon[i];
-         G4LorentzVector tmp=aNucleon->Get4Momentum()-DeltaPResidualNucleus;
+//         G4LorentzVector tmp=aNucleon->Get4Momentum()-DeltaPResidualNucleus;
+         G4LorentzVector tmp=-DeltaPResidualNucleus;
          aNucleon->SetMomentum(tmp);
          aNucleon->SetBindingEnergy(DeltaExcitationE);
@@ -837 +695 @@
 {            //  @@ this method is used in FTFModel as well. Should go somewhere common!
         
-        G4double Pt2;
-        Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() * 
+        G4double Pt2(0.);
+        if(AveragePt2 <= 0.) {Pt2=0.;}
+        else
+        {
+         Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() * 
                 (std::exp(-maxPtSquare/AveragePt2)-1.)); 
-        
+        }
         G4double Pt=std::sqrt(Pt2);
         G4double phi=G4UniformRand() * twopi;

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

@@ -1 @@
-//********************************************************************
+//
+// ********************************************************************
 // * License and Disclaimer                                           *
 // *                                                                  *
@@ -24 +25 @@
 //
 //
-// $Id: G4FTFParameters.cc,v 1.11 2009/10/25 10:50:54 vuzhinsk Exp $
-// GEANT4 tag $Name: geant4-09-03-cand-01 $
+// $Id: G4FTFParameters.cc,v 1.13 2009/12/16 17:51:15 gunter Exp $
+// GEANT4 tag $Name: geant4-09-03 $
 //
 
@@ -31 +32 @@
 
 #include "G4ios.hh"
-#include <utility>                                        // Uzhi 29.03.08
+#include <utility>                                        
 
 G4FTFParameters::G4FTFParameters()
@@ -231 +232 @@
               SetProjMinNonDiffMass(0.3);                 // GeV
               SetProbabilityOfProjDiff(0.*0.62*std::pow(s/GeV/GeV,-0.51)); // 40/32 X-dif/X-inel
-//G4cout<<"Params "<<0.6*0.62*std::pow(s/GeV/GeV,-0.51)<<G4endl;
+
               SetTarMinDiffMass(1.1);                     // GeV
               SetTarMinNonDiffMass(1.1);                  // GeV
-//G4cout<<"       "<<2.7*0.62*std::pow(s/GeV/GeV,-0.51)<<G4endl; // was 2
-//G4int Uzhi; G4cin>>Uzhi;
+
               SetProbabilityOfTarDiff(2.*0.62*std::pow(s/GeV/GeV,-0.51)); // 40/32 X-dif/X-inel
 
-/*
-SetProjMinDiffMass(0.5);
-SetProjMinNonDiffMass(0.3);   // Uzhi 12.06.08
-SetProbabilityOfProjDiff(0.05); 
-SetProbabilityOfTarDiff(0.05);
-*/
               SetAveragePt2(0.3);                         // GeV^2
              }

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

@@ -25 +25 @@
 //
 //
-// $Id: G4FTFParticipants.cc,v 1.15 2009/10/25 10:50:54 vuzhinsk Exp $
-// GEANT4 tag $Name: geant4-09-03-cand-01 $
+// $Id: G4FTFParticipants.cc,v 1.16 2009/11/25 09:14:03 vuzhinsk Exp $
+// GEANT4 tag $Name: geant4-09-03 $
 //
 // ------------------------------------------------------------
@@ -72 +72 @@
 
 void G4FTFParticipants::GetList(const G4ReactionProduct  &thePrimary,
-                                      G4FTFParameters    *theParameters) // Uzhi 29.03.08
+                                      G4FTFParameters    *theParameters) 
 {
     
@@ -96 +96 @@
         G4double impactY = theImpactParameter.second;
 
-        G4ThreeVector thePosition(impactX, impactY, -DBL_MAX);     //Uzhi 24 July 09 
-        primarySplitable->SetPosition(thePosition);                //Uzhi 24 July 09
+        G4ThreeVector thePosition(impactX, impactY, -DBL_MAX);     
+        primarySplitable->SetPosition(thePosition);                
 
         theNucleus->StartLoop();
         G4Nucleon * nucleon;
-//G4int InterNumber=0;           // Uzhi
-G4int NucleonNumber(0);        // Uzhi
-//while ( (nucleon=theNucleus->GetNextNucleon())&& (InterNumber < 1) ) // Uzhi
-        while ( (nucleon=theNucleus->GetNextNucleon()) ) // Uzhi
+
+        while ( (nucleon=theNucleus->GetNextNucleon()) ) 
         {
-//G4cout<<"Nucl# "<<NucleonNumber<<G4endl; // Vova
            G4double impact2= sqr(impactX - nucleon->GetPosition().x()) +
                              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 
+           if ( theParameters->GetProbabilityOfInteraction(impact2/fermi/fermi) 
                 > G4UniformRand() )
            {
-//InterNumber++;
-/*                                                    // Uzhi 20 July 2009
-                if ( nucleusNeedsShift ) 
-                {                       // on the first hit, shift nucleus 
-                     nucleusNeedsShift = false;
-                     theNucleus->DoTranslation(G4ThreeVector(-1*impactX,-1*impactY,0.));
-                     impactX=0;
-                     impactY=0;
-                }
-*/                                                    // Uzhi 20 July 2009
-//G4cout<<"          Interact"<<G4endl;
                 primarySplitable->SetStatus(1);        // It takes part in the interaction
 
@@ -130 +115 @@
                 if ( ! nucleon->AreYouHit() )
                 {
-//G4cout<<"Part "<<nucleon->Get4Momentum()<<G4endl;
                     targetSplitable= new G4DiffractiveSplitableHadron(*nucleon);
                     nucleon->Hit(targetSplitable);
-                    nucleon->SetBindingEnergy(3.*nucleon->GetBindingEnergy()); // Uzhi 5.10.09
+                    nucleon->SetBindingEnergy(3.*nucleon->GetBindingEnergy()); 
                     targetSplitable->SetStatus(1);     // It takes part in the interaction
                 }
@@ -142 +126 @@
                 theInteractions.push_back(aInteraction);
            }
-NucleonNumber++; // Uzhi
         }    
-// // Uzhi
+
 //      G4cout << "Number of Hit nucleons " << theInteractions.size()<<G4endl; //  entries() 
 //              << "\t" << impactX/fermi << "\t"<<impactY/fermi
 //              << "\t" << std::sqrt(sqr(impactX)+sqr(impactY))/fermi <<G4endl;
-// // Uzhi
+
     }