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

Last change on this file since 1071 was 1055, checked in by garnier, 15 years ago

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