BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
DQA_MDC.cxx
Go to the documentation of this file.
1#include <vector>
2
3#include "CLHEP/Vector/LorentzVector.h"
4#include "CLHEP/Vector/ThreeVector.h"
5#include "GaudiKernel/Algorithm.h"
6#include "GaudiKernel/Bootstrap.h"
7#include "GaudiKernel/INTupleSvc.h"
8#include "GaudiKernel/ITHistSvc.h"
9#include "GaudiKernel/MsgStream.h"
10#include "GaudiKernel/NTuple.h"
11#include "GaudiKernel/PropertyMgr.h"
12#include "GaudiKernel/SmartDataPtr.h"
13
14#include "EventModel/Event.h"
15#include "EventModel/EventHeader.h"
16#include "EventModel/EventModel.h"
17#include "MdcRecEvent/RecMdcHit.h"
18#include "MdcRecEvent/RecMdcKalTrack.h"
19#include "MdcRecEvent/RecMdcTrack.h"
20
21#include "DQAEvent/DQAEvent.h"
22#include "DQA_MDC.h"
23
24using CLHEP::Hep3Vector;
25/////////////////////////////////////////////////////////////////////////////
26DECLARE_COMPONENT( DQA_MDC )
27DQA_MDC::DQA_MDC( const std::string& name, ISvcLocator* pSvcLocator )
28 : Algorithm( name, pSvcLocator ) {
29
30 // Declare the properties
31}
32
33// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
34StatusCode DQA_MDC::initialize() {
35 MsgStream log( msgSvc(), name() );
36
37 log << MSG::INFO << "in initialize()" << endmsg;
38 StatusCode status;
39
40 // Call Histogram service
41 if ( service( "THistSvc", m_thsvc ).isFailure() )
42 {
43 log << MSG::ERROR << "Couldn't get THistSvc" << endmsg;
44 return StatusCode::FAILURE;
45 }
46
47 // Resolution histograms
48 m_hresAllIncBb = new TH1F( "HResAllIncBb", "", 200, -1.0, 1.0 );
49 if ( m_thsvc->regHist( "/DQAHist/MDC/hresAllIncBb", m_hresAllIncBb ).isFailure() )
50 { log << MSG::ERROR << "Couldn't register HResAllIncBb" << endmsg; }
51
52 m_hresAllExcBb = new TH1F( "HResAllExcBb", "", 200, -1.0, 1.0 );
53 if ( m_thsvc->regHist( "/DQAHist/MDC/hresAllExcBb", m_hresAllExcBb ).isFailure() )
54 { log << MSG::ERROR << "Couldn't register HResAllExcBb" << endmsg; }
55
56 m_hresAllEvaBb = new TH1F( "HResAllEvaBb", "", 200, -1.0, 1.0 );
57 if ( m_thsvc->regHist( "/DQAHist/MDC/hresAllEvaBb", m_hresAllEvaBb ).isFailure() )
58 { log << MSG::ERROR << "Couldn't register HResAllEvaBb" << endmsg; }
59
60 m_ppLabBb = new TH1F( "PpLabBb", "", 800, 0, 3 );
61 if ( m_thsvc->regHist( "/DQAHist/MDC/hppLabBb", m_ppLabBb ).isFailure() )
62 { log << MSG::ERROR << "Couldn't register PpLabBb" << endmsg; }
63 m_pmLabBb = new TH1F( "PmLabBb", "", 800, 0, 3 );
64 if ( m_thsvc->regHist( "/DQAHist/MDC/hpmLabBb", m_pmLabBb ).isFailure() )
65 { log << MSG::ERROR << "Couldn't register PmLabBb" << endmsg; }
66
67 m_ppCmsBb = new TH1F( "PpCmsBb", "", 800, 0, 3 );
68 if ( m_thsvc->regHist( "/DQAHist/MDC/hppCmsBb", m_ppCmsBb ).isFailure() )
69 { log << MSG::ERROR << "Couldn't register PpCmsBb" << endmsg; }
70 m_pmCmsBb = new TH1F( "PmCmsBb", "", 800, 0, 3 );
71 if ( m_thsvc->regHist( "/DQAHist/MDC/hpmCmsBb", m_pmCmsBb ).isFailure() )
72 { log << MSG::ERROR << "Couldn't register PmCmsBb" << endmsg; }
73
74 m_pTotLabBb = new TH1F( "PTotLabBb", "", 800, 0, 3 );
75 if ( m_thsvc->regHist( "/DQAHist/MDC/hpTotLabBb", m_pTotLabBb ).isFailure() )
76 { log << MSG::ERROR << "Couldn't register PTotLabBb" << endmsg; }
77 m_pTotCmsBb = new TH1F( "PTotCmsBb", "", 800, 0, 3 );
78 if ( m_thsvc->regHist( "/DQAHist/MDC/hpTotCmsBb", m_pTotCmsBb ).isFailure() )
79 { log << MSG::ERROR << "Couldn't register PTotCmsBb" << endmsg; }
80
81 m_hresAllIncDimu = new TH1F( "HResAllIncDimu", "", 200, -1.0, 1.0 );
82 if ( m_thsvc->regHist( "/DQAHist/MDC/hresAllIncDimu", m_hresAllIncDimu ).isFailure() )
83 { log << MSG::ERROR << "Couldn't register HResAllIncDimu" << endmsg; }
84
85 m_hresAllExcDimu = new TH1F( "HResAllExcDimu", "", 200, -1.0, 1.0 );
86 if ( m_thsvc->regHist( "/DQAHist/MDC/hresAllExcDimu", m_hresAllExcDimu ).isFailure() )
87 { log << MSG::ERROR << "Couldn't register HResAllExcDimu" << endmsg; }
88
89 m_hresAllEvaDimu = new TH1F( "HResAllEvaDimu", "", 200, -1.0, 1.0 );
90 if ( m_thsvc->regHist( "/DQAHist/MDC/hresAllEvaDimu", m_hresAllEvaDimu ).isFailure() )
91 { log << MSG::ERROR << "Couldn't register HResAllEvaDimu" << endmsg; }
92
93 m_ppLabDimu = new TH1F( "PpLabDimu", "", 800, 0, 3 );
94 if ( m_thsvc->regHist( "/DQAHist/MDC/hppLabDimu", m_ppLabDimu ).isFailure() )
95 { log << MSG::ERROR << "Couldn't register PpLabDimu" << endmsg; }
96 m_pmLabDimu = new TH1F( "PmLabDimu", "", 800, 0, 3 );
97 if ( m_thsvc->regHist( "/DQAHist/MDC/hpmLabDimu", m_pmLabDimu ).isFailure() )
98 { log << MSG::ERROR << "Couldn't register PmLabDimu" << endmsg; }
99
100 m_ppCmsDimu = new TH1F( "PpCmsDimu", "", 800, 0, 3 );
101 if ( m_thsvc->regHist( "/DQAHist/MDC/hppCmsDimu", m_ppCmsDimu ).isFailure() )
102 { log << MSG::ERROR << "Couldn't register PpCmsDimu" << endmsg; }
103 m_pmCmsDimu = new TH1F( "PmCmsDimu", "", 800, 0, 3 );
104 if ( m_thsvc->regHist( "/DQAHist/MDC/hpmCmsDimu", m_pmCmsDimu ).isFailure() )
105 { log << MSG::ERROR << "Couldn't register PmCmsDimu" << endmsg; }
106
107 m_pTotLabDimu = new TH1F( "PTotLabDimu", "", 800, 0, 3 );
108 if ( m_thsvc->regHist( "/DQAHist/MDC/hpTotLabDimu", m_pTotLabDimu ).isFailure() )
109 { log << MSG::ERROR << "Couldn't register PTotLabDimu" << endmsg; }
110 m_pTotCmsDimu = new TH1F( "PTotCmsDimu", "", 800, 0, 3 );
111 if ( m_thsvc->regHist( "/DQAHist/MDC/hpTotCmsDimu", m_pTotCmsDimu ).isFailure() )
112 { log << MSG::ERROR << "Couldn't register PTotCmsDimu" << endmsg; }
113
114 log << MSG::INFO << "Initialize done!" << endmsg;
115 return StatusCode::SUCCESS;
116}
117
118// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
119StatusCode DQA_MDC::execute() {
120
121 MsgStream log( msgSvc(), name() );
122 log << MSG::INFO << "in execute()" << endmsg;
123
124 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
125 m_run = eventHeader->runNumber();
126 m_event = eventHeader->eventNumber();
127 log << MSG::DEBUG << "Run " << m_run << "\tEvent " << m_event << endmsg;
128
129 SmartDataPtr<DQAEvent::DQAEvent> dqaevt( eventSvc(), "/Event/DQATag" );
130 if ( dqaevt ) { log << MSG::INFO << "success get DQAEvent" << endmsg; }
131 else
132 {
133 log << MSG::ERROR << "Error accessing DQAEvent" << endmsg;
134 return StatusCode::FAILURE;
135 }
136 log << MSG::DEBUG << "DQA event tag = " << dqaevt->EventTag() << endmsg;
137
138 double p;
139 double p_cms;
140 bool isDimu = false;
141 if ( dqaevt->Dimu() ) isDimu = true;
142 log << MSG::INFO << "DQADimuTag:\t" << dqaevt->Dimu() << endmsg;
143 // cout << "DQADimuTag:\t" << dqaevt->Dimu() << endl;
144
145 bool isBb = false;
146 if ( dqaevt->Bhabha() ) isBb = true;
147 log << MSG::INFO << "DQABbTag:\t" << dqaevt->Bhabha() << endmsg;
148 // cout << "DQABbTag:\t" << dqaevt->Bhabha() << endl;
149
150 SmartDataPtr<RecMdcKalTrackCol> kaltrkCol( eventSvc(), "/Event/Recon/RecMdcKalTrackCol" );
151 if ( !kaltrkCol )
152 {
153 log << MSG::FATAL << "Could not find RecMdcKalTrackCol" << endmsg;
154 return StatusCode::FAILURE;
155 }
156 RecMdcKalTrackCol::iterator iter_trk = kaltrkCol->begin();
157 for ( ; iter_trk != kaltrkCol->end(); iter_trk++ )
158 {
159 if ( isBb ) RecMdcKalTrack::setPidType( RecMdcKalTrack::electron );
160 else if ( isDimu ) RecMdcKalTrack::setPidType( RecMdcKalTrack::muon );
161 // else return StatusCode::SUCCESS;
162
163 double dr = ( *iter_trk )->dr();
164 double phi0 = ( *iter_trk )->fi0();
165 double kappa = ( *iter_trk )->kappa();
166 double dz = ( *iter_trk )->dz();
167 double tanl = ( *iter_trk )->tanl();
168 double pt = 1.0 / kappa;
169 HepLorentzVector p4;
170 p4.setPx( -sin( phi0 ) / fabs( kappa ) );
171 p4.setPy( cos( phi0 ) / fabs( kappa ) );
172 p4.setPz( tanl / fabs( kappa ) );
173
174 double p3 = p4.mag();
175 double mass = 0.0;
176 if ( isBb ) mass = 0.000511;
177 else if ( isDimu ) mass = 0.105658;
178
179 p4.setE( sqrt( p3 * p3 + mass * mass ) );
180 p = p4.rho();
181
182 double ecm = 3.68632;
183 HepLorentzVector psip( 0.011 * ecm, 0, 0.0075, ecm );
184 Hep3Vector boostv = psip.boostVector();
185 p4.boost( -boostv );
186 p_cms = p4.rho();
187 TH1* hmom( 0 );
188 if ( isBb )
189 {
190 if ( m_thsvc->getHist( "/DQAHist/MDC/hpTotLabBb", hmom ).isSuccess() )
191 { hmom->Fill( p ); }
192 else { log << MSG::ERROR << "Couldn't retrieve hpTotLabBb" << endmsg; }
193 if ( m_thsvc->getHist( "/DQAHist/MDC/hpTotCmsBb", hmom ).isSuccess() )
194 { hmom->Fill( p_cms ); }
195 else { log << MSG::ERROR << "Couldn't retrieve hpTotCmsBb" << endmsg; }
196 if ( pt > 0 )
197 {
198 if ( m_thsvc->getHist( "/DQAHist/MDC/hppLabBb", hmom ).isSuccess() )
199 { hmom->Fill( p ); }
200 else { log << MSG::ERROR << "Couldn't retrieve hppLabBb" << endmsg; }
201 if ( m_thsvc->getHist( "/DQAHist/MDC/hppCmsBb", hmom ).isSuccess() )
202 { hmom->Fill( p ); }
203 else { log << MSG::ERROR << "Couldn't retrieve hppCmsBb" << endmsg; }
204 }
205 else if ( pt < 0 )
206 {
207 if ( m_thsvc->getHist( "/DQAHist/MDC/hpmLabBb", hmom ).isSuccess() )
208 { hmom->Fill( p ); }
209 else { log << MSG::ERROR << "Couldn't retrieve hpmLabBb" << endmsg; }
210 if ( m_thsvc->getHist( "/DQAHist/MDC/hpmCmsBb", hmom ).isSuccess() )
211 { hmom->Fill( p ); }
212 else { log << MSG::ERROR << "Couldn't retrieve hpmCmsBb" << endmsg; }
213 }
214 }
215 else if ( isDimu )
216 {
217 if ( m_thsvc->getHist( "/DQAHist/MDC/hpTotLabDimu", hmom ).isSuccess() )
218 { hmom->Fill( p ); }
219 else { log << MSG::ERROR << "Couldn't retrieve hpTotLabDimu" << endmsg; }
220 if ( m_thsvc->getHist( "/DQAHist/MDC/hpTotCmsDimu", hmom ).isSuccess() )
221 { hmom->Fill( p_cms ); }
222 else { log << MSG::ERROR << "Couldn't retrieve hpTotCmsDimu" << endmsg; }
223 if ( pt > 0 )
224 {
225 if ( m_thsvc->getHist( "/DQAHist/MDC/hppLabDimu", hmom ).isSuccess() )
226 { hmom->Fill( p ); }
227 else { log << MSG::ERROR << "Couldn't retrieve hppLabDimu" << endmsg; }
228 if ( m_thsvc->getHist( "/DQAHist/MDC/hppCmsDimu", hmom ).isSuccess() )
229 { hmom->Fill( p ); }
230 else { log << MSG::ERROR << "Couldn't retrieve hppCmsDimu" << endmsg; }
231 }
232 else if ( pt < 0 )
233 {
234 if ( m_thsvc->getHist( "/DQAHist/MDC/hpmLabDimu", hmom ).isSuccess() )
235 { hmom->Fill( p ); }
236 else { log << MSG::ERROR << "Couldn't retrieve hpmLabDimu" << endmsg; }
237 if ( m_thsvc->getHist( "/DQAHist/MDC/hpmCmsDimu", hmom ).isSuccess() )
238 { hmom->Fill( p ); }
239 else { log << MSG::ERROR << "Couldn't retrieve hpmCmsDimu" << endmsg; }
240 }
241 }
242 }
243
244 SmartDataPtr<RecMdcTrackCol> newtrkCol( eventSvc(), "/Event/Recon/RecMdcTrackCol" );
245 if ( !newtrkCol )
246 {
247 log << MSG::ERROR << "Could not find RecMdcTrackCol" << endmsg;
248 return ( StatusCode::FAILURE );
249 }
250
251 RecMdcTrackCol::iterator it_trk = newtrkCol->begin();
252 for ( ; it_trk != newtrkCol->end(); it_trk++ )
253 {
254 HitRefVec gothits = ( *it_trk )->getVecHits();
255 HitRefVec::iterator it_hit = gothits.begin();
256 for ( ; it_hit != gothits.end(); it_hit++ )
257 {
258
259 // cout << __LINE__ << endl;
260 // cout << __LINE__ << endl;
261 double lr = ( *it_hit )->getFlagLR();
262 if ( -1 == lr ) lr = 0; // definition not same as MdcRecHit
263 double resilrInc;
264 double resilrExc;
265 // for boss before 650
266 double docaInc = ( *it_hit )->getDoca();
267 double docaExc = docaInc;
268 // boss650
269 // double docaInc = (*it_hit) -> getDocaIncl();
270 // double docaExc = (*it_hit) -> getDocaExcl();
271 double dmeas;
272 if ( 0 == lr ) { dmeas = ( *it_hit )->getDriftDistLeft(); }
273 else { dmeas = ( *it_hit )->getDriftDistRight(); }
274 // the following is for cm to mm
275 double df = 10.0;
276 docaInc *= df;
277 docaExc *= df;
278 dmeas *= df;
279 double resiInc = fabs( dmeas ) - fabs( docaInc );
280 if ( 0 == lr ) resilrInc = -1.0 * resiInc;
281 else resilrInc = resiInc;
282
283 double resiExc = fabs( dmeas ) - fabs( docaExc );
284 if ( 0 == lr ) resilrExc = -1.0 * resiExc;
285 else resilrExc = resiExc;
286
287 double resilrEva = 0.5 * ( resilrInc + resilrExc );
288 // cout << "resilrInc = " << resilrInc << ", resilrExc = " << resilrExc << ", resilrEva =
289 // " << resilrEva << endl;
290 TH1* htmp( 0 );
291 if ( isBb )
292 {
293 if ( m_thsvc->getHist( "/DQAHist/MDC/hresAllIncBb", htmp ).isSuccess() )
294 { htmp->Fill( resilrInc ); }
295 else { log << MSG::ERROR << "Couldn't retrieve hresAllIncBb" << endmsg; }
296 if ( m_thsvc->getHist( "/DQAHist/MDC/hresAllExcBb", htmp ).isSuccess() )
297 { htmp->Fill( resilrExc ); }
298 else { log << MSG::ERROR << "Couldn't retrieve hresAllExcBb" << endmsg; }
299 if ( m_thsvc->getHist( "/DQAHist/MDC/hresAllEvaBb", htmp ).isSuccess() )
300 { htmp->Fill( resilrEva ); }
301 else { log << MSG::ERROR << "Couldn't retrieve hresAllEvaBb" << endmsg; }
302 }
303 else if ( isDimu )
304 {
305 if ( m_thsvc->getHist( "/DQAHist/MDC/hresAllIncDimu", htmp ).isSuccess() )
306 { htmp->Fill( resilrInc ); }
307 else { log << MSG::ERROR << "Couldn't retrieve hresAllIncDimu" << endmsg; }
308 if ( m_thsvc->getHist( "/DQAHist/MDC/hresAllExcDimu", htmp ).isSuccess() )
309 { htmp->Fill( resilrExc ); }
310 else { log << MSG::ERROR << "Couldn't retrieve hresAllExcDimu" << endmsg; }
311 if ( m_thsvc->getHist( "/DQAHist/MDC/hresAllEvaDimu", htmp ).isSuccess() )
312 { htmp->Fill( resilrEva ); }
313 else { log << MSG::ERROR << "Couldn't retrieve hresAllEvaDimu" << endmsg; }
314 }
315 }
316 }
317
318 return StatusCode::SUCCESS;
319}
320
321// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
322StatusCode DQA_MDC::finalize() {
323
324 MsgStream log( msgSvc(), name() );
325 log << MSG::INFO << "in finalize()" << endmsg;
326 log << MSG::INFO << "Finalize successfully! " << endmsg;
327
328 return StatusCode::SUCCESS;
329}
DECLARE_COMPONENT
DECLARE_COMPONENT(BesBdkRc)
mass
double mass
Definition CosmicGenerator.cxx:128
p_cms
const HepLorentzVector p_cms(0.034067, 0.0, 0.0, 3.097)
DQA_MDC.h
sin
double sin(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:185
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187
HitRefVec
SmartRefVector< RecMdcHit > HitRefVec
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcRecEvent/RecMdcTrack.h:22
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:12
DQA_MDC
Definition DQA_MDC.h:11
DQA_MDC::execute
StatusCode execute()
Definition DQA_MDC.cxx:119
DQA_MDC::DQA_MDC
DQA_MDC(const std::string &name, ISvcLocator *pSvcLocator)
Definition DQA_MDC.cxx:27
DQA_MDC::initialize
StatusCode initialize()
Definition DQA_MDC.cxx:34
DQA_MDC::finalize
StatusCode finalize()
Definition DQA_MDC.cxx:322