Line Code
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177

/*
 *  See header file for a description of this class.
 *
 *  \author S. Bolognesi - INFN Torino
 */

#include "DTVDriftAnalyzer.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "CondFormats/DTObjects/interface/DTMtime.h"
#include "CondFormats/DataRecord/interface/DTMtimeRcd.h"
#include "CondFormats/DTObjects/interface/DTRecoConditions.h"
#include "CondFormats/DataRecord/interface/DTRecoConditionsVdriftRcd.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include <iostream>
#include "TFile.h"
#include "TH1D.h"
#include "TString.h"

using namespace edm;
using namespace std;

DTVDriftAnalyzer::DTVDriftAnalyzer(const ParameterSet& pset)
    : readLegacyVDriftDB(pset.getParameter<bool>("readLegacyVDriftDB")) {
  // The root file which will contain the histos
  string rootFileName = pset.getUntrackedParameter<string>("rootFileName");
  theFile = new TFile(rootFileName.c_str(), "RECREATE");
  theFile->cd();
  mTimeMapToken_ = esConsumes<edm::Transition::BeginRun>();
  vDriftMapToken_ = esConsumes<edm::Transition::BeginRun>();
}

DTVDriftAnalyzer::~DTVDriftAnalyzer() { theFile->Close(); }

void DTVDriftAnalyzer::beginRun(const edm::Run& run, const edm::EventSetup& eventSetup) {
  if (readLegacyVDriftDB) {
    ESHandle<DTMtime> mTime = eventSetup.getHandle(mTimeMapToken_);
    mTimeMap = &*mTime;
    vDriftMap_ = nullptr;
    edm::LogVerbatim("DTVDriftAnalyzer") << "[DTVDriftAnalyzer] MTime version: " << mTime->version() << endl;
  } else {
    ESHandle<DTRecoConditions> hVdrift = eventSetup.getHandle(vDriftMapToken_);
    vDriftMap_ = &*hVdrift;
    mTimeMap = nullptr;
    // Consistency check: no parametrization is implemented for the time being
    int version = vDriftMap_->version();
    if (version != 1) {
      throw cms::Exception("Configuration") << "only version 1 is presently supported for VDriftDB";
    }
  }
}

void DTVDriftAnalyzer::endJob() {
  // Loop over DB entries

  map<uint32_t, pair<float, float>> values;

  if (readLegacyVDriftDB) {
    for (DTMtime::const_iterator mtime = mTimeMap->begin(); mtime != mTimeMap->end(); ++mtime) {
      DTWireId wireId(
          (*mtime).first.wheelId, (*mtime).first.stationId, (*mtime).first.sectorId, (*mtime).first.slId, 0, 0);
      float vdrift;
      float reso;
      DetId detId(wireId.rawId());
      // vdrift is cm/ns , resolution is cm
      mTimeMap->get(detId, vdrift, reso, DTVelocityUnits::cm_per_ns);
      values[wireId.rawId()] = make_pair(vdrift, reso);
    }
  } else {
    for (DTRecoConditions::const_iterator vd = vDriftMap_->begin(); vd != vDriftMap_->end(); ++vd) {
      DTWireId wireId(vd->first);
      float vdrift = vDriftMap_->get(wireId);
      values[vd->first] = make_pair(vdrift, 0.f);
    }
  }

  for (map<uint32_t, pair<float, float>>::const_iterator it = values.begin(); it != values.end(); ++it) {
    float vdrift = it->second.first;
    float reso = it->second.second;
    DTWireId wireId(it->first);
    // vdrift is cm/ns , resolution is cm
    edm::LogVerbatim("DTVDriftAnalyzer") << "Wire: " << wireId << endl
                                         << " vdrift (cm/ns): " << vdrift << endl
                                         << " reso (cm): " << reso << endl;

    //Define an histo for each wheel and each superlayer type
    TH1D* hVDriftHisto = theVDriftHistoMap[make_pair(wireId.wheel(), wireId.superlayer())];
    TH1D* hResoHisto = theResoHistoMap[make_pair(wireId.wheel(), wireId.superlayer())];
    if (hVDriftHisto == 0) {
      theFile->cd();
      TString name = getHistoName(wireId).c_str();
      if (wireId.superlayer() != 2) {
        hVDriftHisto = new TH1D(name + "_VDrift", "VDrift calibrated from MT per superlayer", 50, 0, 50);
        hResoHisto = new TH1D(name + "_Reso", "Reso calibrated from MT per superlayer", 50, 0, 50);
      } else {
        hVDriftHisto = new TH1D(name + "_VDrift", "VDrift calibrated from MT per superlayer", 36, 0, 36);
        hResoHisto = new TH1D(name + "_Reso", "Reso calibrated from MT per superlayer", 36, 0, 36);
      }
      theVDriftHistoMap[make_pair(wireId.wheel(), wireId.superlayer())] = hVDriftHisto;
      theResoHistoMap[make_pair(wireId.wheel(), wireId.superlayer())] = hResoHisto;
    }

    //Fill the histos and set the bin label
    int binNumber = wireId.sector() + 12 * (wireId.station() - 1);
    hVDriftHisto->SetBinContent(binNumber, vdrift);
    hResoHisto->SetBinContent(binNumber, reso);
    string labelName;
    stringstream theStream;
    if (wireId.sector() == 1)
      theStream << "MB" << wireId.station() << "_Sec" << wireId.sector();
    else
      theStream << "Sec" << wireId.sector();
    theStream >> labelName;
    hVDriftHisto->GetXaxis()->SetBinLabel(binNumber, labelName.c_str());
    hResoHisto->GetXaxis()->SetBinLabel(binNumber, labelName.c_str());

    //Define a distribution for each wheel,station and each superlayer type
    vector<int> Wh_St_SL;
    Wh_St_SL.push_back(wireId.wheel());
    Wh_St_SL.push_back(wireId.station());
    Wh_St_SL.push_back(wireId.superlayer());
    TH1D* hVDriftDistrib = theVDriftDistribMap[Wh_St_SL];
    TH1D* hResoDistrib = theResoDistribMap[Wh_St_SL];
    if (hVDriftDistrib == 0) {
      theFile->cd();
      TString name = getDistribName(wireId).c_str();
      hVDriftDistrib = new TH1D(name + "_VDrift", "VDrift calibrated from MT per superlayer", 100, 0.00530, 0.00580);
      hResoDistrib = new TH1D(name + "_Reso", "Reso calibrated from MT per superlayer", 300, 0.015, 0.045);
      theVDriftDistribMap[Wh_St_SL] = hVDriftDistrib;
      theResoDistribMap[Wh_St_SL] = hResoDistrib;
    }
    //Fill the distributions
    hVDriftDistrib->Fill(vdrift);
    hResoDistrib->Fill(reso);
  }

  //Write histos in a .root file
  theFile->cd();
  for (map<pair<int, int>, TH1D*>::const_iterator lHisto = theVDriftHistoMap.begin(); lHisto != theVDriftHistoMap.end();
       ++lHisto) {
    (*lHisto).second->GetXaxis()->LabelsOption("v");
    (*lHisto).second->Write();
  }
  for (map<pair<int, int>, TH1D*>::const_iterator lHisto = theResoHistoMap.begin(); lHisto != theResoHistoMap.end();
       ++lHisto) {
    (*lHisto).second->GetXaxis()->LabelsOption("v");
    (*lHisto).second->Write();
  }
  for (map<vector<int>, TH1D*>::const_iterator lDistrib = theVDriftDistribMap.begin();
       lDistrib != theVDriftDistribMap.end();
       ++lDistrib) {
    (*lDistrib).second->Write();
  }
  for (map<vector<int>, TH1D*>::const_iterator lDistrib = theResoDistribMap.begin();
       lDistrib != theResoDistribMap.end();
       ++lDistrib) {
    (*lDistrib).second->Write();
  }
}

string DTVDriftAnalyzer::getHistoName(const DTWireId& wId) const {
  string histoName;
  stringstream theStream;
  theStream << "Wheel" << wId.wheel() << "_SL" << wId.superlayer();
  theStream >> histoName;
  return histoName;
}

string DTVDriftAnalyzer::getDistribName(const DTWireId& wId) const {
  string histoName;
  stringstream theStream;
  theStream << "Wheel" << wId.wheel() << "_Station" << wId.station() << "_SL" << wId.superlayer();
  theStream >> histoName;
  return histoName;
}