loon: AlgShowerSSList.cxx Source File

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