SculptedAbsorberCompress.cc

Go to the documentation of this file.

#include <stdlib.h>
#include <string>
#include <fstream>
#include <sstream>
#include <vector>
#include <iostream>
//
// This is a translated and compiled version of SculptedAbosrber R script. It write a txt file on 
// /lbne/data, which has the Slice -2 and slice 30 Ntuple. Sile 30: there will two files: all the muon, and only 
// the muons that are not seen at slice 2 (created in the spoile, or downstream of it. 
//
// July 14, cross-checks with the new flag to signal the intercept with the slice in fron of HA.
// This is obsolete..  Since the G4 tracks appear in the ntuple out of order, too complicate and unnecessary to fix.. 
//
// However, our colleagues from Coloroda are trying to implement a similar algorythm, from their root ntuple. 
// So, let us fix this 
//
class SCHAmuon {
  public:
     int evt, tr, id, slice;
     double x, y, z, px, py, pz, ek, w, g2;
     inline void out(std::ofstream &fOut) const {
       fOut << " " << evt << " " << tr << " " << id << " ";
       fOut << " " << (x/10.) << " " << (y/10.) << " " << (z/10.);
       fOut << " " << (px) << " " << (py) << " " << (pz);
       fOut << " " << ek << " " << w << " " << g2 << std::endl;
     }
};

int main (int argc, char **argv) {

  const std::string dirOutName("/lbne/data/users/lebrun/SculptedAbsorber/");
//  const std::string dirPnfs("/pnfs/lbne/scratch/users/lebrun/LBNFSculptHA1f/");
//  const std::string tokenDir("LBNFSculptHA1g"); // 50 mRad, X offset -1.8 mRad 
  const std::string tokenDir("LBNFSculptHA2a"); // V2, nominal geometry 
//  const std::string tokenDir("LBNEAbsorbV1a"); // Straight, but default LBNE absorber.
  bool doLBNE =  (tokenDir.find("LBNE") != std::string::npos);
  std::string dirPnfs("/pnfs/lbne/scratch/users/lebrun/");
  dirPnfs += tokenDir + std::string("/");
  const std::string dirScratch("/scratch/lbne/lebrun/");
  const std::string jtitleHA("MuonFluxAtSculptedAbsorber_SculptedAbsorberProd1_");  
//  const std::string jtitleHA("MuonFluxAtSculptedAbsorber_SculptedAbsorberProd3_");  
//  const std::string jtitleHA("MuonFluxAtSculptedAbsorber_LBNEAbsorberProd1_");  
//  std::string header(" evt tr id slice x y z px py pz ek w "); 
  std::string header(" evt tr id x y z px py pz ek w g2 "); // skip the slice, file specific... 

//  const int jCluster = 2048727; // for LBNFSculptHA1f, Horizontal staggering of the cracks. 
//  const int jCluster = 2071904; // for LBNFSculptHA1f, Horizontal staggering of the cracks, + X angle of 0.25 mRad
//  const int jCluster = 2015073; // Cracks are aligned.. LBNFSculptHA1f
//  const int jCluster = 2109228; // Cracks are aligned.. LBNFSculptHA1f
//  const int jCluster = 2109228; // Cracks are not aligned.. LBNFSculptHA1g
  const int jCluster = 2623210; //LBNEAbsorbV1a 
  
  std::vector<SCHAmuon> datSl30Tmp;
  std::vector<SCHAmuon> datSlM2Tmp;
  std::vector<SCHAmuon> datSl30Cr;
  std::vector<SCHAmuon> datOneTr;
  
  std::ostringstream fOutSlM2StrStr; 
  fOutSlM2StrStr << dirOutName << "MuonFluxAtSculptedAbsorber_" << tokenDir << "_V2Slm2_" << jCluster << ".txt";
  std::string fOutSlM2Str(fOutSlM2StrStr.str());
  std::ofstream fOutSlM2(fOutSlM2Str.c_str());
  fOutSlM2 << header << std::endl;
  fOutSlM2.precision(6);
  
  
  std::ostringstream fOutSl30StrStr; 
  fOutSl30StrStr << dirOutName << "MuonFluxAtSculptedAbsorber_" << tokenDir << "_V2Sl30_" << jCluster << ".txt";
  std::string fOutSl30Str(fOutSl30StrStr.str());
  std::ofstream fOutSl30(fOutSl30Str.c_str());
  fOutSl30 << header << std::endl;
  fOutSl30.precision(6);
  
//  const size_t numFiles=50; // run e anf g 
  const size_t numFiles=200; //run a, with default absorber.. (and run g..) 
  int nBadFiles = 0;
//  const size_t numFiles=2; //run g
  size_t nSlM2=0; size_t nSl30=0;
//  for (size_t kf=0; kf != numFiles; kf++) {
  for (size_t kf=0; kf != 3; kf++) {
//    std::ostringstream cmdStrStr; 
//    cmdStrStr << "ifdh cp -D "<< dirPnfs << jtitleHA << jCluster << "_" << kf << ".txt " <<  dirScratch;
//    std::string cmdStr(cmdStrStr.str());
//    system(cmdStr.c_str());
    std::ostringstream aFNameStrStr;
    aFNameStrStr << dirScratch << jtitleHA << jCluster << "_" << kf << ".txt";
    std::string aFNameStr(aFNameStrStr.str());
    std::ifstream fIn(aFNameStr.c_str());
    if (!fIn.is_open()) {
      std::cerr << " Can't open file " << aFNameStr << std::endl;
//      return EXIT_FAILURE;
      nBadFiles++;
      continue; 
    } else {
      std::cerr << " Got file " << aFNameStr << std::endl;
    }
    size_t k=0;
    int prevEvt = -1;
    int prevTr = -1;
    while (fIn.good() &&(!fIn.eof())) {
      k++;
      char line[1024];
      fIn.getline(line, 1023);
      if (k == 1) continue;
      std::string linec(line);
      if (fIn.eof()) break;
      if (linec.length() < 5) continue;
      SCHAmuon a;
      std::istringstream sIn(linec);
      sIn >> a.evt;  sIn >> a.tr; sIn >> a.id; sIn >> a.slice;
      if ((a.slice != -2) && (a.slice != 30)) continue;
      sIn >> a.x; sIn >> a.y; sIn >> a.z;
      sIn >> a.px; sIn >> a.py; sIn >> a.pz;
      sIn >> a.ek; sIn >> a.w ; sIn >> a.g2;
      
      if (a.slice == -2) {
         a.out(fOutSlM2);
         datSlM2Tmp.push_back(a);
         nSlM2++;
      } 
      if (a.slice == 30) {
         if ((!doLBNE) || (doLBNE && (a.z < 236000.))) { // Take only the first detector plane, for sake of consistency with the 
                                                         // sculpted absorber...
           a.out(fOutSl30); 
           datSl30Tmp.push_back(a);
           nSl30++;
         }
      }
      if ((prevEvt !=  a.evt) || (prevTr != a.tr)) { // New track. Analyze the previous track. 
        size_t nSlPrev = datOneTr.size();
        if (nSlPrev != 0) {       
          if ((datOneTr[nSlPrev-1].slice == 30) && (datOneTr[0].slice != -2))  { 
             datSl30Cr.push_back(datOneTr[nSlPrev-1]);
          }
        }
        datOneTr.clear();
        datOneTr.push_back(a);
        prevEvt = a.evt;
        prevTr = a.tr;  
      } else {
        datOneTr.push_back(a);
      }
    }
    if (k < 2) {
      std::cerr << " One empty file " << aFNameStr << std::endl;
      nBadFiles++;
      continue; 
    }
//    std::ostringstream cmdDelStrStr; 
//    cmdDelStrStr << " rm -f  " << dirScratch << jtitleHA << jCluster << "_" << kf << ".txt";
//    std::string cmdDelStr(cmdDelStrStr.str());
//    system(cmdDelStr.c_str());
 //
  // Now pick the muons that are detected in slice 30, but not in slice 2. 
  // We have to do this on g4lbne run (or Dk2nu file) basis... 
  //
    std::vector<SCHAmuon>::const_iterator itm2 = datSlM2Tmp.begin();
    if (itm2 == datSlM2Tmp.end()) continue;
    std::cerr << " Number of muon in slice 30, all " << datSl30Tmp.size() << std::endl;
    size_t nSl30G2G = 0;
    size_t nSl30G2B = 0;
    for(std::vector<SCHAmuon>::const_iterator it30 = datSl30Tmp.begin(); it30 != datSl30Tmp.end(); it30++) {
      if (it30->g2 == 0) nSl30G2B++;
      else nSl30G2G++;
    }
    size_t nSl30CrG2G = 0;
    size_t nSl30CrG2B = 0;
    for(std::vector<SCHAmuon>::const_iterator it30 = datSl30Cr.begin(); it30 != datSl30Cr.end(); it30++) {
      if (it30->g2 == 0) nSl30CrG2B++;
      else nSl30CrG2G++;
    }
      
    std::cerr << " ... nSl30G2G " <<  nSl30G2G << " nSl30G2B " << nSl30G2B << std::endl;
    std::cerr << " ... nSl30CrG2G " <<  nSl30CrG2G << " nSl30CrG2B " << nSl30CrG2B << std::endl;
  } // On files.. 
  //
  fOutSlM2.close();
  fOutSl30.close();
  std::cerr << " Num Sl2 " << nSlM2 << " Slice 30 " << nSl30 << " Sl30 Cr " << datSl30Cr.size() << std::endl;

  std::ostringstream fOutSl30CrStrStr; 
  fOutSl30CrStrStr << dirOutName << "MuonFluxAtSculptedAbsorber_" << tokenDir << "_Sl30Cr_" << jCluster << ".txt";
  std::string fOutSl30CrStr(fOutSl30CrStrStr.str());
  std::ofstream fOutSl30Cr(fOutSl30CrStr.c_str());
  fOutSl30Cr.precision(6);
  fOutSl30Cr << header << std::endl;
  
  for (std::vector<SCHAmuon>::const_iterator it30Cr = datSl30Cr.begin(); it30Cr != datSl30Cr.end(); it30Cr++) {
    it30Cr->out(fOutSl30Cr);
  }
  fOutSl30Cr.close();
    std::cerr << " Number of empty or missing files: " << nBadFiles << std::endl;
}