TestMixedSource.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    TestMixedSource
// Class:      TestMixedSource
// 
//
// Original Author:  Emilia Lubenova Becheva
//         Created:  Wed May 20 16:46:58 CEST 2009
// $Id: TestMixedSource.cc,v 1.5 2011/11/15 21:57:47 gowdy Exp $
//
//


// system include files
#include <memory>

// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EDAnalyzer.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h"
#include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
#include "SimDataFormats/TrackingHit/interface/PSimHit.h"

#include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
#include "SimDataFormats/Track/interface/SimTrackContainer.h"
#include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"

#include "TH1I.h"
#include "TFile.h"

#include "TestMixedSource.h"

#include <iostream>
#include <fstream>


//
// constructors and destructor
//
namespace edm
{
TestMixedSource::TestMixedSource(const edm::ParameterSet& iConfig)
: fileName_(iConfig.getParameter<std::string>("fileName")), minbunch_(iConfig.getParameter<int>("minBunch")),maxbunch_(iConfig.getParameter<int>("maxBunch"))
{ 
  
  histTrack_bunchSignal_ = new TH1I("histoTrackSignal","Bunchcrossings",maxbunch_-minbunch_+1,minbunch_,maxbunch_+1);
  histTrack_bunchPileups_ = new TH1I("histoTrackPileups","Bunchcrossings",maxbunch_-minbunch_+1,minbunch_,maxbunch_+1);

  // Vertex
  histVertex_bunch_ = new TH1I("histoVertex","Bunchcrossings",maxbunch_-minbunch_+1,minbunch_,maxbunch_+1);
  
  // PCaloHit
  histPCaloHit_bunch_ = new TH1I("histoPCaloHit","Bunchcrossings",maxbunch_-minbunch_+1,minbunch_,maxbunch_+1);
  
  // PSimHit
  histPSimHit_bunchSignal_TrackerHitsTECHighTof_ = new TH1I("histoPSimHitTrackerHitsTECHighTofSignal","Bunchcrossings",maxbunch_-minbunch_+1,minbunch_,maxbunch_+1);
  histPSimHit_bunchPileups_TrackerHitsTECHighTof_ = new TH1I("histoPSimHitTrackerHitsTECHighTofPileups","Bunchcrossings",maxbunch_-minbunch_+1,minbunch_,maxbunch_+1);
  histPSimHit_bunchSignal_MuonCSCHits_ = new TH1I("histoPSimHitMuonCSCHitsSignal","Bunchcrossings",maxbunch_-minbunch_+1,minbunch_,maxbunch_+1);
  histPSimHit_bunchPileups_MuonCSCHits_ = new TH1I("histoPSimHitMuonCSCHitsPileups","Bunchcrossings",maxbunch_-minbunch_+1,minbunch_,maxbunch_+1);

  int bsp = 25;//bunchspace
  tofhist_ = new TH1I ("TrackerHit_Tof_bcr","TrackerHit_ToF",100,float(bsp*minbunch_),float(bsp*maxbunch_)+50.);
  tofhist_sig_ = new TH1I("SignalTrackerHit_Tof_bcr","TrackerHit_ToF",100,float(bsp*minbunch_),float(bsp*maxbunch_)+50.);


  // HepMCProduct
  histHepMCProduct_bunch_ = new TH1I("histoHepMCProduct","Bunchcrossings",maxbunch_-minbunch_+1,minbunch_,maxbunch_+1);

}


TestMixedSource::~TestMixedSource()
{ std::cout << " Destructor TestMixed"  << std::endl;

}


//
// member functions
//

// ------------ method called to for each event  ------------
void
TestMixedSource::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{ 
  // test SimTracks
  //----------------------
  edm::Handle<CrossingFrame<SimTrack> > cf_simtrack;
  bool gotTracks = iEvent.getByLabel("mix","g4SimHits",cf_simtrack);
  if (!gotTracks)  outputFile<<" Could not read SimTracks!!!!" 
           << " Please, check if the object SimTracks has been declared in the"
           << " MixingModule configuration file."<<std::endl;
  
  if (gotTracks) {
    outputFile<<"\n=================== Starting SimTrack access ==================="<<std::endl;

    std::auto_ptr<MixCollection<SimTrack> > col1(new MixCollection<SimTrack>(cf_simtrack.product()));
    MixCollection<SimTrack>::iterator cfi1;
    int count1=0;
    std::cout <<" \nWe got "<<col1->sizeSignal()<<" signal tracks and "<<col1->sizePileup()<<" pileup tracks, total: "<<col1->size()<<std::endl;
    for (cfi1=col1->begin(); cfi1!=col1->end();cfi1++) {
      //std::cout << " BUNCH cfi1.bunch() = " << cfi1.bunch() << std::endl;
      if (cfi1.getTrigger()==0){
      histTrack_bunchPileups_->Fill(cfi1.bunch());
      }
      
      if (cfi1.getTrigger()==1){
      histTrack_bunchSignal_->Fill(cfi1.bunch());
      }

      
      int a = count1%4;
      if (a==3){
      outputFile<<" SimTrack "<<count1<<" has genpart index  "<<cfi1->genpartIndex()<<" vertex Index "<<cfi1->vertIndex() <<" bunchcr "<<cfi1.bunch()<<" trigger "<<cfi1.getTrigger()<<", from EncodedEventId: "<<cfi1->eventId().bunchCrossing() <<" "<<cfi1->eventId().event() <<std::endl;
      }
      count1++; 
    }
  }

  
  // test SimVertices
  //---------------------
  edm::Handle<CrossingFrame<SimVertex> > cf_simvtx;
  bool gotSimVertex = iEvent.getByLabel("mix","g4SimHits",cf_simvtx);
  if (!gotSimVertex) outputFile<<" Could not read Simvertices !!!!"<<std::endl;
  else {
    outputFile<<"\n=================== Starting SimVertex access ==================="<<std::endl;
    std::auto_ptr<MixCollection<SimVertex> > col2(new MixCollection<SimVertex>(cf_simvtx.product()));
    MixCollection<SimVertex>::iterator cfi2;
    int count2=0;
    outputFile <<" \nWe got "<<col2->sizeSignal()<<" signal vertices and "<<col2->sizePileup()<<" pileup vertices, total: "<<col2->size()<<std::endl;
    for (cfi2=col2->begin(); cfi2!=col2->end();cfi2++) {
      histVertex_bunch_->Fill(cfi2.bunch());
      int b = count2%4;
      if (count2 == 0 || b==3){
      //outputFile<<" SimVertex "<<count2<<" has parent index  "<<cfi2->parentIndex()<<" bunchcr "<<cfi2.bunch()<<" trigger "<<cfi2.getTrigger()<<", from EncodedEventId: "<<cfi2->eventId().bunchCrossing() <<" "<<cfi2->eventId().event() <<std::endl;
      }
      //SimVertex myvtx=(*cfi2);
      //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;
      count2++; 
    }
  }
  
  
  //test HepMCProducts
  //------------------------------------
  edm::Handle<CrossingFrame<edm::HepMCProduct> > cf_hepmc;
  bool gotHepMCP = iEvent.getByLabel("mix","generator",cf_hepmc);
  if (!gotHepMCP) std::cout<<" Could not read HepMCProducts!!!!"<<std::endl;
  else {
    outputFile<<"\n=================== Starting HepMCProduct access ==================="<<std::endl;
    std::auto_ptr<MixCollection<edm::HepMCProduct> > colhepmc(new MixCollection<edm::HepMCProduct>(cf_hepmc.product()));
    MixCollection<edm::HepMCProduct>::iterator cfihepmc;
    
    int count3=0;
    outputFile <<" \nWe got "<<colhepmc->sizeSignal()<<" signal hepmc products and "<<colhepmc->sizePileup()<<" pileup hepmcs, total: "<<colhepmc->size()<<std::endl;
    for (cfihepmc=colhepmc->begin(); cfihepmc!=colhepmc->end();cfihepmc++) {
      histHepMCProduct_bunch_->Fill(cfihepmc.bunch());
      int c = count3%4;
      if (count3==0 || c==3){
      //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;
      }
      HepMCProduct myprod=colhepmc->getObject(count3);
      //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;
      count3++;
    }
  }


  // test CaloHits
  //--------------------------  
  const std::string subdetcalo("g4SimHitsEcalHitsEB");
  edm::Handle<CrossingFrame<PCaloHit> > cf_calo;
  bool gotPCaloHit = iEvent.getByLabel("mix",subdetcalo,cf_calo);
  if (!gotPCaloHit) outputFile<<" Could not read CaloHits with label "<<subdetcalo<<"!!!!"<<std::endl;
  else {
    outputFile<<"\n\n=================== Starting CaloHit access, subdet "<<subdetcalo<<"  ==================="<<std::endl;
    std::auto_ptr<MixCollection<PCaloHit> > colcalo(new MixCollection<PCaloHit>(cf_calo.product()));
    //outputFile<<*(colcalo.get())<<std::endl;
    MixCollection<PCaloHit>::iterator cficalo;
    int count4=0;
    for (cficalo=colcalo->begin(); cficalo!=colcalo->end();cficalo++) {
      histPCaloHit_bunch_->Fill(cficalo.bunch());
      int d = count4%4;
      if (count4==0 || d==3){
      //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;
      }
      count4++;
    }
  }

  
  // test PSimHit for one particular subdet
  //---------------------------------------
  
  
  const std::string subdet("g4SimHitsTrackerHitsTECHighTof");
  edm::Handle<CrossingFrame<PSimHit> > cf_simhit;
  bool gotPSimHit = iEvent.getByLabel("mix",subdet,cf_simhit);
  if (!gotPSimHit) outputFile<<" Could not read SimHits with label "<<subdet<<"!!!!"<<std::endl;
  else {
    outputFile<<"\n\n=================== Starting SimHit access, subdet "<<subdet<<"  ==================="<<std::endl;

    std::auto_ptr<MixCollection<PSimHit> > col(new MixCollection<PSimHit>(cf_simhit.product()));
    //outputFile<<*(col.get())<<std::endl;
    MixCollection<PSimHit>::iterator cfi;
    int count5=0;
    for (cfi=col->begin(); cfi!=col->end();cfi++) {
      
      // Signal
      if (cfi.getTrigger()==1){ 
        histPSimHit_bunchSignal_TrackerHitsTECHighTof_->Fill(cfi.bunch());
        tofhist_sig_->Fill(cfi->timeOfFlight());
      }
      
      // Pileups
      if (cfi.getTrigger()==0){
       histPSimHit_bunchPileups_TrackerHitsTECHighTof_->Fill(cfi.bunch());
       std::cout << " cfi->timeOfFlight() = " << cfi->timeOfFlight() << std::endl;
       tofhist_->Fill(cfi->timeOfFlight()); 
      }
      
      int e = count5%4;
      if (e==3){
      outputFile<<" Hit "<<count5<<" has tof "<<cfi->timeOfFlight()<<" trackid "<<cfi->trackId() <<" bunchcr "<<cfi.bunch()<<" trigger "<<cfi.getTrigger()<<", from EncodedEventId: "<<cfi->eventId().bunchCrossing()<<" " <<cfi->eventId().event() <<" bcr from MixCol "<<cfi.bunch()<<std::endl;
      }
      count5++;
    }
  }
  
  const std::string subdet1("g4SimHitsMuonCSCHits");
  edm::Handle<CrossingFrame<PSimHit> > cf_simhit1;
  bool gotPSimHit1 = iEvent.getByLabel("mix",subdet1,cf_simhit1);
  if (!gotPSimHit1) outputFile<<" Could not read SimHits with label "<<subdet1<<"!!!!"<<std::endl;
  else {
    outputFile<<"\n\n=================== Starting SimHit access, subdet "<<subdet1<<"  ==================="<<std::endl;

    std::auto_ptr<MixCollection<PSimHit> > col(new MixCollection<PSimHit>(cf_simhit1.product()));
    //outputFile<<*(col.get())<<std::endl;
    MixCollection<PSimHit>::iterator cfi;
    int count5=0;
    for (cfi=col->begin(); cfi!=col->end();cfi++) {
      
      if (cfi.getTrigger()==1) histPSimHit_bunchSignal_MuonCSCHits_->Fill(cfi.bunch());
      
      if (cfi.getTrigger()==0) histPSimHit_bunchPileups_MuonCSCHits_->Fill(cfi.bunch());
      
      int e = count5%4;
      if (e==3){
      outputFile<<" Hit "<<count5<<" has tof "<<cfi->timeOfFlight()<<" trackid "<<cfi->trackId() <<" bunchcr "<<cfi.bunch()<<" trigger "<<cfi.getTrigger()<<", from EncodedEventId: "<<cfi->eventId().bunchCrossing()<<" " <<cfi->eventId().event() <<" bcr from MixCol "<<cfi.bunch()<<std::endl;
      }
      count5++;
    }
  }
  


}


// ------------ method called once each job just before starting event loop  ------------
void 
TestMixedSource::beginJob()
{ 
   outputFile.open("test.log");
   if (!outputFile.is_open())
   {
     std::cout << "Unable to open file!" << std::endl;
   } 
   
}

// ------------ method called once each job just after ending the event loop  ------------
void 
TestMixedSource::endJob() {
  
  outputFile.close();
  
  char t[30];
  sprintf(t,"%s",fileName_.c_str());
  std::cout << " fileName = " << t << std::endl; 
  histFile_=new TFile(t,"RECREATE");
  histTrack_bunchPileups_->Write();
  histTrack_bunchSignal_->Write();
  histVertex_bunch_->Write();
  histPCaloHit_bunch_->Write();
  histPSimHit_bunchPileups_TrackerHitsTECHighTof_->Write();
  histPSimHit_bunchSignal_TrackerHitsTECHighTof_->Write();
  tofhist_sig_->Write();
  tofhist_->Write();
  histPSimHit_bunchSignal_MuonCSCHits_->Write();
  histPSimHit_bunchPileups_MuonCSCHits_->Write();
  histHepMCProduct_bunch_->Write();
  histFile_->Write();
  histFile_->Close();

}
}//edm