BOSS 7.0.1
BESIII Offline Software System
Loading...
Searching...
No Matches
DQAJpsi2PPbarAlg.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
8#include "EventModel/EventModel.h"
9#include "EventModel/Event.h"
10#include "EventModel/EventHeader.h"
11
12#include "EvtRecEvent/EvtRecEvent.h"
13#include "EvtRecEvent/EvtRecTrack.h"
14#include "DstEvent/TofHitStatus.h"
15
16// #include "TMath.h"
17#include "TH1F.h"
18#include "VertexFit/KinematicFit.h"
19#include "VertexFit/VertexFit.h"
20#include "VertexFit/SecondVertexFit.h"
21#include "VertexFit/IVertexDbSvc.h"
22
23#include "ParticleID/ParticleID.h"
24
25
26#include "TMath.h"
27#include "GaudiKernel/INTupleSvc.h"
28#include "GaudiKernel/NTuple.h"
29#include "GaudiKernel/Bootstrap.h"
30#include "GaudiKernel/ISvcLocator.h"
31#include "GaudiKernel/ITHistSvc.h"
32//#include "GaudiKernel/IHistogramSvc.h"
33#include "CLHEP/Vector/ThreeVector.h"
34#include "CLHEP/Vector/LorentzVector.h"
35#include "CLHEP/Vector/TwoVector.h"
36using CLHEP::Hep3Vector;
37using CLHEP::Hep2Vector;
38using CLHEP::HepLorentzVector;
39#include "CLHEP/Geometry/Point3D.h"
40#ifndef ENABLE_BACKWARDS_COMPATIBILITY
41 typedef HepGeom::Point3D<double> HepPoint3D;
42#endif
43#include "DQAJpsi2PPbarAlg/DQAJpsi2PPbarAlg.h"
44
45#include "VertexFit/KinematicFit.h"
46#include "VertexFit/VertexFit.h"
47#include "ParticleID/ParticleID.h"
48
49#include <vector>
50//const double twopi = 6.2831853;
51//const double pi = 3.1415927;
52const double mpi0 = 0.134977;
53const double mks0 = 0.497648;
54const double xmass[5] = {0.000511, 0.105658, 0.139570,0.493677, 0.938272};
55//const double velc = 29.9792458; tof_path unit in cm.
56const double velc = 299.792458; // tof path unit in mm
57typedef std::vector<int> Vint;
58typedef std::vector<HepLorentzVector> Vp4;
59
60//boosted
61//const HepLorentzVector p_cms(0.0407, 0.0, 0.0, 3.686);
62const HepLorentzVector p_cms(0.034067,0.0,0.0,3.097);
63const Hep3Vector u_cms = -p_cms.boostVector();
64
65
66//declare one counter
67static int counter[10]={0,0,0,0,0,0,0,0,0,0};
68/////////////////////////////////////////////////////////////////////////////
69
70DQAJpsi2PPbarAlg::DQAJpsi2PPbarAlg(const std::string& name, ISvcLocator* pSvcLocator) :
71 Algorithm(name, pSvcLocator) {
72
73 //Declare the properties
74 declareProperty("Vr0cut", m_vr0cut=5.0);
75 declareProperty("Vz0cut", m_vz0cut=20.0);
76 declareProperty("Vr1cut", m_vr1cut=1.0);
77 declareProperty("Vz1cut", m_vz1cut=10.0);
78 declareProperty("Vctcut", m_cthcut=0.93);
79 declareProperty("EnergyThreshold", m_energyThreshold=0.04);
80 declareProperty("GammaPhiCut", m_gammaPhiCut=20.0);
81 declareProperty("GammaThetaCut", m_gammaThetaCut=20.0);
82 declareProperty("Test4C", m_test4C = 1);
83 declareProperty("Test5C", m_test5C = 1);
84 declareProperty("CheckDedx", m_checkDedx = 1);
85 declareProperty("CheckTof", m_checkTof = 1);
86 declareProperty("GammaAngCut", m_gammaAngCut=20.0);
87}
88
89// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
90StatusCode DQAJpsi2PPbarAlg::initialize(){
91 MsgStream log(msgSvc(), name());
92
93 log << MSG::INFO << "in initialize()" << endmsg;
94
95 StatusCode status;
96
97
98
99 if(service("THistSvc", m_thsvc).isFailure()) {
100 log << MSG::ERROR << "Couldn't get THistSvc" << endreq;
101 return StatusCode::FAILURE;
102 }
103
104 // "DQAHist" is fixed
105 // "DQAJpsi2PPbar" is control sample name, just as ntuple case.
106 TH1F* hbst_p = new TH1F("bst_p", "bst_p", 80, 1.15, 1.31);
107 if(m_thsvc->regHist("/DQAHist/DQAJpsi2PPbar/hbst_p", hbst_p).isFailure()) {
108 log << MSG::ERROR << "Couldn't register bst_p" << endreq;
109 }
110
111 TH1F* hbst_cos = new TH1F("bst_cos", "bst_cos", 20, -1.0, 1.0);
112 if(m_thsvc->regHist("/DQAHist/DQAJpsi2PPbar/hbst_cos", hbst_cos).isFailure()) {
113 log << MSG::ERROR << "Couldn't register bst_cos" << endreq;
114 }
115
116 TH1F* hmpp = new TH1F("mpp", "mpp", 100, 3.05, 3.15);
117 if(m_thsvc->regHist("/DQAHist/DQAJpsi2PPbar/hmpp", hmpp).isFailure()) {
118 log << MSG::ERROR << "Couldn't register mpp" << endreq;
119 }
120
121 TH1F* hangle = new TH1F("angle", "angle", 50, 175.0, 180.);
122 if(m_thsvc->regHist("/DQAHist/DQAJpsi2PPbar/hangle", hangle).isFailure()) {
123 log << MSG::ERROR << "Couldn't register angle" << endreq;
124 }
125
126 TH1F* hdeltatof = new TH1F("deltatof", "deltatof", 50, -10., 10.);
127 if(m_thsvc->regHist("/DQAHist/DQAJpsi2PPbar/hdeltatof", hdeltatof).isFailure()) {
128 log << MSG::ERROR << "Couldn't register deltatof" << endreq;
129 }
130
131 TH1F* he_emc1 = new TH1F("e_emc1", "e_emc1", 100, 0.0, 2.0);
132 if(m_thsvc->regHist("/DQAHist/DQAJpsi2PPbar/he_emc1", he_emc1).isFailure()) {
133 log << MSG::ERROR << "Couldn't register e_emc1" << endreq;
134 }
135
136 TH1F* he_emc2 = new TH1F("e_emc2", "e_emc2", 100, 0.0, 2.0);
137 if(m_thsvc->regHist("/DQAHist/DQAJpsi2PPbar/he_emc2", he_emc2).isFailure()) {
138 log << MSG::ERROR << "Couldn't register e_emc2" << endreq;
139 }
140
141 // DQA
142 // The first directory specifier must be "DQAFILE"
143 // The second is the control sample name, the first letter is in upper format. eg. "DQAJpsi2PPbar"
144
145 NTuplePtr nt1(ntupleSvc(), "DQAFILE/DQAJpsi2PPbar");
146 if ( nt1 ) m_tuple = nt1;
147 else {
148 m_tuple = ntupleSvc()->book ("DQAFILE/DQAJpsi2PPbar", CLID_ColumnWiseTuple, "N-Tuple");
149 if ( m_tuple ) {
150 status = m_tuple->addItem ("runNo", m_runNo);
151 status = m_tuple->addItem ("event", m_event);
152 status = m_tuple->addItem ("Nchrg", m_nchrg);
153 status = m_tuple->addItem ("Nneu", m_nneu);
154 status = m_tuple->addItem ("ngch", m_ngch, 0, 10);
155
156 status = m_tuple->addIndexedItem ("charge",m_ngch, m_charge);
157 status = m_tuple->addIndexedItem ("vx0", m_ngch, m_vx0);
158 status = m_tuple->addIndexedItem ("vy0", m_ngch, m_vy0);
159 status = m_tuple->addIndexedItem ("vz0", m_ngch, m_vz0);
160 status = m_tuple->addIndexedItem ("vr0", m_ngch, m_vr0);
161
162 status = m_tuple->addIndexedItem ("vx", m_ngch, m_vx);
163 status = m_tuple->addIndexedItem ("vy", m_ngch, m_vy);
164 status = m_tuple->addIndexedItem ("vz", m_ngch, m_vz);
165 status = m_tuple->addIndexedItem ("vr", m_ngch, m_vr);
166
167 status = m_tuple->addIndexedItem ("px", m_ngch, m_px);
168 status = m_tuple->addIndexedItem ("py", m_ngch, m_py);
169 status = m_tuple->addIndexedItem ("pz", m_ngch, m_pz);
170 status = m_tuple->addIndexedItem ("p", m_ngch, m_p );
171 status = m_tuple->addIndexedItem ("cos", m_ngch, m_cos);
172
173 status = m_tuple->addIndexedItem ("bst_px", m_ngch, m_bst_px) ;
174 status = m_tuple->addIndexedItem ("bst_py", m_ngch, m_bst_py) ;
175 status = m_tuple->addIndexedItem ("bst_pz", m_ngch, m_bst_pz) ;
176 status = m_tuple->addIndexedItem ("bst_p", m_ngch, m_bst_p) ;
177 status = m_tuple->addIndexedItem ("bst_cos", m_ngch, m_bst_cos);
178
179
180 status = m_tuple->addIndexedItem ("chie" , m_ngch, m_chie) ;
181 status = m_tuple->addIndexedItem ("chimu" , m_ngch, m_chimu) ;
182 status = m_tuple->addIndexedItem ("chipi" , m_ngch, m_chipi) ;
183 status = m_tuple->addIndexedItem ("chik" , m_ngch, m_chik) ;
184 status = m_tuple->addIndexedItem ("chip" , m_ngch, m_chip) ;
185 status = m_tuple->addIndexedItem ("ghit" , m_ngch, m_ghit) ;
186 status = m_tuple->addIndexedItem ("thit" , m_ngch, m_thit) ;
187 status = m_tuple->addIndexedItem ("probPH" , m_ngch, m_probPH) ;
188 status = m_tuple->addIndexedItem ("normPH" , m_ngch, m_normPH) ;
189
190 status = m_tuple->addIndexedItem ("e_emc" , m_ngch, m_e_emc) ;
191
192
193 status = m_tuple->addIndexedItem ("qual_etof" , m_ngch, m_qual_etof );
194 status = m_tuple->addIndexedItem ("tof_etof" , m_ngch, m_tof_etof );
195 status = m_tuple->addIndexedItem ("te_etof" , m_ngch, m_te_etof );
196 status = m_tuple->addIndexedItem ("tmu_etof" , m_ngch, m_tmu_etof );
197 status = m_tuple->addIndexedItem ("tpi_etof" , m_ngch, m_tpi_etof );
198 status = m_tuple->addIndexedItem ("tk_etof" , m_ngch, m_tk_etof );
199 status = m_tuple->addIndexedItem ("tp_etof" , m_ngch, m_tp_etof );
200
201 status = m_tuple->addIndexedItem ("qual_btof1", m_ngch, m_qual_btof1 );
202 status = m_tuple->addIndexedItem ("tof_btof1" , m_ngch, m_tof_btof1 );
203 status = m_tuple->addIndexedItem ("te_btof1" , m_ngch, m_te_btof1 );
204 status = m_tuple->addIndexedItem ("tmu_btof1" , m_ngch, m_tmu_btof1 );
205 status = m_tuple->addIndexedItem ("tpi_btof1" , m_ngch, m_tpi_btof1 );
206 status = m_tuple->addIndexedItem ("tk_btof1" , m_ngch, m_tk_btof1 );
207 status = m_tuple->addIndexedItem ("tp_btof1" , m_ngch, m_tp_btof1 );
208
209 status = m_tuple->addIndexedItem ("qual_btof2", m_ngch, m_qual_btof2 );
210 status = m_tuple->addIndexedItem ("tof_btof2" , m_ngch, m_tof_btof2 );
211 status = m_tuple->addIndexedItem ("te_btof2" , m_ngch, m_te_btof2 );
212 status = m_tuple->addIndexedItem ("tmu_btof2" , m_ngch, m_tmu_btof2 );
213 status = m_tuple->addIndexedItem ("tpi_btof2" , m_ngch, m_tpi_btof2 );
214 status = m_tuple->addIndexedItem ("tk_btof2" , m_ngch, m_tk_btof2 );
215 status = m_tuple->addIndexedItem ("tp_btof2" , m_ngch, m_tp_btof2 );
216
217 status = m_tuple->addIndexedItem ("dedx_pid", m_ngch, m_dedx_pid);
218 status = m_tuple->addIndexedItem ("tof1_pid", m_ngch, m_tof1_pid);
219 status = m_tuple->addIndexedItem ("tof2_pid", m_ngch, m_tof2_pid);
220 status = m_tuple->addIndexedItem ("prob_pi", m_ngch, m_prob_pi );
221 status = m_tuple->addIndexedItem ("prob_k", m_ngch, m_prob_k );
222 status = m_tuple->addIndexedItem ("prob_p", m_ngch, m_prob_p );
223
224 status = m_tuple->addItem ( "np", m_np );
225 status = m_tuple->addItem ( "npb", m_npb );
226
227 status = m_tuple->addItem ( "m2p", m_m2p );
228 status = m_tuple->addItem ( "angle", m_angle );
229 status = m_tuple->addItem ( "deltatof", m_deltatof );
230
231 status = m_tuple->addItem ( "vtx_m2p", m_vtx_m2p );
232 status = m_tuple->addItem ( "vtx_angle", m_vtx_angle );
233
234 status = m_tuple->addItem ( "m_chi2_4c", m_chi2_4c );
235 status = m_tuple->addItem ( "m_m2p_4c", m_m2p_4c );
236 status = m_tuple->addItem ( "m_angle_4c", m_angle_4c );
237
238 }
239 else {
240 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple) << endmsg;
241 return StatusCode::FAILURE;
242 }
243 }
244
245 //
246 //--------end of book--------
247 //
248
249 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
250 return StatusCode::SUCCESS;
251
252}
253
254// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
255StatusCode DQAJpsi2PPbarAlg::execute() {
256
257
258 MsgStream log(msgSvc(), name());
259 log << MSG::INFO << "in execute()" << endreq;
260
261 setFilterPassed(false);
262
263 counter[0]++;
264 log << MSG::INFO << "counter[0]=" << counter[0]<< endreq;
265
266 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
267 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
268 m_runNo = eventHeader->runNumber();
269 m_event = eventHeader->eventNumber();
270 m_nchrg = evtRecEvent->totalCharged();
271 m_nneu = evtRecEvent->totalNeutral();
272
273 log << MSG::INFO <<"ncharg, nneu, tottks = "
274 << evtRecEvent->totalCharged() << " , "
275 << evtRecEvent->totalNeutral() << " , "
276 << evtRecEvent->totalTracks() <<endreq;
277
278 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
279 //
280 // check x0, y0, z0, r0
281 // suggest cut: |z0|<10 && r0<1
282 //
283
284 Vint iGood, iptrk, imtrk;
285 iGood.clear();
286 iptrk.clear();
287 imtrk.clear();
288 int nCharge = 0;
289 Hep3Vector xorigin(0,0,0);
290
291 //if (m_reader.isRunNumberValid(runNo)) {
292 IVertexDbSvc* vtxsvc;
293 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
294 if (vtxsvc->isVertexValid()){
295 double* dbv = vtxsvc->PrimaryVertex();
296 double* vv = vtxsvc->SigmaPrimaryVertex();
297 // HepVector dbv = m_reader.PrimaryVertex(runNo);
298 // HepVector vv = m_reader.SigmaPrimaryVertex(runNo);
299 xorigin.setX(dbv[0]);
300 xorigin.setY(dbv[1]);
301 xorigin.setZ(dbv[2]);
302 log << MSG::INFO
303 <<"xorigin.x="<<xorigin.x()<<", "
304 <<"xorigin.y="<<xorigin.y()<<", "
305 <<"xorigin.z="<<xorigin.z()<<", "
306 <<endreq ;
307 }
308
309 for (int i = 0; i < evtRecEvent->totalCharged(); i++){
310 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
311 if(!(*itTrk)->isMdcTrackValid()) continue;
312 if(!(*itTrk)->isMdcKalTrackValid()) continue;
313 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
314 double x0=mdcTrk->x();
315 double y0=mdcTrk->y();
316 double z0=mdcTrk->z();
317 double phi0=mdcTrk->helix(1);
318 double xv=xorigin.x();
319 double yv=xorigin.y();
320 double zv=xorigin.z();
321 double rv=(x0-xv)*cos(phi0)+(y0-yv)*sin(phi0);
322 double cost = cos(mdcTrk->theta());
323 if(fabs(z0-zv) >= m_vz1cut) continue;
324 if(fabs(rv) >= m_vr1cut) continue;
325 //if(fabs(cost) >= m_cthcut ) continue;
326
327 iGood.push_back(i);
328 nCharge += mdcTrk->charge();
329
330 if (mdcTrk->charge() > 0) {
331 iptrk.push_back(i);
332 } else {
333 imtrk.push_back(i);
334 }
335 }
336 //
337 // Finish Good Charged Track Selection
338 //
339 int nGood = iGood.size();
340 m_ngch = iGood.size();
341 log << MSG::INFO << "ngood, totcharge = " << nGood << " , " << nCharge << endreq;
342 if((nGood != 2) || (nCharge!=0) ){
343 return StatusCode::SUCCESS;
344 }
345 counter[1]++;
346 log << MSG::INFO << "counter[1]=" << counter[1]<< endreq;
347
348 /////////////////////////////////////////////////
349 // PID
350 /////////////////////////////////////////////////
351
352 Vint ipp, ipm;
353 ipp.clear();
354 ipm.clear();
355
356 Vp4 p_pp, p_pm;
357 p_pp.clear();
358 p_pm.clear();
359
360 int ii=-1 ;
361
362 ParticleID *pid = ParticleID::instance();
363 for(int i = 0; i < nGood; i++) {
364
365 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
366 if(!(*itTrk)->isMdcTrackValid()) continue;
367 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
368 if(!(*itTrk)->isMdcKalTrackValid()) continue;
369 RecMdcKalTrack* mdcKalTrk = (*itTrk)->mdcKalTrack();
370
371 // if(pid) delete pid;
372 pid->init();
373 pid->setMethod(pid->methodProbability());
374 pid->setChiMinCut(4);
375 pid->setRecTrack(*itTrk);
376 // pid->usePidSys(pid->useDedx() | pid->useTof1() | pid->useTof2() | pid->useTofE() | pid->useTofQ());
377 // use PID sub-system
378 pid->usePidSys(pid->useDedx() | pid->useTof1() | pid->useTof2()); // use PID sub-system
379 // pid->identify(pid->onlyProton());
380 pid->identify(pid->onlyPionKaonProton()); // seperater Pion/Kaon/Proton
381 // pid->identify(pid->onlyPion() | pid->onlyKaon()); // seperater Pion/Kaon
382 // pid->identify(pid->onlyPion());
383 // pid->identify(pid->onlyKaon());
384 pid->calculate();
385 if(!(pid->IsPidInfoValid())) continue;
386
387 double prob_pi = pid->probPion();
388 double prob_k = pid->probKaon();
389 double prob_p = pid->probProton();
390
391 // if(pid->probPion() < 0.001 || (pid->probPion() < pid->probKaon())) continue;
392 // if(pid->probPion() < 0.001) continue;
393 if (prob_p > prob_pi && prob_p > prob_k) {
394
395 HepLorentzVector pkaltrk;
396
397 mdcKalTrk->setPidType(RecMdcKalTrack::proton);
398 pkaltrk.setPx(mdcKalTrk->px());
399 pkaltrk.setPy(mdcKalTrk->py());
400 pkaltrk.setPz(mdcKalTrk->pz());
401 double p3 = pkaltrk.mag();
402 pkaltrk.setE(sqrt(p3*p3+xmass[4]*xmass[4]));
403
404 if(mdcTrk->charge() >0) {
405 ii = 0;
406 ipp.push_back(iGood[i]);
407 p_pp.push_back(pkaltrk);
408 } else {
409 ii = 1 ;
410 ipm.push_back(iGood[i]);
411 p_pm.push_back(pkaltrk);
412 }
413
414 m_dedx_pid[ii] = pid->chiDedx(2);
415 m_tof1_pid[ii] = pid->chiTof1(2);
416 m_tof2_pid[ii] = pid->chiTof2(2);
417 m_prob_pi[ii] = pid->probPion();
418 m_prob_k[ii] = pid->probKaon();
419 m_prob_p[ii] = pid->probProton();
420
421 }
422 } // with PID
423
424 m_np = ipp.size();
425 m_npb = ipm.size();
426
427 //if(m_np*m_npb != 1) return SUCCESS;
428
429 counter[2]++;
430 log << MSG::INFO << "counter[2]=" << counter[2]<< endreq;
431
432 ///////////////////////////////////////////////
433 // read information for good charged tracks
434 ///////////////////////////////////////////////
435 Vp4 p_ptrk, p_mtrk;
436 p_ptrk.clear(), p_mtrk.clear();
437 RecMdcKalTrack *ppKalTrk = 0 , *pmKalTrk = 0 ;
438
439 ii = -1 ;
440 for(int i = 0; i < m_ngch; i++ ){
441 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + iGood[i];
442 if(!(*itTrk)->isMdcTrackValid()) continue;
443 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
444 if(!(*itTrk)->isMdcKalTrackValid()) continue;
445 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
446
447 mdcKalTrk->setPidType(RecMdcKalTrack::proton);
448
449 if (mdcTrk->charge() > 0 ) {
450 ii = 0 ;
451 ppKalTrk = mdcKalTrk ;
452 }else{
453 ii = 1 ;
454 pmKalTrk = mdcKalTrk ;
455 }
456
457 m_charge[ii] = mdcTrk->charge();
458
459 double x0=mdcKalTrk->x();
460 double y0=mdcKalTrk->y();
461 double z0=mdcKalTrk->z();
462 double phi0=mdcTrk->helix(1);
463 double xv=xorigin.x();
464 double yv=xorigin.y();
465 double zv=xorigin.z();
466 double rv=(x0-xv)*cos(phi0)+(y0-yv)*sin(phi0);
467
468 m_vx0[ii] = x0 ;
469 m_vy0[ii] = y0 ;
470 m_vz0[ii] = z0 ;
471 m_vr0[ii] = sqrt((x0*x0)+(y0*y0)) ;
472
473 m_vx[ii] = x0-xv ;
474 m_vy[ii] = y0-yv ;
475 m_vz[ii] = fabs(z0-zv) ;
476 m_vr[ii] = fabs(rv) ;
477
478 mdcKalTrk->setPidType(RecMdcKalTrack::proton);
479 m_px[ii] = mdcKalTrk->px();
480 m_py[ii] = mdcKalTrk->py();
481 m_pz[ii] = mdcKalTrk->pz();
482 m_p[ii] = mdcKalTrk->p();
483 m_cos[ii] = mdcKalTrk->pz()/mdcKalTrk->p();
484 double temp_e = sqrt(m_p[ii]*m_p[ii] + xmass[4]*xmass[4]);
485 HepLorentzVector temp_p(m_px[ii],m_py[ii],m_pz[ii],temp_e);
486 if (i==0){
487 p_ptrk.push_back(temp_p);
488 }else{
489 p_mtrk.push_back(temp_p);
490 }
491
492
493 double ptrk = mdcKalTrk->p() ;
494 if ((*itTrk)->isMdcDedxValid()) { //Dedx information
495 RecMdcDedx* dedxTrk = (*itTrk)->mdcDedx();
496 m_chie[ii] = dedxTrk->chiE();
497 m_chimu[ii] = dedxTrk->chiMu();
498 m_chipi[ii] = dedxTrk->chiPi();
499 m_chik[ii] = dedxTrk->chiK();
500 m_chip[ii] = dedxTrk->chiP();
501 m_ghit[ii] = dedxTrk->numGoodHits();
502 m_thit[ii] = dedxTrk->numTotalHits();
503 m_probPH[ii]= dedxTrk->probPH();
504 m_normPH[ii]= dedxTrk->normPH();
505
506 }
507
508 if((*itTrk)->isEmcShowerValid()) {
509 RecEmcShower *emcTrk = (*itTrk)->emcShower();
510 m_e_emc[ii] = emcTrk->energy();
511 }
512
513
514 if((*itTrk)->isTofTrackValid()) {
515
516 SmartRefVector<RecTofTrack> tofTrkCol = (*itTrk)->tofTrack();
517
518 SmartRefVector<RecTofTrack>::iterator iter_tof = tofTrkCol.begin();
519 for(;iter_tof != tofTrkCol.end(); iter_tof++ ) {
520 TofHitStatus *status = new TofHitStatus;
521 status->setStatus((*iter_tof)->status());
522
523 if(!(status->is_barrel())){//endcap
524 if( !(status->is_counter()) ) continue; // ?
525 if( status->layer()!=0 ) continue;//layer1
526 double path=(*iter_tof)->path(); // ?
527 double tof = (*iter_tof)->tof();
528 double ph = (*iter_tof)->ph();
529 double rhit = (*iter_tof)->zrhit();
530 double qual = 0.0 + (*iter_tof)->quality();
531 double cntr = 0.0 + (*iter_tof)->tofID();
532 double texp[5];
533 for(int j = 0; j < 5; j++) {
534 double gb = ptrk/xmass[j];
535 double beta = gb/sqrt(1+gb*gb);
536 texp[j] = path /beta/velc;
537 }
538
539 m_qual_etof[ii] = qual;
540 m_tof_etof[ii] = tof ;
541 m_te_etof[ii] = tof - texp[0];
542 m_tmu_etof[ii] = tof - texp[1];
543 m_tpi_etof[ii] = tof - texp[2];
544 m_tk_etof[ii] = tof - texp[3];
545 m_tp_etof[ii] = tof - texp[4];
546 }
547 else {//barrel
548 if( !(status->is_counter()) ) continue; // ?
549 if(status->layer()==1){ //layer1
550 double path=(*iter_tof)->path(); // ?
551 double tof = (*iter_tof)->tof();
552 double ph = (*iter_tof)->ph();
553 double rhit = (*iter_tof)->zrhit();
554 double qual = 0.0 + (*iter_tof)->quality();
555 double cntr = 0.0 + (*iter_tof)->tofID();
556 double texp[5];
557 for(int j = 0; j < 5; j++) {
558 double gb = ptrk/xmass[j];
559 double beta = gb/sqrt(1+gb*gb);
560 texp[j] = path /beta/velc;
561 }
562
563 m_qual_btof1[ii] = qual;
564 m_tof_btof1[ii] = tof ;
565 m_te_btof1[ii] = tof - texp[0];
566 m_tmu_btof1[ii] = tof - texp[1];
567 m_tpi_btof1[ii] = tof - texp[2];
568 m_tk_btof1[ii] = tof - texp[3];
569 m_tp_btof1[ii] = tof - texp[4];
570
571 }
572
573 if(status->layer()==2){//layer2
574 double path=(*iter_tof)->path(); // ?
575 double tof = (*iter_tof)->tof();
576 double ph = (*iter_tof)->ph();
577 double rhit = (*iter_tof)->zrhit();
578 double qual = 0.0 + (*iter_tof)->quality();
579 double cntr = 0.0 + (*iter_tof)->tofID();
580 double texp[5];
581 for(int j = 0; j < 5; j++) {
582 double gb = ptrk/xmass[j];
583 double beta = gb/sqrt(1+gb*gb);
584 texp[j] = path /beta/velc;
585 }
586
587 m_qual_btof2[ii] = qual;
588 m_tof_btof2[ii] = tof ;
589 m_te_btof2[ii] = tof - texp[0];
590 m_tmu_btof2[ii] = tof - texp[1];
591 m_tpi_btof2[ii] = tof - texp[2];
592 m_tk_btof2[ii] = tof - texp[3];
593 m_tp_btof2[ii] = tof - texp[4];
594 }
595 }
596 delete status;
597 }
598 }
599 }
600 counter[3]++;
601 log << MSG::INFO << "counter[3]=" << counter[3]<< endreq;
602
603
604 //boosted mdcKalTrk
605 //cout <<"before p_ptrk[0]="<<p_ptrk[0]<<endl;
606 //cout <<"before p_mtrk[0]="<<p_mtrk[0]<<endl;
607 //cout <<"_m2p = "<<(p_ptrk[0] + p_mtrk[0]).m()<<endl ;
608 p_ptrk[0].boost(u_cms);
609 p_mtrk[0].boost(u_cms);
610 for (int i=0; i<=1; i++ ) {
611 HepLorentzVector p;
612 if (i==0) p = p_ptrk[0];
613 if (i==1) p = p_mtrk[0];
614
615 m_bst_px[i] = p.px();
616 m_bst_py[i] = p.py();
617 m_bst_pz[i] = p.pz();
618 m_bst_p[i] = p.rho();
619 m_bst_cos[i]= p.cosTheta();
620 }
621
622 m_m2p = (p_ptrk[0] + p_mtrk[0]).m();
623 //cout <<"after p_ptrk[0]="<<p_ptrk[0]<<endl;
624 //cout <<"after p_mtrk[0]="<<p_mtrk[0]<<endl;
625 //cout <<"_m2p = "<<(p_ptrk[0] + p_mtrk[0]).m()<<endl ;
626 m_angle = (p_ptrk[0].vect()).angle((p_mtrk[0]).vect()) * 180.0/(CLHEP::pi) ;
627
628 //cout <<"m_angle="<<m_angle<<endl;
629 //cout <<"m_e_emc="<<m_e_emc[0]<<endl;
630
631 double deltatof = 0.0 ;
632 if(m_tof_btof1[0]*m_tof_btof1[1] != 0.0) deltatof+=m_tof_btof1[0]-m_tof_btof1[1];
633 if(m_tof_btof2[0]*m_tof_btof2[1] != 0.0) deltatof+=m_tof_btof2[0]-m_tof_btof2[1];
634 m_deltatof = deltatof ;
635
636
637
638 // Vertex Fit
639
640 // Kinamatic Fit
641
642 // finale selection
643 if (fabs(m_bst_p[0]-1.232)>0.02 ) {return StatusCode::SUCCESS ; }
644 if (fabs(m_bst_p[1]-1.232)>0.02 ) {return StatusCode::SUCCESS ; }
645 if (fabs(m_deltatof)>4.0) {return StatusCode::SUCCESS ; }
646 if (m_angle<178.0) {return StatusCode::SUCCESS ; }
647 if (m_e_emc[0]>0.7) {return StatusCode::SUCCESS ; }
648
649 counter[4]++;
650 log << MSG::INFO << "counter[4]=" << counter[4]<< endreq;
651
652 // DQA
653
654 (*(evtRecTrkCol->begin()+iptrk[0]))->tagProton();
655 (*(evtRecTrkCol->begin()+imtrk[0]))->tagProton();
656
657 // Quality: defined by whether dE/dx or TOF is used to identify particle
658 // 0 - no dE/dx, no TOF (can be used for dE/dx and TOF calibration)
659 // 1 - only dE/dx (can be used for TOF calibration)
660 // 2 - only TOF (can be used for dE/dx calibration)
661 // 3 - Both dE/dx and TOF
662 (*(evtRecTrkCol->begin()+iptrk[0]))->setQuality(0);
663 (*(evtRecTrkCol->begin()+imtrk[0]))->setQuality(0);
664
665 setFilterPassed(true);
666
667 TH1 *h(0);
668 if (m_thsvc->getHist("/DQAHist/DQAJpsi2PPbar/hbst_p", h).isSuccess()) {
669 h->Fill(m_bst_p[0]);
670 } else {
671 log << MSG::ERROR << "Couldn't retrieve hbst_p" << endreq;
672 }
673
674 if (m_thsvc->getHist("/DQAHist/DQAJpsi2PPbar/hbst_cos", h).isSuccess()) {
675 h->Fill(m_bst_cos[0]);
676 } else {
677 log << MSG::ERROR << "Couldn't retrieve hbst_cos" << endreq;
678 }
679
680 if (m_thsvc->getHist("/DQAHist/DQAJpsi2PPbar/hmpp", h).isSuccess()) {
681 h->Fill(m_m2p);
682 } else {
683 log << MSG::ERROR << "Couldn't retrieve hmpp" << endreq;
684 }
685
686 if (m_thsvc->getHist("/DQAHist/DQAJpsi2PPbar/hangle", h).isSuccess()) {
687 h->Fill(m_angle);
688 } else {
689 log << MSG::ERROR << "Couldn't retrieve hangle" << endreq;
690 }
691
692 if (m_thsvc->getHist("/DQAHist/DQAJpsi2PPbar/hdeltatof", h).isSuccess()) {
693 h->Fill(m_deltatof);
694 } else {
695 log << MSG::ERROR << "Couldn't retrieve hdeltatof" << endreq;
696 }
697
698 if (m_thsvc->getHist("/DQAHist/DQAJpsi2PPbar/he_emc1", h).isSuccess()) {
699 h->Fill(m_e_emc[0]);
700 } else {
701 log << MSG::ERROR << "Couldn't retrieve he_emc1" << endreq;
702 }
703
704 if (m_thsvc->getHist("/DQAHist/DQAJpsi2PPbar/he_emc2", h).isSuccess()) {
705 h->Fill(m_e_emc[1]);
706 } else {
707 log << MSG::ERROR << "Couldn't retrieve he_emc2" << endreq;
708 }
709
710
711
712
713
714 m_tuple -> write();
715
716 return StatusCode::SUCCESS;
717}
718
719
720// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
721StatusCode DQAJpsi2PPbarAlg::finalize() {
722
723 MsgStream log(msgSvc(), name());
724 log << MSG::INFO << "in finalize()" << endmsg;
725
726 std::cout << "*********** Singal.cxx *****************" << std::endl;
727 std::cout << "the totale events is " << counter[0] << std::endl;
728 std::cout << "select good charged track " << counter[1] << std::endl;
729 std::cout << "PID " << counter[2] << std::endl;
730 std::cout << "inform. for charged track " << counter[3] << std::endl;
731 std::cout << "after all selections " << counter[4] << std::endl;
732 std::cout << "****************************************" << std::endl;
733
734 return StatusCode::SUCCESS;
735}
736
u_cms
const Hep3Vector u_cms
Definition: DQADtagAlg.cxx:62
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition: DQAJpsi2PPbarAlg.cxx:41
mks0
const double mks0
Definition: DQAJpsi2PPbarAlg.cxx:53
u_cms
const Hep3Vector u_cms
Definition: DQAJpsi2PPbarAlg.cxx:63
Vp4
std::vector< HepLorentzVector > Vp4
Definition: DQAJpsi2PPbarAlg.cxx:58
xmass
const double xmass[5]
Definition: DQAJpsi2PPbarAlg.cxx:54
velc
const double velc
Definition: DQAJpsi2PPbarAlg.cxx:56
mpi0
const double mpi0
Definition: DQAJpsi2PPbarAlg.cxx:52
Vint
std::vector< int > Vint
Definition: DQAJpsi2PPbarAlg.cxx:57
p_cms
const HepLorentzVector p_cms(0.034067, 0.0, 0.0, 3.097)
Vp4
std::vector< HepLorentzVector > Vp4
Definition: Gam4pikp.cxx:53
xmass
const double xmass[5]
Definition: Gam4pikp.cxx:50
velc
const double velc
Definition: Gam4pikp.cxx:51
Vint
std::vector< int > Vint
Definition: Gam4pikp.cxx:52
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
DQAJpsi2PPbarAlg::initialize
StatusCode initialize()
Definition: DQAJpsi2PPbarAlg.cxx:90
DQAJpsi2PPbarAlg::DQAJpsi2PPbarAlg
DQAJpsi2PPbarAlg(const std::string &name, ISvcLocator *pSvcLocator)
Definition: DQAJpsi2PPbarAlg.cxx:70
DstMdcTrack::helix
const HepVector helix() const
......
Definition: DstMdcTrack.cxx:151
IVertexDbSvc::isVertexValid
virtual bool isVertexValid()=0
IVertexDbSvc::SigmaPrimaryVertex
virtual double * SigmaPrimaryVertex()=0
IVertexDbSvc::PrimaryVertex
virtual double * PrimaryVertex()=0
ParticleID::chiTof2
double chiTof2(int n) const
ParticleID::instance
static ParticleID * instance()
Definition: ParticleID.cxx:22
ParticleID::IsPidInfoValid
bool IsPidInfoValid() const
ParticleID::chiTof1
double chiTof1(int n) const
ParticleID::calculate
void calculate()
Definition: ParticleID.cxx:97
ParticleID::init
void init()
Definition: ParticleID.cxx:27
ParticleID::chiDedx
double chiDedx(int n) const
TofHitStatus::setStatus
void setStatus(unsigned int status)
Definition: TofHitStatus.cxx:134