Ignore:
Timestamp:
Sep 30, 2010, 2:47:17 PM (14 years ago)
Author:
garnier
Message:

tag geant4.9.4 beta 1 + modifs locales

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/models/parton_string/hadronization/src/G4LundStringFragmentation.cc

    r1228 r1337  
    2525//
    2626//
    27 // $Id: G4LundStringFragmentation.cc,v 1.18 2009/10/05 12:52:48 vuzhinsk Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $ 1.8
     27// $Id: G4LundStringFragmentation.cc,v 1.20 2010/06/21 17:50:48 vuzhinsk Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-01 $ 1.8
    2929//
    3030// -----------------------------------------------------------------------------
    … …  
    5858    SmoothParam  = 0.2;                   
    5959
    60 //    SetStringTensionParameter(0.25);                           // Uzhi 20 June 08
    61     SetStringTensionParameter(1.);                           // Uzhi 20 June 09
     60//    SetStringTensionParameter(0.25);                           
     61    SetStringTensionParameter(1.);                         
    6262   }
    6363
    … …  
    9090
    9191//--------------------------------------------------------------------------------------
    92 void G4LundStringFragmentation::SetMinimalStringMass(const G4FragmentingString  * const string)  // Uzhi
    93 {
    94 /*
    95 G4cout<<"In SetMinMass -------------------"<<std::sqrt(string->Mass2())<<G4endl;
    96 G4cout<<string->GetLeftParton()->GetPDGEncoding()<<" "<<
    97         string->GetRightParton()->GetPDGEncoding()<<G4endl;
    98 */
     92void G4LundStringFragmentation::SetMinimalStringMass(const G4FragmentingString  * const string) 
     93{
    9994  G4double EstimatedMass=0.; 
    10095  G4int Number_of_quarks=0;
    … …  
    150145  MinimalStringMass=EstimatedMass;
    151146  SetMinimalStringMass2(EstimatedMass);
    152 //G4cout<<"Out SetMinimalStringMass "<<MinimalStringMass<<G4endl;
    153147}
    154148
    … …  
    164158                                const G4ExcitedString& theString)
    165159{
    166 
    167 //G4cout<<"In FragmentString"<<G4endl;
    168 
    169160//    Can no longer modify Parameters for Fragmentation.
    170161        PastInitPhase=true;
    171        
     162
    172163//      check if string has enough mass to fragment...
    173164       
    … …  
    180171G4cout<<"FragmentString Momentum"<<theString.Get4Momentum()<<theString.Get4Momentum().mag()<<G4endl;
    181172*/
     173        G4FragmentingString aString(theString);
     174        SetMinimalStringMass(&aString);
     175
     176       if((MinimalStringMass+WminLUND)*(1.-SmoothParam) > theString.Get4Momentum().mag())
     177       {SetMassCut(1000.*MeV);}
     178// V.U. 20.06.10 in order to put un correspondence LightFragTest and MinStrMass
     179
    182180        G4KineticTrackVector * LeftVector=LightFragmentationTest(&theString);
     181
    183182        if ( LeftVector != 0 ) {
    184 //G4cout<<"Return single hadron insted of string"<<G4endl;
    185183// Uzhi insert 6.05.08 start
    186184          if(LeftVector->size() == 1){
    … …  
    197195            LeftVector->operator[](0)->SetPosition(aPosition);
    198196*/           
    199 //G4cout<<"Single hadron "<<LeftVector->operator[](0)->Get4Momentum()<<" "<<LeftVector->operator[](0)->Get4Momentum().mag()<<G4endl;
    200197          } else {    // 2 hadrons created from qq-qqbar are stored
     198
    201199            LeftVector->operator[](0)->SetFormationTime(theString.GetTimeOfCreation());
    202200            LeftVector->operator[](0)->SetPosition(theString.GetPosition());
    … …  
    204202            LeftVector->operator[](1)->SetPosition(theString.GetPosition());
    205203          }
    206 // Uzhi insert 6.05.08 end
    207204          return LeftVector;
    208205        }
    … …  
    210207//--------------------- The string can fragment -------------------------------
    211208//--------------- At least two particles can be produced ----------------------
    212 
    213209                               LeftVector =new G4KineticTrackVector;
    214210        G4KineticTrackVector * RightVector=new G4KineticTrackVector;
    … …  
    223219                G4FragmentingString *currentString=new G4FragmentingString(*theStringInCMS);
    224220
    225 //G4cout<<"FragmentString cur M2  "<<std::sqrt(currentString->Mass2())<<G4endl;
    226221          // Cleaning up the previously produced hadrons ------------------------------
    227222                std::for_each(LeftVector->begin() , LeftVector->end() , DeleteKineticTrack());
    … …  
    237232                        G4FragmentingString *newString=0;  // used as output from SplitUp...
    238233
    239 //G4cout<<"FragmentString to Splitup ===================================="<<G4endl;
    240 //G4cout<<"++++++++++++++++++++++++++ Enter num--------------------------"<<G4endl;
    241 //G4int Uzhi; G4cin>>Uzhi;                         // Uzhi
    242 
    243234                        G4KineticTrack * Hadron=Splitup(currentString,newString);
    244 //G4cout<<" Hadron "<<Hadron<<G4endl;
    245235
    246236                        if ( Hadron != 0 )  // Store the hadron                               
    … …  
    255245                };
    256246        // Split remaining string into 2 final Hadrons ------------------------
    257 //G4cout<<"FragmentString to SplitLast if inner_sucess#0"<<inner_sucess<<G4endl;
    258247                if ( inner_sucess &&                                       
    259248                     SplitLast(currentString,LeftVector, RightVector) )
    … …  
    283272        delete RightVector;
    284273
    285 //G4cout<<"CalculateHadronTimePosition"<<G4endl;
    286 
    287274        CalculateHadronTimePosition(theString.Get4Momentum().mag(), LeftVector);
    288275
    … …  
    294281        G4double      TimeOftheStringCreation=theString.GetTimeOfCreation();
    295282        G4ThreeVector PositionOftheStringCreation(theString.GetPosition());
    296 // Uzhi 11.07.09
    297 //G4cout<<" String T Pos"<<TimeOftheStringCreation<<' '<<PositionOftheStringCreation<<G4endl;
    298 /*  // For large formation time open *
    299         G4double      TimeOftheStringCreation=theString.GetTimeOfCreation()+100*fermi;
    300         G4ThreeVector PositionOftheStringCreation(theString.GetPosition().x(),
    301                                                   theString.GetPosition().y(),
    302                                                   theString.GetPosition().z()+100*fermi);
    303 */
    304 
    305 /*
    306         if(theString.GetPosition().y() > 100.*fermi){   
    307 // It is a projectile-like string -------------------------------------
    308           G4double Zmin=theString.GetPosition().y()-1000.*fermi;
    309           G4double Zmax=theString.GetPosition().z();
    310           TimeOftheStringCreation=
    311           (Zmax-Zmin)*theString.Get4Momentum().e()/theString.Get4Momentum().z();
    312 
    313           G4ThreeVector aPosition(0.,0.,Zmax);
    314           PositionOftheStringCreation=aPosition;
    315         }
    316 */
    317 
    318 //G4cout<<"Lund Frag # hadrons"<<LeftVector->size()<<G4endl;       // Uzhi 11.07.09
     283
    319284        for(size_t C1 = 0; C1 < LeftVector->size(); C1++)
    320285        {
    … …  
    326291           G4LorentzVector Coordinate(Hadron->GetPosition(), Hadron->GetFormationTime());
    327292           Momentum = toObserverFrame*Coordinate;
    328            Hadron->SetFormationTime(TimeOftheStringCreation+Momentum.e()    // Uzhi 11.07.09
    329                                                            -fermi/c_light); // Uzhi 11.07.09
     293           Hadron->SetFormationTime(TimeOftheStringCreation+Momentum.e()   
     294                                                           -fermi/c_light);
    330295           G4ThreeVector aPosition(Momentum.vect());
    331296//         Hadron->SetPosition(theString.GetPosition()+aPosition);
    332297           Hadron->SetPosition(PositionOftheStringCreation+aPosition);
    333 //G4cout<<"Hadron "<<C1<<" "<<Hadron->GetPosition()/fermi<<" "<<Hadron->GetFormationTime()/fermi<<G4endl;
    334 /* // Uzhi 11.07.09
    335 G4cout<<C1<<' '<<Hadron->GetDefinition()->GetParticleName()<<G4endl;
    336 G4cout<<Hadron->GetDefinition()->GetPDGMass()<<' '
    337 <<Hadron->Get4Momentum()<<G4endl;
    338 G4cout<<Hadron->GetFormationTime()<<' '
    339 <<Hadron->GetPosition()<<' '
    340 <<Hadron->GetPosition().z()/fermi<<G4endl;
    341 */  // Uzhi
    342298        };
    343299
    344 //G4cout<<"Out FragmentString"<<G4endl;
    345300        return LeftVector;
    346301               
    … …  
    350305G4bool G4LundStringFragmentation::IsFragmentable(const G4FragmentingString * const string)
    351306{
    352 //G4cout<<"In IsFragmentable"<<G4endl;
    353   SetMinimalStringMass(string);                                                     // Uzhi
    354 //G4cout<<"Out IsFragmentable MinMass"<<MinimalStringMass<<" String Mass"<<std::sqrt(string->Get4Momentum().mag2())<<G4endl;
    355   return sqr(MinimalStringMass + WminLUND) < string->Get4Momentum().mag2();         // Uzhi
    356 
    357 //      return sqr(FragmentationMass(string)+MassCut) <                             // Uzhi
    358 //                      string->Mass2();                                            // Uzhi
     307  SetMinimalStringMass(string);                                                     
     308  return sqr(MinimalStringMass + WminLUND) < string->Get4Momentum().mag2();         
    359309}
    360310
    … …  
    362312G4bool G4LundStringFragmentation::StopFragmenting(const G4FragmentingString * const string)
    363313{
    364 //G4cout<<"StopFragmenting"<<G4endl;
    365 
    366314  SetMinimalStringMass(string);                                           
    367315
    368316//G4cout<<"StopFragm MinMass "<<MinimalStringMass<<" String Mass "<<std::sqrt(string->Get4Momentum().mag2())<<G4endl;
    369 
     317//G4cout<<"WminLUND "<<WminLUND<<" SmoothParam "<<SmoothParam<<" "<<string->Mass()<<G4endl;
    370318  return (MinimalStringMass + WminLUND)*
    371319             (1 + SmoothParam * (1.-2*G4UniformRand())) >               
    … …  
    379327{
    380328    //... perform last cluster decay
    381 /*
    382 G4cout<<"SplitLast String mass "<<string->Mass()<<G4endl;
    383 G4cout<<string->GetLeftParton()->GetPDGEncoding()<<" "<<G4endl;
    384 G4cout<<string->GetRightParton()->GetPDGEncoding()<<" "<<G4endl;
    385 */
     329
    386330    G4LorentzVector Str4Mom=string->Get4Momentum();
    387 //G4cout<<"String 4 momentum "<<Str4Mom<<G4endl;
     331
    388332    G4ThreeVector ClusterVel =string->Get4Momentum().boostVector();
    389333
    … …  
    395339    do
    396340    {
    397 //G4cout<<" Cicle "<<cClusterInterrupt<<" "<< ClusterLoopInterrupt<<G4endl;
    398 
    399341        if (cClusterInterrupt++ >= ClusterLoopInterrupt)
    400342        {
    … …  
    454396           }
    455397        }
    456 /*
    457 G4cout<<"SplitLast Left Right hadrons "<<LeftHadron->GetPDGEncoding()<<" "<<RightHadron->GetPDGEncoding()<<G4endl;
    458 G4cout<<"SplitLast Left Right hadrons "<<LeftHadron->GetPDGMass()<<" "<<RightHadron->GetPDGMass()<<G4endl;
    459 G4cout<<"Sum H mass Str Mass "<<LeftHadron->GetPDGMass() + RightHadron->GetPDGMass()<<" "<<ResidualMass<<G4endl;
    460 */
     398
    461399       //... repeat procedure, if mass of cluster is too low to produce hadrons
    462400       //... ClusterMassCut = 0.15*GeV model parameter
    463 
    464401    }
    465     while (ResidualMass <= LeftHadron->GetPDGMass() + RightHadron->GetPDGMass());// UzhiVOVA
     402    while (ResidualMass <= LeftHadron->GetPDGMass() + RightHadron->GetPDGMass());
    466403
    467404    //... compute hadron momenta and energies   
    … …  
    480417    RightVector->push_back(new G4KineticTrack(RightHadron, 0, Pos, RightMom));
    481418
    482 //G4cout<<"Out SplitLast "<<G4endl;
    483419    return true;
    484420
    … …  
    493429     G4double MassMt2, AntiMassMt2;                                         
    494430     G4double AvailablePz, AvailablePz2;                                   
    495 
    496 //G4cout<<"Sample4Momentum "<<G4endl;
    497 //G4cout<<"Sample4Momentum Mass"<<Mass<<" "<<AntiMass<<" "<<InitialMass<<G4endl;
    498431                                                                           
    499432    if(Mass > 930. || AntiMass > 930.)              // If there is a baryon
    … …  
    528461         G4double termN = 2*termD + 4*Mass*Mass + 4*AntiMass*AntiMass;
    529462         G4double pt2max=(termD*termD - termab )/ termN ;
    530          
    531 //       G4cout << " termD, ab, N " << termD << "  " << termab << "  " << termN
    532 //              <<  "   pt2max= " << pt2max ;
    533                                                                                          
     463                                               
    534464         Pt=SampleQuarkPt(std::sqrt(pt2max)); Pt.setZ(0); G4double Pt2=Pt.mag2();
    535 
    536          
    537 //       G4cout << " sampled Pt2 = " << Pt2 << "  " << pt2max-Pt2 << G4endl;             
    538          // end.. GF
    539 
    540 //G4cout<<"Sample4Momentum Pt x y "<<Pt.getX()<<" "<<Pt.getY()<<G4endl;
    541465
    542466         MassMt2    =     Mass *     Mass + Pt2;                             
    543467         AntiMassMt2= AntiMass * AntiMass + Pt2;                             
    544 
    545 //G4cout<<"Mts "<<MassMt2<<" "<<AntiMassMt2<<" "<<InitialMass*InitialMass<<G4endl;
    546468
    547469         AvailablePz2= sqr(InitialMass*InitialMass - MassMt2 - AntiMassMt2) -
    … …  
    552474      AvailablePz2 /=(4.*InitialMass*InitialMass);                           
    553475      AvailablePz = std::sqrt(AvailablePz2);                               
    554 
    555 //G4cout<<"AvailablePz "<<AvailablePz<<G4endl;
    556476 
    557477      G4double Px=Pt.getX();                                                 
    558478      G4double Py=Pt.getY();                                           
    559479
    560 //if(Mass > AntiMass){AvailablePz=-AvailablePz;}  // May30                   // Uzhi
    561 
    562480      Mom->setPx(Px); Mom->setPy(Py); Mom->setPz(AvailablePz);             
    563481      Mom->setE(std::sqrt(MassMt2+AvailablePz2));                         
    564482
    565 //G4cout<<" 1 part "<<Px<<"  "<<Py<<" "<<AvailablePz<<" "<<std::sqrt(MassMt2+AvailablePz2)<<G4endl;
    566 
    567483     AntiMom->setPx(-Px); AntiMom->setPy(-Py); AntiMom->setPz(-AvailablePz);
    568484     AntiMom->setE (std::sqrt(AntiMassMt2+AvailablePz2));                   
    569 
    570 //G4cout<<" 2 part "<<-Px<<"  "<<-Py<<" "<<-AvailablePz<<" "<<std::sqrt(AntiMassMt2+AvailablePz2)<<G4endl;
    571485    }
    572 
    573 //G4cout<<"Out Sample4Momentum "<<G4endl;
    574486  }
    575487
    … …  
    579491        G4FragmentingString * string, G4FragmentingString * newString)
    580492{
    581 /*     
    582 G4cout<<"SplitEandP "<<G4endl;
    583 G4cout<<"SplitEandP string mass "<<string->Mass()<<G4endl;   
    584 G4cout<<string->GetLeftParton()->GetPDGEncoding()<<" "
    585       <<string->GetRightParton()->GetPDGEncoding()<<" "<<G4endl;
    586 G4cout<<G4endl;
    587 G4cout<<newString->GetLeftParton()->GetPDGEncoding()<<" "<<G4endl;
    588 G4cout<<newString->GetRightParton()->GetPDGEncoding()<<" "<<G4endl;
    589 */
    590493       G4LorentzVector String4Momentum=string->Get4Momentum();
    591494       G4double StringMT2=string->Get4Momentum().mt2();
    592 //G4cout<<"StringMt2 "<<StringMT2<<G4endl;
    593495
    594496       G4double HadronMass = pHadron->GetPDGMass();
    595 //G4cout<<"Hadron mass "<<HadronMass<<" "<<pHadron->GetParticleName()<<G4endl;   
    596497
    597498       SetMinimalStringMass(newString);                                       
    598499       String4Momentum.setPz(0.);
    599500       G4ThreeVector StringPt=String4Momentum.vect();
    600 //G4cout<<"StringPt "<<StringPt<<G4endl<<G4endl;
    601 //G4cout<<"Min string mass "<<MinimalStringMass<<G4endl;
    602501
    603502// calculate and assign hadron transverse momentum component HadronPx and HadronPy
    … …  
    607506       G4ThreeVector HadronPt = thePt +string->DecayPt();
    608507       HadronPt.setZ(0);
    609 //G4cout<<"Hadron Pt"<<HadronPt<<G4endl;
    610508
    611509       G4ThreeVector RemSysPt = StringPt - HadronPt;
    612 //G4cout<<"RemSys Pt"<<RemSysPt<<G4endl;
    613510
    614511       //...  sample z to define hadron longitudinal momentum and energy
    … …  
    618515       G4double ResidualMassT2=sqr(MinimalStringMass) + RemSysPt.mag2();
    619516
    620 //G4cout<<"Mt h res str "<<std::sqrt(HadronMassT2)<<" "<<std::sqrt(ResidualMassT2)<<" srt mass"<<StringMT2<<G4endl;
    621 
    622517        G4double Pz2 = (sqr(StringMT2 - HadronMassT2 - ResidualMassT2) -             
    623518                                      4*HadronMassT2 * ResidualMassT2)/4./StringMT2;
    624 
    625 //G4cout<<"Pz**2 "<<Pz2<<G4endl;
    626519
    627520        if(Pz2 < 0 ) {return 0;}          // have to start all over!           
    … …  
    632525       G4double zMin = (std::sqrt(HadronMassT2+Pz2) - Pz)/std::sqrt(StringMT2); 
    633526       G4double zMax = (std::sqrt(HadronMassT2+Pz2) + Pz)/std::sqrt(StringMT2); 
    634 
    635 //G4cout<<"Zmin max "<<zMin<<" "<<zMax<<G4endl;               // Uzhi
    636527
    637528       if (zMin >= zMax) return 0;              // have to start all over!
    … …  
    652543
    653544       G4LorentzVector * a4Momentum= new G4LorentzVector(HadronPt,HadronE);
    654 
    655 //G4cout<<"Hadron Pt"<<HadronPt<<G4endl;
    656 //G4cout<<"Out of SplitEandP Pz E "<<HadronPt.getZ()<<" "<<HadronE<<G4endl;
    657545
    658546       return a4Momentum;
    … …  
    694582    {                                                           
    695583     // ---------------- Di-quark fragmentation ----------------------
    696 //G4cout<<"Di-quark"<<G4endl; // Vova
    697584       alund=0.7/GeV/GeV;    // 0.7 2.0
    698585       G4double zOfMaxyf=alund*Mt2/(alund*Mt2 + 1.);
    … …  
    707594    };                                                           
    708595
    709 //G4cout<<" test z "<<std::pow(2.,3.)<<" "<<z<<G4endl;
    710596    return z;
    711597    }
Note: See TracChangeset for help on using the changeset viewer.