LIPmt.cxx

Go to the documentation of this file.

// Program name: LIPmt.cxx                
//
// Package: LISummary          
//
// Coded by Jeff Hartnell Nov/2003            
//
// Purpose: 
//
// Contact: jeffrey.hartnell@physics.ox.ac.uk                         

#include "MessageService/MsgService.h"

#include "LISummary/LIPmt.h"

#include <cmath>

ClassImp(LIPmt)
  
CVSID("$Id: LIPmt.cxx,v 1.7 2007/01/15 20:12:15 rhatcher Exp $");

//......................................................................

LIPmt::LIPmt()
{
  MSG("LIPmt", Msg::kDebug) 
    <<"Running LIPmt default constructor..."<<endl;
  
  fInitialised=false;

  MSG("LIPmt", Msg::kDebug) 
    <<"Finished LIPmt default constructor"<<endl;
}

//......................................................................

LIPmt::LIPmt(Int_t numPixels,Int_t numPixelSpots)
{
  MSG("LIPmt", Msg::kDebug) 
    <<"Running LIPmt constructor..."<<endl;
  
  this->InitialiseVariables(numPixels,numPixelSpots);  
  fInitialised=true;

  MSG("LIPmt", Msg::kDebug) 
    <<"Finished LIPmt constructor"<<endl;
}

//......................................................................

LIPmt::~LIPmt()
{
  MSG("LIPmt", Msg::kDebug) 
    <<"Running LIPmt destructor..."<<endl;
  
  
  MSG("LIPmt", Msg::kDebug) 
    <<"Finished LIPmt destructor"<<endl;
}

//......................................................................

bool operator==(LIPmt a, LIPmt b)
{
  return a.GetGains()==b.GetGains();
}

//......................................................................

bool operator<(LIPmt a, LIPmt b)
{
  return a.GetGains().size()<b.GetGains().size();
}

//......................................................................

void LIPmt::Initialise(Int_t numPixels,Int_t numPixelSpots)
{
  MSG("LIPmt",Msg::kVerbose) 
    <<"Running Initialise method..."<<endl;

  if (!fInitialised){
    this->InitialiseVariables(numPixels,numPixelSpots);
    fInitialised=true;
  }

  MSG("LIPmt",Msg::kVerbose) 
    <<"Initialise method finished"<<endl;
}

//......................................................................

void LIPmt::InitialiseVariables(Int_t numPixels,Int_t numPixelSpots)
{
  MSG("LIPmt",Msg::kVerbose) 
    <<"Running InitialiseVariables method..."<<endl;

  MSG("LIPmt",Msg::kDebug)
    <<"Initialised variables, fInitialised="<<fInitialised
    <<", numPixels="<<numPixels<<", numPixelSpots="<<numPixelSpots 
    <<endl;

  fNumPixels=numPixels;
  fNumPixelSpots=numPixelSpots;
  fGainsSum=0;
  fGainsNum=0;
  fS="";

  //initialise all the vectors
  for (Int_t spot=0;spot<fNumPixels*fNumPixelSpots;spot++){
    fGains.push_back(0);
    fCounter.push_back(0);
  }

  if (numPixelSpots>1){
    //initialise the vector to hold the average gain for a pixel
    for (Int_t pix=0;pix<fNumPixels;pix++){
      fPixelGainsSum.push_back(0);
      fPixelGainsNum.push_back(0);
    }
  }

  MSG("LIPmt",Msg::kVerbose) 
    <<"InitialiseVariables method finished"<<endl;
}

//......................................................................

Bool_t LIPmt::IsEmpty() const
{
  if (fGainsNum>0) return false;
  else return true;
}

//......................................................................

void LIPmt::Print(Option_t* /* option */) const
{
  if (!fInitialised) return;

  //check that you know the pmt type
  if (fNumPixels<=0){
    MSG("LIPmt",Msg::kWarning)
      <<"Num of pixels <=0, fNumPixels="<<fNumPixels
      <<", fNumPixelSpots="<<fNumPixelSpots<<endl
      <<"The object must be initialised"<<endl;
    return;
  }

  Int_t pmtWidth=static_cast<Int_t>
    (sqrt(static_cast<Float_t>(fNumPixels)));

  for(int row=0;row<pmtWidth;row++){
    for(int p=0;p<pmtWidth;p++){
      
      MSG("LIPmt",Msg::kInfo)
        <<static_cast<Int_t>(this->GetPixelGain(row*pmtWidth+p))<<" ";
    }
    //leave a line between rows
    MSG("LIPmt",Msg::kInfo)<<endl;
  }
}

//......................................................................

TH2F* LIPmt::GetPmtFaceMap(string sTitle) const
{
  if (!fInitialised) return 0;

  //check that you know the pmt type
  if (fNumPixels<=0){
    MSG("LIPmt",Msg::kWarning)
      <<"Num of pixels <=0, fNumPixels="<<fNumPixels
      <<", fNumPixelSpots="<<fNumPixelSpots<<endl
      <<"The object must be initialised"<<endl;
    return 0;
  }

  Int_t pmtWidth=static_cast<Int_t>
    (sqrt(static_cast<Float_t>(fNumPixels)));

  TH2F* hPmtFace=new TH2F(sTitle.c_str(),sTitle.c_str(),
                          pmtWidth+2,-1,pmtWidth+1,
                          pmtWidth+2,-1,pmtWidth+1);
  hPmtFace->SetTitle(sTitle.c_str());
  hPmtFace->SetFillColor(2);  
  //hPmtFace->SetBit(TH1::kCanRebin);

  for(int row=0;row<pmtWidth;row++){
    for(int p=0;p<pmtWidth;p++){

      hPmtFace->Fill(p,row,this->GetPixelGain(row*pmtWidth+p));

    }
  }
  
  return hPmtFace;
}

//......................................................................

void LIPmt::SetPoint(Int_t pixel,Int_t pixelSpot,Double_t gain)
{
  MSG("LIPmt",Msg::kVerbose)<<"Running SetPoint method..."<<endl;

  //make pixelSpot start at zero
  pixelSpot--;

  MSG("LIPmt",Msg::kDebug)
    <<"("<<pixel<<","<<pixelSpot<<")"
    <<" Setting point: gain="<<gain<<endl;

  if (pixel>=0 && pixel<fNumPixels){
    if (pixelSpot>=0 && pixelSpot<fNumPixelSpots){

      Int_t spotIndex=pixel*fNumPixelSpots+pixelSpot;
      
      //work out the running total
      fGainsSum+=gain;
      fGainsNum++;

      //store the individual gain
      fGains[spotIndex]=gain;

      //sanity check
      if (fCounter[spotIndex]!=0){
        MSG("LIPmt",Msg::kWarning) 
          <<"This pixel and pixelSpot have already been set!"<<endl
          <<"Pixel="<<pixel<<", pixelSpot="<<pixelSpot
          <<", counter="<<fCounter[spotIndex]<<endl;
      }
      fCounter[spotIndex]++;

      //calculate and store the running total of gains on a pixel spot
      if (fNumPixelSpots>1){
        fPixelGainsSum[pixel]+=gain;
        fPixelGainsNum[pixel]++;
      }

    }
    else{
      MSG("LIPmt",Msg::kWarning)
        <<"Pixel Spot is wrong, pixel="<<pixel
        <<", pixelSpot="<<pixelSpot<<", gain="<<gain<<endl;
    }
  }
  else{
    MSG("LIPmt",Msg::kWarning)
      <<"Pixel is wrong, pixel="<<pixel
      <<", pixelSpot="<<pixelSpot<<", gain="<<gain<<endl;
  }
  
  MSG("LIPmt",Msg::kVerbose)<<"SetPoint method finished"<<endl;
}

//......................................................................

void LIPmt::AddMultiPoint(Int_t pixel,Int_t pixelSpot,Double_t gain)
{

  MSG("LIPmt",Msg::kVerbose)<<"Running AddMultiPoint method..."<<endl;

  //make pixelSpot start at zero
  pixelSpot--;

  MSG("LIPmt",Msg::kDebug)
    <<"("<<pixel<<","<<pixelSpot<<")"
    <<" Adding new point: gain="<<gain<<endl;

  if (pixel>=0 && pixel<fNumPixels){
    if (pixelSpot>=0 && pixelSpot<fNumPixelSpots){

      Int_t spotIndex=pixel*fNumPixelSpots+pixelSpot;
      
      //work out the running total
      fGainsSum+=gain;
      fGainsNum++;

      //sum the individual gain
      fGains[spotIndex]+=gain;
      fCounter[spotIndex]++;

      //calculate and store the running total of gains on a pixel
      if (fNumPixelSpots>1){
        fPixelGainsSum[pixel]+=gain;
        fPixelGainsNum[pixel]++;
      }

    }
    else{
      MSG("LIPmt",Msg::kWarning)
        <<"Pixel Spot is wrong, pixel="<<pixel
        <<", pixelSpot="<<pixelSpot<<", gain="<<gain<<endl;
    }
  }
  else{
    MSG("LIPmt",Msg::kWarning)
      <<"Pixel is wrong, pixel="<<pixel
      <<", pixelSpot="<<pixelSpot<<", gain="<<gain<<endl;
  }
  
  MSG("LIPmt",Msg::kVerbose)<<"AddMultiPoint method finished"<<endl;
}

//......................................................................

Double_t LIPmt::GetGain(Int_t pixel,Int_t pixelSpot) const
{
  //make pixelSpot start at zero
  pixelSpot--;

  if (pixel>=0 && pixel<fNumPixels){
    if (pixelSpot>=0 && pixelSpot<fNumPixelSpots){
      Int_t spotIndex=pixel*fNumPixelSpots+pixelSpot;

      //protect against fpe
      if (fCounter[spotIndex]>0){
        return fGains[spotIndex]/fCounter[spotIndex];
      }
    }
  }
  return -1;
} 

//......................................................................

const vector<Double_t>& LIPmt::GetGains() const
{
  MSG("LIPmt",Msg::kVerbose) 
    <<"Running GetGain method..."<<endl;

  MSG("LIPmt",Msg::kVerbose) 
    <<"GetGain method finished"<<endl;
  return fGains;
}

//......................................................................

Double_t LIPmt::GetAvPmtGain() const
{
  MSG("LIPmt",Msg::kVerbose) 
    <<"Running GetAvPmtGain method..."<<endl;

  Double_t avGain=-1;
  if (fGainsNum>0) avGain=fGainsSum/fGainsNum;

  MSG("LIPmt",Msg::kVerbose) 
    <<"GetAvPmtGain method finished"<<endl;
  return avGain;
}

//......................................................................

Double_t LIPmt::GetPixelGain(Int_t pixel) const
{
  MSG("LIPmt",Msg::kVerbose) 
    <<"Running GetPixelGain method..."<<endl;

  Double_t pixelGain=-1;
  
  if (pixel>=0 && pixel<fNumPixels){
    
    //calculate the running total of gains on a pixel spot
    if (fNumPixelSpots>1){
      //protect against fpe
      if (fPixelGainsNum[pixel]>0){
        pixelGain=fPixelGainsSum[pixel]/fPixelGainsNum[pixel];
      }
    }
    else {//only one spot per pixel , e.g. M64
      //there is no average gain and only the gain of that pixel
      if (fCounter[pixel]>0) pixelGain=fGains[pixel]/fCounter[pixel];
    }
  }
  else {
    MSG("LIPmt",Msg::kWarning)<<"Pixel is wrong, pixel="<<pixel<<endl;
  }

  MSG("LIPmt",Msg::kVerbose) 
    <<"GetPixelGain method finished"<<endl;
  return pixelGain;
}

//......................................................................

Int_t LIPmt::GetPixelNum(Int_t pixel) const
{
  MSG("LIPmt",Msg::kVerbose) 
    <<"Running GetPixelNum method..."<<endl;

  Int_t pixelNum=0;

  if (pixel>=0 && pixel<fNumPixels){
    
    //calculate the running total of gains on a pixel spot
    if (fNumPixelSpots>1){
      pixelNum=fPixelGainsNum[pixel];
    }
    else {//only one spot per pixel
      //there is no average gain and only the gain of that pixel
      pixelNum=fCounter[pixel];
    }
  }
  else {
    MSG("LIPmt",Msg::kWarning)<<"Pixel is wrong, pixel="<<pixel<<endl;
  }

  MSG("LIPmt",Msg::kVerbose) 
    <<"GetPixelNum method finished"<<endl;
  return pixelNum;
}

//......................................................................

vector<Int_t> LIPmt::GetNearestNeighbours(Int_t pixel) const
{
  MSG("LIPmt",Msg::kVerbose)
    <<"Running GetNearestNeighbours method..."<<endl;

  vector<Int_t> nearestNeighbours;

  //check that you know the pmt type
  if (fNumPixels<=0){
    MSG("LIPmt",Msg::kWarning)
      <<"Num of pixels <=0, fNumPixels="<<fNumPixels
      <<", fNumPixelSpots="<<fNumPixelSpots<<endl
      <<"The object must be initialised, returning empty vector"<<endl;
    return nearestNeighbours;
  }

  if (pixel>=0 && pixel<fNumPixels){
    
    Int_t pmtWidth=static_cast
      <Int_t>(sqrt(static_cast<Float_t>(fNumPixels)));

    //first do the corners
    if (pixel==0 || pixel==pmtWidth-1 || 
        pixel==fNumPixels-pmtWidth || pixel==fNumPixels-1){
      //calculate each of the 4 corners separately
      if (pixel==0){
        nearestNeighbours.push_back(pixel+1);
        nearestNeighbours.push_back(pixel+pmtWidth);
      }
      else if (pixel==pmtWidth-1){
        nearestNeighbours.push_back(pixel-1);
        nearestNeighbours.push_back(pixel+pmtWidth);
      }
      else if (pixel==fNumPixels-pmtWidth){
        nearestNeighbours.push_back(pixel-pmtWidth);
        nearestNeighbours.push_back(pixel+1);
      }
      else if (pixel==fNumPixels-1){
        nearestNeighbours.push_back(pixel-pmtWidth);
        nearestNeighbours.push_back(pixel-1);
      }
    }
    //find the edges (corners are already found though)
    else if (pixel<pmtWidth || pixel>=fNumPixels-pmtWidth ||
             pixel%pmtWidth==0 || pixel%pmtWidth==pmtWidth-1){
      //bottom edge
      if (pixel<pmtWidth){
        nearestNeighbours.push_back(pixel-1);
        nearestNeighbours.push_back(pixel+1);
        nearestNeighbours.push_back(pixel+pmtWidth);
      }
      //left edge
      else if (pixel%pmtWidth==0){
        nearestNeighbours.push_back(pixel-pmtWidth);
        nearestNeighbours.push_back(pixel+1);
        nearestNeighbours.push_back(pixel+pmtWidth);
      }
      else if (pixel%pmtWidth==pmtWidth-1){
        nearestNeighbours.push_back(pixel-pmtWidth);
        nearestNeighbours.push_back(pixel-1);
        nearestNeighbours.push_back(pixel+pmtWidth);
      }
      else if (pixel>=fNumPixels-pmtWidth){
        nearestNeighbours.push_back(pixel-pmtWidth);
        nearestNeighbours.push_back(pixel-1);
        nearestNeighbours.push_back(pixel+1);
      }
    }
    //in the middle somewhere
    else {
      nearestNeighbours.push_back(pixel-pmtWidth);
      nearestNeighbours.push_back(pixel-1);
      nearestNeighbours.push_back(pixel+1);
      nearestNeighbours.push_back(pixel+pmtWidth);
    }
  }
  else {
    MSG("LIPmt",Msg::kWarning)<<"Pixel is wrong, pixel="<<pixel<<endl;
  }

  MSG("LIPmt",Msg::kVerbose) 
    <<"Finished GetNearestNeighbours method"<<endl;
  return nearestNeighbours;
}

//......................................................................

Double_t LIPmt::GetNearestNeighboursAvGain(Int_t pixel) const
{
  MSG("LIPmt",Msg::kVerbose) 
    <<"Running GetNearestNeighboursAvGain method..."<<endl;

  Double_t pixelGain=0;
  Int_t counter=0;
  
  if (pixel>=0 && pixel<fNumPixels){

    vector<Int_t> nearestNeighbours=this->GetNearestNeighbours(pixel);
    vector<Int_t>::iterator nnPixel=nearestNeighbours.begin();
    
    //loop over the nearestNeighbour (nn) pixels
    while (nnPixel!=nearestNeighbours.end()){

      //calculate the running total of gains on a pixel spot
      if (fNumPixelSpots>1){
        //protect against fpe
        if (fPixelGainsNum[*nnPixel]>0){
          pixelGain+=fPixelGainsSum[*nnPixel]/fPixelGainsNum[*nnPixel];
          counter++;
        }
      }
      //only one spot per pixel , e.g. M64
      else {
        //there is no average of spots so only the gain of that pixel
        if (fCounter[*nnPixel]>0){
          pixelGain+=fGains[*nnPixel]/fCounter[*nnPixel];
          counter++;
        }
      }
      nnPixel++;
    }

    //calc average
    if (counter>0) pixelGain/=counter;
    else pixelGain=-1;
  }
  else {
    pixelGain=-1;
    MSG("LIPmt",Msg::kWarning)<<"Pixel is wrong, pixel="<<pixel<<endl;
  }

  MSG("LIPmt",Msg::kVerbose) 
    <<"Finished GetNearestNeighboursAvGain method"<<endl;
  return pixelGain;
}

//......................................................................

Int_t LIPmt::GetNearestNeighboursNum(Int_t pixel) const
{
  MSG("LIPmt",Msg::kVerbose) 
    <<"Running GetNearestNeighboursNum method..."<<endl;

  Int_t totalNum=0;
  Int_t counter=0;

  if (pixel>=0 && pixel<fNumPixels){

    vector<Int_t> nearestNeighbours=this->GetNearestNeighbours(pixel);
    vector<Int_t>::iterator nnPixel=nearestNeighbours.begin();
    
    //loop over the nearestNeighbour (nn) pixels
    while (nnPixel!=nearestNeighbours.end()){

      //calculate the running total of gains on a pixel spot
      if (fNumPixelSpots>1){
        //protect against fpe
        if (fPixelGainsNum[*nnPixel]>0){
          totalNum+=fPixelGainsNum[*nnPixel];
          counter++;
        }
      }
      //only one spot per pixel , e.g. M64
      else {
        //there is no average of spots so only use pixel indexed vector
        if (fCounter[*nnPixel]>0){
          totalNum+=fCounter[*nnPixel];
          counter++;
        }
      }
      nnPixel++;
    }

    //calc average
    if (counter<=0) totalNum=-1;
  }
  else {
    totalNum=-1;
    MSG("LIPmt",Msg::kWarning)<<"Pixel is wrong, pixel="<<pixel<<endl;
  }

  MSG("LIPmt",Msg::kVerbose) 
    <<"Finished GetNearestNeighboursNum method"<<endl;
  return totalNum;
}

//......................................................................

void LIPmt::PrintNearestNeighbours() const
{
  MSG("LIPmt",Msg::kInfo)
    <<"Running PrintNearestNeighbours method..."<<endl;

  //check that you know the pmt type
  if (fNumPixels<=0){
    MSG("LIPmt",Msg::kWarning)
      <<"Num of pixels <=0, fNumPixels="<<fNumPixels
      <<", fNumPixelSpots="<<fNumPixelSpots<<endl
      <<"The object must be initialised"<<endl;
  }

  for (Int_t i=0;i<fNumPixels;i++){
    vector<Int_t> nearestNeighbours=this->GetNearestNeighbours(i);
    
    string sNn=lookup.GetVectorAsString(nearestNeighbours);
    
    MSG("LIPmt",Msg::kInfo) 
      <<"Pixel="<<i<<" has nearest neighbours: "<<sNn;
  }

  MSG("LIPmt",Msg::kInfo) 
    <<"Finished PrintNearestNeighbours method"<<endl;
}

//......................................................................

void LIPmt::FillGainsHisto(TH1F* hGainAvDet) const
{
  MSG("LIPmt",Msg::kDebug)<<"Running FillHisto method..."<<endl;

  //check histo and gains vector exists
  if (hGainAvDet && fGains.size()>0){

    vector<Double_t>::const_iterator gain=fGains.begin();
    vector<Int_t>::const_iterator counter=fCounter.begin();

    while(gain!=fGains.end()){
      //protect against fpe
      if (*counter>0){
        //fill histo
        hGainAvDet->Fill(*gain/(*counter));
      }
      gain++;
      counter++;
    }
  }
  
  MSG("LIPmt",Msg::kDebug)<<"Finished FillHisto method"<<endl;
}

//......................................................................

---