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