BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
DiGam.cxx
Go to the documentation of this file.
1#include "CLHEP/Geometry/Point3D.h"
2#include "CLHEP/Vector/LorentzVector.h"
3#include "CLHEP/Vector/ThreeVector.h"
4#include "CLHEP/Vector/TwoVector.h"
5#include "GaudiKernel/Bootstrap.h"
6#include "GaudiKernel/IDataProviderSvc.h"
7#include "GaudiKernel/ISvcLocator.h"
8#include "GaudiKernel/ITHistSvc.h"
9#include "GaudiKernel/MsgStream.h"
10#include "GaudiKernel/NTuple.h"
11#include "GaudiKernel/SmartDataPtr.h"
12#include "TMath.h"
13
14#include "DstEvent/TofHitStatus.h"
15#include "EventModel/Event.h"
16#include "EventModel/EventHeader.h"
17#include "EventModel/EventModel.h"
18#include "EvtRecEvent/EvtRecEvent.h"
19#include "EvtRecEvent/EvtRecTrack.h"
20#include "McTruth/McParticle.h"
21#include "ParticleID/ParticleID.h"
22#include "VertexFit/KinematicFit.h"
23#include "VertexFit/VertexFit.h"
24
25#include "DiGam.h"
26#include "Parameter.h"
27#include <vector>
28
29const double xmass[5] = { 0.000511, 0.105658, 0.139570, 0.493677, 0.938272 };
30const double velc = 299.792458;
31typedef std::vector<int> Vint;
32typedef std::vector<HepLorentzVector> Vp4;
33
34const double pai = 3.1415926;
35int N0;
36int N1;
37int N2;
38int N3;
39int N4;
40DECLARE_COMPONENT( DiGam )
41DiGam::DiGam( const std::string& name, ISvcLocator* pSvcLocator )
42 : Algorithm( name, pSvcLocator ) {
43 declareProperty( "jpsiCrossSection", jpsiCrossSection = 249.7 );
44 declareProperty( "jpsiMCEff", jpsiMCEff = 0.07145 );
45 declareProperty( "jpsiMCEffBoost", jpsiMCEffBoost = 0.07099 );
46
47 declareProperty( "psi2sCrossSection", psi2sCrossSection = 180.8 );
48 declareProperty( "psi2sMCEff", psi2sMCEff = 0.07159 );
49 declareProperty( "psi2sMCEffBoost", psi2sMCEffBoost = 0.07123 );
50
51 declareProperty( "psi3770CrossSection", psi3770CrossSection = 173.4 );
52 declareProperty( "psi3770MCEff", psi3770MCEff = 0.071725 );
53 declareProperty( "psi3770MCEffBoost", psi3770MCEffBoost = 0.0713125 );
54
55 declareProperty( "CrossSection", CrossSection );
56 declareProperty( "MCEff", MCEff );
57 declareProperty( "MCEffBoost", MCEffBoost );
58
59 declareProperty( "boostPhimin", boostPhimin = 2.5 );
60 declareProperty( "boostPhimax", boostPhimax = 5 );
61 declareProperty( "boostMinEmin", boostMinEmin );
62 declareProperty( "boostMinEmax", boostMinEmax );
63
64 declareProperty( "RunModel", RunModel = 2000 );
65 declareProperty( "dPhiMin", dPhiMin = -7 );
66 declareProperty( "dPhiMax", dPhiMax = 5 );
67 declareProperty( "dPhiMinSig", dPhiMinSig = -4 );
68 declareProperty( "dPhiMaxSig", dPhiMaxSig = 2 );
69
70 declareProperty( "flag", flag = true );
71 declareProperty( "E_cms", E_cms );
72}
73//***************************************************************************************
74StatusCode DiGam::initialize() {
75 MsgStream log( msgSvc(), name() );
76 log << MSG::INFO << "in initialize()" << endmsg;
77
78 if ( RunModel == 1 )
79 { // jpsi
80 CrossSection = jpsiCrossSection;
81 MCEff = jpsiMCEff;
82 MCEffBoost = jpsiMCEffBoost;
83 boostPhimin = 2.5;
84 boostPhimax = 5;
85 boostMinEmin = 1.2;
86 boostMinEmax = 1.6;
87 }
88 else if ( RunModel == 2 )
89 { // psi2s
90 CrossSection = psi2sCrossSection;
91 MCEff = psi2sMCEff;
92 MCEffBoost = psi2sMCEffBoost;
93 boostPhimin = 2.5;
94 boostPhimax = 5;
95 boostMinEmin = 1.4;
96 boostMinEmax = 1.9;
97 }
98 else if ( RunModel == 3 )
99 { // psi3770
100 CrossSection = psi3770CrossSection;
101 MCEff = psi3770MCEff;
102 MCEffBoost = psi3770MCEffBoost;
103 boostPhimin = 2.5;
104 boostPhimax = 5;
105 boostMinEmin = 1.4;
106 boostMinEmax = 2;
107 }
108
109 tot = 0;
110 signal = 0;
111 boost_signal = 0;
112 boost_tot = 0;
113 StatusCode ssc = service( "THistSvc", m_thsvc );
114 if ( ssc.isFailure() )
115 {
116 log << MSG::FATAL << "DiGamAlg:Couldn't get THistSvc" << endmsg;
117 return ssc;
118 }
119 h_lum = new TH1F( "lum", "lum", 4, 0.5, 4.5 );
120
121 StatusCode scBeamEnergy = service( "BeamEnergySvc", m_BeamEnergySvc );
122
123 if ( scBeamEnergy.isFailure() )
124 {
125 log << MSG::FATAL << "can not use BeamEnergySvc" << endmsg;
126 return scBeamEnergy;
127 }
128
129 StatusCode status;
130
131 NTuplePtr nt2( ntupleSvc(), "DQAFILE/zhsDiGam" );
132 if ( nt2 ) m_tuple2 = nt2;
133 else
134 {
135 m_tuple2 =
136 ntupleSvc()->book( "DQAFILE/zhsDiGam", CLID_ColumnWiseTuple, "ks N-Tuple example" );
137 if ( m_tuple2 )
138 {
139 status = m_tuple2->addItem( "ETol", m_tot );
140 status = m_tuple2->addItem( "maxE", m_maxE );
141 status = m_tuple2->addItem( "minE", m_minE );
142 status = m_tuple2->addItem( "maxTheta", m_maxTheta );
143 status = m_tuple2->addItem( "minTheta", m_minTheta );
144 status = m_tuple2->addItem( "maxPhi", m_maxPhi );
145 status = m_tuple2->addItem( "minPhi", m_minPhi );
146 status = m_tuple2->addItem( "delTheta", m_delTheta );
147 status = m_tuple2->addItem( "delPhi", m_delPhi );
148 status = m_tuple2->addItem( "angle", m_angle );
149 status = m_tuple2->addItem( "boostAngle", m_boostAngle );
150 status = m_tuple2->addItem( "boostMaxE", m_boostMaxE );
151 status = m_tuple2->addItem( "boostMinE", m_boostMinE );
152 status = m_tuple2->addItem( "boostDelPhi", m_boostDelPhi );
153 status = m_tuple2->addItem( "boostDelTheta", m_boostDelTheta );
154 status = m_tuple2->addItem( "boostM", m_boostM );
155 status = m_tuple2->addItem( "boostIM", m_boostIM );
156 status = m_tuple2->addItem( "m", m_m );
157 status = m_tuple2->addItem( "IM", m_IM );
158
159 status = m_tuple2->addItem( "run", m_run );
160 status = m_tuple2->addItem( "rec", m_rec );
161 status = m_tuple2->addItem( "indexmc", m_idxmc, 0, 100 );
162 status = m_tuple2->addIndexedItem( "pdgid", m_idxmc, m_pdgid );
163 status = m_tuple2->addIndexedItem( "motheridx", m_idxmc, m_motheridx );
164 }
165 else
166 {
167 log << MSG::FATAL << "Cannot book N-tuple:" << long( m_tuple2 ) << endmsg;
168 return StatusCode::FAILURE;
169 }
170 }
171
172 return StatusCode::SUCCESS;
173}
174//******************************************************************************
175StatusCode DiGam::execute() {
176 MsgStream log( msgSvc(), name() );
177 log << MSG::INFO << "in execute()" << endmsg;
178 N0++;
179 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
180 runNo = eventHeader->runNumber();
181 int event = eventHeader->eventNumber();
182 log << MSG::INFO << "run,evtnum=" << runNo << "," << event << endmsg;
183 m_run = eventHeader->runNumber();
184 m_rec = eventHeader->eventNumber();
185
186 if ( flag == true )
187 {
188 flag = false;
189 log << MSG::DEBUG << "setting parameter" << endmsg;
190
191 m_BeamEnergySvc->getBeamEnergyInfo();
192 if ( !m_BeamEnergySvc->isRunValid() ) return StatusCode::FAILURE;
193 E_cms = 2 * m_BeamEnergySvc->getbeamE();
194 log << MSG::DEBUG << "cms energy is " << E_cms << endmsg;
195
196 Parameter* func = new Parameter();
197 func->parameters( E_cms );
198
199 CrossSection = func->m_CrossSection;
200 MCEff = func->m_MCEff;
201 MCEffBoost = func->m_MCEffBoost;
202 boostMinEmin = func->m_boostMinEmin;
203 boostMinEmax = func->m_boostMinEmax;
204
205 log << MSG::DEBUG << "parameters: " << CrossSection << " " << MCEff << endmsg;
206 }
207
208 SmartDataPtr<EvtRecEvent> evtRecEvent( eventSvc(), EventModel::EvtRec::EvtRecEvent );
209 log << MSG::INFO << "ncharg,nneu,tottks=" << evtRecEvent->totalCharged() << ","
210 << evtRecEvent->totalNeutral() << "," << evtRecEvent->totalTracks() << endmsg;
211 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol( eventSvc(), EventModel::EvtRec::EvtRecTrackCol );
212 std::vector<int> iGam;
213 iGam.clear();
214 N1++;
215 for ( int i = evtRecEvent->totalCharged(); i < evtRecEvent->totalTracks(); i++ )
216 {
217 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
218 if ( !( *itTrk )->isEmcShowerValid() ) continue;
219 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
220 if ( emcTrk->energy() < 0.6 ) continue;
221 if ( fabs( cos( emcTrk->theta() ) ) > 0.83 ) continue;
222 iGam.push_back( ( *itTrk )->trackId() );
223 }
224 if ( iGam.size() < 2 ) return StatusCode::SUCCESS;
225 N2++;
226 double energy1 = 0.5;
227 double energy2 = 0.5;
228 int gam1 = -9999;
229 int gam2 = -9999;
230 for ( int i = 0; i < iGam.size(); i++ )
231 {
232 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGam[i];
233 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
234 // std::cout<<emcTrk->energy()<<std::endl;
235 if ( emcTrk->energy() > energy1 )
236 { // max emc
237 energy2 = energy1;
238 gam2 = gam1;
239 energy1 = emcTrk->energy();
240 gam1 = iGam[i];
241 }
242 else if ( emcTrk->energy() > energy2 )
243 { // second max emc
244 energy2 = emcTrk->energy();
245 gam2 = iGam[i];
246 }
247 }
248 if ( gam1 == -9999 || gam2 == -9999 ) return StatusCode::SUCCESS;
249 N3++;
250 m_tot = energy1 + energy2;
251 RecEmcShower* maxEmc = ( *( evtRecTrkCol->begin() + gam1 ) )->emcShower();
252 RecEmcShower* minEmc = ( *( evtRecTrkCol->begin() + gam2 ) )->emcShower();
253 m_maxE = maxEmc->energy();
254 m_minE = minEmc->energy();
255 m_maxTheta = maxEmc->theta();
256 m_minTheta = minEmc->theta();
257 m_maxPhi = maxEmc->phi();
258 m_minPhi = minEmc->phi();
259 m_delPhi = ( fabs( m_maxPhi - m_minPhi ) - pai ) * 180. / pai;
260 m_delTheta = ( fabs( m_maxTheta + m_minTheta ) - pai ) * 180. / pai;
261
262 HepLorentzVector max, min;
263 max.setPx( m_maxE * sin( m_maxTheta ) * cos( m_maxPhi ) );
264 max.setPy( m_maxE * sin( m_maxTheta ) * sin( m_maxPhi ) );
265 max.setPz( m_maxE * cos( m_maxTheta ) );
266 max.setE( m_maxE );
267 min.setPx( m_minE * sin( m_minTheta ) * cos( m_minPhi ) );
268 min.setPy( m_minE * sin( m_minTheta ) * sin( m_minPhi ) );
269 min.setPz( m_minE * cos( m_minTheta ) );
270 min.setE( m_minE );
271 m_angle = max.angle( min.vect() ) * 180. / pai;
272 m_m = ( max + min ).m();
273 m_IM = max.invariantMass( min );
274 // select signal and sideband to get lum;
275 if ( m_minE >= boostMinEmin && m_delPhi > dPhiMin && m_delPhi < dPhiMax &&
276 m_minE <= boostMinEmax )
277 tot++;
278 if ( m_minE >= boostMinEmin && m_delPhi > dPhiMinSig && m_delPhi < dPhiMaxSig &&
279 m_minE <= boostMinEmax )
280 signal++;
281
282 // boost
283 // HepLorentzVector p_cms(0.034067,0.0,0.0,3.097);
284 HepLorentzVector boost1 = max.boost( -0.011, 0, 0 );
285 HepLorentzVector boost2 = min.boost( -0.011, 0, 0 );
286 m_boostAngle = boost1.angle( boost2.vect() ) * 180. / pai;
287 m_boostMaxE = boost1.e();
288 m_boostMinE = boost2.e();
289 m_boostDelPhi = ( fabs( boost1.phi() - boost2.phi() ) - pai ) * 180. / pai;
290 m_boostDelTheta = ( fabs( boost1.theta() + boost2.theta() ) - pai ) * 180. / pai;
291 m_boostM = ( boost1 + boost2 ).m();
292 m_boostIM = boost1.invariantMass( boost2 );
293 log << MSG::INFO << "num Good Photon " << iGam.size() << endmsg;
294 if ( m_boostMinE >= boostMinEmin && m_boostMinE <= boostMinEmax &&
295 fabs( m_boostDelPhi ) <= boostPhimin )
296 boost_signal++;
297 if ( m_boostMinE >= boostMinEmin && m_boostMinE <= boostMinEmax &&
298 fabs( m_boostDelPhi ) <= boostPhimax )
299 boost_tot++;
300 // mcTruth
301 if ( eventHeader->runNumber() < 0 )
302 {
303 SmartDataPtr<Event::McParticleCol> mcParticleCol( eventSvc(), "/Event/MC/McParticleCol" );
304 int m_numParticle = 0;
305 if ( !mcParticleCol )
306 {
307 // std::cout << "Could not retrieve McParticelCol" << std::endl;
308 return StatusCode::FAILURE;
309 }
310 else
311 {
312 bool jpsipDecay = false;
313 int rootIndex = -1;
314 Event::McParticleCol::iterator iter_mc = mcParticleCol->begin();
315 for ( ; iter_mc != mcParticleCol->end(); iter_mc++ )
316 {
317 if ( ( *iter_mc )->primaryParticle() ) continue;
318 if ( !( *iter_mc )->decayFromGenerator() ) continue;
319
320 if ( ( *iter_mc )->particleProperty() == 443 )
321 {
322 jpsipDecay = true;
323 rootIndex = ( *iter_mc )->trackIndex();
324 }
325 if ( !jpsipDecay ) continue;
326 int mcidx = ( ( *iter_mc )->mother() ).trackIndex() - rootIndex;
327 int pdgid = ( *iter_mc )->particleProperty();
328 m_pdgid[m_numParticle] = pdgid;
329 m_motheridx[m_numParticle] = mcidx;
330 m_numParticle += 1;
331 }
332 }
333 m_idxmc = m_numParticle;
334 }
335 N4++;
336 m_tuple2->write().ignore();
337 return StatusCode::SUCCESS;
338}
339// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
340StatusCode DiGam::finalize() {
341 MsgStream log( msgSvc(), name() );
342 log << MSG::INFO << "in finalize()" << endmsg;
343 // get intLumEndcapEE
344 double ssg = 0.;
345 double lum = 0.;
346 double boost_ssg = 0.;
347 double boost_lum = 0.;
348 // protection if execute() is not executed
349 if ( flag == false && ( CrossSection * MCEff == 0 || CrossSection * MCEffBoost == 0 ) )
350 {
351 log << MSG::FATAL << "DiGam Error: cross_section or MC_eff is not correct!" << endmsg;
352 exit( 1 );
353 }
354 if ( flag == false )
355 {
356 ssg = ( signal - ( tot - signal ) ) * 1.0;
357 // boss 650 mc:7.135%
358 lum = ( ssg ) / ( CrossSection * MCEff );
359 // boost LUM
360 boost_ssg = ( boost_signal - ( boost_tot - boost_signal ) ) * 1.0;
361 // boost 650 mc:7.097%
362 boost_lum = ( boost_ssg ) / ( CrossSection * MCEffBoost );
363 }
364
365 // std::cout<<"flag = "<<flag<<std::endl;
366 log << MSG::DEBUG << "E_cms = " << E_cms << endmsg;
367
368 log << MSG::DEBUG << "N0 = " << N0 << endmsg;
369 log << MSG::DEBUG << "minE, phi in:" << boostMinEmin << " ~ " << boostMinEmax << ", "
370 << boostPhimin << " ~ " << boostPhimax << endmsg;
371 log << MSG::DEBUG << "sigma, eff = " << CrossSection << "," << MCEff << endmsg;
372 log << MSG::DEBUG << "sigma, effBoost = " << CrossSection << ", " << MCEffBoost << endmsg;
373 log << MSG::DEBUG << "Nsignal, lum, Nsig_boost, boost_lum = " << ssg << ", " << lum << ", "
374 << boost_ssg << ", " << boost_lum << endmsg;
375 h_lum->SetBinContent( 1, lum );
376 h_lum->SetBinContent( 2, ssg );
377 h_lum->SetBinContent( 3, boost_lum );
378 h_lum->SetBinContent( 4, boost_ssg );
379 if ( m_thsvc->regHist( "/DQAHist/zhsLUM/lum", h_lum ).isFailure() )
380 {
381 log << MSG::FATAL << "DiGam Error:can't write data to LUM!" << endmsg;
382 exit( 1 );
383 }
384 return StatusCode::SUCCESS;
385}
pai
const double pai
Definition BbEmc.cxx:47
DECLARE_COMPONENT
DECLARE_COMPONENT(BesBdkRc)
lum
Double_t lum[10]
Definition DataBase/tau_mode.c:35
N0
int N0
Definition DiGam.cxx:35
N3
int N3
Definition DiGam.cxx:38
N1
int N1
Definition DiGam.cxx:36
N4
int N4
Definition DiGam.cxx:39
N2
int N2
Definition DiGam.cxx:37
min
#define min(a, b)
Definition Eepipi/src/ee2eepp/basesv5.1/f2c.h:157
max
#define max(a, b)
Definition Eepipi/src/ee2eepp/basesv5.1/f2c.h:158
EvtRecTrackIterator
EvtRecTrackCol::iterator EvtRecTrackIterator
Definition Event/EvtRecEvent/include/EvtRecEvent/EvtRecTrack.h:106
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
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
DiGam
Definition DiGam.h:16
DiGam::execute
StatusCode execute()
Definition DiGam.cxx:175
DiGam::DiGam
DiGam(const std::string &name, ISvcLocator *pSvcLocator)
Definition DiGam.cxx:41
DiGam::finalize
StatusCode finalize()
Definition DiGam.cxx:340
DiGam::initialize
StatusCode initialize()
Definition DiGam.cxx:74
Parameter::m_MCEffBoost
double m_MCEffBoost
Definition Parameter.h:11
Parameter::m_boostMinEmax
double m_boostMinEmax
Definition Parameter.h:13
Parameter::parameters
void parameters(double E_cms)
Definition Parameter.cxx:8
Parameter::m_boostMinEmin
double m_boostMinEmin
Definition Parameter.h:12
Parameter::m_MCEff
double m_MCEff
Definition Parameter.h:10
Parameter::m_CrossSection
double m_CrossSection
Definition Parameter.h:9