EcalTools/plugins/EcalMipGraphs.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:   EcalMipGraphs
0004 // Class:     EcalMipGraphs
0005 //
0006 /**\class EcalMipGraphs EcalMipGraphs.cc
0007
0008  Description: <one line class summary>
0009
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Seth COOPER
0015 //         Created:  Th Nov 22 5:46:22 CEST 2007
0016 //
0017 //
0018
0019 #include "CaloOnlineTools/EcalTools/plugins/EcalMipGraphs.h"
0020 #include "RecoCaloTools/Navigation/interface/CaloNavigator.h"
0021 #include "TCanvas.h"
0022 #include <utility>
0023
0024 using namespace edm;
0025 using namespace std;
0026
0027 //
0028 // constants, enums and typedefs
0029 //
0030
0031 //
0032 // static data member definitions
0033 //
0034 float EcalMipGraphs::gainRatio[3] = {1., 2., 12.};
0035 edm::Service<TFileService> EcalMipGraphs::fileService;
0036
0037 //
0038 // constructors and destructor
0039 //
0040 EcalMipGraphs::EcalMipGraphs(const edm::ParameterSet& iConfig)
0041     : EBRecHitCollection_(iConfig.getParameter<edm::InputTag>("EcalRecHitCollectionEB")),
0042       EERecHitCollection_(iConfig.getParameter<edm::InputTag>("EcalRecHitCollectionEE")),
0043       EBDigis_(iConfig.getParameter<edm::InputTag>("EBDigiCollection")),
0044       EEDigis_(iConfig.getParameter<edm::InputTag>("EEDigiCollection")),
0045       headerProducer_(iConfig.getParameter<edm::InputTag>("headerProducer")),
0046       rawDataToken_(consumes<EcalRawDataCollection>(headerProducer_)),
0047       ebRecHitToken_(consumes<EcalRecHitCollection>(EBRecHitCollection_)),
0048       eeRecHitToken_(consumes<EcalRecHitCollection>(EERecHitCollection_)),
0049       ebDigiToken_(consumes<EBDigiCollection>(EBDigis_)),
0050       eeDigiToken_(consumes<EEDigiCollection>(EEDigis_)),
0051       ecalMappingToken_(esConsumes<edm::Transition::BeginRun>()),
0052       topologyToken_(esConsumes()),
0053       runNum_(-1),
0054       side_(iConfig.getUntrackedParameter<int>("side", 3)),
0055       threshold_(iConfig.getUntrackedParameter<double>("amplitudeThreshold", 12.0)),
0056       minTimingAmp_(iConfig.getUntrackedParameter<double>("minimumTimingAmplitude", 0.100)) {
0057   vector<int> listDefaults;
0058   listDefaults.push_back(-1);
0059
0060   maskedChannels_ = iConfig.getUntrackedParameter<vector<int> >("maskedChannels", listDefaults);
0061   maskedFEDs_ = iConfig.getUntrackedParameter<vector<int> >("maskedFEDs", listDefaults);
0062   seedCrys_ = iConfig.getUntrackedParameter<vector<int> >("seedCrys", vector<int>());
0063
0064   vector<string> defaultMaskedEBs;
0065   defaultMaskedEBs.push_back("none");
0066   maskedEBs_ = iConfig.getUntrackedParameter<vector<string> >("maskedEBs", defaultMaskedEBs);
0067
0068   fedMap_ = new EcalFedMap();
0069
0070   string title1 = "Jitter for all FEDs";
0071   string name1 = "JitterAllFEDs";
0072   allFedsTimingHist_ = fileService->make<TH1F>(name1.c_str(), title1.c_str(), 150, -7, 7);
0073
0074   // load up the maskedFED list with the proper FEDids
0075   if (maskedFEDs_[0] == -1) {
0076     //if "actual" EB id given, then convert to FEDid and put in listFEDs_
0077     if (maskedEBs_[0] != "none") {
0078       maskedFEDs_.clear();
0079       for (vector<string>::const_iterator ebItr = maskedEBs_.begin(); ebItr != maskedEBs_.end(); ++ebItr) {
0080         maskedFEDs_.push_back(fedMap_->getFedFromSlice(*ebItr));
0081       }
0082     }
0083   }
0084
0085   for (int i = 0; i < 10; i++)
0086     abscissa[i] = i;
0087
0088   naiveEvtNum_ = 0;
0089 }
0090
0091 EcalMipGraphs::~EcalMipGraphs() {}
0092
0093 //
0094 // member functions
0095 //
0096
0097 // ------------ method called to for each event  ------------
0098 void EcalMipGraphs::analyze(edm::Event const& iEvent, edm::EventSetup const& iSetup) {
0099   // get the headers
0100   // (one header for each supermodule)
0101   edm::Handle<EcalRawDataCollection> DCCHeaders;
0102   iEvent.getByToken(rawDataToken_, DCCHeaders);
0103
0104   for (EcalRawDataCollection::const_iterator headerItr = DCCHeaders->begin(); headerItr != DCCHeaders->end();
0105        ++headerItr) {
0106     FEDsAndDCCHeaders_[headerItr->id() + 600] = *headerItr;
0107   }
0108
0109   int ievt = iEvent.id().event();
0110   naiveEvtNum_++;
0111
0112   if (runNum_ == -1) {
0113     runNum_ = iEvent.id().run();
0114     canvasNames_ = fileService->make<TTree>("canvasNames", "Names of written canvases");
0115     names = new std::vector<string>();
0116     canvasNames_->Branch("names", "vector<string>", &names);
0117   }
0118
0119   //We only want the 3x3's for this event...
0120   listEBChannels.clear();
0121   listEEChannels.clear();
0122   caloTopo_ = &iSetup.getData(topologyToken_);
0123
0124   Handle<EcalRecHitCollection> EBhits;
0125   Handle<EcalRecHitCollection> EEhits;
0126   iEvent.getByToken(ebRecHitToken_, EBhits);
0127   iEvent.getByToken(eeRecHitToken_, EEhits);
0128   // Now, retrieve the crystal digi from the event
0129   iEvent.getByToken(ebDigiToken_, EBdigisHandle);
0130   iEvent.getByToken(eeDigiToken_, EEdigisHandle);
0131
0132   selectHits(EBhits, ievt);
0133   selectHits(EEhits, ievt);
0134 }
0135
0136 TGraph* EcalMipGraphs::selectDigi(DetId thisDet, int ievt) {
0137   int emptyY[10];
0138   for (int i = 0; i < 10; i++)
0139     emptyY[i] = 0;
0140   TGraph* emptyGraph = fileService->make<TGraph>(10, abscissa, emptyY);
0141   emptyGraph->SetTitle("NOT ECAL");
0142
0143   //If the DetId is not from Ecal, return
0144   if (thisDet.det() != DetId::Ecal)
0145     return emptyGraph;
0146
0147   emptyGraph->SetTitle("NO DIGIS");
0148   //find digi we need  -- can't get find() to work; need DataFrame(DetId det) to work?
0149   EcalElectronicsId elecId = ecalElectronicsMap_->getElectronicsId(thisDet);
0150   int FEDid = 600 + elecId.dccId();
0151   bool isBarrel = true;
0152   if (FEDid < 610 || FEDid > 645)
0153     isBarrel = false;
0154   int cryIndex = isBarrel ? ((EBDetId)thisDet).hashedIndex() : getEEIndex(elecId);
0155   int ic = isBarrel ? ((EBDetId)thisDet).ic() : cryIndex;
0156
0157   string sliceName = fedMap_->getSliceFromFed(FEDid);
0158   EcalDataFrame df;
0159   if (isBarrel) {
0160     EBDigiCollection::const_iterator digiItr = EBdigisHandle->begin();
0161     while (digiItr != EBdigisHandle->end() && ((*digiItr).id() != (EBDetId)thisDet)) {
0162       ++digiItr;
0163     }
0164     if (digiItr == EBdigisHandle->end()) {
0165       //LogWarning("EcalMipGraphs") << "Cannot find digi for ic:" << ic
0166       //  << " FED:" << FEDid << " evt:" << naiveEvtNum_;
0167       return emptyGraph;
0168     } else
0169       df = *digiItr;
0170   } else {
0171     EEDigiCollection::const_iterator digiItr = EEdigisHandle->begin();
0172     while (digiItr != EEdigisHandle->end() && ((*digiItr).id() != (EEDetId)thisDet)) {
0173       ++digiItr;
0174     }
0175     if (digiItr == EEdigisHandle->end()) {
0176       //LogWarning("EcalMipGraphs") << "Cannot find digi for ic:" << ic
0177       //  << " FED:" << FEDid << " evt:" << naiveEvtNum_;
0178       return emptyGraph;
0179     } else
0180       df = *digiItr;
0181   }
0182
0183   int gainId = FEDsAndDCCHeaders_[FEDid].getMgpaGain();
0184   int gainHuman;
0185   if (gainId == 1)
0186     gainHuman = 12;
0187   else if (gainId == 2)
0188     gainHuman = 6;
0189   else if (gainId == 3)
0190     gainHuman = 1;
0191   else
0192     gainHuman = -1;
0193
0194   double pedestal = 200;
0195
0196   emptyGraph->SetTitle("FIRST TWO SAMPLES NOT GAIN12");
0197   if (df.sample(0).gainId() != 1 || df.sample(1).gainId() != 1)
0198     return emptyGraph;  //goes to the next digi
0199   else {
0200     ordinate[0] = df.sample(0).adc();
0201     ordinate[1] = df.sample(1).adc();
0202     pedestal = (double)(ordinate[0] + ordinate[1]) / (double)2;
0203   }
0204
0205   for (int i = 0; i < df.size(); ++i) {
0206     if (df.sample(i).gainId() != 0)
0207       ordinate[i] = (int)(pedestal + (df.sample(i).adc() - pedestal) * gainRatio[df.sample(i).gainId() - 1]);
0208     else
0209       ordinate[i] = 49152;  //Saturation of gain1
0210   }
0211
0212   TGraph* oneGraph = fileService->make<TGraph>(10, abscissa, ordinate);
0213   string name = "Graph_ev" + intToString(naiveEvtNum_) + "_ic" + intToString(ic) + "_FED" + intToString(FEDid);
0214   string gainString = (gainId == 1) ? "Free" : intToString(gainHuman);
0215   string title = "Event" + intToString(naiveEvtNum_) + "_lv1a" + intToString(ievt) + "_ic" + intToString(ic) + "_" +
0216                  sliceName + "_gain" + gainString;
0217
0218   float energy = 0;
0219   map<int, float>::const_iterator itr;
0220   itr = crysAndAmplitudesMap_.find(cryIndex);
0221   if (itr != crysAndAmplitudesMap_.end()) {
0222     //edm::LogWarning("EcalMipGraphs")<< "itr->second(ampli)="<< itr->second;
0223     energy = (float)itr->second;
0224   }
0225   //else
0226   //edm::LogWarning("EcalMipGraphs") << "cry " << ic << "not found in ampMap";
0227
0228   title += "_Energy" + floatToString(round(energy * 1000));
0229
0230   oneGraph->SetTitle(title.c_str());
0231   oneGraph->SetName(name.c_str());
0232   oneGraph->GetXaxis()->SetTitle("sample");
0233   oneGraph->GetYaxis()->SetTitle("ADC");
0234   return oneGraph;
0235 }
0236
0237 void EcalMipGraphs::selectHits(Handle<EcalRecHitCollection> hits, int ievt) {
0238   for (EcalRecHitCollection::const_iterator hitItr = hits->begin(); hitItr != hits->end(); ++hitItr) {
0239     EcalRecHit hit = (*hitItr);
0240     DetId det = hit.id();
0241     EcalElectronicsId elecId = ecalElectronicsMap_->getElectronicsId(det);
0242     int FEDid = 600 + elecId.dccId();
0243     bool isBarrel = true;
0244     if (FEDid < 610 || FEDid > 645)
0245       isBarrel = false;
0246     int cryIndex = isBarrel ? ((EBDetId)det).hashedIndex() : ((EEDetId)det).hashedIndex();
0247     int ic = isBarrel ? ((EBDetId)det).ic() : getEEIndex(elecId);
0248
0249     float ampli = hit.energy();
0250
0251     vector<int>::iterator result;
0252     result = find(maskedFEDs_.begin(), maskedFEDs_.end(), FEDid);
0253     if (result != maskedFEDs_.end()) {
0254       //LogWarning("EcalMipGraphs") << "skipping uncalRecHit for FED " << FEDid << " ; amplitude " << ampli;
0255       continue;
0256     }
0257     result = find(maskedChannels_.begin(), maskedChannels_.end(), cryIndex);
0258     if (result != maskedChannels_.end()) {
0259       //LogWarning("EcalMipGraphs") << "skipping uncalRecHit for channel: " << cryIndex << " in fed: " << FEDid << " with amplitude " << ampli ;
0260       continue;
0261     }
0262     bool cryIsInList = false;
0263     result = find(seedCrys_.begin(), seedCrys_.end(), cryIndex);
0264     if (result != seedCrys_.end())
0265       cryIsInList = true;
0266
0267     // Either we must have a user-requested cry (in which case there is no amplitude selection)
0268     // Or we pick all crys that pass the amplitude cut (in which case there is no fixed crystal selection)
0269     if (cryIsInList || (seedCrys_.empty() && ampli > threshold_)) {
0270       // We have a winner!
0271       crysAndAmplitudesMap_[cryIndex] = ampli;
0272       string name = "Event" + intToString(naiveEvtNum_) + "_ic" + intToString(ic) + "_FED" + intToString(FEDid);
0273       string title = "Digis";
0274       string seed = "ic" + intToString(ic) + "_FED" + intToString(FEDid);
0275       int freq = 1;
0276       pair<map<string, int>::iterator, bool> pair = seedFrequencyMap_.insert(make_pair(seed, freq));
0277       if (!pair.second) {
0278         ++(pair.first->second);
0279       }
0280
0281       TCanvas can(name.c_str(), title.c_str(), 200, 50, 900, 900);
0282       can.Divide(side_, side_);
0283       TGraph* myGraph;
0284       int canvasNum = 1;
0285
0286       CaloNavigator<DetId> cursor = CaloNavigator<DetId>(det, caloTopo_->getSubdetectorTopology(det));
0287       //Now put each graph in one by one
0288       for (int j = side_ / 2; j >= -side_ / 2; --j) {
0289         for (int i = -side_ / 2; i <= side_ / 2; ++i) {
0290           cursor.home();
0291           cursor.offsetBy(i, j);
0292           can.cd(canvasNum);
0293           myGraph = selectDigi(*cursor, ievt);
0294           myGraph->Draw("A*");
0295           canvasNum++;
0296         }
0297       }
0298       can.Write();
0299       names->push_back(name);
0300     }
0301
0302     TH1F* timingHist = FEDsAndTimingHists_[FEDid];
0303     if (timingHist == nullptr) {
0304       initHists(FEDid);
0305       timingHist = FEDsAndTimingHists_[FEDid];
0306     }
0307     if (ampli > minTimingAmp_) {
0308       timingHist->Fill(hit.time());
0309       allFedsTimingHist_->Fill(hit.time());
0310     }
0311   }
0312 }
0313
0314 int EcalMipGraphs::getEEIndex(EcalElectronicsId elecId) {
0315   int FEDid = 600 + elecId.dccId();
0316   return 10000 * FEDid + 100 * elecId.towerId() + 5 * (elecId.stripId() - 1) + elecId.xtalId();
0317 }
0318
0319 // insert the hist map into the map keyed by FED number
0320 void EcalMipGraphs::initHists(int FED) {
0321   using namespace std;
0322
0323   string title1 = "Jitter for ";
0324   title1.append(fedMap_->getSliceFromFed(FED));
0325   string name1 = "JitterFED";
0326   name1.append(intToString(FED));
0327   TH1F* timingHist = fileService->make<TH1F>(name1.c_str(), title1.c_str(), 150, -7, 7);
0328   FEDsAndTimingHists_[FED] = timingHist;
0329 }
0330
0331 // ------------ method called once each job just before starting event loop  ------------
0332 void EcalMipGraphs::beginRun(edm::Run const&, edm::EventSetup const& c) {
0333   ecalElectronicsMap_ = &c.getData(ecalMappingToken_);
0334 }
0335
0336 void EcalMipGraphs::endRun(edm::Run const&, edm::EventSetup const& c) {}
0337
0338 // ------------ method called once each job just after ending the event loop  ------------
0339 void EcalMipGraphs::endJob() {
0340   canvasNames_->Fill();
0341
0342   string frequencies = "";
0343   for (std::map<std::string, int>::const_iterator itr = seedFrequencyMap_.begin(); itr != seedFrequencyMap_.end();
0344        ++itr) {
0345     if (itr->second > 1) {
0346       frequencies += itr->first;
0347       frequencies += " Frequency: ";
0348       frequencies += intToString(itr->second);
0349       frequencies += "\n";
0350     }
0351   }
0352   LogWarning("EcalMipGraphs") << "Found seeds with frequency > 1: "
0353                               << "\n\n"
0354                               << frequencies;
0355
0356   std::string channels;
0357   for (std::vector<int>::const_iterator itr = maskedChannels_.begin(); itr != maskedChannels_.end(); ++itr) {
0358     channels += intToString(*itr);
0359     channels += ",";
0360   }
0361
0362   std::string feds;
0363   for (std::vector<int>::const_iterator itr = maskedFEDs_.begin(); itr != maskedFEDs_.end(); ++itr) {
0364     feds += intToString(*itr);
0365     feds += ",";
0366   }
0367
0368   LogWarning("EcalMipGraphs") << "Masked channels are: " << channels;
0369   LogWarning("EcalMipGraphs") << "Masked FEDs are: " << feds << " and that is all!";
0370 }
0371
0372 std::string EcalMipGraphs::intToString(int num) {
0373   using namespace std;
0374   ostringstream myStream;
0375   myStream << num << flush;
0376   return (myStream.str());  //returns the string form of the stringstream object
0377 }
0378
0379 std::string EcalMipGraphs::floatToString(float num) {
0380   using namespace std;
0381   ostringstream myStream;
0382   myStream << num << flush;
0383   return (myStream.str());  //returns the string form of the stringstream object
0384 }