ShowerLibraryProducer/bin/CastorShowerLibraryMerger.cpp

0001 //
0002 //
0003 //
0004 // Package:        CastorShowerLibraryProducer
0005 // Program:        CastorShowerLibraryMerger
0006 //
0007 // Implementation
0008 //
0009 // Original Author: Maria Elena Pol (polme@mail.cern.ch)
0010 //                  Luiz Mundim     (mundim@mail.cern.ch)
0011 //         Created: 14/Jun/2010
0012 //
0013 /////////////////////////////////////////////////////////////////////
0014 //
0015 #include "FWCore/Utilities/interface/EDMException.h"
0016 #include "FWCore/Utilities/interface/Exception.h"
0017
0018 #include "TROOT.h"
0019 #include "TTree.h"
0020 #include "TFile.h"
0021 #include "TBranchObject.h"
0022 #include "TMath.h"
0023 #include "TLorentzVector.h"
0024 #include "TF1.h"
0025 #include <iostream>
0026 #include <sstream>
0027 #include <iomanip>
0028 #include <cstdlib>
0029 #include <cassert>
0030 #include <unistd.h>
0031 #include <string>
0032 #include <vector>
0033
0034 // Classes for shower library Root file
0035 #include "SimDataFormats/CaloHit/interface/CastorShowerLibraryInfo.h"
0036 #include "SimDataFormats/CaloHit/interface/CastorShowerEvent.h"
0037
0038 void Usage();
0039
0040 int main(int argc, char *argv[]) {
0041   std::vector<std::string> FilesToMerge;
0042   std::string OutPutFile;
0043   if (argc < 3)
0044     Usage();
0045   for (int i = 1; i < argc - 1; i++) {
0046     if (access(argv[i], R_OK) == -1) {
0047       std::cout << "CastorShowerLibraryMerger: File not readable. Existing." << std::endl;
0048       std::cout << "                          " << argv[i] << std::endl;
0049       exit(1);
0050     }
0051     FilesToMerge.push_back(std::string(argv[i]));
0052   }
0053   if (access(argv[argc - 1], R_OK) != -1) {
0054     std::cout << "CastorShowerLibraryMerger: Output File already exists. Exiting." << std::endl;
0055     std::cout << "                           argv[i]" << std::endl;
0056     exit(1);
0057   }
0058   OutPutFile = std::string(argv[argc - 1]);
0059
0060   TFile *of = new TFile(OutPutFile.c_str(), "RECREATE");
0061   // Check that TFile has been successfully opened
0062   if (!of->IsOpen()) {
0063     //edm::LogError("CastorShowerLibraryMerger") << "Opening " << OutPutFile << " failed";
0064     throw cms::Exception("Unknown", "CastorShowerLibraryMerger") << "Opening of " << OutPutFile << " fails\n";
0065   } else {
0066     /*edm::LogInfo*/ std::cout << ("CastorShowerLibraryMerger") << "Successfully openned " << OutPutFile
0067                                << " for output.\n";
0068   }
0069   // Create the output tree
0070   Int_t split = 1;
0071   Int_t bsize = 64000;
0072   TTree *outTree = new TTree("CastorCherenkovPhotons", "Cherenkov Photons");
0073   CastorShowerLibraryInfo *emInfo_out = new CastorShowerLibraryInfo();
0074   CastorShowerEvent *emShower = new CastorShowerEvent();
0075   CastorShowerLibraryInfo *hadInfo_out = new CastorShowerLibraryInfo();
0076   CastorShowerEvent *hadShower = new CastorShowerEvent();
0077   outTree->Branch("emShowerLibInfo.", "CastorShowerLibraryInfo", &emInfo_out, bsize, split);
0078   outTree->Branch("emParticles.", "CastorShowerEvent", &emShower, bsize, split);
0079   outTree->Branch("hadShowerLibInfo.", "CastorShowerLibraryInfo", &hadInfo_out, bsize, split);
0080   outTree->Branch("hadParticles.", "CastorShowerEvent", &hadShower, bsize, split);
0081
0082   // rebuild info branch
0083   CastorShowerLibraryInfo *emInfo = new CastorShowerLibraryInfo();
0084   CastorShowerLibraryInfo *hadInfo = new CastorShowerLibraryInfo();
0085   // Check for the TBranch holding EventInfo in "Events" TTree
0086   std::vector<double> ebin_em;
0087   std::vector<double> etabin_em;
0088   std::vector<double> phibin_em;
0089   std::vector<double> ebin_had;
0090   std::vector<double> etabin_had;
0091   std::vector<double> phibin_had;
0092   int nevt_perbin_e_em = 0;
0093   int nevt_tot_em = 0;
0094   int nevt_perbin_e_had = 0;
0095   int nevt_tot_had = 0;
0096   for (int i = 0; i < int(FilesToMerge.size()); i++) {
0097     TFile in(FilesToMerge.at(i).c_str(), "r");
0098     TTree *event = (TTree *)in.Get("CastorCherenkovPhotons");
0099     TBranchObject *emInfo_b = (TBranchObject *)event->GetBranch("emShowerLibInfo.");
0100     TBranchObject *hadInfo_b = (TBranchObject *)event->GetBranch("hadShowerLibInfo.");
0101     if (emInfo_b) {
0102       emInfo_b->SetAddress(&emInfo);
0103       emInfo_b->GetEntry(0);
0104       ebin_em.insert(ebin_em.end(), emInfo->Energy.getBin().begin(), emInfo->Energy.getBin().end());
0105       if (nevt_perbin_e_em > 0) {
0106         if (nevt_perbin_e_em != int(emInfo->Energy.getNEvtPerBin())) {
0107           std::cout << "CastorShowerLibraryMerger: ERROR: Number of events per energy bin not the same. Exiting."
0108                     << std::endl;
0109           exit(1);
0110         }
0111       } else
0112         nevt_perbin_e_em = emInfo->Energy.getNEvtPerBin();
0113       nevt_tot_em += emInfo->Energy.getNEvts();
0114
0115       if (!emInfo_out->Eta.getBin().empty()) {
0116         if (emInfo_out->Eta.getBin() != emInfo->Eta.getBin()) {
0117           std::cout << "CastorShowerLibraryMerger: ERROR: Eta bins not the same in all files. Exiting." << std::endl;
0118           exit(1);
0119         }
0120       } else {
0121         emInfo_out->Eta.setBin(emInfo->Eta.getBin());
0122         emInfo_out->Eta.setNEvtPerBin(emInfo->Eta.getNEvtPerBin());
0123         emInfo_out->Eta.setNBins(emInfo->Eta.getNBins());
0124         emInfo_out->Eta.setNEvts(emInfo->Eta.getNEvts());
0125       }
0126
0127       if (!emInfo_out->Phi.getBin().empty()) {
0128         if (emInfo_out->Phi.getBin() != emInfo->Phi.getBin()) {
0129           std::cout << "CastorShowerLibraryMerger: ERROR: Phi bins not the same in all files. Exiting." << std::endl;
0130           exit(1);
0131         }
0132       } else {
0133         emInfo_out->Phi.setBin(emInfo->Phi.getBin());
0134         emInfo_out->Phi.setNEvtPerBin(emInfo->Phi.getNEvtPerBin());
0135         emInfo_out->Phi.setNBins(emInfo->Phi.getNBins());
0136         emInfo_out->Phi.setNEvts(emInfo->Phi.getNEvts());
0137       }
0138     }
0139     if (hadInfo_b) {
0140       hadInfo_b->SetAddress(&hadInfo);
0141       hadInfo_b->GetEntry(0);
0142       ebin_had.insert(ebin_had.end(), hadInfo->Energy.getBin().begin(), hadInfo->Energy.getBin().end());
0143       if (nevt_perbin_e_had > 0) {
0144         if (nevt_perbin_e_had != int(hadInfo->Energy.getNEvtPerBin())) {
0145           std::cout << "CastorShowerLibraryMerger: ERROR: Number of events per energy bin not the same. Exiting."
0146                     << std::endl;
0147           exit(1);
0148         }
0149       } else
0150         nevt_perbin_e_had = hadInfo->Energy.getNEvtPerBin();
0151       nevt_tot_had += hadInfo->Energy.getNEvts();
0152
0153       if (!hadInfo_out->Eta.getBin().empty()) {
0154         if (hadInfo_out->Eta.getBin() != hadInfo->Eta.getBin()) {
0155           std::cout << "CastorShowerLibraryMerger: ERROR: Eta bins not the same in all files. Exiting." << std::endl;
0156           exit(1);
0157         }
0158       } else {
0159         hadInfo_out->Eta.setBin(hadInfo->Eta.getBin());
0160         hadInfo_out->Eta.setNEvtPerBin(hadInfo->Eta.getNEvtPerBin());
0161         hadInfo_out->Eta.setNBins(hadInfo->Eta.getNBins());
0162         hadInfo_out->Eta.setNEvts(hadInfo->Eta.getNEvts());
0163       }
0164       if (!hadInfo_out->Phi.getBin().empty()) {
0165         if (hadInfo_out->Phi.getBin() != hadInfo->Phi.getBin()) {
0166           std::cout << "CastorShowerLibraryMerger: ERROR: Phi bins not the same in all files. Exiting." << std::endl;
0167           exit(1);
0168         }
0169       } else {
0170         hadInfo_out->Phi.setBin(hadInfo->Phi.getBin());
0171         hadInfo_out->Phi.setNEvtPerBin(hadInfo->Phi.getNEvtPerBin());
0172         hadInfo_out->Phi.setNBins(hadInfo->Phi.getNBins());
0173         hadInfo_out->Phi.setNEvts(hadInfo->Phi.getNEvts());
0174       }
0175     }
0176     // put ther new info data into emInfo_out and hadInfo_out
0177     in.Close();
0178   }
0179   emInfo_out->Energy.setBin(ebin_em);
0180   emInfo_out->Energy.setNEvtPerBin(nevt_perbin_e_em);
0181   emInfo_out->Energy.setNEvts(nevt_tot_em);
0182   emInfo_out->Energy.setNBins(ebin_em.size());
0183   hadInfo_out->Energy.setBin(ebin_had);
0184   hadInfo_out->Energy.setNEvtPerBin(nevt_perbin_e_had);
0185   hadInfo_out->Energy.setNEvts(nevt_tot_had);
0186   hadInfo_out->Energy.setNBins(ebin_had.size());
0187
0188   // Loop over events from input files merging them into a new file, sequentially
0189   for (int i = 0; i < int(FilesToMerge.size()); i++) {
0190     TFile in(FilesToMerge.at(i).c_str());
0191     TTree *event = (TTree *)in.Get("CastorCherenkovPhotons");
0192     TBranchObject *emShower_b = (TBranchObject *)event->GetBranch("emParticles.");
0193     TBranchObject *hadShower_b = (TBranchObject *)event->GetBranch("hadParticles.");
0194     emShower_b->SetAddress(&emShower);
0195     hadShower_b->SetAddress(&hadShower);
0196     int nevents = event->GetEntries();
0197     for (int n = 0; n < nevents; n++) {
0198       event->GetEntry(n);
0199       outTree->Fill();
0200     }
0201   }
0202   of->cd();
0203   outTree->Write();
0204   of->Close();
0205 }
0206
0207 void Usage() {
0208   std::cout << "Usage: CastorShowerLibraryMerger input_file1 input_file2 [...] output_file" << std::endl;
0209   exit(1);
0210 }