AnalysisRootpleProducer.cc

Go to the documentation of this file.

// Authors: F. Ambroglini, L. Fano'
#include <QCDAnalysis/UEAnalysis/interface/AnalysisRootpleProducer.h>
#include "FWCore/Common/interface/TriggerNames.h"

using namespace edm;
using namespace std;
using namespace reco;

class GreaterPt{
public:
  bool operator()( const math::XYZTLorentzVector& a, const math::XYZTLorentzVector& b) {
    return a.pt() > b.pt();
  }
};

class GenJetSort{
public:
  bool operator()(const GenJet& a, const GenJet& b) {
    return a.pt() > b.pt();
  }
};

class BasicJetSort{
public:
  bool operator()(const BasicJet& a, const BasicJet& b) {
    return a.pt() > b.pt();
  }
};

class CaloJetSort{
public:
  bool operator()(const CaloJet& a, const CaloJet& b) {
    return a.pt() > b.pt();
  }
};


void AnalysisRootpleProducer::store(){

  AnalysisTree->Fill();

  NumberMCParticles=0;
  NumberTracks=0;
  NumberInclusiveJet=0;
  NumberChargedJet=0;
  NumberTracksJet=0;
  NumberCaloJet=0;
}

void AnalysisRootpleProducer::fillEventInfo(int e){
  EventKind = e;
}

void AnalysisRootpleProducer::fillMCParticles(float p, float pt, float eta, float phi){
  MomentumMC[NumberMCParticles]=p;
  TransverseMomentumMC[NumberMCParticles]=pt;
  EtaMC[NumberMCParticles]=eta;
  PhiMC[NumberMCParticles]=phi;
  NumberMCParticles++;
}

void AnalysisRootpleProducer::fillTracks(float p, float pt, float eta, float phi){
  MomentumTK[NumberTracks]=p;
  TransverseMomentumTK[NumberTracks]=pt;
  EtaTK[NumberTracks]=eta;
  PhiTK[NumberTracks]=phi;
  NumberTracks++;
}

void AnalysisRootpleProducer::fillInclusiveJet(float p, float pt, float eta, float phi){
  MomentumIJ[NumberInclusiveJet]=p;
  TransverseMomentumIJ[NumberInclusiveJet]=pt;
  EtaIJ[NumberInclusiveJet]=eta;
  PhiIJ[NumberInclusiveJet]=phi;
  NumberInclusiveJet++;
}

void AnalysisRootpleProducer::fillChargedJet(float p, float pt, float eta, float phi){
  MomentumCJ[NumberChargedJet]=p;
  TransverseMomentumCJ[NumberChargedJet]=pt;
  EtaCJ[NumberChargedJet]=eta;
  PhiCJ[NumberChargedJet]=phi;
  NumberChargedJet++;
}

void AnalysisRootpleProducer::fillTracksJet(float p, float pt, float eta, float phi){
  MomentumTJ[NumberTracksJet]=p;
  TransverseMomentumTJ[NumberTracksJet]=pt;
  EtaTJ[NumberTracksJet]=eta;
  PhiTJ[NumberTracksJet]=phi;
  NumberTracksJet++;
}

void AnalysisRootpleProducer::fillCaloJet(float p, float pt, float eta, float phi){
  MomentumEHJ[NumberCaloJet]=p;
  TransverseMomentumEHJ[NumberCaloJet]=pt;
  EtaEHJ[NumberCaloJet]=eta;
  PhiEHJ[NumberCaloJet]=phi;
  NumberCaloJet++;
}

AnalysisRootpleProducer::AnalysisRootpleProducer( const ParameterSet& pset )
{
  // flag to ignore gen-level analysis
  onlyRECO = pset.getUntrackedParameter<bool>("OnlyRECO",false);

  // particle, track and jet collections
  mcEventToken = mayConsume<edm::HepMCProduct>(pset.getUntrackedParameter<InputTag>("MCEvent",std::string("")));
  genJetCollToken = mayConsume<reco::GenJetCollection>(pset.getUntrackedParameter<InputTag>("GenJetCollectionName",std::string("")));
  chgJetCollToken = mayConsume<reco::GenJetCollection>(pset.getUntrackedParameter<InputTag>("ChgGenJetCollectionName",std::string("")));
  tracksJetCollToken = consumes<reco::BasicJetCollection>(pset.getUntrackedParameter<InputTag>("TracksJetCollectionName",std::string("")));
  recoCaloJetCollToken = consumes<reco::CaloJetCollection>(pset.getUntrackedParameter<InputTag>("RecoCaloJetCollectionName",std::string("")));
  chgGenPartCollToken = mayConsume<std::vector<reco::GenParticle> >(pset.getUntrackedParameter<InputTag>("ChgGenPartCollectionName",std::string("")));
  tracksCollToken = consumes<reco::CandidateCollection>(pset.getUntrackedParameter<InputTag>("TracksCollectionName",std::string("")));

  // trigger results
  triggerResultsToken = consumes<edm::TriggerResults>(pset.getParameter<InputTag>("triggerResults"));
  //   hltFilterTag      = pset.getParameter<InputTag>("hltFilter");
  //   triggerName       = pset.getParameter<InputTag>("triggerName");

  piG = acos(-1.);
  NumberMCParticles=0;
  NumberTracks=0;
  NumberInclusiveJet=0;
  NumberChargedJet=0;
  NumberTracksJet=0;
  NumberCaloJet=0;
}

void AnalysisRootpleProducer::beginJob()
{

  // use TFileService for output to root file
  AnalysisTree = fs->make<TTree>("AnalysisTree","MBUE Analysis Tree ");

  AnalysisTree->Branch("EventKind",&EventKind,"EventKind/I");

  // store p, pt, eta, phi for particles and jets

  // GenParticles at hadron level
  AnalysisTree->Branch("NumberMCParticles",&NumberMCParticles,"NumberMCParticles/I");
  AnalysisTree->Branch("MomentumMC",MomentumMC,"MomentumMC[NumberMCParticles]/F");
  AnalysisTree->Branch("TransverseMomentumMC",TransverseMomentumMC,"TransverseMomentumMC[NumberMCParticles]/F");
  AnalysisTree->Branch("EtaMC",EtaMC,"EtaMC[NumberMCParticles]/F");
  AnalysisTree->Branch("PhiMC",PhiMC,"PhiMC[NumberMCParticles]/F");

  // tracks
  AnalysisTree->Branch("NumberTracks",&NumberTracks,"NumberTracks/I");
  AnalysisTree->Branch("MomentumTK",MomentumTK,"MomentumTK[NumberTracks]/F");
  AnalysisTree->Branch("TrasverseMomentumTK",TransverseMomentumTK,"TransverseMomentumTK[NumberTracks]/F");
  AnalysisTree->Branch("EtaTK",EtaTK,"EtaTK[NumberTracks]/F");
  AnalysisTree->Branch("PhiTK",PhiTK,"PhiTK[NumberTracks]/F");

  // GenJets
  AnalysisTree->Branch("NumberInclusiveJet",&NumberInclusiveJet,"NumberInclusiveJet/I");
  AnalysisTree->Branch("MomentumIJ",MomentumIJ,"MomentumIJ[NumberInclusiveJet]/F");
  AnalysisTree->Branch("TrasverseMomentumIJ",TransverseMomentumIJ,"TransverseMomentumIJ[NumberInclusiveJet]/F");
  AnalysisTree->Branch("EtaIJ",EtaIJ,"EtaIJ[NumberInclusiveJet]/F");
  AnalysisTree->Branch("PhiIJ",PhiIJ,"PhiIJ[NumberInclusiveJet]/F");

  // jets from charged GenParticles
  AnalysisTree->Branch("NumberChargedJet",&NumberChargedJet,"NumberChargedJet/I");
  AnalysisTree->Branch("MomentumCJ",MomentumCJ,"MomentumCJ[NumberChargedJet]/F");
  AnalysisTree->Branch("TrasverseMomentumCJ",TransverseMomentumCJ,"TransverseMomentumCJ[NumberChargedJet]/F");
  AnalysisTree->Branch("EtaCJ",EtaCJ,"EtaCJ[NumberChargedJet]/F");
  AnalysisTree->Branch("PhiCJ",PhiCJ,"PhiCJ[NumberChargedJet]/F");

  // jets from tracks
  AnalysisTree->Branch("NumberTracksJet",&NumberTracksJet,"NumberTracksJet/I");
  AnalysisTree->Branch("MomentumTJ",MomentumTJ,"MomentumTJ[NumberTracksJet]/F");
  AnalysisTree->Branch("TrasverseMomentumTJ",TransverseMomentumTJ,"TransverseMomentumTJ[NumberTracksJet]/F");
  AnalysisTree->Branch("EtaTJ",EtaTJ,"EtaTJ[NumberTracksJet]/F");
  AnalysisTree->Branch("PhiTJ",PhiTJ,"PhiTJ[NumberTracksJet]/F");

  // jets from calorimeter towers
  AnalysisTree->Branch("NumberCaloJet",&NumberCaloJet,"NumberCaloJet/I");
  AnalysisTree->Branch("MomentumEHJ",MomentumEHJ,"MomentumEHJ[NumberCaloJet]/F");
  AnalysisTree->Branch("TrasverseMomentumEHJ",TransverseMomentumEHJ,"TransverseMomentumEHJ[NumberCaloJet]/F");
  AnalysisTree->Branch("EtaEHJ",EtaEHJ,"EtaEHJ[NumberCaloJet]/F");
  AnalysisTree->Branch("PhiEHJ",PhiEHJ,"PhiEHJ[NumberCaloJet]/F");


  // alternative storage method:
  // save TClonesArrays of TLorentzVectors
  // i.e. store 4-vectors of particles and jets

  MonteCarlo = new TClonesArray("TLorentzVector", 10000);
  AnalysisTree->Branch("MonteCarlo", "TClonesArray", &MonteCarlo, 128000, 0);

  Track = new TClonesArray("TLorentzVector", 10000);
  AnalysisTree->Branch("Track", "TClonesArray", &Track, 128000, 0);

  InclusiveJet = new TClonesArray("TLorentzVector", 10000);
  AnalysisTree->Branch("InclusiveJet", "TClonesArray", &InclusiveJet, 128000, 0);

  ChargedJet = new TClonesArray("TLorentzVector", 10000);
  AnalysisTree->Branch("ChargedJet", "TClonesArray", &ChargedJet, 128000, 0);

  TracksJet = new TClonesArray("TLorentzVector", 10000);
  AnalysisTree->Branch("TracksJet", "TClonesArray", &TracksJet, 128000, 0);

  CalorimeterJet = new TClonesArray("TLorentzVector", 10000);
  AnalysisTree->Branch("CalorimeterJet", "TClonesArray", &CalorimeterJet, 128000, 0);

  acceptedTriggers = new TClonesArray("TObjString", 10000);
  AnalysisTree->Branch("acceptedTriggers", "TClonesArray", &acceptedTriggers, 128000, 0);

}


void AnalysisRootpleProducer::analyze( const Event& e, const EventSetup& )
{

  e.getByToken( triggerResultsToken, triggerResults );
  const edm::TriggerNames & triggerNames = e.triggerNames(*triggerResults);

  acceptedTriggers->Clear();
  unsigned int iAcceptedTriggers( 0 );
  if ( triggerResults.product()->wasrun() )
    {
      //cout << "at least one path out of " << triggerResults.product()->size() << " ran? " << triggerResults.product()->wasrun() << endl;

      if ( triggerResults.product()->accept() )
    {
      //cout << endl << "at least one path accepted? " << triggerResults.product()->accept() << endl;

      const unsigned int n_TriggerResults( triggerResults.product()->size() );
      for ( unsigned int itrig( 0 ); itrig < n_TriggerResults; ++itrig )
        {
          if ( triggerResults.product()->accept( itrig ) )
        {
          //cout << "path " << triggerNames.triggerName( itrig );
          //cout << ", module index " << triggerResults.product()->index( itrig );
          //cout << ", state (Ready = 0, Pass = 1, Fail = 2, Exception = 3) " << triggerResults.product()->state( itrig );
          //cout << ", accept " << triggerResults.product()->accept( itrig );
          //cout << endl;

          // save name of accepted trigger path
          new((*acceptedTriggers)[iAcceptedTriggers]) TObjString( (triggerNames.triggerName( itrig )).c_str() );
          ++iAcceptedTriggers;
        }
        }
    }
    }

  // gen level analysis
  // skipped, if onlyRECO flag set to true

  if(!onlyRECO){

    e.getByToken( mcEventToken       , EvtHandle        );
    e.getByToken( chgGenPartCollToken, CandHandleMC     );
    e.getByToken( chgJetCollToken    , ChgGenJetsHandle );
    e.getByToken( genJetCollToken    , GenJetsHandle    );

    const HepMC::GenEvent* Evt = EvtHandle->GetEvent() ;

    EventKind = Evt->signal_process_id();

    std::vector<math::XYZTLorentzVector> GenPart;
    std::vector<GenJet> ChgGenJetContainer;
    std::vector<GenJet> GenJetContainer;

    GenPart.clear();
    ChgGenJetContainer.clear();
    GenJetContainer.clear();
    MonteCarlo->Clear();
    InclusiveJet->Clear();
    ChargedJet->Clear();

    // jets from charged particles at hadron level
    if (!ChgGenJetsHandle->empty()){

      for ( GenJetCollection::const_iterator it(ChgGenJetsHandle->begin()), itEnd(ChgGenJetsHandle->end());
        it!=itEnd; ++it)
    {
      ChgGenJetContainer.push_back(*it);
    }

      std::stable_sort(ChgGenJetContainer.begin(),ChgGenJetContainer.end(),GenJetSort());

      std::vector<GenJet>::const_iterator it(ChgGenJetContainer.begin()), itEnd(ChgGenJetContainer.end());
      for ( int iChargedJet(0); it != itEnd; ++it, ++iChargedJet)
    {
      fillChargedJet(it->p(),it->pt(),it->eta(),it->phi());
      new((*ChargedJet)[iChargedJet]) TLorentzVector(it->px(), it->py(), it->pz(), it->energy());
    }
    }


    // GenJets
    if (!GenJetsHandle->empty()){

      for ( GenJetCollection::const_iterator it(GenJetsHandle->begin()), itEnd(GenJetsHandle->end());
        it!=itEnd; ++it )
    {
      GenJetContainer.push_back(*it);
    }

      std::stable_sort(GenJetContainer.begin(),GenJetContainer.end(),GenJetSort());

      std::vector<GenJet>::const_iterator it(GenJetContainer.begin()), itEnd(GenJetContainer.end());
      for ( int iInclusiveJet(0); it != itEnd; ++it, ++iInclusiveJet)
    {
      fillInclusiveJet(it->p(),it->pt(),it->eta(),it->phi());
      new((*InclusiveJet)[iInclusiveJet]) TLorentzVector(it->px(), it->py(), it->pz(), it->energy());
    }
    }


    // hadron level particles
    if (!CandHandleMC->empty()){

      for (vector<GenParticle>::const_iterator it(CandHandleMC->begin()), itEnd(CandHandleMC->end());
       it != itEnd;it++)
    {
      GenPart.push_back(it->p4());
    }

      std::stable_sort(GenPart.begin(),GenPart.end(),GreaterPt());

      std::vector<math::XYZTLorentzVector>::const_iterator it(GenPart.begin()), itEnd(GenPart.end());
      for( int iMonteCarlo(0); it != itEnd; ++it, ++iMonteCarlo )
    {
      fillMCParticles(it->P(),it->Pt(),it->Eta(),it->Phi());
      new((*MonteCarlo)[iMonteCarlo]) TLorentzVector(it->Px(), it->Py(), it->Pz(), it->E());
    }
    }

  }


  // reco level analysis

  e.getByToken( tracksCollToken     , CandHandleRECO     );
  e.getByToken( recoCaloJetCollToken, RecoCaloJetsHandle );
  e.getByToken( tracksJetCollToken  , TracksJetsHandle   );

  std::vector<math::XYZTLorentzVector> Tracks;
  std::vector<BasicJet> TracksJetContainer;
  std::vector<CaloJet> RecoCaloJetContainer;

  Tracks.clear();
  TracksJetContainer.clear();
  RecoCaloJetContainer.clear();

  Track->Clear();
  TracksJet->Clear();
  CalorimeterJet->Clear();

  if(!RecoCaloJetsHandle->empty())
    {
    for(CaloJetCollection::const_iterator it(RecoCaloJetsHandle->begin()), itEnd(RecoCaloJetsHandle->end());
    it!=itEnd;++it)
      {
    RecoCaloJetContainer.push_back(*it);
      }
    std::stable_sort(RecoCaloJetContainer.begin(),RecoCaloJetContainer.end(),CaloJetSort());

    std::vector<CaloJet>::const_iterator it(RecoCaloJetContainer.begin()), itEnd(RecoCaloJetContainer.end());
    for( int iCalorimeterJet(0); it != itEnd; ++it, ++iCalorimeterJet)
      {
    fillCaloJet(it->p(),it->pt(),it->eta(),it->phi());
    new((*CalorimeterJet)[iCalorimeterJet]) TLorentzVector(it->px(), it->py(), it->pz(), it->energy());
      }
    }

  if(!TracksJetsHandle->empty())
    {
      for(BasicJetCollection::const_iterator it(TracksJetsHandle->begin()), itEnd(TracksJetsHandle->end());
      it!=itEnd;++it)
    {
      TracksJetContainer.push_back(*it);
    }
      std::stable_sort(TracksJetContainer.begin(),TracksJetContainer.end(),BasicJetSort());

      std::vector<BasicJet>::const_iterator it(TracksJetContainer.begin()), itEnd(TracksJetContainer.end());
      for(int iTracksJet(0); it != itEnd; ++it, ++iTracksJet)
    {
      fillTracksJet(it->p(),it->pt(),it->eta(),it->phi());
      new((*TracksJet)[iTracksJet]) TLorentzVector(it->px(), it->py(), it->pz(), it->energy());
    }
    }

  if(!CandHandleRECO->empty())
    {
      for(CandidateCollection::const_iterator it(CandHandleRECO->begin()), itEnd(CandHandleRECO->end());
      it!=itEnd;++it)
    {
      Tracks.push_back(it->p4());
    }
      std::stable_sort(Tracks.begin(),Tracks.end(),GreaterPt());

      std::vector<math::XYZTLorentzVector>::const_iterator it( Tracks.begin()), itEnd(Tracks.end());
      for(int iTracks(0); it != itEnd; ++it, ++iTracks)
    {
      fillTracks(it->P(),it->Pt(),it->Eta(),it->Phi());
      new ((*Track)[iTracks]) TLorentzVector(it->Px(), it->Py(), it->Pz(), it->E());
    }
    }

  store();
}

void AnalysisRootpleProducer::endJob()
{
}