MCCheater_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// \file    MCCheater.cxx
// \brief
// \version $Id: MCCheater.cxx,v 1.16 2012-09-20 22:01:30 greenc Exp $
// \author  $Author: greenc $
////////////////////////////////////////////////////////////////////////
// C++ includes
#include <functional>
#include <cmath>
#include <iostream>
#include <map>
#include <string>
#include <vector>

// ROOT includes
#include "TArray.h"
#include "TArrayD.h"
#include "TH1.h"
#include "TH2.h"
#include "TMatrixD.h"
#include "TObject.h"

// ART includes
#include "art/Framework/Principal/Event.h" 
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h" 
#include "art/Framework/Services/Optional/TFileDirectory.h" 
#include "art/Framework/Services/Optional/TFileService.h" 
#include "art/Framework/Services/Registry/ServiceHandle.h" 
#include "art/Framework/Principal/Handle.h" 
#include "canvas/Persistency/Common/Ptr.h"
#include "canvas/Persistency/Common/PtrVector.h"
#include "fhiclcpp/ParameterSet.h" 
#include "messagefacility/MessageLogger/MessageLogger.h" 

// NOvA includes
#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "MCCheater/BackTracker.h"
#include "MCCheater/SimHit.h" 
#include "MCCheater/SimTrack.h"
#include "Simulation/FLSHit.h"
#include "Simulation/FLSHitList.h"
#include "Simulation/Particle.h"
#include "Simulation/ParticleNavigator.h"

class TH1F;
class TH2F;

namespace cheat 
{

  class MCCheater : public art::EDProducer {

    public:

       explicit MCCheater(fhicl::ParameterSet const &pset);  // the explicit tag tells the compiler that GENIEGen is not an art::ParameterSet
       virtual ~MCCheater();

      // Optional, read/write access to event       
       void produce(art::Event& evt); 

      // Optional if you want to be able to configure from event display, for example
//      void reconfigure(std::istream &is,
//                     std::ostream &os,
//                     fhicl::ParameterSet pset);

      // Optional use if you have histograms, ntuples, etc you want around for every event
       void beginJob();
       
    private:
      
      //double         prod[6][6];
      //double         tran[7][7];
      int            ntest;

      std::string      fMCTruthListLabel;
      std::string      fFLSHitListLabel;
//      std::string      fSimHitListLabel;
//      std::string      fSimTrackListLabel;

      //  histograms
      TH1F*       fMuE;
      TH1F*       fElE;
      TH1F*       felL;
      TH1F*       fmuL;
      TH1F*       ftotalE;
      TH1F*       ftotEXV;
      TH1F*       ftotEYV;
      TH1F*       fdcosX;
      TH1F*       fdcosY;
      TH1F*       fdcosZ;
    //TH1F*       fdelE;
      TH1F*       fsimhitE;
    //TH1F*       fEpar;
      TH1F*       fFLSE;
      TH1F*       fnumsh;
      TH1F*       fNumFLS;
      TH1F*       ftlen;
    //TH1F*       fshlen;
      TH1F*       fthreshE;
      TH1F*       fnplanes;
    //TH2F*       fVertx;
    //TH2F*       fFidVert;
    //TH1F*       fzvtx;
      TH1F*       fparte;
      TH1F*       fnpart;
      TH1F*       fntrajpts;
    //TH1F*       fXRange;
    //TH1F*       fYRange;
    //TH1F*       fZRange;
    //TH1F*       fnshptr;
      TH1F*       fcorx;
      TH1F*       fcory;
    //TH2F*       fstsxy;
    //TH2F*       fstexy;
      TH1F*       fclhit;
      TH1F*       fzstart;

    };
////////////////////////////////////////////////////////////////////////
  MCCheater::MCCheater(fhicl::ParameterSet const& pset) :
    fMCTruthListLabel  (pset.get< std::string >("MCTruthListLabel" )),
    fFLSHitListLabel   (pset.get< std::string >("FLSHitListLabel"  ))
  {
    // tell module what its making
    produces< std::vector<cheat::SimHit> >();
    produces< std::vector<cheat::SimTrack>  >();
    //  produces< std::vector<cheat::SimVertex>  >();
    //  produces< std::vector<cheat::SimCluster>  >();

    // Be noisy to demonstrate what's happening
    mf::LogInfo("MCCheater") << " MCCheater::MCCheater()\n";

    // print this to keep track of version number
    mf::LogInfo("MCCheater") << "          ****************************************\n"
                             << "          *     MCCheater Version Number 4       *\n"
                             << "          ****************************************\n";

  }

  //..............................................................
  MCCheater::~MCCheater()
  {
    //======================================================================
    // Clean up any memory allocated by your module
    //======================================================================
  }

  //......................................................................
  void MCCheater::beginJob()
  {

    //
    // Book histograms, ntuples, initialize counts etc., etc., ...
    //
    art::ServiceHandle<art::TFileService> tfs;

    fMuE      = tfs->make<TH1F>("fMuE", ";True Muon Energy (GeV)", 50, 0., 10.);
    fElE      = tfs->make<TH1F>("fElE", ";True Electron Energy (GeV)", 50, 0., 10.);
    felL      = tfs->make<TH1F>("felL", ";length(cm)", 60, 0., 1200.);
    fmuL      = tfs->make<TH1F>("fmuL", ";length(cm)", 60, 0., 1200.);
    ftotalE   = tfs->make<TH1F>("ftotalE", ";E of SimHit (MeV)", 100, 0., 50.);
    ftotEXV   = tfs->make<TH1F>("ftotEXV", ";E of Xview SimHit (MeV)", 100, 0., 50.);
    ftotEYV   = tfs->make<TH1F>("ftotEYV", ";E of Yview SimHit (MeV)", 100, 0., 50.);
    fdcosX    = tfs->make<TH1F>("fdcosX", ";Dir cos X", 100, -1., 1.);
    fdcosY    = tfs->make<TH1F>("fdcosY", ";Dir cos Y", 100, -1., 1.);
    fdcosZ    = tfs->make<TH1F>("fdcosZ", ";Dir cos Z", 100, -1., 1.);
    fthreshE  = tfs->make<TH1F>("fthreshE", ";E of threshold SimHit (MeV)", 100, 0., 50.);
    fsimhitE  = tfs->make<TH1F>("fsimhitE", ";Total SimHit Energy per Event (GeV)", 200, 0., 20.);
    //      fshlen = tfs->make<TH1F>("fshlen", "Shower Length (cm)", 250, 0., 1000.);
    ftlen     = tfs->make<TH1F>("ftlen", ";Track Length (cm)", 250, 0., 500.);
    //      fEpar = tfs->make<TH1F>("fEpar", ";SimHit E per Event Tr 0 (GeV)", 250, 0., 5.);
    fFLSE     = tfs->make<TH1F>("fFLSE", ";Total FLSHit Energy per Event (GeV)", 200, 0., 20.);
    fnumsh    = tfs->make<TH1F>("fnumsh", ";Number of SimHits per Event", 50, -0.5, 999.5);
    fnplanes  = tfs->make<TH1F>("fnplanes", ";Number of Planes in Track", 100, 0., 200.);
    //      fVertx = tfs->make<TH2F>("fVertx", ";x (cm);y (cm)", 400, -800., 800., 400, -800., 800.);
    //      fstsxy = tfs->make<TH2F>("fstsxy", ";SimTrack start y vs x", 150, -150., 150., 240, -220., 220.);
    //      fstexy = tfs->make<TH2F>("fstexy", ";SimTrack end y vs x", 150, -150., 150., 240, -220., 220.);
    //      fFidVert = tfs->make<TH2F>("fFidfVert", ";x (cm);y (cm)", 400, -800., 800., 400, -800., 800.);
    //      fzvtx = tfs->make<TH1F>("fzvtx", ";Zvtx (cm)", 350, 0., 7000.);
    fNumFLS   = tfs->make<TH1F>("fNumFLS", ";Number of FLSHits in Event", 100, -0.5, 1999.5);
    fnpart    = tfs->make<TH1F>("fnpart", ";Number of Particles in Event", 50, -0.5, 49.5);
    fparte    = tfs->make<TH1F>("fparte", ";Total Energy of Particles in Event (GeV)", 200, 0., 20.);
    fntrajpts = tfs->make<TH1F>("fntrajpts", ";Number of Trajectory Points per Particle", 50, 0., 50.);
    //      fXRange = tfs->make<TH1F>("fXRange", ";Range of SimHit X Coords", 120, -210., 210.);
    //      fYRange = tfs->make<TH1F>("fYRange", ";Range of SimHit Y Coords", 120, -210., 210.);
    //      fZRange = tfs->make<TH1F>("fZRange", ";Range of SimHit Z Coords", 100, 0., 1000.);
    //      fnshptr = tfs->make<TH1F>("fnshptr", ";Number of SimHits per Track", 50, 0., 100.);
    fcorx     = tfs->make<TH1F>("fcorx", ";X Coord of xview plane", 120, -810., 810.);
    fcory     = tfs->make<TH1F>("fcory", ";Y Coord of yview plane", 120, -810., 810.);
    fclhit    = tfs->make<TH1F>("fclhit", ";Penetration Depth with SimHits", 125, 0., 250.);
    fzstart   = tfs->make<TH1F>("fzstart", ";Plane of particle start", 100, -1.0, 199.0);

  }

  void MCCheater::produce(art::Event& evt)
  {
    //======================================================================
    // Called for every event. "Reco" implies that you are adding
    // information to the event
    //======================================================================

    mf::LogInfo("MCCheater") << "MCCheater::produce() " << " New Reco Event\n";

    // increment counter for M matrix
    ntest++;

    /*
    //  get MCTruth info - should loop over this list since there may be more
    //  than 1 interaction, a neutrino event and a cosmic event for example
    //  other than that, don't need this list for the products of the cheater
    art::Handle< std::vector<simb::MCTruth> > mclist;
    evt.getByLabel(fMCTruthListLabel, mclist);
    if (mclist->empty())
    {
    mf::LogError("MCCheater") << "No MCTruth List - bummer\n";
    return;
    }

    // for now, only have 1 MCTruth per event, so only need mclist[0]
    //      art::PtrVector<simb::MCTruth> mctruth;  don't need this either?
    //      for (unsigned int i=0; i<mclist->size(); ++i)
    //      {
    art::Ptr<simb::MCTruth> mct(mclist, 0);
    std::cout << " Number of Particles from MCTruth = " << mct->NParticles() << "\n";
    //      }
      
    mf::LogInfo("MCCheater") << "Number of MCTruth lists = "<< mclist->size() << "\n";
      
      
    // fill a histogram with the size of each MCTruth 
    // for now, there is only 1 MCTruth list per MC event, so only need to use i=0
      
    for(unsigned int i=0; i < mclist->size(); ++i)
    {
    std::cout << " Number of particles from MCTruth = " << mctruth[i]->NParticles() << "\n";
    }
    */

    //  PARTICLE LIST NEEDED FOR CHEATER
    //  get Particle info - needed for Sim Objects
    art::ServiceHandle<cheat::BackTracker> bt;
    const sim::ParticleNavigator &pnav = bt->ParticleNavigator();

    //  First list needed for Cheater - the particle list
    //  using the particle list, make some plots, 
    //  make a map linking trackID with particles
    //  and a map linking trackID with its particleID - this is all that's needed for cheater
    std::map<int, const sim::Particle*> trpartmap;  // dont really need this map
    std::map<int, int> trpidmap;  // TrackID ParticleID map
    double PartE = 0.;
    unsigned int NPart = 0;
    int npar = 0;
    int pcount = 0;
    for (sim::ParticleNavigator::const_iterator j = pnav.begin(); j != pnav.end(); ++j){
      
      const sim::Particle p = *((*j).second);
      NPart += 1;
      npar += 1;

      int pid = p.PdgCode();  // what be this particle
      // energy of particle
      PartE += p.E();
      unsigned int ntrajpts = p.NumberTrajectoryPoints();
      fntrajpts->Fill(ntrajpts);
      //  for now, just check length of particle
      //  double plength = sqrt((enx-stx)*(enx-stx)+(eny-sty)*(eny-sty)+(enz-stz)*(enz-stz));
      if (abs(pid)==13) fMuE->Fill(p.E());  //  true muon energy for some reason?
      if (abs(pid)==11) fElE->Fill(p.E());  //  true electron energy
      int tid = p.TrackId();
      //  reduce number of particles with some obvious? cuts, no neutrinos, pi0s, photons, nuclei
      if (abs(pid)==12 || abs(pid)==14) continue;  // no neutrinos is only cut so far
      pcount++;
      trpartmap[tid] = &p;
      trpidmap[tid]=pid;
    }  // end loop over particles in this list
    mf::LogInfo("MCCheater") << "Final particle count after cuts in particle list = " << pcount << "\n";
    fparte->Fill(PartE);
    fnpart->Fill(NPart);
  
    //      mf::LogInfo("MCCheater") << "TrackID Particle Map size  " << trpartmap.size() << "\n";
    //      mf::LogInfo("MCCheater") << "TrackID ParticleID Map size " << trpidmap.size() << "\n";
    

    //  FLSHITS ALSO NEEDED FOR CHEATER
    // need FLSHits also - get FLSHit info
    art::Handle< std::vector<sim::FLSHitList>  > hitlist;
    evt.getByLabel(fFLSHitListLabel, hitlist);
    if (hitlist->empty()){
      mf::LogError("MCCheater") << "No FLSHit Lists - bummer\n";
      return;
    }
    //      mf::LogInfo("MCCheater") << "Number of FLSHit hitlists = " << hitlist->size() << "\n";
  
    //  Make maps if needed - these maps not necessary for the Cheater
    //  std::map<unsigned int, std::pair<int, int> > tidmgmmap;  // map linking track ID to its mother and grandmother
    std::map<unsigned int, std::vector<std::pair<unsigned int, unsigned int> > > tidpcpmap;
    //  std::map<std::pair<unsigned int, unsigned int>, std::vector<int> > pcptidmap;  // map linking pair plane cell to vector of TrackIDs
    //  std::map<std::pair<unsigned int, unsigned int>, std::vector<int> > pcpmotmap;  // map linking pair plane cell to vector of TrackID mothers
    //  std::map<std::pair<unsigned int, unsigned int>, std::vector<int> > pcppidmap;  // map linking pair plane cell to vector of particleIDs
    //  std::map<std::pair<unsigned int, unsigned int>, std::vector<double> > pcpendmap;  // map linking pair plane cell to vector of Edeps
    //  std::map<std::pair<unsigned int, unsigned int>, std::vector<double> > pcptimmap;  // map linking pair plane cell to vector of times
  
    // FLSHits are formed for each particle that deposits energy in a cell, so there
    // can be more than 1 FLSHit in a particular cell corresponding to each particle
    // that contributed to the energy.  And the cell IDs in FLSHits are not unique, when
    // combined with the PlaneId, then a unique channel ID for the whole detector can be determined.
    int nFLShits = 0;
    double flshitE = 0.;
    //      geo::Geometry&  geo = util::EDMUtils::GetGeometry(evt);
    art::ServiceHandle<geo::Geometry> fgeo;
    // Make a map linking FLSHits to TrackID
    std::map<int, std::vector<const sim::FLSHit*> > trhitmap;  // map linking TrackID to FLSHits
    // can define a map with a pair of keys and a value, so can make a map with plane and cell as key and
    // the value a vector of eg, FLSHits
    std::map<std::pair<unsigned int, unsigned int>, std::vector<const sim::FLSHit*> > pcpflsmap;  // map linking plane cell pair to FLShits
  
    for (unsigned int i=0; i<hitlist->size(); ++i)  // loop over FLSHit Lists - 1 per process interaction
      {
      
        art::Ptr<sim::FLSHitList> hlist(hitlist, i);
        //         mf::LogInfo("MCCheater") << "Num of FLSHits in list = " << hlist->fHits.size() << "\n";
      
        const std::vector<sim::FLSHit>& hits = hlist->fHits;
      
        nFLShits += hits.size();
        //  loop over FLSHits in this list
        for (unsigned int j=0; j<hits.size(); ++j)
          {
            int tid = hits[j].GetTrackID();
            if (tid<0) continue;  // just protection - shouldn't be here this is incoming particle (neutrino)
            // here, loop over particle list to get mother of this trackID
            int tmoth = -2;  // set mother to -2 for now
            int tgmoth = -2;  // set grandmother to -2 for now
            int tpid = 0;  
            // need to use only the particle id from particle list for this track - FLSHit pid can be different one
            for (sim::ParticleNavigator::const_iterator k = pnav.begin(); k != pnav.end(); ++k)
              {
                const sim::Particle part = *((*k).second);
                if (part.TrackId()==tid)
                  {
                    tpid = part.PdgCode();
                    tmoth = part.Mother();
                  }
                //  have to find grandma myself - 
                tgmoth = pnav.EveId(part.TrackId());
              }
            int smoth = tmoth;
            if (smoth > -2) smoth = 1;
            if (tmoth==-2) mf::LogInfo("MCCheater") << "                    THIS TrackID NOT IN PARTICLE LIST = " << tid << " " << hits[i].GetPDG() <<"\n";
            //  make a pair out of tmoth and tgmoth
            std::pair<int, int> mgmpair(tmoth, tgmoth);
            //              if (tpid>6000) continue;
            unsigned int planeid = hits[j].GetPlaneID();
            unsigned int cellid = hits[j].GetCellID();
            //  define the key as a plane cell pair here
            std::pair<unsigned int, unsigned int> pcpair(planeid, cellid);
            //  this is pair of trackID and particleID
            std::pair<unsigned int, unsigned int> tppair(tid, tpid);
            double xyzloc[3];
            double xyzwor[3];
            xyzloc[0] = 0.5*(hits[j].GetEntryX()+hits[j].GetExitX());
            xyzloc[1] = 0.5*(hits[j].GetEntryY()+hits[j].GetExitY());
            xyzloc[2] = 0.5*(hits[j].GetEntryZ()+hits[j].GetExitZ());
            fgeo->Plane(planeid)->Cell(cellid)->LocalToWorld(xyzloc,xyzwor);
            //              if (h.GetEdep()>0.001) trhitmap[tid].push_back(& h);
            // no selection - fill maps with everything
            trhitmap[tid].push_back(& hits[j]);
            pcpflsmap[pcpair].push_back(& hits[j]);
            double end = hits[j].GetEdep();
            flshitE += end;
            //  this is time, but don't use it yet
            //  double tim = hits[j].fT;
            //              geo::CellUniqueId cid = hits[j].fId;
            //  fill maps which will be used to make SimHits, so don't make any cuts yet
            //  if tpid is 22, a photon, don't use
          
            if (end>0.000000000001)
              {
                tidpcpmap[tid].push_back(pcpair);
                //                  tidmgmmap[tid].push_back(mgmpair);
                //                  tidmgmmap[tid]=mgmpair;
                //                  pcptidmap[pcpair].push_back(tid);
                //                  pcpmotmap[pcpair].push_back(tmoth);
                //                  pcppidmap[pcpair].push_back(tpid);
                //                  pcpendmap[pcpair].push_back(end);
                //                  pcptimmap[pcpair].push_back(tim);
              }
          
          }  // end of FLSHit hits loop
      }  // end of FLSHit lists loop
    fFLSE->Fill(flshitE);  // histogram of total FLSHit E
    fNumFLS->Fill(nFLShits);  // histogram of number of FLSHits
    /*
      std::cout << "   TrackID FLSHit Map size  " << trhitmap.size() << "\n";
      std::cout << "   TrackID PlaneCell pair Map size " << tidpcpmap.size() << "\n";
      std::cout << "   TrackID MotherGrandmother pair Map size = " << tidmgmmap.size() << "\n";
      std::cout << "   PlaneCell pair FLSHit Map size " << pcpflsmap.size() << "\n";
      std::cout << "   PlaneCell pair TrackID Map size " << pcptidmap.size() << "\n";
      std::cout << "   PlaneCell pair TID Mother Map size " << pcpmotmap.size() << "\n";
      std::cout << "   PlaneCell pair ParticleID Map size  " << pcppidmap.size() << "\n";
      std::cout << "   PlaneCell pair Edep Map size  " << pcpendmap.size() << "\n";
      std::cout << "   PlaneCell pair Tdep Map size  " << pcptimmap.size() << "\n";
    */ 
    // test some maps
  
    //      mf::LogInfo("MCCheater") << "Testing PlaneCell pair FLSHit map\n";
    /*
      for (std::map<std::pair<unsigned int, unsigned int>, std::vector<const sim::FLSHit*> >::iterator itm = pcpflsmap.begin(); itm != pcpflsmap.end(); ++itm)
      {
      std::pair<unsigned int, unsigned int> planecell = itm->first;
      std::vector<const sim::FLSHit*> harr = itm->second;
      std::cout <<"     Number of FLSHits = " << harr.size() << " For Plane and Cell = " << planecell.first <<  " " << planecell.second << "\n";
      }
    */
  
    // NEED THIS MAP FOR SIMTRACKS
    // find unique planecells for trackIDs, so can tell whether a track is longer than n cells
    // make a map linking each trackID with its unique planecells and the size of the value
    // planecell array is the number of unique planecells on this track
    std::map<unsigned int, std::vector<std::pair<unsigned int, unsigned int> > > tidupcpmap;
    for (std::map<unsigned int, std::vector<std::pair<unsigned int, unsigned int> > >::iterator itm = tidpcpmap.begin(); itm != tidpcpmap.end(); ++itm)
      {
        unsigned int trackid = itm->first;
        std::vector<std::pair<unsigned int, unsigned int> > pcarray = itm->second;
      
        //  eliminate duplicates, resize array with erase
        //  stl unique moves duplicates to the end of the array, doesn't eliminate
        //  erase deletes the items after the new location from unique
        //  this is part of that magic of C++ (R Rechenmacher)
        std::sort(pcarray.begin(), pcarray.end());
        std::vector<std::pair<unsigned int, unsigned int> >::iterator endloc;
        endloc = std::unique(pcarray.begin(), pcarray.end());
        pcarray.erase(endloc, pcarray.end());
        for (unsigned int i=0; i<pcarray.size(); ++i)
          {
            tidupcpmap[trackid].push_back(pcarray[i]);
          }
      }
  
  
    //  MAKE SIMHITS HERE
    //  make SimHits using the pc pair, FLSHit map
    double Etotev = 0.;
    unsigned int NumSH = 0;
    std::vector<cheat::SimHit> simhits;
    for (std::map<std::pair<unsigned int, unsigned int>, std::vector<const sim::FLSHit*> >::iterator itm = pcpflsmap.begin(); itm != pcpflsmap.end(); ++itm)
      {
        //          cheat::SimHit* hit = new cheat::SimHit();
        cheat::SimHit hit;
        std::pair<unsigned int, unsigned int> pcp = itm->first;  // plane-cell pair of this hit
        std::vector<const sim::FLSHit*> fh = itm->second;  // get the FLSHits in this pc pair
        unsigned int hsize = fh.size();  // number of FLSHits in this cell
      
        // can set plane cell and view for this SimHit now
        hit.SetPlane(pcp.first);                         // set plane of this hit
        unsigned int planen = pcp.first;                  // can use this later to set fiducial cut on fd
        hit.SetCell(pcp.second);                         // set cell of this hit
        unsigned int celln = pcp.second;                  // can use this later for fiducial cuts on fd
        // next 2 lines set view for this plane
        geo::View_t view = fgeo->Plane(pcp.first)->View();
        hit.SetView(view);                               // set view of this hit
        uint32_t ichan = 1;   //  have to find this somewhere?
        hit.SetChannel(ichan);
      
        //  make the arrays of E, T, PID, and TrID
        std::vector<double> earr;
        std::vector<double> carr;
        std::vector<int> tarr;
        std::vector<int> parr;
      
        std::vector<std::pair<int, int> > tparr;  // trackID, PID pair
        std::vector<std::pair<int, double> > tearr;
        std::vector<std::pair<int, double> > tcarr;
        double ene = 0.;
        double coorx = 0.;
        double coory = 0.;
        double coorz = 0.;
      
        double_t xxyyzz[3];
        fgeo->Plane(planen)->Cell(celln)->GetCenter(xxyyzz);
        coorz = xxyyzz[2];
        if (view==geo::kX)
          {
            coorx = xxyyzz[0];
            fcorx->Fill(coorx);
            //  need to get y from flshit
          } else if (view==geo::kY)
          {
            coory = xxyyzz[1];
            fcory->Fill(coory);
            //  need to get x from flshit
          }

        for (unsigned int i = 0; i<hsize; ++i)
          {
            earr.push_back(fh[i]->GetEdep());  //energy of this hit
            std::pair<int, double> tep(fh[i]->GetTrackID(), fh[i]->GetEdep());  // pair of trackID, hitE
            tearr.push_back(tep);
            carr.push_back(fh[i]->GetEntryT());  // time of this hit
            std::pair<int, double> ttp(fh[i]->GetTrackID(), fh[i]->GetEntryT());  // pair of trackID, time
            tcarr.push_back(ttp);

            tarr.push_back(fh[i]->GetTrackID());  // trackID of this hit
            //  get particle ID from TrackID ParticleID map, not from FLSHit since its different!
            for (std::map<int, int>::iterator itm = trpidmap.begin(); itm != trpidmap.end(); ++itm)
              {
                int trid = itm->first;
                if (trid == fh[i]->GetTrackID())
                  {
                    parr.push_back(itm->second);  // particle ID of this trackID
                    std::pair<int, int> tpp(fh[i]->GetTrackID(), itm->second);  // pair of trackID, particleID
                    tparr.push_back(tpp);
                  }
              }
            ene += fh[i]->GetEdep();  // energy sum of all contributions to this hit

            //  this uses first FLS contribution to get the non-view coordinate for a SimHit
            //  replace this with an array of the actual values for each contribution
            if (i==0)
              {
                double loc0[3];
                double wor0[3];
                double zexo = fh[0]->GetExitZ();
                double zeno = fh[0]->GetEntryZ();
                loc0[2] = 0.5*(zexo+zeno);
                double xexo = fh[0]->GetExitX();
                double xeno = fh[0]->GetEntryX();
                loc0[0] = 0.5*(xexo+xeno);
                double yexo = fh[0]->GetExitY();
                double yeno = fh[0]->GetEntryY();
                loc0[1] = 0.5*(yexo+yeno);
                //  change from local to global
                fgeo->Plane(planen)->Cell(celln)->LocalToWorld(loc0,wor0);
                if (view==geo::kX) coory = wor0[1];
                if (view==geo::kY) coorx = wor0[0];
              }
          }
        hit.SetTotE(ene);  // set E of this hit
        ftotalE->Fill(ene * 1000.);  // fill this hist also in units of MeV
        hit.SetPE(ene*5000.);  // took a guess to set the number of photoelectrons
        //         hit.SetPECorr(ene*5000.);  // took another guess to seet the number of corrected photoelectrons
        hit.SetCoorX(coorx);
        hit.SetCoorY(coory);
        hit.SetCoorZ(coorz);

        //  now need to finish definition based on number of FLSHits in this cell
        if (hsize == 1)  // there is only 1 hit in this cell, can make SimHit directly
          {
            hit.SetNP(1);  // just 1 here
            hit.SetEnD(earr); 
            hit.SetTim(carr);
            hit.SetTrID(tarr);
            hit.SetPID(parr);
          }
        else   // size greater than 1 combine if duplicates or mother-daughters
          {
            //  without sorting, can do with remove - then can adjust other arrays as well
            std::vector<int>::iterator startloc(tarr.begin());
            std::vector<int>::iterator endloc(tarr.end());
            int value(*startloc);

            while (++startloc < endloc)
              {
                endloc = std::remove(startloc, endloc, value);
                value = *startloc;
              }
            tarr.erase(endloc,tarr.end());          
            if (tarr.size()==1)
              {
                hit.SetNP(1);
                hit.SetTrID(tarr);
                std::vector<int> mparr (1,parr[0]);
                hit.SetPID(mparr);
                std::vector<double> mcarr (1,carr[0]);
                hit.SetTim(mcarr);
                double esum = 0.;
                for (unsigned int i=0; i<hsize; ++i)
                  {
                    esum += earr[i];
                  }
                std::vector<double> mearr (1,esum);
                hit.SetEnD(mearr);
              }
            else  // tarr size greater than 1
              {
                // if the size of the tarr array is the same as the size of tpp, for now, all hits are
                //  valid and unique and the arrays should be used to make a SimHit
                if (hsize==tarr.size())
                  {
                    hit.SetNP(hsize);
                    hit.SetPID(parr);
                    hit.SetEnD(earr);
                    hit.SetTrID(tarr);
                    hit.SetTim(carr);
                  }
                else
                  {
                    //  tarr size less than hsize so must analyze
                    //  can set the size to tarr size
                    hit.SetNP(tarr.size());
                    hit.SetTrID(tarr);

                    //  use pairs to re-generate PID, Energy, time arrays
                    std::vector<int> pvec(tarr.size());
                    std::vector<double> evec(tarr.size());
                    std::vector<double> tvec(tarr.size());
                    for (unsigned int i=0; i<tarr.size(); ++i)
                      {
                        for (unsigned int j=0; j<tparr.size(); ++j)
                          {
                            double lowt = 9999999999.;
                            if (tparr[j].first == tarr[i])
                              {
                                pvec[i] = tparr[j].second;
                                evec[i] += tearr[j].second;
                                if (tcarr[j].second > 0. && tcarr[j].second < lowt)
                                  {
                                    tvec[i] = tcarr[j].second;
                                    lowt = tcarr[j].second;
                                  }
                              }
                          }
                      }  // end of tarr loop
                    hit.SetPID(pvec);
                    hit.SetEnD(evec);
                    hit.SetTim(tvec);
                  }
              }
          }

        //  can put threshold and timing cuts here if wanted before making hit list
        //  0.0005 is half MeV threshold cut
        if( hit.fTotE > 0.00000005)  //  for now, just use a threshold cut
          {
            //  finally, put this hit into list
            simhits.push_back(hit);
            Etotev += hit.fTotE;
            NumSH++;
          }
      }

    fsimhitE->Fill(Etotev);  //fill histogram of total SimHit E in event
    fnumsh->Fill(NumSH);  // fill histogram of number of SimHits in Event
    std::cout << "                    Total E (GeV) in SimHits = " << Etotev << "\n";

    // check some SimHit stuff
    std::cout << "Size of SimHit List = " << simhits.size() << "\n";
    // loop over all SimHits to get Info
    for (unsigned int i=0; i<simhits.size(); ++i)
      {
        cheat::SimHit h=simhits[i];
        //        std::cout <<"  Num of particles for "<< i << " th hit = " << h.fNP << "\n"
        //                    <<"     Plane and Cell of this hit = " << h.fPlane << "  " << h.fCell << "\n"
        //                    <<"     View of this hit = " << h.fView << "\n"
        //                    <<"     Total E of this hit = "<< h.fTotE << "\n";
        // fill some histograms
        fthreshE->Fill(h.fTotE * 1000.);  // this is energy per hit in MeV
        if (h.View()==geo::kX) ftotEXV->Fill(h.fTotE * 1000.);  // this is energy per hit in MeV per X view
        if (h.View()==geo::kY) ftotEYV->Fill(h.fTotE * 1000.);  // this is energy per hti in MeV per Y view
        //          for (unsigned int j=0; j<h.fNP; ++j)
        //          {
        //            std::cout <<"          Particle Edep, ID, TrID, Time = "<< h.fEnD[j] << " " << h.fPID[j] << " " << h.fTrID[j] << " " << h.fTim[j] << "\n";
        //          }
      }
     
    //  MAKE SIMTRACKS HERE
    //  to make SimTracks, use TrackID unique PlaneCell pair map - if a TrackID has more than
    //  2 unique planecell pairs, both views represented, its a SimTrack
    //  will need a map linking trackID to SimHits for SimTracks
    //      std::map<unsigned int, art::Ptr<cheat::SimHit> > simtrhitmap;
    std::map<unsigned int, std::vector<cheat::SimHit*> > simtrhitmap;
    std::vector<cheat::SimTrack> simtracks;  // initialize SimTracks
    for (std::map<unsigned int, std::vector<std::pair<unsigned int, unsigned int> > >::iterator itm = tidupcpmap.begin(); itm != tidupcpmap.end(); ++itm)
      {
        unsigned int tid = itm->first;
        unsigned int upcsize = itm->second.size();
        //  if only 1 or 2 unique pc pairs, this is not reconstructable as a track
        //  also eventually probably want to add view requirement for both views
        if (upcsize>2)
          {
            //  this is a SimTrack
            cheat::SimTrack strack;
            //  first, get start and end coords from true particle length using trackID particle map
            int trpid = 0;
            // make some variables to hold the trajectory information
            TVector3 scoord; // Track start point
            bool foundStart = false;
            unsigned int startpoint = 0;
            TVector3 ecoord; // Track end point
            bool foundEnd = false;
            unsigned int endpoint = 0;
            std::vector<TLorentzVector> trpoints; // vector holding the particle trajectory points between the start and end pont.


            std::map<int, const sim::Particle*>::const_iterator it;
            for (it = trpartmap.begin(); it != trpartmap.end(); ++it)
              {
                unsigned int trid = it->first;
                const sim::Particle*  p=it->second;
                if (trid==tid)
                  {

                    trpid = p->PdgCode();

                    // Get the starting position of the track as the first trajectory point inside the detector volume
                    unsigned int trajPoint = 0;
                    while(!foundStart && trajPoint<p->NumberTrajectoryPoints()){
                      TLorentzVector start = p->Position(trajPoint);
                      if(TMath::Abs(start.X())<fgeo->DetHalfWidth() && TMath::Abs(start.Y())<fgeo->DetHalfHeight() && start.Z()>0.0 && start.Z()<fgeo->DetLength()){
                        scoord.SetXYZ(start.X(),start.Y(),start.Z());
                        foundStart = true;
                        startpoint = trajPoint;
                        //                        std::cout<<"Found the start point: "<<startpoint<<std::endl;
                      }
                      else{
                        ++trajPoint;
                      }
                    }
                    //initiate endpoint to same as startpoint
                    endpoint = startpoint;
                    trajPoint = p->NumberTrajectoryPoints()-1;
                    // Get the endpoint as the last trajectory point that is inside the detector volume
                    while(!foundEnd && trajPoint>=startpoint){
                      TLorentzVector end = p->Position(trajPoint);
                      if(TMath::Abs(end.X())<fgeo->DetHalfWidth() && TMath::Abs(end.Y())<fgeo->DetHalfHeight() && end.Z()>0.0 && end.Z()<fgeo->DetLength()){
                        ecoord.SetXYZ(end.X(),end.Y(),end.Z());
                        foundEnd = true;
                        endpoint = trajPoint;
                        //                        std::cout<<"Found the end point: "<<endpoint<<std::endl;
                      }
                      else{
                        --trajPoint;
                      }
                    }
                    // Save all the trajectory points after start point up to the end point.
                    if(foundStart && foundEnd && endpoint>startpoint){
                      for(unsigned int trpoint = startpoint+1; trpoint <= endpoint;++trpoint){
                        //                        std::cout<<"Start point: "<<startpoint<<" End point: "<<endpoint<<std::endl;
                        //                        std::cout<<"Putting trajectory point into track list: "<<trpoint<<std::endl;
                        trpoints.push_back(p->Position(trpoint));
                      }
                    }

                  }
              }
            // now only make a simtrack object if the track start and stop where found
            if(foundStart && foundEnd && endpoint>=startpoint){
              cheat::SimTrack strack;
              
              // calculate direction cosines from start and end of this MC track
              strack.SetPID(trpid);
              strack.SetTrackID(tid);
              strack.SetStart(scoord);
              for(unsigned int trpoint = 0; trpoint < trpoints.size();++trpoint){
                TVector3 pos(trpoints[trpoint].X(),trpoints[trpoint].Y(),trpoints[trpoint].Z());
                strack.AppendTrajectoryPoint(pos);
                // Get initial track direction
                if(trpoint == 0){
                  // calculate direction cosines from the first and second trajectory points of this MC track

                  double dist = sqrt((ecoord.X()-trpoints[trpoint].X())*(ecoord.X()-trpoints[trpoint].X())+
                                     (ecoord.Y()-trpoints[trpoint].Y())*(ecoord.Y()-trpoints[trpoint].Y())+
                                     (ecoord.Z()-trpoints[trpoint].Z())*(ecoord.Z()-trpoints[trpoint].Z()));
                  TVector3 dir((ecoord.X()-trpoints[trpoint].X())/dist,(ecoord.Y()-trpoints[trpoint].Y())/dist,(ecoord.Z()-trpoints[trpoint].Z())/dist);
                  strack.SetDir(dir);

                }
              }
              // add to SimTracks
              simtracks.push_back(strack);
            }

            //  loop over SimHits, adding to the TrackID, SimHit map
            //              simtrhitmap[tid].resize(0);  // Jon suggested this to fix map?
            for (unsigned int i=0; i<simhits.size(); ++i)
              {
                bool ismatch = false;
                //                  const cheat::SimHit h=simhits[i];
                cheat::SimHit h=simhits[i];
                for (unsigned int j=0; j<h.fNP; ++j)
                  {
                    if (h.fTrID[j]==(int) tid) ismatch = true;
                  }
                if (ismatch) simtrhitmap[tid].push_back(& h);
              }
          }  // end of upc greater than 2 loop
      }  // end of loop over trackID planecell pair map
    
    //  write SimTracks to event
    std::unique_ptr<std::vector<cheat::SimTrack> > simtrackcol(new std::vector<cheat::SimTrack>(simtracks) );
    evt.put(std::move(simtrackcol));

    /*
    //  **********************  Diagnostics for MCCheater  ***********************
    //  setup a test for single particles - if trajectory along z axis, test variable is true
    bool zdir = false;
    //  report some SimTrack stuff
    std::cout << "Size of SimTrack List = " << simtracks.size() << "\n";
    double tlensum = 0.;
    //  define min and max coords for NDOS
    double xmin = -130.1;
    double xmax = 130.1;
    double ymin = -195.1;
    double ymax = 195.1;
    double zmin = 0.;
    double zmax = 1200.;
    for (unsigned int i=0; i<simtracks.size(); ++i)
    {
    fzstart->Fill(simtracks[i].StartPos()[2]/6.0);
    double xlen = simtracks[i].EndPos()[0]-simtracks[i].StartPos()[0];
    double ylen = simtracks[i].EndPos()[1]-simtracks[i].StartPos()[1];
    double zlen = simtracks[i].EndPos()[2]-simtracks[i].StartPos()[2];
    fnplanes->Fill(zlen/6.);
    double tlength = sqrt(xlen*xlen+ylen*ylen+zlen*zlen);
    tlensum += tlength;
    fdcosX->Fill(simtracks[i].Dir()[0]);
    fdcosY->Fill(simtracks[i].Dir()[1]);
    fdcosZ->Fill(simtracks[i].Dir()[2]);
    if (simtracks[i].Dir()[2]>0.5) zdir = true;
    if (zdir && simtracks[i].TrackID()==0) ftlen->Fill(tlength);
    if (simtracks[i].Dir()[2]>0.5 && simtracks[i].PID()==13) fmuL->Fill(tlength);
    if (simtracks[i].Dir()[2]>0.5 && simtracks[i].PID()==11) felL->Fill(tlength);
    std::cout << "     Track Number " << i << "\n";
    std::cout << "     Track ID = " << simtracks[i].TrackID() << " Pdg ID = " << simtracks[i].PID() << "\n";
    std::cout << "     Start Coord = " << simtracks[i].StartPos()[0] << " " << simtracks[i].StartPos()[1] << " " << simtracks[i].StartPos()[2] << "\n";
    std::cout << "     End Coord = " << simtracks[i].EndPos()[0] << " " << simtracks[i].EndPos()[1] << " " << simtracks[i].EndPos()[2] << "\n";
    std::cout << "     Direction = " << simtracks[i].Dir()[0] << " " << simtracks[i].Dir()[1] << " " << simtracks[i].Dir()[2] << "\n";
    //  find point of intersection with plane of NDOS
    //  first, use top plane
    double ttop = (195.-simtracks[i].StartPos()[1])/ylen;
    double xint = simtracks[i].StartPos()[0]+ttop*(xlen);
    double yint = simtracks[i].StartPos()[1]+ttop*(ylen);
    double zint = simtracks[i].StartPos()[2]+ttop*(zlen);
    std::cout << "          TP ttrack = " << ttop << " Coords = " << xint << " " << yint << " " << zint << "\n";
    if (xint>xmin && xint<xmax && yint>ymin && yint<ymax && zint>zmin && zint<zmax)
    {
    //  loop over all simhits to find closest to this point
    double sdist = 999.;
    for (unsigned int j=0; j<simhits.size(); ++j)
    {
    cheat::SimHit h=simhits[j];
    //  does this hit contain contribution from this track?
    bool htid = false;
    for (unsigned int k=0; k<h.fNP; ++k)
    {
    if (h.fTrID[k]==simtracks[i].TrackID()) htid = true;
    }
    if (htid)
    {
    double hdist = sqrt((h.fCoorX-xint)*(h.fCoorX-xint)+(h.fCoorY-yint)*(h.fCoorY-yint)+(h.fCoorZ-zint)*(h.fCoorZ-zint));
    //                      std::cout << "    Hit distance = " << hdist << "\n";
    if (hdist<sdist) sdist = hdist;
    }  // satisfies trackID
    }  // end of SimHit loop
    std::cout << "                    Closest distance to edge = " << sdist << " for TrackID = " << simtracks[i].TrackID() << "\n";
    fclhit->Fill(sdist);
    }
    }
    */
      
    /*
      double E0SHTot = 0.;
      if (zdir)
      {
      for (unsigned int i=0; i<simhits.size(); ++i)
      {
      //              cheat::SimHit hit=simhits[i];
      E0SHTot += simhits[i].fTotE;
      }
      fEpar->Fill(E0SHTot); // histogram of forward-going energy sum
      }
    */

    //  check map of SimTrack TrackID SimHit
    //  this map doesn't work?
    /*
      for (std::map<unsigned int, std::vector<cheat::SimHit*> >::iterator itm = simtrhitmap.begin(); itm != simtrhitmap.end(); ++itm)
      {
      std::vector<cheat::SimHit*> sharr = itm->second;
      mf::LogInfo("MCCheater") <<"     Number of SimHits = " << sharr.size() << " For SimTrackID = " << itm->first << "\n";
      fnshptr->Fill(sharr.size());
      double stc[3];
      double enc[3];
      stc[1] = -210.;
      enc[1] = 210.;

      for (unsigned int i=0; i<sharr.size(); ++i)
      {
      cheat::SimHit* h=sharr[i];
      double xcor = h->fCoorX;
      double ycor = h->fCoorY;
      double zcor = h->fCoorZ;
      // histogram these coordinates
      fXRange->Fill(xcor);
      fYRange->Fill(ycor);
      fZRange->Fill(zcor);

      if (ycor > stc[1])
      {
      stc[1] = ycor;
      stc[2] = zcor;
      stc[0] = xcor;
      continue;
      } else if (ycor < stc[1])
      {
      if (ycor < enc[1])
      {
      enc[1] = ycor;
      enc[2] = zcor;
      enc[0] = xcor;
      continue;
      }
      }
      }
      fstsxy->Fill(stc[0],stc[1]);
      fstexy->Fill(enc[0],enc[1]);
      }
    */

    // write SimHits to event
    std::unique_ptr<std::vector<cheat::SimHit> > simhitcol(new std::vector<cheat::SimHit>(simhits) );
    std::cout << "          Wrote SimHits to event of size = " << simhits.size() << "\n";
    evt.put(std::move(simhitcol));

    return;
  }

   DEFINE_ART_MODULE(MCCheater)
}  // end namespace cheat
////////////////////////////////////////////////////////////////////////