File indexing completed on 2024-04-06 11:55:57
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022 #include <fstream>
0023 #include <sstream>
0024
0025
0026 #include "FWCore/Framework/interface/Frameworkfwd.h"
0027 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0028 #include "FWCore/Framework/interface/EventSetup.h"
0029 #include "FWCore/Framework/interface/Event.h"
0030 #include "FWCore/Framework/interface/MakerMacros.h"
0031
0032 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0033 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0034 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0035
0036 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0037 #include "FWCore/Utilities/interface/EDMException.h"
0038
0039 #include "CLHEP/Matrix/SymMatrix.h"
0040
0041 #include "CondFormats/AlignmentRecord/interface/TrackerAlignmentErrorExtendedRcd.h"
0042 #include "CondFormats/Alignment/interface/AlignmentErrorsExtended.h"
0043
0044
0045 #include "Alignment/APEEstimation/interface/TrackerSectorStruct.h"
0046 #include "Alignment/APEEstimation/interface/ReducedTrackerTreeVariables.h"
0047
0048 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0049 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0050 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
0051 #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
0052 #include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
0053
0054 #include "TString.h"
0055 #include "TFile.h"
0056 #include "TDirectory.h"
0057 #include "TTree.h"
0058 #include "TMath.h"
0059
0060
0061
0062
0063 class ApeTreeCreateDefault : public edm::one::EDAnalyzer<> {
0064 public:
0065 explicit ApeTreeCreateDefault(const edm::ParameterSet&);
0066 ~ApeTreeCreateDefault() override;
0067
0068 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0069
0070 private:
0071 void beginJob() override;
0072 void analyze(const edm::Event&, const edm::EventSetup&) override;
0073 void endJob() override;
0074
0075 void sectorBuilder();
0076 bool checkIntervalsForSectors(const unsigned int sectorCounter, const std::vector<double>&) const;
0077 bool checkModuleIds(const unsigned int, const std::vector<unsigned int>&) const;
0078 bool checkModuleBools(const bool, const std::vector<unsigned int>&) const;
0079 bool checkModuleDirections(const int, const std::vector<int>&) const;
0080 bool checkModulePositions(const float, const std::vector<double>&) const;
0081
0082
0083 const edm::ESGetToken<AlignmentErrorsExtended, TrackerAlignmentErrorExtendedRcd> alignmentErrorToken_;
0084 const std::string resultFile_;
0085 const std::string trackerTreeFile_;
0086 const std::vector<edm::ParameterSet> sectors_;
0087
0088 std::map<unsigned int, TrackerSectorStruct> m_tkSector_;
0089 std::map<unsigned int, ReducedTrackerTreeVariables> m_tkTreeVar_;
0090 unsigned int noSectors;
0091 };
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104 ApeTreeCreateDefault::ApeTreeCreateDefault(const edm::ParameterSet& iConfig)
0105 : alignmentErrorToken_(esConsumes()),
0106 resultFile_(iConfig.getParameter<std::string>("resultFile")),
0107 trackerTreeFile_(iConfig.getParameter<std::string>("trackerTreeFile")),
0108 sectors_(iConfig.getParameter<std::vector<edm::ParameterSet>>("sectors")) {}
0109
0110 ApeTreeCreateDefault::~ApeTreeCreateDefault() {}
0111
0112
0113
0114
0115 void ApeTreeCreateDefault::sectorBuilder() {
0116
0117 TFile* tkTreeFile(TFile::Open(trackerTreeFile_.c_str()));
0118 if (tkTreeFile) {
0119 edm::LogInfo("SectorBuilder") << "TrackerTreeFile OK";
0120 } else {
0121 edm::LogError("SectorBuilder") << "TrackerTreeFile not found";
0122 return;
0123 }
0124 TTree* tkTree(nullptr);
0125 tkTreeFile->GetObject("TrackerTreeGenerator/TrackerTree/TrackerTree", tkTree);
0126 if (tkTree) {
0127 edm::LogInfo("SectorBuilder") << "TrackerTree OK";
0128 } else {
0129 edm::LogError("SectorBuilder") << "TrackerTree not found in file";
0130 return;
0131 }
0132
0133 unsigned int rawId(999), subdetId(999), layer(999), side(999), half(999), rod(999), ring(999), petal(999), blade(999),
0134 panel(999), outerInner(999), module(999), nStrips(999);
0135 bool isDoubleSide(false), isRPhi(false), isStereo(false);
0136 int uDirection(999), vDirection(999), wDirection(999);
0137 float posR(999.F), posPhi(999.F), posEta(999.F), posX(999.F), posY(999.F), posZ(999.F);
0138
0139 tkTree->SetBranchAddress("RawId", &rawId);
0140 tkTree->SetBranchAddress("SubdetId", &subdetId);
0141 tkTree->SetBranchAddress("Layer", &layer);
0142 tkTree->SetBranchAddress("Side", &side);
0143 tkTree->SetBranchAddress("Half", &half);
0144 tkTree->SetBranchAddress("Rod", &rod);
0145 tkTree->SetBranchAddress("Ring", &ring);
0146 tkTree->SetBranchAddress("Petal", &petal);
0147 tkTree->SetBranchAddress("Blade", &blade);
0148 tkTree->SetBranchAddress("Panel", &panel);
0149 tkTree->SetBranchAddress("OuterInner", &outerInner);
0150 tkTree->SetBranchAddress("Module", &module);
0151 tkTree->SetBranchAddress("NStrips", &nStrips);
0152 tkTree->SetBranchAddress("IsDoubleSide", &isDoubleSide);
0153 tkTree->SetBranchAddress("IsRPhi", &isRPhi);
0154 tkTree->SetBranchAddress("IsStereo", &isStereo);
0155 tkTree->SetBranchAddress("UDirection", &uDirection);
0156 tkTree->SetBranchAddress("VDirection", &vDirection);
0157 tkTree->SetBranchAddress("WDirection", &wDirection);
0158 tkTree->SetBranchAddress("PosR", &posR);
0159 tkTree->SetBranchAddress("PosPhi", &posPhi);
0160 tkTree->SetBranchAddress("PosEta", &posEta);
0161 tkTree->SetBranchAddress("PosX", &posX);
0162 tkTree->SetBranchAddress("PosY", &posY);
0163 tkTree->SetBranchAddress("PosZ", &posZ);
0164
0165 int nModules(tkTree->GetEntries());
0166 TrackerSectorStruct allSectors;
0167
0168
0169 unsigned int sectorCounter(0);
0170 std::vector<edm::ParameterSet> v_sectorDef(sectors_);
0171 edm::LogInfo("SectorBuilder") << "There are " << v_sectorDef.size() << " Sectors defined";
0172
0173 for (auto const& parSet : v_sectorDef) {
0174 ++sectorCounter;
0175 const std::string& sectorName(parSet.getParameter<std::string>("name"));
0176 std::vector<unsigned int> v_rawId(parSet.getParameter<std::vector<unsigned int>>("rawId")),
0177 v_subdetId(parSet.getParameter<std::vector<unsigned int>>("subdetId")),
0178 v_layer(parSet.getParameter<std::vector<unsigned int>>("layer")),
0179 v_side(parSet.getParameter<std::vector<unsigned int>>("side")),
0180 v_half(parSet.getParameter<std::vector<unsigned int>>("half")),
0181 v_rod(parSet.getParameter<std::vector<unsigned int>>("rod")),
0182 v_ring(parSet.getParameter<std::vector<unsigned int>>("ring")),
0183 v_petal(parSet.getParameter<std::vector<unsigned int>>("petal")),
0184 v_blade(parSet.getParameter<std::vector<unsigned int>>("blade")),
0185 v_panel(parSet.getParameter<std::vector<unsigned int>>("panel")),
0186 v_outerInner(parSet.getParameter<std::vector<unsigned int>>("outerInner")),
0187 v_module(parSet.getParameter<std::vector<unsigned int>>("module")),
0188 v_nStrips(parSet.getParameter<std::vector<unsigned int>>("nStrips")),
0189 v_isDoubleSide(parSet.getParameter<std::vector<unsigned int>>("isDoubleSide")),
0190 v_isRPhi(parSet.getParameter<std::vector<unsigned int>>("isRPhi")),
0191 v_isStereo(parSet.getParameter<std::vector<unsigned int>>("isStereo"));
0192 std::vector<int> v_uDirection(parSet.getParameter<std::vector<int>>("uDirection")),
0193 v_vDirection(parSet.getParameter<std::vector<int>>("vDirection")),
0194 v_wDirection(parSet.getParameter<std::vector<int>>("wDirection"));
0195 std::vector<double> v_posR(parSet.getParameter<std::vector<double>>("posR")),
0196 v_posPhi(parSet.getParameter<std::vector<double>>("posPhi")),
0197 v_posEta(parSet.getParameter<std::vector<double>>("posEta")),
0198 v_posX(parSet.getParameter<std::vector<double>>("posX")),
0199 v_posY(parSet.getParameter<std::vector<double>>("posY")),
0200 v_posZ(parSet.getParameter<std::vector<double>>("posZ"));
0201
0202 if (!this->checkIntervalsForSectors(sectorCounter, v_posR) ||
0203 !this->checkIntervalsForSectors(sectorCounter, v_posPhi) ||
0204 !this->checkIntervalsForSectors(sectorCounter, v_posEta) ||
0205 !this->checkIntervalsForSectors(sectorCounter, v_posX) ||
0206 !this->checkIntervalsForSectors(sectorCounter, v_posY) ||
0207 !this->checkIntervalsForSectors(sectorCounter, v_posZ)) {
0208 continue;
0209 }
0210
0211 TrackerSectorStruct tkSector;
0212 tkSector.name = sectorName;
0213
0214 ReducedTrackerTreeVariables tkTreeVar;
0215
0216
0217 for (int module = 0; module < nModules; ++module) {
0218 tkTree->GetEntry(module);
0219
0220 if (sectorCounter == 1) {
0221 tkTreeVar.subdetId = subdetId;
0222 tkTreeVar.nStrips = nStrips;
0223 tkTreeVar.uDirection = uDirection;
0224 tkTreeVar.vDirection = vDirection;
0225 tkTreeVar.wDirection = wDirection;
0226 m_tkTreeVar_[rawId] = tkTreeVar;
0227 }
0228
0229 if (!this->checkModuleIds(rawId, v_rawId))
0230 continue;
0231 if (!this->checkModuleIds(subdetId, v_subdetId))
0232 continue;
0233 if (!this->checkModuleIds(layer, v_layer))
0234 continue;
0235 if (!this->checkModuleIds(side, v_side))
0236 continue;
0237 if (!this->checkModuleIds(half, v_half))
0238 continue;
0239 if (!this->checkModuleIds(rod, v_rod))
0240 continue;
0241 if (!this->checkModuleIds(ring, v_ring))
0242 continue;
0243 if (!this->checkModuleIds(petal, v_petal))
0244 continue;
0245 if (!this->checkModuleIds(blade, v_blade))
0246 continue;
0247 if (!this->checkModuleIds(panel, v_panel))
0248 continue;
0249 if (!this->checkModuleIds(outerInner, v_outerInner))
0250 continue;
0251 if (!this->checkModuleIds(module, v_module))
0252 continue;
0253 if (!this->checkModuleIds(nStrips, v_nStrips))
0254 continue;
0255 if (!this->checkModuleBools(isDoubleSide, v_isDoubleSide))
0256 continue;
0257 if (!this->checkModuleBools(isRPhi, v_isRPhi))
0258 continue;
0259 if (!this->checkModuleBools(isStereo, v_isStereo))
0260 continue;
0261 if (!this->checkModuleDirections(uDirection, v_uDirection))
0262 continue;
0263 if (!this->checkModuleDirections(vDirection, v_vDirection))
0264 continue;
0265 if (!this->checkModuleDirections(wDirection, v_wDirection))
0266 continue;
0267 if (!this->checkModulePositions(posR, v_posR))
0268 continue;
0269 if (!this->checkModulePositions(posPhi, v_posPhi))
0270 continue;
0271 if (!this->checkModulePositions(posEta, v_posEta))
0272 continue;
0273 if (!this->checkModulePositions(posX, v_posX))
0274 continue;
0275 if (!this->checkModulePositions(posY, v_posY))
0276 continue;
0277 if (!this->checkModulePositions(posZ, v_posZ))
0278 continue;
0279
0280 tkSector.v_rawId.push_back(rawId);
0281 bool moduleSelected(false);
0282 for (auto const& i_rawId : allSectors.v_rawId) {
0283 if (rawId == i_rawId)
0284 moduleSelected = true;
0285 }
0286 if (!moduleSelected)
0287 allSectors.v_rawId.push_back(rawId);
0288 }
0289
0290
0291 bool isPixel(false);
0292 bool isStrip(false);
0293 for (auto const& i_rawId : tkSector.v_rawId) {
0294 switch (m_tkTreeVar_[i_rawId].subdetId) {
0295 case PixelSubdetector::PixelBarrel:
0296 case PixelSubdetector::PixelEndcap:
0297 isPixel = true;
0298 break;
0299 case StripSubdetector::TIB:
0300 case StripSubdetector::TOB:
0301 case StripSubdetector::TID:
0302 case StripSubdetector::TEC:
0303 isStrip = true;
0304 break;
0305 }
0306 }
0307
0308 if (isPixel && isStrip) {
0309 edm::LogError("SectorBuilder")
0310 << "Incorrect Sector Definition: there are pixel and strip modules within one sector"
0311 << "\n... sector selection is not applied, sector " << sectorCounter << " is not built";
0312 continue;
0313 }
0314 tkSector.isPixel = isPixel;
0315
0316 m_tkSector_[sectorCounter] = tkSector;
0317 edm::LogInfo("SectorBuilder") << "There are " << tkSector.v_rawId.size() << " Modules in Sector " << sectorCounter;
0318 }
0319 noSectors = sectorCounter;
0320 return;
0321 }
0322
0323
0324 bool ApeTreeCreateDefault::checkIntervalsForSectors(const unsigned int sectorCounter,
0325 const std::vector<double>& v_id) const {
0326 if (v_id.empty())
0327 return true;
0328 if (v_id.size() % 2 == 1) {
0329 edm::LogError("SectorBuilder")
0330 << "Incorrect Sector Definition: Position Vectors need even number of arguments (Intervals)"
0331 << "\n... sector selection is not applied, sector " << sectorCounter << " is not built";
0332 return false;
0333 }
0334 int entry(0);
0335 double intervalBegin(999.);
0336 for (auto const& i_id : v_id) {
0337 ++entry;
0338 if (entry % 2 == 1)
0339 intervalBegin = i_id;
0340 if (entry % 2 == 0 && intervalBegin > i_id) {
0341 edm::LogError("SectorBuilder") << "Incorrect Sector Definition (Position Vector Intervals): \t" << intervalBegin
0342 << " is bigger than " << i_id << " but is expected to be smaller"
0343 << "\n... sector selection is not applied, sector " << sectorCounter
0344 << " is not built";
0345 return false;
0346 }
0347 }
0348 return true;
0349 }
0350
0351 bool ApeTreeCreateDefault::checkModuleIds(const unsigned int id, const std::vector<unsigned int>& v_id) const {
0352 if (v_id.empty())
0353 return true;
0354 for (auto const& i_id : v_id) {
0355 if (id == i_id)
0356 return true;
0357 }
0358 return false;
0359 }
0360
0361 bool ApeTreeCreateDefault::checkModuleBools(const bool id, const std::vector<unsigned int>& v_id) const {
0362 if (v_id.empty())
0363 return true;
0364 for (auto const& i_id : v_id) {
0365 if (1 == i_id && id)
0366 return true;
0367 if (2 == i_id && !id)
0368 return true;
0369 }
0370 return false;
0371 }
0372
0373 bool ApeTreeCreateDefault::checkModuleDirections(const int id, const std::vector<int>& v_id) const {
0374 if (v_id.empty())
0375 return true;
0376 for (auto const& i_id : v_id) {
0377 if (id == i_id)
0378 return true;
0379 }
0380 return false;
0381 }
0382
0383 bool ApeTreeCreateDefault::checkModulePositions(const float id, const std::vector<double>& v_id) const {
0384 if (v_id.empty())
0385 return true;
0386 int entry(0);
0387 double intervalBegin(999.);
0388 for (auto const& i_id : v_id) {
0389 ++entry;
0390 if (entry % 2 == 1)
0391 intervalBegin = i_id;
0392 if (entry % 2 == 0 && id >= intervalBegin && id < i_id)
0393 return true;
0394 }
0395 return false;
0396 }
0397
0398
0399 void ApeTreeCreateDefault::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0400
0401
0402
0403 const AlignmentErrorsExtended* alignmentErrors = &iSetup.getData(alignmentErrorToken_);
0404
0405
0406 const std::string defaultFileName(resultFile_);
0407 TFile* defaultFile = new TFile(defaultFileName.c_str(), "RECREATE");
0408
0409
0410 TTree* defaultTreeX(nullptr);
0411 TTree* defaultTreeY(nullptr);
0412 defaultFile->GetObject("iterTreeX;1", defaultTreeX);
0413 defaultFile->GetObject("iterTreeY;1", defaultTreeY);
0414
0415 TTree* sectorNameTree(nullptr);
0416 defaultFile->GetObject("nameTree;1", sectorNameTree);
0417
0418 edm::LogInfo("DefaultAPETree") << "APE Tree is being created";
0419 defaultTreeX = new TTree("iterTreeX", "Tree for default APE x values from GT");
0420 defaultTreeY = new TTree("iterTreeY", "Tree for default APE y values from GT");
0421 sectorNameTree = new TTree("nameTree", "Tree with names of sectors");
0422
0423
0424 std::vector<double*> a_defaultSectorX;
0425 std::vector<double*> a_defaultSectorY;
0426
0427 std::vector<std::string*> a_sectorName;
0428 for (auto const& i_sector : m_tkSector_) {
0429 const unsigned int iSector(i_sector.first);
0430 const bool pixelSector(i_sector.second.isPixel);
0431
0432 a_defaultSectorX.push_back(new double(-99.));
0433 a_defaultSectorY.push_back(new double(-99.));
0434 a_sectorName.push_back(new std::string(i_sector.second.name));
0435
0436 std::stringstream ss_sector;
0437 std::stringstream ss_sectorSuffixed;
0438 ss_sector << "Ape_Sector_" << iSector;
0439
0440 ss_sectorSuffixed << ss_sector.str() << "/D";
0441 defaultTreeX->Branch(ss_sector.str().c_str(), &(*a_defaultSectorX[iSector - 1]), ss_sectorSuffixed.str().c_str());
0442
0443 if (pixelSector) {
0444 defaultTreeY->Branch(ss_sector.str().c_str(), &(*a_defaultSectorY[iSector - 1]), ss_sectorSuffixed.str().c_str());
0445 }
0446 sectorNameTree->Branch(ss_sector.str().c_str(), &(*a_sectorName[iSector - 1]), 32000, 00);
0447 }
0448
0449
0450
0451 for (auto& i_sector : m_tkSector_) {
0452 double defaultApeX(0.);
0453 double defaultApeY(0.);
0454 unsigned int nModules(0);
0455 for (auto const& i_rawId : i_sector.second.v_rawId) {
0456 std::vector<AlignTransformErrorExtended> alignErrors = alignmentErrors->m_alignError;
0457 for (auto const& i_alignError : alignErrors) {
0458 if (i_rawId == i_alignError.rawId()) {
0459 CLHEP::HepSymMatrix errMatrix = i_alignError.matrix();
0460 defaultApeX += errMatrix[0][0];
0461 defaultApeY += errMatrix[1][1];
0462 nModules++;
0463 }
0464 }
0465 }
0466 *a_defaultSectorX[i_sector.first - 1] = defaultApeX / nModules;
0467 *a_defaultSectorY[i_sector.first - 1] = defaultApeY / nModules;
0468 }
0469
0470 sectorNameTree->Fill();
0471 sectorNameTree->Write("nameTree");
0472 defaultTreeX->Fill();
0473 defaultTreeX->Write("iterTreeX");
0474 defaultTreeY->Fill();
0475 defaultTreeY->Write("iterTreeY");
0476
0477 defaultFile->Close();
0478 delete defaultFile;
0479 for (unsigned int i = 0; i < a_defaultSectorX.size(); i++) {
0480 delete a_defaultSectorX[i];
0481 delete a_defaultSectorY[i];
0482 delete a_sectorName[i];
0483 }
0484 }
0485
0486
0487 void ApeTreeCreateDefault::beginJob() { this->sectorBuilder(); }
0488
0489
0490 void ApeTreeCreateDefault::endJob() {}
0491
0492 void ApeTreeCreateDefault::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0493 edm::ParameterSetDescription desc;
0494 edm::ParameterSetDescription sector;
0495
0496 std::vector<unsigned> emptyUnsignedIntVector;
0497 std::vector<int> emptyIntVector;
0498 std::vector<double> emptyDoubleVector;
0499 sector.add<std::string>("name", "default");
0500 sector.add<std::vector<unsigned>>("rawId", emptyUnsignedIntVector);
0501 sector.add<std::vector<unsigned>>("subdetId", emptyUnsignedIntVector);
0502 sector.add<std::vector<unsigned>>("layer", emptyUnsignedIntVector);
0503 sector.add<std::vector<unsigned>>("side", emptyUnsignedIntVector);
0504 sector.add<std::vector<unsigned>>("half", emptyUnsignedIntVector);
0505 sector.add<std::vector<unsigned>>("rod", emptyUnsignedIntVector);
0506 sector.add<std::vector<unsigned>>("ring", emptyUnsignedIntVector);
0507 sector.add<std::vector<unsigned>>("petal", emptyUnsignedIntVector);
0508 sector.add<std::vector<unsigned>>("blade", emptyUnsignedIntVector);
0509 sector.add<std::vector<unsigned>>("panel", emptyUnsignedIntVector);
0510 sector.add<std::vector<unsigned>>("outerInner", emptyUnsignedIntVector);
0511 sector.add<std::vector<unsigned>>("module", emptyUnsignedIntVector);
0512 sector.add<std::vector<unsigned>>("nStrips", emptyUnsignedIntVector);
0513 sector.add<std::vector<unsigned>>("isDoubleSide", emptyUnsignedIntVector);
0514 sector.add<std::vector<unsigned>>("isRPhi", emptyUnsignedIntVector);
0515 sector.add<std::vector<unsigned>>("isStereo", emptyUnsignedIntVector);
0516 sector.add<std::vector<int>>("uDirection", emptyIntVector);
0517 sector.add<std::vector<int>>("vDirection", emptyIntVector);
0518 sector.add<std::vector<int>>("wDirection", emptyIntVector);
0519 sector.add<std::vector<double>>("posR", emptyDoubleVector);
0520 sector.add<std::vector<double>>("posPhi", emptyDoubleVector);
0521 sector.add<std::vector<double>>("posEta", emptyDoubleVector);
0522 sector.add<std::vector<double>>("posX", emptyDoubleVector);
0523 sector.add<std::vector<double>>("posY", emptyDoubleVector);
0524 sector.add<std::vector<double>>("posZ", emptyDoubleVector);
0525
0526 desc.add<std::string>("resultFile", "defaultAPE.root");
0527 desc.add<std::string>("trackerTreeFile");
0528 desc.addVPSet("sectors", sector);
0529
0530 descriptions.add("apeTreeCreateDefault", desc);
0531 }
0532
0533
0534 DEFINE_FWK_MODULE(ApeTreeCreateDefault);