| 722 | |
| 723 | == Optimize Disp (on Crab Data) == |
| 724 | |
| 725 | This is the same analysis using the Excess-Error on Crab Data. |
| 726 | |
| 727 | {{{#!Spoiler |
| 728 | {{{#!cpp |
| 729 | #include <iostream> |
| 730 | #include <iomanip> |
| 731 | |
| 732 | #include <TMath.h> |
| 733 | #include <TH2.h> |
| 734 | #include <TChain.h> |
| 735 | #include <TCanvas.h> |
| 736 | #include <TGraph.h> |
| 737 | #include <TStopwatch.h> |
| 738 | |
| 739 | Double_t LiMa(Double_t s, Double_t b, Double_t alpha=0.2) |
| 740 | { |
| 741 | const Double_t sum = s+b; |
| 742 | |
| 743 | if (s<0 || b<0 || alpha<=0) |
| 744 | return -1; |
| 745 | |
| 746 | const Double_t l = s==0 ? 0 : s*TMath::Log(s/sum*(alpha+1)/alpha); |
| 747 | const Double_t m = b==0 ? 0 : b*TMath::Log(b/sum*(alpha+1) ); |
| 748 | |
| 749 | return l+m<0 ? 0 : TMath::Sqrt((l+m)*2); |
| 750 | } |
| 751 | |
| 752 | Double_t Error(Double_t s, Double_t b, Double_t alpha=0.2) |
| 753 | { |
| 754 | const Double_t sN = s + alpha*alpha*b; |
| 755 | |
| 756 | const Double_t error = sN<0 ? 0 : TMath::Sqrt(sN); |
| 757 | if (error==0) |
| 758 | return 0; |
| 759 | |
| 760 | const Double_t Ns = s - alpha*b; |
| 761 | |
| 762 | return Ns<0 ? 0 : Ns/error; |
| 763 | } |
| 764 | |
| 765 | |
| 766 | double ganymed2(double xi0, double slope0) |
| 767 | { |
| 768 | // Create chain for the tree Result |
| 769 | // This is just easier than using TFile/TTree |
| 770 | TChain c("Result"); |
| 771 | |
| 772 | // Add the input file to the |
| 773 | c.AddFile("crab-data-only.root"); |
| 774 | |
| 775 | // Define variables for all leaves to be accessed |
| 776 | // By definition rootifysql writes only doubles |
| 777 | double FileId, EvtNumber, X, Y, MeanX, MeanY, Width, Length, CosDelta, SinDelta, |
| 778 | M3Long, SlopeLong, Leakage1, SlopeSpreadWeighted, Size, |
| 779 | ConcCore, ConcCOG, NumIslands, NumUsedPixels; |
| 780 | |
| 781 | // Connect the variables to the cordesponding leaves |
| 782 | c.SetBranchAddress("FileId", &FileId); |
| 783 | c.SetBranchAddress("EvtNumber", &EvtNumber); |
| 784 | c.SetBranchAddress("X", &X); |
| 785 | c.SetBranchAddress("Y", &Y); |
| 786 | c.SetBranchAddress("MeanX", &MeanX); |
| 787 | c.SetBranchAddress("MeanY", &MeanY); |
| 788 | c.SetBranchAddress("Width", &Width); |
| 789 | c.SetBranchAddress("Length", &Length); |
| 790 | c.SetBranchAddress("CosDelta", &CosDelta); |
| 791 | c.SetBranchAddress("SinDelta", &SinDelta); |
| 792 | c.SetBranchAddress("M3Long", &M3Long); |
| 793 | c.SetBranchAddress("SlopeLong", &SlopeLong); |
| 794 | c.SetBranchAddress("Leakage1", &Leakage1); |
| 795 | c.SetBranchAddress("NumIslands", &NumIslands); |
| 796 | c.SetBranchAddress("NumUsedPixels", &NumUsedPixels); |
| 797 | //c.SetBranchAddress("SlopeSpreadWeighted", &SlopeSpreadWeighted); |
| 798 | c.SetBranchAddress("Size", &Size); |
| 799 | //c.SetBranchAddress("ConcCOG", &ConcCOG); |
| 800 | //c.SetBranchAddress("ConcCore", &ConcCore); |
| 801 | |
| 802 | // Set some constants (they could be included in the database |
| 803 | // in the future) |
| 804 | double mm2deg = +0.0117193246260285378; |
| 805 | //double abberation = 1.02; |
| 806 | |
| 807 | // -------------------- Source dependent parameter calculation ------------------- |
| 808 | |
| 809 | // Create a histogram for on- and off-data |
| 810 | TH1F hon("on", "", 55, 0, 1); |
| 811 | TH1F hoff("off", "", 55, 0, 1); |
| 812 | |
| 813 | long cnt_on = 0; |
| 814 | long cnt_off = 0; |
| 815 | |
| 816 | for (int i=0; i<c.GetEntries(); i++) |
| 817 | { |
| 818 | // read the i-th event from the file |
| 819 | c.GetEntry(i); |
| 820 | |
| 821 | // First calculate all cuts to speedup the analysis |
| 822 | double area = TMath::Pi()*Width*Length; |
| 823 | |
| 824 | bool cutq = NumIslands<3.5 && NumUsedPixels>5.5 && Leakage1==0; |
| 825 | if (!cutq) |
| 826 | continue; |
| 827 | |
| 828 | bool cut0 = area < log10(Size)*898-1535; |
| 829 | if (!cut0) |
| 830 | continue; |
| 831 | |
| 832 | for (int angle=0; angle<360; angle+=60) |
| 833 | { |
| 834 | // -------------------- Source dependent parameter calculation ------------------- |
| 835 | |
| 836 | double cr = cos(angle*TMath::DegToRad()); |
| 837 | double sr = sin(angle*TMath::DegToRad()); |
| 838 | |
| 839 | double px = cr*X-sr*Y; |
| 840 | double py = cr*Y+sr*X; |
| 841 | |
| 842 | // blau rot |
| 843 | double dx = MeanX - px*1.02; |
| 844 | double dy = MeanY - py*1.02; |
| 845 | |
| 846 | double norm = sqrt(dx*dx + dy*dy); |
| 847 | double dist = norm*mm2deg; |
| 848 | |
| 849 | double lx = min(max((CosDelta*dy - SinDelta*dx)/norm, -1.), 1.); |
| 850 | double ly = min(max((CosDelta*dx + SinDelta*dy)/norm, -1.), 1.); |
| 851 | |
| 852 | double alpha = asin(lx); |
| 853 | double sgn = TMath::Sign(1., ly); |
| 854 | |
| 855 | // ------------------------------- Application ---------------------------------- |
| 856 | |
| 857 | double m3l = M3Long*sgn*mm2deg; |
| 858 | double slope = SlopeLong*sgn/mm2deg; |
| 859 | |
| 860 | // --------------------------------- Analysis ----------------------------------- |
| 861 | |
| 862 | double xi = xi0 + slope0 *slope;// + 1.67972 *(1-1/(1+4.86232*Leakage1)); |
| 863 | |
| 864 | double sign1 = m3l+0.07; |
| 865 | double sign2 = (dist-0.5)*7.2-slope; |
| 866 | |
| 867 | double disp = (sign1<0 || sign2<0 ? -xi : xi)*(1-Width/Length); |
| 868 | |
| 869 | //double thetasq = (disp*disp + dist*dist - 2*disp*dist*cos(alpha))*mm2deg*mm2deg; |
| 870 | double thetasq = disp*disp + dist*dist - 2*disp*dist*sqrt(1-lx*lx); |
| 871 | |
| 872 | // Fill the on- and off-histograms |
| 873 | if (thetasq<0.024) |
| 874 | angle==0 ? cnt_on++ : cnt_off++; |
| 875 | } |
| 876 | } |
| 877 | |
| 878 | cout << cnt_on << "\t" << cnt_off << "\t" << cnt_on-cnt_off*0.2 << endl; |
| 879 | return Error(cnt_on, cnt_off); |
| 880 | } |
| 881 | |
| 882 | |
| 883 | void ganymed_optimdisp() |
| 884 | { |
| 885 | TH2F hist("", "", |
| 886 | 21, 1.200-0.0050, 1.410-0.0050, |
| 887 | 13, 0.010-0.0025, 0.140-0.0025); |
| 888 | |
| 889 | for (float xi0=1.20; xi0<1.41; xi0+=0.01) |
| 890 | for (float slope0=0.01; slope0<0.14; slope0+=0.01) |
| 891 | { |
| 892 | cout << xi0 << "\t" << slope0 << "\t"; |
| 893 | hist.Fill(xi0, slope0, ganymed2(xi0, slope0)); |
| 894 | } |
| 895 | |
| 896 | TCanvas *canv = new TCanvas; |
| 897 | canv->SetTopMargin(0.01); |
| 898 | |
| 899 | hist.SetStats(kFALSE); |
| 900 | hist.SetContour(100); |
| 901 | hist.SetXTitle("Xi_{0} / deg"); |
| 902 | hist.SetYTitle("Slope_{0} #dot ns/deg^{2} "); |
| 903 | hist.GetYaxis()->SetTitleOffset(1.2); |
| 904 | hist.DrawCopy("colz cont2"); |
| 905 | } |
| 906 | }}} |
| 907 | }}} |
| 908 | |
| 909 | Which gives the following result |
| 910 | |
| 911 | [[Image(Result3.png)]] |
| 912 | |
| 913 | |