ShowerMaker.cc

Go to the documentation of this file.

#include <fwk/CentralConfig.h>
#include <fwk/LocalCoordinateSystem.h>
#include <fwk/CoordinateSystemRegistry.h>
#include <utl/Point.h>
#include <evt/Header.h>
#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerSRecData.h>
#include <utl/TimeStamp.h>
#include <utl/AugerUnits.h>
#include <det/Detector.h>
#include <rdet/RDetector.h> 
#include <boost/tokenizer.hpp>

#include "ShowerMaker.h"


//Just for Ben 

#include <rdet/Station.h>
#include <TMath.h>
using namespace std;
using namespace evt;
using namespace utl;
using namespace fwk;
using namespace det;
typedef boost::tokenizer<boost::char_separator<char> > mytok;

namespace ShowerMaker { 

  ShowerMaker::ShowerMaker() : 
    fHeraldFilePath("/localdata/melissas/NewSimParams/shower_parameter_test_file.txt"),
    fHeraldStream(NULL),
    frand2(4173),
    // DEBUG
    foutfile(NULL),
    fH2_GroundMap(NULL),
    fDetGraph(NULL),
    fH2_Skymap(NULL),
    fH_energy(NULL),
    fH_theta(NULL),
    fH_phi(NULL)
    {
    }
  ShowerMaker::~ShowerMaker() { 
    if (fHeraldStream!=NULL) delete fHeraldStream;
  //for debug 
  /*
    if (foutfile!=NULL) delete foutfile;
    if (fH2_GroundMap!=NULL) delete fH2_GroundMap;
    if (fH2_Skymap!=NULL) delete fH2_Skymap;
    if (fH_energy!=NULL) delete fH_energy;*/
    //
    }

  fwk::VModule::ResultFlag ShowerMaker::Init() {
//Debug 
    foutfile= new TFile("out.root","RECREATE");
    fH2_GroundMap = new TH2F("GroundMap","GroundMap",700,-500,6500,350,-200,+3300);
    fH2_Skymap = new TH2F("skymap","sky",180,-2*TMath::Pi(),2*TMath::Pi(),180,-2*TMath::Pi(),2*TMath::Pi());
    fH_phi = new TH1F("phi","phi",180,0,360);
    fH_theta = new TH1F("theta","theta",45,0,90);
    fH_energy = new TH1F("nrj","nrj",30,17,20);
    fHeraldStream = new ifstream(fHeraldFilePath.c_str(), ios_base::in);
    string index("");
    boost::char_separator<char> sep(" ");
    getline(*fHeraldStream,index);
    mytok tok(index,sep);
    int ii=0;
    for(mytok::const_iterator titer=tok.begin();titer!=tok.end();++titer) { 
      if (*titer=="#") continue;
      fm_index[string(*titer)]=ii;
      ii++;
      }
    unsigned int seed=46150;
    frand2.SetSeed(seed);
    return VModule::eSuccess;
      }
    
    fwk::VModule::ResultFlag ShowerMaker::Run(evt::Event& theevent ) {
    if (fHeraldStream->eof()) {
      INFO("End Of File Reached"); 
      return eBreakLoop;
    }
    string zeline("");
    getline(*fHeraldStream,zeline);
//    event->MakeSEvent(); 
    theevent.MakeRecShower();
    ShowerRecData& Show= theevent.GetRecShower();
    Show.MakeSRecShower();
    ShowerSRecData& SdShow = Show.GetSRecShower();
    boost::char_separator<char> sep(" ");
    mytok tok(zeline,sep);
    vector<string> zetokens;
    for(mytok::const_iterator titer=tok.begin();titer!=tok.end();++titer) { 
      zetokens.push_back(*titer);
      }
//    for (int i=0;i<zetokens.size();++i) { 
  //    cout <<i<< "\t" << zetokens[i]<<endl;
//      }
    if (zetokens.size()<12)  { 
    
      for (unsigned int i=0;i<zetokens.size();++i) cout << " " << zetokens[i];
      cout <<endl;
      return eContinueLoop;}
    Detector& Det = Detector::GetInstance();
    Det.Update(utl::GetCurrentSystemTime());
    const rdet::RDetector& RadioDet = Det.GetRDetector();
//    Point refpos=RadioDet.GetStation(115).GetPosition(); // For AERA Ideal = AERA_13
      Point refpos=RadioDet.GetStation(113).GetPosition(); // For AERA Real = AERA_13
//    Point refpos=RadioDet.GetStation(10).GetPosition(); // For MAXIMA = M4
      
    // get reference coordinate system from detector (usually PampaAmarilla)      
    const CoordinateSystemPtr coordsys= Det.GetReferenceCoordinateSystem();
    float coreX=atof(zetokens[fm_index["fXCore"]].c_str())*m;
    float coreY=atof(zetokens[fm_index["fYCore"]].c_str())*m;
    float coreZ=refpos.GetZ(coordsys);
//    float coreZ=0;
    Point core(coreX,coreY,coreZ,coordsys);
    cout << "X : " << core.GetX(coordsys)/m << " Y: " << core.GetY(coordsys)/m << " Z : " << core.GetZ(coordsys)/m << endl;
    SdShow.SetCorePosition(core);
    //Get ID;
    //int eventid=0;
    //eventid=atoi(zetokens[fm_index["AugerId"]].c_str());
    
    //This piece of code is used to keep a monitoring plot of core positions
    //double x=0;
    //double y=0;
    Point refst153=RadioDet.GetStation(115).GetPosition();
    const CoordinateSystemPtr refCS= fwk::LocalCoordinateSystem::Create(refst153);
    if (fDetGraph==NULL) {
      const vector<int>& v_stid= RadioDet.GetFullStationList();
      fDetGraph = new TGraph(v_stid.size());
      int icount=0;
      for (vector<int>::const_iterator stiter=v_stid.begin();stiter!=v_stid.end();++stiter) { 
        Point pos=RadioDet.GetStation(*stiter).GetPosition();
        fDetGraph->SetPoint(icount,pos.GetX(refCS)*m,pos.GetY(refCS)*m);
        ++icount;
        }
      }
    /*************************************************************************************
     *    This bunch of code can generate a random core position in some position        *
     *                Please uncomment the case you want to generate                     *
     *************************************************************************************

      //Circle around AERA_11 
      Point statpos=RadioDet.GetStation(113).GetPosition();
      const CoordinateSystemPtr statCS = fwk::LocalCoordinateSystem::Create(statpos);
      double radius=550;
      this->Diskgen(frand2,x,y,radius);
      x*=m;
      y*=m;
      Point newcore(x,y,0,statCS);
      SdShow.SetCorePosition(newcore);
      fH2_GroundMap->Fill(newcore.GetX(refCS),newcore.GetY(refCS));
      //Circle around AERA_28 
      Point statpos=RadioDet.GetStation(128).GetPosition(); //AERA_28
      const CoordinateSystemPtr statCS = fwk::LocalCoordinateSystem::Create(statpos);
      double radius=1000;
      this->Diskgen(frand2,x,y,radius);
//      frand2.Circle(x,y,radius);
      x*=m;
      y*=m;
      Point newcore(x,y,0,statCS);
      SdShow.SetCorePosition(newcore);
      fH2_GroundMap->Fill(newcore.GetX(refCS),newcore.GetY(refCS));

      //Random in the AERA_ARRAY 
      x=0;y=0;
      while (1==1)   { 
        x=-200+frand2.Rndm()*6550;
        y=-250+frand2.Rndm()*3500;
        if (x<=3800 && y<0.51*x+1200) break; 
        if (x>3800 && x<=5200) break; 
        if (x>5200 && y<2000) break;

        x*=m;y*=m;
      Point newcore(x,y,0,refCS);
      fH2_GroundMap->Fill(x,y);
      SdShow.SetCorePosition(newcore);
      fH2_GroundMap->Fill(newcore.GetX(refCS),newcore.GetY(refCS));
      }
      ************************************************************************************
      *        End of the Random core setting                                            *
      ***********************************************************************************/

    //float energy=atof(zetokens[fm_index["Energy(CIC+FD_cal)"]].c_str())*EeV;
    float energy=atof(zetokens[fm_index["Energy(CICINFILL)"]].c_str())*EeV;
    fH_energy->Fill(TMath::Log10(energy/eV));
    SdShow.SetEnergy(energy,energy/4.);
    Show.SetEnergy(energy,energy/4.);
//    cout << " MM :: Energy ==" << zetokens[fm_index["Energy(CIC+FD_cal)"]] << "UTCTime " <<  zetokens[fm_index["UTCTime"]].c_str() <<endl ;
    cout << " MM :: Energy ==" << zetokens[fm_index["Energy(CICINFILL)"]] << "UTCTime " <<  zetokens[fm_index["UTCTime"]].c_str() <<endl ;
    cout << " MM :: " << fm_index["AugerId"] << " " << fm_index["NrStations"] << " " << fm_index["UTCTime"] << endl; 
//    Show.SetAzimuthAngle(atof(zetokens[11].c_str()));
//    Show.SetZenithAngle(atof(zetokens[10].c_str()));
//    CoordinateSystemPtr coresys=core.GetCoordinateSystem();
    CoordinateSystemPtr coresys=SdShow.GetCorePosition().GetCoordinateSystem();
//    utl::Vector axis(1/tan(atof(zetokens[11].c_str())*degree),1,cos(atof(zetokens[10].c_str())*degree),coresys);

    float theta=atof(zetokens[fm_index["Theta"]].c_str())*degree;
    float phi=atof(zetokens[fm_index["Phi"]].c_str())*degree;
//    float theta=atof(zetokens[fm_index["LocalTheta"]].c_str())*degree;
//    float phi=atof(zetokens[fm_index["LocalPhi"]].c_str())*degree;
    cout << " MM :: " << theta/degree << " ph "<< phi/degree <<endl;
    fH2_Skymap->Fill(theta,phi);
    fH_theta->Fill(theta/degree);
    fH_phi->Fill(phi/degree);
//    utl::Vector axis(cos(phi)*sin(theta),sin(phi)*sin(theta),cos(theta),coresys);
    utl::Vector axis(cos(phi)*sin(theta),sin(phi)*sin(theta),cos(theta),coordsys);
    SdShow.SetAxis(axis); 
    Header& head=theevent.GetHeader();
    head.SetId(zetokens[fm_index["AugerId"]]);
    TimeStamp ts(atoi(zetokens[fm_index["UTCTime"]].c_str()),atoi(zetokens[fm_index["fRec.fT0"]].c_str()));
    head.SetTime(ts);
    det::Detector::GetInstance().Update(head.GetTime());
    return VModule::eSuccess;
    }
  VModule::ResultFlag ShowerMaker::Finish() { 
    fHeraldStream->close();
    fH2_GroundMap->Write();
    fH2_Skymap->Write();
    fH_energy->Write();
    fDetGraph->Write();
//    fH_phi->Write();
    fH_theta->Write();
    foutfile->Close();
    return VModule::eSuccess;
    }
  void ShowerMaker::Diskgen(TRandom2& rnd, double& x, double& y, double r) { 
     double phi = rnd.Rndm()*TMath::TwoPi();
     double radius=rnd.Rndm()*r;
     x = radius*TMath::Cos(phi);
     y= radius*TMath::Sin(phi);
   } 
  } // namespace;