RegionalCaloTrigger/plugins/MaskedRctInputDigiProducer.cc

0001 #include "L1Trigger/RegionalCaloTrigger/interface/MaskedRctInputDigiProducer.h"
0002
0003 #include <fstream>
0004 #include <iostream>
0005 #include <string>
0006 #include <vector>
0007 using std::endl;
0008 using std::ifstream;
0009 using std::vector;
0010
0011 //
0012 // constructors and destructor
0013 //
0014
0015 MaskedRctInputDigiProducer::MaskedRctInputDigiProducer(const edm::ParameterSet &iConfig)
0016     : useEcal_(iConfig.getParameter<bool>("useEcal")),
0017       useHcal_(iConfig.getParameter<bool>("useHcal")),
0018       maskFile_(iConfig.getParameter<edm::FileInPath>("maskFile")) {
0019   if (useEcal_) {
0020     ecalDigisToken_ = consumes(iConfig.getParameter<edm::InputTag>("ecalDigisLabel"));
0021   }
0022   if (useHcal_) {
0023     hcalDigisToken_ = consumes(iConfig.getParameter<edm::InputTag>("hcalDigisLabel"));
0024   }
0025   // register your products
0026   /* Examples
0027      produces<ExampleData2>();
0028
0029      //if do put with a label
0030      produces<ExampleData2>("label");
0031   */
0032
0033   produces<EcalTrigPrimDigiCollection>();
0034   produces<HcalTrigPrimDigiCollection>();
0035
0036   // now do what ever other initialization is needed
0037 }
0038
0039 MaskedRctInputDigiProducer::~MaskedRctInputDigiProducer() {
0040   // do anything here that needs to be done at desctruction time
0041   // (e.g. close files, deallocate resources etc.)
0042 }
0043
0044 //
0045 // member functions
0046 //
0047
0048 // ------------ method called to produce the data  ------------
0049 void MaskedRctInputDigiProducer::produce(edm::StreamID, edm::Event &iEvent, const edm::EventSetup &iSetup) const {
0050   using namespace edm;
0051
0052   edm::Handle<EcalTrigPrimDigiCollection> ecal;
0053   edm::Handle<HcalTrigPrimDigiCollection> hcal;
0054
0055   if (useEcal_) {
0056     ecal = iEvent.getHandle(ecalDigisToken_);
0057   }
0058   if (useHcal_) {
0059     hcal = iEvent.getHandle(hcalDigisToken_);
0060   }
0061
0062   EcalTrigPrimDigiCollection ecalColl;
0063   HcalTrigPrimDigiCollection hcalColl;
0064   if (ecal.isValid()) {
0065     ecalColl = *ecal;
0066   }
0067   if (hcal.isValid()) {
0068     hcalColl = *hcal;
0069   }
0070
0071   ifstream maskFileStream(maskFile_.fullPath().c_str());
0072   if (!maskFileStream.is_open()) {
0073     throw cms::Exception("FileInPathError") << "RCT mask file not opened" << std::endl;
0074     ;
0075   }
0076
0077   // read and process file (transform mask to new position in phi to match
0078   // digis)
0079   // char junk[256];
0080   std::string junk;
0081   // char junk[32];
0082   do {
0083     // maskFileStream.getline(junk, 256);
0084     maskFileStream >> junk;
0085   } while (junk != "ECAL:");
0086
0087   //   std::vector<std::vector<std::vector<unsigned short> > >
0088   //   temp(2,std::vector<std::vector<unsigned short> >(72,std::vector<unsigned
0089   //   short>(28,1)));
0090   std::vector<std::vector<std::vector<unsigned short>>> ecalMask(
0091       2, std::vector<std::vector<unsigned short>>(72, std::vector<unsigned short>(28, 1)));
0092   std::vector<std::vector<std::vector<unsigned short>>> hcalMask(
0093       2, std::vector<std::vector<unsigned short>>(72, std::vector<unsigned short>(28, 1)));
0094   std::vector<std::vector<std::vector<unsigned short>>> hfMask(
0095       2, std::vector<std::vector<unsigned short>>(18, std::vector<unsigned short>(8, 1)));
0096
0097   // read ECAL mask first
0098   // loop through rct iphi
0099   for (int i = 1; i <= 72; i++) {
0100     int phi_index = (72 + 18 - i) % 72;  // transfrom from RCT coords
0101     if (phi_index == 0) {
0102       phi_index = 72;
0103     }
0104     // std::cout << "ecal phi index is " << phi_index << std::endl;
0105     for (int j = 28; j >= 1; j--) {
0106       maskFileStream >> junk;
0107       if (junk == "-") {
0108       } else if ((junk == "X") || (junk == "x")) {
0109         ecalMask.at(0).at(phi_index - 1).at(j - 1) = 0;
0110       } else {
0111         std::cerr << "RCT mask producer: error -- unrecognized character" << std::endl;
0112       }
0113     }
0114     for (int j = 1; j <= 28; j++) {
0115       maskFileStream >> junk;
0116       if (junk == "-") {
0117       } else if ((junk == "X") || (junk == "x")) {
0118         ecalMask.at(1).at(phi_index - 1).at(j - 1) = 0;
0119       } else {
0120         std::cerr << "RCT mask producer: error -- unrecognized character" << std::endl;
0121       }
0122     }
0123   }
0124   // std::cout << "done reading ECAL" << std::endl;
0125
0126   maskFileStream >> junk;
0127   if (junk != "HCAL:") {
0128     throw cms::Exception("FileInPathError") << "RCT mask producer: error reading ECAL mask" << std::endl;
0129     ;
0130   }
0131
0132   // std::cout << "Starting HCAL read" << std::endl;
0133
0134   // now read HCAL mask
0135   // loop through rct iphi
0136   for (int i = 1; i <= 72; i++) {
0137     int phi_index = (72 + 18 - i) % 72;  // transfrom from RCT coords
0138     if (phi_index == 0) {
0139       phi_index = 72;
0140     }
0141     // std::cout << "hcal phi index is " << phi_index << std::endl;
0142     for (int j = 28; j >= 1; j--) {
0143       maskFileStream >> junk;
0144       if (junk == "-") {
0145       } else if ((junk == "X") || (junk == "x")) {
0146         hcalMask.at(0).at(phi_index - 1).at(j - 1) = 0;
0147       } else {
0148         std::cerr << "RCT mask producer: error -- unrecognized character" << std::endl;
0149       }
0150     }
0151     for (int j = 1; j <= 28; j++) {
0152       maskFileStream >> junk;
0153       if (junk == "-") {
0154       } else if ((junk == "X") || (junk == "x")) {
0155         hcalMask.at(1).at(phi_index - 1).at(j - 1) = 0;
0156       } else {
0157         std::cerr << "RCT mask producer: error -- unrecognized character" << std::endl;
0158       }
0159     }
0160   }
0161
0162   maskFileStream >> junk;
0163   if (junk != "HF:") {
0164     throw cms::Exception("FileInPathError") << "RCT mask producer: error reading HCAL mask" << std::endl;
0165     ;
0166   }
0167
0168   // loop through the hf phi columns in file
0169   for (int i = 0; i < 18; i++) {
0170     // std::cout << "iphi for hf file read is " << i << std::endl;
0171     for (int j = 4; j >= 1; j--) {
0172       if (maskFileStream >> junk) {
0173       } else {
0174         std::cerr << "RCT mask producer: error reading HF mask" << std::endl;
0175       }
0176       if (junk == "-") {
0177       } else if ((junk == "X") || (junk == "x")) {
0178         hfMask.at(0).at(i).at(j - 1) = 0;  // just save iphi as 0-17, transform later
0179       } else {
0180         std::cerr << "RCT mask producer: error -- unrecognized character" << std::endl;
0181       }
0182     }
0183     for (int j = 1; j <= 4; j++) {
0184       if (maskFileStream >> junk) {
0185       } else {
0186         std::cerr << "RCT mask producer: error reading HF mask" << std::endl;
0187       }
0188       if (junk == "-") {
0189       } else if ((junk == "X") || (junk == "x")) {
0190         hfMask.at(1).at(i).at(j - 1) = 0;
0191       } else {
0192         std::cerr << "RCT mask producer: error -- unrecognized character" << std::endl;
0193       }
0194     }
0195   }
0196
0197   // std::cout << "HF read done" << std::endl;
0198
0199   maskFileStream.close();
0200
0201   // std::cout << "file closed" << std::endl;
0202
0203   // apply mask
0204
0205   std::unique_ptr<EcalTrigPrimDigiCollection> maskedEcalTPs(new EcalTrigPrimDigiCollection());
0206   std::unique_ptr<HcalTrigPrimDigiCollection> maskedHcalTPs(new HcalTrigPrimDigiCollection());
0207   maskedEcalTPs->reserve(56 * 72);
0208   maskedHcalTPs->reserve(56 * 72 + 18 * 8);
0209   int nEcalSamples = 0;
0210   int nHcalSamples = 0;
0211
0212   for (unsigned int i = 0; i < ecalColl.size(); i++) {
0213     nEcalSamples = ecalColl[i].size();
0214     short ieta = (short)ecalColl[i].id().ieta();
0215     unsigned short absIeta = (unsigned short)abs(ieta);
0216     int sign = ieta / absIeta;
0217     short iphi = (unsigned short)ecalColl[i].id().iphi();
0218     // if (i < 20) {std::cout << "ieta is " << ieta << ", absIeta is " <<
0219     // absIeta
0220     //            << ", iphi is " << iphi << std::endl;}
0221
0222     EcalTriggerPrimitiveDigi ecalDigi(EcalTrigTowerDetId(sign, EcalTriggerTower, absIeta, iphi));
0223     ecalDigi.setSize(nEcalSamples);
0224
0225     for (int nSample = 0; nSample < nEcalSamples; nSample++) {
0226       int energy = 0;
0227       bool fineGrain = false;
0228
0229       if (sign < 0) {
0230         // std::cout << "eta-: mask is " <<
0231         // ecalMask.at(0).at(iphi-1).at(absIeta-1) << std::endl;
0232         energy = ecalMask.at(0).at(iphi - 1).at(absIeta - 1) * ecalColl[i].sample(nSample).compressedEt();
0233         fineGrain = (ecalMask.at(0).at(iphi - 1).at(absIeta - 1) != 0) && ecalColl[i].sample(nSample).fineGrain();
0234       } else if (sign > 0) {
0235         // std::cout << "eta+: mask is " <<
0236         // ecalMask.at(1).at(iphi-1).at(absIeta-1) << std::endl;
0237         energy = ecalMask.at(1).at(iphi - 1).at(absIeta - 1) * ecalColl[i].sample(nSample).compressedEt();
0238         fineGrain = (ecalMask.at(1).at(iphi - 1).at(absIeta - 1) != 0) && ecalColl[i].sample(nSample).fineGrain();
0239       }
0240
0241       ecalDigi.setSample(nSample, EcalTriggerPrimitiveSample(energy, fineGrain, 0));
0242     }
0243     maskedEcalTPs->push_back(ecalDigi);
0244   }
0245   // std::cout << "End of ecal digi masking" << std::endl;
0246
0247   // std::cout << "nHcalDigis is " << hcalColl.size() << std::endl;
0248   for (unsigned int i = 0; i < hcalColl.size(); i++) {
0249     nHcalSamples = hcalColl[i].size();
0250     // if ((i % 100) == 0 ) {std::cout << "Loop " << i << std::endl;}
0251     short ieta = (short)hcalColl[i].id().ieta();
0252     unsigned short absIeta = (unsigned short)abs(ieta);
0253     int sign = ieta / absIeta;
0254     short iphi = (unsigned short)hcalColl[i].id().iphi();
0255     // if (i < 20) {std::cout << "ieta is " << ieta << ", absIeta is " <<
0256     // absIeta
0257     //      << ", iphi is " << iphi << std::endl;}
0258     /*
0259     if (hcalColl[i].SOI_compressedEt() != 0)
0260       {
0261         std::cout << "original et " << hcalColl[i].SOI_compressedEt()
0262                   << "  fg " << hcalColl[i].SOI_fineGrain() << "  iphi "
0263                   << iphi << "  ieta " << ieta << std::endl;
0264       }
0265     */
0266     HcalTriggerPrimitiveDigi hcalDigi(HcalTrigTowerDetId(ieta, iphi));
0267     hcalDigi.setSize(nHcalSamples);
0268     hcalDigi.setPresamples(hcalColl[i].presamples());
0269
0270     for (int nSample = 0; nSample < nHcalSamples; nSample++) {
0271       int energy = 0;
0272       bool fineGrain = false;
0273
0274       if (absIeta < 29) {
0275         if (sign < 0) {
0276           energy = hcalMask.at(0).at(iphi - 1).at(absIeta - 1) * hcalColl[i].sample(nSample).compressedEt();
0277           fineGrain = (hcalMask.at(0).at(iphi - 1).at(absIeta - 1) != 0) && hcalColl[i].sample(nSample).fineGrain();
0278         } else if (sign > 0) {
0279           energy = hcalMask.at(1).at(iphi - 1).at(absIeta - 1) * hcalColl[i].sample(nSample).compressedEt();
0280           fineGrain = (hcalMask.at(1).at(iphi - 1).at(absIeta - 1) != 0) && hcalColl[i].sample(nSample).fineGrain();
0281         }
0282       } else if ((absIeta >= 29) && (absIeta <= 32)) {
0283         // std::cout << "hf iphi: " << iphi << std::endl;
0284         short hf_phi_index = iphi / 4;
0285         // iphi = iphi/4;  // change from 1,5,9, etc to access vector positions
0286         // std::cout << "hf phi index: " << hf_phi_index << std::endl;
0287         if (sign < 0) {
0288           // std::cout << "ieta is " << ieta << ", absIeta is " << absIeta << ",
0289           // iphi is " << iphi << std::endl; std::cout << "eta-: mask is " <<
0290           // hfMask.at(0).at(hf_phi_index).at(absIeta-29) << std::endl; // hf
0291           // ieta 0-3
0292           energy = hfMask.at(0).at(hf_phi_index).at(absIeta - 29) *
0293                    hcalColl[i].sample(nSample).compressedEt();  // for hf, iphi starts at 0
0294           fineGrain = (hfMask.at(0).at(hf_phi_index).at(absIeta - 29) != 0) && hcalColl[i].sample(nSample).fineGrain();
0295         } else if (sign > 0) {
0296           // std::cout << "ieta is " << ieta << ", absIeta is " << absIeta << ",
0297           // iphi is " << iphi << std::endl; std::cout << "eta+: mask is " <<
0298           // hfMask.at(1).at(hf_phi_index).at(absIeta-29) << std::endl;
0299           energy = hfMask.at(1).at(hf_phi_index).at(absIeta - 29) * hcalColl[i].sample(nSample).compressedEt();
0300           fineGrain = (hfMask.at(1).at(hf_phi_index).at(absIeta - 29) != 0) && hcalColl[i].sample(nSample).fineGrain();
0301         }
0302         // iphi = iphi*4 + 1; // change back to original
0303         // std::cout << "New hf iphi = " << iphi << std::endl;
0304       }
0305
0306       hcalDigi.setSample(nSample, HcalTriggerPrimitiveSample(energy, fineGrain, 0, 0));
0307
0308       // if (hcalDigi.SOI_compressedEt() != 0)
0309       //{
0310       //  std::cout << "et " << hcalDigi.SOI_compressedEt()
0311       //     << "fg " << hcalDigi.SOI_fineGrain() << std::endl;
0312       //}
0313     }
0314     maskedHcalTPs->push_back(hcalDigi);
0315   }
0316   // std::cout << "End of hcal digi masking" << std::endl;
0317
0318   // put new data into event
0319
0320   iEvent.put(std::move(maskedEcalTPs));
0321   iEvent.put(std::move(maskedHcalTPs));
0322 }
0323
0324 // define this as a plug-in
0325 DEFINE_FWK_MODULE(MaskedRctInputDigiProducer);