TestOutliers.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    TestOutliers
// Class:      TestOutliers
//
//
// Original Author:  Giuseppe Cerati
//         Created:  Mon Sep 17 10:31:30 CEST 2007
//
//

// system include files
#include <memory>
#include <vector>

// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/one/EDAnalyzer.h"

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

#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/InputTag.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "SimDataFormats/Associations/interface/TrackToTrackingParticleAssociator.h"
#include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h"
#include <TH1.h>
#include <TH2.h>
#include <TFile.h>
#include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
#include "TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h"
#include "MagneticField/Engine/interface/MagneticField.h"
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
#include "Geometry/CommonDetUnit/interface/GeomDet.h"
#include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
#include "Geometry/Records/interface/TrackerTopologyRcd.h"
#include "DataFormats/GeometryCommonDetAlgo/interface/ErrorFrameTransformer.h"
#include "DataFormats/BeamSpot/interface/BeamSpot.h"

//
// class decleration
//

class TestOutliers : public edm::one::EDAnalyzer<edm::one::WatchRuns> {
public:
  explicit TestOutliers(const edm::ParameterSet &);
  ~TestOutliers() override;

private:
  void beginRun(edm::Run const &run, const edm::EventSetup &) override;
  void analyze(const edm::Event &, const edm::EventSetup &) override;
  void endRun(edm::Run const &run, const edm::EventSetup &) override {}
  void endJob() override;

  // ----------member data ---------------------------
  edm::InputTag trackTagsOut_;  //used to select what tracks to read from configuration file
  edm::InputTag trackTagsOld_;  //used to select what tracks to read from configuration file
  edm::InputTag tpTags_;        //used to select what tracks to read from configuration file
  TrackerHitAssociator::Config trackerHitAssociatorConfig_;
  edm::EDGetTokenT<reco::TrackToTrackingParticleAssociator> theAssociatorOldToken;
  edm::EDGetTokenT<reco::TrackToTrackingParticleAssociator> theAssociatorOutToken;
  const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> theGToken;
  edm::ESHandle<TrackerGeometry> theG;
  const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> theTopoToken;
  const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> theMFToken;
  std::string out;
  TFile *file;
  TH1F *histoPtOut, *histoPtOld;
  TH1F *histoQoverpOut, *histoQoverpOld;
  TH1F *histoPhiOut, *histoPhiOld;
  TH1F *histoD0Out, *histoD0Old;
  TH1F *histoDzOut, *histoDzOld;
  TH1F *histoLambdaOut, *histoLambdaOld;
  TH1F *deltahits, *deltahitsAssocGained, *deltahitsAssocLost, *hitsPerTrackLost, *hitsPerTrackAssocLost,
      *hitsPerTrackGained, *hitsPerTrackAssocGained, *deltahitsOK, *deltahitsNO;
  TH1F *okcutsOut, *okcutsOld;
  TH1F *tracks, *goodbadhits, *process, *goodcluster, *goodprocess, *badcluster, *badprocess;
  TH1F *goodhittype, *goodlayer, *goodhittype_clgteq4, *goodlayer_clgteq4, *goodhittype_cllt4, *goodlayer_cllt4,
      *badhittype, *badlayer;
  TH2F *posxy, *poszr;
  TH1F *goodpixgteq4_simvecsize, *goodpixlt4_simvecsize;
  TH1F *goodst1gteq4_simvecsize, *goodst1lt4_simvecsize;
  TH1F *goodst2gteq4_simvecsize, *goodst2lt4_simvecsize;
  TH1F *goodprjgteq4_simvecsize, *goodprjlt4_simvecsize;
  TH1F *goodpix_clustersize;
  TH1F *goodst1_clustersize;
  TH1F *goodst2_clustersize;
  TH1F *goodprj_clustersize;
  TH1F *goodpix_simvecsize;
  TH1F *goodst1_simvecsize;
  TH1F *goodst2_simvecsize;
  TH1F *goodprj_simvecsize;
  TH1F *goodhittype_simvecsmall, *goodlayer_simvecsmall, *goodhittype_simvecbig, *goodlayer_simvecbig, *goodbadmerged,
      *goodbadmergedLost, *goodbadmergedGained;
  TH1F *energyLoss, *energyLossMax, *energyLossRatio, *nOfTrackIds, *mergedPull, *mergedlayer, *mergedcluster,
      *mergedhittype;
  TH1F *sizeOut, *sizeOld, *sizeOutT, *sizeOldT;
  TH1F *countOutA, *countOutT, *countOldT;
  TH1F *gainedhits, *gainedhits2;
  TH1F *probXgood, *probXbad, *probXdelta, *probXshared, *probXnoshare;
  TH1F *probYgood, *probYbad, *probYdelta, *probYshared, *probYnoshare;
};
//
// constants, enums and typedefs
//

//
// static data member definitions
//

using namespace edm;

//
// constructors and destructor
//
TestOutliers::TestOutliers(const edm::ParameterSet &iConfig)
    : trackTagsOut_(iConfig.getUntrackedParameter<edm::InputTag>("tracksOut")),
      trackTagsOld_(iConfig.getUntrackedParameter<edm::InputTag>("tracksOld")),
      tpTags_(iConfig.getUntrackedParameter<edm::InputTag>("tp")),
      trackerHitAssociatorConfig_(consumesCollector()),
      theAssociatorOldToken(consumes<reco::TrackToTrackingParticleAssociator>(
          iConfig.getUntrackedParameter<edm::InputTag>("TrackAssociatorByHitsOld"))),
      theAssociatorOutToken(consumes<reco::TrackToTrackingParticleAssociator>(
          iConfig.getUntrackedParameter<edm::InputTag>("TrackAssociatorByHitsOut"))),
      theGToken(esConsumes<edm::Transition::BeginRun>()),
      theTopoToken(esConsumes()),
      theMFToken(esConsumes()),
      out(iConfig.getParameter<std::string>("out")) {
  LogTrace("TestOutliers") << "constructor";
  //   ParameterSet cuts = iConfig.getParameter<ParameterSet>("RecoTracksCuts");
  //   selectRecoTracks = RecoTrackSelector(cuts.getParameter<double>("ptMin"),
  //                       cuts.getParameter<double>("minRapidity"),
  //                       cuts.getParameter<double>("maxRapidity"),
  //                       cuts.getParameter<double>("tip"),
  //                       cuts.getParameter<double>("lip"),
  //                       cuts.getParameter<int>("minHit"),
  //                       cuts.getParameter<double>("maxChi2"));

  LogTrace("TestOutliers") << "end constructor";
}

TestOutliers::~TestOutliers() {
  // 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 TestOutliers::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
  //Retrieve tracker topology from geometry
  edm::ESHandle<TrackerTopology> tTopo = iSetup.getHandle(theTopoToken);

  theG = iSetup.getHandle(theGToken);
  edm::ESHandle<MagneticField> theMF = iSetup.getHandle(theMFToken);

  using namespace edm;
  using namespace std;
  using reco::TrackCollection;

  LogTrace("TestOutliers") << "analyze event #" << iEvent.id();

  edm::Handle<edm::View<reco::Track> > tracksOut;
  iEvent.getByLabel(trackTagsOut_, tracksOut);
  edm::Handle<edm::View<reco::Track> > tracksOld;
  iEvent.getByLabel(trackTagsOld_, tracksOld);
  Handle<TrackingParticleCollection> tps;
  iEvent.getByLabel(tpTags_, tps);
  edm::Handle<reco::BeamSpot> beamSpot;
  iEvent.getByLabel("offlineBeamSpot", beamSpot);

  TrackerHitAssociator hitAssociator(iEvent, trackerHitAssociatorConfig_);

  edm::Handle<reco::TrackToTrackingParticleAssociator> hAssociatorOld;
  iEvent.getByToken(theAssociatorOldToken, hAssociatorOld);
  const reco::TrackToTrackingParticleAssociator *theAssociatorOld = hAssociatorOld.product();

  edm::Handle<reco::TrackToTrackingParticleAssociator> hAssociatorOut;
  iEvent.getByToken(theAssociatorOutToken, hAssociatorOut);
  const reco::TrackToTrackingParticleAssociator *theAssociatorOut = hAssociatorOut.product();

  reco::RecoToSimCollection recSimCollOut = theAssociatorOut->associateRecoToSim(tracksOut, tps);
  reco::RecoToSimCollection recSimCollOld = theAssociatorOld->associateRecoToSim(tracksOld, tps);
  sizeOut->Fill(recSimCollOut.size());
  sizeOld->Fill(recSimCollOld.size());
  sizeOutT->Fill(tracksOut->size());
  sizeOldT->Fill(tracksOld->size());

  LogTrace("TestOutliers") << "tOld size=" << tracksOld->size() << " tOut size=" << tracksOut->size()
                           << " aOld size=" << recSimCollOld.size() << " aOut size=" << recSimCollOut.size();

#if 0
  LogTrace("TestOutliers") << "recSimCollOld.size()=" << recSimCollOld.size() ;
  for(reco::TrackCollection::size_type j=0; j<tracksOld->size(); ++j){
    reco::TrackRef trackOld(tracksOld, j);
    //if ( !selectRecoTracks( *trackOld,beamSpot.product() ) ) continue;
    LogTrace("TestOutliers") << "trackOld->pt()=" << trackOld->pt() << " trackOld->numberOfValidHits()=" << trackOld->numberOfValidHits();
    std::vector<pair<TrackingParticleRef, double> > tpOld;
    if(recSimCollOld.find(trackOld) != recSimCollOld.end()){
      tpOld = recSimCollOld[trackOld];
      if (tpOld.size()!=0) LogTrace("TestOutliers") << " associated" ;
      else LogTrace("TestOutliers") << " NOT associated" ;
    } else LogTrace("TestOutliers") << " NOT associated" ;
  }
  LogTrace("TestOutliers") << "recSimCollOut.size()=" << recSimCollOut.size() ;
  for(reco::TrackCollection::size_type j=0; j<tracksOut->size(); ++j){
    reco::TrackRef trackOut(tracksOut, j);
    //if ( !selectRecoTracks( *trackOut,beamSpot.product() ) ) continue;
    LogTrace("TestOutliers") << "trackOut->pt()=" << trackOut->pt() << " trackOut->numberOfValidHits()=" << trackOut->numberOfValidHits();
    std::vector<pair<TrackingParticleRef, double> > tpOut;
    if(recSimCollOut.find(trackOut) != recSimCollOut.end()){
      tpOut = recSimCollOut[trackOut];
      if (tpOut.size()!=0) LogTrace("TestOutliers") << " associated" ;
      else LogTrace("TestOutliers") << " NOT associated" ;
    } else LogTrace("TestOutliers") << " NOT associated" ;
  }
#endif

  LogTrace("TestOutliers") << "tracksOut->size()=" << tracksOut->size();
  LogTrace("TestOutliers") << "tracksOld->size()=" << tracksOld->size();

  std::vector<unsigned int> outused;
  for (unsigned int j = 0; j < tracksOld->size(); ++j) {
    countOldT->Fill(1);
    edm::RefToBase<reco::Track> trackOld(tracksOld, j);
    LogTrace("TestOutliers") << "now track old with id=" << j << " seed ref=" << trackOld->seedRef().get()
                             << " pt=" << trackOld->pt() << " eta=" << trackOld->eta()
                             << " chi2=" << trackOld->normalizedChi2()
                             << " tip= " << fabs(trackOld->dxy(beamSpot->position()))
                             << " lip= " << fabs(trackOld->dsz(beamSpot->position()));

    //     if (i==tracksOut->size()) {
    //       if(recSimCollOld.find(trackOld) != recSimCollOld.end()){
    //  vector<pair<TrackingParticleRef, double> > tpOld;
    //  tpOld = recSimCollOld[trackOld];
    //  if (tpOld.size()!=0) {
    //    LogTrace("TestOutliers") <<"no match: old associated and out lost! old has #hits=" << trackOld->numberOfValidHits()
    //                 << " and fraction=" << tpOld.begin()->second;
    //    if (tpOld.begin()->second>0.5) hitsPerTrackLost->Fill(trackOld->numberOfValidHits());
    //  }
    //       }
    //       continue;
    //     }

    std::vector<unsigned int> outtest;  //all the out tracks with the same seed ref
    for (unsigned int k = 0; k < tracksOut->size(); ++k) {
      edm::RefToBase<reco::Track> tmpOut = edm::RefToBase<reco::Track>(tracksOut, k);
      if (tmpOut->seedRef() == trackOld->seedRef()) {
        outtest.push_back(k);
      }
    }

    edm::RefToBase<reco::Track> trackOut;
    if (outtest.size() == 1) {  // if only one that's it
      trackOut = edm::RefToBase<reco::Track>(tracksOut, outtest[0]);
      LogTrace("TestOutliers") << "now track out with id=" << outtest[0] << " seed ref=" << trackOut->seedRef().get()
                               << " pt=" << trackOut->pt();
      outused.push_back(outtest[0]);
    } else if (outtest.size() > 1) {  //if > 1 take the one that shares all the hits with the old track
      for (unsigned int p = 0; p < outtest.size(); ++p) {
        edm::RefToBase<reco::Track> tmpOut = edm::RefToBase<reco::Track>(tracksOut, outtest[p]);
        bool allhits = true;
        for (trackingRecHit_iterator itOut = tmpOut->recHitsBegin(); itOut != tmpOut->recHitsEnd(); itOut++) {
          if ((*itOut)->isValid()) {
            bool thishit = false;
            for (trackingRecHit_iterator itOld = trackOld->recHitsBegin(); itOld != trackOld->recHitsEnd(); itOld++) {
              if ((*itOld)->isValid()) {
                const TrackingRecHit *kt = &(**itOld);
                if ((*itOut)->sharesInput(kt, TrackingRecHit::some)) {
                  thishit = true;
                  break;
                }
              }
            }
            if (!thishit)
              allhits = false;
          }
        }
        if (allhits) {
          trackOut = edm::RefToBase<reco::Track>(tracksOut, outtest[p]);
          LogTrace("TestOutliers") << "now track out with id=" << outtest[p]
                                   << " seed ref=" << trackOut->seedRef().get() << " pt=" << trackOut->pt();
          outused.push_back(outtest[p]);
        }
      }
    }

    if (outtest.empty() || trackOut.get() == nullptr) {  //no out track found for the old track
      if (recSimCollOld.find(trackOld) != recSimCollOld.end()) {
        vector<pair<TrackingParticleRef, double> > tpOld;
        tpOld = recSimCollOld[trackOld];
        if (!tpOld.empty()) {
          LogTrace("TestOutliers") << "no match: old associated and out lost! old has #hits="
                                   << trackOld->numberOfValidHits() << " and fraction=" << tpOld.begin()->second;
          if (tpOld.begin()->second > 0.5)
            hitsPerTrackLost->Fill(trackOld->numberOfValidHits());
        }
      }
      LogTrace("TestOutliers") << "...skip to next old track";
      continue;
    }

    //look if old and out are associated
    LogTrace("TestOutliers") << "trackOut->seedRef()=" << trackOut->seedRef().get()
                             << " trackOld->seedRef()=" << trackOld->seedRef().get();
    bool oldAssoc = recSimCollOld.find(trackOld) != recSimCollOld.end();
    bool outAssoc = recSimCollOut.find(trackOut) != recSimCollOut.end();
    LogTrace("TestOutliers") << "outAssoc=" << outAssoc << " oldAssoc=" << oldAssoc;

    //     if ( trackOut->seedRef()!= trackOld->seedRef() ||
    //   (trackOut->seedRef() == trackOld->seedRef() && trackOut->numberOfValidHits()>trackOld->numberOfValidHits()) ) {
    //       LogTrace("TestOutliers") <<"out and old tracks does not match...";
    //       LogTrace("TestOutliers") <<"old has #hits=" << trackOld->numberOfValidHits();
    //       std::vector<pair<TrackingParticleRef, double> > tpOld;
    //       if(recSimCollOld.find(trackOld) != recSimCollOld.end()) {
    //  tpOld = recSimCollOld[trackOld];
    //  if (tpOld.size()!=0) {
    //    LogTrace("TestOutliers") <<"old was associated with fraction=" << tpOld.begin()->second;
    //    if (tpOld.begin()->second>0.5) hitsPerTrackLost->Fill(trackOld->numberOfValidHits());
    //  }
    //       }
    //       LogTrace("TestOutliers") <<"...skip to next old track";
    //       continue;
    //     }
    //     //++i;
    countOutT->Fill(1);

    //if ( !selectRecoTracks( *trackOld,beamSpot.product() ) ) continue;//no more cuts

    tracks->Fill(0);  //FIXME

    TrackingParticleRef tpr;
    TrackingParticleRef tprOut;
    TrackingParticleRef tprOld;
    double fracOut;
    std::vector<unsigned int> tpids;
    std::vector<std::pair<TrackingParticleRef, double> > tpOut;
    std::vector<pair<TrackingParticleRef, double> > tpOld;
    //contare outliers delle tracce che erano associate e non lo sono piu!!!!

    if (outAssoc) {  //save the ids od the tp associate to the out track
      tpOut = recSimCollOut[trackOut];
      if (!tpOut.empty()) {
        countOutA->Fill(1);
        tprOut = tpOut.begin()->first;
        fracOut = tpOut.begin()->second;
        for (TrackingParticle::g4t_iterator g4T = tprOut->g4Track_begin(); g4T != tprOut->g4Track_end(); ++g4T) {
          tpids.push_back(g4T->trackId());
        }
      }
    }

    if (oldAssoc) {  //save the ids od the tp associate to the old track
      tpOld = recSimCollOld[trackOld];
      if (!tpOld.empty()) {
        tprOld = tpOld.begin()->first;
        //  LogTrace("TestOutliers") <<"old associated and out not! old has #hits=" << trackOld->numberOfValidHits()
        //               << " and fraction=" << tpOld.begin()->second;
        //  if (tpOld.begin()->second>0.5) hitsPerTrackLost->Fill(trackOld->numberOfValidHits());//deve essere un plot diverso tipo LostAssoc
        if (tpOut.empty()) {
          for (TrackingParticle::g4t_iterator g4T = tprOld->g4Track_begin(); g4T != tprOld->g4Track_end(); ++g4T) {
            tpids.push_back(g4T->trackId());
          }
        }
      }
    }

    if (tprOut.get() != nullptr || tprOld.get() != nullptr) {  //at least one of the tracks has to be associated

      tpr = tprOut.get() != nullptr ? tprOut : tprOld;

      const SimTrack *assocTrack = &(*tpr->g4Track_begin());

      //if ( trackOut->numberOfValidHits() < trackOld->numberOfValidHits() ) {
      if (trackOut->numberOfValidHits() != trackOld->numberOfValidHits() ||
          !(*trackOut->recHitsBegin())->sharesInput((*trackOld->recHitsBegin()), TrackingRecHit::some) ||
          !(*(trackOut->recHitsEnd() - 1))
               ->sharesInput(
                   (*(trackOld->recHitsEnd() - 1)),
                   TrackingRecHit::
                       some)) {  //there are outliers if the number of valid hits is != or if the first and last hit does not match
        LogTrace("TestOutliers") << "outliers for track with #hits=" << trackOut->numberOfValidHits();
        tracks->Fill(1);
        LogTrace("TestOutliers")
            << "Out->pt=" << trackOut->pt() << " Old->pt=" << trackOld->pt() << " tp->pt="
            << sqrt(tpr->momentum().perp2())
            //<< " trackOut->ptError()=" << trackOut->ptError() << " trackOld->ptError()=" << trackOld->ptError()
            << " Old->validHits=" << trackOld->numberOfValidHits() << " Out->validHits="
            << trackOut->numberOfValidHits()
            /*<< " fracOld=" << fracOld*/
            << " fracOut=" << fracOut << " deltaHits=" << trackOld->numberOfValidHits() - trackOut->numberOfValidHits();

        //compute all the track parameters
        double PtPullOut = (trackOut->pt() - sqrt(tpr->momentum().perp2())) / trackOut->ptError();
        double PtPullOld = (trackOld->pt() - sqrt(tpr->momentum().perp2())) / trackOld->ptError();
        histoPtOut->Fill(PtPullOut);
        histoPtOld->Fill(PtPullOld);

        //LogTrace("TestOutliers") << "MagneticField";
        FreeTrajectoryState ftsAtProduction(
            GlobalPoint(tpr->vertex().x(), tpr->vertex().y(), tpr->vertex().z()),
            GlobalVector(assocTrack->momentum().x(), assocTrack->momentum().y(), assocTrack->momentum().z()),
            TrackCharge(trackOld->charge()),
            theMF.product());
        TSCPBuilderNoMaterial tscpBuilder;
        TrajectoryStateClosestToPoint tsAtClosestApproach =
            tscpBuilder(ftsAtProduction, GlobalPoint(0, 0, 0));  //as in TrackProducerAlgorithm
        GlobalPoint v = tsAtClosestApproach.theState().position();
        GlobalVector p = tsAtClosestApproach.theState().momentum();

        //LogTrace("TestOutliers") << "qoverpSim";
        double qoverpSim = tsAtClosestApproach.charge() / p.mag();
        double lambdaSim = M_PI / 2 - p.theta();
        double phiSim = p.phi();
        double dxySim = (-v.x() * sin(p.phi()) + v.y() * cos(p.phi()));
        double dszSim = v.z() * p.perp() / p.mag() - (v.x() * p.x() + v.y() * p.y()) / p.perp() * p.z() / p.mag();
        double d0Sim = -dxySim;
        double dzSim = dszSim * p.mag() / p.perp();

        //LogTrace("TestOutliers") << "qoverpPullOut";
        double qoverpPullOut = (trackOut->qoverp() - qoverpSim) / trackOut->qoverpError();
        double qoverpPullOld = (trackOld->qoverp() - qoverpSim) / trackOld->qoverpError();
        double lambdaPullOut = (trackOut->lambda() - lambdaSim) / trackOut->thetaError();
        double lambdaPullOld = (trackOld->lambda() - lambdaSim) / trackOld->thetaError();
        double phi0PullOut = (trackOut->phi() - phiSim) / trackOut->phiError();
        double phi0PullOld = (trackOld->phi() - phiSim) / trackOld->phiError();
        double d0PullOut = (trackOut->d0() - d0Sim) / trackOut->d0Error();
        double d0PullOld = (trackOld->d0() - d0Sim) / trackOld->d0Error();
        double dzPullOut = (trackOut->dz() - dzSim) / trackOut->dzError();
        double dzPullOld = (trackOld->dz() - dzSim) / trackOld->dzError();

        //LogTrace("TestOutliers") << "histoQoverpOut";
        histoQoverpOut->Fill(qoverpPullOut);
        histoQoverpOld->Fill(qoverpPullOld);
        histoPhiOut->Fill(phi0PullOut);
        histoPhiOld->Fill(phi0PullOld);
        histoD0Out->Fill(d0PullOut);
        histoD0Old->Fill(d0PullOld);
        histoDzOut->Fill(dzPullOut);
        histoDzOld->Fill(dzPullOld);
        histoLambdaOut->Fill(lambdaPullOut);
        histoLambdaOld->Fill(lambdaPullOld);

        //delta number of valid hits
        LogTrace("TestOutliers") << "deltahits=" << trackOld->numberOfValidHits() - trackOut->numberOfValidHits();
        deltahits->Fill(trackOld->numberOfValidHits() - trackOut->numberOfValidHits());

        if (tprOut.get() != nullptr && tprOld.get() == nullptr) {  //out associated and old not: gained track
          if (!tpOld.empty() && tpOld.begin()->second <= 0.5) {
            deltahitsAssocGained->Fill(trackOld->numberOfValidHits() - trackOut->numberOfValidHits());
            hitsPerTrackAssocGained->Fill(trackOut->numberOfValidHits());
            LogTrace("TestOutliers") << "a) gained (assoc) track out #hits==" << trackOut->numberOfValidHits()
                                     << " old #hits=" << trackOld->numberOfValidHits();
          } else {
            deltahitsAssocGained->Fill(trackOld->numberOfValidHits() - trackOut->numberOfValidHits());
            hitsPerTrackAssocGained->Fill(trackOut->numberOfValidHits());
            LogTrace("TestOutliers") << "b) gained (assoc) track out #hits==" << trackOut->numberOfValidHits()
                                     << " old #hits=" << trackOld->numberOfValidHits();
          }
        } else if (tprOut.get() == nullptr && tprOld.get() != nullptr) {  //old associated and out not: lost track
          LogTrace("TestOutliers") << "old associated and out not! old has #hits=" << trackOld->numberOfValidHits()
                                   << " and fraction=" << tpOld.begin()->second;
          if (tpOld.begin()->second > 0.5) {
            hitsPerTrackAssocLost->Fill(trackOld->numberOfValidHits());
            deltahitsAssocLost->Fill(trackOld->numberOfValidHits() - trackOut->numberOfValidHits());
          }
        }

        if (fabs(PtPullOut) < fabs(PtPullOld))
          deltahitsOK->Fill(trackOld->numberOfValidHits() - trackOut->numberOfValidHits());
        else
          deltahitsNO->Fill(trackOld->numberOfValidHits() - trackOut->numberOfValidHits());

        //LogTrace("TestOutliers") << "RecoTrackSelector";
        //if (selectRecoTracks(*trackOut,beamSpot.product())) okcutsOut->Fill(1); else okcutsOut->Fill(0);
        //if (selectRecoTracks(*trackOld,beamSpot.product())) okcutsOld->Fill(1); else okcutsOld->Fill(0);

        LogTrace("TestOutliers") << "track old";
        for (trackingRecHit_iterator itOld = trackOld->recHitsBegin(); itOld != trackOld->recHitsEnd(); itOld++) {
          LogTrace("TestOutliers") << (*itOld)->isValid() << " " << (*itOld)->geographicalId().rawId();
        }
        LogTrace("TestOutliers") << "track out";
        for (trackingRecHit_iterator itOut = trackOut->recHitsBegin(); itOut != trackOut->recHitsEnd(); itOut++) {
          LogTrace("TestOutliers") << (*itOut)->isValid() << " " << (*itOut)->geographicalId().rawId();
        }
        //LogTrace("TestOutliers") << "itOut";

        vector<pair<int, trackingRecHit_iterator> > gainedlostoutliers;
        //look for gained hits
        for (trackingRecHit_iterator itOut = trackOut->recHitsBegin(); itOut != trackOut->recHitsEnd(); itOut++) {
          bool gained = true;
          if ((*itOut)->isValid()) {
            for (trackingRecHit_iterator itOld = trackOld->recHitsBegin(); itOld != trackOld->recHitsEnd(); itOld++) {
              if ((*itOld)->geographicalId().rawId() == (*itOut)->geographicalId().rawId())
                gained = false;
            }
            if (gained) {
              gainedlostoutliers.push_back(pair<int, trackingRecHit_iterator>(1, itOut));
              LogTrace("TestOutliers") << "broken trajectory during old fit... gained hit "
                                       << (*itOut)->geographicalId().rawId();
              gainedhits->Fill(1);
            }
          }
        }

        //look for outliers and lost hits
        for (trackingRecHit_iterator itOld = trackOld->recHitsBegin(); itOld != trackOld->recHitsEnd(); itOld++) {
          bool outlier = false;
          bool lost = true;

          for (trackingRecHit_iterator itOut = trackOut->recHitsBegin(); itOut != trackOut->recHitsEnd(); itOut++) {
            if ((*itOld)->geographicalId().rawId() == (*itOut)->geographicalId().rawId()) {
              lost = false;
              if ((*itOld)->isValid() && !(*itOut)->isValid() &&
                  (*itOld)->geographicalId().rawId() == (*itOut)->geographicalId().rawId()) {
                LogTrace("TestOutliers") << (*itOld)->isValid() << " " << (*itOut)->isValid() << " "
                                         << (*itOld)->geographicalId().rawId() << " "
                                         << (*itOut)->geographicalId().rawId();
                outlier = true;
              }
            }
          }
          if (lost)
            gainedlostoutliers.push_back(pair<int, trackingRecHit_iterator>(2, itOld));
          if (lost)
            LogTrace("TestOutliers") << "lost";
          else if (outlier)
            gainedlostoutliers.push_back(pair<int, trackingRecHit_iterator>(3, itOld));
        }

        for (std::vector<pair<int, trackingRecHit_iterator> >::iterator it = gainedlostoutliers.begin();
             it != gainedlostoutliers.end();
             ++it) {
          LogTrace("TestOutliers") << "type of processed hit:" << it->first;
          trackingRecHit_iterator itHit = it->second;
          bool gained = false;
          bool lost = false;
          bool outlier = false;
          if (it->first == 1)
            gained = true;
          else if (it->first == 2)
            lost = true;
          else if (it->first == 3)
            outlier = true;

          if (outlier || lost || gained) {
            //if (1) {

            if (lost && (*itHit)->isValid() == false) {
              goodbadmergedLost->Fill(0);
              LogTrace("TestOutliers") << "lost invalid";
              continue;
            } else if (gained && (*itHit)->isValid() == false) {
              goodbadmergedGained->Fill(0);
              LogTrace("TestOutliers") << "gained invalid";
              continue;
            }

            //LogTrace("TestOutliers") << "vector<SimHitIdpr>";
            //look if the hit comes from a correct sim track
            std::vector<SimHitIdpr> simTrackIds = hitAssociator.associateHitId(**itHit);
            bool goodhit = false;
            for (size_t j = 0; j < simTrackIds.size(); j++) {
              for (size_t jj = 0; jj < tpids.size(); jj++) {
                if (simTrackIds[j].first == tpids[jj])
                  goodhit = true;
                break;
              }
            }

            //find what kind of hit is it
            int clustersize = 0;
            int hittypeval = 0;
            int layerval = 0;
            if (dynamic_cast<const SiPixelRecHit *>(&**itHit)) {
              LogTrace("TestOutliers") << "SiPixelRecHit";
              clustersize = ((const SiPixelRecHit *)(&**itHit))->cluster()->size();
              hittypeval = 1;
            } else if (dynamic_cast<const SiStripRecHit2D *>(&**itHit)) {
              LogTrace("TestOutliers") << "SiStripRecHit2D";
              clustersize = ((const SiStripRecHit2D *)(&**itHit))->cluster()->amplitudes().size();
              hittypeval = 2;
            } else if (dynamic_cast<const SiStripMatchedRecHit2D *>(&**itHit)) {
              LogTrace("TestOutliers") << "SiStripMatchedRecHit2D";
              int clsize1 = ((const SiStripMatchedRecHit2D *)(&**itHit))->monoCluster().amplitudes().size();
              int clsize2 = ((const SiStripMatchedRecHit2D *)(&**itHit))->stereoCluster().amplitudes().size();
              if (clsize1 > clsize2)
                clustersize = clsize1;
              else
                clustersize = clsize2;
              hittypeval = 3;
            } else if (dynamic_cast<const ProjectedSiStripRecHit2D *>(&**itHit)) {
              LogTrace("TestOutliers") << "ProjectedSiStripRecHit2D";
              clustersize =
                  ((const ProjectedSiStripRecHit2D *)(&**itHit))->originalHit().cluster()->amplitudes().size();
              hittypeval = 4;
            }

            //find the layer of the hit
            int subdetId = (*itHit)->geographicalId().subdetId();
            DetId id = (*itHit)->geographicalId();
            int layerId = tTopo->layer(id);
            layerval = subdetId * 10 + layerId;

            //LogTrace("TestOutliers") << "gpos";
            GlobalPoint gpos = theG->idToDet((*itHit)->geographicalId())->surface().toGlobal((*itHit)->localPosition());

            //get the vector of sim hit associated and choose the one with the largest energy loss
            //double delta = 99999;
            //LocalPoint rhitLPv = (*itHit)->localPosition();
            //vector<PSimHit> assSimHits = hitAssociator.associateHit(**itHit);
            //if (assSimHits.size()==0) continue;
            //PSimHit shit;
            //for(std::vector<PSimHit>::const_iterator m=assSimHits.begin(); m<assSimHits.end(); m++){
            //if ((m->localPosition()-rhitLPv).mag()<delta) {
            //  shit=*m;
            //  delta = (m->localPosition()-rhitLPv).mag();
            // }
            //}
            //LogTrace("TestOutliers") << "energyLoss_";
            double energyLoss_ = 0.;
            unsigned int monoId = 0;
            std::vector<double> energyLossM;
            std::vector<double> energyLossS;
            std::vector<PSimHit> assSimHits = hitAssociator.associateHit(**itHit);
            if (assSimHits.empty())
              continue;
            PSimHit shit;
            std::vector<unsigned int> trackIds;
            energyLossS.clear();
            energyLossM.clear();
            //LogTrace("TestOutliers") << "energyLossM";
            for (std::vector<PSimHit>::const_iterator m = assSimHits.begin(); m < assSimHits.end(); m++) {
              if (outlier)
                energyLoss->Fill(m->energyLoss());
              unsigned int tId = m->trackId();
              if (find(trackIds.begin(), trackIds.end(), tId) == trackIds.end())
                trackIds.push_back(tId);
              LogTrace("TestOutliers") << "id=" << tId;
              if (m->energyLoss() > energyLoss_) {
                shit = *m;
                energyLoss_ = m->energyLoss();
              }
              if (hittypeval == 3) {
                if (monoId == 0)
                  monoId = m->detUnitId();
                if (monoId == m->detUnitId()) {
                  energyLossM.push_back(m->energyLoss());
                } else {
                  energyLossS.push_back(m->energyLoss());
                }
                //std::cout << "detUnitId="  << m->detUnitId() << " trackId=" << m->trackId() << " energyLoss=" << m->energyLoss() << std::endl;
              } else {
                energyLossM.push_back(m->energyLoss());
              }
            }
            unsigned int nIds = trackIds.size();

            if (outlier) {
              goodbadhits->Fill(goodhit);
              posxy->Fill(fabs(gpos.x()), fabs(gpos.y()));
              poszr->Fill(fabs(gpos.z()), sqrt(gpos.x() * gpos.x() + gpos.y() * gpos.y()));
              process->Fill(shit.processType());
              energyLossMax->Fill(energyLoss_);
            }

            //look if the hit is shared and if is produced only by ionization processes
            bool shared = true;
            bool ioniOnly = true;
            unsigned int idc = 0;
            for (size_t jj = 0; jj < tpids.size(); jj++) {
              idc += std::count(trackIds.begin(), trackIds.end(), tpids[jj]);
            }
            if (idc == trackIds.size()) {
              shared = false;
            }
            for (std::vector<PSimHit>::const_iterator m = assSimHits.begin() + 1; m < assSimHits.end(); m++) {
              if ((m->processType() != 7 && m->processType() != 8 && m->processType() != 9) &&
                  abs(m->particleType()) != 11) {
                ioniOnly = false;
                break;
              }
            }
            if (ioniOnly && !shared) {
              LogTrace("TestOutliers") << "delta";
            }

            if (goodhit) {
              if (outlier) {
                goodprocess->Fill(shit.processType());
                if (clustersize >= 4) {
                  goodhittype_clgteq4->Fill(hittypeval);
                  goodlayer_clgteq4->Fill(layerval);
                } else {
                  goodhittype_cllt4->Fill(hittypeval);
                  goodlayer_cllt4->Fill(layerval);
                }
                //LogTrace("TestOutliers") << "hittypeval && clustersize";
                if (hittypeval == 1 && clustersize >= 4)
                  goodpixgteq4_simvecsize->Fill(assSimHits.size());
                if (hittypeval == 1 && clustersize < 4)
                  goodpixlt4_simvecsize->Fill(assSimHits.size());
                if (hittypeval == 2 && clustersize >= 4)
                  goodst1gteq4_simvecsize->Fill(assSimHits.size());
                if (hittypeval == 2 && clustersize < 4)
                  goodst1lt4_simvecsize->Fill(assSimHits.size());
                if (hittypeval == 3 && clustersize >= 4)
                  goodst2gteq4_simvecsize->Fill(assSimHits.size());
                if (hittypeval == 3 && clustersize < 4)
                  goodst2lt4_simvecsize->Fill(assSimHits.size());
                if (hittypeval == 4 && clustersize >= 4)
                  goodprjgteq4_simvecsize->Fill(assSimHits.size());
                if (hittypeval == 4 && clustersize < 4)
                  goodprjlt4_simvecsize->Fill(assSimHits.size());

                //LogTrace("TestOutliers") << "hittypeval";
                if (hittypeval == 1)
                  goodpix_clustersize->Fill(clustersize);
                if (hittypeval == 2)
                  goodst1_clustersize->Fill(clustersize);
                if (hittypeval == 3)
                  goodst2_clustersize->Fill(clustersize);
                if (hittypeval == 4)
                  goodprj_clustersize->Fill(clustersize);
                if (hittypeval == 1)
                  goodpix_simvecsize->Fill(assSimHits.size());
                if (hittypeval == 2)
                  goodst1_simvecsize->Fill(assSimHits.size());
                if (hittypeval == 3)
                  goodst2_simvecsize->Fill(assSimHits.size());
                if (hittypeval == 4)
                  goodprj_simvecsize->Fill(assSimHits.size());

                //LogTrace("TestOutliers") << "nOfTrackIds";
                nOfTrackIds->Fill(nIds);
                if (hittypeval != 3) {
                  if (energyLossM.size() > 1) {
                    sort(energyLossM.begin(), energyLossM.end(), greater<double>());
                    energyLossRatio->Fill(energyLossM[1] / energyLossM[0]);
                  }
                } else {
                  //LogTrace("TestOutliers") << "hittypeval==3";
                  if (energyLossM.size() > 1 && energyLossS.size() <= 1) {
                    sort(energyLossM.begin(), energyLossM.end(), greater<double>());
                    energyLossRatio->Fill(energyLossM[1] / energyLossM[0]);
                  } else if (energyLossS.size() > 1 && energyLossM.size() <= 1) {
                    sort(energyLossS.begin(), energyLossS.end(), greater<double>());
                    energyLossRatio->Fill(energyLossS[1] / energyLossS[0]);
                  } else if (energyLossS.size() > 1 && energyLossM.size() > 1) {
                    sort(energyLossM.begin(), energyLossM.end(), greater<double>());
                    sort(energyLossS.begin(), energyLossS.end(), greater<double>());
                    energyLossM[1] / energyLossM[0] > energyLossS[1] / energyLossS[0]
                        ? energyLossRatio->Fill(energyLossM[1] / energyLossM[0])
                        : energyLossRatio->Fill(energyLossS[1] / energyLossS[0]);
                  }
                }

                LogTrace("TestOutliers") << "before merged";
                const SiStripMatchedRecHit2D *tmp = dynamic_cast<const SiStripMatchedRecHit2D *>(&**itHit);
                LogTrace("TestOutliers") << "tmp=" << tmp;
                LogTrace("TestOutliers") << "assSimHits.size()=" << assSimHits.size();
                if ((assSimHits.size() > 1 && tmp == nullptr) || (assSimHits.size() > 2 && tmp != nullptr)) {
                  //std::cout << "MERGED HIT" << std::endl;
                  //LogTrace("TestOutliers") << "merged";
                  mergedlayer->Fill(layerval);
                  mergedcluster->Fill(clustersize);
                  mergedhittype->Fill(hittypeval);

                  for (std::vector<PSimHit>::const_iterator m = assSimHits.begin(); m < assSimHits.end(); m++) {
                    unsigned int tId = m->trackId();
                    LogTrace("TestOutliers") << "component with id=" << tId << " eLoss=" << m->energyLoss()
                                             << " proc=" << m->processType() << " part=" << m->particleType();
                    if (find(tpids.begin(), tpids.end(), tId) == tpids.end())
                      continue;
                    if (m->processType() == 2) {
                      //GlobalPoint gpos = theG->idToDet((*itHit)->geographicalId())->surface().toGlobal((*itHit)->localPosition());
                      //GlobalPoint gpr = rhit->globalPosition();
                      AlgebraicSymMatrix33 ger = ErrorFrameTransformer()
                                                     .transform((*itHit)->localPositionError(),
                                                                theG->idToDet((*itHit)->geographicalId())->surface())
                                                     .matrix();
                      //AlgebraicSymMatrix ger = rhit->globalPositionError().matrix();
                      GlobalPoint gps =
                          theG->idToDet((*itHit)->geographicalId())->surface().toGlobal(m->localPosition());
                      //LogTrace("TestOutliers") << gpr << " " << gps << " " << ger;
                      ROOT::Math::SVector<double, 3> delta;
                      delta[0] = gpos.x() - gps.x();
                      delta[1] = gpos.y() - gps.y();
                      delta[2] = gpos.z() - gps.z();
                      LogTrace("TestOutliers") << delta << " " << ger;
                      double mpull = sqrt(delta[0] * delta[0] / ger[0][0] + delta[1] * delta[1] / ger[1][1] +
                                          delta[2] * delta[2] / ger[2][2]);
                      LogTrace("TestOutliers") << "hit pull=" << mpull;  //ger.similarity(delta);
                      mergedPull->Fill(mpull);
                    }
                  }
                  LogTrace("TestOutliers") << "end merged";
                } else {  //not merged=>good
                  //LogTrace("TestOutliers") << "goodlayer";
                  goodlayer->Fill(layerval);
                  goodcluster->Fill(clustersize);
                  goodhittype->Fill(hittypeval);
                }
              }  //if (outlier)

              const SiPixelRecHit *pix = dynamic_cast<const SiPixelRecHit *>(&**itHit);
              if ((hittypeval != 3 && assSimHits.size() < 2) || (hittypeval == 3 && assSimHits.size() < 3)) {
                if (outlier) {
                  goodhittype_simvecsmall->Fill(hittypeval);
                  goodlayer_simvecsmall->Fill(layerval);
                  goodbadmerged->Fill(1);
                  if (pix) {
                    probXgood->Fill(pix->probabilityX());
                    probYgood->Fill(pix->probabilityY());
                  }
                  LogTrace("TestOutliers") << "out good";
                } else if (lost) {
                  goodbadmergedLost->Fill(1);
                  LogTrace("TestOutliers") << "lost good";
                } else if (gained) {
                  goodbadmergedGained->Fill(1);
                  LogTrace("TestOutliers") << "gained good";
                }
                LogTrace("TestOutliers") << "good";
              } else {
                if (outlier) {
                  goodhittype_simvecbig->Fill(hittypeval);
                  goodlayer_simvecbig->Fill(layerval);
                  if (ioniOnly && !shared) {
                    goodbadmerged->Fill(3);
                    if (pix) {
                      probXdelta->Fill(pix->probabilityX());
                      probYdelta->Fill(pix->probabilityY());
                    }
                  } else if (!ioniOnly && !shared) {
                    goodbadmerged->Fill(4);
                    if (pix) {
                      probXnoshare->Fill(pix->probabilityX());
                      probYnoshare->Fill(pix->probabilityY());
                    }
                  } else {
                    goodbadmerged->Fill(5);
                    if (pix) {
                      probXshared->Fill(pix->probabilityX());
                      probYshared->Fill(pix->probabilityY());
                    }
                  }
                  LogTrace("TestOutliers") << "out merged, ioniOnly=" << ioniOnly << " shared=" << shared;
                } else if (lost) {
                  if (ioniOnly && !shared)
                    goodbadmergedLost->Fill(3);
                  else if (!ioniOnly && !shared)
                    goodbadmergedLost->Fill(4);
                  else
                    goodbadmergedLost->Fill(5);
                  LogTrace("TestOutliers") << "lost merged, ioniOnly=" << ioniOnly << " shared=" << shared;
                } else if (gained) {
                  if (ioniOnly && !shared)
                    goodbadmergedGained->Fill(3);
                  else if (!ioniOnly && !shared)
                    goodbadmergedGained->Fill(4);
                  else
                    goodbadmergedGained->Fill(5);
                  LogTrace("TestOutliers") << "gained merged, ioniOnly=" << ioniOnly << " shared=" << shared;
                }
                LogTrace("TestOutliers") << "merged, ioniOnly=" << ioniOnly << " shared=" << shared;
              }
            }       //if (goodhit)
            else {  //badhit
              //LogTrace("TestOutliers") << "badhit";
              if (outlier) {
                badcluster->Fill(clustersize);
                badhittype->Fill(hittypeval);
                badlayer->Fill(layerval);
                badprocess->Fill(shit.processType());
                goodbadmerged->Fill(2);
                const SiPixelRecHit *pix = dynamic_cast<const SiPixelRecHit *>(&**itHit);
                if (pix) {
                  probXbad->Fill(pix->probabilityX());
                  probYbad->Fill(pix->probabilityY());
                }
                LogTrace("TestOutliers") << "out bad";
              } else if (lost) {
                goodbadmergedLost->Fill(2);
                LogTrace("TestOutliers") << "lost bad";
              } else if (gained) {
                goodbadmergedGained->Fill(2);
                LogTrace("TestOutliers") << "gained bad";
              }
              LogTrace("TestOutliers") << "bad";
            }
          }
        }
      }
      //else if ( trackOut->numberOfValidHits() > trackOld->numberOfValidHits() ) {
      else if (false) {
        LogTrace("TestOutliers") << "outliers for track with #hits=" << trackOut->numberOfValidHits();
        tracks->Fill(1);
        LogTrace("TestOutliers")
            << "Out->pt=" << trackOut->pt() << " Old->pt=" << trackOld->pt() << " tp->pt="
            << sqrt(tpr->momentum().perp2())
            //<< " trackOut->ptError()=" << trackOut->ptError() << " trackOld->ptError()=" << trackOld->ptError()
            << " Old->validHits=" << trackOld->numberOfValidHits() << " Out->validHits="
            << trackOut->numberOfValidHits()
            /*<< " fracOld=" << fracOld*/
            << " fracOut=" << fracOut << " deltaHits=" << trackOld->numberOfValidHits() - trackOut->numberOfValidHits();
        LogTrace("TestOutliers") << "track with gained hits";
        gainedhits2->Fill(trackOut->numberOfValidHits() - trackOld->numberOfValidHits());
      } else {
        LogTrace("TestOutliers") << "no outliers for track with #hits=" << trackOut->numberOfValidHits();
      }
    }
    LogTrace("TestOutliers") << "end track old #" << j;
  }

  for (unsigned int k = 0; k < tracksOut->size(); ++k) {
    if (find(outused.begin(), outused.end(), k) == outused.end()) {
      edm::RefToBase<reco::Track> trackOut(tracksOut, k);
      bool outAssoc = recSimCollOut.find(trackOut) != recSimCollOut.end();
      if (outAssoc) {
        hitsPerTrackGained->Fill(trackOut->numberOfValidHits());
        LogTrace("TestOutliers") << "gained track out id=" << k << " #hits==" << trackOut->numberOfValidHits();
      }
    }
  }
}

// ------------ method called once each job just before starting event loop  ------------
void TestOutliers::beginRun(edm::Run const &run, const edm::EventSetup &es) {
  const bool oldAddDir = TH1::AddDirectoryStatus();
  TH1::AddDirectory(true);
  file = new TFile(out.c_str(), "recreate");
  histoPtOut = new TH1F("histoPtOut", "histoPtOut", 100, -10, 10);
  histoPtOld = new TH1F("histoPtOld", "histoPtOld", 100, -10, 10);
  histoQoverpOut = new TH1F("histoQoverpOut", "histoQoverpOut", 250, -25, 25);
  histoQoverpOld = new TH1F("histoQoverpOld", "histoQoverpOld", 250, -25, 25);
  histoPhiOut = new TH1F("histoPhiOut", "histoPhiOut", 250, -25, 25);
  histoPhiOld = new TH1F("histoPhiOld", "histoPhiOld", 250, -25, 25);
  histoD0Out = new TH1F("histoD0Out", "histoD0Out", 250, -25, 25);
  histoD0Old = new TH1F("histoD0Old", "histoD0Old", 250, -25, 25);
  histoDzOut = new TH1F("histoDzOut", "histoDzOut", 250, -25, 25);
  histoDzOld = new TH1F("histoDzOld", "histoDzOld", 250, -25, 25);
  histoLambdaOut = new TH1F("histoLambdaOut", "histoLambdaOut", 250, -25, 25);
  histoLambdaOld = new TH1F("histoLambdaOld", "histoLambdaOld", 250, -25, 25);
  deltahits = new TH1F("deltahits", "deltahits", 80, -40, 40);
  deltahitsOK = new TH1F("deltahitsOK", "deltahitsOK", 20, 0, 20);
  deltahitsNO = new TH1F("deltahitsNO", "deltahitsNO", 20, 0, 20);
  okcutsOut = new TH1F("okcutsOut", "okcutsOut", 2, -0.5, 1.5);
  okcutsOld = new TH1F("okcutsOld", "okcutsOld", 2, -0.5, 1.5);
  tracks = new TH1F("tracks_", "tracks_", 2, -0.5, 1.5);
  goodbadhits = new TH1F("goodbadhits", "goodbadhits", 2, -0.5, 1.5);
  process = new TH1F("process", "process", 20, -0.5, 19.5);
  goodcluster = new TH1F("goodcluster", "goodcluster", 40, -0.5, 39.5);
  goodprocess = new TH1F("goodprocess", "goodprocess", 20, -0.5, 19.5);
  badcluster = new TH1F("badcluster", "badcluster", 40, -0.5, 39.5);
  badprocess = new TH1F("badprocess", "badprocess", 20, -0.5, 19.5);
  goodhittype = new TH1F("goodhittype", "goodhittype", 5, -0.5, 4.5);
  goodlayer = new TH1F("goodlayer", "goodlayer", 70, -0.5, 69.5);
  goodhittype_clgteq4 = new TH1F("goodhittype_clgteq4", "goodhittype_clgteq4", 5, -0.5, 4.5);
  goodlayer_clgteq4 = new TH1F("goodlayer_clgteq4", "goodlayer_clgteq4", 70, -0.5, 69.5);
  goodhittype_cllt4 = new TH1F("goodhittype_cllt4", "goodhittype_cllt4", 5, -0.5, 4.5);
  goodlayer_cllt4 = new TH1F("goodlayer_cllt4", "goodlayer_cllt4", 70, -0.5, 69.5);
  badhittype = new TH1F("badhittype", "badhittype", 5, -0.5, 4.5);
  badlayer = new TH1F("badlayer", "badlayer", 70, -0.5, 69.5);
  posxy = new TH2F("posxy", "posxy", 1200, 0, 120, 1200, 0, 120);
  poszr = new TH2F("poszr", "poszr", 3000, 0, 300, 1200, 0, 120);
  goodpixgteq4_simvecsize = new TH1F("goodpixgteq4_simvecsize", "goodpixgteq4_simvecsize", 40, -0.5, 39.5);
  goodpixlt4_simvecsize = new TH1F("goodpixlt4_simvecsize", "goodpixlt4_simvecsize", 40, -0.5, 39.5);
  goodst1gteq4_simvecsize = new TH1F("goodst1gteq4_simvecsize", "goodst1gteq4_simvecsize", 40, -0.5, 39.5);
  goodst1lt4_simvecsize = new TH1F("goodst1lt4_simvecsize", "goodst1lt4_simvecsize", 40, -0.5, 39.5);
  goodst2gteq4_simvecsize = new TH1F("goodst2gteq4_simvecsize", "goodst2gteq4_simvecsize", 40, -0.5, 39.5);
  goodst2lt4_simvecsize = new TH1F("goodst2lt4_simvecsize", "goodst2lt4_simvecsize", 40, -0.5, 39.5);
  goodprjgteq4_simvecsize = new TH1F("goodprjgteq4_simvecsize", "goodprjgteq4_simvecsize", 40, -0.5, 39.5);
  goodprjlt4_simvecsize = new TH1F("goodprjlt4_simvecsize", "goodprjlt4_simvecsize", 40, -0.5, 39.5);
  goodpix_clustersize = new TH1F("goodpix_clustersize", "goodpix_clustersize", 40, -0.5, 39.5);
  goodst1_clustersize = new TH1F("goodst1_clustersize", "goodst1_clustersize", 40, -0.5, 39.5);
  goodst2_clustersize = new TH1F("goodst2_clustersize", "goodst2_clustersize", 40, -0.5, 39.5);
  goodprj_clustersize = new TH1F("goodprj_clustersize", "goodprj_clustersize", 40, -0.5, 39.5);
  goodpix_simvecsize = new TH1F("goodpix_simvecsize", "goodpix_simvecsize", 40, -0.5, 39.5);
  goodst1_simvecsize = new TH1F("goodst1_simvecsize", "goodst1_simvecsize", 40, -0.5, 39.5);
  goodst2_simvecsize = new TH1F("goodst2_simvecsize", "goodst2_simvecsize", 40, -0.5, 39.5);
  goodprj_simvecsize = new TH1F("goodprj_simvecsize", "goodprj_simvecsize", 40, -0.5, 39.5);
  goodhittype_simvecsmall = new TH1F("goodhittype_simvecsmall", "goodhittype_simvecsmall", 5, -0.5, 4.5);
  goodlayer_simvecsmall = new TH1F("goodlayer_simvecsmall", "goodlayer_simvecsmall", 70, -0.5, 69.5);
  goodhittype_simvecbig = new TH1F("goodhittype_simvecbig", "goodhittype_simvecbig", 5, -0.5, 4.5);
  goodlayer_simvecbig = new TH1F("goodlayer_simvecbig", "goodlayer_simvecbig", 70, -0.5, 69.5);
  goodbadmerged = new TH1F("goodbadmerged", "goodbadmerged", 5, 0.5, 5.5);
  goodbadmergedLost = new TH1F("goodbadmergedLost", "goodbadmergedLost", 5, 0.5, 5.5);
  goodbadmergedGained = new TH1F("goodbadmergedGained", "goodbadmergedGained", 5, 0.5, 5.5);
  energyLoss = new TH1F("energyLoss", "energyLoss", 1000, 0, 0.1);
  energyLossMax = new TH1F("energyLossMax", "energyLossMax", 1000, 0, 0.1);
  energyLossRatio = new TH1F("energyLossRatio", "energyLossRatio", 100, 0, 1);
  nOfTrackIds = new TH1F("nOfTrackIds", "nOfTrackIds", 10, 0, 10);
  mergedPull = new TH1F("mergedPull", "mergedPull", 100, 0, 10);
  mergedlayer = new TH1F("mergedlayer", "mergedlayer", 70, -0.5, 69.5);
  mergedhittype = new TH1F("mergedhittype", "mergedhittype", 5, -0.5, 4.5);
  mergedcluster = new TH1F("mergedcluster", "mergedcluster", 40, -0.5, 39.5);
  deltahitsAssocGained = new TH1F("deltahitsAssocGained", "deltahitsAssocGained", 80, -40, 40);
  deltahitsAssocLost = new TH1F("deltahitsAssocLost", "deltahitsAssocLost", 80, -40, 40);
  hitsPerTrackLost = new TH1F("hitsPerTrackLost", "hitsPerTrackLost", 40, -0.5, 39.5);
  hitsPerTrackAssocLost = new TH1F("hitsPerTrackAssocLost", "hitsPerTrackAssocLost", 40, -0.5, 39.5);
  hitsPerTrackGained = new TH1F("hitsPerTrackGained", "hitsPerTrackGained", 40, -0.5, 39.5);
  hitsPerTrackAssocGained = new TH1F("hitsPerTrackAssocGained", "hitsPerTrackAssocGained", 40, -0.5, 39.5);
  sizeOut = new TH1F("sizeOut", "sizeOut", 900, -0.5, 899.5);
  sizeOld = new TH1F("sizeOld", "sizeOld", 900, -0.5, 899.5);
  sizeOutT = new TH1F("sizeOutT", "sizeOutT", 900, -0.5, 899.5);
  sizeOldT = new TH1F("sizeOldT", "sizeOldT", 900, -0.5, 899.5);
  countOutA = new TH1F("countOutA", "countOutA", 2, 0, 2);
  countOutT = new TH1F("countOutT", "countOutT", 2, 0, 2);
  countOldT = new TH1F("countOldT", "countOldT", 2, 0, 2);
  gainedhits = new TH1F("gainedhits", "gainedhits", 2, 0, 2);
  gainedhits2 = new TH1F("gainedhits2", "gainedhits2", 30, -0.5, 29.5);
  probXgood = new TH1F("probXgood", "probXgood", 110, 0, 1.1);
  probXbad = new TH1F("probXbad", "probXbad", 110, 0, 1.1);
  probXdelta = new TH1F("probXdelta", "probXdelta", 110, 0, 1.1);
  probXshared = new TH1F("probXshared", "probXshared", 110, 0, 1.1);
  probXnoshare = new TH1F("probXnoshare", "probXnoshare", 110, 0, 1.1);
  probYgood = new TH1F("probYgood", "probYgood", 110, 0, 1.1);
  probYbad = new TH1F("probYbad", "probYbad", 110, 0, 1.1);
  probYdelta = new TH1F("probYdelta", "probYdelta", 110, 0, 1.1);
  probYshared = new TH1F("probYshared", "probYshared", 110, 0, 1.1);
  probYnoshare = new TH1F("probYnoshare", "probYnoshare", 110, 0, 1.1);
  TH1::AddDirectory(oldAddDir);
}
// ------------ method called once each job just after ending the event loop  ------------
void TestOutliers::endJob() {
  LogTrace("TestOutliers") << "TestOutliers::endJob";
  file->Write();
  LogTrace("TestOutliers") << "outfile written";
  file->Close();
  LogTrace("TestOutliers") << "oufile closed";
  LogTrace("TestOutliers") << "exiting TestOutliers::endJob";
}

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