BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
DQASelHadron.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 "EmcRecEventModel/RecEmcEventModel.h"
16#include "EmcRecGeoSvc/EmcRecBarrelGeo.h"
17#include "EmcRecGeoSvc/IEmcRecGeoSvc.h"
18#include "EvTimeEvent/RecEsTime.h"
19#include "EventModel/Event.h"
20#include "EventModel/EventHeader.h"
21#include "EventModel/EventModel.h"
22#include "EvtRecEvent/EvtRecEvent.h"
23#include "EvtRecEvent/EvtRecTrack.h"
24#include "McTruth/McParticle.h"
25#include "ParticleID/ParticleID.h"
26#include "VertexDbSvc/IVertexDbSvc.h"
27#include "VertexFit/FastVertexFit.h"
28#include "VertexFit/Helix.h"
29#include "VertexFit/KinematicFit.h"
30#include "VertexFit/SecondVertexFit.h"
31#include "VertexFit/VertexFit.h"
32
33using CLHEP::Hep2Vector;
34using CLHEP::Hep3Vector;
35using CLHEP::HepLorentzVector;
36#include "CLHEP/Geometry/Point3D.h"
37#ifndef ENABLE_BACKWARDS_COMPATIBILITY
38typedef HepGeom::Point3D<double> HepPoint3D;
39#endif
40
41#include "DQASelHadron.h"
42
43#ifndef ENABLE_BACKWARDS_COMPATIBILITY
44typedef HepGeom::Point3D<double> HepPoint3D;
45#endif
46using CLHEP::HepLorentzVector;
47
48const double mpi = 0.13957;
49const double mk = 0.493677;
50const double xmass[5] = { 0.000511, 0.105658, 0.139570, 0.493677, 0.938272 };
51const double velc = 299.792458; // tof path unit in mm
52typedef std::vector<int> Vint;
53typedef std::vector<HepLorentzVector> Vp4;
54// declare one counter
55static int counter[10] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
56static int nhadron = 0;
57static int n0prong = 0;
58static int n2prong = 0;
59static int n4prong = 0;
60static int ng4prong = 0;
61
62/////////////////////////////////////////////////////////////////////////////
63DECLARE_COMPONENT( DQASelHadron )
64DQASelHadron::DQASelHadron( const std::string& name, ISvcLocator* pSvcLocator )
65 : Algorithm( name, pSvcLocator ) {
66
67 // Declare the properties
68 declareProperty( "writentuple", m_writentuple = false );
69 declareProperty( "useVertexDB", m_useVertexDB = false );
70 declareProperty( "ecms", m_ecms = 3.097 );
71 declareProperty( "beamangle", m_beamangle = 0.022 );
72 declareProperty( "Vr0cut", m_vr0cut = 1.0 );
73 declareProperty( "Vz0cut", m_vz0cut = 10.0 );
74 declareProperty( "Coscut", m_coscut = 0.93 );
75
76 declareProperty( "EnergyThreshold", m_energyThreshold = 0.04 );
77 declareProperty( "GammaPhiCut", m_gammaPhiCut = 20.0 );
78 declareProperty( "GammaThetaCut", m_gammaThetaCut = 20.0 );
79 declareProperty( "GammaTrkCut", m_gammaTrkCut = 20.0 );
80 declareProperty( "GammaTLCut", m_gammatlCut = 0 );
81 declareProperty( "GammaTHCut", m_gammathCut = 60 );
82
83 declareProperty( "acoll_h_cut", m_acoll_h_cut = 10. );
84 declareProperty( "poeb_h_cut", m_poeb_h_cut = 0.2 );
85 declareProperty( "dtof_h_cut", m_dtof_h_cut = 6. );
86 declareProperty( "eop_h_cut", m_eop_h_cut = 0.2 );
87 declareProperty( "etotal_h_cut", m_etotal_h_cut = 0.2 );
88 declareProperty( "ngam_h_cut", m_ngam_h_cut = 2 );
89 declareProperty( "br_h_cut", m_br_h_cut = 0.65 );
90 declareProperty( "bz_h_cut", m_bz_h_cut = 0.7 );
91 declareProperty( "thr_h_cut", m_thr_h_cut = 0.5 );
92
93 // normally, MDC+EMC, otherwise EMC only
94 declareProperty( "m_useEMConly", m_useEMConly = false );
95 declareProperty( "m_usePID", m_usePID = false ); // sub-system is under study
96 declareProperty( "m_useMDC", m_useMDC = true );
97 declareProperty( "m_useDEDX", m_useDEDX = false ); // not used
98 declareProperty( "m_useTOF", m_useTOF = false ); // sub-system is under study
99 declareProperty( "m_useEMC", m_useEMC = true );
100 declareProperty( "m_useMUC", m_useMUC = false ); // efficiency
101}
102
103// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
104StatusCode DQASelHadron::initialize() {
105 MsgStream log( msgSvc(), name() );
106
107 log << MSG::INFO << "in initialize()" << endmsg;
108 StatusCode status;
109 status = service( "THistSvc", m_thistsvc );
110 if ( status.isFailure() )
111 {
112 log << MSG::INFO << "Unable to retrieve pointer to THistSvc" << endmsg;
113 return status;
114 }
115
116 m_ha_costheta = new TH1F( "PHY_HAD_SUM_costheta", "PHY_HAD_SUM_costheta", 100, -1, 1 );
117 status = m_thistsvc->regHist( "/DQAHist/Hadron/ha_costheta", m_ha_costheta );
118 m_ha_phi = new TH1F( "PHY_HAD_SUM_phi", "PHY_HAD_SUM_phi", 128, -3.2, 3.2 );
119 status = m_thistsvc->regHist( "/DQAHist/Hadron/ha_phi", m_ha_phi );
120 m_ha_pmax = new TH1F( "PHY_HAD_SUM_pmax", "PHY_HAD_SUM_pmax", 100, 0, 2 );
121 status = m_thistsvc->regHist( "/DQAHist/Hadron/ha_pmax", m_ha_pmax );
122 m_ha_emax = new TH1F( "PHY_HAD_SUM_emax", "PHY_HAD_SUM_emax", 100, 0, 2 );
123 status = m_thistsvc->regHist( "/DQAHist/Hadron/ha_emax", m_ha_emax );
124 m_ha_etot = new TH1F( "PHY_HAD_SUM_etot", "PHY_HAD_SUM_etot", 100, 0, 4 );
125 status = m_thistsvc->regHist( "/DQAHist/Hadron/ha_etot", m_ha_etot );
126 m_ha_br = new TH1F( "PHY_HAD_SUM_br", "PHY_HAD_SUM_br", 100, 0, 2 );
127 status = m_thistsvc->regHist( "/DQAHist/Hadron/ha_br", m_ha_br );
128 m_ha_bz = new TH1F( "PHY_HAD_SUM_bz", "PHY_HAD_SUM_bz", 100, 0, 2 );
129 status = m_thistsvc->regHist( "/DQAHist/Hadron/ha_bz", m_ha_bz );
130 m_ha_nneu = new TH1I( "PHY_HAD_SUM_nneu", "PHY_HAD_SUM_nneu", 20, 0, 20 );
131 status = m_thistsvc->regHist( "/DQAHist/Hadron/ha_nneu", m_ha_nneu );
132 m_ha_nchg = new TH1I( "PHY_HAD_SUM_nchg", "PHY_HAD_SUM_nchg", 20, 0, 20 );
133 status = m_thistsvc->regHist( "/DQAHist/Hadron/ha_nchg", m_ha_nchg );
134
135 m_ha_vx = new TH1F( "PHY_HAD_FLS_vx", "PHY_HAD_FLS_vx", 100, -1., 1. );
136 status = m_thistsvc->regHist( "/DQAHist/Hadron/ha_vx", m_ha_vx );
137 m_ha_vy = new TH1F( "PHY_HAD_FLS_vy", "PHY_HAD_FLS_vy", 100, -1., 1. );
138 status = m_thistsvc->regHist( "/DQAHist/Hadron/ha_vy", m_ha_vy );
139 m_ha_vz = new TH1F( "PHY_HAD_FLS_vz", "PHY_HAD_FLS_vz", 100, -10.0, 10. );
140 status = m_thistsvc->regHist( "/DQAHist/Hadron/ha_vz", m_ha_vz );
141
142 NTuplePtr nt1( ntupleSvc(), "DQAFILE/Hadron" );
143 if ( nt1 ) m_tuple1 = nt1;
144 else
145 {
146 m_tuple1 = ntupleSvc()->book( "DQAFILE/Hadron", CLID_ColumnWiseTuple, "N-Tuple" );
147 if ( m_tuple1 )
148 {
149 status = m_tuple1->addItem( "run", m_run );
150 status = m_tuple1->addItem( "rec", m_rec );
151 status = m_tuple1->addItem( "Nchrg", m_ncharg );
152 status = m_tuple1->addItem( "Nneu", m_nneu, 0, 40 );
153 status = m_tuple1->addItem( "NGch", m_ngch, 0, 40 );
154 status = m_tuple1->addItem( "NGam", m_nGam );
155
156 status = m_tuple1->addItem( "hadrontag", m_hadrontag );
157
158 status = m_tuple1->addItem( "br", m_br );
159 status = m_tuple1->addItem( "bz", m_bz );
160 status = m_tuple1->addItem( "evis", m_evis );
161 status = m_tuple1->addItem( "thr", m_thr );
162
163 status = m_tuple1->addItem( "acoll", m_acoll );
164 status = m_tuple1->addItem( "acopl", m_acopl );
165 status = m_tuple1->addItem( "deltatof", m_deltatof );
166 status = m_tuple1->addItem( "eop1", m_eop1 );
167 status = m_tuple1->addItem( "eop2", m_eop2 );
168 status = m_tuple1->addItem( "eoeb1", m_eoeb1 );
169 status = m_tuple1->addItem( "eoeb2", m_eoeb2 );
170 status = m_tuple1->addItem( "poeb1", m_poeb1 );
171 status = m_tuple1->addItem( "poeb2", m_poeb2 );
172 status = m_tuple1->addItem( "etoeb1", m_etoeb1 );
173 status = m_tuple1->addItem( "etoeb2", m_etoeb2 );
174 status = m_tuple1->addItem( "mucinfo1", m_mucinfo1 );
175 status = m_tuple1->addItem( "mucinfo2", m_mucinfo2 );
176
177 status = m_tuple1->addIndexedItem( "delang", m_nneu, m_delang );
178 status = m_tuple1->addIndexedItem( "delphi", m_nneu, m_delphi );
179 status = m_tuple1->addIndexedItem( "delthe", m_nneu, m_delthe );
180 status = m_tuple1->addIndexedItem( "npart", m_nneu, m_npart );
181 status = m_tuple1->addIndexedItem( "nemchits", m_nneu, m_nemchits );
182 status = m_tuple1->addIndexedItem( "module", m_nneu, m_module );
183 status = m_tuple1->addIndexedItem( "x", m_nneu, m_x );
184 status = m_tuple1->addIndexedItem( "y", m_nneu, m_y );
185 status = m_tuple1->addIndexedItem( "z", m_nneu, m_z );
186 // status = m_tuple1->addIndexedItem ("px",m_nneu, m_px);
187 // status = m_tuple1->addIndexedItem ("py",m_nneu, m_py);
188 // status = m_tuple1->addIndexedItem ("pz",m_nneu, m_pz);
189 status = m_tuple1->addIndexedItem( "theta", m_nneu, m_theta );
190 status = m_tuple1->addIndexedItem( "phi", m_nneu, m_phi );
191 status = m_tuple1->addIndexedItem( "dx", m_nneu, m_dx );
192 status = m_tuple1->addIndexedItem( "dy", m_nneu, m_dy );
193 status = m_tuple1->addIndexedItem( "dz", m_nneu, m_dz );
194 status = m_tuple1->addIndexedItem( "dtheta", m_nneu, m_dtheta );
195 status = m_tuple1->addIndexedItem( "dphi", m_nneu, m_dphi );
196 status = m_tuple1->addIndexedItem( "energy", m_nneu, m_energy );
197 status = m_tuple1->addIndexedItem( "dE", m_nneu, m_dE );
198 status = m_tuple1->addIndexedItem( "eSeed", m_nneu, m_eSeed );
199 status = m_tuple1->addIndexedItem( "nSeed", m_nneu, m_nSeed );
200 status = m_tuple1->addIndexedItem( "e3x3", m_nneu, m_e3x3 );
201 status = m_tuple1->addIndexedItem( "e5x5", m_nneu, m_e5x5 );
202 status = m_tuple1->addIndexedItem( "secondMoment", m_nneu, m_secondMoment );
203 status = m_tuple1->addIndexedItem( "latMoment", m_nneu, m_latMoment );
204 status = m_tuple1->addIndexedItem( "a20Moment", m_nneu, m_a20Moment );
205 status = m_tuple1->addIndexedItem( "a42Moment", m_nneu, m_a42Moment );
206 status = m_tuple1->addIndexedItem( "getTime", m_nneu, m_getTime );
207 status = m_tuple1->addIndexedItem( "getEAll", m_nneu, m_getEAll );
208
209 status = m_tuple1->addIndexedItem( "charge", m_ngch, m_charge );
210 status = m_tuple1->addIndexedItem( "vx", m_ngch, m_vx0 );
211 status = m_tuple1->addIndexedItem( "vy", m_ngch, m_vy0 );
212 status = m_tuple1->addIndexedItem( "vz", m_ngch, m_vz0 );
213
214 status = m_tuple1->addIndexedItem( "px", m_ngch, m_px );
215 status = m_tuple1->addIndexedItem( "py", m_ngch, m_py );
216 status = m_tuple1->addIndexedItem( "pz", m_ngch, m_pz );
217 status = m_tuple1->addIndexedItem( "p", m_ngch, m_p );
218
219 status = m_tuple1->addIndexedItem( "kal_vx", m_ngch, m_kal_vx0 );
220 status = m_tuple1->addIndexedItem( "kal_vy", m_ngch, m_kal_vy0 );
221 status = m_tuple1->addIndexedItem( "kal_vz", m_ngch, m_kal_vz0 );
222
223 status = m_tuple1->addIndexedItem( "kal_px", m_ngch, m_kal_px );
224 status = m_tuple1->addIndexedItem( "kal_py", m_ngch, m_kal_py );
225 status = m_tuple1->addIndexedItem( "kal_pz", m_ngch, m_kal_pz );
226 status = m_tuple1->addIndexedItem( "kal_p", m_ngch, m_kal_p );
227
228 status = m_tuple1->addIndexedItem( "probPH", m_ngch, m_probPH );
229 status = m_tuple1->addIndexedItem( "normPH", m_ngch, m_normPH );
230 status = m_tuple1->addIndexedItem( "chie", m_ngch, m_chie );
231 status = m_tuple1->addIndexedItem( "chimu", m_ngch, m_chimu );
232 status = m_tuple1->addIndexedItem( "chipi", m_ngch, m_chipi );
233 status = m_tuple1->addIndexedItem( "chik", m_ngch, m_chik );
234 status = m_tuple1->addIndexedItem( "chip", m_ngch, m_chip );
235 status = m_tuple1->addIndexedItem( "ghit", m_ngch, m_ghit );
236 status = m_tuple1->addIndexedItem( "thit", m_ngch, m_thit );
237
238 status = m_tuple1->addIndexedItem( "e_emc", m_ngch, m_e_emc );
239 status = m_tuple1->addIndexedItem( "phi_emc", m_ngch, m_phi_emc );
240 status = m_tuple1->addIndexedItem( "theta_emc", m_ngch, m_theta_emc );
241
242 status = m_tuple1->addIndexedItem( "nhit_muc", m_ngch, m_nhit_muc );
243 status = m_tuple1->addIndexedItem( "nlay_muc", m_ngch, m_nlay_muc );
244 status = m_tuple1->addIndexedItem( "t_btof", m_ngch, m_t_btof );
245 status = m_tuple1->addIndexedItem( "t_etof", m_ngch, m_t_etof );
246 status = m_tuple1->addIndexedItem( "qual_etof", m_ngch, m_qual_etof );
247 status = m_tuple1->addIndexedItem( "tof_etof", m_ngch, m_tof_etof );
248 status = m_tuple1->addIndexedItem( "te_etof", m_ngch, m_te_etof );
249 status = m_tuple1->addIndexedItem( "tmu_etof", m_ngch, m_tmu_etof );
250 status = m_tuple1->addIndexedItem( "tpi_etof", m_ngch, m_tpi_etof );
251 status = m_tuple1->addIndexedItem( "tk_etof", m_ngch, m_tk_etof );
252 status = m_tuple1->addIndexedItem( "tp_etof", m_ngch, m_tp_etof );
253
254 status = m_tuple1->addIndexedItem( "qual_btof1", m_ngch, m_qual_btof1 );
255 status = m_tuple1->addIndexedItem( "tof_btof1", m_ngch, m_tof_btof1 );
256 status = m_tuple1->addIndexedItem( "te_btof1", m_ngch, m_te_btof1 );
257 status = m_tuple1->addIndexedItem( "tmu_btof1", m_ngch, m_tmu_btof1 );
258 status = m_tuple1->addIndexedItem( "tpi_btof1", m_ngch, m_tpi_btof1 );
259 status = m_tuple1->addIndexedItem( "tk_btof1", m_ngch, m_tk_btof1 );
260 status = m_tuple1->addIndexedItem( "tp_btof1", m_ngch, m_tp_btof1 );
261
262 status = m_tuple1->addIndexedItem( "qual_btof2", m_ngch, m_qual_btof2 );
263 status = m_tuple1->addIndexedItem( "tof_btof2", m_ngch, m_tof_btof2 );
264 status = m_tuple1->addIndexedItem( "te_btof2", m_ngch, m_te_btof2 );
265 status = m_tuple1->addIndexedItem( "tmu_btof2", m_ngch, m_tmu_btof2 );
266 status = m_tuple1->addIndexedItem( "tpi_btof2", m_ngch, m_tpi_btof2 );
267 status = m_tuple1->addIndexedItem( "tk_btof2", m_ngch, m_tk_btof2 );
268 status = m_tuple1->addIndexedItem( "tp_btof2", m_ngch, m_tp_btof2 );
269 status = m_tuple1->addIndexedItem( "pidcode", m_ngch, m_pidcode );
270 status = m_tuple1->addIndexedItem( "pidprob", m_ngch, m_pidprob );
271 status = m_tuple1->addIndexedItem( "pidchiDedx", m_ngch, m_pidchiDedx );
272 status = m_tuple1->addIndexedItem( "pidchiTof1", m_ngch, m_pidchiTof1 );
273 status = m_tuple1->addIndexedItem( "pidchiTof2", m_ngch, m_pidchiTof2 );
274
275 status = m_tuple1->addItem( "px_cms_ep", m_px_cms_ep ); // momentum of electron+
276 status = m_tuple1->addItem( "py_cms_ep", m_py_cms_ep ); // momentum of electron+
277 status = m_tuple1->addItem( "pz_cms_ep", m_pz_cms_ep ); // momentum of electron+
278 status = m_tuple1->addItem( "e_cms_ep", m_e_cms_ep ); // momentum of electron+
279 status = m_tuple1->addItem( "cos_ep", m_cos_ep ); // momentum of electron+
280 status = m_tuple1->addItem( "px_cms_em", m_px_cms_em ); // momentum of electron-
281 status = m_tuple1->addItem( "py_cms_em", m_py_cms_em ); // momentum of electron-
282 status = m_tuple1->addItem( "pz_cms_em", m_pz_cms_em ); // momentum of electron-
283 status = m_tuple1->addItem( "e_cms_em", m_e_cms_em ); // momentum of electron-
284 status = m_tuple1->addItem( "cos_em", m_cos_em ); // momentum of electron-
285 status = m_tuple1->addItem( "mass_ee", m_mass_ee ); //
286 status = m_tuple1->addItem( "emax", m_emax ); //
287 status = m_tuple1->addItem( "esum", m_esum ); //
288 status = m_tuple1->addItem( "npip", m_npip );
289 status = m_tuple1->addItem( "npim", m_npim );
290 status = m_tuple1->addItem( "nkp", m_nkp );
291 status = m_tuple1->addItem( "nkm", m_nkm );
292 status = m_tuple1->addItem( "np", m_np );
293 status = m_tuple1->addItem( "npb", m_npb );
294
295 status = m_tuple1->addItem( "nep", m_nep );
296 status = m_tuple1->addItem( "nem", m_nem );
297 status = m_tuple1->addItem( "nmup", m_nmup );
298 status = m_tuple1->addItem( "nmum", m_nmum );
299 }
300 else
301 {
302 log << MSG::ERROR << " Cannot book N-tuple:" << long( m_tuple1 ) << endmsg;
303 return StatusCode::FAILURE;
304 }
305 }
306
307 //
308 //--------end of book--------
309 //
310
311 log << MSG::INFO << "successfully return from initialize()" << endmsg;
312 return StatusCode::SUCCESS;
313}
314
315// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
316StatusCode DQASelHadron::execute() {
317 setFilterPassed( false );
318 const double beamEnergy = m_ecms / 2.;
319 const HepLorentzVector p_cms( m_ecms * sin( m_beamangle * 0.5 ), 0.0, 0.0, m_ecms );
320 const Hep3Vector u_cms = -p_cms.boostVector();
321 MsgStream log( msgSvc(), name() );
322 log << MSG::INFO << "in execute()" << endmsg;
323
324 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
325 if ( !eventHeader )
326 {
327 log << MSG::FATAL << "Could not find Event Header" << endmsg;
328 return StatusCode::SUCCESS;
329 }
330
331 m_run = eventHeader->runNumber();
332 m_rec = eventHeader->eventNumber();
333
334 SmartDataPtr<EvtRecEvent> evtRecEvent( eventSvc(), EventModel::EvtRec::EvtRecEvent );
335 if ( !evtRecEvent )
336 {
337 log << MSG::FATAL << "Could not find EvtRecEvent" << endmsg;
338 return StatusCode::SUCCESS;
339 }
340 log << MSG::INFO << "ncharg, nneu, tottks = " << evtRecEvent->totalCharged() << " , "
341 << evtRecEvent->totalNeutral() << " , " << evtRecEvent->totalTracks() << endmsg;
342 // if(evtRecEvent->totalNeutral()>30)return sc;
343 m_ncharg = evtRecEvent->totalCharged();
344
345 m_nneu = evtRecEvent->totalNeutral();
346
347 HepPoint3D vx( 0., 0., 0. );
348 HepSymMatrix Evx( 3, 0 );
349 if ( m_useVertexDB )
350 {
351 IVertexDbSvc* vtxsvc;
352 Gaudi::svcLocator()->service( "VertexDbSvc", vtxsvc ).ignore();
353 if ( vtxsvc->isVertexValid() )
354 {
355 double* dbv = vtxsvc->PrimaryVertex();
356 double* vv = vtxsvc->SigmaPrimaryVertex();
357 // if (m_reader.isRunNumberValid( m_run)) {
358 // HepVector dbv = m_reader.PrimaryVertex( m_run);
359 // HepVector vv = m_reader.SigmaPrimaryVertex( m_run);
360 vx.setX( dbv[0] );
361 vx.setY( dbv[1] );
362 vx.setZ( dbv[2] );
363 Evx[0][0] = vv[0] * vv[0];
364 Evx[0][1] = vv[0] * vv[1];
365 Evx[1][1] = vv[1] * vv[1];
366 Evx[1][2] = vv[1] * vv[2];
367 Evx[2][2] = vv[2] * vv[2];
368 }
369 }
370
371 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol( eventSvc(), EventModel::EvtRec::EvtRecTrackCol );
372 if ( !evtRecTrkCol )
373 {
374 log << MSG::FATAL << "Could not find EvtRecTrackCol" << endmsg;
375 return StatusCode::SUCCESS;
376 }
377 Vint iGood;
378 iGood.clear();
379
380 int nCharge = 0;
381
382 for ( int i = 0; i < evtRecEvent->totalCharged(); i++ )
383 {
384 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
385 if ( !( *itTrk )->isMdcTrackValid() ) continue;
386 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
387
388 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
389 double pch = mdcTrk->p();
390 double x0 = mdcTrk->x();
391 double y0 = mdcTrk->y();
392 double z0 = mdcTrk->z();
393 // double phi0=mdcTrk->helix(1);
394 // double xv=vx.x();
395 // double yv=vx.y();
396 // double zv=vx.z();
397 // double Rxy=(x0-xv)*cos(phi0)+(y0-yv)*sin(phi0);
398 // double m_vx0 = x0;
399 // double m_vy0 = y0;
400 // double m_vz0 = z0;
401 // double m_vr0 = Rxy;
402 // if(fabs(z0) >= m_vz0cut) continue;
403 // if(fabs(Rxy) >= m_vr0cut) continue;
404
405 // if(fabs(m_vz0) >= m_vz0cut) continue;
406 // if(m_vr0 >= m_vr0cut) continue;
407
408 HepVector a = mdcTrk->helix();
409 HepSymMatrix Ea = mdcTrk->err();
410 HepPoint3D point0( 0., 0., 0. );
411 HepPoint3D IP( vx[0], vx[1], vx[2] );
412 VFHelix helixip( point0, a, Ea );
413 helixip.pivot( IP );
414 HepVector vecipa = helixip.a();
415 double Rvxy0 = fabs( vecipa[0] ); // the distance to IP in xy plane
416 double Rvz0 = vecipa[3]; // the distance to IP in z direction
417 double Rvphi0 = vecipa[1];
418 if ( fabs( Rvz0 ) >= m_vz0cut ) continue;
419 if ( fabs( Rvxy0 ) >= m_vr0cut ) continue;
420
421 // double cost = cos(mdcTrk->theta());
422 // if(fabs(cost) >= m_coscut ) continue;
423 // iGood.push_back((*itTrk)->trackId());
424 iGood.push_back( i );
425 nCharge += mdcTrk->charge();
426 }
427
428 //
429 // Finish Good Charged Track Selection
430 //
431 int nGood = iGood.size();
432 m_ngch = nGood;
433 log << MSG::DEBUG << "ngood, totcharge = " << nGood << " , " << nCharge << endmsg;
434
435 if ( m_ngch < 2 || m_ngch > 20 || ( nCharge > 4 ) ) { return StatusCode::SUCCESS; }
436 counter[1]++;
437
438 //
439 // Particle ID
440 //
441 Vint ipip, ipim, iep, iem, imup, imum;
442 ipip.clear();
443 ipim.clear();
444 iep.clear();
445 iem.clear();
446 imup.clear();
447 imum.clear();
448
449 if ( m_usePID )
450 {
451 ParticleID* pid = ParticleID::instance();
452 for ( int i = 0; i < m_ngch; i++ )
453 {
454 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
455 // if(pid) delete pid;
456 pid->init();
457 pid->setMethod( pid->methodProbability() );
458 pid->setChiMinCut( 4 );
459 pid->setRecTrack( *itTrk );
460 pid->usePidSys( pid->useDedx() | pid->useTof1() |
461 pid->useTof2() ); //|pid->useEmc()|pid->useMuc()); // use PID sub-system
462 pid->identify( pid->onlyElectron() | pid->onlyMuon() |
463 pid->onlyPion() ); // seperater Pion/Kaon/Proton
464 pid->calculate();
465 if ( !( pid->IsPidInfoValid() ) ) continue;
466 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
467 /// RecMdcKalTrack* mdcKalTrk = 0 ;
468 /// if((*itTrk)->isMdcKalTrackValid()) mdcKalTrk = (*itTrk)->mdcKalTrack();
469 double prob_pi = pid->probPion();
470 double prob_K = pid->probKaon();
471 double prob_p = pid->probProton();
472 double prob_e = pid->probElectron();
473 double prob_mu = pid->probMuon();
474 // std::cout << "prob "<< prob_pi << ", "<< prob_K << ", "<< prob_p << std::endl;
475 HepLorentzVector ptrk;
476 ptrk.setPx( mdcTrk->px() );
477 ptrk.setPy( mdcTrk->py() );
478 ptrk.setPz( mdcTrk->pz() );
479 double p3 = ptrk.mag();
480
481 m_pidcode[i] = 1;
482 m_pidprob[i] = pid->prob( 1 );
483 m_pidchiDedx[i] = pid->chiDedx( 1 );
484 m_pidchiTof1[i] = pid->chiTof1( 1 );
485 m_pidchiTof2[i] = pid->chiTof2( 1 );
486 if ( mdcTrk->charge() > 0 ) { imup.push_back( iGood[i] ); }
487 if ( mdcTrk->charge() < 0 ) { imum.push_back( iGood[i] ); }
488 }
489 }
490 m_nep = iep.size();
491 m_nem = iem.size();
492 m_nmup = imup.size();
493 m_nmum = imum.size();
494
495 counter[2]++;
496
497 //
498 // Good neutral track selection
499 //
500 Vint iGam;
501 iGam.clear();
502 int iphoton = 0;
503 for ( int i = evtRecEvent->totalCharged(); i < evtRecEvent->totalTracks(); i++ )
504 {
505 if ( i >= evtRecTrkCol->size() ) break;
506 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
507 if ( !( *itTrk )->isEmcShowerValid() ) continue;
508 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
509 Hep3Vector emcpos( emcTrk->x(), emcTrk->y(), emcTrk->z() );
510
511 RecEmcID showerId = emcTrk->getShowerId();
512 unsigned int npart = EmcID::barrel_ec( showerId );
513 int n = emcTrk->numHits();
514 int module = emcTrk->module();
515 double x = emcTrk->x();
516 double y = emcTrk->y();
517 double z = emcTrk->z();
518 double dx = emcTrk->dx();
519 double dy = emcTrk->dy();
520 double dth = emcTrk->dtheta();
521 double dph = emcTrk->dphi();
522 double dz = emcTrk->dz();
523 double energy = emcTrk->energy();
524 double dE = emcTrk->dE();
525 double eSeed = emcTrk->eSeed();
526 double e3x3 = emcTrk->e3x3();
527 double e5x5 = emcTrk->e5x5();
528 double secondMoment = emcTrk->secondMoment();
529 double latMoment = emcTrk->latMoment();
530 double getTime = emcTrk->time();
531 double getEAll = emcTrk->getEAll();
532 double a20Moment = emcTrk->a20Moment();
533 double a42Moment = emcTrk->a42Moment();
534 // int phigap=emcTrk->PhiGap();
535 // int thetagap=emcTrk->ThetaGap();
536 // double getETof2x1 = emcTrk->getETof2x1();
537 // double getETof2x3 = emcTrk->getETof2x3();
538 // double getELepton = emcTrk->getELepton();
539 double nseed = 0; //(emcTrk->getCluster() )->getSeedSize() ;
540 HepPoint3D EmcPos( x, y, z );
541 m_nemchits[iphoton] = n;
542 m_npart[iphoton] = npart;
543 m_module[iphoton] = module;
544 m_theta[iphoton] = EmcPos.theta();
545 m_phi[iphoton] = EmcPos.phi();
546 m_x[iphoton] = x;
547 m_y[iphoton] = y;
548 m_z[iphoton] = z;
549 m_dx[iphoton] = dx;
550 m_dy[iphoton] = dy;
551 m_dz[iphoton] = dz;
552 m_dtheta[iphoton] = dth;
553 m_dphi[iphoton] = dph;
554 m_energy[iphoton] = energy;
555 m_dE[iphoton] = dE;
556 m_eSeed[iphoton] = eSeed;
557 m_nSeed[iphoton] = nseed;
558 m_e3x3[iphoton] = e3x3;
559 m_e5x5[iphoton] = e5x5;
560 m_secondMoment[iphoton] = secondMoment;
561 m_latMoment[iphoton] = latMoment;
562 m_getTime[iphoton] = getTime;
563 m_getEAll[iphoton] = getEAll;
564 m_a20Moment[iphoton] = a20Moment;
565 m_a42Moment[iphoton] = a42Moment;
566
567 // m_getELepton[iphoton]=getELepton;
568 // m_getETof2x1[iphoton]=getETof2x1;
569 // m_getETof2x3[iphoton]=getETof2x3;
570 // m_PhiGap[iphoton]=phigap;
571 // m_ThetaGap[iphoton]=thetagap;
572 double dthe = 200.;
573 double dphi = 200.;
574 double dang = 200.;
575
576 // find the nearest charged track
577 for ( int j = 0; j < nGood; j++ )
578 {
579
580 EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() + iGood[j];
581 if ( !( *jtTrk )->isMdcTrackValid() ) continue;
582 RecMdcTrack* jtmdcTrk = ( *jtTrk )->mdcTrack();
583 double jtcharge = jtmdcTrk->charge();
584 if ( !( *jtTrk )->isExtTrackValid() ) continue;
585 RecExtTrack* extTrk = ( *jtTrk )->extTrack();
586 if ( extTrk->emcVolumeNumber() == -1 ) continue;
587 Hep3Vector extpos = extTrk->emcPosition();
588 // double ctht = extpos.cosTheta(emcpos);
589 double angd = extpos.angle( emcpos );
590 double thed = extpos.theta() - emcpos.theta();
591 double phid = extpos.deltaPhi( emcpos );
592 thed = fmod( thed + CLHEP::twopi + CLHEP::twopi + pi, CLHEP::twopi ) - CLHEP::pi;
593 phid = fmod( phid + CLHEP::twopi + CLHEP::twopi + pi, CLHEP::twopi ) - CLHEP::pi;
594
595 if ( fabs( thed ) < fabs( dthe ) ) dthe = thed;
596 if ( fabs( phid ) < fabs( dphi ) ) dphi = phid;
597 if ( angd < dang ) dang = angd;
598 }
599
600 //
601 // good photon cut will be set here
602 //
603
604 dthe = dthe * 180 / ( CLHEP::pi );
605 dphi = dphi * 180 / ( CLHEP::pi );
606 dang = dang * 180 / ( CLHEP::pi );
607 double eraw = emcTrk->energy();
608 double phi = emcTrk->phi();
609 double the = emcTrk->theta();
610
611 m_delphi[iphoton] = dphi;
612 m_delthe[iphoton] = dthe;
613 m_delang[iphoton] = dang;
614 if ( energy < m_energyThreshold ) continue;
615 if ( getTime > m_gammathCut || getTime < m_gammatlCut ) continue;
616 // if((fabs(dthe) < m_gammaThetaCut) && (fabs(dphi)<m_gammaPhiCut) ) continue;
617 if ( dang < m_gammaTrkCut ) continue;
618 iphoton++;
619 iGam.push_back( i );
620 if ( iphoton >= 40 ) return StatusCode::SUCCESS;
621 }
622
623 int nGam = iGam.size();
624 m_nGam = nGam;
625 // std::cout << "num Good Photon " << m_nGam << " , "
626 // <<evtRecEvent->totalNeutral()<<std::endl;
627 // std::cout<<"dbg_4"<<std::endl;
628 counter[3]++;
629
630 double egam_ext = 0;
631 double ex_gam = 0;
632 double ey_gam = 0;
633 double ez_gam = 0;
634 double et_gam = 0;
635 double e_gam = 0;
636 for ( int i = 0; i < m_nGam; i++ )
637 {
638 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGam[i];
639 if ( !( *itTrk )->isEmcShowerValid() ) continue;
640 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
641 double eraw = emcTrk->energy();
642 double phi = emcTrk->phi();
643 double the = emcTrk->theta();
644 HepLorentzVector ptrk;
645 ex_gam += eraw * sin( the ) * cos( phi );
646 ey_gam += eraw * sin( the ) * sin( phi );
647 ez_gam += eraw * cos( the );
648 et_gam += eraw * sin( the );
649 e_gam += eraw;
650 if ( eraw >= egam_ext ) { egam_ext = eraw; }
651 }
652
653 double px_had = 0;
654 double py_had = 0;
655 double pz_had = 0;
656 double t_pxy2 = 0;
657 double pt_had = 0;
658 double p_had = 0;
659 double e_had = 0;
660 double p_max = 0.;
661 double e_max = 0.;
662 //
663 // check good charged track's infomation
664 //
665
666 double ctrk_theta = -10;
667 double ctrk_phi = -10;
668 for ( int i = 0; i < m_ngch; i++ )
669 {
670
671 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
672
673 if ( !( *itTrk )->isMdcTrackValid() ) continue; // MDC information
674 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
675
676 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
677 RecMdcKalTrack* mdcKalTrk = ( *itTrk )->mdcKalTrack();
678
679 // if ( m_ngch==2 &&mdcTrk->charge()>0) i = 0 ;
680 // if ( m_ngch==2 &&mdcTrk->charge()<0) i = 1 ;
681 m_charge[i] = mdcTrk->charge();
682 m_vx0[i] = mdcTrk->x();
683 m_vy0[i] = mdcTrk->y();
684 m_vz0[i] = mdcTrk->z();
685 m_px[i] = mdcTrk->px();
686 m_py[i] = mdcTrk->py();
687 m_pz[i] = mdcTrk->pz();
688 m_p[i] = mdcTrk->p();
689 ctrk_theta = mdcTrk->theta();
690 ctrk_phi = mdcTrk->phi();
691 mdcKalTrk->setPidType( RecMdcKalTrack::pion );
692 // double ptrk = mdcKalTrk->p() ;
693 m_kal_vx0[i] = mdcKalTrk->x();
694 m_kal_vy0[i] = mdcKalTrk->y();
695 m_kal_vz0[i] = mdcKalTrk->z();
696
697 m_kal_px[i] = mdcKalTrk->px();
698 m_kal_py[i] = mdcKalTrk->py();
699 m_kal_pz[i] = mdcKalTrk->pz();
700 // m_kal_p[i] = mdcKalTrk->p(); // pxy() and p() are not filled in the
701 // reconstruction algorithm
702 t_pxy2 = m_kal_px[i] * m_kal_px[i] + m_kal_py[i] * m_kal_py[i];
703 m_kal_p[i] = sqrt( t_pxy2 + m_kal_pz[i] * m_kal_pz[i] );
704 double ptrk = m_kal_p[i];
705 px_had += m_kal_px[i];
706 py_had += m_kal_py[i];
707 pz_had += m_kal_pz[i];
708 pt_had += sqrt( t_pxy2 );
709 p_had += m_kal_p[i];
710 e_had += sqrt( m_kal_p[i] * m_kal_p[i] + xmass[2] * xmass[2] );
711 if ( m_useDEDX && ( *itTrk )->isMdcDedxValid() )
712 { // DEDX information
713
714 RecMdcDedx* dedxTrk = ( *itTrk )->mdcDedx();
715 m_probPH[i] = dedxTrk->probPH();
716 m_normPH[i] = dedxTrk->normPH();
717
718 m_chie[i] = dedxTrk->chiE();
719 m_chimu[i] = dedxTrk->chiMu();
720 m_chipi[i] = dedxTrk->chiPi();
721 m_chik[i] = dedxTrk->chiK();
722 m_chip[i] = dedxTrk->chiP();
723 m_ghit[i] = dedxTrk->numGoodHits();
724 m_thit[i] = dedxTrk->numTotalHits();
725 }
726
727 if ( ( *itTrk )->isEmcShowerValid() )
728 {
729
730 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
731 m_e_emc[i] = emcTrk->energy();
732 m_phi_emc[i] = emcTrk->phi();
733 m_theta_emc[i] = emcTrk->theta();
734 if ( m_e_emc[i] > e_max )
735 {
736 p_max = m_p[i];
737 e_max = m_e_emc[i];
738 }
739 }
740 else
741 {
742 m_e_emc[i] = 0;
743 m_phi_emc[i] = -10;
744 m_theta_emc[i] = -10;
745 }
746
747 if ( m_useMUC && ( *itTrk )->isMucTrackValid() )
748 {
749
750 RecMucTrack* mucTrk = ( *itTrk )->mucTrack();
751 m_nhit_muc[i] = mucTrk->numHits();
752 m_nlay_muc[i] = mucTrk->numLayers();
753 }
754
755 if ( m_useTOF && ( *itTrk )->isTofTrackValid() )
756 { // TOF information
757
758 SmartRefVector<RecTofTrack> tofTrkCol = ( *itTrk )->tofTrack();
759
760 SmartRefVector<RecTofTrack>::iterator iter_tof = tofTrkCol.begin();
761
762 for ( ; iter_tof != tofTrkCol.end(); iter_tof++ )
763 {
764 TofHitStatus* status = new TofHitStatus;
765 status->setStatus( ( *iter_tof )->status() );
766
767 if ( !( status->is_barrel() ) )
768 { // endcap
769 if ( ( status->is_cluster() ) ) m_t_etof[i] = ( *iter_tof )->tof();
770 if ( !( status->is_counter() ) )
771 {
772 if ( status ) delete status;
773 continue;
774 } // ?
775 if ( status->layer() != 0 )
776 {
777 if ( status ) delete status;
778 continue;
779 } // layer1
780 double path = ( *iter_tof )->path(); // ?
781 double tof = ( *iter_tof )->tof();
782 double ph = ( *iter_tof )->ph();
783 double rhit = ( *iter_tof )->zrhit();
784 double qual = 0.0 + ( *iter_tof )->quality();
785 double cntr = 0.0 + ( *iter_tof )->tofID();
786 double texp[5];
787 for ( int j = 0; j < 5; j++ )
788 {
789 double gb = ptrk / xmass[j];
790 double beta = gb / sqrt( 1 + gb * gb );
791 texp[j] = path / beta / velc;
792 }
793
794 m_qual_etof[i] = qual;
795 m_tof_etof[i] = tof;
796 }
797 else
798 { // barrel
799 if ( ( status->is_cluster() ) ) m_t_btof[i] = ( *iter_tof )->tof();
800 if ( !( status->is_counter() ) )
801 {
802 if ( status ) delete status;
803 continue;
804 } // ?
805 if ( status->layer() == 1 )
806 { // layer1
807 double path = ( *iter_tof )->path(); // ?
808 double tof = ( *iter_tof )->tof();
809 double ph = ( *iter_tof )->ph();
810 double rhit = ( *iter_tof )->zrhit();
811 double qual = 0.0 + ( *iter_tof )->quality();
812 double cntr = 0.0 + ( *iter_tof )->tofID();
813 double texp[5];
814 for ( int j = 0; j < 5; j++ )
815 {
816 double gb = ptrk / xmass[j];
817 double beta = gb / sqrt( 1 + gb * gb );
818 texp[j] = path / beta / velc;
819 }
820
821 m_qual_btof1[i] = qual;
822 m_tof_btof1[i] = tof;
823 }
824
825 if ( status->layer() == 2 )
826 { // layer2
827 double path = ( *iter_tof )->path(); // ?
828 double tof = ( *iter_tof )->tof();
829 double ph = ( *iter_tof )->ph();
830 double rhit = ( *iter_tof )->zrhit();
831 double qual = 0.0 + ( *iter_tof )->quality();
832 double cntr = 0.0 + ( *iter_tof )->tofID();
833 double texp[5];
834 for ( int j = 0; j < 5; j++ )
835 {
836 double gb = ptrk / xmass[j];
837 double beta = gb / sqrt( 1 + gb * gb );
838 texp[j] = path / beta / velc;
839 }
840
841 m_qual_btof2[i] = qual;
842 m_tof_btof2[i] = tof;
843 }
844 }
845 if ( status ) delete status;
846 }
847 }
848 }
849
850 // tag
851
852 m_hadrontag = 0;
853 FastVertexFit* fvtxfit = FastVertexFit::instance();
854 if ( m_ngch != 2 ) m_hadrontag = 11111;
855 else if ( m_ngch == 2 && nCharge == 0 )
856 {
857 EvtRecTrackIterator itTrk1;
858
859 EvtRecTrackIterator itTrk2;
860
861 RecMdcKalTrack* mdcKalTrk1{ nullptr };
862
863 RecMdcKalTrack* mdcKalTrk2{ nullptr };
864
865 HepLorentzVector p41e, p42e, p4le;
866 Hep3Vector p31e, p32e, p3le;
867 HepLorentzVector p41m, p42m, p4lm;
868 Hep3Vector p31m, p32m, p3lm;
869 HepLorentzVector p41h, p42h, p4lh;
870 Hep3Vector p31h, p32h, p3lh;
871 WTrackParameter w1_ini, w1_ve, w1_vmu;
872 WTrackParameter w2_ini, w2_ve, w2_vmu;
873 int iip = 0;
874 int iim = 0;
875 for ( int i = 0; i < m_ngch; i++ )
876 {
877 if ( m_charge[i] > 0 ) itTrk1 = evtRecTrkCol->begin() + iGood[i];
878 if ( m_charge[i] < 0 ) itTrk2 = evtRecTrkCol->begin() + iGood[i];
879 if ( m_charge[i] > 0 ) mdcKalTrk1 = ( *itTrk1 )->mdcKalTrack();
880 if ( m_charge[i] < 0 ) mdcKalTrk2 = ( *itTrk2 )->mdcKalTrack();
881 if ( m_charge[i] > 0 ) iip = i;
882 if ( m_charge[i] < 0 ) iim = i;
883
884 if ( m_charge[i] > 0 )
885 w1_ini = WTrackParameter( xmass[2], mdcKalTrk1->helix(), mdcKalTrk1->err() );
886 if ( m_charge[i] < 0 )
887 w2_ini = WTrackParameter( xmass[2], mdcKalTrk2->helix(), mdcKalTrk2->err() );
888 if ( m_charge[i] > 0 ) p41h = w1_ini.p();
889 if ( m_charge[i] < 0 ) p42h = w2_ini.p();
890 if ( m_charge[i] > 0 ) p41h.boost( u_cms );
891 if ( m_charge[i] < 0 ) p42h.boost( u_cms );
892 if ( m_charge[i] > 0 ) p31h = p41h.vect();
893 if ( m_charge[i] < 0 ) p32h = p42h.vect();
894
895 if ( m_charge[i] > 0 ) p41e = w1_ini.p();
896 if ( m_charge[i] < 0 ) p42e = w2_ini.p();
897 if ( m_charge[i] > 0 ) p41e.boost( u_cms );
898 if ( m_charge[i] < 0 ) p42e.boost( u_cms );
899 if ( m_charge[i] > 0 ) p31e = p41e.vect();
900 if ( m_charge[i] < 0 ) p32e = p42e.vect();
901
902 if ( m_charge[i] > 0 )
903 {
904 m_px_cms_ep = p41h.px();
905 m_py_cms_ep = p41h.py();
906 m_pz_cms_ep = p41h.pz();
907 m_e_cms_ep = p41h.e();
908 }
909 if ( m_charge[i] < 0 )
910 {
911 m_px_cms_em = p42h.px();
912 m_py_cms_em = p42h.py();
913 m_pz_cms_em = p42h.pz();
914 m_e_cms_em = p42h.e();
915 }
916 }
917 double e01 = m_e_emc[iip]; // m_e_cms_ep;
918 double e02 = m_e_emc[iim]; // m_e_cms_em;
919
920 int ilarge = ( e01 > e02 ) ? iip : iim;
921
922 p4lh = ( e01 > e02 ) ? p41h : p42h;
923
924 p3lh = ( e01 > e02 ) ? p31h : p32h;
925
926 double acollh = 180. - p31h.angle( p32h ) * 180.0 / CLHEP::pi;
927 double acoplh = 180. - ( p31h.perpPart() ).angle( p32h.perpPart() ) * 180.0 / CLHEP::pi;
928 double poeb1h = p41h.rho() / beamEnergy;
929 double poeb2h = p42h.rho() / beamEnergy;
930 double poeblh = p4lh.rho() / beamEnergy;
931
932 double eoeb1 = m_e_emc[iip] / beamEnergy;
933 double eoeb2 = m_e_emc[iim] / beamEnergy;
934 double eop1 = 0;
935 if ( p41h.rho() > 0 ) eop1 = m_e_emc[iip] / p41h.rho();
936 double eop2 = 0;
937 if ( p42h.rho() > 0 ) eop2 = m_e_emc[iim] / p42h.rho();
938
939 double eope1 = 0;
940 if ( p41e.rho() > 0 ) eope1 = m_e_emc[iip] / p41e.rho();
941 double eope2 = 0;
942 if ( p42e.rho() > 0 ) eope2 = m_e_emc[iim] / p42e.rho();
943 double eopm1 = 0;
944 if ( p41m.rho() > 0 ) eopm1 = m_e_emc[iip] / p41m.rho();
945 double eopm2 = 0;
946 if ( p42m.rho() > 0 ) eopm2 = m_e_emc[iim] / p42m.rho();
947
948 double exoeb1 =
949 m_e_emc[iip] * sin( m_theta_emc[iip] ) * cos( m_phi_emc[iip] ) / beamEnergy;
950 double eyoeb1 =
951 m_e_emc[iip] * sin( m_theta_emc[iip] ) * sin( m_phi_emc[iip] ) / beamEnergy;
952 double ezoeb1 = m_e_emc[iip] * cos( m_theta_emc[iip] ) / beamEnergy;
953 double etoeb1 = m_e_emc[iip] * sin( m_theta_emc[iip] ) / beamEnergy;
954
955 double exoeb2 =
956 m_e_emc[iim] * sin( m_theta_emc[iim] ) * cos( m_phi_emc[iim] ) / beamEnergy;
957 double eyoeb2 =
958 m_e_emc[iim] * sin( m_theta_emc[iim] ) * sin( m_phi_emc[iim] ) / beamEnergy;
959 double ezoeb2 = m_e_emc[iim] * cos( m_theta_emc[iim] ) / beamEnergy;
960 double etoeb2 = m_e_emc[iim] * sin( m_theta_emc[iim] ) / beamEnergy;
961
962 double eoebl = m_e_emc[ilarge] / beamEnergy;
963
964 double eopl = 0;
965 if ( p4lh.rho() > 0 ) eopl = m_e_emc[ilarge] / p4lh.rho();
966
967 double exoebl =
968 m_e_emc[ilarge] * sin( m_theta_emc[ilarge] ) * cos( m_phi_emc[ilarge] ) / beamEnergy;
969 double eyoebl =
970 m_e_emc[ilarge] * sin( m_theta_emc[ilarge] ) * sin( m_phi_emc[ilarge] ) / beamEnergy;
971 double ezoebl = m_e_emc[ilarge] * cos( m_theta_emc[ilarge] ) / beamEnergy;
972 double etoebl = m_e_emc[ilarge] * sin( m_theta_emc[ilarge] ) / beamEnergy;
973
974 int mucinfo1 = ( m_nhit_muc[iip] >= 2 && m_nlay_muc[iip] >= 2 ) ? 1 : 0;
975 int mucinfo2 = ( m_nhit_muc[iim] >= 2 && m_nlay_muc[iim] >= 2 ) ? 1 : 0;
976 int mucinfol = ( m_nhit_muc[ilarge] >= 2 && m_nlay_muc[ilarge] >= 2 ) ? 1 : 0;
977 int pidel = ( e01 > e02 ) ? m_nep : m_nem;
978 int pidmul = ( e01 > e02 ) ? m_nmup : m_nmum;
979 double deltatof = 0.0;
980
981 // if(m_tof_btof2[iip]*m_tof_btof2[iim]!=0)
982 // deltatof+=fabs(m_tof_btof2[iip]-m_tof_btof2[iim]);
983 // if(m_tof_btof1[iip]*m_tof_btof1[iim]!=0)deltatof+=fabs(m_tof_btof1[iip]-m_tof_btof1[iim]);
984 // if(m_tof_etof[iip]*m_tof_etof[iim]!=0)deltatof+=fabs(m_tof_etof[iip]-m_tof_etof[iim]);
985
986 // if(!m_endcap) {
987 if ( m_t_btof[iip] * m_t_btof[iim] != 0 ) deltatof = fabs( m_t_btof[iip] - m_t_btof[iim] );
988 // }
989 // else {
990 // if(m_t_etof[iip]*m_t_etof[iim]!=0)deltatof=fabs(m_t_etof[iip]-m_t_etof[iim]);
991 // }
992
993 // if (acollh>m_acoll_h_cut)m_hadrontag+=1;
994 if ( ( acollh > m_acoll_h_cut ) || ( !m_useEMC || m_nGam >= m_ngam_h_cut ) )
995 m_hadrontag += 11;
996 if ( !m_useTOF || ( deltatof < m_dtof_h_cut ) ) m_hadrontag += 100;
997 if ( !m_useMUC || ( mucinfo1 == 0 || mucinfo2 == 0 ) ) m_hadrontag += 1000;
998 if ( !m_useEMC || ( fabs( eope1 - 1 ) > m_eop_h_cut && fabs( eope2 - 1 ) > m_eop_h_cut ) )
999 m_hadrontag += 10000;
1000
1001 m_deltatof = deltatof;
1002 m_eop1 = eope1;
1003 m_eop2 = eope2;
1004 m_eoeb1 = eoeb1;
1005 m_eoeb2 = eoeb2;
1006
1007 m_etoeb1 = etoeb1;
1008 m_etoeb2 = etoeb2;
1009 m_mucinfo1 = mucinfo1;
1010 m_mucinfo2 = mucinfo2;
1011
1012 m_acoll = acollh;
1013 m_acopl = acoplh;
1014 m_poeb1 = poeb1h;
1015 m_poeb2 = poeb2h;
1016 m_cos_ep = p41h.cosTheta();
1017 m_cos_em = p42h.cosTheta();
1018 m_mass_ee = ( p41h + p42h ).m();
1019 }
1020 double br = 0;
1021 double bz = 0;
1022 double thr = 0;
1023 double evis = 0;
1024 WTrackParameter w1_vh, w2_vh, w3_vh;
1025
1026 br = sqrt( ( px_had + ex_gam ) * ( px_had + ex_gam ) +
1027 ( py_had + ey_gam ) * ( py_had + ey_gam ) ) /
1028 ( pt_had + et_gam );
1029 bz = fabs( pz_had + ez_gam ) / ( p_had + e_gam );
1030 thr = p_had + e_gam;
1031 evis = e_had + e_gam;
1032 log << MSG::DEBUG << "hadron " << px_had << " " << py_had << " " << pz_had << " " << pt_had
1033 << " " << p_had << " " << br << " " << bz << endmsg;
1034 log << MSG::DEBUG << "Gamma " << ex_gam << " " << ey_gam << " " << ez_gam << " " << e_gam
1035 << " " << thr / beamEnergy << endmsg;
1036 if ( !m_useEMC || ( ( br < m_br_h_cut ) && ( bz < m_bz_h_cut ) ) ) m_hadrontag += 100000;
1037 if ( !m_useEMC || thr / beamEnergy > m_thr_h_cut ) m_hadrontag += 1000000;
1038 m_emax = egam_ext;
1039 m_esum = e_gam;
1040 m_br = br;
1041 m_bz = bz;
1042 m_thr = thr;
1043 m_evis = evis;
1044 log << MSG::INFO << "m_hadrontag= " << m_hadrontag << endmsg;
1045 // std::cout<<"m_sbhabhatag= "<<m_sbhabhatag<<std::endl;
1046 // std::cout<<"m_sdimutag= "<<m_sdimutag<<std::endl;
1047 // std::cout<<"m_hadrontag= "<<m_hadrontag<<std::endl;
1048 if ( m_hadrontag == 1111111 )
1049 {
1050 nhadron++;
1051 for ( int i = 0; i < m_ngch; i++ )
1052 {
1053 // m_ha_costheta->Fill(cos(m_theta[i]));
1054 // m_ha_phi->Fill(m_phi[i]);
1055 m_ha_costheta->Fill( ctrk_theta );
1056 m_ha_phi->Fill( ctrk_phi );
1057 }
1058 if ( m_ngch >= 3 )
1059 {
1060 RecMdcKalTrack* ktrk0 = ( *( evtRecTrkCol->begin() + iGood[0] ) )->mdcKalTrack();
1061 RecMdcKalTrack* ktrk1 = ( *( evtRecTrkCol->begin() + iGood[1] ) )->mdcKalTrack();
1062 RecMdcKalTrack* ktrk2 = ( *( evtRecTrkCol->begin() + iGood[2] ) )->mdcKalTrack();
1063 // w1_vh=WTrackParameter (xmass[2],ktrk0->getZHelix(),ktrk0->getZError());
1064 // w2_vh=WTrackParameter (xmass[2],ktrk1->getZHelix(),ktrk1->getZError());
1065 // w3_vh=WTrackParameter (xmass[2],ktrk2->getZHelix(),ktrk2->getZError());
1066 // vtxfit->init();
1067 // vtxfit->AddTrack(0, w1_vh);
1068 // vtxfit->AddTrack(1, w2_vh);
1069 // vtxfit->AddTrack(2, w3_vh);
1070 // vtxfit->AddVertex(0, vxpar,0, 1);
1071 // if(vtxfit->Fit(0)) {
1072 // m_ha_vx->Fill((vtxfit->vx(0)).x());
1073 // m_ha_vy->Fill((vtxfit->vx(0)).y());
1074 // m_ha_vz->Fill((vtxfit->vx(0)).z());
1075 // }
1076
1077 // ktrk0->setPidType(RecMdcKalTrack::pion);
1078 // ktrk1->setPidType(RecMdcKalTrack::pion);
1079 // ktrk2->setPidType(RecMdcKalTrack::pion);
1080 fvtxfit->init();
1081 fvtxfit->addTrack( 0, ktrk0->helix(), ktrk0->err() );
1082 fvtxfit->addTrack( 1, ktrk1->helix(), ktrk1->err() );
1083 fvtxfit->addTrack( 2, ktrk2->helix(), ktrk2->err() );
1084 if ( fvtxfit->Fit() )
1085 {
1086 m_ha_vx->Fill( ( fvtxfit->Vx() )[0] );
1087 m_ha_vy->Fill( ( fvtxfit->Vx() )[1] );
1088 m_ha_vz->Fill( ( fvtxfit->Vx() )[2] );
1089 }
1090 }
1091
1092 m_ha_br->Fill( br );
1093 m_ha_bz->Fill( bz );
1094 m_ha_pmax->Fill( p_max );
1095 m_ha_emax->Fill( e_max );
1096 m_ha_etot->Fill( evis );
1097 m_ha_nneu->Fill( nGam );
1098 m_ha_nchg->Fill( m_ngch );
1099 if ( m_ngch == 0 ) { n0prong++; }
1100 if ( m_ngch == 2 && nCharge == 0 ) { n2prong++; }
1101 if ( m_ngch == 4 && nCharge == 0 ) { n4prong++; }
1102
1103 if ( m_ngch > 4 ) { ng4prong++; }
1104 if ( m_writentuple ) m_tuple1->write().ignore();
1105 setFilterPassed( true );
1106 }
1107
1108 return StatusCode::SUCCESS;
1109}
1110
1111// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
1112StatusCode DQASelHadron::finalize() {
1113
1114 MsgStream log( msgSvc(), name() );
1115 log << MSG::INFO << "in finalize()" << endmsg;
1116 log << MSG::INFO << counter << endmsg;
1117 log << MSG::ALWAYS << "nhadron = " << nhadron << endmsg;
1118 log << MSG::INFO << "n0prong = " << n0prong << endmsg;
1119 log << MSG::INFO << "n2prong = " << n2prong << endmsg;
1120 log << MSG::INFO << "n4prong = " << n4prong << endmsg;
1121 log << MSG::INFO << "ng4prong = " << ng4prong << endmsg;
1122 return StatusCode::SUCCESS;
1123}
DECLARE_COMPONENT
DECLARE_COMPONENT(BesBdkRc)
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition CosmicGenerator.h:36
u_cms
const Hep3Vector u_cms
Definition DQADtagAlg.cxx:56
p_cms
const HepLorentzVector p_cms(0.034067, 0.0, 0.0, 3.097)
n
const Int_t n
Definition DataBase/tau_mode.c:57
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
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
energy
************Class m_ypar INTEGER m_KeyWgt INTEGER m_KeyIHVP INTEGER m_KeyGPS INTEGER m_IsBeamPolarized INTEGER m_EvtGenInterface DOUBLE PRECISION m_Emin DOUBLE PRECISION m_sphot DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_q2 DOUBLE PRECISION m_PolBeam2 DOUBLE PRECISION m_xErrPb *COMMON c_KK2f $ !CMS energy average $ !Spin Polarization vector first beam $ !Spin Polarization vector second beam $ !Beam energy spread[GeV] $ !minimum hadronization energy[GeV] $ !input READ never touch them !$ !debug facility $ !maximum weight $ !inverse alfaQED $ !minimum real photon energy
Definition KK2f.h:50
ntupleSvc
INTupleSvc * ntupleSvc()
Definition ServiceAccessor.h:51
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:12
DQASelHadron::initialize
StatusCode initialize()
Definition DQASelHadron.cxx:104
DQASelHadron::finalize
StatusCode finalize()
Definition DQASelHadron.cxx:1112
DQASelHadron::execute
StatusCode execute()
Definition DQASelHadron.cxx:316
DQASelHadron::DQASelHadron
DQASelHadron(const std::string &name, ISvcLocator *pSvcLocator)
Definition DQASelHadron.cxx:64
DstEmcShower::dy
double dy() const
Definition DstEmcShower.cxx:10
DstEmcShower::dz
double dz() const
Definition DstEmcShower.cxx:17
DstEmcShower::dx
double dx() const
Definition DstEmcShower.cxx:3
DstMdcTrack::err
const HepSymMatrix err() const
Definition DstMdcTrack.cxx:149
DstMdcTrack::helix
const HepVector helix() const
......
Definition DstMdcTrack.cxx:147
EmcID::barrel_ec
static unsigned int barrel_ec(const Identifier &id)
Values of different levels (failure returns 0).
Definition EmcID.cxx:36
FastVertexFit::instance
static FastVertexFit * instance()
Definition FastVertexFit.cxx:9
FastVertexFit::addTrack
void addTrack(const int number, const HepVector helix, const HepSymMatrix err)
Definition FastVertexFit.cxx:29
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
VFHelix::pivot
const HepPoint3D & pivot(void) const
returns pivot position.
Definition Analysis/VertexFit/include/VertexFit/Helix.h:177
VFHelix::a
const HepVector & a(void) const
returns helix parameters.
Definition Analysis/VertexFit/include/VertexFit/Helix.h:195