TrackerRemapper/plugins/SiPixelPhase1Analyzer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    DQM/TrackerRemapper
0004 // Class:      SiPixelPhase1Analyzer
0005 //
0006
0007 #include "DQM/TrackerRemapper/interface/SiPixelPhase1Analyzer.h"
0008
0009 using namespace std;
0010 using namespace edm;
0011
0012 SiPixelPhase1Analyzer::SiPixelPhase1Analyzer(const edm::ParameterSet& iConfig)
0013     : geomToken_(esConsumes()),
0014       topoToken_(esConsumes()),
0015       opMode(static_cast<OperationMode>(iConfig.getUntrackedParameter<unsigned int>("opMode"))),
0016       debugFileName(iConfig.getUntrackedParameter<string>("debugFileName")),
0017       firstEvent(true),
0018       rootFileHandle(nullptr),
0019       isBarrelSource(iConfig.getUntrackedParameter<vector<unsigned>>("isBarrelSource")),
0020       analazedRootFileName(iConfig.getUntrackedParameter<vector<string>>("remapRootFileName")),
0021       pathToHistograms(iConfig.getUntrackedParameter<vector<string>>("pathToHistograms")),
0022       baseHistogramName(iConfig.getUntrackedParameter<vector<string>>("baseHistogramName")) {
0023 #ifdef DEBUG_MODE
0024   debugFile = std::ofstream(debugFileName.c_str(), std::ofstream::out);
0025 #endif
0026   usesResource("TFileService");
0027
0028   orthoProjectionMatrix.BuildOrthographicMatrix(1.0f, -1.0f, 1.0f, -1.0f, -10.0f, 10.0f);
0029
0030   switch (opMode) {
0031     case MODE_ANALYZE:
0032
0033       tracksToken = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("src"));
0034
0035       analazedRootFileName.clear();
0036
0037       pathToHistograms.clear();
0038       pathToHistograms.push_back("RecHits/");
0039
0040       baseHistogramName.clear();
0041       baseHistogramName.push_back("RecHits");
0042
0043       break;
0044     case MODE_REMAP:
0045       break;
0046     default:
0047       break;
0048   }
0049 }
0050
0051 SiPixelPhase1Analyzer::~SiPixelPhase1Analyzer() {
0052   for (auto& i : bins) {
0053     delete i.second;
0054   }
0055
0056   for (auto& i : binsSummary) {
0057     delete i.second;
0058   }
0059
0060 #ifdef DEBUG_MODE
0061   debugFile.close();
0062 #endif
0063 }
0064
0065 // ------------ method called for each event  ------------
0066 void SiPixelPhase1Analyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0067   const auto& theTrackerGeometry = iSetup.getData(geomToken_);
0068   const auto& tt = &iSetup.getData(topoToken_);
0069
0070   if (firstEvent) {
0071     /////////////////////////////////
0072     BookHistograms();
0073     /////////////////////////////////
0074     BookBins(theTrackerGeometry, tt);
0075     /////////////////////////////////
0076     if (opMode == MODE_REMAP) {
0077       FillBins(nullptr, theTrackerGeometry, tt);
0078     }
0079     /////////////////////////////////
0080     firstEvent = false;
0081   }
0082   if (opMode == MODE_ANALYZE) {
0083     edm::Handle<reco::TrackCollection> tracks;
0084     iEvent.getByToken(tracksToken, tracks);
0085     if (!tracks.isValid()) {
0086       LogInfo("Analyzer") << "reco::TrackCollection not found... Aborting...\n";
0087       return;
0088     }
0089     FillBins(&tracks, theTrackerGeometry, tt);
0090   }
0091   // debugFile << "SiPixelPhase1Analyzer::analyze() - Event " << iEvent.run() << "/" << iEvent.id().event() << endl;
0092 }
0093
0094 void SiPixelPhase1Analyzer::BookHistograms() {
0095   // ghost object <initializes> FileService (without it, it crashes and creation of directories would not be possible)
0096   TObject* ghostObj = fs->make<TH2Poly>("ghost", "ghost", -1, 1, -1, 1);
0097
0098   TDirectory* topDir = fs->getBareDirectory();
0099   topDir->cd();
0100
0101 #ifdef DEBUG_MODE
0102   debugFile << "Full path: " << fs->fullPath() << endl << endl;
0103 #endif
0104   string histName;
0105   for (unsigned j = 0; j < baseHistogramName.size(); ++j) {
0106     string currentHistoName = baseHistogramName[j];
0107
0108     TDirectory* currentDir = topDir->mkdir(currentHistoName.c_str());
0109     currentDir->cd();
0110
0111     if (opMode == MODE_REMAP) {
0112       if (isBarrelSource[j]) {
0113         BookBarrelHistograms(currentDir, currentHistoName);
0114       } else {
0115         BookForwardHistograms(currentDir, currentHistoName);
0116       }
0117     } else {
0118       BookBarrelHistograms(currentDir, currentHistoName);
0119       BookForwardHistograms(currentDir, currentHistoName);
0120     }
0121
0122     topDir->cd();
0123   }
0124
0125   ghostObj->Delete();  //not needed any more
0126 }
0127
0128 void SiPixelPhase1Analyzer::BookBarrelHistograms(TDirectory* currentDir, const string& currentHistoName) {
0129   string histName;
0130   TH2Poly* th2p;
0131
0132 #ifdef DEBUG_MODE
0133   TH2* th2;
0134 #endif
0135
0136   for (unsigned i = 0; i < 4; ++i) {
0137     histName = "barrel_layer_";
0138
0139     th2p = new TH2Poly((histName + std::to_string(i + 1)).c_str(), "PXBMap", -15.0, 15.0, 0.0, 5.0);
0140
0141     th2p->SetFloat();
0142
0143     th2p->GetXaxis()->SetTitle("z [cm]");
0144     th2p->GetYaxis()->SetTitle("ladder");
0145
0146 #ifdef DEBUG_MODE
0147     th2p->SetOption("COLZ 0 TEXT");
0148 #else
0149     th2p->SetOption("COLZ L");
0150 #endif
0151
0152     currentDir->Add(th2p);
0153     th2PolyBarrel[currentHistoName].push_back(th2p);
0154
0155 #ifdef DEBUG_MODE
0156     if (opMode == MODE_ANALYZE) {
0157       th2 = new TH2I((histName + std::to_string(i + 1) + "_DEBUG").c_str(),
0158                      "position",
0159                      3000,
0160                      -30.0f,
0161                      30.0f,
0162                      1000,
0163                      -4.5f * (i + 1),
0164                      4.5f * (i + 1));
0165
0166       th2->GetXaxis()->SetTitle("z [cm]");
0167       th2->GetYaxis()->SetTitle("-x [?]");
0168
0169       th2->SetOption("COLZ 0 TEXT");
0170
0171       currentDir->Add(th2);
0172       th2PolyBarrelDebug[currentHistoName].push_back(th2);
0173     }
0174 #endif
0175   }
0176
0177   th2p = new TH2Poly("barrel_summary", "PXBMap", -5.0, 5.0, 0.0, 5.0);
0178   th2p->SetFloat();
0179
0180   th2p->GetXaxis()->SetTitle("");
0181   th2p->GetYaxis()->SetTitle("~ladder");
0182
0183   th2p->SetOption("COLZ L");
0184
0185   currentDir->Add(th2p);
0186   th2PolyBarrelSummary[currentHistoName] = th2p;
0187 }
0188
0189 void SiPixelPhase1Analyzer::BookForwardHistograms(TDirectory* currentDir, const string& currentHistoName) {
0190   string histName;
0191   TH2Poly* th2p;
0192 #ifdef DEBUG_MODE
0193   TH2* th2;
0194 #endif
0195
0196   for (unsigned side = 1; side <= 2; ++side) {
0197     for (unsigned disk = 1; disk <= 3; ++disk) {
0198       histName = "forward_disk_";
0199
0200       th2p = new TH2Poly((histName + std::to_string((side == 1 ? -(int(disk)) : (int)disk))).c_str(),
0201                          "PXFMap",
0202                          -15.0,
0203                          15.0,
0204                          -15.0,
0205                          15.0);
0206
0207       th2p->SetFloat();
0208
0209       th2p->GetXaxis()->SetTitle("x [cm]");
0210       th2p->GetYaxis()->SetTitle("y [cm]");
0211
0212 #ifdef DEBUG_MODE
0213       th2p->SetOption("COLZ 0 TEXT");
0214 #else
0215       th2p->SetOption("COLZ L");
0216 #endif
0217       currentDir->Add(th2p);
0218       pxfTh2PolyForward[currentHistoName].push_back(th2p);
0219
0220 #ifdef DEBUG_MODE
0221       if (opMode == MODE_ANALYZE) {
0222         th2 = new TH2I((histName + std::to_string((side == 1 ? -(int(disk)) : (int)disk)) + "_DEBUG").c_str(),
0223                        "position",
0224                        1000,
0225                        -15.0f,
0226                        15.0f,
0227                        1000,
0228                        -15.0f,
0229                        15.0f);
0230
0231         th2->GetXaxis()->SetTitle("x [cm]");
0232         th2->GetYaxis()->SetTitle("y [cm]");
0233
0234         th2->SetOption("COLZ 0 TEXT");
0235
0236         currentDir->Add(th2);
0237         pxfTh2PolyForwardDebug[currentHistoName].push_back(th2);
0238       }
0239 #endif
0240     }
0241   }
0242
0243   th2p = new TH2Poly("forward_summary", "PXFMap", -40.0, 50.0, -20.0, 90.0);
0244   th2p->SetFloat();
0245
0246   th2p->GetXaxis()->SetTitle("");
0247   th2p->GetYaxis()->SetTitle("");
0248
0249 #ifdef DEBUG_MODE
0250   th2p->SetOption("COLZ 0 TEXT");
0251 #else
0252   th2p->SetOption("COLZ L");
0253 #endif
0254
0255   currentDir->Add(th2p);
0256   pxfTh2PolyForwardSummary[currentHistoName] = th2p;
0257 }
0258
0259 void SiPixelPhase1Analyzer::BookBins(const TrackerGeometry& theTrackerGeometry, const TrackerTopology* tt) {
0260   BookBarrelBins(theTrackerGeometry, tt);
0261   BookForwardBins(theTrackerGeometry, tt);
0262
0263 #ifdef DEBUG_MODE
0264   SaveDetectorVertices(tt);
0265 #endif
0266 }
0267
0268 void SiPixelPhase1Analyzer::BookBarrelBins(const TrackerGeometry& theTrackerGeometry, const TrackerTopology* tt) {
0269   TrackingGeometry::DetContainer pxb = theTrackerGeometry.detsPXB();
0270 #ifdef DEBUG_MODE
0271   debugFile << "There are " << pxb.size() << " detector elements in the PXB." << endl;
0272 #endif
0273   for (auto& i : pxb) {
0274     const GeomDet* det = i;
0275
0276     PXBDetId id = det->geographicalId();
0277
0278     Local2DPoint origin;
0279     GlobalPoint p = det->surface().toGlobal(origin);
0280
0281     int layer = tt->pxbLayer(id);
0282     int ladder = tt->pxbLadder(id);
0283
0284 #ifdef DEBUG_MODE
0285     int module = tt->pxbModule(id);
0286     PixelBarrelName pixelBarrelName(id, tt, true);
0287     SaveDetectorData(
0288         true, id.rawId(), pixelBarrelName.shell(), pixelBarrelName.layerName(), pixelBarrelName.ladderName());
0289 #endif
0290
0291 #ifdef DEBUG_MODE
0292     float r = sqrt(p.x() * p.x() + p.y() * p.y());
0293
0294     debugFile << "Layer: " << layer << "\tLadder: " << ladder << "\tModule: " << module << "\t(x, y, z, r2): (" << p.x()
0295               << ", " << p.y() << ", " << p.z() << ", " << r << ")" << endl;
0296 #endif
0297
0298     const Bounds& b = (det->surface().bounds());
0299     float bl = b.length();
0300
0301 #ifdef DEBUG_MODE
0302     float bw = b.width();
0303     float bt = b.thickness();
0304
0305     debugFile << "Length: " << bl << "\tWidth: " << bw << "\tThickness: " << bt << endl;
0306 #endif
0307
0308     float vertX[] = {p.z() - bl * 0.5f, p.z() + bl * 0.5f, p.z() + bl * 0.5f, p.z() - bl * 0.5f, p.z() - bl * 0.5f};
0309     float vertY[] = {(ladder - 1.0f), (ladder - 1.0f), (float)ladder, (float)ladder, (ladder - 1.0f)};
0310
0311     bins[id.rawId()] = new TGraph(5, vertX, vertY);
0312     bins[id.rawId()]->SetName(TString::Format("%u", id.rawId()));
0313
0314     // Summary plot
0315     for (unsigned k = 0; k < 5; ++k) {
0316       vertX[k] += ((layer == 2 || layer == 3) ? 0.0f : -60.0f);
0317       vertY[k] += ((layer > 2) ? 30.0f : 0.0f);
0318     }
0319
0320     binsSummary[id.rawId()] = new TGraph(5, vertX, vertY);
0321     binsSummary[id.rawId()]->SetName(TString::Format("%u", id.rawId()));
0322
0323     for (unsigned nameNum = 0; nameNum < baseHistogramName.size(); ++nameNum) {
0324       if (isBarrelSource[nameNum] || opMode == MODE_ANALYZE) {
0325         const string& strName = baseHistogramName[nameNum];
0326         th2PolyBarrel[strName][layer - 1]->AddBin(bins[id.rawId()]->Clone());
0327         th2PolyBarrelSummary[strName]->AddBin(binsSummary[id.rawId()]->Clone());
0328       }
0329     }
0330   }
0331 }
0332
0333 void SiPixelPhase1Analyzer::BookForwardBins(const TrackerGeometry& theTrackerGeometry, const TrackerTopology* tt) {
0334   TrackingGeometry::DetContainer pxf = theTrackerGeometry.detsPXF();
0335 #ifdef DEBUG_MODE
0336   debugFile << "There are " << pxf.size() << " detector elements in the PXF." << endl;
0337 #endif
0338   bool firstForwardElem = true;
0339
0340   float elemWidth = 1.0f;
0341   float elemLength = 1.0f;
0342
0343   // FIRST PASS -> MAKE MAP OF CORRESPONDING ELEMENTS (BLADES ON BOTH PANELS)
0344   for (auto& i : pxf) {
0345     const GeomDet* det = i;
0346
0347     PXFDetId id = det->geographicalId();
0348
0349     Local2DPoint origin;
0350     GlobalPoint p = det->surface().toGlobal(origin);
0351
0352     int panel = tt->pxfPanel(id);
0353     int side = tt->side(id);   //tt->pxfSide(id);
0354     int disk = tt->layer(id);  //tt->pxfDisk(id);
0355     int blade = tt->pxfBlade(id);
0356
0357 #ifdef DEBUG_MODE
0358     int module = tt->module(id);  //tt->pxfModule(id);
0359     PixelEndcapName pixelEndcapName(id, tt, true);
0360     SaveDetectorData(
0361         false, id.rawId(), pixelEndcapName.halfCylinder(), pixelEndcapName.diskName(), pixelEndcapName.bladeName());
0362 #endif
0363
0364 #ifdef DEBUG_MODE
0365     float r = sqrt(p.x() * p.x() + p.y() * p.y());
0366
0367     debugFile << "Panel: " << panel << "\tSide: " << side << "\tDisk: " << disk << "\tBlade: " << blade
0368               << "\tModule: " << module << "\t(x, y, z, r): (" << p.x() << ", " << p.y() << ", " << p.z() << ", " << r
0369               << ")" << endl;
0370 #endif
0371     if (firstForwardElem) {
0372       const Bounds& b = det->surface().bounds();  //const RectangularPlaneBounds& b
0373
0374       elemLength = b.length();
0375       elemWidth = b.width();
0376
0377       firstForwardElem = false;
0378     }
0379
0380     const auto& rot = det->rotation();
0381
0382     mat4 transMat(
0383         rot.xx(), rot.xy(), rot.xz(), rot.yx(), rot.yy(), rot.yz(), rot.zx(), rot.zy(), rot.zz(), p.x(), p.y(), p.z());
0384
0385     mapOfComplementaryElements[CODE_FORWARD(side, disk, blade)].mat[panel - 1] = transMat;
0386     mapOfComplementaryElements[CODE_FORWARD(side, disk, blade)].rawId[panel - 1] = id.rawId();
0387   }
0388
0389   // SECOND PASS -> USE INFORMATION FROM MAP TO MAKE GEOMETRIC OBJECTS (BINS)
0390   for (auto& i : mapOfComplementaryElements) {
0391     // decode side&disk from the map key
0392     unsigned side = i.first & 0xF;
0393     unsigned disk = (i.first & 0xF0) >> 4;
0394     // unsigned blade = (i.first & 0xFF00) >> 8;
0395
0396     unsigned mapIdx = disk + (side - 1) * 3 - 1;
0397
0398     // normal vectors of elements point to the (almost) opposite direction, so correction is needed before interploation (probably not 100% correct but fast)
0399     i.second.mat[1].data[0] = -i.second.mat[1].data[0];
0400     i.second.mat[1].data[1] = -i.second.mat[1].data[1];
0401     i.second.mat[1].data[2] = -i.second.mat[1].data[2];
0402
0403     i.second.mat[1].data[6] = -i.second.mat[1].data[6];
0404     i.second.mat[1].data[7] = -i.second.mat[1].data[7];
0405     i.second.mat[1].data[8] = -i.second.mat[1].data[8];
0406
0407     mat4 meanTransform = (i.second.mat[0] + i.second.mat[1]) * 0.5f;
0408     // mat4 meanTransform = i.second.mat[0];
0409
0410     static const float baseVertX[4] = {-elemWidth * 0.8f, -elemWidth * 0.5f, elemWidth * 0.8f, elemWidth * 0.5f};
0411     static const float baseVertY[4] = {
0412         elemLength * 0.38f, -elemLength * 0.38f, elemLength * 0.38f, -elemLength * 0.38f};
0413
0414     float vertXPanel[2][4], vertYPanel[2][4];
0415     float vertIn[3], vertOut[3];
0416
0417     /*
0418         (1)  __________________ (3)
0419              \             /
0420               \           /
0421                \             /
0422                 \       /
0423               (2)\_________/(4)
0424
0425              - division line: (2) - (3)
0426
0427              - panel 1: triangle of lower area (2, 3, 4)
0428              - panel 2: triangle of bigger area (1, 2, 3)
0429     */
0430
0431     // obtain transformed vertices
0432     for (unsigned j = 0; j < 4; ++j) {
0433       vertIn[0] = baseVertX[j];
0434       vertIn[1] = baseVertY[j];
0435       vertIn[2] = 0.0f;
0436
0437       meanTransform.MulVec(vertIn, vertOut);
0438       std::swap(vertIn, vertOut);
0439       orthoProjectionMatrix.MulVec(vertIn, vertOut);
0440
0441       // vertical flip
0442       vertOut[0] = -vertOut[0];  // so that inner elements have positive x-coordinate
0443
0444       if (j > 0) {
0445         vertXPanel[0][j - 1] = vertOut[0];
0446         vertYPanel[0][j - 1] = vertOut[1];  // for panel 2
0447       }
0448       if (j < 3) {
0449         vertXPanel[1][j] = vertOut[0];
0450         vertYPanel[1][j] = vertOut[1];  // for panel 1
0451       }
0452     }
0453
0454     for (unsigned j = 0; j < 2; ++j) {
0455       vertXPanel[j][3] = vertXPanel[j][0];
0456       vertYPanel[j][3] = vertYPanel[j][0];
0457
0458       bins[i.second.rawId[j]] = new TGraph(4, vertXPanel[j], vertYPanel[j]);
0459       bins[i.second.rawId[j]]->SetName(TString::Format("%u", i.second.rawId[j]));
0460
0461       // for (auto strName: baseHistogramName)
0462       // {
0463       // pxfTh2PolyForward[strName][mapIdx]->AddBin(bins[i.second.rawId[j]]->Clone());
0464       // }
0465
0466       // Summary plot
0467       for (unsigned k = 0; k < 4; ++k) {
0468         vertXPanel[j][k] += (float(side) - 1.5f) * 40.0f;
0469         vertYPanel[j][k] += (disk - 1) * 35.0f;
0470       }
0471
0472       binsSummary[i.second.rawId[j]] = new TGraph(4, vertXPanel[j], vertYPanel[j]);
0473       binsSummary[i.second.rawId[j]]->SetName(TString::Format("%u", i.second.rawId[j]));
0474
0475       for (unsigned nameNum = 0; nameNum < baseHistogramName.size(); ++nameNum) {
0476         if (!isBarrelSource[nameNum] || opMode == MODE_ANALYZE) {
0477           const string& strName = baseHistogramName[nameNum];
0478           pxfTh2PolyForward[strName][mapIdx]->AddBin(bins[i.second.rawId[j]]->Clone());
0479           pxfTh2PolyForwardSummary[strName]->AddBin(binsSummary[i.second.rawId[j]]->Clone());
0480         }
0481       }
0482     }
0483   }
0484 }
0485
0486 void SiPixelPhase1Analyzer::SaveDetectorVertices(const TrackerTopology* tt) {
0487   vector<std::ofstream*> verticesFiles[2];
0488   for (unsigned i = 0; i < 4; ++i) {
0489     std::ofstream* f = new std::ofstream(("vertices_barrel_" + std::to_string(i + 1)).c_str(), std::ofstream::out);
0490
0491     verticesFiles[0].push_back(f);
0492   }
0493
0494   for (unsigned side = 1; side <= 2; ++side) {
0495     for (unsigned disk = 1; disk <= 3; ++disk) {
0496       std::ofstream* f = new std::ofstream(
0497           ("vertices_forward_" + std::to_string((side == 1 ? -(int(disk)) : (int)disk))).c_str(), std::ofstream::out);
0498
0499       verticesFiles[1].push_back(f);
0500     }
0501   }
0502
0503   for (auto& bin : bins) {
0504     unsigned rawId = bin.first;
0505     DetId id(rawId);
0506     unsigned subdetId = id.subdetId();
0507
0508     if (subdetId != PixelSubdetector::PixelBarrel && subdetId != PixelSubdetector::PixelEndcap)
0509       continue;
0510
0511     double* vertX = bin.second->GetX();
0512     double* vertY = bin.second->GetY();
0513
0514     if (subdetId == PixelSubdetector::PixelBarrel) {
0515       PXBDetId pxbId(rawId);
0516       unsigned layer = tt->pxbLayer(pxbId);
0517       string onlineName = PixelBarrelName(pxbId, tt, true).name();
0518
0519       *(verticesFiles[0][layer - 1]) << rawId << " " << onlineName << " \"";
0520       for (unsigned i = 0; i < 4; ++i) {
0521         *(verticesFiles[0][layer - 1]) << vertX[i] << "," << vertY[i];
0522         if (i == 3)
0523           *(verticesFiles[0][layer - 1]) << "\"\n";
0524         else
0525           *(verticesFiles[0][layer - 1]) << " ";
0526       }
0527     } else {
0528       PXFDetId pxfId(rawId);
0529       unsigned side = tt->pxfSide(pxfId);
0530       unsigned disk = tt->pxfDisk(pxfId);
0531       string onlineName = PixelEndcapName(pxfId, tt, true).name();
0532       unsigned mapIdx = disk + (side - 1) * 3 - 1;
0533
0534       *(verticesFiles[1][mapIdx]) << rawId << " " << onlineName << " \"";
0535       for (unsigned i = 0; i < 3; ++i) {
0536         *(verticesFiles[1][mapIdx]) << vertX[i] << "," << vertY[i];
0537
0538         if (i == 2)
0539           *(verticesFiles[1][mapIdx]) << "\"\n";
0540         else
0541           *(verticesFiles[1][mapIdx]) << " ";
0542       }
0543     }
0544   }
0545
0546   for (unsigned i = 0; i < 2; ++i) {
0547     for (auto& j : verticesFiles[i]) {
0548       j->close();
0549       delete j;
0550     }
0551   }
0552 }
0553
0554 void SiPixelPhase1Analyzer::FillBins(edm::Handle<reco::TrackCollection>* tracks,
0555                                      const TrackerGeometry& theTrackerGeometry,
0556                                      const TrackerTopology* tt) {
0557   switch (opMode) {
0558     case MODE_ANALYZE:
0559       for (auto const& track : *(*tracks)) {
0560         auto recHitsBegin = track.recHitsBegin();
0561         for (unsigned i = 0; i < track.recHitsSize(); ++i) {
0562           auto recHit = *(recHitsBegin + i);
0563           if (!recHit->isValid())
0564             continue;
0565
0566           DetId id = recHit->geographicalId();
0567           unsigned subdetId = id.subdetId();
0568
0569           if (subdetId != PixelSubdetector::PixelBarrel && subdetId != PixelSubdetector::PixelEndcap)
0570             continue;
0571
0572           const PixelGeomDetUnit* geomdetunit =
0573               dynamic_cast<const PixelGeomDetUnit*>(theTrackerGeometry.idToDet(id));  // theTrackerGeometry ?????
0574           //const PixelTopology& topol = geomdetunit->specificTopology();
0575
0576           LocalPoint localPoint = recHit->localPosition();
0577           GlobalPoint globalPoint = geomdetunit->surface().toGlobal(localPoint);
0578
0579           if (subdetId == PixelSubdetector::PixelBarrel)
0580             FillBarrelBinsAnalyze(theTrackerGeometry, tt, id.rawId(), globalPoint);
0581           else
0582             FillForwardBinsAnalyze(theTrackerGeometry, tt, id.rawId(), globalPoint);
0583         }
0584       }
0585       break;
0586     case MODE_REMAP:
0587       FillBarrelBinsRemap(theTrackerGeometry, tt);
0588       FillForwardBinsRemap(theTrackerGeometry, tt);
0589       break;
0590     default:
0591       break;
0592   }
0593 }
0594
0595 void SiPixelPhase1Analyzer::FillBarrelBinsAnalyze(const TrackerGeometry& theTrackerGeometry,
0596                                                   const TrackerTopology* tt,
0597                                                   unsigned rawId,
0598                                                   const GlobalPoint& globalPoint) {
0599   for (unsigned nameNum = 0; nameNum < baseHistogramName.size(); ++nameNum) {
0600     string strName = baseHistogramName[nameNum];
0601
0602     PXBDetId id(rawId);
0603
0604     int layer = tt->pxbLayer(id);
0605
0606     th2PolyBarrel[strName][layer - 1]->Fill(TString::Format("%u", rawId), 1);
0607     th2PolyBarrelSummary[strName]->Fill(TString::Format("%u", rawId), 1);
0608 #ifdef DEBUG_MODE
0609     th2PolyBarrelDebug[strName][layer - 1]->Fill((globalPoint.y() < 0 ? globalPoint.z() + 0.5f : globalPoint.z()),
0610                                                  -globalPoint.x(),
0611                                                  (globalPoint.y() < 0 ? -1 : 1));
0612 #endif
0613   }
0614 }
0615
0616 void SiPixelPhase1Analyzer::FillForwardBinsAnalyze(const TrackerGeometry& theTrackerGeometry,
0617                                                    const TrackerTopology* tt,
0618                                                    unsigned rawId,
0619                                                    const GlobalPoint& globalPoint) {
0620   for (unsigned nameNum = 0; nameNum < baseHistogramName.size(); ++nameNum) {
0621     string strName = baseHistogramName[nameNum];
0622
0623     PXFDetId id(rawId);
0624
0625     int side = tt->side(id);   //tt->pxfSide(id);
0626     int disk = tt->layer(id);  //tt->pxfDisk(id);
0627     // int blade = tt->pxfBlade(id);
0628     unsigned mapIdx = disk + (side - 1) * 3 - 1;
0629
0630     pxfTh2PolyForward[strName][mapIdx]->Fill(TString::Format("%u", rawId), 1);
0631     pxfTh2PolyForwardSummary[strName]->Fill(TString::Format("%u", rawId), 1);
0632
0633 #ifdef DEBUG_MODE
0634     pxfTh2PolyForwardDebug[strName][mapIdx]->Fill(globalPoint.x(), globalPoint.y(), 1);
0635 #endif
0636   }
0637 }
0638
0639 void SiPixelPhase1Analyzer::FillBarrelBinsRemap(const TrackerGeometry& theTrackerGeometry, const TrackerTopology* tt) {
0640   rootFileHandle = new TFile(analazedRootFileName[0].c_str());
0641
0642   if (!rootFileHandle) {
0643     LogInfo("Analyzer") << "Could not open file: " << analazedRootFileName[0] << "..." << endl;
0644     return;
0645   }
0646
0647 #ifdef DEBUG_MODE
0648   rootFileHandle->ls();
0649   LogInfo("Analyzer") << "\n\n";
0650   rootFileHandle->pwd();
0651   LogInfo("Analyzer") << "\n\n";
0652 #endif
0653
0654   for (unsigned nameNum = 0; nameNum < baseHistogramName.size(); ++nameNum) {
0655     if (!isBarrelSource[nameNum])
0656       continue;
0657
0658     // if (pathToHistograms[nameNum][pathToHistograms.size() - 1] != '/') pathToHistograms[nameNum] += "/";
0659     string baseHistogramNameWithPath = pathToHistograms[nameNum] + baseHistogramName[nameNum] + "_";
0660
0661     const TProfile2D* handles[4];
0662 #ifndef DEBUG_MODE
0663     const TProfile2D* h;
0664 #endif
0665     bool problemWithHandles = false;
0666
0667     for (unsigned i = 0; i < 4; ++i) {
0668       string fullFileName = (baseHistogramNameWithPath + std::to_string(i + 1) + ";1 ");
0669       handles[i] = (TProfile2D*)rootFileHandle->Get(fullFileName.c_str());
0670       if (!handles[i]) {
0671         problemWithHandles = true;
0672         LogInfo("Analyzer") << "Histogram: " << fullFileName << " does not exist!\n";
0673
0674         break;
0675       }
0676     }
0677
0678     if (!problemWithHandles) {
0679       LogInfo("Analyzer") << "\nInput histograms: " << baseHistogramNameWithPath << " opened successfully\n";
0680
0681       //Add original histograms to this file
0682
0683       TDirectory* currDir = fs->getBareDirectory()->GetDirectory(baseHistogramName[nameNum].c_str());
0684       currDir->cd();
0685
0686       for (unsigned i = 0; i < 4; ++i) {
0687         currDir->Add(handles[i]->Clone());
0688       }
0689
0690       TrackingGeometry::DetContainer pxb = theTrackerGeometry.detsPXB();
0691 #ifdef DEBUG_MODE
0692       debugFile << "There are " << pxb.size() << " detector elements in the PXB." << endl;
0693 #endif
0694       for (auto& i : pxb) {
0695         const GeomDet* det = i;
0696
0697         PXBDetId id = det->geographicalId();
0698         unsigned rawId = id.rawId();
0699
0700         int module = tt->pxbModule(id);
0701         //int ladder = tt->pxbLadder(id);
0702         int layer = tt->pxbLayer(id);
0703
0704         int signedOnlineModule = module - 4;
0705         if (signedOnlineModule <= 0)
0706           --signedOnlineModule;
0707
0708         PixelBarrelName pixelBarrelName = PixelBarrelName(id, tt, true);
0709         int onlineShell = pixelBarrelName.shell();
0710
0711         int signedOnlineLadder = ((onlineShell & 1) ? -pixelBarrelName.ladderName() : pixelBarrelName.ladderName());
0712         string strName = baseHistogramName[nameNum];
0713
0714 #ifdef DEBUG_MODE
0715         th2PolyBarrel[strName][layer - 1]->Fill(TString::Format("%u", rawId), signedOnlineLadder);
0716         th2PolyBarrelSummary[strName]->Fill(TString::Format("%u", rawId), signedOnlineLadder);
0717 #else
0718         h = handles[layer - 1];
0719         unsigned nx = h->GetNbinsX();
0720         unsigned ny = h->GetNbinsY();
0721         unsigned binX = signedOnlineModule + ((nx + 1) >> 1);
0722         unsigned binY = (signedOnlineLadder) + ((ny + 1) >> 1);
0723         double val = h->GetBinContent(binX, binY);
0724         th2PolyBarrel[strName][layer - 1]->Fill(TString::Format("%u", rawId), val);
0725         th2PolyBarrelSummary[strName]->Fill(TString::Format("%u", rawId), val);
0726 #endif
0727       }
0728     }
0729   }
0730
0731   rootFileHandle->Close();
0732   delete rootFileHandle;
0733 }
0734
0735 void SiPixelPhase1Analyzer::FillForwardBinsRemap(const TrackerGeometry& theTrackerGeometry, const TrackerTopology* tt) {
0736   rootFileHandle = new TFile(analazedRootFileName[0].c_str());
0737
0738   if (!rootFileHandle) {
0739     return;
0740   }
0741
0742   TrackingGeometry::DetContainer pxf = theTrackerGeometry.detsPXF();
0743
0744 #ifdef DEBUG_MODE
0745   rootFileHandle->ls();
0746   LogInfo("Analyzer") << "\n\n";
0747   rootFileHandle->pwd();
0748   LogInfo("Analyzer") << "\n\n";
0749 #endif
0750
0751   for (unsigned nameNum = 0; nameNum < baseHistogramName.size(); ++nameNum) {
0752     if (isBarrelSource[nameNum])
0753       continue;
0754
0755     string baseHistogramNameWithPath = pathToHistograms[nameNum] + baseHistogramName[nameNum] + "_";
0756
0757     const TProfile2D* h_1 = (TProfile2D*)rootFileHandle->Get((baseHistogramNameWithPath + "1;1 ").c_str());
0758     const TProfile2D* h_2 = (TProfile2D*)rootFileHandle->Get((baseHistogramNameWithPath + "2;1 ").c_str());
0759 #ifndef DEBUG_MODE
0760     const TProfile2D* h;
0761 #endif
0762     if (h_2 && h_1) {
0763       LogInfo("Analyzer") << "\nInput histograms: " << baseHistogramNameWithPath << " opened successfully\n";
0764
0765       //Add original histograms to this file
0766       TDirectory* currDir = fs->getBareDirectory()->GetDirectory(baseHistogramName[nameNum].c_str());
0767       currDir->cd();
0768       currDir->Add(h_1->Clone());
0769       currDir->Add(h_2->Clone());
0770
0771       for (auto& i : pxf) {
0772         const GeomDet* det = i;
0773
0774         PXFDetId id = det->geographicalId();
0775
0776         int side = tt->side(id);   //tt->pxfSide(id);
0777         int disk = tt->layer(id);  //tt->pxfDisk(id);
0778
0779         unsigned rawId = id.rawId();
0780         PixelEndcapName pixelEndcapName = PixelEndcapName(PXFDetId(rawId), tt, true);
0781
0782 #ifdef DEBUG_MODE
0783         int blade = tt->pxfBlade(id);
0784 #else
0785         int onlineBlade = pixelEndcapName.bladeName();
0786         bool isInnerOnlineBlade = !(pixelEndcapName.halfCylinder() & 1);  // inner -> blade > 0 (?)
0787
0788         int signedOnlineBlade = (isInnerOnlineBlade) ? onlineBlade : -onlineBlade;
0789         int signedDisk = (side == 2) ? disk : -disk;
0790
0791         int pannel = pixelEndcapName.pannelName() - 1;
0792
0793 #endif
0794         unsigned mapIdx = disk + (side - 1) * 3 - 1;
0795         string strName = baseHistogramName[nameNum];
0796
0797 #ifdef DEBUG_MODE
0798         pxfTh2PolyForward[strName][mapIdx]->Fill(TString::Format("%u", rawId), blade);
0799         pxfTh2PolyForwardSummary[strName]->Fill(TString::Format("%u", rawId), blade);
0800 #else
0801         if (pixelEndcapName.ringName() == 1)
0802           h = h_1;
0803         else
0804           h = h_2;
0805         // ---- REMAP (Online -> Offline)
0806         unsigned nx = h->GetNbinsX();
0807         unsigned ny = h->GetNbinsY();
0808         unsigned binX = signedDisk + ((nx + 1) >> 1);
0809         unsigned binY = (signedOnlineBlade * 2) + (ny >> 1);
0810         double val = h->GetBinContent(binX, binY + pannel);
0811         pxfTh2PolyForward[strName][mapIdx]->Fill(TString::Format("%u", rawId), val);
0812         pxfTh2PolyForwardSummary[strName]->Fill(TString::Format("%u", rawId), val);
0813 #endif
0814       }
0815     }
0816   }
0817
0818   rootFileHandle->Close();
0819   delete rootFileHandle;
0820 }
0821
0822 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0823 void SiPixelPhase1Analyzer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0824   //The following says we do not know what parameters are allowed so do no validation
0825   // Please change this to state exactly what you do use, even if it is no parameters
0826   edm::ParameterSetDescription desc;
0827   desc.setComment(
0828       "Creates TH2Poly Pixel Tracker maps by either analyzing the event or remapping exising DQM historams");
0829   desc.add<edm::InputTag>("src", edm::InputTag("generalTracks"));
0830   desc.addUntracked<unsigned int>("opMode", 1);
0831   desc.addUntracked<std::string>("debugFileName", "debug.txt");
0832   desc.addUntracked<std::vector<unsigned int>>("isBarrelSource", {0, 0, 1});
0833   desc.addUntracked<std::vector<std::string>>("remapRootFileName", {"dqmFile.root"});
0834   desc.addUntracked<std::vector<std::string>>(
0835       "pathToHistograms",
0836       {"DQMData/Run 1/PixelPhase1/Run summary/Phase1_MechanicalView/PXForward/",
0837        "DQMData/Run 1/PixelPhase1/Run summary/Phase1_MechanicalView/PXForward/",
0838        "DQMData/Run 1/PixelPhase1/Run summary/Phase1_MechanicalView/PXBarrel/"});
0839   desc.addUntracked<std::vector<std::string>>("baseHistogramName",
0840                                               {"num_clusters_per_PXDisk_per_SignedBladePanel_PXRing",
0841                                                "num_digis_per_PXDisk_per_SignedBladePanel_PXRing",
0842                                                "num_digis_per_SignedModule_per_SignedLadder_PXLayer"});
0843   descriptions.addWithDefaultLabel(desc);
0844 }
0845
0846 //define this as a plug-in
0847 DEFINE_FWK_MODULE(SiPixelPhase1Analyzer);