CGEM BOSS 6.6.5.i
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDecay Class Reference

#include <EvtDecay.h>

+ Inheritance diagram for EvtDecay:

Public Member Functions

 EvtDecay (const string &name, ISvcLocator *pSvcLocator)
 
StatusCode initialize ()
 
StatusCode execute ()
 
StatusCode finalize ()
 

Public Attributes

IDataInfoSvctmpInfoSvc
 
DataInfoSvcdataInfoSvc
 

Detailed Description

Definition at line 64 of file EvtDecay.h.

Constructor & Destructor Documentation

◆ EvtDecay()

EvtDecay::EvtDecay ( const string & name,
ISvcLocator * pSvcLocator )

Definition at line 75 of file EvtDecay.cxx.

75 :Algorithm( name, pSvcLocator )
76{
77
78
79
80 // these can be used to specify alternative locations or filenames
81 // for the EvtGen particle and channel definition files.
82
83 declareProperty("DecayDecDir", m_DecayDec = "");
84 declareProperty("PdtTableDir", m_PdtTable = "");
85 declareProperty("userDecayTableName", userDecFileName = "NONE");
86 declareProperty("DecayTopology", m_DecayTop = ""); // output decay topology to a file specified by user
87 declareProperty("DecayRecord", m_DecayRec = ""); // output decay topology of one particle specified by user to a file
88 declareProperty("ParentParticle", m_ParentPart = "NONE");// Mothor particle name in pdt.table for only running BesEvtGen
89 declareProperty("Multiplicity", m_Ncharge = false); // output ncharge number of an event
90 declareProperty("NutupleFile",m_NtupleFile = false); // output Ntuple file
91 declareProperty("mDIY",_mDIY= false); // mDIY model
92 declareProperty("FDPdata",m_SB3File = ""); // Fit Dalitz Plot (FDP) data file name (*.root)
93 declareProperty("FDPHisT",m_SB3HT = ""); // histogram title of Dalitz plot in data file
94 declareProperty("FDPpart",m_FDPparticle =""); //to assign the FDP parent particle name (string)
95 declareProperty("TruthFile",m_truthFile ="");
96 declareProperty("TruthPart",m_truthPart ="");
97 declareProperty("Psi3SopenCharm",m_Psi4040OpenCharm=false);
98 declareProperty("Psi2openCharm", m_Psi2openCharm=false);
99 declareProperty("statDecays",m_statDecays=false);
100 declareProperty("TagLundCharmModel", m_tagLundModel=false);
101 //for polarized charmonium production
102 m_polarization.clear();
103 declareProperty("polarization", m_polarization);
104}

Member Function Documentation

◆ execute()

StatusCode EvtDecay::execute ( )

Definition at line 265 of file EvtDecay.cxx.

266{
267
268 MsgStream log(messageService(), name());
269 // log << MSG::INFO << "EvtDecay executing" << endreq;
270 //------------------
271 string key = "GEN_EVENT";
272 // retrieve event from Transient Store (Storegate)
273 SmartDataPtr<McGenEventCol> McEvtColl(eventSvc(), "/Event/Gen");
274
275 m_numberEvent += 1;
276 writeFlag = 0;
277 //std::cout<<"EvtNumber= "<<m_numberEvent<<std::endl;
278 if ( McEvtColl == 0 )
279 {
280 HepMC::GenEvent* evt=new GenEvent();
281 evt->set_event_number(m_numberEvent);
282 callBesEvtGen( evt );
283 if(!(m_DecayTop=="")) evt->print(outfile);
284
285 //create a Transient store
286 McGenEventCol *mcColl = new McGenEventCol;
287 McGenEvent* mcEvent = new McGenEvent(evt);
288 mcColl->push_back(mcEvent);
289 StatusCode sc = eventSvc()->registerObject("/Event/Gen",mcColl);
290 if(sc != SUCCESS) return StatusCode::FAILURE;
291
292 } else // (McEvtColl != 0)
293 {
294 McGenEventCol::iterator mcItr;
295 for( mcItr = McEvtColl->begin(); mcItr != McEvtColl->end(); mcItr++ )
296 {
297 HepMC::GenEvent* hEvt = (*mcItr)->getGenEvt();
298 // MeVToGeV( hEvt );
299
300 callEvtGen( hEvt );
301
302 if(!(m_DecayTop=="")) hEvt->print(outfile);
303 // GeVToMeV( hEvt );
304 // if(!(m_DecayRec=="")) outfile2<<std::endl;
305 if(!(m_DecayRec=="")) std::cout<<std::endl;
306 };
307 }
308 if( m_NtupleFile ){ m_tuple->write();}
309 return StatusCode::SUCCESS;
310}
ObjectVector< McGenEvent > McGenEventCol
Definition McGenEvent.h:39
*************DOUBLE PRECISION m_pi *DOUBLE PRECISION m_HvecTau2 DOUBLE PRECISION m_HvClone2 DOUBLE PRECISION m_gamma1 DOUBLE PRECISION m_gamma2 DOUBLE PRECISION m_thet1 DOUBLE PRECISION m_thet2 INTEGER m_IFPHOT *COMMON c_Taupair $ !Spin Polarimeter vector first Tau $ !Spin Polarimeter vector second Tau $ !Clone Spin Polarimeter vector first Tau $ !Clone Spin Polarimeter vector second Tau $ !Random Euler angle for cloning st tau $ !Random Euler angle for cloning st tau $ !Random Euler angle for cloning st tau $ !Random Euler angle for cloning nd tau $ !Random Euler angle for cloning nd tau $ !Random Euler angle for cloning nd tau $ !phi of HvecTau1 $ !theta of HvecTau1 $ !phi of HvecTau2 $ !theta of HvecTau2 $ !super key
Definition Taupair.h:42

◆ finalize()

StatusCode EvtDecay::finalize ( )

Definition at line 811 of file EvtDecay.cxx.

812{
813 MsgStream log(messageService(), name());
814 truth.close();
815 log << MSG::INFO << "EvtDecay finalized" << endreq;
816 fstream Forfile;
817 Forfile.open("fort.16",ios::in);
818 char delcmd[300];
819 strcpy(delcmd,"rm -f ");
820 strcat(delcmd,"fort.16");
821 // if(Forfile)system(delcmd);
822
823 fstream Forfile2;
824 Forfile2.open("fort.22",ios::in);
825 strcpy(delcmd,"rm -f ");
826 strcat(delcmd,"fort.22");
827 // if(Forfile2)system(delcmd);
828
829 // To get the branching fraction of the specified mode in Lundcharm Model
830 /*
831 EvtLundCharm lundcharmEvt;
832 int nevt=lundcharmEvt.getTotalEvt();
833 std::cout<<"The total number of the specified mode is "<<nevt<<std::endl;
834 */
835 int totalEvt=0;
836 if(!(m_SB3File=="" || m_SB3HT=="")){
837 for(int i=0;i<500;i++){totalEvt=totalEvt+br[i];}
838 for(int i=0;i<500;i++){
839 double bi=1.0*br[i]/totalEvt;
840 std::cout<<"Branching fraction of channel "<<i<<" is: "<<bi<<std::endl;
841 if(br[i]==0) break;
842 }
843 }
844
845 if(m_statDecays){
846 int totalEvt=0;
847 for(int i=0;i<=totalChannels;i++){totalEvt=totalEvt+br[i];}
848 std::cout<<"==================Statistical first chain decay ==============================="<<std::endl;
849 std::cout<<"i-th channel Events Generated Branching fraction generated "<<std::endl;
850 for(int i=0;i<=totalChannels;i++){
851 std::cout<<"| "<<i<<" "<<br[i]<<" "<<br[i]*1.0/totalEvt<<std::endl;
852
853 }
854 std::cout<<"--------------------------------------------------------------------------------"<<std::endl;
855 std::cout<<" sum: "<<totalEvt<<" "<<"1"<<std::endl;
856 std::cout<<"================== End of statistical first chain decay ========================"<<std::endl;
857 }
858
859 return StatusCode::SUCCESS;
860}

◆ initialize()

StatusCode EvtDecay::initialize ( )

Definition at line 107 of file EvtDecay.cxx.

107 {
108
109 MsgStream log(messageService(), name());
110 if(m_truthFile!=""){truth.open(m_truthFile.c_str());}
111 log << MSG::INFO << "EvtDecay initialize" << endreq;
112 static const bool CREATEIFNOTTHERE(true);
113 StatusCode RndmStatus = service("BesRndmGenSvc", p_BesRndmGenSvc, CREATEIFNOTTHERE);
114 if (!RndmStatus.isSuccess() || 0 == p_BesRndmGenSvc)
115 {
116 log << MSG::ERROR << " Could not initialize Random Number Service" << endreq;
117 return RndmStatus;
118 }
119
120 m_numberEvent=0;
121 first = true;
122 m_ampscalflag = true;
123 // Save the random number seeds in the event
124 CLHEP::HepRandomEngine* engine = p_BesRndmGenSvc->GetEngine("EVTGEN");
125 const long s = engine->getSeed();
126 p_BesRndmGenSvc->setGenseed(s+1);
127
128 m_seeds.clear();
129 m_seeds.push_back(s);
130 totalChannels = 0;
131
132 // open an interface to EvtGen
133
134 if ( m_DecayDec == "") { //use default DECAY.DEC and pdt.table if not defined by user
135 m_DecayDec=getenv("BESEVTGENROOT");
136 m_DecayDec +="/share/DECAY.DEC";
137 }
138
139 if ( m_PdtTable == "") {
140 m_PdtTable =getenv("BESEVTGENROOT");
141 m_PdtTable +="/share/pdt.table";
142 }
143
144 char catcmd[300]; //output user decay cards to log file
145 std::cout<<"===================== user decay card ========================"<<std::endl;
146 strcpy(catcmd, "cat ");
147 strcat(catcmd, userDecFileName.c_str());
148 system(catcmd);
149
152
153 // write decay cards to the root file: pingrg@2009-09-09
154 StatusCode status;
155 status = service("DataInfoSvc",tmpInfoSvc);
156 if (status.isSuccess()) {
157 log << MSG::INFO << "got the DataInfoSvc" << endreq;
158 dataInfoSvc=dynamic_cast<DataInfoSvc *>(tmpInfoSvc);
159 dataInfoSvc->setDecayCard(userDecFileName);
160 } else {
161 log << MSG::ERROR << "could not get the DataInfoSvc" << endreq;
162 return StatusCode::FAILURE;
163 }
164
165
166
167 m_RandomEngine = new EvtBesRandom(engine);
168 m_Gen = new EvtGen( m_DecayDec.c_str(), m_PdtTable.c_str(), m_RandomEngine );
169
170 if(userDecFileName =="NONE") log << MSG::INFO << "EvtDecay User did not define his Decay table EvtGen will use standart" << endreq;
171 if(!(userDecFileName =="NONE")) m_Gen->readUDecay(userDecFileName.c_str());
172
173 if(!(m_DecayTop==" ")) {
174 outfile.open(m_DecayTop.c_str());
175 }
176
177 //for output Ntuple file, pingrg-2009-09-07
178
179
180 if(m_NtupleFile) {
181 NTuplePtr nt1(ntupleSvc(),"MYROOTFILE/Trk");
182 if(nt1) {m_tuple=nt1;}
183 else {
184 m_tuple = ntupleSvc()->book ("MYROOTFILE/Trk", CLID_ColumnWiseTuple, "Generator-trk-Ntuple");
185 if(m_tuple){
186 status = m_tuple->addItem ("TotNumTrk", TotNumTrk, 0,1000);
187 status = m_tuple->addIndexedItem ("Trk_index", TotNumTrk, m_Trk_index);
188 status = m_tuple->addIndexedItem ("m_px_trk", TotNumTrk, m_px_trk);
189 status = m_tuple->addIndexedItem ("m_py_trk", TotNumTrk, m_py_trk);
190 status = m_tuple->addIndexedItem ("m_pz_trk", TotNumTrk, m_pz_trk);
191 status = m_tuple->addIndexedItem ("m_en_trk", TotNumTrk, m_en_trk);
192 status = m_tuple->addIndexedItem ("FST", TotNumTrk, m_fst);
193
194 status = m_tuple->addItem ("nchr", m_nchr);
195 status = m_tuple->addItem ("nchr_e", m_nchr_e);
196 status = m_tuple->addItem ("nchr_mu", m_nchr_mu);
197 status = m_tuple->addItem ("nchr_pi", m_nchr_pi);
198 status = m_tuple->addItem ("nchr_k", m_nchr_k);
199 status = m_tuple->addItem ("nchr_p", m_nchr_p);
200 status = m_tuple->addItem ("N_gamma", m_gamma);
201 status = m_tuple->addItem ("N_gammaFSR", m_gammaFSR);
202 } else {
203 log << MSG::ERROR << " Cannot book N-tuple:"<< long(m_tuple) << endmsg;
204 return StatusCode::FAILURE;
205 }
206 }
207
208 NTuplePtr nt2(ntupleSvc(),"MYROOTFILE/mass");
209 if(nt2) {mass_tuple=nt2;}
210 else {
211 mass_tuple = ntupleSvc()->book ("MYROOTFILE/mass", CLID_ColumnWiseTuple, "Generator-mass-Ntuple");
212 if(mass_tuple){
213 status = mass_tuple->addItem ("m12", m_m12);
214 status = mass_tuple->addItem ("m13", m_m13);
215 status = mass_tuple->addItem ("m23", m_m23);
216 status = mass_tuple->addItem ("m1", m_m1);
217 status = mass_tuple->addItem ("m2", m_m2);
218 status = mass_tuple->addItem ("m3", m_m3);
219 status = mass_tuple->addItem ("costheta1", m_cos1);
220 status = mass_tuple->addItem ("costheta2", m_cos2);
221 status = mass_tuple->addItem ("costheta3", m_cos3);
222 status = mass_tuple->addItem ("ichannel", m_ich);
223 } else {
224 log << MSG::ERROR << " Cannot book N-tuple:"<< long(m_tuple) << endmsg;
225 return StatusCode::FAILURE;
226 }
227 }
228
229 NTuplePtr nt3(ntupleSvc(),"MYROOTFILE/massGen");
230 if(nt3) {massgen_tuple=nt3;}
231 else {
232 massgen_tuple = ntupleSvc()->book ("MYROOTFILE/massGen", CLID_ColumnWiseTuple, "Generator-mass-Ntuple");
233 if(massgen_tuple){
234 status = massgen_tuple->addItem ("m12", _m12);
235 status = massgen_tuple->addItem ("m13", _m13);
236 status = massgen_tuple->addItem ("m23", _m23);
237 status = massgen_tuple->addItem ("m1", _m1);
238 status = massgen_tuple->addItem ("m2", _m2);
239 status = massgen_tuple->addItem ("m3", _m3);
240 status = massgen_tuple->addItem ("costheta1", _cos1);
241 status = massgen_tuple->addItem ("costheta2", _cos2);
242 status = massgen_tuple->addItem ("costheta3", _cos3);
243 status = massgen_tuple->addItem ("ichannel", _ich);
244 } else {
245 log << MSG::ERROR << " Cannot book N-tuple:"<< long(m_tuple) << endmsg;
246 return StatusCode::FAILURE;
247 }
248 }
249
250
251 }
252 for(int i=0;i<500;i++){br[i]=0;}
253 if(!(m_SB3File=="" && m_SB3HT=="")) {
254 const char *filename, *histitle;
255 filename=(char*)m_SB3File.c_str();
256 histitle=(char*)m_SB3HT.c_str();
257 SuperBody3decay.setFile(filename);
258 SuperBody3decay.setHTitle(histitle);
259 SuperBody3decay.init();
260 }
261
262 return StatusCode::SUCCESS;
263}
XmlRpcServer s
INTupleSvc * ntupleSvc()
void setDecayCard(string card)
DataInfoSvc * dataInfoSvc
Definition EvtDecay.h:74
IDataInfoSvc * tmpInfoSvc
Definition EvtDecay.h:73
void readUDecay(const char *const udecay_name)
Definition EvtGen.cc:126
void init()
Definition EvtHis2F.cc:98
void setHTitle(const char *htitle)
Definition EvtHis2F.cc:40
void setFile(const char *dtfile)
Definition EvtHis2F.cc:35
virtual void setGenseed(long)=0
virtual CLHEP::HepRandomEngine * GetEngine(const std::string &StreamName)=0
Interface to the CLHEP engine.

Member Data Documentation

◆ dataInfoSvc

DataInfoSvc* EvtDecay::dataInfoSvc

Definition at line 74 of file EvtDecay.h.

Referenced by initialize().

◆ tmpInfoSvc

IDataInfoSvc* EvtDecay::tmpInfoSvc

Definition at line 73 of file EvtDecay.h.

Referenced by initialize().


The documentation for this class was generated from the following files: