---

VisMuon.cc

Go to the documentation of this file.

//<<<<<< INCLUDES                                                       >>>>>>

#include "VisReco/Analyzer/interface/VisMuon.h"
#include "VisReco/Analyzer/interface/IguanaService.h"
#include "DataFormats/GeometrySurface/interface/PlaneBuilder.h"
#include "DataFormats/MuonReco/interface/Muon.h"
#include "DataFormats/MuonReco/interface/MuonFwd.h"
#include "DataFormats/MuonReco/interface/MuonChamberMatch.h"
#include "DataFormats/MuonReco/interface/MuonSelectors.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackBase.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "FWCore/Utilities/interface/Exception.h"
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
#include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
#include "TrackingTools/TransientTrack/interface/TransientTrack.h"
#include "Iguana/Framework/interface/IgCollection.h"
#include <iostream>
#include <sstream>

//<<<<<< PRIVATE DEFINES                                                >>>>>>
//<<<<<< PRIVATE CONSTANTS                                              >>>>>>
//<<<<<< PRIVATE TYPES                                                  >>>>>>
//<<<<<< PRIVATE VARIABLE DEFINITIONS                                   >>>>>>
//<<<<<< PUBLIC VARIABLE DEFINITIONS                                    >>>>>>
//<<<<<< CLASS STRUCTURE INITIALIZATION                                 >>>>>>
//<<<<<< PRIVATE FUNCTION DEFINITIONS                                   >>>>>>
//<<<<<< PUBLIC FUNCTION DEFINITIONS                                    >>>>>>
//<<<<<< MEMBER FUNCTION DEFINITIONS                                    >>>>>>

using namespace edm::service;

VisMuon::VisMuon (const edm::ParameterSet& iConfig)
    : inputTag_ (iConfig.getParameter<edm::InputTag>("visMuonTag")),
      in_ (iConfig.getUntrackedParameter<double>( "propagatorIn", 0.0)),
      out_ (iConfig.getUntrackedParameter<double>( "propagatorOut", 0.0)),
      step_ (iConfig.getUntrackedParameter<double>( "propagatorStep", 0.05))
{}

void 
VisMuon::analyze (const edm::Event& event, const edm::EventSetup& eventSetup)
{
    edm::Service<IguanaService> config;
    if (! config.isAvailable ()) 
    {
        throw cms::Exception ("Configuration")
            << "VisMuon requires the IguanaService\n"
            "which is not present in the configuration file.\n"
            "You must add the service in the configuration file\n"
            "or remove the module that requires it";
    }
    
    edm::Handle<reco::MuonCollection> collection;

    eventSetup.get<IdealMagneticFieldRecord> ().get (field_);
    event.getByLabel (inputTag_, collection);

    if (collection.isValid ())
    {       
        storage_ = config->storage ();
        IgCollection &icollection = storage_->getCollection ("Muons_V1");
        IgProperty PT = icollection.addProperty ("pt", 0.0);
        IgProperty CHARGE = icollection.addProperty ("charge", 0.0);
        IgProperty RP = icollection.addProperty ("rp", IgV3d ());
        IgProperty PHI = icollection.addProperty ("phi", 0.0);
        IgProperty ETA = icollection.addProperty ("eta", 0.0);
        IgProperty T_PT = icollection.addProperty ("tracker_pt", 0.0);
        IgProperty S_PT = icollection.addProperty ("standalone_pt", 0.0);
        IgProperty G_PT = icollection.addProperty ("global_pt", 0.0);
        IgProperty CALO_E = icollection.addProperty ("calo_energy", 0.0);

        IgCollection &points = storage_->getCollection("Points_V1");
        IgProperty POS = points.addProperty("pos", IgV3d());

        IgCollection &detIds = storage_->getCollection("DetIds_V1");
        IgProperty DETID = detIds.addProperty ("detid", int (0));

        IgAssociationSet &muonDetIds = storage_->getAssociationSet("MuonDetIds_V1");

        reco::MuonCollection::const_iterator it = collection->begin ();
        reco::MuonCollection::const_iterator end = collection->end ();
        for (; it != end; ++it) 
        {
            IgCollectionItem imuon = icollection.create ();
            float charge = (*it).charge ();
            float pt = (*it).pt();
            imuon [PT] = static_cast<double>(pt);
            imuon [CHARGE] = static_cast<double>(charge);
            imuon [RP] = IgV3d (static_cast<double>((*it).vx ()/100.), static_cast<double>((*it).vy()/100.), static_cast<double>((*it).vz()/100.));
            imuon [PHI] = static_cast<double>((*it).phi ());
            imuon [ETA] = static_cast<double>((*it).eta ());
            
            if ((*it).isMatchesValid ())
            {               
                const std::vector<reco::MuonChamberMatch> &dets = (*it).matches ();
                for (std::vector<reco::MuonChamberMatch>::const_iterator dit = dets.begin (), ditEnd = dets.end (); dit != ditEnd; ++dit)
                {                   
                    IgCollectionItem idetId = detIds.create ();
                    idetId[DETID] = static_cast<int>((*dit).id);
                    muonDetIds.associate (imuon, idetId);
                }
            }
            if ((*it).innerTrack ().isNonnull ()) // Tracker
            {
                double pt = (*it).innerTrack ()->pt ();
                imuon[T_PT] = static_cast<double>(pt);
                IgAssociationSet &muonTrackerPoints = storage_->getAssociationSet("MuonTrackerPoints_V1");      
                refitTrack (imuon, muonTrackerPoints, (*it).innerTrack (), in_, out_, step_);           
            }
                
            if ((*it).outerTrack ().isNonnull ()) // Standalone
            {
                double pt = (*it).outerTrack ()->pt ();             
                imuon[S_PT] = static_cast<double>(pt);
                IgAssociationSet &muonStandalonePoints = storage_->getAssociationSet("MuonStandalonePoints_V1");        
                refitTrack (imuon, muonStandalonePoints, (*it).outerTrack (), in_, out_, step_);                
            }
                
            if ((*it).globalTrack ().isNonnull ()) // Global
            {
                double pt = (*it).globalTrack ()->pt ();
                imuon[G_PT] = static_cast<double>(pt);
                IgAssociationSet &muonGlobalPoints = storage_->getAssociationSet("MuonGlobalPoints_V1");        
                refitTrack (imuon, muonGlobalPoints, (*it).globalTrack (), in_, out_, step_);           
            }
                
            if ((*it).isEnergyValid ()) // CaloTower
            {
                float e = (*it).calEnergy ().tower;
                imuon[CALO_E] = static_cast<double>(e);         
            }
        }
    }
    else 
    {
        // friendlyName:moduleLabel:instanceName:processName
        std::string error = "### Error: Muons "
                            + edm::TypeID (typeid (reco::MuonCollection)).friendlyClassName () + ":" 
                            + inputTag_.label() + ":"
                            + inputTag_.instance() + ":" 
                            + inputTag_.process() + " are not found.";

        IgDataStorage *storage = config->storage ();
        IgCollection &collection = storage->getCollection ("Errors_V1");
        IgProperty errorMsg = collection.addProperty ("Error", std::string ());
        IgCollectionItem item = collection.create ();
        item ["Error"] = error;
    }
}

void
VisMuon::refitTrack (IgCollectionItem &imuon, IgAssociationSet &association, reco::TrackRef track, double in, double out, double step) 
{
    if (track.isNonnull () && field_.isValid ())
    {
        IgCollection &points = storage_->getCollection("Points_V1");

        reco::TransientTrack aRealTrack (track, &*field_);

        //Get a propagator and inner/outer state
        SteppingHelixPropagator propagator (&*field_, alongMomentum);
        SteppingHelixPropagator reversePropagator (&*field_, oppositeToMomentum);
        GlobalPoint gPin ((*track).innerPosition ().x (),
                          (*track).innerPosition ().y (),
                          (*track).innerPosition ().z ());
        GlobalPoint gPout ((*track).outerPosition ().x (),
                           (*track).outerPosition ().y (),
                           (*track).outerPosition ().z ());
        GlobalVector InOutVector = (gPout - gPin);

        // Define rotation for plane perpendicular to InOutVector
        // z axis coincides with perp
        GlobalVector zAxis = InOutVector.unit();
        // x axis has no global Z component
        GlobalVector xAxis;
        if (zAxis.x() != 0 || zAxis.y() != 0) 
        {
            // precision is not an issue here, just protect against divizion by zero
            xAxis = GlobalVector (-zAxis.y(), zAxis.x(), 0).unit();
        }
        else { // perp coincides with global Z
            xAxis = GlobalVector( 1, 0, 0);
        }
        // y axis obtained by cross product
        GlobalVector yAxis( zAxis.cross( xAxis));
        Surface::RotationType rot( xAxis.x(), xAxis.y(), xAxis.z(),
                                   yAxis.x(), yAxis.y(), yAxis.z(),
                                   zAxis.x(), zAxis.y(), zAxis.z());

        // Define step size and number of extra steps on inside and outside
        int nSteps = 0;
        int nStepsInside = 0;
        int nStepsOutside = 0;
        GlobalVector StepVector;
        if( step > 0.01 ) {
            StepVector = InOutVector * step;
            nSteps = int (0.5 + 1.0/step);
            nStepsInside = int (0.5 + in/step);
            nStepsOutside = int (0.5 + out/step);
        } else {
            StepVector = InOutVector * 0.01;
            nSteps = 100;
            nStepsInside = int (0.5 + in*100);
            nStepsOutside = int (0.5 + out*100);          
        } 

        GlobalVector gVin ((*track).innerMomentum ().x (),
                           (*track).innerMomentum ().y (),
                           (*track).innerMomentum ().z ());
        GlobalVector gVout ((*track).outerMomentum ().x (),
                            (*track).outerMomentum ().y (),
                            (*track).outerMomentum ().z ());
        GlobalTrajectoryParameters GTPin (gPin, gVin, aRealTrack.impactPointTSCP ().charge (), &*field_);
        FreeTrajectoryState FTSin (GTPin);
        GlobalTrajectoryParameters GTPout (gPout, gVout, aRealTrack.impactPointTSCP ().charge (), &*field_);
        FreeTrajectoryState FTSout (GTPout);

        int nGood = 0;
        int nBad = 0;

        // Do nStepsInside propagations on inside to plot track
        GlobalPoint GP = gPin; 
        GP -= nStepsInside * StepVector;
        Basic3DVector<float> Basic3DV (StepVector.x (), StepVector.y (), StepVector.z ());
        for (int istep = 0; istep < nStepsInside; ++istep) 
        {
            GP += StepVector;
            Surface::PositionType pos (GP.x (), GP.y (), GP.z ());
            PlaneBuilder::ReturnType SteppingPlane = PlaneBuilder ().plane (pos, rot);
            TrajectoryStateOnSurface trj =  reversePropagator.propagate (FTSin, *SteppingPlane);
            if (trj.isValid ())
            {
                float x = trj.globalPosition ().x () / 100.0;
                float y = trj.globalPosition ().y () / 100.0;
                float z = trj.globalPosition ().z () / 100.0;              
                IgCollectionItem ipoint = points.create ();
                ipoint["pos"] = IgV3d (static_cast<double>(x), static_cast<double>(y), static_cast<double>(z));
                association.associate (imuon, ipoint);          
                nGood++;
            } 
            else 
            {
                nBad++;
            }
        }
        // Do nStep propagations from track Inner state to track Outer state

        GP = gPin - StepVector;
        float w = 0;
        for (int istep = 0; istep < nSteps+1 ; ++istep) 
        { // from innerPosition to outerPosition
            GP += StepVector;
            Surface::PositionType pos (GP.x (), GP.y (), GP.z ());
            PlaneBuilder::ReturnType SteppingPlane = PlaneBuilder ().plane (pos, rot);
            TrajectoryStateOnSurface trj_in =  propagator.propagate (FTSin, *SteppingPlane);
            TrajectoryStateOnSurface trj_out =  reversePropagator.propagate (FTSout, *SteppingPlane);
            if (trj_in.isValid () && trj_out.isValid ()) 
            {
                float x1 = trj_in.globalPosition ().x () / 100.0;
                float y1 = trj_in.globalPosition ().y () / 100.0;
                float z1 = trj_in.globalPosition ().z () / 100.0;                  

                float x2 = trj_out.globalPosition ().x () / 100.0;
                float y2 = trj_out.globalPosition ().y () / 100.0;
                float z2 = trj_out.globalPosition ().z () / 100.0;                 

                float ww = w<0.4999 ? ww = w : ww=1.0-w;
                float w2 = 0.5*ww*ww/((1.0-ww)*(1.0-ww));
                if(w>0.4999) w2=1.0-w2;  

                float x = (1.0-w2)*x1 + w2*x2;
                float y = (1.0-w2)*y1 + w2*y2;
                float z = (1.0-w2)*z1 + w2*z2;

                IgCollectionItem ipoint = points.create ();
                ipoint["pos"] = IgV3d (static_cast<double>(x), static_cast<double>(y), static_cast<double>(z));
                association.associate (imuon, ipoint);          
                nGood++;
            } 
            else 
            {
                nBad++;
            }
            w += 1.0/float(nSteps);
        }
        // Do nStepsInside propagations on Outside to plot track
        GP = gPout;
        for (int istep = 0; istep < nStepsOutside; ++istep) 
        {
            GP += StepVector;
            Surface::PositionType pos (GP.x (), GP.y (), GP.z ());
            PlaneBuilder::ReturnType SteppingPlane = PlaneBuilder ().plane (pos, rot);
            TrajectoryStateOnSurface trj =  propagator.propagate (FTSout, *SteppingPlane);
            if (trj.isValid ())
            {
                float x = trj.globalPosition ().x () / 100.0;
                float y = trj.globalPosition ().y () / 100.0;
                float z = trj.globalPosition ().z () / 100.0;              
                IgCollectionItem ipoint = points.create ();
                ipoint["pos"] = IgV3d (static_cast<double>(x), static_cast<double>(y), static_cast<double>(z));
                association.associate (imuon, ipoint);          
                nGood++;
            } 
            else 
            {
                nBad++;
            }
        }
    }
}

DEFINE_FWK_MODULE(VisMuon);

---