G4ElasticHNScattering.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 //      ---------------- G4ElasticHNScattering --------------
00032 //                   by V. Uzhinsky, March 2008.
00033 //             elastic scattering used by Fritiof model
00034 //                 Take a projectile and a target
00035 //                 scatter the projectile and target
00036 // ---------------------------------------------------------------------
00037 
00038 
00039 #include "globals.hh"
00040 #include "Randomize.hh"
00041 #include "G4PhysicalConstants.hh"
00042 
00043 #include "G4ElasticHNScattering.hh"
00044 #include "G4LorentzRotation.hh"
00045 #include "G4ThreeVector.hh"
00046 #include "G4ParticleDefinition.hh"
00047 #include "G4VSplitableHadron.hh"
00048 #include "G4ExcitedString.hh"
00049 #include "G4FTFParameters.hh"
00050 //#include "G4ios.hh"
00051 
00052 G4ElasticHNScattering::G4ElasticHNScattering()          
00053 {
00054 }
00055 
00056 G4bool G4ElasticHNScattering::
00057            ElasticScattering (G4VSplitableHadron *projectile, 
00058                               G4VSplitableHadron *target,
00059                               G4FTFParameters    *theParameters) const
00060 {
00061 // -------------------- Projectile parameters -----------------------------------
00062            G4LorentzVector Pprojectile=projectile->Get4Momentum();
00063 
00064            if(Pprojectile.z() < 0.)
00065            {
00066             target->SetStatus(2);
00067             return false;
00068            } 
00069 
00070            G4bool PutOnMassShell(false);
00071 
00072            G4double M0projectile = Pprojectile.mag();                        
00073            if(M0projectile < projectile->GetDefinition()->GetPDGMass()) 
00074            {
00075               PutOnMassShell=true;
00076               M0projectile=projectile->GetDefinition()->GetPDGMass();
00077            }
00078 
00079            G4double Mprojectile2 = M0projectile * M0projectile;
00080 
00081            G4double AveragePt2=theParameters->GetAvaragePt2ofElasticScattering();
00082 
00083 // -------------------- Target parameters ----------------------------------------------
00084 
00085            G4LorentzVector Ptarget=target->Get4Momentum();
00086 
00087            G4double M0target = Ptarget.mag();
00088 
00089            if(M0target < target->GetDefinition()->GetPDGMass()) 
00090            {
00091               PutOnMassShell=true;
00092               M0target=target->GetDefinition()->GetPDGMass();
00093            }
00094      
00095            G4double Mtarget2 = M0target * M0target;                      
00096 
00097 // Transform momenta to cms and then rotate parallel to z axis;
00098 
00099            G4LorentzVector Psum;
00100            Psum=Pprojectile+Ptarget;
00101 
00102            G4LorentzRotation toCms(-1*Psum.boostVector());
00103 
00104            G4LorentzVector Ptmp=toCms*Pprojectile;
00105 
00106            if ( Ptmp.pz() <= 0. )                                
00107            {
00108            // "String" moving backwards in  CMS, abort collision !!
00109            //G4cout << " abort Collision!! " << G4endl;
00110                    target->SetStatus(2);
00111                    return false; 
00112            }
00113                            
00114            toCms.rotateZ(-1*Ptmp.phi());
00115            toCms.rotateY(-1*Ptmp.theta());
00116         
00117            G4LorentzRotation toLab(toCms.inverse());
00118 
00119            Pprojectile.transform(toCms);
00120            Ptarget.transform(toCms);
00121 
00122 // ---------------------- Putting on mass-on-shell, if needed ------------------------
00123            G4double PZcms2, PZcms;                                          
00124 
00125            G4double S=Psum.mag2();                                          
00126 //         G4double SqrtS=std::sqrt(S);                                     
00127 
00128            PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2-
00129                                  2*S*Mprojectile2-2*S*Mtarget2-2*Mprojectile2*Mtarget2)/4./S;
00130 
00131            if(PZcms2 < 0.)
00132            {  // It can be in an interaction with off-shell nuclear nucleon
00133             if(M0projectile > projectile->GetDefinition()->GetPDGMass()) 
00134             {  // An attempt to de-excite the projectile
00135                // It is assumed that the target is in the ground state
00136              M0projectile = projectile->GetDefinition()->GetPDGMass();
00137              Mprojectile2=M0projectile*M0projectile;
00138              PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2-
00139                     2*S*Mprojectile2 - 2*S*Mtarget2 - 2*Mprojectile2*Mtarget2)
00140                     /4./S;
00141 
00142              if(PZcms2 < 0.){ return false;} // Non succesful attempt after the de-excitation
00143             }
00144             else // if(M0projectile > projectile->GetDefinition()->GetPDGMass())
00145             {
00146              target->SetStatus(2);                                   
00147              return false;                   // The projectile was not excited,
00148                                              // but the energy was too low to put
00149                                              // the target nucleon on mass-shell
00150             }   // end of if(M0projectile > projectile->GetDefinition()->GetPDGMass())
00151            }    // end of if(PZcms2 < 0.)
00152 
00153            PZcms = std::sqrt(PZcms2);
00154 
00155            if(PutOnMassShell)
00156            {
00157               if(Pprojectile.z() > 0.)
00158               {
00159                  Pprojectile.setPz( PZcms);
00160                  Ptarget.setPz(    -PZcms);
00161               }
00162               else  // if(Pprojectile.z() > 0.)
00163               {
00164                  Pprojectile.setPz(-PZcms);
00165                  Ptarget.setPz(     PZcms);
00166               };
00167 
00168               Pprojectile.setE(std::sqrt(Mprojectile2+
00169                                                       Pprojectile.x()*Pprojectile.x()+
00170                                                       Pprojectile.y()*Pprojectile.y()+
00171                                                       PZcms2));
00172               Ptarget.setE(std::sqrt(    Mtarget2    +
00173                                                       Ptarget.x()*Ptarget.x()+
00174                                                       Ptarget.y()*Ptarget.y()+
00175                                                       PZcms2));
00176            }  // end of if(PutOnMassShell)
00177 
00178            G4double maxPtSquare = PZcms2;
00179 
00180 // ------ Now we can calculate the transfered Pt --------------------------
00181            G4double Pt2;                                                    
00182            G4double ProjMassT2; //, ProjMassT;                                  
00183            G4double TargMassT2; //, TargMassT;
00184 
00185            G4LorentzVector Qmomentum;
00186            Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
00187 
00188            Pt2=G4ThreeVector(Qmomentum.vect()).mag2();                  
00189 
00190            ProjMassT2=Mprojectile2+Pt2;                           
00191 //           ProjMassT =std::sqrt(ProjMassT2);                            
00192 
00193            TargMassT2=Mtarget2+Pt2;                               
00194 //           TargMassT =std::sqrt(TargMassT2);                            
00195 
00196            PZcms2=(S*S+ProjMassT2*ProjMassT2+                           
00197                        TargMassT2*TargMassT2-                           
00198                     2.*S*ProjMassT2-2.*S*TargMassT2-                 
00199                     2.*ProjMassT2*TargMassT2)/4./S;                  
00200            if(PZcms2 < 0 ) {PZcms2=0;};// to avoid the exactness problem
00201            PZcms =std::sqrt(PZcms2);                                    
00202 
00203            Pprojectile.setPz( PZcms);  
00204            Ptarget.setPz(    -PZcms); 
00205 
00206            Pprojectile += Qmomentum;
00207            Ptarget     -= Qmomentum;
00208 
00209 // Transform back and update SplitableHadron Participant.
00210            Pprojectile.transform(toLab);
00211            Ptarget.transform(toLab);
00212 /*  // Maybe it will be needed for an exact calculations--------------------
00213            G4double TargetMomentum=std::sqrt(Ptarget.x()*Ptarget.x()+
00214                                              Ptarget.y()*Ptarget.y()+
00215                                              Ptarget.z()*Ptarget.z());
00216 */
00217 
00218 // Calculation of the creation time ---------------------
00219       projectile->SetTimeOfCreation(target->GetTimeOfCreation());
00220       projectile->SetPosition(target->GetPosition());
00221 // Creation time and position of target nucleon were determined at
00222 // ReggeonCascade() of G4FTFModel
00223 // ------------------------------------------------------
00224 
00225            projectile->Set4Momentum(Pprojectile);
00226            target->Set4Momentum(Ptarget);
00227 
00228            projectile->IncrementCollisionCount(1);
00229            target->IncrementCollisionCount(1);
00230 
00231            return true;
00232 }
00233 
00234 
00235 // --------- private methods ----------------------
00236 
00237 G4ThreeVector G4ElasticHNScattering::GaussianPt(G4double AveragePt2, G4double maxPtSquare) const
00238 {            //  @@ this method is used in FTFModel as well. Should go somewhere common!
00239         
00240         G4double Pt2(0.);
00241         if(AveragePt2 <= 0.) {Pt2=0.;}
00242         else
00243         {
00244          Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() * 
00245                 (std::exp(-maxPtSquare/AveragePt2)-1.)); 
00246         }
00247         G4double Pt=std::sqrt(Pt2);
00248         G4double phi=G4UniformRand() * twopi;
00249         
00250         return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);    
00251 }
00252 
00253 G4ElasticHNScattering::G4ElasticHNScattering(const G4ElasticHNScattering &)
00254 {
00255         throw G4HadronicException(__FILE__, __LINE__, "G4ElasticHNScattering copy contructor not meant to be called");
00256 }
00257 
00258 
00259 G4ElasticHNScattering::~G4ElasticHNScattering()
00260 {
00261 }
00262 
00263 
00264 const G4ElasticHNScattering & G4ElasticHNScattering::operator=(const G4ElasticHNScattering &)
00265 {
00266         throw G4HadronicException(__FILE__, __LINE__, "G4ElasticHNScattering = operator not meant to be called");
00267         //return *this;  //A.R. 25-Jul-2012 : fix Coverity
00268 }
00269 
00270 
00271 int G4ElasticHNScattering::operator==(const G4ElasticHNScattering &) const
00272 {
00273         throw G4HadronicException(__FILE__, __LINE__, "G4ElasticHNScattering == operator not meant to be called");
00274         //return false;  //A.R. 25-Jul-2012 : fix Coverity
00275 }
00276 
00277 int G4ElasticHNScattering::operator!=(const G4ElasticHNScattering &) const
00278 {
00279         throw G4HadronicException(__FILE__, __LINE__, "G4ElasticHNScattering != operator not meant to be called");
00280         //return true;  //A.R. 25-Jul-2012 : fix Coverity
00281 }

Generated on Mon May 27 17:48:06 2013 for Geant4 by  doxygen 1.4.7