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"
8 #include "HepMC3/ReaderRootTree.h"
9 #include "HepMC3/WriterRootTree.h"
10 #include "HepMC3/ReaderAsciiHepMC2.h"
11 #include "HepMC3/WriterAsciiHepMC2.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
Init(TChain * tree)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 }
SomeAnalysis(const std::string & file)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
main()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 }
184