Ignore:
Timestamp:
Jun 18, 2010, 11:42:07 AM (14 years ago)
Author:
garnier
Message:

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08

Location:
trunk/source/processes/hadronic/models/chiral_inv_phase_space/fragmentation/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/models/chiral_inv_phase_space/fragmentation/src/G4QFragmentation.cc

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4QFragmentation.cc,v 1.6 2009/12/16 17:51:03 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4QFragmentation.cc,v 1.14 2010/06/11 10:35:18 mkossov Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030// -----------------------------------------------------------------------------
    … …  
    5454// Definition of static parameters
    5555G4int    G4QFragmentation::nCutMax=7;
    56 G4double G4QFragmentation::stringTension=1.5*GeV/fermi;
    57 G4double G4QFragmentation::tubeDensity  =1.5/fermi;
     56G4double G4QFragmentation::stringTension=GeV/fermi;
     57//G4double G4QFragmentation::tubeDensity  =1./fermi;
     58G4double G4QFragmentation::tubeDensity  =0./fermi;
    5859// Parameters of diffractional fragmentation (was .72*)
    5960G4double G4QFragmentation::widthOfPtSquare=-GeV*GeV;// pt -width2 forStringExcitation
    … …  
    473474        aTarget=theInteractions[0]->GetTarget();
    474475        aProjectile=theInteractions[0]->GetProjectile();
     476      delete theInteractions[0];
    475477      theInteractions.clear();
    476       delete theInteractions[0];
    477478    }
    478479    else                                             // Create a new target nucleon
    … …  
    555556  G4cout<<"G4QFragmentation::Construct: Creation ofSoftCollisionPartonPair STARTS"<<G4endl;
    556557#endif
    557   for(it = theInteractions.begin(); it != theInteractions.end(); ++it)   
    558   {
     558  G4bool rep=true;
     559  while(rep && theInteractions.size())
     560  {
     561   for(it = theInteractions.begin(); it != theInteractions.end(); ++it)
     562   {
    559563    G4QInteraction* anIniteraction = *it;
    560564    G4QPartonPair*  aPair=0;
    … …  
    583587      delete *it;
    584588      it=theInteractions.erase(it);      // Soft interactions are converted & erased
    585       it--;
    586     }
     589      if( it != theInteractions.begin() )// To avoid going below begin() (for Windows)
     590      {
     591        it--;
     592        rep=false;
     593#ifdef debug
     594        G4cout<<"G4QFragmentation::Construct: *** Decremented ***"<<G4endl;
     595#endif
     596      }
     597      else
     598      {
     599        rep=true;
     600#ifdef debug
     601        G4cout<<"G4QFragmentation::Construct: *** IT Begin ***"<<G4endl;
     602#endif
     603        break;
     604      }
     605    }
     606    else rep=false;
     607#ifdef debug
     608    G4cout<<"G4QFragmentation::Construct: #0fSC="<<nSoftCollisions<<", r="<<rep<<G4endl;
     609#endif
     610   }
     611#ifdef debug
     612   G4cout<<"G4QFragmentation::Construct: *** IT While *** , r="<<rep<<G4endl;
     613#endif
    587614  }
    588615#ifdef debug
    … …  
    748775  G4int problem=0;                                   // 0="no problem", incremented by ASIS
    749776  G4QStringVector::iterator ist;
    750   for(ist = strings.begin(); ist < strings.end(); ist++)
    751   {
     777  G4bool con=true;
     778  while(con && strings.size())
     779  {
     780   for(ist = strings.begin(); ist < strings.end(); ++ist)
     781   {
    752782    G4bool bad=true;
    753783    G4LorentzVector cS4M=(*ist)->Get4Momentum();
    … …  
    15801610        delete (*ist);
    15811611        strings.erase(ist);
    1582         ist--;
    15831612#ifdef debug
    15841613        G4LorentzVector ss4M=pL4M+pR4M;
    15851614        G4cout<<"G4QFragmentation::Construct:Created,4M="<<ss4M<<",m2="<<ss4M.m2()<<G4endl;
    15861615#endif
     1616        if( ist != strings.begin() ) // To avoid going below begin() (for Windows)
     1617        {
     1618          ist--;
     1619          con=false;
     1620#ifdef debug
     1621          G4cout<<"G4QFragmentation::Construct: *** IST Decremented ***"<<G4endl;
     1622#endif
     1623        }
     1624        else
     1625        {
     1626          con=true;
     1627#ifdef debug
     1628          G4cout<<"G4QFragmentation::Construct: *** IST Begin ***"<<G4endl;
     1629#endif
     1630          break;
     1631        }
    15871632      } // End of the IF(the best partnerString candidate was found)
    15881633      else
    … …  
    15921637#endif
    15931638        ++problem;
     1639        con=false;
    15941640      }
    1595     }
    1596   }
     1641    } // End of IF
     1642    else con=false;
     1643   } // End of loop over ist iterator
     1644#ifdef debug
     1645   G4cout<<"G4QFragmentation::Construct: *** IST While *** , con="<<con<<G4endl;
     1646#endif
     1647  } // End of "con" while
    15971648#ifdef edebug
    15981649  // This print has meaning only if something appear between it and the StringFragmLOOP
    … …  
    18261877      G4int miPDG=qsumQC.GetSPDGCode();                     // PDG of minM of hadron/fragm.
    18271878      G4double gsM=0.;                                      // Proto minM of had/frag forQC
     1879#ifdef debug
     1880      G4cout<<"G4QFragmentation::Fragment: QC="<<qsumQC<<",PDG="<<miPDG<<G4endl;
     1881#endif
    18281882      if(miPDG == 10)
    18291883      {
    … …  
    18331887        //    theWorld->GetQParticle(QCh.GetQPDG2())->MinMassOfFragm();
    18341888      }
    1835       else if(miPDG>80000000)                               // Compound Nucleus
    1836       {
    1837         G4QNucleus rtN(qsumQC);                  // Create PseudoNucl for totCompound
    1838         gsM=rtN.GetGSMass(); // MinMass of residQ+(Env-ParC) syst.      }
    1839       }
    1840       else if(std::abs(miPDG)%10 > 2)
     1889      // @@ it is not clear, why it does not work ?!
     1890      //else if(miPDG>80000000)                             // Compound Nucleus
     1891      //{
     1892      //  G4QNucleus rtN(qsumQC);                           // CreatePseudoNucl for totComp
     1893      //  gsM=rtN.GetGSMass();                              // MinMass of residQ+(Env-ParC)
     1894      //}
     1895      else if(miPDG < 80000000 && std::abs(miPDG)%10 > 2)
    18411896                           gsM=theWorld->GetQParticle(G4QPDGCode(miPDG))->MinMassOfFragm();
    1842       else gsM=G4QPDGCode(miPDG).GetMass();      // minM of hadron/fragm. for QC
     1897      else gsM=G4QPDGCode(miPDG).GetMass();                 // minM of hadron/fragm. for QC
    18431898      G4double reM=qsum4M.m();                              // real mass of the compound
    18441899#ifdef debug
    … …  
    18611916      else
    18621917      {
    1863         G4cerr<<"***G4QFragmentation::Fragm:PDG="<<miPDG<<",M="<<reM<<",GSM="<<gsM<<G4endl;
     1918        G4cerr<<"*G4QFr::Fr:PDG="<<miPDG<<",M="<<reM<<",GSM="<<gsM<<",QC="<<qsumQC<<G4endl;
    18641919        G4Exception("G4QFragmentation::Fragment:","27",FatalException,"Can't recover GSM");
    18651920      }
    … …  
    22772332    G4QString* curString=strings[astring];
    22782333    if(!curString->GetDirection()) continue;  // Historic for the dead strings: DoesNotWork
    2279 #ifdef edebug
    22802334    G4int curStrChg = curString->GetCharge();
    22812335    G4int curStrBaN = curString->GetBaryonNumber();
    2282 #endif
    22832336    G4LorentzVector curString4M = curString->Get4Momentum();
    22842337#ifdef debug
    … …  
    28162869#endif
    28172870                G4LorentzVector p4M=selHP->Get4Momentum();
     2871                //
    28182872                curString4M+=p4M;
    2819 #ifdef edebug
    28202873                G4int Chg=selHP->GetCharge();
    28212874                G4int BaN=selHP->GetBaryonNumber();
    28222875                curStrChg+=Chg;
    28232876                curStrBaN+=BaN;
     2877#ifdef edebug
    28242878                G4cout<<"-EMC->>>>G4QFragmentation::Breeder: S+=H, 4M="<<curString4M<<",M="
    28252879                      <<curString4M.m()<<", Charge="<<curStrChg<<", B="<<curStrBaN<<G4endl;
    … …  
    35103564      G4QHadron* curHadron=(*theHadrons)[aTrack];
    35113565      G4int hPDG=curHadron->GetPDGCode();
    3512 #ifdef edebug
    35133566      G4LorentzVector curH4M=curHadron->Get4Momentum();
    35143567      G4int           curHCh=curHadron->GetCharge();
    35153568      G4int           curHBN=curHadron->GetBaryonNumber();
    3516 #endif
    35173569#ifdef debug
    35183570      G4cout<<">>>>>>>>G4QFragmentation::Breeder:S#"<<astring<<",H#"<<aTrack<<",PDG="<<hPDG
    … …  
    35313583        {
    35323584          theResult->push_back((*tmpQHadVec)[aH]);// TheDecayProduct of TheHadron is filled
    3533 #ifdef edebug
     3585          //
    35343586          G4QHadron*   prodH =(*tmpQHadVec)[aH];
    35353587          G4LorentzVector p4M=prodH->Get4Momentum();
    3536           G4int           PDG=prodH->GetPDGCode();
    35373588          G4int           Chg=prodH->GetCharge();
    35383589          G4int           BaN=prodH->GetBaryonNumber();
    3539           curH4M-=p4M;
    35403590          curString4M-=p4M;
    35413591          curStrChg-=Chg;
    35423592          curStrBaN-=BaN;
     3593          curH4M-=p4M;
    35433594          curHCh-=Chg;
    35443595          curHBN-=BaN;
     3596#ifdef edebug
     3597          G4int           PDG=prodH->GetPDGCode();
    35453598          G4cout<<"-EMC->>>>G4QFragmentation::Breeder:Str*Filled, 4M="<<p4M<<", PDG="<<PDG
    35463599                <<", Chg="<<Chg<<", BaN="<<BaN<<G4endl;
    … …  
    35563609      {
    35573610        theResult->push_back(curHadron);        // The original hadron is filled
    3558 #ifdef edebug
     3611        //
    35593612        curString4M-=curH4M;
    35603613        G4int curCh=curHadron->GetCharge();
    … …  
    35623615        curStrChg-=curCh;
    35633616        curStrBaN-=curBN;
     3617#ifdef edebug
    35643618        G4cout<<"-EMC->>>>>>G4QFragmentation::Breeder: curH filled 4M="<<curH4M<<",PDG="
    35653619              <<curHadron->GetPDGCode()<<", Chg="<<curCh<<", BaN="<<curBN<<G4endl;
    … …  
    35733627          <<", rChg="<<curStrChg<<", rBaN="<<curStrBaN<<G4endl;
    35743628#endif
     3629    // Trap with the debugging warning --- Starts ---
     3630    if(curStrChg || curStrBaN || curString4M.t() > eps || std::fabs(curString4M.x()) > eps
     3631       || std::fabs(curString4M.y()) > eps || std::fabs(curString4M.z()) > eps )
     3632    {
     3633      G4double dEn=curString4M.t();
     3634      G4double dPx=curString4M.x();
     3635      G4double dPy=curString4M.y();
     3636      G4double dPz=curString4M.z();
     3637      G4int nHadr=theResult->size();
     3638      G4double hEn=0.;
     3639      G4double hPx=0.;
     3640      G4double hPy=0.;
     3641      G4double hPz=0.;
     3642      G4int    hCh=0;
     3643      G4int    hBN=0;
     3644      G4double mEn=0.;
     3645      G4double mPx=0.;
     3646      G4double mPy=0.;
     3647      G4double mPz=0.;
     3648      G4int    mCh=0;
     3649      G4int    mBN=0;
     3650      for(G4int i=0; i<nHadr; i++)
     3651      {
     3652        mEn=hEn;
     3653        mPx=hPx;
     3654        mPy=hPy;
     3655        mPz=hPz;
     3656        mCh=hCh;
     3657        mBN=hBN;
     3658        G4QHadron* curHadr = (*theResult)[i];
     3659        G4LorentzVector hI4M = curHadr->Get4Momentum();
     3660        hEn=hI4M.t();
     3661        hPx=hI4M.x();
     3662        hPy=hI4M.y();
     3663        hPz=hI4M.z();
     3664        hCh=curHadr->GetCharge();
     3665        hBN=curHadr->GetBaryonNumber();
     3666        G4cout<<"G4QFragmentation::Breeder: H#"<<i<<", d4M="<<curString4M+hI4M
     3667              <<",dCh="<<hCh+curStrChg<<",dBN="<<hBN+curStrBaN<<G4endl;
     3668        if( !(hCh+curStrChg) && !(hBN+curStrBaN) && std::fabs(dEn+hEn)<eps &&
     3669            std::fabs(dPx+hPx)<eps && std::fabs(dPy+hPy)<eps && std::fabs(dPz+hPz)<eps )
     3670        {
     3671          G4cout<<"G4QFragmentation::Breeder: ***Cured*** Redundent Hadron # "<<i<<G4endl;
     3672          G4QHadron* theLast = (*theResult)[nHadr-1];
     3673          curHadr->Set4Momentum(theLast->Get4Momentum()); //4-Mom of CurHadr
     3674          G4QPDGCode lQP=theLast->GetQPDG();
     3675          if(lQP.GetPDGCode()!=10) curHadr->SetQPDG(lQP);
     3676          else curHadr->SetQC(theLast->GetQC());
     3677          theResult->pop_back(); // theLastQHadron is excluded from OUTPUT
     3678          delete theLast;        //*!!When kill, delete theLastQHadr as an Instance!*
     3679          break;
     3680        }
     3681        if( !(hCh+mCh+curStrChg) && !(hBN+mBN+curStrBaN) && std::fabs(dEn+hEn+mEn)<eps &&
     3682            std::fabs(dPx+hPx+mPx)<eps && std::fabs(dPy+hPy+mPy)<eps &&
     3683            std::fabs(dPz+hPz+mPz)<eps )
     3684        {
     3685          G4cout<<"G4QFragmentation::Breeder:***Cured*** Redundent 2Hadrons i="<<i<<G4endl;
     3686          G4QHadron* theLast = (*theResult)[nHadr-1];
     3687          curHadr->Set4Momentum(theLast->Get4Momentum()); //4-Mom of CurHadr
     3688          G4QPDGCode lQP=theLast->GetQPDG();
     3689          if(lQP.GetPDGCode()!=10) curHadr->SetQPDG(lQP);
     3690          else curHadr->SetQC(theLast->GetQC());
     3691          theResult->pop_back(); // theLastQHadron is excluded from OUTPUT
     3692          delete theLast;        //*!!When kill, delete theLastQHadr as an Instance!*
     3693          theLast = (*theResult)[nHadr-2];
     3694          curHadr->Set4Momentum(theLast->Get4Momentum()); //4-Mom of CurHadr
     3695          lQP=theLast->GetQPDG();
     3696          if(lQP.GetPDGCode()!=10) curHadr->SetQPDG(lQP);
     3697          else curHadr->SetQC(theLast->GetQC());
     3698          theResult->pop_back(); // theLastQHadron is excluded from OUTPUT
     3699          delete theLast;        //*!!When kill, delete theLastQHadr as an Instance!*
     3700          break;
     3701        }
     3702        // If the redundent particle decay in 3 FS hadrons -> the corresponding Improvement
     3703        G4cout<<"*Warning*G4QFragmentation::Breeder: Nonconservation isn't cured!"<<G4endl;
     3704      }
     3705    }
     3706    // Trap with the debugging warning ^^^  Ends  ^^^
    35753707  } // End of the main LOOP over decaying strings
    35763708  G4LorentzVector r4M=theNucleus.Get4Momentum(); // Nucleus 4-momentum in LS

  • trunk/source/processes/hadronic/models/chiral_inv_phase_space/fragmentation/src/G4QIonIonCollision.cc

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4QIonIonCollision.cc,v 1.2 2009/12/01 09:24:24 mkossov Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4QIonIonCollision.cc,v 1.8 2010/04/01 15:03:35 mkossov Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030// -----------------------------------------------------------------------------
    … …  
    5757G4double G4QIonIonCollision::tubeDensity  =1./fermi;
    5858// Parameters of diffractional fragmentation
    59 G4double G4QIonIonCollision::widthOfPtSquare=-0.72*GeV*GeV; // ptWidth2 forStringExcitation
     59G4double G4QIonIonCollision::widthOfPtSquare=-0.75*GeV*GeV; // ptWidth2 forStringExcitation
    6060
    6161G4QIonIonCollision::G4QIonIonCollision(G4QNucleus &pNucleus, const G4QNucleus &tNucleus)
    … …  
    386386        aTarget=theInteractions[0]->GetTarget();
    387387             aProjectile=theInteractions[0]->GetProjectile();
     388      delete theInteractions[0];
    388389      theInteractions.clear();
    389       delete theInteractions[0];
    390390    }
    391391    else                                             // Create a new target nucleon (?)
    … …  
    458458  G4cout<<"G4QIonIonCollision::Constr: Creation ofSoftCollisionPartonPair STARTS"<<G4endl;
    459459#endif
    460   for(it = theInteractions.begin(); it != theInteractions.end(); ++it)   
    461   {
     460  G4bool rep=true;
     461  while(rep && theInteractions.size())
     462  {
     463   for(it = theInteractions.begin(); it != theInteractions.end(); ++it)   
     464   {
    462465    G4QInteraction* anIniteraction = *it;
    463466    G4QPartonPair*  aPair=0;
    … …  
    475478                                  pProjectile->GetNextAntiParton(),
    476479                                  G4QPartonPair::SOFT, G4QPartonPair::TARGET);
    477         thePartonPairs.push_back(aPair);            // A target pair (Why TAGRET?)
     480        thePartonPairs.push_back(aPair); // A target pair (Why TAGRET?)
    478481        aPair = new G4QPartonPair(pProjectile->GetNextParton(),
    479482                                  pTarget->GetNextAntiParton(),
    480483                                  G4QPartonPair::SOFT, G4QPartonPair::PROJECTILE);
    481         thePartonPairs.push_back(aPair);            // A projectile pair (Why Projectile?)
     484        thePartonPairs.push_back(aPair); // A projectile pair (Why Projectile?)
    482485#ifdef debug
    483486        G4cout<<"--->G4QIonIonCollision::Constr: SOFT, 2 parton pairs are filled"<<G4endl;
    … …  
    485488      } 
    486489      delete *it;
    487       it=theInteractions.erase(it);                 // SoftInteractions're converted&erased
    488       it--;
    489     }
     490      it=theInteractions.erase(it);      // SoftInteractions're converted&erased
     491      if( it != theInteractions.begin() )// To avoid going below begin() (for Windows)
     492      {
     493        it--;
     494        rep=false;
     495      }
     496      else
     497      {
     498        rep=true;
     499        break;
     500      }
     501    }
     502    else rep=false;
     503   }
    490504  }
    491505#ifdef debug
    … …  
    633647  G4int problem=0;                                   // 0="no problem", incremented by ASIS
    634648  G4QStringVector::iterator ist;
    635   for(ist = strings.begin(); ist < strings.end(); ist++)
    636   {
     649  G4bool con=true;
     650  while(con && strings.size())
     651  {
     652   for(ist = strings.begin(); ist < strings.end(); ++ist)
     653   {
    637654    G4bool bad=true;
    638655    G4LorentzVector cS4M=(*ist)->Get4Momentum();
    … …  
    14001417        delete (*ist);
    14011418        strings.erase(ist);
    1402         ist--;
    14031419#ifdef debug
    14041420        G4LorentzVector ss4M=pL4M+pR4M;
    14051421        G4cout<<"G4QIonIonCollision::Constr:Created,4M="<<ss4M<<",m2="<<ss4M.m2()<<G4endl;
    14061422#endif
     1423        if( ist != strings.begin() ) // To avoid going below begin() (for Windows)
     1424        {
     1425          ist--;
     1426          con=false;
     1427#ifdef debug
     1428          G4cout<<"G4QIonIonCollision::Construct: *** IST Decremented ***"<<G4endl;
     1429#endif
     1430        }
     1431        else
     1432        {
     1433          con=true;
     1434#ifdef debug
     1435          G4cout<<"G4QIonIonCollision::Construct: *** IST Begin ***"<<G4endl;
     1436#endif
     1437          break;
     1438        }
    14071439      } // End of the IF(the best partnerString candidate was found)
    14081440      else
    … …  
    14121444#endif
    14131445        ++problem;
     1446        con=false;
    14141447      }
    14151448    }
    1416   }
     1449    else con=false;
     1450   } // End of loop over ist iterator
     1451#ifdef debug
     1452   G4cout<<"G4QIonIonCollision::Construct: *** IST While *** , con="<<con<<G4endl;
     1453#endif
     1454  } // End of "con" while
    14171455#ifdef edebug
    14181456  // This print has meaning only if something appear between it and the StringFragmLOOP
    … …  
    16451683        //    theWorld->GetQParticle(QCh.GetQPDG2())->MinMassOfFragm();
    16461684      }
    1647       else if(miPDG>80000000)                               // Compound Nucleus
    1648       {
    1649         G4QNucleus rtN(qsumQC);                  // Create PseudoNucl for totCompound
    1650         gsM=rtN.GetGSMass(); // MinMass of residQ+(Env-ParC) syst.      }
    1651       }
    1652       else if(std::abs(miPDG)%10 > 2)
     1685      // @@ it is not clear, why it does not work ?!
     1686      //else if(miPDG>80000000)                             // Compound Nucleus
     1687      //{
     1688      //  G4QNucleus rtN(qsumQC);                           // CreatePseudoNucl for totComp
     1689      //  gsM=rtN.GetGSMass();                              // MinMass of residQ+(Env-ParC)
     1690      //}
     1691      else if(miPDG < 80000000 && std::abs(miPDG)%10 > 2)
    16531692                           gsM=theWorld->GetQParticle(G4QPDGCode(miPDG))->MinMassOfFragm();
    1654       else gsM=G4QPDGCode(miPDG).GetMass();      // minM of hadron/fragm. for QC
     1693      else gsM=G4QPDGCode(miPDG).GetMass();                 // minM of hadron/fragm. for QC
    16551694      G4double reM=qsum4M.m();                              // real mass of the compound
    16561695#ifdef pdebug
Note: See TracChangeset for help on using the changeset viewer.