ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/CSCPriEff/src/TPTrackMuonSys.cc
Revision: 1.4
Committed: Wed Mar 20 11:53:08 2013 UTC (12 years, 1 month ago) by zhangjin
Content type: text/plain
Branch: MAIN
Changes since 1.3: +59 -146 lines
Log Message:
ME11AandB distinguishable,searching of chambercandidate improved

File Contents

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