ResidualRefitting.h

Go to the documentation of this file.

#ifndef _ResidualRefitting_h__
#define __ResidualRefitting_h_ (1)


#include <vector>
#include <string>

#include "TFile.h"
#include "TBranch.h"
#include "TTree.h"

#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/MuonReco/interface/Muon.h"
#include "DataFormats/MuonReco/interface/MuonFwd.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"

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

#include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
#include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"

#include "MagneticField/Engine/interface/MagneticField.h"
#include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
#include "TrackingTools/GeomPropagators/interface/Propagator.h"
//#include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
//#include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"

//#include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"

//#include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"

#include "TrackingTools/PatternTools/interface/Trajectory.h"
//#include "TrackingTools/TransientTrack/interface/TransientTrack.h"
//#include "TrackingTools/TrackRefitter/interface/TrackTransformer.h"

class TrackerTopology;

class ResidualRefitting : public edm::EDAnalyzer{

    static const int N_MAX_STORED = 10;
    static const int N_MAX_STORED_HIT = 1000;

    static const int PXB = 1;
    static const int PXF = 2;
    static const int TIB = 3;
    static const int TID = 4;
    static const int TOB = 5;
    static const int TEC = 6;
    

 public:

//  typedef std::pair<const Trajectory*, const reco::Track*> ConstTrajTrackPair;
//  typedef std::vector< ConstTrajTrackPair >  ConstTrajTrackPairCollection;

    typedef struct {
        int evtNum_;
        int runNum_;
    } storage_event;

    ResidualRefitting::storage_event eventInfo_;

 
    typedef struct StorageMuon {

        int n_;
    
        int     charge_     [N_MAX_STORED];
        float   pt_         [N_MAX_STORED];
        float   eta_        [N_MAX_STORED];
        float   p_          [N_MAX_STORED];
        float   phi_        [N_MAX_STORED];
        int     numRecHits_     [N_MAX_STORED];
        float   chiSq_      [N_MAX_STORED];
        float   ndf_        [N_MAX_STORED];
        float   chiSqOvrNdf_    [N_MAX_STORED];

            StorageMuon() :
          n_(0) {
          for (int i=0;i<N_MAX_STORED;++i) {
            charge_[i] = 0;
            pt_[i] = 0.;
            eta_[i] = 0.;
            p_[i] = 0.;
            phi_[i] = 0.;
            numRecHits_[i] = 0;
            chiSq_[i] = 0.;
            ndf_[i] = 0.;
            chiSqOvrNdf_[i] = 0.;
          }
        }
    } storage_muon; // Storage for standard muon information

    typedef struct StorageHit {
    
        int n_;
        int muonLink_           [N_MAX_STORED_HIT];
        
        int system_     [N_MAX_STORED_HIT];
        int endcap_     [N_MAX_STORED_HIT];
        int station_            [N_MAX_STORED_HIT];
        int ring_       [N_MAX_STORED_HIT];
        int chamber_            [N_MAX_STORED_HIT];
        int layer_      [N_MAX_STORED_HIT];
        int superLayer_         [N_MAX_STORED_HIT];
        int wheel_      [N_MAX_STORED_HIT];
        int sector_     [N_MAX_STORED_HIT];
        
        float   gpX_        [N_MAX_STORED_HIT];
        float   gpY_        [N_MAX_STORED_HIT];
        float   gpZ_        [N_MAX_STORED_HIT];
        float   gpEta_      [N_MAX_STORED_HIT];
        float   gpPhi_      [N_MAX_STORED_HIT];
        float   lpX_        [N_MAX_STORED_HIT];
        float   lpY_        [N_MAX_STORED_HIT];
        float   lpZ_        [N_MAX_STORED_HIT];
        
            StorageHit() :
          n_(0) {
          for (int i=0;i<N_MAX_STORED_HIT;++i) {
            muonLink_[i] = 0;
            system_[i] = 0;
            endcap_[i] = 0;
            station_[i] = 0;
            ring_[i] = 0;
            chamber_[i] = 0;
            layer_[i] = 0;
            superLayer_[i] = 0;
            wheel_[i] = 0;
            sector_[i] = 0;
            gpX_[i] = 0.;
            gpY_[i] = 0.;
            gpZ_[i] = 0.;
            gpEta_[i] = 0.;
            gpPhi_[i] = 0.;
            lpX_[i] = 0.;
            lpY_[i] = 0.;
            lpZ_[i] = 0.;           
          }
        }
    } storage_hit;

    typedef struct StorageTrackExtrap {
    
        int n_;
        
        int muonLink_       [N_MAX_STORED_HIT];
        int recLink_        [N_MAX_STORED_HIT];
        float gpX_          [N_MAX_STORED_HIT];
        float gpY_          [N_MAX_STORED_HIT];
        float gpZ_          [N_MAX_STORED_HIT];
        float gpEta_        [N_MAX_STORED_HIT];
        float gpPhi_        [N_MAX_STORED_HIT];
        float lpX_          [N_MAX_STORED_HIT];
        float lpY_          [N_MAX_STORED_HIT];
        float lpZ_          [N_MAX_STORED_HIT];
        float resX_         [N_MAX_STORED_HIT];
        float resY_         [N_MAX_STORED_HIT];
        float resZ_         [N_MAX_STORED_HIT];

            StorageTrackExtrap() :
          n_(0) {
          for (int i=0;i<N_MAX_STORED_HIT;++i) {
            muonLink_[i] = 0;
            recLink_[i] = 0;
            gpX_[i] = 0.;
            gpY_[i] = 0.;
            gpZ_[i] = 0.;
            gpEta_[i] = 0.;
            gpPhi_[i] = 0.;
            lpX_[i] = 0.;
            lpY_[i] = 0.;
            lpZ_[i] = 0.;
            resX_[i] = 0.;
            resY_[i] = 0.;
            resZ_[i] = 0.;
          }
        }
    } storage_trackExtrap;
        
    typedef struct StorageTrackHit {
    
        int n_;
        
        int muonLink_   [N_MAX_STORED_HIT];
        int detector_   [N_MAX_STORED_HIT];
        int subdetector_[N_MAX_STORED_HIT];
        int blade_      [N_MAX_STORED_HIT];
        int disk_       [N_MAX_STORED_HIT];
        int ladder_     [N_MAX_STORED_HIT];
        int layer_      [N_MAX_STORED_HIT];
        int module_     [N_MAX_STORED_HIT];
        int panel_      [N_MAX_STORED_HIT];
        int ring_       [N_MAX_STORED_HIT];
        int side_       [N_MAX_STORED_HIT];
        int wheel_      [N_MAX_STORED_HIT];
                
        float gpX_      [N_MAX_STORED_HIT];
        float gpY_      [N_MAX_STORED_HIT];
        float gpZ_      [N_MAX_STORED_HIT];
        float gpEta_    [N_MAX_STORED_HIT];
        float gpPhi_    [N_MAX_STORED_HIT];
        float lpX_      [N_MAX_STORED_HIT];
        float lpY_      [N_MAX_STORED_HIT];
        float lpZ_      [N_MAX_STORED_HIT];
                
            StorageTrackHit() :
          n_(0) {
          for (int i=0;i<N_MAX_STORED_HIT;++i) {
            muonLink_[i] = 0;
            detector_[i] = 0;
            subdetector_[i] = 0;
            blade_[i] = 0;
            disk_[i] = 0;
            ladder_[i] = 0;
            layer_[i] = 0;
            module_[i] = 0;
            panel_[i] = 0;
            ring_[i] = 0;
            side_[i] = 0;
            wheel_[i] = 0;
            gpX_[i] = 0.;
            gpY_[i] = 0.;
            gpZ_[i] = 0.;
            gpEta_[i] = 0.;
            gpPhi_[i] = 0.;
            lpX_[i] = 0.;
            lpY_[i] = 0.;
            lpZ_[i] = 0.;
          }
        }
    } storage_trackHit;
        
//Standard Muon info storage
    ResidualRefitting::storage_muon
     storageGmrOld_  ,  
     storageGmrNew_  , storageSamNew_   , storageTrkNew_,
     storageGmrNoSt1_, storageSamNoSt1_ ,
     storageGmrNoSt2_, storageSamNoSt2_ ,
     storageGmrNoSt3_, storageSamNoSt3_ ,
     storageGmrNoSt4_, storageSamNoSt4_ ,
     
     storageGmrNoPXBLayer1, storageGmrNoPXBLayer2, storageGmrNoPXBLayer3, 
     storageTrkNoPXBLayer1, storageTrkNoPXBLayer2, storageTrkNoPXBLayer3,       

     storageGmrNoPXF, storageTrkNoPXF,

     storageGmrNoTIBLayer1, storageGmrNoTIBLayer2, storageGmrNoTIBLayer3, storageGmrNoTIBLayer4, 
     storageTrkNoTIBLayer1, storageTrkNoTIBLayer2, storageTrkNoTIBLayer3, storageTrkNoTIBLayer4, 
     storageGmrNoTID, storageTrkNoTID,
     storageGmrNoTOBLayer1, storageGmrNoTOBLayer2, storageGmrNoTOBLayer3, storageGmrNoTOBLayer4, storageGmrNoTOBLayer5, storageGmrNoTOBLayer6,
     storageTrkNoTOBLayer1, storageTrkNoTOBLayer2, storageTrkNoTOBLayer3, storageTrkNoTOBLayer4, storageTrkNoTOBLayer5, storageTrkNoTOBLayer6,
     storageGmrNoTEC, storageTrkNoTEC;  
     
     
     
//Rec hit storage    
    ResidualRefitting::storage_hit          storageRecMuon_;
    ResidualRefitting::storage_trackHit storageTrackHit_;

//Track Extrapolation to Muon System
    ResidualRefitting::storage_trackExtrap      storageTrackExtrapRec_, storageTrackExtrapRecNoSt1_,
                                            storageTrackExtrapRecNoSt2_, storageTrackExtrapRecNoSt3_, storageTrackExtrapRecNoSt4_;

//Track Extrapolation with Cylinder
    ResidualRefitting::storage_trackExtrap      trackExtrap120_, samExtrap120_;
   
//Track Extrapolation to Tracker system
    ResidualRefitting::storage_trackExtrap      storageTrackExtrapTracker_,
            storageTrackNoPXBLayer1, storageTrackNoPXBLayer2, storageTrackNoPXBLayer3,
            storageTrackNoPXF,
            storageTrackNoTIBLayer1, storageTrackNoTIBLayer2, storageTrackNoTIBLayer3,storageTrackNoTIBLayer4,
            storageTrackNoTID,
            storageTrackNoTOBLayer1, storageTrackNoTOBLayer2, storageTrackNoTOBLayer3,
            storageTrackNoTOBLayer4, storageTrackNoTOBLayer5, storageTrackNoTOBLayer6,
            storageTrackNoTEC;
  
//
// Start of the method declarations
//
 
    explicit ResidualRefitting( const edm::ParameterSet & );
    ~ResidualRefitting() override;
    
    void analyze(const edm::Event&, const edm::EventSetup&) override;
    void beginJob() override ;
    void endJob() override ;
//Zero Storage
    void zero_storage();
    void zero_muon(ResidualRefitting::storage_muon* str);
    void zero_trackExtrap(ResidualRefitting::storage_trackExtrap* str);
    void branchMuon(ResidualRefitting::storage_muon& storageTmp, std::string branchName);
    void branchTrackExtrap(ResidualRefitting::storage_trackExtrap& storageTmp, std::string branchName);
    
    

//  void collectTrackRecExtrap(reco::MuonCollection::const_iterator muon, ResidualRefitting::storage_trackExtrap& storeTemp);
    void muonInfo(ResidualRefitting::storage_muon& storeMuon, reco::TrackRef muon, int val);
    
    void CollectTrackHits(edm::Handle<reco::TrackCollection> trackColl, ResidualRefitting::storage_trackExtrap& trackExtrap,
                              const edm::EventSetup& eventSetup);
    void StoreTrackerRecHits(DetId detid, const TrackerTopology* tTopo, int iTrack, int iRec);
    void NewTrackMeasurements(edm::Handle<reco::TrackCollection> trackCollOrig, edm::Handle<reco::TrackCollection> trackColl,
      ResidualRefitting::storage_trackExtrap& trackExtrap);
    int MatchTrackWithRecHits(reco::TrackCollection::const_iterator trackIt, edm::Handle<reco::TrackCollection> ref);
    
    bool IsSameHit(TrackingRecHit const& hit1, TrackingRecHit const& hit2);


    void trkExtrap(const DetId& detid, int iTrkLink, int iTrk, int iRec,
               const FreeTrajectoryState& freeTrajState,
               const LocalPoint& recPoint,
               storage_trackExtrap& storeTemp);
    
    void cylExtrapTrkSam(int recNum, reco::TrackRef track, ResidualRefitting::storage_trackExtrap& storage, double rho);

//Simplifiying functions
    FreeTrajectoryState freeTrajStateMuon(reco::TrackRef muon);//Returns a Free Trajectory State
//Debug Data Dumps
//  void dumpRecoMuonColl(reco::MuonCollection::const_iterator muon); //
//  void dumpRecoTrack(reco::TrackCollection::const_iterator muon);
    void dumpTrackRef(reco::TrackRef muon, std::string str); 
    void dumpTrackExtrap(const ResidualRefitting::storage_trackExtrap& track);
    void dumpTrackHits(const ResidualRefitting::storage_trackHit& hit);
    void dumpMuonRecHits(const ResidualRefitting::storage_hit& hit);

    int ReturnSector(DetId detid);
    int ReturnStation(DetId detid);

    
// Deprecated Functions
    void omitStation(edm::Handle<reco::MuonCollection> funcMuons, edm::Handle<reco::TrackCollection>,
        ResidualRefitting::storage_muon& storeGmr, ResidualRefitting::storage_muon& storeSam,
        ResidualRefitting::storage_trackExtrap& storeExtrap, int omitStation);
    void omitTrackerSystem(edm::Handle<reco::MuonCollection> trkMuons, ResidualRefitting::storage_muon& storeGmr,
        ResidualRefitting::storage_muon& storeTrk,
        ResidualRefitting::storage_trackExtrap& storeExtrap, int omitSystem);
      
    // output histogram file name
    std::string outputFileName_;
    //edm::InputTag PropagatorSource_;      
    std::string PropagatorSource_;      

    // names of product labels
    edm::InputTag tracks_,
     muons_, muonsRemake_,
     muonsNoStation1_, muonsNoStation2_, muonsNoStation3_,muonsNoStation4_,//Global Muon Collections
     muonsNoPXBLayer1_, muonsNoPXBLayer2_, muonsNoPXBLayer3_,
     muonsNoPXF_,
     muonsNoTIBLayer1_, muonsNoTIBLayer2_, muonsNoTIBLayer3_, muonsNoTIBLayer4_,
     muonsNoTID_,
     muonsNoTOBLayer1_, muonsNoTOBLayer2_, muonsNoTOBLayer3_, muonsNoTOBLayer4_, muonsNoTOBLayer5_, muonsNoTOBLayer6_,
     muonsNoTEC_;
//     tjTag; 


    bool debug_;

    // output ROOT file
    TFile * outputFile_ ;

    TTree   *outputTree_;
    TBranch *outputBranch_;

//  unsigned int nBins_;
  
    const MagneticField* theField;
    const edm::ESHandle<GlobalTrackingGeometry> trackingGeometry;
    MuonServiceProxy* theService;
    edm::ESHandle<Propagator> thePropagator;
  
};

#endif