File indexing completed on 2024-04-06 12:27:56
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012 #include "RecoTBCalo/EcalSimpleTBAnalysis/interface/EcalSimpleTBAnalyzer.h"
0013 #include "DataFormats/EcalRecHit/interface/EcalUncalibratedRecHit.h"
0014 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0015 #include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
0016 #include "TBDataFormats/EcalTBObjects/interface/EcalTBHodoscopeRecInfo.h"
0017 #include "TBDataFormats/EcalTBObjects/interface/EcalTBTDCRecInfo.h"
0018 #include "TBDataFormats/EcalTBObjects/interface/EcalTBEventHeader.h"
0019
0020
0021
0022 #include "TFile.h"
0023 #include "TH1.h"
0024 #include "TH2.h"
0025 #include "TF1.h"
0026
0027 #include <iostream>
0028 #include <string>
0029 #include <stdexcept>
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043 EcalSimpleTBAnalyzer::EcalSimpleTBAnalyzer(const edm::ParameterSet& iConfig)
0044 : xtalInBeam_(0)
0045
0046 {
0047
0048 rootfile_ = iConfig.getUntrackedParameter<std::string>("rootfile", "ecalSimpleTBanalysis.root");
0049 digiCollection_ = iConfig.getParameter<std::string>("digiCollection");
0050 digiProducer_ = iConfig.getParameter<std::string>("digiProducer");
0051 hitCollection_ = iConfig.getParameter<std::string>("hitCollection");
0052 hitProducer_ = iConfig.getParameter<std::string>("hitProducer");
0053 hodoRecInfoCollection_ = iConfig.getParameter<std::string>("hodoRecInfoCollection");
0054 hodoRecInfoProducer_ = iConfig.getParameter<std::string>("hodoRecInfoProducer");
0055 tdcRecInfoCollection_ = iConfig.getParameter<std::string>("tdcRecInfoCollection");
0056 tdcRecInfoProducer_ = iConfig.getParameter<std::string>("tdcRecInfoProducer");
0057 eventHeaderCollection_ = iConfig.getParameter<std::string>("eventHeaderCollection");
0058 eventHeaderProducer_ = iConfig.getParameter<std::string>("eventHeaderProducer");
0059
0060 std::cout << "EcalSimpleTBAnalyzer: fetching hitCollection: " << hitCollection_.c_str() << " produced by "
0061 << hitProducer_.c_str() << std::endl;
0062 }
0063
0064
0065 EcalSimpleTBAnalyzer::~EcalSimpleTBAnalyzer()
0066
0067 {
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090 }
0091
0092
0093 void EcalSimpleTBAnalyzer::beginJob() {
0094
0095
0096
0097 h_ampltdc = new TH2F("h_ampltdc", "Max Amplitude vs TDC offset", 100, 0., 1., 1000, 0., 4000.);
0098
0099
0100 h_tableIsMoving = new TH1F("h_tableIsMoving", "TableIsMoving", 100000, 0., 100000.);
0101
0102 h_e1x1 = new TH1F("h_e1x1", "E1x1 energy", 1000, 0., 4000.);
0103 h_e3x3 = new TH1F("h_e3x3", "E3x3 energy", 1000, 0., 4000.);
0104 h_e5x5 = new TH1F("h_e5x5", "E5x5 energy", 1000, 0., 4000.);
0105
0106 h_e1x1_center = new TH1F("h_e1x1_center", "E1x1 energy", 1000, 0., 4000.);
0107 h_e3x3_center = new TH1F("h_e3x3_center", "E3x3 energy", 1000, 0., 4000.);
0108 h_e5x5_center = new TH1F("h_e5x5_center", "E5x5 energy", 1000, 0., 4000.);
0109
0110 h_e1e9 = new TH1F("h_e1e9", "E1/E9 ratio", 600, 0., 1.2);
0111 h_e1e25 = new TH1F("h_e1e25", "E1/E25 ratio", 600, 0., 1.2);
0112 h_e9e25 = new TH1F("h_e9e25", "E9/E25 ratio", 600, 0., 1.2);
0113
0114 h_bprofx = new TH1F("h_bprofx", "Beam Profile X", 100, -20., 20.);
0115 h_bprofy = new TH1F("h_bprofy", "Beam Profile Y", 100, -20., 20.);
0116
0117 h_qualx = new TH1F("h_qualx", "Beam Quality X", 5000, 0., 5.);
0118 h_qualy = new TH1F("h_qualy", "Beam Quality X", 5000, 0., 5.);
0119
0120 h_slopex = new TH1F("h_slopex", "Beam Slope X", 500, -5e-4, 5e-4);
0121 h_slopey = new TH1F("h_slopey", "Beam Slope Y", 500, -5e-4, 5e-4);
0122
0123 char hname[50];
0124 char htitle[50];
0125 for (unsigned int icry = 0; icry < 25; icry++) {
0126 sprintf(hname, "h_mapx_%d", icry);
0127 sprintf(htitle, "Max Amplitude vs X %d", icry);
0128 h_mapx[icry] = new TH2F(hname, htitle, 80, -20, 20, 1000, 0., 4000.);
0129 sprintf(hname, "h_mapy_%d", icry);
0130 sprintf(htitle, "Max Amplitude vs Y %d", icry);
0131 h_mapy[icry] = new TH2F(hname, htitle, 80, -20, 20, 1000, 0., 4000.);
0132 }
0133
0134 h_e1e9_mapx = new TH2F("h_e1e9_mapx", "E1/E9 vs X", 80, -20, 20, 600, 0., 1.2);
0135 h_e1e9_mapy = new TH2F("h_e1e9_mapy", "E1/E9 vs Y", 80, -20, 20, 600, 0., 1.2);
0136
0137 h_e1e25_mapx = new TH2F("h_e1e25_mapx", "E1/E25 vs X", 80, -20, 20, 600, 0., 1.2);
0138 h_e1e25_mapy = new TH2F("h_e1e25_mapy", "E1/E25 vs Y", 80, -20, 20, 600, 0., 1.2);
0139
0140 h_e9e25_mapx = new TH2F("h_e9e25_mapx", "E9/E25 vs X", 80, -20, 20, 600, 0., 1.2);
0141 h_e9e25_mapy = new TH2F("h_e9e25_mapy", "E9/E25 vs Y", 80, -20, 20, 600, 0., 1.2);
0142
0143 h_Shape_ = new TH2F("h_Shape_", "Xtal in Beam Shape", 250, 0, 10, 350, 0, 3500);
0144 }
0145
0146
0147 void EcalSimpleTBAnalyzer::endJob() {
0148
0149
0150 TFile f(rootfile_.c_str(), "RECREATE");
0151
0152
0153 h_ampltdc->Write();
0154
0155
0156 h_e1x1->Write();
0157 h_e3x3->Write();
0158 h_e5x5->Write();
0159
0160 h_e1x1_center->Write();
0161 h_e3x3_center->Write();
0162 h_e5x5_center->Write();
0163
0164 h_e1e9->Write();
0165 h_e1e25->Write();
0166 h_e9e25->Write();
0167
0168 h_bprofx->Write();
0169 h_bprofy->Write();
0170
0171 h_qualx->Write();
0172 h_qualy->Write();
0173
0174 h_slopex->Write();
0175 h_slopey->Write();
0176
0177 h_Shape_->Write();
0178
0179 for (unsigned int icry = 0; icry < 25; icry++) {
0180 h_mapx[icry]->Write();
0181 h_mapy[icry]->Write();
0182 }
0183
0184 h_e1e9_mapx->Write();
0185 h_e1e9_mapy->Write();
0186
0187 h_e1e25_mapx->Write();
0188 h_e1e25_mapy->Write();
0189
0190 h_e9e25_mapx->Write();
0191 h_e9e25_mapy->Write();
0192
0193 h_tableIsMoving->Write();
0194
0195 f.Close();
0196 }
0197
0198
0199
0200
0201
0202
0203 void EcalSimpleTBAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0204
0205
0206 using namespace edm;
0207 using namespace cms;
0208
0209 Handle<EBDigiCollection> pdigis;
0210 const EBDigiCollection* digis = nullptr;
0211
0212 iEvent.getByLabel(digiProducer_, digiCollection_, pdigis);
0213 if (pdigis.isValid()) {
0214 digis = pdigis.product();
0215
0216 } else {
0217 edm::LogError("EcalSimpleTBAnalyzerError") << "Error! can't get the product " << digiCollection_;
0218 }
0219
0220
0221 Handle<EBUncalibratedRecHitCollection> phits;
0222 const EBUncalibratedRecHitCollection* hits = nullptr;
0223
0224 iEvent.getByLabel(hitProducer_, hitCollection_, phits);
0225 if (phits.isValid()) {
0226 hits = phits.product();
0227
0228 } else {
0229 edm::LogError("EcalSimpleTBAnalyzerError") << "Error! can't get the product " << hitCollection_;
0230 }
0231
0232 Handle<EcalTBHodoscopeRecInfo> pHodo;
0233 const EcalTBHodoscopeRecInfo* recHodo = nullptr;
0234
0235 iEvent.getByLabel(hodoRecInfoProducer_, hodoRecInfoCollection_, pHodo);
0236 if (pHodo.isValid()) {
0237 recHodo = pHodo.product();
0238 } else {
0239 edm::LogError("EcalSimpleTBAnalyzerError") << "Error! can't get the product " << hodoRecInfoCollection_;
0240 }
0241
0242 Handle<EcalTBTDCRecInfo> pTDC;
0243 const EcalTBTDCRecInfo* recTDC = nullptr;
0244
0245 iEvent.getByLabel(tdcRecInfoProducer_, tdcRecInfoCollection_, pTDC);
0246 if (pTDC.isValid()) {
0247 recTDC = pTDC.product();
0248 } else {
0249 edm::LogError("EcalSimpleTBAnalyzerError") << "Error! can't get the product " << tdcRecInfoCollection_;
0250 }
0251
0252 Handle<EcalTBEventHeader> pEventHeader;
0253 const EcalTBEventHeader* evtHeader = nullptr;
0254
0255 iEvent.getByLabel(eventHeaderProducer_, pEventHeader);
0256 if (pEventHeader.isValid()) {
0257 evtHeader = pEventHeader.product();
0258 } else {
0259 edm::LogError("EcalSimpleTBAnalyzerError") << "Error! can't get the product " << eventHeaderProducer_;
0260 }
0261
0262 if (!hits)
0263 return;
0264
0265 if (!recTDC)
0266 return;
0267
0268 if (!recHodo)
0269 return;
0270
0271 if (!evtHeader)
0272 return;
0273
0274 if (hits->empty())
0275 return;
0276
0277 if (evtHeader->tableIsMoving())
0278 h_tableIsMoving->Fill(evtHeader->eventNumber());
0279
0280
0281 if (xtalInBeam_.null()) {
0282 xtalInBeam_ = EBDetId(1, evtHeader->crystalInBeam(), EBDetId::SMCRYSTALMODE);
0283 std::cout << "Xtal In Beam is " << xtalInBeam_.ic() << std::endl;
0284 } else if (xtalInBeam_ != EBDetId(1, evtHeader->crystalInBeam(), EBDetId::SMCRYSTALMODE))
0285 return;
0286
0287 if (evtHeader->tableIsMoving())
0288 return;
0289
0290
0291
0292
0293
0294
0295
0296
0297
0298
0299
0300
0301
0302
0303
0304
0305
0306
0307
0308
0309
0310
0311
0312
0313 EBDetId Xtals5x5[25];
0314 for (unsigned int icry = 0; icry < 25; icry++) {
0315 unsigned int row = icry / 5;
0316 unsigned int column = icry % 5;
0317 int ieta = xtalInBeam_.ieta() + column - 2;
0318 int iphi = xtalInBeam_.iphi() + row - 2;
0319 EBDetId tempId(ieta, iphi, EBDetId::ETAPHIMODE);
0320 if (tempId.ism() == 1)
0321 Xtals5x5[icry] = tempId;
0322 else
0323 Xtals5x5[icry] = EBDetId(0);
0324
0325 }
0326
0327 double samples_save[10];
0328 for (int i = 0; i < 10; ++i)
0329 samples_save[i] = 0.0;
0330
0331
0332 EBDigiCollection::const_iterator myDg = digis->find(xtalInBeam_);
0333 double eMax = 0.;
0334 if (myDg != digis->end()) {
0335 EBDataFrame myDigi = (*myDg);
0336 for (int sample = 0; sample < myDigi.size(); ++sample) {
0337 double analogSample = myDigi.sample(sample).adc();
0338 samples_save[sample] = analogSample;
0339
0340 if (eMax < analogSample) {
0341 eMax = analogSample;
0342 }
0343 }
0344
0345 }
0346
0347 for (int i = 0; i < 10; ++i) {
0348 h_Shape_->Fill(double(i) + recTDC->offset(), samples_save[i]);
0349 }
0350
0351 double amplitude[25];
0352
0353 double amplitude3x3 = 0;
0354 double amplitude5x5 = 0;
0355
0356 for (unsigned int icry = 0; icry < 25; icry++) {
0357 if (!Xtals5x5[icry].null()) {
0358 amplitude[icry] = (hits->find(Xtals5x5[icry]))->amplitude();
0359 amplitude5x5 += amplitude[icry];
0360
0361 if (icry == 6 || icry == 7 || icry == 8 || icry == 11 || icry == 12 || icry == 13 || icry == 16 || icry == 17 ||
0362 icry == 18) {
0363 amplitude3x3 += amplitude[icry];
0364 }
0365 }
0366 }
0367
0368 h_e1x1->Fill(amplitude[12]);
0369 h_e3x3->Fill(amplitude3x3);
0370 h_e5x5->Fill(amplitude5x5);
0371
0372 h_e1e9->Fill(amplitude[12] / amplitude3x3);
0373 h_e1e25->Fill(amplitude[12] / amplitude5x5);
0374 h_e9e25->Fill(amplitude3x3 / amplitude5x5);
0375
0376 if (recTDC)
0377 h_ampltdc->Fill(recTDC->offset(), amplitude[12]);
0378
0379 if (recHodo) {
0380 float x = recHodo->posX();
0381 float y = recHodo->posY();
0382 float xslope = recHodo->slopeX();
0383 float yslope = recHodo->slopeY();
0384 float xqual = recHodo->qualX();
0385 float yqual = recHodo->qualY();
0386
0387
0388 h_bprofx->Fill(x);
0389 h_bprofy->Fill(y);
0390 h_qualx->Fill(xqual);
0391 h_qualy->Fill(yqual);
0392 h_slopex->Fill(xslope);
0393 h_slopey->Fill(yslope);
0394
0395
0396
0397 if (fabs(x + 2.5) < 2.5 && fabs(y + 0.5) < 2.5) {
0398 h_e1x1_center->Fill(amplitude[12]);
0399 h_e3x3_center->Fill(amplitude3x3);
0400 h_e5x5_center->Fill(amplitude5x5);
0401
0402 h_e1e9->Fill(amplitude[12] / amplitude3x3);
0403 h_e1e25->Fill(amplitude[12] / amplitude5x5);
0404 h_e9e25->Fill(amplitude3x3 / amplitude5x5);
0405 }
0406
0407 for (unsigned int icry = 0; icry < 25; icry++) {
0408 h_mapx[icry]->Fill(x, amplitude[icry]);
0409 h_mapy[icry]->Fill(y, amplitude[icry]);
0410 }
0411
0412 h_e1e9_mapx->Fill(x, amplitude[12] / amplitude3x3);
0413 h_e1e9_mapy->Fill(y, amplitude[12] / amplitude3x3);
0414
0415 h_e1e25_mapx->Fill(x, amplitude[12] / amplitude5x5);
0416 h_e1e25_mapy->Fill(y, amplitude[12] / amplitude5x5);
0417
0418 h_e9e25_mapx->Fill(x, amplitude3x3 / amplitude5x5);
0419 h_e9e25_mapy->Fill(y, amplitude3x3 / amplitude5x5);
0420 }
0421 }