Reference Manual

Go to the documentation of this file.
00001 /*
00002  * Package:     CommonAlignmentProducer
00003  * Class  :     AlignmentMonitorMuonVsCurvature
00004  *
00005  * Original Author:  Jim Pivarski
00006  *         Created:  Fri Feb 19 21:45:02 CET 2010
00007  *
00008  * $Id: AlignmentMonitorMuonVsCurvature.cc,v 1.6 2011/10/12 22:59:47 khotilov Exp $
00009  */
00010 
00011 #include "Alignment/CommonAlignmentMonitor/interface/AlignmentMonitorPluginFactory.h"
00012 #include "Alignment/CommonAlignmentMonitor/interface/AlignmentMonitorBase.h"
00013 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsFromTrack.h"
00014 
00015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00016 #include "FWCore/Utilities/interface/InputTag.h"
00017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00018 #include "FWCore/Framework/interface/EventSetup.h"
00019 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
00020 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
00021 
00022 #include <sstream>
00023 #include "TProfile.h"
00024 #include "TH2F.h"
00025 #include "TH1F.h"
00026 
00027 
00028 class AlignmentMonitorMuonVsCurvature: public AlignmentMonitorBase
00029 {
00030 public:
00031   AlignmentMonitorMuonVsCurvature(const edm::ParameterSet& cfg);
00032   virtual ~AlignmentMonitorMuonVsCurvature() {}
00033 
00034   void book();
00035 
00036   void event(const edm::Event &iEvent, const edm::EventSetup &iSetup, const ConstTrajTrackPairCollection& iTrajTracks);
00037   void processMuonResidualsFromTrack(MuonResidualsFromTrack &mrft, const Trajectory* traj = NULL);
00038 
00039   void afterAlignment(const edm::EventSetup &iSetup) {}
00040 
00041 private:
00042   
00043   edm::InputTag m_muonCollectionTag;
00044   double m_minTrackPt;
00045   double m_minTrackP;
00046   int m_minTrackerHits;
00047   double m_maxTrackerRedChi2;
00048   bool m_allowTIDTEC;
00049   bool m_minNCrossedChambers;
00050   double m_maxDxy;
00051   int m_minDT13Hits;
00052   int m_minDT2Hits;
00053   int m_minCSCHits;
00054   int m_layer;
00055   std::string m_propagator;
00056   bool m_doDT;
00057   bool m_doCSC;
00058 
00059   enum {
00060     kDeltaX = 0,
00061     kDeltaDxDz,
00062     kNumComponents
00063   };
00064 
00065   // DT [wheel][station][sector]
00066   TH2F *th2f_wheel_st_sector[5][4][14][kNumComponents];
00067   TProfile *tprofile_wheel_st_sector[5][4][14][kNumComponents];
00068   
00069   // CSC [endcap*station][ring][chamber]
00070   TH2F *th2f_st_ring_chamber[8][3][36][kNumComponents];
00071   TProfile *tprofile_st_ring_chamber[8][3][36][kNumComponents];
00072 
00073   TH1F *th1f_trackerRedChi2;
00074   TH1F *th1f_trackerRedChi2Diff;
00075 };
00076 
00077 
00078 AlignmentMonitorMuonVsCurvature::AlignmentMonitorMuonVsCurvature(const edm::ParameterSet& cfg)
00079 : AlignmentMonitorBase(cfg, "AlignmentMonitorMuonVsCurvature")
00080 , m_muonCollectionTag(cfg.getParameter<edm::InputTag>("muonCollectionTag"))
00081 , m_minTrackPt(cfg.getParameter<double>("minTrackPt"))
00082 , m_minTrackP(cfg.getParameter<double>("minTrackP"))
00083 , m_minTrackerHits(cfg.getParameter<int>("minTrackerHits"))
00084 , m_maxTrackerRedChi2(cfg.getParameter<double>("maxTrackerRedChi2"))
00085 , m_allowTIDTEC(cfg.getParameter<bool>("allowTIDTEC"))
00086 , m_minNCrossedChambers(cfg.getParameter<int>("minNCrossedChambers"))
00087 , m_maxDxy(cfg.getParameter<double>("maxDxy"))
00088 , m_minDT13Hits(cfg.getParameter<int>("minDT13Hits"))
00089 , m_minDT2Hits(cfg.getParameter<int>("minDT2Hits"))
00090 , m_minCSCHits(cfg.getParameter<int>("minCSCHits"))
00091 , m_layer(cfg.getParameter<int>("layer"))
00092 , m_propagator(cfg.getParameter<std::string>("propagator"))
00093 , m_doDT(cfg.getParameter<bool>("doDT"))
00094 , m_doCSC(cfg.getParameter<bool>("doCSC"))
00095 {}
00096 
00097 
00098 void AlignmentMonitorMuonVsCurvature::book()
00099 {
00100   // DT
00101   std::string wheelname[5] = {"wheelm2_", "wheelm1_", "wheelz_", "wheelp1_", "wheelp2_"};
00102   if (m_doDT)
00103   for (int wheel = -2;  wheel <=2 ;  wheel++)
00104   for (int station = 1; station <= 4; station++)
00105   for (int sector = 1;  sector <= 14;  sector++)
00106   {
00107     if (station != 4 && sector > 12) continue;
00108     
00109     char stationname[20];
00110     sprintf(stationname,"st%d_", station);
00111 
00112     char sectorname[20];
00113     sprintf(sectorname,"sector%02d_", sector);
00114 
00115     for (int component = 0;  component < kNumComponents;  component++)
00116     {
00117       std::stringstream th2f_name, tprofile_name;
00118       th2f_name << "th2f_" << wheelname[wheel+2] <<stationname<< sectorname;
00119       tprofile_name << "tprofile_" << wheelname[wheel+2] <<stationname<< sectorname;
00120 
00121       double yminmax = 50., xminmax = 0.05;
00122       if (m_minTrackPt>0.) xminmax = 1./m_minTrackPt;
00123       int ynbins = 50;
00124       if (component == kDeltaX) {
00125         th2f_name << "deltax";
00126         tprofile_name << "deltax";
00127       }
00128       else if (component == kDeltaDxDz) {
00129         th2f_name << "deltadxdz";
00130         tprofile_name << "deltadxdz";
00131       }
00132 
00133       th2f_wheel_st_sector[wheel+2][station-1][sector-1][component] =
00134           book2D("/iterN/", th2f_name.str().c_str(), "", 30, -xminmax, xminmax, ynbins, -yminmax, yminmax);
00135       tprofile_wheel_st_sector[wheel+2][station-1][sector-1][component] =
00136           bookProfile("/iterN/", tprofile_name.str().c_str(), "", 30,  -xminmax, xminmax);
00137     }
00138   }
00139 
00140   // CSC
00141   std::string stname[8] = {"Ep_S1_", "Ep_S2_", "Ep_S3_", "Ep_S4_", "Em_S1_", "Em_S2_", "Em_S3_", "Em_S4_"};
00142   if (m_doCSC)
00143   for (int station = 0;  station < 8;  station++)
00144   for (int ring = 1;  ring <= 3;  ring++)
00145   for (int chamber = 1;  chamber <= 36;  chamber++)
00146   {
00147     int st = station%4+1;
00148     if (st > 1 && ring > 2) continue; // only station 1 has more then 2 rings
00149     if (st > 1 && ring == 1 && chamber > 18) continue; // ring 1 stations 1,2,3 have 18 chambers
00150 
00151     char ringname[20];
00152     sprintf(ringname,"R%d_", ring);
00153 
00154     char chname[20];
00155     sprintf(chname,"C%02d_", chamber);
00156 
00157     for (int component = 0;  component < kNumComponents;  component++)
00158     {
00159       std::stringstream componentname;
00160       double yminmax = 50., xminmax = 0.05;
00161       if (m_minTrackPt>0.) xminmax = 1./m_minTrackPt;
00162       if (ring == 1) xminmax *= 0.5;
00163       if (component == kDeltaX) {
00164         componentname << "deltax";
00165       }
00166       else if (component == kDeltaDxDz) {
00167         componentname << "deltadxdz";
00168       }
00169       
00170       std::stringstream th2f_name, tprofile_name;
00171       th2f_name << "th2f_" << stname[station] << ringname << chname << componentname.str();
00172       tprofile_name << "tprofile_" << stname[station] << ringname << chname << componentname.str();
00173 
00174       th2f_st_ring_chamber[station][ring-1][chamber-1][component] =
00175           book2D("/iterN/", th2f_name.str().c_str(), "", 30, -xminmax, xminmax, 100, -yminmax, yminmax);
00176       tprofile_st_ring_chamber[station][ring-1][chamber-1][component] =
00177           bookProfile("/iterN/", tprofile_name.str().c_str(), "", 30, -xminmax, xminmax);
00178     }
00179   }
00180 
00181   th1f_trackerRedChi2 = book1D("/iterN/", "trackerRedChi2", "Refit tracker reduced chi^2", 100, 0., 30.);
00182   th1f_trackerRedChi2Diff = book1D("/iterN/", "trackerRedChi2Diff", "Fit-minus-refit tracker reduced chi^2", 100, -5., 5.);
00183 }
00184 
00185 
00186 void AlignmentMonitorMuonVsCurvature::event(const edm::Event &iEvent, const edm::EventSetup &iSetup, const ConstTrajTrackPairCollection& trajtracks)
00187 {
00188   edm::ESHandle<GlobalTrackingGeometry> globalGeometry;
00189   iSetup.get<GlobalTrackingGeometryRecord>().get(globalGeometry);
00190 
00191   edm::Handle<reco::BeamSpot> beamSpot;
00192   iEvent.getByLabel(m_beamSpotTag, beamSpot);
00193 
00194   if (m_muonCollectionTag.label().empty()) // use trajectories
00195   {
00196     for (ConstTrajTrackPairCollection::const_iterator trajtrack = trajtracks.begin();  trajtrack != trajtracks.end();  ++trajtrack)
00197     {
00198       const Trajectory* traj = (*trajtrack).first;
00199       const reco::Track* track = (*trajtrack).second;
00200 
00201       if (track->pt() > m_minTrackPt  && track->p() > m_minTrackP  &&  fabs(track->dxy(beamSpot->position())) < m_maxDxy )
00202       {
00203         MuonResidualsFromTrack muonResidualsFromTrack(globalGeometry, traj, track, pNavigator(), 1000.);
00204         processMuonResidualsFromTrack(muonResidualsFromTrack, traj );
00205       } // end if track pT is within range
00206     } // end loop over tracks
00207   }
00208   else
00209   {
00210     edm::Handle<reco::MuonCollection> muons;
00211     iEvent.getByLabel(m_muonCollectionTag, muons);
00212 
00213     for (reco::MuonCollection::const_iterator muon = muons->begin();  muon != muons->end();  ++muon)
00214     {
00215       if ( !(muon->isTrackerMuon() && muon->innerTrack().isNonnull() ) ) continue;
00216 
00217       if (m_minTrackPt < muon->pt()  &&  m_minTrackP < muon->p()  &&  fabs(muon->innerTrack()->dxy(beamSpot->position())) < m_maxDxy)
00218       {
00219         MuonResidualsFromTrack muonResidualsFromTrack(globalGeometry, &(*muon), pNavigator(), 100.);
00220         processMuonResidualsFromTrack(muonResidualsFromTrack);
00221       }
00222     }
00223   }
00224 }
00225 
00226 
00227 void AlignmentMonitorMuonVsCurvature::processMuonResidualsFromTrack(MuonResidualsFromTrack &mrft, const Trajectory* traj)
00228 {
00229   if (mrft.trackerNumHits() < m_minTrackerHits) return;
00230   if (!m_allowTIDTEC  && mrft.contains_TIDTEC()) return;
00231   
00232   int nMuChambers = 0;
00233   std::vector<DetId> chamberIds = mrft.chamberIds();
00234   for (unsigned ch=0; ch < chamberIds.size(); ch++)  if (chamberIds[ch].det() == DetId::Muon)  nMuChambers++;
00235   if (nMuChambers < m_minNCrossedChambers ) return;
00236   
00237   th1f_trackerRedChi2->Fill(mrft.trackerRedChi2());
00238   th1f_trackerRedChi2Diff->Fill(mrft.getTrack()->normalizedChi2() - mrft.trackerRedChi2());
00239 
00240   if (mrft.normalizedChi2() > m_maxTrackerRedChi2) return;
00241 
00242   double qoverpt = mrft.getTrack()->charge() / mrft.getTrack()->pt();
00243   double qoverpz = 0.;
00244   if (fabs(mrft.getTrack()->pz()) > 0.01) qoverpz = mrft.getTrack()->charge() / fabs(mrft.getTrack()->pz());
00245 
00246   for (std::vector<DetId>::const_iterator chamberId = chamberIds.begin();  chamberId != chamberIds.end();  ++chamberId)
00247   {
00248     if (chamberId->det() != DetId::Muon  ) continue;
00249     
00250     if (m_doDT  &&  chamberId->subdetId() == MuonSubdetId::DT)
00251     {
00252       DTChamberId dtid(chamberId->rawId());
00253       MuonChamberResidual *dt13 = mrft.chamberResidual(*chamberId, MuonChamberResidual::kDT13);
00254       
00255       if (dt13 != NULL  &&  dt13->numHits() >= m_minDT13Hits)
00256       {
00257         int wheel = dtid.wheel() + 2;
00258         int station = dtid.station() -1;
00259         int sector = dtid.sector() - 1;
00260 
00261         double resid_x = 10. * dt13->global_residual();
00262         double resid_dxdz = 1000. * dt13->global_resslope();
00263 
00264         if (fabs(resid_x) < 100. && fabs(resid_dxdz) < 100.)
00265         {
00266           th2f_wheel_st_sector[wheel][station][sector][kDeltaX]->Fill(qoverpt, resid_x);
00267           tprofile_wheel_st_sector[wheel][station][sector][kDeltaX]->Fill(qoverpt, resid_x);
00268           th2f_wheel_st_sector[wheel][station][sector][kDeltaDxDz]->Fill(qoverpt, resid_dxdz);
00269           tprofile_wheel_st_sector[wheel][station][sector][kDeltaDxDz]->Fill(qoverpt, resid_dxdz);
00270         }
00271       } // if it's a good segment
00272     } // if DT
00273 
00274     if (m_doCSC  &&  chamberId->subdetId() == MuonSubdetId::CSC)
00275     {
00276       CSCDetId cscid(chamberId->rawId());
00277       MuonChamberResidual *csc = mrft.chamberResidual(*chamberId, MuonChamberResidual::kCSC);
00278 
00279       if (csc != NULL  &&  csc->numHits() >= m_minCSCHits)
00280       {
00281         int station = 4*cscid.endcap() + cscid.station() - 5;
00282         int ring = cscid.ring() - 1;
00283         if (cscid.station()==1 && cscid.ring()==4) ring = 0; // join ME1/a to ME1/b
00284         int chamber = cscid.chamber() - 1;
00285 
00286         double resid_x = 10. * csc->global_residual();
00287         double resid_dxdz = 1000. * csc->global_resslope();
00288 
00289         if (fabs(resid_x) < 100. && fabs(resid_dxdz) < 100.)
00290         {
00291           th2f_st_ring_chamber[station][ring][chamber][kDeltaX]->Fill(qoverpz, resid_x);
00292           tprofile_st_ring_chamber[station][ring][chamber][kDeltaX]->Fill(qoverpz, resid_x);
00293           th2f_st_ring_chamber[station][ring][chamber][kDeltaDxDz]->Fill(qoverpz, resid_dxdz);
00294           tprofile_st_ring_chamber[station][ring][chamber][kDeltaDxDz]->Fill(qoverpz, resid_dxdz);
00295         }
00296       } // if it's a good segment
00297     } // if CSC
00298 
00299   } // end loop over chamberIds
00300 }
00301 
00302 DEFINE_EDM_PLUGIN(AlignmentMonitorPluginFactory, AlignmentMonitorMuonVsCurvature, "AlignmentMonitorMuonVsCurvature");
/data/refman/pasoursint/CMSSW_5_3_3/src/Alignment/CommonAlignmentMonitor/plugins/AlignmentMonitorMuonVsCurvature.cc
