Context Navigation

source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/calcul/G4MultyFragmentationTest.cc @ 1196

Last change on this file since 1196 was 968, checked in by garnier, 15 years ago
fichier ajoutes
File size: 14.3 KB

Line
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(mpmp+mombmomb); // @@ 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	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: