LILinResp.cxx

Go to the documentation of this file.

#define LILinResp_cxx

#include <LILinResp.h>
#include <TF1.h>

using namespace std;

std::vector<Float_t> LILinResp::fXrefPINLG1;
std::vector<Float_t> LILinResp::fXrefPINLG2;
std::vector<Float_t> LILinResp::fXrefPINLG3;
std::vector<Float_t> LILinResp::fXrefPINLG4;
std::vector<Float_t> LILinResp::fXrefPINLG5;
std::vector<Float_t> LILinResp::fXrefPINLG6;

std::vector<Float_t> LILinResp::fXrefPINHG1;
std::vector<Float_t> LILinResp::fXrefPINHG2;
std::vector<Float_t> LILinResp::fXrefPINHG3;
std::vector<Float_t> LILinResp::fXrefPINHG4;
std::vector<Float_t> LILinResp::fXrefPINHG5;
std::vector<Float_t> LILinResp::fXrefPINHG6;

CalVaLinearity* LILinResp::myCV = 0;


CVSID("$Id: LILinResp.cxx,v 1.4 2005/04/18 14:24:54 anatael Exp $");

//-------------------------------------------------------------

LILinResp::LILinResp() { 

  MSG("LILinResp",Msg::kDebug) << "Making a LILinResp..." << endl;

  entryHG = 0;
  entryLG = 0;

  addressKey  = -1;

  slopeHG       = -1;
  interceptHG   = -1;
  slopeHG_e     = -1;
  interceptHG_e = -1;

  slopeLG       = -1;
  interceptLG   = -1;
  slopeLG_e     = -1;
  interceptLG_e = -1;

  gain  = -1;
  drift = -1;
  led   = -1;

  //Zeroing all pointers:
  resPlotHG_g = 0;
  resPlotLG_g = 0;
}

//-------------------------------------------------------------

LILinResp::~LILinResp() {

  if(resPlotHG_g) delete resPlotHG_g; 
  if(resPlotLG_g) delete resPlotLG_g; 

  delete myCV;
  myCV = 0;

  MSG("LILinResp",Msg::kInfo) << "Killing a LILinResp" << endl;
}

//-------------------------------------------------------------

Bool_t LILinResp::GetLinearity() {

  Bool_t success = true;

  //
  // Create VA-nonLinearity object -> linearise PINs, which are FD electronis:
  //
  if(!myCV) myCV = new CalVaLinearity( RawChannelId(),
                                       24,
                                       16061.3,
                                         4.63,
                                       13683.9,
                                       6650 );

  //
  // Fill TGraph with ALL the information
  // for the linearity calculation:
  //
  vector<Float_t>::iterator xIter  = fX.begin();
  vector<Float_t>::iterator xeIter = fX_error.begin();

  vector<Float_t>::iterator xRefIterLG;// = NULL;
  vector<Float_t>::iterator xRefIterHG;// = NULL;
  
  //Use the LED number to ID the corresponding PIN (reference X):
  switch (led) {
  case 1:
    xRefIterLG = fXrefPINLG1.begin();
    xRefIterHG = fXrefPINHG1.begin();
    break;
  case 2:
    xRefIterLG = fXrefPINLG2.begin();
    xRefIterHG = fXrefPINHG2.begin();
    break;
  case 3:
    xRefIterLG = fXrefPINLG3.begin();
    xRefIterHG = fXrefPINHG3.begin();
    break;
  case 4:
    xRefIterLG = fXrefPINLG4.begin();
    xRefIterHG = fXrefPINHG4.begin();
    break;
  case 5:
    xRefIterLG = fXrefPINLG5.begin();
    xRefIterHG = fXrefPINHG5.begin();
    break;
  case 6:
    xRefIterLG = fXrefPINLG6.begin();
    xRefIterHG = fXrefPINHG6.begin();
    break;
  }

  //
  //LOOP over vectors:
  //
  while( xIter != fX.end() ) {

    //Avoid using entries with -1's
    if((*xRefIterHG) > 0  && (*xIter) > 2) {
      
      success = this->FillGraphHG( myCV->Linearize( (*xRefIterHG) ),
                                   0.,
                                   (*xIter),
                                   (*xeIter) );  
    }
    
    if((*xRefIterLG) > 0 && (*xIter) > 2) {
      
      success = this->FillGraphLG( myCV->Linearize( (*xRefIterLG) ),
                                   0., 
                                   (*xIter),
                                   (*xeIter) );  
    }

    xRefIterHG++;
    xRefIterLG++;
    xIter++;
    xeIter++;
  }    

  //
  //Get linearity:
  //

  TF1* lin = new TF1("lin","[0]+[1]*x"); //"pol1");


  lin->FixParameter(0,0);
  //lin->SetParLimits(0,0.001,-0.001);
  lin->SetParameter(1,0.3);


  Int_t out = -1;

  //This complication is needed to tune the fit to the lowest possible region, which is LED dependent
  if(typeFEE == ElecType::kVA) {
    switch( led ) {
    case 1:
      //if(resPlotHG_g) out = resPlotHG_g->Fit(lin,"QN","",0,4000);
      break;
    case 2:
      //if(resPlotHG_g) out = resPlotHG_g->Fit(lin,"QN","",0,2200);
      break;
    case 3:
      //if(resPlotHG_g) out = resPlotHG_g->Fit(lin,"QN","",0,5500);
      break;
    case 4:
      if(resPlotHG_g) out = resPlotHG_g->Fit(lin,"QN","",0,3800);
      break;
    case 5:
      if(resPlotHG_g) out = resPlotHG_g->Fit(lin,"QN","",0,5000);
      break;
    case 6:
      if(resPlotHG_g) out = resPlotHG_g->Fit(lin,"QN","",0,3000);
      break;
    }
  }
  else {
    switch( led ) {
    case 1:
      if(resPlotHG_g) out = resPlotHG_g->Fit(lin,"QN","",0,4000);
      break;
    case 2:
      if(resPlotHG_g) out = resPlotHG_g->Fit(lin,"QN","",0,2200);
      break;
    case 3:
      if(resPlotHG_g) out = resPlotHG_g->Fit(lin,"QN","",0,5500);
      break;
    case 4:
      //if(resPlotHG_g) out = resPlotHG_g->Fit(lin,"QN","",0,3500);
      break;
    case 5:
      //if(resPlotHG_g) out = resPlotHG_g->Fit(lin,"QN","",0,5000);
      break;
    case 6:
      //if(resPlotHG_g) out = resPlotHG_g->Fit(lin,"QN","",0,4000);
      break;
    }
  }

  slopeHG       = lin->GetParameter(1);
  slopeHG_e     = lin->GetParError(1);
  interceptHG   = lin->GetParameter(0);
  interceptHG_e = lin->GetParError(0);


  lin->FixParameter(0,0);
  lin->SetParameter(1,0.3);
  

  //This complication is needed to tune the fit to the lowest possible region, which is LED dependent
  if(typeFEE == ElecType::kVA) {
    switch( led ) {
    case 1:
      //if(resPlotLG_g) out = resPlotLG_g->Fit(lin,"QN","",0,500);
      break;
    case 2:
      //if(resPlotLG_g) out = resPlotLG_g->Fit(lin,"QN","",0,500);
      break;
    case 3:
      //if(resPlotLG_g) out = resPlotLG_g->Fit(lin,"QN","",0,500);
      break;
    case 4:
      if(resPlotLG_g) out = resPlotLG_g->Fit(lin,"QN","",0,800);
      break;
    case 5:
      if(resPlotLG_g) out = resPlotLG_g->Fit(lin,"QN","",0,1000);
      break;
    case 6:
      if(resPlotLG_g) out = resPlotLG_g->Fit(lin,"QN","",0,600);
      break;
    }
  }
  else {
    switch( led ) {
    case 1:
      if(resPlotLG_g) out = resPlotLG_g->Fit(lin,"QN","",0,800);//1000);
      break;
    case 2:
      if(resPlotLG_g) out = resPlotLG_g->Fit(lin,"QN","",0,800);//1100);
      break;
    case 3:
      if(resPlotLG_g) out = resPlotLG_g->Fit(lin,"QN","",0,800);//1000,2300); //Structure at low charges
      break;
    case 4:
      if(resPlotLG_g) out = resPlotLG_g->Fit(lin,"QN","",0,800);//900);
      break;
    case 5:
      if(resPlotLG_g) out = resPlotLG_g->Fit(lin,"QN","",0,800);//1200);
      break;
    case 6:
      if(resPlotLG_g) out = resPlotLG_g->Fit(lin,"QN","",0,800);//1200);
      break;
    }
  }


  //Check quality of LG, HG info is not used
  if(slopeLG==0) return false;
  if(slopeLG>0.29) return false;


  slopeLG       = lin->GetParameter(1);
  slopeLG_e     = lin->GetParError(1);
  interceptLG   = lin->GetParameter(0);
  interceptLG_e = lin->GetParError(0);

  delete lin;



  xIter  = fX.begin();
  //xeIter = fX_error.begin();

  //vector<Float_t>::iterator xRefIterLG;
  //vector<Float_t>::iterator xRefIterHG;
  
  //Use the LED number to ID the corresponding PIN (reference X):
  switch (led) {
  case 1:
    xRefIterLG = fXrefPINLG1.begin();
    xRefIterHG = fXrefPINHG1.begin();
    break;
  case 2:
    xRefIterLG = fXrefPINLG2.begin();
    xRefIterHG = fXrefPINHG2.begin();
    break;
  case 3:
    xRefIterLG = fXrefPINLG3.begin();
    xRefIterHG = fXrefPINHG3.begin();
    break;
  case 4:
    xRefIterLG = fXrefPINLG4.begin();
    xRefIterHG = fXrefPINHG4.begin();
    break;
  case 5:
    xRefIterLG = fXrefPINLG5.begin();
    xRefIterHG = fXrefPINHG5.begin();
    break;
  case 6:
    xRefIterLG = fXrefPINLG6.begin();
    xRefIterHG = fXrefPINHG6.begin();
    break;
  }

  //
  //LOOP over the responses of each LILinResp object:
  //
  while( xIter != fX.end() ) {
    
    //If there was bad data (ID by "-1") -> label as bad data.
    if( (*xRefIterLG)==-1 || (*xRefIterHG)==-1 || (*xIter)==-1 ) {
      fResidual.push_back( -1 );
      fResidualPMT.push_back( -1 );
    }
    else {

      //Furthermore, HG suffers from non-linearity of VA: BEWARE 
      //Float_t func = slopeHG * (*xRefIterHG) + interceptHG; 
      Float_t func = slopeLG * (*xRefIterLG) + interceptLG;
      
      //Calculating residuals:
      Float_t res    = (*xIter) - func;
      //% residual
      //Normalise by LIN:
      if( func != 0 ) res /= func;
      //Normalise by NON-LIN:
      //if( (*xIter)!= 0 ) res /= (*xIter);
      else res = -1;

      //Get residuals into vector of residuals:
      fResidual.push_back( res );


      //Remove linearity from PMT response on FD side:
      if(typeFEE == ElecType::kVA) {

        const Float_t linResp = myCV->Linearize( (*xIter)*gain )/gain;
        //Calculating residuals:
        Float_t resPMT = linResp - func;

        //Normalise by LIN:
        if( func != 0 ) resPMT /= func;
        //Normalise by Non-LIN:
        //if( linResp != 0 ) resPMT /= linResp;
        else resPMT = -1;
        fResidualPMT.push_back( resPMT );

      }      
      else 
        fResidualPMT.push_back( -1 ); //For ND.

    }
    
    xIter++;
    xRefIterHG++;
    xRefIterLG++;
  }

  return success;
}

//--------------------------------------------------------

Bool_t LILinResp::FillGraphHG( const Float_t x,
                             const Float_t xe,
                             const Float_t y,
                             const Float_t ye ) {

  //Check whether created or not:
  if(!resPlotHG_g) {

    MSG("LILinResp",Msg::kVerbose) << "Getting TGraph..." << endl;

    resPlotHG_g = new TGraphErrors();

    //Setting Style of Graph:
    resPlotHG_g->SetLineWidth(4);
    resPlotHG_g->SetLineColor(1);
    resPlotHG_g->SetMarkerSize(1);
    resPlotHG_g->SetMarkerStyle(20);
    resPlotHG_g->SetMarkerColor(2);
  }

  if(resPlotHG_g) {

    resPlotHG_g->SetPoint(entryHG,x,y);
    resPlotHG_g->SetPointError(entryHG,xe,ye);

    entryHG++;
    return true;
  }

  else return false;
}

//--------------------------------------------------------

Bool_t LILinResp::FillGraphLG( const Float_t x,
                             const Float_t xe,
                             const Float_t y,
                             const Float_t ye ) {

  //Check whether created or not:
  if(!resPlotLG_g) {

    MSG("LILinResp",Msg::kVerbose) << "Getting TGraph..." << endl;

    resPlotLG_g = new TGraphErrors();

    //Setting Style of Graph:
    resPlotLG_g->SetLineWidth(4);
    resPlotLG_g->SetLineColor(1);
    resPlotLG_g->SetMarkerSize(1);
    resPlotLG_g->SetMarkerStyle(20);
    resPlotLG_g->SetMarkerColor(2);
  }

  if(resPlotLG_g) {

    resPlotLG_g->SetPoint(entryLG,x,y);
    resPlotLG_g->SetPointError(entryLG,xe,ye);

    entryLG++;
    //resPlotLG_g->Print();
    return true;
  }

  else return false;
}

//--------------------------------------------------------

void LILinResp::PrintMe() {

  //cout << gain << "\t" << drift << "\t" << led << "\t" << addressKey << "\t" << entryHG << "\t" << entryLG << "\t" << typeFEE << endl;

  //resPlotHG_g->Print();
  //resPlotLG_g->Print();

  //cout << "X:    " << fX.size() << " (" << (&fX) << ")" << endl;
  //cout << "X_e:  " << fX_error.size()  << " (" << (&fX_error)  << ")" << endl;
  //cout << "Xres: " << fResidual.size() << " (" << (&fResidual) << ")" << endl;

  //cout << LILinResp::fXrefPINHG1.size() << " (" << (&fXrefPINHG1) << ")" << endl;
  //cout << LILinResp::fXrefPINHG2.size() << " (" << (&fXrefPINHG2) << ")" << endl;
  //cout << LILinResp::fXrefPINHG3.size() << " (" << (&fXrefPINHG3) << ")" << endl;
  //cout << LILinResp::fXrefPINHG4.size() << " (" << (&fXrefPINHG4) << ")" << endl;
  //cout << LILinResp::fXrefPINHG5.size() << " (" << (&fXrefPINHG5) << ")" << endl;
  //cout << LILinResp::fXrefPINHG6.size() << " (" << (&fXrefPINHG6) << ")" << endl;

  //cout << LILinResp::fXrefPINLG1.size() << endl;
  //cout << LILinResp::fXrefPINLG2.size() << endl;
  //cout << LILinResp::fXrefPINLG3.size() << endl;
  //cout << LILinResp::fXrefPINLG4.size() << endl;
  //cout << LILinResp::fXrefPINLG5.size() << endl;
  //cout << LILinResp::fXrefPINLG6.size() << endl;

  //cout << "SlopeHG: " << slopeHG << "\t" << interceptHG << endl;
  //cout << "SlopeHG_e: " << slopeHG_e << "\t" << interceptHG_e << endl;

  //cout << "SlopeLG: " << slopeLG << "\t" << interceptLG << endl;
  //cout << "SlopeLG_e: " << slopeLG_e << "\t" << interceptLG_e << endl << endl;


  return;

}

//EOF

---