Ticket #639: Example3.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/Example3.C'("delphes_output.root")'
*/
using namespace TMath;
//------------------------------------------------------------------------------

struct TestPlots
{
  TH1 *fJetPT;
  TH1 *fdeltaR;

  TH1 *fTowerPT;
  TH1 *fTowerE;
  TH1 *fdeltaE;

};

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

class ExRootResult;
class ExRootTreeReader;

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

void BookHistograms(ExRootResult *result, TestPlots *plots)
{
  TLegend *legend;
  TPaveText *comment;

  plots->fTowerPT = result->AddHist1D(
                                      "TowerPt", "Pt",
                                      "Pt", "number of towers",
                                      100, 1, 4000);

  plots->fTowerE = result->AddHist1D(
                                     "TowerE", "Energy",
                                     "Energy", "number of towers",
                                     100, 1, 4000);

  plots->fJetPT = result->AddHist1D(
                                    "jet_pt", "Jet Pt",
                                    "Jet Pt", "number of jets",
                                    100, 2900, 5000);

  plots->fdeltaR=result->AddHist1D(
                                   "deltaR", "Tower Distance from Jet",
                                   "Delta R", "number of towers",
                                   100, .001, 5);

  plots->fdeltaE=result->AddHist1D(
                                   "deltaE", "Energy Difference Between Jet and Towers",
                                   "Delta E", "Frequency",
                                   100, -1000, 4000);

}

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

void AnalyseEvents(ExRootTreeReader *treeReader, TestPlots *plots)
{
  TClonesArray *branchParticle = treeReader->UseBranch("Particle");
  TClonesArray *branchTower = treeReader->UseBranch("Tower");
  TClonesArray *branchJet = treeReader->UseBranch("Jet");

  Long64_t allEntries = treeReader->GetEntries();

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

  GenParticle *particle;
  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);

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

          // Loop over all jet's constituents if jetPt>3000
          if (jet->PT>3000) {

            cout<<"Jet pt: "<<jet->PT<<", eta: "<<jet->Eta<<", phi: "<<jet->Phi<<", Energy:"<<jet->P4().Energy()<<endl;

            double particleE=0;
            double towerE=0;
            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();
                    particleE +=particle->E;

                  }
        

                  
                else if(object->IsA() == Tower::Class())
                  {
                    tower = (Tower*) object;
                    //cout << "    Tower pt: " << tower->ET << ", eta: " << tower->Eta << ", phi: " << tower->Phi << endl;
                    momentum += tower->P4();
                    plots->fTowerPT->Fill(tower->P4().Perp());
                    plots->fTowerE->Fill(tower->E);

                    double dEta=abs(jet->Eta-tower->Eta);
                    double dPhi=abs(jet->Phi-tower->Phi);
                    //if (dPhi>Pi()) { dPhi=Pi()-dPhi; }
                    double dR=Sqrt(Power(dEta,2)+Power(dPhi,2));
                    plots->fdeltaR->Fill(dR);
                    towerE +=tower->E;

                  }

              }
        
          cout<<"Particle Energy Sum: "<<particleE<<endl;
          cout<<"Tower Energy Sum: "<<towerE<<endl;     
          plots->fdeltaE->Fill(jet->P4().Energy()-(towerE+particleE));
         }
          if (jet->PT>3000) {
            plots->fJetPT->Fill(jet->PT);
          }
        }
    }
}

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

void PrintHistograms(ExRootResult *result, TestPlots *plots)
{
  result->Print("png");
}

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

void Example3()
{
  gSystem->Load("libDelphes");

  const char* inputFile="tev14_pythia8_qcd_pt2500_001.root"; 
 
  TChain *chain = new TChain("Delphes");
  chain->Add(inputFile);

  ExRootTreeReader *treeReader = new ExRootTreeReader(chain);
  ExRootResult *result = new ExRootResult();

  TestPlots *plots = new TestPlots;

  BookHistograms(result, plots);


  AnalyseEvents(treeReader, plots);

  PrintHistograms(result, plots);

  result->Write("results.root");

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

  delete plots;
  delete result;
  delete treeReader;
  delete chain;
}

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