ShowerInventorAnalyzer.cc

Go to the documentation of this file.

#include <sstream>

#include "ShowerInventorAnalyzer.h"

#include <utl/ErrorLogger.h>
#include <utl/TimeDistribution.h>
#include <utl/TimeDistributionAlgorithm.h>
#include <utl/Point.h>
#include <utl/Math.h>
#include <utl/Particle.h>

#include <evt/Event.h>

#include <sevt/SEvent.h>
#include <sevt/Station.h>
#include <sevt/PMT.h>
#include <sevt/PMTSimData.h>
#include <sevt/StationSimData.h>

#include <det/Detector.h>
#include <sdet/SDetector.h>
#include <sdet/Station.h>

#include <evt/ShowerSimData.h>

#include <fwk/CentralConfig.h>

using namespace std;
using namespace utl;
using namespace fwk;
using namespace evt;
using namespace sevt;

using namespace ShowerInventorAnalyzerNS;


VModule::ResultFlag
ShowerInventorAnalyzer::Init()
{
  INFO("ShowerInventorAnalyzer::Init()");

  Branch topB =
    CentralConfig::GetInstance()->GetTopBranch("ShowerInventorAnalyzer");

  topB.GetChild("OutputFilename").GetData(fOutFile);

  return eSuccess;
}


VModule::ResultFlag
ShowerInventorAnalyzer::Run(evt::Event& event)
{
  INFO("ShowerInventorAnalyzer::Run()");

  ShowerSimData& shower = event.GetSimShower();
  CoordinateSystemPtr sCs = shower.GetShowerCoordinateSystem();
  CoordinateSystemPtr lCs = shower.GetLocalCoordinateSystem();

  double zen = shower.GetZenith();
  cout << "----> in analyzer: zenith = " << zen/deg << endl;

  // histogram of total number of pe vs radius (in shower plane) for this zenith

  ostringstream peHistTitleSm;
  peHistTitleSm << "integrated pe number for zenith : " << zen/deg;
  string peHistTitle(peHistTitleSm.str());

  TH1F* peHist;
  if (fTankHistos.find(peHistTitle) == fTankHistos.end()) {
    peHist = new TH1F(peHistTitle.c_str(),
                      peHistTitle.c_str(),
                      501, 0*m, 5000*m);
  } else {
    peHist = fTankHistos[peHistTitle];
  }

  // histogram of number of muons vs radius (in shower plane) for this zenith

  ostringstream nMuHistTitleSr;
  nMuHistTitleSr << "muon number for zenith : " << zen/deg;
  string nMuHistTitle(nMuHistTitleSr.str());

  TH1F* nMuHist;
  if (fTankHistos.find(nMuHistTitle) == fTankHistos.end()) {
    nMuHist = new TH1F(nMuHistTitle.c_str(),
                       nMuHistTitle.c_str(),
                       501, 0*m, 5000*m);
  } else {
    nMuHist = fTankHistos[nMuHistTitle];
  }

  const SEvent& sevent = event.GetSEvent();
  for (SEvent::ConstStationIterator sIt = sevent.StationsBegin();
       sIt != sevent.StationsEnd(); ++sIt) {

    if (sIt->HasSimData()) {

      // Check if any PMT has sim data rather than checking to see
      // if there are particles in the station since it is more convenient
      // to not export particles into the offline file

      bool stationIsHit = false;
      const StationSimData& sim = sIt->GetSimData();
      for (Station::ConstPMTIterator pmtIt = sIt->PMTsBegin();
           pmtIt != sIt->PMTsEnd(); ++pmtIt) {
        if (pmtIt->HasSimData())
          stationIsHit = true;
      }

      if (!stationIsHit)
        continue;

      const unsigned int nMu = sim.GetNumberOfMuons();
      cout << "number of muons " << nMu << endl;

      double sumPe = 0;
      for (Station::ConstPMTIterator pmtIt = sIt->PMTsBegin();
           pmtIt != sIt->PMTsEnd(); ++pmtIt) {

        if (!pmtIt->HasSimData())
          continue;

        const PMTSimData& pmtSim = pmtIt->GetSimData();

        if (!pmtSim.HasPETimeDistribution())
          continue;

        const TimeDistributionI peDistro = pmtSim.GetPETimeDistribution();

        sumPe += TimeDistributionAlgorithm::Sum(peDistro, peDistro.GetStart(), peDistro.GetStop());

      } // loop on pmts

      // fill in histo with number pe's observed at distance from
      // shower axis (distance is in the shower frame).

      CoordinateSystemPtr lcs = shower.GetLocalCoordinateSystem();
      const sdet::Station& detStat = det::Detector::GetInstance().GetSDetector().GetStation(*sIt);

      CoordinateSystemPtr scs = shower.GetShowerCoordinateSystem();

      Point::Triple xyz = detStat.GetPosition().GetCoordinates(scs);
      double radius2 = Sqr(xyz.get<0>()) + Sqr(xyz.get<1>());

      peHist->Fill(sqrt(radius2), sumPe);
      nMuHist->Fill(sqrt(radius2), nMu);

    }

  }

  fTankHistos[peHistTitle] = peHist;
  fTankHistos[nMuHistTitle] = nMuHist;

  return eSuccess;
}


VModule::ResultFlag
ShowerInventorAnalyzer::Finish()
{
  INFO("ShowerInventorAnalyzer::Finish()");

  cout << "creating file : " << fOutFile << endl;

  TFile fOutputFile(fOutFile.c_str() ,"recreate");

  for (map<string, TH1F*>::iterator lIt = fTankHistos.begin();
       lIt != fTankHistos.end(); ++lIt) {
    lIt->second->Write();
  }

  fOutputFile.cd();
  fOutputFile.Close();

  return eSuccess;
}