Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

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)

Constructor & Destructor Documentation

AlgChopListSharp2::AlgChopListSharp2  ) 
 

Definition at line 49 of file AlgChopListSharp2.cxx.

00050 {
00051 }

AlgChopListSharp2::~AlgChopListSharp2  )  [virtual]
 

Definition at line 54 of file AlgChopListSharp2.cxx.

00055 {
00056 }


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 116 of file AlgChopListSharp2.cxx.

References CandHandle::AddDaughterLink(), digit(), 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().

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

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

Definition at line 92 of file AlgChopListSharp2.cxx.

References k1pe.

Referenced by RunAlg().

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

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

Reimplemented from AlgBase.

Definition at line 542 of file AlgChopListSharp2.cxx.

00543 {
00544 }


The documentation for this class was generated from the following files:
Generated on Sat Nov 21 22:49:14 2009 for loon by  doxygen 1.3.9.1