00001 #include "NueAna/ParticlePID/ParticleFinder/Managed/HitManager.h"
00002
00003 #include <math.h>
00004 #include <map>
00005
00006 #include "MessageService/MsgService.h"
00007
00008 #include <sstream>
00009
00010 #include <cstdlib>
00011
00012 CVSID("$Id: HitManager.cxx,v 1.3 2009/09/11 05:00:40 gmieg Exp $");
00013
00014
00015
00016 using namespace Managed;
00017
00018 ClassImp(HitManager)
00019
00020 HitManager::HitManager()
00021 {
00022 Reset();
00023 }
00024
00025 HitManager::~HitManager(){}
00026
00027
00028 Managed::ManagedHit * HitManager::FindHit(int id)
00029 {
00030 for(unsigned int i=0;i<hits.size();i++)
00031 if(hits[i].id==id)return &hits[i];
00032
00033 return 0;
00034 }
00035
00036 void HitManager::Reset()
00037 {
00038 hits.clear();
00039 ManagedHit::ResetIDCounter();
00040 }
00041
00042 int HitManager::InsertHit(int view, int plane, int strip, double z, double t, double e)
00043 {
00044 Managed::ManagedHit h(view,plane,strip,z,t,e);
00045 hits.push_back(h);
00046 return h.id;
00047 }
00048
00049 Managed::ManagedHit *HitManager::FindHit(int , int plane, int strip)
00050 {
00051 for(unsigned int i=0;i<hits.size();i++)
00052 if(hits[i].GetPlane()==plane && hits[i].GetStrip()==strip)return &hits[i];
00053 return 0;
00054 }
00055
00056
00057 Managed::ManagedHit *HitManager::FindHit(int , double z, double t)
00058 {
00059 for(unsigned int i=0;i<hits.size();i++)
00060 if(fabs(hits[i].GetZ()-z)<0.01 && fabs(hits[i].GetT()-t)<0.01)return &hits[i];
00061 return 0;
00062 }
00063
00064
00065 std::vector<ManagedHit> HitManager::GetAvailableHits()
00066 {
00067 std::vector<ManagedHit> ret;
00068 for(unsigned int i=0;i<hits.size();i++)
00069 if(hits[i].GetERemaining()>0)ret.push_back(hits[i]);
00070
00071 return ret;
00072
00073 }
00074
00075
00076 void HitManager::ClearXTalk()
00077 {
00078
00079
00080 std::vector<int> view;
00081 std::vector<int> plane;
00082 std::vector<int> strip;
00083 std::vector<double> t;
00084 std::vector<double> z;
00085 std::vector<double> energy;
00086
00087
00088 std::multimap<double,int> sorter_map;
00089
00090 std::map<int, std::map<int, int> > loc_map;
00091
00092
00093 for(unsigned int i=0;i<hits.size();i++)
00094 {
00095 view.push_back(hits[i].GetView());
00096 plane.push_back(hits[i].GetPlane());
00097 strip.push_back(hits[i].GetStrip());
00098 t.push_back(hits[i].GetT());
00099 z.push_back(hits[i].GetZ());
00100 energy.push_back(hits[i].GetERemaining());
00101
00102
00103
00104
00105 }
00106
00107
00108 double thresh = 3;
00109
00110 for(unsigned int i=0;i<plane.size();i++)
00111 {
00112 sorter_map.insert(std::pair <double,int>(energy[i],i));
00113 loc_map[plane[i]][strip[i]]=i;
00114 }
00115
00116 ostringstream os;
00117
00118 double reassignedxtalke=0.0;
00119
00120 int foundxtalk=0;
00121 std::map<int, std::map<int, int> >::iterator p_iter;
00122 for(p_iter=loc_map.begin();p_iter!=loc_map.end(); p_iter++)
00123 {
00124 std::map<int, int>::iterator s_iter;
00125 for(s_iter=p_iter->second.begin();s_iter!=p_iter->second.end(); s_iter++)
00126 {
00127 if (energy[s_iter->second]<thresh)continue;
00128
00129 std::map<int, int>::iterator s_iter2;
00130 for(s_iter2=p_iter->second.begin();s_iter2!=p_iter->second.end(); s_iter2++)
00131 {
00132 if(s_iter==s_iter2)continue;
00133 int sp=abs(s_iter->first-s_iter2->first);
00134 os<<"sep "<<sp<<" view "<<view[s_iter->second]<<" plane "<<plane[s_iter->second]<<" high "<<energy[s_iter->second]<<" low "<<energy[s_iter2->second]<<"\n";
00135
00136 if( sp>9 && sp<14)
00137 {
00138 if(energy[s_iter2->second] < 0.3 * energy[s_iter->second])
00139 {
00140 reassignedxtalke+=energy[s_iter2->second];
00141
00142 energy[s_iter->second]+=energy[s_iter2->second];
00143 energy[s_iter2->second]=0.0;
00144 foundxtalk=1;
00145 }
00146 }
00147
00148 }
00149 }
00150 }
00151
00152
00153
00154
00155
00156 if(foundxtalk)
00157 {
00158 std::vector<int>tplane;
00159 std::vector<int>tstrip;
00160 std::vector<int>tview;
00161 std::vector<double>tt;
00162 std::vector<double>tz;
00163 std::vector<double>tenergy;
00164
00165
00166 for(unsigned int i=0;i<energy.size();i++)
00167 {
00168 if(energy[i]>0.001)
00169 {
00170 tplane.push_back(plane[i]);
00171 tstrip.push_back(strip[i]);
00172 tenergy.push_back(energy[i]);
00173 tt.push_back(t[i]);
00174 tz.push_back(z[i]);
00175 tview.push_back(view[i]);
00176 }
00177 }
00178
00179 plane=tplane;
00180 strip=tstrip;
00181 view=tview;
00182 t=tt;
00183 z=tz;
00184 energy=tenergy;
00185
00186 os<<"Removing XTalk ---- "<<hits.size() <<" hits merged to make ";
00187 Reset();
00188 for(unsigned int i=0;i<plane.size();i++)
00189 {
00190 Managed::ManagedHit h(view[i],plane[i],strip[i],z[i],t[i],energy[i]);
00191 hits.push_back(h);
00192 }
00193 os << (int)hits.size() <<" hits\n";
00194
00195 }
00196
00197
00198 sorter_map.clear();
00199 loc_map.clear();
00200
00201 MSG("HitManager",Msg::kDebug)<<os;
00202
00203 }
00204
00205
00206