BinCorrelationEM.cxx

Go to the documentation of this file.

#include "BinCorrelationEM.h"
#include "MessageService/MsgService.h"
#include <fstream>

using namespace std;

CVSID("$Id: BinCorrelationEM.cxx,v 1.2 2005/10/21 11:22:30 cbs Exp $");

BinCorrelationEM::BinCorrelationEM(Double_t *inputs){

  fErrArray = NULL;
  fCovMat = NULL;
  
  if(inputs){
    fInputParams = new Double_t[6];
    for(Int_t i=0;i<6;i++){
      fInputParams[i] = inputs[i];
    }
    Init();
  }
  else fInputParams = NULL;

}

BinCorrelationEM::~BinCorrelationEM(){
  if(fErrArray) delete [] fErrArray;
  if(fCovMat) delete [] fCovMat;
  if(fInputParams) delete [] fInputParams;
}

void BinCorrelationEM::ReInit(Double_t *inputs){

  //delete old inputs
  Double_t CurrentEnergy = -1;
  if(fInputParams) {
    CurrentEnergy = fInputParams[0];
    delete [] fInputParams;
  }
  
  //read in new ones
  fInputParams = new Double_t[6];
  for(Int_t i=0;i<6;i++){
    fInputParams[i] = inputs[i];
  }
  
  //delete fErrArray whatever happens
  //only delete fCovMat if the energy has changed
  if(fInputParams[0]>0){
    Init();
    delete [] fErrArray;
    fErrArray = NULL;
    if(CurrentEnergy!=fInputParams[0]) {
      delete [] fCovMat;
      fCovMat = NULL;
    }
  }
  else MSG("FitShowerEM",Msg::kWarning)
    <<"Trying to ReInit with energy <=0 ... Probably everything will break..." 
    << endl;
}

void BinCorrelationEM::Init(){

}

Double_t BinCorrelationEM::CalcCorrelation(Double_t /*t*/, Double_t /*r*/){
  
  return 1;

}

Double_t *BinCorrelationEM::GetErrArray(Int_t LowBin_pl,Int_t HighBin_pl,
                                        Int_t LowBin_st,Int_t HighBin_st){

  if(fErrArray) return fErrArray;
  
  int nplane = int(HighBin_pl-LowBin_pl+1);
  int nstrip = int(HighBin_st-LowBin_st+1);
  int totstp = nplane*nstrip;
  
  fErrArray = new double[totstp*totstp];


  if(!fCovMat){ //if it doesn't exist, make it and fill it    
    MSG("FitShowerEM",Msg::kDebug) << "Setting up Correlation Matrix" << endl;
    
    fCovMat = new double[720*720];
    ifstream infile;
    char name[256];
    sprintf(name,"BinErrors/CovMat_%iGeV_NoXTalk.dat",int(fInputParams[0]));
    infile.open(name);
    if(!infile.is_open()) {
      //open default
      sprintf(name,"BinErrors/CovMat_3GeV_NoXTalk.dat");
      infile.open(name);
    }
    
    if(!infile.is_open()) {
      MSG("FitShowerEM",Msg::kDebug) 
        << "No txt file for correlations found - using identity matrix" << endl;
      for(int i=0;i<720;i++){
        for(int j=0;j<720;j++) {
          if(i==j) fCovMat[i*720+j] = 1;
          else fCovMat[i*720+j] = 0;
        }
      }
    }
    else {
      int a,b;
      double c;
      while(infile>>a>>b>>c) {
        if(c>1e6) fCovMat[a*720+b]=c;
        else fCovMat[a*720+b]=0;
      }
      infile.close();    
      MSG("FitShowerEM",Msg::kDebug) << "Read in file " << name << endl;
    }
  }
  

  Double_t SimVtx[3] = {11.5,11.5,0.5};
  //fill fErrArray
  for(int i=0;i<totstp;i++){
    int plane_i = int(LowBin_pl) + i/nstrip;
    int strip_i = int(LowBin_st) + i%nstrip;
    if(plane_i%2==0) {
      strip_i += int(-fInputParams[1] + SimVtx[0]);      
    }
    else {
      strip_i += int(-fInputParams[2] + SimVtx[1]);
    }
    plane_i += int(-fInputParams[3] + SimVtx[2]);
    for(int j=0;j<totstp;j++){
      int plane_j = int(LowBin_pl) + j/nstrip;
      int strip_j = int(LowBin_st) + j%nstrip;      
      if(plane_j%2==0) {
        strip_j = int(-fInputParams[1] + SimVtx[0]) + strip_j;
      }
      else {
        strip_j = int(-fInputParams[2] + SimVtx[1]) + strip_j;
      }
      plane_j += int(-fInputParams[3] + SimVtx[2]);
            
      if(plane_i>=0&&plane_i<30&&plane_j>=0&&plane_j<30
         &&strip_i>=0&&strip_i<24&&strip_j>=0&&strip_j<24){
        fErrArray[i*totstp+j] = fCovMat[(plane_i*24+strip_i)*720 + 
                                        plane_j*24+strip_j];
      }
      else { //out of defined region
        if(i==j) fErrArray[i*totstp+j] = 1.;  //diagonal terms
        else fErrArray[i*totstp+j] = 0;       //off-diagonal
      }

      MSG("FitShowerEM",Msg::kVerbose) << "Plane_i: " << plane_i 
                                       << " Strip_i: " << strip_i 
                                       << "\nPlane_j: " << plane_j 
                                       << " Strip_j: " << strip_j
                                       << "\nErrArray_ij: " 
                                       << fErrArray[i*totstp+j]
                                       << endl;
    }
  }

  return fErrArray;

}

---