G4QGSDiffractiveExcitation Class Reference

#include <G4QGSDiffractiveExcitation.hh>

Inheritance diagram for G4QGSDiffractiveExcitation:

Public Member Functions
	G4QGSDiffractiveExcitation ()
virtual	~G4QGSDiffractiveExcitation ()
virtual G4bool	ExciteParticipants (G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner) const
virtual G4ExcitedString *	String (G4VSplitableHadron *aHadron, G4bool isProjectile) const

Detailed Description

Definition at line 51 of file G4QGSDiffractiveExcitation.hh.

Constructor & Destructor Documentation

G4QGSDiffractiveExcitation::G4QGSDiffractiveExcitation

(

)

Definition at line 59 of file G4QGSDiffractiveExcitation.cc.

{
}

G4QGSDiffractiveExcitation::~G4QGSDiffractiveExcitation ( )

virtual

Definition at line 63 of file G4QGSDiffractiveExcitation.cc.

{
}

Member Function Documentation

G4bool G4QGSDiffractiveExcitation::ExciteParticipants ( G4VSplitableHadron *  aPartner, G4VSplitableHadron *  bPartner  ) const

virtual

Reimplemented in G4SingleDiffractiveExcitation.

Definition at line 69 of file G4QGSDiffractiveExcitation.cc.

References CLHEP::HepLorentzVector::boostVector(), G4VSplitableHadron::Get4Momentum(), G4VSplitableHadron::GetDefinition(), G4ParticleDefinition::GetPDGEncoding(), G4ParticleDefinition::GetPDGMass(), python.hepunit::GeV, CLHEP::HepLorentzVector::mag(), CLHEP::HepLorentzVector::mag2(), python.hepunit::MeV, CLHEP::HepLorentzVector::rotateY(), CLHEP::HepLorentzVector::rotateZ(), G4VSplitableHadron::Set4Momentum(), CLHEP::HepLorentzVector::setE(), CLHEP::HepLorentzVector::setPz(), CLHEP::HepLorentzRotation::transform(), CLHEP::HepLorentzVector::transform(), and CLHEP::HepLorentzVector::vect().

Referenced by G4QGSParticipants::SelectInteractions().

{

    G4LorentzVector Pprojectile=projectile->Get4Momentum();

    // -------------------- Projectile parameters -----------------------------------
    G4bool PutOnMassShell=0;

    //     G4double M0projectile=projectile->GetDefinition()->GetPDGMass();  // With de-excitation
    G4double M0projectile = Pprojectile.mag();                        // Without de-excitation

    if(M0projectile < projectile->GetDefinition()->GetPDGMass())
    {
        PutOnMassShell=1;
        M0projectile=projectile->GetDefinition()->GetPDGMass();
    }

    G4double Mprojectile2 = M0projectile * M0projectile;

    G4int    PDGcode=projectile->GetDefinition()->GetPDGEncoding();
    G4int    absPDGcode=std::abs(PDGcode);
    G4double ProjectileDiffCut;
    G4double AveragePt2;

    if( absPDGcode > 1000 )                        //------Projectile is baryon --------
    {
        ProjectileDiffCut = 1.1;                    // GeV
        AveragePt2 = 0.3;                           // GeV^2
    }
    else if( absPDGcode == 211 || PDGcode ==  111) //------Projectile is Pion -----------
    {
        ProjectileDiffCut = 1.0;                   // GeV
        AveragePt2 = 0.3;                          // GeV^2
    }
    else if( absPDGcode == 321 || PDGcode == -311) //------Projectile is Kaon -----------
    {
        ProjectileDiffCut = 1.1;                    // GeV
        AveragePt2 = 0.3;                           // GeV^2
    }
    else                                           //------Projectile is undefined, Nucleon assumed
    {
        ProjectileDiffCut = 1.1;                    // GeV
        AveragePt2 = 0.3;                           // GeV^2
    };

    ProjectileDiffCut = ProjectileDiffCut * GeV;
    AveragePt2 = AveragePt2 * GeV*GeV;

    // -------------------- Target parameters ----------------------------------------------
    G4LorentzVector Ptarget=target->Get4Momentum();

    G4double M0target = Ptarget.mag();

    if(M0target < target->GetDefinition()->GetPDGMass())
    {
        PutOnMassShell=1;
        M0target=target->GetDefinition()->GetPDGMass();
    }

    G4double Mtarget2 = M0target * M0target;                      //Ptarget.mag2(); // for AA-inter.

    G4double NuclearNucleonDiffCut = 1.1*GeV;

    G4double ProjectileDiffCut2 = ProjectileDiffCut * ProjectileDiffCut;
    G4double NuclearNucleonDiffCut2 = NuclearNucleonDiffCut * NuclearNucleonDiffCut;

    // Transform momenta to cms and then rotate parallel to z axis;

    G4LorentzVector Psum;
    Psum=Pprojectile+Ptarget;

    G4LorentzRotation toCms(-1*Psum.boostVector());

    G4LorentzVector Ptmp=toCms*Pprojectile;

    if ( Ptmp.pz() <= 0. )
    {
        // "String" moving backwards in  CMS, abort collision !!
        //G4cout << " abort Collision!! " << G4endl;
        return false;
    }

    toCms.rotateZ(-1*Ptmp.phi());
    toCms.rotateY(-1*Ptmp.theta());

    G4LorentzRotation toLab(toCms.inverse());

    Pprojectile.transform(toCms);
    Ptarget.transform(toCms);

    G4double Pt2;
    G4double ProjMassT2, ProjMassT;
    G4double TargMassT2, TargMassT;
    G4double PZcms2, PZcms;
    G4double PMinusNew, TPlusNew;

    G4double S=Psum.mag2();
    G4double SqrtS=std::sqrt(S);

    if(SqrtS < 2200*MeV) {return false;}   // The model cannot work for pp-interactions
    // at Plab < 1.3 GeV/c. Uzhi

    PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2-
            2*S*Mprojectile2-2*S*Mtarget2-2*Mprojectile2*Mtarget2)/4./S;
    if(PZcms2 < 0)
    {return false;}   // It can be in an interaction with off-shell nuclear nucleon

    PZcms = std::sqrt(PZcms2);

    if(PutOnMassShell)
    {
        if(Pprojectile.z() > 0.)
        {
            Pprojectile.setPz( PZcms);
            Ptarget.setPz(    -PZcms);
        }
        else
        {
            Pprojectile.setPz(-PZcms);
            Ptarget.setPz(     PZcms);
        };

        Pprojectile.setE(std::sqrt(Mprojectile2+
                Pprojectile.x()*Pprojectile.x()+
                Pprojectile.y()*Pprojectile.y()+
                PZcms2));
        Ptarget.setE(std::sqrt(    Mtarget2    +
                Ptarget.x()*Ptarget.x()+
                Ptarget.y()*Ptarget.y()+
                PZcms2));
    }

    G4double maxPtSquare = PZcms2;

    //G4cout << "Pprojectile aft boost : " << Pprojectile << G4endl;
    //G4cout << "Ptarget aft boost : " << Ptarget << G4endl;
    //     G4cout << "cms aft boost : " << (Pprojectile+ Ptarget) << G4endl;

    //     G4cout << " Projectile Xplus / Xminus : " <<
            //      Pprojectile.plus() << " / " << Pprojectile.minus() << G4endl;
    //     G4cout << " Target Xplus / Xminus : " <<
            //      Ptarget.plus() << " / " << Ptarget.minus() << G4endl;

    G4LorentzVector Qmomentum;
    G4double Qminus, Qplus;

    // /* Vova
    G4int whilecount=0;
    do {
        //  Generate pt

        if (whilecount++ >= 500 && (whilecount%100)==0)
            //       G4cout << "G4QGSDiffractiveExcitation::ExciteParticipants possibly looping"
            //       << ", loop count/ maxPtSquare : "
            //               << whilecount << " / " << maxPtSquare << G4endl;
            if (whilecount > 1000 )
            {
                Qmomentum=G4LorentzVector(0.,0.,0.,0.);
                return false;     //  Ignore this interaction
            }

        Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);

        //G4cout << "generated Pt " << Qmomentum << G4endl;
        //G4cout << "Pprojectile with pt : " << Pprojectile+Qmomentum << G4endl;
        //G4cout << "Ptarget with pt : " << Ptarget-Qmomentum << G4endl;

        //  Momentum transfer
        /*                                                                          // Uzhi
           G4double Xmin = minmass / ( Pprojectile.e() + Ptarget.e() );
           G4double Xmax=1.;
           G4double Xplus =ChooseX(Xmin,Xmax);
           G4double Xminus=ChooseX(Xmin,Xmax);

//         G4cout << " X-plus  " << Xplus << G4endl;
//         G4cout << " X-minus " << Xminus << G4endl;

           G4double pt2=G4ThreeVector(Qmomentum.vect()).mag2();
           G4double Qplus =-1 * pt2 / Xminus/Ptarget.minus();
           G4double Qminus=     pt2 / Xplus /Pprojectile.plus();
         */                                                                          // Uzhi    *

        Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
        ProjMassT2=Mprojectile2+Pt2;
        ProjMassT =std::sqrt(ProjMassT2);

        TargMassT2=Mtarget2+Pt2;
        TargMassT =std::sqrt(TargMassT2);

        PZcms2=(S*S+ProjMassT2*ProjMassT2+
                TargMassT2*TargMassT2-
                2.*S*ProjMassT2-2.*S*TargMassT2-
                2.*ProjMassT2*TargMassT2)/4./S;
        if(PZcms2 < 0 ) {PZcms2=0;};
        PZcms =std::sqrt(PZcms2);

        G4double PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
        G4double PMinusMax=SqrtS-TargMassT;

        PMinusNew=ChooseP(PMinusMin,PMinusMax);
        Qminus=PMinusNew-Pprojectile.minus();

        G4double TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
        G4double TPlusMax=SqrtS-ProjMassT;

        TPlusNew=ChooseP(TPlusMin, TPlusMax);
        Qplus=-(TPlusNew-Ptarget.plus());

        Qmomentum.setPz( (Qplus-Qminus)/2 );
        Qmomentum.setE(  (Qplus+Qminus)/2 );

        //G4cout << "Qplus / Qminus " << Qplus << " / " << Qminus<<G4endl;
        //       G4cout << "pt2" << pt2 << G4endl;
        //       G4cout << "Qmomentum " << Qmomentum << G4endl;
        //       G4cout << " Masses (P/T) : " << (Pprojectile+Qmomentum).mag() <<
                //                 " / " << (Ptarget-Qmomentum).mag() << G4endl;
        /*                                                                                 // Uzhi
       } while ( (Pprojectile+Qmomentum).mag2() <= Mprojectile2 ||
             (Ptarget-Qmomentum).mag2()     <= Mtarget2 );
         */                                                                                 // Uzhi     *


    } while (( (Pprojectile+Qmomentum).mag2() <  Mprojectile2       ||      // Uzhi No without excitation
            (Ptarget    -Qmomentum).mag2() <  Mtarget2             ) ||  // Uzhi
            ( (Pprojectile+Qmomentum).mag2()  <  ProjectileDiffCut2 &&     // Uzhi No double Diffraction
                    (Ptarget    -Qmomentum).mag2()  <  NuclearNucleonDiffCut2) );// Uzhi

    if((Ptarget-Qmomentum).mag2() < NuclearNucleonDiffCut2)                 // Uzhi Projectile diffraction
    {
        G4double TMinusNew=SqrtS-PMinusNew;
        Qminus=Ptarget.minus()-TMinusNew;
        TPlusNew=TargMassT2/TMinusNew;
        Qplus=Ptarget.plus()-TPlusNew;

        Qmomentum.setPz( (Qplus-Qminus)/2 );
        Qmomentum.setE(  (Qplus+Qminus)/2 );
    }
    else if((Pprojectile+Qmomentum).mag2() <  ProjectileDiffCut2)    // Uzhi Target diffraction
    {
        G4double PPlusNew=SqrtS-TPlusNew;
        Qplus=PPlusNew-Pprojectile.plus();
        PMinusNew=ProjMassT2/PPlusNew;
        Qminus=PMinusNew-Pprojectile.minus();

        Qmomentum.setPz( (Qplus-Qminus)/2 );
        Qmomentum.setE(  (Qplus+Qminus)/2 );
    };

    Pprojectile += Qmomentum;
    Ptarget     -= Qmomentum;

    // Vova

    /*
        Pprojectile.setPz(0.);
        Pprojectile.setE(SqrtS-M0target);

        Ptarget.setPz(0.);
        Ptarget.setE(M0target);
     */

    //G4cout << "Pprojectile with Q : " << Pprojectile << G4endl;
    //G4cout << "Ptarget with Q : " << Ptarget << G4endl;

    //     G4cout << "Projectile back: " << toLab * Pprojectile << G4endl;
    //     G4cout << "Target back: " << toLab * Ptarget << G4endl;

    // Transform back and update SplitableHadron Participant.
    Pprojectile.transform(toLab);
    Ptarget.transform(toLab);

    //G4cout << "Pprojectile with Q M: " << Pprojectile<<" "<<  Pprojectile.mag() << G4endl;
    //G4cout << "Ptarget     with Q M: " << Ptarget    <<" "<<  Ptarget.mag()     << G4endl;

    //G4cout << "Target  mass  " <<  Ptarget.mag() << G4endl;

    target->Set4Momentum(Ptarget);

    //G4cout << "Projectile mass  " <<  Pprojectile.mag() << G4endl;

    projectile->Set4Momentum(Pprojectile);

    return true;
}

G4ExcitedString * G4QGSDiffractiveExcitation::String ( G4VSplitableHadron *  aHadron, G4bool  isProjectile  ) const

virtual

Definition at line 356 of file G4QGSDiffractiveExcitation.cc.

References G4cout, G4endl, G4VSplitableHadron::Get4Momentum(), G4Parton::Get4Momentum(), G4VSplitableHadron::GetNextParton(), G4Parton::GetPDGcode(), G4VSplitableHadron::GetPosition(), CLHEP::HepLorentzVector::mag(), G4INCL::Math::max(), CLHEP::HepLorentzVector::minus(), CLHEP::HepLorentzVector::perp2(), CLHEP::HepLorentzVector::plus(), CLHEP::HepLorentzVector::px(), CLHEP::HepLorentzVector::py(), G4Parton::Set4Momentum(), CLHEP::HepLorentzVector::setE(), CLHEP::HepLorentzVector::setPx(), CLHEP::HepLorentzVector::setPy(), CLHEP::HepLorentzVector::setPz(), UUtils::Sign(), G4VSplitableHadron::SplitUp(), sqr(), CLHEP::Hep3Vector::x(), and CLHEP::Hep3Vector::y().

{
    hadron->SplitUp();
    G4Parton *start= hadron->GetNextParton();
    if ( start==NULL) 
    { G4cout << " G4FTFModel::String() Error:No start parton found"<< G4endl;
    return NULL;
    }
    G4Parton *end  = hadron->GetNextParton();
    if ( end==NULL) 
    { G4cout << " G4FTFModel::String() Error:No end parton found"<< G4endl;
    return NULL;
    }

    G4ExcitedString * string;
    if ( isProjectile ) 
    {
        string= new G4ExcitedString(end,start, +1);
    } else {
        string= new G4ExcitedString(start,end, -1);
    }

    string->SetPosition(hadron->GetPosition());

    // momenta of string ends
    G4double ptSquared= hadron->Get4Momentum().perp2();
    G4double transverseMassSquared= hadron->Get4Momentum().plus()
                            *   hadron->Get4Momentum().minus();


    G4double maxAvailMomentumSquared=
            sqr( std::sqrt(transverseMassSquared) - std::sqrt(ptSquared) );

    G4double widthOfPtSquare = 0.25;          // Uzhi <Pt^2>=0.25 ??????????????????
    G4ThreeVector pt=GaussianPt(widthOfPtSquare,maxAvailMomentumSquared);

    G4LorentzVector Pstart(G4LorentzVector(pt,0.));
    G4LorentzVector Pend;
    Pend.setPx(hadron->Get4Momentum().px() - pt.x());
    Pend.setPy(hadron->Get4Momentum().py() - pt.y());

    G4double tm1=hadron->Get4Momentum().minus() +
            ( Pend.perp2()-Pstart.perp2() ) / hadron->Get4Momentum().plus();

    G4double tm2= std::sqrt( std::max(0., sqr(tm1) -
            4. * Pend.perp2() * hadron->Get4Momentum().minus()
            /  hadron->Get4Momentum().plus() ));

    G4int Sign= isProjectile ? -1 : 1;

    G4double endMinus  = 0.5 * (tm1 + Sign*tm2);
    G4double startMinus= hadron->Get4Momentum().minus() - endMinus;

    G4double startPlus= Pstart.perp2() /  startMinus;
    G4double endPlus  = hadron->Get4Momentum().plus() - startPlus;

    Pstart.setPz(0.5*(startPlus - startMinus));
    Pstart.setE(0.5*(startPlus + startMinus));

    Pend.setPz(0.5*(endPlus - endMinus));
    Pend.setE(0.5*(endPlus + endMinus));

    start->Set4Momentum(Pstart);
    end->Set4Momentum(Pend);

    #ifdef G4_FTFDEBUG
        G4cout << " generated string flavors          " << start->GetPDGcode() << " / " << end->GetPDGcode() << G4endl;
        G4cout << " generated string momenta:   quark " << start->Get4Momentum() << "mass : " <<start->Get4Momentum().mag()<< G4endl;
        G4cout << " generated string momenta: Diquark " << end ->Get4Momentum() << "mass : " <<end->Get4Momentum().mag()<< G4endl;
        G4cout << " sum of ends                       " << Pstart+Pend << G4endl;
        G4cout << " Original                          " << hadron->Get4Momentum() << G4endl;
    #endif

    return string;  
}

---

The documentation for this class was generated from the following files:

• G4QGSDiffractiveExcitation.hh
• G4QGSDiffractiveExcitation.cc