File indexing completed on 2023-03-17 11:23:32
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateClosestToPoint.h"
0020
0021
0022 #include <memory>
0023
0024
0025 #include "FWCore/Framework/interface/Frameworkfwd.h"
0026 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/Framework/interface/EventSetup.h"
0029 #include "FWCore/Framework/interface/MakerMacros.h"
0030
0031 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0032 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0033
0034
0035
0036
0037 #include "DataFormats/VertexReco/interface/Vertex.h"
0038
0039 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0040 #include "MagneticField/VolumeBasedEngine/interface/VolumeBasedMagneticField.h"
0041
0042 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
0043 #include "DataFormats/Candidate/interface/VertexCompositeCandidate.h"
0044
0045 #include "SimDataFormats/Track/interface/SimTrack.h"
0046 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0047 #include "SimDataFormats/Vertex/interface/SimVertex.h"
0048 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0049
0050 #include "Geometry/CommonDetUnit/interface/TrackingGeometry.h"
0051 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0052 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0053 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0054 #include "Geometry/CommonDetUnit/interface/GluedGeomDet.h"
0055
0056 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0057
0058 #include "FWCore/ServiceRegistry/interface/Service.h"
0059 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0060
0061 #include "TH1.h"
0062 #include "TH2.h"
0063 #include "TFile.h"
0064 #include <fstream>
0065
0066
0067
0068
0069 class V0Analyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0070 public:
0071 explicit V0Analyzer(const edm::ParameterSet&);
0072 ~V0Analyzer();
0073
0074 private:
0075
0076 virtual void beginJob();
0077 virtual void analyze(const edm::Event&, const edm::EventSetup&);
0078 virtual void endJob();
0079
0080 const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> m_mfToken;
0081 const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> m_geomToken;
0082
0083 std::string algoLabel;
0084 std::string recoAlgoLabel;
0085 std::string SimTkLabel;
0086 std::string SimVtxLabel;
0087 std::string outputFile;
0088 std::string HistoFileName;
0089 std::string V0CollectionName;
0090 std::string cmsswVersion;
0091
0092
0093 TH1D* myKshortMassHisto;
0094 TH1D* myDecayRadiusHisto;
0095 TH1D* myRefitKshortMassHisto;
0096 TH1D* myRefitDecayRadiusHisto;
0097 TH1D* myNativeParticleMassHisto;
0098 TH1D* myEtaEfficiencyHisto;
0099 TH1D* mySimEtaHisto;
0100 TH1D* myRhoEfficiencyHisto;
0101 TH1D* myRhoEfficiencyHisto2;
0102 TH1D* myRhoEfficiencyHisto3;
0103 TH1D* myRhoEfficiencyHisto4;
0104 TH1D* mySimRhoHisto;
0105 TH1D* myKshortPtHisto;
0106 TH1D* myImpactParameterHisto;
0107 TH1D* myImpactParameterHisto2;
0108 TH1D* myNumSimKshortsHisto;
0109 TH1D* myNumRecoKshortsHisto;
0110 TH1D* myInnermostHitDistanceHisto;
0111
0112
0113 TH1D* step1massHisto;
0114 TH1D* step2massHisto;
0115 TH1D* step3massHisto;
0116 TH1D* step4massHisto;
0117
0118 TH1D* vertexChi2Histo;
0119 TH1D* rVtxHisto1;
0120 TH1D* vtxSigHisto1;
0121 TH1D* rVtxHisto2;
0122
0123 TH1D* vtxSigHisto2;
0124
0125 TH1D* rErrorHisto;
0126 TH1D* tkPtHisto;
0127 TH1D* k0sPtHisto;
0128 TH1D* tkChi2Histo;
0129 TH1D* sqrtTkChi2Histo;
0130 TH1D* tkEtaHisto;
0131 TH1D* kShortEtaHisto;
0132 TH1D* numHitsHisto;
0133
0134 TH1D* simTkMpipiHisto;
0135 int numDiff1, numDiff2;
0136
0137 std::ofstream hitsOut;
0138
0139
0140 };
0141
0142
0143
0144
0145
0146 const double piMassSq = 0.019479101;
0147
0148
0149
0150
0151
0152
0153
0154
0155 V0Analyzer::V0Analyzer(const edm::ParameterSet& iConfig)
0156 : m_mfToken(esConsumes()),
0157 m_geomToken(esConsumes()),
0158 algoLabel(iConfig.getUntrackedParameter("recoAlgorithm", std::string("ctfV0Prod"))),
0159 recoAlgoLabel(iConfig.getUntrackedParameter("trackingAlgo", std::string("ctfWithMaterialTracks"))),
0160 SimTkLabel(iConfig.getUntrackedParameter("moduleLabelTk", std::string("g4SimHits"))),
0161 SimVtxLabel(iConfig.getUntrackedParameter("moduleLabelVtx", std::string("g4SimHits"))),
0162 HistoFileName(iConfig.getUntrackedParameter("histoFileName", std::string("vtxAnalyzerHistos.root"))),
0163 V0CollectionName(iConfig.getUntrackedParameter("v0CollectionName", std::string("Kshort"))),
0164 cmsswVersion(iConfig.getUntrackedParameter("cmsswVersion", std::string("200"))) {
0165
0166
0167
0168 usesResource(TFileService::kSharedResource);
0169
0170 numDiff1 = numDiff2 = 0;
0171 }
0172
0173 V0Analyzer::~V0Analyzer() {
0174
0175
0176 }
0177
0178
0179
0180
0181
0182
0183
0184 void V0Analyzer::beginJob() {
0185 using std::string;
0186
0187 edm::Service<TFileService> fs;
0188
0189
0190
0191 string algo(algoLabel);
0192 string vernum = cmsswVersion;
0193 algo += vernum;
0194 string refit("Refitted");
0195 string native("Native");
0196
0197 string massHistoLabelLong(algo + string(" Invariant Mass of K^{0}_{s}"));
0198 string massHistoLabelShort(vernum + string("K0sInvarMass"));
0199
0200 string decRadHistoLabelLong(algo + string(" Distance of Decay Vtx from z-axis"));
0201 string decRadHistoLabelShort(vernum + string("K0sDecayRad"));
0202
0203 string etaEffHistoLabelLong(algo + string(" Reconstruction Efficiency vs. #eta"));
0204 string etaEffHistoLabelShort(vernum + string("EffVsEta"));
0205
0206 string rhoEffHistoLabelLong(algo + string(" Reconstruction Efficiency vs. #rho"));
0207 string rhoEffHistoLabelShort(vernum + string("EffVsRho"));
0208
0209 string rhoEffHistoLabel2Long(algo + string(" 2-Reconstruction Efficiency vs. #rho"));
0210 string rhoEffHistoLabel2Short(vernum + string("EffVsRho2"));
0211
0212 string rhoEffHistoLabel3Long(algo + string(" 3-Reconstruction Efficiency vs. #rho"));
0213 string rhoEffHistoLabel3Short(vernum + string("EffVsRho3"));
0214
0215 string rhoEffHistoLabel4Long(algo + string(" 4-Reconstruction Efficiency vs. #rho"));
0216 string rhoEffHistoLabel4Short(vernum + string("EffVsRho4"));
0217
0218 string simRhoHistoLabelLong(algo + string(" Simulated K^{0}_{s} #rho"));
0219 string simRhoHistoLabelShort(vernum + string("SimRho"));
0220
0221 string simRHistoLabelLong(algo + string(" Simulated K^{0}_{s} R"));
0222 string simRHistoLabelShort(vernum + string("SimR"));
0223
0224 string simEtaHistoLabelLong(algo + string(" Simulated K^{0}_{s} #eta"));
0225 string simEtaHistoLabelShort(vernum + string("SimEta"));
0226
0227 string kShortPtHistoLabelLong(algo + string(" K^{0}_{s} P_{t}"));
0228 string kShortPtHistoLabelShort(vernum + string("K0sPt"));
0229
0230 string impactParamHistoLabelLong(algo + string(" Impact parameter of reco tracks"));
0231 string impactParamHistoLabelShort(vernum + string("d0"));
0232
0233 string impactParamHisto2LabelLong(algo + string(" distance from beam line to reco tracks"));
0234 string impactParamHisto2LabelShort(vernum + string("d0Radial"));
0235
0236 string numSimKshortsHistoLabelLong(algo + string(" number of simulated Kshorts"));
0237 string numSimKshortsHistoLabelShort(vernum + string("numSimK0s"));
0238
0239 string numRecoKshortsHistoLabelLong(algo + string(" number of reconstructed Kshorts"));
0240 string numRecoKshortsHistoLabelShort(vernum + string("numRecoK0s"));
0241
0242 string innermostHitDistanceHistoLabelLong(algo + string(" distance between innermost hits of track pairs"));
0243 string innermostHitDistanceHistoLabelShort(vernum + string("hitDist"));
0244
0245
0246
0247 string s1massHistoLabelLong(algo + string(" m_{#pi #pi}, all tracks"));
0248 string s1massHistoLabelShort(vernum + string("mPiPiS1"));
0249
0250 string s2massHistoLabelLong(algo + string(" m_{#pi #pi}, R_{vtx} > 0.1, #chi^2 < 1"));
0251 string s2massHistoLabelShort(vernum + string("mPiPiS2"));
0252
0253 string s3massHistoLabelLong(algo + string(" m_{#pi #pi}, most cuts implemented"));
0254 string s3massHistoLabelShort(vernum + string("mPiPiS3"));
0255
0256 string s4massHistoLabelLong(algo + string(" m_{#pi #pi}, all cuts implemented"));
0257 string s4massHistoLabelShort(vernum + string("mPiPiS4"));
0258
0259 string vertexChi2HistoLabelLong(algo + string(" vertex #chi^{2}"));
0260 string vertexChi2HistoLabelShort(vernum + string("vtxChi2"));
0261
0262 string rVtxHistoLabelLong(algo + string(" r_{vtx} of V^{0} decay - radial"));
0263 string rVtxHistoLabelShort(vernum + string("rVtxRadial"));
0264
0265 string vtxSigHistoLabelLong(algo + string(" V^{0} vertex significance - radial"));
0266 string vtxSigHistoLabelShort(vernum + string("VtxRadialSig"));
0267
0268 string rVtxHistoLabel2Long(algo + string(" r_{vtx} of V^{0} decay after hit cut"));
0269 string rVtxHistoLabel2Short(vernum + string("rVtxAfterHitCut"));
0270
0271 string vtxSigHistoLabel2Long(algo + string(" V^{0} vertex significance"));
0272 string vtxSigHistoLabel2Short(vernum + string("VtxSphericalSig"));
0273
0274 string rErrorHistoLabelLong(algo + string(" Error in r_{vtx}"));
0275 string rErrorHistoLabelShort(vernum + string("#sigma R"));
0276
0277 string tkPtHistoLabelLong(algo + string(" Track P_{t}"));
0278 string tkPtHistoLabelShort(vernum + string("tkPt"));
0279
0280 string k0sPtHistoLabelLong(algo + string(" Sim K0s P_{t}"));
0281 string k0sPtHistoLabelShort(vernum + string("k0sPt"));
0282
0283 string tkChi2HistoLabelLong(algo + string(" Track #Chi^{2}"));
0284 string tkChi2HistoLabelShort(vernum + string("tkChi2"));
0285
0286 string sqrtTkChi2HistoLabelLong(algo + string(" sqrt of Track #Chi^{2}"));
0287 string sqrtTkChi2HistoLabelShort(vernum + string("sqrtTkChi2"));
0288
0289 string tkEtaHistoLabelLong(algo + string(" Track #eta"));
0290 string tkEtaHistoLabelShort(vernum + string("tkEta"));
0291
0292 string kShortEtaHistoLabelLong(algo + string(" K^{0}_{s} #eta"));
0293 string kShortEtaHistoLabelShort(vernum + string("k0sEta"));
0294
0295 string numHitsHistoLabelLong(algo + string(" Number Of Hits Per Track"));
0296 string numHitsHistoLabelShort(vernum + string("numHits"));
0297
0298 string simTkMpipiHistoLabelLong(algo + string(" Sim Track Pairs m_{#pi #pi}"));
0299 string simTkMpipiHistoLabelShort(vernum + string("simTkMpipi"));
0300
0301 simTkMpipiHisto = fs->make<TH1D>(simTkMpipiHistoLabelShort.c_str(), simTkMpipiHistoLabelLong.c_str(), 100, 0., 2.);
0302
0303 step1massHisto = fs->make<TH1D>(s1massHistoLabelShort.c_str(), s1massHistoLabelLong.c_str(), 100, 0., 2.);
0304 step2massHisto = fs->make<TH1D>(s2massHistoLabelShort.c_str(), s2massHistoLabelLong.c_str(), 100, 0., 2.);
0305 step3massHisto = fs->make<TH1D>(s3massHistoLabelShort.c_str(), s3massHistoLabelLong.c_str(), 100, 0., 2.);
0306 step4massHisto = fs->make<TH1D>(s4massHistoLabelShort.c_str(), s4massHistoLabelLong.c_str(), 100, 0., 2.);
0307
0308 vertexChi2Histo = fs->make<TH1D>(vertexChi2HistoLabelShort.c_str(), vertexChi2HistoLabelLong.c_str(), 100, 0., 20.);
0309 rVtxHisto1 = fs->make<TH1D>(rVtxHistoLabelShort.c_str(), rVtxHistoLabelLong.c_str(), 100, 0., 50.);
0310 vtxSigHisto1 = fs->make<TH1D>(vtxSigHistoLabelShort.c_str(), vtxSigHistoLabelLong.c_str(), 100, 0., 100.);
0311 rVtxHisto2 = fs->make<TH1D>(rVtxHistoLabel2Short.c_str(), rVtxHistoLabel2Long.c_str(), 100, 0., 50.);
0312 vtxSigHisto2 = fs->make<TH1D>(vtxSigHistoLabel2Short.c_str(), vtxSigHistoLabel2Long.c_str(), 100, 0., 100.);
0313 k0sPtHisto = fs->make<TH1D>(k0sPtHistoLabelShort.c_str(), k0sPtHistoLabelLong.c_str(), 100, 0., 20.);
0314
0315 rErrorHisto = fs->make<TH1D>(rErrorHistoLabelShort.c_str(), rErrorHistoLabelLong.c_str(), 1000, 0., 1.);
0316 tkPtHisto = fs->make<TH1D>(tkPtHistoLabelShort.c_str(), tkPtHistoLabelLong.c_str(), 100, 0., 20.);
0317 tkChi2Histo = fs->make<TH1D>(tkChi2HistoLabelShort.c_str(), tkChi2HistoLabelLong.c_str(), 100, 0., 20.);
0318 sqrtTkChi2Histo = fs->make<TH1D>(sqrtTkChi2HistoLabelShort.c_str(), sqrtTkChi2HistoLabelLong.c_str(), 100, 0., 20.);
0319 tkEtaHisto = fs->make<TH1D>(tkEtaHistoLabelShort.c_str(), tkEtaHistoLabelLong.c_str(), 100, -2.5, 2.5);
0320 kShortEtaHisto = fs->make<TH1D>(kShortEtaHistoLabelShort.c_str(), kShortEtaHistoLabelLong.c_str(), 100, -2.5, 2.5);
0321 numHitsHisto = fs->make<TH1D>(numHitsHistoLabelShort.c_str(), numHitsHistoLabelLong.c_str(), 20, 0., 20.);
0322
0323 myKshortMassHisto = fs->make<TH1D>(massHistoLabelShort.c_str(), massHistoLabelLong.c_str(), 100, 0.3, 0.7);
0324 myDecayRadiusHisto = fs->make<TH1D>(decRadHistoLabelShort.c_str(), decRadHistoLabelLong.c_str(), 100, 0., 40.);
0325 myRefitKshortMassHisto =
0326 fs->make<TH1D>((refit + massHistoLabelShort).c_str(), (refit + massHistoLabelLong).c_str(), 100, 0.3, 0.7);
0327 myRefitDecayRadiusHisto =
0328 fs->make<TH1D>((refit + decRadHistoLabelShort).c_str(), (refit + decRadHistoLabelLong).c_str(), 100, 0., 40.);
0329 myNativeParticleMassHisto =
0330 fs->make<TH1D>((native + massHistoLabelShort).c_str(), (refit + massHistoLabelLong).c_str(), 300, 0.3, 1.8);
0331 myEtaEfficiencyHisto = fs->make<TH1D>(etaEffHistoLabelShort.c_str(), etaEffHistoLabelLong.c_str(), 40, -2.5, 2.5);
0332 mySimEtaHisto = fs->make<TH1D>(simEtaHistoLabelShort.c_str(), simEtaHistoLabelLong.c_str(), 40, -2.5, 2.5);
0333 myRhoEfficiencyHisto = fs->make<TH1D>(rhoEffHistoLabelShort.c_str(), rhoEffHistoLabelLong.c_str(), 60, 0., 60.);
0334 myRhoEfficiencyHisto2 = fs->make<TH1D>(rhoEffHistoLabel2Short.c_str(), rhoEffHistoLabel2Long.c_str(), 60, 0., 60.);
0335 myRhoEfficiencyHisto3 = fs->make<TH1D>(rhoEffHistoLabel3Short.c_str(), rhoEffHistoLabel3Long.c_str(), 60, 0., 60.);
0336 myRhoEfficiencyHisto4 = fs->make<TH1D>(rhoEffHistoLabel4Short.c_str(), rhoEffHistoLabel4Long.c_str(), 60, 0., 60.);
0337 mySimRhoHisto = fs->make<TH1D>(simRhoHistoLabelShort.c_str(), simRhoHistoLabelLong.c_str(), 60, 0., 60.);
0338 myKshortPtHisto = fs->make<TH1D>(kShortPtHistoLabelShort.c_str(), kShortPtHistoLabelLong.c_str(), 100, 0., 100.);
0339 myImpactParameterHisto =
0340 fs->make<TH1D>(impactParamHistoLabelShort.c_str(), impactParamHistoLabelLong.c_str(), 100, 0., 10.);
0341 myImpactParameterHisto2 =
0342 fs->make<TH1D>(impactParamHisto2LabelShort.c_str(), impactParamHisto2LabelLong.c_str(), 100, 0., 10.);
0343 myNumSimKshortsHisto =
0344 fs->make<TH1D>(numSimKshortsHistoLabelShort.c_str(), numSimKshortsHistoLabelLong.c_str(), 100, 0., 100.);
0345 myNumRecoKshortsHisto =
0346 fs->make<TH1D>(numRecoKshortsHistoLabelShort.c_str(), numRecoKshortsHistoLabelLong.c_str(), 100, 0., 100.);
0347 myInnermostHitDistanceHisto = fs->make<TH1D>(
0348 innermostHitDistanceHistoLabelShort.c_str(), innermostHitDistanceHistoLabelLong.c_str(), 100, 0., 10.);
0349
0350 myEtaEfficiencyHisto->Sumw2();
0351 mySimEtaHisto->Sumw2();
0352 myRhoEfficiencyHisto->Sumw2();
0353 myRhoEfficiencyHisto2->Sumw2();
0354 myRhoEfficiencyHisto3->Sumw2();
0355 myRhoEfficiencyHisto4->Sumw2();
0356 mySimRhoHisto->Sumw2();
0357
0358 hitsOut.open("hitsOut.txt");
0359 }
0360
0361
0362 void V0Analyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0363
0364
0365 using namespace edm;
0366
0367 Handle<std::vector<reco::Vertex> > theVtxHandle;
0368
0369 Handle<reco::VertexCompositeCandidateCollection> theCandHand;
0370 Handle<SimTrackContainer> SimTk;
0371 Handle<SimVertexContainer> SimVtx;
0372 Handle<reco::TrackCollection> RecoTk;
0373
0374 ESHandle<MagneticField> bFieldHandle = iSetup.getHandle(m_mfToken);
0375 ESHandle<TrackerGeometry> trackerGeomHandle = iSetup.getHandle(m_geomToken);
0376
0377
0378
0379
0380 iEvent.getByLabel(algoLabel, V0CollectionName, theCandHand);
0381 iEvent.getByLabel(SimTkLabel, SimTk);
0382 iEvent.getByLabel(SimVtxLabel, SimVtx);
0383 iEvent.getByLabel(recoAlgoLabel, RecoTk);
0384
0385
0386 std::vector<reco::Track> theRecoTracks;
0387 std::vector<SimVertex> theSimVerts;
0388 std::vector<SimTrack> theSimTracks;
0389
0390 std::vector<reco::VertexCompositeCandidate> theKshorts;
0391 theRecoTracks.insert(theRecoTracks.end(), RecoTk->begin(), RecoTk->end());
0392 theSimVerts.insert(theSimVerts.end(), SimVtx->begin(), SimVtx->end());
0393 theSimTracks.insert(theSimTracks.end(), SimTk->begin(), SimTk->end());
0394 theKshorts.insert(theKshorts.end(), theCandHand->begin(), theCandHand->end());
0395 std::cout << theRecoTracks.size() << " " << theSimVerts.size() << " " << theSimTracks.size() << " "
0396 << theKshorts.size() << std::endl;
0397
0398
0399
0400
0401 double numSimKshorts = 0.;
0402 double simRad = 0.;
0403 for (unsigned int ndx1 = 0; ndx1 < theSimTracks.size(); ndx1++) {
0404 if (theSimTracks[ndx1].type() == 310) {
0405 math::XYZTLorentzVectorD k0sP(theSimTracks[ndx1].momentum());
0406 k0sPtHisto->Fill(sqrt(k0sP.Perp2()), 1.);
0407 simRad = sqrt(theSimVerts[(theSimTracks[ndx1].vertIndex() + 1)].position().perp2());
0408 numSimKshorts += 1.;
0409 }
0410 }
0411
0412 myNumSimKshortsHisto->Fill(numSimKshorts, 1.);
0413
0414
0415
0416
0417
0418 for (unsigned int stkidx1 = 0; stkidx1 < theSimTracks.size(); stkidx1++) {
0419 for (unsigned int stkidx2 = stkidx1 + 1; stkidx2 < theSimTracks.size(); stkidx2++) {
0420 SimTrack* thePosSimTk = 0;
0421 SimTrack* theNegSimTk = 0;
0422 if (theSimTracks[stkidx1].charge() > 0. && theSimTracks[stkidx2].charge() < 0.) {
0423 thePosSimTk = &theSimTracks[stkidx1];
0424 theNegSimTk = &theSimTracks[stkidx2];
0425 } else if (theSimTracks[stkidx1].charge() < 0. && theSimTracks[stkidx2].charge() > 0.) {
0426 theNegSimTk = &theSimTracks[stkidx1];
0427 thePosSimTk = &theSimTracks[stkidx2];
0428 }
0429
0430 if (thePosSimTk && theNegSimTk) {
0431 double posP2 = thePosSimTk->momentum().px() * thePosSimTk->momentum().px() +
0432 thePosSimTk->momentum().py() * thePosSimTk->momentum().py() +
0433 thePosSimTk->momentum().pz() * thePosSimTk->momentum().pz();
0434 double negP2 = theNegSimTk->momentum().px() * theNegSimTk->momentum().px() +
0435 theNegSimTk->momentum().py() * theNegSimTk->momentum().py() +
0436 theNegSimTk->momentum().pz() * theNegSimTk->momentum().pz();
0437 double posE = sqrt(posP2 + piMassSq);
0438 double negE = sqrt(negP2 + piMassSq);
0439 double k0sE = posE + negE;
0440 math::XYZTLorentzVectorD k0sMomentum(thePosSimTk->momentum() + theNegSimTk->momentum());
0441 k0sMomentum.SetE(k0sE);
0442 double k0sInvMass = k0sMomentum.M();
0443 simTkMpipiHisto->Fill(k0sInvMass, 1.);
0444 }
0445 thePosSimTk = theNegSimTk = 0;
0446 }
0447 }
0448
0449
0450
0451 double simCosTheta = theSimVerts[1].position().z() / theSimVerts[1].position().mag();
0452 double simEta = -log(tan(acos(simCosTheta) / 2.));
0453
0454 mySimRhoHisto->Fill(simRad, 1.);
0455 mySimEtaHisto->Fill(simEta, 1.);
0456
0457
0458
0459
0460
0461
0462 if (V0CollectionName == std::string("Kshort")) {
0463 myNumRecoKshortsHisto->Fill((double)theKshorts.size(), 1.);
0464 }
0465
0466
0467 for (unsigned int ksndx_ = 0; ksndx_ < theKshorts.size(); ksndx_++) {
0468 myNativeParticleMassHisto->Fill(theKshorts[ksndx_].mass());
0469 }
0470
0471 for (unsigned int ksndx = 0; ksndx < theKshorts.size(); ksndx++) {
0472 std::vector<reco::RecoChargedCandidate> v0daughters;
0473 std::vector<reco::TrackRef> theDaughterTracks;
0474
0475 for (unsigned int i = 0; i < theKshorts[ksndx].numberOfDaughters(); i++) {
0476 v0daughters.push_back(*(dynamic_cast<reco::RecoChargedCandidate*>(theKshorts[ksndx].daughter(i))));
0477 }
0478
0479 for (unsigned int j = 0; j < v0daughters.size(); j++) {
0480 theDaughterTracks.push_back(v0daughters[j].track());
0481 }
0482
0483 reco::TrackBase::Point beamSpot(0, 0, 0);
0484
0485 for (unsigned int k = 0; k < theDaughterTracks.size(); k++) {
0486 myImpactParameterHisto->Fill(sqrt(theDaughterTracks[k]->dxy(beamSpot) * theDaughterTracks[k]->dxy(beamSpot) +
0487 theDaughterTracks[k]->dsz(beamSpot) * theDaughterTracks[k]->dsz(beamSpot)),
0488 1.);
0489 }
0490
0491
0492 reco::Particle::Point k0s(theKshorts[ksndx].vx(), theKshorts[ksndx].vy(), theKshorts[ksndx].vz());
0493
0494
0495
0496
0497
0498 myKshortPtHisto->Fill(theKshorts[ksndx].pt(), 1.);
0499
0500
0501
0502
0503
0504
0505
0506
0507
0508
0509
0510
0511
0512
0513
0514
0515 kShortEtaHisto->Fill(theKshorts[ksndx].eta(), 1.);
0516
0517
0518
0519 double decayRad = sqrt(k0s.x() * k0s.x() + k0s.y() * k0s.y());
0520 myDecayRadiusHisto->Fill(decayRad);
0521
0522 GlobalPoint vtxPos(k0s.x(), k0s.y(), k0s.z());
0523 double recoRad = vtxPos.perp();
0524 double recoCosTheta = vtxPos.z() / vtxPos.mag();
0525 double recoEta = -log(tan(acos(recoCosTheta) / 2.));
0526
0527 double x_ = k0s.x();
0528 double y_ = k0s.y();
0529 double z_ = k0s.z();
0530
0531
0532
0533
0534
0535
0536 double sig00 = theKshorts[ksndx].vertexCovariance(0, 0);
0537 double sig11 = theKshorts[ksndx].vertexCovariance(1, 1);
0538 double sig22 = theKshorts[ksndx].vertexCovariance(2, 2);
0539 double sig01 = theKshorts[ksndx].vertexCovariance(0, 1);
0540 double sig02 = theKshorts[ksndx].vertexCovariance(0, 2);
0541 double sig12 = theKshorts[ksndx].vertexCovariance(1, 2);
0542
0543 double vtxRSph = vtxPos.mag();
0544 double vtxR = vtxPos.perp();
0545 double vtxChi2 = theKshorts[ksndx].vertexChi2();
0546 double vtxErrorSph = sqrt(sig00 * (x_ * x_) + sig11 * (y_ * y_) + sig22 * (z_ * z_) +
0547 2 * (sig01 * (x_ * y_) + sig02 * (x_ * z_) + sig12 * (y_ * z_))) /
0548 vtxRSph;
0549 double vtxError = sqrt(sig00 * (x_ * x_) + sig11 * (y_ * y_) + 2 * sig01 * (x_ * y_)) / vtxR;
0550 double vtxSigSph = vtxRSph / vtxErrorSph;
0551 double vtxSig = vtxR / vtxError;
0552
0553 rErrorHisto->Fill(vtxError, 1.);
0554
0555 using namespace reco;
0556 std::vector<reco::TrackRef> theVtxTrax;
0557 for (unsigned int i = 0; i < v0daughters.size(); i++) {
0558 theVtxTrax.push_back(v0daughters[i].track());
0559 }
0560
0561 bool hitsOkay2 = true;
0562 bool tkChi2Cut = true;
0563 if (theVtxTrax.size() == 2) {
0564 if (theVtxTrax[0]->normalizedChi2() > 5. || theVtxTrax[1]->normalizedChi2() > 5.) {
0565 tkChi2Cut = false;
0566 }
0567
0568 GlobalPoint tk1hitPos(
0569 theVtxTrax[0]->innerPosition().x(), theVtxTrax[0]->innerPosition().y(), theVtxTrax[0]->innerPosition().z());
0570 GlobalPoint tk2hitPos(
0571 theVtxTrax[1]->innerPosition().x(), theVtxTrax[1]->innerPosition().y(), theVtxTrax[1]->innerPosition().z());
0572
0573 if (tk1hitPos.mag() < (vtxRSph - 4. * vtxErrorSph) && theVtxTrax[0]->innerOk()) {
0574 hitsOkay2 = false;
0575 }
0576 if (tk2hitPos.mag() < (vtxRSph - 4. * vtxErrorSph) && theVtxTrax[1]->innerOk()) {
0577 hitsOkay2 = false;
0578 }
0579 }
0580
0581 bool hitsOkay = true;
0582 bool nHitsCut = true;
0583
0584
0585
0586 if (theVtxTrax.size() == 2) {
0587 if (theVtxTrax[0]->recHitsSize() && theVtxTrax[1]->recHitsSize()) {
0588 double nHits1 = (double)theVtxTrax[0]->numberOfValidHits();
0589 double nHits2 = (double)theVtxTrax[1]->numberOfValidHits();
0590
0591
0592
0593
0594
0595
0596
0597
0598
0599
0600
0601
0602
0603
0604
0605
0606
0607
0608
0609
0610
0611
0612
0613
0614
0615
0616
0617
0618
0619
0620
0621
0622
0623
0624
0625
0626
0627
0628
0629
0630
0631
0632 numHitsHisto->Fill(nHits1, 1.);
0633 numHitsHisto->Fill(nHits2, 1.);
0634 if (nHits1 < 8. || nHits2 < 8.) {
0635 nHitsCut = false;
0636 }
0637 }
0638 }
0639
0640
0641
0642 hitsOut << "hitsOkay=" << hitsOkay << ", hitsOkay2=" << hitsOkay2 << std::endl;
0643 if (!hitsOkay)
0644 ++numDiff1;
0645 if (!hitsOkay2)
0646 ++numDiff2;
0647 for (unsigned int ndx3 = 0; ndx3 < theRecoTracks.size(); ndx3++) {
0648 tkEtaHisto->Fill(theRecoTracks[ndx3].eta(), 1.);
0649
0650 }
0651
0652 if (nHitsCut) {
0653 for (unsigned int ndx3_1 = 0; ndx3_1 < theRecoTracks.size(); ndx3_1++) {
0654
0655 tkChi2Histo->Fill(theRecoTracks[ndx3_1].normalizedChi2(), 1.);
0656 }
0657 }
0658
0659 vertexChi2Histo->Fill(vtxChi2, 1.);
0660
0661
0662
0663 if (vtxChi2 < 7. && nHitsCut && tkChi2Cut) {
0664 step1massHisto->Fill(theKshorts[ksndx].mass(), 1.);
0665 myRhoEfficiencyHisto->Fill(recoRad, 1.);
0666 rVtxHisto1->Fill(vtxR, 1.);
0667
0668
0669 step2massHisto->Fill(theKshorts[ksndx].mass(), 1.);
0670 myRhoEfficiencyHisto2->Fill(recoRad, 1.);
0671 vtxSigHisto1->Fill(vtxSig, 1.);
0672 if (vtxSig > 20.) {
0673 step3massHisto->Fill(theKshorts[ksndx].mass(), 1.);
0674 myRhoEfficiencyHisto3->Fill(recoRad, 1.);
0675 if (hitsOkay) {
0676 step4massHisto->Fill(theKshorts[ksndx].mass(), 1.);
0677 rVtxHisto2->Fill(vtxR, 1.);
0678 myEtaEfficiencyHisto->Fill(recoEta, 1.);
0679 myRhoEfficiencyHisto4->Fill(recoRad, 1.);
0680 }
0681 }
0682
0683 }
0684
0685
0686 if (vtxChi2 < 1.) {
0687
0688 if (vtxRSph > 0.1) {
0689 vtxSigHisto2->Fill(vtxSigSph, 1.);
0690 }
0691 }
0692
0693
0694
0695
0696
0697
0698
0699
0700 std::vector<reco::Track> theRefTracks;
0701 theRefTracks.push_back(*theDaughterTracks[0]);
0702 theRefTracks.push_back(*theDaughterTracks[1]);
0703
0704
0705
0706
0707
0708
0709
0710
0711
0712
0713
0714
0715
0716
0717
0718
0719
0720
0721
0722 reco::TransientTrack refTemp1(theRefTracks[0], &(*bFieldHandle));
0723 reco::TransientTrack refTemp2(theRefTracks[1], &(*bFieldHandle));
0724
0725
0726
0727
0728
0729
0730
0731
0732 reco::TransientTrack temp1(*theDaughterTracks[0], &(*bFieldHandle));
0733 reco::TransientTrack temp2(*theDaughterTracks[1], &(*bFieldHandle));
0734
0735 math::XYZPoint tk1InnerHitPos = temp1.track().innerPosition();
0736 math::XYZPoint tk2InnerHitPos = temp2.track().innerPosition();
0737 math::XYZVector difference = tk2InnerHitPos - tk1InnerHitPos;
0738 double dist =
0739 sqrt(difference.x() * difference.x() + difference.y() * difference.y() + difference.z() * difference.z());
0740 myInnermostHitDistanceHisto->Fill(dist, 1.);
0741
0742 TrajectoryStateClosestToBeamLine tscb1(temp1.stateAtBeamLine());
0743 TrajectoryStateClosestToBeamLine tscb2(temp2.stateAtBeamLine());
0744 myImpactParameterHisto2->Fill(tscb1.transverseImpactParameter().value(), 1.);
0745 myImpactParameterHisto2->Fill(tscb2.transverseImpactParameter().value(), 1.);
0746
0747 TrajectoryStateClosestToPoint tscpRef1(refTemp1.trajectoryStateClosestToPoint(vtxPos));
0748 TrajectoryStateClosestToPoint tscpRef2(refTemp2.trajectoryStateClosestToPoint(vtxPos));
0749 TrajectoryStateClosestToPoint tscp1(temp1.trajectoryStateClosestToPoint(vtxPos));
0750 TrajectoryStateClosestToPoint tscp2(temp2.trajectoryStateClosestToPoint(vtxPos));
0751
0752 GlobalVector refTrack1P(tscpRef1.momentum());
0753 GlobalVector refTrack2P(tscpRef2.momentum());
0754 GlobalVector track1P(tscp1.momentum());
0755 GlobalVector track2P(tscp2.momentum());
0756
0757
0758 GlobalVector refTotalP(refTrack1P + refTrack2P);
0759 GlobalVector totalP(track1P + track2P);
0760 double refPTotMag =
0761 sqrt(refTotalP.x() * refTotalP.x() + refTotalP.y() * refTotalP.y() + refTotalP.z() * refTotalP.z());
0762 double pTotMag = sqrt(totalP.x() * totalP.x() + totalP.y() * totalP.y() + totalP.z() * totalP.z());
0763
0764
0765 double refETot = sqrt(refTrack1P.x() * refTrack1P.x() + refTrack1P.y() * refTrack1P.y() +
0766 refTrack1P.z() * refTrack1P.z() + piMassSq) +
0767 sqrt(refTrack2P.x() * refTrack2P.x() + refTrack2P.y() * refTrack2P.y() +
0768 refTrack2P.z() * refTrack2P.z() + piMassSq);
0769 double eTot = sqrt(track1P.x() * track1P.x() + track1P.y() * track1P.y() + track1P.z() * track1P.z() + piMassSq) +
0770 sqrt(track2P.x() * track2P.x() + track2P.y() * track2P.y() + track2P.z() * track2P.z() + piMassSq);
0771
0772
0773 double refKShortMass = sqrt((refETot + refPTotMag) * (refETot - refPTotMag));
0774 double kShortMass = sqrt((eTot + pTotMag) * (eTot - pTotMag));
0775
0776 myRefitKshortMassHisto->Fill(refKShortMass);
0777 myKshortMassHisto->Fill(kShortMass);
0778 }
0779
0780 }
0781
0782
0783 void V0Analyzer::endJob() {
0784 static int endJcount = 0;
0785
0786 myEtaEfficiencyHisto->Divide(mySimEtaHisto);
0787 myRhoEfficiencyHisto->Divide(mySimRhoHisto);
0788 myRhoEfficiencyHisto2->Divide(mySimRhoHisto);
0789 myRhoEfficiencyHisto3->Divide(mySimRhoHisto);
0790 myRhoEfficiencyHisto4->Divide(mySimRhoHisto);
0791
0792
0793
0794
0795
0796
0797
0798
0799
0800
0801
0802
0803
0804
0805
0806
0807
0808
0809
0810
0811
0812
0813
0814
0815
0816
0817
0818
0819
0820
0821
0822
0823
0824
0825
0826
0827
0828
0829
0830
0831
0832
0833
0834
0835
0836
0837
0838
0839
0840
0841
0842
0843 endJcount++;
0844 std::cout << "ENDJCOUNT: " << endJcount << std::endl;
0845
0846
0847
0848
0849 hitsOut << "numDiff1 = " << numDiff1 << ", numDiff2 = " << numDiff2 << std::endl;
0850
0851 hitsOut.close();
0852 }
0853
0854
0855 #include "FWCore/PluginManager/interface/ModuleDef.h"
0856
0857 DEFINE_FWK_MODULE(V0Analyzer);