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