Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages
NueAna / ParticlePID / ParticleFinder

PrimaryShowerFinder.cxx

Go to the documentation of this file.
00001 #include "NueAna/ParticlePID/ParticleFinder/PrimaryShowerFinder.h"
00002 #include <map>
00003 #include <math.h>
00004 #include <algorithm>
00005 
00006 #include "MessageService/MsgService.h"
00007 
00008 #include "NueAna/ParticlePID/ParticleFinder/Managed/ManagedCluster.h"
00009 #include <sstream>
00010 #include "NueAna/ParticlePID/ParticleFinder/ShwFit.h"
00011 
00012 
00013 ClassImp(PrimaryShowerFinder)
00014 CVSID("$Id: PrimaryShowerFinder.cxx,v 1.3 2009/09/11 05:00:40 gmieg Exp $");
00015 
00016 
00017 TH2D * PrimaryShowerFinder::houghmapU =0;
00018 TH2D * PrimaryShowerFinder::houghmapV =0;
00019 
00020 TH2D * PrimaryShowerFinder::intU=0;
00021 TH2D * PrimaryShowerFinder::intV=0;     
00022 
00023 ChainHelper *PrimaryShowerFinder::chu=0;
00024 ChainHelper *PrimaryShowerFinder::chv=0;
00025 
00026 Particle3D * PrimaryShowerFinder::foundparticle3d=0;
00027 
00028 Chain * PrimaryShowerFinder::chain_u=0;
00029 Chain * PrimaryShowerFinder::chain_v=0;
00030                 
00031 
00032 struct spoint
00033 {
00034         int view;
00035         double z;
00036         double t;
00037         double e;
00038         int chain;
00039         int chainhit;
00040         double rms_t;
00041 };
00042         
00043 bool spointgreaterlmf(spoint p1, spoint p2){return p1.z < p2.z;}
00044 
00045 PrimaryShowerFinder::PrimaryShowerFinder()
00046 {
00047 
00048         Reset();
00049 }
00050 
00051 PrimaryShowerFinder::~PrimaryShowerFinder()
00052 {
00053         Reset();
00054 }
00055 
00056 void PrimaryShowerFinder::Reset(int reset_vertex)
00057 {
00058         ran=0;
00059         if(foundparticle3d){delete foundparticle3d;foundparticle3d=0;}
00060         foundparticle=0;
00061         chain_u=0; 
00062         chain_v=0;
00063         single_view_long_shower=0;
00064 
00065 
00066 //      if(houghmapU)delete houghmapU;
00067 //              houghmapU=0;
00068 //      if(houghmapV)delete houghmapV;
00069 //              houghmapV=0;
00070 
00071 
00072         if(houghmapU)houghmapU->Reset();
00073         if(houghmapV)houghmapV->Reset();
00074 
00075         houghlinesU.clear();
00076         houghlinesV.clear();
00077         houghlineMatch.clear();
00078 
00079         if(intU)delete intU;
00080                 intU=0;
00081         if(intV)delete intV;
00082                 intV=0;
00083 
00084         if(reset_vertex)
00085         {
00086                 vtx_u=0;
00087                 vtx_v=0;
00088                 vtx_z=0;
00089                 foundvertex=0;
00090                 foundvertexU=0;
00091                 foundvertexV=0;
00092         }
00093 }
00094 
00095 void PrimaryShowerFinder::MakeChains(int view)
00096 {
00097 
00098         ChainHelper *ch = view==2?chu:chv;
00099         if(!ch)return;
00100         
00101 
00102         
00103 
00104         
00105         std::vector<HoughLine> * hl = view==2?&houghlinesU:&houghlinesV;
00106         
00107         if(hl->size()<1)
00108         {
00109                 //does the view have only one cluster? if so... make that a chain!
00110                 std::map<double, std::map<double, int> > * cluster_map = cluster_manager->GetClusterMap(view);
00111                 std::map<double, std::map<double, int> >::reverse_iterator p_iterr;
00112         std::map<double, int >::iterator s_iterr;
00113 
00114                 Managed::ManagedCluster *cluster=0;
00115                 int cluster_count=0;
00116                 for(p_iterr=cluster_map->rbegin();p_iterr!=cluster_map->rend(); p_iterr++)
00117         {       
00118         
00119                 for(s_iterr=p_iterr->second.begin();s_iterr!=p_iterr->second.end(); s_iterr++)
00120                         {
00121                                 cluster = cluster_manager->GetCluster(s_iterr->second);
00122                                 cluster_count++;
00123                         }
00124         }
00125         if(cluster_count==1 && cluster)
00126         {
00127                         int nc = ch->NewChain();
00128                         Chain *c = ch->GetChain(nc);
00129                         
00130                 
00131                 //      int newid = cluster_manager->SaveCluster(cluster->id,cluster->e,2);
00132                 //      cluster = cluster_manager->GetCluster(newid);
00133                         
00134                         if(cluster)
00135                         {
00136                                 c->insert(cluster->t, cluster->z, cluster->e, cluster->id);
00137                         }       
00138                         MSG("PrimaryShowerFinder",Msg::kDebug)<<"single cluster chain made in view "<<view<<"\n";
00139         
00140         }
00141                 
00142         
00143                 return;
00144         }
00145         
00146         
00147         //make a copy....
00148         
00149         std::vector<HoughLine> hlcopy = *hl;
00150 
00151         hl=&hlcopy;
00152 
00153         Managed::ClusterManager *mycm=cluster_manager;
00154         
00155 
00156         //probably vertex 
00157         //to check if chains after the first are vertex pointing or other chain pointing
00158         double prob_vtx_z=0;
00159         double prob_vtx_t=0;
00160 
00161         if(view ==2 && foundvertexU || foundvertex)
00162         {
00163                 prob_vtx_z=vtx_z;
00164                 prob_vtx_t=vtx_u;
00165                 
00166                 //printf("USING VERTEX U %f %f\n",prob_vtx_z,prob_vtx_t);
00167         }
00168 
00169         if(view ==3 && foundvertexV || foundvertex)
00170         {
00171                 prob_vtx_z=vtx_z;
00172                 prob_vtx_t=vtx_u;
00173                 
00174                 //printf("USING VERTEX V %f %f\n",prob_vtx_z,prob_vtx_t);
00175         }       
00176 
00178         //uncomment to do a temp        
00179         Managed::ClusterManager cm=*cluster_manager;
00180         cm.SetClusterSaver(0);
00181         cm.SetHitManager(0);
00182         
00183         mycm=&cm;
00184         mycm->ClearInUse();
00185 
00187 
00188         
00189         std::vector<int>foundChains;//store the found chains here as well... so we can split up hough lines as needed
00190         
00191         
00192         //need to load any existing chains from the long muon finder!!!!
00194         for(unsigned int i=0;i<ch->finished.size();i++)
00195         {
00196                 foundChains.push_back(ch->finished[i].myId);
00197         
00198         }
00199         
00200         
00201         
00203         
00204 
00205         //do a full reload...
00206         for(unsigned int k=0;k<hl->size();k++)
00207         {
00208 //              printf("%d\n",k);
00209                 (*hl)[k].ResetHits(0);
00210                 LoadCloseHits(&(*hl)[k],mycm,view,0.0412*1.5);
00211 //              printf("%f %f   %f %f\n",(*hl)[k].start_z,(*hl)[k].start_t,(*hl)[k].end_z,(*hl)[k].end_t);      
00212         }
00213         std::vector<HoughLine> hlcopy2 = *hl;
00214                 
00215 //      
00216 //      printf("!!!!!!!!!!!!!!!!!!!!!!!\nlooking for chains\n");
00217         for(int ccc=0;ccc<10;ccc++)
00218         {
00219         //      printf("try %d with %d houghlines\n",ccc,hl->size());
00220                 for(unsigned int k=0;k<hl->size();k++)
00221                 {
00222         
00223         //              printf("%d\n",k);
00224                         (*hl)[k].ResetHits(1); //keep the bounds
00225                         
00226                 //      printf("reset hl %d\n",k);
00227                         LoadCloseHits(&(*hl)[k],mycm,view,0.0412*1.5,1);
00228         //              printf("current %d,%f %f %f %f\n",(&(*hl)[k])->ncluster,(&(*hl)[k])->start_z,(&(*hl)[k])->start_t,(&(*hl)[k])->end_z,(&(*hl)[k])->end_t);
00229                         
00230                         
00231 //                      printf("loaded %d hits\n",(*hl)[k].ncluster);
00232                         
00233                 
00234                         //does this hough line intersect any found chains?  if so, we should split it!
00235                         for(unsigned int c=0;c<foundChains.size();c++)
00236                         {
00237 //                              printf("c %d of %d\n",c,foundChains.size());
00238 //                              printf("current %d,%f %f %f %f\n",(&(*hl)[k])->ncluster,(&(*hl)[k])->start_z,(&(*hl)[k])->start_t,(&(*hl)[k])->end_z,(&(*hl)[k])->end_t);
00239                                 
00240                                 Chain *chain = ch->GetChain(foundChains[c]);
00241                                 if(!chain)continue;
00242                                 if(chain->entries<2)continue;
00243                                 
00244                                 HoughLine *lnew = SplitHoughLine(chain,&(*hl)[k],mycm);
00245 //                              printf("split  %d,%f %f %f %f\n",(&(*hl)[k])->ncluster,(&(*hl)[k])->start_z,(&(*hl)[k])->start_t,(&(*hl)[k])->end_z,(&(*hl)[k])->end_t);
00246                                 
00247                                                         
00248                                 if(lnew && lnew->ncluster>0)
00249                                 {
00251                                         hl->push_back(*lnew);
00252 //                                      printf("%d %f %f %f %f\n",lnew->ncluster,lnew->start_z,lnew->start_t,lnew->end_z,lnew->end_t);
00253                                 }
00254                                 
00255                                 if(lnew)delete lnew;
00256                                 //delete it if it is empty!
00257                                 if((*hl)[k].ncluster<1)
00258                                 {
00259 //                                      printf("hl size %d hl at position %d is empty...\n",hl->size(),k);
00260                                         
00261                                 //      std::vector<HoughLine> hh;
00262                                 //      for(unsigned int aw=0;aw<hl->size();aw++)
00263                                 //              if(aw!=k)hh.push_back((*hl)[aw]);
00264                                         
00265                                 //      *hl=hh;
00266                                         
00267                                         hl->erase(hl->begin()+k,hl->begin()+k+1);
00268                                         k--;
00269 //                                      printf("%d %d\n",k,hl->size());
00270                                         break;//restart the check now that this hl is split
00271                                 }
00272 
00273                         }
00274 //                      printf("e %d\n",k);
00275                 }
00276         
00277                 CleanHoughLines(0,hl);
00278                 
00279                 //so we know where to attach secondary chains
00280                 std::vector<int>closest_chain;
00281                 std::vector<double>closest_chain_z;
00282                 for(unsigned int k=0;k<hl->size();k++)
00283                 {
00284                         closest_chain.push_back(-1);
00285                         closest_chain_z.push_back(0);
00286                 }
00287 
00288 
00289                 std::vector<int>toskip;
00290 
00291                 if( ( view ==2 && foundvertexU == 1 ) || ( view==3 && foundvertexV==1)  || foundvertex) //if we don't have a vertex, we don't have any chains... so consider all of them and skip this part....
00292                 for(unsigned int k=0;k<hl->size();k++)
00293                 {
00294                         HoughLine * l = &(*hl)[k];
00295 
00296                         int keep=0;
00297                         //check for chain pointing      
00298                         double lastdist=10000;
00299                         for(unsigned int i=0;i<foundChains.size();i++)
00300                         {       
00301                                 Chain *c = ch->GetChain(foundChains[i]);
00302                                 if(!c)continue;
00303                                 if(c->entries<2)continue; //can come from when a chain is split at the first cluster
00304                 
00305 
00306 
00307                                 //measure the distance of intersection
00308                                 double vz = vtx_z;
00309                                 double vt = view==2? vtx_u:vtx_v;
00310                                 
00311                                 
00312                                 double c_start_z, c_start_t, c_end_z, c_end_t =0;
00313                                 
00314                                 if( fabs((c->start_z-vz)*(c->start_z-vz)+(c->start_t-vt)*(c->start_t-vt)) < fabs((c->end_z-vz)*(c->end_z-vz)+(c->end_t-vt)*(c->end_t-vt)) )
00315                                 {
00316                                         c_start_z = c->start_z;
00317                                         c_start_t = c->start_t;
00318                                         c_end_z = c->end_z;
00319                                         c_end_t = c->end_t;     
00320                                 }else{
00321                                         c_start_z = c->end_z;
00322                                         c_start_t = c->end_t;
00323                                         c_end_z = c->start_z;
00324                                         c_end_t = c->start_t;                           
00325                                 }
00326                                 
00327                                  if(!(c_end_z-c_start_z))continue;
00328 
00329                                 double l_start_z, l_start_t, l_end_z, l_end_t =0;
00330                                 
00331                                 if( fabs((l->start_z-c_start_z)*(l->start_z-c_start_z)+(l->start_t-c_start_t)*(l->start_t-c_start_t)) < fabs((l->end_z-c_start_z)*(l->end_z-c_start_z)+(l->end_t-c_start_t)*(l->end_t-c_start_t)) )
00332                                 {
00333                                         l_start_z = l->start_z;
00334                                         l_start_t = l->start_t;
00335                                         l_end_z = l->end_z;
00336                                         l_end_t = l->end_t;     
00337                                 }else{
00338                                         l_start_z = l->end_z;
00339                                         l_start_t = l->end_t;
00340                                         l_end_z = l->start_z;
00341                                         l_end_t = l->start_t;                           
00342                                 }
00343                                 
00344 
00345                                 //find intersection
00346                                 
00347                                 double lslope = (l_end_t-l_start_t) / (l_end_z-l_start_z) ;
00348                                 double loffset = l_end_t - l_end_z*lslope;
00349                                 
00350                                 double cslope = (c_end_t-c_start_t) / (c_end_z-c_start_z) ;
00351                                 double coffset = c_end_t - c_end_z*cslope;
00352                                 
00353                                 double int_z = cslope-lslope ? (loffset - coffset ) / ( cslope - lslope):l_start_z;
00354                                 double int_t = cslope * int_z + coffset;
00355                                                                                                 
00356                                 //require intersection to be behind the start of the chain
00357 
00358                                 //if(fabs((c_start_z-vz)*(c_start_z-vz)+(c_start_t-vt)*(c_start_t-vt)) > fabs((int_z-vz)*(int_z-vz)+(int_t-vt)*(int_t-vt)) )continue;
00359                                 
00360                                 //calculate distance along chain path... where start is 0 and going towards end increases
00361                                 //if its <0  its bad....
00362                                 double r_c_end_z=c_end_z-c_start_z;
00363                                 double r_c_end_t=c_end_t-c_start_t;
00364                                 double r_l_end_z=l_end_z-c_start_z;
00365                                 double r_l_end_t=l_end_t-c_start_t;
00366                                 double r_l_start_z=l_start_z-c_start_z;
00367                                 double r_l_start_t=l_start_t-c_start_t;
00368                                 
00369                                 double theta = atan(r_c_end_t/r_c_end_z);                               
00370                                 double rr_int_z = (int_z-c_start_z) * cos(theta) + (int_t-c_start_t) * sin(theta);                              
00371                                 double rr_int_t = (int_t-c_start_t) * cos(theta) - (int_z-c_start_z) * sin(theta);                                      
00372                                 if(rr_int_z<0.03)continue;// points to far fowards of chain
00373                                 
00374                                 
00375                                 
00376                                 
00377                                 double ltheta = atan( ( l->end_t-l->start_t) / (l->end_z-l->start_z) );                         
00378                                 double rr_l_int_z = (int_z-l->start_z) * cos(ltheta) + (int_t-l->start_t) * sin(ltheta);                                
00379         
00380                                 //printf("in hl coords, int is at z %f\n",rr_l_int_z);
00381         
00382                                 //require the intersection to be forwards of the hl!
00383                                 if(rr_l_int_z>0)continue;                               
00384                                 
00385                                                 
00386                                 
00387                                 //its the intersection past the end of the chain?  if so... check the relative angle between the hl and the chain to keep it reasonable
00388                                 double rr_c_end_z = r_c_end_z * cos(theta) + r_c_end_t * sin(theta);                                    
00389                                 //double rr_c_end_t =  -r_c_end_z * sin(theta) + r_c_end_t * cos(theta);                                
00390 
00391 
00392                                 //rotate the coordinates so the chain is along the x axis
00393 
00394                                 double rr_l_end_z = r_l_end_z * cos(theta) + r_l_end_t * sin(theta);
00395                                 double rr_l_end_t =  -r_l_end_z * sin(theta) + r_l_end_t * cos(theta); 
00396                                 double rr_l_start_z = r_l_start_z * cos(theta) + r_l_start_t * sin(theta);
00397                                 double rr_l_start_t =  -r_l_start_z * sin(theta) + r_l_start_t * cos(theta);                            
00398         
00399                                 
00400                                 if(rr_int_z > rr_c_end_z)
00401                                 {
00402                                         //count the angle from the intersection!
00403                                         double dy = rr_l_start_t - rr_int_t;
00404                                         double dx = rr_l_start_z - rr_int_z;
00405                                         
00406                                         double intangle = dx ? atan(dy/dx) : 0;
00407                                         if(intangle>3.141592/6 || dx==0)continue; //to steep!
00408                                 
00409                                 }
00410                                 
00411                                 //is it the closest?  
00412                                 
00413                                 double dist = sqrt((l_start_z - int_z)*(l_start_z - int_z)+(l_start_t - int_t)*(l_start_t - int_t));
00414                                 
00415                                 //if the intersection is past the end of the chain... need to include the distance from the intersection to the end of the chain
00416                                 if(fabs((c_end_z-vz)*(c_end_z-vz)+(c_end_t-vt)*(c_end_t-vt)) < fabs((int_z-vz)*(int_z-vz)+(int_t-vt)*(int_t-vt)) )
00417                                 {
00418                                         dist+= sqrt((c_end_z - int_z)*(c_end_z - int_z)+(c_end_t - int_t)*(c_end_t - int_t));
00419                                 }
00420 
00421                                 //is the closest ?
00422                                 if(dist > lastdist)continue;
00423                                 //printf("closer to chain %f %f %f %f at z %f\n",c_start_z,c_start_t,c_end_z,c_end_t,int_z);
00424                                 //if so, check for orientation before recording
00425                         
00426                                 if(fabs(rr_l_end_t) <= fabs(rr_l_start_t) && rr_l_start_z < rr_l_end_z ||
00427                                         fabs(rr_l_end_t) >= fabs(rr_l_start_t) && rr_l_start_z > rr_l_end_z 
00428                                  )
00429                                 {
00430                                         //its closer... but no good... we want to invalidate the last guess...
00431                                         closest_chain[k]=-1;
00432                                         closest_chain_z[k]=int_z;
00433                                         lastdist=dist;  
00434                                         keep=0;                         
00435                                         continue;
00436                                 }
00437 
00438                                 //printf("in chain coords line has start %f %f end %f %f\n",rr_l_start_z,rr_l_start_t,rr_l_end_z,rr_l_end_t);
00439                                 //printf("hl start %f %f end %f %f\n",l_start_z,l_start_t,l_end_z,l_end_t);
00440                                 //record it
00441                                 closest_chain[k]=foundChains[i];
00442                                 closest_chain_z[k]=int_z;
00443                                 lastdist=dist;
00444                                 keep=1;
00445                                 
00446                                 //printf("keeping %d it points to chain %f %f %f %f at z %f\n",k,c_start_z,c_start_t,c_end_z,c_end_t,int_z);
00447                                 
00448                         }
00449                         
00450                         
00451                         //dont have anything yet?
00452                         //see if we are vertex pointing
00453                         if(closest_chain[k]<0)
00454                         {
00455                                 double exp_vtx_t = l->GetExpectedT(vtx_z);
00456                                 if((view==2&&foundvertexU || view==3&&foundvertexV) && fabs(exp_vtx_t-(view==2?vtx_u:vtx_v))<2*0.0451) //smaller than 2 strips away?
00457                                 {
00458                                 
00459                                         //need to make sure that the vertex is not within the hl!
00460         
00461                                 
00462                                         double theta = atan( ( l->end_t-l->start_t) / (l->end_z-l->start_z) );                          
00463                                         double rr_end_z = (l->end_z-l->start_z) * cos(theta) + (l->end_t-l->start_t) * sin(theta);                              
00464                                         //double rr_end_t = (l->end_t-l->start_t) * cos(theta) - (l->end_z-l->start_z) * sin(theta);                                    
00465 
00466                                         double rr_vtx_z = (vtx_z-l->start_z) * cos(theta) + (exp_vtx_t-l->start_t) * sin(theta);                                
00467                                         //double rr_vtx_t = (exp_vtx_t-l->start_t) * cos(theta) - (vtx_z-l->start_z) * sin(theta);      
00468         
00469                                         //printf("in hl coords:  hl 0 0 to %f %f  vtx %f %f\n",rr_end_z,rr_end_t,rr_vtx_z,rr_vtx_t);
00470         
00471                                         if(rr_vtx_z< 0 || rr_vtx_z > rr_end_z)
00472                                         {                               
00473                                                 //printf("keeping %d vp estimate vtx %f %f actual %f %f\n",k,vtx_z,l->GetExpectedT(vtx_z),vtx_z,view==2?vtx_u:vtx_v);
00474                                                 keep=1; //we'll keep it
00475                                         }
00476                                 }       
00477                         }
00478                 
00479 
00480                         if(keep)
00481                         {
00482                                 //printf("keeping %d\n",k);
00483                                 continue;
00484                         }
00485                         //if we get here... we don't want this hl!
00486                         //printf("skipping %d\n",k);
00487                         toskip.push_back(k);
00488                 
00489                 }
00490 
00491 
00492 /*
00493 
00494 
00495                 std::vector<int>toskip;
00496                 //make sure that all hough line are vertex pointing or pointing to a another chain and on the correct angle
00497                 if(foundChains.size()>0)
00498                 {
00499                         for(unsigned int k=0;k<hl->size();k++)
00500                         {
00501                                 HoughLine * l = &(*hl)[k];
00502         //                                      printf("current %d,%f %f %f %f\n",(&(*hl)[k])->ncluster,(&(*hl)[k])->start_z,(&(*hl)[k])->start_t,(&(*hl)[k])->end_z,(&(*hl)[k])->end_t);
00503                                                 
00504                                                         
00505                                 int keep=0;                             
00506                                 //is it vertex pointing?
00507                                 if(foundvertex && fabs(l->GetExpectedT(prob_vtx_z)-prob_vtx_t)<2*0.0451) //smaller than 2 strips away?
00508                                 {
00509                                         printf("keeping %d vp\n",k);
00510                                         keep=1; //we'll keep it
00511                                         //we still want to see if it is pointing to another chain so we know if it needs to be a secondary
00512                                 }
00513 
00514                                 //now we need to see if it is pointing to a current chain ... and find the closest one!
00515                                 int closest_pointing=-1;
00516                                 double closest_dist=0;
00517                                 double dist=0;          
00518                                 for(unsigned int i=0;i<foundChains.size();i++)
00519                                 {
00520                                         Chain *c = ch->GetChain(foundChains[i]);
00521                                         double exp_t_start = l->GetExpectedT(c->start_z);
00522                                         double exp_t_end = l->GetExpectedT(c->end_z);
00523 
00524                         
00525                                         if(( exp_t_start < c->start_t && exp_t_end > c->end_t ) || ( exp_t_start > c->start_t && exp_t_end < c->end_t))
00526                                         {
00527                                         //              //check for orientation
00528                                         //              
00529                                         //              if(fabs(c->end_z - l->end_z) > fabs(c->end_z - l->start_z))continue;
00530                                         //              if(fabs(c->start_z - l->start_z) > fabs(c->start_z - l->end_z))continue;
00531                                                         
00532                                                         MSG("PrimaryShowerFinder",Msg::kDebug)<<"Passing orientation on "<<c->start_z<<" "<<c->start_t<<" "<<c->end_z<<" "<<c->end_t<<"\n";
00533 
00534                                                         //check the distance
00535                                                         //the way we found the intersections requires the intersection to be within the chain....
00536                                                         for(unsigned int kk=0;kk<c->z.size()-1;kk++)
00537                                                         {
00538                                                                 double exp_t = l->GetExpectedT(c->z[kk]);
00539                                                                 double exp_t_n = l->GetExpectedT(c->z[kk+1]);
00540                                                                 
00541                                                                 if(( exp_t>=c->t[kk]&&exp_t_n<=c->t[kk+1] ) || ( exp_t<=c->t[kk]&&exp_t_n>=c->t[kk+1] ) )
00542                                                                 {//intersection is bewteen these two points!
00543                                                         
00544                                                                                 //for now, put it close to the more forward hit...
00545                                                                                 //later we can actually extrapolate it
00546                                                                                 
00547                                                                                 double t = c->t[kk];
00548                                                                                 dist = fabs(t-exp_t);
00549                                                                                 closest_chain_z[k]=c->z[kk];
00550                                                                         
00551                                                                 }
00552                                                         
00553                                                         }
00554                                                         
00555                                                         
00556                                                 if(closest_pointing<0)//the first one
00557                                                 {
00558                                                         closest_pointing=foundChains[i];
00559                                                         closest_dist =  dist;           
00560                                                 }else
00561                                                 {
00562                                                         if(dist<closest_dist)
00563                                                         {
00564                                                                 closest_pointing=foundChains[i];
00565                                                                 closest_dist =  dist;                                                           
00566                                                         }
00567                                                 }
00568                                 
00569                                                 
00570 
00571 
00572                                         }
00573                                 }
00574                                 
00575                                 for(unsigned int i=0;i<foundChains.size();i++)
00576                                 {               
00577                                                                         Chain *c = ch->GetChain(foundChains[i]);        
00578                                 
00579         //printf("end chain %f %f  start hl %f %f\n",c->end_z,c->end_t,l->start_z,l->start_t);
00580 
00581                                 
00582                                 
00583                                 if(c->end_z < l->start_z)
00584                                 {
00585                                 
00586                                                         //also check to see if this chain is at the end of another chain
00587                                                         if(
00588                                                                 fabs(c->end_t - l->GetExpectedT(c->end_z))<0.1  
00589                                                                 && l->start_z - c->end_z  <0.07*3)
00590                                                         {
00591                                                         
00592                                                                 closest_chain_z[k]=c->end_z;
00593                                                                 dist = fabs(c->end_t - l->GetExpectedT(c->end_z));
00594 
00595                                                         }
00596 
00597                                                 //      printf("end chain %f %f  start hl %f %f\n",c->end_z,c->end_t,l->start_z,l->start_t);
00598                                 
00599                                 }
00600                                 
00601                                 
00602                                 
00603                                                 if(closest_pointing<0)//the first one
00604                                                 {
00605                                                         closest_pointing=foundChains[i];
00606                                                         closest_dist =  dist;           
00607                                                 }else
00608                                                 {
00609                                                         if(dist<closest_dist)
00610                                                         {
00611                                                                 closest_pointing=foundChains[i];
00612                                                                 closest_dist =  dist;                                                           
00613                                                         }
00614                                                 }
00615                                 
00616                                 
00617                                 
00618                                 }
00619                                 
00620                                 
00621                                 if(closest_pointing>-1)
00622                                 {
00623                                                 Chain *c = ch->GetChain(closest_pointing);
00624                                                 //its pointing... so now check the orientation
00625                                                 double close_z;
00626                                                 double close_t;
00627                                                 printf("prob vtx z %f t %f\n",prob_vtx_z,prob_vtx_t);
00628                                                 printf("chain start z %f t %f, end z %f t %f\n",c->start_z,c->start_t,c->end_z,c->end_t);
00629                                                 if(fabs((c->start_t-prob_vtx_t)*(c->start_t-prob_vtx_t)+(c->start_z-prob_vtx_z)*(c->start_z-prob_vtx_z))
00630                                                         <
00631                                                         fabs((c->end_t-prob_vtx_t)*(c->end_t-prob_vtx_t)+(c->end_z-prob_vtx_z)*(c->end_z-prob_vtx_z)))
00632                                                 {       
00633                                                         close_z=c->start_z;
00634                                                         close_t=c->start_t;
00635                                                 }
00636                                                 else
00637                                                 {
00638                                                         close_z=c->end_z;
00639                                                         close_t=c->end_t;
00640                                                 }
00641                                                 
00642                                                 printf("close z %f t %f\n",close_z,close_t);
00643                                                 double ch_close_z;
00644                                                 double ch_close_t;
00645                                                 double ch_far_z;
00646                                                 double ch_far_t;
00647                                                 
00648                         //                      if(fabs(l->start_z - close_z) < fabs(l->end_z - close_z))
00649                         //                      {
00650                                                 if(fabs((l->start_t-prob_vtx_t)*(l->start_t-prob_vtx_t)+(l->start_z-prob_vtx_z)*(l->start_z-prob_vtx_z))
00651                                                         <
00652                                                         fabs((l->end_t-prob_vtx_t)*(l->end_t-prob_vtx_t)+(l->end_z-prob_vtx_z)*(l->end_z-prob_vtx_z)))
00653                                                 {       
00654                         
00655                         
00656                                                         ch_close_z = l->start_z;
00657                                                         ch_close_t = l->start_t;
00658                                                         ch_far_z = l->end_z;
00659                                                         ch_far_t = l->end_t;
00660                                                 }else{
00661                                                         ch_close_z = l->end_z;
00662                                                         ch_close_t = l->end_t;
00663                                                         ch_far_z = l->start_z;
00664                                                         ch_far_t = l->start_t;
00665                                                 }
00666                                                 
00667                                                 
00668                                                 //this will cause issues with vertical chains!
00669                                                 double dst = fabs(c->interpolate(ch_close_z)-ch_close_t);
00670                                                 double det = fabs(c->interpolate(ch_far_z)-ch_far_t);
00671                                                 printf("!!! %f\n",c->interpolate(ch_far_z));
00672                                                 printf(" close z %f t %f  far z %f t %f  dst %f det %f\n",ch_close_z,ch_close_t,ch_far_z,ch_far_t,dst,det);
00673                                                 
00674                                                 if(dst<det)
00675                                                 {
00676                                                         keep=1;//we'll keep it
00677                                                         printf("keeping %d cp\n",k);
00678                                                         closest_chain[k]=closest_pointing;
00679                                                 }
00680                                                         
00681                                 
00682                                 }
00683                                 
00684                                 if(keep)continue;
00685                                 //if we get here... we don't want this hl!
00686                                 printf("skipping %d\n",k);
00687                                 toskip.push_back(k);
00688         
00689                         }
00690                 }
00691         
00692 
00693 
00694 */
00695 
00696 
00697 
00698 
00699 
00700 
00701 
00702         for(unsigned int i=0;i<hl->size();i++)
00703         {
00704                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"hl "<<i<<" "<<(*hl)[i].start_z<<" "<<(*hl)[i].start_t<<" "<<(*hl)[i].end_z<<" "<<(*hl)[i].end_t<<"\n";
00705         }
00706         
00707         
00708         
00709         std::map<double,int> orderer;
00710 //      printf("\n");
00711         for(int i=0;i<(int)hl->size();i++)
00712         {
00713 
00714                 int save=1;
00715                 for(int j=0;j<(int)toskip.size();j++)
00716                 {
00717                         if(toskip[j]==i)
00718                         {
00719                                 save =0;
00720                                 break;
00721                         }
00722                 }
00723                 if(!save)continue;
00724 
00725                 HoughLine * l = &(*hl)[i];
00726         
00727                 double weight = l->ncluster *cos(l->phi)*l->sum_e*l->sum_e;
00728 
00729                 orderer.insert(make_pair(weight,i));
00730                 
00731 //              printf("%d %f %f | %f\n",l->ncluster,l->phi,l->sum_e,weight);
00732         }               
00733 //      printf("\n");
00734 
00735         std::map<double,int>::reverse_iterator it_orderer;
00736 
00737         HoughLine *h = 0;       
00738 
00739         //int cnt=0;
00740         //for(it_orderer = orderer.rbegin();it_orderer!=orderer.rend();it_orderer++)
00741         //{
00742         
00743         it_orderer = orderer.rbegin();
00744         while(it_orderer!=orderer.rend())
00745         {
00746                 if((*hl)[it_orderer->second].ncluster<2)it_orderer++;
00747                 else break;
00748         }
00749         
00750         if(it_orderer==orderer.rend())break;
00751         
00752                 int i= it_orderer->second;
00753                 
00754                 HoughLine * l = &(*hl)[i];
00755                 h=l;
00756                 
00757                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"chose "<<l->ncluster<<" "<<l->phi<<" "<<l->sum_e<<" | "<<it_orderer->first<<"\n";
00758 
00759         //      cnt++;
00760         //      if(cnt>3)break;
00761 
00762 
00763                 int nc = ch->NewChain();
00764                 Chain *c = ch->GetChain(nc);
00765                 foundChains.push_back(nc);
00766 
00767                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"making chain from houghline "<<i<<" "<<l->start_z<<" "<<l->start_t<<" "<<l->end_z<<" "<<l->end_t<<"\n";
00768                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"STARTING WITH "<< ch->finished.size()<<" CHAINS\n";
00769         
00770                 for(unsigned int i=0;i<h->cluster_id.size();i++)
00771                 {
00772                         Managed::ManagedCluster *cl = mycm->GetCluster(h->cluster_id[i]);
00773                         
00774                         
00775                 //      int newid = mycm->SaveCluster(cl->id,cl->e,2);
00776                 //      cl = mycm->GetCluster(newid);
00777                         
00778                         if(cl)
00779                         {
00780                         
00781                                 
00782                                 
00783                                 
00784                                 c->insert(cl->t, cl->z, cl->e, cl->id);
00785                                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"("<<cl->t<<" "<< cl->z<<" "<< cl->e<<" "<< cl->id<<") s"<<cl->GetStatus()<<" v"<<cl->view<<"\n";
00786                         }
00787 
00788                         mycm->MarkUsed(cl->id);
00789                 
00790                 }
00791                 //attach the chain if needed...
00792                 
00793                 if(closest_chain[i]>-1)
00794                 {
00795 
00796                         Chain *mom = ch->GetChain(closest_chain[i]);
00797                         Chain newdaughter = ch->AttachAt(mom,c,closest_chain_z[i]);
00798                         if(newdaughter.entries>0)ch->AddFinishedChain(newdaughter);
00799                         MSG("PrimaryShowerFinder",Msg::kDebug)<<"kept by closest_chain, attempting to attach to chain "<<mom->start_z<<" "<<mom->start_t<<" "<<mom->end_z<<" "<<mom->end_t<<" at z "<<closest_chain_z[i]<<"\n";
00800                 }
00801                 
00802                 if(foundChains.size()==1 && !foundvertex)
00803                 {
00804                         prob_vtx_z=c->start_z;
00805                         prob_vtx_t=c->start_t;
00806                         
00807                         if(view==2)foundvertexU=1;
00808                         if(view==3)foundvertexV=1;
00809                         if(foundvertexU && foundvertexV) foundvertex=1;
00810                         vtx_z=prob_vtx_z;
00811                         if(view==2)vtx_u=prob_vtx_t;
00812                         if(view==3)vtx_v=prob_vtx_t;
00813                 }
00814         
00815                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"ENDING WITH "<<ch->finished.size()<<" CHAINS\n";
00816                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"done\n";
00817 
00818         //}
00819         
00820         
00821         
00822         }
00823 
00824         MSG("PrimaryShowerFinder",Msg::kDebug)<<"done making chains\n";
00825         
00826         mycm->ClearInUse();
00827 }
00828 
00829 
00830 HoughLine  * PrimaryShowerFinder::SplitHoughLine(Chain *c,HoughLine *hl, Managed::ClusterManager *mycm)
00831 {
00832         //here we want to split the hough line if the chain bisects it... event by the chain projectionion
00833         
00834         if(!mycm)mycm=cluster_manager;
00835         
00836         if(!c || !hl)return 0;
00837         
00838         //check for intersection?
00839         int intersection=0;
00840         double intersection_z=0;
00841         double intersection_t=0;
00842         
00843         
00844         
00845         //consider the various intersections:
00846         
00847 //      printf("looking for intersection chain %f %f %f %f hl %f %f %f %f\n",c->start_z,c->start_t,c->end_z,c->end_t,hl->start_z,hl->start_t,hl->end_z,hl->end_t);
00848         
00849         //cross inside
00850 /*      if(c->z.size()<1)return 0;
00851         for(unsigned int i=0;i<c->z.size()-1;i++)
00852         {
00853                 double et1 = hl->GetExpectedT(c->z[i]);
00854                 double et2 = hl->GetExpectedT(c->z[i+1]);
00855                 
00856                 double max_t = et1>et2?et1:et2;
00857                 double min_t = et1<et2?et1:et2;
00858                 
00859                 double min_c_t = c->t[i]<c->t[i+1]?c->t[i]:c->t[i+1];
00860                 double max_c_t = c->t[i]>c->t[i+1]?c->t[i]:c->t[i+1];
00861                 
00862                 printf("minct %f mint %f maxct %f maxt %f\n",min_c_t,min_t,max_c_t,max_t);
00863                 
00864                 if( (min_c_t<min_t && max_c_t>max_t) || (max_t>max_c_t && min_t < min_c_t))
00865                 //intersection
00866                 {
00867                         intersection=1;
00868                         intersection_z = (c->z[i]+c->z[i+1])/2;
00869                         intersection_t = (min_c_t+max_c_t)/2;
00870                 
00871                 }
00872         }
00873 */      
00874         if(c->z.size()<1)return 0;
00875         
00876         double last_int_z=c->z[0];
00877         double last_int_dt=c->t[0]-hl->GetExpectedT(c->z[0]);
00878         
00879         for(unsigned int i=1;i<c->z.size();i++)
00880         {
00881                 double et1 = hl->GetExpectedT(c->z[i]);
00882 
00883                 double dt=c->t[i]-et1;  
00884                         
00885 //              printf("last z %f dt %f  now z %f dt %f\n",last_int_z,last_int_dt,c->z[i],dt);
00886                 
00887                 //intersection
00888                 if( ( last_int_dt<0 && dt >0 ) || ( last_int_dt >0 && dt < 0) )
00889                 {
00890                         intersection=1;
00891 
00892 /*
00893                         intersection_z = c->z[i];
00894                         intersection_t = c->t[i];
00895 */
00896 
00897                         //find the actual intersection
00898                         
00899                         double cslope=(c->z[i]-c->z[i-1]) ? (c->t[i]-c->t[i-1])/(c->z[i]-c->z[i-1]) : 0;
00900                         double coffset=c->t[i]-cslope*c->z[i];
00901                         
00902                         double ht2 = hl->GetExpectedT(c->z[i]);
00903                         double ht1 = hl->GetExpectedT(c->z[i-1]);
00904 
00905                         double hslope = (c->z[i]-c->z[i-1]) ? (ht2-ht1)/(c->z[i]-c->z[i-1]) : 0;
00906                         double hoffset = ht2 - hslope*c->z[i];
00907                         
00908                         intersection_t = (hslope-cslope) ? (hslope*coffset - cslope*hoffset)/(hslope-cslope) : 0;
00909                         if(cslope)intersection_z = (intersection_t - coffset ) / cslope;
00910                         else intersection_z=0;
00911                         
00912 
00913                         break;
00914                 }
00915                 
00916                 last_int_z=c->z[i];
00917                 last_int_dt=dt;
00918                 
00919         }
00920         
00921         
00922         //cross with projection
00923         
00924         
00925         
00926         if(!intersection)return 0;
00927 
00928 //      printf("intersection found\n");
00929 
00930         double hl_min_z = hl->start_z<hl->end_z?hl->start_z:hl->end_z;
00931         double hl_max_z = hl->start_z>hl->end_z?hl->start_z:hl->end_z;
00932         
00933         if(intersection_z < hl_min_z || intersection_z > hl_max_z)return 0;
00934         
00935         HoughLine *newHL=new HoughLine();
00936         *newHL=*hl;
00937         newHL->ResetHits();
00938         
00939         HoughLine oldHL=*hl;
00940         oldHL.ResetHits();
00941         
00942         //otherwise, split the hough line at the intersection point....
00943         for(unsigned int i=0;i<hl->cluster_id.size();i++)
00944         {
00945                 Managed::ManagedCluster *cl = mycm->GetCluster(hl->cluster_id[i]);
00946                 
00947                 if(cl->z-intersection_z<-0.01)  //require a buffer here in case the chain has two hits in the same z plane which don't have exactly same z values
00948                 {
00949                         oldHL.AddCluster(cl);
00950 //                      printf("adding to old hl  %f %f\n",cl->z,cl->t);
00951                 }
00952                 else
00953                 {
00954                         newHL->AddCluster(cl);
00955 //                      printf("adding to new hl  %f %f\n",cl->z,cl->t);
00956                 }
00957         }
00958         
00959         
00960 
00961         
00962         //its possible that we made a new houghline which has the same hits as the old houghline, but the old line now is empty....
00963         //in such a case... just keep the old hough line
00964         //this is required to prevent infitie and useless recursion
00965         
00966         if(oldHL.ncluster<1 && newHL->ncluster>0)
00967         {
00968                 *hl=*newHL;
00969                 delete newHL;
00970                 return 0;
00971         }
00972         
00973         *hl=oldHL;      
00974         return newHL;
00975 }
00976 
00977 
00978 
00979 int PrimaryShowerFinder::FindFirstIntersection(int view, double &t, double &z)
00980 {
00981         std::vector<HoughLine> * hl;
00982         if(view==2)hl=&houghlinesU;
00983         else if(view==3)hl=&houghlinesV;
00984         else return 0;
00985 
00986         TH2D * inth = view==2?intU:intV;
00987         
00988         if(inth)delete inth;
00989         inth=0;
00990         
00991         double min_t = view==2?cluster_manager->minu:cluster_manager->minv;
00992         double max_t = view==2?cluster_manager->maxu:cluster_manager->maxv;
00993         inth = new TH2D("","",50,cluster_manager->minz,cluster_manager->maxz,50,min_t,max_t);
00994         inth->SetDirectory(0);
00995 
00996         double best_z = 10000;
00997         double best_t=0;
00998 
00999         //iterate over each line
01000         for(unsigned int i=0;i<(*hl).size();i++)
01001         {
01002                 //nested loop
01003                 for(unsigned int j=0;j<(*hl).size();j++)
01004                 {
01005                         if(i==j)continue;
01006                 
01007                         //computer intersection and save if better
01008                         HoughLine *h1 = &(*hl)[i];
01009                         HoughLine *h2 = &(*hl)[j];
01010                         
01011                         
01012                         z=
01013                         (
01014                                 h2->r/sin(h2->theta) + h2->offset_t 
01015                                 - h1->r/sin(h1->theta) - h1->offset_t 
01016                                 +(cos(h2->theta)/sin(h2->theta))*h2->offset_z 
01017                                 -(cos(h1->theta)/sin(h1->theta))*h1->offset_z
01018                         ) 
01019                                 /
01020                         (
01021                                 cos(h2->theta)/sin(h2->theta)
01022                                 -
01023                                 cos(h1->theta)/sin(h1->theta)
01024                         );
01025                         
01026                         if(z<best_z && z>0)
01027                         {
01028                                 best_z=z;
01029                                 best_t= h1->GetExpectedT(best_z);                       
01030                         }
01031                         
01032                         //require the intersection to be within 2 planes of the start/end of either houghline
01033                         if(fabs(h1->start_z-z)>0.07 && fabs(h2->start_z-z)>0.07 && fabs(h1->end_z-z)>0.07 && fabs(h2->end_z-z)>0.07)continue;
01034                         
01035                                 inth->Fill(z,h1->GetExpectedT(z),h1->ncluster+h2->ncluster);
01036         
01037                 }
01038         }
01039         
01040         inth->Draw("colz");
01041         if(best_z >=10000)return 0;
01042 
01043 
01044         t=best_t;
01045         z=best_z;
01046         return 1;
01047 }
01048 
01049 
01050 
01051 
01052 void PrimaryShowerFinder::Bundle(int view)
01053 {
01054         std::vector<HoughLine> * hl;
01055         if(view==2)hl=&houghlinesU;
01056         else if(view==3)hl=&houghlinesV;
01057         else return;
01058 
01059 
01060         std::vector<std::pair<int,int> >to_merge;
01061         
01062         //iterate over each line
01063         for(unsigned int i=0;i<(*hl).size();i++)
01064         {
01065         
01066 //                      printf("houghline with e %f clusters %d (t,z) (%f,%f) (%f,%f) chi2/ndf %f  offz %f expt %f phi %f\n", (*hl)[i].sum_e, (*hl)[i].ncluster, (*hl)[i].start_t, (*hl)[i].start_z, (*hl)[i].end_t, (*hl)[i].end_z, (*hl)[i].chi2/(*hl)[i].ncluster,(*hl)[i].offset_z,(*hl)[i].GetExpectedT((*hl)[i].offset_z),(*hl)[i].phi);
01067                         
01068                 //nested loop
01069                 for(int j=i+1;j<(int)(*hl).size();j++)
01070                 {
01071 
01072                         //if(i==j)continue;
01073                         
01074                         //did we already match this chain?
01075                         int done=0;
01076                         for(int k=0;k<(int)to_merge.size();k++)
01077                         {
01078                                 if(to_merge[k].first==j || to_merge[k].second==j)
01079                                 {
01080                                         done=1;
01081                                         break;
01082                                 }
01083                         }
01084                         if(done)continue;
01085                         
01086                 //
01087 //                      printf("\t %d ccc e %f clusters %d (t,z) (%f,%f) (%f,%f) chi2/ndf %f  offz %f expt %f phi %f\n", j,(*hl)[j].sum_e, (*hl)[j].ncluster, (*hl)[j].start_t, (*hl)[j].start_z, (*hl)[j].end_t, (*hl)[j].end_z, (*hl)[j].chi2/(*hl)[j].ncluster,(*hl)[j].offset_z,(*hl)[j].GetExpectedT((*hl)[j].offset_z),(*hl)[j].phi);
01088                 
01089                         if(fabs((*hl)[i].phi - (*hl)[j].phi)<0.05)
01090                         {
01091                                 if(fabs( (*hl)[i].GetExpectedT((*hl)[i].offset_z) - (*hl)[j].GetExpectedT((*hl)[i].offset_z)) < 0.05)
01092                                 {
01093                                         to_merge.push_back(make_pair(i,j));
01094 //                                      printf("tm %d %d %f %f %f\n",i,j,(*hl)[i].offset_z,(*hl)[i].GetExpectedT((*hl)[i].offset_z),(*hl)[j].GetExpectedT((*hl)[i].offset_z));
01095                                 
01096                                 }
01097                         
01098                         
01099                         }
01100         
01101                 }
01102         }
01103         
01104         if(to_merge.size()<1)return;
01105         
01106         std::vector<HoughLine> toKeep;
01107 
01108         //keep non merged lines
01109         for(unsigned int i=0;i<hl->size();i++)
01110         {
01111                 int keep=1;
01112                 for(unsigned int j=0;j<to_merge.size();j++)
01113                 {
01114                         if(to_merge[j].first ==(int)i || to_merge[j].second ==(int)i)
01115                         {
01116                                 keep=0;
01117                                 break;
01118                         }
01119                 }
01120                 if(keep){
01121                         toKeep.push_back((*hl)[i]);
01122 //                      printf("not merging %d\n",i);
01123                 }
01124         
01125         }
01126 
01127 
01128         for(unsigned int i=0;i<to_merge.size();i++)
01129         {
01130 /*
01131                 //simple for now... just pick one out of the group
01132                 if(i==0 || (i>0 && (to_merge[i].first != to_merge[i-1].first) ) )
01133                 {
01134 
01135                         toKeep.push_back((*hl)[to_merge[i].first]);
01136                         
01137                         printf("%d \n",to_merge[i].first);
01138                 }
01139 */
01140                 std::vector<int> mlist;
01141                 
01142                 mlist.push_back(to_merge[i].first);
01143                 for(unsigned int j=i+1;j<to_merge.size();j++)
01144                 {
01145                         if(to_merge[j].first !=to_merge[i].first)break;
01146                         mlist.push_back(to_merge[j].second);
01147                 }
01148 
01149 
01150                 //now we have a list of houghlines to merge...
01151                 //do a weighted average by energy
01152                 
01153                 double theta=0;
01154                 double r=0;
01155                 double offset_t=0;
01156                 double offset_z=0;
01157                 double sum_e=0;
01158                 
01159                 for(unsigned int j=0;j<mlist.size();j++)
01160                 {
01161                         HoughLine *h = &(*hl)[mlist[j]];
01162                         theta += h->theta*h->sum_e;
01163                         r += h->r*h->sum_e;
01164                         offset_t += h->offset_t*h->sum_e;
01165                         offset_z += h->offset_z*h->sum_e;
01166                         sum_e+=h->sum_e;
01167                 }
01168                 
01169                 if(sum_e<0.1)continue;
01170                 theta/=sum_e;
01171                 r/=sum_e;
01172                 offset_t/=sum_e;
01173                 offset_z/=sum_e;
01174                 
01175                 HoughLine l( theta,  r,  offset_t,  offset_z);
01176                 LoadCloseHits(&l,cluster_manager,view,0.0412);
01177                 toKeep.push_back(l);
01178 
01179         }
01180 
01181 
01182 
01183         
01184         (*hl)=toKeep;
01185         
01186 //      printf("# of hough lines %d\n",hl->size());
01187 
01188 }
01189 
01190 
01191 void PrimaryShowerFinder::CleanHoughLines(int view,std::vector<HoughLine>*hhh=0)
01192 {
01193 
01194         
01195         std::vector<HoughLine> * hl;
01196         if(!hhh)
01197         {
01198                 if(view==2)hl=&houghlinesU;
01199                 else if(view==3)hl=&houghlinesV;
01200                 else return;
01201         }else{
01202                 hl=hhh;
01203         }
01204         
01205         
01206         
01207         std::vector<int>todel;
01208         //remove all but first vertical hough lines in this view
01209         
01210         //first hit z is:
01211         double first_z = 1000;
01212         int keep_v=-1;
01213         for(unsigned int i=0;i<hl->size();i++)
01214         {
01215                 if(fabs((*hl)[i].phi-3.141592/2 )<0.1){
01216                         if((*hl)[i].start_z < first_z)
01217                         {
01218                                 first_z = (*hl)[i].start_z;
01219                                 keep_v=i;
01220                         }
01221                 }
01222         }                       
01223                 
01224         for(int i=0;i<(int)hl->size();i++)
01225         {
01226                 if(fabs((*hl)[i].phi-3.141592/2 )<0.1){
01227                         if(keep_v!=i)todel.push_back(i);
01228                 }
01229         }       
01230         
01231                 
01232                         
01233         //remove hough lines that have the closest distance between a pair of adjacent hits larger than some value
01234         double min_dist = 0.2;
01235         for(unsigned int i=0;i<hl->size();i++)
01236         {
01237         
01238 
01239                 double small_dist=1000;
01240                 int done=0;
01241                 for(unsigned int j=0;j<(*hl)[i].cluster_id.size();j++)
01242                 {
01243                         
01244                         for(unsigned int k=j+1;k<(*hl)[i].cluster_id.size();k++)
01245                         {
01246                                 int id1 = (*hl)[i].cluster_id[j];
01247                                 int id2 = (*hl)[i].cluster_id[k];
01248                                 
01249                                 Managed::ManagedCluster * c1 = cluster_manager->GetCluster(id1);
01250                                 Managed::ManagedCluster * c2 = cluster_manager->GetCluster(id2);
01251                                 if(!c1 || !c2)continue;
01252                                 
01253                         //      double dist = sqrt((c1->t-c2->t)*(c1->t-c2->t) + (c1->z-c2->z)*(c1->z-c2->z));
01254                                 
01255                         //      //adjust the distance based on the angle of the houghline
01256                         //      dist *=cos((*hl)[i].phi);
01257                         
01258                                 double dist = fabs(c1->z-c2->z);
01259                                 
01260                                 if(dist < small_dist)small_dist=dist;           
01261                                 
01262                 //              printf("dist %f\n",dist);               
01263                                 
01264                                 if(small_dist < min_dist)
01265                                 {
01266                                         done=1;
01267                                         break;
01268                                 }
01269                         }
01270                         if(done)break;  
01271                 }
01272                 if(done)continue;
01273                 
01274                 //if we are here.. we want to delete it
01275                 todel.push_back(i);
01276                 
01277         }
01278 
01279 
01280         //hit density
01281         for(unsigned int i=0;i<hl->size();i++)
01282         {
01283                 //double len = sqrt(((*hl)[i].end_z-(*hl)[i].start_z)*((*hl)[i].end_z-(*hl)[i].start_z) + ((*hl)[i].end_t-(*hl)[i].start_t)*((*hl)[i].end_t-(*hl)[i].start_t) );
01284                 
01285                 double len=fabs((*hl)[i].end_z-(*hl)[i].start_z);
01286                 double frac = len>0?(*hl)[i].ncluster/(len/0.07):0;
01287                 
01288                 if(frac <0.5)todel.push_back(i);
01289         }       
01290 
01291 
01292 
01293         //need to save?
01294         if(todel.size()<1)return;
01295         
01296         std::vector<HoughLine> temp;
01297         for(int i=0;i<(int)hl->size();i++)
01298         {
01299                 int save=1;
01300                 for(int j=0;j<(int)todel.size();j++)
01301                 {
01302                         if(todel[j]==i)
01303                         {
01304                                 save=0;
01305                                 break;
01306                         }
01307                 }       
01308                 if(!save)continue;
01309                 temp.push_back((*hl)[i]);
01310         }
01311 
01312         *hl=temp;
01313 
01314 
01315 }
01316 
01317 
01318 Chain * PrimaryShowerFinder::ExtractMuon(Chain *c)
01319 {
01320         if(!c)return 0;
01321         
01322         int required_consecutive_muon_hits=2;
01323         
01324         if(c->entries<required_consecutive_muon_hits)return 0;
01325         
01326         MSG("PrimaryShowerFinder",Msg::kDebug)<<"extract muon\n";
01327         
01328         //do we have a muon in this chain?
01329         
01330         double muoncount=0;
01331         double muonenergy=0;
01332         double sume=0;
01333         
01334         for(int i=0;i<c->entries;i++)
01335         {
01336                 if(c->e[i]<2.5)
01337                 {
01338                         muoncount++;
01339                         muonenergy+=c->e[i];
01340                 }
01341                 sume+=c->e[i];
01342         
01343         }
01344 
01345         double muonfrac = muoncount / c->entries;
01346 
01347         //double efrac=muonenergy/sume;
01348         
01349         MSG("PrimaryShowerFinder",Msg::kDebug)<<"mf "<<muonfrac<<"\n";
01350         if(muonfrac<0.2 ) return 0;//no muon
01351         
01352         
01353         
01354         //find where the muon exists the other particles... require n in a row muon hits
01355         
01356         int startmuon=0;
01357         int nmuon=0;
01358         
01359         for(int i=0;i<c->entries;i++)
01360         {
01361         
01362         //      printf("%d %f\n",i,c->e[i]);
01363                 startmuon=i;
01364                 if(c->e[i]<2.5)
01365                 {
01366                         nmuon++;
01367                         if(nmuon>=required_consecutive_muon_hits)break;  //require n muon hits in a row
01368                 
01369                 }else
01370                 {
01371                         nmuon=0;
01372                 }
01373         }
01374         
01375         if(nmuon<required_consecutive_muon_hits)return 0;//not enough muon like!
01376         //printf("sm %d\n",startmuon);
01377         startmuon=startmuon-nmuon+1;
01378         
01379         MSG("PrimaryShowerFinder",Msg::kDebug)<<"muon starts at "<<startmuon<<"\n";
01380         
01381         //determine average muon energy in this particle
01382         //using the muon hits, but ignoring large ones which could be muon->e shower
01383         
01384         double muone=0;
01385         muoncount=0;
01386         for(int i=startmuon;i<c->entries;i++)
01387         {
01388                 if(c->e[i]<2.5)
01389                 {
01390                         muone+=c->e[i];
01391                         muoncount++;
01392                 }
01393         }
01394         
01395         double adjmuon = muone/muoncount;
01396         
01397         MSG("PrimaryShowerFinder",Msg::kDebug)<<"avg muon e "<<adjmuon<<"\n";
01398         
01399         Chain ctemp;
01400         
01401         std::vector<double> savedE;
01402         for(int i=0;i<startmuon;i++)
01403         {
01404                 if(c->e[i]-adjmuon>0)
01405                 {
01406                         ctemp.insert(c->t[i],c->z[i],c->e[i]-adjmuon,c->cluster_id[i]);
01407                         savedE.push_back(adjmuon);
01408                 }else{
01409                         savedE.push_back(0);
01410                 }
01411                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"AAAAAAAAAAAA "<<c->t[i]<<" "<<c->z[i]<<" e "<<c->e[i]-adjmuon<<"\n";
01412         
01413         }
01414         
01415 
01416         Chain *muonchain = new Chain();
01417         muonchain->level=c->level; //so we get the parend id, etc
01418         muonchain->parentChain=c->parentChain;
01419         muonchain->children=c->children;
01420         
01421         
01422         for(int i=0;i<startmuon;i++)
01423         {
01424                 if(savedE[i]>0)
01425                         muonchain->insert(c->t[i],c->z[i],savedE[i],c->cluster_id[i]);
01426         }       
01427         
01428         for(int i=startmuon;i<c->entries;i++)
01429         {
01430                 muonchain->insert(c->t[i],c->z[i],c->e[i],c->cluster_id[i]);
01431         }               
01432 
01433 
01434 
01435         *c=ctemp;
01436         
01437         return muonchain;
01438 
01439 }
01440 
01441 
01442 
01443 
01444 
01445 
01446 int PrimaryShowerFinder::FindPrimaryShower(Managed::ClusterManager *cm)
01447 {
01448 
01449         MSG("PrimaryShowerFinder",Msg::kDebug)<<"looking for primary shower\n";
01450         Reset(0); //keep any vertex that we might have set
01451 
01452         ran=1;
01453 
01454         cluster_manager=cm;
01455         
01456         
01457         MakeHoughMap(2);
01458         MakeHoughMap(3);
01459 
01460 //DumpHoughLines(2);
01461 //DumpHoughLines(3);
01462 
01463         CleanHoughLines(2);
01464         CleanHoughLines(3);
01465         
01466         Bundle(2);
01467         Bundle(3);
01468 
01469         MakeChains(2);
01470         MakeChains(3);
01471 
01472 
01473 //      FindHoughLineMatches();
01474         
01475         
01476         //if(houghlineMatch.size()<1)
01477                 
01478         
01479         
01480         //try to find a long Shower chain in each view
01481         //chain_u = FindShowerChain(cm,2);
01482         //chain_v = FindShowerChain(cm,3);
01483         
01484         
01485 //      chain_u = MakeShowerChain(2);
01486 //      chain_v = MakeShowerChain(3);
01487 
01488 
01489         chain_u =0;
01490         chain_v=0;
01491         
01492         for(unsigned int i=0;i<chu->finished.size();i++)
01493         {
01494                 if(chu->finished[i].available)
01495                 {
01496                         chain_u=&chu->finished[i];
01497                         break;
01498                 }
01499         }
01500         
01501         for(unsigned int i=0;i<chv->finished.size();i++)
01502         {
01503                 if(chv->finished[i].available)
01504                 {
01505                         chain_v=&chv->finished[i];
01506                         break;
01507                 }
01508         }
01509 
01510         if(!chain_u || !chain_v)return 0;
01511         
01512 
01513 //make temp chains for use in the finder which represent a max path
01514         
01515 
01516 
01517         //printf("building chain\n");
01518         chu->print_finished();  
01519         std::pair< std::vector<int>, double> path;
01520         path = chu->FindMaxPath(chain_u);       
01521         chain_u=new Chain();
01522 
01523 
01524         //printf("chain made from ");
01525 //      for(int i=0 ;i <path.first.size();i++)
01526 //              printf("%d ",path.first[i]);
01527 //      printf("\n");
01528                 
01529 
01530 
01531         for(int i=0 ;i <(int)path.first.size();i++)
01532         {
01533                 Chain *c = chu->GetChain(path.first[i]);
01534                 if(!c)continue;
01535                 for(int j=0;j<c->entries;j++)
01536                 {
01537                         chain_u->insert(c->t[j],c->z[j],c->e[j],c->cluster_id[j]);
01538                         MSG("PrimaryShowerFinder",Msg::kDebug)<<"u "<<c->z[j]<<" "<<c->t[j]<<" "<<c->e[j]<<"\n";
01539                 }
01540         }
01541         std::vector<int> upath=path.first;
01542 
01543         //printf("building chain\n");
01544         chv->print_finished();
01545         path = chv->FindMaxPath(chain_v);       
01546         chain_v=new Chain();
01547         
01548         //printf("chain made from ");
01549         //      for(int i=0 ;i <path.first.size();i++)
01550         //              printf("%d ",path.first[i]);
01551         //      printf("\n");
01552                 
01553 
01554         for(int i=0 ;i <(int)path.first.size();i++)
01555         {
01556                 Chain *c = chv->GetChain(path.first[i]);
01557                 if(!c)continue;
01558 //              printf("chain %d\n",path.first[i]);
01559                 for(int j=0;j<c->entries;j++)
01560                 {
01561 //                      printf("id %d\n",c->cluster_id[j]);
01562                         chain_v->insert(c->t[j],c->z[j],c->e[j],c->cluster_id[j]);
01563                         MSG("PrimaryShowerFinder",Msg::kDebug)<<"v "<<c->z[j]<<" "<<c->t[j]<<" "<<c->e[j]<<"\n";
01564                 }
01565         }
01566         std::vector<int> vpath=path.first;              
01567         
01568 /*      
01569         printf("\n\n---------\nchain\n");
01570         for(int i=0;i<chain_u->entries;i++)
01571         {
01572                 printf("%d %f %f %f\n",chain_u->cluster_id[i],chain_u->z[i],chain_u->t[i],chain_u->e[i]);
01573         }
01574 
01575         printf("\n\n---------\nchain\n");
01576         for(int i=0;i<chain_v->entries;i++)
01577         {
01578                 printf("%d %f %f %f\n",chain_v->cluster_id[i],chain_v->z[i],chain_v->t[i],chain_v->e[i]);
01579         }       
01580 */      
01581         //muon extraction:
01583         Chain * mchainU=0;
01584         Chain * mchainV=0;
01585         mchainU=ExtractMuon(chain_u);
01586         mchainV=ExtractMuon(chain_v);
01588 /*      
01589         printf("\n\n---------\nmuon chain\n");
01590         if(mchainU)for(int i=0;i<mchainU->entries;i++)
01591         {
01592                 printf("%d %f %f %f\n",mchainU->cluster_id[i],mchainU->z[i],mchainU->t[i],mchainU->e[i]);
01593         }
01594 
01595         printf("\n\n---------\nmuon chain\n");
01596         if(mchainV)for(int i=0;i<mchainV->entries;i++)
01597         {
01598                 printf("%d %f %f %f\n",mchainV->cluster_id[i],mchainV->z[i],mchainV->t[i],mchainV->e[i]);
01599         }       
01600         
01601 
01602         printf("\n\n---------\nchain\n");
01603         for(int i=0;i<chain_u->entries;i++)
01604         {
01605                 printf("%d %f %f %f\n",chain_u->cluster_id[i],chain_u->z[i],chain_u->t[i],chain_u->e[i]);
01606         }
01607 
01608         printf("\n\n---------\nchain\n");
01609         for(int i=0;i<chain_v->entries;i++)
01610         {
01611                 printf("%d %f %f %f\n",chain_v->cluster_id[i],chain_v->z[i],chain_v->t[i],chain_v->e[i]);
01612         }       
01613         */
01614         
01615         //double ustart=chain_u->start_z;
01616         //double vstart=chain_v->start_z;
01617         //double uend=chain_u->end_z;
01618         //double vend=chain_v->end_z;
01619 //      ExpandShowerChain(chain_u,2,vstart,vend);
01620 //      ExpandShowerChain(chain_v,3,ustart,uend);
01621         
01622         
01623         int qual_u=0;
01624         int qual_v=0;
01625         
01626         if(chain_u)
01627         {
01628                 chain_u->Recalc();
01629                 qual_u=CheckChainQuality(chain_u);
01630         }
01631         if(chain_v)
01632         {
01633                 chain_v->Recalc();
01634                 qual_v=CheckChainQuality(chain_v);
01635         }
01636 
01637         if(qual_u&&!qual_v)
01638         {
01639                 if(chain_u->end_z - chain_u->start_z > 2)single_view_long_shower=1;
01640         }else if(qual_v&&!qual_u)
01641         {
01642                 if(chain_v->end_z - chain_v->start_z > 2)single_view_long_shower=1;
01643         }
01644         
01645         if(!qual_u || !qual_v)return 0;
01646 
01647 
01648 
01649         //should do something smarter... like try other chains...
01650         
01651         //if we have a quality chain in each view... then we will be making a particle
01652 /*      if(qual_u && qual_v)
01653         {
01654         
01655                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"qqqqq "<<qual_u <<" "<< qual_v<<"\n";
01656                 
01657         
01658                 //look in both views to find at least 3u and 3v consecutive planes with only 1 strip... that is where we say the Shower exits the rest of the shower/partices/etc
01659 
01660                 int qual=CheckChainOverlap(chain_u, chain_v);
01661                 if(qual)
01662                 {
01663                 
01664                 ClearFrontVertex(chain_u, chain_v);
01665 */
01666 
01667                 
01668                 
01669         
01670         
01671                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"making particle\n";
01672                 MakeParticle3D();
01673                 if(foundparticle3d)
01674                 {
01675                         foundparticle=1;
01676                 
01677                         foundparticle3d->particletype=Particle3D::electron;
01678                 
01679                         MSG("PrimaryShowerFinder",Msg::kDebug)<<"particle made\n";
01680                         DumpParticle3D();
01681                         
01682 //                      chain_u->available=0;
01683 //                      chain_v->available=0;
01684 
01685                         for(int i=0 ;i <(int)upath.size();i++)
01686                         {
01687                                 Chain *c = chu->GetChain(upath[i]);
01688                                 
01689                                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"u chain "<<upath[i]<<"\n";
01690                                 
01691                                 c->available=0;
01692                         }
01693                         for(int i=0 ;i <(int)vpath.size();i++)
01694                         {
01695                                 Chain *c = chv->GetChain(vpath[i]);
01696                                 
01697                                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"v chain "<<vpath[i]<<"\n";
01698                                 
01699                                 c->available=0;
01700                         }                       
01701                         
01702                         //do a fit on the particle3d to see if it is like an electron..
01703                         Particle3D *p=foundparticle3d;          
01704                         ShwFit f;
01705                         
01706                         double startz=0;
01707                         
01708                         startz=p->z[0];
01709                         
01710                         for(int i=0;i<p->entries;i++)
01711                         {
01712                                 //double planewidth=0.035;
01713                                 
01714                 //              f.Insert(p->e[i],(int)((p->z[i]-startz)/planewidth));
01715                         
01716                 //              cout <<"inserting "<<p->e[i]<<" at plane "<<(int)((p->z[i]-startz)/planewidth)<<endl;
01717                 
01718                 //for now, do the simple thing and assume each hit is in the next plane...
01719                         
01720                                 //                      f.Insert(p->e[i],i+1);
01721                         
01722                         
01723                                 f.Insert3d(p->e[i],p->u[i],p->v[i],p->z[i]);
01724                         
01725                                 //cout <<"inserting "<<p->e[i]<<" at plane "<<i+1<<endl;
01726                 
01727                         
01728                         }
01729                         
01730                         //f.Fit();
01731                 
01732 
01733                         f.Fit3d(1);
01734 
01735 
01736 
01737 
01738                         MSG("PrimaryShowerFinder",Msg::kDebug)<<"EM FIT! " << f.par_a << " " << f.par_b << " +- "<<f.par_b_err<<" "<<f.par_e0<<endl;
01739                 
01740                         p->emfit_a=f.par_a;
01741                         p->emfit_b=f.par_b;
01742                         p->emfit_e0=f.par_e0;
01743                         p->emfit_a_err=f.par_a_err;
01744                         p->emfit_b_err=f.par_b_err;
01745                         p->emfit_e0_err=f.par_e0_err;
01746                         p->emfit_prob=f.prob;
01747                         p->calibrated_energy=f.par_e0;
01748                         p->emfit_chisq=f.chisq;
01749                         p->emfit_ndf=f.ndf;
01750 
01751                         p->pred_e_a=f.pred_e_a; 
01752                         p->pred_g_a=f.pred_g_a;
01753                         p->pred_b=f.pred_b;
01754                         p->pred_e0=f.pred_e0;
01755                         p->pred_e_chisq=f.pred_e_chisq;
01756                         p->pred_e_ndf=f.pred_e_ndf;
01757                         p->pred_g_chisq=f.pred_g_chisq;
01758                         p->pred_g_ndf=f.pred_g_ndf;                             
01759                         
01760                         p->pre_over=f.pre_over;         
01761                         p->pre_under=f.pre_under;               
01762                         p->post_over=f.post_over;               
01763                         p->post_under=f.post_under;             
01764 
01765                         p->pp_chisq=f.pp_chisq;
01766                         p->pp_ndf=f.pp_ndf;
01767                         p->pp_igood=f.pp_igood;
01768                         p->pp_p=f.pp_p;
01769         
01770                         p->cmp_chisq = f.cmp_chisq;
01771                         p->cmp_ndf = f.cmp_ndf;
01772                         p->peakdiff = f.peakdiff;       
01773                         
01774                         int redochain=0; 
01775                         //do we need to redo the chain because we are shifting the start past the first hits of the chain?
01776                         if(f.zshift)
01777                         {
01778                                 if(f.zshift < 0)redochain=1;
01779                                 p->start_z -=f.zshift;
01780                                 
01781                                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"zshift "<<f.zshift<<"\n";
01782                         }
01783                         
01784                         
01785                         if(!f.conv)p->particletype=Particle3D::other;
01786                         else{
01787                         
01788                         
01789                         
01790                         
01791                         //redo the chain if needed
01792                         if(redochain)
01793                         {
01794 
01795                                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"REDOING CHAINS\n";
01796                         
01797                                 Chain temp_u;
01798                                 for(int i=0;i<chain_u->entries;i++)
01799                                 {
01800                                         if(chain_u->z[i]<p->start_z){
01801                                                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"throwing away chain hit with z "<<chain_u->z[i]<<"\n";
01802                                                 continue;
01803                                         }
01804                                         temp_u.insert(chain_u->t[i],chain_u->z[i],chain_u->e[i],chain_u->cluster_id[i]);
01805                                 }
01806                                 temp_u.myId=chain_u->myId;
01807                                 *chain_u=temp_u;                        
01808 
01809                                 Chain temp_v;
01810                                 for(int i=0;i<chain_v->entries;i++)
01811                                 {
01812                                         if(chain_v->z[i]<p->start_z){
01813                                                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"throwing away chain hit with z "<<chain_v->z[i]<<"\n";
01814                                                 continue;
01815                                         }
01816                                         temp_v.insert(chain_v->t[i],chain_v->z[i],chain_v->e[i],chain_v->cluster_id[i]);
01817                                 }
01818                                 temp_v.myId=chain_v->myId;              
01819                                 *chain_v=temp_v;
01820                                 chain_u->available=0;
01821                                 chain_v->available=0;
01822 
01823                                 delete foundparticle3d; 
01824                                 foundparticle=0;        
01825                                 foundparticle3d=0;
01826                                         
01827                                 //remake particle
01828                                 MakeParticle3D();
01829 
01830                                 if(!foundparticle3d)return 0;
01831 
01832                                 if(foundparticle3d)
01833                                 {
01834                                 foundparticle=1;        
01835                                 p=foundparticle3d;      
01836                                 
01837                                 f.Reset();
01838                                 
01839                                 double startz=0;
01840                         
01841                                 startz=p->z[0];
01842                         
01843                                 for(int i=0;i<p->entries;i++)
01844                                 {
01845                                         //double planewidth=0.035;
01846                                         f.Insert3d(p->e[i],p->u[i],p->v[i],p->z[i]);
01847                                 }
01848                         
01849 
01850                                 
01851                                 foundparticle3d->particletype=Particle3D::electron;
01852                                 f.Fit3d();
01853                         
01854                                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"EM FIT! " << f.par_a << " " << f.par_b << " +- "<<f.par_b_err<<" "<<f.par_e0<<endl;
01855                 
01856                                 p->emfit_a=f.par_a;
01857                                 p->emfit_b=f.par_b;
01858                                 p->emfit_e0=f.par_e0;
01859                                 p->emfit_a_err=f.par_a_err;
01860                                 p->emfit_b_err=f.par_b_err;
01861                                 p->emfit_e0_err=f.par_e0_err;
01862                                 p->emfit_prob=f.prob;
01863                                 p->calibrated_energy=f.par_e0;
01864                                 p->emfit_chisq=f.chisq;
01865                                 p->emfit_ndf=f.ndf;
01866 
01867                                 p->pred_e_a=f.pred_e_a; 
01868                                 p->pred_g_a=f.pred_g_a;
01869                                 p->pred_b=f.pred_b;
01870                                 p->pred_e0=f.pred_e0;
01871                                 p->pred_e_chisq=f.pred_e_chisq;
01872                                 p->pred_e_ndf=f.pred_e_ndf;
01873                                 p->pred_g_chisq=f.pred_g_chisq;
01874                                 p->pred_g_ndf=f.pred_g_ndf;                             
01875                         
01876                                 p->pre_over=f.pre_over;         
01877                                 p->pre_under=f.pre_under;               
01878                                 p->post_over=f.post_over;               
01879                                 p->post_under=f.post_under;     
01880                                 
01881                                 p->pp_chisq=f.pp_chisq;
01882                                 p->pp_ndf=f.pp_ndf;
01883                                 p->pp_igood=f.pp_igood;
01884                                 p->pp_p=f.pp_p;
01885         
01886                                 p->cmp_chisq = f.cmp_chisq;
01887                                 p->cmp_ndf = f.cmp_ndf;
01888                                 p->peakdiff = f.peakdiff;       
01889                                                                 
01890                         
01891                                 if(!f.conv)p->particletype=Particle3D::other;
01892                                 if(f.zshift)
01893                                 {
01894                                         p->start_z -=f.zshift;
01895                                         MSG("PrimaryShowerFinder",Msg::kDebug)<<"zshift "<<f.zshift<<"\n";
01896                                 }
01897                         
01898                                 }
01899                         
01900                         }
01901                         
01902                         
01903                 }       
01904 
01905                         
01906                 //save clusters involved...
01907 
01908                 //don't save until we know that we want it!
01909 
01910                 for(int i=0;i<chain_u->entries;i++)
01911                 {
01912                         int newid = cluster_manager->SaveCluster(chain_u->cluster_id[i],chain_u->e[i],2);
01913                         if(!newid)continue;//if we have an issue saving, such as there is no energy...
01914                         Managed::ManagedCluster * newcluster = cluster_manager->GetCluster(newid);
01915                         if(fabs(chain_u->e[i]-newcluster->e)>0.001)
01916                                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"saving with different e "<<chain_u->e[i]<<" "<<newcluster->e<<"\n";
01917                         chain_u->e[i]=newcluster->e;
01918                         chain_u->cluster_id[i]=newid;
01919                 }
01920 
01921 
01922 
01923                 for(int i=0;i<chain_v->entries;i++)
01924                 {
01925                         int newid = cluster_manager->SaveCluster(chain_v->cluster_id[i],chain_v->e[i],2);
01926                         if(!newid)continue;//if we have an issue saving, such as there is no energy...
01927                         Managed::ManagedCluster * newcluster = cluster_manager->GetCluster(newid);
01928                         
01929                         //printf("saving cluster id %d e %f newcluster ok ? %d\n",chain_v->cluster_id[i],chain_v->e[i],newcluster!=0);
01930                         //if(newcluster)printf("new cluster id %d  actual %d\n",newid,newcluster->id);
01931                         if(fabs(chain_v->e[i]-newcluster->e)>0.001)
01932                                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"saving with different e "<<chain_v->e[i]<<" "<<newcluster->e<<"\n";
01933                         chain_v->e[i]=newcluster->e;
01934                         chain_v->cluster_id[i]=newid;
01935                 }               
01936 
01937                         
01938                 //save chains involved!
01939                 chain_u->available=0;
01940                 chain_v->available=0;
01941 
01942                 //printf("psf wants to add %d %d\n",chain_u->myId,chain_v->myId);
01943                 chu->AddFinishedChain(*chain_u);
01944                 chv->AddFinishedChain(*chain_v);
01945                 
01946 
01947                 
01948                 //did we have muon chains?
01949                 if(mchainU)chu->AddFinishedChain(*mchainU);     
01950                 if(mchainV)chv->AddFinishedChain(*mchainV);             
01951                         
01952                 //we made the electron chain temporary....
01953                 //now remove the actual chains that were used to make the electron and muon extracted chains
01954                 
01955                 
01956                 for(unsigned int i=0;i<upath.size();i++)
01957                 {
01958                         chu->DeleteChain(upath[i]);
01959                 }
01960  
01961                 for(unsigned int i=0;i<vpath.size();i++)
01962                 {
01963                         chv->DeleteChain(vpath[i]);
01964                 }
01965                                 
01966 
01967                 
01970                                                 //remake particle
01971  
01972                                 delete foundparticle3d; 
01973                                 foundparticle=0;        
01974                                 foundparticle3d=0;
01975                                         
01976                                 //remake particle
01977                                 MakeParticle3D();
01978                                 p=foundparticle3d;      
01979                                 if(!foundparticle3d)return 0;
01980 
01981                                 foundparticle=1;
01982                                                                 
01983                                 foundparticle3d->particletype=Particle3D::electron;
01984                                 f.Fit3d(1);
01985                         
01986                                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"EM FIT! " << f.par_a << " " << f.par_b << " +- "<<f.par_b_err<<" "<<f.par_e0<<endl;
01987                 
01988                                 p->emfit_a=f.par_a;
01989                                 p->emfit_b=f.par_b;
01990                                 p->emfit_e0=f.par_e0;
01991                                 p->emfit_a_err=f.par_a_err;
01992                                 p->emfit_b_err=f.par_b_err;
01993                                 p->emfit_e0_err=f.par_e0_err;
01994                                 p->emfit_prob=f.prob;
01995                                 p->calibrated_energy=f.par_e0;
01996                                 p->emfit_chisq=f.chisq;
01997                                 p->emfit_ndf=f.ndf;
01998 
01999                                 p->pred_e_a=f.pred_e_a; 
02000                                 p->pred_g_a=f.pred_g_a;
02001                                 p->pred_b=f.pred_b;
02002                                 p->pred_e0=f.pred_e0;
02003                                 p->pred_e_chisq=f.pred_e_chisq;
02004                                 p->pred_e_ndf=f.pred_e_ndf;
02005                                 p->pred_g_chisq=f.pred_g_chisq;
02006                                 p->pred_g_ndf=f.pred_g_ndf;                             
02007                         
02008                                 p->pre_over=f.pre_over;         
02009                                 p->pre_under=f.pre_under;               
02010                                 p->post_over=f.post_over;               
02011                                 p->post_under=f.post_under;             
02012 
02013                                 p->pp_chisq=f.pp_chisq;
02014                                 p->pp_ndf=f.pp_ndf;
02015                                 p->pp_igood=f.pp_igood;
02016                                 p->pp_p=f.pp_p;
02017         
02018         
02019                                 p->cmp_chisq = f.cmp_chisq;
02020                                 p->cmp_ndf = f.cmp_ndf;
02021                                 p->peakdiff = f.peakdiff;       
02022                                 
02023                         
02024                                 if(!f.conv)p->particletype=Particle3D::other;
02025                 
02027                         
02028                         
02029                         }
02030                         
02031                 
02032 //              }
02033 //      }
02034 
02035 
02036 
02037         MSG("PrimaryShowerFinder",Msg::kDebug)<<"printing long Shower chains\nu\n";
02038         if(chain_u)chain_u->PrintChain();MSG("PrimaryShowerFinder",Msg::kDebug)<<"v\n";
02039         if(chain_v)chain_v->PrintChain();
02040         MSG("PrimaryShowerFinder",Msg::kDebug)<<"done\n";
02041         
02042 
02043         MSG("PrimaryShowerFinder",Msg::kDebug)<<"done looking for long Showers\n";
02044         return foundparticle;
02045 }
02046 
02047 
02048 
02049 //once we make chains in both view, we need to see if one chain is too short compared to the other
02050 //if so, try to extend it by finding near hits
02051 void PrimaryShowerFinder::ExpandShowerChain(Chain * chain,int view, double start_z, double end_z)
02052 {
02053 
02054         if(chain->start_z < start_z && chain->end_z > end_z)return;
02055         
02056 //      printf("expanding shower chain from %f to %f\n",chain->start_z,chain->end_z);
02057                                 
02058         if(chain->start_z > start_z)
02059         {
02060                 double z=chain->start_z-0.06;
02061                 while(z>start_z-0.04)
02062                 {
02063                         //find hits in this z range and add the best one to the chain
02064                         std::vector<int> cs = cluster_manager->FindClustersInZ(z,view);
02065         
02066         //              printf("found %d hits in z %f\n",cs.size(),z);
02067                         
02068                         int closest=-1;
02069                         double closest_d=1000;
02070                         //for now do a simple straight line from start t
02071                         for(unsigned int i=0;i<cs.size();i++)   
02072                         {
02073                                 Managed::ManagedCluster *c = cluster_manager->GetCluster(cs[i]);
02074         //                      printf("id %d z %f t %f e %f\n",cs[i],c->z,c->t,c->e);
02075                         
02076                                 if(fabs(c->t-chain->start_t)<closest_d)
02077                                 {
02078                                         closest_d=fabs(c->t-chain->start_t);
02079                                         closest=c->id;
02080                                 }
02081                         }
02082                         
02083                         if(closest>-1 && closest_d < 0.1)
02084                         {
02085                                 Managed::ManagedCluster *c = cluster_manager->GetCluster(closest);
02086                                 chain->insert(c->t,c->z,c->e,c->id);
02087                         }
02088         
02089                         z-= 0.06; //subtract off a plane
02090                 }
02091         }
02092 
02093 
02094                 
02095         if(chain->end_z < end_z)
02096         {
02097                 double z=chain->end_z+0.06;
02098                 while(z<end_z+0.04)
02099                 {
02100                         //find hits in this z range and add the best one to the chain
02101                         std::vector<int> cs = cluster_manager->FindClustersInZ(z,view);
02102         
02103         //              printf("found %d hits in z %f\n",cs.size(),z);
02104                         
02105                         int closest=-1;
02106                         double closest_d=1000;
02107                         //for now do a simple straight line from start t
02108                         for(unsigned int i=0;i<cs.size();i++)   
02109                         {
02110                                 Managed::ManagedCluster *c = cluster_manager->GetCluster(cs[i]);
02111         //                      printf("id %d z %f t %f e %f\n",cs[i],c->z,c->t,c->e);
02112                         
02113                                 if(fabs(c->t-chain->end_t)<closest_d)
02114                                 {
02115                                         closest_d=fabs(c->t-chain->end_t);
02116                                         closest=c->id;
02117                                 }
02118                         }
02119                         
02120                         if(closest>-1  && closest_d < 0.1)
02121                         {
02122                                 Managed::ManagedCluster *c = cluster_manager->GetCluster(closest);
02123                                 chain->insert(c->t,c->z,c->e,c->id);
02124                         }
02125         
02126                         z+= 0.06; //subtract off a plane
02127                 }
02128         }
02129 
02130 
02131 }
02132 
02133 
02134 
02135 
02136 //make sure that the chains both start within one plane of each other... 
02137 //for now, simply do this by finding the plane below the most upstream hit
02138 void PrimaryShowerFinder::ClearFrontVertex(Chain * chain_u, Chain * chain_v)
02139 {
02140         double start_z_u = chain_u->start_z;
02141         double start_z_v = chain_v->start_z;
02142 
02143         if(fabs(start_z_u-start_z_v)<0.04)return; //close enough
02144         
02145         double start_z = start_z_u < start_z_v ? start_z_v : start_z_u - 0.05;
02146         
02147         Chain * toclear =  start_z_u < start_z_v ? chain_u : chain_v;
02148         //now delete all hits up to start_z
02149         
02150         Chain tmp = *toclear;
02151         
02152         toclear->ClearHits();
02153         
02154         for(int i=0;i<tmp.entries;i++)
02155         {
02156                 if(tmp.z[i]>start_z)
02157                         toclear->insert(tmp.t[i], tmp.z[i], tmp.e[i], tmp.cluster_id[i]);
02158         }
02159         
02160 }
02161 
02162 //make sure that the chains actually overlap in z (so from the same particle)
02163 //and that the chain is of sufficient length
02164 int PrimaryShowerFinder::CheckChainOverlap(Chain * chain_u, Chain * chain_v)
02165 {
02166 
02167         
02168 
02169         
02170         //now make sure that there is sufficient overlap
02171         double require_overlap_distance=1.0;
02172         
02173         double start_z = chain_u->start_z > chain_v->start_z ? chain_u->start_z : chain_v->start_z;
02174         double end_z = chain_u->end_z < chain_v->end_z ? chain_u->end_z : chain_v->end_z;
02175 
02176         MSG("PrimaryShowerFinder",Msg::kDebug)<<"overlap "<<end_z-start_z<<"\n";
02177         
02178         if(end_z-start_z < require_overlap_distance)return 0;
02179         return 1;
02180         
02181         
02182 
02183 }
02184 
02185 
02186 // check to see if this is a Shower-like chain with a sufficient number of hits 
02187 int PrimaryShowerFinder::CheckChainQuality(Chain *c)
02188 {
02189         if(!c)return 0;
02190         if(c->entries<1)return 0;       
02191         
02192         int qual=1;
02193         
02194         //look at Shower-like strip fraction
02195         //and sparsity (# planes with hits/#planes from front to back of Shower)
02196         
02197         int Showerlike=0;
02198         int hitplane=0;
02199         
02200         double last_hitplane=-100;
02201         
02202         double min = 1000000;
02203         double max = 0;
02204         for(int i=0;i<c->entries;i++)
02205         {
02206                 if(c->e[i]>0.5 && c->e[i]<2.5)Showerlike++;
02207                 if(fabs(last_hitplane - c->z[i])>0.03)
02208                 {
02209                         hitplane++;
02210                         last_hitplane=c->z[i];
02211                 }
02212                 
02213                 min = c->z[i]< min? c->z[i]:min;
02214                 max = c->z[i]> max? c->z[i]:max;
02215         }
02216         
02217         double Showerfrac = (double)Showerlike / (double)c->entries;
02218         double sparsity = max-min?(double)hitplane / (max-min) / 7.08 :1; //length/size of u+v plane
02219         
02220 
02221 //definition of shwerfrac is wrong!     
02222 //      if (Showerfrac < 0.5) qual=0;
02223         if (sparsity < 0.8) qual=0;
02224         
02225         if(c->muonfrac>0.5)qual=0;
02226         
02227         
02228         MSG("PrimaryShowerFinder",Msg::kDebug)<<"MUON FRAC "<<c->muonfrac <<"Shower chain check  Showerfrac "<< Showerfrac<<" sparsity "<< sparsity<<"  qual "<< qual<<"\n";
02229 
02230         return qual;
02231 }
02232 
02233 
02234 Chain * PrimaryShowerFinder::FindShowerChain(Managed::ClusterManager *cl, int view)
02235 {
02236 
02237         MSG("PrimaryShowerFinder",Msg::kDebug)<<"\n\nLooking for long Shower Chain in view "<<view<<"\n";
02238 
02239 
02240         std::map<double, std::map<double, int> > * cluster_map = cl->GetClusterMap(view);
02241 
02242 
02243     std::map<double, std::map<double, int> >::reverse_iterator p_iterr;
02244     std::map<double, int >::iterator s_iterr;
02245 
02246 
02247         Chain * c = new Chain();
02248 
02249         Managed::ManagedCluster * largest=0;
02250         double last_plane=100000;
02251         
02252         std::vector<Managed::ManagedCluster *> maxplaneclusters;
02253         //first we want to find maximum hits in each plane
02254         //we will then check these hits to see if they can form a tightly aligned track
02255     for(p_iterr=cluster_map->rbegin();p_iterr!=cluster_map->rend(); p_iterr++)
02256     {   
02257         //new plane and a largest cluster in the previous plane?
02258                 if(fabs(last_plane-p_iterr->first)>0.03 && largest)
02259                 {
02260                         MSG("PrimaryShowerFinder",Msg::kDebug)<<"largest cluster found "<<largest->t<<" "<<largest->z<<" "<<largest->e<<"\n";
02261                         maxplaneclusters.push_back(largest);
02262                         largest=0;
02263                 }
02264         
02265         for(s_iterr=p_iterr->second.begin();s_iterr!=p_iterr->second.end(); s_iterr++)
02266         {
02267                         Managed::ManagedCluster * clus = cl->GetCluster(s_iterr->second);
02268                         if(!largest)
02269                         {
02270                                 largest=clus;
02271                                 continue;
02272                         }
02273                         if(largest->e<clus->e)largest=clus;
02274                 }
02275         }
02276         if(largest)
02277         {
02278                 MSG("PrimaryShowerFinder",Msg::kDebug)<<"largest cluster found "<<largest->t<<" "<<largest->z<<" "<<largest->e<<"\n";
02279                 maxplaneclusters.push_back(largest);
02280                 largest=0;
02281         }
02282 
02283         
02284         for(unsigned int i=0;i<maxplaneclusters.size();i++)
02285                 c->insert(maxplaneclusters[i]->t,maxplaneclusters[i]->z,maxplaneclusters[i]->e,maxplaneclusters[i]->id);                
02286                         
02287         return c;                       
02288 
02289 
02290 }
02291 
02292 
02293 void PrimaryShowerFinder::DumpParticle3D()
02294 {
02295         ostringstream s;
02296         
02297         s << "Long Shower Particle3D found "<<endl;
02298 
02299                 Particle3D * p = foundparticle3d;
02300                 if (p==0)return;
02301         
02302         
02303 
02304         
02305                 s << "\n---Particle3d "  << "---\nstart, stop (u,v,z) (" << p->start_u << ", "<< p->start_v << ", " << p->start_z <<") (" <<p->end_u<<", "<< p->end_v << ", " << p->end_z <<")"<<endl;
02306                 s << "entries "<< p->entries << "  muonlike " << p->muonfrac <<endl;
02307                 
02308 
02309                 
02310                 s<<"types: ";
02311                 for(unsigned int j=0;j<p->types.size();j++)
02312                 {
02313                         switch( p->types[j].type)
02314                         {
02315                                 case    ParticleType::em:
02316                                         s<<"em ";
02317                                         break;
02318                                 case ParticleType::emshort:
02319                                         s<<"emshort  ";
02320                                         break;                          
02321                                 case ParticleType::muon:
02322                                         s<<"muon ";
02323                                         break;
02324                                 case ParticleType::prot:
02325                                         s<<"prot ";
02326                                         break;          
02327                                 case ParticleType::pi0:
02328                                         s<<"pi0 ";
02329                                         break;
02330                                 case ParticleType::uniquemuon:
02331                                         s<<"uniqueShower ";
02332                                         break;  
02333                         
02334                         }
02335                                 
02336                 }
02337                 s<<endl;
02338                 
02339                 s<<"particletype : "<<p->particletype<<endl;
02340                 
02341                 
02342                 s<<"par a " << p->emfit_a << " par b "<< p->emfit_b << " par e0 "<< p->emfit_e0 << " cale "<<p->calibrated_energy<<" chisq " << p->emfit_chisq<<" ndf " << p->emfit_ndf<<endl;
02343                 s<<"emprob " << p->emfit_prob <<" avg rms_t "<<p->avg_rms_t<<endl;
02344                                 
02345                 s << "spoints (u,v,z,e - chain, chainhit, chainview - shared - rms_t - view) : ";
02346                 for(int j=0;j<p->entries;j++)
02347                         s << "(" << p->u[j]<<", "<<p->v[j]<<", "<<p->z[j]<<", "<<p->e[j]<<" - " << p->chain[j] <<", "<<p->chainhit[j]<<", "<<p->view[j]<<" - "<<p->shared[j]<<" - "<< p->rms_t[j]<< " - " << p->view[j]<<") ";
02348         
02349         /*
02350         s << "spoints (e - shared) : ";
02351                 for(int j=0;j<p->entries;j++)
02352                         s << "(" <<p->e[j]<<" - "<<p->shared[j]<<") ";
02353         
02354         */
02355         
02356         
02357 
02358                         
02359                         
02360 
02361                         
02362                         
02363                 s<<endl<<endl;
02364                 
02365                 
02366         
02367         
02368         
02369 
02370         MSG("PrimaryShowerFinder",Msg::kDebug) <<s.str();
02371 
02372         
02373 }
02374 
02375 
02376 
02377 void PrimaryShowerFinder::MakeParticle3D()
02378 {
02379 
02380         //printf("making psf particle cu %d cv %d\n",chain_u->myId,chain_v->myId);
02381 
02382         //verify that the chains in both views have sufficient overlap
02383         
02384         
02385         //make the 3d particle from these chains
02386         if(foundparticle3d)delete foundparticle3d;
02387         foundparticle3d= new Particle3D();
02388 
02389 
02390         
02391         std::vector<spoint> spoints;
02392         
02393         double start =0;
02394         double end  =1000;
02395         
02396         double endu=0;
02397         double endv=0;
02398                 
02399         
02400         
02401         
02402         int type=0;
02403 
02404 
02405                 Chain *c = chain_u;
02406                 
02407                 if(c->particletype)
02408                         type = c->particletype;
02409                 
02410                 for(int j=0;j<c->entries;j++)
02411                 {
02412                         if(c->parentChain==-1 && j==0)
02413                                 start = start < c->z[j]? c->z[j] : start;
02414                         if( j == c->entries-1)
02415                         {
02416                                 end = end < c->z[j] ? end : c->z[j];
02417                         }
02418                         endu=endu<c->z[j]?c->z[j]:endu;
02419 
02420                         
02421                         spoint a;
02422                         a.z=c->z[j];
02423                         a.t=c->t[j];
02424                         a.e=c->e[j];
02425                         a.chain=c->myId;
02426                         a.chainhit=j;
02427                         a.view=2;
02428                         
02429 
02430                         Managed::ManagedCluster * clu = cluster_manager->GetCluster(c->cluster_id[j]);
02431                 
02432                         if(clu)
02433                         {
02434                                 a.rms_t=clu->rms_t;
02435                                 spoints.push_back(a);
02436                         }
02437                 }
02438         
02439         
02440 
02441                 c = chain_v;
02442                 
02443                 if(c->particletype)
02444                         type = c->particletype;
02445                 for(int j=0;j<c->entries;j++)
02446                 {
02447                         if(c->parentChain==-1 && j==0)
02448                                 start = start < c->z[j]? c->z[j] : start;
02449                         if( j == c->entries-1)
02450                         {
02451                                 end = end < c->z[j] ? end : c->z[j];
02452                         }
02453                         endv=endv<c->z[j]?c->z[j]:endv;
02454                         
02455 
02456                         spoint a;
02457                         a.z=c->z[j];
02458                         a.t=c->t[j];
02459                         a.e=c->e[j];
02460                         a.chain=c->myId;
02461                         a.chainhit=j;
02462                         a.view=3;
02463                         
02464 
02465                         Managed::ManagedCluster * clu = cluster_manager->GetCluster(c->cluster_id[j]);
02466                         if(clu)
02467                         {
02468                                 a.rms_t=clu->rms_t;
02469                                 spoints.push_back(a);
02470                         }
02471                 }
02472         
02473         MSG("PrimaryShowerFinder",Msg::kDebug)<<"MAKING PARTICLE WITH "<<spoints.size()<<" POINTS\n";
02474         
02475         //we don't want the particle to extrapolate too far.... ...but now go all the way to the end
02476         double stopend = endu<endv?endv:endu;
02477 
02478                 
02479         sort(spoints.begin(), spoints.end(),spointgreaterlmf);
02480         
02481         for(unsigned int i=0;i<spoints.size();i++)
02482         {
02483 //              if( start - spoints[i].z > 0.045 )continue;
02484 //              if( spoints[i].z - end > 0.045)continue; //within 1 planes, we want the end of the chain
02485         
02486 //      printf("spoint %f %f %f %d\n",spoints[i].z,spoints[i].t,spoints[i].e,spoints[i].view);
02487         
02488                 if(spoints[i].z>stopend)continue;
02489         
02490                 int myview = spoints[i].view;
02491                 int lower=-1;
02492                 int upper=-1;
02493                 for(int j=i-1;j>-1;j--)
02494                 {
02495                         if(spoints[j].view!=myview)
02496                         {
02497                                 lower=j;
02498                                 break;
02499                         }
02500                 }
02501                 for(unsigned int j=i+1;j<spoints.size();j++)
02502                 {
02503                         if(spoints[j].view!=myview)
02504                         {
02505                                 upper=j;
02506                                 break;
02507                         }
02508                 }       
02509                 
02510                 
02511                 
02512                 double u = spoints[i].view==2 ? spoints[i].t : 0;
02513                 double v = spoints[i].view==3 ? spoints[i].t : 0;
02514                 double e = spoints[i].e;
02515                 double z = spoints[i].z;
02516                 int view = spoints[i].view;
02517                 
02518                 double rms_t = spoints[i].rms_t;
02519                 
02520 
02521                 
02522                 if(lower>-1 && upper > -1 )// good we can extrapolate!
02523                 {
02524                         double s = (spoints[upper].t - spoints[lower].t) /  ( spoints[upper].z - spoints[lower].z);
02525                         
02526                         double t = s * (spoints[i].z-spoints[lower].z) + spoints[lower].t;
02527                         
02528                         u = myview == 2 ? u : t;
02529                         v = myview == 3 ? v : t;
02530                         
02531 
02532                 }else if(lower>-1 && upper < 0)  //just user the closest other view t value
02533                 {
02534                 
02535 
02536                         //do we have another lower value?
02537                         int lower2=-1;
02538                         for(int j=lower-1;j>-1;j--)
02539                         {
02540                                 if(spoints[j].view!=myview)
02541                                 {
02542                                         lower2=j;
02543                                         break;
02544                                 }
02545                         }
02546                         
02547                         if(lower2>-1 && fabs( spoints[lower].z - spoints[lower2].z) >0)
02548                         {
02549                                 double s = (spoints[lower].t - spoints[lower2].t) /  ( spoints[lower].z - spoints[lower2].z);
02550                         
02551                                 double t = s * (spoints[i].z-spoints[lower2].z) + spoints[lower2].t;
02552                                 u = myview == 2 ? u : t;
02553                                 v = myview == 3 ? v : t;
02554                         
02555                         }else{
02556                                 u = myview == 2 ? u : spoints[lower].t;
02557                                 v = myview == 3 ? v : spoints[lower].t;
02558                         }
02559                 }
02560                 else if(upper>-1 && lower < 0)   //just user the closest other view t value
02561                 {
02562                 
02563 
02564                 
02565                         //do we have another upper value?
02566                         int upper2=-1;
02567                         for(unsigned int j=upper+1;j<spoints.size();j++)
02568                         {
02569                                 if(spoints[j].view!=myview)
02570                                 {
02571                                         upper2=j;
02572                                         break;
02573                                 }
02574                         }
02575                         
02576                         if(upper2>-1 && fabs( spoints[upper2].z - spoints[upper].z)>0)
02577                         {
02578                                 double s = (spoints[upper2].t - spoints[upper].t) /  ( spoints[upper2].z - spoints[upper].z);
02579                         
02580                                 double t = s * (spoints[i].z-spoints[upper].z) + spoints[upper].t;
02581                                 u = myview == 2 ? u : t;
02582                                 v = myview == 3 ? v : t;
02583                         
02584                         }else{
02585                                 u = myview == 2 ? u : spoints[upper].t;
02586                                 v = myview == 3 ? v : spoints[upper].t;
02587                         }
02588 
02589 
02590 
02591                         //lets use the vertex!
02592 
02593 /*
02594                         if(spoints[upper].view != spoints[i].view)
02595                                 upper=upper2;
02596                                 
02597                         if(spoints[upper].view == spoints[i].view)
02598                         {
02599 */
02600 //dont yet have a vertex!
02601 /*                              double vz =vtx_z;
02602                                 double vt = 0;
02603                                 vt = spoints[upper].view == 2 ? vtx_u : vt;
02604                                 vt = spoints[upper].view == 3 ? vtx_v : vt;
02605                         
02606                                 double t =vt;
02607                                 
02608                                 //if the spoint at z is not the same as z vtx, we need to extrapolate 
02609                                 if(fabs ( spoints[upper].z - vz)>0.001)
02610                                 {                       
02611                         
02612                                         double s = (spoints[upper].t - vt) /  ( spoints[upper].z - vz);
02613                                         t = s * (spoints[i].z-vz) + vt;                 
02614                         
02615                                 }
02616                                 
02617                 //              printf("---view %d u %f v %f t %f\n",myview,u,v,t);
02618                 //              printf("---vtx z %f t%f upper z %f t %f\n",vz,vt,spoints[upper].z,spoints[upper].t);
02619                         
02620 
02621                                 u = myview == 2 ? u : t;
02622                                 v = myview == 3 ? v : t;
02623 */
02624                                 u = myview == 2 ? u : spoints[upper].t;
02625                                 v = myview == 3 ? v : spoints[upper].t;
02626 
02627 //                      }       
02628                                 
02629 
02630                 }
02631                 else if(upper==-1 && lower==-1) //we have an empty view!!!
02632                 {
02633 
02634                         u = myview == 2 ? u : 0;
02635                         v = myview == 3 ? v : 0;
02636                 }
02637                 
02638                 foundparticle3d->add_to_back(u,v,z,e,spoints[i].chain,spoints[i].chainhit,view,rms_t);
02639                 
02640         //      printf("adding 3d spoint to particle %f %f %f --- %f\n",u,v,z,e);
02641         }
02642         
02643         
02644         if(type)
02645                 foundparticle3d->particletype=(Particle3D::EParticle3DType)type;
02646                 
02647         foundparticle3d->finalize();
02648         
02649         if(foundparticle3d->entries<1){delete foundparticle3d;foundparticle3d=0;}
02650 
02651 
02652 }
02653 
02654 
02655 
02656 void PrimaryShowerFinder::LoadCloseHits(HoughLine * hl, Managed::ClusterManager *cm, int view, double max_tdist,int limit_to_current_size)
02657 {
02658         if(!hl || !cm)return;
02659 
02660         std::map<double, std::map<double, int> > * cluster_map = cm->GetClusterMap(view);
02661     std::map<double, std::map<double, int> >::reverse_iterator p_iterr;
02662     std::map<double, int >::iterator s_iterr;
02663 
02664 
02665 //printf("\nloading hit\n");
02666         double last_z=0;
02667         Managed::ManagedCluster *closest_cluster=0;
02668         double dist=1000;
02669     for(p_iterr=cluster_map->rbegin();p_iterr!=cluster_map->rend(); p_iterr++)
02670     {   
02671         if(!last_z)last_z=p_iterr->first;
02672                 
02673   //    printf("plane %f \n",p_iterr->first);
02674         for(s_iterr=p_iterr->second.begin();s_iterr!=p_iterr->second.end(); s_iterr++)
02675                 {
02676                                 Managed::ManagedCluster *cl = cm->GetCluster(s_iterr->second);
02677                                 if(!cl)continue;
02678                                 if(cl->e<0.001)continue;
02679                                 
02680                                 //do we need to make sure that this cluster is contained withing the dimensions of the houghline?
02681                                 if(limit_to_current_size)
02682                                 {
02683                                         double min_z = hl->start_z < hl->end_z?hl->start_z:hl->end_z;
02684                                         double max_z = hl->start_z > hl->end_z?hl->start_z:hl->end_z;
02685                                         if(min_z-cl->z>0.07 || cl->z-max_z>0.07)
02686                                         {
02687                         //                      printf("not adding %f %f to hl with endpts %f %f  %f %f\n",cl->z,cl->t,hl->start_z,hl->start_t,hl->end_z,hl->end_t);
02688                                                 continue;
02689                                         }
02690                                         
02691                                         //if this hl is contrained to a single plane... we need to keep hits within t!
02692                                         double min_t =  hl->start_t < hl->end_t ? hl->start_t:hl->end_t;
02693                                         double max_t =  hl->start_t > hl->end_t ? hl->start_t:hl->end_t;
02694                                         
02695                                         if(max_z-min_z<0.001 && (cl->t-min_t < 0.001 || cl->t - max_t > 0.001))continue;
02696                                         
02697                                         if(fabs(hl->phi-3.141592/2 )<0.1)
02698                                         {
02699 
02700                                                 
02701                                                 if(cl->t<min_t || cl->t>max_t)continue;
02702                                         }
02703                                 
02704                                 }                               
02705                                 
02706                                 double exp_t = hl->GetExpectedT(cl->z);
02707 //                              printf(" cluster z %f\n",cl->z);
02708                                 //is it a vertical line?
02709                                 if(fabs(hl->phi-3.141592/2 )<0.1){
02710                                 
02711                                         //find where t=0
02712                                         double zpos = hl->offset_z + (hl->offset_t + hl->r/sin(hl->theta))*sin(hl->theta)/cos(hl->theta);
02713                         //      printf("verticle match at %f against %f\n",zpos,cl->z); 
02714                                         //give more for the verticle hits       
02715                                         double tmt = max_tdist > 0.1 ? max_tdist:0.1;
02716                                         if(fabs(zpos-cl->z)<tmt)
02717                                         {
02718                                                 hl->AddCluster(cl);
02719                                         }
02720                                 
02721                                 }else{
02722                                         //is it the closest?
02723                                         if(fabs(exp_t-cl->t)<max_tdist && fabs(exp_t-cl->t)<dist)
02724                                         {
02725                                                 closest_cluster=cl;
02726                                                 dist=fabs(exp_t-cl->t);
02727                                         }
02728                                 }       
02729                                         
02730                 }
02731 
02732                 std::map<double, std::map<double, int> >::reverse_iterator next_p_iterr=p_iterr;
02733                 next_p_iterr++;
02734         if(next_p_iterr!=cluster_map->rend())
02735         if(fabs(p_iterr->first-next_p_iterr->first)>0.04 && closest_cluster)
02736         {//new plane.. save the old plane
02737                         hl->AddCluster(closest_cluster);
02738 //                      printf("next plane %f from %f\n",p_iterr->first,last_z);
02739   //            printf("adding %f %f %f\n",closest_cluster->z,closest_cluster->t,closest_cluster->e);
02740                 closest_cluster=0;
02741                 dist=1000;
02742                 
02743                 last_z=p_iterr->first;
02744         }
02745 
02746 
02747         }
02748 
02749         if(closest_cluster)
02750         {//new plane.. save the old plane
02751  //     printf("last_plane %f\n",last_z);
02752  //                     printf("adding %f %f %f\n",closest_cluster->z,closest_cluster->t,closest_cluster->e);
02753                 hl->AddCluster(closest_cluster);
02754         }
02755 }
02756 
02757 
02758 void PrimaryShowerFinder::MakeHoughMap(int view)
02759 {
02760         if(view !=2 && view !=3)return;
02761         std::map<double, std::map<double, int>  >::iterator p_iterr;
02762     std::map<double, int>::iterator s_iterr;
02763         
02764 //      printf("----view %d----\n\n",view);
02765         
02766         if(houghmapU && view==2){houghmapU->Reset();}//delete houghmapU;houghmapU=0;}
02767         if(houghmapV && view==3){houghmapV->Reset();}//delete houghmapV;houghmapV=0;}
02768         
02769 //      if(houghmapU==0)houghmapU= new TH2D("hhu","Hough U", 300, 0,3.141592,150,-1,1);
02770 //      if(houghmapV==0)houghmapV= new TH2D("hhv","Hough V", 300, 0,3.141592,150,-1,1);
02771 
02772         if(houghmapU==0)houghmapU= new TH2D("hhu","Hough U", 100, 0,3.141592,50,-1,1);
02773         if(houghmapV==0)houghmapV= new TH2D("hhv","Hough V", 100, 0,3.141592,50,-1,1);
02774 
02775 
02776                 
02777         TH2D * h = view==2?houghmapU:view==3?houghmapV:0;
02778         if(!h)return;
02779 
02780         //is the map empty?
02781         if(cluster_manager->GetClusterMap(view)->begin()==cluster_manager->GetClusterMap(view)->end())return; 
02782 
02783  //being on the first hit can caus some trouble... so move a bit away from it
02784         double off_z = cluster_manager->GetClusterMap(view)->begin()->first-0.1;
02785         double off_t = cluster_manager->GetClusterMap(view)->begin()->second.begin()->first-0.1;
02786         //Managed::ManagedCluster *largest=0;
02787         //double last_plane=0;
02788   
02789   
02790         TH2D copy;
02791         h->Copy(copy);    
02792 
02793         //find max value!
02794         double maxvale=0;
02795     for(p_iterr=cluster_manager->GetClusterMap(view)->begin();p_iterr!=cluster_manager->GetClusterMap(view)->end(); p_iterr++)
02796     {   
02797         for(s_iterr=p_iterr->second.begin();s_iterr!=p_iterr->second.end(); s_iterr++)
02798         {
02799                         Managed::ManagedCluster * clus = cluster_manager->GetCluster(s_iterr->second);
02800                         if(clus->e>maxvale)
02801                         {
02802                                 maxvale=clus->e;
02803                         }
02804                 }
02805         }
02806 
02807 
02808     for(p_iterr=cluster_manager->GetClusterMap(view)->begin();p_iterr!=cluster_manager->GetClusterMap(view)->end(); p_iterr++)
02809     {   
02810         for(s_iterr=p_iterr->second.begin();s_iterr!=p_iterr->second.end(); s_iterr++)
02811         {
02812                         Managed::ManagedCluster * clus = cluster_manager->GetCluster(s_iterr->second);
02813 //printf("adding with e %f\n",clus->e);
02814                         //don't care about low e hits...
02815                         if(0.2 > clus->e)continue;
02816                         for(double i=0;i<3.141592;i+=3.141592/(h->GetNbinsX()-1))
02817                         {
02818                                 //double val = clus->e;
02819                                 double val=1;
02820                                 h->Fill(i, (clus->z-off_z) * cos(i) + (clus->t-off_t) * sin(i),  val);  
02821                                 copy.Fill(i, (clus->z-off_z) * cos(i) + (clus->t-off_t) * sin(i), clus->e);
02822                         }                       
02823                 }
02824         }
02825         
02826 
02827    TH2D ch;
02828    h->Copy(ch);
02829         
02830 
02831         h->Reset();
02832         
02833         //int nx=ch.GetNbinsX();
02834         //int ny=ch.GetNbinsY();   
02835 
02836         while(ch.GetMaximum()>1)
02837         {
02838                 int x;
02839                 int y;
02840                 int z;
02841                 int m = ch.GetMaximumBin();
02842                 ch.GetBinXYZ(m,x,y,z);
02843 
02844                 SaveHitGroup(&ch,(TH2D*)h, (double)ch.GetBinContent(x,y),(double)ch.GetBinContent(x,y), x,y);
02845     }    
02846 
02847 
02848         multimap<double, std::pair<double,double> > top_map;
02849         int tops=0;     
02850 
02851 
02853 //make a list of all peaks....to speed up the TH1::GetMaximum(Bin) which is painfully slow....
02854 
02855         std::vector<int> peakbinx;
02856         std::vector<int> peakbiny;
02857         std::vector<double> peakx;
02858         std::vector<double> peaky;
02859         std::vector<double> peakvalue;
02860         std::vector<int> peakstillgood;
02861         
02862         for(int i=1;i<h->GetNbinsX()+1;i++)
02863         for(int j=1;j<h->GetNbinsY()+1;j++)
02864         {
02865                 double thisc = h->GetBinContent(i,j);
02866                 if(thisc<0.0001)continue;
02867                 
02868                 if(
02869                         thisc >= h->GetBinContent(i+1,j)
02870                         &&
02871                         thisc >= h->GetBinContent(i,j+1)
02872                         &&
02873                         thisc >= h->GetBinContent(i-1,j)
02874                         &&
02875                         thisc >= h->GetBinContent(i,j-1)
02876                         &&
02877                         thisc >= h->GetBinContent(i+1,j+1)
02878                         &&
02879                         thisc >= h->GetBinContent(i+1,j-1)
02880                         &&
02881                         thisc >= h->GetBinContent(i-1,j+1)
02882                         &&
02883                         thisc >= h->GetBinContent(i-1,j-1)
02884                 )
02885                 {
02886                         peakx.push_back(h->GetXaxis()->GetBinCenter(i));
02887                         peaky.push_back(h->GetYaxis()->GetBinCenter(j));
02888                         peakbinx.push_back(i);
02889                         peakbiny.push_back(j);
02890                         peakvalue.push_back(thisc);     
02891                         peakstillgood.push_back(1);
02892                 }
02893         
02894         }
02895 
02896 /*
02897         int npeaks = peakstillgood.size();
02898         
02899         printf("found %d peaks",npeaks);
02900         
02901         while(npeaks>0)
02902         {
02903         
02904                 double xv=0;
02905                 double yv=0;
02906                 double val=0;
02907                 int cnt =0;
02908                         
02909                 GetPeakAreaAverage(xv, yv,val, peakbinx,peakbiny,peakx,peaky,peakvalue,peakstillgood);
02910 
02911 
02912                 top_map.insert(make_pair((double)val, std::pair<double,double>(xv,yv)) );
02913                 tops++;
02914 
02915                 npeaks=0;
02916                 for(unsigned int i=0;i<peakstillgood.size();i++)        
02917                         npeaks+=peakstillgood[i]==1?1:0;
02918         }
02919         
02920         printf("%d peaks made\n",tops);
02921 */
02923 
02924 //there must be a better way....
02925 
02926         int npeaks=0;
02927         for(unsigned int i=0;i<peakstillgood.size();i++)        
02928                 npeaks+=peakstillgood[i]==1?1:0;
02929         
02930         
02931         
02932         while(npeaks>0)//h->GetMaximum()>0)
02933         {
02934                 double xv=0;
02935                 double yv=0;
02936                 double val=0;
02937                 int cnt =0;
02938                 
02939                 int x;
02940                 int y;
02941                 //int z;
02942 //              int m = h->GetMaximumBin();
02943 //              h->GetBinXYZ(m,x,y,z);
02944 
02945                 double max=0;
02946                 int maxi=-1;
02947                 for(int i=0;i<(int)peakvalue.size();i++)
02948                 {
02949                         if(peakvalue[i]>max && peakstillgood[i])
02950                         {
02951                                 max = peakvalue[i];
02952                                 maxi=i;
02953                         }
02954                 }
02955 
02956                 if(maxi<0)break;
02957                 
02958                 x=peakbinx[maxi];
02959                 y=peakbiny[maxi];
02960                 
02961                 GetPeakAreaAverage(xv, yv,val, cnt, x, y, (TH2D*)h, peakvalue, peakstillgood, peakbinx, peakbiny);
02962                 xv/=(double)cnt;
02963                 yv/=(double)cnt;
02964                 
02965                 val=copy.GetBinContent(copy.FindBin(xv,yv));
02966                 
02967                 top_map.insert(make_pair((double)val, std::pair<double,double>(xv,yv)) );
02968                 tops++;
02969                 
02970                 npeaks=0;
02971                 for(unsigned int i=0;i<peakstillgood.size();i++)        
02972                         npeaks+=peakstillgood[i]==1?1:0;
02973         
02974         }
02975 
02976 
02977 
02978 
02979         
02980         
02981 //      printf("printing top %d\n",tops);
02982         
02983         multimap<double, std::pair<double,double> >::reverse_iterator it = top_map.rbegin();
02984 
02985         
02986         std::vector<HoughLine> *houghlines = view==2?&houghlinesU:&houghlinesV;
02987         
02988         for(int k=0;k<tops;k++)
02989         {
02990                 double theta=it->second.first;
02991                 double r = it->second.second;   
02992                 HoughLine hl(theta, r, off_t, off_z);
02993                 
02994         //      printf("peak at theta %f r %f off_t %f off_z %f\n",theta,r,off_t,off_z);
02995         
02996                 LoadCloseHits(&hl,cluster_manager,view,0.0412);
02997                         
02998                 if(hl.ncluster>1)       
02999                         houghlines->push_back(hl);
03000 
03001                 it++;
03002 
03003                 
03004         }
03005         
03006         
03007         sort(&(*houghlines)[0],&(*houghlines)[houghlines->size()],CompareLength);
03008         
03009         //should remove duplicates --- (caused by the same clusters being matched to different hough map peaks)
03010         std::vector<HoughLine> houghlinestemp;
03011 
03012         houghlinestemp = *houghlines;
03013         houghlines->clear();
03014 
03015         for(int i=0;i<(int)houghlinestemp.size();i++)
03016         {
03017                 int dup=0;
03018                 HoughLine * l = &houghlinestemp[i];
03019                 for(int j=i+1;j<(int)houghlinestemp.size();j++)
03020                 {
03021                         HoughLine *r = &houghlinestemp[j];
03022                         if(l && r && l->start_z == r->start_z && l->start_t == r->start_t && l ->end_z == r->end_z && l->end_t == r->end_t)
03023                         {
03024                                 dup=1;
03025                                 break;
03026                         }
03027 
03028                 }
03029 
03030                 if(!dup && houghlinestemp[i].ncluster>1)houghlines->push_back(houghlinestemp[i]);
03031 
03032         }
03033 
03034 
03035 /* //this code is highly flawed...
03036         for(int i=0;i<houghlines->size();i++)
03037         {
03038                 HoughLine * l = &(*houghlines)[i];
03039                 HoughLine * r=0;
03040                 if(i<houghlines->size()-1)r = &(*houghlines)[i+1];
03041                 if((*houghlines)[i].ncluster>1)houghlinestemp.push_back((*houghlines)[i]);              
03042                 if(l && r && l->start_z == r->start_z && l->start_t == r->start_t && l ->end_z == r->end_z && l->end_t == r->end_t)
03043                 {
03044 
03045                         while (i<houghlines->size())
03046                         {
03047                                 i++;
03048                                 HoughLine * l = &(*houghlines)[i];
03049                                 HoughLine * r = &(*houghlines)[i+1];
03050                                 if(! (l->start_z == r->start_z && l->start_t == r->start_t && l ->end_z == r->end_z && l->end_t == r->end_t)) break;
03051                         }
03052                 }
03053         }
03054 
03055 */
03056 
03057         
03058 
03059         
03060         sort(&(*houghlines)[0],&(*houghlines)[houghlines->size()],CompareForwardAndClusters);
03061         
03062         //int max =-1;
03063         //if(houghlines->size()>10)max = houghlines->size()-10;
03064         for(int i=0;i<(int)houghlines->size();i++)
03065         {
03066                 if((*houghlines)[i].ncluster<1)break;
03067 //              printf("%d houghline with e %f clusters %d (t,z) (%f,%f) (%f,%f) chi2/ndf %f theta %f r %f phiz %f \n", i,(*houghlines)[i].sum_e, (*houghlines)[i].ncluster, (*houghlines)[i].start_t, (*houghlines)[i].start_z, (*houghlines)[i].end_t, (*houghlines)[i].end_z, (*houghlines)[i].chi2/(*houghlines)[i].ncluster,(*houghlines)[i].theta,(*houghlines)[i].r,(*houghlines)[i].phi);
03068                 
03069                                                                         
03070         }
03071         
03072         
03073 
03074 }
03075 
03076 
03077  
03078 void PrimaryShowerFinder::SaveHitGroup(TH2D * his, TH2D * saver, double save_val, double with_val, int curx, int cury)
03079 {
03080         if(curx<1 || cury < 1 || curx> his->GetNbinsX() || cury >his->GetNbinsY())return;
03081         
03082 
03083         
03084         double thisval = his->GetBinContent(curx,cury);
03085 
03086 //      printf("saving %d %d thisval %f saveval %f withval %f\n",curx,cury,thisval,save_val, with_val);
03087         if(thisval==0)return;
03088         if(fabs(thisval-save_val)<0.001)
03089         {
03090                 saver->SetBinContent(curx, cury,thisval);
03091                 //thisval=0;
03092         }
03093         else if(thisval>=with_val)return;
03094 
03095 
03096         his->SetBinContent(curx,cury,0);
03097                         
03098         SaveHitGroup(his,saver,save_val,thisval,curx+1,cury);
03099         SaveHitGroup(his,saver,save_val,thisval,curx,cury+1);
03100         SaveHitGroup(his,saver,save_val,thisval,curx-1,cury);
03101         SaveHitGroup(his,saver,save_val,thisval,curx,cury-1);
03102         SaveHitGroup(his,saver,save_val,thisval,curx+1,cury+1);
03103         SaveHitGroup(his,saver,save_val,thisval,curx-1,cury-1);
03104         SaveHitGroup(his,saver,save_val,thisval,curx+1,cury-1);
03105         SaveHitGroup(his,saver,save_val,thisval,curx-1,cury+1);
03106 //      printf("save recursion done\n");
03107         
03108 
03109         //printf("saving points %d %d %f\n",curx,cury,thisval);
03110 
03111 }
03112 
03113 /*
03114 void PrimaryShowerFinder::GetPeakAreaAverage(double &xv, double &yv, double &val, std::vector<int> &peakbinx, std::vector<int> &peakbiny, std::vector<double> &peakx, std::vector<double> &peaky, std::vector<double> &peakvalue, std::vector<int> &peakstillgood)
03115 {
03116         xv=0;
03117         yv=0;
03118         val=0;
03119         
03120         //goal is to find all ajacent bins and group them into one bin....
03121         
03122         int peak=-1;
03123         
03124         //find the first available peak
03125         for(unsigned int i=0;i<peakstillgood.size();i++)
03126         {
03127                 if(peakstillgood[i]==1)
03128                 {
03129                         peak=i;
03130                         break;
03131                 }
03132         }
03133         
03134         if(peak<0)return;
03135         
03136         val=peakvalue[peak];
03137         
03138         int binx=peakbinx[peak];
03139         int biny=peakbiny[peak];
03140         
03141         int lastaddcount=1;
03142         
03143         std::vector<int> addedx;
03144         std::vector<int> addedy;
03145         
03146         addedx.push_back(binx);
03147         addedy.push_back(biny);
03148         peakstillgood[peak]=0;
03149         
03150         while(lastaddcount>0)
03151         {
03152                 lastaddcount=0;
03153                 for(unsigned int i=0;i<peakstillgood.size();i++)
03154                 {
03155                         if(peakstillgood[i])
03156                         {
03157                                 for(unsigned int j=0;j<addedx.size();j++)
03158                                 {
03159                                 
03160                                         if(peakbinx[i] == peakbinx[addedx[j]]+1  && peakbiny[i] == peakbiny[addedy[j]]
03161                                                 ||
03162                                                 peakbinx[i] == peakbinx[addedx[j]]  && peakbiny[i] == peakbiny[addedy[j]]+1
03163                                                 ||
03164                                                 peakbinx[i] == peakbinx[addedx[j]]-1  && peakbiny[i] == peakbiny[addedy[j]]
03165                                                 ||
03166                                                 peakbinx[i] == peakbinx[addedx[j]]  && peakbiny[i] == peakbiny[addedy[j]]-1
03167                                                 ||
03168                                                 peakbinx[i] == peakbinx[addedx[j]]+1  && peakbiny[i] == peakbiny[addedy[j]]+1
03169                                                 ||
03170                                                 peakbinx[i] == peakbinx[addedx[j]]+1  && peakbiny[i] == peakbiny[addedy[j]]-1
03171                                                 ||
03172                                                 peakbinx[i] == peakbinx[addedx[j]]-1  && peakbiny[i] == peakbiny[addedy[j]]+1
03173                                                 ||
03174                                                 peakbinx[i] == peakbinx[addedx[j]]-1  && peakbiny[i] == peakbiny[addedy[j]]-1
03175                                         ){
03176                                                 peakstillgood[i]=0;
03177                                                 addedx.push_back(peakbinx[i]);
03178                                                 addedy.push_back(peakbiny[i]);
03179                                                 lastaddcount++;
03180                                                 break;
03181                                         }
03182                                         
03183                                         
03184                                         
03185                         
03186                                 }
03187                         }
03188                 
03189                 }       
03190         }
03191         
03192         
03193         //calculate the value
03194         
03195         double cnt=0;
03196         for(unsigned int i=0;i<addedx.size();i++)
03197         {
03198                 xv+=peakx[addedx[i]];
03199                 yv+=peaky[addedy[i]];
03200                 cnt+=1;
03201         }
03202         
03203         xv/=cnt;
03204         yv/=cnt;
03205         
03206 
03207 }*/
03208 
03209 
03210 //assuming that the region is bounded by bins with 0 as their content
03211 void PrimaryShowerFinder::GetPeakAreaAverage(double &x, double &y,double &val, int & cnt, int curx, int cury, TH2D * hist, std::vector<double> & peakvalue, std::vector<int> &peakstillgood, std::vector<int> &peakbinx,std::vector<int> & peakbiny)
03212 {
03213 
03214         double thisval = hist->GetBinContent(curx,cury);
03215 
03216         if(thisval==0)return;
03217         
03218         val+=thisval;
03219         x+=hist->GetXaxis()->GetBinCenter(curx);
03220         y+=hist->GetYaxis()->GetBinCenter(cury);
03221         cnt++;
03222         hist->SetBinContent(curx,cury,0);
03223 
03225         
03226         for(int i=0;i<(int)peakbinx.size();i++)
03227         {
03228                 if(peakbinx[i]==curx && peakbiny[i]==cury)
03229                 {
03230                         peakstillgood[i]=0;
03231                         peakvalue[i]=0;
03232                         break;
03233                 }
03234         }
03235         
03236         
03238 
03239 
03240         
03241         GetPeakAreaAverage(x, y, val,cnt, curx+1, cury, hist, peakvalue, peakstillgood, peakbinx, peakbiny);
03242         GetPeakAreaAverage(x, y, val,cnt, curx, cury+1, hist, peakvalue, peakstillgood, peakbinx, peakbiny);
03243         GetPeakAreaAverage(x, y, val,cnt, curx, cury-1, hist, peakvalue, peakstillgood, peakbinx, peakbiny);
03244         GetPeakAreaAverage(x, y, val,cnt, curx-1, cury, hist, peakvalue, peakstillgood, peakbinx, peakbiny);
03245         GetPeakAreaAverage(x, y, val,cnt, curx+1, cury+1, hist, peakvalue, peakstillgood, peakbinx, peakbiny);
03246         GetPeakAreaAverage(x, y, val,cnt, curx-1, cury-1, hist, peakvalue, peakstillgood, peakbinx, peakbiny);
03247         GetPeakAreaAverage(x, y, val,cnt, curx+1, cury-1, hist, peakvalue, peakstillgood, peakbinx, peakbiny);
03248         GetPeakAreaAverage(x, y, val,cnt, curx-1, cury+1, hist, peakvalue, peakstillgood, peakbinx, peakbiny);
03249 
03250                 
03251 }
03252 
03253 void PrimaryShowerFinder::FindHoughLineMatches()
03254 {
03255         houghlineMatch.clear();
03256 
03257         //need to iterate over the hough lines in both view and look for matches
03258 
03259         for(unsigned int i=0;i<houghlinesU.size();i++)
03260         {
03261                 //don't want to match to verticle lines
03262                 if(fabs(houghlinesU[i].phi -3.141592/2)<0.1)continue;
03263                 if(houghlinesU[i].ncluster<2)continue;
03264                         HoughLine * l = &houghlinesU[i];
03265 
03266                 for(unsigned int j=0;j<houghlinesV.size();j++)
03267                 {
03268                         //don't want to match to verticle lines
03269                         if(fabs(houghlinesV[j].phi - 3.141592/2)<0.1)continue;
03270                         if(houghlinesV[j].ncluster<2)continue;
03271 
03272                         HoughLine * r = &houghlinesV[j];
03273                         
03274                         double overlap = l->end_z<r->end_z?l->end_z:r->end_z;
03275                         overlap -= l->start_z > r->start_z ? l->start_z:r->start_z;
03276                         
03277                         //require overlap to exceed 90% of the length of the smaller line
03278                         double l_z = l->end_z-l->start_z;
03279                         double r_z = r->end_z-r->start_z;
03280                         double smaller_z = l_z < r_z?l_z:r_z;
03281                         if(smaller_z * 0.5 > overlap)
03282                         {
03283                         //      printf("failing %d %d with overlap %f > %f\n",i,j,smaller_z * 0.5,overlap);
03284                                 continue;
03285                         }
03286                         double smaller_e = l->sum_e < r->sum_e ? l->sum_e:r->sum_e;
03287                         double larger_e = l->sum_e > r->sum_e ? l->sum_e:r->sum_e;
03288 
03289                         //require smaller e to be more than 30% of larger e
03290                         if(larger_e*0.3 > smaller_e)
03291                         {
03292                 //              printf("failing %d %d with energy %f > %f\n",i,j,larger_e*0.5 , smaller_e);
03293                                 continue;
03294                         }
03295 
03296 
03297                         //require a peak hit that is at least 20% of total e
03298                 
03299 
03300                         
03301                         double max_hit_e=0;
03302                         double sum_hit_e=0;
03303                         for(unsigned int ka=0;ka<l->cluster_id.size();ka++)
03304                         {
03305                                 Managed::ManagedCluster * clus = cluster_manager->GetCluster(l->cluster_id[ka]);
03306                                 if(!clus){
03307                                         //printf("missing cluster %d\n",l->cluster_id[ka]);
03308                                         continue;
03309                                 }
03310                                 sum_hit_e+=clus->e;
03311                                 
03312                                 if(clus->e>max_hit_e)max_hit_e=clus->e;
03313                         
03314                         }
03315                         for(unsigned int ka=0;ka<r->cluster_id.size();ka++)
03316                         {
03317                                 Managed::ManagedCluster * clus = cluster_manager->GetCluster(r->cluster_id[ka]);
03318                                 if(!clus){
03319                                         //printf("missing cluster %d\n",r->cluster_id[ka]);
03320                                         continue;
03321                                 }                                       
03322                                 sum_hit_e+=clus->e;
03323                                 
03324                                 if(clus->e>max_hit_e)max_hit_e=clus->e;
03325                         
03326                         }       
03327                         if(sum_hit_e*0.2 > max_hit_e)
03328                         {
03329                         //      printf("continuing sum hit*.3 %f > maxhit %f\n",sum_hit_e*0.3, max_hit_e);
03330                                 continue;
03331                         }
03332                         if( max_hit_e<5)
03333                         {
03334                         //      printf("maxhit too small at %f\n", max_hit_e);
03335                                 continue;
03336                         }
03337                         
03338                 //      printf("saving %d %d\n",i,j);
03339                         houghlineMatch.push_back(make_pair(i,j));                                       
03340                 }
03341         }
03342 
03343         std::map<double,int> orderer;
03344         for(unsigned int i=0;i<houghlineMatch.size();i++)
03345         {
03346                 HoughLine * l = &houghlinesU[houghlineMatch[i].first];
03347                 HoughLine * r = &houghlinesV[houghlineMatch[i].second];
03348                 if(l->ncluster<2 || r->ncluster<2)continue;
03349                 double weight = (l->sum_e+r->sum_e)*(cos(l->phi))*(cos(r->phi));
03350                 weight *= l->ncluster / sqrt( (l->end_z-l->start_z)*(l->end_z-l->start_z)+(l->end_t-l->start_t)*(l->end_t-l->start_t));
03351                 weight *= r->ncluster / sqrt( (r->end_z-r->start_z)*(r->end_z-r->start_z)+(r->end_t-r->start_t)*(r->end_t-r->start_t));
03352 
03353                 orderer.insert(make_pair(weight,i));
03354         }               
03355 
03356         std::map<double,int>::iterator it_orderer;
03357 
03358         
03359         //printf("found %d preliminary matches\n",houghlineMatch.size());
03360 
03361         /*for(it_orderer = orderer.begin();it_orderer!=orderer.end();it_orderer++)
03362         //for(unsigned int i=0;i<houghlineMatch.size();i++)
03363         {
03364                 int i= it_orderer->second;
03365                 
03366                 //printf("match %d %d\n",houghlineMatch[i].first,houghlineMatch[i].second);
03367                 //HoughLine * l = &houghlinesU[houghlineMatch[i].first];
03368                 //HoughLine * r = &houghlinesV[houghlineMatch[i].second];
03369         
03370                 //printf("\t U s %f %f e %f %f sume %f clust %d phi %f\n",l->start_z,l->start_t, l->end_z,l->end_t,l->sum_e,l->ncluster,l->phi);
03371                 //printf("\t V s %f %f e %f %f sume %f clust %d phi %f\n",r->start_z,r->start_t, r->end_z,r->end_t,r->sum_e,r->ncluster,r->phi);
03372         }
03373 */
03374 }
03375 
03376 //make a chain using the primary shower
03377 Chain *  PrimaryShowerFinder::MakeShowerChain(int view)
03378 {
03379         if(houghlineMatch.size()<1)return 0;
03380         
03381         Chain * c = new Chain();
03382 
03383         //take the first match as the best for now...
03384         HoughLine * l = 0;
03385         /*if(view ==2)l=&houghlinesU[houghlineMatch[houghlineMatch.size()-1].first];
03386         else if(view ==3)l=&houghlinesV[houghlineMatch[houghlineMatch.size()-1].second];
03387         else return 0;
03388         */
03389 
03390 /*      
03391         if(view ==2)l=&houghlinesU[houghlinesU.size()-1];
03392         else if(view ==3)l=&houghlinesV[houghlinesV.size()-1];
03393         else return 0;
03394 */      
03395 
03396         if(view ==2)l=&houghlinesU[houghlineMatch[houghlineMatch.size()-1].first];
03397         else if(view ==3)l=&houghlinesV[houghlineMatch[houghlineMatch.size()-1].second];
03398         else return 0;
03399         
03400         
03401         l->primary=1;
03402         
03403         //printf("making chain for view %d\n",view);
03404         //printf("look at houghline %d\n",view==2?houghlineMatch[houghlineMatch.size()-1].first:houghlineMatch[houghlineMatch.size()-1].second);
03405         for(int i=0;i<l->ncluster;i++)
03406         {
03407         //      printf("getting cluster %d\n",l->cluster_id[i]);
03408                 Managed::ManagedCluster *clus = cluster_manager->GetCluster(l->cluster_id[i]);
03409                 if(!clus)continue;
03410                 c->insert(clus->t,clus->z,clus->e,clus->id);
03411         //      printf("inserting %f %f %f %d\n",clus->t,clus->z,clus->e,clus->id);     
03412         }
03413 
03414 
03415         return c;
03416 }
03417 
03418 
03419 void PrimaryShowerFinder::DumpHoughLines(int view)
03420 {
03421         
03422         if(!MsgService::Instance()->IsActive("PrimaryShowerFinder",Msg::kDebug))return;
03423         
03424         printf("\nView %d\n",view);
03425         std::vector<HoughLine>* houghlines=0;
03426         if(view==2)houghlines=&houghlinesU;
03427         if(view==3)houghlines=&houghlinesV;
03428         if(!houghlines)return;
03429 
03430         for(int i=(*houghlines).size()-1;i>-1;i--)
03431         {
03432                 if((*houghlines)[i].ncluster<1)continue;
03433                 printf("houghline with e %f clusters %d (t,z) (%f,%f) (%f,%f) chi2/ndf %f  offz %f expt %f phi %f\n", (*houghlines)[i].sum_e, (*houghlines)[i].ncluster, (*houghlines)[i].start_t, (*houghlines)[i].start_z, (*houghlines)[i].end_t, (*houghlines)[i].end_z, (*houghlines)[i].chi2/(*houghlines)[i].ncluster,(*houghlines)[i].offset_z,(*houghlines)[i].GetExpectedT((*houghlines)[i].offset_z),(*houghlines)[i].phi);
03434         }
03435 
03436         printf("\n");
03437 }
Generated on Mon Feb 15 11:07:26 2010 for loon by  doxygen 1.3.9.1