|
|
|
#include <GeneratorInterface/MuEnrichInterface/interface/PythiaMuEnrichSource.h>
Definition at line 71 of file PythiaMuEnrichSource.h.
| edm::PythiaMuEnrichSource::PythiaMuEnrichSource | ( | const ParameterSet & | pset, | |
| const InputSourceDescription & | desc | |||
| ) |
Definition at line 56 of file PythiaMuEnrichSource.cc.
References call_pygive(), call_txgive(), call_txgive_init(), comenergy, edm::errors::Configuration, lat::endl(), firstevent, firstRun, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), i, initdc(), LENGTH, lmu_ptmin, LogTrace, maxeta, maxEventsToPrint_, nbc, nbcsim, nevnt, ngenbfe, ngenbgs, ngenbpf, ngencfe, ngencgs, ngencpf, nrejbc, nselbfe, nselbgs, nselbpf, nselcfe, nselcgs, nselcpf, pars, pydat1, pythiaHepMCVerbosity_, pythiaPylistVerbosity_, RADIUS, sumwg, and sumwg2.
00057 : 00058 GeneratedInputSource(pset, desc), evt(0), 00059 pythiaPylistVerbosity_ (pset.getUntrackedParameter<int>("pythiaPylistVerbosity",0)), 00060 pythiaHepMCVerbosity_ (pset.getUntrackedParameter<bool>("pythiaHepMCVerbosity",false)), 00061 maxEventsToPrint_ (pset.getUntrackedParameter<int>("maxEventsToPrint",1)), 00062 maxeta(pset.getUntrackedParameter<double>("maxeta", 2.4)), 00063 comenergy(pset.getUntrackedParameter<double>("comEnergy",14000.)), 00064 wtype(pset.getUntrackedParameter<int>("wtype", 1)), 00065 nmu_min(pset.getUntrackedParameter<int>("nmu_min", 1)), 00066 wbc(pset.getUntrackedParameter<double>("wbc", 1.)), 00067 wmin(pset.getUntrackedParameter<double>("wmin", 0.1)), 00068 wprintmin(pset.getUntrackedParameter<double>("wprintmin", 0.5)), 00069 lmu_ptmin(pset.getUntrackedParameter<double>("ptmin", 3.)) 00070 { 00071 00072 //now do what ever initialization is needed 00073 HepMC::HEPEVT_Wrapper::set_max_number_entries(4000); 00074 HepMC::HEPEVT_Wrapper::set_sizeof_real(8); 00075 nevnt=0; 00076 sumwg2=0; 00077 sumwg=0; 00078 firstRun=-1; 00079 LogTrace ("PythiaMuEnrichSource") << "PyhtiaMuEnrichSource: initializing Pythia. " << std::endl; 00080 00081 // PYLIST Verbosity Level 00082 // Valid PYLIST arguments are: 1, 2, 3, 5, 7, 11, 12, 13 00083 pythiaPylistVerbosity_ = pset.getUntrackedParameter<int>("pythiaPylistVerbosity",0); 00084 LogTrace ("PythiaMuEnrichSource") << "Pythia PYLIST verbosity level = " << pythiaPylistVerbosity_ << std::endl; 00085 00086 // HepMC event verbosity Level 00087 pythiaHepMCVerbosity_ = pset.getUntrackedParameter<bool>("pythiaHepMCVerbosity",false); 00088 LogTrace ("PythiaMuEnrichSource") << "Pythia HepMC verbosity = " << pythiaHepMCVerbosity_ << std::endl; 00089 00090 //Max number of events printed on verbosity level 00091 maxEventsToPrint_ = pset.getUntrackedParameter<int>("maxEventsToPrint",0); 00092 LogTrace ("PythiaMuEnrichSource") << "Number of events to be printed = " << maxEventsToPrint_ << std::endl; 00093 00094 // Set PYTHIA parameters in a single ParameterSet 00095 ParameterSet pythia_params = 00096 pset.getParameter<ParameterSet>("PythiaParameters") ; 00097 00098 // The parameter sets to be read (default, min bias, user ...) in the 00099 // proper order. 00100 std::vector<std::string> setNames = pythia_params.getParameter<std::vector<std::string> >("parameterSets"); 00101 00102 // Loop over the sets 00103 for ( unsigned i=0; i<setNames.size(); ++i ) { 00104 00105 std::string mySet = setNames[i]; 00106 00107 // Read the PYTHIA parameters for each set of parameters 00108 std::vector<std::string> pars = 00109 pythia_params.getParameter<std::vector<std::string> >(mySet); 00110 00111 if (mySet != "CSAParameters"){ 00112 LogTrace ("PythiaMuEnrichSource") << "----------------------------------------------" << std::endl; 00113 LogTrace ("PythiaMuEnrichSource") << "Read PYTHIA parameter set " << mySet << std::endl; 00114 LogTrace ("PythiaMuEnrichSource") << "----------------------------------------------" << std::endl; 00115 00116 // Loop over all parameters and stop in case of mistake 00117 for( std::vector<std::string>::const_iterator 00118 itPar = pars.begin(); itPar != pars.end(); ++itPar ) { 00119 static std::string sRandomValueSetting("MRPY(1)"); 00120 if( 0 == itPar->compare(0,sRandomValueSetting.size(),sRandomValueSetting) ) { 00121 throw edm::Exception(edm::errors::Configuration,"PythiaError") 00122 <<" attempted to set random number using pythia command 'MRPY(1)' this is not allowed.\n Please use the RandomNumberGeneratorService to set the random number seed."; 00123 } 00124 if( ! call_pygive(*itPar) ) { 00125 throw edm::Exception(edm::errors::Configuration,"PythiaError") 00126 <<" pythia did not accept the following \""<<*itPar<<"\""; 00127 } 00128 } 00129 }else if(mySet == "CSAParameters"){ 00130 00131 // Read CSA parameter 00132 00133 pars = pythia_params.getParameter<std::vector<std::string> >("CSAParameters"); 00134 00135 LogTrace ("PythiaMuEnrichSource") << "----------------------------------------------" << std::endl; 00136 LogTrace ("PythiaMuEnrichSource") << "Reading CSA parameter settings. " << std::endl; 00137 LogTrace ("PythiaMuEnrichSource") << "----------------------------------------------" << std::endl; 00138 call_txgive_init(); 00139 00140 00141 // Loop over all parameters and stop in case of a mistake 00142 for (std::vector<std::string>::const_iterator 00143 itPar = pars.begin(); itPar != pars.end(); ++itPar) { 00144 call_txgive(*itPar); 00145 00146 } 00147 00148 } 00149 } 00150 00151 //pypars.mstp[141]=2; 00152 pydat1.mstj[21]=1; 00153 pydat1.parj[72]=RADIUS; 00154 pydat1.parj[73]=LENGTH; 00155 nbc=nbcsim=nrejbc=0; 00156 ngenbfe=ngenbpf=ngencfe=ngencpf=ngenbgs=ngencgs=0; 00157 nselbfe=nselbpf=nselcfe=nselcpf=nselbgs=nselcgs=0; 00158 firstevent=true; 00159 initdc(&lmu_ptmin,&maxeta); 00160 00161 //In the future, we will get the random number seed on each event and tell 00162 // pythia to use that new seed 00163 LogTrace ("PythiaMuEnrichSource") << "----------------------------------------------" << std::endl; 00164 LogTrace ("PythiaMuEnrichSource") << "Setting Pythia random number seed " << std::endl; 00165 LogTrace ("PythiaMuEnrichSource") << "----------------------------------------------" << std::endl; 00166 edm::Service<RandomNumberGenerator> rng; 00167 uint32_t seed = rng->mySeed(); 00168 std::ostringstream sRandomSet; 00169 sRandomSet <<"MRPY(1)="<<seed; 00170 call_pygive(sRandomSet.str()); 00171 call_pyinit( "CMS", "p", "p", comenergy ); 00172 //******** 00173 00174 produces<HepMCProduct>(); 00175 LogTrace ("PythiaMuEnrichSource") << "PyhtiaMuEnrichSource: starting event generation ... " << std::endl; 00176 }
| edm::PythiaMuEnrichSource::~PythiaMuEnrichSource | ( | ) | [virtual] |
Definition at line 179 of file PythiaMuEnrichSource.cc.
References clear(), lat::endl(), nbcsim, nevnt, ngenbfe, ngenbgs, ngenbpf, ngencfe, ngencgs, ngencpf, norm, nselbfe, nselbgs, nselbpf, nselcfe, nselcgs, nselcpf, edm::ConfigurableInputSource::numberEventsInRun(), pyint5, and edm::ConfigurableInputSource::run().
00180 { 00181 00182 // do anything here that needs to be done at desctruction time 00183 // (e.g. close files, deallocate resources etc.) 00184 LogVerbatim ("PythiaMuEnrichSource") << "PyhtiaMuEnrichSource: event generation done. " << std::endl; 00185 call_pystat(1); 00186 nbcsim=pyint5.ngen[2][0]-nevnt; 00187 double norm=nbcsim*pyint5.xsec[2][0]; 00188 LogVerbatim ("PythiaMuEnrichSource") << "Simulated "<<nbcsim<<" events for bc production, about to end"; 00189 LogVerbatim ("PythiaMuEnrichSource")<<"Production Nb 2 (bc)"; 00190 LogVerbatim ("PythiaMuEnrichSource")<<"Run Number "<<run(); 00191 LogVerbatim ("PythiaMuEnrichSource")<<"Total Number of Generated bc Events "<<ngenbfe+ngenbpf+ngenbgs+ngencfe+ngencpf+ngencgs; 00192 LogVerbatim ("PythiaMuEnrichSource")<<"Generated B (FE) evts,Xsec prod 2 :"<<ngenbfe<<" "<<((double)ngenbfe)/norm; 00193 LogVerbatim ("PythiaMuEnrichSource")<<"Generated B (PF) evts,Xsec prod 2 :"<<ngenbpf<<" "<<((double)ngenbpf)/norm; 00194 LogVerbatim ("PythiaMuEnrichSource")<<"Generated B (GS) evts,Xsec prod 2 :"<<ngenbgs<<" "<<((double)ngenbgs)/norm; 00195 LogVerbatim ("PythiaMuEnrichSource")<<"Generated B, Total Xsec prod 2 :"<<ngenbpf+ngenbgs+ngenbfe<<" "<<((double)(ngenbpf+ngenbgs+ngenbfe))/norm; 00196 LogVerbatim ("PythiaMuEnrichSource")<<"Generated C (FE) evts,Xsec prod 2 :"<<ngencfe<<" "<<((double)ngenbfe)/norm; 00197 LogVerbatim ("PythiaMuEnrichSource")<<"Generated C (PF) evts,Xsec prod 2 :"<<ngencpf<<" "<<((double)ngenbpf)/norm; 00198 LogVerbatim ("PythiaMuEnrichSource")<<"Generated C (GS) evts,Xsec prod 2 :"<<ngencgs<<" "<<((double)ngenbgs)/norm; 00199 LogVerbatim ("PythiaMuEnrichSource")<<"Generated C, Total Xsec prod 2 :"<<ngencpf+ngencgs+ngencfe<<" "<<((double)(ngencpf+ngencgs+ngencfe))/norm; 00200 LogVerbatim ("PythiaMuEnrichSource")<<"Total Number of Selected bc Events "<<(nselbfe+nselbpf+nselbgs+nselcfe+nselcpf+nselcgs); 00201 LogVerbatim ("PythiaMuEnrichSource")<<"Selected B (FE) evts,Xsec prod 2 :"<<nselbfe<<" "<<((double)nselbfe)/norm; 00202 LogVerbatim ("PythiaMuEnrichSource")<<"Selected B (PF) evts,Xsec prod 2 :"<<nselbpf<<" "<<((double)nselbpf)/norm; 00203 LogVerbatim ("PythiaMuEnrichSource")<<"Selected B (GS) evts,Xsec prod 2 :"<<nselbgs<<" "<<((double)nselbgs)/norm; 00204 LogVerbatim ("PythiaMuEnrichSource")<<"Selected B, Total Xsec prod 2 :"<<nselbpf+nselbgs+nselbfe<<" "<<((double)(nselbpf+nselbgs+nselbfe))/norm; 00205 LogVerbatim ("PythiaMuEnrichSource")<<"Selected C (FE) evts,Xsec prod 2 :"<<nselcfe<<" "<<((double)nselbfe)/norm; 00206 LogVerbatim ("PythiaMuEnrichSource")<<"Selected C (PF) evts,Xsec prod 2 :"<<nselcpf<<" "<<((double)nselbpf)/norm; 00207 LogVerbatim ("PythiaMuEnrichSource")<<"Selected C (GS) evts,Xsec prod 2 :"<<nselcgs<<" "<<((double)nselbgs)/norm; 00208 LogVerbatim ("PythiaMuEnrichSource")<<"Selected C, Total Xsec prod 2 :"<<nselcpf+nselcgs+nselcfe<<" "<<((double)(nselcpf+nselcgs+nselcfe))/norm; 00209 LogVerbatim ("PythiaMuEnrichSource")<<"Total Selected events "<<(numberEventsInRun()+nselbfe+nselbpf+nselbgs+nselcfe+nselcpf+nselcgs); 00210 00211 clear(); 00212 }
| bool edm::PythiaMuEnrichSource::call_pygive | ( | const std::string & | iParm | ) | [private] |
Definition at line 306 of file PythiaMuEnrichSource.cc.
References pydat1, and PYGIVE.
Referenced by PythiaMuEnrichSource().
00307 { 00308 00309 int numWarn = pydat1.mstu[26]; //# warnings 00310 int numErr = pydat1.mstu[22];// # errors 00311 00312 //call the fortran routine pygive with a fortran string 00313 PYGIVE( iParm.c_str(), iParm.length() ); 00314 // PYGIVE( iParm ); 00315 //if an error or warning happens it is problem 00316 return pydat1.mstu[26] == numWarn && pydat1.mstu[22] == numErr; 00317 00318 }
| bool edm::PythiaMuEnrichSource::call_txgive | ( | const std::string & | iParm | ) | [private] |
Definition at line 321 of file PythiaMuEnrichSource.cc.
References lat::endl(), LogTrace, and TXGIVE.
Referenced by PythiaMuEnrichSource().
00321 { 00322 00323 TXGIVE( iParm.c_str(), iParm.length() ); 00324 LogTrace ("PythiaMuEnrichSource") << " " << iParm.c_str() << std::endl; 00325 00326 return 1; 00327 }
| bool edm::PythiaMuEnrichSource::call_txgive_init | ( | ) | [private] |
Definition at line 330 of file PythiaMuEnrichSource.cc.
References lat::endl(), LogTrace, and TXGIVE_INIT.
Referenced by PythiaMuEnrichSource().
00330 { 00331 00332 TXGIVE_INIT(); 00333 LogTrace ("PythiaMuEnrichSource") << " Setting CSA defaults. " << std::endl; 00334 00335 return 1; 00336 }
| void edm::PythiaMuEnrichSource::choose_decay | ( | t_decayed_hadron * | hadron, | |
| unsigned int | nmuon | |||
| ) | [private] |
Definition at line 1022 of file PythiaMuEnrichSource.cc.
References choose_nomuon(), d, t_decay::daughter, t_decayed_hadron::decay, t_decay::first, h, i, t_decayed_hadron::index, K, t_decay::last, LogTrace, MAXMUON, t_decayed_hadron::next, t_decay::next, p, t_decayed_hadron::prob, t_decay::prob, prob_nmuon(), pyr, r, and undecay().
Referenced by wgtmu().
01029 { 01030 unsigned int i; 01031 double p; 01032 double r; 01033 double ph[MAXMUON+1],pn[MAXMUON]; 01034 t_decay *d=0; 01035 t_decayed_hadron *h; 01036 01037 if (!nmuon) 01038 { 01039 choose_nomuon(hadron); 01040 return; 01041 } 01042 01043 if (nmuon>MAXMUON) 01044 { 01045 LogTrace ("PythiaMuEnrichSource") << "error in choose_decay: too many muons requested"; 01046 LogTrace ("PythiaMuEnrichSource") << " - won't get more than "<<MAXMUON; 01047 nmuon=MAXMUON; 01048 } 01049 01050 d=hadron->decay; 01051 if (d==0) return; 01052 if (d->next!=0) // more than one decay 01053 { // => decision needed 01054 p=0; 01055 for (; d!=0; d=d->next) 01056 p+= (nmuon==MAXMUON ? d->prob[nmuon-1] : d->prob[nmuon-1]-d->prob[nmuon]); 01057 r=pyr()*p; 01058 for (d=hadron->decay; d!=0; d=d->next) 01059 { 01060 r-= (nmuon==MAXMUON ? d->prob[nmuon-1] : d->prob[nmuon-1]-d->prob[nmuon]); 01061 if (r<0) 01062 { 01063 r+=1.; // was 2 but 1 should be sufficient 01064 hadron->decay=d; 01065 } 01066 else 01067 undecay (d); 01068 } 01069 } 01070 01071 d=hadron->decay; 01072 d->next=0; // only one decay left 01073 01074 // need to update daughter information here which may be 01075 // wrong because of multiple decays of the same particle: 01076 if ( hadron->index > 0 ) { // don't do it for muons: they have index 0 !!! 01077 K(hadron->index,4)=d->first; 01078 K(hadron->index,5)=d->last; 01079 // note: d->last is the last of all particles arising 01080 // from the decay, not only the immediate daughters 01081 } 01082 h=d->daughter; 01083 if (h==0) 01084 { 01085 LogTrace ("PythiaMuEnrichSource") << "decay without products in cbhoose_decay"; 01086 return; 01087 } 01088 while (h->next!=0) // more than one decay product 01089 { // => decision needed 01090 if (nmuon) 01091 { 01092 prob_nmuon(h->next,nmuon,pn); 01093 if (nmuon==MAXMUON) // for nmuon=MAXMUON, use probabilities for 01094 { // >= MAXMUON, but select only MAXMUON 01095 p=ph[0]=((double)1-h->prob[0])*pn[nmuon-1]; 01096 for (i=1; i<nmuon; ++i) 01097 p+=(ph[i]=(h->prob[i-1]-h->prob[i])*(pn[nmuon-i-1])); 01098 p+=(ph[nmuon]=h->prob[nmuon-1]); 01099 } 01100 else 01101 { 01102 p=ph[0]=((double)1-h->prob[0])*(pn[nmuon-1]-pn[nmuon]); 01103 for (i=1; i<nmuon; ++i) 01104 p+=(ph[i]=(h->prob[i-1]-h->prob[i])*(pn[nmuon-i-1]-pn[nmuon-i])); 01105 p+=(ph[nmuon]=(h->prob[nmuon-1]-h->prob[nmuon])*((double)1-pn[0])); 01106 } 01107 01108 r=pyr()*p; // choose number of muons for this hadron 01109 for (i=nmuon; i; --i) 01110 { 01111 r-=ph[i]; 01112 if (r<0) break; 01113 } 01114 01115 nmuon-=i; 01116 } 01117 else i=0; 01118 choose_decay(h,i); 01119 h=h->next; 01120 } 01121 choose_decay(h,nmuon); 01122 }
Definition at line 657 of file PythiaMuEnrichSource.cc.
References funct::exp(), i, K, LENGTH, funct::log(), LogTrace, P, p, pydat2, pyr, r, radius(), RADIUS, funct::sqrt(), v, and V.
Referenced by decay().
00662 { 00663 00664 int i; 00665 double prob; 00666 00667 double lifetime,mass; 00668 double r,x0,y0; 00669 double px,py,pz,pt,p; 00670 double length,radius,v; 00671 00672 lifetime=pydat2.pmas[3][kc-1]; 00673 r=pyr(); 00674 00675 /* get momentum and determine length of track before intersection with 00676 the bounding cylinder */ 00677 px=P(index,1); 00678 py=P(index,2); 00679 pz=P(index,3); 00680 pt=px*px+py*py; 00681 p=sqrt((double)(pt+pz*pz)); 00682 pt=sqrt((double)(pt)); 00683 00684 length=LENGTH-V(index,3)*((pz>0)-(pz<0)); 00685 v=V(index,1)*V(index,1)+V(index,2)*V(index,2); 00686 if (v>RADIUS*RADIUS || length<0) // particle produced outside detector? 00687 return 0.; 00688 if (pt==0) radius=RADIUS; 00689 else 00690 { 00691 radius=V(index,1)*P(index,1)/pt; 00692 radius+=sqrt((double)radius*radius+RADIUS*RADIUS-v); 00693 } 00694 pz=fabs((double)pz); 00695 if (pt*length>pz*radius) x0=radius*p/pt; 00696 else x0=length*p/pz; 00697 00698 mass=P(index,5); 00699 y0=exp(-mass*x0/(lifetime*p)); 00700 00701 prob=1.-y0; 00702 00703 if (prob<=0) 00704 { 00705 if (prob<0) LogTrace ("PythiaMuEnrichSource") << "choose_decay_vertex: prob<0!"; 00706 return 0; 00707 } 00708 00709 lifetime*=(-log((double)(y0*(1-r)+r))); // choose random proper lifetime 00710 00711 for (i=1; i<=4; ++i) V(index,i)+=lifetime*P(index,i)/mass; 00712 00713 K(index,1)=5; // force decay at a given vertex 00714 00715 return prob; 00716 00717 }
| void edm::PythiaMuEnrichSource::choose_nomuon | ( | t_decayed_hadron * | hadron | ) | [private] |
Definition at line 794 of file PythiaMuEnrichSource.cc.
References d, t_decayed_hadron::decay, t_decayed_hadron::index, K, t_decay::next, and undecay().
Referenced by choose_decay().
00800 { 00801 t_decay *d; 00802 for (d=hadron->decay; d!=0; d=d->next) // delete decay products 00803 undecay(d); 00804 K(hadron->index,1)=1; // declare undecayed 00805 }
Definition at line 628 of file PythiaMuEnrichSource.cc.
References funct::abs(), i, IDMU, K, muon_cut(), and n.
Referenced by pyexec_nomu(), and wgtmu().
00632 { 00633 int i,n; 00634 00635 n=0; 00636 for (i=low; i<=high; ++i) 00637 if (abs(K(i,2))==IDMU && K(i,1)>0 && K(i,1)<=10) 00638 if (muon_cut(i)) ++n; 00639 return n; 00640 }
Definition at line 1154 of file PythiaMuEnrichSource.cc.
References funct::abs(), i, K, N, ngenbfe, ngenbgs, ngenbpf, ngencfe, ngencgs, ngencpf, nselbfe, nselbgs, nselbpf, nselcfe, nselcgs, nselcpf, and ptype.
Referenced by produce().
01156 { 01157 int nbhs=0; 01158 int nchs=0; 01159 int nbgs=0; 01160 int ncgs=0; 01161 if (abs(K(7,2)) == 4) ++nchs; 01162 else if (abs(K(7,2)) == 5 ) ++nbhs; 01163 if (abs(K(8,2)) == 4 ) ++nchs; 01164 else if (abs(K(8,2)) == 5 ) ++nbhs; 01165 for ( int i=1; i<=N; ++i){ 01166 if ( abs(K(K(i,3),2)) == 21 && abs(K(i,2)) == 4 ) ++ncgs; 01167 if ( abs(K(K(i,3),2)) == 21 && abs(K(i,2)) == 5 ) ++nbgs; 01168 } 01169 if ( nbhs == 1 ) { 01170 if (sel==1)++ngenbfe; 01171 else ++nselbfe; 01172 ptype=1; 01173 }else if ( nbhs == 2){ 01174 if (sel==1)ngenbpf++; 01175 else ++nselbpf; 01176 ptype=2; 01177 }else if (nchs==1){ 01178 if (sel==1)ngencfe++; 01179 else ++nselcfe; 01180 ptype=11; 01181 }else if (nchs==2){ 01182 if (sel==1)ngencpf++; 01183 else ++nselcpf; 01184 ptype=12; 01185 } else if (ncgs==2 && (nchs==0 && nbhs==0) ) { 01186 if (sel==1)ngenbgs++; 01187 else ++nselbgs; 01188 ptype=3; 01189 }else if (nbgs==2 && (nchs==0 && nbhs==0) ) { 01190 if (sel==1)ngencgs++; 01191 else ++nselcgs; 01192 ptype=13; 01193 } 01194 }
| t_decayed_hadron * edm::PythiaMuEnrichSource::decay | ( | int | index | ) | [private] |
Definition at line 875 of file PythiaMuEnrichSource.cc.
References a, funct::abs(), brat_muon, brat_nonmuon, choose_decay_vertex(), d1, d2, t_decayed_hadron::decay, hadrons, i, IDMU, t_decayed_hadron::index, K, k, last_decay, loop_hadrons(), MAXMUON, muon_cut(), N, n, t_decayed_hadron::next, t_decay::next, nhadrons, P, t_decayed_hadron::prob, t_decay::prob, ptmin2, pv, pycomp, pydat2, pydat3, pydecy_all(), set_muon_decay(), v, and V.
Referenced by loop_hadrons().
00881 { 00882 int kf,kc; 00883 int i,k,kd,a,n; 00884 t_decay *d1=0; 00885 t_decay *d2=0; 00886 double v[6]; 00887 double pv; 00888 double px,py,pt2; 00889 double prob[MAXMUON]; 00890 00891 kf=K(index,2); 00892 00893 if (abs(kf)==IDMU) // muon 00894 { 00895 if (muon_cut(index)) 00896 { 00897 hadrons[nhadrons].prob[0]=1; 00898 for (i=1; i<MAXMUON; ++i) hadrons[nhadrons].prob[i]=0.; 00899 hadrons[nhadrons].index=index; 00900 hadrons[nhadrons].decay=0; 00901 hadrons[nhadrons].next=0; 00902 return (&(hadrons[nhadrons++])); 00903 } 00904 else 00905 return 0; 00906 } 00907 00908 a=abs(kf); 00909 kc=pycomp(&a); 00910 if (last_decay[kc]<0) return 0; // no possible muon parent 00911 00912 if (ptmin2[kc]>0.) // test whether muon is kinematically possible 00913 { 00914 px=P(index,1); // calculate pt of hadron 00915 py=P(index,2); 00916 pt2=px*px+py*py; 00917 if (pt2<ptmin2[kc]) return 0; 00918 } 00919 00920 for (i=0; i<MAXMUON; ++i) prob[i]=0; 00921 pv=1.; 00922 00923 for (i=1; i<=5; ++i) v[i]=V(index,i); // save production vertex position 00924 if (pydat2.pmas[3][kc-1]>100.) // lifetime long (>10cm ?) 00925 { 00926 pv=choose_decay_vertex(index,kc); 00927 if (pv<=0) 00928 { 00929 K(index,1)=4; // mark particle as "could have decayed but didn't" 00930 00931 return 0; 00932 } 00933 } 00934 00935 pydat3.mdcy[0][kc-1]=1; // switch on decay of this particle 00936 00937 k=kd=K(index,1); // remember status code 00938 00939 if (set_muon_decay (1,kc)) // switch on decays into muons 00940 { 00941 n=N+1; 00942 { 00943 pydecy_all(index); 00944 00945 d2=d1=loop_hadrons (n,N); // decay daughters and count muons 00946 00947 if (d1==0) // no muons? 00948 { 00949 K(index,1)=k; // declare particle as undecayed 00950 N=n-1; // forget all decay products 00951 } 00952 else 00953 for (i=0; i<MAXMUON; ++i) prob[i]=(d1->prob[i]*=brat_muon[kc]); 00954 } 00955 } 00956 if (set_muon_decay (2,kc)) // switch on decays without prompt muons 00957 { 00958 n=N+1; 00959 kd=K(index,1); // remember status code 00960 K(index,1)=k; // declare particle undecayed 00961 pydecy_all(index); // have the particle decay again 00962 00963 d2=loop_hadrons (n,N); // loop over decay products 00964 00965 if (d2==0) // no muons? 00966 { 00967 N=n-1; // => forget about decay 00968 K(index,1)=kd; // restore old status code 00969 d2=d1; 00970 } 00971 else 00972 { 00973 for (i=0; i<MAXMUON; ++i) 00974 prob[i]+=(d2->prob[i]*=brat_nonmuon[kc]); // add probabilities 00975 if (prob[0]>0) d2->next=d1; // => keep decays (d1 may be 0) 00976 00977 } 00978 } 00979 pydat3.mdcy[0][kc-1]=0; // switch decay off again 00980 00981 for (i=1; i<=5; ++i) V(index,i)=v[i]; // restore production vertex position 00982 if (K(index,1)==5) K(index,1)=1; 00983 if (prob[0]>0) // any muon finally? 00984 { 00985 for (i=0; i<MAXMUON; ++i) 00986 hadrons[nhadrons].prob[i]=pv*prob[i]; // take decay probability into account 00987 hadrons[nhadrons].index=index; 00988 hadrons[nhadrons].decay=d2; 00989 hadrons[nhadrons].next=0; 00990 return (&(hadrons[nhadrons++])); 00991 } 00992 else 00993 return (0); 00994 }
| void edm::PythiaMuEnrichSource::default_decay_modes | ( | ) | [private] |
Definition at line 997 of file PythiaMuEnrichSource.cc.
References a, funct::abs(), decay_mode_list, first_muon, i, j, last_decay, LENGTH, pycomp, pydat1, pydat3, and RADIUS.
Referenced by wgtmu().
01001 { 01002 int i,j,a,kc,kf; 01003 01004 i=0; 01005 while ((kf=muon_parents[i])) 01006 { 01007 a=abs(kf); 01008 kc=pycomp(&a); 01009 for (j=first_muon[kc]; j<=last_decay[kc]; ++j) 01010 pydat3.mdme[0][decay_mode_list[j]-1]=1; 01011 pydat3.mdcy[0][kc-1]=1; 01012 ++i; 01013 } 01014 01015 pydat1.mstj[21]=4; // limit decay to cylinder 01016 pydat1.parj[72]=RADIUS; 01017 pydat1.parj[73]=LENGTH; 01018 01019 }
| int edm::PythiaMuEnrichSource::initdc | ( | double * | pt, | |
| double * | eta | |||
| ) | [private] |
Definition at line 388 of file PythiaMuEnrichSource.cc.
References a, funct::abs(), brat_muon, brat_nonmuon, decay_mode_list, etamaxmu, etaminmu, edm::first(), first_muon, i, IDMU, j, edm::es::l(), prof2calltree::last, last_decay, last_muon, last_nonmuon, LogTrace, MASSMU2, par_a, par_b, par_m, ptmin, ptmin2, ptminmu, pycomp, pydat2, pydat3, funct::sqrt(), and switch_off_decays().
Referenced by PythiaMuEnrichSource().
00398 : *minpt - minimum pt required for a muon to be accepted */ 00399 /* *etamax - maximum pseudorapidity for accepted muons */ 00400 /* (muons are accepted in the range -etamax .. etamax)*/ 00401 /* */ 00402 /***********************************************************************/ 00403 { 00404 int i,j,l,first,last,a; 00405 int kf,kc,kcd,mu; 00406 int idecay; 00407 double mass,ptmin,pt2minmu; 00408 00409 for (i=0; i<500; ++i) 00410 first_muon[i]=last_muon[i]=last_nonmuon[i]=last_decay[i]=-1; 00411 00412 i=0; 00413 idecay=0; 00414 kf=muon_parents[i]; 00415 while (kf) 00416 { 00417 a=abs(kf); 00418 kc=pycomp(&a); 00419 first_muon[kc]=idecay; 00420 brat_muon[kc]=0.; 00421 brat_nonmuon[kc]=0.; 00422 kcd=kc; 00423 first=pydat3.mdcy[1][kcd-1]; 00424 if (pydat3.mdme[1][first-1]>=84 && pydat3.mdme[1][first-1]<=88) 00425 { 00426 kcd=pydat3.mdme[1][first-1]; // generic heavy quark decay 00427 first=pydat3.mdcy[1][kcd-1]; 00428 } 00429 last=pydat3.mdcy[1][kcd-1]+pydat3.mdcy[2][kcd-1]; 00430 for (j=first; j<last; ++j) // loop over decay modes 00431 { 00432 mu=0; 00433 for (l=1; l<=5; ++l) // loop over decay products 00434 if (abs(pydat3.kfdp[l-1][j-1])==IDMU || // and count muons 00435 (abs(pydat3.kfdp[l-1][j-1])==443 && // include also J/psi to muon channels 00436 a!=30443)) ++mu; // (but not for psi') 00437 if (mu) // find decay into muons 00438 { 00439 if (pydat3.mdme[0][j-1]==1) // decay switched on? 00440 { 00441 decay_mode_list[idecay++]=j; // store decay mode index 00442 brat_muon[kc]+=pydat3.brat[j-1]; // sum branching ratio 00443 } 00444 } 00445 } 00446 last_muon[kc]=idecay-1; 00447 for (j=first; j<last; ++j) // loop again over decay modes 00448 { 00449 mu=0; 00450 for (l=1; l<=5; ++l) // loop over decay products 00451 if (abs(pydat3.kfdp[l-1][j-1])==IDMU) ++mu; // and count muons 00452 if (!mu) // find decay without muons 00453 { 00454 if (pydat3.mdme[0][j-1]==1) // decay switched on? 00455 { 00456 decay_mode_list[idecay++]=j; // store decay mode index 00457 brat_nonmuon[kc]+=pydat3.brat[j-1]; // sum branching ratio 00458 } 00459 } 00460 } 00461 last_decay[kc]=last_nonmuon[kc]=idecay-1; 00462 00463 /* get pt and eta range for muons to select */ 00464 ptminmu=*minpt; 00465 pt2minmu=ptminmu*ptminmu; 00466 etaminmu=-*etamax; 00467 etamaxmu=*etamax; 00468 00469 /* store some parameters of 2-body decay to allow quick estimate whether 00470 a muon coming out of the decay of the particle can possibly have a 00471 large transverse momentum 00472 */ 00473 mass=pydat2.pmas[0][kc-1]; 00474 par_a[kc]=0.5*(1.+MASSMU2/(mass*mass)); 00475 par_b[kc]=1.-par_a[kc]; 00476 par_m[kc]=mass; 00477 ptmin=ptminmu*(-par_b[kc]*mass*mass/ 00478 (pt2minmu*(sqrt(MASSMU2/pt2minmu+1.)+1.)) +1.); 00479 LogTrace ("PythiaMuEnrichSource") <<"kf="<<kf<<" mass="<<mass<<", ptmin="<<ptmin; 00480 if (ptmin>0.) ptmin2[kc]=ptmin*ptmin; 00481 else ptmin2[kc]=0.; 00482 00483 kf=muon_parents[++i]; // get next hadron 00484 } 00485 00486 /* switch off neutral pion decay */ 00487 a=111; 00488 pydat3.mdcy[0][pycomp(&a)-1]=0; 00489 00490 switch_off_decays(); 00491 00492 return idecay; 00493 }
| bool edm::PythiaMuEnrichSource::isbc | ( | ) | [private] |
Definition at line 1139 of file PythiaMuEnrichSource.cc.
References funct::abs(), i, K, and N.
Referenced by produce().
01141 { 01142 01143 if (abs(K(7,2)) == 4 || abs(K(8,2)) == 4 ) return true; 01144 if (abs(K(7,2)) == 5 || abs(K(8,2)) == 5 ) return true; 01145 int ngcs=0; 01146 int ngbs=0; 01147 for ( int i=1; i<=N; ++i){ 01148 if ( abs(K(K(i,3),2)) == 21 && abs(K(i,2)) == 4 ) ++ngcs; 01149 if ( abs(K(K(i,3),2)) == 21 && abs(K(i,2)) == 5 ) ++ngbs; 01150 } 01151 if ( ngcs >1 || ngbs > 1) return true; 01152 return false; 01153 }
Definition at line 530 of file PythiaMuEnrichSource.cc.
References t_decay::daughter, decay(), decays, t_decay::first, h1, i, K, t_decay::last, MAXMUON, N, n, ndecays, t_decayed_hadron::next, t_decay::next, t_decayed_hadron::prob, and prob_nmuon().
Referenced by decay(), and wgtmu().
00537 { 00538 int i,n; 00539 //int mu; 00540 double pn[MAXMUON]; 00541 t_decayed_hadron *h1; 00542 t_decayed_hadron *previous; 00543 // t_decay d; 00544 t_decay *next; 00545 t_decayed_hadron *daughter; 00546 next=0; 00547 daughter=0; 00548 n=N; 00549 previous=0; 00550 //LogTrace ("PythiaMuEnrichSource") << "in Loop_hadrons first="<<first<<" last="<<last; 00551 for (i=first; i<=last; ++i) 00552 { 00553 if (K(i,1)>0 && K(i,1)<=10) // undecayed particle 00554 { 00555 00556 h1=decay(i); 00557 if (h1!=0) 00558 { 00559 pn[0]=h1->prob[0]; 00560 if (pn[0]>0) // any chance to get a muon? 00561 { 00562 if (previous!=0) previous->next=h1; 00563 else daughter=h1; 00564 previous=h1; 00565 } 00566 } 00567 } 00568 } 00569 if (daughter!=0) prob_nmuon(daughter,MAXMUON-1,pn); 00570 else return 0; 00571 00572 if (pn[0]>0) 00573 { 00574 for (i=0; i<MAXMUON; ++i) decays[ndecays].prob[i]=pn[i]; 00575 decays[ndecays].first=first; 00576 decays[ndecays].last=N; 00577 decays[ndecays].daughter=daughter; 00578 decays[ndecays].next=next; 00579 return (&(decays[ndecays++])); 00580 } 00581 else 00582 return 0; 00583 }
Definition at line 586 of file PythiaMuEnrichSource.cc.
References eta, etamaxmu, etaminmu, funct::log(), nmu_min, P, p, ptminmu, funct::sqrt(), and wtype.
Referenced by count_all_muons(), count_muons(), decay(), and pydecy_mu().
00591 { 00592 double px,py,pz,pt,pt2,eta,p; 00593 int iselp,iselm; 00594 iselp=wtype; 00595 iselm=nmu_min; 00596 px=P(index,1); 00597 py=P(index,2); 00598 pz=P(index,3); 00599 pt2=px*px+py*py; 00600 p=sqrt((double)(pt2+pz*pz)); 00601 pt=sqrt((double)(pt2)); 00602 00603 eta=log((double)((p+fabs((double)pz))/pt))*((pz>0)-(pz<0)); 00604 00605 if(iselp == 1 && iselm==1) 00606 00607 // special for single mu, pt1 (wtype=1) production (U.G.) : 00608 { 00609 00610 if( pt > 4.0 ) return 0; 00611 if (eta<etaminmu || eta>etamaxmu) return 0; 00612 if (eta>etaminmu && eta<-1.7 && p<3.5) return 0; 00613 if (eta>-1.7 && eta<-1.2 && pt<1.8) return 0; 00614 if (eta>-1.2 && eta< 1.2 && pt<3.0) return 0; 00615 if (eta> 1.2 && eta< 1.7 && pt<1.8) return 0; 00616 if (eta> 1.7 && eta<etamaxmu && p<3.5) return 0; 00617 } 00618 00619 // for all other productions : 00620 else 00621 { 00622 if (eta<etaminmu || eta>etamaxmu || pt<ptminmu) return 0; 00623 } 00624 return 1; 00625 }
| void edm::PythiaMuEnrichSource::prob_nmuon | ( | t_decayed_hadron * | hadron, | |
| int | nmuon, | |||
| double * | pn | |||
| ) | [private] |
Definition at line 496 of file PythiaMuEnrichSource.cc.
References h, i, j, t_decayed_hadron::next, and t_decayed_hadron::prob.
Referenced by choose_decay(), and loop_hadrons().
00502 : hadron - first hadron in the list, hadron->next is expected */ 00503 /* to point to the next one, the last has next=0 */ 00504 /* hadron->prob[n] is the probability to get >n muons */ 00505 /* from this hadron */ 00506 /* nmuon - maximum number of muons to calculate */ 00507 /* output: pn[] - array of probabilities */ 00508 /* p[0] is the probability to get >0 muons */ 00509 /* p[1] . . . . . >1 muons */ 00510 /* .. . . . . . .. .... */ 00511 /* p[nmuon] . . . . >nmuon muons */ 00512 /* */ 00513 /***********************************************************************/ 00514 { 00515 int i,j; 00516 t_decayed_hadron *h; 00517 00518 for (i=0; i<=nmuon; ++i) pn[i]=0; 00519 00520 for (h=hadron; h!=0; h=h->next) 00521 for (i=nmuon; i>=0; --i) 00522 { 00523 for (j=i-1; j>=0; --j) 00524 pn[i]+=h->prob[i-j-1]*(pn[j]-pn[j+1]); 00525 pn[i]+=h->prob[i]*((double)1-pn[0]); 00526 } 00527 }
Implements edm::ConfigurableInputSource.
Definition at line 222 of file PythiaMuEnrichSource.cc.
References call_pylist(), conv, countbc(), lat::endl(), edm::ConfigurableInputSource::event(), evt, firstevent, firstRun, int, isbc(), LogTrace, maxEventsToPrint_, nbc, nbcsim, nevnt, nmu_min, nrejbc, edm::ConfigurableInputSource::numberEventsInRun(), edm::Event::put(), pyint5, pypars, pythiaHepMCVerbosity_, pythiaPylistVerbosity_, edm::InputSource::remainingEvents(), edm::ConfigurableInputSource::run(), wbc, wgtmu(), and wmin.
00223 { 00224 std::auto_ptr<HepMCProduct> bare_product(new HepMCProduct()); 00225 LogVerbatim ("PythiaMuEnrichSource") <<"=============================================================="; 00226 LogTrace ("PythiaMuEnrichSource") << "PyhtiaMuEnrichSource: Generating event ... " << std::endl; 00227 double wg=0.; 00228 if (firstevent){ 00229 firstRun=run(); 00230 firstevent=false; 00231 } 00232 if ( (int) run() !=firstRun && nbc==0 ) { 00233 call_pystat(1); 00234 nevnt=pyint5.ngen[2][0]; 00235 nbc=(int) (((double)nevnt)/wbc); 00236 LogVerbatim ("PythiaMuEnrichSource")<<"Number of events to generate for b-c events:"<<nbc; 00237 LogVerbatim ("PythiaMuEnrichSource")<<"Number of rejected events in prod ! because of bc:"<<nrejbc; 00238 LogVerbatim ("PythiaMuEnrichSource")<<"Production Nb1 (without bc)"; 00239 LogVerbatim ("PythiaMuEnrichSource")<<"Run Number "<<firstRun; 00240 LogVerbatim ("PythiaMuEnrichSource")<<"Total Number of Generated Events"<<nevnt; 00241 LogVerbatim ("PythiaMuEnrichSource")<<"Cross Section (mb) "<<pyint5.xsec[2][0]; 00242 LogVerbatim ("PythiaMuEnrichSource")<<"Selected events "<<numberEventsInRun(); 00243 double integlumi=pyint5.ngen[2][0]/(1.E+6*pyint5.xsec[2][0]); 00244 LogVerbatim ("PythiaMuEnrichSource")<<"Integrated LHC lumi (nb-1) "<<integlumi; 00245 edm::Service<JobReport> reportSvc; 00246 std::stringstream stmp,srun; 00247 srun << firstRun; 00248 stmp << integlumi/nanobarn; 00249 reportSvc->reportGeneratorInfo("runX",srun.str()); 00250 reportSvc->reportGeneratorInfo("lumiX",stmp.str()); 00251 //pypars.mstp[81]=1; // for speed of production 00252 } 00253 bool bcevent=false; 00254 while ( wg <wmin && nbcsim <= nbc ){ 00255 call_pyevnt(); // generate one event with Pythia 00256 bcevent=isbc(); 00257 if(bcevent) countbc(1); 00258 if( (int) run()!= firstRun) nbcsim++; 00259 if ( (int) run() == firstRun && bcevent ){ 00260 wg=0.; 00261 nrejbc++; 00262 }else if ( (int) run() == firstRun || bcevent )wg=wgtmu(nmu_min,wmin); 00263 else wg=0.; 00264 } 00265 00266 if (wg >= wmin ){ 00267 if (bcevent)countbc(2); 00268 call_pyhepc(1); 00269 00270 00271 //HepMC::GenEvent* evt = conv.getGenEventfromHEPEVT(); 00272 HepMC::GenEvent* evt = conv.read_next_event(); 00273 evt->set_signal_process_id(pypars.msti[0]); 00274 evt->set_event_scale(pypars.pari[16]); 00275 evt->set_event_number(numberEventsInRun() - remainingEvents() - 1); 00276 evt->weights().push_back(wg); 00277 //******** Verbosity ******** 00278 00279 if(event() <= maxEventsToPrint_ && 00280 (pythiaPylistVerbosity_ || pythiaHepMCVerbosity_)) { 00281 00282 // Prints PYLIST info 00283 if(pythiaPylistVerbosity_) { 00284 call_pylist(pythiaPylistVerbosity_); 00285 } 00286 00287 // Prints HepMC event 00288 if(pythiaHepMCVerbosity_) { 00289 LogTrace ("PythiaMuEnrichSource") << "Event process = " << pypars.msti[0] << std::endl 00290 << "----------------------" << std::endl; 00291 evt->print(); 00292 } 00293 } 00294 00295 00296 if(evt) bare_product->addHepMCData(evt ); 00297 00298 iEvent.put(bare_product); 00299 } 00300 if ( nbcsim >= nbc && nbcsim > 1 )return false; 00301 else return true; 00302 00303 00304 }
Definition at line 720 of file PythiaMuEnrichSource.cc.
References funct::abs(), i, K, prof2calltree::last, LogTrace, N, n, pycomp, pydat3, and pydecy.
Referenced by decay(), pydecy_mu(), and pyexec_nomu().
00723 { 00724 int i,n,kf,last; 00725 00726 n=N+1; 00727 pydecy(&index); 00728 last=N; 00729 00730 for (i=n; i<=last; ++i) 00731 { 00732 if (K(i,1)>0 && K(i,1)<=10) 00733 { 00734 kf=abs(K(i,2)); 00735 if ( K(i,1) != 1 ) LogTrace ("PythiaMuEnrichSource") << " in loop pydecy_all i="<<i<<" kf="<<kf<<" status="<<K(i,1)<<" daughters="<<K(i,4)<<"-"<<K(i,5); 00736 if (pydat3.mdcy[0][pycomp(&kf)-1]) pydecy_all(i); 00737 } 00738 } 00739 }
Definition at line 808 of file PythiaMuEnrichSource.cc.
References funct::abs(), i, IDMU, j, K, k, muon_cut(), N, n, pycomp, pydat3, pydecy, and pydecy_all().
00814 { 00815 int i,j,found,k,kd,kf,n,nsuccess; 00816 00817 n=N+1; 00818 00819 nsuccess=0; 00820 k=kd=K(index,1); // remember status code of undecayed particle 00821 for (i=0; i<ndecay; ++i) 00822 { 00823 pydecy(&index); // decay the particle 00824 found=0; 00825 for (j=n; j<=N; ++j) // loop over decay products 00826 if (abs(K(j,2))==IDMU) // find muons 00827 if (muon_cut(j)) 00828 { 00829 found++; 00830 break; 00831 } 00832 if (found) 00833 { 00834 ++nsuccess; // in case of success, 00835 for (j=n; j<=N; ++j) // make sure all daughters are decayed 00836 if (K(j,1)>0 && K(j,1)<=10) 00837 { 00838 kf=abs(K(j,2)); 00839 if (pydat3.mdcy[0][pycomp(&kf)-1]) pydecy_all(j); 00840 } 00841 n=N+1; // keep decay products 00842 kd=K(index,1); // remember decayed status code 00843 break; 00844 } 00845 else // no muon found 00846 { 00847 K(index,1)=k; // forget about decay 00848 N=n-1; 00849 } 00850 } 00851 00852 for (++i; i<ndecay; ++i) // continue the loop 00853 { 00854 pydecy(&index); // decay the particle again 00855 found=0; 00856 for (j=n; j<=N; ++j) // loop over decay products 00857 if (abs(K(j,2))==IDMU) // find muons 00858 if (muon_cut(j)) 00859 { 00860 found++; 00861 break; 00862 } 00863 if (found) 00864 ++nsuccess; // count successes, 00865 K(index,1)=k; // but forget about the decay anyway 00866 N=n-1; 00867 } 00868 00869 K(index,1)=kd; // update status code (possibly decayed) 00870 00871 return nsuccess; 00872 }
| void edm::PythiaMuEnrichSource::pyexec_nomu | ( | ) | [private] |
Definition at line 753 of file PythiaMuEnrichSource.cc.
References funct::abs(), count_muons(), i, K, prof2calltree::last, LogTrace, MAXTRYNOMU, N, n, pycomp, pydat3, and pydecy_all().
Referenced by wgtmu().
00760 { 00761 #define MAXTRYNOMU 50 00762 int i,kf,last,n,tried; 00763 00764 last=N; 00765 00766 for (i=1; i<=last; ++i) 00767 { 00768 kf=abs(K(i,2)); 00769 if (K(i,1)==1) 00770 { 00771 if (pydat3.mdcy[0][pycomp(&kf)-1]) 00772 { 00773 n=N+1; 00774 pydecy_all(i); // decay 00775 tried=0; 00776 while (count_muons(n,N)) 00777 { 00778 if (++tried>=MAXTRYNOMU) 00779 { 00780 LogTrace ("edm::PythiaMuEnrichSource") << "error in pyexec_nomu: unable to get decay without muons after "<<tried<<" attempts"; 00781 N=n-1; 00782 K(i,1)=1; 00783 break; 00784 } 00785 N=n-1; 00786 K(i,1)=1; 00787 pydecy_all(i); 00788 } 00789 } 00790 } 00791 } 00792 }
Definition at line 339 of file PythiaMuEnrichSource.cc.
References decay_mode_list, first_muon, i, j, k, prof2calltree::last, last_decay, last_muon, and pydat3.
Referenced by decay().
00344 : 1 if selected decay has nonzero branching ratio 00345 0 otherwise 00346 */ 00347 { 00348 int i,j,k,last; 00349 00350 k=0; 00351 j=(mode!=2); 00352 last=last_muon[kc]; 00353 for (i=first_muon[kc]; i<=last; ++i) 00354 { 00355 pydat3.mdme[0][decay_mode_list[i]-1]=j; 00356 k|=j; 00357 } 00358 last=last_decay[kc]; 00359 j=(mode!=1); 00360 for (; i<=last; ++i) 00361 { 00362 pydat3.mdme[0][decay_mode_list[i]-1]=j; 00363 k|=j; 00364 } 00365 return(k); 00366 }
| void edm::PythiaMuEnrichSource::switch_off_decays | ( | ) | [private] |
Definition at line 368 of file PythiaMuEnrichSource.cc.
References a, funct::abs(), i, pycomp, pydat1, and pydat3.
Referenced by initdc(), and wgtmu().
00370 { 00371 int i,a,kc,kf; 00372 00373 i=0; 00374 while ((kf=muon_parents[i])) 00375 { 00376 a=abs(kf); 00377 kc=pycomp(&a); 00378 pydat3.mdcy[0][kc-1]=0; 00379 ++i; 00380 } 00381 00382 pydat1.mstj[21]=1; 00383 00384 }
Definition at line 742 of file PythiaMuEnrichSource.cc.
References t_decay::first, i, K, t_decay::last, and prof2calltree::last.
Referenced by choose_decay(), and choose_nomuon().
00746 { 00747 int i,last; 00748 last=decay->last; 00749 for (i=decay->first; i<=last; ++i) K(i,1)=0; 00750 }
| double edm::PythiaMuEnrichSource::wgtmu | ( | int | nmu, | |
| float | minweight | |||
| ) | [private] |
Definition at line 1195 of file PythiaMuEnrichSource.cc.
References funct::abs(), choose_decay(), count_all_muons(), count_muons(), d, t_decayed_hadron::decay, default_decay_modes(), h, i, t_decayed_hadron::index, K, LogTrace, loop_hadrons(), MAXMUON, N, nbcsim, ndecays, t_decayed_hadron::next, nhadrons, t_decayed_hadron::prob, t_decay::prob, pyedit, pyexec, pyexec_nomu(), pyjets, pyr, r, switch_off_decays(), wbc, and weight.
Referenced by produce().
01197 : nmuon minimum number of desired muons 01198 fl minimum desired weight 01199 returned: weight of the event (discard event if wgtmu_=0.) 01200 */ 01201 { 01202 // bc_separately declaration 01203 01204 int idpy,ndpy; 01205 int hard_b,hard_c; 01206 int ngs_b,ngs_c; 01207 int i; 01208 t_decay *d=0; 01209 t_decayed_hadron h; 01210 double r; 01211 double weight,wtxs; 01212 01213 weight=1; 01214 01215 if ((i=count_muons(1,N))) 01216 { 01217 LogTrace ("PythiaMuEnrichSource") << i<< " muons already found before decay\n"; 01218 } 01219 01220 nmuon-=i; 01221 nhadrons=ndecays=0; 01222 if (nmuon<=0) // already enough muons? 01223 { 01224 default_decay_modes(); 01225 pyexec(); // decay the event without constraints 01226 switch_off_decays(); // prepare for next event 01227 } 01228 else 01229 { 01230 d=loop_hadrons (9,N); // decay all hadrons into muons 01231 if (d==0) return 0; // no more muons? 01232 01233 weight=d->prob[nmuon-1]; 01234 if (weight==0) return 0; 01235 01236 r=pyr()*d->prob[nmuon-1]; // choose number of muons to select 01237 01238 for (i=nmuon; i<MAXMUON; ++i) 01239 if (r>d->prob[i]) break; 01240 01241 if (i==MAXMUON) 01242 { 01243 LogTrace ("PythiaMuEnrichSource") << "Warning in WGTMU: Event wants "<<i<<" or more muons,"; 01244 LogTrace ("PythiaMuEnrichSource") << " but will never get more than "<<i; 01245 LogTrace ("PythiaMuEnrichSource") << " If needed, this can be fixed by increasing the parameter MAXMUON"; 01246 } 01247 01248 nmuon=i; 01249 } 01250 01251 01252 // this is the pyevwt of Nicola 01253 /* get PYTHIA weight */ 01254 01255 //double pt2=pyint1.vint[47]; 01256 //double pthat=sqrt(pt2); 01257 //double w1=1; 01258 //if ( wtype == 3 ) w1=0.2*(1+0.0003*pthat*pthat); 01259 //wtxs=w1<1?w1:1; 01260 //wtxs=weight/wtxs; 01261 wtxs=weight; 01262 01263 // ============================================================================ 01264 // identify leading order c-cbar and b-bbar production and flavour excitation 01265 // and reduce the weight of these events by factor wbc 01266 // These events have to generated separately and added in order to recover 01267 // a complete sample! 01268 // b-events are additionally scaled by 500mb/bbxsec(pthat=0,mstp(82)=1)=0.42 -> not done 01269 // following recommendations for LHC b-physics 01270 // ============================================================================ 01271 /* down-scale gluon splitting (GS) events and b,c events (FE + PF) */ 01272 01273 // 01274 ngs_b =0; 01275 ngs_c =0; 01276 hard_b=0; 01277 hard_c=0; 01278 ndpy=pyjets.n; 01279 01280 if (abs(K(7,2))==5 || abs(K(8,2))==5)hard_b=1; 01281 if (abs(K(7,2))==4 || abs(K(8,2))==4)hard_c=1; 01282 01283 for (idpy=1; idpy<=ndpy; ++idpy) 01284 { 01285 if ( abs(K(K(idpy,3),2))==21 && abs(K(idpy,2))==4 ) ++ngs_c; 01286 if ( abs(K(K(idpy,3),2))==21 && abs(K(idpy,2))==5 ) ++ngs_b; 01287 } 01288 //rescale by 10 (wbc) to account for enriched luminosity 01289 01290 if( ngs_b >= 2 || hard_b==1 || ngs_c>=2 || hard_c==1 )wtxs*=wbc; 01291 01292 //additional rescaling for bb xsec 01293 01294 //if( ngs_b >= 2 || hard_b==1 ) wtxs*=0.42; 01295 01296 // ============================================================================ 01297 01298 /* downscale events with low weights */ 01299 if (wtxs>0 && wtxs<fl) 01300 { 01301 r=pyr(); 01302 if (r<wtxs/fl) 01303 { 01304 wtxs=fl; 01305 } 01306 else 01307 { 01308 return 0.; 01309 } 01310 } 01311 else 01312 LogTrace ("PythiaMuEnrichSource") << "Event not downscaled, weight="<<wtxs<<" Number of muons to produce by decay:"<<nmuon; 01313 01314 if (nmuon>0) 01315 { 01316 // create dummy hadron in order to avoid duplicating the loop over hadrons 01317 // the last decay pointer contains all the probabilities although no real decay 01318 h.index=0; 01319 for (i=0; i<MAXMUON; ++i) h.prob[i]=d->prob[i]; 01320 h.decay=d; 01321 h.next=0; 01322 default_decay_modes(); 01323 choose_decay(&h,nmuon); // select a decay mode with nmu muons 01324 switch_off_decays(); // prepare for next event 01325 i=12; 01326 pyedit(&i); // remove deleted particles from event record 01327 LogTrace ("PythiaMuEnrichSource") << "WGTMU: event weight="<<d->prob[0]<<", after downscaling:"<<wtxs; 01328 } 01329 default_decay_modes(); 01330 pyexec_nomu(); // decay all hadrons without any muons 01331 switch_off_decays(); // prepare for next event 01332 LogTrace ("PythiaMuEnrichSource") << "After All N="<<N<<" and number of muons=" << count_all_muons(9,N) <<" nbcsim="<<nbcsim; 01333 return wtxs; 01334 }
| void edm::PythiaMuEnrichSource::zero_hepevnt | ( | ) | [private] |
Definition at line 1124 of file PythiaMuEnrichSource.cc.
References i.
01125 { 01126 for ( int i=0; i < HepMC::HEPEVT_Wrapper::number_entries(); ++i ) { 01127 HepMC::HEPEVT_Wrapper::set_status( i, 0 ); 01128 HepMC::HEPEVT_Wrapper::set_id( i, 0 ); 01129 HepMC::HEPEVT_Wrapper::set_parents( i, 0, 0 ); 01130 HepMC::HEPEVT_Wrapper::set_children( i, 0, 0 ); 01131 HepMC::HEPEVT_Wrapper::set_momentum( i, 0, 0, 0, 0 ); 01132 HepMC::HEPEVT_Wrapper::set_mass( i, 0 ); 01133 HepMC::HEPEVT_Wrapper::set_position( i, 0, 0, 0, 0 ); 01134 } 01135 HepMC::HEPEVT_Wrapper::set_number_entries( 0 ); 01136 01137 }
double edm::PythiaMuEnrichSource::brat_muon[500] [private] |
double edm::PythiaMuEnrichSource::brat_nonmuon[500] [private] |
double edm::PythiaMuEnrichSource::comenergy [private] |
int edm::PythiaMuEnrichSource::decay_mode_list[2000] [private] |
Definition at line 116 of file PythiaMuEnrichSource.h.
Referenced by default_decay_modes(), initdc(), and set_muon_decay().
t_decay edm::PythiaMuEnrichSource::decays[2000] [private] |
double edm::PythiaMuEnrichSource::etamaxmu [private] |
double edm::PythiaMuEnrichSource::etaminmu [private] |
HepMC::GenEvent* edm::PythiaMuEnrichSource::evt [private] |
int edm::PythiaMuEnrichSource::first_muon[500] [private] |
Definition at line 118 of file PythiaMuEnrichSource.h.
Referenced by default_decay_modes(), initdc(), and set_muon_decay().
bool edm::PythiaMuEnrichSource::firstevent [private] |
Definition at line 105 of file PythiaMuEnrichSource.h.
Referenced by produce(), and PythiaMuEnrichSource().
int edm::PythiaMuEnrichSource::firstRun [private] |
Definition at line 103 of file PythiaMuEnrichSource.h.
Referenced by produce(), and PythiaMuEnrichSource().
t_decayed_hadron edm::PythiaMuEnrichSource::hadrons[4000] [private] |
int edm::PythiaMuEnrichSource::last_decay[500] [private] |
Definition at line 121 of file PythiaMuEnrichSource.h.
Referenced by decay(), default_decay_modes(), initdc(), and set_muon_decay().
int edm::PythiaMuEnrichSource::last_muon[500] [private] |
Definition at line 119 of file PythiaMuEnrichSource.h.
Referenced by initdc(), and set_muon_decay().
int edm::PythiaMuEnrichSource::last_nonmuon[500] [private] |
double edm::PythiaMuEnrichSource::lmu_ptmin [private] |
double edm::PythiaMuEnrichSource::maxeta [private] |
unsigned int edm::PythiaMuEnrichSource::maxEventsToPrint_ [private] |
Events to print if verbosity.
Definition at line 95 of file PythiaMuEnrichSource.h.
Referenced by produce(), and PythiaMuEnrichSource().
int edm::PythiaMuEnrichSource::n_decay_mode [private] |
Definition at line 117 of file PythiaMuEnrichSource.h.
int edm::PythiaMuEnrichSource::nbc [private] |
Definition at line 104 of file PythiaMuEnrichSource.h.
Referenced by produce(), and PythiaMuEnrichSource().
int edm::PythiaMuEnrichSource::nbcsim [private] |
Definition at line 104 of file PythiaMuEnrichSource.h.
Referenced by produce(), PythiaMuEnrichSource(), wgtmu(), and ~PythiaMuEnrichSource().
int edm::PythiaMuEnrichSource::ndecays [private] |
int edm::PythiaMuEnrichSource::nevnt [private] |
Definition at line 100 of file PythiaMuEnrichSource.h.
Referenced by produce(), PythiaMuEnrichSource(), and ~PythiaMuEnrichSource().
int edm::PythiaMuEnrichSource::ngenbfe [private] |
Definition at line 128 of file PythiaMuEnrichSource.h.
Referenced by countbc(), PythiaMuEnrichSource(), and ~PythiaMuEnrichSource().
int edm::PythiaMuEnrichSource::ngenbgs [private] |
Definition at line 128 of file PythiaMuEnrichSource.h.
Referenced by countbc(), PythiaMuEnrichSource(), and ~PythiaMuEnrichSource().
int edm::PythiaMuEnrichSource::ngenbpf [private] |
Definition at line 128 of file PythiaMuEnrichSource.h.
Referenced by countbc(), PythiaMuEnrichSource(), and ~PythiaMuEnrichSource().
int edm::PythiaMuEnrichSource::ngencfe [private] |
Definition at line 128 of file PythiaMuEnrichSource.h.
Referenced by countbc(), PythiaMuEnrichSource(), and ~PythiaMuEnrichSource().
int edm::PythiaMuEnrichSource::ngencgs [private] |
Definition at line 128 of file PythiaMuEnrichSource.h.
Referenced by countbc(), PythiaMuEnrichSource(), and ~PythiaMuEnrichSource().
int edm::PythiaMuEnrichSource::ngencpf [private] |
Definition at line 128 of file PythiaMuEnrichSource.h.
Referenced by countbc(), PythiaMuEnrichSource(), and ~PythiaMuEnrichSource().
int edm::PythiaMuEnrichSource::nhadrons [private] |
int edm::PythiaMuEnrichSource::nmu_min [private] |
int edm::PythiaMuEnrichSource::nrejbc [private] |
Definition at line 130 of file PythiaMuEnrichSource.h.
Referenced by produce(), and PythiaMuEnrichSource().
int edm::PythiaMuEnrichSource::nselbfe [private] |
Definition at line 129 of file PythiaMuEnrichSource.h.
Referenced by countbc(), PythiaMuEnrichSource(), and ~PythiaMuEnrichSource().
int edm::PythiaMuEnrichSource::nselbgs [private] |
Definition at line 129 of file PythiaMuEnrichSource.h.
Referenced by countbc(), PythiaMuEnrichSource(), and ~PythiaMuEnrichSource().
int edm::PythiaMuEnrichSource::nselbpf [private] |
Definition at line 129 of file PythiaMuEnrichSource.h.
Referenced by countbc(), PythiaMuEnrichSource(), and ~PythiaMuEnrichSource().
int edm::PythiaMuEnrichSource::nselcfe [private] |
Definition at line 129 of file PythiaMuEnrichSource.h.
Referenced by countbc(), PythiaMuEnrichSource(), and ~PythiaMuEnrichSource().
int edm::PythiaMuEnrichSource::nselcgs [private] |
Definition at line 129 of file PythiaMuEnrichSource.h.
Referenced by countbc(), PythiaMuEnrichSource(), and ~PythiaMuEnrichSource().
int edm::PythiaMuEnrichSource::nselcpf [private] |
Definition at line 129 of file PythiaMuEnrichSource.h.
Referenced by countbc(), PythiaMuEnrichSource(), and ~PythiaMuEnrichSource().
double edm::PythiaMuEnrichSource::par_a[500] [private] |
double edm::PythiaMuEnrichSource::par_b[500] [private] |
double edm::PythiaMuEnrichSource::par_m[500] [private] |
double edm::PythiaMuEnrichSource::ptmin2[500] [private] |
double edm::PythiaMuEnrichSource::ptminmu [private] |
int edm::PythiaMuEnrichSource::ptype [private] |
HepMC verbosity flag.
Definition at line 93 of file PythiaMuEnrichSource.h.
Referenced by produce(), and PythiaMuEnrichSource().
unsigned int edm::PythiaMuEnrichSource::pythiaPylistVerbosity_ [private] |
Pythia PYLIST Verbosity flag.
Definition at line 91 of file PythiaMuEnrichSource.h.
Referenced by produce(), and PythiaMuEnrichSource().
double edm::PythiaMuEnrichSource::sumwg [private] |
double edm::PythiaMuEnrichSource::sumwg2 [private] |
double edm::PythiaMuEnrichSource::wbc [private] |
double edm::PythiaMuEnrichSource::wmin [private] |
double edm::PythiaMuEnrichSource::wprintmin [private] |
Definition at line 99 of file PythiaMuEnrichSource.h.
int edm::PythiaMuEnrichSource::wtype [private] |
1.5.4