CGEM BOSS 6.6.5.h
BESIII Offline Software System
Loading...
Searching...
No Matches
Phokhara.cxx
Go to the documentation of this file.
1//*****************************************************************************
2//
3// Phokhara.cxx
4//
5// Algorithm runs event generator Phokhara (hep-ph/0710.4227v1)
6// and stores output to transient store
7//
8// Nov. 2007 A. Zhemchugov, initial version written for BES3
9//*****************************************************************************
10
11// Header for this module:-
12#include "Phokhara/Phokhara.h"
13#include "Phokhara/PhokharaDef.h"
14
15// Framework Related Headers:-
16#include "HepMC/GenEvent.h"
17#include "HepMC/GenVertex.h"
18#include "HepMC/GenParticle.h"
19#include "CLHEP/Vector/LorentzVector.h"
20
21#include "GaudiKernel/MsgStream.h"
22#include "GaudiKernel/ISvcLocator.h"
23#include "GaudiKernel/AlgFactory.h"
24#include "GaudiKernel/DataSvc.h"
25#include "GaudiKernel/SmartDataPtr.h"
26
27#include "GeneratorObject/McGenEvent.h"
28#include "BesKernel/IBesRndmGenSvc.h"
29#include "CLHEP/Random/RandomEngine.h"
30
31#include "cfortran/cfortran.h"
32
33#include <stdlib.h>
34#include <iostream>
35#include <fstream>
36
37using namespace std;
38using namespace CLHEP;
39
40Phokhara::Phokhara(const string& name, ISvcLocator* pSvcLocator):Algorithm( name, pSvcLocator )
41{
42 m_initSeed = (0);
43
44 m_tagged = 0;
45
46 declareProperty( "InitialSeed", m_initSeed = 0);
47
48 declareProperty( "InitialEvents", m_nm = 50000 ); // # of events to determine the maximum
49 declareProperty( "NLO", m_nlo = 1 ); // Born(0), NLO(1)
50 declareProperty( "SoftPhotonCutoff", m_w = 0.0001 ); // soft photon cutoff
51 declareProperty( "Channel", m_pion = 0 ); // mu+mu-(0),pi+pi-(1),2pi0pi+pi-(2),
52 // 2pi+2pi-(3),ppbar(4),nnbar(5),
53 // K+K-(6),K0K0bar(7),pi+pi-pi0(8),
54 // Lamb Lambbar->pi-pi+ppbar(9)
55 declareProperty( "FSR", m_fsr = 1 ); // ISR only(0), ISR+FSR(1), ISR+INT+FSR(2)
56 declareProperty( "FSRNLO", m_fsrnlo = 1 ); // yes(1), no(0)
57 declareProperty( "NarrowRes", m_NarrowRes = 1 ); // none(0) jpsi(1) psip(2)
58 declareProperty( "KaonFormfactor", m_FF_Kaon = 1 ); // constrained (0), unconstrained (1), Kuhn-Khodjamirian-Bruch (2)
59 declareProperty( "VacuumPolarization", m_ivac = 0 ); // yes(1), no(0)
60 //declareProperty( "TaggedPhotons", m_tagged = 0 ); // tagged photons(0), untagged photons(1)
61 declareProperty( "PionFormfactor", m_FF_Pion = 0 ); // KS Pionformfactor(0), GS Pionformfactor old(1) and new(2)
62 declareProperty( "F0model", m_f0_model = 0 ); // f0+f0(600): KK model(0), no structure(1),
63 // no f0+f0(600)(2), f0 KLOE(3)
64 declareProperty( "Ecm", m_E = 3.097 ); // CMS-energy (GeV)
65 declareProperty( "CutTotalInvMass", m_q2min = 0.0 ); // minimal hadrons(muons)-gamma-inv mass squared (GeV^2)
66 declareProperty( "CutHadInvMass", m_q2_min_c = 0.0447 ); // minimal inv. mass squared of the hadrons(muons)(GeV^2)
67 declareProperty( "MaxHadInvMass", m_q2_max_c = 1.06 ); // maximal inv. mass squared of the hadrons(muons)(GeV^2)
68 declareProperty( "MinPhotonEnergy", m_gmin = 0.05 ); // minimal photon energy/missing energy (GeV)
69 declareProperty( "MinPhotonAngle", m_phot1cut = 0.0 ); // minimal photon angle/missing momentum angle (grad)
70 declareProperty( "MaxPhotonAngle", m_phot2cut = 180.0 ); // maximal photon angle/missing momentum angle (grad)
71 declareProperty( "MinHadronsAngle", m_pi1cut = 0.0 ); // minimal hadrons(muons) angle (grad)
72 declareProperty( "MaxHadronsAngle", m_pi2cut = 180.0 ); // maximal hadrons(muons) angle (grad)
73
74 ievent = 0;
75}
76
77StatusCode Phokhara::initialize(){
78
79 MsgStream log(messageService(), name());
80
81 log << MSG::INFO << "Phokhara initialize" << endreq;
82
83 int i,s_seed[105];
84 long int k,j;
85
86 if( m_initSeed == 0) // if seed is not set in the jobptions file, set it using BesRndmGenSvc
87 {
88 IBesRndmGenSvc* p_BesRndmGenSvc;
89 static const bool CREATEIFNOTTHERE(true);
90 StatusCode RndmStatus = service("BesRndmGenSvc", p_BesRndmGenSvc, CREATEIFNOTTHERE);
91 if (!RndmStatus.isSuccess() || 0 == p_BesRndmGenSvc)
92 {
93 log << MSG::ERROR << " Could not initialize Random Number Service" << endreq;
94 return RndmStatus;
95 }
96
97 // Save the random number seeds in the event
98 HepRandomEngine* engine = p_BesRndmGenSvc->GetEngine("PHOKHARA");
99 const long* s = engine->getSeeds();
100 m_initSeed = s[0];
101 }
102
103 log << MSG::INFO << "Set initial seed " << m_initSeed << endreq;
104
105 // --- initialize the seed ------
106 /* ifstream seeds("seed.dat");
107 if( seeds.is_open() )
108 {
109 int ii=0;
110 while( ! seeds.eof() )
111 seeds >> s_seed[ii++];
112 }
113 else
114 cerr << "Cannot open file seed.dat for reading" << endl;
115 */
116
117 RLXDINIT(1, m_initSeed);
118
119 //RLXDRESETF(s_seed);
120 //rlxd_reset(s_seed);
121
122 // --- input parameters ----------------------------
123 MAXIMA.iprint = 0;
124
125 FLAGS.nlo = m_nlo;
126 FLAGS.pion = m_pion;
127 FLAGS.fsr = m_fsr;
128 FLAGS.fsrnlo = m_fsrnlo;
129 FLAGS.ivac = m_ivac;
130 // FLAGS.tagged = m_tagged;
131 FLAGS.FF_pion = m_FF_Pion;
132 FLAGS.f0_model = m_f0_model;
133 FLAGS.FF_kaon = m_FF_Kaon;
134 FLAGS.narr_res = m_NarrowRes;
135
136 CTES.Sp = m_E*m_E; ;
137
138 CUTS.w = m_w;
139 CUTS.q2min = m_q2min;
140 CUTS.q2_min_c = m_q2_min_c;
141 CUTS.q2_max_c = m_q2_max_c;
142 CUTS.gmin = m_gmin;
143 CUTS.phot1cut = m_phot1cut;
144 CUTS.phot2cut = m_phot2cut;
145 CUTS.pi1cut = m_pi1cut;
146 CUTS.pi2cut = m_pi2cut;
147
148 INPUT();
149 // --- print run data ------------------------------
150 cout << "-------------------------------------------------------------" << endl;
151 if (FLAGS.pion == 0)
152 cout << " PHOKHARA 7.0 : e^+ e^- -> mu^+ mu^- gamma" << endl;
153 else if (FLAGS.pion == 1)
154 cout << " PHOKHARA 7.0: e^+ e^- -> pi^+ pi^- gamma" << endl;
155 else if (FLAGS.pion == 2)
156 cout << " PHOKHARA 7.0: e^+ e^- -> pi^+ pi^- 2pi^0 gamma" << endl;
157 else if (FLAGS.pion == 3)
158 cout << " PHOKHARA 7.0: e^+ e^- -> 2pi^+ 2pi^- gamma" << endl;
159 else if (FLAGS.pion == 4)
160 cout << " PHOKHARA 7.0: e^+ e^- -> p pbar gamma" << endl;
161 else if (FLAGS.pion == 5)
162 cout << " PHOKHARA 7.0: e^+ e^- -> n nbar gamma" << endl;
163 else if (FLAGS.pion == 6)
164 cout << " PHOKHARA 7.0: e^+ e^- -> K^+ K^- gamma" << endl;
165 else if (FLAGS.pion == 7)
166 cout << " PHOKHARA 7.0: e^+ e^- -> K_0 K_0bar gamma" << endl;
167 else if (FLAGS.pion == 8)
168 cout << " PHOKHARA 7.0: e^+ e^- -> pi^+ pi^- pi^0 gamma" << endl;
169 else if (FLAGS.pion == 9) {
170 cout << "PHOKHARA 7.0 : e^+ e^- ->" << endl;
171 cout << " Lambda (-> pi^- p) Lambda bar (-> pi^+ pbar) gamma" << endl;
172 } else
173 cout << " PHOKHARA 7.0: not yet implemented" << endl;
174
175 // --------------------------------
176 cout << "--------------------------------------------------------------" << endl;
177 printf(" %s %f %s\n","cms total energy = ",sqrt(CTES.Sp)," GeV ");
178// if (FLAGS.tagged == 0) {
179 if((CUTS.gmin/2.0/CTES.ebeam) < 0.0098){
180 cerr << " minimal missing energy set to small" << endl;
181 return 0;
182 }
183 printf(" %s %f %s\n","minimal tagged photon energy = ",CUTS.gmin," GeV ");
184 printf(" %s %f,%f\n","angular cuts on tagged photon = ",CUTS.phot1cut,CUTS.phot2cut);
185/* }
186 else {
187 if((CUTS.gmin/2.0/CTES.ebeam) < 0.0098){
188 cerr << " minimal missing energy set to small" << endl;
189 return 0;
190 }
191 printf(" %s %f %s\n","minimal missing energy = ",CUTS.gmin," GeV ");
192 printf(" %s %f,%f\n","angular cuts on missing momentum = ",CUTS.phot1cut, CUTS.phot2cut);
193 }
194*/
195
196 // --------------------------------
197 if (FLAGS.pion == 0)
198 printf(" %s %f,%f\n","angular cuts on muons = ",CUTS.pi1cut,CUTS.pi2cut);
199 else if (FLAGS.pion == 4)
200 printf(" %s %f,%f\n","angular cuts on protons = ",CUTS.pi1cut,CUTS.pi2cut);
201 else if (FLAGS.pion == 5)
202 printf(" %s %f,%f\n","angular cuts on neutrons = ", CUTS.pi1cut,CUTS.pi2cut);
203 else if ((FLAGS.pion == 6)||(FLAGS.pion == 7))
204 printf(" %s %f,%f\n","angular cuts on kaons = ", CUTS.pi1cut,CUTS.pi2cut);
205 else if (FLAGS.pion == 9)
206 printf(" %s %f,%f\n","angular cuts on pions and protons = ", CUTS.pi1cut,CUTS.pi2cut);
207 else
208 printf(" %s %f,%f\n","angular cuts on pions = ", CUTS.pi1cut,CUTS.pi2cut);
209
210// if (FLAGS.tagged == 0) {
211 if (FLAGS.pion == 0)
212 printf(" %s %f %s\n","min. muons-tagged photon inv.mass^2 = ", CUTS.q2min," GeV^2");
213 else if (FLAGS.pion == 4)
214 printf(" %s %f %s\n","min. protons-tagged photon inv.mass^2 = ", CUTS.q2min," GeV^2" );
215 else if (FLAGS.pion == 5)
216 printf(" %s %f %s\n","min. neutrons-tagged photon inv.mass^2 = ", CUTS.q2min," GeV^2" );
217 else if ((FLAGS.pion == 6)||(FLAGS.pion == 7))
218 printf(" %s %f %s\n","min. kaons-tagged photon inv.mass^2 = ", CUTS.q2min," GeV^2" );
219 else if (FLAGS.pion == 9)
220 printf(" %s %f %s\n"," min. lambdas-tagged photon inv.mass^2 = ", CUTS.q2min," GeV^2" );
221 else
222 printf(" %s %f %s\n","min. pions-tagged photon inv.mass^2 = ", CUTS.q2min," GeV^2" );
223// }
224// --- set cuts ------------------------------------
225// if (FLAGS.tagged == 0) {
226 cos1min = cos(CUTS.phot2cut*CTES.pi/180.0); // photon1 angle cuts in the
227 cos1max = cos(CUTS.phot1cut*CTES.pi/180.0); // LAB rest frame
228// } else {
229// cos1min = -1.0;
230// cos1max = 1.0;
231// CUTS.gmin = CUTS.gmin/2.0;
232// }
233 cos2min = -1.0; // photon2 angle limits
234 cos2max = 1.0; //
235 cos3min = -1.0; // hadrons/muons angle limits
236 cos3max = 1.0; // in their rest frame
237 if (FLAGS.pion == 0) // virtual photon energy cut
238 qqmin = 4.0*CTES.mmu*CTES.mmu;
239 else if (FLAGS.pion == 1)
240 qqmin = 4.0*CTES.mpi*CTES.mpi;
241 else if (FLAGS.pion == 2)
242 qqmin = 4.0*(CTES.mpi+CTES.mpi0)*(CTES.mpi+CTES.mpi0);
243 else if (FLAGS.pion == 3)
244 qqmin = 16.0*CTES.mpi*CTES.mpi;
245 else if (FLAGS.pion == 4)
246 qqmin = 4.0*CTES.mp*CTES.mp;
247 else if (FLAGS.pion == 5)
248 qqmin = 4.0*CTES.mnt*CTES.mnt;
249 else if (FLAGS.pion == 6)
250 qqmin = 4.0*CTES.mKp*CTES.mKp;
251 else if (FLAGS.pion == 7)
252 qqmin = 4.0*CTES.mKn*CTES.mKn;
253 else if (FLAGS.pion == 8)
254 qqmin = (2.0*CTES.mpi+CTES.mpi0)*(2.0*CTES.mpi+CTES.mpi0);
255 else if (FLAGS.pion == 9)
256 qqmin = 4.0*CTES.mlamb*CTES.mlamb;
257 // else
258 // continue;
259
260 qqmax = CTES.Sp-2.0*sqrt(CTES.Sp)*CUTS.gmin; // if only one photon
261 if (CUTS.q2_max_c < qqmax)
262 qqmax=CUTS.q2_max_c; // external cuts
263
264 // -------------------
265 if ( (CUTS.q2_min_c > qqmin) && (CUTS.q2_min_c < (CTES.Sp*(1.0-2.0*(CUTS.gmin/sqrt(CTES.Sp)+CUTS.w))) ) )
266 qqmin = CUTS.q2_min_c;
267 else {
268 cerr << "------------------------------" << endl;
269 cerr << " Q^2_min TOO SMALL" << endl;
270 cerr << " Q^2_min CHANGED BY PHOKHARA = " << qqmin << " GeV^2" << endl;
271 cerr << "------------------------------" << endl;
272 }
273 // -------------------
274 if(qqmax <= qqmin){
275 cerr << " Q^2_max to small " << endl;
276 cerr << " Q^2_max = " << qqmax << endl;
277 cerr << " Q^2_min = " << qqmin << endl;
278 return 0;
279 }
280
281 // -------------------
282 if (FLAGS.pion == 0) {
283 printf(" %s %f %s\n", "minimal muon-pair invariant mass^2 = ", qqmin," GeV^2");
284 printf(" %s %f %s\n", "maximal muon-pair invariant mass^2 = ", qqmax," GeV^2");
285 } else if (FLAGS.pion == 1) {
286 printf(" %s %f %s\n", "minimal pion-pair invariant mass^2 = ", qqmin," GeV^2");
287 printf(" %s %f %s\n", "maximal pion-pair invariant mass^2 = ", qqmax," GeV^2");
288 } else if (FLAGS.pion == 4) {
289 printf(" %s %f %s\n", "minimal proton-pair invariant mass^2 = ", qqmin," GeV^2");
290 printf(" %s %f %s\n", "maximal proton-pair invariant mass^2 = ", qqmax," GeV^2");
291 } else if (FLAGS.pion == 5) {
292 printf(" %s %f %s\n", "minimal neutron-pair invariant mass^2 = ", qqmin," GeV^2");
293 printf(" %s %f %s\n", "maximal neutron-pair invariant mass^2 = ", qqmax," GeV^2");
294 } else if ((FLAGS.pion == 6) || (FLAGS.pion == 7)) {
295 printf(" %s %f %s\n", "minimal kaon-pair invariant mass^2 = ", qqmin," GeV^2");
296 printf(" %s %f %s\n", "maximal kaon-pair invariant mass^2 = ", qqmax," GeV^2");
297 } else if(FLAGS.pion == 8){
298 printf(" %s %f %s\n", "minimal three-pion invariant mass^2 = ", qqmin," GeV^2");
299 printf(" %s %f %s\n", "maximal three-pion invariant mass^2 = ", qqmax," GeV^2");
300 } else if(FLAGS.pion == 9){
301 printf(" %s %f %s\n", "minimal lambda-pair invariant mass^2 = ", qqmin," GeV^2");
302 printf(" %s %f %s\n", "maximal lambda-pair invariant mass^2 = ", qqmax," GeV^2");
303 } else {
304 printf(" %s %f %s\n", "minimal four-pion invariant mass^2 = ", qqmin," GeV^2" );
305 printf(" %s %f %s\n", "maximal four-pion invariant mass^2 = ", qqmax," GeV^2");
306 }
307 // -------------------
308 if (FLAGS.nlo == 0) {
309 cout << "Born" << endl;
310 if(FLAGS.fsrnlo != 0){
311 cerr << "WRONG FSRNLO flag - only fsrnlo=0 allowed for Born" << endl;
312 return 0;
313 }
314 }
315 // -------------------
316 if((FLAGS.pion == 9) && (FLAGS.nlo != 0)) {
317 cerr << "WRONG NLO flag - only Born allowed for Lambdas"<< endl;
318 cerr << "If you feel that you need NLO"<< endl;
319 cerr << "please contact the authors"<< endl;
320 return 0;
321 }
322 // -------------------
323 if (FLAGS.nlo == 1) printf(" %s %f\n", "NLO: soft photon cutoff w = ",CUTS.w);
324 if ((FLAGS.pion <= 1) || (FLAGS.pion == 6))
325 {
326
327 if( ! ((FLAGS.fsr == 1) || (FLAGS.fsr == 2) || (FLAGS.fsrnlo == 0)
328 || (FLAGS.fsr == 1) || (FLAGS.fsrnlo == 1)
329 || (FLAGS.fsr == 0) || (FLAGS.fsrnlo == 0))) {
330 cerr << "WRONG combination of FSR, FSRNLO flags" <<endl;
331 return 0;
332 }
333
334 // ------------------
335 if (FLAGS.fsr == 0)
336 cout << "ISR only" << endl;
337 else if (FLAGS.fsr == 1)
338 cout << "ISR+FSR" << endl;
339 else if (FLAGS.fsr == 2) {
340 if (FLAGS.nlo == 0)
341 cout << "ISR+INT+FSR" << endl;
342 else {
343 cerr << "WRONG FSR flag: interference is included only for nlo=0" << endl;
344 return 0;
345 }
346 }
347 else {
348 cerr << "WRONG FSR flag" << FLAGS.fsr << endl;
349 return 0;
350 }
351
352 if(FLAGS.fsrnlo == 1)
353 cout << "IFSNLO included" << endl;
354 }
355 else
356 {
357 if((FLAGS.fsr == 0) && (FLAGS.fsrnlo == 0))
358 cout << "ISR only" << endl;
359 else {
360 cerr << "FSR is implemented only for pi+pi-, mu+mu- and K+K- modes" << endl;
361 return 0;
362 }
363 }
364
365 // ------------------
366 if(FLAGS.ivac == 0){
367 cout << "Vacuum polarization is NOT included" << endl;}
368 else if(FLAGS.ivac == 1){
369 cout << "Vacuum polarization by Fred Jegerlehner (http://www-com.physik.hu-berlin.de/fjeger/alphaQEDn.uu)" << endl;}
370 else if(FLAGS.ivac == 2){
371 cout << "Vacuum polarization (VP_HLMNT_v1_3nonr) by Daisuke Nomura and Thomas Teubner" << endl;}
372 else {
373 cout << "WRONG vacuum polarization switch" << endl;
374 return 0;
375 }
376
377// -----------------
378 if(FLAGS.pion == 1){
379 if(FLAGS.FF_pion == 0)
380 cout << "Kuhn-Santamaria PionFormFactor" << endl;
381 else if(FLAGS.FF_pion == 1)
382 cout << "Gounaris-Sakurai PionFormFactor old" << endl;
383 else if(FLAGS.FF_pion == 2)
384 cout << "Gounaris-Sakurai PionFormFactor new" << endl;
385 else {
386 cout << "WRONG PionFormFactor switch" << endl;
387 return 0;
388 }
389
390 if(FLAGS.fsr != 0){
391 if(FLAGS.f0_model == 0)
392 cout << "f0+f0(600): K+K- model" << endl;
393 else if(FLAGS.f0_model == 1)
394 cout << "f0+f0(600): \"no structure\" model" << endl;
395 else if(FLAGS.f0_model == 2)
396 cout << "NO f0+f0(600)" << endl;
397 else if(FLAGS.f0_model == 3)
398 cout << "only f0, KLOE: Cesare Bini-private communication" << endl;
399 else {
400 cout << "WRONG f0+f0(600) switch" << endl;
401 return 0;
402 }
403 }
404 }
405
406//-------
407 if((FLAGS.pion == 6) || (FLAGS.pion==7)){
408 if(FLAGS.FF_kaon == 0)
409 cout << "constrained KaonFormFactor" << endl;
410 else if(FLAGS.FF_kaon == 1) {
411 cout << "unconstrained KaonFormFactor" << endl;}
412 else if(FLAGS.FF_kaon == 2) {
413 cout << "Kuhn-Khodjamirian-Bruch KaonFormFactor" << endl;}
414 else{
415 cout << "WRONG KaonFormFactor switch" << endl;
416 return 0;
417 }
418 }
419// --------------------
420 if((FLAGS.pion == 0) || (FLAGS.pion==1) || (FLAGS.pion == 6) || (FLAGS.pion == 7)){
421 if(FLAGS.narr_res == 1){
422 cout << "THE NARROW RESONANCE J/PSI INCLUDED" << endl;}
423 else if(FLAGS.narr_res == 2){
424 cout << "THE NARROW RESONANCE PSI(2S) INCLUDED" << endl;}
425 else if(FLAGS.narr_res != 0){
426 cout << "wrong flag narr_res: only 0, 1 or 2 allowed" << endl;
427 return 0;
428 }
429 }
430// ------
431
432 cout << "--------------------------------------------------------------" << endl;
433//
434 // =================================================
435// --- finding the maximum -------------------------
436 for( i = 0; i<2; i++)
437 {
438 MAXIMA.Mmax[i] = 1.0;
439 MAXIMA.gross[i] = 0.0;
440 MAXIMA.klein[i] = 0.0;
441 }
442
443 if (FLAGS.nlo == 0)
444 MAXIMA.Mmax[1]=0.0; // only 1 photon events generated
445
446 MAXIMA.tr[0] = 0.0;
447 MAXIMA.tr[1] = 0.0;
448 MAXIMA.count[0] = 0.0;
449 MAXIMA.count[1] = 0.0;
450
451 // =================================================
452 // --- beginning the MC loop event generation ------
453 for(j = 1; j <= m_nm; j++)
454 {
455 RANLXDF(Ar_r,1);
456 Ar[1] = Ar_r[0];
457 if (Ar[1] <= (MAXIMA.Mmax[0]/(MAXIMA.Mmax[0]+MAXIMA.Mmax[1]))) {
458 MAXIMA.count[0] = MAXIMA.count[0]+1.0;
459 GEN_1PH(1,qqmin,qqmax,cos1min,cos1max,cos3min,cos3max);
460 }
461 else {
462 MAXIMA.count[1] = MAXIMA.count[1]+1.0;
463 GEN_2PH(1,qqmin,cos1min,cos1max,cos2min,cos2max,cos3min,cos3max);
464 }
465 }
466 // --- end of the MC loop --------------------------
467 // =================================================
468 // --- for the second run ---
469 MAXIMA.Mmax[0] = MAXIMA.gross[0]+.05*sqrt(MAXIMA.gross[0]*MAXIMA.gross[0]);
470 MAXIMA.Mmax[1] = MAXIMA.gross[1]+(.03+.02*CTES.Sp)*sqrt(MAXIMA.gross[1]*MAXIMA.gross[1]);
471 if((FLAGS.pion == 1) && (FLAGS.fsrnlo == 1))
472 MAXIMA.Mmax[1]=MAXIMA.Mmax[1]*1.5;
473 if((FLAGS.pion == 0) && (FLAGS.fsrnlo == 1))
474 MAXIMA.Mmax[1]=MAXIMA.Mmax[1]*1.5;
475
476 if((FLAGS.pion == 0) && (FLAGS.fsr == 1) && (FLAGS.fsrnlo == 0))
477 MAXIMA.Mmax[1]=MAXIMA.Mmax[1]*1.2;
478
479 if((FLAGS.pion == 2) || (FLAGS.pion == 3)){
480 MAXIMA.Mmax[0]=MAXIMA.Mmax[0]*1.1;
481 MAXIMA.Mmax[1]=MAXIMA.Mmax[1]*1.1;
482 }
483
484 if(FLAGS.pion == 8){
485 MAXIMA.Mmax[0]=MAXIMA.Mmax[0]*1.08;
486 MAXIMA.Mmax[1]=MAXIMA.Mmax[1]*1.1;
487 }
488// --- end of the second run -----------------------
489
490 MAXIMA.tr[0] = 0.0;
491 MAXIMA.tr[1] = 0.0;
492 MAXIMA.count[0] = 0.0;
493 MAXIMA.count[1] = 0.0;
494
495 return StatusCode::SUCCESS;
496}
497
498StatusCode Phokhara::execute()
499{
500 MsgStream log(messageService(), name());
501
502 int ntrials = 0;
503 int tr_old[2];
504 tr_old[0] = (int)MAXIMA.tr[0];
505 tr_old[1] = (int)MAXIMA.tr[1];
506
507 while( ntrials < 10000)
508 {
509 ievent++;
510 RANLXDF(Ar_r,1);
511 Ar[1] = Ar_r[0];
512
513 if (Ar[1] <= (MAXIMA.Mmax[0]/(MAXIMA.Mmax[0]+MAXIMA.Mmax[1]))) {
514 MAXIMA.count[0] = MAXIMA.count[0]+1.0;
515 GEN_1PH(2,qqmin,qqmax,cos1min,cos1max,cos3min,cos3max);
516 }
517 else {
518 MAXIMA.count[1] = MAXIMA.count[1]+1.0;
519 GEN_2PH(2,qqmin,cos1min,cos1max,cos2min,cos2max,cos3min,cos3max);
520 }
521
522 if( ((int)MAXIMA.tr[0]+(int)MAXIMA.tr[1]) > (tr_old[0]+tr_old[1]) ) // event accepted after cuts
523 {
524 storeParticles();
525 return StatusCode::SUCCESS;
526 }
527 ntrials ++;
528 }
529
530 log << MSG::FATAL << "Could not satisfy cuts after " << ntrials << "trials. Terminate." << endreq;
531
532 return StatusCode::FAILURE;
533}
534
535
536StatusCode Phokhara::storeParticles()
537{
538 MsgStream log(messageService(), name());
539 int npart = 0;
540
541 // Fill event information
542 GenEvent* evt = new GenEvent(1,1);
543
544 GenVertex* prod_vtx = new GenVertex();
545// prod_vtx->add_particle_out( p );
546 evt->add_vertex( prod_vtx );
547
548 // incoming beam e+
549 GenParticle* p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][0], CTES.momenta[2][0], CTES.momenta[3][0], CTES.momenta[0][0] ), -11, 3);
550 p->suggest_barcode( ++npart );
551 prod_vtx->add_particle_in(p);
552
553 // incoming beam e-
554 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][1], CTES.momenta[2][1], CTES.momenta[3][1], CTES.momenta[0][1] ), 11, 3);
555 p->suggest_barcode( ++npart );
556 prod_vtx->add_particle_in(p);
557
558
559 // gamma
560 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][2], CTES.momenta[2][2], CTES.momenta[3][2], CTES.momenta[0][2] ), 22, 1);
561 p->suggest_barcode( ++npart );
562 prod_vtx->add_particle_out(p);
563
564 if( CTES.momenta[0][3] != 0 ) // second photon exists
565 {
566 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][3], CTES.momenta[2][3], CTES.momenta[3][3], CTES.momenta[0][3] ), 22, 1);
567 p->suggest_barcode( ++npart );
568 prod_vtx->add_particle_out(p);
569 }
570
571 // other products depending on channel
572 if (FLAGS.pion == 0) { // mu+ mu-
573 // mu+
574 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][5], CTES.momenta[2][5], CTES.momenta[3][5], CTES.momenta[0][5] ), -13, 1);
575 p->suggest_barcode( ++npart );
576 prod_vtx->add_particle_out(p);
577 // mu -
578 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][6], CTES.momenta[2][6], CTES.momenta[3][6], CTES.momenta[0][6] ), 13, 1);
579 p->suggest_barcode( ++npart );
580 prod_vtx->add_particle_out(p);
581 }
582
583 if (FLAGS.pion == 1) { // pi+ pi-
584 // pi+
585 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][5], CTES.momenta[2][5], CTES.momenta[3][5], CTES.momenta[0][5] ),211, 1);
586 p->suggest_barcode( ++npart );
587 prod_vtx->add_particle_out(p);
588 // pi-
589 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][6], CTES.momenta[2][6], CTES.momenta[3][6], CTES.momenta[0][6] ), -211, 1);
590 p->suggest_barcode( ++npart );
591 prod_vtx->add_particle_out(p);
592 }
593
594 if (FLAGS.pion == 2) { // pi+ pi- pi0 pi0
595 // pi0
596 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][5], CTES.momenta[2][5], CTES.momenta[3][5], CTES.momenta[0][5] ),111, 1);
597 p->suggest_barcode( ++npart );
598 prod_vtx->add_particle_out(p);
599 // pi0
600 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][6], CTES.momenta[2][6], CTES.momenta[3][6], CTES.momenta[0][6] ), 111, 1);
601 p->suggest_barcode( ++npart );
602 prod_vtx->add_particle_out(p);
603 // pi-
604 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][7], CTES.momenta[2][7], CTES.momenta[3][7], CTES.momenta[0][7] ),-211, 1);
605 p->suggest_barcode( ++npart );
606 prod_vtx->add_particle_out(p);
607 // pi+
608 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][8], CTES.momenta[2][8], CTES.momenta[3][8], CTES.momenta[0][8] ), 211, 1);
609 p->suggest_barcode( ++npart );
610 prod_vtx->add_particle_out(p);
611 }
612
613 if (FLAGS.pion == 3) { // pi+ pi- pi+ pi-
614 // pi+
615 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][5], CTES.momenta[2][5], CTES.momenta[3][5], CTES.momenta[0][5] ),211, 1);
616 p->suggest_barcode( ++npart );
617 prod_vtx->add_particle_out(p);
618 // pi-
619 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][6], CTES.momenta[2][6], CTES.momenta[3][6], CTES.momenta[0][6] ), -211, 1);
620 p->suggest_barcode( ++npart );
621 prod_vtx->add_particle_out(p);
622 // pi-
623 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][7], CTES.momenta[2][7], CTES.momenta[3][7], CTES.momenta[0][7] ),-211, 1);
624 p->suggest_barcode( ++npart );
625 prod_vtx->add_particle_out(p);
626 // pi+
627 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][8], CTES.momenta[2][8], CTES.momenta[3][8], CTES.momenta[0][8] ), 211, 1);
628 p->suggest_barcode( ++npart );
629 prod_vtx->add_particle_out(p);
630 }
631
632 if (FLAGS.pion == 4) { // p pbar
633 // pbar
634 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][5], CTES.momenta[2][5], CTES.momenta[3][5], CTES.momenta[0][5] ),-2212, 1);
635 p->suggest_barcode( ++npart );
636 prod_vtx->add_particle_out(p);
637 // p
638 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][6], CTES.momenta[2][6], CTES.momenta[3][6], CTES.momenta[0][6] ), 2212, 1);
639 p->suggest_barcode( ++npart );
640 prod_vtx->add_particle_out(p);
641 }
642
643 if (FLAGS.pion == 5) { // n nbar
644 // nbar
645 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][5], CTES.momenta[2][5], CTES.momenta[3][5], CTES.momenta[0][5] ),-2112, 1);
646 p->suggest_barcode( ++npart );
647 prod_vtx->add_particle_out(p);
648 // n
649 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][6], CTES.momenta[2][6], CTES.momenta[3][6], CTES.momenta[0][6] ), 2112, 1);
650 p->suggest_barcode( ++npart );
651 prod_vtx->add_particle_out(p);
652 }
653
654 if (FLAGS.pion == 6) { // K+ K-
655 // pbar
656 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][5], CTES.momenta[2][5], CTES.momenta[3][5], CTES.momenta[0][5] ),321, 1);
657 p->suggest_barcode( ++npart );
658 prod_vtx->add_particle_out(p);
659 // p
660 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][6], CTES.momenta[2][6], CTES.momenta[3][6], CTES.momenta[0][6] ), -321, 1);
661 p->suggest_barcode( ++npart );
662 prod_vtx->add_particle_out(p);
663 }
664
665 if (FLAGS.pion == 7) { // K0 K0bar
666 // pbar
667 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][5], CTES.momenta[2][5], CTES.momenta[3][5], CTES.momenta[0][5] ), 311, 1);
668 p->suggest_barcode( ++npart );
669 prod_vtx->add_particle_out(p);
670 // p
671 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][6], CTES.momenta[2][6], CTES.momenta[3][6], CTES.momenta[0][6] ), -311, 1);
672 p->suggest_barcode( ++npart );
673 prod_vtx->add_particle_out(p);
674 }
675
676 if (FLAGS.pion == 8) { // pi+ pi- pi0
677 // pi+
678 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][5], CTES.momenta[2][5], CTES.momenta[3][5], CTES.momenta[0][5] ), 211, 1);
679 p->suggest_barcode( ++npart );
680 prod_vtx->add_particle_out(p);
681 // pi-
682 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][6], CTES.momenta[2][6], CTES.momenta[3][6], CTES.momenta[0][6] ), -211, 1);
683 p->suggest_barcode( ++npart );
684 prod_vtx->add_particle_out(p);
685 // pi0
686 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][7], CTES.momenta[2][7], CTES.momenta[3][7], CTES.momenta[0][7] ), 111, 1);
687 p->suggest_barcode( ++npart );
688 prod_vtx->add_particle_out(p);
689 }
690
691 if (FLAGS.pion == 9) { // Lambda Lambdabar Pi+ Pi- P Pbar
692 // Lambdabar
693 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][5], CTES.momenta[2][5], CTES.momenta[3][5], CTES.momenta[0][5] ), -3122, 2);
694 p->suggest_barcode( ++npart );
695 prod_vtx->add_particle_out(p);
696 // Lambda
697 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][6], CTES.momenta[2][6], CTES.momenta[3][6], CTES.momenta[0][6] ), 3122, 2);
698 p->suggest_barcode( ++npart );
699 prod_vtx->add_particle_out(p);
700 // pi+
701 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][7], CTES.momenta[2][7], CTES.momenta[3][7], CTES.momenta[0][7] ), 211, 1);
702 p->suggest_barcode( ++npart );
703 prod_vtx->add_particle_out(p);
704 // pbar
705 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][8], CTES.momenta[2][8], CTES.momenta[3][8], CTES.momenta[0][8] ), -2212, 1);
706 p->suggest_barcode( ++npart );
707 prod_vtx->add_particle_out(p);
708 // pi-
709 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][9], CTES.momenta[2][9], CTES.momenta[3][9], CTES.momenta[0][9] ), -211, 1);
710 p->suggest_barcode( ++npart );
711 prod_vtx->add_particle_out(p);
712 // p
713 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][10], CTES.momenta[2][10], CTES.momenta[3][10], CTES.momenta[0][10] ), 2212, 1);
714 p->suggest_barcode( ++npart );
715 prod_vtx->add_particle_out(p);
716 }
717
718 if( log.level() < MSG::INFO )
719 {
720 evt->print();
721 }
722
723 // Check if the McCollection already exists
724 SmartDataPtr<McGenEventCol> anMcCol(eventSvc(), "/Event/Gen");
725 if (anMcCol!=0)
726 {
727 // Add event to existing collection
728 MsgStream log(messageService(), name());
729 log << MSG::INFO << "Add McGenEvent to existing collection" << endreq;
730 McGenEvent* mcEvent = new McGenEvent(evt);
731 anMcCol->push_back(mcEvent);
732 }
733 else
734 {
735 // Create Collection and add to the transient store
736 McGenEventCol *mcColl = new McGenEventCol;
737 McGenEvent* mcEvent = new McGenEvent(evt);
738 mcColl->push_back(mcEvent);
739 StatusCode sc = eventSvc()->registerObject("/Event/Gen",mcColl);
740 if (sc != StatusCode::SUCCESS)
741 {
742 log << MSG::ERROR << "Could not register McGenEvent" << endreq;
743 delete mcColl;
744 delete evt;
745 delete mcEvent;
746 return StatusCode::FAILURE;
747 }
748 else
749 {
750 log << MSG::INFO << "McGenEventCol created and " << npart <<" particles stored in McGenEvent" << endreq;
751 }
752 }
753
754 return StatusCode::SUCCESS;
755
756}
757
758StatusCode Phokhara::finalize()
759{
760 MsgStream log(messageService(), name());
761 // =================================================
762 if(FLAGS.pion == 9)
763 MAXIMA.Mmax[1] = MAXIMA.Mmax[1] * (1.0 + LAMBDA_PAR.alpha_lamb)*(1.0 + LAMBDA_PAR.alpha_lamb)
764 * LAMBDA_PAR.ratio_lamb*LAMBDA_PAR.ratio_lamb;
765
766 // --- value of the cross section ------------------
767 if (FLAGS.nlo == 0) {
768 sigma = MAXIMA.Mmax[0]/MAXIMA.count[0]*MAXIMA.tr[0];
769 dsigm = MAXIMA.Mmax[0]*sqrt((MAXIMA.tr[0]/MAXIMA.count[0]-(MAXIMA.tr[0]/MAXIMA.count[0])*(MAXIMA.tr[0]/MAXIMA.count[0]))/MAXIMA.count[0]);
770 } else {
771 sigma1 = MAXIMA.Mmax[0]/MAXIMA.count[0]*MAXIMA.tr[0];
772 sigma2 = MAXIMA.Mmax[1]/MAXIMA.count[1]*MAXIMA.tr[1];
773 dsigm1 = MAXIMA.Mmax[0]*sqrt((MAXIMA.tr[0]/MAXIMA.count[0]-(MAXIMA.tr[0]/MAXIMA.count[0])*(MAXIMA.tr[0]/MAXIMA.count[0]))/MAXIMA.count[0]);
774 dsigm2 = MAXIMA.Mmax[1]*sqrt((MAXIMA.tr[1]/MAXIMA.count[1]-(MAXIMA.tr[1]/MAXIMA.count[1])*(MAXIMA.tr[1]/MAXIMA.count[1]))/MAXIMA.count[1]);
775
776 sigma = sigma1+sigma2;
777 dsigm = sqrt(dsigm1*dsigm1 + dsigm2*dsigm2);
778 }
779
780 // --- output --------------------------------------
781 cout << "-------------------------------------------------------------" << endl;
782 cout << " PHOKHARA 7.0 Final Statistics " << endl;
783 cout << "-------------------------------------------------------------" << endl;
784 cout << int(MAXIMA.tr[0]+MAXIMA.tr[1]) << " total events accepted of " << endl;
785 cout << int(ievent) << " total events generated" << endl;
786 cout << int(MAXIMA.tr[0]) << " one photon events accepted of " << endl;
787 cout << int(MAXIMA.count[0]) << " events generated" << endl;
788 cout << int(MAXIMA.tr[1]) << " two photon events accepted of " << endl;
789 cout << int(MAXIMA.count[1]) << " events generated" << endl;
790 cout << endl;
791 if (FLAGS.nlo != 0) cout << "sigma1(nbarn) = " << sigma1 << " +- " << dsigm1 << endl;
792 if (FLAGS.nlo != 0) cout << "sigma2(nbarn) = " << sigma2 << " +- " << dsigm2 << endl;
793 cout << "sigma (nbarn) = " << sigma << " +- " << dsigm << endl;
794 cout << endl;
795 cout << "maximum1 = " << MAXIMA.gross[0] << " minimum1 = " << MAXIMA.klein[0] << endl;
796 cout << "Mmax1 = " << MAXIMA.Mmax[0] << endl;
797 cout << "maximum2 = " << MAXIMA.gross[1] << " minimum2 = " << MAXIMA.klein[1] << endl;
798 cout << "Mmax2 = " << MAXIMA.Mmax[1] << endl;
799 cout << "-------------------------------------------------------------" << endl;
800
801 log << MSG::INFO << "Phokhara finalized" << endreq;
802
803 return StatusCode::SUCCESS;
804}
cos
double cos(const BesAngle a)
Definition BesAngle.h:213
INPUT
#define INPUT
Definition BesBdkRc.cxx:35
FLAGS
#define FLAGS
Definition BesBdkRc.cxx:82
s
XmlRpcServer s
Definition HelloServer.cpp:11
McGenEventCol
ObjectVector< McGenEvent > McGenEventCol
Definition McGenEvent.h:39
RLXDINIT
#define RLXDINIT(LUXURY, SEED)
Definition PhokharaDef.h:97
GEN_1PH
#define GEN_1PH(I, QQMIN, QQMAX, COS1MIN, COS1MAX, COS3MIN, COS3MAX)
Definition PhokharaDef.h:121
RANLXDF
#define RANLXDF(AR, VAL)
Definition PhokharaDef.h:118
CTES
#define CTES
Definition PhokharaDef.h:20
LAMBDA_PAR
#define LAMBDA_PAR
Definition PhokharaDef.h:61
MAXIMA
#define MAXIMA
Definition PhokharaDef.h:84
CUTS
#define CUTS
Definition PhokharaDef.h:30
GEN_2PH
#define GEN_2PH(I, QQMIN, COS1MIN, COS1MAX, COS2MIN, COS2MAX, COS3MIN, COS3MAX)
Definition PhokharaDef.h:124
IBesRndmGenSvc
manage multiple CLHEP random engines as named streams
Definition IBesRndmGenSvc.h:24
IBesRndmGenSvc::GetEngine
virtual CLHEP::HepRandomEngine * GetEngine(const std::string &StreamName)=0
Interface to the CLHEP engine.
Phokhara::Phokhara
Phokhara(const string &name, ISvcLocator *pSvcLocator)
Definition Phokhara.cxx:40
Phokhara::initialize
StatusCode initialize()
Definition Phokhara.cxx:77
Phokhara::storeParticles
StatusCode storeParticles()
Definition Phokhara.cxx:536
Phokhara::finalize
StatusCode finalize()
Definition Phokhara.cxx:758
Phokhara::execute
StatusCode execute()
Definition Phokhara.cxx:498