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