FillShortEvent.cxx

Go to the documentation of this file.

// $Id: FillShortEvent.cxx,v 1.3 2010/01/05 23:35:18 jyuko Exp $

// C/C++
#include <cmath>
#include <cassert>
#include <iostream>
#include <sstream>

// ROOT
#include "TDirectory.h"
#include "TH1.h"
#include "TH2.h"

// MINOS
#include "Registry/Registry.h"

// Local
#include "PhysicsNtuple/Default.h"
#include "PhysicsNtuple/Factory.h"
#include "PhysicsNtuple/Hist/HistMan.h"
#include "PhysicsNtuple/Record.h"

// Local
#include "FillShortEvent.h"

REGISTER_ANP_OBJECT(AlgSnarl,FillShortEvent)

using namespace std;

//----------------------------------------------------------------------
Anp::FillShortEvent::FillShortEvent()
  :fRecoHighE(3),
   fNEndPlanes(5),
   fMinADCHit(200),
   fNPlaneSkip(4),
   fTimeWindow((2.0*(18.83 + 0.1)*1.0e-9)),
   fStripWindow(4),
   fKeyBase(19000),
   fMinDLength(4),
   fCutOff(0.5),
   fKeyPass(-1),
   fKeyPass2(-1),
   fDirName("shortevent"),
   fDebug(false),
   fErase(false)
{
}

//----------------------------------------------------------------------
Anp::FillShortEvent::~FillShortEvent()
{
}
//----------------------------------------------------------------------
bool Anp::FillShortEvent::Run(Record &record)
{ 


  EventIterator ievent = record.EventBegIterator();
  while(ievent != record.EventEndIterator())
  { 
    //
    // Fill truth variable used by training algorithm(s)
    //

    if(!record.GetHeader().IsData() && record.TruthBeg() != record.TruthEnd())
      {
        const TruthIter truth = record.FindTruth(*ievent);
        if(truth == record.TruthEnd()) return true;
          
        double iscc= 0;
        if(truth-> IsCC()) iscc = 1;
            
        if(!(ievent->Add(fKeyBase, iscc)))
          {
            cerr << "FillShortEvent::Run - failed to insert value at " << fKeyBase+300 << " key" << endl;
          }
      }

    TrackIter track_it = LongestTrack(*ievent,record);
    if(track_it==record.TrackEnd()||
       (fKeyPass>0 && fKeyPass2>0 &&!track_it->KeyExists(fKeyPass) &&   !track_it->KeyExists(fKeyPass2)))
      {
      if(!FillShortEvent::Explore(*ievent, record) && fErase)
       {
        ievent = record.Erase(ievent);
        continue;
      }
    }

    ievent++;
  }
  if(fErase)  Anp::CleanRecord(record.EventBeg(), record.EventEnd(), record);
  return true;
}
//---------------------------------------------------------------------
void Anp::FillShortEvent::Config(const Registry &reg)
{
  reg.Get("FillShortEventMaxRecoE",  fRecoHighE);
  reg.Get("FillShortEventMinNPlane",  fNEndPlanes);
  reg.Get("FillShortEventMinADCHit",  fMinADCHit);
  reg.Get("FillShortEventKeyBase",    fKeyBase);
  reg.Get("FillShortEventTimeWindow", fTimeWindow);
  reg.Get("FillShortEventStripWindow", fStripWindow);
  reg.Get("FillShortEventNPlaneSkip", fNPlaneSkip);
  reg.Get("FillShortEventMinBVDLength",  fMinDLength);
  reg.Get("FillShortEventCutOff", fCutOff);
  reg.Get("FillShortEventKeyPass", fKeyPass);
  reg.Get("FillShortEventKeyPass2", fKeyPass2);

  //
  // Read bool and string variables
  //
  Anp::Read(reg, "FillShortEventDebug",    fDebug);
  Anp::Read(reg, "FillShortEventErase",    fErase);


  if(reg.KeyExists("PrintConfig"))
   {
     cout << "FillShortEvent::Config" << endl;
     cout << "KeyBase = " << fKeyBase<< endl;
     cout << "BDV Length = " << fMinDLength<< endl;
     cout << "TimeWindow = "<< fTimeWindow<<endl;
     cout << "KeyPass = "<<fKeyPass<<endl;
     cout << "KeyPass2 = "<<fKeyPass2<<endl;
     cout<<  "Erase ="<< fErase<<endl;
   }
}
//--------------------------------------------------------------------
bool Anp::FillShortEvent::Init(const Header& ) 
{
  fTooShort=0;
  fBVDmiss =0;
  return true;
}
//--------------------------------------------------------------------                                                                               
void Anp::FillShortEvent::End(const DataBlock &)
{
  cout<<" Nmissed (too short)      "<< fTooShort<<endl;
  cout<<" Nmissed (bvd too short)  "<< fBVDmiss<<endl;
}
//-----------------------------------------------------------------------
bool Anp::FillShortEvent::Explore(Event & event, const Record& record){

  PlaneMap mapu;
  PlaneMap mapv;
  const Basic eventb= event.GetBasic();

  if(eventb.NUPlane()<fNEndPlanes &&eventb.NVPlane()<fNEndPlanes){
    fTooShort++;
    return false;
  }
  if(!FillShortEvent::GetPlaneStrips(mapu, mapv, event, record))   return false;
 
  std::map<int, double> ansv;  std::map<int, double> ansu;

  ansu[0]=(double)mapu.size();
  ansv[0]=(double)mapv.size();
  
  bool ugood = false; 
  bool vgood = false;
  if(fDebug) cout<<" FillShortEvent::Explore - find event quantities "<<endl;
  if(mapu.size()>=fNEndPlanes  &&FillShortEvent::GetEventQuantities(mapu,&ansu))    ugood=true;
  if(mapv.size()>=fNEndPlanes && FillShortEvent::GetEventQuantities(mapv,&ansv))    vgood=true;
  
  std::map< int, double> ans;

  if(!(ugood) && !(vgood) )
  {
    fBVDmiss ++;
    return false;
  }
  else if(!ugood)    ans = ansv;
  else if(!vgood)    ans = ansu;
  else if(ansu[3] >= ansv[3]) ans = ansu;
  else ans = ansv;
  if(fDebug) cout<<" choose line with size "<<ans[3]<<endl;

  ans[12]= event.Gev();

  for(unsigned int ii=0; ii<ans.size(); ii++)    event.Add(fKeyBase+ii+1, ans[ii]);

  return true;
}

//---------------------------------------------------------------------
bool Anp::FillShortEvent::GetPlaneStrips(PlaneMap& plane_mapu, PlaneMap& plane_mapv, const Event ievent, const Record& record)
{
  for(StripIter istrip = record.StripBeg(); istrip != record.StripEnd(); ++istrip)
    {
      if(!istrip -> MatchEvt(ievent.EventIndex()))
        {
          continue;
        }

      if(record.GetHeader().IsNear() && istrip -> GetPlane() > 120)
        {
          return false;
        }

      if(istrip -> SigCor() < fMinADCHit)
        {
          continue;
        }

      if(istrip->IsUview()){
        vector<StripIter> &svec = plane_mapu[istrip -> GetPlane()];
        svec.push_back(istrip);
      }
      else if(istrip->IsVview()){
        vector<StripIter> &svec = plane_mapv[istrip -> GetPlane()];
        svec.push_back(istrip);
      }
    }
  return true;
}
//--------------------------------------------------------------------------
bool Anp::FillShortEvent::GetEventQuantities(PlaneMap plane_m, std::map<int,double>* vect)
{
  //
  // Get Some basic information (# hits, # planes)
  //
  double NHits=0.0; //(number of hits)                                
  int count=0;      // count of planes                                                                
  double maxT=-10000;
  double minT= 10000;

  PlaneMap::iterator last_pit = plane_m.begin();

  for (PlaneMap::iterator pit = plane_m.begin(); pit!=plane_m.end();++pit){
     vector<StripIter> &svec  = pit -> second;
    if(count>2)
      {
      NHits+=svec.size();
  
      // count hits for all but first two planes (can be noisy)                   
      for(vector<StripIter>::iterator it = svec.begin(); it!=svec.end(); it++)
      {    
        if( (*it)->TPos()> maxT ) maxT = (*it)->TPos();
        else if((*it)->TPos()<minT) minT= (*it)->TPos();
      }
    }
    count++;
  }

  if(count<2 || NHits<=0)
  {
    if(fDebug) cout<<" FillShortEvent:: GetEventQuantities - return false for too small vector"<<endl;
    return false; 
  }

 //
 // Get the BVD best line 
 //

 PlaneMap dv_map =plane_m;
 double dijkstra = BVD_BestLine(dv_map);
 double size = dv_map.size();
 if(dijkstra<fCutOff ||size<0)
 {
   if(fDebug) cout<<" dijkstra best line not found. Size = "<<size<<endl;
   return false;
 }

 if(count>2*size) cout<<" Number of planes is much greater than found path "<< count<<" "<<size<<endl;

 if(size>=fMinDLength){

   double sumph =-1; double sum_near_ph=-1;
   double sigph =-1;
   double lastph=-1;
   vector<double> time;
   
   int c=0;
   //
   // first loop
   //
   for(PlaneMap::reverse_iterator si = dv_map.rbegin(); si!=dv_map.rend(); si++){
     if(sumph==-1){
       sumph=0;
       lastph=0;
     }
     const vector<StripIter> &svec  = si -> second;
     sumph += (svec[0])->SigCor();
     time.push_back(svec[0]->Time());
     if(c<3)lastph+=(svec[0])->SigCor();
     c++;
   }
   
   std::sort(time.begin(), time.end());
   //
   // second loop
   //
   int nearstrips=0;
   for(PlaneMap::iterator si = dv_map.begin(); si!=dv_map.end(); si++){
     if(sum_near_ph<0) sum_near_ph=0;
     vector<StripIter> &svec  = si -> second;
     int minstrip = svec[0]->GetStrip();
     int maxstrip = svec[0]->GetStrip();
     for( vector<StripIter>::iterator its = svec.begin(); its!=svec.end();its++){
       if((*its)->GetStrip()< minstrip) minstrip = (*its)->GetStrip();
       if((*its)->GetStrip()> maxstrip) maxstrip = (*its)->GetStrip();
     }
     vector<StripIter> pmit_vec = plane_m[si->first];
     for(vector<StripIter>::iterator it = pmit_vec.begin(); it!=pmit_vec.end(); it++){
       int its=(*it)->GetStrip();
       
       if( ( its> minstrip-fStripWindow && its<maxstrip+fStripWindow)
           && (*it)->Time()< time.back()+fTimeWindow && (*it)->Time()>time.front()-fTimeWindow){
         sum_near_ph+= (*it)->SigCor(); 
         nearstrips++;
       }
     }
     
     if(sigph==-1) sigph =0;
     sigph +=  ((svec[0])->SigCor() -(sumph/size))*((svec[0])->SigCor() -(sumph/size))/size;
   }
   vect->insert( pair<int,double> (1,(double)NHits/(count-2)));  // 1
   vect->insert( pair<int,double> (2,(maxT-minT)/(count-2)));     // 2
   vect->insert( pair<int,double> (3,size));                     // 3  
   vect->insert( pair<int,double> (4,dijkstra));                 // 4    
   vect->insert( pair<int,double> (5,sumph/size));               // 5    
   vect->insert( pair<int,double> (6,sqrt(sigph)));              // 6
   vect->insert( pair<int,double> (7,lastph));                   // 7
   vect->insert( pair<int,double> (8,(sumph+sum_near_ph)/size)); // 8
   vect->insert( pair<int,double> (9,size/count));               // 9
   vect->insert( pair<int,double> (10,(time.back()-time.front())
                                   /(1.0e-9*time.size())));      // 10
   vect->insert( pair<int,double> (11,nearstrips/size));         // 11
  
   return true;
 }

 if(fDebug) cout<<"FillShortEvent::GetEventQuantities  Return false because size is less than minimum length "<<size<<endl;
 return false;

}
//-----------------------------------------------------------------------------------
double Anp::FillShortEvent::BVD_BestLine(PlaneMap& map)
{
  if(fDebug) cout<<" FillShortEvent:: BVD_BestLine---------------------------------------------"<<endl;
  //,,,,,,,,,,,,,,,,,,,,,,,                                                                                                                     
  // Initialize the vector                                                                                               
  //'''''''''''''''''''''''         
  if(map.size() ==0 )return -1;
  PlaneMapList map_list;
  for(PlaneMap::iterator mapit = map.begin(); mapit!=map.end();mapit++)
  {
    vector<Anp::FillShortEvent::StripLinkPtr> map_vec;
    for(vector<StripIter>::iterator it =(mapit->second).begin(); it!=(mapit->second).end();it++)
      {
        StripLink* str_link_it= new StripLink();
        if(!FillLink(*str_link_it, 0,0,0,(*it))){
          if(fDebug) cout<<" StripLink did not fill"<<endl;
          continue;
        }
        map_vec.push_back(str_link_it);
      }
    // For each plane get a vector of Strip 'Links'
    map_list[mapit->first] =map_vec;       
  }

  //......................
  // Sift through and link
  //''''''''''''''''''''''

  vector<Anp::FillShortEvent::StripLinkPtr> last_vec = map_list.begin()->second;

  // Loop over each plane

  for(PlaneMapList::reverse_iterator list_it = map_list.rbegin(); list_it!=map_list.rend(); list_it++){
    if((list_it == map_list.rbegin()))
    {
      last_vec = list_it->second;
      continue;
    }

    // loop over the vector of strips for this plane
    for(vector<StripLinkPtr>::iterator iter = (list_it->second).begin(); iter!=(list_it->second).end(); iter++)
      {
       
        StripLink*  striplink  = (*iter);
        striplink->ext_value=0.0;
        // loop over the vector of strips for the last plane 
        for(vector<StripLinkPtr>::iterator lp_it = last_vec.begin(); lp_it!=last_vec.end(); lp_it++)
          {

            // go through the linked list of all previous strips connected to this particlar strip
            // if any get a higher probability of being 'connected',despite the skipped plane this becomes the new strip to beat
            // This is now the prbability for "striplink" to be connected to "link"         

            StripLink*  link  = (*lp_it); // first strip in the chain
            double linkp = GetProbability(striplink, link);

            StripLink* linkn= link->laststriplink; // strip connected to this strip
            int skipper=0;
            while(linkn!=NULL)
            {                                                                                            
              if(skipper> 2) break;//fNPlaneSkip?
              double c_prob= GetProbability(striplink, linkn);
         
              if(c_prob > linkp)
              {
                link = linkn;
                linkp = c_prob;
              }
      
              linkn= linkn->laststriplink;            
              skipper++;
              if( c_prob< 0)
              {
                if(fDebug) cout<<" probabiltiy is less than zero"<<endl;
                break;
              }
            }
            if(fDebug) cout<<" FillShortEvent::BVD goodlength has  probability  : "<< linkp<<endl;

            double dtpos =((striplink)->strip)->TPos()-((link)->strip)->TPos();
            double dzpos= ((striplink)->strip)->ZPos()-((link)->strip)->ZPos();
            double ddx_value =0;
            if(link->link_number>2) ddx_value =((link->ddx_value)*(link->link_number-3) +((dtpos/dzpos- link->ins_value)/dzpos))/(link->link_number-2);

            if(dzpos==0)
            {
                if(fDebug) cout<<" dzpos is zero"<<endl;
                continue;
            }
            else if(linkp> striplink->ext_value)// If you have a larger probability
            {
              FillLink(*striplink, dtpos/dzpos, linkp, ddx_value, striplink->strip, link);           
            }
            else if(linkp== striplink->ext_value &&  (striplink->laststriplink==NULL  ||
                                                      sqrt(dtpos*dtpos +dzpos*dzpos) < sqrt(pow(striplink->strip->TPos() - striplink->laststriplink->strip->TPos(),2) + 
                                                                                            pow(striplink->strip->ZPos()- striplink->laststriplink->strip->ZPos(),2))))
            {                 
              FillLink(*striplink, dtpos/dzpos, linkp, ddx_value, striplink->strip, link);            
            }
          }     
        if( striplink->laststriplink!=NULL) ((striplink)->laststriplink)->end_link = false;     
      }
    last_vec = list_it->second;
  }
  

  int cnt=0;
  double sumpercnt =0.0;
  StripLinkPtr slink =NULL;

    for(PlaneMapList::reverse_iterator list_it = map_list.rbegin(); list_it!=map_list.rend(); list_it++){
      vector<Anp::FillShortEvent::StripLinkPtr> link_vec = list_it->second;
      cnt++;
      if(cnt<fMinDLength)       continue; 

      for(vector<StripLinkPtr>::iterator iter = link_vec.begin(); iter!=link_vec.end(); iter++)  {
        int link_n =(int) (*iter)->link_number;
        if(fDebug) cout<<" -- "<< link_n<<endl;
        if(link_n+fNPlaneSkip< map.size()) continue;
          if(pow((*iter)->ext_value,1.0/((double)link_n))*pow((double)link_n,2)  > sumpercnt &&link_n >=fMinDLength )  {         
          slink= (*iter);
          sumpercnt = pow((*iter)->ext_value,1.0/((double)link_n))*pow((double)link_n,2);           
        }
      }      

    }

    double size=0;
    if(slink!=NULL) size = slink->link_number;
  PlaneMap mapout;
  while(slink!=NULL){
    int plane = (slink->strip)->GetPlane();
    mapout[plane].push_back(slink->strip);    
    slink= slink->laststriplink;
  }
  if(!(size==mapout.size()) ) cout<<" FillShortEvent::BVD_bestline sizes don't match"<< mapout.size()<<" "<<size<<endl;
  if(fDebug) cout<<" FillShortEvent::BVD_bestline End -------   found line of length "<< mapout.size()<<" from "<<map.size() <<" planes ("<<size<<")"<< endl;

  map = mapout;
  if(mapout.size()==0) return 0;
  return  sumpercnt/pow((double)size,2.0);
  }


//-----------------------------------------------------------------------------------
double Anp::FillShortEvent::GetProbability(StripLink* cur, StripLink* cmp){
  
  //
  // Penalty for missing a plane
  //
  double penalty = 0.728;
  // 
  // Get the derivative - change in position
  //
  double dtpos =((cur)->strip)->TPos()-((cmp)->strip)->TPos();
  double dzpos= ((cur)->strip)->ZPos()-((cmp)->strip)->ZPos();

  //
  // Get the last extremum value, instantaneous (derivative) value
  //
  double lsum =(cmp)->ext_value;
  double lins =(cmp)->ins_value;
  //
  // Get the number of skipped planes
  //
  double pndist = (cmp->strip->GetPlane() - cur->strip->GetPlane())/2 -1;
  if(pndist<0 || dzpos>=0){
    if(fDebug) cout<<" Number of skipped planes is negative! "<<pndist<<" "<<dzpos<<endl;
    return -1;
  }

  double pred_dtdz =0; //predicted change in step       
  double dtpred =0;
  double costheta=0;
  //
  // Get cosine (delta theta)  
  //where delta theta is the angle between pervious derivative and current
  if((cmp)->laststriplink!=NULL)
  {
    if(fDebug) cout<<" cmp->laststriplink ! = NULL "<<endl;
    // not comparing to first link
    if(cmp->laststriplink->laststriplink!=NULL){
      pred_dtdz=(cmp)->ddx_value;
      if(pred_dtdz*dzpos > 10){
        //      cout<<" Very large change in pred_dtdz"<< pred_dtdz<<endl;
        pred_dtdz =0.0;
      }
    }
      dtpred = (lins + pred_dtdz*dzpos)*dzpos; // predicted dt
      costheta= (dtpos*dtpred+dzpos*dzpos)/sqrt((dtpos*dtpos +dzpos*dzpos)*(dtpred*dtpred+dzpos*dzpos));
  }
  else{  
    if(fDebug) cout<<" cmp-> laststriplink == NULL"<<endl;
    return pow(penalty,pndist);   
  }
    //
  // Return a 'probability' 
  //related to cos(theta) (penalty)^(skipped_planes)  (previous probabilities)
  if(fDebug) cout<<" End of Get Probability "<<endl;
  return (double)pow((costheta+1)/2,2)*pow(penalty, pndist)*lsum;
}


//-----------------------------------------------------------------------------
bool Anp::FillShortEvent::FillLink(StripLink& striplink, double ins,double ext, double ddx, StripIter strip, StripLink* link){
  (striplink).ins_value = ins;
  (striplink).ext_value = ext;
  (striplink).link_number = link->link_number+1;
  (striplink).ddx_value = ddx;
  (striplink).strip=strip;
  (striplink).laststriplink = link;

  return true;
}

                                                            
bool Anp::FillShortEvent::FillLink(StripLink& striplink, double ins,double ext, double ddx, StripIter strip){

  (striplink).ins_value = ins;
  (striplink).ext_value = ext;
  (striplink).link_number = 1;
  (striplink).ddx_value = ddx;
  (striplink).strip=strip;
  (striplink).laststriplink = NULL;

  return true;
}

---