CTPPS/plugins/CTPPSAcceptancePlotter.cc

0001 /****************************************************************************
0002  * Authors:
0003  *   Jan Kašpar
0004  ****************************************************************************/
0005
0006 #include "FWCore/Framework/interface/Frameworkfwd.h"
0007 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0008
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011
0012 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0013
0014 #include "DataFormats/CTPPSDetId/interface/CTPPSDetId.h"
0015 #include "DataFormats/CTPPSReco/interface/CTPPSLocalTrackLite.h"
0016 #include "DataFormats/CTPPSReco/interface/CTPPSLocalTrackLiteFwd.h"
0017
0018 #include "TFile.h"
0019 #include "TH1D.h"
0020 #include "TH2D.h"
0021
0022 #include <map>
0023 #include <set>
0024 #include <string>
0025
0026 //----------------------------------------------------------------------------------------------------
0027
0028 class CTPPSAcceptancePlotter : public edm::one::EDAnalyzer<> {
0029 public:
0030   explicit CTPPSAcceptancePlotter(const edm::ParameterSet &);
0031
0032 private:
0033   void beginJob() override {}
0034   void analyze(const edm::Event &, const edm::EventSetup &) override;
0035   void endJob() override;
0036
0037   edm::EDGetTokenT<edm::HepMCProduct> tokenHepMC_;
0038   edm::EDGetTokenT<CTPPSLocalTrackLiteCollection> tokenTracks_;
0039
0040   unsigned int rpId_45_N_, rpId_45_F_, rpId_56_N_, rpId_56_F_;
0041
0042   std::string outputFile_;
0043
0044   struct SingleArmPlots {
0045     std::unique_ptr<TH1D> h_xi_all, h_xi_acc;
0046     SingleArmPlots() : h_xi_all(new TH1D("", ";#xi", 100, 0., 0.25)), h_xi_acc(new TH1D("", ";#xi", 100, 0., 0.25)) {}
0047
0048     void fill(double xi, bool acc) {
0049       h_xi_all->Fill(xi);
0050       if (acc)
0051         h_xi_acc->Fill(xi);
0052     }
0053
0054     void write() const {
0055       h_xi_all->Write("h_xi_all");
0056       h_xi_acc->Write("h_xi_acc");
0057
0058       auto h_xi_rat = std::make_unique<TH1D>(*h_xi_acc);
0059       h_xi_rat->Divide(h_xi_all.get());
0060       h_xi_rat->Write("h_xi_rat");
0061     }
0062   };
0063
0064   std::vector<std::set<unsigned int>> singleArmConfigurations;
0065   std::map<std::set<unsigned int>, SingleArmPlots> singleArmPlots;
0066
0067   struct DoubleArmPlots {
0068     std::unique_ptr<TH1D> h_m_all, h_m_acc;
0069     std::unique_ptr<TH2D> h2_xi_45_vs_xi_56_all, h2_xi_45_vs_xi_56_acc;
0070     std::unique_ptr<TH2D> h2_y_vs_m_all, h2_y_vs_m_acc;
0071
0072     DoubleArmPlots()
0073         : h_m_all(new TH1D("", ";m   (GeV)", 100, 0., 2500.)),
0074           h_m_acc(new TH1D("", ";m   (GeV)", 100, 0., 2500.)),
0075           h2_xi_45_vs_xi_56_all(new TH2D("", ";xi_56;xi_45", 25, 0., 0.25, 25, 0., 0.25)),
0076           h2_xi_45_vs_xi_56_acc(new TH2D("", ";xi_56;xi_45", 25, 0., 0.25, 25, 0., 0.25)),
0077           h2_y_vs_m_all(new TH2D("", ";m   (GeV);y", 25, 0., 2500., 25, -1.5, +1.5)),
0078           h2_y_vs_m_acc(new TH2D("", ";m   (GeV);y", 25, 0., 2500., 25, -1.5, +1.5)) {}
0079
0080     void fill(double xi_45, double xi_56, bool acc) {
0081       const double p_nom = 6500.;
0082       const double m = 2. * p_nom * sqrt(xi_45 * xi_56);
0083       const double y = log(xi_45 / xi_56) / 2.;
0084
0085       h_m_all->Fill(m);
0086       h2_xi_45_vs_xi_56_all->Fill(xi_56, xi_45);
0087       h2_y_vs_m_all->Fill(m, y);
0088
0089       if (acc) {
0090         h_m_acc->Fill(m);
0091         h2_xi_45_vs_xi_56_acc->Fill(xi_56, xi_45);
0092         h2_y_vs_m_acc->Fill(m, y);
0093       }
0094     }
0095
0096     void write() const {
0097       h_m_all->Write("h_m_all");
0098       h_m_acc->Write("h_m_acc");
0099
0100       auto h_m_rat = std::make_unique<TH1D>(*h_m_acc);
0101       h_m_rat->Divide(h_m_all.get());
0102       h_m_rat->Write("h_m_rat");
0103
0104       h2_xi_45_vs_xi_56_all->Write("h2_xi_45_vs_xi_56_all");
0105       h2_xi_45_vs_xi_56_acc->Write("h2_xi_45_vs_xi_56_acc");
0106
0107       auto h2_xi_45_vs_xi_56_rat = std::make_unique<TH2D>(*h2_xi_45_vs_xi_56_acc);
0108       h2_xi_45_vs_xi_56_rat->Divide(h2_xi_45_vs_xi_56_all.get());
0109       h2_xi_45_vs_xi_56_rat->Write("h2_xi_45_vs_xi_56_rat");
0110
0111       h2_y_vs_m_all->Write("h2_y_vs_m_all");
0112       h2_y_vs_m_acc->Write("h2_y_vs_m_acc");
0113
0114       auto h2_y_vs_m_rat = std::make_unique<TH2D>(*h2_y_vs_m_acc);
0115       h2_y_vs_m_rat->Divide(h2_y_vs_m_all.get());
0116       h2_y_vs_m_rat->Write("h2_y_vs_m_rat");
0117     }
0118   };
0119
0120   std::vector<std::set<unsigned int>> doubleArmConfigurations;
0121   std::map<std::set<unsigned int>, DoubleArmPlots> doubleArmPlots;
0122 };
0123
0124 //----------------------------------------------------------------------------------------------------
0125
0126 using namespace std;
0127 using namespace edm;
0128 using namespace HepMC;
0129
0130 //----------------------------------------------------------------------------------------------------
0131
0132 CTPPSAcceptancePlotter::CTPPSAcceptancePlotter(const edm::ParameterSet &iConfig)
0133     : tokenHepMC_(consumes<edm::HepMCProduct>(iConfig.getParameter<edm::InputTag>("tagHepMC"))),
0134       tokenTracks_(consumes<CTPPSLocalTrackLiteCollection>(iConfig.getParameter<edm::InputTag>("tagTracks"))),
0135       rpId_45_N_(iConfig.getParameter<unsigned int>("rpId_45_N")),
0136       rpId_45_F_(iConfig.getParameter<unsigned int>("rpId_45_F")),
0137       rpId_56_N_(iConfig.getParameter<unsigned int>("rpId_56_N")),
0138       rpId_56_F_(iConfig.getParameter<unsigned int>("rpId_56_F")),
0139       outputFile_(iConfig.getParameter<string>("outputFile")) {
0140   singleArmConfigurations = {
0141       {rpId_45_N_},
0142       {rpId_45_F_},
0143       {rpId_56_N_},
0144       {rpId_56_F_},
0145       {rpId_45_N_, rpId_45_F_},
0146       {rpId_56_N_, rpId_56_F_},
0147   };
0148
0149   doubleArmConfigurations = {
0150       {rpId_45_N_, rpId_56_N_},
0151       {rpId_45_F_, rpId_56_F_},
0152       {rpId_45_N_, rpId_45_F_, rpId_56_N_, rpId_56_F_},
0153   };
0154 }
0155
0156 //----------------------------------------------------------------------------------------------------
0157
0158 void CTPPSAcceptancePlotter::analyze(const edm::Event &iEvent, const edm::EventSetup &) {
0159   // get input
0160   edm::Handle<edm::HepMCProduct> hHepMC;
0161   iEvent.getByToken(tokenHepMC_, hHepMC);
0162   HepMC::GenEvent *hepMCEvent = (HepMC::GenEvent *)hHepMC->GetEvent();
0163
0164   edm::Handle<CTPPSLocalTrackLiteCollection> hTracks;
0165   iEvent.getByToken(tokenTracks_, hTracks);
0166
0167   // extract protons
0168   bool proton_45_set = false;
0169   bool proton_56_set = false;
0170   FourVector mom_45, mom_56;
0171
0172   for (auto it = hepMCEvent->particles_begin(); it != hepMCEvent->particles_end(); ++it) {
0173     const auto &part = *it;
0174
0175     // accept only stable non-beam protons
0176     if (part->pdg_id() != 2212)
0177       continue;
0178
0179     if (part->status() != 1)
0180       continue;
0181
0182     if (part->is_beam())
0183       continue;
0184
0185     const auto &mom = part->momentum();
0186
0187     if (mom.e() < 4500.)
0188       continue;
0189
0190     if (mom.z() > 0) {
0191       // 45
0192       if (proton_45_set) {
0193         LogError("CTPPSAcceptancePlotter") << "Multiple protons found in sector 45.";
0194         return;
0195       }
0196
0197       proton_45_set = true;
0198       mom_45 = mom;
0199     } else {
0200       // 56
0201       if (proton_56_set) {
0202         LogError("CTPPSAcceptancePlotter") << "Multiple protons found in sector 56.";
0203         return;
0204       }
0205
0206       proton_56_set = true;
0207       mom_56 = mom;
0208     }
0209   }
0210
0211   // stop if protons missing
0212   if (!proton_45_set || !proton_56_set)
0213     return;
0214
0215   // calculate xi's
0216   const double p_nom = 6500.;
0217   const double xi_45 = (p_nom - mom_45.e()) / p_nom;
0218   const double xi_56 = (p_nom - mom_56.e()) / p_nom;
0219
0220   // process tracks
0221   map<unsigned int, bool> trackPresent;
0222   for (const auto &trk : *hTracks) {
0223     CTPPSDetId rpId(trk.rpId());
0224     unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
0225     trackPresent[rpDecId] = true;
0226   }
0227
0228   // update plots
0229   for (const auto &rpIds : singleArmConfigurations) {
0230     bool acc = true;
0231     signed int arm = -1;
0232     for (const auto rpId : rpIds) {
0233       acc &= trackPresent[rpId];
0234       arm = rpId / 100;
0235     }
0236
0237     if (arm < 0)
0238       continue;
0239
0240     const double xi = (arm == 0) ? xi_45 : xi_56;
0241
0242     singleArmPlots[rpIds].fill(xi, acc);
0243   }
0244
0245   for (const auto &rpIds : doubleArmConfigurations) {
0246     bool acc = true;
0247     for (const auto rpId : rpIds)
0248       acc &= trackPresent[rpId];
0249
0250     doubleArmPlots[rpIds].fill(xi_45, xi_56, acc);
0251   }
0252 }
0253
0254 //----------------------------------------------------------------------------------------------------
0255
0256 void CTPPSAcceptancePlotter::endJob() {
0257   auto f_out = std::make_unique<TFile>(outputFile_.c_str(), "recreate");
0258
0259   for (const auto &p : singleArmPlots) {
0260     string dirName;
0261     for (const auto &rpId : p.first) {
0262       if (!dirName.empty())
0263         dirName += ",";
0264       dirName += Form("%u", rpId);
0265     }
0266
0267     gDirectory = f_out->mkdir(dirName.c_str());
0268     p.second.write();
0269   }
0270
0271   for (const auto &p : doubleArmPlots) {
0272     string dirName;
0273     for (const auto &rpId : p.first) {
0274       if (!dirName.empty())
0275         dirName += ",";
0276       dirName += Form("%u", rpId);
0277     }
0278
0279     gDirectory = f_out->mkdir(dirName.c_str());
0280     p.second.write();
0281   }
0282 }
0283
0284 //----------------------------------------------------------------------------------------------------
0285
0286 DEFINE_FWK_MODULE(CTPPSAcceptancePlotter);