/* ======================================================================== *\
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cherenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!   Author(s): Daniel Mazin, 05/2004 <mailto:mazin@imppmu.mpg.de>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

// Daniel Mazin 17.05.2004  mazin@mppmu.mpg.de
// **********************************************************************************
// this macro is used to produce false source plots
// and it is a part from the macro falsesourcemethod.C.
// See falsesourcemethod.C for more explanation
// input: hillas parameters as vectors
// output: several 2D plots in a root file +
// ascii file with parameters like significance, Nexcess etc.
// **********************************************************************************


void gridloop()
{

    printf("inside of gridloop() ... \n");

/* emty histograms  */

    histskyplot.Reset();
    histNexOnOff.Reset();
    histChi2Off.Reset();
    histskyplotOnOff.Reset();
    histNex.Reset();
    histChi2.Reset();
    histskyplot4.Reset();
    histskyLiMa.Reset();
    histsign.Reset();
    histsign4.Reset();
    histBerlin.Reset();

/* end emty histograms  */


    TH1F histalphagrid;
    histalphagrid.SetName("Alpha");
    histalphagrid.SetTitle("Alpha");
    histalphagrid.SetXTitle("Alpha (deg)");
    histalphagrid.SetYTitle("Counts");
    histalphagrid.SetDirectory(NULL);
    histalphagrid.SetFillStyle(4000);
    histalphagrid.UseCurrentStyle();
  
    bins.SetEdges(18, 0.0, 90.);
    bins.Apply(histalphagrid);

    Double_t xcorr, ycorr, xtest, ytest, xtilde, ytilde, fconc1;

    const Float_t binwidth1 = histalphagrid.GetBinWidth(1);
    if (TMath::Abs( ALPHAMAX/binwidth1 - TMath::Nint(ALPHAMAX/binwidth1)) > TOLERANCE)
    {
         cout << "wrong binning for histogram: histalphagrid \n SORRY ... ABORT" << endl;
         exit(1);
    }
    const Int_t numbbinOnMax = TMath::Nint(ALPHAMAX / binwidth1);  // number bins in ON region

  /*   fit functions  */
    TF1 * fitsig = new TF1("fsig", "[0]", 0., ALPHAMAX);
    fitsig->SetLineColor(4); 
    TF1 * fitbgpar = new TF1("fbgpar", "[0]*x*x + [1]", 30., 90.);
    fitbgpar->SetLineColor(2); 

    TF1 * fitbg4 = new TF1("fbg4", "[0]*x*x*x*x + [1]*x*x + [2]", 30., 90.);
    fitbg4->SetLineColor(3);
  /*   end fit functions  */

    Float_t Nex, Non, Noff, Sign;
    Float_t xsource, ysource, xsournew, ysournew;

    Double_t Nex4, Noff4, Sign4; 
    Double_t alpha, SignLiMa;
 
    Double_t chisquarefit;

    Double_t integon, integoff, normf, NoffOFF, SignOnOff, NexOnOff;
    Float_t logsize, lgsize, lgsize2, tanbeta, beta;
    const Float_t LOG3000 = log(3000.);



    Double_t chi2par[GRIDBINS], apar[GRIDBINS], bpar[GRIDBINS]; // fit parameters from OFF sample
    Float_t offbin1[GRIDBINS], offbin2[GRIDBINS], offbin3[GRIDBINS]; // Number of events in first 3 bins
    Double_t xpos, ypos, binwidth2;
    Int_t  numgrid;


    Int_t numbinsoff = 12;  // FIXME: 12 should be changed according to the binning
    Int_t nbinx, nbiny, gridpoint = 0;

    Char_t psname[120], titelname[120], plot1name[120], plot2name[120], stringsig[120];
 
    FILE * fp;

    if (TYPEOPTION == kTRUE)  // read fit parameters from the OFF sample
    {
      Char_t * datin = DATPARAMIN;
              
      fp = fopen(datin, "r");

      for (Int_t i=0; i<GRIDBINS; i++)
      {
        fscanf(fp,"%d %lf %lf %lf %lf %lf %lf %f %f %f",
                 &numgrid, &xpos, &ypos,  &apar[i], &bpar[i], 
                 &chi2par[i], &binwidth2,
                 &offbin1[i], &offbin2[i], &offbin3[i]);
//  cout << "numgrid: " << numgrid << ",  apar " << apar[i] << endl;
      }

       fclose(fp);
    }


    Char_t * chout = DATPARAMOUT;  // write parameters into file

    fp = fopen(chout, "w");


/* begin of the loop over all bins in the sky map */

   for(Int_t gridx = 0; gridx < NUMSTEPS; gridx++)   
   { 
    xsource = MINXGRID + gridx * STEPGRID;
//cout << xsource << endl;

    for(Int_t gridy = 0; gridy < NUMSTEPS; gridy++)  
    {
      ysource = MINYGRID + gridy * STEPGRID;
//cout << ysource << endl;
    
// empty histogram:
//     nbinx = histalphagrid.GetNbinsX();
//     for(Int_t ii = 0; ii< (nbinx+1); ii++) histalphagrid.SetBinContent(ii, 0.);
//     histalphagrid.SetEntries(0);
     histalphagrid.Reset();


//      Bool_t anglephi = kFALSE;

    /* loop over all events */
      for(Int_t k=0; k< imnum; k++)
      {

          if (ROTOPTION == kTRUE) // derotate into sky coordinates
          {
      /*     derotation    : correct sky coordinates into camera coordinates */
             xsournew = cosgam(k) * xsource - singam(k) * ysource;
             ysournew = singam(k) * xsource + cosgam(k) * ysource;
       /*    end derotatiom    */
          }
          else // do not derotate, plot into camera coordinates
          {
             xsournew = xsource;
             ysournew = ysource;
          }


       /* calculate ALPHA und DIST acording to the source position */
           tanbeta = (ypar(k) - ysournew) / (xpar(k) - xsournew);
           beta = TMath::ATan(tanbeta);
           alphapar(k) = (deltapar(k) - beta) * 180./ TMath::Pi();
           distpar(k) = sqrt((ypar(k) - ysournew) * (ypar(k) - ysournew) +
                             (xpar(k) - xsournew) * (xpar(k) - xsournew));



          if(alphapar(k) > 90.) alphapar(k) -= 180.;
          if(alphapar(k) < -90.) alphapar(k) += 180.;

//          histalpha.Fill(alphapar(k),1.);
//          histdist.Fill( distpar(k),1.);

          if(CUTOPTION==kFALSE)
          { 
/* ******************************************************   */
      /*   static   cuts   */
                if (lengthpar(k) > LENGTHMIN && lengthpar(k) < LENGTHMAX)
                if (widthpar(k) > WIDTHMIN && widthpar(k) < WIDTHMAX)
                if (sizepar(k) > SIZEMIN)
                if (distpar(k) > DISTMIN && distpar(k) < DISTMAX)
                if (leak1(k) < LEAKMAX)
/* ******************************************************   */
                     {
                         alphapar(k) = TMath::Abs(alphapar(k));
//printf("alpha passed the cuts: %5.2f \n", alphapar(k));
//printf("histalphagrid:  %5.2f \n", histalphagrid.GetNbinsX());
                         histalphagrid.Fill(alphapar(k),1.);

                     }

           }
/* ******************************************************   */
/*            dynamical cuts                                */
           else     // if(CUTOPTION==kTRUE)
           {
                 logsize = log(sizepar(k));
                 lgsize = logsize-LOG3000;
                 lgsize2 = lgsize*lgsize;
                 if ( sizepar(k) > SIZEMIN )
                 if ( leak1(k) < LEAKMAX )
                 if ( lengthpar(k) > (LENGTHMINParA + LENGTHMINParB*lgsize + LENGTHMINParC*lgsize2) &&
                      lengthpar(k) < (LENGTHMAXParA + LENGTHMAXParB*lgsize + LENGTHMAXParC*lgsize2))
                 if ( widthpar(k) > (WIDTHMINParA + WIDTHMINParB*lgsize + WIDTHMINParC*lgsize2)  &&
                      widthpar(k) < (WIDTHMAXParA + WIDTHMAXParB*lgsize + WIDTHMAXParC*lgsize2) )
                 if ( distpar(k) > (DISTMINParA + DISTMINParB*lgsize + DISTMINParC*lgsize2)  &&
                      distpar(k) < (DISTMAXParA + DISTMAXParB*lgsize + DISTMAXParC*lgsize2) )
                 {
                         alphapar(k) = TMath::Abs(alphapar(k));
                         histalphagrid.Fill(alphapar(k),1.);
                         if(alphapar(k) < 8.)   histBerlin.Fill(xsournew, ysournew);
                 }
           }
         
        }
//    histalphagrid.Fit("fsig","NR");
//printf("content bin 1: %4.2f\n", histalphagrid.GetBinContent(1));
    histalphagrid.Fit(fitsig,"NWR");
    histalphagrid.Fit("fbg4","+NWR");
    histalphagrid.Fit("fbgpar","+NWR");

// cout << "OK" << endl;
// calculate significance:
    Non = 0.;
    for(Int_t i=1; i<=numbbinOnMax;i++) Non += histalphagrid.GetBinContent(i);
//    Non = (fitsig->GetParameter(0)) * numbbinOnMax;
//cout << "Non : " << Non << endl;
//    histalphagrid.Draw();

    chisquarefit = fitbgpar->GetChisquare(); 
    Noff = (1./3. * (fitbgpar->GetParameter(0)) * pow(ALPHAMAX,3.) + 
           (fitbgpar->GetParameter(1)) * ALPHAMAX) /  binwidth1;
    
//cout << "Noff : " << Noff << endl;
    
    Nex = Non - Noff;

    Noff4 = (1./5. * fitbg4->GetParameter(0) *  pow(ALPHAMAX,5.) +
            1./3. * (fitbg4->GetParameter(1)) * pow(ALPHAMAX,3.) +
           fitbg4->GetParameter(2) * ALPHAMAX ) / binwidth1;
 
//cout << "Noff4 : " << Noff4 << endl;
 
    Nex4 = Non - Noff4;
 
// calculate significance: 

//    Sign4 = Nex4 / sqrt(Nex4 + 2.* Noff4);
    if (Noff4<0.) Sign4 = 0.;
    else Sign4 = LiMa17(Non,Noff4,1.);

//cout << "Sign4 : " << Sign4 << endl;
//    Sign = Nex / sqrt(Nex + 2.* Noff);
    if (Noff<0.) Sign = 0.;
    else Sign = LiMa17(Non,Noff,1.);
//cout << "Sign LiMa : " << Sign << endl;

   /*  calculate Noff from the OFF data  */
  if (TYPEOPTION == kTRUE)
  {
   	// ON: 

   	integon = 0.;  // number of events between 30 and 90 degrees

	   for (Int_t ik = 0; ik < numbinsoff; ik++)
	   {
	    integon += histalphagrid.GetBinContent(7+ik);  // FIXME: 7 should be changed according to the binning
	   }
	   // OFF: 

	   integoff = ((1./3. * apar[gridpoint] * pow(90.,3.) + bpar[gridpoint] * 90.) -
        	      (1./3. * apar[gridpoint] * pow(30.,3.) + bpar[gridpoint] * 30.)) / binwidth2;
	
	   normf = integoff / integon;

	   NoffOFF = (1./3. * apar[gridpoint] * pow(ALPHAMAX,3.) +
        	   (bpar[gridpoint] * ALPHAMAX)) /  binwidth2 / normf;

	   NexOnOff = Non - NoffOFF;

  	   SignOnOff = NexOnOff / sqrt(NexOnOff + 2.* NoffOFF);

	// calculate according to Li Ma:
           if (NoffOFF<0.) SignLiMa=0.;
	   else SignLiMa = LiMa17(Non,NoffOFF*normf,1./normf);

   }

    cout << endl << "   xsource: " << xsource <<  "   ysource: " << ysource << endl;

    cout << " Non : " << Non << " Noff : " << Noff << " Nex : " << Nex << endl;

 if (TYPEOPTION == kTRUE) 
    cout << " ON and OFF : Non : " << Non << " Noff (estimated) : " << NoffOFF << " Nex : " << NexOnOff << endl;

    cout << endl;
    cout << " significance (parab): " << Sign << " sigma" << endl;
    cout << " significance (4th order): " << Sign4 << " sigma" << endl;
 if (TYPEOPTION == kTRUE)
  {
    cout << " significance (using ON and OFF): " << SignOnOff << " sigma" << endl;
    cout << " significance (LiMa 17): " << SignLiMa << " sigma" << endl;
  }


    if (TYPEOPTION == kTRUE)
      fprintf(fp,"%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf \n", 
                xsource, ysource, Sign, Non, Noff, Nex, SignLiMa, 
                SignOnOff, NoffOFF, NexOnOff, chi2par[gridpoint]);
    else 
         fprintf(fp,"%d %lf %lf %lf %lf %lf %lf %lf %lf %lf\n",
            gridpoint+1, xsource, ysource,
            fitbgpar->GetParameter(0), fitbgpar->GetParameter(1),
            fitbgpar->GetChisquare(), histalphagrid.GetBinWidth(1),
            histalphagrid.GetBinContent(1), histalphagrid.GetBinContent(2),
            histalphagrid.GetBinContent(3));


    histskyplot.SetBinContent(gridx+1, gridy+1, Sign);
    histNex.SetBinContent(gridx+1, gridy+1, Nex);
    histChi2.SetBinContent(gridx+1, gridy+1, fitbgpar->GetChisquare() / 10.);

    if (TYPEOPTION == kTRUE) histskyplotOnOff.SetBinContent(gridx+1, gridy+1, SignOnOff);
    if (TYPEOPTION == kTRUE) histNexOnOff.SetBinContent(gridx+1, gridy+1, NexOnOff);
    if (TYPEOPTION == kTRUE) histChi2Off.SetBinContent(gridx+1, gridy+1, chi2par[gridpoint] / 16.);

    if (TYPEOPTION == kTRUE) histskyLiMa.SetBinContent(gridx+1, gridy+1, SignLiMa);
    histskyplot4.SetBinContent(gridx+1, gridy+1, Sign4);
    histsign.Fill(Sign,1.);
    histsign4.Fill(Sign4,1.);
//    break;

    if (TYPEOPTION == kTRUE)
    {
	    TF1 * bgoff = new TF1("bgoff", parabfunc, 0., 90., 3);
	    bgoff->SetParameters(apar[gridpoint], bpar[gridpoint], normf);
	    bgoff->FixParameter(0, apar[gridpoint]);
	    bgoff->FixParameter(1, bpar[gridpoint]);
	    bgoff->FixParameter(2, normf);
	    bgoff->SetLineColor(8);
    }

//    sprintf(stringsig,"S = %.2f sigma", Sign);
    if (TYPEOPTION == kTRUE) sprintf(stringsig,"S = %.2f sigma", SignOnOff);
/*
    sprintf(psname,"/.magic/magicserv01/scratch/Daniel/plots/RootPlots/040215Mrk421_misp/all/alphaX%.2fY%.2fON.root", xsource, ysource);
    sprintf(titelname,"alpha plot. assumed source position: x = %.2f y = %.2f", xsource, ysource);
    TCanvas * cplot  = new TCanvas();
    histalphagrid.SetTitle(titelname);
    histalphagrid.Draw();
    bgoff->Draw("same");
    leg = new TLegend(0.1,0.15,0.52,0.35);
    leg->Draw();
    leg->AddEntry(fitsig,"fit for ON region (0-20 deg)","l");
    leg->AddEntry(fitbgpar,"fit for OFF region (30-90 deg)","l");
    leg->AddEntry(bgoff,"fit from OFF data (0-90 deg)","l");
    leg->SetHeader("Legend");
    leg->SetFillColor(19);
    leg->Draw();

    text = new TPaveText(50,20,90,80);
    text->AddText(0.5, 0.8, "calculated significance");
    text->AddText(0.5, 0.5, "using ON and OFF:");
    text->AddText(1, 0.2, stringsig);
    text->SetTextSize(0.04);
    text->Draw();

    cplot->Modified();
    cplot->Update();

    cplot->SaveAs(psname);
    delete cplot;
*/
    if (TYPEOPTION == kTRUE) delete  bgoff;

    gridpoint++;

    }
   }

  /* close file  */
    fclose(fp);  
    

    TFile rootfile(ROOTPLOTNAME, "RECREATE",
                   "sky plots in all variations");


//    cout << " end of one of the grids 1 " << endl;

     histskyplot.Write();

     histNex.Write();

     histChi2.Write();

     histsign.Write();

     histskyplot4.Write();

     histsign4.Write();

     histBerlin.Write();

   if (TYPEOPTION == kTRUE)
   {
     histNexOnOff.Write();
     histChi2Off.Write();
     histskyplotOnOff.Write();
     histskyLiMa.Write();
    }

   rootfile.Close();

   delete fitsig;
   delete fitbgpar;
   delete fitbg4;

  cout << " end of one of the grids " << endl;
}