AlgSubShowerSRList Class Reference

#include <AlgSubShowerSRList.h>

Inheritance diagram for AlgSubShowerSRList:
AlgBase

List of all members.

Public Member Functions

 AlgSubShowerSRList ()
virtual ~AlgSubShowerSRList ()
virtual void RunAlg (AlgConfig &ac, CandHandle &ch, CandContext &cx)
virtual void Trace (const char *c) const

Static Public Member Functions

static NavKey StripKeyFromPlane (const CandStripHandle *)

Private Member Functions

std::vector< TObjArray * > SplitSubShower (const TObjArray *subShower)
TObjArray * MakeClusters (std::vector< CandStripHandle * > &remainingHits, TObjArray *shower)
Bool_t AreNeighbours (const CandStripHandle *stp1, const CandStripHandle *stp2)
Bool_t FindCluster (TObjArray *, TObjArray *, Int_t)
Bool_t TestOverLap (TObjArray *, CandHandle &, const CandSliceHandle *)
Bool_t TimeTest (TObjArray *, TObjArray *allStrips=NULL)
std::vector< windowTransCluster (std::vector< window >, Int_t)
std::vector< windowHoughTransCluster (std::vector< window >, Int_t)
Bool_t CleanUp (TObjArray *, CandHandle &, Int_t, const CandSliceHandle *cslice=NULL)
Bool_t FormHalo (TObjArray *, CandHandle &, CandContext &, AlgHandle &, const CandSliceListHandle *, Int_t)
std::vector< Double_t > BestHough (std::vector< window >, Bool_t)

Private Attributes

Double_t fLongWindowMipCut
Int_t fLongWindowGapPlaneCut
Int_t fMinHoughPlanes
Double_t fMinHoughPH
Double_t fBestHoughMaxPHFrac
Double_t fCleanUpTimeCut
Int_t fCleanUpPlaneCut
Double_t fCleanUpStripCut
Int_t fFormHalo
Int_t fHaloMaxPlaneGap
Int_t fHaloNPeakFindingBins
Double_t fMinStripPE
Double_t fMaxShowerGap
Int_t fSplitShowers
TH1F * vtxHist
TH2F * houghHist
TH2F * phHoughHist
TF1 * ssPol1

Detailed Description

Definition at line 35 of file AlgSubShowerSRList.h.


Constructor & Destructor Documentation

AlgSubShowerSRList::AlgSubShowerSRList (  ) 

Definition at line 64 of file AlgSubShowerSRList.cxx.

00064                                        :
00065   fLongWindowMipCut(0.0),fLongWindowGapPlaneCut(4),
00066   fMinHoughPlanes(3),fMinHoughPH(1.0), fBestHoughMaxPHFrac(0.75),
00067   fCleanUpTimeCut(50.0),fCleanUpPlaneCut(2),fCleanUpStripCut(2.5),
00068   fFormHalo(1),fHaloMaxPlaneGap(2),fHaloNPeakFindingBins(64),
00069   fMinStripPE(2)
00070 {
00071   vtxHist = new TH1F("vtxHist","vtxHist",1,-4.,4.);
00072   houghHist = new TH2F("houghHist","Hough Histogram",1,-1,1,1,-1,1);
00073   phHoughHist = new TH2F("phHoughHist","PH Weighted Hough Histogram",1,-1,1,1,-1,1);
00074   ssPol1 = new TF1("ssPol1","pol1",0,500);
00075 }

AlgSubShowerSRList::~AlgSubShowerSRList (  )  [virtual]

Definition at line 78 of file AlgSubShowerSRList.cxx.

References houghHist, phHoughHist, ssPol1, and vtxHist.

00079 {
00080   delete vtxHist;
00081   delete houghHist;
00082   delete phHoughHist;
00083   delete ssPol1;
00084 }


Member Function Documentation

Bool_t AlgSubShowerSRList::AreNeighbours ( const CandStripHandle stp1,
const CandStripHandle stp2 
) [private]

Definition at line 3160 of file AlgSubShowerSRList.cxx.

References fMaxShowerGap, CandStripHandle::GetTPos(), and CandStripHandle::GetZPos().

Referenced by MakeClusters().

03161 {
03162   //relax position cuts if in fully-instrumented region
03163   Double_t clusterRadius = fMaxShowerGap;
03164 
03165   //check for neighbouring hits in a 2D ellipsoid
03166   return (pow((stp1->GetZPos() - stp2->GetZPos())/clusterRadius,2) +
03167           pow((stp1->GetTPos() - stp2->GetTPos())/clusterRadius,2)) < 1;
03168 }

std::vector< Double_t > AlgSubShowerSRList::BestHough ( std::vector< window winVec,
Bool_t  makeEPS 
) [private]

Definition at line 2823 of file AlgSubShowerSRList.cxx.

References fBestHoughMaxPHFrac, houghHist, Msg::kDebug, Msg::kVerbose, Munits::m, MSG, phHoughHist, and STRIPWIDTHINMETRES.

Referenced by HoughTransCluster(), and TransCluster().

02825 {
02826 
02827   std::vector<window>::iterator winIterBeg = winVec.begin();
02828   std::vector<window>::iterator winIterEnd = winVec.end();
02829   std::vector<Double_t> MCV(14,0.0);
02830   for(int i=0;i<14;i++) MCV[i] = 0.0;
02831 
02832   //get extrema of windows:
02833   Int_t MaxPlane = 0;
02834   Int_t MinPlane = 500;
02835   Double_t MaxStrip = -5.;
02836   Double_t MinStrip = 5.;
02837   //total PH:
02838   Double_t totPH = 0;
02839   //PH of window with max PH:
02840   Double_t maxWinPH = 0;
02841   //plane with maximum PH, and value of max PH:
02842   Int_t maxPlanePos = -1;
02843   Double_t maxPlanePosPH = 0;
02844   //plane by plane sum PH:
02845   Double_t PlanePH[500] = {0};
02846   while(winIterBeg!=winIterEnd){
02847     if(MaxPlane<winIterBeg->plane) MaxPlane = winIterBeg->plane;
02848     if(MinPlane>winIterBeg->plane) MinPlane = winIterBeg->plane;
02849     if(MaxStrip<winIterBeg->upper_tpos) MaxStrip = winIterBeg->upper_tpos;
02850     if(MinStrip>winIterBeg->lower_tpos) MinStrip = winIterBeg->lower_tpos;
02851     if(maxWinPH<winIterBeg->ph) maxWinPH = winIterBeg->ph;
02852     totPH+=winIterBeg->ph;
02853     PlanePH[winIterBeg->plane] += winIterBeg->ph;
02854     if(PlanePH[winIterBeg->plane]>maxPlanePosPH) {
02855       maxPlanePos = winIterBeg->plane;
02856       maxPlanePosPH = PlanePH[maxPlanePos];
02857     }
02858     winIterBeg++;
02859   }
02860   
02861   //window tpos position with max PH, and value of max PH
02862   Double_t maxWindowPos = 0.;
02863   Double_t maxWindowPH = 0.;
02864   winIterBeg = winVec.begin();
02865   while(winIterBeg!=winIterEnd){
02866     if(winIterBeg->plane==maxPlanePos){
02867       if(winIterBeg->ph>maxWindowPH){
02868         maxWindowPos = winIterBeg->phpos;
02869         maxWindowPH = winIterBeg->ph;
02870       }
02871     }
02872     winIterBeg++;
02873   }
02874 
02875   if(winVec.size()<=1 || MinPlane==MaxPlane){
02876     MCV[0] = 0;
02877     MCV[1] = maxWindowPos;
02878     MCV[2] = maxPlanePos;
02879     MCV[3] = 0;
02880     MCV[4] = 0;
02881     MCV[5] = 0;
02882     MCV[6] = maxWindowPos;
02883     MCV[7] = maxPlanePos;
02884     MCV[8] = 0;
02885     MCV[9] = 0;
02886     MCV[10] = 0;
02887     MCV[11] = maxWindowPos;
02888     MCV[12] = 0;
02889     MCV[13] = maxWindowPos;
02890     return MCV;
02891   }
02892 
02893   Double_t TPosWin = MaxStrip-MinStrip+STRIPWIDTHINMETRES;
02894 
02895   MSG("SubShowerSR",Msg::kDebug) << "Hough extrema: plane " << MinPlane 
02896                                  << " to " << MaxPlane
02897                                  << " strip tpos " << MinStrip 
02898                                  << " to " << MaxStrip << endl;
02899   
02900   //define "x=0" to be at MinPlane-2, 
02901   //so y-intercept is strip value at MinPlane-2
02902   //gradient is in units of Delta(tpos) per plane
02903 
02904   houghHist->Reset();
02905   houghHist->SetBins(Int_t(TPosWin/STRIPWIDTHINMETRES) + 
02906                      MaxPlane-MinPlane+2,
02907                      -2.*TPosWin/TSpecReal_t(MaxPlane-MinPlane+2),
02908                      2.*TPosWin/TSpecReal_t(MaxPlane-MinPlane+2),
02909                      Int_t(2.*TPosWin/STRIPWIDTHINMETRES),
02910                      MinStrip-TPosWin,
02911                      MaxStrip+TPosWin);
02912 
02913   phHoughHist->Reset();
02914   phHoughHist->SetBins(Int_t(TPosWin/STRIPWIDTHINMETRES) + 
02915                        MaxPlane-MinPlane+2,
02916                        -2.*TPosWin/TSpecReal_t(MaxPlane-MinPlane+2),
02917                        2.*TPosWin/TSpecReal_t(MaxPlane-MinPlane+2),
02918                        Int_t(2.*TPosWin/STRIPWIDTHINMETRES),
02919                        MinStrip-TPosWin,
02920                        MaxStrip+TPosWin);
02921   
02922   winIterBeg = winVec.begin();
02923   while(winIterBeg!=winIterEnd){
02924     Int_t plane = winIterBeg->plane;
02925     Double_t strip = winIterBeg->phpos;
02926     Double_t ph = winIterBeg->ph;
02927     Int_t lastBin = -1;
02928     for(Double_t c=MinStrip-TPosWin; c<=MaxStrip+TPosWin;
02929         c+=STRIPWIDTHINMETRES){
02930       Double_t m = (strip - c)/TSpecReal_t(plane-MinPlane+2);
02931       MSG("SubShowerSR",Msg::kVerbose) << "Hough entries: " << m << " " 
02932                                        << c << " " << ph << endl;      
02933       if(houghHist->FindBin(m,c)!=lastBin){
02934         houghHist->Fill(m,c,1);
02935         phHoughHist->Fill(m,c,TMath::Sqrt(ph));
02936         lastBin = houghHist->FindBin(m,c);
02937       }
02938     }
02939     winIterBeg++;
02940   }
02941 
02942   Int_t maxBin = houghHist->GetMaximumBin();
02943   TSpecReal_t maxVal = houghHist->GetMaximum();
02944   Int_t xbin = maxBin%(houghHist->GetNbinsX()+2);
02945   Int_t ybin = maxBin/(houghHist->GetNbinsX()+2);
02946   Int_t phMaxBin = phHoughHist->GetMaximumBin();
02947   TSpecReal_t phMaxVal = phHoughHist->GetMaximum();
02948   Int_t phXbin = phMaxBin%(phHoughHist->GetNbinsX()+2);
02949   Int_t phYbin = phMaxBin/(phHoughHist->GetNbinsX()+2);
02950 
02951   Double_t totx = 0;
02952   Double_t sumx=0,sumx2=0;
02953   Double_t ptotx = 0;
02954   Double_t psumx=0,psumx2=0;
02955 
02956   Double_t toty=0;
02957   Double_t sumy=0,sumy2=0;
02958   Double_t ptoty = 0;
02959   Double_t psumy=0,psumy2=0;
02960 
02961   for(int xx = 1;xx<=houghHist->GetNbinsX();xx++){
02962     for(int yy = 1;yy<=houghHist->GetNbinsY();yy++){
02963 
02964       if(houghHist->GetBinContent(xx,yy)>maxVal*fBestHoughMaxPHFrac) {
02965         
02966         totx += houghHist->GetBinContent(xx,yy);
02967         sumx += ( houghHist->GetBinContent(xx,yy) * 
02968                   houghHist->GetXaxis()->GetBinCenter(xx)  );
02969         sumx2 += ( houghHist->GetBinContent(xx,yy) * 
02970                    houghHist->GetXaxis()->GetBinCenter(xx)  *
02971                    houghHist->GetXaxis()->GetBinCenter(xx)  );
02972         
02973         toty += houghHist->GetBinContent(xx,yy);
02974         sumy += ( houghHist->GetBinContent(xx,yy) * 
02975                   houghHist->GetYaxis()->GetBinCenter(yy)  );
02976         sumy2 += ( houghHist->GetBinContent(xx,yy) * 
02977                    houghHist->GetYaxis()->GetBinCenter(yy)  *
02978                    houghHist->GetYaxis()->GetBinCenter(yy)  );
02979       }
02980 
02981       if(phHoughHist->GetBinContent(xx,yy)>phMaxVal*fBestHoughMaxPHFrac) {
02982 
02983         ptotx += phHoughHist->GetBinContent(xx,yy);
02984         psumx += ( phHoughHist->GetBinContent(xx,yy) * 
02985                    phHoughHist->GetXaxis()->GetBinCenter(xx)  );
02986         psumx2 += ( phHoughHist->GetBinContent(xx,yy) * 
02987                     phHoughHist->GetXaxis()->GetBinCenter(xx)  *
02988                     phHoughHist->GetXaxis()->GetBinCenter(xx)  );
02989         
02990         ptoty += phHoughHist->GetBinContent(xx,yy);
02991         psumy += ( phHoughHist->GetBinContent(xx,yy) * 
02992                    phHoughHist->GetYaxis()->GetBinCenter(yy)  );
02993         psumy2 += ( phHoughHist->GetBinContent(xx,yy) * 
02994                     phHoughHist->GetYaxis()->GetBinCenter(yy)  *
02995                     phHoughHist->GetYaxis()->GetBinCenter(yy)  );
02996       }
02997 
02998     }
02999   }
03000 
03001   MCV[2] = MinPlane - 2;
03002   if(totx>0) {
03003     MCV[0] = sumx/totx;
03004     MCV[3] = (sumx2/totx - TMath::Power(MCV[0],2))/totx;
03005     if(MCV[3]>0) MCV[3] = TMath::Sqrt(MCV[3]);
03006     else MCV[3] = 0;
03007   }
03008   if(toty>0) {
03009     MCV[1] = sumy/toty;
03010     MCV[4] = (sumy2/toty - TMath::Power(MCV[1],2))/toty;
03011     if(MCV[4]>0) MCV[4] = TMath::Sqrt(MCV[4]);
03012     else MCV[4] = 0;
03013   }
03014   if(MCV[3]<houghHist->GetXaxis()->GetBinWidth(1))
03015     MCV[3] = houghHist->GetXaxis()->GetBinWidth(1);
03016   if(MCV[4]<houghHist->GetYaxis()->GetBinWidth(1))
03017     MCV[4] = houghHist->GetYaxis()->GetBinWidth(1);
03018 
03019   MCV[7] = MinPlane - 2;
03020   if(ptotx>0) {
03021     MCV[5] = psumx/ptotx;
03022     MCV[8] = (psumx2/ptotx - TMath::Power(MCV[5],2))/ptotx;
03023     if(MCV[8]>0) MCV[8] = TMath::Sqrt(MCV[8]);
03024     else MCV[8] = 0;
03025   }
03026   if(ptoty>0) {
03027     MCV[6] = psumy/ptoty;
03028     MCV[9] = (psumy2/ptoty - TMath::Power(MCV[6],2))/ptoty;
03029     if(MCV[9]>0) MCV[9] = TMath::Sqrt(MCV[9]);
03030     else MCV[9] = 0;  
03031   }
03032   if(MCV[8]<phHoughHist->GetXaxis()->GetBinWidth(1))
03033     MCV[8] = phHoughHist->GetXaxis()->GetBinWidth(1);  
03034   if(MCV[9]<phHoughHist->GetYaxis()->GetBinWidth(1))
03035     MCV[9] = phHoughHist->GetYaxis()->GetBinWidth(1);
03036 
03037   MCV[10] = houghHist->GetXaxis()->GetBinCenter(xbin);
03038   MCV[11] = houghHist->GetYaxis()->GetBinCenter(ybin);
03039   MCV[12] = phHoughHist->GetXaxis()->GetBinCenter(phXbin);
03040   MCV[13] = phHoughHist->GetYaxis()->GetBinCenter(phYbin);
03041 
03042   if(makeEPS) {
03043     gStyle->SetOptStat(0);
03044     TCanvas *can = new TCanvas();
03045     can->Divide(2,2);
03046     can->cd(1);
03047     houghHist->Draw("COLZ");
03048     can->cd(2);
03049     phHoughHist->Draw("COLZ");
03050     can->cd(3);
03051     houghHist->Draw("LEGO2");
03052     can->cd(4);
03053     phHoughHist->Draw("LEGO2");
03054     Int_t countr = 0;
03055     Bool_t noPlot = true;
03056     while(noPlot){
03057       char name[256];
03058       sprintf(name,"best_%i.root",countr);
03059       std::ifstream Test(name);
03060       if(!Test) { //if files does not exist:
03061         can->Print(name);
03062         noPlot = false;
03063       }
03064       countr+=1;
03065     }
03066     delete can;
03067   }
03068 
03069   MSG("SubShowerSR",Msg::kDebug) 
03070     << "Before Cuts in BestHough:"
03071     << "\nHough gradient, intercept, plane vertex: " 
03072     << MCV[0] << "+/-" << MCV[3] << " " 
03073     << MCV[1] << "+/-" << MCV[4] << " " << MCV[2]
03074     << "\nPH Hough gradient, intercept, plane vertex: " 
03075     << MCV[5] << "+/-" << MCV[8] << " " 
03076     << MCV[6] << "+/-" << MCV[9] << " " << MCV[7] 
03077     << "\n Hough Max Bin Gradient and Intercept Values: "
03078     << MCV[10] << " " << MCV[11] << " "
03079     << MCV[12] << " " << MCV[13] << endl;
03080 
03081   //check that there is a coincidence of 3 or more windows
03082   //and gradient is bigger than error
03083   if(maxVal < 3 || 
03084      (TMath::Abs(MCV[0])-MCV[3]<0 && 
03085       TMath::Abs(MCV[10])-MCV[3]<0)) {
03086     MCV[0] = 0;
03087     MCV[1] = maxWindowPos;
03088     MCV[2] = maxPlanePos;
03089     MCV[3] = 0;
03090     MCV[4] = 0; 
03091     MCV[10] = 0;
03092     MCV[11] = maxWindowPos;
03093   }
03094   
03095   if((phMaxVal - maxWindowPH)/(totPH - maxWindowPH) < 0.1 || 
03096      (TMath::Abs(MCV[5])-MCV[8]<0 ) &&
03097      (TMath::Abs(MCV[12])-MCV[8]<0 )) {
03098     MCV[5] = 0;
03099     MCV[6] = maxWindowPos;
03100     MCV[7] = maxPlanePos;
03101     MCV[8] = 0;
03102     MCV[9] = 0;
03103     MCV[12] = 0;
03104     MCV[13] = maxWindowPos;
03105   }
03106   return MCV;
03107   
03108 }

Bool_t AlgSubShowerSRList::CleanUp ( TObjArray *  allStrips,
CandHandle ch,
Int_t  det,
const CandSliceHandle cslice = NULL 
) [private]

Definition at line 2210 of file AlgSubShowerSRList.cxx.

References fCleanUpPlaneCut, fCleanUpStripCut, fCleanUpTimeCut, CandStripHandle::GetCharge(), CandHandle::GetDaughterIterator(), Calibrator::GetMIP(), CandStripHandle::GetPlane(), CandStripHandle::GetStrip(), CandStripHandle::GetTime(), CandStripHandle::GetTPos(), CandStripHandle::GetZPos(), Calibrator::Instance(), Msg::kDebug, CalDigitType::kSigCorr, CalStripType::kSigCorr, PlaneView::kU, PlaneView::kV, Msg::kVerbose, PlaneView::kX, PlaneView::kY, MSG, and STRIPWIDTHINMETRES.

Referenced by FormHalo(), and RunAlg().

02212 {
02213   MSG("SubShowerSR", Msg::kDebug) << "Running CleanUp()" << endl;
02214 
02215   Int_t totUnassigned = allStrips->GetLast()+1;
02216   if(totUnassigned==0) return false;
02217 
02218   Calibrator& calibrator=Calibrator::Instance();
02219 
02220   assert(ch.InheritsFrom("CandSubShowerSRListHandle"));
02221   CandSubShowerSRListHandle &csslh = 
02222     dynamic_cast<CandSubShowerSRListHandle &>(ch);
02223 
02224   //Get planeview for this cleanup
02225   TObject *firstObj = allStrips->At(0);
02226   assert(firstObj->InheritsFrom("CandStripHandle"));
02227   CandStripHandle *firstStp = dynamic_cast<CandStripHandle*>(firstObj);
02228   PlaneView::PlaneView_t pv = firstStp->GetPlaneView();
02229   MSG("SubShowerSR", Msg::kDebug) << "PlaneView = " << pv << endl;
02230 
02231   //Any changes made?
02232   Bool_t anyChanges = false;
02233   //once strips have been assigned do not try to assign again:
02234   std::vector<Bool_t> useStrip(totUnassigned,0);
02235   std::vector<Bool_t> assignedStrip(totUnassigned,0);
02236   
02237   for(int i=0;i<totUnassigned;i++) {
02238     assignedStrip[i] = false;
02239     //if cslice defined, check that strip is part of slice:
02240     if(cslice){
02241       useStrip[i] = false;
02242       TObject *tobj = allStrips->At(i);
02243       CandStripHandle *unassigned = dynamic_cast<CandStripHandle*>(tobj);
02244       CandStripHandleItr stpSliceItr(cslice->GetDaughterIterator());
02245       while(CandStripHandle *stp = stpSliceItr()) {
02246         if(*stp==*unassigned) {
02247           useStrip[i] = true;
02248           break;
02249         }
02250       }
02251     }
02252     else useStrip[i] = true;
02253   }
02254   MSG("SubShowerSR", Msg::kDebug) << "Checking " << totUnassigned 
02255                                   << " strips" << endl;
02256 
02257   //going to check if any isolated strips are nearby within +/- 2 strip
02258   for (Int_t i=0; i<totUnassigned; i++) {
02259     
02260     MSG("SubShowerSR", Msg::kVerbose) << "Checking unassigned strip " 
02261                                       << i << endl;
02262 
02263     TObject *tobj = allStrips->At(i);
02264     if(!useStrip[i]) continue;
02265     MSG("SubShowerSR", Msg::kVerbose) << "Can still use unassigned strip " 
02266                                       << i << endl;
02267 
02268     CandStripHandle *unassigned = dynamic_cast<CandStripHandle*>(tobj);
02269     std::vector<TSpecReal_t> StripDistance(csslh.GetNDaughters(),0.0);
02270     TSpecReal_t theCharge = 
02271       calibrator.GetMIP(unassigned->GetCharge(CalDigitType::kSigCorr));
02272 
02273     //loop through existing subshowers
02274     Int_t ssCount = 0;
02275     CandSubShowerSRHandleItr subshowerItr(csslh.GetDaughterIterator());
02276     while (CandSubShowerSRHandle *subshower = subshowerItr()) {       
02277 
02278       StripDistance[ssCount] = 99999; //set to 999999 as a flag
02279 
02280       //increment early because of continue statements, later use ssCount-1
02281       ssCount++; 
02282 
02283       //if a slice is given in argument list
02284       //demand that slices are the same
02285       if(cslice){
02286         if(*subshower->GetCandSlice()!=*cslice) continue;
02287       }
02288 
02289       //demand that any hit added to subshower is within +/- 50 ns
02290       if(TMath::Abs(subshower->GetAveTime() - 
02291                     unassigned->GetTime())*1.0e9 > fCleanUpTimeCut) continue;
02292       
02293       //only look at the right planeview subshowers
02294       if(subshower->GetPlaneView()!=pv) continue; 
02295       
02296       MSG("SubShowerSR", Msg::kVerbose) << "SubShower " << ssCount-1 
02297                                         << " in the same planeview" << endl;
02298 
02299       // if plane of strip is more than fCleanUpPlaneCut planes away
02300       // don't try to add it to this subshower 
02301       // (also beware of SM boundary in FarDet)
02302       Int_t thePlane = unassigned->GetPlane();
02303       if(thePlane<subshower->GetBegPlane()-fCleanUpPlaneCut || 
02304          thePlane>subshower->GetEndPlane()+fCleanUpPlaneCut || 
02305          (det==1 && 
02306           ( (subshower->GetBegPlane()>249 && thePlane<249) || 
02307             (subshower->GetEndPlane()<249 && thePlane>249) ))) continue;
02308 
02309       MSG("SubShowerSR", Msg::kVerbose) << "SubShower " << ssCount-1 
02310                                         << " in the same plane range" << endl;
02311  
02312       //loop through strips in subshower
02313       //check for proximity to other strips in subshower
02314       //(on a plane-by-plane basis)
02315       Bool_t checkProximity = false;
02316       TSpecReal_t theStrip_tpos = unassigned->GetTPos();
02317       CandStripHandleItr stripItr(subshower->GetDaughterIterator());
02318       while(CandStripHandle *stp = stripItr()){
02319         //check we're on the same plane 
02320         //or on the begin plane if the unassigned plane < begin plane
02321         //or on the end plane if the unassigned plane > end plane
02322         if((stp->GetPlane() == thePlane) ||
02323            (thePlane < subshower->GetBegPlane() && 
02324             stp->GetPlane() == subshower->GetBegPlane()) ||
02325             (thePlane > subshower->GetEndPlane() && 
02326             stp->GetPlane() == subshower->GetEndPlane())) {
02327           //allow +/-fCleanUpStripCut strips:
02328           if(TMath::Abs(theStrip_tpos-stp->GetTPos()) <=
02329              fCleanUpStripCut*STRIPWIDTHINMETRES) { 
02330             MSG("SubShowerSR",Msg::kVerbose) 
02331               << "TMath::Abs(theStrip_tpos-stp->GetTPos())=" 
02332               << TMath::Abs(theStrip_tpos-stp->GetTPos()) 
02333               << " stp->Plane()=" << stp->GetPlane()
02334               << " stp->Strip()=" << stp->GetStrip()
02335               << endl;
02336             checkProximity = true;
02337           }
02338         }
02339       }
02340       if(!checkProximity) continue;  //not close enough to this subshower
02341 
02342       MSG("SubShowerSR", Msg::kVerbose) 
02343         << "SubShower " << ssCount-1 
02344         << " has strips in the same tpos range" 
02345         << endl;
02346 
02347       //calculate pH weighted distance from cluster centre
02348       //get subshower vertex and angle:
02349       Double_t tposVtx = 0;
02350       if(pv==PlaneView::kU || 
02351          pv==PlaneView::kX) tposVtx = subshower->GetVtxU();
02352       else if(pv==PlaneView::kV || 
02353               pv==PlaneView::kY) tposVtx = subshower->GetVtxV();      
02354       Double_t zposVtx = subshower->GetVtxZ();
02355       TSpecReal_t slope = subshower->GetSlope(); 
02356       TSpecReal_t avgdev = subshower->GetAvgDev();
02357       //get strip position and charge:
02358       TSpecReal_t tpos = unassigned->GetTPos();
02359       TSpecReal_t zpos = unassigned->GetZPos();
02360       if(avgdev!=0.) { 
02361         StripDistance[ssCount-1] = 
02362           (TMath::Abs(tpos - (tposVtx + slope*(zpos - zposVtx))) * 
02363            TMath::Cos(TMath::ATan(slope)))*theCharge/avgdev; 
02364       }
02365       else {
02366         MSG("SubShowerSR",Msg::kVerbose) << "SubShower " << ssCount-1 
02367                                          << " has AvgDev = 0!" << endl;
02368         StripDistance[ssCount-1] = 99998;
02369       }
02370 
02371       MSG("SubShowerSR", Msg::kVerbose) 
02372         << "Unassigned strip " << i 
02373         << " and SubShower " << ssCount-1 
02374         << " are a good match with dev(strip)/avgdev(subshower) = " 
02375         << StripDistance[ssCount-1] << endl;
02376     }
02377 
02378     //make decision here about closest subshower
02379     TSpecReal_t minSD = 99999;
02380     Int_t bestSS = -1;
02381     for(int j=0;j<csslh.GetNDaughters();j++){
02382       if(StripDistance[j]<minSD) {
02383         minSD = StripDistance[j];
02384         bestSS = j;
02385       }
02386     }
02387     if(bestSS!=-1){
02388       MSG("SubShowerSR", Msg::kVerbose) << "Best match for strip " << i 
02389                                         << ": SubShower " 
02390                                         << bestSS << " Adding plane: " 
02391                                         << unassigned->GetPlane() 
02392                                         << " strip: " 
02393                                         << unassigned->GetStrip() << endl; 
02394       
02395       CandSubShowerSRHandleItr subshowerItr(csslh.GetDaughterIterator());
02396       ssCount = 0;
02397       while(CandSubShowerSRHandle *cssh = subshowerItr()){
02398         if(ssCount==bestSS){
02399           MSG("SubShowerSR",Msg::kVerbose) << "Before: cssh->GetNStrip() = " 
02400                                            <<  cssh->GetNStrip() 
02401                                            << " cssh->GetEnergy() = " 
02402                                            << cssh->GetEnergy() << endl; 
02403                                          
02404           cssh->AddDaughterLink(*unassigned);
02405           Double_t mip2gev = cssh->GetEnergy() / 
02406             calibrator.GetMIP(cssh->GetCharge(CalStripType::kSigCorr));
02407           cssh->SetEnergy(cssh->GetEnergy()+mip2gev*theCharge);
02408           
02409           MSG("SubShowerSR",Msg::kVerbose) << "After: cssh->GetNStrip() = " 
02410                                            << cssh->GetNStrip() 
02411                                            << " cssh->GetEnergy() = " 
02412                                            << cssh->GetEnergy() << endl;
02413           
02414           assignedStrip[i] = true;
02415           anyChanges = true;
02416           break;
02417         }
02418         ssCount+=1;
02419       }
02420     }
02421     else {
02422       MSG("SubShowerSR", Msg::kVerbose) 
02423         << "No matching subshower for unassigned strip " 
02424         << i << endl;
02425     }
02426   }
02427   
02428   //remove assigned strips
02429   Int_t assigned = 0;
02430   for(int i=0;i<totUnassigned;i++) {
02431     if(assignedStrip[i]) {
02432       allStrips->RemoveAt(i);
02433       assigned+=1;
02434     }  
02435   }
02436   allStrips->Compress();
02437   
02438   if(anyChanges) {
02439     MSG("SubShowerSR", Msg::kDebug) 
02440       << assigned << "/" << totUnassigned
02441       << " strips assigned to SubShowers" << endl;
02442   }
02443   else {
02444     MSG("SubShowerSR", Msg::kDebug) 
02445       << "No Changes made to SubShowers" << endl;
02446   }
02447   return anyChanges;
02448   
02449 }

Bool_t AlgSubShowerSRList::FindCluster ( TObjArray *  allStrips,
TObjArray *  newSubShower,
Int_t  det 
) [private]

find maximum PH plane within selected longitudinal window:

Definition at line 499 of file AlgSubShowerSRList.cxx.

References window::centroid, fLongWindowGapPlaneCut, fLongWindowMipCut, Calibrator::GetMIP(), CandHandle::GetVldContext(), Calibrator::Instance(), Msg::kDebug, Msg::kError, CalDigitType::kPE, CalDigitType::kSigCorr, Msg::kVerbose, Msg::kWarning, window::lower, window::lower_tpos, MSG, n, window::ph, window::phpos, window::phwidth, window::plane, STRIPWIDTHINMETRES, TransCluster(), window::type, window::upper, window::upper_tpos, and win.

Referenced by RunAlg().

00501 {
00502 
00503   Int_t nstp = allStrips->GetLast()+1;
00504   if(nstp==0){
00505     MSG("SubShower", Msg::kWarning) << "No strips present!" << endl;
00506     return false;
00507   }
00508 
00509   //Get useful info from candidates for clustering
00510   //Read in info:
00511   std::vector<Int_t> plane(nstp,0);
00512   std::vector<Int_t> strip(nstp,0);
00513   std::vector<Double_t> ph(nstp,0.0);
00514   std::vector<Double_t> pherr(nstp,0.0);
00515   std::vector<Double_t> z(nstp,0.0);
00516   std::vector<Double_t> tpos(nstp,0.0);
00517 
00518   //neighbouring strip info
00519   std::vector<Double_t> TransGradientPlus(nstp,0.0);
00520   std::vector<Double_t> TransGradientMinus(nstp,0.0);
00521   std::vector<Double_t> TransGradientErrorPlus(nstp,0.0);
00522   std::vector<Double_t> TransGradientErrorMinus(nstp,0.0);
00523 
00524   CandStripHandle *usefulStrip = 0;
00525   int begPlane = 500;
00526   int endPlane = 0;
00527   
00528   Int_t doubleCounter = 0;
00529   Int_t alreadyGot[500][192] = {{0}};
00530   Calibrator& calibrator=Calibrator::Instance();
00531   for (Int_t i=0; i<=allStrips->GetLast(); i++) {
00532     TObject *tobj = allStrips->At(i);
00533     assert(tobj->InheritsFrom("CandStripHandle"));
00534     CandStripHandle *stp = dynamic_cast<CandStripHandle*>(tobj);
00535     if(alreadyGot[stp->GetPlane()][stp->GetStrip()]==1) {
00536       for(int j=0;j<i-doubleCounter;j++){
00537         if(plane[j]==stp->GetPlane() && strip[j]==stp->GetStrip()){
00538           ph[j] += calibrator.GetMIP(stp->GetCharge(CalDigitType::kSigCorr));
00539           TransGradientPlus[j] = ph[j];
00540           TransGradientMinus[j] = ph[j];
00541           break;
00542         }
00543       }
00544       doubleCounter+=1;
00545       nstp -= 1;
00546     }
00547     else {
00548       alreadyGot[stp->GetPlane()][stp->GetStrip()] = 1;
00549       plane[i-doubleCounter] = stp->GetPlane();
00550       if (plane[i-doubleCounter]<begPlane) begPlane = plane[i-doubleCounter];
00551       if (plane[i-doubleCounter]>endPlane) endPlane = plane[i-doubleCounter];
00552       strip[i-doubleCounter] = stp->GetStrip();
00553       
00554       ph[i-doubleCounter] = 
00555         calibrator.GetMIP(stp->GetCharge(CalDigitType::kSigCorr));
00556       pherr[i-doubleCounter] = ph[i-doubleCounter] /
00557         TMath::Sqrt(stp->GetCharge(CalDigitType::kPE)); //pe stats err
00558       
00559       z[i-doubleCounter] = stp->GetZPos();
00560       tpos[i-doubleCounter] = stp->GetTPos();
00561       
00562       TransGradientPlus[i-doubleCounter] = ph[i-doubleCounter];
00563       TransGradientErrorPlus[i-doubleCounter] = 
00564         TMath::Power(pherr[i-doubleCounter],2);
00565       
00566       TransGradientMinus[i-doubleCounter] = ph[i-doubleCounter];
00567       TransGradientErrorMinus[i-doubleCounter] = 
00568         TMath::Power(pherr[i-doubleCounter],2);
00569       
00570       usefulStrip = stp;
00571     }
00572   }
00573 
00575   //find transverse neighbour strip index and PH gradient for each strip
00576   for(Int_t i=0;i<nstp;i++){
00577     for(Int_t j=0;j<nstp;j++){
00578       if (plane[j]==plane[i]){
00579         if(strip[j]==strip[i]-1) {
00580           TransGradientMinus[i] = TransGradientMinus[i] - ph[j];
00581           TransGradientErrorMinus[i] = TMath::Sqrt(TransGradientErrorMinus[i]+ 
00582                                                    TMath::Power(pherr[j],2));
00583         }
00584         if(strip[j]==strip[i]+1) {
00585           TransGradientPlus[i] = TransGradientPlus[i] - ph[j];
00586           TransGradientErrorPlus[i] = TMath::Sqrt(TransGradientErrorPlus[i]+ 
00587                                                   TMath::Power(pherr[j],2));
00588           TransGradientErrorMinus[i] = 0;
00589           TransGradientErrorPlus[i] = 0;
00590         }
00591       }
00592     }
00593   }
00594 
00596   //find per plane PH
00597   int npls = (endPlane-begPlane)/2+1;
00598   if(det==1&&(begPlane-249)*(endPlane-249)<0) npls = (endPlane-begPlane-1)/2+1;
00599   std::vector<Int_t> allpl(npls,0);
00600   std::vector<Double_t> allplph(npls,0.0);
00601   std::vector<Int_t> allstp(npls,0);
00602 
00603   for(Int_t i=0;i<npls;i++){
00604     allpl[i] = 2*i+begPlane;
00605     if(det==1&&(begPlane-249)*(endPlane-249)<0&&(2*i+begPlane)>249)
00606       allpl[i] = 2*i+begPlane+1;
00607     allplph[i] = 0.;
00608     allstp[i] = 0;
00609     for(Int_t j=0;j<nstp;j++){
00610       if (plane[j]==allpl[i]){
00611         allplph[i] += ph[j];
00612         allstp[i] += 1;
00613       }
00614     }
00615   }
00616   UgliGeomHandle ugh(*usefulStrip->GetVldContext());
00617 
00619   //Start Clustering: 
00621   //Find contiguous longitudinal windows
00622   MSG("SubShowerSR",Msg::kDebug) <<"start longitudinal clustering"<<endl;
00623 
00624   std::vector<int> zeroPlanes;  //vector of planes with PH = 0
00625   int start = begPlane-2;
00626   if(det==1&&(begPlane-2-249)*(begPlane-249)<0) start = begPlane-3;
00627   zeroPlanes.push_back(start-2);
00628   zeroPlanes.push_back(start);
00629   for(Int_t i=0;i<npls;i++){
00630     if (allplph[i]<=fLongWindowMipCut) zeroPlanes.push_back(allpl[i]);
00631   }
00632   int over = endPlane+2;
00633   if(det==1&&(endPlane+2-249)*(endPlane-249)<0) over = endPlane+3;
00634   zeroPlanes.push_back(over);
00635   zeroPlanes.push_back(over+2);
00636   
00637   MSG("SubShowerSR",Msg::kVerbose) << "Number of Zero Planes = " 
00638                                    << zeroPlanes.size()
00639                                    << " First Zero Plane = " << start 
00640                                    << " Last Zero Plane =  "<< over << endl;
00641   
00643   //Find "Zero PH" windows
00644 
00645   //position of first plane in a longitudinal gap window:
00646   std::vector<int> firstGapPlane;
00647   //position of last plane in a longitudinal gap window:
00648   std::vector<int> lastGapPlane;
00649   int begGap = 0;
00650   int endGap = -5;
00651   //filling gap window vectors:
00652   for (UInt_t i = 0; i<zeroPlanes.size(); i++){
00653     bool foundgap = false;
00654     if(zeroPlanes.at(i)>endGap) {
00655       begGap = zeroPlanes.at(i);
00656       endGap = zeroPlanes.at(i);
00657       int j = 0;
00658       while ((i+j+1)<zeroPlanes.size() && 
00659              zeroPlanes.at(i+j+1)-zeroPlanes.at(i+j) < 
00660              fLongWindowGapPlaneCut) {
00661         foundgap = true;
00662         endGap = zeroPlanes.at(i+j+1);
00663         j++;
00664       }
00665       if (foundgap){
00666         firstGapPlane.push_back(begGap);
00667         lastGapPlane.push_back(endGap);
00668       }
00669     }
00670   }
00671 
00672   //Find highest strip density longitudinal window
00673   Int_t maxWinBegPlane = 0;  //first boundary plane of max window
00674   Int_t maxWinEndPlane = 0;
00675   TSpecReal_t maxDensity = 0;
00676   for (UInt_t i = 0; i<firstGapPlane.size(); i++){
00677     if((i+1)<firstGapPlane.size()){
00678       TSpecReal_t density = 0;
00679       Int_t bpln = lastGapPlane.at(i);
00680       Int_t epln = firstGapPlane.at(i+1);
00681       for(int j=0;j<npls;j++){
00682         if(allpl[j]>=bpln && allpl[j]<=epln){
00683           density += TSpecReal_t(allstp[j]);
00684         }
00685       }
00686       density/=TSpecReal_t(epln-bpln+2);
00687       MSG("SubShowerSR",Msg::kDebug) << "Proto-long window [" << bpln
00688                                      << "," << epln << "] density = " 
00689                                      << density << endl;
00690       if(density > maxDensity) {
00691         maxWinBegPlane = bpln;
00692         maxWinEndPlane = epln;
00693         maxDensity = density;
00694       }
00695     }
00696   }
00697 
00698   MSG("SubShowerSR",Msg::kDebug) << "Using Longitudinal Window: "
00699                                  << maxWinBegPlane << "-"
00700                                  << maxWinEndPlane << endl;
00701   
00704   Int_t maxPlane = -1;
00705   Double_t maxPlanePH = 0;
00706   for(int i=0;i<npls;i++){
00707     if(allpl[i]<maxWinEndPlane && allpl[i]>maxWinBegPlane) {
00708       if(allplph[i]>maxPlanePH) {
00709         maxPlane = allpl[i];
00710         maxPlanePH = allplph[i];
00711       }
00712     }
00713   }
00714   MSG("SubShowerSR",Msg::kVerbose) 
00715     << "Maximum PH plane within Long. window = " 
00716     << maxPlane << endl;
00717 
00718   Int_t numPreMaxPlanes = maxPlane - maxWinBegPlane;
00719   Int_t numPostMaxPlanes = maxWinEndPlane - maxPlane;
00720   Int_t numLongPlanes = numPostMaxPlanes + numPreMaxPlanes - 1;
00721 
00723   //Start Transverse Clustering
00724   
00725   //window info:
00726   std::vector<window> winvec;  //vector containing all transverse windows
00727 
00729   //Loop over all planes in long. window
00730   //Find transverse windows per plane
00731 
00732   for (int i=0;i<numLongPlanes;i++){  //loop over numLongPlanes
00733  
00734     //find next appropriate plane to look at starting from max plane:
00735     int pl = -1;
00736     if(numPreMaxPlanes>numPostMaxPlanes){
00737       if (i<numPreMaxPlanes) pl = maxPlane-i;
00738       else pl = maxPlane + i - numPreMaxPlanes + 1; 
00739       //add 1 to prevent re-doing max plane
00740     }
00741     else{
00742       if (i<numPostMaxPlanes) pl = maxPlane+i;
00743       else pl = maxPlane - i + numPostMaxPlanes - 1;
00744       //subtract 1 to prevent re-doing max plane
00745     }
00746 
00747     std::vector<Int_t> winu;  // vector of upper dip
00748     std::vector<Int_t> wind;  // vector of lower dip
00749     std::vector<Double_t> winu_tpos;  // vector of tpos of upper dip
00750     std::vector<Double_t> wind_tpos;  // vector of tpos of lower dip
00751     std::vector<Int_t> wintu; // vector of upper dip type
00752     std::vector<Int_t> wintd; // vector of lower dip type
00753 
00754     //find dips and categorize them
00755     MSG("SubShowerSR",Msg::kVerbose) << "Find Upper/Lower boundaries in plane = "
00756                                      << pl << endl;
00757     for(Int_t i=0;i<nstp;i++){
00758       if(plane[i]==pl){
00759         //using neighbouring strip info from above:
00760         //(only called dips if delta PH is significant given pe stats)
00761         if(TransGradientPlus[i]<0 /*-(TransGradientErrorPlus[i])*/ &&
00762            TransGradientMinus[i]<0 /*-(TransGradientErrorMinus[i])*/){
00763           winu.push_back(strip[i]);
00764           winu_tpos.push_back(tpos[i]);
00765           wintu.push_back(0);      // type-0 (this is a dip) 
00766           MSG("SubShowerSR",Msg::kVerbose) 
00767             << "up dip strip = " << strip[i] << " ph = " << ph[i] 
00768             << " tgp = " << TransGradientPlus[i]
00769             << " tgm = " << TransGradientMinus[i] << endl;
00770         }
00771         else if(TransGradientPlus[i]==ph[i] && 
00772                 TransGradientMinus[i]<=ph[i]){
00773           winu.push_back(strip[i]);
00774           winu_tpos.push_back(tpos[i]);
00775           MSG("SubShowerSR",Msg::kVerbose) 
00776             << "up edge/iso strip = " << strip[i] << " ph = " << ph[i] 
00777             << " tgp = " << TransGradientPlus[i]
00778             << " tgm = " << TransGradientMinus[i] << endl;
00779           if (TransGradientMinus[i]<ph[i]) 
00780             wintu.push_back(1);     // type-1 (this is an edge strip) 
00781           else wintu.push_back(2);  // type-2 (this is an isolated strip)
00782         }
00783         if(TransGradientPlus[i]<0 /*-(TransGradientErrorPlus[i])*/ &&
00784            TransGradientMinus[i]<0 /*-(TransGradientErrorMinus[i])*/ ){
00785           wind.push_back(strip[i]);
00786           wind_tpos.push_back(tpos[i]);
00787           wintd.push_back(0);      // type-0
00788           MSG("SubShowerSR",Msg::kVerbose) 
00789             << "down dip strip = " << strip[i] << " ph = " << ph[i]
00790             << " tgp = " << TransGradientPlus[i]
00791             << " tgm = " << TransGradientMinus[i] << endl;
00792         }
00793         else if(TransGradientMinus[i]==ph[i] && 
00794                 TransGradientPlus[i]<=ph[i]){
00795           wind.push_back(strip[i]);
00796           wind_tpos.push_back(tpos[i]);
00797           MSG("SubShowerSR",Msg::kVerbose) 
00798             << "down edge/iso strip = "  << strip[i] << " ph = " << ph[i]
00799             << " tgp = "  << TransGradientPlus[i]
00800             << " tgm = "  << TransGradientMinus[i] << endl;
00801           if (TransGradientPlus[i]<ph[i]) 
00802             wintd.push_back(1);     //type-1
00803           else wintd.push_back(2);  //type-2
00804         }
00805       }
00806     }
00807 
00809     //use dips to form windows                                   
00810     UInt_t wins = winu.size();
00811     if(wins!=wind.size()) {
00812       MSG("SubShowerSR",Msg::kError) 
00813         << "Number of Upper Transverse Window Boundaries does not equal " 
00814         << "number of Lower Transverse Window Boundaries... "
00815         << "Leaving FindCluster()" << endl;
00816       MSG("SubShowerSR",Msg::kError) << "winu.size() = " << winu.size() << endl;
00817       MSG("SubShowerSR",Msg::kError) << "wind.size() = " << wind.size() << endl;
00818       MSG("SubShowerSR",Msg::kError) << "wins = " << wins << endl;
00819       
00820       return false;
00821     }
00822 
00823     //temporarily hold window info
00824     std::vector<Int_t> win0(wins,0);             // strip of upper dip
00825     std::vector<Int_t> win1(wins,0);             // strip of lower dip
00826     std::vector<TSpecReal_t> win0_tpos(wins,0.0);      // tpos of strip of upper dip
00827     std::vector<TSpecReal_t> win1_tpos(wins,0.0);      // tpos of strip of lower dip
00828     std::vector<TSpecReal_t> win2(wins,0.0);           // PH of window
00829     std::vector<Int_t> win3(wins,0);             // ID of window:
00830                                                  //   0-both dips; 
00831                                                  //   1 or 10-one dip,one gap; 
00832                                                  //   11-both gaps; 
00833                                                  //   22-isolated strip
00834     std::vector<TSpecReal_t> win4(wins,0.0);           // Centroid of window in tpos
00835     std::vector<TSpecReal_t> win5(wins,0.0);           // PH weighted tpos of window
00836     std::vector<TSpecReal_t> win6(wins,0.0);           // PH weighted width of window
00837 
00838     //match up dips to form windows
00839     for (UInt_t w = 0; w<wins; w++){
00840       win0[w] = -1;
00841       win1[w] = -1;
00842       win0_tpos[w] = 0;
00843       win1_tpos[w] = 0;
00844       win2[w] = 0;
00845       win3[w] = -1;
00846       win4[w] = 0;
00847       win5[w] = 0;
00848       win6[w] = 0;
00849       int winw = 200;
00850       for (UInt_t s = 0; s<wind.size(); s++){
00851         int diff = winu.at(w)-wind.at(s);
00852         if ( ( (diff>0  && wintu.at(w)<2  && wintd.at(s)<2) ||
00853                (diff==0 && wintu.at(w)==2 && wintd.at(s)==2) ) &&
00854              diff<winw) {
00855           winw = diff;
00856           win0[w] = winu.at(w);
00857           win1[w] = wind.at(s);
00858           win0_tpos[w] = winu_tpos.at(w);
00859           win1_tpos[w] = wind_tpos.at(s);
00860           win3[w] = wintu.at(w)*10+wintd.at(s);
00861         }
00862       }
00863     }
00864 
00866     //now get PH of windows
00867     //and push_back windows into vectors
00868     for (UInt_t w = 0; w<wins; w++){            
00869       Double_t biggestStripPH = 0;
00870       Int_t n = 0;
00871       for(Int_t i=0;i<nstp;i++){
00872         if (plane[i]==pl && strip[i]<=win0[w] && strip[i]>=win1[w]) {
00873           win2[w] += ph[i];
00874           if(ph[i]>biggestStripPH) {
00875             win4[w] = tpos[i];
00876             biggestStripPH = ph[i];
00877           }
00878           win5[w] += ph[i]*tpos[i];
00879           win6[w] += ph[i]*tpos[i]*tpos[i];
00880           n++;
00881         }
00882       }
00883       win5[w] /= win2[w];
00884       if(n!=1){
00885         win6[w] /= win2[w];
00886         win6[w] -= win5[w]*win5[w];
00887         if(win6[w]>0) win6[w] = TMath::Sqrt(win6[w]);
00888         else win6[w] = STRIPWIDTHINMETRES/TMath::Sqrt(12.); 
00889      }
00890       else win6[w] = STRIPWIDTHINMETRES/TMath::Sqrt(12.);
00891       
00892       window win;
00893       win.plane = pl;
00894       win.upper = win0[w];
00895       win.lower = win1[w];
00896       win.upper_tpos = win0_tpos[w];
00897       win.lower_tpos = win1_tpos[w];
00898       win.ph = win2[w];
00899       win.type = win3[w];
00900       win.centroid = win4[w];
00901       win.phpos = win5[w];
00902       win.phwidth = win6[w];
00903       winvec.push_back(win);
00904 
00905     }
00906   }
00907 
00909   //Perform clustering on windows:
00910 
00911   // if this longitudinal window only contains a single 
00912   // active plane then just return largest window
00913   if(numLongPlanes/2<2){ 
00914     std::vector<window>::iterator winIter = winvec.begin();
00915     std::vector<window>::iterator winEnd = winvec.end();
00916     Int_t best_upper_strip = 0;
00917     Int_t best_lower_strip = 200;
00918     Double_t best_pulse_height = 0;
00919     while(winIter!=winEnd){
00920       if(winIter->ph>best_pulse_height){
00921         best_lower_strip = winIter->lower;
00922         best_upper_strip = winIter->upper;
00923         best_pulse_height = winIter->ph;
00924       }
00925       winIter++;
00926     }
00927     for (Int_t i=0; i<=allStrips->GetLast(); i++) {
00928       TObject *tobj = allStrips->At(i);
00929       assert(tobj->InheritsFrom("CandStripHandle"));
00930       CandStripHandle *stp = dynamic_cast<CandStripHandle*>(tobj);
00931       if(stp->GetPlane()>=maxWinBegPlane && stp->GetPlane()<=maxWinEndPlane &&
00932          stp->GetStrip()>=best_lower_strip && stp->GetStrip()<=best_upper_strip){
00933         newSubShower->Add(stp);
00934       }
00935     }
00936     return true;
00937   }
00938 
00939   //call transverse clustering algorithms
00940   std::vector<window> clusterWinVec = TransCluster(winvec,det);
00941   
00942   //no windows added to cluster:
00943   if(clusterWinVec.size()==0) return false;
00944 
00945   //Clean up proto-subshower
00946   //Check for longitudinal gaps within the proto-subshower
00947   //Remove windows after first gap on both sides of shower max plane
00948   maxPlane = -1;  //recalculate maxPlane
00949   maxPlanePH = 0.;
00950   Double_t PlanePH[500] = {0.};
00951   Int_t non0Win = 0;
00952   std::vector<window>::iterator winIter = clusterWinVec.begin();
00953   std::vector<window>::iterator winEnd = clusterWinVec.end();
00954   while(winIter!=winEnd){
00955     if(winIter->plane>0){
00956       PlanePH[winIter->plane] += winIter->ph;
00957       if(PlanePH[winIter->plane]>0) non0Win++;
00958       if(PlanePH[winIter->plane]>maxPlanePH) {
00959         maxPlane = winIter->plane;
00960         maxPlanePH = PlanePH[maxPlane];
00961       }
00962     }
00963     winIter++;
00964   }
00965   
00966   if(non0Win==0){
00967     MSG("SubShowerSR",Msg::kDebug) <<"No window solutions found"<<endl;
00968     return false;
00969   }
00970 
00971   Int_t maxLongPlane = -1;
00972   Int_t minLongPlane = -1;
00973   Int_t counter = 0;
00974   Int_t delta = 2;
00975   while((maxPlane+counter)>=0&&(PlanePH[maxPlane+counter]>0||(maxPlane+counter)==249)) {
00976     maxLongPlane = maxPlane+counter;
00977     delta = 2;
00978     if(det==1&&(maxLongPlane-249)*(maxLongPlane+delta-249)<0)
00979       delta = 3;
00980     counter+=delta;
00981   }
00982   counter = 0;
00983   while((maxPlane+counter)>=0&&(PlanePH[maxPlane+counter]>0||(maxPlane+counter)==249)) {
00984     minLongPlane = maxPlane+counter;
00985     delta = 2;
00986     if(det==1&&(minLongPlane-249)*(minLongPlane-delta-249)<0)
00987       delta = 3;
00988     counter-=delta;
00989   }
00990   MSG("SubShowerSR",Msg::kDebug) <<"New Longitudinal Window: "<< minLongPlane
00991                                  << "-" << maxLongPlane << endl;
00992 
00993   for (Int_t i=0; i<=allStrips->GetLast(); i++) {
00994     TObject *tobj = allStrips->At(i);
00995     assert(tobj->InheritsFrom("CandStripHandle"));
00996     CandStripHandle *stp = dynamic_cast<CandStripHandle*>(tobj);
00997     if(stp->GetPlane()>=minLongPlane && stp->GetPlane()<=maxLongPlane){
00998       winIter = clusterWinVec.begin();
00999       while(winIter!=winEnd){
01000         if(winIter->plane==stp->GetPlane()){
01001           Int_t strp = stp->GetStrip();
01002           if(strp<=winIter->upper && strp>=winIter->lower) {
01003             newSubShower->Add(stp);
01004             break;
01005           }
01006         }
01007         winIter++;
01008       }
01009     }
01010   }
01011 
01012   return true;
01013 }

Bool_t AlgSubShowerSRList::FormHalo ( TObjArray *  allStrips,
CandHandle ch,
CandContext cxx,
AlgHandle ah,
const CandSliceListHandle slicelist,
Int_t  det 
) [private]

Definition at line 2452 of file AlgSubShowerSRList.cxx.

References CandHandle::AddDaughterLink(), CleanUp(), fHaloMaxPlaneGap, fHaloNPeakFindingBins, fMinStripPE, fSplitShowers, CandHandle::GetDaughterIterator(), CandHandle::GetNDaughters(), CandRecoHandle::GetNStrip(), Msg::kDebug, Msg::kError, ClusterType::kHalo, CalDigitType::kSigCorr, PlaneView::kU, PlaneView::kV, PlaneView::kX, ClusterType::kXTalk, PlaneView::kY, CandSubShowerSR::MakeCandidate(), MSG, CandRecoHandle::SetCandSlice(), CandSubShowerSRHandle::SetClusterID(), CandContext::SetDataIn(), CandSubShowerSRHandle::SetMinStripPE(), SplitSubShower(), TimeTest(), and vtxHist.

Referenced by RunAlg().

02456 {
02457 
02458   MSG("SubShowerSR", Msg::kDebug) << "Running FormHalo()" << endl;
02459   
02460   assert(ch.InheritsFrom("CandSubShowerSRListHandle"));
02461 
02462   Int_t totUnassigned = allStrips->GetLast()+1;
02463   if(totUnassigned==0) return false;
02464   
02465   MSG("SubShowerSR", Msg::kDebug) << "Forming slice-wise, " 
02466                                   << "plane-wise Halo clusters from " 
02467                                   << totUnassigned 
02468                                   << " unassigned strips" << endl;
02469   
02470   //loop over slices
02471   Int_t nslice = 0;
02472   CandSliceHandleItr sliceItr(slicelist->GetDaughterIterator());
02473   while (CandSliceHandle *slice = sliceItr()) {
02474     MSG("SubShowerSR",Msg::kDebug) << "Slice " << nslice << endl;
02475     
02476     //get valid plane ranges for this slice:
02477     CandSubShowerSRListHandle &subshowerlist = 
02478       dynamic_cast<CandSubShowerSRListHandle&>(ch);
02479     CandSubShowerSRHandleItr subshowerItr(subshowerlist.GetDaughterIterator());
02480     Int_t plane[500] = {0};
02481     Int_t firstPlane = 500;
02482     Int_t lastPlane = 0;
02483     while (CandSubShowerSRHandle *subshower = subshowerItr()) {
02484       if (*subshower->GetCandSlice()==*slice) { 
02485         for(int i=subshower->GetBegPlane();i<=subshower->GetEndPlane();i++){
02486           plane[i] = 1;
02487           if(i<firstPlane) firstPlane = i;
02488           if(i>lastPlane) lastPlane = i;
02489         }
02490       }
02491     }
02492 
02493     //loop over valid plane ranges and look for gaps > fHaloMaxPlaneGap
02494     while(firstPlane<=lastPlane){
02495       Int_t begPlane = 500;
02496       Int_t endPlane = 0;
02497       Int_t nGaps = 0;
02498       for(int i=firstPlane;i<=lastPlane;i++){
02499         if(plane[i]==0) {
02500           if(begPlane!=500) nGaps+=1;
02501         }
02502         else {
02503           if(begPlane==500) begPlane = i;
02504           endPlane = i;
02505         }
02506         plane[i] = 0;
02507         if(nGaps>fHaloMaxPlaneGap) {
02508           break;
02509         }
02510       }
02511       firstPlane = endPlane+1;
02512 
02513       //use plane range between begPlane and endPlane to form subshowers
02514       //loop over views:
02515       for(Int_t pv=0;pv<2;pv++){
02516 
02517         MSG("SubShowerSR",Msg::kDebug) << "View : " << pv << endl;
02518 
02519         //first loop through subshower list again and histogram tposvtx
02520         //in order to look for peaks
02521         vtxHist->Reset();
02522         vtxHist->SetBins(fHaloNPeakFindingBins+2,
02523                  -4. - 8./TSpecReal_t(fHaloNPeakFindingBins),
02524                  +4. + 8./TSpecReal_t(fHaloNPeakFindingBins));
02525 
02526         subshowerItr.ResetFirst();
02527         while (CandSubShowerSRHandle *subshower = subshowerItr()) {
02528           if (*subshower->GetCandSlice()==*slice && 
02529               subshower->GetClusterID()!=ClusterType::kXTalk) {
02530             if(pv==0 && 
02531                (subshower->GetPlaneView()==PlaneView::kV || 
02532                 subshower->GetPlaneView()==PlaneView::kY)) continue;
02533             else if(pv==1 && 
02534                     (subshower->GetPlaneView()==PlaneView::kU || 
02535                      subshower->GetPlaneView()==PlaneView::kX)) continue;
02536             if(subshower->GetBegPlane()>=begPlane && 
02537                subshower->GetBegPlane()<=endPlane){
02538               if(pv==0) {
02539                 MSG("SubShowerSR",Msg::kDebug) 
02540                   << "Filling vtxHist with Subshower "
02541                   << " with tposvtv = " << subshower->GetVtxU() 
02542                   << " and energy = " << subshower->GetEnergy() 
02543                   << " and slope = " << subshower->GetSlope()
02544                   << " and AvgDev = " << subshower->GetAvgDev()
02545                   << endl;
02546                 CandStripHandleItr stripssItr(subshower->GetDaughterIterator());
02547                 while (CandStripHandle *stp = stripssItr()) {
02548                   //fill with max charge for this tpos
02549                   TSpecReal_t val = stp->GetCharge(CalDigitType::kSigCorr) - 
02550                     vtxHist->GetBinContent(vtxHist->FindBin(stp->GetTPos()));
02551                   if(val>0) vtxHist->Fill(stp->GetTPos(),val);
02552                 }
02553               }
02554               else if(pv==1) {
02555                 MSG("SubShowerSR",Msg::kDebug) 
02556                   << "Filling vtxHist with Subshower "
02557                   << " with tposvtv = " << subshower->GetVtxV() 
02558                   << " and energy = " << subshower->GetEnergy() 
02559                   << " and slope = " << subshower->GetSlope()
02560                   << " and AvgDev = " << subshower->GetAvgDev()
02561                   << endl;
02562                 CandStripHandleItr stripssItr(subshower->GetDaughterIterator());
02563                 while (CandStripHandle *stp = stripssItr()) {
02564                   //fill with max charge for this tpos
02565                   TSpecReal_t val = stp->GetCharge(CalDigitType::kSigCorr) - 
02566                     vtxHist->GetBinContent(vtxHist->FindBin(stp->GetTPos()));
02567                   if(val>0) vtxHist->Fill(stp->GetTPos(),val);            
02568                 }
02569               }
02570             }
02571           }
02572         }
02573 
02574         if(false){
02575           TCanvas *can = new TCanvas();
02576           can->cd(1);
02577           vtxHist->Draw();
02578           char nomus[256];
02579           sprintf(nomus,"peakPlot_FH_view%i_slice%i.eps",pv,nslice);
02580           can->Print(nomus);
02581           delete can;
02582         }
02583         
02584         //width of shower in units of number of bins
02585         //assume shower typically has a width of ~4 strips
02586         Double_t sigma = (4.*0.0417/8.)*(fHaloNPeakFindingBins);
02587         //don't let width be <=1 bin to prevent
02588         //against clipping error messages
02589         if(sigma<=1.0) sigma = 1.01;
02590 
02591         Int_t nPeaks = 0;
02592         TSpecReal_t peakPos[10] = {};
02593         TSpecReal_t peakVal[10] = {};
02594         TSpecReal_t tmp_peakPos = vtxHist->GetBinCenter(vtxHist->GetMaximumBin());
02595         TSpecReal_t tmp_peakVal = vtxHist->GetMaximum();
02596         for(int bins = 1;bins<=vtxHist->GetNbinsX();bins++){
02597           if(vtxHist->GetBinContent(bins)>0) nPeaks++;
02598         }
02599         if(nPeaks==1) {
02600           peakPos[0] = tmp_peakPos;
02601           peakVal[0] = tmp_peakVal;
02602         }
02603         else if(nPeaks>0) {
02604           TSpectrum spec(10);
02605           TSpecReal_t *source = new TSpecReal_t[vtxHist->GetNbinsX()];
02606           TSpecReal_t *dest   = new TSpecReal_t[vtxHist->GetNbinsX()];
02607           for(int h_loop=0;h_loop<vtxHist->GetNbinsX();h_loop++){
02608             source[h_loop] = vtxHist->GetBinContent(h_loop+1);
02609             dest[h_loop] = 0;
02610           }
02611 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,9,1)
02612           nPeaks = spec.SearchHighRes(source,dest,
02613                                       vtxHist->GetNbinsX(),
02614                                       sigma,5,0,1,0,1);
02615 #else
02616           nPeaks = spec.Search1HighRes(source,dest,
02617                                       vtxHist->GetNbinsX(),
02618                                       sigma,5,0,1,0,1);
02619 #endif
02620           TSpecReal_t *posX = spec.GetPositionX();
02621           for(Int_t p_loop=0;p_loop<nPeaks;p_loop++){
02622             Int_t bin = 1+Int_t(posX[p_loop]+0.5);
02623             peakPos[p_loop] = vtxHist->GetBinCenter(bin);
02624             peakVal[p_loop] = vtxHist->GetBinContent(bin);
02625           }
02626           delete [] source;
02627           delete [] dest;
02628           //spec.Search(vtxHist,sigma);
02629           //nPeaks = spec.GetNPeaks();
02630         }
02631 
02632         //if no peaks are found
02633         //then there's no reconstructed subshower 
02634         //in this slice + plane range + view
02635         if(nPeaks==0) {
02636           MSG("SubShowerSR",Msg::kDebug) << "No peaks found" << endl;
02637           //so take everything and put it in a subshower
02638           CandStripHandleItr stripslcItr(slice->GetDaughterIterator());
02639           TObjArray newSubShower;
02640           while (CandStripHandle *stp = stripslcItr()) {
02641             if(pv==0 && (stp->GetPlaneView()==PlaneView::kV || 
02642                          stp->GetPlaneView()==PlaneView::kY)) continue;
02643             else if(pv==1 && (stp->GetPlaneView()==PlaneView::kU || 
02644                               stp->GetPlaneView()==PlaneView::kX)) continue;
02645             //is this strip in the subshower plane range?
02646             //allow +/- 2 planes in this window
02647             if(stp->GetPlane()>endPlane+2 || 
02648                stp->GetPlane()<begPlane-2) continue;
02649             //loop over unassigned strips, and see if it's there
02650             for (Int_t i=0; i<allStrips->GetLast()+1; i++) {
02651               TObject *tobj = allStrips->At(i);
02652               CandStripHandle *unassigned = 
02653                 dynamic_cast<CandStripHandle*>(tobj);
02654               if(*unassigned == *stp) {
02655                 newSubShower.Add(unassigned);
02656                 allStrips->RemoveAt(i);
02657               }
02658             }
02659             allStrips->Compress();
02660           }
02661           //make a new subshower with the unassigned strips in this slice
02662           if(newSubShower.GetLast()+1>0){
02663             TimeTest(&newSubShower,allStrips);
02664             std::vector<TObjArray*> showers;
02665             //if splitShowers is turned on, create a vector of split subshowers
02666             //otherwise, just put the full subshower in the vector
02667             if (fSplitShowers) showers = this->SplitSubShower(&newSubShower);
02668             else showers.push_back(&newSubShower);
02669             for (UInt_t i = 0; i<showers.size(); i++){
02670               TObjArray* splitSubShower = showers.at(i);
02671               cxx.SetDataIn(splitSubShower);
02672               MSG("SubShowerSR",Msg::kDebug) << "forming SubShowerSR\n";
02673               CandSubShowerSRHandle subshowersrhandle = 
02674                 CandSubShowerSR::MakeCandidate(ah,cxx);
02675               // Set the strip PE cut - this gets passed through the subshower
02676               // to the shower, so that it is accessible from the truthhelper
02677               // when it fills the purity and completeness variables
02678               subshowersrhandle.SetMinStripPE(fMinStripPE);
02679               
02680               if(subshowersrhandle.GetNStrip()!=0) {
02681                 MSG("SubShowerSR",Msg::kDebug) << "Formed SubShowerSR " 
02682                                                << ch.GetNDaughters()
02683                                                << " with " 
02684                                                << subshowersrhandle.GetNStrip()
02685                                                << " strips" << endl;
02686                 subshowersrhandle.SetClusterID(ClusterType::kHalo);
02687                 subshowersrhandle.SetCandSlice(slice);
02688                 ch.AddDaughterLink(subshowersrhandle);        
02689               }
02690             }
02691           }
02692           continue;
02693         }
02694       
02695         //have one or more peak => maybe one or 
02696         //more shower in this plane range
02697         //first find acceptable limits in tpos for a halo
02698         //subshower to be formed:
02699         std::vector<TSpecReal_t> peakPosLow(nPeaks,0.0);
02700         std::vector<TSpecReal_t> peakPosUpp(nPeaks,0.0);
02701         MSG("SubShowerSR",Msg::kDebug) << "Number of Peaks found = "
02702                                        << nPeaks 
02703                                        << " with positions, heights "
02704                                        << "and bounds: "  << endl;
02705         for (int i=0;i<nPeaks;i++){
02706           if(i==0) peakPosLow[i] = -4;
02707           else peakPosLow[i] = peakPosUpp[i-1];
02708           if(i==nPeaks-1) peakPosUpp[i] = 4;
02709           else {
02710             if(peakVal[i]+peakVal[i+1] == 0) {
02711               MSG("SubShowerSR",Msg::kError) << "peakVal[" << i << "]=" << peakVal[i] 
02712                                              << " peakVal[" << i+1 << "]=" << peakVal[i+1]
02713                                              << " - Setting nPeaks=1!" << endl;
02714               nPeaks=1;
02715               peakPosUpp[0] = 4;
02716               break;
02717             }
02718             peakPosUpp[i] = peakPos[i] + 
02719               (peakVal[i])*(peakPos[i+1] - peakPos[i]) / 
02720               ((peakVal[i]) + (peakVal[i+1]));
02721           }
02722           MSG("SubShowerSR",Msg::kDebug) << peakPos[i] << " " 
02723                                          << peakVal[i] << " ["
02724                                          << peakPosLow[i] << "," 
02725                                          << peakPosUpp[i] << "]" 
02726                                          << endl;
02727         }
02728 
02729         //loop through peaks and assign strips to form halo subshowers
02730         //according to acceptable tpos and plane ranges
02731         for(int i=0;i<nPeaks;i++){
02732           //loop through slice strips
02733           CandStripHandleItr stripslcItr(slice->GetDaughterIterator());
02734           TObjArray newSubShower;
02735         
02736           //loop over strips in slices in this view:
02737           // (have to loop through slice strips and not just
02738           // unassigned strips because you cannot get slice number
02739           // from strip handle)
02740           while (CandStripHandle *stp = stripslcItr()) {
02741           
02742             //only look at one view at a time:
02743             if(pv==0 && (stp->GetPlaneView()==PlaneView::kV || 
02744                          stp->GetPlaneView()==PlaneView::kY)) continue;
02745             else if(pv==1 && (stp->GetPlaneView()==PlaneView::kU || 
02746                               stp->GetPlaneView()==PlaneView::kX)) continue;
02747           
02748             //is this strip in the subshower plane range?
02749             //allow +/- 2 planes in this window
02750             if(stp->GetPlane()>endPlane+2 || 
02751                stp->GetPlane()<begPlane-2) continue;
02752             
02753             //is this strip in the tpos range
02754             if(stp->GetTPos()>peakPosUpp[i] ||
02755                stp->GetTPos()<=peakPosLow[i]) continue;
02756 
02757             //loop over unassigned strips, and see if it's there
02758             for (Int_t i=0; i<allStrips->GetLast()+1; i++) {
02759               TObject *tobj = allStrips->At(i);
02760               CandStripHandle *unassigned = 
02761                 dynamic_cast<CandStripHandle*>(tobj);
02762               if(*unassigned == *stp) {
02763                 newSubShower.Add(unassigned);
02764                 allStrips->RemoveAt(i);
02765               }
02766             }
02767             allStrips->Compress();
02768           }
02769           
02770           //make a new subshower with the unassigned strips in this slice
02771           if(newSubShower.GetLast()+1>0){  
02772             TimeTest(&newSubShower,allStrips); 
02773             std::vector<TObjArray*> showers;
02774             //if splitShowers is turned on, create a vector of split subshowers
02775             //otherwise, just put the full subshower in the vector
02776             if (fSplitShowers) showers = this->SplitSubShower(&newSubShower);
02777             else showers.push_back(&newSubShower);
02778             for (UInt_t i = 0; i<showers.size(); i++){
02779               TObjArray* splitSubShower = showers.at(i);
02780               
02781               cxx.SetDataIn(splitSubShower);
02782               MSG("SubShowerSR",Msg::kDebug) << "forming SubShowerSR\n";
02783               CandSubShowerSRHandle subshowersrhandle = 
02784                 CandSubShowerSR::MakeCandidate(ah,cxx);
02785               // Set the strip PE cut - this gets passed through the subshower
02786               // to the shower, so that it is accessible from the truthhelper
02787               // when it fills the purity and completeness variables
02788               subshowersrhandle.SetMinStripPE(fMinStripPE);
02789               
02790               if(subshowersrhandle.GetNStrip()!=0) {
02791                 MSG("SubShowerSR",Msg::kDebug) << "Formed SubShowerSR with " 
02792                                                << subshowersrhandle.GetNStrip()
02793                                                << " strips" << endl;
02794                 subshowersrhandle.SetClusterID(ClusterType::kHalo);
02795                 subshowersrhandle.SetCandSlice(slice);
02796                 ch.AddDaughterLink(subshowersrhandle); 
02797               }
02798             }
02799           }
02800         }
02801       }
02802     }
02803     //slice-wise clean up
02804     Int_t numberOfCleanUps = 0;
02805     while(CleanUp(allStrips,ch,det,slice)) {
02806       numberOfCleanUps+=1;
02807       //prevent against infinite loop! (which should NEVER happen)
02808       if(numberOfCleanUps>1000) break;
02809     }
02810     nslice++;
02811   }
02812   //snarl-wide clean up
02813   Int_t numberOfCleanUps = 0;
02814   while(CleanUp(allStrips,ch,det)) {
02815     numberOfCleanUps+=1;
02816     //prevent against infinite loop! (which should NEVER happen)
02817     if(numberOfCleanUps>1000) break;
02818   }
02819   return true;
02820 }

std::vector< window > AlgSubShowerSRList::HoughTransCluster ( std::vector< window win,
Int_t  det 
) [private]

Definition at line 1686 of file AlgSubShowerSRList.cxx.

References BestHough(), fMinHoughPH, fMinHoughPlanes, Msg::kDebug, MAXHOUGHITER, MSG, size, and STRIPWIDTHINMETRES.

Referenced by TransCluster().

01688 {
01689 
01690   MSG("SubShowerSR",Msg::kDebug) << "In HoughTransCluster()" << endl;
01691   
01692   std::vector<window> cluster;
01693   std::vector<window> nonCluster;
01694   std::vector<window> emptyCluster;
01695   
01696   //copy all windows into nonCluster:
01697   //also get best vertex estimate in first plane
01698   Int_t vertexPlane = 500;
01699   Double_t vertexPlanePH = 0;
01700   Double_t vertexTPos = 0;
01701   Double_t vertexTPosWidth = 0;
01702 
01703   //use these to calculate the density of hits in the long. window
01704   //this can be used to decide which Hough method to use
01705   Double_t totStrips = 0; 
01706   Double_t maxTpos = -5.0;
01707   Double_t minTpos = 5.0;
01708   Int_t lastPlane = 0;
01709 
01710   std::vector<window>::iterator wIt = win.begin();
01711   std::vector<window>::iterator wEn = win.end();
01712   while(wIt!=wEn){
01713     nonCluster.push_back(*wIt);
01714     if(wIt->plane<vertexPlane) {
01715       if(wIt->ph>1.0) {
01716         vertexPlane = wIt->plane;
01717         vertexPlanePH = wIt->ph;
01718         vertexTPos = wIt->phpos;
01719         vertexTPosWidth = wIt->phwidth;
01720         if(vertexTPosWidth<STRIPWIDTHINMETRES) vertexTPosWidth=STRIPWIDTHINMETRES;
01721       }
01722       else if(vertexPlanePH<1.0){
01723         vertexPlane = wIt->plane;
01724         vertexPlanePH = wIt->ph;
01725         vertexTPos = wIt->phpos;
01726         vertexTPosWidth = wIt->phwidth;
01727         if(vertexTPosWidth<STRIPWIDTHINMETRES) vertexTPosWidth=STRIPWIDTHINMETRES;
01728       }
01729     }
01730     else if(wIt->plane==vertexPlane){
01731       if(wIt->ph>vertexPlanePH){
01732         vertexPlanePH = wIt->ph;
01733         vertexTPos = wIt->phpos;
01734         vertexTPosWidth = wIt->phwidth;
01735         if(vertexTPosWidth<STRIPWIDTHINMETRES) vertexTPosWidth=STRIPWIDTHINMETRES;
01736       }
01737     }
01738     if(wIt->ph>1.0) {
01739       if(wIt->plane>lastPlane) lastPlane = wIt->plane;
01740       if(wIt->upper_tpos>maxTpos) maxTpos = wIt->upper_tpos;
01741       if(wIt->lower_tpos<minTpos) minTpos = wIt->lower_tpos;      
01742       totStrips += Double_t(wIt->upper - wIt->lower + 1);
01743     }
01744     wIt++;
01745   }
01746 
01747   //keep track of some quantities for quality checking at the end:
01748   Int_t   lowPlane  = 500;
01749   Int_t   uppPlane  = 0;
01750   Double_t nStrips  = 0;
01751   Double_t   avePH  = 0;
01752   Double_t aveWidth = 0;
01753   
01754   //keep looping until we have all windows
01755   Bool_t keepGoing = true;
01756   Int_t nLoops = 0;
01757   while(keepGoing && nLoops<MAXHOUGHITER) {
01758     
01759     MSG("SubShowerSR",Msg::kDebug) << "Running BestHough(): Loop number " 
01760                                    << nLoops << endl;
01761 
01762     //remember number of strips in cluster after last loop
01763     //once this value does not change, then stop looping
01764     Double_t oldNStrips = nStrips;
01765 
01766     //Get Hough m,c:
01767     std::vector<Double_t> MCV(14,0.0);
01768     Int_t badGradient[4] = {0};
01769 
01770     //If we are on the first loop, use all windows to get angle
01771     //If we are not on first loop, just use previously 
01772     //selected windows to get better determination of angle
01773     if(cluster.size()==0){
01774       MSG("SubShowerSR",Msg::kDebug) << "cluster.size()=0" << endl;
01775       MCV = BestHough(win,false);
01776 
01777       //verify that best hough passes reasonably close to vertex:
01778       if(true){
01779         MSG("SubShowerSR",Msg::kDebug) << "vertexPlane = " 
01780                                        << vertexPlane << endl;
01781         MSG("SubShowerSR",Msg::kDebug) << "vertexTPos = " 
01782                                        << vertexTPos << endl;
01783         MSG("SubShowerSR",Msg::kDebug) << "vertexTPosWidth = " 
01784                                        << vertexTPosWidth << endl;
01785         
01786         double vertexY_simple = (vertexPlane - MCV[2])*MCV[0] + MCV[1];
01787         double vertexY_simple_peak = (vertexPlane - MCV[2])*MCV[10] + MCV[11];
01788         
01789         MSG("SubShowerSR",Msg::kDebug) << "vertexY_simple = " 
01790                                        << vertexY_simple << endl;
01791         
01792         if( vertexY_simple < vertexTPos - vertexTPosWidth ||
01793             vertexY_simple > vertexTPos + vertexTPosWidth ||
01794             MCV[0]==0 ) badGradient[0] = 1;
01795         if( vertexY_simple_peak < vertexTPos - vertexTPosWidth || 
01796             vertexY_simple_peak > vertexTPos + vertexTPosWidth || 
01797             MCV[10]==0 ) badGradient[1] = 1;
01798     
01799         double vertexY_ph      = (vertexPlane - MCV[7])*MCV[5] + MCV[6];
01800         double vertexY_ph_peak = (vertexPlane - MCV[7])*MCV[12] + MCV[13];
01801         MSG("SubShowerSR",Msg::kDebug) << "vertexY_ph = " 
01802                                        << vertexY_ph << endl;
01803         if( vertexY_ph < vertexTPos - vertexTPosWidth || 
01804             vertexY_ph > vertexTPos + vertexTPosWidth ||
01805             MCV[5]==0 ) badGradient[2] = 1;
01806         if( vertexY_ph_peak < vertexTPos - vertexTPosWidth || 
01807             vertexY_ph_peak > vertexTPos + vertexTPosWidth ||
01808             MCV[12]==0 ) badGradient[3] = 1;
01809       }
01810     }
01811     else {
01812       MSG("SubShowerSR",Msg::kDebug) << "cluster.size()=" 
01813                                      << cluster.size() << endl;
01814       MCV = BestHough(cluster,false);
01815     }
01816 
01817     MSG("SubShowerSR",Msg::kDebug) 
01818       << "\nHough gradient, intercept, plane vertex: " 
01819       << MCV[0] << "+/-" << MCV[3] << " " 
01820       << MCV[1] << "+/-" << MCV[4] << " " << MCV[2]
01821       << "\nPH Hough gradient, intercept, plane vertex: " 
01822       << MCV[5] << "+/-" << MCV[8] << " " 
01823       << MCV[6] << "+/-" << MCV[9] << " " << MCV[7]
01824       << "\n Hough Max Bin Gradient and Intercept Values: "
01825       << MCV[10] << " " << MCV[11] << " "
01826       << MCV[12] << " " << MCV[13] << endl;
01827     
01828     //if we've learnt nothing from first run of BestHough
01829     //stop now and try other clustering algorithm
01830     if(badGradient[0]==1 && badGradient[1]==1 && 
01831        badGradient[2]==1 && badGradient[3]==1 && 
01832        cluster.size()==0) return emptyCluster;
01833 
01834     cluster.clear();
01835     avePH = 0;
01836     aveWidth = 0;
01837     nStrips = 0;
01838     lowPlane = 500;
01839     uppPlane = 0;
01840     std::vector<int> tempcluster[4];
01841     std::vector<window>::iterator winIter = nonCluster.begin();
01842     std::vector<window>::iterator winEnd  = nonCluster.end();
01843     Int_t counter = 0;
01844     while(winIter!=winEnd){
01845       double y[4] = {};
01846       double y_upp[4] = {};
01847       double y_low[4] = {};
01848       y[0] = (winIter->plane - MCV[2])*MCV[0]  + MCV[1];
01849       y[1] = (winIter->plane - MCV[2])*MCV[10] + MCV[11];
01850       y[2] = (winIter->plane - MCV[7])*MCV[5]  + MCV[6];
01851       y[3] = (winIter->plane - MCV[7])*MCV[12] + MCV[13];
01852       
01853       if(cluster.size()==0) {
01854         for(int ii=0;ii<4;ii++) {
01855           y_upp[ii] = y[ii] + STRIPWIDTHINMETRES; 
01856           y_low[ii] = y[ii] - STRIPWIDTHINMETRES;
01857         }
01858       }
01859       else {
01860         for(int ii=0;ii<4;ii++) {
01861           y_upp[ii] = y[ii] + MCV[4+int(ii/2)*5];
01862           y_low[ii] = y[ii] - MCV[4+int(ii/2)*5];
01863         }
01864       }
01865       
01866       for(int ii=0;ii<4;ii++){
01867         if( ( winIter->upper_tpos>=y[ii] && winIter->lower_tpos<=y[ii] ) ||
01868             ( winIter->phpos<=y_upp[ii]  && winIter->phpos>=y_low[ii]  ) ) {
01869           if(badGradient[ii]==0) {
01870             tempcluster[ii].push_back(counter);
01871           }
01872         }
01873       }
01874       winIter++;
01875       counter++;
01876     }
01877     
01878     Int_t best_cluster = -1;
01879     UInt_t max_size = 0;
01880     for(int ii=0;ii<4;ii++){
01881       if(tempcluster[ii].size()>max_size) {
01882         max_size = tempcluster[ii].size();
01883         best_cluster = ii;
01884       }
01885     }
01886     if(best_cluster>=0 && best_cluster<=3) {
01887       std::vector<window>::iterator winIter = nonCluster.begin();
01888       std::vector<window>::iterator winEnd  = nonCluster.end();
01889       std::vector<int>::iterator tempIter = tempcluster[best_cluster].begin();
01890       std::vector<int>::iterator tempEnd  = tempcluster[best_cluster].end();
01891       counter = 0;
01892       while(winIter!=winEnd){
01893         if(tempIter!=tempEnd && (*tempIter) == counter) {
01894           cluster.push_back(*winIter);
01895           avePH += winIter->ph;
01896           aveWidth += winIter->upper_tpos - winIter->lower_tpos;
01897           nStrips += Double_t(winIter->upper - winIter->lower + 1);
01898           if(winIter->plane<lowPlane) lowPlane = winIter->plane;
01899           if(winIter->plane>uppPlane) uppPlane = winIter->plane;
01900           tempIter++;
01901         }
01902         winIter++;
01903         counter++;
01904       }
01905     }
01906     
01907     MSG("SubShowerSR",Msg::kDebug) << "Proto-SubShower currently has " 
01908                                    << (uppPlane-lowPlane+2)/2
01909                                    << " planes and " << nStrips 
01910                                    << " strips" << endl;
01911 
01912     if(nStrips == oldNStrips) keepGoing = false;
01913     nLoops+=1;
01914   }
01915 
01916   if(nStrips<=0) return emptyCluster;
01917     
01918   MSG("SubShowerSR",Msg::kDebug) 
01919     << "Returning Proto-SubShower from HoughTransCluster()" << endl;
01920   
01921   if( det==1 && uppPlane>249 && lowPlane<249 ) return cluster;
01922   avePH/=nStrips;
01923   if( ((uppPlane - lowPlane + 2)/2)<fMinHoughPlanes && 
01924       avePH<fMinHoughPH ) return emptyCluster;
01925   
01926   return cluster;
01927 
01928 }

TObjArray * AlgSubShowerSRList::MakeClusters ( std::vector< CandStripHandle * > &  remainingHits,
TObjArray *  shower 
) [private]

Definition at line 3141 of file AlgSubShowerSRList.cxx.

References AreNeighbours().

Referenced by SplitSubShower().

03142 {
03143   for (int j = 0; j<shower->GetEntries(); j++) {
03144     //loop over all strips currently in the shower
03145     //each time a new one is added it is checked
03146     CandStripHandle* stp = static_cast<CandStripHandle*>(shower->At(j));
03147     for(int i = remainingHits.size()-1; i >= 0 ; --i){
03148       //check strip for proximity for all other unassigned strips
03149       //remove strip from unassigned list if added to shower
03150       CandStripHandle* otherStp = remainingHits.at(i);
03151       if (AreNeighbours(stp, otherStp)) {
03152         shower->Add(otherStp);
03153         remainingHits.erase(remainingHits.begin()+i);
03154       }
03155     }
03156   }
03157   return shower;
03158 }

void AlgSubShowerSRList::RunAlg ( AlgConfig ac,
CandHandle ch,
CandContext cx 
) [virtual]

Implements AlgBase.

Definition at line 87 of file AlgSubShowerSRList.cxx.

References CandHandle::AddDaughterLink(), CleanUp(), fBestHoughMaxPHFrac, fCleanUpPlaneCut, fCleanUpStripCut, fCleanUpTimeCut, fFormHalo, fHaloMaxPlaneGap, fHaloNPeakFindingBins, FindCluster(), FIRSTNDSPECPLANE, fLongWindowGapPlaneCut, fLongWindowMipCut, fMaxShowerGap, fMinHoughPH, fMinHoughPlanes, fMinStripPE, FormHalo(), fSplitShowers, Registry::Get(), AlgFactory::GetAlgHandle(), CandContext::GetCandRecord(), CandContext::GetDataIn(), CandHandle::GetDaughterIterator(), VldContext::GetDetector(), Registry::GetDouble(), AlgFactory::GetInstance(), Registry::GetInt(), Calibrator::GetMIP(), CandContext::GetMom(), CandRecoHandle::GetNStrip(), RecMinos::GetVldContext(), Calibrator::Instance(), Msg::kDebug, Msg::kError, Detector::kFar, Msg::kInfo, CalDigitType::kPE, CalDigitType::kSigCorr, PlaneView::kU, PlaneView::kV, Msg::kVerbose, PlaneView::kX, PlaneView::kY, CandSubShowerSR::MakeCandidate(), MSG, CandRecoHandle::SetCandSlice(), CandSubShowerSRHandle::SetMinStripPE(), SplitSubShower(), TestOverLap(), and TimeTest().

00088 {
00089 
00090   MSG("SubShowerSR", Msg::kDebug) << "Starting AlgSubShowerSRList::RunAlg()" 
00091                                   << endl;
00092   assert(cx.GetDataIn());
00093   if (!(cx.GetDataIn()->InheritsFrom("TObjArray"))) {
00094     return;
00095   }
00096   
00097   const CandSliceListHandle *slicelist = 0;
00098   const CandClusterListHandle *clusterlist = 0;
00099   const CandTrackListHandle *tracklist = 0;
00100   const TObjArray *cxin = dynamic_cast<const TObjArray *>(cx.GetDataIn());
00101   for (Int_t i=0; i<=cxin->GetLast(); i++) {
00102     TObject *tobj = cxin->At(i);
00103     MSG("SubShowerSR", Msg::kDebug) << "cxin->At(" << i << ")"<< endl;
00104     if (tobj->InheritsFrom("CandSliceListHandle")) {
00105       slicelist = dynamic_cast<CandSliceListHandle*>(tobj);
00106       MSG("SubShowerSR", Msg::kDebug) << "Got SliceList" << endl;
00107     }
00108     if (tobj->InheritsFrom("CandClusterListHandle")) {
00109       clusterlist = dynamic_cast<CandClusterListHandle*>(tobj);
00110       MSG("SubShowerSR", Msg::kDebug) << "Got ClusterList" << endl;
00111     }
00112     if (tobj->InheritsFrom("CandTrackListHandle")) {
00113       tracklist = dynamic_cast<CandTrackListHandle*>(tobj);
00114       MSG("SubShowerSR", Msg::kDebug) << "Got TrackList" << endl;
00115     }
00116   }
00117   if (!slicelist || !clusterlist || !tracklist) {
00118     MSG("SubShowerSR",Msg::kInfo) << "CandSliceListHandle or " 
00119                                   << "CandClusterListHandle or " 
00120                                   << "CandTrackListHandle missing\n";
00121     return;
00122   }
00123 
00124   MSG("SubShowerSR", Msg::kDebug) << "Creating CandContext" << endl;
00125   // Create Candcontext
00126   CandContext cxx(this,cx.GetMom());
00127 
00128   //get config for AlgSubShowerSRList
00129   const char *charSubShowerSRAlgConfig = 0;
00130   ac.Get("SubShowerSRAlgConfig",charSubShowerSRAlgConfig);
00131 
00132   Int_t minNTrkOnlyPlanesToUseTrackStrips = ac.GetInt("MinNTrkOnlyPlanesToUseTrackStrips");
00133   Int_t maxPlanesFromVtxToUseTrackStrips  = ac.GetInt("MaxPlanesFromVertexToUseTrackStrips");
00134   Double_t minMIPsToUseTrackStrips        = ac.GetDouble("MinMIPsToUseTrackStrips");
00135 
00136   fLongWindowMipCut      = ac.GetDouble("LongWindowMipCut");
00137   fLongWindowGapPlaneCut = ac.GetInt("LongWindowGapPlaneCut");
00138   fMinHoughPlanes        = ac.GetInt("MinHoughPlanes");
00139   fMinHoughPH            = ac.GetDouble("MinHoughPH");
00140   fBestHoughMaxPHFrac    = ac.GetDouble("BestHoughMaxPHFrac");
00141   fCleanUpTimeCut        = ac.GetDouble("CleanUpTimeCut");
00142   fCleanUpPlaneCut       = ac.GetInt("CleanUpPlaneCut");
00143   fCleanUpStripCut       = ac.GetDouble("CleanUpStripCut");
00144   fFormHalo              = ac.GetInt("FormHalo");
00145   fHaloMaxPlaneGap       = ac.GetInt("HaloMaxPlaneGap");
00146   fHaloNPeakFindingBins  = ac.GetInt("HaloNPeakFindingBins");
00147   fMinStripPE            = ac.GetDouble("MinStripPE");
00148   fMaxShowerGap          = ac.GetDouble("MaxShowerGap");  
00149   fSplitShowers          = ac.GetInt("SplitShowers");  
00150 
00151   //Get singleton instance of AlgFactory
00152   AlgFactory &af = AlgFactory::GetInstance();
00153   AlgHandle ah = af.GetAlgHandle("AlgSubShowerSR",charSubShowerSRAlgConfig);
00154                                                                    
00155   const CandRecord *candrec = cx.GetCandRecord();
00156   assert(candrec);
00157   const VldContext *vldcptr = candrec->GetVldContext();
00158   assert(vldcptr);
00159   VldContext vldc = *vldcptr;                                         
00160   UgliGeomHandle ugh(vldc);
00161   Int_t detector = 0;
00162   if(vldc.GetDetector()==Detector::kFar) detector = 1;
00163 
00164   Calibrator& calibrator = Calibrator::Instance();
00165 
00166   //keep track of all unassigned strips in all slices and try to add them later
00167   TObjArray unassignedStrips;
00168 
00169   //loop over all slices
00170   Int_t sliceCounter = 0;
00171   CandSliceHandleItr sliceItr(slicelist->GetDaughterIterator());
00172   while (CandSliceHandle *slice = sliceItr()) {
00173 
00174     MSG("SubShowerSR",Msg::kVerbose) << "Slice " <<  sliceCounter << endl;
00175 
00176     //also keep track of cluster planes and track planes
00177     Int_t cluPln[500] = {0};
00178     Int_t trkPln[500] = {0};
00179     Int_t nTrkOnlyPlanes = 0;
00180     if(tracklist!=0) {
00181       //loop through all tracks in track list
00182       CandTrackHandleItr trackItr(tracklist->GetDaughterIterator());
00183       while (CandTrackHandle *track = trackItr()){
00184         //if this track is in the current slice then...
00185         if (*track->GetCandSlice()==*slice) {
00186           for(Int_t i=track->GetBegPlane();i<=track->GetEndPlane();i++){
00187             trkPln[i] += 1;   //if trkPln > 1 => multiple tracks in plane
00188             nTrkOnlyPlanes += 1;
00189           }
00190         }
00191       }
00192     }
00193     if(clusterlist!=0) {
00194       //loop through all clusters in cluster list
00195       CandClusterHandleItr clusterItr(clusterlist->GetDaughterIterator());
00196       while (CandClusterHandle *cluster = clusterItr()){
00197         //if this cluster is in the current slice then...
00198         if (*cluster->GetCandSlice()==*slice) {
00199           for(Int_t i=cluster->GetBegPlane();i<=cluster->GetEndPlane();i+=2){
00200             cluPln[i] += 1;   //if cluPln > 1 => multiple clusters in plane
00201             if(trkPln[i]>0) nTrkOnlyPlanes-=1;
00202             if(detector==1){//account for SM gap
00203               if(i==247) i = 250; 
00204               if(i==248) i = 251;
00205             }
00206           }
00207         }
00208       }
00209     }
00210 
00211     //for each slice, make an Object Array to hold all the 
00212     //non track strips in the slice
00213     TObjArray allUStrips;
00214     TObjArray allVStrips;
00215     
00216     //loop through the strips in the slice
00217     CandStripHandleItr stripslcItr(slice->GetDaughterIterator());    
00218     MSG("SubShowerSR",Msg::kVerbose) << "About to loop through slice strips" 
00219                                      << endl;
00220     
00221     while (CandStripHandle *stp =stripslcItr()) { //while over slice stp
00222 
00223       // Only consider strips above the required minimum
00224       if (stp->GetCharge(CalDigitType::kPE) < fMinStripPE) continue;
00225 
00226       bool gotStrip = false;
00227       MSG("SubShowerSR",Msg::kVerbose) << "PlaneView " 
00228                                        << stp->GetPlaneView() << endl;
00229       
00230       CandClusterHandleItr clusterItr(clusterlist->GetDaughterIterator());
00231       Int_t clusterCounter = 0;
00232       while(CandClusterHandle *cluster = clusterItr()){ 
00233         MSG("SubShowerSR",Msg::kVerbose) 
00234           << "Looking for stp match in cluster " << clusterCounter << endl;
00235         if(*cluster->GetCandSlice()==*slice){
00236           CandStripHandleItr stripcluItr(cluster->GetDaughterIterator());
00237           while (CandStripHandle *stpclu =stripcluItr()){
00238             if(*stp==*stpclu){
00239               gotStrip = true;
00240               //add strip to object array then break out of loop
00241               if(stp->GetPlaneView()==PlaneView::kU ||
00242                  stp->GetPlaneView()==PlaneView::kX) {
00243                 allUStrips.Add(stp);
00244                 MSG("SubShowerSR",Msg::kVerbose) << "Found match! U strip" 
00245                                                  << endl;             
00246               }
00247               else if(stp->GetPlaneView()==PlaneView::kV ||
00248                       stp->GetPlaneView()==PlaneView::kY) {
00249                 allVStrips.Add(stp);
00250                 MSG("SubShowerSR",Msg::kVerbose) << "Found match! V strip" 
00251                                                  << endl;
00252               }
00253               else MSG("SubShowerSR",Msg::kError) 
00254                 << "Unknown PlaneView! Not adding strip to initial ObjArray" 
00255                 << endl;
00256               break;
00257             }
00258           }
00259         }
00260         if(gotStrip) break;
00261         clusterCounter++;
00262       }
00263       
00264       //if this is an ND strip in spectrometer, don't add to array
00265       if(!gotStrip && detector==0 && stp->GetPlane()>=FIRSTNDSPECPLANE) gotStrip=true;
00266 
00267       if(!gotStrip){  //if !gotStrip
00268         bool inTrack = false;
00269         Int_t nPlanesFromBegPlane = 0;
00270         if(tracklist!=0) {
00271           MSG("SubShowerSR",Msg::kVerbose) 
00272             << "No matched strip found in clusters, looking in tracks"
00273             << endl;
00274           
00275           //loop through all tracks in track list
00276           CandTrackHandleItr trackItr(tracklist->GetDaughterIterator());
00277           Int_t trackCounter = 0;
00278           while (CandTrackHandle *track = trackItr()){
00279             MSG("SubShowerSR",Msg::kVerbose) 
00280               << "Looking for match of slc stp with trk stp in track "
00281               << trackCounter << endl;
00282             //if this track is in the current slice then...
00283             if (*track->GetCandSlice()==*slice) {
00284               //...loop over all strips in this track
00285               CandStripHandleItr striptrkItr(track->GetDaughterIterator());
00286               while (CandStripHandle *stptrk =striptrkItr()){
00287                 //if this strip is in this track then...
00288                 if (*stp == *stptrk) {
00289                   MSG("SubShowerSR",Msg::kVerbose) 
00290                     << "Found match! This strip is in a track!" << endl;
00291                   //...break and move on to the next strip in the slice
00292                   inTrack = true;
00293                   nPlanesFromBegPlane = TMath::Abs(stp->GetPlane()-track->GetBegPlane());
00294                   break;
00295                 }
00296               }
00297             }
00298             if(inTrack) break;
00299             trackCounter++;
00300           }
00301         }
00302         
00303         MSG("SubShowerSR",Msg::kVerbose) << "nTrkOnlyPlanes = " 
00304                                          << nTrkOnlyPlanes << endl;
00305         if(inTrack && nTrkOnlyPlanes<minNTrkOnlyPlanesToUseTrackStrips) {
00306           MSG("SubShowerSR",Msg::kDebug) 
00307             << "nTrkOnlyPlanes<" << minNTrkOnlyPlanesToUseTrackStrips 
00308             << " - going to include this strip " 
00309             << "even though it's in a track." << endl;
00310           inTrack = false; //ignore tracks that are buried in showers
00311         }
00312         else if(inTrack && nPlanesFromBegPlane <
00313                 maxPlanesFromVtxToUseTrackStrips) { //near vertex
00314           if(calibrator.GetMIP(stp->GetCharge(CalDigitType::kSigCorr)) >
00315              minMIPsToUseTrackStrips) {
00316             MSG("SubShowerSR",Msg::kDebug) 
00317               << "nPlanesFromBegPlane < "
00318               << maxPlanesFromVtxToUseTrackStrips 
00319               << " for this strip and ph = "
00320               << calibrator.GetMIP(stp->GetCharge(CalDigitType::kSigCorr))
00321               << " mips. "
00322               << "Going to include it even though it's in a track." 
00323               << endl;
00324             inTrack = false;
00325           }
00326         }
00327 
00328         //then: add it to our array
00329         if(!inTrack){
00330           MSG("SubShowerSR",Msg::kVerbose) << "Strip not found in track! ";
00331           if(stp->GetPlaneView()==PlaneView::kU ||
00332              stp->GetPlaneView()==PlaneView::kX) {
00333             allUStrips.Add(stp);
00334             MSG("SubShowerSR",Msg::kVerbose) << "Adding U strip" << endl;
00335           }
00336           else if(stp->GetPlaneView()==PlaneView::kV ||
00337                   stp->GetPlaneView()==PlaneView::kY) {
00338             allVStrips.Add(stp);
00339             MSG("SubShowerSR",Msg::kVerbose) << "Adding V strip" << endl;
00340           }
00341           else MSG("SubShowerSR",Msg::kError) 
00342             << "Unknown PlaneView! Not adding strip to initial ObjArray" 
00343             << endl;
00344         }
00345       }
00346     }
00347 
00348     MSG("SubShowerSR",Msg::kDebug) << "Found total of " 
00349                                    << allUStrips.GetLast()+1 << " U strips" 
00350                                    << " and " << allVStrips.GetLast()+1 
00351                                    << " V strips" << endl;
00352     
00353     //In here put code to pull out the strips belonging to the SubShowers:
00354     bool keepGoingU = true;
00355     bool keepGoingV = true;
00356     while(keepGoingU || keepGoingV){
00357       if(keepGoingU) 
00358         MSG("SubShowerSR",Msg::kDebug) << "Still have " 
00359                                        << allUStrips.GetLast()+1
00360                                        << " U strips to process" << endl;
00361       
00362       else if(keepGoingV) 
00363         MSG("SubShowerSR",Msg::kDebug) << "Still have " 
00364                                        << allVStrips.GetLast()+1
00365                                        << " V strips to process" << endl;
00366 
00367       if(allUStrips.GetLast()+1<=0) keepGoingU = false;
00368       if(allVStrips.GetLast()+1<=0) keepGoingV = false;
00369 
00370       TObjArray newSubShower;
00371       if(keepGoingU){ //do U first
00372         this->FindCluster(&allUStrips,&newSubShower,detector);
00373         this->TimeTest(&newSubShower);
00374         MSG("SubShowerSR",Msg::kDebug) << "Added " << newSubShower.GetLast()+1 
00375                                        << " U strips to newSubShower" << endl;
00376       }
00377       else if(keepGoingV) { //then V strips
00378         this->FindCluster(&allVStrips,&newSubShower,detector);
00379         this->TimeTest(&newSubShower);
00380         MSG("SubShowerSR",Msg::kDebug) << "Added " << newSubShower.GetLast()+1 
00381                                        << " V strips to newSubShower" << endl;
00382      }
00383 
00384       //Everytime you have a newSubShower with all the strips 
00385       //you want to add do this
00386       //if splitShowers is turned on, create a vector of split subshowers
00387       //otherwise, just put the full subshower in the vector
00388       std::vector<TObjArray*> showers;
00389       if (fSplitShowers) showers = this->SplitSubShower(&newSubShower);
00390       else showers.push_back(&newSubShower);
00391       
00392       for (UInt_t i = 0; i<showers.size(); i++){
00393         TObjArray* splitSubShower = showers.at(i);
00394         while(this->TestOverLap(splitSubShower,ch,slice)){};
00395         cxx.SetDataIn(splitSubShower);
00396         MSG("SubShowerSR",Msg::kDebug) << "forming SubShowerSR\n";
00397         CandSubShowerSRHandle subshowersrhandle = CandSubShowerSR::MakeCandidate(ah,cxx);
00398         // Set the strip PE cut - this gets passed through the subshower
00399         // to the shower, so that it is accessible from the truthhelper
00400         // when it fills the purity and completeness variables
00401         subshowersrhandle.SetMinStripPE(fMinStripPE);
00402         
00403         if(subshowersrhandle.GetNStrip()!=0) {
00404           
00405           MSG("SubShowerSR",Msg::kDebug) << "Formed SubShowerSR with " 
00406                                          << subshowersrhandle.GetNStrip() 
00407                                          << " strips" << endl;
00408           
00409           subshowersrhandle.SetCandSlice(slice);
00410           ch.AddDaughterLink(subshowersrhandle);
00411           
00412           CandStripHandleItr stripItr(subshowersrhandle.GetDaughterIterator());
00413           
00414           std::vector<Int_t> tmpUArray(allUStrips.GetLast()+1,0);
00415           for(int i=0;i<allUStrips.GetLast()+1;i++) tmpUArray[i] = 0;
00416           std::vector<Int_t> tmpVArray(allVStrips.GetLast()+1,0);
00417           for(int i=0;i<allVStrips.GetLast()+1;i++) tmpVArray[i] = 0;   
00418         
00419           while(CandStripHandle *stp = stripItr()){
00420             if(stp->GetPlaneView()==PlaneView::kU ||
00421                stp->GetPlaneView()==PlaneView::kX){
00422               for (int i=0; i<=allUStrips.GetLast(); i++) {
00423                 TObject *tobj = allUStrips.At(i);
00424                 assert(tobj->InheritsFrom("CandStripHandle"));
00425                 CandStripHandle *striphandle = 
00426                   dynamic_cast<CandStripHandle*>(tobj);
00427                 if(*stp==*striphandle) {
00428                   tmpUArray[i] = 1;
00429                   break;
00430                 }
00431               }
00432             }
00433             else if(stp->GetPlaneView()==PlaneView::kV ||
00434                     stp->GetPlaneView()==PlaneView::kY){
00435               for (int i=0; i<=allVStrips.GetLast(); i++) {
00436                 TObject *tobj = allVStrips.At(i);
00437                 assert(tobj->InheritsFrom("CandStripHandle"));
00438                 CandStripHandle *striphandle = 
00439                   dynamic_cast<CandStripHandle*>(tobj);
00440                 if(*stp==*striphandle) {
00441                   tmpVArray[i] = 1;
00442                   break;
00443                 }
00444               }
00445             }
00446           }
00447           
00448           //remove used strips from arrays
00449           for(int i=0;i<allUStrips.GetLast()+1;i++) {
00450             if(tmpUArray[i]==1) allUStrips.RemoveAt(i);
00451           }
00452           allUStrips.Compress();
00453           if(allUStrips.GetLast()==-1) keepGoingU = false;
00454           
00455           for(int i=0;i<allVStrips.GetLast()+1;i++) {
00456             if(tmpVArray[i]==1) allVStrips.RemoveAt(i);
00457           }
00458           allVStrips.Compress();
00459           if(allVStrips.GetLast()==-1) keepGoingV = false;
00460           
00461         }
00462         else {
00463           MSG("SubShowerSR",Msg::kDebug) 
00464             << "SubShowerSR formed with 0 strips - not adding to list" 
00465             << endl;
00466           if(keepGoingU) {
00467             MSG("SubShowerSR",Msg::kDebug) << "Running CleanUp() on U Strips" 
00468                                            << endl;     
00469             keepGoingU = CleanUp(&allUStrips,ch,detector,slice);
00470             if(!keepGoingU) {
00471               for(int i=0;i<allUStrips.GetLast()+1;i++) {
00472                 unassignedStrips.Add(allUStrips.At(i));
00473               }
00474             }
00475           }
00476           else {
00477             MSG("SubShowerSR",Msg::kDebug) << "Running CleanUp() on V Strips" 
00478                                            << endl;
00479             keepGoingV = CleanUp(&allVStrips,ch,detector,slice);
00480             if(!keepGoingV) {
00481               for(int i=0;i<allVStrips.GetLast()+1;i++) {
00482                 unassignedStrips.Add(allVStrips.At(i));
00483               }
00484             }
00485           }
00486         }
00487       }
00488     }
00489     sliceCounter++;
00490   }
00491   
00492   //slice-wise cleanup
00493   if(fFormHalo&&FormHalo(&unassignedStrips,ch,cxx,ah,slicelist,detector)) 
00494     MSG("SubShowerSR",Msg::kDebug) << "Halo SubShowers formed" << endl;  
00495 
00496 }

std::vector< TObjArray * > AlgSubShowerSRList::SplitSubShower ( const TObjArray *  subShower  )  [private]

Definition at line 3110 of file AlgSubShowerSRList.cxx.

References MakeClusters().

Referenced by FormHalo(), and RunAlg().

03111 {
03112   std::vector<TObjArray*> splitSubShowers;
03113   Int_t nStrips = subShower->GetEntries();
03114 
03115   if (nStrips<=1){ //in case subshower has only one or zero hits
03116     splitSubShowers.push_back(new TObjArray(*subShower));
03117     return splitSubShowers;
03118   }
03119 
03120   else{
03121     //copy subshower elements into a vector of hits
03122     std::vector<CandStripHandle*> remainingHits;
03123     for (int i = 0; i<nStrips; i++){
03124       remainingHits.push_back(static_cast<CandStripHandle*>(subShower->At(i)));
03125     }
03126 
03127     //now cluster from the hits in this vector, removing as we go
03128     while(remainingHits.size()>0){
03129       TObjArray* newSubShower = new TObjArray();
03130       Int_t endIndex = remainingHits.size()-1;
03131       CandStripHandle* stp = remainingHits.at(endIndex);
03132       newSubShower->Add(stp);
03133       remainingHits.erase(remainingHits.begin()+endIndex);
03134       newSubShower = MakeClusters(remainingHits, newSubShower);
03135       splitSubShowers.push_back(newSubShower);
03136     }
03137     return splitSubShowers;
03138   }
03139 }

NavKey AlgSubShowerSRList::StripKeyFromPlane ( const CandStripHandle strip  )  [static]

Definition at line 3172 of file AlgSubShowerSRList.cxx.

References VHS::GetPlane().

03173 {
03174   return const_cast<CandStripHandle *>(strip)->GetPlane();
03175 }

Bool_t AlgSubShowerSRList::TestOverLap ( TObjArray *  newSubShower,
CandHandle ch,
const CandSliceHandle slice 
) [private]

Definition at line 1931 of file AlgSubShowerSRList.cxx.

References CandStripHandle::GetCharge(), VHS::GetPlane(), CandStripHandle::GetPlane(), CandStripHandle::GetTPos(), CandStripHandle::GetZPos(), Msg::kDebug, CandSubShowerSRHandle::KeyFromViewEnergy(), ClusterType::kHalo, CalDigitType::kSigCorr, PlaneView::kU, PlaneView::kV, PlaneView::kX, ClusterType::kXTalk, PlaneView::kY, and MSG.

Referenced by RunAlg().

01933 {
01934   //Going to check that proto-subshower does not go through an existing subshower  
01935   //if it does, take largest, plane-wise contiguous set of strips
01936 
01937   MSG("SubShowerSR", Msg::kDebug) << "Running TestOverLap()" << endl;
01938 
01939   //if no strips return false. Empty subshower will be caught later.
01940   Int_t totStp = newSubShower->GetLast()+1;  
01941   if(totStp<=1) return false;
01942 
01943   assert(ch.InheritsFrom("CandSubShowerSRListHandle"));
01944   CandSubShowerSRListHandle &csslh = dynamic_cast<CandSubShowerSRListHandle &>(ch);
01945   //if this is the first subshower, just return false:
01946   if(csslh.GetNDaughters()==0) return false;
01947 
01948   //Get plane view:
01949   TObject *firstObj = newSubShower->At(0);
01950   assert(firstObj->InheritsFrom("CandStripHandle"));
01951   CandStripHandle *firstStp = dynamic_cast<CandStripHandle*>(firstObj);
01952   PlaneView::PlaneView_t pv = firstStp->GetPlaneView();
01953   MSG("SubShowerSR", Msg::kDebug) << "PlaneView = " << pv << endl;
01954 
01955   //loop through subshowers in view + energy order
01956   //so that newSubShower is tested against largest subshower first
01957   CandSubShowerSRHandleItr subshowerItr(csslh.GetDaughterIterator());
01958   CandSubShowerSRHandleKeyFunc *engKF = subshowerItr.CreateKeyFunc();
01959   engKF->SetFun(CandSubShowerSRHandle::KeyFromViewEnergy);
01960   subshowerItr.GetSet()->AdoptSortKeyFunc(engKF);
01961   while (CandSubShowerSRHandle *subshower = subshowerItr()) {
01962     if((*subshower->GetCandSlice())!=(*slice)) continue;
01963     if(subshower->GetPlaneView()!=pv) continue;
01964     if(subshower->GetClusterID()==ClusterType::kXTalk ||
01965        subshower->GetClusterID()==ClusterType::kHalo) continue;
01966     //Get vertex and angle of subshower:
01967     Double_t tpos_vtx = 0;
01968     Double_t z_vtx = subshower->GetVtxZ();
01969     if((subshower->GetPlaneView()==PlaneView::kU || 
01970         subshower->GetPlaneView()==PlaneView::kX)) 
01971       tpos_vtx = subshower->GetVtxU();
01972     if((subshower->GetPlaneView()==PlaneView::kV || 
01973         subshower->GetPlaneView()==PlaneView::kY)) {
01974       tpos_vtx = subshower->GetVtxV();
01975     }
01976     Double_t slope = subshower->GetSlope();
01977     Int_t beg_plane = subshower->GetBegPlane();
01978     Int_t end_plane = subshower->GetEndPlane();
01979 
01980     //for calculating asymmetry:
01981     Double_t n_upp = 0;
01982     Double_t n_low = 0;
01983     Double_t ph_upp = 0;
01984     Double_t ph_low = 0;
01985     Double_t n_plane_upp[500] = {0};
01986     Double_t n_plane_low[500] = {0};
01987     Double_t ph_plane_upp[500] = {0};
01988     Double_t ph_plane_low[500] = {0};
01989 
01990     Int_t firstPlane = 500;
01991     Int_t lastPlane = 0;
01992 
01993     //calculate number of strips above and below axis of subshower
01994     for(int i=0;i<totStp;i++){
01995       CandStripHandle *stp = dynamic_cast<CandStripHandle*>(newSubShower->At(i));
01996       if(firstPlane>stp->GetPlane()) firstPlane = stp->GetPlane();
01997       if(lastPlane<stp->GetPlane()) lastPlane = stp->GetPlane();
01998       if(stp->GetTPos() > (slope*(stp->GetZPos()-z_vtx) + tpos_vtx)) {
01999         n_plane_upp[stp->GetPlane()] += 1;
02000         ph_plane_upp[stp->GetPlane()] += stp->GetCharge(CalDigitType::kSigCorr);
02001         if(stp->GetPlane()>=beg_plane && stp->GetPlane()<=end_plane) {
02002           n_upp += 1;
02003           ph_upp += stp->GetCharge(CalDigitType::kSigCorr);
02004         }
02005       }
02006       else {
02007         n_plane_low[stp->GetPlane()] += 1;
02008         ph_plane_low[stp->GetPlane()] += stp->GetCharge(CalDigitType::kSigCorr);
02009         if(stp->GetPlane()>=beg_plane && stp->GetPlane()<=end_plane) {  
02010           n_low += 1;
02011           ph_low += stp->GetCharge(CalDigitType::kSigCorr);
02012         }
02013       }
02014     }
02015 
02016     //if asymmetry < 95% (should be 100%!) in the plane 
02017     //range of the current subshower then only keep the 
02018     //larger part of the new subshower:
02019     Double_t asymmetry = 1;
02020     Double_t asymmetryPH = 1;
02021     if( n_upp + n_low > 1 ) asymmetry = TMath::Abs(n_upp-n_low)/(n_upp+n_low);
02022     if( ph_upp + ph_low > 1 ) asymmetryPH = TMath::Abs(ph_upp-ph_low)/(ph_upp+ph_low);
02023     MSG("SubShowerSR", Msg::kDebug) << "Asymmetry = " << asymmetry
02024                                     << " AsymmetryPH = " << asymmetryPH
02025                                     << " TotStp = " << totStp
02026                                     << " n_upp = " << n_upp
02027                                     << " n_low = " << n_low
02028                                     << " ph_upp = " << ph_upp
02029                                     << " ph_low = " << ph_low
02030                                     << endl;
02031     if(asymmetry < 0.95) {
02032       //we know now that this proto-subshower is cutting through another subshower
02033       //find the greatest number of contiguous planes above or below slope
02034       //and take those strips only
02035 
02036       Int_t bestAsymFirstPlane = firstPlane;
02037       Int_t bestAsymLastPlane = firstPlane;
02038       Int_t AsymFirstPlane = firstPlane;
02039       Int_t currentSense = 0;
02040       Int_t bestSense = 0;
02041 
02042       //loop through all planes and find longest segments
02043       for(int i=firstPlane;i<=lastPlane;i+=2){
02044         //find longest set of contiguous planes with same asym sign:
02045         if(n_plane_upp[i]>n_plane_low[i] || 
02046            ( n_plane_upp[i]==n_plane_low[i] && 
02047              ph_plane_upp[i]>ph_plane_low[i]) ) {
02048           if(currentSense==0) {
02049             AsymFirstPlane = i;
02050             currentSense = 1;
02051           }
02052           else if(currentSense>0) currentSense += 1;
02053           else if(currentSense<0) { //sense has changed:
02054             if(TMath::Abs(currentSense)>TMath::Abs(bestSense)){
02055               bestSense = currentSense;
02056               bestAsymFirstPlane = AsymFirstPlane;
02057               bestAsymLastPlane = i-1;
02058             }
02059             currentSense = 1;
02060             AsymFirstPlane = i;
02061           }
02062         }
02063         else if(n_plane_low[i]>n_plane_upp[i] || 
02064                 ( n_plane_upp[i]==n_plane_low[i] && 
02065                   ph_plane_low[i]>ph_plane_upp[i]) ){
02066           if(currentSense==0) {
02067             AsymFirstPlane = i;
02068             currentSense = -1;
02069           }
02070           else if(currentSense<0) currentSense -= 1;
02071           else if(currentSense>0) {
02072             if(TMath::Abs(currentSense)>TMath::Abs(bestSense)){
02073               bestSense = currentSense;
02074               bestAsymFirstPlane = AsymFirstPlane;
02075               bestAsymLastPlane = i-1;
02076             }
02077             currentSense = -1;
02078             AsymFirstPlane = i;
02079           }
02080         }
02081       }
02082 
02083       if(TMath::Abs(currentSense)>TMath::Abs(bestSense)){
02084         bestSense = currentSense;
02085         bestAsymFirstPlane = AsymFirstPlane;
02086         bestAsymLastPlane = lastPlane;
02087       }
02088 
02089       for(int i=0;i<totStp;i++){
02090         CandStripHandle *stp = dynamic_cast<CandStripHandle*>(newSubShower->At(i));
02091         if(stp->GetPlane()<bestAsymFirstPlane || stp->GetPlane()>bestAsymLastPlane){
02092           newSubShower->RemoveAt(i);
02093         }
02094       }
02095       newSubShower->Compress();
02096       if(newSubShower->GetLast()+1<totStp) {
02097         return true;
02098       }
02099     }
02100   }
02101   return false;
02102 }

Bool_t AlgSubShowerSRList::TimeTest ( TObjArray *  newSubShower,
TObjArray *  allStrips = NULL 
) [private]

Definition at line 2105 of file AlgSubShowerSRList.cxx.

References fCleanUpTimeCut, CandStripHandle::GetCharge(), CandStripHandle::GetPlane(), CandStripHandle::GetPlaneView(), CandStripHandle::GetStrip(), CandStripHandle::GetTime(), Msg::kDebug, CalDigitType::kSigCorr, and MSG.

Referenced by FormHalo(), and RunAlg().

02106 {
02107   //Going to check that proto-subshower does not have hits with extreme times
02108   MSG("SubShowerSR", Msg::kDebug) << "Running TestTime()" << endl;
02109 
02110   //if no strips return false. Empty subshower will be caught later.
02111   Int_t totStp = newSubShower->GetLast()+1;
02112   if(totStp<=1) return false;
02113   
02114   //Look for >50ns gaps in time structure:
02115   std::vector<Double_t> stripTime(totStp,0.0);
02116   std::vector<TSpecReal_t> stripCharge(totStp,0.0);
02117   for(int i=0;i<totStp;i++){
02118     CandStripHandle *stp = dynamic_cast<CandStripHandle*>(newSubShower->At(i));
02119     Int_t nEarlier = 0;
02120     Int_t nRepeats = 0;
02121     for(int j=0;j<totStp;j++){
02122       if(i==j) continue;
02123       CandStripHandle *stp2 = dynamic_cast<CandStripHandle*>(newSubShower->At(j));    
02124       if(stp->GetTime()>stp2->GetTime()) nEarlier+=1;
02125       else if(stp->GetTime()==stp2->GetTime()) {
02126         if(stp->GetPlane()>stp2->GetPlane()) nEarlier+=1;
02127         else if(stp->GetPlane()==stp2->GetPlane()) {
02128           if(stp->GetStrip()>stp2->GetStrip()) nEarlier+=1;
02129           else if(stp->GetStrip()==stp2->GetStrip()) {
02130             if(stp->GetPlaneView()>stp2->GetPlaneView()) nEarlier+=1;
02131             else if(stp->GetPlaneView()==stp2->GetPlaneView()) nRepeats += 1;         
02132           }
02133         }
02134       }
02135     }
02136     //Fill vectors with time and charge information in time order (earliest first)
02137     stripTime[nEarlier] = stp->GetTime()*1e9;
02138     stripCharge[nEarlier] = stp->GetCharge(CalDigitType::kSigCorr);
02139     for(int j=1;j<=nRepeats;j++) {
02140       stripTime[nEarlier+j] = stp->GetTime()*1e9;
02141       stripCharge[nEarlier+j] = stp->GetCharge(CalDigitType::kSigCorr);
02142     }
02143   }
02144 
02145   //Loop through time vector to find gaps greater than 50ns
02146   //count the number of gaps
02147   Int_t nGaps = 0;
02148   for(int i=0;i<totStp-1;i++){
02149     if(stripTime[i+1] - stripTime[i] > fCleanUpTimeCut) {
02150       nGaps+=1;
02151     }
02152   }
02153   //if no gaps, return now
02154   if(nGaps==0) return false;
02155   //otherwise find the largest chunk of charge between the gaps
02156   std::vector<TSpecReal_t> totalCharge(nGaps+1,0.0);
02157   nGaps = 0;
02158   for(int i=0;i<totStp-1;i++) {
02159     totalCharge[nGaps] += stripCharge[i];
02160     if(stripTime[i+1] - stripTime[i] > fCleanUpTimeCut) {
02161       nGaps+=1;
02162     }
02163   }
02164   //don't forget to add in the last strip!
02165   totalCharge[nGaps] += stripCharge[totStp-1];
02166 
02167   //find the biggest chunk:
02168   Int_t biggestChunk = -1;
02169   Double_t biggestCharge = 0;
02170   for(int i=0;i<nGaps+1;i++) {
02171     if(totalCharge[i] > biggestCharge){
02172       biggestChunk = i;
02173       biggestCharge = totalCharge[i];
02174     }
02175   }
02176 
02177   //now go back to time vector to get the min, max times to use
02178   Double_t minTime = -1;
02179   Double_t maxTime = 0;
02180   nGaps = 0;
02181   for(int i=0;i<totStp-1;i++){
02182     if(nGaps>biggestChunk) break;
02183     if(nGaps == biggestChunk) {
02184       if(minTime<0) minTime = stripTime[i];
02185       maxTime = stripTime[i];
02186     }
02187     if(stripTime[i+1] - stripTime[i] > fCleanUpTimeCut) {
02188       nGaps+=1;
02189     }
02190   }
02191 
02192   //loop through strip list again and remove strips outside of time range
02193   for(int i=0;i<totStp;i++){
02194     CandStripHandle *stp = dynamic_cast<CandStripHandle*>(newSubShower->At(i));
02195     if(stp->GetTime()*1e9<minTime - fCleanUpTimeCut || 
02196        stp->GetTime()*1e9>maxTime + fCleanUpTimeCut) {
02197       newSubShower->RemoveAt(i);
02198       if(allStrips) allStrips->Add(stp);
02199       MSG("SubShowerSR",Msg::kDebug) << "Removing Strip From SubShower" << endl;
02200     }
02201   }
02202   newSubShower->Compress();
02203 
02204   if(newSubShower->GetLast()+1<totStp) return true;
02205   return false;
02206 
02207 }

void AlgSubShowerSRList::Trace ( const char *  c  )  const [virtual]

Reimplemented from AlgBase.

Definition at line 3178 of file AlgSubShowerSRList.cxx.

03179 {
03180 }

std::vector< window > AlgSubShowerSRList::TransCluster ( std::vector< window win,
Int_t  det 
) [private]

Definition at line 1016 of file AlgSubShowerSRList.cxx.

References BestHough(), window::centroid, HoughTransCluster(), Msg::kDebug, window::lower, window::lower_tpos, Munits::m, MSG, window::ph, window::phpos, window::phwidth, window::plane, ssPol1, STRIPWIDTHINMETRES, window::type, window::upper, and window::upper_tpos.

Referenced by FindCluster().

01018 {
01019   std::vector<window> cluster;
01020 
01021   Int_t maxPlane = -1;  //calculate maxPlane
01022   Double_t maxPlanePH = 0;
01023   Double_t PlanePH[500] = {0};
01024   std::vector<window>::iterator winIter = win.begin();
01025   std::vector<window>::iterator winEnd = win.end();
01026   while(winIter!=winEnd){
01027     PlanePH[winIter->plane] += winIter->ph;
01028     if(PlanePH[winIter->plane]>maxPlanePH) {
01029       maxPlane = winIter->plane;
01030       maxPlanePH = PlanePH[maxPlane];
01031     }
01032     winIter++;
01033   }
01034 
01035   Int_t maxLongPlane = -1;
01036   Int_t minLongPlane = -1;
01037   Int_t counter = 0;
01038   Int_t delta = 2;
01039   while(PlanePH[maxPlane+counter]>0||(maxPlane+counter)==249) {
01040     maxLongPlane = maxPlane+counter;
01041     delta = 2;
01042     if(det==1&&(maxLongPlane-249)*(maxLongPlane+delta-249)<0)
01043       delta = 3;
01044     counter+=delta;
01045   }
01046   counter = 0;
01047   while(PlanePH[maxPlane+counter]<=0||(maxPlane+counter)==249) {
01048     minLongPlane = maxPlane+counter;
01049     delta = 2;
01050     if(det==1&&(minLongPlane-249)*(minLongPlane-delta-249)<0)
01051       delta = 3;
01052     counter-=delta;
01053   }
01054   Int_t numPreMaxPlanes = maxPlane-minLongPlane+1;
01055   Int_t numPostMaxPlanes = maxLongPlane-maxPlane+1;
01056   Int_t numLongPlanes = numPostMaxPlanes + numPreMaxPlanes - 1;
01057 
01058   //match up windows plane by plane
01059   Double_t shwmid = 0;    // centroid of cluster transversely
01060   Double_t shwwid = 0;    // trans width of widest window in proto-cluster
01061   Double_t maxwid = 0;
01062   Double_t clusterph = 0; // total PH of cluster
01063   Double_t refp = 0;      // reference plane  }  essentially pivot point for 
01064   Double_t refs = 0;      // reference strip  }  cluster
01065   
01066   //long. window info:  
01067   std::vector<Int_t> winpl(numLongPlanes,0);          // plane number 
01068   std::vector<Int_t> winpl0(numLongPlanes,0);         // upper bound
01069   std::vector<Int_t> winpl1(numLongPlanes,0);         // lower bound
01070   std::vector<TSpecReal_t> winpl0_tpos(numLongPlanes,0.0);  // upper bound in tpos
01071   std::vector<TSpecReal_t> winpl1_tpos(numLongPlanes,0.0);  // lower bound in tpos
01072   std::vector<Double_t> winpl2(numLongPlanes,0.0);      // PH
01073   std::vector<Int_t> winpl3(numLongPlanes,0);         // window type
01074   std::vector<Double_t> winpl4(numLongPlanes,0.0);      // width
01075   std::vector<Double_t> winpl5(numLongPlanes,0.0);      // centroid
01076 
01077   Double_t slopen = 0;
01078   Double_t eslopen = 0;
01079   Bool_t trustslopen = false;
01080 
01081   std::vector<window> houghcluster = HoughTransCluster(win,det);
01082 
01083   if(houghcluster.size()==0) { //houghcluster.size()=0 try usual tracking
01084     for (int d = 0;d<numLongPlanes;d++){
01085       int pl = -1;
01086       if(numPreMaxPlanes>numPostMaxPlanes){
01087         if (d<numPreMaxPlanes) pl = maxPlane-d;
01088         else pl = maxPlane+d-numPreMaxPlanes+1;
01089       }
01090       else{
01091         if (d<numPostMaxPlanes) pl = maxPlane+d;
01092         else pl = maxPlane-d+numPostMaxPlanes-1;
01093       }
01094       winpl[d] = 0;
01095       winpl0[d] = -1;
01096       winpl1[d] = -1;
01097       winpl0_tpos[d] = -10;
01098       winpl1_tpos[d] = -10;
01099       winpl2[d] = 0;
01100       winpl3[d] = 0;
01101       winpl4[d] = -1;
01102       winpl5[d] = -1;
01103       // PH for current best matched window:
01104       double weight = 0.;
01105       //deviation of best matched window from prediction of 
01106       //centroid for this plane:
01107       double weight1 = 10000.;
01108       //loop keeps track of # windows already checked on this plane:
01109       int loop = 0;
01110       // counts number of windows with a good match:
01111       int within = 0;
01112       //prediction of centroid for this plane (not using slope in info):
01113       double mid0 = 0;
01114       //prediction of centroid for this plane using slope info where available:
01115       double mid = 0;
01116       
01117       //loop over all in windows
01118       winIter = win.begin();
01119       while (winIter!=winEnd){
01120         if (winIter->plane==pl){  //if it's in current plane:
01121           loop++;
01122           double dev = shwwid;
01123           if (dev<(winIter->upper_tpos-winIter->lower_tpos+STRIPWIDTHINMETRES))
01124             dev = winIter->upper_tpos-winIter->lower_tpos+STRIPWIDTHINMETRES;
01125           if (dev<2.*STRIPWIDTHINMETRES) dev = 2.*STRIPWIDTHINMETRES;
01126           if (shwmid!=0) {
01127             mid0 = shwmid;
01128             mid = shwmid;
01129           }
01130           else {
01131             mid = (winIter->upper+winIter->lower)/2.;
01132             mid0 = mid;
01133           }
01134           
01135           // check if this window is consistent with previously selected 
01136           // windows and has a higher PH than current best match for this
01137           // plane or else if it's shower max plane
01138           if((shwmid==0 ||
01139               (TMath::Abs(winIter->upper_tpos-mid)<dev &&
01140                TMath::Abs(winIter->lower_tpos-mid)<dev))
01141              && winIter->ph>weight){
01142             
01144             //loop through planes in long. window again to do a local 
01145             //scan to ensure best match is found.
01146             //step forward and backward once each if possible
01147           
01148             int loopp = 0;
01149             int loopm = 0;
01150             int withinp = 0;
01151             int withinm = 0;
01152             //double midp0 = 0;
01153             //double midm0 = 0;
01154             double midp = 0;
01155             double midm = 0;
01156             
01157             // step forward:
01158             std::vector<window>::iterator winIter1 = win.begin();
01159             while(winIter1!=winEnd){
01160               if ((((det==1&&(pl-249)*(pl+2-249)>0)||det==0)
01161                    &&winIter1->plane==pl+2)
01162                   ||(det==1&&(pl-249)*(pl+2-249)<0
01163                      &&winIter1->plane==pl+3)){   //go to next same view plane
01164                 loopp++;
01165                 if(shwmid!=0){  //not shower max plane
01166                   midp = shwmid;
01167                   if(TMath::Abs(winIter1->upper_tpos-midp)<dev/2.
01168                      ||TMath::Abs(winIter1->lower_tpos-midp)<dev/2.){  //good match
01169                     withinp++;
01170                   }
01171                 }
01172                 else if(shwmid==0){  //in shower max plane
01173                   midp = (winIter->upper_tpos+winIter->lower_tpos)/2.;
01174                   if(TMath::Abs(winIter1->upper_tpos-midp)<dev/2.
01175                      ||TMath::Abs(winIter1->lower_tpos-midp)<dev/2.){  //good match
01176                     withinp++;  
01177                   }
01178                 }
01179               }
01180               
01181               //step backwards:
01182               if ((((det==1&&(pl-249)*(pl-2-249)>0)||det==0)
01183                    &&winIter1->plane==pl-2)
01184                   ||(det==1&&(pl-249)*(pl-2-249)<0
01185                      &&winIter1->plane==pl-3)){   
01186                 //go to previous same view plane                
01187                 loopm++;
01188                 if(shwmid!=0){   // not shower max plane
01189                   midm = shwmid;
01190                   if(TMath::Abs(winIter1->upper_tpos-midm)<dev
01191                      ||TMath::Abs(winIter1->lower_tpos-midm)<dev){  //good match
01192                     withinm++;
01193                   }
01194                 }
01195                 else if(shwmid==0){  //shower max plane
01196                   midm = (winIter->upper_tpos+winIter->lower_tpos)/2.;
01197                   if(TMath::Abs(winIter1->upper_tpos-midm)<dev
01198                      ||TMath::Abs(winIter1->lower_tpos-midm)<dev){  //good match
01199                     withinm++;
01200                   }
01201                 }
01202               }
01203               winIter1++;
01204             }
01205             if((withinp>0 && withinm>0) || (loopp==0 && withinm>0) ||
01206                (loopm==0 && withinp>0)){  
01207               //success! current window is a good solution
01208               within++;
01209               // write out solution: 
01210               winpl[d] = winIter->plane;
01211               winpl0[d] = winIter->upper;
01212               winpl1[d] = winIter->lower;
01213               winpl0_tpos[d] = winIter->upper_tpos;
01214               winpl1_tpos[d] = winIter->lower_tpos;
01215               winpl2[d] = winIter->ph;
01216               winpl3[d] = winIter->type;
01217               weight = winIter->ph;
01218               weight1 = (winIter->upper_tpos+winIter->lower_tpos)/2.-mid0;
01219               //these will be overwritten if there is another window 
01220               //which is also a good match, but has more PH, 
01221               //ensuring we get best possible match
01222             }
01223           }
01224         }
01225         winIter++;
01226       }
01227       
01228       if(loop*within>0){   //got a matched window
01229         if (shwmid==0) {   // shower max plane
01230           refp = winpl[d];  // set reference point
01231           refs = (winpl0_tpos[d]+winpl1_tpos[d])/2.;
01232         }
01233         //determine widest window:
01234         if(shwwid<(winpl0_tpos[d]-winpl1_tpos[d]+STRIPWIDTHINMETRES)) 
01235           shwwid = winpl0_tpos[d]-winpl1_tpos[d]+STRIPWIDTHINMETRES; 
01236         // average centroid of event so far:
01237         shwmid = (shwmid*clusterph+(winpl0_tpos[d]+winpl1_tpos[d]) * 
01238                   winpl2[d]/2.)/(clusterph+winpl2[d]); 
01239         clusterph+=winpl2[d];
01240         winpl4[d] = shwwid;
01241         winpl5[d] = shwmid;
01242       }
01243     }
01244   
01245     //Fitting the boundaries of windows matched
01246     Double_t slopetn = 0;
01247     Double_t slopebn = 0;
01248     Double_t eslopetn = 0;
01249     Double_t eslopebn = 0;
01250     Double_t chi2tn = 0;
01251     Double_t chi2bn = 0;
01252     int nt = 0;
01253     int nb = 0;
01254 
01255     for (int d = 0;d<numLongPlanes;d++){
01256       int pl = -1;
01257       if(numPreMaxPlanes>numPostMaxPlanes){
01258         if (d<numPreMaxPlanes) pl = maxPlane-d;
01259         else pl = maxPlane+d-numPreMaxPlanes+1;
01260       }
01261       else{
01262         if (d<numPostMaxPlanes) pl = maxPlane+d;
01263         else pl = maxPlane-d+numPostMaxPlanes-1;
01264       }
01265       if(winpl2[d]>0) {nt++; nb++;}
01266       if (maxwid<winpl4[d]) maxwid = winpl4[d];
01267       if (maxwid<2.*STRIPWIDTHINMETRES) maxwid = 2.*STRIPWIDTHINMETRES;
01268     }
01269     if (nt>2 && nb>2) {
01270       TGraphErrors grat(nt);
01271       TGraphErrors grab(nb);
01272       int t = 0;
01273       int b = 0;
01274       for (int d = 0;d<numLongPlanes;d++){
01275         int pl = -1;
01276         if(numPreMaxPlanes>numPostMaxPlanes){
01277           if (d<numPreMaxPlanes) pl = maxPlane-d;
01278           else pl = maxPlane+d-numPreMaxPlanes+1;
01279         }
01280         else{
01281           if (d<numPostMaxPlanes) pl = maxPlane+d;
01282           else pl = maxPlane-d+numPostMaxPlanes-1;
01283         }
01284         if(winpl2[d]>0) {
01285           grat.SetPoint(t,winpl[d],winpl0_tpos[d]);
01286           grat.SetPointError(t,0,1./TMath::Sqrt(winpl2[d]));
01287           t++;
01288           grab.SetPoint(b,winpl[d],winpl1_tpos[d]);
01289           grab.SetPointError(b,0,1./TMath::Sqrt(winpl2[d]));
01290           b++;
01291         }
01292       }
01293 
01294       grat.Fit("ssPol1","NNQ");
01295       slopetn = ssPol1->GetParameter(1);
01296       eslopetn = ssPol1->GetParError(1);
01297       chi2tn = ssPol1->GetChisquare();
01298       grab.Fit("ssPol1","NNQ");
01299       slopebn = ssPol1->GetParameter(1);
01300       eslopebn = ssPol1->GetParError(1);
01301       chi2bn = ssPol1->GetChisquare();
01302       
01303       if ((slopetn+slopebn)!=0. && slopetn*slopebn!=0.){
01304         if(chi2tn/nt<100. && chi2bn/nb<100. &&
01305            ((((eslopetn/TMath::Abs(slopetn)<0.1 && eslopebn/TMath::Abs(slopebn)<0.1
01306                &&chi2tn/nt>0.1 && chi2bn/nb>0.1)) &&
01307              (TMath::Abs(slopetn-slopebn)/TMath::Abs(slopetn+slopebn)<0.25)) ||
01308             TMath::Abs(slopetn-slopebn)==0)) {
01309           slopen = (slopetn+slopebn)/2.;
01310           eslopen = TMath::Sqrt(TMath::Power(eslopetn,2)+TMath::Power(eslopebn,2))/2.;
01311           if(TMath::Abs(slopetn-slopebn)>eslopen) eslopen = TMath::Abs(slopetn-slopebn);
01312           trustslopen = true;
01313         }
01314       }
01315       MSG("SubShowerSR",Msg::kDebug) <<"fit window "<<slopetn<<" "<<slopebn
01316                                      <<" "<<slopen<<" "<<eslopetn<<" "
01317                                      <<eslopebn<<" "<<eslopen<<" "
01318                                      <<chi2tn/nt<<" "<<chi2bn/nb<<" "
01319                                      <<nt<<" "<<nb<<endl;
01320     }
01321     
01322     //check if solution is good:
01323     double usewid = maxwid;
01324     for (Int_t d = 0;d<numLongPlanes;d++){
01325       int pl = -1;
01326       if(numPreMaxPlanes>numPostMaxPlanes){
01327         if (d<numPreMaxPlanes) pl = maxPlane-d;
01328         else pl = maxPlane+d-numPreMaxPlanes+1;
01329       }
01330       else{
01331         if (d<numPostMaxPlanes) pl = maxPlane+d;
01332         else pl = maxPlane-d+numPostMaxPlanes-1;
01333       }
01334     
01335       if(d>0 && (winpl0[d]>=0||winpl1[d]>=0) ){ //if there is a solution
01336         if(!trustslopen && //don't trust new slope
01337            //window is too far away from centroid
01338            (winpl0_tpos[d]-shwmid)*(winpl1_tpos[d]-shwmid) > 
01339            (maxwid/2.)*(maxwid/2.) ){ 
01340           MSG("SubShowerSR",Msg::kDebug) << "removed plane " << winpl[d] 
01341                                          << " " << winpl0[d] << " " 
01342                                          << winpl1[d] << " " 
01343                                          << winpl0_tpos[d] << " " 
01344                                          << winpl1_tpos[d] << endl;
01345           shwmid = (shwmid*clusterph-(winpl0_tpos[d]+winpl1_tpos[d]) * 
01346                     winpl2[d]/2.)/(clusterph-winpl2[d]);
01347           winpl0[d] = -1;
01348           winpl1[d] = -1;
01349           winpl0_tpos[d] = -10;
01350           winpl1_tpos[d] = -10;
01351           winpl2[d] = 0;
01352           winpl3[d] = 0;
01353           winpl4[d] = -1;
01354           winpl5[d] = -1;
01355         }
01356         else if(trustslopen){ 
01357           //this time account for slope on window boundaries:
01358           //try to find another window as a reference point       
01359           //to compare the window in this plane with
01360           int i = 0;  
01361           while(winpl2[i]==0&&i<=d) i++;
01362           if(i==d){
01363             while(winpl2[i]==0&&i<numLongPlanes) i++;
01364           }
01365           if(usewid<TMath::Abs(winpl[d]-winpl[i])*eslopen)  
01366             //find appropriate cut value on window width
01367             usewid = TMath::Abs(winpl[d]-winpl[i])*eslopen; 
01368           //expected centroid for plane:
01369           double expmid = (winpl0_tpos[i] + 
01370                            winpl1_tpos[i])/2.+(winpl[d]-winpl[i])*slopen; 
01371           if((winpl0_tpos[d]-expmid)*(winpl1_tpos[d]-expmid) >
01372              (usewid/2.)*(usewid/2.)){ 
01373             MSG("SubShowerSR",Msg::kDebug) <<"removed plane using slope "
01374                                            <<winpl[d]<<" "<<winpl0[d]<<" "
01375                                            <<winpl1[d]<<" with pl "
01376                                            <<maxPlane<<" "<<winpl[i]
01377                                            <<" "<<expmid<<" "<<usewid<<endl;
01378             shwmid = (shwmid*clusterph-(winpl0_tpos[d]+winpl1_tpos[d]) * 
01379                       winpl2[d]/2.)/(clusterph-winpl2[d]);
01380             clusterph-=winpl2[d];
01381             winpl0[d] = -1;
01382             winpl1[d] = -1;
01383             winpl0_tpos[d] = -10;
01384             winpl1_tpos[d] = -10;
01385             winpl2[d] = 0;
01386             winpl3[d] = 0;
01387             winpl4[d] = -1;
01388             winpl5[d] = -1;
01389           }
01390         }
01391       }
01392     }
01393   }
01394   else { //use hough solution as a starting point:
01395     for (int d = 0;d<numLongPlanes;d++){  
01396       int pl = -1;
01397       if(numPreMaxPlanes>numPostMaxPlanes){
01398         if (d<numPreMaxPlanes) pl = maxPlane-d;
01399         else pl = maxPlane+d-numPreMaxPlanes+1;
01400       }
01401       else{
01402         if (d<numPostMaxPlanes) pl = maxPlane+d;
01403         else pl = maxPlane-d+numPostMaxPlanes-1;
01404       }
01405       winpl[d] = 0;
01406       winpl0[d] = -1;
01407       winpl1[d] = -1;
01408       winpl0_tpos[d] = -10;
01409       winpl1_tpos[d] = -10;
01410       winpl2[d] = 0;
01411       winpl3[d] = 0;
01412       winpl4[d] = -1;
01413       winpl5[d] = -1;
01414 
01415       //loop over all in windows
01416       std::vector<window>::iterator hcIter = houghcluster.begin();
01417       std::vector<window>::iterator hcEnd = houghcluster.end();      
01418       while (hcIter!=hcEnd){
01419         if (hcIter->plane==pl){  //if it's in current plane:
01420           winpl[d] = hcIter->plane;
01421           winpl0[d] = hcIter->upper;
01422           winpl1[d] = hcIter->lower;
01423           winpl0_tpos[d] = hcIter->upper_tpos;
01424           winpl1_tpos[d] = hcIter->lower_tpos;
01425           winpl2[d] = hcIter->ph;
01426           winpl3[d] = hcIter->type;
01427                 
01428           if (shwmid==0) {   // shower max plane
01429             refp = winpl[d];  // set reference point
01430             refs = (winpl0_tpos[d]+winpl1_tpos[d])/2.;
01431           }
01432           //determine widest window:
01433           if(shwwid<(winpl0_tpos[d]-winpl1_tpos[d]+STRIPWIDTHINMETRES)) 
01434             shwwid = winpl0_tpos[d]-winpl1_tpos[d]+STRIPWIDTHINMETRES; 
01435           // average centroid of event so far:
01436           shwmid = (shwmid*clusterph+(winpl0_tpos[d]+winpl1_tpos[d]) * 
01437                     winpl2[d]/2.)/(clusterph+winpl2[d]); 
01438           // calculate average distance to centroid:
01439           clusterph += winpl2[d];
01440           winpl4[d] = shwwid;
01441           winpl5[d] = shwmid;
01442           if (maxwid<winpl4[d]) maxwid = winpl4[d];
01443           if (maxwid<2.*STRIPWIDTHINMETRES) maxwid = 2.*STRIPWIDTHINMETRES;
01444         }
01445         hcIter++;
01446       }
01447     }
01448     std::vector<Double_t> MCV = BestHough(houghcluster,false);
01449     if(MCV[5]!=0){
01450       slopen = MCV[5];
01451       eslopen = MCV[8];
01452       trustslopen = true;
01453     }
01454     else if(MCV[0]!=0){
01455       slopen = MCV[0];
01456       eslopen = MCV[3];
01457       trustslopen = true;
01458     }
01459     else {
01460       slopen = 0;
01461       eslopen = 0;
01462       trustslopen = false;
01463     }
01464   }
01465   
01467   //Try to fill planes with no windows
01468   
01469   //First find centroid and width for +/- 1 planes
01470   
01471   for (Int_t d = 0;d<numLongPlanes;d++){
01472     int pl = -1;
01473     if(numPreMaxPlanes>numPostMaxPlanes){
01474       if (d<numPreMaxPlanes) pl = maxPlane-d;
01475       else pl = maxPlane+d-numPreMaxPlanes+1;
01476     }
01477     else{
01478       if (d<numPostMaxPlanes) pl = maxPlane+d;
01479       else pl = maxPlane-d+numPostMaxPlanes-1;
01480     }
01481     
01482     double minwei = 200.;   //minimum distance to centroid (weight)
01483     double minwei_tpos = 10.;   //minimum distance to centroid (weight)
01484     int newwinpl0 = -1;     //new window upper and lower strip bounds
01485     int newwinpl1 = -1;
01486     TSpecReal_t newwinpl0_tpos = -10; //new window upp and low strip tpos bounds
01487     TSpecReal_t newwinpl1_tpos = -10;
01488     double newwinpl2 = 0.;
01489     int newwinpl3 = 0;
01490 
01491     if(winpl0[d]==-1&&winpl1[d]==-1){  //does this plane have a solution?
01492       
01493       int p = 0;
01494       int m = 0;
01495       if (det==0||(det==1&&(pl-249)*(pl+2-249)>0)) p = 2;
01496       else p = 3;
01497       if (det==0||(det==1&&(pl-249)*(pl-2-249)>0)) m = 2;
01498       else m = 3;
01499       
01500       double pmid = -1;   //centroid for next plane (plus)
01501       double mmid = -1;   //centroid for prev plane (minus)
01502       double pwid = -1;   //width
01503       double mwid = -1;   
01504       int pwt = -1;      // window type  plus
01505       int mwt = -1;     //               minus
01506       bool useslopep = false; 
01507       bool useslopem = false;
01508       
01509       for (Int_t c = 0;c<numLongPlanes;c++){
01510         if (winpl[c]==pl+p) {
01511           pmid = winpl5[c];
01512           if (trustslopen) useslopep = true;
01513           pwid = winpl4[c];
01514           pwt = winpl3[c];
01515         }
01516         if (winpl[c]==pl-m) {
01517           mmid = winpl5[c];      
01518           if (trustslopen) useslopem = true;
01519           mwid = winpl4[c];
01520           mwt = winpl3[c];
01521         }
01522       }
01523       
01524       //Find compatible windows
01525       
01526       double usewid = maxwid;
01527       if(maxwid>0){
01528         if((mmid==-1&&pmid==-1)||
01529            (!useslopem&&mmid>=0&&TMath::Abs(shwmid-mmid)>maxwid)||
01530            (!useslopep&&pmid>=0&&TMath::Abs(shwmid-pmid)>maxwid)||
01531            (useslopem&&mmid>=0&&
01532               (TMath::Abs(shwmid-mmid-(maxPlane-pl+m)*slopen)>maxwid)||
01533             (TMath::Abs(shwmid-mmid-(maxPlane-pl+m)*slopen) > 
01534              TMath::Abs(maxPlane-pl+m)*eslopen))||
01535            (useslopep&&pmid>=0&&
01536               (TMath::Abs(shwmid-pmid-(maxPlane-pl-p)*slopen)>maxwid)||
01537             (TMath::Abs(shwmid-pmid-(maxPlane-pl-p)*slopen) > 
01538              TMath::Abs(maxPlane-pl-p)*eslopen))){
01539           minwei = 200.;
01540           minwei_tpos = 10.;
01541           newwinpl0 = -1;
01542           newwinpl1 = -1;
01543           newwinpl0_tpos = -10;
01544           newwinpl1_tpos = -10;
01545           newwinpl2 = 0.;
01546           winIter = win.begin();
01547           while(winIter!=winEnd){
01548             if(winIter->plane==pl){
01549               if (//shwmid>=0&&
01550                   ((winIter->upper_tpos-shwmid) * 
01551                    (winIter->lower_tpos-shwmid)<=
01552                    (maxwid/2.)*(maxwid/2.)) && 
01553                   TMath::Abs((winIter->upper_tpos + 
01554                               winIter->lower_tpos)/2.-shwmid)<minwei){
01555                 minwei = TMath::Abs((winIter->upper_tpos + 
01556                                      winIter->lower_tpos)/2.-shwmid);
01557                 newwinpl0 = winIter->upper;
01558                 newwinpl1 = winIter->lower;
01559                 newwinpl0_tpos = winIter->upper_tpos;
01560                 newwinpl1_tpos = winIter->lower_tpos;
01561                 newwinpl2 = winIter->ph;
01562                 newwinpl3 = winIter->type;
01563               }
01564             }
01565             winIter++;
01566           }
01567         }
01568         else{
01569           minwei = 200.;
01570           minwei_tpos = 10.;
01571           newwinpl0 = -1;
01572           newwinpl1 = -1;
01573           newwinpl0_tpos = -10.;
01574           newwinpl1_tpos = -10.;
01575           newwinpl2 = 0.;
01576           winIter = win.begin();
01577           while(winIter!=winEnd){
01578             if(winIter->plane==pl){
01579               if (!useslopem&&mmid>=0
01580                   &&((winIter->upper_tpos-mmid) * 
01581                      (winIter->lower_tpos-mmid) <=
01582                      (maxwid/2.)*(maxwid/2.)) &&
01583                   TMath::Abs((winIter->upper_tpos + 
01584                               winIter->lower_tpos)/2.-mmid)<minwei){
01585                 minwei = TMath::Abs((winIter->upper_tpos + 
01586                                      winIter->lower_tpos)/2.-mmid);
01587                 newwinpl0 = winIter->upper;
01588                 newwinpl1 = winIter->lower;
01589                 newwinpl0_tpos = winIter->upper_tpos;
01590                 newwinpl1_tpos = winIter->lower_tpos;
01591                 newwinpl2 = winIter->ph;
01592                 newwinpl3 = winIter->type;
01593               }
01594               else if (!useslopep&&pmid>=0
01595                        &&((winIter->upper_tpos-pmid) * 
01596                           (winIter->lower_tpos-pmid) <= 
01597                           (maxwid/2.)*(maxwid/2.))
01598                        && TMath::Abs((winIter->upper_tpos + 
01599                                       winIter->lower_tpos)/2.-pmid)<minwei){
01600                 minwei = TMath::Abs((winIter->upper_tpos + 
01601                                      winIter->lower_tpos)/2.-pmid);
01602                 newwinpl0 = winIter->upper;
01603                 newwinpl1 = winIter->lower;
01604                 newwinpl0_tpos = winIter->upper_tpos;
01605                 newwinpl1_tpos = winIter->lower_tpos;
01606                 newwinpl2 = winIter->ph;
01607                 newwinpl3 = winIter->type;
01608               }
01609               else if (useslopem&&mmid>=0){
01610                 if(usewid<m*eslopen) usewid = m*eslopen;
01611                 else usewid = maxwid;
01612                 if((winIter->upper_tpos-mmid-m*slopen) * 
01613                    (winIter->lower_tpos-mmid-m*slopen) <= 
01614                    (usewid/2.)*(usewid/2.)
01615                    &&TMath::Abs((winIter->upper_tpos + 
01616                                  winIter->lower_tpos)/2.-mmid-m*slopen)<minwei){
01617                   minwei = TMath::Abs((winIter->upper_tpos + 
01618                                        winIter->lower_tpos)/2.-mmid-m*slopen);
01619                   newwinpl0 = winIter->upper;
01620                   newwinpl1 = winIter->lower;
01621                   newwinpl0_tpos = winIter->upper_tpos;
01622                   newwinpl1_tpos = winIter->lower_tpos;
01623                   newwinpl2 = winIter->ph;
01624                   newwinpl3 = winIter->type;
01625                 }
01626               }
01627               else if (useslopep&&pmid>=0){
01628                 if(usewid<p*eslopen) usewid = p*eslopen;
01629                 else usewid = maxwid; 
01630                 if((winIter->upper_tpos-pmid+p*slopen) * 
01631                    (winIter->lower_tpos-pmid+p*slopen) <= 
01632                    (usewid/2.)*(usewid/2.)
01633                    &&TMath::Abs((winIter->upper_tpos + 
01634                                  winIter->lower_tpos)/2.-pmid+p*slopen)<minwei){
01635                   minwei = TMath::Abs((winIter->upper_tpos + 
01636                                        winIter->lower_tpos)/2.-pmid+p*slopen);
01637                   newwinpl0 = winIter->upper;
01638                   newwinpl1 = winIter->lower;
01639                   newwinpl0_tpos = winIter->upper_tpos;
01640                   newwinpl1_tpos = winIter->lower_tpos;
01641                   newwinpl2 = winIter->ph;
01642                   newwinpl3 = winIter->type;
01643                 }
01644               }
01645             }
01646             winIter++;
01647           }
01648         }
01649       }
01650       if(newwinpl0>=0) {
01651         winpl0[d] = newwinpl0;
01652         winpl1[d] = newwinpl1;
01653         winpl0_tpos[d] = newwinpl0_tpos;
01654         winpl1_tpos[d] = newwinpl1_tpos;
01655         winpl2[d] = newwinpl2;
01656         winpl3[d] = newwinpl3;
01657         shwmid = (shwmid*clusterph + 
01658                   (winpl0_tpos[d]+winpl1_tpos[d]) * 
01659                   winpl2[d]/2.)/(clusterph+winpl2[d]);
01660         clusterph+=winpl2[d];
01661         MSG("SubShowerSR",Msg::kDebug) <<"find missed windows "<<pl
01662                                        <<" "<<newwinpl0<<" "
01663                                        <<newwinpl1<<" "<<minwei<<endl;
01664       }
01665     }
01666     
01667     window winSolution;
01668     winSolution.plane = pl;
01669     winSolution.upper = winpl0[d];
01670     winSolution.lower = winpl1[d];
01671     winSolution.upper_tpos = winpl0_tpos[d];
01672     winSolution.lower_tpos = winpl1_tpos[d];
01673     winSolution.ph = winpl2[d];
01674     winSolution.type = winpl3[d];
01675     winSolution.centroid = 0.;
01676     winSolution.phpos = 0.;
01677     winSolution.phwidth = 0.;
01678     cluster.push_back(winSolution);
01679   }
01680 
01681   return cluster;
01682 }


Member Data Documentation

Definition at line 66 of file AlgSubShowerSRList.h.

Referenced by BestHough(), and RunAlg().

Definition at line 68 of file AlgSubShowerSRList.h.

Referenced by CleanUp(), and RunAlg().

Definition at line 69 of file AlgSubShowerSRList.h.

Referenced by CleanUp(), and RunAlg().

Definition at line 67 of file AlgSubShowerSRList.h.

Referenced by CleanUp(), RunAlg(), and TimeTest().

Definition at line 70 of file AlgSubShowerSRList.h.

Referenced by RunAlg().

Definition at line 71 of file AlgSubShowerSRList.h.

Referenced by FormHalo(), and RunAlg().

Definition at line 72 of file AlgSubShowerSRList.h.

Referenced by FormHalo(), and RunAlg().

Definition at line 63 of file AlgSubShowerSRList.h.

Referenced by FindCluster(), and RunAlg().

Definition at line 62 of file AlgSubShowerSRList.h.

Referenced by FindCluster(), and RunAlg().

Definition at line 74 of file AlgSubShowerSRList.h.

Referenced by AreNeighbours(), and RunAlg().

Double_t AlgSubShowerSRList::fMinHoughPH [private]

Definition at line 65 of file AlgSubShowerSRList.h.

Referenced by HoughTransCluster(), and RunAlg().

Definition at line 64 of file AlgSubShowerSRList.h.

Referenced by HoughTransCluster(), and RunAlg().

Double_t AlgSubShowerSRList::fMinStripPE [private]

Definition at line 73 of file AlgSubShowerSRList.h.

Referenced by FormHalo(), and RunAlg().

Definition at line 75 of file AlgSubShowerSRList.h.

Referenced by FormHalo(), and RunAlg().

Definition at line 78 of file AlgSubShowerSRList.h.

Referenced by BestHough(), and ~AlgSubShowerSRList().

Definition at line 79 of file AlgSubShowerSRList.h.

Referenced by BestHough(), and ~AlgSubShowerSRList().

Definition at line 80 of file AlgSubShowerSRList.h.

Referenced by TransCluster(), and ~AlgSubShowerSRList().

TH1F* AlgSubShowerSRList::vtxHist [private]

Definition at line 77 of file AlgSubShowerSRList.h.

Referenced by FormHalo(), and ~AlgSubShowerSRList().


The documentation for this class was generated from the following files:

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1