PhotonCalibratedBlueComputer.cxx

Go to the documentation of this file.

#include "PhotonCalibratedBlueComputer.h"
#include "PhotonTransportMaker.h"
#include "PhotonUtil.h"
#include "MessageService/MsgService.h"
#include "UgliGeometry/UgliGeomHandle.h"
#include "Calibrator/Calibrator.h"

CVSID( " $Id: PhotonCalibratedBlueComputer.cxx,v 1.5 2006/12/01 20:12:11 rhatcher Exp $" );

ClassImp(PhotonCalibratedBlueComputer)
PhotonTransportMakerProxy<PhotonCalibratedBlueComputer>
        gPhotonCalibratedBlueComputerProxy("photonCalibratedBlueComputer"); 


// Constructor
PhotonCalibratedBlueComputer::PhotonCalibratedBlueComputer() 
{
  Configure(PhotonConfiguration());
}

// Destructor
PhotonCalibratedBlueComputer::~PhotonCalibratedBlueComputer()
{
 
}

// Configure the model.
void 
PhotonCalibratedBlueComputer::Configure( const Registry& r )
{
  MSG("Photon",Msg::kDebug) << "Configuring CalibratedBlueComputer." << std::endl;

  r.Get("OverallLightOutput",fOverallLightOutput);
  r.Get("BirksConstant",fBirksConstant);
  r.Get("ClearFibreAttenN1",fClearFibreAttenN1);
  r.Get("ClearFibreAttenN2",fClearFibreAttenN2);
  r.Get("ClearFibreAttenLength1",fClearFibreAttenLength1);
  r.Get("ClearFibreAttenLength2",fClearFibreAttenLength2);
 }

// Reset the model.
void 
PhotonCalibratedBlueComputer::Reset( const VldContext&  )
{
}


// Computes the number of photons we'll need to track for the 
// hit. See the PhotonCount class for the various ways this
// can be calculated.

Bool_t 
PhotonCalibratedBlueComputer::ComputePhotons( const DigiScintHit* inHit,
                                              Double_t inBluePrescale,
                                              Double_t inWlsPrescale,
                                              Double_t inGreenPrescale,
                                              PhotonCount& outCount )
{
  // Number of photons made:
  if(!inHit) return false;
  double dE = inHit->DE();
  double dx = inHit->DS();
  double dEdx = dE/dx;
  
  double lightYield = ( 1/100.  // photons/eV
                        * 1e9 // GeV/ev
                        * inBluePrescale
                        * inWlsPrescale
                        * inGreenPrescale
                        * .0431 // fudge factor to tune to other methods.
                        );
  

  PlexStripEndId pos = inHit->StripEndId(); pos.SetEnd(StripEnd::kPositive);
  PlexStripEndId neg = inHit->StripEndId(); neg.SetEnd(StripEnd::kNegative);
        
  // Energy->MIP
  // Apply normalization and birks' constant.
  double mblue =  
    fOverallLightOutput             // Tuning parameter.
    * dE/(1. + fBirksConstant*dEdx) // Birk's law
    * lightYield;                   // photons/GeV

  const Calibrator& cal = Calibrator::Instance();

  // Now we need to adjust only for the strip-to-strip variation.
  // To get this, we ask for MIP correction, strip correction, 
  // and attenuation correction, but corrected to zero green fibre length
  // and zero clear fibre length.
  
  UgliGeomHandle ugli(fContext);
  UgliStripHandle ustrip = ugli.GetStripHandle(inHit->StripEndId());
  float greenLen[2];
  greenLen[0] = ustrip.GetHalfLength() + ustrip.WlsPigtail(StripEnd::kNegative);
  greenLen[1] = ustrip.GetHalfLength() + ustrip.WlsPigtail(StripEnd::kPositive);

  float clearLen[2];
  clearLen[0] = ustrip.ClearFiber(StripEnd::kNegative);
  clearLen[1] = ustrip.ClearFiber(StripEnd::kPositive);

  
  // Correct for MIPs (overall detector scale)
  double sigmap[2];
  sigmap[0] = cal.DecalMIP( mblue * 0.5, neg); // Send half the light left
  sigmap[1] = cal.DecalMIP( mblue * 0.5, pos); // and half to the right

  // Get the scintillator light output.
  
  // The following does a 'no correction' correction: it
  // tells you now much to correct for a hit at the very end of the WLS fibre.
  float sigcorr[2];
  sigcorr[0] = cal.DecalAttenCorrected(sigmap[0],-greenLen[0],neg);
  sigcorr[1] = cal.DecalAttenCorrected(sigmap[1], greenLen[1],pos);


  // Correct for clear fibre attenuation.. i.e. what if ethere were none.
  sigcorr[0]= sigcorr[0] / PhotonUtil::DoubleExp( clearLen[0], 
                                                  fClearFibreAttenN1,      fClearFibreAttenN2,
                                                  fClearFibreAttenLength1, fClearFibreAttenLength2 );
  sigcorr[1]= sigcorr[1] / PhotonUtil::DoubleExp( clearLen[1], 
                                                  fClearFibreAttenN1,      fClearFibreAttenN2,
                                                  fClearFibreAttenLength1, fClearFibreAttenLength2 );
  
  // SigCorr->SigLin (siglin = ADCs, because no drift, no linearity)
  double siglin[2];
  siglin[0] = cal.DecalStripToStrip(sigcorr[0],neg);
  siglin[1] = cal.DecalStripToStrip(sigcorr[1],pos);
  
  // Remove gain part of the strip-to-strip calibration.
  // These aren't really photo-electrons.
  double blue[2];
  blue[0] = cal.GetPhotoElectrons(siglin[0],neg);
  blue[1] = cal.GetPhotoElectrons(siglin[1],pos);

  // Half the photons go each way.
  if(isnan(blue[0])||isnan(blue[1])) {
    // Feh? How did this happen?
    MSG("Photon",Msg::kFatal) << "Got an expectation value of NaN photons in CalibratedBlueComputer.\n";
    outCount.SetPeNegPos(0,0); // Explicitly do nothing.  
    return false;
  }
  
  // Set to a poisson number about this mean.
  outCount.SetBlueNegPos( (int)fRandom->PoissonD(blue[0]),
                          (int)fRandom->PoissonD(blue[1])
                        );
    
  MSG("Photon",Msg::kInfo) << "CalibratedBlueComputer: " 
                           << dE/Munits::MeV  << " MeV,  "
                           << mblue << " Blue nominal,   "
                           << sigmap[0] + sigmap[1] << " SigMap, "
                           << sigcorr[0] + sigcorr[1] << " SigCorr, "
                           << siglin[0] + siglin[1] << " SigLin, "
                           << blue[0] + blue[1] << " PE" << endl;
  MSG("Photon",Msg::kInfo) << "                      " 
                           << outCount.GetBlue_Neg() 
                           << " (" << blue[0] << ") pe Neg" 
                           << outCount.GetBlue_Pos() 
                           << " (" << blue[1] << ") pe Pos" 
                           << endl;
  return true;

}

---