BOSS 7.1.0
BESIII Offline Software System
Loading...
Searching...
No Matches
DimuPreSelect/DimuPreSelect-00-01-00/src/DimuPreSelect.cxx
Go to the documentation of this file.
1
2#include "GaudiKernel/MsgStream.h"
3#include "GaudiKernel/AlgFactory.h"
4#include "GaudiKernel/ISvcLocator.h"
5#include "GaudiKernel/SmartDataPtr.h"
6#include "GaudiKernel/IDataProviderSvc.h"
7#include "GaudiKernel/PropertyMgr.h"
8#include "EventModel/EventModel.h"
9#include "EventModel/Event.h"
10#include "EventModel/EventHeader.h"
11#include "EvtRecEvent/EvtRecEvent.h"
12#include "EvtRecEvent/EvtRecTrack.h"
13
14#include "MdcRecEvent/RecMdcTrack.h"
15#include "EmcRecEventModel/RecEmcEventModel.h"
16#include "TofRecEvent/RecTofTrack.h"
17#include "DstEvent/TofHitStatus.h"
18#include "MucRawEvent/MucDigi.h"
19#include "MucRecEvent/RecMucTrack.h"
20
21#include "TMath.h"
22#include "GaudiKernel/INTupleSvc.h"
23#include "GaudiKernel/NTuple.h"
24#include "GaudiKernel/Bootstrap.h"
25#include "GaudiKernel/IHistogramSvc.h"
26#include "CLHEP/Vector/ThreeVector.h"
27using CLHEP::Hep3Vector;
28#include "CLHEP/Vector/LorentzVector.h"
29using CLHEP::HepLorentzVector;
30#include "CLHEP/Vector/TwoVector.h"
31using CLHEP::Hep2Vector;
32#include "CLHEP/Geometry/Point3D.h"
33#ifndef ENABLE_BACKWARDS_COMPATIBILITY
34typedef HepGeom::Point3D<double> HepPoint3D;
35#endif
36#include "VertexFit/KinematicFit.h"
37#include "VertexFit/VertexFit.h"
38#include "VertexFit/SecondVertexFit.h"
39#include "DimuPreSelect/DimuPreSelect.h"
40#include <vector>
41typedef std::vector<int> Vint;
42typedef std::vector<HepLorentzVector> Vp4;
43const double ECMS = 3.097;
44const double VELC = 299.792458;
45const double PI = 3.14159265;
46const double XMASS[5] = { 0.000511, 0.105658, 0.139570, 0.493677, 0.938272 };
47DECLARE_COMPONENT(DimuPreSelect)
48////////////////////////////////////////////////////////////////////
49DimuPreSelect::DimuPreSelect(const std::string& name, ISvcLocator* pSvcLocator) :
50 Algorithm(name, pSvcLocator)
51{
52 declareProperty("CmsEnergy", m_ecms = 3.686);
53 declareProperty("Vr0Cut", m_vr0cut = 1.0);
54 declareProperty("Vz0Cut", m_vz0cut = 5.0);
55 declareProperty("PUpCut", m_pcut_up = 2.0);
56 declareProperty("PDownCut", m_pcut_down = 0.5);
57 declareProperty("PSymCut", m_psymcut = 0.5);
58 declareProperty("TCut", m_tcut = 4);
59 declareProperty("EUpCut", m_ecut_up = 1.0);
60 declareProperty("EDownCut", m_ecut_down = 0.1);
61 declareProperty("DThetaCut", m_dthetacut = 0.05);
62 declareProperty("DPhiCut", m_dphicut = 0.4);
63 declareProperty("PartSelect", m_partselect= 0);
64 declareProperty("MuDigiCut", m_mudigicut = 6);
65 declareProperty("MuTrkCut", m_mutrkcut = 1);
66 declareProperty("MuDepthCut", m_depthcut= 30);
67 declareProperty("UseMDC", m_useFlag[0]= 1);
68 declareProperty("UseTOF", m_useFlag[1]= 1);
69 declareProperty("UseEMC", m_useFlag[2]= 1);
70 declareProperty("UseMUC", m_useFlag[3]= 1);
71 declareProperty ("Output", m_output = false);
72
73 declareProperty ("SelectDimu",m_selectFlag= true);
74 //Declare the properties
75}
76
77// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
78StatusCode DimuPreSelect::initialize()
79{
80 MsgStream log(msgSvc(), name());
81 log << MSG::INFO << "in initialize()" << endmsg;
82
83 m_totevent = m_currun = m_curevent = 0;
84 for(int i = 0; i < 20; i++) m_cutpass[i] = 0;
85 m_subpass[0] = m_subpass[1] = m_subpass[2] = m_subpass[3] = 0;
86 m_totalpass = 0;
87
88 if(m_output)
89 {
90 StatusCode status;
91 NTuplePtr nt(ntupleSvc(), "FILE1/dimu");
92 if ( nt ) m_passtuple = nt;
93 else {
94 m_passtuple = ntupleSvc()->book ("FILE1/dimu", CLID_ColumnWiseTuple, "DimuPreSelect N-Tuple example");
95 if ( m_passtuple )
96 {
97 status =m_passtuple->addItem ("run", m_run);
98 status =m_passtuple->addItem ("event", m_event);
99 status =m_passtuple->addItem ("part", m_part);
100 status =m_passtuple->addItem ("c1", m_c1);
101 status =m_passtuple->addItem ("c2", m_c2);
102 status =m_passtuple->addItem ("r1", m_r1);
103 status =m_passtuple->addItem ("r2", m_r2);
104 status =m_passtuple->addItem ("z1", m_z1);
105 status =m_passtuple->addItem ("z2", m_z2);
106 status =m_passtuple->addItem ("p1", m_p1);
107 status =m_passtuple->addItem ("p2", m_p2);
108 status =m_passtuple->addItem ("t1", m_t1);
109 status =m_passtuple->addItem ("t2", m_t2);
110 status =m_passtuple->addItem ("e1", m_e1);
111 status =m_passtuple->addItem ("e2", m_e2);
112 status =m_passtuple->addItem ("theta1",m_theta1);
113 status =m_passtuple->addItem ("theta2",m_theta2);
114 status =m_passtuple->addItem ("phi1", m_phi1);
115 status =m_passtuple->addItem ("phi2", m_phi2);
116 status =m_passtuple->addItem ("mudigi",m_mudigi);
117 status =m_passtuple->addItem ("mutrk", m_mutrk);
118 status =m_passtuple->addItem ("depth", m_depth);
119
120 status =m_passtuple->addItem ("zeroC", m_zeroCFlag);
121 status =m_passtuple->addItem ("vtRZ", m_vtRZFlag);
122 status =m_passtuple->addItem ("pLim", m_pLimFlag);
123 status =m_passtuple->addItem ("pSym", m_pSymFlag);
124 status =m_passtuple->addItem ("tLim", m_tLimFlag);
125 status =m_passtuple->addItem ("eLim", m_eLimFlag);
126 status =m_passtuple->addItem ("eBB", m_eBBFlag);
127 status =m_passtuple->addItem ("partSlt",m_partFlag);
128 status =m_passtuple->addItem ("muDigiN",m_mudigiFlag);
129 status =m_passtuple->addItem ("muTrkN", m_mutrkFlag);
130
131 status =m_passtuple->addItem ("mdc", m_mdcFlag);
132 status =m_passtuple->addItem ("tof", m_tofFlag);
133 status =m_passtuple->addItem ("emc", m_emcFlag);
134 status =m_passtuple->addItem ("muc", m_mucFlag);
135
136 }
137 else
138 {
139 log << MSG::ERROR << "Cannot book N-tuple:"<<long(m_passtuple)<<endmsg;
140 return StatusCode::FAILURE;
141 }
142 }
143 }
144
145 log << MSG::INFO << "Initialize done!" << endmsg;
146
147 return StatusCode::SUCCESS;
148}
149
150//*************************************************
151StatusCode DimuPreSelect::execute()
152{
153 if(m_selectFlag==false) return( StatusCode::SUCCESS );
154
155 MsgStream log(msgSvc(),name());
156 log<<MSG::INFO<<"in execute()"<<endreq;
157
158 m_totevent ++;
159 if(m_totevent%1000==0) std::cout << m_totevent << "\tdone!" <<std::endl;
160
161 setFilterPassed(false);
162 m_mdcPass = m_tofPass = m_emcPass = m_mucPass = false;
163 for(int i=0; i<4; i++) m_passFlag[i] = false;
164
165 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
166 m_currun = eventHeader->runNumber();
167 m_curevent = eventHeader->eventNumber();
168 if( m_output ) { m_run = m_currun; m_event = m_curevent; }
169
170 // Select by MDC Info
171 SmartDataPtr<RecMdcTrackCol> mdcTrackCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
172 if(!mdcTrackCol) {
173 log << MSG::FATAL << "Could not find RecMdcTrackCol!" << endreq;
174 return( StatusCode::FAILURE);
175 }
176 log << MSG::INFO << "MDC tracks:\t" << mdcTrackCol->size() << endreq;
177
178 if( mdcTrackCol->size() != 2 ) return ( StatusCode::SUCCESS );
179
180 m_cutpass[0] += 1;
181
182 double c1, c2, r1, r2, z1, z2, p1, p2;
183 c1 = c2 = r1 = r2 = z1 = z2 = p1 = p2 = 0.;
184
185 c1 = (*mdcTrackCol)[0]->charge(); c2 = (*mdcTrackCol)[1]->charge();
186 r1 = (*mdcTrackCol)[0]->r(); r2 = (*mdcTrackCol)[1]->r();
187 z1 = (*mdcTrackCol)[0]->z(); z2 = (*mdcTrackCol)[1]->z();
188 p1 = (*mdcTrackCol)[0]->p(); p2 = (*mdcTrackCol)[1]->p();
189
190 if( m_output ) { m_c1 = c1; m_c2 = c2; m_r1 = r1; m_r2 = r2; m_z1 = z1; m_z2 = z2; m_p1 = p1; m_p2 = p2; }
191
192 bool mdcflag1 = c1 + c2 == 0;
193 bool mdcflag2 = fabs(r1)<=m_vr0cut && fabs(z1)<m_vz0cut;
194 bool mdcflag3 = fabs(r2)<=m_vr0cut && fabs(z2)<m_vz0cut;
195 bool mdcflag4 = p1<m_pcut_up && p2<m_pcut_up;
196 //bool mdcflag4 = p1<2 && p1>1;
197 //bool mdcflag5 = p2<2 && p2>1;
198 bool mdcflag5 = fabs( p1-p2 )/( p1+p2 ) < m_psymcut;
199
200 log << MSG::INFO << "r1:\t" << r1 << "\tz1:" << z1 << endreq;
201 log << MSG::INFO << "r2:\t" << r2 << "\tz2:" << z2 << endreq;
202 log << MSG::INFO << "p1:\t" << p1 << "\tp2:" << p2 << endreq;
203
204 if( mdcflag1 ) { m_cutpass[1] += 1; if(m_output) m_zeroCFlag = 1; }
205 if( mdcflag2 && mdcflag3 ) { m_cutpass[2] += 1; if(m_output) m_vtRZFlag = 1; }
206 if( mdcflag4 ) { m_cutpass[3] += 1; if(m_output) m_pLimFlag = 1; }
207 if( mdcflag5 ) { m_cutpass[4] += 1; if(m_output) m_pSymFlag = 1; }
208 if( mdcflag1 && mdcflag2 && mdcflag3 && mdcflag4 && mdcflag5 ) { m_mdcPass = true; m_subpass[0] += 1; }
209 if( !m_useFlag[0] ) m_passFlag[0] = true;
210 else m_passFlag[0] = m_mdcPass;
211 log << MSG::INFO << "MDC selection done!" << endreq;
212
213 // Select by TOF Info
214 SmartDataPtr<RecTofTrackCol> tofTrackCol(eventSvc(),"/Event/Recon/RecTofTrackCol");
215 if(!tofTrackCol) {
216 log << MSG::FATAL << "Could not find RecTofTrackCol!" << endreq;
217 return( StatusCode::FAILURE);
218 }
219 log << MSG::INFO << "TOF tracks:\t" << tofTrackCol->size() << endreq;
220
221 double t1, t2;
222 t1 = 0., t2 = 1000;
223 // if(tofTrackCol->size() < 8 && tofTrackCol->size() > 20)
224 if(tofTrackCol->size() > 7 && tofTrackCol->size() < 21)
225 {
226 int goodtof = 0;
227 for(int itof = 0; itof < tofTrackCol->size(); itof++)
228 {
229 TofHitStatus *status = new TofHitStatus;
230 status->setStatus((*tofTrackCol)[itof]->status());
231
232 if( !(status->is_cluster()) ) { delete status; continue; }
233 if(goodtof==0) t1 = (*tofTrackCol)[itof]->tof();
234 if(goodtof==1) t2 = (*tofTrackCol)[itof]->tof();
235
236 goodtof++;
237 delete status;
238 }
239 }
240
241 if( m_output ) { m_t1 = t1; m_t2 = t2; }
242
243 bool tofflag1 = fabs( t1-t2 ) < m_tcut; // ns
244 log << MSG::INFO << "t1:\t" << t1 << "\tt2:" << t2 << "dt:\t" << fabs(t1-t2) << endreq;
245 if( tofflag1 ) { m_cutpass[5] += 1; if(m_output) m_tLimFlag = 1; m_tofPass = true; m_subpass[1] += 1; }
246 if( !m_useFlag[1] ) m_passFlag[1] = true;
247 else m_passFlag[1] = m_tofPass;
248 log << MSG::INFO << "TOF selection done!" << endreq;
249
250 // Select by EMC Info
251 SmartDataPtr<RecEmcShowerCol> emcShowerCol(eventSvc(),"/Event/Recon/RecEmcShowerCol");
252 if (!emcShowerCol) {
253 log << MSG::FATAL << "Could not find RecEmcShowerCol!" << endreq;
254 return( StatusCode::FAILURE);
255 }
256 log << MSG::INFO << "EMC showers:\t" << emcShowerCol->size() << endreq;
257
258 if( emcShowerCol->size() < 2 ) return ( StatusCode::SUCCESS );
259
260 double e1, e2, theta1, theta2, phi1, phi2;
261 int part;
262
263 e1 = (*emcShowerCol)[0]->energy(); e2 = (*emcShowerCol)[1]->energy();
264 theta1 = (*emcShowerCol)[0]->theta(); theta2 = (*emcShowerCol)[1]->theta();
265 phi1 = (*emcShowerCol)[0]->phi(); phi2 = (*emcShowerCol)[1]->phi();
266 part = (*emcShowerCol)[0]->module();
267
268 if( m_output ) { m_e1 = e1; m_e2 = e2; m_theta1 = theta1; m_theta2 = theta2; m_phi1 = phi1; m_phi2 = phi2;
269 m_part = part;
270 }
271
272 bool emcFlag1 = e1 < m_ecut_up && e1 > m_ecut_down;
273 bool emcFlag2 = e2 < m_ecut_up && e2 > m_ecut_down;
274 bool emcFlag3 = fabs(theta1 + theta2 - PI) < m_dthetacut;
275 bool emcFlag4 = fabs(phi1 - phi2) - PI < m_dphicut;
276 bool emcFlag5 = !m_partselect || (m_partselect==1&&part==1) || (m_partselect==2&&part!=1);
277
278 log << MSG::INFO << "e1:\t" << e1 << "\te2:\t" << e2 << endreq;
279 log << MSG::INFO << "theta1:\t" << theta1 << "\ttheta2:\t" << theta2 << endreq;
280 log << MSG::INFO << "phi1:\t" << phi1 << "\tphi2:\t" << phi2 << endreq;
281 log << MSG::INFO << "part:\t" << part << "\tpartFlag:\t" << emcFlag5 << endreq;
282
283 if( emcFlag1 && emcFlag2 ) { m_cutpass[6] += 1; if(m_output) m_eLimFlag = 1; }
284 if( emcFlag3 && emcFlag4 ) { m_cutpass[7] += 1; if(m_output) m_eBBFlag = 1; }
285 if( emcFlag5 ) { m_cutpass[8] += 1; if(m_output) m_partFlag = 1; }
286 if( emcFlag1 && emcFlag2 && emcFlag3 && emcFlag4 && emcFlag5 ) { m_emcPass = true; m_subpass[2] += 1; }
287 if( !m_useFlag[2] ) m_passFlag[2] = true;
288 else m_passFlag[2] = m_emcPass;
289 log << MSG::INFO << "EMC selection done!" << endreq;
290
291 // Select by MUC Info
292 SmartDataPtr<MucDigiCol> mucDigiCol(eventSvc(),"/Event/Digi/MucDigiCol");
293 if(!mucDigiCol) {
294 log << MSG::FATAL << "Could not find MucDigiCol!" << endreq;
295 return( StatusCode::FAILURE);
296 }
297 SmartDataPtr<RecMucTrackCol> mucTrackCol(eventSvc(),"/Event/Recon/RecMucTrackCol");
298 if (!mucTrackCol) {
299 log << MSG::FATAL << "Could not find RecMucTrackCol" << endreq;
300 return( StatusCode::FAILURE);
301 }
302
303 int mudigiNum, mutrkNum;
304 double depth = 0;
305 mudigiNum = mutrkNum = 0;
306 mudigiNum = mucDigiCol->size(); mutrkNum = mucTrackCol->size();
307 RecMucTrackCol::iterator mutrkIter = mucTrackCol->begin();
308 for(; mutrkIter != mucTrackCol->end(); mutrkIter++) {
309 if( 0 == (*mutrkIter)->trackId() || 1 == (*mutrkIter)->trackId() ) depth += (*mutrkIter)->depth();
310 }
311
312 if(m_output) { m_mudigi = mudigiNum; m_mutrk = mutrkNum; m_depth = depth; }
313
314 bool mucflag1 = mudigiNum >= m_mudigicut;
315 bool mucflag2 = mutrkNum >= m_mutrkcut;
316 bool mucflag3 = depth > m_depthcut;
317
318 log << MSG::INFO << "MUC digi:\t" << mudigiNum << "\tMUC track:\t" << mutrkNum << endreq;
319
320 if( mucflag1 ) { m_cutpass[9] += 1; if(m_output) m_mudigiFlag = 1; }
321 if( mucflag2 ) { m_cutpass[10] += 1; if(m_output) m_mutrkFlag = 1; }
322 if( mucflag3 ) { m_cutpass[11] += 1; if(m_output) m_depthFlag = 1; }
323 if( mucflag1 && mucflag2 && mucflag3 ) { m_mucPass = true; m_subpass[3] += 1; }
324 if( !m_useFlag[3] ) m_passFlag[3] = true;
325 else m_passFlag[3] = m_mucPass;
326 log << MSG::INFO << "MUC selection done!" << endreq;
327
328
329 // All selections
330 if( m_passFlag[0] && m_passFlag[1] && m_passFlag[2] && m_passFlag[3] )
331 {
332 m_totalpass += 1;
333 setFilterPassed(true);
334 }
335 log << MSG::INFO << "Set filter passed!" << endreq;
336
337 if(m_output) { m_mdcFlag = m_mdcPass; m_tofFlag = m_tofPass; m_emcFlag = m_emcPass; m_mucFlag = m_mucPass; }
338
339 if( m_output ) m_passtuple->write();
340
341 return StatusCode::SUCCESS;
342}
343
344StatusCode DimuPreSelect::finalize()
345{
346 if(m_selectFlag==false) return StatusCode::SUCCESS;
347
348 MsgStream log(msgSvc(), name());
349 log << MSG::INFO << "in finalize()" << endmsg;
350 const string str[3] = {"all","barrel","endcap"};
351
352 cout << "pass 0 - 2 MDC tracks : " << m_cutpass[0] << endl;
353 cout << "pass 1 - 0 charge : " << m_cutpass[1] << endl;
354 cout << "pass 2 - |r|<"<<m_vr0cut<<" && |z|<"<<m_vz0cut<<" : " << m_cutpass[2] << endl;
355 cout << "pass 3 - p < "<<m_pcut_up<<" GeV/c : " << m_cutpass[3] << endl;
356 cout << "pass 4 - p_sym < "<<m_psymcut<<" : " << m_cutpass[4] << endl;
357 cout << "pass 5 - |t1-t1| < "<<m_tcut<<" ns : " << m_cutpass[5] << endl;
358 cout << "pass 6 - "<<m_ecut_down<<" < e < "<<m_ecut_up<<" : " << m_cutpass[6] << endl;
359 cout << "pass 7 - |dth|<"<<m_dthetacut<<" && |dphi|<"<<m_dphicut<<": " << m_cutpass[7] << endl;
360 cout << "pass 8 - "<< str[(int)m_partselect] <<" part is selected : " << m_cutpass[8] << endl;
361 cout << "pass 9 - mudigi >= "<<m_mudigicut<<" : " << m_cutpass[9] << endl;
362 cout << "pass 10- mutrk >= "<<m_mutrkcut<<" : " << m_cutpass[10]<< endl;
363 cout << "pass 11- mudepth>= "<<m_depthcut<<" : " << m_cutpass[11]<< endl;
364
365 cout << "pass MDC : " << m_subpass[0] << "\tUsed: " << (m_useFlag[0]?"Yes":"No") << endl;
366 cout << "pass TOF : " << m_subpass[1] << "\tUsed: " << (m_useFlag[1]?"Yes":"No") << endl;
367 cout << "pass EMC : " << m_subpass[2] << "\tUsed: " << (m_useFlag[2]?"Yes":"No") << endl;
368 cout << "pass MUC : " << m_subpass[3] << "\tUsed: " << (m_useFlag[3]?"Yes":"No") << endl;
369
370
371 cout<<"Total event: " << m_totevent <<endl;
372 cout<<"Dimu event: " << m_totalpass <<endl;
373
374 return StatusCode::SUCCESS;
375}
phi2
Double_t phi2
Definition: DataBase/tau_mode.c:7
phi1
Double_t phi1
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
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition: DimuPreSelect/DimuPreSelect-00-01-00/src/DimuPreSelect.cxx:34
Vp4
std::vector< HepLorentzVector > Vp4
Definition: DimuPreSelect/DimuPreSelect-00-01-00/src/DimuPreSelect.cxx:42
PI
const double PI
Definition: PipiJpsi.cxx:55
ntupleSvc
INTupleSvc * ntupleSvc()
Definition: ServiceAccessor.h:55
msgSvc
IMessageSvc * msgSvc()
Definition: ServiceAccessor.h:13
theta2
@ theta2
Definition: TrkKalDeriv.h:24
theta1
@ theta1
Definition: TrkKalDeriv.h:24
TofHitStatus::is_cluster
bool is_cluster() const
Definition: TofHitStatus.h:25
TofHitStatus::setStatus
void setStatus(unsigned int status)
Definition: TofHitStatus.cxx:134
p2
double * p2
Definition: qcdloop1.h:76