Ticket #1369: jet_constits.C

/*
This macro shows how to access the particle-level reference for reconstructed objects.
It is also shown how to loop over the jet constituents.

root -l examples/jet_constits.C'("delphes_output.root")'
*/

#ifdef __CLING__
R__LOAD_LIBRARY(libDelphes)
#include "classes/DelphesClasses.h"
#include "external/ExRootAnalysis/ExRootTreeReader.h"
#include "external/ExRootAnalysis/ExRootResult.h"
#else
class ExRootTreeReader;
class ExRootResult;
#endif

void AnalyseEvents(ExRootTreeReader *treeReader)
{
    {//output files have a scope?!?
    ofstream jetOut("jetPts.csv");
    ofstream trackOut("trackPts.csv");
    ofstream towerOut("towerPts.csv");
    ofstream genOut("genPts.csv");
    TClonesArray *branchParticle = treeReader->UseBranch("Particle");
    //TClonesArray *branchElectron = treeReader->UseBranch("Electron");
    //TClonesArray *branchPhoton = treeReader->UseBranch("Photon");
    //TClonesArray *branchMuon = treeReader->UseBranch("Muon");
    //~~~~~ Uncommenting this give tracks
    TClonesArray *branchEFlowTrack = treeReader->UseBranch("EFlowTrack");
    //~~~~~ Uncommenting this give towers
    TClonesArray *branchEFlowPhoton = treeReader->UseBranch("EFlowPhoton");
    //~~~~~ Uncommenting this give towers
    TClonesArray *branchEFlowNeutralHadron = treeReader->UseBranch("EFlowNeutralHadron");
    //~~~~~~~ Why are the existance of unused pointers changing the output?
    TClonesArray *branchJet = treeReader->UseBranch("Jet");

    Long64_t allEntries = treeReader->GetEntries();

    cout << "** Chain contains " << allEntries << " events" << endl;

    GenParticle *particle;
    //Electron *electron;
    //Photon *photon;
    //Muon *muon;

    Track *track;
    Tower *tower;

    Jet *jet;
    TObject *object;

    TLorentzVector momentum;

    //Float_t Eem, Ehad;
    //Bool_t skip;

    Long64_t entry;

    Int_t i, j; //, pdgCode;

    // Loop over all events
    for(entry = 0; entry < allEntries; ++entry)
    {
      // Load selected branches with data from specified event
      treeReader->ReadEntry(entry);

      // // Loop over all electrons in event
      // for(i = 0; i < branchElectron->GetEntriesFast(); ++i)
      // {
      //   electron = (Electron*) branchElectron->At(i);
      //   particle = (GenParticle*) electron->Particle.GetObject();

      // }

      // // Loop over all photons in event
      // for(i = 0; i < branchPhoton->GetEntriesFast(); ++i)
      // {
      //   photon = (Photon*) branchPhoton->At(i);

      //   // skip photons with references to multiple particles
      //   if(photon->Particles.GetEntriesFast() != 1) continue;

      //   particle = (GenParticle*) photon->Particles.At(0);
      //       
      // }

      // // Loop over all muons in event
      // for(i = 0; i < branchMuon->GetEntriesFast(); ++i)
      // {
      //   muon = (Muon*) branchMuon->At(i);
      //   particle = (GenParticle*) muon->Particle.GetObject();

      // }

      cout << "--  New event -- " << endl;

      // Loop over all jets in event
      for(i = 0; i < branchJet->GetEntriesFast(); ++i)
      {
        jet = (Jet*) branchJet->At(i);

        momentum.SetPxPyPzE(0.0, 0.0, 0.0, 0.0);

        cout<<"Looping over jet constituents. Jet pt: "<<jet->PT<<", eta: "<<jet->Eta<<", phi: "<<jet->Phi<<endl;
        jetOut << jet->PT;

        // Loop over all jet's constituents
        for(j = 0; j < jet->Constituents.GetEntriesFast(); ++j)
        {
          object = jet->Constituents.At(j);

          // Check if the constituent is accessible
          if(object == 0) continue;

          if(object->IsA() == GenParticle::Class())
          {
            particle = (GenParticle*) object;
            cout << "    GenPart pt: " << particle->PT << ", eta: " << particle->Eta << ", phi: " << particle->Phi << endl;
            momentum += particle->P4();
            genOut << particle->PT << ",";
          }
          else if(object->IsA() == Track::Class())
          {
            track = (Track*) object;
            cout << "    Track pt: " << track->PT << ", eta: " << track->Eta << ", phi: " << track->Phi << endl;
            momentum += track->P4();
            trackOut << track->PT << ",";
          }
          else if(object->IsA() == Tower::Class())
          {
            tower = (Tower*) object;
            cout << "    Tower pt: " << tower->ET << ", eta: " << tower->Eta << ", phi: " << tower->Phi << endl;
            momentum += tower->P4();
            towerOut << tower->ET << ",";
          }
        }
        jetOut << endl;
        towerOut << endl;
        genOut << endl;
        trackOut << endl;
      }
    }
    }
}

//------------------------------------------------------------------------------

void jet_constits(const char *inputFile)
{
    gSystem->Load("libDelphes");

    TChain *chain = new TChain("Delphes");
    chain->Add(inputFile);

    ExRootTreeReader *treeReader = new ExRootTreeReader(chain);

    AnalyseEvents(treeReader);


    cout << "** Exiting..." << endl;

    delete treeReader;
    delete chain;
}

//------------------------------------------------------------------------------