Changeset 1340 for trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4FTFModel.cc

Timestamp:

Nov 5, 2010, 3:45:55 PM (14 years ago)

Author:

garnier

Message:

update ti head

File:

• trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4FTFModel.cc (modified) (32 diffs)

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

@@ -25 +25 @@
 //
 //
-// $Id: G4FTFModel.cc,v 1.34 2009/12/15 19:14:31 vuzhinsk Exp $
-// GEANT4 tag $Name: geant4-09-04-beta-01 $
+// $Id: G4FTFModel.cc,v 1.36 2010/09/20 15:50:46 vuzhinsk Exp $
+// GEANT4 tag $Name: geant4-09-03-ref-09 $
 //
 
@@ -99 +99 @@
 {
         theProjectile = aProjectile;  
-        theParticipants.Init(aNucleus.GetN(),aNucleus.GetZ()); 
+
+        theParticipants.Init(aNucleus.GetA_asInt(),aNucleus.GetZ_asInt()); 
+
 // ----------- N-mass number Z-charge -------------------------
 
@@ -113 +115 @@
 
 //theParameters->SetProbabilityOfElasticScatt(0.); 
+//G4cout<<theParameters->GetProbabilityOfElasticScatt()<<G4endl;
+//G4int Uzhi; G4cin>>Uzhi;
 // To turn on/off (1/0) elastic scattering
 
@@ -125 +129 @@
 { 
         G4ExcitedStringVector * theStrings(0);
-
+//G4cout<<"GetString"<<G4endl;
         theParticipants.GetList(theProjectile,theParameters);
-
+//G4cout<<"Reggeon"<<G4endl;
         ReggeonCascade(); 
 
@@ -133 +137 @@
         if( PutOnMassShell() )
         {
+//G4cout<<"PutOn mass Shell OK"<<G4endl;
          if( ExciteParticipants() )
          {
+//G4cout<<"Excite partic OK"<<G4endl;
           theStrings = BuildStrings();
-
+//G4cout<<"Build String OK"<<G4endl;
           GetResidualNucleus();
 
@@ -194 +200 @@
            TheInvolvedNucleon[NumberOfInvolvedNucleon]=TargetNucleon;
            NumberOfInvolvedNucleon++;
-
+//G4cout<<"Prim NumberOfInvolvedNucleon "<<NumberOfInvolvedNucleon<<G4endl;
            G4double XofWoundedNucleon = TargetNucleon->GetPosition().x();
            G4double YofWoundedNucleon = TargetNucleon->GetPosition().y();
@@ -214 +220 @@
               TheInvolvedNucleon[NumberOfInvolvedNucleon]=Neighbour;
               NumberOfInvolvedNucleon++;
+//G4cout<<"Seco NumberOfInvolvedNucleon "<<NumberOfInvolvedNucleon<<G4endl;
 
               G4VSplitableHadron *targetSplitable; 
@@ -274 +281 @@
         G4LorentzVector Pprojectile=primary->Get4Momentum();
 
+//G4cout<<"Pprojectile "<<Pprojectile<<G4endl;
 // To get original projectile particle
 
@@ -284 +292 @@
         G4double        SumMasses = Mprojectile + 20.*MeV; // 13.12.09
                                                // Separation energy for projectile
-
+//G4cout<<"SumMasses Pr "<<SumMasses<<G4endl;
 //--------------- Target nucleus ------------------------------
         G4V3DNucleus *theNucleus = GetWoundedNucleus();
@@ -305 +313 @@
           SumMasses += aNucleon->GetDefinition()->GetPDGMass();
           SumMasses += 20.*MeV;   // 13.12.09 Separation energy for a nucleon
+//G4cout<<"SumMasses Tr "<<SumMasses<<G4endl;
           ResidualMassNumber--;
           ResidualCharge-=(G4int) aNucleon->GetDefinition()->GetPDGCharge();
@@ -316 +325 @@
 
         Psum += PnuclearResidual;
-
+//G4cout<<"ResidualCharge ,ResidualMassNumber "<<ResidualCharge<<" "<<ResidualMassNumber<<G4endl;
         G4double ResidualMass(0.);
         if(ResidualMassNumber == 0)
@@ -331 +340 @@
  
 //      ResidualMass +=ResidualExcitationEnergy; // Will be given after checks
+//G4cout<<"SumMasses End ResidualMass "<<SumMasses<<" "<<ResidualMass<<G4endl;
         SumMasses += ResidualMass;
-
+//G4cout<<"SumMasses + ResM "<<SumMasses<<G4endl;
+//G4cout<<"Psum "<<Psum<<G4endl;
 //-------------------------------------------------------------
 
@@ -338 +349 @@
         G4double     S=Psum.mag2();
 
+//G4cout<<"SqrtS < SumMasses "<<SqrtS<<" "<<SumMasses<<G4endl;
         if(SqrtS < SumMasses)      {return false;} // It is impossible to simulate
                                                    // after putting nuclear nucleons
@@ -346 +358 @@
         ResidualMass +=ResidualExcitationEnergy;
         SumMasses    +=ResidualExcitationEnergy;
-
+//G4cout<<"ResidualMass "<<ResidualMass<<" "<<SumMasses<<G4endl;
 //-------------------------------------------------------------
 // Sampling of nucleons what are transfered to delta-isobars --
@@ -373 +385 @@
         }   // end of if(theNucleus.GetMassNumber() != 1)
 //-------------------------------------------------------------
+
         G4LorentzRotation toCms(-1*Psum.boostVector());
         G4LorentzVector Ptmp=toCms*Pprojectile;
@@ -393 +406 @@
         G4double AveragePt2  = theParameters->GetPt2ofNuclearDestruction();
         G4double maxPtSquare = theParameters->GetMaxPt2ofNuclearDestruction();
-
+//G4cout<<"Dcor "<<Dcor<<" AveragePt2 "<<AveragePt2<<G4endl;
         G4double M2target(0.);
         G4double WminusTarget(0.);
@@ -409 +422 @@
 
             NumberOfTries++;
+//G4cout<<"NumberOfTries "<<NumberOfTries<<G4endl;
             if(NumberOfTries == 100*(NumberOfTries/100))   // 100
             { // At large number of tries it would be better to reduce the values
@@ -439 +453 @@
 
                G4LorentzVector tmp(tmpPt.x(),tmpPt.y(),Xminus,0.);
+//G4cout<<"Inv i mom "<<i<<" "<<tmp<<G4endl;
                aNucleon->SetMomentum(tmp);
              }   // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
@@ -447 +462 @@
              G4double DeltaXminus(0.);
 
+//G4cout<<"ResidualMassNumber "<<ResidualMassNumber<<" "<<PtSum<<G4endl;
              if(ResidualMassNumber == 0)
              {
@@ -457 +473 @@
               DeltaXminus = -1./theNucleus->GetMassNumber();
              }
-
+//G4cout<<"Dx y xmin "<<DeltaX<<" "<<DeltaY<<" "<<DeltaXminus<<G4endl;
              XminusSum=1.;
              M2target =0.;
@@ -467 +483 @@
                Xminus = aNucleon->Get4Momentum().pz() - DeltaXminus;
                XminusSum-=Xminus;               
-
-               if((Xminus <= 0.)   || (Xminus > 1.) || 
-                  (XminusSum <=0.) || (XminusSum > 1.)) {InerSuccess=false; break;}
- 
+//G4cout<<" i X-sum "<<i<<" "<<Xminus<<" "<<XminusSum<<G4endl;
+               if(ResidualMassNumber == 0)                // Uzhi 5.07.10
+               {
+                if((Xminus <= 0.)   || (Xminus > 1.))    {InerSuccess=false; break;}
+               } else
+               {
+                if((Xminus <= 0.)   || (Xminus > 1.) || 
+                   (XminusSum <=0.) || (XminusSum > 1.)) {InerSuccess=false; break;}
+               }                                          // Uzhi 5.07.10
+
                G4double Px=aNucleon->Get4Momentum().px() - DeltaX;
                G4double Py=aNucleon->Get4Momentum().py() - DeltaY;
@@ -482 +504 @@
              }   // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
 
+//G4cout<<"Rescale O.K."<<G4endl;
+
              if(InerSuccess && (ResidualMassNumber != 0))
              {
               M2target +=(ResidualMass*ResidualMass + PtSum.mag2())/XminusSum;
              }
+//G4cout<<"InerSuccess "<<InerSuccess<<G4endl;
+//G4int Uzhi;G4cin>>Uzhi;
             } while(!InerSuccess);
           } while (SqrtS < Mprojectile + std::sqrt(M2target));
@@ -495 +521 @@
           WminusTarget=(S-M2projectile+M2target+std::sqrt(DecayMomentum2))/2./SqrtS;
           WplusProjectile=SqrtS - M2target/WminusTarget;
+//G4cout<<"DecayMomentum2 "<<DecayMomentum2<<G4endl;
 //-------------------------------------------------------------
           for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
@@ -511 +538 @@
            if( E+Pz > WplusProjectile ){OuterSuccess=false; break;}
           }   // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
+//G4int Uzhi;G4cin>>Uzhi;
         } while(!OuterSuccess);
 
@@ -575 +603 @@
         Successfull=false;
         theParticipants.StartLoop();
+
+G4int MaxNumOfInelCollisions=G4int(theParameters->GetMaxNumberOfCollisions());
+G4double NumberOfInel(0.);
+//
+if(MaxNumOfInelCollisions > 0)  
+{   //  Plab > Pbound, Normal application of FTF is possible
+ G4double ProbMaxNumber=theParameters->GetMaxNumberOfCollisions()-MaxNumOfInelCollisions;
+ if(G4UniformRand() < ProbMaxNumber) {MaxNumOfInelCollisions++;}
+ NumberOfInel=MaxNumOfInelCollisions;
+} else
+{   //  Plab < Pbound, Normal application of FTF is impossible, low energy corrections
+ if(theParticipants.theNucleus->GetMassNumber() > 1)
+ {
+  NumberOfInel = theParameters->GetProbOfInteraction();
+  MaxNumOfInelCollisions = 1;
+ } else
+ { // Special case for hadron-nucleon interactions
+  NumberOfInel = 1.;
+  MaxNumOfInelCollisions = 1;
+ }
+}  // end of if(MaxNumOfInelCollisions > 0)
+//
         while (theParticipants.Next())
         {          
@@ -581 +631 @@
            G4VSplitableHadron * projectile=collision.GetProjectile();
            G4VSplitableHadron * target=collision.GetTarget();
-
+//G4cout<<"ProbabilityOfElasticScatt "<<theParameters->GetProbabilityOfElasticScatt()<<G4endl;
            if(G4UniformRand()< theParameters->GetProbabilityOfElasticScatt())
            { //   Elastic scattering -------------------------
+//G4cout<<"Elastic FTF"<<G4endl;
             if(theElastic->ElasticScattering(projectile, target, theParameters))
             {
@@ -595 +646 @@
            else
            { //   Inelastic scattering ---------------------- 
+/*
             if(theExcitation->ExciteParticipants(projectile, target, 
                                                  theParameters, theElastic))
@@ -604 +656 @@
              target->SetStatus(2);
             }
+*/
+//G4cout<<"InElastic FTF"<<G4endl;
+            if(G4UniformRand()< NumberOfInel/MaxNumOfInelCollisions)
+            {
+             if(theExcitation->ExciteParticipants(projectile, target, 
+                                                 theParameters, theElastic))
+             {
+              Successfull = Successfull || true; 
+NumberOfInel--;
+             } else
+             {
+              Successfull = Successfull || false;
+              target->SetStatus(2);
+             }
+            } else // If NumOfInel
+            {
+             if(theElastic->ElasticScattering(projectile, target, theParameters))
+             {
+              Successfull = Successfull || true;
+             } else
+             {
+              Successfull = Successfull || false;
+              target->SetStatus(2);
+             }
+            }   // end if NumOfInel
            }
         }       // end of while (theParticipants.Next())
@@ -624 +701 @@
 
         theParticipants.StartLoop();    // restart a loop 
-
+//
         while ( theParticipants.Next() ) 
         {
@@ -634 +711 @@
                 primaries.push_back(interaction.GetProjectile());     
         }
-            
+
         unsigned int ahadron;
         for ( ahadron=0; ahadron < primaries.size() ; ahadron++)
@@ -650 +727 @@
             if(SecondString != 0) strings->push_back(SecondString);
         }
-
+//
         for (G4int ahadron=0; ahadron < NumberOfInvolvedNucleon ; ahadron++)
         {
@@ -669 +746 @@
         std::for_each(primaries.begin(), primaries.end(), DeleteVSplitableHadron());
         primaries.clear();
-
         return strings;
 }