AttenCorr.cxx

Go to the documentation of this file.

// 
// October 9 2007
//
// Author: Pedro O. and Mark T. 
// 
// Calculates attenuation corrections for strips

#include <iostream>
#include <cmath>
#include "TMath.h"
#include "TClonesArray.h"
#include "AttenCorr.h"
#include "StandardNtuple/NtpStRecord.h"
#include "CandNtupleSR/NtpSRStrip.h"
#include "UgliGeometry/UgliGeomHandle.h"

using std::cout;
using std::endl;

static const float  MEAN_ATTEN = 0.326; //mean attenuation correction of near and far dets
static const float  ND_CORRECTION = 1.;
static const float  REFLECTIVITY = 0.8;

AttenCorr::AttenCorr()
{  
  Reset();
}

AttenCorr::~AttenCorr()
{
}

// ***************** METHODS ***********************

void AttenCorr::Reset(){

  qpeE=-100;
  qpeW=-100;
  lE=-1;
  lW=-1;
  halflength=-1;
  wlspigtail_w=-1;
  wlspigtail_e=-1;
  clearfiber_w=-1;
  clearfiber_e=-1;
  centeru=-1;
  centerv=-1;
  qpeWcor=-100;
  qpeEcor=-100;;
  attenE=-1;
  attenW=-1;

}

// -----------------------------------------------------
void AttenCorr::Setup(const NtpStRecord* record,const NtpSRStrip* strip,unsigned short qtype){

  Reset();

  //Geometry handles and such
  UgliGeomHandle ugh(record->GetHeader().GetVldContext()); 
  detector = record->GetHeader().GetVldContext().GetDetector();  
  PlexStripEndId pseid(record->GetHeader().GetVldContext().GetDetector()
                       ,strip->plane, strip->strip);
  UgliStripHandle stripgeom(ugh.GetStripHandle(pseid));

  //Necessary quantities
  if (qtype == 0) {
    qpeE=strip->ph0.pe;
    qpeW=strip->ph1.pe;
  }
  if (qtype == 5) {
    qpeE = strip->ph0.sigcor;
    qpeW = strip->ph1.sigcor;
  }
  halflength = stripgeom.GetHalfLength();
  wlspigtail_w = stripgeom.WlsPigtail(StripEnd::kWest);
  wlspigtail_e = stripgeom.WlsPigtail(StripEnd::kEast);
  clearfiber_w = stripgeom.ClearFiber(StripEnd::kWest);
  clearfiber_e = stripgeom.ClearFiber(StripEnd::kEast);
  const TVector3 globalpos(stripgeom.GlobalPos(0));
  centerv = .7071067812*(globalpos.Y()-globalpos.X());
  centeru = .7071067812*(globalpos.Y()+globalpos.X());

  /*cout << "The wlspigtails " << endl;//<--tmp
  cout << wlspigtail_w << "," << wlspigtail_e << endl;
  cout << "The clearfibers " << endl;
  cout << clearfiber_w << "," << clearfiber_e << endl;
  cout << "The centeru and centerv " << endl;
  cout << centeru << "," << centerv << endl;*/

}

// --------------------------------------------------
void AttenCorr::CalcPECorr(const int){

  float lwlsE = 0;
  float lwlsW = 0;
  float fE = 0;
  float fW = 0;

  if(detector==0x02 || detector==0x04){//<-- Far Detector OR Caldet

    lwlsE = wlspigtail_e + lE;
    lwlsW = wlspigtail_w + lW;
    fE = 1./( 0.666*exp(-lwlsE/7.05)+ 0.333*exp(-lwlsE/1.05));
    fW = 1./( 0.666*exp(-lwlsW/7.05)+ 0.333*exp(-lwlsW/1.05));
    fE = fE*exp(+(clearfiber_e/10.));
    fW = fW*exp(+(clearfiber_w/10.));
    qpeEcor = qpeE*fE*MEAN_ATTEN;
    qpeWcor = qpeW*fW*MEAN_ATTEN;
    qpecor = (qpeEcor+qpeWcor)*ND_CORRECTION;  

    attenE=fE;
    attenW=fW;
    
  } else if(detector==0x01){//<--Near Detector (only W readout; reflection in E)

    lwlsW = wlspigtail_w + lW;
    lwlsE = 4*halflength - fabs(lW) + wlspigtail_w;
      
    fW = 1.*(0.666*exp(-lwlsW/7.05)+ 0.333*exp(-lwlsW/1.05));//Caution: it's not 1/fW as in FarDet; same for fE below;

    //A la Pedro
    fE = 1.*(0.666*exp(-lwlsE/7.05)+ 0.333*exp(-lwlsE/1.05))*REFLECTIVITY;
    float fCor = 2./((fE+fW)*exp(-clearfiber_w/10.));
    fCor=fCor*1.587;//determined from muons, on Nov-23-2007, for cedar_phy
    
    //A la Caius
    //float fCor = 1./(fW*exp(-clearfiber_w/10.));//<--tmp
    //fCor=fCor*1.9454;

    qpeWcor = qpeW*fCor*MEAN_ATTEN;
    qpecor = qpeWcor*ND_CORRECTION;
       
    attenE=0;
    attenW=fCor;
    
  }


}


//---------------------------------------------------
void AttenCorr::CalcCorr(const NtpStRecord* record,const NtpSRStrip* strip,Float_t meanUTPos, Float_t meanVTPos, unsigned short qtype){

  Setup(record,strip,qtype);

  if(strip->planeview==2){ //U view
    //location in V of E end:
    lE = centerv + halflength - meanVTPos; 
    lW = meanVTPos - (centerv - halflength);
  }else{  //V view
    lE = meanUTPos - (centeru - halflength);
    lW = centeru+halflength -meanUTPos;
  }

  //cout << "The lE, lW and halflength " << endl;//<--tmp
  //cout << lE << "," << lW << "," << halflength << endl;

  CalcPECorr(detector);

  //MHO0509
  //cout << "MHO - AttenCorr found attenE*MEAN_ATTEN " << attenE*MEAN_ATTEN << " and attenW*MEAN_ATTEN " << attenW*MEAN_ATTEN << " in strip " << strip->strip << " and plane " << strip->plane << endl;
  
}

---