BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
DQARhopi.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/MsgStream.h"
8#include "GaudiKernel/NTuple.h"
9#include "GaudiKernel/SmartDataPtr.h"
10#include "TH1F.h"
11#include "TMath.h"
12#include <vector>
13
14#include "DstEvent/TofHitStatus.h"
15#include "EventModel/Event.h"
16#include "EventModel/EventHeader.h"
17#include "EventModel/EventModel.h"
18#include "EvtRecEvent/EvtRecEvent.h"
19#include "EvtRecEvent/EvtRecTrack.h"
20#include "McTruth/McParticle.h"
21#include "ParticleID/ParticleID.h"
22#include "VertexDbSvc/IVertexDbSvc.h"
23#include "VertexFit/Helix.h"
24#include "VertexFit/KinematicFit.h"
25#include "VertexFit/SecondVertexFit.h"
26#include "VertexFit/VertexFit.h"
27
28using CLHEP::Hep2Vector;
29using CLHEP::Hep3Vector;
30using CLHEP::HepLorentzVector;
31#include "CLHEP/Geometry/Point3D.h"
32#ifndef ENABLE_BACKWARDS_COMPATIBILITY
33typedef HepGeom::Point3D<double> HepPoint3D;
34#endif
35
36#include "DQARhopi.h"
37
38using namespace Event;
39// const double twopi = 6.2831853;
40// const double pi = 3.1415927;
41const double mpi = 0.13957;
42const double mk = 0.493677;
43const double xmass[5] = { 0.000511, 0.105658, 0.139570, 0.493677, 0.938272 };
44// const double velc = 29.9792458; tof_path unit in cm.
45const double velc = 299.792458; // tof path unit in mm
46typedef std::vector<int> Vint;
47typedef std::vector<HepLorentzVector> Vp4;
48
49const HepLorentzVector ecms( 0.034, 0, 0, 3.097 );
50const Hep3Vector u_cms = -ecms.boostVector();
51
52/////////////////////////////////////////////////////////////////////////////
53
54DQARhopi::DQARhopi( const std::string& name, ISvcLocator* pSvcLocator )
55 : Algorithm( name, pSvcLocator ) {
56
57 // Declare the properties
58 declareProperty( "Vr0cut", m_vr0cut = 1.0 );
59 declareProperty( "Vz0cut", m_vz0cut = 5.0 );
60 declareProperty( "Vctcut", m_cthcut = 0.93 );
61 declareProperty( "EnergyThreshold", m_energyThreshold = 0.05 );
62 declareProperty( "GammaAngCut", m_gammaAngCut = 25.0 );
63 declareProperty( "Test4C", m_test4C = 1 );
64 declareProperty( "Test5C", m_test5C = 1 );
65 declareProperty( "CheckDedx", m_checkDedx = 1 );
66 declareProperty( "CheckTof", m_checkTof = 1 );
67}
68
69// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
70StatusCode DQARhopi::initialize() {
71 MsgStream log( msgSvc(), name() );
72
73 log << MSG::INFO << "in initialize()" << endmsg;
74
75 Ncut0 = Ncut1 = Ncut2 = Ncut3 = Ncut4 = Ncut5 = Ncut6 = Ncut7 = Ncut8 = Ncut9 = Ncut10 = 0;
76
77 StatusCode status;
78
79 NTuplePtr nt4( ntupleSvc(), "DQAFILE/Rhopi" );
80 if ( nt4 ) m_tuple4 = nt4;
81 else
82 {
83 m_tuple4 =
84 ntupleSvc()->book( "DQAFILE/Rhopi", CLID_ColumnWiseTuple, "ks N-Tuple example" );
85 if ( m_tuple4 )
86 {
87 status = m_tuple4->addItem( "run", m_run );
88 status = m_tuple4->addItem( "rec", m_rec );
89 status = m_tuple4->addItem( "nch", m_nch );
90 status = m_tuple4->addItem( "nneu", m_nneu );
91 status = m_tuple4->addItem( "chi1", m_chi1 );
92 status = m_tuple4->addItem( "mpi0", m_mpi0 );
93 status = m_tuple4->addItem( "mprho0", m_prho0 );
94 status = m_tuple4->addItem( "mprhop", m_prhop );
95 status = m_tuple4->addItem( "mprhom", m_prhom );
96 status = m_tuple4->addItem( "mgood", m_good );
97 status = m_tuple4->addItem( "mgam", m_gam );
98 status = m_tuple4->addItem( "mpip", m_pip );
99 status = m_tuple4->addItem( "mpim", m_pim );
100 status = m_tuple4->addItem( "m2gam", m_2gam );
101 status = m_tuple4->addItem( "ngch", m_ngch, 0, 10 );
102 status = m_tuple4->addItem( "nggneu", m_nggneu, 0, 10 );
103 status = m_tuple4->addItem( "moutpi0", m_outpi0 );
104 status = m_tuple4->addItem( "cosang", m_cosang );
105 status = m_tuple4->addItem( "moutpip", m_outpip );
106 status = m_tuple4->addItem( "moutpim", m_outpim );
107 status = m_tuple4->addItem( "menpip", m_enpip );
108 status = m_tuple4->addItem( "mdcpip", m_dcpip );
109 status = m_tuple4->addItem( "menpim", m_enpim );
110 status = m_tuple4->addItem( "mdcpim", m_dcpim );
111 status = m_tuple4->addItem( "mpipf", m_pipf );
112 status = m_tuple4->addItem( "mpimf", m_pimf );
113 status = m_tuple4->addItem( "mpi0f", m_pi0f );
114
115 status = m_tuple4->addItem( "mpmax", m_pmax );
116 status = m_tuple4->addItem( "mppx", m_ppx );
117 status = m_tuple4->addItem( "mppy", m_ppy );
118 status = m_tuple4->addItem( "mppz", m_ppz );
119 status = m_tuple4->addItem( "mcosthep", m_costhep );
120 status = m_tuple4->addItem( "mppxkal", m_ppxkal );
121 status = m_tuple4->addItem( "mppykal", m_ppykal );
122 status = m_tuple4->addItem( "mppzkal", m_ppzkal );
123 status = m_tuple4->addItem( "mmpx", m_mpx );
124 status = m_tuple4->addItem( "mmpy", m_mpy );
125 status = m_tuple4->addItem( "mmpz", m_mpz );
126 status = m_tuple4->addItem( "mcosthem", m_costhem );
127 status = m_tuple4->addItem( "mmpxkal", m_mpxkal );
128 status = m_tuple4->addItem( "mmpykal", m_mpykal );
129 status = m_tuple4->addItem( "mmpzkal", m_mpzkal );
130
131 status = m_tuple4->addItem( "mvxpin", m_vxpin );
132 status = m_tuple4->addItem( "mvypin", m_vypin );
133 status = m_tuple4->addItem( "mvzpin", m_vzpin );
134 status = m_tuple4->addItem( "mvrpin", m_vrpin );
135 status = m_tuple4->addItem( "mcosthepin", m_costhepin );
136 status = m_tuple4->addItem( "mvxmin", m_vxmin );
137 status = m_tuple4->addItem( "mvymin", m_vymin );
138 status = m_tuple4->addItem( "mvzmin", m_vzmin );
139 status = m_tuple4->addItem( "mvrmin", m_vrmin );
140 status = m_tuple4->addItem( "mcosthemin", m_costhemin );
141 status = m_tuple4->addItem( "mvxp", m_vxp );
142 status = m_tuple4->addItem( "mvyp", m_vyp );
143 status = m_tuple4->addItem( "mvzp", m_vzp );
144 status = m_tuple4->addItem( "mvrp", m_vrp );
145 status = m_tuple4->addItem( "mvxm", m_vxm );
146 status = m_tuple4->addItem( "mvym", m_vym );
147 status = m_tuple4->addItem( "mvzm", m_vzm );
148 status = m_tuple4->addItem( "mvrm", m_vrm );
149 status = m_tuple4->addItem( "nangecc", m_nangecc, 0, 10 );
150 status = m_tuple4->addIndexedItem( "mdthec", m_nangecc, m_dthec );
151 status = m_tuple4->addIndexedItem( "mdphic", m_nangecc, m_dphic );
152 status = m_tuple4->addIndexedItem( "mdangc", m_nangecc, m_dangc );
153 status = m_tuple4->addIndexedItem( "mspippim", m_nangecc, m_mspippim );
154
155 status = m_tuple4->addItem( "angsg", dangsg );
156 status = m_tuple4->addItem( "thesg", dthesg );
157 status = m_tuple4->addItem( "phisg", dphisg );
158 status = m_tuple4->addItem( "cosgth1", cosgth1 );
159 status = m_tuple4->addItem( "cosgth2", cosgth2 );
160
161 status = m_tuple4->addItem( "mchi5", m_chi5 );
162 status = m_tuple4->addItem( "mkpi0", m_kpi0 );
163 status = m_tuple4->addItem( "mkpkm", m_kpkm );
164 status = m_tuple4->addItem( "mkpp0", m_kpp0 );
165 status = m_tuple4->addItem( "mkmp0", m_kmp0 );
166 status = m_tuple4->addItem( "mpgam2pi1", m_pgam2pi1 );
167 status = m_tuple4->addItem( "mpgam2pi2", m_pgam2pi2 );
168 status = m_tuple4->addItem( "cosva1", cosva1 );
169 status = m_tuple4->addItem( "cosva2", cosva2 );
170 status = m_tuple4->addItem( "laypi1", m_laypi1 );
171 status = m_tuple4->addItem( "hit1", m_hit1 );
172 status = m_tuple4->addItem( "laypi2", m_laypi2 );
173 status = m_tuple4->addItem( "hit2", m_hit2 );
174 status = m_tuple4->addItem( "anglepm", m_anglepm );
175
176 status = m_tuple4->addIndexedItem( "mptrk", m_ngch, m_ptrk );
177 status = m_tuple4->addIndexedItem( "chie", m_ngch, m_chie );
178 status = m_tuple4->addIndexedItem( "chimu", m_ngch, m_chimu );
179 status = m_tuple4->addIndexedItem( "chipi", m_ngch, m_chipi );
180 status = m_tuple4->addIndexedItem( "chik", m_ngch, m_chik );
181 status = m_tuple4->addIndexedItem( "chip", m_ngch, m_chip );
182 status = m_tuple4->addIndexedItem( "probPH", m_ngch, m_probPH );
183 status = m_tuple4->addIndexedItem( "normPH", m_ngch, m_normPH );
184 status = m_tuple4->addIndexedItem( "ghit", m_ngch, m_ghit );
185 status = m_tuple4->addIndexedItem( "thit", m_ngch, m_thit );
186
187 status = m_tuple4->addIndexedItem( "ptot_etof", m_ngch, m_ptot_etof );
188 status = m_tuple4->addIndexedItem( "cntr_etof", m_ngch, m_cntr_etof );
189 status = m_tuple4->addIndexedItem( "te_etof", m_ngch, m_te_etof );
190 status = m_tuple4->addIndexedItem( "tmu_etof", m_ngch, m_tmu_etof );
191 status = m_tuple4->addIndexedItem( "tpi_etof", m_ngch, m_tpi_etof );
192 status = m_tuple4->addIndexedItem( "tk_etof", m_ngch, m_tk_etof );
193 status = m_tuple4->addIndexedItem( "tp_etof", m_ngch, m_tp_etof );
194 status = m_tuple4->addIndexedItem( "ph_etof", m_ngch, m_ph_etof );
195 status = m_tuple4->addIndexedItem( "rhit_etof", m_ngch, m_rhit_etof );
196 status = m_tuple4->addIndexedItem( "qual_etof", m_ngch, m_qual_etof );
197 status = m_tuple4->addIndexedItem( "ec_toff_e", m_ngch, m_ec_toff_e );
198 status = m_tuple4->addIndexedItem( "ec_toff_mu", m_ngch, m_ec_toff_mu );
199 status = m_tuple4->addIndexedItem( "ec_toff_pi", m_ngch, m_ec_toff_pi );
200 status = m_tuple4->addIndexedItem( "ec_toff_k", m_ngch, m_ec_toff_k );
201 status = m_tuple4->addIndexedItem( "ec_toff_p", m_ngch, m_ec_toff_p );
202 status = m_tuple4->addIndexedItem( "ec_tsig_e", m_ngch, m_ec_tsig_e );
203 status = m_tuple4->addIndexedItem( "ec_tsig_mu", m_ngch, m_ec_tsig_mu );
204 status = m_tuple4->addIndexedItem( "ec_tsig_pi", m_ngch, m_ec_tsig_pi );
205 status = m_tuple4->addIndexedItem( "ec_tsig_k", m_ngch, m_ec_tsig_k );
206 status = m_tuple4->addIndexedItem( "ec_tsig_p", m_ngch, m_ec_tsig_p );
207 status = m_tuple4->addIndexedItem( "ec_tof", m_ngch, m_ec_tof );
208 status = m_tuple4->addIndexedItem( "ptot_btof1", m_ngch, m_ptot_btof1 );
209 status = m_tuple4->addIndexedItem( "cntr_btof1", m_ngch, m_cntr_btof1 );
210 status = m_tuple4->addIndexedItem( "te_btof1", m_ngch, m_te_btof1 );
211 status = m_tuple4->addIndexedItem( "tmu_btof1", m_ngch, m_tmu_btof1 );
212 status = m_tuple4->addIndexedItem( "tpi_btof1", m_ngch, m_tpi_btof1 );
213 status = m_tuple4->addIndexedItem( "tk_btof1", m_ngch, m_tk_btof1 );
214 status = m_tuple4->addIndexedItem( "tp_btof1", m_ngch, m_tp_btof1 );
215 status = m_tuple4->addIndexedItem( "ph_btof1", m_ngch, m_ph_btof1 );
216 status = m_tuple4->addIndexedItem( "zhit_btof1", m_ngch, m_zhit_btof1 );
217 status = m_tuple4->addIndexedItem( "qual_btof1", m_ngch, m_qual_btof1 );
218 status = m_tuple4->addIndexedItem( "b1_toff_e", m_ngch, m_b1_toff_e );
219 status = m_tuple4->addIndexedItem( "b1_toff_mu", m_ngch, m_b1_toff_mu );
220 status = m_tuple4->addIndexedItem( "b1_toff_pi", m_ngch, m_b1_toff_pi );
221 status = m_tuple4->addIndexedItem( "b1_toff_k", m_ngch, m_b1_toff_k );
222 status = m_tuple4->addIndexedItem( "b1_toff_p", m_ngch, m_b1_toff_p );
223 status = m_tuple4->addIndexedItem( "b1_tsig_e", m_ngch, m_b1_tsig_e );
224 status = m_tuple4->addIndexedItem( "b1_tsig_mu", m_ngch, m_b1_tsig_mu );
225 status = m_tuple4->addIndexedItem( "b1_tsig_pi", m_ngch, m_b1_tsig_pi );
226 status = m_tuple4->addIndexedItem( "b1_tsig_k", m_ngch, m_b1_tsig_k );
227 status = m_tuple4->addIndexedItem( "b1_tsig_p", m_ngch, m_b1_tsig_p );
228 status = m_tuple4->addIndexedItem( "b1_tof", m_ngch, m_b1_tof );
229
230 status = m_tuple4->addIndexedItem( "mdedx_pid", m_ngch, m_dedx_pid );
231 status = m_tuple4->addIndexedItem( "mtof1_pid", m_ngch, m_tof1_pid );
232 status = m_tuple4->addIndexedItem( "mtof2_pid", m_ngch, m_tof2_pid );
233 status = m_tuple4->addIndexedItem( "mprob_pid", m_ngch, m_prob_pid );
234 status = m_tuple4->addIndexedItem( "mptrk_pid", m_ngch, m_ptrk_pid );
235 status = m_tuple4->addIndexedItem( "mcost_pid", m_ngch, m_cost_pid );
236 status = m_tuple4->addItem( "mpnp", m_pnp );
237 status = m_tuple4->addItem( "mpnm", m_pnm );
238
239 status =
240 m_tuple4->addIndexedItem( "numHits", m_nggneu, m_numHits ); // Total number of hits
241 status = m_tuple4->addIndexedItem( "secondmoment", m_nggneu, m_secondmoment );
242 status =
243 m_tuple4->addIndexedItem( "mx", m_nggneu, m_x ); // Shower coordinates and errors
244 status = m_tuple4->addIndexedItem( "my", m_nggneu, m_y );
245 status = m_tuple4->addIndexedItem( "mz", m_nggneu, m_z );
246 status = m_tuple4->addIndexedItem( "cosemc", m_nggneu,
247 m_cosemc ); // Shower Counter angles and errors
248 status = m_tuple4->addIndexedItem( "phiemc", m_nggneu, m_phiemc );
249 status = m_tuple4->addIndexedItem( "energy", m_nggneu,
250 m_energy ); // Total energy observed in Emc
251 status = m_tuple4->addIndexedItem( "eseed", m_nggneu, m_eSeed );
252 status = m_tuple4->addIndexedItem( "me9", m_nggneu, m_e3x3 );
253 status = m_tuple4->addIndexedItem( "me25", m_nggneu, m_e5x5 );
254 status = m_tuple4->addIndexedItem( "mlat", m_nggneu, m_lat );
255 status = m_tuple4->addIndexedItem( "ma20", m_nggneu, m_a20 );
256 status = m_tuple4->addIndexedItem( "ma42", m_nggneu, m_a42 );
257 }
258 else
259 {
260 log << MSG::ERROR << " Cannot book N-tuple:" << long( m_tuple4 ) << endmsg;
261 return StatusCode::FAILURE;
262 }
263 }
264
265 if ( service( "THistSvc", m_thsvc ).isFailure() )
266 {
267 log << MSG::ERROR << "Couldn't get THistSvc" << endmsg;
268 return StatusCode::FAILURE;
269 }
270
271 // "DQAHist" is fixed
272 // "Rhopi" is control sample name, just as ntuple case.
273
274 TH1F* mrho0 = new TH1F( "mrho0", "mass for #rho^{0}->#pi^{+}#pi^{-}", 160, 0.0, 3.2 );
275 if ( m_thsvc->regHist( "/DQAHist/Rhopi/mrho0", mrho0 ).isFailure() )
276 { log << MSG::ERROR << "Couldn't register mrho0" << endmsg; }
277
278 TH1F* mrhop = new TH1F( "mrhop", "mass for #rho^{+}->#pi^{+}#pi^{0}", 160, 0.0, 3.2 );
279 if ( m_thsvc->regHist( "/DQAHist/Rhopi/mrhop", mrhop ).isFailure() )
280 { log << MSG::ERROR << "Couldn't register mrhop" << endmsg; }
281
282 TH1F* mrhom = new TH1F( "mrhom", "mass for #rho^{-}->#pi^{-}#pi^{0}", 160, 0.0, 3.2 );
283 if ( m_thsvc->regHist( "/DQAHist/Rhopi/mrhom", mrhom ).isFailure() )
284 { log << MSG::ERROR << "Couldn't register mrhom" << endmsg; }
285
286 TH1F* mpi0 = new TH1F( "mpi0", "mass for #pi^{0}->#gamma #gamma", 50, 0.08, 0.18 );
287 if ( m_thsvc->regHist( "/DQAHist/Rhopi/mpi0", mpi0 ).isFailure() )
288 { log << MSG::ERROR << "Couldn't register mpi0" << endmsg; }
289
290 //
291 //--------end of book--------
292 //
293
294 log << MSG::INFO << "successfully return from initialize()" << endmsg;
295 return StatusCode::SUCCESS;
296}
297
298// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
299StatusCode DQARhopi::execute() {
300
301 // std::cout << "execute()" << std::endl;
302
303 MsgStream log( msgSvc(), name() );
304 log << MSG::INFO << "in execute()" << endmsg;
305
306 setFilterPassed( false );
307
308 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
309 int run = eventHeader->runNumber();
310 int event = eventHeader->eventNumber();
311 log << MSG::DEBUG << "run, evtnum = " << run << " , " << event << endmsg;
312
313 m_run = eventHeader->runNumber();
314 m_rec = eventHeader->eventNumber();
315
316 // cout<<"event "<<event<<endl;
317 Ncut0++;
318 SmartDataPtr<EvtRecEvent> evtRecEvent( eventSvc(), EventModel::EvtRec::EvtRecEvent );
319 log << MSG::INFO << "get event tag OK" << endmsg;
320 log << MSG::DEBUG << "ncharg, nneu, tottks = " << evtRecEvent->totalCharged() << " , "
321 << evtRecEvent->totalNeutral() << " , " << evtRecEvent->totalTracks() << endmsg;
322
323 m_nch = evtRecEvent->totalCharged();
324 m_nneu = evtRecEvent->totalNeutral();
325
326 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol( eventSvc(), EventModel::EvtRec::EvtRecTrackCol );
327
328 //
329 // check x0, y0, z0, r0
330 // suggest cut: |z0|<5 && r0<1
331 //
332 Vint iGood, ipip, ipim, ipnp, ipnm;
333 iGood.clear();
334 ipip.clear();
335 ipim.clear();
336 ipnp.clear();
337 ipnm.clear();
338 Vp4 ppip, ppim, ppnp, ppnm;
339 ppip.clear();
340 ppim.clear();
341 ppnp.clear();
342 ppnm.clear();
343
344 Hep3Vector xorigin( 0, 0, 0 );
345
346 // if (m_reader.isRunNumberValid(runNo)) {
347 IVertexDbSvc* vtxsvc;
348 Gaudi::svcLocator()->service( "VertexDbSvc", vtxsvc ).ignore();
349 if ( vtxsvc->isVertexValid() )
350 {
351 double* dbv = vtxsvc->PrimaryVertex();
352 double* vv = vtxsvc->SigmaPrimaryVertex();
353 xorigin.setX( dbv[0] );
354 xorigin.setY( dbv[1] );
355 xorigin.setZ( dbv[2] );
356 }
357
358 int nCharge = 0;
359 for ( int i = 0; i < evtRecEvent->totalCharged(); i++ )
360 {
361 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
362 if ( !( *itTrk )->isMdcTrackValid() ) continue;
363 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
364 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
365 RecMdcKalTrack* mdcKalTrk = ( *itTrk )->mdcKalTrack();
366
367 double pch = mdcTrk->p();
368 double x0 = mdcTrk->x();
369 double y0 = mdcTrk->y();
370 double z0 = mdcTrk->z();
371 double phi0 = mdcTrk->helix( 1 );
372 double xv = xorigin.x();
373 double yv = xorigin.y();
374 double Rxy = fabs( ( x0 - xv ) * cos( phi0 ) + ( y0 - yv ) * sin( phi0 ) );
375 double m_vx0 = x0;
376 double m_vy0 = y0;
377 double m_vz0 = z0 - xorigin.z();
378 double m_vr0 = Rxy;
379 double m_Vct = cos( mdcTrk->theta() );
380 // m_tuple1->write();
381 if ( fabs( m_vz0 ) >= m_vz0cut ) continue;
382 if ( m_vr0 >= m_vr0cut ) continue;
383 if ( fabs( m_Vct ) >= m_cthcut ) continue;
384
385 iGood.push_back( ( *itTrk )->trackId() );
386 nCharge += mdcTrk->charge();
387 }
388
389 //
390 // Finish Good Charged Track Selection
391 //
392 int nGood = iGood.size();
393 log << MSG::DEBUG << "ngood, totcharge = " << nGood << " , " << nCharge << endmsg;
394 if ( ( nGood != 2 ) || ( nCharge != 0 ) ) { return StatusCode::SUCCESS; }
395 Ncut1++;
396
397 Vint iGam;
398 iGam.clear();
399 for ( int i = evtRecEvent->totalCharged(); i < evtRecEvent->totalTracks(); i++ )
400 {
401 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
402 if ( !( *itTrk )->isEmcShowerValid() ) continue;
403 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
404 Hep3Vector emcpos( emcTrk->x(), emcTrk->y(), emcTrk->z() );
405 // find the nearest charged track
406 double dthe = 200.;
407 double dphi = 200.;
408 double dang = 200.;
409 for ( int j = 0; j < evtRecEvent->totalCharged(); j++ )
410 {
411 EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() + j;
412 if ( !( *jtTrk )->isExtTrackValid() ) continue;
413 RecExtTrack* extTrk = ( *jtTrk )->extTrack();
414 if ( extTrk->emcVolumeNumber() == -1 ) continue;
415 Hep3Vector extpos = extTrk->emcPosition();
416 // double ctht = extpos.cosTheta(emcpos);
417 double angd = extpos.angle( emcpos );
418 double thed = extpos.theta() - emcpos.theta();
419 double phid = extpos.deltaPhi( emcpos );
420 thed = fmod( thed + CLHEP::twopi + CLHEP::twopi + pi, CLHEP::twopi ) - CLHEP::pi;
421 phid = fmod( phid + CLHEP::twopi + CLHEP::twopi + pi, CLHEP::twopi ) - CLHEP::pi;
422
423 if ( fabs( thed ) < fabs( dthe ) ) dthe = thed;
424 if ( fabs( phid ) < fabs( dphi ) ) dphi = phid;
425 if ( angd < dang ) dang = angd;
426 }
427 if ( dang >= 200 ) continue;
428 double eraw = emcTrk->energy();
429 double theta1 = emcTrk->theta();
430 double time = emcTrk->time();
431 dthe = dthe * 180 / ( CLHEP::pi );
432 dphi = dphi * 180 / ( CLHEP::pi );
433 dang = dang * 180 / ( CLHEP::pi );
434 double m_dthe = dthe;
435 double m_dphi = dphi;
436 double m_dang = dang;
437 double m_eraw = eraw;
438 // m_tuple2->write();
439 double fc_theta = fabs( cos( theta1 ) );
440 if ( fc_theta < 0.80 )
441 { // Barrel EMC
442 if ( eraw <= m_energyThreshold / 2 ) continue;
443 }
444 else if ( fc_theta > 0.86 && fc_theta < 0.92 )
445 { // Endcap EMC
446 if ( eraw <= m_energyThreshold ) continue;
447 }
448 else continue;
449
450 if ( time < 0 || time > 14 ) continue;
451 if ( dang < m_gammaAngCut ) continue;
452 //
453 // good photon cut will be set here
454 //
455 iGam.push_back( ( *itTrk )->trackId() );
456 }
457
458 //
459 // Finish Good Photon Selection
460 //
461 int nGam = iGam.size();
462
463 log << MSG::DEBUG << "num Good Photon " << nGam << " , " << evtRecEvent->totalNeutral()
464 << endmsg;
465 if ( nGam < 2 ) { return StatusCode::SUCCESS; }
466 Ncut2++;
467
468 //
469 // Assign 4-momentum to each photon
470 //
471
472 Vp4 pGam;
473 pGam.clear();
474 for ( int i = 0; i < nGam; i++ )
475 {
476 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGam[i];
477 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
478 double eraw = emcTrk->energy();
479 double phi = emcTrk->phi();
480 double the = emcTrk->theta();
481 HepLorentzVector ptrk;
482 ptrk.setPx( eraw * sin( the ) * cos( phi ) );
483 ptrk.setPy( eraw * sin( the ) * sin( phi ) );
484 ptrk.setPz( eraw * cos( the ) );
485 ptrk.setE( eraw );
486
487 // ptrk = ptrk.boost(-0.011,0,0);// boost to cms
488
489 pGam.push_back( ptrk );
490 }
491
492 log << MSG::DEBUG << "liuf1 Good Photon " << endmsg;
493
494 for ( int i = 0; i < nGood; i++ )
495 { // for rhopi without PID
496 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
497 if ( !( *itTrk )->isMdcTrackValid() ) continue;
498 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
499 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
500 RecMdcKalTrack* mdcKalTrk = ( *itTrk )->mdcKalTrack();
501 RecMdcKalTrack::setPidType( RecMdcKalTrack::pion );
502
503 if ( mdcTrk->charge() > 0 )
504 {
505 ipip.push_back( iGood[i] );
506 HepLorentzVector ptrk;
507 ptrk.setPx( mdcKalTrk->px() );
508 ptrk.setPy( mdcKalTrk->py() );
509 ptrk.setPz( mdcKalTrk->pz() );
510 double p3 = ptrk.mag();
511 ptrk.setE( sqrt( p3 * p3 + mpi * mpi ) );
512 ppip.push_back( ptrk );
513 }
514 else
515 {
516 ipim.push_back( iGood[i] );
517 HepLorentzVector ptrk;
518 ptrk.setPx( mdcKalTrk->px() );
519 ptrk.setPy( mdcKalTrk->py() );
520 ptrk.setPz( mdcKalTrk->pz() );
521 double p3 = ptrk.mag();
522 ptrk.setE( sqrt( p3 * p3 + mpi * mpi ) );
523 ppim.push_back( ptrk );
524 }
525 } // without PID
526
527 int npip = ipip.size();
528 int npim = ipim.size();
529 log << MSG::DEBUG << "num of pion " << ipip.size() << "," << ipim.size() << endmsg;
530 Ncut3++;
531
532 //***********************************************//
533 // the angle between the two charged tracks //
534 //***********************************************//
535
536 int langcc = 0;
537 double dthec = 200.;
538 double dphic = 200.;
539 double dangc = 200.;
540 for ( int i = 0; i < npip; i++ )
541 {
542 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + ipip[i];
543 RecMdcTrack* mdcTrk1 = ( *itTrk )->mdcTrack();
544 RecMdcKalTrack* mdcKalTrk1 = ( *itTrk )->mdcKalTrack();
545 Hep3Vector emcpos( ppip[i].px(), ppip[i].py(), ppip[i].pz() );
546
547 for ( int j = 0; j < npim; j++ )
548 {
549 EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() + ipim[j];
550 RecMdcTrack* mdcTrk2 = ( *jtTrk )->mdcTrack();
551 RecMdcKalTrack* mdcKalTrk2 = ( *jtTrk )->mdcKalTrack();
552
553 HepLorentzVector pippim = ppip[i] + ppim[j];
554
555 Hep3Vector extpos( ppim[j].px(), ppim[j].py(), ppim[j].pz() );
556
557 double angd = extpos.angle( emcpos );
558 double thed = extpos.theta() - emcpos.theta();
559 double phid = extpos.deltaPhi( emcpos );
560 thed = fmod( thed + CLHEP::twopi + CLHEP::twopi + pi, CLHEP::twopi ) - CLHEP::pi;
561 phid = fmod( phid + CLHEP::twopi + CLHEP::twopi + pi, CLHEP::twopi ) - CLHEP::pi;
562
563 m_dthec[langcc] = thed * 180 / ( CLHEP::pi );
564 m_dphic[langcc] = phid * 180 / ( CLHEP::pi );
565 m_dangc[langcc] = angd * 180 / ( CLHEP::pi );
566 m_mspippim[langcc] = pippim.m();
567 langcc++;
568 }
569 }
570 m_nangecc = langcc;
571
572 //
573 // Loop each gamma pair, check ppi0 and pTot
574 //
575
576 double m_m2gg, m_momentpi0;
577 HepLorentzVector pTot, p2g;
578
579 log << MSG::DEBUG << "liuf2 Good Photon " << langcc << endmsg;
580 //******************************************************//
581 // asign the momentum of MDC and KALFIT //
582 // the deposite energy of EMC for charged tracks //
583 //******************************************************//
584 double m_momentpip, m_momentpim, extmot1, extmot2;
585 double emcTg1 = 0.0;
586 double emcTg2 = 0.0;
587 double nlaypi1 = 0;
588 double nhit1 = 0;
589 double nlaypi2 = 0;
590 double nhit2 = 0;
591
592 EvtRecTrackIterator itTrk11 = evtRecTrkCol->begin() + ipip[0];
593 RecMdcTrack* mdcTrk11 = ( *itTrk11 )->mdcTrack();
594 RecMdcKalTrack* mdcKalTrk11 = ( *itTrk11 )->mdcKalTrack();
595 RecEmcShower* emcTrk11 = ( *itTrk11 )->emcShower();
596 RecMucTrack* mucTrk11 = ( *itTrk11 )->mucTrack();
597 double phi01 = mdcTrk11->helix( 1 );
598
599 EvtRecTrackIterator itTrk12 = evtRecTrkCol->begin() + ipim[0];
600 RecMdcTrack* mdcTrk12 = ( *itTrk12 )->mdcTrack();
601 RecMdcKalTrack* mdcKalTrk12 = ( *itTrk12 )->mdcKalTrack();
602 RecEmcShower* emcTrk12 = ( *itTrk12 )->emcShower();
603 RecMucTrack* mucTrk12 = ( *itTrk12 )->mucTrack();
604 double phi02 = mdcTrk12->helix( 1 );
605
606 m_vxpin = mdcTrk11->x();
607 m_vypin = mdcTrk11->y();
608 m_vzpin = mdcTrk11->z() - xorigin.z();
609 m_vrpin = fabs( ( mdcTrk11->x() - xorigin.x() ) * cos( phi01 ) +
610 ( mdcTrk11->y() - xorigin.y() ) * sin( phi01 ) );
611 m_costhepin = cos( mdcTrk11->theta() );
612
613 m_momentpip = mdcTrk11->p();
614 m_ppx = mdcTrk11->px();
615 m_ppy = mdcTrk11->py();
616 m_ppz = mdcTrk11->pz();
617
618 m_vxp = mdcKalTrk11->x();
619 m_vyp = mdcKalTrk11->y();
620 m_vzp = mdcKalTrk11->z() - xorigin.z();
621 m_vrp = fabs( ( mdcKalTrk11->x() - xorigin.x() ) * cos( phi01 ) +
622 ( mdcKalTrk11->y() - xorigin.y() ) * sin( phi01 ) );
623 m_costhep = cos( mdcKalTrk11->theta() );
624
625 // extmot1=mdcKalTrk11->p(); // only has value when read data from dst file
626 m_ppxkal = mdcKalTrk11->px();
627 m_ppykal = mdcKalTrk11->py();
628 m_ppzkal = mdcKalTrk11->pz();
629 extmot1 = sqrt( m_ppxkal * m_ppxkal + m_ppykal * m_ppykal + m_ppzkal * m_ppzkal );
630
631 m_vxmin = mdcTrk12->x();
632 m_vymin = mdcTrk12->y();
633 m_vzmin = mdcTrk12->z();
634 m_vrmin = fabs( ( mdcTrk12->x() - xorigin.x() ) * cos( phi02 ) +
635 ( mdcTrk12->y() - xorigin.y() ) * sin( phi02 ) );
636 m_costhemin = cos( mdcTrk12->theta() );
637
638 m_momentpim = mdcTrk12->p();
639 m_mpx = mdcTrk12->px();
640 m_mpy = mdcTrk12->py();
641 m_mpz = mdcTrk12->pz();
642
643 m_vxm = mdcKalTrk12->x();
644 m_vym = mdcKalTrk12->y();
645 m_vzm = mdcKalTrk12->z();
646 m_vrm = fabs( ( mdcKalTrk12->x() - xorigin.x() ) * cos( phi02 ) +
647 ( mdcKalTrk12->y() - xorigin.y() ) * sin( phi02 ) );
648 m_costhem = cos( mdcKalTrk12->theta() );
649
650 // extmot2 =mdcKalTrk12->p();
651 m_mpxkal = mdcKalTrk12->px();
652 m_mpykal = mdcKalTrk12->py();
653 m_mpzkal = mdcKalTrk12->pz();
654 extmot2 = sqrt( m_mpxkal * m_mpxkal + m_mpykal * m_mpykal + m_mpzkal * m_mpzkal );
655
656 if ( ( *itTrk11 )->isEmcShowerValid() ) { emcTg1 = emcTrk11->energy(); }
657 if ( ( *itTrk12 )->isEmcShowerValid() ) { emcTg2 = emcTrk12->energy(); }
658 if ( ( *itTrk11 )->isMucTrackValid() )
659 {
660 nlaypi1 = mucTrk11->numLayers();
661 nhit1 = mucTrk11->numHits();
662 }
663 if ( ( *itTrk12 )->isMucTrackValid() )
664 {
665 nlaypi2 = mucTrk12->numLayers();
666 nhit2 = mucTrk12->numHits();
667 }
668
669 m_laypi1 = nlaypi1;
670 m_hit1 = nhit1;
671 m_laypi2 = nlaypi2;
672 m_hit2 = nhit2;
673
674 log << MSG::DEBUG << "liuf3 Good Photon " << ipim[0] << endmsg;
675
676 RecMdcKalTrack* pipTrk = ( *( evtRecTrkCol->begin() + ipip[0] ) )->mdcKalTrack();
677 RecMdcKalTrack* pimTrk = ( *( evtRecTrkCol->begin() + ipim[0] ) )->mdcKalTrack();
678
679 log << MSG::DEBUG << "liuf4 Good Photon " << endmsg;
680
681 WTrackParameter wvpipTrk, wvpimTrk, wvkpTrk, wvkmTrk;
682 wvpipTrk = WTrackParameter( mpi, pipTrk->getZHelix(), pipTrk->getZError() );
683 wvpimTrk = WTrackParameter( mpi, pimTrk->getZHelix(), pimTrk->getZError() );
684
685 wvkpTrk = WTrackParameter( mk, pipTrk->getZHelixK(), pipTrk->getZErrorK() ); // kaon
686 wvkmTrk = WTrackParameter( mk, pimTrk->getZHelixK(), pimTrk->getZErrorK() ); // kaon
687
688 /* Default is pion, for other particles:
689 wvppTrk = WTrackParameter(mp, pipTrk->getZHelixP(), pipTrk->getZErrorP());//proton
690 wvmupTrk = WTrackParameter(mmu, pipTrk->getZHelixMu(), pipTrk->getZErrorMu());//muon
691 wvepTrk = WTrackParameter(me, pipTrk->getZHelixE(), pipTrk->getZErrorE());//electron
692 wvkpTrk = WTrackParameter(mk, pipTrk->getZHelixK(), pipTrk->getZErrorK());//kaon
693 */
694 //
695 // Test vertex fit
696 //
697
698 HepPoint3D vx( 0., 0., 0. );
699 HepSymMatrix Evx( 3, 0 );
700 double bx = 1E+6;
701 double by = 1E+6;
702 double bz = 1E+6;
703 Evx[0][0] = bx * bx;
704 Evx[1][1] = by * by;
705 Evx[2][2] = bz * bz;
706
707 VertexParameter vxpar;
708 vxpar.setVx( vx );
709 vxpar.setEvx( Evx );
710
711 //****************************************************//
712 // Test vertex fit //
713 //****************************************************//
714
715 VertexFit* vtxfit = VertexFit::instance();
716
717 //****************************************************//
718 // if the charged particle is kaon //
719 //****************************************************//
720 double chi5 = 9999.0;
721 double jkpi0 = -0.5;
722 double jkpkm = 0.0;
723 double jkpp0 = 0.0;
724 double jkmp0 = 0.0;
725 vtxfit->init();
726 vtxfit->AddTrack( 0, wvkpTrk );
727 vtxfit->AddTrack( 1, wvkmTrk );
728 vtxfit->AddVertex( 0, vxpar, 0, 1 );
729 if ( vtxfit->Fit( 0 ) )
730 {
731 vtxfit->Swim( 0 );
732 WTrackParameter wkp = vtxfit->wtrk( 0 );
733 WTrackParameter wkm = vtxfit->wtrk( 1 );
734
735 KinematicFit* kmfit = KinematicFit::instance();
736
737 //
738 // Apply Kinematic 4C fit
739 //
740
741 double chisq = 9999.;
742 for ( int i = 0; i < nGam - 1; i++ )
743 {
744 RecEmcShower* g1Trk = ( *( evtRecTrkCol->begin() + iGam[i] ) )->emcShower();
745 for ( int j = i + 1; j < nGam; j++ )
746 {
747 RecEmcShower* g2Trk = ( *( evtRecTrkCol->begin() + iGam[j] ) )->emcShower();
748 kmfit->init();
749 // kmfit->setDynamicerror(1);
750 kmfit->AddTrack( 0, wkp );
751 kmfit->AddTrack( 1, wkm );
752 kmfit->AddTrack( 2, 0.0, g1Trk );
753 kmfit->AddTrack( 3, 0.0, g2Trk );
754 kmfit->AddFourMomentum( 0, ecms );
755 bool oksq = kmfit->Fit();
756
757 if ( oksq )
758 {
759 double chi2 = kmfit->chisq();
760 if ( chi2 < chi5 )
761 {
762 HepLorentzVector kpi0 = kmfit->pfit( 2 ) + kmfit->pfit( 3 );
763 HepLorentzVector kpkm = kmfit->pfit( 0 ) + kmfit->pfit( 1 );
764 HepLorentzVector kpp0 = kmfit->pfit( 0 ) + kpi0;
765 HepLorentzVector kmp0 = kmfit->pfit( 1 ) + kpi0;
766 chi5 = kmfit->chisq();
767 jkpi0 = kpi0.m();
768 jkpkm = kpkm.m();
769 jkpp0 = kpp0.m();
770 jkmp0 = kmp0.m();
771 }
772 }
773 }
774 }
775 } // vetex is true
776
777 //****************************************************//
778 // if the charged particles are pions for real //
779 //****************************************************//
780
781 vtxfit->init();
782 vtxfit->AddTrack( 0, wvpipTrk );
783 vtxfit->AddTrack( 1, wvpimTrk );
784 vtxfit->AddVertex( 0, vxpar, 0, 1 );
785 if ( !vtxfit->Fit( 0 ) ) return StatusCode::SUCCESS;
786 vtxfit->Swim( 0 );
787
788 WTrackParameter wpip = vtxfit->wtrk( 0 );
789 WTrackParameter wpim = vtxfit->wtrk( 1 );
790
791 log << MSG::DEBUG << "liuf5 Good Photon " << endmsg;
792
793 KinematicFit* kmfit = KinematicFit::instance();
794
795 //
796 // Apply Kinematic 4C fit
797 //
798
799 double chisq = 9999.;
800 int ig1 = -1;
801 int ig2 = -1;
802 HepLorentzVector gg1, gg2;
803 for ( int i = 0; i < nGam - 1; i++ )
804 {
805 RecEmcShower* g1Trk = ( *( evtRecTrkCol->begin() + iGam[i] ) )->emcShower();
806 for ( int j = i + 1; j < nGam; j++ )
807 {
808 RecEmcShower* g2Trk = ( *( evtRecTrkCol->begin() + iGam[j] ) )->emcShower();
809 kmfit->init();
810 // kmfit->setDynamicerror(1);
811 kmfit->AddTrack( 0, wpip );
812 kmfit->AddTrack( 1, wpim );
813 kmfit->AddTrack( 2, 0.0, g1Trk );
814 kmfit->AddTrack( 3, 0.0, g2Trk );
815 kmfit->AddFourMomentum( 0, ecms );
816 bool oksq = kmfit->Fit();
817 if ( oksq )
818 {
819 double chi2 = kmfit->chisq();
820 if ( chi2 < chisq )
821 {
822 chisq = chi2;
823 ig1 = iGam[i];
824 ig2 = iGam[j];
825 gg1 = pGam[i];
826 gg2 = pGam[j];
827 }
828 }
829 }
830 }
831
832 p2g = gg1 + gg2;
833 m_pmax = gg1.e() > gg2.e() ? gg1.e() : gg2.e();
834 m_m2gg = p2g.m();
835 m_cosang = ( gg1.e() - gg2.e() ) / p2g.rho();
836 m_momentpi0 = sqrt( p2g.px() * p2g.px() + p2g.py() * p2g.py() + p2g.pz() * p2g.pz() );
837 log << MSG::DEBUG << " chisq for 4c " << chisq << endmsg;
838 if ( chisq > 200 ) { return StatusCode::SUCCESS; }
839
840 // select charge track and nneu track
841 Vint jGood;
842 jGood.clear();
843 jGood.push_back( ipip[0] );
844 jGood.push_back( ipim[0] );
845 m_ngch = jGood.size();
846
847 Vint jGgam;
848 jGgam.clear();
849 jGgam.push_back( ig1 );
850 jGgam.push_back( ig2 );
851 m_nggneu = jGgam.size();
852
853 HepLorentzVector ppip1 = ppip[0];
854 HepLorentzVector ppim1 = ppim[0];
855
856 HepLorentzVector Ppipboost = ppip1.boost( u_cms );
857 HepLorentzVector Ppimboost = ppim1.boost( u_cms );
858 Hep3Vector p3pi1 = Ppipboost.vect(); // pi1
859 Hep3Vector p3pi2 = Ppimboost.vect(); // pi2
860 m_anglepm = ( p3pi1.angle( p3pi2 ) ) * 180 / ( CLHEP::pi );
861
862 //*******************************************************//
863
864 RecEmcShower* g1Trk = ( *( evtRecTrkCol->begin() + ig1 ) )->emcShower();
865 RecEmcShower* g2Trk = ( *( evtRecTrkCol->begin() + ig2 ) )->emcShower();
866 kmfit->init();
867 // kmfit->setDynamicerror(1);
868 kmfit->AddTrack( 0, wpip );
869 kmfit->AddTrack( 1, wpim );
870 kmfit->AddTrack( 2, 0.0, g1Trk );
871 kmfit->AddTrack( 3, 0.0, g2Trk );
872 kmfit->AddFourMomentum( 0, ecms );
873 bool oksq = kmfit->Fit();
874 if ( !oksq ) return StatusCode::SUCCESS;
875
876 HepLorentzVector ppi0 = kmfit->pfit( 2 ) + kmfit->pfit( 3 );
877 HepLorentzVector prho0 = kmfit->pfit( 0 ) + kmfit->pfit( 1 );
878 HepLorentzVector prhop = kmfit->pfit( 0 ) + ppi0;
879 HepLorentzVector prhom = kmfit->pfit( 1 ) + ppi0;
880 HepLorentzVector pgam2pi1 = prho0 + kmfit->pfit( 2 );
881 HepLorentzVector pgam2pi2 = prho0 + kmfit->pfit( 3 );
882 HepLorentzVector pgam11 = kmfit->pfit( 2 );
883 HepLorentzVector pgam12 = kmfit->pfit( 3 );
884
885 /*
886 HepLorentzVector ppi0_a=ppi0;
887 HepLorentzVector pgam11_a =pgam11;
888 HepLorentzVector pgam12_a =pgam12;
889
890 Hep3Vector pi0_cms = -ppi0_a.boostVector();
891 HepLorentzVector pgam1 = pgam11_a.boost(pi0_cms);
892 HepLorentzVector pgam2 = pgam12_a.boost(pi0_cms);
893 */
894 m_chi1 = kmfit->chisq();
895 m_mpi0 = ppi0.m();
896 m_prho0 = prho0.m();
897 m_prhop = prhop.m();
898 m_prhom = prhom.m();
899 m_good = nGood;
900 m_gam = nGam;
901
902 m_pip = npip;
903 m_pim = npim;
904 m_2gam = m_m2gg;
905 m_outpi0 = m_momentpi0;
906 m_outpip = m_momentpip;
907 m_outpim = m_momentpim;
908 m_enpip = emcTg1;
909 m_dcpip = extmot1;
910 m_enpim = emcTg2;
911 m_dcpim = extmot2;
912 m_pipf = kmfit->pfit( 0 ).rho();
913 m_pimf = kmfit->pfit( 1 ).rho();
914 m_pi0f = ppi0.rho();
915
916 m_chi5 = chi5;
917 m_kpi0 = jkpi0;
918 m_kpkm = jkpkm;
919 m_kpp0 = jkpp0;
920 m_kmp0 = jkmp0;
921 m_pgam2pi1 = pgam2pi1.m();
922 m_pgam2pi2 = pgam2pi2.m();
923 cosva1 = kmfit->pfit( 2 ).rho();
924 cosva2 = kmfit->pfit( 3 ).rho();
925 // m_pe1 =pe1;
926 // m_pe2 =pe2;
927
928 //
929 // fill information of dedx and tof
930 //
931
932 //
933 // check dedx infomation
934 //
935
936 for ( int ii = 0; ii < m_ngch; ii++ )
937 {
938 // dedx
939 m_ptrk[ii] = 9999.0;
940 m_chie[ii] = 9999.0;
941 m_chimu[ii] = 9999.0;
942 m_chipi[ii] = 9999.0;
943 m_chik[ii] = 9999.0;
944 m_chip[ii] = 9999.0;
945 m_ghit[ii] = 9999.0;
946 m_thit[ii] = 9999.0;
947 m_probPH[ii] = 9999.0;
948 m_normPH[ii] = 9999.0;
949
950 // endtof
951 m_cntr_etof[ii] = 9999.0;
952 m_ptot_etof[ii] = 9999.0;
953 m_ph_etof[ii] = 9999.0;
954 m_rhit_etof[ii] = 9999.0;
955 m_qual_etof[ii] = 9999.0;
956 m_te_etof[ii] = 9999.0;
957 m_tmu_etof[ii] = 9999.0;
958 m_tpi_etof[ii] = 9999.0;
959 m_tk_etof[ii] = 9999.0;
960 m_tp_etof[ii] = 9999.0;
961 m_ec_tof[ii] = 9999.0;
962 m_ec_toff_e[ii] = 9999.0;
963 m_ec_toff_mu[ii] = 9999.0;
964 m_ec_toff_pi[ii] = 9999.0;
965 m_ec_toff_k[ii] = 9999.0;
966 m_ec_toff_p[ii] = 9999.0;
967 m_ec_tsig_e[ii] = 9999.0;
968 m_ec_tsig_mu[ii] = 9999.0;
969 m_ec_tsig_pi[ii] = 9999.0;
970 m_ec_tsig_k[ii] = 9999.0;
971 m_ec_tsig_p[ii] = 9999.0;
972
973 // barrel tof
974 m_cntr_btof1[ii] = 9999.0;
975 m_ptot_btof1[ii] = 9999.0;
976 m_ph_btof1[ii] = 9999.0;
977 m_zhit_btof1[ii] = 9999.0;
978 m_qual_btof1[ii] = 9999.0;
979 m_te_btof1[ii] = 9999.0;
980 m_tmu_btof1[ii] = 9999.0;
981 m_tpi_btof1[ii] = 9999.0;
982 m_tk_btof1[ii] = 9999.0;
983 m_tp_btof1[ii] = 9999.0;
984 m_b1_tof[ii] = 9999.0;
985 m_b1_toff_e[ii] = 9999.0;
986 m_b1_toff_mu[ii] = 9999.0;
987 m_b1_toff_pi[ii] = 9999.0;
988 m_b1_toff_k[ii] = 9999.0;
989 m_b1_toff_p[ii] = 9999.0;
990 m_b1_tsig_e[ii] = 9999.0;
991 m_b1_tsig_mu[ii] = 9999.0;
992 m_b1_tsig_pi[ii] = 9999.0;
993 m_b1_tsig_k[ii] = 9999.0;
994 m_b1_tsig_p[ii] = 9999.0;
995 // pid
996 m_dedx_pid[ii] = 9999.0;
997 m_tof1_pid[ii] = 9999.0;
998 m_tof2_pid[ii] = 9999.0;
999 m_prob_pid[ii] = 9999.0;
1000 m_ptrk_pid[ii] = 9999.0;
1001 m_cost_pid[ii] = 9999.0;
1002 }
1003
1004 int indx0 = 0;
1005 for ( int i = 0; i < m_ngch; i++ )
1006 {
1007 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + jGood[i];
1008 if ( !( *itTrk )->isMdcTrackValid() ) continue;
1009 if ( !( *itTrk )->isMdcDedxValid() ) continue;
1010 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
1011 RecMdcDedx* dedxTrk = ( *itTrk )->mdcDedx();
1012 m_ptrk[indx0] = mdcTrk->p();
1013
1014 m_chie[indx0] = dedxTrk->chiE();
1015 m_chimu[indx0] = dedxTrk->chiMu();
1016 m_chipi[indx0] = dedxTrk->chiPi();
1017 m_chik[indx0] = dedxTrk->chiK();
1018 m_chip[indx0] = dedxTrk->chiP();
1019 m_ghit[indx0] = dedxTrk->numGoodHits();
1020 m_thit[indx0] = dedxTrk->numTotalHits();
1021 m_probPH[indx0] = dedxTrk->probPH();
1022 m_normPH[indx0] = dedxTrk->normPH();
1023 indx0++;
1024 }
1025
1026 //
1027 // check TOF infomation
1028 //
1029
1030 int indx1 = 0;
1031 for ( int i = 0; i < m_ngch; i++ )
1032 {
1033 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + jGood[i];
1034 if ( !( *itTrk )->isMdcTrackValid() ) continue;
1035 if ( !( *itTrk )->isTofTrackValid() ) continue;
1036
1037 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
1038 SmartRefVector<RecTofTrack> tofTrkCol = ( *itTrk )->tofTrack();
1039
1040 double ptrk = mdcTrk->p();
1041 SmartRefVector<RecTofTrack>::iterator iter_tof = tofTrkCol.begin();
1042 for ( ; iter_tof != tofTrkCol.end(); iter_tof++ )
1043 {
1044 TofHitStatus* status = new TofHitStatus;
1045 status->setStatus( ( *iter_tof )->status() );
1046 if ( !( status->is_barrel() ) )
1047 { // endcap
1048 if ( !( status->is_counter() ) ) continue; // ?
1049 if ( status->layer() != 1 ) continue; // layer1
1050 double path = ( *iter_tof )->path(); // ?
1051 double tof = ( *iter_tof )->tof();
1052 double ph = ( *iter_tof )->ph();
1053 double rhit = ( *iter_tof )->zrhit();
1054 double qual = 0.0 + ( *iter_tof )->quality();
1055 double cntr = 0.0 + ( *iter_tof )->tofID();
1056 double texp[5];
1057 double tsig[5];
1058 for ( int j = 0; j < 5; j++ )
1059 { // 0 e, 1 mu, 2 pi, 3 K, 4 p
1060 texp[j] = ( *iter_tof )->texp( j );
1061 // tsig[j] = (*iter_tof)->sigma(j);
1062 // toffset[j] = (*iter_tof)->offset(j);
1063 }
1064 m_cntr_etof[indx1] = cntr;
1065 m_ptot_etof[indx1] = ptrk;
1066 m_ph_etof[indx1] = ph;
1067 m_rhit_etof[indx1] = rhit;
1068 m_qual_etof[indx1] = qual;
1069 m_te_etof[indx1] = tof - texp[0];
1070 m_tmu_etof[indx1] = tof - texp[1];
1071 m_tpi_etof[indx1] = tof - texp[2];
1072 m_tk_etof[indx1] = tof - texp[3];
1073 m_tp_etof[indx1] = tof - texp[4];
1074
1075 m_ec_tof[indx1] = tof;
1076
1077 m_ec_toff_e[indx1] = ( *iter_tof )->toffset( 0 );
1078 m_ec_toff_mu[indx1] = ( *iter_tof )->toffset( 1 );
1079 m_ec_toff_pi[indx1] = ( *iter_tof )->toffset( 2 );
1080 m_ec_toff_k[indx1] = ( *iter_tof )->toffset( 3 );
1081 m_ec_toff_p[indx1] = ( *iter_tof )->toffset( 4 );
1082
1083 m_ec_tsig_e[indx1] = ( *iter_tof )->sigma( 0 );
1084 m_ec_tsig_mu[indx1] = ( *iter_tof )->sigma( 1 );
1085 m_ec_tsig_pi[indx1] = ( *iter_tof )->sigma( 2 );
1086 m_ec_tsig_k[indx1] = ( *iter_tof )->sigma( 3 );
1087 m_ec_tsig_p[indx1] = ( *iter_tof )->sigma( 4 );
1088 }
1089 else
1090 { // barrel
1091 if ( !( status->is_cluster() ) ) continue; // ?
1092 double path = ( *iter_tof )->path(); // ?
1093 double tof = ( *iter_tof )->tof();
1094 double ph = ( *iter_tof )->ph();
1095 double rhit = ( *iter_tof )->zrhit();
1096 double qual = 0.0 + ( *iter_tof )->quality();
1097 double cntr = 0.0 + ( *iter_tof )->tofID();
1098 double texp[5];
1099 for ( int j = 0; j < 5; j++ ) { texp[j] = ( *iter_tof )->texp( j ); }
1100 m_cntr_btof1[indx1] = cntr;
1101 m_ptot_btof1[indx1] = ptrk;
1102 m_ph_btof1[indx1] = ph;
1103 m_zhit_btof1[indx1] = rhit;
1104 m_qual_btof1[indx1] = qual;
1105 m_te_btof1[indx1] = tof - texp[0];
1106 m_tmu_btof1[indx1] = tof - texp[1];
1107 m_tpi_btof1[indx1] = tof - texp[2];
1108 m_tk_btof1[indx1] = tof - texp[3];
1109 m_tp_btof1[indx1] = tof - texp[4];
1110
1111 m_b1_tof[indx1] = tof;
1112
1113 m_b1_toff_e[indx1] = ( *iter_tof )->toffset( 0 );
1114 m_b1_toff_mu[indx1] = ( *iter_tof )->toffset( 1 );
1115 m_b1_toff_pi[indx1] = ( *iter_tof )->toffset( 2 );
1116 m_b1_toff_k[indx1] = ( *iter_tof )->toffset( 3 );
1117 m_b1_toff_p[indx1] = ( *iter_tof )->toffset( 4 );
1118
1119 m_b1_tsig_e[indx1] = ( *iter_tof )->sigma( 0 );
1120 m_b1_tsig_mu[indx1] = ( *iter_tof )->sigma( 1 );
1121 m_b1_tsig_pi[indx1] = ( *iter_tof )->sigma( 2 );
1122 m_b1_tsig_k[indx1] = ( *iter_tof )->sigma( 3 );
1123 m_b1_tsig_p[indx1] = ( *iter_tof )->sigma( 4 );
1124 }
1125 delete status;
1126 }
1127 indx1++;
1128 } // loop all charged track
1129
1130 //
1131 // Assign 4-momentum to each charged track
1132 //
1133 int indx2 = 0;
1134 ParticleID* pid = ParticleID::instance();
1135 for ( int i = 0; i < m_ngch; i++ )
1136 {
1137 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + jGood[i];
1138 // if(pid) delete pid;
1139 pid->init();
1140 pid->setMethod( pid->methodProbability() );
1141 // pid->setMethod(pid->methodLikelihood()); //for Likelihood Method
1142
1143 pid->setChiMinCut( 4 );
1144 pid->setRecTrack( *itTrk );
1145 pid->usePidSys( pid->useDedx() | pid->useTof1() | pid->useTof2() ); // use PID sub-system
1146 pid->identify( pid->onlyPion() | pid->onlyKaon() ); // seperater Pion/Kaon
1147 // pid->identify(pid->onlyPion());
1148 // pid->identify(pid->onlyKaon());
1149 pid->calculate();
1150 if ( !( pid->IsPidInfoValid() ) ) continue;
1151 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
1152
1153 m_dedx_pid[indx2] = pid->chiDedx( 2 );
1154 m_tof1_pid[indx2] = pid->chiTof1( 2 );
1155 m_tof2_pid[indx2] = pid->chiTof2( 2 );
1156 m_prob_pid[indx2] = pid->probPion();
1157
1158 // if(pid->probPion() < 0.001 || (pid->probPion() < pid->probKaon())) continue;
1159 // if(pid->probPion() < 0.001) continue;
1160 // if(pid->pdf(2)<pid->pdf(3)) continue; // for Likelihood Method(0=electron 1=muon
1161 // 2=pion 3=kaon 4=proton)
1162
1163 RecMdcKalTrack* mdcKalTrk =
1164 ( *itTrk )->mdcKalTrack(); // After ParticleID, use RecMdcKalTrack substitute
1165 // RecMdcTrack
1166 RecMdcKalTrack::setPidType( RecMdcKalTrack::pion ); // PID can set to electron, muon, pion,
1167 // kaon and proton;The default setting
1168 // is pion
1169
1170 // m_ptrk_pid[indx2] = mdcKalTrk->p();
1171 m_cost_pid[indx2] = cos( mdcTrk->theta() );
1172
1173 HepLorentzVector ptrk;
1174 ptrk.setPx( mdcKalTrk->px() );
1175 ptrk.setPy( mdcKalTrk->py() );
1176 ptrk.setPz( mdcKalTrk->pz() );
1177 double p3 = ptrk.mag();
1178 ptrk.setE( sqrt( p3 * p3 + mpi * mpi ) );
1179 // ptrk = ptrk.boost(-0.011,0,0);//boost to cms
1180 m_ptrk_pid[indx2] = p3;
1181
1182 if ( mdcTrk->charge() > 0 && pid->probPion() > pid->probKaon() )
1183 {
1184 ipnp.push_back( jGood[i] );
1185 ppnp.push_back( ptrk );
1186 } // plus charge with with PID
1187 if ( mdcTrk->charge() < 0 && pid->probPion() > pid->probKaon() )
1188 {
1189 ipnm.push_back( jGood[i] );
1190 ppnm.push_back( ptrk );
1191 } // minus charge with with PID
1192 }
1193 int npnp = ipnp.size();
1194 int npnm = ipnm.size();
1195
1196 m_pnp = npnp;
1197 m_pnm = npnm;
1198
1199 int iphoton = 0;
1200 for ( int i = 0; i < m_nggneu; i++ )
1201 {
1202 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + jGgam[i];
1203 if ( !( *itTrk )->isEmcShowerValid() ) continue;
1204 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
1205 m_numHits[iphoton] = emcTrk->numHits();
1206 m_secondmoment[iphoton] = emcTrk->secondMoment();
1207 m_x[iphoton] = emcTrk->x();
1208 m_y[iphoton] = emcTrk->y();
1209 m_z[iphoton] = emcTrk->z();
1210 m_cosemc[iphoton] = cos( emcTrk->theta() );
1211 m_phiemc[iphoton] = emcTrk->phi();
1212 m_energy[iphoton] = emcTrk->energy();
1213 m_eSeed[iphoton] = emcTrk->eSeed();
1214 m_e3x3[iphoton] = emcTrk->e3x3();
1215 m_e5x5[iphoton] = emcTrk->e5x5();
1216 m_lat[iphoton] = emcTrk->latMoment();
1217 m_a20[iphoton] = emcTrk->a20Moment();
1218 m_a42[iphoton] = emcTrk->a42Moment();
1219 iphoton++;
1220 }
1221 m_tuple4->write().ignore();
1222 Ncut4++;
1223
1224 if ( kmfit->chisq() >= chi5 ) return StatusCode::SUCCESS;
1225 if ( pgam2pi2.m() <= 1.0 ) return StatusCode::SUCCESS;
1226 if ( pgam2pi1.m() <= 1.0 ) return StatusCode::SUCCESS;
1227 if ( nGam != 2 ) return StatusCode::SUCCESS;
1228 if ( emcTg1 / extmot1 >= 0.8 ) return StatusCode::SUCCESS;
1229 if ( emcTg2 / extmot2 >= 0.8 ) return StatusCode::SUCCESS;
1230 Ncut5++;
1231
1232 // DQA
1233 TH1* h( 0 );
1234 if ( m_thsvc->getHist( "/DQAHist/Rhopi/mpi0", h ).isSuccess() ) { h->Fill( ppi0.m() ); }
1235 else { log << MSG::ERROR << "Couldn't retrieve mpi0" << endmsg; }
1236
1237 if ( fabs( ppi0.m() - 0.135 ) >= 0.015 ) return StatusCode::SUCCESS;
1238 Ncut6++;
1239
1240 if ( m_thsvc->getHist( "/DQAHist/Rhopi/mrho0", h ).isSuccess() ) { h->Fill( prho0.m() ); }
1241 else { log << MSG::ERROR << "Couldn't retrieve mrho0" << endmsg; }
1242
1243 if ( m_thsvc->getHist( "/DQAHist/Rhopi/mrhop", h ).isSuccess() ) { h->Fill( prhop.m() ); }
1244 else { log << MSG::ERROR << "Couldn't retrieve mrhop" << endmsg; }
1245
1246 if ( m_thsvc->getHist( "/DQAHist/Rhopi/mrhom", h ).isSuccess() ) { h->Fill( prhom.m() ); }
1247 else { log << MSG::ERROR << "Couldn't retrieve mrhom" << endmsg; }
1248
1249 // Pid: 1 - electron, 2 - muon, 3 - pion, 4 - kaon, 5 - proton
1250
1251 ( *( evtRecTrkCol->begin() + ipip[0] ) )->setPartId( 3 );
1252 ( *( evtRecTrkCol->begin() + ipim[0] ) )->setPartId( 3 );
1253
1254 // Quality: defined by whether dE/dx or TOF is used to identify particle
1255 // 0 - no dE/dx, no TOF (can be used for dE/dx and TOF calibration)
1256 // 1 - only dE/dx (can be used for TOF calibration)
1257 // 2 - only TOF (can be used for dE/dx calibration)
1258 // 3 - Both dE/dx and TOF
1259
1260 ( *( evtRecTrkCol->begin() + ipip[0] ) )->setQuality( 0 );
1261 ( *( evtRecTrkCol->begin() + ipim[0] ) )->setQuality( 0 );
1262
1263 setFilterPassed( true );
1264
1265 return StatusCode::SUCCESS;
1266}
1267
1268// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
1269StatusCode DQARhopi::finalize() {
1270 cout << "total number: " << Ncut0 << endl;
1271 cout << "nGood==2, nCharge==0: " << Ncut1 << endl;
1272 cout << "nGam>=2: " << Ncut2 << endl;
1273 cout << "Pass Pid: " << Ncut3 << endl;
1274 cout << "Pass 4C: " << Ncut4 << endl;
1275 cout << "final cut without pi0:" << Ncut5 << endl;
1276 cout << "final cut with pi0: " << Ncut6 << endl;
1277 MsgStream log( msgSvc(), name() );
1278 log << MSG::INFO << "in finalize()" << endmsg;
1279 return StatusCode::SUCCESS;
1280}
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition CosmicGenerator.h:36
u_cms
const Hep3Vector u_cms
Definition DQADtagAlg.cxx:56
time
Double_t time
Definition DataBase/tau_mode.c:7
EvtRecTrackIterator
EvtRecTrackCol::iterator EvtRecTrackIterator
Definition Event/EvtRecEvent/include/EvtRecEvent/EvtRecTrack.h:106
mpi
double mpi
Definition EvtPhokharaDef.hh:17
pi
double pi
Definition EvtPhokharaDef.hh:17
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
mk
const double mk
Definition Gam4pikp.cxx:33
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
theta1
@ theta1
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrkBase/TrkKalDeriv.h:24
ntupleSvc
INTupleSvc * ntupleSvc()
Definition ServiceAccessor.h:51
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:12
DQARhopi::initialize
StatusCode initialize()
Definition DQARhopi.cxx:70
DQARhopi::execute
StatusCode execute()
Definition DQARhopi.cxx:299
DQARhopi::DQARhopi
DQARhopi(const std::string &name, ISvcLocator *pSvcLocator)
Definition DQARhopi.cxx:54
DQARhopi::finalize
StatusCode finalize()
Definition DQARhopi.cxx:1269
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
KinematicFit::instance
static KinematicFit * instance()
Definition KinematicFit.cxx:18
KinematicFit::AddFourMomentum
void AddFourMomentum(int number, HepLorentzVector p4)
Definition KinematicFit.cxx:341
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
TrackPool::AddTrack
void AddTrack(const int number, const double mass, const RecMdcTrack *trk)
Definition TrackPool.cxx:21
VertexFit::wtrk
WTrackParameter wtrk(int n) const
Definition Analysis/VertexFit/include/VertexFit/VertexFit.h:91
VertexFit::init
void init()
Definition VertexFit.cxx:27
VertexFit::AddVertex
void AddVertex(int number, VertexParameter vpar, std::vector< int > lis)
Definition VertexFit.cxx:85
VertexFit::instance
static VertexFit * instance()
Definition VertexFit.cxx:15
VertexFit::Fit
bool Fit()
Definition VertexFit.cxx:288
ecms
const double ecms
Definition inclkstar.cxx:26