Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-10-27 04:18:22

0001 /****************************************************************************
0002  *
0003  * This is a part of CTPPS validation software
0004  * Authors:
0005  *   Jan Kašpar
0006  *   Laurent Forthomme
0007  *
0008  ****************************************************************************/
0009 
0010 #include "FWCore/Framework/interface/Frameworkfwd.h"
0011 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0012 
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/Framework/interface/MakerMacros.h"
0015 
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017 
0018 #include "DataFormats/CTPPSDetId/interface/CTPPSDetId.h"
0019 #include "DataFormats/CTPPSReco/interface/CTPPSLocalTrackLite.h"
0020 #include "DataFormats/CTPPSReco/interface/CTPPSLocalTrackLiteFwd.h"
0021 
0022 #include "TFile.h"
0023 #include "TH2D.h"
0024 #include "TH1D.h"
0025 #include "TProfile.h"
0026 #include "TProfile2D.h"
0027 #include "TGraph.h"
0028 
0029 #include <map>
0030 #include <memory>
0031 
0032 //----------------------------------------------------------------------------------------------------
0033 
0034 class CTPPSTrackDistributionPlotter : public edm::one::EDAnalyzer<> {
0035 public:
0036   explicit CTPPSTrackDistributionPlotter(const edm::ParameterSet&);
0037 
0038   ~CTPPSTrackDistributionPlotter() override {}
0039 
0040 private:
0041   void analyze(const edm::Event&, const edm::EventSetup&) override;
0042   void endJob() override;
0043 
0044   edm::EDGetTokenT<CTPPSLocalTrackLiteCollection> tracksToken_;
0045 
0046   double x_pitch_pixels_;
0047 
0048   std::string outputFile_;
0049 
0050   unsigned int events_total_;
0051   std::map<unsigned int, unsigned int> events_per_arm_;
0052 
0053   struct RPPlots {
0054     bool initialized;
0055 
0056     std::unique_ptr<TH2D> h2_y_vs_x;
0057     std::unique_ptr<TProfile> p_y_vs_x;
0058     std::unique_ptr<TH1D> h_x;
0059     std::unique_ptr<TH1D> h_y;
0060     std::unique_ptr<TH1D> h_time;
0061 
0062     std::unique_ptr<TH2D> h2_de_x_vs_x;  // "delta" refers to distance between two tracks in the same RP
0063     std::unique_ptr<TH2D> h2_de_x_vs_y;
0064     std::unique_ptr<TH2D> h2_de_y_vs_x;
0065     std::unique_ptr<TH2D> h2_de_y_vs_y;
0066 
0067     RPPlots() : initialized(false) {}
0068 
0069     void init(bool pixel, double pitch) {
0070       const double bin_size_x = (pixel) ? pitch * cos(18.4 / 180. * M_PI) : 100E-3;
0071 
0072       h2_y_vs_x = std::make_unique<TH2D>("", "", 300, -10., +70., 600, -30., +30.);
0073       p_y_vs_x = std::make_unique<TProfile>("", "", 300, -10., +70.);
0074 
0075       int n_mi = ceil(10. / bin_size_x);
0076       int n_pl = ceil(70. / bin_size_x);
0077 
0078       h_x = std::make_unique<TH1D>("", "", n_mi + n_pl, -n_mi * bin_size_x, +n_pl * bin_size_x);
0079 
0080       h_y = std::make_unique<TH1D>("", "", 300, -15., +15.);
0081 
0082       h_time = std::make_unique<TH1D>("", ";time", 500, -50., +50.);
0083 
0084       h2_de_x_vs_x =
0085           std::make_unique<TH2D>("h2_de_x_vs_x", "h2_de_x_vs_x;x;distance in x axis", 300, -30., +30., 300, -3., +3.);
0086       h2_de_x_vs_y =
0087           std::make_unique<TH2D>("h2_de_x_vs_y", "h2_de_x_vs_y;y;distance in x axis", 300, -30., +30., 300, -3., +3.);
0088       h2_de_y_vs_x =
0089           std::make_unique<TH2D>("h2_de_y_vs_x", "h2_de_y_vs_x;x;distance in y axis", 300, -30., +30., 300, -3., +3.);
0090       h2_de_y_vs_y =
0091           std::make_unique<TH2D>("h2_de_y_vs_y", "h2_de_y_vs_y;y;distance in y axis", 300, -30., +30., 300, -3., +3.);
0092 
0093       initialized = true;
0094     }
0095 
0096     void fillOneTrack(double x, double y, double time) {
0097       h2_y_vs_x->Fill(x, y);
0098       p_y_vs_x->Fill(x, y);
0099       h_x->Fill(x);
0100       h_y->Fill(y);
0101       h_time->Fill(time);
0102     }
0103 
0104     void write() const {
0105       h2_y_vs_x->Write("h2_y_vs_x");
0106       p_y_vs_x->Write("p_y_vs_x");
0107       h_x->Write("h_x");
0108       h_y->Write("h_y");
0109       h_time->Write("h_time");
0110 
0111       h2_de_x_vs_x->Write("h2_de_x_vs_x");
0112       h2_de_x_vs_y->Write("h2_de_x_vs_y");
0113       h2_de_y_vs_x->Write("h2_de_y_vs_x");
0114       h2_de_y_vs_y->Write("h2_de_y_vs_y");
0115     }
0116   };
0117 
0118   std::map<unsigned int, RPPlots> rpPlots;
0119 
0120   struct ArmPlots {
0121     unsigned int rpId_N, rpId_F;
0122 
0123     std::unique_ptr<TH1D> h_de_x, h_de_y;
0124     std::unique_ptr<TProfile> p_de_x_vs_x, p_de_y_vs_x;
0125     std::unique_ptr<TProfile> p_de_x_vs_y, p_de_y_vs_y;
0126     std::unique_ptr<TProfile2D> p2_de_x_vs_x_y, p2_de_y_vs_x_y;
0127     std::unique_ptr<TH2D> h2_de_x_vs_x, h2_de_y_vs_x;
0128     std::unique_ptr<TH2D> h2_de_y_vs_de_x;
0129 
0130     struct Stat {
0131       unsigned int sN = 0, s1 = 0;
0132     };
0133 
0134     std::map<unsigned int, std::map<unsigned int, Stat>> m_stat;
0135 
0136     ArmPlots()
0137         : h_de_x(new TH1D("", ";x^{F} - x^{N}", 100, -1., +1.)),
0138           h_de_y(new TH1D("", ";y^{F} - y^{N}", 100, -1., +1.)),
0139 
0140           p_de_x_vs_x(new TProfile("", ";x^{N};x^{F} - x^{N}", 80, 0., 40., "s")),
0141           p_de_y_vs_x(new TProfile("", ";x^{N};y^{F} - y^{N}", 80, 0., 40., "s")),
0142 
0143           p_de_x_vs_y(new TProfile("", ";y^{N};x^{F} - x^{N}", 80, -20., +20., "s")),
0144           p_de_y_vs_y(new TProfile("", ";y^{N};y^{F} - y^{N}", 80, -20., +20., "s")),
0145 
0146           p2_de_x_vs_x_y(new TProfile2D("", ";x;y;x^{F} - x^{N}", 80, 0., 40., 80, -20., +20., "s")),
0147           p2_de_y_vs_x_y(new TProfile2D("", ";x;y;y^{F} - y^{N}", 80, 0., 40., 80, -20., +20., "s")),
0148 
0149           h2_de_x_vs_x(new TH2D("", ";x^{N};x^{F} - x^{N}", 80, 0., 40., 100, -1., +1.)),
0150           h2_de_y_vs_x(new TH2D("", ";x^{N};y^{F} - y^{N}", 80, 0., 40., 100, -1., +1.)),
0151           h2_de_y_vs_de_x(new TH2D("", ";x^{F} - x^{N};y^{F} - y^{N}", 100, -1., +1., 100, -1., +1.)) {}
0152 
0153     void fill(double x_N, double y_N, double x_F, double y_F) {
0154       h_de_x->Fill(x_F - x_N);
0155       h_de_y->Fill(y_F - y_N);
0156 
0157       p_de_x_vs_x->Fill(x_N, x_F - x_N);
0158       p_de_y_vs_x->Fill(x_N, y_F - y_N);
0159 
0160       p_de_x_vs_y->Fill(y_N, x_F - x_N);
0161       p_de_y_vs_y->Fill(y_N, y_F - y_N);
0162 
0163       p2_de_x_vs_x_y->Fill(x_N, y_N, x_F - x_N);
0164       p2_de_y_vs_x_y->Fill(x_N, y_N, y_F - y_N);
0165 
0166       h2_de_x_vs_x->Fill(x_N, x_F - x_N);
0167       h2_de_y_vs_x->Fill(x_N, y_F - y_N);
0168 
0169       h2_de_y_vs_de_x->Fill(x_F - x_N, y_F - y_N);
0170     }
0171 
0172     void write() const {
0173       h_de_x->Write("h_de_x");
0174       h_de_y->Write("h_de_y");
0175 
0176       p_de_x_vs_x->Write("p_de_x_vs_x");
0177       buildRMSHistogram(p_de_x_vs_x)->Write("h_rms_de_x_vs_x");
0178       p_de_y_vs_x->Write("p_de_y_vs_x");
0179       buildRMSHistogram(p_de_y_vs_x)->Write("h_rms_de_y_vs_x");
0180 
0181       p_de_x_vs_y->Write("p_de_x_vs_y");
0182       buildRMSHistogram(p_de_x_vs_y)->Write("h_rms_de_x_vs_y");
0183       p_de_y_vs_y->Write("p_de_y_vs_y");
0184       buildRMSHistogram(p_de_y_vs_y)->Write("h_rms_de_y_vs_y");
0185 
0186       p2_de_x_vs_x_y->Write("p2_de_x_vs_x_y");
0187       buildRMSHistogram(p2_de_x_vs_x_y)->Write("h2_rms_de_x_vs_x_y");
0188       p2_de_y_vs_x_y->Write("p2_de_y_vs_x_y");
0189       buildRMSHistogram(p2_de_y_vs_x_y)->Write("h2_rms_de_y_vs_x_y");
0190 
0191       h2_de_x_vs_x->Write("h2_de_x_vs_x");
0192       h2_de_y_vs_x->Write("h2_de_y_vs_x");
0193 
0194       h2_de_y_vs_de_x->Write("h2_de_y_vs_de_x");
0195 
0196       for (const auto& rp : m_stat) {
0197         TGraph* g = new TGraph();
0198 
0199         char buf[100];
0200         sprintf(buf, "g_mean_track_mult_run_%u", rp.first);
0201         g->SetName(buf);
0202 
0203         for (const auto& lsp : rp.second) {
0204           const int idx = g->GetN();
0205           const double m = (lsp.second.s1 > 0) ? double(lsp.second.sN) / lsp.second.s1 : 0.;
0206           g->SetPoint(idx, lsp.first, m);
0207         }
0208 
0209         g->Write();
0210       }
0211     }
0212 
0213     static std::unique_ptr<TH1D> buildRMSHistogram(const std::unique_ptr<TProfile>& p) {
0214       std::unique_ptr<TH1D> output =
0215           std::make_unique<TH1D>("", p->GetTitle(), p->GetNbinsX(), p->GetXaxis()->GetXmin(), p->GetXaxis()->GetXmax());
0216 
0217       for (int bi = 1; bi <= output->GetNbinsX(); ++bi)
0218         output->SetBinContent(bi, p->GetBinError(bi));
0219 
0220       return output;
0221     }
0222 
0223     static std::unique_ptr<TH2D> buildRMSHistogram(const std::unique_ptr<TProfile2D>& p) {
0224       std::unique_ptr<TH2D> output = std::make_unique<TH2D>("",
0225                                                             p->GetTitle(),
0226                                                             p->GetNbinsX(),
0227                                                             p->GetXaxis()->GetXmin(),
0228                                                             p->GetXaxis()->GetXmax(),
0229                                                             p->GetNbinsY(),
0230                                                             p->GetYaxis()->GetXmin(),
0231                                                             p->GetYaxis()->GetXmax());
0232 
0233       for (int bi_x = 1; bi_x <= output->GetNbinsX(); ++bi_x)
0234         for (int bi_y = 1; bi_y <= output->GetNbinsY(); ++bi_y)
0235           output->SetBinContent(bi_x, bi_y, p->GetBinError(bi_x, bi_y));
0236 
0237       return output;
0238     }
0239   };
0240 
0241   std::map<unsigned int, ArmPlots> armPlots;
0242 };
0243 
0244 //----------------------------------------------------------------------------------------------------
0245 
0246 CTPPSTrackDistributionPlotter::CTPPSTrackDistributionPlotter(const edm::ParameterSet& iConfig)
0247     : tracksToken_(consumes<CTPPSLocalTrackLiteCollection>(iConfig.getParameter<edm::InputTag>("tagTracks"))),
0248       x_pitch_pixels_(iConfig.getUntrackedParameter<double>("x_pitch_pixels", 150E-3)),
0249       outputFile_(iConfig.getParameter<std::string>("outputFile")),
0250       events_total_(0) {
0251   armPlots[0].rpId_N = iConfig.getParameter<unsigned int>("rpId_45_N");
0252   armPlots[0].rpId_F = iConfig.getParameter<unsigned int>("rpId_45_F");
0253 
0254   armPlots[1].rpId_N = iConfig.getParameter<unsigned int>("rpId_56_N");
0255   armPlots[1].rpId_F = iConfig.getParameter<unsigned int>("rpId_56_F");
0256 }
0257 
0258 //----------------------------------------------------------------------------------------------------
0259 
0260 void CTPPSTrackDistributionPlotter::analyze(const edm::Event& iEvent, const edm::EventSetup&) {
0261   // get input
0262   edm::Handle<CTPPSLocalTrackLiteCollection> tracks;
0263   iEvent.getByToken(tracksToken_, tracks);
0264 
0265   // process tracks
0266   std::map<unsigned int, unsigned int> m_mult;
0267 
0268   for (const auto& trk : *tracks) {
0269     CTPPSDetId rpId(trk.rpId());
0270     unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
0271     bool rpPixel = (rpId.subdetId() == CTPPSDetId::sdTrackingPixel);
0272 
0273     auto& pl = rpPlots[rpDecId];
0274     if (!pl.initialized)
0275       pl.init(rpPixel, x_pitch_pixels_);
0276 
0277     pl.fillOneTrack(trk.x(), trk.y(), trk.time());
0278 
0279     m_mult[rpId.arm()]++;
0280   }
0281 
0282   for (unsigned int arm = 0; arm < 2; ++arm) {
0283     auto& st = armPlots[arm].m_stat[iEvent.id().run()][iEvent.id().luminosityBlock()];
0284     st.s1++;
0285     st.sN += m_mult[arm];
0286   }
0287 
0288   for (const auto& t1 : *tracks) {
0289     const CTPPSDetId rpId1(t1.rpId());
0290 
0291     for (const auto& t2 : *tracks) {
0292       const CTPPSDetId rpId2(t2.rpId());
0293 
0294       if (rpId1.arm() != rpId2.arm())
0295         continue;
0296 
0297       auto& ap = armPlots[rpId1.arm()];
0298 
0299       const unsigned int rpDecId1 = rpId1.arm() * 100 + rpId1.station() * 10 + rpId1.rp();
0300       const unsigned int rpDecId2 = rpId2.arm() * 100 + rpId2.station() * 10 + rpId2.rp();
0301 
0302       if (rpDecId1 == ap.rpId_N && rpDecId2 == ap.rpId_F)
0303         ap.fill(t1.x(), t1.y(), t2.x(), t2.y());
0304     }
0305   }
0306 
0307   for (unsigned int i = 0; i < tracks->size(); ++i) {
0308     for (unsigned int j = 0; j < tracks->size(); ++j) {
0309       if (i == j)
0310         continue;
0311 
0312       const auto& tr_i = tracks->at(i);
0313       const auto& tr_j = tracks->at(j);
0314 
0315       if (tr_i.rpId() != tr_j.rpId())
0316         continue;
0317 
0318       CTPPSDetId rpId(tr_i.rpId());
0319       unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
0320 
0321       auto& pl = rpPlots[rpDecId];
0322 
0323       pl.h2_de_x_vs_x->Fill(tr_i.x(), tr_j.x() - tr_i.x());
0324       pl.h2_de_x_vs_y->Fill(tr_i.y(), tr_j.x() - tr_i.x());
0325       pl.h2_de_y_vs_x->Fill(tr_i.x(), tr_j.y() - tr_i.y());
0326       pl.h2_de_y_vs_y->Fill(tr_i.y(), tr_j.y() - tr_i.y());
0327     }
0328   }
0329 
0330   // update counters
0331   events_total_++;
0332 
0333   if (m_mult[0] > 0)
0334     events_per_arm_[0]++;
0335   if (m_mult[1] > 0)
0336     events_per_arm_[1]++;
0337 }
0338 
0339 //----------------------------------------------------------------------------------------------------
0340 
0341 void CTPPSTrackDistributionPlotter::endJob() {
0342   edm::LogInfo("CTPPSTrackDistributionPlotter")
0343       << "    events processed: " << events_total_ << " (" << std::scientific << double(events_total_) << ")\n"
0344       << "    events with tracks in sector 45: " << events_per_arm_[0] << " (" << double(events_per_arm_[0]) << ")\n"
0345       << "    events with tracks in sector 56: " << events_per_arm_[1] << " (" << double(events_per_arm_[1]) << ")";
0346 
0347   auto f_out = std::make_unique<TFile>(outputFile_.c_str(), "recreate");
0348 
0349   for (const auto& it : rpPlots) {
0350     gDirectory = f_out->mkdir(Form("RP %u", it.first));
0351     it.second.write();
0352   }
0353 
0354   for (const auto& it : armPlots) {
0355     gDirectory = f_out->mkdir(Form("arm %u", it.first));
0356     it.second.write();
0357   }
0358 }
0359 
0360 //----------------------------------------------------------------------------------------------------
0361 
0362 DEFINE_FWK_MODULE(CTPPSTrackDistributionPlotter);