EcalTools/plugins/EcalPulseShapeGrapher.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    EcalPulseShapeGrapher
0004 // Class:      EcalPulseShapeGrapher
0005 //
0006 /**\class EcalPulseShapeGrapher EcalPulseShapeGrapher.cc Analyzers/EcalPulseShapeGrapher/src/EcalPulseShapeGrapher.cc
0007
0008  Description: <one line class summary>
0009
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Seth Cooper
0015 //         Created:  Tue Feb  5 11:35:45 CST 2008
0016 //
0017 //
0018
0019 #include "CaloOnlineTools/EcalTools/plugins/EcalPulseShapeGrapher.h"
0020
0021 //
0022 // constants, enums and typedefs
0023 //
0024
0025 //
0026 // static data member definitions
0027 //
0028
0029 //
0030 // constructors and destructor
0031 //
0032 EcalPulseShapeGrapher::EcalPulseShapeGrapher(const edm::ParameterSet& iConfig)
0033     : EBUncalibratedRecHitCollection_(consumes<EcalUncalibratedRecHitCollection>(
0034           iConfig.getParameter<edm::InputTag>("EBUncalibratedRecHitCollection"))),
0035       EBDigis_(consumes<EBDigiCollection>(iConfig.getParameter<edm::InputTag>("EBDigiCollection"))),
0036       EEUncalibratedRecHitCollection_(consumes<EcalUncalibratedRecHitCollection>(
0037           iConfig.getParameter<edm::InputTag>("EEUncalibratedRecHitCollection"))),
0038       EEDigis_(consumes<EEDigiCollection>(iConfig.getParameter<edm::InputTag>("EEDigiCollection"))),
0039       ampCut_(iConfig.getUntrackedParameter<int>("AmplitudeCutADC", 13)),
0040       rootFilename_(iConfig.getUntrackedParameter<std::string>("rootFilename", "pulseShapeGrapher")) {
0041   //now do what ever initialization is needed
0042
0043   std::vector<int> listDefaults;
0044   listDefaults.push_back(-1);
0045   listChannels_ = iConfig.getUntrackedParameter<std::vector<int> >("listChannels", listDefaults);
0046
0047   for (std::vector<int>::const_iterator itr = listChannels_.begin(); itr != listChannels_.end(); ++itr) {
0048     std::string title = "Amplitude of cry " + intToString(*itr);
0049     std::string name = "ampOfCry" + intToString(*itr);
0050     ampHistMap_[*itr] = new TH1F(name.c_str(), title.c_str(), 100, 0, 100);
0051     ampHistMap_[*itr]->GetXaxis()->SetTitle("ADC");
0052
0053     title = "Amplitude (over 13 ADC) of cry " + intToString(*itr);
0054     name = "cutAmpOfCry" + intToString(*itr);
0055     cutAmpHistMap_[*itr] = new TH1F(name.c_str(), title.c_str(), 100, 0, 100);
0056     cutAmpHistMap_[*itr]->GetXaxis()->SetTitle("ADC");
0057
0058     title = "Pulse shape of cry " + intToString(*itr);
0059     name = "PulseShapeCry" + intToString(*itr);
0060     pulseShapeHistMap_[*itr] = new TH2F(name.c_str(), title.c_str(), 10, 0, 10, 220, -20, 2);
0061     pulseShapeHistMap_[*itr]->GetXaxis()->SetTitle("sample");
0062
0063     title = "Raw Pulse shape of cry " + intToString(*itr);
0064     name = "RawPulseShapeCry" + intToString(*itr);
0065     rawPulseShapeHistMap_[*itr] = new TH2F(name.c_str(), title.c_str(), 10, 0, 10, 500, 0, 500);
0066     rawPulseShapeHistMap_[*itr]->GetXaxis()->SetTitle("sample");
0067
0068     title = "Amplitude of first sample, cry " + intToString(*itr);
0069     name = "AmpOfFirstSampleCry" + intToString(*itr);
0070     firstSampleHistMap_[*itr] = new TH1F(name.c_str(), title.c_str(), 300, 100, 400);
0071     firstSampleHistMap_[*itr]->GetXaxis()->SetTitle("ADC");
0072   }
0073
0074   fedMap_ = new EcalFedMap();
0075
0076   for (int i = 0; i < 10; i++)
0077     abscissa[i] = i;
0078 }
0079
0080 //
0081 // member functions
0082 //
0083
0084 // ------------ method called to for each event  ------------
0085 void EcalPulseShapeGrapher::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0086   using namespace edm;
0087   using namespace std;
0088
0089   int numHitsWithActivity = 0;
0090
0091   const Handle<EcalUncalibratedRecHitCollection>& EBHits = iEvent.getHandle(EBUncalibratedRecHitCollection_);
0092   const Handle<EcalUncalibratedRecHitCollection>& EEHits = iEvent.getHandle(EEUncalibratedRecHitCollection_);
0093   const Handle<EBDigiCollection>& EBdigis = iEvent.getHandle(EBDigis_);
0094   const Handle<EEDigiCollection>& EEdigis = iEvent.getHandle(EEDigis_);
0095
0096   unique_ptr<EcalElectronicsMapping> ecalElectronicsMap(new EcalElectronicsMapping);
0097
0098   // Loop over the hits
0099   for (EcalUncalibratedRecHitCollection::const_iterator hitItr = EBHits->begin(); hitItr != EBHits->end(); ++hitItr) {
0100     EcalUncalibratedRecHit hit = (*hitItr);
0101     float amplitude = hit.amplitude();
0102     EBDetId hitDetId = hit.id();
0103
0104     // Get the Fedid
0105     EcalElectronicsId elecId = ecalElectronicsMap->getElectronicsId(hitDetId);
0106     int FEDid = 600 + elecId.dccId();
0107     string SMname = fedMap_->getSliceFromFed(FEDid);
0108
0109     vector<int>::const_iterator itr = listChannels_.begin();
0110     while (itr != listChannels_.end() && (*itr) != hitDetId.hashedIndex()) {
0111       itr++;
0112     }
0113     if (itr == listChannels_.end())
0114       continue;
0115
0116     ampHistMap_[hitDetId.hashedIndex()]->Fill(amplitude);
0117     if (amplitude < ampCut_)
0118       continue;
0119
0120     cutAmpHistMap_[hitDetId.hashedIndex()]->Fill(amplitude);
0121     numHitsWithActivity++;
0122     EBDigiCollection::const_iterator digiItr = EBdigis->begin();
0123     while (digiItr != EBdigis->end() && digiItr->id() != hitItr->id()) {
0124       digiItr++;
0125     }
0126     if (digiItr == EBdigis->end())
0127       continue;
0128
0129     double sampleADC[10];
0130     EBDataFrame df(*digiItr);
0131     double pedestal = 200;
0132
0133     if (df.sample(0).gainId() != 1 || df.sample(1).gainId() != 1)
0134       continue;  //goes to the next digi
0135     else {
0136       sampleADC[0] = df.sample(0).adc();
0137       sampleADC[1] = df.sample(1).adc();
0138       pedestal = (double)(sampleADC[0] + sampleADC[1]) / (double)2;
0139     }
0140
0141     for (int i = 0; (unsigned int)i < digiItr->size(); ++i) {
0142       EBDataFrame df(*digiItr);
0143       double gain = 12.;
0144       if (df.sample(i).gainId() == 1)
0145         gain = 1.;
0146       else if (df.sample(i).gainId() == 2)
0147         gain = 2.;
0148       sampleADC[i] = pedestal + (df.sample(i).adc() - pedestal) * gain;
0149     }
0150
0151     for (int i = 0; i < 10; ++i) {
0152       pulseShapeHistMap_[hitDetId.hashedIndex()]->Fill(i, (float)(sampleADC[i] - pedestal) / amplitude);
0153       rawPulseShapeHistMap_[hitDetId.hashedIndex()]->Fill(i, (float)(sampleADC[i]));
0154     }
0155     firstSampleHistMap_[hitDetId.hashedIndex()]->Fill(sampleADC[0]);
0156   }
0157
0158   // Now do the same for the EE hits
0159   for (EcalUncalibratedRecHitCollection::const_iterator hitItr = EEHits->begin(); hitItr != EEHits->end(); ++hitItr) {
0160     EcalUncalibratedRecHit hit = (*hitItr);
0161     float amplitude = hit.amplitude();
0162     EEDetId hitDetId = hit.id();
0163
0164     // Get the Fedid
0165     EcalElectronicsId elecId = ecalElectronicsMap->getElectronicsId(hitDetId);
0166     int FEDid = 600 + elecId.dccId();
0167     string SMname = fedMap_->getSliceFromFed(FEDid);
0168
0169     vector<int>::const_iterator itr = listChannels_.begin();
0170     while (itr != listChannels_.end() && (*itr) != hitDetId.hashedIndex()) {
0171       itr++;
0172     }
0173     if (itr == listChannels_.end())
0174       continue;
0175
0176     ampHistMap_[hitDetId.hashedIndex()]->Fill(amplitude);
0177     if (amplitude < ampCut_)
0178       continue;
0179
0180     cutAmpHistMap_[hitDetId.hashedIndex()]->Fill(amplitude);
0181     numHitsWithActivity++;
0182     EEDigiCollection::const_iterator digiItr = EEdigis->begin();
0183     while (digiItr != EEdigis->end() && digiItr->id() != hitItr->id()) {
0184       digiItr++;
0185     }
0186     if (digiItr == EEdigis->end())
0187       continue;
0188
0189     double sampleADC[10];
0190     EEDataFrame df(*digiItr);
0191     double pedestal = 200;
0192
0193     if (df.sample(0).gainId() != 1 || df.sample(1).gainId() != 1)
0194       continue;  //goes to the next digi
0195     else {
0196       sampleADC[0] = df.sample(0).adc();
0197       sampleADC[1] = df.sample(1).adc();
0198       pedestal = (double)(sampleADC[0] + sampleADC[1]) / (double)2;
0199     }
0200
0201     for (int i = 0; (unsigned int)i < digiItr->size(); ++i) {
0202       EEDataFrame df(*digiItr);
0203       double gain = 12.;
0204       if (df.sample(i).gainId() == 1)
0205         gain = 1.;
0206       else if (df.sample(i).gainId() == 2)
0207         gain = 2.;
0208       sampleADC[i] = pedestal + (df.sample(i).adc() - pedestal) * gain;
0209     }
0210
0211     for (int i = 0; i < 10; ++i) {
0212       pulseShapeHistMap_[hitDetId.hashedIndex()]->Fill(i, (float)(sampleADC[i] - pedestal) / amplitude);
0213       rawPulseShapeHistMap_[hitDetId.hashedIndex()]->Fill(i, (float)(sampleADC[i]));
0214     }
0215     firstSampleHistMap_[hitDetId.hashedIndex()]->Fill(sampleADC[0]);
0216   }
0217 }
0218
0219 // ------------ method called once each job just after ending the event loop  ------------
0220 void EcalPulseShapeGrapher::endJob() {
0221   rootFilename_ += ".root";
0222   file_ = new TFile(rootFilename_.c_str(), "RECREATE");
0223   TH1::AddDirectory(false);
0224
0225   for (std::vector<int>::const_iterator itr = listChannels_.begin(); itr != listChannels_.end(); ++itr) {
0226     ampHistMap_[*itr]->Write();
0227     cutAmpHistMap_[*itr]->Write();
0228     firstSampleHistMap_[*itr]->Write();
0229
0230     rawPulseShapeHistMap_[*itr]->Write();
0231     TProfile* t2 = (TProfile*)(rawPulseShapeHistMap_[*itr]->ProfileX());
0232     t2->Write();
0233     //TODO: fix the normalization so these are correct
0234   }
0235
0236   file_->Write();
0237   file_->Close();
0238 }
0239
0240 std::string EcalPulseShapeGrapher::intToString(int num) {
0241   using namespace std;
0242   ostringstream myStream;
0243   myStream << num << flush;
0244   return (myStream.str());  //returns the string form of the stringstream object
0245 }