LISieve.cxx

Go to the documentation of this file.

// Coded by Jeff Hartnell and Phill Litchfield Feb/2006 onwards
//
// Contact: j.j.hartnell@rl.ac.uk

#include <set>

#include "TClonesArray.h"

#include "MessageService/MsgService.h"
#include "CandNtupleSR/NtpSRStrip.h"
#include "StandardNtuple/NtpStRecord.h"

#include "NtupleUtils/LISieve.h"

CVSID("$Id: LISieve.cxx,v 1.5 2009/11/16 14:02:16 jdejong Exp $");
// Last modified by Jeff de Jong 10 Nov,2009
//......................................................................
//namespace LISieve {
Bool_t LISieve::IsLIforNC(const NtpStRecord& ntp)
{

  //only analyse if FD
  if (ntp.GetHeader().GetVldContext().GetDetector()!=
      Detector::kFar) return false;

  //variable to turn on all the useful messages if required
  Msg::LogLevel_t logLevel=Msg::kVerbose;
    
  //this variable can turn off the noise removal in the event
  const Bool_t cleanTheEvent=true;
  
  //store the set of times for the snarl
  multiset<Double_t> times;
  
  Int_t numStrips=ntp.stp->GetEntries();
  if (numStrips<=0) return false;//pass the event, not LI
  //loop over the strips in the snarl
  for (Int_t hit=0;hit<numStrips;hit++){
    NtpSRStrip* strip = dynamic_cast<NtpSRStrip*>((*ntp.stp)[hit]);

    Double_t time=strip->time0;
    if (time<0 || time>1) time=strip->time1;

    if (time>0 && time<=1) {
      times.insert(time);
    }
    //else just don't put the time in the map - clearly rubbish
  }//end of loop over strips

  //get the median time from the map
  multiset<Double_t>::const_iterator it=times.begin();
  advance(it,times.size()/2);
  Double_t medianTime=*it;
  MAXMSG("LISieve",Msg::kDebug,100)
    <<"LISieve: Median time="<<medianTime<<endl;

  map<Int_t,Int_t> hpp;  // <plane, # of hits>
  Float_t sumSigCorEast=0;
  Float_t sumSigCorWest=0;
  //second loop
  for (Int_t hit=0;hit<numStrips;hit++){
    NtpSRStrip* strip = dynamic_cast<NtpSRStrip*>((*ntp.stp)[hit]);

    if (cleanTheEvent) {
      if (strip->ph0.sigcor > 0){
        if (strip->time0<medianTime-200e-9 || 
            strip->time0>medianTime+200e-9) continue;
      }
      if (strip->ph1.sigcor > 0){
        if (strip->time1<medianTime-200e-9 || 
            strip->time1>medianTime+200e-9) continue;
      }
    } // if (cleanTheEvent)
    
    sumSigCorEast+=strip->ph0.sigcor;
    sumSigCorWest+=strip->ph1.sigcor;

    hpp[strip->plane] += 1;
  }//end of loop over strips

  Float_t avHpp=0;
  vector<Float_t> pb(8);//to store fractions hit in each pb region
  vector<Float_t> pbStrips(8);//to store fractionof strip hit in each pb region

  //loop over all the planes flashed
  for (map<Int_t,Int_t>::iterator hppIt=hpp.begin();
       hppIt!=hpp.end();++hppIt){
    avHpp+=hppIt->second;

    Int_t pl=hppIt->first;

    if (pl>=1 && pl<=64) pb[0]+=1./64;//0=bookend
    else if (pl>=65 && pl<=128) pb[1]+=1./64;
    else if (pl>=129 && pl<=192) pb[2]+=1./64;
    else if (pl>=193 && pl<=248) pb[3]+=1./56;//249=bookend
    else if (pl>=250 && pl<=313) pb[4]+=1./64;
    else if (pl>=314 && pl<=377) pb[5]+=1./64;
    else if (pl>=378 && pl<=441) pb[6]+=1./64;
    else if (pl>=442 && pl<=485) pb[7]+=1./44;


    if      (pl>=1 && pl<=64)    pbStrips[0]+=avHpp*1./(64*192.);//0=bookend
    else if (pl>=65 && pl<=128)  pbStrips[1]+=avHpp*1./(64*192.);
    else if (pl>=129 && pl<=192) pbStrips[2]+=avHpp*1./(64*192.);
    else if (pl>=193 && pl<=248) pbStrips[3]+=avHpp*1./(56*192.);//249=bookend
    else if (pl>=250 && pl<=313) pbStrips[4]+=avHpp*1./(64*192.);
    else if (pl>=314 && pl<=377) pbStrips[5]+=avHpp*1./(64*192.);
    else if (pl>=378 && pl<=441) pbStrips[6]+=avHpp*1./(64*192.);
    else if (pl>=442 && pl<=485) pbStrips[7]+=avHpp*1./(44*192.);



    MAXMSG("LISieve",Msg::kVerbose,2000)
      <<"plane="<<pl<<", hpp="<<hppIt->second<<endl;
  }
  if (hpp.size()>0) avHpp/=hpp.size();
  Float_t sumSig=sumSigCorEast+sumSigCorWest;
  Float_t minSigCorEW = TMath::Min(sumSigCorEast, sumSigCorWest);
  Float_t asym=0;
  
  if (sumSig) asym = (sumSig - 2 * minSigCorEW) / sumSig;

  // Look for pulser box fractions based on planes only

  Float_t fractFlashed1=0;
  Float_t fractFlashed2=0;

  for (vector<Float_t>::iterator it=pb.begin();
       it!=pb.end();++it)
    {
      Float_t fract=*it;
      if (fract>fractFlashed1){
        //copy the 1 into 2
        fractFlashed2=fractFlashed1;
        //then set new highest fractFlashed
        fractFlashed1=fract;
      }
      //make sure to get fractFlashed2 if fractFlashed2<fract<frFla1
      else if (fract>fractFlashed2) {
        fractFlashed2=fract;
      }
      
      MAXMSG("LISieve",Msg::kVerbose,200)
        <<"fract="<<fract<<", f1="<<fractFlashed1
        <<", f2="<<fractFlashed2<<endl;
    }
  
  Float_t ratioFlashed=0;
  if (fractFlashed1>0) ratioFlashed=fractFlashed2/fractFlashed1;
  // --------------------------
  // Same as above but for strip use fraction
  // added Nov 10,2009 by Jeff de Jong
  // Using strips offers abetter way to discriminate
  // LI from signal as LI lights up everything

  Float_t fractStripsFlashed1=0;
  Float_t fractStripsFlashed2=0;

  for (vector<Float_t>::iterator it=pbStrips.begin();
       it!=pbStrips.end();++it)
    {
      Float_t fract=*it;
      if (fract>fractStripsFlashed1){
        //copy the 1 into 2
        fractStripsFlashed2=fractStripsFlashed1;
        //then set new highest fractFlashed
        fractStripsFlashed1=fract;
      }
      //make sure to get fractFlashed2 if fractFlashed2<fract<frFla1
      else if (fract>fractStripsFlashed2) {
        fractStripsFlashed2=fract;
      }
      
      MAXMSG("LISieve",Msg::kVerbose,200)
        <<"fractStrip="<<fract<<", f1="<<fractStripsFlashed1
        <<", f2="<<fractStripsFlashed2<<endl;
    }
  
  Float_t ratioStripsFlashed=0;
  if (fractStripsFlashed1>0) 
    ratioStripsFlashed=fractStripsFlashed2/fractStripsFlashed1;

  //-------------------------


  Bool_t li=false;

  // check if LI
  // don't check the asymmetry if the VA is saturating
//  if ( (asym>0.5 || minSigCorEW>1.7e6) && 
//       avHpp>3 && ratioFlashed<0.05 &&
//       fractFlashed1>0.85 && fractFlashed2<0.1 )
// New cuts added on 10 Nov,2009 by J. de Jong
  if ( (asym>0.55 || minSigCorEW>1.7e6) && 
       ratioStripsFlashed<0.006 &&
       fractStripsFlashed1>0.02 )
    {
      li=true;
      
      MAXMSG("LISieve",logLevel,100)
        <<"LI Event:"<<endl
        <<"  Asymmetry="<<asym<<endl
        <<"  Av. HPP  ="<<avHpp<<endl
        <<"  Fract Fl1="<<fractFlashed1
        <<", Fract Fl2="<<fractFlashed2
        <<"  Ratio Fla="<<ratioFlashed
        <<"  Strips Fl1="<<fractStripsFlashed1
        <<", StripsFl2="<<fractStripsFlashed2
        <<"  Ratio Strips="<<ratioStripsFlashed<<endl;
    }
  else
    {
      MAXMSG("LISieve",logLevel,10)
        <<"Non-LI Event:"<<endl
        <<"  Asymmetry="<<asym<<endl
        <<"  Av. HPP  ="<<avHpp<<endl
        <<"  Fract Fl1="<<fractFlashed1
        <<", Fract Fl2="<<fractFlashed2
        <<"  Ratio Fla="<<ratioFlashed
        <<"  Strips Fl1="<<fractStripsFlashed1
        <<", StripsFl2="<<fractStripsFlashed2
        <<"  Ratio Strips="<<ratioStripsFlashed<<endl;
    }
  
  return li;
}
//}
Bool_t LISieve::IsLI(const NtpStRecord& ntp)
{

  //only analyse if FD
  if (ntp.GetHeader().GetVldContext().GetDetector()!=
      Detector::kFar) return false;

  //variable to turn on all the useful messages if required
  Msg::LogLevel_t logLevel=Msg::kVerbose;
    
  //this variable can turn off the noise removal in the event
  const Bool_t cleanTheEvent=true;
  
  //store the set of times for the snarl
  multiset<Double_t> times;
  
  Int_t numStrips=ntp.stp->GetEntries();
  if (numStrips<=0) return false;//pass the event, not LI
  //loop over the strips in the snarl
  for (Int_t hit=0;hit<numStrips;hit++){
    NtpSRStrip* strip = dynamic_cast<NtpSRStrip*>((*ntp.stp)[hit]);

    Double_t time=strip->time0;
    if (time<0 || time>1) time=strip->time1;

    if (time>0 && time<=1) {
      times.insert(time);
    }
    //else just don't put the time in the map - clearly rubbish
  }//end of loop over strips

  //get the median time from the map
  multiset<Double_t>::const_iterator it=times.begin();
  advance(it,times.size()/2);
  Double_t medianTime=*it;
  MAXMSG("LISieve",Msg::kDebug,100)
    <<"LISieve: Median time="<<medianTime<<endl;

  map<Int_t,Int_t> hpp;  // <plane, # of hits>
  Float_t sumSigCorEast=0;
  Float_t sumSigCorWest=0;
  //second loop
  for (Int_t hit=0;hit<numStrips;hit++){
    NtpSRStrip* strip = dynamic_cast<NtpSRStrip*>((*ntp.stp)[hit]);

    if (cleanTheEvent) {
      if (strip->ph0.sigcor > 0){
        if (strip->time0<medianTime-200e-9 || 
            strip->time0>medianTime+200e-9) continue;
      }
      if (strip->ph1.sigcor > 0){
        if (strip->time1<medianTime-200e-9 || 
            strip->time1>medianTime+200e-9) continue;
      }
    } // if (cleanTheEvent)
    
    sumSigCorEast+=strip->ph0.sigcor;
    sumSigCorWest+=strip->ph1.sigcor;

    hpp[strip->plane] += 1;
  }//end of loop over strips

  Float_t avHpp=0;
  vector<Float_t> pb(8);//to store fractions hit in each pb region

  //loop over all the planes flashed
  for (map<Int_t,Int_t>::iterator hppIt=hpp.begin();
       hppIt!=hpp.end();++hppIt){
    avHpp+=hppIt->second;

    Int_t pl=hppIt->first;

    if (pl>=1 && pl<=64) pb[0]+=1./64;//0=bookend
    else if (pl>=65 && pl<=128) pb[1]+=1./64;
    else if (pl>=129 && pl<=192) pb[2]+=1./64;
    else if (pl>=193 && pl<=248) pb[3]+=1./56;//249=bookend
    else if (pl>=250 && pl<=313) pb[4]+=1./64;
    else if (pl>=314 && pl<=377) pb[5]+=1./64;
    else if (pl>=378 && pl<=441) pb[6]+=1./64;
    else if (pl>=442 && pl<=485) pb[7]+=1./44;

    MAXMSG("LISieve",Msg::kVerbose,2000)
      <<"plane="<<pl<<", hpp="<<hppIt->second<<endl;
  }
  if (hpp.size()>0) avHpp/=hpp.size();
  Float_t sumSig=sumSigCorEast+sumSigCorWest;
  Float_t minSigCorEW = TMath::Min(sumSigCorEast, sumSigCorWest);
  Float_t asym=0;
  
  if (sumSig) asym = (sumSig - 2 * minSigCorEW) / sumSig;

  // Look for pulser box fractions based on planes only

  Float_t fractFlashed1=0;
  Float_t fractFlashed2=0;

  for (vector<Float_t>::iterator it=pb.begin();
       it!=pb.end();++it)
    {
      Float_t fract=*it;
      if (fract>fractFlashed1){
        //copy the 1 into 2
        fractFlashed2=fractFlashed1;
        //then set new highest fractFlashed
        fractFlashed1=fract;
      }
      //make sure to get fractFlashed2 if fractFlashed2<fract<frFla1
      else if (fract>fractFlashed2) {
        fractFlashed2=fract;
      }
      
      MAXMSG("LISieve",Msg::kVerbose,200)
        <<"fract="<<fract<<", f1="<<fractFlashed1
        <<", f2="<<fractFlashed2<<endl;
    }
  
  Float_t ratioFlashed=0;
  if (fractFlashed1>0) ratioFlashed=fractFlashed2/fractFlashed1;


  Bool_t li=false;

  // check if LI
  // don't check the asymmetry if the VA is saturating
  if ( (asym>0.5 || minSigCorEW>1.7e6) && 
       avHpp>3 && ratioFlashed<0.05 &&
       fractFlashed1>0.85 && fractFlashed2<0.1 )
    {
      li=true;
      
      MAXMSG("LISieve",logLevel,100)
        <<"LI Event:"<<endl
        <<"  Asymmetry="<<asym<<endl
        <<"  Av. HPP  ="<<avHpp<<endl
        <<"  Fract Fl1="<<fractFlashed1
        <<", Fract Fl2="<<fractFlashed2
        <<"  Ratio Fla="<<ratioFlashed << endl;
    }
  else
    {
      MAXMSG("LISieve",logLevel,10)
        <<"Non-LI Event:"<<endl
        <<"  Asymmetry="<<asym<<endl
        <<"  Av. HPP  ="<<avHpp<<endl
        <<"  Fract Fl1="<<fractFlashed1
        <<", Fract Fl2="<<fractFlashed2
        <<"  Ratio Fla="<<ratioFlashed <<endl;
    }
  
  return li;
}


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

---