MixingModule/plugins/TestMixedSource.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    TestMixedSource
0004 // Class:      TestMixedSource
0005 //
0006 /**\class TestMixedSource TestMixedSource.cc TestMixed/TestMixedSource/src/TestMixedSource.cc
0007
0008  Description: <one line class summary>
0009
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Emilia Lubenova Becheva
0015 //         Created:  Wed May 20 16:46:58 CEST 2009
0016 //
0017 //
0018
0019 // system include files
0020 #include <memory>
0021
0022 // user include files
0023 #include "FWCore/Framework/interface/Frameworkfwd.h"
0024
0025 #include "FWCore/Framework/interface/Event.h"
0026 #include "FWCore/Framework/interface/MakerMacros.h"
0027
0028 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0029
0030 #include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h"
0031 #include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
0032 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0033
0034 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0035 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0036 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0037 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0038
0039 #include "TH1I.h"
0040 #include "TFile.h"
0041
0042 #include "TestMixedSource.h"
0043
0044 #include <iostream>
0045 #include <fstream>
0046
0047 //
0048 // constructors and destructor
0049 //
0050 namespace edm {
0051   TestMixedSource::TestMixedSource(const edm::ParameterSet& iConfig)
0052       : fileName_(iConfig.getParameter<std::string>("fileName")),
0053         minbunch_(iConfig.getParameter<int>("minBunch")),
0054         maxbunch_(iConfig.getParameter<int>("maxBunch")) {
0055     histTrack_bunchSignal_ =
0056         new TH1I("histoTrackSignal", "Bunchcrossings", maxbunch_ - minbunch_ + 1, minbunch_, maxbunch_ + 1);
0057     histTrack_bunchPileups_ =
0058         new TH1I("histoTrackPileups", "Bunchcrossings", maxbunch_ - minbunch_ + 1, minbunch_, maxbunch_ + 1);
0059
0060     // Vertex
0061     histVertex_bunch_ = new TH1I("histoVertex", "Bunchcrossings", maxbunch_ - minbunch_ + 1, minbunch_, maxbunch_ + 1);
0062
0063     // PCaloHit
0064     histPCaloHit_bunch_ =
0065         new TH1I("histoPCaloHit", "Bunchcrossings", maxbunch_ - minbunch_ + 1, minbunch_, maxbunch_ + 1);
0066
0067     // PSimHit
0068     histPSimHit_bunchSignal_TrackerHitsTECHighTof_ = new TH1I("histoPSimHitTrackerHitsTECHighTofSignal",
0069                                                               "Bunchcrossings",
0070                                                               maxbunch_ - minbunch_ + 1,
0071                                                               minbunch_,
0072                                                               maxbunch_ + 1);
0073     histPSimHit_bunchPileups_TrackerHitsTECHighTof_ = new TH1I("histoPSimHitTrackerHitsTECHighTofPileups",
0074                                                                "Bunchcrossings",
0075                                                                maxbunch_ - minbunch_ + 1,
0076                                                                minbunch_,
0077                                                                maxbunch_ + 1);
0078     histPSimHit_bunchSignal_MuonCSCHits_ = new TH1I(
0079         "histoPSimHitMuonCSCHitsSignal", "Bunchcrossings", maxbunch_ - minbunch_ + 1, minbunch_, maxbunch_ + 1);
0080     histPSimHit_bunchPileups_MuonCSCHits_ = new TH1I(
0081         "histoPSimHitMuonCSCHitsPileups", "Bunchcrossings", maxbunch_ - minbunch_ + 1, minbunch_, maxbunch_ + 1);
0082
0083     int bsp = 25;  //bunchspace
0084     tofhist_ =
0085         new TH1I("TrackerHit_Tof_bcr", "TrackerHit_ToF", 100, float(bsp * minbunch_), float(bsp * maxbunch_) + 50.);
0086     tofhist_sig_ = new TH1I(
0087         "SignalTrackerHit_Tof_bcr", "TrackerHit_ToF", 100, float(bsp * minbunch_), float(bsp * maxbunch_) + 50.);
0088
0089     // HepMCProduct
0090     histHepMCProduct_bunch_ =
0091         new TH1I("histoHepMCProduct", "Bunchcrossings", maxbunch_ - minbunch_ + 1, minbunch_, maxbunch_ + 1);
0092
0093     // Tokens
0094
0095     edm::InputTag tag = edm::InputTag("mix", "g4SimHits");
0096
0097     SimTrackToken_ = consumes<CrossingFrame<SimTrack>>(tag);
0098     SimVertexToken_ = consumes<CrossingFrame<SimVertex>>(tag);
0099
0100     tag = edm::InputTag("mix", "g4SimHitsTrackerHitsTECHighTof");
0101     TrackerToken0_ = consumes<CrossingFrame<PSimHit>>(tag);
0102
0103     tag = edm::InputTag("mix", "g4SimHitsEcalHitsEB");
0104     CaloToken1_ = consumes<CrossingFrame<PCaloHit>>(tag);
0105
0106     tag = edm::InputTag("mix", "g4SimHitsMuonCSCHits");
0107     MuonToken_ = consumes<CrossingFrame<PSimHit>>(tag);
0108
0109     tag = edm::InputTag("mix", "generatorSmeared");
0110     HepMCToken_ = consumes<CrossingFrame<edm::HepMCProduct>>(tag);
0111   }
0112
0113   TestMixedSource::~TestMixedSource() { std::cout << " Destructor TestMixed" << std::endl; }
0114
0115   //
0116   // member functions
0117   //
0118
0119   // ------------ method called to for each event  ------------
0120   void TestMixedSource::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0121     // test SimTracks
0122     //----------------------
0123     edm::Handle<CrossingFrame<SimTrack>> cf_simtrack;
0124     bool gotTracks = iEvent.getByToken(SimTrackToken_, cf_simtrack);
0125     if (!gotTracks)
0126       outputFile << " Could not read SimTracks!!!!"
0127                  << " Please, check if the object SimTracks has been declared in the"
0128                  << " MixingModule configuration file." << std::endl;
0129
0130     if (gotTracks) {
0131       outputFile << "\n=================== Starting SimTrack access ===================" << std::endl;
0132
0133       std::unique_ptr<MixCollection<SimTrack>> col1(new MixCollection<SimTrack>(cf_simtrack.product()));
0134       MixCollection<SimTrack>::iterator cfi1;
0135       int count1 = 0;
0136       std::cout << " \nWe got " << col1->sizeSignal() << " signal tracks and " << col1->sizePileup()
0137                 << " pileup tracks, total: " << col1->size() << std::endl;
0138       for (cfi1 = col1->begin(); cfi1 != col1->end(); cfi1++) {
0139         //std::cout << " BUNCH cfi1.bunch() = " << cfi1.bunch() << std::endl;
0140         if (cfi1.getTrigger() == 0) {
0141           histTrack_bunchPileups_->Fill(cfi1.bunch());
0142         }
0143
0144         if (cfi1.getTrigger() == 1) {
0145           histTrack_bunchSignal_->Fill(cfi1.bunch());
0146         }
0147
0148         int a = count1 % 4;
0149         if (a == 3) {
0150           outputFile << " SimTrack " << count1 << " has genpart index  " << cfi1->genpartIndex() << " vertex Index "
0151                      << cfi1->vertIndex() << " bunchcr " << cfi1.bunch() << " trigger " << cfi1.getTrigger()
0152                      << ", from EncodedEventId: " << cfi1->eventId().bunchCrossing() << " " << cfi1->eventId().event()
0153                      << std::endl;
0154         }
0155         count1++;
0156       }
0157     }
0158
0159     // test SimVertices
0160     //---------------------
0161     edm::Handle<CrossingFrame<SimVertex>> cf_simvtx;
0162     bool gotSimVertex = iEvent.getByToken(SimVertexToken_, cf_simvtx);
0163     if (!gotSimVertex)
0164       outputFile << " Could not read Simvertices !!!!" << std::endl;
0165     else {
0166       outputFile << "\n=================== Starting SimVertex access ===================" << std::endl;
0167       std::unique_ptr<MixCollection<SimVertex>> col2(new MixCollection<SimVertex>(cf_simvtx.product()));
0168       MixCollection<SimVertex>::iterator cfi2;
0169       int count2 = 0;
0170       outputFile << " \nWe got " << col2->sizeSignal() << " signal vertices and " << col2->sizePileup()
0171                  << " pileup vertices, total: " << col2->size() << std::endl;
0172       for (cfi2 = col2->begin(); cfi2 != col2->end(); cfi2++) {
0173         histVertex_bunch_->Fill(cfi2.bunch());
0174         int b = count2 % 4;
0175         if (count2 == 0 || b == 3) {
0176           //outputFile<<" SimVertex "<<count2<<" has parent index  "<<cfi2->parentIndex()<<" bunchcr "<<cfi2.bunch()<<" trigger "<<cfi2.getTrigger()<<", from EncodedEventId: "<<cfi2->eventId().bunchCrossing() <<" "<<cfi2->eventId().event() <<std::endl;
0177         }
0178         //SimVertex myvtx=(*cfi2);
0179         //outputFile<<"Same with op*: "<<count2<<" has parent index  "<<myvtx.parentIndex()<<" bunchcr "<<cfi2.bunch()<<" trigger "<<cfi2.getTrigger()<<", from EncodedEventId: "<<myvtx.eventId().bunchCrossing() <<" "<<myvtx.eventId().event() <<std::endl;
0180         count2++;
0181       }
0182     }
0183
0184     //test HepMCProducts
0185     //------------------------------------
0186     edm::Handle<CrossingFrame<edm::HepMCProduct>> cf_hepmc;
0187     bool gotHepMCP = iEvent.getByToken(HepMCToken_, cf_hepmc);
0188     if (!gotHepMCP)
0189       std::cout << " Could not read HepMCProducts!!!!" << std::endl;
0190     else {
0191       outputFile << "\n=================== Starting HepMCProduct access ===================" << std::endl;
0192       std::unique_ptr<MixCollection<edm::HepMCProduct>> colhepmc(
0193           new MixCollection<edm::HepMCProduct>(cf_hepmc.product()));
0194       MixCollection<edm::HepMCProduct>::iterator cfihepmc;
0195
0196       int count3 = 0;
0197       outputFile << " \nWe got " << colhepmc->sizeSignal() << " signal hepmc products and " << colhepmc->sizePileup()
0198                  << " pileup hepmcs, total: " << colhepmc->size() << std::endl;
0199       for (cfihepmc = colhepmc->begin(); cfihepmc != colhepmc->end(); cfihepmc++) {
0200         histHepMCProduct_bunch_->Fill(cfihepmc.bunch());
0201         int c = count3 % 4;
0202         if (count3 == 0 || c == 3) {
0203           //outputFile<<" edm::HepMCProduct "<<count3<<" has event number "<<cfihepmc->GetEvent()->event_number()<<", "<< cfihepmc->GetEvent()->particles_size()<<" particles and "<<cfihepmc->GetEvent()->vertices_size()<<" vertices,  bunchcr= "<<cfihepmc.bunch()<<" trigger= "<<cfihepmc.getTrigger() <<" sourcetype= "<<cfihepmc.getSourceType()<<std::endl;
0204         }
0205         HepMCProduct myprod = colhepmc->getObject(count3);
0206         //outputFile<<"same with getObject:hepmc product   "<<count3<<" has event number "<<myprod.GetEvent()->event_number()<<", "<<myprod.GetEvent()->particles_size()<<" particles and "<<myprod.GetEvent()->vertices_size()<<" vertices"<<std::endl;
0207         count3++;
0208       }
0209     }
0210
0211     // test CaloHits
0212     //--------------------------
0213     const std::string subdetcalo("g4SimHitsEcalHitsEB");
0214     edm::Handle<CrossingFrame<PCaloHit>> cf_calo;
0215     bool gotPCaloHit = iEvent.getByToken(CaloToken1_, cf_calo);
0216     if (!gotPCaloHit)
0217       outputFile << " Could not read CaloHits with label " << subdetcalo << "!!!!" << std::endl;
0218     else {
0219       outputFile << "\n\n=================== Starting CaloHit access, subdet " << subdetcalo
0220                  << "  ===================" << std::endl;
0221       std::unique_ptr<MixCollection<PCaloHit>> colcalo(new MixCollection<PCaloHit>(cf_calo.product()));
0222       //outputFile<<*(colcalo.get())<<std::endl;
0223       MixCollection<PCaloHit>::iterator cficalo;
0224       int count4 = 0;
0225       for (cficalo = colcalo->begin(); cficalo != colcalo->end(); cficalo++) {
0226         histPCaloHit_bunch_->Fill(cficalo.bunch());
0227         int d = count4 % 4;
0228         if (count4 == 0 || d == 3) {
0229           //outputFile<<" CaloHit "<<count4<<" has tof "<<cficalo->time()<<" trackid "<<cficalo->geantTrackId() <<" bunchcr "<<cficalo.bunch()<<" trigger "<<cficalo.getTrigger()<<", from EncodedEventId: "<<cficalo->eventId().bunchCrossing()<<" " <<cficalo->eventId().event() <<std::endl;
0230         }
0231         count4++;
0232       }
0233     }
0234
0235     // test PSimHit for one particular subdet
0236     //---------------------------------------
0237
0238     const std::string subdet("g4SimHitsTrackerHitsTECHighTof");
0239     edm::Handle<CrossingFrame<PSimHit>> cf_simhit;
0240     bool gotPSimHit = iEvent.getByToken(TrackerToken0_, cf_simhit);
0241     if (!gotPSimHit)
0242       outputFile << " Could not read SimHits with label " << subdet << "!!!!" << std::endl;
0243     else {
0244       outputFile << "\n\n=================== Starting SimHit access, subdet " << subdet
0245                  << "  ===================" << std::endl;
0246
0247       std::unique_ptr<MixCollection<PSimHit>> col(new MixCollection<PSimHit>(cf_simhit.product()));
0248       //outputFile<<*(col.get())<<std::endl;
0249       MixCollection<PSimHit>::iterator cfi;
0250       int count5 = 0;
0251       for (cfi = col->begin(); cfi != col->end(); cfi++) {
0252         // Signal
0253         if (cfi.getTrigger() == 1) {
0254           histPSimHit_bunchSignal_TrackerHitsTECHighTof_->Fill(cfi.bunch());
0255           tofhist_sig_->Fill(cfi->timeOfFlight());
0256         }
0257
0258         // Pileups
0259         if (cfi.getTrigger() == 0) {
0260           histPSimHit_bunchPileups_TrackerHitsTECHighTof_->Fill(cfi.bunch());
0261           std::cout << " cfi->timeOfFlight() = " << cfi->timeOfFlight() << std::endl;
0262           tofhist_->Fill(cfi->timeOfFlight());
0263         }
0264
0265         int e = count5 % 4;
0266         if (e == 3) {
0267           outputFile << " Hit " << count5 << " has tof " << cfi->timeOfFlight() << " trackid " << cfi->trackId()
0268                      << " bunchcr " << cfi.bunch() << " trigger " << cfi.getTrigger()
0269                      << ", from EncodedEventId: " << cfi->eventId().bunchCrossing() << " " << cfi->eventId().event()
0270                      << " bcr from MixCol " << cfi.bunch() << std::endl;
0271         }
0272         count5++;
0273       }
0274     }
0275
0276     const std::string subdet1("g4SimHitsMuonCSCHits");
0277     edm::Handle<CrossingFrame<PSimHit>> cf_simhit1;
0278     bool gotPSimHit1 = iEvent.getByToken(MuonToken_, cf_simhit1);
0279     if (!gotPSimHit1)
0280       outputFile << " Could not read SimHits with label " << subdet1 << "!!!!" << std::endl;
0281     else {
0282       outputFile << "\n\n=================== Starting SimHit access, subdet " << subdet1
0283                  << "  ===================" << std::endl;
0284
0285       std::unique_ptr<MixCollection<PSimHit>> col(new MixCollection<PSimHit>(cf_simhit1.product()));
0286       //outputFile<<*(col.get())<<std::endl;
0287       MixCollection<PSimHit>::iterator cfi;
0288       int count5 = 0;
0289       for (cfi = col->begin(); cfi != col->end(); cfi++) {
0290         if (cfi.getTrigger() == 1)
0291           histPSimHit_bunchSignal_MuonCSCHits_->Fill(cfi.bunch());
0292
0293         if (cfi.getTrigger() == 0)
0294           histPSimHit_bunchPileups_MuonCSCHits_->Fill(cfi.bunch());
0295
0296         int e = count5 % 4;
0297         if (e == 3) {
0298           outputFile << " Hit " << count5 << " has tof " << cfi->timeOfFlight() << " trackid " << cfi->trackId()
0299                      << " bunchcr " << cfi.bunch() << " trigger " << cfi.getTrigger()
0300                      << ", from EncodedEventId: " << cfi->eventId().bunchCrossing() << " " << cfi->eventId().event()
0301                      << " bcr from MixCol " << cfi.bunch() << std::endl;
0302         }
0303         count5++;
0304       }
0305     }
0306   }
0307
0308   // ------------ method called once each job just before starting event loop  ------------
0309   void TestMixedSource::beginJob() {
0310     outputFile.open("test.log");
0311     if (!outputFile.is_open()) {
0312       std::cout << "Unable to open file!" << std::endl;
0313     }
0314   }
0315
0316   // ------------ method called once each job just after ending the event loop  ------------
0317   void TestMixedSource::endJob() {
0318     outputFile.close();
0319
0320     char t[30];
0321     sprintf(t, "%s", fileName_.c_str());
0322     std::cout << " fileName = " << t << std::endl;
0323     histFile_ = new TFile(t, "RECREATE");
0324     histTrack_bunchPileups_->Write();
0325     histTrack_bunchSignal_->Write();
0326     histVertex_bunch_->Write();
0327     histPCaloHit_bunch_->Write();
0328     histPSimHit_bunchPileups_TrackerHitsTECHighTof_->Write();
0329     histPSimHit_bunchSignal_TrackerHitsTECHighTof_->Write();
0330     tofhist_sig_->Write();
0331     tofhist_->Write();
0332     histPSimHit_bunchSignal_MuonCSCHits_->Write();
0333     histPSimHit_bunchPileups_MuonCSCHits_->Write();
0334     histHepMCProduct_bunch_->Write();
0335     histFile_->Write();
0336     histFile_->Close();
0337   }
0338 }  // namespace edm