G4QGSDiffractiveExcitation.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id$
00028 // ------------------------------------------------------------
00029 //      GEANT 4 class implemetation file
00030 //
00031 //      ---------------- G4QGSDiffractiveExcitation --------------
00032 //             by Gunter Folger, October 1998.
00033 //         diffractive Excitation used by strings models
00034 //         Take a projectile and a target
00035 //         excite the projectile and target
00036 //  Essential changed by V. Uzhinsky in November - December 2006
00037 //  in order to put it in a correspondence with original FRITIOF
00038 //  model. Variant of FRITIOF with nucleon de-excitation is implemented.  
00039 // ---------------------------------------------------------------------
00040 
00041 // Modified:
00042 //  25-05-07 : G.Folger
00043 //       move from management/G4DiffractiveExcitation to to qgsm/G4QGSDiffractiveExcitation
00044 //                  
00045 
00046 #include "globals.hh"
00047 #include "G4PhysicalConstants.hh"
00048 #include "G4SystemOfUnits.hh"
00049 #include "Randomize.hh"
00050 
00051 #include "G4QGSDiffractiveExcitation.hh"
00052 #include "G4LorentzRotation.hh"
00053 #include "G4ThreeVector.hh"
00054 #include "G4ParticleDefinition.hh"
00055 #include "G4VSplitableHadron.hh"
00056 #include "G4ExcitedString.hh"
00057 //#include "G4ios.hh"
00058 
00059 G4QGSDiffractiveExcitation::G4QGSDiffractiveExcitation()                                     // Uzhi
00060 {
00061 }
00062 
00063 G4QGSDiffractiveExcitation::~G4QGSDiffractiveExcitation()
00064 {
00065 }
00066 
00067 
00068 G4bool G4QGSDiffractiveExcitation::
00069 ExciteParticipants(G4VSplitableHadron *projectile, G4VSplitableHadron *target) const
00070 {
00071 
00072         G4LorentzVector Pprojectile=projectile->Get4Momentum();
00073 
00074         // -------------------- Projectile parameters -----------------------------------
00075         G4bool PutOnMassShell=0;
00076 
00077         //         G4double M0projectile=projectile->GetDefinition()->GetPDGMass();  // With de-excitation
00078         G4double M0projectile = Pprojectile.mag();                        // Without de-excitation
00079 
00080         if(M0projectile < projectile->GetDefinition()->GetPDGMass())
00081         {
00082                 PutOnMassShell=1;
00083                 M0projectile=projectile->GetDefinition()->GetPDGMass();
00084         }
00085 
00086         G4double Mprojectile2 = M0projectile * M0projectile;
00087 
00088         G4int    PDGcode=projectile->GetDefinition()->GetPDGEncoding();
00089         G4int    absPDGcode=std::abs(PDGcode);
00090         G4double ProjectileDiffCut;
00091         G4double AveragePt2;
00092 
00093         if( absPDGcode > 1000 )                        //------Projectile is baryon --------
00094         {
00095                 ProjectileDiffCut = 1.1;                    // GeV
00096                 AveragePt2 = 0.3;                           // GeV^2
00097         }
00098         else if( absPDGcode == 211 || PDGcode ==  111) //------Projectile is Pion -----------
00099         {
00100                 ProjectileDiffCut = 1.0;                   // GeV
00101                 AveragePt2 = 0.3;                          // GeV^2
00102         }
00103         else if( absPDGcode == 321 || PDGcode == -311) //------Projectile is Kaon -----------
00104         {
00105                 ProjectileDiffCut = 1.1;                    // GeV
00106                 AveragePt2 = 0.3;                           // GeV^2
00107         }
00108         else                                           //------Projectile is undefined, Nucleon assumed
00109         {
00110                 ProjectileDiffCut = 1.1;                    // GeV
00111                 AveragePt2 = 0.3;                           // GeV^2
00112         };
00113 
00114         ProjectileDiffCut = ProjectileDiffCut * GeV;
00115         AveragePt2 = AveragePt2 * GeV*GeV;
00116 
00117         // -------------------- Target parameters ----------------------------------------------
00118         G4LorentzVector Ptarget=target->Get4Momentum();
00119 
00120         G4double M0target = Ptarget.mag();
00121 
00122         if(M0target < target->GetDefinition()->GetPDGMass())
00123         {
00124                 PutOnMassShell=1;
00125                 M0target=target->GetDefinition()->GetPDGMass();
00126         }
00127 
00128         G4double Mtarget2 = M0target * M0target;                      //Ptarget.mag2(); // for AA-inter.
00129 
00130         G4double NuclearNucleonDiffCut = 1.1*GeV;
00131 
00132         G4double ProjectileDiffCut2 = ProjectileDiffCut * ProjectileDiffCut;
00133         G4double NuclearNucleonDiffCut2 = NuclearNucleonDiffCut * NuclearNucleonDiffCut;
00134 
00135         // Transform momenta to cms and then rotate parallel to z axis;
00136 
00137         G4LorentzVector Psum;
00138         Psum=Pprojectile+Ptarget;
00139 
00140         G4LorentzRotation toCms(-1*Psum.boostVector());
00141 
00142         G4LorentzVector Ptmp=toCms*Pprojectile;
00143 
00144         if ( Ptmp.pz() <= 0. )
00145         {
00146                 // "String" moving backwards in  CMS, abort collision !!
00147                 //G4cout << " abort Collision!! " << G4endl;
00148                 return false;
00149         }
00150 
00151         toCms.rotateZ(-1*Ptmp.phi());
00152         toCms.rotateY(-1*Ptmp.theta());
00153 
00154         G4LorentzRotation toLab(toCms.inverse());
00155 
00156         Pprojectile.transform(toCms);
00157         Ptarget.transform(toCms);
00158 
00159         G4double Pt2;
00160         G4double ProjMassT2, ProjMassT;
00161         G4double TargMassT2, TargMassT;
00162         G4double PZcms2, PZcms;
00163         G4double PMinusNew, TPlusNew;
00164 
00165         G4double S=Psum.mag2();
00166         G4double SqrtS=std::sqrt(S);
00167 
00168         if(SqrtS < 2200*MeV) {return false;}   // The model cannot work for pp-interactions
00169         // at Plab < 1.3 GeV/c. Uzhi
00170 
00171         PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2-
00172                         2*S*Mprojectile2-2*S*Mtarget2-2*Mprojectile2*Mtarget2)/4./S;
00173         if(PZcms2 < 0)
00174         {return false;}   // It can be in an interaction with off-shell nuclear nucleon
00175 
00176         PZcms = std::sqrt(PZcms2);
00177 
00178         if(PutOnMassShell)
00179         {
00180                 if(Pprojectile.z() > 0.)
00181                 {
00182                         Pprojectile.setPz( PZcms);
00183                         Ptarget.setPz(    -PZcms);
00184                 }
00185                 else
00186                 {
00187                         Pprojectile.setPz(-PZcms);
00188                         Ptarget.setPz(     PZcms);
00189                 };
00190 
00191                 Pprojectile.setE(std::sqrt(Mprojectile2+
00192                                 Pprojectile.x()*Pprojectile.x()+
00193                                 Pprojectile.y()*Pprojectile.y()+
00194                                 PZcms2));
00195                 Ptarget.setE(std::sqrt(    Mtarget2    +
00196                                 Ptarget.x()*Ptarget.x()+
00197                                 Ptarget.y()*Ptarget.y()+
00198                                 PZcms2));
00199         }
00200 
00201         G4double maxPtSquare = PZcms2;
00202 
00203         //G4cout << "Pprojectile aft boost : " << Pprojectile << G4endl;
00204         //G4cout << "Ptarget aft boost : " << Ptarget << G4endl;
00205         //         G4cout << "cms aft boost : " << (Pprojectile+ Ptarget) << G4endl;
00206 
00207         //         G4cout << " Projectile Xplus / Xminus : " <<
00208                         //              Pprojectile.plus() << " / " << Pprojectile.minus() << G4endl;
00209         //         G4cout << " Target Xplus / Xminus : " <<
00210                         //              Ptarget.plus() << " / " << Ptarget.minus() << G4endl;
00211 
00212         G4LorentzVector Qmomentum;
00213         G4double Qminus, Qplus;
00214 
00215         // /* Vova
00216         G4int whilecount=0;
00217         do {
00218                 //  Generate pt
00219 
00220                 if (whilecount++ >= 500 && (whilecount%100)==0)
00221                         //               G4cout << "G4QGSDiffractiveExcitation::ExciteParticipants possibly looping"
00222                         //               << ", loop count/ maxPtSquare : "
00223                         //               << whilecount << " / " << maxPtSquare << G4endl;
00224                         if (whilecount > 1000 )
00225                         {
00226                                 Qmomentum=G4LorentzVector(0.,0.,0.,0.);
00227                                 return false;     //  Ignore this interaction
00228                         }
00229 
00230                 Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
00231 
00232                 //G4cout << "generated Pt " << Qmomentum << G4endl;
00233                 //G4cout << "Pprojectile with pt : " << Pprojectile+Qmomentum << G4endl;
00234                 //G4cout << "Ptarget with pt : " << Ptarget-Qmomentum << G4endl;
00235 
00236                 //  Momentum transfer
00237                 /*                                                                          // Uzhi
00238                G4double Xmin = minmass / ( Pprojectile.e() + Ptarget.e() );
00239                G4double Xmax=1.;
00240                G4double Xplus =ChooseX(Xmin,Xmax);
00241                G4double Xminus=ChooseX(Xmin,Xmax);
00242 
00243 //             G4cout << " X-plus  " << Xplus << G4endl;
00244 //             G4cout << " X-minus " << Xminus << G4endl;
00245 
00246                G4double pt2=G4ThreeVector(Qmomentum.vect()).mag2();
00247                G4double Qplus =-1 * pt2 / Xminus/Ptarget.minus();
00248                G4double Qminus=     pt2 / Xplus /Pprojectile.plus();
00249                  */                                                                          // Uzhi    *
00250 
00251                 Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
00252                 ProjMassT2=Mprojectile2+Pt2;
00253                 ProjMassT =std::sqrt(ProjMassT2);
00254 
00255                 TargMassT2=Mtarget2+Pt2;
00256                 TargMassT =std::sqrt(TargMassT2);
00257 
00258                 PZcms2=(S*S+ProjMassT2*ProjMassT2+
00259                                 TargMassT2*TargMassT2-
00260                                 2.*S*ProjMassT2-2.*S*TargMassT2-
00261                                 2.*ProjMassT2*TargMassT2)/4./S;
00262                 if(PZcms2 < 0 ) {PZcms2=0;};
00263                 PZcms =std::sqrt(PZcms2);
00264 
00265                 G4double PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
00266                 G4double PMinusMax=SqrtS-TargMassT;
00267 
00268                 PMinusNew=ChooseP(PMinusMin,PMinusMax);
00269                 Qminus=PMinusNew-Pprojectile.minus();
00270 
00271                 G4double TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
00272                 G4double TPlusMax=SqrtS-ProjMassT;
00273 
00274                 TPlusNew=ChooseP(TPlusMin, TPlusMax);
00275                 Qplus=-(TPlusNew-Ptarget.plus());
00276 
00277                 Qmomentum.setPz( (Qplus-Qminus)/2 );
00278                 Qmomentum.setE(  (Qplus+Qminus)/2 );
00279 
00280                 //G4cout << "Qplus / Qminus " << Qplus << " / " << Qminus<<G4endl;
00281                 //           G4cout << "pt2" << pt2 << G4endl;
00282                 //           G4cout << "Qmomentum " << Qmomentum << G4endl;
00283                 //           G4cout << " Masses (P/T) : " << (Pprojectile+Qmomentum).mag() <<
00284                                 //                             " / " << (Ptarget-Qmomentum).mag() << G4endl;
00285                 /*                                                                                 // Uzhi
00286            } while ( (Pprojectile+Qmomentum).mag2() <= Mprojectile2 ||
00287                      (Ptarget-Qmomentum).mag2()     <= Mtarget2 );
00288                  */                                                                                 // Uzhi     *
00289 
00290 
00291         } while (( (Pprojectile+Qmomentum).mag2() <  Mprojectile2       ||      // Uzhi No without excitation
00292                         (Ptarget    -Qmomentum).mag2() <  Mtarget2             ) ||  // Uzhi
00293                         ( (Pprojectile+Qmomentum).mag2()  <  ProjectileDiffCut2 &&     // Uzhi No double Diffraction
00294                                         (Ptarget    -Qmomentum).mag2()  <  NuclearNucleonDiffCut2) );// Uzhi
00295 
00296         if((Ptarget-Qmomentum).mag2() < NuclearNucleonDiffCut2)                 // Uzhi Projectile diffraction
00297         {
00298                 G4double TMinusNew=SqrtS-PMinusNew;
00299                 Qminus=Ptarget.minus()-TMinusNew;
00300                 TPlusNew=TargMassT2/TMinusNew;
00301                 Qplus=Ptarget.plus()-TPlusNew;
00302 
00303                 Qmomentum.setPz( (Qplus-Qminus)/2 );
00304                 Qmomentum.setE(  (Qplus+Qminus)/2 );
00305         }
00306         else if((Pprojectile+Qmomentum).mag2() <  ProjectileDiffCut2)    // Uzhi Target diffraction
00307         {
00308                 G4double PPlusNew=SqrtS-TPlusNew;
00309                 Qplus=PPlusNew-Pprojectile.plus();
00310                 PMinusNew=ProjMassT2/PPlusNew;
00311                 Qminus=PMinusNew-Pprojectile.minus();
00312 
00313                 Qmomentum.setPz( (Qplus-Qminus)/2 );
00314                 Qmomentum.setE(  (Qplus+Qminus)/2 );
00315         };
00316 
00317         Pprojectile += Qmomentum;
00318         Ptarget     -= Qmomentum;
00319 
00320         // Vova
00321 
00322         /*
00323                 Pprojectile.setPz(0.);
00324                 Pprojectile.setE(SqrtS-M0target);
00325 
00326                 Ptarget.setPz(0.);
00327                 Ptarget.setE(M0target);
00328          */
00329 
00330         //G4cout << "Pprojectile with Q : " << Pprojectile << G4endl;
00331         //G4cout << "Ptarget with Q : " << Ptarget << G4endl;
00332 
00333         //         G4cout << "Projectile back: " << toLab * Pprojectile << G4endl;
00334         //         G4cout << "Target back: " << toLab * Ptarget << G4endl;
00335 
00336         // Transform back and update SplitableHadron Participant.
00337         Pprojectile.transform(toLab);
00338         Ptarget.transform(toLab);
00339 
00340         //G4cout << "Pprojectile with Q M: " << Pprojectile<<" "<<  Pprojectile.mag() << G4endl;
00341         //G4cout << "Ptarget     with Q M: " << Ptarget    <<" "<<  Ptarget.mag()     << G4endl;
00342 
00343         //G4cout << "Target      mass  " <<  Ptarget.mag() << G4endl;
00344 
00345         target->Set4Momentum(Ptarget);
00346 
00347         //G4cout << "Projectile mass  " <<  Pprojectile.mag() << G4endl;
00348 
00349         projectile->Set4Momentum(Pprojectile);
00350 
00351         return true;
00352 }
00353 
00354 
00355 G4ExcitedString * G4QGSDiffractiveExcitation::
00356 String(G4VSplitableHadron * hadron, G4bool isProjectile) const
00357 {
00358         hadron->SplitUp();
00359         G4Parton *start= hadron->GetNextParton();
00360         if ( start==NULL) 
00361         { G4cout << " G4FTFModel::String() Error:No start parton found"<< G4endl;
00362         return NULL;
00363         }
00364         G4Parton *end  = hadron->GetNextParton();
00365         if ( end==NULL) 
00366         { G4cout << " G4FTFModel::String() Error:No end parton found"<< G4endl;
00367         return NULL;
00368         }
00369 
00370         G4ExcitedString * string;
00371         if ( isProjectile ) 
00372         {
00373                 string= new G4ExcitedString(end,start, +1);
00374         } else {
00375                 string= new G4ExcitedString(start,end, -1);
00376         }
00377 
00378         string->SetPosition(hadron->GetPosition());
00379 
00380         // momenta of string ends
00381         G4double ptSquared= hadron->Get4Momentum().perp2();
00382         G4double transverseMassSquared= hadron->Get4Momentum().plus()
00383                                                 *       hadron->Get4Momentum().minus();
00384 
00385 
00386         G4double maxAvailMomentumSquared=
00387                         sqr( std::sqrt(transverseMassSquared) - std::sqrt(ptSquared) );
00388 
00389         G4double widthOfPtSquare = 0.25;          // Uzhi <Pt^2>=0.25 ??????????????????
00390         G4ThreeVector pt=GaussianPt(widthOfPtSquare,maxAvailMomentumSquared);
00391 
00392         G4LorentzVector Pstart(G4LorentzVector(pt,0.));
00393         G4LorentzVector Pend;
00394         Pend.setPx(hadron->Get4Momentum().px() - pt.x());
00395         Pend.setPy(hadron->Get4Momentum().py() - pt.y());
00396 
00397         G4double tm1=hadron->Get4Momentum().minus() +
00398                         ( Pend.perp2()-Pstart.perp2() ) / hadron->Get4Momentum().plus();
00399 
00400         G4double tm2= std::sqrt( std::max(0., sqr(tm1) -
00401                         4. * Pend.perp2() * hadron->Get4Momentum().minus()
00402                         /  hadron->Get4Momentum().plus() ));
00403 
00404         G4int Sign= isProjectile ? -1 : 1;
00405 
00406         G4double endMinus  = 0.5 * (tm1 + Sign*tm2);
00407         G4double startMinus= hadron->Get4Momentum().minus() - endMinus;
00408 
00409         G4double startPlus= Pstart.perp2() /  startMinus;
00410         G4double endPlus  = hadron->Get4Momentum().plus() - startPlus;
00411 
00412         Pstart.setPz(0.5*(startPlus - startMinus));
00413         Pstart.setE(0.5*(startPlus + startMinus));
00414 
00415         Pend.setPz(0.5*(endPlus - endMinus));
00416         Pend.setE(0.5*(endPlus + endMinus));
00417 
00418         start->Set4Momentum(Pstart);
00419         end->Set4Momentum(Pend);
00420 
00421         #ifdef G4_FTFDEBUG
00422                 G4cout << " generated string flavors          " << start->GetPDGcode() << " / " << end->GetPDGcode() << G4endl;
00423                 G4cout << " generated string momenta:   quark " << start->Get4Momentum() << "mass : " <<start->Get4Momentum().mag()<< G4endl;
00424                 G4cout << " generated string momenta: Diquark " << end ->Get4Momentum() << "mass : " <<end->Get4Momentum().mag()<< G4endl;
00425                 G4cout << " sum of ends                       " << Pstart+Pend << G4endl;
00426                 G4cout << " Original                          " << hadron->Get4Momentum() << G4endl;
00427         #endif
00428 
00429         return string;  
00430 }
00431 
00432 
00433 // --------- private methods ----------------------
00434 
00435 G4double G4QGSDiffractiveExcitation::ChooseP(G4double Pmin, G4double Pmax) const // Uzhi
00436 {
00437         // choose an x between Xmin and Xmax with P(x) ~ 1/x
00438 
00439         //  to be improved...
00440 
00441         G4double range=Pmax-Pmin;                                         // Uzhi
00442 
00443         if ( Pmin <= 0. || range <=0. ) 
00444         {
00445                 G4cout << " Pmin, range : " << Pmin << " , " << range << G4endl;
00446                 throw G4HadronicException(__FILE__, __LINE__, "G4QGSDiffractiveExcitation::ChooseP : Invalid arguments ");
00447         }
00448 
00449         G4double P;
00450         /*                                                                          // Uzhi
00451         do {
00452             x=Xmin + G4UniformRand() * range;
00453         }  while ( Xmin/x < G4UniformRand() );
00454          */                                                                          // Uzhi
00455 
00456         P=Pmin * std::pow(Pmax/Pmin,G4UniformRand());                       // Uzhi
00457 
00458         //debug-hpw     cout << "DiffractiveX "<<x<<G4endl;
00459         return P;
00460 }
00461 
00462 G4ThreeVector G4QGSDiffractiveExcitation::GaussianPt(G4double AveragePt2, G4double maxPtSquare) const // Uzhi
00463 {            //  @@ this method is used in FTFModel as well. Should go somewhere common!
00464 
00465         G4double Pt2;
00466         /*                                                                          // Uzhi
00467         do {
00468             pt2=widthSquare * std::log( G4UniformRand() );
00469         } while ( pt2 > maxPtSquare);
00470          */                                                                          // Uzhi
00471 
00472         Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() * (std::exp(-maxPtSquare/AveragePt2)-1.));// Uzhi
00473 
00474         G4double Pt=std::sqrt(Pt2);
00475 
00476         G4double phi=G4UniformRand() * twopi;
00477 
00478         return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);    
00479 }

Generated on Mon May 27 17:49:34 2013 for Geant4 by  doxygen 1.4.7