AlgShowerSSList Class Reference

#include <AlgShowerSSList.h>

Inheritance diagram for AlgShowerSSList:

AlgShowerSRList AlgBase List of all members.

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

Detailed Description

Definition at line 15 of file AlgShowerSSList.h.


Constructor & Destructor Documentation

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 }


Member Function Documentation

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]

Reimplemented from AlgShowerSRList.

Definition at line 1288 of file AlgShowerSSList.cxx.

01289 {
01290 }


Member Data Documentation

TH1F* AlgShowerSSList::vtxUHist [private]

Definition at line 27 of file AlgShowerSSList.h.

Referenced by RunAlg(), and ~AlgShowerSSList().

TH1F* AlgShowerSSList::vtxUHistAll [private]

Definition at line 26 of file AlgShowerSSList.h.

Referenced by RunAlg(), and ~AlgShowerSSList().

TH1F* AlgShowerSSList::vtxVHist [private]

Definition at line 29 of file AlgShowerSSList.h.

Referenced by RunAlg(), and ~AlgShowerSSList().

TH1F* AlgShowerSSList::vtxVHistAll [private]

Definition at line 28 of file AlgShowerSSList.h.

Referenced by RunAlg(), and ~AlgShowerSSList().


The documentation for this class was generated from the following files:
Generated on Thu Jul 10 22:52:18 2014 for loon by  doxygen 1.4.7