CSCCalibration/plugins/CSCGainsConditions.cc

0001 #include <fstream>
0002 #include <iostream>
0003
0004 #include "CalibMuon/CSCCalibration/interface/CSCGainsConditions.h"
0005
0006 CSCGains *CSCGainsConditions::prefillGains() {
0007   float mean, min, minchi;
0008   int seed;
0009   int old_chamber_id, old_strip, new_chamber_id, new_strip;
0010   float old_gainslope, old_intercpt, old_chisq;
0011   std::vector<int> old_cham_id;
0012   std::vector<int> old_strips;
0013   std::vector<float> old_slope;
0014   std::vector<float> old_intercept;
0015   std::vector<float> old_chi2;
0016   float new_gainslope, new_intercpt, new_chisq;
0017   std::vector<int> new_cham_id;
0018   std::vector<int> new_strips;
0019   std::vector<float> new_slope;
0020   std::vector<float> new_intercept;
0021   std::vector<float> new_chi2;
0022
0023   const CSCDetId &detId = CSCDetId();
0024   CSCGains *cngains = new CSCGains();
0025
0026   int max_istrip, id_layer, max_ring, max_cham;
0027   unsigned int old_nrlines = 0;
0028   unsigned int new_nrlines = 0;
0029   seed = 10000;
0030   srand(seed);
0031   mean = 6.8, min = -10.0, minchi = 1.0;
0032
0033   std::ifstream olddata;
0034   olddata.open("old_gains.dat", std::ios::in);
0035   if (!olddata) {
0036     std::cerr << "Error: old_gains.dat -> no such file!" << std::endl;
0037     exit(1);
0038   }
0039
0040   while (!olddata.eof()) {
0041     olddata >> old_chamber_id >> old_strip >> old_gainslope >> old_intercpt >> old_chisq;
0042     old_cham_id.push_back(old_chamber_id);
0043     old_strips.push_back(old_strip);
0044     old_slope.push_back(old_gainslope);
0045     old_intercept.push_back(old_intercpt);
0046     old_chi2.push_back(old_chisq);
0047     old_nrlines++;
0048   }
0049   olddata.close();
0050
0051   std::ifstream newdata;
0052   newdata.open("new_gains.txt", std::ios::in);
0053   if (!newdata) {
0054     std::cerr << "Error: new_gains.txt -> no such file!" << std::endl;
0055     exit(1);
0056   }
0057
0058   while (!newdata.eof()) {
0059     newdata >> new_chamber_id >> new_strip >> new_gainslope >> new_intercpt >> new_chisq;
0060     new_cham_id.push_back(new_chamber_id);
0061     new_strips.push_back(new_strip);
0062     new_slope.push_back(new_gainslope);
0063     new_intercept.push_back(new_intercpt);
0064     new_chi2.push_back(new_chisq);
0065     new_nrlines++;
0066   }
0067   newdata.close();
0068
0069   // endcap=1 to 2,station=1 to 4, ring=1 to 4,chamber=1 to 36,layer=1 to 6
0070   for (int iendcap = detId.minEndcapId(); iendcap <= detId.maxEndcapId(); iendcap++) {
0071     for (int istation = detId.minStationId(); istation <= detId.maxStationId(); istation++) {
0072       max_ring = detId.maxRingId();
0073       // station 4 ring 4 not there(36 chambers*2 missing)
0074       // 3 rings max this way of counting (ME1a & b)
0075       if (istation == 1)
0076         max_ring = 3;
0077       if (istation == 2)
0078         max_ring = 2;
0079       if (istation == 3)
0080         max_ring = 2;
0081       if (istation == 4)
0082         max_ring = 1;
0083
0084       for (int iring = detId.minRingId(); iring <= max_ring; iring++) {
0085         max_istrip = 80;
0086         max_cham = detId.maxChamberId();
0087         if (istation == 1 && iring == 1)
0088           max_cham = 36;
0089         if (istation == 1 && iring == 2)
0090           max_cham = 36;
0091         if (istation == 1 && iring == 3)
0092           max_cham = 36;
0093         if (istation == 2 && iring == 1)
0094           max_cham = 18;
0095         if (istation == 2 && iring == 2)
0096           max_cham = 36;
0097         if (istation == 3 && iring == 1)
0098           max_cham = 18;
0099         if (istation == 3 && iring == 2)
0100           max_cham = 36;
0101         if (istation == 4 && iring == 1)
0102           max_cham = 18;
0103
0104         for (int ichamber = detId.minChamberId(); ichamber <= max_cham; ichamber++) {
0105           for (int ilayer = detId.minLayerId(); ilayer <= detId.maxLayerId(); ilayer++) {
0106             // station 1 ring 3 has 64 strips per layer instead of 80
0107             if (istation == 1 && iring == 3)
0108               max_istrip = 64;
0109
0110             std::vector<CSCGains::Item> itemvector;
0111             itemvector.resize(max_istrip);
0112             id_layer = 100000 * iendcap + 10000 * istation + 1000 * iring + 10 * ichamber + ilayer;
0113
0114             for (int istrip = 0; istrip < max_istrip; istrip++) {
0115               itemvector[istrip].gain_slope = ((double)rand() / ((double)(RAND_MAX) + (double)(1))) + mean;
0116               itemvector[istrip].gain_intercept = ((double)rand() / ((double)(RAND_MAX) + (double)(1))) + min;
0117               itemvector[istrip].gain_chi2 = ((double)rand() / ((double)(RAND_MAX) + (double)(1))) + minchi;
0118               cngains->gains[id_layer] = itemvector;
0119             }
0120           }
0121         }
0122       }
0123     }
0124   }
0125
0126   // overwrite fakes with old values from DB
0127   int istrip = 0;
0128   std::vector<CSCGains::Item> itemvector;
0129   itemvector.resize(80);
0130
0131   for (unsigned int mystrip = 0; mystrip < old_nrlines - 1; mystrip++) {
0132     if (old_strips[mystrip] == 0)
0133       istrip = 0;
0134     itemvector[istrip].gain_slope = old_slope[mystrip];
0135     itemvector[istrip].gain_intercept = old_intercept[mystrip];
0136     itemvector[istrip].gain_chi2 = old_chi2[mystrip];
0137     cngains->gains[old_cham_id[mystrip]] = itemvector;
0138     istrip++;
0139   }
0140
0141   itemvector.resize(64);
0142   for (unsigned int mystrip = 0; mystrip < old_nrlines - 1; mystrip++) {
0143     if (old_strips[mystrip] == 0)
0144       istrip = 0;
0145     if (old_cham_id[mystrip] >= 113000 && old_cham_id[mystrip] <= 113999) {
0146       itemvector[istrip].gain_slope = old_slope[mystrip];
0147       itemvector[istrip].gain_intercept = old_intercept[mystrip];
0148       itemvector[istrip].gain_chi2 = old_chi2[mystrip];
0149       cngains->gains[old_cham_id[mystrip]] = itemvector;
0150       istrip++;
0151     }
0152   }
0153
0154   itemvector.resize(64);
0155   for (unsigned int mystrip = 0; mystrip < old_nrlines - 1; mystrip++) {
0156     if (old_strips[mystrip] == 0)
0157       istrip = 0;
0158     if (old_cham_id[mystrip] >= 213000 && old_cham_id[mystrip] <= 213999) {
0159       itemvector[istrip].gain_slope = old_slope[mystrip];
0160       itemvector[istrip].gain_intercept = old_intercept[mystrip];
0161       itemvector[istrip].gain_chi2 = old_chi2[mystrip];
0162       cngains->gains[old_cham_id[mystrip]] = itemvector;
0163       istrip++;
0164     }
0165   }
0166
0167   // overwrite old values with ones from new runs
0168   itemvector.resize(80);
0169   for (unsigned int mystrip = 0; mystrip < new_nrlines - 1; mystrip++) {
0170     if (new_strips[mystrip] == 0)
0171       istrip = 0;
0172     itemvector[istrip].gain_slope = new_slope[mystrip];
0173     itemvector[istrip].gain_intercept = new_intercept[mystrip];
0174     itemvector[istrip].gain_chi2 = new_chi2[mystrip];
0175     cngains->gains[new_cham_id[mystrip]] = itemvector;
0176     istrip++;
0177   }
0178
0179   itemvector.resize(64);
0180   for (unsigned int mystrip = 0; mystrip < new_nrlines - 1; mystrip++) {
0181     if (new_strips[mystrip] == 0)
0182       istrip = 0;
0183     if (new_cham_id[mystrip] >= 113000 && new_cham_id[mystrip] <= 113999) {
0184       itemvector[istrip].gain_slope = new_slope[mystrip];
0185       itemvector[istrip].gain_intercept = new_intercept[mystrip];
0186       itemvector[istrip].gain_chi2 = new_chi2[mystrip];
0187       cngains->gains[new_cham_id[mystrip]] = itemvector;
0188       istrip++;
0189     }
0190   }
0191
0192   itemvector.resize(64);
0193   for (unsigned int mystrip = 0; mystrip < new_nrlines - 1; mystrip++) {
0194     if (new_strips[mystrip] == 0)
0195       istrip = 0;
0196     if (new_cham_id[mystrip] >= 213000 && new_cham_id[mystrip] <= 213999) {
0197       itemvector[istrip].gain_slope = new_slope[mystrip];
0198       itemvector[istrip].gain_intercept = new_intercept[mystrip];
0199       itemvector[istrip].gain_chi2 = new_chi2[mystrip];
0200       cngains->gains[new_cham_id[mystrip]] = itemvector;
0201       istrip++;
0202     }
0203   }
0204   return cngains;
0205 }
0206
0207 CSCGainsConditions::CSCGainsConditions(const edm::ParameterSet &iConfig) {
0208   // the following line is needed to tell the framework what
0209   // added by Zhen (changed since 1_2_0)
0210   setWhatProduced(this, &CSCGainsConditions::produceGains);
0211   findingRecord<CSCGainsRcd>();
0212   // now do what ever other initialization is needed
0213 }
0214
0215 CSCGainsConditions::~CSCGainsConditions() {
0216   // do anything here that needs to be done at desctruction time
0217   // (e.g. close files, deallocate resources etc.)
0218 }
0219
0220 //
0221 // member functions
0222 //
0223
0224 // ------------ method called to produce the data  ------------
0225 CSCGainsConditions::ReturnType CSCGainsConditions::produceGains(const CSCGainsRcd &iRecord) {
0226   // Added by Zhen, need a new object so to not be deleted at exit
0227   return CSCGainsConditions::ReturnType(prefillGains());
0228 }
0229
0230 void CSCGainsConditions::setIntervalFor(const edm::eventsetup::EventSetupRecordKey &,
0231                                         const edm::IOVSyncValue &,
0232                                         edm::ValidityInterval &oValidity) {
0233   oValidity = edm::ValidityInterval(edm::IOVSyncValue::beginOfTime(), edm::IOVSyncValue::endOfTime());
0234 }