#include <ChopHelper.h>
Public Member Functions | |
| ChopHelper (const MomNavigator *mom=0) | |
| virtual | ~ChopHelper () |
| ChopHelp * | GetChopHelp (std::vector< CandHandle > candidates) |
| ChopHelp * | GetChopHelp (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 MomNavigator * | fMom |
|
|
Definition at line 46 of file ChopHelper.cxx. 00046 : fMom(mom) 00047 { 00048 }
|
|
|
Definition at line 51 of file ChopHelper.cxx. 00052 {
00053 }
|
|
||||||||||||||||||||||||
|
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 }
|
|
|
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(), done(), ChopHelp::estcompleteness, ChopHelp::estpurity, DigitVector::FillFrom(), PlexSEIdAltL::GetBestSEId(), CandDigitHandle::GetCharge(), VldContext::GetDetector(), PlexPlaneId::GetPlane(), CandDigitHandle::GetPlexSEIdAltL(), Calibrator::GetTDCFromTime(), CandDigitHandle::GetTime(), Calibrator::Instance(), kQieRcid, PlexPlaneId::LastPlaneNearCalor(), ChopHelp::matchmatrix, MSG, and ShouldSplit(). 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 }
|
|
||||||||||||||||
|
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 }
|
|
|
Definition at line 36 of file ChopHelper.h. |
1.3.9.1