Ignore:
Timestamp:
Apr 6, 2009, 12:30:29 PM (15 years ago)
Author:
garnier
Message:

update processes

Location:
trunk/source/processes/hadronic/models/parton_string/diffraction/src
Files:
5 edited

Legend:

Unmodified
Added
Removed

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

    r819 r962  
    2525//
    2626//
    27 // $Id: G4DiffractiveExcitation.cc,v 1.1 2007/05/25 06:56:53 vuzhinsk Exp $
     27// $Id: G4DiffractiveExcitation.cc,v 1.7 2008/12/18 13:01:58 gunter Exp $
    2828// ------------------------------------------------------------
    2929//      GEANT 4 class implemetation file
    … …  
    3131//      ---------------- G4DiffractiveExcitation --------------
    3232//             by Gunter Folger, October 1998.
    33 //         diffractive Excitation used by strings models
    34 //         Take a projectile and a target
    35 //         excite the projectile and target
     33//        diffractive Excitation used by strings models
     34//             Take a projectile and a target
     35//             excite the projectile and target
    3636//  Essential changed by V. Uzhinsky in November - December 2006
    3737//  in order to put it in a correspondence with original FRITIOF
    3838//  model. Variant of FRITIOF with nucleon de-excitation is implemented.
    3939//  Other changes by V.Uzhinsky in May 2007 were introduced to fit
    40 //  meson-nucleon interactions 
     40//  meson-nucleon interactions. Additional changes by V. Uzhinsky
     41//  were introduced in December 2006. They treat diffraction dissociation
     42//  processes more exactly.
    4143// ---------------------------------------------------------------------
    4244
    … …  
    5153#include "G4VSplitableHadron.hh"
    5254#include "G4ExcitedString.hh"
     55#include "G4FTFParameters.hh"                            // Uzhi 19.04.08
    5356//#include "G4ios.hh"
    54 
    55 G4DiffractiveExcitation::G4DiffractiveExcitation()                         // Uzhi
    56 {
    57 }
    58 
     57//#include "UZHI_diffraction.hh"
     58
     59G4DiffractiveExcitation::G4DiffractiveExcitation()
     60{
     61}
     62
     63// ---------------------------------------------------------------------
    5964G4bool G4DiffractiveExcitation::
    60   ExciteParticipants(G4VSplitableHadron *projectile, G4VSplitableHadron *target) const
    61 {
    62 
    63            G4LorentzVector Pprojectile=projectile->Get4Momentum();
    64 
    65 // -------------------- Projectile parameters -----------------------------------
    66            G4bool PutOnMassShell=0;
    67 
    68 // With de-excitation
    69 //         G4double M0projectile=projectile->GetDefinition()->GetPDGMass(); 
    70 // Without de-excitation
    71            G4double M0projectile = Pprojectile.mag();                       
    72 
    73            if(M0projectile < projectile->GetDefinition()->GetPDGMass())
    74              {
    75               PutOnMassShell=1;
    76               M0projectile=projectile->GetDefinition()->GetPDGMass();
    77              }
    78 
    79            G4double Mprojectile2 = M0projectile * M0projectile;
    80 
    81            G4int    PDGcode=projectile->GetDefinition()->GetPDGEncoding();
    82            G4int    absPDGcode=std::abs(PDGcode);
    83            G4double ProjectileDiffCut;
    84            G4double ProjectileBackgroundWeight;           // Uzhi May 2007
    85 
    86            G4double TargetBackgroundWeight;               // Uzhi May 2007
    87 
    88            G4double AveragePt2;
    89 
    90            if( absPDGcode > 1000 )                        //------Projectile is baryon --------
    91              {
    92               ProjectileDiffCut = 1.1;                    // GeV
    93               ProjectileBackgroundWeight=0.;              // Uzhi May 2007
    94               AveragePt2 = 0.3;                           // GeV^2
    95              }
    96            else if( absPDGcode == 211 || PDGcode ==  111) //------Projectile is Pion -----------
    97              {
    98               ProjectileDiffCut = 0.6;                   // GeV Uzhi May 2007 1.0->0.6
    99               ProjectileBackgroundWeight=0.9;            // Uzhi May 2007
    100               AveragePt2 = 0.3;                          // GeV^2
    101              }
    102            else if( absPDGcode == 321 || PDGcode == -311) //------Projectile is Kaon -----------
    103              {
    104               ProjectileDiffCut = 0.7;                    // GeV 1.1
    105               ProjectileBackgroundWeight=0.5;             // Uzhi May 2007   ???
    106               AveragePt2 = 0.3;                           // GeV^2
    107              }
    108            else                                           //------Projectile is undefined,
    109                                                           //------Nucleon assumed
    110              {
    111               ProjectileDiffCut = 1.1;                    // GeV
    112               ProjectileBackgroundWeight=0.;              // Uzhi May 2007
    113               AveragePt2 = 0.3;                           // GeV^2
    114              };
    115 
    116            ProjectileDiffCut = ProjectileDiffCut * GeV;
    117            AveragePt2 = AveragePt2 * GeV*GeV;
    118 
    119 // -------------------- Target parameters ----------------------------------------------
    120            G4LorentzVector Ptarget=target->Get4Momentum();
    121 
    122            G4double M0target = Ptarget.mag();
    123 
    124            if(M0target < target->GetDefinition()->GetPDGMass())
    125              {
    126               PutOnMassShell=1;
    127               M0target=target->GetDefinition()->GetPDGMass();
    128              }
    129      
    130            G4double Mtarget2 = M0target * M0target;                      //Ptarget.mag2();
    131                                                                          // for AA-inter.
    132 
    133            G4double NuclearNucleonDiffCut = 1.1*GeV;       
    134            TargetBackgroundWeight=0.;                           // Uzhi May 2007
    135 
    136            G4double ProjectileDiffCut2 = ProjectileDiffCut * ProjectileDiffCut;
    137            G4double NuclearNucleonDiffCut2 = NuclearNucleonDiffCut * NuclearNucleonDiffCut;
     65  ExciteParticipants(G4VSplitableHadron *projectile,
     66                     G4VSplitableHadron *target,
     67                     G4FTFParameters    *theParameters) const
     68{
     69     G4bool PutOnMassShell=0;
     70
     71// -------------------- Projectile parameters -----------------------
     72
     73     G4LorentzVector Pprojectile=projectile->Get4Momentum();
     74//   G4double M0projectile=projectile->GetDefinition()->GetPDGMass(); // With de-excitation
     75     G4double M0projectile = Pprojectile.mag();                       // Without de-excitation
     76/*
     77G4cout<<"ExciteParticipants-------------------"<<G4endl;
     78G4cout<<"Mom "<<Pprojectile<<" mass "<<M0projectile<<G4endl;
     79*/
     80     if(M0projectile < projectile->GetDefinition()->GetPDGMass())
     81     {
     82      PutOnMassShell=1;
     83      M0projectile=projectile->GetDefinition()->GetPDGMass();
     84     }
     85
     86     G4double M0projectile2 = M0projectile * M0projectile;
     87
     88     G4int    PDGcode=projectile->GetDefinition()->GetPDGEncoding();
     89     G4int    absPDGcode=std::abs(PDGcode);
     90
     91     G4double ProjectileDiffStateMinMass=theParameters->GetProjMinDiffMass();
     92     G4double ProjectileNonDiffStateMinMass=theParameters->GetProjMinNonDiffMass();
     93     G4double ProbProjectileDiffraction=theParameters->GetProbabilityOfProjDiff();
     94/*
     95G4cout<<ProjectileDiffStateMinMass<<" "<<ProjectileNonDiffStateMinMass<<" "<<ProbProjectileDiffraction<<G4endl;
     96*/
     97// -------------------- Target paraExciteParticipantsmeters -------------------------
     98     G4LorentzVector Ptarget=target->Get4Momentum();
     99     G4double M0target = Ptarget.mag();
     100
     101//G4cout<<"Mom "<<Ptarget<<" mass "<<M0target<<G4endl;
     102
     103     if(M0target < target->GetDefinition()->GetPDGMass())
     104     {
     105      PutOnMassShell=1;
     106      M0target=target->GetDefinition()->GetPDGMass();
     107     }
     108
     109     G4double M0target2 = M0target * M0target;             //Ptarget.mag2();
     110                                                           // for AA-inter.
     111     G4double TargetDiffStateMinMass=theParameters->GetTarMinDiffMass();   
     112     G4double TargetNonDiffStateMinMass=theParameters->GetTarMinNonDiffMass();   
     113     G4double ProbTargetDiffraction=theParameters->GetProbabilityOfTarDiff();
     114/*
     115G4cout<<TargetDiffStateMinMass<<" "<<TargetNonDiffStateMinMass<<" "<<ProbTargetDiffraction<<G4endl;
     116*/
     117     G4double AveragePt2=theParameters->GetAveragePt2();
     118
     119// Kinematical properties of the interactions --------------
     120     G4LorentzVector Psum;      // 4-momentum in CMS
     121     Psum=Pprojectile+Ptarget;
     122     G4double S=Psum.mag2();
     123
     124//G4cout<<" sqrt(s) "<<std::sqrt(S)<<G4endl;
     125
     126// ------------------------------------------------------------------
     127
     128//ProbProjectileDiffraction=1.;
     129//ProbTargetDiffraction    =1.;
     130     G4double ProbOfDiffraction=ProbProjectileDiffraction +
     131                                ProbTargetDiffraction;
     132
     133     if(ProbOfDiffraction!=0.)
     134     {
     135      ProbProjectileDiffraction/=ProbOfDiffraction;
     136     }
     137     else
     138     {
     139      ProbProjectileDiffraction=0.;
     140     }
     141//   ProbTargetDiffraction    /=ProbOfDiffraction;
     142
     143//G4cout<<"ProbOfDiffraction "<<ProbOfDiffraction<<"ProbProjectileDiffraction "<<ProbProjectileDiffraction<<G4endl;   // Vova
     144
     145     G4double ProjectileDiffStateMinMass2    = ProjectileDiffStateMinMass    *
     146                                               ProjectileDiffStateMinMass;
     147     G4double ProjectileNonDiffStateMinMass2 = ProjectileNonDiffStateMinMass *
     148                                               ProjectileNonDiffStateMinMass;
     149
     150     G4double TargetDiffStateMinMass2        = TargetDiffStateMinMass        *
     151                                               TargetDiffStateMinMass;
     152     G4double TargetNonDiffStateMinMass2     = TargetNonDiffStateMinMass     *
     153                                               TargetNonDiffStateMinMass;
    138154
    139155// Transform momenta to cms and then rotate parallel to z axis;
    140156
    141            G4LorentzVector Psum;
    142            Psum=Pprojectile+Ptarget;
    143 
    144            G4LorentzRotation toCms(-1*Psum.boostVector());
    145 
    146            G4LorentzVector Ptmp=toCms*Pprojectile;
    147 
    148            if ( Ptmp.pz() <= 0. )
    149            {
     157//         G4LorentzVector Psum;
     158//         Psum=Pprojectile+Ptarget;
     159
     160     G4LorentzRotation toCms(-1*Psum.boostVector());
     161
     162     G4LorentzVector Ptmp=toCms*Pprojectile;
     163     if ( Ptmp.pz() <= 0. )
     164        {
    150165           // "String" moving backwards in  CMS, abort collision !!
    151166           //G4cout << " abort Collision!! " << G4endl;
    152                    return false;
    153            }
    154                            
    155            toCms.rotateZ(-1*Ptmp.phi());
    156            toCms.rotateY(-1*Ptmp.theta());
    157        
    158            G4LorentzRotation toLab(toCms.inverse());
    159 
    160            Pprojectile.transform(toCms);
    161            Ptarget.transform(toCms);
    162 
    163            G4double Pt2;                                                   
    164            G4double ProjMassT2, ProjMassT;                                 
    165            G4double TargMassT2, TargMassT;                                 
    166            G4double PZcms2, PZcms;                                         
    167            G4double PMinusNew, TPlusNew;                                   
    168 
    169            G4double S=Psum.mag2();                                         
    170            G4double SqrtS=std::sqrt(S);                                     
    171 
    172            if(absPDGcode > 1000 && SqrtS < 2200*MeV)
    173              {return false;}                         // The model cannot work for
    174                                                      // p+p-interactions
    175                                                      // at Plab < 1.3 GeV/c.
    176 
    177            if(( absPDGcode == 211 || PDGcode ==  111) && SqrtS < 1600*MeV)
    178              {return false;}                         // The model cannot work for
    179                                                      // Pi+p-interactions
    180                                                      // at Plab < 1. GeV/c. 
    181 
    182            if(( absPDGcode == 321 || PDGcode == -311) && SqrtS < 1600*MeV)
    183              {return false;}                         // The model cannot work for
    184                                                      // K+p-interactions
    185                                                      // at Plab < ??? GeV/c.  ???
    186 
    187            PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2-
    188                                  2*S*Mprojectile2-2*S*Mtarget2-2*Mprojectile2*Mtarget2)/4./S;
    189            if(PZcms2 < 0)
    190              {return false;}   // It can be in an interaction with off-shell nuclear nucleon
    191 
    192            PZcms = std::sqrt(PZcms2);
    193 
    194            if(PutOnMassShell)
     167          return false;
     168         }
     169
     170     toCms.rotateZ(-1*Ptmp.phi());
     171     toCms.rotateY(-1*Ptmp.theta());
     172
     173     G4LorentzRotation toLab(toCms.inverse());
     174
     175     Pprojectile.transform(toCms);
     176     Ptarget.transform(toCms);
     177
     178     G4double Pt2;
     179     G4double ProjMassT2, ProjMassT;
     180     G4double TargMassT2, TargMassT;
     181     G4double PZcms2, PZcms;
     182     G4double PMinusMin, PMinusMax;
     183//     G4double PPlusMin , PPlusMax;
     184     G4double TPlusMin , TPlusMax;
     185     G4double PMinusNew, PPlusNew, TPlusNew, TMinusNew;
     186
     187//   G4double S=Psum.mag2();
     188     G4double SqrtS=std::sqrt(S);
     189
     190     if(absPDGcode > 1000 && SqrtS < 2200*MeV)
     191     {return false;}                         // The model cannot work for
     192                                             // p+p-interactions
     193                                             // at Plab < 1.3 GeV/c.
     194
     195     if(( absPDGcode == 211 || PDGcode ==  111) && SqrtS < 1600*MeV)
     196     {return false;}                         // The model cannot work for
     197                                             // Pi+p-interactions
     198                                             // at Plab < 1. GeV/c.
     199
     200     if(( absPDGcode == 321 || PDGcode == -311) && SqrtS < 1600*MeV)
     201     {return false;}                         // The model cannot work for
     202                                             // K+p-interactions
     203                                             // at Plab < ??? GeV/c.  ???
     204
     205     PZcms2=(S*S+M0projectile2*M0projectile2+M0target2*M0target2-
     206             2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2)
     207             /4./S;
     208
     209     if(PZcms2 < 0)
     210     {return false;}   // It can be in an interaction with off-shell nuclear nucleon
     211
     212     PZcms = std::sqrt(PZcms2);
     213
     214     if(PutOnMassShell)
     215     {
     216      if(Pprojectile.z() > 0.)
     217      {
     218       Pprojectile.setPz( PZcms);
     219       Ptarget.setPz(    -PZcms);
     220      }
     221      else
     222      {
     223       Pprojectile.setPz(-PZcms);
     224       Ptarget.setPz(     PZcms);
     225      };
     226
     227      Pprojectile.setE(std::sqrt(M0projectile2                  +
     228                                 Pprojectile.x()*Pprojectile.x()+
     229                                 Pprojectile.y()*Pprojectile.y()+
     230                                 PZcms2));
     231      Ptarget.setE(std::sqrt(M0target2              +
     232                             Ptarget.x()*Ptarget.x()+
     233                             Ptarget.y()*Ptarget.y()+
     234                             PZcms2));
     235     }
     236
     237     G4double maxPtSquare; // = PZcms2;
     238/*
     239G4cout << "Pprojectile aft boost : " << Pprojectile <<" "<<Pprojectile.mag()<< G4endl;
     240G4cout << "Ptarget aft boost : " << Ptarget <<" "<<Ptarget.mag()<< G4endl;
     241G4cout << "cms aft boost : " << (Pprojectile+ Ptarget) << G4endl;
     242G4cout << " Projectile Xplus / Xminus : " <<
     243            Pprojectile.plus() << " / " << Pprojectile.minus() << G4endl;
     244G4cout << " Target Xplus / Xminus : " <<           Ptarget.plus() << " / " << Ptarget.minus() << G4endl;
     245G4cout<<"maxPtSquare "<<maxPtSquare<<G4endl;
     246*/
     247     G4LorentzVector Qmomentum;
     248     G4double Qminus, Qplus;
     249
     250     G4int whilecount=0;
     251//  Choose a process
     252
     253     if(G4UniformRand() < ProbOfDiffraction)
     254       {
     255        if(G4UniformRand() < ProbProjectileDiffraction)
     256        { //-------- projectile diffraction ---------------
     257//G4cout<<"   Projectile diffraction"<<G4endl;
     258//Uzhi_projectilediffraction++;
     259         do {
     260//             Generate pt
     261//             if (whilecount++ >= 500 && (whilecount%100)==0)
     262//               G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping"
     263//               << ", loop count/ maxPtSquare : "
     264//               << whilecount << " / " << maxPtSquare << G4endl;
     265             if (whilecount > 1000 )
    195266             {
    196               if(Pprojectile.z() > 0.)
    197                 {
    198                  Pprojectile.setPz( PZcms);
    199                  Ptarget.setPz(    -PZcms);
    200                 }
    201               else
    202                  {
    203                  Pprojectile.setPz(-PZcms);
    204                  Ptarget.setPz(     PZcms);
    205                 };
    206 
    207                Pprojectile.setE(std::sqrt(Mprojectile2+
    208                                                        Pprojectile.x()*Pprojectile.x()+
    209                                                        Pprojectile.y()*Pprojectile.y()+
    210                                                        PZcms2));
    211                Ptarget.setE(std::sqrt(    Mtarget2    +
    212                                                        Ptarget.x()*Ptarget.x()+
    213                                                        Ptarget.y()*Ptarget.y()+
    214                                                        PZcms2));
    215              }
    216 
    217            G4double maxPtSquare = PZcms2;
    218 
    219 //G4cout << "Pprojectile aft boost : " << Pprojectile << G4endl;
    220 //G4cout << "Ptarget aft boost : " << Ptarget << G4endl;
    221 //         G4cout << "cms aft boost : " << (Pprojectile+ Ptarget) << G4endl;
    222 
    223 //         G4cout << " Projectile Xplus / Xminus : " <<
    224 //              Pprojectile.plus() << " / " << Pprojectile.minus() << G4endl;
    225 //         G4cout << " Target Xplus / Xminus : " <<
    226 //              Ptarget.plus() << " / " << Ptarget.minus() << G4endl;
    227 
    228            G4LorentzVector Qmomentum;
    229            G4double Qminus, Qplus;                                         
    230 
    231            G4int whilecount=0;
    232            do {
    233 //  Generate pt         
    234 
    235                if (whilecount++ >= 500 && (whilecount%100)==0)
     267              Qmomentum=G4LorentzVector(0.,0.,0.,0.);
     268              return false;    //  Ignore this interaction
     269             };
     270// --------------- Check that the interaction is possible -----------
     271             ProjMassT2=ProjectileDiffStateMinMass2;
     272             ProjMassT =ProjectileDiffStateMinMass;
     273
     274             TargMassT2=M0target2;
     275             TargMassT =M0target;
     276
     277             PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
     278                     2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
     279                    /4./S;
     280//G4cout<<" Pt2 Mpt Mtt Pz2 "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl;
     281
     282if(PZcms2 < 0 )
     283{
     284/*
     285G4cout<<"whilecount "<<whilecount<<" "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl;
     286G4int Uzhi; G4cin>>Uzhi;
     287*/
     288return false;
     289};
     290             maxPtSquare=PZcms2;
     291
     292             Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
     293             Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
     294
     295             ProjMassT2=ProjectileDiffStateMinMass2+Pt2;
     296             ProjMassT =std::sqrt(ProjMassT2);
     297
     298             TargMassT2=M0target2+Pt2;
     299             TargMassT =std::sqrt(TargMassT2);
     300
     301             PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
     302                     2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
     303                    /4./S;
     304//G4cout<<" Pt2 Mpt Mtt Pz2 "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl;
     305
     306//             if(PZcms2 < 0 ) {PZcms2=0;};
     307             if(PZcms2 < 0 ) continue;
     308             PZcms =std::sqrt(PZcms2);
     309
     310             PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
     311             PMinusMax=SqrtS-TargMassT;
     312//G4cout<<" SqrtS P+mim max "<<SqrtS<<" "<<PMinusMin<<" "<<PMinusMax<<G4endl;
     313
     314             PMinusNew=ChooseP(PMinusMin, PMinusMax);
     315// PMinusNew=1./sqrt(1./PMinusMin-G4UniformRand()*(1./PMinusMin-1./PMinusMax));
     316
     317             TMinusNew=SqrtS-PMinusNew;
     318             Qminus=Ptarget.minus()-TMinusNew;
     319             TPlusNew=TargMassT2/TMinusNew;
     320             Qplus=Ptarget.plus()-TPlusNew;
     321
     322             Qmomentum.setPz( (Qplus-Qminus)/2 );
     323             Qmomentum.setE(  (Qplus+Qminus)/2 );
     324          } while (
     325((Pprojectile+Qmomentum).mag2() <  ProjectileDiffStateMinMass2) ||  //No without excitation
     326((Ptarget    -Qmomentum).mag2() <  M0target2                  ));
     327        }
     328        else
     329        { // -------------- Target diffraction ----------------
     330//G4cout<<"   Target difraction"<<G4endl;
     331//Uzhi_targetdiffraction++;
     332         do {
     333//             Generate pt
     334//             if (whilecount++ >= 500 && (whilecount%100)==0)
    236335//               G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping"
    237 //               << ", loop count/ maxPtSquare : " 
     336//               << ", loop count/ maxPtSquare : "
    238337//               << whilecount << " / " << maxPtSquare << G4endl;
    239                if (whilecount > 1000 )
    240                {
    241                  Qmomentum=G4LorentzVector(0.,0.,0.,0.);
    242                  return false;    //  Ignore this interaction
    243                }
    244 
    245                Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
    246 
    247 //G4cout << "generated Pt " << Qmomentum << G4endl;
    248 //G4cout << "Pprojectile with pt : " << Pprojectile+Qmomentum << G4endl;
    249 //G4cout << "Ptarget with pt : " << Ptarget-Qmomentum << G4endl;
    250 
    251 //  Momentum transfer
    252 /*                                                                          // Uzhi
    253                G4double Xmin = minmass / ( Pprojectile.e() + Ptarget.e() );
    254                G4double Xmax=1.;
    255                G4double Xplus =ChooseX(Xmin,Xmax);
    256                G4double Xminus=ChooseX(Xmin,Xmax);
    257 
    258 //             G4cout << " X-plus  " << Xplus << G4endl;
    259 //             G4cout << " X-minus " << Xminus << G4endl;
    260                
    261                G4double pt2=G4ThreeVector(Qmomentum.vect()).mag2();
    262                G4double Qplus =-1 * pt2 / Xminus/Ptarget.minus();
    263                G4double Qminus=     pt2 / Xplus /Pprojectile.plus();
    264 */                                                                          // Uzhi    *
    265 
    266                Pt2=G4ThreeVector(Qmomentum.vect()).mag2();                 
    267                ProjMassT2=Mprojectile2+Pt2;                           
    268                ProjMassT =std::sqrt(ProjMassT2);                           
    269 
    270                TargMassT2=Mtarget2+Pt2;                               
    271                TargMassT =std::sqrt(TargMassT2);                           
    272 
    273                PZcms2=(S*S+ProjMassT2*ProjMassT2+                           
    274                            TargMassT2*TargMassT2-                           
    275                            2.*S*ProjMassT2-2.*S*TargMassT2-                 
    276                            2.*ProjMassT2*TargMassT2)/4./S;                 
    277                if(PZcms2 < 0 ) {PZcms2=0;};                                 
    278                PZcms =std::sqrt(PZcms2);                                   
    279 
    280                G4double PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;       
    281                G4double PMinusMax=SqrtS-TargMassT;                         
    282 
    283                if(G4UniformRand() < ProjectileBackgroundWeight)              // Uzhi May 2007
    284                  {
    285                   PMinusNew=PMinusMax-(PMinusMax-PMinusMin)*G4UniformRand(); // Uzhi May 2007
    286                  }                                                           // Uzhi May 2007
    287                else                                                          // Uzhi May 2007
    288                  {                                                           // Uzhi May 2007
    289                   PMinusNew=ChooseP(PMinusMin, PMinusMax);                   // Uzhi May 2007
    290                  };                                                          // Uzhi May 2007
    291 //               PMinusNew=ChooseP(PMinusMin, PMinusMax);                    // Uzhi May 2007
    292                Qminus=PMinusNew-Pprojectile.minus();                       
    293 
    294                G4double TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;       
    295                G4double TPlusMax=SqrtS-PMinusNew;                            // Uzhi May 2007
    296 //               G4double TPlusMax=SqrtS-ProjMassT;                          // Uzhi May 2007
    297 
    298                if(G4UniformRand() < TargetBackgroundWeight)                  // Uzhi May 2007
    299                  {
    300                   TPlusNew=TPlusMax-(TPlusMax-TPlusMin)*G4UniformRand();     // Uzhi May 2007
    301                  }                                                           // Uzhi May 2007
    302                else                                                          // Uzhi May 2007
    303                  {                                                           // Uzhi May 2007
    304                   TPlusNew=ChooseP(TPlusMin, TPlusMax);                      // Uzhi May 2007
    305                  };                                                          // Uzhi May 2007
    306 
    307 //               TPlusNew=ChooseP(TPlusMin, TPlusMax);                       // Uzhi May 2007
    308                Qplus=-(TPlusNew-Ptarget.plus());                           
    309 
    310                Qmomentum.setPz( (Qplus-Qminus)/2 );
    311                Qmomentum.setE(  (Qplus+Qminus)/2 );
    312 
    313 //G4cout << "Qplus / Qminus " << Qplus << " / " << Qminus<<G4endl;
    314 //           G4cout << "pt2" << pt2 << G4endl;
    315 //           G4cout << "Qmomentum " << Qmomentum << G4endl;
    316 //           G4cout << " Masses (P/T) : " << (Pprojectile+Qmomentum).mag() <<
    317 //                             " / " << (Ptarget-Qmomentum).mag() << G4endl;
    318 
    319            } while (
    320 ( (Pprojectile+Qmomentum).mag2() <  Mprojectile2       ||      // Uzhi No without excitation
    321   (Ptarget    -Qmomentum).mag2() <  Mtarget2             ) ||  // Uzhi   
    322 ( (Pprojectile+Qmomentum).mag2() <  ProjectileDiffCut2 &&      // Uzhi No double Diffraction
    323   (Ptarget    -Qmomentum).mag2()  <  NuclearNucleonDiffCut2) );// Uzhi
    324 
    325            if((Ptarget-Qmomentum).mag2() < NuclearNucleonDiffCut2)           // Uzhi
    326 //Projectile diffraction
     338             if (whilecount > 1000 )
    327339             {
    328               G4double TMinusNew=SqrtS-PMinusNew;
    329               Qminus=Ptarget.minus()-TMinusNew;
    330               TPlusNew=TargMassT2/TMinusNew;
    331               Qplus=Ptarget.plus()-TPlusNew;
    332 
    333               Qmomentum.setPz( (Qplus-Qminus)/2 );                         
    334               Qmomentum.setE(  (Qplus+Qminus)/2 );                         
    335              }
    336            else if((Pprojectile+Qmomentum).mag2() <  ProjectileDiffCut2)    // Uzhi
    337 //Target diffraction
     340              Qmomentum=G4LorentzVector(0.,0.,0.,0.);
     341              return false;    //  Ignore this interaction
     342             };
     343// --------------- Check that the interaction is possible -----------
     344             ProjMassT2=M0projectile2;
     345             ProjMassT =M0projectile;
     346
     347             TargMassT2=TargetDiffStateMinMass2;
     348             TargMassT =TargetDiffStateMinMass;
     349
     350             PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
     351                     2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
     352                    /4./S;
     353//G4cout<<" Pt2 Mpt Mtt Pz2 "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl;
     354
     355if(PZcms2 < 0 )
     356{
     357/*
     358G4cout<<"whilecount "<<whilecount<<" "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl;
     359G4int Uzhi; G4cin>>Uzhi;
     360*/
     361return false;
     362};
     363             maxPtSquare=PZcms2;
     364
     365             Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
     366             Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
     367
     368             ProjMassT2=M0projectile2+Pt2;
     369             ProjMassT =std::sqrt(ProjMassT2);
     370
     371             TargMassT2=TargetDiffStateMinMass2+Pt2;
     372             TargMassT =std::sqrt(TargMassT2);
     373
     374             PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
     375                     2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
     376                    /4./S;
     377/*
     378if(PZcms2 < 0 )
     379{
     380G4cout<<"whilecount "<<whilecount<<" "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl;
     381G4int Uzhi; G4cin>>Uzhi;
     382return false;
     383};
     384*/
     385             if(PZcms2 < 0 ) continue;
     386             PZcms =std::sqrt(PZcms2);
     387
     388             TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
     389             TPlusMax=SqrtS-ProjMassT;
     390
     391//G4cout<<" Tmin max "<<TPlusMin<<" "<<TPlusMax<<G4endl;
     392
     393             TPlusNew=ChooseP(TPlusMin, TPlusMax);
     394
     395//TPlusNew=TPlusMax;
     396//G4cout<<"T+new "<<TPlusNew<<G4endl;
     397
     398             PPlusNew=SqrtS-TPlusNew;
     399             Qplus=PPlusNew-Pprojectile.plus();
     400             PMinusNew=ProjMassT2/PPlusNew;
     401             Qminus=PMinusNew-Pprojectile.minus();
     402
     403             Qmomentum.setPz( (Qplus-Qminus)/2 );
     404             Qmomentum.setE(  (Qplus+Qminus)/2 );
     405
     406          } while (
     407 ((Pprojectile+Qmomentum).mag2() <  M0projectile2          ) ||  //No without excitation
     408 ((Ptarget    -Qmomentum).mag2() <  TargetDiffStateMinMass2));
     409         }
     410        }
     411        else  //----------- Non-diffraction process ------------
     412        {
     413//G4cout<<"   Non-difraction"<<G4endl;
     414         do {
     415//             Generate pt
     416//             if (whilecount++ >= 500 && (whilecount%100)==0)
     417//               G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping"
     418//               << ", loop count/ maxPtSquare : "
     419//               << whilecount << " / " << maxPtSquare << G4endl;
     420             if (whilecount > 1000 )
    338421             {
    339               G4double PPlusNew=SqrtS-TPlusNew;
    340               Qplus=PPlusNew-Pprojectile.plus();
    341               PMinusNew=ProjMassT2/PPlusNew;
    342               Qminus=PMinusNew-Pprojectile.minus();
    343 
    344               Qmomentum.setPz( (Qplus-Qminus)/2 );                         
    345               Qmomentum.setE(  (Qplus+Qminus)/2 );                         
     422              Qmomentum=G4LorentzVector(0.,0.,0.,0.);
     423              return false;    //  Ignore this interaction
    346424             };
     425// --------------- Check that the interaction is possible -----------
     426             ProjMassT2=ProjectileNonDiffStateMinMass2;
     427             ProjMassT =ProjectileNonDiffStateMinMass;
     428
     429             TargMassT2=TargetNonDiffStateMinMass2;
     430             TargMassT =TargetNonDiffStateMinMass;
     431
     432             PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
     433                    2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
     434                   /4./S;
     435//G4cout<<" Pt2 Mpt Mtt Pz2 "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl;
     436
     437if(PZcms2 < 0 )
     438{
     439/*
     440G4cout<<"whilecount "<<whilecount<<" "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl;
     441G4int Uzhi; G4cin>>Uzhi;
     442*/
     443return false;
     444};
     445             maxPtSquare=PZcms2;
     446             Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
     447             Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
     448
     449             ProjMassT2=ProjectileNonDiffStateMinMass2+Pt2;
     450             ProjMassT =std::sqrt(ProjMassT2);
     451
     452             TargMassT2=TargetNonDiffStateMinMass2+Pt2;
     453             TargMassT =std::sqrt(TargMassT2);
     454
     455             PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
     456                    2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
     457                   /4./S;
     458/*
     459G4cout<<"ProjectileNonDiffStateMinMass2 "<<ProjectileNonDiffStateMinMass2<<G4endl;
     460G4cout<<"TargetNonDiffStateMinMass2     "<<TargetNonDiffStateMinMass2<<G4endl;
     461G4cout<<"Mt "<<ProjMassT<<" "<<TargMassT<<" "<<Pt2<<" "<<PZcms2<<G4endl<<G4endl;
     462*/
     463//             if(PZcms2 < 0 ) {PZcms2=0;};
     464             if(PZcms2 < 0 ) continue;
     465             PZcms =std::sqrt(PZcms2);
     466
     467             PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
     468             PMinusMax=SqrtS-TargMassT;
     469
     470             PMinusNew=ChooseP(PMinusMin, PMinusMax);
     471// PMinusNew=1./sqrt(1./PMinusMin-G4UniformRand()*(1./PMinusMin-1./PMinusMax));
     472
     473//G4cout<<"Proj "<<PMinusMin<<" "<<PMinusMax<<" "<<PMinusNew<<G4endl;
     474
     475//PMinusNew=PMinusMax; //+++++++++++++++++++++++++++++++++++ Vova
     476
     477             Qminus=PMinusNew-Pprojectile.minus();
     478
     479             TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
     480//             TPlusMax=SqrtS-PMinusNew;                      // Vova
     481             TPlusMax=SqrtS-ProjMassT;                        // Vova
     482
     483             TPlusNew=ChooseP(TPlusMin, TPlusMax);
     484
     485//G4cout<<"Targ "<<TPlusMin<<" "<<TPlusMax<<" "<<TPlusNew<<G4endl;
     486//G4cout<<PMinusNew<<" "<<TPlusNew<<G4endl;
     487
     488             Qplus=-(TPlusNew-Ptarget.plus());
     489
     490             Qmomentum.setPz( (Qplus-Qminus)/2 );
     491             Qmomentum.setE(  (Qplus+Qminus)/2 );
     492/*
     493G4cout << "Qplus / Qminus " << Qplus << " / " << Qminus<<G4endl;
     494G4cout << "pt2" << pt2 << G4endl;
     495G4cout << "Qmomentum " << Qmomentum << G4endl;
     496G4cout << " Masses (P/T) : " << (Pprojectile+Qmomentum).mag() <<
     497           " / " << (Ptarget-Qmomentum).mag() << G4endl;   // mag()
     498G4cout<<"Mprojectile "<<std::sqrt(M0projectile2)<<G4endl;
     499G4cout<<"Mtarget     "<<std::sqrt(M0target2    )<<G4endl;
     500G4cout<<"ProjectileDiffStateMinMass "<<std::sqrt(ProjectileDiffStateMinMass2)<<G4endl;
     501G4cout<<"TargetDiffStateMinMass "<<std::sqrt(TargetDiffStateMinMass2)<<G4endl;
     502*/
     503       } while (
     504 ((Pprojectile+Qmomentum).mag2() <  ProjectileNonDiffStateMinMass2) || //No double Diffraction
     505 ((Ptarget    -Qmomentum).mag2() <  TargetNonDiffStateMinMass2    ));
     506         }
     507
     508//G4int Uzhiinp; G4cin>>Uzhiinp;    // Vova
    347509
    348510           Pprojectile += Qmomentum;
    349511           Ptarget     -= Qmomentum;
    350 
    351 /*
    352 Pprojectile.setPz(0.);
    353 Pprojectile.setE(SqrtS-M0target);
    354 
    355 Ptarget.setPz(0.);
    356 Ptarget.setE(M0target);
    357 */
    358 
    359 //G4cout << "Pprojectile with Q : " << Pprojectile << G4endl;
    360 //G4cout << "Ptarget with Q : " << Ptarget << G4endl;
    361        
    362 //         G4cout << "Projectile back: " << toLab * Pprojectile << G4endl;
    363 //         G4cout << "Target back: " << toLab * Ptarget << G4endl;
     512/*
     513G4cout << "Pprojectile with Q : " << Pprojectile << G4endl;
     514G4cout << "Ptarget with Q : " << Ptarget << G4endl;
     515G4cout << "Target        mass  " <<  Ptarget.mag() << G4endl;
     516G4cout << "Projectile mass  " <<  Pprojectile.mag() << G4endl;
     517//
     518//G4cout << "Projectile back: " << toLab * Pprojectile << G4endl;
     519//G4cout << "Target back: " << toLab * Ptarget << G4endl;
     520*/
     521//-------------- Flip if projectale moves in backward direction ------------
     522//G4bool Flip=Pprojectile.pz()< 0.;
     523
    364524
    365525// Transform back and update SplitableHadron Participant.
    … …  
    369529//G4cout << "Pprojectile with Q M: " << Pprojectile<<" "<<  Pprojectile.mag() << G4endl;
    370530//G4cout << "Ptarget     with Q M: " << Ptarget    <<" "<<  Ptarget.mag()     << G4endl;
    371 
    372 //G4cout << "Target      mass  " <<  Ptarget.mag() << G4endl;     
    373                            
     531//G4cout << "Target      mass  " <<  Ptarget.mag() << G4endl;
     532//G4cout << "Projectile mass  " <<  Pprojectile.mag() << G4endl;
     533
     534/*
     535if(!Flip){
     536           projectile->Set4Momentum(Pprojectile);
    374537           target->Set4Momentum(Ptarget);
    375 
    376 //G4cout << "Projectile mass  " <<  Pprojectile.mag() << G4endl;
     538         }
     539else     {
     540           G4ParticleDefinition * t_Definition=projectile->GetDefinition();
     541           projectile->SetDefinition(target->GetDefinition());
     542           projectile->Set4Momentum(Ptarget);
     543           target->SetDefinition(t_Definition);
     544           target->Set4Momentum(Pprojectile);
     545          }
     546*/
     547//
     548/*
     549if(G4UniformRand() < 1.) {
     550           G4ParticleDefinition * t_Definition=projectile->GetDefinition();
     551           projectile->SetDefinition(target->GetDefinition());
     552           target->SetDefinition(t_Definition);
     553}
     554*/ // For flip, for HARP
     555
     556           G4double ZcoordinateOfCurrentInteraction = target->GetPosition().z();
     557// It is assumed that nucleon z-coordinates are ordered on increasing -----------
     558
     559           G4double betta_z=projectile->Get4Momentum().pz()/projectile->Get4Momentum().e();
     560
     561           G4double ZcoordinateOfPreviousCollision=projectile->GetPosition().z();
     562           if(projectile->GetSoftCollisionCount()==0) {
     563              projectile->SetTimeOfCreation(0.);
     564              target->SetTimeOfCreation(0.);
     565              ZcoordinateOfPreviousCollision=ZcoordinateOfCurrentInteraction;
     566           }
     567           
     568           G4ThreeVector thePosition(projectile->GetPosition().x(),
     569                                     projectile->GetPosition().y(),
     570                                     ZcoordinateOfCurrentInteraction);
     571           projectile->SetPosition(thePosition);
     572
     573           G4double TimeOfPreviousCollision=projectile->GetTimeOfCreation();
     574           G4double TimeOfCurrentCollision=TimeOfPreviousCollision+
     575                    (ZcoordinateOfCurrentInteraction-ZcoordinateOfPreviousCollision)/betta_z;
     576
     577           projectile->SetTimeOfCreation(TimeOfCurrentCollision);
     578           target->SetTimeOfCreation(TimeOfCurrentCollision);
    377579
    378580           projectile->Set4Momentum(Pprojectile);
     581           target->Set4Momentum(Ptarget);
     582
     583           projectile->IncrementCollisionCount(1);
     584           target->IncrementCollisionCount(1);
     585
     586//
     587//G4cout<<"Out of Excitation --------------------"<<G4endl;
     588//G4int Uzhiinp; G4cin>>Uzhiinp;    // Vova
    379589
    380590           return true;
    381591}
    382592
    383 
     593// ---------------------------------------------------------------------
    384594G4ExcitedString * G4DiffractiveExcitation::
    385595           String(G4VSplitableHadron * hadron, G4bool isProjectile) const
    386596{
     597
     598//G4cout<<"G4DiffractiveExcitation::String isProj"<<isProjectile<<G4endl;
     599
    387600        hadron->SplitUp();
    388601        G4Parton *start= hadron->GetNextParton();
    389         if ( start==NULL) 
     602        if ( start==NULL)
    390603        { G4cout << " G4FTFModel::String() Error:No start parton found"<< G4endl;
    391604          return NULL;
    392605        }
    393606        G4Parton *end  = hadron->GetNextParton();
    394         if ( end==NULL) 
     607        if ( end==NULL)
    395608        { G4cout << " G4FTFModel::String() Error:No end parton found"<< G4endl;
    396609          return NULL;
    397610        }
    398        
     611
    399612        G4ExcitedString * string;
    400         if ( isProjectile ) 
     613        if ( isProjectile )
    401614        {
    402615                string= new G4ExcitedString(end,start, +1);
    … …  
    404617                string= new G4ExcitedString(start,end, -1);
    405618        }
    406        
     619// Uzhi
     620//G4cout<<"G4ExcitedString * G4DiffractiveExcitation::String"<<G4endl;
     621//G4cout<<hadron->GetTimeOfCreation()<<" "<<hadron->GetPosition()/fermi<<G4endl;
     622
     623        string->SetTimeOfCreation(hadron->GetTimeOfCreation());
    407624        string->SetPosition(hadron->GetPosition());
    408625
    409626// momenta of string ends
     627//
     628        G4double Momentum=hadron->Get4Momentum().vect().mag();
     629        G4double  Plus=hadron->Get4Momentum().e() + Momentum;
     630        G4double Minus=hadron->Get4Momentum().e() - Momentum;
     631
     632        G4ThreeVector tmp;
     633        if(Momentum > 0.)
     634        {
     635         tmp.set(hadron->Get4Momentum().px(),
     636                 hadron->Get4Momentum().py(),
     637                 hadron->Get4Momentum().pz());
     638         tmp/=Momentum;
     639         }
     640         else
     641         {
     642          tmp.set(0.,0.,1.);
     643         };
     644
     645        G4LorentzVector Pstart(tmp,0.);
     646        G4LorentzVector   Pend(tmp,0.);
     647
     648        if(isProjectile)
     649        {
     650         Pstart*=(-1.)*Minus/2.;
     651         Pend  *=(+1.)*Plus /2.;
     652         }
     653        else
     654         {
     655          Pstart*=(+1.)*Plus/2.;
     656          Pend  *=(-1.)*Minus/2.;
     657         };
     658
     659         Momentum=-Pstart.mag();
     660         Pstart.setT(Momentum);  // It is assumed that quark has m=0.
     661
     662         Momentum=-Pend.mag();
     663         Pend.setT(Momentum);    // It is assumed that di-quark has m=0.
     664//
     665/* Uzhi
    410666        G4double ptSquared= hadron->Get4Momentum().perp2();
    411667        G4double transverseMassSquared= hadron->Get4Momentum().plus()
    … …  
    416672                 sqr( std::sqrt(transverseMassSquared) - std::sqrt(ptSquared) );
    417673
    418         G4double widthOfPtSquare = 0.25;          // Uzhi <Pt^2>=0.25 ??????????????????
     674        G4double widthOfPtSquare = 0.25*GeV*GeV;       // Uzhi 11.07 <Pt^2>=0.25 ??????????????????
    419675        G4ThreeVector pt=GaussianPt(widthOfPtSquare,maxAvailMomentumSquared);
    420676
    … …  
    432688
    433689        G4int Sign= isProjectile ? -1 : 1;
    434        
     690
    435691        G4double endMinus  = 0.5 * (tm1 + Sign*tm2);
    436692        G4double startMinus= hadron->Get4Momentum().minus() - endMinus;
    … …  
    444700        Pend.setPz(0.5*(endPlus - endMinus));
    445701        Pend.setE(0.5*(endPlus + endMinus));
    446 
     702*/ // Uzhi
    447703        start->Set4Momentum(Pstart);
    448704        end->Set4Momentum(Pend);
    449        
     705/*
     706G4cout<<"G4DiffractiveExcitation::String hadro"<<hadron->Get4Momentum()<<" "<<hadron->Get4Momentum().mag2()<<G4endl;
     707
     708G4cout<<"G4DiffractiveExcitation::String start"<<start->Get4Momentum()<<" "<<start->GetPDGcode()<<G4endl;
     709
     710G4cout<<"G4DiffractiveExcitation::String end  "<<  end->Get4Momentum()<<" "<<  end->GetPDGcode()<<G4endl;
     711G4int Uzhi; G4cin>>Uzhi;
     712*/
    450713#ifdef G4_FTFDEBUG
    451         G4cout << " generated string flavors          " << start->GetPDGcode() << " / " << end->GetPDGcode() << G4endl;
    452         G4cout << " generated string momenta:   quark " << start->Get4Momentum() << "mass : " <<start->Get4Momentum().mag()<< G4endl;
    453         G4cout << " generated string momenta: Diquark " << end ->Get4Momentum() << "mass : " <<end->Get4Momentum().mag()<< G4endl;
     714        G4cout << " generated string flavors          "
     715               << start->GetPDGcode() << " / "
     716               << end->GetPDGcode()                    << G4endl;
     717        G4cout << " generated string momenta:   quark "
     718               << start->Get4Momentum() << "mass : "
     719               <<start->Get4Momentum().mag()           << G4endl;
     720        G4cout << " generated string momenta: Diquark "
     721               << end ->Get4Momentum()
     722               << "mass : " <<end->Get4Momentum().mag()<< G4endl;
    454723        G4cout << " sum of ends                       " << Pstart+Pend << G4endl;
    455724        G4cout << " Original                          " << hadron->Get4Momentum() << G4endl;
    456725#endif
    457        
    458         return string; 
     726
     727        return string;
    459728}
    460729
    … …  
    462731// --------- private methods ----------------------
    463732
     733// ---------------------------------------------------------------------
    464734G4double G4DiffractiveExcitation::ChooseP(G4double Pmin, G4double Pmax) const // Uzhi
    465735{
    466736// choose an x between Xmin and Xmax with P(x) ~ 1/x
    467 
    468737//  to be improved...
    469738
    470739        G4double range=Pmax-Pmin;                                         // Uzhi
    471        
    472         if ( Pmin <= 0. || range <=0. ) 
     740
     741        if ( Pmin <= 0. || range <=0. )
    473742        {
    474743                G4cout << " Pmin, range : " << Pmin << " , " << range << G4endl;
    … …  
    489758}
    490759
    491 G4ThreeVector G4DiffractiveExcitation::GaussianPt(G4double AveragePt2, G4double maxPtSquare) const // Uzhi
     760// ---------------------------------------------------------------------
     761G4ThreeVector G4DiffractiveExcitation::GaussianPt(G4double AveragePt2,
     762                                                  G4double maxPtSquare) const // Uzhi
    492763{            //  @@ this method is used in FTFModel as well. Should go somewhere common!
    493        
     764
    494765        G4double Pt2;
    495766/*                                                                          // Uzhi
    … …  
    499770*/                                                                          // Uzhi
    500771
    501         Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() * (std::exp(-maxPtSquare/AveragePt2)-1.));// Uzhi
    502        
     772        Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() *
     773                           (std::exp(-maxPtSquare/AveragePt2)-1.));// Uzhi
     774
    503775        G4double Pt=std::sqrt(Pt2);
    504776
    505777        G4double phi=G4UniformRand() * twopi;
    506        
    507         return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);   
    508 }
    509 
     778
     779        return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);
     780}
     781
     782// ---------------------------------------------------------------------
    510783G4DiffractiveExcitation::G4DiffractiveExcitation(const G4DiffractiveExcitation &)
    511784{

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

    r819 r962  
    3030#include "G4KineticTrack.hh"
    3131#include "G4DiffractiveSplitableHadron.hh"
     32#include "G4FTFParameters.hh"                            // Uzhi 21.04.08
    3233
    3334G4DiffractiveHHScatterer::G4DiffractiveHHScatterer()
    … …  
    3738{}
    3839
     40// -------------------------------------------------------------------
    3941G4KineticTrackVector * G4DiffractiveHHScatterer::
    4042Scatter(const G4KineticTrack & aTrack, const G4KineticTrack & bTrack)
    … …  
    4446  G4DiffractiveSplitableHadron aHadron(& aTrack);
    4547  G4DiffractiveSplitableHadron bHadron(& bTrack);
    46   if ( ! theExcitation->ExciteParticipants(& aHadron, & bHadron))
     48  theParameters = new G4FTFParameters(aHadron.GetDefinition(), // -------- Uzhi 21.04.08
     49                                          1.,1., 100.);
     50                                          //s);// ------------------------- Uzhi 21.04.08
     51  if ( ! theExcitation->ExciteParticipants(& aHadron,
     52                                           & bHadron,
     53                                           theParameters))     // -------- Uzhi 21.04.08
    4754  {
    4855        return NULL;

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

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

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

    r819 r962  
    2525//
    2626//
    27 // $Id: G4FTFModel.cc,v 1.7 2007/04/24 10:32:59 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     27// $Id: G4FTFModel.cc,v 1.13 2008/12/09 10:40:52 vuzhinsk Exp $
     28// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2929//
    3030
    … …  
    3838
    3939#include "G4FTFModel.hh"
     40#include "G4FTFParameters.hh"                            // Uzhi 29.03.08
    4041#include "G4FTFParticipants.hh"
    4142#include "G4InteractionContent.hh"
    … …  
    4344#include "G4ParticleDefinition.hh"
    4445#include "G4ios.hh"
     46#include <utility>                                        // Uzhi 29.03.08
    4547
    4648// Class G4FTFModel
    4749
    48 G4FTFModel::G4FTFModel():theExcitation(new G4DiffractiveExcitation()) // Uzhi
     50G4FTFModel::G4FTFModel():theExcitation(new G4DiffractiveExcitation()),
     51                         theElastic(new G4ElasticHNScattering()) // Uzhi 29.03.08
    4952{
    5053        G4VPartonStringModel::SetThisPointer(this);
    51 }
    52 
    53 G4FTFModel::G4FTFModel(G4double a, G4double b, G4double c):theExcitation(new G4DiffractiveExcitation())
     54        theParameters=0;                                         // Uzhi 9.12.08
     55}
     56
     57/*
     58G4FTFModel::G4FTFModel(G4double , G4double , G4double ):theExcitation(new // Uzhi 9.12.08 G4DiffractiveExcitation())
    5459{
    5560        G4VPartonStringModel::SetThisPointer(this);
    … …  
    6267        G4VPartonStringModel::SetThisPointer(this);
    6368}
    64 
     69*/
    6570
    6671
    6772G4FTFModel::~G4FTFModel()
    68 {}
     73{
     74   if( theParameters != 0 ) delete theParameters;             // Uzhi 5.12.08
     75// Because FTF model can be called for various particles
     76// theParameters must be erased at the end of each call.
     77// Thus the delete is olso in G4FTFModel::GetStrings() method
     78   if( theExcitation != 0 ) delete theExcitation;             // Uzhi 5.12.08
     79   if( theElastic    != 0 ) delete theElastic;                // Uzhi 5.12.08
     80}
    6981
    7082
    … …  
    8698}
    8799
     100// ------------------------------------------------------------
    88101void G4FTFModel::Init(const G4Nucleus & aNucleus, const G4DynamicParticle & aProjectile)
    89102{
    90         theParticipants.Init(aNucleus.GetN(),aNucleus.GetZ()); // Uzhi N-mass number Z-charge
    91  
    92103        theProjectile = aProjectile; 
    93 }
    94 
     104//G4cout<<"G4FTFModel::Init "<<aNucleus.GetN()<<" "<<aNucleus.GetZ()<<G4endl;
     105        theParticipants.Init(aNucleus.GetN(),aNucleus.GetZ());
     106// Uzhi N-mass number Z-charge ------------------------- Uzhi 29.03.08
     107
     108// --- cms energy
     109
     110        G4double s = sqr( theProjectile.GetMass() ) +
     111                     sqr( G4Proton::Proton()->GetPDGMass() ) +
     112                     2*theProjectile.GetTotalEnergy()*G4Proton::Proton()->GetPDGMass();
     113/*
     114G4cout << " primary Total E (GeV): " << theProjectile.GetTotalEnergy()/GeV << G4endl;
     115G4cout << " primary Mass    (GeV): " << theProjectile.GetMass() /GeV << G4endl;
     116G4cout << "cms std::sqrt(s) (GeV) = " << std::sqrt(s) / GeV << G4endl;
     117*/
     118
     119      if( theParameters != 0 ) delete theParameters;                    // Uzhi 9.12.08
     120      theParameters = new G4FTFParameters(theProjectile.GetDefinition(),
     121                                          aNucleus.GetN(),aNucleus.GetZ(),
     122                                          s);// ------------------------- Uzhi 19.04.08
     123//theParameters->SetProbabilityOfElasticScatt(0.); // To turn on/off (1/0) elastic scattering
     124}
     125
     126// ------------------------------------------------------------
    95127G4ExcitedStringVector * G4FTFModel::GetStrings()
    96128{
    97         theParticipants.BuildInteractions(theProjectile);
    98 
     129//G4cout<<"theParticipants.GetList"<<G4endl;
     130        theParticipants.GetList(theProjectile,theParameters);
     131//G4cout<<"ExciteParticipants()"<<G4endl;
    99132        if (! ExciteParticipants()) return NULL;;
    100 
     133//G4cout<<"theStrings = BuildStrings()"<<G4endl;
    101134        G4ExcitedStringVector * theStrings = BuildStrings();
    102 
     135//G4cout<<"Return to theStrings "<<G4endl;
     136        if( theParameters != 0 )                              // Uzhi 9.12.08
     137        {                                                     // Uzhi 9.12.08
     138          delete theParameters;                               // Uzhi 9.12.08
     139          theParameters=0;                                    // Uzhi 9.12.08
     140        }                                                     // Uzhi 9.12.08
    103141        return theStrings;
    104142}
    105143
     144// ------------------------------------------------------------
    106145struct DeleteVSplitableHadron { void operator()(G4VSplitableHadron * aH){delete aH;} };
    107146
    108 G4ExcitedStringVector * G4FTFModel::BuildStrings()
    109 {       
    110 
    111 // Loop over all collisions; find all primaries, and all target ( targets may
    112 //  be duplicate in the List ( to unique G4VSplitableHadrons)
    113 
    114         G4ExcitedStringVector * strings;
    115         strings = new G4ExcitedStringVector();
    116        
    117         std::vector<G4VSplitableHadron *> primaries;
    118         std::vector<G4VSplitableHadron *> targets;
    119        
    120         theParticipants.StartLoop();    // restart a loop
    121         while ( theParticipants.Next() )
    122         {
    123             const G4InteractionContent & interaction=theParticipants.GetInteraction();
    124                  //  do not allow for duplicates ...
    125             if ( primaries.end() == std::find(primaries.begin(), primaries.end(), interaction.GetProjectile()) )
    126                 primaries.push_back(interaction.GetProjectile());
    127                
    128             if ( targets.end() == std::find(targets.begin(), targets.end(),interaction.GetTarget()) )
    129                 targets.push_back(interaction.GetTarget());
    130         }
    131            
    132        
    133 //      G4cout << "BuildStrings prim/targ " << primaries.entries() << " , " <<
    134 //                                           targets.entries() << G4endl;
    135 
    136 
    137         unsigned int ahadron;
    138         for ( ahadron=0; ahadron < primaries.size() ; ahadron++)
    139         {
    140             G4bool isProjectile=true;
    141             strings->push_back(theExcitation->String(primaries[ahadron], isProjectile));
    142         }
    143         for ( ahadron=0; ahadron < targets.size() ; ahadron++)
    144         {
    145             G4bool isProjectile=false;
    146             strings->push_back(theExcitation->String(targets[ahadron], isProjectile));
    147         }
    148 
    149         std::for_each(primaries.begin(), primaries.end(), DeleteVSplitableHadron());
    150         primaries.clear();
    151         std::for_each(targets.begin(), targets.end(), DeleteVSplitableHadron());
    152         targets.clear();
    153        
    154         return strings;
    155 }
    156 
     147// ------------------------------------------------------------
    157148G4bool G4FTFModel::ExciteParticipants()
    158149{
    159        
     150/*    // Uzhi 29.03.08                     For elastic Scatt.
     151G4cout<<"  In ExciteParticipants() "<<theParticipants.theInteractions.size()<<G4endl;
     152G4cout<<" test Params Tot "<<theParameters->GetTotalCrossSection()<<G4endl;
     153G4cout<<" test Params Ela "<<theParameters->GetElasticCrossSection()<<G4endl;
     154       
     155G4int counter=0;
     156*/   // Uzhi 29.03.08
     157
     158
     159
    160160        while (theParticipants.Next())
    161161        {         
    162162           const G4InteractionContent & collision=theParticipants.GetInteraction();
    163 
    164 //G4cout << " soft colls : " << collision.GetNumberOfSoftCollisions() << G4endl; // Uzhi no match
     163/*
     164counter++;
     165G4cout<<" Inter # "<<counter<<G4endl;
     166*/
    165167           G4VSplitableHadron * projectile=collision.GetProjectile();
    166168           G4VSplitableHadron * target=collision.GetTarget();
    167169
    168            if ( ! theExcitation->ExciteParticipants(projectile, target) )
     170//   // Uzhi 29.03.08
     171           G4bool Successfull;
     172           if(G4UniformRand()< theParameters->GetProbabilityOfElasticScatt())
     173           {
     174//G4cout<<"Elastic"<<G4endl;
     175            Successfull=theElastic->ElasticScattering(projectile, target, theParameters);
     176           }
     177           else
     178           {
     179//G4cout<<"Inelastic"<<G4endl;
     180            Successfull=theExcitation->ExciteParticipants(projectile, target, theParameters);
     181           }
     182//           if(!Successfull)
     183// // Uzhi 29.03.08
     184
     185//         if ( ! theExcitation->ExciteParticipants(projectile, target) )
     186           if(!Successfull)
    169187           {
    170188//           give up, clean up
    … …  
    176194                    const G4InteractionContent & interaction=theParticipants.GetInteraction();
    177195                         //  do not allow for duplicates ...
    178                     if ( primaries.end() == std::find(primaries.begin(), primaries.end(), interaction.GetProjectile()) )
     196                    if ( primaries.end() == std::find(primaries.begin(), primaries.end(),
     197                                                      interaction.GetProjectile()) )
    179198                        primaries.push_back(interaction.GetProjectile());
    180199
    181                     if ( targets.end() == std::find(targets.begin(), targets.end(),interaction.GetTarget()) )
     200                    if ( targets.end()   == std::find(targets.begin(), targets.end(),
     201                                                      interaction.GetTarget()) )
    182202                        targets.push_back(interaction.GetTarget());
    183203                }
    … …  
    188208
    189209                return false;
    190            }
    191 
     210           }  // End of the loop Uzhi
     211        }
     212        return true;
     213}
     214// ------------------------------------------------------------
     215G4ExcitedStringVector * G4FTFModel::BuildStrings()
     216{       
     217// Loop over all collisions; find all primaries, and all target ( targets may
     218//  be duplicate in the List ( to unique G4VSplitableHadrons)
     219
     220        G4ExcitedStringVector * strings;
     221        strings = new G4ExcitedStringVector();
     222       
     223        std::vector<G4VSplitableHadron *> primaries;
     224        std::vector<G4VSplitableHadron *> targets;
     225       
     226        theParticipants.StartLoop();    // restart a loop
     227        while ( theParticipants.Next() )
     228        {
     229            const G4InteractionContent & interaction=theParticipants.GetInteraction();
     230                 //  do not allow for duplicates ...
     231            if ( primaries.end() == std::find(primaries.begin(), primaries.end(),
     232                                                interaction.GetProjectile()) )
     233                primaries.push_back(interaction.GetProjectile());
     234               
     235            if ( targets.end()   == std::find(targets.begin(), targets.end(),
     236                                                interaction.GetTarget()) )
     237                targets.push_back(interaction.GetTarget());
    192238        }
    193         return true;
    194 }
     239           
     240       
     241//      G4cout << "BuildStrings prim/targ " << primaries.size() << " , " <<
     242//                                             targets.size() << G4endl;
     243
     244        unsigned int ahadron;
     245// Only for hA-interactions Uzhi -------------------------------------
     246        for ( ahadron=0; ahadron < primaries.size() ; ahadron++)
     247        {
     248//G4ThreeVector aPosition=primaries[ahadron]->GetPosition();
     249//G4cout<<"Proj Build "<<aPosition<<" "<<primaries[ahadron]->GetTimeOfCreation()<<G4endl;
     250            G4bool isProjectile=true;
     251            strings->push_back(theExcitation->String(primaries[ahadron], isProjectile));
     252        }
     253
     254        for ( ahadron=0; ahadron < targets.size() ; ahadron++)
     255        {
     256//G4ThreeVector aPosition=targets[ahadron]->GetPosition();
     257//G4cout<<"Targ Build "<<aPosition<<" "<<targets[ahadron]->GetTimeOfCreation()<<G4endl;
     258            G4bool isProjectile=false;
     259            strings->push_back(theExcitation->String(targets[ahadron], isProjectile));
     260        }
     261
     262        std::for_each(primaries.begin(), primaries.end(), DeleteVSplitableHadron());
     263        primaries.clear();
     264        std::for_each(targets.begin(), targets.end(), DeleteVSplitableHadron());
     265        targets.clear();
     266       
     267        return strings;
     268}
     269// ------------------------------------------------------------

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

    r819 r962  
    2525//
    2626//
    27 // $Id: G4FTFParticipants.cc,v 1.7 2007/04/24 10:33:00 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     27// $Id: G4FTFParticipants.cc,v 1.9 2008/06/13 12:49:23 vuzhinsk Exp $
     28// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2929//
    3030// ------------------------------------------------------------
    … …  
    3838// ------------------------------------------------------------
    3939
     40#include "G4FTFParameters.hh"                            // Uzhi 29.03.08
    4041#include "G4FTFParticipants.hh"
    4142#include "G4DiffractiveSplitableHadron.hh"
    4243#include "G4VSplitableHadron.hh"
    43 //#include "G4PomeronCrossSection.hh"                      // Uzhi
    44 #include "G4FTFCrossSection.hh"                            // Uzhi
    4544#include "Randomize.hh"
    46 #include <utility>
    47 
     45#include <utility>                                        // Uzhi 29.03.08
    4846
    4947// Class G4FTFParticipants
    50 
    51 
    5248
    5349G4FTFParticipants::G4FTFParticipants()
    … …  
    7571//{}
    7672
    77 void G4FTFParticipants::BuildInteractions(const G4ReactionProduct  &thePrimary)
     73void G4FTFParticipants::GetList(const G4ReactionProduct  &thePrimary,
     74                                      G4FTFParameters    *theParameters) // Uzhi 29.03.08
    7875{
    7976   
    … …  
    8380    theInteractions.clear();
    8481
    85 // --- cms energy
    86 
    87     G4double s = sqr( thePrimary.GetMass() ) +
    88                  sqr( G4Proton::Proton()->GetPDGMass() ) +
    89                  2*thePrimary.GetTotalEnergy()*G4Proton::Proton()->GetPDGMass();
    90 
    91 //    G4cout << " primary Total E (GeV): " << thePrimary.GetTotalEnergy()/GeV << G4endl;
    92 //    G4cout << " primary Mass    (GeV): " << thePrimary.GetMass() /GeV << G4endl;
    93 //    G4cout << "cms std::sqrt(s) (GeV) = " << std::sqrt(s) / GeV << G4endl;
    94 
    95 //    G4PomeronCrossSection theCrossSection(thePrimary.GetDefinition());      // Uzhi
    96       G4FTFCrossSection theCrossSection(thePrimary.GetDefinition(),s);        // Uzhi
    97    
    98 
    99     G4double deltaxy=2 * fermi;
     82    G4double deltaxy=2 * fermi;                       // Extra nuclear radius
    10083
    10184    G4VSplitableHadron * primarySplitable=new G4DiffractiveSplitableHadron(thePrimary);
    10285
    103     G4double xyradius;
    104     xyradius =theNucleus->GetOuterRadius() + deltaxy;
    105 
     86    G4double xyradius;                         
     87    xyradius =theNucleus->GetOuterRadius() + deltaxy; // Impact parameter sampling
     88                                                      // radius
    10689    G4bool nucleusNeedsShift = true;
    10790   
    … …  
    11598        theNucleus->StartLoop();
    11699        G4Nucleon * nucleon;
    117         while ( (nucleon=theNucleus->GetNextNucleon()) )
     100//G4int InterNumber=0;           // Uzhi
     101//while ( (nucleon=theNucleus->GetNextNucleon())&& (InterNumber < 1) ) // Uzhi
     102        while ( (nucleon=theNucleus->GetNextNucleon()) ) // Uzhi
    118103        {
    119104           G4double impact2= sqr(impactX - nucleon->GetPosition().x()) +
    120                     sqr(impactY - nucleon->GetPosition().y());
    121 //         if ( theCrossSection.GetInelasticProbability(s,impact2)                       // Uzhi
    122            if ( theCrossSection.GetInelasticProbability(  impact2/fermi/fermi)           // Uzhi
     105                             sqr(impactY - nucleon->GetPosition().y());
     106
     107//         if ( theParameters->GetInelasticProbability(impact2/fermi/fermi) // Uzhi 29.03.08
     108           if ( theParameters->GetProbabilityOfInteraction(impact2/fermi/fermi) // Uzhi 29.03.08
    123109                > G4UniformRand() )
    124110           {
     111//InterNumber++;
    125112                if ( nucleusNeedsShift )
    126113                {                       // on the first hit, shift nucleus
    … …  
    136123                    nucleon->Hit(targetSplitable);
    137124                }
    138                 G4InteractionContent * aInteraction = new G4InteractionContent(primarySplitable);
     125                G4InteractionContent * aInteraction =
     126                                       new G4InteractionContent(primarySplitable);
    139127                aInteraction->SetTarget(targetSplitable);
    140128                theInteractions.push_back(aInteraction);
    … …  
    152140
    153141// Implementation (private) methods
    154 
    155 
    156 
    157 
    158 
Note: See TracChangeset for help on using the changeset viewer.