PhotonCalibratedPeComputer.cxx

Go to the documentation of this file.

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

CVSID( " $Id: PhotonCalibratedPeComputer.cxx,v 1.23 2008/04/04 20:55:32 rhatcher Exp $" );

ClassImp(PhotonCalibratedPeComputer)
PhotonTransportMakerProxy<PhotonCalibratedPeComputer>
        gPhotonCalibratedPeComputerProxy("photonCalibratedPeComputer"); 


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

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

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

  bool ok = true;
  ok = ok && r.Get("OverallLightOutput",fOverallLightOutput);
  ok = ok && r.Get("BirksConstant",fBirksConstant);
  ok = ok && r.Get("GeVPerMip",fGeVPerMip);
  ok = ok && r.Get("ScintillatorClipping",fScintillatorClipping);
  ok = ok && r.Get("MaxPePerHit",fMaxPePerHit);
  if (!ok) MSG("Photon",Msg::kError)
             << "Trouble configuring PhotonCalibratedPeComputer." << endl;
}

// Reset the model.
void 
PhotonCalibratedPeComputer::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 
PhotonCalibratedPeComputer::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;

  // Get position.
  double x = 0.5*(inHit->X1() + inHit->X2());
  
  double clippedfrac = 1.0;

  UgliGeomHandle ugli(fContext);
  UgliStripHandle ustrip=ugli.GetStripHandle(inHit->StripEndId());

  if (fScintillatorClipping) {
    // Clip the photon response by shaving off photons lost out the
    // ends of the scintillator either at the planed edges or the bore hole.

    // To do this, we cheat a bit: we assume that ALL strip ends play a role
    // in diminishing the return. 
   
    double distToEnd;
    
    distToEnd = ustrip.GetHalfLength()-x; // +ve end.
    clippedfrac *= ClipFraction( distToEnd ); // End1

    distToEnd = x+ustrip.GetHalfLength(); // -ve end.
    clippedfrac *= ClipFraction( distToEnd ); // End2

    if (ustrip.WlsBypass() > 0.) {
      // There is a coil bypass, so look at coil stripends.
      // Where is the coil center?
      double coilpos = 0.5*(ustrip.PartialLength(StripEnd::kNegative) - ustrip.PartialLength(StripEnd::kPositive));
      if (x<coilpos) {
        distToEnd = (-x) - ( ustrip.GetHalfLength() - ustrip.PartialLength(StripEnd::kNegative));
        clippedfrac *= ClipFraction( distToEnd );

      } else {
        distToEnd = (x) - ( ustrip.GetHalfLength() - ustrip.PartialLength(StripEnd::kPositive));
        clippedfrac *= ClipFraction( distToEnd );
      }
    }
  }

  if (clippedfrac < 0.49) {
    MSG("Photon",Msg::kError) << "Error: saw clipfrac less than 0.5: " << clippedfrac<< endl;
    clippedfrac = 0.5;
  }
 
  PlexStripEndId end[2];
  end[0] = inHit->StripEndId(); end[0].SetEnd(StripEnd::kNegative);
  end[1] = inHit->StripEndId(); end[1].SetEnd(StripEnd::kPositive);

  // Energy->MIP
  // Apply normalization and birks' constant.
  double emips = 
    clippedfrac                       // Light that didnt' leak out ends of strip                   
    * fOverallLightOutput             // Tuning parameter.
    * dE/(1. + fBirksConstant*dEdx)   // Birk's law
    / fGeVPerMip;                     // Convert from energy to mips

  
  // Useful for deugging:
  //MSG("Photon",Msg::kInfo) << "dE: " << dE << "  fGeVPerMip: " << fGeVPerMip << "  Birk suppress: " << 1.0/(1. + fBirksConstant*dEdx) << "  clippfrac: " << clippedfrac
  //     << " gives emips = " << emips << endl;

  const Calibrator& cal = Calibrator::Instance();
  // MIP->SigMap.  These two values should be the same (in principle)
  // but we may as well be pedantic.
  double sigmap[2] = {0,0};
  double sigcorr[2]= {0,0};
  double siglin[2] = {0,0};
  double adc[2]    = {0,0};
  double meanpe[2] = {0,0};
  double pe[2]     = {0,0};
  int npe[2]       = {0,0};

  if (ustrip.IsMirrored(StripEnd::kNegative) || ustrip.IsMirrored(StripEnd::kPositive) ) {
    // All of the energy goes to one end.
  } else {
    // The energy is split so that light goes both ways.
    emips *= 0.5;
  }

  for (int i=0;i<2;i++) {
    if ( ustrip.IsMirrored(end[i].GetEnd()) ) {
      // This end doesn't read out. Set the output to zero.
      sigmap[i] = sigcorr[i] = siglin[i] = pe[i] = 0;
      npe[i] = 0;
    } else {
      // Go through the decalibrtation chain:

      // Mips -> sigmaps.
      sigmap[i] =  cal.DecalMIP( emips , end[i]);

      // SigMap->SigCorr
      sigcorr[i] = cal.DecalAttenCorrected(sigmap[i],x,end[i]);

      // SigCorr->SigLin (siglin = ADCs for no drift, no linearity)
      siglin[i] = cal.DecalStripToStrip(sigcorr[i],end[i]);

      // SigLin->ADCs (Correct for drift. Linearity is applied during PMT simulation.)
      adc[i] = cal.DecalDrift(siglin[i],end[i]);

      // SigLin->Pe
      meanpe[i] = cal.GetPhotoElectrons(adc[i],end[i]);

      if (isnan(meanpe[i])) {
        MSG("Photon",Msg::kError) 
          << "Got an expectation value of NaN photons in CalibratedPeComputer.\n"
          << "\tHit on stripend " << end[i].AsString() << endl
          << "\tPosition: " << x << " m (in strip coords) " << endl
          << "\t" << dE/Munits::MeV  << " MeV,  "<< endl
          << "\t" << emips      << " MIP,   "<< endl
          << "\t" << sigmap[0]  << " / " << sigmap[1]  << " SigMap, "<< endl
          << "\t" << sigcorr[0] << " / " << sigcorr[1] << " SigCorr, "<< endl
          << "\t" << siglin[0]  << " / " << siglin[1]  << " SigLin, "<< endl
          << "\t" << adc[0]  << " / " << adc[1]  << " ADC, "<< endl
          << "\t" << meanpe[0]  << " / " << meanpe[1]  << " PE" << endl;
        outCount.SetPeNegPos(0,0); // Explicitly do nothing.  
        return false;              // and give up on this hit.
      }

      // Set to a poisson number about this mean.
      pe[i] = fRandom->PoissonD(meanpe[i]);

      // Handle some error conditions:
      if (isnan(pe[i])) {
        MSG("Photon",Msg::kError)
          << "Got a random poisson value of NaN photons in CalibratedPeComputer.\n"
          << "\tHit on stripend " << end[i].AsString() << endl
          << "\tPosition: " << x << " m (in strip coords) " << endl
          << "\t" << dE/Munits::MeV  << " MeV,  "<< endl
          << "\t" << emips      << " MIP,   "<< endl
          << "\t" << sigmap[0]  << " / " << sigmap[1]  << " SigMap, "<< endl
          << "\t" << sigcorr[0] << " / " << sigcorr[1] << " SigCorr, "<< endl
          << "\t" << siglin[0]  << " / " << siglin[1]  << " SigLin, "<< endl
          << "\t" << adc[0]  << " / " << adc[1]  << " ADC, "<< endl
          << "\t" << meanpe[0]  << " / " << meanpe[1]  << " PE" << endl;
        outCount.SetPeNegPos(0,0); // Explicitly do nothing.  
        return false;              // and give up on this hit.
      }
        
      if (pe[i] < 0) {
        MSG("Photon",Msg::kError) 
          << "Got -ve PEs on end " << i << endl
          << "\tHit on stripend " << end[i].AsString() << endl
          << "\tPosition: " << x << " m (in strip coords) " << endl
          << "\t" << dE/Munits::MeV  << " MeV,  "<< endl
          << "\t" << emips      << " MIP,   "<< endl
          << "\t" << sigmap[0]  << " / " << sigmap[1]  << " SigMap, "<< endl
          << "\t" << sigcorr[0] << " / " << sigcorr[1] << " SigCorr, "<< endl
          << "\t" << siglin[0]  << " / " << siglin[1]  << " SigLin, "<< endl
          << "\t" << adc[0]  << " / " << adc[1]  << " ADC, "<< endl
          << "\t" << meanpe[0]  << " / " << meanpe[1]  << " PE" << endl;
        pe[i] = 0;
      }
  
      if (pe[i] > fMaxPePerHit) {
        MSG("Photon",Msg::kError) 
          << "Got huge pulse on end " << i << endl
          << "\tHit on stripend " << end[i].AsString() << endl
          << "\tPosition: " << x << " m (strip coords) " << endl
          << "\t" << dE/Munits::MeV  << " MeV,  "<< endl
          << "\t" << emips      << " MIP,   "<< endl
          << "\t" << sigmap[0]  << " / " << sigmap[1]  << " SigMap, "<< endl
          << "\t" << sigcorr[0] << " / " << sigcorr[1] << " SigCorr, "<< endl
          << "\t" << siglin[0]  << " / " << siglin[1]  << " SigLin, "<< endl
          << "\t" << adc[0]  << " / " << adc[1]  << " ADC, "<< endl
          << "\t" << meanpe[0]  << " / " << meanpe[1]  << " PE" << endl
          << "\t ... Truncating to " << fMaxPePerHit << " pes." << endl;
        pe[i] = fMaxPePerHit;
      }
      npe[i] = (int) pe[i];
    }
  }


  // Set the output.
  outCount.SetPeNegPos( npe[0], npe[1] );

  /*
   // Consistency check. Disabled in production.
  float recomip[2];
  for (int i=0;i<2;i++) {
    float gain, width;
    cal.DecalGainAndWidth(gain, width, end[i]);
    float q = (float)(meanpe[i])*gain;
    q = cal.GetStripToStripCorrected(q,end[i]);
    q = cal.GetAttenCorrected(q,x,end[i]);
    q = cal.GetMIP(q,end[i]);
    recomip[i] = q;
  }
  MSG("Photon",Msg::kInfo) 
    << "Input mip: " << emips << "  Output: " << recomip[0] << "  " << recomip[1] << endl;
  */

  MSG("Photon",Msg::kDebug) << "CalibratedPeComputer: " 
                           << dE/Munits::MeV  << " MeV,  "
                           << emips << " MIP,   "
                           << sigmap[0] + sigmap[1] << " SigMap, "
                           << sigcorr[0] + sigcorr[1] << " SigCorr, "
                           << siglin[0] + siglin[1] << " SigLin, "
                           << meanpe[0] + meanpe[1] << " PE" 
                           << endl;
  MSG("Photon",Msg::kDebug) << "                      " 
                           << outCount.GetPe_Neg() 
                           << " (" << pe[0] << ") pe Neg   " 
                           << outCount.GetPe_Pos() 
                           << " (" << pe[1] << ") pe Pos" 
                           << endl;

  return true;
}

Double_t PhotonCalibratedPeComputer::ClipFraction(double x) const
{
  // Returns the integral of the double-lorentzian blue photon profile
  // from -inf to x.
  
  // Same parameters as in PhotonFastBlueModel.

  const double kLorentzWidth1 = 0.017853;
  const double kArea1 = 8.178052;
  const double kLorentzWidth2 = 0.001309;
  const double kArea2 = 0.864215;

  const double kArea = kArea1+kArea2;

  double frac = kArea1 * atan(x/kLorentzWidth1)
    + kArea2 * atan(x/kLorentzWidth2);
  frac = frac / (kArea * TMath::Pi());
  frac += 0.5;

  // Bonus extra fudge-factor:
  // We don't want this effect to be seen outside ~20 cm 
  // (where the effect is already down around 1%)
  // So, supress by an additional amount:
  if (x>0.2) frac = 1-(1-frac)*exp(-(x-0.2)*20.);

  return frac;
}

---