BOSS 7.0.9
BESIII Offline Software System
Loading...
Searching...
No Matches
BesEvtGen-00-04-08/src/EvtDecay.cxx
Go to the documentation of this file.
1//*****************************************************************************
2//
3// EvtDecay.cxx
4//
5// EvtDecay algorithm takes generated event from transient store, and decay
6// particles, including weak decays of its secondary particles.
7// EvtDecay can be used as a TopAlg in conjunction with algorithms Pythia,
8// KKMC or a SingleParticleGun.
9//
10// History:
11// Original LHCb code by Witold Pokorski and Patric Robbe.
12// August 2002: Malte Muller, adopted for ATHENA to work inside algorithm PythiaBModule (only private version within 3.0.0)
13// Sept 2003: James Catmore, adopted for 6.0.3 (not reeased in 6.0.3 ony private version) to work inside PythiaBModule.
14// Nov 2003: M.Smizanska, rewritten a) to be standalone EvtDecay algorithm so it can be combined
15// with any type of Pythia generator b) to match to new LHCb(CDF) code dedicated to LHC, c) to work in 7.3.0.
16// EvtDecay first time released in 7.3.0 20.Nov. 2003
17// Dec 2003: M.Smizanska: define user's input - decay table
18// Feb 2004 J Catmore, addition of random seed facility. TEMPORARY FIX
19// Oct 2005 A. Zhemchugov, adapted for BES3
20// May 2006 K.L He, set spin density matrix for e+e- -> V
21// Sept.2007, R.G.Ping, change to the BesEvtGen
22// Jan. 2008, R.G.Ping, to print McTruth table to the file DecayTopology when only doing generation, not simulation
23// Feb. 2008, R.G.Ping, add an option to only run the BesEvtGen
24// Mar. 2008, R.G. Ping, user options "DecayDecDir" and "PdtTableDir" are changed.
25// Mar. 2008-03, R.G. Ping, dump the user decay card to the bosslog file
26//*****************************************************************************
27
28// Header for this module:-
29#include "EvtDecay.h"
30#include "ReadME.h"
31#include "EvtGenBase/EvtVector4R.hh"
32#include "EvtGenBase/EvtParticle.hh"
33#include "EvtGenBase/EvtParticleFactory.hh"
34#include "EvtGenBase/EvtDecayTag.hh"
35#include "EvtGen.hh"
36#include "EvtGenBase/EvtRandomEngine.hh"
37#include "EvtGenBase/EvtDecayTable.hh"
38#include "EvtGenModels/EvtLundCharm.hh"
39#include "EvtGenModels/EvtPsi3Sdecay.hh"
40#include "EvtGlobalSet.hh"
41#include "EvtGenBase/EvtHelSys.hh"
42#include "EvtGenModels/EvtCalHelAmp.hh"
43#include "EvtGenModels/EvtConExc.hh"
44
45// Framework Related Headers:-
46#include "HepMC/GenEvent.h"
47#include "HepMC/GenVertex.h"
48#include "HepMC/GenParticle.h"
49#include "EventModel/EventHeader.h"
50#include "GaudiKernel/SmartDataPtr.h"
51#include "DataInfoSvc/IDataInfoSvc.h"
52#include "DataInfoSvc/DataInfoSvc.h"
53#include "GaudiKernel/MsgStream.h"
54#include "GaudiKernel/ISvcLocator.h"
55#include "GaudiKernel/AlgFactory.h"
56#include "GaudiKernel/DataSvc.h"
57#include "GaudiKernel/SmartDataPtr.h"
58#include "GaudiKernel/IPartPropSvc.h"
59
60#include "BesKernel/IBesRndmGenSvc.h"
61#include "GeneratorObject/McGenEvent.h"
62#include "CLHEP/Random/Ranlux64Engine.h"
63#include "CLHEP/Random/RandFlat.h"
64#include "EvtGenBase/EvtEulerAngles.hh"
65
66#include <stdlib.h>
67#include <string.h>
68
69static string mydecayrec;
70static double _amps_max;
71int nchr = 0;
72int nchr_e = 0;
73int nchr_mu= 0;
74int nchr_pi= 0;
75int nchr_k = 0;
76int nchr_p = 0;
77int N_gamma=0;
78int N_gammaFSR = 0;
79int fst[50];
80int EvtDecay::m_runNo=0;
81double EvtConExc::SetMthr=0;
82
83EvtDecay::EvtDecay(const string& name, ISvcLocator* pSvcLocator):Algorithm( name, pSvcLocator )
84{
85
86
87
88 // these can be used to specify alternative locations or filenames
89 // for the EvtGen particle and channel definition files.
90
91 declareProperty("DecayDecDir", m_DecayDec = "");
92 declareProperty("PdtTableDir", m_PdtTable = "");
93 declareProperty("userDecayTableName", userDecFileName = "NONE");
94 declareProperty("DecayTopology", m_DecayTop = ""); // output decay topology to a file specified by user
95 declareProperty("DecayRecord", m_DecayRec = ""); // output decay topology of one particle specified by user to a file
96 declareProperty("ParentParticle", m_ParentPart = "NONE");// Mothor particle name in pdt.table for only running BesEvtGen
97 declareProperty("Multiplicity", m_Ncharge = false); // output ncharge number of an event
98 declareProperty("NutupleFile",m_NtupleFile = false); // output Ntuple file
99 declareProperty("mDIY",_mDIY= false); // mDIY model
100 declareProperty("FDPdata",m_SB3File = ""); // Fit Dalitz Plot (FDP) data file name (*.root)
101 declareProperty("FDPHisT",m_SB3HT = ""); // histogram title of Dalitz plot in data file
102 declareProperty("FDPpart",m_FDPparticle =""); //to assign the FDP parent particle name (string)
103 declareProperty("TruthFile",m_truthFile ="");
104 declareProperty("TruthPart",m_truthPart ="");
105 declareProperty("Psi3SopenCharm",m_Psi4040OpenCharm=false);
106 declareProperty("Psi2openCharm", m_Psi2openCharm=false);
107 declareProperty("SetMthrForConExc",m_SetMthr=0.0);
108 declareProperty("statDecays",m_statDecays=false);
109 declareProperty("TagLundCharmModel", m_tagLundModel=false);
110 m_mystring.clear();
111 declareProperty("FileForTrackGen", m_mystring);
112 //for polarized charmonium production
113 m_polarization.clear();//= diag{1+Pe, 0, 1-Pe} with electron polarization Pe
114 declareProperty("polarization", m_polarization);
115 declareProperty("eBeamPolarization", m_eBeamPolarization=-1);
116 m_cluster0.clear();m_wind0.clear();
117 m_cluster1.clear();m_wind1.clear();
118 m_cluster2.clear();m_wind2.clear();
119 declareProperty("cluster0",m_cluster0);
120 declareProperty("cluster1",m_cluster1);
121 declareProperty("cluster2",m_cluster2);
122 declareProperty("masswin0",m_wind0);
123 declareProperty("masswin1",m_wind1);
124 declareProperty("masswin2",m_wind2);
125 declareProperty("OutputP4",m_outputp4="");
126 declareProperty("ReadMeasuredEcms", m_RdMeasuredEcms = false);
127 declareProperty("beamEnergySpread", m_beamEnergySpread = 0);
128 m_ReadTruth.clear();
129 declareProperty("ReadTruthAtCM",_truthAtCM=false);
130 declareProperty("ReadTruth",m_ReadTruth);
131 //ReadTruth={{"ParentName"},{"i0","i1","i2"},{"j0","j1","j2","j3"}}, where the first part. is Parent->getDaug(i0)->getDaug(i1)->getDaug(i2),
132 //and second particle is Parent ->getDaug(j0)->getDaug(j1)->getDaug(j2)->getDaug(j3);
133 declareProperty("RvalueTag",_RvalueTag=false);
134 declareProperty("PrintFSFor",m_printfs="");
135}
136
137
138StatusCode EvtDecay::initialize(){
139
140 MsgStream log(messageService(), name());
141 system("cat $BESEVTGENROOT/share/phokhara_9.1.fferr>phokhara_9.1.fferr");
142 system("cat $BESEVTGENROOT/share/phokhara_9.1.ffwarn>phokhara_9.1.ffwarn");
143
144 if(m_truthFile!=""){truth.open(m_truthFile.c_str());}
145 log << MSG::INFO << "EvtDecay initialize" << endreq;
146 static const bool CREATEIFNOTTHERE(true);
147 StatusCode RndmStatus = service("BesRndmGenSvc", p_BesRndmGenSvc, CREATEIFNOTTHERE);
148 if (!RndmStatus.isSuccess() || 0 == p_BesRndmGenSvc)
149 {
150 log << MSG::ERROR << " Could not initialize Random Number Service" << endreq;
151 return RndmStatus;
152 }
153
154 EvtGlobalSet::SV.clear();
155 for(int i=0;i<m_mystring.size();i++){EvtGlobalSet::SV.push_back(m_mystring[i]);}
156
157 EvtGlobalSet::iVV.clear();
158 EvtGlobalSet::dVV.clear();
159 std::vector<std::vector<int> >myivv;
160 std::vector<std::vector<double> >mydvv;
161
162 EvtConExc::SetMthr=m_SetMthr;
163
164 myivv.push_back(m_cluster0);
165 myivv.push_back(m_cluster1);
166 myivv.push_back(m_cluster2);
167
168 mydvv.push_back(m_wind0);
169 mydvv.push_back(m_wind1);
170 mydvv.push_back(m_wind2);
171
172 for(int i=0;i<myivv.size();i++){
173 std::vector<int> theivv;
174 for(int j=0;j<myivv[i].size();j++){
175 theivv.push_back(myivv[i][j]);
176 }
177 if(theivv.size()) EvtGlobalSet::iVV.push_back(theivv);
178 }
179
180 for(int i=0;i<mydvv.size();i++){
181 std::vector<double> thedvv;
182 for(int j=0;j<mydvv[i].size();j++){
183 thedvv.push_back(mydvv[i][j]);
184 }
185 if(thedvv.size()) EvtGlobalSet::dVV.push_back(thedvv);
186 }
187
188
189 if(m_eBeamPolarization>1) {std::cout<<"The e-beam polaziation should be in 0 to 1"<<std::endl;abort();}
190 m_numberEvent=0;
191 first = true;
192 m_ampscalflag = true;
193 // Save the random number seeds in the event
194 CLHEP::HepRandomEngine* engine = p_BesRndmGenSvc->GetEngine("EVTGEN");
195 const long s = engine->getSeed();
196 p_BesRndmGenSvc->setGenseed(s+1);
197
198 m_seeds.clear();
199 m_seeds.push_back(s);
200 totalChannels = 0;
201
202 // open an interface to EvtGen
203
204 if ( m_DecayDec == "") { //use default DECAY.DEC and pdt.table if not defined by user
205 m_DecayDec=getenv("BESEVTGENROOT");
206 m_DecayDec +="/share/DECAY.DEC";
207 }
208
209 if ( m_PdtTable == "") {
210 m_PdtTable =getenv("BESEVTGENROOT");
211 m_PdtTable +="/share/pdt.table";
212 }
213
214 char catcmd[300]; //output user decay cards to log file
215 std::cout<<"===================== user decay card ========================"<<std::endl;
216 strcpy(catcmd, "cat ");
217 strcat(catcmd, userDecFileName.c_str());
218 system(catcmd);
219
220 IDataInfoSvc *tmpInfoSvc;
221 DataInfoSvc* dataInfoSvc;
222
223 // write decay cards to the root file: pingrg@2009-09-09
224 StatusCode status;
225 status = service("DataInfoSvc",tmpInfoSvc);
226 if (status.isSuccess()) {
227 log << MSG::INFO << "got the DataInfoSvc" << endreq;
228 dataInfoSvc=dynamic_cast<DataInfoSvc *>(tmpInfoSvc);
229 dataInfoSvc->setDecayCard(userDecFileName);
230 } else {
231 log << MSG::ERROR << "could not get the DataInfoSvc" << endreq;
232 return StatusCode::FAILURE;
233 }
234
235
236
237 m_RandomEngine = new EvtBesRandom(engine);
238 m_Gen = new EvtGen( m_DecayDec.c_str(), m_PdtTable.c_str(), m_RandomEngine );
239
240 if(userDecFileName =="NONE") log << MSG::INFO << "EvtDecay User did not define his Decay table EvtGen will use standart" << endreq;
241 if(!(userDecFileName =="NONE")) m_Gen->readUDecay(userDecFileName.c_str());
242
243 if(!(m_DecayTop==" ")) {
244 outfile.open(m_DecayTop.c_str());
245 }
246
247 //for output Ntuple file, pingrg-2009-09-07
248
249
250 if(m_NtupleFile) {
251 NTuplePtr nt1(ntupleSvc(),"MYROOTFILE/Trk");
252 if(nt1) {m_tuple=nt1;}
253 else {
254 m_tuple = ntupleSvc()->book ("MYROOTFILE/Trk", CLID_ColumnWiseTuple, "Generator-trk-Ntuple");
255 if(m_tuple){
256 status = m_tuple->addItem ("TotNumTrk", TotNumTrk, 0,100);
257 status = m_tuple->addIndexedItem ("Trk_index", TotNumTrk, m_Trk_index);
258 status = m_tuple->addIndexedItem ("m_px_trk", TotNumTrk, m_px_trk);
259 status = m_tuple->addIndexedItem ("m_py_trk", TotNumTrk, m_py_trk);
260 status = m_tuple->addIndexedItem ("m_pz_trk", TotNumTrk, m_pz_trk);
261 status = m_tuple->addIndexedItem ("m_en_trk", TotNumTrk, m_en_trk);
262 status = m_tuple->addIndexedItem ("FST", TotNumTrk, m_fst);
263
264 status = m_tuple->addItem ("nchr", m_nchr);
265 status = m_tuple->addItem ("nchr_e", m_nchr_e);
266 status = m_tuple->addItem ("nchr_mu", m_nchr_mu);
267 status = m_tuple->addItem ("nchr_pi", m_nchr_pi);
268 status = m_tuple->addItem ("nchr_k", m_nchr_k);
269 status = m_tuple->addItem ("nchr_p", m_nchr_p);
270 status = m_tuple->addItem ("N_gamma", m_gamma);
271 status = m_tuple->addItem ("N_gammaFSR", m_gammaFSR);
272 status = m_tuple->addItem ("Flag1", m_flag1);
273 } else {
274 log << MSG::ERROR << " Cannot book N-tuple:"<< long(m_tuple) << endmsg;
275 return StatusCode::FAILURE;
276 }
277 }
278
279 NTuplePtr nt2(ntupleSvc(),"MYROOTFILE/mass");
280 if(nt2) {mass_tuple=nt2;}
281 else {
282 mass_tuple = ntupleSvc()->book ("MYROOTFILE/mass", CLID_ColumnWiseTuple, "Generator-mass-Ntuple");
283 if(mass_tuple){
284 status = mass_tuple->addItem ("m12", m_m12);
285 status = mass_tuple->addItem ("m13", m_m13);
286 status = mass_tuple->addItem ("m23", m_m23);
287 status = mass_tuple->addItem ("m1", m_m1);
288 status = mass_tuple->addItem ("m2", m_m2);
289 status = mass_tuple->addItem ("m3", m_m3);
290 status = mass_tuple->addItem ("costheta1", m_cos1);
291 status = mass_tuple->addItem ("costheta2", m_cos2);
292 status = mass_tuple->addItem ("costheta3", m_cos3);
293 status = mass_tuple->addItem ("ichannel", m_ich);
294 } else {
295 log << MSG::ERROR << " Cannot book N-tuple:"<< long(m_tuple) << endmsg;
296 return StatusCode::FAILURE;
297 }
298 }
299
300 NTuplePtr nt3(ntupleSvc(),"MYROOTFILE/massGen");
301 if(nt3) {massgen_tuple=nt3;}
302 else {
303 massgen_tuple = ntupleSvc()->book ("MYROOTFILE/massGen", CLID_ColumnWiseTuple, "Generator-mass-Ntuple");
304 if(massgen_tuple){
305 status = massgen_tuple->addItem ("m12", _m12);
306 status = massgen_tuple->addItem ("m13", _m13);
307 status = massgen_tuple->addItem ("m23", _m23);
308 status = massgen_tuple->addItem ("m1", _m1);
309 status = massgen_tuple->addItem ("m2", _m2);
310 status = massgen_tuple->addItem ("m3", _m3);
311 status = massgen_tuple->addItem ("costheta1", _cos1);
312 status = massgen_tuple->addItem ("costheta2", _cos2);
313 status = massgen_tuple->addItem ("costheta3", _cos3);
314 status = massgen_tuple->addItem ("ichannel", _ich);
315 } else {
316 log << MSG::ERROR << " Cannot book N-tuple:"<< long(m_tuple) << endmsg;
317 return StatusCode::FAILURE;
318 }
319 }
320
321
322 }
323 for(int i=0;i<500;i++){br[i]=0;}
324 if(!(m_SB3File=="" && m_SB3HT=="")) {
325 const char *filename, *histitle;
326 filename=(char*)m_SB3File.c_str();
327 histitle=(char*)m_SB3HT.c_str();
328 SuperBody3decay.setFile(filename);
329 SuperBody3decay.setHTitle(histitle);
330 SuperBody3decay.init();
331 }
332
333 return StatusCode::SUCCESS;
334}
335
336StatusCode EvtDecay::execute()
337{
338
339 MsgStream log(messageService(), name());
340 // log << MSG::INFO << "EvtDecay executing" << endreq;
341 //------------------
342 string key = "GEN_EVENT";
343 // retrieve event from Transient Store (Storegate)
344 SmartDataPtr<McGenEventCol> McEvtColl(eventSvc(), "/Event/Gen");
345
346 m_numberEvent += 1;
347 writeFlag = 0;
348 //std::cout<<"EvtNumber= "<<m_numberEvent<<std::endl;
349 if ( McEvtColl == 0 )
350 {
351 HepMC::GenEvent* evt=new GenEvent();
352 evt->set_event_number(m_numberEvent);
353
354 //read Ecms from the database
355 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
356 int runNo=eventHeader->runNumber();
357 int event=eventHeader->eventNumber();
358 bool newRunFlag=0;
359 if(runNo != 0 && runNo != m_runNo){m_runNo=runNo;newRunFlag = true;}else{newRunFlag=false;}
360 if(m_RdMeasuredEcms&& newRunFlag){// using cms energy of beam read from database
361 runNo = std::abs(runNo);
362 ReadME theME(runNo);
363 if(theME.isRunNoValid()){
364 dbEcms=theME.getEcms();
365 std::cout<<"Read Ecms= "<<dbEcms<<std::endl;
366 if(dbEcms!=0){
367 std::cout<<"INFO: Read the MeasuredEcms"<<std::endl;
368 }
369 else{
370 std::cout<<"ERROR: Can not read the MeasuredEcms, try to turn off the ReadEcmsDB"<<std::endl;
371 return StatusCode::FAILURE;
372 }
373 }
374 }
375 //end of read Ecms
376
377 callBesEvtGen( evt );
378 if(!(m_DecayTop=="")) evt->print(outfile);
379
380 //create a Transient store
381 McGenEventCol *mcColl = new McGenEventCol;
382 McGenEvent* mcEvent = new McGenEvent(evt);
383 mcColl->push_back(mcEvent);
384 StatusCode sc = eventSvc()->registerObject("/Event/Gen",mcColl);
385 if(sc != SUCCESS) return StatusCode::FAILURE;
386
387 } else // (McEvtColl != 0)
388 {
389 McGenEventCol::iterator mcItr;
390 for( mcItr = McEvtColl->begin(); mcItr != McEvtColl->end(); mcItr++ )
391 {
392 HepMC::GenEvent* hEvt = (*mcItr)->getGenEvt();
393 // MeVToGeV( hEvt );
394
395 callEvtGen( hEvt );
396
397 if(!(m_DecayTop=="")) hEvt->print(outfile);
398 // GeVToMeV( hEvt );
399 // if(!(m_DecayRec=="")) outfile2<<std::endl;
400 if(!(m_DecayRec=="")) std::cout<<std::endl;
401 };
402 }
403 if( m_NtupleFile ){ m_tuple->write();}
404
405 if(_mDIY){
406 std::cout<<"helangs ";
407 for(int i=0;i<m_vangs.size();i++) std::cout<<m_vangs[i]<<" ";
408 std::cout<<std::endl;
409 }
410 return StatusCode::SUCCESS;
411}
412
413// main procedure, loops over all particles in given event,
414// converts them to EvtGenParticles,
415// feeds them to EvtGen,
416// converts back to HepMC particles and puts them in the event.
417
418StatusCode EvtDecay::callEvtGen( HepMC::GenEvent* hepMCevt )
419{
420 MsgStream log(messageService(), name());
421 nchr = 0;
422 nchr_e = 0;
423 nchr_mu= 0;
424 nchr_pi= 0;
425 nchr_k = 0;
426 nchr_p = 0;
427 N_gamma=0;
428 N_gammaFSR = 0;
429 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
430
431 for (int i=0;i<13;i++) fst[i]=0;
432 HepMC::GenEvent::particle_const_iterator itp;
433 AllTrk_index=0;
434 for( itp=hepMCevt->particles_begin(); itp!=hepMCevt->particles_end(); ++itp )
435 {
436 // This is the main loop.
437 // We iterate over particles and we look for ones that
438 // should be decayed with EvtGen.
439 //
440 // status == 1 - undecayed Pythia particle (also for particles that are not supposed to decay)
441 // status == 999 - particle already decayed by EvtGen
442 // status == 899 - particle was supposed to decay by EvtGen - but found no model
443 // this may be also intentional - such events are then excluded from
444 // being written into persistent output.
445 // ISTHEP(IHEP):
446 // status code for entry IHEP, with the following meanings:
447 //
448 // = 0 : null entry.
449 // = 1 : an existing entry, which has not decayed or fragmented. This is the main class of entries, which represents the `final state' given by the generator.
450 // = 2 : an entry which has decayed or fragmented and is therefore not appearing in the final state, but is retained for event history information.
451 // = 3 : a documentation line, defined separately from the event history. This could include the two incoming reacting particles, etc.
452 // = 4 - 10 : undefined, but reserved for future standards.
453 // = 11 - 200 : at the disposal of each model builder for constructs specific to his program, but equivalent to a null line in the context of any other program.
454 // = 201 - : at the disposal of users, in particular for event tracking in the detector.
455
456 HepMC::GenParticle* hepMCpart=*itp;
457 int stat = hepMCpart->status();
458
459
460 if ( stat != 1 ) continue;
461
462 loop_to_select_eventsB:
463 int id = hepMCpart->pdg_id();
464 parentPDGcode=id;
465 hepMCpart->set_status(899);
466 EvtId eid=EvtPDL::evtIdFromStdHep(id);
467 string parentName=EvtPDL::name(eid);
468
469 double en =(hepMCpart->momentum()).e();
470 double pex=(hepMCpart->momentum()).px();
471 double pey=(hepMCpart->momentum()).py();
472 double pez=(hepMCpart->momentum()).pz();
473
474 EvtVector4R p_init(en,pex,pey,pez);
475 parentMass=p_init.mass();
476 EvtPDL::reSetMass(eid,parentMass);
477
478 EvtParticle* part=EvtParticleFactory::particleFactory(eid,p_init);
479
480 // set spin density matrix for e+ e- -> V
481 bool parentFlag=false;
482 if( writeFlag==0 && part->getSpinType() == EvtSpinType::VECTOR) {
483 if(m_polarization.size()>0) {
484 part->setPolarizedSpinDensity(m_polarization[0],m_polarization[1],m_polarization[2]);
485 }else if(m_eBeamPolarization>0){
486 part->setPolarizedSpinDensity(1+ m_eBeamPolarization,0,1- m_eBeamPolarization);
487 } else{
488 part->setVectorSpinDensity();
489 }
490 parentFlag=true;
491 writeFlag++;
492 } else {
493 int id = hepMCpart->pdg_id();
494 if( id == 22 ) {N_gamma ++;fst[0]++;} //fst[0]:photon
495 if( id == -22 ){N_gammaFSR ++;}
496 if( id == 11 ) {fst[1]++;} else if(id == -11){fst[2]++;} // fst[1]: electron ; fst[]: positron
497 else if(id == 13){fst[3]++;} // fst[3]: mu-
498 else if(id == -13){fst[4]++;} // fst[4]: mu+
499 else if(id == 211){fst[5]++;} // fst[5]: pi+
500 else if(id ==-211){fst[6]++;} // fst[6]: pi-
501 else if(id == 321){fst[7]++;} // fst[7]: K+
502 else if(id ==-321){fst[8]++;} // fst[8]: K-
503 else if(id ==2212){fst[9]++;} // fst[9]: pronton
504 else if(id ==-2212){fst[10]++;} // fst[10]: anti-proton
505 else if(id == 2112){fst[11]++;} // fst[11]: nucleon
506 else if(id ==-2112){fst[12]++;} // fst[12]: ant-nucleon
507 if( fabs(id) == 11) {nchr_e++;}
508 if( fabs(id) == 13) {nchr_mu++;}
509 if( fabs(id) == 211) {nchr_pi++;}
510 if( fabs(id) == 321) {nchr_k++;}
511 if( fabs(id) == 2212) {nchr_p++;}
512
513 //std::cout<<"id= "<<id<<" AllTrk_index= "<<AllTrk_index<<std::endl;
514 if( m_NtupleFile==true ){
515 m_Trk_index[AllTrk_index]=id;
516 m_px_trk[AllTrk_index]=pex;
517 m_py_trk[AllTrk_index]=pey;
518 m_pz_trk[AllTrk_index]=pez;
519 m_en_trk[AllTrk_index]=en ;
520
521 AllTrk_index++;
522 Trk_index[AllTrk_index]=0;
523 }
524
525
526 }
527
528 // for SuperBody3decay generator
529 // EvtVector4R pd1,pd2,pd3;
530 EvtVector4R fdp_init;
531 EvtId fdp_ppid;
532 if(m_FDPparticle!=""){
533 findPart(part);
534 fdp_ppid = FDP_id;
535 fdp_init = FDP_init;
536 } else{
537 fdp_ppid = eid;
538 fdp_init = p_init;
539 }
540
541 if(!(m_SB3File=="" || m_SB3HT=="") && parentFlag && first ){
542 SuperBody3decay_make(fdp_ppid, fdp_init );
543 }
544
545 loop_to_select_eventsA:
546 m_Gen->generateDecay(part);
547 if(m_truthFile!="" && m_truthPart!=""){
548 if(EvtPDL::name(part->getId())==m_truthPart ){
549 int ndaug = part->getNDaug();
550 for(int id=0;id<ndaug;id++){
551 EvtParticle* sonp=part->getDaug(id);
552 EvtVector4R son=sonp->getP4();
553 truth<<son.get(0)<<" "<<son.get(1)<<" "<<son.get(2)<<" "<<son.get(3)<<std::endl;
554 }
555 }
556 }
557 double ratio,rdm;
558 if(_mDIY && parentFlag && m_ampscalflag ) _amps_max=CalAmpsMax(part);
559 if(_mDIY && parentFlag) {
560 rdm=EvtRandom::Flat(0.0, 1.0);
561 double amps=CalAmpsMDIY( part );
562 ratio=amps/_amps_max;
563 }
564 if(_mDIY && parentFlag && ratio<=rdm){ goto loop_to_select_eventsA;}
565
566 //std::cout<<"m_vangs.size= "<<m_vangs.size()<<std::endl;
567 //-------- For superbody3------------------------------
568 bool accept;
569 if(m_FDPparticle==""){FDP_part=part;}
570 if(!(m_SB3File=="" || m_SB3HT=="") && parentFlag ){ accept=SuperBody3decay_judge( FDP_part); }
571 if(!(m_SB3File=="" || m_SB3HT=="") && parentFlag && !accept){
572 part->deleteTree();
573 writeFlag=0;
574 goto loop_to_select_eventsB;
575 }else if(!(m_SB3File=="" || m_SB3HT=="") && parentFlag && accept){
576 countChannel(FDP_part);
577 }
578 //-------- for psi4040 open charm inclusive decay
579 if((m_Psi2openCharm || m_Psi4040OpenCharm) && parentFlag ){
580 if(getModel(part) == "OPENCHARM"){
581 std::cout<<"OPENCHARM model need to set Psi2openCharm and Psi3SopenCharm to be false!"<<std::endl;
582 abort();
583 }
584 EvtPsi3Sdecay opencharm(en, part);
585 bool theDecay = opencharm.choseDecay();
586 if(!theDecay ) {
587 part->deleteTree();
588 writeFlag=0;
589 goto loop_to_select_eventsB;
590 }
591 }
592 //------------- dump the decay tree to the event model
593 // if(Nchannel>=0) part->dumpTree();
594 // part->printTree();
595
596 if(parentFlag){
597 EvtDecayTag theTag(part);
598 unsigned int modeTag, multiplicityTag;
599 modeTag = theTag.getModeTag();
600 if( getModel(part) == "OPENCHARM"||getModel(part) == "LUNDCHARM" && m_tagLundModel){
601 multiplicityTag = decayType(part);
602 //debugging
603 // std::cout<<"Opencharm Mode: "<<decayType(part)<<std::endl;
604 } else {
605 multiplicityTag = theTag.getMultTag() + decayType(part);
606 }
607
608 //std::cout<<"StdHep= "<<EvtPDL::getStdHep(part->getId())<<" getChannel "<<part->getChannel()<<" flag1,flag2= "<<modeTag<<" "<<multiplicityTag<<std::endl;
609 if(eventHeader == 0) std::cout << "event header unavailable: " << std::endl;
610 if (eventHeader != 0) {
611 eventHeader->setFlag1(modeTag);
612 eventHeader->setFlag2(multiplicityTag);
613 //std::cout<<"flag1,flag2= "<<modeTag<<" "<<multiplicityTag<<std::endl;
614 }
615 }
616
617 if(!(m_DecayRec=="")) mydecayrec=part->writeTreeRec(m_DecayRec);
618 //-- decay statistics
619 if(m_statDecays && parentFlag ) countChannel(part);
620 // ----------- pingrg 2008-03-25 ---------
621
622 if ( log.level() <= MSG::DEBUG ) part->printTree() ;
623 log << MSG::DEBUG << "Converting particles to HepMC" << endreq;
624
625 makeHepMC(part, hepMCevt, hepMCpart);
626 if(part->getNDaug()!=0)
627 hepMCpart->set_status(999);
628
629 //debug
630
631 if(part->getId()==EvtPDL::getId(m_outputp4)){
632 int nds=part->getNDaug();
633 for(int i=0;i<nds;i++){
634 if(EvtPDL::name(part->getDaug(i)->getId())=="gammaFSR") continue;
635 EvtVector4R vp1=part->getDaug(i)->getP4Lab();
636 std::cout<<"vvpp "<<vp1.get(0)<<" "<<vp1.get(1)<<" "<<vp1.get(2)<<" "<<vp1.get(3)<<std::endl;
637 }
638 }
639
640 if(m_ReadTruth.size())ReadTruth(part,m_ReadTruth);
641 //if(part->getId()==EvtPDL::getId("psi(4260)")) std::cout<<"mydebug "<<part->getP4().mass()<<std::endl;
642 part->deleteTree();
643 };
644
645 for( itp=hepMCevt->particles_begin(); itp!=hepMCevt->particles_end(); ++itp )
646 {
647 if (!(*itp)->end_vertex() && ( (*itp)->status()==899 || (*itp)->status()==999 ))
648 (*itp)->set_status(1);
649 }
650 //nchr = nchr_e + nchr_mu + nchr_pi + nchr_k +nchr_p;
651 nchr = nchr_pi + nchr_k +nchr_p;
652 if(m_Ncharge == true ) {std::cout<<"Multiplicity: TotoalCharge_e_mu_pi_K_p_gamma_gammaFSR: "
653 <<nchr <<" "
654 <<nchr_e <<" "
655 << nchr_mu <<" "
656 << nchr_pi <<" "
657 << nchr_k <<" "
658 <<nchr_p <<" "
659 <<N_gamma <<" "
660 <<N_gammaFSR<<endl;}
661 if(m_Ncharge == true ){std::cout<<"FinalStates: "
662 <<fst[0]<<" "
663 <<fst[1]<<" "
664 <<fst[2]<<" "
665 <<fst[3]<<" "
666 <<fst[4]<<" "
667 <<fst[5]<<" "
668 <<fst[6]<<" "
669 <<fst[7]<<" "
670 <<fst[8]<<" "
671 <<fst[9]<<" "
672 <<fst[10]<<" "
673 <<fst[11]<<" "
674 <<fst[12]<<std::endl;}
675 if( m_NtupleFile ){
676 // TotNumTrk=500;
677 m_nchr = nchr;
678 m_nchr_e = nchr_e;
679 m_nchr_mu = nchr_mu;
680 m_nchr_pi = nchr_pi;
681 m_nchr_k = nchr_k;
682 m_nchr_p = nchr_p;
683 m_gamma=N_gamma;
684 m_gammaFSR=N_gammaFSR;
685 TotNumTrk = AllTrk_index;// nchr + N_gamma + N_gammaFSR;
686 }
687 // std::cout<<"Flag= "<<eventHeader->flag1()<<" "<<eventHeader->flag2()<<std::endl;
688
689 return StatusCode::SUCCESS;
690}
691
692//****** CallBesEvtGen ************
693// This part is responsible for only ruuing BesEvtGen. //pingrg Feb. 3,2008
694//***************************************************
695StatusCode EvtDecay::callBesEvtGen( HepMC::GenEvent* hepMCevt )
696{
697 MsgStream log(messageService(), name());
698
699 nchr = -1;
700 nchr_e = 0;
701 nchr_mu= 0;
702 nchr_pi= 0;
703 nchr_k = 0;
704 nchr_p = 0;
705 N_gamma=0;
706 N_gammaFSR = 0;
707 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
708 for (int i=0;i<13;i++) fst[i]=0;
709 HepMC::GenEvent::particle_const_iterator itp;
710 // This is the main loop.
711 // We iterate over particles and we look for ones that
712 // should be decayed with EvtGen.
713 //
714 // status == 1 - undecayed Pythia particle (also for particles that are not supposed to decay)
715 // status == 999 - particle already decayed by EvtGen
716 // status == 899 - particle was supposed to decay by EvtGen - but found no model
717 // this may be also intentional - such events are then excluded from
718 // being written into persistent output.
719 // ISTHEP(IHEP):
720 // status code for entry IHEP, with the following meanings:
721 //
722 // = 0 : null entry.
723 // = 1 : an existing entry, which has not decayed or fragmented. This is the main class of entries, which represents the `final state' given by the generator.
724 // = 2 : an entry which has decayed or fragmented and is therefore not appearing in the final state, but is retained for event history information.
725 // = 3 : a documentation line, defined separately from the event history. This could include the two incoming reacting particles, etc.
726 // = 4 - 10 : undefined, but reserved for future standards.
727 // = 11 - 200 : at the disposal of each model builder for constructs specific to his program, but equivalent to a null line in the context of any other program.
728 // = 201 - : at the disposal of users, in particular for event tracking in the detector.
729
730 loop_to_select_eventsB:
731 EvtId ppid=EvtPDL::getId(m_ParentPart); //parent particle name
732 double ppmass=0;
733 if(m_RdMeasuredEcms && m_beamEnergySpread==0) { //read Ecms, no energy spread
734 EvtPDL::reSetMass(ppid,dbEcms);
735 ppmass=EvtPDL::getMass(ppid);
736 }else if(m_RdMeasuredEcms && m_beamEnergySpread!=0){ //read Ecms, with nergy spread
737 double cmssig=sqrt(2.)*m_beamEnergySpread;
738 double myppmass = energySpread(dbEcms, cmssig);
739 EvtPDL::reSetMass(ppid,myppmass);
740 ppmass=EvtPDL::getMass(ppid);
741 }else if( (!m_RdMeasuredEcms) && m_beamEnergySpread!=0){//not read Ecms, with energy spread
742 if(m_numberEvent==1 )dbEcms = EvtPDL::getMass(ppid);
743 double cmssig=sqrt(2.)*m_beamEnergySpread;
744 double myppmass = energySpread(dbEcms, cmssig);
745 EvtPDL::reSetMass(ppid,myppmass);
746 ppmass=EvtPDL::getMass(ppid);
747 }else if((!m_RdMeasuredEcms) && m_beamEnergySpread==0 ){//not read Ecms, no energy spread
748 ppmass=EvtPDL::getMass(ppid);
749 }else {std::cout<<"EvtDecay::callBesEvtGen:: bad option"<<std::endl; abort();}
750
751 // std::cout<<"testMass "<<ppmass<<std::endl;
752 //using Ecms from data base, for XYZ data sets
753 parentMass=ppmass;
754 int pid=EvtPDL::getStdHep(ppid);
755 parentPDGcode=pid;
756
757 EvtVector4R p_init(ppmass,0.0,0.0,0.0);
758
759 EvtParticle* part=EvtParticleFactory::particleFactory(ppid,p_init);
760 HepMC::GenParticle* hepMCpart=new GenParticle(HepLorentzVector(0,0,0,ppmass),pid,3);
761
762 bool parentFlag = false;
763 // set spin density matrix for e+ e- -> V
764 if(writeFlag ==0 && part->getSpinType() == EvtSpinType::VECTOR) {
765 if(m_polarization.size()>0) {
766 part->setPolarizedSpinDensity(m_polarization[0],m_polarization[1],m_polarization[2]);
767 } else if(m_eBeamPolarization>0){
768 part->setPolarizedSpinDensity(1+ m_eBeamPolarization,0,1- m_eBeamPolarization);
769 } else{
770 part->setVectorSpinDensity();
771 }
772 parentFlag = true;
773 writeFlag++;
774 }
775
776 // for SuperBody3decay generator
777 EvtVector4R fdp_init;
778 EvtId fdp_ppid;
779 if(m_FDPparticle!=""){
780 findPart(part);
781 fdp_ppid = FDP_id;
782 fdp_init = FDP_init;
783 } else{
784 fdp_ppid = ppid;
785 fdp_init = p_init;
786 }
787
788 if(!(m_SB3File=="" || m_SB3HT=="") && first ){ SuperBody3decay_make( fdp_ppid, fdp_init);}
789 // -----------------------------------
790
791 loop_to_select_events:
792 m_Gen->generateDecay(part); if(m_numberEvent<=1){ std::cout<<"after m_Gen->decay "<<std::endl; part->printTree();}
793
794 double ratio,rdm;
795 if(_mDIY && m_ampscalflag ) _amps_max=CalAmpsMax(part);
796
797 if(_mDIY ) {
798 for(int myloop=0;myloop<100;myloop++){
799 m_Gen->generateDecay(part);
800 double amps=CalAmpsMDIY( part);
801 ratio=amps/_amps_max;
802 rdm=EvtRandom::Flat(0.0, 1.0);
803 if( !isNumber(amps) || !isNumber(_amps_max) || ratio<=0 ) {
804 part->deleteTree();
805 part=EvtParticleFactory::particleFactory(ppid,p_init);
806 continue;
807 }else if( rdm<=ratio) {
808 break;
809 }
810 }
811 }
812
813 if(m_truthFile!="" && m_truthPart!=""){
814 if(EvtPDL::name(part->getId())==m_truthPart ){
815 int ndaug = part->getNDaug();
816 for(int id=0;id<ndaug;id++){
817 EvtParticle* sonp=part->getDaug(id);
818 EvtVector4R son=sonp->getP4();
819 truth<<son.get(0)<<" "<<son.get(1)<<" "<<son.get(2)<<" "<<son.get(3)<<std::endl;
820 }
821 }
822 }
823//--------- For superbody3
824 bool accept;
825 if(m_FDPparticle==""){FDP_part=part;}
826 if(!(m_SB3File=="" || m_SB3HT=="")){
827 accept=SuperBody3decay_judge( FDP_part);
828 }
829 if(!(m_SB3File=="" || m_SB3HT=="") && !accept) {
830 delete hepMCpart;
831 part->deleteTree();
832 goto loop_to_select_eventsB;
833 }else if(!(m_SB3File=="" || m_SB3HT=="") && parentFlag && accept){
834 countChannel(FDP_part);
835 }
836
837 int Nchannel=part->getChannel();
838
839 if(m_statDecays && parentFlag) countChannel(part);
840 //------------- dump the decay tree to the event model
841 // int Nchannel=part->getChannel();
842 // if(Nchannel>=0) part->dumpTree();
843
844 if(parentFlag){
845 EvtDecayTag theTag(part);
846 unsigned int modeTag, multiplicityTag;
847 modeTag = theTag.getModeTag();
848 if( getModel(part) == "OPENCHARM"|| getModel(part) == "LUNDCHARM" && m_tagLundModel ){
849 multiplicityTag = decayType(part);
850 //debugging
851 // std::cout<<"Opencharm Mode: "<<decayType(part)<<std::endl;
852 }else if(_RvalueTag){
853 multiplicityTag = decayType(part);
854 } else {
855 multiplicityTag = theTag.getMultTag() + decayType(part);
856 }
857 if(eventHeader == 0) std::cout << "event header unavailable: " << std::endl;
858 if (eventHeader != 0) {
859 eventHeader->setFlag1(modeTag);
860 eventHeader->setFlag2(multiplicityTag);
861 }
862 //debugg
863 // std::cout<<modeTag<<" "<<multiplicityTag<<std::endl;
864 if(m_NtupleFile)m_flag1 = modeTag;
865 //std::cout<<"Flag1 "<<modeTag<<std::endl;
866 }
867
868 if(!(m_DecayRec=="")) {mydecayrec=part->writeTreeRec(m_DecayRec);std::cout<<std::endl;};
869 // ----------- pingrg 2008-03-25 ---------
870
871 if ( log.level() <= MSG::DEBUG ) part->printTree() ;
872 log << MSG::DEBUG << "Converting particles to HepMC" << endreq;
873
874 makeHepMC(part, hepMCevt, hepMCpart);
875
876 if(part->getNDaug()!=0) hepMCpart->set_status(999);
877
878 //-------------
879 if(EvtPDL::getId(m_printfs)==part->getId()) {
880 int myd=part->getNDaug();
881 std::cout<<"_FinalState "<<myd;
882 for(int ii=0;ii<myd;ii++) {
883 std::cout<<" "<<EvtPDL::name(part->getDaug(ii)->getId());
884 }
885 std::cout<<std::endl;
886 }
887
888 AllTrk_index=0;
889 for( itp=hepMCevt->particles_begin(); itp!=hepMCevt->particles_end(); ++itp )
890 {
891 if (!(*itp)->end_vertex() && ( (*itp)->status()==899 || (*itp)->status()==999 ))
892 (*itp)->set_status(1);
893
894 HepMC::GenParticle* hepMCpart=*itp;
895 int stat = hepMCpart->status();
896 int id = hepMCpart->pdg_id();
897 if(abs(id)==411 ||abs(id)==413 ) { (*itp)->set_status(999);stat=999;}
898
899 if ( stat != 1 ) continue;
900 if( id == 22 ) {N_gamma ++;fst[0]++;}
901 if( id == -22 ){N_gammaFSR ++;}
902 if(id == 11 ) {fst[1]++;}
903 else if(id == -11) {fst[2]++;}
904 else if(id == 13 ) {fst[3]++;}
905 else if(id ==-13) {fst[4]++;}
906 else if(id==211) {fst[5]++;}
907 else if(id==-211) {fst[6]++;}
908 else if(id== 321) {fst[7]++;}
909 else if(id ==-321) {fst[8]++;}
910 else if(id ==2212) {fst[9]++;}
911 else if(id==-2212) {fst[10]++;}
912 else if(id == 2112){fst[11]++;}
913 else if(id==-2112) {fst[12]++;}
914 if( fabs(id) == 11) {nchr_e++;}
915 if( fabs(id) == 13) {nchr_mu++;}
916 if( fabs(id) == 211) {nchr_pi++;}
917 if( fabs(id) == 321) {nchr_k++;}
918 if( fabs(id) == 2212) {nchr_p++;}
919
920 //for Nuple
921 double en =(hepMCpart->momentum()).e();
922 double pex=(hepMCpart->momentum()).px();
923 double pey=(hepMCpart->momentum()).py();
924 double pez=(hepMCpart->momentum()).pz();
925
926 if( m_NtupleFile==true && id!=0){
927 m_Trk_index[AllTrk_index]=id;
928 m_px_trk[AllTrk_index]=pex;
929 m_py_trk[AllTrk_index]=pey;
930 m_pz_trk[AllTrk_index]=pez;
931 m_en_trk[AllTrk_index]=en ; /*
932 std::cout<<"trk# " <<AllTrk_index
933 <<" PID:" <<id
934 <<" px: " <<pex
935 <<" py: " <<pey
936 <<" pz: " <<pez
937 <<" E: " <<en<<std::endl; */
938 AllTrk_index++;
939 }
940
941 }
942
943 itp=hepMCevt->particles_begin(); // parent decaying particle status set
944 (*itp)->set_status(2);
945
946 //nchr = nchr_e + nchr_mu + nchr_pi + nchr_k +nchr_p;
947 nchr = nchr_pi + nchr_k +nchr_p;
948 if(m_Ncharge == true ) {std::cout<<"Multiplicity: TotoalCharge_e_mu_pi_K_p_gamma_gammaFSR: "
949 <<nchr<<" "
950 <<nchr_e<<" "
951 << nchr_mu <<" "
952 << nchr_pi <<" "
953 << nchr_k <<" "
954 <<nchr_p<<" "
955 <<N_gamma<<" "
956 <<N_gammaFSR<<endl;}
957 if(m_Ncharge == true ){std::cout<<"FinalStates: "
958 <<fst[0]<<" "
959 <<fst[1]<<" "
960 <<fst[2]<<" "
961 <<fst[3]<<" "
962 <<fst[4]<<" "
963 <<fst[5]<<" "
964 <<fst[6]<<" "
965 <<fst[7]<<" "
966 <<fst[8]<<" "
967 <<fst[9]<<" "
968 <<fst[10]<<" "
969 <<fst[11]<<" "
970 <<fst[12]<<std::endl;}
971
972 //if(nchr==3) part->printTree();
973 if( m_NtupleFile ){
974 m_nchr = nchr;
975 m_nchr_e = nchr_e;
976 m_nchr_mu = nchr_mu;
977 m_nchr_pi = nchr_pi;
978 m_nchr_k = nchr_k;
979 m_nchr_p = nchr_p;
980 m_gamma=N_gamma;
981 m_gammaFSR=N_gammaFSR;
982 TotNumTrk = AllTrk_index; //nchr + N_gamma + N_gammaFSR;
983 }
984
985 //debug
986 if(m_ReadTruth.size())ReadTruth(part,m_ReadTruth);
987
988
989 //if(part->getId()==EvtPDL::getId("psi(4260)")) std::cout<<"mydebug "<<part->getP4().mass()<<std::endl;
990
991 part->deleteTree();
992 return StatusCode::SUCCESS;
993}
994
995//************ end of callBesEvtGen *********************
996StatusCode EvtDecay::finalize()
997{
998
999 if(EvtCalHelAmp::nevt>0){
1000 double H2=EvtCalHelAmp::_H2;
1001 double H1=EvtCalHelAmp::_H1;
1002 double H0=EvtCalHelAmp::_H0;
1003 double H2err=EvtCalHelAmp::_H2err;
1004 double H1err=EvtCalHelAmp::_H1err;
1005 double H0err=EvtCalHelAmp::_H0err;
1006 int nd = EvtCalHelAmp::nevt;
1007 std::cout<<"Calculated helicity amplitude sqaured are:"<<std::endl;
1008 std::cout<<"|H_{2}|^2=|H_{-2}|^2= "<<H2/nd<<" +/- "<<sqrt(H2err/nd)<<endl;
1009 std::cout<<"|H_{1}|^2=|H_{-1}|^2= "<<H1/nd<<" +/- "<<sqrt(H1err/nd)<<endl;
1010 //std::cout<<"|H_{0}|^2= "<<H0/nd <<" +/- "<<sqrt(H0err/nd)<<endl;
1011 }
1012 MsgStream log(messageService(), name());
1013 truth.close();
1014 log << MSG::INFO << "EvtDecay finalized" << endreq;
1015 fstream Forfile;
1016 Forfile.open("fort.16",ios::in);
1017 char delcmd[300];
1018 strcpy(delcmd,"rm -f ");
1019 strcat(delcmd,"fort.16");
1020 // if(Forfile)system(delcmd);
1021
1022 fstream Forfile2;
1023 Forfile2.open("fort.22",ios::in);
1024 strcpy(delcmd,"rm -f ");
1025 strcat(delcmd,"fort.22");
1026 // if(Forfile2)system(delcmd);
1027
1028 // To get the branching fraction of the specified mode in Lundcharm Model
1029 /*
1030 EvtLundCharm lundcharmEvt;
1031 int nevt=lundcharmEvt.getTotalEvt();
1032 std::cout<<"The total number of the specified mode is "<<nevt<<std::endl;
1033 */
1034 int totalEvt=0;
1035 if(!(m_SB3File=="" || m_SB3HT=="")){
1036 for(int i=0;i<500;i++){totalEvt=totalEvt+br[i];}
1037 for(int i=0;i<500;i++){
1038 double bi=1.0*br[i]/totalEvt;
1039 std::cout<<"Branching fraction of channel "<<i<<" is: "<<bi<<std::endl;
1040 if(br[i]==0) break;
1041 }
1042 }
1043
1044 if(m_statDecays){
1045 int totalEvt=0;
1046 for(int i=0;i<=totalChannels;i++){totalEvt=totalEvt+br[i];}
1047 std::cout<<"==================Statistical first chain decay ==============================="<<std::endl;
1048 std::cout<<"i-th channel Events Generated Branching fraction generated "<<std::endl;
1049 for(int i=0;i<=totalChannels;i++){
1050 std::cout<<"| "<<i<<" "<<br[i]<<" "<<br[i]*1.0/totalEvt<<std::endl;
1051
1052 }
1053 std::cout<<"--------------------------------------------------------------------------------"<<std::endl;
1054 std::cout<<" sum: "<<totalEvt<<" "<<"1"<<std::endl;
1055 std::cout<<"================== End of statistical first chain decay ========================"<<std::endl;
1056 }
1057
1058 return StatusCode::SUCCESS;
1059}
1060
1061
1062StatusCode EvtDecay::makeHepMC(EvtParticle* part, HepMC::GenEvent* hEvt,
1063 HepMC::GenParticle* hPart)
1064{
1065 MsgStream log(messageService(), name());
1066
1067 if(part->getNDaug()!=0)
1068 {
1069 double ct=(part->getDaug(0)->get4Pos()).get(0);
1070 double x=(part->getDaug(0)->get4Pos()).get(1);
1071 double y=(part->getDaug(0)->get4Pos()).get(2);
1072 double z=(part->getDaug(0)->get4Pos()).get(3);
1073
1074 HepMC::GenVertex* end_vtx = new HepMC::GenVertex(HepLorentzVector(x,y,z,ct));
1075 hEvt->add_vertex( end_vtx );
1076 end_vtx->add_particle_in(hPart);
1077
1078 int ndaug=part->getNDaug();
1079
1080 for(int it=0;it<ndaug;it++)
1081 {
1082
1083 double e=(part->getDaug(it)->getP4Lab()).get(0);
1084 double px=(part->getDaug(it)->getP4Lab()).get(1);
1085 double py=(part->getDaug(it)->getP4Lab()).get(2);
1086 double pz=(part->getDaug(it)->getP4Lab()).get(3);
1087 int id=EvtPDL::getStdHep(part->getDaug(it)->getId());
1088 int status=1;
1089
1090 if(part->getDaug(it)->getNDaug()!=0) status=999;
1091
1092 HepMC::GenParticle* prod_part = new HepMC::GenParticle(HepLorentzVector(px,py,pz,e),id,status);
1093 end_vtx->add_particle_out(prod_part);
1094 makeHepMC(part->getDaug(it),hEvt,prod_part);
1095
1096
1097
1098 if( log.level()<MSG::INFO )
1099 prod_part->print();
1100 }
1101 }
1102
1103 return StatusCode::SUCCESS;
1104}
1105
1106
1107void EvtDecay::MeVToGeV (HepMC::GenEvent* evt)
1108{
1109 for ( HepMC::GenEvent::particle_iterator p = evt->particles_begin(); p != evt->particles_end(); ++p ) {
1110 // cout << " PDG, BAR " << (*p)->pdg_id() << " " << (*p)->barcode() << endl;
1111 // (*p)->set_momentum( (*p)->momentum() / 1000. );
1112
1113 HepMC::FourVector tmpfv((*p)->momentum().x()/1000.,
1114 (*p)->momentum().y()/1000.,
1115 (*p)->momentum().z()/1000.,
1116 (*p)->momentum().t()/1000.
1117 );
1118
1119 (*p)->set_momentum( tmpfv );
1120 }
1121}
1122
1123void EvtDecay::GeVToMeV (HepMC::GenEvent* evt)
1124{
1125 for ( HepMC::GenEvent::particle_iterator p = evt->particles_begin(); p != evt->particles_end(); ++p ) {
1126 // cout << " PDG, BAR " << (*p)->pdg_id() << " " << (*p)->barcode() << endl;
1127 // (*p)->set_momentum( (*p)->momentum() * 1000. );
1128 HepMC::FourVector tmpfv((*p)->momentum().x()*1000.,
1129 (*p)->momentum().y()*1000.,
1130 (*p)->momentum().z()*1000.,
1131 (*p)->momentum().t()*1000.
1132 );
1133
1134 (*p)->set_momentum( tmpfv );
1135 }
1136}
1137
1138
1139double EvtDecay::CalAmpsMax( EvtParticle* part )
1140{
1141 double amps_max=0;
1142 for(int i=0;i<100000;i++){
1143 m_Gen->generateDecay(part);
1144 double amps=CalAmpsMDIY(part);
1145 if(amps > amps_max) amps_max=amps;
1146 }
1147 amps_max=amps_max*1.05;
1148 m_ampscalflag = false;
1149 return amps_max;
1150}
1151
1152// EvtBesRandom class implementation, basically BesRandomSvc -> EvtGen random engine interface
1153
1154EvtBesRandom::EvtBesRandom(HepRandomEngine* engine)
1155{
1156 m_engine = engine;
1157 m_engine->showStatus();
1158}
1159
1160EvtBesRandom::~EvtBesRandom()
1161{}
1162
1163double EvtBesRandom::random()
1164{
1165 //double rdm=EvtRandom::Flat(0.0, 1.0);;
1166 double rdm= RandFlat::shoot(m_engine);
1167 // std::cout<<"rdm= "<<rdm<<std::endl;
1168 return rdm;
1169 //return RandFlat::shoot(m_engine);
1170}
1171
1172
1173void EvtDecay::SuperBody3decay_make(EvtId ppid, EvtVector4R p_init ){
1174 double xmass2,ymass2;
1175 bool accept;
1176 EvtVector4R pd1,pd2,pd3;
1177
1178 //---- find out the pdg codes of final state and count the identical particles
1179 FinalState_make( ppid, p_init);
1180
1181 //begin to loop with events
1182 for(int i=0;i<100000;i++){
1183 regen:
1184 EvtParticle* part=EvtParticleFactory::particleFactory(ppid,p_init);
1185 HepMC::GenParticle* hepMCpart=new GenParticle(HepLorentzVector(0,0,0,parentMass),parentPDGcode,3); //this line make the decay to select different channel
1186
1187 if(part->getSpinType() == EvtSpinType::VECTOR) {
1188 part->setVectorSpinDensity();
1189 }
1190 m_Gen->generateDecay(part);
1191
1192
1193 FinalState_sort(part);
1194 pd1=son0;
1195 pd2=son1;
1196 pd3=son2;
1197
1198 // std::cout<<"multi, pdg, pi= "<<multi<<" "<<pdg0<<" "<<pdg1<<" "<<pdg2<<" "<<son0<<son1<<son2<<std::endl;
1199
1200 xmass2=(pd1+pd2).mass2(); // Dalitz plot m12^2 ~ m13^2
1201 ymass2=(pd1+pd3).mass2();
1202 double xlow=SuperBody3decay.getXlow();
1203 double xup=SuperBody3decay.getXup();
1204 double ylow=SuperBody3decay.getYlow();
1205 double yup=SuperBody3decay.getYup();
1206 if(xmass2<xlow || xmass2>xup || ymass2<ylow || ymass2>yup) { part->deleteTree();goto regen;}
1207 SuperBody3decay.HFill(xmass2,ymass2);
1208
1209 if( m_NtupleFile ){
1210 m_ich=part->getChannel();
1211 m_m1=pd1.mass();
1212 m_m2=pd2.mass();
1213 m_m3=pd3.mass();
1214 m_m12=(pd1+pd2).mass();
1215 m_m23=(pd2+pd3).mass();
1216 m_m13=(pd1+pd3).mass();
1217 m_cos1=pd1.get(3)/pd1.d3mag();
1218 m_cos2=pd2.get(3)/pd2.d3mag();
1219 m_cos3=pd3.get(3)/pd3.d3mag();
1220 mass_tuple->write();
1221 }
1222 double m1=pd1.mass(); double m2=pd2.mass();double m3=pd3.mass();
1223 double mj=(pd2+pd3).mass();
1224 // std::cout<<"mass 1 2 3 chicj <<=="<<m1<<" "<<m2<<" "<<m3<<" "<<mj<<std::endl;
1225 // delete hepMCpart;
1226 }
1227
1228 SuperBody3decay.HReweight(); //reweighting Dalitz plot
1229 SuperBody3decay.setZmax(); // find out the maximum value after reweighting
1230 first=false;
1231}
1232
1233bool EvtDecay::SuperBody3decay_judge(EvtParticle* part){
1234 double xmass2,ymass2;
1235 bool accept;
1236 EvtVector4R pd1,pd2,pd3;
1237
1238
1239 FinalState_sort( part);
1240
1241 pd1=son0;
1242 pd2=son1;
1243 pd3=son2;
1244
1245
1246 xmass2=(pd1+pd2).mass2(); // Dalitz plot m12^2 ~ m13^2
1247 ymass2=(pd1+pd3).mass2();
1248 accept=SuperBody3decay.AR(xmass2,ymass2);
1249
1250 //debugging
1251 if(accept && m_NtupleFile) {
1252 _ich=part->getChannel();
1253 _m1=pd1.mass();
1254 _m2=pd2.mass();
1255 _m3=pd3.mass();
1256 _m12=(pd1+pd2).mass();
1257 _m23=(pd2+pd3).mass();
1258 _m13=(pd1+pd3).mass();
1259 _cos1=pd1.get(3)/pd1.d3mag();
1260 _cos2=pd2.get(3)/pd2.d3mag();
1261 _cos3=pd3.get(3)/pd3.d3mag();
1262 massgen_tuple->write();
1263 }
1264 // std::cout<<"mass 1 2 3 chicj>> "<<_m1<<" "<<_m2<<" "<<_m3<<" "<<_m23<<std::endl;
1265
1266 return accept;
1267}
1268
1269
1270void EvtDecay:: FinalState_make(EvtId ppid, EvtVector4R p_init){
1271 // get the final state pdg cond and count the identicle particle
1272 multi=1;
1273 identical_flag=true;
1274
1275 for(int i=1;i<10000;i++){ //sigle out the final state
1276 EvtParticle* part=EvtParticleFactory::particleFactory(ppid,p_init);
1277 HepMC::GenParticle* hepMCpart=new GenParticle(HepLorentzVector(0,0,0,parentMass),parentPDGcode,3);
1278
1279 m_Gen->generateDecay(part);
1280 int nson=part->getNDaug();
1281
1282 if(nson == 2) {continue;} else
1283 if(nson ==3){
1284 EvtId xid0,xid1,xid2;
1285 xid0=part->getDaug(0)->getId();
1286 xid1=part->getDaug(1)->getId();
1287 xid2=part->getDaug(2)->getId();
1288
1289//-- pdg code
1290 pdg0=EvtPDL::getStdHep(xid0);
1291 pdg1=EvtPDL::getStdHep(xid1);
1292 pdg2=EvtPDL::getStdHep(xid2);
1293
1294 if(pdg0==pdg1 && pdg1==pdg2) {multi=3;}
1295 else if(pdg0==pdg1 ){multi=2;} // two identical particle list as 0,1 index
1296 else if(pdg0==pdg2 ){multi=2;pdg=pdg1; pdg1=pdg2; pdg2=pdg;}
1297 else if(pdg1==pdg2) {multi=2;pdg=pdg0; pdg0=pdg1; pdg1=pdg2; pdg2=pdg;}
1298 else {multi=1;}
1299 break;
1300 } else{
1301 std::cout<<"No direct three body decay channel found in decay card, I stop run"<<std::endl;
1302 ::abort();
1303 }
1304 }
1305 // std::cout<<"pdg_i "<<pdg0<<" "<<pdg1<<" "<<pdg2<<std::endl;
1306 // std::cout<<"identical particle "<<multi<<std::endl;
1307}
1308
1309void EvtDecay::FinalState_sort(EvtParticle* part){
1310 // sort the momentum to aon0, son1, son2
1311 EvtVector4R pd0,pd1,pd2;
1312 EvtId xid0,xid1,xid2;
1313 int id0,id1,id2;
1314
1315 int nson=part->getNDaug();
1316 if(nson==2){
1317 pd0=part->getDaug(0)->getP4Lab();
1318 pd1=part->getDaug(1)->getDaug(0)->getP4Lab();
1319 pd2=part->getDaug(1)->getDaug(1)->getP4Lab();
1320
1321 xid0=part->getDaug(0)->getId();
1322 xid1=part->getDaug(1)->getDaug(0)->getId();
1323 xid2=part->getDaug(1)->getDaug(1)->getId();
1324
1325 } else if(nson==3){
1326 pd0=part->getDaug(0)->getP4Lab();
1327 pd1=part->getDaug(1)->getP4Lab();
1328 pd2=part->getDaug(2)->getP4Lab();
1329
1330 xid0=part->getDaug(0)->getId();
1331 xid1=part->getDaug(1)->getId();
1332 xid2=part->getDaug(2)->getId();
1333 }
1334 //-- pdg code
1335 id0=EvtPDL::getStdHep(xid0);
1336 id1=EvtPDL::getStdHep(xid1);
1337 id2=EvtPDL::getStdHep(xid2);
1338
1339 // std::cout<<"pid before match: "<<id0<<" "<<id1<<" "<<id2<<std::endl;
1340 //-- assign sons momentum
1341 if(multi==1){
1342 assign_momentum(id0,pd0);
1343 assign_momentum(id1,pd1);
1344 assign_momentum(id2,pd2);
1345 } else if(multi==2){
1346 assign_momentum2(id0,pd0);
1347 assign_momentum2(id1,pd1);
1348 assign_momentum2(id2,pd2);
1349 if(son0.get(0)>son1.get(0)){son=son0;son0=son1;son1=son;} // change into E_0 < E_1
1350 } else if(multi==3){ // sort sons according to energy E_0 < E_1 < E_2
1351 double En0=pd0.get(0);
1352 double En1=pd1.get(0);
1353 double En2=pd2.get(0);
1354 if(En0 < En1 && En1 < En2) {son0=pd0;son1=pd1;son2=pd2;}
1355 if(En0 < En2 && En2 < En1) {son0=pd0;son1=pd2;son2=pd1;}
1356 if(En1 < En0 && En0 < En2) {son0=pd1;son1=pd0;son2=pd2;}
1357 if(En1 < En2 && En2 < En0) {son0=pd1;son1=pd2;son2=pd0;}
1358 if(En2 < En0 && En0 < En1) {son0=pd2;son1=pd0;son2=pd1;}
1359 if(En2 < En1 && En1 < En0) {son0=pd2;son1=pd1;son2=pd0;}
1360 }
1361
1362}
1363
1364
1365void EvtDecay::assign_momentum(int id0, EvtVector4R pd0){
1366 if(id0==pdg0){son0=pd0;}
1367 else if(id0==pdg1){son1=pd0;}
1368 else if(id0==pdg2){son2=pd0;}
1369 else {std::cout<<"PID= "<<id0<<" not machted, I stop"<<std::endl; ::abort();}
1370}
1371
1372void EvtDecay::assign_momentum2(int id0, EvtVector4R pd0){ // for two identicle particle case
1373 if(id0==pdg0 && identical_flag){son0=pd0;identical_flag=false;}
1374 else if(id0==pdg1){son1=pd0;identical_flag=true;}
1375 else if(id0==pdg2){son2=pd0;}
1376 else {std::cout<<"PID not machted, I stop"<<std::endl; ::abort();}
1377}
1378
1379void EvtDecay::findPart(EvtParticle* part){
1380 FDP_id=EvtPDL::getId(m_FDPparticle);
1381 int FDP_pdg=EvtPDL::getStdHep(FDP_id);
1382
1383 m_Gen->generateDecay(part);
1384 int dn=part->getNDaug();
1385
1386 if(dn!=0){
1387 for(int i=0;i<dn;i++){
1388 EvtParticle* part1=part->getDaug(i);
1389 EvtId sonid=part1->getId();
1390 int son_pdg = EvtPDL::getStdHep(sonid);
1391 if(son_pdg == FDP_pdg){
1392 FDP_init=part1->getP4Lab();
1393 FDP_part=part1;
1394 return;
1395 } else{
1396 findPart(part1);
1397 }
1398 }
1399 } //if (dn block
1400
1401}
1402
1403void EvtDecay::countChannel(EvtParticle* part){
1404 int ich=part->getChannel();
1405
1406 //debuging
1407 //if(getModel(part,ich)=="OPENCHARM" &&EvtPDL::name( part->getId() )=="psi(4260)") ich = part->getGeneratorFlag();
1408 //std::cout<<"the channel is "<<ich<<std::endl;
1409
1410 br[ich]++;
1411 if(ich>totalChannels) totalChannels = ich;
1412 //std::cout<<"EVENT IN br_i "<<br[ich]<<std::endl;
1413}
1414
1415bool EvtDecay::isCharmonium(EvtId xid){
1416EvtId psi4415 = EvtPDL::getId(std::string("psi(4415)"));
1417EvtId psi4160 = EvtPDL::getId(std::string("psi(4160)"));
1418EvtId psi4040 = EvtPDL::getId(std::string("psi(4040)"));
1419EvtId psi3770 = EvtPDL::getId(std::string("psi(3770)"));
1420EvtId psiprim = EvtPDL::getId(std::string("psi(2S)"));
1421EvtId jpsi = EvtPDL::getId(std::string("J/psi"));
1422EvtId etac = EvtPDL::getId(std::string("eta_c"));
1423EvtId etac2s = EvtPDL::getId(std::string("eta_c(2S)"));
1424EvtId hc = EvtPDL::getId(std::string("h_c"));
1425EvtId chic0 = EvtPDL::getId(std::string("chi_c0"));
1426EvtId chic1 = EvtPDL::getId(std::string("chi_c1"));
1427EvtId chic2 = EvtPDL::getId(std::string("chi_c2"));
1428EvtId chic0p = EvtPDL::getId(std::string("chi_c0p"));
1429EvtId chic1p = EvtPDL::getId(std::string("chi_c1p"));
1430EvtId chic2p = EvtPDL::getId(std::string("chi_c1p"));
1431 std::vector<EvtId> Vid; Vid.clear();
1432 Vid.push_back(psi4415);
1433 Vid.push_back(psi4160);
1434 Vid.push_back(psi4040);
1435 Vid.push_back(psi3770);
1436 Vid.push_back(psiprim);
1437 Vid.push_back(jpsi);
1438 Vid.push_back(etac);
1439 Vid.push_back(etac2s);
1440 Vid.push_back(hc);
1441 Vid.push_back(chic0);
1442 Vid.push_back(chic1);
1443 Vid.push_back(chic2);
1444 Vid.push_back(chic0p);
1445 Vid.push_back(chic1p);
1446 Vid.push_back(chic2p);
1447
1448 bool flag=true;
1449 for(int i=0;i<Vid.size();i++){ if(xid == Vid[i]) return flag;}
1450 return false;
1451}
1452
1453
1454bool EvtDecay::isCharm(EvtId xid){
1455EvtId d0 = EvtPDL::getId(std::string("D0"));
1456EvtId d0bar = EvtPDL::getId(std::string("anti-D0"));
1457EvtId dp = EvtPDL::getId(std::string("D+"));
1458EvtId dm = EvtPDL::getId(std::string("D-"));
1459EvtId d0h = EvtPDL::getId(std::string("D0H"));
1460EvtId d0l = EvtPDL::getId(std::string("D0L"));
1461EvtId dstp = EvtPDL::getId(std::string("D*+"));
1462EvtId dstm = EvtPDL::getId(std::string("D*-"));
1463EvtId ds0 = EvtPDL::getId(std::string("D*0"));
1464EvtId ds0bar = EvtPDL::getId(std::string("anti-D*0"));
1465EvtId dsp = EvtPDL::getId(std::string("D_s+"));
1466EvtId dsm = EvtPDL::getId(std::string("D_s-"));
1467EvtId dsstp = EvtPDL::getId(std::string("D_s*+"));
1468EvtId dsstm = EvtPDL::getId(std::string("D_s*-"));
1469EvtId ds0stp = EvtPDL::getId(std::string("D_s0*+"));
1470EvtId ds0stm = EvtPDL::getId(std::string("D_s0*-"));
1471
1472 std::vector<EvtId> Vid; Vid.clear();
1473 Vid.push_back(d0);
1474 Vid.push_back(d0bar);
1475 Vid.push_back(dp);
1476 Vid.push_back(dm);
1477 Vid.push_back(d0h);
1478 Vid.push_back(d0l);
1479 Vid.push_back(dstp);
1480 Vid.push_back(dstm);
1481 Vid.push_back(ds0);
1482 Vid.push_back(ds0bar );
1483 Vid.push_back(dsp );
1484 Vid.push_back(dsm );
1485 Vid.push_back(dsstp );
1486 Vid.push_back(dsstm );
1487 Vid.push_back(ds0stp );
1488 Vid.push_back(ds0stm );
1489
1490 bool flag=true;
1491 for(int i=0;i<Vid.size();i++){ if(xid == Vid[i]) return flag;}
1492 return false;
1493}
1494
1495bool EvtDecay::isRadecay(EvtParticle *par){
1496 int ndg = par->getNDaug();
1497 for(int i=0;i<ndg;i++){
1498 EvtId xid = par->getDaug(i)->getId();
1499 if(EvtPDL::getStdHep(xid) == 22){return true;}
1500 }
1501 return false;
1502}
1503
1504int EvtDecay::decayType(EvtParticle *par){
1505 /*************************************
1506 * 0: gamma light_hadrons
1507 * 1: gamma charmonium
1508 * 2: DDbar (D0 D0bar or D+D-)
1509 * 3: ligh_hadrons
1510 * 4: others
1511 *************************************/
1512 int Nchannel=par->getChannel();
1513 int ndg = par->getNDaug();
1514 int nfsr=0;
1515
1516 if(_RvalueTag){return Nchannel;}
1517 // std::cout<<"decayType= "<<Nchannel<<std::endl;
1518
1519 std::string model = getModel(par,Nchannel);
1520 if( model == "OPENCHARM" || model == "LUNDCHARM" && m_tagLundModel){ // lundcharm model tag
1521 int ldm = par->getGeneratorFlag();
1522 // std::cout<<"LundCharmFlag = "<<ldm<<std::endl;
1523 return ldm;
1524 //return 9;
1525 }
1526
1527 for(int i=0;i<ndg;i++){ //remove the FSR photon
1528 EvtId xid =par->getDaug(i)->getId();
1529 if(EvtPDL::getStdHep(xid) == -22 ){nfsr++;}
1530 }
1531
1532 if( isRadecay(par) ){
1533 for(int i=0;i<ndg;i++){
1534 EvtId xid = par->getDaug(i)->getId();
1535 if( isCharmonium(xid) ) return 1;
1536 }
1537 return 0;
1538 } else{
1539 if(ndg-nfsr ==2 ){
1540 int ndg = par->getNDaug();
1541 EvtId xd1 = par->getDaug(0)->getId();
1542 EvtId xd2 = par->getDaug(1)->getId();
1543 if( isCharm(xd1) && isCharm(xd2) ) {return 2;} else
1544 if( !isCharmonium(xd1) && !isCharmonium(xd2) && !isCharm(xd1) && !isCharm(xd2) ){
1545 return 3;
1546 }
1547 } else{ // ndg>=3
1548 bool flag = false;
1549 for(int i=0;i<ndg;i++){
1550 EvtId xid = par->getDaug(i)->getId();
1551 if( isCharmonium(xid) ) {flag = true; break;}
1552 if( isCharm(xid) ) {flag = true; break;}
1553 } // for_i block
1554 if( !flag ){return 3;} else {return 4;}
1555 }// if_ndg block
1556 }
1557
1558}
1559
1560
1561std::string EvtDecay::getModel(EvtParticle *par, int mode){
1562 EvtDecayBase* thedecay = EvtDecayTable::gettheDecay(par->getId(),mode);
1563 return thedecay->getModelName();
1564}
1565
1566std::string EvtDecay::getModel(EvtParticle *par){
1567 int mode = par->getChannel();
1568 EvtDecayBase* thedecay = EvtDecayTable::gettheDecay(par->getId(),mode);
1569 return thedecay->getModelName();
1570}
1571
1572
1573void EvtDecay::ReadTruth(EvtParticle* part,std::vector<std::vector<string> > mylist){
1574 if(mylist.size()<2) {std::cout<<" Empty list"<<std::endl;abort();}
1575 EvtVector4R vp1;
1576 for(int k=0;k<mylist[0].size();k++){
1577 if(part->getId()==EvtPDL::getId(mylist[0][k])){
1578 if(part->getDaug(1)->getId()==EvtPDL::getId("vhdr")){part=part->getDaug(1);continue;}
1579 if(part->getDaug(0)->getId()==EvtPDL::getId("vgam")){part=part->getDaug(1);continue;}
1580 for(int i=1;i<mylist.size();i++){
1581 EvtParticle* mypart=part;
1582 for(int j=0;j<mylist[i].size();j++){
1583 mypart=mypart->getDaug(atoi(mylist[i][j].c_str()));
1584 //std::cout<<atoi(mylist[i][j].c_str());
1585 }
1586 //std::cout<<std::endl;
1587 std::cout<<EvtPDL::name(mypart->getId())<<std::endl;
1588 if(_truthAtCM){
1589 vp1=mypart->getP4();
1590 }else{
1591 vp1=mypart->getP4Lab();
1592 }
1593 if( !isNumber(vp1.get(0)) ) {part->printTree(); abort();}
1594 //std::cout<<"truth "<<vp1.get(0)<<" "<<vp1.get(1)<<" "<<vp1.get(2)<<" "<<vp1.get(3)<<std::endl;
1595 printf("truth %.12lf %.12lf %.12lf %.12lf \n",vp1.get(0),vp1.get(1),vp1.get(2),vp1.get(3) );
1596 }
1597 }
1598 }
1599}
1600
1601int EvtDecay::isNumber(double d){return (d==d);}//if d=nan, return 0, otherwise, return 1
1602
1603
1604double EvtDecay::energySpread(double mu,double sigma){
1605 //mu: mean value in Gaussian
1606 //sigma: variance in Gaussian
1607 rloop:
1608 int n=12;
1609 double ri=0;
1610 for(int i=0;i<n;i++){
1611 double pm= EvtRandom::Flat(0.,1);
1612 ri += pm;
1613 }
1614 double eta=sqrt(n*12.0)*(ri/12-0.5);
1615 double xsig= sigma*eta+mu;
1616 double mx0=EvtConExc::mlow;
1617 double mx1=EvtConExc::mup;
1618 if(xsig<mx0 || xsig>mx1) goto rloop;
1619 return xsig;
1620}
1621
1622#include "../user/Lenu.inc"
mass
double mass
Definition: CosmicGenerator.cxx:138
sigma
TTree * sigma
Definition: CurConst_dEdXCal.cxx:6
runNo
int runNo
Definition: DQA_TO_DB.cxx:12
n
const Int_t n
Definition: DataBase/tau_mode.c:65
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57
s
XmlRpcServer s
Definition: HelloServer.cpp:11
McGenEventCol
ObjectVector< McGenEvent > McGenEventCol
Definition: McGenEvent.h:39
ntupleSvc
INTupleSvc * ntupleSvc()
Definition: ServiceAccessor.h:55
key
*************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
pd0
double complex pa0 double complex pb0ij double complex pc0ijk double complex pd0
Definition: bak-BesEvtGen-00-04-08/src/phokhara/eemmg-lib/eemmg5/pjfry.h:5
DataInfoSvc::setDecayCard
void setDecayCard(string card)
Definition: DataInfoSvc.cxx:50
EvtBesRandom::EvtBesRandom
EvtBesRandom(HepRandomEngine *engine)
Definition: bak-BesEvtGen-00-04-08/src/EvtDecay.cxx:1154
EvtDecayTable::gettheDecay
static EvtDecayBase * gettheDecay(EvtId parent, int imode)
Definition: bak-BesEvtGen-00-04-08/src/EvtGen/EvtGenBase/EvtDecayTable.cc:57
EvtDecay::EvtDecay
EvtDecay(const string &name, ISvcLocator *pSvcLocator)
Definition: bak-BesEvtGen-00-04-08/src/EvtDecay.cxx:83
EvtDecay::energySpread
double energySpread(double mu, double sigma)
Definition: bak-BesEvtGen-00-04-08/src/EvtDecay.cxx:1604
EvtGen::readUDecay
void readUDecay(const char *const udecay_name)
Definition: bak-BesEvtGen-00-04-08/src/EvtGen/EvtGen.cc:126
EvtGen::generateDecay
void generateDecay(int stdhepid, EvtVector4R P, EvtVector4R D, EvtStdHep *evtStdHep, EvtSpinDensity *spinDensity=0)
Definition: bak-BesEvtGen-00-04-08/src/EvtGen/EvtGen.cc:152
EvtGlobalSet::dVV
static std::vector< std::vector< double > > dVV
Definition: bak-BesEvtGen-00-04-08/src/EvtGen/EvtGenModels/EvtGlobalSet.hh:21
EvtGlobalSet::iVV
static std::vector< std::vector< int > > iVV
Definition: bak-BesEvtGen-00-04-08/src/EvtGen/EvtGenModels/EvtGlobalSet.hh:22
EvtHis2F::AR
bool AR(double xmass2, double ymass2)
Definition: bak-BesEvtGen-00-04-08/src/EvtGen/EvtGenBase/EvtHis2F.cc:205
EvtHis2F::HFill
void HFill(double xmass2, double ymass2)
Definition: bak-BesEvtGen-00-04-08/src/EvtGen/EvtGenBase/EvtHis2F.cc:129
EvtPDL::reSetMass
static void reSetMass(EvtId i, double mass)
Definition: bak-BesEvtGen-00-04-08/src/EvtGen/EvtGenBase/EvtPDL.hh:68
EvtPDL::getId
static EvtId getId(const std::string &name)
Definition: bak-BesEvtGen-00-04-08/src/EvtGen/EvtGenBase/EvtPDL.cc:287
EvtParticleFactory::particleFactory
static EvtParticle * particleFactory(EvtSpinType::spintype spinType)
Definition: bak-BesEvtGen-00-04-08/src/EvtGen/EvtGenBase/EvtParticleFactory.cc:46
EvtParticle::setPolarizedSpinDensity
void setPolarizedSpinDensity(double r00, double r11, double r22)
Definition: bak-BesEvtGen-00-04-08/src/EvtGen/EvtGenBase/EvtParticle.cc:157
IBesRndmGenSvc::setGenseed
virtual void setGenseed(long)=0
IBesRndmGenSvc::GetEngine
virtual CLHEP::HepRandomEngine * GetEngine(const std::string &StreamName)=0
Interface to the CLHEP engine.
y
double y[1000]
Definition: draw_charge_space_coarse_II.cxx:8
rb::hc
const double hc
Definition: TConstant.h:41