File indexing completed on 2021-12-08 08:15:44
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 #include <memory>
0021 #include <iostream>
0022 #include <fstream>
0023
0024
0025 #include "FWCore/Framework/interface/Frameworkfwd.h"
0026 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0027
0028 #include "FWCore/Framework/interface/Event.h"
0029 #include "FWCore/Framework/interface/MakerMacros.h"
0030
0031 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0032 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0033
0034 #include "DataFormats/EcalRecHit/interface/EcalUncalibratedRecHit.h"
0035 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0036 #include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
0037 #include "DataFormats/EcalDigi/interface/EBDataFrame.h"
0038 #include "DataFormats/EcalRawData/interface/EcalRawDataCollections.h"
0039
0040 #include "TFile.h"
0041 #include "TH1F.h"
0042 #include "TH2F.h"
0043 #include "TROOT.h"
0044 #include "TF1.h"
0045 #include "TNtuple.h"
0046 #include "TDirectory.h"
0047
0048
0049
0050
0051
0052 class EcalLaserAnalyzerYousi : public edm::one::EDAnalyzer<> {
0053 public:
0054 explicit EcalLaserAnalyzerYousi(const edm::ParameterSet &);
0055 ~EcalLaserAnalyzerYousi() override;
0056
0057 private:
0058 void beginJob() override;
0059 void analyze(const edm::Event &, const edm::EventSetup &) override;
0060 void endJob() override;
0061
0062
0063
0064 TH1F *C_APD[1700];
0065 TH1F *C_APDPN[1700];
0066 TH1F *C_PN[1700];
0067 TH1F *C_J[1700];
0068 TH2F *C_APDPN_J[1700];
0069
0070 TH1F *peakAPD[2];
0071 TH1F *peakAPDPN[2];
0072 TH1F *APD_LM[9];
0073 TH1F *APDPN_LM[9];
0074 TH2F *APDPN_J_LM[9];
0075 TH2F *APDPN_J_H[2];
0076
0077
0078 TH2F *APD;
0079 TH2F *APD_RMS;
0080 TH2F *APDPN;
0081 TH2F *APDPN_RMS;
0082 TH2F *PN;
0083 TH2F *APDPN_J;
0084 TH2F *APDPN_C;
0085
0086 TH1F *FitHist;
0087 TH2F *Count;
0088
0089 TFile *fPN;
0090 TFile *fAPD;
0091 TFile *fROOT;
0092
0093 TNtuple *C_Tree[1700];
0094
0095
0096
0097 std::string hitCollection_;
0098 std::string hitProducer_;
0099
0100
0101 std::string outFileName_;
0102 std::string SM_;
0103 std::string Run_;
0104 std::string digiProducer_;
0105 std::string PNdigiCollection_;
0106 };
0107
0108
0109
0110
0111
0112
0113
0114
0115
0116
0117
0118
0119 EcalLaserAnalyzerYousi::EcalLaserAnalyzerYousi(const edm::ParameterSet &iConfig) {
0120
0121
0122
0123
0124 hitCollection_ = iConfig.getUntrackedParameter<std::string>("hitCollection");
0125 hitProducer_ = iConfig.getUntrackedParameter<std::string>("hitProducer");
0126
0127
0128 outFileName_ = iConfig.getUntrackedParameter<std::string>("outFileName");
0129 SM_ = iConfig.getUntrackedParameter<std::string>("SM");
0130 Run_ = iConfig.getUntrackedParameter<std::string>("Run");
0131 digiProducer_ = iConfig.getUntrackedParameter<std::string>("digiProducer");
0132 PNdigiCollection_ = iConfig.getUntrackedParameter<std::string>("PNdigiCollection");
0133 }
0134
0135 EcalLaserAnalyzerYousi::~EcalLaserAnalyzerYousi() {
0136
0137
0138 }
0139
0140
0141
0142
0143
0144
0145 void EcalLaserAnalyzerYousi::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0146 using namespace edm;
0147
0148
0149
0150 edm::Handle<EcalRawDataCollection> DCCHeaders;
0151 iEvent.getByLabel(digiProducer_, DCCHeaders);
0152
0153 EcalDCCHeaderBlock::EcalDCCEventSettings settings = DCCHeaders->begin()->getEventSettings();
0154
0155 int wavelength = settings.wavelength;
0156
0157
0158
0159 if (wavelength != 0) {
0160 return;
0161 }
0162
0163 edm::Handle<EBUncalibratedRecHitCollection> hits;
0164
0165 try {
0166 iEvent.getByLabel(hitProducer_, hitCollection_, hits);
0167
0168 } catch (std::exception &ex) {
0169 LogError("EcalLaserAnalyzerYousi") << "Cannot get product: EBRecHitCollection from: " << hitCollection_
0170 << " - returning.\n\n";
0171
0172 }
0173
0174 edm::Handle<EcalPnDiodeDigiCollection> pndigis;
0175
0176 try {
0177
0178 iEvent.getByLabel(digiProducer_, PNdigiCollection_, pndigis);
0179
0180 } catch (std::exception &ex) {
0181 LogError("EcalLaserAnalyzerYousi") << "Cannot get product: EBdigiCollection from: "
0182 << "getbytype"
0183 << " - returning.\n\n";
0184
0185 }
0186
0187 Float_t PN_amp[5];
0188
0189
0190 for (int j = 0; j < 5; ++j) {
0191 PN_amp[j] = 0;
0192 for (int z = 0; z < 2; ++z) {
0193 FitHist->Reset();
0194 TF1 peakFit("peakFit", "[0] +[1]*x +[2]*x^2", 30, 50);
0195 TF1 pedFit("pedFit", "[0]", 0, 5);
0196
0197 for (int k = 0; k < 50; k++) {
0198 FitHist->SetBinContent(k, (*pndigis)[j + z * 5].sample(k).adc());
0199 }
0200 pedFit.SetParameter(0, 750);
0201 FitHist->Fit(&pedFit, "RQI");
0202 Float_t ped = pedFit.GetParameter(0);
0203
0204 Int_t maxbin = FitHist->GetMaximumBin();
0205 peakFit.SetRange(FitHist->GetBinCenter(maxbin) - 4 * FitHist->GetBinWidth(maxbin),
0206 FitHist->GetBinCenter(maxbin) + 4 * FitHist->GetBinWidth(maxbin));
0207 peakFit.SetParameters(750, 4, -.05);
0208 FitHist->Fit(&peakFit, "RQI");
0209 Float_t max = peakFit.Eval(-peakFit.GetParameter(1) / (2 * peakFit.GetParameter(2)));
0210 if (ped != max) {
0211 PN_amp[j] = PN_amp[j] + max - ped;
0212 } else {
0213 PN_amp[j] = PN_amp[j] + max;
0214 }
0215
0216 }
0217 PN_amp[j] = PN_amp[j] / 2.0;
0218 }
0219
0220
0221
0222
0223
0224
0225
0226
0227
0228 Float_t fTree[7];
0229
0230
0231 EBDetId ID;
0232 Float_t theAPD;
0233 Float_t thePN;
0234 Float_t Jitter;
0235 Float_t Chi2;
0236 Int_t CN = hits->size();
0237
0238 for (int j = 0; j < CN; j++) {
0239 ID = (*hits)[j].id();
0240 theAPD = (*hits)[j].amplitude();
0241 Jitter = (*hits)[j].jitter();
0242 Chi2 = (*hits)[j].chi2();
0243 thePN = PN_amp[(ID.ic() + 299) / 400];
0244
0245
0246
0247
0248 C_APD[ID.ic() - 1]->Fill(theAPD);
0249 C_APDPN[ID.ic() - 1]->Fill(theAPD / thePN);
0250 C_PN[ID.ic() - 1]->Fill(thePN);
0251 C_J[ID.ic() - 1]->Fill(Jitter);
0252 C_APDPN_J[ID.ic() - 1]->Fill(Jitter, theAPD / thePN);
0253 APDPN_J->Fill(Jitter, theAPD / thePN);
0254 APDPN_C->Fill(Chi2, theAPD / thePN);
0255 fTree[0] = theAPD;
0256 fTree[1] = thePN;
0257 fTree[2] = theAPD / thePN;
0258 fTree[3] = Jitter;
0259 fTree[4] = Chi2;
0260 fTree[5] = (*hits)[j].pedestal();
0261 fTree[6] = iEvent.id().event();
0262 C_Tree[ID.ic() - 1]->Fill(fTree);
0263 if (((ID.ic() - 1) % 20 > 9) || ((ID.ic() - 1) < 100)) {
0264 peakAPD[0]->Fill(theAPD);
0265 peakAPDPN[0]->Fill(theAPD / thePN);
0266 APDPN_J_H[0]->Fill(Jitter, theAPD / thePN);
0267 } else {
0268 peakAPD[1]->Fill(theAPD);
0269 peakAPDPN[1]->Fill(theAPD / thePN);
0270 APDPN_J_H[1]->Fill(Jitter, theAPD / thePN);
0271 }
0272 if ((ID.ic() - 1) < 100) {
0273 APD_LM[0]->Fill(theAPD);
0274 APDPN_LM[0]->Fill(theAPD / thePN);
0275 APDPN_J_LM[0]->Fill(Jitter, theAPD / thePN);
0276 } else {
0277 Int_t index;
0278 if (((ID.ic() - 1) % 20) < 10) {
0279 index = ((ID.ic() - 101) / 400) * 2 + 1;
0280 APD_LM[index]->Fill(theAPD);
0281 APDPN_LM[index]->Fill(theAPD / thePN);
0282 APDPN_J_LM[index]->Fill(Jitter, theAPD / thePN);
0283 } else {
0284 index = ((ID.ic() - 101) / 400) * 2 + 2;
0285 APD_LM[index]->Fill(theAPD);
0286 APDPN_LM[index]->Fill(theAPD / thePN);
0287 APDPN_J_LM[index]->Fill(Jitter, theAPD / thePN);
0288 }
0289 }
0290 }
0291
0292
0293
0294
0295
0296
0297
0298
0299
0300
0301 }
0302
0303
0304 void EcalLaserAnalyzerYousi::beginJob() {
0305 edm::LogInfo("EcalLaserAnalyzerYousi") << "running laser analyzer \n\n";
0306
0307 fROOT = new TFile(outFileName_.c_str(), "RECREATE");
0308 fROOT->cd();
0309
0310
0311 APD = new TH2F("APD", "APD", 85, 0., 85., 20, 0., 20.);
0312 APD_RMS = new TH2F("APD_RMS", "APD_RMS", 85, 0., 85., 20, 0., 20.);
0313 APDPN = new TH2F("APDPN", "APDPN", 85, 0., 85., 20, 0., 20.);
0314 APDPN_RMS = new TH2F("APDPN_RMS", "APDPN_RMS", 85, 0., 85., 20, 0., 20.);
0315 PN = new TH2F("PN", "PN", 85, 0., 85., 20, 0., 20.);
0316 APDPN_J = new TH2F("JittervAPDPN", "JittervAPDPN", 250, 3., 7., 250, 1., 2.);
0317 APDPN_C = new TH2F("Chi2vAPDPN", "Chi2vAPDPN", 250, 0., 50., 250, 0., 5.0);
0318 FitHist = new TH1F("FitHist", "FitHist", 50, 0, 50);
0319 Count = new TH2F("Count", "Count", 85, 0., 1., 20, 0., 1.);
0320
0321 for (int i = 0; i < 1700; i++) {
0322 std::ostringstream name_1;
0323 std::ostringstream name_2;
0324 std::ostringstream name_3;
0325 std::ostringstream name_4;
0326 std::ostringstream name_5;
0327 name_1 << "C_APD_" << i + 1;
0328 name_2 << "C_APDPN_" << i + 1;
0329 name_3 << "C_PN_" << i + 1;
0330 name_4 << "C_J_" << i + 1;
0331 name_5 << "C_APDPN_J_" << i + 1;
0332 C_APD[i] = new TH1F(name_1.str().c_str(), name_1.str().c_str(), 2500, 0., 5000.);
0333 C_APDPN[i] = new TH1F(name_2.str().c_str(), name_2.str().c_str(), 20000, 0., 25.);
0334 C_PN[i] = new TH1F(name_3.str().c_str(), name_3.str().c_str(), 1000, 0., 4000.);
0335 C_J[i] = new TH1F(name_4.str().c_str(), name_4.str().c_str(), 250, 3.0, 7.);
0336 C_APDPN_J[i] = new TH2F(name_5.str().c_str(), name_5.str().c_str(), 250, 3.0, 6., 250, 1., 2.2);
0337 }
0338
0339 for (int i = 0; i < 2; i++) {
0340 std::ostringstream aname_1;
0341 std::ostringstream aname_2;
0342 std::ostringstream aname_3;
0343 aname_1 << "peakAPD_" << i;
0344 aname_2 << "peakAPDPN_" << i;
0345 aname_3 << "JittervAPDPN_Half_" << i;
0346 peakAPD[i] = new TH1F(aname_1.str().c_str(), aname_1.str().c_str(), 1000, 0., 5000.);
0347 peakAPDPN[i] = new TH1F(aname_2.str().c_str(), aname_2.str().c_str(), 1000, 0., 8.);
0348 APDPN_J_H[i] = new TH2F(aname_3.str().c_str(), aname_3.str().c_str(), 250, 3., 7., 250, 1., 2.2);
0349 }
0350
0351 for (int i = 0; i < 9; i++) {
0352 std::ostringstream bname_1;
0353 std::ostringstream bname_2;
0354 std::ostringstream bname_3;
0355 bname_1 << "APD_LM_" << i;
0356 bname_2 << "APDPN_LM_" << i;
0357 bname_3 << "APDPN_J_LM_" << i;
0358 APD_LM[i] = new TH1F(bname_1.str().c_str(), bname_1.str().c_str(), 500, 0., 5000.);
0359 APDPN_LM[i] = new TH1F(bname_2.str().c_str(), bname_2.str().c_str(), 500, 0., 8.);
0360 APDPN_J_LM[i] = new TH2F(bname_3.str().c_str(), bname_3.str().c_str(), 250, 3., 7., 250, 1., 2.2);
0361 }
0362
0363
0364
0365
0366
0367
0368 std::ostringstream varlist;
0369 varlist << "APD:PN:APDPN:Jitter:Chi2:ped:EVT";
0370 for (int i = 0; i < 1700; i++) {
0371 std::ostringstream name;
0372 name << "C_Tree_" << i + 1;
0373 C_Tree[i] = (TNtuple *)new TNtuple(name.str().c_str(), name.str().c_str(), varlist.str().c_str());
0374 }
0375 }
0376
0377
0378 void EcalLaserAnalyzerYousi::endJob() {
0379
0380 TFile *fROOT = (TFile *)new TFile(outFileName_.c_str(), "RECREATE");
0381
0382
0383 TDirectory *DIR;
0384
0385 DIR = fROOT->mkdir(Run_.c_str());
0386
0387 DIR->cd();
0388 for (int j = 0; j < 1700; j++) {
0389 Float_t min_r, max_r;
0390 Float_t RMS, Sigma, K;
0391 Int_t iCount;
0392 TF1 *gs3;
0393
0394 RMS = C_APD[j]->GetRMS();
0395 APD_RMS->SetBinContent(85 - (j / 20), 20 - (j % 20), RMS);
0396 Sigma = 999999;
0397 K = 2.5;
0398 iCount = 0;
0399 while (Sigma > RMS) {
0400 min_r = C_APD[j]->GetBinCenter(C_APD[j]->GetMaximumBin()) - K * RMS;
0401 max_r = C_APD[j]->GetBinCenter(C_APD[j]->GetMaximumBin()) + K * RMS;
0402 TF1 gs1("gs1", "gaus", min_r, max_r);
0403 C_APD[j]->Fit(&gs1, "RQI");
0404 Sigma = gs1.GetParameter(2);
0405 K = K * 1.5;
0406 iCount++;
0407 if (iCount > 2) {
0408 C_APD[j]->Fit("gaus", "QI");
0409 break;
0410 }
0411 }
0412
0413 RMS = C_APDPN[j]->GetRMS();
0414 APDPN_RMS->SetBinContent(85 - (j / 20), 20 - (j % 20), RMS);
0415 Sigma = 999999;
0416 K = 2.5;
0417 iCount = 0;
0418 while (Sigma > RMS) {
0419 min_r = C_APDPN[j]->GetBinCenter(C_APDPN[j]->GetMaximumBin()) - K * RMS;
0420 max_r = C_APDPN[j]->GetBinCenter(C_APDPN[j]->GetMaximumBin()) + K * RMS;
0421 TF1 gs2("gs2", "gaus", min_r, max_r);
0422 C_APDPN[j]->Fit(&gs2, "RQI");
0423 Sigma = gs2.GetParameter(2);
0424 K = K * 1.5;
0425 iCount++;
0426 if (iCount > 2) {
0427 C_APDPN[j]->Fit("gaus", "QI");
0428 break;
0429 }
0430 }
0431
0432 TF1 *newgs1;
0433 TF1 *newgs2;
0434
0435 C_PN[j]->Fit("gaus", "Q");
0436 C_APD[j]->Fit("gaus", "QI");
0437 C_APDPN[j]->Fit("gaus", "QI");
0438 C_APD[j]->Write("", TObject::kOverwrite);
0439 C_APDPN[j]->Write("", TObject::kOverwrite);
0440 C_PN[j]->Write("", TObject::kOverwrite);
0441 C_J[j]->Write("", TObject::kOverwrite);
0442 C_APDPN_J[j]->Write("", TObject::kOverwrite);
0443 newgs1 = C_APD[j]->GetFunction("gaus");
0444 newgs2 = C_APDPN[j]->GetFunction("gaus");
0445 gs3 = C_PN[j]->GetFunction("gaus");
0446 Float_t theAPD = newgs1->GetParameter(1);
0447 APD->SetBinContent(85 - (j / 20), 20 - (j % 20), theAPD);
0448 Float_t theAPDPN = newgs2->GetParameter(1);
0449 APDPN->SetBinContent(85 - (j / 20), 20 - (j % 20), theAPDPN);
0450 Float_t thePN = gs3->GetParameter(1);
0451
0452 PN->SetBinContent(85 - (j / 20), 20 - (j % 20), thePN);
0453 C_Tree[j]->Write("", TObject::kOverwrite);
0454 }
0455
0456 for (int i = 0; i < 9; i++) {
0457 APD_LM[i]->Write("", TObject::kOverwrite);
0458 APDPN_LM[i]->Write("", TObject::kOverwrite);
0459 APDPN_J_LM[i]->Write("", TObject::kOverwrite);
0460 }
0461
0462
0463 APD->Write("", TObject::kOverwrite);
0464 APD_RMS->Write("", TObject::kOverwrite);
0465 APDPN_RMS->Write("", TObject::kOverwrite);
0466 APDPN->Write("", TObject::kOverwrite);
0467 APDPN_J->Write("", TObject::kOverwrite);
0468 APDPN_C->Write("", TObject::kOverwrite);
0469 PN->Write("", TObject::kOverwrite);
0470 peakAPD[0]->Write("", TObject::kOverwrite);
0471 peakAPD[1]->Write("", TObject::kOverwrite);
0472 peakAPDPN[0]->Write("", TObject::kOverwrite);
0473 peakAPDPN[1]->Write("", TObject::kOverwrite);
0474 APDPN_J_H[0]->Write("", TObject::kOverwrite);
0475 APDPN_J_H[1]->Write("", TObject::kOverwrite);
0476
0477
0478
0479
0480
0481
0482
0483 fROOT->Write();
0484
0485 }
0486
0487
0488 DEFINE_FWK_MODULE(EcalLaserAnalyzerYousi);