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