PrimaryShowerFinder Class Reference

#include <PrimaryShowerFinder.h>

List of all members.

Public Member Functions

 PrimaryShowerFinder ()
 ~PrimaryShowerFinder ()
int FindPrimaryShower (Managed::ClusterManager *cm)
Particle3DGetParticle3D ()
ChainGetChain (int view)
int FoundSingleViewPrimaryShower ()
TH2D * GetHoughMap (int view)
std::vector< HoughLine > * GetHoughLines (int view)
int FindFirstIntersection (int view, double &t, double &z)
void Bundle (int view)
void SetVertex (double u, double v, double z)

Static Public Member Functions

static void LoadCloseHits (HoughLine *hl, Managed::ClusterManager *cm, int view, double max_tdist, int limit_to_current_size=0)

Public Attributes

int ran

Static Public Attributes

static TH2D * intU = 0
static TH2D * intV = 0
static ChainHelperchu = 0
static ChainHelperchv = 0

Private Member Functions

int CheckChainQuality (Chain *c)
ChainFindShowerChain (Managed::ClusterManager *cl, int view)
void Reset (int reset_vertex=1)
void MakeParticle3D ()
void MakeHoughMap (int view)
void SaveHitGroup (TH2D *his, TH2D *saver, double save_val, double with_val, int curx, int cury)
void GetPeakAreaAverage (double &x, double &y, double &val, int &cnt, int curx, int cury, TH2D *hist, std::vector< double > &peakvalue, std::vector< int > &peakstillgood, std::vector< int > &peakbinx, std::vector< int > &peakbiny)
void ExpandShowerChain (Chain *chain, int view, double start_z, double end_z)
ChainExtractMuon (Chain *c)
void DumpParticle3D ()
void DumpHoughLines (int view)
void MakeChains (int view)
ChainMakeShowerChain (int view)
void MergeChainClusters (Chain *ch, std::vector< int > *clusters)
int CheckChainOverlap (Chain *chain_u, Chain *chain_v)
void ClearFrontVertex (Chain *chain_u, Chain *chain_v)
HoughLineSplitHoughLine (Chain *c, HoughLine *hl, Managed::ClusterManager *mycm=0)
void FindHoughLineMatches ()
void CleanHoughLines (int view, std::vector< HoughLine > *hhh)
 ClassDef (PrimaryShowerFinder, 1)

Private Attributes

int foundparticle
Managed::ClusterManagercluster_manager
int single_view_long_shower
std::vector< HoughLinehoughlinesU
std::vector< HoughLinehoughlinesV
std::vector< std::pair< int,
int > > 
houghlineMatch
double vtx_u
double vtx_v
double vtx_z
int foundvertex
int foundvertexU
int foundvertexV

Static Private Attributes

static Particle3Dfoundparticle3d = 0
static Chainchain_u = 0
static Chainchain_v = 0
static TH2D * houghmapU = 0
static TH2D * houghmapV = 0

Detailed Description

Definition at line 17 of file PrimaryShowerFinder.h.


Constructor & Destructor Documentation

PrimaryShowerFinder::PrimaryShowerFinder (  ) 

Definition at line 45 of file PrimaryShowerFinder.cxx.

References Reset().

00046 {
00047 
00048         Reset();
00049 }

PrimaryShowerFinder::~PrimaryShowerFinder (  ) 

Definition at line 51 of file PrimaryShowerFinder.cxx.

References Reset().

00052 {
00053         Reset();
00054 }


Member Function Documentation

void PrimaryShowerFinder::Bundle ( int  view  ) 

Definition at line 1052 of file PrimaryShowerFinder.cxx.

References cluster_manager, done(), houghlinesU, houghlinesV, LoadCloseHits(), HoughLine::offset_t, HoughLine::offset_z, HoughLine::r, Munits::second, HoughLine::sum_e, and HoughLine::theta.

Referenced by FindPrimaryShower().

01053 {
01054         std::vector<HoughLine> * hl;
01055         if(view==2)hl=&houghlinesU;
01056         else if(view==3)hl=&houghlinesV;
01057         else return;
01058 
01059 
01060         std::vector<std::pair<int,int> >to_merge;
01061         
01062         //iterate over each line
01063         for(unsigned int i=0;i<(*hl).size();i++)
01064         {
01065         
01066 //                      printf("houghline with e %f clusters %d (t,z) (%f,%f) (%f,%f) chi2/ndf %f  offz %f expt %f phi %f\n", (*hl)[i].sum_e, (*hl)[i].ncluster, (*hl)[i].start_t, (*hl)[i].start_z, (*hl)[i].end_t, (*hl)[i].end_z, (*hl)[i].chi2/(*hl)[i].ncluster,(*hl)[i].offset_z,(*hl)[i].GetExpectedT((*hl)[i].offset_z),(*hl)[i].phi);
01067                         
01068                 //nested loop
01069                 for(int j=i+1;j<(int)(*hl).size();j++)
01070                 {
01071 
01072                         //if(i==j)continue;
01073                         
01074                         //did we already match this chain?
01075                         int done=0;
01076                         for(int k=0;k<(int)to_merge.size();k++)
01077                         {
01078                                 if(to_merge[k].first==j || to_merge[k].second==j)
01079                                 {
01080                                         done=1;
01081                                         break;
01082                                 }
01083                         }
01084                         if(done)continue;
01085                         
01086                 //
01087 //                      printf("\t %d ccc e %f clusters %d (t,z) (%f,%f) (%f,%f) chi2/ndf %f  offz %f expt %f phi %f\n", j,(*hl)[j].sum_e, (*hl)[j].ncluster, (*hl)[j].start_t, (*hl)[j].start_z, (*hl)[j].end_t, (*hl)[j].end_z, (*hl)[j].chi2/(*hl)[j].ncluster,(*hl)[j].offset_z,(*hl)[j].GetExpectedT((*hl)[j].offset_z),(*hl)[j].phi);
01088                 
01089                         if(fabs((*hl)[i].phi - (*hl)[j].phi)<0.05)
01090                         {
01091                                 if(fabs( (*hl)[i].GetExpectedT((*hl)[i].offset_z) - (*hl)[j].GetExpectedT((*hl)[i].offset_z)) < 0.05)
01092                                 {
01093                                         to_merge.push_back(make_pair(i,j));
01094 //                                      printf("tm %d %d %f %f %f\n",i,j,(*hl)[i].offset_z,(*hl)[i].GetExpectedT((*hl)[i].offset_z),(*hl)[j].GetExpectedT((*hl)[i].offset_z));
01095                                 
01096                                 }
01097                         
01098                         
01099                         }
01100         
01101                 }
01102         }
01103         
01104         if(to_merge.size()<1)return;
01105         
01106         std::vector<HoughLine> toKeep;
01107 
01108         //keep non merged lines
01109         for(unsigned int i=0;i<hl->size();i++)
01110         {
01111                 int keep=1;
01112                 for(unsigned int j=0;j<to_merge.size();j++)
01113                 {
01114                         if(to_merge[j].first ==(int)i || to_merge[j].second ==(int)i)
01115                         {
01116                                 keep=0;
01117                                 break;
01118                         }
01119                 }
01120                 if(keep){
01121                         toKeep.push_back((*hl)[i]);
01122 //                      printf("not merging %d\n",i);
01123                 }
01124         
01125         }
01126 
01127 
01128         for(unsigned int i=0;i<to_merge.size();i++)
01129         {
01130 /*
01131                 //simple for now... just pick one out of the group
01132                 if(i==0 || (i>0 && (to_merge[i].first != to_merge[i-1].first) ) )
01133                 {
01134 
01135                         toKeep.push_back((*hl)[to_merge[i].first]);
01136                         
01137                         printf("%d \n",to_merge[i].first);
01138                 }
01139 */
01140                 std::vector<int> mlist;
01141                 
01142                 mlist.push_back(to_merge[i].first);
01143                 for(unsigned int j=i+1;j<to_merge.size();j++)
01144                 {
01145                         if(to_merge[j].first !=to_merge[i].first)break;
01146                         mlist.push_back(to_merge[j].second);
01147                 }
01148 
01149 
01150                 //now we have a list of houghlines to merge...
01151                 //do a weighted average by energy
01152                 
01153                 double theta=0;
01154                 double r=0;
01155                 double offset_t=0;
01156                 double offset_z=0;
01157                 double sum_e=0;
01158                 
01159                 for(unsigned int j=0;j<mlist.size();j++)
01160                 {
01161                         HoughLine *h = &(*hl)[mlist[j]];
01162                         theta += h->theta*h->sum_e;
01163                         r += h->r*h->sum_e;
01164                         offset_t += h->offset_t*h->sum_e;
01165                         offset_z += h->offset_z*h->sum_e;
01166                         sum_e+=h->sum_e;
01167                 }
01168                 
01169                 if(sum_e<0.1)continue;
01170                 theta/=sum_e;
01171                 r/=sum_e;
01172                 offset_t/=sum_e;
01173                 offset_z/=sum_e;
01174                 
01175                 HoughLine l( theta,  r,  offset_t,  offset_z);
01176                 LoadCloseHits(&l,cluster_manager,view,0.0412);
01177                 toKeep.push_back(l);
01178 
01179         }
01180 
01181 
01182 
01183         
01184         (*hl)=toKeep;
01185         
01186 //      printf("# of hough lines %d\n",hl->size());
01187 
01188 }

int PrimaryShowerFinder::CheckChainOverlap ( Chain chain_u,
Chain chain_v 
) [private]

Definition at line 2164 of file PrimaryShowerFinder.cxx.

References Chain::end_z, Msg::kDebug, MSG, and Chain::start_z.

02165 {
02166 
02167         
02168 
02169         
02170         //now make sure that there is sufficient overlap
02171         double require_overlap_distance=1.0;
02172         
02173         double start_z = chain_u->start_z > chain_v->start_z ? chain_u->start_z : chain_v->start_z;
02174         double end_z = chain_u->end_z < chain_v->end_z ? chain_u->end_z : chain_v->end_z;
02175 
02176         MSG("PrimaryShowerFinder",Msg::kDebug)<<"overlap "<<end_z-start_z<<"\n";
02177         
02178         if(end_z-start_z < require_overlap_distance)return 0;
02179         return 1;
02180         
02181         
02182 
02183 }

int PrimaryShowerFinder::CheckChainQuality ( Chain c  )  [private]

Definition at line 2187 of file PrimaryShowerFinder.cxx.

References Chain::e, Chain::entries, Msg::kDebug, max, min, MSG, Chain::muonfrac, and Chain::z.

Referenced by FindPrimaryShower().

02188 {
02189         if(!c)return 0;
02190         if(c->entries<1)return 0;       
02191         
02192         int qual=1;
02193         
02194         //look at Shower-like strip fraction
02195         //and sparsity (# planes with hits/#planes from front to back of Shower)
02196         
02197         int Showerlike=0;
02198         int hitplane=0;
02199         
02200         double last_hitplane=-100;
02201         
02202         double min = 1000000;
02203         double max = 0;
02204         for(int i=0;i<c->entries;i++)
02205         {
02206                 if(c->e[i]>0.5 && c->e[i]<2.5)Showerlike++;
02207                 if(fabs(last_hitplane - c->z[i])>0.03)
02208                 {
02209                         hitplane++;
02210                         last_hitplane=c->z[i];
02211                 }
02212                 
02213                 min = c->z[i]< min? c->z[i]:min;
02214                 max = c->z[i]> max? c->z[i]:max;
02215         }
02216         
02217         double Showerfrac = (double)Showerlike / (double)c->entries;
02218         double sparsity = max-min?(double)hitplane / (max-min) / 7.08 :1; //length/size of u+v plane
02219         
02220 
02221 //definition of shwerfrac is wrong!     
02222 //      if (Showerfrac < 0.5) qual=0;
02223         if (sparsity < 0.8) qual=0;
02224         
02225         if(c->muonfrac>0.5)qual=0;
02226         
02227         
02228         MSG("PrimaryShowerFinder",Msg::kDebug)<<"MUON FRAC "<<c->muonfrac <<"Shower chain check  Showerfrac "<< Showerfrac<<" sparsity "<< sparsity<<"  qual "<< qual<<"\n";
02229 
02230         return qual;
02231 }

PrimaryShowerFinder::ClassDef ( PrimaryShowerFinder  ,
 
) [private]
void PrimaryShowerFinder::CleanHoughLines ( int  view,
std::vector< HoughLine > *  hhh = 0 
) [private]

Definition at line 1191 of file PrimaryShowerFinder.cxx.

References cluster_manager, done(), Managed::ClusterManager::GetCluster(), houghlinesU, houghlinesV, len, and Managed::ManagedCluster::z.

Referenced by FindPrimaryShower(), and MakeChains().

01192 {
01193 
01194         
01195         std::vector<HoughLine> * hl;
01196         if(!hhh)
01197         {
01198                 if(view==2)hl=&houghlinesU;
01199                 else if(view==3)hl=&houghlinesV;
01200                 else return;
01201         }else{
01202                 hl=hhh;
01203         }
01204         
01205         
01206         
01207         std::vector<int>todel;
01208         //remove all but first vertical hough lines in this view
01209         
01210         //first hit z is:
01211         double first_z = 1000;
01212         int keep_v=-1;
01213         for(unsigned int i=0;i<hl->size();i++)
01214         {
01215                 if(fabs((*hl)[i].phi-3.141592/2 )<0.1){
01216                         if((*hl)[i].start_z < first_z)
01217                         {
01218                                 first_z = (*hl)[i].start_z;
01219                                 keep_v=i;
01220                         }
01221                 }
01222         }                       
01223                 
01224         for(int i=0;i<(int)hl->size();i++)
01225         {
01226                 if(fabs((*hl)[i].phi-3.141592/2 )<0.1){
01227                         if(keep_v!=i)todel.push_back(i);
01228                 }
01229         }       
01230         
01231                 
01232                         
01233         //remove hough lines that have the closest distance between a pair of adjacent hits larger than some value
01234         double min_dist = 0.2;
01235         for(unsigned int i=0;i<hl->size();i++)
01236         {
01237         
01238 
01239                 double small_dist=1000;
01240                 int done=0;
01241                 for(unsigned int j=0;j<(*hl)[i].cluster_id.size();j++)
01242                 {
01243                         
01244                         for(unsigned int k=j+1;k<(*hl)[i].cluster_id.size();k++)
01245                         {
01246                                 int id1 = (*hl)[i].cluster_id[j];
01247                                 int id2 = (*hl)[i].cluster_id[k];
01248                                 
01249                                 Managed::ManagedCluster * c1 = cluster_manager->GetCluster(id1);
01250                                 Managed::ManagedCluster * c2 = cluster_manager->GetCluster(id2);
01251                                 if(!c1 || !c2)continue;
01252                                 
01253                         //      double dist = sqrt((c1->t-c2->t)*(c1->t-c2->t) + (c1->z-c2->z)*(c1->z-c2->z));
01254                                 
01255                         //      //adjust the distance based on the angle of the houghline
01256                         //      dist *=cos((*hl)[i].phi);
01257                         
01258                                 double dist = fabs(c1->z-c2->z);
01259                                 
01260                                 if(dist < small_dist)small_dist=dist;           
01261                                 
01262                 //              printf("dist %f\n",dist);               
01263                                 
01264                                 if(small_dist < min_dist)
01265                                 {
01266                                         done=1;
01267                                         break;
01268                                 }
01269                         }
01270                         if(done)break;  
01271                 }
01272                 if(done)continue;
01273                 
01274                 //if we are here.. we want to delete it
01275                 todel.push_back(i);
01276                 
01277         }
01278 
01279 
01280         //hit density
01281         for(unsigned int i=0;i<hl->size();i++)
01282         {
01283                 //double len = sqrt(((*hl)[i].end_z-(*hl)[i].start_z)*((*hl)[i].end_z-(*hl)[i].start_z) + ((*hl)[i].end_t-(*hl)[i].start_t)*((*hl)[i].end_t-(*hl)[i].start_t) );
01284                 
01285                 double len=fabs((*hl)[i].end_z-(*hl)[i].start_z);
01286                 double frac = len>0?(*hl)[i].ncluster/(len/0.07):0;
01287                 
01288                 if(frac <0.5)todel.push_back(i);
01289         }       
01290 
01291 
01292 
01293         //need to save?
01294         if(todel.size()<1)return;
01295         
01296         std::vector<HoughLine> temp;
01297         for(int i=0;i<(int)hl->size();i++)
01298         {
01299                 int save=1;
01300                 for(int j=0;j<(int)todel.size();j++)
01301                 {
01302                         if(todel[j]==i)
01303                         {
01304                                 save=0;
01305                                 break;
01306                         }
01307                 }       
01308                 if(!save)continue;
01309                 temp.push_back((*hl)[i]);
01310         }
01311 
01312         *hl=temp;
01313 
01314 
01315 }

void PrimaryShowerFinder::ClearFrontVertex ( Chain chain_u,
Chain chain_v 
) [private]

Definition at line 2138 of file PrimaryShowerFinder.cxx.

References Chain::ClearHits(), Chain::cluster_id, Chain::e, Chain::entries, Chain::insert(), Chain::start_z, Chain::t, and Chain::z.

02139 {
02140         double start_z_u = chain_u->start_z;
02141         double start_z_v = chain_v->start_z;
02142 
02143         if(fabs(start_z_u-start_z_v)<0.04)return; //close enough
02144         
02145         double start_z = start_z_u < start_z_v ? start_z_v : start_z_u - 0.05;
02146         
02147         Chain * toclear =  start_z_u < start_z_v ? chain_u : chain_v;
02148         //now delete all hits up to start_z
02149         
02150         Chain tmp = *toclear;
02151         
02152         toclear->ClearHits();
02153         
02154         for(int i=0;i<tmp.entries;i++)
02155         {
02156                 if(tmp.z[i]>start_z)
02157                         toclear->insert(tmp.t[i], tmp.z[i], tmp.e[i], tmp.cluster_id[i]);
02158         }
02159         
02160 }

void PrimaryShowerFinder::DumpHoughLines ( int  view  )  [private]

Definition at line 3419 of file PrimaryShowerFinder.cxx.

References houghlinesU, houghlinesV, MsgService::Instance(), and Msg::kDebug.

03420 {
03421         
03422         if(!MsgService::Instance()->IsActive("PrimaryShowerFinder",Msg::kDebug))return;
03423         
03424         printf("\nView %d\n",view);
03425         std::vector<HoughLine>* houghlines=0;
03426         if(view==2)houghlines=&houghlinesU;
03427         if(view==3)houghlines=&houghlinesV;
03428         if(!houghlines)return;
03429 
03430         for(int i=(*houghlines).size()-1;i>-1;i--)
03431         {
03432                 if((*houghlines)[i].ncluster<1)continue;
03433                 printf("houghline with e %f clusters %d (t,z) (%f,%f) (%f,%f) chi2/ndf %f  offz %f expt %f phi %f\n", (*houghlines)[i].sum_e, (*houghlines)[i].ncluster, (*houghlines)[i].start_t, (*houghlines)[i].start_z, (*houghlines)[i].end_t, (*houghlines)[i].end_z, (*houghlines)[i].chi2/(*houghlines)[i].ncluster,(*houghlines)[i].offset_z,(*houghlines)[i].GetExpectedT((*houghlines)[i].offset_z),(*houghlines)[i].phi);
03434         }
03435 
03436         printf("\n");
03437 }

void PrimaryShowerFinder::DumpParticle3D (  )  [private]

Definition at line 2293 of file PrimaryShowerFinder.cxx.

References Particle3D::avg_rms_t, Particle3D::calibrated_energy, Particle3D::chain, Particle3D::chainhit, Particle3D::e, ParticleType::em, Particle3D::emfit_a, Particle3D::emfit_b, Particle3D::emfit_chisq, Particle3D::emfit_e0, Particle3D::emfit_ndf, Particle3D::emfit_prob, ParticleType::emshort, Particle3D::end_u, Particle3D::end_v, Particle3D::end_z, Particle3D::entries, foundparticle3d, Msg::kDebug, MSG, ParticleType::muon, Particle3D::muonfrac, Particle3D::particletype, ParticleType::pi0, ParticleType::prot, Particle3D::rms_t, Particle3D::shared, Particle3D::start_u, Particle3D::start_v, Particle3D::start_z, Particle3D::types, Particle3D::u, ParticleType::uniquemuon, Particle3D::v, Particle3D::view, and Particle3D::z.

Referenced by FindPrimaryShower().

02294 {
02295         ostringstream s;
02296         
02297         s << "Long Shower Particle3D found "<<endl;
02298 
02299                 Particle3D * p = foundparticle3d;
02300                 if (p==0)return;
02301         
02302         
02303 
02304         
02305                 s << "\n---Particle3d "  << "---\nstart, stop (u,v,z) (" << p->start_u << ", "<< p->start_v << ", " << p->start_z <<") (" <<p->end_u<<", "<< p->end_v << ", " << p->end_z <<")"<<endl;
02306                 s << "entries "<< p->entries << "  muonlike " << p->muonfrac <<endl;
02307                 
02308 
02309                 
02310                 s<<"types: ";
02311                 for(unsigned int j=0;j<p->types.size();j++)
02312                 {
02313                         switch( p->types[j].type)
02314                         {
02315                                 case    ParticleType::em:
02316                                         s<<"em ";
02317                                         break;
02318                                 case ParticleType::emshort:
02319                                         s<<"emshort  ";
02320                                         break;                          
02321                                 case ParticleType::muon:
02322                                         s<<"muon ";
02323                                         break;
02324                                 case ParticleType::prot:
02325                                         s<<"prot ";
02326                                         break;          
02327                                 case ParticleType::pi0:
02328                                         s<<"pi0 ";
02329                                         break;
02330                                 case ParticleType::uniquemuon:
02331                                         s<<"uniqueShower ";
02332                                         break;  
02333                         
02334                         }
02335                                 
02336                 }
02337                 s<<endl;
02338                 
02339                 s<<"particletype : "<<p->particletype<<endl;
02340                 
02341                 
02342                 s<<"par a " << p->emfit_a << " par b "<< p->emfit_b << " par e0 "<< p->emfit_e0 << " cale "<<p->calibrated_energy<<" chisq " << p->emfit_chisq<<" ndf " << p->emfit_ndf<<endl;
02343                 s<<"emprob " << p->emfit_prob <<" avg rms_t "<<p->avg_rms_t<<endl;
02344                                 
02345                 s << "spoints (u,v,z,e - chain, chainhit, chainview - shared - rms_t - view) : ";
02346                 for(int j=0;j<p->entries;j++)
02347                         s << "(" << p->u[j]<<", "<<p->v[j]<<", "<<p->z[j]<<", "<<p->e[j]<<" - " << p->chain[j] <<", "<<p->chainhit[j]<<", "<<p->view[j]<<" - "<<p->shared[j]<<" - "<< p->rms_t[j]<< " - " << p->view[j]<<") ";
02348         
02349         /*
02350         s << "spoints (e - shared) : ";
02351                 for(int j=0;j<p->entries;j++)
02352                         s << "(" <<p->e[j]<<" - "<<p->shared[j]<<") ";
02353         
02354         */
02355         
02356         
02357 
02358                         
02359                         
02360 
02361                         
02362                         
02363                 s<<endl<<endl;
02364                 
02365                 
02366         
02367         
02368         
02369 
02370         MSG("PrimaryShowerFinder",Msg::kDebug) <<s.str();
02371 
02372         
02373 }

void PrimaryShowerFinder::ExpandShowerChain ( Chain chain,
int  view,
double  start_z,
double  end_z 
) [private]

Definition at line 2051 of file PrimaryShowerFinder.cxx.

References cluster_manager, Managed::ManagedCluster::e, Chain::end_t, Chain::end_z, Managed::ClusterManager::FindClustersInZ(), Managed::ClusterManager::GetCluster(), Managed::ManagedCluster::id, Chain::insert(), Chain::start_t, Chain::start_z, Managed::ManagedCluster::t, and Managed::ManagedCluster::z.

02052 {
02053 
02054         if(chain->start_z < start_z && chain->end_z > end_z)return;
02055         
02056 //      printf("expanding shower chain from %f to %f\n",chain->start_z,chain->end_z);
02057                                 
02058         if(chain->start_z > start_z)
02059         {
02060                 double z=chain->start_z-0.06;
02061                 while(z>start_z-0.04)
02062                 {
02063                         //find hits in this z range and add the best one to the chain
02064                         std::vector<int> cs = cluster_manager->FindClustersInZ(z,view);
02065         
02066         //              printf("found %d hits in z %f\n",cs.size(),z);
02067                         
02068                         int closest=-1;
02069                         double closest_d=1000;
02070                         //for now do a simple straight line from start t
02071                         for(unsigned int i=0;i<cs.size();i++)   
02072                         {
02073                                 Managed::ManagedCluster *c = cluster_manager->GetCluster(cs[i]);
02074         //                      printf("id %d z %f t %f e %f\n",cs[i],c->z,c->t,c->e);
02075                         
02076                                 if(fabs(c->t-chain->start_t)<closest_d)
02077                                 {
02078                                         closest_d=fabs(c->t-chain->start_t);
02079                                         closest=c->id;
02080                                 }
02081                         }
02082                         
02083                         if(closest>-1 && closest_d < 0.1)
02084                         {
02085                                 Managed::ManagedCluster *c = cluster_manager->GetCluster(closest);
02086                                 chain->insert(c->t,c->z,c->e,c->id);
02087                         }
02088         
02089                         z-= 0.06; //subtract off a plane
02090                 }
02091         }
02092 
02093 
02094                 
02095         if(chain->end_z < end_z)
02096         {
02097                 double z=chain->end_z+0.06;
02098                 while(z<end_z+0.04)
02099                 {
02100                         //find hits in this z range and add the best one to the chain
02101                         std::vector<int> cs = cluster_manager->FindClustersInZ(z,view);
02102         
02103         //              printf("found %d hits in z %f\n",cs.size(),z);
02104                         
02105                         int closest=-1;
02106                         double closest_d=1000;
02107                         //for now do a simple straight line from start t
02108                         for(unsigned int i=0;i<cs.size();i++)   
02109                         {
02110                                 Managed::ManagedCluster *c = cluster_manager->GetCluster(cs[i]);
02111         //                      printf("id %d z %f t %f e %f\n",cs[i],c->z,c->t,c->e);
02112                         
02113                                 if(fabs(c->t-chain->end_t)<closest_d)
02114                                 {
02115                                         closest_d=fabs(c->t-chain->end_t);
02116                                         closest=c->id;
02117                                 }
02118                         }
02119                         
02120                         if(closest>-1  && closest_d < 0.1)
02121                         {
02122                                 Managed::ManagedCluster *c = cluster_manager->GetCluster(closest);
02123                                 chain->insert(c->t,c->z,c->e,c->id);
02124                         }
02125         
02126                         z+= 0.06; //subtract off a plane
02127                 }
02128         }
02129 
02130 
02131 }

Chain * PrimaryShowerFinder::ExtractMuon ( Chain c  )  [private]

Definition at line 1318 of file PrimaryShowerFinder.cxx.

References Chain::children, Chain::cluster_id, Chain::e, Chain::entries, Chain::insert(), Msg::kDebug, Chain::level, MSG, Chain::parentChain, Chain::t, and Chain::z.

Referenced by FindPrimaryShower().

01319 {
01320         if(!c)return 0;
01321         
01322         int required_consecutive_muon_hits=2;
01323         
01324         if(c->entries<required_consecutive_muon_hits)return 0;
01325         
01326         MSG("PrimaryShowerFinder",Msg::kDebug)<<"extract muon\n";
01327         
01328         //do we have a muon in this chain?
01329         
01330         double muoncount=0;
01331         double muonenergy=0;
01332         double sume=0;
01333         
01334         for(int i=0;i<c->entries;i++)
01335         {
01336                 if(c->e[i]<2.5)
01337                 {
01338                         muoncount++;
01339                         muonenergy+=c->e[i];
01340                 }
01341                 sume+=c->e[i];
01342         
01343         }
01344 
01345         double muonfrac = muoncount / c->entries;
01346 
01347         //double efrac=muonenergy/sume;
01348         
01349         MSG("PrimaryShowerFinder",Msg::kDebug)<<"mf "<<muonfrac<<"\n";
01350         if(muonfrac<0.2 ) return 0;//no muon
01351         
01352         
01353         
01354         //find where the muon exists the other particles... require n in a row muon hits
01355         
01356         int startmuon=0;
01357         int nmuon=0;
01358         
01359         for(int i=0;i<c->entries;i++)
01360         {
01361         
01362         //      printf("%d %f\n",i,c->e[i]);
01363                 startmuon=i;
01364                 if(c->e[i]<2.5)
01365                 {
01366                         nmuon++;
01367                         if(nmuon>=required_consecutive_muon_hits)break;  //require n muon hits in a row
01368                 
01369                 }else
01370                 {
01371                         nmuon=0;
01372                 }
01373         }
01374         
01375         if(nmuon<required_consecutive_muon_hits)return 0;//not enough muon like!
01376         //printf("sm %d\n",startmuon);
01377         startmuon=startmuon-nmuon+1;
01378         
01379         MSG("PrimaryShowerFinder",Msg::kDebug)<<"muon starts at "<<startmuon<<"\n";
01380         
01381         //determine average muon energy in this particle
01382         //using the muon hits, but ignoring large ones which could be muon->e shower
01383         
01384         double muone=0;
01385         muoncount=0;
01386         for(int i=startmuon;i<c->entries;i++)
01387         {
01388                 if(c->e[i]<2.5)
01389                 {
01390                         muone+=c->e[i];
01391                         muoncount++;
01392                 }
01393         }
01394         
01395         double adjmuon = muone/muoncount;
01396         
01397         MSG("PrimaryShowerFinder",Msg::kDebug)<<"avg muon e "<<adjmuon<<"\n";
01398         
01399         Chain ctemp;
01400         
01401         std::vector<double> savedE;
01402         for(int i=0;i<startmuon;i++)
01403         {
01404                 if(c->e[i]-adjmuon>0)
01405                 {
01406                         ctemp.insert(c->t[i],c->z[i],c->e[i]-adjmuon,c->cluster_id[i]);
01407                         savedE.push_back(adjmuon);
01408                 }else{
01409                         savedE.push_back(0);
01410                 }
01411                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"AAAAAAAAAAAA "<<c->t[i]<<" "<<c->z[i]<<" e "<<c->e[i]-adjmuon<<"\n";
01412         
01413         }
01414         
01415 
01416         Chain *muonchain = new Chain();
01417         muonchain->level=c->level; //so we get the parend id, etc
01418         muonchain->parentChain=c->parentChain;
01419         muonchain->children=c->children;
01420         
01421         
01422         for(int i=0;i<startmuon;i++)
01423         {
01424                 if(savedE[i]>0)
01425                         muonchain->insert(c->t[i],c->z[i],savedE[i],c->cluster_id[i]);
01426         }       
01427         
01428         for(int i=startmuon;i<c->entries;i++)
01429         {
01430                 muonchain->insert(c->t[i],c->z[i],c->e[i],c->cluster_id[i]);
01431         }               
01432 
01433 
01434 
01435         *c=ctemp;
01436         
01437         return muonchain;
01438 
01439 }

int PrimaryShowerFinder::FindFirstIntersection ( int  view,
double &  t,
double &  z 
)

Definition at line 979 of file PrimaryShowerFinder.cxx.

References cluster_manager, HoughLine::end_z, HoughLine::GetExpectedT(), houghlinesU, houghlinesV, intU, intV, Managed::ClusterManager::maxu, Managed::ClusterManager::maxv, Managed::ClusterManager::maxz, Managed::ClusterManager::minu, Managed::ClusterManager::minv, Managed::ClusterManager::minz, HoughLine::ncluster, HoughLine::offset_t, HoughLine::offset_z, HoughLine::r, HoughLine::start_z, and HoughLine::theta.

Referenced by HoughView::DrawHough().

00980 {
00981         std::vector<HoughLine> * hl;
00982         if(view==2)hl=&houghlinesU;
00983         else if(view==3)hl=&houghlinesV;
00984         else return 0;
00985 
00986         TH2D * inth = view==2?intU:intV;
00987         
00988         if(inth)delete inth;
00989         inth=0;
00990         
00991         double min_t = view==2?cluster_manager->minu:cluster_manager->minv;
00992         double max_t = view==2?cluster_manager->maxu:cluster_manager->maxv;
00993         inth = new TH2D("","",50,cluster_manager->minz,cluster_manager->maxz,50,min_t,max_t);
00994         inth->SetDirectory(0);
00995 
00996         double best_z = 10000;
00997         double best_t=0;
00998 
00999         //iterate over each line
01000         for(unsigned int i=0;i<(*hl).size();i++)
01001         {
01002                 //nested loop
01003                 for(unsigned int j=0;j<(*hl).size();j++)
01004                 {
01005                         if(i==j)continue;
01006                 
01007                         //computer intersection and save if better
01008                         HoughLine *h1 = &(*hl)[i];
01009                         HoughLine *h2 = &(*hl)[j];
01010                         
01011                         
01012                         z=
01013                         (
01014                                 h2->r/sin(h2->theta) + h2->offset_t 
01015                                 - h1->r/sin(h1->theta) - h1->offset_t 
01016                                 +(cos(h2->theta)/sin(h2->theta))*h2->offset_z 
01017                                 -(cos(h1->theta)/sin(h1->theta))*h1->offset_z
01018                         ) 
01019                                 /
01020                         (
01021                                 cos(h2->theta)/sin(h2->theta)
01022                                 -
01023                                 cos(h1->theta)/sin(h1->theta)
01024                         );
01025                         
01026                         if(z<best_z && z>0)
01027                         {
01028                                 best_z=z;
01029                                 best_t= h1->GetExpectedT(best_z);                       
01030                         }
01031                         
01032                         //require the intersection to be within 2 planes of the start/end of either houghline
01033                         if(fabs(h1->start_z-z)>0.07 && fabs(h2->start_z-z)>0.07 && fabs(h1->end_z-z)>0.07 && fabs(h2->end_z-z)>0.07)continue;
01034                         
01035                                 inth->Fill(z,h1->GetExpectedT(z),h1->ncluster+h2->ncluster);
01036         
01037                 }
01038         }
01039         
01040         inth->Draw("colz");
01041         if(best_z >=10000)return 0;
01042 
01043 
01044         t=best_t;
01045         z=best_z;
01046         return 1;
01047 }

void PrimaryShowerFinder::FindHoughLineMatches (  )  [private]

Definition at line 3253 of file PrimaryShowerFinder.cxx.

References HoughLine::cluster_id, cluster_manager, Managed::ManagedCluster::e, HoughLine::end_t, HoughLine::end_z, Managed::ClusterManager::GetCluster(), houghlineMatch, houghlinesU, houghlinesV, HoughLine::ncluster, HoughLine::phi, HoughLine::start_t, HoughLine::start_z, and HoughLine::sum_e.

03254 {
03255         houghlineMatch.clear();
03256 
03257         //need to iterate over the hough lines in both view and look for matches
03258 
03259         for(unsigned int i=0;i<houghlinesU.size();i++)
03260         {
03261                 //don't want to match to verticle lines
03262                 if(fabs(houghlinesU[i].phi -3.141592/2)<0.1)continue;
03263                 if(houghlinesU[i].ncluster<2)continue;
03264                         HoughLine * l = &houghlinesU[i];
03265 
03266                 for(unsigned int j=0;j<houghlinesV.size();j++)
03267                 {
03268                         //don't want to match to verticle lines
03269                         if(fabs(houghlinesV[j].phi - 3.141592/2)<0.1)continue;
03270                         if(houghlinesV[j].ncluster<2)continue;
03271 
03272                         HoughLine * r = &houghlinesV[j];
03273                         
03274                         double overlap = l->end_z<r->end_z?l->end_z:r->end_z;
03275                         overlap -= l->start_z > r->start_z ? l->start_z:r->start_z;
03276                         
03277                         //require overlap to exceed 90% of the length of the smaller line
03278                         double l_z = l->end_z-l->start_z;
03279                         double r_z = r->end_z-r->start_z;
03280                         double smaller_z = l_z < r_z?l_z:r_z;
03281                         if(smaller_z * 0.5 > overlap)
03282                         {
03283                         //      printf("failing %d %d with overlap %f > %f\n",i,j,smaller_z * 0.5,overlap);
03284                                 continue;
03285                         }
03286                         double smaller_e = l->sum_e < r->sum_e ? l->sum_e:r->sum_e;
03287                         double larger_e = l->sum_e > r->sum_e ? l->sum_e:r->sum_e;
03288 
03289                         //require smaller e to be more than 30% of larger e
03290                         if(larger_e*0.3 > smaller_e)
03291                         {
03292                 //              printf("failing %d %d with energy %f > %f\n",i,j,larger_e*0.5 , smaller_e);
03293                                 continue;
03294                         }
03295 
03296 
03297                         //require a peak hit that is at least 20% of total e
03298                 
03299 
03300                         
03301                         double max_hit_e=0;
03302                         double sum_hit_e=0;
03303                         for(unsigned int ka=0;ka<l->cluster_id.size();ka++)
03304                         {
03305                                 Managed::ManagedCluster * clus = cluster_manager->GetCluster(l->cluster_id[ka]);
03306                                 if(!clus){
03307                                         //printf("missing cluster %d\n",l->cluster_id[ka]);
03308                                         continue;
03309                                 }
03310                                 sum_hit_e+=clus->e;
03311                                 
03312                                 if(clus->e>max_hit_e)max_hit_e=clus->e;
03313                         
03314                         }
03315                         for(unsigned int ka=0;ka<r->cluster_id.size();ka++)
03316                         {
03317                                 Managed::ManagedCluster * clus = cluster_manager->GetCluster(r->cluster_id[ka]);
03318                                 if(!clus){
03319                                         //printf("missing cluster %d\n",r->cluster_id[ka]);
03320                                         continue;
03321                                 }                                       
03322                                 sum_hit_e+=clus->e;
03323                                 
03324                                 if(clus->e>max_hit_e)max_hit_e=clus->e;
03325                         
03326                         }       
03327                         if(sum_hit_e*0.2 > max_hit_e)
03328                         {
03329                         //      printf("continuing sum hit*.3 %f > maxhit %f\n",sum_hit_e*0.3, max_hit_e);
03330                                 continue;
03331                         }
03332                         if( max_hit_e<5)
03333                         {
03334                         //      printf("maxhit too small at %f\n", max_hit_e);
03335                                 continue;
03336                         }
03337                         
03338                 //      printf("saving %d %d\n",i,j);
03339                         houghlineMatch.push_back(make_pair(i,j));                                       
03340                 }
03341         }
03342 
03343         std::map<double,int> orderer;
03344         for(unsigned int i=0;i<houghlineMatch.size();i++)
03345         {
03346                 HoughLine * l = &houghlinesU[houghlineMatch[i].first];
03347                 HoughLine * r = &houghlinesV[houghlineMatch[i].second];
03348                 if(l->ncluster<2 || r->ncluster<2)continue;
03349                 double weight = (l->sum_e+r->sum_e)*(cos(l->phi))*(cos(r->phi));
03350                 weight *= l->ncluster / sqrt( (l->end_z-l->start_z)*(l->end_z-l->start_z)+(l->end_t-l->start_t)*(l->end_t-l->start_t));
03351                 weight *= r->ncluster / sqrt( (r->end_z-r->start_z)*(r->end_z-r->start_z)+(r->end_t-r->start_t)*(r->end_t-r->start_t));
03352 
03353                 orderer.insert(make_pair(weight,i));
03354         }               
03355 
03356         std::map<double,int>::iterator it_orderer;
03357 
03358         
03359         //printf("found %d preliminary matches\n",houghlineMatch.size());
03360 
03361         /*for(it_orderer = orderer.begin();it_orderer!=orderer.end();it_orderer++)
03362         //for(unsigned int i=0;i<houghlineMatch.size();i++)
03363         {
03364                 int i= it_orderer->second;
03365                 
03366                 //printf("match %d %d\n",houghlineMatch[i].first,houghlineMatch[i].second);
03367                 //HoughLine * l = &houghlinesU[houghlineMatch[i].first];
03368                 //HoughLine * r = &houghlinesV[houghlineMatch[i].second];
03369         
03370                 //printf("\t U s %f %f e %f %f sume %f clust %d phi %f\n",l->start_z,l->start_t, l->end_z,l->end_t,l->sum_e,l->ncluster,l->phi);
03371                 //printf("\t V s %f %f e %f %f sume %f clust %d phi %f\n",r->start_z,r->start_t, r->end_z,r->end_t,r->sum_e,r->ncluster,r->phi);
03372         }
03373 */
03374 }

int PrimaryShowerFinder::FindPrimaryShower ( Managed::ClusterManager cm  ) 

!!! we also need to remake the particle here after adjusting chain energy! ///////////

Definition at line 1446 of file PrimaryShowerFinder.cxx.

References ChainHelper::AddFinishedChain(), Chain::available, Bundle(), Particle3D::calibrated_energy, chain_u, chain_v, CheckChainQuality(), ShwFit::chisq, chu, chv, CleanHoughLines(), Chain::cluster_id, cluster_manager, Particle3D::cmp_chisq, ShwFit::cmp_chisq, ShwFit::cmp_ndf, Particle3D::cmp_ndf, ShwFit::conv, ChainHelper::DeleteChain(), DumpParticle3D(), Managed::ManagedCluster::e, Chain::e, Particle3D::e, Particle3D::electron, Particle3D::emfit_a, Particle3D::emfit_a_err, Particle3D::emfit_b, Particle3D::emfit_b_err, Particle3D::emfit_chisq, Particle3D::emfit_e0, Particle3D::emfit_e0_err, Particle3D::emfit_ndf, Particle3D::emfit_prob, Chain::end_z, Chain::entries, Particle3D::entries, ExtractMuon(), ChainHelper::FindMaxPath(), ChainHelper::finished, ShwFit::Fit3d(), foundparticle, foundparticle3d, ChainHelper::GetChain(), Managed::ClusterManager::GetCluster(), Chain::insert(), ShwFit::Insert3d(), Msg::kDebug, MakeChains(), MakeHoughMap(), MakeParticle3D(), MSG, Chain::myId, ShwFit::ndf, Particle3D::other, ShwFit::par_a, ShwFit::par_a_err, ShwFit::par_b, ShwFit::par_b_err, ShwFit::par_e0, ShwFit::par_e0_err, Particle3D::particletype, Particle3D::peakdiff, ShwFit::peakdiff, ShwFit::post_over, Particle3D::post_over, Particle3D::post_under, ShwFit::post_under, ShwFit::pp_chisq, Particle3D::pp_chisq, ShwFit::pp_igood, Particle3D::pp_igood, Particle3D::pp_ndf, ShwFit::pp_ndf, ShwFit::pp_p, Particle3D::pp_p, Particle3D::pre_over, ShwFit::pre_over, ShwFit::pre_under, Particle3D::pre_under, Particle3D::pred_b, ShwFit::pred_b, Particle3D::pred_e0, ShwFit::pred_e0, Particle3D::pred_e_a, ShwFit::pred_e_a, ShwFit::pred_e_chisq, Particle3D::pred_e_chisq, ShwFit::pred_e_ndf, Particle3D::pred_e_ndf, Particle3D::pred_g_a, ShwFit::pred_g_a, Particle3D::pred_g_chisq, ShwFit::pred_g_chisq, Particle3D::pred_g_ndf, ShwFit::pred_g_ndf, ChainHelper::print_finished(), Chain::PrintChain(), ShwFit::prob, ran, Chain::Recalc(), Reset(), ShwFit::Reset(), Managed::ClusterManager::SaveCluster(), single_view_long_shower, Chain::start_z, Particle3D::start_z, Chain::t, Particle3D::u, Particle3D::v, Particle3D::z, Chain::z, and ShwFit::zshift.

Referenced by Finder::Process().

01447 {
01448 
01449         MSG("PrimaryShowerFinder",Msg::kDebug)<<"looking for primary shower\n";
01450         Reset(0); //keep any vertex that we might have set
01451 
01452         ran=1;
01453 
01454         cluster_manager=cm;
01455         
01456         
01457         MakeHoughMap(2);
01458         MakeHoughMap(3);
01459 
01460 //DumpHoughLines(2);
01461 //DumpHoughLines(3);
01462 
01463         CleanHoughLines(2);
01464         CleanHoughLines(3);
01465         
01466         Bundle(2);
01467         Bundle(3);
01468 
01469         MakeChains(2);
01470         MakeChains(3);
01471 
01472 
01473 //      FindHoughLineMatches();
01474         
01475         
01476         //if(houghlineMatch.size()<1)
01477                 
01478         
01479         
01480         //try to find a long Shower chain in each view
01481         //chain_u = FindShowerChain(cm,2);
01482         //chain_v = FindShowerChain(cm,3);
01483         
01484         
01485 //      chain_u = MakeShowerChain(2);
01486 //      chain_v = MakeShowerChain(3);
01487 
01488 
01489         chain_u =0;
01490         chain_v=0;
01491         
01492         for(unsigned int i=0;i<chu->finished.size();i++)
01493         {
01494                 if(chu->finished[i].available)
01495                 {
01496                         chain_u=&chu->finished[i];
01497                         break;
01498                 }
01499         }
01500         
01501         for(unsigned int i=0;i<chv->finished.size();i++)
01502         {
01503                 if(chv->finished[i].available)
01504                 {
01505                         chain_v=&chv->finished[i];
01506                         break;
01507                 }
01508         }
01509 
01510         if(!chain_u || !chain_v)return 0;
01511         
01512 
01513 //make temp chains for use in the finder which represent a max path
01514         
01515 
01516 
01517         //printf("building chain\n");
01518         chu->print_finished();  
01519         std::pair< std::vector<int>, double> path;
01520         path = chu->FindMaxPath(chain_u);       
01521         chain_u=new Chain();
01522 
01523 
01524         //printf("chain made from ");
01525 //      for(int i=0 ;i <path.first.size();i++)
01526 //              printf("%d ",path.first[i]);
01527 //      printf("\n");
01528                 
01529 
01530 
01531         for(int i=0 ;i <(int)path.first.size();i++)
01532         {
01533                 Chain *c = chu->GetChain(path.first[i]);
01534                 if(!c)continue;
01535                 for(int j=0;j<c->entries;j++)
01536                 {
01537                         chain_u->insert(c->t[j],c->z[j],c->e[j],c->cluster_id[j]);
01538                         MSG("PrimaryShowerFinder",Msg::kDebug)<<"u "<<c->z[j]<<" "<<c->t[j]<<" "<<c->e[j]<<"\n";
01539                 }
01540         }
01541         std::vector<int> upath=path.first;
01542 
01543         //printf("building chain\n");
01544         chv->print_finished();
01545         path = chv->FindMaxPath(chain_v);       
01546         chain_v=new Chain();
01547         
01548         //printf("chain made from ");
01549         //      for(int i=0 ;i <path.first.size();i++)
01550         //              printf("%d ",path.first[i]);
01551         //      printf("\n");
01552                 
01553 
01554         for(int i=0 ;i <(int)path.first.size();i++)
01555         {
01556                 Chain *c = chv->GetChain(path.first[i]);
01557                 if(!c)continue;
01558 //              printf("chain %d\n",path.first[i]);
01559                 for(int j=0;j<c->entries;j++)
01560                 {
01561 //                      printf("id %d\n",c->cluster_id[j]);
01562                         chain_v->insert(c->t[j],c->z[j],c->e[j],c->cluster_id[j]);
01563                         MSG("PrimaryShowerFinder",Msg::kDebug)<<"v "<<c->z[j]<<" "<<c->t[j]<<" "<<c->e[j]<<"\n";
01564                 }
01565         }
01566         std::vector<int> vpath=path.first;              
01567         
01568 /*      
01569         printf("\n\n---------\nchain\n");
01570         for(int i=0;i<chain_u->entries;i++)
01571         {
01572                 printf("%d %f %f %f\n",chain_u->cluster_id[i],chain_u->z[i],chain_u->t[i],chain_u->e[i]);
01573         }
01574 
01575         printf("\n\n---------\nchain\n");
01576         for(int i=0;i<chain_v->entries;i++)
01577         {
01578                 printf("%d %f %f %f\n",chain_v->cluster_id[i],chain_v->z[i],chain_v->t[i],chain_v->e[i]);
01579         }       
01580 */      
01581         //muon extraction:
01583         Chain * mchainU=0;
01584         Chain * mchainV=0;
01585         mchainU=ExtractMuon(chain_u);
01586         mchainV=ExtractMuon(chain_v);
01588 /*      
01589         printf("\n\n---------\nmuon chain\n");
01590         if(mchainU)for(int i=0;i<mchainU->entries;i++)
01591         {
01592                 printf("%d %f %f %f\n",mchainU->cluster_id[i],mchainU->z[i],mchainU->t[i],mchainU->e[i]);
01593         }
01594 
01595         printf("\n\n---------\nmuon chain\n");
01596         if(mchainV)for(int i=0;i<mchainV->entries;i++)
01597         {
01598                 printf("%d %f %f %f\n",mchainV->cluster_id[i],mchainV->z[i],mchainV->t[i],mchainV->e[i]);
01599         }       
01600         
01601 
01602         printf("\n\n---------\nchain\n");
01603         for(int i=0;i<chain_u->entries;i++)
01604         {
01605                 printf("%d %f %f %f\n",chain_u->cluster_id[i],chain_u->z[i],chain_u->t[i],chain_u->e[i]);
01606         }
01607 
01608         printf("\n\n---------\nchain\n");
01609         for(int i=0;i<chain_v->entries;i++)
01610         {
01611                 printf("%d %f %f %f\n",chain_v->cluster_id[i],chain_v->z[i],chain_v->t[i],chain_v->e[i]);
01612         }       
01613         */
01614         
01615         //double ustart=chain_u->start_z;
01616         //double vstart=chain_v->start_z;
01617         //double uend=chain_u->end_z;
01618         //double vend=chain_v->end_z;
01619 //      ExpandShowerChain(chain_u,2,vstart,vend);
01620 //      ExpandShowerChain(chain_v,3,ustart,uend);
01621         
01622         
01623         int qual_u=0;
01624         int qual_v=0;
01625         
01626         if(chain_u)
01627         {
01628                 chain_u->Recalc();
01629                 qual_u=CheckChainQuality(chain_u);
01630         }
01631         if(chain_v)
01632         {
01633                 chain_v->Recalc();
01634                 qual_v=CheckChainQuality(chain_v);
01635         }
01636 
01637         if(qual_u&&!qual_v)
01638         {
01639                 if(chain_u->end_z - chain_u->start_z > 2)single_view_long_shower=1;
01640         }else if(qual_v&&!qual_u)
01641         {
01642                 if(chain_v->end_z - chain_v->start_z > 2)single_view_long_shower=1;
01643         }
01644         
01645         if(!qual_u || !qual_v)return 0;
01646 
01647 
01648 
01649         //should do something smarter... like try other chains...
01650         
01651         //if we have a quality chain in each view... then we will be making a particle
01652 /*      if(qual_u && qual_v)
01653         {
01654         
01655                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"qqqqq "<<qual_u <<" "<< qual_v<<"\n";
01656                 
01657         
01658                 //look in both views to find at least 3u and 3v consecutive planes with only 1 strip... that is where we say the Shower exits the rest of the shower/partices/etc
01659 
01660                 int qual=CheckChainOverlap(chain_u, chain_v);
01661                 if(qual)
01662                 {
01663                 
01664                 ClearFrontVertex(chain_u, chain_v);
01665 */
01666 
01667                 
01668                 
01669         
01670         
01671                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"making particle\n";
01672                 MakeParticle3D();
01673                 if(foundparticle3d)
01674                 {
01675                         foundparticle=1;
01676                 
01677                         foundparticle3d->particletype=Particle3D::electron;
01678                 
01679                         MSG("PrimaryShowerFinder",Msg::kDebug)<<"particle made\n";
01680                         DumpParticle3D();
01681                         
01682 //                      chain_u->available=0;
01683 //                      chain_v->available=0;
01684 
01685                         for(int i=0 ;i <(int)upath.size();i++)
01686                         {
01687                                 Chain *c = chu->GetChain(upath[i]);
01688                                 
01689                                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"u chain "<<upath[i]<<"\n";
01690                                 
01691                                 c->available=0;
01692                         }
01693                         for(int i=0 ;i <(int)vpath.size();i++)
01694                         {
01695                                 Chain *c = chv->GetChain(vpath[i]);
01696                                 
01697                                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"v chain "<<vpath[i]<<"\n";
01698                                 
01699                                 c->available=0;
01700                         }                       
01701                         
01702                         //do a fit on the particle3d to see if it is like an electron..
01703                         Particle3D *p=foundparticle3d;          
01704                         ShwFit f;
01705                         
01706                         double startz=0;
01707                         
01708                         startz=p->z[0];
01709                         
01710                         for(int i=0;i<p->entries;i++)
01711                         {
01712                                 //double planewidth=0.035;
01713                                 
01714                 //              f.Insert(p->e[i],(int)((p->z[i]-startz)/planewidth));
01715                         
01716                 //              cout <<"inserting "<<p->e[i]<<" at plane "<<(int)((p->z[i]-startz)/planewidth)<<endl;
01717                 
01718                 //for now, do the simple thing and assume each hit is in the next plane...
01719                         
01720                                 //                      f.Insert(p->e[i],i+1);
01721                         
01722                         
01723                                 f.Insert3d(p->e[i],p->u[i],p->v[i],p->z[i]);
01724                         
01725                                 //cout <<"inserting "<<p->e[i]<<" at plane "<<i+1<<endl;
01726                 
01727                         
01728                         }
01729                         
01730                         //f.Fit();
01731                 
01732 
01733                         f.Fit3d(1);
01734 
01735 
01736 
01737 
01738                         MSG("PrimaryShowerFinder",Msg::kDebug)<<"EM FIT! " << f.par_a << " " << f.par_b << " +- "<<f.par_b_err<<" "<<f.par_e0<<endl;
01739                 
01740                         p->emfit_a=f.par_a;
01741                         p->emfit_b=f.par_b;
01742                         p->emfit_e0=f.par_e0;
01743                         p->emfit_a_err=f.par_a_err;
01744                         p->emfit_b_err=f.par_b_err;
01745                         p->emfit_e0_err=f.par_e0_err;
01746                         p->emfit_prob=f.prob;
01747                         p->calibrated_energy=f.par_e0;
01748                         p->emfit_chisq=f.chisq;
01749                         p->emfit_ndf=f.ndf;
01750 
01751                         p->pred_e_a=f.pred_e_a; 
01752                         p->pred_g_a=f.pred_g_a;
01753                         p->pred_b=f.pred_b;
01754                         p->pred_e0=f.pred_e0;
01755                         p->pred_e_chisq=f.pred_e_chisq;
01756                         p->pred_e_ndf=f.pred_e_ndf;
01757                         p->pred_g_chisq=f.pred_g_chisq;
01758                         p->pred_g_ndf=f.pred_g_ndf;                             
01759                         
01760                         p->pre_over=f.pre_over;         
01761                         p->pre_under=f.pre_under;               
01762                         p->post_over=f.post_over;               
01763                         p->post_under=f.post_under;             
01764 
01765                         p->pp_chisq=f.pp_chisq;
01766                         p->pp_ndf=f.pp_ndf;
01767                         p->pp_igood=f.pp_igood;
01768                         p->pp_p=f.pp_p;
01769         
01770                         p->cmp_chisq = f.cmp_chisq;
01771                         p->cmp_ndf = f.cmp_ndf;
01772                         p->peakdiff = f.peakdiff;       
01773                         
01774                         int redochain=0; 
01775                         //do we need to redo the chain because we are shifting the start past the first hits of the chain?
01776                         if(f.zshift)
01777                         {
01778                                 if(f.zshift < 0)redochain=1;
01779                                 p->start_z -=f.zshift;
01780                                 
01781                                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"zshift "<<f.zshift<<"\n";
01782                         }
01783                         
01784                         
01785                         if(!f.conv)p->particletype=Particle3D::other;
01786                         else{
01787                         
01788                         
01789                         
01790                         
01791                         //redo the chain if needed
01792                         if(redochain)
01793                         {
01794 
01795                                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"REDOING CHAINS\n";
01796                         
01797                                 Chain temp_u;
01798                                 for(int i=0;i<chain_u->entries;i++)
01799                                 {
01800                                         if(chain_u->z[i]<p->start_z){
01801                                                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"throwing away chain hit with z "<<chain_u->z[i]<<"\n";
01802                                                 continue;
01803                                         }
01804                                         temp_u.insert(chain_u->t[i],chain_u->z[i],chain_u->e[i],chain_u->cluster_id[i]);
01805                                 }
01806                                 temp_u.myId=chain_u->myId;
01807                                 *chain_u=temp_u;                        
01808 
01809                                 Chain temp_v;
01810                                 for(int i=0;i<chain_v->entries;i++)
01811                                 {
01812                                         if(chain_v->z[i]<p->start_z){
01813                                                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"throwing away chain hit with z "<<chain_v->z[i]<<"\n";
01814                                                 continue;
01815                                         }
01816                                         temp_v.insert(chain_v->t[i],chain_v->z[i],chain_v->e[i],chain_v->cluster_id[i]);
01817                                 }
01818                                 temp_v.myId=chain_v->myId;              
01819                                 *chain_v=temp_v;
01820                                 chain_u->available=0;
01821                                 chain_v->available=0;
01822 
01823                                 delete foundparticle3d; 
01824                                 foundparticle=0;        
01825                                 foundparticle3d=0;
01826                                         
01827                                 //remake particle
01828                                 MakeParticle3D();
01829 
01830                                 if(!foundparticle3d)return 0;
01831 
01832                                 if(foundparticle3d)
01833                                 {
01834                                 foundparticle=1;        
01835                                 p=foundparticle3d;      
01836                                 
01837                                 f.Reset();
01838                                 
01839                                 double startz=0;
01840                         
01841                                 startz=p->z[0];
01842                         
01843                                 for(int i=0;i<p->entries;i++)
01844                                 {
01845                                         //double planewidth=0.035;
01846                                         f.Insert3d(p->e[i],p->u[i],p->v[i],p->z[i]);
01847                                 }
01848                         
01849 
01850                                 
01851                                 foundparticle3d->particletype=Particle3D::electron;
01852                                 f.Fit3d();
01853                         
01854                                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"EM FIT! " << f.par_a << " " << f.par_b << " +- "<<f.par_b_err<<" "<<f.par_e0<<endl;
01855                 
01856                                 p->emfit_a=f.par_a;
01857                                 p->emfit_b=f.par_b;
01858                                 p->emfit_e0=f.par_e0;
01859                                 p->emfit_a_err=f.par_a_err;
01860                                 p->emfit_b_err=f.par_b_err;
01861                                 p->emfit_e0_err=f.par_e0_err;
01862                                 p->emfit_prob=f.prob;
01863                                 p->calibrated_energy=f.par_e0;
01864                                 p->emfit_chisq=f.chisq;
01865                                 p->emfit_ndf=f.ndf;
01866 
01867                                 p->pred_e_a=f.pred_e_a; 
01868                                 p->pred_g_a=f.pred_g_a;
01869                                 p->pred_b=f.pred_b;
01870                                 p->pred_e0=f.pred_e0;
01871                                 p->pred_e_chisq=f.pred_e_chisq;
01872                                 p->pred_e_ndf=f.pred_e_ndf;
01873                                 p->pred_g_chisq=f.pred_g_chisq;
01874                                 p->pred_g_ndf=f.pred_g_ndf;                             
01875                         
01876                                 p->pre_over=f.pre_over;         
01877                                 p->pre_under=f.pre_under;               
01878                                 p->post_over=f.post_over;               
01879                                 p->post_under=f.post_under;     
01880                                 
01881                                 p->pp_chisq=f.pp_chisq;
01882                                 p->pp_ndf=f.pp_ndf;
01883                                 p->pp_igood=f.pp_igood;
01884                                 p->pp_p=f.pp_p;
01885         
01886                                 p->cmp_chisq = f.cmp_chisq;
01887                                 p->cmp_ndf = f.cmp_ndf;
01888                                 p->peakdiff = f.peakdiff;       
01889                                                                 
01890                         
01891                                 if(!f.conv)p->particletype=Particle3D::other;
01892                                 if(f.zshift)
01893                                 {
01894                                         p->start_z -=f.zshift;
01895                                         MSG("PrimaryShowerFinder",Msg::kDebug)<<"zshift "<<f.zshift<<"\n";
01896                                 }
01897                         
01898                                 }
01899                         
01900                         }
01901                         
01902                         
01903                 }       
01904 
01905                         
01906                 //save clusters involved...
01907 
01908                 //don't save until we know that we want it!
01909 
01910                 for(int i=0;i<chain_u->entries;i++)
01911                 {
01912                         int newid = cluster_manager->SaveCluster(chain_u->cluster_id[i],chain_u->e[i],2);
01913                         if(!newid)continue;//if we have an issue saving, such as there is no energy...
01914                         Managed::ManagedCluster * newcluster = cluster_manager->GetCluster(newid);
01915                         if(fabs(chain_u->e[i]-newcluster->e)>0.001)
01916                                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"saving with different e "<<chain_u->e[i]<<" "<<newcluster->e<<"\n";
01917                         chain_u->e[i]=newcluster->e;
01918                         chain_u->cluster_id[i]=newid;
01919                 }
01920 
01921 
01922 
01923                 for(int i=0;i<chain_v->entries;i++)
01924                 {
01925                         int newid = cluster_manager->SaveCluster(chain_v->cluster_id[i],chain_v->e[i],2);
01926                         if(!newid)continue;//if we have an issue saving, such as there is no energy...
01927                         Managed::ManagedCluster * newcluster = cluster_manager->GetCluster(newid);
01928                         
01929                         //printf("saving cluster id %d e %f newcluster ok ? %d\n",chain_v->cluster_id[i],chain_v->e[i],newcluster!=0);
01930                         //if(newcluster)printf("new cluster id %d  actual %d\n",newid,newcluster->id);
01931                         if(fabs(chain_v->e[i]-newcluster->e)>0.001)
01932                                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"saving with different e "<<chain_v->e[i]<<" "<<newcluster->e<<"\n";
01933                         chain_v->e[i]=newcluster->e;
01934                         chain_v->cluster_id[i]=newid;
01935                 }               
01936 
01937                         
01938                 //save chains involved!
01939                 chain_u->available=0;
01940                 chain_v->available=0;
01941 
01942                 //printf("psf wants to add %d %d\n",chain_u->myId,chain_v->myId);
01943                 chu->AddFinishedChain(*chain_u);
01944                 chv->AddFinishedChain(*chain_v);
01945                 
01946 
01947                 
01948                 //did we have muon chains?
01949                 if(mchainU)chu->AddFinishedChain(*mchainU);     
01950                 if(mchainV)chv->AddFinishedChain(*mchainV);             
01951                         
01952                 //we made the electron chain temporary....
01953                 //now remove the actual chains that were used to make the electron and muon extracted chains
01954                 
01955                 
01956                 for(unsigned int i=0;i<upath.size();i++)
01957                 {
01958                         chu->DeleteChain(upath[i]);
01959                 }
01960  
01961                 for(unsigned int i=0;i<vpath.size();i++)
01962                 {
01963                         chv->DeleteChain(vpath[i]);
01964                 }
01965                                 
01966 
01967                 
01970                                                 //remake particle
01971  
01972                                 delete foundparticle3d; 
01973                                 foundparticle=0;        
01974                                 foundparticle3d=0;
01975                                         
01976                                 //remake particle
01977                                 MakeParticle3D();
01978                                 p=foundparticle3d;      
01979                                 if(!foundparticle3d)return 0;
01980 
01981                                 foundparticle=1;
01982                                                                 
01983                                 foundparticle3d->particletype=Particle3D::electron;
01984                                 f.Fit3d(1);
01985                         
01986                                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"EM FIT! " << f.par_a << " " << f.par_b << " +- "<<f.par_b_err<<" "<<f.par_e0<<endl;
01987                 
01988                                 p->emfit_a=f.par_a;
01989                                 p->emfit_b=f.par_b;
01990                                 p->emfit_e0=f.par_e0;
01991                                 p->emfit_a_err=f.par_a_err;
01992                                 p->emfit_b_err=f.par_b_err;
01993                                 p->emfit_e0_err=f.par_e0_err;
01994                                 p->emfit_prob=f.prob;
01995                                 p->calibrated_energy=f.par_e0;
01996                                 p->emfit_chisq=f.chisq;
01997                                 p->emfit_ndf=f.ndf;
01998 
01999                                 p->pred_e_a=f.pred_e_a; 
02000                                 p->pred_g_a=f.pred_g_a;
02001                                 p->pred_b=f.pred_b;
02002                                 p->pred_e0=f.pred_e0;
02003                                 p->pred_e_chisq=f.pred_e_chisq;
02004                                 p->pred_e_ndf=f.pred_e_ndf;
02005                                 p->pred_g_chisq=f.pred_g_chisq;
02006                                 p->pred_g_ndf=f.pred_g_ndf;                             
02007                         
02008                                 p->pre_over=f.pre_over;         
02009                                 p->pre_under=f.pre_under;               
02010                                 p->post_over=f.post_over;               
02011                                 p->post_under=f.post_under;             
02012 
02013                                 p->pp_chisq=f.pp_chisq;
02014                                 p->pp_ndf=f.pp_ndf;
02015                                 p->pp_igood=f.pp_igood;
02016                                 p->pp_p=f.pp_p;
02017         
02018         
02019                                 p->cmp_chisq = f.cmp_chisq;
02020                                 p->cmp_ndf = f.cmp_ndf;
02021                                 p->peakdiff = f.peakdiff;       
02022                                 
02023                         
02024                                 if(!f.conv)p->particletype=Particle3D::other;
02025                 
02027                         
02028                         
02029                         }
02030                         
02031                 
02032 //              }
02033 //      }
02034 
02035 
02036 
02037         MSG("PrimaryShowerFinder",Msg::kDebug)<<"printing long Shower chains\nu\n";
02038         if(chain_u)chain_u->PrintChain();MSG("PrimaryShowerFinder",Msg::kDebug)<<"v\n";
02039         if(chain_v)chain_v->PrintChain();
02040         MSG("PrimaryShowerFinder",Msg::kDebug)<<"done\n";
02041         
02042 
02043         MSG("PrimaryShowerFinder",Msg::kDebug)<<"done looking for long Showers\n";
02044         return foundparticle;
02045 }

Chain * PrimaryShowerFinder::FindShowerChain ( Managed::ClusterManager cl,
int  view 
) [private]

Definition at line 2234 of file PrimaryShowerFinder.cxx.

References Managed::ManagedCluster::e, Managed::ClusterManager::GetCluster(), Managed::ClusterManager::GetClusterMap(), Chain::insert(), Msg::kDebug, MSG, Managed::ManagedCluster::t, and Managed::ManagedCluster::z.

02235 {
02236 
02237         MSG("PrimaryShowerFinder",Msg::kDebug)<<"\n\nLooking for long Shower Chain in view "<<view<<"\n";
02238 
02239 
02240         std::map<double, std::map<double, int> > * cluster_map = cl->GetClusterMap(view);
02241 
02242 
02243     std::map<double, std::map<double, int> >::reverse_iterator p_iterr;
02244     std::map<double, int >::iterator s_iterr;
02245 
02246 
02247         Chain * c = new Chain();
02248 
02249         Managed::ManagedCluster * largest=0;
02250         double last_plane=100000;
02251         
02252         std::vector<Managed::ManagedCluster *> maxplaneclusters;
02253         //first we want to find maximum hits in each plane
02254         //we will then check these hits to see if they can form a tightly aligned track
02255     for(p_iterr=cluster_map->rbegin();p_iterr!=cluster_map->rend(); p_iterr++)
02256     {   
02257         //new plane and a largest cluster in the previous plane?
02258                 if(fabs(last_plane-p_iterr->first)>0.03 && largest)
02259                 {
02260                         MSG("PrimaryShowerFinder",Msg::kDebug)<<"largest cluster found "<<largest->t<<" "<<largest->z<<" "<<largest->e<<"\n";
02261                         maxplaneclusters.push_back(largest);
02262                         largest=0;
02263                 }
02264         
02265         for(s_iterr=p_iterr->second.begin();s_iterr!=p_iterr->second.end(); s_iterr++)
02266         {
02267                         Managed::ManagedCluster * clus = cl->GetCluster(s_iterr->second);
02268                         if(!largest)
02269                         {
02270                                 largest=clus;
02271                                 continue;
02272                         }
02273                         if(largest->e<clus->e)largest=clus;
02274                 }
02275         }
02276         if(largest)
02277         {
02278                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"largest cluster found "<<largest->t<<" "<<largest->z<<" "<<largest->e<<"\n";
02279                 maxplaneclusters.push_back(largest);
02280                 largest=0;
02281         }
02282 
02283         
02284         for(unsigned int i=0;i<maxplaneclusters.size();i++)
02285                 c->insert(maxplaneclusters[i]->t,maxplaneclusters[i]->z,maxplaneclusters[i]->e,maxplaneclusters[i]->id);                
02286                         
02287         return c;                       
02288 
02289 
02290 }

int PrimaryShowerFinder::FoundSingleViewPrimaryShower (  )  [inline]

Definition at line 32 of file PrimaryShowerFinder.h.

References single_view_long_shower.

Referenced by Finder::Process().

00032 {return single_view_long_shower;};

Chain* PrimaryShowerFinder::GetChain ( int  view  )  [inline]

Definition at line 29 of file PrimaryShowerFinder.h.

References chain_u, chain_v, and foundparticle.

00029 {if(foundparticle){if(view==2)return chain_u;if(view==3)return chain_v;}return 0;};

std::vector<HoughLine>* PrimaryShowerFinder::GetHoughLines ( int  view  )  [inline]

Definition at line 37 of file PrimaryShowerFinder.h.

References houghlinesU, and houghlinesV.

Referenced by HoughView::DrawHough().

00037 {if(view==2)return &houghlinesU;if(view==3)return &houghlinesV; return 0;};

TH2D* PrimaryShowerFinder::GetHoughMap ( int  view  )  [inline]

Definition at line 36 of file PrimaryShowerFinder.h.

References houghmapU, and houghmapV.

00036 {if(view==2)return houghmapU; if(view==3)return houghmapV; return 0;};

Particle3D* PrimaryShowerFinder::GetParticle3D (  )  [inline]

Definition at line 28 of file PrimaryShowerFinder.h.

References foundparticle, and foundparticle3d.

Referenced by Finder::Process().

00028 {if(foundparticle)return foundparticle3d;return 0;};

void PrimaryShowerFinder::GetPeakAreaAverage ( double &  x,
double &  y,
double &  val,
int &  cnt,
int  curx,
int  cury,
TH2D *  hist,
std::vector< double > &  peakvalue,
std::vector< int > &  peakstillgood,
std::vector< int > &  peakbinx,
std::vector< int > &  peakbiny 
) [private]

Definition at line 3211 of file PrimaryShowerFinder.cxx.

Referenced by MakeHoughMap().

03212 {
03213 
03214         double thisval = hist->GetBinContent(curx,cury);
03215 
03216         if(thisval==0)return;
03217         
03218         val+=thisval;
03219         x+=hist->GetXaxis()->GetBinCenter(curx);
03220         y+=hist->GetYaxis()->GetBinCenter(cury);
03221         cnt++;
03222         hist->SetBinContent(curx,cury,0);
03223 
03225         
03226         for(int i=0;i<(int)peakbinx.size();i++)
03227         {
03228                 if(peakbinx[i]==curx && peakbiny[i]==cury)
03229                 {
03230                         peakstillgood[i]=0;
03231                         peakvalue[i]=0;
03232                         break;
03233                 }
03234         }
03235         
03236         
03238 
03239 
03240         
03241         GetPeakAreaAverage(x, y, val,cnt, curx+1, cury, hist, peakvalue, peakstillgood, peakbinx, peakbiny);
03242         GetPeakAreaAverage(x, y, val,cnt, curx, cury+1, hist, peakvalue, peakstillgood, peakbinx, peakbiny);
03243         GetPeakAreaAverage(x, y, val,cnt, curx, cury-1, hist, peakvalue, peakstillgood, peakbinx, peakbiny);
03244         GetPeakAreaAverage(x, y, val,cnt, curx-1, cury, hist, peakvalue, peakstillgood, peakbinx, peakbiny);
03245         GetPeakAreaAverage(x, y, val,cnt, curx+1, cury+1, hist, peakvalue, peakstillgood, peakbinx, peakbiny);
03246         GetPeakAreaAverage(x, y, val,cnt, curx-1, cury-1, hist, peakvalue, peakstillgood, peakbinx, peakbiny);
03247         GetPeakAreaAverage(x, y, val,cnt, curx+1, cury-1, hist, peakvalue, peakstillgood, peakbinx, peakbiny);
03248         GetPeakAreaAverage(x, y, val,cnt, curx-1, cury+1, hist, peakvalue, peakstillgood, peakbinx, peakbiny);
03249 
03250                 
03251 }

void PrimaryShowerFinder::LoadCloseHits ( HoughLine hl,
Managed::ClusterManager cm,
int  view,
double  max_tdist,
int  limit_to_current_size = 0 
) [static]

Definition at line 2656 of file PrimaryShowerFinder.cxx.

References HoughLine::AddCluster(), Managed::ManagedCluster::e, HoughLine::end_t, HoughLine::end_z, Managed::ClusterManager::GetCluster(), Managed::ClusterManager::GetClusterMap(), HoughLine::GetExpectedT(), HoughLine::offset_t, HoughLine::offset_z, HoughLine::phi, HoughLine::r, HoughLine::start_t, HoughLine::start_z, Managed::ManagedCluster::t, HoughLine::theta, and Managed::ManagedCluster::z.

Referenced by Bundle(), MakeChains(), and MakeHoughMap().

02657 {
02658         if(!hl || !cm)return;
02659 
02660         std::map<double, std::map<double, int> > * cluster_map = cm->GetClusterMap(view);
02661     std::map<double, std::map<double, int> >::reverse_iterator p_iterr;
02662     std::map<double, int >::iterator s_iterr;
02663 
02664 
02665 //printf("\nloading hit\n");
02666         double last_z=0;
02667         Managed::ManagedCluster *closest_cluster=0;
02668         double dist=1000;
02669     for(p_iterr=cluster_map->rbegin();p_iterr!=cluster_map->rend(); p_iterr++)
02670     {   
02671         if(!last_z)last_z=p_iterr->first;
02672                 
02673   //    printf("plane %f \n",p_iterr->first);
02674         for(s_iterr=p_iterr->second.begin();s_iterr!=p_iterr->second.end(); s_iterr++)
02675                 {
02676                                 Managed::ManagedCluster *cl = cm->GetCluster(s_iterr->second);
02677                                 if(!cl)continue;
02678                                 if(cl->e<0.001)continue;
02679                                 
02680                                 //do we need to make sure that this cluster is contained withing the dimensions of the houghline?
02681                                 if(limit_to_current_size)
02682                                 {
02683                                         double min_z = hl->start_z < hl->end_z?hl->start_z:hl->end_z;
02684                                         double max_z = hl->start_z > hl->end_z?hl->start_z:hl->end_z;
02685                                         if(min_z-cl->z>0.07 || cl->z-max_z>0.07)
02686                                         {
02687                         //                      printf("not adding %f %f to hl with endpts %f %f  %f %f\n",cl->z,cl->t,hl->start_z,hl->start_t,hl->end_z,hl->end_t);
02688                                                 continue;
02689                                         }
02690                                         
02691                                         //if this hl is contrained to a single plane... we need to keep hits within t!
02692                                         double min_t =  hl->start_t < hl->end_t ? hl->start_t:hl->end_t;
02693                                         double max_t =  hl->start_t > hl->end_t ? hl->start_t:hl->end_t;
02694                                         
02695                                         if(max_z-min_z<0.001 && (cl->t-min_t < 0.001 || cl->t - max_t > 0.001))continue;
02696                                         
02697                                         if(fabs(hl->phi-3.141592/2 )<0.1)
02698                                         {
02699 
02700                                                 
02701                                                 if(cl->t<min_t || cl->t>max_t)continue;
02702                                         }
02703                                 
02704                                 }                               
02705                                 
02706                                 double exp_t = hl->GetExpectedT(cl->z);
02707 //                              printf(" cluster z %f\n",cl->z);
02708                                 //is it a vertical line?
02709                                 if(fabs(hl->phi-3.141592/2 )<0.1){
02710                                 
02711                                         //find where t=0
02712                                         double zpos = hl->offset_z + (hl->offset_t + hl->r/sin(hl->theta))*sin(hl->theta)/cos(hl->theta);
02713                         //      printf("verticle match at %f against %f\n",zpos,cl->z); 
02714                                         //give more for the verticle hits       
02715                                         double tmt = max_tdist > 0.1 ? max_tdist:0.1;
02716                                         if(fabs(zpos-cl->z)<tmt)
02717                                         {
02718                                                 hl->AddCluster(cl);
02719                                         }
02720                                 
02721                                 }else{
02722                                         //is it the closest?
02723                                         if(fabs(exp_t-cl->t)<max_tdist && fabs(exp_t-cl->t)<dist)
02724                                         {
02725                                                 closest_cluster=cl;
02726                                                 dist=fabs(exp_t-cl->t);
02727                                         }
02728                                 }       
02729                                         
02730                 }
02731 
02732                 std::map<double, std::map<double, int> >::reverse_iterator next_p_iterr=p_iterr;
02733                 next_p_iterr++;
02734         if(next_p_iterr!=cluster_map->rend())
02735         if(fabs(p_iterr->first-next_p_iterr->first)>0.04 && closest_cluster)
02736         {//new plane.. save the old plane
02737                         hl->AddCluster(closest_cluster);
02738 //                      printf("next plane %f from %f\n",p_iterr->first,last_z);
02739   //            printf("adding %f %f %f\n",closest_cluster->z,closest_cluster->t,closest_cluster->e);
02740                 closest_cluster=0;
02741                 dist=1000;
02742                 
02743                 last_z=p_iterr->first;
02744         }
02745 
02746 
02747         }
02748 
02749         if(closest_cluster)
02750         {//new plane.. save the old plane
02751  //     printf("last_plane %f\n",last_z);
02752  //                     printf("adding %f %f %f\n",closest_cluster->z,closest_cluster->t,closest_cluster->e);
02753                 hl->AddCluster(closest_cluster);
02754         }
02755 }

void PrimaryShowerFinder::MakeChains ( int  view  )  [private]

printf("new from split\n");

Definition at line 95 of file PrimaryShowerFinder.cxx.

References ChainHelper::AddFinishedChain(), ChainHelper::AttachAt(), chu, chv, CleanHoughLines(), Managed::ClusterManager::ClearInUse(), HoughLine::cluster_id, cluster_manager, Munits::cm, Managed::ManagedCluster::e, HoughLine::end_t, Chain::end_t, HoughLine::end_z, Chain::end_z, Chain::entries, ChainHelper::finished, foundvertex, foundvertexU, foundvertexV, ChainHelper::GetChain(), Managed::ClusterManager::GetCluster(), Managed::ClusterManager::GetClusterMap(), HoughLine::GetExpectedT(), Managed::ManagedCluster::GetStatus(), houghlinesU, houghlinesV, Managed::ManagedCluster::id, Chain::insert(), Msg::kDebug, LoadCloseHits(), Managed::ClusterManager::MarkUsed(), MSG, HoughLine::ncluster, ChainHelper::NewChain(), HoughLine::phi, Managed::ClusterManager::SetClusterSaver(), Managed::ClusterManager::SetHitManager(), SplitHoughLine(), Chain::start_t, HoughLine::start_t, HoughLine::start_z, Chain::start_z, HoughLine::sum_e, Managed::ManagedCluster::t, Managed::ManagedCluster::view, vtx_u, vtx_v, vtx_z, and Managed::ManagedCluster::z.

Referenced by FindPrimaryShower().

00096 {
00097 
00098         ChainHelper *ch = view==2?chu:chv;
00099         if(!ch)return;
00100         
00101 
00102         
00103 
00104         
00105         std::vector<HoughLine> * hl = view==2?&houghlinesU:&houghlinesV;
00106         
00107         if(hl->size()<1)
00108         {
00109                 //does the view have only one cluster? if so... make that a chain!
00110                 std::map<double, std::map<double, int> > * cluster_map = cluster_manager->GetClusterMap(view);
00111                 std::map<double, std::map<double, int> >::reverse_iterator p_iterr;
00112         std::map<double, int >::iterator s_iterr;
00113 
00114                 Managed::ManagedCluster *cluster=0;
00115                 int cluster_count=0;
00116                 for(p_iterr=cluster_map->rbegin();p_iterr!=cluster_map->rend(); p_iterr++)
00117         {       
00118         
00119                 for(s_iterr=p_iterr->second.begin();s_iterr!=p_iterr->second.end(); s_iterr++)
00120                         {
00121                                 cluster = cluster_manager->GetCluster(s_iterr->second);
00122                                 cluster_count++;
00123                         }
00124         }
00125         if(cluster_count==1 && cluster)
00126         {
00127                         int nc = ch->NewChain();
00128                         Chain *c = ch->GetChain(nc);
00129                         
00130                 
00131                 //      int newid = cluster_manager->SaveCluster(cluster->id,cluster->e,2);
00132                 //      cluster = cluster_manager->GetCluster(newid);
00133                         
00134                         if(cluster)
00135                         {
00136                                 c->insert(cluster->t, cluster->z, cluster->e, cluster->id);
00137                         }       
00138                         MSG("PrimaryShowerFinder",Msg::kDebug)<<"single cluster chain made in view "<<view<<"\n";
00139         
00140         }
00141                 
00142         
00143                 return;
00144         }
00145         
00146         
00147         //make a copy....
00148         
00149         std::vector<HoughLine> hlcopy = *hl;
00150 
00151         hl=&hlcopy;
00152 
00153         Managed::ClusterManager *mycm=cluster_manager;
00154         
00155 
00156         //probably vertex 
00157         //to check if chains after the first are vertex pointing or other chain pointing
00158         double prob_vtx_z=0;
00159         double prob_vtx_t=0;
00160 
00161         if(view ==2 && foundvertexU || foundvertex)
00162         {
00163                 prob_vtx_z=vtx_z;
00164                 prob_vtx_t=vtx_u;
00165                 
00166                 //printf("USING VERTEX U %f %f\n",prob_vtx_z,prob_vtx_t);
00167         }
00168 
00169         if(view ==3 && foundvertexV || foundvertex)
00170         {
00171                 prob_vtx_z=vtx_z;
00172                 prob_vtx_t=vtx_u;
00173                 
00174                 //printf("USING VERTEX V %f %f\n",prob_vtx_z,prob_vtx_t);
00175         }       
00176 
00178         //uncomment to do a temp        
00179         Managed::ClusterManager cm=*cluster_manager;
00180         cm.SetClusterSaver(0);
00181         cm.SetHitManager(0);
00182         
00183         mycm=&cm;
00184         mycm->ClearInUse();
00185 
00187 
00188         
00189         std::vector<int>foundChains;//store the found chains here as well... so we can split up hough lines as needed
00190         
00191         
00192         //need to load any existing chains from the long muon finder!!!!
00194         for(unsigned int i=0;i<ch->finished.size();i++)
00195         {
00196                 foundChains.push_back(ch->finished[i].myId);
00197         
00198         }
00199         
00200         
00201         
00203         
00204 
00205         //do a full reload...
00206         for(unsigned int k=0;k<hl->size();k++)
00207         {
00208 //              printf("%d\n",k);
00209                 (*hl)[k].ResetHits(0);
00210                 LoadCloseHits(&(*hl)[k],mycm,view,0.0412*1.5);
00211 //              printf("%f %f   %f %f\n",(*hl)[k].start_z,(*hl)[k].start_t,(*hl)[k].end_z,(*hl)[k].end_t);      
00212         }
00213         std::vector<HoughLine> hlcopy2 = *hl;
00214                 
00215 //      
00216 //      printf("!!!!!!!!!!!!!!!!!!!!!!!\nlooking for chains\n");
00217         for(int ccc=0;ccc<10;ccc++)
00218         {
00219         //      printf("try %d with %d houghlines\n",ccc,hl->size());
00220                 for(unsigned int k=0;k<hl->size();k++)
00221                 {
00222         
00223         //              printf("%d\n",k);
00224                         (*hl)[k].ResetHits(1); //keep the bounds
00225                         
00226                 //      printf("reset hl %d\n",k);
00227                         LoadCloseHits(&(*hl)[k],mycm,view,0.0412*1.5,1);
00228         //              printf("current %d,%f %f %f %f\n",(&(*hl)[k])->ncluster,(&(*hl)[k])->start_z,(&(*hl)[k])->start_t,(&(*hl)[k])->end_z,(&(*hl)[k])->end_t);
00229                         
00230                         
00231 //                      printf("loaded %d hits\n",(*hl)[k].ncluster);
00232                         
00233                 
00234                         //does this hough line intersect any found chains?  if so, we should split it!
00235                         for(unsigned int c=0;c<foundChains.size();c++)
00236                         {
00237 //                              printf("c %d of %d\n",c,foundChains.size());
00238 //                              printf("current %d,%f %f %f %f\n",(&(*hl)[k])->ncluster,(&(*hl)[k])->start_z,(&(*hl)[k])->start_t,(&(*hl)[k])->end_z,(&(*hl)[k])->end_t);
00239                                 
00240                                 Chain *chain = ch->GetChain(foundChains[c]);
00241                                 if(!chain)continue;
00242                                 if(chain->entries<2)continue;
00243                                 
00244                                 HoughLine *lnew = SplitHoughLine(chain,&(*hl)[k],mycm);
00245 //                              printf("split  %d,%f %f %f %f\n",(&(*hl)[k])->ncluster,(&(*hl)[k])->start_z,(&(*hl)[k])->start_t,(&(*hl)[k])->end_z,(&(*hl)[k])->end_t);
00246                                 
00247                                                         
00248                                 if(lnew && lnew->ncluster>0)
00249                                 {
00251                                         hl->push_back(*lnew);
00252 //                                      printf("%d %f %f %f %f\n",lnew->ncluster,lnew->start_z,lnew->start_t,lnew->end_z,lnew->end_t);
00253                                 }
00254                                 
00255                                 if(lnew)delete lnew;
00256                                 //delete it if it is empty!
00257                                 if((*hl)[k].ncluster<1)
00258                                 {
00259 //                                      printf("hl size %d hl at position %d is empty...\n",hl->size(),k);
00260                                         
00261                                 //      std::vector<HoughLine> hh;
00262                                 //      for(unsigned int aw=0;aw<hl->size();aw++)
00263                                 //              if(aw!=k)hh.push_back((*hl)[aw]);
00264                                         
00265                                 //      *hl=hh;
00266                                         
00267                                         hl->erase(hl->begin()+k,hl->begin()+k+1);
00268                                         k--;
00269 //                                      printf("%d %d\n",k,hl->size());
00270                                         break;//restart the check now that this hl is split
00271                                 }
00272 
00273                         }
00274 //                      printf("e %d\n",k);
00275                 }
00276         
00277                 CleanHoughLines(0,hl);
00278                 
00279                 //so we know where to attach secondary chains
00280                 std::vector<int>closest_chain;
00281                 std::vector<double>closest_chain_z;
00282                 for(unsigned int k=0;k<hl->size();k++)
00283                 {
00284                         closest_chain.push_back(-1);
00285                         closest_chain_z.push_back(0);
00286                 }
00287 
00288 
00289                 std::vector<int>toskip;
00290 
00291                 if( ( view ==2 && foundvertexU == 1 ) || ( view==3 && foundvertexV==1)  || foundvertex) //if we don't have a vertex, we don't have any chains... so consider all of them and skip this part....
00292                 for(unsigned int k=0;k<hl->size();k++)
00293                 {
00294                         HoughLine * l = &(*hl)[k];
00295 
00296                         int keep=0;
00297                         //check for chain pointing      
00298                         double lastdist=10000;
00299                         for(unsigned int i=0;i<foundChains.size();i++)
00300                         {       
00301                                 Chain *c = ch->GetChain(foundChains[i]);
00302                                 if(!c)continue;
00303                                 if(c->entries<2)continue; //can come from when a chain is split at the first cluster
00304                 
00305 
00306 
00307                                 //measure the distance of intersection
00308                                 double vz = vtx_z;
00309                                 double vt = view==2? vtx_u:vtx_v;
00310                                 
00311                                 
00312                                 double c_start_z, c_start_t, c_end_z, c_end_t =0;
00313                                 
00314                                 if( fabs((c->start_z-vz)*(c->start_z-vz)+(c->start_t-vt)*(c->start_t-vt)) < fabs((c->end_z-vz)*(c->end_z-vz)+(c->end_t-vt)*(c->end_t-vt)) )
00315                                 {
00316                                         c_start_z = c->start_z;
00317                                         c_start_t = c->start_t;
00318                                         c_end_z = c->end_z;
00319                                         c_end_t = c->end_t;     
00320                                 }else{
00321                                         c_start_z = c->end_z;
00322                                         c_start_t = c->end_t;
00323                                         c_end_z = c->start_z;
00324                                         c_end_t = c->start_t;                           
00325                                 }
00326                                 
00327                                  if(!(c_end_z-c_start_z))continue;
00328 
00329                                 double l_start_z, l_start_t, l_end_z, l_end_t =0;
00330                                 
00331                                 if( fabs((l->start_z-c_start_z)*(l->start_z-c_start_z)+(l->start_t-c_start_t)*(l->start_t-c_start_t)) < fabs((l->end_z-c_start_z)*(l->end_z-c_start_z)+(l->end_t-c_start_t)*(l->end_t-c_start_t)) )
00332                                 {
00333                                         l_start_z = l->start_z;
00334                                         l_start_t = l->start_t;
00335                                         l_end_z = l->end_z;
00336                                         l_end_t = l->end_t;     
00337                                 }else{
00338                                         l_start_z = l->end_z;
00339                                         l_start_t = l->end_t;
00340                                         l_end_z = l->start_z;
00341                                         l_end_t = l->start_t;                           
00342                                 }
00343                                 
00344 
00345                                 //find intersection
00346                                 
00347                                 double lslope = (l_end_t-l_start_t) / (l_end_z-l_start_z) ;
00348                                 double loffset = l_end_t - l_end_z*lslope;
00349                                 
00350                                 double cslope = (c_end_t-c_start_t) / (c_end_z-c_start_z) ;
00351                                 double coffset = c_end_t - c_end_z*cslope;
00352                                 
00353                                 double int_z = cslope-lslope ? (loffset - coffset ) / ( cslope - lslope):l_start_z;
00354                                 double int_t = cslope * int_z + coffset;
00355                                                                                                 
00356                                 //require intersection to be behind the start of the chain
00357 
00358                                 //if(fabs((c_start_z-vz)*(c_start_z-vz)+(c_start_t-vt)*(c_start_t-vt)) > fabs((int_z-vz)*(int_z-vz)+(int_t-vt)*(int_t-vt)) )continue;
00359                                 
00360                                 //calculate distance along chain path... where start is 0 and going towards end increases
00361                                 //if its <0  its bad....
00362                                 double r_c_end_z=c_end_z-c_start_z;
00363                                 double r_c_end_t=c_end_t-c_start_t;
00364                                 double r_l_end_z=l_end_z-c_start_z;
00365                                 double r_l_end_t=l_end_t-c_start_t;
00366                                 double r_l_start_z=l_start_z-c_start_z;
00367                                 double r_l_start_t=l_start_t-c_start_t;
00368                                 
00369                                 double theta = atan(r_c_end_t/r_c_end_z);                               
00370                                 double rr_int_z = (int_z-c_start_z) * cos(theta) + (int_t-c_start_t) * sin(theta);                              
00371                                 double rr_int_t = (int_t-c_start_t) * cos(theta) - (int_z-c_start_z) * sin(theta);                                      
00372                                 if(rr_int_z<0.03)continue;// points to far fowards of chain
00373                                 
00374                                 
00375                                 
00376                                 
00377                                 double ltheta = atan( ( l->end_t-l->start_t) / (l->end_z-l->start_z) );                         
00378                                 double rr_l_int_z = (int_z-l->start_z) * cos(ltheta) + (int_t-l->start_t) * sin(ltheta);                                
00379         
00380                                 //printf("in hl coords, int is at z %f\n",rr_l_int_z);
00381         
00382                                 //require the intersection to be forwards of the hl!
00383                                 if(rr_l_int_z>0)continue;                               
00384                                 
00385                                                 
00386                                 
00387                                 //its the intersection past the end of the chain?  if so... check the relative angle between the hl and the chain to keep it reasonable
00388                                 double rr_c_end_z = r_c_end_z * cos(theta) + r_c_end_t * sin(theta);                                    
00389                                 //double rr_c_end_t =  -r_c_end_z * sin(theta) + r_c_end_t * cos(theta);                                
00390 
00391 
00392                                 //rotate the coordinates so the chain is along the x axis
00393 
00394                                 double rr_l_end_z = r_l_end_z * cos(theta) + r_l_end_t * sin(theta);
00395                                 double rr_l_end_t =  -r_l_end_z * sin(theta) + r_l_end_t * cos(theta); 
00396                                 double rr_l_start_z = r_l_start_z * cos(theta) + r_l_start_t * sin(theta);
00397                                 double rr_l_start_t =  -r_l_start_z * sin(theta) + r_l_start_t * cos(theta);                            
00398         
00399                                 
00400                                 if(rr_int_z > rr_c_end_z)
00401                                 {
00402                                         //count the angle from the intersection!
00403                                         double dy = rr_l_start_t - rr_int_t;
00404                                         double dx = rr_l_start_z - rr_int_z;
00405                                         
00406                                         double intangle = dx ? atan(dy/dx) : 0;
00407                                         if(intangle>3.141592/6 || dx==0)continue; //to steep!
00408                                 
00409                                 }
00410                                 
00411                                 //is it the closest?  
00412                                 
00413                                 double dist = sqrt((l_start_z - int_z)*(l_start_z - int_z)+(l_start_t - int_t)*(l_start_t - int_t));
00414                                 
00415                                 //if the intersection is past the end of the chain... need to include the distance from the intersection to the end of the chain
00416                                 if(fabs((c_end_z-vz)*(c_end_z-vz)+(c_end_t-vt)*(c_end_t-vt)) < fabs((int_z-vz)*(int_z-vz)+(int_t-vt)*(int_t-vt)) )
00417                                 {
00418                                         dist+= sqrt((c_end_z - int_z)*(c_end_z - int_z)+(c_end_t - int_t)*(c_end_t - int_t));
00419                                 }
00420 
00421                                 //is the closest ?
00422                                 if(dist > lastdist)continue;
00423                                 //printf("closer to chain %f %f %f %f at z %f\n",c_start_z,c_start_t,c_end_z,c_end_t,int_z);
00424                                 //if so, check for orientation before recording
00425                         
00426                                 if(fabs(rr_l_end_t) <= fabs(rr_l_start_t) && rr_l_start_z < rr_l_end_z ||
00427                                         fabs(rr_l_end_t) >= fabs(rr_l_start_t) && rr_l_start_z > rr_l_end_z 
00428                                  )
00429                                 {
00430                                         //its closer... but no good... we want to invalidate the last guess...
00431                                         closest_chain[k]=-1;
00432                                         closest_chain_z[k]=int_z;
00433                                         lastdist=dist;  
00434                                         keep=0;                         
00435                                         continue;
00436                                 }
00437 
00438                                 //printf("in chain coords line has start %f %f end %f %f\n",rr_l_start_z,rr_l_start_t,rr_l_end_z,rr_l_end_t);
00439                                 //printf("hl start %f %f end %f %f\n",l_start_z,l_start_t,l_end_z,l_end_t);
00440                                 //record it
00441                                 closest_chain[k]=foundChains[i];
00442                                 closest_chain_z[k]=int_z;
00443                                 lastdist=dist;
00444                                 keep=1;
00445                                 
00446                                 //printf("keeping %d it points to chain %f %f %f %f at z %f\n",k,c_start_z,c_start_t,c_end_z,c_end_t,int_z);
00447                                 
00448                         }
00449                         
00450                         
00451                         //dont have anything yet?
00452                         //see if we are vertex pointing
00453                         if(closest_chain[k]<0)
00454                         {
00455                                 double exp_vtx_t = l->GetExpectedT(vtx_z);
00456                                 if((view==2&&foundvertexU || view==3&&foundvertexV) && fabs(exp_vtx_t-(view==2?vtx_u:vtx_v))<2*0.0451) //smaller than 2 strips away?
00457                                 {
00458                                 
00459                                         //need to make sure that the vertex is not within the hl!
00460         
00461                                 
00462                                         double theta = atan( ( l->end_t-l->start_t) / (l->end_z-l->start_z) );                          
00463                                         double rr_end_z = (l->end_z-l->start_z) * cos(theta) + (l->end_t-l->start_t) * sin(theta);                              
00464                                         //double rr_end_t = (l->end_t-l->start_t) * cos(theta) - (l->end_z-l->start_z) * sin(theta);                                    
00465 
00466                                         double rr_vtx_z = (vtx_z-l->start_z) * cos(theta) + (exp_vtx_t-l->start_t) * sin(theta);                                
00467                                         //double rr_vtx_t = (exp_vtx_t-l->start_t) * cos(theta) - (vtx_z-l->start_z) * sin(theta);      
00468         
00469                                         //printf("in hl coords:  hl 0 0 to %f %f  vtx %f %f\n",rr_end_z,rr_end_t,rr_vtx_z,rr_vtx_t);
00470         
00471                                         if(rr_vtx_z< 0 || rr_vtx_z > rr_end_z)
00472                                         {                               
00473                                                 //printf("keeping %d vp estimate vtx %f %f actual %f %f\n",k,vtx_z,l->GetExpectedT(vtx_z),vtx_z,view==2?vtx_u:vtx_v);
00474                                                 keep=1; //we'll keep it
00475                                         }
00476                                 }       
00477                         }
00478                 
00479 
00480                         if(keep)
00481                         {
00482                                 //printf("keeping %d\n",k);
00483                                 continue;
00484                         }
00485                         //if we get here... we don't want this hl!
00486                         //printf("skipping %d\n",k);
00487                         toskip.push_back(k);
00488                 
00489                 }
00490 
00491 
00492 /*
00493 
00494 
00495                 std::vector<int>toskip;
00496                 //make sure that all hough line are vertex pointing or pointing to a another chain and on the correct angle
00497                 if(foundChains.size()>0)
00498                 {
00499                         for(unsigned int k=0;k<hl->size();k++)
00500                         {
00501                                 HoughLine * l = &(*hl)[k];
00502         //                                      printf("current %d,%f %f %f %f\n",(&(*hl)[k])->ncluster,(&(*hl)[k])->start_z,(&(*hl)[k])->start_t,(&(*hl)[k])->end_z,(&(*hl)[k])->end_t);
00503                                                 
00504                                                         
00505                                 int keep=0;                             
00506                                 //is it vertex pointing?
00507                                 if(foundvertex && fabs(l->GetExpectedT(prob_vtx_z)-prob_vtx_t)<2*0.0451) //smaller than 2 strips away?
00508                                 {
00509                                         printf("keeping %d vp\n",k);
00510                                         keep=1; //we'll keep it
00511                                         //we still want to see if it is pointing to another chain so we know if it needs to be a secondary
00512                                 }
00513 
00514                                 //now we need to see if it is pointing to a current chain ... and find the closest one!
00515                                 int closest_pointing=-1;
00516                                 double closest_dist=0;
00517                                 double dist=0;          
00518                                 for(unsigned int i=0;i<foundChains.size();i++)
00519                                 {
00520                                         Chain *c = ch->GetChain(foundChains[i]);
00521                                         double exp_t_start = l->GetExpectedT(c->start_z);
00522                                         double exp_t_end = l->GetExpectedT(c->end_z);
00523 
00524                         
00525                                         if(( exp_t_start < c->start_t && exp_t_end > c->end_t ) || ( exp_t_start > c->start_t && exp_t_end < c->end_t))
00526                                         {
00527                                         //              //check for orientation
00528                                         //              
00529                                         //              if(fabs(c->end_z - l->end_z) > fabs(c->end_z - l->start_z))continue;
00530                                         //              if(fabs(c->start_z - l->start_z) > fabs(c->start_z - l->end_z))continue;
00531                                                         
00532                                                         MSG("PrimaryShowerFinder",Msg::kDebug)<<"Passing orientation on "<<c->start_z<<" "<<c->start_t<<" "<<c->end_z<<" "<<c->end_t<<"\n";
00533 
00534                                                         //check the distance
00535                                                         //the way we found the intersections requires the intersection to be within the chain....
00536                                                         for(unsigned int kk=0;kk<c->z.size()-1;kk++)
00537                                                         {
00538                                                                 double exp_t = l->GetExpectedT(c->z[kk]);
00539                                                                 double exp_t_n = l->GetExpectedT(c->z[kk+1]);
00540                                                                 
00541                                                                 if(( exp_t>=c->t[kk]&&exp_t_n<=c->t[kk+1] ) || ( exp_t<=c->t[kk]&&exp_t_n>=c->t[kk+1] ) )
00542                                                                 {//intersection is bewteen these two points!
00543                                                         
00544                                                                                 //for now, put it close to the more forward hit...
00545                                                                                 //later we can actually extrapolate it
00546                                                                                 
00547                                                                                 double t = c->t[kk];
00548                                                                                 dist = fabs(t-exp_t);
00549                                                                                 closest_chain_z[k]=c->z[kk];
00550                                                                         
00551                                                                 }
00552                                                         
00553                                                         }
00554                                                         
00555                                                         
00556                                                 if(closest_pointing<0)//the first one
00557                                                 {
00558                                                         closest_pointing=foundChains[i];
00559                                                         closest_dist =  dist;           
00560                                                 }else
00561                                                 {
00562                                                         if(dist<closest_dist)
00563                                                         {
00564                                                                 closest_pointing=foundChains[i];
00565                                                                 closest_dist =  dist;                                                           
00566                                                         }
00567                                                 }
00568                                 
00569                                                 
00570 
00571 
00572                                         }
00573                                 }
00574                                 
00575                                 for(unsigned int i=0;i<foundChains.size();i++)
00576                                 {               
00577                                                                         Chain *c = ch->GetChain(foundChains[i]);        
00578                                 
00579         //printf("end chain %f %f  start hl %f %f\n",c->end_z,c->end_t,l->start_z,l->start_t);
00580 
00581                                 
00582                                 
00583                                 if(c->end_z < l->start_z)
00584                                 {
00585                                 
00586                                                         //also check to see if this chain is at the end of another chain
00587                                                         if(
00588                                                                 fabs(c->end_t - l->GetExpectedT(c->end_z))<0.1  
00589                                                                 && l->start_z - c->end_z  <0.07*3)
00590                                                         {
00591                                                         
00592                                                                 closest_chain_z[k]=c->end_z;
00593                                                                 dist = fabs(c->end_t - l->GetExpectedT(c->end_z));
00594 
00595                                                         }
00596 
00597                                                 //      printf("end chain %f %f  start hl %f %f\n",c->end_z,c->end_t,l->start_z,l->start_t);
00598                                 
00599                                 }
00600                                 
00601                                 
00602                                 
00603                                                 if(closest_pointing<0)//the first one
00604                                                 {
00605                                                         closest_pointing=foundChains[i];
00606                                                         closest_dist =  dist;           
00607                                                 }else
00608                                                 {
00609                                                         if(dist<closest_dist)
00610                                                         {
00611                                                                 closest_pointing=foundChains[i];
00612                                                                 closest_dist =  dist;                                                           
00613                                                         }
00614                                                 }
00615                                 
00616                                 
00617                                 
00618                                 }
00619                                 
00620                                 
00621                                 if(closest_pointing>-1)
00622                                 {
00623                                                 Chain *c = ch->GetChain(closest_pointing);
00624                                                 //its pointing... so now check the orientation
00625                                                 double close_z;
00626                                                 double close_t;
00627                                                 printf("prob vtx z %f t %f\n",prob_vtx_z,prob_vtx_t);
00628                                                 printf("chain start z %f t %f, end z %f t %f\n",c->start_z,c->start_t,c->end_z,c->end_t);
00629                                                 if(fabs((c->start_t-prob_vtx_t)*(c->start_t-prob_vtx_t)+(c->start_z-prob_vtx_z)*(c->start_z-prob_vtx_z))
00630                                                         <
00631                                                         fabs((c->end_t-prob_vtx_t)*(c->end_t-prob_vtx_t)+(c->end_z-prob_vtx_z)*(c->end_z-prob_vtx_z)))
00632                                                 {       
00633                                                         close_z=c->start_z;
00634                                                         close_t=c->start_t;
00635                                                 }
00636                                                 else
00637                                                 {
00638                                                         close_z=c->end_z;
00639                                                         close_t=c->end_t;
00640                                                 }
00641                                                 
00642                                                 printf("close z %f t %f\n",close_z,close_t);
00643                                                 double ch_close_z;
00644                                                 double ch_close_t;
00645                                                 double ch_far_z;
00646                                                 double ch_far_t;
00647                                                 
00648                         //                      if(fabs(l->start_z - close_z) < fabs(l->end_z - close_z))
00649                         //                      {
00650                                                 if(fabs((l->start_t-prob_vtx_t)*(l->start_t-prob_vtx_t)+(l->start_z-prob_vtx_z)*(l->start_z-prob_vtx_z))
00651                                                         <
00652                                                         fabs((l->end_t-prob_vtx_t)*(l->end_t-prob_vtx_t)+(l->end_z-prob_vtx_z)*(l->end_z-prob_vtx_z)))
00653                                                 {       
00654                         
00655                         
00656                                                         ch_close_z = l->start_z;
00657                                                         ch_close_t = l->start_t;
00658                                                         ch_far_z = l->end_z;
00659                                                         ch_far_t = l->end_t;
00660                                                 }else{
00661                                                         ch_close_z = l->end_z;
00662                                                         ch_close_t = l->end_t;
00663                                                         ch_far_z = l->start_z;
00664                                                         ch_far_t = l->start_t;
00665                                                 }
00666                                                 
00667                                                 
00668                                                 //this will cause issues with vertical chains!
00669                                                 double dst = fabs(c->interpolate(ch_close_z)-ch_close_t);
00670                                                 double det = fabs(c->interpolate(ch_far_z)-ch_far_t);
00671                                                 printf("!!! %f\n",c->interpolate(ch_far_z));
00672                                                 printf(" close z %f t %f  far z %f t %f  dst %f det %f\n",ch_close_z,ch_close_t,ch_far_z,ch_far_t,dst,det);
00673                                                 
00674                                                 if(dst<det)
00675                                                 {
00676                                                         keep=1;//we'll keep it
00677                                                         printf("keeping %d cp\n",k);
00678                                                         closest_chain[k]=closest_pointing;
00679                                                 }
00680                                                         
00681                                 
00682                                 }
00683                                 
00684                                 if(keep)continue;
00685                                 //if we get here... we don't want this hl!
00686                                 printf("skipping %d\n",k);
00687                                 toskip.push_back(k);
00688         
00689                         }
00690                 }
00691         
00692 
00693 
00694 */
00695 
00696 
00697 
00698 
00699 
00700 
00701 
00702         for(unsigned int i=0;i<hl->size();i++)
00703         {
00704                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"hl "<<i<<" "<<(*hl)[i].start_z<<" "<<(*hl)[i].start_t<<" "<<(*hl)[i].end_z<<" "<<(*hl)[i].end_t<<"\n";
00705         }
00706         
00707         
00708         
00709         std::map<double,int> orderer;
00710 //      printf("\n");
00711         for(int i=0;i<(int)hl->size();i++)
00712         {
00713 
00714                 int save=1;
00715                 for(int j=0;j<(int)toskip.size();j++)
00716                 {
00717                         if(toskip[j]==i)
00718                         {
00719                                 save =0;
00720                                 break;
00721                         }
00722                 }
00723                 if(!save)continue;
00724 
00725                 HoughLine * l = &(*hl)[i];
00726         
00727                 double weight = l->ncluster *cos(l->phi)*l->sum_e*l->sum_e;
00728 
00729                 orderer.insert(make_pair(weight,i));
00730                 
00731 //              printf("%d %f %f | %f\n",l->ncluster,l->phi,l->sum_e,weight);
00732         }               
00733 //      printf("\n");
00734 
00735         std::map<double,int>::reverse_iterator it_orderer;
00736 
00737         HoughLine *h = 0;       
00738 
00739         //int cnt=0;
00740         //for(it_orderer = orderer.rbegin();it_orderer!=orderer.rend();it_orderer++)
00741         //{
00742         
00743         it_orderer = orderer.rbegin();
00744         while(it_orderer!=orderer.rend())
00745         {
00746                 if((*hl)[it_orderer->second].ncluster<2)it_orderer++;
00747                 else break;
00748         }
00749         
00750         if(it_orderer==orderer.rend())break;
00751         
00752                 int i= it_orderer->second;
00753                 
00754                 HoughLine * l = &(*hl)[i];
00755                 h=l;
00756                 
00757                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"chose "<<l->ncluster<<" "<<l->phi<<" "<<l->sum_e<<" | "<<it_orderer->first<<"\n";
00758 
00759         //      cnt++;
00760         //      if(cnt>3)break;
00761 
00762 
00763                 int nc = ch->NewChain();
00764                 Chain *c = ch->GetChain(nc);
00765                 foundChains.push_back(nc);
00766 
00767                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"making chain from houghline "<<i<<" "<<l->start_z<<" "<<l->start_t<<" "<<l->end_z<<" "<<l->end_t<<"\n";
00768                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"STARTING WITH "<< ch->finished.size()<<" CHAINS\n";
00769         
00770                 for(unsigned int i=0;i<h->cluster_id.size();i++)
00771                 {
00772                         Managed::ManagedCluster *cl = mycm->GetCluster(h->cluster_id[i]);
00773                         
00774                         
00775                 //      int newid = mycm->SaveCluster(cl->id,cl->e,2);
00776                 //      cl = mycm->GetCluster(newid);
00777                         
00778                         if(cl)
00779                         {
00780                         
00781                                 
00782                                 
00783                                 
00784                                 c->insert(cl->t, cl->z, cl->e, cl->id);
00785                                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"("<<cl->t<<" "<< cl->z<<" "<< cl->e<<" "<< cl->id<<") s"<<cl->GetStatus()<<" v"<<cl->view<<"\n";
00786                         }
00787 
00788                         mycm->MarkUsed(cl->id);
00789                 
00790                 }
00791                 //attach the chain if needed...
00792                 
00793                 if(closest_chain[i]>-1)
00794                 {
00795 
00796                         Chain *mom = ch->GetChain(closest_chain[i]);
00797                         Chain newdaughter = ch->AttachAt(mom,c,closest_chain_z[i]);
00798                         if(newdaughter.entries>0)ch->AddFinishedChain(newdaughter);
00799                         MSG("PrimaryShowerFinder",Msg::kDebug)<<"kept by closest_chain, attempting to attach to chain "<<mom->start_z<<" "<<mom->start_t<<" "<<mom->end_z<<" "<<mom->end_t<<" at z "<<closest_chain_z[i]<<"\n";
00800                 }
00801                 
00802                 if(foundChains.size()==1 && !foundvertex)
00803                 {
00804                         prob_vtx_z=c->start_z;
00805                         prob_vtx_t=c->start_t;
00806                         
00807                         if(view==2)foundvertexU=1;
00808                         if(view==3)foundvertexV=1;
00809                         if(foundvertexU && foundvertexV) foundvertex=1;
00810                         vtx_z=prob_vtx_z;
00811                         if(view==2)vtx_u=prob_vtx_t;
00812                         if(view==3)vtx_v=prob_vtx_t;
00813                 }
00814         
00815                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"ENDING WITH "<<ch->finished.size()<<" CHAINS\n";
00816                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"done\n";
00817 
00818         //}
00819         
00820         
00821         
00822         }
00823 
00824         MSG("PrimaryShowerFinder",Msg::kDebug)<<"done making chains\n";
00825         
00826         mycm->ClearInUse();
00827 }

void PrimaryShowerFinder::MakeHoughMap ( int  view  )  [private]

Definition at line 2758 of file PrimaryShowerFinder.cxx.

References cluster_manager, CompareForwardAndClusters(), CompareLength(), copy(), Managed::ManagedCluster::e, HoughLine::end_t, HoughLine::end_z, Managed::ClusterManager::GetCluster(), Managed::ClusterManager::GetClusterMap(), GetPeakAreaAverage(), houghlinesU, houghlinesV, houghmapU, houghmapV, it, LoadCloseHits(), Munits::m, max, HoughLine::ncluster, Managed::ManagedCluster::Reset(), SaveHitGroup(), HoughLine::start_t, HoughLine::start_z, Managed::ManagedCluster::t, and Managed::ManagedCluster::z.

Referenced by FindPrimaryShower().

02759 {
02760         if(view !=2 && view !=3)return;
02761         std::map<double, std::map<double, int>  >::iterator p_iterr;
02762     std::map<double, int>::iterator s_iterr;
02763         
02764 //      printf("----view %d----\n\n",view);
02765         
02766         if(houghmapU && view==2){houghmapU->Reset();}//delete houghmapU;houghmapU=0;}
02767         if(houghmapV && view==3){houghmapV->Reset();}//delete houghmapV;houghmapV=0;}
02768         
02769 //      if(houghmapU==0)houghmapU= new TH2D("hhu","Hough U", 300, 0,3.141592,150,-1,1);
02770 //      if(houghmapV==0)houghmapV= new TH2D("hhv","Hough V", 300, 0,3.141592,150,-1,1);
02771 
02772         if(houghmapU==0)houghmapU= new TH2D("hhu","Hough U", 100, 0,3.141592,50,-1,1);
02773         if(houghmapV==0)houghmapV= new TH2D("hhv","Hough V", 100, 0,3.141592,50,-1,1);
02774 
02775 
02776                 
02777         TH2D * h = view==2?houghmapU:view==3?houghmapV:0;
02778         if(!h)return;
02779 
02780         //is the map empty?
02781         if(cluster_manager->GetClusterMap(view)->begin()==cluster_manager->GetClusterMap(view)->end())return; 
02782 
02783  //being on the first hit can caus some trouble... so move a bit away from it
02784         double off_z = cluster_manager->GetClusterMap(view)->begin()->first-0.1;
02785         double off_t = cluster_manager->GetClusterMap(view)->begin()->second.begin()->first-0.1;
02786         //Managed::ManagedCluster *largest=0;
02787         //double last_plane=0;
02788   
02789   
02790         TH2D copy;
02791         h->Copy(copy);    
02792 
02793         //find max value!
02794         double maxvale=0;
02795     for(p_iterr=cluster_manager->GetClusterMap(view)->begin();p_iterr!=cluster_manager->GetClusterMap(view)->end(); p_iterr++)
02796     {   
02797         for(s_iterr=p_iterr->second.begin();s_iterr!=p_iterr->second.end(); s_iterr++)
02798         {
02799                         Managed::ManagedCluster * clus = cluster_manager->GetCluster(s_iterr->second);
02800                         if(clus->e>maxvale)
02801                         {
02802                                 maxvale=clus->e;
02803                         }
02804                 }
02805         }
02806 
02807 
02808     for(p_iterr=cluster_manager->GetClusterMap(view)->begin();p_iterr!=cluster_manager->GetClusterMap(view)->end(); p_iterr++)
02809     {   
02810         for(s_iterr=p_iterr->second.begin();s_iterr!=p_iterr->second.end(); s_iterr++)
02811         {
02812                         Managed::ManagedCluster * clus = cluster_manager->GetCluster(s_iterr->second);
02813 //printf("adding with e %f\n",clus->e);
02814                         //don't care about low e hits...
02815                         if(0.2 > clus->e)continue;
02816                         for(double i=0;i<3.141592;i+=3.141592/(h->GetNbinsX()-1))
02817                         {
02818                                 //double val = clus->e;
02819                                 double val=1;
02820                                 h->Fill(i, (clus->z-off_z) * cos(i) + (clus->t-off_t) * sin(i),  val);  
02821                                 copy.Fill(i, (clus->z-off_z) * cos(i) + (clus->t-off_t) * sin(i), clus->e);
02822                         }                       
02823                 }
02824         }
02825         
02826 
02827    TH2D ch;
02828    h->Copy(ch);
02829         
02830 
02831         h->Reset();
02832         
02833         //int nx=ch.GetNbinsX();
02834         //int ny=ch.GetNbinsY();   
02835 
02836         while(ch.GetMaximum()>1)
02837         {
02838                 int x;
02839                 int y;
02840                 int z;
02841                 int m = ch.GetMaximumBin();
02842                 ch.GetBinXYZ(m,x,y,z);
02843 
02844                 SaveHitGroup(&ch,(TH2D*)h, (double)ch.GetBinContent(x,y),(double)ch.GetBinContent(x,y), x,y);
02845     }    
02846 
02847 
02848         multimap<double, std::pair<double,double> > top_map;
02849         int tops=0;     
02850 
02851 
02853 //make a list of all peaks....to speed up the TH1::GetMaximum(Bin) which is painfully slow....
02854 
02855         std::vector<int> peakbinx;
02856         std::vector<int> peakbiny;
02857         std::vector<double> peakx;
02858         std::vector<double> peaky;
02859         std::vector<double> peakvalue;
02860         std::vector<int> peakstillgood;
02861         
02862         for(int i=1;i<h->GetNbinsX()+1;i++)
02863         for(int j=1;j<h->GetNbinsY()+1;j++)
02864         {
02865                 double thisc = h->GetBinContent(i,j);
02866                 if(thisc<0.0001)continue;
02867                 
02868                 if(
02869                         thisc >= h->GetBinContent(i+1,j)
02870                         &&
02871                         thisc >= h->GetBinContent(i,j+1)
02872                         &&
02873                         thisc >= h->GetBinContent(i-1,j)
02874                         &&
02875                         thisc >= h->GetBinContent(i,j-1)
02876                         &&
02877                         thisc >= h->GetBinContent(i+1,j+1)
02878                         &&
02879                         thisc >= h->GetBinContent(i+1,j-1)
02880                         &&
02881                         thisc >= h->GetBinContent(i-1,j+1)
02882                         &&
02883                         thisc >= h->GetBinContent(i-1,j-1)
02884                 )
02885                 {
02886                         peakx.push_back(h->GetXaxis()->GetBinCenter(i));
02887                         peaky.push_back(h->GetYaxis()->GetBinCenter(j));
02888                         peakbinx.push_back(i);
02889                         peakbiny.push_back(j);
02890                         peakvalue.push_back(thisc);     
02891                         peakstillgood.push_back(1);
02892                 }
02893         
02894         }
02895 
02896 /*
02897         int npeaks = peakstillgood.size();
02898         
02899         printf("found %d peaks",npeaks);
02900         
02901         while(npeaks>0)
02902         {
02903         
02904                 double xv=0;
02905                 double yv=0;
02906                 double val=0;
02907                 int cnt =0;
02908                         
02909                 GetPeakAreaAverage(xv, yv,val, peakbinx,peakbiny,peakx,peaky,peakvalue,peakstillgood);
02910 
02911 
02912                 top_map.insert(make_pair((double)val, std::pair<double,double>(xv,yv)) );
02913                 tops++;
02914 
02915                 npeaks=0;
02916                 for(unsigned int i=0;i<peakstillgood.size();i++)        
02917                         npeaks+=peakstillgood[i]==1?1:0;
02918         }
02919         
02920         printf("%d peaks made\n",tops);
02921 */
02923 
02924 //there must be a better way....
02925 
02926         int npeaks=0;
02927         for(unsigned int i=0;i<peakstillgood.size();i++)        
02928                 npeaks+=peakstillgood[i]==1?1:0;
02929         
02930         
02931         
02932         while(npeaks>0)//h->GetMaximum()>0)
02933         {
02934                 double xv=0;
02935                 double yv=0;
02936                 double val=0;
02937                 int cnt =0;
02938                 
02939                 int x;
02940                 int y;
02941                 //int z;
02942 //              int m = h->GetMaximumBin();
02943 //              h->GetBinXYZ(m,x,y,z);
02944 
02945                 double max=0;
02946                 int maxi=-1;
02947                 for(int i=0;i<(int)peakvalue.size();i++)
02948                 {
02949                         if(peakvalue[i]>max && peakstillgood[i])
02950                         {
02951                                 max = peakvalue[i];
02952                                 maxi=i;
02953                         }
02954                 }
02955 
02956                 if(maxi<0)break;
02957                 
02958                 x=peakbinx[maxi];
02959                 y=peakbiny[maxi];
02960                 
02961                 GetPeakAreaAverage(xv, yv,val, cnt, x, y, (TH2D*)h, peakvalue, peakstillgood, peakbinx, peakbiny);
02962                 xv/=(double)cnt;
02963                 yv/=(double)cnt;
02964                 
02965                 val=copy.GetBinContent(copy.FindBin(xv,yv));
02966                 
02967                 top_map.insert(make_pair((double)val, std::pair<double,double>(xv,yv)) );
02968                 tops++;
02969                 
02970                 npeaks=0;
02971                 for(unsigned int i=0;i<peakstillgood.size();i++)        
02972                         npeaks+=peakstillgood[i]==1?1:0;
02973         
02974         }
02975 
02976 
02977 
02978 
02979         
02980         
02981 //      printf("printing top %d\n",tops);
02982         
02983         multimap<double, std::pair<double,double> >::reverse_iterator it = top_map.rbegin();
02984 
02985         
02986         std::vector<HoughLine> *houghlines = view==2?&houghlinesU:&houghlinesV;
02987         
02988         for(int k=0;k<tops;k++)
02989         {
02990                 double theta=it->second.first;
02991                 double r = it->second.second;   
02992                 HoughLine hl(theta, r, off_t, off_z);
02993                 
02994         //      printf("peak at theta %f r %f off_t %f off_z %f\n",theta,r,off_t,off_z);
02995         
02996                 LoadCloseHits(&hl,cluster_manager,view,0.0412);
02997                         
02998                 if(hl.ncluster>1)       
02999                         houghlines->push_back(hl);
03000 
03001                 it++;
03002 
03003                 
03004         }
03005         
03006         
03007         sort(&(*houghlines)[0],&(*houghlines)[houghlines->size()],CompareLength);
03008         
03009         //should remove duplicates --- (caused by the same clusters being matched to different hough map peaks)
03010         std::vector<HoughLine> houghlinestemp;
03011 
03012         houghlinestemp = *houghlines;
03013         houghlines->clear();
03014 
03015         for(int i=0;i<(int)houghlinestemp.size();i++)
03016         {
03017                 int dup=0;
03018                 HoughLine * l = &houghlinestemp[i];
03019                 for(int j=i+1;j<(int)houghlinestemp.size();j++)
03020                 {
03021                         HoughLine *r = &houghlinestemp[j];
03022                         if(l && r && l->start_z == r->start_z && l->start_t == r->start_t && l ->end_z == r->end_z && l->end_t == r->end_t)
03023                         {
03024                                 dup=1;
03025                                 break;
03026                         }
03027 
03028                 }
03029 
03030                 if(!dup && houghlinestemp[i].ncluster>1)houghlines->push_back(houghlinestemp[i]);
03031 
03032         }
03033 
03034 
03035 /* //this code is highly flawed...
03036         for(int i=0;i<houghlines->size();i++)
03037         {
03038                 HoughLine * l = &(*houghlines)[i];
03039                 HoughLine * r=0;
03040                 if(i<houghlines->size()-1)r = &(*houghlines)[i+1];
03041                 if((*houghlines)[i].ncluster>1)houghlinestemp.push_back((*houghlines)[i]);              
03042                 if(l && r && l->start_z == r->start_z && l->start_t == r->start_t && l ->end_z == r->end_z && l->end_t == r->end_t)
03043                 {
03044 
03045                         while (i<houghlines->size())
03046                         {
03047                                 i++;
03048                                 HoughLine * l = &(*houghlines)[i];
03049                                 HoughLine * r = &(*houghlines)[i+1];
03050                                 if(! (l->start_z == r->start_z && l->start_t == r->start_t && l ->end_z == r->end_z && l->end_t == r->end_t)) break;
03051                         }
03052                 }
03053         }
03054 
03055 */
03056 
03057         
03058 
03059         
03060         sort(&(*houghlines)[0],&(*houghlines)[houghlines->size()],CompareForwardAndClusters);
03061         
03062         //int max =-1;
03063         //if(houghlines->size()>10)max = houghlines->size()-10;
03064         for(int i=0;i<(int)houghlines->size();i++)
03065         {
03066                 if((*houghlines)[i].ncluster<1)break;
03067 //              printf("%d houghline with e %f clusters %d (t,z) (%f,%f) (%f,%f) chi2/ndf %f theta %f r %f phiz %f \n", i,(*houghlines)[i].sum_e, (*houghlines)[i].ncluster, (*houghlines)[i].start_t, (*houghlines)[i].start_z, (*houghlines)[i].end_t, (*houghlines)[i].end_z, (*houghlines)[i].chi2/(*houghlines)[i].ncluster,(*houghlines)[i].theta,(*houghlines)[i].r,(*houghlines)[i].phi);
03068                 
03069                                                                         
03070         }
03071         
03072         
03073 
03074 }

void PrimaryShowerFinder::MakeParticle3D (  )  [private]

Definition at line 2377 of file PrimaryShowerFinder.cxx.

References MuELoss::a, Particle3D::add_to_back(), spoint::chain, chain_u, chain_v, spoint::chainhit, Chain::cluster_id, cluster_manager, Chain::e, MuELoss::e, spoint::e, Chain::entries, Particle3D::entries, Particle3D::finalize(), foundparticle3d, Managed::ClusterManager::GetCluster(), Msg::kDebug, MSG, Chain::myId, Chain::parentChain, Particle3D::particletype, Chain::particletype, spoint::rms_t, Managed::ManagedCluster::rms_t, spointgreaterlmf(), Chain::t, spoint::t, spoint::view, spoint::z, and Chain::z.

Referenced by FindPrimaryShower().

02378 {
02379 
02380         //printf("making psf particle cu %d cv %d\n",chain_u->myId,chain_v->myId);
02381 
02382         //verify that the chains in both views have sufficient overlap
02383         
02384         
02385         //make the 3d particle from these chains
02386         if(foundparticle3d)delete foundparticle3d;
02387         foundparticle3d= new Particle3D();
02388 
02389 
02390         
02391         std::vector<spoint> spoints;
02392         
02393         double start =0;
02394         double end  =1000;
02395         
02396         double endu=0;
02397         double endv=0;
02398                 
02399         
02400         
02401         
02402         int type=0;
02403 
02404 
02405                 Chain *c = chain_u;
02406                 
02407                 if(c->particletype)
02408                         type = c->particletype;
02409                 
02410                 for(int j=0;j<c->entries;j++)
02411                 {
02412                         if(c->parentChain==-1 && j==0)
02413                                 start = start < c->z[j]? c->z[j] : start;
02414                         if( j == c->entries-1)
02415                         {
02416                                 end = end < c->z[j] ? end : c->z[j];
02417                         }
02418                         endu=endu<c->z[j]?c->z[j]:endu;
02419 
02420                         
02421                         spoint a;
02422                         a.z=c->z[j];
02423                         a.t=c->t[j];
02424                         a.e=c->e[j];
02425                         a.chain=c->myId;
02426                         a.chainhit=j;
02427                         a.view=2;
02428                         
02429 
02430                         Managed::ManagedCluster * clu = cluster_manager->GetCluster(c->cluster_id[j]);
02431                 
02432                         if(clu)
02433                         {
02434                                 a.rms_t=clu->rms_t;
02435                                 spoints.push_back(a);
02436                         }
02437                 }
02438         
02439         
02440 
02441                 c = chain_v;
02442                 
02443                 if(c->particletype)
02444                         type = c->particletype;
02445                 for(int j=0;j<c->entries;j++)
02446                 {
02447                         if(c->parentChain==-1 && j==0)
02448                                 start = start < c->z[j]? c->z[j] : start;
02449                         if( j == c->entries-1)
02450                         {
02451                                 end = end < c->z[j] ? end : c->z[j];
02452                         }
02453                         endv=endv<c->z[j]?c->z[j]:endv;
02454                         
02455 
02456                         spoint a;
02457                         a.z=c->z[j];
02458                         a.t=c->t[j];
02459                         a.e=c->e[j];
02460                         a.chain=c->myId;
02461                         a.chainhit=j;
02462                         a.view=3;
02463                         
02464 
02465                         Managed::ManagedCluster * clu = cluster_manager->GetCluster(c->cluster_id[j]);
02466                         if(clu)
02467                         {
02468                                 a.rms_t=clu->rms_t;
02469                                 spoints.push_back(a);
02470                         }
02471                 }
02472         
02473         MSG("PrimaryShowerFinder",Msg::kDebug)<<"MAKING PARTICLE WITH "<<spoints.size()<<" POINTS\n";
02474         
02475         //we don't want the particle to extrapolate too far.... ...but now go all the way to the end
02476         double stopend = endu<endv?endv:endu;
02477 
02478                 
02479         sort(spoints.begin(), spoints.end(),spointgreaterlmf);
02480         
02481         for(unsigned int i=0;i<spoints.size();i++)
02482         {
02483 //              if( start - spoints[i].z > 0.045 )continue;
02484 //              if( spoints[i].z - end > 0.045)continue; //within 1 planes, we want the end of the chain
02485         
02486 //      printf("spoint %f %f %f %d\n",spoints[i].z,spoints[i].t,spoints[i].e,spoints[i].view);
02487         
02488                 if(spoints[i].z>stopend)continue;
02489         
02490                 int myview = spoints[i].view;
02491                 int lower=-1;
02492                 int upper=-1;
02493                 for(int j=i-1;j>-1;j--)
02494                 {
02495                         if(spoints[j].view!=myview)
02496                         {
02497                                 lower=j;
02498                                 break;
02499                         }
02500                 }
02501                 for(unsigned int j=i+1;j<spoints.size();j++)
02502                 {
02503                         if(spoints[j].view!=myview)
02504                         {
02505                                 upper=j;
02506                                 break;
02507                         }
02508                 }       
02509                 
02510                 
02511                 
02512                 double u = spoints[i].view==2 ? spoints[i].t : 0;
02513                 double v = spoints[i].view==3 ? spoints[i].t : 0;
02514                 double e = spoints[i].e;
02515                 double z = spoints[i].z;
02516                 int view = spoints[i].view;
02517                 
02518                 double rms_t = spoints[i].rms_t;
02519                 
02520 
02521                 
02522                 if(lower>-1 && upper > -1 )// good we can extrapolate!
02523                 {
02524                         double s = (spoints[upper].t - spoints[lower].t) /  ( spoints[upper].z - spoints[lower].z);
02525                         
02526                         double t = s * (spoints[i].z-spoints[lower].z) + spoints[lower].t;
02527                         
02528                         u = myview == 2 ? u : t;
02529                         v = myview == 3 ? v : t;
02530                         
02531 
02532                 }else if(lower>-1 && upper < 0)  //just user the closest other view t value
02533                 {
02534                 
02535 
02536                         //do we have another lower value?
02537                         int lower2=-1;
02538                         for(int j=lower-1;j>-1;j--)
02539                         {
02540                                 if(spoints[j].view!=myview)
02541                                 {
02542                                         lower2=j;
02543                                         break;
02544                                 }
02545                         }
02546                         
02547                         if(lower2>-1 && fabs( spoints[lower].z - spoints[lower2].z) >0)
02548                         {
02549                                 double s = (spoints[lower].t - spoints[lower2].t) /  ( spoints[lower].z - spoints[lower2].z);
02550                         
02551                                 double t = s * (spoints[i].z-spoints[lower2].z) + spoints[lower2].t;
02552                                 u = myview == 2 ? u : t;
02553                                 v = myview == 3 ? v : t;
02554                         
02555                         }else{
02556                                 u = myview == 2 ? u : spoints[lower].t;
02557                                 v = myview == 3 ? v : spoints[lower].t;
02558                         }
02559                 }
02560                 else if(upper>-1 && lower < 0)   //just user the closest other view t value
02561                 {
02562                 
02563 
02564                 
02565                         //do we have another upper value?
02566                         int upper2=-1;
02567                         for(unsigned int j=upper+1;j<spoints.size();j++)
02568                         {
02569                                 if(spoints[j].view!=myview)
02570                                 {
02571                                         upper2=j;
02572                                         break;
02573                                 }
02574                         }
02575                         
02576                         if(upper2>-1 && fabs( spoints[upper2].z - spoints[upper].z)>0)
02577                         {
02578                                 double s = (spoints[upper2].t - spoints[upper].t) /  ( spoints[upper2].z - spoints[upper].z);
02579                         
02580                                 double t = s * (spoints[i].z-spoints[upper].z) + spoints[upper].t;
02581                                 u = myview == 2 ? u : t;
02582                                 v = myview == 3 ? v : t;
02583                         
02584                         }else{
02585                                 u = myview == 2 ? u : spoints[upper].t;
02586                                 v = myview == 3 ? v : spoints[upper].t;
02587                         }
02588 
02589 
02590 
02591                         //lets use the vertex!
02592 
02593 /*
02594                         if(spoints[upper].view != spoints[i].view)
02595                                 upper=upper2;
02596                                 
02597                         if(spoints[upper].view == spoints[i].view)
02598                         {
02599 */
02600 //dont yet have a vertex!
02601 /*                              double vz =vtx_z;
02602                                 double vt = 0;
02603                                 vt = spoints[upper].view == 2 ? vtx_u : vt;
02604                                 vt = spoints[upper].view == 3 ? vtx_v : vt;
02605                         
02606                                 double t =vt;
02607                                 
02608                                 //if the spoint at z is not the same as z vtx, we need to extrapolate 
02609                                 if(fabs ( spoints[upper].z - vz)>0.001)
02610                                 {                       
02611                         
02612                                         double s = (spoints[upper].t - vt) /  ( spoints[upper].z - vz);
02613                                         t = s * (spoints[i].z-vz) + vt;                 
02614                         
02615                                 }
02616                                 
02617                 //              printf("---view %d u %f v %f t %f\n",myview,u,v,t);
02618                 //              printf("---vtx z %f t%f upper z %f t %f\n",vz,vt,spoints[upper].z,spoints[upper].t);
02619                         
02620 
02621                                 u = myview == 2 ? u : t;
02622                                 v = myview == 3 ? v : t;
02623 */
02624                                 u = myview == 2 ? u : spoints[upper].t;
02625                                 v = myview == 3 ? v : spoints[upper].t;
02626 
02627 //                      }       
02628                                 
02629 
02630                 }
02631                 else if(upper==-1 && lower==-1) //we have an empty view!!!
02632                 {
02633 
02634                         u = myview == 2 ? u : 0;
02635                         v = myview == 3 ? v : 0;
02636                 }
02637                 
02638                 foundparticle3d->add_to_back(u,v,z,e,spoints[i].chain,spoints[i].chainhit,view,rms_t);
02639                 
02640         //      printf("adding 3d spoint to particle %f %f %f --- %f\n",u,v,z,e);
02641         }
02642         
02643         
02644         if(type)
02645                 foundparticle3d->particletype=(Particle3D::EParticle3DType)type;
02646                 
02647         foundparticle3d->finalize();
02648         
02649         if(foundparticle3d->entries<1){delete foundparticle3d;foundparticle3d=0;}
02650 
02651 
02652 }

Chain * PrimaryShowerFinder::MakeShowerChain ( int  view  )  [private]

Definition at line 3377 of file PrimaryShowerFinder.cxx.

References HoughLine::cluster_id, cluster_manager, Managed::ManagedCluster::e, Managed::ClusterManager::GetCluster(), houghlineMatch, houghlinesU, houghlinesV, Managed::ManagedCluster::id, Chain::insert(), HoughLine::ncluster, HoughLine::primary, Managed::ManagedCluster::t, and Managed::ManagedCluster::z.

03378 {
03379         if(houghlineMatch.size()<1)return 0;
03380         
03381         Chain * c = new Chain();
03382 
03383         //take the first match as the best for now...
03384         HoughLine * l = 0;
03385         /*if(view ==2)l=&houghlinesU[houghlineMatch[houghlineMatch.size()-1].first];
03386         else if(view ==3)l=&houghlinesV[houghlineMatch[houghlineMatch.size()-1].second];
03387         else return 0;
03388         */
03389 
03390 /*      
03391         if(view ==2)l=&houghlinesU[houghlinesU.size()-1];
03392         else if(view ==3)l=&houghlinesV[houghlinesV.size()-1];
03393         else return 0;
03394 */      
03395 
03396         if(view ==2)l=&houghlinesU[houghlineMatch[houghlineMatch.size()-1].first];
03397         else if(view ==3)l=&houghlinesV[houghlineMatch[houghlineMatch.size()-1].second];
03398         else return 0;
03399         
03400         
03401         l->primary=1;
03402         
03403         //printf("making chain for view %d\n",view);
03404         //printf("look at houghline %d\n",view==2?houghlineMatch[houghlineMatch.size()-1].first:houghlineMatch[houghlineMatch.size()-1].second);
03405         for(int i=0;i<l->ncluster;i++)
03406         {
03407         //      printf("getting cluster %d\n",l->cluster_id[i]);
03408                 Managed::ManagedCluster *clus = cluster_manager->GetCluster(l->cluster_id[i]);
03409                 if(!clus)continue;
03410                 c->insert(clus->t,clus->z,clus->e,clus->id);
03411         //      printf("inserting %f %f %f %d\n",clus->t,clus->z,clus->e,clus->id);     
03412         }
03413 
03414 
03415         return c;
03416 }

void PrimaryShowerFinder::MergeChainClusters ( Chain ch,
std::vector< int > *  clusters 
) [private]
void PrimaryShowerFinder::Reset ( int  reset_vertex = 1  )  [private]

Definition at line 56 of file PrimaryShowerFinder.cxx.

References chain_u, chain_v, foundparticle, foundparticle3d, foundvertex, foundvertexU, foundvertexV, houghlineMatch, houghlinesU, houghlinesV, houghmapU, houghmapV, intU, intV, ran, single_view_long_shower, vtx_u, vtx_v, and vtx_z.

Referenced by FindPrimaryShower(), PrimaryShowerFinder(), and ~PrimaryShowerFinder().

00057 {
00058         ran=0;
00059         if(foundparticle3d){delete foundparticle3d;foundparticle3d=0;}
00060         foundparticle=0;
00061         chain_u=0; 
00062         chain_v=0;
00063         single_view_long_shower=0;
00064 
00065 
00066 //      if(houghmapU)delete houghmapU;
00067 //              houghmapU=0;
00068 //      if(houghmapV)delete houghmapV;
00069 //              houghmapV=0;
00070 
00071 
00072         if(houghmapU)houghmapU->Reset();
00073         if(houghmapV)houghmapV->Reset();
00074 
00075         houghlinesU.clear();
00076         houghlinesV.clear();
00077         houghlineMatch.clear();
00078 
00079         if(intU)delete intU;
00080                 intU=0;
00081         if(intV)delete intV;
00082                 intV=0;
00083 
00084         if(reset_vertex)
00085         {
00086                 vtx_u=0;
00087                 vtx_v=0;
00088                 vtx_z=0;
00089                 foundvertex=0;
00090                 foundvertexU=0;
00091                 foundvertexV=0;
00092         }
00093 }

void PrimaryShowerFinder::SaveHitGroup ( TH2D *  his,
TH2D *  saver,
double  save_val,
double  with_val,
int  curx,
int  cury 
) [private]

Definition at line 3078 of file PrimaryShowerFinder.cxx.

Referenced by MakeHoughMap().

03079 {
03080         if(curx<1 || cury < 1 || curx> his->GetNbinsX() || cury >his->GetNbinsY())return;
03081         
03082 
03083         
03084         double thisval = his->GetBinContent(curx,cury);
03085 
03086 //      printf("saving %d %d thisval %f saveval %f withval %f\n",curx,cury,thisval,save_val, with_val);
03087         if(thisval==0)return;
03088         if(fabs(thisval-save_val)<0.001)
03089         {
03090                 saver->SetBinContent(curx, cury,thisval);
03091                 //thisval=0;
03092         }
03093         else if(thisval>=with_val)return;
03094 
03095 
03096         his->SetBinContent(curx,cury,0);
03097                         
03098         SaveHitGroup(his,saver,save_val,thisval,curx+1,cury);
03099         SaveHitGroup(his,saver,save_val,thisval,curx,cury+1);
03100         SaveHitGroup(his,saver,save_val,thisval,curx-1,cury);
03101         SaveHitGroup(his,saver,save_val,thisval,curx,cury-1);
03102         SaveHitGroup(his,saver,save_val,thisval,curx+1,cury+1);
03103         SaveHitGroup(his,saver,save_val,thisval,curx-1,cury-1);
03104         SaveHitGroup(his,saver,save_val,thisval,curx+1,cury-1);
03105         SaveHitGroup(his,saver,save_val,thisval,curx-1,cury+1);
03106 //      printf("save recursion done\n");
03107         
03108 
03109         //printf("saving points %d %d %f\n",curx,cury,thisval);
03110 
03111 }

void PrimaryShowerFinder::SetVertex ( double  u,
double  v,
double  z 
) [inline]

Definition at line 56 of file PrimaryShowerFinder.h.

References foundvertex, foundvertexU, foundvertexV, vtx_u, vtx_v, and vtx_z.

Referenced by Finder::Process().

HoughLine * PrimaryShowerFinder::SplitHoughLine ( Chain c,
HoughLine hl,
Managed::ClusterManager mycm = 0 
) [private]

Definition at line 830 of file PrimaryShowerFinder.cxx.

References HoughLine::AddCluster(), HoughLine::cluster_id, cluster_manager, HoughLine::end_z, Managed::ClusterManager::GetCluster(), HoughLine::GetExpectedT(), HoughLine::ncluster, HoughLine::ResetHits(), HoughLine::start_z, Chain::t, Managed::ManagedCluster::z, and Chain::z.

Referenced by MakeChains().

00831 {
00832         //here we want to split the hough line if the chain bisects it... event by the chain projectionion
00833         
00834         if(!mycm)mycm=cluster_manager;
00835         
00836         if(!c || !hl)return 0;
00837         
00838         //check for intersection?
00839         int intersection=0;
00840         double intersection_z=0;
00841         double intersection_t=0;
00842         
00843         
00844         
00845         //consider the various intersections:
00846         
00847 //      printf("looking for intersection chain %f %f %f %f hl %f %f %f %f\n",c->start_z,c->start_t,c->end_z,c->end_t,hl->start_z,hl->start_t,hl->end_z,hl->end_t);
00848         
00849         //cross inside
00850 /*      if(c->z.size()<1)return 0;
00851         for(unsigned int i=0;i<c->z.size()-1;i++)
00852         {
00853                 double et1 = hl->GetExpectedT(c->z[i]);
00854                 double et2 = hl->GetExpectedT(c->z[i+1]);
00855                 
00856                 double max_t = et1>et2?et1:et2;
00857                 double min_t = et1<et2?et1:et2;
00858                 
00859                 double min_c_t = c->t[i]<c->t[i+1]?c->t[i]:c->t[i+1];
00860                 double max_c_t = c->t[i]>c->t[i+1]?c->t[i]:c->t[i+1];
00861                 
00862                 printf("minct %f mint %f maxct %f maxt %f\n",min_c_t,min_t,max_c_t,max_t);
00863                 
00864                 if( (min_c_t<min_t && max_c_t>max_t) || (max_t>max_c_t && min_t < min_c_t))
00865                 //intersection
00866                 {
00867                         intersection=1;
00868                         intersection_z = (c->z[i]+c->z[i+1])/2;
00869                         intersection_t = (min_c_t+max_c_t)/2;
00870                 
00871                 }
00872         }
00873 */      
00874         if(c->z.size()<1)return 0;
00875         
00876         double last_int_z=c->z[0];
00877         double last_int_dt=c->t[0]-hl->GetExpectedT(c->z[0]);
00878         
00879         for(unsigned int i=1;i<c->z.size();i++)
00880         {
00881                 double et1 = hl->GetExpectedT(c->z[i]);
00882 
00883                 double dt=c->t[i]-et1;  
00884                         
00885 //              printf("last z %f dt %f  now z %f dt %f\n",last_int_z,last_int_dt,c->z[i],dt);
00886                 
00887                 //intersection
00888                 if( ( last_int_dt<0 && dt >0 ) || ( last_int_dt >0 && dt < 0) )
00889                 {
00890                         intersection=1;
00891 
00892 /*
00893                         intersection_z = c->z[i];
00894                         intersection_t = c->t[i];
00895 */
00896 
00897                         //find the actual intersection
00898                         
00899                         double cslope=(c->z[i]-c->z[i-1]) ? (c->t[i]-c->t[i-1])/(c->z[i]-c->z[i-1]) : 0;
00900                         double coffset=c->t[i]-cslope*c->z[i];
00901                         
00902                         double ht2 = hl->GetExpectedT(c->z[i]);
00903                         double ht1 = hl->GetExpectedT(c->z[i-1]);
00904 
00905                         double hslope = (c->z[i]-c->z[i-1]) ? (ht2-ht1)/(c->z[i]-c->z[i-1]) : 0;
00906                         double hoffset = ht2 - hslope*c->z[i];
00907                         
00908                         intersection_t = (hslope-cslope) ? (hslope*coffset - cslope*hoffset)/(hslope-cslope) : 0;
00909                         if(cslope)intersection_z = (intersection_t - coffset ) / cslope;
00910                         else intersection_z=0;
00911                         
00912 
00913                         break;
00914                 }
00915                 
00916                 last_int_z=c->z[i];
00917                 last_int_dt=dt;
00918                 
00919         }
00920         
00921         
00922         //cross with projection
00923         
00924         
00925         
00926         if(!intersection)return 0;
00927 
00928 //      printf("intersection found\n");
00929 
00930         double hl_min_z = hl->start_z<hl->end_z?hl->start_z:hl->end_z;
00931         double hl_max_z = hl->start_z>hl->end_z?hl->start_z:hl->end_z;
00932         
00933         if(intersection_z < hl_min_z || intersection_z > hl_max_z)return 0;
00934         
00935         HoughLine *newHL=new HoughLine();
00936         *newHL=*hl;
00937         newHL->ResetHits();
00938         
00939         HoughLine oldHL=*hl;
00940         oldHL.ResetHits();
00941         
00942         //otherwise, split the hough line at the intersection point....
00943         for(unsigned int i=0;i<hl->cluster_id.size();i++)
00944         {
00945                 Managed::ManagedCluster *cl = mycm->GetCluster(hl->cluster_id[i]);
00946                 
00947                 if(cl->z-intersection_z<-0.01)  //require a buffer here in case the chain has two hits in the same z plane which don't have exactly same z values
00948                 {
00949                         oldHL.AddCluster(cl);
00950 //                      printf("adding to old hl  %f %f\n",cl->z,cl->t);
00951                 }
00952                 else
00953                 {
00954                         newHL->AddCluster(cl);
00955 //                      printf("adding to new hl  %f %f\n",cl->z,cl->t);
00956                 }
00957         }
00958         
00959         
00960 
00961         
00962         //its possible that we made a new houghline which has the same hits as the old houghline, but the old line now is empty....
00963         //in such a case... just keep the old hough line
00964         //this is required to prevent infitie and useless recursion
00965         
00966         if(oldHL.ncluster<1 && newHL->ncluster>0)
00967         {
00968                 *hl=*newHL;
00969                 delete newHL;
00970                 return 0;
00971         }
00972         
00973         *hl=oldHL;      
00974         return newHL;
00975 }


Member Data Documentation

Chain * PrimaryShowerFinder::chain_u = 0 [static, private]

Definition at line 85 of file PrimaryShowerFinder.h.

Referenced by FindPrimaryShower(), GetChain(), MakeParticle3D(), and Reset().

Chain * PrimaryShowerFinder::chain_v = 0 [static, private]

Definition at line 86 of file PrimaryShowerFinder.h.

Referenced by FindPrimaryShower(), GetChain(), MakeParticle3D(), and Reset().

Definition at line 50 of file PrimaryShowerFinder.h.

Referenced by FindPrimaryShower(), MakeChains(), and Finder::Process().

Definition at line 51 of file PrimaryShowerFinder.h.

Referenced by FindPrimaryShower(), MakeChains(), and Finder::Process().

Definition at line 84 of file PrimaryShowerFinder.h.

Referenced by FindPrimaryShower(), GetChain(), GetParticle3D(), and Reset().

Definition at line 127 of file PrimaryShowerFinder.h.

Referenced by MakeChains(), Reset(), and SetVertex().

Definition at line 128 of file PrimaryShowerFinder.h.

Referenced by MakeChains(), Reset(), and SetVertex().

Definition at line 129 of file PrimaryShowerFinder.h.

Referenced by MakeChains(), Reset(), and SetVertex().

std::vector<std::pair<int,int> > PrimaryShowerFinder::houghlineMatch [private]

Definition at line 116 of file PrimaryShowerFinder.h.

Referenced by FindHoughLineMatches(), MakeShowerChain(), and Reset().

TH2D * PrimaryShowerFinder::houghmapU = 0 [static, private]

Definition at line 107 of file PrimaryShowerFinder.h.

Referenced by GetHoughMap(), MakeHoughMap(), and Reset().

TH2D * PrimaryShowerFinder::houghmapV = 0 [static, private]

Definition at line 108 of file PrimaryShowerFinder.h.

Referenced by GetHoughMap(), MakeHoughMap(), and Reset().

TH2D * PrimaryShowerFinder::intU = 0 [static]

Definition at line 47 of file PrimaryShowerFinder.h.

Referenced by HoughView::DrawHough(), FindFirstIntersection(), and Reset().

TH2D * PrimaryShowerFinder::intV = 0 [static]

Definition at line 48 of file PrimaryShowerFinder.h.

Referenced by HoughView::DrawHough(), FindFirstIntersection(), and Reset().

Definition at line 53 of file PrimaryShowerFinder.h.

Referenced by HoughView::DrawHough(), FindPrimaryShower(), and Reset().

Definition at line 97 of file PrimaryShowerFinder.h.

Referenced by FindPrimaryShower(), FoundSingleViewPrimaryShower(), and Reset().

double PrimaryShowerFinder::vtx_u [private]

Definition at line 124 of file PrimaryShowerFinder.h.

Referenced by MakeChains(), Reset(), and SetVertex().

double PrimaryShowerFinder::vtx_v [private]

Definition at line 125 of file PrimaryShowerFinder.h.

Referenced by MakeChains(), Reset(), and SetVertex().

double PrimaryShowerFinder::vtx_z [private]

Definition at line 126 of file PrimaryShowerFinder.h.

Referenced by MakeChains(), Reset(), and SetVertex().


The documentation for this class was generated from the following files:

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1