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