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 | } |
---|