AlgChopListSharp2 Class Reference

#include <AlgChopListSharp2.h>

Inheritance diagram for AlgChopListSharp2:

AlgBase 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)

Detailed Description

Definition at line 13 of file AlgChopListSharp2.h.


Constructor & Destructor Documentation

AlgChopListSharp2::AlgChopListSharp2 (  ) 

Definition at line 50 of file AlgChopListSharp2.cxx.

00051 {
00052 }

AlgChopListSharp2::~AlgChopListSharp2 (  )  [virtual]

Definition at line 55 of file AlgChopListSharp2.cxx.

00056 {
00057 }


Member Function Documentation

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

Implements AlgBase.

Definition at line 117 of file AlgChopListSharp2.cxx.

References digit(), digits(), done(), Form(), AlgFactory::GetAlgHandle(), CandContext::GetCandRecord(), OscFit::GetCharge(), CandContext::GetDataIn(), AlgFactory::GetInstance(), CandContext::GetMom(), PlexPlaneId::GetPlane(), Calibrator::GetTDCFromTime(), RecMinos::GetVldContext(), Calibrator::Instance(), Msg::kDebug, Msg::kError, Detector::kNear, CalTimeType::kNone, kQieRcid, CalDigitType::kSigCorr, Msg::kWarning, PlexPlaneId::LastPlaneNearCalor(), CandDigitList::MakeCandidate(), MSG, CandHandle::SetName(), and ShouldSplit().

00120 {
00124 
00125   // Config.
00126   bool cDoRecombobulation = false;
00127 
00128   assert(candHandle.InheritsFrom("CandChopListHandle"));
00129   CandChopListHandle &chopList = dynamic_cast<CandChopListHandle &>(candHandle);
00130 
00131    assert(candContext.GetDataIn());
00132    // Check for CandDigitListHandle input
00133    if (!(candContext.GetDataIn()->InheritsFrom("CandDigitListHandle"))) {
00134      MSG("Chop",Msg::kWarning ) << "Data into AlgChopListSharp2 is not a digit list." << std::endl;
00135    }
00136    
00137    const CandDigitListHandle *cdlh_ptr = 
00138      dynamic_cast<const CandDigitListHandle*>(candContext.GetDataIn());
00139    
00140    const MomNavigator *mom = candContext.GetMom();
00141    const RawDigitDataBlock* rddb = DataUtil::GetRawBlock<RawDigitDataBlock>(mom);
00142    if (!rddb) {
00143      MSG("Chop", Msg::kWarning) << "No RawDigitDataBlock in RawRecord." << endl;
00144      return;
00145    }
00146    
00147    // Get setup for the DigitList maker algorithm
00148    AlgFactory &af = AlgFactory::GetInstance();
00149    AlgHandle ah = af.GetAlgHandle("AlgChop","default");
00150    CandContext cxx(this,candContext.GetMom());
00151 
00152    const VldContext &context = *(candContext.GetCandRecord()->GetVldContext());
00153    if(context.GetDetector() != Detector::kNear) 
00154      MSG("Chop",Msg::kError) << "Running the Sharp2 algorithm on FD data is a no-no!" << endl;
00155 
00156    Calibrator& cal = Calibrator::Instance();
00157    UgliGeomHandle ugli(context);
00158   
00159    // Now do the actual slicing.
00160 
00161    // First, make a nice stl vector of the digits.
00162    DigitVector digits(cdlh_ptr);
00163 
00164    UInt_t ndigits = digits.size();
00165    UInt_t nchop = 0;
00166     
00167    // Sort the list by time.
00168    // Not neccessary for this algorithm.
00169    // std::sort(digits.begin(), digits.end(), compareDigitTimes());
00170 
00171 
00172    // Also, I want some other pieces of info:
00173    std::vector<int>    digit_tdc(ndigits);
00174    std::vector<UInt_t> digit_plane(ndigits);
00175    //std::vector<float>  digit_tpos(ndigits);
00176    for(UInt_t i=0;i<ndigits;i++) {
00177      digit_tdc[i] = (cal.GetTDCFromTime(digits[i].GetTime(CalTimeType::kNone), kQieRcid));
00178      digit_plane[i] = digits[i].GetPlexSEIdAltL().GetPlane(); 
00179      //if(digit_plane[i]<=PlexPlaneId::LastPlaneNearCalor())
00180      //  digit_tpos[i]  = ugli.GetStripHandle(digits[i].GetPlexSEIdAltL().GetBestSEId()).GetTPos(); 
00181      //else 
00182      //  digit_tpos[i]  = -999;
00183    }
00184 
00185    // Find first and last times. Add some padding so sertain operations are easier to code.
00186    Int_t tfirst = digit_tdc[0];
00187    Int_t tlast  = digit_tdc[0];
00188    for(UInt_t i=0;i<ndigits;i++) {
00189      if(digit_tdc[i] < tfirst) tfirst = digit_tdc[i];
00190      if(digit_tdc[i] > tlast ) tlast  = digit_tdc[i];
00191    }
00192    tfirst-=5;
00193    tlast +=5;
00194 
00195 
00196    // Make the energy histogram.
00197    MSG("Chop",Msg::kDebug) << "Running Chop_Sharp2" << endl;
00198       
00199    UInt_t numBins = tlast-tfirst;
00200 
00201    // Create the energy-time profile.
00202    std::vector<float> energyVsTime(numBins,0.);
00203    
00204    for(UInt_t idig = 0; idig < ndigits; idig++ ) {
00205      float sigcor = digits[idig].GetCharge(CalDigitType::kSigCorr);
00206      int tdcbin = digit_tdc[idig]-tfirst;
00207      if((tdcbin<0) || ((int)numBins<=tdcbin)) MSG("Chop",Msg::kDebug) << "Whups!" << endl;
00208      else if(digit_plane[idig]<=PlexPlaneId::LastPlaneNearCalor()) {
00209        energyVsTime[digit_tdc[idig]-tfirst] += sigcor; 
00210      }
00211    }
00212 
00213 
00214    
00215    std::vector<int> cuts;
00216    /*
00217    cuts.push_back(0); // cut at start.   
00218    int peak_bin = 0;
00219    bool rising = true;
00220    for(UInt_t bin=0;bin<numBins-1;bin++) {
00221      float dq = energyVsTime[bin+1] - energyVsTime[bin];
00222      float max_climb = 2.5*sqrt(fabs(energyVsTime[bin])/k1pe)*k1pe;
00223      if(max_climb < 4.*k1pe) max_climb = max_climb*2.;
00224 
00225      bool falling = !rising;
00226 
00227      if(falling && ((bin-peak_bin)<8)) 
00228        if(max_climb<5.*k1pe) max_climb = 5.*k1pe;
00229 
00230      if(falling) {
00231        if(dq > max_climb){
00232          cuts.push_back(bin);
00233          rising = true;
00234        }
00235      } 
00236 
00237      if(rising) {
00238        if(dq< 0){
00239          rising = false;
00240          peak_bin = bin;
00241        }
00242      }
00243    }
00244    
00245    cuts.push_back(numBins-1); // cut at the end.
00246    */
00247 
00248    // Used bins:
00249    std::vector<char> binsUsed(numBins,0);
00250    binsUsed[0]=1;
00251    binsUsed[numBins-1]=1;
00252 
00253    cuts.push_back(0);
00254    cuts.push_back(numBins-1);
00255 
00256    
00257 
00258    do {
00259      // Look for biggest peak. 
00260      UInt_t biggest_bin  = 99999;
00261      float  biggest_size = 0;
00262      for(UInt_t i=0;i<numBins;i++) {
00263        if(binsUsed[i]==0) 
00264          if(energyVsTime[i]>biggest_size) {
00265            biggest_size = energyVsTime[i];
00266            biggest_bin = i;
00267          }
00268      }
00269      
00270      if(biggest_bin==99999) break; // We've gone through all of them.
00271      if(biggest_size<=0.) break; // We've hit rock bottom.
00272      
00273      // Collect the start and stop time for this chop.
00274      // Start 1 bin before the peak, and at least 1 bin after the peak.
00275      UInt_t bin_start = biggest_bin;
00276      UInt_t bin_stop =  biggest_bin;
00277      
00278      if(binsUsed[bin_start-1]==0) bin_start--;
00279      if(binsUsed[bin_stop +1]==0) bin_stop++;
00280      
00281      bool done = false;
00282      while(!done) {
00283 
00284        // Stop if there's a rise. If so, mark as a contested bin.
00285        if(ShouldSplit(energyVsTime[bin_start], 
00286                       energyVsTime[bin_start-1],
00287                       bin_start-1 - biggest_bin ) ) {
00288          done = true;
00289          cuts.push_back(bin_start);
00290        }
00291        
00292        // Stop if we've hit another chop.
00293        if(binsUsed[bin_start-1]) {
00294          done = true;
00295        }
00296        
00297        if(!done) {
00298          bin_start--;
00299        }
00300      };
00301      
00302      // Expand forwards until the energy starts climbing.
00303      // But, allow small pulses in for the first 5 buckets.
00304      done = false;
00305      while(!done) {
00306        
00307        // Allow 5 buckets worth of small stuff:
00308        if((energyVsTime[bin_stop+1] < 500.) && (bin_stop+1 < biggest_bin+7)) {
00309          // keep going
00310        } else {
00311          if(ShouldSplit(energyVsTime[bin_stop],
00312                         energyVsTime[bin_stop+1],
00313                         bin_stop+1 - biggest_bin) ) {
00314            done = true;
00315            cuts.push_back(bin_stop); // mark this one as contested.
00316          }
00317        }
00318        
00319        // Stop if we hit another chop.
00320        if(binsUsed[bin_stop+1]) {
00321          //MSG("Chop",Msg::kDebug) << "Didn't move forward; hit another chop." << endl;
00322          done = true;
00323        }
00324        
00325        // If we're ok, increment and continue.
00326        if(!done) bin_stop++;
00327      }
00328      
00329      // Zero out these buckets so they won't be caught again.
00330      for(UInt_t i=bin_start; i<=bin_stop; i++) {
00331        binsUsed[i] = 1;
00332      }
00333      
00334    } while(true); 
00335 
00336 
00337    // Time-order the cuts.
00338    std::sort(cuts.begin(), cuts.end());
00339 
00340    // One chop for every cut, corresponding to the cut that ends the chop.
00341    std::vector<Chop> chops(cuts.size());
00342 
00343    // Build initial cuts, excluding hits in contested bins.
00344    // loop cuts, 1 to n
00345    for(UInt_t icut = 1; icut<= cuts.size(); icut++) {
00346      int bin_start = cuts[icut-1]+1; // Excluding last contested bin
00347      int bin_stop  = cuts[icut]-1;   // Excluding current contested bin.
00348 
00349      int tdc_start = bin_start + tfirst;
00350      int tdc_stop  = bin_stop  + tfirst;
00351      
00352      for(UInt_t idig = 0; idig < ndigits; idig++ ) {
00353        int tdc = digit_tdc[idig];
00354        if((tdc>=(int)tdc_start)&&(tdc<=(int)tdc_stop)) {
00355          chops[icut].Add(digits[idig]);
00356        }
00357      }     
00358    }
00359 
00360    // Now loop through again and deal with any contested bins.
00361    // Ignore the first and last cuts, since they are at the start
00362    // and end of the event.
00363    for(UInt_t icut = 1; icut< cuts.size()-1; icut++) {
00364      int contested_bin = cuts[icut];
00365      int contested_tdc = contested_bin + tfirst;
00366      
00367      // Make a list of digits in the contested bin.
00368      DigitVector contested_digits;
00369      for(UInt_t idig = 0; idig < ndigits; idig++ ) {
00370        int tdc = digit_tdc[idig];
00371        if(tdc==contested_tdc) {
00372          contested_digits.push_back(digits[idig]);
00373        }
00374      }
00375 
00376      if(contested_digits.size()>0) {
00377        // only bother if there's actually some digits in the contested bin.
00378 
00379        int chop1 = icut;
00380        int chop2 = icut+1;
00381 
00382        // Tell the fighting chops to build maps.
00383        chops[chop1].BuildMaps();
00384        chops[chop2].  BuildMaps();
00385 
00386        MSG("Chop",Msg::kDebug) << "Dealing with contested bin " << cuts[icut] 
00387             << " with " << contested_digits.size() << " contested digits." << endl;
00388 
00389        int assign_strip1 = 0;
00390        int assign_strip2 = 0;
00391        int assign_plane1 = 0;
00392        int assign_plane2 = 0;
00393        int assign_total1 = 0;
00394        int assign_total2 = 0;
00395 
00396        
00397        // Now loop through contested digits and assign them:
00398        for(UInt_t ic=0; ic<contested_digits.size(); ic++) {
00399          // Assign to chop with most energy on same strip:
00400          CandDigitHandle digit = contested_digits[ic];
00401          PlexStripEndId seid =digit.GetPlexSEIdAltL().GetBestSEId();
00402          float chop1stripE = chops[chop1].fStripEnergy[seid];
00403          float chop2stripE  =chops[chop2  ].fStripEnergy[seid];
00404 
00405          /*
00406          MSG("Chop",Msg::kDebug) << "  Consider " << seid.AsString() << " q=" << digit.GetCharge(CalDigitType::kSigCorr)
00407               << "  chop1(s/p/t) = " 
00408               << chops[chop1].fStripEnergy[seid] << "/"
00409               << chops[chop1].fPlaneEnergy[seid.GetPlane()] << "/"
00410               << chops[chop1].fTotalEnergy 
00411               << "  chop2(s/p/t) = " 
00412               << chops[chop2].fStripEnergy[seid] << "/"
00413               << chops[chop2].fPlaneEnergy[seid.GetPlane()] << "/"
00414               << chops[chop2].fTotalEnergy 
00415               << endl;
00416          */
00417 
00418          if(chop1stripE > chop2stripE){
00419            chops[chop1].Add(digit);
00420            assign_strip1++;
00421          } else if(chop2stripE > chop1stripE){
00422            chops[chop2].Add(digit);
00423            assign_strip2++;
00424          } else {
00425            // Darn. The strip content is the same in both. (probably 0)
00426            // So compare on the plane level:
00427            
00428            float chop1planeE = chops[chop1].fPlaneEnergy[seid.GetPlane()];
00429            float chop2planeE = chops[chop2].fPlaneEnergy[seid.GetPlane()];
00430            if(chop1planeE > chop2planeE){
00431              chops[chop1].Add(digit);
00432              assign_plane1++;
00433            } else  if(chop2planeE > chop1planeE){
00434              chops[chop2].Add(digit);
00435              assign_plane2++;      
00436            } else {
00437 
00438              // Ok, so there's nothing in either plane. Assign it to the biggest chop,
00439              // or the earlier one if tied.
00440              if(chops[chop1].fTotalEnergy>=chops[chop2].fTotalEnergy) {
00441                chops[chop1].Add(digit);
00442                assign_total1++;
00443              } else {
00444                chops[chop2].Add(digit);
00445                assign_total2++;
00446              }
00447            }
00448          }
00449        } //uuuuuuugly!
00450        
00451       
00452        MSG("Chop",Msg::kDebug) << "  Assigned by strip: " << assign_strip1 << "/" << assign_strip2
00453             << "   by plane: "         << assign_plane1 << "/" << assign_plane2
00454             << "   by total: "         << assign_total1 << "/" << assign_total2
00455             << endl;
00456      }
00457     
00458 
00459    }
00460    
00462    // And now...
00463    // Recombobulate!
00464    if(cDoRecombobulation) {
00465 
00466    for(UInt_t ichop=1;ichop<chops.size();ichop++) chops[ichop].BuildMaps();
00467    
00468    for(UInt_t icut=1;icut<cuts.size()-1;icut++) {
00469      for(int ilr = 0; ilr<2; ilr++) { //loop left and right
00470        int contested_tdc;
00471        int chopfrom;
00472        int chopto;
00473        if(ilr==0) {
00474          contested_tdc = cuts[icut]-1 + tfirst;
00475          chopfrom      = icut;
00476          chopto        = icut+1;
00477        } else {
00478          contested_tdc = cuts[icut]+1 + tfirst;
00479          chopfrom      = icut+1;
00480          chopto        = icut;
00481        }
00482 
00483        int moved = 0;
00484        
00485        for(UInt_t idig = 0; idig<chops[chopfrom].size(); idig++) {
00486          CandDigitHandle digit = chops[chopfrom][idig];
00487          int tdc = (cal.GetTDCFromTime(digit.GetTime(CalTimeType::kNone), kQieRcid));
00488          if(tdc == contested_tdc) {
00489            PlexStripEndId seid =digit.GetPlexSEIdAltL().GetBestSEId();
00490            float sigcor = digits[idig].GetCharge(CalDigitType::kSigCorr);
00491            
00492            // If the chop that owns the digit has no other digits on that strip:
00493            if(chops[chopfrom].fStripEnergy[seid] <= sigcor) {
00494              // If the neigboring chop HAS got some energy on that strip
00495              if(chops[chopto].fStripEnergy[seid] > 0) {
00496                // Then move the digit.
00497                chops[chopto].push_back(digit);
00498                chops[chopfrom].erase(chops[chopfrom].begin() + idig);
00499                idig--; // set so we look at the next element instead of skipping one.
00500                moved++;
00501              }
00502            }
00503          }
00504        }
00505        
00506        MSG("Chop",Msg::kDebug) << "Recombobulated " << moved << " digits in bin " << contested_tdc-tfirst
00507             << " from chop " << chopfrom << " to chop " << chopto << endl;
00508      } // End left/right loop
00509    } // End loop over cuts
00510 
00511    } // end Recombobulation
00512 
00513 
00515    // Finally, store the chops.
00516    
00517    for(UInt_t ichop=1;ichop<chops.size();ichop++) {
00518 
00519      if(chops[ichop].size()>0) {
00520    
00521        cxx.SetDataIn(&(chops[ichop]));
00522        CandDigitListHandle chopHandle = CandDigitList::MakeCandidate(ah,cxx);
00523        chopHandle.SetName(Form("Chop %d",ichop));
00524        chopList.AddDaughterLink(chopHandle);
00525        
00526        MSG("Chop",Msg::kDebug) << "Creating chop. "
00527             << "  Start: " << cuts[ichop-1] << "   Stop: " << cuts[ichop] 
00528             << "  Digits: " << chops[ichop].size() 
00529             << endl;
00530 
00531        nchop++;
00532      }
00533    } 
00534 }

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

Definition at line 93 of file AlgChopListSharp2.cxx.

References k1pe.

Referenced by RunAlg().

00097 {
00098   float climb = next_ph - this_ph;
00099   
00100   // the maximum delta that the algorithm will climb before making a new chop
00101   float max_climb = 2.5*sqrt(fabs(this_ph)/k1pe)*k1pe;
00102   if(max_climb < (6 * k1pe)) max_climb=max_climb*2;
00103 
00104   // the maximum pulse height in this bin if we're making a new chop.
00105   //const float size_limit = 20 * k1pe; 
00106   //const float min_time = 2;
00107 
00108   //if(d_tmax  < min_time) return false;
00109   //if(this_ph < size_limit) return false;
00110 
00111   if( climb >= max_climb ) return true;
00112   else return false;    
00113 }

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

Reimplemented from AlgBase.

Definition at line 537 of file AlgChopListSharp2.cxx.

00538 {
00539 }


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