/search.css" rel="stylesheet" type="text/css"/> /search.js">
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

In This Package:

Classes | Public Member Functions | Static Public Member Functions | Private Attributes
RollingGain Class Reference

#include <RollingGain.h>

Collaboration diagram for RollingGain:
Collaboration graph
[legend]

List of all members.

Classes

struct  DetProp
struct  PmtProp

Public Member Functions

 RollingGain (const std::string &name, ISvcLocator *pSvcLocator)
virtual ~RollingGain ()
virtual StatusCode initialize ()
virtual StatusCode execute ()
virtual StatusCode finalize ()
void Fit (Detector det)
int InitPmtProp (const ServiceMode &svc)
int ResetPmtProp (Detector det)

Static Public Member Functions

static double poisson (double u, int n)
static double gn (double x, double n, double q0, double sigma0, double q1, double sigma1)
static double NIMmodel (double *xpar, double *par)
static double Pedestal (double *xpar, double *par)

Private Attributes

ICableSvcm_cableSvc
 Cable mapping service.
IStatisticsSvc * m_statsSvc
IFloatingFeePedestalSvcm_pedesSvc
 Average Pedestal.
TF1 * f1
TF1 * f2
TFile * m_rootfile
TTree * m_tree
FILE * m_masterfile
FILE * m_outfile
string m_fileName
map< Detector, DetPropm_DetPropMap
int m_Fix
int m_output
int m_NStop
int m_lastFit
int m_site
int m_simFlag
TimeStamp m_currentTime
int m_thisRun
double m_fitPeriod
 Both m_fitPeriod and m_fitStatistic must be satisfied.
double m_fitStatistic
list __all__ = [ "run" ]

Detailed Description

Definition at line 66 of file RollingGain.h.


Constructor & Destructor Documentation

RollingGain::RollingGain ( const std::string &  name,
ISvcLocator *  pSvcLocator 
)

Definition at line 15 of file RollingGain.cc.

                                                                       :
  GaudiAlgorithm(name, pSvcLocator),
  m_statsSvc(0)
{
  declareProperty("FileName", m_fileName = "AdcFit.root", "Root file name");  
  declareProperty("FitPeriod", m_fitPeriod = 60*60, "Fit period in second, Default is 60 minutes");
  declareProperty("FitStatistic", m_fitStatistic = 3000, "Fit period in statistic, average statistic. Set to <=0 to ignore");
  declareProperty("Fix", m_Fix = 1, "a temp fix for 40M noise");
  declareProperty("RGOutput",m_output=1,"output rolling gain file");
  declareProperty("LastFit", m_lastFit = 1, "Fit the last set of data. Enforced output");
}
RollingGain::~RollingGain ( ) [virtual]

Definition at line 26 of file RollingGain.cc.

{
}

Member Function Documentation

StatusCode RollingGain::initialize ( ) [virtual]

Cable service

Statistics service

FloatingFeePedestalSvc

Definition at line 590 of file RollingGain.cc.

{
  StatusCode sc;
  info()<<"FileName "<<m_fileName<<endreq;
  info()<<"FitPeriod "<<m_fitPeriod<<"s"<<endreq;
  info()<<"FitStatistic "<<m_fitStatistic<<endreq;
  info()<<"Fix "<<m_Fix<<endreq;
  info()<<"RGOutput "<<m_output<<endreq;
  info()<<"LastFit "<<m_lastFit<<endreq;
  info()<<"Random Trigger added for gain fitting"<<endreq;

  sc = this->GaudiAlgorithm::initialize();
  if( sc.isFailure() ) {
    error() << "Base class initialization error" << endreq;
    return sc;
  }

  sc = service("CableSvc",m_cableSvc);
  if( sc.isFailure() ) {
    error() << "Can't get CableSvc" << endreq;
    return sc;
  }

  sc = service("StatisticsSvc",m_statsSvc,true);                        
  if(sc.isFailure()){                                                                    
    error() << "Failed to get StatisticsSvc" << endreq;                                  
    return sc;
  } 

  sc = service("FloatingDaqFeePedestalSvc",m_pedesSvc);
  if( sc.isFailure() ) {
    error() << "Can't get FloatingDaqFeePedestalSvc" << endreq;
    return sc;
  }

  f1 = new TF1("f1","gaus",20,250);
  //f1 = new TF1("f1",Pedestal,0,150,5);
  f2 = new TF1("f2",NIMmodel,20,300,8);

  if(m_output){
    m_rootfile = new TFile(m_fileName.c_str(),"recreate");
    m_masterfile = fopen("Master.pmtCalibMap.txt", "a");
    fprintf(m_masterfile, "# [StartRun] [EndRun] [Sites] [SimFlag] [StartTimeSec] [StartTimeNanoSec] [EndTimeSec] [EndTimeNanoSec] [map file]\n");
    fclose(m_masterfile);
  }
  else{
    //book the tree for DQ summary value
    m_tree = new TTree("DQTree", "Data Quality Summary tree");
    //m_tree->SetMaxTreeSize(200000000);  // 200M
  }

  return sc;
}
StatusCode RollingGain::execute ( ) [virtual]

Get ReadoutHeader

for a pedestal reference

Confirm detector type

check AD or water pool

Get the information for current detector

Checking the crate info

Cut: Block ESum-only trigger, always noisy

Cut: Delta T cut (Save time. Singles rate may be wrong. But that needs to decoupled from this later.)

Count number of available triggers

Loop over channel data

Loop over channel data; calculate NoiseMulti (Multiplicity in noise window)

get the first tdc

Loop over channel data; hit filling

Remove 40MHz noise (DocDB 5529)

1. Including NHit trigger

Excluding Random trigger

Out of main event trigger window

Noise Mulitiplicity <3. Remove below trigger events

2. Just for random trigger

Other cuts: Only one hit per channel

High Gain i.e. Fine Range

Remove irregular hits with nonstandard peak cycle (DocDB 6819 & 6973)

Remove events with bad baseline

Remove 40MHz noise (DocDB 5527)

Definition at line 647 of file RollingGain.cc.

{
  StatusCode sc;
  
  ReadoutHeader* readoutHdr = get<ReadoutHeader>( ReadoutHeader::defaultLocation() );
  if(!readoutHdr) {
    error()<<"Failed to get readout header"<<endreq;
    return StatusCode::FAILURE;
  }

  sc = m_pedesSvc->cumulate( readoutHdr );

  const DaqCrate* daqCrate = readoutHdr->daqCrate();
  if(!daqCrate) {
    error()<<"Failed to get daqCrate from header"<<endreq;
    return StatusCode::FAILURE;
  }
  Detector det = daqCrate->detector();
  if( !(det.isAD() ||
        det.isWaterShield()) ) {
    //info()<<"this is not PMT data, skip to next"<<endreq;
    return StatusCode::SUCCESS;
  }

  static bool first = true;
  if( first ) {
    ServiceMode svc(readoutHdr->context(), 0);
    InitPmtProp(svc); // Initializes all detectors at site

    Context ctx = readoutHdr->context();
    m_simFlag=ctx.GetSimFlag(); 

    first = false;
  }
  
  m_currentTime = daqCrate->triggerTime();

  DetProp& currDet = m_DetPropMap[det];
  if( currDet.FirstTrigTime == TimeStamp::GetBOT() ) {
    currDet.FirstTrigTime = m_currentTime;
    currDet.LastTrigTime  = m_currentTime;
  }
  TimeStamp PrevTrigTime = currDet.LastTrigTime;
  currDet.LastTrigTime = m_currentTime;

  const DaqPmtCrate* pmtCrate = daqCrate->asPmtCrate();
  if (!pmtCrate) {
    error()<<"Failed to get DaqPmtCrate"<<endreq;
    return StatusCode::FAILURE;
  }

  currDet.LastRun = pmtCrate->runNumber();
  ServiceMode svcMode( readoutHdr->context(),0 );
  
  Trigger::TriggerType_t trigType = daqCrate->triggerType();
  if( trigType == Trigger::kESumAll )   return StatusCode::SUCCESS;

  if( (m_currentTime - PrevTrigTime).GetSeconds() < 20e-6 )  return StatusCode::SUCCESS;

  currDet.NTrigger++;

  const DaqPmtCrate::PmtChannelPtrList& channels = pmtCrate->pmtChannelReadouts();
  DaqPmtCrate::PmtChannelPtrList::const_iterator ci, ci_end = channels.end();

  int NoiseMulti=0;
  for(ci = channels.begin(); ci!=ci_end; ci++) {
    const DaqPmtChannel* channel = *ci;
    // number of hits in this channel
    unsigned int multi = channel->hitCount();
    int tdc;
    if(multi>=1)  {
      tdc = channel->tdc(0);
      if(tdc>1070)  NoiseMulti += 1;
    }
  }
  
  for(ci = channels.begin(); ci!=ci_end; ci++) {
    const DaqPmtChannel* channel = *ci;
    FeeChannelId channelId = channel->channelId();
    
    // number of hits in this channel
    unsigned int multi = channel->hitCount();
    int id = 0;
    if(det.isAD()){
      AdPmtSensor pmtId = m_cableSvc->adPmtSensor( channelId, svcMode );
      id = pmtId.fullPackedData();
    }
    else{
      PoolPmtSensor pmtId = m_cableSvc->poolPmtSensor( channelId, svcMode );
      id = pmtId.fullPackedData();
    }
    map<int, PmtProp>::iterator mapit = currDet.PmtPropMap.find(id);
    if(mapit==currDet.PmtPropMap.end()) continue; //skip the one not in map
    PmtProp& curPmt = mapit->second;

    /* Add Channel ID */
    curPmt.ChannelId = channelId.fullPackedData();

    for(unsigned int nr = 0; nr<multi; ++nr) {
      int tdc = channel->tdc(nr);
      int adc = channel->adc(nr);
      float preAdc = channel->preAdcAvg(nr);
      int gain = channel->isHighGainAdc(nr) ? FeeGain::kHigh
                                            : FeeGain::kLow ;
      int peakCycle = channel->peakCycle(nr);
      double avePed = m_pedesSvc->pedestal( channelId, gain );

      // TDC
      curPmt.h_Tdc->Fill(tdc);

      // Noise
      if(tdc>1070) {
        if(m_Fix) {
          int residue = (tdc-15)%16;                
          if(residue >= 4) {
            curPmt.DarkRate++;
          }
          if(residue == 1) {
            curPmt.ElecNoiseRate++;
          }
        }else{
          curPmt.DarkRate++;
        }
      }

      // Accumulate hits for gain fitting
      // Incorporate regular triggers and random triggers
      if( ( ((trigType & Trigger::kMult) == Trigger::kMult) &&            
            ((trigType & Trigger::kRandom) != Trigger::kRandom) &&        
            (tdc>1070) &&                                                 
            (NoiseMulti<3) ) ||                                           
          
          ( ((trigType & Trigger::kRandom) == Trigger::kRandom) ) ) {     
        
        if( multi==1 &&                               
            gain==FeeGain::kHigh &&                   
            peakCycle>=4 && peakCycle<=6 &&           
            fabs( preAdc-avePed ) <20                 
            ) {
          if(m_Fix){                //  special treat for run 5773 and later runs, lose one fourth data.
            int residue = (tdc-15)%16;                
            if(residue>=4){
              curPmt.h_Adc->Fill(adc-preAdc);
              curPmt.h_PreAdc->Fill(preAdc);        
              currDet.Statistics++;
            }
          }
          else{
            curPmt.h_Adc->Fill(adc-preAdc);
            curPmt.h_PreAdc->Fill(preAdc);
            currDet.Statistics++;
          }
        }
      }
    }
  } // end of channels
    
  //-------------------------------------
  //fit the data
  int NChannels=currDet.PmtPropMap.size();
  int TotalHits=currDet.Statistics;
  double aveSta= NChannels>0 ? TotalHits/NChannels : 0 ;
  
  if( aveSta>m_fitStatistic && 
      (m_currentTime - currDet.FirstTrigTime).GetSeconds() > m_fitPeriod )  { // Fit period test and statistic test
    //Fit
    Fit(det);
    //after fit
    ResetPmtProp(det);
    m_NStop++;
  }
  
  return StatusCode::SUCCESS;
}
StatusCode RollingGain::finalize ( ) [virtual]

Definition at line 835 of file RollingGain.cc.

{
  map<Detector, DetProp>::iterator it_det, it_det_end=m_DetPropMap.end();
  for( it_det=m_DetPropMap.begin(); it_det!=it_det_end; it_det++ ) {
    DetProp & currDet = it_det->second;
    Detector det = it_det->first;

    if( m_lastFit ) {
      if( currDet.Statistics>0 ) {
        Fit(det);
        m_NStop++;
      }
    }
        
    map<int/*PmtId*/,PmtProp>::iterator it,idend=currDet.PmtPropMap.end();
    for(it=currDet.PmtPropMap.begin(); it!=idend; it++) {
      PmtProp& curr = it->second;
      delete curr.h_Adc;
      delete curr.h_PreAdc;
      delete curr.h_Tdc;
    }
    currDet.PmtPropMap.clear();
  }

  if(m_output){
    m_rootfile->Write();
    m_rootfile->Close();
    delete m_rootfile;
  }
  else{
    //Construct histogram path in statistics service
    std::stringstream path;
    path << "/file1/diagnostics/run_" << std::setfill('0') << std::setw(7)
         << m_thisRun << "/DQTree";
    m_statsSvc->put( path.str(), m_tree );
  }
  delete f1;
  delete f2;
  //delete m_tree;

  StatusCode sc;
  sc = this->GaudiAlgorithm::finalize();
  return sc;
}
double RollingGain::poisson ( double  u,
int  n 
) [static]

Definition at line 30 of file RollingGain.cc.

                                          {
  return pow(u,n)*exp(-u)/TMath::Gamma(n+1);
}
double RollingGain::gn ( double  x,
double  n,
double  q0,
double  sigma0,
double  q1,
double  sigma1 
) [static]

Definition at line 34 of file RollingGain.cc.

                                   {
  double sigman = sqrt(sigma0*sigma0+n*sigma1*sigma1);
  double xn = x-q0-n*q1;
  double tmp = pow(xn,2)/(2*sigman*sigman);
  return exp(-tmp)/(sigman*sqrt(2*TMath::Pi()));
}
double RollingGain::NIMmodel ( double *  xpar,
double *  par 
) [static]

Definition at line 42 of file RollingGain.cc.

                                                    {
  double x =      xpar[0];
  double norm =   par[0];
  double q0 =     par[1];
  double sigma0 = par[2];
  double q1 =     par[3];  
  double sigma1 = par[4];
  double w =      par[5];
  double a =      par[6];
  double u =      par[7];
  int nmin =   1;
  int nmax =   6;

  double total = 0;
  for(int i=nmin; i<nmax; i++){
    int n = i;
    double tmp = poisson(u,n);
    double tmp2 = gn(x,n,q0,sigma0,q1,sigma1);
    double qn = q0+n*q1;
    double sigman = sqrt(sigma0*sigma0+n*sigma1*sigma1);

    double tmp3 = exp(-a*(x-qn-a*sigman*sigman/2.))*a/2.;
    double tmp4 = TMath::Erf(fabs(q0-qn-a*sigman*sigman)/(sigman*sqrt(2.)));
    double sign = 1.0;
    if((x-qn-sigman*sigman*a)<0) sign = -1.0;
    tmp4 += sign*TMath::Erf(fabs(x-qn-sigman*sigman*a)/(sigman*sqrt(2.)));
    tmp3 *= tmp4;
    
    double result = tmp*((1-w)*tmp2+w*tmp3);
    total += result;
  }
  return total*norm;
}
double RollingGain::Pedestal ( double *  xpar,
double *  par 
) [static]

Definition at line 76 of file RollingGain.cc.

                                                    {
  double x =      xpar[0];
  double norm =   par[0];
  double q0 =     par[1];
  double sigma0 = par[2];
  double w =      par[3];
  double a =      par[4];

  double tmp = exp(-pow(x-q0,2)/(2*pow(sigma0,2)));    
  double result = (1-w)*tmp/(sigma0*sqrt(2*TMath::Pi()));
  if(x>q0) result += w*a*exp(-a*(x-q0));
  
  return result*norm;
}
void RollingGain::Fit ( Detector  det)

Pick out the detector need to fit

For TDC [1070, 1085], 8 40MHz clock cycles

skip empty channels

Default values are ok to store.

skip low statisticd tubes

start a real fit

Step2: The gain fitting ^^^^^^^^^^^^^^^^^^^^^^^ Do a smarter fitting to automatically adjust some parameters until get a good fit So far it is found that lowerbound is the most sensitive parameter.

lowerbound trial is ordered by the chance a good fit can be reached.

uppperTrial is for the simple gauss fit only

Start right above pedestal then include the long tail.

Start right above pedestal then include the long tail.

Step3: Let's see if this channel can be bailed out by a simple gauss fit ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Start right above pedestal then include the long tail.

Last step: Fit result summary ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Definition at line 91 of file RollingGain.cc.

                                 {
    
  DetProp & currDet = m_DetPropMap[det];
  TimeStamp FirstTime = currDet.FirstTrigTime.GetSeconds();
  TimeStamp LastTime = currDet.LastTrigTime.GetSeconds();

  double gain,gainErr,darkRate,darkRateErr,elecNoiseRate,elecNoiseRateErr,preADC,preADCErr,chiSquare;
  int pmtID,fitStatus;

  if(m_output){
    ostringstream tmp;
    tmp <<setw(2) <<setfill('0') << m_NStop;
    string filename = "pmtDataTable_"+tmp.str()+".txt";
    m_outfile = fopen(filename.c_str(), "w+");
  
    fprintf(m_outfile, "# First Trigger time: (TrigSecond, TrigNanoSec) = (%d,%d).\n", (int)FirstTime.GetSec(), (int)FirstTime.GetNanoSec());
    fprintf(m_outfile, "# Last Trigger time: (TrigSecond, TrigNanoSec) = (%d,%d).\n", (int)LastTime.GetSec(), (int)LastTime.GetNanoSec()+1);
    fprintf(m_outfile, "# [pmtID] [channelID] [description] [status] [nhits] [fitStat] [chi2ndf] [speHigh] [gainErr] [sigmaSpe] [speLow] [timeOffset] [timeSpread] [efficiency] [prePulse] [afterPulse] [darkRate] [darkRateErr] [elecNoiseRate] [elecNoiseRateErr] [preAdc] [preAdcErr]\n");

    m_rootfile->cd();
    string newdir = tmp.str();
    m_rootfile->mkdir(newdir.c_str());
    m_rootfile->cd(newdir.c_str());

    m_masterfile = fopen("Master.pmtCalibMap.txt", "a");
    fprintf(m_masterfile, "0000000 0000000 %4.0d %4.0d %10.0d %10.0d %10.0d %10.0d %s\n",
            (currDet.Det.fullPackedData()>>16), m_simFlag, 
            (int)FirstTime.GetSec(), (int)FirstTime.GetNanoSec(),
            (int)LastTime.GetSec(),  (int)LastTime.GetNanoSec()+1, filename.c_str());
    fclose(m_masterfile);
  }
  else{
    //set branch for DQTree
    if(!m_tree->GetBranch("gain")){
      m_tree->Branch("gain",&gain,"gain/D");
      m_tree->Branch("gainErr",&gainErr,"gainErr/D");
      m_tree->Branch("darkRate",&darkRate,"darkRate/D");
      m_tree->Branch("darkRateErr",&darkRateErr,"darkRateErr/D");
      m_tree->Branch("elecNoiseRate",&elecNoiseRate,"elecNoiseRate/D");
      m_tree->Branch("elecNoiseRateErr",&elecNoiseRateErr,"elecNoiseRateErr/D");
      m_tree->Branch("preADC",&preADC,"preADC/D");
      m_tree->Branch("preADCErr",&preADCErr,"preADCErr/D");
      m_tree->Branch("pmtID",&pmtID,"pmtID/I");
      m_tree->Branch("fitStatus",&fitStatus,"fitStatus/I");
      m_tree->Branch("chiSquare",&chiSquare,"chiSquare/D");
    }
  }

  map<int/*PmtId*/,PmtProp>::iterator it,idend=currDet.PmtPropMap.end();
  for(it=currDet.PmtPropMap.begin(); it!=idend; it++) {
    PmtProp& curr = it->second;

    //calculate dark rate
    double tlength = 
      (double)currDet.NTrigger
      *115                   // For the range of TDC [1070, 1085], 115 TDCs
      *1/TdcFrequencyHz;     // 1 tdc corresponding 1.5625ns

    if(m_Fix) // 40 MHz noise fix applied, 25% hits not counted.
      {
        tlength*=0.75;
      }
    curr.DarkRateErr = sqrt(curr.DarkRate)/tlength;
    curr.DarkRate = curr.DarkRate/tlength;

    //calculate ElecNoiseRate
    if(m_Fix) { // 40MHz noise fix applied
      double noiselength = 
        (double)currDet.NTrigger
        *8                    
        *1/TdcFrequencyHz;

      // In principle the electronic noise is added to the normal PMT dark rate. Substract them.
      curr.ElecNoiseRateErr = sqrt(curr.ElecNoiseRate)/noiselength;
      curr.ElecNoiseRateErr = sqrt(curr.ElecNoiseRateErr*curr.ElecNoiseRateErr + curr.DarkRateErr*curr.DarkRateErr );
      curr.ElecNoiseRate = curr.ElecNoiseRate/noiselength;
      curr.ElecNoiseRate = curr.ElecNoiseRate - curr.DarkRate;
    } else {
      curr.ElecNoiseRateErr = 0;
      curr.ElecNoiseRate = 0;
    }

    if( curr.h_PreAdc->GetEntries() == 0 ) {          

      curr.Status=PmtCalibData::kUnknown;
    } else if ( (curr.h_PreAdc->GetEntries() > 0 &&           
                 curr.h_PreAdc->GetEntries() < 100) ||
                (curr.h_Adc->GetEntries() > 0 &&
                 curr.h_Adc->GetEntries() < 100) ) {

      curr.Status = PmtCalibData::kBadEfficiency;
      curr.NHits = (int)curr.h_Adc->GetEntries();
      //curr.FitStatus = 
      //curr.Chi2 = 
      curr.Gain = curr.h_Adc->GetMean();
      curr.Sigma = curr.h_Adc->GetRMS();
      curr.GainErr = curr.Sigma/TMath::Sqrt(curr.NHits);
      curr.Pedestal = curr.h_PreAdc->GetMean();
      curr.PedestalErr = (curr.h_PreAdc->GetRMS())/TMath::Sqrt(curr.h_PreAdc->GetEntries());
    } else {                                          
      // Step1: fit pedestal
      // ^^^^^^^^^^^^^^^^^^^
      ostringstream tmp3;
      tmp3<<setw(2)<<setfill('0')<< m_NStop;
      string histname = "PreAdc_"+curr.PmtLabel+"_"+tmp3.str();
      string histname2 = "Adc_"+curr.PmtLabel+"_"+tmp3.str();
      
      double pedarea = curr.h_PreAdc->GetEntries();
      double pedmean = curr.h_PreAdc->GetMean();
      double pedsigm = curr.h_PreAdc->GetRMS();
      f1->SetParameters(pedarea,pedmean,pedsigm);
      f1->SetRange(pedmean-3*pedsigm,pedmean+3*pedsigm);
      int status1 = curr.h_PreAdc->Fit("f1","RQ");
      
      double q0 = f1->GetParameter(1);
      double q0err = f1->GetParError(1);
      double sigma0 = f1->GetParameter(2);
      //double w = f1->GetParameter(3);
      //double a = f1->GetParameter(4);
      if(f1->GetChisquare()/f1->GetNDF()>10 || status1!=0){
        warning()<< "Batch "<<tmp3.str()<<" bad pedestal fit for PMT: "
                 << curr.PmtLabel <<endreq;;
        q0 = curr.h_PreAdc->GetMean();
        q0err = 0;
        sigma0 = 1.9;
        curr.Status |= PmtCalibData::kBadPedestal;
      }
      curr.Pedestal = q0;
      //curr.h_PreAdc->Write(histname.c_str());
      //----------------------------------------
      //Adc is pedestal subtracted, set q0=0
      q0 = 0;
      //----------------------------------------
      double AreaGuess, Q1MeanGuess, Q1SigmaGuess;
      AreaGuess = curr.h_Adc->GetEntries();
      Q1MeanGuess = curr.h_Adc->GetMean();
      Q1SigmaGuess = curr.h_Adc->GetRMS();
      
      double lowerbound, upperbound;
      double lowerTrial[9]={0.65*Q1MeanGuess,0.625*Q1MeanGuess,0.6*Q1MeanGuess,
                            0.55*Q1MeanGuess,0.5*Q1MeanGuess,0.45*Q1MeanGuess,
                            0.4*Q1MeanGuess,0.35*Q1MeanGuess,0.3*Q1MeanGuess};
      
      double upperTrial[9]={1.65*Q1MeanGuess,1.55*Q1MeanGuess,2.00*Q1MeanGuess,
                            1.65*Q1MeanGuess,1.55*Q1MeanGuess,2.00*Q1MeanGuess,
                            1.65*Q1MeanGuess,1.55*Q1MeanGuess,2.00*Q1MeanGuess};

      bool GotGoodFit = false;

      int status2;      // Fit status
      double q1;        // Gain
      double sigma1;    // Sigma
      double q1err;     // Gain Err
      double sigma1err; // Sigma Err
      double avechi2;   // Chi2/ndf

      for( unsigned int tidx=0; tidx<9; tidx++ ) {
        lowerbound = q0+lowerTrial[tidx];    
        upperbound = q0+Q1MeanGuess+5*Q1SigmaGuess;

        // Initialize fitting parameters
        f2->SetParameters(AreaGuess*10, 0, 0, Q1MeanGuess, Q1SigmaGuess, 0.01, 0.03, 0.1);
        f2->SetParNames("N","Q_{0}","#sigma_{0}","Q_{1}","#sigma_{1}","w","#alpha","#mu");
      
        f2->FixParameter(1,q0);
        f2->FixParameter(2,sigma0);
        f2->FixParameter(5,0.);
        f2->FixParameter(6,0.);
        //f2->FixParameter(5,w);
        //f2->FixParameter(6,a);
      
        f2->SetRange(lowerbound, upperbound); 
        status2 = curr.h_Adc->Fit("f2","RQ");
        
        q1 = f2->GetParameter(3);
        sigma1 = f2->GetParameter(4);
        q1err = f2->GetParError(3);
        sigma1err = f2->GetParError(4);
        avechi2 = f2->GetChisquare()/f2->GetNDF();
        if(sigma1<0) sigma1 = -sigma1;//sometimes fit gives a negative sigma
        
        // Keep trying until get a good fit.
        // Define a good fit
        // ^^^^^^^^^^^^^^^^^
        if( status2==0 &&                         // Root TH1 fit status 0 if ok (http://root.cern.ch/root/html/TH1.html#TH1:Fit)
            avechi2<10 &&                         // Good Chi2/ndf
            q1>6 && q1<80 &&                      // Too low or too high gain
            q1>lowerbound && q1<upperbound &&     // Within lower and upper bounds otherwise it always has big error
            q1err<10 &&                           // Reasonable err of Gain
            sigma1<q1 &&                          // Gain sigma is usually 1/3 of Gain
            sigma1>1 && sigma1<30                 // Too low or too high sigma
            ) {
          // Escape from the loop
          GotGoodFit = true;
          break;
        }
      }

      if( !GotGoodFit ) {
        info()<<"Batch "<<tmp3.str()<<" bad gain fit for PMT: "<<
          curr.PmtLabel<<", try a simple Gaussian fit." <<endreq;
        
        for( unsigned int tidx=0; tidx<9; tidx++ ) {
          lowerbound = q0+lowerTrial[tidx];    
          upperbound = q0+upperTrial[tidx];

          f1->SetParameters(AreaGuess,Q1MeanGuess,7.);
          f1->SetRange(lowerbound, upperbound);
          status2 = curr.h_Adc->Fit("f1","RQ");
          
          q1 = f1->GetParameter(1)-q0;
          sigma1 = f1->GetParameter(2);
          q1err = f1->GetParError(1);
          sigma1err = f1->GetParError(2);
          avechi2 = f1->GetChisquare()/f1->GetNDF();
          
          // Define a good fit
          // ^^^^^^^^^^^^^^^^^
          if( status2==0 &&                         // Root TH1 fit status 0 if ok (http://root.cern.ch/root/html/TH1.html#TH1:Fit)
              avechi2<10 &&                         // Good Chi2/ndf
              q1>6 && q1<80 &&                      // Too low or too high gain
              q1>lowerbound && q1<upperbound &&     // Within lower and upper bounds otherwise it always has big error
              q1err<10 &&                           // Reasonable err of Gain
              sigma1<q1 &&                          // Gain sigma is usually 1/3 of Gain
              sigma1>1 && sigma1<40                 // Too low or too high sigma
              ) {
            // Escape from the loop
            GotGoodFit = true;
            break;
          }
        }
      }

      if( !GotGoodFit ) {
        warning()<< "Batch "<<tmp3.str()<<" gaussian fit for PMT: "
                 << curr.PmtLabel
                 <<" is not good either. Set the gain with mean value of h_Adc."<<endreq;
        q1 = curr.h_Adc->GetMean();
        q1err = curr.h_Adc->GetMean() / TMath::Sqrt( curr.h_Adc->GetEntries() );
        sigma1 = 0;
        sigma1err = 0;
      }

      //q0 is  always 0 here
      //curr.Pedestal = q0; 
      curr.PedestalErr = q0err;
      curr.Gain = q1;
      curr.Sigma = sigma1;
      curr.GainErr = q1err;
      curr.SigmaErr = sigma1err;
      curr.NHits = (int)AreaGuess;
      curr.FitStatus = status2;
      curr.Chi2 = avechi2;
      if( !GotGoodFit ) {
        curr.Status |= PmtCalibData::kBadGain;
      }
      if( curr.Sigma/curr.Gain > 1 ) {
        curr.Status |= PmtCalibData::kBadGainWidth;
      }
      if( curr.Status == PmtCalibData::kUnknown) {
        curr.Status = PmtCalibData::kGood;
      }
    }
    
    
    if(curr.h_Adc->GetEntries()>0){//skip 0s    
      if(m_output){
        fprintf(m_outfile, "%d %d %s %2d %6d %2d %6.2f %8.3f %8.3f  %7.3f %.1f %.1f %.1f %.1f %.1f %.1f %9.1f %9.1f %8.3f %8.3f %8.3f %8.3f\n", 
                curr.PmtId, curr.ChannelId, curr.PmtLabel.c_str(), curr.Status, curr.NHits, curr.FitStatus,
                curr.Chi2, curr.Gain, curr.GainErr, curr.Sigma, curr.Gain/19 ,0.,0.,0.,0.,0.,curr.DarkRate, curr.DarkRateErr, curr.ElecNoiseRate, curr.ElecNoiseRateErr, curr.Pedestal, curr.PedestalErr);
      
        curr.h_PreAdc->Write();
        curr.h_Adc->Write();
        curr.h_Tdc->Write();
      }
      else{
        //fill DQTree
        pmtID = curr.PmtId;
        fitStatus = curr.FitStatus;
        chiSquare = curr.Chi2;
        gain = curr.Gain;
        gainErr = curr.GainErr;
        darkRate = curr.DarkRate;
        darkRateErr = curr.DarkRateErr;
        elecNoiseRate = curr.ElecNoiseRate;
        elecNoiseRateErr = curr.ElecNoiseRateErr;
        preADC = curr.Pedestal;
        preADCErr = curr.PedestalErr;
        m_tree->Fill();
      }
    }
  }
 

  if(m_output){
    fclose(m_outfile);
  }
  else{ 
    // Construct histogram path in statistics service
    //std::stringstream path;
    //path << "/file1/diagnostics/run_" << std::setfill('0') << std::setw(7)
    //     << currDet.LastRun << "/DQTree";
    //m_statsSvc->put( path.str() ,    m_tree );
    m_thisRun = currDet.LastRun;
  }
}
int RollingGain::InitPmtProp ( const ServiceMode svc)

Default pmt properties

Default detector maps

Start

Locate this detector

Get the PMT map

Locate this detector

Get the PMT map

Definition at line 409 of file RollingGain.cc.

{

  vector<AdPmtSensor> adPmts;
  vector<PoolPmtSensor> poolPmts;
  
  // Loop over all detectors at site, and add to PMT lists
  for(DetectorId::DetectorId_t detId=DetectorId::kUnknown;
      detId<=DetectorId::kAll;){
    Context context = svc.context();
    context.SetDetId(detId);
    ServiceMode svcByDet(context,svc.task());
    if(detId==DetectorId::kAD1 || detId==DetectorId::kAD2 
       || detId==DetectorId::kAD3 || detId==DetectorId::kAD4){
      const vector<AdPmtSensor>& adPmtsByDet = 
        m_cableSvc->adPmtSensors(svcByDet);
      adPmts.insert(adPmts.end(), adPmtsByDet.begin(), adPmtsByDet.end());
    }else if(detId==DetectorId::kIWS || detId==DetectorId::kOWS){
      const vector<PoolPmtSensor>& poolPmtsByDet = 
        m_cableSvc->poolPmtSensors(svcByDet);
      poolPmts.insert(poolPmts.end(), poolPmtsByDet.begin(), poolPmtsByDet.end());
    }
    detId = (DetectorId::DetectorId_t)(1+(int)detId);
  }

  info()<<"Initializing " << adPmts.size() << " AD PMTs." << endreq;
  info()<<"Initializing " << poolPmts.size() << " Pool PMTs." << endreq;

  PmtProp defPmtProp;
  defPmtProp.h_Adc=0;
  defPmtProp.h_Tdc=0;
  defPmtProp.h_PreAdc=0;
  defPmtProp.Pedestal=0;
  defPmtProp.PedestalErr=0;
  defPmtProp.DarkRate=0;
  defPmtProp.DarkRateErr=0;
  defPmtProp.Gain=0;
  defPmtProp.GainErr=0;
  defPmtProp.Sigma=0;
  defPmtProp.SigmaErr=0;
  defPmtProp.Status=PmtCalibData::kUnknown;
  defPmtProp.NHits=0;
  defPmtProp.FitStatus=-1;
  defPmtProp.Chi2=-1;
  defPmtProp.PmtLabel="";
  defPmtProp.PmtId=0;
  defPmtProp.ChannelId=0;

  DetProp defDetProp;
  defDetProp.Det.set(Site::kUnknown, DetectorId::kUnknown);
  defDetProp.PmtPropMap.clear();
  defDetProp.FirstTrigTime=TimeStamp::GetBOT();
  defDetProp.LastTrigTime=TimeStamp::GetBOT();
  defDetProp.NTrigger=0;
  defDetProp.Statistics=0;
  defDetProp.FirstRun=-1;
  defDetProp.LastRun=-1;

  m_DetPropMap.clear();
  map<Detector, DetProp>::iterator it;

  //define histograms
  //  AD PMTs 
  for( unsigned int ct=0; ct<adPmts.size(); ct++ ) {
    AdPmtSensor sensor = adPmts[ct];
    
    Detector det(sensor.site(), sensor.detectorId());
    it = m_DetPropMap.find( det );
    if( it == m_DetPropMap.end() ) {
      m_DetPropMap[det]=defDetProp;
    }
    m_DetPropMap[det].Det = det;

    map<int/*PmtId*/,PmtProp> & pmtMap = m_DetPropMap[det].PmtPropMap;
    
    int id = sensor.fullPackedData();
    pmtMap[id] = defPmtProp;
    // get label
    ostringstream sring, scolumn;
    sring<<setw(2)<<setfill('0')<<sensor.ring();
    scolumn<<setw(2)<<setfill('0')<<sensor.column();
    string label = sensor.detName()+"_R"+sring.str()+"_C"+scolumn.str();
    //
    string hname = "h_Adc_"+label;
    string hname2 = "h_PreAdc_"+label;
    string hname3 = "h_Tdc_"+label;
    pmtMap[id].h_Adc = new TH1F(hname.c_str(),hname.c_str(),350,-50,300);
    pmtMap[id].h_PreAdc = new TH1F(hname2.c_str(),hname2.c_str(),350,-50,300);
    pmtMap[id].h_Tdc = new TH1F(hname3.c_str(),hname3.c_str(),400,400,1200);
    pmtMap[id].PmtLabel=label;
    pmtMap[id].PmtId=id;
  }
  //  Water Poll PMTs
  for( unsigned int ct=0; ct<poolPmts.size(); ct++ ) {
    PoolPmtSensor sensor = poolPmts[ct];

    Detector det(sensor.site(), sensor.detectorId());
    it = m_DetPropMap.find( det );
    if( it == m_DetPropMap.end() ) {
      m_DetPropMap[det]=defDetProp;
    }
    m_DetPropMap[det].Det = det;

    map<int/*PmtId*/,PmtProp> & pmtMap = m_DetPropMap[det].PmtPropMap;

    int id = sensor.fullPackedData();
    pmtMap[id] = defPmtProp;
    // get label
    ostringstream swall, sspot,sfacing;
    swall<<setw(2)<<setfill('0')<<sensor.wallNumber();
    sspot<<setw(2)<<setfill('0')<<sensor.wallSpot();
    sfacing<<setw(2)<<setfill('0')<<sensor.inwardFacing();
    string label = sensor.detName()+"_W"+swall.str()+"_S"+sspot.str()+"_F"+sfacing.str();
    // 
    string hname = "h_Adc_"+label;
    string hname2 = "h_PreAdc_"+label;
    string hname3 = "h_Tdc_"+label;
    pmtMap[id].h_Adc = new TH1F(hname.c_str(),hname.c_str(),350,-50,300);
    pmtMap[id].h_PreAdc = new TH1F(hname2.c_str(),hname2.c_str(),350,-50,300);
    pmtMap[id].h_Tdc = new TH1F(hname3.c_str(),hname3.c_str(),400,400,1200);
    pmtMap[id].PmtLabel=label;
    pmtMap[id].PmtId=id;
  }
  

  map<Detector, DetProp>::iterator itdet,itdetend=m_DetPropMap.end();
  for(itdet = m_DetPropMap.begin(); itdet!=itdetend; itdet++) {
    Detector det = itdet->first;
    ResetPmtProp(det);//reset
  }

  m_NStop = 1;
  
  return 1;
}
int RollingGain::ResetPmtProp ( Detector  det)

Definition at line 552 of file RollingGain.cc.

{
  DetProp& currDet = m_DetPropMap[det];
  currDet.FirstTrigTime = TimeStamp::GetBOT();
  currDet.LastTrigTime = TimeStamp::GetBOT();
  currDet.FirstRun = -1;
  currDet.LastRun = -1;
  currDet.NTrigger=0;
  currDet.Statistics=0;

  map<int/*PmtId*/,PmtProp> & pmtMap = currDet.PmtPropMap;

  map<int/*PmtId*/,PmtProp>::iterator it,idend=pmtMap.end();
  for(it=pmtMap.begin(); it!=idend; it++) {
    
    PmtProp& curr = it->second;    
    curr.h_Adc->Reset();
    curr.h_PreAdc->Reset();
    curr.h_Tdc->Reset();
    curr.Pedestal=0;
    curr.PedestalErr=0;
    curr.DarkRate=0;
    curr.DarkRateErr=0;
    curr.ElecNoiseRate=0;
    curr.ElecNoiseRateErr=0;
    curr.Gain=0;
    curr.GainErr=0;
    curr.Sigma=0;
    curr.SigmaErr=0;
    curr.Status=PmtCalibData::kUnknown;
    curr.NHits=0;
    curr.FitStatus=-1;
    curr.Chi2=-1;
  }

  return 1;
}

Member Data Documentation

Cable mapping service.

Definition at line 90 of file RollingGain.h.

IStatisticsSvc* RollingGain::m_statsSvc [private]

Definition at line 91 of file RollingGain.h.

Average Pedestal.

Definition at line 93 of file RollingGain.h.

TF1* RollingGain::f1 [private]

Definition at line 95 of file RollingGain.h.

TF1* RollingGain::f2 [private]

Definition at line 96 of file RollingGain.h.

TFile* RollingGain::m_rootfile [private]

Definition at line 97 of file RollingGain.h.

TTree* RollingGain::m_tree [private]

Definition at line 98 of file RollingGain.h.

FILE* RollingGain::m_masterfile [private]

Definition at line 99 of file RollingGain.h.

FILE* RollingGain::m_outfile [private]

Definition at line 100 of file RollingGain.h.

string RollingGain::m_fileName [private]

Definition at line 101 of file RollingGain.h.

Definition at line 139 of file RollingGain.h.

int RollingGain::m_Fix [private]

Definition at line 141 of file RollingGain.h.

int RollingGain::m_output [private]

Definition at line 142 of file RollingGain.h.

int RollingGain::m_NStop [private]

Definition at line 143 of file RollingGain.h.

int RollingGain::m_lastFit [private]

Definition at line 144 of file RollingGain.h.

int RollingGain::m_site [private]

Definition at line 146 of file RollingGain.h.

int RollingGain::m_simFlag [private]

Definition at line 147 of file RollingGain.h.

Definition at line 148 of file RollingGain.h.

int RollingGain::m_thisRun [private]

Definition at line 149 of file RollingGain.h.

double RollingGain::m_fitPeriod [private]

Both m_fitPeriod and m_fitStatistic must be satisfied.

Definition at line 151 of file RollingGain.h.

double RollingGain::m_fitStatistic [private]

Definition at line 152 of file RollingGain.h.

list RollingGain::__all__ = [ "run" ] [private]

Definition at line 3 of file __init__.py.


The documentation for this class was generated from the following files:
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:04:18 for RollingGain by doxygen 1.7.4