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