Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-01-13 23:20:59

0001 // -*- C++ -*-
0002 //
0003 // Package:    TestOutliers
0004 // Class:      TestOutliers
0005 //
0006 /**\class TestOutliers TestOutliers.cc RecoTracker/DebugTools/plugins/TestOutliers.cc
0007 
0008    Description: <one line class summary>
0009 
0010    Implementation:
0011    <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Giuseppe Cerati
0015 //         Created:  Mon Sep 17 10:31:30 CEST 2007
0016 //
0017 //
0018 
0019 // system include files
0020 #include <memory>
0021 #include <vector>
0022 
0023 // user include files
0024 #include "FWCore/Framework/interface/Frameworkfwd.h"
0025 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0026 
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/Framework/interface/MakerMacros.h"
0029 
0030 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0031 #include "FWCore/Utilities/interface/InputTag.h"
0032 #include "DataFormats/TrackReco/interface/Track.h"
0033 #include "SimDataFormats/Associations/interface/TrackToTrackingParticleAssociator.h"
0034 #include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h"
0035 #include <TH1.h>
0036 #include <TH2.h>
0037 #include <TFile.h>
0038 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0039 #include "TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h"
0040 #include "MagneticField/Engine/interface/MagneticField.h"
0041 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0042 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0043 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0044 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0045 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0046 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0047 #include "DataFormats/GeometryCommonDetAlgo/interface/ErrorFrameTransformer.h"
0048 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0049 
0050 //
0051 // class decleration
0052 //
0053 
0054 class TestOutliers : public edm::one::EDAnalyzer<edm::one::WatchRuns> {
0055 public:
0056   explicit TestOutliers(const edm::ParameterSet &);
0057   ~TestOutliers() override;
0058 
0059 private:
0060   void beginRun(edm::Run const &run, const edm::EventSetup &) override;
0061   void analyze(const edm::Event &, const edm::EventSetup &) override;
0062   void endRun(edm::Run const &run, const edm::EventSetup &) override {}
0063   void endJob() override;
0064 
0065   // ----------member data ---------------------------
0066   edm::InputTag trackTagsOut_;  //used to select what tracks to read from configuration file
0067   edm::InputTag trackTagsOld_;  //used to select what tracks to read from configuration file
0068   edm::InputTag tpTags_;        //used to select what tracks to read from configuration file
0069   TrackerHitAssociator::Config trackerHitAssociatorConfig_;
0070   edm::EDGetTokenT<reco::TrackToTrackingParticleAssociator> theAssociatorOldToken;
0071   edm::EDGetTokenT<reco::TrackToTrackingParticleAssociator> theAssociatorOutToken;
0072   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> theGToken;
0073   edm::ESHandle<TrackerGeometry> theG;
0074   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> theTopoToken;
0075   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> theMFToken;
0076   std::string out;
0077   TFile *file;
0078   TH1F *histoPtOut, *histoPtOld;
0079   TH1F *histoQoverpOut, *histoQoverpOld;
0080   TH1F *histoPhiOut, *histoPhiOld;
0081   TH1F *histoD0Out, *histoD0Old;
0082   TH1F *histoDzOut, *histoDzOld;
0083   TH1F *histoLambdaOut, *histoLambdaOld;
0084   TH1F *deltahits, *deltahitsAssocGained, *deltahitsAssocLost, *hitsPerTrackLost, *hitsPerTrackAssocLost,
0085       *hitsPerTrackGained, *hitsPerTrackAssocGained, *deltahitsOK, *deltahitsNO;
0086   TH1F *okcutsOut, *okcutsOld;
0087   TH1F *tracks, *goodbadhits, *process, *goodcluster, *goodprocess, *badcluster, *badprocess;
0088   TH1F *goodhittype, *goodlayer, *goodhittype_clgteq4, *goodlayer_clgteq4, *goodhittype_cllt4, *goodlayer_cllt4,
0089       *badhittype, *badlayer;
0090   TH2F *posxy, *poszr;
0091   TH1F *goodpixgteq4_simvecsize, *goodpixlt4_simvecsize;
0092   TH1F *goodst1gteq4_simvecsize, *goodst1lt4_simvecsize;
0093   TH1F *goodst2gteq4_simvecsize, *goodst2lt4_simvecsize;
0094   TH1F *goodprjgteq4_simvecsize, *goodprjlt4_simvecsize;
0095   TH1F *goodpix_clustersize;
0096   TH1F *goodst1_clustersize;
0097   TH1F *goodst2_clustersize;
0098   TH1F *goodprj_clustersize;
0099   TH1F *goodpix_simvecsize;
0100   TH1F *goodst1_simvecsize;
0101   TH1F *goodst2_simvecsize;
0102   TH1F *goodprj_simvecsize;
0103   TH1F *goodhittype_simvecsmall, *goodlayer_simvecsmall, *goodhittype_simvecbig, *goodlayer_simvecbig, *goodbadmerged,
0104       *goodbadmergedLost, *goodbadmergedGained;
0105   TH1F *energyLoss, *energyLossMax, *energyLossRatio, *nOfTrackIds, *mergedPull, *mergedlayer, *mergedcluster,
0106       *mergedhittype;
0107   TH1F *sizeOut, *sizeOld, *sizeOutT, *sizeOldT;
0108   TH1F *countOutA, *countOutT, *countOldT;
0109   TH1F *gainedhits, *gainedhits2;
0110   TH1F *probXgood, *probXbad, *probXdelta, *probXshared, *probXnoshare;
0111   TH1F *probYgood, *probYbad, *probYdelta, *probYshared, *probYnoshare;
0112 };
0113 //
0114 // constants, enums and typedefs
0115 //
0116 
0117 //
0118 // static data member definitions
0119 //
0120 
0121 using namespace edm;
0122 
0123 //
0124 // constructors and destructor
0125 //
0126 TestOutliers::TestOutliers(const edm::ParameterSet &iConfig)
0127     : trackTagsOut_(iConfig.getUntrackedParameter<edm::InputTag>("tracksOut")),
0128       trackTagsOld_(iConfig.getUntrackedParameter<edm::InputTag>("tracksOld")),
0129       tpTags_(iConfig.getUntrackedParameter<edm::InputTag>("tp")),
0130       trackerHitAssociatorConfig_(consumesCollector()),
0131       theAssociatorOldToken(consumes<reco::TrackToTrackingParticleAssociator>(
0132           iConfig.getUntrackedParameter<edm::InputTag>("TrackAssociatorByHitsOld"))),
0133       theAssociatorOutToken(consumes<reco::TrackToTrackingParticleAssociator>(
0134           iConfig.getUntrackedParameter<edm::InputTag>("TrackAssociatorByHitsOut"))),
0135       theGToken(esConsumes<edm::Transition::BeginRun>()),
0136       theTopoToken(esConsumes()),
0137       theMFToken(esConsumes()),
0138       out(iConfig.getParameter<std::string>("out")) {
0139   LogTrace("TestOutliers") << "constructor";
0140   //   ParameterSet cuts = iConfig.getParameter<ParameterSet>("RecoTracksCuts");
0141   //   selectRecoTracks = RecoTrackSelector(cuts.getParameter<double>("ptMin"),
0142   //                       cuts.getParameter<double>("minRapidity"),
0143   //                       cuts.getParameter<double>("maxRapidity"),
0144   //                       cuts.getParameter<double>("tip"),
0145   //                       cuts.getParameter<double>("lip"),
0146   //                       cuts.getParameter<int>("minHit"),
0147   //                       cuts.getParameter<double>("maxChi2"));
0148 
0149   LogTrace("TestOutliers") << "end constructor";
0150 }
0151 
0152 TestOutliers::~TestOutliers() {
0153   // do anything here that needs to be done at desctruction time
0154   // (e.g. close files, deallocate resources etc.)
0155 }
0156 
0157 //
0158 // member functions
0159 //
0160 
0161 // ------------ method called to for each event  ------------
0162 void TestOutliers::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0163   //Retrieve tracker topology from geometry
0164   edm::ESHandle<TrackerTopology> tTopo = iSetup.getHandle(theTopoToken);
0165 
0166   theG = iSetup.getHandle(theGToken);
0167   edm::ESHandle<MagneticField> theMF = iSetup.getHandle(theMFToken);
0168 
0169   using namespace edm;
0170   using namespace std;
0171   using reco::TrackCollection;
0172 
0173   LogTrace("TestOutliers") << "analyze event #" << iEvent.id();
0174 
0175   edm::Handle<edm::View<reco::Track> > tracksOut;
0176   iEvent.getByLabel(trackTagsOut_, tracksOut);
0177   edm::Handle<edm::View<reco::Track> > tracksOld;
0178   iEvent.getByLabel(trackTagsOld_, tracksOld);
0179   Handle<TrackingParticleCollection> tps;
0180   iEvent.getByLabel(tpTags_, tps);
0181   edm::Handle<reco::BeamSpot> beamSpot;
0182   iEvent.getByLabel("offlineBeamSpot", beamSpot);
0183 
0184   TrackerHitAssociator hitAssociator(iEvent, trackerHitAssociatorConfig_);
0185 
0186   edm::Handle<reco::TrackToTrackingParticleAssociator> hAssociatorOld;
0187   iEvent.getByToken(theAssociatorOldToken, hAssociatorOld);
0188   const reco::TrackToTrackingParticleAssociator *theAssociatorOld = hAssociatorOld.product();
0189 
0190   edm::Handle<reco::TrackToTrackingParticleAssociator> hAssociatorOut;
0191   iEvent.getByToken(theAssociatorOutToken, hAssociatorOut);
0192   const reco::TrackToTrackingParticleAssociator *theAssociatorOut = hAssociatorOut.product();
0193 
0194   reco::RecoToSimCollection recSimCollOut = theAssociatorOut->associateRecoToSim(tracksOut, tps);
0195   reco::RecoToSimCollection recSimCollOld = theAssociatorOld->associateRecoToSim(tracksOld, tps);
0196   sizeOut->Fill(recSimCollOut.size());
0197   sizeOld->Fill(recSimCollOld.size());
0198   sizeOutT->Fill(tracksOut->size());
0199   sizeOldT->Fill(tracksOld->size());
0200 
0201   LogTrace("TestOutliers") << "tOld size=" << tracksOld->size() << " tOut size=" << tracksOut->size()
0202                            << " aOld size=" << recSimCollOld.size() << " aOut size=" << recSimCollOut.size();
0203 
0204 #if 0
0205   LogTrace("TestOutliers") << "recSimCollOld.size()=" << recSimCollOld.size() ;
0206   for(reco::TrackCollection::size_type j=0; j<tracksOld->size(); ++j){
0207     reco::TrackRef trackOld(tracksOld, j);
0208     //if ( !selectRecoTracks( *trackOld,beamSpot.product() ) ) continue;
0209     LogTrace("TestOutliers") << "trackOld->pt()=" << trackOld->pt() << " trackOld->numberOfValidHits()=" << trackOld->numberOfValidHits();
0210     std::vector<pair<TrackingParticleRef, double> > tpOld;
0211     if(recSimCollOld.find(trackOld) != recSimCollOld.end()){
0212       tpOld = recSimCollOld[trackOld];
0213       if (tpOld.size()!=0) LogTrace("TestOutliers") << " associated" ;
0214       else LogTrace("TestOutliers") << " NOT associated" ;
0215     } else LogTrace("TestOutliers") << " NOT associated" ;
0216   }
0217   LogTrace("TestOutliers") << "recSimCollOut.size()=" << recSimCollOut.size() ;
0218   for(reco::TrackCollection::size_type j=0; j<tracksOut->size(); ++j){
0219     reco::TrackRef trackOut(tracksOut, j);
0220     //if ( !selectRecoTracks( *trackOut,beamSpot.product() ) ) continue;
0221     LogTrace("TestOutliers") << "trackOut->pt()=" << trackOut->pt() << " trackOut->numberOfValidHits()=" << trackOut->numberOfValidHits();
0222     std::vector<pair<TrackingParticleRef, double> > tpOut;
0223     if(recSimCollOut.find(trackOut) != recSimCollOut.end()){
0224       tpOut = recSimCollOut[trackOut];
0225       if (tpOut.size()!=0) LogTrace("TestOutliers") << " associated" ;
0226       else LogTrace("TestOutliers") << " NOT associated" ;
0227     } else LogTrace("TestOutliers") << " NOT associated" ;
0228   }
0229 #endif
0230 
0231   LogTrace("TestOutliers") << "tracksOut->size()=" << tracksOut->size();
0232   LogTrace("TestOutliers") << "tracksOld->size()=" << tracksOld->size();
0233 
0234   std::vector<unsigned int> outused;
0235   for (unsigned int j = 0; j < tracksOld->size(); ++j) {
0236     countOldT->Fill(1);
0237     edm::RefToBase<reco::Track> trackOld(tracksOld, j);
0238     LogTrace("TestOutliers") << "now track old with id=" << j << " seed ref=" << trackOld->seedRef().get()
0239                              << " pt=" << trackOld->pt() << " eta=" << trackOld->eta()
0240                              << " chi2=" << trackOld->normalizedChi2()
0241                              << " tip= " << fabs(trackOld->dxy(beamSpot->position()))
0242                              << " lip= " << fabs(trackOld->dsz(beamSpot->position()));
0243 
0244     //     if (i==tracksOut->size()) {
0245     //       if(recSimCollOld.find(trackOld) != recSimCollOld.end()){
0246     //  vector<pair<TrackingParticleRef, double> > tpOld;
0247     //  tpOld = recSimCollOld[trackOld];
0248     //  if (tpOld.size()!=0) {
0249     //    LogTrace("TestOutliers") <<"no match: old associated and out lost! old has #hits=" << trackOld->numberOfValidHits()
0250     //                 << " and fraction=" << tpOld.begin()->second;
0251     //    if (tpOld.begin()->second>0.5) hitsPerTrackLost->Fill(trackOld->numberOfValidHits());
0252     //  }
0253     //       }
0254     //       continue;
0255     //     }
0256 
0257     std::vector<unsigned int> outtest;  //all the out tracks with the same seed ref
0258     for (unsigned int k = 0; k < tracksOut->size(); ++k) {
0259       edm::RefToBase<reco::Track> tmpOut = edm::RefToBase<reco::Track>(tracksOut, k);
0260       if (tmpOut->seedRef() == trackOld->seedRef()) {
0261         outtest.push_back(k);
0262       }
0263     }
0264 
0265     edm::RefToBase<reco::Track> trackOut;
0266     if (outtest.size() == 1) {  // if only one that's it
0267       trackOut = edm::RefToBase<reco::Track>(tracksOut, outtest[0]);
0268       LogTrace("TestOutliers") << "now track out with id=" << outtest[0] << " seed ref=" << trackOut->seedRef().get()
0269                                << " pt=" << trackOut->pt();
0270       outused.push_back(outtest[0]);
0271     } else if (outtest.size() > 1) {  //if > 1 take the one that shares all the hits with the old track
0272       for (unsigned int p = 0; p < outtest.size(); ++p) {
0273         edm::RefToBase<reco::Track> tmpOut = edm::RefToBase<reco::Track>(tracksOut, outtest[p]);
0274         bool allhits = true;
0275         for (trackingRecHit_iterator itOut = tmpOut->recHitsBegin(); itOut != tmpOut->recHitsEnd(); itOut++) {
0276           if ((*itOut)->isValid()) {
0277             bool thishit = false;
0278             for (trackingRecHit_iterator itOld = trackOld->recHitsBegin(); itOld != trackOld->recHitsEnd(); itOld++) {
0279               if ((*itOld)->isValid()) {
0280                 const TrackingRecHit *kt = &(**itOld);
0281                 if ((*itOut)->sharesInput(kt, TrackingRecHit::some)) {
0282                   thishit = true;
0283                   break;
0284                 }
0285               }
0286             }
0287             if (!thishit)
0288               allhits = false;
0289           }
0290         }
0291         if (allhits) {
0292           trackOut = edm::RefToBase<reco::Track>(tracksOut, outtest[p]);
0293           LogTrace("TestOutliers") << "now track out with id=" << outtest[p]
0294                                    << " seed ref=" << trackOut->seedRef().get() << " pt=" << trackOut->pt();
0295           outused.push_back(outtest[p]);
0296         }
0297       }
0298     }
0299 
0300     if (outtest.empty() || trackOut.get() == nullptr) {  //no out track found for the old track
0301       if (recSimCollOld.find(trackOld) != recSimCollOld.end()) {
0302         vector<pair<TrackingParticleRef, double> > tpOld;
0303         tpOld = recSimCollOld[trackOld];
0304         if (!tpOld.empty()) {
0305           LogTrace("TestOutliers") << "no match: old associated and out lost! old has #hits="
0306                                    << trackOld->numberOfValidHits() << " and fraction=" << tpOld.begin()->second;
0307           if (tpOld.begin()->second > 0.5)
0308             hitsPerTrackLost->Fill(trackOld->numberOfValidHits());
0309         }
0310       }
0311       LogTrace("TestOutliers") << "...skip to next old track";
0312       continue;
0313     }
0314 
0315     //look if old and out are associated
0316     LogTrace("TestOutliers") << "trackOut->seedRef()=" << trackOut->seedRef().get()
0317                              << " trackOld->seedRef()=" << trackOld->seedRef().get();
0318     bool oldAssoc = recSimCollOld.find(trackOld) != recSimCollOld.end();
0319     bool outAssoc = recSimCollOut.find(trackOut) != recSimCollOut.end();
0320     LogTrace("TestOutliers") << "outAssoc=" << outAssoc << " oldAssoc=" << oldAssoc;
0321 
0322     //     if ( trackOut->seedRef()!= trackOld->seedRef() ||
0323     //   (trackOut->seedRef() == trackOld->seedRef() && trackOut->numberOfValidHits()>trackOld->numberOfValidHits()) ) {
0324     //       LogTrace("TestOutliers") <<"out and old tracks does not match...";
0325     //       LogTrace("TestOutliers") <<"old has #hits=" << trackOld->numberOfValidHits();
0326     //       std::vector<pair<TrackingParticleRef, double> > tpOld;
0327     //       if(recSimCollOld.find(trackOld) != recSimCollOld.end()) {
0328     //  tpOld = recSimCollOld[trackOld];
0329     //  if (tpOld.size()!=0) {
0330     //    LogTrace("TestOutliers") <<"old was associated with fraction=" << tpOld.begin()->second;
0331     //    if (tpOld.begin()->second>0.5) hitsPerTrackLost->Fill(trackOld->numberOfValidHits());
0332     //  }
0333     //       }
0334     //       LogTrace("TestOutliers") <<"...skip to next old track";
0335     //       continue;
0336     //     }
0337     //     //++i;
0338     countOutT->Fill(1);
0339 
0340     //if ( !selectRecoTracks( *trackOld,beamSpot.product() ) ) continue;//no more cuts
0341 
0342     tracks->Fill(0);  //FIXME
0343 
0344     TrackingParticleRef tpr;
0345     TrackingParticleRef tprOut;
0346     TrackingParticleRef tprOld;
0347     double fracOut;
0348     std::vector<unsigned int> tpids;
0349     std::vector<std::pair<TrackingParticleRef, double> > tpOut;
0350     std::vector<pair<TrackingParticleRef, double> > tpOld;
0351     //contare outliers delle tracce che erano associate e non lo sono piu!!!!
0352 
0353     if (outAssoc) {  //save the ids od the tp associate to the out track
0354       tpOut = recSimCollOut[trackOut];
0355       if (!tpOut.empty()) {
0356         countOutA->Fill(1);
0357         tprOut = tpOut.begin()->first;
0358         fracOut = tpOut.begin()->second;
0359         for (TrackingParticle::g4t_iterator g4T = tprOut->g4Track_begin(); g4T != tprOut->g4Track_end(); ++g4T) {
0360           tpids.push_back(g4T->trackId());
0361         }
0362       }
0363     }
0364 
0365     if (oldAssoc) {  //save the ids od the tp associate to the old track
0366       tpOld = recSimCollOld[trackOld];
0367       if (!tpOld.empty()) {
0368         tprOld = tpOld.begin()->first;
0369         //  LogTrace("TestOutliers") <<"old associated and out not! old has #hits=" << trackOld->numberOfValidHits()
0370         //               << " and fraction=" << tpOld.begin()->second;
0371         //  if (tpOld.begin()->second>0.5) hitsPerTrackLost->Fill(trackOld->numberOfValidHits());//deve essere un plot diverso tipo LostAssoc
0372         if (tpOut.empty()) {
0373           for (TrackingParticle::g4t_iterator g4T = tprOld->g4Track_begin(); g4T != tprOld->g4Track_end(); ++g4T) {
0374             tpids.push_back(g4T->trackId());
0375           }
0376         }
0377       }
0378     }
0379 
0380     if (tprOut.get() != nullptr || tprOld.get() != nullptr) {  //at least one of the tracks has to be associated
0381 
0382       tpr = tprOut.get() != nullptr ? tprOut : tprOld;
0383 
0384       const SimTrack *assocTrack = &(*tpr->g4Track_begin());
0385 
0386       //if ( trackOut->numberOfValidHits() < trackOld->numberOfValidHits() ) {
0387       if (trackOut->numberOfValidHits() != trackOld->numberOfValidHits() ||
0388           !(*trackOut->recHitsBegin())->sharesInput((*trackOld->recHitsBegin()), TrackingRecHit::some) ||
0389           !(*(trackOut->recHitsEnd() - 1))
0390                ->sharesInput(
0391                    (*(trackOld->recHitsEnd() - 1)),
0392                    TrackingRecHit::
0393                        some)) {  //there are outliers if the number of valid hits is != or if the first and last hit does not match
0394         LogTrace("TestOutliers") << "outliers for track with #hits=" << trackOut->numberOfValidHits();
0395         tracks->Fill(1);
0396         LogTrace("TestOutliers")
0397             << "Out->pt=" << trackOut->pt() << " Old->pt=" << trackOld->pt() << " tp->pt="
0398             << sqrt(tpr->momentum().perp2())
0399             //<< " trackOut->ptError()=" << trackOut->ptError() << " trackOld->ptError()=" << trackOld->ptError()
0400             << " Old->validHits=" << trackOld->numberOfValidHits() << " Out->validHits="
0401             << trackOut->numberOfValidHits()
0402             /*<< " fracOld=" << fracOld*/
0403             << " fracOut=" << fracOut << " deltaHits=" << trackOld->numberOfValidHits() - trackOut->numberOfValidHits();
0404 
0405         //compute all the track parameters
0406         double PtPullOut = (trackOut->pt() - sqrt(tpr->momentum().perp2())) / trackOut->ptError();
0407         double PtPullOld = (trackOld->pt() - sqrt(tpr->momentum().perp2())) / trackOld->ptError();
0408         histoPtOut->Fill(PtPullOut);
0409         histoPtOld->Fill(PtPullOld);
0410 
0411         //LogTrace("TestOutliers") << "MagneticField";
0412         FreeTrajectoryState ftsAtProduction(
0413             GlobalPoint(tpr->vertex().x(), tpr->vertex().y(), tpr->vertex().z()),
0414             GlobalVector(assocTrack->momentum().x(), assocTrack->momentum().y(), assocTrack->momentum().z()),
0415             TrackCharge(trackOld->charge()),
0416             theMF.product());
0417         TSCPBuilderNoMaterial tscpBuilder;
0418         TrajectoryStateClosestToPoint tsAtClosestApproach =
0419             tscpBuilder(ftsAtProduction, GlobalPoint(0, 0, 0));  //as in TrackProducerAlgorithm
0420         GlobalPoint v = tsAtClosestApproach.theState().position();
0421         GlobalVector p = tsAtClosestApproach.theState().momentum();
0422 
0423         //LogTrace("TestOutliers") << "qoverpSim";
0424         double qoverpSim = tsAtClosestApproach.charge() / p.mag();
0425         double lambdaSim = M_PI / 2 - p.theta();
0426         double phiSim = p.phi();
0427         double dxySim = (-v.x() * sin(p.phi()) + v.y() * cos(p.phi()));
0428         double dszSim = v.z() * p.perp() / p.mag() - (v.x() * p.x() + v.y() * p.y()) / p.perp() * p.z() / p.mag();
0429         double d0Sim = -dxySim;
0430         double dzSim = dszSim * p.mag() / p.perp();
0431 
0432         //LogTrace("TestOutliers") << "qoverpPullOut";
0433         double qoverpPullOut = (trackOut->qoverp() - qoverpSim) / trackOut->qoverpError();
0434         double qoverpPullOld = (trackOld->qoverp() - qoverpSim) / trackOld->qoverpError();
0435         double lambdaPullOut = (trackOut->lambda() - lambdaSim) / trackOut->thetaError();
0436         double lambdaPullOld = (trackOld->lambda() - lambdaSim) / trackOld->thetaError();
0437         double phi0PullOut = (trackOut->phi() - phiSim) / trackOut->phiError();
0438         double phi0PullOld = (trackOld->phi() - phiSim) / trackOld->phiError();
0439         double d0PullOut = (trackOut->d0() - d0Sim) / trackOut->d0Error();
0440         double d0PullOld = (trackOld->d0() - d0Sim) / trackOld->d0Error();
0441         double dzPullOut = (trackOut->dz() - dzSim) / trackOut->dzError();
0442         double dzPullOld = (trackOld->dz() - dzSim) / trackOld->dzError();
0443 
0444         //LogTrace("TestOutliers") << "histoQoverpOut";
0445         histoQoverpOut->Fill(qoverpPullOut);
0446         histoQoverpOld->Fill(qoverpPullOld);
0447         histoPhiOut->Fill(phi0PullOut);
0448         histoPhiOld->Fill(phi0PullOld);
0449         histoD0Out->Fill(d0PullOut);
0450         histoD0Old->Fill(d0PullOld);
0451         histoDzOut->Fill(dzPullOut);
0452         histoDzOld->Fill(dzPullOld);
0453         histoLambdaOut->Fill(lambdaPullOut);
0454         histoLambdaOld->Fill(lambdaPullOld);
0455 
0456         //delta number of valid hits
0457         LogTrace("TestOutliers") << "deltahits=" << trackOld->numberOfValidHits() - trackOut->numberOfValidHits();
0458         deltahits->Fill(trackOld->numberOfValidHits() - trackOut->numberOfValidHits());
0459 
0460         if (tprOut.get() != nullptr && tprOld.get() == nullptr) {  //out associated and old not: gained track
0461           if (!tpOld.empty() && tpOld.begin()->second <= 0.5) {
0462             deltahitsAssocGained->Fill(trackOld->numberOfValidHits() - trackOut->numberOfValidHits());
0463             hitsPerTrackAssocGained->Fill(trackOut->numberOfValidHits());
0464             LogTrace("TestOutliers") << "a) gained (assoc) track out #hits==" << trackOut->numberOfValidHits()
0465                                      << " old #hits=" << trackOld->numberOfValidHits();
0466           } else {
0467             deltahitsAssocGained->Fill(trackOld->numberOfValidHits() - trackOut->numberOfValidHits());
0468             hitsPerTrackAssocGained->Fill(trackOut->numberOfValidHits());
0469             LogTrace("TestOutliers") << "b) gained (assoc) track out #hits==" << trackOut->numberOfValidHits()
0470                                      << " old #hits=" << trackOld->numberOfValidHits();
0471           }
0472         } else if (tprOut.get() == nullptr && tprOld.get() != nullptr) {  //old associated and out not: lost track
0473           LogTrace("TestOutliers") << "old associated and out not! old has #hits=" << trackOld->numberOfValidHits()
0474                                    << " and fraction=" << tpOld.begin()->second;
0475           if (tpOld.begin()->second > 0.5) {
0476             hitsPerTrackAssocLost->Fill(trackOld->numberOfValidHits());
0477             deltahitsAssocLost->Fill(trackOld->numberOfValidHits() - trackOut->numberOfValidHits());
0478           }
0479         }
0480 
0481         if (fabs(PtPullOut) < fabs(PtPullOld))
0482           deltahitsOK->Fill(trackOld->numberOfValidHits() - trackOut->numberOfValidHits());
0483         else
0484           deltahitsNO->Fill(trackOld->numberOfValidHits() - trackOut->numberOfValidHits());
0485 
0486         //LogTrace("TestOutliers") << "RecoTrackSelector";
0487         //if (selectRecoTracks(*trackOut,beamSpot.product())) okcutsOut->Fill(1); else okcutsOut->Fill(0);
0488         //if (selectRecoTracks(*trackOld,beamSpot.product())) okcutsOld->Fill(1); else okcutsOld->Fill(0);
0489 
0490         LogTrace("TestOutliers") << "track old";
0491         for (trackingRecHit_iterator itOld = trackOld->recHitsBegin(); itOld != trackOld->recHitsEnd(); itOld++) {
0492           LogTrace("TestOutliers") << (*itOld)->isValid() << " " << (*itOld)->geographicalId().rawId();
0493         }
0494         LogTrace("TestOutliers") << "track out";
0495         for (trackingRecHit_iterator itOut = trackOut->recHitsBegin(); itOut != trackOut->recHitsEnd(); itOut++) {
0496           LogTrace("TestOutliers") << (*itOut)->isValid() << " " << (*itOut)->geographicalId().rawId();
0497         }
0498         //LogTrace("TestOutliers") << "itOut";
0499 
0500         vector<pair<int, trackingRecHit_iterator> > gainedlostoutliers;
0501         //look for gained hits
0502         for (trackingRecHit_iterator itOut = trackOut->recHitsBegin(); itOut != trackOut->recHitsEnd(); itOut++) {
0503           bool gained = true;
0504           if ((*itOut)->isValid()) {
0505             for (trackingRecHit_iterator itOld = trackOld->recHitsBegin(); itOld != trackOld->recHitsEnd(); itOld++) {
0506               if ((*itOld)->geographicalId().rawId() == (*itOut)->geographicalId().rawId())
0507                 gained = false;
0508             }
0509             if (gained) {
0510               gainedlostoutliers.push_back(pair<int, trackingRecHit_iterator>(1, itOut));
0511               LogTrace("TestOutliers") << "broken trajectory during old fit... gained hit "
0512                                        << (*itOut)->geographicalId().rawId();
0513               gainedhits->Fill(1);
0514             }
0515           }
0516         }
0517 
0518         //look for outliers and lost hits
0519         for (trackingRecHit_iterator itOld = trackOld->recHitsBegin(); itOld != trackOld->recHitsEnd(); itOld++) {
0520           bool outlier = false;
0521           bool lost = true;
0522 
0523           for (trackingRecHit_iterator itOut = trackOut->recHitsBegin(); itOut != trackOut->recHitsEnd(); itOut++) {
0524             if ((*itOld)->geographicalId().rawId() == (*itOut)->geographicalId().rawId()) {
0525               lost = false;
0526               if ((*itOld)->isValid() && !(*itOut)->isValid() &&
0527                   (*itOld)->geographicalId().rawId() == (*itOut)->geographicalId().rawId()) {
0528                 LogTrace("TestOutliers") << (*itOld)->isValid() << " " << (*itOut)->isValid() << " "
0529                                          << (*itOld)->geographicalId().rawId() << " "
0530                                          << (*itOut)->geographicalId().rawId();
0531                 outlier = true;
0532               }
0533             }
0534           }
0535           if (lost)
0536             gainedlostoutliers.push_back(pair<int, trackingRecHit_iterator>(2, itOld));
0537           if (lost)
0538             LogTrace("TestOutliers") << "lost";
0539           else if (outlier)
0540             gainedlostoutliers.push_back(pair<int, trackingRecHit_iterator>(3, itOld));
0541         }
0542 
0543         for (std::vector<pair<int, trackingRecHit_iterator> >::iterator it = gainedlostoutliers.begin();
0544              it != gainedlostoutliers.end();
0545              ++it) {
0546           LogTrace("TestOutliers") << "type of processed hit:" << it->first;
0547           trackingRecHit_iterator itHit = it->second;
0548           bool gained = false;
0549           bool lost = false;
0550           bool outlier = false;
0551           if (it->first == 1)
0552             gained = true;
0553           else if (it->first == 2)
0554             lost = true;
0555           else if (it->first == 3)
0556             outlier = true;
0557 
0558           if (outlier || lost || gained) {
0559             //if (1) {
0560 
0561             if (lost && (*itHit)->isValid() == false) {
0562               goodbadmergedLost->Fill(0);
0563               LogTrace("TestOutliers") << "lost invalid";
0564               continue;
0565             } else if (gained && (*itHit)->isValid() == false) {
0566               goodbadmergedGained->Fill(0);
0567               LogTrace("TestOutliers") << "gained invalid";
0568               continue;
0569             }
0570 
0571             //LogTrace("TestOutliers") << "vector<SimHitIdpr>";
0572             //look if the hit comes from a correct sim track
0573             std::vector<SimHitIdpr> simTrackIds = hitAssociator.associateHitId(**itHit);
0574             bool goodhit = false;
0575             for (size_t j = 0; j < simTrackIds.size(); j++) {
0576               for (size_t jj = 0; jj < tpids.size(); jj++) {
0577                 if (simTrackIds[j].first == tpids[jj])
0578                   goodhit = true;
0579                 break;
0580               }
0581             }
0582 
0583             //find what kind of hit is it
0584             int clustersize = 0;
0585             int hittypeval = 0;
0586             int layerval = 0;
0587             if (dynamic_cast<const SiPixelRecHit *>(&**itHit)) {
0588               LogTrace("TestOutliers") << "SiPixelRecHit";
0589               clustersize = ((const SiPixelRecHit *)(&**itHit))->cluster()->size();
0590               hittypeval = 1;
0591             } else if (dynamic_cast<const SiStripRecHit2D *>(&**itHit)) {
0592               LogTrace("TestOutliers") << "SiStripRecHit2D";
0593               clustersize = ((const SiStripRecHit2D *)(&**itHit))->cluster()->amplitudes().size();
0594               hittypeval = 2;
0595             } else if (dynamic_cast<const SiStripMatchedRecHit2D *>(&**itHit)) {
0596               LogTrace("TestOutliers") << "SiStripMatchedRecHit2D";
0597               int clsize1 = ((const SiStripMatchedRecHit2D *)(&**itHit))->monoCluster().amplitudes().size();
0598               int clsize2 = ((const SiStripMatchedRecHit2D *)(&**itHit))->stereoCluster().amplitudes().size();
0599               if (clsize1 > clsize2)
0600                 clustersize = clsize1;
0601               else
0602                 clustersize = clsize2;
0603               hittypeval = 3;
0604             } else if (dynamic_cast<const ProjectedSiStripRecHit2D *>(&**itHit)) {
0605               LogTrace("TestOutliers") << "ProjectedSiStripRecHit2D";
0606               clustersize =
0607                   ((const ProjectedSiStripRecHit2D *)(&**itHit))->originalHit().cluster()->amplitudes().size();
0608               hittypeval = 4;
0609             }
0610 
0611             //find the layer of the hit
0612             int subdetId = (*itHit)->geographicalId().subdetId();
0613             DetId id = (*itHit)->geographicalId();
0614             int layerId = tTopo->layer(id);
0615             layerval = subdetId * 10 + layerId;
0616 
0617             //LogTrace("TestOutliers") << "gpos";
0618             GlobalPoint gpos = theG->idToDet((*itHit)->geographicalId())->surface().toGlobal((*itHit)->localPosition());
0619 
0620             //get the vector of sim hit associated and choose the one with the largest energy loss
0621             //double delta = 99999;
0622             //LocalPoint rhitLPv = (*itHit)->localPosition();
0623             //vector<PSimHit> assSimHits = hitAssociator.associateHit(**itHit);
0624             //if (assSimHits.size()==0) continue;
0625             //PSimHit shit;
0626             //for(std::vector<PSimHit>::const_iterator m=assSimHits.begin(); m<assSimHits.end(); m++){
0627             //if ((m->localPosition()-rhitLPv).mag()<delta) {
0628             //  shit=*m;
0629             //  delta = (m->localPosition()-rhitLPv).mag();
0630             // }
0631             //}
0632             //LogTrace("TestOutliers") << "energyLoss_";
0633             double energyLoss_ = 0.;
0634             unsigned int monoId = 0;
0635             std::vector<double> energyLossM;
0636             std::vector<double> energyLossS;
0637             std::vector<PSimHit> assSimHits = hitAssociator.associateHit(**itHit);
0638             if (assSimHits.empty())
0639               continue;
0640             PSimHit shit;
0641             std::vector<unsigned int> trackIds;
0642             energyLossS.clear();
0643             energyLossM.clear();
0644             //LogTrace("TestOutliers") << "energyLossM";
0645             for (std::vector<PSimHit>::const_iterator m = assSimHits.begin(); m < assSimHits.end(); m++) {
0646               if (outlier)
0647                 energyLoss->Fill(m->energyLoss());
0648               unsigned int tId = m->trackId();
0649               if (find(trackIds.begin(), trackIds.end(), tId) == trackIds.end())
0650                 trackIds.push_back(tId);
0651               LogTrace("TestOutliers") << "id=" << tId;
0652               if (m->energyLoss() > energyLoss_) {
0653                 shit = *m;
0654                 energyLoss_ = m->energyLoss();
0655               }
0656               if (hittypeval == 3) {
0657                 if (monoId == 0)
0658                   monoId = m->detUnitId();
0659                 if (monoId == m->detUnitId()) {
0660                   energyLossM.push_back(m->energyLoss());
0661                 } else {
0662                   energyLossS.push_back(m->energyLoss());
0663                 }
0664                 //std::cout << "detUnitId="  << m->detUnitId() << " trackId=" << m->trackId() << " energyLoss=" << m->energyLoss() << std::endl;
0665               } else {
0666                 energyLossM.push_back(m->energyLoss());
0667               }
0668             }
0669             unsigned int nIds = trackIds.size();
0670 
0671             if (outlier) {
0672               goodbadhits->Fill(goodhit);
0673               posxy->Fill(fabs(gpos.x()), fabs(gpos.y()));
0674               poszr->Fill(fabs(gpos.z()), sqrt(gpos.x() * gpos.x() + gpos.y() * gpos.y()));
0675               process->Fill(shit.processType());
0676               energyLossMax->Fill(energyLoss_);
0677             }
0678 
0679             //look if the hit is shared and if is produced only by ionization processes
0680             bool shared = true;
0681             bool ioniOnly = true;
0682             unsigned int idc = 0;
0683             for (size_t jj = 0; jj < tpids.size(); jj++) {
0684               idc += std::count(trackIds.begin(), trackIds.end(), tpids[jj]);
0685             }
0686             if (idc == trackIds.size()) {
0687               shared = false;
0688             }
0689             for (std::vector<PSimHit>::const_iterator m = assSimHits.begin() + 1; m < assSimHits.end(); m++) {
0690               if ((m->processType() != 7 && m->processType() != 8 && m->processType() != 9) &&
0691                   abs(m->particleType()) != 11) {
0692                 ioniOnly = false;
0693                 break;
0694               }
0695             }
0696             if (ioniOnly && !shared) {
0697               LogTrace("TestOutliers") << "delta";
0698             }
0699 
0700             if (goodhit) {
0701               if (outlier) {
0702                 goodprocess->Fill(shit.processType());
0703                 if (clustersize >= 4) {
0704                   goodhittype_clgteq4->Fill(hittypeval);
0705                   goodlayer_clgteq4->Fill(layerval);
0706                 } else {
0707                   goodhittype_cllt4->Fill(hittypeval);
0708                   goodlayer_cllt4->Fill(layerval);
0709                 }
0710                 //LogTrace("TestOutliers") << "hittypeval && clustersize";
0711                 if (hittypeval == 1 && clustersize >= 4)
0712                   goodpixgteq4_simvecsize->Fill(assSimHits.size());
0713                 if (hittypeval == 1 && clustersize < 4)
0714                   goodpixlt4_simvecsize->Fill(assSimHits.size());
0715                 if (hittypeval == 2 && clustersize >= 4)
0716                   goodst1gteq4_simvecsize->Fill(assSimHits.size());
0717                 if (hittypeval == 2 && clustersize < 4)
0718                   goodst1lt4_simvecsize->Fill(assSimHits.size());
0719                 if (hittypeval == 3 && clustersize >= 4)
0720                   goodst2gteq4_simvecsize->Fill(assSimHits.size());
0721                 if (hittypeval == 3 && clustersize < 4)
0722                   goodst2lt4_simvecsize->Fill(assSimHits.size());
0723                 if (hittypeval == 4 && clustersize >= 4)
0724                   goodprjgteq4_simvecsize->Fill(assSimHits.size());
0725                 if (hittypeval == 4 && clustersize < 4)
0726                   goodprjlt4_simvecsize->Fill(assSimHits.size());
0727 
0728                 //LogTrace("TestOutliers") << "hittypeval";
0729                 if (hittypeval == 1)
0730                   goodpix_clustersize->Fill(clustersize);
0731                 if (hittypeval == 2)
0732                   goodst1_clustersize->Fill(clustersize);
0733                 if (hittypeval == 3)
0734                   goodst2_clustersize->Fill(clustersize);
0735                 if (hittypeval == 4)
0736                   goodprj_clustersize->Fill(clustersize);
0737                 if (hittypeval == 1)
0738                   goodpix_simvecsize->Fill(assSimHits.size());
0739                 if (hittypeval == 2)
0740                   goodst1_simvecsize->Fill(assSimHits.size());
0741                 if (hittypeval == 3)
0742                   goodst2_simvecsize->Fill(assSimHits.size());
0743                 if (hittypeval == 4)
0744                   goodprj_simvecsize->Fill(assSimHits.size());
0745 
0746                 //LogTrace("TestOutliers") << "nOfTrackIds";
0747                 nOfTrackIds->Fill(nIds);
0748                 if (hittypeval != 3) {
0749                   if (energyLossM.size() > 1) {
0750                     sort(energyLossM.begin(), energyLossM.end(), greater<double>());
0751                     energyLossRatio->Fill(energyLossM[1] / energyLossM[0]);
0752                   }
0753                 } else {
0754                   //LogTrace("TestOutliers") << "hittypeval==3";
0755                   if (energyLossM.size() > 1 && energyLossS.size() <= 1) {
0756                     sort(energyLossM.begin(), energyLossM.end(), greater<double>());
0757                     energyLossRatio->Fill(energyLossM[1] / energyLossM[0]);
0758                   } else if (energyLossS.size() > 1 && energyLossM.size() <= 1) {
0759                     sort(energyLossS.begin(), energyLossS.end(), greater<double>());
0760                     energyLossRatio->Fill(energyLossS[1] / energyLossS[0]);
0761                   } else if (energyLossS.size() > 1 && energyLossM.size() > 1) {
0762                     sort(energyLossM.begin(), energyLossM.end(), greater<double>());
0763                     sort(energyLossS.begin(), energyLossS.end(), greater<double>());
0764                     energyLossM[1] / energyLossM[0] > energyLossS[1] / energyLossS[0]
0765                         ? energyLossRatio->Fill(energyLossM[1] / energyLossM[0])
0766                         : energyLossRatio->Fill(energyLossS[1] / energyLossS[0]);
0767                   }
0768                 }
0769 
0770                 LogTrace("TestOutliers") << "before merged";
0771                 const SiStripMatchedRecHit2D *tmp = dynamic_cast<const SiStripMatchedRecHit2D *>(&**itHit);
0772                 LogTrace("TestOutliers") << "tmp=" << tmp;
0773                 LogTrace("TestOutliers") << "assSimHits.size()=" << assSimHits.size();
0774                 if ((assSimHits.size() > 1 && tmp == nullptr) || (assSimHits.size() > 2 && tmp != nullptr)) {
0775                   //std::cout << "MERGED HIT" << std::endl;
0776                   //LogTrace("TestOutliers") << "merged";
0777                   mergedlayer->Fill(layerval);
0778                   mergedcluster->Fill(clustersize);
0779                   mergedhittype->Fill(hittypeval);
0780 
0781                   for (std::vector<PSimHit>::const_iterator m = assSimHits.begin(); m < assSimHits.end(); m++) {
0782                     unsigned int tId = m->trackId();
0783                     LogTrace("TestOutliers") << "component with id=" << tId << " eLoss=" << m->energyLoss()
0784                                              << " proc=" << m->processType() << " part=" << m->particleType();
0785                     if (find(tpids.begin(), tpids.end(), tId) == tpids.end())
0786                       continue;
0787                     if (m->processType() == 2) {
0788                       //GlobalPoint gpos = theG->idToDet((*itHit)->geographicalId())->surface().toGlobal((*itHit)->localPosition());
0789                       //GlobalPoint gpr = rhit->globalPosition();
0790                       AlgebraicSymMatrix33 ger = ErrorFrameTransformer()
0791                                                      .transform((*itHit)->localPositionError(),
0792                                                                 theG->idToDet((*itHit)->geographicalId())->surface())
0793                                                      .matrix();
0794                       //AlgebraicSymMatrix ger = rhit->globalPositionError().matrix();
0795                       GlobalPoint gps =
0796                           theG->idToDet((*itHit)->geographicalId())->surface().toGlobal(m->localPosition());
0797                       //LogTrace("TestOutliers") << gpr << " " << gps << " " << ger;
0798                       ROOT::Math::SVector<double, 3> delta;
0799                       delta[0] = gpos.x() - gps.x();
0800                       delta[1] = gpos.y() - gps.y();
0801                       delta[2] = gpos.z() - gps.z();
0802                       LogTrace("TestOutliers") << delta << " " << ger;
0803                       double mpull = sqrt(delta[0] * delta[0] / ger[0][0] + delta[1] * delta[1] / ger[1][1] +
0804                                           delta[2] * delta[2] / ger[2][2]);
0805                       LogTrace("TestOutliers") << "hit pull=" << mpull;  //ger.similarity(delta);
0806                       mergedPull->Fill(mpull);
0807                     }
0808                   }
0809                   LogTrace("TestOutliers") << "end merged";
0810                 } else {  //not merged=>good
0811                   //LogTrace("TestOutliers") << "goodlayer";
0812                   goodlayer->Fill(layerval);
0813                   goodcluster->Fill(clustersize);
0814                   goodhittype->Fill(hittypeval);
0815                 }
0816               }  //if (outlier)
0817 
0818               const SiPixelRecHit *pix = dynamic_cast<const SiPixelRecHit *>(&**itHit);
0819               if ((hittypeval != 3 && assSimHits.size() < 2) || (hittypeval == 3 && assSimHits.size() < 3)) {
0820                 if (outlier) {
0821                   goodhittype_simvecsmall->Fill(hittypeval);
0822                   goodlayer_simvecsmall->Fill(layerval);
0823                   goodbadmerged->Fill(1);
0824                   if (pix) {
0825                     probXgood->Fill(pix->probabilityX());
0826                     probYgood->Fill(pix->probabilityY());
0827                   }
0828                   LogTrace("TestOutliers") << "out good";
0829                 } else if (lost) {
0830                   goodbadmergedLost->Fill(1);
0831                   LogTrace("TestOutliers") << "lost good";
0832                 } else if (gained) {
0833                   goodbadmergedGained->Fill(1);
0834                   LogTrace("TestOutliers") << "gained good";
0835                 }
0836                 LogTrace("TestOutliers") << "good";
0837               } else {
0838                 if (outlier) {
0839                   goodhittype_simvecbig->Fill(hittypeval);
0840                   goodlayer_simvecbig->Fill(layerval);
0841                   if (ioniOnly && !shared) {
0842                     goodbadmerged->Fill(3);
0843                     if (pix) {
0844                       probXdelta->Fill(pix->probabilityX());
0845                       probYdelta->Fill(pix->probabilityY());
0846                     }
0847                   } else if (!ioniOnly && !shared) {
0848                     goodbadmerged->Fill(4);
0849                     if (pix) {
0850                       probXnoshare->Fill(pix->probabilityX());
0851                       probYnoshare->Fill(pix->probabilityY());
0852                     }
0853                   } else {
0854                     goodbadmerged->Fill(5);
0855                     if (pix) {
0856                       probXshared->Fill(pix->probabilityX());
0857                       probYshared->Fill(pix->probabilityY());
0858                     }
0859                   }
0860                   LogTrace("TestOutliers") << "out merged, ioniOnly=" << ioniOnly << " shared=" << shared;
0861                 } else if (lost) {
0862                   if (ioniOnly && !shared)
0863                     goodbadmergedLost->Fill(3);
0864                   else if (!ioniOnly && !shared)
0865                     goodbadmergedLost->Fill(4);
0866                   else
0867                     goodbadmergedLost->Fill(5);
0868                   LogTrace("TestOutliers") << "lost merged, ioniOnly=" << ioniOnly << " shared=" << shared;
0869                 } else if (gained) {
0870                   if (ioniOnly && !shared)
0871                     goodbadmergedGained->Fill(3);
0872                   else if (!ioniOnly && !shared)
0873                     goodbadmergedGained->Fill(4);
0874                   else
0875                     goodbadmergedGained->Fill(5);
0876                   LogTrace("TestOutliers") << "gained merged, ioniOnly=" << ioniOnly << " shared=" << shared;
0877                 }
0878                 LogTrace("TestOutliers") << "merged, ioniOnly=" << ioniOnly << " shared=" << shared;
0879               }
0880             }       //if (goodhit)
0881             else {  //badhit
0882               //LogTrace("TestOutliers") << "badhit";
0883               if (outlier) {
0884                 badcluster->Fill(clustersize);
0885                 badhittype->Fill(hittypeval);
0886                 badlayer->Fill(layerval);
0887                 badprocess->Fill(shit.processType());
0888                 goodbadmerged->Fill(2);
0889                 const SiPixelRecHit *pix = dynamic_cast<const SiPixelRecHit *>(&**itHit);
0890                 if (pix) {
0891                   probXbad->Fill(pix->probabilityX());
0892                   probYbad->Fill(pix->probabilityY());
0893                 }
0894                 LogTrace("TestOutliers") << "out bad";
0895               } else if (lost) {
0896                 goodbadmergedLost->Fill(2);
0897                 LogTrace("TestOutliers") << "lost bad";
0898               } else if (gained) {
0899                 goodbadmergedGained->Fill(2);
0900                 LogTrace("TestOutliers") << "gained bad";
0901               }
0902               LogTrace("TestOutliers") << "bad";
0903             }
0904           }
0905         }
0906       }
0907       //else if ( trackOut->numberOfValidHits() > trackOld->numberOfValidHits() ) {
0908       else if (false) {
0909         LogTrace("TestOutliers") << "outliers for track with #hits=" << trackOut->numberOfValidHits();
0910         tracks->Fill(1);
0911         LogTrace("TestOutliers")
0912             << "Out->pt=" << trackOut->pt() << " Old->pt=" << trackOld->pt() << " tp->pt="
0913             << sqrt(tpr->momentum().perp2())
0914             //<< " trackOut->ptError()=" << trackOut->ptError() << " trackOld->ptError()=" << trackOld->ptError()
0915             << " Old->validHits=" << trackOld->numberOfValidHits() << " Out->validHits="
0916             << trackOut->numberOfValidHits()
0917             /*<< " fracOld=" << fracOld*/
0918             << " fracOut=" << fracOut << " deltaHits=" << trackOld->numberOfValidHits() - trackOut->numberOfValidHits();
0919         LogTrace("TestOutliers") << "track with gained hits";
0920         gainedhits2->Fill(trackOut->numberOfValidHits() - trackOld->numberOfValidHits());
0921       } else {
0922         LogTrace("TestOutliers") << "no outliers for track with #hits=" << trackOut->numberOfValidHits();
0923       }
0924     }
0925     LogTrace("TestOutliers") << "end track old #" << j;
0926   }
0927 
0928   for (unsigned int k = 0; k < tracksOut->size(); ++k) {
0929     if (find(outused.begin(), outused.end(), k) == outused.end()) {
0930       edm::RefToBase<reco::Track> trackOut(tracksOut, k);
0931       bool outAssoc = recSimCollOut.find(trackOut) != recSimCollOut.end();
0932       if (outAssoc) {
0933         hitsPerTrackGained->Fill(trackOut->numberOfValidHits());
0934         LogTrace("TestOutliers") << "gained track out id=" << k << " #hits==" << trackOut->numberOfValidHits();
0935       }
0936     }
0937   }
0938 }
0939 
0940 // ------------ method called once each job just before starting event loop  ------------
0941 void TestOutliers::beginRun(edm::Run const &run, const edm::EventSetup &es) {
0942   const bool oldAddDir = TH1::AddDirectoryStatus();
0943   TH1::AddDirectory(true);
0944   file = new TFile(out.c_str(), "recreate");
0945   histoPtOut = new TH1F("histoPtOut", "histoPtOut", 100, -10, 10);
0946   histoPtOld = new TH1F("histoPtOld", "histoPtOld", 100, -10, 10);
0947   histoQoverpOut = new TH1F("histoQoverpOut", "histoQoverpOut", 250, -25, 25);
0948   histoQoverpOld = new TH1F("histoQoverpOld", "histoQoverpOld", 250, -25, 25);
0949   histoPhiOut = new TH1F("histoPhiOut", "histoPhiOut", 250, -25, 25);
0950   histoPhiOld = new TH1F("histoPhiOld", "histoPhiOld", 250, -25, 25);
0951   histoD0Out = new TH1F("histoD0Out", "histoD0Out", 250, -25, 25);
0952   histoD0Old = new TH1F("histoD0Old", "histoD0Old", 250, -25, 25);
0953   histoDzOut = new TH1F("histoDzOut", "histoDzOut", 250, -25, 25);
0954   histoDzOld = new TH1F("histoDzOld", "histoDzOld", 250, -25, 25);
0955   histoLambdaOut = new TH1F("histoLambdaOut", "histoLambdaOut", 250, -25, 25);
0956   histoLambdaOld = new TH1F("histoLambdaOld", "histoLambdaOld", 250, -25, 25);
0957   deltahits = new TH1F("deltahits", "deltahits", 80, -40, 40);
0958   deltahitsOK = new TH1F("deltahitsOK", "deltahitsOK", 20, 0, 20);
0959   deltahitsNO = new TH1F("deltahitsNO", "deltahitsNO", 20, 0, 20);
0960   okcutsOut = new TH1F("okcutsOut", "okcutsOut", 2, -0.5, 1.5);
0961   okcutsOld = new TH1F("okcutsOld", "okcutsOld", 2, -0.5, 1.5);
0962   tracks = new TH1F("tracks_", "tracks_", 2, -0.5, 1.5);
0963   goodbadhits = new TH1F("goodbadhits", "goodbadhits", 2, -0.5, 1.5);
0964   process = new TH1F("process", "process", 20, -0.5, 19.5);
0965   goodcluster = new TH1F("goodcluster", "goodcluster", 40, -0.5, 39.5);
0966   goodprocess = new TH1F("goodprocess", "goodprocess", 20, -0.5, 19.5);
0967   badcluster = new TH1F("badcluster", "badcluster", 40, -0.5, 39.5);
0968   badprocess = new TH1F("badprocess", "badprocess", 20, -0.5, 19.5);
0969   goodhittype = new TH1F("goodhittype", "goodhittype", 5, -0.5, 4.5);
0970   goodlayer = new TH1F("goodlayer", "goodlayer", 70, -0.5, 69.5);
0971   goodhittype_clgteq4 = new TH1F("goodhittype_clgteq4", "goodhittype_clgteq4", 5, -0.5, 4.5);
0972   goodlayer_clgteq4 = new TH1F("goodlayer_clgteq4", "goodlayer_clgteq4", 70, -0.5, 69.5);
0973   goodhittype_cllt4 = new TH1F("goodhittype_cllt4", "goodhittype_cllt4", 5, -0.5, 4.5);
0974   goodlayer_cllt4 = new TH1F("goodlayer_cllt4", "goodlayer_cllt4", 70, -0.5, 69.5);
0975   badhittype = new TH1F("badhittype", "badhittype", 5, -0.5, 4.5);
0976   badlayer = new TH1F("badlayer", "badlayer", 70, -0.5, 69.5);
0977   posxy = new TH2F("posxy", "posxy", 1200, 0, 120, 1200, 0, 120);
0978   poszr = new TH2F("poszr", "poszr", 3000, 0, 300, 1200, 0, 120);
0979   goodpixgteq4_simvecsize = new TH1F("goodpixgteq4_simvecsize", "goodpixgteq4_simvecsize", 40, -0.5, 39.5);
0980   goodpixlt4_simvecsize = new TH1F("goodpixlt4_simvecsize", "goodpixlt4_simvecsize", 40, -0.5, 39.5);
0981   goodst1gteq4_simvecsize = new TH1F("goodst1gteq4_simvecsize", "goodst1gteq4_simvecsize", 40, -0.5, 39.5);
0982   goodst1lt4_simvecsize = new TH1F("goodst1lt4_simvecsize", "goodst1lt4_simvecsize", 40, -0.5, 39.5);
0983   goodst2gteq4_simvecsize = new TH1F("goodst2gteq4_simvecsize", "goodst2gteq4_simvecsize", 40, -0.5, 39.5);
0984   goodst2lt4_simvecsize = new TH1F("goodst2lt4_simvecsize", "goodst2lt4_simvecsize", 40, -0.5, 39.5);
0985   goodprjgteq4_simvecsize = new TH1F("goodprjgteq4_simvecsize", "goodprjgteq4_simvecsize", 40, -0.5, 39.5);
0986   goodprjlt4_simvecsize = new TH1F("goodprjlt4_simvecsize", "goodprjlt4_simvecsize", 40, -0.5, 39.5);
0987   goodpix_clustersize = new TH1F("goodpix_clustersize", "goodpix_clustersize", 40, -0.5, 39.5);
0988   goodst1_clustersize = new TH1F("goodst1_clustersize", "goodst1_clustersize", 40, -0.5, 39.5);
0989   goodst2_clustersize = new TH1F("goodst2_clustersize", "goodst2_clustersize", 40, -0.5, 39.5);
0990   goodprj_clustersize = new TH1F("goodprj_clustersize", "goodprj_clustersize", 40, -0.5, 39.5);
0991   goodpix_simvecsize = new TH1F("goodpix_simvecsize", "goodpix_simvecsize", 40, -0.5, 39.5);
0992   goodst1_simvecsize = new TH1F("goodst1_simvecsize", "goodst1_simvecsize", 40, -0.5, 39.5);
0993   goodst2_simvecsize = new TH1F("goodst2_simvecsize", "goodst2_simvecsize", 40, -0.5, 39.5);
0994   goodprj_simvecsize = new TH1F("goodprj_simvecsize", "goodprj_simvecsize", 40, -0.5, 39.5);
0995   goodhittype_simvecsmall = new TH1F("goodhittype_simvecsmall", "goodhittype_simvecsmall", 5, -0.5, 4.5);
0996   goodlayer_simvecsmall = new TH1F("goodlayer_simvecsmall", "goodlayer_simvecsmall", 70, -0.5, 69.5);
0997   goodhittype_simvecbig = new TH1F("goodhittype_simvecbig", "goodhittype_simvecbig", 5, -0.5, 4.5);
0998   goodlayer_simvecbig = new TH1F("goodlayer_simvecbig", "goodlayer_simvecbig", 70, -0.5, 69.5);
0999   goodbadmerged = new TH1F("goodbadmerged", "goodbadmerged", 5, 0.5, 5.5);
1000   goodbadmergedLost = new TH1F("goodbadmergedLost", "goodbadmergedLost", 5, 0.5, 5.5);
1001   goodbadmergedGained = new TH1F("goodbadmergedGained", "goodbadmergedGained", 5, 0.5, 5.5);
1002   energyLoss = new TH1F("energyLoss", "energyLoss", 1000, 0, 0.1);
1003   energyLossMax = new TH1F("energyLossMax", "energyLossMax", 1000, 0, 0.1);
1004   energyLossRatio = new TH1F("energyLossRatio", "energyLossRatio", 100, 0, 1);
1005   nOfTrackIds = new TH1F("nOfTrackIds", "nOfTrackIds", 10, 0, 10);
1006   mergedPull = new TH1F("mergedPull", "mergedPull", 100, 0, 10);
1007   mergedlayer = new TH1F("mergedlayer", "mergedlayer", 70, -0.5, 69.5);
1008   mergedhittype = new TH1F("mergedhittype", "mergedhittype", 5, -0.5, 4.5);
1009   mergedcluster = new TH1F("mergedcluster", "mergedcluster", 40, -0.5, 39.5);
1010   deltahitsAssocGained = new TH1F("deltahitsAssocGained", "deltahitsAssocGained", 80, -40, 40);
1011   deltahitsAssocLost = new TH1F("deltahitsAssocLost", "deltahitsAssocLost", 80, -40, 40);
1012   hitsPerTrackLost = new TH1F("hitsPerTrackLost", "hitsPerTrackLost", 40, -0.5, 39.5);
1013   hitsPerTrackAssocLost = new TH1F("hitsPerTrackAssocLost", "hitsPerTrackAssocLost", 40, -0.5, 39.5);
1014   hitsPerTrackGained = new TH1F("hitsPerTrackGained", "hitsPerTrackGained", 40, -0.5, 39.5);
1015   hitsPerTrackAssocGained = new TH1F("hitsPerTrackAssocGained", "hitsPerTrackAssocGained", 40, -0.5, 39.5);
1016   sizeOut = new TH1F("sizeOut", "sizeOut", 900, -0.5, 899.5);
1017   sizeOld = new TH1F("sizeOld", "sizeOld", 900, -0.5, 899.5);
1018   sizeOutT = new TH1F("sizeOutT", "sizeOutT", 900, -0.5, 899.5);
1019   sizeOldT = new TH1F("sizeOldT", "sizeOldT", 900, -0.5, 899.5);
1020   countOutA = new TH1F("countOutA", "countOutA", 2, 0, 2);
1021   countOutT = new TH1F("countOutT", "countOutT", 2, 0, 2);
1022   countOldT = new TH1F("countOldT", "countOldT", 2, 0, 2);
1023   gainedhits = new TH1F("gainedhits", "gainedhits", 2, 0, 2);
1024   gainedhits2 = new TH1F("gainedhits2", "gainedhits2", 30, -0.5, 29.5);
1025   probXgood = new TH1F("probXgood", "probXgood", 110, 0, 1.1);
1026   probXbad = new TH1F("probXbad", "probXbad", 110, 0, 1.1);
1027   probXdelta = new TH1F("probXdelta", "probXdelta", 110, 0, 1.1);
1028   probXshared = new TH1F("probXshared", "probXshared", 110, 0, 1.1);
1029   probXnoshare = new TH1F("probXnoshare", "probXnoshare", 110, 0, 1.1);
1030   probYgood = new TH1F("probYgood", "probYgood", 110, 0, 1.1);
1031   probYbad = new TH1F("probYbad", "probYbad", 110, 0, 1.1);
1032   probYdelta = new TH1F("probYdelta", "probYdelta", 110, 0, 1.1);
1033   probYshared = new TH1F("probYshared", "probYshared", 110, 0, 1.1);
1034   probYnoshare = new TH1F("probYnoshare", "probYnoshare", 110, 0, 1.1);
1035   TH1::AddDirectory(oldAddDir);
1036 }
1037 // ------------ method called once each job just after ending the event loop  ------------
1038 void TestOutliers::endJob() {
1039   LogTrace("TestOutliers") << "TestOutliers::endJob";
1040   file->Write();
1041   LogTrace("TestOutliers") << "outfile written";
1042   file->Close();
1043   LogTrace("TestOutliers") << "oufile closed";
1044   LogTrace("TestOutliers") << "exiting TestOutliers::endJob";
1045 }
1046 
1047 //define this as a plug-in
1048 DEFINE_FWK_MODULE(TestOutliers);