| 1 |
// -*- C++ -*-
|
| 2 |
//
|
| 3 |
// Package: TPTrackMuonSys
|
| 4 |
// Class: TPTrackMuonSys
|
| 5 |
//
|
| 6 |
/**\class TPTrackMuonSys TPTrackMuonSys.cc RPhysicsStudies/TPTrackMuonSys/src/TPTrackMuonSys.cc
|
| 7 |
|
| 8 |
Description: [one line class summary]
|
| 9 |
|
| 10 |
Implementation:
|
| 11 |
[Notes on implementation]
|
| 12 |
*/
|
| 13 |
//
|
| 14 |
// Original Author: Chaouki Boulahouache,32 4-A05,+41227674763;
|
| 15 |
// Jinzhong Zhang,32 4-C21,+41227671337,
|
| 16 |
// Created: Wed May 12 15:52:13 CEST 2010
|
| 17 |
// $Id$
|
| 18 |
//
|
| 19 |
// system include files
|
| 20 |
#include <memory>
|
| 21 |
|
| 22 |
// user include files
|
| 23 |
#include "../interface/TPTrackMuonSys.h"
|
| 24 |
|
| 25 |
TPTrackMuonSys::TPTrackMuonSys(const edm::ParameterSet& Conf) : theDbConditions( Conf ) {
|
| 26 |
//now do what ever initialization is needed
|
| 27 |
m_rootFileName = Conf.getUntrackedParameter<std::string>("rootFileName","PhyTree.root");
|
| 28 |
///////////////////////////////
|
| 29 |
// Various input parameters.
|
| 30 |
//////////////////////////////
|
| 31 |
|
| 32 |
// A debug flag, here I can set different debugging steps based on the m_debug value,
|
| 33 |
m_gTracksTag = Conf.getUntrackedParameter<edm::InputTag>("gTracksTag"); // "generalTracks"
|
| 34 |
m_dEdxDiscrimTag = Conf.getUntrackedParameter<edm::InputTag>("dedxTag"); //dedxHarmonic2
|
| 35 |
trackProducer = Conf.getParameter<edm::InputTag>("trackProducer");
|
| 36 |
m_vertexSrc = Conf.getParameter<edm::InputTag>("vertexCollection");
|
| 37 |
|
| 38 |
m_hlt = Conf.getUntrackedParameter<edm::InputTag>("hltTag"); //TriggerResults::HLT
|
| 39 |
m_hltTrgEv = Conf.getUntrackedParameter<edm::InputTag>("hltEvTag"); //TriggerResults::HLT
|
| 40 |
m_L1extraTag = Conf.getUntrackedParameter<edm::InputTag>("L1extraTag");
|
| 41 |
|
| 42 |
m_HepMCTag = Conf.getUntrackedParameter<string>("HepMCTag","generator");
|
| 43 |
|
| 44 |
m_HLTMuTrgNames = Conf.getParameter< vector<string> >("HLTMuTrgNames");
|
| 45 |
m_HLTDiMuTrgName = Conf.getParameter< string >("HLTDiMuTrgName");
|
| 46 |
|
| 47 |
// TrackAssociator parameters
|
| 48 |
edm::ParameterSet parameters = Conf.getParameter<edm::ParameterSet>("TrackAssociatorParameters");
|
| 49 |
trackExtractorPSet_ = Conf.getParameter<edm::ParameterSet>("TrackExtractor");
|
| 50 |
parameters_.loadParameters( parameters );
|
| 51 |
trackAssociator_.useDefaultPropagator();
|
| 52 |
|
| 53 |
/// Now the MC specific information... if available...
|
| 54 |
//
|
| 55 |
m_isMC = Conf.getUntrackedParameter<bool>("isMC");
|
| 56 |
m_mcTag = Conf.getUntrackedParameter<edm::InputTag>("mcTag"); //
|
| 57 |
|
| 58 |
// flags to switch on/off individual modules
|
| 59 |
// set counters to zero
|
| 60 |
nEventsAnalyzed = 0;
|
| 61 |
treeCount = 0;
|
| 62 |
|
| 63 |
// Create the root file for the histograms
|
| 64 |
theFile = new TFile(m_rootFileName.c_str(), "RECREATE");
|
| 65 |
|
| 66 |
|
| 67 |
//////////////////// for LUTs
|
| 68 |
/*
|
| 69 |
bzero(srLUTs_, sizeof(srLUTs_));
|
| 70 |
//Int_t endcap =1, sector =1;
|
| 71 |
Bool_t TMB07=true;
|
| 72 |
edm::ParameterSet srLUTset;
|
| 73 |
srLUTset.addUntrackedParameter<bool>("ReadLUTs", false);
|
| 74 |
srLUTset.addUntrackedParameter<bool>("Binary", false);
|
| 75 |
srLUTset.addUntrackedParameter<std::string>("LUTPath", "./");
|
| 76 |
for(Int_t endcapItr = CSCDetId::minEndcapId(); endcapItr <= CSCDetId::maxEndcapId(); endcapItr++){
|
| 77 |
for(Int_t sectorItr = CSCTriggerNumbering::minTriggerSectorId();sectorItr <= CSCTriggerNumbering::maxTriggerSectorId();sectorItr++){
|
| 78 |
for(Int_t stationItr = 1; stationItr <= 4; stationItr++){
|
| 79 |
if(stationItr == 1){
|
| 80 |
for(Int_t subsectorItr = 0; subsectorItr < 2; subsectorItr++){
|
| 81 |
srLUTs_[endcapItr-1][sectorItr-1][subsectorItr] = new CSCSectorReceiverLUT(endcapItr, sectorItr, subsectorItr+1, stationItr, srLUTset, TMB07);
|
| 82 |
}
|
| 83 |
} else {
|
| 84 |
srLUTs_[endcapItr-1][sectorItr-1][stationItr] = new CSCSectorReceiverLUT(endcapItr, sectorItr, 0, stationItr, srLUTset, TMB07);
|
| 85 |
} //if for station 1 or 234
|
| 86 |
} // stationItr loop
|
| 87 |
} // sectorItr loop
|
| 88 |
} // endcapItr loop
|
| 89 |
*/
|
| 90 |
// ------------------------------
|
| 91 |
// --- Summary ntuple booking ---
|
| 92 |
// ------------------------------
|
| 93 |
|
| 94 |
fractNtuple = new TTree("Fraction","/Fraction");
|
| 95 |
fractNtuple->Branch("run_number", &run_number, "run_number/I") ;
|
| 96 |
fractNtuple->Branch("event_number", &event_number, "event_number/I") ;
|
| 97 |
fractNtuple->Branch("LumiBlock", &LumiBlock,"LumiBlock/I") ;
|
| 98 |
fractNtuple->Branch("bunchX", &bunchX, "bunchX/I") ;
|
| 99 |
fractNtuple->Branch("orbitNumb", &orbitNumb,"orbitNumb/I") ;
|
| 100 |
fractNtuple->Branch("mcweight", &mcweight,"mcweight/F") ;
|
| 101 |
fractNtuple->Branch("numberOfPUVertices", &numberOfPUVertices, "numberOfPUVertices/i");
|
| 102 |
fractNtuple->Branch("numberOfPUVerticesMixingTruth", &numberOfPUVerticesMixingTruth, "numberOfPUVerticesMixingTruth/F");
|
| 103 |
fractNtuple->Branch("numberOfPUVerticesTot", &numberOfPUVerticesTot, "numberOfPUVerticesTot/i");
|
| 104 |
fractNtuple->Branch("numberOfPrimaryVertices" , &numberOfPrimaryVertices , "numberOfPrimaryVertices/i");
|
| 105 |
fractNtuple->Branch("trgSingle", &trgSingle,"trgSingle/O") ;
|
| 106 |
fractNtuple->Branch("nTrkCountCSCSeg", &nTrkCountCSCSeg, "nTrkCountCSCSeg/I") ;
|
| 107 |
//fractNtuple->Branch("nTrkCountRPCE", &nTrkCountRPCE, "nTrkCountRPCE/I") ;
|
| 108 |
//fractNtuple->Branch("nTrkCountEC", &nTrkCountEC, "nTrkCountEC/I") ;
|
| 109 |
//fractNtuple->Branch("nPosTrk", &nPosTrk, "nPosTrk/I") ;
|
| 110 |
//bin/fractNtuple->Branch("nNegTrk", &nNegTrk, "nNegTrk/I") ;
|
| 111 |
fractNtuple->Branch("nTotalTrks", &nTotalTrks, "nTotalTrks/I") ;
|
| 112 |
fractNtuple->Branch("trackVeto_strict", &trackVeto_strict, "trackVeto_strict/O") ;//strick track veto: if two tracks are expected to pass one chamber, veto both (value=true).
|
| 113 |
fractNtuple->Branch("trackVeto_isClosestToLCT", &trackVeto_isClosestToLCT, "trackVeto_isClosestToLCT/O") ;//if there is a track closer to the LCT of the chamber, veto the further ones (value=true).
|
| 114 |
fractNtuple->Branch("myRegion", &myRegion, "myRegion/I") ;
|
| 115 |
fractNtuple->Branch("MuTagPt", &MuTagPt, "MuTagPt/F") ;
|
| 116 |
fractNtuple->Branch("MuTagEta", &MuTagEta, "MuTagEta/F") ;
|
| 117 |
fractNtuple->Branch("MuTagPhi", &MuTagPhi, "MuTagPhi/F") ;
|
| 118 |
fractNtuple->Branch("MuTagIsoR03Ratio", &MuTagIsoR03Ratio, "MuTagIsoR03Ratio/F") ;
|
| 119 |
fractNtuple->Branch("MuTagIsoR05Ratio", &MuTagIsoR05Ratio, "MuTagIsoR05Ratio/F") ;
|
| 120 |
fractNtuple->Branch("MuTagPromt", &MuTagPromt, "MuTagPromt/I");
|
| 121 |
fractNtuple->Branch("MuTagnSegTrkArb", &MuTagnSegTrkArb, "MuTagnSegTrkArb/I");
|
| 122 |
fractNtuple->Branch("MuTagCaloL", &MuTagCaloL, "MuTagCaloL/O") ;
|
| 123 |
fractNtuple->Branch("MuTagCaloT", &MuTagCaloT, "MuTagCaloT/O") ;
|
| 124 |
fractNtuple->Branch("MuTagtracktruth_pt", &MuTagtracktruth_pt, "MuTagtracktruth_pt/F") ;
|
| 125 |
fractNtuple->Branch("MuTagtracktruth_p", &MuTagtracktruth_p, "MuTagtracktruth_p/F") ;
|
| 126 |
fractNtuple->Branch("MuTagtracktruth_id", &MuTagtracktruth_id, "MuTagtracktruth_id/F") ;
|
| 127 |
fractNtuple->Branch("MuTagtracktruth_type", &MuTagtracktruth_type, "MuTagtracktruth_type/l") ;
|
| 128 |
fractNtuple->Branch("MuTagtracktruth_isPileup", &MuTagtracktruth_isPileup, "MuTagtracktruth_isPileup/O") ;
|
| 129 |
fractNtuple->Branch("MuTagtracktruth_thesamewith", &MuTagtracktruth_thesamewith, "MuTagtracktruth_thesamewith/I") ;
|
| 130 |
fractNtuple->Branch("vtx_r", &vtx_r, "vtx_r/F") ;
|
| 131 |
fractNtuple->Branch("vtx_z", &vtx_z, "vtx_z/F") ;
|
| 132 |
fractNtuple->Branch("vtx_rError", &vtx_rError, "vtx_rError/F") ;
|
| 133 |
fractNtuple->Branch("vtx_zError", &vtx_zError, "vtx_zError/F") ;
|
| 134 |
fractNtuple->Branch("vtx_normChi2", &vtx_normChi2, "vtx_normChi2/F") ;
|
| 135 |
fractNtuple->Branch("vtx_size", &vtx_size, "vtx_size/I") ;
|
| 136 |
fractNtuple->Branch("iSameVtx", &iSameVtx, "iSameVtx/O") ;
|
| 137 |
fractNtuple->Branch("invMass", &invMass, "invMass/F");
|
| 138 |
fractNtuple->Branch("tracks_e", &tracks_e, "tracks_e/F");
|
| 139 |
fractNtuple->Branch("tracks_pt", &tracks_pt, "tracks_pt/F");
|
| 140 |
fractNtuple->Branch("tracks_eta", &tracks_eta, "tracks_eta/F");
|
| 141 |
fractNtuple->Branch("tracks_phi", &tracks_phi, "tracks_phi/F");
|
| 142 |
fractNtuple->Branch("tracks_charge", &tracks_charge, "tracks_charge/I");
|
| 143 |
fractNtuple->Branch("tracks_id", &tracks_id, "tracks_id/I");
|
| 144 |
fractNtuple->Branch("tracks_normchi2", &tracks_chi2, "tracks_normchi2/F");
|
| 145 |
fractNtuple->Branch("tracks_dxy", &tracks_dxy, "tracks_dxy/F");
|
| 146 |
fractNtuple->Branch("tracks_dz", &tracks_dz, "tracks_dz/F");
|
| 147 |
fractNtuple->Branch("tracks_vx", &tracks_vx, "tracks_vx/F");
|
| 148 |
fractNtuple->Branch("tracks_vy", &tracks_vy, "tracks_vy/F");
|
| 149 |
fractNtuple->Branch("tracks_vz", &tracks_vz, "tracks_vz/F");
|
| 150 |
fractNtuple->Branch("tracks_IsoR03Ratio", &tracks_IsoR03Ratio, "tracks_IsoR03Ratio/F");
|
| 151 |
fractNtuple->Branch("tracks_IsoR05Ratio", &tracks_IsoR05Ratio, "tracks_IsoR05Ratio/F");
|
| 152 |
fractNtuple->Branch("tracks_qoverp", &tracks_qoverp, "tracks_qoverp/F");
|
| 153 |
fractNtuple->Branch("tracks_lambda", &tracks_lambda, "tracks_lambda/F");
|
| 154 |
fractNtuple->Branch("tracks_recHitsSize", &tracks_recHitsSize, "tracks_recHitsSize/I");
|
| 155 |
fractNtuple->Branch("tracks_numberOfValidHits", &tracks_numberOfValidHits, "tracks_numberOfValidHits/I");
|
| 156 |
fractNtuple->Branch("tracks_numberOfLostHits", &tracks_numberOfLostHits, "tracks_numberOfLostHits/I");
|
| 157 |
fractNtuple->Branch("tracks_qoverpError", &tracks_qoverpError, "tracks_qoverpError/F");
|
| 158 |
fractNtuple->Branch("tracks_ptError", &tracks_ptError, "tracks_ptError/F");
|
| 159 |
fractNtuple->Branch("tracks_thetaError", &tracks_thetaError, "tracks_thetaError/F");
|
| 160 |
fractNtuple->Branch("tracks_lambdaError", &tracks_lambdaError, "tracks_lambdaError/F");
|
| 161 |
fractNtuple->Branch("tracks_etaError", &tracks_etaError, "tracks_etaError/F");
|
| 162 |
fractNtuple->Branch("tracks_phiError", &tracks_phiError, "tracks_phiError/F");
|
| 163 |
fractNtuple->Branch("tracks_dxyError", &tracks_dxyError, "tracks_dxyError/F");
|
| 164 |
fractNtuple->Branch("tracks_d0Error", &tracks_d0Error, "tracks_d0Error/F");
|
| 165 |
fractNtuple->Branch("tracks_dszError", &tracks_dszError, "tracks_dszError/F");
|
| 166 |
fractNtuple->Branch("tracks_dzError", &tracks_dzError, "tracks_dzError/F");
|
| 167 |
fractNtuple->Branch("tracks_isCaloMuTrk", &tracks_isCaloMuTrk, "tracks_isCaloMuTrk/O");
|
| 168 |
fractNtuple->Branch("tracks_isTrackerMuTrk", &tracks_isTrackerMuTrk, "tracks_isTrackerMuTrk/O");
|
| 169 |
// fractNtuple->Branch("tracks_dedx", &tracks_dedx, "tracks_dedx/F") ;
|
| 170 |
// fractNtuple->Branch("tracks_dedxErr", &tracks_dedxErr, "tracks_dedxErr/F") ;
|
| 171 |
// fractNtuple->Branch("tracks_nSatMeas", &tracks_nSatMeas, "tracks_nSatMeas/F") ;
|
| 172 |
//fractNtuple->Branch("tracks_nMeas", &tracks_nMeas, "tracks_nMeas/F") ;
|
| 173 |
fractNtuple->Branch("tracktruth_pt", &tracktruth_pt, "tracktruth_pt/F") ;
|
| 174 |
fractNtuple->Branch("tracktruth_p", &tracktruth_p, "tracktruth_p/F") ;
|
| 175 |
fractNtuple->Branch("tracktruth_e", &tracktruth_e, "tracktruth_e/F") ;
|
| 176 |
fractNtuple->Branch("tracktruth_id", &tracktruth_id, "tracktruth_id/F") ;
|
| 177 |
fractNtuple->Branch("tracktruth_type", &tracktruth_type, "tracktruth_type/l") ;
|
| 178 |
fractNtuple->Branch("tracktruth_isPileup", &tracktruth_isPileup, "tracktruth_isPileup/O") ;
|
| 179 |
fractNtuple->Branch("tracktruth_thesamewith", &tracktruth_thesamewith, "tracktruth_thesamewith/I") ;
|
| 180 |
|
| 181 |
fractNtuple->Branch("CSCEndCapPlus", &CSCEndCapPlus, "CSCEndCapPlus/O") ;
|
| 182 |
|
| 183 |
#define MakeBranchAllSt(Name,Type,Var) sprintf(BranchName,"%s%d",Name,st+1); \
|
| 184 |
sprintf(BranchTitle,"%s%d/%s",Name,st+1,Type); \
|
| 185 |
fractNtuple->Branch(BranchName, &Var[st],BranchTitle)
|
| 186 |
|
| 187 |
for (Short_t st=0;st<4;st++) {
|
| 188 |
char BranchName[30],BranchTitle[30];
|
| 189 |
/*CSC Chamber Candidates in each station*/
|
| 190 |
MakeBranchAllSt("CSCRg","b",CSCRg);
|
| 191 |
MakeBranchAllSt("CSCCh","b",CSCChCand);
|
| 192 |
MakeBranchAllSt("CSCCBad","O",CSCChBad);
|
| 193 |
|
| 194 |
/*Extrapolated Tracks on CSC Chamber Candidates in each station*/
|
| 195 |
MakeBranchAllSt("CSCDyProjHVGap","F",CSCDyProjHVGap);
|
| 196 |
MakeBranchAllSt("CSCDyErrProjHVGap","F",CSCDyErrProjHVGap);
|
| 197 |
MakeBranchAllSt("CSCProjDistEdge","F",CSCProjEdgeDist);
|
| 198 |
MakeBranchAllSt("CSCProjDistErrEdge","F",CSCProjEdgeDistErr);
|
| 199 |
|
| 200 |
/*Segments characteristics*/
|
| 201 |
MakeBranchAllSt("CSCSegxLc","F",CSCSegxLc);
|
| 202 |
MakeBranchAllSt("CSCSegyLc","F",CSCSegyLc);
|
| 203 |
MakeBranchAllSt("CSCSegxErrLc","F",CSCSegxErrLc);
|
| 204 |
MakeBranchAllSt("CSCSegyErrLc","F",CSCSegyErrLc);
|
| 205 |
MakeBranchAllSt("CSCdXdZTTSeg","F",CSCdXdZTTSeg);
|
| 206 |
MakeBranchAllSt("CSCdYdZTTSeg","F",CSCdYdZTTSeg);
|
| 207 |
MakeBranchAllSt("CSCSegChisqProb","F",CSCSegChisqProb);
|
| 208 |
MakeBranchAllSt("CSCnSegHits","I",CSCnSegHits);
|
| 209 |
|
| 210 |
/*Distance from the Extrapolated Tracks to CSC Segments, 99999. for no CSC segment found*/
|
| 211 |
MakeBranchAllSt("CSCDxTTSeg","F",CSCDxTTSeg);
|
| 212 |
MakeBranchAllSt("CSCDxErrTTSeg","F",CSCDxErrTTSeg);
|
| 213 |
MakeBranchAllSt("CSCDyTTSeg","F",CSCDyTTSeg);
|
| 214 |
MakeBranchAllSt("CSCDyErrTTSeg","F",CSCDyErrTTSeg);
|
| 215 |
MakeBranchAllSt("CSCDxyTTSeg","F",CSCDxyTTSeg);
|
| 216 |
MakeBranchAllSt("CSCDxyErrTTSeg","F",CSCDxyErrTTSeg);
|
| 217 |
|
| 218 |
/*LCT characteristics*/
|
| 219 |
MakeBranchAllSt("CSCLCTxLc","F",CSCLCTxLc);
|
| 220 |
MakeBranchAllSt("CSCLCTyLc","F",CSCLCTyLc);
|
| 221 |
MakeBranchAllSt("CSCLCTbx","I",CSCLCTbx);
|
| 222 |
|
| 223 |
/*Distance from the Extrapolated Tracks to LCT, 99999. for no LCT found*/
|
| 224 |
MakeBranchAllSt("CSCDxTTLCT","F",CSCDxTTLCT);
|
| 225 |
MakeBranchAllSt("CSCDxErrTTLCT","F",CSCDxErrTTLCT);
|
| 226 |
MakeBranchAllSt("CSCDyTTLCT","F",CSCDyTTLCT);
|
| 227 |
MakeBranchAllSt("CSCDxErrTTLCT","F",CSCDxErrTTLCT);
|
| 228 |
MakeBranchAllSt("CSCDxyTTLCT","F",CSCDxyTTLCT);
|
| 229 |
MakeBranchAllSt("CSCDxyErrTTLCT","F",CSCDxyErrTTLCT);
|
| 230 |
|
| 231 |
/*DetlaR between the extrapolated tracker track on muon system and the tagged muon*/
|
| 232 |
MakeBranchAllSt("dRTkMu","F",dRTkMu);
|
| 233 |
|
| 234 |
/*Default decision of whether a segment or LCT is found*/
|
| 235 |
// MakeBranchAllSt("segSt","I",segSt);
|
| 236 |
// MakeBranchAllSt("lctSt","I",lctSt);
|
| 237 |
}
|
| 238 |
|
| 239 |
HLTMuAcceptance=new vector<Bool_t>();
|
| 240 |
fractNtuple->Branch("HLTMuAcceptance", &HLTMuAcceptance);
|
| 241 |
fractNtuple->Branch("HLTDiMuAcceptance", &HLTDiMuAcceptance, "HLTDiMuAcceptance/F") ;
|
| 242 |
minDRHLTMu=new vector<Float_t>();
|
| 243 |
fractNtuple->Branch("minDRHLTMu", &minDRHLTMu) ;
|
| 244 |
fractNtuple->Branch("minDRHLTDiMu", &minDRHLTDiMu, "minDRHLTDiMu/F") ;
|
| 245 |
fractNtuple->Branch("minDRHLTAllSingleMu", &minDRHLTAllSingleMu, "minDRHLTAllSingleMu/F") ;
|
| 246 |
|
| 247 |
RunInfo = new TTree("RunInfo","RunInfo");
|
| 248 |
RunInfo->Branch("RunNumber",&run_number,"RunNum/i");
|
| 249 |
RunInfo->Branch("TableName",&HLTTableName);
|
| 250 |
HLTMuNames=new vector<string>();
|
| 251 |
RunInfo->Branch("HLTMuNames",&HLTMuNames);
|
| 252 |
HLTMuObjModuleNames=new vector<string>();
|
| 253 |
RunInfo->Branch("HLTMuObjModuleNames",&HLTMuObjModuleNames);
|
| 254 |
RunInfo->Branch("HLTDiMuName",&HLTDiMuName);
|
| 255 |
RunInfo->Branch("HLTDiMuObjModuleName",&HLTDiMuObjModuleName);
|
| 256 |
RunInfo->Branch("BadCSCChamberIndexList",&badChambersIndices);
|
| 257 |
}
|
| 258 |
|
| 259 |
TPTrackMuonSys::~TPTrackMuonSys(){
|
| 260 |
// do anything here that needs to be done at desctruction time
|
| 261 |
// (e.g. close files, deallocate resources etc.)
|
| 262 |
}
|
| 263 |
|
| 264 |
// ------------ method called to for each event ------------
|
| 265 |
void
|
| 266 |
TPTrackMuonSys::analyze(const edm::Event& event, const edm::EventSetup& setup){
|
| 267 |
//cout << "HERE *************************" << endl;
|
| 268 |
|
| 269 |
nEventsAnalyzed++;
|
| 270 |
// access conditions data for this event
|
| 271 |
theDbConditions.initializeEvent( setup ); //initializeEvent(const edm::EventSetup & es);// fetch the maps from the database
|
| 272 |
|
| 273 |
//////////////////////////////////////////////////////////////////////////
|
| 274 |
reco::BeamSpot beamSpot;
|
| 275 |
edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
|
| 276 |
Bool_t beamExists = false;
|
| 277 |
try{
|
| 278 |
event.getByLabel(m_beamSpot,recoBeamSpotHandle);
|
| 279 |
if ( recoBeamSpotHandle.isValid()) {
|
| 280 |
beamSpot = *recoBeamSpotHandle;
|
| 281 |
beamExists = true;
|
| 282 |
}
|
| 283 |
}catch (cms::Exception){
|
| 284 |
edm::LogError("")<< "Error! Can't get m_beamSpot by label. ";
|
| 285 |
}
|
| 286 |
event_number = event.id().event();
|
| 287 |
LumiBlock = (int)event.luminosityBlock();
|
| 288 |
expType = (int)event.experimentType();
|
| 289 |
bunchX = (int)event.bunchCrossing();
|
| 290 |
orbitNumb = (int)event.orbitNumber();
|
| 291 |
Bool_t isRData = event.isRealData(); //( L1Decision ? "passed" : "failed")
|
| 292 |
|
| 293 |
//Get the Magnetic field from the setup
|
| 294 |
setup.get<IdealMagneticFieldRecord>().get(theBField);
|
| 295 |
// Get the GlobalTrackingGeometry from the setup
|
| 296 |
// setup.get<GlobalTrackingGeometryRecord>().get(theTrackingGeometry);
|
| 297 |
// Get the DT Geometry from the setup
|
| 298 |
// setup.get<MuonGeometryRecord>().get(dtGeom);
|
| 299 |
|
| 300 |
setup.get<MuonGeometryRecord>().get(cscGeom);
|
| 301 |
//#define m_debug
|
| 302 |
#ifdef m_debug
|
| 303 |
//check the chamber active region geometry and numbering --- something good to know
|
| 304 |
if (nEventsAnalyzed==1) {
|
| 305 |
UChar_t rings[]={11,12,21,31,41,22};//,14,12,13,21,31,41,22};
|
| 306 |
for (UChar_t iendcap=1;iendcap<3;iendcap++)
|
| 307 |
for (UChar_t iring=0;iring<sizeof(rings)/sizeof(UChar_t);iring++) {
|
| 308 |
CSCDetId id=CSCDetId(iendcap, rings[iring]/10, rings[iring]%10, 1, 1);
|
| 309 |
printf("============ ME%c%d chamber outline=============\n",iendcap==1?'+':'-',rings[iring]);
|
| 310 |
const CSCLayerGeometry *layerGeom = cscGeom->chamber(id)->layer(1)->geometry ();
|
| 311 |
const UChar_t NStrips=layerGeom->numberOfStrips(),NWires=layerGeom->numberOfWireGroups();
|
| 312 |
LocalPoint interSect_ = layerGeom->stripWireGroupIntersection(1, 5);
|
| 313 |
printf("(strip 1, wire group 5) at (%.3f,%.3f)\n",interSect_.x(),interSect_.y());
|
| 314 |
interSect_ = layerGeom->stripWireGroupIntersection(NStrips, 5);
|
| 315 |
printf("(strip %d, wire group 5) at (%.3f,%.3f)\n",NStrips,interSect_.x(),interSect_.y());
|
| 316 |
interSect_ = layerGeom->stripWireGroupIntersection(NStrips, NWires);
|
| 317 |
printf("(strip %d, wire group %d) at (%.3f,%.3f)\n",NStrips, NWires,interSect_.x(),interSect_.y());
|
| 318 |
interSect_ = layerGeom->stripWireGroupIntersection(1, NWires);
|
| 319 |
printf("(strip 1, wire group %d) at (%.3f,%.3f)\n", NWires,interSect_.x(),interSect_.y());
|
| 320 |
printf(" ======== middle x of wires(check)==========\n");
|
| 321 |
const CSCWireTopology* wireTopology = layerGeom->wireTopology();
|
| 322 |
for (UChar_t wireGroup_id=10;wireGroup_id<=15;wireGroup_id++) {
|
| 323 |
Float_t wideWidth = wireTopology->wideWidthOfPlane();
|
| 324 |
Float_t narrowWidth = wireTopology->narrowWidthOfPlane();
|
| 325 |
Float_t length = wireTopology->lengthOfPlane();
|
| 326 |
Float_t tangent = (wideWidth-narrowWidth)/(2.*length);
|
| 327 |
Float_t wireangle= wireTopology->wireAngle();
|
| 328 |
std::vector<float> wirecenters=wireTopology->wireValues(wireTopology->middleWireOfGroup(wireGroup_id));
|
| 329 |
printf("x center of wire %d is at %.3f=%.3f\n",wireGroup_id,wireTopology->wireValues(wireTopology->middleWireOfGroup(wireGroup_id))[0],wirecenters[2]*sin(wireangle)*tangent/2);
|
| 330 |
}
|
| 331 |
printf("y center of the wires is at %.3f\n",wireTopology->yOfWire(1)+0.5*wireTopology->lengthOfPlane());
|
| 332 |
}
|
| 333 |
}
|
| 334 |
#endif
|
| 335 |
// Get the propagators
|
| 336 |
setup.get<TrackingComponentsRecord>().get("SmartPropagatorAnyRK", propagatorAlong );
|
| 337 |
setup.get<TrackingComponentsRecord>().get("SmartPropagatorAnyOpposite", propagatorOpposite);
|
| 338 |
|
| 339 |
//vertices
|
| 340 |
edm::Handle<reco::VertexCollection> pvHandle;
|
| 341 |
try{
|
| 342 |
event.getByLabel(m_vertexSrc, pvHandle);
|
| 343 |
}catch (cms::Exception){
|
| 344 |
edm::LogError("")<< "Error! Can't get m_vertexSrc by label. ";
|
| 345 |
}
|
| 346 |
const reco::VertexCollection & vertices = *pvHandle.product();
|
| 347 |
numberOfPrimaryVertices =pvHandle ->size();
|
| 348 |
|
| 349 |
// Get the RecHits collection :
|
| 350 |
Handle<CSCRecHit2DCollection> recHits;
|
| 351 |
event.getByLabel("csc2DRecHits",recHits);
|
| 352 |
|
| 353 |
// get CSC segment collection
|
| 354 |
event.getByLabel("cscSegments", cscSegments);
|
| 355 |
|
| 356 |
edm::Handle<CSCCorrelatedLCTDigiCollection> mpclcts;
|
| 357 |
try{
|
| 358 |
event.getByLabel("csctfunpacker","", mpclcts);
|
| 359 |
}catch (cms::Exception){
|
| 360 |
edm::LogError("")<< "Error! Can't get m_gTracksTag by label. ";
|
| 361 |
}
|
| 362 |
|
| 363 |
/*
|
| 364 |
// get RPC and DT
|
| 365 |
setup.get<MuonGeometryRecord>().get(rpcGeo);
|
| 366 |
|
| 367 |
edm::Handle<RPCRecHitCollection> rpcRecHits;
|
| 368 |
event.getByLabel("rpcRecHits",rpcRecHits);
|
| 369 |
|
| 370 |
edm::Handle<DTRecSegment4DCollection> dtSegments;
|
| 371 |
event.getByLabel( "dt4DSegments", dtSegments );
|
| 372 |
*/
|
| 373 |
|
| 374 |
// muon trigger scales
|
| 375 |
edm::ESHandle<L1MuTriggerScales> trigscales_h;
|
| 376 |
setup.get<L1MuTriggerScalesRcd> ().get(trigscales_h);
|
| 377 |
theTriggerScales = trigscales_h.product();
|
| 378 |
|
| 379 |
//if (! dtSegments.isValid()) throw cms::Exception("FatalError") << "Unable to find DTRecSegment4DCollection in event!\n";
|
| 380 |
|
| 381 |
std::vector<DetId> chamberIds;
|
| 382 |
std::vector<DetId>::const_iterator chamberIdIt;
|
| 383 |
|
| 384 |
Handle<reco::TrackCollection> gTracks;
|
| 385 |
// Bool_t trkBool[2]; trkBool[0]= true;
|
| 386 |
try{
|
| 387 |
event.getByLabel(m_gTracksTag, gTracks);
|
| 388 |
}catch (cms::Exception){
|
| 389 |
edm::LogError("")<< "Error! Can't get m_gTracksTag by label. ";
|
| 390 |
// trkBool[0] = false;
|
| 391 |
}
|
| 392 |
|
| 393 |
try{
|
| 394 |
event.getByLabel("muons", muons);
|
| 395 |
}catch (cms::Exception){
|
| 396 |
edm::LogError("")<< "Error! Can't get muons ";
|
| 397 |
}
|
| 398 |
Handle<ValueMap<DeDxData> > dEdxTrackHandle;
|
| 399 |
event.getByLabel(m_dEdxDiscrimTag, dEdxTrackHandle);
|
| 400 |
const ValueMap<DeDxData> dEdxTrack = *dEdxTrackHandle.product();
|
| 401 |
|
| 402 |
// calo geometry
|
| 403 |
edm::ESHandle<CaloGeometry> calogeo;
|
| 404 |
setup.get<CaloGeometryRecord>().get(calogeo);
|
| 405 |
if (! calogeo.isValid()) throw cms::Exception("FatalError") << "Unable to find CaloGeometryRecord in event!\n";
|
| 406 |
|
| 407 |
//=======Generation and Simulation information========
|
| 408 |
Handle< View<Track> > trackCollectionHV;
|
| 409 |
event.getByLabel(m_gTracksTag,trackCollectionHV);
|
| 410 |
//Simulated Vertices: the vertexId() is just the index
|
| 411 |
Handle<SimVertexContainer> SVCollectionH;
|
| 412 |
HepMC::GenEvent * HepGenEvent=NULL;
|
| 413 |
|
| 414 |
//TrackingParticles
|
| 415 |
Handle<TrackingParticleCollection> TPCollectionH ;
|
| 416 |
RecoToSimCollection RecoToSimByHits;
|
| 417 |
ESHandle<TrackAssociatorBase> AssociatorByHits;
|
| 418 |
|
| 419 |
vector< vector< vector<Int_t> > > TracksSimChains, MuTracksSimChains;
|
| 420 |
|
| 421 |
numberOfPUVertices=0; // the number of pileup interactions that have been added to the event from this bunch crossing
|
| 422 |
numberOfPUVerticesMixingTruth=0;// the *true* mean number of pileup interactions for this event from which each bunch crossing has been sampled; same for all bunch crossings in an event
|
| 423 |
numberOfPUVerticesTot=0;
|
| 424 |
mcweight = 1.;
|
| 425 |
Handle<edm::HepMCProduct> HepMCH;
|
| 426 |
if (m_isMC) {
|
| 427 |
event.getByLabel("g4SimHits", SVCollectionH);
|
| 428 |
SVC = *SVCollectionH.product();
|
| 429 |
event.getByLabel(m_HepMCTag, HepMCH);
|
| 430 |
HepGenEvent=const_cast<HepMC::GenEvent *>( HepMCH->GetEvent() );
|
| 431 |
Handle<GenEventInfoProduct> hEvtInfo;
|
| 432 |
event.getByLabel("generator", hEvtInfo);
|
| 433 |
mcweight = hEvtInfo->weight();
|
| 434 |
event.getByLabel("mergedtruth",TPCollectionH);
|
| 435 |
//event.getByType(TPCollectionH);
|
| 436 |
setup.get<TrackAssociatorRecord>().get("TrackAssociatorByHits", AssociatorByHits);
|
| 437 |
RecoToSimByHits = AssociatorByHits->associateRecoToSim(trackCollectionHV,TPCollectionH,&event,&setup);
|
| 438 |
MCParticlesList.clear();
|
| 439 |
SavedSimTrk.clear();
|
| 440 |
SavedHepPar.clear();
|
| 441 |
/*Start of getting pileup information*/
|
| 442 |
edm::InputTag PileupSrc_("addPileupInfo");
|
| 443 |
Handle<std::vector< PileupSummaryInfo > > puInfo;
|
| 444 |
if ( event.getByLabel(PileupSrc_, puInfo) ) {
|
| 445 |
#ifdef GetPUFromEarlierThan_4_1_2
|
| 446 |
numberOfPUVertices = (*puInfo)[0].getPU_NumInteractions();
|
| 447 |
numberOfPUVerticesMixingTruth=numberOfPUVertices;
|
| 448 |
numberOfPUVerticesTot=numberOfPUVertices;
|
| 449 |
#else
|
| 450 |
for(std::vector<PileupSummaryInfo>::const_iterator PVI = puInfo->begin(); PVI != puInfo->end(); ++PVI) {
|
| 451 |
Int_t BX = PVI->getBunchCrossing();
|
| 452 |
if(BX == 0) {
|
| 453 |
numberOfPUVertices = PVI->getPU_NumInteractions();
|
| 454 |
#ifdef GetPUFromEarlierThan_4_4_0
|
| 455 |
numberOfPUVerticesMixingTruth = -1.;//not available
|
| 456 |
#else
|
| 457 |
numberOfPUVerticesMixingTruth = PVI->getTrueNumInteractions();
|
| 458 |
#endif
|
| 459 |
}
|
| 460 |
numberOfPUVerticesTot+=numberOfPUVertices;
|
| 461 |
}
|
| 462 |
#endif
|
| 463 |
}//====== End of getting pileup information =============
|
| 464 |
}//======= End of if m_isMC========
|
| 465 |
|
| 466 |
//////////////////////////////////////////////
|
| 467 |
edm::ESHandle<L1GtTriggerMenu> menuRcd;
|
| 468 |
setup.get<L1GtTriggerMenuRcd>().get(menuRcd) ;
|
| 469 |
const L1GtTriggerMenu* menu = menuRcd.product();
|
| 470 |
Handle<L1GlobalTriggerReadoutRecord> gtBeforeMaskRecord;
|
| 471 |
event.getByLabel("gtDigis", gtBeforeMaskRecord);
|
| 472 |
trgSingle = false;
|
| 473 |
if (!gtBeforeMaskRecord.isValid()) {
|
| 474 |
LogInfo("PhysicsTrees: ") << " Error: no L1GlobalTriggerReadoutRecord found with input tag --> gtDigis" ;
|
| 475 |
} else {
|
| 476 |
const TechnicalTriggerWord techDecisionWord = gtBeforeMaskRecord->technicalTriggerWord(0);
|
| 477 |
const DecisionWord gtDecisionWordBeforeMask = gtBeforeMaskRecord->decisionWord(0);
|
| 478 |
#ifdef m_debug
|
| 479 |
Int_t nl1s = 0;
|
| 480 |
#endif
|
| 481 |
for (CItAlgo algo = menu->gtAlgorithmMap().begin(); algo!=menu->gtAlgorithmMap().end(); ++algo) {
|
| 482 |
#ifdef m_debug
|
| 483 |
if (nEventsAnalyzed <= 1) {
|
| 484 |
cout <<" L1Menu Item m1 = " << nl1s
|
| 485 |
<<" corresponds to AlgoName = "
|
| 486 |
<< (algo->second).algoName() << " Alias: " << (algo->second).algoAlias()
|
| 487 |
<< " and the result= " << menu->gtAlgorithmResult((algo->second).algoAlias(), gtDecisionWordBeforeMask)
|
| 488 |
<< endl;
|
| 489 |
}
|
| 490 |
nl1s = nl1s + 1;
|
| 491 |
#endif
|
| 492 |
if((algo->second).algoName() == "L1_SingleMuOpen") trgSingle = true;
|
| 493 |
if((algo->second).algoName() == "L1_SingleMu5") trgSingle=false;
|
| 494 |
}
|
| 495 |
}
|
| 496 |
|
| 497 |
/*-----------Start getting HLT results------------*/
|
| 498 |
Handle<TriggerResults> triggerResults;
|
| 499 |
Handle< trigger::TriggerEvent > handleTriggerEvent;
|
| 500 |
HLTDiMuAcceptance=false;
|
| 501 |
HLTMuAcceptance->clear();
|
| 502 |
try {
|
| 503 |
event.getByLabel(m_hlt, triggerResults);
|
| 504 |
if ( triggerResults.product()->wasrun() ){
|
| 505 |
LogInfo("")<<" At least one path out of " << triggerResults.product()->size() << " ran? " << triggerResults.product()->wasrun() << endl;
|
| 506 |
if ( triggerResults.product()->accept() ) {
|
| 507 |
for (vector<Int_t>::const_iterator iter=m_HLTMuTrgBit.begin();iter<m_HLTMuTrgBit.end();iter++)
|
| 508 |
HLTMuAcceptance->push_back( triggerResults.product()->accept(*iter) );
|
| 509 |
if ( m_HLTDiMuTrgBit>-1 ) HLTDiMuAcceptance=triggerResults.product()->accept( m_HLTDiMuTrgBit );
|
| 510 |
} // end if at least one triggerResult accepted
|
| 511 |
}// end if wasRun
|
| 512 |
} catch (...) {// some old codes, it is possibly not going to happen in >CMSSW_4_4_X versions
|
| 513 |
LogWarning("")<<"Not found HLT result, is HLT process name correct?"<< endl;
|
| 514 |
if(isRData == 0){
|
| 515 |
if(event.getByLabel( m_hltTrgEv, handleTriggerEvent )){
|
| 516 |
const trigger::TriggerObjectCollection & toc(handleTriggerEvent->getObjects());
|
| 517 |
for ( size_t ia = 0; ia < handleTriggerEvent->sizeFilters(); ++ ia) {
|
| 518 |
std::string fullname = handleTriggerEvent->filterTag(ia).encode();
|
| 519 |
std::string name;
|
| 520 |
size_t p = fullname.find_first_of(':');
|
| 521 |
if ( p != std::string::npos) name = fullname.substr(0, p);
|
| 522 |
else name = fullname;
|
| 523 |
if ( &toc !=0 ) {
|
| 524 |
HLTMuAcceptance->assign(m_HLTMuTrgBit.size(),false);
|
| 525 |
for (vector<string>::const_iterator iter=HLTMuObjModuleNames->begin(); iter<HLTMuObjModuleNames->end();iter++ )
|
| 526 |
if ( name == *iter ) HLTMuAcceptance->at( iter-HLTMuObjModuleNames->begin() )=true;
|
| 527 |
if ( name == *HLTDiMuObjModuleName ) HLTDiMuAcceptance=true;
|
| 528 |
}
|
| 529 |
}
|
| 530 |
}
|
| 531 |
}
|
| 532 |
}
|
| 533 |
Bool_t m_GotTrgObj=event.getByLabel( m_hltTrgEv, handleTriggerEvent );
|
| 534 |
/*-----------End getting HLT results------------*/
|
| 535 |
|
| 536 |
if(!gTracks.isValid() ) return;
|
| 537 |
if(!muons.isValid() ) return;
|
| 538 |
if (gTracks->size() > MAXNTRACKS) return;
|
| 539 |
|
| 540 |
#ifdef m_debug
|
| 541 |
cout << "Run: " << run_number << " Event: " << event_number << " LumiBlock: " << LumiBlock << " BX: " << bunchX
|
| 542 |
<< " isRealDAT: " << isRData
|
| 543 |
<< endl;
|
| 544 |
#endif
|
| 545 |
|
| 546 |
Int_t trackNumber[2][4][4][36];
|
| 547 |
Int_t BestTrkWithLCT[2][4][4][36];
|
| 548 |
Float_t CloestDisTrkLCT[2][4][4][36];
|
| 549 |
Bool_t trkVeto[MAXNTRACKS];
|
| 550 |
|
| 551 |
for(UChar_t i1 = 0; i1 < 2; i1++)
|
| 552 |
for(UChar_t i2 = 0; i2 < 4; i2++)
|
| 553 |
for(UChar_t i3 = 0; i3 < 4; i3++)
|
| 554 |
for(UChar_t i4 = 0; i4 < 36; i4++) {
|
| 555 |
trackNumber[i1][i2][i3][i4] = -1;
|
| 556 |
BestTrkWithLCT[i1][i2][i3][i4] = -1;
|
| 557 |
CloestDisTrkLCT[i1][i2][i3][i4] = 9999.;
|
| 558 |
}
|
| 559 |
|
| 560 |
nTotalTrks = gTracks->size();
|
| 561 |
for(reco::TrackCollection::const_iterator itTrack = gTracks->begin(); itTrack != gTracks->end(); itTrack++){//start loop tracks
|
| 562 |
UInt_t itrk=itTrack - gTracks->begin();
|
| 563 |
trkVeto[ itrk ]=false;
|
| 564 |
if ( itTrack->charge() == 0 ) continue;
|
| 565 |
if ( itTrack->p() < 3.0) continue;
|
| 566 |
|
| 567 |
reco::TrackRef trackRef(gTracks, itrk );
|
| 568 |
Float_t trk_eta = itTrack->eta(),
|
| 569 |
trk_phi = itTrack->phi(),
|
| 570 |
trk_dxy = itTrack->dxy(),
|
| 571 |
trk_dz = itTrack->dz(),
|
| 572 |
trk_chi2 = 9999.;
|
| 573 |
while (trk_phi < 0 ) trk_phi = trk_phi + 2*M_PI;
|
| 574 |
if(itTrack->ndof() != 0) trk_chi2 = itTrack->chi2()/itTrack->ndof();
|
| 575 |
if(beamExists){
|
| 576 |
trk_dxy = itTrack->dxy(beamSpot.position());
|
| 577 |
trk_dz = itTrack->dz(beamSpot.position());
|
| 578 |
}
|
| 579 |
if ( ! (itTrack->numberOfValidHits()>7 && fabs(trk_dz)<24 && fabs(trk_dxy)< 2.0 && trk_chi2>0 && trk_chi2<4) ) continue;
|
| 580 |
|
| 581 |
Bool_t ec= ( trk_eta > 0 );
|
| 582 |
// if(fabs(trk_eta) < 1.2 ) trackDT = true;
|
| 583 |
if(fabs(trk_eta) > 0.9 ) {//if it is a CSCTrack
|
| 584 |
for(UChar_t j =0; j < 6; j++){//start looping stations
|
| 585 |
TrajectoryStateOnSurface tsos = surfExtrapTrkSam(trackRef, ec?MEZ[j]:-MEZ[j]);
|
| 586 |
if (!tsos.isValid()) continue;
|
| 587 |
Float_t trkEta = tsos.globalPosition().eta(), trkPhi = tsos.globalPosition().phi();
|
| 588 |
UChar_t endcapCSC=ec?0:1,
|
| 589 |
stationCSC = j>2?j-2:0,
|
| 590 |
ringCSC = ringCandidate(stationCSC+1, trkEta, trkPhi);
|
| 591 |
if (!ringCSC) continue;
|
| 592 |
UChar_t chamberCSC = thisChamberCandidate(stationCSC+1, ringCSC, trkPhi)-1;
|
| 593 |
if ( trackNumber[endcapCSC][stationCSC][ringCSC-1][chamberCSC] < 0 )
|
| 594 |
trackNumber[endcapCSC][stationCSC][ringCSC-1][chamberCSC] = itrk;
|
| 595 |
else {
|
| 596 |
trkVeto[ trackNumber[endcapCSC][stationCSC][ringCSC-1][chamberCSC] ]=true;
|
| 597 |
trkVeto[ itrk ]=true;
|
| 598 |
}
|
| 599 |
|
| 600 |
CSCDetId Layer3id = CSCDetId( endcapCSC+1, stationCSC+1, ringCSC, chamberCSC+1, 3 );//go to layer 3 that corresponds to the LCTPos
|
| 601 |
const BoundPlane &Layer3Surface=cscGeom->idToDet(Layer3id)->surface();
|
| 602 |
tsos=surfExtrapTrkSam(trackRef, Layer3Surface.position().z());
|
| 603 |
LocalPoint localL3pCSC = Layer3Surface.toLocal(tsos.freeState()->position());
|
| 604 |
Float_t dRTrkLCT=9999.;
|
| 605 |
Int_t dummy_bx=0.;
|
| 606 |
LocalPoint *LCTPos=matchTTwithLCTs( localL3pCSC.x(), localL3pCSC.y(), endcapCSC+1, stationCSC+1, ringCSC, chamberCSC+1, mpclcts, dRTrkLCT, dummy_bx);
|
| 607 |
if (LCTPos!=NULL) {
|
| 608 |
if (dRTrkLCT<CloestDisTrkLCT[endcapCSC][stationCSC][ringCSC-1][chamberCSC]) {
|
| 609 |
CloestDisTrkLCT[endcapCSC][stationCSC][ringCSC-1][chamberCSC] = dRTrkLCT;
|
| 610 |
BestTrkWithLCT[endcapCSC][stationCSC][ringCSC-1][chamberCSC] = itrk;
|
| 611 |
}
|
| 612 |
delete LCTPos;
|
| 613 |
}
|
| 614 |
}//end looping stations
|
| 615 |
}//end of track loop
|
| 616 |
}//end of if CSC track
|
| 617 |
|
| 618 |
std::string trackExtractorName = trackExtractorPSet_.getParameter<std::string>("ComponentName");
|
| 619 |
reco::isodeposit::IsoDepositExtractor* muIsoExtractorTrack_ = IsoDepositExtractorFactory::get()->create( trackExtractorName, trackExtractorPSet_);
|
| 620 |
for (reco::MuonCollection::const_iterator muIter1 = muons->begin(); muIter1 != muons->end(); ++muIter1) {
|
| 621 |
if (!muIter1->isTrackerMuon() ) continue;
|
| 622 |
if(!muIter1->track().isNonnull()) continue;
|
| 623 |
MuTagPt = muIter1->track()->pt() ;
|
| 624 |
if (MuTagPt < 5.) continue;
|
| 625 |
|
| 626 |
Float_t mu1dxy = 1000., mu1dz = 1000.;
|
| 627 |
|
| 628 |
if(beamExists){
|
| 629 |
mu1dxy = muIter1->track()->dxy(beamSpot.position());
|
| 630 |
mu1dz = muIter1->track()->dz(beamSpot.position());
|
| 631 |
}else{
|
| 632 |
mu1dxy = muIter1->track()->dxy();
|
| 633 |
mu1dz = muIter1->track()->dz();
|
| 634 |
}
|
| 635 |
Float_t mu1Chi2 = muIter1->track()->normalizedChi2();
|
| 636 |
MuTagHitsTrkSys = muIter1->track()->hitPattern().numberOfValidTrackerHits();
|
| 637 |
Bool_t goodTrack = (fabs(mu1dxy) < 2.0 && fabs(mu1dz) < 24.0 && fabs(mu1Chi2) < 4.0 && MuTagHitsTrkSys > 7 ); //&& MuTagHitsMuSys > 3
|
| 638 |
if(!goodTrack)continue;
|
| 639 |
|
| 640 |
MuTagnSegTrkArb = muIter1->numberOfMatches(); // get number of chambers with matched segments
|
| 641 |
if (MuTagnSegTrkArb < 2) continue;
|
| 642 |
|
| 643 |
MuTagPx = muIter1->track()->px() ;
|
| 644 |
MuTagPy = muIter1->track()->py() ;
|
| 645 |
MuTagPz = muIter1->track()->pz() ;
|
| 646 |
|
| 647 |
MuTagE = muIter1->track()->p() ;
|
| 648 |
MuTagEta = muIter1->track()->eta() ;
|
| 649 |
MuTagPhi = muIter1->track()->phi() ;
|
| 650 |
|
| 651 |
Int_t MuTagcharge = muIter1->track()->charge() ;
|
| 652 |
|
| 653 |
// MuTagHitsPixSys = muIter1->track()->hitPattern().numberOfValidPixelHits();
|
| 654 |
// MuTagHitsRPCSys = muIter1->track()->hitPattern().numberOfValidMuonRPCHits();
|
| 655 |
|
| 656 |
/// check for HLT matching
|
| 657 |
minDRHLTDiMu=100.; minDRHLTAllSingleMu=100.;
|
| 658 |
minDRHLTMu->assign(HLTMuObjModuleNames->size(),100.);
|
| 659 |
Float_t PhiTemp1 = MuTagPhi<0?MuTagPhi + 2*M_PI:MuTagPhi;
|
| 660 |
if ( m_GotTrgObj ){
|
| 661 |
const trigger::TriggerObjectCollection & toc(handleTriggerEvent->getObjects());
|
| 662 |
for ( size_t ia = 0; ia < handleTriggerEvent->sizeFilters(); ++ ia) {
|
| 663 |
std::string fullname = handleTriggerEvent->filterTag(ia).encode();
|
| 664 |
std::string name;
|
| 665 |
size_t p = fullname.find_first_of(':');
|
| 666 |
if ( p != std::string::npos) name = fullname.substr(0, p);
|
| 667 |
else name = fullname;
|
| 668 |
if ( &toc ==0) continue;
|
| 669 |
const trigger::Keys & k = handleTriggerEvent->filterKeys(ia);
|
| 670 |
for (trigger::Keys::const_iterator ki = k.begin(); ki !=k.end(); ++ki ) {
|
| 671 |
double l3phi = toc[*ki].phi();
|
| 672 |
double l3eta = toc[*ki].eta();
|
| 673 |
|
| 674 |
if(l3phi < 0 )l3phi = l3phi + 2*M_PI;
|
| 675 |
|
| 676 |
Float_t deltaR_TTHLT = sqrt((l3eta - MuTagEta) * (l3eta - MuTagEta) + (l3phi - PhiTemp1) *(l3phi - PhiTemp1));
|
| 677 |
if ( deltaR_TTHLT > 0.4 ) continue;
|
| 678 |
if ( name == *HLTDiMuObjModuleName ) {
|
| 679 |
if ( deltaR_TTHLT < minDRHLTDiMu ){
|
| 680 |
minDRHLTDiMu = deltaR_TTHLT;
|
| 681 |
}
|
| 682 |
}
|
| 683 |
UChar_t idx=0;
|
| 684 |
for (vector<string>::const_iterator iter=HLTMuObjModuleNames->begin(); iter<HLTMuObjModuleNames->end();iter++,idx++ )
|
| 685 |
if ( name == *iter ) {
|
| 686 |
if ( deltaR_TTHLT < minDRHLTMu->at(idx) ){
|
| 687 |
minDRHLTMu->at(idx) = deltaR_TTHLT;
|
| 688 |
}
|
| 689 |
if ( deltaR_TTHLT < minDRHLTAllSingleMu ){
|
| 690 |
minDRHLTAllSingleMu = deltaR_TTHLT;
|
| 691 |
}
|
| 692 |
break;
|
| 693 |
}
|
| 694 |
}//end of loop ki
|
| 695 |
}//end of loop ia
|
| 696 |
}//end of if m_GotTrgObj
|
| 697 |
MuTagIsoR03Ratio = MuTagPt>0?muIter1->isolationR03().sumPt/MuTagPt:9999.;
|
| 698 |
MuTagIsoR05Ratio = MuTagPt>0?muIter1->isolationR05().sumPt/MuTagPt:9999.;
|
| 699 |
MuTagCaloL= muon::isGoodMuon(*muIter1,muon::TM2DCompatibilityLoose);
|
| 700 |
MuTagCaloT= muon::isGoodMuon(*muIter1,muon::TM2DCompatibilityTight);
|
| 701 |
|
| 702 |
/*------------------------Start getting the Monte Carlo Truth--------------------------*/
|
| 703 |
MuTagtracktruth_pt = -9999.;
|
| 704 |
MuTagtracktruth_p = -9999.;
|
| 705 |
MuTagtracktruth_id = -9999.;
|
| 706 |
MuTagtracktruth_type = 0;
|
| 707 |
MuTagtracktruth_isPileup=false;
|
| 708 |
MuTagtracktruth_thesamewith=-1;
|
| 709 |
if ( m_isMC) {
|
| 710 |
//TrackingParticles
|
| 711 |
MuTracksSimChains.push_back(* new vector< vector<Int_t> >() );
|
| 712 |
if ( TPCollectionH.isValid() ) {
|
| 713 |
//SimToReco Tracks Association
|
| 714 |
reco::TrackBaseRef trackBaseRef(muIter1->track());
|
| 715 |
if ( RecoToSimByHits.find(trackBaseRef) != RecoToSimByHits.end() ) {
|
| 716 |
pair<TrackingParticleRef, double> BestMatch=RecoToSimByHits[trackBaseRef].front();
|
| 717 |
//vector<pair<TrackingParticleRef, double> > TPCByChi2=RecoToSimByChi2[trk];
|
| 718 |
TrackingParticleRef tpr = BestMatch.first;
|
| 719 |
//Simulated Tracks
|
| 720 |
MuTagtracktruth_pt=tpr->pt();
|
| 721 |
MuTagtracktruth_p=tpr->p();
|
| 722 |
MuTagtracktruth_id=tpr->pdgId();
|
| 723 |
MuTagtracktruth_isPileup=GetDecayChains(tpr,HepGenEvent,MuTagtracktruth_type,MuTagtracktruth_thesamewith, MuTracksSimChains);
|
| 724 |
}
|
| 725 |
}//end of TPCollection.isValid
|
| 726 |
}//end of m_isMC
|
| 727 |
/*------------------------End of getting the Monte Carlo Truth--------------------------*/
|
| 728 |
|
| 729 |
CLHEP::Hep3Vector Mu1Vector(MuTagPx, MuTagPy, MuTagPz);
|
| 730 |
|
| 731 |
Float_t etaMuVec[4] ={-9999., -9999., -9999., -9999.};
|
| 732 |
Float_t phiMuVec[4] ={-9999., -9999., -9999., -9999.};
|
| 733 |
|
| 734 |
Float_t allCSC_mu = MEZ[0];
|
| 735 |
if(fabs(MuTagEta) < 1.5)allCSC_mu = MEZ[1];
|
| 736 |
if(MuTagEta < 0) allCSC_mu = -allCSC_mu;
|
| 737 |
|
| 738 |
reco::TrackRef trackMuRef = muIter1->track();
|
| 739 |
TrajectoryStateOnSurface tsos = surfExtrapTrkSam(trackMuRef, allCSC_mu);
|
| 740 |
if(tsos.isValid()){
|
| 741 |
etaMuVec[0] = tsos.globalPosition().eta();
|
| 742 |
phiMuVec[0] = tsos.globalPosition().phi();
|
| 743 |
}
|
| 744 |
allCSC_mu = MEZ[3];
|
| 745 |
if(MuTagEta < 0) allCSC_mu = -allCSC_mu;
|
| 746 |
tsos = surfExtrapTrkSam(trackMuRef, allCSC_mu);
|
| 747 |
if(tsos.isValid()){
|
| 748 |
etaMuVec[1] = tsos.globalPosition().eta();
|
| 749 |
phiMuVec[1] = tsos.globalPosition().phi();
|
| 750 |
}
|
| 751 |
|
| 752 |
allCSC_mu = MEZ[4];
|
| 753 |
if(MuTagEta < 0) allCSC_mu = -allCSC_mu;
|
| 754 |
tsos = surfExtrapTrkSam(trackMuRef, allCSC_mu);
|
| 755 |
if(tsos.isValid()){
|
| 756 |
etaMuVec[2] = tsos.globalPosition().eta();
|
| 757 |
phiMuVec[2] = tsos.globalPosition().phi();
|
| 758 |
}
|
| 759 |
|
| 760 |
allCSC_mu = MEZ[5];
|
| 761 |
if(MuTagEta < 0) allCSC_mu = -allCSC_mu;
|
| 762 |
tsos = surfExtrapTrkSam(trackMuRef, allCSC_mu);
|
| 763 |
if(tsos.isValid()){
|
| 764 |
etaMuVec[3] = tsos.globalPosition().eta();
|
| 765 |
phiMuVec[3] = tsos.globalPosition().phi();
|
| 766 |
}
|
| 767 |
|
| 768 |
Bool_t firsttrackmatchingtoMuTag=true;
|
| 769 |
|
| 770 |
for(reco::TrackCollection::const_iterator itTrack = gTracks->begin(); itTrack != gTracks->end(); itTrack++){
|
| 771 |
UInt_t itrk = itTrack - gTracks->begin();
|
| 772 |
if(itTrack->charge() == 0) continue;
|
| 773 |
if(itTrack->p() < 3.0) continue;
|
| 774 |
trackVeto_strict = trkVeto[itrk];
|
| 775 |
reco::TrackRef trackRef(gTracks, itrk );
|
| 776 |
|
| 777 |
if(itTrack->charge()*MuTagcharge != -1) continue;
|
| 778 |
|
| 779 |
tracks_eta = itTrack->eta();
|
| 780 |
|
| 781 |
if(fabs(tracks_eta) > 0.9 ) {
|
| 782 |
if(fabs(tracks_eta) < 1.2 ) myRegion = 2;
|
| 783 |
else myRegion = 3;
|
| 784 |
}
|
| 785 |
else myRegion = 1;
|
| 786 |
if (myRegion == 1) continue;//currently, only CSC
|
| 787 |
|
| 788 |
tracks_phi = itTrack->phi();
|
| 789 |
if(tracks_phi < 0 )tracks_phi = tracks_phi + 2*M_PI;
|
| 790 |
|
| 791 |
tracks_chi2 = -9999.;
|
| 792 |
if(itTrack->ndof() != 0)tracks_chi2 = itTrack->chi2()/itTrack->ndof();
|
| 793 |
tracks_charge = itTrack->charge();
|
| 794 |
|
| 795 |
tracks_dxy = itTrack->dxy();
|
| 796 |
tracks_dz = itTrack->dz();
|
| 797 |
if(beamExists){
|
| 798 |
tracks_dxy = itTrack->dxy(beamSpot.position());
|
| 799 |
tracks_dz = itTrack->dz(beamSpot.position());
|
| 800 |
}
|
| 801 |
|
| 802 |
MuProbenHitsTrkSys = itTrack->hitPattern().numberOfValidTrackerHits();
|
| 803 |
MuProbenHitsPixSys = itTrack->hitPattern().numberOfValidPixelHits();
|
| 804 |
|
| 805 |
goodTrack = ( fabs(itTrack->eta())< 2.4 && fabs(tracks_dz)< 24.0 &&
|
| 806 |
fabs(tracks_dxy)< 2.0 && tracks_chi2> 0.0 && tracks_chi2< 4.0 && MuProbenHitsTrkSys > 7 );
|
| 807 |
|
| 808 |
if (!goodTrack) continue;
|
| 809 |
if(tracks_eta > 0) CSCEndCapPlus = true;
|
| 810 |
else CSCEndCapPlus = false;
|
| 811 |
|
| 812 |
tracks_pt = itTrack->pt();
|
| 813 |
tracks_e = itTrack->p();
|
| 814 |
tracks_vx = itTrack->vx(); // x coordinate of the reference point on track
|
| 815 |
tracks_vy = itTrack->vy(); // y coordinate of the reference point on track
|
| 816 |
tracks_vz = itTrack->vz(); // z coordinate of the reference point on track
|
| 817 |
reco::IsoDeposit depTrk = muIsoExtractorTrack_->deposit(event, setup, *itTrack );
|
| 818 |
tracks_IsoR03Ratio = tracks_pt>0?depTrk.depositWithin(0.3)/tracks_pt:9999.;
|
| 819 |
tracks_IsoR05Ratio = tracks_pt>0?depTrk.depositWithin(0.5)/tracks_pt:9999.;
|
| 820 |
tracks_qoverp = itTrack->qoverp(); // q/p
|
| 821 |
tracks_lambda = itTrack->lambda();
|
| 822 |
tracks_recHitsSize= itTrack->recHitsSize();
|
| 823 |
tracks_qoverpError = itTrack->qoverpError();//TrackAlgorithm tAlg = itTrack->algo();
|
| 824 |
tracks_ptError = itTrack->ptError();/// error on Pt (set to 1000 TeV if charge==0 for safety)
|
| 825 |
tracks_thetaError = itTrack->thetaError();/// error on theta
|
| 826 |
tracks_lambdaError = itTrack->lambdaError(); /// error on lambda
|
| 827 |
tracks_etaError = itTrack->etaError(); tracks_phiError = itTrack->phiError(); tracks_dxyError = itTrack->dxyError();/// error on dxy
|
| 828 |
tracks_d0Error = itTrack->d0Error(); tracks_dszError = itTrack->dszError();tracks_dzError = itTrack->dzError();/// error on dz
|
| 829 |
tracks_numberOfValidHits = itTrack->numberOfValidHits(); ///unsigned short number of valid hits found
|
| 830 |
tracks_numberOfLostHits = itTrack->numberOfLostHits(); ///unsigned short number of cases where track crossed a layer without getting a hit.
|
| 831 |
if(tracks_phi < 0)tracks_phi = tracks_phi + 2*M_PI;
|
| 832 |
reco::MuonCollection::const_iterator matchedMu=matchTTwithMT(itTrack);
|
| 833 |
if ( matchedMu!=muons->end() ) {
|
| 834 |
tracks_isCaloMuTrk=matchedMu->isCaloMuon();
|
| 835 |
tracks_isTrackerMuTrk=matchedMu->isTrackerMuon();
|
| 836 |
}
|
| 837 |
else {
|
| 838 |
tracks_isCaloMuTrk=false;
|
| 839 |
tracks_isTrackerMuTrk=false;
|
| 840 |
}
|
| 841 |
|
| 842 |
Bool_t trQuality = (fabs(MuTagEta) < 2.4 && fabs(tracks_eta) > 0.9 && fabs(tracks_eta) < 2.4
|
| 843 |
&& tracks_etaError < 0.003 && tracks_phiError < 0.003
|
| 844 |
&& tracks_ptError/tracks_pt < 0.05 && tracks_numberOfValidHits >= 10); // cuts removed from the SkimDPG.C file and put here...
|
| 845 |
|
| 846 |
if(!trQuality)continue;
|
| 847 |
|
| 848 |
Float_t mMu = 0.1134289256;
|
| 849 |
invMass = pow( ( sqrt(pow(itTrack->p(),2)+ mMu*mMu) + sqrt(pow(muIter1->track()->p(),2)+ mMu*mMu) ) ,2 ) -
|
| 850 |
(
|
| 851 |
pow((itTrack->px() + muIter1->track()->px()),2) +
|
| 852 |
pow((itTrack->py() + muIter1->track()->py()),2) +
|
| 853 |
pow((itTrack->pz() + muIter1->track()->pz()),2)
|
| 854 |
);
|
| 855 |
|
| 856 |
if(invMass < 0) continue;
|
| 857 |
invMass = sqrt(invMass);
|
| 858 |
Bool_t gotMass = ((invMass > 2.5 && invMass < 3.6) || (invMass > 75.)) ;
|
| 859 |
if(!gotMass)continue;
|
| 860 |
|
| 861 |
/*------------------------Start getting the Monte Carlo Truth--------------------------*/
|
| 862 |
tracktruth_pt = -9999.;
|
| 863 |
tracktruth_e = -9999.;
|
| 864 |
tracktruth_p = -9999.;
|
| 865 |
tracktruth_id = -9999.;
|
| 866 |
tracktruth_type = 0;
|
| 867 |
tracktruth_isPileup=false;
|
| 868 |
tracktruth_thesamewith=-1;
|
| 869 |
|
| 870 |
if ( m_isMC) {
|
| 871 |
//TrackingParticles
|
| 872 |
if ( TPCollectionH.isValid() ) {
|
| 873 |
#ifdef jz_debug
|
| 874 |
// cout<<"track "<<i<<" has "<<tracks_numberOfValidHits << " hits" <<endl;
|
| 875 |
#endif
|
| 876 |
if (!firsttrackmatchingtoMuTag) MuTracksSimChains.push_back( MuTracksSimChains.back() );
|
| 877 |
TracksSimChains.push_back(* new vector< vector<Int_t> >() );
|
| 878 |
//SimToReco Tracks Association
|
| 879 |
reco::TrackBaseRef trackBaseRef(trackRef);
|
| 880 |
if ( RecoToSimByHits.find(trackBaseRef) != RecoToSimByHits.end() ) {
|
| 881 |
pair<TrackingParticleRef, double> BestMatch=RecoToSimByHits[trackBaseRef].front();
|
| 882 |
//vector<pair<TrackingParticleRef, double> > TPCByChi2=RecoToSimByChi2[trk];
|
| 883 |
TrackingParticleRef tpr = BestMatch.first;
|
| 884 |
//Simulated Tracks
|
| 885 |
tracktruth_pt=tpr->pt();
|
| 886 |
tracktruth_e=tpr->energy();
|
| 887 |
tracktruth_p=tpr->p();
|
| 888 |
tracktruth_id=tpr->pdgId();
|
| 889 |
tracktruth_isPileup=GetDecayChains(tpr,HepGenEvent,tracktruth_type,tracktruth_thesamewith,TracksSimChains);
|
| 890 |
}
|
| 891 |
}//end of TPCollection.isValid
|
| 892 |
}//end of m_isMC
|
| 893 |
/*------------------------End of getting the Monte Carlo Truth--------------------------*/
|
| 894 |
|
| 895 |
CLHEP::Hep3Vector trackVector(itTrack->px(), itTrack->py(), itTrack->pz());
|
| 896 |
//////////// Check if they have the same good vertex...
|
| 897 |
|
| 898 |
Bool_t gotVertexTrk1 = false, gotVertexTrk2 = false;
|
| 899 |
|
| 900 |
iSameVtx = false;
|
| 901 |
vtx_r = -9999.; vtx_z = -9999.;
|
| 902 |
vtx_rError = -9999.; vtx_zError = -9999.;
|
| 903 |
vtx_normChi2 = 100000.; vtx_size = -1;
|
| 904 |
|
| 905 |
for(reco::VertexCollection::const_iterator it=vertices.begin() ; it!=vertices.end() ; ++it) {
|
| 906 |
if(!it->isValid())continue;
|
| 907 |
if(it->isFake())continue;
|
| 908 |
for(std::vector<TrackBaseRef>::const_iterator itr=it->tracks_begin() ; itr!=it->tracks_end() ; ++itr) {
|
| 909 |
if ( & (*muIter1->track().get() ) == &( * const_cast<reco::Track *>(itr->get()) ) ) {
|
| 910 |
// cerr<<"Got MuTag trk matched"<<endl;
|
| 911 |
gotVertexTrk1 = true;
|
| 912 |
if (gotVertexTrk2) break;
|
| 913 |
}
|
| 914 |
if ( & (* itTrack) == & ( * (itr->get()) ) ) {
|
| 915 |
// cerr<<"Got TT trk matched"<<endl;
|
| 916 |
gotVertexTrk2 = true;
|
| 917 |
if (gotVertexTrk1) break;
|
| 918 |
}
|
| 919 |
}// tracks within vertex
|
| 920 |
|
| 921 |
if(gotVertexTrk1 && gotVertexTrk2){
|
| 922 |
vtx_size = it->tracksSize();
|
| 923 |
vtx_r = sqrt(pow(it->x(),2)+ pow(it->y(),2)); vtx_z = it->z();
|
| 924 |
if(it->x() != 0 && it->y() != 0) vtx_rError = vtx_r*sqrt(pow((it->xError()/it->x()),2) + pow((it->yError()/it->y()),2));
|
| 925 |
if(it->x() != 0 && it->y() == 0) vtx_rError = it->xError();
|
| 926 |
if(it->x() == 0 && it->y() != 0) vtx_rError = it->yError();
|
| 927 |
vtx_zError = it->zError();
|
| 928 |
if(it->ndof() != 0) vtx_normChi2 = it->chi2()/it->ndof();
|
| 929 |
iSameVtx = true;
|
| 930 |
}
|
| 931 |
}// for loop over vertices..
|
| 932 |
|
| 933 |
/////////// finished checking the vertex..
|
| 934 |
|
| 935 |
nTrkCountCSCSeg=0;
|
| 936 |
for(UChar_t j =0; j < 4; j++){
|
| 937 |
/*CSC Chamber Candidates in each station*/
|
| 938 |
CSCRg[j]= 0;
|
| 939 |
CSCChCand[j]= 0;
|
| 940 |
CSCChBad[j] = false;
|
| 941 |
/*Extrapolated Tracks on CSC Chamber Candidates in each station*/
|
| 942 |
CSCDyProjHVGap[j]=9999.;
|
| 943 |
CSCDyErrProjHVGap[j]=-9999.;
|
| 944 |
CSCProjEdgeDist[j]=-9999.;
|
| 945 |
CSCProjEdgeDistErr[j]=-9999.;
|
| 946 |
/*Segments characteristics*/
|
| 947 |
CSCSegxLc[j]=-9999.;
|
| 948 |
CSCSegyLc[j]=-9999.;
|
| 949 |
CSCSegxErrLc[j]=-9999.;
|
| 950 |
CSCSegyErrLc[j]=-9999.;
|
| 951 |
CSCSegChisqProb[j]=-9999.;
|
| 952 |
CSCdXdZTTSeg[j]=-9999.;
|
| 953 |
CSCdYdZTTSeg[j]=-9999.;
|
| 954 |
CSCnSegHits[j]=-9999;
|
| 955 |
/*Distance from the Extrapolated Tracks to CSC Segments, 9999. for no CSC segment found*/
|
| 956 |
CSCDxTTSeg[j] = -9999.;
|
| 957 |
CSCDxErrTTSeg[j] = -9999.;
|
| 958 |
CSCDyTTSeg[j] = -9999.;
|
| 959 |
CSCDyErrTTSeg[j] = -9999.;
|
| 960 |
CSCDxyTTSeg[j] = -9999.;
|
| 961 |
CSCDxyErrTTSeg[j] = -9999.;
|
| 962 |
/*LCT characteristics*/
|
| 963 |
CSCLCTxLc[j] = -9999.;
|
| 964 |
CSCLCTyLc[j] = -9999.;
|
| 965 |
CSCLCTbx[j] = -9999;
|
| 966 |
/*Distance from the Extrapolated Tracks to LCT, 9999. for no LCT found*/
|
| 967 |
CSCDxTTLCT[j] = -9999.;
|
| 968 |
CSCDxErrTTLCT[j] = -9999.;
|
| 969 |
CSCDyTTLCT[j] = -9999.;
|
| 970 |
CSCDyErrTTLCT[j] = -9999.;
|
| 971 |
CSCDxyTTLCT[j] = -9999.;
|
| 972 |
CSCDxyErrTTLCT[j] = -9999.;
|
| 973 |
/*DetlaR between the extrapolated tracker track on muon system and the tagged muon*/
|
| 974 |
dRTkMu[j] = -9999.;
|
| 975 |
/*Default decision of whether a segment or LCT is found*/
|
| 976 |
segSt[j]=false;
|
| 977 |
lctSt[j]=false;
|
| 978 |
}
|
| 979 |
|
| 980 |
for(UInt_t j =0; j < 6; j++) {
|
| 981 |
|
| 982 |
if( MEZ[j] == 0)continue;
|
| 983 |
|
| 984 |
Float_t zzPlaneME = MEZ[j];
|
| 985 |
if(!CSCEndCapPlus) zzPlaneME = -MEZ[j];
|
| 986 |
tsos = surfExtrapTrkSam(trackRef, zzPlaneME);
|
| 987 |
if (!tsos.isValid()) continue;
|
| 988 |
|
| 989 |
UChar_t st = j>2?j-2:0, ec=CSCEndCapPlus?1:2; // determine the station number...
|
| 990 |
Float_t trkEta = tsos.globalPosition().eta(), trkPhi = tsos.globalPosition().phi();
|
| 991 |
|
| 992 |
UChar_t rg = ringCandidate(st+1, trkEta, trkPhi);
|
| 993 |
if (!rg) continue;
|
| 994 |
if (!(
|
| 995 |
(j==0&&(rg==1||rg==4)) ||
|
| 996 |
(j==1&&rg==2) ||
|
| 997 |
(j==2&&rg==3) ||
|
| 998 |
(j>2)
|
| 999 |
)) continue;//check if it is runing on the right ring for a certain Z
|
| 1000 |
else CSCRg[st]=rg;
|
| 1001 |
|
| 1002 |
//for chamber overlap region
|
| 1003 |
/* This part of code is useless. Just directly extrapolate to chamber
|
| 1004 |
cerr<< "CHK"<<thisChamberCandidate(st+1, rg, trkPhi);
|
| 1005 |
if ( thisChamberCandidate(st+1, rg, trkPhi)%2 == 0 ) zzPlaneME = MEZEven[j];
|
| 1006 |
else zzPlaneME = MEZOdd[j];
|
| 1007 |
|
| 1008 |
if(!CSCEndCapPlus)zzPlaneME = -zzPlaneME;
|
| 1009 |
|
| 1010 |
tsos = surfExtrapTrkSam(trackRef, zzPlaneME);
|
| 1011 |
if (!tsos.isValid()) continue;
|
| 1012 |
trkEta = tsos.globalPosition().eta(); trkPhi = tsos.globalPosition().phi();
|
| 1013 |
rg = ringCandidate(st+1, trkEta, trkPhi);
|
| 1014 |
if ( rg==-9999 ) continue;
|
| 1015 |
*/
|
| 1016 |
|
| 1017 |
//Calculate the solid angle between the muon and this extrapolated track
|
| 1018 |
if(trkPhi < 0) trkPhi += 2*M_PI;
|
| 1019 |
/*dR between muon and track*/
|
| 1020 |
Float_t MPhi = phiMuVec[st];
|
| 1021 |
if (MPhi < 0 ) MPhi += 2*M_PI;
|
| 1022 |
Float_t TPhi = tsos.globalPosition().phi();
|
| 1023 |
if (TPhi < 0 ) TPhi += 2*M_PI;
|
| 1024 |
dRTkMu[st] = deltaR( etaMuVec[st], MPhi, tsos.globalPosition().eta() , TPhi);
|
| 1025 |
// printf("Mu(%f,%f),Trk(%f,%f)-->dR(%f)",etaMuVec[st], MPhi, tsos.globalPosition().eta() , TPhi,dRTkMu[st]);
|
| 1026 |
//Find the chamber candidate. If two chamber overlaps, simply divide them from middle. e.g. if chamber1 is from phi=1-3, chamber 2 is from phi=2-4. Then phi=1-2.5-->chamber 1; phi=2.5-4-->chamber 2. Since the distribution over phi is uniform, it won't bias the result.
|
| 1027 |
CSCChCand[st] = thisChamberCandidate(st+1, rg, trkPhi);
|
| 1028 |
CSCDetId Layer0Id=CSCDetId(ec, st+1, rg, CSCChCand[st], 0);//layer 0 is the mid point of the chamber. It is not a real layer.
|
| 1029 |
CSCChBad[st] = badChambers_->isInBadChamber( Layer0Id );
|
| 1030 |
//skip not-existing ME42 chambers
|
| 1031 |
if (CSCChBad[st]&&st==3&&rg==2) continue;
|
| 1032 |
#ifdef jz_debug
|
| 1033 |
if (CSCChBad[st]) cerr<<(CSCEndCapPlus?"ME+":"ME-")<<st+1<<"/"<<rg<<"/"<<CSCChBad[st]<<" is a dead chamber."<<endl;
|
| 1034 |
#endif
|
| 1035 |
/*Not necessary. Chamber search using phi is accurate enough.
|
| 1036 |
CSCDetId Layerid = CSCDetId( ec, st+1, rg, CSCChCand[st], 1 );
|
| 1037 |
vector<Float_t> EdgeAndDistToGap( GetEdgeAndDistToGap(trackRef,Layerid) );//values: 1-edge;2-err of edge;3-disttogap;4-err of dist to gap
|
| 1038 |
if (EdgeAndDistToGap[0]>0.||EdgeAndDistToGap[0]<-9990.) {
|
| 1039 |
cerr<<"ch"<<CSCChCand[st]<<", oldedge="<<EdgeAndDistToGap[0]<<"cm --> ";
|
| 1040 |
//try neighborhood chambers
|
| 1041 |
CSCDetId Layerid_plusone = CSCDetId( ec, st+1, rg, CSCChCand[st]+1, 1 );
|
| 1042 |
vector<Float_t> EdgeAndDistToGap_plusone( GetEdgeAndDistToGap(trackRef,Layerid_plusone) );
|
| 1043 |
CSCDetId Layerid_minusone = CSCDetId( ec, st+1, rg, CSCChCand[st]-1, 1 );
|
| 1044 |
vector<Float_t> EdgeAndDistToGap_minusone( GetEdgeAndDistToGap(trackRef,Layerid_minusone) );
|
| 1045 |
if (EdgeAndDistToGap_plusone[0]<EdgeAndDistToGap[0]&&EdgeAndDistToGap_plusone[0]>-9990.) {
|
| 1046 |
if (EdgeAndDistToGap_minusone[0]>EdgeAndDistToGap_plusone[0]) CSCChCand[st]+=1;
|
| 1047 |
else { if (EdgeAndDistToGap_minusone[0]>-9990.) CSCChCand[st]-=1;}
|
| 1048 |
}
|
| 1049 |
else { if (EdgeAndDistToGap_minusone[0]<EdgeAndDistToGap[0]&&EdgeAndDistToGap_minusone[0]>-9990.) CSCChCand[st]-=1;}
|
| 1050 |
cerr<<"ch"<<CSCChCand[st]<<", plusoneedge="<<EdgeAndDistToGap_plusone[0]<<"cm,"<<", minusoneedge="<<EdgeAndDistToGap_minusone[0]<<"cm."<<endl;
|
| 1051 |
}*/
|
| 1052 |
for (Int_t ly=1;ly<7;ly++) {
|
| 1053 |
CSCDetId Layerid = CSCDetId( ec, st+1, rg, CSCChCand[st], ly );
|
| 1054 |
vector<Float_t> EdgeAndDistToGap( GetEdgeAndDistToGap(trackRef,Layerid) );//values: 1-edge;2-err of edge;3-disttogap;4-err of dist to gap
|
| 1055 |
if (EdgeAndDistToGap[0]>CSCProjEdgeDist[st]) {
|
| 1056 |
CSCProjEdgeDist[st]=EdgeAndDistToGap[0];
|
| 1057 |
CSCProjEdgeDistErr[st]=EdgeAndDistToGap[1];
|
| 1058 |
}
|
| 1059 |
if (EdgeAndDistToGap[2]<CSCDyProjHVGap[st]) {
|
| 1060 |
CSCDyProjHVGap[st]=EdgeAndDistToGap[2];
|
| 1061 |
CSCDyErrProjHVGap[st]=EdgeAndDistToGap[3];
|
| 1062 |
}
|
| 1063 |
}
|
| 1064 |
//cerr<<"To Edge:"<<CSCProjEdgeDist[st]<<"; To HVGap:"<<CSCDyProjHVGap[st]<<endl;
|
| 1065 |
// Check the CSC segments in that region..
|
| 1066 |
CSCSegmentCollection::const_iterator cscSegOut;
|
| 1067 |
TrajectoryStateOnSurface *TrajToSeg = matchTTwithCSCSeg( trackRef, cscSegments, cscSegOut, Layer0Id);//update tsos with segment Z position, the segment z has little difference with layer 0 Z after the second or the third decimal place (in cm).
|
| 1068 |
if ( TrajToSeg!=NULL ) {//start saving matching information with segments
|
| 1069 |
/*
|
| 1070 |
const GeomDet* tmp_gdet=cscGeom->idToDet(cscSegOut->cscDetId());
|
| 1071 |
const CSCChamber* cscchamber = cscGeom->chamber(cscSegOut->cscDetId());
|
| 1072 |
cerr<<"CSCSEGDetID:"<<cscSegOut->cscDetId()<<" -- posZ: "<<tmp_gdet->surface().position().z()
|
| 1073 |
<<" -- CSCSeg posZ: "<< cscchamber->toGlobal(cscSegOut->localPosition()).z()<<endl;
|
| 1074 |
*/
|
| 1075 |
/* Save the chamber ID */
|
| 1076 |
CSCDetId id = (CSCDetId) cscSegOut->cscDetId();
|
| 1077 |
/* Save the segment postion and direction */
|
| 1078 |
LocalPoint localSegPos = (*cscSegOut).localPosition();
|
| 1079 |
//GlobalPoint globalSegPosition = cscchamber->toGlobal( localSegPos );
|
| 1080 |
CSCSegxLc[st] = localSegPos.x();
|
| 1081 |
CSCSegyLc[st] = localSegPos.y();
|
| 1082 |
LocalError localSegErr = (*cscSegOut).localPositionError();
|
| 1083 |
CSCSegxErrLc[st] = sqrt(localSegErr.xx());
|
| 1084 |
CSCSegyErrLc[st] = sqrt(localSegErr.yy());
|
| 1085 |
//cerr<<"CSCSegxLc:"<<localSegPos.x()<<"+-"<<CSCSegxErrLc[st]<<endl;
|
| 1086 |
//cerr<<"CSCSegyLc:"<<localSegPos.y()<<"+-"<<CSCSegyErrLc[st]<<endl;
|
| 1087 |
|
| 1088 |
/* Save the segment quality */
|
| 1089 |
CSCnSegHits[st] = cscSegOut->specificRecHits().size();
|
| 1090 |
Int_t nDOFCSC = 2*CSCnSegHits[st]-4;
|
| 1091 |
CSCSegChisqProb[st] = ChiSquaredProbability( double( (*cscSegOut).chi2() ), nDOFCSC );
|
| 1092 |
|
| 1093 |
/* Save the difference between the ex-tracker track and the segment */
|
| 1094 |
const GeomDet* gdet=cscGeom->idToDet(id);
|
| 1095 |
LocalPoint localpCSC = gdet->surface().toLocal(TrajToSeg->freeState()->position());
|
| 1096 |
CSCDxTTSeg[st] = CSCSegxLc[st] - localpCSC.x();
|
| 1097 |
CSCDyTTSeg[st] = CSCSegyLc[st] - localpCSC.y();
|
| 1098 |
CSCDxyTTSeg[st] = sqrt(pow(CSCDxTTSeg[st],2)+pow(CSCDyTTSeg[st],2));
|
| 1099 |
|
| 1100 |
LocalError localTTErr =TrajToSeg->localError().positionError();
|
| 1101 |
CSCDxErrTTSeg[st] = localSegErr.xx() + localTTErr.xx(); //the x error of the distance between segment and track squared
|
| 1102 |
CSCDyErrTTSeg[st] = localSegErr.yy() + localTTErr.yy(); //the y error of the distance between segment and track squared
|
| 1103 |
CSCDxyErrTTSeg[st] = sqrt(pow(CSCDxTTSeg[st],2)*CSCDxErrTTSeg[st]+pow(CSCDyTTSeg[st],2)*CSCDyErrTTSeg[st])/CSCDxyTTSeg[st];
|
| 1104 |
CSCDxErrTTSeg[st] = sqrt(CSCDxErrTTSeg[st]); CSCDyErrTTSeg[st] = sqrt(CSCDyErrTTSeg[st]);
|
| 1105 |
|
| 1106 |
LocalVector trackLocalDir = TrajToSeg->localDirection();
|
| 1107 |
LocalVector segDir = (*cscSegOut).localDirection();
|
| 1108 |
if ( trackLocalDir.z()!=0. && segDir.z()!=0.) {
|
| 1109 |
Float_t dxdz_trk = trackLocalDir.x()/trackLocalDir.z();
|
| 1110 |
Float_t dydz_trk = trackLocalDir.y()/trackLocalDir.z();
|
| 1111 |
Float_t dxdz_seg = segDir.x()/segDir.z();
|
| 1112 |
Float_t dydz_seg = segDir.y()/segDir.z();
|
| 1113 |
if(fabs(id.station()) == 3 || fabs(id.station()) == 4) dydz_seg = -dydz_seg;
|
| 1114 |
CSCdXdZTTSeg[st] = dxdz_trk - dxdz_seg;
|
| 1115 |
CSCdYdZTTSeg[st] = dydz_trk - dydz_seg;
|
| 1116 |
}
|
| 1117 |
nTrkCountCSCSeg++;
|
| 1118 |
delete TrajToSeg;
|
| 1119 |
}//end saving matching information with segments
|
| 1120 |
|
| 1121 |
////// Loop over MPC infromation to look for LCTs....
|
| 1122 |
trackVeto_isClosestToLCT=(BestTrkWithLCT[ec-1][st][rg-1][CSCChCand[st]-1]!=-1&&(signed)itrk!=BestTrkWithLCT[ec-1][st][rg-1][CSCChCand[st]-1]);
|
| 1123 |
CSCDetId Layer3id = CSCDetId( ec, st+1, rg, CSCChCand[st], 3 );//go to layer 3 that corresponds to the LCTPos
|
| 1124 |
const BoundPlane &Layer3Surface=cscGeom->idToDet(Layer3id)->surface();
|
| 1125 |
tsos=surfExtrapTrkSam(trackRef, Layer3Surface.position().z());
|
| 1126 |
if (tsos.isValid()) {//start matching with LCTs
|
| 1127 |
LocalPoint localL3pCSC = Layer3Surface.toLocal(tsos.freeState()->position());
|
| 1128 |
LocalPoint *LCTPos=matchTTwithLCTs( localL3pCSC.x(), localL3pCSC.y(), CSCEndCapPlus?1:2, st+1, CSCRg[st], CSCChCand[st], mpclcts, CSCDxyTTLCT[st], CSCLCTbx[st]);
|
| 1129 |
if (LCTPos!=NULL) {
|
| 1130 |
CSCLCTxLc[st]=LCTPos->x();
|
| 1131 |
CSCLCTyLc[st]=LCTPos->y();
|
| 1132 |
CSCDxTTLCT[st]=CSCLCTxLc[st]-localL3pCSC.x();
|
| 1133 |
CSCDyTTLCT[st]=CSCLCTyLc[st]-localL3pCSC.y();
|
| 1134 |
LocalError localTTErr =tsos.localError().positionError();
|
| 1135 |
CSCDxErrTTLCT[st] = sqrt(localTTErr.xx()); CSCDyErrTTLCT[st] = sqrt(localTTErr.yy());
|
| 1136 |
CSCDxyErrTTLCT[st] =sqrt(pow(CSCDxTTLCT[st],2)*localTTErr.xx() + pow(CSCDyTTLCT[st],2)*localTTErr.yy())/CSCDxyTTLCT[st];
|
| 1137 |
delete LCTPos;
|
| 1138 |
//cerr<<"segZ_TTy,LCTZ_TTy:"<<CSCSegyLc[st]-CSCDyTTSeg[st]<<","<<localL3pCSC.y()<<";"<<endl;
|
| 1139 |
}//end of if found LCT
|
| 1140 |
}//end matching with LCTs
|
| 1141 |
} // for loop for the stations -- j
|
| 1142 |
/*
|
| 1143 |
for (UChar_t st=0;st<4;st++) {
|
| 1144 |
lctSt[st] = ( CSCDxyTTLCT[st] >0. && CSCDxyTTLCT[st] < 40 )?1:0;
|
| 1145 |
segSt[st] = ( CSCDxyTTSeg[st] >0. && CSCDxyTTSeg[st] < 40)?1:0;
|
| 1146 |
}*/
|
| 1147 |
//cerr<<"eta="<<tracks_eta<<":ME1"<<CSCRg[0]<<":TTLCTx,TTLCTy,TTLCTxy:"<<CSCDxTTLCT[0]<<","<<CSCDyTTLCT[0]<<","<<CSCDxyTTLCT[0]<<endl;
|
| 1148 |
fractNtuple->Fill();
|
| 1149 |
firsttrackmatchingtoMuTag=false;
|
| 1150 |
} // loop over tracks...
|
| 1151 |
} // loop over muon tracks
|
| 1152 |
}
|
| 1153 |
|
| 1154 |
// ------------ method called once each job just before starting event loop ------------
|
| 1155 |
void TPTrackMuonSys::beginJob() {}
|
| 1156 |
// ------------ method called once each job just after ending the event loop ------------
|
| 1157 |
void TPTrackMuonSys::endJob() {
|
| 1158 |
theFile->Write();
|
| 1159 |
theFile->Close();
|
| 1160 |
std::cout<<"Events in "<<nEventsAnalyzed<<std::endl;
|
| 1161 |
}
|
| 1162 |
|
| 1163 |
//copied and modified from http://www.codeproject.com/KB/string/wildcmp.aspx
|
| 1164 |
//It does not belong to any class
|
| 1165 |
Bool_t wildcmp(const char *wild, const char *string) {
|
| 1166 |
// Written by Jack Handy - jakkhandy@hotmail.com
|
| 1167 |
const char *cp = NULL, *mp = NULL;
|
| 1168 |
|
| 1169 |
while ((*string) && (*wild != '*')) {
|
| 1170 |
if ((*wild != *string) && (*wild != '?')) {
|
| 1171 |
return false;
|
| 1172 |
}
|
| 1173 |
wild++;
|
| 1174 |
string++;
|
| 1175 |
}
|
| 1176 |
|
| 1177 |
while (*string) {
|
| 1178 |
if (*wild == '*') {
|
| 1179 |
if (!*++wild) {
|
| 1180 |
return true;
|
| 1181 |
}
|
| 1182 |
mp = wild;
|
| 1183 |
cp = string+1;
|
| 1184 |
} else if ((*wild == *string) || (*wild == '?')) {
|
| 1185 |
wild++;
|
| 1186 |
string++;
|
| 1187 |
} else {
|
| 1188 |
wild = mp;
|
| 1189 |
string = cp++;
|
| 1190 |
}
|
| 1191 |
}
|
| 1192 |
|
| 1193 |
while (*wild == '*') {
|
| 1194 |
wild++;
|
| 1195 |
}
|
| 1196 |
return !*wild;
|
| 1197 |
}
|
| 1198 |
|
| 1199 |
// ------------ method called in the beginning of each run ------------
|
| 1200 |
void TPTrackMuonSys::beginRun(const Run& r, const EventSetup& iSet)
|
| 1201 |
{
|
| 1202 |
run_number = r.runAuxiliary().run();
|
| 1203 |
iSet.get<CSCBadChambersRcd>().get(pBad);
|
| 1204 |
badChambers_=const_cast<CSCBadChambers*>(pBad.product());
|
| 1205 |
badChambersIndices=new vector<Int_t>( badChambers_->container() );
|
| 1206 |
|
| 1207 |
char plotname[100],plottitle[100];
|
| 1208 |
sprintf(plotname,"Run%d_BadChambers",run_number);
|
| 1209 |
sprintf(plottitle,"Known BadChambers in Run %d; chamber number",run_number);
|
| 1210 |
/* Draw a badchamber plot and save it to the ntuple*/
|
| 1211 |
TH2F *TH2F_BadChambers=new TH2F(plotname,plottitle,36,1,37,18,-9,9);
|
| 1212 |
const char *chambers[36] = {"01","02","03","04","05","06","07","08","09","10","11","12","13","14","15","16","17","18","19","20","21","22","23","24","25","26","27","28","29","30","31","32","33","34","35","36"};
|
| 1213 |
const char *rings[18] = {"ME-42","ME-41","ME-32","ME-31","ME-22","ME-21","ME-13","ME-12","ME-11","ME+11","ME+12","ME+13","ME+21","ME+22","ME+31","ME+32","ME+41","ME+42"};
|
| 1214 |
for (UChar_t i=0;i<36;i++)
|
| 1215 |
TH2F_BadChambers->GetXaxis()->SetBinLabel(i+1,chambers[i]);
|
| 1216 |
for (UChar_t i=0;i<18;i++)
|
| 1217 |
TH2F_BadChambers->GetYaxis()->SetBinLabel(i+1,rings[i]);
|
| 1218 |
for( Short_t indexc = 1; indexc<=540; ++indexc ) {// chamber indices are in range 1-468 (CSCs 2008) or 469-540 (ME42)
|
| 1219 |
CSCDetId id = CSCIndexer().detIdFromChamberIndex( indexc );
|
| 1220 |
if ( !badChambers_->isInBadChamber( id ) ) continue;
|
| 1221 |
UChar_t ring=id.station()*10+id.ring();
|
| 1222 |
Float_t fillY;
|
| 1223 |
switch ( ring )
|
| 1224 |
{
|
| 1225 |
case 14:
|
| 1226 |
fillY=0.5;
|
| 1227 |
break;
|
| 1228 |
case 11:
|
| 1229 |
fillY=0.5;
|
| 1230 |
break;
|
| 1231 |
case 12:
|
| 1232 |
fillY=1.5;
|
| 1233 |
break;
|
| 1234 |
case 13:
|
| 1235 |
fillY=2.5;
|
| 1236 |
break;
|
| 1237 |
case 21:
|
| 1238 |
fillY=3.5;
|
| 1239 |
break;
|
| 1240 |
case 22:
|
| 1241 |
fillY=4.5;
|
| 1242 |
break;
|
| 1243 |
case 31:
|
| 1244 |
fillY=5.5;
|
| 1245 |
break;
|
| 1246 |
case 32:
|
| 1247 |
fillY=6.5;
|
| 1248 |
break;
|
| 1249 |
case 41:
|
| 1250 |
fillY=7.5;
|
| 1251 |
break;
|
| 1252 |
case 42:
|
| 1253 |
fillY=8.5;
|
| 1254 |
break;
|
| 1255 |
default:
|
| 1256 |
printf("Unexpected ring number: %d",ring);
|
| 1257 |
fillY=9.5;
|
| 1258 |
}
|
| 1259 |
if (id.endcap()==2) fillY*=-1;
|
| 1260 |
TH2F_BadChambers->Fill(id.chamber()+0.5,fillY);
|
| 1261 |
#ifdef jz_debug
|
| 1262 |
cerr<<(id.endcap()==1?"ME+":"ME-")<<Int_t(ring)<<"/"<<id.chamber()<<"("<<indexc<<")";
|
| 1263 |
#endif
|
| 1264 |
}
|
| 1265 |
// TCanvas *BadChambersView=new TCanvas("badch","badch",1200,1000);
|
| 1266 |
// BadChambersView->SetGrid();
|
| 1267 |
TH2F_BadChambers->SetStats(0);
|
| 1268 |
TH2F_BadChambers->SetMinimum(0);
|
| 1269 |
TH2F_BadChambers->SetMaximum(0.8);
|
| 1270 |
TH2F_BadChambers->Draw("colz");
|
| 1271 |
TH2F_BadChambers->SetLabelSize(0.035,"X");
|
| 1272 |
#ifndef GetCSCHitsBefore500
|
| 1273 |
theFile->cd();
|
| 1274 |
#endif
|
| 1275 |
TH2F_BadChambers->Write();
|
| 1276 |
// BadChambersView->Write("BadChambersPlot");
|
| 1277 |
/* End of drawing the badchamber plot*/
|
| 1278 |
|
| 1279 |
Bool_t isConfigChanged = false;
|
| 1280 |
m_HLTMuTrgBit.clear();
|
| 1281 |
HLTMuNames->clear();
|
| 1282 |
HLTMuObjModuleNames->clear();
|
| 1283 |
m_HLTDiMuTrgBit=-1;
|
| 1284 |
HLTDiMuName=new string();
|
| 1285 |
HLTDiMuObjModuleName=new string();
|
| 1286 |
if ( hltConfigProvider_.init( r, iSet, m_hlt.process() , isConfigChanged ) ) {
|
| 1287 |
HLTTableName=new string( hltConfigProvider_.tableName() );
|
| 1288 |
const vector<string> & HLTNamesSet_=hltConfigProvider_.triggerNames();
|
| 1289 |
UInt_t idx=0;
|
| 1290 |
for ( vector<string>::const_iterator itertable=HLTNamesSet_.begin();itertable != HLTNamesSet_.end();itertable++,idx++ ) {
|
| 1291 |
#ifdef jz_debug
|
| 1292 |
cout<<endl<<idx<<":"<<*itertable;
|
| 1293 |
#endif
|
| 1294 |
Bool_t selected=false;
|
| 1295 |
for ( vector<string>::const_iterator iter=m_HLTMuTrgNames.begin();iter!=m_HLTMuTrgNames.end();iter++ ) {
|
| 1296 |
if ( wildcmp(iter->c_str(),itertable->c_str()) ) {
|
| 1297 |
selected=true;
|
| 1298 |
break;
|
| 1299 |
}
|
| 1300 |
}
|
| 1301 |
if (selected) {
|
| 1302 |
m_HLTMuTrgBit.push_back(idx);
|
| 1303 |
HLTMuNames->push_back(*itertable);
|
| 1304 |
HLTMuObjModuleNames->push_back( *(hltConfigProvider_.moduleLabels(idx).end()-2) );
|
| 1305 |
#ifdef jz_debug
|
| 1306 |
cout<<", saved for singlemu";
|
| 1307 |
#endif
|
| 1308 |
}
|
| 1309 |
if ( m_HLTDiMuTrgBit<0 && wildcmp(m_HLTDiMuTrgName.c_str(),itertable->c_str()) ) {
|
| 1310 |
m_HLTDiMuTrgBit=idx;
|
| 1311 |
HLTDiMuName->assign(*itertable);
|
| 1312 |
HLTDiMuObjModuleName->assign( *(hltConfigProvider_.moduleLabels(idx).end()-2) );
|
| 1313 |
#ifdef jz_debug
|
| 1314 |
cout<<", saved for doublemu";
|
| 1315 |
#endif
|
| 1316 |
}
|
| 1317 |
}//end of iter HLT table
|
| 1318 |
RunInfo->Fill();
|
| 1319 |
} else {
|
| 1320 |
// if init returns FALSE, initialisation has NOT succeeded, which indicates a problem
|
| 1321 |
// with the file and/or code and needs to be investigated!
|
| 1322 |
LogWarning("DataLost") << " HLT config extraction failure with process name " <<m_hlt.process();
|
| 1323 |
// In this case, all access methods will return empty values!
|
| 1324 |
}
|
| 1325 |
}
|
| 1326 |
|
| 1327 |
vector<Float_t> TPTrackMuonSys::GetEdgeAndDistToGap(reco::TrackRef trackRef, CSCDetId & detid) {
|
| 1328 |
vector<Float_t> result(4,9999.);
|
| 1329 |
result[3]=-9999;
|
| 1330 |
const GeomDet* gdet=cscGeom->idToDet( detid );
|
| 1331 |
TrajectoryStateOnSurface tsos=surfExtrapTrkSam(trackRef, gdet->surface().position().z());
|
| 1332 |
if (!tsos.isValid()) return result;
|
| 1333 |
LocalPoint localTTPos = gdet->surface().toLocal(tsos.freeState()->position());
|
| 1334 |
const CSCWireTopology* wireTopology = cscGeom->layer(detid)->geometry()->wireTopology();
|
| 1335 |
Float_t wideWidth = wireTopology->wideWidthOfPlane();
|
| 1336 |
Float_t narrowWidth = wireTopology->narrowWidthOfPlane();
|
| 1337 |
Float_t length = wireTopology->lengthOfPlane();
|
| 1338 |
// If slanted, there is no y offset between local origin and symmetry center of wire plane
|
| 1339 |
Float_t yOfFirstWire = fabs(wireTopology->wireAngle())>1.E-06 ? -0.5*length : wireTopology->yOfWire(1);
|
| 1340 |
// y offset between local origin and symmetry center of wire plane
|
| 1341 |
Float_t yCOWPOffset = yOfFirstWire+0.5*length;
|
| 1342 |
// tangent of the incline angle from inside the trapezoid
|
| 1343 |
Float_t tangent = (wideWidth-narrowWidth)/(2.*length);
|
| 1344 |
// y position wrt bottom of trapezoid
|
| 1345 |
Float_t yPrime = localTTPos.y()+fabs(yOfFirstWire);
|
| 1346 |
// half trapezoid width at y' is 0.5 * narrowWidth + x side of triangle with the above tangent and side y'
|
| 1347 |
Float_t halfWidthAtYPrime = 0.5*narrowWidth+yPrime*tangent;
|
| 1348 |
// x offset between local origin and symmetry center of wire plane is zero
|
| 1349 |
// x offset of ME11s is also zero. x center of wire groups is not at zero, because it is not parallel to x. The wire groups of ME11s have a complex geometry, see the code in m_debug.
|
| 1350 |
Float_t edgex = fabs(localTTPos.x()) - halfWidthAtYPrime;
|
| 1351 |
Float_t edgey = fabs(localTTPos.y()-yCOWPOffset) - 0.5*length;
|
| 1352 |
LocalError localTTErr = tsos.localError().positionError();
|
| 1353 |
// cerr<<"x err = "<<sqrt(localTTErr.xx())<<" ; y err = "<<sqrt(localTTErr.yy())<<endl;
|
| 1354 |
if ( edgex > edgey ) {
|
| 1355 |
result[0] = edgex;
|
| 1356 |
result[1] = sqrt( localTTErr.xx() );
|
| 1357 |
//result[1] = sqrt(tsos.cartesianError().position().cxx());
|
| 1358 |
}
|
| 1359 |
else {
|
| 1360 |
result[0] = edgey;
|
| 1361 |
result[1] = sqrt( localTTErr.yy() );
|
| 1362 |
//result[1] = sqrt(tsos.cartesianError().position().cyy());
|
| 1363 |
}
|
| 1364 |
result[2]=YDistToHVDeadZone(localTTPos.y(), detid.station()*10+detid.ring());
|
| 1365 |
result[3]=sqrt( localTTErr.yy() );
|
| 1366 |
return result;//return values: 1-edge;2-err of edge;3-disttogap;4-err of dist to gap
|
| 1367 |
}
|
| 1368 |
|
| 1369 |
////////////////////////////////////////////
|
| 1370 |
reco::MuonCollection::const_iterator TPTrackMuonSys::matchTTwithMT(reco::TrackCollection::const_iterator &itrack)
|
| 1371 |
{
|
| 1372 |
for (reco::MuonCollection::const_iterator muIter = muons->begin(); muIter != muons->end(); muIter++) {
|
| 1373 |
if (muIter->combinedMuon().isNonnull() )
|
| 1374 |
if ( &( * muIter->combinedMuon().get() ) == & (*itrack) ) return muIter;
|
| 1375 |
if (muIter->track().isNonnull() )
|
| 1376 |
if ( &(* muIter->track().get() ) == & (*itrack) ) return muIter;
|
| 1377 |
if (muIter->standAloneMuon().isNonnull() )
|
| 1378 |
if ( &( *muIter->standAloneMuon().get() ) == & (*itrack) ) return muIter;
|
| 1379 |
}
|
| 1380 |
return muons->end();
|
| 1381 |
}
|
| 1382 |
|
| 1383 |
////// Match TT with RPCEndCapHit
|
| 1384 |
Bool_t TPTrackMuonSys::matchTTwithCSCRecHit(Bool_t trackDir,
|
| 1385 |
Int_t j,
|
| 1386 |
reco::TrackRef trackRef,
|
| 1387 |
edm::Handle<CSCRecHit2DCollection> recHits,
|
| 1388 |
//std::vector<CSCRecHit2DCollection> recHitOut,
|
| 1389 |
std::vector<CSCRecHit2D> & recHitOut,
|
| 1390 |
std::vector<Int_t >& deltaRecHitX,
|
| 1391 |
std::vector<Int_t >& deltaRecHitY)
|
| 1392 |
{
|
| 1393 |
|
| 1394 |
Bool_t aMatch = false;
|
| 1395 |
CSCRecHit2D tmp[6];
|
| 1396 |
|
| 1397 |
Float_t minDX[6]={999., 999., 999., 999., 999., 999.};
|
| 1398 |
Float_t minDY[6]={999., 999., 999., 999., 999., 999.};
|
| 1399 |
Float_t minDR[6]={999., 999., 999., 999., 999., 999.};
|
| 1400 |
recHitOut.clear();
|
| 1401 |
|
| 1402 |
Float_t rCut = 50.;
|
| 1403 |
|
| 1404 |
for (CSCRecHit2DCollection::const_iterator recIt = recHits->begin(); recIt != recHits->end(); recIt++) {
|
| 1405 |
CSCDetId id = (CSCDetId)(*recIt).cscDetId();
|
| 1406 |
|
| 1407 |
if((trackDir && id.endcap() == 2) || (!trackDir && id.endcap() == 1)) continue;
|
| 1408 |
if(
|
| 1409 |
(j == 0 && id.station() == 1 && ((id.ring() == 1) || (id.ring() == 4 ))) || // ME1/1,4
|
| 1410 |
(j == 1 && id.station() == 1 && (id.ring() == 2)) || // ME1/2
|
| 1411 |
(j == 2 && id.station() == 1 && (id.ring() == 3)) || // ME1/3
|
| 1412 |
(j == 3 && id.station() == 2) || // ME2/1,2
|
| 1413 |
(j == 4 && id.station() == 3) || // ME3/1,2
|
| 1414 |
(j == 5 && id.station() == 4) // ME4/1
|
| 1415 |
){
|
| 1416 |
|
| 1417 |
// Get pointer to the layer:
|
| 1418 |
const CSCLayer* csclayer = cscGeom->layer( id );
|
| 1419 |
if (!csclayer)continue;
|
| 1420 |
|
| 1421 |
// Transform hit position from local chamber geometry to global CMS geom
|
| 1422 |
LocalPoint rhitlocal = (*recIt).localPosition();
|
| 1423 |
GlobalPoint rhitglobal = csclayer->toGlobal(rhitlocal);
|
| 1424 |
Float_t RecHit_grecx = rhitglobal.x();
|
| 1425 |
Float_t RecHit_grecy = rhitglobal.y();
|
| 1426 |
TrajectoryStateOnSurface tsos = surfExtrapTrkSam(trackRef, rhitglobal.z());
|
| 1427 |
if (!tsos.isValid()) continue;
|
| 1428 |
Float_t cscdeltaX = RecHit_grecx - tsos.globalPosition().x();
|
| 1429 |
Float_t cscdeltaY = RecHit_grecy - tsos.globalPosition().y();
|
| 1430 |
Float_t cscdeltaR = sqrt(pow(cscdeltaX,2)+pow(cscdeltaY,2));
|
| 1431 |
|
| 1432 |
Int_t mlayer = id.layer() - 1;
|
| 1433 |
if(fabs(cscdeltaR) < fabs(minDR[mlayer]) && fabs(cscdeltaR) < rCut){
|
| 1434 |
aMatch = true;
|
| 1435 |
minDX[mlayer] = cscdeltaX;
|
| 1436 |
minDY[mlayer] = cscdeltaY;
|
| 1437 |
minDR[mlayer] = cscdeltaR;
|
| 1438 |
tmp[mlayer] = (*recIt);
|
| 1439 |
}
|
| 1440 |
}
|
| 1441 |
}
|
| 1442 |
|
| 1443 |
for (Int_t i = 0; i < 6; i++) {
|
| 1444 |
Float_t cscdeltaR = sqrt(pow(minDX[i],2)+pow(minDY[i],2));
|
| 1445 |
if(cscdeltaR < rCut) recHitOut.push_back(tmp[i]);
|
| 1446 |
}
|
| 1447 |
for (Int_t i = 0; i < 6; i++) {
|
| 1448 |
deltaRecHitX.push_back(minDX[i]);
|
| 1449 |
deltaRecHitY.push_back(minDY[i]);
|
| 1450 |
}
|
| 1451 |
return aMatch;
|
| 1452 |
}
|
| 1453 |
|
| 1454 |
////// Match TT with RPCEndCapHit
|
| 1455 |
/*
|
| 1456 |
Bool_t TPTrackMuonSys::matchTTwithRPCEChit(Bool_t trackDir,
|
| 1457 |
Int_t j,
|
| 1458 |
reco::TrackRef trackRef,
|
| 1459 |
edm::Handle<RPCRecHitCollection> rpcRecHits,
|
| 1460 |
RPCRecHitCollection::const_iterator &rpcHitOut)
|
| 1461 |
// Float_t deltaRPCX, Float_t deltaRPCY)
|
| 1462 |
{
|
| 1463 |
Bool_t aMatch = false;
|
| 1464 |
|
| 1465 |
//Float_t deltaRPCX = 999.0, deltaRPCY = 999.0;
|
| 1466 |
Float_t deltaRPCR = 999.0;// deltaRPCPhi = 999.0;
|
| 1467 |
|
| 1468 |
for (RPCRecHitCollection::const_iterator rpcIt = rpcRecHits->begin(); rpcIt != rpcRecHits->end(); rpcIt++) {
|
| 1469 |
RPCDetId id = (RPCDetId)(*rpcIt).rpcId();
|
| 1470 |
const GeomDet* gdet=rpcGeo->idToDet(id);
|
| 1471 |
const BoundPlane & surface = gdet->surface();
|
| 1472 |
if(id.region() == 0)continue;
|
| 1473 |
if((trackDir && id.region() == -1) || (!trackDir && id.region() == 1) ) continue;
|
| 1474 |
if(
|
| 1475 |
(j >= 1 && j <= 2 && id.station() == 1) ||
|
| 1476 |
(j == 3 && id.station() == 2) ||
|
| 1477 |
(j == 4 && id.station() == 3)
|
| 1478 |
){
|
| 1479 |
|
| 1480 |
LocalPoint localPos = (*rpcIt).localPosition();
|
| 1481 |
GlobalPoint globalPosition=surface.toGlobal(localPos);
|
| 1482 |
TrajectoryStateOnSurface tsos = surfExtrapTrkSam(trackRef, globalPosition.z());
|
| 1483 |
if (!tsos.isValid()) continue;
|
| 1484 |
|
| 1485 |
Float_t RPCEdeltaXMin = globalPosition.x() - tsos.globalPosition().x();
|
| 1486 |
Float_t RPCEdeltaYMin = globalPosition.y() - tsos.globalPosition().y();
|
| 1487 |
//Float_t RPCEdeltaRMin = sqrt(pow(RPCEdeltaXMin,2)+pow(RPCEdeltaYMin,2));
|
| 1488 |
|
| 1489 |
Float_t RPCEPhiProj = tsos.globalPosition().phi();
|
| 1490 |
Float_t RPCEPhiPos = globalPosition.phi();
|
| 1491 |
|
| 1492 |
////// New for RPCEC
|
| 1493 |
if(RPCEPhiProj < 0) RPCEPhiProj = RPCEPhiProj + 2*M_PI;
|
| 1494 |
if(RPCEPhiPos < 0) RPCEPhiPos = RPCEPhiPos + 2*M_PI;
|
| 1495 |
|
| 1496 |
// Float_t RPCEdeltaPhiMin = fabs(RPCEPhiProj - RPCEPhiPos);//sqrt(pow(RPCEdeltaXMin,2)+pow(RPCEdeltaYMin,2));
|
| 1497 |
|
| 1498 |
Float_t RPCEdeltaRMin = sqrt(pow(RPCEdeltaXMin,2)+pow(RPCEdeltaYMin,2));
|
| 1499 |
if(fabs(RPCEdeltaRMin) < fabs(deltaRPCR) && fabs(RPCEdeltaRMin) < 50. ){
|
| 1500 |
aMatch = true;
|
| 1501 |
// deltaRPCPhi = RPCEdeltaPhiMin;
|
| 1502 |
deltaRPCR = RPCEdeltaRMin;
|
| 1503 |
rpcHitOut = rpcIt;
|
| 1504 |
}
|
| 1505 |
}
|
| 1506 |
} // loop over RPCE rechits
|
| 1507 |
|
| 1508 |
return aMatch;
|
| 1509 |
}
|
| 1510 |
*/
|
| 1511 |
////////////// Get the matching with LCTs...
|
| 1512 |
LocalPoint * TPTrackMuonSys::matchTTwithLCTs(Float_t xPos, Float_t yPos, UChar_t ec, UChar_t st, UChar_t &rg, UChar_t cham,
|
| 1513 |
edm::Handle<CSCCorrelatedLCTDigiCollection> mpclcts, Float_t &dRTrkLCT, Int_t &lctBX ) {
|
| 1514 |
LocalPoint *interSect=NULL;
|
| 1515 |
|
| 1516 |
for (CSCCorrelatedLCTDigiCollection::DigiRangeIterator detMPCUnitIt = mpclcts->begin();
|
| 1517 |
detMPCUnitIt != mpclcts->end(); detMPCUnitIt++) {
|
| 1518 |
CSCDetId id = (*detMPCUnitIt).first;
|
| 1519 |
|
| 1520 |
if(ec != id.endcap())continue;
|
| 1521 |
if(st != id.station())continue;
|
| 1522 |
if(cham != id.chamber())continue;
|
| 1523 |
|
| 1524 |
Bool_t ed1 = (st == 1) && ((rg == 1 || rg == 4) && (id.ring() == 1 || id.ring() == 4));
|
| 1525 |
Bool_t ed2 = (st == 1) && ((rg == 2 && id.ring() == 2) || (rg == 3 && id.ring() == 3));
|
| 1526 |
Bool_t ed3 = (st != 1) && (rg == id.ring());
|
| 1527 |
if ( !(ed1 || ed2 || ed3) ) continue;
|
| 1528 |
const CSCCorrelatedLCTDigiCollection::Range& MPCrange = (*detMPCUnitIt).second;
|
| 1529 |
for (CSCCorrelatedLCTDigiCollection::const_iterator mpcIt = MPCrange.first; mpcIt != MPCrange.second; mpcIt++) {
|
| 1530 |
Bool_t lct_valid = (*mpcIt).isValid();
|
| 1531 |
if(!lct_valid)continue;
|
| 1532 |
//In CSC offline/hlt software in general, is COUNT FROM ONE, such as CSCGeometry.
|
| 1533 |
//However, the LCT software counts from zero. strip_id is from 0 to 159. wireGroup_id is from 0 to 31(ME13),47(ME11),63(ME1234/2),95(ME31,41),111(ME21). So here we need to plus one.
|
| 1534 |
UChar_t wireGroup_id = (*mpcIt).getKeyWG()+1;
|
| 1535 |
UChar_t strip_id=(*mpcIt).getStrip()/2+1;
|
| 1536 |
//if ( id.ring()==4 ) cerr<<"Alert id.ring()==4"<<endl;
|
| 1537 |
Bool_t me11=(st == 1) && (id.ring() == 1 || id.ring() == 4),
|
| 1538 |
me11a = me11 && strip_id>64;
|
| 1539 |
//me11b = me11 && strip_id<=64;
|
| 1540 |
//http://cmssdt.cern.ch/SDT/lxr/source/CondFormats/CSCObjects/src/CSCChannelTranslator.cc
|
| 1541 |
//Translate a raw strip channel in range 1-80, iraw, into corresponding geometry-oriented channel in which increasing
|
| 1542 |
//channel number <-> strip number increasing with +ve local x.
|
| 1543 |
//However, the LCT firmwire has already done the number flipping but it numbers ME11B strips from 0 to 63, ME11A strips from 64-79. (confirmed by Slava)
|
| 1544 |
//That is the reason we commented out the following number flipping codes.
|
| 1545 |
if ( me11a ) {
|
| 1546 |
strip_id-=64;
|
| 1547 |
// The CSCCorrelatedLCTDigi DetId does NOT distinguish ME11A and B. All of the DetIDs are labelled as ME11B (all ME11, none ME14)
|
| 1548 |
// However, stripWireGroupIntersection must know that since ME11A has 48 strips and ME11B has 64.
|
| 1549 |
id=CSCDetId(ec, 1, 4, cham, 3);
|
| 1550 |
//if ( id.endcap()==1 ) strip_id=17-strip_id;
|
| 1551 |
}
|
| 1552 |
//if ( me11b && id.endcap()!=1 ) strip_id=65-strip_id;
|
| 1553 |
const CSCLayerGeometry *layerGeom = cscGeom->chamber(id)->layer (3)->geometry ();
|
| 1554 |
for(UChar_t ii = 0; ii < 3; ii++){
|
| 1555 |
// if ( strip_id>64 ) LogWarning("Strip_id") << "Got "<<strip_id<<", but there are "<< Nstrips <<" strips in total." <<m_hlt.process();
|
| 1556 |
LocalPoint interSect_ = layerGeom->stripWireGroupIntersection(strip_id, wireGroup_id);
|
| 1557 |
// printf( "ME%d/%d: %.1f/%d, %d/%d: xLCT-xTT=%.2f-%.2f; yLCT-yTT=%.2f-%.2f \n",st,id.ring(),strip_id,Nstrips,wireGroup_id,layerGeom->numberOfWireGroups(),interSect_.x(),xPos,interSect_.y(),yPos);
|
| 1558 |
Float_t DeltaR_ = sqrt(pow((interSect_.x()-xPos),2) + pow((interSect_.y()-yPos),2));
|
| 1559 |
if( DeltaR_ < fabs(dRTrkLCT) ) {
|
| 1560 |
delete interSect;
|
| 1561 |
interSect=new LocalPoint(interSect_);
|
| 1562 |
dRTrkLCT = DeltaR_ ;
|
| 1563 |
lctBX = (*mpcIt).getBX();
|
| 1564 |
if (me11a) rg=4;
|
| 1565 |
else rg=id.ring();
|
| 1566 |
//cout << "1: BX = " << (*mpcIt).getBX() << " BX0 = " << (*mpcIt).getBX0() << std::endl;
|
| 1567 |
} // for the matching if statement...
|
| 1568 |
if (me11a) strip_id+=16;
|
| 1569 |
else break;
|
| 1570 |
}// end iteration over of ME11A triplet
|
| 1571 |
}// end iteration over lcts_mpc_data in one chamber
|
| 1572 |
}// end iteration over lcts_mpc_data
|
| 1573 |
return interSect;
|
| 1574 |
}
|
| 1575 |
|
| 1576 |
inline Float_t TPTrackMuonSys::TrajectoryDistToSeg( TrajectoryStateOnSurface *TrajSuf, CSCSegmentCollection::const_iterator segIt) {
|
| 1577 |
if (!TrajSuf->isValid()) return 9999.;
|
| 1578 |
const GeomDet* gdet=cscGeom->idToDet( (CSCDetId)(*segIt).cscDetId() );
|
| 1579 |
LocalPoint localTTPos = gdet->surface().toLocal(TrajSuf->freeState()->position());
|
| 1580 |
LocalPoint localSegPos = (*segIt).localPosition();
|
| 1581 |
Float_t CSCdeltaX = localSegPos.x() - localTTPos.x();
|
| 1582 |
Float_t CSCdeltaY = localSegPos.y() - localTTPos.y();
|
| 1583 |
return sqrt(pow(CSCdeltaX,2)+pow(CSCdeltaY,2));
|
| 1584 |
}
|
| 1585 |
|
| 1586 |
TrajectoryStateOnSurface *TPTrackMuonSys::matchTTwithCSCSeg( reco::TrackRef trackRef, edm::Handle<CSCSegmentCollection> cscSegments,
|
| 1587 |
CSCSegmentCollection::const_iterator &cscSegOut, CSCDetId & idCSC ) {
|
| 1588 |
TrajectoryStateOnSurface *TrajSuf=NULL;
|
| 1589 |
Float_t deltaCSCR = 9999.;
|
| 1590 |
for(CSCSegmentCollection::const_iterator segIt=cscSegments->begin(); segIt != cscSegments->end(); segIt++) {
|
| 1591 |
CSCDetId id = (CSCDetId)(*segIt).cscDetId();
|
| 1592 |
if(idCSC.endcap() != id.endcap())continue;
|
| 1593 |
if(idCSC.station() != id.station())continue;
|
| 1594 |
if(idCSC.chamber() != id.chamber())continue;
|
| 1595 |
|
| 1596 |
Bool_t ed1 = (idCSC.station() == 1) && ((idCSC.ring() == 1 || idCSC.ring() == 4) && (id.ring() == 1 || id.ring() == 4));
|
| 1597 |
Bool_t ed2 = (idCSC.station() == 1) && ((idCSC.ring() == 2 && id.ring() == 2) || (idCSC.ring() == 3 && id.ring() == 3));
|
| 1598 |
Bool_t ed3 = (idCSC.station() != 1) && (idCSC.ring() == id.ring());
|
| 1599 |
Bool_t TMCSCMatch = (ed1 || ed2 || ed3);
|
| 1600 |
if(! TMCSCMatch)continue;
|
| 1601 |
|
| 1602 |
const CSCChamber* cscchamber = cscGeom->chamber(id);
|
| 1603 |
if (!cscchamber) continue;
|
| 1604 |
TrajectoryStateOnSurface TrajSuf_ = surfExtrapTrkSam(trackRef, cscchamber->toGlobal( (*segIt).localPosition() ).z());
|
| 1605 |
|
| 1606 |
Float_t dR_= fabs( TrajectoryDistToSeg( &TrajSuf_, segIt ) );
|
| 1607 |
if ( dR_ < deltaCSCR ){
|
| 1608 |
delete TrajSuf;
|
| 1609 |
TrajSuf=new TrajectoryStateOnSurface(TrajSuf_);
|
| 1610 |
deltaCSCR = dR_;
|
| 1611 |
cscSegOut = segIt;
|
| 1612 |
}
|
| 1613 |
}//loop over segments
|
| 1614 |
|
| 1615 |
return TrajSuf;
|
| 1616 |
}
|
| 1617 |
////////////// Get the matching with CSC-sgements...
|
| 1618 |
Bool_t TPTrackMuonSys::matchTTwithCSCSeg(Bool_t trackDir, Int_t st_, reco::TrackRef trackRef, edm::Handle<CSCSegmentCollection> cscSegments,
|
| 1619 |
CSCSegmentCollection::const_iterator &cscSegOut ) {
|
| 1620 |
Bool_t aMatch = false;
|
| 1621 |
// Float_t deltaCSCX = 9999.0, deltaCSCY = 9999.0;
|
| 1622 |
Float_t deltaCSCR = 50.0;
|
| 1623 |
for(CSCSegmentCollection::const_iterator segIt=cscSegments->begin(); segIt != cscSegments->end(); segIt++) {
|
| 1624 |
CSCDetId id = (CSCDetId)(*segIt).cscDetId();
|
| 1625 |
if((trackDir && id.endcap() == 2) || (!trackDir && id.endcap() == 1)) continue;
|
| 1626 |
if(
|
| 1627 |
(st_ == 0 && id.station() == 1 && ((id.ring() == 1) || (id.ring() == 4 ))) || // ME1/1,4
|
| 1628 |
(st_ == 1 && id.station() == 1 && (id.ring() == 2)) || // ME1/2
|
| 1629 |
(st_ == 2 && id.station() == 1 && (id.ring() == 3)) || // ME1/3
|
| 1630 |
(st_ == 3 && id.station() == 2) || // ME2/1,2
|
| 1631 |
(st_ == 4 && id.station() == 3) || // ME3/1,2
|
| 1632 |
(st_ == 5 && id.station() == 4) // ME4/1
|
| 1633 |
) {
|
| 1634 |
const CSCChamber* cscchamber = cscGeom->chamber(id);
|
| 1635 |
if (!cscchamber) continue;
|
| 1636 |
TrajectoryStateOnSurface TrajSuf = surfExtrapTrkSam(trackRef, cscchamber->toGlobal( (*segIt).localPosition() ).z());
|
| 1637 |
|
| 1638 |
Float_t dR_= fabs( TrajectoryDistToSeg( &TrajSuf, segIt ) );
|
| 1639 |
|
| 1640 |
if ( dR_ < deltaCSCR ){
|
| 1641 |
aMatch = true;
|
| 1642 |
deltaCSCR = dR_;
|
| 1643 |
cscSegOut = segIt;
|
| 1644 |
}
|
| 1645 |
}
|
| 1646 |
} // loop over segments
|
| 1647 |
return aMatch;
|
| 1648 |
}
|
| 1649 |
|
| 1650 |
////
|
| 1651 |
Int_t TPTrackMuonSys::getNLayerMatchedCSCSeg(CSCSegmentCollection::const_iterator &cscSegMatch,
|
| 1652 |
edm::Handle<CSCRecHit2DCollection> recHits,
|
| 1653 |
Float_t *delRecSegX,
|
| 1654 |
Float_t *delRecSegY,
|
| 1655 |
Int_t &nGhits
|
| 1656 |
)
|
| 1657 |
{
|
| 1658 |
Int_t nhits = 0;
|
| 1659 |
CSCDetId idCSC = (CSCDetId)(*cscSegMatch).cscDetId();
|
| 1660 |
Int_t ec = idCSC.endcap(), st = idCSC.station(), rg = idCSC.ring(), cham = idCSC.chamber();
|
| 1661 |
LocalPoint localPos = (*cscSegMatch).localPosition();
|
| 1662 |
|
| 1663 |
Float_t cscx = localPos.x();
|
| 1664 |
Float_t cscy = localPos.y();
|
| 1665 |
|
| 1666 |
for(Int_t i = 0; i < 6; i++){
|
| 1667 |
delRecSegX[i] = 999.;
|
| 1668 |
delRecSegY[i] = 999.;
|
| 1669 |
}
|
| 1670 |
Float_t minDeltaR = 999.;
|
| 1671 |
|
| 1672 |
for (CSCRecHit2DCollection::const_iterator recIt = recHits->begin(); recIt != recHits->end(); recIt++) {
|
| 1673 |
CSCDetId id = (CSCDetId)(*recIt).cscDetId();
|
| 1674 |
if(ec != id.endcap())continue;
|
| 1675 |
if(st != id.station())continue;
|
| 1676 |
if(rg != id.ring())continue;
|
| 1677 |
if(cham != id.chamber())continue;
|
| 1678 |
LocalPoint rhitlocal = (*recIt).localPosition();
|
| 1679 |
Float_t recx = rhitlocal.x();
|
| 1680 |
Float_t recy = rhitlocal.y();
|
| 1681 |
Float_t myDeltaR = sqrt(pow((cscx - recx),2)+ pow((cscy - recy),2));
|
| 1682 |
|
| 1683 |
Int_t layer = id.layer();
|
| 1684 |
|
| 1685 |
if(myDeltaR < minDeltaR){
|
| 1686 |
minDeltaR = myDeltaR;
|
| 1687 |
delRecSegX[layer-1] = cscx - recx;
|
| 1688 |
delRecSegY[layer-1] = cscy - recy;
|
| 1689 |
}
|
| 1690 |
nhits = nhits + 1;
|
| 1691 |
}
|
| 1692 |
|
| 1693 |
nGhits = 0;
|
| 1694 |
for(Int_t i = 0; i < 6; i++){
|
| 1695 |
Float_t myDeltaR = sqrt(pow(delRecSegX[i],2)+ pow(delRecSegY[i],2));
|
| 1696 |
if(myDeltaR < 15.)nGhits = nGhits + 1; /// look for this hit is within 15 cm.
|
| 1697 |
}
|
| 1698 |
|
| 1699 |
|
| 1700 |
return nhits; // this is probably not the right number of hits.
|
| 1701 |
}
|
| 1702 |
|
| 1703 |
|
| 1704 |
|
| 1705 |
|
| 1706 |
#ifdef GetCSCHitsBefore500
|
| 1707 |
void TPTrackMuonSys::getCSCSegWkeyHalfStrip(const std::vector<CSCRecHit2D> &theseRecHits, Float_t &cStrp, Float_t &ckWG){
|
| 1708 |
|
| 1709 |
#ifdef m_debug
|
| 1710 |
cout << " Start of the function " << endl;
|
| 1711 |
#endif
|
| 1712 |
|
| 1713 |
Float_t m_cStrp = 999.0, m_ckWG = 999.;
|
| 1714 |
|
| 1715 |
Float_t strpSum =0.0, wkeySum = 0.0;
|
| 1716 |
Int_t nhits =0 ;
|
| 1717 |
Bool_t layer3 = false;
|
| 1718 |
Bool_t me11a = false;
|
| 1719 |
|
| 1720 |
for ( std::vector<CSCRecHit2D>::const_iterator iRH = theseRecHits.begin(); iRH != theseRecHits.end(); iRH++) {
|
| 1721 |
CSCDetId idRH = (CSCDetId)(*iRH).cscDetId();
|
| 1722 |
Int_t mx_layer = idRH.layer();
|
| 1723 |
|
| 1724 |
if (idRH.ring() == 4) me11a = true;
|
| 1725 |
|
| 1726 |
/////////////////////////
|
| 1727 |
//
|
| 1728 |
// Get the strip number... ------> m_cStrp
|
| 1729 |
//
|
| 1730 |
/////////////////////////
|
| 1731 |
LocalPoint rhitlocal = (*iRH).localPosition();
|
| 1732 |
Float_t pWS = (*iRH).positionWithinStrip();
|
| 1733 |
// Find the strip containing this hit
|
| 1734 |
CSCRecHit2D::ChannelContainer hitstrips = (*iRH).channels();
|
| 1735 |
|
| 1736 |
Int_t nStrips = hitstrips.size();
|
| 1737 |
// Int_t centerid = nStrips/2 + 1;
|
| 1738 |
|
| 1739 |
if(nStrips == 3) m_cStrp = (int) (2.0*(hitstrips[1] + pWS - 0.5 ));
|
| 1740 |
else m_cStrp = (int) ( 2.0*(hitstrips[0] - pWS - 0.5 ));
|
| 1741 |
|
| 1742 |
Bool_t evenLayer = (mx_layer % 2 == 0);
|
| 1743 |
if ( evenLayer )m_cStrp -= 1;
|
| 1744 |
if ( (idRH.station() == 1) && (idRH.layer() != 3) && evenLayer ) m_cStrp += 1;
|
| 1745 |
//m_cStrp = hitstrips[centerid - 1];
|
| 1746 |
/////////////////////////
|
| 1747 |
//
|
| 1748 |
// Get the wireKey number...------> m_ckWG
|
| 1749 |
//
|
| 1750 |
/////////////////////////
|
| 1751 |
|
| 1752 |
CSCRecHit2D::ChannelContainer wiresg = (*iRH).wgroups();
|
| 1753 |
// m_ckWG = wiresg[0]-1;
|
| 1754 |
m_ckWG = wiresg[wiresg.size()/2]-1;//corrected from to wiresg[0] to [wiresg.size()/2]
|
| 1755 |
|
| 1756 |
if(mx_layer == 3){
|
| 1757 |
cStrp = m_cStrp;
|
| 1758 |
ckWG = m_ckWG;
|
| 1759 |
layer3 = true;
|
| 1760 |
}
|
| 1761 |
strpSum = strpSum + m_cStrp;
|
| 1762 |
wkeySum = wkeySum + m_ckWG;
|
| 1763 |
nhits = nhits + 1;
|
| 1764 |
}
|
| 1765 |
if(! layer3 && nhits != 0){
|
| 1766 |
cStrp = strpSum/nhits;
|
| 1767 |
ckWG = wkeySum/nhits;
|
| 1768 |
}
|
| 1769 |
if (me11a)cStrp = ((int)cStrp - 1)%32 + 1;
|
| 1770 |
}
|
| 1771 |
#else
|
| 1772 |
void TPTrackMuonSys::getCSCSegWkeyHalfStrip(const vector<CSCRecHit2D> &theseRecHits, Float_t &cStrp, Float_t &ckWG){
|
| 1773 |
Float_t strpSum = 0., wkeySum = 0.;
|
| 1774 |
Int_t nhits = 0; Bool_t me11a=false;
|
| 1775 |
for ( vector<CSCRecHit2D>::const_iterator iRH = theseRecHits.begin(); iRH != theseRecHits.end(); iRH++,nhits++ ) {
|
| 1776 |
CSCDetId idRH = (CSCDetId)(*iRH).cscDetId();
|
| 1777 |
if (idRH.ring() == 4) me11a = true;
|
| 1778 |
UChar_t layer = idRH.layer();
|
| 1779 |
Float_t pWS = (*iRH).positionWithinStrip();
|
| 1780 |
Int_t m_cStrp = (Int_t) (2.0*(iRH->channels(1) + pWS - 0.5 ));
|
| 1781 |
if ( layer%2!=0 &&
|
| 1782 |
!( (idRH.station() == 1) && (idRH.layer() != 3) ) ) m_cStrp--;
|
| 1783 |
|
| 1784 |
if( layer == 3 ) {
|
| 1785 |
cStrp=m_cStrp;
|
| 1786 |
ckWG=iRH->hitWire()-1;
|
| 1787 |
nhits=1;
|
| 1788 |
break;
|
| 1789 |
}
|
| 1790 |
else {
|
| 1791 |
strpSum += m_cStrp;
|
| 1792 |
wkeySum += iRH->hitWire()-1;
|
| 1793 |
}
|
| 1794 |
#ifdef m_debug
|
| 1795 |
clog << "ME" << idRH.station() << idRH.ring() <<":" << endl
|
| 1796 |
<< "cstrips ("<< iRH->nStrips() <<"): " //should be always 3
|
| 1797 |
<< iRH->channels(0) << "," << iRH->channels(1) << "," << iRH->channels(2)
|
| 1798 |
<< "; cwire:" << iRH->hitWire() << endl;
|
| 1799 |
#endif
|
| 1800 |
}
|
| 1801 |
cStrp = strpSum/nhits;
|
| 1802 |
ckWG = wkeySum/nhits;
|
| 1803 |
if (me11a) cStrp = ( (Int_t) cStrp - 1)%32 + 1;
|
| 1804 |
#ifdef m_debug
|
| 1805 |
clog << "cStrip: "<< cStrp << "; cWire: " << ckWG << endl;
|
| 1806 |
#endif
|
| 1807 |
}
|
| 1808 |
#endif// end of GetCSCHitsBefore500
|
| 1809 |
/*
|
| 1810 |
Bool_t TPTrackMuonSys::matchCSCSegWithLCT(edm::Handle<CSCCorrelatedLCTDigiCollection> mpclcts,
|
| 1811 |
CSCDetId & idCSC,
|
| 1812 |
Int_t TT,
|
| 1813 |
Float_t TrkPhi, Float_t TrkEta,
|
| 1814 |
Float_t c1, Float_t w1,
|
| 1815 |
CSCCorrelatedLCTDigiCollection::const_iterator &mpcItOut,
|
| 1816 |
CSCCorrelatedLCTDigiCollection::const_iterator &mpcHsWkOut,
|
| 1817 |
Bool_t *xMatch,
|
| 1818 |
Float_t *mDAngle,
|
| 1819 |
Float_t *diffTrkEta,
|
| 1820 |
Float_t *diffTrkPhi,
|
| 1821 |
Float_t *delHStrp,
|
| 1822 |
Float_t *delWkey
|
| 1823 |
)
|
| 1824 |
{
|
| 1825 |
/// Input
|
| 1826 |
// idCSC, halfstrip and wireKey
|
| 1827 |
|
| 1828 |
/// Output
|
| 1829 |
/// Whether there is a match based on halfstrip and wireKey
|
| 1830 |
/// or based on phi.
|
| 1831 |
/// also: minPhi, minDAngle;
|
| 1832 |
|
| 1833 |
for(Int_t i = 0; i < 2; i++){
|
| 1834 |
xMatch[i] = false;
|
| 1835 |
mDAngle[i] = 999.;
|
| 1836 |
diffTrkEta[i] = 999.;
|
| 1837 |
diffTrkPhi[i] = 999.;
|
| 1838 |
delHStrp[i] = 99999.;
|
| 1839 |
delWkey[i] = 99999.;
|
| 1840 |
}
|
| 1841 |
|
| 1842 |
Int_t ec = idCSC.endcap(), st = idCSC.station(), rg = idCSC.ring(), cham = idCSC.chamber();
|
| 1843 |
Bool_t me11a = false;
|
| 1844 |
if (rg == 4) me11a = true;
|
| 1845 |
|
| 1846 |
//// For the minimization:
|
| 1847 |
|
| 1848 |
for (CSCCorrelatedLCTDigiCollection::DigiRangeIterator detMPCUnitIt = mpclcts->begin(); detMPCUnitIt != mpclcts->end(); detMPCUnitIt++) {
|
| 1849 |
const CSCDetId& id = (*detMPCUnitIt).first;
|
| 1850 |
|
| 1851 |
if(ec != id.endcap())continue;
|
| 1852 |
if(st != id.station())continue;
|
| 1853 |
if(cham != id.chamber())continue;
|
| 1854 |
|
| 1855 |
Bool_t ed1 = (st == 1) && ((rg == 1 || rg == 4) && (id.ring() == 1 || id.ring() == 4));
|
| 1856 |
Bool_t ed2 = (st == 1) && ((rg == 2 && id.ring() == 2) || (rg == 3 && id.ring() == 3));
|
| 1857 |
Bool_t ed3 = (st != 1) && (rg == id.ring());
|
| 1858 |
Bool_t TMCSCMatch = (ed1 || ed2 || ed3);
|
| 1859 |
if(! TMCSCMatch)continue;
|
| 1860 |
|
| 1861 |
const CSCCorrelatedLCTDigiCollection::Range& MPCrange = (*detMPCUnitIt).second;
|
| 1862 |
for (CSCCorrelatedLCTDigiCollection::const_iterator mpcIt = MPCrange.first; mpcIt != MPCrange.second; mpcIt++) {
|
| 1863 |
Bool_t lct_valid = (*mpcIt).isValid();
|
| 1864 |
if(!lct_valid)continue;
|
| 1865 |
|
| 1866 |
Int_t station = id.station()-1;
|
| 1867 |
Int_t sector = id.triggerSector()-1;
|
| 1868 |
Int_t subSector = CSCTriggerNumbering::triggerSubSectorFromLabels(id);
|
| 1869 |
Int_t fpga = ( subSector ? subSector-1 : station+1 );
|
| 1870 |
Int_t endcap = id.endcap()-1;
|
| 1871 |
lclphidat lclPhi = srLUTs_[endcap][sector][fpga]->localPhi((*mpcIt).getStrip(),(*mpcIt).getPattern(),(*mpcIt).getQuality(),(*mpcIt).getBend() );
|
| 1872 |
gblphidat gblPhi = srLUTs_[endcap][sector][fpga]->globalPhiME( lclPhi.phi_local, (*mpcIt).getKeyWG(), id.triggerCscId() );
|
| 1873 |
gbletadat gblEta = srLUTs_[endcap][sector][fpga]->globalEtaME(lclPhi.phi_bend_local, lclPhi.phi_local,(*mpcIt).getKeyWG(), id.triggerCscId() );
|
| 1874 |
|
| 1875 |
UInt_t mpcphi = gblPhi.global_phi;
|
| 1876 |
UInt_t mpceta = gblEta.global_eta;
|
| 1877 |
double geta = theTriggerScales->getRegionalEtaScale(2)->getCenter(mpceta); //Type 2 is CSC
|
| 1878 |
double radphi = ((mpcphi)/4096.)*(62./180.)*M_PI; // cscphi in rad
|
| 1879 |
|
| 1880 |
double gphi = radphi + sector*M_PI/3.+ (14*M_PI/180.);//0.24434609;// Global Phi with respect to 0 degrees.
|
| 1881 |
if(sector == 5){
|
| 1882 |
if(radphi <=45*(M_PI/180.)) gphi = radphi + 315*M_PI/180.;
|
| 1883 |
else gphi = radphi - 45*M_PI/180.;
|
| 1884 |
}
|
| 1885 |
if(TrkPhi < 0)TrkPhi = TrkPhi + 2*M_PI;
|
| 1886 |
if(gphi < 0)gphi = gphi + 2*M_PI;
|
| 1887 |
Float_t DeltaAngle = sqrt(pow((TrkEta-geta),2) + pow((TrkPhi-gphi),2));
|
| 1888 |
|
| 1889 |
Int_t m_strip = (*mpcIt).getStrip(); // halfstrip that goes from 0 to 31
|
| 1890 |
if (me11a)m_strip = (m_strip-1)%32+1;
|
| 1891 |
Float_t m_keyWG = (*mpcIt).getKeyWG(); //
|
| 1892 |
Float_t m_delHStrp = c1 - m_strip;
|
| 1893 |
Float_t m_delWkey = w1 - m_keyWG;
|
| 1894 |
|
| 1895 |
if (me11a){ // the ganging of ME11a causes wraparound effects at the boundaries for delta strip
|
| 1896 |
if (m_delHStrp > 16) m_delHStrp -= 32;
|
| 1897 |
if (m_delHStrp < -16)m_delHStrp += 32;
|
| 1898 |
}
|
| 1899 |
|
| 1900 |
Float_t phiDiff = TrkPhi-gphi;
|
| 1901 |
if( fabs(phiDiff) < fabs(diffTrkPhi[0]) ){
|
| 1902 |
xMatch[0] = true;
|
| 1903 |
mDAngle[0] = DeltaAngle;
|
| 1904 |
diffTrkEta[0] = TrkEta-geta;
|
| 1905 |
diffTrkPhi[0] = phiDiff;
|
| 1906 |
delHStrp[0] = m_delHStrp;
|
| 1907 |
delWkey[0] = m_delWkey;
|
| 1908 |
mpcItOut = mpcIt;
|
| 1909 |
} // for the matching if statement...
|
| 1910 |
|
| 1911 |
if( (fabs(m_delHStrp) < fabs(delHStrp[1])) && (fabs(m_delWkey) < fabs(delWkey[1])) ){ // match strips at <=10 and wirekeys at <=5
|
| 1912 |
xMatch[1] = true;
|
| 1913 |
mDAngle[1] = DeltaAngle;
|
| 1914 |
diffTrkEta[1] = TrkEta-geta;
|
| 1915 |
diffTrkPhi[1] = TrkPhi-gphi;
|
| 1916 |
delHStrp[1] = m_delHStrp;
|
| 1917 |
delWkey[1] = m_delWkey;
|
| 1918 |
mpcHsWkOut = mpcIt;
|
| 1919 |
} // for the matching if statement...
|
| 1920 |
|
| 1921 |
#ifdef m_debug
|
| 1922 |
std::cout << "MPC E: " << id.endcap() << " R:" << id.ring() << " S: " << id.station() << " C: " << id.chamber()
|
| 1923 |
<< std::endl;
|
| 1924 |
#endif
|
| 1925 |
}
|
| 1926 |
}// end iteration over lcts_mpc_data
|
| 1927 |
return (xMatch[0] || xMatch[1]);
|
| 1928 |
}
|
| 1929 |
*/
|
| 1930 |
void TPTrackMuonSys::chamberCandidates(Int_t station, Float_t feta, Float_t phi, std::vector <int> &coupleOfChambers){
|
| 1931 |
// yeah, hardcoded again...
|
| 1932 |
coupleOfChambers.clear();
|
| 1933 |
|
| 1934 |
UChar_t ring = ringCandidate(station, feta, phi);
|
| 1935 |
if(ring){
|
| 1936 |
Float_t phi_zero = 0.;// check! the phi at the "edge" of Ch 1
|
| 1937 |
Float_t phi_const = 2.*M_PI/36.;
|
| 1938 |
Int_t first_chamber = 1, last_chamber = 36;
|
| 1939 |
if(1 != station && 1==ring){ // 18 chambers in the ring
|
| 1940 |
phi_const*=2;
|
| 1941 |
last_chamber /= 2;
|
| 1942 |
}
|
| 1943 |
if(phi < 0.) phi += 2*M_PI;
|
| 1944 |
|
| 1945 |
Float_t chamber_float = (phi - phi_zero)/phi_const;
|
| 1946 |
Int_t chamber_int = int(chamber_float);
|
| 1947 |
if (chamber_float - Float_t(chamber_int) -0.5 <0.){
|
| 1948 |
|
| 1949 |
if(0!=chamber_int ) coupleOfChambers.push_back(chamber_int);
|
| 1950 |
else coupleOfChambers.push_back(last_chamber);
|
| 1951 |
coupleOfChambers.push_back(chamber_int+1);
|
| 1952 |
|
| 1953 |
}else{
|
| 1954 |
coupleOfChambers.push_back(chamber_int+1);
|
| 1955 |
if(last_chamber!=chamber_int+1) coupleOfChambers.push_back(chamber_int+2);
|
| 1956 |
else coupleOfChambers.push_back(first_chamber);
|
| 1957 |
}
|
| 1958 |
}
|
| 1959 |
}
|
| 1960 |
|
| 1961 |
|
| 1962 |
Int_t TPTrackMuonSys::ringCandidate(Int_t station, Float_t feta, Float_t phi){
|
| 1963 |
switch (station){
|
| 1964 |
case 1:
|
| 1965 |
if(fabs(feta)>=0.85 && fabs(feta)<1.18){//ME13
|
| 1966 |
return 3;
|
| 1967 |
}
|
| 1968 |
if(fabs(feta)>=1.18 && fabs(feta)<=1.5){//ME12 if(fabs(feta)>1.18 && fabs(feta)<1.7){//ME12
|
| 1969 |
return 2;
|
| 1970 |
}
|
| 1971 |
if(fabs(feta)>1.5 && fabs(feta)<2.1){//ME11
|
| 1972 |
return 1;
|
| 1973 |
}
|
| 1974 |
if(fabs(feta)>=2.1 && fabs(feta)<2.45){//ME11
|
| 1975 |
return 4;
|
| 1976 |
}
|
| 1977 |
break;
|
| 1978 |
case 2:
|
| 1979 |
if(fabs(feta)>0.95 && fabs(feta)<1.6){//ME22
|
| 1980 |
return 2;
|
| 1981 |
}
|
| 1982 |
if(fabs(feta)>1.55 && fabs(feta)<2.45){//ME21
|
| 1983 |
return 1;
|
| 1984 |
}
|
| 1985 |
break;
|
| 1986 |
case 3:
|
| 1987 |
if(fabs(feta)>1.08 && fabs(feta)<1.72){//ME32
|
| 1988 |
return 2;
|
| 1989 |
}
|
| 1990 |
if(fabs(feta)>1.69 && fabs(feta)<2.45){//ME31
|
| 1991 |
return 1;
|
| 1992 |
}
|
| 1993 |
break;
|
| 1994 |
case 4:
|
| 1995 |
if(fabs(feta)>1.78 && fabs(feta) <2.45){//ME41
|
| 1996 |
return 1;
|
| 1997 |
}
|
| 1998 |
if(fabs(feta)>1.15 && fabs(feta) <=1.78){//ME42
|
| 1999 |
return 2;
|
| 2000 |
}
|
| 2001 |
break;
|
| 2002 |
default:
|
| 2003 |
LogError("")<<"Invalid station: "<<station<<endl;
|
| 2004 |
break;
|
| 2005 |
}
|
| 2006 |
return 0;
|
| 2007 |
}
|
| 2008 |
|
| 2009 |
UChar_t TPTrackMuonSys::thisChamberCandidate(UChar_t station, UChar_t ring, Float_t phi){
|
| 2010 |
//search for chamber candidate based on CMS IN-2007/024
|
| 2011 |
//10 deg chambers are ME1,ME22,ME32,ME42 chambers; 20 deg chambers are ME21,31,41 chambers
|
| 2012 |
//Chambers one always starts from approx -5 deg.
|
| 2013 |
const UChar_t nVal = (station>1 && ring == 1)?18:36;
|
| 2014 |
const Float_t ChamberSpan=2*M_PI/nVal;
|
| 2015 |
Float_t dphi = phi + M_PI/36;
|
| 2016 |
while (dphi >= 2*M_PI) dphi -= 2*M_PI;
|
| 2017 |
while (dphi < 0) dphi += 2*M_PI;
|
| 2018 |
UChar_t ChCand=floor(dphi/ChamberSpan)+1;
|
| 2019 |
return ChCand>nVal?nVal:ChCand;
|
| 2020 |
}
|
| 2021 |
|
| 2022 |
///////////////////////
|
| 2023 |
// to get the track position info at a particular rho
|
| 2024 |
TrajectoryStateOnSurface TPTrackMuonSys::cylExtrapTrkSam(reco::TrackRef track, double rho)
|
| 2025 |
{
|
| 2026 |
Cylinder::PositionType pos(0, 0, 0);
|
| 2027 |
Cylinder::RotationType rot;
|
| 2028 |
Cylinder::CylinderPointer myCylinder = Cylinder::build(pos, rot, rho);
|
| 2029 |
try{
|
| 2030 |
FreeTrajectoryState recoStart = freeTrajStateMuon(track);
|
| 2031 |
TrajectoryStateOnSurface recoProp = propagatorAlong->propagate(recoStart, *myCylinder);
|
| 2032 |
if (!recoProp.isValid()) {
|
| 2033 |
recoProp = propagatorOpposite->propagate(recoStart, *myCylinder);
|
| 2034 |
}
|
| 2035 |
return recoProp;
|
| 2036 |
}
|
| 2037 |
catch(cms::Exception){
|
| 2038 |
edm::LogError("")<<"invalid track extrapolation to cylinder"<<endl;
|
| 2039 |
TrajectoryStateOnSurface recoProp;
|
| 2040 |
return recoProp;
|
| 2041 |
}
|
| 2042 |
}
|
| 2043 |
|
| 2044 |
// to get track position at a particular (xy) plane given its z
|
| 2045 |
TrajectoryStateOnSurface TPTrackMuonSys::surfExtrapTrkSam(reco::TrackRef track, double z)
|
| 2046 |
{
|
| 2047 |
Plane::PositionType pos(0, 0, z);
|
| 2048 |
Plane::RotationType rot;
|
| 2049 |
Plane::PlanePointer myPlane = Plane::build(pos, rot);
|
| 2050 |
try{
|
| 2051 |
FreeTrajectoryState recoStart = freeTrajStateMuon(track);
|
| 2052 |
TrajectoryStateOnSurface recoProp = propagatorAlong->propagate(recoStart, *myPlane);
|
| 2053 |
if (!recoProp.isValid()) {
|
| 2054 |
recoProp = propagatorOpposite->propagate(recoStart, *myPlane);
|
| 2055 |
}
|
| 2056 |
return recoProp;
|
| 2057 |
}
|
| 2058 |
catch(cms::Exception){
|
| 2059 |
edm::LogError("")<<"invalid track extrapolation to plane"<<endl;
|
| 2060 |
TrajectoryStateOnSurface recoProp;
|
| 2061 |
return recoProp;
|
| 2062 |
}
|
| 2063 |
}
|
| 2064 |
|
| 2065 |
FreeTrajectoryState TPTrackMuonSys::freeTrajStateMuon(reco::TrackRef track)
|
| 2066 |
{
|
| 2067 |
GlobalPoint innerPoint(track->innerPosition().x(), track->innerPosition().y(), track->innerPosition().z());
|
| 2068 |
GlobalVector innerVec (track->innerMomentum().x(), track->innerMomentum().y(), track->innerMomentum().z());
|
| 2069 |
|
| 2070 |
GlobalTrajectoryParameters gtPars(innerPoint, innerVec, track->charge(), &*theBField);
|
| 2071 |
|
| 2072 |
//FreeTrajectoryState recoStart(innerPoint, innerVec, track->charge(), &*theBField);
|
| 2073 |
//return recoStart;
|
| 2074 |
|
| 2075 |
AlgebraicSymMatrix66 cov;
|
| 2076 |
cov *= 1e-20;
|
| 2077 |
|
| 2078 |
CartesianTrajectoryError tCov(cov);
|
| 2079 |
|
| 2080 |
return (cov.kRows == 6 ? FreeTrajectoryState(gtPars, tCov) : FreeTrajectoryState(gtPars)) ;
|
| 2081 |
}
|
| 2082 |
|
| 2083 |
|
| 2084 |
Bool_t TPTrackMuonSys::GetDecayChains(TrackingParticleRef tpr, HepMC::GenEvent *HepGenEvent, ULong64_t &truth_type, Int_t & truth_thesamewith, vector<vector< vector<Int_t> > > & ExistingSimChains) {
|
| 2085 |
Bool_t IsPileup=false;
|
| 2086 |
truth_thesamewith=-1;
|
| 2087 |
vector<TheTrackType> types;
|
| 2088 |
for (vector<SimTrack>::const_iterator g4Track_iter = tpr->g4Track_begin(); g4Track_iter != tpr->g4Track_end(); ++g4Track_iter ) {
|
| 2089 |
//g4Track loop begin
|
| 2090 |
DChain.clear();SimChains.clear();
|
| 2091 |
const SimTrack *thisTrk=&(*g4Track_iter);
|
| 2092 |
SimTrackDaughtersTree( thisTrk,tpr );
|
| 2093 |
Bool_t ChainEnd; TrackingParticleRef tpr_tmp=tpr;
|
| 2094 |
do {
|
| 2095 |
ChainEnd=true;
|
| 2096 |
if ( !thisTrk->noVertex() ) {
|
| 2097 |
TrackingVertexRef tvr=tpr_tmp->parentVertex();
|
| 2098 |
for ( TrackingParticleRefVector::iterator parenttp=tvr->sourceTracks_begin();parenttp!=tvr->sourceTracks_end();parenttp++ ) {
|
| 2099 |
for (vector<SimTrack>::const_iterator g4Trk_iter = (*parenttp)->g4Track_begin() ; g4Trk_iter != (*parenttp)->g4Track_end(); ++g4Trk_iter )
|
| 2100 |
if ( SVC[thisTrk->vertIndex()].parentIndex()==Int_t(g4Trk_iter->trackId()) && g4Trk_iter->eventId().rawId()==thisTrk->eventId().rawId()) {
|
| 2101 |
thisTrk=&(*g4Trk_iter);tpr_tmp=*parenttp;
|
| 2102 |
vector<const SimTrack *>::iterator SavedSimTrk_iter=find(SavedSimTrk.begin(),SavedSimTrk.end(),thisTrk);
|
| 2103 |
if (SavedSimTrk_iter==SavedSimTrk.end()) {
|
| 2104 |
DChain.push_back( MCParticlesList.size() );
|
| 2105 |
SavedSimTrk.push_back(thisTrk);
|
| 2106 |
MCParticlesList.push_back( MCParticleInfo_Creator(thisTrk,tpr_tmp) );
|
| 2107 |
}
|
| 2108 |
else DChain.push_back( FindSimTrackInMCParticlesList(SavedSimTrk_iter-SavedSimTrk.begin()));
|
| 2109 |
ChainEnd=false;
|
| 2110 |
break;
|
| 2111 |
}
|
| 2112 |
if (!ChainEnd) break;
|
| 2113 |
}
|
| 2114 |
}
|
| 2115 |
}while (!ChainEnd);
|
| 2116 |
for (vector< vector<Int_t> >::iterator SimChains_iter = SimChains.begin(); SimChains_iter != SimChains.end(); ++SimChains_iter )
|
| 2117 |
SimChains_iter->insert(SimChains_iter->begin(),DChain.rbegin(),DChain.rend());
|
| 2118 |
DChain.clear();
|
| 2119 |
|
| 2120 |
//HepMC Particles
|
| 2121 |
HepMCChains.clear();
|
| 2122 |
if (!thisTrk->noGenpart() ) {
|
| 2123 |
HepMC::GenParticle *genPar=HepGenEvent->barcode_to_particle(thisTrk->genpartIndex());
|
| 2124 |
if (genPar!=NULL) {
|
| 2125 |
vector<const SimTrack *>::iterator SavedSimTrk_iter=find(SavedSimTrk.begin(),SavedSimTrk.end(),thisTrk);
|
| 2126 |
if (SavedSimTrk_iter!=SavedSimTrk.end()) MCParticlesList[FindSimTrackInMCParticlesList(SavedSimTrk_iter-SavedSimTrk.begin())].IsParticleFromGenerator=true;
|
| 2127 |
else LogError("CodeWrong")<<"Cannot find the simulated track in saved sim tracks.";
|
| 2128 |
HepMCParentTree( genPar );
|
| 2129 |
}
|
| 2130 |
else LogWarning("RefNull")<<"Either SimTrack::genpartIndex() or HepMC::GenParticle::barcode is wrong. Pilup track?";
|
| 2131 |
}
|
| 2132 |
|
| 2133 |
//merge the HepMC and SimTrack Decay Chains
|
| 2134 |
for (vector< vector<Int_t> >::iterator SimChains_iter = SimChains.begin(); SimChains_iter != SimChains.end(); ++SimChains_iter )
|
| 2135 |
if ( !HepMCChains.empty() )
|
| 2136 |
for (vector< vector<Int_t> >::iterator HepMCChains_iter = HepMCChains.begin(); HepMCChains_iter != HepMCChains.end(); ++HepMCChains_iter ) {
|
| 2137 |
vector<Int_t> thisChain(HepMCChains_iter->rbegin(),HepMCChains_iter->rend());
|
| 2138 |
thisChain.insert(thisChain.end(),SimChains_iter->begin(),SimChains_iter->end());
|
| 2139 |
IsPileup=SaveAndClassify(thisChain,types, truth_thesamewith, ExistingSimChains);
|
| 2140 |
}
|
| 2141 |
else IsPileup=SaveAndClassify( *SimChains_iter,types, truth_thesamewith, ExistingSimChains );
|
| 2142 |
}//g4Track loop end
|
| 2143 |
for (vector <TheTrackType>::iterator type_iter=types.begin();type_iter!=types.end();type_iter++)
|
| 2144 |
truth_type=truth_type*100+Long64_t(*type_iter);
|
| 2145 |
return IsPileup;
|
| 2146 |
}
|
| 2147 |
|
| 2148 |
void TPTrackMuonSys::SimTrackDaughtersTree(const SimTrack * thisTrk, TrackingParticleRef tpr)
|
| 2149 |
{
|
| 2150 |
// Find MC Truth Segment - the offical one use chi2 to match simtrk (MuonIdTruthInfo.cc) and it won't know the decay in flight segment truth
|
| 2151 |
// The particle type of the hit may differ from the particle type of the SimTrack with id trackId().
|
| 2152 |
// This happends if the hit was created by a secondary track(e.g. a delta ray) originating from the trackId() and not existing as aseparate SimTrack.
|
| 2153 |
// ( particle type match notice is from haiyun.teng@cern.ch )
|
| 2154 |
Bool_t ChainEnd=true;
|
| 2155 |
//To avoid duplicate particle saving
|
| 2156 |
vector<const SimTrack *>::iterator SavedSimTrk_iter=find(SavedSimTrk.begin(),SavedSimTrk.end(),thisTrk);
|
| 2157 |
if (SavedSimTrk_iter==SavedSimTrk.end()) {
|
| 2158 |
DChain.push_back( MCParticlesList.size() );
|
| 2159 |
SavedSimTrk.push_back(thisTrk);
|
| 2160 |
MCParticlesList.push_back( MCParticleInfo_Creator(thisTrk,tpr) );
|
| 2161 |
}
|
| 2162 |
else DChain.push_back( FindSimTrackInMCParticlesList(SavedSimTrk_iter-SavedSimTrk.begin()) );
|
| 2163 |
|
| 2164 |
for ( TrackingVertexRefVector::iterator tvr=tpr->decayVertices().begin();tvr!=tpr->decayVertices().end();tvr++ )
|
| 2165 |
for ( TrackingParticleRefVector::iterator daughtertp=(*tvr)->daughterTracks_begin();daughtertp!=(*tvr)->daughterTracks_end();daughtertp++ )
|
| 2166 |
for (vector<SimTrack>::const_iterator g4Trk_iter = (*daughtertp)->g4Track_begin() ; g4Trk_iter != (*daughtertp)->g4Track_end(); ++g4Trk_iter )
|
| 2167 |
if ( SVC[g4Trk_iter->vertIndex()].parentIndex()==Int_t(thisTrk->trackId()) && g4Trk_iter->eventId().rawId()==thisTrk->eventId().rawId()) {
|
| 2168 |
ChainEnd=false;
|
| 2169 |
SimTrackDaughtersTree( &(*g4Trk_iter), *daughtertp );
|
| 2170 |
}
|
| 2171 |
if (ChainEnd) SimChains.push_back(DChain);
|
| 2172 |
DChain.pop_back();
|
| 2173 |
}
|
| 2174 |
|
| 2175 |
void TPTrackMuonSys::HepMCParentTree(HepMC::GenParticle *genPar) {
|
| 2176 |
HepMC::GenVertex *thisVtx = genPar->production_vertex();
|
| 2177 |
Bool_t ChainEnd=true;
|
| 2178 |
if (thisVtx) {
|
| 2179 |
for (HepMC::GenVertex::particles_in_const_iterator pgenD = thisVtx->particles_in_const_begin(); pgenD != thisVtx->particles_in_const_end(); ++pgenD)
|
| 2180 |
if ((*pgenD)->pdg_id()!=92) {//Pythia special code for string, we only care about the particles after hadronization
|
| 2181 |
ChainEnd=false;
|
| 2182 |
vector<HepMC::GenParticle *>::iterator SavedHepPar_iter=find(SavedHepPar.begin(),SavedHepPar.end(),*pgenD);
|
| 2183 |
if (SavedHepPar_iter==SavedHepPar.end())
|
| 2184 |
{
|
| 2185 |
DChain.push_back(MCParticlesList.size());
|
| 2186 |
SavedHepPar.push_back(*pgenD);
|
| 2187 |
MCParticlesList.push_back( MCParticleInfo_Creator( (*pgenD) ) );
|
| 2188 |
}
|
| 2189 |
else DChain.push_back(FindHepMCInMCParticlesList(SavedHepPar_iter-SavedHepPar.begin()));
|
| 2190 |
HepMCParentTree(*pgenD);
|
| 2191 |
DChain.pop_back();
|
| 2192 |
}
|
| 2193 |
}
|
| 2194 |
if (ChainEnd) HepMCChains.push_back(DChain);
|
| 2195 |
}
|
| 2196 |
|
| 2197 |
Bool_t TPTrackMuonSys::SaveAndClassify(vector<Int_t> &Chain, vector<TheTrackType> &types, Int_t &truth_thesamewith, vector< vector< vector<Int_t> > > &ExistingSimChains) {
|
| 2198 |
//Find out if another track is just a decay product or mother particle of this track
|
| 2199 |
truth_thesamewith=-1;
|
| 2200 |
Bool_t truth_isPileup=false;
|
| 2201 |
for ( vector< vector< vector<Int_t> > >::const_iterator atrack=ExistingSimChains.begin(); atrack!=ExistingSimChains.end(); atrack++ ) {
|
| 2202 |
for ( vector< vector<Int_t> >::const_iterator adecaychain=atrack->begin(); adecaychain!=atrack->end(); adecaychain++ )
|
| 2203 |
if ( IstheSameDChain(Chain,*adecaychain) ) {
|
| 2204 |
truth_thesamewith=atrack-ExistingSimChains.begin();
|
| 2205 |
break;
|
| 2206 |
}
|
| 2207 |
if ( truth_thesamewith>-1 ) break;
|
| 2208 |
}
|
| 2209 |
|
| 2210 |
//Save
|
| 2211 |
ExistingSimChains.back().push_back(Chain);
|
| 2212 |
|
| 2213 |
//Classification
|
| 2214 |
TheTrackType type=Classify(Chain);
|
| 2215 |
if ( !MCParticlesList[Chain[0]].IsParticleFromGenerator ) truth_isPileup=true;
|
| 2216 |
#ifdef jz_debug
|
| 2217 |
// if ( Int_t(type)<2||Int_t(type)>100|| truth_isPileup ) {
|
| 2218 |
cerr<<"======================================================="<<endl;
|
| 2219 |
for (vector<Int_t>::reverse_iterator iter = Chain.rbegin(); iter != Chain.rend(); iter++)
|
| 2220 |
cerr<<"("<<MCParticlesList[*iter].pdgId<<","<<MCParticlesList[*iter].DoesParticleHaveMuonHit<<") <--|";
|
| 2221 |
cerr<<endl;
|
| 2222 |
cerr<<"MuTagpdgID:"<<MuTagtracktruth_id<<"; pdgID:"<<tracktruth_id;
|
| 2223 |
cerr<<"; type is "<<Int_t(type)<<"; pileup: "<<Int_t(truth_isPileup)<<endl
|
| 2224 |
<<"======================================================="<<endl;
|
| 2225 |
// }
|
| 2226 |
#endif
|
| 2227 |
if ( types.empty() ) types.push_back(type);
|
| 2228 |
else if (type!=Others) {
|
| 2229 |
if ( find(types.begin(),types.end(),type)==types.end() ) {
|
| 2230 |
vector<TheTrackType>::iterator DeterminedAsOthers=find(types.begin(),types.end(),Others);
|
| 2231 |
if ( DeterminedAsOthers!=types.end() ) *DeterminedAsOthers=type;
|
| 2232 |
else {
|
| 2233 |
vector<TheTrackType>::iterator DeterminedAsNoMuSysHit=find(types.begin(),types.end(),Others);
|
| 2234 |
if ( DeterminedAsNoMuSysHit!=types.end() ) *DeterminedAsOthers=type;
|
| 2235 |
else types.push_back(type);
|
| 2236 |
}
|
| 2237 |
}
|
| 2238 |
}
|
| 2239 |
return truth_isPileup;
|
| 2240 |
}
|
| 2241 |
|
| 2242 |
inline TPTrackMuonSys::ParticleType TPTrackMuonSys::ParticleCata(Int_t pid)
|
| 2243 |
{
|
| 2244 |
pid=pid>0?pid:pid*-1;
|
| 2245 |
if (pid==13) return Muon;
|
| 2246 |
if (pid==443) return JPsi;
|
| 2247 |
if (pid==23) return Z;
|
| 2248 |
if (pid==24) return W;
|
| 2249 |
if ((pid>110&&pid<900)||(pid>10000))
|
| 2250 |
{
|
| 2251 |
pid=pid%1000;
|
| 2252 |
if (pid<400) return LightMeson;
|
| 2253 |
if (pid>400&&pid<440) return CharmedMeson;
|
| 2254 |
if (pid>440&&pid<500) return ccbarMeson;
|
| 2255 |
if (pid>500&&pid<550) return BottomMeson;
|
| 2256 |
if (pid>550&&pid<600) return bbarMeson;
|
| 2257 |
}
|
| 2258 |
if (pid>1111&&pid<3350&&(pid%100)/10>0) return LightBaryon;
|
| 2259 |
if (pid>4111&&pid<4450&&(pid%100)/10>0) return CharmedBaryon;
|
| 2260 |
if (pid>5111&&pid<5560&&(pid%100)/10>0) return BottomBaryon;
|
| 2261 |
if (pid>1000&&pid<6000&&(pid%100)/10==0) return DiQuarks;
|
| 2262 |
if (pid>10&&pid<19) return Lepton;
|
| 2263 |
return Other;
|
| 2264 |
}
|
| 2265 |
|
| 2266 |
TPTrackMuonSys::TheTrackType TPTrackMuonSys::Classify(vector<Int_t> &Chain) {
|
| 2267 |
Int_t MuPos=-1;
|
| 2268 |
for (vector<Int_t>::reverse_iterator iter = Chain.rbegin(); iter != Chain.rend(); iter++) {
|
| 2269 |
ParticleType ParType=ParticleCata(MCParticlesList[*iter].pdgId);
|
| 2270 |
if (ParType==Muon&&MuPos<0) {//only care of the last muon
|
| 2271 |
MuPos=*iter;
|
| 2272 |
if ( !MCParticlesList[MuPos].DoesParticleHaveMuonHit ) return NoMuSysHit;
|
| 2273 |
continue;
|
| 2274 |
}
|
| 2275 |
if (MuPos>=0&&(ParType==W||ParType==Z)) {
|
| 2276 |
if ( MCParticlesList[MuPos].IsParticleFromGenerator ) {
|
| 2277 |
if (ParType==W) return PromptMuFromW;
|
| 2278 |
else return PromptMuFromZ;
|
| 2279 |
}
|
| 2280 |
else return NotPromptMufromWZ;//should never happen
|
| 2281 |
}
|
| 2282 |
if (MuPos>=0&&(ParType==JPsi)) return PromptMuFromJPsi;
|
| 2283 |
if (MuPos<0&&ParType>=LightMeson&&ParType<=DiQuarks) {
|
| 2284 |
if ( MCParticlesList[*iter].DoesParticleHaveMuonHit ) return PunchThrough;
|
| 2285 |
else return NoMuSysHit;
|
| 2286 |
}
|
| 2287 |
if (MuPos>=0&&ParType>=LightMeson&&ParType<=BottomBaryon) {
|
| 2288 |
if ( MCParticlesList[MuPos].IsParticleFromGenerator||MCParticlesList[MuPos].IsParticleBornInsideOfARegion ) {//inside pixel detector
|
| 2289 |
if (ParType==LightMeson) return PromptMuFromLightMeson;
|
| 2290 |
if (ParType>LightMeson&&ParType<LightBaryon) return PromptMuFromHeavyMeson;
|
| 2291 |
if (ParType==LightBaryon) return PromptMuFromLightBaryon;
|
| 2292 |
if (ParType>LightBaryon&&ParType<=BottomBaryon) return PromptMuFromHeavyBaryon;
|
| 2293 |
}
|
| 2294 |
if ( MCParticlesList[*iter].DoesParticleHaveMuonHit ) return PunchThroughAndDecayInFlight;
|
| 2295 |
else {
|
| 2296 |
if (ParType==LightMeson) return DecayInFlightFromLightMeson;
|
| 2297 |
if (ParType>LightMeson&&ParType<LightBaryon) return DecayInFlightFromHeavyMeson;
|
| 2298 |
if (ParType==LightBaryon) return DecayInFlightFromLightBaryon;
|
| 2299 |
if (ParType>LightBaryon&&ParType<=BottomBaryon) return DecayInFlightFromHeavyBaryon;
|
| 2300 |
}
|
| 2301 |
}
|
| 2302 |
}
|
| 2303 |
if (MuPos>=0) return PromptMuFromOthers;
|
| 2304 |
else return Others;
|
| 2305 |
}
|
| 2306 |
|
| 2307 |
Bool_t TPTrackMuonSys::IstheSameDChain(const vector<Int_t> &Chain1,const vector<Int_t> &Chain2) {//whether two chains include each other
|
| 2308 |
Bool_t ChainIncluded=false;
|
| 2309 |
vector<Int_t>::const_iterator Chain1_Particle = Chain1.begin(),Chain2_Particle = Chain2.begin();
|
| 2310 |
for (; Chain2_Particle != Chain2.end()&&Chain1_Particle != Chain1.end(); Chain2_Particle++) {
|
| 2311 |
if (ChainIncluded&&*Chain1_Particle!=*Chain2_Particle) ChainIncluded=false;
|
| 2312 |
if (Chain1.front() == *Chain2_Particle) {
|
| 2313 |
ChainIncluded=true;
|
| 2314 |
Chain1_Particle = Chain1.begin();
|
| 2315 |
}
|
| 2316 |
if (ChainIncluded) Chain1_Particle++;
|
| 2317 |
}
|
| 2318 |
if (!ChainIncluded) {
|
| 2319 |
Chain2_Particle = Chain2.begin(); Chain1_Particle = Chain1.begin();
|
| 2320 |
for (; Chain2_Particle != Chain2.end()&&Chain1_Particle != Chain1.end(); Chain1_Particle++) {
|
| 2321 |
if (ChainIncluded&&*Chain1_Particle!=*Chain2_Particle) ChainIncluded=false;
|
| 2322 |
if (Chain2.front() == *Chain1_Particle) {
|
| 2323 |
Chain2_Particle = Chain2.begin();
|
| 2324 |
ChainIncluded=true;
|
| 2325 |
}
|
| 2326 |
if (ChainIncluded) Chain2_Particle++;
|
| 2327 |
}
|
| 2328 |
}
|
| 2329 |
return ChainIncluded;
|
| 2330 |
}
|
| 2331 |
|
| 2332 |
TPTrackMuonSys::MCParticleInfo TPTrackMuonSys::MCParticleInfo_Creator( const SimTrack * thisTrk, TrackingParticleRef tpr ) {
|
| 2333 |
MCParticleInfo TBA;
|
| 2334 |
TBA.IsThisFromSimTrk=true;
|
| 2335 |
TBA.IsParticleFromGenerator=false;TBA.DoesParticleHaveMuonHit=false;
|
| 2336 |
TBA.IsParticleBornInsideOfARegion=false;
|
| 2337 |
TBA.IsPileup=false;
|
| 2338 |
TBA.pdgId=thisTrk->type();
|
| 2339 |
for ( vector<PSimHit>::const_iterator g4Hit_iter=tpr->pSimHit_begin();g4Hit_iter!=tpr->pSimHit_end();g4Hit_iter++ )
|
| 2340 |
if ( g4Hit_iter->trackId()==thisTrk->trackId() && g4Hit_iter->eventId().rawId()==thisTrk->eventId().rawId() && g4Hit_iter->particleType() == thisTrk->type() ) {
|
| 2341 |
DetId DetectorId( g4Hit_iter->detUnitId() );
|
| 2342 |
if ( DetectorId.det() == DetId::Muon ) {
|
| 2343 |
TBA.DoesParticleHaveMuonHit=true;
|
| 2344 |
break;
|
| 2345 |
}
|
| 2346 |
}
|
| 2347 |
if (!thisTrk->noVertex()) {
|
| 2348 |
SimVertex thisVtx=SVC[thisTrk->vertIndex()];
|
| 2349 |
//if ( (thisVtx.position().Perp2()<161604.&&deltaZ<568.)||(deltaRPhi2>82024.96&&deltaRPhi2<161604.&&deltaZ<666.)) //inside HCAL
|
| 2350 |
if ( thisVtx.position().Perp2()<100&&
|
| 2351 |
fabs(thisVtx.position().z())<30 ) //inside pixel detector
|
| 2352 |
TBA.IsParticleBornInsideOfARegion=true;
|
| 2353 |
}
|
| 2354 |
else TBA.IsPileup=true;
|
| 2355 |
return TBA;
|
| 2356 |
}
|
| 2357 |
|
| 2358 |
//deadzone center is according to http://cmssdt.cern.ch/SDT/lxr/source/RecoLocalMuon/CSCEfficiency/src/CSCEfficiency.cc#605
|
| 2359 |
//wire spacing is according to CSCTDR
|
| 2360 |
Float_t TPTrackMuonSys::YDistToHVDeadZone(Float_t yLocal, Int_t StationAndRing){
|
| 2361 |
//the ME11 wires are not parallel to x, but no gap
|
| 2362 |
//chamber edges are not included.
|
| 2363 |
const Float_t deadZoneCenterME1_2[2] = {-32.88305,32.867423};
|
| 2364 |
const Float_t deadZoneCenterME1_3[2] = {-22.7401,27.86665};
|
| 2365 |
const Float_t deadZoneCenterME2_1[2] = {-27.47,33.67};
|
| 2366 |
const Float_t deadZoneCenterME3_1[2] = {-36.21,23.68};
|
| 2367 |
const Float_t deadZoneCenterME4_1[2] = {-26.14,23.85};
|
| 2368 |
const Float_t deadZoneCenterME234_2[4] = {-81.8744,-21.18165,39.51105,100.2939};
|
| 2369 |
const Float_t *deadZoneCenter;
|
| 2370 |
Float_t deadZoneHeightHalf=0.32*7/2;// wire spacing * (wires missing + 1)/2
|
| 2371 |
Float_t minY=999999.;
|
| 2372 |
UChar_t nGaps=2;
|
| 2373 |
switch (abs(StationAndRing)) {
|
| 2374 |
case 11:
|
| 2375 |
case 14:
|
| 2376 |
return 162;//the height of ME11
|
| 2377 |
break;
|
| 2378 |
case 12:
|
| 2379 |
deadZoneCenter=deadZoneCenterME1_2;
|
| 2380 |
break;
|
| 2381 |
case 13:
|
| 2382 |
deadZoneCenter=deadZoneCenterME1_3;
|
| 2383 |
break;
|
| 2384 |
case 21:
|
| 2385 |
deadZoneCenter=deadZoneCenterME2_1;
|
| 2386 |
break;
|
| 2387 |
case 31:
|
| 2388 |
deadZoneCenter=deadZoneCenterME3_1;
|
| 2389 |
break;
|
| 2390 |
case 41:
|
| 2391 |
deadZoneCenter=deadZoneCenterME4_1;
|
| 2392 |
break;
|
| 2393 |
default:
|
| 2394 |
deadZoneCenter=deadZoneCenterME234_2;
|
| 2395 |
nGaps=4;
|
| 2396 |
}
|
| 2397 |
for ( UChar_t iGap=0;iGap<nGaps;iGap++ ) {
|
| 2398 |
Float_t newMinY=yLocal<deadZoneCenter[iGap]?deadZoneCenter[iGap]-deadZoneHeightHalf-yLocal:yLocal-(deadZoneCenter[iGap]+deadZoneHeightHalf);
|
| 2399 |
if ( newMinY<minY ) minY=newMinY;
|
| 2400 |
}
|
| 2401 |
return minY;
|
| 2402 |
}
|
| 2403 |
|
| 2404 |
//define this as a plug-in
|
| 2405 |
DEFINE_FWK_MODULE(TPTrackMuonSys);
|
| 2406 |
|
| 2407 |
|
