G4StringChipsInterface.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // !!! Was used in QBBC PL, NOW it is not. Must be absolete !!!
00027 // =----------------------------------------------------------=
00028 
00029 //#define debug
00030 //#define pdebug
00031 
00032 #include <utility>
00033 #include <list>
00034 
00035 #include "G4StringChipsInterface.hh"
00036 #include "globals.hh"
00037 #include "G4SystemOfUnits.hh"
00038 #include "G4KineticTrackVector.hh"
00039 #include "G4Nucleon.hh"
00040 #include "G4LorentzRotation.hh"
00041 
00042 G4StringChipsInterface::G4StringChipsInterface()
00043 {
00044 #ifdef CHIPSdebug
00045   G4cout << "Please enter the energy loss per fermi in GeV"<<G4endl;
00046   G4cin >> theEnergyLossPerFermi;
00047 #endif
00048 #ifdef debug
00049   G4cout<<"G4StringChipsInterface::Constructor is called"<<G4endl;
00050 #endif
00051   theEnergyLossPerFermi = 0.5*GeV;
00052   // theEnergyLossPerFermi = 1.*GeV;
00053 }
00054 
00055 G4HadFinalState* G4StringChipsInterface::
00056 ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& theNucleus)
00057 {
00058 #ifdef debug
00059   G4cout<<"G4StringChipsInterface::ApplyYourself is called"<<G4endl;
00060 #endif
00061   return theModel.ApplyYourself(aTrack, theNucleus);
00062 }
00063 
00064 G4ReactionProductVector* G4StringChipsInterface::
00065 Propagate(G4KineticTrackVector* theSecondaries, G4V3DNucleus* theNucleus)
00066 {
00067 #ifdef debug
00068   G4cout<<"G4StringChipsInterface::Propagate is called"<<G4endl;
00069 #endif
00070   // Protection for non physical conditions
00071   
00072   if(theSecondaries->size() == 1) 
00073     throw G4HadronicException(__FILE__, __LINE__, "G4StringChipsInterface: Only one particle from String models!");
00074   
00075   // Calculate the mean energy lost
00076   std::pair<G4double, G4double> theImpact = theNucleus->RefetchImpactXandY();
00077   G4double impactX = theImpact.first;
00078   G4double impactY = theImpact.second;
00079   G4double inpactPar2 = impactX*impactX + impactY*impactY;
00080   
00081   G4double radius2 = theNucleus->GetNuclearRadius(5*perCent);
00082   radius2 *= radius2;
00083   G4double pathlength = 0;
00084   if(radius2 - inpactPar2>0) pathlength = 2.*std::sqrt(radius2 - inpactPar2);
00085   G4double theEnergyLostInFragmentation = theEnergyLossPerFermi*pathlength/fermi;
00086   
00087   // now select all particles in range
00088   std::list<std::pair<G4double, G4KineticTrack *> > theSorted;
00089   std::list<std::pair<G4double, G4KineticTrack *> >::iterator current;
00090   for(unsigned int secondary = 0; secondary<theSecondaries->size(); secondary++)
00091   {
00092     G4LorentzVector a4Mom = theSecondaries->operator[](secondary)->Get4Momentum();
00093 
00094 #ifdef CHIPSdebug
00095     G4cout <<"ALL STRING particles "<<theSecondaries->operator[](secondary)->GetDefinition()->GetPDGCharge()<<" "
00096            << theSecondaries->operator[](secondary)->GetDefinition()->GetPDGEncoding()<<" "
00097            << a4Mom <<G4endl; 
00098 #endif
00099 
00100     G4double toSort = a4Mom.rapidity();
00101     std::pair<G4double, G4KineticTrack *> it;
00102     it.first = toSort;
00103     it.second = theSecondaries->operator[](secondary);
00104     G4bool inserted = false;
00105     for(current = theSorted.begin(); current!=theSorted.end(); current++)
00106     {
00107       if((*current).first > toSort)
00108       {
00109                theSorted.insert(current, it);
00110                inserted = true;
00111                break;
00112       }
00113     }
00114     if(!inserted)
00115     {
00116       theSorted.push_front(it);
00117     }
00118   }
00119   
00120   G4LorentzVector proj4Mom(0.,0.,0.,0.);
00121   // @@ Use the G4QContent class, which is exactly (nD,nU,nS,nAD,nAU,nAS) !
00122   // The G4QContent class is a basic clas in CHIPS (not PDG Code as in GEANT4),
00123   // so in CHIPS on can has a hadronic obgect (Quasmon) with any Quark Content.
00124   // As a simple extantion for the hadron (which is a special case for Hadron)
00125   // there is a clas G4QChipolino, which is a Quasmon, which can decay in two
00126   // hadrons. In future the three-hadron class can be added...
00127   G4int nD  = 0;
00128   G4int nU  = 0;
00129   G4int nS  = 0;
00130   G4int nAD = 0;
00131   G4int nAU = 0;
00132   G4int nAS = 0;
00133   std::list<std::pair<G4double, G4KineticTrack *> >::iterator firstEscaping = theSorted.begin();
00134   G4double runningEnergy = 0;
00135   G4int particleCount = 0;
00136   G4LorentzVector theLow = (*(theSorted.begin())).second->Get4Momentum();
00137   G4LorentzVector theHigh;
00138 
00139 #ifdef CHIPSdebug
00140   G4cout << "CHIPS ENERGY LOST "<<theEnergyLostInFragmentation<<G4endl;
00141   G4cout << "sorted rapidities event start"<<G4endl;
00142 #endif 
00143 
00144   G4ReactionProductVector * theResult = new G4ReactionProductVector;
00145   G4ReactionProduct * theSec;
00146   G4KineticTrackVector * secondaries;
00147 #ifdef pdebug
00148   G4cout<<"G4StringChipsInterface::Propagate: Absorption"<<G4endl; 
00149 #endif
00150   
00151   // first decay and add all escaping particles.
00152   for(current = theSorted.begin(); current!=theSorted.end(); current++)
00153   {
00154     firstEscaping = current;
00155     if((*current).second->GetDefinition()->GetQuarkContent(3)!=0 ||
00156        (*current).second->GetDefinition()->GetAntiQuarkContent(3) !=0)
00157     {
00158       G4KineticTrack * aResult = (*current).second;
00159       G4ParticleDefinition * pdef=aResult->GetDefinition();
00160       secondaries = NULL;
00161       if ( pdef->GetPDGWidth() > 0 && pdef->GetPDGLifeTime() < 5E-17*s )
00162       {
00163         secondaries = aResult->Decay();
00164       }
00165       if ( secondaries == NULL )
00166       {
00167         theSec = new G4ReactionProduct(aResult->GetDefinition());
00168         G4LorentzVector current4Mom = aResult->Get4Momentum();
00169         theSec->SetTotalEnergy(current4Mom.t());
00170         theSec->SetMomentum(current4Mom.vect());
00171         theResult->push_back(theSec);
00172       } 
00173       else
00174       {
00175         for (unsigned int aSecondary=0; aSecondary<secondaries->size(); aSecondary++)
00176         {
00177           theSec = new G4ReactionProduct(secondaries->operator[](aSecondary)->GetDefinition());
00178           G4LorentzVector current4Mom = secondaries->operator[](aSecondary)->Get4Momentum();
00179           theSec->SetTotalEnergy(current4Mom.t());
00180           theSec->SetMomentum(current4Mom.vect());
00181           theResult->push_back(theSec);
00182         }
00183         std::for_each(secondaries->begin(), secondaries->end(), DeleteKineticTrack());
00184         delete secondaries;
00185       }
00186       continue;
00187     }
00188     runningEnergy += (*current).second->Get4Momentum().t();
00189     if(runningEnergy > theEnergyLostInFragmentation) break;
00190     
00191 #ifdef CHIPSdebug
00192     G4cout <<"ABSORBED STRING particles "<<current->second->GetDefinition()->GetPDGCharge()<<" "
00193            << current->second->GetDefinition()->GetPDGEncoding()<<" "
00194                   << current->second->Get4Momentum() <<G4endl; 
00195 #endif
00196 #ifdef pdebug
00197     G4cout<<"G4StringChipsInterface::Propagate:C="
00198           <<current->second->GetDefinition()->GetPDGCharge()<<", PDG="
00199           << current->second->GetDefinition()->GetPDGEncoding()<<", 4M="
00200                  << current->second->Get4Momentum() <<G4endl; 
00201 #endif
00202 
00203     // projectile 4-momentum needed in constructor of quasmon
00204     particleCount++;
00205     theHigh = (*current).second->Get4Momentum(); 
00206     proj4Mom += (*current).second->Get4Momentum(); 
00207 
00208 #ifdef CHIPSdebug
00209     G4cout << "sorted rapidities "<<current->second->Get4Momentum().rapidity()<<G4endl;
00210 #endif
00211     
00212      // projectile quark contents needed for G4QContent construction (@@ ? -> G4QContent class)
00213     nD  += (*current).second->GetDefinition()->GetQuarkContent(1);
00214     nU  += (*current).second->GetDefinition()->GetQuarkContent(2);
00215     nS  += (*current).second->GetDefinition()->GetQuarkContent(3);
00216     nAD += (*current).second->GetDefinition()->GetAntiQuarkContent(1);
00217     nAU += (*current).second->GetDefinition()->GetAntiQuarkContent(2);
00218     nAS += (*current).second->GetDefinition()->GetAntiQuarkContent(3);
00219   }
00220   // construct G4QContent
00221 
00222 #ifdef CHIPSdebug
00223   G4cout << "Quark content: d="<<nD<<", u="<<nU<< ", s="<< nS << G4endl;
00224   G4cout << "Anti-quark content: anit-d="<<nAD<<", anti-u="<<nAU<< ", anti-s="<< nAS << G4endl;
00225 #endif
00226 
00227   G4QContent theProjectiles(nD, nU, nS, nAD, nAU, nAS);
00228 
00229 #ifdef CHIPSdebug
00230   G4cout << "G4QContent is constructed"<<endl;
00231 #endif
00232   
00233   // target properties needed in constructor of quasmon
00234   // remove all hit nucleons to get Target code
00235   theNucleus->StartLoop();
00236   G4Nucleon * aNucleon;
00237   G4int resA = 0;
00238   G4int resZ = 0;
00239   G4ThreeVector hitMomentum(0,0,0);
00240   G4double hitMass = 0;
00241   G4int hitCount = 0;
00242   while((aNucleon = theNucleus->GetNextNucleon()))
00243   {
00244     if(!aNucleon->AreYouHit())
00245     {
00246       resA++;
00247       resZ+=G4int (aNucleon->GetDefinition()->GetPDGCharge());
00248     }
00249     else
00250     {
00251       hitMomentum += aNucleon->GetMomentum().vect();
00252       hitMass += aNucleon->GetMomentum().m();
00253       hitCount ++;
00254     }
00255   }
00256   G4int targetPDGCode = 90000000 + 1000*resZ + (resA-resZ);
00257   G4double targetMass = theNucleus->GetMass();
00258   targetMass -= hitMass;
00259   G4double targetEnergy = std::sqrt(hitMomentum.mag2()+targetMass*targetMass);
00260   // !! @@ Target should be at rest: hitMomentum=(0,0,0) @@ !! M.K.
00261   G4LorentzVector targ4Mom(-1.*hitMomentum, targetEnergy);
00262   
00263   // construct the quasmon
00264   G4int nop = 85; // clusters up to Alpha cluster (Reduced)
00265   //  G4int nop = 122; // clusters up to Alpha cluster
00266   //  G4int nop = 152; // not reduced upto Li6
00267   G4double fractionOfSingleQuasiFreeNucleons = 0.5; // It is A-dependent (C=.85, U=.40)
00268   G4double fractionOfPairedQuasiFreeNucleons = 0.05;
00269   G4double clusteringCoefficient = 5.;
00270   G4double temperature = 180.;
00271   G4double halfTheStrangenessOfSee = 0.3; // = s/d = s/u
00272   G4double etaToEtaPrime = 0.3;
00273 
00274   G4QNucleus::SetParameters(fractionOfSingleQuasiFreeNucleons,
00275                             fractionOfPairedQuasiFreeNucleons,
00276                             clusteringCoefficient);
00277   G4Quasmon::SetParameters(temperature,
00278                            halfTheStrangenessOfSee,
00279                            etaToEtaPrime);
00280 
00281 #ifdef CHIPSdebug
00282   G4cout << "G4QNucleus parameters "<< fractionOfSingleQuasiFreeNucleons << " "
00283          << fractionOfPairedQuasiFreeNucleons << " "<< clusteringCoefficient << G4endl;
00284   G4cout << "G4Quasmon parameters "<< temperature << " "<< halfTheStrangenessOfSee << " "
00285          <<etaToEtaPrime << G4endl;
00286   G4cout << "The Target PDG code = "<<targetPDGCode<<G4endl;
00287   G4cout << "The projectile momentum = "<<1./MeV*proj4Mom<<G4endl;
00288   G4cout << "The target momentum = "<<1./MeV*targ4Mom<<G4endl;
00289 #endif
00290 
00291   
00292   // Chips expects all in target rest frame, along z.
00293   // G4QCHIPSWorld aWorld(nop);              // Create CHIPS World of nop particles
00294   G4QCHIPSWorld::Get()->GetParticles(nop);
00295   G4QHadronVector projHV;
00296   // target rest frame
00297   proj4Mom.boost(-1.*targ4Mom.boostVector());
00298   // now go along z
00299   G4LorentzRotation toZ;
00300   toZ.rotateZ(-1*proj4Mom.phi());
00301   toZ.rotateY(-1*proj4Mom.theta());
00302   proj4Mom = toZ*proj4Mom;
00303   G4LorentzRotation toLab(toZ.inverse());
00304 
00305 #ifdef CHIPSdebug
00306   G4cout << "a Boosted projectile vector along z"<<proj4Mom<<" "<<proj4Mom.mag()<<G4endl;
00307 #endif
00308   
00309   G4QHadron* iH = new G4QHadron(theProjectiles, 1./MeV*proj4Mom);
00310   projHV.push_back(iH);
00311 
00312   // now call chips with this info in place
00313   G4QHadronVector * output = 0;
00314   if (particleCount!=0)
00315   {
00316     G4QEnvironment* pan= new G4QEnvironment(projHV, targetPDGCode);
00317     try
00318     {
00319       output = pan->Fragment();
00320     }
00321     catch(G4HadronicException & aR)
00322     {
00323       G4cerr << "Exception thrown passing through G4ChiralInvariantPhaseSpace "<<G4endl;
00324       G4cerr << " targetPDGCode = "<< targetPDGCode <<G4endl;
00325       G4cerr << " Dumping the information in the pojectile list"<<G4endl;
00326       for(size_t i=0; i< projHV.size(); i++)
00327       {
00328                G4cerr <<"  Incoming 4-momentum and PDG code of "<<i<<"'th hadron: "
00329                <<" "<< projHV[i]->Get4Momentum()<<" "<<projHV[i]->GetPDGCode()<<G4endl;
00330       }
00331       throw;
00332     }
00333     std::for_each(projHV.begin(), projHV.end(), DeleteQHadron());
00334     projHV.clear();
00335     delete pan;
00336   }
00337   else output = new G4QHadronVector;
00338    
00339   // Fill the result.
00340 #ifdef CHIPSdebug
00341   G4cout << "NEXT EVENT"<<endl;
00342 #endif
00343   
00344   // first decay and add all escaping particles.
00345   for(current = firstEscaping; current!=theSorted.end(); current++)
00346   {
00347     G4KineticTrack * aResult = (*current).second;
00348     G4ParticleDefinition * pdef=aResult->GetDefinition();
00349     secondaries = NULL;
00350     if ( pdef->GetPDGWidth() > 0 && pdef->GetPDGLifeTime() < 5E-17*s )
00351     {
00352       secondaries = aResult->Decay();
00353     }
00354     if ( secondaries == NULL )
00355     {
00356       theSec = new G4ReactionProduct(aResult->GetDefinition());
00357       G4LorentzVector current4Mom = aResult->Get4Momentum();
00358       current4Mom = toLab*current4Mom;
00359       current4Mom.boost(targ4Mom.boostVector());
00360       theSec->SetTotalEnergy(current4Mom.t());
00361       theSec->SetMomentum(current4Mom.vect());
00362       theResult->push_back(theSec);
00363     } 
00364     else
00365     {
00366       for (unsigned int aSecondary=0; aSecondary<secondaries->size(); aSecondary++)
00367       {
00368         theSec = new G4ReactionProduct(secondaries->operator[](aSecondary)->GetDefinition());
00369         G4LorentzVector current4Mom = secondaries->operator[](aSecondary)->Get4Momentum();
00370         current4Mom = toLab*current4Mom;
00371         current4Mom.boost(targ4Mom.boostVector());
00372         theSec->SetTotalEnergy(current4Mom.t());
00373         theSec->SetMomentum(current4Mom.vect());
00374         theResult->push_back(theSec);
00375       }
00376       std::for_each(secondaries->begin(), secondaries->end(), DeleteKineticTrack());
00377       delete secondaries;
00378     }
00379   }
00380   std::for_each(theSecondaries->begin(), theSecondaries->end(), DeleteKineticTrack());
00381   delete theSecondaries;
00382     
00383   // now add the quasmon output
00384 #ifdef CHIPSdebug
00385   G4cout << "Number of particles from string"<<theResult->size()<<G4endl;
00386   G4cout << "Number of particles from chips"<<output->size()<<G4endl;
00387 #endif
00388 
00389   for(unsigned int particle = 0; particle < output->size(); particle++)
00390   {
00391     if(output->operator[](particle)->GetNFragments() != 0) 
00392     {
00393       delete output->operator[](particle);
00394       continue;
00395     }
00396     //  theSec = new G4ReactionProduct;   // JA - not used, and memory leaked (Coverity)
00397     G4int pdgCode = output->operator[](particle)->GetPDGCode();
00398     G4ParticleDefinition * theDefinition;
00399     // Note that I still have to take care of strange nuclei
00400     // For this I need the mass calculation, and a changed interface
00401     // for ion-table ==> work for Hisaya @@@@@@@
00402     // Then I can sort out the pdgCode. I also need a decau process 
00403     // for strange nuclei; may be another chips interface
00404     if(pdgCode>90000000) 
00405     {
00406       G4int aZ = (pdgCode-90000000)/1000;
00407       if (aZ>1000) aZ=aZ%1000;  // patch for strange nuclei, to be repaired @@@@
00408       G4int anN = pdgCode-90000000-1000*aZ;
00409       if(anN>1000) anN=anN%1000; // patch for strange nuclei, to be repaired @@@@
00410       if(pdgCode==91000000) theDefinition = G4Lambda::LambdaDefinition();
00411       else if(pdgCode==92000000) theDefinition = G4Lambda::LambdaDefinition();
00412       else if(pdgCode==93000000) theDefinition = G4Lambda::LambdaDefinition();
00413       else if(pdgCode==94000000) theDefinition = G4Lambda::LambdaDefinition();
00414       else if(pdgCode==95000000) theDefinition = G4Lambda::LambdaDefinition();
00415       else if(pdgCode==96000000) theDefinition = G4Lambda::LambdaDefinition();
00416       else if(pdgCode==97000000) theDefinition = G4Lambda::LambdaDefinition();
00417       else if(pdgCode==98000000) theDefinition = G4Lambda::LambdaDefinition();
00418       else if(aZ == 0 && anN == 1) theDefinition = G4Neutron::Neutron();
00419       else theDefinition = G4ParticleTable::GetParticleTable()->FindIon(aZ,anN+aZ,0,aZ);
00420     }    
00421     else
00422     {
00423       theDefinition = G4ParticleTable::GetParticleTable()->FindParticle(output->operator[](particle)->GetPDGCode());
00424     }
00425 #ifdef CHIPSdebug
00426     G4cout << "Particle code produced = "<< pdgCode <<G4endl;
00427 #endif
00428 
00429     theSec = new G4ReactionProduct(theDefinition);
00430     G4LorentzVector current4Mom = output->operator[](particle)->Get4Momentum();
00431     current4Mom = toLab*current4Mom;
00432     current4Mom.boost(targ4Mom.boostVector());
00433     theSec->SetTotalEnergy(current4Mom.t());
00434     theSec->SetMomentum(current4Mom.vect());
00435     theResult->push_back(theSec);
00436     
00437 #ifdef CHIPSdebug
00438     G4cout <<"CHIPS particles "<<theDefinition->GetPDGCharge()<<" "
00439            << theDefinition->GetPDGEncoding()<<" "
00440            << current4Mom <<G4endl; 
00441 #endif
00442 
00443     delete output->operator[](particle);
00444   }
00445   delete output;
00446 #ifdef CHIPSdebug
00447   G4cout << "Number of particles"<<theResult->size()<<G4endl;
00448   G4cout << G4endl;
00449   // @@ G4QContent has even the out option!
00450   G4cout << "QUASMON preparation info "
00451          << 1./MeV*proj4Mom<<" "
00452          << 1./MeV*targ4Mom<<" "
00453          << nD<<" "<<nU<<" "<<nS<<" "<<nAD<<" "<<nAU<<" "<<nAS<<" "
00454          << hitCount<<" "
00455          << particleCount<<" "
00456          << theLow<<" "
00457          << theHigh<<" "
00458          << G4endl;
00459 #endif
00460 
00461   return theResult;
00462 } 

Generated on Mon May 27 17:49:55 2013 for Geant4 by  doxygen 1.4.7