CSCCalibration/plugins/CSCCrosstalkConditions.cc

0001 #include <fstream>
0002 #include <iostream>
0003
0004 #include "CalibMuon/CSCCalibration/interface/CSCCrosstalkConditions.h"
0005
0006 CSCcrosstalk *CSCCrosstalkConditions::prefillCrosstalk() {
0007   float mean, min, minchi;
0008   int seed;
0009   int old_chamber_id, old_strip, new_chamber_id, new_strip;
0010   float old_slope_right, old_slope_left, old_intercept_right;
0011   float old_intercept_left, old_chi2_right, old_chi2_left;
0012   std::vector<int> old_cham_id;
0013   std::vector<int> old_strips;
0014   std::vector<float> old_slope_r;
0015   std::vector<float> old_intercept_r;
0016   std::vector<float> old_chi2_r;
0017   std::vector<float> old_slope_l;
0018   std::vector<float> old_intercept_l;
0019   std::vector<float> old_chi2_l;
0020   float new_slope_right, new_slope_left, new_intercept_right;
0021   float new_intercept_left, new_chi2_right, new_chi2_left;
0022   std::vector<int> new_cham_id;
0023   std::vector<int> new_strips;
0024   std::vector<float> new_slope_r;
0025   std::vector<float> new_intercept_r;
0026   std::vector<float> new_chi2_r;
0027   std::vector<float> new_slope_l;
0028   std::vector<float> new_intercept_l;
0029   std::vector<float> new_chi2_l;
0030
0031   const CSCDetId &detId = CSCDetId();
0032   CSCcrosstalk *cncrosstalk = new CSCcrosstalk();
0033
0034   int max_istrip, id_layer, max_ring, max_cham;
0035   unsigned int old_nrlines = 0;
0036   unsigned int new_nrlines = 0;
0037   seed = 10000;
0038   srand(seed);
0039   mean = -0.0009, min = 0.035, minchi = 1.5;
0040
0041   // endcap=1 to 2,station=1 to 4, ring=1 to 4,chamber=1 to 36,layer=1 to 6
0042   std::ifstream olddata;
0043   olddata.open("old_xtalk.dat", std::ios::in);
0044   if (!olddata) {
0045     std::cerr << "Error: old_xtalk.dat -> no such file!" << std::endl;
0046     exit(1);
0047   }
0048
0049   while (!olddata.eof()) {
0050     olddata >> old_chamber_id >> old_strip >> old_slope_right >> old_intercept_right >> old_chi2_right >>
0051         old_slope_left >> old_intercept_left >> old_chi2_left;
0052     old_cham_id.push_back(old_chamber_id);
0053     old_strips.push_back(old_strip);
0054     old_slope_r.push_back(old_slope_right);
0055     old_slope_l.push_back(old_slope_left);
0056     old_intercept_r.push_back(old_intercept_right);
0057     old_intercept_l.push_back(old_intercept_left);
0058     old_chi2_r.push_back(old_chi2_right);
0059     old_chi2_l.push_back(old_chi2_left);
0060     old_nrlines++;
0061   }
0062   olddata.close();
0063
0064   std::ifstream newdata;
0065   newdata.open("new_xtalk.txt", std::ios::in);
0066   if (!newdata) {
0067     std::cerr << "Error: new_xtalk.dat -> no such file!" << std::endl;
0068     exit(1);
0069   }
0070
0071   while (!newdata.eof()) {
0072     newdata >> new_chamber_id >> new_strip >> new_slope_right >> new_intercept_right >> new_chi2_right >>
0073         new_slope_left >> new_intercept_left >> new_chi2_left;
0074     new_cham_id.push_back(new_chamber_id);
0075     new_strips.push_back(new_strip);
0076     new_slope_r.push_back(new_slope_right);
0077     new_slope_l.push_back(new_slope_left);
0078     new_intercept_r.push_back(new_intercept_right);
0079     new_intercept_l.push_back(new_intercept_left);
0080     new_chi2_r.push_back(new_chi2_right);
0081     new_chi2_l.push_back(new_chi2_left);
0082     new_nrlines++;
0083   }
0084   newdata.close();
0085
0086   for (int iendcap = detId.minEndcapId(); iendcap <= detId.maxEndcapId(); iendcap++) {
0087     for (int istation = detId.minStationId(); istation <= detId.maxStationId(); istation++) {
0088       max_ring = detId.maxRingId();
0089       // station 4 ring 4 not there(36 chambers*2 missing)
0090       // 3 rings max this way of counting (ME1a & b)
0091       if (istation == 1)
0092         max_ring = 3;
0093       if (istation == 2)
0094         max_ring = 2;
0095       if (istation == 3)
0096         max_ring = 2;
0097       if (istation == 4)
0098         max_ring = 1;
0099
0100       for (int iring = detId.minRingId(); iring <= max_ring; iring++) {
0101         max_istrip = 80;
0102         max_cham = detId.maxChamberId();
0103         if (istation == 1 && iring == 1)
0104           max_cham = 36;
0105         if (istation == 1 && iring == 2)
0106           max_cham = 36;
0107         if (istation == 1 && iring == 3)
0108           max_cham = 36;
0109         if (istation == 2 && iring == 1)
0110           max_cham = 18;
0111         if (istation == 2 && iring == 2)
0112           max_cham = 36;
0113         if (istation == 3 && iring == 1)
0114           max_cham = 18;
0115         if (istation == 3 && iring == 2)
0116           max_cham = 36;
0117         if (istation == 4 && iring == 1)
0118           max_cham = 18;
0119         // station 1 ring 3 has 64 strips per layer instead of 80(minus & plus
0120         // side!!!)
0121
0122         for (int ichamber = detId.minChamberId(); ichamber <= max_cham; ichamber++) {
0123           for (int ilayer = detId.minLayerId(); ilayer <= detId.maxLayerId(); ilayer++) {
0124             // station 1 ring 3 has 64 strips per layer instead of 80
0125             if (istation == 1 && iring == 3)
0126               max_istrip = 64;
0127
0128             std::vector<CSCcrosstalk::Item> itemvector;
0129             itemvector.resize(max_istrip);
0130             id_layer = 100000 * iendcap + 10000 * istation + 1000 * iring + 10 * ichamber + ilayer;
0131
0132             for (int istrip = 0; istrip < max_istrip; istrip++) {
0133               // create fake values
0134               itemvector[istrip].xtalk_slope_right =
0135                   -((double)rand() / ((double)(RAND_MAX) + (double)(1))) / 10000 + mean;
0136               itemvector[istrip].xtalk_intercept_right =
0137                   ((double)rand() / ((double)(RAND_MAX) + (double)(1))) / 100 + min;
0138               itemvector[istrip].xtalk_chi2_right = ((double)rand() / ((double)(RAND_MAX) + (double)(1))) + minchi;
0139               itemvector[istrip].xtalk_slope_left =
0140                   -((double)rand() / ((double)(RAND_MAX) + (double)(1))) / 10000 + mean;
0141               itemvector[istrip].xtalk_intercept_left =
0142                   ((double)rand() / ((double)(RAND_MAX) + (double)(1))) / 100 + min;
0143               itemvector[istrip].xtalk_chi2_left = ((double)rand() / ((double)(RAND_MAX) + (double)(1))) + minchi;
0144               cncrosstalk->crosstalk[id_layer] = itemvector;
0145
0146               if (istrip == 0) {
0147                 itemvector[istrip].xtalk_slope_right =
0148                     -((double)rand() / ((double)(RAND_MAX) + (double)(1))) / 10000 + mean;
0149                 itemvector[istrip].xtalk_intercept_right =
0150                     ((double)rand() / ((double)(RAND_MAX) + (double)(1))) / 100 + min;
0151                 itemvector[istrip].xtalk_chi2_right = ((double)rand() / ((double)(RAND_MAX) + (double)(1))) + minchi;
0152                 itemvector[istrip].xtalk_slope_left = 0.0;
0153                 itemvector[istrip].xtalk_intercept_left = 0.0;
0154                 itemvector[istrip].xtalk_chi2_left = 0.0;
0155                 cncrosstalk->crosstalk[id_layer] = itemvector;
0156               }
0157
0158               if (istrip == 79) {
0159                 itemvector[istrip].xtalk_slope_right = 0.0;
0160                 itemvector[istrip].xtalk_intercept_right = 0.0;
0161                 itemvector[istrip].xtalk_chi2_right = 0.0;
0162                 itemvector[istrip].xtalk_slope_left =
0163                     -((double)rand() / ((double)(RAND_MAX) + (double)(1))) / 10000 + mean;
0164                 itemvector[istrip].xtalk_intercept_left =
0165                     ((double)rand() / ((double)(RAND_MAX) + (double)(1))) / 100 + min;
0166                 itemvector[istrip].xtalk_chi2_left = ((double)rand() / ((double)(RAND_MAX) + (double)(1))) + minchi;
0167                 cncrosstalk->crosstalk[id_layer] = itemvector;
0168               }
0169             }
0170           }
0171         }
0172       }
0173     }
0174
0175     // overwrite fakes with old values from DB
0176     int istrip = 0;
0177     std::vector<CSCcrosstalk::Item> itemvector;
0178     itemvector.resize(80);
0179
0180     for (unsigned int mystrip = 0; mystrip < old_nrlines - 1; mystrip++) {
0181       if (old_strips[mystrip] == 0)
0182         istrip = 0;
0183       itemvector[istrip].xtalk_slope_right = old_slope_r[mystrip];
0184       itemvector[istrip].xtalk_intercept_right = old_intercept_r[mystrip];
0185       itemvector[istrip].xtalk_chi2_right = old_chi2_r[mystrip];
0186       itemvector[istrip].xtalk_slope_left = old_slope_l[mystrip];
0187       itemvector[istrip].xtalk_intercept_left = old_intercept_l[mystrip];
0188       itemvector[istrip].xtalk_chi2_left = old_chi2_l[mystrip];
0189       cncrosstalk->crosstalk[old_cham_id[mystrip]] = itemvector;
0190       istrip++;
0191     }
0192
0193     itemvector.resize(64);
0194     for (unsigned int mystrip = 0; mystrip < old_nrlines - 1; mystrip++) {
0195       if (old_strips[mystrip] == 0)
0196         istrip = 0;
0197       if (old_cham_id[mystrip] >= 113000 && old_cham_id[mystrip] <= 113999) {
0198         itemvector[istrip].xtalk_slope_right = old_slope_r[mystrip];
0199         itemvector[istrip].xtalk_intercept_right = old_intercept_r[mystrip];
0200         itemvector[istrip].xtalk_chi2_right = old_chi2_r[mystrip];
0201         itemvector[istrip].xtalk_slope_left = old_slope_l[mystrip];
0202         itemvector[istrip].xtalk_intercept_left = old_intercept_l[mystrip];
0203         itemvector[istrip].xtalk_chi2_left = old_chi2_l[mystrip];
0204         cncrosstalk->crosstalk[old_cham_id[mystrip]] = itemvector;
0205         istrip++;
0206       }
0207     }
0208
0209     itemvector.resize(64);
0210     for (unsigned int mystrip = 0; mystrip < old_nrlines - 1; mystrip++) {
0211       if (old_strips[mystrip] == 0)
0212         istrip = 0;
0213       if (old_cham_id[mystrip] >= 213000 && old_cham_id[mystrip] <= 213999) {
0214         itemvector[istrip].xtalk_slope_right = old_slope_r[mystrip];
0215         itemvector[istrip].xtalk_intercept_right = old_intercept_r[mystrip];
0216         itemvector[istrip].xtalk_chi2_right = old_chi2_r[mystrip];
0217         itemvector[istrip].xtalk_slope_left = old_slope_l[mystrip];
0218         itemvector[istrip].xtalk_intercept_left = old_intercept_l[mystrip];
0219         itemvector[istrip].xtalk_chi2_left = old_chi2_l[mystrip];
0220         cncrosstalk->crosstalk[old_cham_id[mystrip]] = itemvector;
0221         istrip++;
0222       }
0223     }
0224
0225     // overwrite old values with ones from new runs
0226     itemvector.resize(80);
0227     for (unsigned int mystrip = 0; mystrip < new_nrlines - 1; mystrip++) {
0228       if (new_strips[mystrip] == 0)
0229         istrip = 0;
0230       itemvector[istrip].xtalk_slope_right = new_slope_r[mystrip];
0231       itemvector[istrip].xtalk_intercept_right = new_intercept_r[mystrip];
0232       itemvector[istrip].xtalk_chi2_right = new_chi2_r[mystrip];
0233       itemvector[istrip].xtalk_slope_left = new_slope_l[mystrip];
0234       itemvector[istrip].xtalk_intercept_left = new_intercept_l[mystrip];
0235       itemvector[istrip].xtalk_chi2_left = new_chi2_l[mystrip];
0236       cncrosstalk->crosstalk[new_cham_id[mystrip]] = itemvector;
0237       istrip++;
0238     }
0239
0240     itemvector.resize(64);
0241     for (unsigned int mystrip = 0; mystrip < new_nrlines - 1; mystrip++) {
0242       if (new_strips[mystrip] == 0)
0243         istrip = 0;
0244       if (new_cham_id[mystrip] >= 113000 && new_cham_id[mystrip] <= 113999) {
0245         itemvector[istrip].xtalk_slope_right = new_slope_r[mystrip];
0246         itemvector[istrip].xtalk_intercept_right = new_intercept_r[mystrip];
0247         itemvector[istrip].xtalk_chi2_right = new_chi2_r[mystrip];
0248         itemvector[istrip].xtalk_slope_left = new_slope_l[mystrip];
0249         itemvector[istrip].xtalk_intercept_left = new_intercept_l[mystrip];
0250         itemvector[istrip].xtalk_chi2_left = new_chi2_l[mystrip];
0251         cncrosstalk->crosstalk[new_cham_id[mystrip]] = itemvector;
0252         istrip++;
0253       }
0254     }
0255
0256     itemvector.resize(64);
0257     for (unsigned int mystrip = 0; mystrip < new_nrlines - 1; mystrip++) {
0258       if (new_strips[mystrip] == 0)
0259         istrip = 0;
0260       if (new_cham_id[mystrip] >= 213000 && new_cham_id[mystrip] <= 213999) {
0261         itemvector[istrip].xtalk_slope_right = new_slope_r[mystrip];
0262         itemvector[istrip].xtalk_intercept_right = new_intercept_r[mystrip];
0263         itemvector[istrip].xtalk_chi2_right = new_chi2_r[mystrip];
0264         itemvector[istrip].xtalk_slope_left = new_slope_l[mystrip];
0265         itemvector[istrip].xtalk_intercept_left = new_intercept_l[mystrip];
0266         itemvector[istrip].xtalk_chi2_left = new_chi2_l[mystrip];
0267         cncrosstalk->crosstalk[new_cham_id[mystrip]] = itemvector;
0268         istrip++;
0269       }
0270     }
0271   }
0272
0273   return cncrosstalk;
0274 }
0275
0276 CSCCrosstalkConditions::CSCCrosstalkConditions(const edm::ParameterSet &iConfig) {
0277   // the following line is needed to tell the framework what
0278   // data is being produced
0279   // added by Zhen (changed since 1_2_0)
0280   setWhatProduced(this, &CSCCrosstalkConditions::produceCrosstalk);
0281   findingRecord<CSCcrosstalkRcd>();
0282   // now do what ever other initialization is needed
0283 }
0284
0285 CSCCrosstalkConditions::~CSCCrosstalkConditions() {
0286   // do anything here that needs to be done at desctruction time
0287   // (e.g. close files, deallocate resources etc.)
0288 }
0289
0290 //
0291 // member functions
0292 //
0293
0294 // ------------ method called to produce the data  ------------
0295 CSCCrosstalkConditions::ReturnType CSCCrosstalkConditions::produceCrosstalk(const CSCcrosstalkRcd &iRecord) {
0296   // Added by Zhen, need a new object so to not be deleted at exit
0297   return CSCCrosstalkConditions::ReturnType(prefillCrosstalk());
0298 }
0299
0300 void CSCCrosstalkConditions::setIntervalFor(const edm::eventsetup::EventSetupRecordKey &,
0301                                             const edm::IOVSyncValue &,
0302                                             edm::ValidityInterval &oValidity) {
0303   oValidity = edm::ValidityInterval(edm::IOVSyncValue::beginOfTime(), edm::IOVSyncValue::endOfTime());
0304 }