Blinder.cxx

Go to the documentation of this file.

// $Id: Blinder.cxx,v 1.12 2007/11/11 07:29:53 rhatcher Exp $
//
// Module to apply blindness for physics analysis.
//
// Blinder::Ana returns kPassed for events that should be in the open dataset
// and kFailed for events that should be only in the blind dataset.
//
// It returns kPassed for all events in the ND, and uses
// the algorithm presented at the Jan 2005 collab meeting to
// select events from the FD.
//
// Nathaniel Tagg
// n.tagg1@physics.ox.ac.uk
#include "Filtration/Blinder.h"
#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "JobControl/JobCModuleRegistry.h" // For JOBMODULE macro
#include "RawData/RawRecord.h"
#include "RawData/RawDaqSnarlHeader.h"
#include "RawData/RawDigitDataBlock.h"
#include "RawData/RawDigit.h"
#include "Plex/PlexHandle.h"

#include "Record/SimSnarlRecord.h"
#include "REROOT_Classes/REROOT_NeuKin.h"

#include <TFile.h>
#include <TNtuple.h>
#include <TTree.h>
#include <TRandom3.h>
#include <TClonesArray.h>

#include <vector>
#include <cassert>
#include <cmath>

JOBMODULE(Blinder, "Blinder",
          "Blinds events away");


CVSID(" $Id: Blinder.cxx,v 1.12 2007/11/11 07:29:53 rhatcher Exp $ ");


Blinder::BlindingFunction::BlindingFunction()
{
  fLFreq  = 0;
  fEFreq  = 0;
  fLPhase = 0;
  fEPhase = 0;
}

void
Blinder::BlindingFunction::SetParameters( TRandom* paramChooser, Int_t seed )
{
  assert(paramChooser);
  paramChooser->SetSeed(seed);
  
  fLFreq  = paramChooser->Uniform(1.0/15.0, 1.0/500.0); // Lphase. 5-200 planes.
  fEFreq  = paramChooser->Uniform(1.0/0.8, 1.0/2.0) / 10000; // efreq. 10000 ~= 1GeV
  fLPhase = paramChooser->Uniform(0,3.14); // phase
  fEPhase = 3.14-fLPhase; 
}

Float_t 
Blinder::BlindingFunction::EventBlindProb( Float_t length, 
                                          Float_t energy )
{
  float b1 = 0.5*(sin(length*fLFreq + fLPhase)+1.);
  float b2 = 0.5*(sin(energy*fEFreq + fEPhase)+1.);
  float f = 0.25 + b1*b2;
  if(f>1.0) f=1.0;
  return f; 

}


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

Blinder::Blinder() 
{
  fSeed = 0;
  fFile = 0;
  fTree = 0;

  fRandom = new TRandom3;
  fParameterChooser = new TRandom3;
  fBlindingFunction.SetParameters(fParameterChooser,fSeed);
}

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

Blinder::~Blinder()
{
  if(fFile) {
    fFile->cd();
    fTree->Write();
    fFile->Write();
    fFile->Close();
  }
  if(fRandom) delete fRandom;
  if(fParameterChooser) delete fParameterChooser;
}

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

JobCResult Blinder::Ana(const MomNavigator* mom)
{
  Int_t run;
  Int_t snarl;
  Int_t epoch;
  Float_t length;
  Float_t energy;


  if(GetInfo(mom,run,snarl,epoch,length,energy)) {
    // Set the parameters from the epoch:
    fBlindingFunction.SetParameters(fParameterChooser,fSeed + 3 * epoch);
    
    // Find the probability of blinding this event:
    Float_t prob = fBlindingFunction.EventBlindProb(length,energy);

    // Get a pseudo-random number specific for this event:
    fRandom->SetSeed(run*1000 + snarl);
    Float_t r = fRandom->Rndm();
    
    // Allow event to be open if r<prob
    if(r<=prob) return JobCResult::kPassed;
    else        return JobCResult::kFailed;
  };

  // No blindness info obtained; return default. (This is probably a monitoring block.)
  return JobCResult::kPassed;
}


Bool_t Blinder::GetInfo(const MomNavigator* mom,
                        Int_t &run,
                        Int_t &snarl,
                        Int_t &epoch,
                        Float_t &length,
                        Float_t &energy)
{
  //
  // This function pulls the neccessary info from the current Mom to find
  // the run, snarl, event length (with a basic clustering algorithm) and 
  // event energy (in raw ADC counts).
  //
  run = 0; 
  snarl = 0;
  length = 0;
  energy = 0;
  epoch = 0;
  
  // Find raw data.
  const RawRecord* rr = dynamic_cast<const RawRecord*>(mom->GetFragment("RawRecord"));
  if(rr==0) return false;

  const RawDaqSnarlHeader* snarlHead= dynamic_cast<const RawDaqSnarlHeader*>(rr->GetRawHeader());
  if(snarlHead==0)  return false;


  TIter rdbit = rr->GetRawBlockIter();
  TObject *tob;
  const RawDigitDataBlock *rddb = 0;
  while ((tob = rdbit())) {
    if (tob->InheritsFrom("RawDigitDataBlock")) {
      rddb = (const RawDigitDataBlock *) tob;
    }
  }
  if(rddb==0) return false;

  if(snarlHead->GetVldContext().GetDetector()!=Detector::kFar) return false;

  if(snarlHead->GetVldContext().GetSimFlag()!=SimFlag::kData) epoch = 10;

  // Put epoch here:
  // These lines change the blinding function at certain
  // watershed dates.
  const VldTimeStamp kShutdown2006(2006,4,15,0,0,0); // April 15, 2006, midnight
  if(snarlHead->GetVldContext().GetTimeStamp() > kShutdown2006) epoch +=1;

  MAXMSG("Blinder",Msg::kInfo,1) << "Blinding FD data. Epoch = " << epoch << endl;
  
  // Get plex info.
  PlexHandle plex(snarlHead->GetVldContext());

  // Get the header data.
  run = snarlHead->GetRun();
  snarl = snarlHead->GetSnarl();
 
  // Map of the planes.
  std::vector<int> planes(500,0);

  // Loop the digits.
  TIter rdit = rddb->GetDatumIter();
  const RawDigit *digit;
  Int_t rawdigitindex = -1;
  while ((digit = (const RawDigit *) rdit())) {
    ++rawdigitindex;
    
    if(plex.GetReadoutType(digit->GetChannel()) != ReadoutType::kScintStrip) 
      continue;

    // Find plane.
    PlexSEIdAltL altl = plex.GetSEIdAltL(digit->GetChannel());
    Int_t plane = altl.GetPlane();

    // Find energy.
    Int_t adc = digit->GetADC();

    // add it up.
    if((plane>=0)&&(plane<490))
      planes[plane] ++;

    energy+=(float)adc;
  }

  // do the number-of-planes matching.  Look for the largest 'cluster' 
  int start = 0;
  int end = 0;
  int longest = 0;
  while(start<500) {
    if(planes[start]==0) {
      start++;
      continue;
    };
    
    // Found a starting place.
    end = start;
    int gap = 0;
    int thislen = 0;
    while(end<490) {
      end++;
      if(planes[end]==0) gap++;
      else thislen = end-start;
      if(gap>2) break;
    }
    
    if(thislen > longest) longest = thislen;

    start=end;    
  }
  length = longest;

  return true;
}


//......................................................................
// The stuff below is used only for studies of the blindness itself, not
// for the actual blinding routine.
//......................................................................

JobCResult Blinder::Reco(MomNavigator* mom)
{
  Int_t run;
  Int_t snarl;
  Int_t epoch;
  Float_t length;
  Float_t energy;

  if(GetInfo(mom,run,epoch,snarl,length,energy)) {
    float prob = fBlindingFunction.EventBlindProb(length,energy);

    fRandom->SetSeed(run*1000 + snarl);
    Float_t r = fRandom->Rndm();

    float blind = 0;
    if(r<=prob) blind = 1;

    float enu=0;
    float evis=0;
    float cc=0;
    float emu=0;
    float esh=0;
    float nu=0;

    const REROOT_NeuKin* truth = GetTruth(mom);
    if(truth) {
      enu = fabs(truth->P4Neu()[3]);
      cc  = truth->IAction();
      emu = fabs(truth->P4Mu1()[3]) + fabs(truth->P4Mu2()[3]);
      esh = fabs(truth->P4Shw()[3]);
      evis = emu + esh 
        + fabs(truth->P4El1()[3])
        + fabs(truth->P4El2()[3])
        + fabs(truth->P4Tau()[3]);      
      nu = truth->INu();
    }

    if(!fTree) {
      fFile = new TFile("blind.root","RECREATE");
      fTree = new TNtuple("blind","blind",
                          "run:snarl:len:energy"
                          ":enu:evis:cc:nu:emu:esh"
                          ":prob:blind"
                          );      
    }

    if(fTree) {
      fTree->Fill(run,snarl,length,energy
                  ,enu,evis,cc,nu,emu,esh
                  ,prob,blind
                  );      
    }
  };
  
  return JobCResult::kPassed; // kNoDecision, kFailed, etc.
}

const REROOT_NeuKin* Blinder::GetTruth(const MomNavigator* mom)
{
    // Look for the simSnarl in Mom.
  SimSnarlRecord* simsnarl = 0;
  TObject* tobj;
  TIter    fragiter = mom->FragmentIter();
  while( ( tobj = fragiter.Next() ) ) {
    simsnarl = dynamic_cast<SimSnarlRecord*>(tobj);
    if(simsnarl) break;
  }

  if(simsnarl==0) return 0;

  // Get the NeuKin list.
  const TClonesArray* neuKins = dynamic_cast<const TClonesArray*>(simsnarl->FindComponent("TClonesArray","NeuKinList"));
  if(neuKins==0) return 0;

  return dynamic_cast<const REROOT_NeuKin*>((*neuKins)[0]); 
}




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

const Registry& Blinder::DefaultConfig() const
{
//======================================================================
// Supply the default configuration for the module
//======================================================================
  static Registry r; // Default configuration for module

  // Set name of config
  std::string name = this->GetName();
  name += ".config.default";
  r.SetName(name.c_str());

  // Set values in configuration
  r.UnLockValues();
  r.Set("Seed", 10127);  // Seed with Nathaniel's fnal ID number.. that's pretty random.
  r.LockValues();

  return r;
}

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

void Blinder::Config(const Registry& r)
{
  //======================================================================
  // Configure the module given the Registry r
  //======================================================================
  
  r.Get("Seed",fSeed);
  
  // reset the paramaters:
  fBlindingFunction.SetParameters(fParameterChooser,fSeed);
}

---