convert.C
Go to the documentation of this file.
1 #include "TFile.h"
2 #include "TTree.h"
3 
4 #include <iostream>
5 
6 #include "sys/stat.h"
7 #include "wordexp.h"
8 
9 std::vector<std::string> wildcard(const std::string& wildcardString)
10 {
11  // Expand environment variables and wildcards like the shell would
12  wordexp_t p;
13  const int status = wordexp(wildcardString.c_str(), &p, WRDE_SHOWERR);
14 
15  if(status != 0){
16  std::cerr << "Wildcard string '" << wildcardString
17  << "' returned error " << status << " from wordexp()."
18  << std::endl;
19  return {};
20  }
21 
22  std::vector<std::string> fileList;
23 
24  for(unsigned int i = 0; i < p.we_wordc; ++i){
25  // Check the file exists before adding it
26  struct stat sb;
27  if(stat(p.we_wordv[i], &sb) == 0)
28  fileList.push_back(p.we_wordv[i]);
29  }
30 
31  wordfree(&p);
32 
33  return fileList;
34 }
35 
36 int gibuu_to_pdg(int id, int charge)
37 {
38  // See https://gibuu.hepforge.org/trac/wiki/ParticleIDs
39 
40  if(id == 1 && charge == -1) return -2212; // pbar
41  if(id == 1 && charge == 0) return 2112; // n
42  if(id == 1 && charge == +1) return +2212; // p
43 
44  if(id == 32 && charge == 0) return 3122; // Lambda0
45 
46  if(id == 33 && charge == -1) return 3112; // Sigma-
47  if(id == 33 && charge == 0) return 3212; // Sigma0
48  if(id == 33 && charge == +1) return 3222; // Sigma+
49 
50  if(id == 53 && charge == -1) return +3312; // Xi-
51  if(id == 53 && charge == 0) return 3322; // Xi0
52  if(id == 53 && charge == +1) return -3312; // Xi+
53 
54  if(id == 56 && charge == -1) return -4122; // Lambda_c-
55  if(id == 56 && charge == +1) return +4122; // Lambda_c+
56 
57  if(id == 57 && charge == -2) return -4222; // Sigma_c--
58  if(id == 57 && charge == -1) return -4212; // Sigma_c-
59  if(id == 57 && charge == 0) return 4112; // Sigma_c0
60  if(id == 57 && charge == +1) return 4212; // Sigma_c+
61  if(id == 57 && charge == +2) return 4222; // Sigma_c++
62 
63  if(id == 59 && charge == -1) return -4232;// Xi_c-
64  if(id == 59 && charge == 0) return 4132;// Xi_c0
65  if(id == 59 && charge == +1) return +4232;// Xi_c+
66 
67  if(id == 101 && charge == -1) return -211; // pi-
68  if(id == 101 && charge == 0) return 111; // pi0
69  if(id == 101 && charge == +1) return +211; // pi+
70 
71  if(id == 110 && charge == 0) return +311; // K0
72  if(id == 110 && charge == +1) return +321; // K+
73 
74  if(id == 111 && charge == -1) return -321; // K-
75  if(id == 111 && charge == 0) return -311; // K0bar
76 
77  if(id == 114 && charge == 0) return +421; // D0
78  if(id == 114 && charge == +1) return +411; // D+
79 
80  if(id == 115 && charge == -1) return -411; // D-
81  if(id == 115 && charge == 0) return -421; // D0bar
82 
83  if(id == 118 && charge == +1) return +431; // D_s+
84  if(id == 119 && charge == -1) return -431; // D_s-
85 
86  if(id == 901 && charge == -1) return +11; // e-
87  if(id == 901 && charge == +1) return -11; // e+
88  if(id == 902 && charge == -1) return +13; // mu-
89  if(id == 902 && charge == +1) return -13; // mu+
90  if(id == 903 && charge == -1) return +15; // tau-
91  if(id == 903 && charge == +1) return -15; // tau+
92 
93  if(id == 911 && charge == 0) return +12; // nu_e
94  if(id == 912 && charge == 0) return +14; // nu_mu
95  if(id == 913 && charge == 0) return +16; // nu_tau
96  if(id == -911 && charge == 0) return -12; // nu_ebar
97  if(id == -912 && charge == 0) return -14; // nu_mubar
98  if(id == -913 && charge == 0) return -16; // nu_taubar
99 
100 
101  std::cout << "Unknown particle: ID=" << id << ", charge=" << charge
102  << ". Please implement." << std::endl;
103  // abort();
104  return 0;
105 }
106 
107 void convert(std::string dir = "cc_numu/C12")
108 {
109  TFile* fOut = new TFile((dir+"/records.root").c_str(), "RECREATE");
110  TTree* trOut = new TTree("tr", "tr");
111  float Enu, weight;
112  trOut->Branch("Enu", &Enu);
113  trOut->Branch("weight", &weight);
114  int prod_id;
115  trOut->Branch("prod_id", &prod_id);
116  int nparts = 0;
117  int pdg[1000];
118  float E[1000];
119  float px[1000];
120  float py[1000];
121  float pz[1000];
122  trOut->Branch("nparts", &nparts);
123  trOut->Branch("pdg", &pdg, "pdg[nparts]/I");
124  trOut->Branch("E", &E, "E[nparts]/F");
125  trOut->Branch("px", &px, "px[nparts]/F");
126  trOut->Branch("py", &py, "py[nparts]/F");
127  trOut->Branch("pz", &pz, "pz[nparts]/F");
128 
129  int prevEvt = -1;
130 
131  const std::vector<std::string> fnames = wildcard(dir+"/*/FinalEvents.dat");
132 
133  for(unsigned int i = 0; i < fnames.size(); ++i){
134  const std::string fname = fnames[i];
135  std::cout << "Adding " << fname << std::endl;
136 
137  FILE* f = fopen(fname.c_str(), "r");
138  char header[4096];
139  fgets(header, 4096, f);
140 
141  int run, event, id, charge;
142  float perweight, pos1, pos2, pos3, p0, p1, p2, p3;
143  int history, prod_id_in;
144  float Enu_in;
145 
146  while(true){
147  int ret = fscanf(f, "%d %d %d %d %f %f %f %f %f %f %f %f %d %d %f",
148  &run, &event, &id, &charge,
149  &perweight, &pos1, &pos2, &pos3, &p0, &p1, &p2, &p3,
150  &history, &prod_id_in,
151  &Enu_in);
152 
153  if(ret != 15){
154  fclose(f);
155  if(nparts > 0) trOut->Fill(); // Catch the last event
156  nparts = 0; // Reset to the start of the particle list
157  break;
158  }
159 
160  if(prevEvt == -1) prevEvt = event;
161 
162  // If it's a new event, write out the old one first
163  if(event != prevEvt){
164  prevEvt = event;
165  trOut->Fill();
166  // Reset to the start of the particle list
167  nparts = 0;
168  }
169 
170  Enu = Enu_in;
171  prod_id = prod_id_in;
172 
173  // Skip particles with zero weight (TODO?)
174  if(perweight > 0){
175  pdg[nparts] = gibuu_to_pdg(id, charge);
176  E[nparts] = p0;
177  px[nparts] = p1;
178  py[nparts] = p2;
179  pz[nparts] = p3;
180  ++nparts;
181 
182  weight = perweight/fnames.size();
183  }
184  } // end while
185  }
186 
187  fOut->cd();
188  trOut->Write();
189 
190  std::cout << "Wrote combination of " << fnames.size()
191  << " files with " << trOut->GetEntries() << " events" << std::endl;
192 }
::xsd::cxx::tree::id< char, ncname > id
Definition: Database.h:165
def history()
Definition: g4zmq.py:84
int status
Definition: fabricate.py:1613
const Var weight
const char * p
Definition: xmltok.h:285
OStream cerr
Definition: OStream.cxx:7
int gibuu_to_pdg(int id, int charge)
Definition: convert.C:36
list fileList
Definition: cafExposure.py:30
void convert(std::string dir="cc_numu/C12")
Definition: convert.C:107
fclose(fg1)
Float_t E
Definition: plot.C:20
std::vector< std::string > wildcard(const std::string &wildcardString)
Definition: convert.C:9
OStream cout
Definition: OStream.cxx:6
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
TDirectory * dir
Definition: macro.C:5