LongMuonFinder Class Reference

#include <LongMuonFinder.h>

List of all members.

Public Member Functions

 LongMuonFinder (Detector::Detector_t d)
 ~LongMuonFinder ()
int FindLongMuon (Managed::ClusterManager *cm)
Particle3DGetParticle3D ()
ChainGetChain (int view)
int FoundSingleViewLongMuon ()

Private Member Functions

int CheckChainQuality (Chain *c, int view, int partialcount)
ChainFindMuonChain (Managed::ClusterManager *cl, int view)
void Reset ()
void MakeParticle3D ()
void AbsorbMuonClusters (Chain *c, int view, double past_z)
void RemoveNonMuonEnergy (Chain *c, int view, double past_z)
std::pair< int, int > CountInPartiallyInstrumentedRegion (Chain *viewU, Chain *viewV)
void DumpParticle3D ()
double FindIsolationZ ()
void MergeChainClusters (Chain *ch, std::vector< int > *clusters)
int CheckChainOverlap (Chain *chain_u, Chain *chain_v, double isolation_z)
void ClearFrontVertex (Chain *chain_u, Chain *chain_v)
int IsPartiallyInstrumented (double t, double z, int view)

Private Attributes

int foundparticle
Managed::ClusterManagercluster_manager
int single_view_long_muon
Detector::Detector_t detector

Static Private Attributes

static Particle3Dfoundparticle3d = 0
static Chainchain_u = 0
static Chainchain_v = 0

Detailed Description

Definition at line 10 of file LongMuonFinder.h.


Constructor & Destructor Documentation

LongMuonFinder::LongMuonFinder ( Detector::Detector_t  d  ) 

Definition at line 194 of file LongMuonFinder.cxx.

References detector, and Reset().

00195 {
00196         detector=d;
00197         Reset();
00198 }

LongMuonFinder::~LongMuonFinder (  ) 

Definition at line 200 of file LongMuonFinder.cxx.

00201 {}


Member Function Documentation

void LongMuonFinder::AbsorbMuonClusters ( Chain c,
int  view,
double  past_z 
) [private]

Definition at line 450 of file LongMuonFinder.cxx.

References cluster_manager, Managed::ClusterManager::GetClusterMap(), and MergeChainClusters().

Referenced by FindLongMuon().

00451 {
00452         //for each cluster in the chain, merge it with all other clusters in the plane past z distance 
00453 
00454         if(!cluster_manager)return;
00455         
00456         std::map<double, std::map<double, int> > * cluster_map = cluster_manager->GetClusterMap(view);
00457 
00458     std::map<double, std::map<double, int> >::iterator p_iterr;
00459     std::map<double, int >::iterator s_iterr;
00460 
00461         std::vector<int> clusters;
00462 
00463         double last_z=0;
00464         for(p_iterr=cluster_map->begin();p_iterr!=cluster_map->end(); p_iterr++)
00465     {   
00466         if(p_iterr->first-past_z<0.02)continue;
00467                 
00468         if(fabs(last_z-p_iterr->first)>0.01)
00469                 {
00470                         if(clusters.size()>0)MergeChainClusters(ch,&clusters);
00471                         clusters.clear();
00472                         last_z=p_iterr->first;
00473                 }
00474         
00475         
00476         
00477         for(s_iterr=p_iterr->second.begin();s_iterr!=p_iterr->second.end(); s_iterr++)
00478         {
00479                         clusters.push_back(s_iterr->second);
00480                 }
00481         }
00482         
00483         if(clusters.size()>0)MergeChainClusters(ch,&clusters);
00484         clusters.clear();
00485 }

int LongMuonFinder::CheckChainOverlap ( Chain chain_u,
Chain chain_v,
double  isolation_z 
) [private]

Definition at line 365 of file LongMuonFinder.cxx.

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

Referenced by FindLongMuon().

00366 {
00367 
00368         MSG("LongMuonFinder",Msg::kDebug)<<"checking chain overlap clear_z at "<<clear_z<<"\n";
00369         MSG("LongMuonFinder",Msg::kDebug)<<"isolation distance u "<<chain_u->end_z-clear_z<<"\n";
00370         MSG("LongMuonFinder",Msg::kDebug)<<"isolation distance v "<<chain_v->end_z-clear_z<<"\n";
00371 
00372 
00373         //first check that the muon is sufficiently long past the isolation point!
00374         double require_isolation_distance=1.0;
00375         if(clear_z <1e-9)return 0; //no valid z clearance measured!
00376         if(chain_u->end_z-clear_z<require_isolation_distance && 
00377                 fabs(chain_u->end_z-chain_u->start_z)*0.7 > clear_z)return 0;
00378         if(chain_v->end_z-clear_z<require_isolation_distance && 
00379                 fabs(chain_v->end_z-chain_v->start_z)*0.7 > clear_z)return 0;
00380         
00381         //now make sure that there is sufficient overlap
00382         double require_overlap_distance=1.0;
00383         
00384         double start_z = chain_u->start_z > chain_v->start_z ? chain_u->start_z : chain_v->start_z;
00385         double end_z = chain_u->end_z < chain_v->end_z ? chain_u->end_z : chain_v->end_z;
00386 
00387         MSG("LongMuonFinder",Msg::kDebug)<<"overlap "<<end_z-start_z<<"\n";
00388         
00389         if(end_z-start_z < require_overlap_distance)return 0;
00390         return 1;
00391         
00392         
00393 
00394 }

int LongMuonFinder::CheckChainQuality ( Chain c,
int  view,
int  partialcount 
) [private]

Definition at line 633 of file LongMuonFinder.cxx.

References Chain::cluster_id, cluster_manager, Chain::e, Chain::entries, Managed::ClusterManager::GetCluster(), Managed::ManagedCluster::hitplane, Msg::kDebug, max, min, MSG, and Chain::z.

Referenced by FindLongMuon().

00634 {
00635         if(!c)return 0;
00636         if(c->entries<3)return 0;       
00637         
00638         int qual=1;
00639         
00640         //look at muon-like strip fraction
00641         //and sparsity (# planes with hits/#planes from front to back of muon)
00642         
00643         int muonlike=0;
00644         double hitplane=0.;
00645         
00646         double last_hitplane=-100.;
00647         
00648         double min = 1000000.;
00649         double max = 0.;
00650         
00651         double partial=(double)partialcount;
00652         
00653         for(int i=0;i<c->entries;i++)
00654         {
00655                 Managed::ManagedCluster *clus = cluster_manager->GetCluster(c->cluster_id[i]);
00656         
00657                 if(c->e[i]>0.5 && c->e[i]<2.5)muonlike++;
00658                 if(fabs(last_hitplane - c->z[i])>0.03)
00659                 {
00660                         hitplane+=1;
00661                         last_hitplane=c->z[i];
00662                 }
00663                 
00664                 min = clus->hitplane[0]< min? clus->hitplane[0]:min;
00665                 max = clus->hitplane[0]> max? clus->hitplane[0]:max;
00666         
00667 /*              int part=IsPartiallyInstrumented(c->t[i], c->z[i], view);
00668                 if(part)
00669                 {
00670                         //if(part==1) //its in z, so no question!
00671                                 partial+=1;
00672                         //if(part==2)
00674                         //so we should just err on the side of caution and count all failures...
00675                         //this will result in sparsities > 1 !
00676                         //this will also result in under counted sparsities!
00677                 }
00678 */      }
00679         
00680         double muonfrac = (double)muonlike / (double)c->entries;
00681         
00682         double nplanes = (max-min)/2.0 ; //get the total number of planes in this view that the track could go through// / 0.0708  ;
00683         
00684         //double detscale=1;
00685         //if(partial>0)detscale=4.*partial/hitplane + 1.*(1-partial/hitplane);
00686         
00687         
00688         //printf("partially instrumented planes %f max %f min %f estplanes %f hit planes %f adj planes %f\n",partial,max,min, (max-min)/0.0708, hitplane,hitplane+partial*8.0);
00689         
00690         //adjust number of hitplanes for partially instrumented regions
00691         hitplane += 8*partial;
00692         
00693         double sparsity = nplanes<=0? 1 : (double)hitplane / nplanes ; //length/size of u+v plane
00694         
00695         double sparsity_cut=0.8;
00696         
00697         //we don't know about the other view... so if we have some partial hits....make a large guess...count the rest as partial as well!
00698         if(partial/(hitplane-8*partial)>0.2 && sparsity<1)sparsity=(hitplane+8*(hitplane-partial))/nplanes;
00699         
00700         
00701         if (muonfrac < 0.5) qual=0;
00702         if (sparsity < sparsity_cut) qual=0;
00703         
00704         MSG("LongMuonFinder",Msg::kDebug)<<"muon chain check  muonfrac "<< muonfrac<<" sparsity "<< sparsity<<"  qual "<< qual<<" hitplane "<<hitplane<< " nplane "<<nplanes <<" partialhits "<<partial<<"\n";
00705 
00706         return qual;
00707 }

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

Definition at line 339 of file LongMuonFinder.cxx.

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

Referenced by FindLongMuon().

00340 {
00341         double start_z_u = chain_u->start_z;
00342         double start_z_v = chain_v->start_z;
00343 
00344         if(fabs(start_z_u-start_z_v)<0.04)return; //close enough
00345         
00346         double start_z = start_z_u < start_z_v ? start_z_v : start_z_u - 0.05;
00347         
00348         Chain * toclear =  start_z_u < start_z_v ? chain_u : chain_v;
00349         //now delete all hits up to start_z
00350         
00351         Chain tmp = *toclear;
00352         
00353         toclear->ClearHits();
00354         
00355         for(int i=0;i<tmp.entries;i++)
00356         {
00357                 if(tmp.z[i]>start_z)
00358                         toclear->insert(tmp.t[i], tmp.z[i], tmp.e[i], tmp.cluster_id[i]);
00359         }
00360         
00361 }

std::pair< int, int > LongMuonFinder::CountInPartiallyInstrumentedRegion ( Chain viewU,
Chain viewV 
) [private]

Definition at line 65 of file LongMuonFinder.cxx.

References detector, IsPartiallyInstrumented(), it, Detector::kNear, Chain::t, and Chain::z.

Referenced by FindLongMuon().

00066 {
00067         int ucount=0;
00068         int vcount=0;
00069         
00070         if(detector!=Detector::kNear)return std::pair<int,int>(ucount,vcount);
00071 
00072         std::map<double,std::pair<double,int> > zorder;// z, <t, view>
00073         
00074         for(unsigned int i=0;i<viewU->z.size();i++)
00075                 zorder.insert(std::pair<double,std::pair<double,int> >(viewU->z[i],std::pair<double,int>(viewU->t[i],2)));
00076                 
00077         for(unsigned int i=0;i<viewV->z.size();i++)
00078                 zorder.insert(std::pair<double,std::pair<double,int> >(viewV->z[i],std::pair<double,int>(viewV->t[i],3)));
00079 
00080         
00081 
00082         //count the number in each plane that are in the partially instrumented region
00083         
00084         std::map<double,std::pair<double,int> >::iterator it;
00085         
00086         for(it=zorder.begin();it!=zorder.end();it++)
00087         {
00088                 //values in the other views for extrapolation!
00089                 double prevz=0;
00090                 double prevt=0;
00091                 double nextz=0;
00092                 double nextt=0;
00093 
00094                 double prevz2=0;
00095                 double prevt2=0;
00096                 double nextz2=0;
00097                 double nextt2=0;
00098                 
00099                 int thisview = it->second.second;
00100                 double thist = it->second.first;
00101                 double thisz = it->first;
00102                 
00103                 if(thisview !=2 && thisview !=3)
00104                 {
00105                         printf("unknown view\n");
00106                         abort();
00107                 }
00108                 
00109                 std::map<double,std::pair<double,int> >::iterator itforward;
00110                 std::map<double,std::pair<double,int> >::iterator itbackward;
00111                 
00112                 for(itforward=it;itforward!=zorder.end();itforward++)
00113                 {
00114                         if(itforward->second.first!=thisview )
00115                         {
00116                                 if(!nextz)
00117                                 {
00118                                         nextz=itforward->first;
00119                                         nextt=itforward->second.second;
00120                                 }else if(fabs(itforward->first - nextz)>0.04){  //make sure we are in a different plane
00121                                         nextz2=itforward->first;
00122                                         nextt2=itforward->second.second;                                
00123                                         break;
00124                                 }
00125                         }       
00126                 }
00127                 for(itbackward=it;itbackward!=zorder.begin();itbackward--)
00128                 {
00129                         if(itbackward->second.first!=thisview)
00130                         {
00131                                 if(!nextz)
00132                                 {
00133                                         prevz=itbackward->first;
00134                                         prevt=itbackward->second.second;
00135                                 }else if(fabs(itbackward->first - prevz)>0.04){  //make sure we are in a different plane
00136                                         prevz2=itbackward->first;
00137                                         prevt2=itbackward->second.second;                               
00138                                         break;
00139                                 }
00140                         }
00141                 }
00142                 
00143         //      printf("n %f %f p %f %f\n",nextz,nextz2,prevz,prevz2);
00144                 
00145                 if((nextz?1:0)+ (nextz2?1:0)+ (prevz?1:0)+ (prevz2?1:0) <2)break;//we can't do anything with information in only one view! need at least two points in the other view!
00146                 
00147                 //we need an estimation of the other view t position at the point in question...
00148                 
00149                 double est_other_t=0;
00150                 
00151                 double slope=0;
00152                 double offset=0;
00153                 
00154                 if(nextz && prevz)
00155                 {
00156                         slope = (nextt-prevt)/(nextz-prevz);
00157                         offset = prevt-slope*prevz;
00158                 }else if(nextz && nextz2)
00159                 {
00160                         slope = (nextt-nextt2)/(nextz-nextz2);
00161                         offset = nextt2-slope*nextz2;
00162                 }else if(prevz && prevz2)
00163                 {
00164                         slope = (prevt-prevt2)/(prevz-prevz2);
00165                         offset = prevt2-slope*prevz2;
00166                 }
00167                 
00168                 est_other_t = thisz*slope+offset;
00169                 
00170         //      printf("checking view %d z %f thist %f othert %f\n",thisview,thisz,thist,est_other_t);
00171                 
00172                 int partialviewthis = IsPartiallyInstrumented(thist,thisz,thisview);
00173                 int partialviewother = IsPartiallyInstrumented(est_other_t,thisz,thisview==2?3:2);
00174                 
00175                 if(partialviewthis || partialviewother)
00176                 {
00177                         if(thisview==2)ucount++;
00178                         else vcount++;
00179                 }
00180                 
00181         }
00182         
00183         
00184         
00185         
00186         return std::pair<int,int>(ucount,vcount);
00187         
00188                                 
00189 
00190 }

void LongMuonFinder::DumpParticle3D (  )  [private]

Definition at line 1499 of file LongMuonFinder.cxx.

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

Referenced by FindLongMuon().

01500 {
01501         ostringstream s;
01502         
01503         s << "Long Muon Particle3D found "<<endl;
01504 
01505                 Particle3D * p = foundparticle3d;
01506                 if (p==0)return;
01507         
01508         
01509 
01510         
01511                 s << "\n---Particle3d "  << "---\nstart, stop (u,v,z) (" << p->start_u << ", "<< p->start_v << ", " << p->start_z <<") (" <<p->end_u<<", "<< p->end_v << ", " << p->end_z <<")"<<endl;
01512                 s << "entries "<< p->entries << " muonlike " << p->muonfrac <<endl;
01513                 
01514 
01515                 
01516                 s<<"types: ";
01517                 for(unsigned int j=0;j<p->types.size();j++)
01518                 {
01519                         switch( p->types[j].type)
01520                         {
01521                                 case    ParticleType::em:
01522                                         s<<"em ";
01523                                         break;
01524                                 case ParticleType::emshort:
01525                                         s<<"emshort  ";
01526                                         break;                          
01527                                 case ParticleType::muon:
01528                                         s<<"muon ";
01529                                         break;
01530                                 case ParticleType::prot:
01531                                         s<<"prot ";
01532                                         break;          
01533                                 case ParticleType::pi0:
01534                                         s<<"pi0 ";
01535                                         break;
01536                                 case ParticleType::uniquemuon:
01537                                         s<<"uniquemuon ";
01538                                         break;  
01539                         
01540                         }
01541                                 
01542                 }
01543                 s<<endl;
01544                 
01545                 s<<"particletype : "<<p->particletype<<endl;
01546                 
01547                 
01548                 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;
01549                 s<<"emprob " << p->emfit_prob <<" avg rms_t "<<p->avg_rms_t<<endl;
01550                                 
01551                 s << "points (u,v,z,e - chain, chainhit, chainview - shared - rms_t - view) : ";
01552                 for(int j=0;j<p->entries;j++)
01553                         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]<<") ";
01554         
01555         /*
01556         s << "points (e - shared) : ";
01557                 for(int j=0;j<p->entries;j++)
01558                         s << "(" <<p->e[j]<<" - "<<p->shared[j]<<") ";
01559         
01560         */
01561         
01562         
01563 
01564                         
01565                         
01566 
01567                         
01568                         
01569                 s<<endl<<endl;
01570                 
01571                 
01572         
01573         
01574         
01575 
01576         MSG("LongMuonFinder",Msg::kDebug) <<s.str();
01577 
01578         
01579 }

double LongMuonFinder::FindIsolationZ (  )  [private]

Definition at line 552 of file LongMuonFinder.cxx.

References cluster_manager, Managed::ManagedCluster::e, Managed::ClusterManager::GetCluster(), Managed::ClusterManager::GetClusterMap(), Msg::kDebug, MSG, and Managed::ManagedCluster::view.

Referenced by FindLongMuon().

00553 {
00554         if(!cluster_manager)return 0;
00555         
00556         //require 3 consecutive u and v planes to signify start of muon isolation
00557         std::map<double, std::map<double, int> > * cluster_map = cluster_manager->GetClusterMap();
00558 
00559 
00560     std::map<double, std::map<double, int> >::iterator p_iterr;
00561     std::map<double, int >::iterator s_iterr;
00562 
00563         double start_z[6];
00564         for(int i=0;i<6;i++)start_z[i]=0;
00565         int ucount=0;
00566         int vcount=0;
00567         int cnt=0;
00568         int lastview=0;
00569         
00570         std::vector<double> goodz;
00571         
00572     for(p_iterr=cluster_map->begin();p_iterr!=cluster_map->end(); p_iterr++)
00573     {   
00574         
00575                 if(fabs(start_z[5]-p_iterr->first)>0.03)
00576                 {
00577                 
00578                         
00579                         for(int i=0;i<5;i++)start_z[i]=start_z[i+1];
00580                         start_z[5]=p_iterr->first;
00581                         
00582                         if(cnt<=1)
00583                         {
00584                                 if(lastview==2)ucount++;
00585                                 else vcount++;
00586                                 
00587                                 goodz.push_back(p_iterr->first);
00588                         }else{
00589                                 ucount=0;
00590                                 vcount=0;
00591                                 goodz.clear();
00592                         }
00593                         
00594                         MSG("LongMuonFinder",Msg::kDebug)<<"z "<<start_z[5] <<" lastview "<<lastview <<", ucount "<< ucount<<" vcount "<<vcount <<" pcnt "<< cnt<<"\n";
00595                                 
00596                         cnt=0;
00597                 }
00598 
00599         if(ucount>=3 && vcount>=3)break;
00600 
00601         for(s_iterr=p_iterr->second.begin();s_iterr!=p_iterr->second.end(); s_iterr++)
00602         {
00603                 Managed::ManagedCluster *this_cluster = cluster_manager->GetCluster(s_iterr->second);
00604                         int view=this_cluster->view;
00605                         
00606                 //did we have an empty plane?
00607                 //or do we have a very large, unmuon type hit?
00608                 if( /*(cnt==0 && view==lastview) ||*/ this_cluster->e >3)
00609                 {
00610                                 ucount=0;
00611                                 vcount=0;
00612                                 goodz.clear();
00613                 }
00614                 
00615                 
00616         
00617                 lastview=view;
00618                 cnt++;
00619                 }
00620         
00621         }
00622 
00623         MSG("LongMuonFinder",Msg::kDebug)<<"!!!!! u "<<ucount<<" v "<<vcount<<" z "<<start_z[0]<<" good z "<<goodz[0]<<"\n";
00624         if(ucount>=3 && vcount>=3) return goodz[0]; //start_z[0]
00625         
00626         return 0;       
00627 
00628 }

int LongMuonFinder::FindLongMuon ( Managed::ClusterManager cm  ) 

Definition at line 212 of file LongMuonFinder.cxx.

References AbsorbMuonClusters(), Chain::available, chain_u, chain_v, CheckChainOverlap(), CheckChainQuality(), ClearFrontVertex(), Chain::cluster_id, cluster_manager, CountInPartiallyInstrumentedRegion(), DumpParticle3D(), Managed::ManagedCluster::e, Chain::e, Chain::end_z, Chain::entries, FindIsolationZ(), FindMuonChain(), foundparticle, foundparticle3d, Managed::ClusterManager::GetSavedCluster(), Msg::kDebug, MakeParticle3D(), MSG, Particle3D::muon, Particle3D::particletype, Chain::PrintChain(), Chain::Recalc(), RemoveNonMuonEnergy(), Reset(), Managed::ClusterManager::SaveCluster(), single_view_long_muon, and Chain::start_z.

00213 {
00214 
00215         MSG("LongMuonFinder",Msg::kDebug)<<"looking for long muons\n";
00216         Reset();
00217 
00218         cluster_manager=cm;
00219         
00220         //try to find a long muon chain in each view
00221         if(chain_u){delete chain_u;chain_u=0;}
00222         if(chain_v){delete chain_v;chain_v=0;}
00223         chain_u = FindMuonChain(cm,2);
00224         chain_v = FindMuonChain(cm,3);
00225         
00226         int qual_u=0;
00227         int qual_v=0;
00228         
00229         
00230         
00231         if(chain_u)chain_u->Recalc();
00232         if(chain_v)chain_v->Recalc();
00233         
00234         std::pair<int,int> partialcount = CountInPartiallyInstrumentedRegion(chain_u, chain_v);
00235         
00236         //printf("partial count %d %d\n",partialcount.first,partialcount.second);
00237         
00238         if(chain_u)qual_u=CheckChainQuality(chain_u,2,partialcount.first);
00239         if(chain_v)qual_v=CheckChainQuality(chain_v,3,partialcount.second);
00240         
00241 
00242         if(qual_u&&!qual_v)
00243         {
00244                 if(chain_u->end_z - chain_u->start_z > 2)single_view_long_muon=1;
00245         }else if(qual_v&&!qual_u)
00246         {
00247                 if(chain_v->end_z - chain_v->start_z > 2)single_view_long_muon=1;
00248         }
00249         
00250         //if we have a quality chain in each view... then we will be making a particle
00251         if(qual_u && qual_v)
00252         {
00253         
00254                 MSG("LongMuonFinder",Msg::kDebug)<<"qqqqq "<<qual_u <<" "<< qual_v<<"\n";
00255                 
00256         
00257                 //look in both views to find at least 3u and 3v consecutive planes with only 1 strip... that is where we say the muon exits the rest of the shower/partices/etc
00258                 double isolation_z = FindIsolationZ();
00259                 
00260 
00261                 int qual=CheckChainOverlap(chain_u, chain_v, isolation_z);
00262                 if(qual)
00263                 {
00264                 
00265                 ClearFrontVertex(chain_u, chain_v);
00266 
00267                 if(isolation_z>0)
00268                 {               
00269                         //need to remove energy from clusters forward of that point which are too high for a muon track
00270                         
00271                         RemoveNonMuonEnergy(chain_u,2,isolation_z);
00272                         RemoveNonMuonEnergy(chain_v,3,isolation_z);
00273                         
00274                 
00275                         //need to absorb clusters along the chain which came off the muon (brehms, etc)
00276                         AbsorbMuonClusters(chain_u,2,isolation_z);
00277                         AbsorbMuonClusters(chain_v,3,isolation_z);
00278                 }
00279         
00280 
00281                 //save clusters involved...
00282 
00283                 for(int i=0;i<chain_u->entries;i++)
00284                 {
00285                         int newid = cluster_manager->SaveCluster(chain_u->cluster_id[i],chain_u->e[i],1);
00286                         Managed::ManagedCluster * newcluster = cluster_manager->GetSavedCluster(newid);
00287                         if(chain_u->e[i]!=newcluster->e)
00288                                 MSG("LongMuonFinder",Msg::kDebug)<<"saving with different e "<<chain_u->e[i]<<" "<<newcluster->e<<"\n";
00289                         chain_u->e[i]=newcluster->e;
00290                         chain_u->cluster_id[i]=newid;
00291                 }
00292 
00293                 for(int i=0;i<chain_v->entries;i++)
00294                 {
00295                         int newid = cluster_manager->SaveCluster(chain_v->cluster_id[i],chain_v->e[i],1);
00296                         Managed::ManagedCluster * newcluster = cluster_manager->GetSavedCluster(newid);
00297                         if(chain_v->e[i]!=newcluster->e)
00298                                 MSG("LongMuonFinder",Msg::kDebug)<<"saving with different e "<<chain_v->e[i]<<" "<<newcluster->e<<"\n";
00299                         chain_v->e[i]=newcluster->e;
00300                         chain_v->cluster_id[i]=newid;
00301                 }               
00302                 
00303                 
00304                 
00305                 
00306         
00307         
00308                 MSG("LongMuonFinder",Msg::kDebug)<<"making particle\n";
00309                 MakeParticle3D();
00310                 if(foundparticle3d)
00311                 {
00312                         foundparticle=1;
00313                 
00314                         foundparticle3d->particletype=Particle3D::muon;
00315                 
00316                         MSG("LongMuonFinder",Msg::kDebug)<<"particle made\n";
00317                         DumpParticle3D();
00318                         
00319                         chain_u->available=0;
00320                         chain_v->available=0;
00321                 }
00322                 }
00323         }
00324 
00325 
00326         MSG("LongMuonFinder",Msg::kDebug)<<"printing long muon chains\nu\n";
00327         if(chain_u)chain_u->PrintChain();MSG("LongMuonFinder",Msg::kDebug)<<"v\n";
00328         if(chain_v)chain_v->PrintChain();
00329         MSG("LongMuonFinder",Msg::kDebug)<<"done\n";
00330         
00331 
00332         MSG("LongMuonFinder",Msg::kDebug)<<"done looking for long muons\n";
00333         return foundparticle;
00334 }

Chain * LongMuonFinder::FindMuonChain ( Managed::ClusterManager cl,
int  view 
) [private]

Definition at line 983 of file LongMuonFinder.cxx.

References Managed::ManagedCluster::e, Chain::entries, Chain::front_offset, Chain::front_slope, Managed::ClusterManager::GetCluster(), Managed::ClusterManager::GetClusterMap(), Managed::ManagedCluster::id, Chain::insert(), Chain::interpolate(), Msg::kDebug, MSG, n, Chain::start_t, Chain::start_z, Managed::ManagedCluster::t, Chain::t, and Managed::ManagedCluster::z.

Referenced by FindLongMuon().

00984 {
00985 
00986         MSG("LongMuonFinder",Msg::kDebug)<<"\n\nLooking for long muon Chain in view "<<view<<"\n";
00987 
00988 
00989         std::map<double, std::map<double, int> > * cluster_map = cl->GetClusterMap(view);
00990 
00991 
00992     std::map<double, std::map<double, int> >::reverse_iterator p_iterr;
00993     std::map<double, int >::iterator s_iterr;
00994 
00995 
00996         Chain * c = new Chain();
00997 
00998         double last_t=0;
00999         double last_d=100000;
01000         Managed::ManagedCluster*closest=0;
01001         double last_z=100000;
01002         double last_plane=100000;
01003         Managed::ManagedCluster*last_closest=0;
01004         int notadded=0;
01005         int planecount=0;
01006         int lastplanecount=0;
01007         Managed::ManagedCluster*last_added=0;
01008         
01009         Managed::ManagedCluster * closest_proj=0;
01010         double last_d_proj=100000;
01011         
01012         double proj_t=0;
01013         double lin_proj_t=0;
01014         
01015         int first=1;
01016         
01017         
01018 /*      
01019         //for the first one......want to make sure we are choosing the hit closest to the track
01020         //go 5 planes and see where the track actually is...
01021 
01022         int npln=0;
01023         int ncnt=0;
01024         double mean_t=0;
01025     for(p_iterr=cluster_map->rbegin();p_iterr!=cluster_map->rend(); p_iterr++)
01026     {   
01027                 
01028 
01029                 if(fabs(last_plane-p_iterr->first)>0.03)npln++;
01030                 if(npln>5)break;
01031 
01032          for(s_iterr=p_iterr->second.begin();s_iterr!=p_iterr->second.end(); s_iterr++)
01033          {
01034                                 Managed::ManagedCluster * clus = cl->GetCluster(s_iterr->second);
01035                                 ncnt++;
01036                                 mean_t += clus->t;
01037         
01038                         if(npln==1)last_z=clus->z;
01039                 
01040          }
01041                 last_plane=p_iterr->first+0.04;
01042                 
01043                 
01044         }
01045         
01046         if(ncnt>0)mean_t/=ncnt;
01047         
01048         printf("last 5 planes mean t %f\n",mean_t);
01049 
01050         //find the closest cluster to the mean t within 2 strips
01051 
01052         first=1;
01053         last_t=mean_t;
01054 */
01055 
01056 
01057 
01058         
01059         
01060         //for the first one......want to make sure we are choosing the hit closest to the track
01061         //go 5 planes and see where the track actually is...
01062         //look for sigle hits in each view and extrapolate to end plane
01063 
01064         int npln=0;
01065         int ncnt=0;
01066 
01067 
01068                 double sum_z_t=0;
01069                 double sum_z=0;
01070                 double sum_t=0;
01071                 double sum_z_z=0;
01072                 
01073                 
01074         double lasts_t=0;
01075         double lasts_z=0;
01076         int n=0;
01077 
01078         double end_z=0;
01079 
01080     for(p_iterr=cluster_map->rbegin();p_iterr!=cluster_map->rend(); p_iterr++)
01081     {   
01082                 
01083 
01084                 if(fabs(last_plane-p_iterr->first)>0.03)
01085                 {
01086                         MSG("LongMuonFinder",Msg::kDebug)<<"newplane "<<npln<<" last plane count "<<ncnt<<"\n";
01087                         if(ncnt>0)
01088                         {
01089                                 lasts_z/=ncnt;
01090                                 lasts_t/=ncnt;
01091                         
01092                                 sum_z+=lasts_z;
01093                                 sum_t+=lasts_t;
01094                                 sum_z_t+=lasts_z*lasts_t;
01095                                 sum_z_z+=lasts_z*lasts_z;
01096                                 n++;
01097                         }
01098                 
01099                         npln++;
01100                         last_plane=p_iterr->first;
01101                         ncnt=0;
01102                         lasts_z=0;
01103                         lasts_t=0;
01104                 }
01105                 if(npln>5)break;
01106 
01107          for(s_iterr=p_iterr->second.begin();s_iterr!=p_iterr->second.end(); s_iterr++)
01108          {
01109                                 Managed::ManagedCluster * clus = cl->GetCluster(s_iterr->second);
01110                                 if(!clus)continue;
01111                                 ncnt++;
01112                 //      printf("cluster at z t e %f %f %f \n",clus->z,clus->t,clus->e);
01113                         lasts_z+=clus->z;
01114                         lasts_t+=clus->t;
01115                         if(clus->z>end_z)end_z=clus->z;
01116          }
01117                 
01118                 
01119                 
01120         }
01121 
01122                 
01123         double back_slope=0;
01124         double back_offset=0;
01125         last_t=0;
01126 
01127         if(n>1)
01128         {
01129                 back_slope=(n*sum_z_t-sum_z*sum_t)/(n*sum_z_z-(sum_z)*(sum_z));
01130                 back_offset=(sum_t*sum_z_z-sum_z*sum_z_t)/(n*sum_z_z-(sum_z)*(sum_z));
01131                 last_t=back_offset + back_slope*(end_z + 0.07); //want the proj hit in the plane past what we have in the data
01132                 last_z=last_z+0.07;
01133                 MSG("LongMuonFinder",Msg::kDebug)<<"!!! found back slope "<<back_slope<<" offset "<<back_offset<<" endz "<<end_z<<" proj "<<last_t<<" \n";
01134         
01135         }
01136         
01137         
01138 
01139         MSG("LongMuonFinder",Msg::kDebug)<<"last 5 planes proj t "<<last_t<<"\n";
01140 
01141         //find the closest cluster to the mean t within 2 strips
01142 
01143         first=1;
01144         last_plane=1000000;
01145 last_d=100000;
01146 
01147 
01148 /*      
01149         
01150         for(p_iterr=cluster_map->rbegin();p_iterr!=cluster_map->rend(); p_iterr++)
01151     {   
01152         if(fabs(last_plane-p_iterr->first)>0.03)npln++;
01153                 if(npln>5)break;
01154         for(s_iterr=p_iterr->second.begin();s_iterr!=p_iterr->second.end(); s_iterr++)
01155         {
01156                                 Managed::ManagedCluster * clus = cl->GetCluster(s_iterr->second);
01157                                 if(fabs(clus->t-mean_t)<0.05)
01158                                 {
01159                                         last_d=0;
01160                                         closest=clus;
01161                                 
01162                                 
01163                                         last_d_proj=0;
01164                                         closest_proj=clus;
01165                                         last_z=clus->z;
01166                                 
01167                                         printf("in first plane check %f %f  dist %f lastdist %f  lastdistproj %f\n",clus->z, clus->t,fabs(clus->t-lin_proj_t),last_d,last_d_proj);
01168                                                         
01169                                 planecount++;
01170                                 
01171                                         first=0;
01172                                 }
01173         }
01174                 last_plane=p_iterr->first;
01175                 
01176                 if(first)
01177                 {
01178                         p_iterr++;
01179                         break;
01180                 }
01181         }
01182 */      
01183         
01184         
01185         
01186         //end check
01187         
01188                         double last_notadded_z=0;
01189                 double last_notadded_t=0;
01190         
01191         
01192         int punch_through=0;
01193         double punch_through_t=0;
01194         double punch_through_z=0;
01195         double closest_punch_through=10000;
01196         
01197         //start in the plane after where we left off from starting hit check
01198     for(p_iterr=cluster_map->rbegin();p_iterr!=cluster_map->rend(); p_iterr++)
01199     {
01200 
01201                         
01202 
01203         //std::map<double, std::map<double, int> >::reverse_iterator checker = p_iterr;
01204         //checker ++;
01205         
01206         int nextplane=fabs(last_plane-p_iterr->first)>0.03;
01207 
01208         //if(checker !=cluster_map->rend()) if(fabs(checker->first -p_iterr->first)<0.03)nextplane=0;
01209 
01210 
01211 
01212         //are we now in a different plane and do we have something to add?
01213         if(nextplane && closest)
01214         {
01215 
01216                         //see if we didn't add the last guy and it is the only one in the plane and lies upon 
01217                         //the line made by the last added point and the current point being considered
01218                         
01219                         MSG("LongMuonFinder",Msg::kDebug)<<"notadded "<<notadded<<", planecount "<<planecount<<", lastplanecount "<<lastplanecount<<", last_added "<<(last_added!=0)<<", last_closest "<< (last_closest!=0)<<"\n";
01220                         if(notadded&&planecount==1&&lastplanecount==1&&last_added&&last_closest)
01221                         {
01222                                 double slope = (last_added->t-closest->t)/(last_added->z-closest->z);
01223                                 double p  = slope * (last_closest->z-closest->z)+closest->t;
01224                                 
01225                                 MSG("LongMuonFinder",Msg::kDebug)<<"expecting not added point at "<<p<<" and it is at "<<last_closest->t<<" with z "<<last_closest->z<<"\n";
01226                                                 
01227                                 if(fabs(p-last_closest->t) < 0.05 +(c->front_slope==0?0.05:0))
01228                                 {
01229                                         c->insert(last_closest->t,last_closest->z,last_closest->e,last_closest->id);
01230                                         last_closest=0;
01231                                         notadded=0;
01232                                 }
01233                         }
01234 
01235 
01236 
01237                 double prev_dt=0;
01238                 if(c->entries>1)prev_dt=fabs(c->t[0]-c->t[1]);
01239 //              if(closest &&( (prev_dt<0.001?0.025:prev_dt)*2+0.05 +0.02*(fabs(last_z-closest->z)/0.035)> last_d || first))
01240                 if(closest &&( /*(prev_dt<0.025?0.025:prev_dt)*2*/ 0.025 + 0.01*(fabs(last_z-closest->z)/0.035) + (lastplanecount==1?0.002:0)> last_d || first))
01241                         {
01242                                 c->insert(closest->t,closest->z,closest->e,closest->id);
01243                                 first=0;
01244                                 last_t=closest->t;
01245                         last_z=closest->z;
01246                         MSG("LongMuonFinder",Msg::kDebug)<<"adding "<<last_z<<" "<<last_t<<"\n";
01247                         notadded=0;
01248                         last_added=closest;
01249                         lastplanecount=planecount;
01250                         }else{
01251                                 MSG("LongMuonFinder",Msg::kDebug)<<"failed projection to z "<<closest->z<<" t "<<closest->t<<" gap "<<last_d<<"\n";
01252                                 
01253                                 if(closest_proj)
01254                                         MSG("LongMuonFinder",Msg::kDebug)<<"checking other extrap prev_dt "<<prev_dt<<" last_z "<<last_z<<" proj_z "<<closest_proj->z<<" lpc "<<lastplanecount<<" "<<(prev_dt<0.025?0.025:prev_dt)*2 +0.2+0.002*(fabs(last_z-closest_proj->z)/0.035) + (lastplanecount==1?0.02:0)<<" > "<< last_d_proj<<"\n";
01255                                 if(closest_proj &&( /*(prev_dt<0.025?0.025:prev_dt)*2*/ 0.025 +0.01*(fabs(last_z-closest_proj->z)/0.035) + (lastplanecount==1?0.002:0)> last_d_proj || first))
01256                                 {
01257                                         c->insert(closest_proj->t,closest_proj->z,closest_proj->e,closest_proj->id);
01258                                         first=0;
01259                                         last_t=closest_proj->t;
01260                                 last_z=closest_proj->z;
01261                                         MSG("LongMuonFinder",Msg::kDebug)<<"adding "<<last_z<<" "<<last_t<<"\n";
01262                                         notadded=0;
01263                                 last_added=closest_proj;
01264                                 lastplanecount=planecount;
01265                                 }else{
01266                                 
01267                                         lastplanecount=planecount;
01268                                         last_closest=closest;
01269                                         closest=0;
01270                                         closest_proj=0;
01271                                         notadded=1;
01272                                 }
01273                         }
01274                 
01275 
01276 
01277                         
01278         
01279                 closest=0;
01280                 last_d=100000;
01281                         closest_proj=0;
01282                         last_d_proj=100000;
01283                 }
01284                 
01285                 if(nextplane)
01286                 {
01287                         MSG("LongMuonFinder",Msg::kDebug)<<"\n\n-------\nnext plane\n";
01288                 MSG("LongMuonFinder",Msg::kDebug)<<"last "<<last_plane<<" cur "<<p_iterr->first<<"\n";
01289         }
01290     
01291                 
01292 
01293                 if(nextplane)
01294                 {
01295                         planecount=0;
01296                 
01297                         proj_t=last_t;
01298                         if(c->entries>3) proj_t = c->interpolate(p_iterr->first);
01299                         
01300                         lin_proj_t = c->front_slope*c->start_z+c->front_offset;
01301                         if(lin_proj_t==0)lin_proj_t=last_t;
01302                         MSG("LongMuonFinder",Msg::kDebug)<<"quad proj "<<proj_t<<" front proj "<<lin_proj_t<<"\n";
01303                         
01304 
01305                         //punch-through....
01306                         //check next plane --- only  1 hit?
01307                         //extrap to next plane
01308                         //is the hit in the next plane sufficiently close? <1 strip?
01309                         //set a flag....
01310                         //then the closest in this current plane is the closest to the punchthrough line
01311 
01312                         punch_through=0;
01313                         punch_through_t=0;
01314                         punch_through_z=0;
01315                         closest_punch_through=10000;    
01316                         std::map<double, std::map<double, int> >::reverse_iterator z_it = p_iterr;
01317                         z_it++;
01318                         std::map<double, int>::iterator t_it;
01319                         
01320                         //max number of planes to look forward
01321                         int maxcheck=4;
01322                         
01323 //                      printf("next plane\n");
01324                         while(maxcheck>0 && z_it !=cluster_map->rend())
01325                         {
01326 
01327                         
01328                         punch_through=0;
01329                         punch_through_t=0;
01330                         punch_through_z=0;
01331                         double last_z = z_it->first;
01332                         int cnt=0;
01333                 //      printf("\n\nlast_z %f\n",last_z);
01334                         for(;z_it!=cluster_map->rend();z_it++)
01335                         {
01336                                 if(fabs(last_z - z_it->first)>0.03)break;
01337                                 for(t_it=z_it->second.begin();t_it!=z_it->second.end();t_it++)
01338                                 {
01339                                         cnt++;
01340                                         if(cnt>1)break;
01341                                         punch_through_z=z_it->first;
01342                                         punch_through_t=t_it->first;
01343                 //                      printf("%f %f   \n",punch_through_z,punch_through_t);
01344                                 }
01345                                 if(cnt>1)break;
01346                         }
01347                 //      printf("cnt %d\n",cnt);
01348                         if(cnt==1)
01349                         {
01350                                 //is that point close to where we expect it?
01351                                 double exp  = c->front_slope*(punch_through_z) +c->front_offset ;
01352                                 
01353                 //              printf("fs %f fo %f cz %f ct %f pz %f pt %f\n",c->front_slope,c->front_offset,c->start_z,c->start_t,punch_through_z,punch_through_t);
01354                                 
01355                 //              printf("ext %f  diff %f\n",exp,fabs(exp-punch_through_t));
01356                                 if(fabs(exp-punch_through_t)<0.02)
01357                                 {
01358                                         punch_through=1;
01359                                         //save the expected point in the plane currently under consideration!
01360                                         double curz=p_iterr->first;
01361                                         double dz = (c->start_z - punch_through_z);
01362                                         double dt = (c->start_t - punch_through_t);
01363                                         
01364                         //              printf("pt %f pz %f dz %f dt %f sl %f\n",punch_through_t,punch_through_z,dz,dt,dt/dz);
01365                                         
01366 
01367                                         punch_through_t=dz ? dt/dz*(curz-punch_through_z)+punch_through_t : punch_through_t;
01368                                         punch_through_z=curz;
01369                                 }
01370                         }
01371 
01372                 //      if(punch_through)printf("punch through found at z %f t%f\n",punch_through_z,punch_through_t);
01373                         if(punch_through)break;
01374                         
01375                         z_it++;
01376                         maxcheck--;                     
01377                         }
01378 
01379                 }       
01380                 
01381 
01382          for(s_iterr=p_iterr->second.begin();s_iterr!=p_iterr->second.end(); s_iterr++)
01383          {
01384                                 Managed::ManagedCluster * clus = cl->GetCluster(s_iterr->second);
01385 
01386 
01387                                 if(fabs(clus->t-lin_proj_t)<last_d)
01388                                 {
01389                                         last_d=fabs(clus->t-lin_proj_t);
01390                                         closest=clus;
01391                                 }
01392                                 
01393                                 if(fabs(clus->t-proj_t)<last_d_proj)
01394                                 {
01395                                         last_d_proj=fabs(clus->t-proj_t);
01396                                         closest_proj=clus;
01397                                 }
01398                                 
01399                                 //also consider a hit at the same dt
01400                                 if(last_d < 0.05 && fabs(clus->t-last_t)<last_d)
01401                                 {
01402                                         last_d=fabs(clus->t-last_t);
01403                                         closest=clus;
01404                                 }       
01405                                 
01406                                 //also consider if the last hit was not added... if we added that last signle hit, would this hit then fit better?
01407                                 if(notadded && lastplanecount==1)
01408                                 {
01409                                         double slope = (last_notadded_t-last_t)/(last_notadded_z-last_z);
01410                                         double proj = last_t - (last_z-clus->z) * slope;
01411                                         MSG("LongMuonFinder",Msg::kDebug)<<"notadded fit projects to z "<<clus->z<<" t "<<proj<<"\n";
01412                                         MSG("LongMuonFinder",Msg::kDebug)<<"slope from past points z t "<<last_notadded_z<<" "<<last_notadded_t<<"    "<<last_z<<" "<<last_t<<"\n";
01413                                         if(fabs(proj - clus->t)<0.15)
01414                                         {
01415                                                 closest=clus;
01416                                                 last_d =  fabs(proj - clus->t);
01417                                                 MSG("LongMuonFinder",Msg::kDebug)<<"closest!\n";
01418                                         } 
01419                                 }
01420                                 
01421                                 //give priority to the punchthrough
01422                                 if(punch_through)
01423                                 {
01424                                         if(fabs(clus->t-punch_through_t)<closest_punch_through)
01425                                         {
01426                                                 closest_punch_through=fabs(clus->t-punch_through_t);
01427                                                 closest=clus;
01428                                         }       
01429                                 }
01430                                 
01431                                 MSG("LongMuonFinder",Msg::kDebug)<<clus->z<<" "<< clus->t<< " dist "<<fabs(clus->t-lin_proj_t)<<" lastdist "<<last_d<<" lastdistproj "<<last_d_proj<<"\n";
01432                                                         
01433                 planecount++;
01434          }
01435                 last_plane=p_iterr->first;
01436                 
01437                 //pick up the last hit details... only use if you know there is only one hit in this plane!
01438                 for(s_iterr=p_iterr->second.begin();s_iterr!=p_iterr->second.end(); s_iterr++)
01439         {
01440                         Managed::ManagedCluster * clus = cl->GetCluster(s_iterr->second);
01441                         last_notadded_z=clus->z;
01442                         last_notadded_t=clus->t;
01443                 }
01444         }
01445 
01446 
01447                 double prev_dt=0;
01448                 if(c->entries>1)prev_dt=fabs(c->t[0]-c->t[1]);
01449 //              if(closest &&( (prev_dt<0.001?0.025:prev_dt)*2+0.05 +0.02*(fabs(last_z-closest->z)/0.035)> last_d || first))
01450                 if(closest &&( /*(prev_dt<0.025?0.025:prev_dt)*2*/ 0.025 + 0.01*(fabs(last_z-closest->z)/0.035) + (lastplanecount==1?0.002:0)> last_d || first))
01451                         {
01452                                 c->insert(closest->t,closest->z,closest->e,closest->id);
01453                                 first=0;
01454                                 last_t=closest->t;
01455                         last_z=closest->z;
01456                         MSG("LongMuonFinder",Msg::kDebug)<<"adding "<<last_z<<" "<<last_t<<"\n";
01457                         notadded=0;
01458                         last_added=closest;
01459                         lastplanecount=planecount;
01460                         }else{
01461                                 if(closest)MSG("LongMuonFinder",Msg::kDebug)<<"failed projection to z "<<closest->z<<" t "<<closest->t<<" gap is "<<last_d<<"\n";
01462                                 else MSG("LongMuonFinder",Msg::kDebug)<<"no closest\n";
01463                                 
01464                                 if(closest_proj)
01465                                         MSG("LongMuonFinder",Msg::kDebug)<<"checking other extrap prev_dt "<<prev_dt<<" last_z "<<last_z<<" proj_z "<<closest_proj->z<<" lpc "<<lastplanecount<<" "<<(prev_dt<0.025?0.025:prev_dt)*2 +0.2+0.002*(fabs(last_z-closest_proj->z)/0.035) + (lastplanecount==1?0.02:0)<<" > "<< last_d_proj<<"\n";
01466                                 if(closest_proj &&( /*(prev_dt<0.025?0.025:prev_dt)*2*/ 0.025 +0.01*(fabs(last_z-closest_proj->z)/0.035) + (lastplanecount==1?0.002:0)> last_d_proj || first))
01467                                 {
01468                                         c->insert(closest_proj->t,closest_proj->z,closest_proj->e,closest_proj->id);
01469                                         first=0;
01470                                         last_t=closest_proj->t;
01471                                 last_z=closest_proj->z;
01472                                         MSG("LongMuonFinder",Msg::kDebug)<<"adding "<<last_z<<" "<<last_t<<"\n";
01473                                         notadded=0;
01474                                 last_added=closest_proj;
01475                                 lastplanecount=planecount;
01476                                 }else{
01477                                 
01478                                         lastplanecount=planecount;
01479                                         last_closest=closest;
01480                                         closest=0;
01481                                         closest_proj=0;
01482                                         notadded=1;
01483                                 }
01484                         }
01485                 
01486 
01487 
01488                         
01489         
01490    
01491                         
01492                         
01493 return c;                       
01494 
01495 
01496 }

int LongMuonFinder::FoundSingleViewLongMuon (  )  [inline]

Definition at line 25 of file LongMuonFinder.h.

References single_view_long_muon.

00025 {return single_view_long_muon;};

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

Definition at line 22 of file LongMuonFinder.h.

References chain_u, chain_v, and foundparticle.

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

Particle3D* LongMuonFinder::GetParticle3D (  )  [inline]

Definition at line 21 of file LongMuonFinder.h.

References foundparticle, and foundparticle3d.

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

int LongMuonFinder::IsPartiallyInstrumented ( double  t,
double  z,
int  view 
) [private]

Definition at line 38 of file LongMuonFinder.cxx.

References detector, Detector::kFar, and Detector::kNear.

Referenced by CountInPartiallyInstrumentedRegion().

00039 {
00040         if(detector==Detector::kFar)return 0;
00041 
00042         if(detector==Detector::kNear)
00043         {
00044                 if(view==2)
00045                 {
00046                         if(z<0 ||z > 7)return 1;
00047                         if(t<-0.2 || t>2.3)return 2;
00048                         return 0;
00049                 
00050                 }else if(view==3)
00051                 {
00052                         if(z<0 ||z > 7)return 1;
00053                         if(t<-2.3 || t>0.2)return 2;
00054                         return 0;                       
00055                 }else
00056 //              printf("unknown view in LongMuonFinder::IsPartiallyInstrumented\n");
00057                 abort();
00058         }
00059 
00060 
00061 //      printf("unknown detector in LongMuonFinder::IsPartiallyInstrumented\n");
00062         abort();
00063 }

void LongMuonFinder::MakeParticle3D (  )  [private]

Definition at line 709 of file LongMuonFinder.cxx.

References MuELoss::a, Particle3D::add_to_back(), point::chain, chain_u, chain_v, point::chainhit, Chain::cluster_id, cluster_manager, point::e, Chain::e, MuELoss::e, Chain::entries, Particle3D::entries, Particle3D::finalize(), foundparticle3d, Managed::ClusterSaver::GetCluster(), Managed::ClusterManager::GetClusterSaver(), Chain::myId, Chain::parentChain, Particle3D::particletype, Chain::particletype, pointgreaterlmf(), point::rms_t, Managed::ManagedCluster::rms_t, point::t, Chain::t, point::view, point::z, and Chain::z.

Referenced by FindLongMuon().

00710 {
00711 
00712         //verify that the chains in both views have sufficient overlap
00713         
00714         
00715         //make the 3d particle from these chains
00716         if(foundparticle3d){delete foundparticle3d;foundparticle3d=0;};
00717         foundparticle3d= new Particle3D();
00718 
00719 
00720         
00721         std::vector<point> points;
00722         
00723         double start =0;
00724         double end  =1000;
00725         
00726         double endu=0;
00727         double endv=0;
00728                 
00729         
00730         
00731         
00732         int type=0;
00733 
00734 
00735                 Chain *c = chain_u;
00736                 
00737                 if(c->particletype)
00738                         type = c->particletype;
00739                 
00740                 for(int j=0;j<c->entries;j++)
00741                 {
00742                         if(c->parentChain==-1 && j==0)
00743                                 start = start < c->z[j]? c->z[j] : start;
00744                         if( j == c->entries-1)
00745                         {
00746                                 end = end < c->z[j] ? end : c->z[j];
00747                         }
00748                         endu=endu<c->z[j]?c->z[j]:endu;
00749 
00750                         
00751                         point a;
00752                         a.z=c->z[j];
00753                         a.t=c->t[j];
00754                         a.e=c->e[j];
00755                         a.chain=c->myId;
00756                         a.chainhit=j;
00757                         a.view=2;
00758                         
00759                         Managed::ManagedCluster * clu = cluster_manager->GetClusterSaver()->GetCluster(c->cluster_id[j]);
00760                         if(clu)
00761                         {
00762                                 a.rms_t=clu->rms_t;
00763                                 points.push_back(a);
00764                         }
00765                 }
00766         
00767         
00768 
00769                 c = chain_v;
00770                 
00771                 if(c->particletype)
00772                         type = c->particletype;
00773                 for(int j=0;j<c->entries;j++)
00774                 {
00775                         if(c->parentChain==-1 && j==0)
00776                                 start = start < c->z[j]? c->z[j] : start;
00777                         if( j == c->entries-1)
00778                         {
00779                                 end = end < c->z[j] ? end : c->z[j];
00780                         }
00781                         endv=endv<c->z[j]?c->z[j]:endv;
00782                         
00783 
00784                         point a;
00785                         a.z=c->z[j];
00786                         a.t=c->t[j];
00787                         a.e=c->e[j];
00788                         a.chain=c->myId;
00789                         a.chainhit=j;
00790                         a.view=3;
00791                         
00792                         Managed::ManagedCluster * clu =cluster_manager->GetClusterSaver()->GetCluster(c->cluster_id[j]);
00793                         if(clu)
00794                         {
00795                                 a.rms_t=clu->rms_t;
00796                                 points.push_back(a);
00797                         }
00798                 }
00799         
00800         
00801         
00802         //we don't want the particle to extrapolate too far.... ...but now go all the way to the end
00803         double stopend = endu<endv?endv:endu;
00804 
00805                 
00806         sort(points.begin(), points.end(),pointgreaterlmf);
00807         
00808         for(unsigned int i=0;i<points.size();i++)
00809         {
00810 //              if( start - points[i].z > 0.045 )continue;
00811 //              if( points[i].z - end > 0.045)continue; //within 1 planes, we want the end of the chain
00812         
00813 //      printf("point %f %f %f %d\n",points[i].z,points[i].t,points[i].e,points[i].view);
00814         
00815                 if(points[i].z>stopend)continue;
00816         
00817                 int myview = points[i].view;
00818                 int lower=-1;
00819                 int upper=-1;
00820                 for(int j=i-1;j>-1;j--)
00821                 {
00822                         if(points[j].view!=myview)
00823                         {
00824                                 lower=j;
00825                                 break;
00826                         }
00827                 }
00828                 for(unsigned int j=i+1;j<points.size();j++)
00829                 {
00830                         if(points[j].view!=myview)
00831                         {
00832                                 upper=j;
00833                                 break;
00834                         }
00835                 }       
00836                 
00837                 
00838                 
00839                 double u = points[i].view==2 ? points[i].t : 0;
00840                 double v = points[i].view==3 ? points[i].t : 0;
00841                 double e = points[i].e;
00842                 double z = points[i].z;
00843                 int view = points[i].view;
00844                 
00845                 double rms_t = points[i].rms_t;
00846                 
00847 
00848                 
00849                 if(lower>-1 && upper > -1 )// good we can extrapolate!
00850                 {
00851                         double s = (points[upper].t - points[lower].t) /  ( points[upper].z - points[lower].z);
00852                         
00853                         double t = s * (points[i].z-points[lower].z) + points[lower].t;
00854                         
00855                         u = myview == 2 ? u : t;
00856                         v = myview == 3 ? v : t;
00857                         
00858 
00859                 }else if(lower>-1 && upper < 0)  //just user the closest other view t value
00860                 {
00861                 
00862 
00863                         //do we have another lower value?
00864                         int lower2=-1;
00865                         for(int j=lower-1;j>-1;j--)
00866                         {
00867                                 if(points[j].view!=myview && fabs(points[lower].z-points[j].z)>0.04)
00868                                 {
00869                                         lower2=j;
00870                                         break;
00871                                 }
00872                         }
00873                         
00874                         if(lower2>-1 && fabs( points[lower].z - points[lower2].z) >0)
00875                         {
00876                                 double s = (points[lower].t - points[lower2].t) /  ( points[lower].z - points[lower2].z);
00877                         
00878                                 double t = s * (points[i].z-points[lower2].z) + points[lower2].t;
00879                                 u = myview == 2 ? u : t;
00880                                 v = myview == 3 ? v : t;
00881                         
00882                         }else{
00883                                 u = myview == 2 ? u : points[lower].t;
00884                                 v = myview == 3 ? v : points[lower].t;
00885                         }
00886                 }
00887                 else if(upper>-1 && lower < 0)   //just user the closest other view t value
00888                 {
00889                 
00890 
00891                 
00892                         //do we have another upper value?
00893                         int upper2=-1;
00894                         for(unsigned int j=upper+1;j<points.size();j++)
00895                         {
00896                                 if(points[j].view!=myview && fabs(points[upper].z-points[j].z)>0.04)
00897                                 {
00898                                         upper2=j;
00899                                         break;
00900                                 }
00901                         }
00902                         
00903                         if(upper2>-1 && fabs( points[upper2].z - points[upper].z)>0)
00904                         {
00905                                 double s = (points[upper2].t - points[upper].t) /  ( points[upper2].z - points[upper].z);
00906                         
00907                                 double t = s * (points[i].z-points[upper].z) + points[upper].t;
00908                                 u = myview == 2 ? u : t;
00909                                 v = myview == 3 ? v : t;
00910                         
00911                         }else{
00912                                 u = myview == 2 ? u : points[upper].t;
00913                                 v = myview == 3 ? v : points[upper].t;
00914                         }
00915 
00916 
00917 
00918                         //lets use the vertex!
00919 
00920 /*
00921                         if(points[upper].view != points[i].view)
00922                                 upper=upper2;
00923                                 
00924                         if(points[upper].view == points[i].view)
00925                         {
00926 */
00927 //dont yet have a vertex!
00928 /*                              double vz =vtx_z;
00929                                 double vt = 0;
00930                                 vt = points[upper].view == 2 ? vtx_u : vt;
00931                                 vt = points[upper].view == 3 ? vtx_v : vt;
00932                         
00933                                 double t =vt;
00934                                 
00935                                 //if the point at z is not the same as z vtx, we need to extrapolate 
00936                                 if(fabs ( points[upper].z - vz)>0.001)
00937                                 {                       
00938                         
00939                                         double s = (points[upper].t - vt) /  ( points[upper].z - vz);
00940                                         t = s * (points[i].z-vz) + vt;                  
00941                         
00942                                 }
00943                                 
00944                 //              printf("---view %d u %f v %f t %f\n",myview,u,v,t);
00945                 //              printf("---vtx z %f t%f upper z %f t %f\n",vz,vt,points[upper].z,points[upper].t);
00946                         
00947 
00948                                 u = myview == 2 ? u : t;
00949                                 v = myview == 3 ? v : t;
00950 */
00951                                 u = myview == 2 ? u : points[upper].t;
00952                                 v = myview == 3 ? v : points[upper].t;
00953 
00954 //                      }       
00955                                 
00956 
00957                 }
00958                 else if(upper==-1 && lower==-1) //we have an empty view!!!
00959                 {
00960 
00961                         u = myview == 2 ? u : 0;
00962                         v = myview == 3 ? v : 0;
00963                 }
00964                 
00965                 foundparticle3d->add_to_back(u,v,z,e,points[i].chain,points[i].chainhit,view,rms_t);
00966                 
00967         //      printf("adding 3d point to particle %f %f %f --- %f\n",u,v,z,e);
00968         }
00969         
00970         
00971         if(type)
00972                 foundparticle3d->particletype=(Particle3D::EParticle3DType)type;
00973                 
00974         foundparticle3d->finalize();
00975         
00976         if(foundparticle3d->entries<1){delete foundparticle3d;foundparticle3d=0;}
00977 
00978 
00979 }

void LongMuonFinder::MergeChainClusters ( Chain ch,
std::vector< int > *  clusters 
) [private]

Definition at line 488 of file LongMuonFinder.cxx.

References Chain::cluster_id, cluster_manager, Managed::ManagedCluster::e, Chain::e, Managed::ClusterManager::GetCluster(), Managed::ManagedCluster::id, Chain::insert(), Managed::ClusterManager::MergeClusters(), Managed::ManagedCluster::t, Managed::ManagedCluster::z, and Chain::z.

Referenced by AbsorbMuonClusters().

00489 {
00490         if(clusters->size()<1)return;
00491 
00492 /*
00493         //find the appropriate spot in the chain for these clusters... 
00494         Managed::ManagedCluster *current=0;
00495         
00496         unsigned int idx;
00497         for(idx=0;idx<ch->z.size();idx++)
00498         {
00499                 if(fabs(ch->z[idx]-(*clusters)[0]->z)<0.01)
00500                 {
00501                         current = cluster_manager->GetCluster(ch->cluster_id[idx]);
00502                         break;
00503                 }
00504         }
00505                 
00506         int add_to_chain=0;     
00507         if(!current)
00508         {       
00509                 add_to_chain=1;
00510                 current=(*clusters)[0];
00511         }
00512 */
00513 
00514         
00515         int mergeid=(*clusters)[0];
00516         for(unsigned int i=1;i<(*clusters).size();i++)
00517         {
00518                 mergeid=cluster_manager->MergeClusters(mergeid,(*clusters)[i]);
00519         }
00520 
00521 
00522         Managed::ManagedCluster *current = cluster_manager->GetCluster(mergeid);
00523 
00524         //do we already have a cluster for that plane?
00525         int add_to_chain=1;     
00526         int idx=-1;
00527         for(idx=0;idx<(int)ch->z.size();idx++)
00528         {
00529                 if(fabs(ch->z[idx]-current->z)<0.01)
00530                 {
00531                         add_to_chain=0;
00532                         break;
00533                 }
00534         }
00535 
00536 
00537         if(add_to_chain)
00538         {
00539                 ch->insert(current->t, current->z, current->e, current->id);
00540         }else{
00541                 ch->e[idx]=current->e;
00542                 //for our purposes... we just want to account from brehm like energy... we dont want to move the point of the muon track hit
00543         //      ch->t[idx]=current->t;
00544         //      ch->z[idx]=current->z;
00545                 ch->cluster_id[idx]=current->id;
00546         }
00547 
00548 }

void LongMuonFinder::RemoveNonMuonEnergy ( Chain c,
int  view,
double  past_z 
) [private]

Definition at line 396 of file LongMuonFinder.cxx.

References Chain::cluster_id, cluster_manager, Managed::ManagedCluster::e, Chain::e, MuELoss::e, Managed::ClusterManager::GetCluster(), Managed::ClusterManager::GetClusterMap(), Managed::ManagedCluster::id, Msg::kDebug, and MSG.

Referenced by FindLongMuon().

00397 {
00398 
00399 
00400         if(!cluster_manager)return;
00401         
00402         std::map<double, std::map<double, int> > * cluster_map = cluster_manager->GetClusterMap(view);
00403     std::map<double, std::map<double, int> >::iterator p_iterr;
00404     std::map<double, int >::iterator s_iterr;
00405 
00406         //first determine the average energy deposition
00407         double sum_e=0;
00408         int cnt_e=0;
00409         for(p_iterr=cluster_map->begin();p_iterr!=cluster_map->end(); p_iterr++)
00410     {   
00411         if(p_iterr->first-past_z<0.01)continue;
00412     
00413         for(s_iterr=p_iterr->second.begin();s_iterr!=p_iterr->second.end(); s_iterr++)
00414         {
00415                 double e =cluster_manager->GetCluster(s_iterr->second)->e;
00416                 if(e<0.1)continue;
00417                 cnt_e++;
00418                 sum_e+=e;
00419         }
00420     }
00421 
00422         if(cnt_e<1)return;
00423         double avg_e=sum_e/(double)cnt_e;
00424         
00425 
00426         for(p_iterr=cluster_map->begin();p_iterr!=cluster_map->end(); p_iterr++)
00427     {   
00428         if(p_iterr->first-past_z>0.01)return;
00429         for(s_iterr=p_iterr->second.begin();s_iterr!=p_iterr->second.end(); s_iterr++)
00430         {
00431                 Managed::ManagedCluster *c = cluster_manager->GetCluster(s_iterr->second);
00432                 if(c->e < avg_e * 1.2)continue;
00433                 for(unsigned int i=0;i<ch->cluster_id.size();i++)
00434                 {
00435                         if(ch->cluster_id[i]==c->id)
00436                         {       
00437                                 MSG("LongMuonFinder",Msg::kDebug)<<"changing cluster energy in chain from "<<ch->e[i]<<" to "<<avg_e<<"\n";
00438                                 ch->e[i]=avg_e;
00439 
00440                                 break;
00441                         }
00442                 }
00443                         
00444                 }
00445         }
00446         
00447 
00448 }

void LongMuonFinder::Reset (  )  [private]

Definition at line 203 of file LongMuonFinder.cxx.

References chain_u, chain_v, foundparticle, foundparticle3d, and single_view_long_muon.

Referenced by FindLongMuon(), and LongMuonFinder().

00204 {
00205         if(foundparticle3d){delete foundparticle3d; foundparticle3d=0;}
00206         foundparticle=0;
00207         if(chain_u){delete chain_u; chain_u=0;}
00208         if(chain_v){delete chain_v; chain_v=0;}
00209         single_view_long_muon=0;
00210 }


Member Data Documentation

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

Definition at line 48 of file LongMuonFinder.h.

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

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

Definition at line 49 of file LongMuonFinder.h.

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

Definition at line 47 of file LongMuonFinder.h.

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

Definition at line 46 of file LongMuonFinder.h.

Referenced by DumpParticle3D(), FindLongMuon(), GetParticle3D(), MakeParticle3D(), and Reset().

Definition at line 56 of file LongMuonFinder.h.

Referenced by FindLongMuon(), FoundSingleViewLongMuon(), and Reset().


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

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1