BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
FarmMonitorAlg.cxx
Go to the documentation of this file.
1#include "GaudiKernel/IDataProviderSvc.h"
2#include "GaudiKernel/MsgStream.h"
3#include "GaudiKernel/PropertyMgr.h"
4#include "GaudiKernel/SmartDataPtr.h"
5
6#include "EventModel/Event.h"
7#include "EventModel/EventHeader.h"
8#include "EventModel/EventModel.h"
9#include "EvtRecEvent/EvtRecEvent.h"
10#include "EvtRecEvent/EvtRecTrack.h"
11#include "EvtRecEvent/EvtRecVeeVertex.h"
12
13#include "GaudiKernel/Bootstrap.h"
14#include "GaudiKernel/IHistogramSvc.h"
15#include "GaudiKernel/INTupleSvc.h"
16#include "GaudiKernel/NTuple.h"
17#include "TMath.h"
18
19#include "DstEvent/TofHitStatus.h"
20
21#include "CLHEP/Vector/LorentzVector.h"
22#include "CLHEP/Vector/ThreeVector.h"
23
24#include "FarmMonitorAlg.h"
25
26using namespace std;
27using namespace Event;
28
29double mPi = 0.13957;
30double mKs = 0.49767;
31double mLambda = 1.11568;
32
33/////////////////////////////////////////////////////////////////////////////
34DECLARE_COMPONENT( FarmMonitorAlg )
35FarmMonitorAlg::FarmMonitorAlg( const std::string& name, ISvcLocator* pSvcLocator )
36 : Algorithm( name, pSvcLocator ) {
37
38 // Declare the properties
39 declareProperty( "PrintRunEventFreq", m_RunEventFreq = 10 );
40
41 declareProperty( "Ecm", m_ecm = 3.686 );
42 declareProperty( "Vr0cut", m_vr0cut = 1.0 );
43 declareProperty( "Vz0cut", m_vz0cut = 5.0 );
44}
45
46// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
47StatusCode FarmMonitorAlg::initialize() {
48 MsgStream log( msgSvc(), name() );
49
50 log << MSG::INFO << "in initialize()" << endmsg;
51
52 StatusCode status;
53
54 // histograms
55
56 /// Total energies
57 h_eVisibleDivEcm = histoSvc()->book( "eVisibleDivEcm", 1, "eVisible/Ecm", 150, 0.0, 1.5 );
58 h_eEMCDivEcm = histoSvc()->book( "eEMCDivEcm", 1, "eEMC/Ecm", 150, 0.0, 1.5 );
59 h_eNeutralDivEcm = histoSvc()->book( "eNeutralDivEcm", 1, "eNeutral/Ecm", 150, 0.0, 1.5 );
60 h_eChargedDivEcm = histoSvc()->book( "eChargedDivEcm", 1, "eCharged/Ecm", 150, 0.0, 1.5 );
61
62 /// Net momenta
63 h_netMomentumDivEcm_AllChargedTracks = histoSvc()->book(
64 "netMomentumDivEcm_AllChargedTracks", 1, "#Sigma p/Ecm (Charged)", 100, 0.0, 1.0 );
65 h_netTransMomentumDivEcm_AllChargedTracks =
66 histoSvc()->book( "netTransMomentumDivEcm_AllChargedTracks", 1,
67 "#Sigma ptr/Ecm (Charged)", 100, 0.0, 1.0 );
68 h_cosTheta_AllChargedTracks = histoSvc()->book( "cosTheta_AllChargedTracks", 1,
69 "cos#theta (Charged)", 100, -1.0, 1.0 );
70
71 h_netMomentumDivEcm_AllNeutralTracks = histoSvc()->book(
72 "netMomentumDivEcm_AllNeutralTracks", 1, "#Sigma p/Ecm (Neutral)", 100, 0.0, 1.0 );
73 h_netTransMomentumDivEcm_AllNeutralTracks =
74 histoSvc()->book( "netTransMomentumDivEcm_AllNeutralTracks", 1,
75 "#Sigma ptr/Ecm (Neutral)", 100, 0.0, 1.0 );
76 h_cosTheta_AllNeutralTracks = histoSvc()->book( "cosTheta_AllNeutralTracks", 1,
77 "cos#theta (Neutral)", 100, -1.0, 1.0 );
78
79 h_netMomentumDivEcm_AllTracks =
80 histoSvc()->book( "netMomentumDivEcm_AllTracks", 1, "#Sigma p/Ecm", 100, 0.0, 1.0 );
81 h_netTransMomentumDivEcm_AllTracks = histoSvc()->book( "netTransMomentumDivEcm_AllTracks", 1,
82 "#Sigma ptr/Ecm", 100, 0.0, 1.0 );
83 h_cosTheta_AllTracks =
84 histoSvc()->book( "cosTheta_AllTracks", 1, "cos#theta", 100, -1.0, 1.0 );
85
86 /// Charged Tracks
87 h_trackR0 = histoSvc()->book( "trackR0", 1, "R0 (cm)", 100, 0.0, 2.0 );
88 h_trackZ0 = histoSvc()->book( "trackZ0", 1, "Z0 (cm)", 100, -20.0, 20.0 );
89
90 h_nChargedTracks = histoSvc()->book( "nChargedTracks", 1, "NChargedTracks", 17, -0.5, 16.5 );
91 h_nChargedTracksIP =
92 histoSvc()->book( "nChargedTracksIP", 1, "NChargedTracks from IP", 17, -0.5, 16.5 );
93
94 h_netCharge = histoSvc()->book( "netCharge", 1, "Net Charge", 11, -5.5, 5.5 );
95 h_netChargeIP =
96 histoSvc()->book( "netChargeIP", 1, "Net Charge with IP tracks", 11, -5.5, 5.5 );
97
98 /// 2 highest momentum charged tracks
99 h_pIPTrack1DivEb = histoSvc()->book( "pIPTrack1DivEb", 1, "pTrack1/Ebeam", 150, 0.0, 1.5 );
100 h_pIPTrack2DivEb = histoSvc()->book( "pIPTrack2DivEb", 1, "pTrack2/Ebeam", 150, 0.0, 1.5 );
101
102 h_eEMCIPTrack1DivEb =
103 histoSvc()->book( "eEMCIPTrack1DivEb", 1, "eEMCTrack1/Ebeam", 150, 0.0, 1.5 );
104 h_eEMCIPTrack2DivEb =
105 histoSvc()->book( "eEMCIPTrack2DivEb", 1, "eEMCTrack2/Ebeam", 150, 0.0, 1.5 );
106
107 h_acoplanarity_2HighestPIPTracks =
108 histoSvc()->book( "acoplanarity_2HighestPIPTracks", 1, "acoplanarity", 158, 0.0, 3.16 );
109
110 /// Neutral Tracks
111 h_nNeutralTracks = histoSvc()->book( "nNeutralTracks", 1, "NNeutralTracks", 31, -0.5, 30.5 );
112 h_nNeutralTracksGT30MeV =
113 histoSvc()->book( "nNeutralTracksGT30MeV", 1, "NNeutralTracksGT30MeV", 31, -0.5, 30.5 );
114
115 h_eEMCShower1DivEb =
116 histoSvc()->book( "eEMCShower1DivEb", 1, "eEMCShower1/Ebeam", 150, 0.0, 1.5 );
117 h_eEMCShower2DivEb =
118 histoSvc()->book( "eEMCShower2DivEb", 1, "eEMCShower2/Ebeam", 150, 0.0, 1.5 );
119 h_eEMCShower3DivEb =
120 histoSvc()->book( "eEMCShower3DivEb", 1, "eEMCShower3/Ebeam", 150, 0.0, 1.5 );
121
122 /// MUC information
123 h_mucDepth = histoSvc()->book( "mucDepth", 1, "mucDepth", 200, 0.0, 200.0 );
124 h_mucDepthVsCosTheta = histoSvc()->book( "mucDepthVsCosTheta", 2, "mucDepthVsCosTheta", 100,
125 -1.0, 1.0, 200, 0.0, 200.0 );
126 h_mucDepthVsPhi = histoSvc()->book( "mucDepthVsPhi", 2, "mucDepthVsPhi", 180, -3.15, 3.15,
127 200, 0.0, 200.0 );
128
129 /// PID (dE/dx) information
130 h_dedxTotalHitsIP =
131 histoSvc()->book( "dedxTotalHitsIP", 1, "dedxTotalHitsIP", 70, 0.0, 70.0 );
132 h_dedxGoodHitsIP = histoSvc()->book( "dedxGoodHitsIP", 1, "dedxGoodHitsIP", 70, 0.0, 70.0 );
133 h_dedxElecChiIP = histoSvc()->book( "dedxElecChiIP", 1, "dedxElecChiIP", 120, -6.0, 6.0 );
134 h_dedxMuonChiIP = histoSvc()->book( "dedxMuonChiIP", 1, "dedxMuonChiIP", 120, -6.0, 6.0 );
135 h_dedxPionChiIP = histoSvc()->book( "dedxPionChiIP", 1, "dedxPionChiIP", 120, -6.0, 6.0 );
136 h_dedxKaonChiIP = histoSvc()->book( "dedxKaonChiIP", 1, "dedxKaonChiIP", 120, -6.0, 6.0 );
137 h_dedxProtonChiIP =
138 histoSvc()->book( "dedxProtonChiIP", 1, "dedxProtonChiIP", 120, -6.0, 6.0 );
139 h_dedxProbPHIP = histoSvc()->book( "dedxProbPHIP", 1, "dedxProbPHIP", 200, 0.0, 2000.0 );
140 h_dedxProbPHVsMomentumIP = histoSvc()->book(
141 "dedxProbPHVsMomentumIP", 2, "dedxProbPHVsMomentumIP", 100, 0.0, 2.0, 200, 0.0, 2000.0 );
142
143 /// PID (TOF) information
144 h_tofPHIP_BarrelLayer1 =
145 histoSvc()->book( "tofPHIP_BarrelLayer1", 1, "tofPHIP_BarrelLayer1", 250, 0.0, 5000.0 );
146 h_tofPHIP_BarrelLayer2 =
147 histoSvc()->book( "tofPHIP_BarrelLayer2", 1, "tofPHIP_BarrelLayer2", 250, 0.0, 5000.0 );
148 h_tofPHIP_EastEndcap =
149 histoSvc()->book( "tofPHIP_EastEndcap", 1, "tofPHIP_EastEndcap", 200, 0.0, 4000.0 );
150 h_tofPHIP_WestEndcap =
151 histoSvc()->book( "tofPHIP_WestEndcap", 1, "tofPHIP_WestEndcap", 200, 0.0, 4000.0 );
152 h_tofIP_BarrelLayer1 =
153 histoSvc()->book( "tofIP_BarrelLayer1", 1, "tofIP_BarrelLayer1", 100, 0.0, 20.0 );
154 h_tofIP_BarrelLayer2 =
155 histoSvc()->book( "tofIP_BarrelLayer2", 1, "tofIP_BarrelLayer1", 100, 0.0, 20.0 );
156 h_tofIP_EastEndcap =
157 histoSvc()->book( "tofIP_EastEndcap", 1, "tofIP_EastEndcap", 100, 0.0, 20.0 );
158 h_tofIP_WestEndcap =
159 histoSvc()->book( "tofIP_WestEndcap", 1, "tofIP_WestEndcap", 100, 0.0, 20.0 );
160 h_tofElecIP_Barrel =
161 histoSvc()->book( "tofElecIP_Barrel", 1, "tofElecIP_Barrel", 120, -6.0, 6.0 );
162 h_tofMuonIP_Barrel =
163 histoSvc()->book( "tofMuonIP_Barrel", 1, "tofMuonIP_Barrel", 120, -6.0, 6.0 );
164 h_tofPionIP_Barrel =
165 histoSvc()->book( "tofPionIP_Barrel", 1, "tofPionIP_Barrel", 120, -6.0, 6.0 );
166 h_tofKaonIP_Barrel =
167 histoSvc()->book( "tofKaonIP_Barrel", 1, "tofKaonIP_Barrel", 120, -6.0, 6.0 );
168 h_tofProtonIP_Barrel =
169 histoSvc()->book( "tofProtonIP_Barrel", 1, "tofProtonIP_Barrel", 120, -6.0, 6.0 );
170 h_tofElecIP_Endcap =
171 histoSvc()->book( "tofElecIP_Endcap", 1, "tofElecIP_Endcap", 120, -6.0, 6.0 );
172 h_tofMuonIP_Endcap =
173 histoSvc()->book( "tofMuonIP_Endcap", 1, "tofMuonIP_Endcap", 120, -6.0, 6.0 );
174 h_tofPionIP_Endcap =
175 histoSvc()->book( "tofPionIP_Endcap", 1, "tofPionIP_Endcap", 120, -6.0, 6.0 );
176 h_tofKaonIP_Endcap =
177 histoSvc()->book( "tofKaonIP_Endcap", 1, "tofKaonIP_Endcap", 120, -6.0, 6.0 );
178 h_tofProtonIP_Endcap =
179 histoSvc()->book( "tofProtonIP_Endcap", 1, "tofProtonIP_Endcap", 120, -6.0, 6.0 );
180 h_tofVsMomentumIP = histoSvc()->book( "tofVsMomentumIP", 2, "tofVsMomentumIP", 100, 0.0, 2.0,
181 100, 0.0, 20.0 );
182
183 /// VeeVertex information
184 h_nKs = histoSvc()->book( "nKs", 1, "nKs", 9, -0.5, 8.5 );
185 h_ksMass = histoSvc()->book( "ksMass", 1, "ksMass", 70, 0.480, 0.515 );
186
187 h_nLambda = histoSvc()->book( "nLambda", 1, "nLambda", 9, -0.5, 8.5 );
188 h_lambdaMass = histoSvc()->book( "lambdaMass", 1, "lambdaMass", 76, 1.100, 1.133 );
189
190 // Add Pi0 histograms when official pi0 list is ready
191
192 log << MSG::INFO << "successfully return from initialize()" << endmsg;
193 return StatusCode::SUCCESS;
194}
195
196// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
197StatusCode FarmMonitorAlg::execute() {
198
199 double eBeam = m_ecm / 2;
200
201 MsgStream log( msgSvc(), name() );
202 log << MSG::INFO << "in execute()" << endmsg;
203
204 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
205 int run = eventHeader->runNumber();
206 int event = eventHeader->eventNumber();
207 if ( event % m_RunEventFreq == 0 )
208 std::cout << "Run " << run << ", event " << event << std::endl;
209
210 SmartDataPtr<EvtRecEvent> evtRecEvent( eventSvc(), EventModel::EvtRec::EvtRecEvent );
211 log << MSG::DEBUG << "ncharg, nneu, tottks = " << evtRecEvent->totalCharged() << " , "
212 << evtRecEvent->totalNeutral() << " , " << evtRecEvent->totalTracks() << endmsg;
213
214 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol( eventSvc(), EventModel::EvtRec::EvtRecTrackCol );
215
216 ////////////////////////////////////////
217 /// Variables for charged track FOR loop
218 int nChargedTracks = 0, nChargedTracksIP = 0;
219 int nCharge = 0, nChargeIP = 0;
220 double totalVisibleEnergy = 0, totalChargedEnergy = 0, totalEMCEnergy = 0;
221 double totalChargedPX = 0, totalChargedPY = 0, totalChargedPZ = 0;
222
223 double highestIPTrackP = -1, secondHighestIPTrackP = -2;
224 int highestPIPTrackId = -1, secondHighestPIPTrackId = -1;
225
226 /// Charged track FOR loop
227 for ( int i = 0; i < evtRecEvent->totalCharged(); i++ )
228 {
229 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
230 if ( !( *itTrk )->isMdcTrackValid() ) continue;
231 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
232
233 int trackId = mdcTrk->trackId();
234 int charge = mdcTrk->charge();
235 double r0 = mdcTrk->r();
236 double z0 = mdcTrk->z();
237 h_trackR0->fill( r0 );
238 h_trackZ0->fill( z0 );
239
240 nChargedTracks++;
241 nCharge += charge;
242
243 double pX = mdcTrk->px();
244 double pY = mdcTrk->py();
245 double pZ = mdcTrk->pz();
246 double pMag = mdcTrk->p();
247 double cosTheta = cos( mdcTrk->theta() );
248 double phi = mdcTrk->phi();
249
250 double chargedEnergy = sqrt( pMag * pMag + mPi * mPi );
251 totalVisibleEnergy += chargedEnergy;
252 totalChargedEnergy += chargedEnergy;
253
254 totalChargedPX += pX;
255 totalChargedPY += pY;
256 totalChargedPZ += pZ;
257
258 /// EMC energy associated with track
259 if ( ( *itTrk )->isEmcShowerValid() )
260 {
261 RecEmcShower* emcChargedTrk = ( *itTrk )->emcShower();
262 totalEMCEnergy += emcChargedTrk->energy();
263 }
264
265 /// MUC information
266 if ( ( *itTrk )->isMucTrackValid() )
267 {
268 RecMucTrack* mucTrk = ( *itTrk )->mucTrack();
269 double mucDepth = mucTrk->depth();
270 if ( mucDepth > 0 )
271 {
272 h_mucDepth->fill( mucDepth );
273 h_mucDepthVsCosTheta->fill( cosTheta, mucDepth );
274 h_mucDepthVsPhi->fill( phi, mucDepth );
275 }
276 } // End of "isMucShowerValid()" IF
277
278 ////////////////////////
279 /// Tracks after IP cuts
280 if ( fabs( z0 ) >= m_vz0cut ) continue;
281 if ( r0 >= m_vr0cut ) continue;
282
283 nChargedTracksIP++;
284 nChargeIP += charge;
285
286 /// dE/dx information
287 double dedxProbPH = -10;
288 int dedxNumTotalHits = -10;
289 int dedxNumGoodHits = -10;
290
291 if ( ( *itTrk )->isMdcDedxValid() )
292 {
293 RecMdcDedx* dedxTrk = ( *itTrk )->mdcDedx();
294
295 dedxNumTotalHits = dedxTrk->numTotalHits();
296 dedxNumGoodHits = dedxTrk->numGoodHits();
297 h_dedxTotalHitsIP->fill( dedxNumTotalHits );
298 h_dedxGoodHitsIP->fill( dedxNumGoodHits );
299
300 h_dedxElecChiIP->fill( dedxTrk->chiE() );
301 h_dedxMuonChiIP->fill( dedxTrk->chiMu() );
302 h_dedxPionChiIP->fill( dedxTrk->chiPi() );
303 h_dedxKaonChiIP->fill( dedxTrk->chiK() );
304 h_dedxProtonChiIP->fill( dedxTrk->chiP() );
305
306 dedxProbPH = dedxTrk->probPH();
307 h_dedxProbPHIP->fill( dedxProbPH );
308 h_dedxProbPHVsMomentumIP->fill( pMag, dedxProbPH );
309
310 } // End of "isMdcDedxValid()" IF
311
312 /// TOF information
313 if ( ( *itTrk )->isTofTrackValid() )
314 {
315 SmartRefVector<RecTofTrack> tofTrkCol = ( *itTrk )->tofTrack();
316 SmartRefVector<RecTofTrack>::iterator iter_tof = tofTrkCol.begin();
317
318 for ( ; iter_tof != tofTrkCol.end(); iter_tof++ )
319 {
320 TofHitStatus* status = new TofHitStatus;
321 status->setStatus( ( *iter_tof )->status() );
322
323 if ( status->is_barrel() )
324 {
325 if ( !( status->is_counter() ) ) continue; // ?
326
327 double tofPH = ( *iter_tof )->ph();
328 double tof = ( *iter_tof )->tof();
329
330 h_tofElecIP_Barrel->fill( tof - ( *iter_tof )->texpElectron() );
331 h_tofMuonIP_Barrel->fill( tof - ( *iter_tof )->texpMuon() );
332 h_tofPionIP_Barrel->fill( tof - ( *iter_tof )->texpPion() );
333 h_tofKaonIP_Barrel->fill( tof - ( *iter_tof )->texpKaon() );
334 h_tofProtonIP_Barrel->fill( tof - ( *iter_tof )->texpProton() );
335 h_tofVsMomentumIP->fill( pMag, tof );
336
337 if ( status->layer() == 1 )
338 {
339 h_tofPHIP_BarrelLayer1->fill( tofPH );
340 h_tofIP_BarrelLayer1->fill( tof );
341 }
342 if ( status->layer() == 2 )
343 {
344 h_tofPHIP_BarrelLayer2->fill( tofPH );
345 h_tofIP_BarrelLayer2->fill( tof );
346 }
347 } // End of TOF barrel IF
348
349 else
350 {
351 if ( !( status->is_counter() ) ) continue; // ?
352
353 double tofPH = ( *iter_tof )->ph();
354 double tof = ( *iter_tof )->tof();
355
356 h_tofElecIP_Endcap->fill( tof - ( *iter_tof )->texpElectron() );
357 h_tofMuonIP_Endcap->fill( tof - ( *iter_tof )->texpMuon() );
358 h_tofPionIP_Endcap->fill( tof - ( *iter_tof )->texpPion() );
359 h_tofKaonIP_Endcap->fill( tof - ( *iter_tof )->texpKaon() );
360 h_tofProtonIP_Endcap->fill( tof - ( *iter_tof )->texpProton() );
361 h_tofVsMomentumIP->fill( pMag, tof );
362
363 if ( status->is_east() )
364 {
365 h_tofPHIP_EastEndcap->fill( tofPH );
366 h_tofIP_EastEndcap->fill( tof );
367 }
368 else
369 {
370 h_tofPHIP_WestEndcap->fill( tofPH );
371 h_tofIP_WestEndcap->fill( tof );
372 }
373 } // End of TOF endcap IF
374
375 } // End of "iter_tof" FOR
376 } // End of "isTofTrackValid()" IF
377
378 /// For the 2 highest momentum charged tracks
379 if ( pMag > highestIPTrackP )
380 {
381 secondHighestPIPTrackId = highestPIPTrackId;
382 secondHighestIPTrackP = highestIPTrackP;
383 highestPIPTrackId = trackId;
384 highestIPTrackP = pMag;
385 }
386 if ( ( pMag > secondHighestIPTrackP ) && ( pMag < highestIPTrackP ) )
387 {
388 secondHighestPIPTrackId = trackId;
389 secondHighestIPTrackP = pMag;
390 }
391
392 } // End of charged track FOR
393
394 ///////////////////////////////////////
395 /// If the event has 1 IP charged track
396 if ( nChargedTracksIP > 0 )
397 {
398 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + highestPIPTrackId;
399 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
400 double highestPPhi0 = mdcTrk->phi();
401 h_pIPTrack1DivEb->fill( mdcTrk->p() / eBeam );
402
403 if ( ( *itTrk )->isEmcShowerValid() )
404 {
405 RecEmcShower* emcChargedTrk = ( *itTrk )->emcShower();
406 h_eEMCIPTrack1DivEb->fill( emcChargedTrk->energy() / eBeam );
407 }
408
409 /// If the event has 2 IP charged tracks
410 if ( nChargedTracksIP > 1 )
411 {
412 itTrk = evtRecTrkCol->begin() + secondHighestPIPTrackId;
413 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
414 double secondHighestPPhi0 = mdcTrk->phi();
415 h_pIPTrack2DivEb->fill( mdcTrk->p() / eBeam );
416
417 if ( ( *itTrk )->isEmcShowerValid() )
418 {
419 RecEmcShower* emcChargedTrk = ( *itTrk )->emcShower();
420 h_eEMCIPTrack2DivEb->fill( emcChargedTrk->energy() / eBeam );
421 }
422
423 h_acoplanarity_2HighestPIPTracks->fill(
424 fabs( CLHEP::pi - fabs( highestPPhi0 - secondHighestPPhi0 ) ) );
425 } // End of "nChargedTracksIP > 1" IF
426 } // End of "nChargedTracksIP > 0" IF
427
428 ///////////////////////////////////////
429 /// Shower (aka Neutral track) FOR loop
430 int nNeutralTracks = 0, nNeutralTracksGT30MeV = 0;
431 double totalNeutralEnergy = 0;
432 double totalNeutralPX = 0, totalNeutralPY = 0, totalNeutralPZ = 0;
433
434 for ( int i = evtRecEvent->totalCharged(); i < evtRecEvent->totalTracks(); i++ )
435 {
436 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
437 if ( !( *itTrk )->isEmcShowerValid() ) continue;
438 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
439
440 nNeutralTracks++;
441 double eraw = emcTrk->energy();
442 if ( eraw > 0.030 ) nNeutralTracksGT30MeV++;
443
444 totalVisibleEnergy += eraw;
445 totalEMCEnergy += eraw;
446 totalNeutralEnergy += eraw;
447
448 double theta = emcTrk->theta();
449 double phi = emcTrk->phi();
450
451 double pX = eraw * cos( phi ) * sin( theta );
452 double pY = eraw * sin( phi ) * sin( theta );
453 double pZ = eraw * cos( theta );
454
455 totalNeutralPX += pX;
456 totalNeutralPY += pY;
457 totalNeutralPZ += pZ;
458
459 /// Energy of most energetic showers
460 if ( nNeutralTracks == 1 ) h_eEMCShower1DivEb->fill( eraw / eBeam );
461 if ( nNeutralTracks == 2 ) h_eEMCShower2DivEb->fill( eraw / eBeam );
462 if ( nNeutralTracks == 3 ) h_eEMCShower3DivEb->fill( eraw / eBeam );
463
464 } // End of neutral track FOR
465
466 ///////////////////////////////
467 /// Histograms filled per event
468
469 h_nChargedTracks->fill( nChargedTracks );
470 h_nChargedTracksIP->fill( nChargedTracksIP );
471
472 h_netCharge->fill( nCharge );
473 h_netChargeIP->fill( nChargeIP );
474
475 h_nNeutralTracks->fill( nNeutralTracks );
476 h_nNeutralTracksGT30MeV->fill( nNeutralTracksGT30MeV );
477
478 /// Total Energy Histograms
479 h_eVisibleDivEcm->fill( totalVisibleEnergy / m_ecm );
480 h_eEMCDivEcm->fill( totalEMCEnergy / m_ecm );
481 h_eNeutralDivEcm->fill( totalNeutralEnergy / m_ecm );
482 h_eChargedDivEcm->fill( totalChargedEnergy / m_ecm );
483
484 /// Total Charged Momentum Histograms
485 double totalChargedPXY =
486 sqrt( totalChargedPX * totalChargedPX + totalChargedPY * totalChargedPY );
487 double totalChargedP =
488 sqrt( totalChargedPXY * totalChargedPXY + totalChargedPZ * totalChargedPZ );
489
490 h_netMomentumDivEcm_AllChargedTracks->fill( totalChargedP / m_ecm );
491 h_netTransMomentumDivEcm_AllChargedTracks->fill( totalChargedPXY / m_ecm );
492 if ( totalChargedP > 0 )
493 { h_cosTheta_AllChargedTracks->fill( totalChargedPZ / totalChargedP ); }
494 else
495 {
496 if ( nChargedTracks > 0 )
497 {
498 log << MSG::INFO << "Run " << run << ", event " << event << ": totalChargedP <= 0! "
499 << endmsg;
500 }
501 }
502
503 /// Total Neutral Momentum Histograms
504 double totalNeutralPXY =
505 sqrt( totalNeutralPX * totalNeutralPX + totalNeutralPY * totalNeutralPY );
506 double totalNeutralP =
507 sqrt( totalNeutralPXY * totalNeutralPXY + totalNeutralPZ * totalNeutralPZ );
508
509 h_netMomentumDivEcm_AllNeutralTracks->fill( totalNeutralP / m_ecm );
510 h_netTransMomentumDivEcm_AllNeutralTracks->fill( totalNeutralPXY / m_ecm );
511 if ( totalNeutralP > 0 )
512 { h_cosTheta_AllNeutralTracks->fill( totalNeutralPZ / totalNeutralP ); }
513 else
514 {
515 if ( nNeutralTracks > 0 )
516 {
517 log << MSG::INFO << "Run " << run << ", event " << event << ": totalNeutralP <= 0! "
518 << endmsg;
519 }
520 }
521
522 /// Total Momentum Histograms
523 double totalEventPX = totalChargedPX + totalNeutralPX;
524 double totalEventPY = totalChargedPY + totalNeutralPY;
525 double totalEventPZ = totalChargedPZ + totalNeutralPZ;
526
527 double totalEventPXY = sqrt( totalEventPX * totalEventPX + totalEventPY * totalEventPY );
528 double totalEventP = sqrt( totalEventPXY * totalEventPXY + totalEventPZ * totalEventPZ );
529
530 h_netMomentumDivEcm_AllTracks->fill( totalEventP / m_ecm );
531 h_netTransMomentumDivEcm_AllTracks->fill( totalEventPXY / m_ecm );
532 if ( totalEventP > 0 ) { h_cosTheta_AllTracks->fill( totalEventPZ / totalEventP ); }
533 else
534 {
535 log << MSG::INFO << "Run " << run << ", event " << event << ": totalEventP <= 0! "
536 << endmsg;
537 }
538
539 /// VeeVertex information
540 SmartDataPtr<EvtRecVeeVertexCol> evtRecVeeVertexCol( eventSvc(),
541 "/Event/EvtRec/EvtRecVeeVertexCol" );
542 if ( !evtRecVeeVertexCol )
543 {
544 log << MSG::FATAL << "Could not find EvtRecVeeVertexCol" << endmsg;
545 return StatusCode::FAILURE;
546 }
547
548 /// Loop over VeeVertex candidates
549 int numKs = 0, numLambda = 0;
550 for ( EvtRecVeeVertexCol::iterator veeItr = evtRecVeeVertexCol->begin();
551 veeItr < evtRecVeeVertexCol->end(); veeItr++ )
552 {
553
554 h_ksMass->fill( ( *veeItr )->mass() );
555 if ( fabs( ( *veeItr )->mass() - mKs ) < 0.010 ) ++numKs;
556
557 h_lambdaMass->fill( ( *veeItr )->mass() );
558 if ( fabs( ( *veeItr )->mass() - mLambda ) < 0.010 ) ++numLambda;
559
560 } // End of "evtRecVeeVertexCol" FOR
561
562 h_nKs->fill( numKs );
563 h_nLambda->fill( numLambda );
564
565 return StatusCode::SUCCESS;
566}
567
568// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
569StatusCode FarmMonitorAlg::finalize() {
570
571 MsgStream log( msgSvc(), name() );
572 log << MSG::INFO << "in finalize()" << endmsg;
573
574 return StatusCode::SUCCESS;
575}
DECLARE_COMPONENT
DECLARE_COMPONENT(BesBdkRc)
EvtRecTrackIterator
EvtRecTrackCol::iterator EvtRecTrackIterator
Definition Event/EvtRecEvent/include/EvtRecEvent/EvtRecTrack.h:106
mKs
double mKs
Definition FarmMonitorAlg.cxx:30
mPi
double mPi
Definition FarmMonitorAlg.cxx:29
mLambda
double mLambda
Definition FarmMonitorAlg.cxx:31
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
histoSvc
IHistogramSvc * histoSvc()
Definition ServiceAccessor.h:23
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:12
FarmMonitorAlg::finalize
StatusCode finalize()
Definition FarmMonitorAlg.cxx:569
FarmMonitorAlg::execute
StatusCode execute()
Definition FarmMonitorAlg.cxx:197
FarmMonitorAlg::FarmMonitorAlg
FarmMonitorAlg(const std::string &name, ISvcLocator *pSvcLocator)
Definition FarmMonitorAlg.cxx:35
FarmMonitorAlg::initialize
StatusCode initialize()
Definition FarmMonitorAlg.cxx:47
TofHitStatus::setStatus
void setStatus(unsigned int status)
Definition TofHitStatus.cxx:115