SimPmtM64ToyNL.cxx

Go to the documentation of this file.

#include "SimPmtM64ToyNL.h"
#include "SimPmtM64CrosstalkTable.h"
#include "SimPmtMaker.h"
#include "SimPixelTimeBucket.h"
#include "MessageService/MsgService.h"
//#include "MessageService/MsgFormat.h"
#include "TMath.h"
#include <cassert>
#include <iostream>

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

CVSID("$Id: SimPmtM64ToyNL.cxx,v 1.1 2005/06/06 20:36:52 litchfld Exp $");

SimPmtMakerProxy<SimPmtM64ToyNL> gSimPmtM64ToyNLProxy("SimPmtM64ToyNL");

ClassImp(SimPmtM64ToyNL)

//STACICS
Float_t SimPmtM64ToyNL::fsNlSlope = 4.32e-4; //from fit
Float_t SimPmtM64ToyNL::fsNlTurnOn = 15.8; //from fit
Float_t SimPmtM64ToyNL::fsNlSaturate = 500;   //a lot of pes


//CONSTRUCTOR

SimPmtM64ToyNL::SimPmtM64ToyNL(PlexPixelSpotId tube, 
                                 VldContext context,
                                 TRandom* random ): 
  SimPmtM64Full(tube,context,random)
{
  if (fNPixels > 64) {
    // If you want to do exotic things like making M64s with more than 64 
    // pixels then sorry, but you're on your own.
    MSG("Detsim",Msg::kFatal) 
      << "M64 requested with " << fNPixels << "pixels. What?" << endl;
    assert(0);  
    }
}


//......................................................................
void SimPmtM64ToyNL::SimulateNonlinearity()
{
  Float_t totalQ[65];  //nasty, hardcoded value of fNPixels+1.
  
  for (UInt_t pix = 1 ; pix <= 64 ; ++pix ) totalQ[pix] = 0.0;

  //iterate over time buckets
  for (SimPmtBucketIterator it(*this); !it.End(); it.Next()) {
    SimPmtTimeBucket& bucket = it.Bucket();

  //iterate over pixels
    for (Int_t pix =1; pix <= fNPixels; pix++) {
      SimPixelTimeBucket& pixelBucket = bucket.GetPixelBucket(pix);
      
      //add charge in this pixelBucket to the charge on the pixel
      totalQ[pix] += pixelBucket.GetCharge();
      
    } //pixels
  } //buckets
  
  // calculate nonlinearity factor for each pixel, based on 
  // the summed charge (do it outside the bucket loop)
  
  //cout << totalQ[22] / (0.8 *1e6 *Munits::e_SI) << endl;
 Float_t ratioNL[65];
 for (Int_t pix = 1 ; pix <= fNPixels ; ++pix ) {
   ratioNL[pix] = GenNonlinearRatio(pix, totalQ[pix]);
  }
  
  //iterate over timebuckets > pixels again
  for (SimPmtBucketIterator it(*this); !it.End(); it.Next()) {
    SimPmtTimeBucket& bucket = it.Bucket();

    for (Int_t pix =1; pix <= fNPixels; pix++) {
      SimPixelTimeBucket& pixelBucket = bucket.GetPixelBucket(pix);

      // Reduce the charge inthis pixelBucket by the corresponding
      // nonlinearity ratio
      pixelBucket.SetCharge( ratioNL[pix] * pixelBucket.GetCharge() );
      
    } //pixels
  } //buckets
  
  return;
}



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

Float_t SimPmtM64ToyNL::GenNonlinearCharge( Int_t pixel,
                                             Float_t linCharge) 
{ 
  // Override of method in SimPmt. Using it is dangerous---does 
  // linCharge mean the charge in the bucket or the total
  // charge on the pixel? More useful is GenNonlinearRatio.
 
  return linCharge * GenNonlinearRatio(pixel, linCharge);
}
//.............................................................

Float_t SimPmtM64ToyNL::GenNonlinearRatio( Int_t /* pixel */,
                                            Float_t linCharge) 
{
  Float_t q_in_pe = linCharge / (0.8 * 1e6 * Munits::e_SI); 
  
 //Parametrisation NL of ND readout system with GC-LI

  //New parameterisation
  // - no nonlinearity below fsNlTurnOn
  // - saturating (1/Q) above fsNlSaturate
  // - linear in the intevening range

  if ( q_in_pe > fsNlTurnOn ){
    if (q_in_pe <= fsNlSaturate){
      // use form y = 1-m(x-t)
      return 1 - fsNlSlope * (q_in_pe - fsNlTurnOn);
    
    } 
    else {
      // use form y = a/(q-b) 
      // matching f & df/dq at the boundary [Upper case]
      // gives 
      //   a = F^2 / m   and   b = Q - F/m
      Float_t F = 1 - fsNlSlope * (fsNlSaturate - fsNlTurnOn);
      Float_t a = F*F / fsNlSlope;
      Float_t b = fsNlSaturate - F / fsNlSlope;
      return a / (q_in_pe - b);
      
    } // else (above saturation)
  } //if(above turnOn)
  
  return 1;

  
  //---Old parameteriastion--- 
  //  ~15% wrong
  //  unnessecary use of 5th order polynomial

  //   const double pol1[6]={  0,  //go to zero nonlinearity at zero pe      
  //                      3.71E-4,
  //                      -1.86E-5,
  //                      1.69E-7,
  //                      -6.54E-10,
  //                      9.06E-13 };
  //
  //
  //   const double pol2[2]={  0.011,
  //                      -4.65e-4};
  //
  //   Float_t fracLost = 0;
  //
  //      if(q_in_pe>=500) {
  //        //Flatten at some high value so charge dosen't go -ive
  //   fracLost =  pol2[0] 
  //       + pol2[1] * 500. ;
  //      }
  //      else if(q_in_pe>100 && q_in_pe < 500){
  //        fracLost = pol2[0] 
  //            + pol2[1] * q_in_pe;
  //      }
  //      else if(q_in_pe < 100){
  //        fracLost = pol1[0] 
  //            + pol1[1] * q_in_pe
  //            + pol1[2] * TMath::Power(q_in_pe,2)
  //            + pol1[3] * TMath::Power(q_in_pe,3)
  //            + pol1[4] * TMath::Power(q_in_pe,4)
  //            + pol1[5] * TMath::Power(q_in_pe,5);
  //      }
  //
  //     return (1 + fracLost);
  // ---======----


}

---