| | 576 | == Optim Disp Plot1 == |
| | 577 | |
| | 578 | {{{#!Spoiler |
| | 579 | {{{#!cpp |
| | 580 | #include <iostream> |
| | 581 | #include <iomanip> |
| | 582 | |
| | 583 | #include <TMath.h> |
| | 584 | #include <TH2.h> |
| | 585 | #include <TChain.h> |
| | 586 | #include <TCanvas.h> |
| | 587 | |
| | 588 | void optimdisp2() |
| | 589 | { |
| | 590 | // Create chain for the tree Result |
| | 591 | // This is just easier than using TFile/TTree |
| | 592 | TChain c("Result"); |
| | 593 | |
| | 594 | // Add the input file to the |
| | 595 | c.AddFile("simulation.root"); |
| | 596 | |
| | 597 | // Define variables for all leaves to be accessed |
| | 598 | // By definition rootifysql writes only doubles |
| | 599 | double FileId, EvtNumber, X, Y, MeanX, MeanY, Width, Length, CosDelta, SinDelta, |
| | 600 | M3Long, SlopeLong, Leakage1, SlopeSpreadWeighted, Size, |
| | 601 | ConcCore, ConcCOG, NumIslands, NumUsedPixels, Zd; |
| | 602 | |
| | 603 | // Connect the variables to the cordesponding leaves |
| | 604 | c.SetBranchAddress("FileId", &FileId); |
| | 605 | c.SetBranchAddress("EvtNumber", &EvtNumber); |
| | 606 | c.SetBranchAddress("X", &X); |
| | 607 | c.SetBranchAddress("Y", &Y); |
| | 608 | c.SetBranchAddress("MeanX", &MeanX); |
| | 609 | c.SetBranchAddress("MeanY", &MeanY); |
| | 610 | c.SetBranchAddress("Width", &Width); |
| | 611 | c.SetBranchAddress("Length", &Length); |
| | 612 | c.SetBranchAddress("CosDelta", &CosDelta); |
| | 613 | c.SetBranchAddress("SinDelta", &SinDelta); |
| | 614 | c.SetBranchAddress("M3Long", &M3Long); |
| | 615 | c.SetBranchAddress("SlopeLong", &SlopeLong); |
| | 616 | c.SetBranchAddress("Leakage1", &Leakage1); |
| | 617 | c.SetBranchAddress("NumIslands", &NumIslands); |
| | 618 | c.SetBranchAddress("NumUsedPixels", &NumUsedPixels); |
| | 619 | c.SetBranchAddress("Size", &Size); |
| | 620 | |
| | 621 | // Set some constants (they could be included in the database |
| | 622 | // in the future) |
| | 623 | double mm2deg = +0.0117193246260285378; |
| | 624 | |
| | 625 | // -------------------- Source dependent parameter calculation ------------------- |
| | 626 | |
| | 627 | TH2F hist("", "", |
| | 628 | 26, 1.300-0.0025, 1.430-0.0025, |
| | 629 | 20, 0.040-0.0025, 0.140-0.0025); |
| | 630 | |
| | 631 | for (float xi0=1.30; xi0<1.43; xi0+=0.005) |
| | 632 | for (float slope0=0.04; slope0<0.14; slope0+=0.005) |
| | 633 | { |
| | 634 | int cnt_on = 0; |
| | 635 | |
| | 636 | for (int i=0; i<c.GetEntries(); i++) |
| | 637 | { |
| | 638 | // read the i-th event from the file |
| | 639 | c.GetEntry(i); |
| | 640 | |
| | 641 | // First calculate all cuts to speedup the analysis |
| | 642 | double area = TMath::Pi()*Width*Length; |
| | 643 | |
| | 644 | // The abberation correction does increase also Width and Length by 1.02 |
| | 645 | |
| | 646 | bool cutq = NumIslands<3.5 && NumUsedPixels>5.5 && Leakage1==0; |
| | 647 | if (!cutq) |
| | 648 | continue; |
| | 649 | |
| | 650 | bool cut0 = area < log10(Size)*898-1535; |
| | 651 | if (!cut0) |
| | 652 | continue; |
| | 653 | |
| | 654 | int angle = 0; |
| | 655 | |
| | 656 | // -------------------- Source dependent parameter calculation ------------------- |
| | 657 | |
| | 658 | double cr = cos(angle*TMath::DegToRad()); |
| | 659 | double sr = sin(angle*TMath::DegToRad()); |
| | 660 | |
| | 661 | double px = cr*X-sr*Y; |
| | 662 | double py = cr*Y+sr*X; |
| | 663 | |
| | 664 | double dx = MeanX - px*1.02; |
| | 665 | double dy = MeanY - py*1.02; |
| | 666 | |
| | 667 | double norm = sqrt(dx*dx + dy*dy); |
| | 668 | double dist = norm*mm2deg; |
| | 669 | |
| | 670 | double lx = min(max((CosDelta*dy - SinDelta*dx)/norm, -1.), 1.); |
| | 671 | double ly = min(max((CosDelta*dx + SinDelta*dy)/norm, -1.), 1.); |
| | 672 | |
| | 673 | double alpha = asin(lx); |
| | 674 | double sgn = TMath::Sign(1., ly); |
| | 675 | |
| | 676 | // ------------------------------- Application ---------------------------------- |
| | 677 | |
| | 678 | double m3l = M3Long*sgn*mm2deg; |
| | 679 | double slope = SlopeLong*sgn/mm2deg; |
| | 680 | |
| | 681 | // --------------------------------- Analysis ----------------------------------- |
| | 682 | |
| | 683 | double xi = xi0 + slope0 *slope; |
| | 684 | |
| | 685 | double sign1 = m3l+0.07; |
| | 686 | double sign2 = (dist-0.5)*7.2-slope; |
| | 687 | |
| | 688 | double disp = (sign1<0 || sign2<0 ? -xi : xi)*(1-Width/Length); |
| | 689 | |
| | 690 | double thetasq = disp*disp + dist*dist - 2*disp*dist*sqrt(1-lx*lx); |
| | 691 | |
| | 692 | if (thetasq<0.024) |
| | 693 | cnt_on++; |
| | 694 | } |
| | 695 | |
| | 696 | hist.Fill(xi0, slope0, cnt_on); |
| | 697 | |
| | 698 | cout << xi0 << " " << slope0 << " " << cnt_on << endl; |
| | 699 | } |
| | 700 | |
| | 701 | TCanvas *c = new TCanvas; |
| | 702 | c->SetTopMargin(0.01); |
| | 703 | |
| | 704 | hist.SetStats(kFALSE); |
| | 705 | hist.SetContour(100); |
| | 706 | hist.SetXTitle("Xi_{0} / deg"); |
| | 707 | hist.SetYTitle("Slope_{0} #dot ns/deg^{2} "); |
| | 708 | hist.GetYaxis()->SetTitleOffset(1.2); |
| | 709 | hist.DrawCopy("colz cont2"); |
| | 710 | |
| | 711 | } |
| | 712 | }}} |
| | 713 | |
| | 714 | Results in |
| | 715 | |
| | 716 | [[Image(Result.png)]] |
| | 717 | |