ApeTreeCreateDefault.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    DefaultApeTree
// Class:      DefaultApeTree
// 
//
// Original Author:  Marius Teroerde (code from ApeEstimator.cc and
//                   ApeEstimatorSummary.cc)
//         Created:  Tue Nov 14 11:43 CET 2017
//
//


// system include files
#include <fstream>
#include <sstream>

// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/one/EDAnalyzer.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"

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

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Utilities/interface/EDMException.h"

#include "CLHEP/Matrix/SymMatrix.h"

#include "CondFormats/AlignmentRecord/interface/TrackerAlignmentErrorExtendedRcd.h"
#include "CondFormats/Alignment/interface/AlignmentErrorsExtended.h"

//...............
#include "Alignment/APEEstimation/interface/TrackerSectorStruct.h"
#include "Alignment/APEEstimation/interface/ReducedTrackerTreeVariables.h"

#include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
#include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
#include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
#include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
#include "DataFormats/SiStripCluster/interface/SiStripCluster.h"


#include "TString.h"
#include "TFile.h"
#include "TDirectory.h"
#include "TTree.h"
#include "TMath.h"
//
// class declaration
//

class ApeTreeCreateDefault : public edm::one::EDAnalyzer<> {
  public:
    explicit ApeTreeCreateDefault(const edm::ParameterSet&);
    ~ApeTreeCreateDefault() override;
    
    static void fillDescriptions(edm::ConfigurationDescriptions & descriptions);
  
  private:
    void beginJob() override ;
    void analyze(const edm::Event&, const edm::EventSetup&) override;
    void endJob() override ;
    
    void sectorBuilder();
    bool checkIntervalsForSectors(const unsigned int sectorCounter, const std::vector<double>&)const;
    bool checkModuleIds(const unsigned int, const std::vector<unsigned int>&)const;
    bool checkModuleBools(const bool, const std::vector<unsigned int>&)const;
    bool checkModuleDirections(const int, const std::vector<int>&)const;
    bool checkModulePositions(const float, const std::vector<double>&)const;
    
    // ----------member data ---------------------------
    const std::string resultFile_;
    const std::string trackerTreeFile_;
    const std::vector<edm::ParameterSet> sectors_;
    
    std::map<unsigned int, TrackerSectorStruct> m_tkSector_;
    std::map<unsigned int, ReducedTrackerTreeVariables> m_tkTreeVar_;
    unsigned int noSectors;
};              

//
// constants, enums and typedefs
//

//
// static data member definitions
//

//
// constructors and destructor
//
ApeTreeCreateDefault::ApeTreeCreateDefault(const edm::ParameterSet& iConfig):
  resultFile_(iConfig.getParameter<std::string>("resultFile")),
  trackerTreeFile_(iConfig.getParameter<std::string>("trackerTreeFile")),
  sectors_(iConfig.getParameter<std::vector<edm::ParameterSet> >("sectors"))
{
}


ApeTreeCreateDefault::~ApeTreeCreateDefault()
{
}


//
// member functions
//
void
ApeTreeCreateDefault::sectorBuilder()
{
  // Same procedure as in ApeEstimator.cc
  TFile* tkTreeFile(TFile::Open(trackerTreeFile_.c_str()));
  if(tkTreeFile){
    edm::LogInfo("SectorBuilder")<<"TrackerTreeFile OK";
  }else{
    edm::LogError("SectorBuilder")<<"TrackerTreeFile not found";
    return;
  }
  TTree* tkTree(nullptr);
  tkTreeFile->GetObject("TrackerTreeGenerator/TrackerTree/TrackerTree",tkTree);
  if(tkTree){
    edm::LogInfo("SectorBuilder")<<"TrackerTree OK";
  }else{
    edm::LogError("SectorBuilder")<<"TrackerTree not found in file";
    return;
  }
  
  unsigned int rawId(999), subdetId(999), layer(999), side(999), half(999), rod(999), ring(999), petal(999),
         blade(999), panel(999), outerInner(999), module(999), nStrips(999);
  bool isDoubleSide(false), isRPhi(false), isStereo(false);
  int uDirection(999), vDirection(999), wDirection(999);
  float posR(999.F), posPhi(999.F), posEta(999.F), posX(999.F), posY(999.F), posZ(999.F); 
  
  tkTree->SetBranchAddress("RawId", &rawId);
  tkTree->SetBranchAddress("SubdetId", &subdetId);
  tkTree->SetBranchAddress("Layer", &layer);
  tkTree->SetBranchAddress("Side", &side);
  tkTree->SetBranchAddress("Half", &half);
  tkTree->SetBranchAddress("Rod", &rod);
  tkTree->SetBranchAddress("Ring", &ring);
  tkTree->SetBranchAddress("Petal", &petal);
  tkTree->SetBranchAddress("Blade", &blade);
  tkTree->SetBranchAddress("Panel", &panel);
  tkTree->SetBranchAddress("OuterInner", &outerInner);
  tkTree->SetBranchAddress("Module", &module);
  tkTree->SetBranchAddress("NStrips", &nStrips);
  tkTree->SetBranchAddress("IsDoubleSide", &isDoubleSide);
  tkTree->SetBranchAddress("IsRPhi", &isRPhi);
  tkTree->SetBranchAddress("IsStereo", &isStereo);
  tkTree->SetBranchAddress("UDirection", &uDirection);
  tkTree->SetBranchAddress("VDirection", &vDirection);
  tkTree->SetBranchAddress("WDirection", &wDirection);
  tkTree->SetBranchAddress("PosR", &posR);
  tkTree->SetBranchAddress("PosPhi", &posPhi);
  tkTree->SetBranchAddress("PosEta", &posEta);
  tkTree->SetBranchAddress("PosX", &posX);
  tkTree->SetBranchAddress("PosY", &posY);
  tkTree->SetBranchAddress("PosZ", &posZ);
  
  int nModules(tkTree->GetEntries());
  TrackerSectorStruct allSectors;
  
  //Loop over all Sectors
  unsigned int sectorCounter(0);
  std::vector<edm::ParameterSet> v_sectorDef(sectors_);
  edm::LogInfo("SectorBuilder")<<"There are "<<v_sectorDef.size()<<" Sectors defined";

  for(auto const & parSet : v_sectorDef){
    ++sectorCounter;
    const std::string& sectorName(parSet.getParameter<std::string>("name"));
    std::vector<unsigned int> v_rawId(parSet.getParameter<std::vector<unsigned int> >("rawId")),
                              v_subdetId(parSet.getParameter<std::vector<unsigned int> >("subdetId")),
                              v_layer(parSet.getParameter<std::vector<unsigned int> >("layer")),
                              v_side(parSet.getParameter<std::vector<unsigned int> >("side")),
                              v_half(parSet.getParameter<std::vector<unsigned int> >("half")),
                              v_rod(parSet.getParameter<std::vector<unsigned int> >("rod")),
                              v_ring(parSet.getParameter<std::vector<unsigned int> >("ring")),
                              v_petal(parSet.getParameter<std::vector<unsigned int> >("petal")),
                              v_blade(parSet.getParameter<std::vector<unsigned int> >("blade")),
                              v_panel(parSet.getParameter<std::vector<unsigned int> >("panel")),
                              v_outerInner(parSet.getParameter<std::vector<unsigned int> >("outerInner")),
                              v_module(parSet.getParameter<std::vector<unsigned int> >("module")),
                              v_nStrips(parSet.getParameter<std::vector<unsigned int> >("nStrips")),
                              v_isDoubleSide(parSet.getParameter<std::vector<unsigned int> >("isDoubleSide")),
                              v_isRPhi(parSet.getParameter<std::vector<unsigned int> >("isRPhi")),
                              v_isStereo(parSet.getParameter<std::vector<unsigned int> >("isStereo"));
    std::vector<int>  v_uDirection(parSet.getParameter<std::vector<int> >("uDirection")),
                      v_vDirection(parSet.getParameter<std::vector<int> >("vDirection")),
                      v_wDirection(parSet.getParameter<std::vector<int> >("wDirection"));
    std::vector<double> v_posR(parSet.getParameter<std::vector<double> >("posR")),
                        v_posPhi(parSet.getParameter<std::vector<double> >("posPhi")),
                        v_posEta(parSet.getParameter<std::vector<double> >("posEta")),
                        v_posX(parSet.getParameter<std::vector<double> >("posX")),
                        v_posY(parSet.getParameter<std::vector<double> >("posY")),
                        v_posZ(parSet.getParameter<std::vector<double> >("posZ"));
    
    if(!this->checkIntervalsForSectors(sectorCounter,v_posR) || !this->checkIntervalsForSectors(sectorCounter,v_posPhi) ||
       !this->checkIntervalsForSectors(sectorCounter,v_posEta) || !this->checkIntervalsForSectors(sectorCounter,v_posX) ||
       !this->checkIntervalsForSectors(sectorCounter,v_posY)   || !this->checkIntervalsForSectors(sectorCounter,v_posZ)){
      continue;
    }
    
    TrackerSectorStruct tkSector;
    tkSector.name = sectorName;
    
    ReducedTrackerTreeVariables tkTreeVar;
    
    //Loop over all Modules
    for(int module = 0; module < nModules; ++module){
      tkTree->GetEntry(module);
      
      if(sectorCounter==1){
        tkTreeVar.subdetId = subdetId;
        tkTreeVar.nStrips = nStrips;
        tkTreeVar.uDirection = uDirection;
        tkTreeVar.vDirection = vDirection;
        tkTreeVar.wDirection = wDirection;
        m_tkTreeVar_[rawId] = tkTreeVar;
      }
      //Check if modules from Sector builder equal those from TrackerTree
      if(!this->checkModuleIds(rawId,v_rawId))continue;
      if(!this->checkModuleIds(subdetId,v_subdetId))continue;
      if(!this->checkModuleIds(layer,v_layer))continue;
      if(!this->checkModuleIds(side,v_side))continue;
      if(!this->checkModuleIds(half,v_half))continue;
      if(!this->checkModuleIds(rod,v_rod))continue;
      if(!this->checkModuleIds(ring,v_ring))continue;
      if(!this->checkModuleIds(petal,v_petal))continue;
      if(!this->checkModuleIds(blade,v_blade))continue;
      if(!this->checkModuleIds(panel,v_panel))continue;
      if(!this->checkModuleIds(outerInner,v_outerInner))continue;
      if(!this->checkModuleIds(module,v_module))continue;
      if(!this->checkModuleIds(nStrips,v_nStrips))continue;
      if(!this->checkModuleBools(isDoubleSide,v_isDoubleSide))continue;
      if(!this->checkModuleBools(isRPhi,v_isRPhi))continue;
      if(!this->checkModuleBools(isStereo,v_isStereo))continue;
      if(!this->checkModuleDirections(uDirection,v_uDirection))continue;
      if(!this->checkModuleDirections(vDirection,v_vDirection))continue;
      if(!this->checkModuleDirections(wDirection,v_wDirection))continue;
      if(!this->checkModulePositions(posR,v_posR))continue;
      if(!this->checkModulePositions(posPhi,v_posPhi))continue;
      if(!this->checkModulePositions(posEta,v_posEta))continue;
      if(!this->checkModulePositions(posX,v_posX))continue;
      if(!this->checkModulePositions(posY,v_posY))continue;
      if(!this->checkModulePositions(posZ,v_posZ))continue;
      
      tkSector.v_rawId.push_back(rawId);
      bool moduleSelected(false);
      for(auto const & i_rawId : allSectors.v_rawId){
        if(rawId == i_rawId)moduleSelected = true;
      }
      if(!moduleSelected)allSectors.v_rawId.push_back(rawId);
    }
    
    // Stops you from combining pixel and strip detector into one sector
    bool isPixel(false);
    bool isStrip(false);
    for(auto const & i_rawId : tkSector.v_rawId){
      switch (m_tkTreeVar_[i_rawId].subdetId) {
        case PixelSubdetector::PixelBarrel:
        case PixelSubdetector::PixelEndcap:
          isPixel = true;
          break;
        case StripSubdetector::TIB:
        case StripSubdetector::TOB:
        case StripSubdetector::TID:
        case StripSubdetector::TEC:
          isStrip = true;
          break;
      }
    }

    if(isPixel && isStrip){
      edm::LogError("SectorBuilder")<<"Incorrect Sector Definition: there are pixel and strip modules within one sector"
                                     <<"\n... sector selection is not applied, sector "<<sectorCounter<<" is not built";
      continue;
    }
    tkSector.isPixel = isPixel;
    
    m_tkSector_[sectorCounter] = tkSector;
    edm::LogInfo("SectorBuilder")<<"There are "<<tkSector.v_rawId.size()<<" Modules in Sector "<<sectorCounter;
  }
  noSectors = sectorCounter;
  return;
}


// Checking methods copied from ApeEstimator.cc
bool
ApeTreeCreateDefault::checkIntervalsForSectors(const unsigned int sectorCounter, const std::vector<double>& v_id)const
{
    
  if(v_id.empty())return true;
  if(v_id.size()%2==1){
    edm::LogError("SectorBuilder")<<"Incorrect Sector Definition: Position Vectors need even number of arguments (Intervals)"
                                     <<"\n... sector selection is not applied, sector "<<sectorCounter<<" is not built";
    return false;
  }
  int entry(0); double intervalBegin(999.);
  for(auto const & i_id : v_id){
    ++entry;
    if(entry%2==1) intervalBegin = i_id;
    if(entry%2==0 && intervalBegin > i_id){
      edm::LogError("SectorBuilder")<<"Incorrect Sector Definition (Position Vector Intervals): \t"
                                    <<intervalBegin<<" is bigger than "<<i_id<<" but is expected to be smaller"
                                    <<"\n... sector selection is not applied, sector "<<sectorCounter<<" is not built";
      return false;
    }
  }
  return true;
}

bool
ApeTreeCreateDefault::checkModuleIds(const unsigned int id, const std::vector<unsigned int>& v_id)const
{
    
  if(v_id.empty())return true;
  for(auto const & i_id : v_id){
    if(id==i_id)return true;
  }
  return false;
}

bool
ApeTreeCreateDefault::checkModuleBools(const bool id, const std::vector<unsigned int>& v_id)const
{
    
  if(v_id.empty())return true;
  for(auto const & i_id : v_id){
    if(1==i_id && id)return true;
    if(2==i_id && !id)return true;
  }
  return false;
}

bool
ApeTreeCreateDefault::checkModuleDirections(const int id, const std::vector<int>& v_id)const
{
    
  if(v_id.empty())return true;
  for(auto const & i_id : v_id){
    if(id==i_id)return true;
  }
  return false;
}

bool
ApeTreeCreateDefault::checkModulePositions(const float id, const std::vector<double>& v_id)const
{
    
  if(v_id.empty())return true;
  int entry(0); double intervalBegin(999.);
  for(auto const & i_id : v_id){
    ++entry;
    if(entry%2==1)intervalBegin = i_id;
    if(entry%2==0 && id>=intervalBegin && id<i_id)return true;
  }
  return false;
}



// ------------ method called to for each event  ------------
void
ApeTreeCreateDefault::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{
  // Same procedure as in ApeEstimatorSummary.cc minus reading of baseline tree
  
  // Load APEs from the GT and write them to root files similar to the ones from calculateAPE()
  edm::ESHandle<AlignmentErrorsExtended> alignmentErrors;
  iSetup.get<TrackerAlignmentErrorExtendedRcd>().get(alignmentErrors);  
  
  // Set up root file for default APE values
  const std::string defaultFileName(resultFile_);
  TFile* defaultFile = new TFile(defaultFileName.c_str(),"RECREATE");
  
  // Naming in the root files has to be iterTreeX to be consistent for the plotting tool
  TTree* defaultTreeX(nullptr);
  TTree* defaultTreeY(nullptr);
  defaultFile->GetObject("iterTreeX;1",defaultTreeX);
  defaultFile->GetObject("iterTreeY;1",defaultTreeY);
  // The same for TTree containing the names of the sectors (no additional check, since always handled exactly as defaultTree)
  TTree* sectorNameTree(nullptr);
  defaultFile->GetObject("nameTree;1",sectorNameTree);
  
  edm::LogInfo("DefaultAPETree")<<"APE Tree is being created";
  defaultTreeX = new TTree("iterTreeX","Tree for default APE x values from GT");
  defaultTreeY = new TTree("iterTreeY","Tree for default APE y values from GT");     
  sectorNameTree = new TTree("nameTree","Tree with names of sectors");
  
    
  // Assign the information stored in the trees to arrays
  std::vector<double*> a_defaultSectorX;
  std::vector<double*> a_defaultSectorY;
  
  std::vector<std::string*> a_sectorName;
  for(auto const & i_sector : m_tkSector_){
    const unsigned int iSector(i_sector.first);
    const bool pixelSector(i_sector.second.isPixel);
    
    a_defaultSectorX.push_back(new double(-99.));
    a_defaultSectorY.push_back(new double(-99.));
    a_sectorName.push_back(new std::string(i_sector.second.name));
    
    std::stringstream ss_sector;
    std::stringstream ss_sectorSuffixed;
    ss_sector << "Ape_Sector_" << iSector;
  
    ss_sectorSuffixed << ss_sector.str() << "/D";
    defaultTreeX->Branch(ss_sector.str().c_str(), &(*a_defaultSectorX[iSector-1]), ss_sectorSuffixed.str().c_str());
    
    if(pixelSector){
      defaultTreeY->Branch(ss_sector.str().c_str(), &(*a_defaultSectorY[iSector-1]), ss_sectorSuffixed.str().c_str());
    }
    sectorNameTree->Branch(ss_sector.str().c_str(), &(*a_sectorName[iSector-1]), 32000, 00);  
  }

  
  // Loop over sectors for getting default APE
  
  for(auto & i_sector : m_tkSector_){
    
    double defaultApeX(0.);
    double defaultApeY(0.);
    unsigned int nModules(0);
    for(auto const & i_rawId : i_sector.second.v_rawId){
      std::vector<AlignTransformErrorExtended> alignErrors = alignmentErrors->m_alignError;
      for(auto const & i_alignError : alignErrors){
        if(i_rawId ==  i_alignError.rawId()){
          CLHEP::HepSymMatrix errMatrix = i_alignError.matrix();
          defaultApeX += errMatrix[0][0];
          defaultApeY += errMatrix[1][1];
          nModules++;
        }
      }
  }
  *a_defaultSectorX[i_sector.first-1] = defaultApeX/nModules;
  *a_defaultSectorY[i_sector.first-1] = defaultApeY/nModules;
  }
  
  
  sectorNameTree->Fill();
  sectorNameTree->Write("nameTree");
  defaultTreeX->Fill();
  defaultTreeX->Write("iterTreeX");
  defaultTreeY->Fill();
  defaultTreeY->Write("iterTreeY");
  
  defaultFile->Close(); 
  delete defaultFile;
  for(unsigned int i = 0; i < a_defaultSectorX.size(); i++){
    delete a_defaultSectorX[i];
    delete a_defaultSectorY[i];
    delete a_sectorName[i];
  }
}
  



// ------------ method called once each job just before starting event loop  ------------
void 
ApeTreeCreateDefault::beginJob()
{  
  this->sectorBuilder();
}

// ------------ method called once each job just after ending the event loop  ------------
void 
ApeTreeCreateDefault::endJob() 
{   
}

void
ApeTreeCreateDefault::fillDescriptions(edm::ConfigurationDescriptions & descriptions)
{
  edm::ParameterSetDescription desc;
  edm::ParameterSetDescription sector;
  
  std::vector<unsigned> emptyUnsignedIntVector;
  std::vector<int> emptyIntVector;
  std::vector<double> emptyDoubleVector;
  sector.add<std::string>("name", "default");
  sector.add<std::vector<unsigned>>("rawId", emptyUnsignedIntVector);
  sector.add<std::vector<unsigned>>("subdetId", emptyUnsignedIntVector);
  sector.add<std::vector<unsigned>>("layer", emptyUnsignedIntVector);
  sector.add<std::vector<unsigned>>("side", emptyUnsignedIntVector);
  sector.add<std::vector<unsigned>>("half", emptyUnsignedIntVector);
  sector.add<std::vector<unsigned>>("rod", emptyUnsignedIntVector);
  sector.add<std::vector<unsigned>>("ring", emptyUnsignedIntVector);
  sector.add<std::vector<unsigned>>("petal", emptyUnsignedIntVector);
  sector.add<std::vector<unsigned>>("blade", emptyUnsignedIntVector);
  sector.add<std::vector<unsigned>>("panel", emptyUnsignedIntVector);
  sector.add<std::vector<unsigned>>("outerInner", emptyUnsignedIntVector);
  sector.add<std::vector<unsigned>>("module", emptyUnsignedIntVector);
  sector.add<std::vector<unsigned>>("nStrips", emptyUnsignedIntVector);
  sector.add<std::vector<unsigned>>("isDoubleSide", emptyUnsignedIntVector);
  sector.add<std::vector<unsigned>>("isRPhi", emptyUnsignedIntVector);
  sector.add<std::vector<unsigned>>("isStereo", emptyUnsignedIntVector);
  sector.add<std::vector<int>>("uDirection", emptyIntVector);
  sector.add<std::vector<int>>("vDirection", emptyIntVector);
  sector.add<std::vector<int>>("wDirection", emptyIntVector);
  sector.add<std::vector<double>>("posR", emptyDoubleVector);
  sector.add<std::vector<double>>("posPhi", emptyDoubleVector);
  sector.add<std::vector<double>>("posEta", emptyDoubleVector);
  sector.add<std::vector<double>>("posX", emptyDoubleVector);
  sector.add<std::vector<double>>("posY", emptyDoubleVector);
  sector.add<std::vector<double>>("posZ", emptyDoubleVector);

  desc.add<std::string>("resultFile", "defaultAPE.root");
  desc.add<std::string>("trackerTreeFile");
  desc.addVPSet("sectors", sector);
 
  descriptions.add("apeTreeCreateDefault", desc);
}

//define this as a plug-in
DEFINE_FWK_MODULE(ApeTreeCreateDefault);