loon: Anp::FillShortEvent Class Reference

00034   :fRecoHighE(3),
00035    fNEndPlanes(5),
00036    fMinADCHit(200),
00037    fNPlaneSkip(4),
00038    fTimeWindow((2.0*(18.83 + 0.1)*1.0e-9)),
00039    fStripWindow(4),
00040    fKeyBase(19000),
00041    fMinDLength(4),
00042    fCutOff(0.5),
00043    fKeyPass(-1),
00044    fKeyPass2(-1),
00045    fDirName("shortevent"),
00046    fDebug(false),
00047    fErase(false)
00048 {
00049 }

00341 {
00342   if(fDebug) cout<<" FillShortEvent:: BVD_BestLine---------------------------------------------"<<endl;
00343   //,,,,,,,,,,,,,,,,,,,,,,,                                                                                                                     
00344   // Initialize the vector                                                                                               
00345   //'''''''''''''''''''''''         
00346   if(map.size() ==0 )return -1;
00347   PlaneMapList map_list;
00348   for(PlaneMap::iterator mapit = map.begin(); mapit!=map.end();mapit++)
00349   {
00350     vector<Anp::FillShortEvent::StripLinkPtr> map_vec;
00351     for(vector<StripIter>::iterator it =(mapit->second).begin(); it!=(mapit->second).end();it++)
00352       {
00353         StripLink* str_link_it= new StripLink();
00354         if(!FillLink(*str_link_it, 0,0,0,(*it))){
00355           if(fDebug) cout<<" StripLink did not fill"<<endl;
00356           continue;
00357         }
00358         map_vec.push_back(str_link_it);
00359       }
00360     // For each plane get a vector of Strip 'Links'
00361     map_list[mapit->first] =map_vec;       
00362   }
00363 
00364   //......................
00365   // Sift through and link
00366   //''''''''''''''''''''''
00367 
00368   vector<Anp::FillShortEvent::StripLinkPtr> last_vec = map_list.begin()->second;
00369 
00370   // Loop over each plane
00371 
00372   for(PlaneMapList::reverse_iterator list_it = map_list.rbegin(); list_it!=map_list.rend(); list_it++){
00373     if((list_it == map_list.rbegin()))
00374     {
00375       last_vec = list_it->second;
00376       continue;
00377     }
00378 
00379     // loop over the vector of strips for this plane
00380     for(vector<StripLinkPtr>::iterator iter = (list_it->second).begin(); iter!=(list_it->second).end(); iter++)
00381       {
00382        
00383         StripLink*  striplink  = (*iter);
00384         striplink->ext_value=0.0;
00385         // loop over the vector of strips for the last plane 
00386         for(vector<StripLinkPtr>::iterator lp_it = last_vec.begin(); lp_it!=last_vec.end(); lp_it++)
00387           {
00388 
00389             // go through the linked list of all previous strips connected to this particlar strip
00390             // if any get a higher probability of being 'connected',despite the skipped plane this becomes the new strip to beat
00391             // This is now the prbability for "striplink" to be connected to "link"         
00392 
00393             StripLink*  link  = (*lp_it); // first strip in the chain
00394             double linkp = GetProbability(striplink, link);
00395 
00396             StripLink* linkn= link->laststriplink; // strip connected to this strip
00397             int skipper=0;
00398             while(linkn!=NULL)
00399             {                                                                                            
00400               if(skipper> 2) break;//fNPlaneSkip?
00401               double c_prob= GetProbability(striplink, linkn);
00402          
00403               if(c_prob > linkp)
00404               {
00405                 link = linkn;
00406                 linkp = c_prob;
00407               }
00408       
00409               linkn= linkn->laststriplink;            
00410               skipper++;
00411               if( c_prob< 0)
00412               {
00413                 if(fDebug) cout<<" probabiltiy is less than zero"<<endl;
00414                 break;
00415               }
00416             }
00417             if(fDebug) cout<<" FillShortEvent::BVD goodlength has  probability  : "<< linkp<<endl;
00418 
00419             double dtpos =((striplink)->strip)->TPos()-((link)->strip)->TPos();
00420             double dzpos= ((striplink)->strip)->ZPos()-((link)->strip)->ZPos();
00421             double ddx_value =0;
00422             if(link->link_number>2) ddx_value =((link->ddx_value)*(link->link_number-3) +((dtpos/dzpos- link->ins_value)/dzpos))/(link->link_number-2);
00423 
00424             if(dzpos==0)
00425             {
00426                 if(fDebug) cout<<" dzpos is zero"<<endl;
00427                 continue;
00428             }
00429             else if(linkp> striplink->ext_value)// If you have a larger probability
00430             {
00431               FillLink(*striplink, dtpos/dzpos, linkp, ddx_value, striplink->strip, link);           
00432             }
00433             else if(linkp== striplink->ext_value &&  (striplink->laststriplink==NULL  ||
00434                                                       sqrt(dtpos*dtpos +dzpos*dzpos) < sqrt(pow(striplink->strip->TPos() - striplink->laststriplink->strip->TPos(),2) + 
00435                                                                                             pow(striplink->strip->ZPos()- striplink->laststriplink->strip->ZPos(),2))))
00436             {                 
00437               FillLink(*striplink, dtpos/dzpos, linkp, ddx_value, striplink->strip, link);            
00438             }
00439           }     
00440         if( striplink->laststriplink!=NULL) ((striplink)->laststriplink)->end_link = false;     
00441       }
00442     last_vec = list_it->second;
00443   }
00444   
00445 
00446   int cnt=0;
00447   double sumpercnt =0.0;
00448   StripLinkPtr slink =NULL;
00449 
00450     for(PlaneMapList::reverse_iterator list_it = map_list.rbegin(); list_it!=map_list.rend(); list_it++){
00451       vector<Anp::FillShortEvent::StripLinkPtr> link_vec = list_it->second;
00452       cnt++;
00453       if(cnt<fMinDLength)       continue; 
00454 
00455       for(vector<StripLinkPtr>::iterator iter = link_vec.begin(); iter!=link_vec.end(); iter++)  {
00456         int link_n =(int) (*iter)->link_number;
00457         if(fDebug) cout<<" -- "<< link_n<<endl;
00458         if(link_n+fNPlaneSkip< map.size()) continue;
00459           if(pow((*iter)->ext_value,1.0/((double)link_n))*pow((double)link_n,2)  > sumpercnt &&link_n >=fMinDLength )  {         
00460           slink= (*iter);
00461           sumpercnt = pow((*iter)->ext_value,1.0/((double)link_n))*pow((double)link_n,2);           
00462         }
00463       }      
00464 
00465     }
00466 
00467     double size=0;
00468     if(slink!=NULL) size = slink->link_number;
00469   PlaneMap mapout;
00470   while(slink!=NULL){
00471     int plane = (slink->strip)->GetPlane();
00472     mapout[plane].push_back(slink->strip);    
00473     slink= slink->laststriplink;
00474   }
00475   if(!(size==mapout.size()) ) cout<<" FillShortEvent::BVD_bestline sizes don't match"<< mapout.size()<<" "<<size<<endl;
00476   if(fDebug) cout<<" FillShortEvent::BVD_bestline End -------   found line of length "<< mapout.size()<<" from "<<map.size() <<" planes ("<<size<<")"<< endl;
00477 
00478   map = mapout;
00479   if(mapout.size()==0) return 0;
00480   return  sumpercnt/pow((double)size,2.0);
00481   }

00099 {
00100   reg.Get("FillShortEventMaxRecoE",  fRecoHighE);
00101   reg.Get("FillShortEventMinNPlane",  fNEndPlanes);
00102   reg.Get("FillShortEventMinADCHit",  fMinADCHit);
00103   reg.Get("FillShortEventKeyBase",    fKeyBase);
00104   reg.Get("FillShortEventTimeWindow", fTimeWindow);
00105   reg.Get("FillShortEventStripWindow", fStripWindow);
00106   reg.Get("FillShortEventNPlaneSkip", fNPlaneSkip);
00107   reg.Get("FillShortEventMinBVDLength",  fMinDLength);
00108   reg.Get("FillShortEventCutOff", fCutOff);
00109   reg.Get("FillShortEventKeyPass", fKeyPass);
00110   reg.Get("FillShortEventKeyPass2", fKeyPass2);
00111 
00112   //
00113   // Read bool and string variables
00114   //
00115   Anp::Read(reg, "FillShortEventDebug",    fDebug);
00116   Anp::Read(reg, "FillShortEventErase",    fErase);
00117 
00118 
00119   if(reg.KeyExists("PrintConfig"))
00120    {
00121      cout << "FillShortEvent::Config" << endl;
00122      cout << "KeyBase = " << fKeyBase<< endl;
00123      cout << "BDV Length = " << fMinDLength<< endl;
00124      cout << "TimeWindow = "<< fTimeWindow<<endl;
00125      cout << "KeyPass = "<<fKeyPass<<endl;
00126      cout << "KeyPass2 = "<<fKeyPass2<<endl;
00127      cout<<  "Erase ="<< fErase<<endl;
00128    }
00129 }

00139 {
00140   cout<<" Nmissed (too short)      "<< fTooShort<<endl;
00141   cout<<" Nmissed (bvd too short)  "<< fBVDmiss<<endl;
00142 }

00144                                                                   {
00145 
00146   PlaneMap mapu;
00147   PlaneMap mapv;
00148   const Basic eventb= event.GetBasic();
00149 
00150   if(eventb.NUPlane()<fNEndPlanes &&eventb.NVPlane()<fNEndPlanes){
00151     fTooShort++;
00152     return false;
00153   }
00154   if(!FillShortEvent::GetPlaneStrips(mapu, mapv, event, record))   return false;
00155  
00156   std::map<int, double> ansv;  std::map<int, double> ansu;
00157 
00158   ansu[0]=(double)mapu.size();
00159   ansv[0]=(double)mapv.size();
00160   
00161   bool ugood = false; 
00162   bool vgood = false;
00163   if(fDebug) cout<<" FillShortEvent::Explore - find event quantities "<<endl;
00164   if(mapu.size()>=fNEndPlanes  &&FillShortEvent::GetEventQuantities(mapu,&ansu))    ugood=true;
00165   if(mapv.size()>=fNEndPlanes && FillShortEvent::GetEventQuantities(mapv,&ansv))    vgood=true;
00166   
00167   std::map< int, double> ans;
00168 
00169   if(!(ugood) && !(vgood) )
00170   {
00171     fBVDmiss ++;
00172     return false;
00173   }
00174   else if(!ugood)    ans = ansv;
00175   else if(!vgood)    ans = ansu;
00176   else if(ansu[3] >= ansv[3]) ans = ansu;
00177   else ans = ansv;
00178   if(fDebug) cout<<" choose line with size "<<ans[3]<<endl;
00179 
00180   ans[12]= event.Gev();
00181 
00182   for(unsigned int ii=0; ii<ans.size(); ii++)    event.Add(fKeyBase+ii+1, ans[ii]);
00183 
00184   return true;
00185 }

00556                                                                                                         {
00557 
00558   (striplink).ins_value = ins;
00559   (striplink).ext_value = ext;
00560   (striplink).link_number = 1;
00561   (striplink).ddx_value = ddx;
00562   (striplink).strip=strip;
00563   (striplink).laststriplink = NULL;
00564 
00565   return true;
00566 }

00544                                                                                                                          {
00545   (striplink).ins_value = ins;
00546   (striplink).ext_value = ext;
00547   (striplink).link_number = link->link_number+1;
00548   (striplink).ddx_value = ddx;
00549   (striplink).strip=strip;
00550   (striplink).laststriplink = link;
00551 
00552   return true;
00553 }

00220 {
00221   //
00222   // Get Some basic information (# hits, # planes)
00223   //
00224   double NHits=0.0; //(number of hits)                                
00225   int count=0;      // count of planes                                                                
00226   double maxT=-10000;
00227   double minT= 10000;
00228 
00229   PlaneMap::iterator last_pit = plane_m.begin();
00230 
00231   for (PlaneMap::iterator pit = plane_m.begin(); pit!=plane_m.end();++pit){
00232      vector<StripIter> &svec  = pit -> second;
00233     if(count>2)
00234       {
00235       NHits+=svec.size();
00236   
00237       // count hits for all but first two planes (can be noisy)                   
00238       for(vector<StripIter>::iterator it = svec.begin(); it!=svec.end(); it++)
00239       {    
00240         if( (*it)->TPos()> maxT ) maxT = (*it)->TPos();
00241         else if((*it)->TPos()<minT) minT= (*it)->TPos();
00242       }
00243     }
00244     count++;
00245   }
00246 
00247   if(count<2 || NHits<=0)
00248   {
00249     if(fDebug) cout<<" FillShortEvent:: GetEventQuantities - return false for too small vector"<<endl;
00250     return false; 
00251   }
00252 
00253  //
00254  // Get the BVD best line 
00255  //
00256 
00257  PlaneMap dv_map =plane_m;
00258  double dijkstra = BVD_BestLine(dv_map);
00259  double size = dv_map.size();
00260  if(dijkstra<fCutOff ||size<0)
00261  {
00262    if(fDebug) cout<<" dijkstra best line not found. Size = "<<size<<endl;
00263    return false;
00264  }
00265 
00266  if(count>2*size) cout<<" Number of planes is much greater than found path "<< count<<" "<<size<<endl;
00267 
00268  if(size>=fMinDLength){
00269 
00270    double sumph =-1; double sum_near_ph=-1;
00271    double sigph =-1;
00272    double lastph=-1;
00273    vector<double> time;
00274    
00275    int c=0;
00276    //
00277    // first loop
00278    //
00279    for(PlaneMap::reverse_iterator si = dv_map.rbegin(); si!=dv_map.rend(); si++){
00280      if(sumph==-1){
00281        sumph=0;
00282        lastph=0;
00283      }
00284      const vector<StripIter> &svec  = si -> second;
00285      sumph += (svec[0])->SigCor();
00286      time.push_back(svec[0]->Time());
00287      if(c<3)lastph+=(svec[0])->SigCor();
00288      c++;
00289    }
00290    
00291    std::sort(time.begin(), time.end());
00292    //
00293    // second loop
00294    //
00295    int nearstrips=0;
00296    for(PlaneMap::iterator si = dv_map.begin(); si!=dv_map.end(); si++){
00297      if(sum_near_ph<0) sum_near_ph=0;
00298      vector<StripIter> &svec  = si -> second;
00299      int minstrip = svec[0]->GetStrip();
00300      int maxstrip = svec[0]->GetStrip();
00301      for( vector<StripIter>::iterator its = svec.begin(); its!=svec.end();its++){
00302        if((*its)->GetStrip()< minstrip) minstrip = (*its)->GetStrip();
00303        if((*its)->GetStrip()> maxstrip) maxstrip = (*its)->GetStrip();
00304      }
00305      vector<StripIter> pmit_vec = plane_m[si->first];
00306      for(vector<StripIter>::iterator it = pmit_vec.begin(); it!=pmit_vec.end(); it++){
00307        int its=(*it)->GetStrip();
00308        
00309        if( ( its> minstrip-fStripWindow && its<maxstrip+fStripWindow)
00310            && (*it)->Time()< time.back()+fTimeWindow && (*it)->Time()>time.front()-fTimeWindow){
00311          sum_near_ph+= (*it)->SigCor(); 
00312          nearstrips++;
00313        }
00314      }
00315      
00316      if(sigph==-1) sigph =0;
00317      sigph +=  ((svec[0])->SigCor() -(sumph/size))*((svec[0])->SigCor() -(sumph/size))/size;
00318    }
00319    vect->insert( pair<int,double> (1,(double)NHits/(count-2)));  // 1
00320    vect->insert( pair<int,double> (2,(maxT-minT)/(count-2)));     // 2
00321    vect->insert( pair<int,double> (3,size));                     // 3  
00322    vect->insert( pair<int,double> (4,dijkstra));                 // 4    
00323    vect->insert( pair<int,double> (5,sumph/size));               // 5    
00324    vect->insert( pair<int,double> (6,sqrt(sigph)));              // 6
00325    vect->insert( pair<int,double> (7,lastph));                   // 7
00326    vect->insert( pair<int,double> (8,(sumph+sum_near_ph)/size)); // 8
00327    vect->insert( pair<int,double> (9,size/count));               // 9
00328    vect->insert( pair<int,double> (10,(time.back()-time.front())
00329                                    /(1.0e-9*time.size())));      // 10
00330    vect->insert( pair<int,double> (11,nearstrips/size));         // 11
00331   
00332    return true;
00333  }
00334 
00335  if(fDebug) cout<<"FillShortEvent::GetEventQuantities  Return false because size is less than minimum length "<<size<<endl;
00336  return false;
00337 
00338 }

00189 {
00190   for(StripIter istrip = record.StripBeg(); istrip != record.StripEnd(); ++istrip)
00191     {
00192       if(!istrip -> MatchEvt(ievent.EventIndex()))
00193         {
00194           continue;
00195         }
00196 
00197       if(record.GetHeader().IsNear() && istrip -> GetPlane() > 120)
00198         {
00199           return false;
00200         }
00201 
00202       if(istrip -> SigCor() < fMinADCHit)
00203         {
00204           continue;
00205         }
00206 
00207       if(istrip->IsUview()){
00208         vector<StripIter> &svec = plane_mapu[istrip -> GetPlane()];
00209         svec.push_back(istrip);
00210       }
00211       else if(istrip->IsVview()){
00212         vector<StripIter> &svec = plane_mapv[istrip -> GetPlane()];
00213         svec.push_back(istrip);
00214       }
00215     }
00216   return true;
00217 }

00485                                                                       {
00486   
00487   //
00488   // Penalty for missing a plane
00489   //
00490   double penalty = 0.728;
00491   // 
00492   // Get the derivative - change in position
00493   //
00494   double dtpos =((cur)->strip)->TPos()-((cmp)->strip)->TPos();
00495   double dzpos= ((cur)->strip)->ZPos()-((cmp)->strip)->ZPos();
00496 
00497   //
00498   // Get the last extremum value, instantaneous (derivative) value
00499   //
00500   double lsum =(cmp)->ext_value;
00501   double lins =(cmp)->ins_value;
00502   //
00503   // Get the number of skipped planes
00504   //
00505   double pndist = (cmp->strip->GetPlane() - cur->strip->GetPlane())/2 -1;
00506   if(pndist<0 || dzpos>=0){
00507     if(fDebug) cout<<" Number of skipped planes is negative! "<<pndist<<" "<<dzpos<<endl;
00508     return -1;
00509   }
00510 
00511   double pred_dtdz =0; //predicted change in step       
00512   double dtpred =0;
00513   double costheta=0;
00514   //
00515   // Get cosine (delta theta)  
00516   //where delta theta is the angle between pervious derivative and current
00517   if((cmp)->laststriplink!=NULL)
00518   {
00519     if(fDebug) cout<<" cmp->laststriplink ! = NULL "<<endl;
00520     // not comparing to first link
00521     if(cmp->laststriplink->laststriplink!=NULL){
00522       pred_dtdz=(cmp)->ddx_value;
00523       if(pred_dtdz*dzpos > 10){
00524         //      cout<<" Very large change in pred_dtdz"<< pred_dtdz<<endl;
00525         pred_dtdz =0.0;
00526       }
00527     }
00528       dtpred = (lins + pred_dtdz*dzpos)*dzpos; // predicted dt
00529       costheta= (dtpos*dtpred+dzpos*dzpos)/sqrt((dtpos*dtpos +dzpos*dzpos)*(dtpred*dtpred+dzpos*dzpos));
00530   }
00531   else{  
00532     if(fDebug) cout<<" cmp-> laststriplink == NULL"<<endl;
00533     return pow(penalty,pndist);   
00534   }
00535     //
00536   // Return a 'probability' 
00537   //related to cos(theta) (penalty)^(skipped_planes)  (previous probabilities)
00538   if(fDebug) cout<<" End of Get Probability "<<endl;
00539   return (double)pow((costheta+1)/2,2)*pow(penalty, pndist)*lsum;
00540 }

00132 {
00133   fTooShort=0;
00134   fBVDmiss =0;
00135   return true;
00136 }

00057 { 
00058 
00059 
00060   EventIterator ievent = record.EventBegIterator();
00061   while(ievent != record.EventEndIterator())
00062   { 
00063     //
00064     // Fill truth variable used by training algorithm(s)
00065     //
00066 
00067     if(!record.GetHeader().IsData() && record.TruthBeg() != record.TruthEnd())
00068       {
00069         const TruthIter truth = record.FindTruth(*ievent);
00070         if(truth == record.TruthEnd()) return true;
00071           
00072         double iscc= 0;
00073         if(truth-> IsCC()) iscc = 1;
00074             
00075         if(!(ievent->Add(fKeyBase, iscc)))
00076           {
00077             cerr << "FillShortEvent::Run - failed to insert value at " << fKeyBase+300 << " key" << endl;
00078           }
00079       }
00080 
00081     TrackIter track_it = LongestTrack(*ievent,record);
00082     if(track_it==record.TrackEnd()||
00083        (fKeyPass>0 && fKeyPass2>0 &&!track_it->KeyExists(fKeyPass) &&   !track_it->KeyExists(fKeyPass2)))
00084       {
00085       if(!FillShortEvent::Explore(*ievent, record) && fErase)
00086        {
00087         ievent = record.Erase(ievent);
00088         continue;
00089       }
00090     }
00091 
00092     ievent++;
00093   }
00094   if(fErase)  Anp::CleanRecord(record.EventBeg(), record.EventEnd(), record);
00095   return true;
00096 }

Anp::FillShortEvent Class Reference

Public Member Functions

Private Types

Private Member Functions

Private Attributes

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation


Public Member Functions
	FillShortEvent ()
virtual	~FillShortEvent ()
bool	Run (Record &record)
void	Config (const Registry &reg)
bool	Init (const Header &)
void	End (const DataBlock &)
bool	Explore (Event &ievent, const Record &record)
Private Types
typedef StripLink *	StripLinkPtr
typedef std::map< unsigned short, std::vector< Anp::StripIter > >	PlaneMap
typedef std::map< unsigned short, std::vector< StripLinkPtr > >	PlaneMapList
Private Member Functions
bool	GetPlaneStrips (PlaneMap &plane_mapu, PlaneMap &plane_mapv, const Event ievent, const Record &record)
bool	GetEventQuantities (PlaneMap plane_m, std::map< int, double > *vect)
double	BVD_BestLine (PlaneMap &map)
double	GetProbability (StripLink cur, StripLink cmp)
bool	FillLink (StripLink &striplink, double ins, double ext, double ddx, StripIter strip, StripLink *link)
bool	FillLink (StripLink &striplink, double ins, double ext, double ddx, StripIter strip)
Private Attributes
double	fRecoHighE
double	fNEndPlanes
double	fMinADCHit
double	fNPlaneSkip
double	fTimeWindow
int	fStripWindow
int	fKeyBase
int	fMinDLength
double	fCutOff
int	fKeyPass
int	fKeyPass2
std::string	fDirName
TDirectory *	fDir
bool	fDebug
bool	fErase
bool	fPlot
int	fNMiss
int	fTooShort
int	fBVDmiss