V0Producer/test/V0Analyzer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    V0Analyzer
0004 // Class:      V0Analyzer
0005 //
0006 /**\class V0Analyzer V0Analyzer.cc RecoVertex/V0Producer/test/V0Analyzer.cc
0007
0008  Description: <one line class summary>
0009
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Brian Drell
0015 //         Created:  Tue May 22 23:54:16 CEST 2007
0016 //
0017 //
0018
0019 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateClosestToPoint.h"
0020
0021 // system include files
0022 #include <memory>
0023
0024 // user include files
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 //#include "DataFormats/TrackReco/interface/Track.h"
0035 //#include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0036 //#include "DataFormats/TrackReco/interface/TrackBase.h"
0037 #include "DataFormats/VertexReco/interface/Vertex.h"
0038
0039 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0040 #include "MagneticField/VolumeBasedEngine/interface/VolumeBasedMagneticField.h"
0041 //#include "DataFormats/V0Candidate/interface/V0Candidate.h"
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 // class declaration
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   //virtual void beginJob() ;
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   //TFile* theHistoFile;
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   // Histograms for figuring out the best cuts
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   //  TH1D* simRHisto;
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   // ----------member data ---------------------------
0140 };
0141
0142 //
0143 // constants, enums and typedefs
0144 //
0145
0146 const double piMassSq = 0.019479101;
0147
0148 //
0149 // static data member definitions
0150 //
0151
0152 //
0153 // constructors and destructor
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   //now do what ever initialization is needed
0166   //theHistoFile = 0;
0167
0168   usesResource(TFileService::kSharedResource);
0169
0170   numDiff1 = numDiff2 = 0;
0171 }
0172
0173 V0Analyzer::~V0Analyzer() {
0174   // do anything here that needs to be done at desctruction time
0175   // (e.g. close files, deallocate resources etc.)
0176 }
0177
0178 //
0179 // member functions
0180 //
0181
0182 // ------------ method called once each job just before starting event loop  ------------
0183 //void V0Analyzer::beginJob() {
0184 void V0Analyzer::beginJob() {
0185   using std::string;
0186
0187   edm::Service<TFileService> fs;
0188   //theHistoFile = new TFile(HistoFileName.c_str(), "RECREATE");
0189   //theHistoFile->Open();
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 // ------------ method called to for each event  ------------
0362 void V0Analyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0363   //std::cout << std::endl << "@@@In module..." << std::endl;
0364   //std::cout << "Creating HANDLES..." << std::endl;
0365   using namespace edm;
0366   //Handle<reco::VertexCollection> theVtxHandle;
0367   Handle<std::vector<reco::Vertex> > theVtxHandle;
0368   //Handle< std::vector<reco::V0Candidate> > theCandHand;
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   //const TrackerGeometry* trackerGeom = trackerGeomHandle.product();
0378
0379   //std::cout << "Getting by label..." << std::endl;
0380   iEvent.getByLabel(algoLabel, V0CollectionName, theCandHand);
0381   iEvent.getByLabel(SimTkLabel, SimTk);
0382   iEvent.getByLabel(SimVtxLabel, SimVtx);
0383   iEvent.getByLabel(recoAlgoLabel, RecoTk);
0384
0385   //std::cout << "Creating and filling vectors..." << std::endl;
0386   std::vector<reco::Track> theRecoTracks;
0387   std::vector<SimVertex> theSimVerts;
0388   std::vector<SimTrack> theSimTracks;
0389   //std::vector<reco::V0Candidate> theKshorts;
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   // Done getting and filling
0399
0400   // Count how many simulated K0s we have
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   //std::cout << "Found " << theKshorts.size() << " K0s events."
0415   // << std::endl;
0416
0417   // loop over sim tracks and calculate m_pipi, and histogram it
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   // Calculate sim k0s parameters from the sim decay vertex
0450   //double simRad = sqrt(theSimVerts[1].position().perp2());
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   //bool hasKshort = false;
0458   //bool hasLambda = false;
0459   //bool hasLambdaBar = false;
0460   //std::cout << theKshorts.size() << std::endl;
0461
0462   if (V0CollectionName == std::string("Kshort")) {
0463     myNumRecoKshortsHisto->Fill((double)theKshorts.size(), 1.);
0464   }
0465
0466   // Histogram the invariant mass (retrieved from the V0Candidate)
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     //reco::Vertex k0s = theKshorts[ksndx].vertex();
0492     reco::Particle::Point k0s(theKshorts[ksndx].vx(), theKshorts[ksndx].vy(), theKshorts[ksndx].vz());
0493
0494     /*myKshortPtHisto->Fill(sqrt(theKshorts[ksndx].px() *
0495                    theKshorts[ksndx].px() +
0496                    theKshorts[ksndx].py() *
0497                    theKshorts[ksndx].py()), 1.);*/
0498     myKshortPtHisto->Fill(theKshorts[ksndx].pt(), 1.);
0499
0500     /*
0501     if(theKshorts[ksndx].pdgId() == 310) {
0502       hasKshort = true;
0503     }
0504     if(theKshorts[ksndx].pdgId() == 3122) {
0505       hasLambda = true;
0506     }
0507     if(theKshorts[ksndx].pdgId() == -3122) {
0508       hasLambdaBar = true;
0509     }
0510     std::cout << "$#@#$: " << hasKshort << " " << hasLambda << " "
0511               << hasLambdaBar << " size: " << theKshorts.size() << std::endl;
0512     std::cout << "@@@ MASS: " << theKshorts[ksndx].mass() << std::endl;
0513     */
0514
0515     kShortEtaHisto->Fill(theKshorts[ksndx].eta(), 1.);
0516
0517     //std::cout << "tracksSize()=" << k0s.tracksSize() << std::endl;
0518     //std::cout << "hasRefittedTracks()="<< k0s.hasRefittedTracks() << std::endl;
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     //double sig00 = k0s.covariance(0,0);
0531     //double sig11 = k0s.covariance(1,1);
0532     //double sig22 = k0s.covariance(2,2);
0533     //double sig01 = k0s.covariance(0,1);
0534     //double sig02 = k0s.covariance(0,2);
0535     //double sig12 = k0s.covariance(1,2);
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       //if( tk1hitPos.perp() < (vtxR - 5.*vtxError)
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     //std::cout << "theVtxTrax.size = " << theVtxTrax.size() << std::endl;
0584
0585     //hitsOut << "theVtxTrax.size = " << theVtxTrax.size() << std::endl;
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     double nHits1 = 0.;
0593     double nHits2 = 0.;
0594     for(auto const& tk1Hit : theVtxTrax[0]->recHits()) {
0595       const TrackingRecHit* tk1HitPtr = tk1Hit.get();
0596       if( tk1Hit->isValid() ) nHits1 += 1.;
0597       if( tk1Hit->isValid() && hitsOkay) {
0598         GlobalPoint tk1HitPosition
0599           = trackerGeom->idToDet(tk1HitPtr->
0600                      geographicalId())->
0601           surface().toGlobal(tk1HitPtr->localPosition());
0602         //std::cout << typeid(*tk1HitPtr).name();<--This is how
0603         //                               we can access the hit type.
0604
0605         //std::cout << ":::" << tk1HitPosition.mag() << ", "
0606         //<< vtxRSph - 4.*vtxErrorSph << std::endl;
0607
0608         if( tk1HitPosition.mag() < (vtxRSph - 4.*vtxErrorSph) ) {
0609           hitsOkay = false;
0610           //std::cout << "Flagged on track 1." << std::endl;
0611         }
0612       }
0613     }
0614
0615     for(auto const& tk2Hit : theVtxTrax[1]->recHits()) {
0616       const TrackingRecHit* tk2HitPtr = tk2Hit.get();
0617       if( tk2Hit->isValid() ) nHits2 += 1.;
0618       if( tk2Hit->isValid() && hitsOkay) {
0619         GlobalPoint tk2HitPosition
0620           = trackerGeom->idToDet(tk2HitPtr->
0621                      geographicalId())->
0622           surface().toGlobal(tk2HitPtr->localPosition());
0623         //std::cout << ":::" << tk2HitPosition.mag() << ", "
0624         //<< vtxRSph - 4.*vtxErrorSph << std::endl;
0625
0626         if( tk2HitPosition.mag() < (vtxRSph - 4.*vtxErrorSph) ) {
0627           hitsOkay = false;
0628           //std::cout << "Flagged on track 2." << std::endl;
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     //std::cout << "hitsOkay=" << hitsOkay << ", hitsOkay2=" << hitsOkay2
0641     //<< std::endl;
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       //tkChi2Histo->Fill(theRecoTracks[ndx3].normalizedChi2(), 1.);
0650     }
0651
0652     if (nHitsCut) {
0653       for (unsigned int ndx3_1 = 0; ndx3_1 < theRecoTracks.size(); ndx3_1++) {
0654         //tkEtaHisto->Fill(theRecoTracks[ndx3_1].eta(), 1.);
0655         tkChi2Histo->Fill(theRecoTracks[ndx3_1].normalizedChi2(), 1.);
0656       }
0657     }
0658
0659     vertexChi2Histo->Fill(vtxChi2, 1.);
0660
0661     //step1massHisto->Fill(theKshorts[ksndx].mass(), 1.);
0662     //if(hitsOkay2) {
0663     if (vtxChi2 < 7. && nHitsCut && tkChi2Cut) {
0664       step1massHisto->Fill(theKshorts[ksndx].mass(), 1.);
0665       myRhoEfficiencyHisto->Fill(recoRad, 1.);
0666       rVtxHisto1->Fill(vtxR, 1.);
0667       //if(vtxR > 0.1) {
0668       //if(vtxR > 1.) {
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       //rVtxHisto2->Fill(vtxRSph, 1.);
0688       if (vtxRSph > 0.1) {
0689         vtxSigHisto2->Fill(vtxSigSph, 1.);
0690       }
0691     }
0692
0693     //theKshorts[ksndx]
0694
0695     //if(theKshorts.size() < 2) {
0696
0697     //}
0698
0699     //    std::vector<reco::Track> theRefTracks = k0s.refittedTracks();
0700     std::vector<reco::Track> theRefTracks;
0701     theRefTracks.push_back(*theDaughterTracks[0]);
0702     theRefTracks.push_back(*theDaughterTracks[1]);
0703
0704     /*for (unsigned int ndx = 0; ndx < theRefTracks.size(); ndx++) {
0705
0706       std::cout << "IN LOOP, THIS ISN'T THE PROBLEM. "
0707         << theRefTracks.size() << ": " << std::endl;
0708
0709       //reco::Track tmpTk( *(k0s.originalTrack( theRefTracks[ndx] )) );
0710       reco::TrackBaseRef tmpTkRef = k0s.originalTrack( theRefTracks[ndx] );
0711       std::cout << "Got the TrackBaseRef." << std::endl;
0712       reco::TrackRef tmpTkRef2 = tmpTkRef.castTo<reco::TrackRef>();
0713       std::cout << "Created the TrackRef." << std::endl;
0714       reco::Track tmpTk( *(tmpTkRef2) );
0715       std::cout << "Did I even get here???" << std::endl;
0716       myImpactParameterHisto->Fill( sqrt( tmpTk.dxy( beamSpot )
0717                       * tmpTk.dxy( beamSpot )
0718                       + tmpTk.dsz( beamSpot )
0719                       * tmpTk.dsz( beamSpot ) ), 1.);
0720                       }*/
0721
0722     reco::TransientTrack refTemp1(theRefTracks[0], &(*bFieldHandle));
0723     reco::TransientTrack refTemp2(theRefTracks[1], &(*bFieldHandle));
0724
0725     /*reco::TrackBaseRef tkBRef1 = *k0s.tracks_begin();
0726     reco::TrackBaseRef tkBRef2 = *(++k0s.tracks_begin());//tkBIt2;
0727
0728     reco::TransientTrack temp1( *( tkBRef1.castTo<reco::TrackRef>() ),
0729                 &( *bFieldHandle ) );
0730     reco::TransientTrack temp2( *( tkBRef2.castTo<reco::TrackRef>() ),
0731     &( *bFieldHandle ) );*/
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     // Find the total momentum of the tracks and its magnitude
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     // Find the total energy of the tracks
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     // Calculate the invariant mass
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   //std::cout << "At end of analyze() function" << std::endl;
0780 }
0781
0782 // ------------ method called once each job just after ending the event loop  ------------
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   /*theHistoFile->cd();
0793
0794   myKshortMassHisto->Write();
0795   myDecayRadiusHisto->Write();
0796   myRefitKshortMassHisto->Write();
0797   myRefitDecayRadiusHisto->Write();
0798   myNativeParticleMassHisto->Write();
0799   myEtaEfficiencyHisto->Write();
0800   mySimEtaHisto->Write();
0801   myRhoEfficiencyHisto->Write();
0802   myRhoEfficiencyHisto2->Write();
0803   myRhoEfficiencyHisto3->Write();
0804   myRhoEfficiencyHisto4->Write();
0805   mySimRhoHisto->Write();
0806   myKshortPtHisto->Write();
0807   myImpactParameterHisto->Write();
0808   myImpactParameterHisto2->Write();
0809   myNumSimKshortsHisto->Write();
0810   myNumRecoKshortsHisto->Write();
0811   myInnermostHitDistanceHisto->Write();
0812
0813   // Histograms for figuring out the best cuts
0814   step1massHisto->Write();
0815   step2massHisto->Write();
0816   step3massHisto->Write();
0817   step4massHisto->Write();
0818
0819   vertexChi2Histo->Write();
0820   rVtxHisto1->Write();
0821   vtxSigHisto1->Write();
0822   rVtxHisto2->Write();
0823   //simRHisto->Write();
0824   vtxSigHisto2->Write();
0825
0826   rErrorHisto->Write();
0827   tkPtHisto->Write();
0828   k0sPtHisto->Write();
0829   tkChi2Histo->Write();
0830   sqrtTkChi2Histo->Write();
0831   tkEtaHisto->Write();
0832   kShortEtaHisto->Write();
0833   numHitsHisto->Write();
0834
0835   std::cout << "Writing out histogram file." << std::endl;
0836
0837   simTkMpipiHisto->Write();
0838   //std::cout << "Pointer address: " << theHistoFile << std::endl;
0839   theHistoFile->Write();
0840   theHistoFile->Close();
0841   delete theHistoFile;
0842   theHistoFile=0;*/
0843   endJcount++;
0844   std::cout << "ENDJCOUNT: " << endJcount << std::endl;
0845
0846   //std::cout << "numDiff1 = " << numDiff1 << ", numDiff2 = "
0847   //<< numDiff2 << std::endl;
0848
0849   hitsOut << "numDiff1 = " << numDiff1 << ", numDiff2 = " << numDiff2 << std::endl;
0850
0851   hitsOut.close();
0852 }
0853
0854 //define this as a plug-in
0855 #include "FWCore/PluginManager/interface/ModuleDef.h"
0856
0857 DEFINE_FWK_MODULE(V0Analyzer);