#include <LongMuonFinder.h>
|
|
Definition at line 194 of file LongMuonFinder.cxx. References detector, and Reset().
|
|
|
Definition at line 200 of file LongMuonFinder.cxx. 00201 {}
|
|
||||||||||||||||
|
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 }
|
|
||||||||||||||||
|
Definition at line 365 of file LongMuonFinder.cxx. References Chain::end_z, 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 }
|
|
||||||||||||||||
|
Definition at line 633 of file LongMuonFinder.cxx. References Chain::cluster_id, cluster_manager, Chain::e, Chain::entries, Managed::ClusterManager::GetCluster(), Managed::ManagedCluster::hitplane, 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 }
|
|
||||||||||||
|
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 }
|
|
||||||||||||
|
Definition at line 65 of file LongMuonFinder.cxx. References detector, IsPartiallyInstrumented(), 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 }
|
|
|
Definition at line 1499 of file LongMuonFinder.cxx. References Particle3D::avg_rms_t, Particle3D::calibrated_energy, Particle3D::chain, Particle3D::chainhit, Particle3D::e, Particle3D::emfit_a, Particle3D::emfit_b, Particle3D::emfit_chisq, Particle3D::emfit_e0, Particle3D::emfit_ndf, Particle3D::emfit_prob, Particle3D::end_u, Particle3D::end_v, Particle3D::end_z, Particle3D::entries, MSG, Particle3D::muonfrac, Particle3D::particletype, Particle3D::rms_t, s(), Particle3D::shared, Particle3D::start_u, Particle3D::start_v, Particle3D::start_z, Particle3D::types, Particle3D::u, 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 }
|
|
|
Definition at line 552 of file LongMuonFinder.cxx. References cluster_manager, Managed::ManagedCluster::e, Managed::ClusterManager::GetCluster(), Managed::ClusterManager::GetClusterMap(), 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 }
|
|
|
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(), MakeParticle3D(), MSG, 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 }
|
|
||||||||||||
|
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, Chain::start_t, Chain::start_z, Chain::t, Managed::ManagedCluster::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 }
|
|
|
Definition at line 25 of file LongMuonFinder.h. 00025 {return single_view_long_muon;};
|
|
|
Definition at line 22 of file LongMuonFinder.h. 00022 {if(foundparticle){if(view==2)return chain_u;if(view==3)return chain_v;}return 0;};
|
|
|
Definition at line 21 of file LongMuonFinder.h. 00021 {if(foundparticle)return foundparticle3d;return 0;};
|
|
||||||||||||||||
|
Definition at line 38 of file LongMuonFinder.cxx. References detector. 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 }
|
|
|
Definition at line 709 of file LongMuonFinder.cxx. References Particle3D::add_to_back(), point::chain, point::chainhit, Chain::cluster_id, cluster_manager, Managed::ManagedCluster::e, Chain::e, point::e, Particle3D::entries, Chain::entries, Particle3D::finalize(), foundparticle3d, Managed::ClusterSaver::GetCluster(), Managed::ClusterManager::GetClusterSaver(), Chain::myId, Chain::parentChain, Particle3D::particletype, Chain::particletype, Managed::ManagedCluster::rms_t, point::rms_t, s(), Managed::ManagedCluster::t, Chain::t, point::t, Managed::ManagedCluster::view, point::view, Managed::ManagedCluster::z, 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 }
|
|
||||||||||||
|
Definition at line 488 of file LongMuonFinder.cxx. References Chain::cluster_id, cluster_manager, Chain::e, Managed::ManagedCluster::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 }
|
|
||||||||||||||||
|
Definition at line 396 of file LongMuonFinder.cxx. References Chain::cluster_id, cluster_manager, Chain::e, Managed::ManagedCluster::e, Managed::ClusterManager::GetCluster(), Managed::ClusterManager::GetClusterMap(), Managed::ManagedCluster::id, 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 }
|
|
|
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 }
|
|
|
Definition at line 18 of file LongMuonFinder.cxx. Referenced by FindLongMuon(), and Reset(). |
|
|
Definition at line 19 of file LongMuonFinder.cxx. Referenced by FindLongMuon(), and Reset(). |
|
|
Definition at line 53 of file LongMuonFinder.h. Referenced by AbsorbMuonClusters(), CheckChainQuality(), FindIsolationZ(), FindLongMuon(), MakeParticle3D(), MergeChainClusters(), and RemoveNonMuonEnergy(). |
|
|
Definition at line 65 of file LongMuonFinder.h. Referenced by CountInPartiallyInstrumentedRegion(), IsPartiallyInstrumented(), and LongMuonFinder(). |
|
|
Definition at line 47 of file LongMuonFinder.h. Referenced by FindLongMuon(), and Reset(). |
|
|
Definition at line 17 of file LongMuonFinder.cxx. Referenced by FindLongMuon(), MakeParticle3D(), and Reset(). |
|
|
Definition at line 56 of file LongMuonFinder.h. Referenced by FindLongMuon(), and Reset(). |
1.3.9.1