File indexing completed on 2024-04-06 11:57:26
0001 #include <string>
0002 #include <fstream>
0003 #include <map>
0004 #include <iostream>
0005
0006 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0007 #include "FWCore/Framework/interface/EventPrincipal.h"
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/Framework/interface/EventSetup.h"
0010 #include "FWCore/Framework/interface/ESHandle.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0014 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0015 #include "FWCore/Utilities/interface/InputTag.h"
0016
0017 #include "TFile.h"
0018 #include "TString.h"
0019 #include "TChain.h"
0020 #include "TStopwatch.h"
0021
0022 class TkAlCaSkimTreeMerger : public edm::one::EDAnalyzer<> {
0023 public:
0024 TkAlCaSkimTreeMerger(const edm::ParameterSet &iConfig);
0025 ~TkAlCaSkimTreeMerger() override;
0026 void beginJob() override;
0027 void endJob() override;
0028 void analyze(const edm::Event &, const edm::EventSetup &) override;
0029 static void fillDescriptions(edm::ConfigurationDescriptions &);
0030
0031 private:
0032 TTree *out_;
0033 TTree *firsttree_;
0034 TChain *ch_;
0035 std::string filelist_;
0036 std::string firstfilename_;
0037 std::string treename_;
0038 std::string outfilename_;
0039
0040
0041 typedef std::map<uint32_t, uint32_t> DetHitMap;
0042 DetHitMap hitmap_;
0043 DetHitMap overlapmap_;
0044 int maxhits_;
0045 edm::ParameterSet maxhitsSet_;
0046 int maxPXBhits_, maxPXFhits_, maxTIBhits_, maxTIDhits_, maxTOBhits_, maxTEChits_;
0047
0048 TStopwatch myclock;
0049 };
0050
0051 TkAlCaSkimTreeMerger::TkAlCaSkimTreeMerger(const edm::ParameterSet &iConfig)
0052 : filelist_(iConfig.getParameter<std::string>("FileList")),
0053 treename_(iConfig.getParameter<std::string>("TreeName")),
0054 outfilename_(iConfig.getParameter<std::string>("OutputFile")),
0055 maxhits_(iConfig.getParameter<int32_t>("NhitsMaxLimit")),
0056 maxhitsSet_(iConfig.getParameter<edm::ParameterSet>("NhitsMaxSet")) {
0057 maxPXBhits_ = maxhitsSet_.getParameter<int32_t>("PXBmaxhits");
0058 maxPXFhits_ = maxhitsSet_.getParameter<int32_t>("PXFmaxhits");
0059 maxTIBhits_ = maxhitsSet_.getParameter<int32_t>("TIBmaxhits");
0060 maxTIDhits_ = maxhitsSet_.getParameter<int32_t>("TIDmaxhits");
0061 maxTOBhits_ = maxhitsSet_.getParameter<int32_t>("TOBmaxhits");
0062 maxTEChits_ = maxhitsSet_.getParameter<int32_t>("TECmaxhits");
0063
0064 LogDebug("TkAlCaSkimTreeMerger") << "\n\n*** MAX N HITS = " << maxhits_ << std::endl << std::endl;
0065 out_ = nullptr;
0066 firsttree_ = nullptr;
0067 ch_ = nullptr;
0068 }
0069
0070 TkAlCaSkimTreeMerger::~TkAlCaSkimTreeMerger() {
0071 delete ch_;
0072 LogDebug("TkAlCaSkimTreeMerger") << "finished." << std::endl;
0073 }
0074
0075
0076 void TkAlCaSkimTreeMerger::beginJob() {
0077 myclock.Start();
0078
0079
0080 ch_ = new TChain(treename_.c_str());
0081 LogDebug("TkAlCaSkimTreeMerger") << "The chain contains " << ch_->GetNtrees() << " trees" << std::endl;
0082
0083
0084 std::ifstream flist(filelist_.c_str(), std::ios::in);
0085 std::string filename;
0086 std::string firstfilename;
0087 bool first = true;
0088 while (!flist.eof()) {
0089 filename = "";
0090 flist >> filename;
0091 if (filename.empty())
0092 continue;
0093 LogDebug("TkAlCaSkimTreeMerger") << "Adding " << filename << std::endl;
0094 ch_->Add(filename.c_str());
0095 if (first) {
0096 firstfilename_ = filename;
0097 first = false;
0098 }
0099 }
0100 LogDebug("TkAlCaSkimTreeMerger") << "Now the chain contains " << ch_->GetNtrees() << " trees (" << ch_->GetEntries()
0101 << " entries)" << std::endl;
0102
0103 unsigned int id_ch = 0;
0104 uint32_t nhits_ch = 0, noverlaps_ch = 0;
0105 ch_->SetBranchAddress("DetId", &id_ch);
0106 ch_->SetBranchAddress("Nhits", &nhits_ch);
0107 ch_->SetBranchAddress("Noverlaps", &noverlaps_ch);
0108
0109 ch_->SetBranchStatus("SubDet", false);
0110 ch_->SetBranchStatus("Layer", false);
0111 ch_->SetBranchStatus("is2D", false);
0112 ch_->SetBranchStatus("isStereo", false);
0113 ch_->SetBranchStatus("posX", false);
0114 ch_->SetBranchStatus("posY", false);
0115 ch_->SetBranchStatus("posZ", false);
0116 ch_->SetBranchStatus("posR", false);
0117 ch_->SetBranchStatus("posEta", false);
0118 ch_->SetBranchStatus("posPhi", false);
0119
0120 int totnhits(0), totnoverlaps(0);
0121
0122
0123 DetHitMap::iterator mapiter;
0124 DetHitMap::iterator overlapiter;
0125
0126 for (int ent = 0; ent < ch_->GetEntries(); ++ent) {
0127
0128 ch_->GetEntry(ent);
0129 totnhits += nhits_ch;
0130 totnoverlaps += noverlaps_ch;
0131
0132 mapiter = hitmap_.find(id_ch);
0133 if (mapiter != hitmap_.end()) {
0134 hitmap_[id_ch] = hitmap_[id_ch] + nhits_ch;
0135 } else {
0136 hitmap_.insert(std::pair<uint32_t, uint32_t>(id_ch, nhits_ch));
0137 }
0138
0139 overlapiter = overlapmap_.find(id_ch);
0140 if (overlapiter != overlapmap_.end()) {
0141 overlapmap_[id_ch] = overlapmap_[id_ch] + noverlaps_ch;
0142 } else {
0143 overlapmap_.insert(std::pair<uint32_t, uint32_t>(id_ch, noverlaps_ch));
0144 }
0145
0146 }
0147
0148 LogDebug("TkAlCaSkimTreeMerger") << "Nhits in the chain: " << totnhits << std::endl;
0149 LogDebug("TkAlCaSkimTreeMerger") << "NOverlaps in the chain: " << totnoverlaps << std::endl;
0150
0151 myclock.Stop();
0152 LogDebug("TkAlCaSkimTreeMerger") << "Finished beginJob after " << myclock.RealTime() << " s (real time) / "
0153 << myclock.CpuTime() << " s (cpu time)" << std::endl;
0154 myclock.Continue();
0155 }
0156
0157
0158 void TkAlCaSkimTreeMerger::analyze(const edm::Event &, const edm::EventSetup &) {
0159 LogDebug("TkAlCaSkimTreeMerger") << firsttree_->GetEntries() << std::endl;
0160 }
0161
0162
0163 void TkAlCaSkimTreeMerger::endJob() {
0164
0165 TFile *firstfile = new TFile(firstfilename_.c_str(), "READ");
0166 firsttree_ = (TTree *)firstfile->Get(treename_.c_str());
0167 LogDebug("TkAlCaSkimTreeMerger") << "the first tree has " << firsttree_->GetEntries() << " entries" << std::endl;
0168 unsigned int id = 0;
0169 uint32_t nhits = 0, noverlaps = 0;
0170 float posX(-99999.0), posY(-77777.0), posZ(-88888.0);
0171 float posEta(-6666.0), posPhi(-5555.0), posR(-4444.0);
0172 int subdet = 0;
0173 unsigned int layer = 0;
0174
0175 firsttree_->SetBranchAddress("DetId", &id);
0176 firsttree_->SetBranchAddress("Nhits", &nhits);
0177 firsttree_->SetBranchAddress("Noverlaps", &noverlaps);
0178 firsttree_->SetBranchAddress("SubDet", &subdet);
0179 firsttree_->SetBranchAddress("Layer", &layer);
0180
0181
0182 firsttree_->SetBranchAddress("posX", &posX);
0183 firsttree_->SetBranchAddress("posY", &posY);
0184 firsttree_->SetBranchAddress("posZ", &posZ);
0185 firsttree_->SetBranchAddress("posR", &posR);
0186 firsttree_->SetBranchAddress("posEta", &posEta);
0187 firsttree_->SetBranchAddress("posPhi", &posPhi);
0188
0189
0190
0191 TFile *outfile = new TFile(outfilename_.c_str(), "RECREATE");
0192 out_ = new TTree(treename_.c_str(), "AlignmentHitMapsTOTAL");
0193 unsigned int id_out = 0;
0194 uint32_t nhits_out = 0, noverlaps_out = 0;
0195 float posX_out(-99999.0), posY_out(-77777.0), posZ_out(-88888.0);
0196 float posEta_out(-6666.0), posPhi_out(-5555.0), posR_out(-4444.0);
0197 int subdet_out = 0;
0198 unsigned int layer_out = 0;
0199 bool is2D_out = false, isStereo_out = false;
0200 float prescfact_out = 1.0;
0201 float prescfact_overlap_out = 1.0;
0202
0203 out_->Branch("DetId", &id_out, "DetId/i");
0204 out_->Branch("Nhits", &nhits_out, "Nhits/i");
0205 out_->Branch("Noverlaps", &noverlaps_out, "Noverlaps/i");
0206 out_->Branch("SubDet", &subdet_out, "SubDet/I");
0207 out_->Branch("Layer", &layer_out, "Layer/i");
0208 out_->Branch("is2D", &is2D_out, "is2D/B");
0209 out_->Branch("isStereo", &isStereo_out, "isStereo/B");
0210 out_->Branch("posX", &posX_out, "posX/F");
0211 out_->Branch("posY", &posY_out, "posY/F");
0212 out_->Branch("posZ", &posZ_out, "posZ/F");
0213 out_->Branch("posR", &posR_out, "posR/F");
0214 out_->Branch("posEta", &posEta_out, "posEta/F");
0215 out_->Branch("posPhi", &posPhi_out, "posPhi/F");
0216 out_->Branch("PrescaleFactor", &prescfact_out, "PrescaleFact/F");
0217 out_->Branch("PrescaleFactorOverlap", &prescfact_overlap_out, "PrescaleFactOverlap/F");
0218
0219
0220 DetHitMap::iterator mapiter;
0221
0222 for (int mod = 0; mod < firsttree_->GetEntries(); mod++) {
0223
0224
0225
0226
0227 firsttree_->GetEntry(mod);
0228 nhits_out = hitmap_[id];
0229 noverlaps_out = overlapmap_[id];
0230
0231 id_out = id;
0232 subdet_out = subdet;
0233 layer_out = layer;
0234 posX_out = posX;
0235 posY_out = posY;
0236 posZ_out = posZ;
0237 posR_out = posR;
0238 posEta_out = posEta;
0239 posPhi_out = posPhi;
0240
0241
0242 int subdetmax = -1;
0243 if (subdet_out == 1)
0244 subdetmax = maxPXBhits_;
0245 else if (subdet_out == 2)
0246 subdetmax = maxPXFhits_;
0247 else if (subdet_out == 3)
0248 subdetmax = maxTIBhits_;
0249 else if (subdet_out == 4)
0250 subdetmax = maxTIDhits_;
0251 else if (subdet_out == 5)
0252 subdetmax = maxTOBhits_;
0253 else if (subdet_out == 6)
0254 subdetmax = maxTEChits_;
0255 else
0256 subdetmax = -9999;
0257
0258 if (maxhits_ > -1) {
0259 if (int(nhits_out) > maxhits_) {
0260 prescfact_out = float(maxhits_) / float(nhits_out);
0261 }
0262 if (int(noverlaps_out) > maxhits_) {
0263 prescfact_overlap_out = float(maxhits_) / float(noverlaps_out);
0264 }
0265 } else if (subdetmax > 0) {
0266 if (int(nhits_out) > subdetmax) {
0267 prescfact_out = float(subdetmax / nhits_out);
0268 }
0269 if (int(noverlaps_out) > subdetmax) {
0270 prescfact_overlap_out = float(subdetmax) / float(noverlaps_out);
0271 }
0272 } else {
0273 prescfact_out = 1.0;
0274 prescfact_overlap_out = 1.0;
0275 }
0276 out_->Fill();
0277
0278 }
0279
0280 myclock.Stop();
0281 LogDebug("TkAlCaSkimTreeMerger") << "Finished endJob after " << myclock.RealTime() << " s (real time) / "
0282 << myclock.CpuTime() << " s (cpu time)" << std::endl;
0283 LogDebug("TkAlCaSkimTreeMerger") << "Ending the tree merging." << std::endl;
0284 out_->Write();
0285 LogDebug("TkAlCaSkimTreeMerger") << "Deleting..." << std::flush;
0286 delete firstfile;
0287 delete outfile;
0288 }
0289
0290 void TkAlCaSkimTreeMerger::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0291 edm::ParameterSetDescription desc;
0292 desc.setComment("Merger of TkAlCaSkim Trees");
0293 desc.add<std::string>("FileList", "DQMHitMapsList.txt");
0294 desc.add<std::string>("TreeName", "AlignmentHitMaps");
0295 desc.add<std::string>("OutputFile", "AlignmentHitMapsMerged.root");
0296 desc.add<int>("NhitsMaxLimit", 0);
0297
0298 edm::ParameterSetDescription tkAlCaSkimTreeParamsDesc;
0299 std::vector<edm::ParameterSet> tkAlCaSkimDefaults(1);
0300 tkAlCaSkimDefaults[0].addParameter("PXBmaxhits", -1);
0301 tkAlCaSkimDefaults[0].addParameter("PXFmaxhits", -1);
0302 tkAlCaSkimDefaults[0].addParameter("TIBmaxhits", -1);
0303 tkAlCaSkimDefaults[0].addParameter("TIDmaxhits", -1);
0304 tkAlCaSkimDefaults[0].addParameter("TOBmaxhits", -1);
0305 tkAlCaSkimDefaults[0].addParameter("TECmaxhits", -1);
0306 desc.addVPSet("NhitsMaxSet", tkAlCaSkimTreeParamsDesc, tkAlCaSkimDefaults);
0307 descriptions.addWithDefaultLabel(desc);
0308 }
0309
0310
0311 #include "FWCore/PluginManager/interface/ModuleDef.h"
0312 #include "FWCore/Framework/interface/MakerMacros.h"
0313 DEFINE_FWK_MODULE(TkAlCaSkimTreeMerger);