G4DiffractiveExcitation Class Reference

#include <G4DiffractiveExcitation.hh>


Public Member Functions

 G4DiffractiveExcitation ()
virtual ~G4DiffractiveExcitation ()
virtual G4bool ExciteParticipants (G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4FTFParameters *theParameters, G4ElasticHNScattering *theElastic) const
virtual void CreateStrings (G4VSplitableHadron *aHadron, G4bool isProjectile, G4ExcitedString *&FirstString, G4ExcitedString *&SecondString, G4FTFParameters *theParameters) const


Detailed Description

Definition at line 48 of file G4DiffractiveExcitation.hh.


Constructor & Destructor Documentation

G4DiffractiveExcitation::G4DiffractiveExcitation (  ) 

Definition at line 68 of file G4DiffractiveExcitation.cc.

00069 {
00070 }

G4DiffractiveExcitation::~G4DiffractiveExcitation (  )  [virtual]

Definition at line 1554 of file G4DiffractiveExcitation.cc.

01555 {
01556 }


Member Function Documentation

void G4DiffractiveExcitation::CreateStrings ( G4VSplitableHadron aHadron,
G4bool  isProjectile,
G4ExcitedString *&  FirstString,
G4ExcitedString *&  SecondString,
G4FTFParameters theParameters 
) const [virtual]

Definition at line 1076 of file G4DiffractiveExcitation.cc.

References G4cout, G4endl, G4UniformRand, G4Parton::Get4Momentum(), G4VSplitableHadron::Get4Momentum(), G4Parton::GetDefinition(), G4VSplitableHadron::GetNextParton(), G4ParticleDefinition::GetParticleSubType(), G4ParticleDefinition::GetParticleType(), G4Parton::GetPDGcode(), G4ParticleDefinition::GetPDGEncoding(), G4VSplitableHadron::GetPosition(), G4FTFParameters::GetProjMinDiffMass(), G4FTFParameters::GetPt2Kink(), G4FTFParameters::GetQuarkProbabilitiesAtGluonSplitUp(), G4VSplitableHadron::GetStatus(), G4FTFParameters::GetTarMinDiffMass(), G4VSplitableHadron::GetTimeOfCreation(), G4INCL::Math::pi, G4Parton::Set4Momentum(), G4ExcitedString::SetPosition(), G4ExcitedString::SetTimeOfCreation(), G4VSplitableHadron::SplitUp(), and sqr().

01081 {
01082 /*
01083 G4cout<<"Create Strings SplitUp "<<hadron<<G4endl;
01084 G4cout<<"Defin "<<hadron->GetDefinition()<<G4endl;
01085 G4cout<<"Defin "<<hadron->GetDefinition()->GetPDGEncoding()<<G4endl;
01086 */
01087         hadron->SplitUp();
01088         G4Parton *start= hadron->GetNextParton();
01089 
01090         if ( start==NULL)
01091         { G4cout << " G4FTFModel::String() Error:No start parton found"<< G4endl;
01092           FirstString=0; SecondString=0;
01093           return;
01094         }
01095         G4Parton *end  = hadron->GetNextParton();
01096         if ( end==NULL)
01097         { G4cout << " G4FTFModel::String() Error:No end parton found"<< G4endl;
01098           FirstString=0; SecondString=0;
01099           return;
01100         }
01101 //G4cout<<start<<" "<<start->GetPDGcode()<<" "<<end<<" "<<end->GetPDGcode()<<G4endl;
01102 //G4cout<<"Create string "<<start->GetPDGcode()<<" "<<end->GetPDGcode()<<G4endl;
01103         G4LorentzVector Phadron=hadron->Get4Momentum();
01104 //G4cout<<"String mom "<<Phadron<<G4endl;
01105         G4LorentzVector Pstart(0.,0.,0.,0.);
01106         G4LorentzVector Pend(0.,0.,0.,0.);
01107         G4LorentzVector Pkink(0.,0.,0.,0.);
01108         G4LorentzVector PkinkQ1(0.,0.,0.,0.);
01109         G4LorentzVector PkinkQ2(0.,0.,0.,0.);
01110 
01111         G4int PDGcode_startQ = std::abs(start->GetDefinition()->GetPDGEncoding());
01112         G4int PDGcode_endQ   = std::abs(  end->GetDefinition()->GetPDGEncoding());
01113 //G4cout<<"PDGcode_startQ "<<PDGcode_startQ<<" PDGcode_endQ   "<<PDGcode_endQ  <<G4endl;
01114 
01115 //--------------------------------------------------------------------------------        
01116         G4double Wmin(0.);
01117         if(isProjectile)
01118         {
01119           Wmin=theParameters->GetProjMinDiffMass();
01120         } else
01121         {
01122           Wmin=theParameters->GetTarMinDiffMass();
01123         } // end of if(isProjectile)
01124 
01125         G4double W = hadron->Get4Momentum().mag();
01126 //G4cout<<"Wmin W "<<Wmin<<" "<<W<<G4endl;
01127         G4double W2=W*W;
01128 
01129         G4double Pt(0.), x1(0.), x3(0.); // x2(0.), 
01130         G4bool Kink=false;
01131 
01132        if(!(((start->GetDefinition()->GetParticleSubType() == "di_quark") &&
01133              (  end->GetDefinition()->GetParticleSubType() == "di_quark")   ) ||
01134             ((start->GetDefinition()->GetParticleSubType() == "quark")    &&
01135              (  end->GetDefinition()->GetParticleSubType() == "quark")      )))
01136        {   // Kinky strings are allowed only for qq-q strings
01137            // Kinky strings are impossible for other systems (qq-qqbar, q-qbar)
01138            // according to the analysis of Pbar P interactions
01139 //G4cout<<G4endl<<"Check for Kink!##############"<<G4endl<<G4endl;
01140         if(W > Wmin)
01141         {                                        // Kink is possible
01142           if(hadron->GetStatus() == 0)           // VU 10.04.2012
01143           {
01144            G4double Pt2kink=theParameters->GetPt2Kink(); // For non-diffractive
01145            Pt = std::sqrt(Pt2kink*(std::pow(W2/16./Pt2kink+1.,G4UniformRand()) - 1.));
01146           }
01147           else
01148           {Pt=0.;}                                       // For diffractive
01149 
01150           if(Pt > 500.*MeV)
01151           {
01152              G4double Ymax = std::log(W/2./Pt + std::sqrt(W2/4./Pt/Pt - 1.));
01153              G4double Y=Ymax*(1.- 2.*G4UniformRand());
01154 
01155              x1=1.-Pt/W*std::exp( Y);
01156              x3=1.-Pt/W*std::exp(-Y);
01157 //             x2=2.-x1-x3;
01158 
01159              G4double Mass_startQ = 650.*MeV;
01160              if(PDGcode_startQ <  3) Mass_startQ =  325.*MeV;
01161              if(PDGcode_startQ == 3) Mass_startQ =  500.*MeV;
01162              if(PDGcode_startQ == 4) Mass_startQ = 1600.*MeV;
01163 
01164              G4double Mass_endQ = 650.*MeV;
01165              if(PDGcode_endQ <  3) Mass_endQ =  325.*MeV;
01166              if(PDGcode_endQ == 3) Mass_endQ =  500.*MeV;
01167              if(PDGcode_endQ == 4) Mass_endQ = 1600.*MeV;
01168 
01169              G4double P2_1=W2*x1*x1/4.-Mass_endQ  *Mass_endQ;
01170              G4double P2_3=W2*x3*x3/4.-Mass_startQ*Mass_startQ;
01171      
01172              G4double P2_2=sqr((2.-x1-x3)*W/2.);
01173 
01174              if((P2_1 <= 0.) || (P2_3 <= 0.))
01175              { Kink=false;}
01176              else
01177              {
01178                G4double P_1=std::sqrt(P2_1);
01179                G4double P_2=std::sqrt(P2_2);
01180                G4double P_3=std::sqrt(P2_3);
01181 
01182                G4double CosT12=(P2_3-P2_1-P2_2)/(2.*P_1*P_2);
01183                G4double CosT13=(P2_2-P2_1-P2_3)/(2.*P_1*P_3);
01184 //             Pt=P_2*std::sqrt(1.-CosT12*CosT12);  // because system was rotated 11.12.09
01185 
01186                if((std::abs(CosT12) >1.) || (std::abs(CosT13) > 1.)) 
01187                { Kink=false;}
01188                else
01189                { 
01190                  Kink=true; 
01191                  Pt=P_2*std::sqrt(1.-CosT12*CosT12);  // because system was rotated 11.12.09
01192                  Pstart.setPx(-Pt); Pstart.setPy(0.); Pstart.setPz(P_3*CosT13); 
01193                  Pend.setPx(   0.); Pend.setPy(  0.); Pend.setPz(          P_1); 
01194                  Pkink.setPx(  Pt); Pkink.setPy( 0.); Pkink.setPz(  P_2*CosT12);
01195                  Pstart.setE(x3*W/2.);                
01196                  Pkink.setE(Pkink.vect().mag());
01197                  Pend.setE(x1*W/2.);
01198 
01199                  G4double XkQ=GetQuarkFractionOfKink(0.,1.);
01200                  if(Pkink.getZ() > 0.) 
01201                  {
01202                    if(XkQ > 0.5) {PkinkQ1=XkQ*Pkink;} else {PkinkQ1=(1.-XkQ)*Pkink;}
01203                  } else {
01204                    if(XkQ > 0.5) {PkinkQ1=(1.-XkQ)*Pkink;} else {PkinkQ1=XkQ*Pkink;}
01205                  }
01206 
01207                  PkinkQ2=Pkink - PkinkQ1;
01208 //------------------------- Minimizing Pt1^2+Pt3^2 ------------------------------
01209 
01210                  G4double Cos2Psi=(sqr(x1) -sqr(x3)+2.*sqr(x3*CosT13))/
01211                           std::sqrt(sqr(sqr(x1)-sqr(x3)) + sqr(2.*x1*x3*CosT13));
01212                  G4double Psi=std::acos(Cos2Psi);
01213 
01214 G4LorentzRotation Rotate;
01215 if(isProjectile) {Rotate.rotateY(Psi);}
01216 else             {Rotate.rotateY(pi-Psi);}                   
01217 Rotate.rotateZ(twopi*G4UniformRand());
01218 
01219 Pstart*=Rotate;
01220 Pkink*=Rotate;
01221 PkinkQ1*=Rotate;
01222 PkinkQ2*=Rotate;
01223 Pend*=Rotate;
01224 
01225                }
01226              }      // end of if((P2_1 < 0.) || (P2_3 <0.))
01227           }         // end of if(Pt > 500.*MeV)
01228         }           // end of if(W > Wmin) Check for a kink
01229        }            // end of qq-q string selection
01230 //--------------------------------------------------------------------------------
01231 /*
01232 G4cout<<"Kink "<<Kink<<" "
01233 <<start->GetDefinition()->GetParticleSubType()<<" "
01234 <<  end->GetDefinition()->GetParticleSubType() <<G4endl;
01235 
01236 G4cout<<"Kink "<<Kink<<" "
01237 <<start->GetDefinition()->GetPDGEncoding()<<" "
01238 <<  end->GetDefinition()->GetPDGEncoding() <<G4endl;
01239 G4int Uzhi; G4cin>>Uzhi;
01240 */
01241 
01242         if(Kink)
01243         {                                        // Kink is possible
01244 //G4cout<<"Kink is sampled!"<<G4endl;
01245           std::vector<G4double> QuarkProbabilitiesAtGluonSplitUp =
01246               theParameters->GetQuarkProbabilitiesAtGluonSplitUp();
01247 
01248           G4int QuarkInGluon(1); G4double Ksi=G4UniformRand();
01249           for(unsigned int Iq=0; Iq <3; Iq++)
01250           {
01251 //G4cout<<"Iq "<<Iq<<G4endl;
01252 
01253 if(Ksi > QuarkProbabilitiesAtGluonSplitUp[Iq]) QuarkInGluon++;}
01254 
01255 //G4cout<<"Last Iq "<<QuarkInGluon<<G4endl;
01256           G4Parton * Gquark = new G4Parton(QuarkInGluon);
01257           G4Parton * Ganti_quark = new G4Parton(-QuarkInGluon);
01258 //G4cout<<"Lorentz "<<G4endl;
01259 
01260 //-------------------------------------------------------------------------------
01261           G4LorentzRotation toCMS(-1*Phadron.boostVector());
01262 
01263           G4LorentzRotation toLab(toCMS.inverse());
01264 //G4cout<<"Pstart "<<Pstart<<G4endl;
01265 //G4cout<<"Pend   "<<Pend<<G4endl;
01266           Pstart.transform(toLab);  start->Set4Momentum(Pstart);
01267           PkinkQ1.transform(toLab);
01268           PkinkQ2.transform(toLab);
01269           Pend.transform(toLab);    end->Set4Momentum(Pend);
01270 //G4cout<<"Pstart "<<Pstart<<G4endl;
01271 //G4cout<<"Pend   "<<Pend<<G4endl;
01272 //G4cout<<"Defin "<<hadron->GetDefinition()<<G4endl;
01273 //G4cout<<"Defin "<<hadron->GetDefinition()->GetPDGEncoding()<<G4endl;
01274 
01275 //        G4int absPDGcode=std::abs(hadron->GetDefinition()->GetPDGEncoding());
01276           G4int absPDGcode=1500;  // 23 Dec
01277 if((start->GetDefinition()->GetParticleSubType() == "quark") &&
01278    (  end->GetDefinition()->GetParticleSubType() == "quark")  )
01279   absPDGcode=110;
01280 
01281 //G4cout<<"absPDGcode "<<absPDGcode<<G4endl;
01282 
01283           if(absPDGcode < 1000)
01284           {                                // meson
01285             if ( isProjectile )
01286             {                              // Projectile
01287               if(end->GetDefinition()->GetPDGEncoding() > 0 )  // A quark on the end
01288               {                            // Quark on the end
01289                 FirstString = new G4ExcitedString(end   ,Ganti_quark, +1);
01290                 SecondString= new G4ExcitedString(Gquark,start      ,+1);
01291                 Ganti_quark->Set4Momentum(PkinkQ1);
01292                 Gquark->Set4Momentum(PkinkQ2);
01293 
01294               } else
01295               {                            // Anti_Quark on the end
01296                 FirstString = new G4ExcitedString(end        ,Gquark, +1);
01297                 SecondString= new G4ExcitedString(Ganti_quark,start ,+1);
01298                 Gquark->Set4Momentum(PkinkQ1);
01299                 Ganti_quark->Set4Momentum(PkinkQ2);
01300 
01301               }   // end of if(end->GetPDGcode() > 0)
01302             } else {                      // Target
01303               if(end->GetDefinition()->GetPDGEncoding() > 0 )  // A quark on the end
01304               {                           // Quark on the end
01305                 FirstString = new G4ExcitedString(Ganti_quark,end   ,-1);
01306                 SecondString= new G4ExcitedString(start      ,Gquark,-1);
01307                 Ganti_quark->Set4Momentum(PkinkQ2);
01308                 Gquark->Set4Momentum(PkinkQ1);
01309 
01310               } else
01311               {                            // Anti_Quark on the end
01312                 FirstString = new G4ExcitedString(Gquark,end        ,-1);
01313                 SecondString= new G4ExcitedString(start ,Ganti_quark,-1);
01314                 Gquark->Set4Momentum(PkinkQ2);
01315                 Ganti_quark->Set4Momentum(PkinkQ1);
01316 
01317               }   // end of if(end->GetPDGcode() > 0)
01318             }     // end of if ( isProjectile )
01319           } else  // if(absPDGCode < 1000)
01320           {                             // Baryon/AntiBaryon
01321             if ( isProjectile )
01322             {                              // Projectile
01323               if((end->GetDefinition()->GetParticleType() == "diquarks") &&
01324                  (end->GetDefinition()->GetPDGEncoding() > 0           )   ) 
01325               {                            // DiQuark on the end
01326                 FirstString = new G4ExcitedString(end        ,Gquark, +1);
01327                 SecondString= new G4ExcitedString(Ganti_quark,start ,+1);
01328                 Gquark->Set4Momentum(PkinkQ1);
01329                 Ganti_quark->Set4Momentum(PkinkQ2);
01330 
01331               } else
01332               {                            // Anti_DiQuark on the end or quark
01333                 FirstString = new G4ExcitedString(end   ,Ganti_quark, +1);
01334                 SecondString= new G4ExcitedString(Gquark,start      ,+1);
01335                 Ganti_quark->Set4Momentum(PkinkQ1);
01336                 Gquark->Set4Momentum(PkinkQ2);
01337 
01338               }   // end of if(end->GetPDGcode() > 0)
01339             } else {                      // Target
01340 
01341               if((end->GetDefinition()->GetParticleType() == "diquarks") &&
01342                  (end->GetDefinition()->GetPDGEncoding() > 0           )   ) 
01343               {                            // DiQuark on the end
01344                 FirstString = new G4ExcitedString(Gquark,end        ,-1);
01345 
01346                 SecondString= new G4ExcitedString(start ,Ganti_quark,-1);
01347                 Gquark->Set4Momentum(PkinkQ1);
01348                 Ganti_quark->Set4Momentum(PkinkQ2);
01349 
01350               } else
01351               {                            // Anti_DiQuark on the end or Q
01352                 FirstString = new G4ExcitedString(Ganti_quark,end   ,-1);
01353                 SecondString= new G4ExcitedString(start      ,Gquark,-1);
01354                 Gquark->Set4Momentum(PkinkQ2);
01355                 Ganti_quark->Set4Momentum(PkinkQ1);
01356 
01357               }   // end of if(end->GetPDGcode() > 0)
01358             }     // end of if ( isProjectile )
01359           }  // end of if(absPDGcode < 1000)
01360 
01361           FirstString->SetTimeOfCreation(hadron->GetTimeOfCreation());
01362           FirstString->SetPosition(hadron->GetPosition());
01363 
01364           SecondString->SetTimeOfCreation(hadron->GetTimeOfCreation());
01365           SecondString->SetPosition(hadron->GetPosition());
01366 
01367 // -------------------------------------------------------------------------  
01368         } else                                   // End of kink is possible
01369         {                                        // Kink is impossible
01370 //G4cout<<start<<" "<<start->GetPDGcode()<<" "<<end<<" "<<end->GetPDGcode()<<G4endl;
01371           if ( isProjectile )
01372           {
01373                 FirstString= new G4ExcitedString(end,start, +1);
01374           } else {
01375                 FirstString= new G4ExcitedString(start,end, -1);
01376           }
01377 
01378           SecondString=0;
01379 
01380           FirstString->SetTimeOfCreation(hadron->GetTimeOfCreation());
01381           FirstString->SetPosition(hadron->GetPosition());
01382 
01383 // momenta of string ends
01384 //
01385           G4double Momentum=hadron->Get4Momentum().vect().mag();
01386           G4double  Plus=hadron->Get4Momentum().e() + Momentum;
01387           G4double Minus=hadron->Get4Momentum().e() - Momentum;
01388 
01389           G4ThreeVector tmp;
01390           if(Momentum > 0.) 
01391           {
01392            tmp.set(hadron->Get4Momentum().px(),
01393                    hadron->Get4Momentum().py(),
01394                    hadron->Get4Momentum().pz());
01395            tmp/=Momentum;
01396           }
01397           else
01398           {
01399            tmp.set(0.,0.,1.);
01400           }
01401 
01402           G4LorentzVector Pstart1(tmp,0.);
01403           G4LorentzVector   Pend1(tmp,0.);
01404 
01405           if(isProjectile)
01406           {
01407            Pstart1*=(-1.)*Minus/2.;
01408            Pend1  *=(+1.)*Plus /2.;
01409           } 
01410           else
01411           {
01412            Pstart1*=(+1.)*Plus/2.;
01413            Pend1  *=(-1.)*Minus/2.;
01414           }
01415 
01416           Momentum=-Pstart1.mag();
01417           Pstart1.setT(Momentum);  // It is assumed that quark has m=0.
01418 
01419           Momentum=-Pend1.mag();
01420           Pend1.setT(Momentum);    // It is assumed that di-quark has m=0.
01421 
01422           start->Set4Momentum(Pstart1);
01423           end->Set4Momentum(Pend1);
01424           SecondString=0;
01425         }            // End of kink is impossible 
01426 
01427 //G4cout<<"Quarks in the string at creation"<<FirstString->GetRightParton()->GetPDGcode()<<" "<<FirstString->GetLeftParton()->GetPDGcode()<<G4endl;
01428 //G4cout<<FirstString<<" "<<SecondString<<G4endl;
01429 
01430 #ifdef G4_FTFDEBUG
01431           G4cout << " generated string flavors          " 
01432                  << start->GetPDGcode() << " / " 
01433                  << end->GetPDGcode()                    << G4endl;
01434           G4cout << " generated string momenta:   quark " 
01435                  << start->Get4Momentum() << "mass : " 
01436                  <<start->Get4Momentum().mag()           << G4endl;
01437           G4cout << " generated string momenta: Diquark " 
01438                  << end ->Get4Momentum() 
01439                  << "mass : " <<end->Get4Momentum().mag()<< G4endl;
01440           G4cout << " sum of ends                       " << Pstart+Pend << G4endl;
01441           G4cout << " Original                          " << hadron->Get4Momentum() << G4endl;
01442 #endif
01443 
01444         return;
01445    
01446 }

G4bool G4DiffractiveExcitation::ExciteParticipants ( G4VSplitableHadron aPartner,
G4VSplitableHadron bPartner,
G4FTFParameters theParameters,
G4ElasticHNScattering theElastic 
) const [virtual]

Definition at line 74 of file G4DiffractiveExcitation.cc.

References G4ElasticHNScattering::ElasticScattering(), G4ParticleTable::FindParticle(), G4UniformRand, G4VSplitableHadron::Get4Momentum(), G4FTFParameters::GetAveragePt2(), G4VSplitableHadron::GetDefinition(), G4FTFParameters::GetDeltaProbAtQuarkExchange(), G4FTFParameters::GetMagQuarkExchange(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetPDGEncoding(), G4ParticleDefinition::GetPDGiIsospin(), G4ParticleDefinition::GetPDGMass(), G4VSplitableHadron::GetPosition(), G4FTFParameters::GetProbabilityOfProjDiff(), G4FTFParameters::GetProbabilityOfTarDiff(), G4FTFParameters::GetProbLogDistr(), G4FTFParameters::GetProbOfSameQuarkExchange(), G4FTFParameters::GetProjMinDiffMass(), G4FTFParameters::GetProjMinNonDiffMass(), G4FTFParameters::GetSlopeQuarkExchange(), G4VSplitableHadron::GetStatus(), G4FTFParameters::GetTarMinDiffMass(), G4FTFParameters::GetTarMinNonDiffMass(), G4VSplitableHadron::GetTimeOfCreation(), G4VSplitableHadron::IncrementCollisionCount(), G4VSplitableHadron::Set4Momentum(), G4VSplitableHadron::SetDefinition(), G4VSplitableHadron::SetPosition(), G4VSplitableHadron::SetStatus(), and G4VSplitableHadron::SetTimeOfCreation().

00078 {
00079 //G4cout<<G4endl<<"ExciteParticipants --------------"<<G4endl;
00080 // -------------------- Projectile parameters -----------------------
00081      G4LorentzVector Pprojectile=projectile->Get4Momentum();
00082 
00083      if(Pprojectile.z() < 0.)
00084      {
00085        target->SetStatus(2);
00086        return false;
00087      } 
00088 
00089      G4double ProjectileRapidity = Pprojectile.rapidity();
00090 
00091      G4int    ProjectilePDGcode=projectile->GetDefinition()->GetPDGEncoding();
00092      G4int    absProjectilePDGcode=std::abs(ProjectilePDGcode);
00093 
00094      G4bool PutOnMassShell(false);
00095 //   G4double M0projectile=projectile->GetDefinition()->GetPDGMass(); // With de-excitation
00096      G4double M0projectile = Pprojectile.mag();                       // Without de-excitation
00097 //G4cout<<"M0projectile "<<M0projectile<<" "<<ProjectileRapidity<<G4endl;
00098 
00099      if(M0projectile < projectile->GetDefinition()->GetPDGMass())
00100      {
00101       PutOnMassShell=true;
00102       M0projectile=projectile->GetDefinition()->GetPDGMass();
00103      }
00104 
00105      G4double M0projectile2 = M0projectile * M0projectile;
00106 
00107      G4double ProjectileDiffStateMinMass=theParameters->GetProjMinDiffMass();
00108      G4double ProjectileNonDiffStateMinMass=theParameters->GetProjMinNonDiffMass();
00109      G4double ProbProjectileDiffraction=theParameters->GetProbabilityOfProjDiff();
00110 
00111 // -------------------- Target parameters -------------------------
00112      G4int    TargetPDGcode=target->GetDefinition()->GetPDGEncoding();
00113      G4int    absTargetPDGcode=std::abs(TargetPDGcode);
00114 //G4cout<<"Entry to QE or Excit "<<ProjectilePDGcode<<" "<<TargetPDGcode<<G4endl;
00115 
00116      G4LorentzVector Ptarget=target->Get4Momentum();
00117 //G4cout<<"Pproj "<<Pprojectile<<G4endl;
00118 //G4cout<<"Ptarget "<<Ptarget<<G4endl;
00119      G4double M0target = Ptarget.mag();
00120 
00121 //   G4double TargetRapidity = Ptarget.rapidity();
00122 //G4cout<<"M0target "<<M0target<<" "<<TargetRapidity<<G4endl;
00123      if(M0target < target->GetDefinition()->GetPDGMass())
00124      {
00125       PutOnMassShell=true;
00126       M0target=target->GetDefinition()->GetPDGMass();
00127      }
00128 
00129      G4double M0target2 = M0target * M0target; 
00130  
00131      G4double TargetDiffStateMinMass=theParameters->GetTarMinDiffMass();    
00132      G4double TargetNonDiffStateMinMass=theParameters->GetTarMinNonDiffMass();    
00133      G4double ProbTargetDiffraction=theParameters->GetProbabilityOfTarDiff();
00134 
00135      G4double AveragePt2=theParameters->GetAveragePt2();
00136      G4double ProbLogDistr=theParameters->GetProbLogDistr();         // Uzhi 21.05.2012
00137 
00138 //     G4double ProbOfDiffraction=ProbProjectileDiffraction +
00139 //                                ProbTargetDiffraction;
00140 
00141 //     G4double SumMasses=M0projectile+M0target+220.*MeV; // 200->220 7 June 2011
00142      G4double SumMasses=M0projectile+M0target+220.*MeV; // 200->220 7 June 2011
00143 
00144 // Kinematical properties of the interactions --------------
00145      G4LorentzVector Psum;      // 4-momentum in CMS
00146      Psum=Pprojectile+Ptarget;
00147      G4double S=Psum.mag2(); 
00148 
00149 //Uzhi_SqrtS=std::sqrt(S);
00150 
00151 // Transform momenta to cms and then rotate parallel to z axis;
00152      G4LorentzRotation toCms(-1*Psum.boostVector());
00153 
00154      G4LorentzVector Ptmp=toCms*Pprojectile;
00155 
00156      if ( Ptmp.pz() <= 0. )
00157      {
00158        target->SetStatus(2); 
00159    // "String" moving backwards in  CMS, abort collision !!
00160        return false;
00161      }
00162 
00163      toCms.rotateZ(-1*Ptmp.phi());
00164      toCms.rotateY(-1*Ptmp.theta());
00165 
00166      G4LorentzRotation toLab(toCms.inverse());
00167 
00168      Pprojectile.transform(toCms);
00169      Ptarget.transform(toCms);
00170 
00171      G4double PZcms2, PZcms;
00172 
00173      G4double SqrtS=std::sqrt(S);
00174 
00175 /*
00176 G4cout<<"Proj "<<Pprojectile<<G4endl;
00177 G4cout<<"Targ "<<Ptarget<<G4endl;
00178 G4cout<<"SqrtS     "<<SqrtS<<G4endl;
00179 G4cout<<"M0pr M0tr "<<M0projectile<<" "<<M0target<<" "<<SumMasses<<G4endl;
00180 */
00181 //G4cout<<"SqrtS < 2300*MeV Bary "<<SqrtS<<G4endl;
00182 //   if(absProjectilePDGcode > 1000 && (SqrtS < 2300*MeV || SqrtS < SumMasses)) // 7.06.11
00183      if(absProjectilePDGcode > 1000 && (SqrtS < SumMasses))
00184      {target->SetStatus(2);  return false;}  // The model cannot work for
00185                                              // p+p-interactions
00186                                              // at Plab < 1.62 GeV/c.
00187 
00188 //G4cout<<"SqrtS < 1600*MeV Pion "<<SqrtS<<G4endl;
00189 
00190 //   if(( absProjectilePDGcode == 211 || ProjectilePDGcode ==  111) &&        // 7.06.11
00191 //     ((SqrtS < 1600*MeV) || (SqrtS < SumMasses)))
00192      if(( absProjectilePDGcode == 211 || ProjectilePDGcode ==  111) && 
00193        (SqrtS < SumMasses))
00194      {target->SetStatus(2);  return false;}  // The model cannot work for
00195                                              // Pi+p-interactions
00196                                              // at Plab < 1. GeV/c.
00197 //G4cout<<"SqrtS < 1600*MeV "<<SqrtS<<G4endl;
00198 //SumMasses=M0projectile+M0target+20.*MeV;
00199      if(( (absProjectilePDGcode == 321) || (ProjectilePDGcode == -311)   ||
00200           (absProjectilePDGcode == 311) || (absProjectilePDGcode == 130) ||
00201           (absProjectilePDGcode == 310)) && (SqrtS < SumMasses)) 
00202 //        (absProjectilePDGcode == 310)) && ((SqrtS < 1600*MeV) || (SqrtS < SumMasses))) 
00203      {target->SetStatus(2);  return false;}  // The model cannot work for
00204                                              // K+p-interactions
00205                                              // at Plab < ??? GeV/c.  ???
00206 
00207      PZcms2=(S*S+M0projectile2*M0projectile2+M0target2*M0target2-
00208              2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2)
00209              /4./S;
00210 //G4cout<<"PZcms2 "<<PZcms2<<G4endl;
00211      if(PZcms2 < 0)
00212      {target->SetStatus(2);  return false;}   // It can be in an interaction with 
00213                                               // off-shell nuclear nucleon
00214 
00215      PZcms = std::sqrt(PZcms2);
00216 
00217      if(PutOnMassShell)
00218      {
00219       if(Pprojectile.z() > 0.)
00220       {
00221        Pprojectile.setPz( PZcms);
00222        Ptarget.setPz(    -PZcms);
00223       }
00224       else
00225       {
00226        Pprojectile.setPz(-PZcms);
00227        Ptarget.setPz(     PZcms);
00228       };
00229 
00230       Pprojectile.setE(std::sqrt(M0projectile2                  +
00231                                  Pprojectile.x()*Pprojectile.x()+
00232                                  Pprojectile.y()*Pprojectile.y()+
00233                                  PZcms2));
00234       Ptarget.setE(std::sqrt(M0target2              +
00235                              Ptarget.x()*Ptarget.x()+
00236                              Ptarget.y()*Ptarget.y()+
00237                              PZcms2));
00238      }
00239 
00240 
00241 //G4cout<<"Proj "<<Pprojectile<<G4endl;
00242 //G4cout<<"Targ "<<Ptarget<<G4endl;
00243      G4double maxPtSquare; // = PZcms2;
00244 /*
00245 Uzhi_targetdiffraction=0;
00246 Uzhi_projectilediffraction=0;
00247 Uzhi_Mx2=1.0;
00248 Uzhi_modT=0.;
00249 G4int Uzhi_QE=0;
00250 */
00251 /*
00252 G4cout<<"Start --------------------"<<G4endl;
00253 G4cout<<"Proj "<<M0projectile<<" "<<ProjectileDiffStateMinMass<<"  "<<ProjectileNonDiffStateMinMass<<G4endl;
00254 G4cout<<"Targ "<<M0target    <<" "<<TargetDiffStateMinMass    <<" "<<TargetNonDiffStateMinMass<<G4endl;
00255 G4cout<<"SqrtS "<<SqrtS<<G4endl;
00256 G4cout<<"Rapid "<<ProjectileRapidity<<G4endl; //" "<<TargetRapidity<<G4endl;
00257 */
00258 // Charge exchange can be possible for baryons -----------------
00259 
00260 // Getting the values needed for exchange ----------------------
00261      G4double MagQuarkExchange        =theParameters->GetMagQuarkExchange();//0.045; //
00262      G4double SlopeQuarkExchange      =theParameters->GetSlopeQuarkExchange();//0.; //
00263 // 777777777777777777777777777777777777777777777777777777777777777777777777777777
00264      G4double DeltaProbAtQuarkExchange=theParameters->GetDeltaProbAtQuarkExchange();
00265 
00266 //     G4double NucleonMass=
00267 //              (G4ParticleTable::GetParticleTable()->FindParticle(2112))->GetPDGMass();     
00268      G4double DeltaMass=
00269               (G4ParticleTable::GetParticleTable()->FindParticle(2224))->GetPDGMass();
00270 
00271 //G4double TargetRapidity(0.);
00272 //G4cout<<"Prob Q Exch "<<MagQuarkExchange*std::exp(-SlopeQuarkExchange*ProjectileRapidity)<<G4endl;
00273 
00274 //G4cout<<"Q exc Mag Slop Wdelta"<<MagQuarkExchange<<" "<<SlopeQuarkExchange<<" "<<DeltaProbAtQuarkExchange<<G4endl;
00275 //G4cout<<"ProjectileRapidity "<<ProjectileRapidity<<G4endl;
00276 //G4cout<<"Prob Exc "<<MagQuarkExchange*std::exp(-SlopeQuarkExchange*(ProjectileRapidity))<<G4endl;
00277 //G4int Uzhi; G4cin>>Uzhi;
00278 // Check for possible quark exchange -----------------------------------
00279 
00280      if(G4UniformRand() < MagQuarkExchange*
00281         std::exp(-SlopeQuarkExchange*ProjectileRapidity))  //TargetRapidity))) 1.45
00282      {    
00283 //        std::exp(-SlopeQuarkExchange*(ProjectileRapidity - 1.36)))  //TargetRapidity))) 1.45
00284 //G4cout<<"Q exchange"<<G4endl;
00285 //Uzhi_QE=1;
00286       G4int NewProjCode(0), NewTargCode(0);
00287 
00288       G4int ProjQ1(0), ProjQ2(0), ProjQ3(0);
00289 
00290 //  Projectile unpacking --------------------------
00291       if(absProjectilePDGcode < 1000 )
00292       {    // projectile is meson ----------------- 
00293        UnpackMeson(ProjectilePDGcode, ProjQ1, ProjQ2);  
00294       } else 
00295       {    // projectile is baryon ----------------
00296        UnpackBaryon(ProjectilePDGcode, ProjQ1, ProjQ2, ProjQ3);
00297       } // End of the hadron's unpacking ----------
00298 
00299 //  Target unpacking ------------------------------
00300       G4int TargQ1(0), TargQ2(0), TargQ3(0);
00301       UnpackBaryon(TargetPDGcode, TargQ1, TargQ2, TargQ3); 
00302 
00303 //G4cout<<ProjQ1<<" "<<ProjQ2<<" "<<ProjQ3<<G4endl;
00304 //G4cout<<TargQ1<<" "<<TargQ2<<" "<<TargQ3<<G4endl;
00305 // Sampling of exchanged quarks -------------------
00306       G4int ProjExchangeQ(0);
00307       G4int TargExchangeQ(0);
00308 
00309       if(absProjectilePDGcode < 1000 )
00310       {    // projectile is meson ----------------- 
00311 
00312        if(ProjQ1 > 0 ) // ProjQ1 is quark
00313        {  
00314         G4int Navailable=0;
00315         ProjExchangeQ = ProjQ1;
00316         if(ProjExchangeQ != TargQ1) Navailable++;
00317         if(ProjExchangeQ != TargQ2) Navailable++;
00318         if(ProjExchangeQ != TargQ3) Navailable++;
00319 
00320         G4int Nsampled=CLHEP::RandFlat::shootInt(G4long(Navailable))+1;
00321 //G4cout<<"Navailable Nsampled "<<Navailable<<" "<<Nsampled<<G4endl;
00322         Navailable=0;
00323         if(ProjExchangeQ != TargQ1) 
00324         {
00325          Navailable++;
00326          if(Navailable == Nsampled) 
00327          {TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjQ1=TargExchangeQ;}
00328         }
00329 
00330         if(ProjExchangeQ != TargQ2)
00331         {
00332          Navailable++;
00333          if(Navailable == Nsampled)
00334          {TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjQ1=TargExchangeQ;}
00335         }          
00336 
00337         if(ProjExchangeQ != TargQ3)
00338         {
00339          Navailable++;
00340          if(Navailable == Nsampled)
00341          {TargExchangeQ = TargQ3;  TargQ3=ProjExchangeQ; ProjQ1=TargExchangeQ;}
00342         }
00343        } else          // ProjQ2 is quark
00344        {  
00345         G4int Navailable=0;
00346         ProjExchangeQ = ProjQ2;
00347         if(ProjExchangeQ != TargQ1) Navailable++;
00348         if(ProjExchangeQ != TargQ2) Navailable++;
00349         if(ProjExchangeQ != TargQ3) Navailable++;
00350 
00351         G4int Nsampled=CLHEP::RandFlat::shootInt(G4long(Navailable))+1;
00352 //G4cout<<"Navailable Nsampled "<<Navailable<<" "<<Nsampled<<G4endl;
00353         Navailable=0;
00354         if(ProjExchangeQ != TargQ1) 
00355         {
00356          Navailable++;
00357          if(Navailable == Nsampled) 
00358          {TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjQ2=TargExchangeQ;}
00359         }
00360 
00361         if(ProjExchangeQ != TargQ2)
00362         {
00363          Navailable++;
00364          if(Navailable == Nsampled)
00365          {TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjQ2=TargExchangeQ;}
00366         }
00367 
00368         if(ProjExchangeQ != TargQ3)
00369         {
00370          Navailable++;
00371          if(Navailable == Nsampled)
00372          {TargExchangeQ = TargQ3;  TargQ3=ProjExchangeQ; ProjQ2=TargExchangeQ;}
00373         }
00374        } // End of if(ProjQ1 > 0 ) // ProjQ1 is quark
00375 
00376 //G4cout<<"Exch Pr Tr "<<ProjExchangeQ<<" "<<TargExchangeQ<<G4endl;
00377 //G4cout<<ProjQ1<<" "<<ProjQ2<<" "<<ProjQ3<<G4endl;
00378 //G4cout<<TargQ1<<" "<<TargQ2<<" "<<TargQ3<<G4endl;
00379 
00380        G4int aProjQ1=std::abs(ProjQ1);
00381        G4int aProjQ2=std::abs(ProjQ2);
00382        if(aProjQ1 == aProjQ2)          {NewProjCode = 111;} // Pi0-meson 
00383        else  // |ProjQ1| # |ProjQ2|
00384        {
00385         if(aProjQ1 > aProjQ2)          {NewProjCode = aProjQ1*100+aProjQ2*10+1;}
00386         else                           {NewProjCode = aProjQ2*100+aProjQ1*10+1;}
00387 //      NewProjCode *=(ProjectilePDGcode/absProjectilePDGcode);
00388        }
00389 
00390 G4bool ProjExcited=false;
00391 //G4cout<<"NewProjCode "<<NewProjCode<<G4endl;
00392         if(G4UniformRand() < 0.5) 
00393         {
00394          NewProjCode +=2; // Excited Pi0-meson 
00395          ProjExcited=true;
00396         }
00397         if(aProjQ1 != aProjQ2) NewProjCode *=(ProjectilePDGcode/absProjectilePDGcode); // Uzhi
00398 //G4cout<<"NewProjCode +2 or 0 "<<NewProjCode<<G4endl;
00399 
00400 G4ParticleDefinition* TestParticle=0;
00401 TestParticle=G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode);
00402 //G4cout<<"TestParticle ? "<<TestParticle<<G4endl;
00403 
00404 if(TestParticle) 
00405 {
00406  G4double MtestPart=                             // 31.05.2012
00407  (G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode))->GetPDGMass();
00408 /*
00409 G4cout<<"TestParticle Name "<<NewProjCode<<" "<<TestParticle->GetParticleName()<<G4endl;
00410 G4cout<<"MtestPart M0projectile projectile->GetDefinition()->GetPDGMass() "<<MtestPart<<" "<<M0projectile<<" "<<projectile->GetDefinition()->GetPDGMass()<<G4endl;
00411 G4bool Test =M0projectile <= projectile->GetDefinition()->GetPDGMass(); 
00412 G4cout<<"M0projectile <= projectile->GetDefinition()->GetPDGMass() "<<Test<<G4endl;
00413 */
00414 
00415   if(MtestPart > M0projectile) 
00416   {M0projectile = MtestPart;}
00417   else 
00418   {
00419    if(std::abs(M0projectile - projectile->GetDefinition()->GetPDGMass()) < 140.*MeV)
00420    {M0projectile = MtestPart;}
00421   }
00422 //G4cout<<"M0projectile After check "<<M0projectile<<G4endl;
00423   M0projectile2 = M0projectile * M0projectile;
00424 
00425  ProjectileDiffStateMinMass   =M0projectile+210.*MeV; //210 MeV=m_pi+70 MeV 
00426  ProjectileNonDiffStateMinMass=M0projectile+210.*MeV; //210 MeV=m_pi+70 MeV
00427 } else
00428 {return false;}
00429 
00430 //G4cout<<"New TrQ "<<TargQ1<<" "<<TargQ2<<" "<<TargQ3<<G4endl;
00431        NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3);
00432 //G4cout<<"NewTargCode "<<NewTargCode<<G4endl;
00433 
00434 //       if( (TargQ1 != TargQ2) && (TargQ1 != TargQ3) && (TargQ2 != TargQ3) // Lambda or Sigma0
00435 //       {if(G4UniformRand() < 0.5) NewTargCode=
00436 
00437 
00438        if( (TargQ1 == TargQ2) && (TargQ1 == TargQ3) &&
00439            (SqrtS > M0projectile+DeltaMass))           {NewTargCode +=2;  //Create Delta isobar
00440                                                         ProjExcited=true;}
00441        else if( target->GetDefinition()->GetPDGiIsospin() == 3 )          //Delta was the target
00442        { if(G4UniformRand() > DeltaProbAtQuarkExchange){NewTargCode +=2;  //Save   Delta isobar
00443                                                         ProjExcited=true;}
00444          else                                          {}     // De-excite initial Delta isobar
00445        }
00446 
00447 //       else if((!CreateDelta)                               &&
00448        else if((!ProjExcited)                               &&
00449                (G4UniformRand() < DeltaProbAtQuarkExchange) &&         //Nucleon was the target
00450                (SqrtS > M0projectile+DeltaMass))      {NewTargCode +=2;}  //Create Delta isobar
00451 //       else if( CreateDelta)                          {NewTargCode +=2;}
00452        else                                           {}                 //Save initial nucleon
00453 
00454 //       target->SetDefinition(                                          // Fix 15.12.09
00455 //       G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode));// Fix 15.12.09 
00456 
00457 TestParticle=G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode);
00458 //G4cout<<"New targ "<<NewTargCode<<" "<<TestParticle->GetParticleName()<<G4endl;
00459 if(TestParticle) 
00460 {
00461  G4double MtestPart=                             // 31.05.2012
00462  (G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode))->GetPDGMass();
00463 
00464  if(MtestPart > M0target)
00465  {M0target=MtestPart;}
00466  else
00467  {
00468   if(std::abs(M0target - target->GetDefinition()->GetPDGMass()) < 140.*MeV)
00469   {M0target=MtestPart;}
00470  }
00471 
00472  TargetDiffStateMinMass   =M0target+220.*MeV;         //220 MeV=m_pi+80 MeV;    
00473  TargetNonDiffStateMinMass=M0target+220.*MeV;         //220 MeV=m_pi+80 MeV; 
00474 } else
00475 {return false;}
00476       } else 
00477       {    // projectile is baryon ----------------
00478 //=========================================================================
00479        G4double Same=theParameters->GetProbOfSameQuarkExchange(); //0.3; //0.5; 0.
00480        G4bool ProjDeltaHasCreated(false);
00481        G4bool TargDeltaHasCreated(false);
00482  
00483        G4double Ksi=G4UniformRand();
00484        if(G4UniformRand() < 0.5)     // Sampling exchange quark from proj. or targ.
00485        {   // Sampling exchanged quark from the projectile ---
00486         if( Ksi < 0.333333 ) 
00487         {ProjExchangeQ = ProjQ1;}
00488         else if( (0.333333 <= Ksi) && (Ksi < 0.666667))
00489         {ProjExchangeQ = ProjQ2;}
00490         else
00491         {ProjExchangeQ = ProjQ3;}
00492 
00493 //G4cout<<"ProjExchangeQ "<<ProjExchangeQ<<G4endl;
00494         if((ProjExchangeQ != TargQ1)||(G4UniformRand()<Same)) 
00495         {
00496          TargExchangeQ = TargQ1; TargQ1=ProjExchangeQ; ProjExchangeQ=TargExchangeQ;
00497         } else
00498         if((ProjExchangeQ != TargQ2)||(G4UniformRand()<Same)) 
00499         {
00500          TargExchangeQ = TargQ2; TargQ2=ProjExchangeQ; ProjExchangeQ=TargExchangeQ;
00501         } else 
00502         {
00503          TargExchangeQ = TargQ3;  TargQ3=ProjExchangeQ; ProjExchangeQ=TargExchangeQ;
00504         }
00505 
00506 //G4cout<<"ProjExchangeQ "<<ProjExchangeQ<<G4endl;
00507 //G4cout<<"TargExchangeQ "<<TargExchangeQ<<G4endl;
00508         if( Ksi < 0.333333 ) 
00509         {ProjQ1=ProjExchangeQ;}
00510         else if( (0.333333 <= Ksi) && (Ksi < 0.666667))
00511         {ProjQ2=ProjExchangeQ;}
00512         else
00513         {ProjQ3=ProjExchangeQ;}
00514 
00515        } else
00516        {   // Sampling exchanged quark from the target -------
00517         if( Ksi < 0.333333 ) 
00518         {TargExchangeQ = TargQ1;}
00519         else if( (0.333333 <= Ksi) && (Ksi < 0.666667)) 
00520         {TargExchangeQ = TargQ2;}
00521         else
00522         {TargExchangeQ = TargQ3;}
00523 
00524         if((TargExchangeQ != ProjQ1)||(G4UniformRand()<Same)) 
00525         {
00526          ProjExchangeQ = ProjQ1; ProjQ1=TargExchangeQ; TargExchangeQ=ProjExchangeQ;
00527         } else
00528         if((TargExchangeQ != ProjQ2)||(G4UniformRand()<Same)) 
00529         {
00530          ProjExchangeQ = ProjQ2; ProjQ2=TargExchangeQ; TargExchangeQ=ProjExchangeQ;
00531         } else 
00532         {
00533          ProjExchangeQ = ProjQ3;  ProjQ3=TargExchangeQ; TargExchangeQ=ProjExchangeQ;
00534         }
00535 
00536         if( Ksi < 0.333333 ) 
00537         {TargQ1=TargExchangeQ;}
00538         else if( (0.333333 <= Ksi) && (Ksi < 0.666667)) 
00539         {TargQ2=TargExchangeQ;}
00540         else
00541         {TargQ3=TargExchangeQ;}
00542 
00543        } // End of sampling baryon
00544 
00545        NewProjCode = NewNucleonId(ProjQ1, ProjQ2, ProjQ3); // *****************************
00546 
00547        if((ProjQ1==ProjQ2) && (ProjQ1==ProjQ3)) {NewProjCode +=2; ProjDeltaHasCreated=true;}
00548        else if(projectile->GetDefinition()->GetPDGiIsospin() == 3)// Projectile was Delta
00549        { if(G4UniformRand() > DeltaProbAtQuarkExchange)
00550                                                 {NewProjCode +=2; ProjDeltaHasCreated=true;} 
00551          else                                   {NewProjCode +=0; ProjDeltaHasCreated=false;}
00552        }
00553        else                                                       // Projectile was Nucleon
00554        {
00555         if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS > DeltaMass+M0target)) 
00556                                                 {NewProjCode +=2; ProjDeltaHasCreated=true;}
00557         else                                    {NewProjCode +=0; ProjDeltaHasCreated=false;}
00558        } 
00559 
00560 //G4ParticleDefinition* NewTestParticle=
00561 //                      G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode);
00562 //G4cout<<"TestParticleMass NewTestParticle->GetPDGMass() "<<TestParticleMass<<" "<< NewTestParticle->GetPDGMass()<<G4endl;
00563 //if(TestParticleMass < NewTestParticle->GetPDGMass()) {NewProjCode=TestParticleID;}
00564  
00565 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++=
00566 
00567        NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3); // *****************************
00568 
00569 //G4cout<<"TargQ1, TargQ2, TargQ3 "<<TargQ1<<" "<<TargQ2<<" "<<TargQ3<<" "<<NewTargCode<<G4endl;
00570 
00571 //TestParticleID=NewTargCode;
00572 //TestParticleMass=DBL_MAX;
00573 
00574 //TestParticle=G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode);
00575 //if(TestParticle) TestParticleMass=TestParticle->GetPDGMass(); 
00576 
00577        if((TargQ1==TargQ2) && (TargQ1==TargQ3)) {NewTargCode +=2; TargDeltaHasCreated=true;}  
00578        else if(target->GetDefinition()->GetPDGiIsospin() == 3)    // Target was Delta
00579        { if(G4UniformRand() > DeltaProbAtQuarkExchange)
00580                                                 {NewTargCode +=2; TargDeltaHasCreated=true;} 
00581          else                                   {NewTargCode +=0; TargDeltaHasCreated=false;}
00582        }
00583        else                                                       // Target was Nucleon
00584        {
00585         if((G4UniformRand() < DeltaProbAtQuarkExchange) && (SqrtS > M0projectile+DeltaMass)) 
00586                                                 {NewTargCode +=2; TargDeltaHasCreated=true;}
00587         else                                    {NewTargCode +=0; TargDeltaHasCreated=false;}
00588        }         
00589 
00590 //NewTestParticle=G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode);
00591 //G4cout<<"TestParticleMass NewTestParticle->GetPDGMass() "<<TestParticleMass<<" "<< NewTestParticle->GetPDGMass()<<G4endl;
00592 //if(TestParticleMass < NewTestParticle->GetPDGMass()) {NewTargCode=TestParticleID;}
00593 
00594 //G4cout<<"NewProjCode NewTargCode "<<NewProjCode<<" "<<NewTargCode<<G4endl;
00595 //G4int Uzhi; G4cin>>Uzhi;
00596 
00597        if((absProjectilePDGcode == NewProjCode) && (absTargetPDGcode == NewTargCode))
00598        { // Nothing was changed! It is not right!?
00599        }
00600 // Forming baryons --------------------------------------------------
00601 if(ProjDeltaHasCreated) {ProbProjectileDiffraction=1.; ProbTargetDiffraction=0.;}
00602 if(TargDeltaHasCreated) {ProbProjectileDiffraction=0.; ProbTargetDiffraction=1.;}
00603        if(ProjDeltaHasCreated) 
00604        {
00605         G4double MtestPart=                             // 31.05.2012
00606           (G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode))->GetPDGMass();
00607 
00608         if(MtestPart >= M0projectile)                   // 31.05.2012
00609         {                                               // 31.05.2012
00610          M0projectile = MtestPart;                      // 31.05.2012
00611          M0projectile2 = M0projectile * M0projectile;   // 31.05.2012
00612         }                                               // 31.05.2012
00613 
00614         ProjectileDiffStateMinMass   =M0projectile+210.*MeV; //210 MeV=m_pi+70 MeV
00615         ProjectileNonDiffStateMinMass=M0projectile+210.*MeV; //210 MeV=m_pi+70 MeV
00616        }
00617 
00618 //      if(M0target < 
00619 //         (G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode))->GetPDGMass())
00620        if(TargDeltaHasCreated)
00621        {
00622         G4double MtestPart=                             // 31.05.2012
00623           (G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode))->GetPDGMass();
00624 
00625         if(MtestPart >=M0target)                        // 31.05.2012
00626         {                                               // 31.05.2012
00627          M0target=MtestPart;                            // 31.05.2012
00628          M0target2 = M0target * M0target;               // 31.05.2012
00629         }                                               // 31.05.2012
00630 
00631         TargetDiffStateMinMass   =M0target+210.*MeV;         //210 MeV=m_pi+70 MeV;    
00632         TargetNonDiffStateMinMass=M0target+210.*MeV;         //210 MeV=m_pi+70 MeV; 
00633        }
00634       } // End of if projectile is baryon ---------------------------
00635 
00636 //G4cout<<"At end// NewProjCode "<<NewProjCode<<G4endl;
00637 //G4cout<<"At end// NewTargCode "<<NewTargCode<<G4endl;
00638 
00639 // If we assume that final state hadrons after the charge exchange will be
00640 // in the ground states, we have to put ----------------------------------
00641 //G4cout<<"M0pr M0tr SqS "<<M0projectile<<" "<<M0target<<" "<<SqrtS<<G4endl;
00642 
00643       PZcms2=(S*S+M0projectile2*M0projectile2+M0target2*M0target2-
00644              2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2)
00645              /4./S;
00646 //G4cout<<"PZcms2 1 "<<PZcms2<<G4endl<<G4endl;
00647       if(PZcms2 < 0) {return false;}  // It can be if energy is not sufficient for Delta
00648 //----------------------------------------------------------
00649       projectile->SetDefinition(
00650                   G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode)); 
00651 
00652       target->SetDefinition(
00653                   G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode)); 
00654 //----------------------------------------------------------
00655       PZcms = std::sqrt(PZcms2);
00656 
00657       Pprojectile.setPz( PZcms);
00658       Pprojectile.setE(std::sqrt(M0projectile2+PZcms2));
00659 
00660       Ptarget.setPz(    -PZcms);
00661       Ptarget.setE(std::sqrt(M0target2+PZcms2));
00662 
00663 // ----------------------------------------------------------
00664 
00665       if(absProjectilePDGcode < 1000)
00666       { // For projectile meson
00667        G4double Wexcit=1.-2.256*std::exp(-0.6*ProjectileRapidity);
00668        Wexcit=0.;
00669        if(G4UniformRand() > Wexcit)
00670        {                             // Make elastic scattering
00671 //G4cout<<"Make elastic scattering of new hadrons"<<G4endl;
00672         Pprojectile.transform(toLab);
00673         Ptarget.transform(toLab);
00674 
00675         projectile->SetTimeOfCreation(target->GetTimeOfCreation());
00676         projectile->SetPosition(target->GetPosition());
00677 
00678         projectile->Set4Momentum(Pprojectile);
00679         target->Set4Momentum(Ptarget);
00680 
00681         G4bool Result= theElastic->ElasticScattering (projectile,target,theParameters);
00682 //G4cout<<"Result of el. scatt "<<Result<<G4endl;
00683         return Result;
00684        } // end of if(Make elastic scattering for projectile meson?)
00685       } else
00686       { // For projectile baryon
00687        G4double Wexcit=1.-2.256*std::exp(-0.6*ProjectileRapidity);
00688        //Wexcit=0.;
00689        if(G4UniformRand() > Wexcit)
00690        {                             // Make elastic scattering
00691 //G4cout<<"Make elastic scattering of new hadrons"<<G4endl;
00692         Pprojectile.transform(toLab);
00693         Ptarget.transform(toLab);
00694 
00695         projectile->SetTimeOfCreation(target->GetTimeOfCreation());
00696         projectile->SetPosition(target->GetPosition());
00697 
00698         projectile->Set4Momentum(Pprojectile);
00699         target->Set4Momentum(Ptarget);
00700 
00701         G4bool Result= theElastic->ElasticScattering (projectile,target,theParameters);
00702         return Result;
00703        } // end of if(Make elastic scattering for projectile baryon?)
00704       }
00705 //G4cout<<"Make excitation of new hadrons"<<G4endl;
00706      }  // End of charge exchange part ------------------------------
00707 
00708 // ------------------------------------------------------------------
00709      G4double ProbOfDiffraction=ProbProjectileDiffraction + ProbTargetDiffraction;
00710 /*
00711 G4cout<<"Excite --------------------"<<G4endl;
00712 G4cout<<"Proj "<<M0projectile<<" "<<ProjectileDiffStateMinMass<<"  "<<ProjectileNonDiffStateMinMass<<G4endl;
00713 G4cout<<"Targ "<<M0target    <<" "<<TargetDiffStateMinMass    <<" "<<TargetNonDiffStateMinMass<<G4endl;
00714 G4cout<<"SqrtS "<<SqrtS<<G4endl;
00715 
00716 G4cout<<"Prob ProjDiff TargDiff "<<ProbProjectileDiffraction<<" "<<ProbTargetDiffraction<<" "<<ProbOfDiffraction<<G4endl;
00717 G4cout<<"Pr Y "<<Pprojectile.rapidity()<<" Tr Y "<<Ptarget.rapidity()<<G4endl;
00718 //G4int Uzhi; G4cin>>Uzhi;
00719 */
00720 /*
00721      if(ProjectileNonDiffStateMinMass + TargetNonDiffStateMinMass > SqrtS) // 24.07.10
00722      {
00723       if(ProbOfDiffraction!=0.)
00724       {
00725        ProbProjectileDiffraction/=ProbOfDiffraction;
00726        ProbOfDiffraction=1.;
00727       } else {return false;}      
00728      } 
00729 
00730 */
00731 
00732      if(ProbOfDiffraction!=0.)
00733      {
00734       ProbProjectileDiffraction/=ProbOfDiffraction;
00735      }
00736      else
00737      {
00738       ProbProjectileDiffraction=0.;
00739      }
00740 
00741 //G4cout<<"Prob ProjDiff TargDiff "<<ProbProjectileDiffraction<<" "<<ProbTargetDiffraction<<" "<<ProbOfDiffraction<<G4endl;
00742 
00743      G4double ProjectileDiffStateMinMass2    = ProjectileDiffStateMinMass    *
00744                                                ProjectileDiffStateMinMass;
00745      G4double ProjectileNonDiffStateMinMass2 = ProjectileNonDiffStateMinMass *
00746                                                ProjectileNonDiffStateMinMass;
00747 
00748      G4double TargetDiffStateMinMass2        = TargetDiffStateMinMass        *
00749                                                TargetDiffStateMinMass;
00750      G4double TargetNonDiffStateMinMass2     = TargetNonDiffStateMinMass     *
00751                                                TargetNonDiffStateMinMass;
00752 
00753      G4double Pt2;
00754      G4double ProjMassT2, ProjMassT;
00755      G4double TargMassT2, TargMassT;
00756      G4double PMinusMin, PMinusMax;
00757 //   G4double PPlusMin , PPlusMax;
00758      G4double TPlusMin , TPlusMax;
00759      G4double PMinusNew, PPlusNew, TPlusNew, TMinusNew;
00760 
00761      G4LorentzVector Qmomentum;
00762      G4double Qminus, Qplus;
00763 
00764      G4int whilecount=0;
00765 
00766 //   Choose a process ---------------------------
00767 //ProbOfDiffraction=1.;                 // Uzhi Difr
00768 //ProbProjectileDiffraction=1.;         // Uzhi 
00769      if(G4UniformRand() < ProbOfDiffraction)
00770        {
00771         if(G4UniformRand() < ProbProjectileDiffraction)
00772         { //-------- projectile diffraction ---------------
00773 //G4cout<<"projectile diffraction"<<G4endl;
00774 
00775          do {
00776 /*
00777 Uzhi_projectilediffraction=1;
00778 Uzhi_targetdiffraction=0;
00779 Uzhi_Mx2=1.;
00780 */
00781 //             Generate pt
00782 //             if (whilecount++ >= 500 && (whilecount%100)==0)
00783 //               G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping"
00784 //               << ", loop count/ maxPtSquare : "
00785 //               << whilecount << " / " << maxPtSquare << G4endl;
00786 
00787 //             whilecount++;
00788              if (whilecount > 1000 )
00789              {
00790               Qmomentum=G4LorentzVector(0.,0.,0.,0.);
00791               target->SetStatus(2);  return false;    //  Ignore this interaction
00792              };
00793 
00794 // --------------- Check that the interaction is possible -----------
00795              ProjMassT2=ProjectileDiffStateMinMass2;
00796              ProjMassT =ProjectileDiffStateMinMass;
00797 
00798              TargMassT2=M0target2;
00799              TargMassT =M0target;
00800 //G4cout<<"Masses "<<ProjMassT<<" "<<TargMassT<<" "<<SqrtS<<" "<<ProjMassT+TargMassT<<G4endl;
00801              PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
00802                      2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
00803                     /4./S;
00804 
00805 //G4cout<<"PZcms2 PrD"<<PZcms2<<G4endl;
00806              if(PZcms2 < 0 ) 
00807              {
00808                target->SetStatus(2);  
00809                return false;
00810              }
00811              maxPtSquare=PZcms2;
00812 
00813              Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
00814              Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
00815 
00816              ProjMassT2=ProjectileDiffStateMinMass2+Pt2;
00817              ProjMassT =std::sqrt(ProjMassT2);
00818 
00819              TargMassT2=M0target2+Pt2;
00820              TargMassT =std::sqrt(TargMassT2);
00821 
00822 //G4cout<<"Masses "<<ProjMassT<<" "<<TargMassT<<" "<<SqrtS<<" "<<ProjMassT+TargMassT<<G4endl;
00823 
00824              PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
00825                      2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
00826                     /4./S;
00827 //G4cout<<"PZcms2 PrD"<<PZcms2<<G4endl;
00828              if(PZcms2 < 0 ) continue;
00829              PZcms =std::sqrt(PZcms2);
00830 
00831              PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
00832              PMinusMax=SqrtS-TargMassT;
00833 
00834              PMinusNew=ChooseP(PMinusMin, PMinusMax);
00835 // PMinusNew=1./sqrt(1./PMinusMin-G4UniformRand()*(1./PMinusMin-1./PMinusMax));
00836 //PMinusNew=1./sqr(1./std::sqrt(PMinusMin)-G4UniformRand()*(1./std::sqrt(PMinusMin)-1./std::sqrt(PMinusMax)));
00837 
00838              TMinusNew=SqrtS-PMinusNew;
00839              Qminus=Ptarget.minus()-TMinusNew;
00840              TPlusNew=TargMassT2/TMinusNew;
00841              Qplus=Ptarget.plus()-TPlusNew;
00842 
00843              Qmomentum.setPz( (Qplus-Qminus)/2 );
00844              Qmomentum.setE(  (Qplus+Qminus)/2 );
00845 
00846           } while ((Pprojectile+Qmomentum).mag2() <  ProjectileDiffStateMinMass2); //||  
00847                   //Repeat the sampling because there was not any excitation
00848 //((Ptarget    -Qmomentum).mag2() <  M0target2                  )) );
00849           projectile->SetStatus(1*projectile->GetStatus());          // VU 10.04.2012
00850         }
00851         else
00852         { // -------------- Target diffraction ----------------
00853 
00854 //G4cout<<"Target diffraction"<<G4endl;
00855          do {
00856 /*
00857 Uzhi_projectilediffraction=0;
00858 Uzhi_targetdiffraction=1;
00859 Uzhi_Mx2=1.;
00860 */
00861 //             Generate pt
00862 //             if (whilecount++ >= 500 && (whilecount%100)==0)
00863 //               G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping"
00864 //               << ", loop count/ maxPtSquare : "
00865 //               << whilecount << " / " << maxPtSquare << G4endl;
00866 
00867 //             whilecount++;
00868              if (whilecount > 1000 )
00869              {
00870               Qmomentum=G4LorentzVector(0.,0.,0.,0.);
00871               target->SetStatus(2);  return false;    //  Ignore this interaction
00872              };
00873 //G4cout<<"Qm while "<<Qmomentum<<" "<<whilecount<<G4endl;
00874 // --------------- Check that the interaction is possible -----------
00875              ProjMassT2=M0projectile2;
00876              ProjMassT =M0projectile;
00877 
00878              TargMassT2=TargetDiffStateMinMass2;
00879              TargMassT =TargetDiffStateMinMass;
00880 
00881              PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
00882                      2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
00883                     /4./S;
00884 
00885 //G4cout<<"PZcms2 TrD <0 "<<PZcms2<<" return"<<G4endl;
00886              if(PZcms2 < 0 ) 
00887              {
00888                target->SetStatus(2);  
00889                return false;
00890              }
00891              maxPtSquare=PZcms2;
00892 
00893              Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
00894 
00895 //G4cout<<"Qm while "<<Qmomentum<<" "<<whilecount<<G4endl;
00896              Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
00897 
00898              ProjMassT2=M0projectile2+Pt2;
00899              ProjMassT =std::sqrt(ProjMassT2);
00900 
00901              TargMassT2=TargetDiffStateMinMass2+Pt2;
00902              TargMassT =std::sqrt(TargMassT2);
00903 
00904              PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
00905                      2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
00906                     /4./S;
00907 
00908 //G4cout<<"PZcms2 <0 "<<PZcms2<<" continue"<<G4endl;
00909              if(PZcms2 < 0 ) continue;
00910              PZcms =std::sqrt(PZcms2);
00911 
00912              TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
00913              TPlusMax=SqrtS-ProjMassT;
00914 
00915              TPlusNew=ChooseP(TPlusMin, TPlusMax);
00916 //TPlusNew=1./sqr(1./std::sqrt(TPlusMin)-G4UniformRand()*(1./std::sqrt(TPlusMin)-1./std::sqrt(TPlusMax)));
00917 
00918 //TPlusNew=TPlusMax;
00919 
00920              PPlusNew=SqrtS-TPlusNew;
00921              Qplus=PPlusNew-Pprojectile.plus();
00922              PMinusNew=ProjMassT2/PPlusNew;
00923              Qminus=PMinusNew-Pprojectile.minus();
00924 
00925              Qmomentum.setPz( (Qplus-Qminus)/2 );
00926              Qmomentum.setE(  (Qplus+Qminus)/2 );
00927 
00928 /*
00929 G4cout<<(Pprojectile+Qmomentum).mag()<<" "<<M0projectile<<G4endl;
00930 G4bool First=(Pprojectile+Qmomentum).mag2() <  M0projectile2;
00931 G4cout<<First<<G4endl;
00932 
00933 G4cout<<(Ptarget    -Qmomentum).mag()<<" "<<TargetDiffStateMinMass<<" "<<TargetDiffStateMinMass2<<G4endl;
00934 G4bool Seco=(Ptarget    -Qmomentum).mag2() < TargetDiffStateMinMass2;
00935 G4cout<<Seco<<G4endl;
00936 */
00937 
00938          } while ((Ptarget    -Qmomentum).mag2() <  TargetDiffStateMinMass2);
00939                  // Repeat the sampling because there was not any excitation
00940 // (((Pprojectile+Qmomentum).mag2() <  M0projectile2          ) ||  //No without excitation
00941 //  ((Ptarget    -Qmomentum).mag2() <  TargetDiffStateMinMass2)) );
00942 //G4cout<<"Go out"<<G4endl;
00943           target->SetStatus(1*target->GetStatus());     // VU 10.04.2012
00944          } // End of if(G4UniformRand() < ProbProjectileDiffraction)
00945         }
00946         else  //----------- Non-diffraction process ------------
00947         {
00948 
00949 //G4cout<<"Non-diffraction process"<<G4endl;
00950          do {
00951 /*
00952 Uzhi_projectilediffraction=0;
00953 Uzhi_targetdiffraction=0;
00954 Uzhi_Mx2=1.;
00955 */
00956 //             Generate pt
00957 //             if (whilecount++ >= 500 && (whilecount%100)==0)
00958 //               G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping"
00959 //               << ", loop count/ maxPtSquare : "
00960 //               << whilecount << " / " << maxPtSquare << G4endl;
00961 
00962 //             whilecount++;
00963              if (whilecount > 1000 )
00964              {
00965               Qmomentum=G4LorentzVector(0.,0.,0.,0.);
00966               target->SetStatus(2);  return false;    //  Ignore this interaction
00967              };
00968 // --------------- Check that the interaction is possible -----------
00969              ProjMassT2=ProjectileNonDiffStateMinMass2;
00970              ProjMassT =ProjectileNonDiffStateMinMass;
00971 
00972              TargMassT2=TargetNonDiffStateMinMass2;
00973              TargMassT =TargetNonDiffStateMinMass;
00974 
00975              PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
00976                     2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
00977                    /4./S;
00978 
00979              if(PZcms2 < 0 ) 
00980              {
00981                target->SetStatus(2);  
00982                return false;
00983              }
00984              maxPtSquare=PZcms2;
00985              Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
00986              Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
00987 
00988              ProjMassT2=ProjectileNonDiffStateMinMass2+Pt2;
00989              ProjMassT =std::sqrt(ProjMassT2);
00990 
00991              TargMassT2=TargetNonDiffStateMinMass2+Pt2;
00992              TargMassT =std::sqrt(TargMassT2);
00993 
00994              PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
00995                     2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
00996                    /4./S;
00997 //G4cout<<"PZcms2 ND"<<PZcms2<<G4endl;
00998 
00999              if(PZcms2 < 0 ) continue;
01000              PZcms =std::sqrt(PZcms2);
01001 
01002              PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
01003              PMinusMax=SqrtS-TargMassT;
01004 
01005              if(G4UniformRand() < ProbLogDistr)                       // Uzhi 25.04.2012
01006              {     PMinusNew=ChooseP(PMinusMin, PMinusMax);}               // 12.06.11
01007              else {PMinusNew=(PMinusMax-PMinusMin)*G4UniformRand() + PMinusMin;}
01008              Qminus=PMinusNew-Pprojectile.minus();
01009 
01010              TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
01011              TPlusMax=SqrtS-PMinusNew;                   // Why is it closed???
01012 //             TPlusMax=SqrtS-ProjMassT;      
01013 
01014              if(G4UniformRand() < 0.5) //ProbLogDistr)                 // Uzhi 29.05.2012 0.5)
01015              {     TPlusNew=ChooseP(TPlusMin, TPlusMax);}                   // 12.06.11
01016              else {TPlusNew=(TPlusMax-TPlusMin)*G4UniformRand() +TPlusMin;} 
01017 
01018              Qplus=-(TPlusNew-Ptarget.plus());
01019 
01020              Qmomentum.setPz( (Qplus-Qminus)/2 );
01021              Qmomentum.setE(  (Qplus+Qminus)/2 );
01022 /*
01023 G4cout<<(Pprojectile+Qmomentum).mag2()<<" "<<ProjectileNonDiffStateMinMass2<<G4endl;
01024 G4cout<<(Ptarget    -Qmomentum).mag2()<<" "<<TargetNonDiffStateMinMass2<<G4endl;
01025 G4int Uzhi; G4cin>>Uzhi;
01026 */
01027        } while ( 
01028  ((Pprojectile+Qmomentum).mag2() <  ProjectileNonDiffStateMinMass2) || //No double Diffraction
01029  ((Ptarget    -Qmomentum).mag2() <  TargetNonDiffStateMinMass2    ));
01030 
01031       projectile->SetStatus(0*projectile->GetStatus());     // VU 10.04.2012
01032       target->SetStatus(0*target->GetStatus());             // VU 10.04.2012
01033      }
01034 
01035      Pprojectile += Qmomentum;
01036      Ptarget     -= Qmomentum;
01037 
01038 //G4cout<<"Pr Y "<<Pprojectile.rapidity()<<" Tr Y "<<Ptarget.rapidity()<<G4endl;
01039 
01040 // Transform back and update SplitableHadron Participant.
01041      Pprojectile.transform(toLab);
01042      Ptarget.transform(toLab);
01043 
01044 // Calculation of the creation time ---------------------
01045      projectile->SetTimeOfCreation(target->GetTimeOfCreation());
01046      projectile->SetPosition(target->GetPosition());
01047 // Creation time and position of target nucleon were determined at
01048 // ReggeonCascade() of G4FTFModel
01049 // ------------------------------------------------------
01050 /*
01051 if(Uzhi_projectilediffraction != 0) 
01052 {Uzhi_Mx2=Pprojectile.mag2(); Uzhi_modT=(target->Get4Momentum()-Ptarget).mag2();}
01053 
01054 if(Uzhi_targetdiffraction     != 0) 
01055 {Uzhi_Mx2=Ptarget.mag2(); Uzhi_modT=(projectile->Get4Momentum()-Pprojectile).mag2();}
01056 
01057 if(Uzhi_QE!= 0) 
01058 {
01059  Uzhi_projectilediffraction=0;
01060  Uzhi_targetdiffraction    =0;
01061  Uzhi_Mx2                  =1.;
01062 }
01063 */
01064 //G4cout<<"Mproj "<<Pprojectile.mag()<<G4endl;
01065 //G4cout<<"Mtarg "<<Ptarget.mag()<<G4endl;
01066      projectile->Set4Momentum(Pprojectile);
01067      target->Set4Momentum(Ptarget);
01068 
01069      projectile->IncrementCollisionCount(1);
01070      target->IncrementCollisionCount(1);
01071 
01072      return true;
01073 }


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:51:46 2013 for Geant4 by  doxygen 1.4.7