EcalMipGraphs.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:   EcalMipGraphs
// Class:     EcalMipGraphs
//
//
// Original Author:  Seth COOPER
//         Created:  Th Nov 22 5:46:22 CEST 2007
//
//
 
#include "CaloOnlineTools/EcalTools/plugins/EcalMipGraphs.h"
#include "RecoCaloTools/Navigation/interface/CaloNavigator.h"
#include "TCanvas.h"
#include <utility>
#include <string>
#include <vector>
 
using namespace edm;
using namespace std;
 
//
// constants, enums and typedefs
//
 
//
// static data member definitions
//
float EcalMipGraphs::gainRatio[3] = {1., 2., 12.};
edm::Service<TFileService> EcalMipGraphs::fileService;
 
//
// constructors and destructor
//
EcalMipGraphs::EcalMipGraphs(const edm::ParameterSet& iConfig)
    : EBRecHitCollection_(iConfig.getParameter<edm::InputTag>("EcalRecHitCollectionEB")),
      EERecHitCollection_(iConfig.getParameter<edm::InputTag>("EcalRecHitCollectionEE")),
      EBDigis_(iConfig.getParameter<edm::InputTag>("EBDigiCollection")),
      EEDigis_(iConfig.getParameter<edm::InputTag>("EEDigiCollection")),
      headerProducer_(iConfig.getParameter<edm::InputTag>("headerProducer")),
      runNum_(-1),
      side_(iConfig.getUntrackedParameter<int>("side", 3)),
      threshold_(iConfig.getUntrackedParameter<double>("amplitudeThreshold", 12.0)),
      minTimingAmp_(iConfig.getUntrackedParameter<double>("minimumTimingAmplitude", 0.100)) {
  vector<int> listDefaults;
  listDefaults.push_back(-1);
 
  maskedChannels_ = iConfig.getUntrackedParameter<vector<int> >("maskedChannels", listDefaults);
  maskedFEDs_ = iConfig.getUntrackedParameter<vector<int> >("maskedFEDs", listDefaults);
  seedCrys_ = iConfig.getUntrackedParameter<vector<int> >("seedCrys", vector<int>());
 
  vector<string> defaultMaskedEBs;
  defaultMaskedEBs.push_back("none");
  maskedEBs_ = iConfig.getUntrackedParameter<vector<string> >("maskedEBs", defaultMaskedEBs);
 
  fedMap_ = new EcalFedMap();
 
  string title1 = "Jitter for all FEDs";
  string name1 = "JitterAllFEDs";
  allFedsTimingHist_ = fileService->make<TH1F>(name1.c_str(), title1.c_str(), 150, -7, 7);
 
  // load up the maskedFED list with the proper FEDids
  if (maskedFEDs_[0] == -1) {
    //if "actual" EB id given, then convert to FEDid and put in listFEDs_
    if (maskedEBs_[0] != "none") {
      maskedFEDs_.clear();
      for (vector<string>::const_iterator ebItr = maskedEBs_.begin(); ebItr != maskedEBs_.end(); ++ebItr) {
        maskedFEDs_.push_back(fedMap_->getFedFromSlice(*ebItr));
      }
    }
  }
 
  for (int i = 0; i < 10; i++)
    abscissa[i] = i;
 
  naiveEvtNum_ = 0;
}
 
EcalMipGraphs::~EcalMipGraphs() {}
 
//
// member functions
//
 
// ------------ method called to for each event  ------------
void EcalMipGraphs::analyze(edm::Event const& iEvent, edm::EventSetup const& iSetup) {
  // get the headers
  // (one header for each supermodule)
  edm::Handle<EcalRawDataCollection> DCCHeaders;
  iEvent.getByLabel(headerProducer_, DCCHeaders);
 
  for (EcalRawDataCollection::const_iterator headerItr = DCCHeaders->begin(); headerItr != DCCHeaders->end();
       ++headerItr) {
    FEDsAndDCCHeaders_[headerItr->id() + 600] = *headerItr;
  }
 
  int ievt = iEvent.id().event();
  naiveEvtNum_++;
 
  if (runNum_ == -1) {
    runNum_ = iEvent.id().run();
    canvasNames_ = fileService->make<TTree>("canvasNames", "Names of written canvases");
    names = new std::vector<string>();
    canvasNames_->Branch("names", "vector<string>", &names);
  }
 
  //We only want the 3x3's for this event...
  listEBChannels.clear();
  listEEChannels.clear();
  Handle<EcalRecHitCollection> EBhits;
  Handle<EcalRecHitCollection> EEhits;
  ESHandle<CaloTopology> caloTopo;
  iSetup.get<CaloTopologyRecord>().get(caloTopo);
  iEvent.getByLabel(EBRecHitCollection_, EBhits);
  iEvent.getByLabel(EERecHitCollection_, EEhits);
  // Now, retrieve the crystal digi from the event
  iEvent.getByLabel(EBDigis_, EBdigisHandle);
  iEvent.getByLabel(EEDigis_, EEdigisHandle);
  //debug
  //LogWarning("EcalMipGraphs") << "event " << ievt << " EBhits collection size " << EBhits->size();
  //LogWarning("EcalMipGraphs") << "event " << ievt << " EEhits collection size " << EEhits->size();
 
  selectHits(EBhits, ievt, caloTopo);
  selectHits(EEhits, ievt, caloTopo);
}
 
TGraph* EcalMipGraphs::selectDigi(DetId thisDet, int ievt) {
  int emptyY[10];
  for (int i = 0; i < 10; i++)
    emptyY[i] = 0;
  TGraph* emptyGraph = fileService->make<TGraph>(10, abscissa, emptyY);
  emptyGraph->SetTitle("NOT ECAL");
 
  //If the DetId is not from Ecal, return
  if (thisDet.det() != DetId::Ecal)
    return emptyGraph;
 
  emptyGraph->SetTitle("NO DIGIS");
  //find digi we need  -- can't get find() to work; need DataFrame(DetId det) to work?
  EcalElectronicsId elecId = ecalElectronicsMap_->getElectronicsId(thisDet);
  int FEDid = 600 + elecId.dccId();
  bool isBarrel = true;
  if (FEDid < 610 || FEDid > 645)
    isBarrel = false;
  int cryIndex = isBarrel ? ((EBDetId)thisDet).hashedIndex() : getEEIndex(elecId);
  int ic = isBarrel ? ((EBDetId)thisDet).ic() : cryIndex;
 
  string sliceName = fedMap_->getSliceFromFed(FEDid);
  EcalDataFrame df;
  if (isBarrel) {
    EBDigiCollection::const_iterator digiItr = EBdigisHandle->begin();
    while (digiItr != EBdigisHandle->end() && ((*digiItr).id() != (EBDetId)thisDet)) {
      ++digiItr;
    }
    if (digiItr == EBdigisHandle->end()) {
      //LogWarning("EcalMipGraphs") << "Cannot find digi for ic:" << ic
      //  << " FED:" << FEDid << " evt:" << naiveEvtNum_;
      return emptyGraph;
    } else
      df = *digiItr;
  } else {
    EEDigiCollection::const_iterator digiItr = EEdigisHandle->begin();
    while (digiItr != EEdigisHandle->end() && ((*digiItr).id() != (EEDetId)thisDet)) {
      ++digiItr;
    }
    if (digiItr == EEdigisHandle->end()) {
      //LogWarning("EcalMipGraphs") << "Cannot find digi for ic:" << ic
      //  << " FED:" << FEDid << " evt:" << naiveEvtNum_;
      return emptyGraph;
    } else
      df = *digiItr;
  }
 
  int gainId = FEDsAndDCCHeaders_[FEDid].getMgpaGain();
  int gainHuman;
  if (gainId == 1)
    gainHuman = 12;
  else if (gainId == 2)
    gainHuman = 6;
  else if (gainId == 3)
    gainHuman = 1;
  else
    gainHuman = -1;
 
  double pedestal = 200;
 
  emptyGraph->SetTitle("FIRST TWO SAMPLES NOT GAIN12");
  if (df.sample(0).gainId() != 1 || df.sample(1).gainId() != 1)
    return emptyGraph;  //goes to the next digi
  else {
    ordinate[0] = df.sample(0).adc();
    ordinate[1] = df.sample(1).adc();
    pedestal = (double)(ordinate[0] + ordinate[1]) / (double)2;
  }
 
  for (int i = 0; i < df.size(); ++i) {
    if (df.sample(i).gainId() != 0)
      ordinate[i] = (int)(pedestal + (df.sample(i).adc() - pedestal) * gainRatio[df.sample(i).gainId() - 1]);
    else
      ordinate[i] = 49152;  //Saturation of gain1
  }
 
  TGraph* oneGraph = fileService->make<TGraph>(10, abscissa, ordinate);
  string name = "Graph_ev" + intToString(naiveEvtNum_) + "_ic" + intToString(ic) + "_FED" + intToString(FEDid);
  string gainString = (gainId == 1) ? "Free" : intToString(gainHuman);
  string title = "Event" + intToString(naiveEvtNum_) + "_lv1a" + intToString(ievt) + "_ic" + intToString(ic) + "_" +
                 sliceName + "_gain" + gainString;
 
  float energy = 0;
  map<int, float>::const_iterator itr;
  itr = crysAndAmplitudesMap_.find(cryIndex);
  if (itr != crysAndAmplitudesMap_.end()) {
    //edm::LogWarning("EcalMipGraphs")<< "itr->second(ampli)="<< itr->second;
    energy = (float)itr->second;
  }
  //else
  //edm::LogWarning("EcalMipGraphs") << "cry " << ic << "not found in ampMap";
 
  title += "_Energy" + floatToString(round(energy * 1000));
 
  oneGraph->SetTitle(title.c_str());
  oneGraph->SetName(name.c_str());
  oneGraph->GetXaxis()->SetTitle("sample");
  oneGraph->GetYaxis()->SetTitle("ADC");
  return oneGraph;
}
 
void EcalMipGraphs::selectHits(Handle<EcalRecHitCollection> hits, int ievt, ESHandle<CaloTopology> caloTopo) {
  for (EcalRecHitCollection::const_iterator hitItr = hits->begin(); hitItr != hits->end(); ++hitItr) {
    EcalRecHit hit = (*hitItr);
    DetId det = hit.id();
    EcalElectronicsId elecId = ecalElectronicsMap_->getElectronicsId(det);
    int FEDid = 600 + elecId.dccId();
    bool isBarrel = true;
    if (FEDid < 610 || FEDid > 645)
      isBarrel = false;
    int cryIndex = isBarrel ? ((EBDetId)det).hashedIndex() : ((EEDetId)det).hashedIndex();
    int ic = isBarrel ? ((EBDetId)det).ic() : getEEIndex(elecId);
 
    float ampli = hit.energy();
 
    vector<int>::iterator result;
    result = find(maskedFEDs_.begin(), maskedFEDs_.end(), FEDid);
    if (result != maskedFEDs_.end()) {
      //LogWarning("EcalMipGraphs") << "skipping uncalRecHit for FED " << FEDid << " ; amplitude " << ampli;
      continue;
    }
    result = find(maskedChannels_.begin(), maskedChannels_.end(), cryIndex);
    if (result != maskedChannels_.end()) {
      //LogWarning("EcalMipGraphs") << "skipping uncalRecHit for channel: " << cryIndex << " in fed: " << FEDid << " with amplitude " << ampli ;
      continue;
    }
    bool cryIsInList = false;
    result = find(seedCrys_.begin(), seedCrys_.end(), cryIndex);
    if (result != seedCrys_.end())
      cryIsInList = true;
 
    // Either we must have a user-requested cry (in which case there is no amplitude selection)
    // Or we pick all crys that pass the amplitude cut (in which case there is no fixed crystal selection)
    if (cryIsInList || (seedCrys_.empty() && ampli > threshold_)) {
      // We have a winner!
      crysAndAmplitudesMap_[cryIndex] = ampli;
      string name = "Event" + intToString(naiveEvtNum_) + "_ic" + intToString(ic) + "_FED" + intToString(FEDid);
      string title = "Digis";
      string seed = "ic" + intToString(ic) + "_FED" + intToString(FEDid);
      int freq = 1;
      pair<map<string, int>::iterator, bool> pair = seedFrequencyMap_.insert(make_pair(seed, freq));
      if (!pair.second) {
        ++(pair.first->second);
      }
 
      TCanvas can(name.c_str(), title.c_str(), 200, 50, 900, 900);
      can.Divide(side_, side_);
      TGraph* myGraph;
      int canvasNum = 1;
 
      CaloNavigator<DetId> cursor = CaloNavigator<DetId>(det, caloTopo->getSubdetectorTopology(det));
      //Now put each graph in one by one
      for (int j = side_ / 2; j >= -side_ / 2; --j) {
        for (int i = -side_ / 2; i <= side_ / 2; ++i) {
          cursor.home();
          cursor.offsetBy(i, j);
          can.cd(canvasNum);
          myGraph = selectDigi(*cursor, ievt);
          myGraph->Draw("A*");
          canvasNum++;
        }
      }
      can.Write();
      names->push_back(name);
    }
 
    TH1F* timingHist = FEDsAndTimingHists_[FEDid];
    if (timingHist == nullptr) {
      initHists(FEDid);
      timingHist = FEDsAndTimingHists_[FEDid];
    }
    if (ampli > minTimingAmp_) {
      timingHist->Fill(hit.time());
      allFedsTimingHist_->Fill(hit.time());
    }
  }
}
 
int EcalMipGraphs::getEEIndex(EcalElectronicsId elecId) {
  int FEDid = 600 + elecId.dccId();
  return 10000 * FEDid + 100 * elecId.towerId() + 5 * (elecId.stripId() - 1) + elecId.xtalId();
}
 
// insert the hist map into the map keyed by FED number
void EcalMipGraphs::initHists(int FED) {
  using namespace std;
 
  string title1 = "Jitter for ";
  title1.append(fedMap_->getSliceFromFed(FED));
  string name1 = "JitterFED";
  name1.append(intToString(FED));
  TH1F* timingHist = fileService->make<TH1F>(name1.c_str(), title1.c_str(), 150, -7, 7);
  FEDsAndTimingHists_[FED] = timingHist;
}
 
// ------------ method called once each job just before starting event loop  ------------
void EcalMipGraphs::beginRun(edm::Run const&, edm::EventSetup const& c) {
  edm::ESHandle<EcalElectronicsMapping> handle;
  c.get<EcalMappingRcd>().get(handle);
  ecalElectronicsMap_ = handle.product();
}
 
// ------------ method called once each job just after ending the event loop  ------------
void EcalMipGraphs::endJob() {
  canvasNames_->Fill();
 
  string frequencies = "";
  for (std::map<std::string, int>::const_iterator itr = seedFrequencyMap_.begin(); itr != seedFrequencyMap_.end();
       ++itr) {
    if (itr->second > 1) {
      frequencies += itr->first;
      frequencies += " Frequency: ";
      frequencies += intToString(itr->second);
      frequencies += "\n";
    }
  }
  LogWarning("EcalMipGraphs") << "Found seeds with frequency > 1: "
                              << "\n\n"
                              << frequencies;
 
  std::string channels;
  for (std::vector<int>::const_iterator itr = maskedChannels_.begin(); itr != maskedChannels_.end(); ++itr) {
    channels += intToString(*itr);
    channels += ",";
  }
 
  std::string feds;
  for (std::vector<int>::const_iterator itr = maskedFEDs_.begin(); itr != maskedFEDs_.end(); ++itr) {
    feds += intToString(*itr);
    feds += ",";
  }
 
  LogWarning("EcalMipGraphs") << "Masked channels are: " << channels;
  LogWarning("EcalMipGraphs") << "Masked FEDs are: " << feds << " and that is all!";
}
 
std::string EcalMipGraphs::intToString(int num) {
  using namespace std;
  ostringstream myStream;
  myStream << num << flush;
  return (myStream.str());  //returns the string form of the stringstream object
}
 
std::string EcalMipGraphs::floatToString(float num) {
  using namespace std;
  ostringstream myStream;
  myStream << num << flush;
  return (myStream.str());  //returns the string form of the stringstream object
}