BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
HistSample.cxx
Go to the documentation of this file.
1#include <fstream>
2
3#include <math.h>
4// #include <string.h>
5
6#include "HistSample.h"
7// #include "GeneratorObject/McEventCollection.h"
8
9#include "GaudiKernel/MsgStream.h"
10
11#include "GaudiKernel/DataSvc.h"
12#include "GaudiKernel/SmartDataPtr.h"
13
14#include "GaudiKernel/IHistogramSvc.h"
15
16#include "AIDA/IHistogram1D.h"
17#include "AIDA/IHistogram2D.h"
18
19#include "GaudiKernel/IDataProviderSvc.h"
20#include "GaudiKernel/INTupleSvc.h"
21#include "GaudiKernel/ISvcLocator.h"
22#include "GaudiKernel/MsgStream.h"
23#include "GaudiKernel/ObjectList.h"
24#include "GaudiKernel/PropertyMgr.h"
25
26// #include "StoreGate/StoreGateSvc.h"
27
28#include "HepMC/GenEvent.h"
29#include "HepMC/GenParticle.h"
30#include "HepMC/GenVertex.h"
31// #include "CLHEP/HepMC/GenEvent.h"
32// #include "CLHEP/HepMC/GenParticle.h"
33// #include "CLHEP/HepMC/GenVertex.h"
34// #include "CLHEP/HepMC/ParticleDataTableConfig.h"
35
36// #include "PartPropSvc/PartPropSvc.h"
37#include "HepPDT/ParticleData.hh"
38
39#include "GeneratorModule/GeneratorName.h"
40
41#
42// static const AlgFactory<HistSample> Factory;
43// const IAlgFactory& HistSampleFactory = Factory;
44
45HistSample::HistSample( const std::string& name, ISvcLocator* pSvcLocator )
46 : Algorithm( name, pSvcLocator )
47 , m_hgenerated( 0 )
48 , m_hfinal( 0 )
49 , m_ncharged( 0 )
50 , m_hChargedPt( 0 )
51 , m_hChargedEta( 0 )
52 , m_hZPtall( 0 )
53 , m_hZPt( 0 )
54 , m_hZPte( 0 )
55 , m_hZPtm( 0 )
56 , m_hZPtt( 0 )
57 , m_massZall( 0 )
58 , m_massZ( 0 )
59 , m_massZe( 0 )
60 , m_massZm( 0 )
61 , m_massZt( 0 )
62 , m_hPtPaire( 0 )
63 , m_hPtPairm( 0 )
64 , m_hPtPairt( 0 )
65 , m_massPaire( 0 )
66 , m_massPairm( 0 )
67 , m_massPairt( 0 )
68 , m_rapidity( 0 )
69 , m_pseudorapidity( 0 )
70 , m_hpte() {
71 // Declare the algorithm's properties
72 declareProperty( "HistogramFlag", m_produceHistogram = true );
73 // declareProperty("HistogramFlag", m_produceHistogram = false );
74}
75
76StatusCode HistSample::initialize() {
77
78 StatusCode result = StatusCode::SUCCESS;
79 MsgStream msglog( msgSvc(), name() );
80 msglog << MSG::INFO << ">>> HistSample from Initialize" << endmsg;
81 // cout << "----- HistSample World From initialize" << endl;
82
83 /*
84 StatusCode sc = service("StoreGateSvc", m_sgSvc);
85 if (sc.isFailure()) {
86 msglog << MSG::ERROR << "Could not find StoreGateSvc" << endmsg;
87 return sc;
88 }
89 */
90
91 /*
92 // Get the Particle Properties Service
93 IPartPropSvc* p_PartPropSvc;
94 static const bool CREATEIFNOTTHERE(true);
95 StatusCode PartPropStatus = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
96 if (!PartPropStatus.isSuccess() || 0 == p_PartPropSvc) {
97 msglog << MSG::ERROR << " Could not initialize Particle Properties Service" << endmsg;
98 return PartPropStatus;
99 }
100
101 m_particleTable = p_PartPropSvc->PDT();
102 */
103
104 // Register the histograms
105 // m_hgenerated = histoSvc()->book("/stat/1Dhist/1","Generated",100,0,1000);
106 // if( 0 == m_hgenerated ) {
107 // msglog << MSG::ERROR << "Cannot register histo Generated" << endmsg;
108 // }
109
110 m_hgenerated = histoSvc()->book( "/stat/1Dhist/1", "Generated", 100, 0, 1200 );
111 if ( 0 == m_hgenerated )
112 {
113 msglog << MSG::ERROR << " ERROR booking histogram" << endmsg;
114 result = StatusCode::FAILURE;
115 }
116 m_hfinal = histoSvc()->book( "/stat/1Dhist/2", "Final state", 100, 0, 800 );
117 m_ncharged = histoSvc()->book( "/stat/1Dhist/3", "number of charged tracks", 100, 0, 500 );
118 m_hChargedPt = histoSvc()->book( "/stat/1Dhist/4", "Pt of charged tracks", 100, 0, 100 );
119 m_hChargedEta =
120 histoSvc()->book( "/stat/1Dhist/5", "Pseudorapidity of charged tracks", 100, -12, 12 );
121 m_hZPtall = histoSvc()->book( "/stat/1Dhist/6", "ZPt, all", 100, 0, 200 );
122 m_hZPt = histoSvc()->book( "/stat/1Dhist/7", "ZPt, charged leptons", 100, 0, 200 );
123 m_hZPte = histoSvc()->book( "/stat/1Dhist/8", "ZPt electrons", 100, 0, 200 );
124 m_hZPtm = histoSvc()->book( "/stat/1Dhist/9", "ZPt muons", 100, 0, 200 );
125 m_hZPtt = histoSvc()->book( "/stat/1Dhist/10", "ZPt taus", 100, 0, 200 );
126 m_massZall = histoSvc()->book( "/stat/1Dhist/11", "Z mass, all", 40, 70, 110 );
127 m_massZ = histoSvc()->book( "/stat/1Dhist/12", "Z mass, charged leptons", 40, 70, 110 );
128 m_massZe = histoSvc()->book( "/stat/1Dhist/13", "Z mass, electrons", 40, 70, 110 );
129 m_massZm = histoSvc()->book( "/stat/1Dhist/14", "Z mass, muons", 40, 70, 110 );
130 m_massZt = histoSvc()->book( "/stat/1Dhist/15", "Z mass, taus", 40, 70, 110 );
131 m_hPtPaire = histoSvc()->book( "/stat/1Dhist/16", "Pt electron pairs", 100, 0, 200 );
132 m_hPtPairm = histoSvc()->book( "/stat/1Dhist/17", "Pt muon pairs", 100, 0, 200 );
133 m_hPtPairt = histoSvc()->book( "/stat/1Dhist/18", "Pt tau pairs", 100, 0, 200 );
134 m_massPaire = histoSvc()->book( "/stat/1Dhist/19", "mass, electron pairs", 40, 70, 110 );
135 m_massPairm = histoSvc()->book( "/stat/1Dhist/20", "mass, muon pairs", 40, 70, 110 );
136 m_massPairt = histoSvc()->book( "/stat/1Dhist/21", "mass, tau pairs", 40, 70, 110 );
137 m_rapidity =
138 histoSvc()->book( "/stat/1Dhist/22", "Rapidity of charged tracks", 100, -12, 12 );
139 m_pseudorapidity =
140 histoSvc()->book( "/stat/1Dhist/23", "Pseudorapidity of charged tracks", 100, -12, 12 );
141 m_hpte = histoSvc()->book( "/stat/1Dhist/24", "electon pt", 100, 0, 100 );
142
143 // Initialization terminated
144 // return StatusCode::SUCCESS;
145 return result;
146}
147StatusCode HistSample::execute() {
148
149 MsgStream msglog( msgSvc(), name() );
150 msglog << MSG::INFO << ">>> HistSample from execute" << endmsg;
151 // cout << "----- HistSample World From execute" << endl;
152
153 // Read Data from Transient Store
154 /*
155 const McEventCollection* mcCollptr;
156 if ( m_sgSvc->retrieve(mcCollptr).isFailure() ) {
157 msglog << MSG::ERROR << "Could not retrieve McEventCollection"
158 << endmsg;
159 return StatusCode::FAILURE;
160 }
161
162 // cout << " ---- Begin Iterating Over McEventCollection --- " << endl;
163 McEventCollection::const_iterator it;
164 for(it=mcCollptr->begin(); it!=mcCollptr->end(); it++) {
165 // cout << " Generator: " << (*it)->generatorName() << endl;
166 //// Dump on screen
167 // (*it)->pGenEvt()->print();
168
169 // fix the STDHEP mess for the Z status
170 int properStatus = 2;
171 // switch ( (*it)->generatorName()[0] ) {
172 // case 'P':
173 // properStatus = 2;
174 // break;
175 // case 'I':
176 // properStatus = 12;
177 // break;
178 // case 'H':
179 // properStatus = 195;
180 // break;
181 // default:
182 // properStatus = 2;
183 // break;
184 // }
185 HepMC::GenEvent* genEvt = (*it);
186 int g_id = genEvt->signal_process_id();
187 switch ( first_generator(g_id) ) {
188 case PYTHIA:
189 properStatus = 2;
190 break;
191 case ISAJET:
192 properStatus = 12;
193 break;
194 case HERWIG:
195 properStatus = 195;
196 break;
197 default:
198 properStatus = 2;
199 break;
200 }
201
202
203 int number_particles = 0;
204 int final_state = 0;
205 int number_charged = 0;
206 // Iterate over MC particles
207 for(HepMC::GenEvent::particle_iterator pitr=genEvt->particles_begin();
208 pitr!=genEvt->particles_end(); ++pitr ){
209 ++number_particles;
210 // cout << "particle " << ((*pitr)->pdg_id()) << " status " <<
211 ((*pitr)->status()) << endl;
212
213 // Z decays
214 if( ((*pitr)->pdg_id()) == 23 && ((*pitr)->status()) == properStatus ){
215 HepMC::GenVertex::particle_iterator firstChild =
216 (*pitr)->end_vertex()->particles_begin(HepMC::children);
217 HepMC::GenVertex::particle_iterator endChild =
218 (*pitr)->end_vertex()->particles_end(HepMC::children);
219 // plot Pt and mass of the Z (generator values)
220 if( ((*firstChild)->pdg_id()) != 23 ){
221 double ZPt = (*pitr)->momentum().perp();
222 m_hZPtall->fill( ZPt, 1.);
223 double Zmass = (*pitr)->momentum().m();
224 m_massZall->fill(Zmass, 1.);
225 }
226 // Z decays to l+l-
227 double sumx = 0.0;
228 double sumy = 0.0;
229 double sumz = 0.0;
230 double sume = 0.0;
231 int nelds = 0;
232 int nmuds = 0;
233 int ntauds = 0;
234 int Zchild = 0;
235 HepMC::GenVertex::particle_iterator thisChild = firstChild;
236 for(; thisChild != endChild++; ++thisChild){
237 if( ((*thisChild)->pdg_id()) != 23 ){
238 // cout << "Zchild id " << (*thisChild)->pdg_id() << endl;
239 ++Zchild;
240 if( abs((*thisChild)->pdg_id()) == 11 || abs((*thisChild)->pdg_id()) ==
241 13 || abs((*thisChild)->pdg_id()) == 15 ){ if( abs((*thisChild)->pdg_id()) == 11 )++nelds;
242 if( abs((*thisChild)->pdg_id()) == 13 )++nmuds;
243 if( abs((*thisChild)->pdg_id()) == 15 )++ntauds;
244 sumx += (*thisChild)->momentum().x();
245 sumy += (*thisChild)->momentum().y();
246 sumz += (*thisChild)->momentum().z();
247 sume += (*thisChild)->momentum().e();
248 }
249 }
250 }
251 // cout << "Zchildren " << Zchild << " nelds " << nelds << " nmuds " << nmuds
252 << " ntauds " << ntauds << endl; if( Zchild != 0 ){ double PtZ2 = sumx*sumx + sumy*sumy;
253 double PtZ = 0.;
254 if(PtZ2 > 0.) PtZ = sqrt(PtZ2);
255 if(PtZ != 0.)m_hZPt->fill( PtZ, 1.);
256 double massZ2 = sume*sume - sumx*sumx - sumy*sumy - sumz*sumz;
257 double massZ = 0.;
258 if(massZ2 > 0.) massZ = sqrt(massZ2);
259 if(massZ != 0.)m_massZ->fill( massZ, 1.);
260 if(nelds == 2){
261 m_hZPte->fill( PtZ, 1.);
262 m_massZe->fill( massZ, 1.);
263 }
264 if(nmuds == 2){
265 m_hZPtm->fill( PtZ, 1.);
266 m_massZm->fill( massZ, 1.);
267 }
268 if(ntauds == 2){
269 m_hZPtt->fill( PtZ, 1.);
270 m_massZt->fill( massZ, 1.);
271 }
272 }
273 }
274 // all decays to l+l- pairs in the event
275 if( (*pitr)->end_vertex() ){
276 HepMC::GenVertex::particle_iterator fstChild =
277 (*pitr)->end_vertex()->particles_begin(HepMC::children);
278 HepMC::GenVertex::particle_iterator lstChild =
279 (*pitr)->end_vertex()->particles_end(HepMC::children);
280 double sumpx = 0.0;
281 double sumpy = 0.0;
282 double sumpz = 0.0;
283 double sumpe = 0.0;
284 int nel = 0;
285 int nmu = 0;
286 int ntau = 0;
287 HepMC::GenVertex::particle_iterator aChild = fstChild;
288 for(; aChild != lstChild++; ++aChild){
289 if( abs((*aChild)->pdg_id()) == 11 || abs((*aChild)->pdg_id()) == 13 ||
290 abs((*aChild)->pdg_id()) == 15 ){
291 if( abs((*aChild)->pdg_id()) == 11 )++nel;
292 if( abs((*aChild)->pdg_id()) == 13 )++nmu;
293 if( abs((*aChild)->pdg_id()) == 15 )++ntau;
294 sumpx += (*aChild)->momentum().x();
295 sumpy += (*aChild)->momentum().y();
296 sumpz += (*aChild)->momentum().z();
297 sumpe += (*aChild)->momentum().e();
298 }
299 }
300 if(nel == 2 || nmu == 2 || ntau == 2 ){
301 double PtPair2 = sumpx*sumpx + sumpy*sumpy;
302 double PtPair = 0.;
303 if(PtPair2 > 0.) PtPair = sqrt(PtPair2);
304 double massPair2 = sumpe*sumpe - sumpx*sumpx - sumpy*sumpy - sumpz*sumpz;
305 double massPair = 0.;
306 if(massPair2 > 0.) massPair = sqrt(massPair2);
307 if(nel == 2){
308 m_hPtPaire->fill( PtPair, 1.);
309 m_massPaire->fill( massPair, 1.);
310 }
311 if(nmu == 2){
312 m_hPtPairm->fill( PtPair, 1.);
313 m_massPairm->fill( massPair, 1.);
314 }
315 if(ntau == 2){
316 m_hPtPairt->fill( PtPair, 1.);
317 m_massPairt->fill( massPair, 1.);
318 }
319 }
320 }
321 // charged tracks
322 double c = 0.;
323 HepPDT::ParticleData* ap = m_particleTable->particle( abs( (*pitr)->pdg_id() )
324 );
325
326 if(!ap){
327 // std::cout << "ChargeService WARNING: abs(id) "
328 // << abs((*pitr)->pdg_id())
329 // << " is not in particle data table" << std::endl;
330 }
331 else
332 {
333 c = ap->charge();
334 }
335
336 if( ((*pitr)->status() == 1) && ( ! (*pitr)->end_vertex()) ) {
337 ++final_state;
338 if( ((*pitr)->pdg_id()) == 11 ){
339 double ePt = (*pitr)->momentum().perp();
340 m_hpte->fill( ePt, 1.);
341 }
342 if(c!=0.){
343 ++number_charged;
344 double chargedPt = (*pitr)->momentum().perp();
345 double chargedEta = (*pitr)->momentum().pseudoRapidity();
346 m_hChargedPt->fill( chargedPt, 1.);
347 m_hChargedEta->fill( chargedEta, 1.);
348
349 double px = (*pitr)->momentum().x();
350 double py = (*pitr)->momentum().y();
351 double pz = (*pitr)->momentum().z();
352 double pe = (*pitr)->momentum().e();
353 double pp2 = px*px + py*py + pz*pz;
354 double pp3 = 0.;
355 if(pp2 > 0.) pp3 = sqrt(pp2);
356 double y = -999.;
357 if( (pe-pz) != 0. && (pe+pz)/(pe-pz) > 0. ) y =
358 (1./2.)*log((pe+pz)/(pe-pz)); double eta = -999.; if( (pp3-pz) != 0. && (pp3+pz)/(pp3-pz) >
359 0. ) eta = (1./2.)*log((pp3+pz)/(pp3-pz)); m_rapidity->fill( y, 1.); m_pseudorapidity->fill(
360 eta, 1.);
361 }
362 }
363 }
364 m_hgenerated->fill( number_particles, 1.);
365 m_hfinal->fill( final_state, 1.);
366 m_ncharged->fill( number_charged, 1.);
367
368 }
369 */
370 // End of execution for each event
371 return StatusCode::SUCCESS;
372}
373StatusCode HistSample::finalize() {
374
375 MsgStream msglog( msgSvc(), name() );
376 msglog << MSG::INFO << ">>> HistSample from finalize" << endmsg;
377 // cout << "----- HistSample World From finalize" << endl;
378
379 // End of finalization step
380 return StatusCode::SUCCESS;
381}
histoSvc
IHistogramSvc * histoSvc()
Definition ServiceAccessor.h:23
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:12
HistSample::execute
StatusCode execute()
Definition HistSample.cxx:147
HistSample::HistSample
HistSample(const std::string &name, ISvcLocator *pSvcLocator)
Definition HistSample.cxx:45
HistSample::finalize
StatusCode finalize()
Definition HistSample.cxx:373
HistSample::initialize
StatusCode initialize()
Definition HistSample.cxx:76