#include <AlgShowerSSList.h>
Inheritance diagram for AlgShowerSSList:
Public Member Functions | |
| AlgShowerSSList () | |
| virtual | ~AlgShowerSSList () |
| virtual void | RunAlg (AlgConfig &ac, CandHandle &ch, CandContext &cx) |
| virtual void | Trace (const char *c) const |
Private Attributes | |
| TH1F * | vtxUHistAll |
| TH1F * | vtxUHist |
| TH1F * | vtxVHistAll |
| TH1F * | vtxVHist |
Definition at line 15 of file AlgShowerSSList.h.
| AlgShowerSSList::AlgShowerSSList | ( | ) |
Definition at line 57 of file AlgShowerSSList.cxx.
00058 { 00059 00060 vtxUHistAll = new TH1F("vtxUHistAll","vtxUHistAll",1,-4.,+4); 00061 vtxUHist = new TH1F("vtxUHist","vtxUHist",1,-4.,+4); 00062 vtxVHistAll = new TH1F("vtxVHistAll","vtxVHistAll",1,-4.,+4); 00063 vtxVHist = new TH1F("vtxVHist","vtxVHist",1,-4.,+4); 00064 00065 }
| AlgShowerSSList::~AlgShowerSSList | ( | ) | [virtual] |
Definition at line 68 of file AlgShowerSSList.cxx.
References vtxUHist, vtxUHistAll, vtxVHist, and vtxVHistAll.
00069 { 00070 00071 delete vtxUHistAll; 00072 delete vtxUHist; 00073 delete vtxVHistAll; 00074 delete vtxVHist; 00075 00076 }
| void AlgShowerSSList::RunAlg | ( | AlgConfig & | ac, | |
| CandHandle & | ch, | |||
| CandContext & | cx | |||
| ) | [virtual] |
We begin with a nested loop over all CandSlices, and over all CandClusters within each CandSlice, generating list of all clusters in the U and V planes. Next, we execute a nested loop over U-view and V-view CandClusters, checking for matching 2D clusters. The following criteria are used: Beginning planes for the two 2D 'long' CandClusters (defined by PlaneScale) must be within DiffViewPlaneMatch., or DiffViewPlaneMatchShort if one or more clusters is short. Beginning times for the two 2D CandClusters must be within DiffViewTimeMatch. The fractional overlap of U and V clusters must exceed DiffViewPlaneCoverage, if both clusters are long. The total pulse heights of the U and V clusters must agree to within DiffViewPulseHeightCut if at least one of the clusters is long. (Note: This should probably be changed to a charge ratio) If a U/V cluster set is found which satisfies these requirements, a second loop over U clusters is performed, searching for alternative U/V pairs involving the original V cluster satisfying the requirements above, but higher in total energy. If no better alternative is found a CandShower is constructed. In this way, a given CandCluster will only be used once, in the shower with the highest total energy.
Reimplemented from AlgShowerSRList.
Definition at line 80 of file AlgShowerSSList.cxx.
References CandHandle::AddDaughterLink(), Mphysical::c_light, MuELoss::e, Registry::Get(), AlgFactory::GetAlgHandle(), CandContext::GetCandRecord(), CandContext::GetDataIn(), CandHandle::GetDaughterIterator(), CandShowerHandle::GetEnergy(), CandSubShowerSRHandle::GetEnergy(), AlgFactory::GetInstance(), CandSubShowerSRHandle::GetMaxStpTime(), CandSubShowerSRHandle::GetMinStpTime(), CandContext::GetMom(), CandHandle::GetUidInt(), RecMinos::GetVldContext(), Msg::kDebug, ClusterType::kEMLike, ClusterType::kHadLike, CalDigitType::kSigCorr, ClusterType::kTrkLike, PlaneView::kU, PlaneView::kV, Msg::kWarning, PlaneView::kX, PlaneView::kY, CandShowerSR::MakeCandidate(), MSG, PITCHINMETRES, CandRecoHandle::SetCandSlice(), vtxUHist, vtxUHistAll, vtxVHist, and vtxVHistAll.
00081 { 00082 00083 assert(cx.GetDataIn()); 00084 if (!(cx.GetDataIn()->InheritsFrom("TObjArray"))) { 00085 return; 00086 } 00087 00088 const CandSliceListHandle *slicelist = 0; 00089 const CandClusterListHandle *clusterlist = 0; 00090 const CandSubShowerSRListHandle *subshowerlist = 0; 00091 const TObjArray *cxin = dynamic_cast<const TObjArray *>(cx.GetDataIn()); 00092 for (Int_t i=0; i<=cxin->GetLast(); i++) { 00093 TObject *tobj = cxin->At(i); 00094 if (tobj->InheritsFrom("CandSliceListHandle")) { 00095 slicelist = dynamic_cast<CandSliceListHandle*>(tobj); 00096 } 00097 if (tobj->InheritsFrom("CandClusterListHandle")) { 00098 clusterlist = dynamic_cast<CandClusterListHandle*>(tobj); 00099 } 00100 if (tobj->InheritsFrom("CandSubShowerSRListHandle")) { 00101 subshowerlist = dynamic_cast<CandSubShowerSRListHandle*>(tobj); 00102 } 00103 } 00104 if (!slicelist || !clusterlist || !subshowerlist) { 00105 MSG("AlgShowerSS",Msg::kWarning) << 00106 "CandSliceListHandle, CandClusterListHandle or CandSubShowerSRListHandle missing\n"; 00107 return; 00108 } 00109 00110 00111 CandContext cxx(this,cx.GetMom()); 00112 00113 Double_t MinNonPhysicsEnergy = 0.2; //in GeV 00114 Int_t MaxPlaneGap = 2; 00115 Int_t NPeakFindingBins = 32; 00116 Double_t SubShowerTimingCut = 50.0; //in ns 00117 Double_t EnergyAsymmetryCut = 0.3; 00118 Double_t PeakMergeTimeCut = 10.0; //in ns 00119 //get config for AlgSubShowerSS 00120 const char *charShowerSSAlgConfig = 0; 00121 ac.Get("ShowerSSAlgConfig",charShowerSSAlgConfig); 00122 ac.Get("MaxPlaneGap",MaxPlaneGap); 00123 ac.Get("NPeakFindingBins",NPeakFindingBins); 00124 ac.Get("MinNonPhysicsEnergy",MinNonPhysicsEnergy); 00125 ac.Get("SubShowerTimingCut",SubShowerTimingCut); 00126 ac.Get("EnergyAsymmetryCut",EnergyAsymmetryCut); 00127 ac.Get("PeakMergeTimeCut",PeakMergeTimeCut); 00128 00129 //reset all hists and set bins 00130 // (add an extra bin at each end outside of expected range, 00131 // in order to be able to detect peaks at (far) detector edges) 00132 00133 vtxUHistAll->Reset(); 00134 vtxUHistAll->SetBins(NPeakFindingBins+2, 00135 -4. - 8./Double_t(NPeakFindingBins), 00136 +4. + 8./Double_t(NPeakFindingBins)); 00137 vtxVHistAll->Reset(); 00138 vtxVHistAll->SetBins(NPeakFindingBins+2, 00139 -4. - 8./Double_t(NPeakFindingBins), 00140 +4. + 8./Double_t(NPeakFindingBins)); 00141 vtxUHist->Reset(); 00142 vtxUHist->SetBins(NPeakFindingBins+2, 00143 -4. - 8./Double_t(NPeakFindingBins), 00144 +4. + 8./Double_t(NPeakFindingBins)); 00145 vtxVHist->Reset(); 00146 vtxVHist->SetBins(NPeakFindingBins+2, 00147 -4. - 8./Double_t(NPeakFindingBins), 00148 +4. + 8./Double_t(NPeakFindingBins)); 00149 00150 //Get singleton instance of AlgFactory 00151 AlgFactory &af = AlgFactory::GetInstance(); 00152 AlgHandle ah = af.GetAlgHandle("AlgShowerSS",charShowerSSAlgConfig); 00153 00154 const CandRecord *candrec = cx.GetCandRecord(); 00155 assert(candrec); 00156 const VldContext *vldcptr = candrec->GetVldContext(); 00157 assert(vldcptr); 00158 VldContext vldc = *vldcptr; 00159 00160 UgliGeomHandle ugh(vldc); 00161 00162 CandSliceHandleItr sliceItr(slicelist->GetDaughterIterator()); 00163 Int_t nslice = 0; 00164 Int_t nsubshower = 0; 00165 00166 while (CandSliceHandle *slice = sliceItr()) { //slice loop 00167 00168 MSG("AlgShowerSS",Msg::kDebug) << "Slice: " << nslice << endl; 00169 00170 // Select SubShowers within this slice 00171 CandSubShowerSRHandleItr subshowerItr(subshowerlist->GetDaughterIterator()); 00172 00173 // iterate over subshowers, and make a shower out of 00174 // all subshowers in each slice, in contiguous planes 00175 int plane[500] = {0}; 00176 int firstPlane = 500; 00177 int lastPlane = 0; 00178 nsubshower = 0; 00179 while (CandSubShowerSRHandle *subshower = subshowerItr()) { 00180 if (*subshower->GetCandSlice()==*slice) { 00181 MSG("AlgShowerSS",Msg::kDebug) << "SubShower: " << nsubshower 00182 << " in slice " << nslice << endl; 00183 for(int i=subshower->GetBegPlane();i<=subshower->GetEndPlane();i++){ 00184 MSG("AlgShowerSS",Msg::kDebug) << "Plane: " << i << endl; 00185 plane[i] = 1; 00186 if(i<firstPlane) firstPlane = i; 00187 if(i>lastPlane) lastPlane = i; 00188 } 00189 } 00190 nsubshower++; 00191 } 00192 00193 MSG("AlgShowerSS",Msg::kDebug) << "First Plane = " << firstPlane << endl; 00194 MSG("AlgShowerSS",Msg::kDebug) << "Last Plane = " << lastPlane << endl; 00195 00196 //separate up into plane chunks according to MaxPlaneGap 00197 Int_t nshower = 0; 00198 while(firstPlane<=lastPlane){ //while loop over planes 00199 00200 Int_t begPlane = 500; 00201 Int_t endPlane = 0; 00202 Int_t nGaps = 0; 00203 00204 for(int i=firstPlane;i<=lastPlane;i++){ 00205 if(plane[i]==0) { 00206 if(begPlane!=500) nGaps+=1; 00207 } 00208 else { 00209 if(begPlane==500) begPlane = i; 00210 endPlane = i; 00211 } 00212 plane[i] = 0; 00213 if(nGaps>MaxPlaneGap) { 00214 MSG("AlgShowerSS",Msg::kDebug) << "begPlane = " << begPlane << endl; 00215 MSG("AlgShowerSS",Msg::kDebug) << "endPlane = " << endPlane << endl; 00216 break; 00217 } 00218 } 00219 firstPlane = endPlane+1; 00220 00221 //could have multiple showers per plane chunk 00222 //pull out all subshowers in plane chunk and look for TPos peaks: 00223 00224 //loop through subshower list again and histogram tposvtx 00225 //in order to look for peaks 00226 00227 vtxUHistAll->Reset(); 00228 vtxVHistAll->Reset(); 00229 vtxUHist->Reset(); 00230 vtxVHist->Reset(); 00231 00232 subshowerItr.ResetFirst(); 00233 nsubshower = 0; 00234 while (CandSubShowerSRHandle *subshower = subshowerItr()) { 00235 if (*subshower->GetCandSlice()==*slice) { 00236 if(subshower->GetBegPlane()>=begPlane && 00237 subshower->GetBegPlane()<=endPlane){ 00238 if(subshower->GetPlaneView()==PlaneView::kU || 00239 subshower->GetPlaneView()==PlaneView::kX) { 00240 MSG("AlgShowerSS",Msg::kDebug) 00241 << "Filling vtxUHist with Subshower " << nsubshower 00242 << " with vtvU = " << subshower->GetVtxU() 00243 << " and energy = " << subshower->GetEnergy() 00244 << " and slope = " << subshower->GetSlope() 00245 << " and AvgDev = " << subshower->GetAvgDev() 00246 << endl; 00247 CandStripHandleItr stripssItr(subshower->GetDaughterIterator()); 00248 while (CandStripHandle *stp = stripssItr()) { 00249 //fill with max charge for this tpos 00250 Double_t val = stp->GetCharge(CalDigitType::kSigCorr) - 00251 vtxUHistAll->GetBinContent(vtxUHistAll->FindBin(stp->GetTPos())); 00252 if(val>0) vtxUHistAll->Fill(stp->GetTPos(),val); 00253 if( subshower->GetClusterID()==ClusterType::kEMLike || 00254 subshower->GetClusterID()==ClusterType::kHadLike || 00255 subshower->GetClusterID()==ClusterType::kTrkLike){ 00256 val = stp->GetCharge(CalDigitType::kSigCorr) - 00257 vtxUHist->GetBinContent(vtxUHist->FindBin(stp->GetTPos())); 00258 if(val>0) vtxUHist->Fill(stp->GetTPos(),val); 00259 } 00260 } 00261 } 00262 else if(subshower->GetPlaneView()==PlaneView::kV || 00263 subshower->GetPlaneView()==PlaneView::kY){ 00264 MSG("AlgShowerSS",Msg::kDebug) 00265 << "Filling vtxVHist with Subshower " << nsubshower 00266 << " with vtvV = " << subshower->GetVtxV() 00267 << " and energy = " << subshower->GetEnergy() 00268 << " and slope = " << subshower->GetSlope() 00269 << " and AvgDev = " << subshower->GetAvgDev() 00270 << endl; 00271 CandStripHandleItr stripssItr(subshower->GetDaughterIterator()); 00272 while (CandStripHandle *stp = stripssItr()) { 00273 //fill with max charge for this tpos 00274 Double_t val = stp->GetCharge(CalDigitType::kSigCorr) - 00275 vtxVHistAll->GetBinContent(vtxVHistAll->FindBin(stp->GetTPos())); 00276 if(val>0) vtxVHistAll->Fill(stp->GetTPos(),val); 00277 if( subshower->GetClusterID()==ClusterType::kEMLike || 00278 subshower->GetClusterID()==ClusterType::kHadLike || 00279 subshower->GetClusterID()==ClusterType::kTrkLike){ 00280 val = stp->GetCharge(CalDigitType::kSigCorr) - 00281 vtxVHist->GetBinContent(vtxVHist->FindBin(stp->GetTPos())); 00282 if(val>0) vtxVHist->Fill(stp->GetTPos(),val); 00283 } 00284 } 00285 } 00286 } 00287 } 00288 nsubshower++; 00289 } 00290 00291 if(false){ 00292 TCanvas *can = new TCanvas(); 00293 can->Divide(2,2); 00294 can->cd(1); 00295 vtxUHistAll->Draw(); 00296 can->cd(2); 00297 vtxVHistAll->Draw(); 00298 can->cd(3); 00299 vtxUHist->Draw(); 00300 can->cd(4); 00301 vtxVHist->Draw(); 00302 char nomus[256]; 00303 sprintf(nomus,"peakPlot_slice%i.eps",nslice); 00304 can->Print(nomus); 00305 delete can; 00306 } 00307 00308 //width of shower in units of number of bins 00309 //assume shower typically has a width of ~4 strips 00310 Double_t sigma = (4.*0.0417/8.)*(NPeakFindingBins); 00311 //don't let width be <=1 bin to prevent 00312 //against clipping error messages 00313 if(sigma<=1.0) sigma = 1.01; 00314 00315 Int_t nUPeaks = 0; 00316 Int_t nUPeaks_save = 0; 00317 Int_t nUPeaksAll = 0; 00318 TSpecReal_t peakUPos[10] = {0}; 00319 TSpecReal_t peakUVal[10] = {0}; 00320 TSpecReal_t tmp_peakUPos = vtxUHist->GetBinCenter(vtxUHist->GetMaximumBin()); 00321 TSpecReal_t tmp_peakUVal = vtxUHist->GetMaximum(); 00322 MSG("AlgShowerSS",Msg::kDebug) << "U Hist Integrals: " 00323 << vtxUHist->Integral() << " " 00324 << vtxUHistAll->Integral() << endl; 00325 for(int bins = 1;bins<=vtxUHist->GetNbinsX();bins++){ 00326 if(vtxUHist->GetBinContent(bins)>0) nUPeaks++; 00327 if(vtxUHistAll->GetBinContent(bins)>0) nUPeaksAll++; 00328 } 00329 if(nUPeaks==1) { 00330 peakUPos[0] = tmp_peakUPos; 00331 peakUVal[0] = tmp_peakUVal; 00332 } 00333 else if(nUPeaks>0) { 00334 TSpectrum spec(10); 00335 TSpecReal_t *source = new TSpecReal_t[vtxUHist->GetNbinsX()]; 00336 TSpecReal_t *dest = new TSpecReal_t[vtxUHist->GetNbinsX()]; 00337 for(int h_loop=0;h_loop<vtxUHist->GetNbinsX();h_loop++){ 00338 source[h_loop] = vtxUHist->GetBinContent(h_loop+1); 00339 dest[h_loop] = 0; 00340 } 00341 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,9,1) 00342 nUPeaks = spec.SearchHighRes(source,dest, 00343 vtxUHist->GetNbinsX(), 00344 sigma,5,0,1,0,1); 00345 #else 00346 nUPeaks = spec.Search1HighRes(source,dest, 00347 vtxUHist->GetNbinsX(), 00348 sigma,5,0,1,0,1); 00349 #endif 00350 00351 TSpecReal_t *posX = spec.GetPositionX(); 00352 TSpecReal_t *posX_copy = new TSpecReal_t[nUPeaks]; 00353 for(int p_loop=0;p_loop<nUPeaks;p_loop++) posX_copy[p_loop] = posX[p_loop]; 00354 //make sure peaks are in ascending order 00355 for(int p_loop=0;p_loop<nUPeaks-1;p_loop++) { 00356 if(posX_copy[p_loop+1]<posX_copy[p_loop]) { 00357 TSpecReal_t temp = posX_copy[p_loop]; 00358 posX_copy[p_loop] = posX_copy[p_loop+1]; 00359 posX_copy[p_loop+1] = temp; 00360 p_loop = -1; 00361 } 00362 } 00363 for(int p_loop=0;p_loop<nUPeaks;p_loop++){ 00364 Int_t bin = 1+Int_t(posX_copy[p_loop]+0.5); 00365 peakUPos[p_loop] = vtxUHist->GetBinCenter(bin); 00366 peakUVal[p_loop] = vtxUHist->GetBinContent(bin); 00367 } 00368 delete [] posX_copy; 00369 delete [] source; 00370 delete [] dest; 00371 } 00372 if(nUPeaksAll>1) { 00373 TSpectrum spec(10); 00374 TSpecReal_t *source = new TSpecReal_t[vtxUHistAll->GetNbinsX()]; 00375 TSpecReal_t *dest = new TSpecReal_t[vtxUHistAll->GetNbinsX()]; 00376 for(int h_loop=0;h_loop<vtxUHistAll->GetNbinsX();h_loop++){ 00377 source[h_loop] = vtxUHistAll->GetBinContent(h_loop+1); 00378 dest[h_loop] = 0; 00379 } 00380 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,9,1) 00381 nUPeaksAll = spec.SearchHighRes(source,dest, 00382 vtxUHistAll->GetNbinsX(), 00383 sigma,5,0,1,0,1); 00384 #else 00385 nUPeaksAll = spec.Search1HighRes(source,dest, 00386 vtxUHistAll->GetNbinsX(), 00387 sigma,5,0,1,0,1); 00388 #endif 00389 delete [] source; 00390 delete [] dest; 00391 } 00392 nUPeaks_save = nUPeaks; 00393 if(nUPeaksAll==0 && nUPeaks>0) nUPeaksAll = nUPeaks; 00394 00395 Int_t nVPeaks = 0; 00396 Int_t nVPeaks_save = 0; 00397 Int_t nVPeaksAll = 0; 00398 TSpecReal_t peakVPos[10] = {0}; 00399 TSpecReal_t peakVVal[10] = {0}; 00400 TSpecReal_t tmp_peakVPos = vtxVHist->GetBinCenter(vtxVHist->GetMaximumBin()); 00401 TSpecReal_t tmp_peakVVal = vtxVHist->GetMaximum(); 00402 MSG("AlgShowerSS",Msg::kDebug) << "V Hist Integrals: " 00403 << vtxVHist->Integral() << " " 00404 << vtxVHistAll->Integral() << endl; 00405 for(int bins = 1;bins<=vtxVHist->GetNbinsX();bins++){ 00406 if(vtxVHist->GetBinContent(bins)>0) nVPeaks++; 00407 if(vtxVHistAll->GetBinContent(bins)>0) nVPeaksAll++; 00408 } 00409 if(nVPeaks==1) { 00410 peakVPos[0] = tmp_peakVPos; 00411 peakVVal[0] = tmp_peakVVal; 00412 } 00413 else if(nVPeaks>0) { 00414 TSpectrum spec(10); 00415 TSpecReal_t *source = new TSpecReal_t[vtxVHist->GetNbinsX()]; 00416 TSpecReal_t *dest = new TSpecReal_t[vtxVHist->GetNbinsX()]; 00417 for(int h_loop=0;h_loop<vtxVHist->GetNbinsX();h_loop++){ 00418 source[h_loop] = vtxVHist->GetBinContent(h_loop+1); 00419 dest[h_loop] = 0; 00420 } 00421 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,9,1) 00422 nVPeaks = spec.SearchHighRes(source,dest, 00423 vtxVHist->GetNbinsX(), 00424 sigma,5,0,1,0,1); 00425 #else 00426 nVPeaks = spec.Search1HighRes(source,dest, 00427 vtxVHist->GetNbinsX(), 00428 sigma,5,0,1,0,1); 00429 #endif 00430 00431 TSpecReal_t *posX = spec.GetPositionX(); 00432 TSpecReal_t *posX_copy = new TSpecReal_t[nVPeaks]; 00433 for(int p_loop=0;p_loop<nVPeaks;p_loop++) posX_copy[p_loop] = posX[p_loop]; 00434 //make sure peaks are in ascending order 00435 for(int p_loop=0;p_loop<nVPeaks-1;p_loop++) { 00436 if(posX_copy[p_loop+1]<posX_copy[p_loop]) { 00437 TSpecReal_t temp = posX_copy[p_loop]; 00438 posX_copy[p_loop] = posX_copy[p_loop+1]; 00439 posX_copy[p_loop+1] = temp; 00440 p_loop = -1; 00441 } 00442 } 00443 for(int p_loop=0;p_loop<nVPeaks;p_loop++){ 00444 Int_t bin = 1+Int_t(posX_copy[p_loop]+0.5); 00445 peakVPos[p_loop] = vtxVHist->GetBinCenter(bin); 00446 peakVVal[p_loop] = vtxVHist->GetBinContent(bin); 00447 } 00448 delete [] posX_copy; 00449 delete [] source; 00450 delete [] dest; 00451 } 00452 if(nVPeaksAll>1) { 00453 TSpectrum spec(10); 00454 TSpecReal_t *source = new TSpecReal_t[vtxVHistAll->GetNbinsX()]; 00455 TSpecReal_t *dest = new TSpecReal_t[vtxVHistAll->GetNbinsX()]; 00456 for(int h_loop=0;h_loop<vtxVHistAll->GetNbinsX();h_loop++){ 00457 source[h_loop] = vtxVHistAll->GetBinContent(h_loop+1); 00458 dest[h_loop] = 0; 00459 } 00460 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,9,1) 00461 nVPeaksAll = spec.SearchHighRes(source,dest, 00462 vtxVHistAll->GetNbinsX(), 00463 sigma,5,0,1,0,1); 00464 #else 00465 nVPeaksAll = spec.Search1HighRes(source,dest, 00466 vtxVHistAll->GetNbinsX(), 00467 sigma,5,0,1,0,1); 00468 #endif 00469 delete [] source; 00470 delete [] dest; 00471 } 00472 nVPeaks_save = nVPeaks; 00473 if(nVPeaksAll==0 && nVPeaks>0) nVPeaksAll = nVPeaks; 00474 00475 //if no peaks in either view 00476 if(nUPeaksAll==0 || nVPeaksAll==0) { 00477 //can't form a 3D shower in these conditions! 00478 MSG("AlgShowerSS",Msg::kDebug) 00479 << "nUPeaks = " << nUPeaks 00480 << " nVPeaks = " << nVPeaks 00481 << " Can't form a 3D shower" << endl; 00482 continue; 00483 } 00484 00485 //one or both views have no good "physics" subshower: 00486 if(nUPeaks==0 || nVPeaks==0){ 00487 MSG("AlgShowerSS",Msg::kDebug) << "Have something in both views but " 00488 << "nUPeaks = " << nUPeaks 00489 << "and nVPeaks = " << nVPeaks 00490 << endl; 00491 //try to form anyway with whatever's there 00492 if(nUPeaks==0) nUPeaks = 1; 00493 if(nVPeaks==0) nVPeaks = 1; 00494 } 00495 00496 //Try to form some 3D showers, yeah! 00497 00498 //have one or more peak => maybe one or 00499 //more shower in this plane range 00500 00501 //get boundaries that define each shower in U view: 00502 std::vector<TSpecReal_t> peakUPosLow(nUPeaks,0); 00503 std::vector<TSpecReal_t> peakUPosUpp(nUPeaks,0); 00504 00505 MSG("AlgShowerSS",Msg::kDebug) << "Number of U Peaks found = " 00506 << nUPeaks << " with positions, heights " 00507 << "and bounds: " << endl; 00508 for (int i=0;i<nUPeaks;i++){ 00509 if(i==0) peakUPosLow[i] = -5; 00510 else peakUPosLow[i] = peakUPosUpp[i-1]; 00511 if(i==nUPeaks-1) peakUPosUpp[i] = 5; 00512 else { 00513 Int_t bin = vtxUHist->GetXaxis()->FindBin(peakUPos[i]); 00514 Int_t bestBin = vtxUHist->GetNbinsX(); 00515 TSpecReal_t threeBinAve = 1000000; //sigcor! 00516 while(bin<vtxUHist->GetNbinsX() && 00517 vtxUHist->GetBinCenter(bin)<peakUPos[i+1]) { 00518 TSpecReal_t ave = (vtxUHist->GetBinContent(bin-1) + 00519 vtxUHist->GetBinContent(bin) + 00520 vtxUHist->GetBinContent(bin+1)); 00521 if(ave<threeBinAve) { 00522 threeBinAve = ave; 00523 bestBin = bin; 00524 } 00525 bin++; 00526 } 00527 peakUPosUpp[i] = vtxUHist->GetBinCenter(bestBin); 00528 } 00529 MSG("AlgShowerSS",Msg::kDebug) << peakUPos[i] << " " 00530 << peakUVal[i] << " [" 00531 << peakUPosLow[i] << "," 00532 << peakUPosUpp[i] << "]" 00533 << endl; 00534 } 00535 00536 //get boundaries that define each shower in V view: 00537 std::vector<TSpecReal_t> peakVPosLow(nVPeaks,0); 00538 std::vector<TSpecReal_t> peakVPosUpp(nVPeaks,0); 00539 MSG("AlgShowerSS",Msg::kDebug) << "Number of V Peaks found = " 00540 << nVPeaks 00541 << " with positions, heights " 00542 << "and bounds: " << endl; 00543 for (int i=0;i<nVPeaks;i++){ 00544 if(i==0) peakVPosLow[i] = -5; 00545 else peakVPosLow[i] = peakVPosUpp[i-1]; 00546 if(i==nVPeaks-1) peakVPosUpp[i] = 5; 00547 else { 00548 Int_t bin = vtxVHist->GetXaxis()->FindBin(peakVPos[i]); 00549 Int_t bestBin = vtxVHist->GetNbinsX(); 00550 TSpecReal_t threeBinAve = 1000000; //gev! 00551 while(bin<vtxVHist->GetNbinsX() && 00552 vtxVHist->GetBinCenter(bin)<peakVPos[i+1]) { 00553 TSpecReal_t ave = (vtxVHist->GetBinContent(bin-1) + 00554 vtxVHist->GetBinContent(bin) + 00555 vtxVHist->GetBinContent(bin+1)); 00556 if(ave<threeBinAve) { 00557 threeBinAve = ave; 00558 bestBin = bin; 00559 } 00560 bin++; 00561 } 00562 peakVPosUpp[i] = vtxVHist->GetBinCenter(bestBin); 00563 } 00564 MSG("AlgShowerSS",Msg::kDebug) << peakVPos[i] << " " 00565 << peakVVal[i] << " [" 00566 << peakVPosLow[i] << "," 00567 << peakVPosUpp[i] << "]" 00568 << endl; 00569 } 00570 00571 //for matching up U/V views of showers: 00572 std::vector<TSpecReal_t> totUEnergy(nUPeaks,0); 00573 std::vector<Double_t> aveUTime(nUPeaks,0); 00574 std::vector<TSpecReal_t> approxVfromTime(nUPeaks,0); 00575 std::vector<TSpecReal_t> maxPhysU(nUPeaks,0); 00576 std::vector<TSpecReal_t> minPhysU(nUPeaks,0); 00577 std::vector<Int_t> maxPlaneU(nUPeaks,0); 00578 std::vector<Int_t> minPlaneU(nUPeaks,0); 00579 std::vector<Int_t> nSubShowersU(nUPeaks,0); 00580 for(int i=0;i<nUPeaks;i++) { 00581 totUEnergy[i] = 0; 00582 aveUTime[i] = 0; 00583 approxVfromTime[i] = 0; 00584 maxPhysU[i] = -1; maxPlaneU[i] = -1; 00585 minPhysU[i] = 999; minPlaneU[i] = -1; 00586 nSubShowersU[i] = 0; 00587 } 00588 00589 std::vector<TSpecReal_t> totVEnergy(nVPeaks,0); 00590 std::vector<Double_t> aveVTime(nVPeaks,0); 00591 std::vector<TSpecReal_t> approxUfromTime(nVPeaks,0); 00592 std::vector<TSpecReal_t> maxPhysV(nVPeaks,0); 00593 std::vector<TSpecReal_t> minPhysV(nVPeaks,0); 00594 std::vector<Int_t> maxPlaneV(nVPeaks,0); 00595 std::vector<Int_t> minPlaneV(nVPeaks,0); 00596 std::vector<Int_t> nSubShowersV(nVPeaks,0); 00597 for(int i=0;i<nVPeaks;i++) { 00598 totVEnergy[i] = 0; 00599 aveVTime[i] = 0; 00600 approxUfromTime[i] = 0; 00601 maxPhysV[i] = -1; maxPlaneV[i] = -1; 00602 minPhysV[i] = 999; minPlaneV[i] = -1; 00603 nSubShowersV[i] = 0; 00604 } 00605 00606 //loop through subshowerlist and add up energy of showers 00607 subshowerItr.ResetFirst(); 00608 while (CandSubShowerSRHandle *subshower = subshowerItr()) { 00609 if (*subshower->GetCandSlice()==*slice) { 00610 if(subshower->GetBegPlane()>=begPlane && 00611 subshower->GetEndPlane()<=endPlane) { 00612 if(subshower->GetPlaneView()==PlaneView::kU || 00613 subshower->GetPlaneView()==PlaneView::kX){ 00614 for(int i=0;i<nUPeaks;i++){ 00615 TSpecReal_t fracPHInRange = 0; 00616 TSpecReal_t totPHInSS = 0; 00617 CandStripHandleItr stripssItr(subshower->GetDaughterIterator()); 00618 while (CandStripHandle *stp = stripssItr()) { 00619 totPHInSS += stp->GetCharge(CalDigitType::kSigCorr); 00620 if(stp->GetTPos()>peakUPosLow[i] && 00621 stp->GetTPos()<=peakUPosUpp[i]) 00622 fracPHInRange+=stp->GetCharge(CalDigitType::kSigCorr); 00623 } 00624 if(totPHInSS>=0 && fracPHInRange/totPHInSS >= 0.5) { 00625 totUEnergy[i] += subshower->GetEnergy(); 00626 aveUTime[i] += 1e9*subshower->GetAveTime()*subshower->GetEnergy(); 00627 if(nUPeaks>1 && 00628 subshower->GetClusterID()==ClusterType::kEMLike || 00629 subshower->GetClusterID()==ClusterType::kHadLike || 00630 subshower->GetClusterID()==ClusterType::kTrkLike){ 00631 approxVfromTime[i] += subshower->GetVtxV(); 00632 nSubShowersU[i] += 1; 00633 for(int j=subshower->GetBegPlane();j<=subshower->GetEndPlane();j++){ 00634 TSpecReal_t testMinU = subshower->GetMinU(j); 00635 TSpecReal_t testMaxU = subshower->GetMaxU(j); 00636 if(testMaxU>maxPhysU[i]) { 00637 maxPhysU[i] = testMaxU; 00638 maxPlaneU[i] = j; 00639 } 00640 if(testMinU<minPhysU[i]) { 00641 minPhysU[i] = testMinU; 00642 minPlaneU[i] = j; 00643 } 00644 } 00645 } 00646 } 00647 } 00648 } 00649 else if(subshower->GetPlaneView()==PlaneView::kV || 00650 subshower->GetPlaneView()==PlaneView::kY){ 00651 for(int i=0;i<nVPeaks;i++){ 00652 TSpecReal_t fracPHInRange = 0; 00653 TSpecReal_t totPHInSS = 0; 00654 CandStripHandleItr stripssItr(subshower->GetDaughterIterator()); 00655 while (CandStripHandle *stp = stripssItr()) { 00656 totPHInSS += stp->GetCharge(CalDigitType::kSigCorr); 00657 if(stp->GetTPos()>peakVPosLow[i] && 00658 stp->GetTPos()<=peakVPosUpp[i]) 00659 fracPHInRange+=stp->GetCharge(CalDigitType::kSigCorr); 00660 } 00661 if(totPHInSS>=0 && fracPHInRange/totPHInSS >= 0.5) { 00662 totVEnergy[i] += subshower->GetEnergy(); 00663 aveVTime[i] += 1e9*subshower->GetAveTime()*subshower->GetEnergy(); 00664 if(nVPeaks>1 && 00665 subshower->GetClusterID()==ClusterType::kEMLike || 00666 subshower->GetClusterID()==ClusterType::kHadLike || 00667 subshower->GetClusterID()==ClusterType::kTrkLike){ 00668 approxUfromTime[i] += subshower->GetVtxU(); 00669 nSubShowersV[i] += 1; 00670 for(int j=subshower->GetBegPlane();j<=subshower->GetEndPlane();j++){ 00671 TSpecReal_t testMinV = subshower->GetMinV(j); 00672 TSpecReal_t testMaxV = subshower->GetMaxV(j); 00673 if(testMaxV>maxPhysV[i]) { 00674 maxPhysV[i] = testMaxV; 00675 maxPlaneV[i] = j; 00676 } 00677 if(testMinV<minPhysV[i]) { 00678 minPhysV[i] = testMinV; 00679 minPlaneV[i] = j; 00680 } 00681 } 00682 } 00683 } 00684 } 00685 } 00686 } 00687 } 00688 } 00689 00690 00691 //define number of showers there will be in this plane range 00692 Int_t nShowers = nUPeaks; 00693 if(nShowers>nVPeaks) nShowers = nVPeaks; 00694 MSG("AlgShowerSS",Msg::kDebug) << "nShowers = " << nShowers << endl; 00695 00696 00697 //Combine peaks if there is a U/V mismatch in number of peaks 00698 //U planes: 00699 Int_t sanity_counter = 0; //just in case it decides to keep on looping... 00700 Int_t nUPeaks_tmp = nUPeaks; 00701 while(nUPeaks_tmp!=nShowers && sanity_counter<3){ 00702 Int_t merge_peak_1 = -1; //lowest merge peak 00703 Int_t merge_peak_2 = -1; //highest merge peak 00704 TSpecReal_t best_merge_metric = 999999; 00705 for(int i=0;i<nUPeaks-1;i++){ 00706 //don't care about empty peak regions: 00707 if(totUEnergy[i]<=0) continue; 00708 MSG("AlgShowerSS",Msg::kDebug) << "Peak " << i << " is not empty" << endl; 00709 //find next peak region that isn't empty: 00710 Int_t comp_index = 1; 00711 for (; (i+comp_index)<nUPeaks; ++comp_index) 00712 if (totUEnergy[i+comp_index]>0) break; 00713 /* 00714 while( totUEnergy[i+comp_index]<=0 && i+comp_index<nUPeaks) { 00715 comp_index += 1; 00716 } 00717 */ 00718 00719 if(i+comp_index>=nUPeaks) break; //no more comparisons available 00720 MSG("AlgShowerSS",Msg::kDebug) << "Compare with peak " 00721 << i+comp_index << endl; 00722 00723 TSpecReal_t merge_metric = 999999; 00724 TSpecReal_t time_diff = 999999; 00725 if(totUEnergy[i+comp_index]>0 && minPlaneU[i+comp_index]!=-1) { 00726 //get distance of closest approach in metres: 00727 merge_metric = TMath::Power(minPhysU[i+comp_index] - maxPhysU[i],2); 00728 merge_metric += PITCHINMETRES*TMath::Power(minPlaneU[i+comp_index] - 00729 maxPlaneU[i],2); 00730 merge_metric = TMath::Sqrt(merge_metric); 00731 MSG("AlgShowerSS",Msg::kDebug) << "Distance metric = " 00732 << merge_metric << endl; 00733 00734 //add in effective distance based on timing: 00735 merge_metric += Mphysical::c_light*1e-9*TMath::Abs(aveUTime[i+comp_index] / 00736 totUEnergy[i+comp_index] - 00737 aveUTime[i]/totUEnergy[i]); 00738 MSG("AlgShowerSS",Msg::kDebug) 00739 << "Time metric = " << Mphysical::c_light*1e-9* 00740 TMath::Abs(aveUTime[i+comp_index] / totUEnergy[i+comp_index] - 00741 aveUTime[i] / totUEnergy[i]) << endl; 00742 time_diff = ( TMath::Abs(aveUTime[i+comp_index] / 00743 totUEnergy[i+comp_index] - 00744 aveUTime[i] / totUEnergy[i] ) ); 00745 } 00746 else if(totUEnergy[i+comp_index]>0) { 00747 // just check timing since this is a non-physics subshower cluster 00748 // this will lead to a smaller merge_metric value in general 00749 // but it is probably better to merge the non-physics clusters in any case 00750 merge_metric = Mphysical::c_light*1e-9*TMath::Abs(aveUTime[i+comp_index] / 00751 totUEnergy[i+comp_index] - 00752 aveUTime[i]/totUEnergy[i]); 00753 time_diff = ( TMath::Abs(aveUTime[i+comp_index] / 00754 totUEnergy[i+comp_index] - 00755 aveUTime[i] / totUEnergy[i] ) ); 00756 } 00757 00758 MSG("AlgShowerSS",Msg::kDebug) << "Final metric = " << merge_metric << endl; 00759 if(merge_metric<best_merge_metric && time_diff<PeakMergeTimeCut){ 00760 best_merge_metric = merge_metric; 00761 merge_peak_1 = i; 00762 merge_peak_2 = i+comp_index; 00763 MSG("AlgShowerSS",Msg::kDebug) << "New best merge peak!\n" 00764 << "Merge peak " << merge_peak_1 00765 << " with " << merge_peak_2 00766 << " based on metric value " 00767 << best_merge_metric << endl; 00768 } 00769 } 00770 00771 if(best_merge_metric<999999 && merge_peak_1!=-1 && merge_peak_2!=-1){ 00772 //add energy to second peak and 00773 //remove energy from first peak 00774 totUEnergy[merge_peak_2] += totUEnergy[merge_peak_1]; 00775 totUEnergy[merge_peak_1] = 0; 00776 aveUTime[merge_peak_2] += aveUTime[merge_peak_1]; 00777 aveUTime[merge_peak_1] = 0; 00778 //add other-view vertex estimate sum 00779 approxVfromTime[merge_peak_2] += approxVfromTime[merge_peak_1]; 00780 approxVfromTime[merge_peak_1] = 0; 00781 00782 if(nSubShowersU[merge_peak_1]>=0 && nSubShowersU[merge_peak_2]>=0) { 00783 //shift min U of second peak down and 00784 //set min U plane of second peak 00785 minPhysU[merge_peak_2] = minPhysU[merge_peak_1]; 00786 minPlaneU[merge_peak_2] = minPlaneU[merge_peak_1]; 00787 //set min/max limits of first peak to be the same 00788 maxPhysU[merge_peak_1] = minPhysU[merge_peak_1]; 00789 maxPlaneU[merge_peak_1] = minPlaneU[merge_peak_1]; 00790 //set peak lower limit of second peak to that of first 00791 peakUPosLow[merge_peak_2] = peakUPosLow[merge_peak_1]; 00792 //set peak lower and upper limit of first peak to be the same 00793 peakUPosUpp[merge_peak_1] = peakUPosLow[merge_peak_1]; 00794 } 00795 else if(nSubShowersU[merge_peak_1]==-1) { 00796 //first peak is non-physical so don't change anything for merge_peak_2 00797 //set min/max limits of first peak to be the same 00798 maxPhysU[merge_peak_1] = minPhysU[merge_peak_1]; 00799 maxPlaneU[merge_peak_1] = minPlaneU[merge_peak_1]; 00800 //set peak lower and upper limit of first peak to be the same 00801 peakUPosUpp[merge_peak_1] = peakUPosLow[merge_peak_1]; 00802 } 00803 else if(nSubShowersU[merge_peak_2]==-1) { 00804 //set min/max limits of second peak to be the same as the first peak 00805 minPhysU[merge_peak_2] = minPhysU[merge_peak_1]; 00806 minPlaneU[merge_peak_2] = minPlaneU[merge_peak_1]; 00807 maxPhysU[merge_peak_2] = maxPhysU[merge_peak_1]; 00808 maxPlaneU[merge_peak_2] = maxPlaneU[merge_peak_1]; 00809 //set min/max limits of first peak to be the same 00810 maxPhysU[merge_peak_1] = minPhysU[merge_peak_1]; 00811 maxPlaneU[merge_peak_1] = minPlaneU[merge_peak_1]; 00812 //set peak lower limit of second peak to that of first 00813 peakUPosLow[merge_peak_2] = peakUPosLow[merge_peak_1]; 00814 //set peak lower and upper limit of first peak to be the same 00815 peakUPosUpp[merge_peak_1] = peakUPosLow[merge_peak_1]; 00816 00817 } 00818 nSubShowersU[merge_peak_2] += nSubShowersU[merge_peak_1]; 00819 nSubShowersU[merge_peak_1] = 0; 00820 nUPeaks_tmp -= 1; 00821 } 00822 else { 00823 MSG("AlgShowerSS",Msg::kDebug) 00824 << "Incrementing U view sanity_counter:\n" 00825 << " nUPeaks_tmp = " << nUPeaks_tmp 00826 << " best_merge_metric = " << best_merge_metric 00827 << " merge_peak_1 = " << merge_peak_1 00828 << " merge_peak_2 = " << merge_peak_2 00829 << endl; 00830 sanity_counter+=1; 00831 } 00832 } 00833 00834 //V planes: 00835 Int_t nVPeaks_tmp = nVPeaks; 00836 while(nVPeaks_tmp!=nShowers && sanity_counter<5){ 00837 Int_t merge_peak_1 = -1; //lowest merge peak 00838 Int_t merge_peak_2 = -1; //highest merge peak 00839 TSpecReal_t best_merge_metric = 999999; 00840 for(int i=0;i<nVPeaks-1;i++){ 00841 //don't care about empty peak regions: 00842 if(totVEnergy[i]<=0) continue; 00843 MSG("AlgShowerSS",Msg::kDebug) << "Peak " << i << " is not empty" << endl; 00844 //find next peak region that isn't empty: 00845 Int_t comp_index = 1; 00846 for (; (i+comp_index)<nVPeaks; ++comp_index) 00847 if (totVEnergy[i+comp_index]>0) break; 00848 /* 00849 while( totVEnergy[i+comp_index]<=0 && i+comp_index<nVPeaks) { 00850 comp_index += 1; 00851 } 00852 */ 00853 if(i+comp_index>=nVPeaks) break; //no more comparisons available 00854 MSG("AlgShowerSS",Msg::kDebug) << "Compare with peak " 00855 << i+comp_index << endl; 00856 00857 TSpecReal_t merge_metric = 999999; 00858 if(totVEnergy[i+comp_index]>0 && minPlaneV[i+comp_index]!=-1) { 00859 //get distance of closest approach in metres: 00860 merge_metric = TMath::Power(minPhysV[i+comp_index] - maxPhysV[i],2); 00861 merge_metric += PITCHINMETRES*TMath::Power(minPlaneV[i+comp_index] - 00862 maxPlaneV[i],2); 00863 merge_metric = TMath::Sqrt(merge_metric); 00864 MSG("AlgShowerSS",Msg::kDebug) << "Distance metric = " 00865 << merge_metric << endl; 00866 00867 //add in effective distance based on timing: 00868 merge_metric += Mphysical::c_light*1e-9*TMath::Abs(aveVTime[i+comp_index] / 00869 totVEnergy[i+comp_index] - 00870 aveVTime[i]/totVEnergy[i]); 00871 MSG("AlgShowerSS",Msg::kDebug) 00872 << "Time metric = " << Mphysical::c_light*1e-9* 00873 TMath::Abs(aveVTime[i+comp_index] / totVEnergy[i+comp_index] - 00874 aveVTime[i] / totVEnergy[i]) << endl; 00875 } 00876 else if(totVEnergy[i+comp_index]>0) { 00877 // just check timing since this is a non-physics subshower cluster 00878 // this will lead to a smaller merge_metric value in general 00879 // but it is probably better to merge the non-physics clusters in any case 00880 merge_metric = Mphysical::c_light*1e-9*TMath::Abs(aveVTime[i+comp_index] / 00881 totVEnergy[i+comp_index] - 00882 aveVTime[i]/totVEnergy[i]); 00883 } 00884 00885 MSG("AlgShowerSS",Msg::kDebug) << "Final metric = " << merge_metric << endl; 00886 if(merge_metric<best_merge_metric){ 00887 best_merge_metric = merge_metric; 00888 merge_peak_1 = i; 00889 merge_peak_2 = i+comp_index; 00890 MSG("AlgShowerSS",Msg::kDebug) << "New best merge peak!\n" 00891 << "Merge peak " << merge_peak_1 00892 << " with " << merge_peak_2 00893 << " based on metric value " 00894 << best_merge_metric << endl; 00895 } 00896 } 00897 00898 if(best_merge_metric<999999 && merge_peak_1!=-1 && merge_peak_2!=-1){ 00899 //add energy to second peak and 00900 //remove energy from first peak 00901 totVEnergy[merge_peak_2] += totVEnergy[merge_peak_1]; 00902 totVEnergy[merge_peak_1] = 0; 00903 aveVTime[merge_peak_2] += aveVTime[merge_peak_1]; 00904 aveVTime[merge_peak_1] = 0; 00905 //add other-view vertex estimate sum 00906 approxUfromTime[merge_peak_2] += approxUfromTime[merge_peak_1]; 00907 approxUfromTime[merge_peak_1] = 0; 00908 00909 if(nSubShowersV[merge_peak_1]>=0 && nSubShowersV[merge_peak_2]>=0) { 00910 //shift min V of second peak down and 00911 //set min V plane of second peak 00912 minPhysV[merge_peak_2] = minPhysV[merge_peak_1]; 00913 minPlaneV[merge_peak_2] = minPlaneV[merge_peak_1]; 00914 //set min/max limits of first peak to be the same 00915 maxPhysV[merge_peak_1] = minPhysV[merge_peak_1]; 00916 maxPlaneV[merge_peak_1] = minPlaneV[merge_peak_1]; 00917 //set peak lower limit of second peak to that of first 00918 peakVPosLow[merge_peak_2] = peakVPosLow[merge_peak_1]; 00919 //set peak lower and upper limit of first peak to be the same 00920 peakVPosUpp[merge_peak_1] = peakVPosLow[merge_peak_1]; 00921 } 00922 else if(nSubShowersV[merge_peak_1]==-1) { 00923 //first peak is non-physical so don't change anything for merge_peak_2 00924 //set min/max limits of first peak to be the same 00925 maxPhysV[merge_peak_1] = minPhysV[merge_peak_1]; 00926 maxPlaneV[merge_peak_1] = minPlaneV[merge_peak_1]; 00927 //set peak lower and upper limit of first peak to be the same 00928 peakVPosUpp[merge_peak_1] = peakVPosLow[merge_peak_1]; 00929 } 00930 else if(nSubShowersV[merge_peak_2]==-1) { 00931 //set min/max limits of second peak to be the same as the first peak 00932 minPhysV[merge_peak_2] = minPhysV[merge_peak_1]; 00933 minPlaneV[merge_peak_2] = minPlaneV[merge_peak_1]; 00934 maxPhysV[merge_peak_2] = maxPhysV[merge_peak_1]; 00935 maxPlaneV[merge_peak_2] = maxPlaneV[merge_peak_1]; 00936 //set min/max limits of first peak to be the same 00937 maxPhysV[merge_peak_1] = minPhysV[merge_peak_1]; 00938 maxPlaneV[merge_peak_1] = minPlaneV[merge_peak_1]; 00939 //set peak lower limit of second peak to that of first 00940 peakVPosLow[merge_peak_2] = peakVPosLow[merge_peak_1]; 00941 //set peak lower and upper limit of first peak to be the same 00942 peakVPosUpp[merge_peak_1] = peakVPosLow[merge_peak_1]; 00943 00944 } 00945 nSubShowersV[merge_peak_2] += nSubShowersV[merge_peak_1]; 00946 nSubShowersV[merge_peak_1] = 0; 00947 nVPeaks_tmp -= 1; 00948 } 00949 else { 00950 MSG("AlgShowerSS",Msg::kDebug) 00951 << "Incrementing V view sanity_counter:\n" 00952 << " nVPeaks_tmp = " << nVPeaks_tmp 00953 << " best_merge_metric = " << best_merge_metric 00954 << " merge_peak_1 = " << merge_peak_1 00955 << " merge_peak_2 = " << merge_peak_2 00956 << endl; 00957 sanity_counter+=1; 00958 } 00959 } 00960 00961 for(int i=0;i<nShowers;i++){ 00962 MSG("AlgShowerSS",Msg::kDebug) << "Positions, peaks and limits after merging:" 00963 << endl; 00964 MSG("AlgShowerSS",Msg::kDebug) << peakUPos[i] << " " 00965 << peakUVal[i] << " [" 00966 << peakUPosLow[i] << "," 00967 << peakUPosUpp[i] << "]" 00968 << endl; 00969 MSG("AlgShowerSS",Msg::kDebug) << peakVPos[i] << " " 00970 << peakVVal[i] << " [" 00971 << peakVPosLow[i] << "," 00972 << peakVPosUpp[i] << "]" 00973 << endl; 00974 } 00975 00976 00977 //get best order for matching subshower groups in U/V 00978 std::vector<Int_t> bestOrderU(nUPeaks,0); 00979 std::vector<Int_t> bestOrderV(nVPeaks,0); 00980 for(int i=0;i<nUPeaks;i++){ 00981 if(totUEnergy[i]>0){ 00982 Int_t nbigger = 0; 00983 for(int j=0;j<nUPeaks;j++){ 00984 if(totUEnergy[i]<totUEnergy[j]) nbigger+=1; 00985 } 00986 bestOrderU[nbigger] = i; 00987 } 00988 else { // if there is no energy in this "peak" 00989 Int_t nsmaller = 0; // put zero peaks in reverse order 00990 for(int j=0;j<i;j++){ // purely to ensure that all indices 00991 if(totUEnergy[j]==0) nsmaller+=1; // of bestOrderX are filled 00992 } // with legitimate values (prevents 00993 bestOrderU[nUPeaks-1-nsmaller] = i; // against a rare case seg fault) 00994 } 00995 } 00996 for(int i=0;i<nVPeaks;i++){ 00997 if(totVEnergy[i]>0){ 00998 Int_t nbigger = 0; 00999 for(int j=0;j<nVPeaks;j++){ 01000 if(totVEnergy[i]<totVEnergy[j]) nbigger+=1; 01001 } 01002 bestOrderV[nbigger] = i; 01003 } 01004 else { 01005 Int_t nsmaller = 0; 01006 for(int j=0;j<i;j++){ 01007 if(totVEnergy[j]==0) nsmaller+=1; 01008 } 01009 bestOrderV[nVPeaks-1-nsmaller] = i; 01010 } 01011 } 01012 01013 //check energy ordering is ok: 01014 for(int i=0;i<nShowers-1;i++){ //loop over showers 01015 TSpecReal_t EUi = totUEnergy[bestOrderU[i]]; 01016 TSpecReal_t EVi = totVEnergy[bestOrderV[i]]; 01017 TSpecReal_t EUii = totUEnergy[bestOrderU[i+1]]; 01018 TSpecReal_t EVii = totVEnergy[bestOrderV[i+1]]; 01019 if(EUi==0 || EUii==0 || EVi==0 || EVii==0) continue; 01020 //can we swap the order and still get two good energy matches? 01021 TSpecReal_t asym = TMath::Abs(EUi - EVi)/(EUi + EVi) + 01022 TMath::Abs(EUii - EVii)/(EUii + EVii); 01023 TSpecReal_t asymSwap = TMath::Abs(EUi - EVii)/(EUi + EVii) + 01024 TMath::Abs(EUii - EVi)/(EUii + EVi); 01025 //if difference between two asymmetry estimates is <30% 01026 //see if vertex estimate from timing 01027 //or if plane range information helps 01028 if( TMath::Abs(asym - asymSwap)/asym < EnergyAsymmetryCut || 01029 TMath::Abs(aveUTime[bestOrderU[i]]/EUi - 01030 aveVTime[bestOrderV[i]]/EVi) > SubShowerTimingCut ) { 01031 if(nSubShowersU[bestOrderU[i]]==0 || 01032 nSubShowersV[bestOrderV[i]]==0 || 01033 nSubShowersU[bestOrderU[i+1]]==0 || 01034 nSubShowersV[bestOrderV[i+1]]==0) continue; 01035 01036 Double_t timediff = 01037 TMath::Abs( (aveUTime[bestOrderU[i]]/EUi) - 01038 (aveVTime[bestOrderV[i]]/EVi) ) + 01039 TMath::Abs( (aveUTime[bestOrderU[i+1]]/EUii) - 01040 (aveVTime[bestOrderV[i+1]]/EVii) ); 01041 Double_t timediff_swap = 01042 TMath::Abs( (aveUTime[bestOrderU[i]]/EUi) - 01043 (aveVTime[bestOrderV[i+1]]/EVii) ) + 01044 TMath::Abs( (aveUTime[bestOrderU[i+1]]/EUii) - 01045 (aveVTime[bestOrderV[i]]/EVi) ); 01046 01047 if(timediff_swap<timediff && 01048 1.0e9*timediff_swap<SubShowerTimingCut) { //check that timing diff is reasonable 01049 //swap! 01050 //arbitrary choice which view to switch: 01051 Int_t tempIndex = bestOrderU[i]; 01052 bestOrderU[i] = bestOrderU[i+1]; 01053 bestOrderU[i+1] = tempIndex; 01054 //start checks again so that clusters 01055 //can propagate if better matches are found 01056 i=0; 01057 continue; 01058 } 01059 01060 //now check plane range agreement: 01061 Int_t planediff = 01062 TMath::Abs(minPlaneU[bestOrderU[i]] - minPlaneV[bestOrderV[i]]) + 01063 TMath::Abs(maxPlaneU[bestOrderU[i]] - maxPlaneV[bestOrderV[i]]) + 01064 TMath::Abs(minPlaneU[bestOrderU[i+1]] - minPlaneV[bestOrderV[i+1]]) + 01065 TMath::Abs(maxPlaneU[bestOrderU[i+1]] - maxPlaneV[bestOrderV[i+1]]); 01066 01067 Int_t planediff_swap = 01068 TMath::Abs(minPlaneU[bestOrderU[i]] - minPlaneV[bestOrderV[i+1]]) + 01069 TMath::Abs(maxPlaneU[bestOrderU[i]] - maxPlaneV[bestOrderV[i+1]]) + 01070 TMath::Abs(minPlaneU[bestOrderU[i+1]] - minPlaneV[bestOrderV[i]]) + 01071 TMath::Abs(maxPlaneU[bestOrderU[i+1]] - maxPlaneV[bestOrderV[i]]); 01072 01073 if(planediff_swap<planediff && 01074 1.0e9*timediff_swap<SubShowerTimingCut) { 01075 //swap! 01076 //arbitrary choice which view to switch: 01077 Int_t tempIndex = bestOrderU[i]; 01078 bestOrderU[i] = bestOrderU[i+1]; 01079 bestOrderU[i+1] = tempIndex; 01080 //start checks again so that clusters 01081 //can propagate if better matches are found 01082 i=0; 01083 continue; 01084 } 01085 01086 if(false){ //space + time - not implemented yet 01087 approxVfromTime[bestOrderU[i]] /= nSubShowersU[bestOrderU[i]]; 01088 approxUfromTime[bestOrderV[i]] /= nSubShowersV[bestOrderV[i]]; 01089 approxVfromTime[bestOrderU[i+1]] /= nSubShowersU[bestOrderU[i+1]]; 01090 approxUfromTime[bestOrderV[i+1]] /= nSubShowersV[bestOrderV[i+1]]; 01091 01092 //sum distance from both showers 01093 TSpecReal_t vertexDistance = 01094 TMath::Sqrt(TMath::Power(approxUfromTime[bestOrderV[i]] - 01095 peakUPos[bestOrderU[i]],2.) + 01096 TMath::Power(approxVfromTime[bestOrderU[i]] - 01097 peakVPos[bestOrderV[i]],2)) + 01098 TMath::Sqrt(TMath::Power(approxUfromTime[bestOrderV[i+1]] - 01099 peakUPos[bestOrderU[i+1]],2.) + 01100 TMath::Power(approxVfromTime[bestOrderU[i+1]] - 01101 peakVPos[bestOrderV[i+1]],2)); 01102 TSpecReal_t vertexDistanceSwap = 01103 TMath::Sqrt(TMath::Power(approxUfromTime[bestOrderV[i+1]] - 01104 peakUPos[bestOrderU[i]],2.) + 01105 TMath::Power(approxVfromTime[bestOrderU[i]] - 01106 peakVPos[bestOrderV[i+1]],2)) + 01107 TMath::Sqrt(TMath::Power(approxUfromTime[bestOrderV[i]] - 01108 peakUPos[bestOrderU[i+1]],2.) + 01109 TMath::Power(approxVfromTime[bestOrderU[i+1]] - 01110 peakVPos[bestOrderV[i]],2)); 01111 01112 if(vertexDistanceSwap<vertexDistance && 01113 vertexDistanceSwap<0.5) { //check that best vertex is reasonable 01114 //swap! 01115 //arbitrary choice which view to switch: 01116 Int_t tempIndex = bestOrderU[i]; 01117 bestOrderU[i] = bestOrderU[i+1]; 01118 bestOrderU[i+1] = tempIndex; 01119 //start checks again so that clusters 01120 //can propagate if better matches are found 01121 i=0; 01122 } 01123 } 01124 01125 } 01126 } 01127 01128 //loop over number of expected showers: 01129 for(int i=0;i<nShowers;i++){ 01130 TObjArray newshower; 01131 Int_t nsubshower = 0; 01132 subshowerItr.ResetFirst(); 01133 while (CandSubShowerSRHandle *subshower = subshowerItr()) { 01134 if (*subshower->GetCandSlice()==*slice) { 01135 if(subshower->GetBegPlane()>=begPlane && 01136 subshower->GetEndPlane()<=endPlane) { 01137 if(subshower->GetPlaneView()==PlaneView::kU || 01138 subshower->GetPlaneView()==PlaneView::kX){ 01139 TSpecReal_t fracPHInRange = 0; 01140 TSpecReal_t totPHInSS = 0; 01141 CandStripHandleItr stripssItr(subshower->GetDaughterIterator()); 01142 while (CandStripHandle *stp = stripssItr()) { 01143 totPHInSS += stp->GetCharge(CalDigitType::kSigCorr); 01144 if(stp->GetTPos()>peakUPosLow[bestOrderU[i]] && 01145 stp->GetTPos()<=peakUPosUpp[bestOrderU[i]]) 01146 fracPHInRange+=stp->GetCharge(CalDigitType::kSigCorr); 01147 } 01148 if(totPHInSS>=0 && fracPHInRange/totPHInSS >= 0.5) { 01149 newshower.Add(subshower); 01150 MSG("AlgShowerSS",Msg::kDebug) 01151 << "Adding SubShower: " 01152 << nsubshower << " in slice " << nslice 01153 << " to newshower " << i << endl; 01154 } 01155 } 01156 if(subshower->GetPlaneView()==PlaneView::kV || 01157 subshower->GetPlaneView()==PlaneView::kY){ 01158 TSpecReal_t fracPHInRange = 0; 01159 TSpecReal_t totPHInSS = 0; 01160 CandStripHandleItr stripssItr(subshower->GetDaughterIterator()); 01161 while (CandStripHandle *stp = stripssItr()) { 01162 totPHInSS += stp->GetCharge(CalDigitType::kSigCorr); 01163 if(stp->GetTPos()>peakVPosLow[bestOrderV[i]] && 01164 stp->GetTPos()<=peakVPosUpp[bestOrderV[i]]) 01165 fracPHInRange+=stp->GetCharge(CalDigitType::kSigCorr); 01166 } 01167 if(totPHInSS>=0 && fracPHInRange/totPHInSS >= 0.5) { 01168 newshower.Add(subshower); 01169 MSG("AlgShowerSS",Msg::kDebug) 01170 << "Adding SubShower: " 01171 << nsubshower << " in slice " << nslice 01172 << " to newshower " << i << endl; 01173 } 01174 } 01175 } 01176 } 01177 nsubshower+=1; 01178 } 01179 01180 //now check that max/min strip times from subshowers do not have gaps 01181 //greater than 50ns 01182 //make vectors to hold min/max times per subshower 01183 //then loop through other subshowers looking for gaps>100ns 01184 //find "contiguous" time range and add up the energy 01185 //will end up with duplicate entries in the vectors since most subshowers 01186 //should be contiguous 01187 Int_t nAddedSS = newshower.GetLast()+1; 01188 vector<Double_t> contiguousMin(nAddedSS,0); 01189 vector<Double_t> contiguousMax(nAddedSS,0); 01190 vector<Double_t> totalEnergyInRange(nAddedSS,0); 01191 for(int j=0;j<nAddedSS;j++){ 01192 CandSubShowerSRHandle *subshower = 01193 dynamic_cast<CandSubShowerSRHandle*>(newshower.At(j)); 01194 contiguousMin[j] = 1e9*subshower->GetMinStpTime(); 01195 contiguousMax[j] = 1e9*subshower->GetMaxStpTime(); 01196 MSG("AlgShowerSS",Msg::kDebug) << "Initial " << j << " " 01197 << contiguousMin[j] 01198 << " " << contiguousMax[j] << endl; 01199 totalEnergyInRange[j] = subshower->GetEnergy(); 01200 for(int k=0;k<nAddedSS;k++){ 01201 if(j==k) continue; 01202 CandSubShowerSRHandle *subshower2 = 01203 dynamic_cast<CandSubShowerSRHandle*>(newshower.At(k)); 01204 Double_t max2 = 1e9*subshower2->GetMaxStpTime(); 01205 Double_t min2 = 1e9*subshower2->GetMinStpTime(); 01206 if(min2 < contiguousMin[j] && max2+SubShowerTimingCut < contiguousMin[j]) continue; 01207 if(min2-SubShowerTimingCut > contiguousMax[j] && max2 > contiguousMax[j]) continue; 01208 if(min2 > contiguousMin[j] && max2 < contiguousMax[j]) continue; 01209 if(min2 < contiguousMin[j] && max2 > contiguousMax[j]){ 01210 contiguousMin[j] = min2; 01211 contiguousMax[j] = max2; 01212 totalEnergyInRange[j] += subshower2->GetEnergy(); 01213 continue; 01214 } 01215 if(max2 <= contiguousMax[j] && 01216 max2+SubShowerTimingCut >= contiguousMin[j] && 01217 min2 < contiguousMin[j]) { 01218 contiguousMin[j] = min2; 01219 totalEnergyInRange[j] += subshower2->GetEnergy(); 01220 continue; 01221 } 01222 if(min2 >= contiguousMin[j] && 01223 min2-SubShowerTimingCut <= contiguousMax[j] && 01224 max2 > contiguousMax[j]) { 01225 contiguousMax[j] = max2; 01226 totalEnergyInRange[j] += subshower2->GetEnergy(); 01227 continue; 01228 } 01229 } 01230 MSG("AlgShowerSS",Msg::kDebug) << "Final " << j << " " 01231 << contiguousMin[j] 01232 << " " << contiguousMax[j] << endl; 01233 } 01234 01235 //find largest energy range: 01236 Int_t rangeIndex = -1; 01237 Double_t largestEnergy = 0; 01238 for(int j=0;j<nAddedSS;j++){ 01239 if(totalEnergyInRange[j]>largestEnergy){ 01240 rangeIndex = j; 01241 largestEnergy = totalEnergyInRange[j]; 01242 } 01243 } 01244 if(nAddedSS>0) MSG("AlgShowerSS",Msg::kDebug) << "Largest " << rangeIndex 01245 << " " << largestEnergy << " " 01246 << contiguousMin[rangeIndex] 01247 << " " << contiguousMax[rangeIndex] 01248 << endl; 01249 01250 //remove subshowers out of this range: 01251 for(int j=0;j<nAddedSS;j++){ 01252 CandSubShowerSRHandle *subshower = 01253 dynamic_cast<CandSubShowerSRHandle*>(newshower.At(j)); 01254 if(1e9*subshower->GetMinStpTime()<(contiguousMin[rangeIndex]-0.00001) || 01255 1e9*subshower->GetMaxStpTime()>(contiguousMax[rangeIndex]+0.00001)) { 01256 newshower.RemoveAt(j); 01257 MSG("AlgShowerSS",Msg::kDebug) << "Removing SubShower: " 01258 << j << " in slice " << nslice 01259 << " from newshower " << i << endl; 01260 } 01261 } 01262 newshower.Compress(); 01263 01264 cxx.SetDataIn(&newshower); 01265 CandShowerSRHandle showerhandle = CandShowerSR::MakeCandidate(ah,cxx); 01266 if((showerhandle.GetEnergy()>0 && nUPeaks_save>0 && nVPeaks_save>0) || 01267 showerhandle.GetEnergy()>MinNonPhysicsEnergy) { 01268 showerhandle.SetCandSlice(slice); 01269 ch.AddDaughterLink(showerhandle); 01270 MSG("AlgShowerSS",Msg::kDebug) << "Slice " << nslice << " Shower "<< nshower 01271 << " uid = " << showerhandle.GetUidInt() << endl; 01272 } 01273 else { 01274 MSG("AlgShowerSS",Msg::kDebug) << "Shower " << nshower 01275 << " has zero energy." 01276 << " - Not adding shower to list." 01277 << endl; 01278 } 01279 } 01280 } 01281 nslice++; 01282 } 01283 }
| void AlgShowerSSList::Trace | ( | const char * | c | ) | const [virtual] |
TH1F* AlgShowerSSList::vtxUHist [private] |
TH1F* AlgShowerSSList::vtxUHistAll [private] |
TH1F* AlgShowerSSList::vtxVHist [private] |
TH1F* AlgShowerSSList::vtxVHistAll [private] |
1.4.7