Finder Class Reference

#include <Finder.h>

List of all members.

Public Member Functions

 Finder ()
 ~Finder ()
void Reset ()
void AddStrip (int myplane, int mystrip, double myenergy, double st, double sz, int view)
int GetStrips ()
void RecordLostHits (std::vector< Particle3D * > p3d)
void SetPOH (ParticleObjectHolder *p)
void Process (ParticleObjectHolder &p)
void FindIsolatedHits ()
void FindMuons (ChainHelper *ch)
void FindNeutrons (std::vector< Particle3D * > &pout)
void MergeChains (ChainHelper *ch)
void SetTrueVertex (double myu, double myv, double myz)
void FindVertex (ChainHelper *cu, ChainHelper *cv)
void MakeChains (Managed::ClusterManager *cl, ChainHelper *ch, int view)
void RemoveNonVertexPointingChains (ChainHelper *ch, int view)
void Weave (ChainHelper *chu, ChainHelper *chv)
void ClearXTalk ()
std::vector< Particle3D * > ProcessParticle3D (Particle3D *p3d)
std::vector< Particle3D * > ProcessParticle3D1 (Particle3D *p3d, int check_unique_muon=1)
std::pair< Particle3D
*, Particle3D * > 
StripMuon1 (Particle3D *p3d)
void RemoveSharedHits (std::vector< Particle3D * > pv)
std::vector< Particle3D * > SetShared (std::vector< Particle3D * > p3v)
void ShareHit (int view, int chain, int chainhit, std::vector< Particle3D * > shared)
std::vector< std::pair< int,
int > > 
matchViews (std::vector< foundpath > pu, std::vector< foundpath > pv)
std::vector< Particle3DshareEnergy (std::vector< Particle3D > p3v)
std::pair< Particle3D
*, Particle3D * > 
StripMuon (Particle3D *p3d)
void FinalizeParticles3D (std::vector< Particle3D * >pout)
void SetClusterManagerU (Managed::ClusterManager &m)
void SetClusterManagerV (Managed::ClusterManager &m)
void SetMRCC (int i)
void SetMEUperGeV (double d)

Public Attributes

double vtx_u
double vtx_v
double vtx_z
Managed::ClusterManager clustermanager_u
Managed::ClusterManager clustermanager_v

Private Member Functions

std::vector< foundpathGetPaths (ChainHelper *ch)
void DumpPaths (std::vector< foundpath > a)
Particle3D Make3DParticle (std::vector< int >upath, std::vector< int >vpath, ChainHelper *chu, ChainHelper *chv, int multu, int multv)
void DumpParticle3D (std::vector< Particle3D * > p3d)
void SetStatus (Chain *c, int status)

Private Attributes

ParticleObjectHoldermypoh
std::vector< int > view
std::vector< int > plane
std::vector< int > strip
std::vector< double > t
std::vector< double > z
std::vector< double > energy
std::vector< Particleparticles
std::multimap< double, int > sorter_map
std::map< int, std::map< int,
int > > 
loc_map
double true_vtx_u
double true_vtx_v
double true_vtx_z
double meupergev
ShareHolder shareholder
int DoMRCC

Detailed Description

Definition at line 23 of file Finder.h.


Constructor & Destructor Documentation

Finder::Finder (  ) 

Definition at line 35 of file Finder.cxx.

References DoMRCC, meupergev, and sorter_map.

00035               :  vtx_u(0.0), vtx_v(0.0), vtx_z(0.0)
00036 {
00037 
00038         sorter_map.clear();
00039         DoMRCC=0;
00040 
00041 
00042 
00043 
00044 
00045         
00046         //cluster_map.clear();
00047         meupergev=0;
00048 }

Finder::~Finder (  ) 

Definition at line 50 of file Finder.cxx.

References sorter_map.

00051 {
00052 
00053         sorter_map.clear();
00054         //cluster_map.clear();
00055 
00056 }


Member Function Documentation

void Finder::AddStrip ( int  myplane,
int  mystrip,
double  myenergy,
double  st,
double  sz,
int  view 
)

Definition at line 85 of file Finder.cxx.

References StripHolder::AddStrip(), energy, ParticleObjectHolder::event, ParticleEvent::large_maxu, ParticleEvent::large_maxv, ParticleEvent::large_maxz, ParticleEvent::large_minu, ParticleEvent::large_minv, ParticleEvent::large_minz, loc_map, mypoh, ParticleEvent::nstrips, plane, sorter_map, strip, ParticleObjectHolder::strips, t, view, and z.

Referenced by ParticleFinder::Reco().

00086 {
00087         mypoh->event.nstrips++;
00088 
00089 
00090         mypoh->strips.AddStrip(myplane,mystrip,myenergy,st,sz,myview);
00091 
00092         plane.push_back(myplane); 
00093         strip.push_back(mystrip); 
00094         energy.push_back(myenergy); 
00095         t.push_back(st); 
00096         z.push_back(sz);
00097         view.push_back(myview);
00098 
00099         mypoh->event.large_minz = (mypoh->event.large_minz < sz) ? mypoh->event.large_minz : sz;
00100         mypoh->event.large_maxz = (mypoh->event.large_maxz > sz) ? mypoh->event.large_maxz : sz;
00101         if(myview==2)
00102         {
00103                 mypoh->event.large_minu = (mypoh->event.large_minu < st) ? mypoh->event.large_minu : st;
00104                 mypoh->event.large_maxu = (mypoh->event.large_maxu > st) ? mypoh->event.large_maxu : st;
00105         }else if(myview==3){
00106                 mypoh->event.large_minv = (mypoh->event.large_minv < st) ? mypoh->event.large_minv : st;
00107                 mypoh->event.large_maxv = (mypoh->event.large_maxv > st) ? mypoh->event.large_maxv : st;        
00108         }
00109 
00110 
00111         //printf("adding %f %f %d\n",st,sz,myview);
00112     sorter_map.insert(std::pair <double,int>(myenergy, (int)energy.size()-1))
00113 ;
00114         
00115     loc_map[myplane][mystrip]=plane.size()-1;
00116 
00117 
00118 
00119 
00120 }

void Finder::ClearXTalk (  ) 

Definition at line 123 of file Finder.cxx.

References energy, exit(), loc_map, plane, sorter_map, strip, t, view, and z.

00124 {
00125         printf("DONT USE Finder::ClearXTalk()!\n");
00126         exit(1);
00127 
00128         double thresh = 3; //number of mips below which we don't care about looking for xtalk....
00129 
00130         for(unsigned int i=0;i<plane.size();i++)
00131         {
00132                 sorter_map.insert(std::pair <double,int>(energy[i],i));
00133                 loc_map[plane[i]][strip[i]]=i;
00134         }
00135         
00136 
00137         
00138         double reassignedxtalke=0.0;
00139         
00140         int foundxtalk=0;
00141         std::map<int, std::map<int, int> >::iterator p_iter;
00142         for(p_iter=loc_map.begin();p_iter!=loc_map.end(); p_iter++)
00143     {
00144         std::map<int, int>::iterator s_iter;
00145         for(s_iter=p_iter->second.begin();s_iter!=p_iter->second.end(); s_iter++)
00146         {
00147                         if (energy[s_iter->second]<thresh)continue;
00148                         
00149                         std::map<int, int>::iterator s_iter2;
00150                         for(s_iter2=p_iter->second.begin();s_iter2!=p_iter->second.end(); s_iter2++)
00151                 {
00152                                 if(s_iter==s_iter2)continue;
00153                                 int sp=abs(s_iter->first-s_iter2->first);
00154         //                      printf("sep %d view %d plane %d high %f low %f\n",sp,view[s_iter->second],plane[s_iter->second],energy[s_iter->second],energy[s_iter2->second]);
00155                                 
00156                                 if( sp>9 && sp<14)//sp == 13)
00157                                 {
00158                                         if(energy[s_iter2->second] < 0.3 * energy[s_iter->second]) //suspected crosstalk hit is < 30% of the high hit
00159                                         {
00160                                                 reassignedxtalke+=energy[s_iter2->second];
00161                                                 
00162                                                 energy[s_iter->second]+=energy[s_iter2->second];
00163                                                 energy[s_iter2->second]=0.0;
00164                                                 foundxtalk=1;
00165                                         }
00166                                 }
00167                         
00168                         }
00169                 }
00170         }
00171 
00172 
00173 //      if(reassignedxtalke>0)printf("XTALK FILTER --- %f reassigned!\n",reassignedxtalke);
00174 
00175 
00176         if(foundxtalk)
00177         {
00178                 std::vector<int>tplane;
00179                 std::vector<int>tstrip;
00180                 std::vector<int>tview;
00181                 std::vector<double>tt;
00182                 std::vector<double>tz;
00183                 std::vector<double>tenergy;
00184         
00185         
00186                 for(unsigned int i=0;i<energy.size();i++)
00187                 {
00188                         if(energy[i]>0.001)
00189                         {
00190                                 tplane.push_back(plane[i]); 
00191                                 tstrip.push_back(strip[i]); 
00192                                 tenergy.push_back(energy[i]); 
00193                                 tt.push_back(t[i]); 
00194                                 tz.push_back(z[i]);
00195                                 tview.push_back(view[i]);
00196                         }
00197                 }
00198                 
00199                 plane=tplane;
00200                 strip=tstrip;
00201                 view=tview;
00202                 t=tt;
00203                 z=tz;
00204                 energy=tenergy;
00205                 
00206         }
00207 
00208 
00209         sorter_map.clear();
00210         loc_map.clear();
00211 
00212 }

void Finder::DumpParticle3D ( std::vector< Particle3D * >  p3d  )  [private]

Definition at line 3570 of file Finder.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, ShareHolder::GetTotRemaining(), Msg::kDebug, MSG, ParticleType::muon, Particle3D::muonfrac, Particle3D::particletype, ParticleType::pi0, ParticleType::prot, Particle3D::rms_t, Particle3D::shared, shareholder, Particle3D::start_u, Particle3D::start_v, Particle3D::start_z, Particle3D::sum_e, Particle3D::types, Particle3D::u, ParticleType::uniquemuon, Particle3D::v, Particle3D::view, and Particle3D::z.

Referenced by Process(), and Weave().

03571 {
03572         ostringstream s;
03573         
03574         s << "Dump of Particle3D's "<<endl;
03575         s << "Number of Particle3D : "<<p3d.size()<< endl<<endl;
03576                 
03577                 
03578         for (unsigned int i=0;i<p3d.size();i++)
03579         {
03580                 //ostringstream s;
03581         
03582                 Particle3D * p = p3d[i];
03583                 if (p==0)continue;
03584         
03585         
03586 
03587         
03588                 s << "\n---Particle3d " << i << "---\nstart, stop (u,v,z) (" << p->start_u << ", "<< p->start_v << ", " << p->start_z <<") (" <<p->end_u<<", "<< p->end_v << ", " << p->end_z <<")"<<endl;
03589                 s << "entries "<< p->entries << " muonlike " << p->muonfrac <<endl;
03590                 
03591                 
03592 
03593                 
03594                 s<<"types: ";
03595                 for(unsigned int j=0;j<p->types.size();j++)
03596                 {
03597                         switch( p->types[j].type)
03598                         {
03599                                 case    ParticleType::em:
03600                                         s<<"em ";
03601                                         break;
03602                                 case ParticleType::emshort:
03603                                         s<<"emshort  ";
03604                                         break;                          
03605                                 case ParticleType::muon:
03606                                         s<<"muon ";
03607                                         break;
03608                                 case ParticleType::prot:
03609                                         s<<"prot ";
03610                                         break;          
03611                                 case ParticleType::pi0:
03612                                         s<<"pi0 ";
03613                                         break;
03614                                 case ParticleType::uniquemuon:
03615                                         s<<"uniquemuon ";
03616                                         break;  
03617                         
03618                         }
03619                                 
03620                 }
03621                 s<<endl;
03622                 
03623                 s<<"particletype : "<<p->particletype<<endl;
03624                 
03625                 
03626                 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;
03627                 s<<"emprob " << p->emfit_prob <<" avg rms_t "<<p->avg_rms_t<<endl;
03628                 s<<"vise "<<p->sum_e<<endl;
03629                                 
03630                 s << "points (u,v,z,e - chain, chainhit, chainview - shared - rms_t - view) : ";
03631                 for(int j=0;j<p->entries;j++)
03632                         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]<<") ";
03633         
03634         /*
03635         s << "points (e - shared) : ";
03636                 for(int j=0;j<p->entries;j++)
03637                         s << "(" <<p->e[j]<<" - "<<p->shared[j]<<") ";
03638         
03639         */
03640         
03641         
03642 
03643                         
03644                         
03645 
03646                         
03647                         
03648                 s<<endl<<endl;
03649         
03650                         //check to see if it is a good particle........
03651                 for(unsigned int q=0;q<p->u.size();q++)
03652                 {
03653                   if (p->e[q]<0 ) { MSG("Finder",Msg::kDebug) <<"BAD E "<<p->e[q]<<" at "<<q<<"\n"; }
03654                   if (p->z[q]<0  ||p->z[q]>30 || isnan(p->z[q])) { MSG("Finder",Msg::kDebug) <<"BAD Z "<<p->z[q]<<" at "<<q<<"\n"; }
03655                   
03656                   if (p->u[q]<-10||p->u[q]>10 || isnan(p->u[q])) { MSG("Finder",Msg::kDebug) <<"BAD U "<<p->u[q]<<" at "<<q<<"\n"; }
03657                   if (p->v[q]<-10||p->v[q]>10 || isnan(p->v[q])) { MSG("Finder",Msg::kDebug) <<"BAD U "<<p->v[q]<<" at "<<q<<"\n"; }
03658                 }
03659 
03660         
03661                 
03662                 
03663         }
03664         
03665         s << "\n!!! shareholder has "<<shareholder.GetTotRemaining() << " hits still shared!\n";
03666 
03667 
03668         MSG("Finder",Msg::kDebug) <<s.str();
03669 
03670         
03671 
03672 }

void Finder::DumpPaths ( std::vector< foundpath a  )  [private]

Definition at line 3186 of file Finder.cxx.

References Msg::kDebug, and MSG.

Referenced by Weave().

03187 {
03188 
03189         ostringstream os;
03190         os << "Dumping Paths\n";
03191         for(unsigned int i=0;i<a.size();i++)
03192         {
03193                 os << "Path "<<i<<" :  ";
03194                 os << " start t,z "<<a[i].start_t<<", "<<a[i].start_z;
03195                 os << " end t,z "<<a[i].end_t<<", "<<a[i].end_z;
03196                 os << " entries " <<a[i].entries << " energy " << a[i].energy << " muonlike " << a[i].muonlike;
03197                 os << "\n\t Chain Path : ";
03198                 for(unsigned int j=0;j<a[i].path.size();j++)
03199                         os<<a[i].path[j]<<" ";
03200                 os <<"\n\n";
03201         }
03202         
03203         MSG("Finder",Msg::kDebug) << os.str();  
03204 
03205 }

void Finder::FinalizeParticles3D ( std::vector< Particle3D * >  pout  ) 

Finalize each Particle3D

Finalize each Particle3D by determining final type and calculating final values based on that type

Definition at line 2694 of file Finder.cxx.

References Particle3D::calibrated_energy, ShwFit::chisq, Particle3D::Clean(), Particle3D::cmp_chisq, ShwFit::cmp_chisq, ShwFit::cmp_ndf, Particle3D::cmp_ndf, ShwFit::conv, MuELoss::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, Particle3D::entries, ShwFit::Fit3d(), ShwFit::Insert3d(), Msg::kDebug, meupergev, MSG, Particle3D::muon, n, 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, Particle3D::phi, 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, ShwFit::prob, ParticleType::prot, Particle3D::proton, Particle3D::sum_e, Particle3D::theta, Particle3D::types, Particle3D::u, Particle3D::v, vtx_u, vtx_v, vtx_z, and Particle3D::z.

Referenced by Process(), and Weave().

02695 {
02696         for(unsigned int i=0;i<pout.size();i++)
02697         {
02698                 Particle3D * p = pout[i];
02699                 if(!p)continue;
02700         
02701                 p->Clean();
02702         
02703                 
02704                 if(p->entries>0)
02705                 {               
02706                         double sum_z=0;
02707                         double sum_z_z=0;
02708                         double sum_u=0;
02709                         double sum_z_u=0;
02710                         double sum_v=0;
02711                         double sum_z_v=0;
02712                                 
02713 
02714                         double n = 5 <p->entries? 5:p->entries;
02715                         if(n<1)continue;
02716                         for(int j=0;j<n;j++)
02717                         {
02718                 
02719                                 sum_z+=p->z[j];
02720                                 sum_z_z+=p->z[j]*p->z[j];
02721                                 sum_u+=p->u[j];
02722                                 sum_z_u+=p->z[j]*p->u[j];
02723                                 sum_v+=p->v[j];
02724                                 sum_z_v+=p->z[j]*p->v[j];
02725         
02726                         }
02727                 
02728                         if(n==1)//add the vertex as a point...
02729                         {
02730                         
02731                                 sum_z+=vtx_z;
02732                                 sum_z_z+=vtx_z*vtx_z;
02733                                 sum_u+=vtx_u;
02734                                 sum_z_u+=vtx_z*vtx_u;
02735                                 sum_v+=vtx_v;
02736                                 sum_z_v+=vtx_z*vtx_v;
02737                                 n++;
02738                         }
02739                 
02740                         double denom = (n*sum_z_z-(sum_z)*(sum_z)) < 1e-10 ? 0 :  (n*sum_z_z-(sum_z)*(sum_z));
02741                         double au=0;
02742                         double bu=0;
02743                         double av=0;
02744                         double bv=0;
02745                                 
02746                         if(denom==0)            
02747                         {
02748                           au=std::numeric_limits<double>::infinity();
02749                           bu=std::numeric_limits<double>::infinity();
02750                           av=std::numeric_limits<double>::infinity();
02751                           bv=std::numeric_limits<double>::infinity();
02752 
02753 
02754 
02755                         }else{
02756                                 au=(sum_u*sum_z_z-sum_z*sum_z_u)/denom;
02757                                 bu=(n*sum_z_u-sum_z*sum_u)/denom;
02758                                 av=(sum_v*sum_z_z-sum_z*sum_z_v)/denom;
02759                                 bv=(n*sum_z_v-sum_z*sum_v)/denom;
02760                         }
02761 
02762 
02763 
02764 
02765 
02766                 //      avg_slope = b;
02767                 //      avg_offset = a;
02768                 
02769                         double theta = TMath::ATan2(bv,bu);
02770                         double phi = TMath::ATan(bu*bu+bv*bv);
02771                 
02772                 
02773                         p->theta=theta;
02774                         p->phi=phi;
02775                 
02776                 }
02777                 
02778                 
02779                 //for now... if everything else is muon energy like...
02780                 p->calibrated_energy = p->sum_e * (1.46676) / 0.0241;
02781                 
02782                 
02783                 //attempt em fit!
02784                 if(p->particletype!=Particle3D::muon && p->entries>2 && p->emfit_a==0 /*already done?*/)
02785                 {
02786                 
02787                         ShwFit f;
02788                         
02789                         double startz=0;
02790                         
02791                         startz=p->z[0];
02792                         
02793                         for(int i=0;i<p->entries;i++)
02794                         {
02795                                 //double planewidth=0.035;
02796                                 
02797                 //              f.Insert(p->e[i],(int)((p->z[i]-startz)/planewidth));
02798                         
02799                 //              cout <<"inserting "<<p->e[i]<<" at plane "<<(int)((p->z[i]-startz)/planewidth)<<endl;
02800                 
02801                 //for now, do the simple thing and assume each hit is in the next plane...
02802                         
02803                                 //                      f.Insert(p->e[i],i+1);
02804                         
02805                         
02806                                 f.Insert3d(p->e[i],p->u[i],p->v[i],p->z[i]);
02807                         
02808                                 //cout <<"inserting "<<p->e[i]<<" at plane "<<i+1<<endl;
02809                 
02810                         
02811                         }
02812                         
02813                         //f.Fit();
02814                         
02815                         f.Fit3d();
02816                         
02817                         MSG("Finder",Msg::kDebug) <<"EM FIT! " << f.par_a << " " << f.par_b << " +- "<<f.par_b_err<<" "<<f.par_e0<<endl;
02818                 
02819                         p->emfit_a=f.par_a;
02820                         p->emfit_b=f.par_b;
02821                         p->emfit_e0=f.par_e0;
02822                         p->emfit_a_err=f.par_a_err;
02823                         p->emfit_b_err=f.par_b_err;
02824                         p->emfit_e0_err=f.par_e0_err;
02825                         p->emfit_prob=f.prob;
02826                         p->calibrated_energy=f.par_e0;
02827                         p->emfit_chisq=f.chisq;
02828                         p->emfit_ndf=f.ndf;
02829                         
02830 
02831                         p->pred_e_a=f.pred_e_a; 
02832                         p->pred_g_a=f.pred_g_a;
02833                         p->pred_b=f.pred_b;
02834                         p->pred_e0=f.pred_e0;
02835                         p->pred_e_chisq=f.pred_e_chisq;
02836                         p->pred_e_ndf=f.pred_e_ndf;
02837                         p->pred_g_chisq=f.pred_g_chisq;
02838                         p->pred_g_ndf=f.pred_g_ndf;                             
02839                                                 
02840                         p->pre_over=f.pre_over;         
02841                         p->pre_under=f.pre_under;               
02842                         p->post_over=f.post_over;               
02843                         p->post_under=f.post_under;     
02844                         
02845                         p->pp_chisq=f.pp_chisq;
02846                         p->pp_ndf=f.pp_ndf;
02847                         p->pp_igood=f.pp_igood;
02848                         p->pp_p=f.pp_p;
02849                         
02850                         p->cmp_chisq = f.cmp_chisq;
02851                         p->cmp_ndf = f.cmp_ndf;
02852                         p->peakdiff = f.peakdiff;
02853         
02854 //      if(p)           printf("finder-- %f %d %f\n",p->cmp_chisq,p->cmp_ndf, p->peakdiff);             
02855                 
02856                 //      if(f.par_b<0.31 || f.par_b > 0.69) //not too emlike!
02857                 //              if(p->particletype==Particle3D::electron)p->particletype=Particle3D::other;
02858                 //      if(f.par_b>0.31 && f.par_b < 0.69) //force it to be emlike if not already
02859                 //              p->particletype=Particle3D::electron;
02860                                 
02861                 //      if(p->particletype==Particle3D::electron && f.par_e0/25/60. < 0.8) p->particletype=Particle3D::other; //<0.8gev elec is hard to find... so thats probably not it!
02862                 /*
02863                         if(TMath::Prob(p->emfit_chisq,p->emfit_chisq/p->emfit_chisq_ndf)>0.8)
02864                         {
02865                                 p->particletype=Particle3D::electron;
02866                         }
02867                         else
02868                         {
02869                                 p->particletype=Particle3D::other;
02870                         }
02871                 */
02872                 
02873                 
02874                 //      if(f.prob>0.9)  //good fit
02875                         if(f.conv && 
02876                         f.par_b - f.par_b_err*2.0 < 0.6 && f.par_b + f.par_b_err*2.0 > 0.4 //2.0 sigma from correct b value     range
02877                         && f.par_b_err < f.par_b //< 0.4
02878                         && f.par_e0/60./meupergev>0.25) //at least 0.25 gev cal e  to trust it!
02879                         {
02880                                 p->particletype=Particle3D::electron;           
02881                         }
02882                         else
02883                         {
02884                                 p->particletype=Particle3D::other;
02885                                 
02886                 //              f.Fit3dx2();
02887                         }
02888                 
02889                 }
02890         
02891                 
02892                 //proton filter
02893                 for(unsigned int i=0;i<p->types.size();i++)
02894                 {
02895                         if(p->types[i].type==ParticleType::prot)
02896                         {
02897                                 //if 80% of the particle's energy is in the last two hits....
02898                                 
02899                                 if(p->sum_e<0.001)continue;
02900                                 
02901                                 if(p->entries<3)continue;
02902                                 
02903                                 double efrac =( p->e[p->entries-1] + p->e[p->entries-2] ) / p->sum_e;
02904                                 
02905                                 if(efrac > 0.6)
02906                                 {
02907                                         p->particletype=Particle3D::proton;
02908                                         
02909                                         p->calibrated_energy=p->sum_e*60.;
02910                                         break;
02911                                 }
02912                         }
02913                 
02914                 }
02915         
02916         
02917         }
02918 }

void Finder::FindIsolatedHits (  ) 

Definition at line 1635 of file Finder.cxx.

References Particle::begPlane, Particle::begStrip, Particle::begt, Particle::begz, Particle::endPlane, Particle::endStrip, Particle::endt, Particle::endz, Particle::energy, energy, exit(), Msg::kDebug, MSG, particles, Particle::plane, plane, sorter_map, Particle::strip, strip, Particle::stripEnergy, t, Particle::type, and z.

01636 {
01637 
01638         printf("dont use Finder::FindIsolatedHits()\n");
01639         exit(1);
01640 
01641         double minthreshold = 0.5;  //peak must be at least this big
01642         double peakfraction = 0.85;  //fraction of energy in peak in this single strip
01643 
01644         //iterate until all isolated peaks are exhaused....
01645 
01646 
01647       
01648 
01649         for(std::multimap<double,int>::reverse_iterator iter = sorter_map.rbegin(); iter!=sorter_map.rend();iter++)
01650         {
01651 
01652         //for now, just look for max
01653 
01654         int maxi=-1;
01655         double maxe=0;
01656 
01657         int p=0;
01658         int c=0;
01659 
01660         double st=0;
01661         double sz=0;
01662 
01663 
01664 
01665         double locale=iter->first;
01666 
01667         maxi = iter->second;
01668         int i = maxi;
01669                         maxe=energy[i];
01670                         p=plane[i];
01671                         c=strip[i];
01672                         st=t[i];
01673                         sz=z[i];
01674         MSG("Finder",Msg::kDebug) <<"isolated hit peak at plane, cell " << p << " "<< c << " with energy " << maxe <<endl;
01675 
01676 
01677 
01678         if(maxe<minthreshold)break; //since iteration is in e order, all further entries will also be below threshold
01679 
01680 
01681 
01682 
01683         for(unsigned int i=0;i<plane.size();i++)
01684         {
01685 //              if(plane[i]==p && abs(strip[i]-c)==1)locale+=energy[i];
01686 //              if(strip[i]==c && abs(plane[i]-p)==1)locale+=energy[i];
01687 
01688                 if(plane[i]==p && abs(strip[i]-c)==1)locale+=energy[i];
01689 
01690 
01691         }
01692 
01693 
01694 
01695         MSG("Finder",Msg::kDebug) << "peak "<< maxe/locale << " "<< maxe << " "<< p << " "<< sz << " " << sz/(1.0*p) <<endl;
01696         if(maxe/locale>peakfraction) //add particle
01697         {
01698                 Particle ptemp;
01699                 ptemp.begPlane=p;
01700                 ptemp.endPlane=p;
01701                 ptemp.begStrip=c;
01702                 ptemp.endStrip=c;
01703                 ptemp.stripEnergy=maxe - (locale-maxe) /4;  //assign the energy above the local average
01704                 ptemp.type=2112; //2112 neut 2212 prot
01705                 ptemp.endz=sz;
01706                 ptemp.begz=sz;
01707                 ptemp.endt=st;
01708                 ptemp.begt=st;
01709 
01710                 ptemp.plane.push_back(p);
01711                 ptemp.strip.push_back(c);
01712                 ptemp.energy.push_back(ptemp.stripEnergy);
01713                 
01714 
01715 
01716                 MSG("Finder",Msg::kDebug) <<"p at "<< p <<" "<< c <<" "<< st <<" "<< sz<<endl;
01717 
01718                 particles.push_back(ptemp);
01719 
01720         
01721                 //adjust the energy in the peak map, removing the value assigned to this hit, and replacing it with remaining value
01722 //              sorter_map.erase((--iter).base());
01723 //              sorter_map.insert(std::pair <double,int>((locale-maxe) /4, (int)energy.size()));        
01724 
01725 
01726 
01727 
01728 
01729         }else{continue;}
01730 
01731                 
01732         }
01733 
01734 
01735 }

void Finder::FindMuons ( ChainHelper ch  ) 

Definition at line 1033 of file Finder.cxx.

References Chain::add_to_back(), Chain::e, MuELoss::e, Chain::entries, ChainHelper::FindMaxPath(), ChainHelper::GetChain(), Msg::kDebug, Msg::kError, MSG, Particle3D::muon, ChainHelper::NewChain(), ChainHelper::particles, Chain::particletype, Chain::Recalc(), Munits::second, Chain::t, and Chain::z.

01034 {
01035 
01036 
01037 //travel over chain - if we have at least two muon like hits, make a new muon chain going out to last muon hit
01038         //this new chain starts at first hit in the chain, and is not part of the primary chain
01039         //all muon like hits have all energy moved to new chain
01040         //larger hits have some average muon energy removed!
01041         
01042         
01043 
01044         std::vector<int> path = ch->FindMaxPath();
01045         if(path.size()<1)return;
01046 
01047 
01048         std::vector< std::pair<int,int> > muonlike;
01049         std::vector<double> muonlikee;
01050         
01051         int lastconsec=0;
01052         
01053         double avgmuon=0;
01054         int amuon=0;
01055         
01056         for(int i=0;i<(int)path.size();i++)
01057         {
01058         
01059                 MSG("Finder",Msg::kDebug) <<"path chain "<<path[i]<<endl;
01060                 int head=(i>0)?1:0;
01061                 for(int j=head;j<ch->GetChain(path[i])->entries;j++)
01062                 {
01063                         double e = ch->GetChain(path[i])->e[j];
01064                         if(e<2 && e>0.5){
01065                                 if((i!=(int)path.size()-1 || j!=ch->GetChain(path[i])->entries-1) || amuon>0)  
01066                                 //don't add if the only muonlike hit is at the end of a long chain
01067                                 //it could be a low energy shower hit!
01068                                 {
01069                                         muonlike.push_back(std::make_pair(path[i],j));
01070                                         muonlikee.push_back(e);
01071                                         avgmuon+=e;
01072                                         amuon++;
01073                                         lastconsec++;
01074                                 }
01075                         
01076                         }else{
01077                         
01078                           /*
01079                                 if(i==(int)path.size()-1 && j>0)
01080                                 if(e<6.0*ch->GetChain(path[i])->e[j-1])  //try not to take overlapping large hits
01081                                 {
01082                                         if(lastconsec>0)
01083                                         {
01084                                                 muonlike.push_back(std::make_pair(path[i],j));
01085                                                 muonlikee.push_back(e);
01086                                         }
01087                                 }else{
01088                                         lastconsec=0;
01089                                 }       
01090                           */
01091                           // RWH 2014-02-17 Above is what used to be here ...
01092                           // but the compiler is completely confused about dangling else 
01093                           // conditions and can't decide which to attach to which if
01094                           // this is current best guess of what (might have been)/was intended
01095                           if (i==(int)path.size()-1 && j>0) {
01096                             if (e<6.0*ch->GetChain(path[i])->e[j-1]) {
01097                               //try not to take overlapping large hits
01098                               
01099                               if (lastconsec>0) {
01100                                 muonlike.push_back(std::make_pair(path[i],j));
01101                                 muonlikee.push_back(e);
01102                               }
01103                             } else {
01104                               lastconsec=0;
01105                             }       
01106                           }
01107                           // RWH end of problematic area
01108 
01109                         }
01110                 }
01111         }
01112         
01113         if(amuon==0)return; //no muon hits found!
01114         
01115         
01116         avgmuon/=amuon;
01117         
01118 //      double avgmuon;
01119 //      for(int i=0;i<muonlikee.size();i++)
01120 //              avgmuon+=muonlikee[i];
01121 //      avgmuon/=muonlikee.size();
01122         
01123         //subtract out found muon hits
01124 //      for(int i=0;i<muonlike.size();i++)
01125 //              ch->GetChain(muonlike[i].first)->e[muonlike[i].second]=0.0;
01126         
01127         
01128         //subtract out avgmuon from all hits greater than avgmuon
01129         //make new chain, with found levels, and avglevels for rest
01130         int cid = ch->NewChain();
01131         Chain * c = ch->GetChain(cid);
01132 
01133 
01134 
01135 //      Chain c;
01136 
01137         for(unsigned int i=0;i<path.size();i++)
01138         {
01139                 int head=i>0?1:0;
01140                 
01141                 Chain * d = ch->GetChain(path[i]);
01142                 
01143                 if(d==0)
01144                 {
01145                         MSG("Finder",Msg::kError) <<"Bad path chain used - chain id "<<path[i]<<endl;
01146                         continue;
01147                 
01148                 }
01149                 
01150                 for(int j=head;j<d->entries;j++)
01151                 {
01152         //              double e = 0;
01153 /*              
01154                         if(d->e[j]>avgmuon)
01155                         {       d->e[j]-=avgmuon;
01156                                 e=avgmuon;
01157                         }else{
01158                                 e=d->e[j];
01159                                 
01160                         
01161                         
01162                         
01163                 //      double e = d->e[j];
01164                         for(int k=0;k<muonlike.size();k++)
01165                         {
01166                                 if(muonlike[k].first==path[i] && muonlike[k].second==j)
01167                                 {
01168                                         e= e > muonlikee[k] ? muonlikee[k] : e;
01169                                 }
01170                         }
01171                                 d->e[j]-=e;
01172                         }
01173 */
01174 
01175                         
01176                         //if(d->e[j]>avgmuon)d->e[j]-=avgmuon;
01177                         //double e = d->e[j];
01178                         
01179                         int found =0;
01180                         double e=0;
01181                         for(unsigned int k=0;k<muonlike.size();k++)
01182                         {
01183                                 if(muonlike[k].first==path[i] && muonlike[k].second==j)
01184                                 {
01185                                         e= muonlikee[k];
01186                                         found=1;
01187                                         break;
01188                                 }
01189                         }
01190                         if(found) //we already taged it as muon like - so move the hit
01191                         {
01192                                 //e=d->e[j];
01193                                 d->e[j]=0;
01194                         }else{
01195                                 e=avgmuon;
01196                                 d->e[j]-=avgmuon;
01197                         }
01198 
01199 
01200                 
01201                         if(e>0)
01202                         {       c->add_to_back(d->t[j],d->z[j],e,-1);
01203                                 MSG("Finder",Msg::kDebug) <<"muon adding "<< path[i] <<", "<< j<<" energy "<<e <<endl;
01204                         }
01205                 }
01206                 
01207                 
01208                 d->Recalc(); //remove the hits that we set e to 0 from the chain 
01209         }
01210         
01211         
01212 //print out chain ids....
01213 
01214 //      int cid = ch->NewChain();
01215   //  Chain * e =ch->GetChain(cid);
01216   //  *e = c;
01217         
01218         
01219         //attach new chain to current head....
01220 // //   c.parentChain=ch->GetChain(path[0]).myId;
01221 
01222         //ch->AttachAsChild(ch->GetChain(path[0]),c);
01223         
01224         c->particletype=Particle3D::muon; //say that we think it is a muon
01225         
01226         //add new chain to finished
01227 //      Chain d = *c;
01228         ch->particles.push_back(*c);    
01229         delete c;       
01230 
01231 }

void Finder::FindNeutrons ( std::vector< Particle3D * > &  pout  ) 

Definition at line 2383 of file Finder.cxx.

References Particle3D::add_to_back(), ParticleObjectHolder::chu, ParticleObjectHolder::chv, Chain::cluster_id, ParticleObjectHolder::clusterer, Managed::ClusterManager::clusters, Managed::ManagedCluster::e, ChainHelper::GetChain(), Managed::ClusterManager::GetCluster(), Managed::ClusterManager::GetViewIndex(), Managed::ManagedCluster::inuse, mypoh, Particle3D::neutron, ParticleObjectHolder::particles3d, Particle3D::particletype, Managed::ManagedCluster::t, vtx_u, vtx_v, and Managed::ManagedCluster::z.

Referenced by Weave().

02384 {
02385         std::vector<Particle3D*> pout1=pout;
02386         //mark all clusters as unused
02387         Managed::ClusterManager  *clhu = &mypoh->clusterer;
02388         Managed::ClusterManager  *clhv = &mypoh->clusterer;
02389 
02390 /*
02391         ChainHelper *cht = &mypoh->chu;
02392         for(int i=0;i<cht->finished.size();i++)
02393         {
02394                 printf("chain %d \n",cht->finished[i].myId);
02395                 for(int j=0;j<cht->finished[i].cluster_id.size();j++)
02396                         printf("%d ",cht->finished[i].cluster_id[j]);
02397                 printf("\n");
02398         
02399         }
02400 
02401 
02402         cht =& mypoh->chv;
02403         for(int i=0;i<cht->finished.size();i++)
02404         {
02405                 printf("chain %d \n",cht->finished[i].myId);
02406                 for(int j=0;j<cht->finished[i].cluster_id.size();j++)
02407                         printf("%d ",cht->finished[i].cluster_id[j]);
02408                 printf("\n");
02409         
02410         }
02411 */              
02412 
02413 
02414 //cant do this check here.. because other particles might already be in poh and are not in pout 
02415         for(unsigned int i=0;i<clhu->clusters.size();i++)clhu->clusters[i].inuse=0;
02416         for(unsigned int i=0;i<clhv->clusters.size();i++)clhv->clusters[i].inuse=0;
02417         
02418         //go through existing partcle3d's marking used clusters as used
02419 
02420         for(unsigned int i=0;i<mypoh->particles3d.size();i++)pout1.push_back(&(mypoh->particles3d[i]));
02421 
02422         for(unsigned int i=0;i<pout1.size();i++)
02423         {
02424                 for(int j=0;j<pout1[i]->entries;j++)
02425                 {
02426                         ChainHelper *ch = pout1[i]->view[j]==2? &mypoh->chu : pout1[i]->view[j]==3?& mypoh->chv:0;
02427                         if(!ch)continue;
02428                         
02429                         //printf("c %d ch %d\n",pout1[i]->chain[j],pout1[i]->chainhit[j]);
02430                         
02431                         Chain * c = ch->GetChain(pout1[i]->chain[j]);
02432                         if(!c)
02433                         {
02434                                 //cout<<"chain information lost\n";
02435                                 continue;
02436                         }
02437                         
02438                         if(pout1[i]->chainhit[j]<0)continue;  //not available
02439                         
02440                         int id = c->cluster_id[pout1[i]->chainhit[j]];
02441                         
02442                         //cout <<"chain/hit "<<pout1[i]->chain[j] << " " <<pout1[i]->chainhit[j]<<endl;
02443                         
02444                         Managed::ClusterManager  * clh = pout1[i]->view[j]==2? clhu : pout1[i]->view[j]==3? clhv:0;
02445                         if(!clh)continue;
02446                         
02447                         Managed::ManagedCluster *cluster=clh->GetCluster(id);
02448                         if(cluster)cluster->inuse=1;
02449                         //cout<<"in use\n";
02450                 
02451                 }
02452         
02453         }
02454 
02455 /*
02456 
02457         for(unsigned int i=0;i<mypoh->particles3d.size();i++)pout1.push_back(&(mypoh->particles3d[i]));
02458 
02459         for(unsigned int i=0;i<pout1.size();i++)
02460         {
02461                 for(unsigned int j=0;j<pout1[i]->entries;j++)
02462                 {
02463                         ChainHelper *ch = pout1[i]->view[j]==2? &mypoh->chu : pout1[i]->view[j]==3?& mypoh->chv:0;
02464                         if(!ch)continue;
02465                         
02466                         printf("c %d ch %d\n",pout1[i]->chain[j],pout1[i]->chainhit[j]);
02467                         
02468                         Chain * c = ch->GetChain(pout1[i]->chain[j]);
02469                         if(!c)
02470                         {
02471                                 //cout<<"chain information lost\n";
02472                                 continue;
02473                         }
02474                         
02475                         if(pout1[i]->chainhit[j]<0)continue;  //not available
02476                         
02477                         int id = c->cluster_id[pout1[i]->chainhit[j]];
02478                         
02479                         cout <<"chain/hit "<<pout1[i]->chain[j] << " " <<pout1[i]->chainhit[j]<<endl;
02480                         
02481                         Managed::ClusterManager  * clh = pout1[i]->view[j]==2? clhu : pout1[i]->view[j]==3? clhv:0;
02482                         if(!clh)continue;
02483                         
02484                         Managed::ManagedCluster *cluster=clh->GetCluster(id);
02485                         if(cluster)cluster->inuse=1;
02486                         cout<<"in use\n";
02487                 
02488                 }
02489         
02490         }
02491 
02492 
02493 */
02494 
02495         //find unused clusters with energy > some thresh
02496         //make these into particles...
02497 /*      
02498         cout << "//////////////////////////////////////////////\nunused clutsers"<<endl;
02499         cout <<"U\n";
02500         
02501         for(unsigned int i=0;i<clhu->clusters.size();i++)
02502         {
02503                 Managed::ManagedCluster *c = & clhu->clusters[i];
02504                 if(c->inuse)continue;
02505                 cout << "Cluster "<<c->id << " z " << c->z << " t " << c->t << " E " << c->e<<endl;
02506         }
02507         cout <<"V\n";
02508         for(unsigned int i=0;i<clhv->clusters.size();i++)
02509         {
02510                 Managed::ManagedCluster *c = & clhv->clusters[i];
02511                 if(c->inuse)continue;
02512                 cout << "Cluster "<<c->id << " z " << c->z << " t " << c->t << " E " << c->e<<endl;
02513         }
02514         
02515         cout << "//////////////////////////////////////////////\n";
02516 
02517 */
02518 /*
02519         //for now just make 2d neutrons... (don't match views..)
02520         
02521         double neut_thresh=3.0;
02522         
02523         for(int q=0;q<2;q++)
02524         {
02525                 Managed::ClusterManager * ch = q==0?clhu:clhv;
02526                 for(unsigned int i=0;i<ch->clusters.size();i++)
02527                 {
02528                         Managed::ManagedCluster *c = & ch->clusters[i];
02529                         if(c->inuse)continue;
02530                         if(c->e<neut_thresh)continue;
02531                         
02532                                 Particle3D * pnew = new Particle3D();
02533                                 pnew->particletype=Particle3D::neutron;
02534                                 double u = q==0? c->t : 0.0;
02535                                 double v = q==1? c->t : 0.0;
02536                                 pnew->add_to_back(u, v, c->z, c->e, -1,-1, q==0?2:3, 0);
02537                                 pout.push_back(pnew);                   
02538                 }
02539         
02540         }
02541 */
02542         //for now just make 2d neutrons... (don't match views..)
02543         
02544         double neut_thresh=3.0;
02545         
02546         for(int q=0;q<2;q++)
02547         {
02548                 //Managed::ClusterManager * ch = q==0?clhu:clhv;
02549                 
02550                 std::vector<int> chidx = clhu->GetViewIndex(q==0?2:3);
02551                 std::vector<int> chidx_other = clhu->GetViewIndex(q==1?2:3);
02552                 
02553                 for(unsigned int i=0;i<chidx.size();i++)
02554                 {
02555                         Managed::ManagedCluster *c = & clhu->clusters[chidx[i]];
02556                         if(c->inuse)continue;
02557                         if(c->e<neut_thresh)continue;
02558                         
02559                         
02560                         //Managed::ClusterManager * ch_other = q==1?clhu:clhv;
02561                         std::vector<int>other_view;
02562                         for(unsigned int j=0;j<chidx_other.size();j++)
02563                         {
02564                                 Managed::ManagedCluster *cother = & clhu->clusters[chidx_other[j]];
02565                                 if(TMath::Abs(cother->z - c->z)<0.06 && !cother->inuse)
02566                                         other_view.push_back(j);
02567                         }       
02568                                 
02569                         if(other_view.size()==0)
02570                         {       
02571                                 Particle3D * pnew = new Particle3D();
02572                                 pnew->particletype=Particle3D::neutron;
02573                                 double u = q==0? c->t : vtx_u;
02574                                 double v = q==1? c->t : vtx_v;
02575                                 pnew->add_to_back(u, v, c->z, c->e, -1,-1, q==0?2:3, 0);
02576                                 pout.push_back(pnew);   
02577                                 c->inuse=1;             
02578                         }else{
02579                                 //don't interpolate....
02580                                 
02581                                 Managed::ManagedCluster *cother0 = & clhu->clusters[chidx_other[other_view[0]]];
02582                                 
02583                                 double o_t = cother0->t;
02584                                 if(other_view.size()==2) 
02585                                 {
02586                                         Managed::ManagedCluster *cother1 = & clhu->clusters[chidx_other[other_view[1]]];
02587 
02588                                         o_t += cother1->t;
02589                                         o_t /=2.0;
02590                                 }
02591                                 
02592                                 
02593                                 double u = q==0? c->t : o_t;
02594                                 double v = q==1? c->t : o_t;                            
02595                         
02596                         
02597                                 Particle3D * pnew = new Particle3D();
02598                                 pnew->particletype=Particle3D::neutron;
02599                                 pnew->add_to_back(u, v, c->z, c->e, -1,-1, q==0?2:3, 0);
02600                                 c->inuse=1;
02601                                 pnew->add_to_back(u, v, cother0->z, cother0->e, -1,-1, q==1?2:3, 0);
02602                                 cother0->inuse=1;                               
02603                                 if(other_view.size()==2)
02604                                 {
02605                                         Managed::ManagedCluster *cother1 = & clhu->clusters[chidx_other[other_view[1]]];
02606                                         pnew->add_to_back(u, v, cother1->z, cother1->e, -1,-1, q==1?2:3, 0);
02607                                         cother1->inuse=1;
02608                                 }
02609                                 
02610                                 pout.push_back(pnew);   
02611                         }
02612         
02613         
02614                 }
02615         
02616         }
02617 
02618 
02619 }

void Finder::FindVertex ( ChainHelper cu,
ChainHelper cv 
)

Definition at line 1290 of file Finder.cxx.

References Chain::children, Chain::e, Chain::end_z, Chain::entries, ChainHelper::FindMaxPath(), ChainHelper::FindPaths(), ChainHelper::finished, ChainHelper::GetChain(), Chain::interpolate(), Msg::kDebug, Chain::level, Munits::m, MSG, Chain::myId, Chain::parentChain, Chain::Reverse(), ChainHelper::SplitAt(), Chain::start_t, Chain::start_z, Chain::t, vtx_u, vtx_v, vtx_z, and Chain::z.

Referenced by Process().

01291 {
01292         std::vector<int> mpu = cu->FindMaxPath();
01293         std::vector<int> mpv = cv->FindMaxPath();
01294         
01295         if(mpu.size() <1 || mpv.size() <1)return;
01296         
01297         
01298         double cu_t = cu->GetChain(mpu[0])->start_t;
01299         double cu_z = cu->GetChain(mpu[0])->start_z;
01300 
01301         double cv_t = cv->GetChain(mpv[0])->start_t;
01302         double cv_z = cv->GetChain(mpv[0])->start_z;    
01303         
01304         
01306         //try iterating through, measuring dt between adjacent hits in chain.
01307         //go through until dt(0,1)-dt(1,2)<thresh
01308         
01309         
01310 /*      int forcedadvancedu=0;
01311         int forcedadvancedv=0;
01312 */
01313         
01314 
01315 MSG("Finder",Msg::kDebug)<< "cuz " << cu_z <<" cvz "<< cv_z <<endl;
01316         
01317         
01319 /*      
01320         if((forcedadvancedu && forcedadvancedv) || (!forcedadvancedu && !forcedadvancedv))
01321         {
01322         
01323                 if(cv_z<cu_z)
01324                         cu_z=cv_z;
01325                 else
01326                         cv_z=cu_z;
01327         }else{
01328                 if(forcedadvancedu)
01329                         cv_z=cu_z;
01330                 if(forcedadvancedv)
01331                         cu_z=cv_z;
01332         }
01333 */              
01334 
01335 
01336 
01337                 if(cv_z<cu_z)
01338                         cu_z=cv_z;
01339                 else
01340                         cv_z=cu_z;
01341                         
01342 
01343                 int chuid=0;
01344                 
01345         //      if(cu->GetChain(mpu[0]).start_z<cu_z)
01346                 for(unsigned int i=0;i<mpu.size();i++)
01347                 {
01348                         if ( ( cu->GetChain(mpu[i])->start_z <= cu_z && cu->GetChain(mpu[i])->end_z >= cu_z )
01349                                  || (cu->GetChain(mpu[i])->start_z >= cu_z && cu->GetChain(mpu[i])->end_z <= cu_z))
01350                         {
01351                                 chuid=i;
01352                                 
01353                         }
01354                 }
01355                 
01356                 if(cu->GetChain(mpu[0])->entries<2 && mpu.size()>1)chuid=1;
01357                 
01358                 
01359                 int chvid=0;
01360                 
01361         //      if(cv->GetChain(mpv[0]).start_z<cv_z)
01362                 for(unsigned int i=0;i<mpv.size();i++)
01363                 {
01364                         if ( ( cv->GetChain(mpv[i])->start_z <= cv_z && cv->GetChain(mpv[i])->end_z >= cv_z )
01365                                  || (cv->GetChain(mpv[i])->start_z >= cv_z && cv->GetChain(mpv[i])->end_z <= cv_z))
01366                         {
01367                                 chvid=i;
01368                                 
01369                         }
01370                 }
01371         
01372                 if(cv->GetChain(mpv[0])->entries<2 && mpv.size()>1)chvid=1;
01373                 //printf("chain has %d entries, path has %d entries\n",cv->GetChain(mpv[0])->entries,mpv.size());
01374 
01375                 /*
01376                 printf("uchain\n");
01377                 for(unsigned int i=0;i<mpu.size();i++)printf("%d\n",mpu[i]);
01378                 printf("vchain\n");
01379                 for(unsigned int i=0;i<mpv.size();i++)printf("%d\n",mpv[i]);
01380         */      
01381                 
01382                         Chain * uuu = cu->GetChain(mpu[chuid]);
01383                         Chain * vvv =  cv->GetChain(mpv[chvid]);
01384                 
01385                 
01386                 
01387                         cu_z = uuu->start_z;
01388                         cv_z = vvv->start_z;
01389                 
01390                         if(cv_z<cu_z)
01391                         {
01392                                 cu_z=cv_z;
01393                         }
01394                         else
01395                         {
01396                                 cv_z=cu_z;
01397 
01398                         }
01399                 
01400                 
01401                         double tmpuz=0;
01402                         double tmpvz=0;
01403                         
01404                 //      double slopeu = cu->GetChain(mpu[chuid])->avg_slope;
01405                 //      double bzu = cu->GetChain(mpu[chuid])->start_z;
01406                 //      double btu = cu->GetChain(mpu[chuid])->start_t;
01407                         //double bt = cu->GetChain(mpu[chuid])->weighted_t;
01408                         
01409                                         //see if there is a big energy difference between the first 2 hits... if so, we probably want the 2nd hit
01410                                         if(uuu->entries>1 && uuu->e[0] < uuu->e[1]*0.2 )
01411                                         {
01412                                                 //bz = uuu->z[1];
01413                                                 //bt = uuu->t[1];
01414                                                 tmpuz=uuu->z[1];
01415                                                 cu_z=uuu->z[1];
01416                                         }
01417                                         else if(uuu->entries>2 && uuu->e[1] < uuu->e[2]*0.2 )
01418                                         {
01419                                                 //bz = uuu->z[2];
01420                                                 //bt = uuu->t[2];
01421                                                 tmpuz=uuu->z[2];
01422                                                 cu_z=uuu->z[1];
01423                                         }
01424                         
01425                         
01426 
01427 
01428         
01429         
01430                         
01431 //                      double slopev = cv->GetChain(mpv[chvid])->avg_slope;
01432 //                      double bzv = cv->GetChain(mpv[chvid])->start_z;
01433 //                      double btv = cv->GetChain(mpv[chvid])->start_t;
01434 
01435                         //double slopev = cv->GetChain(mpv[chvid])->front_slope;
01436                         //double offsetv = cv->GetChain(mpv[chvid])->front_offset;
01437                         //double slopeu = cu->GetChain(mpu[chuid])->front_slope;
01438                         //double offsetu = cu->GetChain(mpu[chuid])->front_offset;
01439 
01440                                                 
01441                                         if(vvv->entries>1 && vvv->e[0] < vvv->e[1]*0.2 )
01442                                         {
01443                                                 //bz = vvv->z[1];
01444                                                 //bt = vvv->t[1];
01445                                                 tmpvz=vvv->z[1];
01446                                                 cv_z=vvv->z[1];
01447                                         }
01448                                         else if(vvv->entries>2 && vvv->e[1] < vvv->e[2]*0.2 )
01449                                         {
01450                                         //      bz = vvv->z[2];
01451                                         //      bt = vvv->t[2];
01452                                                 tmpvz=vvv->z[2];
01453                                                 cv_z=vvv->z[2];
01454                                         }
01455 
01456                         //bt = cv->GetChain(mpv[chvid])->weighted_t;
01457                         
01458                 
01460 
01461                 
01462         vtx_z=cu_z<cv_z?cu_z:cv_z;
01463         
01464 //      if(tmpuz>0 && tmpvz == 0){ vtx_z=tmpuz; cu_z = vtx_z; }
01465 //      if(tmpvz>0 && tmpuz == 0){ vtx_z=tmpvz; cv_z = vtx_z;}
01466 //      if(tmpuz>0 && tmpvz > 0){ vtx_z=tmpuz<tmpvz ? tmpuz:tmpvz;  cu_z=vtx_z; cv_z=vtx_z; }           
01467                 
01468         //vtx_z-= 0.04;         
01469         
01470         
01471                 
01472                         //if we have a hit at the right z spot, use that t, otherwise extrapolate
01473                         if(fabs(vtx_z-uuu->z[0])<0.01)
01474                         {
01475                                 cu_t = uuu->t[0];
01476                         }else{
01477                                 
01478                                 double oldcut=cu_t;                     
01479                         //      cu_t = slopeu*(-bzu+cu_z) + btu;
01481                 
01482                         cu_t = uuu->interpolate(vtx_z);
01483                 
01484                         
01485                                 MSG("Finder",Msg::kDebug) << "using chain "<< mpu[chuid]<<" "<<cu->GetChain(mpu[chuid])->myId <<" for slope, moving ut from "<< oldcut <<" to "<< cu_t <<endl;
01486                         }
01487         
01488 
01489                         if(fabs(vtx_z-vvv->z[0])<0.01)
01490                         {
01491                                 cv_t = vvv->t[0];
01492                         }else{
01493                                                 
01494                                 double oldcut=cv_t;
01495                                 //cv_t = slopev*(-bzv+cv_z) + btv;              
01496 
01497                                 //cv_t = slopev * vtx_z + offsetv;
01498 
01499                         cv_t = vvv->interpolate(vtx_z);
01500 
01501 
01502                                 MSG("Finder",Msg::kDebug) << "using chain "<< mpv[chvid]<<" "<<cv->GetChain(mpv[chvid])->myId <<" for slope, moving ut from "<< oldcut <<" to "<< cv_t <<endl;
01503                         }
01504         
01505         
01506 
01507         
01508         //we should see if there is a chain path that passes close to this vertex
01509         //that is, if the found vertex is between the beginning and end of this chain, the chain should be split at the vertex, have the front part reversed
01510         //and the vertex should be adjusted to fall on this chain
01511         
01512         for(int q=0;q<2;q++)
01513         {
01514         
01515         ChainHelper * ch = q==0? cu:cv;
01516         
01517         for(unsigned int i=0;i<ch->finished.size();i++)
01518         {
01519                 if(ch->finished[i].start_z > vtx_z) continue;
01520                 std::vector<foundpath> paths = ch->FindPaths(ch->GetChain(ch->finished[i].myId));
01521                 
01522                 //find the chain where we would need to split....
01523                 for(unsigned int p =0;p<paths.size();p++)
01524                 {
01525         
01526                         //make sure the path is sufficiently far from the split point
01527                         if(! ( paths[p].start_z < vtx_z - 0.04 && paths[p].end_z > vtx_z + 0.04) )break;
01528         
01529                         int getp=-1;
01530                         for(unsigned int j=0;j<paths[p].path.size();j++)
01531                         {
01532                                 getp = paths[p].path[j];
01533                                 if(ch->GetChain(getp)->start_z < vtx_z -0.03 && ch->GetChain(getp)->end_z >= vtx_z+0.03)break;
01534                                 getp=-1;
01535                         }
01536                         if(getp<0)continue;
01537 
01538                 
01539         //              cout <<"WANT TO SPLIT "<<getp << " AT " <<vtx_z<<endl;
01540                 
01541                         Chain * m = ch->GetChain(getp);
01542                         //if its close to the vertex, ajust vertex position
01543                         
01544                         int pt_bef=-1;
01545                         int pt_aft=-1;
01546                         for(int j=1;j<m->entries;j++)
01547                         {
01548                                 if(m->z[j-1]<=vtx_z && m->z[j]>=vtx_z)
01549                                 {
01550                                         pt_bef=j-1;
01551                                         pt_aft=j;
01552                                         break;
01553                                 }
01554                         }
01555                 
01556 //                      cout<<"recalculating t from idx "<<pt_bef<<" "<<pt_aft<<endl;
01557                 
01558                         if(pt_bef>-1 && pt_aft>-1)
01559                         {
01560                                 double dt = m->t[pt_aft]-m->t[pt_bef];
01561                                 double dz = m->z[pt_aft]-m->z[pt_bef];
01562                         
01563                 //              cout << "dt "<<dt <<" dz "<<dz<<endl;
01564                         
01565                                 double vt=0;
01566                                 if(dt<0.0001)vt=m->t[pt_aft];
01567                                 
01568                                 else vt = (vtx_z-m->z[pt_bef]) * dt/dz + m->t[pt_bef];
01569                                 
01570                                 double oldt = q==0?cu_t:cv_t;
01571                                 
01572                 //              cout <<"old t "<<oldt << " new t "<<vt<<endl;
01573                                 
01574                                 if(fabs(oldt-vt)<0.05)
01575                                 {
01576                                         if(q==0)cu_t=vt;
01577                                         if(q==1)cv_t=vt;
01578                                 }
01579                         }
01580                 
01581                         
01582                         
01583                         Chain * d = ch->SplitAt(m,vtx_z);
01584                         
01585                         if(d->children.size()<1)continue; //unable to split the chain...
01586                         
01587                         //d->parentChain=-1;
01588                         d->level=0;
01589                         for(unsigned int chi =0;chi<d->children.size();chi++)
01590                         {
01591                                 Chain * cld = ch->GetChain(d->children[chi]);
01592                                 if(cld){
01593                                         cld->parentChain=-1;
01594                                 }
01595                         }
01596                         
01597                         Chain*cp = ch->GetChain(d->parentChain);
01598                         if(cp){
01599                                         for(unsigned int i=0;i<cp->children.size();i++)
01600                                         {
01601                                         if(cp->children[i]==d->myId)
01602                                         {
01603                                                 cp->children.erase(cp->children.begin()+i);
01604                         
01605                         //                      cout<<"\t reversal causing removing of child reference from parent "<<d->parentChain<<endl;
01606                         
01607                                                 break;  
01608                                         }
01609                                         }       
01610                         }
01611                         d->parentChain=-1;
01612                         
01613                         d->children.clear();
01614                         d->Reverse();
01615                         
01616                 }
01617         
01618         }
01619         
01620         
01621         }
01622         
01623 
01624         
01625 //      vtx_z-= 0.04;//0.025;  //event is usually before first hit.... make it better later, based on if this looks like a proton (indicating event forward of strip), or if this is a low hit strip(indicating hit well inside steel)
01626         
01627         vtx_u=cu_t;
01628         vtx_v=cv_t;
01629 
01630 }

std::vector< foundpath > Finder::GetPaths ( ChainHelper ch  )  [private]

Definition at line 3166 of file Finder.cxx.

References ChainHelper::FindPaths(), ChainHelper::finished, and t.

Referenced by Weave().

03167 {
03168 
03169 
03170         std::vector<foundpath> p;
03171         for(unsigned int i=0;i<ch->finished.size();i++)
03172         {
03173                 if(ch->finished[i].parentChain<0 )
03174                 {
03175                         std::vector<foundpath> t = ch->FindPaths(&ch->finished[i]);
03176                         //printf("head chain %d has %d paths\n",ch->finished[i].myId, t.size());
03177                         
03178                         for(unsigned int j=0;j<t.size();j++)
03179                         p.push_back(t[j]);
03180         
03181                 }
03182          }
03183          return p;
03184 }

int Finder::GetStrips (  )  [inline]

Definition at line 34 of file Finder.h.

References plane.

Referenced by ParticleFinder::Reco().

00034 {return plane.size();};

Particle3D Finder::Make3DParticle ( std::vector< int >  upath,
std::vector< int >  vpath,
ChainHelper chu,
ChainHelper chv,
int  multu = 0,
int  multv = 0 
) [private]

printf("LARGEST Z %f t %f e %f\n",largest_z,points[largest_idx].t,points[largest_idx].e);

Definition at line 3220 of file Finder.cxx.

References MuELoss::a, Particle3D::add_to_back(), Chain::available, point::chain, point::chainhit, Chain::cluster_id, ParticleObjectHolder::clusterer, done(), Managed::ManagedCluster::e, point::e, MuELoss::e, Chain::entries, Particle3D::finalize(), ChainHelper::GetChain(), Managed::ClusterManager::GetCluster(), Managed::ManagedCluster::id, Chain::myId, mypoh, Chain::parentChain, Particle3D::particletype, Chain::particletype, pointgreater(), point::rms_t, Managed::ManagedCluster::rms_t, point::t, Chain::t, t, view, point::view, vtx_u, vtx_v, vtx_z, z, point::z, and Chain::z.

Referenced by Weave().

03221 {
03222         Particle3D p;
03223 
03224 
03225         
03226         std::vector<point> points;
03227         
03228         double start =0;
03229         double end  =1000;
03230         
03231         double endu=0;
03232         double endv=0;
03233                 
03234         
03235         
03236         
03237         int type=0;
03238         
03239         
03240         for(unsigned int i=0;i<upath.size();i++)
03241         {
03242                 Chain *c = chu->GetChain(upath[i]);
03243                 //cout <<"using chain "<<upath[i]<<"\n";
03244                                 
03245                 if(c->particletype)
03246                         type = c->particletype;
03247                 
03248                 for(int j=0;j<c->entries;j++)
03249                 {
03250                         if(c->parentChain==-1 && j==0)
03251                                 start = start < c->z[j]? c->z[j] : start;
03252                         if(i == upath.size()-1 && j == c->entries-1)
03253                         {
03254                                 end = end < c->z[j] ? end : c->z[j];
03255                         }
03256                         endu=endu<c->z[j]?c->z[j]:endu;
03257 
03258                         //don't double count children head hits
03259                         if(c->parentChain>-1 && j==0)
03260                         {
03261                                 if(i==0)continue;
03262                                 Chain *cpast = chu->GetChain(upath[i-1]);
03263                                 if(cpast->cluster_id[cpast->entries-1] == c->cluster_id[0])
03264                                         continue; 
03265                         }
03266                         Managed::ManagedCluster * clu = mypoh->clusterer.GetCluster(c->cluster_id[j]);  
03267                         if(!clu)continue;               
03268                 
03269                         
03270                         point a;
03271                         a.z=c->z[j];
03272                         a.t=c->t[j];
03273                         if(c->available && clu->id>0)
03274                                 a.e=clu->e;
03275                         else a.e=0;
03276                         a.chain=c->myId;
03277                         if(c->available && clu->id>0)
03278                                 a.chainhit=j;
03279                         else
03280                                 a.chainhit=-1;
03281                         a.view=2;
03282                         
03283 
03284                         if(clu)
03285                         {
03286                                 a.rms_t=clu->rms_t;
03287                                 points.push_back(a);
03288                         }
03289                         
03290                         
03291                 //      printf("adding pt u %f %f %f\n",a.z,a.t,a.e);
03292                         
03293                 }
03294         
03295         }
03296         for(int i=0;i<(int)vpath.size();i++)
03297         {
03298                 Chain *c = chv->GetChain(vpath[i]);
03299                 
03300                 //cout <<"using chain "<<vpath[i]<<"\n";
03301                 
03302                 if(c->particletype)
03303                         type = c->particletype;
03304                 for(int j=0;j<c->entries;j++)
03305                 {
03306                         if(c->parentChain==-1 && j==0)
03307                                 start = start < c->z[j]? c->z[j] : start;
03308                         if(i == (int)vpath.size()-1 && j == c->entries-1)
03309                         {
03310                                 end = end < c->z[j] ? end : c->z[j];
03311                         }
03312                         endv=endv<c->z[j]?c->z[j]:endv;
03313                         
03314                         //don't double count children head hits
03315                         if(c->parentChain>-1 && j==0)
03316                         {
03317                                 if(i==0)continue;
03318                                 Chain *cpast = chv->GetChain(vpath[i-1]);
03319                                 if(cpast->cluster_id[cpast->entries-1] == c->cluster_id[0])
03320                                         continue; 
03321                         }
03322 
03323                         Managed::ManagedCluster * clu = mypoh->clusterer.GetCluster(c->cluster_id[j]);  
03324                         if(!clu)continue;               
03325                         
03326 
03327                         point a;
03328                         a.z=c->z[j];
03329                         a.t=c->t[j];
03330                         if(c->available && clu->id>0)
03331                                 a.e=clu->e;
03332                         else
03333                                 a.e=0;
03334                         a.chain=c->myId;
03335                         if(c->available && clu->id>0)
03336                                 a.chainhit=j;
03337                         else
03338                                 a.chainhit=-1;
03339                         a.view=3;
03340                         
03341 
03342                         if(clu){
03343                                 a.rms_t=clu->rms_t;
03344                                 points.push_back(a);
03345                         }
03346 
03347                 //      printf("adding pt v %f %f %f\n",a.z,a.t,a.e);
03348 
03349                 }
03350         
03351         }
03352         
03353         //we don't want the particle to extrapolate too far....
03354         double stopend = endu<endv?endu:endv;
03355         //but if that chain path is only used once, then we want to go all the way to the end!
03356         if(!multu && multv)stopend = stopend<endu?endu:stopend;
03357         if(!multv && multu)stopend = stopend<endv?endv:stopend;
03358         if(!multv && !multu)stopend = endu<endv?endv:endu;
03359                 
03360         std::sort(points.begin(), points.end(),pointgreater);
03361 
03362 /*
03363 double largest_z=0;
03364 double largest_idx=0;
03365         for(unsigned int i=0;i<points.size();i++)
03366         {
03367                 if(points[i].z>largest_z)
03368                 {
03369                         largest_z=points[i].z;
03370                         largest_idx=i;
03371                 }
03372         }
03373 */      
03375 
03376         
03377         for(unsigned int i=0;i<points.size();i++)
03378         {
03379 //              if( start - points[i].z > 0.045 )continue;
03380 //              if( points[i].z - end > 0.045)continue; //within 1 planes, we want the end of the chain
03381         
03382 
03383         
03384         
03385         
03386         //are we at the end of usuable points?
03387         int done=0;
03388         for(unsigned int j=i;j<points.size();j++)
03389         {
03390                 if(points[j].e>0)break;
03391         
03392                 if(j==points.size()-1)done=1;
03393         }
03394         if(done)break;
03395         
03396                 if(points[i].z>stopend)continue;
03397         
03398                 int myview = points[i].view;
03399                 int lower=-1;
03400                 int upper=-1;
03401                 for(int j=i-1;j>-1;j--)
03402                 {
03403                         if(points[j].view!=myview)
03404                         {
03405                                 lower=j;
03406                                 break;
03407                         }
03408                 }
03409                 for(unsigned int j=i+1;j<points.size();j++)
03410                 {
03411                         if(points[j].view!=myview)
03412                         {
03413                                 upper=j;
03414                                 break;
03415                         }
03416                 }       
03417                 
03418                 
03419                 
03420                 double u = points[i].view==2 ? points[i].t : 0;
03421                 double v = points[i].view==3 ? points[i].t : 0;
03422                 double e = points[i].e;
03423                 double z = points[i].z;
03424                 int view = points[i].view;
03425                 
03426                 double rms_t = points[i].rms_t;
03427                 
03428 
03429                 
03430                 if(lower>-1 && upper > -1 )// good we can extrapolate!
03431                 {
03432                         double s = (points[upper].t - points[lower].t) /  ( points[upper].z - points[lower].z);
03433                         
03434                         double t = s * (points[i].z-points[lower].z) + points[lower].t;
03435                         
03436                         u = myview == 2 ? u : t;
03437                         v = myview == 3 ? v : t;
03438                         
03439 
03440                 }else if(lower>-1 && upper < 0)  //just user the closest other view t value
03441                 {
03442                         
03443                 //      printf("looking for lower value\n");
03444                         //do we have another lower value?
03445                         int lower2=-1;
03446                         for(int j=lower-1;j>-1;j--)
03447                         {
03448                                 if(points[j].view!=myview && fabs(points[lower].z-points[j].z)>0.04)
03449                                 {
03450                                         lower2=j;
03451                 //                      printf("second lower found\n");
03452                                         break;
03453                                 }
03454                         }
03455                         
03456                         
03457                         if(lower2>-1 && fabs( points[lower].z - points[lower2].z) >0)
03458                         {
03459                                 double s = (points[lower].t - points[lower2].t) /  ( points[lower].z - points[lower2].z);
03460                         
03461                                 double off = points[lower].t - points[lower].z*s;
03462                         
03463                                 double t = s * points[i].z + off;
03464                                 u = myview == 2 ? u : t;
03465                                 v = myview == 3 ? v : t;
03466                 //              printf("extrap  tz tz   %f %f    %f %f to %f\n",points[lower].t,points[lower].z,points[lower2].t,points[lower2].z,points[i].z);
03467                         }else{
03468                                 u = myview == 2 ? u : points[lower].t;
03469                                 v = myview == 3 ? v : points[lower].t;
03470         //                      printf("using closer point\n");
03471                         }
03472         //              printf("uv %f %f\n",u,v);
03473                         
03474                 }
03475                 else if(upper>-1 && lower < 0)   //just user the closest other view t value
03476                 {
03477                 
03478 
03479                 
03480                         //do we have another upper value?
03481                         int upper2=-1;
03482                         for(unsigned int j=upper+1;j<points.size();j++)
03483                         {
03484                                 if(points[j].view!=myview && fabs(points[upper].z-points[j].z)>0.04)
03485                                 {
03486                                         upper2=j;
03487                                         break;
03488                                 }
03489                         }
03490                         
03491                         if(upper2>-1 && fabs( points[upper2].z - points[upper].z)>0)
03492                         {
03493                                 double s = (points[upper2].t - points[upper].t) /  ( points[upper2].z - points[upper].z);
03494                                 double off = points[upper].t - points[upper].z*s;
03495                         
03496                                 double t = s * points[i].z + off;
03497                                 u = myview == 2 ? u : t;
03498                                 v = myview == 3 ? v : t;
03499                                                                 
03500                 //              printf("ext next hit with t %f z %f  t %f z %f for ext\n",points[upper].t,points[upper].z,points[upper2].t,points[upper2].z);
03501                         
03502                         }else{
03503                                 u = myview == 2 ? u : points[upper].t;
03504                                 v = myview == 3 ? v : points[upper].t;
03505 
03506                         }
03507 
03508 
03509 
03510                         //lets use the vertex!
03511 
03512 /*
03513                         if(points[upper].view != points[i].view)
03514                                 upper=upper2;
03515                                 
03516                         if(points[upper].view == points[i].view)
03517                         {
03518 */                              double vz =vtx_z;
03519                                 double vt = 0;
03520                                 vt = points[upper].view == 2 ? vtx_u : vt;
03521                                 vt = points[upper].view == 3 ? vtx_v : vt;
03522                         
03523                                 double t =vt;
03524                                 
03525                                 //if the point at z is not the same as z vtx, we need to extrapolate 
03526                                 if(fabs ( points[upper].z - vz)>0.001)
03527                                 {                       
03528                         
03529                                         double s = (points[upper].t - vt) /  ( points[upper].z - vz);
03530                                         t = s * (points[i].z-vz) + vt;                  
03531                         
03532                                 }
03533                                 
03534                         //      printf("---view %d u %f v %f t %f\n",myview,u,v,t);
03535                         //      printf("---vtx z %f t%f upper z %f t %f\n",vz,vt,points[upper].z,points[upper].t);
03536                         
03537 
03538                                 u = myview == 2 ? u : t;
03539                                 v = myview == 3 ? v : t;
03540 //                      }       
03541                                 
03542 
03543                 }
03544                 else if(upper==-1 && lower==-1) //we have an empty view!!!
03545                 {
03546 
03547                         u = myview == 2 ? u : 0;
03548                         v = myview == 3 ? v : 0;
03549                 }
03550                 
03551                 p.add_to_back(u,v,z,e,points[i].chain,points[i].chainhit,view,rms_t);
03552                 
03553                 //printf("adding %f %f %f --- %f\n",u,v,z,e);
03554         }
03555         
03556         
03557         if(type)
03558                 p.particletype=(Particle3D::EParticle3DType)type;
03559         
03560         p.finalize();
03561         
03562         return p;
03563 }

void Finder::MakeChains ( Managed::ClusterManager cl,
ChainHelper ch,
int  view 
)

Definition at line 1235 of file Finder.cxx.

References ChainHelper::add_to_plane(), Managed::ManagedCluster::e, ParticleObjectHolder::event, ChainHelper::FindMaxPath(), ChainHelper::finish(), Managed::ClusterManager::GetCluster(), Managed::ClusterManager::GetClusterMap(), Managed::ManagedCluster::id, mypoh, ChainHelper::print_finished(), ChainHelper::process_plane(), ChainHelper::vtx_t, ParticleEvent::vtx_u, ParticleEvent::vtx_v, ChainHelper::vtx_z, and ParticleEvent::vtx_z.

01236 {
01237 
01238         if(mypoh->event.vtx_z>0)
01239         {
01240                 ch->vtx_z=mypoh->event.vtx_z;
01241                 ch->vtx_t=view==2?mypoh->event.vtx_u:mypoh->event.vtx_v;
01242                 
01243         }
01244 
01245 
01246      std::map<double, std::map<double, int> >::iterator p_iterr;
01247      std::map<double, int >::iterator s_iterr;
01248 
01249      double lastz=-1;
01250      for(p_iterr=cl->GetClusterMap(view)->begin();p_iterr!=cl->GetClusterMap(view)->end(); p_iterr++)
01251      {
01252         //for the first one...
01253                 if(lastz==-1)lastz=p_iterr->first;
01254                 
01255                 
01256                 if(fabs(p_iterr->first - lastz) > 0.01)
01257                 {
01258                         ch->process_plane();
01259                         lastz=p_iterr->first;
01260                         //printf("process plane\n");
01261                 }
01262       
01263 
01264          for(s_iterr=p_iterr->second.begin();s_iterr!=p_iterr->second.end(); s_iterr++)
01265          {
01266 
01267       
01268                 Managed::ManagedCluster * c = cl->GetCluster(s_iterr->second);
01269                 if(!c)
01270                 {
01271                         //printf("!!! bad cluster map! finder.cxx around 897\n");
01272                         continue;
01273                 }
01274                 //printf("t %f z %f e %f id %d\n",s_iterr->first, p_iterr->first, c->e,c->id);
01275                 ch->add_to_plane(s_iterr->first, p_iterr->first, c->e,c->id);   
01276              }
01277 // ch->process_plane();
01278         
01279      }
01280      
01281      ch->process_plane();
01282  
01283         ch->finish();
01284         ch->print_finished();
01285         ch->FindMaxPath();
01286 
01287 }

std::vector< std::pair< int, int > > Finder::matchViews ( std::vector< foundpath pu,
std::vector< foundpath pv 
)

Definition at line 1739 of file Finder.cxx.

References Msg::kDebug, Msg::kWarning, MSG, and Munits::second.

Referenced by Weave().

01740 {
01741         std::vector<std::pair<int,int> >matches;
01742         
01743 //      cout <<"views sizes "<<pu.size()<<" "<<pv.size()<<endl;
01744         for(unsigned int i=0;i<pu.size();i++)
01745         {
01746 
01747                 std::vector<std::pair<int,int> >zmatches;
01748 
01749 MSG("Finder",Msg::kDebug) <<"--- "<<i<<endl;
01750 
01751                 //if(pu[i].entries<2)continue;
01752         
01753                 std::vector<double> zoverlap;   
01754         
01755                 for(unsigned int j=0;j<pv.size();j++)
01756                 {
01757 
01758 MSG("Finder",Msg::kDebug) <<"----- "<<j<<endl;
01759 
01760                 //      if(pv[j].entries<2)continue;
01761                 
01762                         if((pu[i].entries==1 && pv[j].entries>2 )|| (pu[i].entries>2 && pv[j].entries==1 ))
01763                         {
01764                         
01765                                 zoverlap.push_back(0.0);
01766                                 continue; //don't match single hit with path longer than 2
01767                         }
01768                         
01769                         
01770         MSG("Finder",Msg::kDebug)<< "looking at "<<i << "  "<<j<<endl;
01771         MSG("Finder",Msg::kDebug) << " pu ent " << pu[i].entries << " pv ent " << pv[j].entries<<endl;
01772                 
01773                 
01774                         double o=0.0;
01775                         if(pu[i].start_z>pv[j].end_z || pu[i].end_z < pv[j].start_z) //no overlap!
01776                         {
01777                                 zoverlap.push_back(o);
01778                                 
01779                         MSG("Finder",Msg::kDebug) << "no overlap " << pu[i].start_z << " to " << pu[i].end_z << " vs " << pv[j].start_z << " to "<< pv[j].end_z <<endl;
01780                                 
01781                         //      continue;
01782                         }else{
01783                         
01784                         double start = pu[i].start_z < pv[j].start_z ? pv[j].start_z : pu[i].start_z;
01785                         double end = pu[i].end_z< pv[j].end_z ? pu[i].end_z : pv[j].end_z;
01786                         
01787                         o = fabs(end-start);
01788                         
01789                         double sm1= fabs(pu[i].end_z - pu[i].start_z);
01790                         double sm2= fabs(pv[j].end_z - pv[j].start_z);
01791                         double sm=0;
01792                         sm = sm1 > 0 ? sm1 : sm;
01793                         sm = sm2 > 0 && (sm2 < sm || sm ==0)? sm2 : sm;
01794                         
01795                         sm = sm1>sm2?sm1:sm2;
01796                         
01797                 MSG("Finder",Msg::kDebug)<< "sm1 "<<sm1<<" sm2 "<<sm2<<" pl " << o << " sm " << sm <<endl;
01798                         
01799                         if(sm<0.00001)continue;
01800                         
01801                         if(o==0) //special case --- we already made sure that there was overlap ... so this happens when one of the views has 1 hit, so its pathlength is 0
01802                         o=sm=1;
01803                         
01804                         zoverlap.push_back(o / sm);  // recording fraction of overlap over smaller path
01805                         
01806         MSG("Finder",Msg::kDebug)<< "overlap : " << o/sm<< " smallest path length "<<sm <<endl;
01807                         }
01808                 }
01809 
01810 
01811                 //find the best
01812                 double bestz=0.0;
01813                 for(unsigned int j=0;j<zoverlap.size();j++)
01814                         bestz = bestz > zoverlap[j]  ? bestz : zoverlap[j] > 0.999 ? bestz : zoverlap[j]; //if its a complete overlap, we don't want to compare against it...
01815                 
01816                         //require minimum overlap of 50%
01817                         if(bestz < 0.4 && bestz!=0.0) continue;
01818                 
01819                 
01820                 //find all matches within 30% of best
01821                 for(unsigned int j=0;j<zoverlap.size();j++)
01822                 {
01823                         if((bestz - zoverlap[j]) < 0.3 * bestz || zoverlap[j] > 0.999)
01824                         {
01825                                 zmatches.push_back(std::pair<int,int>(i,j));
01826                                 MSG("Finder",Msg::kDebug)<< "\tkeeping "<<i<<" "<<j<<endl;
01827                         }
01828                 }
01829                 
01830                 if(zmatches.size()==0)continue;
01831                 
01832 /*              if(zmatches.size()>1) //we need to pick a single match, so now look at energy and pick the closest!
01833                 {
01834                         double beste = 100000.0; // distance in energy from ith pu
01835                         int besteidx=-1;
01836                         for(unsigned int j=0;j<zmatches.size();j++)
01837                         {
01838                                 double de =  fabs(pu[zmatches[j].first].energy - pv[zmatches[j].second].energy);
01839                                 if(beste > de)
01840                                 {
01841                                         beste = de;
01842                                         besteidx=j;
01843                                 }
01844                         
01845                         }
01846                         
01847                         if(besteidx>-1)
01848                         matches.push_back(zmatches[besteidx]);
01849                         
01850                 }else{
01851                         matches.push_back(zmatches[0]);
01852                 }               
01853 */
01854 
01855 
01856 
01857                 if(zmatches.size()>1) //we need to pick a single match, so now look at energy / # entries and pick the closest!
01858                 {
01859                         int closest=100000;
01860                         int cidx=-1;
01861                         for(unsigned int j=0;j<zmatches.size();j++)
01862                         {
01863                         //      int d = fabs(pu[zmatches[j].first].energy / pu[zmatches[j].first].entries - pv[zmatches[j].first].energy / pv[zmatches[j].second].entries);
01864                                 
01865                         //      d=d/fabs(pu[zmatches[j].first].entries - pv[zmatches[j].second].entries);
01866                                 
01867                                 
01868                                 int d = (int)(TMath::Abs(pu[zmatches[j].first].entries - pv[zmatches[j].second].entries));
01869 
01870                                 if(d<closest)
01871                                 {
01872                                         closest=d;
01873                                         cidx=j;
01874                                 }
01875                         
01876                         }
01877                         
01878                         if(cidx>-1)
01879                         {
01880                                 matches.push_back(zmatches[cidx]);
01881                                 continue;                               
01882                         }
01883 
01884 
01885                 }
01886 
01887 
01888 
01889         
01890 
01891 
01892 
01893                 if (zmatches.size()>1) {
01894                         MSG("Finder",Msg::kWarning) << "Degeneracy not resolved in making 3d particle tracks!"<<endl;
01895                 }
01896                 for(unsigned int j=0;j<zmatches.size();j++)
01897                         matches.push_back(zmatches[j]);
01898 
01899         }
01900 
01901 
01902         return matches;
01903 }

void Finder::MergeChains ( ChainHelper ch  ) 

Definition at line 628 of file Finder.cxx.

References ChainHelper::AttachAsChild(), Chain::back_offset, Chain::back_slope, MuELoss::e, Chain::end_t, Chain::end_z, ChainHelper::finished, Chain::front_offset, Chain::front_slope, ChainHelper::GetChain(), VldContext::GetDetector(), RecRecordImp< T >::GetHeader(), RecHeader::GetVldContext(), Msg::kDebug, Detector::kFar, Detector::kNear, MSG, Chain::myId, mypoh, Munits::second, Chain::start_t, Chain::start_z, t, and z.

Referenced by Process().

00629 {
00630 
00631         std::vector<std::pair<int,int> >match;
00632 
00633         //iter over head chains
00634         for(int i=0;i<(int)ch->finished.size();i++)
00635         {
00636                 if(ch->finished[i].parentChain>-1)continue; //make it a head chain
00637                 if(ch->finished[i].entries<2)continue;
00638                 
00639         
00640                 double bestdist=10000;
00641                 //double maxtdist=0.1;
00642                 int bestid=-1;
00643                 double bestzdist=10000;
00644                 
00645                 //iter over ends of chains
00646                 for(int j=0;j<(int)ch->finished.size();j++)
00647                 {
00648                         if(i==j)continue;
00649                         if(ch->finished[j].children.size()>0)continue; //make it and ending chain
00650                         if(ch->finished[j].entries<2)continue;
00651                         
00652                         Chain *e  = &ch->finished[j]; //end of one chain
00653                         Chain *b = &ch->finished[i];  //beginning of another chain
00654                         
00655                         if(b->start_z-e->end_z < 0.035*2)continue; // require at least 2 planes away...
00656                 //      if(e->start_z - b->start_z< 0.035*2)continue; // require at least 2 planes away...
00657                         
00658                         double d = e->end_t - (e->end_z * b->front_slope + b->front_offset);
00659                 
00660                 //      double d = maxtdist;
00661                         
00662                         //need to make sure that the chain to be attached is vertex pointing...
00663                         
00664 
00665 /*
00667                         double canslope=b->front_slope;
00668                         double canoffset=b->front_offset;
00669                         double vpslope=e->back_slope;
00670                         double vpoffset=e->back_offset;
00671                         
00672                         
00673                         double t = ( vpoffset * canslope - vpslope *canoffset ) / (canslope - vpslope);
00674                         double z = (canoffset - vpoffset) / (vpslope - canslope);                       
00675                         
00676                         
00677                         printf("vp %d a %f b %f  can %d a %f b %f ------ t %f z %f\n",e->myId,e->a, e->b,b->myId,b->a, b->b,t,z);
00678                         printf("D t %f z%f \n",t,z);
00679                                 
00680                         if( z > e->end_z && z < b->start_z) // its a possible match! --- see if its the closest match...
00681                         if( t > e->end_t && t < b->start_t ||  t < e->end_t && t > b->start_t) //make sure its not kinky
00682                         {
00683                         
00684                                 printf("found match at z %f   t %f\n",z,t);
00685                                 
00686                                 d = e->end_t - (e->end_z * b->front_slope + b->front_offset);
00687                         
00688                         }
00690 */              
00691                         
00692                         
00693 //                      if(d<maxtdist)
00694 //                      {
00695                                 double totdist = sqrt( d*d + (e->end_z-b->start_z)*(e->end_z-b->start_z));
00696                         
00697                                 double zdist  = fabs(e->end_z-b->start_z);
00698                                 if(totdist<bestdist || zdist < bestzdist)
00699                                 {
00700                                         int points=0;
00701                                         double canslope=b->front_slope;
00702                                         double canoffset=b->front_offset;
00703                                         double vpslope=e->back_slope;
00704                                         double vpoffset=e->back_offset;
00705                                         double t = 0;
00706                                         double z = 0;
00707 
00708                                         if(canslope-vpslope>1e-10)
00709                                         {
00710                                                 t = ( vpoffset * canslope - vpslope *canoffset ) / (canslope - vpslope);
00711                                                 z = (canoffset - vpoffset) / (vpslope - canslope);
00712                                         }else{
00713                                                 continue;
00714                                         }
00715 
00716                                 //      if( z > e->end_z && z < b->start_z) points=1;
00717                                         
00718                                         if(fabs(e->back_slope - b->front_slope)<0.5)points=1;
00719                                 
00720                         //              cout <<"back slope "<< e->back_slope<<" front slope "<<b->front_slope<<endl;    
00721                                         
00722                                         //between supermodules in far?
00723                                         double maxtotdist=0.5;
00724                                         if(mypoh->GetHeader().GetVldContext().GetDetector() == Detector::kFar && fabs(e->end_z - 14.6)<0.2 && fabs(b->start_z-16)<0.2)maxtotdist+=16-14.6;
00725                                 
00726                                         //in back of near (partially instrumented every fourth plane)?
00727                                         if(mypoh->GetHeader().GetVldContext().GetDetector() == Detector::kNear && e->end_z >6 && b->start_z>6)maxtotdist*=4.;
00728                                 
00729                                 
00730                                         //make sure that the connection would not be kinky
00731                                         if(points)
00732                                         {
00733                                                 if((e->end_z-b->start_z)==0)continue;// that would be kinky anyways...
00734                                                 double connecting_slope = (e->end_t-b->start_t) / (e->end_z-b->start_z);
00735                                                 
00736                                                 if(fabs(connecting_slope - b->front_slope)>=0.5)points=0;
00737                                                 if(fabs(e->back_slope - connecting_slope)>=0.5)points=0;                                        
00738                                         }
00739                                 
00740                                         
00741                                         
00742                                         if(!points /*&& totdist>maxtotdist*/)
00743                                         {
00744                                                 MSG("Finder",Msg::kDebug) << "failing on points "<<e->myId <<" to "<< b->myId<<" totdist "<< totdist<<" zdist "<<zdist <<" tdist "<<d <<" intersection at z t"<<z<<" "<<t<<"\n";
00745                                                 continue;
00746                                         }
00747                                         
00748                                         
00749                                         
00750                                         bestdist=totdist;
00751                                         bestid=e->myId;
00752                                         bestzdist=zdist;
00753                                         
00754                                         MSG("Finder",Msg::kDebug) << "found better match "<<e->myId <<" to "<< b->myId<<" totdist "<< totdist<<" zdist "<<zdist <<" tdist "<<d <<"\n";
00755                                         
00756                                 }
00757                         }       
00758                         
00759 //              }
00760                 
00761                 
00762                 if(bestid>-1)
00763                 {
00764                         match.push_back(std::pair<int,int>(bestid,ch->finished[i].myId));
00765                 
00766                         MSG("Finder",Msg::kDebug) <<"!!!matched "<<bestid<<" "<<ch->finished[i].myId<<"\n";
00767                 
00768                 }
00769         }
00770         
00771         
00772         for(int i=0;i<(int)match.size();i++)
00773         {
00774                 ch->AttachAsChild(ch->GetChain(match[i].first),ch->GetChain(match[i].second));
00775         
00776         }
00777 
00778 
00779 
00780 }

void Finder::Process ( ParticleObjectHolder p  ) 

!!! do not remake clusters unless we are starting the chains from scratch!!!!

Definition at line 227 of file Finder.cxx.

References ParticleObjectHolder::AddParticle3D(), PrimaryShowerFinder::chu, ParticleObjectHolder::chu, ParticleObjectHolder::chv, PrimaryShowerFinder::chv, ParticleObjectHolder::cluster_map, ParticleObjectHolder::cluster_saver, ParticleObjectHolder::clusterer, clustermanager_v, Managed::ClusterManager::clusters, DoMRCC, Managed::ClusterManager::DumpClusters(), DumpParticle3D(), MuELoss::e, ParticleObjectHolder::event, ParticleObjectHolder::eventquality, exit(), Managed::ClusterManager::FillClusterMap(), FinalizeParticles3D(), PrimaryShowerFinder::FindPrimaryShower(), FindVertex(), EventQuality::foundlongmuon, EventQuality::foundprimaryshower, PrimaryShowerFinder::FoundSingleViewPrimaryShower(), ChainHelper::GetChain(), VldContext::GetDetector(), RecRecordImp< T >::GetHeader(), ParticleObjectHolder::GetMRCCInfo(), PrimaryShowerFinder::GetParticle3D(), Managed::ClusterSaver::GetStripEnergy(), RecHeader::GetVldContext(), it, Msg::kError, Msg::kInfo, Managed::ClusterManager::MakeClusters(), Managed::ClusterSaver::maxu, ParticleEvent::maxu, ParticleEvent::maxv, Managed::ClusterSaver::maxv, ParticleEvent::maxz, Managed::ClusterSaver::maxz, MergeChains(), meupergev, ParticleEvent::minu, Managed::ClusterSaver::minu, ParticleEvent::minv, Managed::ClusterSaver::minv, ParticleEvent::minz, Managed::ClusterSaver::minz, MSG, Particle3D::muon, mypoh, Managed::ClusterSaver::nClusters, ParticleEvent::nclusters, ChainHelper::NewChain(), ParticleEvent::nstrips, ParticleObjectHolder::particles3d, ChainHelper::print_finished(), ParticleObjectHolder::psf, Managed::ClusterSaver::recomputeBounds(), MRCCInfo::removedmuon, SetStatus(), PrimaryShowerFinder::SetVertex(), EventQuality::single_view_long_muon, EventQuality::single_view_primary_shower, MRCCInfo::stage, Particle3D::start_u, Particle3D::start_v, Particle3D::start_z, true_vtx_u, true_vtx_v, true_vtx_z, ParticleEvent::visenergy, vtx_u, vtx_v, vtx_z, and Weave().

Referenced by ParticleFinder::Reco().

00228 {
00229         if(!meupergev)
00230         {
00231                 cout <<"At Finder::Process... meupergev has not been set!\n";
00232                 exit(1);
00233         }
00234 
00235         
00236         mypoh=&p;
00237 
00238 //don't clear the xtalk here... use the XTalkFilter module!     
00239 //      if(p.GetHeader().GetVldContext().GetDetector()==Detector::kFar)ClearXTalk(); //only far has xtalk
00240         
00241 
00242 /*      
00243         for(unsigned int i=0;i<plane.size();i++)
00244         {
00245                 sorter_map.insert(std::pair <double,int>(energy[i],i));
00246                 loc_map[plane[i]][strip[i]]=i;
00247         
00248                 Managed::ClusterManager * clusterhelper = (view[i]==2) ? &mypoh->clusterer : (view[i]==3) ? &mypoh->clusterer : 0;
00249                 if(clusterhelper)clusterhelper->AddStrip(plane[i], strip[i], energy[i], t[i], z[i], view[i]);
00250 
00251         }
00252         
00253 */
00254 
00255         
00256 //mypoh->clusterer=clustermanager_u;
00257 mypoh->clusterer=clustermanager_v;
00258         
00259 //      mypoh->clusterer.MakeClusters();
00260 //      mypoh->clusterer.MakeClusters();
00261 
00262         mypoh->event.nclusters =  mypoh->clusterer.clusters.size();
00263 //      mypoh->event.nclusters += mypoh->clusterer.clusters.size();
00264         
00265 
00266         std::vector<Managed::ManagedCluster> clusters = mypoh->clusterer.clusters;
00267 /*      for(int i=0;i<clusters.size();i++)
00268         {
00269                 printf("clu id %d view %d z %f t %f e %f\n",clusters[i].id,clusters[i].view,clusters[i].z,clusters[i].t,clusters[i].e);
00270         
00271         
00272         }
00273 */
00274         //currently required for the display... should clean it up later!
00275         mypoh->clusterer.FillClusterMap(&mypoh->cluster_map);
00276         mypoh->clusterer.FillClusterMap(&mypoh->cluster_map);
00277         
00278 
00279         //look for long muons
00280         mypoh->clusterer.MakeClusters(0.02,0.05,0.01);
00281 
00282         LongMuonFinder lmf(mypoh->GetHeader().GetVldContext().GetDetector());
00283         int foundlongmuon = lmf.FindLongMuon(&mypoh->clusterer);
00284 
00285         int didMRCC=0;
00286 
00287         if(foundlongmuon )
00288         {
00289                         mypoh->eventquality.foundlongmuon=1;
00290                                 Particle3D *lp=lmf.GetParticle3D();
00291                                 
00292                                 if(lp)
00293                                 {
00294                                         std::vector<Particle3D*> pout1;
00295                                         pout1.push_back(lp);    
00296                 
00297                                         FinalizeParticles3D(pout1); //final computation of values, etc...
00298                         
00299                                         if(!DoMRCC)
00300                                         {
00301                                                 for(unsigned int i=0;i<pout1.size();i++)
00302                                                 {
00303                                                         mypoh->AddParticle3D(*(pout1[i]));
00304                                                 }
00305                                         }else{
00306                                                 MRCCInfo * mrccinfo = mypoh->GetMRCCInfo();
00307                                                 mrccinfo->removedmuon = *(pout1[0]); 
00308                                                 mrccinfo->stage=1;
00309                                                 didMRCC=1;
00310                                                 
00311                                         }               
00312                                 }
00313                                 
00314                                 //we need to store the chains to prevent further use...even if doing MRCC!
00315                                 
00316                                 int cu = mypoh->chu.NewChain();
00317                                 Chain * tempcu = mypoh->chu.GetChain(cu);
00318                                 *tempcu = *lmf.GetChain(2);
00319 
00320                                 int cv = mypoh->chv.NewChain();
00321                                 Chain * tempcv = mypoh->chv.GetChain(cv);
00322                                 *tempcv = *lmf.GetChain(3);
00323                                 
00324                                 if(DoMRCC)
00325                                 {
00326                                         SetStatus(tempcu,-10);
00327                                         SetStatus(tempcv,-10);
00328                                 }
00329                         
00330         }
00331 
00332 
00333         //cout <<"*********Currently have "<<mypoh->particles3d.size()<<" particles in POH\n";
00334         if(DoMRCC)foundlongmuon=0;//we don't want to use the muon vertex in the reco!
00335 
00336         if(lmf.FoundSingleViewLongMuon())
00337         {
00338                 mypoh->eventquality.single_view_long_muon=1;
00339                 //printf("!!!! Found a single good muon in only one view!\n");
00340         
00341         }
00342 
00343 
00344         //if we don't have a long muon, then this is probably either a NC or nue event
00345         //so lets look for a primary shower to get a fix on the event vertex
00346         mypoh->clusterer.MakeClusters(0.2,0.05,0.05);
00347 
00348 
00349         mypoh->clusterer.DumpClusters();
00350 
00351         int foundprimaryshower =0;
00352 
00353         PrimaryShowerFinder*psf = &mypoh->psf;  //mypoh has the psf for the display...
00354         psf->chu=&mypoh->chu;
00355         psf->chv=&mypoh->chv;
00356 
00357         //if foundlongmuon
00358         if(foundlongmuon)
00359         {
00360                 Particle3D *lp=lmf.GetParticle3D();
00361                 //printf("FOUND LONG MUON setting psf vertex %f %f %f\n",lp->start_u, lp->start_v, lp->start_z);
00362 
00363 
00364                 psf->SetVertex(lp->start_u, lp->start_v, lp->start_z);
00365         }
00366 
00367         foundprimaryshower = psf->FindPrimaryShower(&mypoh->clusterer);
00368 
00369 
00370         if(foundprimaryshower)
00371         {
00372                 mypoh->eventquality.foundprimaryshower=1;
00373         
00374                                 Particle3D *lp=psf->GetParticle3D();
00375                                 
00376                                 if(lp)
00377                                 {
00378                                         std::vector<Particle3D*> pout1;
00379                                         pout1.push_back(lp);    
00380                 
00381                                         FinalizeParticles3D(pout1); //final computation of values, etc...
00382 
00383                 
00384                                         if(!DoMRCC || (pout1.size()>0 && pout1[0]->particletype != Particle3D::muon) || didMRCC)
00385                                         {
00386                                                 for(unsigned int i=0;i<pout1.size();i++)
00387                                                 {
00388                                                         mypoh->AddParticle3D(*(pout1[i]));
00389                                                 }
00390                                         }else{
00391                                                 MRCCInfo * mrccinfo = mypoh->GetMRCCInfo();
00392                                                 mrccinfo->removedmuon = *(pout1[0]); 
00393                                                 mrccinfo->stage=2;
00394                                                 
00395                                         }               
00396 
00397                                 }
00398                 
00399                 //chains are being added in psf now...          
00400                 /*              
00401                                 int cu = mypoh->chu.NewChain();
00402                                 Chain * tempcu = mypoh->chu.GetChain(cu);
00403                                 *tempcu = *psf->GetChain(2);
00404 
00405                                 int cv = mypoh->chv.NewChain();
00406                                 Chain * tempcv = mypoh->chv.GetChain(cv);
00407                                 *tempcv = *psf->GetChain(3);
00408                 */              
00409                         
00410         }
00411 
00412 
00413         if(psf->FoundSingleViewPrimaryShower())
00414         {
00415                 mypoh->eventquality.single_view_primary_shower=1;
00416                 //printf("!!!! Found a single primary shower in only one view!\n");
00417         
00418         }
00419 
00420 
00421 
00422 
00423         //if foundprimaryshower
00424         if(foundprimaryshower)
00425         {
00426                 Particle3D *lp=psf->GetParticle3D();
00427                 vtx_u = lp->start_u;
00428                 vtx_v = lp->start_v;
00429                 vtx_z = lp->start_z;
00430         }       
00431         
00432         //if foundlongmuon
00433         if(foundlongmuon)
00434         {
00435                 Particle3D *lp=lmf.GetParticle3D();
00436                 
00437                 if(foundprimaryshower)
00438                 {
00439                         Particle3D *lpe=psf->GetParticle3D();
00440                         if(lpe->start_z > lp->start_z)
00441                         {
00442                                 vtx_u = lp->start_u;
00443                                 vtx_v = lp->start_v;
00444                                 vtx_z = lp->start_z;
00445                         }
00446                 }else{
00447                                 vtx_u = lp->start_u;
00448                                 vtx_v = lp->start_v;
00449                                 vtx_z = lp->start_z;            
00450                 
00451                 }
00452         }       
00453         
00454 
00455         //cout <<"*********Currently have "<<mypoh->particles3d.size()<<" particles in POH\n";
00456 
00457 //mypoh->clusterer.GetClusterSaver()->DumpClusters();
00458 
00459         //now look in clusters for tracks
00460         //remember to recluster!
00462 //      mypoh->clusterer.MakeClusters(0.2,0.05,0.05);
00463 
00464 
00465         MSG("Finder",Msg::kInfo) <<"----U----"<<endl;   
00466 
00467 
00468         
00469         ChainHelper * chu= &mypoh->chu;
00470         //Managed::ClusterManager  * clu = &mypoh->clusterer;
00471 /*      chu->max_z_gap=0.3;
00472         chu->max_slope=3;
00473         chu->SetDetector(mypoh->GetHeader().GetVldContext().GetDetector());
00474         MakeChains(clu,chu,2);
00475 */      
00476         MSG("Finder",Msg::kInfo) <<"----V----"<<endl;   
00477         
00478         ChainHelper * chv= &mypoh->chv;
00479         //Managed::ClusterManager  * clv = &mypoh->clusterer;
00480 /*      chv->max_z_gap=0.3;
00481         chv->max_slope=3;
00482         chv->SetDetector(mypoh->GetHeader().GetVldContext().GetDetector());
00483         MakeChains(clv,chv,3);
00484 */
00485 
00486  
00487         //first guess of vertex 
00488 //      FindVertex(chu,chv);
00489 //      MSG("Finder",Msg::kInfo) << "first guess - vertex at (u,v,z)  ("<< vtx_u <<", "<< vtx_v <<", "<< vtx_z <<")"<<endl;
00490 //      MSG("Finder",Msg::kInfo) << "first guess - true vertex at (u,v,z)  ("<<true_vtx_u <<", "<< true_vtx_v <<", "<< true_vtx_z <<")"<<endl;
00491 
00492         chu->print_finished();
00493         chv->print_finished();
00494 
00495 
00496 
00497         MergeChains(chu);
00498         MergeChains(chv);
00499 
00500         
00501 
00502 
00503         FindVertex(chu,chv);
00504         
00505         
00506         //if foundprimaryshower
00507         if(foundprimaryshower)
00508         {
00509                 Particle3D *lp=psf->GetParticle3D();
00510                 vtx_u = lp->start_u;
00511                 vtx_v = lp->start_v;
00512                 vtx_z = lp->start_z;
00513                 
00514         }       
00515         
00516         //if foundlongmuon
00517         if(foundlongmuon)
00518         {
00519                 Particle3D *lp=lmf.GetParticle3D();
00520                 
00521                 if(foundprimaryshower)
00522                 {
00523                         Particle3D *lpe=psf->GetParticle3D();
00524                         if(lpe->start_z > lp->start_z)
00525                         {
00526                                 vtx_u = lp->start_u;
00527                                 vtx_v = lp->start_v;
00528                                 vtx_z = lp->start_z;
00529                         }
00530                 }else{
00531                 
00532                         vtx_u = lp->start_u;
00533                         vtx_v = lp->start_v;
00534                         vtx_z = lp->start_z;            
00535                 }
00536         }       
00537                 
00538         
00539         MSG("Finder",Msg::kInfo) << "vertex at (u,v,z)  ("<< vtx_u <<", "<< vtx_v <<", "<< vtx_z <<")"<<endl;
00540         MSG("Finder",Msg::kInfo) << "true vertex at (u,v,z)  ("<<true_vtx_u <<", "<< true_vtx_v <<", "<< true_vtx_z <<")"<<endl;
00541 
00542 
00543 
00544 
00545 //      RemoveNonVertexPointingChains(chu,2);
00546 //      RemoveNonVertexPointingChains(chv,3);   
00547 
00548 
00549 
00550         chu->print_finished();
00551         chv->print_finished();
00552 
00553 //      mypoh->clusterer.GetClusterSaver()->DumpClusters();
00554 
00555         //cout <<"*********Currently have "<<mypoh->particles3d.size()<<" particles in POH\n";
00556         Weave(chu,chv);
00557         //cout <<"*********Currently have "<<mypoh->particles3d.size()<<" particles in POH\n";
00558 
00559 
00560         std::vector<Particle3D * > p3d;
00561 /*      for(int i=0;i<mypoh->particles3d1->GetEntries();i++)
00562                 p3d.push_back( (Particle3D * )mypoh->particles3d1->At(i) );
00563 */
00564 
00565         p3d.clear();
00566         for(unsigned int i=0;i<mypoh->particles3d.size();i++)
00567                 p3d.push_back(&mypoh->particles3d[i]);
00568         
00569         DumpParticle3D(p3d);
00570         //cout <<"*********Currently have "<<mypoh->particles3d.size()<<" particles in POH\n";
00571                                 
00572 /*      for (unsigned int i=0;i<p3d.size();i++)
00573         {
00574                 //ostringstream s;
00575         
00576                 Particle3D * p = p3d[i];
00577                 if (p==0)continue;
00578         
00579 //              printf("saving with parb %f\n",p->emfit_b);
00580         }
00581 */
00582 
00583 
00584         //and store event params now...
00585         mypoh->cluster_saver.recomputeBounds();
00586         mypoh->event.minu=mypoh->cluster_saver.minu;
00587         mypoh->event.maxu=mypoh->cluster_saver.maxu;
00588         mypoh->event.minv=mypoh->cluster_saver.minv;
00589         mypoh->event.maxv=mypoh->cluster_saver.maxv;
00590 
00591         mypoh->event.minz=mypoh->cluster_saver.minz;
00592         mypoh->event.maxz=mypoh->cluster_saver.maxz;
00593 
00594         //recalculate the number of strips and event energy
00595 
00596         std::map<std::pair<int,int>, double> stripmap =mypoh->cluster_saver.GetStripEnergy();
00597         std::map<std::pair<int,int>, double>::iterator it;
00598         
00599         int stripcount=0;
00600         double  stripe=0;
00601 
00602         for(it = stripmap.begin();it!=stripmap.end();it++)
00603         {
00604                 if(it->second<1e-6)continue;
00605                 stripcount++;
00606                 stripe+=it->second;
00607 //              printf("%d %d %f\n",it->first.first,it->first.second,it->second);
00608         }
00609 
00610         double olde = mypoh->event.visenergy;
00611         mypoh->event.nstrips=stripcount;
00612         mypoh->event.visenergy=stripe;
00613         mypoh->event.nclusters =  mypoh->cluster_saver.nClusters;
00614 
00615 
00616 //      printf("old e %f new e %f strips %d\n",olde,stripe,stripcount);
00617 
00618 
00619         if(stripe - olde > 1)
00620         {
00621                 MSG("Finder",Msg::kError) << "Event Energy MISMATCH! Strip SumE = "<<olde<<" new ClusterSaver Energy "<<stripe<<"\n";   
00622 
00623         }
00624 }

std::vector< Particle3D * > Finder::ProcessParticle3D ( Particle3D p3d  ) 

Definition at line 3698 of file Finder.cxx.

References MuELoss::a, Particle3D::e, Particle3D::electron, ParticleType::em, ParticleType::emshort, Particle3D::entries, Msg::kDebug, MSG, ParticleType::muon, Particle3D::muon, Particle3D::muonfrac, Particle3D::other, Particle3D::particletype, ParticleType::pi0, ParticleType::prot, Particle3D::shared, ParticleType::start, ParticleType::stop, StripMuon(), ParticleType::type, Particle3D::types, and ParticleType::uniquemuon.

03699 {
03700 
03701         std::vector<Particle3D*> pout;
03702 
03703         ostringstream s;
03704         s<<"\nparticle types:\n";
03705 
03706         //std::vector<ptype> types;
03707 
03708         //for the muon checks
03709         double lowm=0.2;
03710         double highm=3;
03711         int minm=2;
03712 
03713 
03714         //em check
03715         int lastemend=-1;
03716         int lastemlow=-1;
03717         if(p3d->entries>3)
03718         for(int i=0;i<p3d->entries-2;i++)
03719         {
03720                 double p0 = p3d->e[i];
03721                 double p1 = p3d->e[i+1];
03722                 double p2 = p3d->e[i+2];
03723                 if(p0<p1 && p1>p2 && p1 > highm) 
03724                 {
03725 
03726                         int low=i;
03727                         int high=i+2;
03728 
03729                         
03730                         for(int h=i-1;h>-1;h--)
03731                         {
03732                                 if(p3d->e[h]<p3d->e[h+1])
03733                                 {
03734                                         low=h;
03735                                 }else{
03736                                         break;
03737                                 }
03738                         }
03739                         for(int h=i+3;h<p3d->entries;h++)
03740                         {
03741                                 if(p3d->e[h]<p3d->e[h-1])
03742                                 {
03743                                         high=h;
03744                                 }else{
03745                                         break;
03746                                 }
03747                         }
03748                         if(low==lastemend)
03749                         {
03750 
03751                                 s<<"pi0 like - found consecutive em like from "<<lastemlow<<" to "<<high<<endl;
03752                                 ParticleType a;
03753                                 a.type = ParticleType::pi0;
03754                                 a.start=lastemlow;
03755                                 a.stop=high;
03756                                 p3d->types.push_back(a);
03757                         }
03758                         s<<"em like from "<< low << " to " <<high<<endl;
03759                         i=high; //start after this found em like structure
03760                         lastemend=high;
03761                         lastemlow=low;
03762                         ParticleType a;
03763                         a.type = ParticleType::em;
03764                         a.start=low;
03765                         a.stop=high;
03766                         p3d->types.push_back(a);        
03767                 }               
03768         }
03769         
03770         //shortem check
03771         if(p3d->entries>3)
03772         {
03773                 double p0 = p3d->e[0];
03774                 double p1 = p3d->e[1];
03775                 double p2 = p3d->e[2];  
03776                 if(p0>p1 && p1>p2 && p0 > highm)
03777                 {
03778 
03779                         s<<"shortem like"<<endl;
03780                         
03781                         int high=2;
03782                         for(int i=2;i<p3d->entries-1;i++)
03783                         {
03784                                 if(p3d->e[i] > p3d->e[i+1])
03785                                 {
03786                                         high=i+1;
03787                                 }       
03788                         }
03789                         ParticleType a;
03790                         a.type = ParticleType::emshort;
03791                         a.start=0;
03792                         a.stop=high;
03793                         p3d->types.push_back(a);        
03794                 }
03795         }       
03796         
03797         //muon check
03798         if(p3d->entries>1)
03799         for(int i=0;i<p3d->entries;i++)
03800         {
03801                 int low=i;
03802                 int high=i;
03803                 if(p3d->e[i] > lowm && p3d->e[i] < highm)
03804                 {
03805                         for(int j=i;j<p3d->entries;j++)
03806                         {
03807                                 if(p3d->e[j] > lowm && p3d->e[j] < highm)
03808                                 {
03809                                         high=j;
03810                                 }else{
03811                                         if(j < p3d->entries-1 && high == j-1 && p3d->e[j] < highm*2 && p3d->e[j+1] > lowm && p3d->e[j+1] < highm)  //allow 1 out of range hit
03812                                         {
03813                                                 high=j+1;
03814                                                 j++;
03815                                                 continue;
03816                                         }
03817                                         i=+1;
03818                                         break;
03819                                 
03820                                 }
03821                         }
03822                 }
03823                 
03824                 if(high-low+1>=minm)
03825                 {
03826 
03827                         s<<"muon like from "<<low<<" to "<<high<<endl;
03828                         ParticleType a;
03829                         a.type = ParticleType::muon;
03830                         a.start=low;
03831                         a.stop=high;
03832                         p3d->types.push_back(a);        
03833                 }
03834                 
03835                 i=high+1;
03836         }
03837 
03838         //unique muon check --- should we be extracting a muon from a particle that has larger shared hits?
03839         double  mc=0;
03840         double me=0;
03841         int c=0;
03842         for(int i=0;i<p3d->entries;i++)
03843         {
03844 
03845                 if(p3d->shared[i]==0)
03846                 {
03847                         c++;
03848                         if(p3d->e[i] > lowm && p3d->e[i] < highm)
03849                         {
03850                                 mc++;
03851                                 me+=p3d->e[i];
03852                         }
03853                 }
03854         }       
03855 
03856         if(mc/c>0.7)
03857         {
03858                 s<<"unique muon found"<<endl;
03859                 ParticleType a;
03860                 a.type = ParticleType::uniquemuon;
03861                 a.start=-1;
03862                 a.stop=-1;
03863                 p3d->types.push_back(a);        
03864         }
03866 
03867 
03868 
03869         
03870         //prot check // really stopping of +muon
03871         int plow=p3d->entries-1;
03872         int phigh=p3d->entries-1;
03873         if(p3d->e[p3d->entries-1]>highm)
03874         for(int i=p3d->entries-2;i>-1;i--)
03875         {
03876                 if(p3d->e[i]<p3d->e[i+1]) plow=i;
03877                 else break;
03878         }
03879         if(phigh-plow>1)
03880         {
03881 
03882                 s<<"proton/stopping u+ found from "<<plow<<" to "<<phigh<<endl;
03883                 ParticleType a;
03884                 a.type = ParticleType::prot;
03885                 a.start=plow;
03886                 a.stop=phigh;
03887                 p3d->types.push_back(a);        
03888         }
03889         
03890         //we now have a list of types.... go through, and split up the Particle3D as needed until all Particle3Ds have a single type
03891         
03892         
03893         int emt=0;
03894         
03895         p3d->particletype=Particle3D::other;    
03896         for(unsigned int i=0;i<p3d->types.size();i++)
03897         {
03898         
03899         
03900                 if(p3d->types[i].type==ParticleType::uniquemuon)
03901                 {
03902                         std::pair<Particle3D*,Particle3D*> a = StripMuon(p3d);
03903                         if(a.first)
03904                         {
03905                                 pout.push_back(a.first); //thats the stripped muon...
03906                                 if(a.second) //make sure we removed something from p3d...
03907                                 {
03908                                         std::vector<Particle3D*> b =ProcessParticle3D(a.second); //the adjusted remnant needs to be rechecked!
03909                                         for(unsigned int j=0;j<b.size();j++)
03910                                                 pout.push_back(b[j]);
03911                                 }
03912                                 return pout;
03913                         }       
03914                 }
03915                 
03916         
03917                 //see if it is definitely a muon!
03918                 
03919                 if(p3d->types[i].type==ParticleType::muon && p3d->types[i].stop-p3d->types[i].start+1 == p3d->entries)
03920                 {
03921                         p3d->particletype=Particle3D::muon;
03922                         break;
03923                 
03924                 }
03925 
03926                 
03927                 if(p3d->types[i].type==ParticleType::em && p3d->types[i].stop-p3d->types[i].start+1 == p3d->entries)
03928                 {
03929                         p3d->particletype=Particle3D::electron;
03930                         break;
03931                 }
03932 
03933                 if(p3d->types[i].type==ParticleType::em && (double)( p3d->types[i].stop-p3d->types[i].start) / (double) p3d->entries > 0.8) //more than half of it is em....
03934                 {
03935                         p3d->particletype=Particle3D::electron;
03936                         break;
03937                 }
03938 
03939                 if(p3d->types[i].type==ParticleType::em )
03940                 {
03941                         emt+= p3d->types[i].stop-p3d->types[i].start+1;
03942                 }
03943         
03944 
03945         
03946         }
03947 
03948         if(p3d->particletype==0 && p3d->muonfrac>0.6)
03949         {
03950                 p3d->particletype=Particle3D::muon;
03951                         
03952         }       
03953 
03954         if(p3d->particletype==0 && (double)emt / (double)p3d->entries > 40)
03955         {
03956                 p3d->particletype=Particle3D::electron;
03957                         
03958         }       
03959 
03960 
03961         MSG("Finder",Msg::kDebug) << s.str();   
03962         
03963         pout.push_back(p3d);
03964         return pout;
03965         
03966 }

std::vector< Particle3D * > Finder::ProcessParticle3D1 ( Particle3D p3d,
int  check_unique_muon = 1 
)

Definition at line 4055 of file Finder.cxx.

References MuELoss::a, Particle3D::e, Particle3D::electron, ParticleType::em, ParticleType::emshort, Particle3D::entries, it, Msg::kDebug, MSG, ParticleType::muon, Particle3D::muon, Particle3D::muonfrac, Particle3D::other, Particle3D::particletype, ParticleType::pi0, ParticleType::prot, ParticleType::start, ParticleType::stop, StripMuon1(), ParticleType::type, Particle3D::types, and ParticleType::uniquemuon.

Referenced by Weave().

04056 {
04057 
04058         std::vector<Particle3D*> pout;
04059 
04060         ostringstream s;
04061         s<<"\nparticle types:\n";
04062 
04063         //std::vector<ptype> types;
04064 
04065         //for the muon checks
04066         double lowm=0.2;
04067         double highm=3;
04068         int minm=2;
04069 
04070 
04071         //em check
04072         int lastemend=-1;
04073         int lastemlow=-1;
04074         if(p3d->entries>3)
04075         for(int i=0;i<p3d->entries-2;i++)
04076         {
04077                 double p0 = p3d->e[i];
04078                 double p1 = p3d->e[i+1];
04079                 double p2 = p3d->e[i+2];
04080                 if(p0<p1 && p1>p2 && p1 > highm) 
04081                 {
04082 
04083                         int low=i;
04084                         int high=i+2;
04085 
04086                         
04087                         for(int h=i-1;h>-1;h--)
04088                         {
04089                                 if(p3d->e[h]<p3d->e[h+1])
04090                                 {
04091                                         low=h;
04092                                 }else{
04093                                         break;
04094                                 }
04095                         }
04096                         for(int h=i+3;h<p3d->entries;h++)
04097                         {
04098                                 if(p3d->e[h]<p3d->e[h-1])
04099                                 {
04100                                         high=h;
04101                                 }else{
04102                                         break;
04103                                 }
04104                         }
04105                         if(low==lastemend)
04106                         {
04107 
04108                                 s<<"pi0 like - found consecutive em like from "<<lastemlow<<" to "<<high<<endl;
04109                                 ParticleType a;
04110                                 a.type = ParticleType::pi0;
04111                                 a.start=lastemlow;
04112                                 a.stop=high;
04113                                 p3d->types.push_back(a);
04114                         }
04115                         s<<"em like from "<< low << " to " <<high<<endl;
04116                         i=high; //start after this found em like structure
04117                         lastemend=high;
04118                         lastemlow=low;
04119                         ParticleType a;
04120                         a.type = ParticleType::em;
04121                         a.start=low;
04122                         a.stop=high;
04123                         p3d->types.push_back(a);        
04124                 }               
04125         }
04126         
04127         //shortem check
04128         if(p3d->entries>3)
04129         {
04130                 double p0 = p3d->e[0];
04131                 double p1 = p3d->e[1];
04132                 double p2 = p3d->e[2];  
04133                 if(p0>p1 && p1>p2 && p0 > highm)
04134                 {
04135 
04136                         s<<"shortem like"<<endl;
04137                         
04138                         int high=2;
04139                         for(int i=2;i<p3d->entries-1;i++)
04140                         {
04141                                 if(p3d->e[i] > p3d->e[i+1])
04142                                 {
04143                                         high=i+1;
04144                                 }       
04145                         }
04146                         ParticleType a;
04147                         a.type = ParticleType::emshort;
04148                         a.start=0;
04149                         a.stop=high;
04150                         p3d->types.push_back(a);        
04151                 }
04152         }       
04153         
04154         //muon check
04155         if(p3d->entries>1)
04156         for(int i=0;i<p3d->entries;i++)
04157         {
04158                 int low=i;
04159                 int high=i;
04160                 if(p3d->e[i] > lowm && p3d->e[i] < highm)
04161                 {
04162                         for(int j=i;j<p3d->entries;j++)
04163                         {
04164                                 if(p3d->e[j] > lowm && p3d->e[j] < highm)
04165                                 {
04166                                         high=j;
04167                                 }else{
04168                                         if(j < p3d->entries-1 && high == j-1 && p3d->e[j] < highm*2 && p3d->e[j+1] > lowm && p3d->e[j+1] < highm)  //allow 1 out of range hit
04169                                         {
04170                                                 high=j+1;
04171                                                 j++;
04172                                                 continue;
04173                                         }
04174                                         i=+1;
04175                                         break;
04176                                 
04177                                 }
04178                         }
04179                 }
04180                 
04181                 if(high-low+1>=minm)
04182                 {
04183 
04184                         s<<"muon like from "<<low<<" to "<<high<<endl;
04185                         ParticleType a;
04186                         a.type = ParticleType::muon;
04187                         a.start=low;
04188                         a.stop=high;
04189                         p3d->types.push_back(a);        
04190                 }
04191                 
04192                 i=high+1;
04193         }
04194 
04195         //unique muon check --- should we be extracting a muon from a particle that has larger shared hits?
04196 
04197         if(check_unique_muon)
04198         {
04199         double  mc=0;
04200         double me=0;
04201         int c=0;
04202         for(int i=0;i<p3d->entries;i++)
04203         {
04204 
04205 //              if(p3d->shared[i]==0)
04206 //              {
04207                         c++;
04208                         if(p3d->e[i] > lowm && p3d->e[i] < highm)
04209                         {
04210                                 mc++;
04211                                 me+=p3d->e[i];
04212                         }
04213 //              }
04214         }       
04215 
04216         if(mc>0.7*c)
04217         {
04218                 s<<"unique muon found muonlike to all hits ratio "<<mc/c<<endl;
04219                 ParticleType a;
04220                 a.type = ParticleType::uniquemuon;
04221                 a.start=-1;
04222                 a.stop=-1;
04223                 p3d->types.push_back(a);        
04224         }
04225         }
04227 
04228 
04229 
04230         
04231         //prot check // can also be a  stopping of +muon
04232         int plow=p3d->entries-1;
04233         int phigh=p3d->entries-1;
04234         if(p3d->e[p3d->entries-1]>highm)
04235         for(int i=p3d->entries-2;i>-1;i--)
04236         {
04237                 if(p3d->e[i]<p3d->e[i+1]) plow=i;
04238                 else break;
04239         }
04240         if(phigh-plow>0)
04241         {
04242 
04243                 s<<"proton/stopping u+ found from "<<plow<<" to "<<phigh<<endl;
04244                 ParticleType a;
04245                 a.type = ParticleType::prot;
04246                 a.start=plow;
04247                 a.stop=phigh;
04248                 p3d->types.push_back(a);        
04249         }
04250         
04251         //we now have a list of types.... go through, and split up the Particle3D as needed until all Particle3Ds have a single type
04252         
04253         
04254         int emt=0;
04255         
04256         //p3d->particletype=0;
04257         p3d->particletype=Particle3D::other;
04258         
04259         
04260         
04261         std::vector<ParticleType>::iterator it;
04262         
04263         
04264         //cout <<"num of types for this particle "<<  p3d->types.size()<<endl;
04265                 
04266         for(it=p3d->types.begin();it!=p3d->types.end();it++)
04267         {
04268         
04269                 
04270                 if(it->type==ParticleType::uniquemuon)
04271                 {
04272                         ostringstream s;
04273                         s<<"ATTEMPTING TO STRIP MUON\n";
04274                         std::pair<Particle3D*,Particle3D*> a = StripMuon1(p3d);
04275                         if(a.first)
04276                         {
04277                                 s << "STRIPPED MUON\n";
04278                                 pout.push_back(a.first); //thats the stripped muon...
04279                                 if(a.second) //make sure we removed something from p3d...
04280                                 {
04281                                         s <<"  OTHERS REMAIN!\n";
04282                                 
04283                                         std::vector<Particle3D*> b =ProcessParticle3D1(a.second); //the adjusted remnant needs to be rechecked!
04284                                         for(unsigned int j=0;j<b.size();j++)
04285                                                 pout.push_back(b[j]);
04286                                 }
04287                                 //return pout;
04288                                 break;
04289                         }else{  //remove the incorrectly assigned unique muon tag
04290                                 p3d->types.erase(it);
04291                         
04292                         }       
04293                 }
04294                 
04295                 MSG("Finder",Msg::kDebug)<<s;
04296         
04297         
04298                 //see if it is definitely a muon!
04299                 
04300                 if(it->type==ParticleType::muon && it->stop - it->start+1 == p3d->entries)
04301                 {
04302                         p3d->particletype=Particle3D::muon;
04303                 //      break;
04304                 
04305                 }
04306 
04307                 
04308                 if(it->type==ParticleType::em && it->stop - it->start+1 == p3d->entries)
04309                 {
04310                         p3d->particletype=Particle3D::electron;
04311                 //      break;
04312                 }
04313 
04314                 if(it->type==ParticleType::em && (double)( it->stop - it->start) / (double) p3d->entries > 0.8) //more than half of it is em....
04315                 {
04316                         p3d->particletype=Particle3D::electron;
04317                 //      break;
04318                 }
04319 
04320                 if(it->type==ParticleType::em )
04321                 {
04322                         emt+= it->stop - it->start+1;
04323                 }
04324         
04325 
04326         
04327         }
04328 
04329         if(p3d->particletype==0 && p3d->muonfrac>0.6)
04330         {
04331                 p3d->particletype=Particle3D::muon;
04332                         
04333         }       
04334 
04335         if(p3d->particletype==0 && (double)emt / (double)p3d->entries > 40)
04336         {
04337                 p3d->particletype=Particle3D::electron;
04338                         
04339         }       
04340 
04341 
04342         MSG("Finder",Msg::kDebug) << s.str()<<endl;     
04343         
04344         pout.push_back(p3d);
04345         return pout;
04346         
04347 }

void Finder::RecordLostHits ( std::vector< Particle3D * >  p3d  ) 

determine what hits are unused

Definition at line 2625 of file Finder.cxx.

References Particle3D::chain, Particle3D::chainhit, ParticleObjectHolder::chu, ParticleObjectHolder::chv, Chain::cluster_id, ParticleObjectHolder::clusterer, Managed::ClusterManager::clusters, ParticleObjectHolder::event, ChainHelper::GetChain(), Managed::ClusterManager::GetCluster(), Managed::ManagedCluster::inuse, mypoh, ParticleEvent::unused_e, ParticleEvent::unused_e_avg, ParticleEvent::unused_e_rms, ParticleEvent::unused_strips, and Particle3D::view.

02626 {
02627         for(int i=0;i<2;i++)
02628         {
02629                 Managed::ClusterManager * ch = i==0? &mypoh->clusterer : &mypoh->clusterer;
02630                 
02631                 for(unsigned int j=0;j<ch->clusters.size();j++)ch->clusters[j].inuse=0;
02632         }
02633         
02634         for(unsigned int i=0;i<p3d.size();i++)
02635         {
02636                 Particle3D * p = p3d[i];
02637                 for(unsigned int j=0;j<p->chain.size();j++)
02638                 {
02639                         Managed::ClusterManager * clh = p->view[j]==2 ? &mypoh->clusterer :  &mypoh->clusterer ;
02640                         ChainHelper * ch = p->view[j]==2 ? &mypoh->chu :  &mypoh->chv ;
02641                         Chain * chain = ch->GetChain(p->chain[j]);
02642                         if(p->chainhit[j]>-1)
02643                                 clh->GetCluster(chain->cluster_id[p->chainhit[j]])->inuse=1;
02644                 }
02645         }
02646         
02647         for(int i=0;i<2;i++)
02648         {
02649                 Managed::ClusterManager * ch = i==0? &mypoh->clusterer : &mypoh->clusterer;
02650                 
02651                 for(unsigned int j=0;j<ch->clusters.size();j++)
02652                 {
02653                         if(ch->clusters[j].inuse)continue;
02654                         for(unsigned int k=0;k<ch->clusters[j].hite.size();k++)
02655                         {
02656                                 mypoh->event.unused_e +=ch->clusters[j].hite[k];
02657                                 mypoh->event.unused_strips++;
02658                         
02659                         }
02660                 }
02661         }
02662 
02663         mypoh->event.unused_e_avg =0;
02664         if(mypoh->event.unused_strips>0)        
02665                 mypoh->event.unused_e_avg = mypoh->event.unused_e / mypoh->event.unused_strips;
02666         
02667         for(int i=0;i<2;i++)
02668         {
02669                 Managed::ClusterManager * ch = i==0? &mypoh->clusterer : &mypoh->clusterer;
02670                 
02671                 for(unsigned int j=0;j<ch->clusters.size();j++)
02672                 {
02673                         if(ch->clusters[j].inuse)continue;
02674                         for(unsigned int k=0;k<ch->clusters[j].hite.size();k++)
02675                         {
02676                                 mypoh->event.unused_e_rms += (ch->clusters[j].hite[k] - mypoh->event.unused_e_avg)*(ch->clusters[j].hite[k] - mypoh->event.unused_e_avg);
02677                         
02678                         }
02679                 }
02680         }
02681         
02682         mypoh->event.unused_e_rms=sqrt(mypoh->event.unused_e_rms);
02683         if(mypoh->event.unused_e_avg>0)mypoh->event.unused_e_rms/=mypoh->event.unused_e_avg;
02684         
02685 }

void Finder::RemoveNonVertexPointingChains ( ChainHelper ch,
int  view 
)

Definition at line 793 of file Finder.cxx.

References Chain::a, ChainHelper::AttachAsChild(), ChainHelper::DeleteChain(), Chain::end_t, Chain::end_z, Chain::entries, ChainHelper::finished, Chain::front_offset, Chain::front_slope, ChainHelper::GetChain(), Msg::kInfo, MSG, Chain::myId, ChainHelper::SplitAt(), Chain::start_t, Chain::start_z, t, vtx_u, vtx_v, vtx_z, and z.

00794 {
00795 
00796         double max_dt = 0.1; // max distance of chain pointing to vertex z from vertex in t
00797 
00798         double vt =0.0;
00799         if(view==2) vt=vtx_u;
00800         if(view==3) vt=vtx_v;
00801 
00802 
00803         std::vector<int> marked;
00804         std::vector<int> vtxpointing;
00805         std::vector<int> head;
00806 
00807 
00808 
00809         for(unsigned int i=0;i<ch->finished.size();i++)
00810         {       
00811                 ostringstream os;
00812                 
00813                 //here is where we can reverse chains if they point the wrong way
00814                 //chain begin should be closer to the vertex than the end....
00815                 
00816                 
00817                 
00818                 
00820 /*              
00821                 os << "Chain id " << ch->finished[i].myId << " entries: " << ch->finished[i].entries<< " extrap posvtx: " <<  ch->finished[i].avg_slope * vtx_z + ch->finished[i].avg_offset << " vtxt " << vt << "  diff " <<fabs(vt - (ch->finished[i].avg_slope * vtx_z + ch->finished[i].avg_offset) ) ;
00822         
00823                 if(ch->finished[i].entries>1)
00824                 if(ch->finished[i].parentChain<0 || (ch->finished[i].level==1 && ch->GetChain(ch->finished[i].parentChain)->entries==1) ) //only look at head chains...         
00825                 if(max_dt < fabs(vt - (ch->finished[i].avg_slope * vtx_z + ch->finished[i].avg_offset) ) )
00826                 {
00827                         //its too far away...
00828                         //throw it out
00829                         marked.push_back(ch->finished[i].myId);
00830                         
00831                         os << "not pointing to vertex ";
00832                         
00833                 }else{
00834                         vtxpointing.push_back(ch->finished[i].myId);
00835                 }
00836                 MSG("Finder",Msg::kInfo) << os.str() << endl; 
00837 */
00838 
00839 
00840                 int close_enough=0;
00841                 double diff = fabs(vt - (ch->finished[i].front_slope * vtx_z + ch->finished[i].front_offset) );
00842                 close_enough = max_dt > diff;
00843                 
00844 os << "Chain id " << ch->finished[i].myId << " entries: " << ch->finished[i].entries<< " extrap posvtx: " <<  ch->finished[i].front_slope * vtx_z + ch->finished[i].front_offset << " vtxt " << vt<<" diff "<<diff<<endl;
00845                                 
00846                 //using interpolation is probably just better for finding hits to attach...
00847                 diff = fabs(vt - ch->finished[i].interpolate(vtx_z) );
00848                 if(!close_enough)close_enough = max_dt > diff;
00849 
00850                 os << "Chain id " << ch->finished[i].myId << " entries: " << ch->finished[i].entries<< " extrap posvtx: " <<  ch->finished[i].interpolate(vtx_z)  << " vtxt " << vt<<" diff "<<diff<<endl;              
00851 
00852                 //using interpolation is probably just better for finding hits to attach...
00853                 diff = fabs(vt - (ch->finished[i].avg_slope * vtx_z + ch->finished[i].avg_offset) );
00854                 if(!close_enough)close_enough = max_dt > diff;
00855                 
00856                 
00857 
00858 
00859                 os << "Chain id " << ch->finished[i].myId << " entries: " << ch->finished[i].entries<< " extrap posvtx: " <<  (ch->finished[i].avg_slope * vtx_z + ch->finished[i].avg_offset)  << " vtxt " << vt <<" diff "<<diff;
00860 
00861                 /*        
00862                 if(ch->finished[i].entries>1)
00863                 if(ch->finished[i].parentChain<0 || (ch->finished[i].level==1 && ch->GetChain(ch->finished[i].parentChain)->entries==1) ) //only look at head chains...         
00864                 if(!close_enough)
00865                 {
00866                         //its too far away...
00867                         //throw it out
00868                         marked.push_back(ch->finished[i].myId);
00869                         
00870                         os << " not pointing to vertex ";
00871                         
00872                 }else{
00873                         vtxpointing.push_back(ch->finished[i].myId);
00874                 }
00875                 */
00876                 // RWH 2014-02-17 Above is what used to be here ...
00877                 // but the compiler is completely confused about dangling else 
00878                 // conditions and can't decide which to attach to which if
00879                 // this is current best guess of what (might have been)/was intended
00880                 if (ch->finished[i].entries>1) {
00881                   if (ch->finished[i].parentChain<0 || 
00882                       (ch->finished[i].level==1 && ch->GetChain(ch->finished[i].parentChain)->entries==1) ) {
00883                     //only look at head chains...         
00884                     if (!close_enough) {
00885                       //its too far away...
00886                       //throw it out
00887                       marked.push_back(ch->finished[i].myId);
00888                       
00889                       os << " not pointing to vertex ";
00890                       
00891                     } else {
00892                       vtxpointing.push_back(ch->finished[i].myId);
00893                     }
00894                   }
00895                 }
00896                 // RWH end of problematic area
00897 
00898                 MSG("Finder",Msg::kInfo) << os.str() << endl;
00899 
00900 
00901         }
00902         
00903         
00904         std::vector<int> todel;
00905         
00906         for (unsigned int i=0;i<marked.size();i++)
00907         {
00908                 //see if the chain is pointing to a vertexpointing chain....  if so, attach it...
00909                 //see if the projections cross somewhere between the start of the vertexpointing chain and the start of the candidate chain
00910                 
00911                 Chain * can = ch->GetChain(marked[i]);
00912                 
00913                 if(can->entries<2)continue;
00914                 
00915                 
00916                 double bestdist = 100000;
00917                 int bestpid = -1;
00918                 double attachpoint=0;
00919                 
00920                 for(unsigned int j=0;j<vtxpointing.size();j++)
00921                 {
00922                         //printf("looking to see if it points to a good chain for attachment....\n");
00923                 
00924                         Chain *vp = ch->GetChain(vtxpointing[j]);
00925                         //printf("A\n");        
00926                         if(vp->entries<2)continue;
00927                 
00928                         //printf("B\n");        
00929                         if(fabs(vp->a - can->a) < 0.000001)continue;  //dont want a fpe
00930                         //printf("C\n");        
00931                         if(can->start_z < vp->start_z)continue;
00932                         
00933         //              double t = ( vp->b * can->a - vp->a * can->b ) / (can->a - vp->a);
00934         //              double z = (can->b - vp->b) / (vp->a - can->a);
00935         
00936                 
00937 //                      double t = ( vp->a * can->b - vp->b * can->a ) / (can->b - vp->b);
00938 //                      double z = (can->a - vp->a) / (vp->b - can->b);
00939                 
00940 
00941                         double canslope=can->front_slope;
00942                         double canoffset=can->front_offset;
00943                         double vpslope=vp->front_slope;
00944                         double vpoffset=vp->front_offset;
00945                         
00946                         
00947                         double t = ( vpoffset * canslope - vpslope *canoffset ) / (canslope - vpslope);
00948                         double z = (canoffset - vpoffset) / (vpslope - canslope);                       
00949                         
00950                         
00951                         //printf("vp a %f b %f  can a %f b %f ------ t %f z %f\n",vp->a, vp->b,can->a, can->b,t,z);
00952                         
00953                         
00954                         
00955                 //      z+=vp->start_z;
00956                         //printf("D t %f z%f \n",t,z);  
00957                         if( z > vp->start_z && z < can->start_z) // its a possible match! --- see if its the closest match...
00958                         {
00959                         
00960                                 //printf("found match at z %f   t %f\n",z,t);
00961                                 
00962                                 double bd = bestdist+1  ;
00963                         
00964                                 attachpoint =z;
00965                         
00966                                 if( z < vp->end_z) // attach to middle of chain...
00967                                 {
00968                                         //printf("!!!! need to attach to middle of chain!!!\n");
00969                                         bd = sqrt ( (can->start_z-z)*(can->start_z-z) + (can->start_t - t)*(can->start_t - t));
00970                                         bestpid = vp->myId;
00971                                         
00972                                 }else{
00973                                         //printf("attach to end of chain\n");
00974                                         bd = sqrt ( (vp->end_z-can->start_z)*(vp->end_z-can->start_z) + (vp->end_t - can->start_t)*(vp->end_t - can->start_t));
00975                                 }
00976                                 if(bd<bestdist)
00977                                 {
00978                                         bestdist=bd;
00979                                         bestpid = vp->myId;
00980                                 }
00981                         }
00982                 }
00983                 
00984                 
00985         
00986         
00987                 //if not, delete it
00988                 if(bestpid>-1)
00989                 {
00990                         if(attachpoint < ch->GetChain(bestpid)->end_z) //need to break apart vertex pointing chain
00991                         {
00992                                 Chain * np = ch->SplitAt(ch->GetChain(bestpid),attachpoint);  //returns the chain that now ends at attachpoint
00993                 
00994                                 ch->AttachAsChild(np,can);
00995                 //              can->parentChain = np->myId;
00996                 //              np->children.push_back(can->myId);
00997                         
00998                         }else{   //just attach to the end
00999                         
01000                                 ch->AttachAsChild(ch->GetChain(bestpid),can);
01001                 //              can->parentChain = bestpid;
01002                 //              ch->GetChain(bestpid)->children.push_back(can->myId);
01003                         }
01004                 }else{
01005                         todel.push_back(marked[i]);
01006                 }
01007         }
01008         
01009 
01010 
01011         
01012         
01013         for(unsigned int i=0;i<todel.size();i++)
01014                 ch->DeleteChain(todel[i]);
01015                 
01016         
01017         
01018         
01019         ostringstream os;
01020         os << "remaing chains: ";
01021         for(unsigned int i=0;i<ch->finished.size();i++)
01022         {
01023                 os << ch->finished[i].myId<< " ";
01024         }
01025         
01026         MSG("Finder",Msg::kInfo) << os.str() << endl; 
01027 }

void Finder::RemoveSharedHits ( std::vector< Particle3D * >  pv  ) 

Definition at line 4472 of file Finder.cxx.

References ShareHolder::dump(), MuELoss::e, Particle3D::electron, ShareHolder::GetEShared(), ShareHolder::GetNumShared(), Msg::kDebug, MSG, Particle3D::muon, p3dsharedsort(), shareholder, and ShareHolder::Take().

04473 {
04474 
04475         MSG("Finder",Msg::kDebug) <<"\ncleaning list\n";
04476         for (unsigned int i=0;i<pv.size();i++)
04477         {
04478                         MSG("Finder",Msg::kDebug) << " removing with "<<pv[i]->numshared<<" shared \n";
04479         
04480                         //first clear out any hits marked as shared that are not!
04481                         for(int j=0;j<pv[i]->entries;j++)
04482                         {
04483                                 if(pv[i]->shared[j])
04484                                 {
04485                                         if(shareholder.GetNumShared(pv[i]->chain[j],pv[i]->chainhit[j])<2)
04486                                                 pv[i]->UnsetShared(pv[i]->chain[j],pv[i]->chainhit[j]);
04487                                         if(shareholder.GetNumShared(pv[i]->chain[j],pv[i]->chainhit[j])==1)
04488                                         {
04489                                                 int c=pv[i]->chain[j];
04490                                                 int ch=pv[i]->chainhit[j];
04491                                                 double e = shareholder.GetEShared(c,ch);
04492                                         
04493                                         
04494                                                 pv[i]->e[j]=e;
04495                                                 shareholder.Take(c,ch,e);
04496                                         }
04497                                 }
04498                         
04499                         }
04500         
04501         }
04502 
04503 
04504         sort(pv.begin(), pv.end(), p3dsharedsort);
04505 
04506         MSG("Finder",Msg::kDebug) <<"\nremoving hits\n";
04507 
04508         for (unsigned int i=0;i<pv.size();i++)
04509         {               
04510                         shareholder.dump();
04511 
04512                         //first clear out any hits marked as shared that are not!
04513                         for(int j=0;j<pv[i]->entries;j++)
04514                         {
04515                                 if(pv[i]->shared[j])
04516                                 {
04517                                         if(shareholder.GetNumShared(pv[i]->chain[j],pv[i]->chainhit[j])<2)
04518                                                 pv[i]->UnsetShared(pv[i]->chain[j],pv[i]->chainhit[j]); 
04519                                         if(shareholder.GetNumShared(pv[i]->chain[j],pv[i]->chainhit[j])==1)
04520                                         {
04521                                                 int c=pv[i]->chain[j];
04522                                                 int ch=pv[i]->chainhit[j];
04523                                                 double e = shareholder.GetEShared(c,ch);
04524                                         
04525                                         
04526                                                 pv[i]->e[j]=e;
04527                                                 shareholder.Take(c,ch,e);
04528                                         }               
04529                                 }
04530                         
04531                         }
04532 
04533 
04534 
04535 
04536 
04537                 if (pv[i]->particletype==Particle3D::muon)
04538                 {
04539                         MSG("Finder",Msg::kDebug)<<"!!muon\n";
04540                 }
04541                 else if (pv[i]->particletype==Particle3D::electron)
04542                 {
04543                         MSG("Finder",Msg::kDebug)<<"!!elec\n";
04544                 }
04545                 else MSG("Finder",Msg::kDebug)<<"!!other "<<pv[i]->particletype<<"\n";
04546                 
04547                 
04548                 //for now, treat all types similary.... 
04549                 //first look for shared hit surrounded by two unshared hits, and set it to the average
04550                 for(int j=1;j<pv[i]->entries-2;j++)
04551                 {
04552                         if(pv[i]->shared[j] && !pv[i]->shared[j-1] && !pv[i]->shared[j+1])
04553                         {
04554                                 int c=pv[i]->chain[j];
04555                                 int ch=pv[i]->chainhit[j];
04556                                 double e = shareholder.GetEShared(c,ch);
04557                                         
04558                                 double totake=(pv[i]->e[j-1]+pv[i]->e[j+1])/2;
04559                                 totake=totake>e?e:totake;
04560                                 
04561                                 pv[i]->e[j]=totake;
04562                                 shareholder.Take(c,ch,totake);
04563                                 pv[i]->UnsetShared(c,ch);
04564                         
04565                         }
04566                 }
04567                 
04568                 //if that doesn't work, try the average of closest unshared hits
04569                 for(int j=0;j<pv[i]->entries;j++)
04570                 {
04571                         int foundupper=-1;
04572                         int foundlower=-1;
04573                         for(int k=j+1;k<pv[i]->entries;k++)
04574                         {
04575                                 //printf("%d(%d) %d(%d)\n",j,pv[i]->shared[j],k,pv[i]->shared[k]);
04576                                 if(pv[i]->shared[j]&&! pv[i]->shared[k])
04577                                 {
04578                                         foundupper=k;
04579                                         break;
04580                                 }
04581                         }
04582                         
04583                         for(int k=j-1;k>-1;k--)
04584                         {
04585                                 if(pv[i]->shared[j]&&! pv[i]->shared[k])
04586                                 {
04587                                         foundlower=k;
04588                                         break;
04589                                 }
04590                         }
04591                         
04592                         double totake=0;
04593                         
04594                         if(foundupper>-1 && foundlower>-1)
04595                         {
04596                                 totake = (pv[i]->e[foundupper]+pv[i]->e[foundlower])/2;
04597                         }
04598                         else if(foundupper>-1)
04599                         {
04600                                 totake = (pv[i]->e[foundupper]);
04601                         }
04602                         else if(foundlower>-1)
04603                         {
04604                                 totake = (pv[i]->e[foundlower]);
04605                         }else continue;
04606                         
04607                         int c=pv[i]->chain[j];
04608                         int ch=pv[i]->chainhit[j];
04609                         double e = shareholder.GetEShared(c,ch);        
04610                         totake=totake>e?e:totake;
04611                         pv[i]->e[j]=totake;
04612                         MSG("Finder",Msg::kDebug)<<"taking from "<<c<<" "<<ch << " with shard "<<shareholder.GetNumShared(c,ch)<<"\n";
04613                         shareholder.Take(c,ch,totake);
04614                         MSG("Finder",Msg::kDebug)<<"now has "<<shareholder.GetNumShared(c,ch)<<" shared\n";
04615                         pv[i]->UnsetShared(c,ch);
04616                         
04617                 }
04618                 
04619                 pv[i]->finalize();
04620                 
04621         }
04622 
04623 
04624 }

void Finder::Reset (  ) 

Definition at line 58 of file Finder.cxx.

References energy, loc_map, particles, plane, ShareHolder::Reset(), shareholder, sorter_map, strip, t, true_vtx_u, true_vtx_v, true_vtx_z, view, vtx_u, vtx_v, vtx_z, and z.

Referenced by ParticleFinder::Reco().

00059 {
00060         sorter_map.clear(); 
00061 
00062 
00063         loc_map.clear(); 
00064         plane.clear(); 
00065         strip.clear(); 
00066         energy.clear(); 
00067         particles.clear(); 
00068         t.clear();
00069         z.clear(); 
00070         view.clear();  
00071         //cluster_map.clear();
00072 
00073         true_vtx_u=0.0;
00074         true_vtx_v=0.0;
00075         true_vtx_z=0.0;
00076 
00077 
00078         vtx_u=0.0;
00079         vtx_v=0.0;
00080         vtx_z=0.0;
00081         shareholder.Reset();
00082 }

void Finder::SetClusterManagerU ( Managed::ClusterManager m  )  [inline]

Definition at line 91 of file Finder.h.

References clustermanager_u.

Referenced by ParticleFinder::Reco().

00091 {clustermanager_u=m;};

void Finder::SetClusterManagerV ( Managed::ClusterManager m  )  [inline]

Definition at line 92 of file Finder.h.

References clustermanager_v.

Referenced by ParticleFinder::Reco().

00092 {clustermanager_v=m;};

void Finder::SetMEUperGeV ( double  d  )  [inline]

Definition at line 96 of file Finder.h.

References meupergev.

Referenced by ParticleFinder::Reco().

00096 {meupergev=d;};

void Finder::SetMRCC ( int  i  )  [inline]

Definition at line 94 of file Finder.h.

References DoMRCC.

Referenced by ParticleFinder::Reco().

00094 {DoMRCC=i;};

void Finder::SetPOH ( ParticleObjectHolder p  )  [inline]

Definition at line 39 of file Finder.h.

References mypoh.

Referenced by ParticleFinder::Reco().

00039 {mypoh = p;};

std::vector< Particle3D * > Finder::SetShared ( std::vector< Particle3D * >  p3v  ) 

Definition at line 3021 of file Finder.cxx.

References view.

03022 {
03023 
03024         //now iterate over found particles
03025         //and try to divide up shared energy!
03026         
03027         
03028         std::multimap< std::pair<int, std::pair<int,int> >, Particle3D* > pmap;
03029         for(int i=0;i<(int)p3v.size();i++)
03030         for(int j=0;j<p3v[i]->entries;j++)
03031         {
03032                 pmap.insert(make_pair( make_pair( p3v[i]->view[j],  make_pair(p3v[i]->chain[j],p3v[i]->chainhit[j])),p3v[i]) );
03033                 p3v[i]->UnsetShared(p3v[i]->chain[j],p3v[i]->chainhit[j]);
03034                 
03035         }
03036         
03037         
03038         std::pair<int,int>last;
03039         int lastcount=0;
03040         int lastview=0; 
03041         std::vector<Particle3D*> shared3ds;
03042         std::multimap< std::pair<int, std::pair<int,int> >, Particle3D* >::reverse_iterator it1;
03043         for(it1=pmap.rbegin();it1!=pmap.rend();it1++)
03044         {
03045         
03046                 if(it1->second->ShareLocked(it1->first.second.first, it1->first.second.second)>0)continue;
03047         
03048                 if(lastcount==0)
03049                 {
03050                         last=it1->first.second;
03051                         lastcount=1;
03052                         lastview=it1->first.first;
03053                         shared3ds.clear();
03054                         shared3ds.push_back(it1->second);
03055                         continue;       
03056                 }
03057                 
03058                 if(lastcount==1 && ( last!=it1->first.second || lastview!=it1->first.first))
03059                 {
03060                         last=it1->first.second;
03061                         lastview=it1->first.first;
03062                         shared3ds.clear();
03063                         shared3ds.push_back(it1->second);
03064                         continue;
03065                 }
03066                 
03067                 if(last==it1->first.second && lastview==it1->first.first )
03068                 {
03069                         lastcount++;
03070                         shared3ds.push_back(it1->second);
03071                         continue;
03072                 }else{
03073                         //do what we need to do with all of the last ones...
03074                         
03075                         //find all with this key...
03076                         
03077                         //then run a function over the list of the effected particle3d's noting the chain and chain id
03078                         //this function will look for adjacent hits in each of the particle3ds and attempt to extrapolate
03079                         //energy will be divided up evenly
03080                         //start at the back (highest z) and work down... as most shared hits should be closer to the vertex
03081                         
03082                         for(unsigned int i=0;i<shared3ds.size();i++)
03083                         {
03084                                 shared3ds[i]->SetShared(last.first, last.second);
03085                         }
03086                         
03087                         //printf("cmp %d %d shared %d times\n",last.first, last.second,lastcount);
03088                         
03089                         
03090                         
03091                         //then start with the new one..
03092                         lastcount=1;
03093                         last=it1->first.second;
03094                         lastview=it1->first.first;
03095                         shared3ds.clear();
03096                         shared3ds.push_back(it1->second);
03097                 }
03098                 
03099                 
03100         }
03101         
03102         if(lastcount>1)
03103         {
03104                 for(unsigned int i=0;i<shared3ds.size();i++)
03105                 {
03106                         shared3ds[i]->SetShared(last.first, last.second);
03107                 }
03108                 //printf("cmp %d %d shared %d times\n",last.first, last.second,lastcount);
03109         }
03110 
03111         return p3v;
03112 
03113 }

void Finder::SetStatus ( Chain c,
int  status 
) [private]

Definition at line 216 of file Finder.cxx.

References Chain::cluster_id, ParticleObjectHolder::clusterer, Managed::ClusterManager::GetCluster(), mypoh, and Managed::ManagedCluster::SetStatus().

Referenced by Process().

00217 {
00218 
00219         for(unsigned int i=0;i<c->cluster_id.size();i++)
00220         {
00221                 mypoh->clusterer.GetCluster(c->cluster_id[i])->SetStatus(status);
00222         //      printf("setting status cluster id %d to %d\n",c->cluster_id[i],status);
00223         }
00224 }

void Finder::SetTrueVertex ( double  myu,
double  myv,
double  myz 
) [inline]

Definition at line 52 of file Finder.h.

References true_vtx_u, true_vtx_v, and true_vtx_z.

Referenced by ParticleFinder::Reco().

00052 {true_vtx_u = myu;true_vtx_v = myv;true_vtx_z = myz;};

std::vector< Particle3D > Finder::shareEnergy ( std::vector< Particle3D p3v  ) 

Definition at line 2923 of file Finder.cxx.

References ShareHit(), and view.

02924 {
02925 
02926         //now iterate over found particles
02927         //and try to divide up shared energy!
02928         
02929         
02930         std::multimap< std::pair<int, std::pair<int,int> >, Particle3D* > pmap;
02931         for(int i=0;i<(int)p3v.size();i++)
02932         for(int j=0;j<p3v[i].entries;j++)
02933         {
02934                 pmap.insert(make_pair( make_pair( p3v[i].view[j],  make_pair(p3v[i].chain[j],p3v[i].chainhit[j])),&(p3v[i])) );
02935         
02936                 
02937         }
02938         
02939         
02940         std::pair<int,int>last;
02941         int lastcount=0;
02942         int lastview=0; 
02943         std::vector<Particle3D*> shared3ds;
02944         std::multimap< std::pair<int, std::pair<int,int> >, Particle3D* >::reverse_iterator it1;
02945         for(it1=pmap.rbegin();it1!=pmap.rend();it1++)
02946         {
02947                 if(lastcount==0)
02948                 {
02949                         last=it1->first.second;
02950                         lastcount=1;
02951                         lastview=it1->first.first;
02952                         shared3ds.clear();
02953                         shared3ds.push_back(it1->second);
02954                         continue;       
02955                 }
02956                 
02957                 if(lastcount==1 && ( last!=it1->first.second || lastview!=it1->first.first))
02958                 {
02959                         last=it1->first.second;
02960                         lastview=it1->first.first;
02961                         shared3ds.clear();
02962                         shared3ds.push_back(it1->second);
02963                         continue;
02964                 }
02965                 
02966                 if(last==it1->first.second && lastview==it1->first.first)
02967                 {
02968                         lastcount++;
02969                         shared3ds.push_back(it1->second);
02970                         continue;
02971                 }else{
02972                         //do what we need to do with all of the last ones...
02973                         
02974                         //find all with this key...
02975                         
02976                         //then run a function over the list of the effected particle3d's noting the chain and chain id
02977                         //this function will look for adjacent hits in each of the particle3ds and attempt to extrapolate
02978                         //energy will be divided up evenly
02979                         //start at the back (highest z) and work down... as most shared hits should be closer to the vertex
02980                         
02981                         for(unsigned int i=0;i<shared3ds.size();i++)
02982                         {
02983                                 shared3ds[i]->SetShared(last.first, last.second);
02984                         }
02985                         
02986                         ShareHit(lastview,last.first,last.second,shared3ds);
02987                         
02988                         
02989                         
02990                         
02991                         //printf("cmp %d %d shared %d times\n",last.first, last.second,lastcount);
02992                         
02993                         
02994                         
02995                         //then start with the new one..
02996                         lastcount=1;
02997                         last=it1->first.second;
02998                         lastview=it1->first.first;
02999                         shared3ds.clear();
03000                         shared3ds.push_back(it1->second);
03001                 }
03002                 
03003                 
03004         }
03005         
03006         if(lastcount>1)
03007         {
03008                 for(unsigned int i=0;i<shared3ds.size();i++)
03009                 {
03010                         shared3ds[i]->SetShared(last.first, last.second);
03011                 }
03012                 ShareHit(lastview,last.first,last.second,shared3ds);
03013                 //printf("cmp %d %d shared %d times\n",last.first, last.second,lastcount);
03014         }
03015 
03016         return p3v;
03017 
03018 }

void Finder::ShareHit ( int  view,
int  chain,
int  chainhit,
std::vector< Particle3D * >  shared 
)

Definition at line 3116 of file Finder.cxx.

References ParticleObjectHolder::chu, ParticleObjectHolder::chv, Chain::e, energy, ChainHelper::GetChain(), mypoh, n, and total().

Referenced by shareEnergy().

03117 {
03118 
03119         //for now, divide the energy based on the average of the next/previous energies
03120         //in a chain as a fraction of the total n/p energy in all of the shared chains....
03121         
03122         
03123         std::vector<double> es;
03124         for(unsigned int i=0;i<shared.size();i++)
03125         {
03126                 int inview = view == 2 ? 3 :2; //look in the opposite view
03127                 
03128                 double n = shared[i]->GetNextEnergy(chain,chainhit,inview);
03129                 double p = shared[i]->GetPreviousEnergy(chain,chainhit,inview);
03130                 //if either is 0 then we are the end of the chain, so count the non zero one twice
03131                 if(n==0 && p >0)n=p;
03132                 if(p==0 && n>0)p=n;
03133                 
03134                 //printf("e %f\n",n+p);
03135                 
03136                 es.push_back(n+p);
03137         }
03138         
03139         double total=0;;
03140         double totaln=0;
03141         for(unsigned int i=0;i<es.size();i++)
03142         {
03143                 totaln++;
03144                 total+=es[i];
03145         }
03146         
03147         ChainHelper * ch = view == 2 ? &mypoh->chu : &mypoh->chv;
03148         
03149         Chain *c =ch->GetChain(chain);
03150         double energy = c->e[chainhit];
03151         
03152         //printf("Sharing %f\n",energy);
03153         
03154         if(total==0)return; //shared hits with no energy!
03155         for(unsigned int i=0;i<es.size();i++)
03156         {
03157                 shared[i]->SetEnergy(chain,chainhit,view,energy*es[i]/total);
03158         }
03159         
03160 
03161 }

std::pair< Particle3D *, Particle3D * > Finder::StripMuon ( Particle3D p3d  ) 

Definition at line 3969 of file Finder.cxx.

References MuELoss::a, Particle3D::add_to_back(), Particle3D::chain, Particle3D::chainhit, MuELoss::e, Particle3D::e, Particle3D::entries, Particle3D::lockshared, Munits::m, Particle3D::muon, Particle3D::particletype, Particle3D::shared, Particle3D::u, Particle3D::v, Particle3D::view, and Particle3D::z.

Referenced by ProcessParticle3D().

03970 {
03971         std::pair<Particle3D*,Particle3D*>  a;
03972         a.first=0;
03973         a.second=p3d;
03974         
03975         //see if all unique hits are muons hits....
03976         double lowm=0.0;
03977         double highm=3.5;
03978         
03979         int u=0;
03980         int m=0;
03981         double me=0;
03982         
03983         for(int i=0;i<p3d->entries;i++)
03984         {
03985                 if(p3d->shared[i]==0)
03986                 {
03987                         u++;
03988                         if(p3d->e[i] > lowm && p3d->e[i]<highm)
03989                         {
03990                                 m++;
03991                                 me+=p3d->e[i];
03992                         }
03993                 }
03994         }
03995         
03996         if(u==m)
03997         { 
03998                 a.second=0;  //releasing shared particles which cant be part of something else...
03999                 //cout <<"!!!---releasing particle"<<endl;
04000         }
04001         
04002         
04003         Particle3D * c = new Particle3D();
04004         a.first=c;
04005         
04006         
04007         double eavg = me/m;
04008         for(int i=0;i<p3d->entries;i++)
04009         {
04010                 double e = p3d->e[i] > lowm && p3d->e[i]<highm ? p3d->e[i] : eavg;
04011                 e = e > p3d->e[i] ? p3d->e[i] : e;
04012         
04013                 c->add_to_back( p3d->u[i],  p3d->v[i],  p3d->z[i], e, p3d->chain[i], p3d->chainhit[i],  p3d->view[i],0);
04014                 c->lockshared[c->entries-1]=1;
04015                 
04016                 if(a.second)
04017                 {
04018                         double reste=p3d->e[i]-e>0 ? p3d->e[i]-e : 0;
04019                         a.second->e[i]= reste;
04020                         
04021                         
04022                 }
04023         }
04024         c->particletype=Particle3D::muon;
04025 
04026 
04027 
04028         if(a.second)
04029         {
04030                 double sume=0;
04031                 int num0=0;
04032                 for(int i=0;i<a.second->entries;i++)
04033                 {
04034                         sume+=a.second->e[i];
04035                         if(a.second->e[i]<0.01)num0++;
04036                 }
04037                 //cout<<"second e"<<sume<<endl;
04038                 
04039                 if(sume<0.0001 || (double)num0/(double)a.second->entries > 0.8)
04040                 {
04041                         a.second=0;  //releasing shared particles which cant be part of something else...
04042                         //cout <<"!!!---releasing particle"<<endl;
04043                 }
04044         }
04045 
04046         return a;
04047 
04048 }

std::pair< Particle3D *, Particle3D * > Finder::StripMuon1 ( Particle3D p3d  ) 

Definition at line 4351 of file Finder.cxx.

References MuELoss::a, Particle3D::add_to_back(), Particle3D::chain, Particle3D::chainhit, MuELoss::e, Particle3D::e, Particle3D::entries, it, Msg::kDebug, Particle3D::lockshared, Munits::m, MSG, ParticleType::muon, Particle3D::muon, Particle3D::particletype, Particle3D::shared, shareholder, ParticleType::start, ParticleType::stop, ShareHolder::Take(), ParticleType::type, Particle3D::types, Particle3D::u, ParticleType::uniquemuon, Particle3D::v, Particle3D::view, and Particle3D::z.

Referenced by ProcessParticle3D1().

04352 {
04353         std::pair<Particle3D*,Particle3D*>  a;
04354         a.first=0;
04355         a.second=p3d;
04356         
04357         //see if all unique hits are muons hits....
04358         double lowm=0.0;
04359         double highm=3.5;
04360         
04361         int u=0;
04362         int m=0;
04363         double me=0;
04364         
04365         for(int i=0;i<p3d->entries;i++)
04366         {
04367                 if(p3d->shared[i]==0)
04368                 {
04369                         u++;
04370                         if(p3d->e[i] > lowm && p3d->e[i]<highm)
04371                         {
04372                                 m++;
04373                                 me+=p3d->e[i];
04374                         }
04375                 }
04376         }
04377         
04378         if(u==m)
04379         { 
04380                 a.second=0;  //releasing shared particles which cant be part of something else...
04381                 //cout <<"!!!---releasing particle"<<endl;
04382         }
04383         
04384         
04385         Particle3D * c = new Particle3D();
04386         a.first=c;
04387         
04388         
04389         double eavg = me/m;
04390         for(int i=0;i<p3d->entries;i++)
04391         {
04392                 double e = p3d->e[i] > lowm && p3d->e[i]<highm ? p3d->e[i] : eavg;
04393                 e = e > p3d->e[i] ? p3d->e[i] : e;
04394         
04395                 c->add_to_back( p3d->u[i],  p3d->v[i],  p3d->z[i], e, p3d->chain[i], p3d->chainhit[i],  p3d->view[i],0);
04396                 c->lockshared[c->entries-1]=1;
04397 
04398                 shareholder.Take(p3d->chain[i],p3d->chainhit[i],e);
04399                 
04400                 if(a.second)
04401                 {
04402                         double reste=p3d->e[i]-e>0 ? p3d->e[i]-e : 0;
04403                         a.second->e[i]= reste;
04404                         
04405                         shareholder.Take(p3d->chain[i],p3d->chainhit[i],e);
04406                         
04407                         
04408                 }
04409         }
04410         c->particletype=Particle3D::muon;
04411 
04412         ParticleType pta;
04413         pta.type = ParticleType::muon;
04414         pta.start=0;
04415         pta.stop=p3d->entries;
04416         c->types.push_back(pta);        
04417 
04418 
04419         if(a.second)
04420         {
04421                 double sume=0;
04422                 int num0=0;
04423                 for(int i=0;i<a.second->entries;i++)
04424                 {
04425                         sume+=a.second->e[i];
04426                         if(a.second->e[i]<0.01)num0++;
04427                 }
04428                 MSG("Finder",Msg::kDebug)<<"second e"<<sume<<endl;
04429                 
04430                 if(sume<0.0001 || (double)num0/(double)a.second->entries > 0.8)
04431                 {
04432                         a.second=0;  //releasing shared particles which cant be part of something else...
04433                         //cout <<"!!!---releasing particle"<<endl;
04434                 }
04435                 else  //we are keeping the particle....
04436                 {
04437 
04438 
04439                         //remove the unique muon indicator....
04440                 
04441                         std::vector<ParticleType>::iterator it;
04442                         for(it=a.second->types.begin();it!=a.second->types.end();it++)
04443                         {
04444                                 if(it->type==ParticleType::uniquemuon)
04445                                 {
04446                                         a.second->types.erase(it);
04447                                         it=a.second->types.begin(); //for now restart at the beginning... not the best way
04448                                         MSG("Finder",Msg::kDebug)<<"removed uniquemuon indicator from particle3d "<<endl;
04449                                 }
04450                         
04451                         }
04452 
04453 
04454                         //clear particle list....
04455                         a.second->types.clear();
04456 
04457                 }
04458         }
04459 
04460 
04461         return a;
04462 
04463 }

void Finder::Weave ( ChainHelper chu,
ChainHelper chv 
)

Definition at line 1906 of file Finder.cxx.

References ParticleObjectHolder::AddParticle3D(), Particle3D::chain, Particle3D::chainhit, ChainHelper::ChangeHitId(), Chain::cluster_id, ParticleObjectHolder::clusterer, DumpParticle3D(), DumpPaths(), Managed::ManagedCluster::e, MuELoss::e, Particle3D::e, Particle3D::entries, ParticleObjectHolder::event, FinalizeParticles3D(), FindNeutrons(), ChainHelper::GetChain(), Managed::ClusterManager::GetCluster(), GetPaths(), Msg::kDebug, Msg::kError, Make3DParticle(), matchViews(), MSG, mypoh, ProcessParticle3D1(), Managed::ClusterManager::SaveCluster(), Munits::second, Particle3D::sum_e, Particle3D::view, and ParticleEvent::visenergy.

Referenced by Process().

01907 {
01908         std::vector<foundpath> pu = GetPaths(chu);
01909         std::vector<foundpath> pv = GetPaths(chv);
01910 
01911         if(pu.size()<1 || pv.size()<1)return;
01912 
01913         
01914         DumpPaths(pu);
01915         DumpPaths(pv);
01916 
01918         //we should clean each found path vector to remove paths which differ only by the last chain
01919         //as happens when a muon brems
01920         
01921         for(int i=0;i<2;i++){
01922         
01923         std::vector<foundpath>::iterator it1;
01924         std::vector<foundpath>::iterator it2;   
01925         std::vector<foundpath> paths=i==0?pu:pv;
01926         if(paths.size()<2)continue;
01927         
01928         double minzlength=1.0; //distance in m before we can delete...
01929         
01930         it1 =paths.begin();
01931         while(it1!=paths.end())
01932         //for(it1 =paths.begin();it1!=paths.end();it1++)
01933         {
01934                 if(it1->path.size()<2)
01935                 {
01936                         it1++;
01937                         continue;
01938                 }
01939                 int diddel=0;
01940                 
01941                 it2=it1;
01942                 it2++;
01943                 while(it2!=paths.end())
01944                 //for(it2 =it1;it2!=paths.end();it2++)
01945                 {
01946                         diddel=0;
01947                         if(it2->path.size()<2)
01948                         {
01949                                 it2++;
01950                                 continue;
01951                         }
01952 
01953                         
01954                         //MSG("Finder",Msg::kDebug)<<"paths start / end /it1 /it2 "<<&(paths.begin())<<" "<<&(paths.end())<<" "<<&it1<<" "<<&it2<<endl;
01955                         
01956         
01957                         
01958                         //MSG("Finder",Msg::kDebug)<<"Comparing with sizes "<<it1->path.size()<<" "<<it2->path.size()<<endl;
01959                         
01960                         int close=1;
01961                         
01962                 //      int larger=it1->path.size()>it2->path.size()?0:1;
01963                         int larger=it1->end_z-it1->start_z > it2->end_z-it2->start_z ? 0:1;
01964                 
01965                 
01966                         std::vector<int> maxpi = larger==0?it1->path:it2->path;
01967                         std::vector<int> maxpj = larger==1?it1->path:it2->path; 
01968                         
01969                         for(unsigned int k=0;k<maxpj.size()-1;k++)
01970                         {
01971                         MSG("Finder",Msg::kDebug)<<"chains "<<maxpi[k]<<" "<<maxpj[k]<<endl;
01972                                 if(maxpi[k]!=maxpj[k])
01973                                 {
01974                                         close=0;
01975                                         break;
01976                                 }
01977                         }
01978                         
01979                         //delete the smaller one?
01980                         std::vector<foundpath>::iterator itsmall=larger==0?it2:it1;
01981                         double dist = itsmall->end_z-itsmall->start_z;
01982                         if(close && dist > minzlength)
01983                         {
01984                         /*MSG("Finder",Msg::kError)*/
01985                         /*cout<<"removing close chain (";
01986                         for(unsigned int q=0;q<it1->path.size();q++)cout <<it1->path[q]<<" ";
01987                         cout<<") (";
01988                         for(unsigned int q=0;q<it2->path.size();q++)cout <<it2->path[q]<<" ";
01989                         cout<<") \n";                   
01990                         */
01991 
01992                                 int both=it1==it2?1:0;
01993                                 if(larger==1)
01994                                 {
01995                                         it1=paths.erase(it1);
01996                                         if(both)it2=it1;
01997                                 }
01998                                 if(larger==0)
01999                                 {
02000                                         it2=paths.erase(it2);
02001                                         if(both)it2=it1;
02002                                 }
02003                                 
02004                                 //something doesnt work right....!
02005                                 //do this sucky hack
02006                                 diddel=1;
02007                                 
02008                                                 //MSG("Finder",Msg::kError)<<"current pointers start / end /it1 /it2 "<<&(paths.begin())<<" "<<&(paths.end())<<" "<<&it1<<" "<<&it2<<endl;
02009                         }
02010                         
02011                         if(diddel)
02012                         {
02013                                 it1=paths.begin();
02014                                 it2=it1;
02015                         }
02016                         it2++;
02017                         
02018                 }       
02019                 it1++;
02020                 
02021         }
02022                 if(i==0)pu=paths;
02023                 if(i==1)pv=paths;
02024         }//loop over views
02026 
02027         if(pu.size()<1 || pv.size()<1)
02028         {
02029                 MSG("Finder",Msg::kError)<<"Removal of similar paths results in empty view!\n"; 
02030                 return;
02031         }
02032 
02033 
02034         //cout <<"*********Currently have "<<mypoh->particles3d.size()<<" particles in POH\n";
02035 
02036         DumpPaths(pu);
02037         DumpPaths(pv);
02038         
02039         
02040         //now try to find best matches!
02041         std::vector<std::pair<int,int> >matches = matchViews(pu,pv);    
02042         for(unsigned int i=0;i<matches.size();i++)
02043         {
02044                 MSG("Finder",Msg::kDebug) << "found pair " << matches[i].first << " " <<matches[i].second<<endl;
02045         }
02046 
02047         //and check the other view first as well!
02048         std::vector<std::pair<int,int> >matcheso = matchViews(pv,pu);
02049         for(unsigned int i=0;i<matcheso.size();i++)
02050         {
02051                 MSG("Finder",Msg::kDebug) << "found pair " << matcheso[i].second << " " <<matcheso[i].first<<endl;
02052         }
02053         for(unsigned int i=0;i<matcheso.size();i++)
02054         {
02055                 int found =0;
02056                 std::pair<int,int> tmp = make_pair(matcheso[i].second, matcheso[i].first); //need to flip the pair
02057                 for(unsigned int j=0;j<matches.size();j++)
02058                 {
02059                         
02060                         if(matches[j]==tmp)
02061                         {
02062                                 found=1;
02063                                 break;
02064                         }
02065                 }
02066                 if(!found)matches.push_back(tmp);
02067         
02068         }
02069         
02070 
02071         
02072 //      printf("!!1 matches size %d\n",matches.size());
02073         
02074         
02075         if(matches.size()<1)return;
02076         
02077         ostringstream os;
02078         os << "Candidate matches on z overlap : ";
02079         for(unsigned int i=0;i<matches.size();i++)
02080         {
02081                 os << matches[i].first <<"-"<<matches[i].second<<" ";
02082         }
02083         MSG("Finder",Msg::kDebug) << os.str()<<endl;    
02084 
02085 
02086 
02087         //when getting rid of a match... make sure its the best match!  (as we can share a chain in one view!)
02088         //need to see if all matches are "fair"
02089         std::vector<std::pair<int,int> >matchestmp;
02090         std::vector<int>multu;
02091         std::vector<int>multv;
02092         for(unsigned int i=0;i<matches.size();i++)
02093         {
02094                 int mult0=0;
02095                 int mult1=0;
02096                 for(unsigned int j=0;j<matches.size();j++)
02097                 {
02098                         if(i==j)continue;
02099                         if(matches[i].first==matches[j].first)mult0++;
02100                         if(matches[i].second==matches[j].second)mult1++;
02101                 }
02102                 if(! ( mult0>0 && mult1>0 )) // no degeneracy
02103                 {
02104                         matchestmp.push_back(matches[i]);
02105                         multu.push_back(mult0);
02106                         multv.push_back(mult1);
02107                 }else if( mult0>0 && mult1>0 ) //degeneracy, but see if we really want it
02108                 {
02109                         //see if the reason why both of these are degenerate is because they are really long
02110                         //for example, this can happen to a long muon in a DIS CC
02111                         //we want each of these to be the longest out of all matches using each of these paths
02112                         int pass=1;
02113                         
02114                         for(unsigned int k=0;k<matches.size();k++)
02115                         {
02116                                 if(matches[i].second!=matches[k].second)continue;
02117                                 if(pu[matches[i].first].end_z-pu[matches[i].first].start_z < pu[matches[k].first].end_z-pu[matches[k].first].start_z)
02118                                 {
02119                                         pass=0;
02120                                         break;
02121                                 }
02122                         }
02123                         for(unsigned int k=0;k<matches.size();k++)
02124                         {
02125                                 if(matches[i].first!=matches[k].first)continue;
02126                                 if(pv[matches[i].second].end_z-pv[matches[i].second].start_z < pv[matches[k].second].end_z-pv[matches[k].second].start_z)
02127                                 {
02128                                         pass=0;
02129                                         break;
02130                                 }
02131                         }
02132                 
02133                         if(pass)
02134                         {
02135                                 matchestmp.push_back(matches[i]);
02136                                 multu.push_back(mult0);
02137                                 multv.push_back(mult1);
02138                         }
02139                 }
02140                 
02141                 MSG("Finder",Msg::kDebug) << "pair " << matches[i].first <<" " << matches[i].second << " matches " <<mult0 <<" "<< mult1<<"\n";
02142         }
02143         matches=matchestmp;
02144 
02145         //cout <<"*********Currently have "<<mypoh->particles3d.size()<<" particles in POH\n";
02146         
02147         //we now have matches...
02148         //for each, construct a 3d path view
02149         
02150         
02151         std::vector<Particle3D> p3v;
02152         for(unsigned int i=0;i<matches.size();i++)
02153         {
02154                 MSG("Finder",Msg::kDebug) << "Making 3d particle with paths "<<matches[i].first<<" "<<matches[i].second<<"\n";
02155         
02156                 Particle3D p3 = Make3DParticle(pu[matches[i].first].path, pv[matches[i].second].path, chu, chv,multu[i],multv[i]) ;
02157                 
02158                 //we also need to save the clusters....
02159                 
02160                 for(unsigned int j=0;j<p3.view.size();j++)
02161                 {
02162                         Managed::ClusterManager * cluster_manager = &mypoh->clusterer;//p3.view[j]==2?&clustermanager_u:&clustermanager_v;
02163                         if(p3.chainhit[j]>-1)
02164                         {
02165                                 ChainHelper *ch = p3.view[j]==2?chu:chv;
02166                                 int oldid = ch->GetChain(p3.chain[j])->cluster_id[p3.chainhit[j]];
02167                                 int newid = cluster_manager->SaveCluster(oldid,p3.e[j],3);
02168                                 
02169                                 Managed::ManagedCluster * cluster = cluster_manager->GetCluster(newid);
02170                                 if(cluster)
02171                                 {
02172                                         //printf("SSSS saved %d  with e %f new id %d with e %f \n",oldid,p3.e[j],newid,cluster->e);
02173                                 
02174                                         p3.e[j]=cluster->e;
02175                                 
02176                                         p3.chainhit[j]=newid;
02177                                         
02178                                         ch->ChangeHitId(oldid,newid);
02179                         
02180                                 }else{
02181                                         //printf("NNNN cluster %d doesn't exist! removing e %f from particle \n",oldid,p3.e[j]);
02182                                         p3.e[j]=0;
02183                                         p3.chainhit[j]=-1;
02184                                 
02185                                         
02186                                 }
02187                         }
02188                         
02189                                 
02190                 
02191                 }
02192                 
02193                 if(p3.entries && p3.sum_e>1e-5)
02194                 {
02195                         //printf("adding p3d with e %f\n",p3.sum_e);
02196                         p3v.push_back(p3);
02197                 }
02198         }
02199 
02200         
02201 //printf("!!2 p3 size %d \n",p3v.size());
02202 
02203         
02204         
02205         
02206         
02207         
02208         
02209 
02210         std::vector<Particle3D* > p3d;  
02211         for(unsigned int i=0;i<p3v.size();i++)
02212                 p3d.push_back(&p3v[i]); 
02213         
02214 //      p3d=SetShared(p3d);     
02215 
02216 //      shareholder.Reset();
02217 //      shareholder.BulkInsert(p3d);
02218 
02219         /*
02220         DumpParticle3D(p3d);
02221 
02222 
02223         std::vector<Particle3D*> pout;
02224         std::vector<Particle3D*> toshare;
02225         for(unsigned int i=0;i<p3d.size();i++)
02226         {
02227                 //cout<<"on particle3d "<<i<<endl;
02228         
02229                 std::vector<Particle3D* > po = ProcessParticle3D1(p3d[i]);
02230                 for(unsigned int j=0;j<po.size();j++)
02231                 {
02232                         //if(po[j]->hasShared()<1)
02233                                 pout.push_back(po[j]);  
02234                         //else toshare.push_back(po[j]);
02235                 }
02236         
02237         }
02238         */
02239         
02240 /*      
02241         cout <<"PROCESSING COMPLETE\n";
02242         cout <<"not shared\n";
02243                 DumpParticle3D(pout);
02244         cout <<"has shared\n";
02245 
02246                 DumpParticle3D(toshare);
02247                 
02248         cout << "UNSHARING ENERGY\n";
02249 */      
02250         
02251 
02252 /*
02253 
02254         int its=5; //max iterations to try to remove shared hits!
02255         while(toshare.size()>0 && its>0)
02256         {
02257                 std::vector<Particle3D*> tmp;
02258                 
02259                 //some code to process them, returning a vector....
02260                 
02261                 for(unsigned int i=0;i<toshare.size();i++)
02262                         tmp.push_back(toshare[i]);
02263                 toshare.clear();
02264                 
02265                 RemoveSharedHits(tmp);
02266                 
02267                 for(unsigned int i=0;i<tmp.size();i++)
02268                 {
02269                         if(tmp[i]->hasShared())
02270                                 toshare.push_back(tmp[i]);
02271                         else pout.push_back(tmp[i]);
02272                 
02273                 }
02274                 
02275                 its--;
02276         }
02277         
02278         if(toshare.size()>0)
02279         MSG("Finder",Msg::kWarning)<<"particles lost .... could not uniquely assign shared energy !"<<endl;
02280         
02281         */
02282 
02283 
02284 /*
02285         std::vector<Particle3D*> pout;
02286         for(unsigned int i=0;i<p3d.size();i++)
02287         {
02288                 std::vector<Particle3D* > po = ProcessParticle3D(p3d[i]);       //returned vector may not have shared hits....
02289                 
02290                 std::vector<Particle3D*> forsharing;
02291                 
02292                 for(unsigned int j=0;j<po.size();j++)
02293                         pout.push_back(po[j]);
02294                 
02295                 for(unsigned int j=0;j<pout.size();j++)
02296                         forsharing.push_back(pout[j]);
02297                 
02298                 for(unsigned int j=i+1;j<p3v.size();j++)
02299                         forsharing.push_back(p3d[j]);
02300 
02301                         
02302                 forsharing=SetShared(forsharing); //to reset if a hit is shared....
02303                 
02304                 
02305         }
02306         
02307         
02308         
02309 */      
02310 /*
02311         cout <<"CLEANING\n";
02312         for(unsigned int i=0;i<pout.size();i++)pout[i]->Clean(); //remove 0 energy hits...
02313         DumpParticle3D(pout);
02314 */      
02315 
02316 
02317 
02318         //we need to recalculate values after sharing hits
02319         //but don't try to find more unique muons
02320         //we just want to make sure that the type hasn't changed!
02321         std::vector<Particle3D*> pout1;
02322         for(unsigned int i=0;i<p3d.size();i++)
02323         {
02324                 std::vector<Particle3D* > po = ProcessParticle3D1(p3d[i],0);
02325                 for(unsigned int j=0;j<po.size();j++)
02326                 {
02327                         pout1.push_back(po[j]); 
02328                 }
02329         
02330         }
02331         
02332         //printf("weave found %d\n",pout1.size());
02333 
02334 //      pout1=p3d;
02335         //look for straglers (large hit, not associated with any chain.... condider these to be neutrons...)
02336         FindNeutrons(pout1);    
02337         //printf("weave found %d\n",pout1.size());
02338         
02339 
02340         FinalizeParticles3D(pout1); //final computation of values, etc...
02341         //printf("weave found %d\n",pout1.size());
02342 
02343 
02344         for(unsigned int i=0;i<pout1.size();i++)
02345         {
02346                 if(pout1[i]->sum_e>0.01)
02347                         mypoh->AddParticle3D(*(pout1[i]));
02348         }
02349         
02350         //printf("weave found %d\n",pout1.size());
02351 
02352                 
02353         DumpParticle3D(pout1);  
02354 
02355                 
02356         //RecordLostHits(pout1);        
02357                 
02358 //              printf("!!3 pout size %d\n",pout.size());
02359                 
02360                 
02361 //              printf("poh entries %d\n", mypoh->particles3d1->GetEntries());
02362                                 
02363 
02364 //      for(unsigned int i=0;i<p3d.size();i++)
02365 //              mypoh->AddParticle3D(*(p3d[i]));        
02366 
02367 
02368 
02369         //do a check to see if the total particle energy is less than the total event energy...
02370         double parte=0;
02371         for(unsigned int i=0;i<p3d.size();i++) parte+=p3d[i]->sum_e;
02372         
02373         if(parte-mypoh->event.visenergy>.001)printf("!!!!!!!!!!!!!!!!!!\n********************\nenergy mismatch!!!! particle sum %f   energy total %f\n!!!!!!!!!!!!!!!!!!\n********************\n",parte,mypoh->event.visenergy);
02374         
02375         
02376 
02377         
02378 }


Member Data Documentation

Definition at line 88 of file Finder.h.

Referenced by SetClusterManagerU().

Definition at line 89 of file Finder.h.

Referenced by Process(), and SetClusterManagerV().

int Finder::DoMRCC [private]

Definition at line 146 of file Finder.h.

Referenced by Finder(), Process(), and SetMRCC().

std::vector<double> Finder::energy [private]

Definition at line 111 of file Finder.h.

Referenced by AddStrip(), ClearXTalk(), FindIsolatedHits(), Reset(), and ShareHit().

std::map<int, std::map<int, int> > Finder::loc_map [private]

Definition at line 119 of file Finder.h.

Referenced by AddStrip(), ClearXTalk(), and Reset().

double Finder::meupergev [private]

Definition at line 143 of file Finder.h.

Referenced by FinalizeParticles3D(), Finder(), Process(), and SetMEUperGeV().

std::vector<Particle> Finder::particles [private]

Definition at line 115 of file Finder.h.

Referenced by FindIsolatedHits(), and Reset().

std::vector<int> Finder::plane [private]

Definition at line 106 of file Finder.h.

Referenced by AddStrip(), ClearXTalk(), FindIsolatedHits(), GetStrips(), and Reset().

Definition at line 145 of file Finder.h.

Referenced by DumpParticle3D(), RemoveSharedHits(), Reset(), and StripMuon1().

std::multimap<double,int> Finder::sorter_map [private]

Definition at line 117 of file Finder.h.

Referenced by AddStrip(), ClearXTalk(), Finder(), FindIsolatedHits(), Reset(), and ~Finder().

std::vector<int> Finder::strip [private]

Definition at line 107 of file Finder.h.

Referenced by AddStrip(), ClearXTalk(), FindIsolatedHits(), and Reset().

std::vector<double> Finder::t [private]
double Finder::true_vtx_u [private]

Definition at line 130 of file Finder.h.

Referenced by Process(), Reset(), and SetTrueVertex().

double Finder::true_vtx_v [private]

Definition at line 131 of file Finder.h.

Referenced by Process(), Reset(), and SetTrueVertex().

double Finder::true_vtx_z [private]

Definition at line 132 of file Finder.h.

Referenced by Process(), Reset(), and SetTrueVertex().

std::vector<int> Finder::view [private]

Definition at line 105 of file Finder.h.

Referenced by AddStrip(), ClearXTalk(), Make3DParticle(), Reset(), SetShared(), and shareEnergy().

double Finder::vtx_u
double Finder::vtx_v
double Finder::vtx_z
std::vector<double> Finder::z [private]

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

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1