File indexing completed on 2024-04-06 11:57:49
0001
0002
0003
0004
0005
0006
0007
0008 #include "TFile.h"
0009 #include "TTree.h"
0010
0011 #include "EcalPerEvtLaserAnalyzer.h"
0012
0013 #include <sstream>
0014 #include <iomanip>
0015 #include <ctime>
0016
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018
0019 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/MEEEGeom.h"
0020 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/MEEBGeom.h"
0021
0022 #include "FWCore/Framework/interface/Event.h"
0023 #include "FWCore/Framework/interface/MakerMacros.h"
0024 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0025 #include "FWCore/Framework/interface/EventSetup.h"
0026
0027 #include "DataFormats/EcalDetId/interface/EcalElectronicsId.h"
0028 #include "DataFormats/EcalDetId/interface/EcalDetIdCollections.h"
0029
0030 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/TPNFit.h"
0031 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/PulseFitWithFunction.h"
0032
0033 using namespace std;
0034
0035
0036 EcalPerEvtLaserAnalyzer::EcalPerEvtLaserAnalyzer(const edm::ParameterSet& iConfig)
0037
0038 : iEvent(0),
0039 eventHeaderCollection_(iConfig.getParameter<std::string>("eventHeaderCollection")),
0040 eventHeaderProducer_(iConfig.getParameter<std::string>("eventHeaderProducer")),
0041 digiCollection_(iConfig.getParameter<std::string>("digiCollection")),
0042 digiProducer_(iConfig.getParameter<std::string>("digiProducer")),
0043 digiPNCollection_(iConfig.getParameter<std::string>("digiPNCollection")),
0044 rawDataToken_(consumes<EcalRawDataCollection>(edm::InputTag(eventHeaderProducer_, eventHeaderCollection_))),
0045 pnDiodeDigiToken_(consumes<EcalPnDiodeDigiCollection>(edm::InputTag(digiProducer_, digiPNCollection_))),
0046 mappingToken_(esConsumes()),
0047
0048 _nsamples(iConfig.getUntrackedParameter<unsigned int>("nSamples", 10)),
0049 _presample(iConfig.getUntrackedParameter<unsigned int>("nPresamples", 3)),
0050 _firstsample(iConfig.getUntrackedParameter<unsigned int>("firstSample", 1)),
0051 _lastsample(iConfig.getUntrackedParameter<unsigned int>("lastSample", 2)),
0052 _nsamplesPN(iConfig.getUntrackedParameter<unsigned int>("nSamplesPN", 50)),
0053 _presamplePN(iConfig.getUntrackedParameter<unsigned int>("nPresamplesPN", 6)),
0054 _firstsamplePN(iConfig.getUntrackedParameter<unsigned int>("firstSamplePN", 7)),
0055 _lastsamplePN(iConfig.getUntrackedParameter<unsigned int>("lastSamplePN", 8)),
0056 _timingcutlow(iConfig.getUntrackedParameter<unsigned int>("timingCutLow", 3)),
0057 _timingcuthigh(iConfig.getUntrackedParameter<unsigned int>("timingCutHigh", 7)),
0058 _niter(iConfig.getUntrackedParameter<unsigned int>("nIter", 3)),
0059 _fedid(iConfig.getUntrackedParameter<unsigned int>("fedID", 0)),
0060 _tower(iConfig.getUntrackedParameter<unsigned int>("tower", 1)),
0061 _channel(iConfig.getUntrackedParameter<unsigned int>("channel", 1)),
0062 _ecalPart(iConfig.getUntrackedParameter<std::string>("ecalPart", "EB")),
0063 resdir_(iConfig.getUntrackedParameter<std::string>("resDir")),
0064 refalphabeta_(iConfig.getUntrackedParameter<std::string>("refAlphaBeta")),
0065 nCrys(NCRYSEB),
0066 IsFileCreated(0),
0067 runType(-1),
0068 runNum(0),
0069 dccID(-1),
0070 lightside(2),
0071 doesRefFileExist(0),
0072 ttMat(-1),
0073 peakMat(-1),
0074 peak(-1),
0075 evtMat(-1),
0076 colMat(-1)
0077
0078
0079 {
0080 if (_ecalPart == "EB") {
0081 ebDigiToken_ = consumes<EBDigiCollection>(edm::InputTag(digiProducer_, digiCollection_));
0082 } else if (_ecalPart == "EE") {
0083 eeDigiToken_ = consumes<EEDigiCollection>(edm::InputTag(digiProducer_, digiCollection_));
0084 }
0085
0086
0087
0088 if (_ecalPart == "EB") {
0089 nCrys = NCRYSEB;
0090 } else {
0091 nCrys = NCRYSEE;
0092 }
0093 }
0094
0095
0096 EcalPerEvtLaserAnalyzer::~EcalPerEvtLaserAnalyzer() {
0097
0098
0099
0100
0101 }
0102
0103
0104 void EcalPerEvtLaserAnalyzer::beginJob() {
0105
0106
0107
0108
0109 stringstream namefile1;
0110 namefile1 << resdir_ << "/ADCSamples.root";
0111
0112 ADCfile = namefile1.str();
0113
0114
0115
0116 ADCFile = new TFile(ADCfile.c_str(), "RECREATE");
0117
0118 stringstream name;
0119 name << "ADCTree";
0120
0121 ADCtrees = new TTree(name.str().c_str(), name.str().c_str());
0122
0123 ADCtrees->Branch("iphi", &phi, "phi/I");
0124 ADCtrees->Branch("ieta", &eta, "eta/I");
0125 ADCtrees->Branch("dccID", &dccID, "dccID/I");
0126 ADCtrees->Branch("event", &event, "event/I");
0127 ADCtrees->Branch("color", &color, "color/I");
0128 ADCtrees->Branch("adc", &adc, "adc[10]/D");
0129 ADCtrees->Branch("ttrigMatacq", &ttrig, "ttrig/D");
0130 ADCtrees->Branch("peakMatacq", &peak, "peak/D");
0131 ADCtrees->Branch("pn0", &pn0, "pn0/D");
0132 ADCtrees->Branch("pn1", &pn1, "pn1/D");
0133
0134 ADCtrees->SetBranchAddress("ieta", &eta);
0135 ADCtrees->SetBranchAddress("iphi", &phi);
0136 ADCtrees->SetBranchAddress("dccID", &dccID);
0137 ADCtrees->SetBranchAddress("event", &event);
0138 ADCtrees->SetBranchAddress("color", &color);
0139 ADCtrees->SetBranchAddress("adc", adc);
0140 ADCtrees->SetBranchAddress("ttrigMatacq", &ttrig);
0141 ADCtrees->SetBranchAddress("peakMatacq", &peak);
0142 ADCtrees->SetBranchAddress("pn0", &pn0);
0143 ADCtrees->SetBranchAddress("pn1", &pn1);
0144
0145 IsFileCreated = 0;
0146 }
0147
0148
0149 void EcalPerEvtLaserAnalyzer::analyze(const edm::Event& e, const edm::EventSetup& c) {
0150
0151
0152 ++iEvent;
0153
0154
0155 edm::Handle<EcalRawDataCollection> pDCCHeader;
0156 const EcalRawDataCollection* DCCHeader = nullptr;
0157 e.getByToken(rawDataToken_, pDCCHeader);
0158 if (!pDCCHeader.isValid()) {
0159 edm::LogError("nodata") << "Error! can't get the product retrieving DCC header" << eventHeaderCollection_.c_str();
0160 } else {
0161 DCCHeader = pDCCHeader.product();
0162 }
0163
0164
0165 edm::Handle<EBDigiCollection> pEBDigi;
0166 const EBDigiCollection* EBDigi = nullptr;
0167 edm::Handle<EEDigiCollection> pEEDigi;
0168 const EEDigiCollection* EEDigi = nullptr;
0169 if (_ecalPart == "EB") {
0170 e.getByToken(ebDigiToken_, pEBDigi);
0171 if (!pEBDigi.isValid()) {
0172 edm::LogError("nodata") << "Error! can't get the product retrieving EB crystal data " << digiCollection_.c_str();
0173 } else {
0174 EBDigi = pEBDigi.product();
0175 }
0176 } else {
0177 e.getByToken(eeDigiToken_, pEEDigi);
0178 if (!pEEDigi.isValid()) {
0179 edm::LogError("nodata") << "Error! can't get the product retrieving EE crystal data " << digiCollection_.c_str();
0180 } else {
0181 EEDigi = pEEDigi.product();
0182 }
0183 }
0184
0185
0186 edm::Handle<EcalPnDiodeDigiCollection> pPNDigi;
0187 const EcalPnDiodeDigiCollection* PNDigi = nullptr;
0188 e.getByToken(pnDiodeDigiToken_, pPNDigi);
0189 if (!pPNDigi.isValid()) {
0190 edm::LogError("nodata") << "Error! can't get the product " << digiPNCollection_.c_str();
0191 } else {
0192 PNDigi = pPNDigi.product();
0193 }
0194
0195
0196 const auto& TheMapping = c.getData(mappingToken_);
0197
0198
0199
0200
0201
0202 for (EcalRawDataCollection::const_iterator headerItr = DCCHeader->begin(); headerItr != DCCHeader->end();
0203 ++headerItr) {
0204
0205
0206 int fed = headerItr->fedId();
0207
0208 if (fed != _fedid && _fedid != -999)
0209 continue;
0210
0211 runType = headerItr->getRunType();
0212 runNum = headerItr->getRunNumber();
0213 event = headerItr->getLV1();
0214 dccID = headerItr->getDccInTCCCommand();
0215 fedID = headerItr->fedId();
0216
0217
0218 if (600 + dccID != fedID)
0219 continue;
0220
0221
0222 if (runType != EcalDCCHeaderBlock::LASER_STD && runType != EcalDCCHeaderBlock::LASER_GAP)
0223 return;
0224
0225
0226
0227 if (IsFileCreated == 0) {
0228 stringstream namefile2;
0229 namefile2 << resdir_ << "/APDAmpl_Run" << runNum << "_" << _fedid << "_" << _tower << "_" << _channel << ".root";
0230 resfile = namefile2.str();
0231
0232
0233
0234 stringstream namefile;
0235 namefile << resdir_ << "/Matacq-Run" << runNum << ".root";
0236
0237 doesRefFileExist = 0;
0238
0239 FILE* test;
0240 test = fopen(namefile.str().c_str(), "r");
0241 if (test)
0242 doesRefFileExist = 1;
0243
0244 if (doesRefFileExist == 1) {
0245 matacqFile = new TFile((namefile.str().c_str()));
0246 matacqTree = (TTree*)matacqFile->Get("MatacqShape");
0247
0248 matacqTree->SetBranchAddress("event", &evtMat);
0249 matacqTree->SetBranchAddress("color", &colMat);
0250 matacqTree->SetBranchAddress("peak", &peakMat);
0251 matacqTree->SetBranchAddress("ttrig", &ttMat);
0252 }
0253
0254 IsFileCreated = 1;
0255 }
0256
0257
0258
0259 EcalDCCHeaderBlock::EcalDCCEventSettings settings = headerItr->getEventSettings();
0260 int color = settings.wavelength;
0261
0262 vector<int>::iterator iter = find(colors.begin(), colors.end(), color);
0263 if (iter == colors.end()) {
0264 colors.push_back(color);
0265 edm::LogVerbatim("EcalPerEvtLaserAnalyzer") << " new color found " << color << " " << colors.size();
0266 }
0267 }
0268
0269
0270
0271 if (fedID != _fedid && _fedid != -999)
0272 return;
0273
0274
0275
0276
0277
0278 TPNFit* pnfit = new TPNFit();
0279 pnfit->init(_nsamplesPN, _firstsamplePN, _lastsamplePN);
0280
0281 double chi2pn = 0;
0282 double ypnrange[50];
0283 double dsum = 0.;
0284 double dsum1 = 0.;
0285 double bl = 0.;
0286 double bl1 = 0.;
0287 double val_max = 0.;
0288 unsigned int samplemax = 0;
0289 unsigned int k;
0290
0291 std::vector<double> allPNAmpl;
0292
0293 for (EcalPnDiodeDigiCollection::const_iterator pnItr = PNDigi->begin(); pnItr != PNDigi->end();
0294 ++pnItr) {
0295
0296 for (int samId = 0; samId < (*pnItr).size(); samId++) {
0297 pn[samId] = (*pnItr).sample(samId).adc();
0298 }
0299
0300 for (dsum = 0., k = 0; k < _presamplePN; k++) {
0301 dsum += pn[k];
0302 }
0303 bl = dsum / ((double)_presamplePN);
0304
0305 for (val_max = 0., k = 0; k < _nsamplesPN; k++) {
0306 ypnrange[k] = pn[k] - bl;
0307
0308 if (ypnrange[k] > val_max) {
0309 val_max = ypnrange[k];
0310 samplemax = k;
0311 }
0312 }
0313
0314 chi2pn = pnfit->doFit(samplemax, &ypnrange[0]);
0315
0316 if (chi2pn == 101 || chi2pn == 102 || chi2pn == 103)
0317 pnAmpl = 0.;
0318 else
0319 pnAmpl = pnfit->getAmpl();
0320
0321 allPNAmpl.push_back(pnAmpl);
0322 }
0323
0324
0325
0326
0327
0328 ttrig = -1.;
0329 peak = -1.;
0330
0331 if (doesRefFileExist == 1) {
0332
0333 if (color == 0)
0334 matacqTree->GetEntry(event - 1);
0335 else if (color == 3)
0336 matacqTree->GetEntry(matacqTree->GetEntries("color==0") + event - 1);
0337 ttrig = ttMat;
0338 peak = peakMat;
0339 }
0340
0341
0342
0343
0344
0345 double yrange[10];
0346 int adcGain = 0;
0347 int side = 0;
0348
0349 if (EBDigi) {
0350 for (EBDigiCollection::const_iterator digiItr = EBDigi->begin(); digiItr != EBDigi->end();
0351 ++digiItr) {
0352
0353 EBDetId id_crystal(digiItr->id());
0354 EBDataFrame df(*digiItr);
0355
0356 int etaG = id_crystal.ieta();
0357 int phiG = id_crystal.iphi();
0358
0359 int etaL;
0360 int phiL;
0361 std::pair<int, int> LocalCoord = MEEBGeom::localCoord(etaG, phiG);
0362
0363 etaL = LocalCoord.first;
0364 phiL = LocalCoord.second;
0365
0366 eta = etaG;
0367 phi = phiG;
0368
0369 side = MEEBGeom::side(etaG, phiG);
0370
0371 EcalElectronicsId elecid_crystal = TheMapping.getElectronicsId(id_crystal);
0372
0373 int towerID = elecid_crystal.towerId();
0374
0375 int strip = elecid_crystal.stripId();
0376 int xtal = elecid_crystal.xtalId();
0377 int channelID = 5 * (strip - 1) + xtal - 1;
0378
0379 int module = MEEBGeom::lmmod(etaG, phiG);
0380
0381 std::pair<int, int> pnpair = MEEBGeom::pn(module);
0382 unsigned int MyPn0 = pnpair.first;
0383 unsigned int MyPn1 = pnpair.second;
0384
0385 unsigned int channel = MEEBGeom::electronic_channel(etaL, phiL);
0386 assert(channel < nCrys);
0387
0388 double adcmax = 0.0;
0389
0390 if (towerID != int(_tower) || channelID != int(_channel) || dccID != int(_fedid - 600))
0391 continue;
0392 else
0393 channelNumber = channel;
0394
0395 for (unsigned int i = 0; i < (*digiItr).size(); ++i) {
0396
0397 EcalMGPASample samp_crystal(df.sample(i));
0398 adc[i] = samp_crystal.adc();
0399 adcG[i] = samp_crystal.gainId();
0400
0401 if (i == 0)
0402 adcGain = adcG[i];
0403 if (i > 0)
0404 adcGain = TMath::Max(adcG[i], adcGain);
0405
0406 if (adc[i] > adcmax) {
0407 adcmax = adc[i];
0408 }
0409 }
0410
0411 for (dsum = 0., dsum1 = 0., k = 0; k < _presample; k++) {
0412 dsum += adc[k];
0413 if (k < _presample - 1)
0414 dsum1 += adc[k];
0415 }
0416
0417 bl = dsum / ((double)_presample);
0418 bl1 = dsum1 / ((double)_presample - 1);
0419
0420 for (val_max = 0., k = 0; k < _nsamples; k++) {
0421 yrange[k] = adc[k] - bl;
0422 if (yrange[k] > val_max) {
0423 val_max = yrange[k];
0424 samplemax = k;
0425 }
0426 }
0427
0428 if (samplemax == 4 || samplemax == 5) {
0429 for (k = 0; k < _nsamples; k++) {
0430 yrange[k] = yrange[k] + bl - bl1;
0431 }
0432 }
0433
0434 for (unsigned int k = 0; k < _nsamples; k++) {
0435 adc[k] = yrange[k];
0436 }
0437
0438 pn0 = allPNAmpl[MyPn0];
0439 pn1 = allPNAmpl[MyPn1];
0440
0441 if (samplemax >= _timingcutlow && samplemax <= _timingcuthigh && lightside == side)
0442 ADCtrees->Fill();
0443 }
0444
0445 } else {
0446 for (EEDigiCollection::const_iterator digiItr = EEDigi->begin(); digiItr != EEDigi->end();
0447 ++digiItr) {
0448
0449 EEDetId id_crystal(digiItr->id());
0450 EEDataFrame df(*digiItr);
0451
0452 phi = id_crystal.ix() - 1;
0453 eta = id_crystal.iy() - 1;
0454 side = 0;
0455
0456
0457
0458 EcalElectronicsId elecid_crystal = TheMapping.getElectronicsId(id_crystal);
0459
0460 int towerID = elecid_crystal.towerId();
0461 int channelID = elecid_crystal.channelId() - 1;
0462
0463 int module = MEEEGeom::lmmod(phi, eta);
0464
0465 std::pair<int, int> pnpair = MEEEGeom::pn(module, _fedid);
0466 unsigned int MyPn0 = pnpair.first;
0467 unsigned int MyPn1 = pnpair.second;
0468
0469 unsigned int channel = MEEEGeom::crystal(phi, eta);
0470 assert(channel < nCrys);
0471
0472 double adcmax = 0.0;
0473
0474 if (towerID != int(_tower) || channelID != int(_channel) || dccID != int(_fedid - 600))
0475 continue;
0476 else
0477 channelNumber = channel;
0478
0479 for (unsigned int i = 0; i < (*digiItr).size(); ++i) {
0480
0481 EcalMGPASample samp_crystal(df.sample(i));
0482 adc[i] = samp_crystal.adc();
0483 adcG[i] = samp_crystal.gainId();
0484
0485 if (i == 0)
0486 adcGain = adcG[i];
0487 if (i > 0)
0488 adcGain = TMath::Max(adcG[i], adcGain);
0489
0490 if (adc[i] > adcmax) {
0491 adcmax = adc[i];
0492 }
0493 }
0494
0495 for (dsum = 0., dsum1 = 0., k = 0; k < _presample; k++) {
0496 dsum += adc[k];
0497 if (k < _presample - 1)
0498 dsum1 += adc[k];
0499 }
0500
0501 bl = dsum / ((double)_presample);
0502 bl1 = dsum1 / ((double)_presample - 1);
0503
0504 for (val_max = 0., k = 0; k < _nsamples; k++) {
0505 yrange[k] = adc[k] - bl;
0506 if (yrange[k] > val_max) {
0507 val_max = yrange[k];
0508 samplemax = k;
0509 }
0510 }
0511
0512 if (samplemax == 4 || samplemax == 5) {
0513 for (k = 0; k < _nsamples; k++) {
0514 yrange[k] = yrange[k] + bl - bl1;
0515 }
0516 }
0517
0518 for (unsigned int k = 0; k < _nsamples; k++) {
0519 adc[k] = yrange[k];
0520 }
0521
0522 pn0 = allPNAmpl[MyPn0];
0523 pn1 = allPNAmpl[MyPn1];
0524
0525 if (samplemax >= _timingcutlow && samplemax <= _timingcuthigh && lightside == side)
0526 ADCtrees->Fill();
0527 }
0528 }
0529
0530 }
0531
0532
0533 void EcalPerEvtLaserAnalyzer::endJob() {
0534
0535
0536 assert(colors.size() <= nColor);
0537 unsigned int nCol = colors.size();
0538
0539 ADCtrees->Write();
0540
0541 edm::LogVerbatim("EcalPerEvtLaserAnalyzer")
0542 << "\n\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+";
0543 edm::LogVerbatim("EcalPerEvtLaserAnalyzer")
0544 << "\t+=+ Analyzing laser data: getting per event +=+";
0545 edm::LogVerbatim("EcalPerEvtLaserAnalyzer")
0546 << "\t+=+ APD Amplitudes and ADC samples +=+";
0547 edm::LogVerbatim("EcalPerEvtLaserAnalyzer")
0548 << "\t+=+ for fed:" << _fedid << ", tower:" << _tower << ", and channel:" << _channel;
0549
0550
0551
0552 APDFile = new TFile(resfile.c_str(), "RECREATE");
0553
0554 int ieta, iphi, channelID, towerID, flag;
0555 double alpha, beta;
0556
0557 colors.push_back(color);
0558
0559 for (unsigned int i = 0; i < nCol; i++) {
0560 stringstream name1;
0561 name1 << "headerCol" << colors[i];
0562
0563 header[i] = new TTree(name1.str().c_str(), name1.str().c_str());
0564
0565 header[i]->Branch("alpha", &alpha, "alpha/D");
0566 header[i]->Branch("beta", &beta, "beta/D");
0567 header[i]->Branch("iphi", &iphi, "iphi/I");
0568 header[i]->Branch("ieta", &ieta, "ieta/I");
0569 header[i]->Branch("dccID", &dccID, "dccID/I");
0570 header[i]->Branch("towerID", &towerID, "towerID/I");
0571 header[i]->Branch("channelID", &channelID, "channelID/I");
0572
0573 header[i]->SetBranchAddress("alpha", &alpha);
0574 header[i]->SetBranchAddress("beta", &beta);
0575 header[i]->SetBranchAddress("iphi", &iphi);
0576 header[i]->SetBranchAddress("ieta", &ieta);
0577 header[i]->SetBranchAddress("dccID", &dccID);
0578 header[i]->SetBranchAddress("towerID", &towerID);
0579 header[i]->SetBranchAddress("channelID", &channelID);
0580 }
0581
0582 stringstream name2;
0583 name2 << "APDTree";
0584 APDtrees = new TTree(name2.str().c_str(), name2.str().c_str());
0585
0586
0587
0588 APDtrees->Branch("event", &event, "event/I");
0589 APDtrees->Branch("color", &color, "color/I");
0590 APDtrees->Branch("adc", &adc, "adc[10]/D");
0591 APDtrees->Branch("peakMatacq", &peak, "peak/D");
0592 APDtrees->Branch("ttrigMatacq", &ttrig, "ttrig/D");
0593 APDtrees->Branch("apdAmpl", &apdAmpl, "apdAmpl/D");
0594 APDtrees->Branch("apdTime", &apdTime, "apdTime/D");
0595 APDtrees->Branch("flag", &flag, "flag/I");
0596 APDtrees->Branch("pn0", &pn0, "pn0/D");
0597 APDtrees->Branch("pn1", &pn1, "pn1/D");
0598
0599 APDtrees->SetBranchAddress("event", &event);
0600 APDtrees->SetBranchAddress("color", &color);
0601 APDtrees->SetBranchAddress("adc", adc);
0602 APDtrees->SetBranchAddress("peakMatacq", &peak);
0603 APDtrees->SetBranchAddress("ttrigMatacq", &ttrig);
0604 APDtrees->SetBranchAddress("apdAmpl", &apdAmpl);
0605 APDtrees->SetBranchAddress("apdTime", &apdTime);
0606 APDtrees->SetBranchAddress("flag", &flag);
0607 APDtrees->SetBranchAddress("pn0", &pn0);
0608 APDtrees->SetBranchAddress("pn1", &pn1);
0609
0610
0611
0612 TFile* alphaFile = new TFile(refalphabeta_.c_str());
0613 TTree* alphaTree[2];
0614
0615 Double_t alphaRun, betaRun;
0616 int ietaRun, iphiRun, channelIDRun, towerIDRun, dccIDRun, flagRun;
0617
0618 for (unsigned int i = 0; i < nCol; i++) {
0619 stringstream name3;
0620 name3 << "ABCol" << i;
0621 alphaTree[i] = (TTree*)alphaFile->Get(name3.str().c_str());
0622 alphaTree[i]->SetBranchStatus("*", false);
0623 alphaTree[i]->SetBranchStatus("alpha", true);
0624 alphaTree[i]->SetBranchStatus("beta", true);
0625 alphaTree[i]->SetBranchStatus("iphi", true);
0626 alphaTree[i]->SetBranchStatus("ieta", true);
0627 alphaTree[i]->SetBranchStatus("dccID", true);
0628 alphaTree[i]->SetBranchStatus("towerID", true);
0629 alphaTree[i]->SetBranchStatus("channelID", true);
0630 alphaTree[i]->SetBranchStatus("flag", true);
0631
0632 alphaTree[i]->SetBranchAddress("alpha", &alphaRun);
0633 alphaTree[i]->SetBranchAddress("beta", &betaRun);
0634 alphaTree[i]->SetBranchAddress("iphi", &iphiRun);
0635 alphaTree[i]->SetBranchAddress("ieta", &ietaRun);
0636 alphaTree[i]->SetBranchAddress("dccID", &dccIDRun);
0637 alphaTree[i]->SetBranchAddress("towerID", &towerIDRun);
0638 alphaTree[i]->SetBranchAddress("channelID", &channelIDRun);
0639 alphaTree[i]->SetBranchAddress("flag", &flagRun);
0640 }
0641
0642 PulseFitWithFunction* pslsfit = new PulseFitWithFunction();
0643
0644 double chi2;
0645
0646 for (unsigned int icol = 0; icol < nCol; icol++) {
0647 IsThereDataADC[icol] = 1;
0648 stringstream cut;
0649 cut << "color==" << colors.at(icol);
0650 if (ADCtrees->GetEntries(cut.str().c_str()) < 10)
0651 IsThereDataADC[icol] = 0;
0652 IsHeaderFilled[icol] = 0;
0653 }
0654
0655
0656
0657 for (Long64_t jentry = 0; jentry < ADCtrees->GetEntriesFast(); jentry++) {
0658 ADCtrees->GetEntry(jentry);
0659
0660 int iCry = channelNumber;
0661
0662
0663
0664 unsigned int iCol = 0;
0665 for (unsigned int i = 0; i < nCol; i++) {
0666 if (color == colors[i]) {
0667 iCol = i;
0668 i = colors.size();
0669 }
0670 }
0671
0672 alphaTree[iCol]->GetEntry(iCry);
0673
0674 flag = flagRun;
0675 iphi = iphiRun;
0676 ieta = ietaRun;
0677 towerID = towerIDRun;
0678 channelID = channelIDRun;
0679 alpha = alphaRun;
0680 beta = betaRun;
0681
0682 if (IsHeaderFilled[iCol] == 0) {
0683 header[iCol]->Fill();
0684 IsHeaderFilled[iCol] = 1;
0685 }
0686
0687
0688 apdAmpl = 0;
0689 apdTime = 0;
0690
0691 pslsfit->init(_nsamples, _firstsample, _lastsample, _niter, alphaRun, betaRun);
0692 chi2 = pslsfit->doFit(&adc[0]);
0693
0694 if (chi2 < 0. || chi2 == 102 || chi2 == 101) {
0695 apdAmpl = 0;
0696 apdTime = 0;
0697
0698 } else {
0699 apdAmpl = pslsfit->getAmpl();
0700 apdTime = pslsfit->getTime();
0701 }
0702
0703 APDtrees->Fill();
0704 }
0705
0706 alphaFile->Close();
0707
0708 ADCFile->Close();
0709
0710 APDFile->Write();
0711 APDFile->Close();
0712
0713
0714
0715 stringstream del;
0716 del << "rm " << ADCfile;
0717 system(del.str().c_str());
0718
0719 edm::LogVerbatim("EcalPerEvtLaserAnalyzer")
0720 << "\t+=+ .................................................. done +=+";
0721 edm::LogVerbatim("EcalPerEvtLaserAnalyzer")
0722 << "\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+";
0723 }
0724
0725 DEFINE_FWK_MODULE(EcalPerEvtLaserAnalyzer);