source: trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4DiffractiveExcitation.cc @ 962

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

update processes
File size: 29.8 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: G4DiffractiveExcitation.cc,v 1.7 2008/12/18 13:01:58 gunter Exp $
28// ------------------------------------------------------------
29//      GEANT 4 class implemetation file
30//
31//      ---------------- G4DiffractiveExcitation --------------
32//             by Gunter Folger, October 1998.
33//        diffractive Excitation used by strings models
34//             Take a projectile and a target
35//             excite the projectile and target
36//  Essential changed by V. Uzhinsky in November - December 2006
37//  in order to put it in a correspondence with original FRITIOF
38//  model. Variant of FRITIOF with nucleon de-excitation is implemented.
39//  Other changes by V.Uzhinsky in May 2007 were introduced to fit
40//  meson-nucleon interactions. Additional changes by V. Uzhinsky
41//  were introduced in December 2006. They treat diffraction dissociation
42//  processes more exactly.
43// ---------------------------------------------------------------------
44
45
46#include "globals.hh"
47#include "Randomize.hh"
48
49#include "G4DiffractiveExcitation.hh"
50#include "G4LorentzRotation.hh"
51#include "G4ThreeVector.hh"
52#include "G4ParticleDefinition.hh"
53#include "G4VSplitableHadron.hh"
54#include "G4ExcitedString.hh"
55#include "G4FTFParameters.hh"                            // Uzhi 19.04.08
56//#include "G4ios.hh"
57//#include "UZHI_diffraction.hh"
58
59G4DiffractiveExcitation::G4DiffractiveExcitation()
60{
61}
62
63// ---------------------------------------------------------------------
64G4bool G4DiffractiveExcitation::
65  ExciteParticipants(G4VSplitableHadron *projectile, 
66                     G4VSplitableHadron *target,
67                     G4FTFParameters    *theParameters) const
68{
69     G4bool PutOnMassShell=0;
70
71// -------------------- Projectile parameters -----------------------
72
73     G4LorentzVector Pprojectile=projectile->Get4Momentum();
74//   G4double M0projectile=projectile->GetDefinition()->GetPDGMass(); // With de-excitation
75     G4double M0projectile = Pprojectile.mag();                       // Without de-excitation
76/*
77G4cout<<"ExciteParticipants-------------------"<<G4endl;
78G4cout<<"Mom "<<Pprojectile<<" mass "<<M0projectile<<G4endl;
79*/
80     if(M0projectile < projectile->GetDefinition()->GetPDGMass())
81     {
82      PutOnMassShell=1;
83      M0projectile=projectile->GetDefinition()->GetPDGMass();
84     }
85
86     G4double M0projectile2 = M0projectile * M0projectile;
87
88     G4int    PDGcode=projectile->GetDefinition()->GetPDGEncoding();
89     G4int    absPDGcode=std::abs(PDGcode);
90
91     G4double ProjectileDiffStateMinMass=theParameters->GetProjMinDiffMass();
92     G4double ProjectileNonDiffStateMinMass=theParameters->GetProjMinNonDiffMass();
93     G4double ProbProjectileDiffraction=theParameters->GetProbabilityOfProjDiff();
94/*
95G4cout<<ProjectileDiffStateMinMass<<" "<<ProjectileNonDiffStateMinMass<<" "<<ProbProjectileDiffraction<<G4endl;
96*/
97// -------------------- Target paraExciteParticipantsmeters -------------------------
98     G4LorentzVector Ptarget=target->Get4Momentum();
99     G4double M0target = Ptarget.mag();
100
101//G4cout<<"Mom "<<Ptarget<<" mass "<<M0target<<G4endl;
102
103     if(M0target < target->GetDefinition()->GetPDGMass())
104     {
105      PutOnMassShell=1;
106      M0target=target->GetDefinition()->GetPDGMass();
107     }
108
109     G4double M0target2 = M0target * M0target;             //Ptarget.mag2();
110                                                           // for AA-inter.
111     G4double TargetDiffStateMinMass=theParameters->GetTarMinDiffMass();   
112     G4double TargetNonDiffStateMinMass=theParameters->GetTarMinNonDiffMass();   
113     G4double ProbTargetDiffraction=theParameters->GetProbabilityOfTarDiff();
114/*
115G4cout<<TargetDiffStateMinMass<<" "<<TargetNonDiffStateMinMass<<" "<<ProbTargetDiffraction<<G4endl;
116*/
117     G4double AveragePt2=theParameters->GetAveragePt2();
118
119// Kinematical properties of the interactions --------------
120     G4LorentzVector Psum;      // 4-momentum in CMS
121     Psum=Pprojectile+Ptarget;
122     G4double S=Psum.mag2(); 
123
124//G4cout<<" sqrt(s) "<<std::sqrt(S)<<G4endl;
125
126// ------------------------------------------------------------------
127
128//ProbProjectileDiffraction=1.;
129//ProbTargetDiffraction    =1.;
130     G4double ProbOfDiffraction=ProbProjectileDiffraction +
131                                ProbTargetDiffraction;
132
133     if(ProbOfDiffraction!=0.)
134     {
135      ProbProjectileDiffraction/=ProbOfDiffraction;
136     }
137     else
138     {
139      ProbProjectileDiffraction=0.;
140     }
141//   ProbTargetDiffraction    /=ProbOfDiffraction;
142
143//G4cout<<"ProbOfDiffraction "<<ProbOfDiffraction<<"ProbProjectileDiffraction "<<ProbProjectileDiffraction<<G4endl;   // Vova
144
145     G4double ProjectileDiffStateMinMass2    = ProjectileDiffStateMinMass    *
146                                               ProjectileDiffStateMinMass;
147     G4double ProjectileNonDiffStateMinMass2 = ProjectileNonDiffStateMinMass *
148                                               ProjectileNonDiffStateMinMass;
149
150     G4double TargetDiffStateMinMass2        = TargetDiffStateMinMass        *
151                                               TargetDiffStateMinMass;
152     G4double TargetNonDiffStateMinMass2     = TargetNonDiffStateMinMass     *
153                                               TargetNonDiffStateMinMass;
154
155// Transform momenta to cms and then rotate parallel to z axis;
156
157//         G4LorentzVector Psum;
158//         Psum=Pprojectile+Ptarget;
159
160     G4LorentzRotation toCms(-1*Psum.boostVector());
161
162     G4LorentzVector Ptmp=toCms*Pprojectile;
163     if ( Ptmp.pz() <= 0. )
164        {
165           // "String" moving backwards in  CMS, abort collision !!
166           //G4cout << " abort Collision!! " << G4endl;
167          return false;
168         }
169
170     toCms.rotateZ(-1*Ptmp.phi());
171     toCms.rotateY(-1*Ptmp.theta());
172
173     G4LorentzRotation toLab(toCms.inverse());
174
175     Pprojectile.transform(toCms);
176     Ptarget.transform(toCms);
177
178     G4double Pt2;
179     G4double ProjMassT2, ProjMassT;
180     G4double TargMassT2, TargMassT;
181     G4double PZcms2, PZcms;
182     G4double PMinusMin, PMinusMax;
183//     G4double PPlusMin , PPlusMax;
184     G4double TPlusMin , TPlusMax;
185     G4double PMinusNew, PPlusNew, TPlusNew, TMinusNew;
186
187//   G4double S=Psum.mag2();
188     G4double SqrtS=std::sqrt(S);
189
190     if(absPDGcode > 1000 && SqrtS < 2200*MeV)
191     {return false;}                         // The model cannot work for
192                                             // p+p-interactions
193                                             // at Plab < 1.3 GeV/c.
194
195     if(( absPDGcode == 211 || PDGcode ==  111) && SqrtS < 1600*MeV)
196     {return false;}                         // The model cannot work for
197                                             // Pi+p-interactions
198                                             // at Plab < 1. GeV/c.
199
200     if(( absPDGcode == 321 || PDGcode == -311) && SqrtS < 1600*MeV)
201     {return false;}                         // The model cannot work for
202                                             // K+p-interactions
203                                             // at Plab < ??? GeV/c.  ???
204
205     PZcms2=(S*S+M0projectile2*M0projectile2+M0target2*M0target2-
206             2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2)
207             /4./S;
208
209     if(PZcms2 < 0)
210     {return false;}   // It can be in an interaction with off-shell nuclear nucleon
211
212     PZcms = std::sqrt(PZcms2);
213
214     if(PutOnMassShell)
215     {
216      if(Pprojectile.z() > 0.)
217      {
218       Pprojectile.setPz( PZcms);
219       Ptarget.setPz(    -PZcms);
220      }
221      else
222      {
223       Pprojectile.setPz(-PZcms);
224       Ptarget.setPz(     PZcms);
225      };
226
227      Pprojectile.setE(std::sqrt(M0projectile2                  +
228                                 Pprojectile.x()*Pprojectile.x()+
229                                 Pprojectile.y()*Pprojectile.y()+
230                                 PZcms2));
231      Ptarget.setE(std::sqrt(M0target2              +
232                             Ptarget.x()*Ptarget.x()+
233                             Ptarget.y()*Ptarget.y()+
234                             PZcms2));
235     }
236
237     G4double maxPtSquare; // = PZcms2;
238/*
239G4cout << "Pprojectile aft boost : " << Pprojectile <<" "<<Pprojectile.mag()<< G4endl;
240G4cout << "Ptarget aft boost : " << Ptarget <<" "<<Ptarget.mag()<< G4endl;
241G4cout << "cms aft boost : " << (Pprojectile+ Ptarget) << G4endl;
242G4cout << " Projectile Xplus / Xminus : " <<
243            Pprojectile.plus() << " / " << Pprojectile.minus() << G4endl;
244G4cout << " Target Xplus / Xminus : " <<           Ptarget.plus() << " / " << Ptarget.minus() << G4endl;
245G4cout<<"maxPtSquare "<<maxPtSquare<<G4endl;
246*/
247     G4LorentzVector Qmomentum;
248     G4double Qminus, Qplus;
249
250     G4int whilecount=0;
251//  Choose a process
252
253     if(G4UniformRand() < ProbOfDiffraction)
254       {
255        if(G4UniformRand() < ProbProjectileDiffraction)
256        { //-------- projectile diffraction ---------------
257//G4cout<<"   Projectile diffraction"<<G4endl;
258//Uzhi_projectilediffraction++;
259         do {
260//             Generate pt
261//             if (whilecount++ >= 500 && (whilecount%100)==0)
262//               G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping"
263//               << ", loop count/ maxPtSquare : "
264//               << whilecount << " / " << maxPtSquare << G4endl;
265             if (whilecount > 1000 )
266             {
267              Qmomentum=G4LorentzVector(0.,0.,0.,0.);
268              return false;    //  Ignore this interaction
269             };
270// --------------- Check that the interaction is possible -----------
271             ProjMassT2=ProjectileDiffStateMinMass2;
272             ProjMassT =ProjectileDiffStateMinMass;
273
274             TargMassT2=M0target2;
275             TargMassT =M0target;
276
277             PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
278                     2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
279                    /4./S;
280//G4cout<<" Pt2 Mpt Mtt Pz2 "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl;
281
282if(PZcms2 < 0 ) 
283{
284/*
285G4cout<<"whilecount "<<whilecount<<" "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl;
286G4int Uzhi; G4cin>>Uzhi;
287*/
288return false;
289};
290             maxPtSquare=PZcms2;
291
292             Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
293             Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
294
295             ProjMassT2=ProjectileDiffStateMinMass2+Pt2;
296             ProjMassT =std::sqrt(ProjMassT2);
297
298             TargMassT2=M0target2+Pt2;
299             TargMassT =std::sqrt(TargMassT2);
300
301             PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
302                     2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
303                    /4./S;
304//G4cout<<" Pt2 Mpt Mtt Pz2 "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl;
305
306//             if(PZcms2 < 0 ) {PZcms2=0;};
307             if(PZcms2 < 0 ) continue;
308             PZcms =std::sqrt(PZcms2);
309
310             PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
311             PMinusMax=SqrtS-TargMassT;
312//G4cout<<" SqrtS P+mim max "<<SqrtS<<" "<<PMinusMin<<" "<<PMinusMax<<G4endl;
313
314             PMinusNew=ChooseP(PMinusMin, PMinusMax);
315// PMinusNew=1./sqrt(1./PMinusMin-G4UniformRand()*(1./PMinusMin-1./PMinusMax));
316
317             TMinusNew=SqrtS-PMinusNew;
318             Qminus=Ptarget.minus()-TMinusNew;
319             TPlusNew=TargMassT2/TMinusNew;
320             Qplus=Ptarget.plus()-TPlusNew;
321
322             Qmomentum.setPz( (Qplus-Qminus)/2 );
323             Qmomentum.setE(  (Qplus+Qminus)/2 );
324          } while (
325((Pprojectile+Qmomentum).mag2() <  ProjectileDiffStateMinMass2) ||  //No without excitation
326((Ptarget    -Qmomentum).mag2() <  M0target2                  ));
327        }
328        else
329        { // -------------- Target diffraction ----------------
330//G4cout<<"   Target difraction"<<G4endl;
331//Uzhi_targetdiffraction++;
332         do {
333//             Generate pt
334//             if (whilecount++ >= 500 && (whilecount%100)==0)
335//               G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping"
336//               << ", loop count/ maxPtSquare : "
337//               << whilecount << " / " << maxPtSquare << G4endl;
338             if (whilecount > 1000 )
339             {
340              Qmomentum=G4LorentzVector(0.,0.,0.,0.);
341              return false;    //  Ignore this interaction
342             };
343// --------------- Check that the interaction is possible -----------
344             ProjMassT2=M0projectile2;
345             ProjMassT =M0projectile;
346
347             TargMassT2=TargetDiffStateMinMass2;
348             TargMassT =TargetDiffStateMinMass;
349
350             PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
351                     2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
352                    /4./S;
353//G4cout<<" Pt2 Mpt Mtt Pz2 "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl;
354
355if(PZcms2 < 0 ) 
356{
357/*
358G4cout<<"whilecount "<<whilecount<<" "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl;
359G4int Uzhi; G4cin>>Uzhi;
360*/
361return false;
362};
363             maxPtSquare=PZcms2;
364
365             Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
366             Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
367
368             ProjMassT2=M0projectile2+Pt2;
369             ProjMassT =std::sqrt(ProjMassT2);
370
371             TargMassT2=TargetDiffStateMinMass2+Pt2;
372             TargMassT =std::sqrt(TargMassT2);
373
374             PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
375                     2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
376                    /4./S;
377/*
378if(PZcms2 < 0 )
379{
380G4cout<<"whilecount "<<whilecount<<" "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl;
381G4int Uzhi; G4cin>>Uzhi;
382return false;
383};
384*/
385             if(PZcms2 < 0 ) continue;
386             PZcms =std::sqrt(PZcms2);
387
388             TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
389             TPlusMax=SqrtS-ProjMassT;
390
391//G4cout<<" Tmin max "<<TPlusMin<<" "<<TPlusMax<<G4endl;
392
393             TPlusNew=ChooseP(TPlusMin, TPlusMax);
394
395//TPlusNew=TPlusMax;
396//G4cout<<"T+new "<<TPlusNew<<G4endl;
397
398             PPlusNew=SqrtS-TPlusNew;
399             Qplus=PPlusNew-Pprojectile.plus();
400             PMinusNew=ProjMassT2/PPlusNew;
401             Qminus=PMinusNew-Pprojectile.minus();
402
403             Qmomentum.setPz( (Qplus-Qminus)/2 );
404             Qmomentum.setE(  (Qplus+Qminus)/2 );
405
406          } while (
407 ((Pprojectile+Qmomentum).mag2() <  M0projectile2          ) ||  //No without excitation
408 ((Ptarget    -Qmomentum).mag2() <  TargetDiffStateMinMass2));
409         }
410        }
411        else  //----------- Non-diffraction process ------------
412        {
413//G4cout<<"   Non-difraction"<<G4endl;
414         do {
415//             Generate pt
416//             if (whilecount++ >= 500 && (whilecount%100)==0)
417//               G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping"
418//               << ", loop count/ maxPtSquare : "
419//               << whilecount << " / " << maxPtSquare << G4endl;
420             if (whilecount > 1000 )
421             {
422              Qmomentum=G4LorentzVector(0.,0.,0.,0.);
423              return false;    //  Ignore this interaction
424             };
425// --------------- Check that the interaction is possible -----------
426             ProjMassT2=ProjectileNonDiffStateMinMass2;
427             ProjMassT =ProjectileNonDiffStateMinMass;
428
429             TargMassT2=TargetNonDiffStateMinMass2;
430             TargMassT =TargetNonDiffStateMinMass;
431
432             PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
433                    2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
434                   /4./S;
435//G4cout<<" Pt2 Mpt Mtt Pz2 "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl;
436
437if(PZcms2 < 0 ) 
438{
439/*
440G4cout<<"whilecount "<<whilecount<<" "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<PZcms2<<G4endl;
441G4int Uzhi; G4cin>>Uzhi;
442*/
443return false;
444};
445             maxPtSquare=PZcms2;
446             Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
447             Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
448
449             ProjMassT2=ProjectileNonDiffStateMinMass2+Pt2;
450             ProjMassT =std::sqrt(ProjMassT2);
451
452             TargMassT2=TargetNonDiffStateMinMass2+Pt2;
453             TargMassT =std::sqrt(TargMassT2);
454
455             PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
456                    2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
457                   /4./S;
458/*
459G4cout<<"ProjectileNonDiffStateMinMass2 "<<ProjectileNonDiffStateMinMass2<<G4endl;
460G4cout<<"TargetNonDiffStateMinMass2     "<<TargetNonDiffStateMinMass2<<G4endl;
461G4cout<<"Mt "<<ProjMassT<<" "<<TargMassT<<" "<<Pt2<<" "<<PZcms2<<G4endl<<G4endl;
462*/
463//             if(PZcms2 < 0 ) {PZcms2=0;};
464             if(PZcms2 < 0 ) continue;
465             PZcms =std::sqrt(PZcms2);
466
467             PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
468             PMinusMax=SqrtS-TargMassT;
469
470             PMinusNew=ChooseP(PMinusMin, PMinusMax);
471// PMinusNew=1./sqrt(1./PMinusMin-G4UniformRand()*(1./PMinusMin-1./PMinusMax));
472
473//G4cout<<"Proj "<<PMinusMin<<" "<<PMinusMax<<" "<<PMinusNew<<G4endl;
474
475//PMinusNew=PMinusMax; //+++++++++++++++++++++++++++++++++++ Vova
476
477             Qminus=PMinusNew-Pprojectile.minus();
478
479             TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
480//             TPlusMax=SqrtS-PMinusNew;                      // Vova
481             TPlusMax=SqrtS-ProjMassT;                        // Vova
482
483             TPlusNew=ChooseP(TPlusMin, TPlusMax);
484
485//G4cout<<"Targ "<<TPlusMin<<" "<<TPlusMax<<" "<<TPlusNew<<G4endl;
486//G4cout<<PMinusNew<<" "<<TPlusNew<<G4endl;
487
488             Qplus=-(TPlusNew-Ptarget.plus());
489
490             Qmomentum.setPz( (Qplus-Qminus)/2 );
491             Qmomentum.setE(  (Qplus+Qminus)/2 );
492/*
493G4cout << "Qplus / Qminus " << Qplus << " / " << Qminus<<G4endl;
494G4cout << "pt2" << pt2 << G4endl;
495G4cout << "Qmomentum " << Qmomentum << G4endl;
496G4cout << " Masses (P/T) : " << (Pprojectile+Qmomentum).mag() <<
497           " / " << (Ptarget-Qmomentum).mag() << G4endl;   // mag()
498G4cout<<"Mprojectile "<<std::sqrt(M0projectile2)<<G4endl;
499G4cout<<"Mtarget     "<<std::sqrt(M0target2    )<<G4endl;
500G4cout<<"ProjectileDiffStateMinMass "<<std::sqrt(ProjectileDiffStateMinMass2)<<G4endl;
501G4cout<<"TargetDiffStateMinMass "<<std::sqrt(TargetDiffStateMinMass2)<<G4endl;
502*/
503       } while (
504 ((Pprojectile+Qmomentum).mag2() <  ProjectileNonDiffStateMinMass2) || //No double Diffraction
505 ((Ptarget    -Qmomentum).mag2() <  TargetNonDiffStateMinMass2    ));
506         }
507
508//G4int Uzhiinp; G4cin>>Uzhiinp;    // Vova
509
510           Pprojectile += Qmomentum;
511           Ptarget     -= Qmomentum;
512/*
513G4cout << "Pprojectile with Q : " << Pprojectile << G4endl;
514G4cout << "Ptarget with Q : " << Ptarget << G4endl;
515G4cout << "Target        mass  " <<  Ptarget.mag() << G4endl;
516G4cout << "Projectile mass  " <<  Pprojectile.mag() << G4endl;
517//
518//G4cout << "Projectile back: " << toLab * Pprojectile << G4endl;
519//G4cout << "Target back: " << toLab * Ptarget << G4endl;
520*/
521//-------------- Flip if projectale moves in backward direction ------------
522//G4bool Flip=Pprojectile.pz()< 0.;
523
524
525// Transform back and update SplitableHadron Participant.
526           Pprojectile.transform(toLab);
527           Ptarget.transform(toLab);
528
529//G4cout << "Pprojectile with Q M: " << Pprojectile<<" "<<  Pprojectile.mag() << G4endl;
530//G4cout << "Ptarget     with Q M: " << Ptarget    <<" "<<  Ptarget.mag()     << G4endl;
531//G4cout << "Target      mass  " <<  Ptarget.mag() << G4endl;
532//G4cout << "Projectile mass  " <<  Pprojectile.mag() << G4endl;
533
534/*
535if(!Flip){
536           projectile->Set4Momentum(Pprojectile);
537           target->Set4Momentum(Ptarget);
538         }
539else     {
540           G4ParticleDefinition * t_Definition=projectile->GetDefinition();
541           projectile->SetDefinition(target->GetDefinition());
542           projectile->Set4Momentum(Ptarget);
543           target->SetDefinition(t_Definition);
544           target->Set4Momentum(Pprojectile);
545          }
546*/
547//
548/*
549if(G4UniformRand() < 1.) {
550           G4ParticleDefinition * t_Definition=projectile->GetDefinition();
551           projectile->SetDefinition(target->GetDefinition());
552           target->SetDefinition(t_Definition);
553}
554*/ // For flip, for HARP
555
556           G4double ZcoordinateOfCurrentInteraction = target->GetPosition().z();
557// It is assumed that nucleon z-coordinates are ordered on increasing -----------
558
559           G4double betta_z=projectile->Get4Momentum().pz()/projectile->Get4Momentum().e();
560
561           G4double ZcoordinateOfPreviousCollision=projectile->GetPosition().z();
562           if(projectile->GetSoftCollisionCount()==0) {
563              projectile->SetTimeOfCreation(0.);
564              target->SetTimeOfCreation(0.);
565              ZcoordinateOfPreviousCollision=ZcoordinateOfCurrentInteraction;
566           }
567           
568           G4ThreeVector thePosition(projectile->GetPosition().x(),
569                                     projectile->GetPosition().y(),
570                                     ZcoordinateOfCurrentInteraction);
571           projectile->SetPosition(thePosition);
572
573           G4double TimeOfPreviousCollision=projectile->GetTimeOfCreation();
574           G4double TimeOfCurrentCollision=TimeOfPreviousCollision+ 
575                    (ZcoordinateOfCurrentInteraction-ZcoordinateOfPreviousCollision)/betta_z; 
576
577           projectile->SetTimeOfCreation(TimeOfCurrentCollision);
578           target->SetTimeOfCreation(TimeOfCurrentCollision);
579
580           projectile->Set4Momentum(Pprojectile);
581           target->Set4Momentum(Ptarget);
582
583           projectile->IncrementCollisionCount(1);
584           target->IncrementCollisionCount(1);
585
586//
587//G4cout<<"Out of Excitation --------------------"<<G4endl;
588//G4int Uzhiinp; G4cin>>Uzhiinp;    // Vova
589
590           return true;
591}
592
593// ---------------------------------------------------------------------
594G4ExcitedString * G4DiffractiveExcitation::
595           String(G4VSplitableHadron * hadron, G4bool isProjectile) const
596{
597
598//G4cout<<"G4DiffractiveExcitation::String isProj"<<isProjectile<<G4endl;
599
600        hadron->SplitUp();
601        G4Parton *start= hadron->GetNextParton();
602        if ( start==NULL)
603        { G4cout << " G4FTFModel::String() Error:No start parton found"<< G4endl;
604          return NULL;
605        }
606        G4Parton *end  = hadron->GetNextParton();
607        if ( end==NULL)
608        { G4cout << " G4FTFModel::String() Error:No end parton found"<< G4endl;
609          return NULL;
610        }
611
612        G4ExcitedString * string;
613        if ( isProjectile )
614        {
615                string= new G4ExcitedString(end,start, +1);
616        } else {
617                string= new G4ExcitedString(start,end, -1);
618        }
619// Uzhi
620//G4cout<<"G4ExcitedString * G4DiffractiveExcitation::String"<<G4endl;
621//G4cout<<hadron->GetTimeOfCreation()<<" "<<hadron->GetPosition()/fermi<<G4endl;
622
623        string->SetTimeOfCreation(hadron->GetTimeOfCreation());
624        string->SetPosition(hadron->GetPosition());
625
626// momenta of string ends
627//
628        G4double Momentum=hadron->Get4Momentum().vect().mag();
629        G4double  Plus=hadron->Get4Momentum().e() + Momentum;
630        G4double Minus=hadron->Get4Momentum().e() - Momentum;
631
632        G4ThreeVector tmp;
633        if(Momentum > 0.) 
634        {
635         tmp.set(hadron->Get4Momentum().px(),
636                 hadron->Get4Momentum().py(),
637                 hadron->Get4Momentum().pz());
638         tmp/=Momentum;
639         }
640         else
641         {
642          tmp.set(0.,0.,1.);
643         };
644
645        G4LorentzVector Pstart(tmp,0.);
646        G4LorentzVector   Pend(tmp,0.);
647
648        if(isProjectile)
649        {
650         Pstart*=(-1.)*Minus/2.;
651         Pend  *=(+1.)*Plus /2.;
652         } 
653        else
654         {
655          Pstart*=(+1.)*Plus/2.;
656          Pend  *=(-1.)*Minus/2.;
657         };
658
659         Momentum=-Pstart.mag();
660         Pstart.setT(Momentum);  // It is assumed that quark has m=0.
661
662         Momentum=-Pend.mag();
663         Pend.setT(Momentum);    // It is assumed that di-quark has m=0.
664//
665/* Uzhi
666        G4double ptSquared= hadron->Get4Momentum().perp2();
667        G4double transverseMassSquared= hadron->Get4Momentum().plus()
668                                    *   hadron->Get4Momentum().minus();
669
670
671        G4double maxAvailMomentumSquared=
672                 sqr( std::sqrt(transverseMassSquared) - std::sqrt(ptSquared) );
673
674        G4double widthOfPtSquare = 0.25*GeV*GeV;       // Uzhi 11.07 <Pt^2>=0.25 ??????????????????
675        G4ThreeVector pt=GaussianPt(widthOfPtSquare,maxAvailMomentumSquared);
676
677        G4LorentzVector Pstart(G4LorentzVector(pt,0.));
678        G4LorentzVector Pend;
679        Pend.setPx(hadron->Get4Momentum().px() - pt.x());
680        Pend.setPy(hadron->Get4Momentum().py() - pt.y());
681
682        G4double tm1=hadron->Get4Momentum().minus() +
683          ( Pend.perp2()-Pstart.perp2() ) / hadron->Get4Momentum().plus();
684
685        G4double tm2= std::sqrt( std::max(0., sqr(tm1) -
686             4. * Pend.perp2() * hadron->Get4Momentum().minus()
687              /  hadron->Get4Momentum().plus() ));
688
689        G4int Sign= isProjectile ? -1 : 1;
690
691        G4double endMinus  = 0.5 * (tm1 + Sign*tm2);
692        G4double startMinus= hadron->Get4Momentum().minus() - endMinus;
693
694        G4double startPlus= Pstart.perp2() /  startMinus;
695        G4double endPlus  = hadron->Get4Momentum().plus() - startPlus;
696
697        Pstart.setPz(0.5*(startPlus - startMinus));
698        Pstart.setE(0.5*(startPlus + startMinus));
699
700        Pend.setPz(0.5*(endPlus - endMinus));
701        Pend.setE(0.5*(endPlus + endMinus));
702*/ // Uzhi
703        start->Set4Momentum(Pstart);
704        end->Set4Momentum(Pend);
705/*
706G4cout<<"G4DiffractiveExcitation::String hadro"<<hadron->Get4Momentum()<<" "<<hadron->Get4Momentum().mag2()<<G4endl;
707
708G4cout<<"G4DiffractiveExcitation::String start"<<start->Get4Momentum()<<" "<<start->GetPDGcode()<<G4endl;
709
710G4cout<<"G4DiffractiveExcitation::String end  "<<  end->Get4Momentum()<<" "<<  end->GetPDGcode()<<G4endl;
711G4int Uzhi; G4cin>>Uzhi;
712*/
713#ifdef G4_FTFDEBUG
714        G4cout << " generated string flavors          " 
715               << start->GetPDGcode() << " / " 
716               << end->GetPDGcode()                    << G4endl;
717        G4cout << " generated string momenta:   quark " 
718               << start->Get4Momentum() << "mass : " 
719               <<start->Get4Momentum().mag()           << G4endl;
720        G4cout << " generated string momenta: Diquark " 
721               << end ->Get4Momentum() 
722               << "mass : " <<end->Get4Momentum().mag()<< G4endl;
723        G4cout << " sum of ends                       " << Pstart+Pend << G4endl;
724        G4cout << " Original                          " << hadron->Get4Momentum() << G4endl;
725#endif
726
727        return string;
728}
729
730
731// --------- private methods ----------------------
732
733// ---------------------------------------------------------------------
734G4double G4DiffractiveExcitation::ChooseP(G4double Pmin, G4double Pmax) const // Uzhi
735{
736// choose an x between Xmin and Xmax with P(x) ~ 1/x
737//  to be improved...
738
739        G4double range=Pmax-Pmin;                                         // Uzhi
740
741        if ( Pmin <= 0. || range <=0. )
742        {
743                G4cout << " Pmin, range : " << Pmin << " , " << range << G4endl;
744                throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation::ChooseP : Invalid arguments ");
745        }
746
747        G4double P;
748/*                                                                          // Uzhi
749        do {
750            x=Xmin + G4UniformRand() * range;
751        }  while ( Xmin/x < G4UniformRand() );
752*/                                                                          // Uzhi
753
754        P=Pmin * std::pow(Pmax/Pmin,G4UniformRand());                       // Uzhi
755
756//debug-hpw     cout << "DiffractiveX "<<x<<G4endl;
757        return P;
758}
759
760// ---------------------------------------------------------------------
761G4ThreeVector G4DiffractiveExcitation::GaussianPt(G4double AveragePt2, 
762                                                  G4double maxPtSquare) const // Uzhi
763{            //  @@ this method is used in FTFModel as well. Should go somewhere common!
764
765        G4double Pt2;
766/*                                                                          // Uzhi
767        do {
768            pt2=widthSquare * std::log( G4UniformRand() );
769        } while ( pt2 > maxPtSquare);
770*/                                                                          // Uzhi
771
772        Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() * 
773                           (std::exp(-maxPtSquare/AveragePt2)-1.));// Uzhi
774
775        G4double Pt=std::sqrt(Pt2);
776
777        G4double phi=G4UniformRand() * twopi;
778
779        return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);
780}
781
782// ---------------------------------------------------------------------
783G4DiffractiveExcitation::G4DiffractiveExcitation(const G4DiffractiveExcitation &)
784{
785        throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation copy contructor not meant to be called");
786}
787
788
789G4DiffractiveExcitation::~G4DiffractiveExcitation()
790{
791}
792
793
794const G4DiffractiveExcitation & G4DiffractiveExcitation::operator=(const G4DiffractiveExcitation &)
795{
796        throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation = operator meant to be called");
797        return *this;
798}
799
800
801int G4DiffractiveExcitation::operator==(const G4DiffractiveExcitation &) const
802{
803        throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation == operator meant to be called");
804        return false;
805}
806
807int G4DiffractiveExcitation::operator!=(const G4DiffractiveExcitation &) const
808{
809        throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation != operator meant to be called");
810        return true;
811}
Note: See TracBrowser for help on using the repository browser.