EcalPedHists.cc

Go to the documentation of this file.

#include "CaloOnlineTools/EcalTools/plugins/EcalPedHists.h"

EcalPedHists::EcalPedHists(const edm::ParameterSet& ps)
    : runNum_(-1),
      fileName_(ps.getUntrackedParameter<std::string>("fileName", std::string("ecalPedDigiDump"))),
      barrelDigiCollection_(ps.getParameter<edm::InputTag>("EBdigiCollection")),
      endcapDigiCollection_(ps.getParameter<edm::InputTag>("EEdigiCollection")),
      headerProducer_(ps.getParameter<edm::InputTag>("headerProducer")),
      rawDataToken_(consumes<EcalRawDataCollection>(headerProducer_)),
      ebDigiToken_(consumes<EBDigiCollection>(barrelDigiCollection_)),
      eeDigiToken_(consumes<EEDigiCollection>(endcapDigiCollection_)),
      ecalMappingToken_(esConsumes<edm::Transition::BeginRun>()) {
  using namespace std;

  fedMap_ = new EcalFedMap();
  histsFilled_ = false;
  //for(int i=601; i<655; ++i)
  //{
  //  listDefaults.push_back(i);
  //}
  listFEDs_ = ps.getUntrackedParameter<vector<int> >("listFEDs");
  listEBs_ = ps.getUntrackedParameter<vector<string> >("listEBs");

  if (listFEDs_.empty()) {
    allFEDsSelected_ = false;
    //if "actual" EB id given, then convert to FEDid and put in listFEDs_
    if (!listEBs_.empty()) {
      listFEDs_.clear();
      for (vector<string>::const_iterator itr = listEBs_.begin(); itr != listEBs_.end(); ++itr) {
        listFEDs_.push_back(fedMap_->getFedFromSlice(*itr));
      }
    }
  } else if (listFEDs_[0] == -1) {
    // Apply no selection if -1 is passed in FED list
    allFEDsSelected_ = true;
    //debug
    //cout << "no selection on FEDs!" << endl;
    //inputIsOk_=false;
    //return;
    //listFEDs_ = listDefaults;
  } else {
    //in this case, listFEDs should be populated
    allFEDsSelected_ = false;
  }

  if (!allFEDsSelected_) {
    // Verify FED numbers are valid
    for (vector<int>::const_iterator intIter = listFEDs_.begin(); intIter != listFEDs_.end(); intIter++) {
      if (((*intIter) < 601) || (654 < (*intIter))) {
        cout << "[EcalPedHists] FED value: " << (*intIter) << " found in listFEDs. "
             << " Valid range is 601-654. Returning." << endl;
        inputIsOk_ = false;
        return;
      } else
        theRealFedSet_.insert(*intIter);
    }
  }

  vector<int> listDefaults = vector<int>();
  listDefaults.clear();
  for (int i = 1; i < 1701; ++i) {
    listDefaults.push_back(i);
  }
  listChannels_ = ps.getUntrackedParameter<vector<int> >("listChannels", listDefaults);
  listDefaults.clear();
  // Get samples to plot (default to 1,2,3)
  listDefaults.push_back(0);
  listDefaults.push_back(1);
  listDefaults.push_back(2);
  listSamples_ = ps.getUntrackedParameter<vector<int> >("listSamples", listDefaults);

  inputIsOk_ = true;
  vector<int>::iterator intIter;

  // Verify crystal numbers are valid
  for (intIter = listChannels_.begin(); intIter != listChannels_.end(); ++intIter) {
    //TODO: Fix crystal index checking?
    //if ( ((*intIter) < 1)||(1700 < (*intIter)) )
    //{
    //  cout << "[EcalPedHists] ic value: " << (*intIter) << " found in listChannels. "
    //        << " Valid range is 1-1700. Returning." << endl;
    //  inputIsOk_ = false;
    //  return;
    //}
  }
  // Verify sample numbers are valid
  for (intIter = listSamples_.begin(); intIter != listSamples_.end(); intIter++) {
    if (((*intIter) < 1) || (10 < (*intIter))) {
      cout << "[EcalPedHists] sample number: " << (*intIter) << " found in listSamples. "
           << " Valid range is 1-10. Returning." << endl;
      inputIsOk_ = false;
      return;
    }
  }
}

EcalPedHists::~EcalPedHists() {}

void EcalPedHists::beginRun(edm::Run const&, edm::EventSetup const& c) {
  ecalElectronicsMap_ = &c.getData(ecalMappingToken_);
}

void EcalPedHists::endRun(edm::Run const&, edm::EventSetup const& c) {}

void EcalPedHists::endJob(void) {
  using namespace std;
  if (inputIsOk_) {
    //debug
    //cout << "endJob:creating root file!" << endl;

    fileName_ += "-" + intToString(runNum_) + ".graph.root";

    TFile root_file_(fileName_.c_str(), "RECREATE");
    //Loop over FEDs first
    for (set<int>::const_iterator FEDitr = theRealFedSet_.begin(); FEDitr != theRealFedSet_.end(); ++FEDitr) {
      if (!histsFilled_)
        break;
      string dir = fedMap_->getSliceFromFed(*FEDitr);
      TDirectory* FEDdir = gDirectory->mkdir(dir.c_str());
      FEDdir->cd();
      //root_file_.mkdir(dir.c_str());
      //root_file_.cd(dir.c_str());
      map<string, TH1F*> mapHistos = FEDsAndHistMaps_[*FEDitr];

      //Loop over channels; write histos and directory structure
      for (vector<int>::const_iterator itr = listChannels_.begin(); itr != listChannels_.end(); itr++) {
        //debug
        //cout << "loop over channels" << endl;

        TH1F* hist = nullptr;
        string chnl = intToString(*itr);
        string name1 = "Cry";
        name1.append(chnl + "Gain1");
        string name2 = "Cry";
        name2.append(chnl + "Gain6");
        string name3 = "Cry";
        name3.append(chnl + "Gain12");
        hist = mapHistos[name1];
        // This is a sanity check only
        if (hist != nullptr) {
          string cryDirName = "Cry_" + chnl;
          TDirectory* cryDir = FEDdir->mkdir(cryDirName.c_str());
          cryDir->cd();
          hist->SetDirectory(cryDir);
          hist->Write();
          hist = mapHistos[name2];
          hist->SetDirectory(cryDir);
          hist->Write();
          hist = mapHistos[name3];
          hist->SetDirectory(cryDir);
          hist->Write();
          //root_file_.cd(dir.c_str());
          root_file_.cd();
        } else {
          cerr << "EcalPedHists: Error: This shouldn't happen!" << endl;
        }
      }
      root_file_.cd();
    }
    root_file_.Close();
  }
}

void EcalPedHists::analyze(const edm::Event& e, const edm::EventSetup& c) {
  using namespace std;
  using namespace edm;

  if (!inputIsOk_)
    return;

  // loop over the headers, this is to detect missing FEDs if all are selected
  if (allFEDsSelected_) {
    edm::Handle<EcalRawDataCollection> DCCHeaders;
    e.getByToken(rawDataToken_, DCCHeaders);
    if (!DCCHeaders.isValid()) {
      edm::LogError("EcalPedHists") << "Error! can't get the product " << headerProducer_;
      return;
    }

    for (EcalRawDataCollection::const_iterator headerItr = DCCHeaders->begin(); headerItr != DCCHeaders->end();
         ++headerItr) {
      int FEDid = 600 + headerItr->id();
      theRealFedSet_.insert(FEDid);
    }
  }

  // loop over fed list and make sure that there are histo maps
  for (set<int>::const_iterator fedItr = theRealFedSet_.begin(); fedItr != theRealFedSet_.end(); ++fedItr) {
    if (FEDsAndHistMaps_.find(*fedItr) == FEDsAndHistMaps_.end())
      initHists(*fedItr);
  }

  //debug
  //cout << "analyze...input is ok? " << inputIsOk_ << endl;

  bool barrelDigisFound = true;
  bool endcapDigisFound = true;
  // get the barrel digis
  // (one digi for each crystal)
  // TODO; SIC: fix this behavior
  Handle<EBDigiCollection> barrelDigis;
  e.getByToken(ebDigiToken_, barrelDigis);
  if (!barrelDigis.isValid()) {
    edm::LogError("EcalPedOffset") << "Error! can't get the product " << barrelDigiCollection_;
    barrelDigisFound = false;
  }
  // get the endcap digis
  // (one digi for each crystal)
  // TODO; SIC: fix this behavior
  Handle<EEDigiCollection> endcapDigis;
  e.getByToken(eeDigiToken_, endcapDigis);
  if (!endcapDigis.isValid()) {
    edm::LogError("EcalPedOffset") << "Error! can't get the product " << endcapDigiCollection_;
    endcapDigisFound = false;
  }

  if (barrelDigisFound)
    readEBdigis(barrelDigis);
  if (endcapDigisFound)
    readEEdigis(endcapDigis);
  if (!barrelDigisFound && !endcapDigisFound)
    edm::LogError("EcalPedOffset") << "No digis found in the event!";

  if (runNum_ == -1)
    runNum_ = e.id().run();
}

// insert the 3-entry hist map into the map keyed by FED number
void EcalPedHists::initHists(int FED) {
  using namespace std;
  //using namespace edm;

  std::map<string, TH1F*> histMap;
  //debug
  //cout << "Initializing map for FED:" << *FEDitr << endl;
  for (vector<int>::const_iterator intIter = listChannels_.begin(); intIter != listChannels_.end(); ++intIter) {
    //Put 3 histos (1 per gain) for the channel into the map
    string FEDid = intToString(FED);
    string chnl = intToString(*intIter);
    string title1 = "Gain1 ADC Counts for channel ";
    title1.append(chnl);
    string name1 = "Cry";
    name1.append(chnl + "Gain1");
    string title2 = "Gain6 ADC Counts for channel ";
    title2.append(chnl);
    string name2 = "Cry";
    name2.append(chnl + "Gain6");
    string title3 = "Gain12 ADC Counts for channel ";
    title3.append(chnl);
    string name3 = "Cry";
    name3.append(chnl + "Gain12");
    histMap.insert(make_pair(name1, new TH1F(name1.c_str(), title1.c_str(), 75, 175.0, 250.0)));
    histMap[name1]->SetDirectory(nullptr);
    histMap.insert(make_pair(name2, new TH1F(name2.c_str(), title2.c_str(), 75, 175.0, 250.0)));
    histMap[name2]->SetDirectory(nullptr);
    histMap.insert(make_pair(name3, new TH1F(name3.c_str(), title3.c_str(), 75, 175.0, 250.0)));
    histMap[name3]->SetDirectory(nullptr);
  }
  FEDsAndHistMaps_.insert(make_pair(FED, histMap));
}

void EcalPedHists::readEBdigis(edm::Handle<EBDigiCollection> digis) {
  using namespace std;
  using namespace edm;
  //debug
  //cout << "readEBdigis" << endl;

  // Loop over digis
  for (EBDigiCollection::const_iterator digiItr = digis->begin(); digiItr != digis->end(); ++digiItr) {
    EBDetId detId = EBDetId(digiItr->id());
    EcalElectronicsId elecId = ecalElectronicsMap_->getElectronicsId(detId);
    int FEDid = 600 + elecId.dccId();
    int crystalId = detId.ic();

    //debug
    //cout << "FEDid:" << FEDid << " cryId:" << crystalId << endl;
    //cout << "FEDid:" << FEDid << endl;
    //Select desired supermodules only
    set<int>::const_iterator fedIter = find(theRealFedSet_.begin(), theRealFedSet_.end(), FEDid);
    if (fedIter == theRealFedSet_.end())
      continue;

    // Select desired channels only
    vector<int>::iterator icIter;
    icIter = find(listChannels_.begin(), listChannels_.end(), crystalId);
    if (icIter == listChannels_.end())
      continue;

    // Get the adc counts from the selected samples and fill the corresponding histogram
    // Must subtract 1 from user-given sample list (e.g., user's sample 1 -> sample 0)
    for (vector<int>::iterator itr = listSamples_.begin(); itr != listSamples_.end(); itr++) {
      histsFilled_ = true;
      map<string, TH1F*> mapHistos = FEDsAndHistMaps_[FEDid];
      string chnl = intToString(crystalId);
      string name1 = "Cry";
      name1.append(chnl + "Gain1");
      string name2 = "Cry";
      name2.append(chnl + "Gain6");
      string name3 = "Cry";
      name3.append(chnl + "Gain12");
      TH1F* hist = nullptr;
      if (((EBDataFrame)(*digiItr)).sample(*itr - 1).gainId() == 3)
        hist = mapHistos[name1];
      if (((EBDataFrame)(*digiItr)).sample(*itr - 1).gainId() == 2)
        hist = mapHistos[name2];
      if (((EBDataFrame)(*digiItr)).sample(*itr - 1).gainId() == 1)
        hist = mapHistos[name3];
      if (hist != nullptr)
        hist->Fill(((EBDataFrame)(*digiItr)).sample(*itr - 1).adc());
      else
        cerr << "EcalPedHistDumper: Error: This shouldn't happen!" << endl;
    }
  }
}

void EcalPedHists::readEEdigis(edm::Handle<EEDigiCollection> digis) {
  using namespace std;
  using namespace edm;
  //debug
  //cout << "readEEdigis" << endl;

  // Loop over digis
  for (EEDigiCollection::const_iterator digiItr = digis->begin(); digiItr != digis->end(); ++digiItr) {
    EEDetId detId = EEDetId(digiItr->id());
    EcalElectronicsId elecId = ecalElectronicsMap_->getElectronicsId(detId);
    int FEDid = 600 + elecId.dccId();
    int crystalId = 10000 * FEDid + 100 * elecId.towerId() + 5 * (elecId.stripId() - 1) + elecId.xtalId();

    //Select desired FEDs only
    set<int>::const_iterator fedIter = find(theRealFedSet_.begin(), theRealFedSet_.end(), FEDid);
    if (fedIter == theRealFedSet_.end())
      continue;

    // Select desired channels only
    vector<int>::iterator icIter;
    icIter = find(listChannels_.begin(), listChannels_.end(), crystalId);
    if (icIter == listChannels_.end())
      continue;

    // Get the adc counts from the selected samples and fill the corresponding histogram
    // Must subtract 1 from user-given sample list (e.g., user's sample 1 -> sample 0)
    for (vector<int>::iterator itr = listSamples_.begin(); itr != listSamples_.end(); itr++) {
      histsFilled_ = true;
      map<string, TH1F*> mapHistos = FEDsAndHistMaps_[FEDid];
      string chnl = intToString(crystalId);
      string name1 = "Cry";
      name1.append(chnl + "Gain1");
      string name2 = "Cry";
      name2.append(chnl + "Gain6");
      string name3 = "Cry";
      name3.append(chnl + "Gain12");
      TH1F* hist = nullptr;
      if (((EBDataFrame)(*digiItr)).sample(*itr - 1).gainId() == 3)
        hist = mapHistos[name1];
      if (((EBDataFrame)(*digiItr)).sample(*itr - 1).gainId() == 2)
        hist = mapHistos[name2];
      if (((EBDataFrame)(*digiItr)).sample(*itr - 1).gainId() == 1)
        hist = mapHistos[name3];
      if (hist != nullptr)
        hist->Fill(((EBDataFrame)(*digiItr)).sample(*itr - 1).adc());
      else
        cerr << "EcalPedHistDumper: Error: This shouldn't happen!" << endl;
    }
  }
}

std::string EcalPedHists::intToString(int num) {
  using namespace std;
  //
  // outputs the number into the string stream and then flushes
  // the buffer (makes sure the output is put into the stream)
  //
  ostringstream myStream;  //creates an ostringstream object
  myStream << num << flush;
  return (myStream.str());  //returns the string form of the stringstream object
}