source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/body/src/G4QFragmentation.cc @ 967

Last change on this file since 967 was 962, checked in by garnier, 15 years ago

update processes
File size: 31.2 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer                                           *
4// *                                                                  *
5// * The  Geant4 software  is  copyright of the Copyright Holders  of *
6// * the Geant4 Collaboration.  It is provided  under  the terms  and *
7// * conditions of the Geant4 Software License,  included in the file *
8// * LICENSE and available at  http://cern.ch/geant4/license .  These *
9// * include a list of copyright holders.                             *
10// *                                                                  *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work  make  any representation or  warranty, express or implied, *
14// * regarding  this  software system or assume any liability for its *
15// * use.  Please see the license in the file  LICENSE  and URL above *
16// * for the full disclaimer and the limitation of liability.         *
17// *                                                                  *
18// * This  code  implementation is the result of  the  scientific and *
19// * technical work of the GEANT4 collaboration.                      *
20// * By using,  copying,  modifying or  distributing the software (or *
21// * any work based  on the software)  you  agree  to acknowledge its *
22// * use  in  resulting  scientific  publications,  and indicate your *
23// * acceptance of all terms of the Geant4 Software license.          *
24// ********************************************************************
25//
26//
27// $Id: G4QFragmentation.cc,v 1.3 2007/05/02 14:59:55 gunter Exp $
28// GEANT4 tag $Name: geant4-09-02-ref-02 $
29//
30// -----------------------------------------------------------------------------
31//      GEANT 4 class header file
32//
33//                 History:
34//     Created by Mikhail Kossov, October 2006
35//     CHIPS QGS fragmentation class
36//     For comparison mirror member functions are taken from G4 classes:
37//     G4QGSParticipants
38//     G4QGSModels
39//     G4ExcitedStringDecay
40// -----------------------------------------------------------------------------
41//
42
43//#define debug
44//#define pdebug
45//#define ppdebug
46
47#include "globals.hh"
48#include "G4QFragmentation.hh"
49#include "G4LorentzVector.hh"
50#include <utility>
51
52// Promoting model parameters from local variables class properties @@(? M.K.)
53
54// Definition of static parameters
55G4int    G4QFragmentation::nCutMax=7; 
56G4double G4QFragmentation::ThersholdParameter=450.*MeV; 
57G4double G4QFragmentation::QGSMThershold=3.*GeV; 
58G4double G4QFragmentation::theNucleonRadius=1.5*fermi;
59 // Parameters of diffractional fragmentation
60G4double G4QFragmentation::widthOfPtSquare=-0.72*GeV*GeV; // pt -width2 forStringExcitation
61G4double G4QFragmentation::minExtraMass=250.*MeV;// minimum excitation mass
62G4double G4QFragmentation::minmass=250.*MeV;        // mean pion transverse mass for Xmin
63
64G4QFragmentation::G4QFragmentation() : theNucleus(0)
65{
66  // Construct Shortlived particles (needed after the 2006 Particles revolution)
67         G4ShortLivedConstructor ShortLived;
68         ShortLived.ConstructParticle();
69}
70
71G4QFragmentation::~G4QFragmentation() {if(theNucleus) delete theNucleus;}
72
73G4QHadronVector* G4QFragmentation::Scatter(const G4QNucleus &aNucleus,
74                                           const G4QHadron &aPrimary)
75{ 
76  G4QStringVector* strings=0;
77
78  G4QHadron thePrimary = aPrimary;
79 
80  G4LorentzRotation toZ;
81  G4LorentzVector Ptmp=thePrimary.Get4Momentum();
82  toZ.rotateZ(-1*Ptmp.phi());
83  toZ.rotateY(-1*Ptmp.theta());
84  thePrimary.Set4Momentum(toZ*Ptmp);
85  G4LorentzRotation toLab(toZ.inverse());
86
87  G4int attempts = 0, maxAttempts=20;
88  while(!strings)
89  {
90         if (attempts++ > maxAttempts ) 
91         {
92                                        G4cerr<<"***G4QFragmentation::Scatter: "<<attempts<<" to create a string ( > max="
93            <<maxAttempts<<") --> try to increase maxAttempts"<<G4endl;
94      G4Exception("G4QFragmentation::Scatter:","72",FatalException,"StringCreation");
95         }
96           InitModel(aNucleus,thePrimary);
97         strings = GetStrings();
98  }
99 
100  G4QHadronVector* theResult = 0;
101  G4double stringEnergy(0);
102  for( unsigned astring=0; astring < strings->size(); astring++)
103  {
104    //    rotate string to lab frame, models have it aligned to z
105    stringEnergy += (*strings)[astring]->GetLeftParton()->Get4Momentum().t();
106    stringEnergy += (*strings)[astring]->GetRightParton()->Get4Momentum().t();
107    (*strings)[astring]->LorentzRotate(toLab);
108  }
109#ifdef debug
110  G4cout<<"G4QFragmentation::Scatter: Total string energy = "<<stringEnergy<<G4endl;
111#endif
112  theResult = FragmentStrings(strings);
113  std::for_each(strings->begin(), strings->end(), DeleteQString() );
114  delete strings;
115
116  return theResult;
117}
118
119void G4QFragmentation::InitModel(const G4QNucleus& aNucleus,const G4QHadron& aProjectile)
120{
121  static const G4double         mProt = G4Proton::Proton()->GetPDGMass(); // Mass of proton
122  Init(aNucleus.GetN(),aNucleus.GetZ());
123  theCurrentVelocity.setX(0);   
124  theCurrentVelocity.setY(0); 
125  // @@ "target Nucleon" == "Proton at rest" (M.K. ?)
126  G4double nCons = 1;                         // 1 or baryon number of the projectile
127  G4int projAbsB=std::abs(aProjectile.GetBaryonNumber()); // Baryon/anti-baryon
128  if(projAbsB>1) nCons = projAbsB;
129  G4LorentzVector proj4M = aProjectile.Get4Momentum();
130  G4double pz_per_projectile = proj4M.pz()/nCons;
131  G4double e_per_projectile = proj4M.e()/nCons+mProt; // @@ Add mass of TargetProtonAtRest
132  //G4double e_per_projectile = aProjectile.Get4Momentum()*aProjectile.Get4Momentum();
133  //e_per_projectile += proj4M.vect()*proj4M.vect();
134  //e_per_projectile /=nCons*nCons;
135         //e_per_projectile = std::sqrt(e_per_projectile);
136         //e_per_projectile += G4Proton::Proton()->GetPDGMass();//@@Add mass of TargetProtonAtRest
137  G4double vz = pz_per_projectile/e_per_projectile;
138#ifdef debug
139                G4cout<<"G4QFragmentation::Init: Projectile4M="<<proj4M<<", vz="<<vz<<G4endl;
140#endif
141  theCurrentVelocity.setZ(vz);
142  DoLorentzBoost(-theCurrentVelocity); 
143  G4LorentzVector Mom = proj4M;
144  Mom.boost(-theCurrentVelocity);
145  G4QHadron theProjectile(aProjectile.GetQPDG(),Mom);
146#ifdef debug
147  G4cout<<"G4QFragmentation::Init: PreInteractionMomentum"<<Mom<<G4endl;
148#endif
149  BuildInteractions(theProjectile);
150  GetWoundedNucleus()->DoLorentzBoost(theCurrentVelocity);
151}
152
153G4QStringVector* G4QFragmentation::GetStrings()
154{
155  G4QPartonPair* aPair;
156  G4QStringVector* theStrings = new G4QStringVector;
157  G4QString* aString=0;
158  while((aPair = GetNextPartonPair())) // @@ At present no difference in stringBuild ? M.K.
159  {
160    if (aPair->GetCollisionType() == G4QPartonPair::DIFFRACTIVE)
161    {
162      aString = BuildString(aPair);           // @@ ?
163#ifdef debug
164      G4cout<<"G4QFragmentation::GetStrings:DifString4M="<<aString->Get4Momentum()<<G4endl;
165#endif
166    }
167    else
168    {
169      aString = BuildString(aPair);           // @@ ?
170#ifdef debug
171      G4cout<<"G4QFragmentation::GetStrings:SftString4M="<<aString->Get4Momentum()<<G4endl;
172#endif
173    }
174    aString->Boost(theCurrentVelocity); 
175    theStrings->push_back(aString);
176    delete aPair;
177  }
178#ifdef debug
179  for(G4int i=0; i<theStrings->size(); i++) G4cout<<"G4QFragmentation::GetStrings:String #"
180                                    <<i<<", 4M="<<(*theStrings)[i]->Get4Momentum()<<G4endl;
181#endif
182  return theStrings;
183}
184
185void G4QFragmentation::BuildInteractions(const G4QHadron &thePrimary) 
186{
187 
188  // Find the collisions and collition conditions
189  G4QHadron* aProjectile = SelectInteractions(thePrimary);
190
191  // now build the parton pairs
192  SplitHadrons();
193 
194  // soft collisions, ordering is vital
195  PerformSoftCollisions();
196   
197  // the rest is diffractive
198  PerformDiffractiveCollisions();
199 
200  // clean-up, if necessary
201  std::for_each(theInteractions.begin(),theInteractions.end(), DeleteQInteraction());
202  theInteractions.clear();
203  std::for_each(theTargets.begin(), theTargets.end(), DeleteQHadron());
204  theTargets.clear();
205  delete aProjectile;
206}
207
208//
209G4QHadron* G4QFragmentation::SelectInteractions(const G4QHadron &thePrimary)
210{
211  G4QHadron* aProjectile = new G4QHadron(thePrimary);
212  G4QPomeron theProbability(thePrimary.GetPDGCode()); // must be data member
213  G4double outerRadius = theNucleus->GetOuterRadius();
214
215  // Check reaction threshold
216  theNucleus->StartLoop();
217  G4QHadron* pNucleon = theNucleus->GetNextNucleon();
218  G4LorentzVector aPrimaryMomentum=thePrimary.Get4Momentum();
219  G4double s = (aPrimaryMomentum + pNucleon->Get4Momentum()).mag2();
220  G4double primaryMass=thePrimary.GetMass();
221  G4double ThresholdMass = primaryMass + pNucleon->GetMass(); 
222  ModelMode = SOFT;
223  if (sqr(ThresholdMass + ThersholdParameter) > s)
224  {
225    G4cerr<<"***G4QFragmentation::SelectInteractions: ThrM="<<ThresholdMass<<" + ThrPa="
226          <<ThersholdParameter<<" = "<<ThresholdMass+ThersholdParameter<<" > std::sqrt(s)="
227          <<std::sqrt(s)<<G4endl;
228    G4Exception("G4QFragmentation::SelectInteractions:","72",FatalException,"LowEnergy");
229  }
230  if (sqr(ThresholdMass + QGSMThershold) > s) // thus only diffractive in cascade!
231  {
232#ifdef debug
233    G4cout<<"G4QFragmentation::SelectInteractions: ThrM="<<ThresholdMass<<" + ThrQGS="
234          <<QGSMThershold<<" = "<<ThresholdMass+QGSMThershold<<" > std::sqrt(s)="
235          <<std::sqrt(s)<<" -> only Diffraction is possible"<<G4endl; // @@ to Quasmon
236#endif
237    ModelMode = DIFFRACTIVE;
238  }
239 
240  // first find the collisions HPW
241  std::for_each(theInteractions.begin(),theInteractions.end(), DeleteQInteraction());
242  theInteractions.clear();
243  G4int totalCuts = 0;
244  G4double impactUsed = 0;
245
246  G4double eKin = aPrimaryMomentum.e()-primaryMass; // Primary kinetic energy ? GeV ? M.K.
247
248  while(theInteractions.size() == 0)
249  {
250    // choose random impact parameter HPW
251    std::pair<G4double, G4double> theImpactParameter;
252    theImpactParameter = theNucleus->ChooseImpactXandY(outerRadius+theNucleonRadius);
253    G4double impactX = theImpactParameter.first; 
254    G4double impactY = theImpactParameter.second;
255   
256    // loop over nuclei to find collissions HPW
257    theNucleus->StartLoop();
258    G4int nucleonCount = 0; // debug
259    //G4QFragmentation_NPart = 0;                        // ? M.K.
260    while( (pNucleon = theNucleus->GetNextNucleon()) )
261    {
262      if(totalCuts>1.5*eKin) break;
263      nucleonCount++; // debug
264      // Needs to be moved to Probability class @@@
265      G4double s = (aPrimaryMomentum + pNucleon->Get4Momentum()).mag2();
266      G4double Distance2 = sqr(impactX - pNucleon->GetPosition().x()) +
267                           sqr(impactY - pNucleon->GetPosition().y());
268      G4double Probability = theProbability.GetInelasticProbability(s, Distance2); 
269      // test for inelastic collision
270      G4double rndNumber = G4UniformRand();
271      // ModelMode = DIFFRACTIVE;
272      if (Probability > rndNumber)
273      {
274#ifdef debug
275        G4cout<<"G4QFragmentation::SelectInteractions: p="<<Probability<<", r="<<rndNumber
276              <<", d="<<std::sqrt(Distance2)<<G4endl;
277#endif
278        G4QHadron* aTarget = new G4QHadron(*pNucleon);
279        //G4QFragmentation_NPart ++;                       // ? M.K.
280               theTargets.push_back(aTarget);
281              pNucleon=aTarget;
282        if((theProbability.GetDiffractiveProbability(s,Distance2)/Probability >
283            G4UniformRand() && (ModelMode==SOFT)) || ModelMode==DIFFRACTIVE)
284               { 
285                 // diffractive interaction occurs
286                 if(IsSingleDiffractive()) ExciteSingDiffParticipants(aProjectile, aTarget);
287                 else                          ExciteDiffParticipants(aProjectile, aTarget);
288          G4QInteraction* aInteraction = new G4QInteraction(aProjectile);
289          aInteraction->SetTarget(aTarget); 
290          theInteractions.push_back(aInteraction);
291                 aInteraction->SetNumberOfDiffractiveCollisions(1);
292          totalCuts += 1;
293               }
294               else
295               {
296                 // nondiffractive soft interaction occurs
297                 // sample nCut+1 (cut Pomerons) pairs of strings can be produced
298          G4int nCut;
299          G4double* running = new G4double[nCutMax];
300          running[0] = 0;
301          for(nCut = 0; nCut < nCutMax; nCut++)
302          {
303                   running[nCut] = theProbability.GetCutPomeronProbability(s, Distance2, nCut+1);
304            if(nCut!=0) running[nCut] += running[nCut-1];
305          }
306                 G4double random = running[nCutMax-1]*G4UniformRand();
307          for(nCut = 0; nCut < nCutMax; nCut++) {if(running[nCut] > random) break;}
308          delete [] running;
309          nCut = 0;
310          aTarget->IncrementCollisionCount(nCut+1);
311          aProjectile->IncrementCollisionCount(nCut+1);
312          G4QInteraction* aInteraction = new G4QInteraction(aProjectile);
313          aInteraction->SetTarget(aTarget);
314          aInteraction->SetNumberOfSoftCollisions(nCut+1);
315          theInteractions.push_back(aInteraction);
316          totalCuts += nCut+1;
317          impactUsed=Distance2;
318        }
319      }
320    }
321#ifdef debug
322    G4cout<<"G4QFragmentation::SelectInteractions: NUCLEONCOUNT="<<nucleonCount<<G4endl;
323#endif
324  }
325#ifdef debug
326  G4cout<<"G4QFragmentation::SelectInteractions: CUTDEBUG="<<totalCuts
327        <<", ImpactParameter="<<impactUsed<<G4endl;
328#endif
329  return aProjectile;
330}
331
332void G4QFragmentation::PerformDiffractiveCollisions()
333{
334  for(unsigned i = 0; i < theInteractions.size(); i++) 
335  {
336    G4QInteraction* anIniteraction = theInteractions[i];
337    G4QHadron* aProjectile = anIniteraction->GetProjectile();
338    G4QParton* aParton = aProjectile->GetNextParton();
339    G4QPartonPair* aPartonPair;
340    // projectile first
341    if (aParton)
342    {
343      aPartonPair = new G4QPartonPair(aParton, aProjectile->GetNextAntiParton(), 
344                                      G4QPartonPair::DIFFRACTIVE, 
345                                      G4QPartonPair::PROJECTILE);
346      thePartonPairs.push_back(aPartonPair);
347    }
348    // then target
349    G4QHadron* aTarget = anIniteraction->GetTarget();
350    aParton = aTarget->GetNextParton();
351    if (aParton)
352    {
353      aPartonPair = new G4QPartonPair(aParton, aTarget->GetNextAntiParton(), 
354                                      G4QPartonPair::DIFFRACTIVE, 
355                                      G4QPartonPair::TARGET);
356      thePartonPairs.push_back(aPartonPair);
357    }
358  }
359}
360
361void G4QFragmentation::PerformSoftCollisions()
362{
363  G4QInteractionVector::iterator i;
364  for(i = theInteractions.begin(); i != theInteractions.end(); i++)   
365  {
366    G4QInteraction* anIniteraction = *i;
367    G4QPartonPair* aPair=0;
368    if (anIniteraction->GetNumberOfSoftCollisions())
369    { 
370      G4QHadron* pProjectile = anIniteraction->GetProjectile();
371      G4QHadron* pTarget     = anIniteraction->GetTarget();
372      for (G4int j = 0; j < anIniteraction->GetNumberOfSoftCollisions(); j++)
373      {
374        aPair= new G4QPartonPair(pTarget->GetNextParton(),pProjectile->GetNextAntiParton(),
375                                 G4QPartonPair::SOFT, G4QPartonPair::TARGET);
376        thePartonPairs.push_back(aPair);
377        aPair= new G4QPartonPair(pProjectile->GetNextParton(),pTarget->GetNextAntiParton(),
378                                G4QPartonPair::SOFT, G4QPartonPair::PROJECTILE);
379        thePartonPairs.push_back(aPair);
380      } 
381      delete *i;
382      i=theInteractions.erase(i);
383      i--;
384    }
385  }
386}
387
388G4QPartonPair* G4QFragmentation::GetNextPartonPair()
389{
390  if(thePartonPairs.empty()) return 0;
391  G4QPartonPair* result = thePartonPairs.back();
392  thePartonPairs.pop_back();
393  return result;
394}
395
396G4QHadronVector* G4QFragmentation::FragmentStrings(const G4QStringVector* theStrings)
397{
398  G4QHadronVector* theResult = new G4QHadronVector;
399
400  G4LorentzVector KTsum(0.,0.,0.);
401  G4bool NeedEnergyCorrector=false;
402 
403  for( unsigned astring=0; astring < theStrings->size(); astring++)
404  {
405    KTsum+= (*theStrings)[astring]->Get4Momentum();
406    if(!(KTsum.e()<1.) && !(KTsum.e()>-1.))
407    {
408      G4cerr<<"***G4QFragmentation::FragmentStrings: KTsum="<<KTsum<<G4endl;
409      G4Exception("G4QFragmentation::FragmentStrings:","72",FatalException,"NANin3Vector");
410    }
411    G4QHadronVector* generatedHadrons = 0;     // A prototype of the string output
412    if( (*theStrings)[astring]->IsExcited() )
413    {
414      generatedHadrons=(*theStrings)[astring]->FragmentString(true);// Fragment QGSM String
415    }
416    else
417    {
418             //generatedHadrons = new G4QHadronVector;
419             //generatedHadrons->push_back((*theStrings)[astring]->GetAsQHadron()); //@@ NotImplem
420    }   
421           if (!generatedHadrons) 
422           {
423                    G4cerr<<"G4QFragmentation::FragmentStrings: No Hadrons produced" << G4endl;
424                    continue;
425           }
426    G4LorentzVector KTsum1(0.,0.,0.,0.);
427    for(unsigned aTrack=0; aTrack<generatedHadrons->size(); aTrack++)
428           {
429                    theResult->push_back((*generatedHadrons)[aTrack]);
430                    KTsum1+= (*generatedHadrons)[aTrack]->Get4Momentum();
431           }
432       
433           if(std::abs((KTsum1.e()-(*theStrings)[astring]->Get4Momentum().e())/KTsum1.e())
434       > perMillion) NeedEnergyCorrector=true;
435    //      clean up
436           delete generatedHadrons;
437  }
438#ifdef debug
439  G4cout<<"G4QFragmentation::FragmentStrings: String4mom="<<KTsum<<endl; 
440#endif
441  if(NeedEnergyCorrector) EnergyAndMomentumCorrector(theResult, KTsum);
442  return theResult;
443}
444
445G4bool G4QFragmentation::EnergyAndMomentumCorrector(G4QHadronVector* Output,
446                                                    G4LorentzVector& TotalCollisionMom)   
447{
448  const int    nAttemptScale = 500;
449  const double ErrLimit = 1.E-5;
450  if (Output->empty()) return TRUE;
451  G4LorentzVector SumMom;
452  G4double        SumMass = 0;     
453  G4double        TotalCollisionMass = TotalCollisionMom.m();
454  if( !(TotalCollisionMass<1) && !(TotalCollisionMass>-1) )
455  {
456    G4cerr<<"***G4QFragmentation::EnergyAndMomentumCorrect:M="<<TotalCollisionMass<<G4endl;
457    G4Exception("G4QFragmentation::EnergyAndMomentumCorr:","72",FatalException,"NAN_totM");
458  }
459  // Calculate sum hadron 4-momenta and summing hadron mass
460  unsigned int cHadron;
461  for(cHadron = 0; cHadron < Output->size(); cHadron++)
462  {
463    SumMom  += Output->operator[](cHadron)->Get4Momentum();
464    if( !(SumMom<1) && !(SumMom>-1) )
465    {
466      G4cerr<<"***G4QFragmentation::EnergyAndMomentumCorrector: SumMom="<<SumMom<<G4endl;
467      G4Exception("G4QFragmentation::EnergyAndMomentumCorr:","72",FatalException,"NANMom");
468    }
469    SumMass += (*Output)[cHadron]->GetMass();
470    if(!(SumMass<1) && !(SumMass>-1))
471    {
472      G4cerr<<"***G4QFragmentation::EnergyAndMomentumCorrector: SumMass="<<SumMass<<G4endl;
473      G4Exception("G4QFragmentation::EnergyAndMomentumCor:","72",FatalException,"NANMass");
474           }
475  }
476  // Cannot correct a single particle
477  if(Output->size() < 2) return FALSE;
478  if (SumMass > TotalCollisionMass) return FALSE;
479  SumMass = SumMom.m2();
480  if (SumMass < 0) return FALSE;
481  SumMass = std::sqrt(SumMass);
482  // Compute c.m.s. hadron velocity and boost KTV to hadron c.m.s.
483  G4ThreeVector Beta = -SumMom.boostVector();
484  G4int nOut=Output->size();
485  if(nOut) for(G4int o=0; o<nOut; o++) (*Output)[o]->Boost(Beta);
486  // Scale total c.m.s. hadron energy (hadron system mass).
487  // It should be equal interaction mass
488  G4double Scale = 1;
489  G4int cAttempt = 0;
490  G4double Sum = 0;
491  G4bool success = false;
492  for(cAttempt = 0; cAttempt < nAttemptScale; cAttempt++)
493  {
494    Sum = 0;
495    for(cHadron = 0; cHadron < Output->size(); cHadron++)
496    {
497      G4LorentzVector HadronMom = Output->operator[](cHadron)->Get4Momentum();
498      HadronMom.setVect(Scale*HadronMom.vect());
499      G4double E = std::sqrt(HadronMom.vect().mag2() + sqr((*Output)[cHadron]->GetMass()));
500      HadronMom.setE(E);
501      Output->operator[](cHadron)->Set4Momentum(HadronMom);
502      Sum += E;
503    }   
504    Scale = TotalCollisionMass/Sum;   
505    if (Scale - 1 <= ErrLimit) 
506    {
507      success = true;
508      break;
509    }
510#ifdef debug
511    G4cout<<"G4QFragmentation::EnergyAndMomentumCorrector: Scale-1="<<Scale-1<<", TotM="
512          <<TotalCollisionMass<<", Sum="<<Sum<<G4endl;
513#endif     
514  }   
515  if(!success)
516  {
517    G4cout<<"***G4QFragmentation::EnergyAndMomentumCorrector: Scale #1 at end of loop M="
518          <<TotalCollisionMass<<", S"<<Sum<<", Sc="<<Scale
519          <<" Increase number of attempts or increase ERRLIMIT"<<G4endl;
520    G4Exception("G4QFragmentation::SelectInteractions:","72",FatalException,"NotCorrect");
521  }
522  // Compute c.m.s. interaction velocity and KTV back boost
523  Beta = TotalCollisionMom.boostVector(); 
524  nOut=Output->size();
525  if(nOut) for(G4int o=0; o<nOut; o++) (*Output)[o]->Boost(Beta);
526  return TRUE;
527}
528
529// Excite double diffractive string
530G4bool G4QFragmentation::ExciteDiffParticipants(G4QHadron* projectile,
531                                                G4QHadron* target) const
532{
533         G4LorentzVector Pprojectile=projectile->Get4Momentum();
534         G4double Mprojectile=projectile->GetMass() + minExtraMass;
535         G4double Mprojectile2=Mprojectile*Mprojectile;
536  G4LorentzVector Ptarget=target->Get4Momentum();
537  G4double Mtarget=target->GetMass() + minExtraMass;
538  G4double Mtarget2=Mtarget*Mtarget;
539#ifdef debug
540         G4cout<<"G4QFragm::ExciteDiffPartici:Ep="<<Pprojectile.e()<<",Et="<<Ptarget.e()<<G4endl;
541#endif
542  // Transform momenta to cms and then rotate parallel to z axis;
543         G4LorentzVector Psum=Pprojectile+Ptarget;
544         G4LorentzRotation toCms(-Psum.boostVector()); // Boost Rotation to CMS
545         G4LorentzVector Ptmp=toCms*Pprojectile;
546         if(Ptmp.pz()<=0.) // "String" moving backwards in CMS, abort collision !! ? M.K.
547         {
548#ifdef debug
549           G4cout<<"G4QFragmentation::ExciteDiffParticipants: *1* abort Collision!! *1*"<<G4endl;
550#endif
551                  return false; 
552         }                         
553         toCms.rotateZ(-Ptmp.phi());
554         toCms.rotateY(-Ptmp.theta());
555#ifdef debug
556  G4cout<<"G4QFragment::ExciteDiffParticipantts: Be4Boost Pproj="<<Pprojectile<<", Ptarg="
557        <<Ptarget<<G4endl;
558#endif
559         G4LorentzRotation toLab(toCms.inverse()); // Boost Rotation to LabSys (LS)
560         Pprojectile.transform(toCms);
561         Ptarget.transform(toCms);
562#ifdef debug
563         G4cout<< "G4QFragment::ExciteDiffParticipantts: AfterBoost Pproj="<<Pprojectile<<"Ptarg="
564        <<Ptarget<<", cms4M="<<Pprojectile+Ptarget<<G4endl;
565         G4cout<<"G4QFragment::ExciteDiffParticipantts: ProjX+="<<Pprojectile.plus()<<", ProjX-="
566        <<Pprojectile.minus()<<", TargX+="<< Ptarget.plus()<<", TargX-="<<Ptarget.minus()
567        <<G4endl;
568#endif
569         G4LorentzVector Qmomentum(0.,0.,0.,0.);
570         G4int whilecount=0;
571         do
572  {
573    //  Generate pt             
574           G4double maxPtSquare=sqr(Ptarget.pz());
575           if(whilecount++>=500 && whilecount%100==0) // @@ M.K. Hardwired limits
576#ifdef debug
577                 G4cout<<"G4QFragmentation::ExciteDiffParticipantts: can loop, loopCount="<<whilecount
578          <<", maxPtSquare="<<maxPtSquare<<G4endl;
579#endif
580    if(whilecount>1000)                        // @@ M.K. Hardwired limits
581    {
582#ifdef debug
583      G4cout<<"G4QFragmentation::ExciteDiffParticipants: *2* abort Loop!! *2*"<<G4endl;
584#endif
585                   return false;          //  Ignore this interaction
586    }
587           Qmomentum=G4LorentzVector(GaussianPt(widthOfPtSquare,maxPtSquare),0);
588#ifdef debug
589    G4cout<<"G4QFragment::ExciteDiffParticipants: generated Pt="<<Qmomentum<<", ProjPt="
590          <<Pprojectile+Qmomentum<<", TargPt="<<Ptarget-Qmomentum<<G4endl;
591#endif
592    //  Momentum transfer
593           G4double Xmin = minmass/(Pprojectile.e() + Ptarget.e());
594           G4double Xmax=1.;
595           G4double Xplus =ChooseX(Xmin,Xmax);
596           G4double Xminus=ChooseX(Xmin,Xmax);
597#ifdef debug
598           G4cout<<"G4QFragment::ExciteDiffParticip: X-plus="<<Xplus<<",X-minus="<<Xminus<<G4endl;
599#endif
600           G4double pt2=G4ThreeVector(Qmomentum.vect()).mag2();
601           G4double Qplus =-pt2/Xminus/Ptarget.minus();
602           G4double Qminus= pt2/Xplus /Pprojectile.plus();
603           Qmomentum.setPz((Qplus-Qminus)/2);
604           Qmomentum.setE( (Qplus+Qminus)/2);
605#ifdef debug
606           G4cout<<"G4QFragment::ExciteDiffParticip: Qplus="<<Qplus<<", Qminus="<<Qminus<<", pt2="
607          <<pt2<<", Qmomentum="<<Qmomentum<<", ProjM="<<(Pprojectile+Qmomentum).mag()
608          <<", TargM="<<(Ptarget-Qmomentum).mag()<<G4endl;
609#endif
610         } while((Pprojectile+Qmomentum).mag2()<=Mprojectile2 ||
611          (Ptarget-Qmomentum).mag2()<=Mtarget2);
612         Pprojectile += Qmomentum;
613         Ptarget     -= Qmomentum;
614#ifdef debug
615         G4cout<<"G4QFragment::ExciteDiffParticipan: Proj(Q)="<<Pprojectile<<", Targ(Q)="<<Ptarget
616               <<", Proj(back)="<<toLab*Pprojectile<<", Targ(bac)="<< toLab*Ptarget << G4endl;
617#endif
618  // Transform back and update SplitableHadron Participant.
619         Pprojectile.transform(toLab);
620         Ptarget.transform(toLab);
621#ifdef debug
622         G4cout<< "G4QFragmentation::ExciteDiffParticipants: TargetMass="<<Ptarget.mag()<<G4endl;
623#endif
624         target->Set4Momentum(Ptarget);         
625#ifdef debug
626         G4cout<<"G4QFragment::ExciteDiffParticipants:ProjectileMass="<<Pprojectile.mag()<<G4endl;
627#endif
628         projectile->Set4Momentum(Pprojectile);
629         return true;
630} // End of ExciteDiffParticipants
631
632
633// Excite single diffractive string
634G4bool G4QFragmentation::ExciteSingDiffParticipants(G4QHadron* projectile,
635                                                    G4QHadron* target) const
636{
637         G4LorentzVector Pprojectile=projectile->Get4Momentum();
638         G4double Mprojectile=projectile->GetMass() + minExtraMass;
639         G4double Mprojectile2=Mprojectile*Mprojectile;
640  G4LorentzVector Ptarget=target->Get4Momentum();
641  G4double Mtarget=target->GetMass() + minExtraMass;
642  G4double Mtarget2=Mtarget*Mtarget;
643#ifdef debug
644         G4cout<<"G4QFragm::ExSingDiffPartici:Ep="<<Pprojectile.e()<<",Et="<<Ptarget.e()<<G4endl;
645#endif
646         G4bool KeepProjectile= G4UniformRand() > 0.5;
647  // Reset minMass of the non diffractive particle to its value, (minus for rounding...)
648         if(KeepProjectile ) 
649         {
650#ifdef debug
651    G4cout<<"--1/2--G4QFragmentation::ExSingDiffParticipants: Projectile is fixed"<<G4endl;
652#endif
653           Mprojectile2 = projectile->GetMass2()*(1.-perCent); // Isn't it too big reduction? M.K.
654         }
655  else
656  {
657#ifdef debug
658    G4cout<<"---1/2---G4QFragmentation::ExSingDiffParticipants: Target is fixed"<<G4endl;
659#endif
660           Mtarget2 = target->GetMass2()*(1.-perCent); // @@ Isn't it too big reduction? M.K.
661         }
662  // @@ From this point it repeats the Diffractional excitation (? Use flag ?)
663  // Transform momenta to cms and then rotate parallel to z axis;
664         G4LorentzVector Psum=Pprojectile+Ptarget;
665         G4LorentzRotation toCms(-Psum.boostVector()); // Boost Rotation to CMS
666         G4LorentzVector Ptmp=toCms*Pprojectile;
667         if(Ptmp.pz()<=0.) // "String" moving backwards in CMS, abort collision !! ? M.K.
668         {
669#ifdef debug
670           G4cout<<"G4QFragment::ExciteSingDiffParticipants: *1* abort Collision!! *1*"<<G4endl;
671#endif
672                  return false; 
673         }                         
674         toCms.rotateZ(-Ptmp.phi());
675         toCms.rotateY(-Ptmp.theta());
676#ifdef debug
677  G4cout<<"G4QFragm::ExciteSingDiffParticipantts: Be4Boost Pproj="<<Pprojectile<<",Ptarg="
678        <<Ptarget<<G4endl;
679#endif
680         G4LorentzRotation toLab(toCms.inverse()); // Boost Rotation to LabSys (LS)
681         Pprojectile.transform(toCms);
682         Ptarget.transform(toCms);
683#ifdef debug
684         G4cout<< "G4QFragment::ExciteDiffParticipantts: AfterBoost Pproj="<<Pprojectile<<"Ptarg="
685        <<Ptarget<<", cms4M="<<Pprojectile+Ptarget<<G4endl;
686
687         G4cout<<"G4QFragment::ExciteDiffParticipantts: ProjX+="<<Pprojectile.plus()<<", ProjX-="
688        <<Pprojectile.minus()<<", TargX+="<< Ptarget.plus()<<", TargX-="<<Ptarget.minus()
689        <<G4endl;
690#endif
691         G4LorentzVector Qmomentum(0.,0.,0.,0.);
692         G4int whilecount=0;
693         do
694  {
695    //  Generate pt             
696           G4double maxPtSquare=sqr(Ptarget.pz());
697           if(whilecount++>=500 && whilecount%100==0) // @@ M.K. Hardwired limits
698#ifdef debug
699                 G4cout<<"G4QFragment::ExciteSingDiffParticipantts: can loop, loopCount="<<whilecount
700          <<", maxPtSquare="<<maxPtSquare<<G4endl;
701#endif
702    if(whilecount>1000)                        // @@ M.K. Hardwired limits
703    {
704#ifdef debug
705      G4cout<<"G4QFragmentation::ExciteSingDiffParticipants: *2* abort Loop!! *2*"<<G4endl;
706#endif
707                   return false;          //  Ignore this interaction
708    }
709           Qmomentum=G4LorentzVector(GaussianPt(widthOfPtSquare,maxPtSquare),0);
710#ifdef debug
711    G4cout<<"G4QFragm::ExciteSingDiffParticipants: generated Pt="<<Qmomentum<<", ProjPt="
712          <<Pprojectile+Qmomentum<<", TargPt="<<Ptarget-Qmomentum<<G4endl;
713#endif
714    //  Momentum transfer
715           G4double Xmin = minmass/(Pprojectile.e() + Ptarget.e());
716           G4double Xmax=1.;
717           G4double Xplus =ChooseX(Xmin,Xmax);
718           G4double Xminus=ChooseX(Xmin,Xmax);
719#ifdef debug
720           G4cout<<"G4QFragm::ExciteSingDiffPartici:X-plus="<<Xplus<<",X-minus="<<Xminus<<G4endl;
721#endif
722           G4double pt2=G4ThreeVector(Qmomentum.vect()).mag2();
723           G4double Qplus =-pt2/Xminus/Ptarget.minus();
724           G4double Qminus= pt2/Xplus /Pprojectile.plus();
725           if (KeepProjectile)
726             Qminus=(projectile->GetMass2()+pt2)/(Pprojectile.plus()+Qplus) - Pprojectile.minus();
727           else Qplus=Ptarget.plus() - (target->GetMass2()+pt2)/(Ptarget.minus()-Qminus);               
728           Qmomentum.setPz((Qplus-Qminus)/2);
729           Qmomentum.setE( (Qplus+Qminus)/2);
730#ifdef debug
731           G4cout<<"G4QFragm::ExciteDiffParticip: Qplus="<<Qplus<<", Qminus="<<Qminus<<", pt2="
732          <<pt2<<", Qmomentum="<<Qmomentum<<", ProjM="<<(Pprojectile+Qmomentum).mag()
733          <<", TargM="<<(Ptarget-Qmomentum).mag()<<G4endl;
734#endif
735    // while is different from the Double Diffractive Excitation (@@ !)
736                  //} while((Pprojectile+Qmomentum).mag2()<= Mprojectile2 ||
737    //        (Ptarget-Qmomentum).mag2()<=Mtarget2);
738         } while((Ptarget-Qmomentum).mag2()<=Mtarget2 ||
739          (Pprojectile+Qmomentum).mag2()<=Mprojectile2 ||
740          (Ptarget-Qmomentum).e() < 0. || (Pprojectile+Qmomentum).e() < 0.);
741         Pprojectile += Qmomentum;
742         Ptarget     -= Qmomentum;
743#ifdef debug
744         G4cout<<"G4QFragmentation::ExciteSingDiffParticipan: Proj(Q)="<<Pprojectile<<"(E="
745        <<Pprojectile.e()<<"), Targ(Q)="<<Ptarget<<"(E="<<Ptarget.e()
746               <<"), Proj(back)="<<toLab*Pprojectile<<", Targ(bac)="<< toLab*Ptarget << G4endl;
747#endif
748  // Transform back and update SplitableHadron Participant.
749         Pprojectile.transform(toLab);
750         Ptarget.transform(toLab);
751#ifdef debug
752         G4cout<< "G4QFragm::ExciteSingDiffParticipants: TargetMass="<<Ptarget.mag()<<G4endl;
753#endif
754         target->Set4Momentum(Ptarget);         
755#ifdef debug
756         G4cout<<"G4QFragm::ExciteDiffParticipants:ProjectileMass="<<Pprojectile.mag()<<G4endl;
757#endif
758         projectile->Set4Momentum(Pprojectile);
759         return true;
760} // End of ExciteSingleDiffParticipants
761
762void G4QFragmentation::SetParameters(G4int nCM, G4double thresh, G4double QGSMth,
763                           G4double radNuc, G4double SigPt, G4double extraM, G4double minM)
764{//  =============================================================================
765  nCutMax            = nCM;            // max number of pomeron cuts
766  ThersholdParameter = thresh;         // internal threshold
767  QGSMThershold      = QGSMth;         // QGSM threshold
768  theNucleonRadius   = radNuc;         // effective radius of the nucleon inside Nucleus
769         widthOfPtSquare    = -2*SigPt*SigPt;   // width^2 of pt for string excitation
770         minExtraMass       = extraM;           // minimum excitation mass
771         minmass            = minM;               // mean pion transverse mass; used for Xmin
772}
773
774G4double G4QFragmentation::ChooseX(G4double Xmin, G4double Xmax) const
775{
776// choose an x between Xmin and Xmax with P(x) ~ 1/x
777//  to be improved...
778         G4double range=Xmax-Xmin;
779         if( Xmin<= 0. || range <=0.) 
780         {
781                  G4cerr<<"***G4QFragmentation::ChooseX: Xmin="<<Xmin<<", Xmax="<<Xmax<< G4endl;
782    G4Exception("G4QFragmentation::ChooseX:","72",FatalException,"BadXRange");
783         }
784         G4double x;
785         do {x=Xmin+G4UniformRand()*range;} while ( Xmin/x < G4UniformRand() );
786#ifdef debug
787  G4cout<<"G4QFragmentation::ChooseX: DiffractiveX="<<x<<G4endl;
788#endif
789         return x;
790} // End of ChooseX
791
792// Pt distribution @@ one can use 1/(1+A*Pt^2)^B
793G4ThreeVector G4QFragmentation::GaussianPt(G4double widthSq, G4double maxPtSquare) const
794{
795         G4double pt2; do{pt2=widthSq*std::log(G4UniformRand());} while (pt2>maxPtSquare);
796                pt2=std::sqrt(pt2);
797         G4double phi=G4UniformRand()*twopi;
798         return G4ThreeVector(pt2*std::cos(phi),pt2*std::sin(phi),0.);   
799} // End of GaussianPt
Note: See TracBrowser for help on using the repository browser.