Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 23:31:32

0001 /****************************************************************************
0002  * Authors:
0003  *   Jan Kašpar
0004  ****************************************************************************/
0005 
0006 #include "FWCore/Framework/interface/Frameworkfwd.h"
0007 #include "FWCore/Framework/interface/EventSetup.h"
0008 #include "FWCore/Framework/interface/ESHandle.h"
0009 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0010 
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/Framework/interface/MakerMacros.h"
0013 
0014 #include "DataFormats/CTPPSDetId/interface/CTPPSDetId.h"
0015 
0016 #include "CondFormats/DataRecord/interface/CTPPSInterpolatedOpticsRcd.h"
0017 #include "CondFormats/PPSObjects/interface/LHCInterpolatedOpticalFunctionsSetCollection.h"
0018 
0019 #include "TFile.h"
0020 #include "TGraph.h"
0021 
0022 #include <map>
0023 #include <string>
0024 
0025 //----------------------------------------------------------------------------------------------------
0026 
0027 class CTPPSOpticsPlotter : public edm::one::EDAnalyzer<> {
0028 public:
0029   explicit CTPPSOpticsPlotter(const edm::ParameterSet&);
0030 
0031 private:
0032   void analyze(const edm::Event&, const edm::EventSetup&) override;
0033   void endJob() override;
0034 
0035   edm::ESGetToken<LHCInterpolatedOpticalFunctionsSetCollection, CTPPSInterpolatedOpticsRcd> opticsESToken_;
0036 
0037   unsigned int rpId_45_N_, rpId_45_F_;
0038   unsigned int rpId_56_N_, rpId_56_F_;
0039 
0040   std::string outputFile_;
0041 
0042   struct RPPlots {
0043     std::unique_ptr<TGraph> g_v_x_vs_xi, g_L_x_vs_xi, g_x_D_vs_xi;
0044     std::unique_ptr<TGraph> g_v_y_vs_xi, g_L_y_vs_xi, g_y_D_vs_xi;
0045     std::unique_ptr<TGraph> h_y_vs_x_disp;
0046 
0047     RPPlots()
0048         : g_v_x_vs_xi(new TGraph),
0049           g_L_x_vs_xi(new TGraph),
0050           g_x_D_vs_xi(new TGraph),
0051           g_v_y_vs_xi(new TGraph),
0052           g_L_y_vs_xi(new TGraph),
0053           g_y_D_vs_xi(new TGraph),
0054           h_y_vs_x_disp(new TGraph) {}
0055 
0056     void write() const {
0057       g_v_x_vs_xi->SetTitle(";xi;v_{x}");
0058       g_v_x_vs_xi->Write("g_v_x_vs_xi");
0059 
0060       g_L_x_vs_xi->SetTitle(";xi;L_{x}   (cm)");
0061       g_L_x_vs_xi->Write("g_L_x_vs_xi");
0062 
0063       g_x_D_vs_xi->SetTitle(";xi;x_{D}   (cm)");
0064       g_x_D_vs_xi->Write("g_x_D_vs_xi");
0065 
0066       g_v_y_vs_xi->SetTitle(";xi;v_{y}");
0067       g_v_y_vs_xi->Write("g_v_y_vs_xi");
0068 
0069       g_L_y_vs_xi->SetTitle(";xi;L_{y}   (cm)");
0070       g_L_y_vs_xi->Write("g_L_y_vs_xi");
0071 
0072       g_y_D_vs_xi->SetTitle(";xi;y_{D}   (cm)");
0073       g_y_D_vs_xi->Write("g_y_D_vs_xi");
0074 
0075       h_y_vs_x_disp->SetTitle(";x   (cm);y   (cm)");
0076       h_y_vs_x_disp->Write("h_y_vs_x_disp");
0077     }
0078   };
0079 
0080   std::map<unsigned int, RPPlots> rp_plots_;
0081 
0082   struct ArmPlots {
0083     unsigned int id_N, id_F;
0084 
0085     std::unique_ptr<TGraph> g_de_x_vs_x_disp, g_de_y_vs_x_disp;
0086 
0087     ArmPlots() : g_de_x_vs_x_disp(new TGraph), g_de_y_vs_x_disp(new TGraph) {}
0088 
0089     void write() const {
0090       g_de_x_vs_x_disp->SetTitle(";x_N   (cm);x_F - x_N   (cm)");
0091       g_de_x_vs_x_disp->Write("g_de_x_vs_x_disp");
0092 
0093       g_de_y_vs_x_disp->SetTitle(";x_N   (cm);y_F - y_N   (cm)");
0094       g_de_y_vs_x_disp->Write("g_de_y_vs_x_disp");
0095     }
0096   };
0097 
0098   std::map<unsigned int, ArmPlots> arm_plots_;
0099 };
0100 
0101 //----------------------------------------------------------------------------------------------------
0102 
0103 using namespace std;
0104 using namespace edm;
0105 
0106 //----------------------------------------------------------------------------------------------------
0107 
0108 CTPPSOpticsPlotter::CTPPSOpticsPlotter(const edm::ParameterSet& iConfig)
0109     : opticsESToken_(esConsumes(ESInputTag("", iConfig.getParameter<std::string>("opticsLabel")))),
0110 
0111       rpId_45_N_(iConfig.getParameter<unsigned int>("rpId_45_N")),
0112       rpId_45_F_(iConfig.getParameter<unsigned int>("rpId_45_F")),
0113       rpId_56_N_(iConfig.getParameter<unsigned int>("rpId_56_N")),
0114       rpId_56_F_(iConfig.getParameter<unsigned int>("rpId_56_F")),
0115 
0116       outputFile_(iConfig.getParameter<string>("outputFile")) {
0117   arm_plots_[0].id_N = rpId_45_N_;
0118   arm_plots_[0].id_F = rpId_45_F_;
0119 
0120   arm_plots_[1].id_N = rpId_56_N_;
0121   arm_plots_[1].id_F = rpId_56_F_;
0122 }
0123 
0124 //----------------------------------------------------------------------------------------------------
0125 
0126 void CTPPSOpticsPlotter::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0127   // stop if plots already made
0128   if (!rp_plots_.empty())
0129     return;
0130 
0131   // get conditions
0132   const auto& opticalFunctions = iSetup.getData(opticsESToken_);
0133 
0134   // stop if conditions invalid
0135   if (opticalFunctions.empty())
0136     return;
0137 
0138   // make per-RP plots
0139   for (const auto& it : opticalFunctions) {
0140     CTPPSDetId rpId(it.first);
0141     unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
0142 
0143     auto& pl = rp_plots_[rpDecId];
0144 
0145     LHCInterpolatedOpticalFunctionsSet::Kinematics k_in_beam = {0., 0., 0., 0., 0.};
0146     LHCInterpolatedOpticalFunctionsSet::Kinematics k_out_beam;
0147     it.second.transport(k_in_beam, k_out_beam);
0148 
0149     const double vtx_ep = 1E-4;  // cm
0150     const double th_ep = 1E-6;   // rad
0151 
0152     for (double xi = 0.; xi < 0.30001; xi += 0.001) {
0153       LHCInterpolatedOpticalFunctionsSet::Kinematics k_in_xi = {0., 0., 0., 0., xi};
0154       LHCInterpolatedOpticalFunctionsSet::Kinematics k_out_xi;
0155       it.second.transport(k_in_xi, k_out_xi);
0156 
0157       LHCInterpolatedOpticalFunctionsSet::Kinematics k_in_xi_vtx_x = {vtx_ep, 0., 0., 0., xi};
0158       LHCInterpolatedOpticalFunctionsSet::Kinematics k_out_xi_vtx_x;
0159       it.second.transport(k_in_xi_vtx_x, k_out_xi_vtx_x);
0160 
0161       LHCInterpolatedOpticalFunctionsSet::Kinematics k_in_xi_th_x = {0., th_ep, 0., 0., xi};
0162       LHCInterpolatedOpticalFunctionsSet::Kinematics k_out_xi_th_x;
0163       it.second.transport(k_in_xi_th_x, k_out_xi_th_x);
0164 
0165       LHCInterpolatedOpticalFunctionsSet::Kinematics k_in_xi_vtx_y = {0., 0., vtx_ep, 0., xi};
0166       LHCInterpolatedOpticalFunctionsSet::Kinematics k_out_xi_vtx_y;
0167       it.second.transport(k_in_xi_vtx_y, k_out_xi_vtx_y);
0168 
0169       LHCInterpolatedOpticalFunctionsSet::Kinematics k_in_xi_th_y = {0., 0., 0., th_ep, xi};
0170       LHCInterpolatedOpticalFunctionsSet::Kinematics k_out_xi_th_y;
0171       it.second.transport(k_in_xi_th_y, k_out_xi_th_y);
0172 
0173       int idx = pl.g_v_x_vs_xi->GetN();
0174 
0175       pl.g_v_x_vs_xi->SetPoint(idx, xi, (k_out_xi_vtx_x.x - k_out_xi.x) / vtx_ep);
0176       pl.g_L_x_vs_xi->SetPoint(idx, xi, (k_out_xi_th_x.x - k_out_xi.x) / th_ep);
0177       pl.g_x_D_vs_xi->SetPoint(idx, xi, k_out_xi.x - k_out_beam.x);
0178 
0179       pl.g_v_y_vs_xi->SetPoint(idx, xi, (k_out_xi_vtx_y.y - k_out_xi.y) / vtx_ep);
0180       pl.g_L_y_vs_xi->SetPoint(idx, xi, (k_out_xi_th_y.y - k_out_xi.y) / th_ep);
0181       pl.g_y_D_vs_xi->SetPoint(idx, xi, k_out_xi.y - k_out_beam.y);
0182 
0183       pl.h_y_vs_x_disp->SetPoint(idx, k_out_xi.x - k_out_beam.x, k_out_xi.y - k_out_beam.y);
0184     }
0185   }
0186 
0187   // make per-arm plots
0188   for (const auto& ap : arm_plots_) {
0189     // find optics objects
0190     const LHCInterpolatedOpticalFunctionsSet *opt_N = nullptr, *opt_F = nullptr;
0191 
0192     for (const auto& it : opticalFunctions) {
0193       CTPPSDetId rpId(it.first);
0194       unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
0195 
0196       if (rpDecId == ap.second.id_N)
0197         opt_N = &it.second;
0198       if (rpDecId == ap.second.id_F)
0199         opt_F = &it.second;
0200     }
0201 
0202     if (!opt_N || !opt_F) {
0203       edm::LogError("CTPPSOpticsPlotter::analyze") << "Cannot find optics objects for arm " << ap.first;
0204       continue;
0205     }
0206 
0207     LHCInterpolatedOpticalFunctionsSet::Kinematics k_in_beam = {0., 0., 0., 0., 0.};
0208 
0209     LHCInterpolatedOpticalFunctionsSet::Kinematics k_out_beam_N, k_out_beam_F;
0210     opt_N->transport(k_in_beam, k_out_beam_N);
0211     opt_F->transport(k_in_beam, k_out_beam_F);
0212 
0213     for (double xi = 0.; xi < 0.30001; xi += 0.001) {
0214       LHCInterpolatedOpticalFunctionsSet::Kinematics k_in_xi = {0., 0., 0., 0., xi};
0215 
0216       LHCInterpolatedOpticalFunctionsSet::Kinematics k_out_xi_N, k_out_xi_F;
0217       opt_N->transport(k_in_xi, k_out_xi_N);
0218       opt_F->transport(k_in_xi, k_out_xi_F);
0219 
0220       int idx = ap.second.g_de_x_vs_x_disp->GetN();
0221 
0222       ap.second.g_de_x_vs_x_disp->SetPoint(
0223           idx, k_out_xi_N.x - k_out_beam_N.x, (k_out_xi_F.x - k_out_beam_F.x) - (k_out_xi_N.x - k_out_beam_N.x));
0224       ap.second.g_de_y_vs_x_disp->SetPoint(
0225           idx, k_out_xi_N.x - k_out_beam_N.x, (k_out_xi_F.y - k_out_beam_F.y) - (k_out_xi_N.y - k_out_beam_N.y));
0226     }
0227   }
0228 }
0229 
0230 //----------------------------------------------------------------------------------------------------
0231 
0232 void CTPPSOpticsPlotter::endJob() {
0233   auto f_out = std::make_unique<TFile>(outputFile_.c_str(), "recreate");
0234 
0235   for (const auto& p : rp_plots_) {
0236     gDirectory = f_out->mkdir(Form("%u", p.first));
0237     p.second.write();
0238   }
0239 
0240   for (const auto& p : arm_plots_) {
0241     gDirectory = f_out->mkdir(Form("arm %u", p.first));
0242     p.second.write();
0243   }
0244 }
0245 
0246 //----------------------------------------------------------------------------------------------------
0247 
0248 DEFINE_FWK_MODULE(CTPPSOpticsPlotter);