AlgChopListSharp2 Class Reference

#include <AlgChopListSharp2.h>

Inheritance diagram for AlgChopListSharp2:

List of all members.

Public Member Functions
	AlgChopListSharp2 ()
virtual	~AlgChopListSharp2 ()
virtual void	RunAlg (AlgConfig &ac, CandHandle &ch, CandContext &cx)
virtual void	Trace (const char *c) const
Private Member Functions
bool	ShouldSplit (float this_ph, float next_ph, float d_tmax)

---

Constructor & Destructor Documentation

AlgChopListSharp2::AlgChopListSharp2 (   )

Definition at line 48 of file AlgChopListSharp2.cxx. { }

AlgChopListSharp2::~AlgChopListSharp2 (   )  [virtual]

Definition at line 53 of file AlgChopListSharp2.cxx. { }

---

Member Function Documentation

void AlgChopListSharp2::RunAlg (  AlgConfig &  ac, CandHandle &  candHandle, CandContext &  candContext )  [virtual]

Algorithm to chop by using the summed energy waveform of the whole calorimeter. Implements AlgBase. Definition at line 115 of file AlgChopListSharp2.cxx. References CandHandle::AddDaughterLink(), digits(), done(), RawRecord::FindRawBlock(), Form(), AlgFactory::GetAlgHandle(), PlexSEIdAltL::GetBestSEId(), CandContext::GetCandRecord(), CandContext::GetDataIn(), VldContext::GetDetector(), MomNavigator::GetFragment(), AlgFactory::GetInstance(), CandContext::GetMom(), PlexPlaneId::GetPlane(), CandDigitHandle::GetPlexSEIdAltL(), Calibrator::GetTDCFromTime(), CandDigitHandle::GetTime(), RecMinos::GetVldContext(), Calibrator::Instance(), kQieRcid, PlexPlaneId::LastPlaneNearCalor(), CandDigitList::MakeCandidate(), MSG, CandHandle::SetName(), and ShouldSplit(). { // Config. bool cDoRecombobulation = false; assert(candHandle.InheritsFrom("CandChopListHandle")); CandChopListHandle &chopList = dynamic_cast<CandChopListHandle &>(candHandle); assert(candContext.GetDataIn()); // Check for CandDigitListHandle input if (!(candContext.GetDataIn()->InheritsFrom("CandDigitListHandle"))) { MSG("Chop",Msg::kWarning ) << "Data into AlgChopListSharp2 is not a digit list." << std::endl; } const CandDigitListHandle *cdlh_ptr = dynamic_cast<const CandDigitListHandle*>(candContext.GetDataIn()); const MomNavigator *mom = candContext.GetMom(); RawRecord *rr = dynamic_cast<RawRecord *>(mom->GetFragment("RawRecord")); if (!rr) { MSG("Chop", Msg::kWarning) << "No RawRecord in MOM." << endl; return; } const RawDigitDataBlock *rddb = dynamic_cast<const RawDigitDataBlock *> (rr->FindRawBlock("RawDigitDataBlock")); if (!rddb) { MSG("Chop", Msg::kWarning) << "No RawDigitDataBlock in RawRecord." << endl; return; } // Get setup for the DigitList maker algorithm AlgFactory &af = AlgFactory::GetInstance(); AlgHandle ah = af.GetAlgHandle("AlgChop","default"); CandContext cxx(this,candContext.GetMom()); const VldContext &context = *(candContext.GetCandRecord()->GetVldContext()); if(context.GetDetector() != DetectorType::kNear) MSG("Chop",Msg::kError) << "Running the Sharp2 algorithm on FD data is a no-no!" << endl; Calibrator& cal = Calibrator::Instance(); UgliGeomHandle ugli(context); // Now do the actual slicing. // First, make a nice stl vector of the digits. DigitVector digits(cdlh_ptr); UInt_t ndigits = digits.size(); UInt_t nchop = 0; // Sort the list by time. // Not neccessary for this algorithm. // std::sort(digits.begin(), digits.end(), compareDigitTimes()); // Also, I want some other pieces of info: std::vector<int> digit_tdc(ndigits); std::vector<UInt_t> digit_plane(ndigits); //std::vector<float> digit_tpos(ndigits); for(UInt_t i=0;i<ndigits;i++) { digit_tdc[i] = (cal.GetTDCFromTime(digits[i].GetTime(CalTimeType::kNone), kQieRcid)); digit_plane[i] = digits[i].GetPlexSEIdAltL().GetPlane(); //if(digit_plane[i]<=PlexPlaneId::LastPlaneNearCalor()) // digit_tpos[i] = ugli.GetStripHandle(digits[i].GetPlexSEIdAltL().GetBestSEId()).GetTPos(); //else // digit_tpos[i] = -999; } // Find first and last times. Add some padding so sertain operations are easier to code. Int_t tfirst = digit_tdc[0]; Int_t tlast = digit_tdc[0]; for(UInt_t i=0;i<ndigits;i++) { if(digit_tdc[i] < tfirst) tfirst = digit_tdc[i]; if(digit_tdc[i] > tlast ) tlast = digit_tdc[i]; } tfirst-=5; tlast +=5; // Make the energy histogram. MSG("Chop",Msg::kDebug) << "Running Chop_Sharp2" << endl; UInt_t numBins = tlast-tfirst; // Create the energy-time profile. std::vector<float> energyVsTime(numBins,0.); for(UInt_t idig = 0; idig < ndigits; idig++ ) { float sigcor = digits[idig].GetCharge(CalDigitType::kSigCorr); int tdcbin = digit_tdc[idig]-tfirst; if((tdcbin<0) || ((int)numBins<=tdcbin)) MSG("Chop",Msg::kDebug) << "Whups!" << endl; else if(digit_plane[idig]<=PlexPlaneId::LastPlaneNearCalor()) { energyVsTime[digit_tdc[idig]-tfirst] += sigcor; } } std::vector<int> cuts; /* cuts.push_back(0); // cut at start. int peak_bin = 0; bool rising = true; for(UInt_t bin=0;bin<numBins-1;bin++) { float dq = energyVsTime[bin+1] - energyVsTime[bin]; float max_climb = 2.5*sqrt(fabs(energyVsTime[bin])/k1pe)*k1pe; if(max_climb < 4.*k1pe) max_climb = max_climb*2.; bool falling = !rising; if(falling && ((bin-peak_bin)<8)) if(max_climb<5.*k1pe) max_climb = 5.*k1pe; if(falling) { if(dq > max_climb){ cuts.push_back(bin); rising = true; } } if(rising) { if(dq< 0){ rising = false; peak_bin = bin; } } } cuts.push_back(numBins-1); // cut at the end. */ // Used bins: std::vector<char> binsUsed(numBins,0); binsUsed[0]=1; binsUsed[numBins-1]=1; cuts.push_back(0); cuts.push_back(numBins-1); do { // Look for biggest peak. UInt_t biggest_bin = 99999; float biggest_size = 0; for(UInt_t i=0;i<numBins;i++) { if(binsUsed[i]==0) if(energyVsTime[i]>biggest_size) { biggest_size = energyVsTime[i]; biggest_bin = i; } } if(biggest_bin==99999) break; // We've gone through all of them. if(biggest_size<=0.) break; // We've hit rock bottom. // Collect the start and stop time for this chop. // Start 1 bin before the peak, and at least 1 bin after the peak. UInt_t bin_start = biggest_bin; UInt_t bin_stop = biggest_bin; if(binsUsed[bin_start-1]==0) bin_start--; if(binsUsed[bin_stop +1]==0) bin_stop++; bool done = false; while(!done) { // Stop if there's a rise. If so, mark as a contested bin. if(ShouldSplit(energyVsTime[bin_start], energyVsTime[bin_start-1], bin_start-1 - biggest_bin ) ) { done = true; cuts.push_back(bin_start); } // Stop if we've hit another chop. if(binsUsed[bin_start-1]) { done = true; } if(!done) { bin_start--; } }; // Expand forwards until the energy starts climbing. // But, allow small pulses in for the first 5 buckets. done = false; while(!done) { // Allow 5 buckets worth of small stuff: if((energyVsTime[bin_stop+1] < 500.) && (bin_stop+1 < biggest_bin+7)) { // keep going } else { if(ShouldSplit(energyVsTime[bin_stop], energyVsTime[bin_stop+1], bin_stop+1 - biggest_bin) ) { done = true; cuts.push_back(bin_stop); // mark this one as contested. } } // Stop if we hit another chop. if(binsUsed[bin_stop+1]) { //MSG("Chop",Msg::kDebug) << "Didn't move forward; hit another chop." << endl; done = true; } // If we're ok, increment and continue. if(!done) bin_stop++; } // Zero out these buckets so they won't be caught again. for(UInt_t i=bin_start; i<=bin_stop; i++) { binsUsed[i] = 1; } } while(true); // Time-order the cuts. std::sort(cuts.begin(), cuts.end()); // One chop for every cut, corresponding to the cut that ends the chop. std::vector<Chop> chops(cuts.size()); // Build initial cuts, excluding hits in contested bins. // loop cuts, 1 to n for(UInt_t icut = 1; icut<= cuts.size(); icut++) { int bin_start = cuts[icut-1]+1; // Excluding last contested bin int bin_stop = cuts[icut]-1; // Excluding current contested bin. int tdc_start = bin_start + tfirst; int tdc_stop = bin_stop + tfirst; for(UInt_t idig = 0; idig < ndigits; idig++ ) { int tdc = digit_tdc[idig]; if((tdc>=(int)tdc_start)&&(tdc<=(int)tdc_stop)) { chops[icut].Add(digits[idig]); } } } // Now loop through again and deal with any contested bins. // Ignore the first and last cuts, since they are at the start // and end of the event. for(UInt_t icut = 1; icut< cuts.size()-1; icut++) { int contested_bin = cuts[icut]; int contested_tdc = contested_bin + tfirst; // Make a list of digits in the contested bin. DigitVector contested_digits; for(UInt_t idig = 0; idig < ndigits; idig++ ) { int tdc = digit_tdc[idig]; if(tdc==contested_tdc) { contested_digits.push_back(digits[idig]); } } if(contested_digits.size()>0) { // only bother if there's actually some digits in the contested bin. int chop1 = icut; int chop2 = icut+1; // Tell the fighting chops to build maps. chops[chop1].BuildMaps(); chops[chop2]. BuildMaps(); MSG("Chop",Msg::kDebug) << "Dealing with contested bin " << cuts[icut] << " with " << contested_digits.size() << " contested digits." << endl; int assign_strip1 = 0; int assign_strip2 = 0; int assign_plane1 = 0; int assign_plane2 = 0; int assign_total1 = 0; int assign_total2 = 0; // Now loop through contested digits and assign them: for(UInt_t ic=0; ic<contested_digits.size(); ic++) { // Assign to chop with most energy on same strip: CandDigitHandle digit = contested_digits[ic]; PlexStripEndId seid =digit.GetPlexSEIdAltL().GetBestSEId(); float chop1stripE = chops[chop1].fStripEnergy[seid]; float chop2stripE =chops[chop2 ].fStripEnergy[seid]; /* MSG("Chop",Msg::kDebug) << " Consider " << seid.AsString() << " q=" << digit.GetCharge(CalDigitType::kSigCorr) << " chop1(s/p/t) = " << chops[chop1].fStripEnergy[seid] << "/" << chops[chop1].fPlaneEnergy[seid.GetPlane()] << "/" << chops[chop1].fTotalEnergy << " chop2(s/p/t) = " << chops[chop2].fStripEnergy[seid] << "/" << chops[chop2].fPlaneEnergy[seid.GetPlane()] << "/" << chops[chop2].fTotalEnergy << endl; */ if(chop1stripE > chop2stripE){ chops[chop1].Add(digit); assign_strip1++; } else if(chop2stripE > chop1stripE){ chops[chop2].Add(digit); assign_strip2++; } else { // Darn. The strip content is the same in both. (probably 0) // So compare on the plane level: float chop1planeE = chops[chop1].fPlaneEnergy[seid.GetPlane()]; float chop2planeE = chops[chop2].fPlaneEnergy[seid.GetPlane()]; if(chop1planeE > chop2planeE){ chops[chop1].Add(digit); assign_plane1++; } else if(chop2planeE > chop1planeE){ chops[chop2].Add(digit); assign_plane2++; } else { // Ok, so there's nothing in either plane. Assign it to the biggest chop, // or the earlier one if tied. if(chops[chop1].fTotalEnergy>=chops[chop2].fTotalEnergy) { chops[chop1].Add(digit); assign_total1++; } else { chops[chop2].Add(digit); assign_total2++; } } } } //uuuuuuugly! MSG("Chop",Msg::kDebug) << " Assigned by strip: " << assign_strip1 << "/" << assign_strip2 << " by plane: " << assign_plane1 << "/" << assign_plane2 << " by total: " << assign_total1 << "/" << assign_total2 << endl; } } // And now... // Recombobulate! if(cDoRecombobulation) { for(UInt_t ichop=1;ichop<chops.size();ichop++) chops[ichop].BuildMaps(); for(UInt_t icut=1;icut<cuts.size()-1;icut++) { for(int ilr = 0; ilr<2; ilr++) { //loop left and right int contested_tdc; int chopfrom; int chopto; if(ilr==0) { contested_tdc = cuts[icut]-1 + tfirst; chopfrom = icut; chopto = icut+1; } else { contested_tdc = cuts[icut]+1 + tfirst; chopfrom = icut+1; chopto = icut; } int moved = 0; for(UInt_t idig = 0; idig<chops[chopfrom].size(); idig++) { CandDigitHandle digit = chops[chopfrom][idig]; int tdc = (cal.GetTDCFromTime(digit.GetTime(CalTimeType::kNone), kQieRcid)); if(tdc == contested_tdc) { PlexStripEndId seid =digit.GetPlexSEIdAltL().GetBestSEId(); float sigcor = digits[idig].GetCharge(CalDigitType::kSigCorr); // If the chop that owns the digit has no other digits on that strip: if(chops[chopfrom].fStripEnergy[seid] <= sigcor) { // If the neigboring chop HAS got some energy on that strip if(chops[chopto].fStripEnergy[seid] > 0) { // Then move the digit. chops[chopto].push_back(digit); chops[chopfrom].erase(chops[chopfrom].begin() + idig); idig--; // set so we look at the next element instead of skipping one. moved++; } } } } MSG("Chop",Msg::kDebug) << "Recombobulated " << moved << " digits in bin " << contested_tdc-tfirst << " from chop " << chopfrom << " to chop " << chopto << endl; } // End left/right loop } // End loop over cuts } // end Recombobulation // Finally, store the chops. for(UInt_t ichop=1;ichop<chops.size();ichop++) { if(chops[ichop].size()>0) { cxx.SetDataIn(&(chops[ichop])); CandDigitListHandle chopHandle = CandDigitList::MakeCandidate(ah,cxx); chopHandle.SetName(Form("Chop %d",ichop)); chopList.AddDaughterLink(chopHandle); MSG("Chop",Msg::kDebug) << "Creating chop. " << " Start: " << cuts[ichop-1] << " Stop: " << cuts[ichop] << " Digits: " << chops[ichop].size() << endl; nchop++; } } }

bool AlgChopListSharp2::ShouldSplit (  float  this_ph, float  next_ph, float  d_tmax )  [private]

Definition at line 91 of file AlgChopListSharp2.cxx. References k1pe. Referenced by RunAlg(). { float climb = next_ph - this_ph; // the maximum delta that the algorithm will climb before making a new chop float max_climb = 2.5*sqrt(fabs(this_ph)/k1pe)*k1pe; if(max_climb < (6 * k1pe)) max_climb=max_climb*2; // the maximum pulse height in this bin if we're making a new chop. //const float size_limit = 20 * k1pe; //const float min_time = 2; //if(d_tmax < min_time) return false; //if(this_ph < size_limit) return false; if( climb >= max_climb ) return true; else return false; }

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

Reimplemented from AlgBase. Definition at line 541 of file AlgChopListSharp2.cxx. { }

---

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

• AlgChopListSharp2.h
• AlgChopListSharp2.cxx

---