Changeset 15368 for fact/tools/rootmacros
- Timestamp:
- 04/19/13 01:37:40 (12 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
fact/tools/rootmacros/zerosearch.C
r15119 r15368 19 19 int VerbosityLevel 20 20 ){ 21 22 21 vector<Region> * ZeroPositions = new vector<Region>; 23 Region currentRegion = {0, 0, 0, 0.0, 0, 0.0, 0.0, 0}; 22 Region currentRegion = {0, 0, 0, 0.0, 0, 0.0, 0.0, 0, 0}; 23 24 24 if (step < 1){ 25 25 if (VerbosityLevel > 0) … … 67 67 } 68 68 } 69 70 69 } 71 70 else // search for zero-x-ings on both esges … … 81 80 } 82 81 } 83 84 82 85 83 return ZeroPositions; … … 172 170 << reg->maxPos << " with Value:" << reg->maxVal << endl; 173 171 } 172 173 if (reg->maxPos <= 0) 174 { 175 reg = regions.erase( reg ) ; 176 --reg; 177 continue; 178 } 174 179 } 175 180 return regions.size(); … … 411 416 */ 412 417 418 //size_t calcAmplitudeWeightedTime( 419 // vector<Region> ®ions, 420 // vector<float> &data, 421 422 // ){ 423 424 //} 425 426 size_t calcCFDPositions( 427 vector<Region> ®ions, 428 vector<float> &cfd_data 429 ) 430 { 431 vector<Region>::iterator reg; 432 for (reg = regions.begin() ; reg < regions.end() ; reg++) 433 { 434 int neg_zero_crossing = -1; 435 int min_in_cfd = -1; 436 int pos_zero_crossing = reg->maxPos; 437 float minVal = 10.0; 438 439 for (int i = pos_zero_crossing; i > pos_zero_crossing - 40; i--) 440 { 441 if (i <= 1 || i >= 1439) break; 442 443 if ( cfd_data[i] < minVal) 444 { 445 minVal = cfd_data[i]; 446 min_in_cfd = i; 447 } 448 449 if (cfd_data[i] > 0 && cfd_data[i+1] <= 0 && cfd_data[i+2] < 0 ) 450 { 451 neg_zero_crossing = i+1; 452 } 453 } 454 455 reg->cfdMinPos = min_in_cfd; 456 reg->cfdNegZerocrossing = neg_zero_crossing; 457 } 458 return regions.size(); 459 } 460 461 413 462 414 463 size_t findTimeOfHalfMaxLeft( 415 vector<Region> ®ions,416 vector<float> &data,417 float baseline,418 int beginRisingEdge,419 int endRisingEdge,420 int VerbosityLevel)421 { 422 float pulse_size = 0; // 423 float thr = 0; 424 464 vector<Region> ®ions, 465 vector<float> &data, 466 float baseline, 467 int beginRisingEdge, 468 int endRisingEdge, 469 int VerbosityLevel 470 ) 471 { 472 // float pulse_size = 0; // 473 float thr = 0; 425 474 //int begin = beginRisingEdge; 426 int end = endRisingEdge;427 int counter= 0;475 int end = 20; 476 int counter = 0; 428 477 // local vars for line fitting 429 float xmean= 0.0;430 float ymean= 0.0;431 float cov = 0.0; // empiric covarianve between x and y432 float var= 0.0; // empiric variance of x478 float xmean = 0.0; 479 float ymean = 0.0; 480 float cov = 0.0; // empiric covariance between x and y 481 float var = 0.0; // empiric variance of x 433 482 // result of line fitting 434 float slope= 0.0;435 float intercept= 0.0;483 float slope = 0.0; 484 float intercept = 0.0; 436 485 // result of this function -- the position of the half rising edge 437 float pos_of_thr_Xing = 0.0; 438 int result = 0; 486 float pos_of_thr_Xing = 0.0; 487 int result = 0; 488 439 489 440 490 vector<Region>::iterator reg; 441 for (reg = regions.begin() ; reg < regions.end() ; reg++) {442 counter = 0;491 for (reg = regions.begin() ; reg < regions.end() ; reg++) 492 { 443 493 // check if linear rising edge is completely in pipeline 444 494 // if not --> delete … … 450 500 continue; 451 501 } 452 453 // The maximum was probably found using smoothed data,454 // so I read the value at the position of the maximum from the data455 // again. So I rely on a well defined maxPos.456 if (VerbosityLevel > 1) cout << "## baseline: " << baseline << endl;457 pulse_size = data[reg->maxPos] - baseline;458 if (VerbosityLevel > 1) cout << "## pulse_size: " << pulse_size << endl;459 thr = pulse_size / 2. + baseline;460 if (VerbosityLevel > 1) cout << "## thr: " << thr << endl;461 502 462 503 // I think 5 slices left of the max there begins the rising edge … … 469 510 // ALARM check for out of range values... 470 511 512 float average_of_max = -1; 513 int stepwidth = 2; 514 515 counter = 0; 516 for ( int i = reg->maxPos -stepwidth; i <= reg->maxPos + stepwidth; i++) 517 { 518 average_of_max += data[i]; 519 counter++; 520 } 521 average_of_max /= counter; 522 523 thr = average_of_max / 2.0; 524 525 float max_slope = 0; 526 int max_slope_pos = -1; 527 int half_height_pos = -1; 528 471 529 bool passed = false; 472 for (int slice=reg->maxPos; 473 slice > 1; --slice) 474 { 475 if ( data[slice] < 0.9*data[reg->maxPos] && !passed) 476 { 477 beginRisingEdge = reg->maxPos - slice; 478 passed = true; 479 } 480 481 if ( data[slice] < 0.1*data[reg->maxPos]) 482 { 483 endRisingEdge = reg->maxPos - slice; 484 passed = false; 485 reg->lengthOfRisingEdge=endRisingEdge - beginRisingEdge; 486 break; 487 } 488 } 489 490 530 for (int slice=reg->maxPos; slice > reg->maxPos - end; --slice) 531 { 532 if (slice <= 0) break; 533 // if ( data[slice] < 0.9*data[reg->maxPos] && !passed) 534 // { 535 // endRisingEdge = slice; 536 // beginRisingEdge = endRisingEdge -3; 537 // passed = true; 538 // } 539 // else if ( data[slice] < 0) //0.4*data[reg->maxPos] ) 540 // { 541 // beginRisingEdge = slice; 542 // reg->lengthOfRisingEdge=endRisingEdge - beginRisingEdge; 543 // break; 544 // } 545 546 //calculate max slope on leading edge 547 if (data[slice] - data[slice-4] > max_slope && slice-4 > 0) 548 { 549 max_slope = data[slice] - data[slice-4]; 550 max_slope_pos = slice-2; 551 } 552 553 if ( data[slice] <= thr && !passed){ 554 half_height_pos = slice; 555 passed = true; 556 } 557 } 558 reg->maxSlope = max_slope/4; 559 reg->maxSlopePos = max_slope_pos; 560 reg->halfRisingEdgePos = half_height_pos; 561 reg->halfHeight = data[half_height_pos]; 562 563 result = half_height_pos; 564 565 endRisingEdge = result +2; 566 beginRisingEdge = result -2; 567 568 counter = 0; 491 569 //calculate mean of x and y coordinate 492 for (int slice=reg->maxPos - beginRisingEdge; 493 slice > reg->maxPos - endRisingEdge; --slice) 494 { 570 for (int slice= beginRisingEdge; slice < endRisingEdge; ++slice) 571 { 572 if (slice <= 0) break; 573 495 574 xmean += slice; 496 575 ymean += data[slice]; 497 576 counter++; 498 577 } 499 // xmean /= beginRisingEdge - endRisingEdge;500 // ymean /= beginRisingEdge - endRisingEdge;501 578 xmean /= counter; 502 579 ymean /= counter; … … 508 585 var = 0.0; 509 586 510 for (int slice=reg->maxPos - beginRisingEdge; 511 slice > reg->maxPos - endRisingEdge; --slice) 512 { 587 //calculate covariance and variance 588 for (int slice=beginRisingEdge; slice < endRisingEdge; ++slice) 589 { 590 if (slice <= 0) break; 591 513 592 cov += (data[slice] - ymean) * (slice - xmean); 514 593 var += (slice - xmean) * (slice - xmean); … … 517 596 if (VerbosityLevel > 2) cout << "## var: " << var << endl; 518 597 519 slope = cov / var; 520 intercept = ymean - slope * xmean; 598 slope = cov / var; 599 intercept = ymean - slope * xmean; 600 601 // The maximum was probably found using smoothed data, 602 // so I read the value at the position of the maximum from the data 603 // again. So I rely on a well defined maxPos. 604 // if (VerbosityLevel > 4) cout << "## baseline: " << baseline << endl; 605 // pulse_size = data[reg->maxPos] - baseline; 606 // if (VerbosityLevel > 4) cout << "## pulse_size: " << pulse_size << endl; 607 // thr = pulse_size / 2. + baseline; 608 // if (VerbosityLevel > 4) cout << "## thr: " << thr << endl; 521 609 522 610 // now calculate, where the fittet line crosses the threshold 523 pos_of_thr_Xing = (thr - intercept) / slope;524 result = (int)(pos_of_thr_Xing + 0.5);611 // pos_of_thr_Xing = (thr - intercept) / slope; 612 // result = (int)(pos_of_thr_Xing + 0.5); 525 613 526 614 … … 537 625 } 538 626 539 reg->halfRisingEdgePos = result; 540 541 reg->distanceEdgeToMax = reg->maxPos - result; 542 reg->interceptRisingEdge = intercept; 543 reg->slopeOfRisingEdge = slope; 627 reg->distanceEdgeToMax = reg->maxPos - result; 628 reg->interceptRisingEdge = reg->maxPos - intercept; 629 reg->slopeOfRisingEdge = slope; 630 631 if (reg->halfRisingEdgePos <= 0 || reg->halfRisingEdgePos >= 1024) 632 { 633 reg = regions.erase( reg ) ; 634 --reg; 635 continue; 636 } 544 637 } 545 638 return regions.size(); … … 570 663 size_t GetMaxPositions( 571 664 vector<Region> ®ions, 572 vector<int> & maxPositions,573 int VerbosityLevel) 574 { 575 maxPositions.clear();665 vector<int> &positions, 666 int VerbosityLevel) 667 { 668 positions.clear(); 576 669 vector<Region>::iterator it = regions.begin(); 577 670 while( it != regions.end() ) 578 671 { 579 maxPositions.push_back(it->maxPos);672 positions.push_back(it->maxPos); 580 673 if (VerbosityLevel > 3){ 581 674 cout << "getting maxPos@ " << it->maxPos << "\t"; … … 590 683 size_t GetEdgePositions( 591 684 vector<Region> ®ions, 592 vector<int> & edgePos,593 int VerbosityLevel) 594 { 595 maxPositions.clear();685 vector<int> &positions, 686 int VerbosityLevel) 687 { 688 positions.clear(); 596 689 vector<Region>::iterator it = regions.begin(); 597 690 while( it != regions.end() ) 598 691 { 599 maxPositions.push_back(it->halfRisingEdgePos);692 positions.push_back(it->halfRisingEdgePos); 600 693 if (VerbosityLevel > 3){ 601 694 cout << "getting maxPos@ " << it->halfRisingEdgePos << "\t";
Note:
See TracChangeset
for help on using the changeset viewer.