BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
SD0Tag Class Reference

#include <SD0Tag.h>

Inheritance diagram for SD0Tag:

Public Member Functions

 SD0Tag (const std::string &name, ISvcLocator *pSvcLocator)
StatusCode initialize ()
StatusCode execute ()
StatusCode finalize ()

Detailed Description

Definition at line 13 of file SD0Tag.h.

Constructor & Destructor Documentation

◆ SD0Tag()

SD0Tag::SD0Tag ( const std::string & name,
ISvcLocator * pSvcLocator )

Definition at line 67 of file SD0Tag.cxx.

68 : Algorithm( name, pSvcLocator ) {
69 // Declare the properties
70 declareProperty( "PID_FLAG", PID_flag = 1 );
71 declareProperty( "MC_sample", m_MC_sample = 1 );
72
73 declareProperty( "m_Seperate_Charge", Seperate_Charge = 1 );
74 declareProperty( "m_Charge_default", Charge_default = 1 );
75}

Referenced by SD0Tag().

Member Function Documentation

◆ execute()

StatusCode SD0Tag::execute ( )

Definition at line 131 of file SD0Tag.cxx.

131 {
132
133 MsgStream log( msgSvc(), name() );
134
135 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
136 int runNo = eventHeader->runNumber();
137 int event = eventHeader->eventNumber();
138
139 SmartDataPtr<EvtRecEvent> evtRecEvent( eventSvc(), EventModel::EvtRec::EvtRecEvent );
140 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol( eventSvc(), EventModel::EvtRec::EvtRecTrackCol );
141 // cout<<"//////////////////////////////////////"<<endl;
142 nD0++;
143
144 //
145 // The default energy is 3.773 GeV for psi(3770) data.
146 // Alayizer can add calibrated energy here.
147 double Ebeam = 3.773 / 2.0;
148
149 if ( runNo >= 11414 && runNo <= 23500 )
150 {
151 for ( int i = 0; i < 3467; i++ )
152 {
153 if ( runNo == p_run[i] ) { Ebeam = p_Ecm[i] / 2.0; }
154 }
155 }
156
157 if ( runNo < 0 )
158 {
159 int irun = abs( runNo );
160 // Ecmset* ECMSET = Ecmset::instance();
161 // ECMSET->EcmSet(irun);
162 // Ebeam = ECMSET->ECM()/2.0;
163 }
164
165 if ( nD0 % 1000 == 0 )
166 cout << "SD0TagAlg-13-11-15pp = " << nD0 << " " << runNo << " " << event << " "
167 << Ebeam * 2 << endl;
168
169 Hep3Vector xorigin( 0, 0, 0 );
170 IVertexDbSvc* vtxsvc;
171 Gaudi::svcLocator()->service( "VertexDbSvc", vtxsvc );
172 if ( vtxsvc->isVertexValid() )
173 {
174 double* dbv = vtxsvc->PrimaryVertex();
175 double* vv = vtxsvc->SigmaPrimaryVertex();
176 xorigin.setX( dbv[0] );
177 xorigin.setY( dbv[1] );
178 xorigin.setZ( dbv[2] );
179 }
180
181 /////////////////////////////////////
182 Vint iCharge_good;
183 iCharge_good.clear();
184 for ( int i = 0; i < evtRecEvent->totalCharged(); i++ )
185 {
186 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
187 if ( !( *itTrk )->isMdcTrackValid() ) continue;
188 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
189 //%%%%%%%%%%%%%%%%%%%%%%%%
190 HepVector a = mdcTrk->helix();
191 HepSymMatrix Ea = mdcTrk->err();
192 HepPoint3D point0( 0., 0., 0. ); // the initial point for MDC recosntruction
193 HepPoint3D IP( xorigin[0], xorigin[1], xorigin[2] );
194 VFHelix helixip( point0, a, Ea );
195 helixip.pivot( IP );
196 HepVector vecipa = helixip.a();
197 double Rvxy0 = fabs( vecipa[0] ); // the nearest distance to IP in xy plane
198 double Rvz0 = vecipa[3]; // the nearest distance to IP in z direction
199 double costheta = cos( mdcTrk->theta() );
200 //%%%%%%%%%%%%%%%%%%%%%%%%
201 if ( fabs( Rvxy0 ) >= 1.0 ) continue;
202 if ( fabs( Rvz0 ) >= 10.0 ) continue;
203 if ( fabs( costheta ) >= 0.930 ) continue;
204 iCharge_good.push_back( i );
205 }
206
207 ///////////////////////////////////////////////////////
208 Vint iGam;
209 iGam.clear();
210 for ( int i = evtRecEvent->totalCharged(); i < evtRecEvent->totalTracks(); i++ )
211 {
212 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
213 if ( !( *itTrk )->isEmcShowerValid() ) continue;
214 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
215 //
216 // We here remove the hot channels of EMC
217 //
218 if ( abs( runNo ) >= 11615 && abs( runNo ) <= 11617 && emcTrk->cellId() == 805438481 )
219 continue;
220 if ( abs( runNo ) >= 11615 && abs( runNo ) <= 11655 && emcTrk->cellId() == 805376008 )
221 continue;
222 if ( abs( runNo ) >= 13629 && abs( runNo ) <= 13669 && emcTrk->cellId() == 805374754 )
223 continue;
224 if ( abs( runNo ) >= 21190 && abs( runNo ) <= 21231 && emcTrk->cellId() == 805375262 )
225 continue;
226
227 int emctdc = emcTrk->time();
228 if ( emctdc > 14 || emctdc < 0 ) continue;
229
230 double eraw = emcTrk->energy();
231 double theta = emcTrk->theta();
232 int Module = 0;
233 if ( fabs( cos( theta ) ) < 0.80 && eraw > 0.025 ) { Module = 1; }
234 if ( fabs( cos( theta ) ) > 0.86 && fabs( cos( theta ) ) < 0.92 && eraw > 0.05 )
235 { Module = 2; }
236 if ( Module == 0 ) continue;
237 iGam.push_back( ( *itTrk )->trackId() );
238 }
239 /////////////////////////////////////////////////////////////////
240 DTagTool trk;
241 bool cosmic_ok = trk.cosmicandleptonVeto();
242 /////////////////////////////////////////////////////////////////
243
244 int ntk;
245 double tagmass, ebeam, tagmode, ksmass, dlte;
246
247 Vint tagtrk;
248 tagtrk.clear();
249 Vint tagGam;
250 tagGam.clear();
251
252 HepLorentzVector tagp;
253
254 double mass_bc_temp, mass_kf_temp, kf_chi2_temp, mks_temp, mpi0_temp, ptag_temp;
255 int Charge_candidate_D = 0;
256 double EGam_max_0 = 0;
257
258 for ( int i1 = 0; i1 < 2; i1++ )
259 {
260 if ( Seperate_Charge == 2 )
261 {
262 Charge_candidate_D = Charge_default;
263 i1 = 2;
264 }
265 if ( Seperate_Charge == 1 ) { Charge_candidate_D = pow( -1, i1 ); }
266 if ( Seperate_Charge == 0 )
267 {
268 Charge_candidate_D = 0;
269 i1 = 2;
270 }
271
272 for ( int i = 0; i < 15; i++ )
273 {
274 EGam_max_0 = 0;
275 int mdslct = pow( 2., i );
276 Sing sing;
277 sing.Mdset( event, evtRecTrkCol, iCharge_good, iGam, mdslct, Ebeam, PID_flag,
278 Charge_candidate_D );
279 bool oktag = sing.Getoktg();
280 if ( oktag )
281 {
282 //
283 // Here analysizer should pick up variables from anti-D0 tags
284 // to do physics analysis
285 //
286 tagtrk = sing.Gettagtrk1();
287 tagmode = sing.Gettagmd();
288 //==========================================================================
289 if ( ( abs( tagmode ) == 11 || abs( tagmode ) == 15 || abs( tagmode ) == 19 ) )
290 {
291 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
292 for ( int i = evtRecEvent->totalCharged(); i < evtRecEvent->totalTracks(); i++ )
293 {
294 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
295 int itrk = ( *itTrk )->trackId();
296 if ( !( *itTrk )->isEmcShowerValid() ) continue;
297 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
298
299 int emctdc = emcTrk->time();
300 if ( emctdc > 14 || emctdc < 0 ) continue;
301
302 Hep3Vector emcpos( emcTrk->x(), emcTrk->y(), emcTrk->z() );
303 double dang = 200.;
304 for ( int j = 0; j < tagtrk.size(); j++ )
305 {
306 EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() + tagtrk[j];
307 if ( !( *jtTrk )->isExtTrackValid() ) continue;
308 RecExtTrack* extTrk = ( *jtTrk )->extTrack();
309 if ( extTrk->emcVolumeNumber() == -1 ) continue;
310 Hep3Vector extpos = extTrk->emcPosition();
311 double angd = extpos.angle( emcpos );
312 if ( angd < dang ) dang = angd;
313 }
314 dang = dang * 180 / ( CLHEP::pi );
315 if ( fabs( dang ) < 20. ) continue;
316
317 double eraw = emcTrk->energy();
318 if ( eraw > EGam_max_0 ) EGam_max_0 = eraw;
319 }
320 }
321 //==========================================================================
322 m_cosmic_ok = cosmic_ok;
323 m_EGam_max_0 = EGam_max_0;
324 m_nGood = iCharge_good.size();
325 m_nGam = iGam.size();
326 m_runNo = runNo;
327 m_event = event;
328
329 m_tagmode = sing.Gettagmd();
330 m_mass_bc = sing.Getmass_bc();
331 m_delE_tag = sing.GetdelE_tag();
332
333 m_tuple1->write();
334
335 // Here one can do double tag analysis based on the variables
336 // from single anti-D0 tag analysis.
337 // To do so, analyzer should make his/her codes
338 //
339 }
340 }
341 } // The end of loopping over mdslct
342
343 return StatusCode::SUCCESS;
344}
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition CosmicGenerator.h:36
runNo
int runNo
Definition DQA_TO_DB.cxx:13
EvtRecTrackIterator
EvtRecTrackCol::iterator EvtRecTrackIterator
Definition Event/EvtRecEvent/include/EvtRecEvent/EvtRecTrack.h:106
ebeam
double ebeam
Definition EvtPhokharaDef.hh:18
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187
nD0
int nD0
Definition SD0Tag.cxx:57
Vint
std::vector< int > Vint
Definition SD0Tag.h:11
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:12
DTagTool::cosmicandleptonVeto
bool cosmicandleptonVeto(bool emc=true)
Definition DTagTool.cxx:1187
DstMdcTrack::err
const HepSymMatrix err() const
Definition DstMdcTrack.cxx:149
DstMdcTrack::helix
const HepVector helix() const
......
Definition DstMdcTrack.cxx:147
IVertexDbSvc::isVertexValid
virtual bool isVertexValid()=0
IVertexDbSvc::SigmaPrimaryVertex
virtual double * SigmaPrimaryVertex()=0
IVertexDbSvc::PrimaryVertex
virtual double * PrimaryVertex()=0
Sing::Gettagmd
double Gettagmd()
Definition Sing.h:27
Sing::GetdelE_tag
double GetdelE_tag()
Definition Sing.h:30
Sing::Getmass_bc
double Getmass_bc()
Definition Sing.h:28
Sing::Mdset
void Mdset(double event, SmartDataPtr< EvtRecTrackCol > evtRecTrkCol, Vint iGood, Vint iGam, int mdset, double Ebeam, int PID_flag, int Charge_candidate_D)
Definition Sing.cxx:51
Sing::Getoktg
bool Getoktg()
Definition Sing.h:26
Sing::Gettagtrk1
Vint Gettagtrk1()
Definition Sing.h:31

◆ finalize()

StatusCode SD0Tag::finalize ( )

Definition at line 350 of file SD0Tag.cxx.

350 {
351 MsgStream log( msgSvc(), name() );
352 cout << "nD0 = " << nD0 << endl;
353 cout << "n1 = " << n1 << endl;
354 cout << "n2 = " << n2 << endl;
355 cout << "ND0 ==== " << ND0 << endl;
356 cout << "NDp ==== " << NDp << endl;
357 log << MSG::INFO << "in finalize()" << endmsg;
358 return StatusCode::SUCCESS;
359}
NDp
int NDp
Definition SD0Tag.cxx:61
ND0
int ND0
Definition SD0Tag.cxx:60
n2
int n2
Definition SD0Tag.cxx:59
n1
int n1
Definition SD0Tag.cxx:58

◆ initialize()

StatusCode SD0Tag::initialize ( )

Definition at line 79 of file SD0Tag.cxx.

79 {
80
81 MsgStream log( msgSvc(), name() );
82
83 log << MSG::INFO << "in initialize()" << endmsg;
84
85 StatusCode status;
86
87 ifstream readrunEcm;
88 //
89 // Read in the energy calibration constant for run-by-run data
90 readrunEcm.open( "../share/psipp_ecms_calrunbyrun_runno_3648runs.dat" );
91 cout << "readrunEcm = " << readrunEcm.is_open() << endl;
92 for ( int i = 0; i < 3467; i++ )
93 {
94 readrunEcm >> p_run[i];
95 readrunEcm >> p_Ecm[i];
96 }
97
98 NTuplePtr nt1( ntupleSvc(), "FILE1/tag" );
99 if ( nt1 ) m_tuple1 = nt1;
100 else
101 {
102 m_tuple1 = ntupleSvc()->book( "FILE1/tag", CLID_ColumnWiseTuple, "tag N-Tuple example" );
103 if ( m_tuple1 )
104 {
105 status = m_tuple1->addItem( "tagmode", m_tagmode );
106 status = m_tuple1->addItem( "mass_bc", m_mass_bc );
107 status = m_tuple1->addItem( "delE_tag", m_delE_tag );
108
109 status = m_tuple1->addItem( "cosmic_ok", m_cosmic_ok );
110 status = m_tuple1->addItem( "EGam_max_0", m_EGam_max_0 );
111 status = m_tuple1->addItem( "nGood", m_nGood );
112 status = m_tuple1->addItem( "nGam", m_nGam );
113 status = m_tuple1->addItem( "runNo", m_runNo );
114 status = m_tuple1->addItem( "event", m_event );
115 }
116 else
117 {
118 log << MSG::ERROR << " Cannot book N-tuple:" << long( m_tuple1 ) << endmsg;
119 return StatusCode::FAILURE;
120 }
121 }
122
123 log << MSG::INFO << "successfully return from initialize()" << endmsg;
124 return StatusCode::SUCCESS;
125}
ntupleSvc
INTupleSvc * ntupleSvc()
Definition ServiceAccessor.h:51
The documentation for this class was generated from the following files: