MomentumScaleCalibration/bin/ZntupleToTreeConverter.cc

0001 #include <memory>
0002 #include <string>
0003 #include <vector>
0004 #include <sstream>
0005 #include <fstream>
0006 #include <iostream>
0007
0008 #include <TH1F.h>
0009 #include <TROOT.h>
0010 #include <TFile.h>
0011 #include <TSystem.h>
0012
0013 #include "DataFormats/Common/interface/Handle.h"
0014 #include "DataFormats/FWLite/interface/Event.h"
0015 #include "DataFormats/PatCandidates/interface/Muon.h"
0016 #include "FWCore/FWLite/interface/FWLiteEnabler.h"
0017 #include "DataFormats/FWLite/interface/Handle.h"
0018 #include "DataFormats/Common/interface/Handle.h"
0019 #include "PhysicsTools/FWLite/interface/TFileService.h"
0020 #include "MuonAnalysis/MomentumScaleCalibration/interface/RootTreeHandler.h"
0021
0022 // Useful function to convert 4-vector coordinates
0023 // -----------------------------------------------
0024 lorentzVector fromPtEtaPhiToPxPyPz(const double* ptEtaPhiE) {
0025   double muMass = 0.105658;
0026   double px = ptEtaPhiE[0] * cos(ptEtaPhiE[2]);
0027   double py = ptEtaPhiE[0] * sin(ptEtaPhiE[2]);
0028   double tmp = 2 * atan(exp(-ptEtaPhiE[1]));
0029   double pz = ptEtaPhiE[0] * cos(tmp) / sin(tmp);
0030   double E = sqrt(px * px + py * py + pz * pz + muMass * muMass);
0031
0032   // lorentzVector corrMu(px,py,pz,E);
0033   // To fix memory leaks, this is to be substituted with
0034   // std::unique_ptr<lorentzVector> corrMu(new lorentzVector(px, py, pz, E));
0035
0036   return lorentzVector(px, py, pz, E);
0037 }
0038
0039 int main(int argc, char* argv[]) {
0040   if (argc != 2) {
0041     std::cout << "Please provide the name of the file with file: or rfio: as needed" << std::endl;
0042     exit(1);
0043   }
0044   std::string fileName(argv[1]);
0045   if (fileName.find("file:") != 0 && fileName.find("rfio:") != 0) {
0046     std::cout << "Please provide the name of the file with file: or rfio: as needed" << std::endl;
0047     exit(1);
0048   }
0049
0050   // ----------------------------------------------------------------------
0051   // First Part:
0052   //
0053   //  * enable FWLite
0054   //  * book the histograms of interest
0055   //  * open the input file
0056   // ----------------------------------------------------------------------
0057
0058   // load framework libraries
0059   gSystem->Load("libFWCoreFWLite");
0060   FWLiteEnabler::enable();
0061
0062   // book a set of histograms
0063   fwlite::TFileService fs = fwlite::TFileService("analyzeBasics.root");
0064   TFileDirectory theDir = fs.mkdir("analyzeBasic");
0065   TH1F* muonPt_ = theDir.make<TH1F>("muonPt", "pt", 100, 0., 300.);
0066   TH1F* muonEta_ = theDir.make<TH1F>("muonEta", "eta", 100, -3., 3.);
0067   TH1F* muonPhi_ = theDir.make<TH1F>("muonPhi", "phi", 100, -5., 5.);
0068
0069   // open input file (can be located on castor)
0070   TFile* inFile = TFile::Open(fileName.c_str());
0071
0072   //   TFile* inFile = TFile::Open("rfio:/castor/cern.ch/user/f/fabozzi/36XSkimData/run_139791-140159/NtupleLoose_139791-140159_v2.root");
0073   //   TFile* inFile = TFile::Open("rfio:/castor/cern.ch/user/f/fabozzi/36XSkimData/run_140160-140182/NtupleLoose_140160-140182.root");
0074   //   TFile* inFile = TFile::Open("rfio:/castor/cern.ch/user/f/fabozzi/36XSkimData/run_140183-140399/NtupleLoose_140183-140399.root");
0075   //   TFile* inFile = TFile::Open("rfio:/castor/cern.ch/user/d/degrutto/36XSkimData/run_140440-141961/NtupleLoose_140440-141961.root");
0076   //   TFile* inFile = TFile::Open("rfio:/castor/cern.ch/user/d/degrutto/36XSkimData/run_142035-142664/NtupleLoose_142035-142664.root");
0077
0078   // ----------------------------------------------------------------------
0079   // Second Part:
0080   //
0081   //  * loop the events in the input file
0082   //  * receive the collections of interest via fwlite::Handle
0083   //  * fill the histograms
0084   //  * after the loop close the input file
0085   // ----------------------------------------------------------------------
0086
0087   // Create the RootTreeHandler to save the events in the root tree
0088   RootTreeHandler treeHandler;
0089   // MuonPairVector pairVector;
0090   std::vector<MuonPair> pairVector;
0091
0092   // loop the events
0093   unsigned int iEvent = 0;
0094   fwlite::Event ev(inFile);
0095   for (ev.toBegin(); !ev.atEnd(); ++ev, ++iEvent) {
0096     // simple event counter
0097     if (iEvent > 0 && iEvent % 100 == 0) {
0098       std::cout << "  processing event: " << iEvent << std::endl;
0099     }
0100
0101     // Handle to the muon collection
0102     fwlite::Handle<std::vector<float> > muon1pt;
0103     fwlite::Handle<std::vector<float> > muon1eta;
0104     fwlite::Handle<std::vector<float> > muon1phi;
0105     fwlite::Handle<std::vector<float> > muon2pt;
0106     fwlite::Handle<std::vector<float> > muon2eta;
0107     fwlite::Handle<std::vector<float> > muon2phi;
0108     muon1pt.getByLabel(ev, "goodZToMuMuEdmNtupleLoose", "zGoldenDau1Pt");
0109     muon1eta.getByLabel(ev, "goodZToMuMuEdmNtupleLoose", "zGoldenDau1Eta");
0110     muon1phi.getByLabel(ev, "goodZToMuMuEdmNtupleLoose", "zGoldenDau1Phi");
0111     muon2pt.getByLabel(ev, "goodZToMuMuEdmNtupleLoose", "zGoldenDau2Pt");
0112     muon2eta.getByLabel(ev, "goodZToMuMuEdmNtupleLoose", "zGoldenDau2Eta");
0113     muon2phi.getByLabel(ev, "goodZToMuMuEdmNtupleLoose", "zGoldenDau2Phi");
0114
0115     if (!muon1pt.isValid())
0116       continue;
0117     if (!muon1eta.isValid())
0118       continue;
0119     if (!muon1phi.isValid())
0120       continue;
0121     if (!muon2pt.isValid())
0122       continue;
0123     if (!muon2eta.isValid())
0124       continue;
0125     if (!muon2phi.isValid())
0126       continue;
0127     // std::cout << "muon1pt = " << muon1pt->size() << std::endl;
0128
0129     // loop muon collection and fill histograms
0130     if (muon1pt->size() != muon2pt->size()) {
0131       std::cout << "Error: size of muon1 and muon2 is different. Skipping event" << std::endl;
0132       continue;
0133     }
0134     for (unsigned i = 0; i < muon1pt->size(); ++i) {
0135       muonPt_->Fill((*muon1pt)[i]);
0136       muonEta_->Fill((*muon1eta)[i]);
0137       muonPhi_->Fill((*muon1phi)[i]);
0138       muonPt_->Fill((*muon2pt)[i]);
0139       muonEta_->Fill((*muon2eta)[i]);
0140       muonPhi_->Fill((*muon2phi)[i]);
0141
0142       double muon1[3] = {(*muon1pt)[i], (*muon1eta)[i], (*muon1phi)[i]};
0143       double muon2[3] = {(*muon2pt)[i], (*muon2eta)[i], (*muon2phi)[i]};
0144
0145       // pairVector.push_back( std::make_pair( fromPtEtaPhiToPxPyPz(muon1), fromPtEtaPhiToPxPyPz(muon2) ) );
0146       pairVector.push_back(
0147           MuonPair(fromPtEtaPhiToPxPyPz(muon1), fromPtEtaPhiToPxPyPz(muon2), MuScleFitEvent(0, 0, 0, 0, 0, 0)));
0148     }
0149   }
0150   size_t namePos = fileName.find_last_of('/');
0151   treeHandler.writeTree(("tree_" + fileName.substr(namePos + 1, fileName.size())).c_str(), &pairVector);
0152
0153   // close input file
0154   inFile->Close();
0155
0156   // ----------------------------------------------------------------------
0157   // Third Part:
0158   //
0159   //  * never forget to free the memory of objects you created
0160   // ----------------------------------------------------------------------
0161
0162   // in this example there is nothing to do
0163
0164   // that's it!
0165   return 0;
0166 }