BremShowerFilter_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
/// \file    BremShowerFilter.cxx
/// \author  H. Duyang 12/18/13
////////////////////////////////////////////////////////////////////////
// ART includes
#include "art/Framework/Principal/Event.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/FindManyP.h"
#include "canvas/Persistency/Common/Ptr.h" 
#include "canvas/Persistency/Common/PtrVector.h"
#include "fhiclcpp/ParameterSet.h"  
#include "messagefacility/MessageLogger/MessageLogger.h" 
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Core/EDFilter.h"

// NOvA includes
#include "Geometry/Geometry.h"
//#include "Geometry/Geo.h"
#include "GeometryObjects/CellGeo.h"
#include "Calibrator/Calibrator.h"
#include "MCCheater/BackTracker.h"
#include "Simulation/ParticleNavigator.h"
#include "nusimdata/SimulationBase/MCParticle.h"
#include "RecoBase/Track.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/RecoHit.h"
#include "RawData/DAQHeader.h"
#include "RawData/RawDigit.h"
#include "Utilities/func/MathUtil.h"
#include "NovaDAQConventions/DAQConventions.h"
#include "DAQDataFormats/NanoSliceVersionConvention.h"

//root includes
#include <vector>
#include "TH1F.h"
#include "TH2F.h"
#include "TH3F.h"
#include "TTree.h"
#include "TLorentzVector.h"

namespace bsf {
  
  class BremShowerFilter : public art::EDFilter {
    
  public:
    explicit BremShowerFilter(fhicl::ParameterSet const& pset);
    virtual ~BremShowerFilter();
    
    virtual bool filter(art::Event & evt);
    void reconfigure(const fhicl::ParameterSet& pset);
    void beginJob();
    void findShowerByReco(art::Ptr<rb::Track> &trk, int &startPlane, int &endPlane);
    void findShowerByTruth(art::Ptr<rb::Track> &trk, int &startPlane, int &endPlane);
    bool inFiducial(double x, double y, double z);
    bool trackStartContained(art::Ptr<rb::Track> &trk);    
    bool trackEndContained(art::Ptr<rb::Track> &trk);    

  private:
    bool dokeep;
    //parameters here

    const art::ProductToken<std::vector<rb::Track>> fTrackToken;
    const art::ProductToken<std::vector<rb::Cluster>> fClusterToken;
    const art::ProductToken<std::vector<rb::CellHit>> fCellHitToken;

    std::string fTrackLabel;   ///< label of module making fls hits

    bool fWriteParentSliceMode;
    bool fApplyFilter;

    //find shower by truth or by reco
    bool fShowerByTruth;//if true, look for only DIF using truth
    bool fShowerByReco;//if true look for all shower using reco
    bool fFindDiF;//if true, look for only DIF using reco
    float fAngleCut;//cosZ cut on muon track
    float fPlaneCut;//how many plane out of all planes we require?
    float fEnergyCut;//requre a minimum shower energy    
    float fFiducialX;
    float fFiducialY;
    float fFiducialZ;
    int fPlaneMin;
    int fPlaneMax;

    //make a tree to save variables for analysis
    TTree   *fTree;
    
    int run;
    int subrun;
    int event; 
    int particleId;//this is the true particle id
    int trackId;//this is the reco track id
    //int process;
    //int decayPlane;
    int startPlane;
    int endPlane;
    //int nplane;
    double cosz;
    double eshower;
    double planeDedx[2];
    double planeDedxMR[2];
    double cellDedx[2];
    double cellDedxMR[2];
    int    fDCSDistance;


    //double E_mu_decay;
    //double ADCTrack[1000][500];
    //double edepTrack[1000][500];
    //double ADCSlice[1000][500];
    //double edepSlice[1000][500];
    //double ADCShower[1000][500];
    //double edepShower[1000][500];
    //double edepEl[1000][500];
    //double edepMu[1000][500];

    //TLorentzVector *v4;         // particle vertex
    //TLorentzVector *p4Mu;         // particle 4 momentum
    //TLorentzVector *p4El;   // mother 4 momentum 
    //TVector3 *dir;
    TVector3 *vtx= new TVector3(0,0,0);
    
    //hists here
    //TH1F *fHEmu[5];
    //TH1F *fHEndPlane;
    //TH1F *fHEndCell;
    TH1F *fHNTrk;
    TH1F *fHCosZ;
    TH1F *fHCosZCut[5];
    TH1F *fHStartContainment;
    TH1F *fHEndContainment;
    TH1F *fHNPlane;
    TH1F *fHEShower;
    //TH3F *fHTrackStart;
    //TH3F *fHTrackEnd;
    TH1F *fHTrackStartZ;
    //TH2F *fHEndPoint;
    TH1F *fHEdep;
    TH1F *fHADC;
    TH1F *fHPE;
    TH1F *fHPECorr;
    
    TH2F *fHHitMap;
    TH2F *fHHitMapRW;

  };
}

namespace bsf {
  //......................................................................
  
  BremShowerFilter::BremShowerFilter(fhicl::ParameterSet const& pset)
    : fTrackToken(consumes<std::vector<rb::Track>>(pset.get<std::string>("TrackLabel"))),
      fClusterToken(consumes<std::vector<rb::Cluster>>(pset.get<std::string>("ClusterLabel"))),
      fCellHitToken(consumes<std::vector<rb::CellHit>>(pset.get<std::string>("CellHitLabel")))
    //v4(), p4Mu(), p4El()
  {
    reconfigure(pset);
    //produces<std::vector<dm::BremShowerInfo> >();
    produces< std::vector<rb::Cluster> >();
    produces<std::vector<rawdata::RawDigit> >();    
    
    //produces< std::vector<rb::Track> >();
  }
  
  //......................................................................
  
  BremShowerFilter::~BremShowerFilter() { }
  
  //......................................................................
  
  
  void BremShowerFilter::beginJob(){
    art::ServiceHandle<art::TFileService> tfs;
    
    fTree = tfs->make<TTree>("brem", "brem");
    fTree->Branch("run", &run, "run/I");
    fTree->Branch("subrun", &subrun, "subrun/I");
    fTree->Branch("event", &event, "event/I");
    fTree->Branch("particleId", &particleId, "particleId/I");
    fTree->Branch("trackId", &trackId, "trackId/I");
    //fTree->Branch("decayPlane", &decayPlane, "decayPlane/I");
    //fTree->Branch("process", &process, "process/I");
    fTree->Branch("startPlane", &startPlane, "startPlane/I");
    fTree->Branch("endPlane", &endPlane, "endPlane/I");
    //fTree->Branch("nPlane", &nPlane, "nPlane/I");
    fTree->Branch("cosz", &cosz, "cosz/D");
    fTree->Branch("eshower", &eshower, "eshower/D");
    fTree->Branch("vtx", "TVector3", &vtx);
    fTree->Branch("planeDedx", &planeDedx, "planeDedx[2]/D");
    fTree->Branch("planeDedxMR", &planeDedxMR, "planeDedxMR[2]/D");
    fTree->Branch("cellDedx", &cellDedx, "cellDedx[2]/D");
    fTree->Branch("cellDedxMR", &cellDedxMR, "cellDedxMR[2]/D");
    //fTree->Branch("E_mu_decay", &E_mu_decay, "E_mu_decay/D");
    //fTree->Branch("v4", "TLorentzVector", &v4);
    //fTree->Branch("p4El", "TLorentzVector", &p4El);
    //fTree->Branch("p4Mu", "TLorentzVector", &p4Mu);
    //fTree->Branch("dir", "TVector3", &dir);

    //fTree->Branch("edepMu", &edepMu, "edepMu[1000][500]/D");
    //fTree->Branch("edepEl", &edepEl, "edepEl[1000][500]/D");
    //fTree->Branch("edepTrack", &edepTrack, "edepTrack[1000][500]/D");
    //fTrackTree->Branch("ADCTrack", &ADCTrack, "ADCTrack[1000][500]/D");
    //fTree->Branch("edepSlice", &edepSlice, "edepSlice[1000][500]/D");
    //fTrackTree->Branch("ADCSlice", &ADCSlice, "ADCSlice[1000][500]/D");
    //fTree->Branch("edepShower", &edepShower, "edepShower[1000][500]/D");
    //fTrackTree->Branch("ADCShower", &ADCShower, "ADCShower[1000][500]/D");
    
    fHNTrk = tfs->make<TH1F>("hNTrk", "", 100, 0, 100);

    fHCosZ = tfs->make<TH1F>("hCosZ", "", 200, -1, 1);
    for(int i=0; i<5; i++)
      fHCosZCut[i] = tfs->make<TH1F>(Form("hCosZCut%d", i), "", 200, -1, 1);
    fHNPlane = tfs->make<TH1F>("hNPlane", "", 1000, 0, 1000);                                                                                                  
    fHStartContainment = tfs->make<TH1F>("hStartContainment", "", 2, 0, 2);    
    fHEndContainment = tfs->make<TH1F>("hEndContainment", "", 2, 0, 2);    
    fHEShower = tfs->make<TH1F>("hEShower", "", 1000, 0, 10);    
    fHEdep = tfs->make<TH1F>("hEdep", "", 1000, 0, 1);    
    fHADC = tfs->make<TH1F>("hADC", "", 1000, 0, 1000);    
    fHPE = tfs->make<TH1F>("hPE", "", 1000, 0, 1000);    
    fHPECorr = tfs->make<TH1F>("hPECorr", "", 1000, 0, 1000);    
    //fHTrackStart = tfs->make<TH3F>("hTrackStart", "", 200, -1000, 1000, 200, -1000, 1000, 6000, 0, 6000);
    //fHTrackEnd = tfs->make<TH3F>("hTrackEnd", "", 200, -1000, 100, 200, -1000, 1000, 6000, 0, 6000);
    fHTrackStartZ = tfs->make<TH1F>("hTrackStartZ", "", 6000, 0, 6000);
    fHHitMap = tfs->make<TH2F>("hHitMap", "", 200, 0, 200, 400, 0, 400);
    fHHitMapRW = tfs->make<TH2F>("hHitMapRW", "", 200, 0, 200, 400, 0, 400);

  }
  
  //......................................................................
  
  void BremShowerFilter::reconfigure(const fhicl::ParameterSet& pset)
  {
    fTrackLabel        = pset.get< std::string >("TrackLabel");
    fWriteParentSliceMode = pset.get< bool     >("WriteParentSliceMode");
    fApplyFilter       = pset.get< bool        >("ApplyFilter");
    fShowerByTruth     = pset.get< bool        >("ShowerByTruth", "false");
    fShowerByReco      = pset.get< bool        >("ShowerByReco", "false");
    fFindDiF           = pset.get< bool        >("FindDiF", "false");
    fAngleCut          = pset.get< float       >("AngleCut");
    fPlaneCut          = pset.get< float       >("PlaneCut");
    fEnergyCut         = pset.get< float       >("EnergyCut");
    fFiducialX         = pset.get< float       >("FiducialX");
    fFiducialY         = pset.get< float       >("FiducialY");
    fFiducialZ         = pset.get< float       >("FiducialZ");
    fPlaneMin          = pset.get< int         >("PlaneMin");
    fPlaneMax          = pset.get< int         >("PlaneMax");
    fDCSDistance       = pset.get< int         >("DCSDistance");
  }
  
  //......................................................................
  bool BremShowerFilter::filter(art::Event & evt)
  { 
    art::ServiceHandle<geo::Geometry> geom;
    art::ServiceHandle<calib::Calibrator> cal;

    //Look at the truth
    //Get the particle navigator
    //art::ServiceHandle<cheat::BackTracker> bt;
    //const sim::ParticleNavigator &pnav = bt->ParticleNavigator();    

    //use reco base
    art::Handle< std::vector<rb::Track> > trackcol;
    evt.getByToken(fTrackToken, trackcol);
    
    art::Handle< std::vector<rb::CellHit> > cellhit;
    evt.getByToken(fCellHitToken, cellhit);
    
    // Get the slices in the event
    art::Handle< std::vector<rb::Cluster> > slicecol;
    evt.getByToken(fClusterToken, slicecol);
    
    art::FindManyP<rb::Track> fmt(slicecol, evt, fTrackLabel);
    
    std::vector< rawdata::RawDigit> newDigits;
    newDigits.reserve(cellhit->size());
    for(const rawdata::RawDigit& digit: *cellhit) newDigits.push_back(digit);
    
    //std::unique_ptr<std::vector<dm::BremShowerInfo> > 
    //BremShower(new std::vector<dm::BremShowerInfo>);
    
    //and the rawdigit for the showers
    std::unique_ptr< std::vector<rawdata::RawDigit> > showerDigits(new std::vector<rawdata::RawDigit>);
    
    //a new cluster to hold the showers
    std::unique_ptr< std::vector<rb::Cluster> > showercol(new std::vector<rb::Cluster>);
    //std::unique_ptr< std::vector<rb::Track> > etrackcol(new std::vector<rb::Track>);

    dokeep = false;    //indicate to save the event or not

    run    = evt.run();
    subrun = evt.subRun();
    event  = evt.id().event();    

    fHNTrk->Fill(trackcol->size());
   
    //prepare a vector of hits that needs to be removed
    std::vector<rb::WeightedHit > trackHits;//while I call this track hits, it is actually all hits except those belong to ther shower.

    //loop over slices here   
    for(unsigned int iSlice = 0; iSlice < slicecol->size(); ++iSlice) {

      art::Ptr<rb::Cluster> slc(slicecol,iSlice);

      bool goodSlice = true;
      if (slc->NCell() <= 0) continue;
      if(slc->MeanTNS () <25000 || slc->MeanTNS () > 475000) {goodSlice = false;} // rejecting the trigger window on edges 

      std::vector< art::Ptr<rb::Track> > tracks = fmt.at(iSlice);
      
      //is there a shower on this track?
      //if yes, where it starts and where it ends?
      //assume there is no more than one shower per track for now
      
      startPlane = -1;
      endPlane = -1;
      trackId = iSlice; //what saving here is actually the slice where the shower from

      bool isShower = false;
      bool isCandTrack = false;
      //bool isBrem = false;
      bool isDiF = false;

      if(!slc->IsNoise() && tracks.size()==1 && tracks[0]->Is3D() && goodSlice){//if there is no track, then no shower at all
        
        bool startContained = false;
        bool endContained = false;

        //first decide is this track an candidate track?
        if(int(tracks[0]->MaxPlane()) > fPlaneMin && int(tracks[0]->MinPlane()) < fPlaneMax){//well this means the track passes the active area of the detector
          startContained = trackStartContained(tracks[0]);
          endContained = trackEndContained(tracks[0]);
          
          fHStartContainment->Fill(startContained);
          fHEndContainment->Fill(endContained);

          //fHTrackStart->Fill(tracks[0]->Start().X(), tracks[0]->Start().Y(), tracks[0]->Start().Z());
          fHTrackStartZ->Fill(tracks[0]->Start().Z());

          int nPlaneTrack = 0;
          
          if(int(tracks[0]->MaxPlane())>fPlaneMax) nPlaneTrack = fPlaneMax - int(tracks[0]->MinPlane());
          else if(int(tracks[0]->MinPlane())<fPlaneMin) nPlaneTrack = int(tracks[0]->MaxPlane()) - fPlaneMin;
          else nPlaneTrack = tracks[0]->MaxPlane() - tracks[0]->MinPlane();              
          
          cosz = tracks[0]->Dir().Z();
          
          if(startContained) isCandTrack = false;//#1 cut on containment
          else{
            fHCosZ->Fill(tracks[0]->Dir().Z());           
            //if(fabs(tracks[0]->Dir().Z())<fAngleCut) isCandTrack = false;//#2 cut on angle    
            if(tracks[0]->Dir().Z()<fAngleCut) isCandTrack = false;//#2 cut on angle, keep only downstream tracks       
            else {
              fHNPlane->Fill(nPlaneTrack);  //#3 cut on nPlane
              if(nPlaneTrack<fPlaneCut) isCandTrack = false;        
              else{
                isCandTrack = true;
                fHCosZCut[0]->Fill(tracks[0]->Dir().Z());         
              }
            }
          }
          
        }//end if track passes the activ area of the detector
        
        if(isCandTrack){//we found the candidate track!
          if(fShowerByReco){//means to find both brem and DiF
            
            findShowerByReco(tracks[0], startPlane, endPlane);    
            if(startPlane>0 && endPlane>0){
              fHCosZCut[1]->Fill(tracks[0]->Dir().Z());           
              if(startPlane>fPlaneMin && endPlane<fPlaneMax){//select shower contained in the active area of the detector
                fHCosZCut[2]->Fill(tracks[0]->Dir().Z());         
                isShower = true; 
                //if(endContained) 
                //isDiF = true;
              }

            }
          }//end if fShowerByReco
          
          //if(fFindDiF){//means to find only DiF
          //findShowerByReco(tracks[0], startPlane, endPlane);    
          //if(!startContained && endContained && startPlane>0 && endPlane<0) {isShower = true; isDiF = true;}
          //}
          
          if(fShowerByTruth){
            findShowerByTruth(tracks[0], startPlane, endPlane); 
            if(!startContained && endContained && startPlane>0 && endPlane>0) {isShower = true; isDiF = true;}
          }
        }//end if is candidate track      

      }//end if track size is 1
      
      //the following part fills a vector of weighted hits trackHits to be removed     
      //consider to use a function for this
      //getTrackHits(tracks[0], startPlane, endPlane, trackHits)

      if(isShower){//we found the shower!
        //in this case, keep the hits in shower region, remove the rest in the slice
        //also save the shower to a new cluster object
        //bool muonRemove = false;
        //if(endPlane>0) muonRemove=true;//we call it brem. This criteria is not a good one of course
        
        //prepare an cluster object to save the shower hits                                                           
        rb::Cluster shower;     

        //shower region should be taken care of first because we want to remove some hit while keep others
        //we need a loop over planes in shower region! 
        
        int minPlane = startPlane;
        int maxPlane;
        if(endPlane>startPlane) maxPlane = endPlane;
        else maxPlane = tracks[0]->MaxPlane();
        
        
        TVector3 startPoint;
        geom->Plane( startPlane )->Cell( 10 )->GetCenter(startPoint); 

        //try to define a shower vertex here
        //loop over trajectory points, find the one on the startPlane, call it the vertex
        //note that it is only true for a down stream shower. 
        int nTrajPt = tracks[0]->NTrajectoryPoints();
        
        for( int iTraj = 0 ; iTraj < nTrajPt ; iTraj++ ){
          //this is to find the vertex
          if( startPoint.Z()-3.3 < tracks[0]->TrajectoryPoint(iTraj).Z() && 
              startPoint.Z()+3.3 > tracks[0]->TrajectoryPoint(iTraj).Z()){
            vtx->SetXYZ(tracks[0]->TrajectoryPoint(iTraj).X(),
                        tracks[0]->TrajectoryPoint(iTraj).Y(),
                        tracks[0]->TrajectoryPoint(iTraj).Z());

          }

          
        }//end loop over trajectory points
        
        //loop over planes in shower region
        for(int iPlane=minPlane; iPlane<=maxPlane; iPlane++){
          std::vector< art::Ptr<rb::CellHit> > planeHits;//prepare a vector of cell hits on this plane
          
          //the cell on this plane closest, and the 2nd closest to track, 
          unsigned int cellClose1=0;
          unsigned int cellClose2=0;
          double cellCloseDist1=10000000;
          double cellCloseDist2=10000000;
          double edepClose1 = 0;
          double edepClose2 = 0;

          //int trajpoint = -1;

          //loop over trajectory points
          //for( int iTraj = 1; iTraj < nTrajPt ; iTraj++){
            
            //const geo::CellUniqueId cid = geom->CellId(tracks[0]->TrajectoryPoint(iTraj).X(), 
            //                                 tracks[0]->TrajectoryPoint(iTraj).Y(), 
            //                                 tracks[0]->TrajectoryPoint(iTraj).Z(),
            //                                 1,1,0, 2);
            //int trajPlane = -1;
            
            //geom->IdToPlane( cid, &trajPlane);

            //if(trajPlane==iPlane){
            //trajpoint = iTraj;
            //break;
            //}
          //}//end loop over trajectory points
          
          //double trajXY = 0;
          //if(trajpoint>0){
          //if(geom->Plane(iPlane)->View() == geo::kX)
          //  trajXY = tracks[0]->TrajectoryPoint(trajpoint).X();
          //else 
          //  trajXY = tracks[0]->TrajectoryPoint(trajpoint).Y();
          //}
          //else continue;//not trajpoint found in this plane; use the second method. (I don't think this will happen)
            
          //this is the energy a muon *should* deposit on this plane
          double mipDeDx  = 0.00157;
          //double mipDeDx  = 0.00135;
          double mipPlane = 2*mipDeDx*geom->Plane(iPlane)->Cell(0)->HalfD()/fabs(tracks[0]->Dir().Z());
          
          //then loop over track hits, for those in shower region, fill a vector with hits on this plane
          //
          for(unsigned int iCell=0; iCell<tracks[0]->NCell(); ++iCell){
            const rb::RecoHit rhit = tracks[0]->RecoHit(tracks[0]->Cell(iCell));
            const rb::CellHit *chit = tracks[0]->Cell(iCell).get();
            
            if(!rhit.IsCalibrated()) continue;
            if(chit->Plane() != iPlane) continue;
            planeHits.push_back(tracks[0]->Cell(iCell));
            
            //decide what is the distance of the hits to track, if too far away, just remove it(will do)
            double xyz[3];
            geom->Plane( chit->Plane() )->Cell( chit->Cell() )->GetCenter(xyz);                                                                                
            
            //well this is a stupid way of calculating the distance
            double cellXY, trackStartXY;
            double cellZ = xyz[2];
            double trackStartZ = tracks[0]->Start().Z();
            double slope = tracks[0]->Dir().X()/tracks[0]->Dir().Z();
            
            if(geom->Plane(chit->Plane())->View() == geo::kX){
              cellXY = xyz[0];
              trackStartXY = tracks[0]->Start().X();
              slope = tracks[0]->Dir().X()/tracks[0]->Dir().Z();
            }
            else{
              cellXY = xyz[1];
              trackStartXY = tracks[0]->Start().Y();
              slope = tracks[0]->Dir().Y()/tracks[0]->Dir().Z();
            }
            
            double cellDist = fabs(cellXY - trackStartXY - slope*(cellZ - trackStartZ));
            //double cellDist1 = fabs(cellXY - trajXY);
            //sort plane hits according to the dist form track.

            if(cellDist<cellCloseDist1){
              cellCloseDist2 = cellCloseDist1;
              cellClose2 = cellClose1;
              edepClose2 = edepClose1;
              cellCloseDist1 = cellDist;
              cellClose1 = chit->Cell();
              edepClose1 = rhit.GeV();
            }
            else if(cellDist<cellCloseDist2){
              cellCloseDist2 = cellDist;
              cellClose2 = chit->Cell();
              edepClose2 = rhit.GeV();      
            }
          }//end loop over track hits
          
          //loop over plane hits... try to define a weight and fill weighted hits 
          for(unsigned int iHit = 0; iHit < planeHits.size(); ++iHit){      
            
            art::Ptr<rb::CellHit> chit = planeHits[iHit];
            //rb::RecoHit rhit = tracks[0]->RecoHit(chit->Cell());
            
            double weight = 0;//the default should be zero, mean don't remove 
            if(!isDiF){//only for brem we need to remove the muon in shower region
              
              if(chit->Cell()==cellClose1){//the cell closest to track  
                if(edepClose1>mipPlane){//if this cell alone has energy > mip
                  //reset the energy to be a mip
                  weight = mipPlane/edepClose1;
                  if(weight<0 || weight>1) {
                    mf::LogWarning("BremShowerFilter") <<"weight = "<<weight<<" "<<mipPlane<<" "<<edepClose1<<" "<<edepClose2;
                  }
                }
                else //if the cell has energy <= mip 
                  weight = 1;//fill this one with weight=1, and remove it later
              }
              else if(chit->Cell()==cellClose2){//the 2nd closest cell  
                //if the sum of the first and the second hits larger than a mip, reset the energy of the second to be what left
                if(edepClose1<mipPlane){
                  if(edepClose1+edepClose2>mipPlane){
                    weight = (-edepClose1 + mipPlane)/edepClose2;
                    if(weight<0 || weight>1) {
                      mf::LogWarning("BremShowerFilter") <<"weight = "<<weight<<" "<<mipPlane<<" "<<edepClose1<<" "<<edepClose2;
                    }
                  }
                  else //if the two closest cells still don't have enough energy to make the mip, remove them all. leave the rest.
                    weight = 1;
                  //otherwise do nothing, weight = 1
                }
                // else continue;//if the closest aleady has > mip energy, don't touch the second
              }
              else weight=0;//keep all other cell hits, means weight should be 0!
              //else continue;
            }//endl if !isDiF
            else weight = 0;//if it is dif then don't remove anything

            if(weight<0 || weight>1) {
              mf::LogWarning("BremShowerFilter") <<"weight = "<<weight;
            }
            
            if(weight>0){//save trackhits
              rb::WeightedHit newHit(chit, weight);
              trackHits.push_back(newHit);
            }
            else if((1-weight)>0)//save shower hits
              shower.Add(chit, (1-weight));

          }//endl loop over plane hits      
        }//end loop over planes 
        
        //loop over track hits again, this time for removing those outside shower region (I know this method sounds stupid) 
        //4/11/2014: change to slice hits
        for(unsigned int iCell=0; iCell<slc->NCell(); ++iCell){
          //art::Ptr<rb::CellHit> chit = slices[i]->Cell(iCell);                                                                                                 
          const rb::CellHit *chit = slc->Cell(iCell).get();
          int plane = chit->Plane();     
          
          if(plane<minPlane || plane>maxPlane){//not in shower region, simplay remve the hits       
            rb::WeightedHit newHit(slc->Cell(iCell), 1);
            trackHits.push_back(newHit);            
          }

        }//end loop other track hits

        if(!fWriteParentSliceMode) showercol->push_back(shower);

        shower.Clear();
        dokeep  = true;//keep this event as long as there is shower
        fTree->Fill();

        if(fWriteParentSliceMode) showercol->push_back(*slc);
      }//end if shower found
      else{//no shower found. remove all hits in this slice
        //loop over slice hits. For now we just use the slice with the same index
        //if(iTrack<slicecol->size())
        for(unsigned int iCell=0; iCell<slc->NCell(); ++iCell){
          //art::Ptr<rb::CellHit> chit = slices[i]->Cell(iCell);
          //const rb::RecoHit rhit = slices[i]->RecoHit(slices[j]->Cell(iCell));
          rb::WeightedHit newHit(slc->Cell(iCell), 1);
          trackHits.push_back(newHit);
        }//endl loop over slice hits
        
      }//end else (no shower found)
      
    }//end loop over slices      

    //maybe it is faster to remove the trackhits all together here, instead in each track loop
    //-----------------------------------------------------------------------
    //prepare all digits in this event
    //now we want to clean up the track hits
    //if(!RemoveMuon) 

    if(dokeep){//only when the event has shower we want to remove the non-shower digits. otherwise the event will not ke saved anyway
    std::map<std::tuple<int, int, int>, rb::WeightedHit> trkMap;
    for(int iHit = 0; iHit<(int)trackHits.size(); iHit++ ){
      const rb::WeightedHit& h = trackHits[iHit];
      trkMap.emplace(std::make_tuple(h.hit->Plane(), h.hit->Cell(), h.hit->TDC()), h);
    }

    for(int iDigit = 0; iDigit<(int)newDigits.size(); iDigit++ ){
     
      const rawdata::RawDigit& rd = newDigits[iDigit];
      
      const rb::CellHit cHit =  cal->MakeCellHit( &rd );

      auto it = trkMap.find(std::make_tuple(cHit.Plane(), cHit.Cell(), cHit.TDC()));
      if(it == trkMap.end()){
        showerDigits->push_back(newDigits[iDigit]);
        continue;
      }

      const float muonWeight = it->second.weight;

          if( muonWeight == 1){
            // Nothing, hit is not used in output
          }
          else{ 
          //Mulitpoint readout begins. Copied from MRCC
          //  newDigits[iDigit].SetADC(0, (1-muonWeight)*cHit.ADC(0) );

            if( newDigits[iDigit].IsMC() && newDigits[iDigit].Version() == 0 && newDigits[iDigit].NADC() == 3){
              newDigits[iDigit].SetADC(0, (1-muonWeight)*cHit.ADC(0) );
              newDigits[iDigit].SetADC(1, (1-muonWeight)*cHit.ADC(1));
              newDigits[iDigit].SetADC(2, (1-muonWeight)*cHit.ADC(2));
            }
            else if(newDigits[iDigit].Version() == 0 )
              newDigits[iDigit].SetADC(0, (1-muonWeight)*cHit.ADC(0) );
            else{

              static daqdataformats::NanoSliceVersionConvention conv;
              unsigned int nPretrig = conv.getNPretriggeredSamples(newDigits[iDigit].Version());
              int16_t baseline = newDigits[iDigit].ADC(nPretrig - fDCSDistance );
              int nadc = newDigits[iDigit].NADC();
 
             for( int iadc = 0; iadc < nadc; ++iadc){
                int16_t newadc = ((newDigits[iDigit].ADC(iadc) - baseline )*( 1-muonWeight) ) + baseline;
                newDigits[iDigit].SetADC( iadc, newadc);
              }
            }

            showerDigits->push_back(newDigits[iDigit]);
          }
         //Mulitpoint readout ends.
         //   newDigits[iDigit].SetADC(0, (1-muonWeight)*cHit.ADC(0) );
         //   newDigits[iDigit].SetADC(1, (1-muonWeight)*cHit.ADC(1));
         //   newDigits[iDigit].SetADC(2, (1-muonWeight)*cHit.ADC(2));
    }// End loop over digits.            
    
    
    //evt.put(std::move(etrackcol));
    }//end if dokeep    
    
    evt.put(std::move(showerDigits));
    evt.put(std::move(showercol));

    if(fApplyFilter) return dokeep; else return true;
    
  }//end filter
  
}

bool bsf::BremShowerFilter::inFiducial(double x, double y, double z){
  art::ServiceHandle<geo::Geometry> geom;

  bool isContained = false;
  
  double xyzMax[3], xyzMin[3];

  //get the position of the boundary plane
  geom->Plane(fPlaneMax)->Cell(0)->GetCenter(xyzMax); 
  geom->Plane(fPlaneMin)->Cell(0)->GetCenter(xyzMin); 

  //now get the boundary
  double zMax = xyzMax[2] + geom->Plane(fPlaneMax)->Cell(0)->HalfD();
  double zMin = xyzMin[2] - geom->Plane(fPlaneMin)->Cell(0)->HalfD();

  if(x < (geom->DetHalfWidth() - fFiducialX) && x > (-geom->DetHalfWidth() + fFiducialX) && 
     y < (geom->DetHalfHeight() - fFiducialY) && y > (-geom->DetHalfHeight() + fFiducialY) && 
     //z < (geom->DetLength() - fFiducialZ) && z > fFiducialZ
     z < (zMax - fFiducialZ) && z > zMin + fFiducialZ
     )
    
    isContained = true;

  return isContained;
  
}

bool bsf::BremShowerFilter::trackStartContained(art::Ptr<rb::Track> &trk){

  bool isContained = inFiducial(trk->Start().X(), trk->Start().Y(), trk->Start().Z()); 

  return isContained;
}

bool bsf::BremShowerFilter::trackEndContained(art::Ptr<rb::Track> &trk){
  //art::ServiceHandle<geo::Geometry> geom;
  //loop over track hits, prepare eplane for future use                                                                                                                                                     
  double ymin = 10000;
  
  bool isContained = false;

  /*

  for(unsigned int iCell=0; iCell<trk->NCell(); ++iCell){

    const rb::CellHit *chit = trk->Cell(iCell).get();
    
    double xyz[3] = {0.};                                                                                                                                                                                   
    geom->Plane(chit->Plane())->Cell(chit->Cell())->GetCenter(xyz);                                                                                                                                      

    if(xyz[2]<ymin){
      isContained = geom->isInsideFiducialVolume(xyz[0], xyz[1], xyz[2]); 
      ymin = xyz[2];
    }
  }//end loop over track hits
  */
  
  int nTrajPt = trk->NTrajectoryPoints();

  for(int iTrajPt=0; iTrajPt<nTrajPt; iTrajPt++){
    double xTrajpt =  trk->TrajectoryPoint(iTrajPt).X();
    double yTrajpt =  trk->TrajectoryPoint(iTrajPt).Y();
    double zTrajpt =  trk->TrajectoryPoint(iTrajPt).Z();
   
    if(yTrajpt<ymin){
      isContained = inFiducial(xTrajpt, yTrajpt, zTrajpt); 
      ymin = yTrajpt;
    }
  }
  
  return isContained;
}

//find shower by reco info. the input is rb::Track, output is startPlane, endPlane. 
//if no shower, startPlane=endPlane=-1
void bsf::BremShowerFilter::findShowerByReco(art::Ptr<rb::Track> &trk, int &startPlane, int &endPlane){
  //for now, simply cut off tracks too short, or to steep.
 
  art::ServiceHandle<geo::Geometry> geom;

  int minPlane = trk->MinPlane();
  int maxPlane = trk->MaxPlane();
  //int nPlane = geom->NPlanes();
  //fHCosZ->Fill(trk->Dir().Z());  
  //if(fabs(trk->Dir().Z())<fAngleCut) return;  
  //fHNPlane->Fill(fabs(maxPlane-minPlane));  
  //if(fabs(maxPlane-minPlane)<fPlaneCut) return;


  if(minPlane<fPlaneMin) minPlane = fPlaneMin;
  if(maxPlane>fPlaneMax) maxPlane = fPlaneMax;

  std::map<unsigned int, double> ePlane;
  std::map<unsigned int, double> ePlaneAve;

  ePlane.clear();

  //loop over track hits, prepare eplane for future use
  for(unsigned int iCell=0; iCell<trk->NCell(); ++iCell){
    const rb::CellHit* cHit = trk->Cell(iCell).get();
    const rb::RecoHit rHit = trk->RecoHit(trk->Cell(iCell));

    if(cHit->Plane()<fPlaneMin || cHit->Plane()>fPlaneMax) continue;

    if( rHit.IsCalibrated() ){                                                                                                                             
      //ADCTrack[chit->Plane()][chit->Cell()] += chit->ADC(0);
      //edepTrack[chit->Plane()][chit->Cell()] += rhit.GeV();
      
      ePlane[cHit->Plane()] += rHit.GeV();
      fHEdep->Fill(rHit.GeV());
      fHADC->Fill(cHit->ADC(0));
      fHPE->Fill(cHit->PE());
      fHPECorr->Fill(rHit.PECorr());
    }
  }//end loop over cells
  
  //start a loop over track planes
  for(int iPlane=minPlane; iPlane<=maxPlane; iPlane++){

    if(iPlane == minPlane) ePlaneAve[iPlane] = 0.5*(ePlane[iPlane] + ePlane[iPlane+1]);
    else if(iPlane == maxPlane) ePlaneAve[iPlane] = 0.5*(ePlane[iPlane-1] + ePlane[iPlane]);
    else ePlaneAve[iPlane] = (ePlane[iPlane-1] + ePlane[iPlane] + ePlane[iPlane+1])/3;//the average of 3 planes
  }//end loop over track planes

  bool isShower = false;
  bool showerStart = false;
  bool showerEnd = false;

  const int nPlaneStart = 5;
  const int nPlaneEnd = 5;

  eshower = 0;

  //start a loop over track planes
  for(int iPlane=minPlane; iPlane<=maxPlane; iPlane++){
    double mipDeDx  = 0.00157;
    double mipPlane = 2*mipDeDx*geom->Plane(iPlane)->Cell(0)->HalfD()/fabs(trk->Dir().Z());
    
    if(!isShower){//if it's not already in a shower, try to find one                                                                                         
      //5 planes in a row with energy > 2*mip, call it shower start                                                                                          
      if(iPlane<=(maxPlane-nPlaneStart))
        for(int k=0; k<nPlaneStart; ++k){
          if(k==0) showerStart= (ePlane[iPlane+k]>2*mipPlane);//use eplane or eplaneAve?
          else showerStart *= (ePlane[iPlane+k]>2*mipPlane);
          if(!showerStart) break;
        }
      else showerStart = false;
      
      if(showerStart){//found the shower!                                                                                                                    
        isShower = true;//                                                                                                                                   
        startPlane = iPlane;
      }//end if shower start                                                                                                                                 
    }//end if not shower               
    else{//if in the shower region find where it ends
      for(int k=0; k<nPlaneEnd; k++)
        if((iPlane+k)<=maxPlane){
          
          if(k==0) showerEnd = (ePlaneAve[iPlane+k]<1.5*mipPlane && ePlaneAve[iPlane+k]>0.5*mipPlane);
          else showerEnd *= (ePlaneAve[iPlane+k]<1.5*mipPlane && ePlaneAve[iPlane+k]>0.5*mipPlane);
          if(!showerEnd) break;
        }
      if(showerEnd){
        isShower = false;
        endPlane = iPlane;
      }
      
    }//end else
    
  }//end loop over track planes

  //ok here is another loop over track planes...
  if(startPlane>0 && endPlane>0 && startPlane>fPlaneMin && endPlane<fPlaneMax){
    for(int iPlane = startPlane; iPlane <= endPlane; iPlane++) {
      double mipDeDx  = 0.00157;
      double mipPlane = 2*mipDeDx*geom->Plane(iPlane)->Cell(0)->HalfD()/fabs(trk->Dir().Z());
      eshower += ePlane[iPlane] - mipPlane;
    }
  
    //nPlane = endPlane - startPlane;

    fHEShower->Fill(eshower);    

    //if(eshower>0.5) fTree->Fill();
  }
 

  if(eshower < fEnergyCut) startPlane = -1;
 
  /*int startPlane;
  if(trk->Dir().Z()>0) startPlane = minPlane;//going right
  else startPlane = maxPlane;//going left

  int iPlane = startPlane;
  while(iPlane!=endPlane){
    

    if(trk->Dir().Z()>0) iPlane++;
    else iPlane--;
  }
  */
  return;  
}//end findShowerByReco

/*
std::vector<rb::WeightedHit> getTrackHits(const rb::Track &trk, const rb::Cluster &slice){
  //prepare an cluster object to save the shower hits     
  //rb::Cluster shower;
  
  //bool isContained = true;
  //bool isShower = false;
  bool hasShower = false;
  
  //int recoDecayPlane;
  //recoDecayPlane = decayPlane;
  
  //loop over slice cell hits
  //if(j>=slices.size()){
  //continue;
  //}

  //get the direction of this track
  TVector3 direction = trk->Dir();

  //double length = trk->TotalLength();                                                                                                                      

  int minPlane = trk->MinPlane();
  int maxPlane = trk->MaxPlane();
  int nPlane = geom->NPlanes();

  if(direction.Z()<0.75 ||
     fabs(maxPlane-minPlane)<0.3*nPlane ||
     hasShower = false)//if this is simply not the track we want, remove all the hits in the slice

  for(unsigned int iCell=0; iCell<slice->NCell(); ++iCell){
    art::Ptr<rb::CellHit> chit = slices[j]->Cell(iCell);
    const rb::RecoHit rhit = slices[j]->RecoHit(slices[j]->Cell(iCell));
    
    //edepSlice[chit->Plane()][chit->Cell()] += rhit.GeV();
    //need to know where is this cell hit
    double xyz[3] = {0.};           
        geom->Plane(chit->Plane())->Cell(chit->Cell())->GetCenter(xyz);     
        
        //use truth for now
        //isShower = (particle->Momentum().Pz()>0)==(xyz[2]>particle->Position().Z());
        isShower = (trk->Dir().Z()>0)==(chit->Plane()>recoDecayPlane);
        if(!isShower) continue;
        hasShower = true;           
        
        isContained*=geom->isInsideFiducialVolume(xyz[0], xyz[1], xyz[2]);
        
        if(!isContained){
          break;            
        }           
        edepShower[chit->Plane()][chit->Cell()] += rhit.GeV();
        shower.Add(chit);             
        
      }//end loop over cell hits
      
      if(!hasShower) continue;
      if(!isContained) continue;
      
      //loop over track cell hits
      for(unsigned int iCell=0; iCell<trk->NCell(); ++iCell){
        art::Ptr<rb::CellHit> chit = trk->Cell(iCell);
        const rb::RecoHit rhit = trk->RecoHit(trk->Cell(iCell));
        
            edepTrack[chit->Plane()][chit->Cell()] += rhit.GeV();
      }//end loop over cell hits 
      
      //use true vertex for now
      TVector3 eStart(particle->Position().X(), particle->Position().Y(), particle->Position().Z());      
      
      rb::Track etrack(shower, eStart, trk->Dir());
      etrack.AppendTrajectoryPoint(trk->Stop());
      
      showercol->push_back(shower);
      etrackcol->push_back(etrack);
      dokeep = true;
      fTree->Fill();
      break;
    }//end loop over tracks
}
*/

//find the DiF shower region according to the truth
void bsf::BremShowerFilter::findShowerByTruth(art::Ptr<rb::Track> &trk, int &startPlane, int &endPlane){

  art::ServiceHandle<cheat::BackTracker> bt;
  const sim::ParticleNavigator &pnav = bt->ParticleNavigator();
  art::ServiceHandle<geo::Geometry> geom;

  if(!(bt->HaveTruthInfo())){
    mf::LogWarning("BremShowerFilter") <<"No truth info found! ";
    return;
  }

  //int nPlane = geom->NPlanes();
  
  if(fabs(trk->Dir().Z())<fAngleCut) return;  
  if(fabs(trk->MaxPlane()-trk->MinPlane())<fPlaneCut) return;

  //loop over true particles to find DiF
  for (sim::ParticleNavigator::const_iterator i = pnav.begin(); i != pnav.end(); ++i) {
    const sim::Particle *particle = (*i).second;
    
    //find what is the ture particle of the reco track by timing
    //if(particle->Position().T()>trk->MinTNS() && 
    // particle->Position().T()<trk->MaxTNS() ){//found it!
    //}
    //else continue;//don't look at the rest!
    
    if(particle->Position().T()<trk->MinTNS() || 
       particle->Position().T()>trk->MaxTNS() )
      continue;
    
    //----------------------------------------------------------
    //is this a DiF
    //bool isDiF = false;
    
    // electrons only
    if ( abs(particle->PdgCode()) != 11 )
      continue;

    // must from decay
    if ( particle->Process() != "Decay" )
      continue;      
    
    // make sure electrons start inside the detector
    if (!geom->isInsideFiducialVolume(particle->Position().X(), 
                                      particle->Position().Y(), 
                                      particle->Position().Z()))
      continue;      
        
    // must decay from muons      
    const sim::Particle *mother = pnav[particle->Mother()];
    if(abs(mother->PdgCode()) != 13) continue;
    
    //particleId = particle->Mother();//this is suppose to be the id of muon
    
    
    // not Michel electrons
    if ( particle->Momentum().E() < 0.06) continue;

    //get flshits of el and mu
    if (bt->ParticleToFLSHit(particle->TrackId()).size()==0) continue;
    //if (bt->ParticleToFLSHit(particle->Mother()).size()==0) continue;
    //we found DiF!
    //isDiF = true;

    if (particle->Momentum().E() < fEnergyCut) continue;


    /*p4Mu->SetPxPyPzE(mother->Momentum().Px(), 
                     mother->Momentum().Py(), 
                     mother->Momentum().Pz(), 
                     mother->Momentum().E());  

    v4->SetXYZT(particle->Position().X(), 
                particle->Position().Y(), 
                particle->Position().Z(), 
                particle->Position().T());
    
    p4El->SetPxPyPzE(particle->Momentum().Px(), 
                     particle->Momentum().Py(), 
                     particle->Momentum().Pz(), 
                     particle->Momentum().E());      
    */

    const std::vector<sim::FLSHit>& flshitsEl
      = bt->ParticleToFLSHit(particle->TrackId());
    
    int minPlane = 1000;
    int maxPlane = 0;
    
    for(unsigned int h = 1; h < flshitsEl.size(); ++h){
      int plane=flshitsEl[h].GetPlaneID();
      if(plane<minPlane) minPlane = plane;
      if(plane>maxPlane) maxPlane = plane;      
    }

    startPlane = minPlane;
    endPlane = maxPlane;        

  }//end the loop over true particles

  return;
}//end findShowerByTruth

namespace bsf{
  DEFINE_ART_MODULE(BremShowerFilter)
}