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