| 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 | |