BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
ValidJpsiRhopi.cxx
Go to the documentation of this file.
1#include "GaudiKernel/IDataProviderSvc.h"
2#include "GaudiKernel/ISvcLocator.h"
3#include "GaudiKernel/MsgStream.h"
4#include "GaudiKernel/PropertyMgr.h"
5#include "GaudiKernel/SmartDataPtr.h"
6
7#include "EventModel/Event.h"
8#include "EventModel/EventModel.h"
9
10#include "DstEvent/TofHitStatus.h"
11#include "EventModel/EventHeader.h"
12#include "EvtRecEvent/EvtRecEvent.h"
13#include "EvtRecEvent/EvtRecTrack.h"
14
15#include "CLHEP/Vector/LorentzVector.h"
16#include "CLHEP/Vector/ThreeVector.h"
17#include "CLHEP/Vector/TwoVector.h"
18#include "GaudiKernel/Bootstrap.h"
19#include "GaudiKernel/IHistogramSvc.h"
20#include "GaudiKernel/INTupleSvc.h"
21#include "GaudiKernel/NTuple.h"
22#include "TMath.h"
23using CLHEP::Hep2Vector;
24using CLHEP::Hep3Vector;
25using CLHEP::HepLorentzVector;
26#include "CLHEP/Geometry/Point3D.h"
27#ifndef ENABLE_BACKWARDS_COMPATIBILITY
28typedef HepGeom::Point3D<double> HepPoint3D;
29#endif
30
31#include "ParticleID/ParticleID.h"
32#include "VertexFit/KinematicFit.h"
33#include "VertexFit/VertexFit.h"
34
35#include "ValidJpsiRhopi.h"
36
37#include <vector>
38
39const double mpi0c = 0.1349766;
40const double mpi = 0.13957;
41const double ecms = 3.097;
42const double mk = 0.493677;
43const double xmass[5] = { 0.000511, 0.105658, 0.139570, 0.493677, 0.938272 };
44// e mu pi K p
45
46typedef std::vector<int> Vint;
47typedef std::vector<HepLorentzVector> Vp4;
48
49/////////////////////////////////////////////////////////////////////////////
50// //
51// J/psi-->rho + pi //
52// |-->pi pi //
53// //
54// Final : (2Gamma + 2pi) xugf 2008.09 //
55/////////////////////////////////////////////////////////////////////////////
56DECLARE_COMPONENT( ValidJpsiRhopi )
57ValidJpsiRhopi::ValidJpsiRhopi( const std::string& name, ISvcLocator* pSvcLocator )
58 : Algorithm( name, pSvcLocator ) {
59
60 m_tuple4 = 0;
61 m_tuple5 = 0;
62
63 for ( int i = 0; i < 12; i++ ) m_pass[i] = 0;
64
65 // Declare the properties
66 declareProperty( "Vr0cut", m_vr0cut = 5.0 );
67 declareProperty( "Vz0cut", m_vz0cut = 15.0 );
68 declareProperty( "EnergyThreshold", m_energyThreshold = 0.03 );
69 declareProperty( "GammaPhiCut", m_gammaPhiCut = 20.0 );
70 declareProperty( "GammaThetaCut", m_gammaThetaCut = 20.0 );
71}
72
73// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
74StatusCode ValidJpsiRhopi::initialize() {
75 MsgStream log( msgSvc(), name() );
76
77 log << MSG::INFO << "in initialize()" << endmsg;
78
79 StatusCode status;
80
81 status = service( "THistSvc", m_thistsvc );
82 if ( status.isFailure() )
83 {
84 log << MSG::INFO << "Unable to retrieve pointer to THistSvc" << endmsg;
85 return status;
86 }
87 m_pi_vx = new TH1F( "pi_vx", "pi_vx", 100, -1.0, 1.0 );
88 status = m_thistsvc->regHist( "/VAL/PHY/pi_vx", m_pi_vx );
89 m_pi_vy = new TH1F( "pi_vy", "pi_vy", 100, -1.0, 1.0 );
90 status = m_thistsvc->regHist( "/VAL/PHY/pi_vy", m_pi_vy );
91 m_pi_vz = new TH1F( "pi_vz", "pi_vz", 100, -6.0, 6.0 );
92 status = m_thistsvc->regHist( "/VAL/PHY/pi_vz", m_pi_vz );
93 m_pip_mom = new TH1F( "pip_mom", "pip_moment", 100, 0.0, 1.6 );
94 status = m_thistsvc->regHist( "/VAL/PHY/pip_mom", m_pip_mom );
95 m_pim_mom = new TH1F( "pim_mom", "pim_moment", 100, 0.0, 1.6 );
96 status = m_thistsvc->regHist( "/VAL/PHY/pim_mom", m_pim_mom );
97 m_rhop_mass = new TH1F( "rhop_mass", "rhop_mass", 100, 0.0, 1.5 );
98 status = m_thistsvc->regHist( "/VAL/PHY/rhop_mass", m_rhop_mass );
99 m_rhom_mass = new TH1F( "rhom_mass", "rhom_mass", 100, 0.0, 1.5 );
100 status = m_thistsvc->regHist( "/VAL/PHY/rhom_mass", m_rhom_mass );
101 m_rho0_mass = new TH1F( "rho0_mass", "rho0_mass", 100, 0.0, 1.5 );
102 status = m_thistsvc->regHist( "/VAL/PHY/rho0_mass", m_rho0_mass );
103 m_pi0_mass = new TH1F( "pi0_mass", "pi0_mass", 100, 0.0, 0.5 );
104 status = m_thistsvc->regHist( "/VAL/PHY/pi0_mass", m_pi0_mass );
105 m_chisq_4c = new TH1F( "chisq_4c", "chisq_4c", 100, 0.0, 150. );
106 status = m_thistsvc->regHist( "/VAL/PHY/chisq_4c", m_chisq_4c );
107 m_chisq_5c = new TH1F( "chisq_5c", "chisq_5c", 100, 0.0, 150. );
108 status = m_thistsvc->regHist( "/VAL/PHY/chisq_5c", m_chisq_5c );
109
110 m_cos_pip = new TH1F( "cos_pip", "cos_pip", 100, -1.0, 1.0 );
111 status = m_thistsvc->regHist( "/VAL/PHY/cos_pip", m_cos_pip );
112 m_cos_pim = new TH1F( "cos_pim", "cos_pim", 100, -1.0, 1.0 );
113 status = m_thistsvc->regHist( "/VAL/PHY/cos_pim", m_cos_pim );
114
115 m_chipi_dedx = new TH1F( "chipi_dedx", "chipi_dedx", 60, -5.0, 10. );
116 status = m_thistsvc->regHist( "/VAL/PHY/chipi_dedx", m_chipi_dedx );
117 m_chie_dedx = new TH1F( "chie_dedx", "chie_dedx", 60, -15.0, 5. );
118 status = m_thistsvc->regHist( "/VAL/PHY/chie_dedx", m_chie_dedx );
119 m_chimu_dedx = new TH1F( "chimu_dedx", "chimu_dedx", 60, -5.0, 10. );
120 status = m_thistsvc->regHist( "/VAL/PHY/chimu_dedx", m_chimu_dedx );
121 m_chik_dedx = new TH1F( "chik_dedx", "chik_dedx", 100, -20.0, 10. );
122 status = m_thistsvc->regHist( "/VAL/PHY/chik_dedx", m_chik_dedx );
123 m_chip_dedx = new TH1F( "chip_dedx", "chip_dedx", 100, -20.0, 5. );
124 status = m_thistsvc->regHist( "/VAL/PHY/chip_dedx", m_chip_dedx );
125
126 NTuplePtr nt4( ntupleSvc(), "FILE1/h4" );
127 if ( nt4 ) m_tuple4 = nt4;
128 else
129 {
130 m_tuple4 = ntupleSvc()->book( "FILE1/h4", CLID_ColumnWiseTuple, "4gam6pi Ntuple" );
131 if ( m_tuple4 )
132 {
133 status = m_tuple4->addItem( "ngam", m_ngoodneu );
134 status = m_tuple4->addItem( "npip", m_npip );
135 status = m_tuple4->addItem( "npim", m_npim );
136 status = m_tuple4->addItem( "ngoodch", m_ngoodch );
137 status = m_tuple4->addItem( "chisq4c", m_chisq4c );
138 status = m_tuple4->addItem( "ppi04c", m_ppi0 );
139 status = m_tuple4->addItem( "mpi04c", m_mpi0 );
140 status = m_tuple4->addItem( "chisq5c", m_chisq5c );
141 status = m_tuple4->addItem( "ppi05c", m_ppi0fit );
142 status = m_tuple4->addItem( "mpi05c", m_mpi0fit );
143 status = m_tuple4->addItem( "g1inpi0the", m_g1inpi0the );
144 status = m_tuple4->addItem( "g2inpi0the", m_g2inpi0the );
145 status = m_tuple4->addItem( "theta2pi", m_theta2pi );
146 status = m_tuple4->addItem( "ppip", m_ppip );
147 status = m_tuple4->addItem( "ppim", m_ppim );
148 status = m_tuple4->addItem( "p2pi", m_p2pi );
149 status = m_tuple4->addItem( "m2pi", m_m2pi );
150 status = m_tuple4->addItem( "ppip0", m_ppip0 );
151 status = m_tuple4->addItem( "mpip0", m_mpip0 );
152 status = m_tuple4->addItem( "ppim0", m_ppim0 );
153 status = m_tuple4->addItem( "mpim0", m_mpim0 );
154 status = m_tuple4->addItem( "eneumiss", m_eneumiss );
155 status = m_tuple4->addItem( "pneumiss", m_pneumiss );
156 status = m_tuple4->addItem( "mneumiss", m_mneumiss );
157 // status = m_tuple4->addIndexedItem ("kal_cos", m_ngoodch , m_kal_cos);
158 }
159 else
160 {
161 log << MSG::ERROR << " Cannot book N-tuple:" << long( m_tuple4 ) << endmsg;
162 return StatusCode::FAILURE;
163 }
164 }
165
166 NTuplePtr nt5( ntupleSvc(), "FILE1/h5" );
167 if ( nt5 ) m_tuple5 = nt5;
168 else
169 {
170 m_tuple5 = ntupleSvc()->book( "FILE1/h5", CLID_ColumnWiseTuple, "4gam6pi Ntuple" );
171 if ( m_tuple5 )
172 {
173 status = m_tuple5->addItem( "m2graw", m_m2graw );
174 status = m_tuple5->addItem( "emiss", m_emiss );
175 status = m_tuple5->addItem( "pmiss", m_pmiss );
176 status = m_tuple5->addItem( "mmiss", m_mmiss );
177 status = m_tuple5->addItem( "prho0", m_prho0 );
178 status = m_tuple5->addItem( "mrho0", m_mrho0 );
179 status = m_tuple5->addItem( "pmrho0", m_pmrho0 );
180 status = m_tuple5->addItem( "mmrho0", m_mmrho0 );
181 status = m_tuple5->addItem( "prhop", m_prhop );
182 status = m_tuple5->addItem( "mrhop", m_mrhop );
183 status = m_tuple5->addItem( "pmrhom", m_pmrhom );
184 status = m_tuple5->addItem( "mmrhom", m_mmrhom );
185 status = m_tuple5->addItem( "prhom", m_prhom );
186 status = m_tuple5->addItem( "ppipraw", m_ppipraw );
187 status = m_tuple5->addItem( "mrhom", m_mrhom );
188 }
189 else
190 {
191 log << MSG::ERROR << " Cannot book N-tuple:" << long( m_tuple4 ) << endmsg;
192 return StatusCode::FAILURE;
193 }
194 }
195
196 //
197 //--------end of book--------
198 //
199
200 log << MSG::INFO << "successfully return from initialize()" << endmsg;
201 return StatusCode::SUCCESS;
202}
203
204// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
205StatusCode ValidJpsiRhopi::execute() {
206
207 StatusCode sc = StatusCode::SUCCESS;
208
209 MsgStream log( msgSvc(), name() );
210 // log << MSG::INFO << "in execute()" << endmsg;
211
212 m_pass[0] += 1;
213
214 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
215 // SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::Recon::EvtRecEvent);
216 SmartDataPtr<EvtRecEvent> evtRecEvent( eventSvc(), EventModel::EvtRec::EvtRecEvent );
217
218 // SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::Recon::EvtRecTrackCol);
219 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol( eventSvc(), EventModel::EvtRec::EvtRecTrackCol );
220 //
221 // check x0, y0, z0, r0
222 // suggest cut: |z0|<10 && r0<5
223 //
224
225 Vint ipip, ipim, iGood;
226 iGood.clear();
227 ipip.clear();
228 ipim.clear();
229
230 Vp4 ppip, ppim;
231 ppip.clear();
232 ppim.clear();
233
234 int TotCharge = 0;
235 for ( int i = 0; i < evtRecEvent->totalCharged(); i++ )
236 {
237 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
238 if ( !( *itTrk )->isMdcTrackValid() ) continue;
239 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
240
241 m_pi_vx->Fill( mdcTrk->x() );
242 m_pi_vy->Fill( mdcTrk->y() );
243 m_pi_vz->Fill( mdcTrk->z() );
244
245 if ( fabs( mdcTrk->z() ) >= m_vz0cut ) continue;
246 if ( mdcTrk->r() >= m_vr0cut ) continue;
247 iGood.push_back( ( *itTrk )->trackId() );
248 TotCharge += mdcTrk->charge();
249 }
250 //
251 // Finish Good Charged Track Selection
252 //
253 int nGood = iGood.size();
254 m_ngoodch = nGood;
255
256 Vint iGam;
257 iGam.clear();
258 for ( int i = evtRecEvent->totalCharged(); i < evtRecEvent->totalTracks(); i++ )
259 {
260 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
261 if ( !( *itTrk )->isEmcShowerValid() ) continue;
262 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
263 Hep3Vector emcpos( emcTrk->x(), emcTrk->y(), emcTrk->z() );
264 // find the nearest charged track
265 double dthe = 200.;
266 double dphi = 200.;
267 double dang = 200.;
268 for ( int j = 0; j < evtRecEvent->totalCharged(); j++ )
269 {
270 EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() + j;
271 if ( !( *jtTrk )->isExtTrackValid() ) continue;
272 RecExtTrack* extTrk = ( *jtTrk )->extTrack();
273 if ( extTrk->emcVolumeNumber() == -1 ) continue;
274 Hep3Vector extpos = extTrk->emcPosition();
275 // double ctht = extpos.cosTheta(emcpos);
276 double angd = extpos.angle( emcpos );
277 double thed = extpos.theta() - emcpos.theta();
278 double phid = extpos.deltaPhi( emcpos );
279 thed = fmod( thed + CLHEP::twopi + CLHEP::twopi + pi, CLHEP::twopi ) - CLHEP::pi;
280 phid = fmod( phid + CLHEP::twopi + CLHEP::twopi + pi, CLHEP::twopi ) - CLHEP::pi;
281 angd = fmod( angd, CLHEP::twopi );
282
283 if ( angd < dang )
284 {
285 dang = angd;
286 dthe = thed;
287 dphi = phid;
288 }
289 }
290
291 if ( dang >= 200. ) continue;
292
293 double eraw = emcTrk->energy();
294 dthe = dthe * 180 / ( CLHEP::pi );
295 dphi = dphi * 180 / ( CLHEP::pi );
296 dang = dang * 180 / ( CLHEP::pi );
297 if ( eraw < m_energyThreshold ) continue;
298 // if((fabs(dthe) < m_gammaThetaCut) && (fabs(dphi)<m_gammaPhiCut) ) continue;
299 //
300 // good photon cut will be set here
301 //
302 iGam.push_back( ( *itTrk )->trackId() );
303 }
304 //
305 // Finish Good Photon Selection
306 //
307 int nGam = iGam.size();
308 m_ngoodneu = nGam;
309 //
310 // Assign 4-momentum to each photon
311 //
312
313 HepLorentzVector ptcharg, ptneu, ptchargp, ptchargm;
314 Vp4 pGam;
315 pGam.clear();
316 for ( int i = 0; i < nGam; i++ )
317 {
318 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGam[i];
319 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
320 double eraw = emcTrk->energy();
321 double phi = emcTrk->phi();
322 double the = emcTrk->theta();
323 HepLorentzVector ptrk;
324 ptrk.setPx( eraw * sin( the ) * cos( phi ) );
325 ptrk.setPy( eraw * sin( the ) * sin( phi ) );
326 ptrk.setPz( eraw * cos( the ) );
327 ptrk.setE( eraw );
328 pGam.push_back( ptrk );
329 }
330
331 //
332 // Dedx Check ================
333 //
334 for ( int i = 0; i < nGood; i++ )
335 {
336 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
337 if ( !( *itTrk )->isMdcTrackValid() ) continue;
338 if ( !( *itTrk )->isMdcDedxValid() ) continue;
339 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
340 RecMdcDedx* dedxTrk = ( *itTrk )->mdcDedx();
341 m_chie_dedx->Fill( dedxTrk->chiE() );
342 m_chimu_dedx->Fill( dedxTrk->chiMu() );
343 m_chipi_dedx->Fill( dedxTrk->chiPi() );
344 m_chik_dedx->Fill( dedxTrk->chiK() );
345 m_chip_dedx->Fill( dedxTrk->chiP() );
346 }
347
348 //
349 // Assign 4-momentum to each charged track
350 //
351
352 int TotQ = 0;
353
354 ParticleID* pid = ParticleID::instance();
355 for ( int i = 0; i < nGood; i++ )
356 {
357 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
358 // if(pid) delete pid;
359 pid->init();
360 pid->setMethod( pid->methodProbability() );
361 pid->setChiMinCut( 4 );
362 pid->setRecTrack( *itTrk );
363 pid->usePidSys( pid->useDedx() | pid->useTof1() | pid->useTof2() ); // use PID sub-system
364 pid->identify( pid->onlyPion() | pid->onlyKaon() ); // seperater Pion/Kaon
365 pid->calculate();
366 if ( !( pid->IsPidInfoValid() ) ) continue;
367
368 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
369 RecMdcKalTrack* mdcKalTrk = ( *itTrk )->mdcKalTrack();
370 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
371
372 TotQ += mdcKalTrk->charge();
373
374 // RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
375
376 // if(pid->probPion() < 0.001 || (pid->probPion() < pid->probKaon())) continue;
377 if ( pid->probPion() < pid->probKaon() ) continue;
378 mdcKalTrk->setPidType( RecMdcKalTrack::pion );
379 HepLorentzVector ptrk;
380 ptrk.setPx( mdcKalTrk->px() );
381 ptrk.setPy( mdcKalTrk->py() );
382 ptrk.setPz( mdcKalTrk->pz() );
383 double p3 = ptrk.mag();
384 ptrk.setE( sqrt( p3 * p3 + mpi * mpi ) );
385 if ( mdcTrk->charge() > 0 )
386 {
387 ipip.push_back( iGood[i] );
388 ppip.push_back( ptrk );
389 }
390 else
391 {
392 ipim.push_back( iGood[i] );
393 ppim.push_back( ptrk );
394 }
395 }
396
397 int npip = ipip.size();
398 int npim = ipim.size();
399 m_npip = npip;
400 m_npim = npim;
401
402 if ( nGood != 2 || TotCharge != 0 ) return sc;
403 m_pass[1] += 1;
404 if ( nGam < 2 ) return sc;
405 m_pass[2] += 1;
406
407 // if(npip != 1) return sc;
408 // if(npip*npim != 1) return sc;
409
410 HepLorentzVector BoostP( 0.0, 0.0, 0.0, ecms );
411 BoostP[0] = 2 * sin( 0.011 ) * ( ecms / 2 );
412 Hep3Vector u_Boost = -BoostP.boostVector();
413
414 if ( npip == 1 )
415 {
416 //
417 //---- Find the Best Pi0 -----
418 //
419 HepLorentzVector p2g, ppi0;
420 double delmpi = 999.;
421 int ig11 = -1, ig21 = -1;
422
423 for ( int i = 0; i < nGam - 1; i++ )
424 {
425 for ( int j = i + 1; j < nGam; j++ )
426 {
427 HepLorentzVector ppi0 = pGam[i] + pGam[j];
428 if ( fabs( ppi0.m() - mpi0c ) > delmpi ) continue;
429 if ( ppi0.m() > 0.2 ) continue;
430 delmpi = fabs( ppi0.m() - mpi0c );
431 ig11 = i;
432 ig21 = j;
433 }
434 }
435 if ( ig11 == -1 || ig21 == -1 ) return sc;
436 p2g = pGam[ig11] + pGam[ig21];
437 m_m2graw = p2g.m();
438
439 HepLorentzVector Ppip1 = ppip[0];
440 Ppip1.boost( u_Boost );
441 p2g.boost( u_Boost );
442 //
443 // ----- Pion Efficiency -----
444 //
445 //
446 //------ Missing pi-
447 //
448 HepLorentzVector Pmiss = -Ppip1 - p2g;
449 m_pmiss = Pmiss.rho();
450 double emiss = ecms - Ppip1.e() - p2g.e();
451 m_emiss = emiss;
452 m_mmiss = sqrt( fabs( emiss * emiss - Pmiss.rho() * Pmiss.rho() ) );
453
454 HepLorentzVector Prho0, Prhop, Prhom, pmisspi, ptot;
455 pmisspi.setPx( Pmiss.px() );
456 pmisspi.setPy( Pmiss.py() );
457 pmisspi.setPz( Pmiss.pz() );
458 pmisspi.setE( emiss );
459
460 Prho0 = Ppip1 + pmisspi;
461 Prhop = Ppip1 + p2g;
462 Prhom = pmisspi + p2g;
463 ptot = Ppip1 + pmisspi + p2g;
464
465 m_pmrho0 = Prho0.rho();
466 m_mmrho0 = Prho0.m();
467 // m_prho0 = (Ppip1 + Ppim1).rho();
468 // m_mrho0 = (Ppip1 + Ppim1).m();
469 m_prhop = Prhop.rho();
470 m_mrhop = Prhop.m();
471 // m_prhom = (p2g + Ppim1).rho();
472 // m_mrhom = (p2g + Ppim1).m();
473 m_pmrhom = Prhom.rho();
474 m_mmrhom = Prhom.m();
475 m_ppipraw = Ppip1.rho();
476
477 m_tuple5->write();
478 }
479
480 if ( npip * npim != 1 ) return sc;
481 m_pass[3] += 1;
482
483 HepLorentzVector Ppim1 = ppim[0];
484 Ppim1.boost( u_Boost );
485 HepLorentzVector Ppip1 = ppip[0];
486 Ppip1.boost( u_Boost );
487
488 HepLorentzVector Pneumiss = -( Ppim1 + Ppip1 );
489 m_pneumiss = Pneumiss.rho();
490 double eneumiss = ecms - Ppim1.e() - Ppip1.e();
491 m_eneumiss = eneumiss;
492 m_mneumiss = sqrt( fabs( eneumiss * eneumiss - Pneumiss.rho() * Pneumiss.rho() ) );
493
494 //
495 // Let's play vertex fit here
496 //
497 HepPoint3D vx( 0., 0., 0. );
498 HepSymMatrix Evx( 3, 0 );
499 double bx = 1E+6;
500 double by = 1E+6;
501 double bz = 1E+6;
502 Evx[0][0] = bx * bx;
503 Evx[1][1] = by * by;
504 Evx[2][2] = bz * bz;
505
506 VertexParameter vxpar;
507 vxpar.setVx( vx );
508 vxpar.setEvx( Evx );
509
510 VertexFit* vtxfit = VertexFit::instance();
511 KinematicFit* kmfit = KinematicFit::instance();
512
513 RecMdcKalTrack* pipTrk = ( *( evtRecTrkCol->begin() + ipip[0] ) )->mdcKalTrack();
514 pipTrk->setPidType( RecMdcKalTrack::pion );
515 WTrackParameter wpip0( mpi, pipTrk->helix(), pipTrk->err() );
516
517 pipTrk->setPidType( RecMdcKalTrack::kaon );
518 WTrackParameter wkp0( mk, pipTrk->helix(), pipTrk->err() );
519
520 RecMdcKalTrack* pimTrk = ( *( evtRecTrkCol->begin() + ipim[0] ) )->mdcKalTrack();
521 pimTrk->setPidType( RecMdcKalTrack::pion );
522 WTrackParameter wpim0( mpi, pimTrk->helix(), pimTrk->err() );
523
524 pimTrk->setPidType( RecMdcKalTrack::kaon );
525 WTrackParameter wkm0( mk, pimTrk->helix(), pimTrk->err() );
526
527 vtxfit->init();
528 vtxfit->AddTrack( 0, wpip0 );
529 vtxfit->AddTrack( 1, wpim0 );
530 vtxfit->AddVertex( 0, vxpar, 0, 1 );
531 if ( !vtxfit->Fit( 0 ) ) return sc;
532
533 m_pass[4] += 1;
534
535 WTrackParameter wpip = vtxfit->wtrk( 0 );
536 WTrackParameter wpim = vtxfit->wtrk( 1 );
537
538 //
539 // Apply Kinematic 4C fit
540 //
541
542 double chisq4c = 9999.;
543
544 int ig1 = -1;
545 int ig2 = -1;
546
547 for ( int i = 0; i < nGam - 1; i++ )
548 {
549 RecEmcShower* g1Trk = ( *( evtRecTrkCol->begin() + iGam[i] ) )->emcShower();
550 for ( int j = i + 1; j < nGam; j++ )
551 {
552 RecEmcShower* g2Trk = ( *( evtRecTrkCol->begin() + iGam[j] ) )->emcShower();
553
554 kmfit->init();
555 kmfit->AddTrack( 0, wpip );
556 kmfit->AddTrack( 1, wpim );
557 kmfit->AddTrack( 2, 0.0, g1Trk );
558 kmfit->AddTrack( 3, 0.0, g2Trk );
559 kmfit->AddFourMomentum( 0, BoostP );
560
561 if ( !kmfit->Fit() ) continue;
562 double chi2 = kmfit->chisq();
563 if ( chi2 > chisq4c ) continue;
564 chisq4c = chi2;
565 ig1 = i;
566 ig2 = j;
567 } // gamma1
568 } // gamma2
569
570 if ( chisq4c > 500 || ig1 < 0 ) return sc;
571 RecEmcShower* g1Trk = ( *( evtRecTrkCol->begin() + iGam[ig1] ) )->emcShower();
572 RecEmcShower* g2Trk = ( *( evtRecTrkCol->begin() + iGam[ig2] ) )->emcShower();
573
574 kmfit->init();
575 kmfit->AddTrack( 0, wpip );
576 kmfit->AddTrack( 1, wpim );
577 kmfit->AddTrack( 2, 0.0, g1Trk );
578 kmfit->AddTrack( 3, 0.0, g2Trk );
579 kmfit->AddFourMomentum( 0, BoostP );
580 if ( !kmfit->Fit() ) return sc;
581 double chi = kmfit->chisq();
582 m_chisq4c = chi;
583 m_chisq_4c->Fill( chi );
584
585 HepLorentzVector p2gam = kmfit->pfit( 2 ) + kmfit->pfit( 3 );
586 m_ppi0 = p2gam.rho();
587 m_mpi0 = p2gam.m();
588
589 m_pi0_mass->Fill( p2gam.m() );
590
591 m_pass[5] += 1;
592
593 //
594 // Apply Kinematic 5C fit
595 //
596 //
597 // 5c for J/psi-->pi+pi-pi0
598 //
599 kmfit->init();
600 kmfit->AddTrack( 0, wpip );
601 kmfit->AddTrack( 1, wpim );
602 kmfit->AddTrack( 2, 0.0, g1Trk );
603 kmfit->AddTrack( 3, 0.0, g2Trk );
604 kmfit->AddResonance( 0, mpi0c, 2, 3 );
605 kmfit->AddFourMomentum( 1, BoostP );
606 if ( !kmfit->Fit() ) return sc;
607 double chis = kmfit->chisq();
608
609 m_pass[6] += 1;
610 //
611 // 5c for J/psi-->K+K-pi0
612 //
613 double chisk5c = 9999.;
614
615 vtxfit->init();
616 vtxfit->AddTrack( 0, wkp0 );
617 vtxfit->AddTrack( 1, wkm0 );
618 vtxfit->AddVertex( 0, vxpar, 0, 1 );
619
620 if ( vtxfit->Fit( 0 ) )
621 {
622
623 WTrackParameter wkp = vtxfit->wtrk( 0 );
624 WTrackParameter wkm = vtxfit->wtrk( 1 );
625
626 kmfit->init();
627 kmfit->AddTrack( 0, wkp );
628 kmfit->AddTrack( 1, wkm );
629 kmfit->AddTrack( 2, 0.0, g1Trk );
630 kmfit->AddTrack( 3, 0.0, g2Trk );
631 kmfit->AddResonance( 0, mpi0c, 2, 3 );
632 kmfit->AddFourMomentum( 1, BoostP );
633 if ( kmfit->Fit() ) chisk5c = kmfit->chisq();
634 }
635 if ( chisk5c < chis ) return sc;
636
637 m_pass[7] += 1;
638 //
639 // 5c Fit again for J/psi-->pi+pi-pi0
640 //
641 kmfit->init();
642 kmfit->AddTrack( 0, wpip );
643 kmfit->AddTrack( 1, wpim );
644 kmfit->AddTrack( 2, 0.0, g1Trk );
645 kmfit->AddTrack( 3, 0.0, g2Trk );
646 kmfit->AddResonance( 0, mpi0c, 2, 3 );
647 kmfit->AddFourMomentum( 1, BoostP );
648 if ( !kmfit->Fit( 0 ) ) return sc;
649 if ( !kmfit->Fit() ) return sc;
650
651 m_pass[8] += 1;
652 m_chisq5c = kmfit->chisq();
653
654 m_chisq_5c->Fill( kmfit->chisq() );
655
656 HepLorentzVector ppi1 = kmfit->pfit( 0 );
657 HepLorentzVector ppi2 = kmfit->pfit( 1 );
658 HepLorentzVector pgam1 = kmfit->pfit( 2 );
659 HepLorentzVector pgam2 = kmfit->pfit( 3 );
660 HepLorentzVector p2gam1 = kmfit->pfit( 2 ) + kmfit->pfit( 3 );
661 HepLorentzVector p2pi = kmfit->pfit( 0 ) + kmfit->pfit( 1 );
662 m_ppi0fit = p2gam1.rho();
663 m_mpi0fit = p2gam1.m();
664 m_ppip = ppi1.rho();
665 m_ppim = ppi2.rho();
666 m_p2pi = p2pi.rho();
667 m_m2pi = p2pi.m();
668 m_ppip0 = ( ppi1 + p2gam1 ).rho();
669 m_mpip0 = ( ppi1 + p2gam1 ).m();
670 m_ppim0 = ( ppi2 + p2gam1 ).rho();
671 m_mpim0 = ( ppi2 + p2gam1 ).m();
672
673 m_pip_mom->Fill( ppi1.rho() );
674 m_pim_mom->Fill( ppi2.rho() );
675 m_rho0_mass->Fill( p2pi.m() );
676 m_rhop_mass->Fill( ( ppi1 + p2gam1 ).m() );
677 m_rhom_mass->Fill( ( ppi2 + p2gam1 ).m() );
678
679 //--------------Angle between two charged pions ---------
680 double theta2pi = ppi1.px() * ppi2.px() + ppi1.py() * ppi2.py() + ppi1.pz() * ppi2.pz();
681 double time2pi = ppi1.rho() * ppi2.rho();
682 if ( time2pi == 0. ) return sc;
683
684 m_pass[9] += 1;
685
686 theta2pi = theta2pi / time2pi;
687
688 // cout <<"time2pi=====" << time2pi <<endl;
689
690 if ( theta2pi <= -1. ) theta2pi = 180.;
691 else if ( theta2pi >= 1. ) theta2pi = 0.;
692 else if ( fabs( theta2pi ) < 1. ) theta2pi = acos( theta2pi ) * 180. / 3.14159;
693 m_theta2pi = theta2pi;
694
695 //----- Theta angle of gamma in pi0 system --------
696 HepLorentzVector pg1inpi0, pg2inpi0;
697
698 // HepLorentzVector betapi0 = p2gam1/(p2gam1.e());
699 double betapi0 = p2gam1.beta();
700 // ==================================
701 Hep3Vector twogam( p2gam.px(), p2gam.py(), p2gam.pz() );
702 // CLHEP::boostOf(pg1inpi0,gam1,betapi0);
703 m_g1inpi0the =
704 ( ( boostOf( pgam1, twogam / p2gam.rho(), betapi0 ) ).theta() ) * 180. / 3.14159;
705 m_g2inpi0the =
706 ( ( boostOf( pgam2, twogam / p2gam.rho(), betapi0 ) ).theta() ) * 180. / 3.14159;
707 /*
708 //
709 // ----- Pion Efficiency -----
710 //
711 HepLorentzVector p2g = pGam[ig1] + pGam[ig2];
712 m_m2graw = p2g.m();
713 p2g.boost(u_Boost);
714 //
715 //------ Missing pi-
716 //
717 HepLorentzVector Pmiss = -Ppip1 - p2g;
718 m_pmiss = Pmiss.rho();
719 double emiss = ecms - Ppip1.e() - p2g.e();
720 m_emiss = emiss;
721 m_mmiss = sqrt(fabs(emiss*emiss - Pmiss.rho()*Pmiss.rho()));
722
723 HepLorentzVector Prho0, Prhop, Prhom, pmisspi, ptot;
724 pmisspi.setPx(Pmiss.px());
725 pmisspi.setPy(Pmiss.py());
726 pmisspi.setPz(Pmiss.pz());
727 pmisspi.setE(emiss);
728
729 Prho0 = Ppip1 + pmisspi;
730 Prhop = Ppip1 + p2g;
731 Prhom = pmisspi + p2g;
732 ptot = Ppip1 + pmisspi + p2g;
733
734 m_pmrho0 = Prho0.rho();
735 m_mmrho0 = Prho0.m();
736 m_prho0 = (Ppip1 + Ppim1).rho();
737 m_mrho0 = (Ppip1 + Ppim1).m();
738 m_prhop = Prhop.rho();
739 m_mrhop = Prhop.m();
740 m_prhom = (p2g + Ppim1).rho();
741 m_mrhom = (p2g + Ppim1).m();
742 m_pmrhom = Prhom.rho();
743 m_mmrhom = Prhom.m();
744 */
745 m_tuple4->write();
746
747 return StatusCode::SUCCESS;
748}
749
750// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
751StatusCode ValidJpsiRhopi::finalize() {
752
753 MsgStream log( msgSvc(), name() );
754 log << MSG::ALWAYS << "Total Entries : " << m_pass[0] << endmsg;
755 log << MSG::ALWAYS << "Ncharge=2 Cut : " << m_pass[1] << endmsg;
756 log << MSG::ALWAYS << "Ngam<2 cut : " << m_pass[2] << endmsg;
757 log << MSG::ALWAYS << "Nch+=1,Nch-=1 cut : " << m_pass[3] << endmsg;
758 log << MSG::ALWAYS << "Vertex Fit : " << m_pass[4] << endmsg;
759 log << MSG::ALWAYS << "4c-Fit : " << m_pass[5] << endmsg;
760 log << MSG::ALWAYS << "5c-Fit : " << m_pass[6] << endmsg;
761 log << MSG::ALWAYS << "chi3pi<chikkpi cut : " << m_pass[7] << endmsg;
762 log << MSG::ALWAYS << "5c-Fit Again: " << m_pass[8] << endmsg;
763 log << MSG::ALWAYS << "Final : " << m_pass[9] << endmsg;
764 return StatusCode::SUCCESS;
765}
DECLARE_COMPONENT
DECLARE_COMPONENT(BesBdkRc)
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition CosmicGenerator.h:36
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
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
ntupleSvc
INTupleSvc * ntupleSvc()
Definition ServiceAccessor.h:51
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:12
mpi0c
const double mpi0c
Definition ValidJpsiRhopi.cxx:39
KinematicFit::instance
static KinematicFit * instance()
Definition KinematicFit.cxx:18
KinematicFit::AddResonance
void AddResonance(int number, double mres, std::vector< int > tlis)
Definition KinematicFit.cxx:137
KinematicFit::AddFourMomentum
void AddFourMomentum(int number, HepLorentzVector p4)
Definition KinematicFit.cxx:341
ParticleID::instance
static ParticleID * instance()
Definition ParticleID.cxx:21
ParticleID::IsPidInfoValid
bool IsPidInfoValid() const
ParticleID::calculate
void calculate()
Definition ParticleID.cxx:97
ParticleID::init
void init()
Definition ParticleID.cxx:26
TrackPool::AddTrack
void AddTrack(const int number, const double mass, const RecMdcTrack *trk)
Definition TrackPool.cxx:21
ValidJpsiRhopi::finalize
StatusCode finalize()
Definition ValidJpsiRhopi.cxx:751
ValidJpsiRhopi::execute
StatusCode execute()
Definition ValidJpsiRhopi.cxx:205
ValidJpsiRhopi::initialize
StatusCode initialize()
Definition ValidJpsiRhopi.cxx:74
ValidJpsiRhopi::ValidJpsiRhopi
ValidJpsiRhopi(const std::string &name, ISvcLocator *pSvcLocator)
Definition ValidJpsiRhopi.cxx:57
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