AnalysisRootpleProducerOnlyMC.cc

Go to the documentation of this file.

// Authors: F. Ambroglini, L. Fano'
#include <QCDAnalysis/UEAnalysis/interface/AnalysisRootpleProducerOnlyMC.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();
  }
};


void AnalysisRootpleProducerOnlyMC::store()
{
  AnalysisTree->Fill();

  NumberMCParticles=0;
  NumberInclusiveJet=0;
  NumberChargedJet=0;
}

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

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

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

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

AnalysisRootpleProducerOnlyMC::AnalysisRootpleProducerOnlyMC( const ParameterSet& pset )
{
  mcEventToken = consumes<edm::HepMCProduct>(pset.getUntrackedParameter<InputTag>("MCEvent",std::string("")));
  genJetCollToken = consumes<reco::GenJetCollection>(pset.getUntrackedParameter<InputTag>("GenJetCollectionName",std::string("")));
  chgJetCollToken = consumes<reco::GenJetCollection>(pset.getUntrackedParameter<InputTag>("ChgGenJetCollectionName",std::string("")));
  chgGenPartCollToken = consumes<std::vector<reco::GenParticle> >(pset.getUntrackedParameter<InputTag>("ChgGenPartCollectionName",std::string("")));

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


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

  // process type
  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");

  // 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");


  // 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);

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

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


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

  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();

  ChargedJet->Clear();
  InclusiveJet->Clear();
  MonteCarlo->Clear();

  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());
      }
  }

  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());
      }
  }

  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());
      }
  }

  store();
}

void AnalysisRootpleProducerOnlyMC::endJob()
{
}