snippets/FullChainBnlFarmScript/example_macro.C

0001 #include "TFile.h"
0002 #include "TH2D.h"
0003 #include "TTreeReader.h"
0004 #include "TTreeReaderArray.h"
0005 #include <cmath>
0006 #include <iostream>
0007
0008 /**
0009  * ROOT macro for analyzing MC particle distributions
0010  * Creates pt vs eta histogram from Monte Carlo particle data
0011  *
0012  * @param inFileName  Input ROOT file containing MC data
0013  * @param outFileName Output ROOT file for histograms
0014  */
0015 void example_macro(TString inFileName = "input.root",
0016                    TString outFileName = "ana.root") {
0017
0018     // Open input file containing the event tree
0019     TFile* file = new TFile(inFileName);
0020     if (!file || file->IsZombie()) {
0021         std::cerr << "Error: Cannot open input file " << inFileName << std::endl;
0022         return;
0023     }
0024
0025     // Setup tree reader for MC particle data
0026     TTreeReader myReader("events", file);
0027
0028     // Define data containers for particle momentum components
0029     TTreeReaderArray<double> px_mc(myReader, "MCParticles.momentum.x");
0030     TTreeReaderArray<double> py_mc(myReader, "MCParticles.momentum.y");
0031     TTreeReaderArray<double> pz_mc(myReader, "MCParticles.momentum.z");
0032     TTreeReaderArray<int>    status(myReader, "MCParticles.generatorStatus");
0033     TTreeReaderArray<int>    pdg(myReader, "MCParticles.PDG");
0034
0035     // Create output file and histogram
0036     TFile* output = new TFile(outFileName, "recreate");
0037     TH2D* hPtEtaMc = new TH2D("hPtEtaMc",
0038                               "MC Particle Distribution;p_{t,gen} (GeV/c);#eta_{gen}",
0039                               1500, 0.0, 10.0,    // pt bins: 1500 bins from 0 to 10 GeV/c
0040                               200, -5.0, 5.0);    // eta bins: 200 bins from -5 to 5
0041
0042     // Event loop
0043     int nEvents = myReader.GetEntries();
0044     std::cout << "Total Events: " << nEvents << std::endl;
0045
0046     for (int iEvent = 0; iEvent < nEvents; ++iEvent) {
0047         myReader.SetEntry(iEvent);
0048
0049         // Progress indicator
0050         if (iEvent % 1000 == 0) {
0051             std::cout << "Processing Event: " << iEvent << std::endl;
0052         }
0053
0054         // Loop over particles in current event
0055         for (int iParticle = 0; iParticle < status.GetSize(); ++iParticle) {
0056
0057             // Only consider stable particles (status == 1)
0058             if (status[iParticle] != 1) {
0059                 continue;
0060             }
0061
0062             // Calculate kinematic variables
0063             double px = px_mc[iParticle];
0064             double py = py_mc[iParticle];
0065             double pz = pz_mc[iParticle];
0066
0067             // Total momentum magnitude
0068             double p_mc = std::sqrt(px*px + py*py + pz*pz);
0069
0070             // Transverse momentum
0071             double pt_mc = std::sqrt(px*px + py*py);
0072
0073             // Pseudorapidity (eta)
0074             double theta_mc = std::acos(pz / p_mc);
0075             double eta_mc = -std::log(std::tan(theta_mc / 2.0));
0076
0077             // Optional: Azimuthal angle (phi) - currently calculated but not used
0078             // double phi_mc = std::atan2(py, px);  // Use atan2 for proper quadrant
0079
0080             // Fill histogram
0081             hPtEtaMc->Fill(pt_mc, eta_mc);
0082
0083         } // End particle loop
0084
0085     } // End event loop
0086
0087     // Write output
0088     output->cd();
0089     hPtEtaMc->Write();
0090     output->Save();
0091     output->Close();
0092
0093     // Clean up
0094     file->Close();
0095     delete file;
0096     delete output;
0097
0098     std::cout << "Analysis complete. Output saved to: " << outFileName << std::endl;
0099 }