File indexing completed on 2021-12-08 08:15:42
0001
0002
0003
0004
0005
0006
0007
0008 #include "TAxis.h"
0009 #include "TH1.h"
0010 #include "TProfile.h"
0011 #include "TTree.h"
0012 #include "TChain.h"
0013 #include "TFile.h"
0014 #include "TMath.h"
0015
0016 #include "CalibCalorimetry/EcalLaserAnalyzer/plugins/EcalLaserAnalyzer.h"
0017
0018 #include <sstream>
0019 #include <fstream>
0020 #include <iomanip>
0021
0022 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0023
0024 #include "FWCore/Framework/interface/EventSetup.h"
0025
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0029 #include "DataFormats/EcalDetId/interface/EcalElectronicsId.h"
0030 #include "DataFormats/EcalDetId/interface/EcalDetIdCollections.h"
0031
0032 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/TShapeAnalysis.h"
0033 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/TPNFit.h"
0034 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/TMom.h"
0035 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/TAPD.h"
0036 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/TAPDPulse.h"
0037 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/TPN.h"
0038 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/TPNPulse.h"
0039 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/TPNCor.h"
0040 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/TMem.h"
0041 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/PulseFitWithFunction.h"
0042 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/ME.h"
0043 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/MEGeom.h"
0044
0045
0046 EcalLaserAnalyzer::EcalLaserAnalyzer(const edm::ParameterSet& iConfig)
0047
0048 : iEvent(0),
0049 eventHeaderCollection_(iConfig.getParameter<std::string>("eventHeaderCollection")),
0050 eventHeaderProducer_(iConfig.getParameter<std::string>("eventHeaderProducer")),
0051 digiCollection_(iConfig.getParameter<std::string>("digiCollection")),
0052 digiProducer_(iConfig.getParameter<std::string>("digiProducer")),
0053 digiPNCollection_(iConfig.getParameter<std::string>("digiPNCollection")),
0054 rawDataToken_(consumes<EcalRawDataCollection>(edm::InputTag(eventHeaderProducer_, eventHeaderCollection_))),
0055 pnDiodeDigiToken_(consumes<EcalPnDiodeDigiCollection>(edm::InputTag(digiProducer_))),
0056 mappingToken_(esConsumes()),
0057
0058 _nsamples(iConfig.getUntrackedParameter<unsigned int>("nSamples", 10)),
0059 _presample(iConfig.getUntrackedParameter<unsigned int>("nPresamples", 2)),
0060 _firstsample(iConfig.getUntrackedParameter<unsigned int>("firstSample", 1)),
0061 _lastsample(iConfig.getUntrackedParameter<unsigned int>("lastSample", 2)),
0062 _nsamplesPN(iConfig.getUntrackedParameter<unsigned int>("nSamplesPN", 50)),
0063 _presamplePN(iConfig.getUntrackedParameter<unsigned int>("nPresamplesPN", 6)),
0064 _firstsamplePN(iConfig.getUntrackedParameter<unsigned int>("firstSamplePN", 7)),
0065 _lastsamplePN(iConfig.getUntrackedParameter<unsigned int>("lastSamplePN", 8)),
0066 _timingcutlow(iConfig.getUntrackedParameter<unsigned int>("timingCutLow", 2)),
0067 _timingcuthigh(iConfig.getUntrackedParameter<unsigned int>("timingCutHigh", 9)),
0068 _timingquallow(iConfig.getUntrackedParameter<unsigned int>("timingQualLow", 3)),
0069 _timingqualhigh(iConfig.getUntrackedParameter<unsigned int>("timingQualHigh", 8)),
0070 _ratiomincutlow(iConfig.getUntrackedParameter<double>("ratioMinCutLow", 0.4)),
0071 _ratiomincuthigh(iConfig.getUntrackedParameter<double>("ratioMinCutHigh", 0.95)),
0072 _ratiomaxcutlow(iConfig.getUntrackedParameter<double>("ratioMaxCutLow", 0.8)),
0073 _presamplecut(iConfig.getUntrackedParameter<double>("presampleCut", 5.0)),
0074 _niter(iConfig.getUntrackedParameter<unsigned int>("nIter", 3)),
0075 _fitab(iConfig.getUntrackedParameter<bool>("fitAB", false)),
0076 _alpha(iConfig.getUntrackedParameter<double>("alpha", 1.5076494)),
0077 _beta(iConfig.getUntrackedParameter<double>("beta", 1.5136036)),
0078 _nevtmax(iConfig.getUntrackedParameter<unsigned int>("nEvtMax", 200)),
0079 _noise(iConfig.getUntrackedParameter<double>("noise", 2.0)),
0080 _chi2cut(iConfig.getUntrackedParameter<double>("chi2cut", 10.0)),
0081 _ecalPart(iConfig.getUntrackedParameter<std::string>("ecalPart", "EB")),
0082 _docorpn(iConfig.getUntrackedParameter<bool>("doCorPN", false)),
0083 _fedid(iConfig.getUntrackedParameter<int>("fedID", -999)),
0084 _saveallevents(iConfig.getUntrackedParameter<bool>("saveAllEvents", false)),
0085 _qualpercent(iConfig.getUntrackedParameter<double>("qualPercent", 0.2)),
0086 _debug(iConfig.getUntrackedParameter<int>("debug", 0)),
0087 resdir_(iConfig.getUntrackedParameter<std::string>("resDir")),
0088 pncorfile_(iConfig.getUntrackedParameter<std::string>("pnCorFile")),
0089 nCrys(NCRYSEB),
0090 nPNPerMod(NPNPERMOD),
0091 nMod(NMODEE),
0092 nSides(NSIDES),
0093 runType(-1),
0094 runNum(0),
0095 fedID(-1),
0096 dccID(-1),
0097 side(2),
0098 lightside(2),
0099 iZ(1),
0100 phi(-1),
0101 eta(-1),
0102 event(0),
0103 color(-1),
0104 pn0(0),
0105 pn1(0),
0106 apdAmpl(0),
0107 apdAmplA(0),
0108 apdAmplB(0),
0109 apdTime(0),
0110 pnAmpl(0),
0111 pnID(-1),
0112 moduleID(-1),
0113 channelIteratorEE(0)
0114
0115
0116
0117 {
0118 if (_ecalPart == "EB") {
0119 ebDigiToken_ = consumes<EBDigiCollection>(edm::InputTag(digiProducer_, digiCollection_));
0120 } else if (_ecalPart == "EE") {
0121 eeDigiToken_ = consumes<EEDigiCollection>(edm::InputTag(digiProducer_, digiCollection_));
0122 }
0123
0124
0125
0126 if (_ecalPart == "EB") {
0127 nCrys = NCRYSEB;
0128 } else {
0129 nCrys = NCRYSEE;
0130 }
0131 iZ = 1;
0132 if (_fedid <= 609)
0133 iZ = -1;
0134 modules = ME::lmmodFromDcc(_fedid);
0135 nMod = modules.size();
0136 nRefChan = NREFCHAN;
0137
0138 for (unsigned int j = 0; j < nCrys; j++) {
0139 iEta[j] = -1;
0140 iPhi[j] = -1;
0141 iModule[j] = 10;
0142 iTowerID[j] = -1;
0143 iChannelID[j] = -1;
0144 idccID[j] = -1;
0145 iside[j] = -1;
0146 wasTimingOK[j] = true;
0147 wasGainOK[j] = true;
0148 }
0149
0150 for (unsigned int j = 0; j < nMod; j++) {
0151 int ii = modules[j];
0152 firstChanMod[ii - 1] = 0;
0153 isFirstChanModFilled[ii - 1] = 0;
0154 }
0155
0156
0157
0158 isGainOK = true;
0159 isTimingOK = true;
0160
0161
0162
0163 pnCorrector = new TPNCor(pncorfile_);
0164
0165
0166
0167 APDPulse = new TAPDPulse(_nsamples,
0168 _presample,
0169 _firstsample,
0170 _lastsample,
0171 _timingcutlow,
0172 _timingcuthigh,
0173 _timingquallow,
0174 _timingqualhigh,
0175 _ratiomincutlow,
0176 _ratiomincuthigh,
0177 _ratiomaxcutlow);
0178 PNPulse = new TPNPulse(_nsamplesPN, _presamplePN);
0179
0180
0181
0182 Mem = new TMem(_fedid);
0183
0184
0185
0186 Delta01 = new TMom();
0187 Delta12 = new TMom();
0188 }
0189
0190
0191 EcalLaserAnalyzer::~EcalLaserAnalyzer() {
0192
0193
0194
0195
0196 }
0197
0198
0199 void EcalLaserAnalyzer::beginJob() {
0200
0201
0202
0203
0204
0205 ADCfile = resdir_;
0206 ADCfile += "/APDSamplesLaser.root";
0207
0208 APDfile = resdir_;
0209 APDfile += "/APDPNLaserAllEvents.root";
0210
0211 ADCFile = new TFile(ADCfile.c_str(), "RECREATE");
0212
0213 for (unsigned int i = 0; i < nCrys; i++) {
0214 std::stringstream name;
0215 name << "ADCTree" << i + 1;
0216 ADCtrees[i] = new TTree(name.str().c_str(), name.str().c_str());
0217
0218 ADCtrees[i]->Branch("ieta", &eta, "eta/I");
0219 ADCtrees[i]->Branch("iphi", &phi, "phi/I");
0220 ADCtrees[i]->Branch("side", &side, "side/I");
0221 ADCtrees[i]->Branch("dccID", &dccID, "dccID/I");
0222 ADCtrees[i]->Branch("towerID", &towerID, "towerID/I");
0223 ADCtrees[i]->Branch("channelID", &channelID, "channelID/I");
0224 ADCtrees[i]->Branch("event", &event, "event/I");
0225 ADCtrees[i]->Branch("color", &color, "color/I");
0226 ADCtrees[i]->Branch("adc", &adc, "adc[10]/D");
0227 ADCtrees[i]->Branch("pn0", &pn0, "pn0/D");
0228 ADCtrees[i]->Branch("pn1", &pn1, "pn1/D");
0229
0230 ADCtrees[i]->SetBranchAddress("ieta", &eta);
0231 ADCtrees[i]->SetBranchAddress("iphi", &phi);
0232 ADCtrees[i]->SetBranchAddress("side", &side);
0233 ADCtrees[i]->SetBranchAddress("dccID", &dccID);
0234 ADCtrees[i]->SetBranchAddress("towerID", &towerID);
0235 ADCtrees[i]->SetBranchAddress("channelID", &channelID);
0236 ADCtrees[i]->SetBranchAddress("event", &event);
0237 ADCtrees[i]->SetBranchAddress("color", &color);
0238 ADCtrees[i]->SetBranchAddress("adc", adc);
0239 ADCtrees[i]->SetBranchAddress("pn0", &pn0);
0240 ADCtrees[i]->SetBranchAddress("pn1", &pn1);
0241
0242 nevtAB[i] = 0;
0243 }
0244
0245
0246
0247
0248
0249
0250 doesABTreeExist = true;
0251
0252 std::stringstream nameabinitfile;
0253 nameabinitfile << resdir_ << "/ABInit.root";
0254 alphainitfile = nameabinitfile.str();
0255
0256 std::stringstream nameabfile;
0257 nameabfile << resdir_ << "/AB.root";
0258 alphafile = nameabfile.str();
0259
0260 FILE* test;
0261 if (_fitab)
0262 test = fopen(alphainitfile.c_str(), "r");
0263 else
0264 test = fopen(alphafile.c_str(), "r");
0265 if (test == nullptr) {
0266 doesABTreeExist = false;
0267 _fitab = true;
0268 };
0269 delete test;
0270
0271 TFile* fAB = nullptr;
0272 TTree* ABInit = nullptr;
0273 if (doesABTreeExist) {
0274 fAB = new TFile(nameabinitfile.str().c_str());
0275 }
0276 if (doesABTreeExist && fAB) {
0277 ABInit = (TTree*)fAB->Get("ABCol0");
0278 }
0279
0280
0281
0282 if (doesABTreeExist && fAB && ABInit && ABInit->GetEntries() != 0) {
0283 shapana = new TShapeAnalysis(ABInit, _alpha, _beta, 5.5, 1.0);
0284 doesABTreeExist = true;
0285 } else {
0286 shapana = new TShapeAnalysis(_alpha, _beta, 5.5, 1.0);
0287 doesABTreeExist = false;
0288 _fitab = true;
0289 }
0290 shapana->set_const(_nsamples, _firstsample, _lastsample, _presample, _nevtmax, _noise, _chi2cut);
0291
0292 if (doesABTreeExist && fAB)
0293 fAB->Close();
0294
0295
0296
0297 std::stringstream nameapdfile;
0298 nameapdfile << resdir_ << "/APDPN_LASER.root";
0299 resfile = nameapdfile.str();
0300
0301
0302
0303 laserEvents = 0;
0304 }
0305
0306
0307 void EcalLaserAnalyzer::analyze(const edm::Event& e, const edm::EventSetup& c) {
0308
0309
0310 ++iEvent;
0311
0312
0313 edm::Handle<EcalRawDataCollection> pDCCHeader;
0314 const EcalRawDataCollection* DCCHeader = nullptr;
0315 e.getByToken(rawDataToken_, pDCCHeader);
0316 if (!pDCCHeader.isValid()) {
0317 edm::LogError("nodata") << "Error! can't get the product retrieving DCC header" << eventHeaderCollection_.c_str()
0318 << " " << eventHeaderProducer_.c_str();
0319 } else {
0320 DCCHeader = pDCCHeader.product();
0321 }
0322
0323
0324 edm::Handle<EBDigiCollection> pEBDigi;
0325 const EBDigiCollection* EBDigi = nullptr;
0326 edm::Handle<EEDigiCollection> pEEDigi;
0327 const EEDigiCollection* EEDigi = nullptr;
0328 if (_ecalPart == "EB") {
0329 e.getByToken(ebDigiToken_, pEBDigi);
0330 if (!pEBDigi.isValid()) {
0331 edm::LogError("nodata") << "Error! can't get the product retrieving EB crystal data " << digiCollection_.c_str();
0332 } else {
0333 EBDigi = pEBDigi.product();
0334 }
0335 } else if (_ecalPart == "EE") {
0336 e.getByToken(eeDigiToken_, pEEDigi);
0337 if (!pEEDigi.isValid()) {
0338 edm::LogError("nodata") << "Error! can't get the product retrieving EE crystal data " << digiCollection_.c_str();
0339 } else {
0340 EEDigi = pEEDigi.product();
0341 }
0342 } else {
0343 edm::LogError("cfg_error") << " Wrong ecalPart in cfg file ";
0344 return;
0345 }
0346
0347
0348 edm::Handle<EcalPnDiodeDigiCollection> pPNDigi;
0349 const EcalPnDiodeDigiCollection* PNDigi = nullptr;
0350 e.getByToken(pnDiodeDigiToken_, pPNDigi);
0351 if (!pPNDigi.isValid()) {
0352 edm::LogError("nodata") << "Error! can't get the product " << digiPNCollection_.c_str();
0353 } else {
0354 PNDigi = pPNDigi.product();
0355 }
0356
0357
0358 const auto& TheMapping = c.getData(mappingToken_);
0359
0360
0361
0362
0363
0364 for (EcalRawDataCollection::const_iterator headerItr = DCCHeader->begin(); headerItr != DCCHeader->end();
0365 ++headerItr) {
0366
0367
0368 int fed = headerItr->fedId();
0369 if (fed != _fedid && _fedid != -999)
0370 continue;
0371
0372 runType = headerItr->getRunType();
0373 runNum = headerItr->getRunNumber();
0374 event = headerItr->getLV1();
0375
0376 dccID = headerItr->getDccInTCCCommand();
0377 fedID = headerItr->fedId();
0378 lightside = headerItr->getRtHalf();
0379
0380
0381
0382 if (600 + dccID != fedID)
0383 continue;
0384
0385
0386
0387 if (runType != EcalDCCHeaderBlock::LASER_STD && runType != EcalDCCHeaderBlock::LASER_GAP &&
0388 runType != EcalDCCHeaderBlock::LASER_POWER_SCAN && runType != EcalDCCHeaderBlock::LASER_DELAY_SCAN)
0389 return;
0390
0391
0392
0393 EcalDCCHeaderBlock::EcalDCCEventSettings settings = headerItr->getEventSettings();
0394 color = settings.wavelength;
0395 if (color < 0)
0396 return;
0397
0398 std::vector<int>::iterator iter = find(colors.begin(), colors.end(), color);
0399 if (iter == colors.end()) {
0400 colors.push_back(color);
0401 edm::LogVerbatim("EcalLaserAnalyzer") << " new color found " << color << " " << colors.size();
0402 }
0403 }
0404
0405
0406
0407 if (fedID != _fedid && _fedid != -999)
0408 return;
0409
0410
0411
0412 laserEvents++;
0413
0414
0415
0416
0417
0418 TPNFit* pnfit = new TPNFit();
0419 pnfit->init(_nsamplesPN, _firstsamplePN, _lastsamplePN);
0420
0421 double chi2pn = 0;
0422 unsigned int samplemax = 0;
0423 int pnGain = 0;
0424
0425 std::map<int, std::vector<double> > allPNAmpl;
0426 std::map<int, std::vector<double> > allPNGain;
0427
0428
0429
0430 for (EcalPnDiodeDigiCollection::const_iterator pnItr = PNDigi->begin(); pnItr != PNDigi->end(); ++pnItr) {
0431 EcalPnDiodeDetId pnDetId = EcalPnDiodeDetId((*pnItr).id());
0432
0433 if (_debug == 1)
0434 edm::LogVerbatim("EcalLaserAnalyzer")
0435 << "-- debug -- Inside PNDigi - pnID=" << pnDetId.iPnId() << ", dccID=" << pnDetId.iDCCId();
0436
0437
0438
0439 bool isMemRelevant = Mem->isMemRelevant(pnDetId.iDCCId());
0440 if (!isMemRelevant)
0441 continue;
0442
0443
0444
0445 for (int samId = 0; samId < (*pnItr).size(); samId++) {
0446 pn[samId] = (*pnItr).sample(samId).adc();
0447 pnG[samId] = (*pnItr).sample(samId).gainId();
0448 if (samId == 0)
0449 pnGain = pnG[samId];
0450 if (samId > 0)
0451 pnGain = int(TMath::Max(pnG[samId], pnGain));
0452 }
0453
0454 if (pnGain != 1)
0455 edm::LogVerbatim("EcalLaserAnalyzer") << "PN gain different from 1";
0456
0457
0458
0459 PNPulse->setPulse(pn);
0460 pnNoPed = PNPulse->getAdcWithoutPedestal();
0461 samplemax = PNPulse->getMaxSample();
0462 chi2pn = pnfit->doFit(samplemax, &pnNoPed[0]);
0463 if (chi2pn == 101 || chi2pn == 102 || chi2pn == 103)
0464 pnAmpl = 0.;
0465 else
0466 pnAmpl = pnfit->getAmpl();
0467
0468
0469
0470 double corr = 1.0;
0471 if (_docorpn)
0472 corr = pnCorrector->getPNCorrectionFactor(pnAmpl, pnGain);
0473 pnAmpl *= corr;
0474
0475
0476
0477 allPNAmpl[pnDetId.iDCCId()].push_back(pnAmpl);
0478
0479 if (_debug == 1)
0480 edm::LogVerbatim("EcalLaserAnalyzer") << "-- debug -- Inside PNDigi - PNampl=" << pnAmpl << ", PNgain=" << pnGain;
0481 }
0482
0483
0484
0485
0486
0487 int adcGain = 0;
0488
0489 if (EBDigi) {
0490
0491
0492
0493 for (EBDigiCollection::const_iterator digiItr = EBDigi->begin(); digiItr != EBDigi->end(); ++digiItr) {
0494
0495
0496
0497 EBDetId id_crystal(digiItr->id());
0498 EBDataFrame df(*digiItr);
0499 EcalElectronicsId elecid_crystal = TheMapping.getElectronicsId(id_crystal);
0500
0501 int etaG = id_crystal.ieta();
0502 int phiG = id_crystal.iphi();
0503
0504 std::pair<int, int> LocalCoord = MEEBGeom::localCoord(etaG, phiG);
0505
0506 int etaL = LocalCoord.first;
0507 int phiL = LocalCoord.second;
0508
0509 int strip = elecid_crystal.stripId();
0510 int xtal = elecid_crystal.xtalId();
0511
0512 int module = MEEBGeom::lmmod(etaG, phiG);
0513 int tower = elecid_crystal.towerId();
0514
0515 int apdRefTT = MEEBGeom::apdRefTower(module);
0516
0517 std::pair<int, int> pnpair = MEEBGeom::pn(module);
0518 unsigned int MyPn0 = pnpair.first;
0519 unsigned int MyPn1 = pnpair.second;
0520
0521 int lmr = MEEBGeom::lmr(etaG, phiG);
0522 unsigned int channel = MEEBGeom::electronic_channel(etaL, phiL);
0523 assert(channel < nCrys);
0524
0525 setGeomEB(etaG, phiG, module, tower, strip, xtal, apdRefTT, channel, lmr);
0526
0527 if (_debug == 1)
0528 edm::LogVerbatim("EcalLaserAnalyzer")
0529 << "-- debug -- Inside EBDigi - towerID:" << towerID << " channelID:" << channelID << " module:" << module
0530 << " modules:" << modules.size();
0531
0532
0533
0534
0535
0536
0537 for (unsigned int i = 0; i < (*digiItr).size(); ++i) {
0538 EcalMGPASample samp_crystal(df.sample(i));
0539 adc[i] = samp_crystal.adc();
0540 adcG[i] = samp_crystal.gainId();
0541 adc[i] *= adcG[i];
0542 if (i == 0)
0543 adcGain = adcG[i];
0544 if (i > 0)
0545 adcGain = TMath::Max(adcG[i], adcGain);
0546 }
0547
0548 APDPulse->setPulse(adc);
0549
0550
0551
0552
0553 if (adcGain != 1)
0554 nEvtBadGain[channel]++;
0555 if (!APDPulse->isTimingQualOK())
0556 nEvtBadTiming[channel]++;
0557 nEvtTot[channel]++;
0558
0559
0560
0561
0562 int mem0 = Mem->Mem(lmr, 0);
0563 int mem1 = Mem->Mem(lmr, 1);
0564
0565 if (allPNAmpl[mem0].size() > MyPn0)
0566 pn0 = allPNAmpl[mem0][MyPn0];
0567 else
0568 pn0 = 0;
0569 if (allPNAmpl[mem1].size() > MyPn1)
0570 pn1 = allPNAmpl[mem1][MyPn1];
0571 else
0572 pn1 = 0;
0573
0574
0575
0576
0577 if (APDPulse->isPulseOK() && lightside == side) {
0578 ADCtrees[channel]->Fill();
0579
0580 Delta01->addEntry(APDPulse->getDelta(0, 1));
0581 Delta12->addEntry(APDPulse->getDelta(1, 2));
0582
0583 if (nevtAB[channel] < _nevtmax && _fitab) {
0584 if (doesABTreeExist)
0585 shapana->putAllVals(channel, adc, eta, phi);
0586 else
0587 shapana->putAllVals(channel, adc, eta, phi, dccID, side, towerID, channelID);
0588 nevtAB[channel]++;
0589 }
0590 }
0591 }
0592
0593 } else if (EEDigi) {
0594
0595
0596
0597 for (EEDigiCollection::const_iterator digiItr = EEDigi->begin(); digiItr != EEDigi->end(); ++digiItr) {
0598
0599
0600
0601 EEDetId id_crystal(digiItr->id());
0602 EEDataFrame df(*digiItr);
0603 EcalElectronicsId elecid_crystal = TheMapping.getElectronicsId(id_crystal);
0604
0605 int etaG = id_crystal.iy();
0606 int phiG = id_crystal.ix();
0607
0608 int iX = (phiG - 1) / 5 + 1;
0609 int iY = (etaG - 1) / 5 + 1;
0610
0611 int tower = elecid_crystal.towerId();
0612 int ch = elecid_crystal.channelId() - 1;
0613
0614 int module = MEEEGeom::lmmod(iX, iY);
0615 if (module >= 18 && side == 1)
0616 module += 2;
0617 int lmr = MEEEGeom::lmr(iX, iY, iZ);
0618 int dee = MEEEGeom::dee(lmr);
0619 int apdRefTT = MEEEGeom::apdRefTower(lmr, module);
0620
0621 std::pair<int, int> pnpair = MEEEGeom::pn(dee, module);
0622 unsigned int MyPn0 = pnpair.first;
0623 unsigned int MyPn1 = pnpair.second;
0624
0625 int hashedIndex = 100000 * eta + phi;
0626 if (channelMapEE.count(hashedIndex) == 0) {
0627 channelMapEE[hashedIndex] = channelIteratorEE;
0628 channelIteratorEE++;
0629 }
0630 unsigned int channel = channelMapEE[hashedIndex];
0631 assert(channel < nCrys);
0632
0633 setGeomEE(etaG, phiG, iX, iY, iZ, module, tower, ch, apdRefTT, channel, lmr);
0634
0635 if (_debug == 1)
0636 edm::LogVerbatim("EcalLaserAnalyzer")
0637 << "-- debug -- Inside EEDigi - towerID:" << towerID << " channelID:" << channelID << " module:" << module
0638 << " modules:" << modules.size();
0639
0640
0641
0642
0643 if ((*digiItr).size() > 10)
0644 edm::LogVerbatim("EcalLaserAnalyzer") << "SAMPLES SIZE > 10!" << (*digiItr).size();
0645
0646
0647
0648 for (unsigned int i = 0; i < (*digiItr).size(); ++i) {
0649 EcalMGPASample samp_crystal(df.sample(i));
0650 adc[i] = samp_crystal.adc();
0651 adcG[i] = samp_crystal.gainId();
0652 adc[i] *= adcG[i];
0653
0654 if (i == 0)
0655 adcGain = adcG[i];
0656 if (i > 0)
0657 adcGain = TMath::Max(adcG[i], adcGain);
0658 }
0659
0660 APDPulse->setPulse(adc);
0661
0662
0663
0664
0665 if (adcGain != 1)
0666 nEvtBadGain[channel]++;
0667 if (!APDPulse->isTimingQualOK())
0668 nEvtBadTiming[channel]++;
0669 nEvtTot[channel]++;
0670
0671
0672
0673
0674 int mem0 = Mem->Mem(lmr, 0);
0675 int mem1 = Mem->Mem(lmr, 1);
0676
0677 if (allPNAmpl[mem0].size() > MyPn0)
0678 pn0 = allPNAmpl[mem0][MyPn0];
0679 else
0680 pn0 = 0;
0681 if (allPNAmpl[mem1].size() > MyPn1)
0682 pn1 = allPNAmpl[mem1][MyPn1];
0683 else
0684 pn1 = 0;
0685
0686
0687
0688
0689 if (APDPulse->isPulseOK() && lightside == side) {
0690 ADCtrees[channel]->Fill();
0691
0692 Delta01->addEntry(APDPulse->getDelta(0, 1));
0693 Delta12->addEntry(APDPulse->getDelta(1, 2));
0694
0695 if (nevtAB[channel] < _nevtmax && _fitab) {
0696 if (doesABTreeExist)
0697 shapana->putAllVals(channel, adc, eta, phi);
0698 else
0699 shapana->putAllVals(channel, adc, eta, phi, dccID, side, towerID, channelID);
0700 nevtAB[channel]++;
0701 }
0702 }
0703 }
0704 }
0705 }
0706
0707
0708
0709 void EcalLaserAnalyzer::endJob() {
0710
0711
0712
0713
0714
0715 if (_ecalPart == "EE") {
0716 nCrys = channelMapEE.size();
0717 shapana->set_nch(nCrys);
0718 }
0719
0720
0721
0722
0723 double delta01 = Delta01->getMean();
0724 double delta12 = Delta12->getMean();
0725 if (delta12 > _presamplecut) {
0726 _presample = 2;
0727 if (delta01 > _presamplecut)
0728 _presample = 1;
0729 }
0730
0731 APDPulse->setPresamples(_presample);
0732 shapana->set_presample(_presample);
0733
0734
0735
0736
0737 if (_fitab) {
0738 edm::LogVerbatim("EcalLaserAnalyzer") << "\n\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+";
0739 edm::LogVerbatim("EcalLaserAnalyzer") << "\t+=+ Analyzing data: getting (alpha, beta) +=+";
0740 TFile* fAB = nullptr;
0741 TTree* ABInit = nullptr;
0742 if (doesABTreeExist) {
0743 fAB = new TFile(alphainitfile.c_str());
0744 }
0745 if (doesABTreeExist && fAB) {
0746 ABInit = (TTree*)fAB->Get("ABCol0");
0747 }
0748 shapana->computeShape(alphafile, ABInit);
0749 edm::LogVerbatim("EcalLaserAnalyzer") << "\t+=+ .................................... done +=+";
0750 edm::LogVerbatim("EcalLaserAnalyzer") << "\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+";
0751 }
0752
0753
0754
0755
0756 if (laserEvents == 0) {
0757 ADCFile->Close();
0758 std::stringstream del;
0759 del << "rm " << ADCfile;
0760 system(del.str().c_str());
0761 edm::LogVerbatim("EcalLaserAnalyzer") << " No Laser Events ";
0762 return;
0763 }
0764
0765
0766
0767
0768 double BadGainEvtPercentage = 0.0;
0769 double BadTimingEvtPercentage = 0.0;
0770
0771 int nChanBadGain = 0;
0772 int nChanBadTiming = 0;
0773
0774 for (unsigned int i = 0; i < nCrys; i++) {
0775 if (nEvtTot[i] != 0) {
0776 BadGainEvtPercentage = double(nEvtBadGain[i]) / double(nEvtTot[i]);
0777 BadTimingEvtPercentage = double(nEvtBadTiming[i]) / double(nEvtTot[i]);
0778 }
0779 if (BadGainEvtPercentage > _qualpercent) {
0780 wasGainOK[i] = false;
0781 nChanBadGain++;
0782 }
0783 if (BadTimingEvtPercentage > _qualpercent) {
0784 wasTimingOK[i] = false;
0785 nChanBadTiming++;
0786 }
0787 }
0788
0789 double BadGainChanPercentage = double(nChanBadGain) / double(nCrys);
0790 double BadTimingChanPercentage = double(nChanBadTiming) / double(nCrys);
0791
0792 if (BadGainChanPercentage > _qualpercent)
0793 isGainOK = false;
0794 if (BadTimingChanPercentage > _qualpercent)
0795 isTimingOK = false;
0796
0797
0798
0799
0800 edm::LogVerbatim("EcalLaserAnalyzer") << "\n\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+";
0801 edm::LogVerbatim("EcalLaserAnalyzer") << "\t+=+ Analyzing laser data: getting APD, PN, APD/PN, PN/PN +=+";
0802
0803 if (!isGainOK)
0804 edm::LogVerbatim("EcalLaserAnalyzer") << "\t+=+ ............................ WARNING! APD GAIN WAS NOT 1 +=+";
0805 if (!isTimingOK)
0806 edm::LogVerbatim("EcalLaserAnalyzer") << "\t+=+ ............................ WARNING! TIMING WAS BAD +=+";
0807
0808 APDFile = new TFile(APDfile.c_str(), "RECREATE");
0809
0810 int ieta, iphi, channelID, towerID, flag;
0811
0812 for (unsigned int i = 0; i < nCrys; i++) {
0813 std::stringstream name;
0814 name << "APDTree" << i + 1;
0815
0816 APDtrees[i] = new TTree(name.str().c_str(), name.str().c_str());
0817
0818
0819
0820 APDtrees[i]->Branch("event", &event, "event/I");
0821 APDtrees[i]->Branch("color", &color, "color/I");
0822 APDtrees[i]->Branch("iphi", &iphi, "iphi/I");
0823 APDtrees[i]->Branch("ieta", &ieta, "ieta/I");
0824 APDtrees[i]->Branch("side", &side, "side/I");
0825 APDtrees[i]->Branch("dccID", &dccID, "dccID/I");
0826 APDtrees[i]->Branch("towerID", &towerID, "towerID/I");
0827 APDtrees[i]->Branch("channelID", &channelID, "channelID/I");
0828 APDtrees[i]->Branch("apdAmpl", &apdAmpl, "apdAmpl/D");
0829 APDtrees[i]->Branch("apdTime", &apdTime, "apdTime/D");
0830 if (_saveallevents)
0831 APDtrees[i]->Branch("adc", &adc, "adc[10]/D");
0832 APDtrees[i]->Branch("flag", &flag, "flag/I");
0833 APDtrees[i]->Branch("flagAB", &flagAB, "flagAB/I");
0834 APDtrees[i]->Branch("pn0", &pn0, "pn0/D");
0835 APDtrees[i]->Branch("pn1", &pn1, "pn1/D");
0836
0837 APDtrees[i]->SetBranchAddress("event", &event);
0838 APDtrees[i]->SetBranchAddress("color", &color);
0839 APDtrees[i]->SetBranchAddress("iphi", &iphi);
0840 APDtrees[i]->SetBranchAddress("ieta", &ieta);
0841 APDtrees[i]->SetBranchAddress("side", &side);
0842 APDtrees[i]->SetBranchAddress("dccID", &dccID);
0843 APDtrees[i]->SetBranchAddress("towerID", &towerID);
0844 APDtrees[i]->SetBranchAddress("channelID", &channelID);
0845 APDtrees[i]->SetBranchAddress("apdAmpl", &apdAmpl);
0846 APDtrees[i]->SetBranchAddress("apdTime", &apdTime);
0847 if (_saveallevents)
0848 APDtrees[i]->SetBranchAddress("adc", adc);
0849 APDtrees[i]->SetBranchAddress("flag", &flag);
0850 APDtrees[i]->SetBranchAddress("flagAB", &flagAB);
0851 APDtrees[i]->SetBranchAddress("pn0", &pn0);
0852 APDtrees[i]->SetBranchAddress("pn1", &pn1);
0853 }
0854
0855 for (unsigned int iref = 0; iref < nRefChan; iref++) {
0856 for (unsigned int imod = 0; imod < nMod; imod++) {
0857 int jmod = modules[imod];
0858
0859 std::stringstream nameref;
0860 nameref << "refAPDTree" << imod << "_" << iref;
0861
0862 RefAPDtrees[iref][jmod] = new TTree(nameref.str().c_str(), nameref.str().c_str());
0863
0864 RefAPDtrees[iref][jmod]->Branch("eventref", &eventref, "eventref/I");
0865 RefAPDtrees[iref][jmod]->Branch("colorref", &colorref, "colorref/I");
0866 if (iref == 0)
0867 RefAPDtrees[iref][jmod]->Branch("apdAmplA", &apdAmplA, "apdAmplA/D");
0868 if (iref == 1)
0869 RefAPDtrees[iref][jmod]->Branch("apdAmplB", &apdAmplB, "apdAmplB/D");
0870
0871 RefAPDtrees[iref][jmod]->SetBranchAddress("eventref", &eventref);
0872 RefAPDtrees[iref][jmod]->SetBranchAddress("colorref", &colorref);
0873 if (iref == 0)
0874 RefAPDtrees[iref][jmod]->SetBranchAddress("apdAmplA", &apdAmplA);
0875 if (iref == 1)
0876 RefAPDtrees[iref][jmod]->SetBranchAddress("apdAmplB", &apdAmplB);
0877 }
0878 }
0879
0880 assert(colors.size() <= nColor);
0881 unsigned int nCol = colors.size();
0882
0883
0884
0885
0886 for (unsigned int iM = 0; iM < nMod; iM++) {
0887 unsigned int iMod = modules[iM] - 1;
0888
0889 for (unsigned int ich = 0; ich < nPNPerMod; ich++) {
0890 for (unsigned int icol = 0; icol < nCol; icol++) {
0891 PNFirstAnal[iMod][ich][icol] = new TPN(ich);
0892 PNAnal[iMod][ich][icol] = new TPN(ich);
0893 }
0894 }
0895 }
0896
0897
0898
0899
0900 PulseFitWithFunction* pslsfit = new PulseFitWithFunction();
0901 double chi2;
0902
0903 for (unsigned int iCry = 0; iCry < nCrys; iCry++) {
0904 for (unsigned int icol = 0; icol < nCol; icol++) {
0905
0906
0907
0908 APDFirstAnal[iCry][icol] = new TAPD();
0909 IsThereDataADC[iCry][icol] = 1;
0910 std::stringstream cut;
0911 cut << "color==" << colors[icol];
0912 if (ADCtrees[iCry]->GetEntries(cut.str().c_str()) < 10)
0913 IsThereDataADC[iCry][icol] = 0;
0914 }
0915
0916 unsigned int iMod = iModule[iCry] - 1;
0917 double alpha, beta;
0918
0919
0920
0921
0922 Long64_t nbytes = 0, nb = 0;
0923 for (Long64_t jentry = 0; jentry < ADCtrees[iCry]->GetEntriesFast(); jentry++) {
0924 nb = ADCtrees[iCry]->GetEntry(jentry);
0925 nbytes += nb;
0926
0927
0928
0929 unsigned int iCol = 0;
0930 for (unsigned int i = 0; i < nCol; i++) {
0931 if (color == colors[i]) {
0932 iCol = i;
0933 i = colors.size();
0934 }
0935 }
0936
0937
0938
0939 std::vector<double> abvals = shapana->getVals(iCry);
0940 alpha = abvals[0];
0941 beta = abvals[1];
0942 flagAB = int(abvals[4]);
0943 iphi = iPhi[iCry];
0944 ieta = iEta[iCry];
0945 towerID = iTowerID[iCry];
0946 channelID = iChannelID[iCry];
0947
0948
0949
0950 APDPulse->setPulse(adc);
0951 adcNoPed = APDPulse->getAdcWithoutPedestal();
0952
0953 apdAmpl = 0;
0954 apdAmplA = 0;
0955 apdAmplB = 0;
0956 apdTime = 0;
0957
0958 if (APDPulse->isPulseOK()) {
0959 pslsfit->init(_nsamples, _firstsample, _lastsample, _niter, alpha, beta);
0960 chi2 = pslsfit->doFit(&adcNoPed[0]);
0961
0962 if (chi2 < 0. || chi2 == 102 || chi2 == 101) {
0963 apdAmpl = 0;
0964 apdTime = 0;
0965 flag = 0;
0966 } else {
0967 apdAmpl = pslsfit->getAmpl();
0968 apdTime = pslsfit->getTime();
0969 flag = 1;
0970 }
0971 } else {
0972 apdAmpl = 0;
0973 apdTime = 0;
0974 flag = 0;
0975 }
0976
0977 if (_debug == 1)
0978 edm::LogVerbatim("EcalLaserAnalyzer") << "-- debug test -- apdAmpl=" << apdAmpl << ", apdTime=" << apdTime;
0979 double pnmean;
0980 if (pn0 < 10 && pn1 > 10) {
0981 pnmean = pn1;
0982 } else if (pn1 < 10 && pn0 > 10) {
0983 pnmean = pn0;
0984 } else
0985 pnmean = 0.5 * (pn0 + pn1);
0986
0987 if (_debug == 1)
0988 edm::LogVerbatim("EcalLaserAnalyzer") << "-- debug test -- pn0=" << pn0 << ", pn1=" << pn1;
0989
0990
0991
0992
0993 if (firstChanMod[iMod] == iCry && IsThereDataADC[iCry][iCol] == 1) {
0994 for (unsigned int ichan = 0; ichan < nPNPerMod; ichan++) {
0995 PNFirstAnal[iMod][ichan][iCol]->addEntry(pnmean, pn0, pn1);
0996 }
0997 }
0998
0999
1000
1001
1002 if (APDPulse->isPulseOK()) {
1003 APDFirstAnal[iCry][iCol]->addEntry(apdAmpl, pnmean, pn0, pn1, apdTime);
1004 APDtrees[iCry]->Fill();
1005
1006
1007
1008
1009 if (apdRefMap[0][iMod + 1] == iCry || apdRefMap[1][iMod + 1] == iCry) {
1010 apdAmplA = 0.0;
1011 apdAmplB = 0.0;
1012 eventref = event;
1013 colorref = color;
1014
1015 for (unsigned int ir = 0; ir < nRefChan; ir++) {
1016 if (apdRefMap[ir][iMod + 1] == iCry) {
1017 if (ir == 0)
1018 apdAmplA = apdAmpl;
1019 else if (ir == 1)
1020 apdAmplB = apdAmpl;
1021 RefAPDtrees[ir][iMod + 1]->Fill();
1022 }
1023 }
1024 }
1025 }
1026 }
1027 }
1028
1029 delete pslsfit;
1030
1031 ADCFile->Close();
1032
1033
1034
1035 std::stringstream del;
1036 del << "rm " << ADCfile;
1037 system(del.str().c_str());
1038
1039
1040
1041
1042 resFile = new TFile(resfile.c_str(), "RECREATE");
1043
1044 for (unsigned int iColor = 0; iColor < nCol; iColor++) {
1045 std::stringstream nametree;
1046 nametree << "APDCol" << colors[iColor];
1047 std::stringstream nametree2;
1048 nametree2 << "PNCol" << colors[iColor];
1049
1050 restrees[iColor] = new TTree(nametree.str().c_str(), nametree.str().c_str());
1051 respntrees[iColor] = new TTree(nametree2.str().c_str(), nametree2.str().c_str());
1052
1053 restrees[iColor]->Branch("iphi", &iphi, "iphi/I");
1054 restrees[iColor]->Branch("ieta", &ieta, "ieta/I");
1055 restrees[iColor]->Branch("side", &side, "side/I");
1056 restrees[iColor]->Branch("dccID", &dccID, "dccID/I");
1057 restrees[iColor]->Branch("moduleID", &moduleID, "moduleID/I");
1058 restrees[iColor]->Branch("towerID", &towerID, "towerID/I");
1059 restrees[iColor]->Branch("channelID", &channelID, "channelID/I");
1060 restrees[iColor]->Branch("APD", &APD, "APD[6]/D");
1061 restrees[iColor]->Branch("Time", &Time, "Time[6]/D");
1062 restrees[iColor]->Branch("APDoPN", &APDoPN, "APDoPN[6]/D");
1063 restrees[iColor]->Branch("APDoPNA", &APDoPNA, "APDoPNA[6]/D");
1064 restrees[iColor]->Branch("APDoPNB", &APDoPNB, "APDoPNB[6]/D");
1065 restrees[iColor]->Branch("APDoAPDA", &APDoAPDA, "APDoAPDA[6]/D");
1066 restrees[iColor]->Branch("APDoAPDB", &APDoAPDB, "APDoAPDB[6]/D");
1067 restrees[iColor]->Branch("flag", &flag, "flag/I");
1068
1069 respntrees[iColor]->Branch("side", &side, "side/I");
1070 respntrees[iColor]->Branch("moduleID", &moduleID, "moduleID/I");
1071 respntrees[iColor]->Branch("pnID", &pnID, "pnID/I");
1072 respntrees[iColor]->Branch("PN", &PN, "PN[6]/D");
1073 respntrees[iColor]->Branch("PNoPN", &PNoPN, "PNoPN[6]/D");
1074 respntrees[iColor]->Branch("PNoPNA", &PNoPNA, "PNoPNA[6]/D");
1075 respntrees[iColor]->Branch("PNoPNB", &PNoPNB, "PNoPNB[6]/D");
1076
1077 restrees[iColor]->SetBranchAddress("iphi", &iphi);
1078 restrees[iColor]->SetBranchAddress("ieta", &ieta);
1079 restrees[iColor]->SetBranchAddress("side", &side);
1080 restrees[iColor]->SetBranchAddress("dccID", &dccID);
1081 restrees[iColor]->SetBranchAddress("moduleID", &moduleID);
1082 restrees[iColor]->SetBranchAddress("towerID", &towerID);
1083 restrees[iColor]->SetBranchAddress("channelID", &channelID);
1084 restrees[iColor]->SetBranchAddress("APD", APD);
1085 restrees[iColor]->SetBranchAddress("Time", Time);
1086 restrees[iColor]->SetBranchAddress("APDoPN", APDoPN);
1087 restrees[iColor]->SetBranchAddress("APDoPNA", APDoPNA);
1088 restrees[iColor]->SetBranchAddress("APDoPNB", APDoPNB);
1089 restrees[iColor]->SetBranchAddress("APDoAPDA", APDoAPDA);
1090 restrees[iColor]->SetBranchAddress("APDoAPDB", APDoAPDB);
1091 restrees[iColor]->SetBranchAddress("flag", &flag);
1092
1093 respntrees[iColor]->SetBranchAddress("side", &side);
1094 respntrees[iColor]->SetBranchAddress("moduleID", &moduleID);
1095 respntrees[iColor]->SetBranchAddress("pnID", &pnID);
1096 respntrees[iColor]->SetBranchAddress("PN", PN);
1097 respntrees[iColor]->SetBranchAddress("PNoPN", PNoPN);
1098 respntrees[iColor]->SetBranchAddress("PNoPNA", PNoPNA);
1099 respntrees[iColor]->SetBranchAddress("PNoPNB", PNoPNB);
1100 }
1101
1102
1103
1104
1105 for (unsigned int iM = 0; iM < nMod; iM++) {
1106 unsigned int iMod = modules[iM] - 1;
1107
1108 for (unsigned int ich = 0; ich < nPNPerMod; ich++) {
1109 for (unsigned int icol = 0; icol < nCol; icol++) {
1110 PNAnal[iMod][ich][icol]->setPNCut(PNFirstAnal[iMod][ich][icol]->getPN().at(0),
1111 PNFirstAnal[iMod][ich][icol]->getPN().at(1));
1112 }
1113 }
1114 }
1115
1116
1117
1118 for (unsigned int imod = 0; imod < nMod; imod++) {
1119 int jmod = modules[imod];
1120 if (RefAPDtrees[0][jmod]->GetEntries() != 0 && RefAPDtrees[1][jmod]->GetEntries() != 0) {
1121 RefAPDtrees[0][jmod]->BuildIndex("eventref");
1122 RefAPDtrees[1][jmod]->BuildIndex("eventref");
1123 }
1124 }
1125
1126
1127
1128
1129 for (unsigned int iCry = 0; iCry < nCrys; iCry++) {
1130 unsigned int iMod = iModule[iCry] - 1;
1131
1132
1133
1134
1135 for (unsigned int iCol = 0; iCol < nCol; iCol++) {
1136 std::vector<double> lowcut;
1137 std::vector<double> highcut;
1138 double cutMin;
1139 double cutMax;
1140
1141 cutMin = APDFirstAnal[iCry][iCol]->getAPD().at(0) - 2.0 * APDFirstAnal[iCry][iCol]->getAPD().at(1);
1142 if (cutMin < 0)
1143 cutMin = 0;
1144 cutMax = APDFirstAnal[iCry][iCol]->getAPD().at(0) + 2.0 * APDFirstAnal[iCry][iCol]->getAPD().at(1);
1145
1146 lowcut.push_back(cutMin);
1147 highcut.push_back(cutMax);
1148
1149 cutMin = APDFirstAnal[iCry][iCol]->getTime().at(0) - 2.0 * APDFirstAnal[iCry][iCol]->getTime().at(1);
1150 cutMax = APDFirstAnal[iCry][iCol]->getTime().at(0) + 2.0 * APDFirstAnal[iCry][iCol]->getTime().at(1);
1151 lowcut.push_back(cutMin);
1152 highcut.push_back(cutMax);
1153
1154 APDAnal[iCry][iCol] = new TAPD();
1155 APDAnal[iCry][iCol]->setAPDCut(APDFirstAnal[iCry][iCol]->getAPD().at(0),
1156 APDFirstAnal[iCry][iCol]->getAPD().at(1));
1157 APDAnal[iCry][iCol]->setAPDoPNCut(APDFirstAnal[iCry][iCol]->getAPDoPN().at(0),
1158 APDFirstAnal[iCry][iCol]->getAPDoPN().at(1));
1159 APDAnal[iCry][iCol]->setAPDoPN0Cut(APDFirstAnal[iCry][iCol]->getAPDoPN0().at(0),
1160 APDFirstAnal[iCry][iCol]->getAPDoPN0().at(1));
1161 APDAnal[iCry][iCol]->setAPDoPN1Cut(APDFirstAnal[iCry][iCol]->getAPDoPN1().at(0),
1162 APDFirstAnal[iCry][iCol]->getAPDoPN1().at(1));
1163 APDAnal[iCry][iCol]->setTimeCut(APDFirstAnal[iCry][iCol]->getTime().at(0),
1164 APDFirstAnal[iCry][iCol]->getTime().at(1));
1165 APDAnal[iCry][iCol]->set2DAPDoAPD0Cut(lowcut, highcut);
1166 APDAnal[iCry][iCol]->set2DAPDoAPD1Cut(lowcut, highcut);
1167 }
1168
1169
1170
1171
1172 Long64_t nbytes = 0, nb = 0;
1173 for (Long64_t jentry = 0; jentry < APDtrees[iCry]->GetEntriesFast(); jentry++) {
1174 nb = APDtrees[iCry]->GetEntry(jentry);
1175 nbytes += nb;
1176
1177 double pnmean;
1178 if (pn0 < 10 && pn1 > 10) {
1179 pnmean = pn1;
1180 } else if (pn1 < 10 && pn0 > 10) {
1181 pnmean = pn0;
1182 } else
1183 pnmean = 0.5 * (pn0 + pn1);
1184
1185
1186
1187
1188 unsigned int iCol = 0;
1189 for (unsigned int i = 0; i < nCol; i++) {
1190 if (color == colors[i]) {
1191 iCol = i;
1192 i = colors.size();
1193 }
1194 }
1195
1196
1197
1198
1199 if (firstChanMod[iMod] == iCry && IsThereDataADC[iCry][iCol] == 1) {
1200 for (unsigned int ichan = 0; ichan < nPNPerMod; ichan++) {
1201 PNAnal[iMod][ichan][iCol]->addEntry(pnmean, pn0, pn1);
1202 }
1203 }
1204
1205
1206
1207
1208 if (_debug == 1)
1209 edm::LogVerbatim("EcalLaserAnalyzer") << "-- debug test -- Last Loop event:" << event << " apdAmpl:" << apdAmpl;
1210 apdAmplA = 0.0;
1211 apdAmplB = 0.0;
1212
1213 for (unsigned int iRef = 0; iRef < nRefChan; iRef++) {
1214 RefAPDtrees[iRef][iMod + 1]->GetEntryWithIndex(event);
1215 }
1216
1217 if (_debug == 1)
1218 edm::LogVerbatim("EcalLaserAnalyzer")
1219 << "-- debug test -- Last Loop apdAmplA:" << apdAmplA << " apdAmplB:" << apdAmplB << ", event:" << event
1220 << ", eventref:" << eventref;
1221
1222
1223
1224
1225 APDAnal[iCry][iCol]->addEntry(apdAmpl, pnmean, pn0, pn1, apdTime, apdAmplA, apdAmplB);
1226 }
1227
1228 moduleID = iMod + 1;
1229
1230 if (moduleID >= 20)
1231 moduleID -= 2;
1232
1233
1234
1235
1236 for (unsigned int iColor = 0; iColor < nCol; iColor++) {
1237 std::vector<double> apdvec = APDAnal[iCry][iColor]->getAPD();
1238 std::vector<double> apdpnvec = APDAnal[iCry][iColor]->getAPDoPN();
1239 std::vector<double> apdpn0vec = APDAnal[iCry][iColor]->getAPDoPN0();
1240 std::vector<double> apdpn1vec = APDAnal[iCry][iColor]->getAPDoPN1();
1241 std::vector<double> timevec = APDAnal[iCry][iColor]->getTime();
1242 std::vector<double> apdapd0vec = APDAnal[iCry][iColor]->getAPDoAPD0();
1243 std::vector<double> apdapd1vec = APDAnal[iCry][iColor]->getAPDoAPD1();
1244
1245 for (unsigned int i = 0; i < apdvec.size(); i++) {
1246 APD[i] = apdvec.at(i);
1247 APDoPN[i] = apdpnvec.at(i);
1248 APDoPNA[i] = apdpn0vec.at(i);
1249 APDoPNB[i] = apdpn1vec.at(i);
1250 APDoAPDA[i] = apdapd0vec.at(i);
1251 APDoAPDB[i] = apdapd1vec.at(i);
1252 Time[i] = timevec.at(i);
1253 }
1254
1255
1256
1257
1258 iphi = iPhi[iCry];
1259 ieta = iEta[iCry];
1260 dccID = idccID[iCry];
1261 side = iside[iCry];
1262 towerID = iTowerID[iCry];
1263 channelID = iChannelID[iCry];
1264
1265 if (!wasGainOK[iCry] || !wasTimingOK[iCry] || IsThereDataADC[iCry][iColor] == 0) {
1266 flag = 0;
1267 } else
1268 flag = 1;
1269
1270 restrees[iColor]->Fill();
1271 }
1272 }
1273
1274
1275
1276
1277 for (unsigned int iM = 0; iM < nMod; iM++) {
1278 unsigned int iMod = modules[iM] - 1;
1279
1280 side = iside[firstChanMod[iMod]];
1281
1282 for (unsigned int ch = 0; ch < nPNPerMod; ch++) {
1283 pnID = ch;
1284 moduleID = iMod + 1;
1285
1286 if (moduleID >= 20)
1287 moduleID -= 2;
1288
1289 for (unsigned int iColor = 0; iColor < nCol; iColor++) {
1290 std::vector<double> pnvec = PNAnal[iMod][ch][iColor]->getPN();
1291 std::vector<double> pnopnvec = PNAnal[iMod][ch][iColor]->getPNoPN();
1292 std::vector<double> pnopn0vec = PNAnal[iMod][ch][iColor]->getPNoPN0();
1293 std::vector<double> pnopn1vec = PNAnal[iMod][ch][iColor]->getPNoPN1();
1294
1295 for (unsigned int i = 0; i < pnvec.size(); i++) {
1296 PN[i] = pnvec.at(i);
1297 PNoPN[i] = pnopnvec.at(i);
1298 PNoPNA[i] = pnopn0vec.at(i);
1299 PNoPNB[i] = pnopn1vec.at(i);
1300 }
1301
1302
1303
1304
1305 respntrees[iColor]->Fill();
1306 }
1307 }
1308 }
1309
1310
1311
1312 if (!_saveallevents) {
1313 APDFile->Close();
1314 std::stringstream del2;
1315 del2 << "rm " << APDfile;
1316 system(del2.str().c_str());
1317
1318 } else {
1319 APDFile->cd();
1320 APDtrees[0]->Write();
1321
1322 APDFile->Close();
1323 resFile->cd();
1324 }
1325
1326
1327
1328
1329 for (unsigned int i = 0; i < nCol; i++) {
1330 restrees[i]->Write();
1331 respntrees[i]->Write();
1332 }
1333
1334 resFile->Close();
1335
1336 edm::LogVerbatim("EcalLaserAnalyzer") << "\t+=+ .................................................. done +=+";
1337 edm::LogVerbatim("EcalLaserAnalyzer") << "\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+";
1338 }
1339
1340 void EcalLaserAnalyzer::setGeomEB(
1341 int etaG, int phiG, int module, int tower, int strip, int xtal, int apdRefTT, int channel, int lmr) {
1342 side = MEEBGeom::side(etaG, phiG);
1343
1344 assert(module >= *min_element(modules.begin(), modules.end()) &&
1345 module <= *max_element(modules.begin(), modules.end()));
1346
1347 eta = etaG;
1348 phi = phiG;
1349 channelID = 5 * (strip - 1) + xtal - 1;
1350 towerID = tower;
1351
1352 std::vector<int> apdRefChan = ME::apdRefChannels(module, lmr);
1353 for (unsigned int iref = 0; iref < nRefChan; iref++) {
1354 if (channelID == apdRefChan[iref] && towerID == apdRefTT && apdRefMap[iref].count(module) == 0) {
1355 apdRefMap[iref][module] = channel;
1356 }
1357 }
1358
1359 if (isFirstChanModFilled[module - 1] == 0) {
1360 firstChanMod[module - 1] = channel;
1361 isFirstChanModFilled[module - 1] = 1;
1362 }
1363
1364 iEta[channel] = eta;
1365 iPhi[channel] = phi;
1366 iModule[channel] = module;
1367 iTowerID[channel] = towerID;
1368 iChannelID[channel] = channelID;
1369 idccID[channel] = dccID;
1370 iside[channel] = side;
1371 }
1372
1373 void EcalLaserAnalyzer::setGeomEE(
1374 int etaG, int phiG, int iX, int iY, int iZ, int module, int tower, int ch, int apdRefTT, int channel, int lmr) {
1375 side = MEEEGeom::side(iX, iY, iZ);
1376
1377 assert(module >= *min_element(modules.begin(), modules.end()) &&
1378 module <= *max_element(modules.begin(), modules.end()));
1379
1380 eta = etaG;
1381 phi = phiG;
1382 channelID = ch;
1383 towerID = tower;
1384
1385 std::vector<int> apdRefChan = ME::apdRefChannels(module, lmr);
1386 for (unsigned int iref = 0; iref < nRefChan; iref++) {
1387 if (channelID == apdRefChan[iref] && towerID == apdRefTT && apdRefMap[iref].count(module) == 0) {
1388 apdRefMap[iref][module] = channel;
1389 }
1390 }
1391
1392 if (isFirstChanModFilled[module - 1] == 0) {
1393 firstChanMod[module - 1] = channel;
1394 isFirstChanModFilled[module - 1] = 1;
1395 }
1396
1397 iEta[channel] = eta;
1398 iPhi[channel] = phi;
1399 iModule[channel] = module;
1400 iTowerID[channel] = towerID;
1401 iChannelID[channel] = channelID;
1402 idccID[channel] = dccID;
1403 iside[channel] = side;
1404 }
1405
1406 DEFINE_FWK_MODULE(EcalLaserAnalyzer);