BOSS 7.0.3
BESIII Offline Software System
Loading...
Searching...
No Matches
TwoGamma.cxx
Go to the documentation of this file.
1#include "GaudiKernel/MsgStream.h"
2#include "GaudiKernel/AlgFactory.h"
3#include "GaudiKernel/ISvcLocator.h"
4#include "GaudiKernel/SmartDataPtr.h"
5#include "GaudiKernel/IDataProviderSvc.h"
6#include "GaudiKernel/PropertyMgr.h"
7#include "EventModel/EventModel.h"
8#include "EventModel/Event.h"
9#include "EvtRecEvent/EvtRecEvent.h"
10#include "EvtRecEvent/EvtRecTrack.h"
11
12#include "DstEvent/TofHitStatus.h"
13#include "EventModel/EventHeader.h"
14//#include "EvTimeEvent/RecEsTime.h"
15
16#include "TMath.h"
17#include "GaudiKernel/INTupleSvc.h"
18#include "GaudiKernel/NTuple.h"
19#include "GaudiKernel/Bootstrap.h"
20#include "GaudiKernel/IHistogramSvc.h"
21#include "CLHEP/Vector/ThreeVector.h"
22using CLHEP::Hep3Vector;
23#include "CLHEP/Vector/LorentzVector.h"
24using CLHEP::HepLorentzVector;
25#include "CLHEP/Vector/TwoVector.h"
26using CLHEP::Hep2Vector;
27#include "CLHEP/Geometry/Point3D.h"
28#ifndef ENABLE_BACKWARDS_COMPATIBILITY
29typedef HepGeom::Point3D<double> HepPoint3D;
30#endif
31#include "VertexFit/KinematicFit.h"
32#include "VertexFit/VertexFit.h"
33#include "VertexFit/SecondVertexFit.h"
34#include "ParticleID/ParticleID.h"
35#include "TwoGammaAlg/TwoGamma.h"
36#include <vector>
37#include <fstream>
38typedef std::vector<int> Vint;
39typedef std::vector<HepLorentzVector> Vp4;
40const double ecms = 3.686;
41const double velc = 299.792458;
42const double xmass[5] = {0.000511, 0.105658, 0.139570, 0.493677, 0.938272};
43const double pai=3.1415926;
44
45////////////////////////////////////////////////////////////////////
46TwoGamma::TwoGamma(const std::string& name, ISvcLocator* pSvcLocator) :
47 Algorithm(name, pSvcLocator)
48{
49 // m_tuple1 = 0;
50 for(int i = 0; i < 10; i++)
51 {
52 m_pass[i] = 0;
53 }
54 m_event = 0;
55
56 Ndata1 = 0;
57 Ndata2 = 0;
58 m_runNo = 0;
59
60 //Declare the properties
61 declareProperty("CmsEnergy", m_ecms = 3.686);
62
63 declareProperty("max1", m_max1 = 1.2);
64 declareProperty("max2", m_max2 = 0.8);
65 declareProperty("costheta", m_costheta = 0.8);
66 declareProperty("dphi1", m_dphi1 = 2.5);
67 declareProperty("dphi2", m_dphi2 = 5);
68 declareProperty("eff", m_eff = 0.07);
69 declareProperty("sec", m_sec = 184.1);
70}
71
72// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
73StatusCode TwoGamma::initialize()
74{
75 MsgStream log(msgSvc(), name());
76
77 log << MSG::INFO << "in initialize()" << endmsg;
78
79 StatusCode status;
80
81 NTuplePtr nt1(ntupleSvc(), "FILE1/DiGam");
82 if ( nt1 ) m_tuple1 = nt1;
83 else
84 {
85 m_tuple1 = ntupleSvc()->book ("FILE1/DiGam", CLID_ColumnWiseTuple, "DiGam N-Tuple signal");
86 if ( m_tuple1 )
87 {
88 status = m_tuple1->addItem ("run", m_run );
89 status = m_tuple1->addItem ("rec", m_rec );
90 status = m_tuple1->addItem ("time", m_time );
91
92 status = m_tuple1->addItem ("ngood", m_ngood);
93 status = m_tuple1->addItem ("nchrg", m_nchrg);
94
95 status = m_tuple1->addItem ("Cms_e1", m_e1);
96 status = m_tuple1->addItem ("Cms_e2", m_e2);
97 status = m_tuple1->addItem ("Cms_e", m_e);
98 status = m_tuple1->addItem ("Cms_costheta1", m_costheta1);
99 status = m_tuple1->addItem ("Cms_costheta2", m_costheta2);
100 status = m_tuple1->addItem ("Cms_dltphi", m_dltphi);
101 status = m_tuple1->addItem ("Cms_dltphi_1", m_dltphi_1);
102 status = m_tuple1->addItem ("Cms_dlttheta", m_dlttheta);
103 status = m_tuple1->addItem ("Cms_phi1", m_phi1);
104 status = m_tuple1->addItem ("Cms_phi2", m_phi2);
105
106 status = m_tuple1->addItem ("Lab_e1", m_e1_lab);
107 status = m_tuple1->addItem ("Lab_e2", m_e2_lab);
108 status = m_tuple1->addItem ("Lab_e", m_e_lab);
109 status = m_tuple1->addItem ("Lab_costheta1", m_costheta1_lab);
110 status = m_tuple1->addItem ("Lab_costheta2", m_costheta2_lab);
111 status = m_tuple1->addItem ("Lab_dltphi", m_dltphi_lab);
112 status = m_tuple1->addItem ("Lab_dlttheta", m_dlttheta_lab);
113 status = m_tuple1->addItem ("Lab_phi1", m_phi1_lab);
114 status = m_tuple1->addItem ("Lab_phi2", m_phi2_lab);
115
116 status = m_tuple1->addItem ("xBoost", m_xBoost);
117 status = m_tuple1->addItem ("yBoost", m_yBoost);
118 status = m_tuple1->addItem ("zBoost", m_zBoost);
119 }
120 else
121 {
122 log << MSG::ERROR << "Cannot book N-tuple1:"<<long(m_tuple1)<<endmsg;
123 return StatusCode::FAILURE;
124 }
125 }
126 Ndata1 = 0;
127 Ndata2 = 0;
128 m_runNo = 0;
129
130 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
131 return StatusCode::SUCCESS;
132}
133
134//*********************************************************************************************************************
135StatusCode TwoGamma::execute()
136{
137 StatusCode sc=StatusCode::SUCCESS;
138 m_event++;
139
140 MsgStream log(msgSvc(),name());
141 log<<MSG::INFO<<"in execute()"<<endreq;
142 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
143 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(),EventModel::EvtRec::EvtRecEvent);
144 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(),EventModel::EvtRec::EvtRecTrackCol);
145
146 m_run = eventHeader->runNumber();
147 m_rec = eventHeader->eventNumber();
148 int runNo = m_run;
149 int event = m_rec;
150 int time = eventHeader->time();
151 if(m_event > 1 && runNo != m_runNo) return sc;
152 m_runNo = runNo;
153 m_time = time;
154
155 if(m_rec%1000==0)cout<<"Run "<<m_run<<" Event "<<m_rec<<endl;
156
157 HepLorentzVector p4psip(0.011*m_ecms, 0.0, 0.005, m_ecms);
158 double psipBeta = (p4psip.vect()).mag()/(p4psip.e());
159
160 m_pass[0]+=1;
161
162 int NCharge = evtRecEvent->totalCharged();
163 if(NCharge != 0)return sc;
164
165 m_pass[1]+=1;
166
167 double Emax1=0.0;
168 double Emax2=0.0;
169 int Imax[2];
170
171 if(((evtRecEvent->totalTracks() - evtRecEvent->totalCharged() )<2) || ((evtRecEvent->totalTracks() - evtRecEvent->totalCharged() )>15)) return sc;
172 m_pass[2]+=1;
173
174 HepLorentzVector p4Gam_1;
175 HepLorentzVector p4Gam_2;
176
177 for(int i=evtRecEvent->totalCharged(); i<evtRecEvent->totalTracks(); i++)
178 {
179 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
180 if(!(*itTrk)->isEmcShowerValid()) continue;
181 RecEmcShower *emcTrk = (*itTrk)->emcShower();
182
183 HepLorentzVector p4Gam;
184 double Phi = emcTrk->phi();
185 double Theta = emcTrk->theta();
186 double Ener=emcTrk->energy();
187
188 p4Gam.setPx(Ener * sin(Theta) * cos(Phi));
189 p4Gam.setPy(Ener * sin(Theta) * sin(Phi));
190 p4Gam.setPz(Ener * cos(Theta));
191 p4Gam.setE(Ener);
192 p4Gam.boost(-0.011, 0 ,-0.005 * 1.0 / 3.686);
193
194 if(Ener>Emax2)
195 {
196 Emax2=Ener;
197 Imax[1]=i;
198 p4Gam_2 = p4Gam;
199 }
200 if(Ener>Emax1)
201 {
202 Emax2=Emax1;
203 p4Gam_2 = p4Gam_1;
204 Imax[1]=Imax[0];
205 Emax1=Ener;
206 p4Gam_1 = p4Gam;
207 Imax[0]=i;
208 }
209 }
210 m_e1_lab = Emax1;
211 m_e2_lab = Emax2;
212 m_e_lab = Emax1 + Emax2;
213 log << MSG::INFO << "Emax1 = " << Emax1 <<"Emax2= "<<Emax2<< endreq;
214 if(Emax1 < m_max1 || Emax2 < m_max2) return sc;
215 m_pass[3]+=1;
216
217 //P4 in Lab
218 double emcphi[2],emctht[2];
219 for(int i = 0;i < 2; i++)
220 {
221 if(i>=evtRecTrkCol->size()) break;
222 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + Imax[i];
223 if(!(*itTrk)->isEmcShowerValid()) continue;
224 RecEmcShower *emcTrk = (*itTrk)->emcShower();
225 emcphi[i]=emcTrk->phi();
226 emctht[i]=emcTrk->theta();
227 }
228 double dltphi_lab = (fabs(emcphi[0] - emcphi[1]) - pai) * 180.0 / pai;
229 double dlttheta_lab = (fabs(emctht[0] + emctht[1]) - pai) * 180.0 / pai;
230 m_costheta1_lab = cos(emctht[0]);
231 m_costheta2_lab = cos(emctht[1]);
232 m_phi1_lab = emcphi[0] * 180.0 / pai;
233 m_phi2_lab = emcphi[1] * 180.0 / pai;
234 m_dlttheta_lab = dlttheta_lab;
235 m_dltphi_lab = dltphi_lab;
236
237 if(fabs(m_costheta1_lab) > m_costheta || fabs(m_costheta2_lab) > m_costheta) return sc;
238 m_pass[4]+=1;
239 //P4 in Lab
240
241 //P4 in Cms
242 double px1 = p4Gam_1.px();
243 double py1 = p4Gam_1.py();
244 double pz1 = p4Gam_1.pz();
245 double pxy1 = sqrt(px1 * px1 + py1 * py1);
246 double e1 = p4Gam_1.e();
247
248 double px2 = p4Gam_2.px();
249 double py2 = p4Gam_2.py();
250 double pz2 = p4Gam_2.pz();
251 double pxy2 = sqrt(px2 * px2 + py2 * py2);
252 double e2 = p4Gam_2.e();
253
254
255 double phi1 = 0;
256 double phi2 = 0;
257 if(atan(py1 * 1.0 / px1) > 0){
258 if(px1 > 0)phi1 = atan(py1 * 1.0 / px1);
259 else phi1 = -(pai - atan(py1 * 1.0 / px1));
260 }
261 else{
262 if(px1 > 0)phi1 = atan(py1 * 1.0 / px1);
263 else phi1 = pai + atan(py1 * 1.0 / px1);
264 }
265
266 if(atan(py2 * 1.0 / px2) > 0){
267 if(px2 > 0)phi2 = atan(py2 * 1.0 / px2);
268 else phi2 = -(pai - atan(py2 * 1.0 / px2));
269 }
270 else{
271 if(px2 > 0)phi2 = atan(py2 * 1.0 / px2);
272 else phi2 = pai + atan(py2 * 1.0 / px2);
273 }
274
275 double dltphi = (fabs(phi1 - phi2) - pai) * 180.0 / pai;
276
277 double theta1 = 0;
278 double theta2 = 0;
279
280 if(pz1 > 0) theta1 = atan(pxy1 * 1.0 / pz1);
281 else theta1 = (pai + atan(pxy1 * 1.0 / pz1));
282
283 if(pz2 > 0) theta2 = atan(pxy2 * 1.0 / pz2);
284 else theta2 = (pai + atan(pxy2 * 1.0 / pz2));
285
286 double dlttheta = (theta1 + theta2 - pai) * 180.0 / pai;
287
288 double xBoost = p4Gam_1.px() + p4Gam_2.px();
289 double yBoost = p4Gam_1.py() + p4Gam_2.py();
290 double zBoost = p4Gam_1.pz() + p4Gam_2.pz();
291 m_xBoost = xBoost;
292 m_yBoost = yBoost;
293 m_zBoost = zBoost;
294
295 m_costheta1 = cos(theta1);
296 m_costheta2 = cos(theta2);
297 m_phi1 = phi1 * 180.0 / pai;
298 m_phi2 = phi2 * 180.0 / pai;
299 m_dlttheta = dlttheta;
300 m_dltphi = dltphi;
301 m_e1 = e1;
302 m_e2 = e2;
303 m_e = e1 + e2;
304 //P4 in Cms
305
306 if(fabs(m_dltphi) < m_dphi1) Ndata1++;
307 if(fabs(m_dltphi) > m_dphi1 && fabs(m_dltphi) < m_dphi2) Ndata2++;
308
309
310 m_tuple1->write();
311// runNo = 0;
312
313 return StatusCode::SUCCESS;
314}
315
316//*************************************************************************************************************
317StatusCode TwoGamma::finalize()
318{
319 char head[100];
320 char foot[100] = ".txt";
321 sprintf(head,"OffLineLum_%04d",m_runNo);
322 strcat(head,foot);
323 ofstream fout(head,ios::out);
324 fout<<"==============================================="<<endl;
325 fout<<"DESIGNED For OffLine Luminosity BY 2Gam Events"<<endl<<endl;
326 fout<<"MANY THANKS to Prof. C.Z Yuan & Z.Y Wang"<<endl<<endl;
327 fout<<" 2009/05"<<endl;
328 // fout<<" Charmonium Group"<<endl;
329 fout<<"==============================================="<<endl<<endl<<endl;
330
331 double lum = (Ndata1 - Ndata2) * 1.0 / (m_eff * m_sec);
332 fout<<" Table List "<<endl<<endl;
333 fout<<"-----------------------------------------------"<<endl;
334 fout<<"Firstly Energy Gamma "<<m_max1<<" Gev"<<endl;
335 fout<<"Secondly Energy Gamma "<<m_max2<<" Gev"<<endl;
336 fout<<"Solid Angle Acceptance "<<m_costheta<<endl;
337 fout<<"QED XSection "<<m_sec<<" NB"<<endl;
338 fout<<"Monte Carto Efficiency "<<m_eff * 100<<"%"<<endl<<endl;
339 fout<<"runNo Luminosity "<<endl<<endl;
340 fout<<m_runNo<<" "<<lum<<" nb^(-1)"<<endl;
341 fout<<"-----------------------------------------------"<<endl;
342 fout.close();
343
344 MsgStream log(msgSvc(), name());
345 cout<< "in finalize()" <<endl;
346 cout<< "all events " <<m_pass[0]<<endl;
347 cout<< "NCharged==0 " <<m_pass[1]<<endl;
348 cout<< "Number of Gam " <<m_pass[2]<<endl;
349 cout<< "N events " <<m_pass[3]<<endl;
350 cout<< "N events 0.8 " <<m_pass[4]<<endl;
351 return StatusCode::SUCCESS;
352}
353
pai
const double pai
Definition: BbEmc.cxx:48
Phi
double Phi(RecMdcKalTrack *trk)
Definition: CalibEventSelect.cxx:99
runNo
int runNo
Definition: DQA_TO_DB.cxx:12
phi2
Double_t phi2
Definition: DataBase/tau_mode.c:7
lum
Double_t lum[10]
Definition: DataBase/tau_mode.c:36
dltphi
Double_t dltphi
Definition: DataBase/tau_mode.c:7
phi1
Double_t phi1
Definition: DataBase/tau_mode.c:7
time
Double_t time
Definition: DataBase/tau_mode.c:7
e1
Double_t e1
Definition: DataBase/tau_mode.c:7
e2
Double_t e2
Definition: DataBase/tau_mode.c:7
sin
double sin(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:213
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition: TwoGamma.cxx:29
Vp4
std::vector< HepLorentzVector > Vp4
Definition: TwoGamma.cxx:39
xmass
const double xmass[5]
Definition: TwoGamma.cxx:42
pai
const double pai
Definition: TwoGamma.cxx:43
velc
const double velc
Definition: TwoGamma.cxx:41
ecms
const double ecms
Definition: TwoGamma.cxx:40
Vint
std::vector< int > Vint
Definition: TwoGamma.cxx:38
TwoGamma::execute
StatusCode execute()
Definition: TwoGamma.cxx:135
TwoGamma::TwoGamma
TwoGamma(const std::string &name, ISvcLocator *pSvcLocator)
Definition: TwoGamma.cxx:46
TwoGamma::finalize
StatusCode finalize()
Definition: TwoGamma.cxx:317
TwoGamma::initialize
StatusCode initialize()
Definition: TwoGamma.cxx:73