MuonSimHitsValidAnalyzer.cc

Go to the documentation of this file.

#include "Validation/MuonHits/plugins/MuonSimHitsValidAnalyzer.h"

#include "TFile.h"
#include "TTree.h"
#include "TBranch.h"
#include "TH1F.h"

#include <iostream>
#include <string>

using namespace edm;
using namespace std;

MuonSimHitsValidAnalyzer::MuonSimHitsValidAnalyzer(const edm::ParameterSet& iPSet)
    : fName(""),
      verbosity(0),
      label(""),
      getAllProvenances(false),
      printProvenanceInfo(false),
      nRawGenPart(0),
      count(0)

{
  fName = iPSet.getUntrackedParameter<std::string>("Name");
  verbosity = iPSet.getUntrackedParameter<int>("Verbosity");
  label = iPSet.getParameter<std::string>("Label");
  edm::ParameterSet m_Prov = iPSet.getParameter<edm::ParameterSet>("ProvenanceLookup");
  getAllProvenances = m_Prov.getUntrackedParameter<bool>("GetAllProvenances");
  printProvenanceInfo = m_Prov.getUntrackedParameter<bool>("PrintProvenanceInfo");

  nRawGenPart = 0;

  DTHitsToken_ = consumes<edm::PSimHitContainer>(iPSet.getParameter<edm::InputTag>("DTHitsSrc"));


  if (verbosity) {
    Labels l;
    labelsForToken(DTHitsToken_, l);
    edm::LogInfo("MuonSimHitsValidAnalyzer::MuonSimHitsValidAnalyzer")
        << "\n===============================\n"
        << "Initialized as EDAnalyzer with parameter values:\n"
        << "    Name      = " << fName << "\n"
        << "    Verbosity = " << verbosity << "\n"
        << "    Label     = " << label << "\n"
        << "    GetProv   = " << getAllProvenances << "\n"
        << "    PrintProv = " << printProvenanceInfo
        << "\n"
        //    << "    CSCHitsSrc=  " <<CSCHitsSrc_.label()
        //    << ":" << CSCHitsSrc_.instance() << "\n"
        << "    DTHitsSrc =  " << l.module << ":" << l.productInstance
        << "\n"
        //     << "    RPCHitsSrc=  " <<RPCHitsSrc_.label()
        //     << ":" << RPCHitsSrc_.instance() << "\n"
        << "===============================\n";
  }

  pow6 = 1000000.0;
  mom4 = 0.;
  mom1 = 0;
  costeta = 0.;
  radius = 0;
  sinteta = 0.;
  globposx = 0.;
  globposy = 0;
  nummu_DT = 0;
  nummu_CSC = 0;
  nummu_RPC = 0;
  geomToken_ = esConsumes<DTGeometry, MuonGeometryRecord>();
}

MuonSimHitsValidAnalyzer::~MuonSimHitsValidAnalyzer() {}

void MuonSimHitsValidAnalyzer::bookHistograms(DQMStore::IBooker& iBooker,
                                              edm::Run const& iRun,
                                              edm::EventSetup const& /* iSetup */) {
  meAllDTHits = nullptr;
  meMuDTHits = nullptr;
  meToF = nullptr;
  meEnergyLoss = nullptr;
  meMomentumMB1 = nullptr;
  meMomentumMB4 = nullptr;
  meLossMomIron = nullptr;
  meLocalXvsZ = nullptr;
  meLocalXvsY = nullptr;
  meGlobalXvsZ = nullptr;
  meGlobalXvsY = nullptr;
  meGlobalXvsZWm2 = nullptr;
  meGlobalXvsZWm1 = nullptr;
  meGlobalXvsZW0 = nullptr;
  meGlobalXvsZWp1 = nullptr;
  meGlobalXvsZWp2 = nullptr;
  meGlobalXvsYWm2 = nullptr;
  meGlobalXvsYWm1 = nullptr;
  meGlobalXvsYW0 = nullptr;
  meGlobalXvsYWp1 = nullptr;
  meGlobalXvsYWp2 = nullptr;
  meWheelOccup = nullptr;
  meStationOccup = nullptr;
  meSectorOccup = nullptr;
  meSuperLOccup = nullptr;
  meLayerOccup = nullptr;
  meWireOccup = nullptr;
  mePathMuon = nullptr;
  meChamberOccup = nullptr;
  meHitRadius = nullptr;
  meCosTheta = nullptr;
  meGlobalEta = nullptr;
  meGlobalPhi = nullptr;

  Char_t histo_n[100];
  Char_t histo_t[100];

  iBooker.setCurrentFolder("MuonDTHitsV/DTHitsValidationTask");

  sprintf(histo_n, "Number_of_all_DT_hits");
  sprintf(histo_t, "Number_of_all_DT_hits");
  meAllDTHits = iBooker.book1D(histo_n, histo_t, 200, 1.0, 201.0);

  sprintf(histo_n, "Number_of_muon_DT_hits");
  sprintf(histo_t, "Number_of_muon_DT_hits");
  meMuDTHits = iBooker.book1D(histo_n, histo_t, 150, 1.0, 151.0);

  sprintf(histo_n, "Tof_of_hits ");
  sprintf(histo_t, "Tof_of_hits ");
  meToF = iBooker.book1D(histo_n, histo_t, 100, -0.5, 50.);

  sprintf(histo_n, "DT_energy_loss_keV");
  sprintf(histo_t, "DT_energy_loss_keV");
  meEnergyLoss = iBooker.book1D(histo_n, histo_t, 100, 0.0, 10.0);

  sprintf(histo_n, "Momentum_at_MB1");
  sprintf(histo_t, "Momentum_at_MB1");
  meMomentumMB1 = iBooker.book1D(histo_n, histo_t, 100, 10.0, 200.0);

  sprintf(histo_n, "Momentum_at_MB4");
  sprintf(histo_t, "Momentum_at_MB4");
  meMomentumMB4 = iBooker.book1D(histo_n, histo_t, 100, 10.0, 200.0);

  sprintf(histo_n, "Loss_of_muon_Momentum_in_Iron");
  sprintf(histo_t, "Loss_of_muon_Momentum_in_Iron");
  meLossMomIron = iBooker.book1D(histo_n, histo_t, 80, 0.0, 40.0);

  sprintf(histo_n, "Local_x-coord_vs_local_z-coord_of_muon_hit");
  sprintf(histo_t, "Local_x-coord_vs_local_z-coord_of_muon_hit");
  meLocalXvsZ = iBooker.book2D(histo_n, histo_t, 100, -150., 150., 100, -0.8, 0.8);

  sprintf(histo_n, "local_x-coord_vs_local_y-coord_of_muon_hit");
  sprintf(histo_t, "local_x-coord_vs_local_y-coord_of_muon_hit");
  meLocalXvsY = iBooker.book2D(histo_n, histo_t, 100, -150., 150., 100, -150., 150.);

  sprintf(histo_n, "Global_x-coord_vs_global_z-coord_of_muon_hit");
  sprintf(histo_t, "Global_x-coord_vs_global_z-coord_of_muon_hit");
  meGlobalXvsZ = iBooker.book2D(histo_n, histo_t, 100, -800., 800., 100, -800., 800.);

  sprintf(histo_n, "Global_x-coord_vs_global_y-coord_of_muon_hit");
  sprintf(histo_t, "Global_x-coord_vs_global_y-coord_of_muon_hit");
  meGlobalXvsY = iBooker.book2D(histo_n, histo_t, 100, -800., 800., 100, -800., 800.);

  // New histos

  sprintf(histo_n, "Global_x-coord_vs_global_z-coord_of_muon_hit_w-2");
  sprintf(histo_t, "Global_x-coord_vs_global_z-coord_of_muon_hit_w-2");
  meGlobalXvsZWm2 = iBooker.book2D(histo_n, histo_t, 100, -800., 800., 100, -800., 800.);

  sprintf(histo_n, "Global_x-coord_vs_global_y-coord_of_muon_hit_w-2");
  sprintf(histo_t, "Global_x-coord_vs_global_y-coord_of_muon_hit_w-2");
  meGlobalXvsYWm2 = iBooker.book2D(histo_n, histo_t, 100, -800., 800., 100, -800., 800.);

  sprintf(histo_n, "Global_x-coord_vs_global_z-coord_of_muon_hit_w-1");
  sprintf(histo_t, "Global_x-coord_vs_global_z-coord_of_muon_hit_w-1");
  meGlobalXvsZWm1 = iBooker.book2D(histo_n, histo_t, 100, -800., 800., 100, -800., 800.);

  sprintf(histo_n, "Global_x-coord_vs_global_y-coord_of_muon_hit_w-1");
  sprintf(histo_t, "Global_x-coord_vs_global_y-coord_of_muon_hit_w-1");
  meGlobalXvsYWm1 = iBooker.book2D(histo_n, histo_t, 100, -800., 800., 100, -800., 800.);

  sprintf(histo_n, "Global_x-coord_vs_global_z-coord_of_muon_hit_w0");
  sprintf(histo_t, "Global_x-coord_vs_global_z-coord_of_muon_hit_w0");
  meGlobalXvsZW0 = iBooker.book2D(histo_n, histo_t, 100, -800., 800., 100, -800., 800.);

  sprintf(histo_n, "Global_x-coord_vs_global_y-coord_of_muon_hit_w0");
  sprintf(histo_t, "Global_x-coord_vs_global_y-coord_of_muon_hit_w0");
  meGlobalXvsYW0 = iBooker.book2D(histo_n, histo_t, 100, -800., 800., 100, -800., 800.);

  sprintf(histo_n, "Global_x-coord_vs_global_z-coord_of_muon_hit_w1");
  sprintf(histo_t, "Global_x-coord_vs_global_z-coord_of_muon_hit_w1");
  meGlobalXvsZWp1 = iBooker.book2D(histo_n, histo_t, 100, -800., 800., 100, -800., 800.);

  sprintf(histo_n, "Global_x-coord_vs_global_y-coord_of_muon_hit_w1");
  sprintf(histo_t, "Global_x-coord_vs_global_y-coord_of_muon_hit_w1");
  meGlobalXvsYWp1 = iBooker.book2D(histo_n, histo_t, 100, -800., 800., 100, -800., 800.);

  sprintf(histo_n, "Global_x-coord_vs_global_z-coord_of_muon_hit_w2");
  sprintf(histo_t, "Global_x-coord_vs_global_z-coord_of_muon_hit_w2");
  meGlobalXvsZWp2 = iBooker.book2D(histo_n, histo_t, 100, -800., 800., 100, -800., 800.);

  sprintf(histo_n, "Global_x-coord_vs_global_y-coord_of_muon_hit_w2");
  sprintf(histo_t, "Global_x-coord_vs_global_y-coord_of_muon_hit_w2");
  meGlobalXvsYWp2 = iBooker.book2D(histo_n, histo_t, 100, -800., 800., 100, -800., 800.);

  //

  sprintf(histo_n, "Wheel_occupancy");
  sprintf(histo_t, "Wheel_occupancy");
  meWheelOccup = iBooker.book1D(histo_n, histo_t, 10, -5.0, 5.0);

  sprintf(histo_n, "Station_occupancy");
  sprintf(histo_t, "Station_occupancy");
  meStationOccup = iBooker.book1D(histo_n, histo_t, 6, 0., 6.0);

  sprintf(histo_n, "Sector_occupancy");
  sprintf(histo_t, "Sector_occupancy");
  meSectorOccup = iBooker.book1D(histo_n, histo_t, 20, 0., 20.);

  sprintf(histo_n, "SuperLayer_occupancy");
  sprintf(histo_t, "SuperLayer_occupancy");
  meSuperLOccup = iBooker.book1D(histo_n, histo_t, 5, 0., 5.);

  sprintf(histo_n, "Layer_occupancy");
  sprintf(histo_t, "Layer_occupancy");
  meLayerOccup = iBooker.book1D(histo_n, histo_t, 6, 0., 6.);

  sprintf(histo_n, "Wire_occupancy");
  sprintf(histo_t, "Wire_occupancy");
  meWireOccup = iBooker.book1D(histo_n, histo_t, 100, 0., 100.);

  sprintf(histo_n, "path_followed_by_muon");
  sprintf(histo_t, "path_followed_by_muon");
  mePathMuon = iBooker.book1D(histo_n, histo_t, 160, 0., 160.);

  sprintf(histo_n, "chamber_occupancy");
  sprintf(histo_t, "chamber_occupancy");
  meChamberOccup = iBooker.book1D(histo_n, histo_t, 251, 0., 251.);

  sprintf(histo_n, "radius_of_hit");
  sprintf(histo_t, "radius_of_hit");
  meHitRadius = iBooker.book1D(histo_n, histo_t, 100, 0., 1200.);

  sprintf(histo_n, "costheta_of_hit");
  sprintf(histo_t, "costheta_of_hit");
  meCosTheta = iBooker.book1D(histo_n, histo_t, 100, -1., 1.);

  sprintf(histo_n, "global_eta_of_hit");
  sprintf(histo_t, "global_eta_of_hit");
  meGlobalEta = iBooker.book1D(histo_n, histo_t, 60, -2.7, 2.7);

  sprintf(histo_n, "global_phi_of_hit");
  sprintf(histo_t, "global_phi_of_hit");
  meGlobalPhi = iBooker.book1D(histo_n, histo_t, 60, -3.14, 3.14);
}

void MuonSimHitsValidAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
  ++count;

  edm::RunNumber_t nrun = iEvent.id().run();
  edm::EventNumber_t nevt = iEvent.id().event();

  if (verbosity > 0) {
    edm::LogInfo("MuonSimHitsValidAnalyzer::analyze") << "Processing run " << nrun << ", event " << nevt;
  }

  if (getAllProvenances) {
    std::vector<const edm::StableProvenance*> AllProv;
    iEvent.getAllStableProvenance(AllProv);

    if (verbosity > 0)
      edm::LogInfo("MuonSimHitsValidAnalyzer::analyze") << "Number of Provenances = " << AllProv.size();

    if (printProvenanceInfo && (verbosity > 0)) {
      TString eventout("\nProvenance info:\n");

      for (unsigned int i = 0; i < AllProv.size(); ++i) {
        eventout += "\n       ******************************";
        eventout += "\n       Module       : ";
        eventout += AllProv[i]->moduleLabel();
        eventout += "\n       ProductID    : ";
        eventout += AllProv[i]->productID().id();
        eventout += "\n       ClassName    : ";
        eventout += AllProv[i]->className();
        eventout += "\n       InstanceName : ";
        eventout += AllProv[i]->productInstanceName();
        eventout += "\n       BranchName   : ";
        eventout += AllProv[i]->branchName();
      }
      eventout += "       ******************************\n";
      edm::LogInfo("MuonSimHitsValidAnalyzer::analyze") << eventout << "\n";
    }
  }

  fillDT(iEvent, iSetup);

  if (verbosity > 0)
    edm::LogInfo("MuonSimHitsValidAnalyzer::analyze") << "Done gathering data from event.";

  return;
}

void MuonSimHitsValidAnalyzer::fillDT(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
  TString eventout;
  if (verbosity > 0)
    eventout = "\nGathering DT info:";

  edm::PSimHitContainer::const_iterator itHit;

  edm::ESHandle<DTGeometry> hGeom = iSetup.getHandle(geomToken_);
  if (!hGeom.isValid()) {
    edm::LogWarning("MuonSimHitsValidAnalyzer::fillDT")
        << "Unable to find MuonGeometryRecord for the DTGeometry in event!";
    return;
  }
  const DTGeometry& theDTMuon(*hGeom);

  edm::Handle<edm::PSimHitContainer> MuonDTContainer;
  iEvent.getByToken(DTHitsToken_, MuonDTContainer);
  if (!MuonDTContainer.isValid()) {
    edm::LogWarning("MuonSimHitsValidAnalyzer::fillDT") << "Unable to find MuonDTHits in event!";
    return;
  }

  touch1 = 0;
  touch4 = 0;
  nummu_DT = 0;

  meAllDTHits->Fill(MuonDTContainer->size());

  int i = 0, j = 0;
  for (itHit = MuonDTContainer->begin(); itHit != MuonDTContainer->end(); ++itHit) {
    ++i;

    DetId theDetUnitId(itHit->detUnitId());
    int detector = theDetUnitId.det();
    int subdetector = theDetUnitId.subdetId();

    if ((detector == dMuon) && (subdetector == sdMuonDT)) {
      const GeomDetUnit* theDet = theDTMuon.idToDetUnit(theDetUnitId);

      if (!theDet) {
        edm::LogWarning("MuonSimHitsValidAnalyzer::fillDT") << "Unable to get GeomDetUnit from theDTMuon for hit " << i;
        continue;
      }

      ++j;

      const BoundPlane& bsurf = theDet->surface();


      if (abs(itHit->particleType()) == 13) {
        nummu_DT++;
        meToF->Fill(itHit->tof());
        meEnergyLoss->Fill(itHit->energyLoss() * pow6);

        iden = itHit->detUnitId();

        wheel = ((iden >> 15) & 0x7) - 3;
        station = ((iden >> 22) & 0x7);
        sector = ((iden >> 18) & 0xf);
        superlayer = ((iden >> 13) & 0x3);
        layer = ((iden >> 10) & 0x7);
        wire = ((iden >> 3) & 0x7f);

        meWheelOccup->Fill((float)wheel);
        meStationOccup->Fill((float)station);
        meSectorOccup->Fill((float)sector);
        meSuperLOccup->Fill((float)superlayer);
        meLayerOccup->Fill((float)layer);
        meWireOccup->Fill((float)wire);

        // Define a quantity to take into account station, splayer and layer being hit.
        path = (station - 1) * 40 + superlayer * 10 + layer;
        mePathMuon->Fill((float)path);

        // Define a quantity to take into chamber being hit.
        pathchamber = (wheel + 2) * 50 + (station - 1) * 12 + sector;
        meChamberOccup->Fill((float)pathchamber);

        if (station == 1) {
          if (touch1 == 0) {
            mom1 = itHit->pabs();
            meMomentumMB1->Fill(mom1);
            touch1 = 1;
          }
        }

        if (station == 4) {
          if (touch4 == 0) {
            mom4 = itHit->pabs();
            touch4 = 1;
            meMomentumMB4->Fill(mom4);
            if (touch1 == 1) {
              meLossMomIron->Fill(mom1 - mom4);
            }
          }
        }

        meLocalXvsZ->Fill(itHit->localPosition().x(), itHit->localPosition().z());

        meLocalXvsY->Fill(itHit->localPosition().x(), itHit->localPosition().y());


        globposz = bsurf.toGlobal(itHit->localPosition()).z();
        globposeta = bsurf.toGlobal(itHit->localPosition()).eta();
        globposphi = bsurf.toGlobal(itHit->localPosition()).phi();

        radius = globposz * (1. + exp(-2. * globposeta)) / (1. - exp(-2. * globposeta));

        costeta = (1. - exp(-2. * globposeta)) / (1. + exp(-2. * globposeta));
        sinteta = 2. * exp(-globposeta) / (1. + exp(-2. * globposeta));

        globposx = radius * sinteta * cos(globposphi);
        globposy = radius * sinteta * sin(globposphi);

        meGlobalXvsZ->Fill(globposz, globposx);
        meGlobalXvsY->Fill(globposx, globposy);

        //  New Histos
        if (wheel == -2) {
          meGlobalXvsZWm2->Fill(globposz, globposx);
          meGlobalXvsYWm2->Fill(globposx, globposy);
        }
        if (wheel == -1) {
          meGlobalXvsZWm1->Fill(globposz, globposx);
          meGlobalXvsYWm1->Fill(globposx, globposy);
        }
        if (wheel == 0) {
          meGlobalXvsZW0->Fill(globposz, globposx);
          meGlobalXvsYW0->Fill(globposx, globposy);
        }
        if (wheel == 1) {
          meGlobalXvsZWp1->Fill(globposz, globposx);
          meGlobalXvsYWp1->Fill(globposx, globposy);
        }
        if (wheel == 2) {
          meGlobalXvsZWp2->Fill(globposz, globposx);
          meGlobalXvsYWp2->Fill(globposx, globposy);
        }
        //
        meHitRadius->Fill(radius);
        meCosTheta->Fill(costeta);
        meGlobalEta->Fill(globposeta);
        meGlobalPhi->Fill(globposphi);
      }
    } else {
      edm::LogWarning("MuonSimHitsValidAnalyzer::fillDT")
          << "MuonDT PSimHit " << i << " is expected to be (det,subdet) = (" << dMuon << "," << sdMuonDT
          << "); value returned is: (" << detector << "," << subdetector << ")";
      continue;
    }
  }

  if (verbosity > 1) {
    eventout += "\n          Number of DT muon Hits collected:......... ";
    eventout += j;
  }
  meMuDTHits->Fill((float)nummu_DT);

  if (verbosity > 0)
    edm::LogInfo("MuonSimHitsValidAnalyzer::fillDT") << eventout << "\n";
  return;
}