#include <TROOT.h>
#include <iostream>
#include <cometEventLoop.hxx>
#include <IG4Trajectory.hxx>
#include <IG4TrajectoryPoint.hxx>
#include <IG4PrimaryVertex.hxx>
#include <ICOMETEvent.hxx>
#include <IG4HitSegment.hxx>
#include <TH1.h>
#include <TH2.h>
#include <TH3.h>
#include <IHandle.hxx>
#include <TCanvas.h>
#include <vector>
#include <iomanip>
#include <TFile.h>
#include <typeinfo>


class IDumpEvent: public COMET::ICOMETEventLoopFunction {

public:

typedef std::vector<COMET::IG4TrajectoryPoint> Points; 

IDumpEvent() {

  fLSOption = "";
  fQuiet = false;
  fPrintList = false;
  showcounts = true;
  lsoption = false;
  trajcount = false;
  otherbyevent = false;

}

virtual ~IDumpEvent() {};

void Usage(void) {

        std::cout << "    -O quiet    Don't call TObject::ls()"
                  << std::endl;
        std::cout << "    -O list     List events in the file"
                  << std::endl;
        std::cout << "    -O <option> Set the option for TObject::ls()"
                  << std::endl;
}

virtual bool SetOption(std::string option,std::string value="") {

  if (value != "") return false;
  if (option == "quiet") fQuiet = true;
  else if (option == "list") {fPrintList = true;}	
  else fLSOption = option;
  return true;

}
void Initialize () {
  
  Histfile=new TFile("dump-event-ana.root","RECREATE");	

  h1a= new TH1F("Gamma Initial Momentum","Gamma Initial Momenta;Momentum (MeV/c);Counts",3000,0,0);
  h1b= new TH1F("Gamma Final Momentum","Gamma Final Momenta;Momentum (MeV/c);Counts",3000,0,0);
  H1a= new TH2F("Gamma Initial Position","Gamma Initial Position;x position;z position",3000,5000,-5000,3000,-5000,5000);
  H1b= new TH2F("Gamma Final Position","Gamma Final Position;x position;z position",3000,5000,-5000,3000,-5000,5000);

  h2a= new TH1F("Electron Initial Momentum","Electron Initial Momenta;Momentum (MeV/c);Counts",3000,0,0);
  h2b= new TH1F("Electron Final Momentum","Electron Final Momenta;Momentum (MeV/c);Counts",3000,0,0);
  H2a= new TH2F("Electron Initial Position","Electron Initial Position;x position;z position",3000,5000,-5000,3000,-5000,5000);
  H2b= new TH2F("Electron Final Position","Electron Final Position;x position;z position",3000,5000,-5000,3000,-5000,5000);

  h3a= new TH1F("Positron Initial Momentum","Positron Initial Momenta;Momentum (MeV/c);Counts",3000,0,0);
  h3b= new TH1F("Positron Final Momentum","Positron Final Momenta;Momentum (MeV/c);Counts",3000,0,0);
  H3a= new TH2F("Positron Initial Position","Positron Initial Position;x position;z position",3000,5000,-5000,3000,-5000,5000);
  H3b= new TH2F("Positron Final Position","Positron Final Position;x position;z position",3000,5000,-5000,3000,-5000,5000);

  h4a= new TH1F("Mouon Initial Momentum","Muon Initial Momenta;Momentum (MeV/c);Counts",3000,0,0);
  h4b= new TH1F("Muon Final Momentum","Muon Final  Momenta;Momentum (MeV/c);Counts",3000,0,0);
  H4a= new TH2F("Muon Initial Position","Muon Initial Position;x position;z position",3000,5000,-5000,3000,-5000,5000);
    H4aa= new TH2F("Muon 29.8 MeV/c Initial Position","Muon 29.8 MeV/c initial positions;;Counts",3000,-5000,5000,3000,-3000,8000);
  H4b= new TH2F("Muon Final Position","Muon Final Position;x position;z position",3000,-5000,5000,3000,-5000,5000);
    H4ba= new TH2F("Muon 29.8 MeV/c Final Position","Muon 29.8 MeV/c final positions;;Counts",3000,-5000,5000,3000,-3000,8000);

  h5a= new TH1F("Antimuon Initial Momentum","Antimuon Initial Momenta;Momentum (MeV/c);Counts",3000,0,0);
  h5b= new TH1F("Antimuon Final Momentum","Antimuon Final Momenta;Momentum (MeV/c);Counts",3000,0,0);
  H5a= new TH2F("Antimuon Initial Position","Antimuon Initial Position;x position;z position",3000,5000,-5000,3000,-5000,5000);
  H5b= new TH2F("Antimuon Final Position","Antimuon Final Position;x position;z position",3000,5000,-5000,3000,-5000,5000);

  h6a= new TH1F("Pion+ Initial Momentum","Pion+ Initial Momenta;Momentum (MeV/c);Counts",3000,0,0);
  h6b= new TH1F("Pion+ Final Momentum","Pion+ Final Momenta;Momentum (MeV/c);Counts",3000,0,0);
  H6a= new TH2F("Pion+ Initial Position","Pion+ Initial Position;x position;z position",3000,5000,-5000,3000,-5000,5000);
  H6b= new TH2F("Pion+ Final Position","Pion+ Final Position;x position;z position",3000,5000,-5000,3000,-5000,5000);

  h7a= new TH1F("Pion- Initial Momentum","Pion- Initial Momenta;Momentum (MeV/c);Counts",3000,0,0);
  h7b= new TH1F("Pion- Final Momentum","Pion- Final Momenta;Momentum (MeV/c);Counts",3000,0,0);
  H7a= new TH2F("Pion- Initial Position","Pion- Initial Position;x position;z position",3000,5000,-5000,3000,-5000,5000);
  H7b= new TH2F("Pion- Final Position","Pion- Final Position;x position;z position",3000,5000,-5000,3000,-5000,5000);

  h8a= new TH1F("Neutrino Inital Momentum","Neutrino Initial Momenta;Momentum (MeV/c);Counts",3000,0,0);
    h8ab= new TH1F("Electron Neutrino Initial Momentum","Electron neutrino initial momentum; Momentum (MeV/c);Counts",3000,0,350);
    h8ac= new TH1F("Electron Antineutrino Initial Momentum","Electron antineutrino initial momentum; Momentum (MeV/c);Counts",3000,0,350);
    h8ad= new TH1F("Muon Neutrino Initial Momentum","Muon neutrino initial momentum; Momentum (MeV/c);Counts",3000,0,350);
    h8ae= new TH1F("Muon Antineutrino Inital Momentum","Muon antineutrino initial momentum; Momentum (MeV/c);Counts",3000,0,350);
  h8b= new TH1F("Neutrino Final Momentum","Neutrino Final Momenta;Momentum (MeV/c);Counts",3000,0,0);
  H8a= new TH2F("Neutrino Inital Position","Neutrino Initial Position;x position;z position",3000,5000,-5000,3000,-5000,5000);
    H8aa= new TH2F("Neutrino 29.8 MeV/c Inital Position","Neutrino 29.8 MeV/c initial positions;;Counts",3000,-5000,5000,3000,-2000,7000);
  H8b= new TH2F("Neutrino Final Position","Neutrino Final Position;x position;z position",3000,5000,-5000,3000,-5000,5000);
    H8ba= new TH2F("Neutrino 29.8 MeV/c Final Position","Neutrino 29.8 MeV/c final positions;;Counts",3000,-5000,5000,3000,-5000,5000);

  h9a= new TH1F("Neutron Initial Momentum","Neutron Initial Momenta;Momentum (MeV/c);Counts",3000,0,0);
  h9b= new TH1F("Neutron Final Momentum","Neutron Final Momenta;Momentum (MeV/c);Counts",3000,0,0);
  H9a= new TH2F("Neutron Initial Position","Neutron Initial Position;x position;z position",3000,5000,-5000,3000,-5000,5000);
  H9b= new TH2F("Neutron Final Position","Neutron Final Position;x position;z position",3000,5000,-5000,3000,-5000,5000);

  h10a= new TH1F("Proton Initial Momentum","Proton Initial Momenta;Moemtnum (MeV/c);Counts",3000,0,0);
  h10b= new TH1F("Proton Final Momentum","Proton Final Momenta;Momentnum (MeV/c);Counts",3000,0,0);
  H10a= new TH2F("Proton Initial Position","Proton Initial Position;x position;z position",3000,5000,-5000,3000,-5000,5000);
  H10b= new TH2F("Proton Final Position","Proton Final Position;x position;z position",3000,5000,-5000,3000,-5000,5000);

  h11a= new TH1F("Other Initial Momentum","Other Initial Momenta;Momentum (MeV/c);Counts",3000,0,0);
  h11b= new TH1F("Other Final Momentum","Other Final Momenta;Momentum (MeV/c);Counts",3000,0,0);
  H11a= new TH2F("Other Initial Position","Other Initial Position;x position;z position",3000,5000,-5000,3000,-5000,5000);
  H11b= new TH2F("Other Final Position","Other Final Position;x position;z position",3000,5000,-5000,3000,-5000,5000);

  H12a= new TH3F("3D All Initial Position","Initial positions",100,-5000,5000,100,-5000,5000,100,2000,7000);
  H12b= new TH3F("3D All Final Position","Final positions",100,-5000,5000,100,-5000,5000,100,2000,7000);
  H12c= new TH3F("3D All Position","All positions",100,-5000,5000,100,-5000,5000,100,2000,7000);
  H12d= new TH3F("3D All Primary Vertices","All primary particle positions",100,-5000,5000,100,-5000,5000,100,5000,-5000);
  
  h13a= new TH1F("PP Gamma Initial Momentum","PP Gamma Initial Momenta;Momentum (MeV/c);Counts",3000,0,0);
  h13b= new TH1F("PP Electron Initial Momentum","PP Electron Initial Momenta;Momentum (MeV/c);Counts",3000,0,0);
  h13c= new TH1F("PP Positron Initial Momentum","PP Positron Initial Momenta;Momentum (MeV/c);Counts",3000,0,0);
  h13d= new TH1F("PP Muon Initial Momentum","PP Muon Initial Momenta;Momentum (MeV/c);Counts",3000,0,0);
  h13e= new TH1F("PP Antimuon Initial Momentum","PP Antimuon Initial Momenta;Momentum (MeV/c);Counts",3000,0,0);
  h13f= new TH1F("PP Pi+ Initial Momentum","PP Pi+ Initial Momenta;Momentum (MeV/c);Counts",3000,0,0);
  h13g= new TH1F("PP Pi- Initial Momentum","PP Pi- Initial Momenta;Momentum (MeV/c);Counts",3000,0,0);
  h13h= new TH1F("PP Neutrino Initial Momentum","PP Neutrino Initial Momenta;Momentum (MeV/c);Counts",3000,0,0);
  h13i= new TH1F("PP Neutron Initial Momentum","PP Neutron Initial Momenta;Momentum (MeV/c);Counts",3000,0,0);
  h13j= new TH1F("PP Proton Initial Momentum","PP Proton Initial Momenta;Momentum (MeV/c);Counts",3000,0,0);
  h13k= new TH1F("PP Other Initial Momentum","PP Other Initial Momenta;Momentum (MeV/c);Counts",3000,0,0);

  G1a= new TH3F("PP Gammma Initial Position","PP Gamma Initial Position;x position;z position",100,-5000,5000,100,-5000,5000,100,-5000,5000);
  G1b= new TH3F("PP Electron Initial Position","PP Electron Initial Position;x position;z position",100,5000,-5000,100,-5000,5000,100,-5000,5000);
  G1c= new TH3F("PP Postitron Initial Position","PP PositronInitial Position;x position;z position",100,5000,-5000,100,-5000,5000,100,-5000,5000);
  G1d= new TH3F("PP Muon Initial Position","PP Muon Initial Position;x position;z position",100,-5000,5000,100,-5000,5000,100,-5000,5000);
  G1e= new TH3F("PP Antimuon Initial Position","PP Antimuon Initial Position;x position;z position",100,5000,-5000,100,-5000,5000,100,-5000,5000);
  G1f= new TH3F("PP Pion+ Initial Position","PP Pion+ Initial Position;x position;z position",100,-5000,5000,100,-5000,5000,100,-5000,5000);
  G1g= new TH3F("PP Pion- Initial Position","PP Pion- Initial Position;x position;z position",100,-5000,5000,100,-5000,5000,100,-5000,5000);
  G1h= new TH3F("PP Neutrino Initial Position","PP Neutrino Initial Position;x position;z position",100,5000,-5000,100,-5000,5000,100,-5000,5000);
  G1i= new TH3F("PP Neutron Initial Position","PP Neutron Initial Position;x position;z position",100,5000,-5000,100,-5000,5000,100,-5000,5000);
  G1j= new TH3F("PP Proton Initial Position","PP Proton Initial Position;x position;z position",100,5000,-5000,100,-5000,5000,100,-5000,5000);
  G1k= new TH3F("PP Other Initial Position","PP Other Initial Position;x position;z position",100,5000,-5000,100,-5000,5000,100,-5000,5000);

  numall=0;
  numnull01=0;
  numnull02=0;
  numvertex=0;

  gammamin=3000.;
  electronmin=3000.;
  positronmin=3000.;
  muonmin=3000.;
  antimuonmin=3000.;
  pionpmin=3000.;
  pionnmin=3000.;
  neutrinomin=3000.;
  neutronmin=3000.;
  protonmin=3000.;
  othermin=3000.;

  ppgammamin=3000.;
  ppelectronmin=3000.;
  pppositronmin=3000.;
  ppmuonmin=3000.;
  ppantimuonmin=3000.;
  pppionpmin=3000.;
  pppionnmin=3000.;
  ppneutrinomin=3000.;
  ppneutronmin=3000.;
  ppprotonmin=3000.;
  ppothermin=3000.;

  gammamax=0.;
  electronmax=0.;
  positronmax=0.;
  muonmax=0.;
  antimuonmax=0.;
  pionpmax=0.;
  pionnmax=0.;
  neutrinomax=0.;
  neutronmax=0.;
  protonmax=0.;
  othermax=0.;

  ppgammamax=0.;
  ppelectronmax=0.;
  pppositronmax=0.;
  ppmuonmax=0.;
  ppantimuonmax=0.;
  pppionpmax=0.;
  pppionnmax=0.;
  ppneutrinomax=0.;
  ppneutronmax=0.;
  ppprotonmax=0.;
  ppothermax=0.;

}

virtual    bool operator () (COMET::ICOMETEvent& event) {

  if (lsoption) {
    if (COMET::ICOMETLog::GetLogLevel()>COMET::ICOMETLog::QuietLevel && !fQuiet) {
      event.ls(fLSOption.c_str());
      }
  }
  if (fPrintList) {
      COMETLog("Event " << event.GetContext());
        COMET::IHandle<COMET::IG4TrajectoryContainer> trajs=event.Get<COMET::IG4TrajectoryContainer>("truth/G4Trajectories");
    if (trajcount) {
      if (trajs) {
        std::cout<<"Number of trajectories "<<trajs->size()<<std::endl;
      }
      else {
        std::cout<<"Trajectories not added "<<std::endl;
        return true;
      }
    }
  if (!trajs) {numnull01++;}
  if (trajs) {
    for (COMET::IG4TrajectoryContainer::const_iterator k=trajs->begin();k!=trajs->end();k++) {
//k is a pair with first= track id, second=trajectory
      Int_t trackid=k->second.GetTrackId();
      COMET::IHandle<COMET::IG4Trajectory> traj=(*trajs).GetTrajectory(trackid);
      counts[(*traj).GetParticleName()]++;
      const Points& points=(*traj).GetTrajectoryPoints();
      const TVector3 initm=points.front().GetMomentum();
      const TVector3 finm=points.back().GetMomentum();
      const TLorentzVector initialm=(*traj).GetInitialMomentum();
      const TLorentzVector initp=(*traj).GetInitialPosition();
      const TLorentzVector finp=(*traj).GetFinalPosition();
      COMET::IHandle<COMET::IG4PrimaryVertexContainer> PVC=event.Get<COMET::IG4PrimaryVertexContainer>("truth/G4PrimVertex00");
      if (PVC) {
        for (COMET::IG4PrimaryVertexContainer::iterator it=PVC->begin();it!=PVC->end();it++) {
          try {
          COMET::TG4PrimaryParticleContainer& PPC=it->GetPrimaryParticles();
          if (&PPC==NULL) {throw 10;}
          if (&PPC!=NULL) {
        for (COMET::TG4PrimaryParticleContainer::iterator iter=PPC.begin();iter!=PPC.end();iter++) {
        Int_t ppPDG=iter->GetPDGCode();
        TLorentzVector ppinitp=it->GetPosition();
        TLorentzVector ppinitm4=iter->GetMomentum();
        TVector3 ppinitm=ppinitm4.Vect();
        if (ppPDG==22) {
          h13a->Fill(ppinitm.Mag());
          G1a->Fill(ppinitp[0],ppinitp[1],ppinitp[2]);
          if (ppinitm.Mag()>ppgammamax) {ppgammamax=ppinitm.Mag();}	
          if (ppinitm.Mag()<ppgammamin) {ppgammamin=ppinitm.Mag();}
          }
        else if (ppPDG==11) {
          h13b->Fill(ppinitm.Mag());
          G1b->Fill(ppinitp[0],ppinitp[1],ppinitp[2]);
          if (ppinitm.Mag()>ppelectronmax) {ppelectronmax=ppinitm.Mag();}
          if (initm.Mag()<ppelectronmin) {ppelectronmin=ppinitm.Mag();}
          }
        else if (ppPDG==-11) {
          h13c->Fill(ppinitm.Mag());
          G1c->Fill(ppinitp[0],ppinitp[1],ppinitp[2]);
          if (ppinitm.Mag()>pppositronmax) {pppositronmax=ppinitm.Mag();}
          if (ppinitm.Mag()<pppositronmin) {pppositronmin=ppinitm.Mag();}
          }
        else if (ppPDG==13) {
          h13d->Fill(ppinitm.Mag());
          G1d->Fill(ppinitp[0],ppinitp[1],ppinitp[2]);
          if (ppinitm.Mag()>ppmuonmax) {ppmuonmax=ppinitm.Mag();}
          if (ppinitm.Mag()<ppmuonmin) {ppmuonmin=ppinitm.Mag();}
          }
        else if (ppPDG==-13) {
          h13e->Fill(ppinitm.Mag());
          G1e->Fill(ppinitp[0],ppinitp[1],ppinitp[2]);
          if (ppinitm.Mag()>ppantimuonmax) {ppantimuonmax=ppinitm.Mag();}
          if (ppinitm.Mag()<ppantimuonmin) {ppantimuonmin=ppinitm.Mag();}
          }
        else if (ppPDG==211) {
          h13f->Fill(ppinitm.Mag());
          G1f->Fill(ppinitp[0],ppinitp[1],ppinitp[2]);
          if (ppinitm.Mag()>pppionpmax) {pppionpmax=ppinitm.Mag();}
          if (ppinitm.Mag()<pppionpmin) {pppionpmin=ppinitm.Mag();}
          }
        else if (ppPDG==-211) {
          h13g->Fill(ppinitm.Mag());
          G1g->Fill(ppinitp[0],ppinitp[1],ppinitp[2]);
          if (ppinitm.Mag()>pppionnmax) {pppionnmax=ppinitm.Mag();}
          if (ppinitm.Mag()<pppionnmin) {pppionnmin=ppinitm.Mag();}
          }
        else if (abs(ppPDG)==12||abs(ppPDG)==14||abs(ppPDG)==16) {
          h13h->Fill(ppinitm.Mag());
          G1h->Fill(ppinitp[0],ppinitp[1],ppinitp[2]);
          if (ppinitm.Mag()>ppneutrinomax) {ppneutrinomax=ppinitm.Mag();}
          if (ppinitm.Mag()<ppneutrinomin) {ppneutrinomin=ppinitm.Mag();}
          }
        else if (ppPDG==2112) {
          h13i->Fill(ppinitm.Mag());
          G1i->Fill(ppinitp[0],ppinitp[1],ppinitp[2]);
          if (ppinitm.Mag()>ppneutronmax) {ppneutronmax=ppinitm.Mag();}
          if (initm.Mag()<ppneutronmin) {ppneutronmin=ppinitm.Mag();}
          }
        else if (ppPDG==2212) {
          h13j->Fill(ppinitm.Mag());
          G1j->Fill(ppinitp[0],ppinitp[1],ppinitp[2]);
          if (ppinitm.Mag()>ppprotonmax) {ppprotonmax=ppinitm.Mag();}
          if (ppinitm.Mag()<ppprotonmin) {ppprotonmin=ppinitm.Mag();}
          }
        else {
          h13k->Fill(ppinitm.Mag());
          G1k->Fill(ppinitp[0],ppinitp[1],ppinitp[2]);
          if (ppinitm.Mag()>ppothermax) {ppothermax=ppinitm.Mag();}
          if (ppinitm.Mag()<ppothermin) {ppothermin=ppinitm.Mag();}
          }
        }
        }
        }
        catch (int x) {std::cout<<"Null PPC handle"<<std::endl;return false;}
      numvertex++;
        }
      }
    if (!PVC) {numnull02++;}
    H12a->Fill(initp[0],initp[1],initp[2]);
    H12b->Fill(finp[0],finp[1],finp[2]);
    for (Points::const_iterator n=points.begin();n!=points.end();n++) {
      const TLorentzVector positions=(*n).GetPosition();
      H12c->Fill(positions.X(),positions.Y(),positions.Z());
    };
    if (k->second.GetPDGEncoding()==22) {
      h1a->Fill(initm.Mag());
      h1b->Fill(finm.Mag());
      H1a->Fill(initp[0],initp[2]);
      H1b->Fill(finp[0],finp[2]);
      if (initm.Mag()>gammamax) {gammamax=initm.Mag();}	
      if (initm.Mag()<gammamin) {gammamin=initm.Mag();}
      }
    else if (k->second.GetPDGEncoding()==11) {
      h2a->Fill(initm.Mag());
      h2b->Fill(finm.Mag());
      H2a->Fill(initp[0],initp[2]);
      H2b->Fill(finp[0],finp[2]);
      if (initm.Mag()>electronmax) {electronmax=initm.Mag();}
      if (initm.Mag()<electronmin) {electronmin=initm.Mag();}
      }
    else if (k->second.GetPDGEncoding()==-11) {
      h3a->Fill(initm.Mag());
      h3b->Fill(finm.Mag());
      H3a->Fill(initp[0],initp[2]);
      H3b->Fill(finp[0],finp[2]);
      if (initm.Mag()>positronmax) {positronmax=initm.Mag();}
      if (initm.Mag()<positronmin) {positronmin=initm.Mag();}
      }
    else if (k->second.GetPDGEncoding()==13) {
      h4a->Fill(initm.Mag());
//Looking at initial and final position distributions of muons with p=29.79 MeV/c
      h4b->Fill(finm.Mag());
      if (initm.Mag()<29.8&&initm.Mag()>29.78) {
        H4aa->Fill(initp[0],initp[2]);
        H4ba->Fill(finp[0],finp[2]);
      }
      H4a->Fill(initp[0],initp[2]);
      H4b->Fill(finp[0],finp[2]);
      if (initm.Mag()>muonmax) {muonmax=initm.Mag();}
      if (initm.Mag()<muonmin) {muonmin=initm.Mag();}
      }
    else if (k->second.GetPDGEncoding()==-13) {
      h5a->Fill(initm.Mag());
      h5b->Fill(finm.Mag());
      H5a->Fill(initp[0],initp[2]);
      H5b->Fill(finp[0],finp[2]);
      if (initm.Mag()>antimuonmax) {antimuonmax=initm.Mag();}
      if (initm.Mag()<antimuonmin) {antimuonmin=initm.Mag();}
      }
    else if (k->second.GetPDGEncoding()==211) {
      h6a->Fill(initm.Mag());
      h6b->Fill(finm.Mag());
      H6a->Fill(initp[0],initp[2]);
      H6b->Fill(finp[0],finp[2]);
      if (initm.Mag()>pionpmax) {pionpmax=initm.Mag();}
      if (initm.Mag()<pionpmin) {pionpmin=initm.Mag();}
      }
    else if (k->second.GetPDGEncoding()==-211) {
      h7a->Fill(initm.Mag());
      h7b->Fill(finm.Mag());
      H7a->Fill(initp[0],initp[2]);
      H7b->Fill(finp[0],finp[2]);
      if (initm.Mag()>pionnmax) {pionnmax=initm.Mag();}
      if (initm.Mag()<pionnmin) {pionnmin=initm.Mag();}
      }
    else if (abs(k->second.GetPDGEncoding())==12||abs(k->second.GetPDGEncoding())==14||abs(k->second.GetPDGEncoding()==16)) {
      h8a->Fill(initm.Mag());
      h8b->Fill(finm.Mag());
      if (k->second.GetPDGEncoding()==12) {h8ab->Fill(initm.Mag());}
      if (k->second.GetPDGEncoding()==-12) {h8ac->Fill(initm.Mag());}
      if (k->second.GetPDGEncoding()==14) {h8ad->Fill(initm.Mag());}
      if (k->second.GetPDGEncoding()==-14) {h8ae->Fill(initm.Mag());}
      H8a->Fill(initp[0],initp[2]);
      H8b->Fill(finp[0],finp[2]);
      if (initm.Mag()<29.8&&initm.Mag()>29.78) {
        H8aa->Fill(initp[0],initp[2]);
        H8ba->Fill(finp[0],finp[2]);
      }
      if (initm.Mag()>neutrinomax) {neutrinomax=initm.Mag();}
      if (initm.Mag()<neutrinomin) {neutrinomin=initm.Mag();}
      }
    else if (k->second.GetPDGEncoding()==2112) {
      h9a->Fill(initm.Mag());
      h9b->Fill(finm.Mag());
      H9a->Fill(initp[0],initp[2]);
      H9b->Fill(finp[0],finp[2]);
      if (initm.Mag()>neutronmax) {neutronmax=initm.Mag();}
      if (initm.Mag()<neutronmin) {neutronmin=initm.Mag();}
      }
    else if (k->second.GetPDGEncoding()==2212) {
      h10a->Fill(initm.Mag());
      h10b->Fill(finm.Mag());
      H10a->Fill(initp[0],initp[2]);
      H10b->Fill(finp[0],finp[2]);
      if (initm.Mag()>protonmax) {protonmax=initm.Mag();}
      if (initm.Mag()<protonmin) {protonmin=initm.Mag();}
      }
    else {
      h11a->Fill(initm.Mag());
      h11b->Fill(finm.Mag());
      H11a->Fill(initp[0],initp[2]);
      H11b->Fill(finp[0],finp[2]);
      if (initm.Mag()>othermax) {othermax=initm.Mag();}
      if (initm.Mag()<othermin) {othermin=initm.Mag();}
      if (otherbyevent) {
        std::cout<<(*traj).GetParticleName()<<" PDG Encoding= "<<(*traj).GetPDGEncoding();
        std::cout<<" From track number "<<(*traj).GetTrackId()<<std::endl;
        }
      }
    }	
  }
  }
return true;
}
// Looking at the hits container
//		COMET::IHandle<COMET::IG4HitContainer> cdcg4hits = event.Get<COMET::IG4HitContainer>("name");
//	if (!cdcg4hits){
//		COMETLog("ICDCDetResp: GenerateHits: No hits found ");
//		return true;
//		}
//		COMETLog("ICDCDetResp: GenerateHits: Number of G4 hits found " << cdcg4hits->size());
// the Hits container
//		std::string containerTitle("CDCDetResponse hits for the CDC detector");
//		COMET::IHitSelection *hits = new COMET::IHitSelection("CDC",containerTitle.c_str());
// the digits container
//		COMET::IDigitContainer* digits = new COMET::IDigitContainer("cdc");
// detector response simulation
// This must be static so we delete the old event data at the beginning of the next event.
//		static std::map<COMET::IChannelId,std::vector<IFreeCarrier*>> freeCarrierList;
//		std::cout << "**** Size of freeCarrierList: " << freeCarrierList.size() << std::endl;
void Finalize(COMET::ICOMETOutput * const output) {


//all commented out canvases are redundant; they show that the particles are stopped (momentum=0) at their final positions

    Histfile->cd();
    Histfile->Write();
    Histfile->Close();

    if (showcounts) {
      for (std::map<std::string,Int_t>::iterator n=counts.begin();n!=counts.end();n++) {
        std::cout<<"Number of "<<n->first<<"\t\t"<<n->second<<std::endl;
        numall+=n->second;
      }
    std::cout<<"Total number of particles= "<<numall<<std::endl;
    }
    std::cout<<"Maximum gamma momentum (MeV/c)=\t\t"<<gammamax<<"\t\tMinimum gamma momentum (MeV/c)=\t\t"<<gammamin<<std::endl;
    std::cout<<"Maximum electron momentum (MeV/c)=\t"<<electronmax<<"\t\tMinimum electron momentum (MeV/c)=\t"<<electronmin<<std::endl;
    std::cout<<"Maximum positron momentum (MeV/c)=\t"<<positronmax<<"\t\tMinimum positron momentum (MeV/c)=\t"<<positronmin<<std::endl;
    std::cout<<"Maximum muon momentum (MeV/c)=\t\t"<<muonmax<<"\t\tMinimum muon momentum (MeV/c)=\t\t"<<muonmin<<std::endl;
    std::cout<<"Maximum antimuon momentum (MeV/c)=\t"<<antimuonmax<<"\t\tMinimum antimuon momentum (MeV/c)=\t"<<antimuonmin<<std::endl;
    std::cout<<"Maximum pion+ momentum (MeV/c)=\t\t"<<pionpmax<<"\t\tMinimum pion+ momentum (MeV/c)=\t\t"<<pionpmin<<std::endl;
    std::cout<<"Maximum pion- momentum (MeV/c)=\t\t"<<pionnmax<<"\t\tMinimum pion- momentum (MeV/c)=\t\t"<<pionnmin<<std::endl;
    std::cout<<"Maximum neutrino momentum (MeV/c)=\t"<<neutrinomax<<"\t\tMinimum neutrino momentum (MeV/c)=\t"<<neutrinomin<<std::endl;
    std::cout<<"Maximum neutron momentum (MeV/c)=\t"<<neutronmax<<"\t\tMinimum neutron momentum (MeV/c)=\t"<<neutronmin<<std::endl;
    std::cout<<"Maximum proton momentum (MeV/c)=\t"<<protonmax<<"\t\tMinimum proton momentum (MeV/c)=\t"<<protonmin<<std::endl;
    std::cout<<"Maximum other momentum (MeV/c)=\t\t"<<othermax<<"\t\tMinimum other momentum (MeV/c)=\t\t"<<othermin<<std::endl;
    std::cout<<"Number of trajs null handles= "<<numnull01<<std::endl;
    std::cout<<"Number of PVC null handles= "<<numnull02<<std::endl;
    std::cout<<"Number of vertices= "<<numvertex<<std::endl;
    std::cout<<"Maximum primary gamma momentum (MeV/c)=\t\t"<<ppgammamax<<"\t\tMinimum primary gamma momentum (MeV/c)=\t\t"<<ppgammamin<<std::endl;
    std::cout<<"Maximum primary electron momentum (MeV/c)=\t"<<ppelectronmax<<"\t\tMinimum primary electron momentum (MeV/c)=\t"<<ppelectronmin<<std::endl;
    std::cout<<"Maximum primary positron momentum (MeV/c)=\t"<<pppositronmax<<"\t\tMinimum primary positron momentum (MeV/c)=\t"<<pppositronmin<<std::endl;
    std::cout<<"Maximum primary muon momentum (MeV/c)=\t\t"<<ppmuonmax<<"\t\tMinimum primary muon momentum (MeV/c)=\t\t"<<ppmuonmin<<std::endl;
    std::cout<<"Maximum primary antimuon momentum (MeV/c)=\t"<<ppantimuonmax<<"\t\tMinimum prmary antimuon momentum (MeV/c)=\t"<<ppantimuonmin<<std::endl;
    std::cout<<"Maximum primary pion+ momentum (MeV/c)=\t\t"<<pppionpmax<<"\t\tMinimum primary pion+ momentum (MeV/c)=\t\t"<<pppionpmin<<std::endl;
    std::cout<<"Maximum primary pion- momentum (MeV/c)=\t\t"<<pppionnmax<<"\t\tMinimum primary pion- momentum (MeV/c)=\t\t"<<pppionnmin<<std::endl;
    std::cout<<"Maximum primary neutrino momentum (MeV/c)=\t"<<ppneutrinomax<<"\t\tMinimum primary neutrino momentum (MeV/c)=\t"<<ppneutrinomin<<std::endl;
    std::cout<<"Maximum primary neutron momentum (MeV/c)=\t"<<ppneutronmax<<"\t\tMinimum primary neutron momentum (MeV/c)=\t"<<ppneutronmin<<std::endl;
    std::cout<<"Maximum primary proton momentum (MeV/c)=\t"<<ppprotonmax<<"\t\tMinimum primary proton momentum (MeV/c)=\t"<<ppprotonmin<<std::endl;
    std::cout<<"Maximum primary other momentum (MeV/c)=\t\t"<<ppothermax<<"\t\tMinimum primary other momentum (MeV/c)=\t\t"<<ppothermin<<std::endl;
    }


private:
  std::string fLSOption;
  bool fQuiet;
  bool fPrintList;
  bool showcounts;
  bool lsoption;
  bool trajcount;
  bool otherbyevent;

  TH1F *h1a;
  TH1F *h1b;
  TH2F *H1a;
  TH2F *H1b;
  TH1F *h2a;
  TH1F *h2b;
  TH2F *H2a;
  TH2F *H2b;
  TH1F *h3a;
  TH1F *h3b;
  TH2F *H3a;
  TH2F *H3b;
  TH1F *h4a;
  TH1F *h4b;
  TH2F *H4a;
  TH2F *H4aa;
  TH2F *H4b;
  TH2F *H4ba;
  TH1F *h5a;
  TH1F *h5b;
  TH2F *H5a;
  TH2F *H5b;
  TH1F *h6a;
  TH1F *h6b;
  TH2F *H6a;
  TH2F *H6b;
  TH1F *h7a;
  TH1F *h7b;
  TH2F *H7a;
  TH2F *H7b;
  TH1F *h8a;
  TH1F *h8aa;
  TH1F *h8ab;
  TH1F *h8ac;
  TH1F *h8ad;
  TH1F *h8ae;
  TH1F *h8af;
  TH1F *h8b;
  TH2F *H8a;
  TH2F *H8aa;
  TH2F *H8ac;
  TH2F *H8ad;
  TH2F *H8ae;
  TH2F *H8af;
  TH2F *H8b;
  TH2F *H8ba;
  TH1F *h9a;
  TH1F *h9b;
  TH2F *H9a;
  TH2F *H9b;
  TH1F *h10a;
  TH1F *h10b;
  TH2F *H10a;
  TH2F *H10b;
  TH1F *h11a;
  TH1F *h11b;
  TH2F *H11a;
  TH2F *H11b;
  TH3F *H12a;
  TH3F *H12b;
  TH3F *H12c;
  TH3F *H12d;
  TH1F *h13a;
  TH1F *h13b;
  TH1F *h13c;
  TH1F *h13d;
  TH1F *h13e;
  TH1F *h13f;
  TH1F *h13g;
  TH1F *h13h;
  TH1F *h13i;
  TH1F *h13j;
  TH1F *h13k;
  TH3F *G1a;
  TH3F *G1b;
  TH3F *G1c;
  TH3F *G1d;
  TH3F *G1e;
  TH3F *G1f;
  TH3F *G1g;
  TH3F *G1h;
  TH3F *G1i;
  TH3F *G1j;
  TH3F *G1k;

  TFile *Histfile;

  Int_t	numall;
  Int_t	numnull01;
  Int_t	numnull02;
  Int_t	numvertex;

  Float_t	gammamin;
  Float_t	electronmin;
  Float_t	positronmin;
  Float_t	muonmin;
  Float_t	antimuonmin;
  Float_t	pionpmin;
  Float_t	pionnmin;
  Float_t	neutrinomin;
  Float_t	neutronmin;
  Float_t	protonmin;
  Float_t	othermin;

  Float_t	gammamax;
  Float_t	electronmax;
  Float_t	positronmax;
  Float_t	muonmax;
  Float_t	antimuonmax;
  Float_t	pionpmax;
  Float_t	pionnmax;
  Float_t	neutrinomax;
  Float_t	neutronmax;
  Float_t	protonmax;
  Float_t	othermax;
  
  Float_t	ppgammamin;
  Float_t	ppelectronmin;
  Float_t	pppositronmin;
  Float_t	ppmuonmin;
  Float_t	ppantimuonmin;
  Float_t	pppionpmin;
  Float_t	pppionnmin;
  Float_t	ppneutrinomin;
  Float_t	ppneutronmin;
  Float_t	ppprotonmin;
  Float_t	ppothermin;

  Float_t	ppgammamax;
  Float_t	ppelectronmax;
  Float_t	pppositronmax;
  Float_t	ppmuonmax;
  Float_t	ppantimuonmax;
  Float_t	pppionpmax;
  Float_t	pppionnmax;
  Float_t	ppneutrinomax;
  Float_t	ppneutronmax;
  Float_t	ppprotonmax;
  Float_t	ppothermax;

  std::map<std::string,Int_t> counts;
};

int main(int argc, char **argv) {
    IDumpEvent userCode;
    return COMET::cometEventLoop(argc,argv,userCode,0);
}