convert.C

Go to the documentation of this file.

#include "TFile.h"
#include "TTree.h"

#include <iostream>

#include "sys/stat.h"
#include "wordexp.h"

std::vector<std::string> wildcard(const std::string& wildcardString)
{
  // Expand environment variables and wildcards like the shell would
  wordexp_t p;
  const int status = wordexp(wildcardString.c_str(), &p, WRDE_SHOWERR);

  if(status != 0){
    std::cerr << "Wildcard string '" << wildcardString
              << "' returned error " << status << " from wordexp()."
              << std::endl;
    return {};
  }

  std::vector<std::string> fileList;

  for(unsigned int i = 0; i < p.we_wordc; ++i){
    // Check the file exists before adding it
    struct stat sb;
    if(stat(p.we_wordv[i], &sb) == 0)
      fileList.push_back(p.we_wordv[i]);
  }

  wordfree(&p);

  return fileList;
}

int gibuu_to_pdg(int id, int charge)
{
  // See https://gibuu.hepforge.org/trac/wiki/ParticleIDs

  if(id ==   1 && charge == -1) return -2212; // pbar
  if(id ==   1 && charge ==  0) return  2112; // n
  if(id ==   1 && charge == +1) return +2212; // p

  if(id ==  32 && charge ==  0) return  3122; // Lambda0

  if(id ==  33 && charge == -1) return  3112; // Sigma-
  if(id ==  33 && charge ==  0) return  3212; // Sigma0
  if(id ==  33 && charge == +1) return  3222; // Sigma+

  if(id ==  53 && charge == -1) return +3312; // Xi-
  if(id ==  53 && charge ==  0) return  3322; // Xi0
  if(id ==  53 && charge == +1) return -3312; // Xi+

  if(id ==  56 && charge == -1) return -4122; // Lambda_c-
  if(id ==  56 && charge == +1) return +4122; // Lambda_c+

  if(id ==  57 && charge == -2) return -4222; // Sigma_c--
  if(id ==  57 && charge == -1) return -4212; // Sigma_c-
  if(id ==  57 && charge ==  0) return  4112; // Sigma_c0
  if(id ==  57 && charge == +1) return  4212; // Sigma_c+
  if(id ==  57 && charge == +2) return  4222; // Sigma_c++

  if(id ==  59 && charge == -1) return -4232;// Xi_c-
  if(id ==  59 && charge ==  0) return  4132;// Xi_c0
  if(id ==  59 && charge == +1) return +4232;// Xi_c+

  if(id == 101 && charge == -1) return  -211; // pi-
  if(id == 101 && charge ==  0) return   111; // pi0
  if(id == 101 && charge == +1) return  +211; // pi+

  if(id == 110 && charge ==  0) return  +311; // K0
  if(id == 110 && charge == +1) return  +321; // K+

  if(id == 111 && charge == -1) return  -321; // K-
  if(id == 111 && charge ==  0) return  -311; // K0bar

  if(id == 114 && charge ==  0) return  +421; // D0
  if(id == 114 && charge == +1) return  +411; // D+

  if(id == 115 && charge == -1) return  -411; // D-
  if(id == 115 && charge ==  0) return  -421; // D0bar

  if(id == 118 && charge == +1) return  +431; // D_s+
  if(id == 119 && charge == -1) return  -431; // D_s-

  if(id == 901 && charge == -1) return   +11; // e-
  if(id == 901 && charge == +1) return   -11; // e+
  if(id == 902 && charge == -1) return   +13; // mu-
  if(id == 902 && charge == +1) return   -13; // mu+
  if(id == 903 && charge == -1) return   +15; // tau-
  if(id == 903 && charge == +1) return   -15; // tau+

  if(id ==  911 && charge == 0) return   +12; // nu_e
  if(id ==  912 && charge == 0) return   +14; // nu_mu
  if(id ==  913 && charge == 0) return   +16; // nu_tau
  if(id == -911 && charge == 0) return   -12; // nu_ebar
  if(id == -912 && charge == 0) return   -14; // nu_mubar
  if(id == -913 && charge == 0) return   -16; // nu_taubar


  std::cout << "Unknown particle: ID=" << id << ", charge=" << charge
            << ". Please implement." << std::endl;
  //  abort();
  return 0;
}

void convert(std::string dir = "cc_numu/C12")
{
  TFile* fOut = new TFile((dir+"/records.root").c_str(), "RECREATE");
  TTree* trOut = new TTree("tr", "tr");
  float Enu, weight;
  trOut->Branch("Enu", &Enu);
  trOut->Branch("weight", &weight);
  int prod_id;
  trOut->Branch("prod_id", &prod_id);
  int nparts = 0;
  int pdg[1000];
  float E[1000];
  float px[1000];
  float py[1000];
  float pz[1000];
  trOut->Branch("nparts", &nparts);
  trOut->Branch("pdg", &pdg, "pdg[nparts]/I");
  trOut->Branch("E", &E, "E[nparts]/F");
  trOut->Branch("px", &px, "px[nparts]/F");
  trOut->Branch("py", &py, "py[nparts]/F");
  trOut->Branch("pz", &pz, "pz[nparts]/F");

  int prevEvt = -1;

  const std::vector<std::string> fnames = wildcard(dir+"/*/FinalEvents.dat");

  for(unsigned int i = 0; i < fnames.size(); ++i){
    const std::string fname = fnames[i];
    std::cout << "Adding " << fname << std::endl;

    FILE* f = fopen(fname.c_str(), "r");
    char header[4096];
    fgets(header, 4096, f);

    int run, event, id, charge;
    float perweight, pos1, pos2, pos3, p0, p1, p2, p3;
    int history, prod_id_in;
    float Enu_in;

    while(true){
      int ret = fscanf(f, "%d %d %d %d %f %f %f %f %f %f %f %f %d %d %f",
                       &run, &event, &id, &charge,
                       &perweight, &pos1, &pos2, &pos3, &p0, &p1, &p2, &p3,
                       &history, &prod_id_in,
                       &Enu_in);

      if(ret != 15){
        fclose(f);
        if(nparts > 0) trOut->Fill(); // Catch the last event
        nparts = 0; // Reset to the start of the particle list
        break;
      }

      if(prevEvt == -1) prevEvt = event;

      // If it's a new event, write out the old one first
      if(event != prevEvt){
        prevEvt = event;
        trOut->Fill();
        // Reset to the start of the particle list
        nparts = 0;
      }

      Enu = Enu_in;
      prod_id = prod_id_in;

      // Skip particles with zero weight (TODO?)
      if(perweight > 0){
        pdg[nparts] = gibuu_to_pdg(id, charge);
        E[nparts] = p0;
        px[nparts] = p1;
        py[nparts] = p2;
        pz[nparts] = p3;
        ++nparts;

        weight = perweight/fnames.size();
      }
    } // end while
  }

  fOut->cd();
  trOut->Write();

  std::cout << "Wrote combination of " << fnames.size()
            << " files with " << trOut->GetEntries() << " events" << std::endl;
}