#include <AlgSliceSRList.h>
Inheritance diagram for AlgSliceSRList:
Public Member Functions | |
| AlgSliceSRList () | |
| virtual | ~AlgSliceSRList () |
| virtual void | RunAlg (AlgConfig &ac, CandHandle &ch, CandContext &cx) |
| virtual void | SlicetheSnarl (AlgConfig &ac, CandHandle &ch, CandContext &cx) |
| virtual void | SlicetheSnarl_ASAP (AlgConfig &ac, CandHandle &ch, CandContext &cx) |
| virtual void | SlicetheSnarl_MST (AlgConfig &ac, CandHandle &ch, CandContext &cx) |
| virtual void | PassAll (AlgConfig &ac, CandHandle &ch, CandContext &cx) |
| virtual void | Trace (const char *c) const |
|
|
Definition at line 58 of file AlgSliceSRList.cxx. 00059 {
00060 }
|
|
|
Definition at line 63 of file AlgSliceSRList.cxx. 00064 {
00065 }
|
|
||||||||||||||||
|
pass all strips in the snarl to a single slice Definition at line 826 of file AlgSliceSRList.cxx. References CandHandle::AddDaughterLink(), AlgFactory::GetAlgHandle(), CandContext::GetDataIn(), CandHandle::GetDaughterIterator(), AlgFactory::GetInstance(), CandContext::GetMom(), and CandSlice::MakeCandidate(). Referenced by RunAlg(). 00827 {
00828
00829 const CandStripListHandle *cslh =
00830 dynamic_cast<const CandStripListHandle*>(cx.GetDataIn());
00831
00832 AlgFactory &af = AlgFactory::GetInstance();
00833 AlgHandle ah = af.GetAlgHandle("AlgSliceSR","default");
00834 CandContext cxx(this,cx.GetMom());
00835
00836 CandStripHandleItr cshItr(cslh->GetDaughterIterator());
00837
00838 // Iterate over CandStripList daughters
00839 TObjArray striparray;
00840 striparray.Clear();
00841 while (CandStripHandle *curr = cshItr()) {
00842 striparray.Add(curr);
00843 }
00844 cxx.SetDataIn(&striparray);
00845 CandSliceHandle slicehandle = CandSlice::MakeCandidate(ah,cxx);
00846 ch.AddDaughterLink(slicehandle);
00847 }
|
|
||||||||||||||||
|
Main driver routine for CandSliceList algorithm. Based on context, call is made to near or far det slicer. For far detector, value of SliceFar is used to determine whether slicing is performed, or entire snarl is passed. Default currently is to bypass slicing in the far detector. Implements AlgBase. Definition at line 71 of file AlgSliceSRList.cxx. References RawRecord::FindRawBlock(), CandContext::GetDataIn(), MomNavigator::GetFragment(), Registry::GetInt(), CandContext::GetMom(), MSG, PassAll(), SlicetheSnarl(), SlicetheSnarl_ASAP(), and SlicetheSnarl_MST(). 00072 {
00073 assert(cx.GetDataIn());
00074
00075 Int_t passall = ac.GetInt("PassAll");
00076 Int_t passasap = ac.GetInt("PassASAP");
00077 Int_t passmst = ac.GetInt("PassMST");
00078
00079 // Check for CandStripListHandle input
00080 if (cx.GetDataIn()->InheritsFrom("CandStripListHandle")) {
00081 const MomNavigator *mom = cx.GetMom();
00082 RawRecord *rr = dynamic_cast<RawRecord *>(mom->GetFragment("RawRecord"));
00083 if (!rr) {
00084 MSG("SliceSR", Msg::kWarning) << "No RawRecord in MOM." << endl;
00085 return;
00086 }
00087 const RawDigitDataBlock *rddb = dynamic_cast<const RawDigitDataBlock *>
00088 (rr->FindRawBlock("RawDigitDataBlock"));
00089 if (!rddb) {
00090 MSG("SliceSR", Msg::kWarning) << "No RawDigitDataBlock in RawRecord." << endl;
00091 return;
00092 }
00093
00094 if(passall){
00095 PassAll(ac,ch,cx);
00096 }
00097 if(passasap){
00098 SlicetheSnarl_ASAP(ac,ch,cx);
00099 }
00100 if(passmst){
00101 SlicetheSnarl_MST(ac,ch,cx);
00102 }
00103 else if(!passall && !passasap && !passmst){
00104 SlicetheSnarl(ac,ch,cx);
00105 }
00106 }
00107 }
|
|
||||||||||||||||
|
Candstrips in the upstream (non-spectrometer) planes are treated first. For these planes the Candstrips are sorted by ToF-corrected time, keeping only strips which return GetCharge()>MinCharge. CandStrips meeting this requirement are placed in a list, and the time of the new CandStrip is checked against that of the previous last entry. If this time difference is greater than MaxTimeGap, and the number of CandStrips is greater than MinStrip, the CandStrips in the list are assumed to constitute a CandSlice, and a CandSlice is constructed from these members. If this time difference is less than MaxTimeGap, but the total list time duration is greater than TimeWindow, we assume that we have missed the true start of this CandSlice. In this case, we first check to see if, by removing the first entry on the list, and adding the new entry, we end up with a shorter time duration in the list. If so, we perform this operation, adding the new CandStrip, and dropping the earliest from the list. In the next phase, we pick up CandStrips in the spectrometer. To do this, we iterate over CandSlices constructed in the previous stage, and determine which remaining CandStrips should be added. For each CandSlice, a time interval bounded by the CandSlice start time minus EarlyTimeDiff and the CandSlice start time plus TimeDiffSpect is defined, and a list of all spectrometer CandStrips with times inside this interval is constructed. This list of CandStrips is now looped over. If the CandStrip is in the spectrometer, then this CandStrip is added to the current CandSlice. In the upstream planes, a more stringent requirement is placed on the addition of the CandStrip. First, the time difference between the CandStrip and CandSlice start must be less than TimeDiff. If this is the case, the location of the CandStrip is checked against the locations of the current CandSlice members. There must be at least one existing CandSlice member with a plane within PlaneDiff of the CandStrip, and a transverse position within TposDiff of the CandStrip. If all conditions are met, the CandStrip is added to the CandSlice. Next, the CandSlice list is iterated over again, this time to pick up the CandStrips with charge < MinCharge. For each CandSlice, a time interval bounded by the CandSlice start time and the CandSlice start time plus TimeWindow is defined, and a list of all CandStrips with times inside this interval is constructed. If the pulse height is less than MinCharge and the CandStrip is not in the spectrometer, the CandStrip is added to the CandSlice. Definition at line 115 of file AlgSliceSRList.cxx. References CandHandle::AddDaughterLink(), Registry::Clear(), AlgFactory::GetAlgHandle(), CandSliceHandle::GetCharge(), CandStripHandle::GetCharge(), CandSliceHandle::GetCorrBegTime(), CandStripHandle::GetCorrBegTime(), CandSliceHandle::GetCorrTime(), CandContext::GetDataIn(), CandHandle::GetDaughterIterator(), VldContext::GetDetector(), Registry::GetDouble(), AlgFactory::GetInstance(), Registry::GetInt(), CandContext::GetMom(), CandStripHandle::GetPlane(), CandStripHandle::GetPlaneView(), CandStripHandle::GetStrip(), CandStripHandle::GetTime(), CandStripHandle::GetTPos(), CandHandle::GetVldContext(), CandStripHandle::GetZPos(), CandSlice::MakeCandidate(), MSG, CandHandle::RemoveDaughter(), SliceSRKeyFromForwardTime(), and StripSRKeyFromCorrTime(). Referenced by RunAlg(). 00116 {
00117
00118 const CandStripListHandle *cslh =
00119 dynamic_cast<const CandStripListHandle*>(cx.GetDataIn());
00120
00121 MSG("AlgSliceSR",Msg::kDebug) << "Starting SR Slicer!" << endl;
00122
00123 // Sort by Plane and Strip
00124 CandStripHandleItr cshItr(cslh->GetDaughterIterator());
00125 CandStripHandleKeyFunc *cshKf = cshItr.CreateKeyFunc();
00126 cshKf->SetFun(StripSRKeyFromCorrTime);
00127 cshItr.GetSet()->AdoptSortKeyFunc(cshKf);
00128 cshKf = 0;
00129
00130 AlgFactory &af = AlgFactory::GetInstance();
00131 AlgHandle ah = af.GetAlgHandle("AlgSliceSR","default");
00132 CandContext cxx(this,cx.GetMom());
00133
00134 Double_t timewindow = ac.GetDouble("TimeWindow");
00135 Int_t minstrip = ac.GetInt("MinStrip");
00136 Double_t mincharge = ac.GetDouble("MinCharge");
00137 Double_t maxtimegap = ac.GetDouble("MaxTimeGap");
00138 Double_t earlytimediff = ac.GetDouble("EarlyTimeDiff");
00139 Double_t timediffspect = ac.GetDouble("TimeDiffSpect");
00140 Double_t minpulseheight = ac.GetDouble("MinPulseHeight");
00141 Double_t zsplitdistance = ac.GetDouble("ZSplitDistance");
00142 Double_t trackiness = ac.GetDouble("Trackiness");
00143 Double_t mincombiningdistance = ac.GetDouble("MinCombiningDistance");
00144 Int_t combiningonoff = ac.GetInt("CombiningOnOff");
00145 earlytimediff = fabs(earlytimediff); // make sure this is a positive quantity
00146
00147 Int_t specplane = 999;
00148 if (cslh->GetVldContext()->GetDetector() == DetectorType::kNear) specplane = 121;
00149
00150 // Iterate over CandStripList daughters
00151 TObjArray striparray;
00152 striparray.Clear();
00153 Int_t ifirst = 0;
00154 Int_t nstrip = 0;
00155 Double_t dtime = 0.;
00156 while (CandStripHandle *curr = cshItr()) {
00157
00158 MSG("AlgSliceSR",Msg::kDebug) << " strip plane " << curr->GetPlane() << " strip " << curr->GetStrip() << " charge " << curr->GetCharge() << endl;
00159
00160 if (curr->GetPlane()<specplane && curr->GetCharge()>=mincharge) {
00161 // Only consider upstream of muon spectrometer region and strips of a minimum pulse height
00162 if (!striparray.First()) {
00163 MSG("AlgSliceSR",Msg::kDebug) << " first strip " << endl;
00164 striparray.Add(curr);
00165 nstrip++;
00166 } else if (curr->GetCorrBegTime()-dynamic_cast<CandStripHandle*> (striparray.Last())->GetCorrBegTime()>maxtimegap) {
00167 // If difference between this TOF corrected time and last TOF corrected
00168 // time in striparray is greater than maxtimegap, form a new candslice
00169 if (striparray.GetLast()+1>=minstrip) {
00170 MSG("AlgSliceSR",Msg::kDebug) << " starting new slice " << endl;
00171 cxx.SetDataIn(&striparray);
00172 CandSliceHandle slicehandle =
00173 CandSlice::MakeCandidate(ah,cxx);
00174 ch.AddDaughterLink(slicehandle);
00175 }
00176 striparray.Clear();
00177 striparray.Add(curr);
00178 ifirst = 0;
00179 nstrip = 1;
00180 } else {
00181 if (curr->GetCorrBegTime()<=dynamic_cast<CandStripHandle*>
00182 (striparray.At(ifirst))->GetCorrBegTime()+timewindow) {
00183 // If difference between this TOF corrected time and first TOF corrected
00184 // time in striparray is less than timewindow, add to striparray
00185 MSG("AlgSliceSR",Msg::kDebug) << " adding strip " << endl;
00186 striparray.Add(curr);
00187 nstrip++;
00188 dtime = curr->GetCorrBegTime()-dynamic_cast<CandStripHandle*>
00189 (striparray.At(0))->GetCorrBegTime();
00190 // dtime is the time extend of striparray
00191 } else {
00192 /* Gap is less than MaxTimeGap, but the total list time duration is greater than TimeWindow. In this case we assume that we have missed the true start of this CandSlice. We now check to see if, by removing the first entry on the list, and adding the new entry, we end up with a shorter time duration in the list.
00193 */
00194 if (striparray.At(ifirst+1) &&
00195 curr->GetCorrBegTime()-dynamic_cast<CandStripHandle*>
00196 (striparray.At(ifirst+1))->GetCorrBegTime()<dtime &&
00197 striparray.GetLast()+1-ifirst>=nstrip) {
00198 MSG("AlgSliceSR",Msg::kDebug) << " removing first/adding " << endl;
00199 // Removing first strip and adding last results in a shorter slice, so assume this is the thing to do
00200 striparray.RemoveAt(0);
00201 striparray.Compress();
00202 striparray.Add(curr);
00203 nstrip = striparray.GetLast()-ifirst;
00204 dtime = curr->GetCorrBegTime()-dynamic_cast<CandStripHandle*> (striparray.At(0))->GetCorrBegTime();
00205 } else {
00206 // Otherwise, make a new slice
00207 if (striparray.GetLast()+1>=minstrip) {
00208 MSG("AlgSliceSR",Msg::kDebug) << " make new slice " << endl;
00209 cxx.SetDataIn(&striparray);
00210 CandSliceHandle slicehandle =
00211 CandSlice::MakeCandidate(ah,cxx);
00212 ch.AddDaughterLink(slicehandle);
00213 }
00214 // Prepare for next slice
00215 striparray.Clear();
00216 striparray.Add(curr);
00217 ifirst = 0;
00218 nstrip = 1;
00219 }
00220 }
00221 }
00222 }
00223 }
00224
00225 // Cleanup - make slice out of leftovers
00226 if (striparray.GetLast()+1>=minstrip) {
00227 MSG("AlgSliceSR",Msg::kDebug) << " final slice make " << endl;
00228 cxx.SetDataIn(&striparray);
00229 CandSliceHandle slicehandle = CandSlice::MakeCandidate(ah,cxx);
00230 ch.AddDaughterLink(slicehandle);
00231 }
00232
00233 /* Fix eaten slices, doing it with Z first (only NEAR DET). Start off with 1st strip hit and find everything within a distance of it, for all hits, put left overs in a separate array, now have 2 slices. In the slice that is further in Z (presumably a track), one checks how many planes have > 2 hits, if it truly is a track then it should have about 1 hit / plane. Counting the number of planes that have at least a hit and those that have > 2 hits, I take the ratio, if it is less than 0.05 split, otherwise do not.
00234 */
00235
00236 if (cslh->GetVldContext()->GetDetector() == DetectorType::kNear){
00237 MSG("AlgSliceSR",Msg::kDebug) << "Z Splitting!" << endl;
00238 CandSliceHandleItr sliceItr(ch.GetDaughterIterator());
00239 CandSliceHandleKeyFunc *sliceKf = sliceItr.CreateKeyFunc();
00240 sliceKf->SetFun(SliceSRKeyFromForwardTime);
00241 sliceItr.GetSet()->AdoptSortKeyFunc(sliceKf);
00242 sliceKf = 0;
00243
00244 sliceItr.ResetLast();
00245 TObjArray z_split;
00246 z_split.Clear();
00247 TObjArray z_split2;
00248 z_split2.Clear();
00249 while (CandSliceHandle *curr = sliceItr()) {
00250 z_split.Clear(); z_split2.Clear();
00251 CandStripHandleItr currstripItr(curr->GetDaughterIterator());
00252 while(CandStripHandle *strip3 = currstripItr()) {
00253 z_split.Add(strip3);
00254 }
00255 CandStripHandle *temp = dynamic_cast<CandStripHandle*>(z_split.At(0));
00256 z_split2.Add(temp); z_split.RemoveAt(0); z_split.Compress();
00257 while(z_split.GetEntries() > 0) {
00258 Int_t sa_size = z_split.GetEntries();
00259 for(int i=0; i<=z_split2.GetLast(); i++) {
00260 CandStripHandle *sa2 = dynamic_cast<CandStripHandle*>(z_split2.At(i));
00261 for(int j=0; j<=z_split.GetLast(); j++) {
00262 CandStripHandle *sa = dynamic_cast<CandStripHandle*>(z_split.At(j));
00263 if(fabs(sa->GetZPos()-sa2->GetZPos()) < zsplitdistance) {
00264 z_split2.Add(sa);
00265 z_split.RemoveAt(j);
00266 z_split.Compress();
00267 }
00268 }
00269 }
00270 Int_t sa_size2 = z_split.GetEntries();
00271 if(sa_size2 == sa_size) {
00272 break;
00273 }
00274 }
00275
00276 //Now need to check if I'm not splitting a shower and a track that went through the coil!
00277 int planer[300] = {0}; int planer2[300] = {0};
00278 int bigplane = 0, bigplane2 = 0;
00279 for (Int_t i=0; i<=z_split.GetLast(); i++) {
00280 CandStripHandle *csh = dynamic_cast<CandStripHandle*>(z_split.At(i));
00281 planer[csh->GetPlane()]++;
00282 if(csh->GetPlane() > bigplane) bigplane = csh->GetPlane();
00283 }
00284 for (Int_t i=0; i<=z_split2.GetLast(); i++) {
00285 CandStripHandle *csh = dynamic_cast<CandStripHandle*>(z_split2.At(i));
00286 planer2[csh->GetPlane()]++;
00287 if(csh->GetPlane() > bigplane2) bigplane2 = csh->GetPlane();
00288 }
00289
00290 int isbad = 0, nump = 0;
00291 if(bigplane > bigplane2) {
00292 for(int i=0; i<300; i++) {
00293 if(planer[i]>0) nump++;
00294 if(planer[i]>2) isbad++;
00295 }
00296 if(nump == 0) {
00297 nump = 100; isbad = 1;
00298 }
00299 } else {
00300 for(int i=0; i<300; i++) {
00301 if(planer2[i]>0) nump++;
00302 if(planer2[i]>2) isbad++;
00303 }
00304 if(nump == 0) {
00305 nump = 100; isbad = 1;
00306 }
00307 }
00308
00309 if((float) isbad / nump > trackiness) {
00310 if(z_split2.GetEntries() >= minstrip) {
00311 striparray.Clear();
00312 for(int i=0; i<=z_split2.GetLast(); i++) {
00313 CandStripHandle *temp = dynamic_cast<CandStripHandle*>(z_split2.At(i));
00314 striparray.Add(temp);
00315 }
00316 cxx.SetDataIn(&striparray);
00317 CandSliceHandle slicehandle = CandSlice::MakeCandidate(ah,cxx);
00318 ch.AddDaughterLink(slicehandle);
00319 }
00320 if(z_split.GetEntries() >= minstrip) {
00321 striparray.Clear();
00322 for(int i=0; i<=z_split.GetLast(); i++) {
00323 CandStripHandle *temp = dynamic_cast<CandStripHandle*>(z_split.At(i));
00324 striparray.Add(temp);
00325 }
00326 cxx.SetDataIn(&striparray);
00327 CandSliceHandle slicehandle = CandSlice::MakeCandidate(ah,cxx);
00328 ch.AddDaughterLink(slicehandle);
00329 }
00330 ch.RemoveDaughter(curr);
00331 }
00332 z_split.Clear(); z_split2.Clear(); striparray.Clear();
00333 }
00334 }
00335
00336 // Prepare slice iterator for next phase
00337 CandSliceHandleItr sliceItr2(ch.GetDaughterIterator());
00338 CandSliceHandleKeyFunc *sliceKf2 = sliceItr2.CreateKeyFunc();
00339 sliceKf2->SetFun(SliceSRKeyFromForwardTime);
00340 sliceItr2.GetSet()->AdoptSortKeyFunc(sliceKf2);
00341 sliceKf2 = 0;
00342
00343 /* We now pick up CandStrips in the spectrometer. To do this, we iterate over CandSlices constructed in the previous stage, and determine which remaining CandStrips should be added. For each CandSlice, a time interval bounded by the CandSlice start time minus EarlyTimeDiff and the CandSlice start time plus TimeDiffSpect is defined, and CandStrips in the spectrometer within this time range are added to the slice */
00344
00345 if (cslh->GetVldContext()->GetDetector() == DetectorType::kNear){
00346 MSG("AlgSliceSR",Msg::kDebug) << "Adding Spectro Hits!" << endl;
00347 sliceItr2.Reset();
00348 while (CandSliceHandle *slice = sliceItr2()) {
00349 cshItr.GetSet()->Slice();
00350 cshItr.GetSet()->
00351 Slice(slice->GetCorrTime()-earlytimediff,
00352 slice->GetCorrTime()+timediffspect);
00353 while (CandStripHandle *strip1 = cshItr()) {
00354 if (strip1->GetPlane()>= specplane){
00355 slice->AddDaughterLink(*strip1);
00356 }
00357 }
00358 }
00359 }
00360
00361 /* we next pick up low pulse height strips. To do this, we iterate over CandSlices constructed in the previous stage, and determine which remaining CandStrips should be added. For each CandSlice, a time interval bounded by the CandSlice start time and the start time of the NEXT CandSlice is defined, and CandStrips with q<mincharge within this time range are added to the slice */
00362
00363 MSG("AlgSliceSR",Msg::kDebug) << "Adding Small Hits" << endl;
00364 Double_t SliceTime[100] = {0.0};
00365 Int_t counter = 0;
00366 sliceItr2.Reset();
00367 while (CandSliceHandle *slice = sliceItr2()) {
00368 SliceTime[counter] = slice->GetCorrBegTime()*1e9;
00369 counter++;
00370 }
00371 SliceTime[counter] = 30000.0;
00372 counter = 0;
00373 sliceItr2.Reset();
00374 while (CandSliceHandle *slice = sliceItr2()) {
00375 MSG("AlgSliceSR",Msg::kDebug) << " ## slice ## " << counter << " time " << SliceTime[counter] << " next time " << SliceTime[counter+1] << endl;
00376 Float_t ticker = fabs(SliceTime[counter+1]/1e9-SliceTime[counter]/1e9);
00377 cshItr.GetSet()->Slice();
00378 cshItr.GetSet()->Slice(slice->GetCorrBegTime()-(earlytimediff/2.0), slice->GetCorrBegTime()+ticker-(earlytimediff/2.0));
00379 while (CandStripHandle *strip = cshItr()) {
00380 MSG("AlgSliceSR",Msg::kDebug) << " strip plane " << strip->GetPlane() << " strip " << strip->GetStrip() << " charge " << strip->GetCharge() << " time " << strip->GetCorrBegTime()*1e9 << endl;
00381 if (strip->GetPlane()<specplane && strip->GetCharge() < mincharge) {
00382 slice->AddDaughterLink(*strip);
00383 MSG("AlgSliceSR",Msg::kDebug) << " adding " << endl;
00384 }
00385 }
00386 counter++;
00387 }
00388
00389 /* Here we get rid of slices with < 2000 ADC counts (mostly junk anyway) */
00390
00391 MSG("AlgSliceSR",Msg::kDebug) << "Killing Small Slices" << endl;
00392 sliceItr2.Reset();
00393 while( CandSliceHandle *slice = sliceItr2()) {
00394 if(slice->GetCharge(CalDigitType::kNone) < minpulseheight) {
00395 ch.RemoveDaughter(slice);
00396 }
00397 }
00398
00399 // Prepare Slice iterator for next phase.
00400 CandSliceHandleItr sliceItr3(ch.GetDaughterIterator());
00401 CandSliceHandleKeyFunc *sliceKf3 = sliceItr3.CreateKeyFunc();
00402 sliceKf3->SetFun(SliceSRKeyFromForwardTime);
00403 sliceItr3.GetSet()->AdoptSortKeyFunc(sliceKf3);
00404 sliceKf3 = 0;
00405
00406 /* Combining Step. Calculate the Center of Pulse height for all slices < 10,000 ADC counts and loop over all slices that occur before said slice in time and have > 10,000 ADC counts. Find the min distance between the small slice and the bigger slice. If the dist in either view is less than a parameter (mincombiningdistance) then combine the 2 slices. This significantly reduces the number of low-completeness slices, while not affecting the high-completeness regime negatively. One can see the improvement on slices, but CandEventSR is at a stage where this doesn't help 'yet'. It can be turned on or off. 12-15-04 */
00407
00408 if(combiningonoff) {
00409 MSG("AlgSliceSR",Msg::kDebug) << "Starting Combining!" << endl;
00410 Int_t badindex = 0;
00411 TObjArray stripu; stripu.Clear();
00412 TObjArray stripv; stripv.Clear();
00413 while(CandSliceHandle *bad = sliceItr3()) {
00414 stripu.Clear(); stripv.Clear();
00415 if(bad->GetCharge(CalDigitType::kNone) < 10000) {
00416 CandStripHandleItr stripitr(bad->GetDaughterIterator());
00417 Float_t BegTimeBad = 30000.0;
00418 while(CandStripHandle *curr = stripitr()) {
00419 if(curr->GetTime()*1e9 < BegTimeBad) BegTimeBad = curr->GetTime()*1e9;
00420 if(curr->GetPlane()<121 && curr->GetCharge()>2.0) {
00421 if(curr->GetPlaneView() == 2) {
00422 stripu.Add(curr);
00423 } else {
00424 stripv.Add(curr);
00425 }
00426 }
00427 }
00428 //Calculating Center of Pulse Height Coordinates
00429 Float_t phtotuz = 0.0, phtotvz = 0.0, phu = 0.0, phv = 0.0, phzu = 0.0, phzv = 0.0;
00430 for(int i=0; i<stripu.GetEntriesFast(); i++) {
00431 CandStripHandle *strip = dynamic_cast<CandStripHandle*>(stripu.At(i));
00432 phtotuz += strip->GetCharge(CalDigitType::kNone);
00433 phu += strip->GetCharge(CalDigitType::kNone)*strip->GetTPos();
00434 phzu += strip->GetCharge(CalDigitType::kNone)*strip->GetZPos();
00435 }
00436 for(int i=0; i<stripv.GetEntriesFast(); i++) {
00437 CandStripHandle *strip = dynamic_cast<CandStripHandle*>(stripv.At(i));
00438 phtotvz += strip->GetCharge(CalDigitType::kNone);
00439 phv += strip->GetCharge(CalDigitType::kNone)*strip->GetTPos();
00440 phzv += strip->GetCharge(CalDigitType::kNone)*strip->GetZPos();
00441 }
00442
00443 Float_t Upos = -10000.0, ZUpos = -10000.0, Vpos = -10000.0, ZVpos = -10000.0;
00444 if(phtotuz > 0) {
00445 Upos = phu / phtotuz;
00446 ZUpos = phzu / phtotuz;
00447 }
00448 if(phtotvz > 0) {
00449 Vpos = phv / phtotvz;
00450 ZVpos = phzv / phtotvz;
00451 }
00452
00453 Float_t minu = 100000.0, minv = 100000.0; Int_t indexu = 0, indexv = 0;
00454 Float_t BegTimeParent = 30000.0;
00455 Int_t goodindex = 0;
00456 Double_t phdadu = -10.0, phdadv = -10.0;
00457
00458 CandSliceHandleItr sliceItr4(ch.GetDaughterIterator());
00459 CandSliceHandleKeyFunc *sliceKf4 = sliceItr4.CreateKeyFunc();
00460 sliceKf4->SetFun(SliceSRKeyFromForwardTime);
00461 sliceItr4.GetSet()->AdoptSortKeyFunc(sliceKf4);
00462 sliceKf4 = 0;
00463
00464 sliceItr4.Reset();
00465 while(CandSliceHandle *parent = sliceItr4()) {
00466 if(goodindex != badindex) {
00467 if(parent->GetCharge(CalDigitType::kNone) >= 10000) {
00468 CandStripHandleItr parentitr(parent->GetDaughterIterator());
00469 while(CandStripHandle *parentstrip = parentitr()) {
00470 if(parentstrip->GetTime()*1e9 < BegTimeParent)
00471 BegTimeParent = parentstrip->GetTime()*1e9;
00472 }
00473 if(BegTimeParent < BegTimeBad) {
00474 stripu.Clear(); stripv.Clear();
00475 parentitr.Reset();
00476 while(CandStripHandle *curr = parentitr()) {
00477 if(curr->GetPlane()<121 && curr->GetCharge()>2.0) {
00478 if(curr->GetPlaneView() == 2) {
00479 stripu.Add(curr);
00480 } else {
00481 stripv.Add(curr);
00482 }
00483 }
00484 }
00485
00486 Float_t tempu = 300000.0, tempv = 300000.0;
00487
00488 for(int i=0; i<stripu.GetEntriesFast(); i++) {
00489 CandStripHandle *strip = dynamic_cast<CandStripHandle*>(stripu.At(i));
00490 if(sqrt(pow(Upos-strip->GetTPos(),2)+pow(ZUpos-strip->GetZPos(),2)) < tempu)
00491 tempu = sqrt(pow(Upos-strip->GetTPos(),2)+pow(ZUpos-strip->GetZPos(),2));
00492 }
00493 for(int i=0; i<stripv.GetEntriesFast(); i++) {
00494 CandStripHandle *strip = dynamic_cast<CandStripHandle*>(stripv.At(i));
00495 if(sqrt(pow(Vpos-strip->GetTPos(),2)+pow(ZVpos-strip->GetZPos(),2)) < tempv)
00496 tempv = sqrt(pow(Vpos-strip->GetTPos(),2)+pow(ZVpos-strip->GetZPos(),2));
00497 }
00498 stripu.Clear(); stripv.Clear();
00499 if(tempu < minu) {
00500 minu = tempu;
00501 indexu = goodindex;
00502 phdadu = parent->GetCharge(CalDigitType::kNone);
00503 }
00504 if(tempv < minv) {
00505 minv = tempv;
00506 indexv = goodindex;
00507 phdadv = parent->GetCharge(CalDigitType::kNone);
00508 }
00509
00510 }
00511 } //Loop over slices with ph > 10k, trying to find parent
00512 }
00513 goodindex++;
00514 } //Loop over all slices
00515
00516 Float_t minofuandv = 0.0; Int_t minofindexuandv = 0;
00517 if(minu < minv) {
00518 minofuandv = minu;
00519 minofindexuandv = indexu;
00520 } else {
00521 minofuandv = minv;
00522 minofindexuandv = indexv;
00523 }
00524
00525 if(minofuandv <= mincombiningdistance) {
00526 CandSliceHandleItr sliceItr5(ch.GetDaughterIterator());
00527 CandSliceHandleKeyFunc *sliceKf5 = sliceItr5.CreateKeyFunc();
00528 sliceKf5->SetFun(SliceSRKeyFromForwardTime);
00529 sliceItr5.GetSet()->AdoptSortKeyFunc(sliceKf5);
00530 sliceKf5 = 0;
00531
00532 goodindex = 0;
00533 sliceItr5.Reset();
00534 while(CandSliceHandle *slice = sliceItr5()) {
00535 if(goodindex == minofindexuandv) {
00536 CandStripHandleItr badstripitr(slice->GetDaughterIterator());
00537 while(CandStripHandle *badstrip = badstripitr()) {
00538 slice->AddDaughterLink(*badstrip);
00539 }
00540 ch.RemoveDaughter(bad);
00541 break;
00542 }
00543 goodindex++;
00544 }
00545 } //Combining Loop
00546 } //Loop over slices with ph < 10k, trying to find low-C slices
00547 badindex++;
00548 } //end loop over all slices
00549 } // Combining On if statement
00550 CandSliceHandleItr sliceItr10(ch.GetDaughterIterator());
00551 Int_t nslice=0;
00552 MSG("AlgSliceSR",Msg::kDebug) << "Final Slice Contents " << endl;
00553 while(CandSliceHandle *slice = sliceItr10()) {
00554 MSG("AlgSliceSR",Msg::kDebug) << " ** Slice ** " << nslice << endl;
00555 CandStripHandleItr liststripitr(slice->GetDaughterIterator());
00556 while(CandStripHandle *strip = liststripitr()) {
00557 MSG("AlgSliceSR",Msg::kDebug) << "strip plane " << strip->GetPlane() << " strip " << strip->GetStrip() << " charge " << strip->GetCharge() <<" time " << strip->GetCorrBegTime()*1e9 << endl;
00558 }
00559 nslice++;
00560 }
00561
00562
00563 } //End Slice the Snarl
|
|
||||||||||||||||
|
Really simple slicing which gives lower efficiency and higher purity than SR. As it is, it is difficult to tune it further to get better results. Sort Strips by BegTime and simply cut when difference between one strip and the previous one is greater than one bycket. Simple and quick way with results slighly lower in efficiency than SR and higher in purity. Overall they are comparable. Definition at line 572 of file AlgSliceSRList.cxx. References CandHandle::AddDaughterLink(), Registry::Clear(), AlgFactory::GetAlgHandle(), CandStripHandle::GetBegTime(), CandStripHandle::GetCharge(), CandContext::GetDataIn(), CandHandle::GetDaughterIterator(), Registry::GetDouble(), AlgFactory::GetInstance(), Registry::GetInt(), CandContext::GetMom(), CandSlice::MakeCandidate(), and StripSRKeyFromBegTime(). Referenced by RunAlg(). 00573 {
00574
00575 cout << " IN SLICE ASAP " <<endl;
00576 const CandStripListHandle *cslh =
00577 dynamic_cast<const CandStripListHandle*>(cx.GetDataIn());
00578
00579 // Sort by Plane and Strip
00580 CandStripHandleItr cshItr(cslh->GetDaughterIterator());
00581 CandStripHandleKeyFunc *cshKf = cshItr.CreateKeyFunc();
00582 cshKf->SetFun(StripSRKeyFromBegTime);
00583 cshItr.GetSet()->AdoptSortKeyFunc(cshKf);
00584 cshKf = 0;
00585
00586 AlgFactory &af = AlgFactory::GetInstance();
00587 AlgHandle ah = af.GetAlgHandle("AlgSliceSR","default");
00588 CandContext cxx(this,cx.GetMom());
00589
00590 Double_t nbuck = ac.GetDouble("NBuck");
00591 Int_t minstrip = ac.GetInt("MinStrip");
00592 Double_t mincharge = ac.GetDouble("MinCharge");
00593
00594 TObjArray striparraysl;
00595 striparraysl.Clear();
00596
00597 Double_t prevtime;
00598 Bool_t first=true;
00599 Double_t timediff=0;
00600 Double_t timeu;
00601 Int_t ncount=-1;
00602 Int_t nentries=-1;
00603 Double_t cuttime=nbuck*19.*Munits::ns;
00604 // NBuck is be number of buckets
00605 while(CandStripHandle *currstrip = cshItr()){
00606 if(currstrip->GetCharge()>=mincharge) nentries=nentries+1;
00607 }
00608
00609 cshItr.Reset();
00610
00611
00612
00613 while(CandStripHandle *currstrip = cshItr()){
00614 if(currstrip->GetCharge()>=mincharge){
00615 ncount=ncount+1;
00616 if(first) {
00617 first=false;
00618 prevtime=currstrip->GetBegTime();
00619 }
00620 timeu=currstrip->GetBegTime();
00621 timediff = timeu-prevtime;
00622 if(timediff<cuttime){
00623 striparraysl.Add(currstrip);
00624 }
00625 if(timediff>=cuttime || nentries==ncount){
00626 if(striparraysl.GetLast()+1>=minstrip){
00627 cxx.SetDataIn(&striparraysl);
00628 CandSliceHandle slicehandle = CandSlice::MakeCandidate(ah,cxx);
00629 ch.AddDaughterLink(slicehandle);
00630 }
00631 striparraysl.Clear();
00632 striparraysl.Add(currstrip);
00633 }
00634 prevtime = timeu;
00635 }
00636 }
00637
00638 }
|
|
||||||||||||||||
|
Slicing using MSTs with a metric definition (d_mst) that is tunable and at the moment takes into account time and z difference of strips. Gives results the same or slightly better than the SR in both purity and completeness and can be further tuned in various ways. Quite sophisticated algorithm using MSTs that performs slicing with results comparable or some times better than the SR. The MaxLen and Zfact are tunable parameters. I (Niki) will soon describe that in a write up the describes the whole method. Definition at line 648 of file AlgSliceSRList.cxx. References CandHandle::AddDaughterLink(), Registry::Clear(), AlgFactory::GetAlgHandle(), CandStripHandle::GetCharge(), CandContext::GetDataIn(), CandHandle::GetDaughterIterator(), Registry::GetDouble(), AlgFactory::GetInstance(), Registry::GetInt(), CandContext::GetMom(), CandSlice::MakeCandidate(), and StripSRKeyFromCorrTime(). Referenced by RunAlg(). 00649 {
00650 cout << " IN SLICE MST " << endl;
00651 const CandStripListHandle *cslh =
00652 dynamic_cast<const CandStripListHandle*>(cx.GetDataIn());
00653 // Sort by Plane and Strip
00654 CandStripHandleItr cshItr(cslh->GetDaughterIterator());
00655 CandStripHandleKeyFunc *cshKf = cshItr.CreateKeyFunc();
00656 cshKf->SetFun(StripSRKeyFromCorrTime);
00657 cshItr.GetSet()->AdoptSortKeyFunc(cshKf);
00658 cshKf = 0;
00659
00660 AlgFactory &af = AlgFactory::GetInstance();
00661 AlgHandle ah = af.GetAlgHandle("AlgSliceSR","default");
00662 CandContext cxx(this,cx.GetMom());
00663
00664 Int_t minstrip = ac.GetInt("MinStrip");
00665 Double_t mincharge = ac.GetDouble("MinCharge");
00666 Double_t zfact = ac.GetDouble("Zfact");
00667 Double_t maxlen = ac.GetDouble("MaxLen");
00668
00669 cout << " zfact " << zfact << " maxlen " << maxlen << endl;
00670 TObjArray striparrayn;
00671 striparrayn.Clear();
00672 Int_t n_strip=-1;
00673 while(CandStripHandle *currstrip = cshItr()){
00674 if(currstrip->GetCharge()>=mincharge){
00675 striparrayn.Add(currstrip);
00676 n_strip++;
00677 }
00678 }
00679
00680 // Create lower triangular matrix of MST distances.
00681
00682 Int_t n=-1;
00683 Double_t cc=3.*1.e8;
00684 Int_t kMax_d=((n_strip+1)*n_strip)/2;
00685 Double_t *d_mst = new Double_t[kMax_d];
00686 for(Int_t i=0;i<=(n_strip-1);i++){ //
00687 for(Int_t j=0;j<=n_strip;j++) {
00688 if((i+1)>j) {
00689 n++;
00690 Double_t zpos1=dynamic_cast<CandStripHandle*>(striparrayn.At(i+1))->GetZPos();
00691 Double_t zpos2=dynamic_cast<CandStripHandle*>(striparrayn.At(j))->GetZPos();
00692 Double_t t1=dynamic_cast<CandStripHandle*>(striparrayn.At(i+1))->GetCorrBegTime();
00693 Double_t t2=dynamic_cast<CandStripHandle*>(striparrayn.At(j))->GetCorrBegTime();
00694 d_mst[n]=sqrt((t1*cc-t2*cc)*(t1*cc-t2*cc))+((zpos1-zpos2)*(zpos1-zpos2)*zfact*zfact);
00695 }
00696 }
00697 }
00698 //
00699 // Create the MST
00700 Bool_t *a_mst = new Bool_t [n_strip+1];
00701 Int_t *b_mst = new Int_t [n_strip+1];
00702 Double_t *c_mst = new Double_t [n_strip+1];
00703 Int_t *hist_mst = new Int_t [n_strip+1];
00704 Int_t *route_mst = new Int_t [n_strip+1];
00705 Int_t *bb_mst = new Int_t [n_strip+1];
00706 Int_t *ii_mst = new Int_t [n_strip+1];
00707 Double_t *cc_mst = new Double_t [n_strip+1];
00708
00709 Double_t dlarge=99999999.;
00710 Double_t am,dist;
00711 Int_t next;
00712 // Initialization
00713 for(Int_t i=1;i<=n_strip;i++){
00714 a_mst[i]=true;
00715 b_mst[i]=0;
00716 c_mst[i]=dlarge;
00717 }
00718 // Make the tree
00719 Int_t j=0;
00720 Int_t l;
00721 for(Int_t i=1;i<=n_strip;i++){
00722 am=dlarge;
00723 for(Int_t k=1;k<=n_strip;k++){
00724 if(a_mst[k]) {
00725 if(j>k) l=(j-0)*(j-1)/2+k;
00726 if(j<=k) l=(k-0)*(k-1)/2+j;
00727 dist=d_mst[l];
00728 if(dist<c_mst[k]) {
00729 c_mst[k]=dist;
00730 b_mst[k]=j;
00731 }
00732 if(am>c_mst[k]){
00733 am=c_mst[k];
00734 next=k;
00735 }
00736 }
00737 }
00738 j=next;
00739 a_mst[j]=false;
00740 }
00741
00742 delete [] d_mst;
00743 delete [] a_mst;
00744
00745 // DO MST REORDERING
00746
00747 for (Int_t i=0;i<=n_strip;i++){
00748 hist_mst[i]=0;
00749 }
00750 for (Int_t i=1;i<=n_strip;i++){
00751 j=b_mst[i];
00752 hist_mst[j]=hist_mst[j]+1;
00753 }
00754 route_mst[0]=0;
00755 j=0;
00756 Int_t k=0;
00757 Int_t kkk=0;
00758 for(Int_t i=1;i<=n_strip;i++){
00759 while(hist_mst[k]==0){
00760 j=j-1;
00761 k=route_mst[j];
00762 }
00763 hist_mst[k]=hist_mst[k]-1;
00764
00765 Int_t m=1;
00766 Bool_t foundk=false;
00767 while(m<=n_strip && !foundk){
00768 if(k==b_mst[m]) {
00769 foundk=true;
00770 }
00771 else m=m+1;
00772 }
00773 kkk=kkk+1;
00774 bb_mst[kkk]=k;
00775 ii_mst[kkk]=m;
00776 cc_mst[kkk]=c_mst[m];
00777
00778 j=j+1;
00779 route_mst[j]=m;
00780 k=m;
00781 b_mst[m]=-b_mst[m];
00782 if(b_mst[m]==0) b_mst[m]=9999999;
00783 }
00784 delete [] b_mst;
00785 delete [] c_mst;
00786 delete [] hist_mst;
00787 delete [] route_mst;
00788
00789 // NOW CUT BRANCHES AND MAKE THE CLUSTERS
00790
00791 TObjArray striparraysl;
00792 striparraysl.Clear();
00793
00794 for(Int_t i=1;i<=n_strip;i++){
00795 if(cc_mst[i]<maxlen) {
00796 if(i==1) {
00797 striparraysl.Add(dynamic_cast<CandStripHandle*>(striparrayn.At(ii_mst[i])));
00798 striparraysl.Add(dynamic_cast<CandStripHandle*>(striparrayn.At(bb_mst[i])));
00799 }
00800 else striparraysl.Add(dynamic_cast<CandStripHandle*>(striparrayn.At(ii_mst[i])));
00801 }
00802 else if(cc_mst[i]>=maxlen || i==n_strip) {
00803 if(i==1) {
00804 striparraysl.Add(dynamic_cast<CandStripHandle*>(striparrayn.At(ii_mst[i])));
00805 }
00806 else {
00807 // new slice
00808 if(striparraysl.GetLast()+1>=minstrip){
00809 cxx.SetDataIn(&striparraysl);
00810 CandSliceHandle slicehandle = CandSlice::MakeCandidate(ah,cxx);
00811 ch.AddDaughterLink(slicehandle);
00812 }
00813 striparraysl.Clear();
00814 striparraysl.Add(dynamic_cast<CandStripHandle*>(striparrayn.At(ii_mst[i])));
00815 }
00816 }
00817 }
00818
00819 delete [] bb_mst;
00820 delete [] cc_mst;
00821 delete [] ii_mst;
00822 }
|
|
|
Reimplemented from AlgBase. Definition at line 851 of file AlgSliceSRList.cxx. 00852 {
00853 }
|
1.3.9.1