[968] | 1 | // |
---|
| 2 | // ******************************************************************** |
---|
| 3 | // * License and Disclaimer * |
---|
| 4 | // * * |
---|
| 5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
| 6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
| 7 | // * conditions of the Geant4 Software License, included in the file * |
---|
| 8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
| 9 | // * include a list of copyright holders. * |
---|
| 10 | // * * |
---|
| 11 | // * Neither the authors of this software system, nor their employing * |
---|
| 12 | // * institutes,nor the agencies providing financial support for this * |
---|
| 13 | // * work make any representation or warranty, express or implied, * |
---|
| 14 | // * regarding this software system or assume any liability for its * |
---|
| 15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
| 16 | // * for the full disclaimer and the limitation of liability. * |
---|
| 17 | // * * |
---|
| 18 | // * This code implementation is the result of the scientific and * |
---|
| 19 | // * technical work of the GEANT4 collaboration. * |
---|
| 20 | // * By using, copying, modifying or distributing the software (or * |
---|
| 21 | // * any work based on the software) you agree to acknowledge its * |
---|
| 22 | // * use in resulting scientific publications, and indicate your * |
---|
| 23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
| 24 | // ******************************************************************** |
---|
| 25 | // |
---|
| 26 | //#define debug |
---|
| 27 | //#define pdebug |
---|
| 28 | |
---|
| 29 | //#include "G4UIterminal.hh" |
---|
| 30 | #include "G4ios.hh" |
---|
| 31 | #include "G4QEnvironment.hh" |
---|
| 32 | |
---|
| 33 | int main() |
---|
| 34 | { |
---|
| 35 | //G4double theEnergyLossPerFermi = 1.*GeV; |
---|
| 36 | G4int nop = 152; // clusters (A<6) |
---|
| 37 | G4double fractionOfSingleQuasiFreeNucleons = 0.5; // It is A-dependent (C=.85, U=.40) |
---|
| 38 | G4double fractionOfPairedQuasiFreeNucleons = 0.05; |
---|
| 39 | G4double clusteringCoefficient = 5.; |
---|
| 40 | G4double temperature = 180.; |
---|
| 41 | G4double halfTheStrangenessOfSee = 0.3; // = s/d = s/u |
---|
| 42 | G4double etaToEtaPrime = 0.3; |
---|
| 43 | G4double fusionToExchange = 100.; |
---|
| 44 | //G4double theInnerCoreDensityCut = 50.; |
---|
| 45 | #ifdef debug |
---|
| 46 | G4cout<<"G4MulFragTst: nop="<<nop<<",fN="<<fractionOfSingleQuasiFreeNucleons |
---|
| 47 | <<",fD="<<fractionOfPairedQuasiFreeNucleons<<",wC="<<clusteringCoefficient |
---|
| 48 | <<",vn="<<fusionToExchange<<",T="<<temperature |
---|
| 49 | <<",ss="<<halfTheStrangenessOfSee<<",se="<<etaToEtaPrime<<G4endl; |
---|
| 50 | #endif |
---|
| 51 | //-------- Initialize CHIPS |
---|
| 52 | G4QCHIPSWorld* theW=G4QCHIPSWorld::Get(); |
---|
| 53 | theW->GetParticles(nop); // Create CHIPS World of nop particles |
---|
| 54 | G4QNucleus::SetParameters(fractionOfSingleQuasiFreeNucleons, |
---|
| 55 | fractionOfPairedQuasiFreeNucleons, |
---|
| 56 | clusteringCoefficient, |
---|
| 57 | fusionToExchange); |
---|
| 58 | G4Quasmon::SetParameters(temperature, halfTheStrangenessOfSee, etaToEtaPrime); |
---|
| 59 | |
---|
| 60 | G4int tPDG = 90006005; // PDG of the Target Nucleus |
---|
| 61 | G4int nEvents= 100000; |
---|
| 62 | //#ifdef debug |
---|
| 63 | G4cout<<"G4MulFragTst: targetPDG="<<tPDG<<", nEv="<<nEvents<<G4endl; |
---|
| 64 | //#endif |
---|
| 65 | G4QHadronVector projHV; |
---|
| 66 | |
---|
| 67 | G4int PDG1=211; |
---|
| 68 | G4LorentzVector LV1(-485.221,-216.252,267.887,611.103); |
---|
| 69 | //#ifdef debug |
---|
| 70 | G4cout<<"G4MulFragTst: #1 projPDG="<<PDG1<<", 4M="<<LV1<<G4endl; |
---|
| 71 | //#endif |
---|
| 72 | |
---|
| 73 | G4int PDG2=-211; |
---|
| 74 | G4LorentzVector LV2(-201.944,475.785,952.224,1092.41); |
---|
| 75 | //#ifdef debug |
---|
| 76 | G4cout<<"G4MulFragTst: #2 projPDG="<<PDG2<<", 4M="<<LV2<<G4endl; |
---|
| 77 | //#endif |
---|
| 78 | |
---|
| 79 | G4int PDG3=111; |
---|
| 80 | G4LorentzVector LV3(213.293,183.683,996.566,1044.32); |
---|
| 81 | //#ifdef debug |
---|
| 82 | G4cout<<"G4MulFragTst: #3 projPDG="<<PDG3<<", 4M="<<LV3<<G4endl; |
---|
| 83 | //#endif |
---|
| 84 | |
---|
| 85 | //std::ifstream inFile("chipstest.in", std::ios::in); |
---|
| 86 | //G4double temperature; |
---|
| 87 | //G4double ssin2g; |
---|
| 88 | //G4double eteps; |
---|
| 89 | //G4double momb; |
---|
| 90 | //G4double enb; |
---|
| 91 | //G4double cP; |
---|
| 92 | //G4double fN; |
---|
| 93 | //G4double fD; |
---|
| 94 | //G4double rM; |
---|
| 95 | //G4double sA; |
---|
| 96 | //G4int nop; |
---|
| 97 | //G4int pPDG; |
---|
| 98 | //G4int tPDG; |
---|
| 99 | //G4int nEvt; |
---|
| 100 | ////G4int nofdecays; |
---|
| 101 | ////G4int decmask=0; |
---|
| 102 | //inFile>>temperature>>ssin2g>>eteps>>nop>>momb>>enb>>pPDG>>tPDG>>nEvt>>fN>>fD>>cP>>rM>>sA; |
---|
| 103 | //G4cout<<"CHIPStest Par's: Temp="<<temperature<<",SSs="<<ssin2g<<",P/S="<<eteps<<",nop=" |
---|
| 104 | // <<nop<<",p="<<momb<<",e="<<enb<<",projCode="<<pPDG<<",targCode="<<tPDG<<",nEv=" |
---|
| 105 | // <<nEvt<<",fN="<<fN<<",fD="<<fD<<",cP="<<cP<<",rM="<<rM<<",sA="<<sA<<G4endl; |
---|
| 106 | //-------- Initialize CHIPS |
---|
| 107 | //G4QCHIPSWorld* theW=G4QCHIPSWorld::Get(); |
---|
| 108 | //theW->GetParticles(nop); // Create CHIPS World of nop particles |
---|
| 109 | //G4Exception("***CHIPStest: TMP"); |
---|
| 110 | //G4QNucleus::SetParameters(fN,fD,cP,rM); |
---|
| 111 | //G4Quasmon::SetParameters(temperature,ssin2g,eteps); |
---|
| 112 | //G4QEnvironment::SetParameters(sA); // SolAngle (pbar-A secondary capture) |
---|
| 113 | //-------- write results onto a file -------- |
---|
| 114 | //std::ofstream outFile( "chipstest.out", std::ios::out); |
---|
| 115 | //outFile.setf( std::ios::scientific, std::ios::floatfield ); |
---|
| 116 | |
---|
| 117 | //--------- Example how to write in the outFile ----------- |
---|
| 118 | //outFile << " " << endl; |
---|
| 119 | |
---|
| 120 | // *********** Now momb is a momentum of the incident particle ************* |
---|
| 121 | //G4double mp=G4QPDGCode(pPDG).GetMass(); // @@ just for the check |
---|
| 122 | //G4QContent pQC=G4QPDGCode(pPDG).GetQuarkContent(); |
---|
| 123 | //G4int cp=pQC.GetCharge(); |
---|
| 124 | //G4int bnp=pQC.GetBaryonNumber(); |
---|
| 125 | //G4double mp1=G4QPDGCode(PDG1).GetMass(); // @@ just for the check |
---|
| 126 | G4QContent pQC1=G4QPDGCode(PDG1).GetQuarkContent(); |
---|
| 127 | G4int cp1=pQC1.GetCharge(); |
---|
| 128 | G4int bnp1=pQC1.GetBaryonNumber(); |
---|
| 129 | //G4double mp2=G4QPDGCode(PDG2).GetMass(); // @@ just for the check |
---|
| 130 | G4QContent pQC2=G4QPDGCode(PDG2).GetQuarkContent(); |
---|
| 131 | G4int cp2=pQC2.GetCharge(); |
---|
| 132 | G4int bnp2=pQC2.GetBaryonNumber(); |
---|
| 133 | //G4double mp3=G4QPDGCode(PDG3).GetMass(); // @@ just for the check |
---|
| 134 | G4QContent pQC3=G4QPDGCode(PDG3).GetQuarkContent(); |
---|
| 135 | G4int cp3=pQC3.GetCharge(); |
---|
| 136 | G4int bnp3=pQC3.GetBaryonNumber(); |
---|
| 137 | //momb=momb; |
---|
| 138 | //G4double ep=sqrt(mp*mp+momb*momb); // @@ just for the check |
---|
| 139 | //if(enb>0.) ep=enb; |
---|
| 140 | G4double mt=G4QPDGCode(tPDG).GetMass(); // @@ just for the check |
---|
| 141 | G4QContent tQC=G4QPDGCode(tPDG).GetQuarkContent(); |
---|
| 142 | G4int ct=tQC.GetCharge(); |
---|
| 143 | G4int bnt=tQC.GetBaryonNumber(); |
---|
| 144 | G4int totC=cp1+cp2+cp3+ct; |
---|
| 145 | G4int totBN=bnp1+bnp2+bnp3+bnt; |
---|
| 146 | G4LorentzVector preSumLV(0.,0.,0.,mt); |
---|
| 147 | preSumLV+=LV1+LV2+LV3; |
---|
| 148 | G4double fEvt=nEvents; |
---|
| 149 | G4double sumE=0.; |
---|
| 150 | G4double sumK=0.; |
---|
| 151 | G4double sumG=0.; |
---|
| 152 | G4double sumT=0.; |
---|
| 153 | G4double sumN=0.; |
---|
| 154 | G4double sum0=0.; |
---|
| 155 | G4double sumP=0.; |
---|
| 156 | G4double sum1N=0.; |
---|
| 157 | G4double sumNN=0.; |
---|
| 158 | G4double sumPP=0.; |
---|
| 159 | G4double sumAL=0.; |
---|
| 160 | G4double time=clock()/CLOCKS_PER_SEC; |
---|
| 161 | // Main LOOP over events ====================================== |
---|
| 162 | for (G4int ir=0; ir<nEvents; ir++) |
---|
| 163 | { |
---|
| 164 | // Randomization loop: cycle random generator, using 2 lower digits in nEvents |
---|
| 165 | G4int iRandCount = nEvents%100; |
---|
| 166 | G4double vRandCount = 0.; |
---|
| 167 | while (iRandCount>0) |
---|
| 168 | { |
---|
| 169 | vRandCount = G4UniformRand(); |
---|
| 170 | iRandCount--; |
---|
| 171 | } |
---|
| 172 | if(!(ir%1000) && ir) G4cout<<"G4MultFragTst: "<<ir<<" events are simulated"<<G4endl; |
---|
| 173 | //G4cout<<"G4MultFragTst: "<<ir<<" events are simulated"<<G4endl; |
---|
| 174 | G4LorentzVector totSum = preSumLV; |
---|
| 175 | G4QHadron* H1 = new G4QHadron(PDG1,LV1); |
---|
| 176 | G4QHadron* H2 = new G4QHadron(PDG2,LV2); |
---|
| 177 | G4QHadron* H3 = new G4QHadron(PDG3,LV3); |
---|
| 178 | G4int totCharge = totC; |
---|
| 179 | G4int totBaryN = totBN; |
---|
| 180 | G4QHadronVector projHV; |
---|
| 181 | projHV.push_back(H1); // DESTROYED over 3 line |
---|
| 182 | projHV.push_back(H2); // DESTROYED over 2 line |
---|
| 183 | projHV.push_back(H3); // DESTROYED over 1 line |
---|
| 184 | G4QEnvironment* pan= new G4QEnvironment(projHV,tPDG); // DELETED over 8 lines |
---|
| 185 | std::for_each(projHV.begin(), projHV.end(), DeleteQHadron()); |
---|
| 186 | projHV.clear(); |
---|
| 187 | #ifdef debug |
---|
| 188 | G4cout<<"CHIPStest:===>>> Now call Fragment (HadronizeQuasmon) function" << G4endl; |
---|
| 189 | #endif |
---|
| 190 | G4QHadronVector* output; // Prototype of the output |
---|
| 191 | try |
---|
| 192 | { |
---|
| 193 | output = pan->Fragment();// DESTROYED in the end of the LOOP work space |
---|
| 194 | } |
---|
| 195 | catch (G4QException& error) |
---|
| 196 | { |
---|
| 197 | #ifdef pdebug |
---|
| 198 | G4cout<<"***CHIPStest: Exception is catched"<<G4endl; |
---|
| 199 | #endif |
---|
| 200 | G4cerr<<"***CHIPStest Abort: "<<error.GetMessage()<<G4endl; |
---|
| 201 | abort(); |
---|
| 202 | } |
---|
| 203 | #ifdef debug |
---|
| 204 | G4cout<<"CHIPStest:--->>>Now come out of Fragment (HadronizeQuasmon) function"<<G4endl; |
---|
| 205 | #endif |
---|
| 206 | delete pan; // Destruct theHadronVector (& theCandidateVector) of the Quasmon |
---|
| 207 | #ifdef debug |
---|
| 208 | G4cout << "CHIPStest: >>> Here the histograms are filled" << G4endl; |
---|
| 209 | #endif |
---|
| 210 | G4int tNH = output->size(); |
---|
| 211 | G4int npt=0; |
---|
| 212 | G4int nGamma=0; |
---|
| 213 | G4double EGamma=0; |
---|
| 214 | G4int nP0=0; |
---|
| 215 | G4int nPP=0; |
---|
| 216 | G4int nPN=0; |
---|
| 217 | G4int nKaons=0; |
---|
| 218 | G4int nEta=0; |
---|
| 219 | G4int nAlphas=0; |
---|
| 220 | G4int nPhotons=0; |
---|
| 221 | G4int nProtons=0; |
---|
| 222 | G4int nNeutrons=0; |
---|
| 223 | G4int nSpNeut=0; |
---|
| 224 | G4int nSpAlph=0; |
---|
| 225 | G4int nOmega=0; |
---|
| 226 | G4int nDec=0; |
---|
| 227 | G4int dirN=0; |
---|
| 228 | #ifdef pdebug |
---|
| 229 | G4cout<<"----------DONE^^^^^^^************^^^^^^^^^^^:ir="<<ir<<": #ofH="<<tNH<<G4endl; |
---|
| 230 | if(!(ir%100)) G4cerr<<"#"<<ir<<G4endl; |
---|
| 231 | #endif |
---|
| 232 | G4bool alarm=false; |
---|
| 233 | G4bool rad=false; |
---|
| 234 | G4bool hyp=false; |
---|
| 235 | G4bool badPDG=false; |
---|
| 236 | for (G4int ind=0; ind<tNH; ind++) |
---|
| 237 | { |
---|
| 238 | G4QHadron* curH=output->operator[](ind); |
---|
| 239 | G4double m = curH->GetMass(); // Mass of the particle |
---|
| 240 | G4LorentzVector lorV = curH->Get4Momentum(); // 4-momentum of the particle |
---|
| 241 | if(std::fabs(m-lorV.m())>.005) |
---|
| 242 | { |
---|
| 243 | G4cerr<<"***CHIPStest: m="<<lorV.m()<<" # "<<m<<", d="<<lorV.m()-m<<G4endl; |
---|
| 244 | alarm=true; |
---|
| 245 | } |
---|
| 246 | if(!(lorV.e()>=0||lorV.e()<0) || !(lorV.px()>=0||lorV.px()<0) || |
---|
| 247 | !(lorV.py()>=0||lorV.py()<0) || !(lorV.pz()>=0||lorV.pz()<0)) |
---|
| 248 | { |
---|
| 249 | G4cerr<<"***CHIPStest: NAN in LorentzVector="<<lorV<<G4endl; |
---|
| 250 | alarm=true; |
---|
| 251 | } |
---|
| 252 | G4int d=curH->GetNFragments(); // In how many particles this particle decayed |
---|
| 253 | G4ThreeVector p = lorV.vect(); // 3-momentum of the particle |
---|
| 254 | G4double e = lorV.e(); // Energy of the particle |
---|
| 255 | G4int c=curH->GetPDGCode(); // PDG Code of the particle |
---|
| 256 | //if(!d&&(c==90000002||c==90002000||c==92000000||c==221||c==331)) |
---|
| 257 | if(!d&&(c==90000002||c==90002000||c==92000000)) |
---|
| 258 | { |
---|
| 259 | //G4cout<<"***CHIPStest:***Dibaryon or Eta*** ind="<<ind<<", PDG="<<c<<G4endl; |
---|
| 260 | G4cout<<"***CHIPStest:***Dibaryon *** ind="<<ind<<", PDG="<<c<<G4endl; |
---|
| 261 | alarm=true; |
---|
| 262 | } |
---|
| 263 | if(!d&&(c==90000003||c==90003000||c==93000000)) |
---|
| 264 | { |
---|
| 265 | G4cout<<"***CHIPStest:***Tribaryon *** ind="<<ind<<", PDG="<<c<<G4endl; |
---|
| 266 | alarm=true; |
---|
| 267 | } |
---|
| 268 | if(!d) npt++; |
---|
| 269 | if(d) nDec+=d; |
---|
| 270 | if(c==223) nOmega++; |
---|
| 271 | if(c==22) nPhotons++; |
---|
| 272 | if(c==311||c==321||c==-311||c==-321) nKaons++; // kaons |
---|
| 273 | if(c==221) nEta++; // etas |
---|
| 274 | if(c==90002002) nAlphas++; // Alphas |
---|
| 275 | if(c==2212 || c==90001000) nProtons++; // Protons |
---|
| 276 | if(c==90000001 || c==90001000) dirN++; // Ditrect nucleons |
---|
| 277 | if(c==2112 || c==90000001) nNeutrons++; // Neutrons |
---|
| 278 | if((c==2112 || c==90000001) && std::fabs(e-1005.)<3.) nSpNeut++;// Dibar-Neutrons |
---|
| 279 | if(!d && c==90002002 && e-m<7.) nSpAlph++; // Special Alphas |
---|
| 280 | if(c==111) nP0++; // Neutral pions |
---|
| 281 | if(c==-211) nPN++; // Negative pions |
---|
| 282 | if(c==211) nPP++; // Positive pions |
---|
| 283 | if(c==22) nGamma++; // Gammas |
---|
| 284 | if(c==22) EGamma+=e; // Energy of gammas |
---|
| 285 | if(!d) totCharge-=curH->GetCharge(); |
---|
| 286 | if(!d) totBaryN-=curH->GetBaryonNumber(); |
---|
| 287 | if(!d) totSum -= lorV; |
---|
| 288 | if(c>80000000 && (c<90000000 || c%1000>500 || c%1000000>500000) || |
---|
| 289 | !(c>=0 || c<0)) |
---|
| 290 | { |
---|
| 291 | G4cout<<"***OUTPUT ERROR*** CHIPStest: bad PDG is found. It is "<<c<<G4endl; |
---|
| 292 | badPDG=true; |
---|
| 293 | } |
---|
| 294 | rad=rad||c==90002000||c==90003000||c==90004000||c==90000002||c==90000003||c==90000004 |
---|
| 295 | ||c==90002003||c==90003002||c==90004002||c==90002005||c==90005002||c==90004004 |
---|
| 296 | ||c==90006002; |
---|
| 297 | hyp = hyp || c>90999999; |
---|
| 298 | #ifdef pdebug |
---|
| 299 | G4cout<<"#"<<ind<<"(d="<<d<<"), PDG="<<c<<",4M="<<lorV<<m<<",T="<<lorV.e()-m<<G4endl; |
---|
| 300 | #endif |
---|
| 301 | } |
---|
| 302 | #ifdef pdebug |
---|
| 303 | G4cout<<"CHECK: 4M="<<totSum<<", Charge="<<totCharge<<", BaryN="<<totBaryN<<G4endl; |
---|
| 304 | #endif |
---|
| 305 | G4double ss=std::fabs(totSum.t())+std::fabs(totSum.x())+std::fabs(totSum.y())+std::fabs(totSum.z()); |
---|
| 306 | if (totCharge ||totBaryN || !(ss<.01) || alarm || nGamma&&!EGamma || badPDG) |
---|
| 307 | //if (totCharge || ss>.01 || alarm || nSpNeut) |
---|
| 308 | //if (totCharge || ss>.01 || alarm || nSpAlph) |
---|
| 309 | { |
---|
| 310 | G4cerr<<"***CHIPStest:#"<<ir<<":n="<<tNH<<",4M="<<totSum<<",Charge="<<totCharge |
---|
| 311 | <<",BaryN="<<totBaryN<<G4endl; |
---|
| 312 | if(nGamma&&!EGamma)G4cerr<<"***CHIPStest: Egamma=0"<<G4endl; |
---|
| 313 | totSum = preSumLV; |
---|
| 314 | for (int indx=0; indx<tNH; indx++) |
---|
| 315 | { |
---|
| 316 | G4QHadron* curH=output->operator[](indx); |
---|
| 317 | G4double m = curH->GetMass(); |
---|
| 318 | G4LorentzVector lorV = curH->Get4Momentum(); |
---|
| 319 | G4int d=curH->GetNFragments(); |
---|
| 320 | G4int c=curH->GetPDGCode(); |
---|
| 321 | if(!d) totSum -= lorV; |
---|
| 322 | G4cerr<<"#"<<indx<<"("<<d<<"), PDG="<<c<<", m/LV="<<m<<lorV<<", T="<<lorV.e()-m |
---|
| 323 | <<", d4M="<<totSum<<G4endl; |
---|
| 324 | } |
---|
| 325 | G4Exception("***CHIPStest: ALARM or charge/energy/momentum is not conserved"); |
---|
| 326 | } |
---|
| 327 | sum0+=nP0; |
---|
| 328 | sumP+=nPP; |
---|
| 329 | sumN+=nPN; |
---|
| 330 | G4int nPions=nP0+nPP+nPN; |
---|
| 331 | sumG+=nGamma; |
---|
| 332 | sumT+=EGamma; |
---|
| 333 | if(nAlphas)sumAL++; |
---|
| 334 | if(nProtons)sumPP++; |
---|
| 335 | if(nNeutrons)sumNN++; |
---|
| 336 | if(nNeutrons&&!nProtons&&!nAlphas) sum1N++; |
---|
| 337 | if(nKaons)sumK++; |
---|
| 338 | if(nEta)sumE++; |
---|
| 339 | if (nPhotons) nPions=11; |
---|
| 340 | if (nKaons==2&&npt==2) nPions=1; |
---|
| 341 | else if (nKaons) nPions=10; |
---|
| 342 | //histPi.fill(nPions); |
---|
| 343 | //histNeut.fill(nNeutrons); |
---|
| 344 | std::for_each(output->begin(), output->end(), DeleteQHadron()); |
---|
| 345 | output->clear(); |
---|
| 346 | delete output; |
---|
| 347 | } |
---|
| 348 | time=(clock()/CLOCKS_PER_SEC-time)/fEvt; |
---|
| 349 | G4cerr<<"CHIPStest::t="<<time<<",Yields:pi-="<<sumN/fEvt<<",pi+="<<sumP/fEvt<<",pi0=" |
---|
| 350 | <<sum0/fEvt<<",K="<<sumK/fEvt<<",eta="<<sumE/fEvt<<",gamma="<<sumG/fEvt<<"(<E>=" |
---|
| 351 | <<sumT/sumG<<"),n="<<sumNN/fEvt<<",p="<<sumPP/fEvt<<",alpha="<<sumAL/fEvt |
---|
| 352 | <<",onlyN="<<sum1N/fEvt<<G4endl; |
---|
| 353 | return EXIT_SUCCESS; |
---|
| 354 | } |
---|