001
002
003
004
005
006
007
008
009
010
011
012
013
014
015
016
017
018
019
020
021
022
023
024
025
026
027
028
029
030
031
032
033
034
035
036 #ifdef __GNUG__
037 #pragma implementation
038 #endif
039
040 #include <sstream>
041 #include "inc/misc.hh"
042 #include "generatorMods/FakeEv.hh"
043 #include "evt/Event.hh"
044 #include "stdhep_i/CdfHepevt.hh"
045 #include "ParticleDB/ParticleDb.hh"
046 #include "Framework/APPFramework.hh"
047
048
049 #include "r_n/CdfRn.hh"
050 #include "CLHEP/Random/RandFlat.h"
051 #include "CLHEP/Random/RandGaussT.h"
052 #include "CLHEP/Random/RandExponential.h"
053 #include "CLHEP/Random/RandomEngine.h"
054
055 const long FakeEv::_defaultRandomSeed1 = 922883591;
056 const long FakeEv::_defaultRandomSeed2 = 109735476;
057 const char* FakeEv::genId = "FAKE_EVENT";
058 HepRandomEngine* FakeEv::fakeEvEngine = 0;
059
060 using std::ostream;
061 using std::cout;
062 using std::endl;
063 using std::ostringstream ;
064
065
066 FakeEv::FakeEv() :
067 AbsGenModule( FakeEv::genId, "Single Particle Gun Module" ),
068 _useCommand("use",this),
069 _genCommand("generate",this),
070 _randomSeed1("RandomSeed1" ,this,FakeEv::_defaultRandomSeed1),
071 _randomSeed2("RandomSeed2" ,this,FakeEv::_defaultRandomSeed2),
072 _antiParticle("antiParticle",this,0),
073 _genMom(50.,0., 1.,100.,0.,1),
074 _genPt(50.,0., 1.,100.,0.,1),
075 _genTheta(90.,0.,-45., 45.,0.,1),
076 _genEta(0.0,0., -1., 1.,0.,2),
077 _genY(0.0,0., -1., 1.,0.,2),
078 _genPhi(7.5,0., 0.,360.,0.,2),
079 _usePtHist("usePtHist",this,false),
080 _lowPt("lowPt",this,0.),
081 _highPt("highPt",this,0.,0.),
082 _nbinPt("nbinPt",this,0),
083 _valPt("valPt",this,1,200)
084 {
085 commands()->append(&_useCommand);
086 commands()->append(&_genCommand);
087 commands()->append(&_antiParticle);
088 _initializeRandomTalkTo();
089 _initializeUsePtHistTalkTo();
090
091
092
093
094
095
096
097
098
099
100
101 NParticles = 1;
102
103 UsePt = 1;
104
105 PolarChoice = NOCHOICE;
106
107 CdfCode = 401;
108 memset(_myValPt,0,200*sizeof(double));
109 }
110
111 FakeEv::~FakeEv() {
112 }
113
114
115
116 static int process_show_command(FakeEv* module, int argc, char** argv) {
117 char cmd[100];
118
119 for (int i=0; i<argc; i++) {
120 for (int j=0; j<=strlen(argv[i]); j++) cmd[j]=toupper(argv[i][j]);
121
122 if (strcmp(cmd,"P") == 0) {
123 module->genMom()->print();
124 }
125 else if (strcmp(cmd,"PT") == 0) {
126 module->genPt()->print();
127 }
128 else if (strcmp(cmd,"THETA") == 0) {
129 module->genTheta()->print();
130 }
131 else if (strcmp(cmd,"ETA") == 0) {
132 module->genEta()->print();
133 }
134 else if (strcmp(cmd,"PHI") == 0) {
135 module->genPhi()->print();
136 }
137 else if (strcmp(cmd,"Y") == 0) {
138 module->genY()->print();
139 }
140 else if (strcmp(cmd,"ALL") == 0) {
141 module->genMom()->print();
142 module->genPt()->print();
143 module->genTheta()->print();
144 module->genEta()->print();
145 module->genPhi()->print();
146 module->genY()->print();
147 std::cout << "CdfCode= " << module->getCdfCode() << std::endl;
148 }
149 }
150 return 1;
151 }
152
153 static int process_set_command(FakeEv* module, int argc, char** argv) {
154
155
156 FakeEv::GenPars* genpar;
157 float* par;
158 int* ipar;
159 char cmd[100];
160 int len_cmd;
161
162 if (argc <= 0) return -1;
163 len_cmd = strlen(argv[0]);
164
165 for (int j=0; j<=len_cmd; j++) cmd[j]=toupper(argv[0][j]);
166
167 if (strcmp(cmd,"P") == 0) {
168
169 genpar = module->genMom();
170 }
171 else if (strcmp(cmd,"PT") == 0) {
172 genpar = module->genPt();
173 }
174 else if (strcmp(cmd,"THETA") == 0) {
175 genpar = module->genTheta();
176 }
177 else if (strcmp(cmd,"Y") == 0) {
178 genpar = module->genY();
179 }
180 else if (strcmp(cmd,"ETA") == 0) {
181 genpar = module->genEta();
182 }
183 else if (strcmp(cmd,"PHI") == 0) {
184 genpar = module->genPhi();
185 }
186 else if (strcmp(cmd,"CDFCODE") == 0) {
187 if(argc>1) {
188 int code = atoi(argv[1]);
189 module->setCdfCode(code);
190 return 1 ;
191 }
192 }
193 else if (strncmp(cmd,"NPARTICLES",len_cmd) == 0) {
194 if(argc>1) {
195 int np = atoi(argv[1]);
196 module->setNParticles(np);
197 return 1;
198 }
199 }
200 else {
201 std::cout << "process_set_command : unknown parameter" << std::endl;
202 return -2;
203 }
204
205 par = (float*) genpar;
206 ipar = (int* ) genpar;
207
208 if(argc<6) {
209 std::cout << " Not all parameters specified: Use show to check status\n"
210 << std::endl;
211 }
212 for (int i=1; i<argc; i++) {
213 if (i < 6) {
214
215
216 par[i-1] = (float) atof(argv[i]);
217 }
218 else if (i == 6) {
219
220
221 ipar[i-1] = atoi(argv[i]);
222 if(ipar[i-1]<1 || ipar[i-1]>5) {
223 std::cout << " illegal mode " << ipar[i-1] << " set to " << 1 << std::endl;
224 ipar[i-1] = 1;
225 }
226 }
227 else {
228 std::cout << "Too many parameters given: only 1st 6 used\n";
229 }
230 }
231 if (par[3] <= par[2] ) {
232 std::cout << " FAKEEV: generated Min value is larger than Max value:"
233 << " min: " << par[2] << " max: " << par[3]
234 << ". Reseting values to min = " << par[2];
235
236 par[3]=par[2]+1.0e-5;
237 std::cout << ", max = " << par[3] << std::endl;
238 }
239 return 1;
240 }
241
242 static void process_use_command(FakeEv* module, int argc, char** argv) {
243 char cmd[100];
244 for (int i=0; i<argc; i++) {
245 for (int j=0; j<=strlen(argv[i]); j++) cmd[j]=toupper(argv[i][j]);
246 if (strcmp(cmd,"PT") == 0) {
247 module->setUsePt(1);
248 }
249 else if (strcmp(cmd,"P") == 0) {
250 module->setUsePt(0);
251 }
252 else if (strcmp(cmd,"ETA") == 0) {
253 if (module->useY()) {
254 cout << " \t \t FAKE_EVENT "
255 << " Using Eta overrides Using Rapidity " << endl;
256 }
257 else if (module->useTheta()) {
258 cout << " \t \t FAKE_EVENT "
259 << " Using Eta overrides Using Theta " << endl;
260 }
261 module->setUseEta();
262 }
263 else if (strcmp(cmd,"THETA") == 0) {
264 if (module->useY()) {
265 cout << " \t \t FAKE_EVENT "
266 << " Using Theta overrides Using Rapidity " << endl;
267 }
268 else if (module->useEta()) {
269 cout << " \t \t FAKE_EVENT "
270 << " Using Theta overrides Using Eta " << endl;
271 }
272 module->setUseTheta();
273 }
274 else if (strcmp(cmd,"Y") == 0) {
275 if (module->useEta()) {
276 cout << " \t \t FAKE_EVENT "
277 << " Using Rapidity overrides Using Eta " << endl;
278 } else if (module->useTheta()) {
279 cout << " \t \t FAKE_EVENT "
280 << " Using Rapidity overrides Using Theta " << endl;
281 }
282 module->setUseY();
283 }
284 }
285 }
286
287 void FakeEv::menuHandler(char* menu, char* command, int argc, char** argv) {
288 char cmd[100];
289 for (int j=0; j<=strlen(command); j++) cmd[j]=toupper(command[j]);
290 if (strcmp(cmd,"SHOW") == 0) {
291 process_show_command(this,argc,argv);
292 }
293 else if (strcmp(cmd,"SET") == 0) {
294 process_set_command(this,argc,argv);
295 }
296 else if (strcmp(cmd,"USE") == 0) {
297 process_use_command(this,argc,argv);
298 }
299 }
300
301 int FakeEv::callGenerator(AbsEvent* event) {
302 double pt, p, eta, theta, phi;
303 double px, py, pz, rn;
304 HepRandomEngine* engine = FakeEv::fakeEvEngine;
305
306 for (int i=0; i<NParticles; i++) {
307
308
309
310
311 int code = CdfCode;
312 if(!ParticleDb::Instance()->Particle(abs(code))){
313 ERRLOG(ELabort,"[FAKEV_BAD_PARTICLE]")
314 << "FakeEv::callGenerator: Particle CdfCode "<<code
315 <<" was not found in database"<<endmsg;
316 return AppResult::ERROR;
317 }
318 if (_antiParticle.value() > 0) {
319
320
321
322 rn = RandFlat::shoot(engine,0.,1.);
323 if (rn < _antiParticle.value()) {
324 code = -CdfCode;
325 }
326 }
327
328 double mass = ParticleDb::Instance()->Mass(code);
329
330 if(useY()) {
331 eta = generateValue(_genY);
332 }
333 else if(useEta()) {
334 eta = generateValue(_genEta);
335 theta = 2.*atan(exp(-eta));
336 }
337 else {
338 theta = DEGRAD*generateValue(_genTheta);
339 }
340
341
342 phi = DEGRAD*generateValue(_genPhi);
343
344 if(UsePt) {
345 if(_usePtHist.value()) {
346 pt = generatePtUsingHist();
347 if (verbose()) Plot("pt", (float)pt, _nbinPt.value(),
348 (float)_lowPt.value(),
349 (float)_highPt.value());
350 }
351 else {
352 pt = generateValue(_genPt); }
353 if (useY()) {
354 double transverseMass = sqrt(pt*pt + mass*mass);
355 theta = atan((pt/(transverseMass*sinh(eta))));
356 }
357 p = pt/sin(theta);
358 }
359 else {
360 p = generateValue(_genMom);
361 if (useY()) {
362 double transverseMass = sqrt(p*p + mass*mass)/cosh(eta);
363 theta = atan( sqrt(transverseMass*transverseMass - mass*mass)/(transverseMass*sinh(eta)));
364 }
365 pt = p*sin(theta);
366 }
367
368
369 px = pt*cos(phi);
370 py = pt*sin(phi);
371 pz = p*cos(theta);
372
373 Hepevt_t* hepevt = CdfHepevt::Instance()->HepevtPtr();
374
375 hepevt->IDHEP[i] = ParticleDb::Instance()->ParticlePdgCode(code);
376
377
378
379 hepevt->ISTHEP[i] = 1;
380 hepevt->JMOHEP[i][0] = 0;
381 hepevt->JMOHEP[i][1] = 0;
382
383 hepevt->JDAHEP[i][0] = 0;
384 hepevt->JDAHEP[i][0] = 0;
385
386 hepevt->PHEP[i][0] = px;
387 hepevt->PHEP[i][1] = py;
388 hepevt->PHEP[i][2] = pz;
389 hepevt->PHEP[i][3] = sqrt( px*px + py*py + pz*pz + mass*mass );
390 hepevt->PHEP[i][4] = mass;
391
392 hepevt->VHEP[i][0] = 0.;
393 hepevt->VHEP[i][1] = 0.;
394 hepevt->VHEP[i][2] = 0.;
395 hepevt->VHEP[i][3] = 0.;
396 }
397
398 CdfHepevt::Instance()->HepevtPtr()->NHEP = NParticles;
399 return 0;
400 }
401
402
403 void FakeEv::fillHepevt() {}
404
405
406 void FakeEv::PrintState(std::ostream& output) const {
407 output << " " << std::endl;
408 output << "----- " << name( ) << " Parameters of Generation" << std::endl;
409 output << " CdfCode " << CdfCode << std::endl;
410 output << " NParticles " << NParticles << std::endl;
411
412 if (UsePt == 1) {
413 output << " Pt "; _genPt.print(output);
414 }
415 else {
416 output << " Momentum "; _genMom.print(output);
417 }
418 if (PolarChoice == USE_ETA){
419 output << " Eta "; _genEta.print(output);
420 }
421 else if (PolarChoice == USE_Y){
422 output << " Y "; _genY.print(output);
423 }
424 else if (PolarChoice == USE_THETA) {
425 output << " Theta "; _genTheta.print(output);
426 }
427 else {
428 output << " No Y/ETA/THETA Choice made. Will use ETA as default " << endl
429 << " Eta "; _genEta.print(output);
430 }
431 output << " Phi "; _genPhi.print(output);
432 }
433
434
435 FakeEvCommand::FakeEvCommand() {
436 }
437
438
439 FakeEvCommand::~FakeEvCommand(){
440 }
441
442
443 FakeEvCommand::FakeEvCommand(const char* const theCommand,
444 FakeEv* theTarget) : APPCommand( theCommand, theTarget)
445 {
446 _thisModule = theTarget;
447 }
448
449
450 int FakeEvCommand::handle(int argc, char* argv[])
451 {
452 int rc = 0;
453 if(strncmp(command(),"generate",8)==0) {
454 target()->menuHandler("FAKE_EVENT","SET",argc-1,argv+1);
455 }
456 else if(strncmp(command(),"use",3)==0) {
457 target()->menuHandler("FAKE_EVENT","USE", argc-1, argv+1);
458 }
459 else {
460 std::cout << " unknown command " << command() << std::endl;
461 }
462
463 return rc;
464
465 }
466
467
468 void FakeEv::Set(char* nm, float mean, float sigma, float pmax, float pmin,
469 float power, int mode) {
470
471 char name[100];
472 GenPars* par;
473
474 for (int i=0; i<=strlen(nm); i++) name[i] = toupper(nm[i]);
475
476 if (strcmp(name,"P") == 0) {
477 par = &_genMom;
478 UsePt = 0;
479 }
480 else if (strcmp(name,"PT") == 0) {
481 par = &_genPt;
482 UsePt = 1;
483 }
484 else if (strcmp(name,"THETA") == 0) {
485 par = &_genTheta;
486 setUseTheta();
487 }
488 else if (strcmp(name,"ETA") == 0) {
489 par = &_genEta;
490 setUseEta();
491 }
492 else if (strcmp(name,"Y") == 0) {
493 par = &_genY;
494 setUseY();
495 }
496 else if (strcmp(name,"PHI") == 0) {
497 par = &_genPhi;
498 }
499 else {
500 std::cout << " FakeEv::Set detected wrong parameter name : " << nm << std::endl;
501 return;
502 }
503 par->init(mean,sigma,pmin,pmax,power,mode);
504 }
505
506
507
508 void FakeEvCommand::show() const {
509
510 if(strncmp(command(),"generate",8)==0) {
511 _thisModule->PrintState();
512 }
513 }
514
515
516 bool FakeEvCommand::isShowable() const {return true;}
517
518
519 string FakeEvCommand::description() const {
520 string retval;
521 if(strncmp(command(),"generate",8)==0) {
522 retval += "possible syntaxes are:: \n";
523 retval += "\t\t generate PARAM_NAME mean sigma pmin pmax power mode\n";
524 retval += "\t\t where PARAM_NAME is one of P, PT, ETA, THETA,Y, PHI\n";
525 retval += "\t\t\t mode=1 means Gaussian (uses mean, sigma) \n";
526 retval += "\t\t\t mode=2 means power law (uses mean, pmin, pmax, power)\n";
527 retval += "\t\t\t mode=4 means exponential(uses mean, pmin, pmax) \n";
528 retval += "\t\t note: PHI is in degrees not radians!!\n";
529 retval += "\t OR: \n";
530 retval += "\t\t generate CDFCODE code\n";
531 retval += "\t OR: \n";
532 retval += "\t\t NPARTICLES number_of_particles";
533 }
534 else if(strncmp(command(),"use",3)==0) {
535 retval += "specify if generation is in P or PT and/or THETA or ETA or Y\n";
536 retval += "\t\t syntax is: use variable\n";
537 retval += "\t\t where variable means P, PT, THETA, ETA, Y\n";
538 }
539 else {
540 std::cout << " unknown command " << command() << std::endl;
541 }
542
543 return retval;
544 }
545
546
547 AppResult FakeEv::genBeginRun(AbsEvent* run) {
548 return AppResult::OK;
549 }
550
551
552 AppResult FakeEv::genEndRun(AbsEvent* run) {
553 return AppResult::OK;
554 }
555
556
557 AppResult FakeEv::genBeginJob() {
558
559 FakeEv::fakeEvEngine = CdfRn::Instance()->GetEngine("FAKE_EVENT");
560
561 CdfRn* rn = CdfRn::Instance();
562 if ( !rn->isReadingFromFile() ) {
563 rn->SetEngineSeeds(_randomSeed1.value(), _randomSeed2.value(),FakeEv::genId);
564 }
565 if (PolarChoice == NOCHOICE) {
566 std::cout << " FAKE_EVENT: No Choice of ETA/THETA/Y"
567 << std::endl << " Setting default to ETA " << std::endl;
568 setUseEta();
569 }
570
571 if (_usePtHist.value() ) {
572 if ( _nbinPt.value() > 200 ) {
573 std::cout << " FAKE_EVENT: usePtHist is used; \n Tcl error: Maximum number of bins is 200. \n";
574 return AppResult::ERROR;
575 }
576 if ( _lowPt.value() >= _highPt.value() ) {
577 std::cout << " FAKE_EVENT: usePtHist is used; \n Tcl error: low Pt value is greater or equal then high Pt .\n";
578 return AppResult::ERROR;
579 }
580 if ( _nbinPt.value() <= 0 ) {
581 std::cout << " FAKE_EVENT: usePtHist is used; \n Tcl error: zero or negative number of pt bins .\n";
582 return AppResult::ERROR;
583 }
584 double min=0,max=0;
585 int j = 0;
586 for (AbsParmList<double>::ConstIterator i=_valPt.begin();
587 i != _valPt.end(); i++) {
588 min=min <= (*i) ? min : (*i);
589 max=max >= (*i) ? max : (*i);
590 _myValPt[j]=(*i);
591 j++;
592 }
593 if (min<0.) {
594 std::cout << " FAKE_EVENT: usePtHist is used; \n Tcl error: negative Pt value in the histogram .\n";
595 return AppResult::ERROR;
596 }
597 if (min==max && max==0.) {
598 std::cout << " FAKE_EVENT: usePtHist is used; \n Tcl error: zero Pt histogram .\n";
599 return AppResult::ERROR;
600 }
601
602 }
603 return AppResult::OK;
604 }
605
606
607 AppResult FakeEv::genEndJob() {
608 return AppResult::OK;
609 }
610
611 float FakeEv::generateValue(const GenPars& pars) const {
612 double a, b, p, tmp, z, val;
613 HepRandomEngine* engine = FakeEv::fakeEvEngine;
614 switch (pars.mode()) {
615 case 1:
616 if (!pars.sigma())
617 return pars.mean();
618 else {
619 tmp = 0.0;
620 do {
621 tmp = RandGaussT::shoot( engine, pars.mean(), pars.sigma() );
622 } while
623 ((pars.min() || pars.max()) && (tmp < pars.min() || tmp > pars.max()));
624 return (float) tmp;
625 }
626 case 2:
627 case 3:
628 if (!pars.power()) {
629 return RandFlat::shoot(engine,pars.min(),pars.max());
630 }
631 else
632 {
633 p = pars.power()+1.0;
634 if(abs(p) < 1.e-8) {
635 std::cerr<< " FakeEv:generateValue() cant generate power=-1" << std::endl;
636 return 0;
637 }
638 a = pow( (double) pars.min(),p);
639 b = pow( (double) pars.max(),p)-a;
640 tmp = RandFlat::shoot(engine,0.,1.);
641 z = a+b*tmp;
642 val = pow(z,1./p);
643 return (float) val;
644 }
645 case 4:
646 {
647 tmp = 0.0;
648 do {
649 tmp = RandExponential::shoot( engine, pars.mean());
650 } while
651 ((pars.min() || pars.max()) && (tmp < pars.min() || tmp > pars.max()));
652 return (float) tmp;
653 }
654 case 5:
655
656
657
658 if (!pars.power()) {
659 return RandFlat::shoot(engine,pars.min(),pars.max());
660 }
661 else {
662 p = pars.power()+1.0;
663 if(abs(p) < 1.e-8) {
664 std::cerr<< " FakeEv:generateValue() cant generate power=-1" << std::endl;
665 return 0;
666 }
667 a = pow((double) (pars.min()-pars.mean()),p);
668 b = pow((double) (pars.max()-pars.mean()),p)-a;
669 tmp = RandFlat::shoot(engine,0.,1.);
670 z = a+b*tmp;
671 val = pars.mean()+pow(z,1./p);
672 return (float) val;
673 }
674 default:
675 std::cerr << "FakeEv::generateValue(): Invalid mode "
676 << pars.mode() << "!" << std::endl;
677 return 0.0;
678 };
679 }
680
681 void FakeEv::_initializeRandomTalkTo(void) {
682
683
684 _randomNumberMenu.initialize("RandomNumberMenu",this);
685 _randomNumberMenu.initTitle("Fake event random number menu");
686 commands()->append(&_randomNumberMenu);
687
688
689 _randomNumberMenu.commands()->append(&_randomSeed1);
690 _randomNumberMenu.commands()->append(&_randomSeed2);
691
692 ostringstream tmpSstream1;
693 ostringstream tmpSstream2;
694
695
696 tmpSstream1 << " \t\t\tSeed #1 for the random number generator"
697 << "\n\t\t\t(default " << _randomSeed1.value() << ").";
698 _randomSeed1.addDescription(tmpSstream1.str());
699 tmpSstream2 << " \t\t\tSeed #2 for the random number generator"
700 << "\n\t\t\t(default " << _randomSeed2.value() << ").";
701 _randomSeed2.addDescription(tmpSstream2.str());
702
703 }
704
705 void FakeEv::_initializeUsePtHistTalkTo(void) {
706
707
708 _PtHistMenu.initialize("PtHistMenu",this);
709 _randomNumberMenu.initTitle("FakeEvent Pt Histogram Menu");
710 commands()->append(&_PtHistMenu);
711
712
713 _PtHistMenu.commands()->append(&_usePtHist);
714 _PtHistMenu.commands()->append(&_lowPt);
715 _PtHistMenu.commands()->append(&_highPt);
716 _PtHistMenu.commands()->append(&_nbinPt);
717 _PtHistMenu.commands()->append(&_valPt);
718 }
719
720 void FakeEv::GenPars::print(std::ostream& os) const {
721 os << " Mean=" << _mean
722 << " Sigma=" << _sigma
723 << " Min=" << _min
724 << " Max=" << _max
725 << " Power=" << _power
726 << " Mode= " << _mode << std::endl;
727 }
728
729 double FakeEv::generatePtUsingHist() {
730
731 double* val = _myValPt;
732 HepRandomEngine* engine = FakeEv::fakeEvEngine;
733 return CdfRn::CdfHrndm1(_lowPt.value(),_highPt.value(),
734 _nbinPt.value(),val,engine);
735 }
Send problems or questions to cdfcode@fnal.gov
