FastMMStudy_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       FastMMStudy
// Module Type: analyzer
// File:        FastMMStudy_module.cc
//
// Created at Nov 22 2015 by Enhao Song 
////////////////////////////////////////////////////////////////////////

#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Principal/Run.h"

#include "MCCheater/BackTracker.h"
#include "Calibrator/Calibrator.h"
#include "Simulation/Particle.h"
#include "Simulation/ParticleNavigator.h"
#include "art/Framework/Principal/SubRun.h"
#include "canvas/Utilities/InputTag.h"
#include "Utilities/func/MathUtil.h"
#include "fhiclcpp/ParameterSet.h"
#include "art/Framework/Services/Optional/TFileService.h"

#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "GeometryObjects/Geo.h"
#include "GeometryObjects/PlaneGeo.h"

#include "NovaDAQConventions/DAQConventions.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/RecoHit.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Track.h"


#include <algorithm> 
#include <utility>
#include <cmath>
#include <iostream>
#include <string>
#include <vector>
#include <numeric>
#include <map>

#include "TF1.h"
#include "TH2F.h"
#include "TMath.h"
#include "TProfile.h"
#include "TGraph.h"
#include "TPad.h"
#include "TCanvas.h"
#include "TROOT.h"
#include "TTree.h"
#include "TVector3.h"
#include "TRandom3.h"

namespace zcl {
  class FastMMStudy;
}
class zcl::FastMMStudy : public art::EDAnalyzer {
  public:
    explicit FastMMStudy(fhicl::ParameterSet const & p);
    FastMMStudy(FastMMStudy const &) = delete;
    FastMMStudy(FastMMStudy &&) = delete;
    FastMMStudy & operator = (FastMMStudy const &) = delete;
    FastMMStudy & operator = (FastMMStudy &&) = delete;
    void analyze(art::Event const & e) override;
    void beginJob() override;

  private:
    std::string _slicelabel;
    std::string _sliceinstance;
    //MC info:
    int    nxhits_MC, nyhits_MC; 
    double ti_MC, tf_MC;
    double beta_MC;

    int _xMin;
    int _xMax;
    int _yMin;
    int _yMax;
    int _zMin;
    int _zMax;
    int _deltx;
    int _delty;
    int _deltz;
    int _sigmaP;
    int _trueHits;
    int    _nSatHits;
    double _stdevcellsperplane;
    double _MaxHitsPerLength;
    double _MaxOffCenter;
    int _MinNumberSurfaceHits;
    int PX_max, PY_max, PX_min, PY_min, CX_max, CX_min, CY_max, CY_min;
    long long  TX_max, TY_max, TX_min, TY_min;
    double branch_default_value_;
    TTree *tree_;
    TTree *tree2_;
    std::map<std::string, double> t_;
    std::map<std::string, double> t2_;
    art::ServiceHandle<art::TFileService> tfs_;

    void tset(std::string const& branch_name, double const& value);
    void tset2(std::string const& branch_name, double const& value);
    void StdevCellsPerPlane(rb::Cluster const &hits, double &stdev_cells_per_plane_xview, double &stdev_cells_per_plane_yview ) const;
    void NumberOfCellsPerLength(rb::Cluster const &hits,double tracklength_of_xview, double tracklength_of_yview, double &number_of_cells_per_length_xview,double &number_of_cells_per_length_yview) const;
    float Distance(unsigned const &xcell1, unsigned const  &ycell1, unsigned const &plane1,unsigned const &xcell2, unsigned const &ycell2 , unsigned const &plane2) const;
    bool SurfAssign(rb::CellHit const & hit) const;
    int GetSurfaceId(rb::CellHit const & hit) const;
    double WeightedCenterCut(rb::Cluster const &hits,float PX_min,float PX_max,float PY_min,float PY_max,float CX_min, float CX_max, float CY_min, float CY_max, double &weighted_off_center_xx,double &weighted_off_center_xz,double &weighted_off_center_yy,double &weighted_off_center_yz) const;
    double NumberOfHitsInOverlapPlanesCut(rb::Cluster const &hits,float PX_min,float PX_max,float PY_min,float PY_max ) const;
    int NumberOfSurfaceHits(rb::Cluster const &hits) const;
    void LinFit(const std::vector<double>& x, const std::vector<double>& y, double *fitpar);
};
zcl::FastMMStudy::FastMMStudy(fhicl::ParameterSet const & p)
  :
    EDAnalyzer(p),
    _slicelabel(p.get<std::string>("slice_label")),
    _sliceinstance(p.get<std::string>("slice_instance")),
    _xMin(p.get<int>("x_min")),
    _xMax(p.get<int>("x_max")),
    _yMin(p.get<int>("y_min")),
    _yMax(p.get<int>("y_max")),
    _zMin(p.get<int>("z_min")),
    _zMax(p.get<int>("z_max")),
    _deltx(p.get<int>("x_delt")),
    _delty(p.get<int>("y_delt")),
    _deltz(p.get<int>("z_delt")),
    _sigmaP(p.get<int>("sigma_P")),
    _stdevcellsperplane(p.get<double>("stdev_cells_per_plane")),
    _MaxHitsPerLength(p.get<double>("max_hits_per_length")),
    _MaxOffCenter(p.get<double>("max_off_center")),
    _MinNumberSurfaceHits(p.get<int>("min_number_surface_hits"))
{
}


void zcl::FastMMStudy::analyze(art::Event const & e)    {
  //  std::unique_ptr<std::vector<novaddt::TriggerDecision>> 
  //    trigger_decisions(new std::vector<novaddt::TriggerDecision>);
  art::ServiceHandle<geo::Geometry> geom;
  art::ServiceHandle<calib::Calibrator> calib;
  art::Handle< std::vector<rb::Cluster> > slices;
  e.getByLabel(_slicelabel,  slices);
  //e.getByLabel(_slicelabel, _sliceinstance, slices);

  //Get MC info first:
  art::ServiceHandle<cheat::BackTracker> bt;
  const sim::ParticleNavigator& pnav = bt->ParticleNavigator();
  for (sim::ParticleNavigator::const_iterator i = pnav.begin(); i != pnav.end(); ++i) {
    nxhits_MC = 0;
    nyhits_MC = 0;
    const sim::Particle *p = (*i).second;
    if (p->PdgCode()!=42 ) continue;
    //  std::cout<<"Found monopole mctruth"<<std::endl;
    const std::vector<sim::FLSHit>& flshits = bt->ParticleToFLSHit(p->TrackId());
    unsigned ntruehits_Tot = flshits.size();
    if ( ntruehits_Tot==0 ) continue;
    //  std::cout<<"Monopole mctruth hits more than 0!"<<std::endl;
    double momentum2 = (p->P())*(p->P());
    double mass2     = (p->Mass())*(p->Mass());
    beta_MC = sqrt(momentum2 / (momentum2+mass2) );
    ti_MC = flshits[0].GetEntryT();
    tf_MC = flshits[0].GetExitT();
    size_t planeID_i = flshits[0].GetPlaneID();
    if (planeID_i%2) ++nxhits_MC;
    else ++nyhits_MC;
    for (unsigned h =1; h!= flshits.size(); ++h) {
      double ti = flshits[h].GetEntryT();
      double tf = flshits[h].GetExitT(); 
      if (ti < ti_MC) {
        ti_MC = ti;
      }
      if (tf > tf_MC) {
        tf_MC = tf;
      }
      if (flshits[h].GetPlaneID() % 2) ++nxhits_MC;
      else ++nyhits_MC;
    }
    tset2("betaMC",beta_MC);
    tset2("run_number", e.run());
    tset2("subrun_number", e.subRun());
    tset2("event_number", e.event());
    tree2_->Fill();
    //1 monopole per event, so you can break out from the loop:
    break;
  }

  ///I want to remove isolated hits from beginning
  //  std::cout<<"number of slices you are getting: "<<slices->size()<<std::endl;
  // std::cout<<"event start time: "<<event_time*15.625<<std::endl;
  // std::cout<<"TDCT0: "<<TDCT0*15.625<<std::endl;
  // std::cout<<"ti MC and tf MC: "<<ti_MC<<"  "<<tf_MC<<std::endl;
  for (unsigned i = 0; i!= slices->size(); ++i) {
    const rb::Cluster& hits = (*slices)[i];
    PX_max = 0;//P is for plane .enhao
    CX_max = 0;
    CY_max = 0;
    PY_max = 0;
    PX_min = 99999;
    CX_min = 99999;
    CY_min = 99999;
    PY_min = 99999;
    TX_max = hits.XCell(0)->TDC();
    TY_max = TX_max;
    TX_min = hits.XCell(0)->TDC();
    TY_min = TX_min;
    _trueHits = 0;
    _nSatHits   = 0;

    //int nb_hits = hits.NCell();
    //    std::cout<<"slice number = "<<i+1<<" number of hits: "<<nb_hits<<std::endl;
    int  entry_id = 0;
    bool PassThrough = false;
    bool Intrusion   = false;

    int nx = 0;
    int ny = 0;
    double sumADC=0;
    double stdevADC=0;
    double sq_sum_ADC=0;

    std::vector<double> x_hit;
    std::vector<double> y_hit;
    std::vector<double> zy_hit;
    std::vector<double> zx_hit;
    //const double peakToBaselineOffset = 96;
    for (auto const hit_ptr: hits.AllCells()) {
      const rb::CellHit hit = *hit_ptr;
      sumADC+=hit.ADC();
      sq_sum_ADC+=hit.ADC()*hit.ADC();
      if (hit.ADC()>3100) ++ _nSatHits;

      //X view hits
      if (hit.View()== geo::kX) {
        //std::cout<<"X hits"<<std::endl;
        ++nx;
        x_hit.push_back(hit.Cell());
        zx_hit.push_back(hit.Plane());
        if (PX_max < hit.Plane())   PX_max  = hit.Plane();
        if (PX_min > hit.Plane())   PX_min  = hit.Plane();
        if (CX_max < hit.Cell())    CX_max  = hit.Cell();
        if (CX_min > hit.Cell())    CX_min  = hit.Cell();
        if (TX_max < hit.TDC())     TX_max  = hit.TDC();
        if (TX_min > hit.TDC())     TX_min  = hit.TDC();  

        if (!Intrusion) {
          if (hit.Cell() < _xMin+_deltx) {
            entry_id  = 1;
            Intrusion = true;
          }
          else if (hit.Cell()> _xMax-_deltx) {
            entry_id  = 2;
            Intrusion = true;
          }
        }
        else if (!PassThrough) {
          if (hit.Cell() < _xMin+_deltx && entry_id != 1) PassThrough = true;
          else if (hit.Cell()> _xMax-_deltx && entry_id != 2) PassThrough = true;
        }
      } 
      //Y view hits
      else {
        //std::cout<<"Y hits"<<std::endl;
        ++ny;
        y_hit.push_back(hit.Cell());
        zy_hit.push_back(hit.Plane());
        if (PY_max < hit.Plane())   PY_max   = hit.Plane();
        if (PY_min > hit.Plane())   PY_min   = hit.Plane();
        if (CY_max < hit.Cell())    CY_max  = hit.Cell();
        if (CY_min > hit.Cell())    CY_min  = hit.Cell();
        if (TY_max < hit.TDC() )   TY_max   = hit.TDC();
        if (TY_min > hit.TDC() )   TY_min   = hit.TDC();

        if (!Intrusion) {
          if (hit.Cell() < _yMin+_delty) {
            entry_id  = 3;
            Intrusion = true;
          }
          else if (hit.Cell()> _yMax-_delty) {
            entry_id  = 4;
            Intrusion = true;
          }
        }
        else if (!PassThrough) {
          if (hit.Cell() < _yMin+_delty && entry_id != 3) PassThrough = true;
          else if (hit.Cell()> _yMax-_delty && entry_id != 4) PassThrough = true;
        }
      }
      //No matter which view:
      if (!Intrusion) {
        if (hit.Plane() < _zMin+_deltz) {
          entry_id  = 5;
          Intrusion = true;
        }
        else if (hit.Plane() > _zMax-_deltz) {
          entry_id  = 6;
          Intrusion = true;
        }
      }
      else if (!PassThrough) {
        if (hit.Plane() < _zMin+_deltz && entry_id != 5) PassThrough=true;
        else if (hit.Plane() > _zMax-_deltz && entry_id != 6) PassThrough=true;
      }//I think there is a bug if this hit is at corner,  it's cell value will set x intrusion true and entry id =1 and then it will pass through if it's plane val is 0 . by enhao
    }

    //Min-hits on both view:
    //       std::cout<<"xhits number : "<<nx<<"   yhits number : "<<ny<<std::endl;
    double preselection=0;
    if (_trueHits>1 || nx<3 || ny<3||!PassThrough||(std::min(PX_max,PY_max) < std::max(PX_min,PY_min)+_sigmaP)||(std::min(TX_max,TY_max) < std::max(TX_min,TY_min))) { //select cosmic I should add a handler in fcl file.
      continue;
    }
    // if (_trueHits<1 || nx<3 || ny<3||!PassThrough||(std::min(PX_max,PY_max) < std::max(PX_min,PY_min)+_sigmaP)||(std::min(TX_max,TY_max) < std::max(TX_min,TY_min))) {//select MC
    //tset("betaMC",beta_MC);
    //tset("trueHits",_trueHits);
    //tset("preselection",preselection);
    //  if(_trueHits>0) { 
    //    tree_->Fill();
    //std::cout<<"event id: "<<e.event()<<" trueHits: "<<_trueHits<<" nx: "<<nx <<" ny: "<< ny<<" Pass Through: "<<PassThrough<<" Pmax: "<<std::min(PX_max,PY_max)<<" Pmin: "<< std::max(PX_min,PY_min)<<" Tmax: "<<std::min(TX_max,TY_max) <<" Tmin:"<< std::max(TX_min,TY_min)<<std::endl;
    //  }
    //continue;
    //}

    //min number of surface hits test:
    //cut    if (NumberOfSurfaceHits(hits)<_MinNumberSurfaceHits) continue;
    float tracklength =Distance(CX_max, CY_max, PX_max, CX_min, CY_min, PX_min); 
    double tracklength_of_xview=Distance(CX_max, 0, PX_max, CX_min, 0, PX_min);
    double tracklength_of_yview=Distance(0, CY_max, PY_max, 0, CY_min, PY_min);
    double number_of_cells_per_length_xview=0;
    double number_of_cells_per_length_yview=0;
    NumberOfCellsPerLength(hits,tracklength_of_xview,tracklength_of_yview,number_of_cells_per_length_xview,number_of_cells_per_length_yview);
    double fitpar[3];
    LinFit(x_hit, zx_hit, fitpar);
    double r2x = fitpar[2];
    LinFit(y_hit, zy_hit, fitpar);
    double r2y = fitpar[2];
    //std::cout<<"hits/length is "<<hits.size()/tracklength<<std::endl;
    //cut    if (hits.size()/tracklength>_MaxHitsPerLength) continue;
    //I didn't add minhitsperlength because there might be noise hits make tracklength too long.

    //cut    if (StdevCellsPerPlane(hits)==false) continue;
    //cut    if(WeightedCenterCut(hits, PX_min, PX_max, PY_min, PY_max, CX_min,  CX_max,  CY_min,  CY_max )==false) continue;
    //    gROOT->SetBatch(kTRUE);
    stdevADC = std::sqrt( sq_sum_ADC*1.0/(nx+ny) - sumADC*1.0/(nx+ny)*sumADC*1.0/(nx+ny));

    // Fill Tree
    preselection=1;
    for (auto & branch : t_)
      branch.second = branch_default_value_;
    double weighted_off_center_xx=0, weighted_off_center_xz=0, weighted_off_center_yy=0, weighted_off_center_yz=0;
    double stdev_cells_per_plane_xview=0, stdev_cells_per_plane_yview=0;
    tset("run_number", e.run());
    tset("subrun_number", e.subRun());
    tset("event_number", e.event());
    tset("slice_number",i);
    //  tset("event_length", event_length);
    tset("number_of_surface_hits",NumberOfSurfaceHits(hits));
    tset("number_of_hits_per_length",hits.NCell()/tracklength);
    tset("number_of_cells_per_length_xview",number_of_cells_per_length_xview);
    tset("number_of_cells_per_length_yview",number_of_cells_per_length_yview);
    tset("track_length",tracklength);
    StdevCellsPerPlane(hits, stdev_cells_per_plane_xview, stdev_cells_per_plane_yview);
    tset("stdev_cells_per_plane_xview",stdev_cells_per_plane_xview);
    tset("stdev_cells_per_plane_yview",stdev_cells_per_plane_yview);
    tset("weighted_off_center",WeightedCenterCut(hits, PX_min, PX_max, PY_min, PY_max, CX_min,  CX_max,  CY_min,  CY_max, weighted_off_center_xx, weighted_off_center_xz, weighted_off_center_yy, weighted_off_center_yz ));
    tset("weighted_off_center_xx",weighted_off_center_xx);
    tset("weighted_off_center_xz",weighted_off_center_xz);
    tset("weighted_off_center_yy",weighted_off_center_yy);
    tset("weighted_off_center_yz",weighted_off_center_yz);
    tset("meanADC",sumADC*1.0/(nx+ny));
    tset("stdevADC",stdevADC);
    tset("nxhits",nx);
    tset("nyhits",ny);
    tset("percent_of_hits_in_overlap_planes",NumberOfHitsInOverlapPlanesCut(hits,PX_min,PX_max,PY_min,PY_max));
    tset("trueHits",_trueHits);
    tset("betaMC",beta_MC);
    tset("xhitsMC",nxhits_MC);
    tset("yhitsMC",nyhits_MC);
    tset("deltaTDC",std::min(TX_max,TY_max) - std::max(TX_min,TY_min));
    tset("deltaP",std::min(PX_max,PY_max) - std::max(PX_min,PY_min));
    tset("preselection",preselection);
    tset("nSatHits",_nSatHits);
    tset("r2x",r2x);
    tset("r2y",r2y);
    tset("dPX",PX_max-PX_min+1);
    tset("dPY",PY_max-PY_min+1);
    tset("dCX",CX_max-CX_min+1);
    tset("dCY",CY_max-CY_min+1);
    //   std::cout<<"event id: "<<e.event()<<" trueHits: "<<_trueHits<<" nx: "<<nx <<" ny: "<< ny<<" Pass Through: "<<PassThrough<<" Pmax: "<<std::min(PX_max,PY_max)<<" Pmin: "<< std::max(PX_min,PY_min)<<" Tmax: "<<std::min(TX_max,TY_max) <<" Tmin:"<< std::max(TX_min,TY_min)<<std::endl;
    tree_->Fill();
    //if(_trueHits>0) { 
    //  tree_->Fill();
    //}
  }
}

void zcl::FastMMStudy::beginJob()
{
  tree_ = tfs_->make<TTree>("Event", "Monopole Event Tree");
  tree2_ = tfs_->make<TTree>("MC", "Monopole MC Tree");
  std::vector<std::string> branch_names =
  { "run_number", "subrun_number", "event_number", "slice_number","number_of_hits_per_length","number_of_cells_per_length_xview","number_of_cells_per_length_yview","track_length","stdev_cells_per_plane_xview","stdev_cells_per_plane_yview", "weighted_off_center", "weighted_off_center_xx", "weighted_off_center_xz", "weighted_off_center_yy", "weighted_off_center_yz", "percent_of_hits_in_overlap_planes","number_of_surface_hits","betaMC","xhitsMC","yhitsMC","trueHits","preselection","meanADC","stdevADC","nxhits","nyhits","deltaTDC","deltaP","nSatHits","r2x","r2y","dPX","dPY","dCX","dCY"};
  for (auto const& name : branch_names)
    tree_->Branch(name.c_str(), &t_[name], (name + "/D").c_str());

  std::vector<std::string> branch_names2 =
  { "run_number", "subrun_number", "event_number", "betaMC"};
  for (auto const& name : branch_names2)
    tree2_->Branch(name.c_str(), &t2_[name], (name + "/D").c_str());
}

void zcl::FastMMStudy::tset(std::string const& branch_name, double const& value)
{
  auto it = t_.find(branch_name);
  if (it == t_.end())
  {
    std::cerr << "\n\t Branch " << branch_name << " does not exist! \n"
      << std::endl;
    assert(false);
  }

  it->second = value;
}

void zcl::FastMMStudy::tset2(std::string const& branch_name, double const& value)
{
  auto it = t2_.find(branch_name);
  if (it == t2_.end())
  {
    std::cerr << "\n\t Branch " << branch_name << " does not exist! \n"
      << std::endl;
    assert(false);
  }

  it->second = value;
}
void zcl::FastMMStudy::NumberOfCellsPerLength(rb::Cluster const &hits,double tracklength_of_xview, double tracklength_of_yview, double &number_of_cells_per_length_xview,double &number_of_cells_per_length_yview) const{
  //Maybe use hits instead of cells?
  std::vector<std::vector<int>> x_hits_plane(PX_max-PX_min+1,std::vector<int>(0));
  std::vector<std::vector<int>> y_hits_plane(PY_max-PY_min+1,std::vector<int>(0)); 
  for (auto const hit_ptr: hits.AllCells()) {
    const rb::CellHit hit = *hit_ptr;
    if (hit.View() == geo::kX) {
      //std::cout<<"X hits"<<std::endl;
      x_hits_plane[hit.Plane()-PX_min].push_back(hit.Cell()); 
    }
    if (hit.View() == geo::kY) {
      //std::cout<<"Y hits"<<std::endl;
      y_hits_plane[hit.Plane()-PY_min].push_back(hit.Cell()); 
    }
  }
  //following is to get unique number of cells hit per plane
  std::vector<int> x_cells_plane(PX_max-PX_min+1,0);
  std::vector<int> y_cells_plane(PY_max-PY_min+1,0); 
  for (unsigned i=0;i!=x_hits_plane.size();i++){
    std::sort(x_hits_plane[i].begin(),x_hits_plane[i].end());
    std::vector<int>::iterator it;
    it = std::unique (x_hits_plane[i].begin(),x_hits_plane[i].end());
    x_hits_plane[i].resize(std::distance(x_hits_plane[i].begin(),it));
    x_cells_plane[i]=x_hits_plane[i].size();
  }
  for (unsigned i=0;i!=y_hits_plane.size();i++){
    std::sort(y_hits_plane[i].begin(),y_hits_plane[i].end());
    std::vector<int>::iterator it;
    it = std::unique (y_hits_plane[i].begin(),y_hits_plane[i].end());
    y_hits_plane[i].resize(std::distance(y_hits_plane[i].begin(),it));
    y_cells_plane[i]=y_hits_plane[i].size();
  }
  //following is to get number of cells per track length in x view and y view
  int sum_of_cells=0;
  for (unsigned i=0;i!=x_cells_plane.size();i++){
    sum_of_cells+=x_cells_plane[i];
  }
  number_of_cells_per_length_xview=sum_of_cells/tracklength_of_xview;
  sum_of_cells=0;
  for (unsigned i=0;i!=y_cells_plane.size();i++){
    sum_of_cells+=y_cells_plane[i];
  }
  number_of_cells_per_length_yview=sum_of_cells/tracklength_of_yview;

}
void zcl::FastMMStudy::StdevCellsPerPlane(rb::Cluster const &hits, double &stdev_cells_per_plane_xview, double &stdev_cells_per_plane_yview ) const{
  //fluctuation of number of cells of each plane check by Enhao
  //I need to use cells not hits, because there might be many hits in same cell, and the cosmic shower has many different cells hits
  //Maybe use hits instead of cells?
  std::vector<std::vector<int>> x_hits_plane(PX_max-PX_min+1,std::vector<int>(0));
  std::vector<std::vector<int>> y_hits_plane(PY_max-PY_min+1,std::vector<int>(0)); 
  for (auto const hit_ptr: hits.AllCells()) {
    const rb::CellHit hit = *hit_ptr;
    if (hit.View() == geo::kX) {
      //std::cout<<"X hits"<<std::endl;
      x_hits_plane[hit.Plane()-PX_min].push_back(hit.Cell()); 
    }
    if (hit.View() == geo::kY) {
      //std::cout<<"Y hits"<<std::endl;
      y_hits_plane[hit.Plane()-PY_min].push_back(hit.Cell()); 
    }
  }
  double sum =0;
  double mean = 0;
  double sq_sum = 0;
  double stdev = 0;
  //following is to get unique number of cells hit per plane
  std::vector<int> x_cells_plane(PX_max-PX_min+1,0);
  std::vector<int> y_cells_plane(PY_max-PY_min+1,0); 
  for (unsigned i=0;i!=x_hits_plane.size();i++){
    std::sort(x_hits_plane[i].begin(),x_hits_plane[i].end());
    std::vector<int>::iterator it;
    it = std::unique (x_hits_plane[i].begin(),x_hits_plane[i].end());
    x_hits_plane[i].resize(std::distance(x_hits_plane[i].begin(),it));
    x_cells_plane[i]=x_hits_plane[i].size();
  }
  for (unsigned i=0;i!=y_hits_plane.size();i++){
    std::sort(y_hits_plane[i].begin(),y_hits_plane[i].end());
    std::vector<int>::iterator it;
    it = std::unique (y_hits_plane[i].begin(),y_hits_plane[i].end());
    y_hits_plane[i].resize(std::distance(y_hits_plane[i].begin(),it));
    y_cells_plane[i]=y_hits_plane[i].size();
  }
  sum = std::accumulate(x_cells_plane.begin(),x_cells_plane.end(),0.0);
  mean = sum / x_cells_plane.size();
  sq_sum = std::inner_product(x_cells_plane.begin(),x_cells_plane.end(),x_cells_plane.begin(),0.0);
  stdev = std::sqrt(sq_sum / x_cells_plane.size()-mean * mean);
  //        std::cout<<"x stdev is "<<stdev<<" "<<stdev/mean<<std::endl;
  //if (stdev/mean > _stdevcellsperplane) return false;
  stdev_cells_per_plane_xview=stdev/mean;
  sum = std::accumulate(y_cells_plane.begin(),y_cells_plane.end(),0.0);
  mean = sum / y_cells_plane.size();
  sq_sum = std::inner_product(y_cells_plane.begin(),y_cells_plane.end(),y_cells_plane.begin(),0.0);
  stdev = std::sqrt(sq_sum / y_cells_plane.size()-mean * mean);
  stdev_cells_per_plane_yview=stdev/mean;
  //        std::cout<<"y stdev is "<<stdev<<" "<<stdev/mean<<std::endl;
  //        std::cout<<"sum of x and y stdev is "<<sum_stdev<<std::endl;
  //if (stdev/mean > _stdevcellsperplane) return false;
}
float zcl::FastMMStudy::Distance( unsigned const &xcell1, unsigned const  &ycell1, unsigned const &plane1,unsigned const &xcell2, unsigned const &ycell2 , unsigned const &plane2) const {
  return std::sqrt((plane1-plane2)*(plane1-plane2)*2.82072+(xcell1-xcell2)*(xcell1-xcell2)+(ycell1-ycell2)*(ycell1-ycell2));
  //This part I should use Martin's smt::Constant::cell width, plane width
}

bool zcl::FastMMStudy::SurfAssign(rb::CellHit const &hit) const{
  //Check if this is close to surface
  if      ( hit.Cell()  > _xMax - _deltx   ) return true;
  else if ( hit.Cell()  < _xMin+_deltx     ) return true;
  else if ( hit.Plane() < _zMin+_deltz     ) return true;
  else if ( hit.Plane() > _zMax - _deltz   ) return true;

  return false;
}


int zcl::FastMMStudy::GetSurfaceId(rb::CellHit const & hit) const{
  int entry_id=0;
  if (hit.View() == geo::kX) {
    if (hit.Cell() < _xMin+_deltx) {
      entry_id  = 1;
    }
    else if (hit.Cell()> _xMax-_deltx) {
      entry_id  = 2;
    }
  } 
  //Y view hits
  else {
    if (hit.Cell() < _yMin+_delty) {
      entry_id  = 3;
    }
    else if (hit.Cell()> _yMax-_delty) {
      entry_id  = 4;
    }
  }
  //No matter which view:
  if (hit.Plane() < _zMin+_deltz) {
    entry_id  = 5;
  }
  else if (hit.Plane() > _zMax-_deltz) {
    entry_id  = 6;
  }
  return entry_id;

}

int zcl::FastMMStudy::NumberOfSurfaceHits(rb::Cluster const &hits) const{
  int number_of_surface_hits=0;
  for(unsigned i=0;i!=hits.NCell();++i){
    if(SurfAssign(*(hits.Cell(i)))==true) number_of_surface_hits++;
  }
  return number_of_surface_hits;
}
double zcl::FastMMStudy::WeightedCenterCut(rb::Cluster const &hits,float PX_min,float PX_max,float PY_min,float PY_max,float CX_min, float CX_max, float CY_min, float CY_max, double &weighted_off_center_xx,double &weighted_off_center_xz,double &weighted_off_center_yy,double &weighted_off_center_yz ) const{
  float Weighted_PX_Center=0.5,Weighted_CX_Center=0.5,Weighted_PY_Center=0.5,Weighted_CY_Center=0.5,SumOfXHitsADC=0,SumOfYHitsADC=0;
  for (auto const hit_ptr: hits.AllCells()) {
    const rb::CellHit hit = *hit_ptr;
    if (hit.View() == geo::kX) {
      SumOfXHitsADC+=hit.ADC();
      Weighted_PX_Center+=hit.ADC()*(hit.Plane()-PX_min);
      Weighted_CX_Center+=hit.ADC()*(hit.Cell()-CX_min);
    }
    if (hit.View() == geo::kY) {
      SumOfYHitsADC+=hit.ADC();
      Weighted_PY_Center+=hit.ADC()*(hit.Plane()-PY_min);
      Weighted_CY_Center+=hit.ADC()*(hit.Cell()-CY_min);
    }
  }
  if(PX_max!=PX_min)  Weighted_PX_Center=Weighted_PX_Center/SumOfXHitsADC/(PX_max-PX_min);
  if(PY_max!=PY_min)  Weighted_PY_Center=Weighted_PY_Center/SumOfYHitsADC/(PY_max-PY_min);
  if(CX_max!=CX_min)  Weighted_CX_Center=Weighted_CX_Center/SumOfXHitsADC/(CX_max-CX_min);
  if(CY_max!=CY_min)  Weighted_CY_Center=Weighted_CY_Center/SumOfYHitsADC/(CY_max-CY_min);
  //  std::cout<<"weighted px, py, cx, cy center is "<<Weighted_PX_Center<<" "<<Weighted_PY_Center<<" "<<Weighted_CX_Center<<" "<<Weighted_CY_Center<<std::endl;
  // if(std::abs(Weighted_PX_Center-0.5)>_MaxOffCenter||std::abs(Weighted_PY_Center-0.5)>_MaxOffCenter||std::abs(Weighted_CX_Center-0.5)>_MaxOffCenter||std::abs(Weighted_CY_Center-0.5)>_MaxOffCenter) return false;
  //else  return true;
  weighted_off_center_xx=std::abs(Weighted_CX_Center-0.5);
  weighted_off_center_xz=std::abs(Weighted_PX_Center-0.5);
  weighted_off_center_yy=std::abs(Weighted_CY_Center-0.5); 
  weighted_off_center_yz=std::abs(Weighted_PY_Center-0.5);
  return std::max(std::max(std::abs(Weighted_PX_Center-0.5),std::abs(Weighted_PY_Center-0.5)),std::max(std::abs(Weighted_CX_Center-0.5),std::abs(Weighted_CY_Center-0.5)));
}

double zcl::FastMMStudy::NumberOfHitsInOverlapPlanesCut(rb::Cluster const &hits,float PX_min,float PX_max,float PY_min,float PY_max ) const{
  int max=std::min(PX_max,PY_max);
  int min=std::max(PX_min,PY_min);
  double number_hits_in_overlap_planes=0;
  for (auto const hit_ptr: hits.AllCells()) {
    const rb::CellHit hit = *hit_ptr;
    if(hit.Plane()<=max && hit.Plane()>=min)    number_hits_in_overlap_planes++;
  }

  return number_hits_in_overlap_planes/hits.NCell();

}

void zcl::FastMMStudy::LinFit(const std::vector<double>& x_val, const std::vector<double>& y_val, double *fitpar){
  // http://en.wikipedia.org/wiki/Simple_linear_regression
  const auto n    = x_val.size();
  const auto x  = (std::accumulate(x_val.begin(), x_val.end(), 0.0))/n;                       // <x>
  const auto y  = (std::accumulate(y_val.begin(), y_val.end(), 0.0))/n;                       // <y>
  const auto xx = (std::inner_product(x_val.begin(), x_val.end(), x_val.begin(), 0.0))/n;     // <xx>
  const auto yy = (std::inner_product(y_val.begin(), y_val.end(), y_val.begin(), 0.0))/n;     // <yy>
  const auto xy = (std::inner_product(x_val.begin(), x_val.end(), y_val.begin(), 0.0))/n;     // <xy>

  const auto b    = (xy - x*y)/(xx - x*x);                                                    // slope
  const auto a    = y - b*x;                                                                  // intercept
  const auto r    = (xy - x*y)/sqrt((xx - x*x)*(yy - y*y));                                   // Rxy - coeffcient of determination
  fitpar[0] = a;
  fitpar[1] = b;
  fitpar[2] = r*r;

}

DEFINE_ART_MODULE(zcl::FastMMStudy)