ChopHelper Class Reference

#include <ChopHelper.h>

List of all members.

Public Member Functions

 ChopHelper (const MomNavigator *mom=0)
virtual ~ChopHelper ()
ChopHelpGetChopHelp (std::vector< CandHandle > candidates)
ChopHelpGetChopHelp (const CandHandle *ev1, const CandHandle *ev2=0, const CandHandle *ev3=0, const CandHandle *ev4=0, const CandHandle *ev5=0)

Private Member Functions

bool ShouldSplit (float this_ph, float next_ph, float d_tmax)

Private Attributes

const MomNavigatorfMom

Detailed Description

Definition at line 18 of file ChopHelper.h.


Constructor & Destructor Documentation

ChopHelper::ChopHelper ( const MomNavigator mom = 0  ) 

Definition at line 46 of file ChopHelper.cxx.

00046                                               : fMom(mom)
00047 {
00048 }

ChopHelper::~ChopHelper (  )  [virtual]

Definition at line 51 of file ChopHelper.cxx.

00052 {
00053 }


Member Function Documentation

ChopHelp * ChopHelper::GetChopHelp ( const CandHandle ev1,
const CandHandle ev2 = 0,
const CandHandle ev3 = 0,
const CandHandle ev4 = 0,
const CandHandle ev5 = 0 
)

Definition at line 113 of file ChopHelper.cxx.

References GetChopHelp().

00118 {
00119   std::vector<CandHandle> v;
00120   if(ev1) v.push_back(*ev1);
00121   if(ev2) v.push_back(*ev2);
00122   if(ev3) v.push_back(*ev3);
00123   if(ev4) v.push_back(*ev4);
00124   if(ev5) v.push_back(*ev5);
00125   
00126   return GetChopHelp(v);
00127 }

ChopHelp * ChopHelper::GetChopHelp ( std::vector< CandHandle candidates  ) 

Algorithm to chop by using the summed energy waveform of the whole calorimeter.

Definition at line 129 of file ChopHelper.cxx.

References ChopHelp::candph, ChopHelp::chopmatch, ChopHelp::chopph, digit(), digits(), done(), ChopHelp::estcompleteness, ChopHelp::estpurity, DigitVector::FillFrom(), fMom, PlexSEIdAltL::GetBestSEId(), CandDigitHandle::GetCharge(), VldContext::GetDetector(), PlexPlaneId::GetPlane(), CandDigitHandle::GetPlexSEIdAltL(), Calibrator::GetTDCFromTime(), CandDigitHandle::GetTime(), Calibrator::Instance(), Msg::kDebug, Msg::kError, Detector::kNear, CalTimeType::kNone, kQieRcid, CalDigitType::kSigCorr, Msg::kWarning, PlexPlaneId::LastPlaneNearCalor(), ChopHelp::matchmatrix, MSG, ShouldSplit(), and size.

Referenced by ChopHelperModule::Ana(), and GetChopHelp().

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

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

Definition at line 89 of file ChopHelper.cxx.

References k1pe.

Referenced by GetChopHelp().

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


Member Data Documentation

const MomNavigator* ChopHelper::fMom [private]

Definition at line 36 of file ChopHelper.h.

Referenced by GetChopHelp().


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

Generated on 22 Nov 2017 for loon by  doxygen 1.6.1