Changeset 1340 for trunk/source/processes/hadronic/models/parton_string/hadronization/src/G4ExcitedStringDecay.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/hadronization/src/G4ExcitedStringDecay.cc (modified) (4 diffs)

trunk/source/processes/hadronic/models/parton_string/hadronization/src/G4ExcitedStringDecay.cc

@@ -61 +61 @@
 }
 
-G4bool G4ExcitedStringDecay::
-EnergyAndMomentumCorrector(G4KineticTrackVector* Output, G4LorentzVector& TotalCollisionMom)   
+G4KineticTrackVector *G4ExcitedStringDecay::FragmentString
+                                (const G4ExcitedString &theString)
+{
+        if ( theStringDecay == NULL ) 
+
+            theStringDecay=new G4LundStringFragmentation();
+    
+        return theStringDecay->FragmentString(theString);
+}
+        
+G4KineticTrackVector *G4ExcitedStringDecay::FragmentStrings
+                                (const G4ExcitedStringVector * theStrings)
+{
+  G4LorentzVector KTsum(0.,0.,0.,0.);
+  for ( unsigned int astring=0; astring < theStrings->size(); astring++)
+  {
+        KTsum+= theStrings->operator[](astring)->Get4Momentum();
+
+        if( !(KTsum.e()<1) && !(KTsum.e()>-1) )
+        {
+          throw G4HadronicException(__FILE__, __LINE__, 
+                                   "G4ExcitedStringDecay::FragmentStrings received nan string...");
+        }
+  }
+
+  G4KineticTrackVector * theResult = new G4KineticTrackVector;
+  G4int attempts(0);
+  G4bool success=false;
+  do {
+
+        std::for_each(theResult->begin() , theResult->end() , DeleteKineticTrack());
+        theResult->clear();
+
+        attempts++;
+        G4LorentzVector KTsecondaries(0.,0.,0.,0.);
+        G4bool NeedEnergyCorrector=false;
+
+        for ( unsigned int astring=0; astring < theStrings->size(); astring++)
+        {
+//G4cout<<G4endl<<"String No "<<astring+1<<" "<<theStrings->operator[](astring)->Get4Momentum().mag()<<G4endl;
+
+          G4KineticTrackVector * generatedKineticTracks = NULL;
+  
+          if ( theStrings->operator[](astring)->IsExcited() )
+          {
+             generatedKineticTracks=FragmentString(*theStrings->operator[](astring));
+          } else {
+             generatedKineticTracks = new G4KineticTrackVector;
+             generatedKineticTracks->push_back(theStrings->operator[](astring)->GetKineticTrack());
+          }    
+
+          if (generatedKineticTracks == NULL) 
+          {
+             G4cerr << "G4VPartonStringModel:No KineticTracks produced" << G4endl;
+             continue;
+          }
+
+          G4LorentzVector KTsum1(0.,0.,0.,0.);
+          for ( unsigned int aTrack=0; aTrack<generatedKineticTracks->size();aTrack++)
+          {
+//G4cout<<" Prod part "<<(*generatedKineticTracks)[aTrack]->GetDefinition()->GetParticleName()<<" "<<(*generatedKineticTracks)[aTrack]->Get4Momentum()<<G4endl;
+             theResult->push_back(generatedKineticTracks->operator[](aTrack));
+             KTsum1+= (*generatedKineticTracks)[aTrack]->Get4Momentum();
+          }
+          KTsecondaries+=KTsum1;
+        
+          if  ( KTsum1.e() > 0 && std::abs((KTsum1.e()-theStrings->operator[](astring)->Get4Momentum().e()) / KTsum1.e()) > perMillion ) 
+          {
+//--debug--           G4cout << "String secondaries(" <<generatedKineticTracks->size()<< ")  momentum: " 
+//--debug--               << theStrings->operator[](astring)->Get4Momentum() << " " << KTsum1 << G4endl;
+            NeedEnergyCorrector=true;
+          }
+
+//        clean up
+          delete generatedKineticTracks;
+        }
+//--DEBUG  G4cout << "Strings/secs total  4 momentum " << KTsum << " " <<KTsecondaries << G4endl;
+
+        success=true;
+        if ( NeedEnergyCorrector ) success=EnergyAndMomentumCorrector(theResult, KTsum);
+
+  } while(!success && (attempts < 100));
+
+#ifdef debug_ExcitedStringDecay
+  G4LorentzVector  KTsum1=0;
+  for ( unsigned int aTrack=0; aTrack<theResult->size();aTrack++)
+  {
+      G4cout << " corrected tracks .. " << (*theResult)[aTrack]->GetDefinition()->GetParticleName()
+      <<"  " << (*theResult)[aTrack]->Get4Momentum() << G4endl;
+      KTsum1+= (*theResult)[aTrack]->Get4Momentum();
+  }
+  G4cout << "Needcorrector/success " << NeedEnergyCorrector << "/" << success << ", Corrected total  4 momentum " << KTsum1  << G4endl;
+  if ( ! success ) G4cout << "failed to correct E/p" << G4endl;  
+#endif
+
+  return theResult;
+}
+
+G4bool G4ExcitedStringDecay::EnergyAndMomentumCorrector
+                (G4KineticTrackVector* Output, G4LorentzVector& TotalCollisionMom)   
   {
     const int    nAttemptScale = 500;
@@ -71 +169 @@
     G4double        SumMass = 0;     
     G4double        TotalCollisionMass = TotalCollisionMom.m();
+
     if( !(TotalCollisionMass<1) && !(TotalCollisionMass>-1) )
     {
@@ -77 +176 @@
     }
 
+//G4cout<<G4endl<<"EnergyAndMomentumCorrector "<<Output->size()<<G4endl;
     // Calculate sum hadron 4-momenta and summing hadron mass
     unsigned int cHadron;
@@ -92 +192 @@
         }
     }
+
+//G4cout<<"SumMass TotalCollisionMass "<<SumMass<<" "<<TotalCollisionMass<<G4endl;
 
     // Cannot correct a single particle