Changeset 5313 for trunk/MagicSoft/Mars/mtemp/mmpi/MSkyPlot.cc

Timestamp:

10/25/04 17:11:26 (20 years ago)

Author:

mazin

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mtemp/mmpi/MSkyPlot.cc (modified) (33 diffs)

trunk/MagicSoft/Mars/mtemp/mmpi/MSkyPlot.cc

@@ -115 +115 @@
 #include "MSkyPlot.h"
 
+#include <vector>
 #include <TF1.h>
 #include <TH2.h>
@@ -190 +191 @@
     fHistOn.SetYTitle("y [\\circ]");
 
+    fHistSignifGaus.SetName("SignifDistrib");
+    fHistSignifGaus.SetTitle("Distribution of the significances from the sky plot");
+    fHistSignifGaus.SetXTitle("significance");
+    fHistSignifGaus.SetYTitle("counts");
+
     // set some values for the sky plot geometry:
     fMinXGrid    =-1.;                  // [deg] 
@@ -274 +280 @@
     fDistLo[7] =  0.0;
     }
+
+    fAlphaHName = "alphaplot.root";
+    fSkyHName   = "skyplot.root";
 }
 
 void MSkyPlot::ReadCuts(const TString parSCinit="OptSCParametersONOFFThetaRange0_1570mRad.root")
 {
+
+cout << " parameters read from file: " << parSCinit << endl;
     //--------------------------------
     // create container for the supercut parameters
@@ -390 +401 @@
     fRepDrive = (MReportDrive*)plist->FindObject(AddSerialNumber("MReportDrive"));
     if (!fRepDrive)
-        *fLog << warn << "MReportDrive not found... no sky map." << endl;
+        *fLog << warn << "MReportDrive not found... could be problem for sky map." << endl;
 
 
@@ -397 +408 @@
         *fLog << warn << "MSrcPosCam not found... no sky map." << endl;
 
+/*
     fPntPosCam = (MSrcPosCam*)plist->FindObject(AddSerialNumber("MPntPosCam"));
     if (!fPntPosCam)
         *fLog << warn << "MPntPosCam not found... no sky map." << endl;
+*/
 
     fPointPos = (MPointingPos*)plist->FindObject(AddSerialNumber("MPointingPos"));
@@ -407 +420 @@
     fTime = (MTime*)plist->FindObject(AddSerialNumber("MTime"));
     if (!fTime)
-        *fLog << warn << "MTime not found... no sky map." << endl;
+        *fLog << warn << "MTime not found... could be problem for sky map." << endl;
 
     fObservatory = (MObservatory*)plist->FindObject(AddSerialNumber("MObservatory"));
@@ -443 +456 @@
     fNumalphahist  =  (int) (fNumStepsX * fNumStepsY  + 0.5);
 
-    fHistSignif.SetName("SPSignif2ndOrder");
-    fHistSignif.SetTitle("Sky Plot of significance (2nd order fit)");
+   // fHistSignif.SetName("SPSignif2ndOrder");
+   // fHistSignif.SetTitle("Sky Plot of significance (2nd order fit)");
     fHistSignif.SetBins(fNumStepsX, fMinXGrid-0.5*fBinStepGrid, fMaxXGrid+0.5*fBinStepGrid, 
                         fNumStepsY, fMinYGrid-0.5*fBinStepGrid, fMaxYGrid+0.5*fBinStepGrid);
-    fHistSignif.SetDirectory(NULL);
+   // fHistSignif.SetDirectory(NULL);
     fHistSignif.SetFillStyle(4000);
     fHistSignif.UseCurrentStyle();
 
-    fHistNexcess.SetName("SPNex2ndOrder");
-    fHistNexcess.SetTitle("Sky Plot of Number of excess events (2nd order fit)");
+   // fHistNexcess.SetName("SPNex2ndOrder");
+   // fHistNexcess.SetTitle("Sky Plot of Number of excess events (2nd order fit)");
     fHistNexcess.SetBins(fNumStepsX, fMinXGrid-0.5*fBinStepGrid, fMaxXGrid+0.5*fBinStepGrid, 
                         fNumStepsY, fMinYGrid-0.5*fBinStepGrid, fMaxYGrid+0.5*fBinStepGrid);
-    fHistNexcess.SetDirectory(NULL);
+   // fHistNexcess.SetDirectory(NULL);
     fHistNexcess.SetFillStyle(4000);
     fHistNexcess.UseCurrentStyle();
 
-    fHistOn.SetName("SPOnCounted");
-    fHistOn.SetTitle("Sky Plot of ON events (counted)");
+   // fHistOn.SetName("SPOnCounted");
+   // fHistOn.SetTitle("Sky Plot of ON events (counted)");
     fHistOn.SetBins(fNumStepsX, fMinXGrid-0.5*fBinStepGrid, fMaxXGrid+0.5*fBinStepGrid, 
                         fNumStepsY, fMinYGrid-0.5*fBinStepGrid, fMaxYGrid+0.5*fBinStepGrid);
-    fHistOn.SetDirectory(NULL);
+   // fHistOn.SetDirectory(NULL);
     fHistOn.SetFillStyle(4000);
     fHistOn.UseCurrentStyle();
 
+    //fHistSignifGaus.SetName("SignifDistrib");
+    fHistSignifGaus.SetTitle("Distribution of the significances from the sky plot");
+    fHistSignifGaus.SetBins(100, -10., 10.); 
 
     // prepare alpha plots
         // temp histogram for an alpha plot
-    TH1I *temp = new TH1I("alphaplot","alphaplot",fNumBinsAlpha,fAlphaLeftEdge,fAlphaRightEdge);
+    TH1D *temp = new TH1D("alphaplot","alphaplot",fNumBinsAlpha,fAlphaLeftEdge,fAlphaRightEdge);
     temp->SetDirectory(NULL);
     fHistAlphaBinWidth = temp->GetBinWidth(1);
@@ -499 +515 @@
   if (fSetCenter==kTRUE)
   {
-     fRa0  = fRepDrive->GetRa();  // [h]
-     fDec0 = fRepDrive->GetDec(); // [deg]
+     if (fRepDrive)
+     {
+         fRa0  = fRepDrive->GetRa();  // [h]
+         fDec0 = fRepDrive->GetDec(); // [deg]
+         if (fRa0 < 0. || fRa0 > 24. || fDec0 < -90. || fDec0 > 90. || (fRa0==0 && fDec0==0)) return kTRUE;  // temp!, should be changed
+     }
+     else
+     {
+        fRa0 = fPointPos->GetRa(); 
+        fDec0 = fPointPos->GetDec();
+     } 
+     *fLog << inf << "Ra (center) = " << fRa0 << ", Dec = " << fDec0 << endl;
      fSetCenter=kFALSE;
   }
@@ -518 +544 @@
     Double_t rho = 0;
     if (fTime && fObservatory && fPointPos)
-{
+    {
         rho = fPointPos->RotationAngle(*fObservatory, *fTime);
-*fLog << inf << " rho = " << rho << ", Zd = " << fPointPos->GetZd() << 
-                ", Az = " << fPointPos->GetAz() << ", Ra = " << fPointPos->GetRa() << ", Dec = " << fPointPos->GetDec() << endl;
-}
-    //if (fPointPos)
-    //    rho = fPointPos->RotationAngle(*fObservatory);
+//*fLog << inf << " rho = " << rho*180./TMath::Pi() << ", Zd = " << fPointPos->GetZd() << 
+//                ", Az = " << fPointPos->GetAz() << ", Ra = " << fPointPos->GetRa() << ", Dec = " << fPointPos->GetDec() << endl;
+
+    // => coord. system: xtilde, ytilde with center in RA_0, DEC_0 
+
+// Get Rot. Angle:
+    cosgam = TMath::Cos(rho);
+    singam = TMath::Sin(rho);
+// Get x_0, y_0 for RA_0,DEC_0 of the current event
+    }
+    else 
+    {
+//      rho = fPointPos->RotationAngle(*fObservatory);
+      Double_t theta, phi, sing, cosg;
+      theta = fPointPos->GetZd()*TMath::Pi()/180.;
+      phi = fPointPos->GetAz()*TMath::Pi()/180.;
+//   printf("theta: %5.3f, phi: %5.3f\n", theta*180./4.1415, phi*180./4.1415);
+      fObservatory->RotationAngle(theta, phi, sing, cosg);        
+      cosgam = cosg;
+      singam = sing;
+    }
+   // if (fPointPos)
+   //    rho = fPointPos->RotationAngle(*fObservatory);
+
 /*
 //TEMP
@@ -556 +601 @@
 */
 
-    // => coord. system: xtilde, ytilde with center in RA_0, DEC_0 
-
-// Get Rot. Angle:
-        cosgam = TMath::Cos(rho);
-        singam = TMath::Sin(rho);
-// Get x_0, y_0 for RA_0,DEC_0 of the current event
-        x_0 = fPntPosCam->GetX()*fMm2Deg; 
-        y_0 = fPntPosCam->GetY()*fMm2Deg; 
+/*
+    x_0 = fPntPosCam->GetX()*fMm2Deg; 
+    y_0 = fPntPosCam->GetY()*fMm2Deg; 
+*/      
+    x_0 = 0.;
+    y_0 = 0.;
 
     Int_t index = 0;  // index for alpha histograms
@@ -598 +641 @@
            lgsize = logsize-log3000;
            lgsize2 = lgsize*lgsize;
-           dist2 = (distpar-fDistOffset)*(distpar-fDistOffset);
+           dist2 = distpar*distpar - fDistOffset*fDistOffset;
 
            if ( (fHillas->GetLength()*fMm2Deg) < CalcLimit(fLengthUp, lgsize, lgsize2, dist2) &&
@@ -609 +652 @@
            {
 // gamma candidates!
+//*fLog <<  "Length : " << fHillas->GetLength()*fMm2Deg << ", Width : " << fHillas->GetWidth()*fMm2Deg << endl;
+//*fLog <<  "distpar: " << distpar << ", Size : " << fHillas->GetSize() << endl;
                 fHistAlpha[index].Fill(alphapar);
-                index++;
            }
+           index++;
         }
      }
@@ -628 +673 @@
     const Int_t onbinsalpha = TMath::Nint(fAlphaONMax/fHistAlphaBinWidth);
 
+// fit with gaus 
+    fHistSignifGaus.Fit("gaus");
+
 // fit function for the background
     TF1 * fitbgpar = new TF1("fbgpar", "[0]*x*x + [1]", fAlphaBgLow, fAlphaBgUp);
@@ -635 +683 @@
     numdegfreed = (fAlphaBgUp - fAlphaBgLow) / fHistAlphaBinWidth - 2.; 
 
-    vector <TH1I>::iterator alpha_iterator;
+    vector <TH1D>::iterator alpha_iterator;
     int index2 = 0;  // index of the TH2F histograms
 
@@ -643 +691 @@
          // find the bin in the 2dim histogram
          nrow    = index2/fHistOn.GetNbinsX() + 1;                    // row of the histogram (y)
-         ncolumn = TMath::Nint(fmod(index2,fHistOn.GetNbinsX()))+1;   // column of the histogram (x)
+         ncolumn = index2%fHistOn.GetNbinsX()+1;   // column of the histogram (x)
+         //ncolumn = TMath::Nint(fmod(index2,fHistOn.GetNbinsX()))+1;   // column of the histogram (x)
 
          // number of ON events below fAlphaONMax
@@ -659 +708 @@
 
          // estimate Noff from the fit:
-         Noff_fit = (1./3. * (fitbgpar->GetParameter(0)) * pow(fAlphaONMax,3.) +
+         Noff_fit = (1./3. * (fitbgpar->GetParameter(0)) * TMath::Power(fAlphaONMax,3.) +
                     (fitbgpar->GetParameter(1)) * fAlphaONMax) /  fHistAlphaBinWidth;
 
@@ -671 +720 @@
 
          fHistSignif.SetBinContent(ncolumn, nrow, Sign);        // fill in the fHistSignif
+         fHistSignifGaus.Fill(Sign);
 
          alpha_iterator++;
@@ -677 +727 @@
 
 //temp
+/*
    Double_t decl, hang;
    MStarCamTrans mstarc(*fGeomCam, *fObservatory);
@@ -683 +734 @@
 *fLog << warn << "MDrive, RA, DEC = " << fRepDrive->GetRa() << ", " << fRepDrive->GetDec() << endl;
 *fLog << warn << "MStarCamTrans, H angle , DEC = " << hang << ", " << decl << endl;
+*/
 //endtemp
 
+
 // save alpha plots:
-  if(kSaveAlphaPlots==kTRUE) SaveAlphaPlots();
+//  TString stri1 = "alphaplots.root";
+  if(kSaveAlphaPlots==kTRUE) SaveAlphaPlots(fAlphaHName);
 
 // save sky plots:
-  if(kSaveSkyPlots==kTRUE) SaveSkyPlots();
+//  TString stri2 = "skyplots.root";
+  if(kSaveSkyPlots==kTRUE) SaveSkyPlots(fSkyHName);
 
   fHistAlpha.clear();
@@ -721 +776 @@
 
 
-void MSkyPlot::SaveSkyPlots(const TString skyplotfilename="skyplots.root")
+void MSkyPlot::SaveSkyPlots(TString skyplotfilename)
 {
    TFile rootfile(skyplotfilename, "RECREATE",
@@ -728 +783 @@
      fHistNexcess.Write();
      fHistOn.Write();
+     fHistSignif.Write();
 
 // from Wolfgang: 
@@ -733 +789 @@
     // Dec-Hour lines
     Double_t rho, sinrho, cosrho;
-//    if (fTime && fObservatory && fPointPos)
+    Double_t theta, phi, sing, cosg;
+// do I need it?
+    if (fTime && fObservatory && fPointPos)
+    {
         rho = fPointPos->RotationAngle(*fObservatory, *fTime);
-    sinrho=TMath::Sin(rho);
-    cosrho=TMath::Cos(rho);
-
-    Double_t theta, phi, sing, cosg;
-    theta = fPointPos->GetZd()*TMath::Pi()/180.;
-    phi = fPointPos->GetAz()*TMath::Pi()/180.;
+        sinrho=TMath::Sin(rho);
+        cosrho=TMath::Cos(rho);
+    }
+
+       theta = fPointPos->GetZd()*TMath::Pi()/180.;
+       phi = fPointPos->GetAz()*TMath::Pi()/180.;
 //   printf("theta: %5.3f, phi: %5.3f\n", theta*180./4.1415, phi*180./4.1415);
-    fObservatory->RotationAngle(theta, phi, sing, cosg);
+       fObservatory->RotationAngle(theta, phi, sing, cosg);
 
 *fLog << "1: sinrho, cosrho = " << sinrho << ", " << cosrho << endl;
 *fLog << "2: sinrho, cosrho = " << sing << ", " << cosg << endl;
+       sinrho=sing;
+       cosrho=cosg;
 
     Double_t fDistCam = fGeomCam->GetCameraDist() * 1000.0;     //  [mm]
@@ -756 +817 @@
     Double_t declin, hangle;  // [deg], [h]
     MStarCamTrans mstarc(*fGeomCam, *fObservatory);
-    mstarc.LocToCel(fRepDrive->GetNominalZd(),fRepDrive->GetNominalAz(),declin, hangle);
-
-    TGraph *graph1;
+    if (fRepDrive) mstarc.LocToCel(fRepDrive->GetNominalZd(),fRepDrive->GetNominalAz(),declin, hangle);
+    else mstarc.LocToCel(theta*180./TMath::Pi(),phi*180./TMath::Pi(),declin, hangle); // NOT GOOD!
+
+//    std::vector <TGraph> graphvec;
     TLatex *pix;
 
@@ -776 +838 @@
     // direction for the center of the camera
     Double_t dec0 = declin;
-    Double_t h0   = hangle*360./24.;
+    Double_t h0   = hangle*360./24.; //deg
+    Double_t RaHOffset = fRepDrive->GetRa() - h0;
     //trans.LocToCel(theta0, phi0, dec0, h0);
 
+    gStyle->SetOptFit(0);
+    gStyle->SetOptStat(0);
+    gStyle->SetPalette(1);
     TCanvas *c1 = new TCanvas("SPlotsRaDecLines","SPlotsRaDecLines", 400, 0, 700, 600);
+    c1->Divide(2,2);
     c1->cd(1);
+    fHistSignif.Draw("colz");
+    c1->cd(2);
+    fHistNexcess.Draw("colz");
+    c1->cd(3);
     fHistOn.Draw("colz");
+    c1->cd(4);
+    gPad->SetLogy();
+    fHistSignifGaus.Draw(); 
 
     //-----   lines for fixed dec   ------------------------------------
@@ -796 +870 @@
       dh = 360.0/((Double_t)(Nh-1));
 
-
+// start to copy
     for (Int_t j=0; j<Ndec; j++)
     {
@@ -810 +884 @@
     Double_t xh[Nh];
     Double_t yh[Nh];
-
-
 
     for (Int_t j=0; j<Ndec; j++)
@@ -839 +911 @@
 
 //      c1->cd(2);
-      graph1 = new TGraph(Nh,xh,yh);
-      graph1->SetLineColor(j+1);
+      TGraph * graph1 = new TGraph(Nh,xh,yh);
+      //graph1->SetLineColor(j+1);
+      graph1->SetLineColor(36);
       graph1->SetLineStyle(1);
-      graph1->SetMarkerColor(j+1);
+      graph1->SetLineWidth(1);
+      //graph1->SetMarkerColor(j+1);
+      graph1->SetMarkerColor(1);
       graph1->SetMarkerSize(.2);
       graph1->SetMarkerStyle(20);
-      graph1->Draw("PL");
+      
+      c1->cd(1);
+      graph1->Draw("L");
+      c1->cd(2);
+      graph1->Draw("L");
+      c1->cd(3);
+      graph1->Draw("L");
       //delete graph1;
+//      graphvec.push_back(*graph1);
+//      graphvec[j].Draw("L");
 
       sprintf(tit,"Dec = %6.2f", dec[j]);
@@ -852 +935 @@
       ytxt = ylo + (yup-ylo)*0.80 - ((Double_t)j) *(yup-ylo)/20.0;
       pix = new TLatex(xtxt, ytxt, tit);
-      pix->SetTextColor(j+1);
+      pix->SetTextColor(36);
       pix->SetTextSize(.03);
-      pix->Draw("");
+      //pix->Draw("");
       //delete pix;
 
     }
-
+//stop copy
     //-----   lines for fixed hour angle   ------------------------------------
 
@@ -872 +955 @@
       dh1 = 360.0/((Double_t)Nh1);
 
-
+// start copy
     for (Int_t j=0; j<Ndec1; j++)
     {
@@ -915 +998 @@
 
 //      c1->cd(2);
-      graph1 = new TGraph(count,xd,yd);
-      graph1->SetLineColor(k+1);
-      graph1->SetLineStyle(2);
-      graph1->SetMarkerColor(k+1);
+      TGraph * graph1 = new TGraph(count,xd,yd);
+      //graph1->SetLineColor(k+1);
+      graph1->SetLineColor(36);
+      graph1->SetLineWidth(2);
+      graph1->SetLineStyle(1);
+      //graph1->SetMarkerColor(k+1);
+      graph1->SetMarkerColor(1);
       graph1->SetMarkerSize(.2);
       graph1->SetMarkerStyle(20);
-      graph1->Draw("PL");
-      //delete graph1;
-
-      sprintf(tit,"H = %6.2f", h1[k]);
+      c1->cd(1);
+      graph1->Draw("L");
+      c1->cd(2);
+      graph1->Draw("L");
+      c1->cd(3);
+      graph1->Draw("L");
+
+      sprintf(tit,"RA = %6.2f", fRa0 + (h0 - h1[k]));
       Double_t xtxt = xlo + (xup-xlo)*0.75;
       Double_t ytxt = ylo + (yup-ylo)*0.80 - ((Double_t)k) *(yup-ylo)/20.0;
       pix = new TLatex(xtxt, ytxt, tit);
-      pix->SetTextColor(k+1);
+      pix->SetTextColor(36);
       pix->SetTextSize(.03);
-      pix->Draw("");
+      //pix->Draw("");
       //delete pix;
 
@@ -936 +1026 @@
 
 //    c1->cd(2);
-    sprintf(tit,"Dec0 = %6.2f    H0 = %6.2f", dec0, h0);
-    xtxt = xlo + (xup-xlo)*0.05;
-    ytxt = ylo + (yup-ylo)*0.92;
+    sprintf(tit,"Dec0 = %6.2f [deg]   Ra0 = %6.2f [h]", dec0, fRa0);
+    xtxt = xlo + (xup-xlo)*0.05 + 0.80;
+    ytxt = ylo + (yup-ylo)*0.75;
     pix = new TLatex(xtxt, ytxt, tit);
     pix->SetTextColor(1);
-    pix->SetTextSize(.06);
+    pix->SetTextSize(.05);
+    c1->cd(1);
+    pix->Draw("");
+    c1->cd(2);
+    pix->Draw("");
+    c1->cd(3);
     pix->Draw("");
     //delete pix;
 
-    sprintf(tit,"              [deg]             [deg]");
-    xtxt = xlo + (xup-xlo)*0.05;
-    ytxt = ylo + (yup-ylo)*0.86;
-    pix = new TLatex(xtxt, ytxt, tit);
-    pix->SetTextColor(1);
-    pix->SetTextSize(.06);
-    pix->Draw("");
-
     c1->Write();
+    // we suppose that the {skyplotfilename} ends with .root
+    Int_t sizeofout = skyplotfilename.Sizeof();
+    TString outps = skyplotfilename.Remove(sizeofout-5,5) + "ps";
+    c1->Print(outps);  // temporary!!!!!
 
     TCanvas *c2 = new TCanvas("SkyPlotsWithRaDecLines","SkyPlotsWithRaDecLines", 0, 0, 300, 600);
@@ -967 +1058 @@
    delete c1;
    delete c2;
-   delete graph1;
    delete pix;
+
+
 }
 
-void MSkyPlot::SaveAlphaPlots(const TString alphaplotfilename="alphaplots.root")
+void MSkyPlot::SaveAlphaPlots(const TString alphaplotfilename)
 {
    TFile rootfile(alphaplotfilename, "RECREATE",
                    "all the alpha plots");
 
-    vector <TH1I>::iterator alpha_iterator;
+    vector <TH1D>::iterator alpha_iterator;
     int index = 0;  // index of the TH2F histograms
     Char_t strtitle[100];
@@ -989 +1081 @@
 
           nrow    = index/fHistOn.GetNbinsX() + 1;                    // row of the histogram (y)
-          ncolumn = TMath::Nint(fmod(index,fHistOn.GetNbinsX()))+1;   // column of the histogram (x)
+          //ncolumn = TMath::Nint(fmod(index,fHistOn.GetNbinsX()))+1;   // column of the histogram (x)
+          ncolumn = index%fHistOn.GetNbinsX()+1;   // column of the histogram (x)
           xpos    = fMinXGrid + fBinStepGrid*(ncolumn-1);
           ypos    = fMinYGrid + fBinStepGrid*(nrow-1);