Ignore:
Timestamp:
Nov 5, 2010, 3:45:55 PM (14 years ago)
Author:
garnier
Message:

update ti head

Location:
trunk/source/processes/hadronic/models/parton_string/diffraction/src
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4DiffractiveExcitation.cc

    r1228 r1340  
    2525//
    2626//
    27 // $Id: G4DiffractiveExcitation.cc,v 1.21 2009/12/15 19:14:31 vuzhinsk Exp $
     27// $Id: G4DiffractiveExcitation.cc,v 1.22 2010/09/20 15:50:46 vuzhinsk Exp $
    2828// ------------------------------------------------------------
    2929//      GEANT 4 class implemetation file
    … …  
    108108     G4int    TargetPDGcode=target->GetDefinition()->GetPDGEncoding();
    109109     G4int    absTargetPDGcode=std::abs(TargetPDGcode);
     110//G4cout<<"Excit "<<ProjectilePDGcode<<" "<<TargetPDGcode<<G4endl;
    110111
    111112     G4LorentzVector Ptarget=target->Get4Momentum();
    … …  
    113114     G4double M0target = Ptarget.mag();
    114115
    115      G4double TargetRapidity = Ptarget.rapidity();
     116//   G4double TargetRapidity = Ptarget.rapidity();
    116117
    117118     if(M0target < target->GetDefinition()->GetPDGMass())
    … …  
    129130     G4double AveragePt2=theParameters->GetAveragePt2();
    130131
    131      G4double ProbOfDiffraction=ProbProjectileDiffraction +
    132                                 ProbTargetDiffraction;
     132//     G4double ProbOfDiffraction=ProbProjectileDiffraction +
     133//                                ProbTargetDiffraction;
    133134
    134135     G4double SumMasses=M0projectile+M0target+200.*MeV;
    … …  
    162163
    163164     G4double SqrtS=std::sqrt(S);
    164               
     165           
    165166     if(absProjectilePDGcode > 1000 && (SqrtS < 2300*MeV || SqrtS < SumMasses))
    166167     {target->SetStatus(2);  return false;}  // The model cannot work for
    … …  
    215216
    216217     G4double maxPtSquare; // = PZcms2;
    217 
     218/*
     219G4cout<<"Start --------------------"<<G4endl;
     220G4cout<<"Proj "<<M0projectile<<" "<<ProjectileDiffStateMinMass<<"  "<<ProjectileNonDiffStateMinMass<<G4endl;
     221G4cout<<"Targ "<<M0target    <<" "<<TargetDiffStateMinMass    <<" "<<TargetNonDiffStateMinMass<<G4endl;
     222G4cout<<"SqrtS "<<SqrtS<<G4endl;
     223G4cout<<"Rapid "<<ProjectileRapidity<<" "<<TargetRapidity<<G4endl;
     224*/
    218225// Charge exchange can be possible for baryons -----------------
    219226
    … …  
    223230     G4double DeltaProbAtQuarkExchange=theParameters->GetDeltaProbAtQuarkExchange();
    224231
     232//G4cout<<"Q exc "<<MagQuarkExchange<<" "<<SlopeQuarkExchange<<" "<<DeltaProbAtQuarkExchange<<G4endl;
    225233//     G4double NucleonMass=
    226234//              (G4ParticleTable::GetParticleTable()->FindParticle(2112))->GetPDGMass();     
    … …  
    228236              (G4ParticleTable::GetParticleTable()->FindParticle(2224))->GetPDGMass();
    229237
    230 // Check for possible quark excjane -----------------------------------
     238//G4cout<<MagQuarkExchange*std::exp(-SlopeQuarkExchange*(ProjectileRapidity - TargetRapidity))<<G4endl;
     239//G4cout<<MagQuarkExchange*std::exp(-SlopeQuarkExchange*(ProjectileRapidity))<<G4endl;
     240//G4cout<<MagQuarkExchange*std::exp(-SlopeQuarkExchange*(ProjectileRapidity - 1.36))<<G4endl;
     241//G4int Uzhi; G4cin>>Uzhi;
     242// Check for possible quark exchange -----------------------------------
     243
    231244     if(G4UniformRand() < MagQuarkExchange*
    232         std::exp(-SlopeQuarkExchange*(ProjectileRapidity - TargetRapidity)))
     245        std::exp(-SlopeQuarkExchange*ProjectileRapidity))  //TargetRapidity))) 1.45
    233246     {   
     247//        std::exp(-SlopeQuarkExchange*(ProjectileRapidity - 1.36)))  //TargetRapidity))) 1.45
     248//G4cout<<"Q exchange"<<G4endl;
    234249      G4int NewProjCode(0), NewTargCode(0);
    235250
    … …  
    249264      UnpackBaryon(TargetPDGcode, TargQ1, TargQ2, TargQ3);
    250265
     266//G4cout<<ProjQ1<<" "<<ProjQ2<<" "<<ProjQ3<<G4endl;
     267//G4cout<<TargQ1<<" "<<TargQ2<<" "<<TargQ3<<G4endl;
    251268// Sampling of exchanged quarks -------------------
    252269      G4int ProjExchangeQ(0);
    … …  
    326343        {ProjExchangeQ = ProjQ3;}
    327344
     345//G4cout<<"ProjExchangeQ "<<ProjExchangeQ<<G4endl;
    328346        if((ProjExchangeQ != TargQ1)||(G4UniformRand()<Same))
    329347        {
    … …  
    338356        }
    339357
     358//G4cout<<"ProjExchangeQ "<<ProjExchangeQ<<G4endl;
     359//G4cout<<"TargExchangeQ "<<TargExchangeQ<<G4endl;
    340360        if( Ksi < 0.333333 )
    341361        {ProjQ1=ProjExchangeQ;}
    … …  
    376396
    377397       NewProjCode = NewNucleonId(ProjQ1, ProjQ2, ProjQ3); // *****************************
     398
     399//G4cout<<"ProjQ1, ProjQ2, ProjQ3 "<<ProjQ1<<" "<<ProjQ2<<" "<<ProjQ3<<" "<<NewProjCode<<G4endl;
     400
     401G4int                 TestParticleID=NewProjCode;
     402G4ParticleDefinition* TestParticle=0;
     403G4double              TestParticleMass=DBL_MAX;
     404
     405TestParticle=G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode);
     406if(TestParticle) TestParticleMass=TestParticle->GetPDGMass();
    378407
    379408       if((ProjQ1==ProjQ2) && (ProjQ1==ProjQ3)) {NewProjCode +=2; ProjDeltaHasCreated=true;}
    … …  
    390419       }
    391420
     421G4ParticleDefinition* NewTestParticle=
     422                      G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode);
     423//G4cout<<"TestParticleMass NewTestParticle->GetPDGMass() "<<TestParticleMass<<" "<< NewTestParticle->GetPDGMass()<<G4endl;
     424//if(TestParticleMass < NewTestParticle->GetPDGMass()) {NewProjCode=TestParticleID;}
     425 
     426//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++=
     427
    392428       NewTargCode = NewNucleonId(TargQ1, TargQ2, TargQ3); // *****************************
     429
     430//G4cout<<"TargQ1, TargQ2, TargQ3 "<<TargQ1<<" "<<TargQ2<<" "<<TargQ3<<" "<<NewTargCode<<G4endl;
     431
     432TestParticleID=NewTargCode;
     433TestParticleMass=DBL_MAX;
     434
     435TestParticle=G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode);
     436if(TestParticle) TestParticleMass=TestParticle->GetPDGMass();
    393437
    394438       if((TargQ1==TargQ2) && (TargQ1==TargQ3)) {NewTargCode +=2; TargDeltaHasCreated=true;} 
    … …  
    405449       }         
    406450
     451NewTestParticle=G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode);
     452//G4cout<<"TestParticleMass NewTestParticle->GetPDGMass() "<<TestParticleMass<<" "<< NewTestParticle->GetPDGMass()<<G4endl;
     453//if(TestParticleMass < NewTestParticle->GetPDGMass()) {NewTargCode=TestParticleID;}
     454
     455//G4cout<<"NewProjCode NewTargCode "<<NewProjCode<<" "<<NewTargCode<<G4endl;
     456//G4int Uzhi; G4cin>>Uzhi;
     457
    407458       if((absProjectilePDGcode == NewProjCode) && (absTargetPDGcode == NewTargCode))
    408459       { // Nothing was changed! It is not right!?
    409460       }
    410461// Forming baryons --------------------------------------------------
    411 
     462if(ProjDeltaHasCreated) {ProbProjectileDiffraction=1.; ProbTargetDiffraction=0.;}
     463if(TargDeltaHasCreated) {ProbProjectileDiffraction=0.; ProbTargetDiffraction=1.;}
     464       if(ProjDeltaHasCreated)
     465       {
     466        M0projectile=
     467          (G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode))->GetPDGMass();
     468        M0projectile2 = M0projectile * M0projectile;
     469
     470        ProjectileDiffStateMinMass   =M0projectile+210.*MeV; //210 MeV=m_pi+70 MeV
     471        ProjectileNonDiffStateMinMass=M0projectile+210.*MeV; //210 MeV=m_pi+70 MeV
     472       }
     473
     474//      if(M0target <
     475//         (G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode))->GetPDGMass())
     476       if(TargDeltaHasCreated)
     477       {
     478        M0target=
     479          (G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode))->GetPDGMass();
     480        M0target2 = M0target * M0target;
     481
     482        TargetDiffStateMinMass   =M0target+210.*MeV;         //210 MeV=m_pi+70 MeV;   
     483        TargetNonDiffStateMinMass=M0target+210.*MeV;         //210 MeV=m_pi+70 MeV;
     484       }
    412485      } // End of if projectile is baryon ---------------------------
    413486
    … …  
    416489// in the ground states, we have to put ----------------------------------
    417490
    418       if(M0projectile <
    419          (G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode))->GetPDGMass())
     491/*
     492//      if(M0projectile <
     493//         (G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode))->GetPDGMass())
     494      if(ProjDeltaHasCreated)
    420495      {
    421496       M0projectile=
    422497         (G4ParticleTable::GetParticleTable()->FindParticle(NewProjCode))->GetPDGMass();
    423498       M0projectile2 = M0projectile * M0projectile;
     499
     500       ProjectileDiffStateMinMass   =M0projectile+160.*MeV; //160 MeV=m_pi+20 MeV
     501       ProjectileNonDiffStateMinMass=M0projectile+160.*MeV; //160 MeV=m_pi+20 MeV
    424502      }
    425503
    426       if(M0target <
    427          (G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode))->GetPDGMass())
     504//      if(M0target <
     505//         (G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode))->GetPDGMass())
     506      if(TargDeltaHasCreated)
    428507      {
    429508       M0target=
    430509         (G4ParticleTable::GetParticleTable()->FindParticle(NewTargCode))->GetPDGMass();
    431510       M0target2 = M0target * M0target;
     511
     512       TargetDiffStateMinMass   =M0target+160.*MeV;         //160 MeV=m_pi+20 MeV;   
     513       TargetNonDiffStateMinMass=M0target+160.*MeV;         //160 MeV=m_pi+20 MeV;
    432514      }
    433 
     515*/
    434516      PZcms2=(S*S+M0projectile2*M0projectile2+M0target2*M0target2-
    435517             2*S*M0projectile2 - 2*S*M0target2 - 2*M0projectile2*M0target2)
    436518             /4./S;
    437 
     519//G4cout<<"PZcms2 1 "<<PZcms2<<G4endl;
    438520      if(PZcms2 < 0) {return false;}  // It can be if energy is not sufficient for Delta
    439521//----------------------------------------------------------
    … …  
    452534      Ptarget.setE(std::sqrt(M0target2+PZcms2));
    453535
    454       {
     536// ----------------------------------------------------------
     537
     538//      G4double Wexcit=1.-1.97*std::exp(-0.5*ProjectileRapidity);
     539      G4double Wexcit=1.-2.256*std::exp(-0.6*ProjectileRapidity);
     540
     541//G4cout<<ProjectileRapidity<<" "<<1.72*std::exp(-0.4*ProjectileRapidity)<<" "<<std::exp(0.4*ProjectileRapidity)<<G4endl;
     542//G4int Uzhi;G4cin>>Uzhi;
     543//Wexcit=0.;
     544      if(G4UniformRand() > Wexcit)
     545      {                             // Make elastic scattering
     546//G4cout<<"Make elastic scattering"<<G4endl;
    455547       Pprojectile.transform(toLab);
    456548       Ptarget.transform(toLab);
    … …  
    463555
    464556       G4bool Result= theElastic->ElasticScattering (projectile,target,theParameters);
    465 
    466557       return Result;
    467       }
     558      } // end of if(G4UniformRand() > Wexcit)
    468559     }  // End of charge exchange part ------------------------------
    469560
    470561// ------------------------------------------------------------------
     562     G4double ProbOfDiffraction=ProbProjectileDiffraction + ProbTargetDiffraction;
     563/*
     564G4cout<<"Excite --------------------"<<G4endl;
     565G4cout<<"Proj "<<M0projectile<<" "<<ProjectileDiffStateMinMass<<"  "<<ProjectileNonDiffStateMinMass<<G4endl;
     566G4cout<<"Targ "<<M0target    <<" "<<TargetDiffStateMinMass    <<" "<<TargetNonDiffStateMinMass<<G4endl;
     567G4cout<<"SqrtS "<<SqrtS<<G4endl;
     568
     569G4cout<<"Prob ProjDiff TargDiff "<<ProbProjectileDiffraction<<" "<<ProbTargetDiffraction<<" "<<ProbOfDiffraction<<G4endl;
     570G4cout<<"Pr Y "<<Pprojectile.rapidity()<<" Tr Y "<<Ptarget.rapidity()<<G4endl;
     571//G4int Uzhi; G4cin>>Uzhi;
     572*/
     573/*
     574     if(ProjectileNonDiffStateMinMass + TargetNonDiffStateMinMass > SqrtS) // 24.07.10
     575     {
     576      if(ProbOfDiffraction!=0.)
     577      {
     578       ProbProjectileDiffraction/=ProbOfDiffraction;
     579       ProbOfDiffraction=1.;
     580      } else {return false;}     
     581     }
     582
     583*/
     584
    471585     if(ProbOfDiffraction!=0.)
    472586     {
    … …  
    478592     }
    479593
     594//G4cout<<"Prob ProjDiff TargDiff "<<ProbProjectileDiffraction<<" "<<ProbTargetDiffraction<<" "<<ProbOfDiffraction<<G4endl;
     595
    480596     G4double ProjectileDiffStateMinMass2    = ProjectileDiffStateMinMass    *
    481597                                               ProjectileDiffStateMinMass;
    … …  
    500616
    501617     G4int whilecount=0;
     618
    502619//   Choose a process ---------------------------
    503620
    … …  
    506623        if(G4UniformRand() < ProbProjectileDiffraction)
    507624        { //-------- projectile diffraction ---------------
     625//G4cout<<"projectile diffraction"<<G4endl;
     626
    508627         do {
    509628//             Generate pt
    … …  
    512631//               << ", loop count/ maxPtSquare : "
    513632//               << whilecount << " / " << maxPtSquare << G4endl;
     633
     634//             whilecount++;
    514635             if (whilecount > 1000 )
    515636             {
    … …  
    517638              target->SetStatus(2);  return false;    //  Ignore this interaction
    518639             };
     640
    519641// --------------- Check that the interaction is possible -----------
    520642             ProjMassT2=ProjectileDiffStateMinMass2;
    … …  
    523645             TargMassT2=M0target2;
    524646             TargMassT =M0target;
    525 
     647//G4cout<<"Masses "<<ProjMassT<<" "<<TargMassT<<" "<<SqrtS<<" "<<ProjMassT+TargMassT<<G4endl;
    526648             PZcms2=(S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2-
    527649                     2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
    528650                    /4./S;
    529651
     652//G4cout<<"PZcms2 PrD"<<PZcms2<<G4endl;
    530653             if(PZcms2 < 0 )
    531654             {
    … …  
    556679             PMinusNew=ChooseP(PMinusMin, PMinusMax);
    557680// PMinusNew=1./sqrt(1./PMinusMin-G4UniformRand()*(1./PMinusMin-1./PMinusMax));
     681//PMinusNew=1./sqr(1./std::sqrt(PMinusMin)-G4UniformRand()*(1./std::sqrt(PMinusMin)-1./std::sqrt(PMinusMax)));
    558682
    559683             TMinusNew=SqrtS-PMinusNew;
    … …  
    564688             Qmomentum.setPz( (Qplus-Qminus)/2 );
    565689             Qmomentum.setE(  (Qplus+Qminus)/2 );
    566           } while (
    567 ((Pprojectile+Qmomentum).mag2() <  ProjectileDiffStateMinMass2) ||  //No without excitation
    568 ((Ptarget    -Qmomentum).mag2() <  M0target2                  ));
     690
     691          } while ((Pprojectile+Qmomentum).mag2() <  ProjectileDiffStateMinMass2); //|| 
     692                  //Repeat the sampling because there was not any excitation
     693//((Ptarget    -Qmomentum).mag2() <  M0target2                  )) );
    569694        }
    570695        else
    571696        { // -------------- Target diffraction ----------------
     697
     698//G4cout<<"Target diffraction"<<G4endl;
    572699         do {
    573700//             Generate pt
    … …  
    576703//               << ", loop count/ maxPtSquare : "
    577704//               << whilecount << " / " << maxPtSquare << G4endl;
     705
     706//             whilecount++;
    578707             if (whilecount > 1000 )
    579708             {
    … …  
    581710              target->SetStatus(2);  return false;    //  Ignore this interaction
    582711             };
     712//G4cout<<"Qm while "<<Qmomentum<<" "<<whilecount<<G4endl;
    583713// --------------- Check that the interaction is possible -----------
    584714             ProjMassT2=M0projectile2;
    … …  
    592722                    /4./S;
    593723
     724//G4cout<<"PZcms2 TrD <0 "<<PZcms2<<" return"<<G4endl;
    594725             if(PZcms2 < 0 )
    595726             {
    … …  
    600731
    601732             Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
     733
     734//G4cout<<"Qm while "<<Qmomentum<<" "<<whilecount<<G4endl;
    602735             Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
    603736
    … …  
    612745                    /4./S;
    613746
     747//G4cout<<"PZcms2 <0 "<<PZcms2<<" continue"<<G4endl;
    614748             if(PZcms2 < 0 ) continue;
    615749             PZcms =std::sqrt(PZcms2);
    … …  
    619753
    620754             TPlusNew=ChooseP(TPlusMin, TPlusMax);
     755//TPlusNew=1./sqr(1./std::sqrt(TPlusMin)-G4UniformRand()*(1./std::sqrt(TPlusMin)-1./std::sqrt(TPlusMax)));
    621756
    622757//TPlusNew=TPlusMax;
    … …  
    630765             Qmomentum.setE(  (Qplus+Qminus)/2 );
    631766
    632           } while (
    633  ((Pprojectile+Qmomentum).mag2() <  M0projectile2          ) ||  //No without excitation
    634  ((Ptarget    -Qmomentum).mag2() <  TargetDiffStateMinMass2));
    635          }
     767/*
     768G4cout<<(Pprojectile+Qmomentum).mag()<<" "<<M0projectile<<G4endl;
     769G4bool First=(Pprojectile+Qmomentum).mag2() <  M0projectile2;
     770G4cout<<First<<G4endl;
     771
     772G4cout<<(Ptarget    -Qmomentum).mag()<<" "<<TargetDiffStateMinMass<<" "<<TargetDiffStateMinMass2<<G4endl;
     773G4bool Seco=(Ptarget    -Qmomentum).mag2() < TargetDiffStateMinMass2;
     774G4cout<<Seco<<G4endl;
     775*/
     776
     777         } while ((Ptarget    -Qmomentum).mag2() <  TargetDiffStateMinMass2);
     778                 // Repeat the sampling because there was not any excitation
     779// (((Pprojectile+Qmomentum).mag2() <  M0projectile2          ) ||  //No without excitation
     780//  ((Ptarget    -Qmomentum).mag2() <  TargetDiffStateMinMass2)) );
     781//G4cout<<"Go out"<<G4endl;
     782         } // End of if(G4UniformRand() < ProbProjectileDiffraction)
    636783        }
    637784        else  //----------- Non-diffraction process ------------
    638785        {
     786
     787//G4cout<<"Non-diffraction process"<<G4endl;
    639788         do {
    640789//             Generate pt
    … …  
    643792//               << ", loop count/ maxPtSquare : "
    644793//               << whilecount << " / " << maxPtSquare << G4endl;
     794
     795//             whilecount++;
    645796             if (whilecount > 1000 )
    646797             {
    … …  
    677828                    2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)
    678829                   /4./S;
     830//G4cout<<"PZcms2 ND"<<PZcms2<<G4endl;
    679831
    680832             if(PZcms2 < 0 ) continue;
    … …  
    698850             Qmomentum.setPz( (Qplus-Qminus)/2 );
    699851             Qmomentum.setE(  (Qplus+Qminus)/2 );
    700 
     852/*
     853G4cout<<(Pprojectile+Qmomentum).mag2()<<" "<<ProjectileNonDiffStateMinMass2<<G4endl;
     854G4cout<<(Ptarget    -Qmomentum).mag2()<<" "<<TargetNonDiffStateMinMass2<<G4endl;
     855G4int Uzhi; G4cin>>Uzhi;
     856*/
    701857       } while (
    702858 ((Pprojectile+Qmomentum).mag2() <  ProjectileNonDiffStateMinMass2) || //No double Diffraction
    703859 ((Ptarget    -Qmomentum).mag2() <  TargetNonDiffStateMinMass2    ));
    704          }
    705 
    706            Pprojectile += Qmomentum;
    707            Ptarget     -= Qmomentum;
     860     }
     861
     862     Pprojectile += Qmomentum;
     863     Ptarget     -= Qmomentum;
     864
     865//G4cout<<"Pr Y "<<Pprojectile.rapidity()<<" Tr Y "<<Ptarget.rapidity()<<G4endl;
    708866
    709867// Transform back and update SplitableHadron Participant.
    710            Pprojectile.transform(toLab);
    711            Ptarget.transform(toLab);
     868     Pprojectile.transform(toLab);
     869     Ptarget.transform(toLab);
    712870
    713871// Calculation of the creation time ---------------------
    714       projectile->SetTimeOfCreation(target->GetTimeOfCreation());
    715       projectile->SetPosition(target->GetPosition());
     872     projectile->SetTimeOfCreation(target->GetTimeOfCreation());
     873     projectile->SetPosition(target->GetPosition());
    716874// Creation time and position of target nucleon were determined at
    717875// ReggeonCascade() of G4FTFModel
    718876// ------------------------------------------------------
    719877
    720            projectile->Set4Momentum(Pprojectile);
    721            target->Set4Momentum(Ptarget);
    722 
    723            projectile->IncrementCollisionCount(1);
    724            target->IncrementCollisionCount(1);
    725 
    726            return true;
     878//G4cout<<"Mproj "<<Pprojectile.mag()<<G4endl;
     879//G4cout<<"Mtarg "<<Ptarget.mag()<<G4endl;
     880     projectile->Set4Momentum(Pprojectile);
     881     target->Set4Momentum(Ptarget);
     882
     883     projectile->IncrementCollisionCount(1);
     884     target->IncrementCollisionCount(1);
     885
     886     return true;
    727887}
    728888

  • trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4DiffractiveSplitableHadron.cc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4DiffractiveSplitableHadron.cc,v 1.8 2009/07/31 11:03:00 vuzhinsk Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4DiffractiveSplitableHadron.cc,v 1.9 2010/09/20 15:50:46 vuzhinsk Exp $
     28// GEANT4 tag $Name: geant4-09-03-ref-09 $
    2929//
    3030

  • trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4FTFModel.cc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4FTFModel.cc,v 1.34 2009/12/15 19:14:31 vuzhinsk Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4FTFModel.cc,v 1.36 2010/09/20 15:50:46 vuzhinsk Exp $
     28// GEANT4 tag $Name: geant4-09-03-ref-09 $
    2929//
    3030
    … …  
    9999{
    100100        theProjectile = aProjectile; 
    101         theParticipants.Init(aNucleus.GetN(),aNucleus.GetZ());
     101
     102        theParticipants.Init(aNucleus.GetA_asInt(),aNucleus.GetZ_asInt());
     103
    102104// ----------- N-mass number Z-charge -------------------------
    103105
    … …  
    113115
    114116//theParameters->SetProbabilityOfElasticScatt(0.);
     117//G4cout<<theParameters->GetProbabilityOfElasticScatt()<<G4endl;
     118//G4int Uzhi; G4cin>>Uzhi;
    115119// To turn on/off (1/0) elastic scattering
    116120
    … …  
    125129{
    126130        G4ExcitedStringVector * theStrings(0);
    127 
     131//G4cout<<"GetString"<<G4endl;
    128132        theParticipants.GetList(theProjectile,theParameters);
    129 
     133//G4cout<<"Reggeon"<<G4endl;
    130134        ReggeonCascade();
    131135
    … …  
    133137        if( PutOnMassShell() )
    134138        {
     139//G4cout<<"PutOn mass Shell OK"<<G4endl;
    135140         if( ExciteParticipants() )
    136141         {
     142//G4cout<<"Excite partic OK"<<G4endl;
    137143          theStrings = BuildStrings();
    138 
     144//G4cout<<"Build String OK"<<G4endl;
    139145          GetResidualNucleus();
    140146
    … …  
    194200           TheInvolvedNucleon[NumberOfInvolvedNucleon]=TargetNucleon;
    195201           NumberOfInvolvedNucleon++;
    196 
     202//G4cout<<"Prim NumberOfInvolvedNucleon "<<NumberOfInvolvedNucleon<<G4endl;
    197203           G4double XofWoundedNucleon = TargetNucleon->GetPosition().x();
    198204           G4double YofWoundedNucleon = TargetNucleon->GetPosition().y();
    … …  
    214220              TheInvolvedNucleon[NumberOfInvolvedNucleon]=Neighbour;
    215221              NumberOfInvolvedNucleon++;
     222//G4cout<<"Seco NumberOfInvolvedNucleon "<<NumberOfInvolvedNucleon<<G4endl;
    216223
    217224              G4VSplitableHadron *targetSplitable;
    … …  
    274281        G4LorentzVector Pprojectile=primary->Get4Momentum();
    275282
     283//G4cout<<"Pprojectile "<<Pprojectile<<G4endl;
    276284// To get original projectile particle
    277285
    … …  
    284292        G4double        SumMasses = Mprojectile + 20.*MeV; // 13.12.09
    285293                                               // Separation energy for projectile
    286 
     294//G4cout<<"SumMasses Pr "<<SumMasses<<G4endl;
    287295//--------------- Target nucleus ------------------------------
    288296        G4V3DNucleus *theNucleus = GetWoundedNucleus();
    … …  
    305313          SumMasses += aNucleon->GetDefinition()->GetPDGMass();
    306314          SumMasses += 20.*MeV;   // 13.12.09 Separation energy for a nucleon
     315//G4cout<<"SumMasses Tr "<<SumMasses<<G4endl;
    307316          ResidualMassNumber--;
    308317          ResidualCharge-=(G4int) aNucleon->GetDefinition()->GetPDGCharge();
    … …  
    316325
    317326        Psum += PnuclearResidual;
    318 
     327//G4cout<<"ResidualCharge ,ResidualMassNumber "<<ResidualCharge<<" "<<ResidualMassNumber<<G4endl;
    319328        G4double ResidualMass(0.);
    320329        if(ResidualMassNumber == 0)
    … …  
    331340 
    332341//      ResidualMass +=ResidualExcitationEnergy; // Will be given after checks
     342//G4cout<<"SumMasses End ResidualMass "<<SumMasses<<" "<<ResidualMass<<G4endl;
    333343        SumMasses += ResidualMass;
    334 
     344//G4cout<<"SumMasses + ResM "<<SumMasses<<G4endl;
     345//G4cout<<"Psum "<<Psum<<G4endl;
    335346//-------------------------------------------------------------
    336347
    … …  
    338349        G4double     S=Psum.mag2();
    339350
     351//G4cout<<"SqrtS < SumMasses "<<SqrtS<<" "<<SumMasses<<G4endl;
    340352        if(SqrtS < SumMasses)      {return false;} // It is impossible to simulate
    341353                                                   // after putting nuclear nucleons
    … …  
    346358        ResidualMass +=ResidualExcitationEnergy;
    347359        SumMasses    +=ResidualExcitationEnergy;
    348 
     360//G4cout<<"ResidualMass "<<ResidualMass<<" "<<SumMasses<<G4endl;
    349361//-------------------------------------------------------------
    350362// Sampling of nucleons what are transfered to delta-isobars --
    … …  
    373385        }   // end of if(theNucleus.GetMassNumber() != 1)
    374386//-------------------------------------------------------------
     387
    375388        G4LorentzRotation toCms(-1*Psum.boostVector());
    376389        G4LorentzVector Ptmp=toCms*Pprojectile;
    … …  
    393406        G4double AveragePt2  = theParameters->GetPt2ofNuclearDestruction();
    394407        G4double maxPtSquare = theParameters->GetMaxPt2ofNuclearDestruction();
    395 
     408//G4cout<<"Dcor "<<Dcor<<" AveragePt2 "<<AveragePt2<<G4endl;
    396409        G4double M2target(0.);
    397410        G4double WminusTarget(0.);
    … …  
    409422
    410423            NumberOfTries++;
     424//G4cout<<"NumberOfTries "<<NumberOfTries<<G4endl;
    411425            if(NumberOfTries == 100*(NumberOfTries/100))   // 100
    412426            { // At large number of tries it would be better to reduce the values
    … …  
    439453
    440454               G4LorentzVector tmp(tmpPt.x(),tmpPt.y(),Xminus,0.);
     455//G4cout<<"Inv i mom "<<i<<" "<<tmp<<G4endl;
    441456               aNucleon->SetMomentum(tmp);
    442457             }   // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
    … …  
    447462             G4double DeltaXminus(0.);
    448463
     464//G4cout<<"ResidualMassNumber "<<ResidualMassNumber<<" "<<PtSum<<G4endl;
    449465             if(ResidualMassNumber == 0)
    450466             {
    … …  
    457473              DeltaXminus = -1./theNucleus->GetMassNumber();
    458474             }
    459 
     475//G4cout<<"Dx y xmin "<<DeltaX<<" "<<DeltaY<<" "<<DeltaXminus<<G4endl;
    460476             XminusSum=1.;
    461477             M2target =0.;
    … …  
    467483               Xminus = aNucleon->Get4Momentum().pz() - DeltaXminus;
    468484               XminusSum-=Xminus;               
    469 
    470                if((Xminus <= 0.)   || (Xminus > 1.) ||
    471                   (XminusSum <=0.) || (XminusSum > 1.)) {InerSuccess=false; break;}
    472  
     485//G4cout<<" i X-sum "<<i<<" "<<Xminus<<" "<<XminusSum<<G4endl;
     486               if(ResidualMassNumber == 0)                // Uzhi 5.07.10
     487               {
     488                if((Xminus <= 0.)   || (Xminus > 1.))    {InerSuccess=false; break;}
     489               } else
     490               {
     491                if((Xminus <= 0.)   || (Xminus > 1.) ||
     492                   (XminusSum <=0.) || (XminusSum > 1.)) {InerSuccess=false; break;}
     493               }                                          // Uzhi 5.07.10
     494
    473495               G4double Px=aNucleon->Get4Momentum().px() - DeltaX;
    474496               G4double Py=aNucleon->Get4Momentum().py() - DeltaY;
    … …  
    482504             }   // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
    483505
     506//G4cout<<"Rescale O.K."<<G4endl;
     507
    484508             if(InerSuccess && (ResidualMassNumber != 0))
    485509             {
    486510              M2target +=(ResidualMass*ResidualMass + PtSum.mag2())/XminusSum;
    487511             }
     512//G4cout<<"InerSuccess "<<InerSuccess<<G4endl;
     513//G4int Uzhi;G4cin>>Uzhi;
    488514            } while(!InerSuccess);
    489515          } while (SqrtS < Mprojectile + std::sqrt(M2target));
    … …  
    495521          WminusTarget=(S-M2projectile+M2target+std::sqrt(DecayMomentum2))/2./SqrtS;
    496522          WplusProjectile=SqrtS - M2target/WminusTarget;
     523//G4cout<<"DecayMomentum2 "<<DecayMomentum2<<G4endl;
    497524//-------------------------------------------------------------
    498525          for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
    … …  
    511538           if( E+Pz > WplusProjectile ){OuterSuccess=false; break;}
    512539          }   // end of for(G4int i=0; i < NumberOfInvolvedNucleon; i++ )
     540//G4int Uzhi;G4cin>>Uzhi;
    513541        } while(!OuterSuccess);
    514542
    … …  
    575603        Successfull=false;
    576604        theParticipants.StartLoop();
     605
     606G4int MaxNumOfInelCollisions=G4int(theParameters->GetMaxNumberOfCollisions());
     607G4double NumberOfInel(0.);
     608//
     609if(MaxNumOfInelCollisions > 0) 
     610{   //  Plab > Pbound, Normal application of FTF is possible
     611 G4double ProbMaxNumber=theParameters->GetMaxNumberOfCollisions()-MaxNumOfInelCollisions;
     612 if(G4UniformRand() < ProbMaxNumber) {MaxNumOfInelCollisions++;}
     613 NumberOfInel=MaxNumOfInelCollisions;
     614} else
     615{   //  Plab < Pbound, Normal application of FTF is impossible, low energy corrections
     616 if(theParticipants.theNucleus->GetMassNumber() > 1)
     617 {
     618  NumberOfInel = theParameters->GetProbOfInteraction();
     619  MaxNumOfInelCollisions = 1;
     620 } else
     621 { // Special case for hadron-nucleon interactions
     622  NumberOfInel = 1.;
     623  MaxNumOfInelCollisions = 1;
     624 }
     625}  // end of if(MaxNumOfInelCollisions > 0)
     626//
    577627        while (theParticipants.Next())
    578628        {         
    … …  
    581631           G4VSplitableHadron * projectile=collision.GetProjectile();
    582632           G4VSplitableHadron * target=collision.GetTarget();
    583 
     633//G4cout<<"ProbabilityOfElasticScatt "<<theParameters->GetProbabilityOfElasticScatt()<<G4endl;
    584634           if(G4UniformRand()< theParameters->GetProbabilityOfElasticScatt())
    585635           { //   Elastic scattering -------------------------
     636//G4cout<<"Elastic FTF"<<G4endl;
    586637            if(theElastic->ElasticScattering(projectile, target, theParameters))
    587638            {
    … …  
    595646           else
    596647           { //   Inelastic scattering ----------------------
     648/*
    597649            if(theExcitation->ExciteParticipants(projectile, target,
    598650                                                 theParameters, theElastic))
    … …  
    604656             target->SetStatus(2);
    605657            }
     658*/
     659//G4cout<<"InElastic FTF"<<G4endl;
     660            if(G4UniformRand()< NumberOfInel/MaxNumOfInelCollisions)
     661            {
     662             if(theExcitation->ExciteParticipants(projectile, target,
     663                                                 theParameters, theElastic))
     664             {
     665              Successfull = Successfull || true;
     666NumberOfInel--;
     667             } else
     668             {
     669              Successfull = Successfull || false;
     670              target->SetStatus(2);
     671             }
     672            } else // If NumOfInel
     673            {
     674             if(theElastic->ElasticScattering(projectile, target, theParameters))
     675             {
     676              Successfull = Successfull || true;
     677             } else
     678             {
     679              Successfull = Successfull || false;
     680              target->SetStatus(2);
     681             }
     682            }   // end if NumOfInel
    606683           }
    607684        }       // end of while (theParticipants.Next())
    … …  
    624701
    625702        theParticipants.StartLoop();    // restart a loop
    626 
     703//
    627704        while ( theParticipants.Next() )
    628705        {
    … …  
    634711                primaries.push_back(interaction.GetProjectile());     
    635712        }
    636            
     713
    637714        unsigned int ahadron;
    638715        for ( ahadron=0; ahadron < primaries.size() ; ahadron++)
    … …  
    650727            if(SecondString != 0) strings->push_back(SecondString);
    651728        }
    652 
     729//
    653730        for (G4int ahadron=0; ahadron < NumberOfInvolvedNucleon ; ahadron++)
    654731        {
    … …  
    669746        std::for_each(primaries.begin(), primaries.end(), DeleteVSplitableHadron());
    670747        primaries.clear();
    671 
    672748        return strings;
    673749}

  • trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4FTFParameters.cc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4FTFParameters.cc,v 1.13 2009/12/16 17:51:15 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4FTFParameters.cc,v 1.14 2010/09/20 15:50:46 vuzhinsk Exp $
     28// GEANT4 tag $Name: geant4-09-03-ref-09 $
    2929//
    3030
    … …  
    5656    G4double Plab = std::sqrt(Elab * Elab - ProjectileMass*ProjectileMass);
    5757
     58    G4double Ylab=0.5*std::log((Elab+Plab)/(Elab-Plab));
     59
     60    Plab/=GeV;                               // Uzhi 8.07.10
    5861    G4double LogPlab = std::log( Plab );
    5962    G4double sqrLogPlab = LogPlab * LogPlab;
    … …  
    186189      SetInelasticCrossSection(Xtotal-Xelastic);
    187190
     191//G4cout<<"Xtotal, Xelastic "<<Xtotal<<" "<<Xelastic<<G4endl;
    188192//  // Interactions with elastic and inelastic collisions
    189193      SetProbabilityOfElasticScatt(Xtotal, Xelastic);
    … …  
    196200
    197201//----------------------------------------------------------------------------------- 
     202
    198203      SetSlope( Xtotal*Xtotal/16./pi/Xelastic/0.3894 ); // Slope parameter of elastic scattering
    199204                                                        //      (GeV/c)^(-2))
    … …  
    209214           if( absPDGcode > 1000 )                        //------Projectile is baryon --------
    210215             {
    211               SetMagQuarkExchange(3.4); //3.8);
    212               SetSlopeQuarkExchange(1.2);
    213               SetDeltaProbAtQuarkExchange(0.1); //(0.1*4.);
    214 
    215               SetProjMinDiffMass(1.1);                    // GeV
    216               SetProjMinNonDiffMass(1.1);                 // GeV
    217               SetProbabilityOfProjDiff(0.76*std::pow(s/GeV/GeV,-0.35));
    218 
    219               SetTarMinDiffMass(1.1);                     // GeV
    220               SetTarMinNonDiffMass(1.1);                  // GeV
    221               SetProbabilityOfTarDiff(0.76*std::pow(s/GeV/GeV,-0.35));
    222 
    223               SetAveragePt2(0.3);                         // GeV^2
     216              SetMagQuarkExchange(1.84);//(3.63);
     217              SetSlopeQuarkExchange(0.7);//(1.2);
     218              SetDeltaProbAtQuarkExchange(0.);
     219
     220              SetProjMinDiffMass(1.16);                   // GeV
     221              SetProjMinNonDiffMass(1.16);                // GeV
     222             
     223SetProbabilityOfProjDiff(0.805*std::exp(-0.35*Ylab));// 0.5
     224
     225              SetTarMinDiffMass(1.16);                    // GeV
     226              SetTarMinNonDiffMass(1.16);                 // GeV
     227             
     228SetProbabilityOfTarDiff(0.805*std::exp(-0.35*Ylab));// 0.5
     229
     230              SetAveragePt2(0.15);                        // 0.15 GeV^2
    224231             }
    225232           else if( absPDGcode == 211 || PDGcode ==  111) //------Projectile is Pion -----------
    … …  
    286293    if( absPDGcode < 1000 )
    287294    {
    288       SetCofNuclearDestruction(1.); //1.0);                 // for meson projectile
    289     } else if( theA > 20. )
     295      SetMaxNumberOfCollisions(1000.,1.); //(Plab,2.); //3.); ##############################
     296
     297      SetCofNuclearDestruction(0.3); //1.0);           // for meson projectile
     298
     299      SetDofNuclearDestruction(0.4);
     300      SetPt2ofNuclearDestruction(0.17*GeV*GeV);
     301      SetMaxPt2ofNuclearDestruction(1.0*GeV*GeV);
     302
     303      SetExcitationEnergyPerWoundedNucleon(100*MeV);
     304    } else                                             // for baryon projectile
    290305    {
    291       SetCofNuclearDestruction(0.2); //2);                 // for baryon projectile and heavy target
    292     } else
    293     {
    294       SetCofNuclearDestruction(0.2); //1.0);                 // for baryon projectile and light target
     306//      SetMaxNumberOfCollisions(Plab,0.1); //6.); // ##############################
     307      SetMaxNumberOfCollisions(Plab,4.); //6.); // ##############################
     308
     309      SetCofNuclearDestruction(0.62*std::exp(4.*(Ylab-2.1))/(1.+std::exp(4.*(Ylab-2.1))));
     310
     311      SetDofNuclearDestruction(0.4);
     312      SetPt2ofNuclearDestruction((0.035+
     313          0.04*std::exp(4.*(Ylab-2.5))/(1.+std::exp(4.*(Ylab-2.5))))*GeV*GeV); //0.09
     314      SetMaxPt2ofNuclearDestruction(1.0*GeV*GeV);
     315
     316      SetExcitationEnergyPerWoundedNucleon(75.*MeV);
    295317    }
    296318
    297319    SetR2ofNuclearDestruction(1.5*fermi*fermi);
    298320
    299     SetExcitationEnergyPerWoundedNucleon(100*MeV);
    300 
    301     SetDofNuclearDestruction(0.4);
    302     SetPt2ofNuclearDestruction(0.17*GeV*GeV);
    303     SetMaxPt2ofNuclearDestruction(1.0*GeV*GeV);
     321//SetCofNuclearDestruction(0.47*std::exp(2.*(Ylab-2.5))/(1.+std::exp(2.*(Ylab-2.5))));
     322//SetPt2ofNuclearDestruction((0.035+0.1*std::exp(4.*(Ylab-3.))/(1.+std::exp(4.*(Ylab-3.))))*GeV*GeV);
     323
     324//SetProbabilityOfElasticScatt(1.,1.); //(Xtotal, Xelastic);
     325//SetProbabilityOfProjDiff(1.*0.62*std::pow(s/GeV/GeV,-0.51)); // 0->1
     326//SetProbabilityOfTarDiff(4.*0.62*std::pow(s/GeV/GeV,-0.51)); // 2->4
     327//SetAveragePt2(0.3);                              //(0.15);
     328//SetAvaragePt2ofElasticScattering(0.);
     329
     330//SetCofNuclearDestruction(0.6); //(0.4);                 
     331SetExcitationEnergyPerWoundedNucleon(75.*MeV); //(75.*MeV);
     332//SetDofNuclearDestruction(0.6); //(0.4);                 
     333//SetPt2ofNuclearDestruction(0.12*GeV*GeV); //(0.168*GeV*GeV);
     334
    304335}
    305336//**********************************************************************************************

  • trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4FTFParticipants.cc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4FTFParticipants.cc,v 1.16 2009/11/25 09:14:03 vuzhinsk Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4FTFParticipants.cc,v 1.17 2010/09/20 15:50:46 vuzhinsk Exp $
     28// GEANT4 tag $Name: geant4-09-03-ref-09 $
    2929//
    3030// ------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.