TrackerAlignment/plugins/TkAlCaSkimTreeMerger.cc

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_;            //TTree containing the merged result
0033   TTree *firsttree_;      //first tree of the list; this gives the structure to all the others
0034   TChain *ch_;            //chain containing all the tree you want to merge
0035   std::string filelist_;  //text file containing the list of input files
0036   std::string firstfilename_;
0037   std::string treename_;     //name of the tree you want to merge (contained in the file)
0038   std::string outfilename_;  //name of the file where you want to save the output
0039
0040   //Hit Population
0041   typedef std::map<uint32_t, uint32_t> DetHitMap;
0042   DetHitMap hitmap_;
0043   DetHitMap overlapmap_;
0044   int maxhits_;  //above this number, the hit population is prescaled. Configurable for each subdet
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   //anything you want to do for initializing
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 // ------------ method called before analyzing the first event  ------------
0076 void TkAlCaSkimTreeMerger::beginJob() {
0077   myclock.Start();
0078
0079   //prepare the chain
0080   ch_ = new TChain(treename_.c_str());
0081   LogDebug("TkAlCaSkimTreeMerger") << "The chain contains " << ch_->GetNtrees() << " trees" << std::endl;
0082
0083   //load the trees into the chain
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);  //now only id, nhits and noverlaps are on...
0119
0120   int totnhits(0), totnoverlaps(0);
0121
0122   //look if you find this detid in the map
0123   DetHitMap::iterator mapiter;
0124   DetHitMap::iterator overlapiter;
0125
0126   for (int ent = 0; ent < ch_->GetEntries(); ++ent) {
0127     //  for(int ent=0;ent<100;++ent){
0128     ch_->GetEntry(ent);
0129     totnhits += nhits_ch;
0130     totnoverlaps += noverlaps_ch;
0131
0132     mapiter = hitmap_.find(id_ch);
0133     if (mapiter != hitmap_.end()) {  //present, increase its value
0134       hitmap_[id_ch] = hitmap_[id_ch] + nhits_ch;
0135     } else {  //not present, let's add this key to the map with value=1
0136       hitmap_.insert(std::pair<uint32_t, uint32_t>(id_ch, nhits_ch));
0137     }
0138     //do the same for overlaps
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   }  //end loop on ent - entries in the chain
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 }  //end beginJob
0156
0157 // ------------ method called to for each event  ------------
0158 void TkAlCaSkimTreeMerger::analyze(const edm::Event &, const edm::EventSetup &) {
0159   LogDebug("TkAlCaSkimTreeMerger") << firsttree_->GetEntries() << std::endl;
0160 }  //end analyze
0161
0162 // ------------ method called after having analyzed all the events  ------------
0163 void TkAlCaSkimTreeMerger::endJob() {
0164   //address variables in the first tree and in the chain
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   // bool is2D=false,isStereo=false;
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   //  firsttree_->SetBranchAddress("is2D" ,    &is2D);
0181   // firsttree_->SetBranchAddress("isStereo", &isStereo);
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   //create and book the output
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   //look if you find this detid in the map
0220   DetHitMap::iterator mapiter;
0221
0222   for (int mod = 0; mod < firsttree_->GetEntries(); mod++) {
0223     //for(int mod=0;mod<100;++mod){
0224     //   nhits_out=0;
0225     // noverlaps_out=0;
0226
0227     firsttree_->GetEntry(mod);
0228     nhits_out = hitmap_[id];
0229     noverlaps_out = overlapmap_[id];
0230     // if(mod<25)LogDebug("TkAlCaSkimTreeMerger")<<"Nhits 1st tree: "<<nhits<<"\tTotal nhits chain: "<<nhits_out<<std::endl;
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     //calculate prescaling factor
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) {  //calculate different prescaling factors for each subdet
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   }  //end loop on mod - first tree modules
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 }  //end endJob
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 // ========= MODULE DEF ==============
0311 #include "FWCore/PluginManager/interface/ModuleDef.h"
0312 #include "FWCore/Framework/interface/MakerMacros.h"
0313 DEFINE_FWK_MODULE(TkAlCaSkimTreeMerger);