---

StandaloneMerger.cc

Go to the documentation of this file.

#include "Alignment/MuonStandaloneAlgorithm/plugins/StandaloneMerger.h"

// Collaborating Class Header
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/ESHandle.h"

#include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
#include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
#include "Geometry/CommonDetUnit/interface/GeomDet.h"
#include "DataFormats/DetId/interface/DetId.h"
#include "DataFormats/MuonDetId/interface/DTChamberId.h"
#include "DataFormats/MuonDetId/interface/CSCDetId.h"

#include "Alignment/MuonStandaloneAlgorithm/interface/AlignmentDetector.h"
#include "Alignment/MuonStandaloneAlgorithm/interface/PointForAlignment.h"
#include "Alignment/MuonStandaloneAlgorithm/interface/TrackForAlignment.h"
#include "Alignment/MuonStandaloneAlgorithm/interface/AlignmentDetectorCollection.h"
#include "Alignment/MuonStandaloneAlgorithm/interface/TrackForAlignmentCollection.h"

#include "FWCore/ServiceRegistry/interface/Service.h"

//#include "CondCore/DBOutputService/interface/PoolDBOutputService.h"

#include "TFile.h"
#include "TH1F.h"
#include <stdio.h>
#include <fstream>

using namespace std;
using namespace edm;



StandaloneMerger::StandaloneMerger(const ParameterSet& pset){
  
  maxFiles = pset.getParameter<int>("maxFiles");
  surveyFileName = pset.getParameter<std::string>("surveyFileName");
  outputFile = pset.getParameter<std::string>("outputFileName");
  submatrices = pset.getParameter<std::string>("subMatrixDir");
}


StandaloneMerger::~StandaloneMerger(){
}

void StandaloneMerger::updateMatrix(TMatrixD *Cm, TMatrixD *bm, TMatrixD *in, TMatrixD *numberOfTracks) {
 

  /*if(Cm->E2Norm() > 1e7 || bm->E2Norm() > 1.0e5) {
    std::cout << "Bad Matrix. Skipping file" << std::endl;
    return;
  }*/

  //Cm->Print("CMatrix");
  //bm->Print("bMatrix");


  TMatrixD realIndex(in->GetNrows(), 1);

  N(0,0) += (*numberOfTracks)(0,0);

  for(int indexCounter = 0; indexCounter < in->GetNrows(); ++indexCounter) {
    bool oldIndex = false;
    for(int theIndexCounter = 0; theIndexCounter < theIndex.GetNrows(); ++theIndexCounter) {
      if(theIndex(theIndexCounter, 0) == (*in)(indexCounter, 0)) {
        realIndex(indexCounter, 0) = theIndexCounter;
        oldIndex = true;
        break;
      }
    }
    if(oldIndex == false) {
      theIndex.ResizeTo(theIndex.GetNrows()+1,1);
      C.ResizeTo(C.GetNrows()+NDOFAlign, C.GetNcols()+NDOFAlign);
      b.ResizeTo(b.GetNrows()+NDOFAlign, 1);
      theIndex(theIndex.GetNrows()-1, 0) = (*in)(indexCounter, 0);
      realIndex(indexCounter, 0) = theIndex.GetNrows()-1;
    }
  }
  
  for(int param = 0; param < NDOFAlign; param++) {
    for(int rowCounter = 0; rowCounter < in->GetNrows(); rowCounter++) {
      for(int param2 = 0; param2 < NDOFAlign; param2++) {
        for(int columnCounter = 0; columnCounter < in->GetNrows(); columnCounter++) {
          C((int)realIndex(rowCounter,0)*NDOFAlign+param, (int)realIndex(columnCounter,0)*NDOFAlign+param2) += (*Cm)(rowCounter*NDOFAlign+param, columnCounter*NDOFAlign+param2);
          if(C((int)realIndex(rowCounter,0)*NDOFAlign+param, (int)realIndex(columnCounter,0)*NDOFAlign+param2) > 1e100) {
            std::cout << "Error. Anormal values of matrix: " << std::endl;
          }
        }
      }
      b((int)realIndex(rowCounter,0)*NDOFAlign+param,0) += (*bm)(rowCounter*NDOFAlign+param, 0);
    }
  }
}


void StandaloneMerger::mergeMatrix() {

  for(int counterOfFiles = 0; counterOfFiles < maxFiles; ++counterOfFiles) {
        
    char name[100];
    sprintf(name, "%s/SubMatrices_%d.root", submatrices.c_str(), counterOfFiles);
    
    f = new TFile(name, "READ");
    if(f->IsZombie())  continue;
    std::cout << "Analyzing: " << name << std::endl;
    TMatrixD *Cm = (TMatrixD *)f->Get("CMatrix");
    TMatrixD *bm = (TMatrixD *)f->Get("bMatrix");
    TMatrixD *in = (TMatrixD *)f->Get("Index");
    TMatrixD *numberOfTracks = (TMatrixD *)f->Get("numberOfTracks");   
    updateMatrix(Cm, bm, in, numberOfTracks);
    
    f->Close();
    delete f;
    delete Cm;
    delete bm;
    delete in;
  }
}


void StandaloneMerger::beginJob(const EventSetup& eventSetup){
  
  
  //Set the size of Matrizes to 0
  C.ResizeTo(0,0);
  b.ResizeTo(0,0);
  theIndex.ResizeTo(0,0);  
  N.ResizeTo(1,1);

  
  mergeMatrix(); // This method reads matrizes in C and alignables vector
  surveyAdd();
  CMatrixInvert.ResizeTo(C.GetNrows(), C.GetNcols());
  CMatrixInvert = TotalMatrix;
  CMatrixInvert.Invert();
  Delta.ResizeTo(C.GetNrows(), 1);
  Delta = CMatrixInvert*TotalVector;  
  TMatrixD unit, unit2;
  unit.ResizeTo(C.GetNrows(), C.GetNcols());
  unit2.ResizeTo(C.GetNrows(), C.GetNcols());
  unit.UnitMatrix();
  unit2 = unit-CMatrixInvert*TotalMatrix;
  std::cout << "Inversion check: " << unit2.E2Norm() << std::endl;
  
}



void StandaloneMerger::BuildDelta() {

 Delta.ResizeTo(C.GetNrows(), 1);
 for(int j = C.GetNrows()-1; j >= 0; --j) {
   double acum = 0.0;
   for(int i = C.GetNrows()-1; i > j; --i) {
      acum += R(j, i)*Delta(i,0);
   }
   Delta(j,0) = (IndependentVector(j,0)-acum)/R(j,j);
 }
}




bool StandaloneMerger::CalculateSolution() {
  
  //Decomp = new TDecompQRH(CMatrixInvert, 1e-30);
  Q.ResizeTo(CMatrixInvert.GetNrows(), CMatrixInvert.GetNcols());
  R.ResizeTo(CMatrixInvert.GetNrows(), CMatrixInvert.GetNcols());
  IndependentVector.ResizeTo(CMatrixInvert.GetNrows(), 1);
  return true;
  /*if(Decomp->Decompose()) {
    Q = Decomp->GetQ();
    R = Decomp->GetR();
    IndependentVector = Q*TotalVector;
    BuildDelta();   
    //delete Decomp;
    return true;
  } else {
    //delete Decomp;
    return false;
  }*/
}


void StandaloneMerger::endJob(){
  
  //Output root file is opened for writting
  
  theFile = TFile::Open(outputFile.c_str(), "RECREATE");
  theFile->cd();
  
  C.Write("CMatrix");
  b.Write("bMatrix");
  Delta.Write("Delta");
  CMatrixInvert.Write("CMatrixInvert");
  theIndex.Write("Index");
  N.Write("numberOfTracks");
  SurveyMatrix.Write("SurveyMatrix");
  SurveyVector.Write("SurveyVector");  
  TotalMatrix.Write("TotalMatrix");
  TotalVector.Write("TotalVector");
  //Q.Write("Q");
  //R.Write("R");
  //IndependentVector.Write("IndependentVector");

  theFile->Close();
  
}



void StandaloneMerger::analyze(const Event & event, const EventSetup& eventSetup){}


void StandaloneMerger::surveyAdd() {
  
  SurveyMatrix.ResizeTo(theIndex.GetNrows()*NDOFAlign, theIndex.GetNrows()*NDOFAlign);
  SurveyVector.ResizeTo(theIndex.GetNrows()*NDOFAlign, 1);
  TotalMatrix.ResizeTo(theIndex.GetNrows()*NDOFAlign, theIndex.GetNrows()*NDOFAlign);
  TotalVector.ResizeTo(theIndex.GetNrows()*NDOFAlign, 1);


  
  TMatrixD SurveyDat(theIndex.GetNrows()*NDOFAlign, 1);

  TFile surveyFile(surveyFileName.c_str());
  TMatrixD *SurveyCov = (TMatrixD *)surveyFile.Get("SurveyCov");
  TMatrixD *SurveyData = (TMatrixD *)surveyFile.Get("SurveyData");
  TMatrixD *SurveyIndex = (TMatrixD *)surveyFile.Get("SurveyIndex");
  


  std::vector<long> realIndex;
  for(int counter = 0; counter < theIndex.GetNrows(); ++counter) {
    for(int counterC = 0; counterC < SurveyIndex->GetNrows(); ++counterC) {
      if(theIndex(counter,0) == (*SurveyIndex)(counterC,0)) {
        realIndex.push_back(counterC);
        break;
      }
    }
  }

  for(int counteri = 0; counteri < theIndex.GetNrows(); ++counteri) {
    for(int alignN = 0; alignN < NDOFAlign; ++alignN) {
      SurveyDat(counteri*NDOFAlign+alignN, 0) = (*SurveyData)(realIndex.at(counteri)*NDOFAlign+alignN, 0);
      for(int alignN2 = 0; alignN2 < NDOFAlign; ++alignN2) {
        for(int counterj = 0; counterj < theIndex.GetNrows(); ++counterj) {
          SurveyMatrix(counteri*NDOFAlign+alignN, counterj*NDOFAlign+alignN2) = (*SurveyCov)(realIndex.at(counteri)*NDOFAlign+alignN, realIndex.at(counterj)*NDOFAlign+alignN2);
        }
      }
    }
  }

  SurveyVector = SurveyMatrix*SurveyDat;
  
  for(int i = 0; i < SurveyMatrix.GetNrows(); ++i) {
    for(int j = 0; j < SurveyMatrix.GetNcols(); ++j) {
      SurveyMatrix(i,j) = SurveyMatrix(i,j);
    }
    SurveyVector(i,0) = SurveyVector(i,0);
  }


  TotalMatrix = SurveyMatrix + C;
  TotalVector = SurveyVector + b;



  surveyFile.Close();

}
  



DEFINE_FWK_MODULE(StandaloneMerger);

---