Calibration_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       Calibration
// Module Type: analyzer
// File:        Calibration_module.cc
//
// Generated at Wed Nov 11 14:47:12 2015 by Andrey Sheshukov using artmod
// from cetpkgsupport v1_08_07.
////////////////////////////////////////////////////////////////////////

#include <iostream>
#include "art/Framework/Core/EDAnalyzer.h"
#include "canvas/Persistency/Common/FindOneP.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 "art/Framework/Principal/SubRun.h"
#include "canvas/Utilities/InputTag.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include <TH2I.h>
#include <TProfile.h>
#include <TMath.h>
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"
#include "DDTBaseDataProducts/BaseProducts.h"
#include "DDTBaseDataProducts/HitList.h"
#include "DDTBaseDataProducts/DAQHit.h"
#include "DDTBaseDataProducts/Track3D.h"
#include "DDTBaseDataProducts/HitSet.h"
#include "OnlineCalibration/ChannelMapService.h"


// This analyzer's goal is to collect data from events into calibration histograms.
// It will collect EventsForCalibration before analizing this data and store for further usage.


namespace novaddt {
  class Calibration;
}


class novaddt::Calibration : public art::EDAnalyzer {
public:
  explicit Calibration(fhicl::ParameterSet const & p);
  // The destructor generated by the compiler is fine for classes
  // without bare pointers or other resource use.

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

  // Required functions.
  void analyze(art::Event const & e) override;
  
  // Selected optional functions.
  void beginJob() override;
  void endJob() override;

private:
    void accumulate(art::Event const& e);
    void recalibrate();
    TH2I* MakeNoiseHisto(const novaddt::calib::FDHitMap &m);
novaddt::HitSet fetchAllTrackHits(const art::Event & event);

private:
  // Declare member data here.
  unsigned _EventCounter = 0;
  unsigned _EventsToCollect = 50;  //number of events to accumulate for calibration
  std::vector<std::string> _hits_id;   //[label,instance] of hits
  std::vector<std::string> _tracks_id; //[label, instance] of tracks
  TProfile* _histAtt; //attenuation study histogram
};



double FindW(const novaddt::Track3D& t,const novaddt::DAQHit& h){
    TVector3 dir=(t.End()-t.Start());
    dir*=1./dir.Z();
    TVector3 pos=t.Start()+dir*(h.Plane().val-t.Start().Z());
/*    std::cout<<"Z: "<<t.Start().Z()<<" - "<<t.End().Z()
             <<"pos["<<h.Plane().val<<"] = "
             <<pos.X()<<", "<<pos.Y()<<", "<<pos.Z()
    <<std::endl;*/
    if(h.View().val==1)return pos.Y();
    if(h.View().val==2)return pos.X();
    return -999;
}
double FindE(const novaddt::Track3D& t,const novaddt::DAQHit& h){
    TVector3 dir=(t.End()-t.Start()).Unit();
    //geometric factor = sin(angle between dir and cell orientation)
    double Factor=1;
    if(h.View().val==1)Factor=TMath::Sqrt(1-dir.X()*dir.X());
    if(h.View().val==2)Factor=TMath::Sqrt(1-dir.Y()*dir.Y());
    return h.ADC().val*Factor;
}

novaddt::Calibration::Calibration(fhicl::ParameterSet const & p):
  EDAnalyzer(p),
  _hits_id(p.get<std::vector<std::string>>("hits_id")),
  _tracks_id(p.get<std::vector<std::string>>("tracks_id"))
{
}

void novaddt::Calibration::analyze(art::Event const & e)
{
    ++_EventCounter;
     // Hitlist from event (input)
    art::Handle<novaddt::HitList> hitsHandle;
    e.getByLabel(_hits_id[0],_hits_id[1], hitsHandle);

    //Create a set of all hits in tracks
    novaddt::HitSet allTrackHits=fetchAllTrackHits(e);
    novaddt::HitSet freeHits(*hitsHandle);
    freeHits.subtract(allTrackHits);

    //get access to the shared hit map
    art::ServiceHandle<ChannelMapService> cms;
    //increment event counter
    std::cout<<"Counter="<<cms->getCounter()<<std::endl;
    cms->addCounter();
    //obtain sharable lock, to indicate we're currently R\W data in map
    auto lock = cms->GetLock();
    auto &NoiseMap = cms->GetWMap();
    //fill the hits map
    for(const auto &hit: freeHits)
        NoiseMap[hit]+=1;
    
    //now try reading the previous map
    if(cms->getCounter()==1){
        MakeNoiseHisto(cms->GetRMap());
    }
    //do the attenuation calibration:
    //read tracks from event
    art::Handle<std::vector<novaddt::Track3D>> tracksHandle;
    e.getByLabel(_tracks_id[0],_tracks_id[1], tracksHandle);
    //get associated hits
    art::FindOneP<novaddt::HitList> fohl(tracksHandle, e, _tracks_id[0]);
    for(size_t nt = 0; nt < tracksHandle->size(); ++nt){
        auto track=tracksHandle->at(nt);
        //skip non-3D tracks
        if(track.Is3D() == false)continue;
        //skip too vertical tracks
        if(std::abs(track.End().Z()-track.Start().Z())<30)continue;
        //track direction cosine
        //double Cz = (track.End()-track.Start()).CosTheta();
        //get all hits for this track
        art::Ptr<novaddt::HitList> trackHits = fohl.at(nt);
        
        for(const auto& hit: *trackHits){
            //get the distance from hit to the center of detector
             double W=FindW(track,hit);
             double E=FindE(track,hit);
          /*   std::cout<<"hit (ADC="<<hit.ADC().val
                      <<", cell="<<hit.Cell().val
                      <<", plane="<<hit.Plane().val
                      <<",view="<<hit.View().val<<")"
                      <<" W = "<<W
                      <<" E = "<<E
                      <<std::endl;*/
             _histAtt->Fill(W,E);
        }

    }
    
}

novaddt::HitSet novaddt::Calibration::fetchAllTrackHits(const art::Event & event){
    //Create a set of all hits, that belong to tracks
    novaddt::HitSet allTrackHits;
    //Fill this set from tracks
    art::Handle<std::vector<novaddt::HitList>> trackHitLists;
    event.getByLabel(_tracks_id[0],_tracks_id[1],trackHitLists);
    for(const auto& trackHits: *trackHitLists)
        allTrackHits.add(trackHits);
    return allTrackHits;
}


TH2I* novaddt::Calibration::MakeNoiseHisto(const novaddt::calib::FDHitMap& m){
   art::ServiceHandle<art::TFileService> tfs;
   auto hist=tfs->make<TH2I>(Form("h_%d",_EventCounter),
                        Form("Noise, event #%d",_EventCounter),
                        m.Nplanes(),0,m.Nplanes(),
                        m.Ncells(),0,m.Ncells());
   uint64_t sum=0;
   for(size_t nc=0;nc<m.Ncells();++nc)
     for(size_t np=0;np<m.Nplanes();++np){
        auto val=m(novaddt::Plane(np),novaddt::Cell(nc));
        sum+=val;
        hist->SetBinContent(np,nc,val);
     }
   std::cout<<"Made hist with "<<sum<<" hits"<<std::endl;
   return hist;
}

void novaddt::Calibration::beginJob()
{
    std::cout<<"---- Calibration analyzer options -----"<<std::endl;
    std::cout<<"\tEventsToCollect = "<<_EventsToCollect <<std::endl;
    std::cout<<"\tInput hits label   = ["<<_hits_id[0]  <<" : "<<_hits_id[1]  <<"]"<<std::endl;
    std::cout<<"\tInput tracks label = ["<<_tracks_id[0]<<" : "<<_tracks_id[1]<<"]"<<std::endl;
    art::ServiceHandle<art::TFileService> tfs;
    _histAtt=tfs->make<TProfile>("hAtt","Attenuation;#cell;ADC_{corr}",
                              200,0,384);
}

void novaddt::Calibration::endJob()
{
  // Implementation of optional member function here.
}

DEFINE_ART_MODULE(novaddt::Calibration)