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