HepMC3 event record library
testIO4.cc
1 // -*- C++ -*-
2 //
3 // This file is part of HepMC
4 // Copyright (C) 2014-2019 The HepMC collaboration (see AUTHORS for details)
5 //
6 #include "HepMC3/GenEvent.h"
7 #include "HepMC3/GenParticle.h"
12 #include "HepMC3TestUtils.h"
13 #include <TChain.h>
14 #include <TFile.h>
15 #include <TROOT.h>
16 #include <TH1D.h>
17 using namespace HepMC3;
18 const Int_t kMaxparticles = 2000;
19 const Int_t kMaxvertices = 2000;
20 #ifndef DOXYGEN_SHOULD_SKIP_THIS
21 class SomeAnalysis
22 {
23 public :
24  TChain *fChain; //!pointer to the analyzed TTree or TChain
25  // Declaration of leaf types
26 //GenEventData *hepmc3_event;
27  Int_t event_number;
28  Int_t momentum_unit;
29  Int_t length_unit;
30  Int_t particles_;
31  Int_t particles_pid[kMaxparticles]; //[particles_]
32  Int_t particles_status[kMaxparticles]; //[particles_]
33  Bool_t particles_is_mass_set[kMaxparticles]; //[particles_]
34  Double_t particles_mass[kMaxparticles]; //[particles_]
35  Double_t particles_momentum_m_v1[kMaxparticles]; //[particles_]
36  Double_t particles_momentum_m_v2[kMaxparticles]; //[particles_]
37  Double_t particles_momentum_m_v3[kMaxparticles]; //[particles_]
38  Double_t particles_momentum_m_v4[kMaxparticles]; //[particles_]
39  Int_t vertices_;
40  Int_t vertices_status[kMaxvertices]; //[vertices_]
41  Double_t vertices_position_m_v1[kMaxvertices]; //[vertices_]
42  Double_t vertices_position_m_v2[kMaxvertices]; //[vertices_]
43  Double_t vertices_position_m_v3[kMaxvertices]; //[vertices_]
44  Double_t vertices_position_m_v4[kMaxvertices]; //[vertices_]
45  vector<double> weights;
46  Double_t event_pos_m_v1;
47  Double_t event_pos_m_v2;
48  Double_t event_pos_m_v3;
49  Double_t event_pos_m_v4;
50  vector<int> links1;
51  vector<int> links2;
52  vector<int> attribute_id;
53  vector<string> attribute_name;
54  vector<string> attribute_string;
55 
56  // List of branches
57  TBranch *b_hepmc3_event_event_number; //!
58  TBranch *b_hepmc3_event_momentum_unit; //!
59  TBranch *b_hepmc3_event_length_unit; //!
60  TBranch *b_hepmc3_event_particles_; //!
61  TBranch *b_particles_pid; //!
62  TBranch *b_particles_status; //!
63  TBranch *b_particles_is_mass_set; //!
64  TBranch *b_particles_mass; //!
65  TBranch *b_particles_momentum_m_v1; //!
66  TBranch *b_particles_momentum_m_v2; //!
67  TBranch *b_particles_momentum_m_v3; //!
68  TBranch *b_particles_momentum_m_v4; //!
69  TBranch *b_hepmc3_event_vertices_; //!
70  TBranch *b_vertices_status; //!
71  TBranch *b_vertices_position_m_v1; //!
72  TBranch *b_vertices_position_m_v2; //!
73  TBranch *b_vertices_position_m_v3; //!
74  TBranch *b_vertices_position_m_v4; //!
75  TBranch *b_hepmc3_event_weights; //!
76  TBranch *b_hepmc3_event_event_pos_m_v1; //!
77  TBranch *b_hepmc3_event_event_pos_m_v2; //!
78  TBranch *b_hepmc3_event_event_pos_m_v3; //!
79  TBranch *b_hepmc3_event_event_pos_m_v4; //!
80  TBranch *b_hepmc3_event_links1; //!
81  TBranch *b_hepmc3_event_links2; //!
82  TBranch *b_hepmc3_event_attribute_id; //!
83  TBranch *b_hepmc3_event_attribute_name; //!
84  TBranch *b_hepmc3_event_attribute_string; //!
85 
86 
87  void Init(TChain *tree)
88  {
89  if (!tree) return;
90  fChain = tree;
91  fChain->SetMakeClass(1);
92 
93  fChain->SetBranchAddress("event_number", &event_number, &b_hepmc3_event_event_number);
94  fChain->SetBranchAddress("momentum_unit", &momentum_unit, &b_hepmc3_event_momentum_unit);
95  fChain->SetBranchAddress("length_unit", &length_unit, &b_hepmc3_event_length_unit);
96  fChain->SetBranchAddress("particles", &particles_, &b_hepmc3_event_particles_);
97  fChain->SetBranchAddress("particles.pid", particles_pid, &b_particles_pid);
98  fChain->SetBranchAddress("particles.status", particles_status, &b_particles_status);
99  fChain->SetBranchAddress("particles.is_mass_set", particles_is_mass_set, &b_particles_is_mass_set);
100  fChain->SetBranchAddress("particles.mass", particles_mass, &b_particles_mass);
101  fChain->SetBranchAddress("particles.momentum.m_v1", particles_momentum_m_v1, &b_particles_momentum_m_v1);
102  fChain->SetBranchAddress("particles.momentum.m_v2", particles_momentum_m_v2, &b_particles_momentum_m_v2);
103  fChain->SetBranchAddress("particles.momentum.m_v3", particles_momentum_m_v3, &b_particles_momentum_m_v3);
104  fChain->SetBranchAddress("particles.momentum.m_v4", particles_momentum_m_v4, &b_particles_momentum_m_v4);
105  fChain->SetBranchAddress("vertices", &vertices_, &b_hepmc3_event_vertices_);
106  fChain->SetBranchAddress("vertices.status", vertices_status, &b_vertices_status);
107  fChain->SetBranchAddress("vertices.position.m_v1", vertices_position_m_v1, &b_vertices_position_m_v1);
108  fChain->SetBranchAddress("vertices.position.m_v2", vertices_position_m_v2, &b_vertices_position_m_v2);
109  fChain->SetBranchAddress("vertices.position.m_v3", vertices_position_m_v3, &b_vertices_position_m_v3);
110  fChain->SetBranchAddress("vertices.position.m_v4", vertices_position_m_v4, &b_vertices_position_m_v4);
111  fChain->SetBranchAddress("weights", &weights, &b_hepmc3_event_weights);
112  fChain->SetBranchAddress("event_pos.m_v1", &event_pos_m_v1, &b_hepmc3_event_event_pos_m_v1);
113  fChain->SetBranchAddress("event_pos.m_v2", &event_pos_m_v2, &b_hepmc3_event_event_pos_m_v2);
114  fChain->SetBranchAddress("event_pos.m_v3", &event_pos_m_v3, &b_hepmc3_event_event_pos_m_v3);
115  fChain->SetBranchAddress("event_pos.m_v4", &event_pos_m_v4, &b_hepmc3_event_event_pos_m_v4);
116  fChain->SetBranchAddress("links1", &links1, &b_hepmc3_event_links1);
117  fChain->SetBranchAddress("links2", &links2, &b_hepmc3_event_links2);
118  fChain->SetBranchAddress("attribute_id", &attribute_id, &b_hepmc3_event_attribute_id);
119  fChain->SetBranchAddress("attribute_name", &attribute_name, &b_hepmc3_event_attribute_name);
120  fChain->SetBranchAddress("attribute_string", &attribute_string, &b_hepmc3_event_attribute_string);
121 
122  fChain->SetBranchStatus("*",0);
123  fChain->SetBranchStatus("particles",1);
124  fChain->SetBranchStatus("particles.pid",1);
125  fChain->SetBranchStatus("particles.status",1);
126  fChain->SetBranchStatus("particles.momentum.m_v1",1);
127  fChain->SetBranchStatus("particles.momentum.m_v2",1);
128 
129  }
130  SomeAnalysis(const std::string& file)
131  {
132  TChain* TempChain= new TChain("hepmc3_tree");
133  TempChain->Add(file.c_str());
134  Init(TempChain);
135  }
136 };
137 #endif
138 
139 int main()
140 {
141 //Plain tree
142  TH1D H1("H1","Pt of pions;Events/100MeV;P_{T},GeV",1000,0,100);
143  SomeAnalysis* A= new SomeAnalysis("inputIO4.root");
144  if (!A->fChain->GetEntries()) return 10001;
145  for (int entry=0; entry<A->fChain->GetEntries(); entry++)
146  {
147  A->fChain->GetEntry(entry);
148  for (int i=0; i<A->particles_; i++)
149  if (A->particles_status[i]==1&&(std::abs(A->particles_pid[i])==211||std::abs(A->particles_pid[i])==11))
150  H1.Fill(std::sqrt(A->particles_momentum_m_v1[i]*A->particles_momentum_m_v1[i]+A->particles_momentum_m_v2[i]*A->particles_momentum_m_v2[i]) );
151  }
152  delete A;
153 //GenEvent
154  TH1D H2("H2","Pt of pions;Events/100MeV;P_{T},GeV",1000,0,100);
155  ReaderRootTree inputA("inputIO4.root");
156  if(inputA.failed()) return 10002;
157  while( !inputA.failed() )
158  {
159  GenEvent evt(Units::GEV,Units::MM);
160  inputA.read_event(evt);
161  if( inputA.failed() ) {
162  printf("End of file reached. Exit.\n");
163  break;
164  }
165  for (ConstGenParticlePtr p: evt.particles())
166  if ( std::abs(p->status()) == 1 && (std::abs(p->pdg_id()) == 211||std::abs(p->pdg_id()) == 11) )
167  H2.Fill( p->momentum().perp());
168  evt.clear();
169  }
170  inputA.close();
171 //Comparison
172  int diff=0;
173  for (int i=0; i<H1.GetNbinsX(); i++)
174  {
175  double eps=std::abs(H1.GetBinContent(i)-H2.GetBinContent(i));
176  if (eps<1e-5) continue;
177  std::cout<<"Bin: "<<i<<" "<<H1.GetBinContent(i)<<" "<<H2.GetBinContent(i)<<std::endl;
178  diff++;
179  }
180  H1.Print("All");
181  H2.Print("All");
182  return diff;
183 }
Definition of class GenParticle.
Definition of class WriterRootTree.
Definition of class ReaderRootTree.
Definition of class ReaderAsciiHepMC2.
Stores event-related information.
Definition: GenEvent.h:41
Definition of class WriterAsciiHepMC2.
GenEvent I/O parsing and serialization for root files based on root TTree.
int main(int argc, char **argv)
Definition of class GenEvent.
Feature< Feature_type > abs(const Feature< Feature_type > &input)
Obtain the absolute value of a Feature. This works as you&#39;d expect. If foo is a valid Feature...
Definition: Feature.h:316