EcalTools/plugins/EcalPedHists.cc

0001 /**
0002  * Module which outputs a root file of ADC counts (all three gains)
0003  *   of user-selected channels (defaults to channel 1) for
0004  *   user-selected samples (defaults to samples 1,2,3) for
0005  *   user-selected supermodules.
0006  *
0007  * \author S. Cooper
0008  *
0009  */
0010
0011 #include "CaloOnlineTools/EcalTools/plugins/EcalPedHists.h"
0012
0013 EcalPedHists::EcalPedHists(const edm::ParameterSet& ps)
0014     : runNum_(-1),
0015       fileName_(ps.getUntrackedParameter<std::string>("fileName", std::string("ecalPedDigiDump"))),
0016       barrelDigiCollection_(ps.getParameter<edm::InputTag>("EBdigiCollection")),
0017       endcapDigiCollection_(ps.getParameter<edm::InputTag>("EEdigiCollection")),
0018       headerProducer_(ps.getParameter<edm::InputTag>("headerProducer")),
0019       rawDataToken_(consumes<EcalRawDataCollection>(headerProducer_)),
0020       ebDigiToken_(consumes<EBDigiCollection>(barrelDigiCollection_)),
0021       eeDigiToken_(consumes<EEDigiCollection>(endcapDigiCollection_)),
0022       ecalMappingToken_(esConsumes<edm::Transition::BeginRun>()) {
0023   using namespace std;
0024
0025   fedMap_ = new EcalFedMap();
0026   histsFilled_ = false;
0027   //for(int i=601; i<655; ++i)
0028   //{
0029   //  listDefaults.push_back(i);
0030   //}
0031   listFEDs_ = ps.getUntrackedParameter<vector<int> >("listFEDs");
0032   listEBs_ = ps.getUntrackedParameter<vector<string> >("listEBs");
0033
0034   if (listFEDs_.empty()) {
0035     allFEDsSelected_ = false;
0036     //if "actual" EB id given, then convert to FEDid and put in listFEDs_
0037     if (!listEBs_.empty()) {
0038       listFEDs_.clear();
0039       for (vector<string>::const_iterator itr = listEBs_.begin(); itr != listEBs_.end(); ++itr) {
0040         listFEDs_.push_back(fedMap_->getFedFromSlice(*itr));
0041       }
0042     }
0043   } else if (listFEDs_[0] == -1) {
0044     // Apply no selection if -1 is passed in FED list
0045     allFEDsSelected_ = true;
0046     //debug
0047     //cout << "no selection on FEDs!" << endl;
0048     //inputIsOk_=false;
0049     //return;
0050     //listFEDs_ = listDefaults;
0051   } else {
0052     //in this case, listFEDs should be populated
0053     allFEDsSelected_ = false;
0054   }
0055
0056   if (!allFEDsSelected_) {
0057     // Verify FED numbers are valid
0058     for (vector<int>::const_iterator intIter = listFEDs_.begin(); intIter != listFEDs_.end(); intIter++) {
0059       if (((*intIter) < 601) || (654 < (*intIter))) {
0060         cout << "[EcalPedHists] FED value: " << (*intIter) << " found in listFEDs. "
0061              << " Valid range is 601-654. Returning." << endl;
0062         inputIsOk_ = false;
0063         return;
0064       } else
0065         theRealFedSet_.insert(*intIter);
0066     }
0067   }
0068
0069   vector<int> listDefaults = vector<int>();
0070   listDefaults.clear();
0071   for (int i = 1; i < 1701; ++i) {
0072     listDefaults.push_back(i);
0073   }
0074   listChannels_ = ps.getUntrackedParameter<vector<int> >("listChannels", listDefaults);
0075   listDefaults.clear();
0076   // Get samples to plot (default to 1,2,3)
0077   listDefaults.push_back(0);
0078   listDefaults.push_back(1);
0079   listDefaults.push_back(2);
0080   listSamples_ = ps.getUntrackedParameter<vector<int> >("listSamples", listDefaults);
0081
0082   inputIsOk_ = true;
0083   vector<int>::iterator intIter;
0084
0085   // Verify crystal numbers are valid
0086   for (intIter = listChannels_.begin(); intIter != listChannels_.end(); ++intIter) {
0087     //TODO: Fix crystal index checking?
0088     //if ( ((*intIter) < 1)||(1700 < (*intIter)) )
0089     //{
0090     //  cout << "[EcalPedHists] ic value: " << (*intIter) << " found in listChannels. "
0091     //        << " Valid range is 1-1700. Returning." << endl;
0092     //  inputIsOk_ = false;
0093     //  return;
0094     //}
0095   }
0096   // Verify sample numbers are valid
0097   for (intIter = listSamples_.begin(); intIter != listSamples_.end(); intIter++) {
0098     if (((*intIter) < 1) || (10 < (*intIter))) {
0099       cout << "[EcalPedHists] sample number: " << (*intIter) << " found in listSamples. "
0100            << " Valid range is 1-10. Returning." << endl;
0101       inputIsOk_ = false;
0102       return;
0103     }
0104   }
0105 }
0106
0107 EcalPedHists::~EcalPedHists() {}
0108
0109 void EcalPedHists::beginRun(edm::Run const&, edm::EventSetup const& c) {
0110   ecalElectronicsMap_ = &c.getData(ecalMappingToken_);
0111 }
0112
0113 void EcalPedHists::endRun(edm::Run const&, edm::EventSetup const& c) {}
0114
0115 void EcalPedHists::endJob(void) {
0116   using namespace std;
0117   if (inputIsOk_) {
0118     //debug
0119     //cout << "endJob:creating root file!" << endl;
0120
0121     fileName_ += "-" + intToString(runNum_) + ".graph.root";
0122
0123     TFile root_file_(fileName_.c_str(), "RECREATE");
0124     //Loop over FEDs first
0125     for (set<int>::const_iterator FEDitr = theRealFedSet_.begin(); FEDitr != theRealFedSet_.end(); ++FEDitr) {
0126       if (!histsFilled_)
0127         break;
0128       string dir = fedMap_->getSliceFromFed(*FEDitr);
0129       TDirectory* FEDdir = gDirectory->mkdir(dir.c_str());
0130       FEDdir->cd();
0131       //root_file_.mkdir(dir.c_str());
0132       //root_file_.cd(dir.c_str());
0133       map<string, TH1F*> mapHistos = FEDsAndHistMaps_[*FEDitr];
0134
0135       //Loop over channels; write histos and directory structure
0136       for (vector<int>::const_iterator itr = listChannels_.begin(); itr != listChannels_.end(); itr++) {
0137         //debug
0138         //cout << "loop over channels" << endl;
0139
0140         TH1F* hist = nullptr;
0141         string chnl = intToString(*itr);
0142         string name1 = "Cry";
0143         name1.append(chnl + "Gain1");
0144         string name2 = "Cry";
0145         name2.append(chnl + "Gain6");
0146         string name3 = "Cry";
0147         name3.append(chnl + "Gain12");
0148         hist = mapHistos[name1];
0149         // This is a sanity check only
0150         if (hist != nullptr) {
0151           string cryDirName = "Cry_" + chnl;
0152           TDirectory* cryDir = FEDdir->mkdir(cryDirName.c_str());
0153           cryDir->cd();
0154           hist->SetDirectory(cryDir);
0155           hist->Write();
0156           hist = mapHistos[name2];
0157           hist->SetDirectory(cryDir);
0158           hist->Write();
0159           hist = mapHistos[name3];
0160           hist->SetDirectory(cryDir);
0161           hist->Write();
0162           //root_file_.cd(dir.c_str());
0163           root_file_.cd();
0164         } else {
0165           cerr << "EcalPedHists: Error: This shouldn't happen!" << endl;
0166         }
0167       }
0168       root_file_.cd();
0169     }
0170     root_file_.Close();
0171   }
0172 }
0173
0174 void EcalPedHists::analyze(const edm::Event& e, const edm::EventSetup& c) {
0175   using namespace std;
0176   using namespace edm;
0177
0178   if (!inputIsOk_)
0179     return;
0180
0181   // loop over the headers, this is to detect missing FEDs if all are selected
0182   if (allFEDsSelected_) {
0183     edm::Handle<EcalRawDataCollection> DCCHeaders;
0184     e.getByToken(rawDataToken_, DCCHeaders);
0185     if (!DCCHeaders.isValid()) {
0186       edm::LogError("EcalPedHists") << "Error! can't get the product " << headerProducer_;
0187       return;
0188     }
0189
0190     for (EcalRawDataCollection::const_iterator headerItr = DCCHeaders->begin(); headerItr != DCCHeaders->end();
0191          ++headerItr) {
0192       int FEDid = 600 + headerItr->id();
0193       theRealFedSet_.insert(FEDid);
0194     }
0195   }
0196
0197   // loop over fed list and make sure that there are histo maps
0198   for (set<int>::const_iterator fedItr = theRealFedSet_.begin(); fedItr != theRealFedSet_.end(); ++fedItr) {
0199     if (FEDsAndHistMaps_.find(*fedItr) == FEDsAndHistMaps_.end())
0200       initHists(*fedItr);
0201   }
0202
0203   //debug
0204   //cout << "analyze...input is ok? " << inputIsOk_ << endl;
0205
0206   bool barrelDigisFound = true;
0207   bool endcapDigisFound = true;
0208   // get the barrel digis
0209   // (one digi for each crystal)
0210   // TODO; SIC: fix this behavior
0211   Handle<EBDigiCollection> barrelDigis;
0212   e.getByToken(ebDigiToken_, barrelDigis);
0213   if (!barrelDigis.isValid()) {
0214     edm::LogError("EcalPedOffset") << "Error! can't get the product " << barrelDigiCollection_;
0215     barrelDigisFound = false;
0216   }
0217   // get the endcap digis
0218   // (one digi for each crystal)
0219   // TODO; SIC: fix this behavior
0220   Handle<EEDigiCollection> endcapDigis;
0221   e.getByToken(eeDigiToken_, endcapDigis);
0222   if (!endcapDigis.isValid()) {
0223     edm::LogError("EcalPedOffset") << "Error! can't get the product " << endcapDigiCollection_;
0224     endcapDigisFound = false;
0225   }
0226
0227   if (barrelDigisFound)
0228     readEBdigis(barrelDigis);
0229   if (endcapDigisFound)
0230     readEEdigis(endcapDigis);
0231   if (!barrelDigisFound && !endcapDigisFound)
0232     edm::LogError("EcalPedOffset") << "No digis found in the event!";
0233
0234   if (runNum_ == -1)
0235     runNum_ = e.id().run();
0236 }
0237
0238 // insert the 3-entry hist map into the map keyed by FED number
0239 void EcalPedHists::initHists(int FED) {
0240   using namespace std;
0241   //using namespace edm;
0242
0243   std::map<string, TH1F*> histMap;
0244   //debug
0245   //cout << "Initializing map for FED:" << *FEDitr << endl;
0246   for (vector<int>::const_iterator intIter = listChannels_.begin(); intIter != listChannels_.end(); ++intIter) {
0247     //Put 3 histos (1 per gain) for the channel into the map
0248     string FEDid = intToString(FED);
0249     string chnl = intToString(*intIter);
0250     string title1 = "Gain1 ADC Counts for channel ";
0251     title1.append(chnl);
0252     string name1 = "Cry";
0253     name1.append(chnl + "Gain1");
0254     string title2 = "Gain6 ADC Counts for channel ";
0255     title2.append(chnl);
0256     string name2 = "Cry";
0257     name2.append(chnl + "Gain6");
0258     string title3 = "Gain12 ADC Counts for channel ";
0259     title3.append(chnl);
0260     string name3 = "Cry";
0261     name3.append(chnl + "Gain12");
0262     histMap.insert(make_pair(name1, new TH1F(name1.c_str(), title1.c_str(), 75, 175.0, 250.0)));
0263     histMap[name1]->SetDirectory(nullptr);
0264     histMap.insert(make_pair(name2, new TH1F(name2.c_str(), title2.c_str(), 75, 175.0, 250.0)));
0265     histMap[name2]->SetDirectory(nullptr);
0266     histMap.insert(make_pair(name3, new TH1F(name3.c_str(), title3.c_str(), 75, 175.0, 250.0)));
0267     histMap[name3]->SetDirectory(nullptr);
0268   }
0269   FEDsAndHistMaps_.insert(make_pair(FED, histMap));
0270 }
0271
0272 void EcalPedHists::readEBdigis(edm::Handle<EBDigiCollection> digis) {
0273   using namespace std;
0274   using namespace edm;
0275   //debug
0276   //cout << "readEBdigis" << endl;
0277
0278   // Loop over digis
0279   for (EBDigiCollection::const_iterator digiItr = digis->begin(); digiItr != digis->end(); ++digiItr) {
0280     EBDetId detId = EBDetId(digiItr->id());
0281     EcalElectronicsId elecId = ecalElectronicsMap_->getElectronicsId(detId);
0282     int FEDid = 600 + elecId.dccId();
0283     int crystalId = detId.ic();
0284
0285     //debug
0286     //cout << "FEDid:" << FEDid << " cryId:" << crystalId << endl;
0287     //cout << "FEDid:" << FEDid << endl;
0288     //Select desired supermodules only
0289     set<int>::const_iterator fedIter = find(theRealFedSet_.begin(), theRealFedSet_.end(), FEDid);
0290     if (fedIter == theRealFedSet_.end())
0291       continue;
0292
0293     // Select desired channels only
0294     vector<int>::iterator icIter;
0295     icIter = find(listChannels_.begin(), listChannels_.end(), crystalId);
0296     if (icIter == listChannels_.end())
0297       continue;
0298
0299     // Get the adc counts from the selected samples and fill the corresponding histogram
0300     // Must subtract 1 from user-given sample list (e.g., user's sample 1 -> sample 0)
0301     for (vector<int>::iterator itr = listSamples_.begin(); itr != listSamples_.end(); itr++) {
0302       histsFilled_ = true;
0303       map<string, TH1F*> mapHistos = FEDsAndHistMaps_[FEDid];
0304       string chnl = intToString(crystalId);
0305       string name1 = "Cry";
0306       name1.append(chnl + "Gain1");
0307       string name2 = "Cry";
0308       name2.append(chnl + "Gain6");
0309       string name3 = "Cry";
0310       name3.append(chnl + "Gain12");
0311       TH1F* hist = nullptr;
0312       if (((EBDataFrame)(*digiItr)).sample(*itr - 1).gainId() == 3)
0313         hist = mapHistos[name1];
0314       if (((EBDataFrame)(*digiItr)).sample(*itr - 1).gainId() == 2)
0315         hist = mapHistos[name2];
0316       if (((EBDataFrame)(*digiItr)).sample(*itr - 1).gainId() == 1)
0317         hist = mapHistos[name3];
0318       if (hist != nullptr)
0319         hist->Fill(((EBDataFrame)(*digiItr)).sample(*itr - 1).adc());
0320       else
0321         cerr << "EcalPedHistDumper: Error: This shouldn't happen!" << endl;
0322     }
0323   }
0324 }
0325
0326 void EcalPedHists::readEEdigis(edm::Handle<EEDigiCollection> digis) {
0327   using namespace std;
0328   using namespace edm;
0329   //debug
0330   //cout << "readEEdigis" << endl;
0331
0332   // Loop over digis
0333   for (EEDigiCollection::const_iterator digiItr = digis->begin(); digiItr != digis->end(); ++digiItr) {
0334     EEDetId detId = EEDetId(digiItr->id());
0335     EcalElectronicsId elecId = ecalElectronicsMap_->getElectronicsId(detId);
0336     int FEDid = 600 + elecId.dccId();
0337     int crystalId = 10000 * FEDid + 100 * elecId.towerId() + 5 * (elecId.stripId() - 1) + elecId.xtalId();
0338
0339     //Select desired FEDs only
0340     set<int>::const_iterator fedIter = find(theRealFedSet_.begin(), theRealFedSet_.end(), FEDid);
0341     if (fedIter == theRealFedSet_.end())
0342       continue;
0343
0344     // Select desired channels only
0345     vector<int>::iterator icIter;
0346     icIter = find(listChannels_.begin(), listChannels_.end(), crystalId);
0347     if (icIter == listChannels_.end())
0348       continue;
0349
0350     // Get the adc counts from the selected samples and fill the corresponding histogram
0351     // Must subtract 1 from user-given sample list (e.g., user's sample 1 -> sample 0)
0352     for (vector<int>::iterator itr = listSamples_.begin(); itr != listSamples_.end(); itr++) {
0353       histsFilled_ = true;
0354       map<string, TH1F*> mapHistos = FEDsAndHistMaps_[FEDid];
0355       string chnl = intToString(crystalId);
0356       string name1 = "Cry";
0357       name1.append(chnl + "Gain1");
0358       string name2 = "Cry";
0359       name2.append(chnl + "Gain6");
0360       string name3 = "Cry";
0361       name3.append(chnl + "Gain12");
0362       TH1F* hist = nullptr;
0363       if (((EBDataFrame)(*digiItr)).sample(*itr - 1).gainId() == 3)
0364         hist = mapHistos[name1];
0365       if (((EBDataFrame)(*digiItr)).sample(*itr - 1).gainId() == 2)
0366         hist = mapHistos[name2];
0367       if (((EBDataFrame)(*digiItr)).sample(*itr - 1).gainId() == 1)
0368         hist = mapHistos[name3];
0369       if (hist != nullptr)
0370         hist->Fill(((EBDataFrame)(*digiItr)).sample(*itr - 1).adc());
0371       else
0372         cerr << "EcalPedHistDumper: Error: This shouldn't happen!" << endl;
0373     }
0374   }
0375 }
0376
0377 std::string EcalPedHists::intToString(int num) {
0378   using namespace std;
0379   //
0380   // outputs the number into the string stream and then flushes
0381   // the buffer (makes sure the output is put into the stream)
0382   //
0383   ostringstream myStream;  //creates an ostringstream object
0384   myStream << num << flush;
0385   return (myStream.str());  //returns the string form of the stringstream object
0386 }