Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:58:28 

0001 
0002 /*
0003  *  See header file for a description of this class.
0004  *
0005  *  \author M. Giunta
0006  */
0007 
0008 #include "CalibMuon/DTCalibration/plugins/DTVDriftCalibration.h"
0009 #include "CalibMuon/DTCalibration/interface/DTMeanTimerFitter.h"
0010 #include "CalibMuon/DTCalibration/interface/DTCalibDBUtils.h"
0011 
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/Framework/interface/ESHandle.h"
0014 #include "FWCore/Framework/interface/ConsumesCollector.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016 
0017 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0018 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0019 
0020 #include "CalibMuon/DTDigiSync/interface/DTTTrigSyncFactory.h"
0021 #include "CalibMuon/DTDigiSync/interface/DTTTrigBaseSync.h"
0022 
0023 #include "CondFormats/DTObjects/interface/DTMtime.h"
0024 #include "CondFormats/DTObjects/interface/DTRecoConditions.h"
0025 
0026 #include "CondFormats/DataRecord/interface/DTStatusFlagRcd.h"
0027 #include "CondFormats/DTObjects/interface/DTStatusFlag.h"
0028 
0029 /* C++ Headers */
0030 #include <map>
0031 #include <iostream>
0032 #include <fstream>
0033 #include <string>
0034 #include <sstream>
0035 #include "TFile.h"
0036 #include "TH1.h"
0037 #include "TF1.h"
0038 #include "TROOT.h"
0039 
0040 using namespace std;
0041 using namespace edm;
0042 using namespace dttmaxenums;
0043 
0044 DTVDriftCalibration::DTVDriftCalibration(const ParameterSet& pset)
0045     :  // Get the synchronizer
0046       theDTGeomToken{esConsumes()},
0047       theSync{DTTTrigSyncFactory::get()->create(pset.getParameter<string>("tTrigMode"),
0048                                                 pset.getParameter<ParameterSet>("tTrigModeConfig"),
0049                                                 consumesCollector())}
0050 
0051 {
0052   edm::ConsumesCollector collector(consumesCollector());
0053   select_ = std::make_unique<DTSegmentSelector>(pset, collector);
0054 
0055   // The name of the 4D rec hits collection
0056   theRecHits4DToken = (consumes<DTRecSegment4DCollection>(pset.getParameter<InputTag>("recHits4DLabel")));
0057 
0058   // The root file which will contain the histos
0059   string rootFileName = pset.getUntrackedParameter<string>("rootFileName");
0060   theFile = new TFile(rootFileName.c_str(), "RECREATE");
0061   theFile->cd();
0062 
0063   debug = pset.getUntrackedParameter<bool>("debug", false);
0064 
0065   theFitter = std::make_unique<DTMeanTimerFitter>(theFile);
0066   if (debug)
0067     theFitter->setVerbosity(1);
0068 
0069   hChi2 = new TH1F("hChi2", "Chi squared tracks", 100, 0, 100);
0070   h2DSegmRPhi = new h2DSegm("SLRPhi");
0071   h2DSegmRZ = new h2DSegm("SLRZ");
0072   h4DSegmAllCh = new h4DSegm("AllCh");
0073   histograms_.bookHistos();
0074 
0075   findVDriftAndT0 = pset.getUntrackedParameter<bool>("fitAndWrite", false);
0076 
0077   // Chamber/s to calibrate
0078   theCalibChamber = pset.getUntrackedParameter<string>("calibChamber", "All");
0079 
0080   // the txt file which will contain the calibrated constants
0081   theVDriftOutputFile = pset.getUntrackedParameter<string>("vDriftFileName");
0082 
0083   // get parameter set for DTCalibrationMap constructor
0084   theCalibFilePar = pset.getUntrackedParameter<ParameterSet>("calibFileConfig");
0085 
0086   // the granularity to be used for tMax
0087   string tMaxGranularity = pset.getUntrackedParameter<string>("tMaxGranularity", "bySL");
0088 
0089   writeLegacyVDriftDB = pset.getParameter<bool>("writeLegacyVDriftDB");
0090 
0091   // Enforce it to be by SL since rest is not implemented
0092   if (tMaxGranularity != "bySL") {
0093     LogError("Calibration") << "[DTVDriftCalibration] tMaxGranularity will be fixed to bySL.";
0094     tMaxGranularity = "bySL";
0095   }
0096   // Initialize correctly the enum which specify the granularity for the calibration
0097   if (tMaxGranularity == "bySL") {
0098     theGranularity = bySL;
0099   } else if (tMaxGranularity == "byChamber") {
0100     theGranularity = byChamber;
0101   } else if (tMaxGranularity == "byPartition") {
0102     theGranularity = byPartition;
0103   } else
0104     throw cms::Exception("Configuration")
0105         << "[DTVDriftCalibration] Check parameter tMaxGranularity: " << tMaxGranularity << " option not available";
0106 
0107   LogVerbatim("Calibration") << "[DTVDriftCalibration]Constructor called!";
0108 }
0109 
0110 DTVDriftCalibration::~DTVDriftCalibration() {
0111   theFile->Close();
0112   LogVerbatim("Calibration") << "[DTVDriftCalibration]Destructor called!";
0113 }
0114 
0115 void DTVDriftCalibration::analyze(const Event& event, const EventSetup& eventSetup) {
0116   LogTrace("Calibration") << "--- [DTVDriftCalibration] Event analysed #Run: " << event.id().run()
0117                           << " #Event: " << event.id().event();
0118   theFile->cd();
0119   DTChamberId chosenChamberId;
0120 
0121   if (theCalibChamber != "All") {
0122     stringstream linestr;
0123     int selWheel, selStation, selSector;
0124     linestr << theCalibChamber;
0125     linestr >> selWheel >> selStation >> selSector;
0126     chosenChamberId = DTChamberId(selWheel, selStation, selSector);
0127     LogTrace("Calibration") << "chosen chamber " << chosenChamberId;
0128   }
0129 
0130   // Get the DT Geometry
0131   const DTGeometry& dtGeom = eventSetup.getData(theDTGeomToken);
0132 
0133   // Get the rechit collection from the event
0134   Handle<DTRecSegment4DCollection> all4DSegments;
0135   event.getByToken(theRecHits4DToken, all4DSegments);
0136 
0137   // Set the event setup in the Synchronizer
0138   theSync->setES(eventSetup);
0139 
0140   // Loop over segments by chamber
0141   DTRecSegment4DCollection::id_iterator chamberIdIt;
0142   for (chamberIdIt = all4DSegments->id_begin(); chamberIdIt != all4DSegments->id_end(); ++chamberIdIt) {
0143     // Get the chamber from the setup
0144     const DTChamber* chamber = dtGeom.chamber(*chamberIdIt);
0145     LogTrace("Calibration") << "Chamber Id: " << *chamberIdIt;
0146 
0147     // Calibrate just the chosen chamber/s
0148     if ((theCalibChamber != "All") && ((*chamberIdIt) != chosenChamberId))
0149       continue;
0150 
0151     // Get the range for the corresponding ChamberId
0152     DTRecSegment4DCollection::range range = all4DSegments->get((*chamberIdIt));
0153 
0154     // Loop over the rechits of this DetUnit
0155     for (DTRecSegment4DCollection::const_iterator segment = range.first; segment != range.second; ++segment) {
0156       if (!(*segment).hasZed() && !(*segment).hasPhi()) {
0157         LogError("Calibration") << "4D segment without Z and Phi segments";
0158         continue;
0159       }
0160 
0161       LogTrace("Calibration") << "Segment local pos (in chamber RF): " << (*segment).localPosition()
0162                               << "\nSegment global pos: " << chamber->toGlobal((*segment).localPosition());
0163 
0164       if (!((*select_)(*segment, event, eventSetup)))
0165         continue;
0166 
0167       LocalPoint phiSeg2DPosInCham;
0168       LocalVector phiSeg2DDirInCham;
0169       map<DTSuperLayerId, vector<DTRecHit1D> > hitsBySLMap;
0170       if ((*segment).hasPhi()) {
0171         const DTChamberRecSegment2D* phiSeg = (*segment).phiSegment();  // phiSeg lives in the chamber RF
0172         phiSeg2DPosInCham = phiSeg->localPosition();
0173         phiSeg2DDirInCham = phiSeg->localDirection();
0174         vector<DTRecHit1D> phiHits = phiSeg->specificRecHits();
0175         for (vector<DTRecHit1D>::const_iterator hit = phiHits.begin(); hit != phiHits.end(); ++hit) {
0176           DTWireId wireId = (*hit).wireId();
0177           DTSuperLayerId slId = wireId.superlayerId();
0178           hitsBySLMap[slId].push_back(*hit);
0179         }
0180       }
0181       // Get the Theta 2D segment and plot the angle in the chamber RF
0182       LocalVector zSeg2DDirInCham;
0183       LocalPoint zSeg2DPosInCham;
0184       if ((*segment).hasZed()) {
0185         const DTSLRecSegment2D* zSeg = (*segment).zSegment();  // zSeg lives in the SL RF
0186         const DTSuperLayer* sl = chamber->superLayer(zSeg->superLayerId());
0187         zSeg2DPosInCham = chamber->toLocal(sl->toGlobal((*zSeg).localPosition()));
0188         zSeg2DDirInCham = chamber->toLocal(sl->toGlobal((*zSeg).localDirection()));
0189         hitsBySLMap[zSeg->superLayerId()] = zSeg->specificRecHits();
0190       }
0191 
0192       LocalPoint segment4DLocalPos = (*segment).localPosition();
0193       LocalVector segment4DLocalDir = (*segment).localDirection();
0194       double chiSquare = ((*segment).chi2() / (*segment).degreesOfFreedom());
0195 
0196       hChi2->Fill(chiSquare);
0197       if ((*segment).hasPhi())
0198         h2DSegmRPhi->Fill(phiSeg2DPosInCham.x(), phiSeg2DDirInCham.x() / phiSeg2DDirInCham.z());
0199       if ((*segment).hasZed())
0200         h2DSegmRZ->Fill(zSeg2DPosInCham.y(), zSeg2DDirInCham.y() / zSeg2DDirInCham.z());
0201 
0202       if ((*segment).hasZed() && (*segment).hasPhi())
0203         h4DSegmAllCh->Fill(segment4DLocalPos.x(),
0204                            segment4DLocalPos.y(),
0205                            atan(segment4DLocalDir.x() / segment4DLocalDir.z()) * 180. / Geom::pi(),
0206                            atan(segment4DLocalDir.y() / segment4DLocalDir.z()) * 180. / Geom::pi(),
0207                            180 - segment4DLocalDir.theta() * 180. / Geom::pi());
0208       else if ((*segment).hasPhi())
0209         h4DSegmAllCh->Fill(segment4DLocalPos.x(),
0210                            atan(segment4DLocalDir.x() / segment4DLocalDir.z()) * 180. / Geom::pi());
0211       else if ((*segment).hasZed())
0212         LogWarning("Calibration") << "4d segment with only Z";
0213 
0214       //loop over the segments
0215       for (map<DTSuperLayerId, vector<DTRecHit1D> >::const_iterator slIdAndHits = hitsBySLMap.begin();
0216            slIdAndHits != hitsBySLMap.end();
0217            ++slIdAndHits) {
0218         if (slIdAndHits->second.size() < 3)
0219           continue;
0220         DTSuperLayerId slId = slIdAndHits->first;
0221 
0222         // Create the DTTMax, that computes the 4 TMax
0223         DTTMax slSeg(slIdAndHits->second,
0224                      *(chamber->superLayer(slIdAndHits->first)),
0225                      chamber->toGlobal((*segment).localDirection()),
0226                      chamber->toGlobal((*segment).localPosition()),
0227                      *theSync,
0228                      histograms_);
0229 
0230         if (theGranularity == bySL) {
0231           vector<const TMax*> tMaxes = slSeg.getTMax(slId);
0232           DTWireId wireId(slId, 0, 0);
0233           theFile->cd();
0234           cellInfo* cell = theWireIdAndCellMap[wireId];
0235           if (cell == nullptr) {
0236             TString name = (((((TString) "TMax" + (long)slId.wheel()) + (long)slId.station()) + (long)slId.sector()) +
0237                             (long)slId.superLayer());
0238             cell = new cellInfo(name);
0239             theWireIdAndCellMap[wireId] = cell;
0240           }
0241           cell->add(tMaxes);
0242           cell->update();  // FIXME to reset the counter to avoid triple counting, which actually is not used...
0243         } else {
0244           LogWarning("Calibration") << "[DTVDriftCalibration] ###Warning: the chosen granularity is not implemented "
0245                                        "yet, only bySL available!";
0246         }
0247         // to be implemented: granularity different from bySL
0248 
0249         //       else if (theGranularity == byPartition) {
0250         //  // Use the custom granularity defined by partition();
0251         //  // in this case, add() should be called once for each Tmax of each layer
0252         //  // and triple counting should be avoided within add()
0253         //  vector<cellInfo*> cells;
0254         //  for (int i=1; i<=4; i++) {
0255         //    const DTTMax::InfoLayer* iLayer = slSeg.getInfoLayer(i);
0256         //    if(iLayer == 0) continue;
0257         //    cellInfo * cell = partition(iLayer->idWire);
0258         //    cells.push_back(cell);
0259         //    vector<const TMax*> tMaxes = slSeg.getTMax(iLayer->idWire);
0260         //    cell->add(tMaxes);
0261         //  }
0262         //  //reset the counter to avoid triple counting
0263         //  for (vector<cellInfo*>::const_iterator i = cells.begin();
0264         //       i!= cells.end(); i++) {
0265         //    (*i)->update();
0266         //  }
0267         //       }
0268       }
0269     }
0270   }
0271 }
0272 
0273 void DTVDriftCalibration::endJob() {
0274   theFile->cd();
0275   gROOT->GetList()->Write();
0276   h2DSegmRPhi->Write();
0277   h2DSegmRZ->Write();
0278   h4DSegmAllCh->Write();
0279   hChi2->Write();
0280   // Instantiate a DTCalibrationMap object if you want to calculate the calibration constants
0281   DTCalibrationMap calibValuesFile(theCalibFilePar);
0282   // Create the object to be written to DB
0283   std::unique_ptr<DTMtime> mTime;
0284   std::unique_ptr<DTRecoConditions> vDrift;
0285   if (writeLegacyVDriftDB) {
0286     mTime = std::make_unique<DTMtime>();
0287   } else {
0288     vDrift = std::make_unique<DTRecoConditions>();
0289     vDrift->setFormulaExpr("[0]");
0290     //vDriftNewMap->setFormulaExpr("[0]*(1-[1]*x)"); // add parametrization for dependency along Y
0291     vDrift->setVersion(1);
0292   }
0293 
0294   // write the TMax histograms of each SL to the root file
0295   if (theGranularity == bySL) {
0296     for (map<DTWireId, cellInfo*>::const_iterator wireCell = theWireIdAndCellMap.begin();
0297          wireCell != theWireIdAndCellMap.end();
0298          ++wireCell) {
0299       cellInfo* cell = theWireIdAndCellMap[(*wireCell).first];
0300       hTMaxCell* cellHists = cell->getHists();
0301       theFile->cd();
0302       cellHists->Write();
0303       if (findVDriftAndT0) {  // if TRUE: evaluate calibration constants from TMax hists filled in this job
0304         // evaluate v_drift and sigma from the TMax histograms
0305         DTWireId wireId = (*wireCell).first;
0306         vector<float> newConstants;
0307         TString N = (((((TString) "TMax" + (long)wireId.wheel()) + (long)wireId.station()) + (long)wireId.sector()) +
0308                      (long)wireId.superLayer());
0309         vector<float> vDriftAndReso = theFitter->evaluateVDriftAndReso(N);
0310 
0311         // Don't write the constants for the SL if the vdrift was not computed
0312         if (vDriftAndReso.front() == -1)
0313           continue;
0314         const DTCalibrationMap::CalibConsts* oldConstants = calibValuesFile.getConsts(wireId);
0315         if (oldConstants != nullptr) {
0316           newConstants.push_back((*oldConstants)[0]);
0317           newConstants.push_back((*oldConstants)[1]);
0318           newConstants.push_back((*oldConstants)[2]);
0319         } else {
0320           newConstants.push_back(-1);
0321           newConstants.push_back(-1);
0322           newConstants.push_back(-1);
0323         }
0324         for (int ivd = 0; ivd <= 5; ivd++) {
0325           // 0=vdrift, 1=reso, 2=(3deltat0-2deltat0), 3=(2deltat0-1deltat0),
0326           //  4=(1deltat0-0deltat0), 5=deltat0 from hists with max entries,
0327           newConstants.push_back(vDriftAndReso[ivd]);
0328         }
0329 
0330         calibValuesFile.addCell(calibValuesFile.getKey(wireId), newConstants);
0331 
0332         // vdrift is cm/ns , resolution is cm
0333         if (writeLegacyVDriftDB) {
0334           mTime->set((wireId.layerId()).superlayerId(), vDriftAndReso[0], vDriftAndReso[1], DTVelocityUnits::cm_per_ns);
0335         } else {
0336           vector<double> params = {vDriftAndReso[0]};
0337           vDrift->set(wireId, params);
0338         }
0339         LogTrace("Calibration") << " SL: " << (wireId.layerId()).superlayerId() << " vDrift = " << vDriftAndReso[0]
0340                                 << " reso = " << vDriftAndReso[1];
0341       }
0342     }
0343   }
0344 
0345   // to be implemented: granularity different from bySL
0346 
0347   //   if(theGranularity == "byChamber") {
0348   //     const vector<DTChamber*> chambers = dMap.chambers();
0349 
0350   //     // Loop over all chambers
0351   //     for(vector<MuBarChamber*>::const_iterator chamber = chambers.begin();
0352   //    chamber != chambers.end(); chamber ++) {
0353   //       MuBarChamberId chamber_id = (*chamber)->id();
0354   //       MuBarDigiParameters::Key wire_id(chamber_id, 0, 0, 0);
0355   //       vector<float> newConstants;
0356   //       vector<float> vDriftAndReso = evaluateVDriftAndReso(wire_id, f);
0357   //       const CalibConsts* oldConstants = digiParams.getConsts(wire_id);
0358   //       if(oldConstants !=0) {
0359   //    newConstants = *oldConstants;
0360   //    newConstants.push_back(vDriftAndReso[0]);
0361   //    newConstants.push_back(vDriftAndReso[1]);
0362   //    newConstants.push_back(vDriftAndReso[2]);
0363   //    newConstants.push_back(vDriftAndReso[3]);
0364   //       } else {
0365   //    newConstants.push_back(-1);
0366   //    newConstants.push_back(-1);
0367   //    newConstants.push_back(vDriftAndReso[0]);
0368   //    newConstants.push_back(vDriftAndReso[1]);
0369   //    newConstants.push_back(vDriftAndReso[2]);
0370   //    newConstants.push_back(vDriftAndReso[3]);
0371   //       }
0372   //       digiParams.addCell(wire_id, newConstants);
0373   //     }
0374   //   }
0375 
0376   // Write values to a table
0377   calibValuesFile.writeConsts(theVDriftOutputFile);
0378 
0379   LogVerbatim("Calibration") << "[DTVDriftCalibration]Writing vdrift object to DB!";
0380 
0381   // Write the vdrift object to DB
0382   if (writeLegacyVDriftDB) {
0383     string record = "DTMtimeRcd";
0384     DTCalibDBUtils::writeToDB<DTMtime>(record, *mTime);
0385   } else {
0386     DTCalibDBUtils::writeToDB<DTRecoConditions>("DTRecoConditionsVdriftRcd", *vDrift);
0387   }
0388 }
0389 
0390 // to be implemented: granularity different from bySL
0391 
0392 // // Create partitions
0393 // DTVDriftCalibration::cellInfo* DTVDriftCalibration::partition(const DTWireId& wireId) {
0394 //   for( map<MuBarWireId, cellInfo*>::const_iterator iter =
0395 //   mapCellTmaxPart.begin(); iter != mapCellTmaxPart.end(); iter++) {
0396 //     // Divide wires per SL (with phi symmetry)
0397 //     if(iter->first.wheel() == wireId.wheel() &&
0398 //        iter->first.station() == wireId.station() &&
0399 //        //       iter->first.sector() == wireId.sector() && // phi symmetry!
0400 //        iter->first.superlayer() == wireId.superlayer()) {
0401 //       return iter->second;
0402 //     }
0403 //   }
0404 //   cellInfo * result = new cellInfo("dummy string"); // FIXME: change constructor; create tree?
0405 //   mapCellTmaxPart.insert(make_pair(wireId, result));
0406 //   return result;
0407 //}
0408 
0409 void DTVDriftCalibration::cellInfo::add(const vector<const TMax*>& _tMaxes) {
0410   vector<const TMax*> tMaxes = _tMaxes;
0411   float tmax123 = -1.;
0412   float tmax124 = -1.;
0413   float tmax134 = -1.;
0414   float tmax234 = -1.;
0415   SigmaFactor s124 = noR;
0416   SigmaFactor s134 = noR;
0417   unsigned t0_123 = 0;
0418   unsigned t0_124 = 0;
0419   unsigned t0_134 = 0;
0420   unsigned t0_234 = 0;
0421   unsigned hSubGroup = 0;
0422   for (vector<const TMax*>::const_iterator it = tMaxes.begin(); it != tMaxes.end(); ++it) {
0423     if (*it == nullptr) {
0424       continue;
0425     } else {
0426       //FIXME check cached,
0427       if (addedCells.size() == 4 || find(addedCells.begin(), addedCells.end(), (*it)->cells) != addedCells.end()) {
0428         continue;
0429       }
0430       addedCells.push_back((*it)->cells);
0431       SigmaFactor sigma = (*it)->sigma;
0432       float t = (*it)->t;
0433       TMaxCells cells = (*it)->cells;
0434       unsigned t0Factor = (*it)->t0Factor;
0435       hSubGroup = (*it)->hSubGroup;
0436       if (t < 0.)
0437         continue;
0438       switch (cells) {
0439         case notInit:
0440           cout << "Error: no cell type assigned to TMax" << endl;
0441           break;
0442         case c123:
0443           tmax123 = t;
0444           t0_123 = t0Factor;
0445           break;
0446         case c124:
0447           tmax124 = t;
0448           s124 = sigma;
0449           t0_124 = t0Factor;
0450           break;
0451         case c134:
0452           tmax134 = t;
0453           s134 = sigma;
0454           t0_134 = t0Factor;
0455           break;
0456         case c234:
0457           tmax234 = t;
0458           t0_234 = t0Factor;
0459           break;
0460       }
0461     }
0462   }
0463   //add entries to the TMax histograms
0464   histos->Fill(tmax123, tmax124, tmax134, tmax234, s124, s134, t0_123, t0_124, t0_134, t0_234, hSubGroup);
0465 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!