File indexing completed on 2023-03-17 11:24:39
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
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
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
0062 if (!of->IsOpen()) {
0063
0064 throw cms::Exception("Unknown", "CastorShowerLibraryMerger") << "Opening of " << OutPutFile << " fails\n";
0065 } else {
0066 std::cout << ("CastorShowerLibraryMerger") << "Successfully openned " << OutPutFile
0067 << " for output.\n";
0068 }
0069
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
0083 CastorShowerLibraryInfo *emInfo = new CastorShowerLibraryInfo();
0084 CastorShowerLibraryInfo *hadInfo = new CastorShowerLibraryInfo();
0085
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
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
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 }