Line Code
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234

/*
 *  See header file for a description of this class.
 *
 *  \author S. Maselli - INFN Torino
 */

#include "CalibMuon/DTCalibration/plugins/DTTTrigCorrectionFirst.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "Geometry/Records/interface/MuonGeometryRecord.h"
#include "Geometry/DTGeometry/interface/DTGeometry.h"
#include "Geometry/DTGeometry/interface/DTSuperLayer.h"

#include "CondFormats/DTObjects/interface/DTTtrig.h"
#include "CondFormats/DataRecord/interface/DTTtrigRcd.h"
#include "CalibMuon/DTCalibration/interface/DTCalibDBUtils.h"

#include <iostream>
#include <fstream>

using namespace edm;
using namespace std;

DTTTrigCorrectionFirst::DTTTrigCorrectionFirst(const ParameterSet& pset)
    : dtGeomToken_(esConsumes<edm::Transition::BeginRun>()),
      ttrigToken_(
          esConsumes<edm::Transition::BeginRun>(edm::ESInputTag("", pset.getUntrackedParameter<string>("dbLabel")))) {
  debug = pset.getUntrackedParameter<bool>("debug", false);
  ttrigMin = pset.getUntrackedParameter<double>("ttrigMin", 0);
  ttrigMax = pset.getUntrackedParameter<double>("ttrigMax", 5000);
  rmsLimit = pset.getUntrackedParameter<double>("rmsLimit", 5.);
}

DTTTrigCorrectionFirst::~DTTTrigCorrectionFirst() {}

void DTTTrigCorrectionFirst::beginRun(const edm::Run& run, const edm::EventSetup& setup) {
  ESHandle<DTTtrig> tTrig;
  tTrig = setup.getHandle(ttrigToken_);
  tTrigMap = &setup.getData(ttrigToken_);
  cout << "[DTTTrigCorrection]: TTrig version: " << tTrig->version() << endl;
  muonGeom = setup.getHandle(dtGeomToken_);
}

void DTTTrigCorrectionFirst::endJob() {
  // Create the object to be written to DB
  DTTtrig tTrigNewMap;
  //Get the superlayers list
  vector<const DTSuperLayer*> dtSupLylist = muonGeom->superLayers();

  //Loop on all superlayers to compute the mean
  double average = 0.;
  double average2 = 0.;
  double rms = 0.;
  double averageSigma = 0.;
  double counter = 0.;
  double averagekfactor = 0;
  float kfactor = 0;

  for (auto sl = dtSupLylist.begin(); sl != dtSupLylist.end(); sl++) {
    float ttrigMean = 0;
    float ttrigSigma = 0;
    tTrigMap->get((*sl)->id(), ttrigMean, ttrigSigma, kfactor, DTTimeUnits::ns);
    if (ttrigMean < ttrigMax && ttrigMean > ttrigMin) {
      average += ttrigMean;
      averageSigma += ttrigSigma;
      if (ttrigMean > 0)
        counter += 1.;
    }
  }  //End of loop on superlayers

  average = average / counter;
  averageSigma = averageSigma / counter;
  averagekfactor = kfactor;

  //  cout << " average counter "<< average << " "<< counter <<endl;

  for (auto sl = dtSupLylist.begin(); sl != dtSupLylist.end(); sl++) {
    float ttrigMean = 0;
    float ttrigSigma = 0;
    float kfactor = 0;
    tTrigMap->get((*sl)->id(), ttrigMean, ttrigSigma, kfactor, DTTimeUnits::ns);
    if (ttrigMean < ttrigMax && ttrigMean > ttrigMin) {
      average2 += (ttrigMean - average) * (ttrigMean - average);
    }
  }  //End of loop on superlayers

  rms = average2 / (counter - 1);
  rms = sqrt(rms);
  cout << "average averageSigma counter rms " << average << " " << averageSigma << " " << counter << " " << rms << endl;

  for (auto sl = dtSupLylist.begin(); sl != dtSupLylist.end(); sl++) {
    //Compute new ttrig
    double newTTrigMean = 0;
    double newTTrigSigma = 0;
    double newkfactor = 0;
    float tempttrigMean = 0;
    float tempttrigSigma = 0;
    float tempkfactor = 0;
    float ttrigMean = 0;
    float ttrigSigma = 0;
    float kfactor = 0;
    tTrigMap->get((*sl)->id(), ttrigMean, ttrigSigma, kfactor, DTTimeUnits::ns);
    int chamber = (*sl)->id().chamberId().station();

    //check if ttrigMean is similar to the mean
    if (abs(ttrigMean - average) < rmsLimit * rms) {
      newTTrigMean = ttrigMean;
      newTTrigSigma = ttrigSigma;
      newkfactor = kfactor;
    } else {
      // do not consider if ttrig == 0
      if (ttrigMean > 0) {
        //cout << "ttrig chamber " << ttrigMean <<" "<<chamber<<endl;
        if (((*sl)->id().superlayer()) == 1) {
          //cout << " superlayer " << ((*sl)->id().superlayer()) <<endl;
          DTSuperLayerId slId((*sl)->id().chamberId(), 3);
          tTrigMap->get(slId, tempttrigMean, tempttrigSigma, tempkfactor, DTTimeUnits::ns);
          if (abs(tempttrigMean - average) < rmsLimit * rms) {
            newTTrigMean = tempttrigMean;
            newTTrigSigma = tempttrigSigma;
            newkfactor = tempkfactor;
            cout << "Chamber " << chamber << " sl " << ((*sl)->id().superlayer()) << "has ttrig " << ttrigMean
                 << " -> takes value of sl 3 " << newTTrigMean << endl;
          } else if (chamber == 4) {
            cout << "No correction possible within same chamber (sl1) " << endl;
            newTTrigMean = average;
            newTTrigSigma = averageSigma;
            newkfactor = averagekfactor;
            cout << "####### Bad SL: " << (*sl)->id() << " from " << ttrigMean << " to " << newTTrigMean << endl;
          } else if (chamber != 4) {
            DTSuperLayerId slId((*sl)->id().chamberId(), 2);
            tTrigMap->get(slId, tempttrigMean, tempttrigSigma, tempkfactor, DTTimeUnits::ns);
            if (abs(tempttrigMean - average) < rmsLimit * rms) {
              newTTrigMean = tempttrigMean;
              newTTrigSigma = tempttrigSigma;
              newkfactor = tempkfactor;
              cout << "Chamber " << chamber << " sl " << ((*sl)->id().superlayer()) << "has ttrig " << ttrigMean
                   << " -> takes value of sl 2 " << newTTrigMean << endl;
            } else {
              cout << "No correction possible within same chamber (sl1) " << endl;
              newTTrigMean = average;
              newTTrigSigma = averageSigma;
              newkfactor = averagekfactor;
              cout << "####### Bad SL: " << (*sl)->id() << " from " << ttrigMean << " to " << newTTrigMean << endl;
            }
          }
        } else if (((*sl)->id().superlayer()) == 2) {
          //cout << " superlayer " << ((*sl)->id().superlayer()) <<endl;
          DTSuperLayerId slId((*sl)->id().chamberId(), 1);
          tTrigMap->get(slId, tempttrigMean, tempttrigSigma, tempkfactor, DTTimeUnits::ns);
          if (abs(tempttrigMean - average) < rmsLimit * rms) {
            newTTrigMean = tempttrigMean;
            newTTrigSigma = tempttrigSigma;
            newkfactor = tempkfactor;
            cout << "Chamber " << chamber << " sl " << ((*sl)->id().superlayer()) << "has ttrig " << ttrigMean
                 << " -> takes value of sl 1 " << newTTrigMean << endl;
          } else {
            DTSuperLayerId slId((*sl)->id().chamberId(), 3);
            tTrigMap->get(slId, tempttrigMean, tempttrigSigma, tempkfactor, DTTimeUnits::ns);
            if (abs(tempttrigMean - average) < rmsLimit * rms) {
              newTTrigMean = tempttrigMean;
              newTTrigSigma = tempttrigSigma;
              newkfactor = tempkfactor;
              cout << "Chamber " << chamber << " sl " << ((*sl)->id().superlayer()) << "has ttrig " << ttrigMean
                   << "-> takes value of sl 3 " << newTTrigMean << endl;
            } else {
              cout << "No correction possible within same chamber (sl2)  " << endl;
              newTTrigMean = average;
              newTTrigSigma = averageSigma;
              newkfactor = averagekfactor;
              cout << "####### Bad SL: " << (*sl)->id() << " from " << ttrigMean << " to " << newTTrigMean << endl;
            }
          }

        } else if (((*sl)->id().superlayer()) == 3) {
          //cout << " superlayer " << ((*sl)->id().superlayer()) <<endl;
          DTSuperLayerId slId((*sl)->id().chamberId(), 1);
          tTrigMap->get(slId, tempttrigMean, tempttrigSigma, tempkfactor, DTTimeUnits::ns);
          if (abs(tempttrigMean - average) < rmsLimit * rms) {
            newTTrigMean = tempttrigMean;
            newTTrigSigma = tempttrigSigma;
            newkfactor = tempkfactor;
            cout << "Chamber " << chamber << " sl " << ((*sl)->id().superlayer()) << "has ttrig " << ttrigMean
                 << " -> takes value of sl 1 " << newTTrigMean << endl;
          } else if (chamber == 4) {
            cout << "No correction possible within same chamber (sl3)" << endl;
            newTTrigMean = average;
            newTTrigSigma = averageSigma;
            newkfactor = averagekfactor;
            cout << "####### Bad SL: " << (*sl)->id() << " from " << ttrigMean << " to " << newTTrigMean << endl;
          } else if (chamber != 4) {
            DTSuperLayerId slId((*sl)->id().chamberId(), 2);
            tTrigMap->get(slId, tempttrigMean, tempttrigSigma, tempkfactor, DTTimeUnits::ns);
            if (abs(tempttrigMean - average) < rmsLimit * rms) {
              newTTrigMean = tempttrigMean;
              newTTrigSigma = tempttrigSigma;
              newkfactor = tempkfactor;
              cout << "Chamber " << chamber << " sl " << ((*sl)->id().superlayer()) << "has ttrig " << ttrigMean
                   << " -> takes value of sl 2 " << newTTrigMean << endl;
            } else {
              cout << "No correction possible within same chamber (sl3) " << endl;
              newTTrigMean = average;
              newTTrigSigma = averageSigma;
              newkfactor = averagekfactor;
              cout << "####### Bad SL: " << (*sl)->id() << " from " << ttrigMean << " to " << newTTrigMean << endl;
            }
          }
        }
      } else {
        //	cout << "SL not present " << ttrigMean << " " <<((*sl)->id().superlayer()) <<endl;
        newTTrigMean = average;
        newTTrigSigma = averageSigma;
        newkfactor = averagekfactor;
        cout << "####### NotPresent SL: " << (*sl)->id() << " from " << ttrigMean << " to " << newTTrigMean << endl;
      }
    }

    //Store new ttrig in the new map
    tTrigNewMap.set((*sl)->id(), newTTrigMean, newTTrigSigma, newkfactor, DTTimeUnits::ns);
    if (debug) {
      cout << "New tTrig: " << (*sl)->id() << " from " << ttrigMean << " to " << newTTrigMean << endl;
      cout << "New tTrigSigma: " << (*sl)->id() << " from " << ttrigSigma << " to " << newTTrigSigma << endl;
    }
  }  //End of loop on superlayers

  //Write object to DB
  cout << "[DTTTrigCorrection]: Writing ttrig object to DB!" << endl;
  string record = "DTTtrigRcd";
  DTCalibDBUtils::writeToDB<DTTtrig>(record, tTrigNewMap);
}