ADST2CORSIKA.cc

Go to the documentation of this file.

#include <TROOT.h>
#include <TFile.h>
#include <TMath.h>
#include <TVector3.h>
#include <TSystem.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TProfile.h>

#include <iomanip>
#include <iostream>
#include <string>
#include <cmath>
#include <sstream>
#include <vector>
#include <algorithm>
#include <cmath>  
#include <cstdio> 
#include <iterator>
#include <fstream>
#include <algorithm> 
#include <map>

#include <RecEvent.h>
#include <RecEventFile.h>
#include <FileInfo.h>
#include <DetectorGeometry.h>
#include <EyeGeometry.h>
#include <FDEvent.h>
#include <EventInfo.h>
#include <Shower.h>
#include <GenShower.h>
#include <FdGenShower.h>
#include <TelescopeGeometry.h>
#include <FdTelescopeData.h>
#include <FdApertureLight.h>
#include <FdGeometry.h>
#include <FdApertureLight.h>

#define BOOST_SYSTEM_NO_DEPRECATED
#include <boost/algorithm/string.hpp> 
#include <boost/lexical_cast.hpp>
#include <boost/filesystem.hpp>

#include <chrono>
#include <thread>

using namespace std;


const bool uniqueFileNames = true; // discard path information


void 
writeSteeringCard(boost::filesystem::path steerdir, 
                  const int iRun, 
                  const double energyGeV, 
                  const double thetaDeg,
                  const double phiDeg,
                  const double arrangDeg,
                  const double obslevCM,
                  const double corexCM,
                  const double coreyCM)
{
  map<string,string> mapping;
  mapping["@RUNNR@"]  = boost::lexical_cast<std::string>(iRun);
  mapping["@ENERGY@"] = boost::lexical_cast<std::string>(energyGeV);
  mapping["@THETA@"]  = boost::lexical_cast<std::string>(thetaDeg);
  mapping["@PHI@"]    = boost::lexical_cast<std::string>(phiDeg);  
  mapping["@SEED1@"]  = boost::lexical_cast<std::string>(iRun*3+1);
  mapping["@SEED2@"]  = boost::lexical_cast<std::string>(iRun*3+2);
  mapping["@SEED3@"]  = boost::lexical_cast<std::string>(iRun*3+3);
  mapping["@ARRANG@"] = boost::lexical_cast<std::string>(arrangDeg);
  mapping["@OBSLEV@"] = boost::lexical_cast<std::string>(obslevCM);
  mapping["@COREX@"]  = boost::lexical_cast<std::string>(corexCM); 
  mapping["@COREY@"]  = boost::lexical_cast<std::string>(coreyCM);    
  ostringstream oname;
  oname << "DAT" << setw(6) << setfill('0') << iRun << ".inp";  
  if (boost::filesystem::exists(steerdir / oname.str())) {
    cout << "output steering card " << (steerdir / oname.str()).string() << " already exists!" << endl;
    exit(1);
  }
  ifstream input("DATxxxxxx_inp.template");
  ofstream output((steerdir / oname.str()).string().c_str());
  while (input.good() && output.good()) {
    string line;
    getline(input, line);
    for (const auto& replace : mapping) 
      boost::replace_all(line, replace.first, replace.second);
    output << line << endl;
  }
  input.close();
  output.close();
}


void 
ADST2CORSIKA(boost::filesystem::path steerdir, 
             const vector<string>& inputs) 
{
  const double degree = 180 / M_PI;
  
  map<string,map<int,int>> mapADST;
  map<int,pair<string,int>> mapCORS;
  {
    cout << "checking lock: " << (steerdir / "ShowerInfo2CORSIKA.txt-LOCK").string() << endl;
    bool isLocked = true;
    do {
      while (boost::filesystem::exists(steerdir / "ShowerInfo2CORSIKA.txt-LOCK")) { std::this_thread::sleep_for(std::chrono::milliseconds(1000)); /* wait */ } 
      std::ofstream lock ((steerdir / "ShowerInfo2CORSIKA.txt-LOCK").string());
      isLocked &= lock.good();
      lock << " ";
      isLocked &= lock.good();
      lock.close();
      isLocked &= lock.good();
    } while(!isLocked);
    cout << "locked: " << (steerdir / "ShowerInfo2CORSIKA.txt-LOCK").string() << endl;
    
    ifstream infile;
    infile.open ((steerdir / "ShowerInfo2CORSIKA.txt").string());
    if (infile.good()) {      
      string line;
      while(std::getline(infile, line)) {
        string adst = line.substr(line.find(" file=")+6, line.find(" ", line.find(" file=")+6)-line.find(" file=")-6);
        string evt_id = line.substr(line.find(" id=")+4, line.find(" ", line.find(" id=")+4)-line.find(" id=")-4);
        int evt = atoi((evt_id.substr(evt_id.find(":Use_")+5, evt_id.length()-evt_id.find(":Use_")-5)).c_str());
        int cors = atoi((line.substr(line.find(" corsikaID=")+11, line.find(" ", line.find(" corsikaID=")+11)-line.find(" corsikaID=")-11)).c_str());
        if (uniqueFileNames) 
          adst = adst.substr(adst.rfind('/')+1, adst.length()-adst.rfind('/')-1);
        mapADST[adst][evt] = cors;
        mapCORS[cors] = std::make_pair(adst, evt);
        // cout << line << " " << adst << " " << evt << " " << cors << endl;
      }
      infile.close();
    }
  }
  
  ofstream outfile;
  outfile.open ((steerdir / "ShowerInfo2CORSIKA.txt").string(), std::ofstream::app);
    
  RecEvent* theRecEvent = new RecEvent();
  DetectorGeometry* theGeometry = 0;
  theGeometry = new DetectorGeometry();
  FileInfo theInfo;

  unsigned int countFile = 0;
  unsigned int countEvent = 0;
  unsigned int countSelected = 0;
  unsigned int countSelectedNew = 0;
  unsigned int countSelectedAndReco = 0;
  unsigned int countHeatFOV = 0;
  unsigned int countMiss = 0; // count missing Cherenkov events
  
  for (const auto& input : inputs) {

    ++countFile;

    //if (countFile>100)
    //break;

    RecEventFile dataFile(input.c_str());
    dataFile.SetBuffers(&theRecEvent);
    dataFile.ReadDetectorGeometry(*theGeometry);
    dataFile.ReadFileInfo(theInfo);  
    
    //const double HeightOfRefCSOriginMeter = 1400; // from Offline xml
    double ObservationLevelMeter = 1300; // CORSIKA observation level [m]
    double diffangleDeg = -85.767; // CORSIKA Array rotation ARRANG [degree]
    
    unsigned int countInFile = 0;
    
    while(dataFile.ReadNextEvent() == RecEventFile::eSuccess){
      
      const string EventID = theRecEvent->GetEventId();         
      const int EventID_int = atoi((EventID.substr(EventID.find(":Use_")+5, EventID.length()-EventID.find(":Use_")-5)).c_str());
      
      ++countEvent;
      ++countInFile;
      
      TVector3 shower_axis = theRecEvent->GetGenShower().GetAxisSiteCS();
      
      // this return core at relative height in reference-CS (ePampaAmarilla)
      //#warning  there was a bug in corsika TELOUT. Check next time that this works:
      TVector3 shower_core = theRecEvent->GetGenShower().GetCoreAtAltitudeSiteCS(ObservationLevelMeter);//-HeightOfRefCSOriginMeter); 
      
      // find hottest Sim tel over all eyes
      
      int hottestEye = -1;
      int hottestTel = -1;
      double hottestTelPhotons = -1;
      for (RecEvent::EyeIterator eyeIt= theRecEvent->EyesBegin(); 
           eyeIt !=theRecEvent->EyesEnd(); 
           ++eyeIt){    
        const int eyeId = eyeIt->GetEyeId();
        //if (eyeId != 5) // only heat
        //continue;
        for (FDEvent::TelescopeDataConstIterator telIt = eyeIt->TelescopeDataBegin(); 
             telIt != eyeIt->TelescopeDataEnd(); ++telIt) {       
          const FdTelescopeData& telData = telIt->second;
          const FdApertureLight& theGenlight = telData.GetGenApertureLight();     
          const vector<double>& vecTotal = theGenlight.GetTotalLightAtAperture();
          double totalPhotonsTel = 0;
          for (const auto& itot : vecTotal) {totalPhotonsTel += itot;}
          if (totalPhotonsTel > hottestTelPhotons) {
            hottestEye = eyeId;
            hottestTel = telIt->first;
            hottestTelPhotons = totalPhotonsTel;
          }
        }
      }
      
      int selectedEye = 0;
      double genCherPhotons = 0;
      double genTotalPhotons = 0;
      double genMeanAngle = -1;
      
      if (hottestTelPhotons>0) { // no HEAT... ?
        
        const FDEvent& eye = theRecEvent->GetEye(hottestEye);
        const FdTelescopeData& telData = eye.GetTelescopeData(hottestTel);
        
        const int eyeId = eye.GetEyeId();
        selectedEye = eyeId;
      
        const FdApertureLight& theGenlight = telData.GetGenApertureLight();
        const vector<double>& vecCher = theGenlight.GetCherLightAtAperture();
        const vector<double>& vecTotal = theGenlight.GetTotalLightAtAperture();
        for (const auto& icher : vecCher) {genCherPhotons += icher;}
        for (const auto& itot : vecTotal) {genTotalPhotons += itot;}

        genMeanAngle = theGenlight.GetMeanAngle() * degree;

        /*
        if (genTotalPhotons>100 && genMeanAngle<40 && eyeId==5) {
          isHeatFOV = true;
        }
        */
      }
      
      //fill a text file with shower info to fed to CORSIKA steering cards. For showers with cherenkov fraction more than 50%
          
      int Prim = theRecEvent->GetGenShower().GetPrimary();
      double ZenithDeg = theRecEvent->GetGenShower().GetAxisSiteCS().Theta() * degree; //in degree
      double EnergyGeV = (theRecEvent->GetGenShower().GetEnergy())/pow(10,9);//crosika use energy in  Gev
      double Azimuth_correctedDeg =  theRecEvent->GetGenShower().GetAxisSiteCS().Phi() * degree; // in degree
      Azimuth_correctedDeg += diffangleDeg + 180;
      if (Azimuth_correctedDeg<0) Azimuth_correctedDeg += 360;
      if (Azimuth_correctedDeg>360) Azimuth_correctedDeg -= 360;
      double x_posCM =(shower_core.X()) * pow(10,2); //corsika use cm
      double y_posCM =(shower_core.Y()) * pow(10,2);
      //vdouble z_pos =(shower_core.Z()) * pow(10,2);     
      
      // DATyyyxxx
      // yyy: file ID
      // xxx: event ID
      int corsikaID = 1; //(countFile*1000 + countSelected%1000) % 1000000; // not above 999999
      // ostringstream << setw(3) << setfill('0') << countFile << setw(3) << setfill('0') << countSelected;
      
      ++countSelected;
      
      // input file can be "presel_adst.root" or just "adst.root"
      string input_test = input, input_test_alt;
      if (input.find("presel") == string::npos) {
        input_test_alt = input.substr(0, input.find("adst.root")) + "presel_adst.root";
      } else {
        input_test_alt = input.substr(0, input.rfind("presel_adst.root")) + "adst.root";
      }
      
      if (uniqueFileNames) {
        input_test = input_test.substr(input_test.rfind('/')+1, input_test.length()-input_test.rfind('/')-1);
        input_test_alt = input_test_alt.substr(input_test_alt.rfind('/')+1, input_test_alt.length()-input_test_alt.rfind('/')-1);
      } 
      
      if (mapADST.count(input_test)>0 && mapADST[input_test].count(EventID_int)>0) {
        
        cout << input << " evt=" << EventID_int << " already converted to CORSIKA steering card (" << mapADST[input_test][EventID_int] << ")" << endl;

      } else if (mapADST.count(input_test_alt)>0 && mapADST[input_test_alt].count(EventID_int)>0) {
        
        cout << input << " evt=" << EventID_int << " already converted to CORSIKA steering card (" << mapADST[input_test_alt][EventID_int] << ")" << endl;
        
      } else {
        
        while(mapCORS.count(corsikaID)>0) { 
          // cout << "found duplicate CORSIKA ID=" << corsikaID << " already. Avoid. Try next. " << endl;
          ++corsikaID; 
        }       
        
        ++countSelectedNew;
        
        if (true) {
          writeSteeringCard(steerdir,
                            corsikaID, 
                            EnergyGeV, 
                            ZenithDeg,
                            Azimuth_correctedDeg,
                            diffangleDeg,
                            ObservationLevelMeter*100,
                            x_posCM,
                            y_posCM);
          
          outfile << Prim << " " << EnergyGeV << " " << ZenithDeg <<" "<< Azimuth_correctedDeg << " "<< x_posCM << " " << y_posCM << " "<< " "<<ObservationLevelMeter*100. <<" "<< diffangleDeg << " file=" << input << " id="<<EventID << " eye=" << selectedEye << " tot=" << genTotalPhotons << " ch=" << genCherPhotons << " ViewAngle=" << genMeanAngle << " corsikaID=" << corsikaID << " EventID_int=" << EventID_int << endl;
        }
        
        mapADST[input][EventID_int] = corsikaID;
        mapCORS[corsikaID] = make_pair(input,EventID_int);      
      } // end new shower

    } // events loop
    
  } // file loop
  
  // ********************************************
  outfile.close();
  
  cout << "Finished: countFile=" << countFile << " countEvent=" << countEvent
       << " countSelected=" << countSelected << " countSelectedNew=" << countSelectedNew
       << " countSelectedAndReco=" << countSelectedAndReco
       << " countHeatFOV=" << countHeatFOV
       << " countMiss=" << countMiss << " (w/ chFrac>0.5)"
       << endl;

  boost::filesystem::remove( steerdir / "ShowerInfo2CORSIKA.txt-LOCK");
  cout << "unlocked: " << (steerdir / "ShowerInfo2CORSIKA.txt-LOCK").string() << endl;
}




int 
main(int argc, char** argv)
{
  if (argc<3) {
    cout << " specify steerdir AND input file" << endl;
    return 1;
  }
  const string steerdir(argv[1]);
  boost::filesystem::path dir(steerdir);
  if (!boost::filesystem::is_directory(dir)) {
    cout << " specify steerdir AND input file" << endl;
    cout << dir << " is not a directory" << endl;
    return 1;
  }
  
  vector<string> inputs;
  for (int i=1; i<argc-1; ++i) 
    inputs.push_back(string(argv[i+1]));  
  ADST2CORSIKA(steerdir, inputs);
  return 0;
}