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_(iConfig.getParameter<edm::InputTag>("EBUncalibratedRecHitCollection")),
0034       EBDigis_(iConfig.getParameter<edm::InputTag>("EBDigiCollection")),
0035       EEUncalibratedRecHitCollection_(iConfig.getParameter<edm::InputTag>("EEUncalibratedRecHitCollection")),
0036       EEDigis_(iConfig.getParameter<edm::InputTag>("EEDigiCollection")),
0037       ampCut_(iConfig.getUntrackedParameter<int>("AmplitudeCutADC", 13)),
0038       rootFilename_(iConfig.getUntrackedParameter<std::string>("rootFilename", "pulseShapeGrapher")) {
0039   //now do what ever initialization is needed
0040
0041   std::vector<int> listDefaults;
0042   listDefaults.push_back(-1);
0043   listChannels_ = iConfig.getUntrackedParameter<std::vector<int> >("listChannels", listDefaults);
0044
0045   for (std::vector<int>::const_iterator itr = listChannels_.begin(); itr != listChannels_.end(); ++itr) {
0046     std::string title = "Amplitude of cry " + intToString(*itr);
0047     std::string name = "ampOfCry" + intToString(*itr);
0048     ampHistMap_[*itr] = new TH1F(name.c_str(), title.c_str(), 100, 0, 100);
0049     ampHistMap_[*itr]->GetXaxis()->SetTitle("ADC");
0050
0051     title = "Amplitude (over 13 ADC) of cry " + intToString(*itr);
0052     name = "cutAmpOfCry" + intToString(*itr);
0053     cutAmpHistMap_[*itr] = new TH1F(name.c_str(), title.c_str(), 100, 0, 100);
0054     cutAmpHistMap_[*itr]->GetXaxis()->SetTitle("ADC");
0055
0056     title = "Pulse shape of cry " + intToString(*itr);
0057     name = "PulseShapeCry" + intToString(*itr);
0058     pulseShapeHistMap_[*itr] = new TH2F(name.c_str(), title.c_str(), 10, 0, 10, 220, -20, 2);
0059     pulseShapeHistMap_[*itr]->GetXaxis()->SetTitle("sample");
0060
0061     title = "Raw Pulse shape of cry " + intToString(*itr);
0062     name = "RawPulseShapeCry" + intToString(*itr);
0063     rawPulseShapeHistMap_[*itr] = new TH2F(name.c_str(), title.c_str(), 10, 0, 10, 500, 0, 500);
0064     rawPulseShapeHistMap_[*itr]->GetXaxis()->SetTitle("sample");
0065
0066     title = "Amplitude of first sample, cry " + intToString(*itr);
0067     name = "AmpOfFirstSampleCry" + intToString(*itr);
0068     firstSampleHistMap_[*itr] = new TH1F(name.c_str(), title.c_str(), 300, 100, 400);
0069     firstSampleHistMap_[*itr]->GetXaxis()->SetTitle("ADC");
0070   }
0071
0072   fedMap_ = new EcalFedMap();
0073
0074   for (int i = 0; i < 10; i++)
0075     abscissa[i] = i;
0076 }
0077
0078 EcalPulseShapeGrapher::~EcalPulseShapeGrapher() {
0079   // do anything here that needs to be done at desctruction time
0080   // (e.g. close files, deallocate resources etc.)
0081 }
0082
0083 //
0084 // member functions
0085 //
0086
0087 // ------------ method called to for each event  ------------
0088 void EcalPulseShapeGrapher::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0089   using namespace edm;
0090   using namespace std;
0091
0092   int numHitsWithActivity = 0;
0093
0094   //int eventNum = iEvent.id().event();
0095   //vector<EcalUncalibratedRecHit> sampleHitsPastCut;
0096
0097   Handle<EcalUncalibratedRecHitCollection> EBHits;
0098   iEvent.getByLabel(EBUncalibratedRecHitCollection_, EBHits);
0099   Handle<EcalUncalibratedRecHitCollection> EEHits;
0100   iEvent.getByLabel(EEUncalibratedRecHitCollection_, EEHits);
0101   //cout << "event: " << eventNum << " sample hits collection size: " << sampleHits->size() << endl;
0102   //Handle<EcalUncalibratedRecHitCollection> fittedHits;
0103   //iEvent.getByLabel(FittedUncalibratedRecHitCollection_, fittedHits);
0104   Handle<EBDigiCollection> EBdigis;
0105   iEvent.getByLabel(EBDigis_, EBdigis);
0106   Handle<EEDigiCollection> EEdigis;
0107   iEvent.getByLabel(EEDigis_, EEdigis);
0108   //cout << "event: " << eventNum << " digi collection size: " << digis->size() << endl;
0109
0110   unique_ptr<EcalElectronicsMapping> ecalElectronicsMap(new EcalElectronicsMapping);
0111
0112   // Loop over the hits
0113   for (EcalUncalibratedRecHitCollection::const_iterator hitItr = EBHits->begin(); hitItr != EBHits->end(); ++hitItr) {
0114     EcalUncalibratedRecHit hit = (*hitItr);
0115     float amplitude = hit.amplitude();
0116     EBDetId hitDetId = hit.id();
0117
0118     // Get the Fedid
0119     EcalElectronicsId elecId = ecalElectronicsMap->getElectronicsId(hitDetId);
0120     int FEDid = 600 + elecId.dccId();
0121     string SMname = fedMap_->getSliceFromFed(FEDid);
0122
0123     vector<int>::const_iterator itr = listChannels_.begin();
0124     while (itr != listChannels_.end() && (*itr) != hitDetId.hashedIndex()) {
0125       itr++;
0126     }
0127     if (itr == listChannels_.end())
0128       continue;
0129
0130     ampHistMap_[hitDetId.hashedIndex()]->Fill(amplitude);
0131     //cout << "Cry hash:" << hitDetId.hashedIndex() << " amplitude: " << amplitude << endl;
0132     if (amplitude < ampCut_)
0133       continue;
0134
0135     cutAmpHistMap_[hitDetId.hashedIndex()]->Fill(amplitude);
0136     numHitsWithActivity++;
0137     EBDigiCollection::const_iterator digiItr = EBdigis->begin();
0138     while (digiItr != EBdigis->end() && digiItr->id() != hitItr->id()) {
0139       digiItr++;
0140     }
0141     if (digiItr == EBdigis->end())
0142       continue;
0143
0144     double sampleADC[10];
0145     EBDataFrame df(*digiItr);
0146     double pedestal = 200;
0147
0148     if (df.sample(0).gainId() != 1 || df.sample(1).gainId() != 1)
0149       continue;  //goes to the next digi
0150     else {
0151       sampleADC[0] = df.sample(0).adc();
0152       sampleADC[1] = df.sample(1).adc();
0153       pedestal = (double)(sampleADC[0] + sampleADC[1]) / (double)2;
0154     }
0155
0156     for (int i = 0; (unsigned int)i < digiItr->size(); ++i) {
0157       EBDataFrame df(*digiItr);
0158       double gain = 12.;
0159       if (df.sample(i).gainId() == 1)
0160         gain = 1.;
0161       else if (df.sample(i).gainId() == 2)
0162         gain = 2.;
0163       sampleADC[i] = pedestal + (df.sample(i).adc() - pedestal) * gain;
0164     }
0165
0166     //cout << "1) maxsample amp:" << maxSampleAmp << " maxSampleIndex:" << maxSampleIndex << endl;
0167     for (int i = 0; i < 10; ++i) {
0168       //cout << "Filling hist for:" << hitDetId.hashedIndex() << " with sample:" << i
0169       //<< " amp:" <<(float)(df.sample(i).adc()-baseline)/(maxSampleAmp-baseline) << endl;
0170       //cout << "ADC of sample:" << df.sample(i).adc() << " baseline:" << baseline << " maxSampleAmp:"
0171       //<< maxSampleAmp << endl << endl;
0172       pulseShapeHistMap_[hitDetId.hashedIndex()]->Fill(i, (float)(sampleADC[i] - pedestal) / amplitude);
0173       rawPulseShapeHistMap_[hitDetId.hashedIndex()]->Fill(i, (float)(sampleADC[i]));
0174     }
0175     firstSampleHistMap_[hitDetId.hashedIndex()]->Fill(sampleADC[0]);
0176   }
0177
0178   // Now do the same for the EE hits
0179   for (EcalUncalibratedRecHitCollection::const_iterator hitItr = EEHits->begin(); hitItr != EEHits->end(); ++hitItr) {
0180     EcalUncalibratedRecHit hit = (*hitItr);
0181     float amplitude = hit.amplitude();
0182     EEDetId hitDetId = hit.id();
0183
0184     // Get the Fedid
0185     EcalElectronicsId elecId = ecalElectronicsMap->getElectronicsId(hitDetId);
0186     int FEDid = 600 + elecId.dccId();
0187     string SMname = fedMap_->getSliceFromFed(FEDid);
0188
0189     vector<int>::const_iterator itr = listChannels_.begin();
0190     while (itr != listChannels_.end() && (*itr) != hitDetId.hashedIndex()) {
0191       itr++;
0192     }
0193     if (itr == listChannels_.end())
0194       continue;
0195
0196     ampHistMap_[hitDetId.hashedIndex()]->Fill(amplitude);
0197     //cout << "Cry hash:" << hitDetId.hashedIndex() << " amplitude: " << amplitude << endl;
0198     if (amplitude < ampCut_)
0199       continue;
0200
0201     cutAmpHistMap_[hitDetId.hashedIndex()]->Fill(amplitude);
0202     numHitsWithActivity++;
0203     EEDigiCollection::const_iterator digiItr = EEdigis->begin();
0204     while (digiItr != EEdigis->end() && digiItr->id() != hitItr->id()) {
0205       digiItr++;
0206     }
0207     if (digiItr == EEdigis->end())
0208       continue;
0209
0210     double sampleADC[10];
0211     EEDataFrame df(*digiItr);
0212     double pedestal = 200;
0213
0214     if (df.sample(0).gainId() != 1 || df.sample(1).gainId() != 1)
0215       continue;  //goes to the next digi
0216     else {
0217       sampleADC[0] = df.sample(0).adc();
0218       sampleADC[1] = df.sample(1).adc();
0219       pedestal = (double)(sampleADC[0] + sampleADC[1]) / (double)2;
0220     }
0221
0222     for (int i = 0; (unsigned int)i < digiItr->size(); ++i) {
0223       EEDataFrame df(*digiItr);
0224       double gain = 12.;
0225       if (df.sample(i).gainId() == 1)
0226         gain = 1.;
0227       else if (df.sample(i).gainId() == 2)
0228         gain = 2.;
0229       sampleADC[i] = pedestal + (df.sample(i).adc() - pedestal) * gain;
0230     }
0231
0232     //cout << "1) maxsample amp:" << maxSampleAmp << " maxSampleIndex:" << maxSampleIndex << endl;
0233     for (int i = 0; i < 10; ++i) {
0234       //cout << "Filling hist for:" << hitDetId.hashedIndex() << " with sample:" << i
0235       //<< " amp:" <<(float)(df.sample(i).adc()-baseline)/(maxSampleAmp-baseline) << endl;
0236       //cout << "ADC of sample:" << df.sample(i).adc() << " baseline:" << baseline << " maxSampleAmp:"
0237       //<< maxSampleAmp << endl << endl;
0238       pulseShapeHistMap_[hitDetId.hashedIndex()]->Fill(i, (float)(sampleADC[i] - pedestal) / amplitude);
0239       rawPulseShapeHistMap_[hitDetId.hashedIndex()]->Fill(i, (float)(sampleADC[i]));
0240     }
0241     firstSampleHistMap_[hitDetId.hashedIndex()]->Fill(sampleADC[0]);
0242   }
0243 }
0244
0245 // ------------ method called once each job just after ending the event loop  ------------
0246 void EcalPulseShapeGrapher::endJob() {
0247   rootFilename_ += ".root";
0248   file_ = new TFile(rootFilename_.c_str(), "RECREATE");
0249   TH1::AddDirectory(false);
0250
0251   for (std::vector<int>::const_iterator itr = listChannels_.begin(); itr != listChannels_.end(); ++itr) {
0252     ampHistMap_[*itr]->Write();
0253     cutAmpHistMap_[*itr]->Write();
0254     firstSampleHistMap_[*itr]->Write();
0255
0256     rawPulseShapeHistMap_[*itr]->Write();
0257     TProfile* t2 = (TProfile*)(rawPulseShapeHistMap_[*itr]->ProfileX());
0258     t2->Write();
0259     //TODO: fix the normalization so these are correct
0260     //pulseShapeHistMap_[*itr]->Write();
0261     //TProfile* t1 = (TProfile*) (pulseShapeHistMap_[*itr]->ProfileX());
0262     //t1->Write();
0263   }
0264
0265   file_->Write();
0266   file_->Close();
0267 }
0268
0269 std::string EcalPulseShapeGrapher::intToString(int num) {
0270   using namespace std;
0271   ostringstream myStream;
0272   myStream << num << flush;
0273   return (myStream.str());  //returns the string form of the stringstream object
0274 }