SeedMultiplicityAnalyzer.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    TrackRecoMonitoring
// Class:      SeedMultiplicityAnalyzer
// 
//
// Original Author:  Andrea Venturi
//         Created:  Mon Oct 27 17:37:53 CET 2008
// $Id: SeedMultiplicityAnalyzer.cc,v 1.8 2012/04/21 14:03:26 venturia Exp $
//
//


// system include files
#include <memory>

// user include files

#include <vector>
#include <string>
#include <numeric>

#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EDAnalyzer.h"
#include "FWCore/Framework/interface/ConsumesCollector.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/Run.h"
#include "FWCore/Framework/interface/MakerMacros.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "FWCore/ServiceRegistry/interface/Service.h"
#include "CommonTools/UtilAlgos/interface/TFileService.h"
#include "TH1F.h"
#include "TH2F.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "FWCore/Utilities/interface/InputTag.h"

#include "DataFormats/DetId/interface/DetId.h"
#include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
#include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"

#include "MagneticField/Engine/interface/MagneticField.h"
#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"

#include "FWCore/Framework/interface/ESHandle.h"
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" 
#include "TrackingTools/Records/interface/TransientRecHitRecord.h" 

#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
#include "TrackingTools/PatternTools/interface/TSCBLBuilderNoMaterial.h"



//
// class decleration
//

class SeedMultiplicityAnalyzer : public edm::EDAnalyzer {
public:
  explicit SeedMultiplicityAnalyzer(const edm::ParameterSet&);
  ~SeedMultiplicityAnalyzer() override;
  
  class FromTrackRefSeedFilter {
  public:
    FromTrackRefSeedFilter();
    FromTrackRefSeedFilter(edm::ConsumesCollector&& iC, const edm::ParameterSet& iConfig);
    const std::string& suffix() const;
    void prepareEvent(const edm::Event& iEvent);
    bool isSelected(const unsigned int iseed) const;

  private:

    std::string m_suffix;
    bool m_passthrough;
    edm::EDGetTokenT<reco::TrackCollection> m_trackcollToken;
    edm::EDGetTokenT<reco::TrackRefVector> m_seltrackrefcollToken;
    edm::Handle<reco::TrackCollection> m_tracks;
    edm::Handle<reco::TrackRefVector> m_seltrackrefs;
    
  };
  
private:
  void analyze(const edm::Event&, const edm::EventSetup&) override;

      // ----------member data ---------------------------

  std::string _buildername;
  std::vector<edm::EDGetTokenT<TrajectorySeedCollection> > _seedcollTokens;
  std::vector<unsigned int> _seedbins;
  std::vector<double> _seedmax;
  std::vector<FromTrackRefSeedFilter> _seedfilters;
  std::vector<edm::EDGetTokenT<std::map<unsigned int, int> > > _multiplicityMapTokens;
  std::vector<std::string> _labels;
  std::vector<unsigned int> _selections;
  std::vector<unsigned int> _binsmult;
  std::vector<unsigned int> _binseta;
  std::vector<double> _maxs;
  std::vector<TH1F*> _hseedmult;
  std::vector<TH1F*> _hseedeta;
  std::vector<TH2F*> _hseedphieta;
  std::vector<TH1F*> _hpixelrhmult;
  std::vector<TH2F*> _hbpixclusleneta;
  std::vector<TH2F*> _hfpixclusleneta;
  std::vector<TH2F*> _hbpixcluslenangle;
  std::vector<TH2F*> _hfpixcluslenangle;
  std::vector<std::vector<TH2F*> > _hseedmult2D;
  std::vector<std::vector<TH2F*> > _hseedeta2D;


};

//
// constants, enums and typedefs
//

//
// static data member definitions
//

//
// constructors and destructor
//
SeedMultiplicityAnalyzer::SeedMultiplicityAnalyzer(const edm::ParameterSet& iConfig):
  _buildername(iConfig.getParameter<std::string>("TTRHBuilder")),
  _seedcollTokens(),_seedbins(),_seedmax(),_seedfilters(),
  _multiplicityMapTokens(),_labels(),_selections(),_binsmult(),_binseta(),_maxs(),_hseedmult2D(),_hseedeta2D()
{
  //now do what ever initialization is needed
  
  //
  
  
  std::vector<edm::ParameterSet> seedCollectionConfigs = 
    iConfig.getParameter<std::vector<edm::ParameterSet> >("seedCollections");
  
  for(std::vector<edm::ParameterSet>::const_iterator scps=seedCollectionConfigs.begin();scps!=seedCollectionConfigs.end();++scps) {
    
    _seedcollTokens.push_back(consumes<TrajectorySeedCollection>(scps->getParameter<edm::InputTag>("src")));
    _seedbins.push_back(scps->getUntrackedParameter<unsigned int>("nBins",1000));
    _seedmax.push_back(scps->getUntrackedParameter<double>("maxValue",100000.));

    if(scps->exists("trackFilter")) {
      _seedfilters.push_back(FromTrackRefSeedFilter(consumesCollector(),scps->getParameter<edm::ParameterSet>("trackFilter")));
    }
    else {
      _seedfilters.push_back(FromTrackRefSeedFilter());
    }
  }
  
  std::vector<edm::ParameterSet> correlationConfigs = 
    iConfig.getParameter<std::vector<edm::ParameterSet> >("multiplicityCorrelations");
  
  for(std::vector<edm::ParameterSet>::const_iterator ps=correlationConfigs.begin();ps!=correlationConfigs.end();++ps) {
    
    _multiplicityMapTokens.push_back(consumes<std::map<unsigned int, int> >(ps->getParameter<edm::InputTag>("multiplicityMap")));
    _labels.push_back(ps->getParameter<std::string>("detLabel"));
    _selections.push_back(ps->getParameter<unsigned int>("detSelection"));
    _binsmult.push_back(ps->getParameter<unsigned int>("nBins"));
    _binseta.push_back(ps->getParameter<unsigned int>("nBinsEta"));
    _maxs.push_back(ps->getParameter<double>("maxValue"));
    
  }
  
  
  
  edm::Service<TFileService> tfserv;
  
  std::vector<unsigned int>::const_iterator nseedbins = _seedbins.begin();
  std::vector<double>::const_iterator seedmax = _seedmax.begin();
  std::vector<FromTrackRefSeedFilter>::const_iterator filter = _seedfilters.begin();

  for(std::vector<edm::ParameterSet>::const_iterator scps=seedCollectionConfigs.begin();scps!=seedCollectionConfigs.end();++scps,++nseedbins,++seedmax,++filter) {
    
    std::string extendedlabel = std::string(scps->getParameter<edm::InputTag>("src").encode()) + filter->suffix();

    std::string hname = extendedlabel + std::string("_mult");
    std::string htitle = extendedlabel + std::string(" seed multiplicity");
    _hseedmult.push_back(tfserv->make<TH1F>(hname.c_str(),htitle.c_str(),*nseedbins+1,0.5-*seedmax/(*nseedbins),*seedmax+0.5));
    _hseedmult[_hseedmult.size()-1]->GetXaxis()->SetTitle("seeds");  _hseedmult[_hseedmult.size()-1]->GetYaxis()->SetTitle("events"); 

    hname = extendedlabel + std::string("_eta");
    htitle = extendedlabel + std::string(" seed pseudorapidity");
    _hseedeta.push_back(tfserv->make<TH1F>(hname.c_str(),htitle.c_str(),80,-4.,4.));
    _hseedeta[_hseedeta.size()-1]->GetXaxis()->SetTitle("#eta");  _hseedeta[_hseedeta.size()-1]->GetYaxis()->SetTitle("seeds"); 

    hname = extendedlabel + std::string("_phieta");
    htitle = extendedlabel + std::string(" seed phi vs pseudorapidity");
    _hseedphieta.push_back(tfserv->make<TH2F>(hname.c_str(),htitle.c_str(),80,-4.,4.,80,-M_PI,M_PI));
    _hseedphieta[_hseedphieta.size()-1]->GetXaxis()->SetTitle("#eta");  _hseedphieta[_hseedphieta.size()-1]->GetYaxis()->SetTitle("#phi"); 

    _hseedmult2D.push_back(std::vector<TH2F*>());
    _hseedeta2D.push_back(std::vector<TH2F*>());
    

    hname = extendedlabel + std::string("_npixelrh");
    htitle = extendedlabel + std::string(" seed SiPixelRecHit multiplicity");
    _hpixelrhmult.push_back(tfserv->make<TH1F>(hname.c_str(),htitle.c_str(),5,-.5,4.5));
    _hpixelrhmult[_hpixelrhmult.size()-1]->GetXaxis()->SetTitle("NRecHits");  _hpixelrhmult[_hpixelrhmult.size()-1]->GetYaxis()->SetTitle("seeds"); 

    hname = extendedlabel + std::string("_bpixleneta");
    htitle = extendedlabel + std::string(" seed BPIX cluster length vs pseudorapidity");
    _hbpixclusleneta.push_back(tfserv->make<TH2F>(hname.c_str(),htitle.c_str(),80,-4.,4.,40,-0.5,39.5));
    _hbpixclusleneta[_hbpixclusleneta.size()-1]->GetXaxis()->SetTitle("#eta");  _hbpixclusleneta[_hbpixclusleneta.size()-1]->GetYaxis()->SetTitle("length"); 

    hname =  extendedlabel + std::string("_fpixleneta");
    htitle = extendedlabel + std::string(" seed FPIX cluster length vs pseudorapidity");
    _hfpixclusleneta.push_back(tfserv->make<TH2F>(hname.c_str(),htitle.c_str(),80,-4.,4.,40,-0.5,39.5));
    _hfpixclusleneta[_hfpixclusleneta.size()-1]->GetXaxis()->SetTitle("#eta");  _hfpixclusleneta[_hfpixclusleneta.size()-1]->GetYaxis()->SetTitle("length"); 

    hname = extendedlabel + std::string("_bpixlenangle");
    htitle = extendedlabel + std::string(" seed BPIX cluster length vs track projection");
    _hbpixcluslenangle.push_back(tfserv->make<TH2F>(hname.c_str(),htitle.c_str(),200,-1.,1.,40,-0.5,39.5));
    _hbpixcluslenangle[_hbpixcluslenangle.size()-1]->GetXaxis()->SetTitle("projection");  _hbpixcluslenangle[_hbpixcluslenangle.size()-1]->GetYaxis()->SetTitle("length"); 

    hname = extendedlabel + std::string("_fpixlenangle");
    htitle = extendedlabel + std::string(" seed FPIX cluster length vs track projection");
    _hfpixcluslenangle.push_back(tfserv->make<TH2F>(hname.c_str(),htitle.c_str(),200,-1.,1.,40,-0.5,39.5));
    _hfpixcluslenangle[_hfpixcluslenangle.size()-1]->GetXaxis()->SetTitle("projection");  _hfpixcluslenangle[_hfpixcluslenangle.size()-1]->GetYaxis()->SetTitle("length"); 


    for(unsigned int i=0;i< _multiplicityMapTokens.size() ; ++i) {
      
      
      std::string hname2D = extendedlabel + _labels[i];
      hname2D += "_mult";
      std::string htitle2D = extendedlabel + " seeds multiplicity"; 
      htitle2D += " vs ";
      htitle2D +=  _labels[i];
      htitle2D += " hits";
      _hseedmult2D[_hseedmult2D.size()-1].push_back(tfserv->make<TH2F>(hname2D.c_str(),htitle2D.c_str(),_binsmult[i],0.,_maxs[i],
                                 *nseedbins+1,0.5-*seedmax/(*nseedbins),*seedmax+0.5));
      _hseedmult2D[_hseedmult2D.size()-1][_hseedmult2D[_hseedmult2D.size()-1].size()-1]->GetXaxis()->SetTitle("hits");  
      _hseedmult2D[_hseedmult2D.size()-1][_hseedmult2D[_hseedmult2D.size()-1].size()-1]->GetYaxis()->SetTitle("seeds");  

      hname2D = extendedlabel + _labels[i];
      hname2D += "_eta";
      htitle2D = extendedlabel + " seeds pseudorapidity"; 
      htitle2D += " vs ";
      htitle2D +=  _labels[i];
      htitle2D += " hits";
      _hseedeta2D[_hseedeta2D.size()-1].push_back(tfserv->make<TH2F>(hname2D.c_str(),htitle2D.c_str(),_binseta[i],0.,_maxs[i],
                                     80,-4.,4.));
      _hseedeta2D[_hseedeta2D.size()-1][_hseedeta2D[_hseedeta2D.size()-1].size()-1]->GetXaxis()->SetTitle("hits");  
      _hseedeta2D[_hseedeta2D.size()-1][_hseedeta2D[_hseedeta2D.size()-1].size()-1]->GetYaxis()->SetTitle("#eta");  


    }

  }

}


SeedMultiplicityAnalyzer::~SeedMultiplicityAnalyzer()
{
 
   // do anything here that needs to be done at desctruction time
   // (e.g. close files, deallocate resources etc.) 

}


//
// member functions
//

// ------------ method called to for each event  ------------
void
SeedMultiplicityAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{
  using namespace edm;

  // compute cluster multiplicities
  
  std::vector<int> tmpmult(_multiplicityMapTokens.size(),-1);
  for(unsigned int i=0;i<_multiplicityMapTokens.size();++i) {
    Handle<std::map<unsigned int, int> > mults;
    iEvent.getByToken(_multiplicityMapTokens[i],mults);
    
    // check if the selection exists
    
    std::map<unsigned int, int>::const_iterator mult = mults->find(_selections[i]);
    
    if(mult!=mults->end()) {
      tmpmult[i]=mult->second;
    }
    else {
      edm::LogWarning("DetSelectionNotFound") << " DetSelection " 
                          << _selections[i] << " not found"; 
    }
  }
  
  // preparation for loop on seeds
  
  //  TrajectoryStateTransform tsTransform;
  TSCBLBuilderNoMaterial tscblBuilder; // I could have used TSCBLBuilderWithPropagator
  
  edm::ESHandle<MagneticField> theMF;
  iSetup.get<IdealMagneticFieldRecord>().get(theMF);  
  
  edm::ESHandle<TransientTrackingRecHitBuilder> theTTRHBuilder;
  iSetup.get<TransientRecHitRecord>().get(_buildername,theTTRHBuilder);
  
  // I need: 
  // - beamspot bs POSTPONED
  
  std::vector<TH1F*>::iterator histomult=_hseedmult.begin();
  std::vector<std::vector<TH2F*> >::iterator histomult2D=_hseedmult2D.begin();
  std::vector<TH1F*>::iterator histoeta=_hseedeta.begin();
  std::vector<TH2F*>::iterator histophieta=_hseedphieta.begin();
  std::vector<std::vector<TH2F*> >::iterator histoeta2D=_hseedeta2D.begin();
  std::vector<TH1F*>::iterator hpixelrhmult=_hpixelrhmult.begin();
  std::vector<TH2F*>::iterator histobpixleneta=_hbpixclusleneta.begin();
  std::vector<TH2F*>::iterator histofpixleneta=_hfpixclusleneta.begin();
  std::vector<TH2F*>::iterator histobpixlenangle=_hbpixcluslenangle.begin();
  std::vector<TH2F*>::iterator histofpixlenangle=_hfpixcluslenangle.begin();
  std::vector<FromTrackRefSeedFilter>::iterator filter=_seedfilters.begin();
  
  // loop on seed collections
  
  for(std::vector<edm::EDGetTokenT<TrajectorySeedCollection>>::const_iterator coll=_seedcollTokens.begin();
      coll!=_seedcollTokens.end() && 
    histomult!=_hseedmult.end() && histomult2D!=_hseedmult2D.end() && 
    histoeta!=_hseedeta.end() && histoeta2D!=_hseedeta2D.end() &&
    histophieta!=_hseedphieta.end() &&
    hpixelrhmult!= _hpixelrhmult.end() &&
    histobpixleneta != _hbpixclusleneta.end() && histofpixleneta != _hfpixclusleneta.end() &&
    histobpixlenangle != _hbpixcluslenangle.end() && histofpixlenangle != _hfpixcluslenangle.end() ;
      ++coll,++histomult,++histomult2D,++histoeta,++histophieta,++histoeta2D,++hpixelrhmult,
    ++histobpixleneta,++histofpixleneta,++histobpixlenangle,++histofpixlenangle,++filter) {
    
    filter->prepareEvent(iEvent);

    Handle<TrajectorySeedCollection> seeds;
    iEvent.getByToken(*coll,seeds);
    
    /*
    (*histomult)->Fill(seeds->size());
    
    for(unsigned int i=0;i<_multiplicityMaps.size();++i) {
      if(tmpmult[i]>=0) (*histomult2D)[i]->Fill(tmpmult[i],seeds->size());
    }
    */    
    
    // loop on seeds
    
    unsigned int nseeds=0;
    unsigned int iseed=0;
    for(TrajectorySeedCollection::const_iterator seed= seeds->begin(); seed!=seeds->end(); ++seed,++iseed) {

      if(filter->isSelected(iseed)) {

    ++nseeds;

    TransientTrackingRecHit::RecHitPointer recHit = theTTRHBuilder->build(&*(seed->recHits().second-1));
    TrajectoryStateOnSurface state = trajectoryStateTransform::transientState( seed->startingState(), recHit->surface(), theMF.product());
    //      TrajectoryStateClosestToBeamLine tsAtClosestApproachSeed = tscblBuilder(*state.freeState(),bs); // here I need them BS
    
    double eta = state.globalMomentum().eta();
    double phi = state.globalMomentum().phi();
    
    (*histoeta)->Fill(eta);
    (*histophieta)->Fill(eta,phi);
    
    for(unsigned int i=0;i<_multiplicityMapTokens.size();++i) {
      if(tmpmult[i]>=0) (*histoeta2D)[i]->Fill(tmpmult[i],eta);
    }
    
    // loop on rec hits
    
    TrajectorySeed::range rechits = seed->recHits();
    
    int npixelrh = 0;
    for(TrajectorySeed::const_iterator hit = rechits.first; hit!= rechits.second; ++hit) {
      const SiPixelRecHit* sphit = dynamic_cast<const SiPixelRecHit*>(&(*hit));
      if(sphit) {
        ++npixelrh;
        // compute state on recHit surface
        TransientTrackingRecHit::RecHitPointer ttrhit = theTTRHBuilder->build(&*hit);
        TrajectoryStateOnSurface tsos = trajectoryStateTransform::transientState( seed->startingState(), ttrhit->surface(), theMF.product());       

        if(sphit->geographicalId().det()==DetId::Tracker && sphit->geographicalId().subdetId()==PixelSubdetector::PixelBarrel ) {
          (*histobpixleneta)->Fill(eta,sphit->cluster()->sizeY());
          if(tsos.isValid()) {
        //      double normdx = sin(atan2(tsos.localMomentum().x(),tsos.localMomentum().z()));
        double normdx = tsos.localMomentum().x()/sqrt(tsos.localMomentum().x()*tsos.localMomentum().x()+
                                  tsos.localMomentum().z()*tsos.localMomentum().z());
        (*histobpixlenangle)->Fill(normdx,sphit->cluster()->sizeY());
          }
        }
        else if(sphit->geographicalId().det()==DetId::Tracker && sphit->geographicalId().subdetId()==PixelSubdetector::PixelEndcap) {
          (*histofpixleneta)->Fill(eta,sphit->cluster()->sizeX());
          if(tsos.isValid()) {
        //      double normdy = sin(atan2(tsos.localMomentum().y(),tsos.localMomentum().z()));
        double normdy = tsos.localMomentum().y()/sqrt(tsos.localMomentum().y()*tsos.localMomentum().y()+
                                  tsos.localMomentum().z()*tsos.localMomentum().z());
        (*histofpixlenangle)->Fill(normdy,sphit->cluster()->sizeX());
          }
        }
        else {
          edm::LogError("InconsistentSiPixelRecHit") << "SiPixelRecHit with a non-pixel DetId " << sphit->geographicalId().rawId();
        }
      }
    }
    (*hpixelrhmult)->Fill(npixelrh);
      }
    }
    (*histomult)->Fill(nseeds);
    
    for(unsigned int i=0;i<_multiplicityMapTokens.size();++i) {
      if(tmpmult[i]>=0) (*histomult2D)[i]->Fill(tmpmult[i],nseeds);
    }

  }
}

 


SeedMultiplicityAnalyzer::FromTrackRefSeedFilter::FromTrackRefSeedFilter():
  m_suffix(""),m_passthrough(true),
  m_trackcollToken(),m_seltrackrefcollToken(),m_tracks(),m_seltrackrefs() { }

SeedMultiplicityAnalyzer::FromTrackRefSeedFilter::FromTrackRefSeedFilter(edm::ConsumesCollector&& iC, const edm::ParameterSet& iConfig):
  m_suffix(iConfig.getParameter<std::string>("suffix")),
  m_passthrough(false),
  m_trackcollToken(iC.consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("trkCollection"))),
  m_seltrackrefcollToken(iC.consumes<reco::TrackRefVector>(iConfig.getParameter<edm::InputTag>("selRefTrkCollection"))),
  m_tracks(),m_seltrackrefs() { }

const std::string&  SeedMultiplicityAnalyzer::FromTrackRefSeedFilter::suffix() const {return m_suffix; }

void SeedMultiplicityAnalyzer::FromTrackRefSeedFilter::prepareEvent(const edm::Event& iEvent) {

  if(!m_passthrough) {

    iEvent.getByToken(m_trackcollToken,m_tracks);
    iEvent.getByToken(m_seltrackrefcollToken,m_seltrackrefs);

  }

  return;
}

bool SeedMultiplicityAnalyzer::FromTrackRefSeedFilter::isSelected(const unsigned int iseed) const {

  if(m_passthrough) {
    return true;
  }
  else {
    // create a reference for the element iseed of the track collection
    const reco::TrackRef trkref(m_tracks,iseed);

    // loop on the selected trackref to check if there is the same track
    for(reco::TrackRefVector::const_iterator seltrkref=m_seltrackrefs->begin();seltrkref!=m_seltrackrefs->end();++seltrkref) {
      if(trkref==*seltrkref) return true;
    }

  }
  return false;
}



//define this as a plug-in
DEFINE_FWK_MODULE(SeedMultiplicityAnalyzer);