Changeset 1055 for trunk/source/processes/hadronic/models/chiral_inv_phase_space/interface/src/G4QStringChipsParticleLevelInterface.cc

Timestamp:

May 28, 2009, 4:26:57 PM (15 years ago)

Author:

garnier

Message:

maj sur la beta de geant 4.9.3

File:

• trunk/source/processes/hadronic/models/chiral_inv_phase_space/interface/src/G4QStringChipsParticleLevelInterface.cc (modified) (19 diffs)

trunk/source/processes/hadronic/models/chiral_inv_phase_space/interface/src/G4QStringChipsParticleLevelInterface.cc

@@ -23 +23 @@
 // * acceptance of all terms of the Geant4 Software license.          *
 // ********************************************************************
-//
+// Short description: Interface of QGSC to CHIPS (Energy Flow  of soft hadrons) 
+// ----------------------------------------------------------------------------
 
 //#define debug
@@ -143 +144 @@
   }
   G4int targetPDGCode = 90000000 + 1000*resZ + (resA-resZ);    // PDG of theResidualNucleus
-  G4double targetMass = theNucleus->GetMass();                 // Its mass
-  targetMass -= hitMass; // subtract masses of knocked out nucleons (binding?! M.K.) E/M
+  G4double targetMass=mNeut;
+  if (!resZ)                                                   // Nucleus of only neutrons
+  {
+    if (resA>1) targetMass*=resA;
+  }
+  else targetMass=G4ParticleTable::GetParticleTable()->FindIon(resZ,resA,0,resZ)
+                                                                            ->GetPDGMass();
   G4double targetEnergy = std::sqrt(hitMomentum.mag2()+targetMass*targetMass);
   // !! @@ Target should be at rest: hitMomentum=(0,0,0) @@ !! M.K. (go to this system)
@@ -153 +159 @@
   G4double impactX = theImpact.first;
   G4double impactY = theImpact.second;
-  G4double inpactPar2 = impactX*impactX + impactY*impactY;
+  G4double impactPar2 = impactX*impactX + impactY*impactY;
   
   G4double radius2 = theNucleus->GetNuclearRadius(theInnerCoreDensityCut*perCent);
@@ -160 +166 @@
 #ifdef pdebug
   G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: r="<<std::sqrt(radius2)/fermi
-        <<", b="<<std::sqrt(inpactPar2)/fermi<<", R="<<theNucleus->GetOuterRadius()/fermi
-        <<", b/r="<<std::sqrt(inpactPar2/radius2)<<G4endl; 
-#endif
-  if(radius2 - inpactPar2>0) pathlength = 2.*std::sqrt(radius2 - inpactPar2);
+        <<", b="<<std::sqrt(impactPar2)/fermi<<", R="<<theNucleus->GetOuterRadius()/fermi
+        <<", b/r="<<std::sqrt(impactPar2/radius2)<<G4endl; 
+#endif
+  if(radius2 - impactPar2>0) pathlength = 2.*std::sqrt(radius2 - impactPar2);
   G4double theEnergyLostInFragmentation = theEnergyLossPerFermi*pathlength/fermi;
   
@@ -176 +182 @@
           << theSecondaries->operator[](secondary)->GetDefinition()->GetPDGCharge()<<" "
           << theSecondaries->operator[](secondary)->GetDefinition()->GetPDGEncoding()<<" "
-          << a4Mom <<G4endl; 
+          << a4Mom <<G4endl; 
 #endif
 #ifdef pdebug
@@ -193 +199 @@
       if((*current).first > toSort)        // The current is smaller then existing
       {
-               theSorted.insert(current, it);     // It shifts the others up
-               inserted = true;
-               break;
+        theSorted.insert(current, it);     // It shifts the others up
+        inserted = true;
+        break;
       }
     }
@@ -234 +240 @@
         <<theSorted.size()<<G4endl; 
 #endif
-  
+  G4bool EscapeExists = false;
   for(current = theSorted.begin(); current!=theSorted.end(); current++)
   {
 #ifdef pdebug
-                                G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: nq="
+    G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: nq="
           <<(*current).second->GetDefinition()->GetQuarkContent(3)<<", naq="
           <<(*current).second->GetDefinition()->GetAntiQuarkContent(3)<<", PDG="
@@ -306 +312 @@
 
 #ifdef pdebug
-                                G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: E="<<runningEnergy<<", EL="
+    G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: E="<<runningEnergy<<", EL="
                                                     <<theEnergyLostInFragmentation<<G4endl;
 #endif
-
-    if(runningEnergy > theEnergyLostInFragmentation) break;
-    
+    if(runningEnergy > theEnergyLostInFragmentation)
+    {
+      EscapeExists = true;
+      break;
+    }
 #ifdef CHIPSdebug
     G4cout <<"G4QStringChipsParticleLevelInterface::Propagate: ABSORBED STRING particles "
            <<(*current).second->GetDefinition()->GetPDGCharge()<<" "
            << (*current).second->GetDefinition()->GetPDGEncoding()<<" "
-                  << (*current).second->Get4Momentum() <<G4endl; 
+           << (*current).second->Get4Momentum() <<G4endl; 
 #endif
 #ifdef pdebug
@@ -322 +330 @@
           <<current->second->GetDefinition()->GetPDGCharge()<<", PDG="
           <<current->second->GetDefinition()->GetPDGEncoding()<<", 4M="
-                 <<current->second->Get4Momentum()<<G4endl; 
+          <<current->second->Get4Momentum()<<G4endl; 
 #endif
 
@@ -382 +390 @@
       {
         theFinalContents[hp] +=theContents[running];
-               *(theFinalMomenta[hp])+=*(theMomenta[running]);
-               running++;
-               if(running == theContents.size()) break;
+        *(theFinalMomenta[hp])+=*(theMomenta[running]);
+        running++;
+        if(running == theContents.size()) break;
       }
     }
@@ -402 +410 @@
   G4QNucleus::SetParameters(fractionOfSingleQuasiFreeNucleons,
                             fractionOfPairedQuasiFreeNucleons,
-                                                 clusteringCoefficient, 
-                                                               fusionToExchange);
+                            clusteringCoefficient, 
+                            fusionToExchange);
   G4Quasmon::SetParameters(temperature, halfTheStrangenessOfSee, etaToEtaPrime);
 
@@ -424 +432 @@
     G4QCHIPSWorld::Get()->GetParticles(nop);
     G4QEnvironment* pan= new G4QEnvironment(projHV, targetPDGCode);
+#ifdef pdebug
+      G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: CHIPS fragmentation, rA="
+            <<resA<<", #AbsPt="<<particleCount<<G4endl; 
+#endif
     try
     {
@@ -438 +450 @@
       {
         G4cerr <<"  Incoming 4-momentum and PDG code of "<<i<<"'th hadron: "
-               <<" "<< projHV[i]->Get4Momentum()<<" "<<projHV[i]->GetPDGCode()<<G4endl;
+        <<" "<< projHV[i]->Get4Momentum()<<" "<<projHV[i]->GetPDGCode()<<G4endl;
       }
       throw;
@@ -447 +459 @@
     delete pan;
   }
-  else output = new G4QHadronVector;
-   
+  else
+  {
+#ifdef pdebug
+    G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: NO CHIPS fragmentation, rA="
+          <<resA<<", #AbsPt="<<particleCount<<G4endl; 
+#endif
+    output = new G4QHadronVector;
+  }   
   // Fill the result.
 #ifdef CHIPSdebug
-  G4cout << "NEXT EVENT"<<endl;
+  G4cout << "NEXT EVENT, EscapeExists="<<EscapeExists<<G4endl;
 #endif
 
   // first decay and add all escaping particles.
-  for(current = firstEscape; current!=theSorted.end(); current++)
+  if (EscapeExists) for (current = firstEscape; current!=theSorted.end(); current++)
   {
     G4KineticTrack* aResult = (*current).second;
@@ -472 +490 @@
       theSec->SetMomentum(current4Mom.vect());
 #ifdef pdebug
-                                  G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: *OUT* QGS stable PDG="
+      G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: *OUT* QGS stable PDG="
             <<aResult->GetDefinition()->GetPDGEncoding()<<",4M="<<current4Mom<<G4endl; 
 #endif
@@ -487 +505 @@
         theSec->SetMomentum(current4Mom.vect());
 #ifdef pdebug
-                                    G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: *OUT* QGS decay PDG="
+        G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: *OUT* QGS decay PDG="
               <<secondaries->operator[](aSecondary)->GetDefinition()->GetPDGEncoding()
               <<",4M="<<current4Mom<<G4endl; 
@@ -561 +579 @@
       theSec->SetMomentum(current4Mom.vect());
 #ifdef pdebug
-                                  G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: *OUT* CHIPS PDG="
+      G4cout<<"G4QStringChipsParticleLevelInterface::Propagate: *OUT* CHIPS PDG="
               <<theDefinition->GetPDGEncoding()<<",4M="<<current4Mom<<G4endl; 
 #endif
@@ -576 +594 @@
     G4cout <<"CHIPS particles "<<theDefinition->GetPDGCharge()<<" "
            << theDefinition->GetPDGEncoding()<<" "
-           << output->operator[](particle)->Get4Momentum() <<G4endl; 
+    << output->operator[](particle)->Get4Momentum() <<G4endl; 
 #endif
 
     delete output->operator[](particle);
+  }
+  else
+  {
+    if(resA>0)
+    {
+      G4ParticleDefinition* theDefinition = G4Neutron::Neutron();
+      if(resA==1) // The residual nucleus at rest must be added to conserve BaryN & Charge
+      {
+        if(resZ == 1) theDefinition = G4Proton::Proton();
+      }
+      else theDefinition = G4ParticleTable::GetParticleTable()->FindIon(resZ,resA,0,resZ);
+      theSec = new G4ReactionProduct(theDefinition);
+      theSec->SetTotalEnergy(theDefinition->GetPDGMass());
+      theSec->SetMomentum(G4ThreeVector(0.,0.,0.));
+      theResult->push_back(theSec);
+      if(!resZ && resA>0) for(G4int ni=1; ni<resA; ni++) 
+      {
+        theSec = new G4ReactionProduct(theDefinition);
+        theSec->SetTotalEnergy(theDefinition->GetPDGMass());
+        theSec->SetMomentum(G4ThreeVector(0.,0.,0.));
+        theResult->push_back(theSec);
+      }
+    }
   }
   delete output;
@@ -588 +629 @@
   G4cout << "QUASMON preparation info "
          << 1./MeV*proj4Mom<<" "
-         << 1./MeV*targ4Mom<<" "
-         << nD<<" "<<nU<<" "<<nS<<" "<<nAD<<" "<<nAU<<" "<<nAS<<" "
-         << hitCount<<" "
-         << particleCount<<" "
-         << theLow<<" "
-         << theHigh<<" "
-         << G4endl;
+  << 1./MeV*targ4Mom<<" "
+  << nD<<" "<<nU<<" "<<nS<<" "<<nAD<<" "<<nAU<<" "<<nAS<<" "
+  << hitCount<<" "
+  << particleCount<<" "
+  << theLow<<" "
+  << theHigh<<" "
+  << G4endl;
 #endif