BOSS 7.1.0
BESIII Offline Software System
Loading...
Searching...
No Matches
Pi0EtaToGGRecAlg.cxx
Go to the documentation of this file.
1#include "GaudiKernel/SmartIF.h"
2#include "GaudiKernel/MsgStream.h"
3#include "GaudiKernel/SmartDataPtr.h"
4#include "GaudiKernel/IDataManagerSvc.h"
5
12
14
15const double xmpi0 = 0.134976; // BOSS 6.4.1 pdt.table value
16const double xmeta = 0.54784; // PDG08 = 0.547853, BOSS 6.4.1 pdt.table = 0.54784
17DECLARE_COMPONENT(Pi0EtaToGGRecAlg)
18//*************************************************************************************
19Pi0EtaToGGRecAlg::Pi0EtaToGGRecAlg(const std::string& name, ISvcLocator* pSvcLocator) :
20 Algorithm(name, pSvcLocator)
21{
22 //Declare the properties
23
24 declareProperty("RejectBothInEndcap", m_rejectEndEnd = false);
25
26 /// pi0 properties
27 declareProperty("Pi0MinMass", m_pi0minmass_cut = 0.08);
28 declareProperty("Pi0MaxMass", m_pi0maxmass_cut = 0.18);
29 declareProperty("Pi0ChisqCut", m_pi0chisq_cut = 2500);
30
31 /// eta properties
32 declareProperty("EtaMinMass", m_etaminmass_cut = 0.4);
33 declareProperty("EtaMaxMass", m_etamaxmass_cut = 0.7);
34 declareProperty("EtaChisqCut", m_etachisq_cut = 2500);
35
36
37 /// Individual photons properties
38 declareProperty("PhotonMinEnergy", m_minEnergy = 0.025);
39
40 declareProperty("PhotonInBarrelOrEndcap", m_useBarrelEndcap = true);
41 declareProperty("PhotonMaxCosThetaBarrel", m_maxCosThetaBarrel = 0.8);
42 declareProperty("PhotonMinCosThetaEndcap", m_minCosThetaEndcap = 0.86);
43 declareProperty("PhotonMaxCosThetaEndcap", m_maxCosThetaEndcap = 0.92);
44 declareProperty("PhotonMinEndcapEnergy", m_minEndcapEnergy = 0.050);
45
46 declareProperty("PhotonApplyTimeCut", m_applyTimeCut = true);
47 declareProperty("PhotonMinTime", m_minTime = 0.);
48 declareProperty("PhotonMaxTime", m_maxTime = 14.);
49 declareProperty("PhotonDeltaTime", m_deltaTime = 10.);
50
51 declareProperty("PhotonApplyDangCut", m_applyDangCut = false);
52 declareProperty("PhotonMinDang", m_minDang = 20.0);
53
54 /// Limits for number of photons, pi0s and etas
55 declareProperty("MaxNGoodPhoton", m_maxNGoodPhoton = 50);
56 declareProperty("MaxNPi0", m_maxNPi0 = 100);
57 //declareProperty("MaxNEta", m_maxNEta = 100);
58
59 declareProperty("DoClear", m_doClear = false);
60}
61
62//******************************************************************************************
64
65 MsgStream log(msgSvc(), name());
66 log << MSG::INFO << "in initialize()" <<endreq;
67
68 return StatusCode::SUCCESS;
69}
70
71
72//*********************************************************************************************
73StatusCode Pi0EtaToGGRecAlg::execute() {
74
75 MsgStream log(msgSvc(), name());
76 log << MSG::INFO << "in execute()" <<endreq;
77
78 // get event header, eventlist and tracklist from TDS
79 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
80 SmartDataPtr<EvtRecEvent> recEvt(eventSvc(), EventModel::EvtRec::EvtRecEvent);
81 SmartDataPtr<EvtRecTrackCol> recTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
82
83 SmartDataPtr<EvtRecPi0Col> recPi0Col(eventSvc(), EventModel::EvtRec::EvtRecPi0Col);
84 if ( !recPi0Col ) {
85 recPi0Col = new EvtRecPi0Col;
86 StatusCode sc = eventSvc()->registerObject(EventModel::EvtRec::EvtRecPi0Col, recPi0Col);
87 if ( sc.isFailure() ) {
88 log << MSG::ERROR << "could not register EvtRecPi0Col in TDS" <<endreq;
89 return StatusCode::FAILURE;
90 }
91 }
92
93 SmartDataPtr<EvtRecEtaToGGCol> recEtaToGGCol(eventSvc(), EventModel::EvtRec::EvtRecEtaToGGCol);
94 if ( !recEtaToGGCol ) {
95 recEtaToGGCol = new EvtRecEtaToGGCol;
96 StatusCode sc = eventSvc()->registerObject(EventModel::EvtRec::EvtRecEtaToGGCol, recEtaToGGCol);
97 if ( sc.isFailure() ) {
98 log << MSG::ERROR << "could not register EvtRecEtaToGGCol in TDS" <<endreq;
99 return StatusCode::FAILURE;
100 }
101 }
102
103 if ( m_doClear ) {
104 SmartIF<IDataManagerSvc> dataManSvc(eventSvc());
105 dataManSvc->clearSubTree(EventModel::EvtRec::EvtRecPi0Col);
106 dataManSvc->clearSubTree(EventModel::EvtRec::EvtRecEtaToGGCol);
107 }
108
109 // select good photons
110 std::vector<EvtRecTrack*> vGoodPhotons;
111 for ( int i = recEvt->totalCharged(); i < recEvt->totalTracks(); i++ ) {
112 EvtRecTrack* gTrk = *(recTrkCol->begin() + i);
113 if ( !validPhoton(gTrk,recEvt,recTrkCol) ) continue;
114 vGoodPhotons.push_back(gTrk);
115 }
116
117 int nGoodPhoton = vGoodPhotons.size();
118 if ( nGoodPhoton > m_maxNGoodPhoton ) return StatusCode::SUCCESS;
119
120 /// Needed for 1C mass fits for pi0 and eta
122
123 //
124 // Form two-photon combinations
125 //
126 for(int i1 = 0; i1 < (nGoodPhoton-1); i1++) {
127 EvtRecTrack* g1Trk = vGoodPhotons[i1];
128 RecEmcShower* g1Shower = g1Trk->emcShower();
129 HepLorentzVector g1P4 = getP4(g1Shower);
130
131 for(int i2 = i1+1; i2 < nGoodPhoton; i2++) {
132 EvtRecTrack* g2Trk = vGoodPhotons[i2];
133 RecEmcShower* g2Shower = g2Trk->emcShower();
134 HepLorentzVector g2P4 = getP4(g2Shower);
135
136 HepLorentzVector p2g = g1P4 + g2P4;
137
138 /// If RejectBothInEndcap=true, reject candidate if both photons are in endcaps
139 if(m_rejectEndEnd&&bothInEndcap(g1P4,g2P4)) continue;
140
141 /// Select pi0 candidates
142 if ((p2g.m() > m_pi0minmass_cut) && (p2g.m() < m_pi0maxmass_cut ))
143 {
144
145 /// Perform pi0 1C KinematicFit
146 kmfit->init();
147 kmfit->setIterNumber(5);
148 kmfit->AddTrack(0, 0.0, g1Shower);
149 kmfit->AddTrack(1, 0.0, g2Shower);
150 kmfit->AddResonance(0, xmpi0, 0, 1);
151
152 kmfit->Fit(0); // Perform fit
153 kmfit->BuildVirtualParticle(0);
154
155 double pi0_chisq = kmfit->chisq(0);
156 if ( pi0_chisq >= m_pi0chisq_cut ) continue;
157
158
159 /// Fill EvtRecPi0
160 EvtRecPi0* recPi0 = new EvtRecPi0();
161 recPi0->setUnconMass(p2g.restMass());
162 recPi0->setChisq(pi0_chisq);
163
164 if ( g1P4.e() >= g2P4.e() ) {
165 recPi0->setHiPfit(kmfit->pfit(0));
166 recPi0->setLoPfit(kmfit->pfit(1));
167 recPi0->setHiEnGamma(g1Trk);
168 recPi0->setLoEnGamma(g2Trk);
169 }
170 else {
171 recPi0->setHiPfit(kmfit->pfit(1));
172 recPi0->setLoPfit(kmfit->pfit(0));
173 recPi0->setHiEnGamma(g2Trk);
174 recPi0->setLoEnGamma(g1Trk);
175 }
176 recPi0Col->push_back(recPi0);
177
178 } // End of pi0 mass window IF
179
180
181 /// Select eta candidates
182 if ((p2g.m() > m_etaminmass_cut) && (p2g.m() < m_etamaxmass_cut ))
183 {
184
185 /// Perform eta 1C KinematicFit
186 kmfit->init();
187 kmfit->setIterNumber(5);
188 kmfit->AddTrack(0, 0.0, g1Shower);
189 kmfit->AddTrack(1, 0.0, g2Shower);
190 kmfit->AddResonance(0, xmeta, 0, 1);
191
192 kmfit->Fit(0); // Perform fit
193 kmfit->BuildVirtualParticle(0);
194
195 double eta_chisq = kmfit->chisq(0);
196 if ( eta_chisq >= m_etachisq_cut ) continue;
197
198
199 /// Fill EvtRecEtaToGG
200 EvtRecEtaToGG* recEtaToGG = new EvtRecEtaToGG();
201 recEtaToGG->setUnconMass(p2g.restMass());
202 recEtaToGG->setChisq(eta_chisq);
203
204 if ( g1P4.e() >= g2P4.e() ) {
205 recEtaToGG->setHiPfit(kmfit->pfit(0));
206 recEtaToGG->setLoPfit(kmfit->pfit(1));
207 recEtaToGG->setHiEnGamma(g1Trk);
208 recEtaToGG->setLoEnGamma(g2Trk);
209 }
210 else {
211 recEtaToGG->setHiPfit(kmfit->pfit(1));
212 recEtaToGG->setLoPfit(kmfit->pfit(0));
213 recEtaToGG->setHiEnGamma(g2Trk);
214 recEtaToGG->setLoEnGamma(g1Trk);
215 }
216 recEtaToGGCol->push_back(recEtaToGG);
217
218 } // End of etaToGG mass window IF
219
220 } // End of g2Trk evtRec FOR
221 } // End of g1Trk evtRec FOR
222
223 // If there are too many pi0s, it's possibly a bad event and all pi0s are abandoned
224 if ( recPi0Col->size() > m_maxNPi0 ) {
225 recPi0Col->clear();
226 }
227
228 return StatusCode::SUCCESS;
229}
230
231
232//********************************************************************************************
233StatusCode Pi0EtaToGGRecAlg::finalize() {
234
235 MsgStream log(msgSvc(), name());
236 log << MSG::INFO << "in finalize()" << endreq;
237
238 return StatusCode::SUCCESS;
239}
240
241
242
243/// ************************** FUNCTIONS *****************************************************
244HepLorentzVector Pi0EtaToGGRecAlg::getP4(RecEmcShower* gTrk){
245
246 double eraw = gTrk->energy();
247 double phi = gTrk->phi();
248 double the = gTrk->theta();
249
250 return HepLorentzVector( eraw * sin(the) * cos(phi),
251 eraw * sin(the) * sin(phi),
252 eraw * cos(the),
253 eraw );
254}
255
256
257bool Pi0EtaToGGRecAlg::validPhoton(EvtRecTrack* recTrk,
258 SmartDataPtr<EvtRecEvent> recEvt,
259 SmartDataPtr<EvtRecTrackCol> recTrkCol)
260{
261 if ( !recTrk->isEmcShowerValid() ) return false;
262
263 RecEmcShower *emcTrk = recTrk->emcShower();
264
265 HepLorentzVector shP4 = getP4(emcTrk);
266 double cosThetaSh = shP4.vect().cosTheta();
267
268
269 /// Minimum energy
270 if (shP4.e() <= m_minEnergy) return false;
271
272
273 /// Barrel/Endcap
274 if ( m_useBarrelEndcap ) {
275 bool inBarrelEndcap = false;
276
277 if(fabs(cosThetaSh) < m_maxCosThetaBarrel) inBarrelEndcap = true;
278
279 if((fabs(cosThetaSh) > m_minCosThetaEndcap)
280 &&(fabs(cosThetaSh) < m_maxCosThetaEndcap)
281 &&(shP4.e() > m_minEndcapEnergy)) inBarrelEndcap = true;
282
283 if ( !inBarrelEndcap ) return false;
284 }
285
286
287 /// Time, only apply timing cuts if "recEvt->totalCharged() > 0"
288 if ( m_applyTimeCut ) {
289 double time = emcTrk->time();
290 if ( recEvt->totalCharged() > 0 ) {
291 if ( time < m_minTime || time > m_maxTime ) return false;
292 }
293 else {
294 RecEmcShower* firstG = (*(recTrkCol->begin()))->emcShower();
295 double deltaTime = fabs(time - firstG->time());
296 if ( deltaTime > 10 ) return false;
297 }
298 }
299
300
301 /// Dang
302 if (m_applyDangCut && recEvt->totalCharged() > 0)
303 {
304 Hep3Vector emcpos(emcTrk->x(), emcTrk->y(), emcTrk->z());
305 double dang = 200.;
306
307 for (int j = 0; j < recEvt->totalCharged(); j++) {
308 EvtRecTrackIterator jtTrk = recTrkCol->begin() + j;
309 if ( !(*jtTrk)->isExtTrackValid() ) continue;
310 RecExtTrack* extTrk = (*jtTrk)->extTrack();
311 if ( extTrk->emcVolumeNumber() == -1 ) continue;
312 Hep3Vector extpos = extTrk->emcPosition();
313 double angd1 = extpos.angle(emcpos);
314 if ( angd1 < dang ) dang = angd1;
315 }
316
317 if ( dang < 200 ) {
318 dang = dang * 180 / (CLHEP::pi);
319 if ( dang <= m_minDang ) return false;
320 }
321 } // End of "recEvt->totalCharged() > 0" IF
322
323 return true;
324}
325
326
327bool Pi0EtaToGGRecAlg::bothInEndcap(HepLorentzVector g1_P4, HepLorentzVector g2_P4){
328
329 double g1_CosTheta = g1_P4.vect().cosTheta();
330 double g2_CosTheta = g2_P4.vect().cosTheta();
331
332 bool g1InEndcap = false;
333 bool g2InEndcap = false;
334
335 if((fabs(g1_CosTheta) > m_minCosThetaEndcap)
336 &&(fabs(g1_CosTheta) < m_maxCosThetaEndcap)) g1InEndcap = true;
337
338 if((fabs(g2_CosTheta) > m_minCosThetaEndcap)
339 &&(fabs(g2_CosTheta) < m_maxCosThetaEndcap)) g2InEndcap = true;
340
341 if(g1InEndcap&&g2InEndcap) return( true );
342
343 return( false );
344}
double sin(const BesAngle a)
Definition: BesAngle.h:210
double cos(const BesAngle a)
Definition: BesAngle.h:213
Double_t time
ObjectVector< EvtRecEtaToGG > EvtRecEtaToGGCol
Definition: EvtRecEtaToGG.h:58
ObjectVector< EvtRecPi0 > EvtRecPi0Col
Definition: EvtRecPi0.h:58
EvtRecTrackCol::iterator EvtRecTrackIterator
Definition: EvtRecTrack.h:111
const double xmeta
const double xmpi0
***************************************************************************************Pseudo Class RRes *****************************************************************************************Parameters and physical constants **Maarten sept ************************************************************************DOUBLE PRECISION xsmu **************************************************************************PARTICLE DATA all others are from PDG *Only resonances with known widths into electron pairs are sept ************************************************************************C Declarations xmb DOUBLE PRECISION xcb DOUBLE PRECISION xmpi0
Definition: RRes.h:32
IMessageSvc * msgSvc()
double theta() const
Definition: DstEmcShower.h:38
double phi() const
Definition: DstEmcShower.h:39
double x() const
Definition: DstEmcShower.h:35
double time() const
Definition: DstEmcShower.h:50
double z() const
Definition: DstEmcShower.h:37
double energy() const
Definition: DstEmcShower.h:45
double y() const
Definition: DstEmcShower.h:36
const Hep3Vector emcPosition() const
Definition: DstExtTrack.h:126
const int emcVolumeNumber() const
Definition: DstExtTrack.h:132
void setLoEnGamma(const EvtRecTrack *trk)
Definition: EvtRecEtaToGG.h:41
void setHiEnGamma(const EvtRecTrack *trk)
Definition: EvtRecEtaToGG.h:40
void setUnconMass(const double unconMass)
Definition: EvtRecEtaToGG.h:34
void setChisq(const double chisq)
Definition: EvtRecEtaToGG.h:35
void setHiPfit(const HepLorentzVector &hiPfit)
Definition: EvtRecEtaToGG.h:37
void setLoPfit(const HepLorentzVector &loPfit)
Definition: EvtRecEtaToGG.h:38
void setChisq(const double chisq)
Definition: EvtRecPi0.h:35
void setLoPfit(const HepLorentzVector &loPfit)
Definition: EvtRecPi0.h:38
void setLoEnGamma(const EvtRecTrack *trk)
Definition: EvtRecPi0.h:41
void setHiPfit(const HepLorentzVector &hiPfit)
Definition: EvtRecPi0.h:37
void setUnconMass(const double unconMass)
Definition: EvtRecPi0.h:34
void setHiEnGamma(const EvtRecTrack *trk)
Definition: EvtRecPi0.h:40
RecEmcShower * emcShower()
Definition: EvtRecTrack.h:58
bool isEmcShowerValid()
Definition: EvtRecTrack.h:47
void setIterNumber(const int niter=5)
HepLorentzVector pfit(int n) const
void BuildVirtualParticle(int number)
void AddResonance(int number, double mres, std::vector< int > tlis)
static KalmanKinematicFit * instance()
Pi0EtaToGGRecAlg(const std::string &name, ISvcLocator *pSvcLocator)
void AddTrack(const int number, const double mass, const RecMdcTrack *trk)
Definition: TrackPool.cxx:22
_EXTERN_ std::string EvtRecPi0Col
Definition: EventModel.h:123
_EXTERN_ std::string EvtRecEvent
Definition: EventModel.h:116
_EXTERN_ std::string EvtRecEtaToGGCol
Definition: EventModel.h:124
_EXTERN_ std::string EvtRecTrackCol
Definition: EventModel.h:117