CreateHybridLibrary_module.cc

Go to the documentation of this file.

// Chris Backhouse, UCL, Nov 2017

#include "larcore/Geometry/Geometry.h"
#include "larcorealg/Geometry/OpDetGeo.h"
#include "larsim/PhotonPropagation/PhotonVisibilityService.h"
#include "larsim/Simulation/PhotonVoxels.h"

#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"

#include "TBranch.h"
#include "TCanvas.h"
#include "TF1.h"
#include "TFile.h"
#include "TGraph.h"
#include "TH1.h"
#include "TStyle.h"
#include "TTree.h"
#include "TVectorD.h"

#include <iostream>

#define PI 3.14159265

namespace fhicl {
  class ParameterSet;
}

namespace phot {
  class CreateHybridLibrary : public art::EDAnalyzer {
  public:
    explicit CreateHybridLibrary(const fhicl::ParameterSet& p);

    // Plugins should not be copied or assigned.
    CreateHybridLibrary(const CreateHybridLibrary&) = delete;
    CreateHybridLibrary(CreateHybridLibrary&&) = delete;
    CreateHybridLibrary& operator=(const CreateHybridLibrary&) = delete;
    CreateHybridLibrary& operator=(CreateHybridLibrary&&) = delete;

    void analyze(const art::Event& e) override;
  };

  //--------------------------------------------------------------------
  CreateHybridLibrary::CreateHybridLibrary(const fhicl::ParameterSet& p) : EDAnalyzer(p)
  {

    art::ServiceHandle<geo::Geometry const> geom;

    art::ServiceHandle<phot::PhotonVisibilityService const> pvs;
    sim::PhotonVoxelDef voxdef = pvs->GetVoxelDef();

    TFile* fout_full = new TFile("full.root", "RECREATE");
    TFile* fout_fit = new TFile("fit.root", "RECREATE");

    std::cout << voxdef.GetNVoxels() << " voxels for each of " << geom->NOpDets() << " OpDets"
              << std::endl;
    std::cout << std::endl;

    //EP = Exception points: the parameterization is not a good description of the visibility and the value of the Photon Library is kept.
    long totExceptions = 0;
    long totPts = 0;

    for (unsigned int opdetIdx = 0; opdetIdx < geom->NOpDets(); ++opdetIdx) {
      std::cout << opdetIdx << " / " << geom->NOpDets() << std::endl;

      TDirectory* d_full = fout_full->mkdir(TString::Format("opdet_%d", opdetIdx).Data());
      TDirectory* d_fit = fout_fit->mkdir(TString::Format("opdet_%d", opdetIdx).Data());

      d_full->cd();
      TTree* tr_full = new TTree("tr", "tr");
      int vox, taken;
      float dist, vis, psi, theta, xpos;
      tr_full->Branch("vox", &vox);
      tr_full->Branch("dist", &dist);
      tr_full->Branch("vis", &vis);
      tr_full->Branch("taken", &taken);
      tr_full->Branch("psi",
                      &psi); //not needed to parameterize the visibilities, useful for tests in SP
      tr_full->Branch("theta",
                      &theta); //not needed to parameterize the visibilities, useful for tests in SP
      tr_full->Branch("xpos",
                      &xpos); //not needed to parameterize the visibilities, useful for tests in SP

      const geo::OpDetGeo& opdet = geom->OpDetGeoFromOpDet(opdetIdx);
      auto const opdetCenter = opdet.GetCenter();
      //std::cout << "OpDet position " << opdetIdx << ": " << opdetCenter << std::endl;

      struct Visibility {
        Visibility(int vx, int t, float d, float v, float p, float th, float xp)
          : vox(vx), taken(t), dist(d), vis(v), psi(p), theta(th), xpos(xp)
        {}
        int vox;
        int taken;
        float dist;
        float vis;
        float psi;
        float theta;
        float xpos;
      };
      TCanvas* c1 = new TCanvas("c1", "c1");
      //c1->SetCanvasSize(1500, 1500);
      c1->SetWindowSize(600, 600);
      //c1->Divide(1,2);
      TGraph g;
      TGraph g2;

      std::vector<Visibility> viss;
      viss.reserve(voxdef.GetNVoxels());

      for (unsigned int voxIdx = 0U; voxIdx < voxdef.GetNVoxels(); ++voxIdx) {

        const auto voxvec = voxdef.GetPhotonVoxel(voxIdx).GetCenter();
        const double xyzvox[] = {voxvec.X(), voxvec.Y(), voxvec.Z()};
        const double fc_y =
          600; //624cm is below the center of the first voxel outside the Field Cage
        const double fc_z = 1394;
        const double fc_x = 350;
        taken = 0;
        //DP does not need variable "taken" because all voxels are inside the Field Cage for the Photon Library created in LightSim.
        //DP taken = 1;
        dist = opdet.DistanceToPoint(voxvec);
        vis = pvs->GetVisibility(xyzvox, opdetIdx); // TODO use Point_t when available
        // all voxels outside the Field Cage would be assigned these values of psi and theta
        psi = 100;
        theta = 200;
        xpos = voxvec.X();

        if ((voxvec.X() - opdetCenter.X()) < 0) {
          psi = atan((voxvec.Y() - opdetCenter.Y()) / (-voxvec.X() + opdetCenter.X()));
        }

        if (voxvec.X() < fc_x && voxvec.X() > -fc_x && voxvec.Y() < fc_y && voxvec.Y() > -fc_y &&
            voxvec.Z() > -9 && voxvec.Z() < fc_z) {
          g.SetPoint(g.GetN(), dist, vis * dist * dist);
          taken = 1;
          psi = atan((voxvec.Y() - opdetCenter.Y()) / (voxvec.X() - opdetCenter.X())) * 180.0 /
                PI; // psi takes values within (-PI/2, PI/2)
          theta = acos((voxvec.Z() - opdetCenter.Z()) / dist) * 180.0 /
                  PI; // theta takes values within (0 (beam direction, z), PI (-beam direction, -z))

          if ((voxvec.X() - opdetCenter.X()) < 0) {
            psi = atan((voxvec.Y() - opdetCenter.Y()) / (-voxvec.X() + opdetCenter.X()));
          }
        }
        tr_full->Fill();
        viss.emplace_back(voxIdx, taken, dist, vis, psi, theta, xpos);
      } // end for voxIdx

      d_full->cd();
      tr_full->Write();
      delete tr_full;

      g.SetMarkerStyle(7);
      c1->cd();
      g.Draw("ap");
      g.GetXaxis()->SetTitle("Distance (cm)");
      g.GetYaxis()->SetTitle("Visibility #times r^{2}");
      g.Fit("expo");
      TF1* fit = g.GetFunction("expo");

      d_fit->cd();
      TVectorD fitres(2);
      fitres[0] = fit->GetParameter(0);
      fitres[1] = fit->GetParameter(1);
      fitres.Write("fit");

      gPad->SetLogy();

      TH1F h("", "", 200, 0, 20);

      d_fit->cd();
      TTree* tr_fit = new TTree("tr", "tr");
      tr_fit->Branch("vox", &vox);
      tr_fit->Branch("dist", &dist);
      tr_fit->Branch("vis", &vis);
      tr_fit->Branch("taken", &taken);
      tr_fit->Branch("psi", &psi);
      tr_fit->Branch("theta", &theta);
      tr_fit->Branch("xpos", &xpos);

      for (const Visibility& v : viss) {
        // 2e-5 is the magic scaling factor to get back to integer photon
        // counts. TODO this will differ for new libraries, should work out a
        // way to communicate it or derive it.
        const double obs = v.vis / 2e-5; //taken from the Photon Library
        const double pred =
          fit->Eval(v.dist) / (v.dist * v.dist) / 2e-5; //calculated with parameterization

        //DP const double obs = v.vis / 1e-7; //magic scaling factor for DP library created in LightSim
        //Minimal amount of detected photons is 50, bc of Landau dustribution
        //Those voxels with detected photons < 50 were set to 0
        //DP const double pred = fit->Eval(v.dist) / (v.dist*v.dist) / 1e-7; //calculated with parametrisation
        //std::cout << "observed = "<<obs<<" predicted (by parametrization) = "<<pred <<std::endl;

        // Log-likelihood ratio for poisson statistics
        double chisq = 2 * (pred - obs);
        if (obs) chisq += 2 * obs * log(obs / pred);

        vox = v.vox;
        dist = v.dist;
        vis = pred * 2e-5;
        psi = v.psi;
        theta = v.theta;
        xpos = v.xpos;
        //DP vis = pred *1e-7;

        if (v.taken == 1) { h.Fill(chisq); }

        if (
          chisq >
          9) { //equivalent to more than 9 chisquare = 3 sigma    //maybe play around with this cutoff
          g2.SetPoint(g2.GetN(), v.dist, v.vis * v.dist * v.dist);
          vis = obs * 2e-5;
          //DP vis = obs *1e-7;
          tr_fit->Fill();
        }
      }

      d_fit->cd();
      tr_fit->Write();
      delete tr_fit;
      g2.SetMarkerSize(1);
      g2.SetLineWidth(3);
      g2.SetMarkerStyle(7);
      g2.SetMarkerColor(kRed);
      gStyle->SetLabelSize(.04, "XY");
      gStyle->SetTitleSize(.04, "XY");
      c1->cd();
      g2.Draw("p same");
      c1->SaveAs(TString::Format("plots/Chris_vis_vs_dist_%d.png", opdetIdx).Data());

      h.GetXaxis()->SetTitle("#chi^{2}");
      h.GetYaxis()->SetTitle("Counts");
      gStyle->SetLabelSize(.04, "XY");
      gStyle->SetTitleSize(.04, "XY");
      h.Draw("hist");
      std::cout << "Integral(90,-1)/Integral(all) = " << h.Integral(90, -1) / h.Integral(0, -1)
                << std::endl;
      std::cout << "*****************************" << std::endl;
      totExceptions += h.Integral(90, -1);
      totPts += h.Integral(0, -1);
      gPad->Print(TString::Format("plots/chisq_opdet_%d.eps", opdetIdx).Data());

      delete c1;
    } // end for opdetIdx

    std::cout << totExceptions << " exceptions from " << totPts
              << " points = " << (100. * totExceptions) / totPts << "%" << std::endl;

    delete fout_full;
    delete fout_fit;
    exit(0); // We're done :)
  }

  //--------------------------------------------------------------------
  void CreateHybridLibrary::analyze(const art::Event&) {}

  DEFINE_ART_MODULE(CreateHybridLibrary)
} // namespace